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.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418384/*1* Normaliz2* Copyright (C) 2007-2014 Winfried Bruns, Bogdan Ichim, Christof Soeger3* This program is free software: you can redistribute it and/or modify4* it under the terms of the GNU General Public License as published by5* the Free Software Foundation, either version 3 of the License, or6* (at your option) any later version.7*8* This program is distributed in the hope that it will be useful,9* but WITHOUT ANY WARRANTY; without even the implied warranty of10* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the11* GNU General Public License for more details.12*13* You should have received a copy of the GNU General Public License14* along with this program. If not, see <http://www.gnu.org/licenses/>.15*16* As an exception, when this program is distributed through (i) the App Store17* by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play18* by Google Inc., then that store may impose any digital rights management,19* device limits and/or redistribution restrictions that are required by its20* terms of service.21*/2223//---------------------------------------------------------------------------2425#include <stdlib.h>26#include <vector>27#include <list>28#include <fstream>29#include <iostream>30#include <string>31#include <algorithm>32#include <iomanip>3334#include "output.h"35#include "libnormaliz/general.h"36#include "libnormaliz/matrix.h"37#include "libnormaliz/vector_operations.h"38#include "libnormaliz/map_operations.h"3940//---------------------------------------------------------------------------4142template<typename Integer>43Output<Integer>::Output(){44out=true;45inv=false;46ext=false;47esp=false;48typ=false;49egn=false;50gen=false;51cst=false;52tri=false;53tgn=false;54ht1=false;55dec=false;56lat=false;57mod=false;58msp=false;59lattice_ideal_input = false;60no_ext_rays_output=false;61}6263//---------------------------------------------------------------------------6465template<typename Integer>66void Output<Integer>::set_lattice_ideal_input(bool value){67lattice_ideal_input=value;68}6970//---------------------------------------------------------------------------7172template<typename Integer>73void Output<Integer>::set_no_ext_rays_output(){74no_ext_rays_output=true;75}7677//---------------------------------------------------------------------------7879template<typename Integer>80void Output<Integer>::read() const{81cout<<"\nname="<<name<<"\n";82cout<<"\nout="<<out<<"\n";83cout<<"\ninv="<<inv<<"\n";84cout<<"\next="<<ext<<"\n";85cout<<"\nesp="<<esp<<"\n";86cout<<"\ntyp="<<typ<<"\n";87cout<<"\negn="<<egn<<"\n";88cout<<"\ngen="<<gen<<"\n";89cout<<"\ncst="<<cst<<"\n";90cout<<"\ntri="<<tri<<"\n";91cout<<"\ntgn="<<tgn<<"\n";92cout<<"\nht1="<<ht1<<"\n";93cout<<"\nResult is:\n";94Result->print();95}9697//---------------------------------------------------------------------------9899template<typename Integer>100void Output<Integer>::set_name(const string& n){101name=n;102}103104//---------------------------------------------------------------------------105106template<typename Integer>107void Output<Integer>::setCone(Cone<Integer> & C) {108this->Result = &C;109dim = Result->getEmbeddingDim();110homogeneous = !Result->isInhomogeneous();111if (homogeneous) {112of_cone = "";113of_monoid = "";114of_polyhedron = "";115} else {116of_cone = " of recession cone";117of_monoid = " of recession monoid";118of_polyhedron = " of polyhedron (homogenized)";119}120}121122//---------------------------------------------------------------------------123124template<typename Integer>125void Output<Integer>::set_write_out(const bool& flag){126out=flag;127}128129//---------------------------------------------------------------------------130131template<typename Integer>132void Output<Integer>::set_write_inv(const bool& flag){133inv=flag;134}135136//---------------------------------------------------------------------------137138template<typename Integer>139void Output<Integer>::set_write_ext(const bool& flag){140ext=flag;141}142143//---------------------------------------------------------------------------144145template<typename Integer>146void Output<Integer>::set_write_esp(const bool& flag){147esp=flag;148}149150//---------------------------------------------------------------------------151152template<typename Integer>153void Output<Integer>::set_write_typ(const bool& flag){154typ=flag;155}156157//---------------------------------------------------------------------------158159template<typename Integer>160void Output<Integer>::set_write_egn(const bool& flag){161egn=flag;162}163164//---------------------------------------------------------------------------165166template<typename Integer>167void Output<Integer>::set_write_gen(const bool& flag){168gen=flag;169}170171//---------------------------------------------------------------------------172173template<typename Integer>174void Output<Integer>::set_write_cst(const bool& flag){175cst=flag;176}177178//---------------------------------------------------------------------------179180template<typename Integer>181void Output<Integer>::set_write_tri(const bool& flag) {182tri=flag;183}184185//---------------------------------------------------------------------------186187template<typename Integer>188void Output<Integer>::set_write_tgn(const bool& flag) {189tgn=flag;190}191192//---------------------------------------------------------------------------193194template<typename Integer>195void Output<Integer>::set_write_ht1(const bool& flag) {196ht1=flag;197}198199//---------------------------------------------------------------------------200201template<typename Integer>202void Output<Integer>::set_write_dec(const bool& flag) {203dec=flag;204}205206//---------------------------------------------------------------------------207208template<typename Integer>209void Output<Integer>::set_write_mod(const bool& flag) {210mod=flag;211}212213//---------------------------------------------------------------------------214215template<typename Integer>216void Output<Integer>::set_write_lat(const bool& flag) {217lat=flag;218}219220221//---------------------------------------------------------------------------222223template<typename Integer>224void Output<Integer>::set_write_msp(const bool& flag) {225msp=flag;226}227//---------------------------------------------------------------------------228229template<typename Integer>230void Output<Integer>::set_write_extra_files(){231out=true;232inv=true;233gen=true;234cst=true;235}236237//---------------------------------------------------------------------------238239template<typename Integer>240void Output<Integer>::set_write_all_files(){241out=true;242inv=true;243ext=true;244esp=true;245typ=true;246egn=true;247gen=true;248cst=true;249ht1=true;250lat=true;251mod=true;252msp=true;253}254255256257//---------------------------------------------------------------------------258259template<typename Integer>260void Output<Integer>::write_matrix_ext(const Matrix<Integer>& M) const{261if (ext==true) {262M.print(name,"ext");263}264}265266//---------------------------------------------------------------------------267268template<typename Integer>269void Output<Integer>::write_matrix_mod(const Matrix<Integer>& M) const{270if (mod==true) {271M.print(name,"mod");272}273}274275276//---------------------------------------------------------------------------277278template<typename Integer>279void Output<Integer>::write_matrix_lat(const Matrix<Integer>& M) const{280if (ext==true) {281M.print(name,"lat");282}283}284285//---------------------------------------------------------------------------286287template<typename Integer>288void Output<Integer>::write_matrix_esp(const Matrix<Integer>& M) const{289if (esp==true) {290M.print(name,"esp");291}292}293294//---------------------------------------------------------------------------295296template<typename Integer>297void Output<Integer>::write_matrix_typ(const Matrix<Integer>& M) const{298if (typ==true) {299M.print(name,"typ");300}301}302303//---------------------------------------------------------------------------304305template<typename Integer>306void Output<Integer>::write_matrix_egn(const Matrix<Integer>& M) const {307if (egn==true) {308M.print(name,"egn");309}310}311312//---------------------------------------------------------------------------313314template<typename Integer>315void Output<Integer>::write_matrix_gen(const Matrix<Integer>& M) const {316if (gen==true) {317M.print(name,"gen");318}319}320//---------------------------------------------------------------------------321322template<typename Integer>323void Output<Integer>::write_matrix_msp(const Matrix<Integer>& M) const {324if (msp==true) {325M.print(name,"msp");326}327}328329//---------------------------------------------------------------------------330331template<typename Integer>332void Output<Integer>::write_tri() const{333if (tri==true) {334string file_name = name+".tri";335ofstream out(file_name.c_str());336337const vector< pair<vector<libnormaliz::key_t>,Integer> >& Tri = Result->getTriangulation();338typename vector< pair<vector<libnormaliz::key_t>,Integer> >::const_iterator tit = Tri.begin();339const vector<vector<bool> >& Dec = Result->isComputed(ConeProperty::ConeDecomposition) ?340Result->getOpenFacets() : vector<vector<bool> >();341typename vector< vector<bool> >::const_iterator idd = Dec.begin();342343out << Tri.size() << endl;344size_t nr_extra_entries=1;345if (Result->isComputed(ConeProperty::ConeDecomposition))346nr_extra_entries+=Result->getSublattice().getRank();347out << Result->getSublattice().getRank()+nr_extra_entries << endl; //works also for empty list348349for(; tit != Tri.end(); ++tit) {350for (size_t i=0; i<tit->first.size(); i++) {351out << tit->first[i] +1 << " ";352}353out << " " << tit->second;354if(Result->isComputed(ConeProperty::ConeDecomposition)){355out << " ";356for (size_t i=0; i<tit->first.size(); i++) {357out << " " << (*idd)[i];358}359idd++;360}361out << endl;362}363if (Result->isTriangulationNested()) out << "nested" << endl;364else out << "plain" << endl;365if (Result->isTriangulationPartial()) out << "partial" << endl;366out.close();367}368}369370//---------------------------------------------------------------------------371372373template<typename Integer>374void Output<Integer>::write_Stanley_dec() const {375if (dec && Result->isComputed(ConeProperty::StanleyDec)) {376ofstream out((name+".dec").c_str());377378if (Result->isComputed(ConeProperty::InclusionExclusionData)) {379const vector< pair<vector<libnormaliz::key_t>, long> >& InExData = Result->getInclusionExclusionData();380out << "in_ex_data" << endl;381out << InExData.size() << endl;382for (size_t i=0; i<InExData.size(); ++i) {383out << InExData[i].first.size() << " ";384for (size_t j=0; j<InExData[i].first.size(); ++j) {385out << InExData[i].first[j]+1 << " ";386}387out << InExData[i].second << endl;388}389}390391out << "Stanley_dec" << endl;392list<STANLEYDATA_int>& StanleyDec = Result->getStanleyDec_mutable();393auto S = StanleyDec.begin();394size_t i;395396out << StanleyDec.size() << endl;397for (; S!=StanleyDec.end(); ++S) {398for (i=0; i < S->key.size(); ++i)399out << S->key[i]+1 <<" ";400out << endl;401S->offsets.print(out);402out << endl;403}404out.close();405}406}407408//---------------------------------------------------------------------------409410template<typename Integer>411void Output<Integer>::write_matrix_ht1(const Matrix<Integer>& M) const{412if (ht1==true) {413M.print(name,"ht1");414}415}416417//---------------------------------------------------------------------------418419string is_maximal(long a, long b) {420return (a == b) ? " (maximal)" : "";421}422423//---------------------------------------------------------------------------424425template<typename Integer>426void Output<Integer>::write_inv_file() const{427if (inv==true) {//printing .inv file428size_t i;429string name_open=name+".inv"; //preparing output files430const char* file=name_open.c_str();431ofstream inv(file);432433if (Result->isComputed(ConeProperty::ModuleGenerators)) {434inv << "integer number_module_generators = "435<< Result->getNrModuleGenerators() << endl;436}437if (Result->isComputed(ConeProperty::HilbertBasis)) {438inv<<"integer hilbert_basis_elements = "<<Result->getNrHilbertBasis()<<endl;439}440441if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {442inv << "integer number_vertices_polyhedron = "443<< Result->getNrVerticesOfPolyhedron() << endl;444}445if (Result->isComputed(ConeProperty::ExtremeRays)) {446size_t nr_ex_rays = Result->getNrExtremeRays();447inv<<"integer number_extreme_rays = "<<nr_ex_rays<<endl;448}449if (Result->isComputed(ConeProperty::MaximalSubspace)) {450size_t dim_max_subspace = Result->getDimMaximalSubspace();451inv<<"integer dim_max_subspace = "<<dim_max_subspace<<endl;452}453if (Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {454inv << "integer number_module_generators_original_monoid = "455<< Result->getNrModuleGeneratorsOverOriginalMonoid() << endl;456}457458inv << "integer embedding_dim = " << dim << endl;459if (homogeneous) {460if (Result->isComputed(ConeProperty::Sublattice)) {461inv << "integer rank = " << Result->getRank() << endl;462inv << "integer external_index = " << Result->getSublattice().getExternalIndex() << endl;463}464} else {465if (Result->isComputed(ConeProperty::AffineDim))466inv << "integer affine_dim_polyhedron = " << Result->getAffineDim() << endl;467if (Result->isComputed(ConeProperty::RecessionRank))468inv << "integer recession_rank = " << Result->getRecessionRank() << endl;469}470471if(Result->isComputed(ConeProperty::OriginalMonoidGenerators)){472inv << "integer internal_index = " << Result->getIndex() << endl;473}474if (Result->isComputed(ConeProperty::SupportHyperplanes)) {475inv<<"integer number_support_hyperplanes = "<<Result->getNrSupportHyperplanes()<<endl;476}477if (Result->isComputed(ConeProperty::TriangulationSize)) {478inv << "integer size_triangulation = " << Result->getTriangulationSize() << endl;479}480if (Result->isComputed(ConeProperty::TriangulationDetSum)) {481inv << "integer sum_dets = " << Result->getTriangulationDetSum() << endl;482}483484if (Result->isComputed(ConeProperty::IsIntegrallyClosed)) {485if(Result->isIntegrallyClosed())486inv << "boolean integrally_closed = true" << endl;487else488inv << "boolean integrally_closed = false" << endl;489}490491if (!Result->isComputed(ConeProperty::Dehomogenization)) {492inv << "boolean inhomogeneous = false" << endl;493}494else {495inv << "boolean inhomogeneous = true" << endl;496vector<Integer> Linear_Form = Result->getDehomogenization();497inv << "vector " << Linear_Form.size()498<< " dehomogenization = " << Linear_Form;499}500if (Result->isComputed(ConeProperty::Grading)==false) {501if (Result->isComputed(ConeProperty::ExtremeRays)) {502inv<<"boolean graded = "<<"false"<<endl;503}504}505else {506inv<<"boolean graded = "<<"true"<<endl;507if (Result->isComputed(ConeProperty::Deg1Elements)) {508inv<<"integer degree_1_elements = "<<Result->getNrDeg1Elements()<<endl;509}510vector<Integer> Linear_Form = Result->getGrading();511inv<<"vector "<<Linear_Form.size()<<" grading = ";512for (i = 0; i < Linear_Form.size(); i++) {513inv<<Linear_Form[i]<<" ";514}515inv<<endl;516inv <<"integer grading_denom = " << Result->getGradingDenom() << endl;517if (Result->isComputed(ConeProperty::Multiplicity)){518mpq_class mult = Result->getMultiplicity();519inv <<"integer multiplicity = " << mult.get_num() << endl;520inv <<"integer multiplicity_denom = "<< mult.get_den() << endl;521}522if (Result->isComputed(ConeProperty::WeightedEhrhartSeries)){523const HilbertSeries& HS=Result->getIntData().getWeightedEhrhartSeries().first;524525inv <<"vector "<< HS.getNum().size() <<" weighted_ehrhart_series_num = " << HS.getNum();526inv <<"integer num_common_denom = " << Result->getIntData().getWeightedEhrhartSeries().second << endl;527inv <<"vector "<< to_vector(HS.getDenom()).size() <<" weighted_ehrhart_series_num = "528<< to_vector(HS.getDenom());529Result->getIntData().computeWeightedEhrhartQuasiPolynomial();530if (Result->getIntData().isWeightedEhrhartQuasiPolynomialComputed()) {531if(HS.get_nr_coeff_quasipol()>=0){532inv << "integer nr_coeff_weightedEhrhart__quasipol = " << HS.get_nr_coeff_quasipol() << endl;533}534vector< vector <mpz_class> > hqp= Result->getIntData().getWeightedEhrhartQuasiPolynomial();535inv << "matrix " << hqp.size() << " " << hqp[0].size()536<< " weighted_ehrhart_quasipolynomial = ";537inv << endl << hqp;538inv << "integer weighted_ehrhart_quasipolynomial_denom = "539<< Result->getIntData().getWeightedEhrhartQuasiPolynomialDenom() << endl;540}541}542543if (Result->isComputed(ConeProperty::VirtualMultiplicity)){544mpq_class mult = Result->getVirtualMultiplicity();545inv <<"integer virtual_multiplicity = " << mult.get_num() << endl;546inv <<"integer virtual_multiplicity_denom = "<< mult.get_den() << endl;547}548if (Result->isComputed(ConeProperty::Integral)){549mpq_class mult = Result->getIntegral();550inv <<"integer integral = " << mult.get_num() << endl;551inv <<"integer integral_denom = "<< mult.get_den() << endl;552}553if (Result->isComputed(ConeProperty::HilbertSeries)) {554const HilbertSeries& HS = Result->getHilbertSeries();555vector<mpz_class> HSnum;556vector<denom_t> HSdenom;557if(Result->isComputed(ConeProperty::HSOP)){558HSnum = HS.getHSOPNum();559HSdenom = to_vector(HS.getHSOPDenom());560561} else {562HSnum = HS.getNum();563HSdenom = to_vector(HS.getDenom());564}565inv <<"vector "<< HSnum.size() <<" hilbert_series_num = ";566inv << HSnum;567inv <<"vector "<< HSdenom.size() <<" hilbert_series_denom = ";568inv << HSdenom;569HS.computeHilbertQuasiPolynomial();570if (HS.isHilbertQuasiPolynomialComputed()) {571if(HS.get_nr_coeff_quasipol()>=0){572inv << "integer nr_coeff_hilbert_quasipol = " << HS.get_nr_coeff_quasipol() << endl;573}574vector< vector <mpz_class> > hqp = HS.getHilbertQuasiPolynomial();575inv << "matrix " << hqp.size() << " " << hqp[0].size()576<< " hilbert_quasipolynomial = ";577inv << endl << hqp;578inv << "integer hilbert_quasipolynomial_denom = "579<< HS.getHilbertQuasiPolynomialDenom() << endl;580}581}582}583if (Result->isComputed(ConeProperty::IsReesPrimary)) {584if (Result->isReesPrimary()) {585inv<<"boolean primary = true"<<endl;586} else {587inv<<"boolean primary = false"<<endl;588}589}590if (Result->isComputed(ConeProperty::ReesPrimaryMultiplicity)) {591inv<<"integer ideal_multiplicity = "<<Result->getReesPrimaryMultiplicity()<<endl;592}593594if (Result->isComputed(ConeProperty::ClassGroup)) {595inv <<"vector "<< Result->getClassGroup().size() <<" class_group = "<< Result->getClassGroup();596}597598if (Result->isComputed(ConeProperty::IsGorenstein)) {599if(Result->isGorenstein()){600inv << "boolean Gorenstein = true" << endl;601inv <<"vector "<< Result->getGeneratorOfInterior().size() <<" generator_of_interior = "<< Result->getGeneratorOfInterior();602}603else604inv << "boolean Gorenstein = false" << endl;605}606607inv.close();608}609}610611//---------------------------------------------------------------------------612613template<typename Integer>614void Output<Integer>::writeWeightedEhrhartSeries(ofstream& out) const{615616HilbertSeries HS=Result->getIntData().getWeightedEhrhartSeries().first;617out << "Weighted Ehrhart series:" << endl;618vector<mpz_class> num( HS.getNum());619for(size_t i=0;i<num.size();++i)620out << num[i] << " ";621out << endl << "Common denominator of coefficients: ";622out << Result->getIntData().getWeightedEhrhartSeries().second << endl;623map<long, long> HS_Denom = HS.getDenom();624long nr_factors = 0;625for (map<long, long>::iterator it = HS_Denom.begin(); it!=HS_Denom.end(); ++it) {626nr_factors += it->second;627}628out << "Series denominator with " << nr_factors << " factors:" << endl;629out << HS.getDenom();630if (HS.getShift() != 0) {631out << "shift = " << HS.getShift() << endl << endl;632}633out << endl;634out << "degree of weighted Ehrhart series as rational function = "635<< HS.getDegreeAsRationalFunction() << endl << endl;636long period = HS.getPeriod();637if (period == 1) {638out << "Weighted Ehrhart polynomial:" << endl;639for(size_t i=0; i<HS.getHilbertQuasiPolynomial()[0].size();++i)640out << HS.getHilbertQuasiPolynomial()[0][i] << " ";641out << endl;642out << "with common denominator: ";643out << HS.getHilbertQuasiPolynomialDenom()*Result->getIntData().getNumeratorCommonDenom();644} else {645// output cyclonomic representation646out << "Weighted Ehrhart series with cyclotomic denominator:" << endl;647num=HS.getCyclotomicNum();648for(size_t i=0;i<num.size();++i)649out << num[i] << " ";650out << endl << "Common denominator of coefficients: ";651out << Result->getIntData().getWeightedEhrhartSeries().second << endl;652out << "Series cyclotomic denominator:" << endl;653out << HS.getCyclotomicDenom();654out << endl;655// Weighted Ehrhart quasi-polynomial656// vector< vector<mpz_class> > hilbert_quasi_poly = HS.getHilbertQuasiPolynomial();657if (HS.isHilbertQuasiPolynomialComputed()) {658out<<"Weighted Ehrhart quasi-polynomial of period " << period << ":" << endl;659if(HS.get_nr_coeff_quasipol()>=0){660out << "only " << HS.get_nr_coeff_quasipol() << " highest coefficients computed" << endl;661out << "their common period is " << HS.getHilbertQuasiPolynomial().size() << "." << endl;662}663Matrix<mpz_class> HQP(HS.getHilbertQuasiPolynomial());664HQP.pretty_print(out,true);665out<<"with common denominator: "666<<Result->getIntData().getWeightedEhrhartQuasiPolynomialDenom() << endl;667}668else{669out<<"Weighted Ehrhart quasi-polynomial has period " << period << endl;670}671}672673out << endl << endl;674675if (HS.isHilbertQuasiPolynomialComputed()){676long deg=HS.getHilbertQuasiPolynomial()[0].size()-1;677out << "Degree of (quasi)polynomial: " << deg << endl;678679long virtDeg=Result->getRank()+Result->getIntData().getDegreeOfPolynomial()-1;680681out << endl << "Expected degree: " << virtDeg << endl;682}683684if(Result->isComputed(ConeProperty::VirtualMultiplicity)){685out << endl << "Virtual multiplicity: ";686out << Result->getIntData().getVirtualMultiplicity() << endl << endl;687}688}689690//---------------------------------------------------------------------------691template<typename Integer>692void Output<Integer>::write_float(ofstream& out, const Matrix<nmz_float>& mat, size_t nr, size_t nc) const{693694for(size_t i=0;i<nr;++i){695for(size_t j=0; j< nc;++j){696out << std::setw(10) << mat[i][j];697}698out << endl;699}700}701702//---------------------------------------------------------------------------703704template<typename Integer>705void Output<Integer>::write_files() const {706size_t i, nr;707vector<libnormaliz::key_t> rees_ideal_key;708709if (esp && Result->isComputed(ConeProperty::SupportHyperplanes) && Result->isComputed(ConeProperty::Sublattice)) {710//write the suport hyperplanes of the full dimensional cone711const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();712Matrix<Integer> Support_Hyperplanes_Full_Cone = BasisChange.to_sublattice_dual(Result->getSupportHyperplanesMatrix());713// Support_Hyperplanes_Full_Cone.print(name,"esp");714string esp_string = name+".esp";715const char* esp_file = esp_string.c_str();716ofstream esp_out(esp_file);717Support_Hyperplanes_Full_Cone.print(esp_out);718esp_out << "inequalities" << endl;719if (Result->isComputed(ConeProperty::Grading)) {720esp_out << 1 << endl << Result->getRank() << endl;721esp_out << BasisChange.to_sublattice_dual(Result->getGrading());722esp_out << "grading" << endl;723}724if (Result->isComputed(ConeProperty::Dehomogenization)) {725esp_out << 1 << endl << Result->getRank() << endl;726esp_out << BasisChange.to_sublattice_dual(Result->getDehomogenization());727esp_out << "dehomogenization" << endl;728}729esp_out.close();730}731if (tgn && Result->isComputed(ConeProperty::Generators))732Result->getGeneratorsMatrix().print(name,"tgn");733if (tri && Result->isComputed(ConeProperty::Triangulation)) { //write triangulation734write_tri();735}736737if (out==true) { //printing .out file738string name_open=name+".out"; //preparing output files739const char* file=name_open.c_str();740ofstream out(file);741if(out.fail()){742errorOutput() << "Cannot write to output file." << endl;743exit(1);744}745746// write "header" of the .out file747size_t nr_orig_gens = 0;748if (lattice_ideal_input) {749nr_orig_gens = Result->getNrOriginalMonoidGenerators();750out << nr_orig_gens <<" original generators of the toric ring"<<endl;751}752if (Result->isComputed(ConeProperty::ModuleGenerators)) {753out << Result->getNrModuleGenerators() <<" module generators" << endl;754}755if (Result->isComputed(ConeProperty::HilbertBasis)) {756out << Result->getNrHilbertBasis() <<" Hilbert basis elements"757<< of_monoid << endl;758}759if (homogeneous && Result->isComputed(ConeProperty::Deg1Elements)) {760out << Result->getNrDeg1Elements() <<" Hilbert basis elements of degree 1"<<endl;761}762if (Result->isComputed(ConeProperty::IsReesPrimary)763&& Result->isComputed(ConeProperty::HilbertBasis)) {764const Matrix<Integer>& Hilbert_Basis = Result->getHilbertBasisMatrix();765nr = Hilbert_Basis.nr_of_rows();766for (i = 0; i < nr; i++) {767if (Hilbert_Basis[i][dim-1]==1) {768rees_ideal_key.push_back(i);769}770}771out << rees_ideal_key.size() <<" generators of integral closure of the ideal"<<endl;772}773if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {774out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron" << endl;775}776if (Result->isComputed(ConeProperty::ExtremeRays)) {777out << Result->getNrExtremeRays() <<" extreme rays" << of_cone << endl;778}779if(Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {780out << Result->getNrModuleGeneratorsOverOriginalMonoid() <<" module generators over original monoid" << endl;781}782if (Result->isComputed(ConeProperty::SupportHyperplanes)) {783out << Result->getNrSupportHyperplanes() <<" support hyperplanes"784<< of_polyhedron << endl;785}786out<<endl;787if (Result->isComputed(ConeProperty::ExcludedFaces)) {788out << Result->getNrExcludedFaces() <<" excluded faces"<<endl;789out << endl;790}791out << "embedding dimension = " << dim << endl;792if (homogeneous) {793if (Result->isComputed(ConeProperty::Sublattice)) {794auto rank = Result->getRank();795out << "rank = "<< rank << is_maximal(rank,dim) << endl;796out << "external index = "<< Result->getSublattice().getExternalIndex() << endl;797}798} else { // now inhomogeneous case799if (Result->isComputed(ConeProperty::AffineDim))800out << "affine dimension of the polyhedron = "801<< Result->getAffineDim() << is_maximal(Result->getAffineDim(),dim-1) << endl;802if (Result->isComputed(ConeProperty::RecessionRank))803out << "rank of recession monoid = " << Result->getRecessionRank() << endl;804}805806if(Result->isComputed(ConeProperty::OriginalMonoidGenerators)){807out << "internal index = " << Result->getIndex() << endl;808}809810if(Result->isComputed(ConeProperty::MaximalSubspace)){811size_t dim_max_subspace=Result->getDimMaximalSubspace();812if(dim_max_subspace>0)813out << "dimension of maximal subspace = " << dim_max_subspace << endl;814}815816817if (homogeneous && Result->isComputed(ConeProperty::IsIntegrallyClosed)) {818if (Result->isIntegrallyClosed()) {819out << "original monoid is integrally closed"<<endl;820} else {821out << "original monoid is not integrally closed"<<endl;822if ( Result->isComputed(ConeProperty::IsIntegrallyClosed)823&& !Result->isComputed(ConeProperty::HilbertBasis)) {824out << "witness for not being integrally closed:" << endl;825out << Result->getWitnessNotIntegrallyClosed();826}827if(Result->getUnitGroupIndex()>1){828out << "unit group index = " << Result->getUnitGroupIndex() << endl;829}830}831}832out << endl;833if (Result->isComputed(ConeProperty::TriangulationSize)) {834out << "size of ";835if (Result->isTriangulationNested()) out << "nested ";836if (Result->isTriangulationPartial()) out << "partial ";837out << "triangulation = " << Result->getTriangulationSize() << endl;838}839if (Result->isComputed(ConeProperty::TriangulationDetSum)) {840out << "resulting sum of |det|s = " << Result->getTriangulationDetSum() << endl;841}842if (Result->isComputed(ConeProperty::TriangulationSize)) {843out << endl;844}845if ( Result->isComputed(ConeProperty::Dehomogenization) ) {846out << "dehomogenization:" << endl847<< Result->getDehomogenization() << endl;848}849if ( Result->isComputed(ConeProperty::Grading) ) {850out << "grading:" << endl851<< Result->getGrading();852Integer denom = Result->getGradingDenom();853if (denom != 1) {854out << "with denominator = " << denom << endl;855}856out << endl;857if (homogeneous && Result->isComputed(ConeProperty::ExtremeRays)) {858out << "degrees of extreme rays:"<<endl;859map<Integer,long> deg_count;860vector<Integer> degs = Result->getExtremeRaysMatrix().MxV(Result->getGrading());861for (i=0; i<degs.size(); ++i) {862deg_count[degs[i]/denom]++;863}864out << deg_count;865}866}867else if (Result->isComputed(ConeProperty::IsDeg1ExtremeRays)) {868if ( !Result->isDeg1ExtremeRays() ) {869out << "No implicit grading found" << endl;870}871}872out<<endl;873if (homogeneous && Result->isComputed(ConeProperty::IsDeg1HilbertBasis)874&& Result->isDeg1ExtremeRays() ) {875if (Result->isDeg1HilbertBasis()) {876out << "Hilbert basis elements are of degree 1";877} else {878out << "Hilbert basis elements are not of degree 1";879}880out<<endl<<endl;881}882if ( Result->isComputed(ConeProperty::ModuleRank) ) {883out << "module rank = "<< Result->getModuleRank() << endl;884}885if ( Result->isComputed(ConeProperty::Multiplicity) ) {886out << "multiplicity = "<< Result->getMultiplicity() << endl;887}888if ( Result->isComputed(ConeProperty::ModuleRank) || Result->isComputed(ConeProperty::Multiplicity)) {889out << endl;890}891892if ( Result->isComputed(ConeProperty::HilbertSeries) ) {893const HilbertSeries& HS = Result->getHilbertSeries();894vector<mpz_class> HS_Num;895map<long, long> HS_Denom;896if ( Result->isComputed(ConeProperty::HSOP) ){897HS_Denom=HS.getHSOPDenom();898HS_Num=HS.getHSOPNum();899out << "Hilbert series (HSOP):" << endl << HS_Num;900} else {901HS_Denom=HS.getDenom();902HS_Num=HS.getNum();903out << "Hilbert series:" << endl << HS_Num;904}905long nr_factors = 0;906for (map<long, long>::iterator it = HS_Denom.begin(); it!=HS_Denom.end(); ++it) {907nr_factors += it->second;908}909out << "denominator with " << nr_factors << " factors:" << endl;910out << HS_Denom;911out << endl;912if (HS.getShift() != 0) {913out << "shift = " << HS.getShift() << endl << endl;914}915916out << "degree of Hilbert Series as rational function = "917<< HS.getDegreeAsRationalFunction() << endl << endl;918if(v_is_symmetric(HS_Num)){919out << "The numerator of the Hilbert Series is symmetric." << endl << endl;920}921long period = HS.getPeriod();922if (period == 1 && (HS_Denom.size() == 0923|| HS_Denom.begin()->first== (long) HS_Denom.size())) {924out << "Hilbert polynomial:" << endl;925out << HS.getHilbertQuasiPolynomial()[0];926out << "with common denominator = ";927out << HS.getHilbertQuasiPolynomialDenom();928out << endl<< endl;929} else {930// output cyclonomic representation931out << "Hilbert series with cyclotomic denominator:" << endl;932out << HS.getCyclotomicNum();933out << "cyclotomic denominator:" << endl;934out << HS.getCyclotomicDenom();935out << endl;936// Hilbert quasi-polynomial937HS.computeHilbertQuasiPolynomial();938if (HS.isHilbertQuasiPolynomialComputed()) {939out<<"Hilbert quasi-polynomial of period " << period << ":" << endl;940if(HS.get_nr_coeff_quasipol()>=0){941out << "only " << HS.get_nr_coeff_quasipol() << " highest coefficients computed" << endl;942out << "their common period is " << HS.getHilbertQuasiPolynomial().size() << "" << endl;943}944Matrix<mpz_class> HQP(HS.getHilbertQuasiPolynomial());945HQP.pretty_print(out,true);946out<<"with common denominator = "<<HS.getHilbertQuasiPolynomialDenom();947}else{948out<<"Hilbert quasi-polynomial has period " << period << endl;949}950out << endl << endl;951}952953}954955if ( Result->isComputed(ConeProperty::WeightedEhrhartSeries) )956writeWeightedEhrhartSeries(out);957958959if ( Result->isComputed(ConeProperty::VirtualMultiplicity)960&& !Result->isComputed(ConeProperty::WeightedEhrhartQuasiPolynomial)) {961out << "virtual multiplicity of weighted Ehrhart series = "<< Result->getVirtualMultiplicity() << endl;962out << endl;963}964965if ( Result->isComputed(ConeProperty::Integral)) {966out << "integral = "<< Result->getIntegral() << endl;967out << endl;968}969970if (Result->isComputed(ConeProperty::IsReesPrimary)) {971if (Result->isReesPrimary()) {972out<<"ideal is primary to the ideal generated by the indeterminates"<<endl;973} else {974out<<"ideal is not primary to the ideal generated by the indeterminates"<<endl;975}976if (Result->isComputed(ConeProperty::ReesPrimaryMultiplicity)) {977out<<"multiplicity of the ideal = "<<Result->getReesPrimaryMultiplicity()<<endl;978}979out << endl;980}981982if(Result->isComputed(ConeProperty::ClassGroup)) {983vector<Integer> ClassGroup=Result->getClassGroup();984out << "rank of class group = " << ClassGroup[0] << endl;985if(ClassGroup.size()==1)986out << "class group is free" << endl << endl;987else{988ClassGroup.erase(ClassGroup.begin());989out << "finite cyclic summands:" << endl;990out << count_in_map<Integer,size_t>(ClassGroup);991out << endl;992}993}994995if(Result->isComputed(ConeProperty::IsGorenstein)) {996if(Result->isGorenstein()){997out << "Monoid is Gorenstein " << endl;998out << "Generator of interior:" << endl;999out << Result->getGeneratorOfInterior();1000}1001else1002out << "Monoid is not Gorenstein " << endl;1003out << endl;1004}10051006out << "***********************************************************************"1007<< endl << endl;100810091010if (lattice_ideal_input) {1011out << nr_orig_gens <<" original generators:"<<endl;1012Result->getOriginalMonoidGeneratorsMatrix().pretty_print(out);1013out << endl;1014}1015if (Result->isComputed(ConeProperty::ModuleGenerators)) {1016out << Result->getNrModuleGenerators() <<" module generators:" << endl;1017Result->getModuleGeneratorsMatrix().pretty_print(out);1018out << endl;1019}10201021if ( Result->isComputed(ConeProperty::Deg1Elements) ) {1022const Matrix<Integer>& Hom = Result->getDeg1ElementsMatrix();1023write_matrix_ht1(Hom);1024nr=Hom.nr_of_rows();1025out<<nr<<" Hilbert basis elements of degree 1:"<<endl;1026Hom.pretty_print(out);1027out << endl;1028}10291030if (Result->isComputed(ConeProperty::HilbertBasis)) {10311032const Matrix<Integer>& Hilbert_Basis = Result->getHilbertBasisMatrix();10331034if(!Result->isComputed(ConeProperty::Deg1Elements)){1035nr=Hilbert_Basis.nr_of_rows();1036out << nr << " Hilbert basis elements" << of_monoid << ":" << endl;1037Hilbert_Basis.pretty_print(out);1038out << endl;1039}1040else {1041nr=Hilbert_Basis.nr_of_rows()-Result->getNrDeg1Elements();1042out << nr << " further Hilbert basis elements" << of_monoid << " of higher degree:" << endl;1043Matrix<Integer> HighDeg(nr,dim);1044for(size_t i=0;i<nr;++i)1045HighDeg[i]=Hilbert_Basis[i+Result->getNrDeg1Elements()];1046HighDeg.pretty_print(out);1047out << endl;1048}1049Matrix<Integer> complete_Hilbert_Basis(0,dim);1050if (gen || egn || typ) {1051// for these files we append the module generators if there are any1052if (Result->isComputed(ConeProperty::ModuleGenerators)) {1053complete_Hilbert_Basis.append(Hilbert_Basis);1054complete_Hilbert_Basis.append(Result->getModuleGeneratorsMatrix());1055write_matrix_gen(complete_Hilbert_Basis);1056} else {1057write_matrix_gen(Hilbert_Basis);1058}1059}1060if ((egn || typ) && Result->isComputed(ConeProperty::Sublattice)) {1061const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();1062Matrix<Integer> Hilbert_Basis_Full_Cone = BasisChange.to_sublattice(Hilbert_Basis);1063if (Result->isComputed(ConeProperty::ModuleGenerators)) {1064Hilbert_Basis_Full_Cone.append(BasisChange.to_sublattice(Result->getModuleGeneratorsMatrix()));1065}1066if (egn) {1067string egn_string = name+".egn";1068const char* egn_file = egn_string.c_str();1069ofstream egn_out(egn_file);1070Hilbert_Basis_Full_Cone.print(egn_out);1071// egn_out<<"cone"<<endl;1072egn_out.close();1073}10741075if (typ && homogeneous) {1076write_matrix_typ(Hilbert_Basis_Full_Cone.multiplication(BasisChange.to_sublattice_dual(Result->getSupportHyperplanesMatrix()).transpose()));1077}1078}10791080if (Result->isComputed(ConeProperty::IsReesPrimary)) {1081out << rees_ideal_key.size() <<" generators of integral closure of the ideal:"<<endl;1082Matrix<Integer> Ideal_Gens = Hilbert_Basis.submatrix(rees_ideal_key);1083Ideal_Gens.resize_columns(dim-1);1084Ideal_Gens.pretty_print(out);1085out << endl;1086}1087}1088if (Result->isComputed(ConeProperty::VerticesOfPolyhedron) && !no_ext_rays_output) {1089out << Result->getNrVerticesOfPolyhedron() <<" vertices of polyhedron:" << endl;1090if(Result->isComputed(ConeProperty::VerticesFloat))1091write_float(out,Result->getVerticesFloatMatrix(),Result->getNrVerticesFloat(),dim);1092else1093Result->getVerticesOfPolyhedronMatrix().pretty_print(out);1094out << endl;1095}1096if (Result->isComputed(ConeProperty::ExtremeRays) && !no_ext_rays_output) {1097out << Result->getNrExtremeRays() << " extreme rays" << of_cone << ":" << endl;1098if(homogeneous && Result->isComputed(ConeProperty::VerticesFloat))1099write_float(out,Result->getVerticesFloatMatrix(),Result->getNrVerticesFloat(),dim);1100else1101Result->getExtremeRaysMatrix().pretty_print(out);1102out << endl;1103}11041105if (Result->isComputed(ConeProperty::ExtremeRays) && ext){1106// for the .gen file we append the vertices of polyhedron if there are any1107if (Result->isComputed(ConeProperty::VerticesOfPolyhedron)) {1108Matrix<Integer> Extreme_Rays(Result->getExtremeRaysMatrix());1109Extreme_Rays.append(Result->getVerticesOfPolyhedronMatrix());1110write_matrix_ext(Extreme_Rays);1111} else {1112write_matrix_ext(Result->getExtremeRaysMatrix());1113}1114}11151116if(Result->isComputed(ConeProperty::MaximalSubspace) && Result->getDimMaximalSubspace()>0){1117out << Result->getDimMaximalSubspace() <<" basis elements of maximal subspace:" << endl;1118Result->getMaximalSubspaceMatrix().pretty_print(out);1119out << endl;1120if(msp)1121write_matrix_msp(Result->getMaximalSubspaceMatrix());1122}11231124if(Result->isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {1125out << Result->getNrModuleGeneratorsOverOriginalMonoid() <<" module generators over original monoid:" << endl;1126Result->getModuleGeneratorsOverOriginalMonoidMatrix().pretty_print(out);1127out << endl;1128if(mod)1129write_matrix_mod(Result->getModuleGeneratorsOverOriginalMonoidMatrix());1130}11311132//write constrains (support hyperplanes, congruences, equations)11331134if (Result->isComputed(ConeProperty::SupportHyperplanes)) {1135const Matrix<Integer>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();1136out << Support_Hyperplanes.nr_of_rows() <<" support hyperplanes"1137<< of_polyhedron << ":" << endl;1138Support_Hyperplanes.pretty_print(out);1139out << endl;1140}1141if (Result->isComputed(ConeProperty::Sublattice)) {1142const Sublattice_Representation<Integer>& BasisChange = Result->getSublattice();1143//equations1144const Matrix<Integer>& Equations = BasisChange.getEquationsMatrix();1145size_t nr_of_equ = Equations.nr_of_rows();1146if (nr_of_equ > 0) {1147out << nr_of_equ <<" equations:" <<endl;1148Equations.pretty_print(out);1149out << endl;1150}11511152//congruences1153const Matrix<Integer>& Congruences = BasisChange.getCongruencesMatrix();1154size_t nr_of_cong = Congruences.nr_of_rows();1155if (nr_of_cong > 0) {1156out << nr_of_cong <<" congruences:" <<endl;1157Congruences.pretty_print(out);1158out << endl;1159}11601161//lattice1162const Matrix<Integer>& LatticeBasis = BasisChange.getEmbeddingMatrix();1163size_t nr_of_latt = LatticeBasis.nr_of_rows();1164if (nr_of_latt < dim || BasisChange.getExternalIndex()!=1) {1165out << nr_of_latt <<" basis elements of lattice:" <<endl;1166LatticeBasis.pretty_print(out);1167out << endl;1168}1169if(lat)1170write_matrix_lat(LatticeBasis);117111721173//excluded faces1174if (Result->isComputed(ConeProperty::ExcludedFaces)) {1175const Matrix<Integer>& ExFaces = Result->getExcludedFacesMatrix();1176out << ExFaces.nr_of_rows() <<" excluded faces:" <<endl;1177ExFaces.pretty_print(out);1178out << endl;1179}11801181if (cst && Result->isComputed(ConeProperty::SupportHyperplanes)) {1182const Matrix<Integer>& Support_Hyperplanes = Result->getSupportHyperplanesMatrix();1183string cst_string = name+".cst";1184const char* cst_file = cst_string.c_str();1185ofstream cst_out(cst_file);11861187Support_Hyperplanes.print(cst_out);1188cst_out<<"inequalities"<<endl;1189Equations.print(cst_out);1190cst_out<<"equations"<<endl;1191Congruences.print(cst_out);1192cst_out<<"congruences"<<endl;1193if (Result->isComputed(ConeProperty::ExcludedFaces)) {1194Result->getExcludedFacesMatrix().print(cst_out);1195cst_out<<"excluded_faces"<<endl;1196}1197if (Result->isComputed(ConeProperty::Grading)) {1198cst_out << 1 << endl << dim << endl;1199cst_out << Result->getGrading();1200cst_out << "grading" << endl;1201}1202if (Result->isComputed(ConeProperty::Dehomogenization)) {1203cst_out << 1 << endl << dim << endl;1204cst_out << Result->getDehomogenization();1205cst_out << "dehomogenization" << endl;1206}1207cst_out.close();1208}1209}12101211out.close();1212}12131214write_inv_file();1215write_Stanley_dec();1216}1217121812191220