CoCalc -- output.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 / source / output.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
#include <iomanip>
34

35
#include "output.h"
36
#include "libnormaliz/general.h"
37
#include "libnormaliz/matrix.h"
38
#include "libnormaliz/vector_operations.h"
39
#include "libnormaliz/map_operations.h"
40

41
//---------------------------------------------------------------------------
42

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

64
//---------------------------------------------------------------------------
65

66
template<typename Integer>
67
void Output<Integer>::set_lattice_ideal_input(bool value){
68
    lattice_ideal_input=value;
69
}
70

71
//---------------------------------------------------------------------------
72

73
template<typename Integer>
74
void Output<Integer>::set_no_ext_rays_output(){
75
    no_ext_rays_output=true;
76
}
77

78
//---------------------------------------------------------------------------
79

80
template<typename Integer>
81
void Output<Integer>::read() const{
82
    cout<<"\nname="<<name<<"\n";
83
    cout<<"\nout="<<out<<"\n";
84
    cout<<"\ninv="<<inv<<"\n";
85
    cout<<"\next="<<ext<<"\n";
86
    cout<<"\nesp="<<esp<<"\n";
87
    cout<<"\ntyp="<<typ<<"\n";
88
    cout<<"\negn="<<egn<<"\n";
89
    cout<<"\ngen="<<gen<<"\n";
90
    cout<<"\ncst="<<cst<<"\n";
91
    cout<<"\ntri="<<tri<<"\n";
92
    cout<<"\ntgn="<<tgn<<"\n";
93
    cout<<"\nht1="<<ht1<<"\n";
94
    cout<<"\nResult is:\n";
95
    Result->print();
96
}
97

98
//---------------------------------------------------------------------------
99

100
template<typename Integer>
101
void Output<Integer>::set_name(const string& n){
102
    name=n;
103
}
104

105
//---------------------------------------------------------------------------
106

107
template<typename Integer>
108
void Output<Integer>::setCone(Cone<Integer> & C) {
109
    this->Result = &C;
110
    dim = Result->getEmbeddingDim();
111
    homogeneous = !Result->isInhomogeneous();
112
    if (homogeneous) {
113
        of_cone       = "";
114
        of_monoid     = "";
115
        of_polyhedron = "";
116
    } else {
117
        of_cone       = " of recession cone";
118
        of_monoid     = " of recession monoid";
119
        of_polyhedron = " of polyhedron (homogenized)";
120
    }
121
}
122

123
//---------------------------------------------------------------------------
124

125
template<typename Integer>
126
void Output<Integer>::set_write_out(const bool& flag){
127
    out=flag;
128
}
129

130
//---------------------------------------------------------------------------
131

132
template<typename Integer>
133
void Output<Integer>::set_write_inv(const bool& flag){
134
    inv=flag;
135
}
136

137
//---------------------------------------------------------------------------
138

139
template<typename Integer>
140
void Output<Integer>::set_write_ext(const bool& flag){
141
    ext=flag;
142
}
143

144
//---------------------------------------------------------------------------
145

146
template<typename Integer>
147
void Output<Integer>::set_write_esp(const bool& flag){
148
    esp=flag;
149
}
150

151
//---------------------------------------------------------------------------
152

153
template<typename Integer>
154
void Output<Integer>::set_write_typ(const bool& flag){
155
    typ=flag;
156
}
157

158
//---------------------------------------------------------------------------
159

160
template<typename Integer>
161
void Output<Integer>::set_write_egn(const bool& flag){
162
    egn=flag;
163
}
164

165
//---------------------------------------------------------------------------
166

167
template<typename Integer>
168
void Output<Integer>::set_write_gen(const bool& flag){
169
    gen=flag;
170
}
171

172
//---------------------------------------------------------------------------
173

174
template<typename Integer>
175
void Output<Integer>::set_write_cst(const bool& flag){
176
    cst=flag;
177
}
178

179
//---------------------------------------------------------------------------
180

181
template<typename Integer>
182
void Output<Integer>::set_write_tri(const bool& flag) {
183
    tri=flag;
184
}
185

186
//---------------------------------------------------------------------------
187

188
template<typename Integer>
189
void Output<Integer>::set_write_tgn(const bool& flag) {
190
    tgn=flag;
191
}
192

193
//---------------------------------------------------------------------------
194

195
template<typename Integer>
196
void Output<Integer>::set_write_ht1(const bool& flag) {
197
    ht1=flag;
198
}
199

200
//---------------------------------------------------------------------------
201

202
template<typename Integer>
203
void Output<Integer>::set_write_dec(const bool& flag) {
204
    dec=flag;
205
}
206

207
//---------------------------------------------------------------------------
208

209
template<typename Integer>
210
void Output<Integer>::set_write_mod(const bool& flag) {
211
    mod=flag;
212
}
213

214
//---------------------------------------------------------------------------
215

216
template<typename Integer>
217
void Output<Integer>::set_write_lat(const bool& flag) {
218
    lat=flag;
219
}
220

221

222
//---------------------------------------------------------------------------
223

224
template<typename Integer>
225
void Output<Integer>::set_write_msp(const bool& flag) {
226
    msp=flag;
227
}
228
//---------------------------------------------------------------------------
229

230
template<typename Integer>
231
void Output<Integer>::set_write_extra_files(){
232
    out=true;
233
    inv=true;
234
    gen=true;
235
    cst=true;
236
}
237

238
//---------------------------------------------------------------------------
239

240
template<typename Integer>
241
void Output<Integer>::set_write_all_files(){
242
    out=true;
243
    inv=true;
244
    ext=true;
245
    esp=true;
246
    typ=true;
247
    egn=true;
248
    gen=true;
249
    cst=true;
250
    ht1=true;
251
    lat=true;
252
    mod=true;
253
    msp=true;
254
}
255

256

257

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

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

267
//---------------------------------------------------------------------------
268

269
template<typename Integer>
270
void Output<Integer>::write_matrix_mod(const Matrix<Integer>& M) const{
271
    if (mod==true) {
272
        M.print(name,"mod");
273
    }
274
}
275

276

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

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

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

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

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

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

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

306
template<typename Integer>
307
void Output<Integer>::write_matrix_egn(const Matrix<Integer>& M) const {
308
    if (egn==true) {
309
        M.print(name,"egn");
310
    }
311
}
312

313
//---------------------------------------------------------------------------
314

315
template<typename Integer>
316
void Output<Integer>::write_matrix_gen(const Matrix<Integer>& M) const {
317
    if (gen==true) {
318
        M.print(name,"gen");
319
    }
320
}
321
//---------------------------------------------------------------------------
322

323
template<typename Integer>
324
void Output<Integer>::write_matrix_msp(const Matrix<Integer>& M) const {
325
    if (msp==true) {
326
        M.print(name,"msp");
327
    }
328
}
329

330
//---------------------------------------------------------------------------
331

332
template<typename Integer>
333
void Output<Integer>::write_tri() const{
334
    if (tri==true) {
335
        string file_name = name+".tri";
336
        ofstream out(file_name.c_str());
337

338
        const vector< pair<vector<libnormaliz::key_t>,Integer> >& Tri = Result->getTriangulation();
339
        typename vector< pair<vector<libnormaliz::key_t>,Integer> >::const_iterator tit = Tri.begin();        
340
        const vector<vector<bool> >& Dec = Result->isComputed(ConeProperty::ConeDecomposition) ?
341
                Result->getOpenFacets() : vector<vector<bool> >();
342
        typename vector< vector<bool> >::const_iterator idd = Dec.begin();
343

344
        out << Tri.size() << endl;
345
        size_t nr_extra_entries=1;
346
        if (Result->isComputed(ConeProperty::ConeDecomposition))
347
            nr_extra_entries+=Result->getSublattice().getRank();
348
        out << Result->getSublattice().getRank()+nr_extra_entries << endl; //works also for empty list
349

350
        for(; tit != Tri.end(); ++tit) {
351
            for (size_t i=0; i<tit->first.size(); i++) {
352
                out << tit->first[i] +1 << " ";
353
            }
354
            out << "   " << tit->second;
355
            if(Result->isComputed(ConeProperty::ConeDecomposition)){
356
                out << "   ";
357
                for (size_t i=0; i<tit->first.size(); i++) {
358
                    out << " " << (*idd)[i];
359
                }                
360
                idd++;
361
            }
362
            out << endl;
363
        }
364
        if (Result->isTriangulationNested()) out << "nested" << endl;
365
        else out << "plain" << endl;
366
        if (Result->isTriangulationPartial()) out << "partial" << endl;
367
        out.close();
368
    }
369
}
370

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

373

374
template<typename Integer>
375
void Output<Integer>::write_Stanley_dec() const {
376
    if (dec && Result->isComputed(ConeProperty::StanleyDec)) {
377
        ofstream out((name+".dec").c_str());
378

379
        if (Result->isComputed(ConeProperty::InclusionExclusionData)) {
380
            const vector< pair<vector<libnormaliz::key_t>, long> >& InExData = Result->getInclusionExclusionData();
381
            out << "in_ex_data" << endl;
382
            out << InExData.size() << endl;
383
            for (size_t i=0; i<InExData.size(); ++i) {
384
                out << InExData[i].first.size() << " ";
385
                for (size_t j=0; j<InExData[i].first.size(); ++j) {
386
                    out << InExData[i].first[j]+1 << " ";
387
                }
388
                out << InExData[i].second << endl;  
389
            }
390
        }
391

392
        out << "Stanley_dec" << endl;
393
        list<STANLEYDATA_int>& StanleyDec = Result->getStanleyDec_mutable();
394
        auto S = StanleyDec.begin();
395
        size_t i;
396

397
        out << StanleyDec.size() << endl; 
398
        for (; S!=StanleyDec.end(); ++S) {
399
            for (i=0; i < S->key.size(); ++i)
400
                out << S->key[i]+1 <<" ";
401
            out << endl;
402
            S->offsets.print(out);
403
            out << endl;
404
        }
405
        out.close();
406
    }
407
}
408

409
//---------------------------------------------------------------------------
410

411
template<typename Integer>
412
void Output<Integer>::write_matrix_ht1(const Matrix<Integer>& M) const{
413
    if (ht1==true) {
414
        M.print(name,"ht1");
415
    }
416
}
417

418
//---------------------------------------------------------------------------
419

420
string is_maximal(long a, long b) {
421
    return (a == b) ? " (maximal)" : "";
422
}
423

424
//---------------------------------------------------------------------------
425

426
template<typename Integer>
427
void Output<Integer>::write_inv_file() const{
428
    if (inv==true) {//printing .inv file
429
        size_t i;
430
        string name_open=name+".inv";                              //preparing output files
431
        const char* file=name_open.c_str();
432
        ofstream inv(file);
433

434
        if (Result->isComputed(ConeProperty::ModuleGenerators)) {
435
            inv << "integer number_module_generators = "
436
                << Result->getNrModuleGenerators() << endl;
437
        }
438
        if (Result->isComputed(ConeProperty::HilbertBasis)) {
439
            inv<<"integer hilbert_basis_elements = "<<Result->getNrHilbertBasis()<<endl;
440
        }
441

442
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
443
            inv << "integer number_vertices_polyhedron = "
444
                << Result->getNrVerticesOfPolyhedron() << endl;
445
        }
446
        if (Result->isComputed(ConeProperty::ExtremeRays)) {
447
            size_t nr_ex_rays = Result->getNrExtremeRays();
448
            inv<<"integer number_extreme_rays = "<<nr_ex_rays<<endl;
449
        }
450
        if (Result->isComputed(ConeProperty::MaximalSubspace)) {
451
            size_t dim_max_subspace = Result->getDimMaximalSubspace();
452
            inv<<"integer dim_max_subspace = "<<dim_max_subspace<<endl;
453
        }
454
        if (Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {
455
            inv << "integer number_module_generators_original_monoid = "
456
                << Result->getNrModuleGeneratorsOverOriginalMonoid() << endl;
457
        }
458

459
        inv << "integer embedding_dim = " << dim << endl;
460
        if (homogeneous) {
461
            if (Result->isComputed(ConeProperty::Sublattice)) {
462
                inv << "integer rank = " << Result->getRank() << endl;
463
                inv << "integer external_index = " << Result->getSublattice().getExternalIndex() << endl;
464
            }
465
        } else {
466
            if (Result->isComputed(ConeProperty::AffineDim))
467
                inv << "integer affine_dim_polyhedron = " << Result->getAffineDim() << endl;
468
            if (Result->isComputed(ConeProperty::RecessionRank))
469
                inv << "integer recession_rank = "  << Result->getRecessionRank() << endl;
470
        }
471

472
        if(Result->isComputed(ConeProperty::OriginalMonoidGenerators)){
473
            inv << "integer internal_index = " << Result->getIndex() << endl;           
474
        }
475
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) { 
476
            inv<<"integer number_support_hyperplanes = "<<Result->getNrSupportHyperplanes()<<endl;
477
        }
478
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
479
            inv << "integer size_triangulation = " << Result->getTriangulationSize() << endl;
480
        }
481
        if (Result->isComputed(ConeProperty::TriangulationDetSum)) {
482
            inv << "integer sum_dets = " << Result->getTriangulationDetSum() << endl;
483
        }
484

485
        if (Result->isComputed(ConeProperty::IsIntegrallyClosed)) {
486
            if(Result->isIntegrallyClosed())
487
                inv << "boolean integrally_closed = true" << endl;
488
            else
489
                inv << "boolean integrally_closed = false" << endl;
490
        }
491
        
492
        if (!Result->isComputed(ConeProperty::Dehomogenization)) {
493
            inv << "boolean inhomogeneous = false" << endl;
494
        }
495
        else {
496
            inv << "boolean inhomogeneous = true" << endl;
497
            vector<Integer> Linear_Form = Result->getDehomogenization();
498
            inv << "vector " << Linear_Form.size()
499
                << " dehomogenization = " << Linear_Form;
500
        }
501
        if (Result->isComputed(ConeProperty::Grading)==false) {
502
            if (Result->isComputed(ConeProperty::ExtremeRays)) {
503
                inv<<"boolean graded = "<<"false"<<endl;
504
            }
505
        }
506
        else {
507
            inv<<"boolean graded = "<<"true"<<endl;
508
            if (Result->isComputed(ConeProperty::Deg1Elements)) {
509
                inv<<"integer degree_1_elements = "<<Result->getNrDeg1Elements()<<endl;
510
            }
511
            vector<Integer> Linear_Form = Result->getGrading();
512
            inv<<"vector "<<Linear_Form.size()<<" grading = ";
513
            for (i = 0; i < Linear_Form.size(); i++) {
514
                inv<<Linear_Form[i]<<" ";
515
            }
516
            inv<<endl;
517
            inv <<"integer grading_denom = " << Result->getGradingDenom() << endl;
518
            if (Result->isComputed(ConeProperty::Multiplicity)){
519
                mpq_class mult = Result->getMultiplicity();
520
                inv <<"integer multiplicity = "      << mult.get_num() << endl;
521
                inv <<"integer multiplicity_denom = "<< mult.get_den() << endl;
522
            }
523
            if (Result->isComputed(ConeProperty::WeightedEhrhartSeries)){
524
                const HilbertSeries& HS=Result->getIntData().getWeightedEhrhartSeries().first;
525
                
526
                inv <<"vector "<< HS.getNum().size() <<" weighted_ehrhart_series_num = " << HS.getNum();
527
                inv <<"integer num_common_denom = " << Result->getIntData().getWeightedEhrhartSeries().second << endl;
528
                inv <<"vector "<< to_vector(HS.getDenom()).size() <<" weighted_ehrhart_series_num = " 
529
                        << to_vector(HS.getDenom());
530
                Result->getIntData().computeWeightedEhrhartQuasiPolynomial();
531
                if (Result->getIntData().isWeightedEhrhartQuasiPolynomialComputed()) {
532
                    if(HS.get_nr_coeff_quasipol()>=0){
533
                        inv << "integer nr_coeff_weightedEhrhart__quasipol = " << HS.get_nr_coeff_quasipol() << endl;                       
534
                    }
535
                    vector< vector <mpz_class> > hqp= Result->getIntData().getWeightedEhrhartQuasiPolynomial();
536
                    inv << "matrix " << hqp.size() << " " << hqp[0].size()
537
                        << " weighted_ehrhart_quasipolynomial = ";
538
                    inv << endl << hqp;
539
                    inv << "integer weighted_ehrhart_quasipolynomial_denom = " 
540
                        << Result->getIntData().getWeightedEhrhartQuasiPolynomialDenom() << endl;
541
                }
542
            }
543
            
544
            if (Result->isComputed(ConeProperty::VirtualMultiplicity)){
545
                mpq_class mult = Result->getVirtualMultiplicity();
546
                inv <<"integer virtual_multiplicity = "      << mult.get_num() << endl;
547
                inv <<"integer virtual_multiplicity_denom = "<< mult.get_den() << endl;
548
            }
549
            if (Result->isComputed(ConeProperty::Integral)){
550
                mpq_class mult = Result->getIntegral();
551
                inv <<"integer integral = "      << mult.get_num() << endl;
552
                inv <<"integer integral_denom = "<< mult.get_den() << endl;
553
            }
554
            if (Result->isComputed(ConeProperty::HilbertSeries)) {
555
                const HilbertSeries& HS = Result->getHilbertSeries();
556
                vector<mpz_class> HSnum;
557
                vector<denom_t> HSdenom;
558
                if(Result->isComputed(ConeProperty::HSOP)){
559
                    HSnum = HS.getHSOPNum();
560
                    HSdenom = to_vector(HS.getHSOPDenom());
561
                    
562
                } else {
563
                    HSnum = HS.getNum();
564
                    HSdenom = to_vector(HS.getDenom());
565
                }
566
                inv <<"vector "<< HSnum.size() <<" hilbert_series_num = ";
567
                inv << HSnum;
568
                inv <<"vector "<< HSdenom.size() <<" hilbert_series_denom = ";
569
                    inv << HSdenom;
570
                HS.computeHilbertQuasiPolynomial();
571
                if (HS.isHilbertQuasiPolynomialComputed()) {
572
                    if(HS.get_nr_coeff_quasipol()>=0){
573
                        inv << "integer nr_coeff_hilbert_quasipol = " << HS.get_nr_coeff_quasipol() << endl;                       
574
                    }
575
                    vector< vector <mpz_class> > hqp = HS.getHilbertQuasiPolynomial();
576
                    inv << "matrix " << hqp.size() << " " << hqp[0].size()
577
                        << " hilbert_quasipolynomial = ";
578
                    inv << endl << hqp;
579
                    inv << "integer hilbert_quasipolynomial_denom = " 
580
                        << HS.getHilbertQuasiPolynomialDenom() << endl;
581
                }
582
            }
583
        }
584
        if (Result->isComputed(ConeProperty::IsReesPrimary)) {
585
            if (Result->isReesPrimary()) {
586
                inv<<"boolean primary = true"<<endl;
587
            } else {
588
                inv<<"boolean primary = false"<<endl;
589
            }
590
        }
591
        if (Result->isComputed(ConeProperty::ReesPrimaryMultiplicity)) {
592
            inv<<"integer ideal_multiplicity = "<<Result->getReesPrimaryMultiplicity()<<endl;
593
        }
594
        
595
        if (Result->isComputed(ConeProperty::ClassGroup)) {
596
            inv <<"vector "<< Result->getClassGroup().size() <<" class_group = "<< Result->getClassGroup();
597
        }
598
        
599
        if (Result->isComputed(ConeProperty::IsGorenstein)) {
600
            if(Result->isGorenstein()){
601
                inv << "boolean Gorenstein = true" << endl;
602
                inv <<"vector "<< Result->getGeneratorOfInterior().size() <<" generator_of_interior = "<< Result->getGeneratorOfInterior();                
603
            }
604
            else
605
                inv << "boolean Gorenstein = false" << endl;
606
        }
607
        
608
        inv.close();
609
    }
610
}
611

612
//---------------------------------------------------------------------------
613

614
template<typename Integer>
615
void Output<Integer>::writeWeightedEhrhartSeries(ofstream& out) const{
616

617
    HilbertSeries HS=Result->getIntData().getWeightedEhrhartSeries().first;
618
    out << "Weighted Ehrhart series:" << endl;
619
    vector<mpz_class> num( HS.getNum());
620
    for(size_t i=0;i<num.size();++i)
621
        out << num[i] << " ";
622
    out << endl << "Common denominator of coefficients: ";
623
    out << Result->getIntData().getWeightedEhrhartSeries().second << endl;
624
    map<long, long> HS_Denom = HS.getDenom();
625
    long nr_factors = 0;
626
    for (map<long, long>::iterator it = HS_Denom.begin(); it!=HS_Denom.end(); ++it) {
627
        nr_factors += it->second;
628
    }
629
    out << "Series denominator with " << nr_factors << " factors:" << endl;
630
    out << HS.getDenom();
631
    if (HS.getShift() != 0) {
632
        out << "shift = " << HS.getShift() << endl << endl;
633
    }
634
    out << endl;
635
    out << "degree of weighted Ehrhart series as rational function = "
636
            << HS.getDegreeAsRationalFunction() << endl << endl;
637
    long period = HS.getPeriod();
638
    if (period == 1) {
639
        out << "Weighted Ehrhart polynomial:" << endl;
640
        for(size_t i=0; i<HS.getHilbertQuasiPolynomial()[0].size();++i)
641
            out << HS.getHilbertQuasiPolynomial()[0][i] << " ";
642
        out << endl;
643
        out << "with common denominator: ";
644
        out << HS.getHilbertQuasiPolynomialDenom()*Result->getIntData().getNumeratorCommonDenom();
645
    } else {
646
        // output cyclonomic representation
647
        out << "Weighted Ehrhart series with cyclotomic denominator:" << endl;
648
        num=HS.getCyclotomicNum();
649
        for(size_t i=0;i<num.size();++i)
650
            out << num[i] << " ";
651
        out << endl << "Common denominator of coefficients: ";
652
        out << Result->getIntData().getWeightedEhrhartSeries().second << endl;
653
        out << "Series cyclotomic denominator:" << endl;
654
        out << HS.getCyclotomicDenom();
655
        out << endl;
656
        // Weighted Ehrhart quasi-polynomial
657
        // vector< vector<mpz_class> > hilbert_quasi_poly = HS.getHilbertQuasiPolynomial();
658
        if (HS.isHilbertQuasiPolynomialComputed()) { 
659
            out<<"Weighted Ehrhart quasi-polynomial of period " << period << ":" << endl;
660
            if(HS.get_nr_coeff_quasipol()>=0){
661
                out << "only " << HS.get_nr_coeff_quasipol() << " highest coefficients computed" << endl;
662
                out << "their common period is " << HS.getHilbertQuasiPolynomial().size() << "." << endl;                        
663
            }
664
            Matrix<mpz_class> HQP(HS.getHilbertQuasiPolynomial());
665
            HQP.pretty_print(out,true);
666
            out<<"with common denominator: "
667
               <<Result->getIntData().getWeightedEhrhartQuasiPolynomialDenom() << endl;
668
        }
669
        else{
670
            out<<"Weighted Ehrhart quasi-polynomial has period " << period << endl;    
671
        }
672
    }
673
    
674
   out  << endl << endl;
675
    
676
   if (HS.isHilbertQuasiPolynomialComputed()){    
677
    long deg=HS.getHilbertQuasiPolynomial()[0].size()-1;
678
        out << "Degree of (quasi)polynomial: " << deg << endl;
679
    
680
    long virtDeg=Result->getRank()+Result->getIntData().getDegreeOfPolynomial()-1;
681

682
    out << endl << "Expected degree: " <<  virtDeg << endl;
683
    }
684
        
685
    if(Result->isComputed(ConeProperty::VirtualMultiplicity)){
686
        out << endl << "Virtual multiplicity: ";
687
        out << Result->getIntData().getVirtualMultiplicity() << endl << endl;
688
    }
689
}
690

691
//---------------------------------------------------------------------------
692
template<typename Integer>
693
void Output<Integer>::write_float(ofstream& out, const Matrix<nmz_float>& mat, size_t nr, size_t nc) const{
694

695
    for(size_t i=0;i<nr;++i){
696
        for(size_t j=0; j< nc;++j){
697
            out << std::setw(10) << mat[i][j];
698
        }
699
        out << endl;
700
    }
701
}
702

703
//---------------------------------------------------------------------------
704

705
template<typename Integer>
706
void Output<Integer>::write_files() const {
707
    size_t i, nr;
708
    vector<libnormaliz::key_t> rees_ideal_key;
709

710
    if (esp && Result->isComputed(ConeProperty::SupportHyperplanes) && Result->isComputed(ConeProperty::Sublattice)) {
711
        //write the suport hyperplanes of the full dimensional cone
712
        const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();
713
        Matrix<Integer> Support_Hyperplanes_Full_Cone = BasisChange.to_sublattice_dual(Result->getSupportHyperplanesMatrix());
714
        // Support_Hyperplanes_Full_Cone.print(name,"esp");
715
        string esp_string = name+".esp";
716
        const char* esp_file = esp_string.c_str();
717
        ofstream esp_out(esp_file);
718
        Support_Hyperplanes_Full_Cone.print(esp_out);
719
        esp_out << "inequalities" << endl;
720
        if (Result->isComputed(ConeProperty::Grading)) {
721
            esp_out << 1 << endl << Result->getRank() << endl;
722
            esp_out << BasisChange.to_sublattice_dual(Result->getGrading());
723
            esp_out << "grading" << endl;
724
        }
725
        if (Result->isComputed(ConeProperty::Dehomogenization)) {
726
            esp_out << 1 << endl << Result->getRank() << endl;
727
            esp_out << BasisChange.to_sublattice_dual(Result->getDehomogenization());
728
            esp_out << "dehomogenization" << endl;
729
        }
730
        esp_out.close();
731
    }
732
    if (tgn && Result->isComputed(ConeProperty::Generators))
733
        Result->getGeneratorsMatrix().print(name,"tgn");
734
    if (tri && Result->isComputed(ConeProperty::Triangulation)) {     //write triangulation
735
        write_tri();
736
    }
737

738
    if (out==true) {  //printing .out file
739
        string name_open=name+".out";                              //preparing output files
740
        const char* file=name_open.c_str();
741
        ofstream out(file);
742
        if(out.fail()){
743
            errorOutput() << "Cannot write to output file." << endl;
744
            exit(1);
745
        }
746

747
        // write "header" of the .out file
748
        size_t nr_orig_gens = 0;
749
        if (lattice_ideal_input) {
750
            nr_orig_gens = Result->getNrOriginalMonoidGenerators();
751
            out << nr_orig_gens <<" original generators of the toric ring"<<endl;
752
        }
753
        if (Result->isComputed(ConeProperty::ModuleGenerators)) {
754
            out << Result->getNrModuleGenerators() <<" module generators" << endl;
755
        }
756
        if (Result->isComputed(ConeProperty::HilbertBasis)) {
757
            out << Result->getNrHilbertBasis() <<" Hilbert basis elements"
758
                << of_monoid << endl;
759
        }
760
        if (homogeneous && Result->isComputed(ConeProperty::Deg1Elements)) {
761
            out << Result->getNrDeg1Elements() <<" Hilbert basis elements of degree 1"<<endl;
762
        }
763
        if (Result->isComputed(ConeProperty::IsReesPrimary)
764
            && Result->isComputed(ConeProperty::HilbertBasis)) {
765
            const Matrix<Integer>& Hilbert_Basis = Result->getHilbertBasisMatrix();
766
            nr = Hilbert_Basis.nr_of_rows();
767
            for (i = 0; i < nr; i++) {
768
                if (Hilbert_Basis[i][dim-1]==1) {
769
                    rees_ideal_key.push_back(i);
770
                 }
771
            }
772
            out << rees_ideal_key.size() <<" generators of integral closure of the ideal"<<endl;
773
        }
774
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
775
            out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron" << endl;
776
        }
777
        if (Result->isComputed(ConeProperty::ExtremeRays)) {
778
            out << Result->getNrExtremeRays() <<" extreme rays" << of_cone << endl;
779
        }
780
        if(Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {
781
            out << Result->getNrModuleGeneratorsOverOriginalMonoid() <<" module generators over original monoid" << endl;    
782
        }
783
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) {
784
            out << Result->getNrSupportHyperplanes() <<" support hyperplanes"
785
                << of_polyhedron << endl;
786
        }
787
        out<<endl;
788
        if (Result->isComputed(ConeProperty::ExcludedFaces)) {
789
            out << Result->getNrExcludedFaces() <<" excluded faces"<<endl;
790
            out << endl;
791
        }
792
        out << "embedding dimension = " << dim << endl;
793
        if (homogeneous) {
794
            if (Result->isComputed(ConeProperty::Sublattice)) {
795
                auto rank = Result->getRank();
796
                out << "rank = "<< rank << is_maximal(rank,dim) << endl;
797
                out << "external index = "<< Result->getSublattice().getExternalIndex() << endl;
798
            }
799
        } else { // now inhomogeneous case
800
            if (Result->isComputed(ConeProperty::AffineDim))
801
                out << "affine dimension of the polyhedron = "
802
                    << Result->getAffineDim() << is_maximal(Result->getAffineDim(),dim-1) << endl;
803
            if (Result->isComputed(ConeProperty::RecessionRank))
804
                out << "rank of recession monoid = "  << Result->getRecessionRank() << endl;
805
        }
806
        
807
        if(Result->isComputed(ConeProperty::OriginalMonoidGenerators)){
808
            out << "internal index = " << Result->getIndex() << endl;           
809
        }
810
        
811
        if(Result->isComputed(ConeProperty::MaximalSubspace)){
812
            size_t dim_max_subspace=Result->getDimMaximalSubspace();
813
            if(dim_max_subspace>0)
814
                out << "dimension of maximal subspace = " << dim_max_subspace << endl;      
815
        }
816
            
817
        
818
        if (homogeneous && Result->isComputed(ConeProperty::IsIntegrallyClosed)) {
819
            if (Result->isIntegrallyClosed()) {
820
                out << "original monoid is integrally closed"<<endl;
821
            } else {
822
                out << "original monoid is not integrally closed"<<endl;
823
                if ( Result->isComputed(ConeProperty::IsIntegrallyClosed)
824
                       && !Result->isComputed(ConeProperty::HilbertBasis)) {
825
                    out << "witness for not being integrally closed:" << endl;
826
                    out << Result->getWitnessNotIntegrallyClosed();
827
                }
828
                if(Result->getUnitGroupIndex()>1){
829
                    out << "unit group index = " << Result->getUnitGroupIndex() << endl;
830
                }
831
            }
832
        }
833
        out << endl;
834
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
835
            out << "size of ";
836
            if (Result->isTriangulationNested()) out << "nested ";
837
            if (Result->isTriangulationPartial()) out << "partial ";
838
            out << "triangulation   = " << Result->getTriangulationSize() << endl;
839
        }
840
        if (Result->isComputed(ConeProperty::TriangulationDetSum)) {
841
            out << "resulting sum of |det|s = " << Result->getTriangulationDetSum() << endl;
842
        }
843
        if (Result->isComputed(ConeProperty::TriangulationSize)) {
844
            out << endl;
845
        }
846
        if ( Result->isComputed(ConeProperty::Dehomogenization) ) {
847
            out << "dehomogenization:" << endl
848
                << Result->getDehomogenization() << endl;
849
        }
850
        if ( Result->isComputed(ConeProperty::Grading) ) {
851
            out << "grading:" << endl
852
                << Result->getGrading();
853
            Integer denom = Result->getGradingDenom();
854
            if (denom != 1) {
855
                out << "with denominator = " << denom << endl;
856
            }
857
            out << endl;
858
            if (homogeneous && Result->isComputed(ConeProperty::ExtremeRays)) {
859
                out << "degrees of extreme rays:"<<endl;
860
                map<Integer,long> deg_count;
861
                vector<Integer> degs = Result->getExtremeRaysMatrix().MxV(Result->getGrading());
862
                for (i=0; i<degs.size(); ++i) {
863
                    deg_count[degs[i]/denom]++;
864
                }
865
                out << deg_count;
866
            }
867
        }
868
        else if (Result->isComputed(ConeProperty::IsDeg1ExtremeRays)) {
869
            if ( !Result->isDeg1ExtremeRays() ) {
870
                out << "No implicit grading found" << endl;
871
            }
872
        }
873
        out<<endl;
874
        if (homogeneous && Result->isComputed(ConeProperty::IsDeg1HilbertBasis)
875
          && Result->isDeg1ExtremeRays() ) {
876
            if (Result->isDeg1HilbertBasis()) {
877
                out << "Hilbert basis elements are of degree 1";
878
            } else {
879
                out << "Hilbert basis elements are not of degree 1";
880
            }
881
            out<<endl<<endl;
882
        }
883
        if ( Result->isComputed(ConeProperty::ModuleRank) ) {
884
            out << "module rank = "<< Result->getModuleRank() << endl;
885
        }
886
        if ( Result->isComputed(ConeProperty::Multiplicity) ) {
887
            out << "multiplicity = "<< Result->getMultiplicity() << endl;
888
        }
889
        if ( Result->isComputed(ConeProperty::ModuleRank) || Result->isComputed(ConeProperty::Multiplicity)) {
890
            out << endl;
891
        }
892
        
893
        if ( Result->isComputed(ConeProperty::HilbertSeries) ) {
894
            const HilbertSeries& HS = Result->getHilbertSeries();
895
            vector<mpz_class> HS_Num;
896
            map<long, long> HS_Denom;
897
            if ( Result->isComputed(ConeProperty::HSOP) ){
898
                    HS_Denom=HS.getHSOPDenom();
899
                    HS_Num=HS.getHSOPNum();
900
                    out << "Hilbert series (HSOP):" << endl << HS_Num;
901
            } else {
902
                    HS_Denom=HS.getDenom();
903
                    HS_Num=HS.getNum();
904
                    out << "Hilbert series:" << endl << HS_Num;
905
            }
906
            long nr_factors = 0;
907
            for (map<long, long>::iterator it = HS_Denom.begin(); it!=HS_Denom.end(); ++it) {
908
                nr_factors += it->second;
909
            }
910
            out << "denominator with " << nr_factors << " factors:" << endl;
911
            out << HS_Denom;
912
            out << endl;
913
            if (HS.getShift() != 0) {
914
                out << "shift = " << HS.getShift() << endl << endl;
915
            }
916
            
917
            out << "degree of Hilbert Series as rational function = "
918
                << HS.getDegreeAsRationalFunction() << endl << endl;
919
            if(v_is_symmetric(HS_Num)){
920
                out << "The numerator of the Hilbert Series is symmetric." << endl << endl;
921
            }
922
            long period = HS.getPeriod();
923
            if (period == 1 && (HS_Denom.size() == 0
924
                                || HS_Denom.begin()->first== (long) HS_Denom.size())) {
925
                out << "Hilbert polynomial:" << endl;
926
                out << HS.getHilbertQuasiPolynomial()[0];
927
                out << "with common denominator = ";
928
                out << HS.getHilbertQuasiPolynomialDenom();
929
                out << endl<< endl;
930
            } else {
931
                // output cyclonomic representation
932
                out << "Hilbert series with cyclotomic denominator:" << endl;
933
                out << HS.getCyclotomicNum();
934
                out << "cyclotomic denominator:" << endl;
935
                out << HS.getCyclotomicDenom();
936
                out << endl;
937
                // Hilbert quasi-polynomial
938
                HS.computeHilbertQuasiPolynomial();
939
                if (HS.isHilbertQuasiPolynomialComputed()) {
940
                    out<<"Hilbert quasi-polynomial of period " << period << ":" << endl;
941
                    if(HS.get_nr_coeff_quasipol()>=0){
942
                        out << "only " << HS.get_nr_coeff_quasipol() << " highest coefficients computed" << endl;
943
                        out << "their common period is " << HS.getHilbertQuasiPolynomial().size() << "" << endl;                        
944
                    }
945
                    Matrix<mpz_class> HQP(HS.getHilbertQuasiPolynomial());
946
                    HQP.pretty_print(out,true);
947
                    out<<"with common denominator = "<<HS.getHilbertQuasiPolynomialDenom();
948
                }else{
949
                    out<<"Hilbert quasi-polynomial has period " << period << endl;    
950
                }
951
                out << endl << endl;
952
            }
953

954
        }
955
 
956
        if ( Result->isComputed(ConeProperty::WeightedEhrhartSeries) )
957
            writeWeightedEhrhartSeries(out);
958
            
959
        
960
        if ( Result->isComputed(ConeProperty::VirtualMultiplicity)
961
                   && !Result->isComputed(ConeProperty::WeightedEhrhartQuasiPolynomial)) {
962
            out << "virtual multiplicity of weighted Ehrhart series = "<< Result->getVirtualMultiplicity() << endl;
963
            out << endl;
964
        }
965
        
966
        if ( Result->isComputed(ConeProperty::Integral)) {
967
            out << "integral  = "<< Result->getIntegral() << endl;
968
            out << endl;
969
        }
970

971
        if (Result->isComputed(ConeProperty::IsReesPrimary)) {
972
            if (Result->isReesPrimary()) {
973
                out<<"ideal is primary to the ideal generated by the indeterminates"<<endl;
974
            } else {
975
                out<<"ideal is not primary to the ideal generated by the indeterminates"<<endl;
976
            }
977
            if (Result->isComputed(ConeProperty::ReesPrimaryMultiplicity)) {
978
                out<<"multiplicity of the ideal = "<<Result->getReesPrimaryMultiplicity()<<endl;
979
            }
980
            out << endl;
981
        }
982
        
983
        if(Result->isComputed(ConeProperty::ClassGroup)) {
984
            vector<Integer> ClassGroup=Result->getClassGroup();
985
            out << "rank of class group = " << ClassGroup[0] << endl;
986
            if(ClassGroup.size()==1)
987
                out << "class group is free" << endl << endl;
988
            else{
989
                ClassGroup.erase(ClassGroup.begin());
990
                out << "finite cyclic summands:" << endl;
991
                out << count_in_map<Integer,size_t>(ClassGroup);
992
                out << endl;
993
            }                   
994
        }
995
        
996
        if(Result->isComputed(ConeProperty::IsGorenstein)) {
997
            if(Result->isGorenstein()){
998
                out << "Monoid is Gorenstein " << endl;
999
                out << "Generator of interior:" << endl;
1000
                out << Result->getGeneratorOfInterior();              
1001
            }
1002
            else
1003
                out << "Monoid is not Gorenstein " << endl;
1004
            out << endl;
1005
        }
1006

1007
        out << "***********************************************************************"
1008
            << endl << endl;
1009

1010

1011
        if (lattice_ideal_input) {
1012
            out << nr_orig_gens <<" original generators:"<<endl;
1013
            Result->getOriginalMonoidGeneratorsMatrix().pretty_print(out);
1014
            out << endl;
1015
        }
1016
        if (Result->isComputed(ConeProperty::ModuleGenerators)) {
1017
            out << Result->getNrModuleGenerators() <<" module generators:" << endl;
1018
            Result->getModuleGeneratorsMatrix().pretty_print(out);
1019
            out << endl;
1020
        }
1021
        
1022
        if ( Result->isComputed(ConeProperty::Deg1Elements) ) {
1023
            const Matrix<Integer>& Hom = Result->getDeg1ElementsMatrix();
1024
            write_matrix_ht1(Hom);
1025
            nr=Hom.nr_of_rows();
1026
            out<<nr<<" Hilbert basis elements of degree 1:"<<endl;
1027
            Hom.pretty_print(out);
1028
            out << endl;
1029
        }
1030
        
1031
        if (Result->isComputed(ConeProperty::HilbertBasis)) {
1032

1033
            const Matrix<Integer>& Hilbert_Basis = Result->getHilbertBasisMatrix();
1034
            
1035
            if(!Result->isComputed(ConeProperty::Deg1Elements)){
1036
                nr=Hilbert_Basis.nr_of_rows();
1037
                out << nr << " Hilbert basis elements" << of_monoid << ":" << endl;
1038
                Hilbert_Basis.pretty_print(out);
1039
                out << endl;
1040
            }
1041
            else {
1042
                nr=Hilbert_Basis.nr_of_rows()-Result->getNrDeg1Elements();
1043
                out << nr << " further Hilbert basis elements" << of_monoid << " of higher degree:" << endl;
1044
                Matrix<Integer> HighDeg(nr,dim);
1045
                for(size_t i=0;i<nr;++i)
1046
                    HighDeg[i]=Hilbert_Basis[i+Result->getNrDeg1Elements()];                               
1047
                HighDeg.pretty_print(out);
1048
                out << endl;                       
1049
            }
1050
            Matrix<Integer> complete_Hilbert_Basis(0,dim);
1051
            if (gen || egn || typ) {
1052
                // for these files we append the module generators if there are any
1053
                if (Result->isComputed(ConeProperty::ModuleGenerators)) {
1054
                    complete_Hilbert_Basis.append(Hilbert_Basis);
1055
                    complete_Hilbert_Basis.append(Result->getModuleGeneratorsMatrix());
1056
                    write_matrix_gen(complete_Hilbert_Basis);
1057
                } else {
1058
                    write_matrix_gen(Hilbert_Basis);
1059
                }                
1060
            }
1061
            if ((egn || typ) && Result->isComputed(ConeProperty::Sublattice)) {
1062
                const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();
1063
                Matrix<Integer> Hilbert_Basis_Full_Cone = BasisChange.to_sublattice(Hilbert_Basis);
1064
                if (Result->isComputed(ConeProperty::ModuleGenerators)) {
1065
                    Hilbert_Basis_Full_Cone.append(BasisChange.to_sublattice(Result->getModuleGeneratorsMatrix()));
1066
                }
1067
                if (egn) {
1068
                    string egn_string = name+".egn";
1069
                    const char* egn_file = egn_string.c_str();
1070
                    ofstream egn_out(egn_file);        
1071
                    Hilbert_Basis_Full_Cone.print(egn_out);
1072
                    // egn_out<<"cone"<<endl;
1073
                    egn_out.close();
1074
                }    
1075

1076
                if (typ && homogeneous) {
1077
                    write_matrix_typ(Hilbert_Basis_Full_Cone.multiplication(BasisChange.to_sublattice_dual(Result->getSupportHyperplanesMatrix()).transpose()));
1078
                }
1079
            }
1080

1081
            if (Result->isComputed(ConeProperty::IsReesPrimary)) {
1082
                out << rees_ideal_key.size() <<" generators of integral closure of the ideal:"<<endl;
1083
                Matrix<Integer> Ideal_Gens = Hilbert_Basis.submatrix(rees_ideal_key);
1084
                Ideal_Gens.resize_columns(dim-1);
1085
                Ideal_Gens.pretty_print(out);
1086
                out << endl;
1087
            }
1088
        }
1089
        if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)  && !no_ext_rays_output) {
1090
            out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron:" << endl;
1091
            if(Result->isComputed(ConeProperty::VerticesFloat))
1092
                write_float(out,Result->getVerticesFloatMatrix(),Result->getNrVerticesFloat(),dim);
1093
            else
1094
                Result->getVerticesOfPolyhedronMatrix().pretty_print(out);
1095
            out << endl;
1096
        }
1097
        if (Result->isComputed(ConeProperty::ExtremeRays) && !no_ext_rays_output) {
1098
            out << Result->getNrExtremeRays() << " extreme rays" << of_cone << ":" << endl;
1099
            if(homogeneous && Result->isComputed(ConeProperty::VerticesFloat))
1100
                write_float(out,Result->getVerticesFloatMatrix(),Result->getNrVerticesFloat(),dim); 
1101
            else
1102
                Result->getExtremeRaysMatrix().pretty_print(out);
1103
            out << endl;
1104
        }
1105
        
1106
        if (Result->isComputed(ConeProperty::ExtremeRays) && ext){
1107
            // for the .gen file we append the vertices of polyhedron if there are any
1108
            if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {
1109
                Matrix<Integer> Extreme_Rays(Result->getExtremeRaysMatrix());
1110
                Extreme_Rays.append(Result->getVerticesOfPolyhedronMatrix());
1111
                write_matrix_ext(Extreme_Rays);
1112
            } else {
1113
                write_matrix_ext(Result->getExtremeRaysMatrix());
1114
            }
1115
        }
1116
        
1117
        if(Result->isComputed(ConeProperty::MaximalSubspace) && Result->getDimMaximalSubspace()>0){
1118
            out << Result->getDimMaximalSubspace() <<" basis elements of maximal subspace:" << endl;
1119
            Result->getMaximalSubspaceMatrix().pretty_print(out);
1120
            out << endl;
1121
            if(msp)
1122
                write_matrix_msp(Result->getMaximalSubspaceMatrix());            
1123
        }
1124
        
1125
        if(Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {
1126
            out << Result->getNrModuleGeneratorsOverOriginalMonoid() <<" module generators over original monoid:" << endl;
1127
            Result->getModuleGeneratorsOverOriginalMonoidMatrix().pretty_print(out);
1128
            out << endl;
1129
            if(mod)
1130
                write_matrix_mod(Result->getModuleGeneratorsOverOriginalMonoidMatrix());
1131
        }
1132

1133
        //write constrains (support hyperplanes, congruences, equations)
1134

1135
        if (Result->isComputed(ConeProperty::SupportHyperplanes)) {
1136
            const Matrix<Integer>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();
1137
            out << Support_Hyperplanes.nr_of_rows() <<" support hyperplanes" 
1138
                << of_polyhedron << ":" << endl;
1139
            Support_Hyperplanes.pretty_print(out);
1140
            out << endl;
1141
        }
1142
        if (Result->isComputed(ConeProperty::Sublattice)) {
1143
            const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();
1144
            //equations
1145
            const Matrix<Integer>& Equations = BasisChange.getEquationsMatrix();
1146
            size_t nr_of_equ = Equations.nr_of_rows();
1147
            if (nr_of_equ > 0) {
1148
                out << nr_of_equ <<" equations:" <<endl;
1149
                Equations.pretty_print(out);
1150
                out << endl;
1151
            }
1152

1153
            //congruences
1154
            const Matrix<Integer>& Congruences = BasisChange.getCongruencesMatrix();
1155
            size_t nr_of_cong = Congruences.nr_of_rows();
1156
            if (nr_of_cong > 0) {
1157
                out << nr_of_cong <<" congruences:" <<endl;
1158
                Congruences.pretty_print(out);
1159
                out << endl;
1160
            }
1161
            
1162
            //lattice
1163
            const Matrix<Integer>& LatticeBasis = BasisChange.getEmbeddingMatrix();
1164
            size_t nr_of_latt = LatticeBasis.nr_of_rows();
1165
            if (nr_of_latt < dim ||  BasisChange.getExternalIndex()!=1) {
1166
                out << nr_of_latt <<" basis elements of lattice:" <<endl;
1167
                LatticeBasis.pretty_print(out);
1168
                out << endl;
1169
            }
1170
            if(lat)
1171
                write_matrix_lat(LatticeBasis);
1172
            
1173

1174
            //excluded faces
1175
            if (Result->isComputed(ConeProperty::ExcludedFaces)) {
1176
                const Matrix<Integer>& ExFaces = Result->getExcludedFacesMatrix();
1177
                out << ExFaces.nr_of_rows() <<" excluded faces:" <<endl;
1178
                ExFaces.pretty_print(out);
1179
                out << endl;
1180
            }
1181

1182
            if (cst && Result->isComputed(ConeProperty::SupportHyperplanes)) {
1183
                const Matrix<Integer>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();
1184
                string cst_string = name+".cst";
1185
                const char* cst_file = cst_string.c_str();
1186
                ofstream cst_out(cst_file);
1187
    
1188
                Support_Hyperplanes.print(cst_out);
1189
                cst_out<<"inequalities"<<endl;
1190
                Equations.print(cst_out);
1191
                cst_out<<"equations"<<endl;
1192
                Congruences.print(cst_out);
1193
                cst_out<<"congruences"<<endl;
1194
                if (Result->isComputed(ConeProperty::ExcludedFaces)) {
1195
                    Result->getExcludedFacesMatrix().print(cst_out);
1196
                    cst_out<<"excluded_faces"<<endl;
1197
                }
1198
                if (Result->isComputed(ConeProperty::Grading)) {
1199
                    cst_out << 1 << endl << dim << endl;
1200
                    cst_out << Result->getGrading();
1201
                    cst_out << "grading" << endl;
1202
                }
1203
                if (Result->isComputed(ConeProperty::Dehomogenization)) {
1204
                    cst_out << 1 << endl << dim << endl;
1205
                    cst_out << Result->getDehomogenization();
1206
                    cst_out << "dehomogenization" << endl;
1207
                }
1208
                cst_out.close();
1209
            }
1210
        }
1211
        
1212
        out.close();
1213
    }
1214

1215
    write_inv_file();
1216
    write_Stanley_dec();
1217
}
1218

1219

1220
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.
