CoCalc -- Qoutput.cpp

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Path: gap4r8 / pkg / NormalizInterface-1.0.2 / Normaliz.git / Qsource / Qoutput.cpp
Views: ⁴¹⁸³⁸⁴
1
/*
2
 * Normaliz
3
 * Copyright (C) 2007-2014  Winfried Bruns, Bogdan Ichim, Christof Soeger
4
 * This program is free software: you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation, either version 3 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
 *
17
 * As an exception, when this program is distributed through (i) the App Store
18
 * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
19
 * by Google Inc., then that store may impose any digital rights management,
20
 * device limits and/or redistribution restrictions that are required by its
21
 * terms of service.
22
 */
23

24
//---------------------------------------------------------------------------
25

26
#include <stdlib.h>
27
#include <vector>
28
#include <list>
29
#include <fstream>
30
#include <iostream>
31
#include <string>
32
#include <algorithm>
33

34
#include "Qoutput.h"
35
#include "libQnormaliz/Qgeneral.h"
36
#include "libQnormaliz/Qmatrix.h"
37
#include "libQnormaliz/Qvector_operations.h"
38
#include "libQnormaliz/Qmap_operations.h"
39

40
//---------------------------------------------------------------------------
41

42
template<typename Number>
43
Output<Number>::Output(){
44
    out=true;
45
    inv=false;
46
    ext=false;
47
    esp=false;
48
    typ=false;
49
    egn=false;
50
    gen=false;
51
    cst=false;
52
    tri=false;
53
    tgn=false;
54
    ht1=false;
55
    dec=false;
56
    lat=false;
57
    mod=false;
58
    msp=false;
59
    lattice_ideal_input = false;
60
}
61

62
//---------------------------------------------------------------------------
63

64
template<typename Number>
65
void Output<Number>::set_lattice_ideal_input(bool value){
66
    lattice_ideal_input=value;
67
}
68

69
//---------------------------------------------------------------------------
70

71
template<typename Number>
72
void Output<Number>::read() const{
73
    cout<<"\nname="<<name<<"\n";
74
    cout<<"\nout="<<out<<"\n";
75
    cout<<"\ninv="<<inv<<"\n";
76
    cout<<"\next="<<ext<<"\n";
77
    cout<<"\nesp="<<esp<<"\n";
78
    cout<<"\ntyp="<<typ<<"\n";
79
    cout<<"\negn="<<egn<<"\n";
80
    cout<<"\ngen="<<gen<<"\n";
81
    cout<<"\ncst="<<cst<<"\n";
82
    cout<<"\ntri="<<tri<<"\n";
83
    cout<<"\ntgn="<<tgn<<"\n";
84
    cout<<"\nht1="<<ht1<<"\n";
85
    cout<<"\nResult is:\n";
86
    Result->print();
87
}
88

89
//---------------------------------------------------------------------------
90

91
template<typename Number>
92
void Output<Number>::set_name(const string& n){
93
    name=n;
94
}
95

96
//---------------------------------------------------------------------------
97

98
template<typename Number>
99
void Output<Number>::setCone(Cone<Number> & C) {
100
    this->Result = &C;
101
    dim = Result->getEmbeddingDim();
102
    homogeneous = !Result->isInhomogeneous();
103
    if (homogeneous) {
104
        of_cone       = "";
105
        of_monoid     = "";
106
        of_polyhedron = "";
107
    } else {
108
        of_cone       = " of recession cone";
109
        of_monoid     = " of recession monoid";
110
        of_polyhedron = " of polyhedron (homogenized)";
111
    }
112
}
113

114
//---------------------------------------------------------------------------
115

116
template<typename Number>
117
void Output<Number>::set_write_out(const bool& flag){
118
    out=flag;
119
}
120

121
//---------------------------------------------------------------------------
122

123
template<typename Number>
124
void Output<Number>::set_write_inv(const bool& flag){
125
    inv=flag;
126
}
127

128
//---------------------------------------------------------------------------
129

130
template<typename Number>
131
void Output<Number>::set_write_ext(const bool& flag){
132
    ext=flag;
133
}
134

135
//---------------------------------------------------------------------------
136

137
template<typename Number>
138
void Output<Number>::set_write_esp(const bool& flag){
139
    esp=flag;
140
}
141

142
//---------------------------------------------------------------------------
143

144
template<typename Number>
145
void Output<Number>::set_write_typ(const bool& flag){
146
    typ=flag;
147
}
148

149
//---------------------------------------------------------------------------
150

151
template<typename Number>
152
void Output<Number>::set_write_egn(const bool& flag){
153
    egn=flag;
154
}
155

156
//---------------------------------------------------------------------------
157

158
template<typename Number>
159
void Output<Number>::set_write_gen(const bool& flag){
160
    gen=flag;
161
}
162

163
//---------------------------------------------------------------------------
164

165
template<typename Number>
166
void Output<Number>::set_write_cst(const bool& flag){
167
    cst=flag;
168
}
169

170
//---------------------------------------------------------------------------
171

172
template<typename Number>
173
void Output<Number>::set_write_tri(const bool& flag) {
174
    tri=flag;
175
}
176

177
//---------------------------------------------------------------------------
178

179
template<typename Number>
180
void Output<Number>::set_write_tgn(const bool& flag) {
181
    tgn=flag;
182
}
183

184
//---------------------------------------------------------------------------
185

186
template<typename Number>
187
void Output<Number>::set_write_ht1(const bool& flag) {
188
    ht1=flag;
189
}
190

191
//---------------------------------------------------------------------------
192

193
template<typename Number>
194
void Output<Number>::set_write_dec(const bool& flag) {
195
    dec=flag;
196
}
197

198
//---------------------------------------------------------------------------
199

200
template<typename Number>
201
void Output<Number>::set_write_mod(const bool& flag) {
202
    mod=flag;
203
}
204

205
//---------------------------------------------------------------------------
206

207
template<typename Number>
208
void Output<Number>::set_write_lat(const bool& flag) {
209
    lat=flag;
210
}
211

212

213
//---------------------------------------------------------------------------
214

215
template<typename Number>
216
void Output<Number>::set_write_msp(const bool& flag) {
217
    msp=flag;
218
}
219
//---------------------------------------------------------------------------
220

221
template<typename Number>
222
void Output<Number>::set_write_extra_files(){
223
    out=true;
224
    inv=true;
225
    gen=true;
226
    cst=true;
227
}
228

229
//---------------------------------------------------------------------------
230

231
template<typename Number>
232
void Output<Number>::set_write_all_files(){
233
    out=true;
234
    inv=true;
235
    ext=true;
236
    esp=true;
237
    typ=true;
238
    egn=true;
239
    gen=true;
240
    cst=true;
241
    ht1=true;
242
    lat=true;
243
    mod=true;
244
    msp=true;
245
}
246

247

248

249
//---------------------------------------------------------------------------
250

251
template<typename Number>
252
void Output<Number>::write_matrix_ext(const Matrix<Number>& M) const{
253
    if (ext==true) {
254
        M.print(name,"ext");
255
    }
256
}
257

258
//---------------------------------------------------------------------------
259

260
template<typename Number>
261
void Output<Number>::write_matrix_mod(const Matrix<Number>& M) const{
262
    if (mod==true) {
263
        M.print(name,"mod");
264
    }
265
}
266

267

268
//---------------------------------------------------------------------------
269

270
template<typename Number>
271
void Output<Number>::write_matrix_lat(const Matrix<Number>& M) const{
272
    if (ext==true) {
273
        M.print(name,"lat");
274
    }
275
}
276

277
//---------------------------------------------------------------------------
278

279
template<typename Number>
280
void Output<Number>::write_matrix_esp(const Matrix<Number>& M) const{
281
    if (esp==true) {
282
        M.print(name,"esp");
283
    }
284
}
285

286
//---------------------------------------------------------------------------
287

288
template<typename Number>
289
void Output<Number>::write_matrix_typ(const Matrix<Number>& M) const{
290
    if (typ==true) {
291
        M.print(name,"typ");
292
    }
293
}
294

295
//---------------------------------------------------------------------------
296

297
template<typename Number>
298
void Output<Number>::write_matrix_egn(const Matrix<Number>& M) const {
299
    if (egn==true) {
300
        M.print(name,"egn");
301
    }
302
}
303

304
//---------------------------------------------------------------------------
305

306
template<typename Number>
307
void Output<Number>::write_matrix_gen(const Matrix<Number>& M) const {
308
    if (gen==true) {
309
        M.print(name,"gen");
310
    }
311
}
312
//---------------------------------------------------------------------------
313

314
template<typename Number>
315
void Output<Number>::write_matrix_msp(const Matrix<Number>& M) const {
316
    if (msp==true) {
317
        M.print(name,"msp");
318
    }
319
}
320

321
//---------------------------------------------------------------------------
322

323
template<typename Number>
324
void Output<Number>::write_tri() const{
325
    if (tri==true) {
326
        string file_name = name+".tri";
327
        ofstream out(file_name.c_str());
328

329
        const vector< pair<vector<libQnormaliz::key_t>,Number> >& Tri = Result->getTriangulation();
330
        typename vector< pair<vector<libQnormaliz::key_t>,Number> >::const_iterator tit = Tri.begin();        
331
        const vector<vector<bool> >& Dec = Result->isComputed(ConeProperty::ConeDecomposition) ?
332
                Result->getOpenFacets() : vector<vector<bool> >();
333
        typename vector< vector<bool> >::const_iterator idd = Dec.begin();
334

335
        out << Tri.size() << endl;
336
        size_t nr_extra_entries=1;
337
        if (Result->isComputed(ConeProperty::ConeDecomposition))
338
            nr_extra_entries+=Result->getSublattice().getRank();
339
        out << Result->getSublattice().getRank()+nr_extra_entries << endl; //works also for empty list
340

341
        for(; tit != Tri.end(); ++tit) {
342
            for (size_t i=0; i<tit->first.size(); i++) {
343
                out << tit->first[i] +1 << " ";
344
            }
345
            out << "   " << tit->second;
346
            if(Result->isComputed(ConeProperty::ConeDecomposition)){
347
                out << "   ";
348
                for (size_t i=0; i<tit->first.size(); i++) {
349
                    out << " " << (*idd)[i];
350
                }                
351
                idd++;
352
            }
353
            out << endl;
354
        }
355
        if (Result->isTriangulationNested()) out << "nested" << endl;
356
        else out << "plain" << endl;
357
        if (Result->isTriangulationPartial()) out << "partial" << endl;
358
        out.close();
359
    }
360
}
361

362
//---------------------------------------------------------------------------
363

364
template<typename Number>
365
void Output<Number>::write_matrix_ht1(const Matrix<Number>& M) const{
366
    if (ht1==true) {
367
        M.print(name,"ht1");
368
    }
369
}
370

371
//---------------------------------------------------------------------------
372

373
string is_maximal(long a, long b) {
374
    return (a == b) ? " (maximal)" : "";
375
}
376

377
//---------------------------------------------------------------------------
378

379
template<typename Number>
380
void Output<Number>::write_inv_file() const{
381
    if (inv==true) {//printing .inv file
382
        
383
        string name_open=name+".inv";                              //preparing output files
384
        const char* file=name_open.c_str();
385
        ofstream inv(file);
386

387
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
388
            inv << "integer number_vertices_polyhedron = "
389
                << Result->getNrVerticesOfPolyhedron() << endl;
390
        }
391
        if (Result->isComputed(ConeProperty::ExtremeRays)) {
392
            size_t nr_ex_rays = Result->getNrExtremeRays();
393
            inv<<"integer number_extreme_rays = "<<nr_ex_rays<<endl;
394
        }
395
        if (Result->isComputed(ConeProperty::MaximalSubspace)) {
396
            size_t dim_max_subspace = Result->getDimMaximalSubspace();
397
            inv<<"integer dim_max_subspace = "<<dim_max_subspace<<endl;
398
        }
399

400
        inv << "integer embedding_dim = " << dim << endl;
401
        if (!homogeneous){
402
            if (Result->isComputed(ConeProperty::AffineDim))
403
                inv << "integer affine_dim_polyhedron = " << Result->getAffineDim() << endl;
404
            if (Result->isComputed(ConeProperty::RecessionRank))
405
                inv << "integer recession_rank = "  << Result->getRecessionRank() << endl;
406
        }
407
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) { 
408
            inv<<"integer number_support_hyperplanes = "<<Result->getNrSupportHyperplanes()<<endl;
409
        }
410
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
411
            inv << "integer size_triangulation = " << Result->getTriangulationSize() << endl;
412
        }
413
        if (Result->isComputed(ConeProperty::TriangulationDetSum)) {
414
            inv << "integer sum_dets = " << Result->getTriangulationDetSum() << endl;
415
        }
416
        
417
        if (!Result->isComputed(ConeProperty::Dehomogenization)) {
418
            inv << "boolean inhomogeneous = false" << endl;
419
        }
420
        else {
421
            inv << "boolean inhomogeneous = true" << endl;
422
            vector<Number> Linear_Form = Result->getDehomogenization();
423
            inv << "vector " << Linear_Form.size()
424
                << " dehomogenization = " << Linear_Form;
425
        }
426
        
427
        inv.close();
428
    }
429
}
430

431

432
//---------------------------------------------------------------------------
433

434
template<typename Number>
435
void Output<Number>::write_files() const {
436
    vector<libQnormaliz::key_t> rees_ideal_key;
437

438
    if (esp && Result->isComputed(ConeProperty::SupportHyperplanes) && Result->isComputed(ConeProperty::Sublattice)) {
439
        //write the suport hyperplanes of the full dimensional cone
440
        const Sublattice_Representation<Number>& BasisChange = Result->getSublattice();
441
        Matrix<Number> Support_Hyperplanes_Full_Cone = BasisChange.to_sublattice_dual(Result->getSupportHyperplanesMatrix());
442
        // Support_Hyperplanes_Full_Cone.print(name,"esp");
443
        string esp_string = name+".esp";
444
        const char* esp_file = esp_string.c_str();
445
        ofstream esp_out(esp_file);
446
        Support_Hyperplanes_Full_Cone.print(esp_out);
447
        esp_out << "inequalities" << endl;
448
        if (Result->isComputed(ConeProperty::Grading)) {
449
            esp_out << 1 << endl << Result->getRank() << endl;
450
        }
451
        if (Result->isComputed(ConeProperty::Dehomogenization)) {
452
            esp_out << 1 << endl << Result->getRank() << endl;
453
            esp_out << BasisChange.to_sublattice_dual(Result->getDehomogenization());
454
            esp_out << "dehomogenization" << endl;
455
        }
456
        esp_out.close();
457
    }
458
    if (tgn && Result->isComputed(ConeProperty::Generators))
459
        Result->getGeneratorsMatrix().print(name,"tgn");
460
    if (tri && Result->isComputed(ConeProperty::Triangulation)) {     //write triangulation
461
        write_tri();
462
    }
463

464
    if (out==true) {  //printing .out file
465
        string name_open=name+".out";                              //preparing output files
466
        const char* file=name_open.c_str();
467
        ofstream out(file);
468
        if(out.fail()){
469
            errorOutput() << "Cannot write to output file." << endl;
470
            exit(1);
471
        }
472

473
        // write "header" of the .out file
474
        size_t nr_orig_gens = 0;
475
        if (lattice_ideal_input) {
476
            nr_orig_gens = Result->getNrOriginalMonoidGenerators();
477
            out << nr_orig_gens <<" original generators of the toric ring"<<endl;
478
        }
479

480
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
481
            out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron" << endl;
482
        }
483
        if (Result->isComputed(ConeProperty::ExtremeRays)) {
484
            out << Result->getNrExtremeRays() <<" extreme rays" << of_cone << endl;
485
        }
486
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) {
487
            out << Result->getNrSupportHyperplanes() <<" support hyperplanes"
488
                << of_polyhedron << endl;
489
        }
490
        out<<endl;
491

492
        out << "embedding dimension = " << dim << endl;
493
        if (homogeneous) {
494
            if (Result->isComputed(ConeProperty::Sublattice)) {
495
                auto rank = Result->getRank();
496
                out << "rank = "<< rank << is_maximal(rank,dim) << endl;
497
            }
498
        } else { // now inhomogeneous case
499
            if (Result->isComputed(ConeProperty::AffineDim))
500
                out << "affine dimension of the polyhedron = "
501
                    << Result->getAffineDim() << is_maximal(Result->getAffineDim(),dim-1) << endl;
502
            if (Result->isComputed(ConeProperty::RecessionRank))
503
                out << "rank of recession monoid = "  << Result->getRecessionRank() << endl;
504
        }
505
        
506
        if(Result->isComputed(ConeProperty::MaximalSubspace)){
507
            size_t dim_max_subspace=Result->getDimMaximalSubspace();
508
            if(dim_max_subspace>0)
509
                out << "dimension of maximal subspace = " << dim_max_subspace << endl;      
510
        }
511
            
512
        out << endl;
513
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
514
            out << "size of ";
515
            if (Result->isTriangulationNested()) out << "nested ";
516
            if (Result->isTriangulationPartial()) out << "partial ";
517
            out << "triangulation   = " << Result->getTriangulationSize() << endl;
518
        }
519
        if (Result->isComputed(ConeProperty::TriangulationDetSum)) {
520
            out << "resulting sum of |det|s = " << Result->getTriangulationDetSum() << endl;
521
        }
522
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
523
            out << endl;
524
        }
525
        if ( Result->isComputed(ConeProperty::Dehomogenization) ) {
526
            out << "dehomogenization:" << endl
527
                << Result->getDehomogenization() << endl;
528
        }
529
 
530
        out << "***********************************************************************"
531
            << endl << endl;
532
        
533
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
534
            out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron:" << endl;
535
            Result->getVerticesOfPolyhedronMatrix().pretty_print(out);
536
            out << endl;
537
        }
538
        if (Result->isComputed(ConeProperty::ExtremeRays)) {
539
            out << Result->getNrExtremeRays() << " extreme rays" << of_cone << ":" << endl;
540
            Result->getExtremeRaysMatrix().pretty_print(out);
541
            out << endl;
542
            if (ext) {
543
                // for the .gen file we append the vertices of polyhedron if there are any
544
                if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
545
                    Matrix<Number> Extreme_Rays(Result->getExtremeRaysMatrix());
546
                    Extreme_Rays.append(Result->getVerticesOfPolyhedronMatrix());
547
                    write_matrix_ext(Extreme_Rays);
548
                } else {
549
                    write_matrix_ext(Result->getExtremeRaysMatrix());
550
                }
551
            }
552
        }
553
        
554
        if(Result->isComputed(ConeProperty::MaximalSubspace) && Result->getDimMaximalSubspace()>0){
555
            out << Result->getDimMaximalSubspace() <<" basis elements of maximal subspace:" << endl;
556
            Result->getMaximalSubspaceMatrix().pretty_print(out);
557
            out << endl;
558
            if(msp)
559
                write_matrix_msp(Result->getMaximalSubspaceMatrix());            
560
        }
561

562
        //write constrains (support hyperplanes, congruences, equations)
563

564
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) {
565
            const Matrix<Number>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();
566
            out << Support_Hyperplanes.nr_of_rows() <<" support hyperplanes" 
567
                << of_polyhedron << ":" << endl;
568
            Support_Hyperplanes.pretty_print(out);
569
            out << endl;
570
        }
571
        if (Result->isComputed(ConeProperty::Sublattice)) {
572
            const Sublattice_Representation<Number>& BasisChange = Result->getSublattice();
573
            //equations
574
            const Matrix<Number>& Equations = BasisChange.getEquationsMatrix();
575
            size_t nr_of_equ = Equations.nr_of_rows();
576
            if (nr_of_equ > 0) {
577
                out << nr_of_equ <<" equations:" <<endl;
578
                Equations.pretty_print(out);
579
                out << endl;
580
            }
581
            
582
            //lattice
583
            const Matrix<Number>& LatticeBasis = BasisChange.getEmbeddingMatrix();
584
            size_t nr_of_latt = LatticeBasis.nr_of_rows();
585
            if (nr_of_latt < dim) {
586
                out << nr_of_latt <<" basis elements of lattice:" <<endl;
587
                LatticeBasis.pretty_print(out);
588
                out << endl;
589
            }
590
            if(lat)
591
                write_matrix_lat(LatticeBasis);
592
            
593

594
            if (cst && Result->isComputed(ConeProperty::SupportHyperplanes)) {
595
                const Matrix<Number>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();
596
                string cst_string = name+".cst";
597
                const char* cst_file = cst_string.c_str();
598
                ofstream cst_out(cst_file);
599
    
600
                Support_Hyperplanes.print(cst_out);
601
                cst_out<<"inequalities"<<endl;
602
                Equations.print(cst_out);
603
                cst_out<<"equations"<<endl;
604

605
                if (Result->isComputed(ConeProperty::Dehomogenization)) {
606
                    cst_out << 1 << endl << dim << endl;
607
                    cst_out << Result->getDehomogenization();
608
                    cst_out << "dehomogenization" << endl;
609
                }
610
                cst_out.close();
611
            }
612
        }
613
        
614
        out.close();
615
    }
616

617
    write_inv_file();
618
}
619

620

621
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place.

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
