CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

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

Views: 418451
Tweet
Share
Share
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
#ifndef CONE_H_

25
#define CONE_H_

26


27
#include <sys/stat.h>

28
#include <vector>

29
#include <map>

30
#include <utility> //for pair

31
#include <boost/dynamic_bitset.hpp>

32


33
#include <libnormaliz/libnormaliz.h>

34
#include <libnormaliz/cone_property.h>

35
#include <libnormaliz/sublattice_representation.h>

36
#include <libnormaliz/matrix.h>

37
#include <libnormaliz/HilbertSeries.h>

38


39
namespace libnormaliz {

40
using std::vector;

41
using std::map;

42
using std::pair;

43


44
template<typename Integer> class Full_Cone;

45
//template<typename Integer> class Matrix;

46


47
// type for simplex, short in contrast to class Simplex

48
template<typename Integer> struct SHORTSIMPLEX {

49
    vector<key_t> key;                // full key of simplex

50
    Integer height;                   // height of last vertex over opposite facet

51
    Integer vol;                      // volume if computed, 0 else

52
    vector<bool> Excluded;           // for disjoint decomposition of cone

53
                                      // true in position i indictate sthat the facet 

54
                                      // opposite of generator i must be excluded

55
};

56


57
template<typename Integer>

58
bool compareKeys(const SHORTSIMPLEX<Integer>& A, const SHORTSIMPLEX<Integer>& B){

59


60
    return(A.key < B.key);

61
}

62


63
struct STANLEYDATA_int { // for internal use

64
    vector<key_t> key;

65
    Matrix<long> offsets;

66
    vector<long> degrees; // degrees and classNr are used in nmz_integral.cpp

67
    size_t classNr;  // number of class of this simplicial cone

68
};

69


70
template<typename Integer> struct STANLEYDATA {

71
    vector<key_t> key;

72
    Matrix<Integer> offsets;

73
};

74


75
template<typename Integer>

76
class Cone {

77


78
//---------------------------------------------------------------------------

79
//                               public methods

80
//---------------------------------------------------------------------------

81
public:

82


83
//---------------------------------------------------------------------------

84
//                    Constructors, they preprocess the input

85
//---------------------------------------------------------------------------

86
    

87
    Cone(); //default constructor

88


89
    /* give up to 3 matrices as input

90
     * the types must be pairwise different

91
     */

92
    Cone(InputType type, const vector< vector<Integer> >& input_data);

93


94
    Cone(InputType type1, const vector< vector<Integer> >& input_data1,

95
         InputType type2, const vector< vector<Integer> >& input_data2);

96


97
    Cone(InputType type1, const vector< vector<Integer> >& input_data1,

98
         InputType type2, const vector< vector<Integer> >& input_data2,

99
         InputType type3, const vector< vector<Integer> >& input_data3);

100


101
    /* give multiple input */

102
    Cone(const map< InputType , vector< vector<Integer> > >& multi_input_data);

103
    

104
//-----------------------------------------------------------------------------

105
// the same for mpq_class

106
    

107
    Cone(InputType type, const vector< vector<mpq_class> >& input_data);

108


109
    Cone(InputType type1, const vector< vector<mpq_class> >& input_data1,

110
         InputType type2, const vector< vector<mpq_class> >& input_data2);

111


112
    Cone(InputType type1, const vector< vector<mpq_class> >& input_data1,

113
         InputType type2, const vector< vector<mpq_class> >& input_data2,

114
         InputType type3, const vector< vector<mpq_class> >& input_data3);

115


116
    /* give multiple input */

117
    Cone(const map< InputType , vector< vector<mpq_class> > >& multi_input_data);

118
    

119
//-----------------------------------------------------------------------------

120
// the same for nmz_float

121
    

122
    Cone(InputType type, const vector< vector<nmz_float> >& input_data);

123


124
    Cone(InputType type1, const vector< vector<nmz_float> >& input_data1,

125
         InputType type2, const vector< vector<nmz_float> >& input_data2);

126


127
    Cone(InputType type1, const vector< vector<nmz_float> >& input_data1,

128
         InputType type2, const vector< vector<nmz_float> >& input_data2,

129
         InputType type3, const vector< vector<nmz_float> >& input_data3);

130


131
    /* give multiple input */

132
    Cone(const map< InputType , vector< vector<nmz_float> > >& multi_input_data);

133


134
//-----------------------------------------------------------------------------

135
// Now with Matrix

136
    

137
    Cone(InputType type, const Matrix<Integer>& input_data);

138


139
    Cone(InputType type1, const Matrix<Integer>& input_data1,

140
         InputType type2, const Matrix<Integer>& input_data2);

141


142
    Cone(InputType type1, const Matrix<Integer>& input_data1,

143
         InputType type2, const Matrix<Integer>& input_data2,

144
         InputType type3, const Matrix<Integer>& input_data3);

145


146
    /* give multiple input */

147
    Cone(const map< InputType , Matrix<Integer> >& multi_input_data);

148
    

149
//-----------------------------------------------------------------------------

150
// Now with Matrix and mpq_class

151
    

152
    Cone(InputType type, const Matrix<mpq_class>& input_data);

153


154
    Cone(InputType type1, const Matrix<mpq_class>& input_data1,

155
         InputType type2, const Matrix<mpq_class>& input_data2);

156


157
    Cone(InputType type1, const Matrix<mpq_class>& input_data1,

158
         InputType type2, const Matrix<mpq_class>& input_data2,

159
         InputType type3, const Matrix<mpq_class>& input_data3);

160


161
    /* give multiple input */

162
    Cone(const map< InputType , Matrix<mpq_class> >& multi_input_data);

163
    

164
//-----------------------------------------------------------------------------

165
// Now with Matrix and nmz_float

166
    

167
    Cone(InputType type, const Matrix<nmz_float>& input_data);

168


169
    Cone(InputType type1, const Matrix<nmz_float>& input_data1,

170
         InputType type2, const Matrix<nmz_float>& input_data2);

171


172
    Cone(InputType type1, const Matrix<nmz_float>& input_data1,

173
         InputType type2, const Matrix<nmz_float>& input_data2,

174
         InputType type3, const Matrix<nmz_float>& input_data3);

175


176
    /* give multiple input */

177
    Cone(const map< InputType , Matrix<nmz_float> >& multi_input_data);

178
    

179
    

180
//---------------------------------------------------------------------------

181
//                                Destructor

182
//---------------------------------------------------------------------------

183


184
    ~Cone();

185


186
//---------------------------------------------------------------------------

187
//                          give additional data

188
//---------------------------------------------------------------------------

189


190
    /* Sets if the Cone prints verbose output.

191
     * The default value for the Cone is the global verbose.

192
     * returns the old value

193
     */

194
    bool setVerbose (bool v);

195


196
    void deactivateChangeOfPrecision();

197


198
//---------------------------------------------------------------------------

199
//                           make computations

200
//---------------------------------------------------------------------------

201


202
    // return what was NOT computed

203
    // ConeProperties compute(ComputationMode mode = Mode::hilbertBasisSeries); //default: everything

204
    ConeProperties compute_inner(ConeProperties ToCompute);

205
    // special case for up to 3 CPs

206
    ConeProperties compute(ConeProperties ToCompute);

207
    ConeProperties compute(ConeProperty::Enum);

208
    ConeProperties compute(ConeProperty::Enum, ConeProperty::Enum);

209
    ConeProperties compute(ConeProperty::Enum, ConeProperty::Enum, ConeProperty::Enum);

210


211
//---------------------------------------------------------------------------

212
//                         check what is computed

213
//---------------------------------------------------------------------------

214


215
    bool isComputed(ConeProperty::Enum prop) const;

216
    //returns true, when ALL properties in CheckComputed are computed

217
    bool isComputed(ConeProperties CheckComputed) const;

218
    

219
    void resetComputed(ConeProperty::Enum prop);

220


221
//---------------------------------------------------------------------------

222
//   get the results, these methods will start a computation if necessary

223
//   throws an NotComputableException if not succesful

224
//---------------------------------------------------------------------------

225


226
    // dimension and rank invariants

227
    size_t getEmbeddingDim() const { return dim; };   // is always known

228
    size_t getRank();                           // depends on ExtremeRays

229
    Integer getIndex(); // depends on OriginalMonoidGenerators

230
    Integer getInternalIndex(); // = getIndex()

231
    Integer getUnitGroupIndex(); // ditto

232
    // only for inhomogeneous case:

233
    size_t getRecessionRank();

234
    long getAffineDim();

235
    size_t getModuleRank();

236
    

237
    Cone<Integer>& getIntegerHullCone() const;

238
    Cone<Integer>& getSymmetrizedCone() const;

239


240
    const Matrix<Integer>& getGeneratorsMatrix();

241
    const vector< vector<Integer> >& getGenerators();

242
    size_t getNrGenerators();

243


244
    const Matrix<Integer>& getExtremeRaysMatrix();

245
    const vector< vector<Integer> >& getExtremeRays();

246
    size_t getNrExtremeRays();

247
    

248
    const Matrix<nmz_float>& getVerticesFloatMatrix();

249
    const vector< vector<nmz_float> >& getVerticesFloat();

250
    size_t getNrVerticesFloat();

251


252
    const Matrix<Integer>& getVerticesOfPolyhedronMatrix();

253
    const vector< vector<Integer> >& getVerticesOfPolyhedron();

254
    size_t getNrVerticesOfPolyhedron();

255


256
    const Matrix<Integer>& getSupportHyperplanesMatrix();

257
    const vector< vector<Integer> >& getSupportHyperplanes();

258
    size_t getNrSupportHyperplanes();

259
    

260
    const Matrix<Integer>& getMaximalSubspaceMatrix();

261
    const vector< vector<Integer> >& getMaximalSubspace();

262
    size_t getDimMaximalSubspace();

263


264
    // depends on the ConeProperty::s SupportHyperplanes and Sublattice

265
    map< InputType, vector< vector<Integer> > > getConstraints();

266


267
    const Matrix<Integer>& getExcludedFacesMatrix();

268
    const vector< vector<Integer> >& getExcludedFaces();

269
    size_t getNrExcludedFaces();

270


271
    size_t getTriangulationSize();

272
    Integer getTriangulationDetSum();

273


274
    vector<Integer> getWitnessNotIntegrallyClosed();

275
    vector<Integer> getGeneratorOfInterior();

276


277
    const Matrix<Integer>& getHilbertBasisMatrix();

278
    const vector< vector<Integer> >& getHilbertBasis();

279
    size_t getNrHilbertBasis();

280


281
    const Matrix<Integer>& getModuleGeneratorsOverOriginalMonoidMatrix();

282
    const vector< vector<Integer> >& getModuleGeneratorsOverOriginalMonoid();

283
    size_t getNrModuleGeneratorsOverOriginalMonoid();

284


285
    const Matrix<Integer>& getModuleGeneratorsMatrix();

286
    const vector< vector<Integer> >& getModuleGenerators();

287
    size_t getNrModuleGenerators();

288


289
    const Matrix<Integer>& getDeg1ElementsMatrix();

290
    const vector< vector<Integer> >& getDeg1Elements();

291
    size_t getNrDeg1Elements();

292


293
    // the actual grading is Grading/GradingDenom

294
    vector<Integer> getGrading();

295
    Integer getGradingDenom();

296


297
    vector<Integer> getDehomogenization();

298


299
    vector<Integer> getClassGroup();

300


301
    mpq_class getMultiplicity();

302
    mpq_class getVirtualMultiplicity();

303
    mpq_class getIntegral();

304
    const pair<HilbertSeries, mpz_class>& getWeightedEhrhartSeries();

305
    

306
    string getPolynomial() const;

307
    

308
    bool inequalities_present;

309


310
    bool isPointed();

311
    bool isInhomogeneous();

312
    bool isDeg1ExtremeRays();

313
    bool isDeg1HilbertBasis();

314
    bool isIntegrallyClosed();

315
    bool isGorenstein();

316
    bool isReesPrimary();

317
    Integer getReesPrimaryMultiplicity();

318
    const Matrix<Integer>& getOriginalMonoidGeneratorsMatrix();

319
    const vector< vector<Integer> >& getOriginalMonoidGenerators();

320
    size_t getNrOriginalMonoidGenerators();

321


322
    const Sublattice_Representation<Integer>& getSublattice();

323
    const HilbertSeries& getHilbertSeries(); //general purpose object

324


325
    // the following 2 methods give information about the last used triangulation

326
    // if no triangulation was computed so far they return false

327
    bool isTriangulationNested();

328
    bool isTriangulationPartial();

329
    const vector< pair<vector<key_t>, Integer> >& getTriangulation();

330
    const vector< vector<bool> >& getOpenFacets();

331
    const vector< pair<vector<key_t>, long> >& getInclusionExclusionData();

332
    const list< STANLEYDATA<Integer> >& getStanleyDec();

333
    list< STANLEYDATA_int >& getStanleyDec_mutable(); //allows us to erase the StanleyDec

334
                             // in order to save memeory for weighted Ehrhart

335
    

336
    void set_project(string name);

337
    void set_nmz_call(const string& path);

338
    void set_output_dir(string name);

339
    

340
    void setPolynomial(string poly);

341
    void setNrCoeffQuasiPol(long nr_coeff);

342
    

343
    bool get_verbose ();

344
    

345
    IntegrationData& getIntData();

346


347
//---------------------------------------------------------------------------

348
//                          private part

349
//---------------------------------------------------------------------------

350


351
private:

352
    

353
    bool already_in_compute; // protection against call of compute within compute

354
                             // such calls must go via recursive_compute.

355
    

356
    string project;

357
    string output_dir;

358
    string nmz_call;

359
    size_t dim;

360
    

361
    // the following three matrices store the constraints of the input

362
    Matrix<Integer> Inequalities;

363
    Matrix<Integer> Equations;

364
    Matrix<Integer> Congruences;

365
    // we must register some information about thew input

366
    bool lattice_ideal_input;

367
    size_t nr_latt_gen, nr_cone_gen; // they count matrices in the input 

368


369
    Sublattice_Representation<Integer> BasisChange;  //always use compose_basis_change() !

370
    Sublattice_Representation<Integer> BasisChangePointed; // to the pointed cone

371
    bool BC_set;

372
    bool verbose;

373
    ConeProperties is_Computed;

374
    // Matrix<Integer> GeneratorsOfToricRing;

375
    Matrix<Integer> OriginalMonoidGenerators;

376
    Matrix<Integer> Generators;

377
    Matrix<Integer> ExtremeRays;

378
    Matrix<nmz_float> VerticesFloat;

379
    vector<bool> ExtremeRaysIndicator;

380
    Matrix<Integer> VerticesOfPolyhedron;

381
    Matrix<Integer> SupportHyperplanes;

382
    Matrix<Integer> ExcludedFaces;

383
    Matrix<Integer> PreComputedSupportHyperplanes;

384
    size_t TriangulationSize;

385
    Integer TriangulationDetSum;

386
    bool triangulation_is_nested;

387
    bool triangulation_is_partial;

388
    vector< pair<vector<key_t>, Integer> > Triangulation;

389
    vector<vector<bool> > OpenFacets;

390
    vector< pair<vector<key_t>, long> > InExData;

391
    list< STANLEYDATA_int > StanleyDec;

392
    list< STANLEYDATA<Integer> > StanleyDec_export;

393
    mpq_class multiplicity;

394
    mpq_class Integral;

395
    mpq_class VirtualMultiplicity;

396
    vector<Integer> WitnessNotIntegrallyClosed;

397
    vector<Integer> GeneratorOfInterior;

398
    Matrix<Integer> HilbertBasis;

399
    Matrix<Integer> BasisMaxSubspace;

400
    Matrix<Integer> ModuleGeneratorsOverOriginalMonoid;

401
    Matrix<Integer> Deg1Elements;

402
    HilbertSeries HSeries;

403
    IntegrationData IntData;

404
    vector<Integer> Grading;

405
    vector<Integer> Dehomogenization;

406
    Integer GradingDenom;

407
    Integer index;  // the internal index

408
    Integer unit_group_index;

409
    

410
    vector<boost::dynamic_bitset<> > Pair; // for indicator vectors in project-and_lift

411
    vector<boost::dynamic_bitset<> > ParaInPair; // if polytope is a parallelotope

412
    bool check_parallelotope();

413


414
    bool pointed;

415
    bool inhomogeneous;

416
    bool gorensetin;

417
    bool deg1_extreme_rays;

418
    bool deg1_hilbert_basis;

419
    bool integrally_closed;

420
    bool Gorenstein;

421
    bool rees_primary;

422
    Integer ReesPrimaryMultiplicity;

423
    int affine_dim; //dimension of polyhedron

424
    size_t recession_rank; // rank of recession monoid

425
    size_t module_rank; // for the inhomogeneous case

426
    Matrix<Integer> ModuleGenerators;

427
    vector<Integer> ClassGroup;

428
    

429
    bool is_approximation;

430
    Cone* ApproximatedCone;

431
    

432
    // some properties of the current computation taken from ToCompute

433
    bool explicit_HilbertSeries; // true = Hilbert series set explicitly and not only via default mode

434
    bool naked_dual; // true = dual mode set, but neither Hilbert basis nor deg 1 points

435
    bool default_mode; // true default mode set

436


437
    Matrix<Integer> WeightsGrad;

438
    vector<bool> GradAbs;

439


440
    bool no_lattice_restriction; // true if cine generators are known to be in the relevant lattice

441
    bool normalization; // true if input type normalization is used

442


443
    // if this is true we allow to change to a smaller integer type in the computation

444
    bool change_integer_type;

445
    

446
    Cone<Integer>* IntHullCone;

447
    Cone<Integer>* SymmCone;

448
    

449
    // In cone based algorithms we use the following information

450
    bool Grading_Is_Coordinate; // indicates that the grading or dehomogenization is a coordinate

451
    key_t GradingCoordinate;  // namely this one

452


453
    void compose_basis_change(const Sublattice_Representation<Integer>& SR); // composes SR

454


455
    // main input processing

456
    void process_multi_input(const map< InputType, vector< vector<Integer> > >& multi_input_data);

457
    void process_multi_input_inner(map< InputType, vector< vector<Integer> > >& multi_input_data);

458
    void process_multi_input(const map< InputType, vector< vector<mpq_class> > >& multi_input_data);

459
    void process_multi_input(const map< InputType, vector< vector<nmz_float> > >& multi_input_data);

460
    

461
    void prepare_input_lattice_ideal(map< InputType, vector< vector<Integer> > >& multi_input_data);

462
    void prepare_input_constraints(const map< InputType, vector< vector<Integer> > >& multi_input_data);

463
    void prepare_input_generators(map< InputType, vector< vector<Integer> > >& multi_input_data,

464
                     Matrix<Integer>& LatticeGenerators);

465
    void homogenize_input(map< InputType, vector< vector<Integer> > >& multi_input_data);

466
    void check_precomputed_support_hyperplanes();

467
    void check_excluded_faces();

468
    

469
    void setGrading (const vector<Integer>& lf);

470
    void setWeights ();

471
    void setDehomogenization (const vector<Integer>& lf);

472
    void checkGrading();

473
    void checkDehomogenization();

474
    void check_vanishing_of_grading_and_dehom();

475
    void process_lattice_data(const Matrix<Integer>& LatticeGenerators, Matrix<Integer>& Congruences, Matrix<Integer>& Equations);

476
    

477
    ConeProperties recursive_compute(ConeProperties ToCompute);

478
    

479
    void try_symmetrization(ConeProperties& ToCompute);

480
    void try_approximation_or_projection (ConeProperties& ToCompute);

481


482
    Matrix<Integer> prepare_input_type_2(const vector< vector<Integer> >& Input);

483
    Matrix<Integer> prepare_input_type_3(const vector< vector<Integer> >& Input);

484
    void prepare_input_type_4(Matrix<Integer>& Inequalities);

485


486
    /* only used by the constructors */

487
    void initialize();

488


489
    template<typename IntegerFC>

490
    void compute_full_cone(ConeProperties& ToCompute);

491


492
    /* compute the generators using the support hyperplanes */

493
    void compute_generators();

494
    template<typename IntegerFC>

495
    void compute_generators_inner();

496


497
    /* compute method for the dual_mode, used in compute(mode) */

498
    void compute_dual(ConeProperties& ToCompute);

499
    template<typename IntegerFC>

500
    void compute_dual_inner(ConeProperties& ToCompute);

501
    

502
    void set_implicit_dual_mode(ConeProperties& ToCompute);

503


504
    /* extract the data from Full_Cone, this may remove data from Full_Cone!*/

505
    template<typename IntegerFC>

506
    void extract_data(Full_Cone<IntegerFC>& FC);

507
    template<typename IntegerFC>

508
    void extract_supphyps(Full_Cone<IntegerFC>& FC);

509
    

510
    void extract_supphyps(Full_Cone<Integer>& FC);

511


512


513
    /* set OriginalMonoidGenerators */

514
    void set_original_monoid_generators(const Matrix<Integer>&);

515


516
    /* set ExtremeRays, in inhomogeneous case also VerticesOfPolyhedron */

517
    void set_extreme_rays(const vector<bool>&);

518


519
    /* If the Hilbert basis and the original monoid generators are computed,

520
     * use them to check whether the original monoid is integrally closed. */

521
    void check_integrally_closed();

522
    void compute_unit_group_index();

523
    /* try to find a witness for not integrally closed in the Hilbert basis */

524
    void find_witness();

525
    

526
    void check_Gorenstein (ConeProperties&  ToCompute);

527


528
    Integer compute_primary_multiplicity();

529
    template<typename IntegerFC>

530
    Integer compute_primary_multiplicity_inner();

531
    

532
    void compute_integer_hull();

533
    void complete_sublattice_comp(ConeProperties& ToCompute); // completes the sublattice computations

534
    void complete_HilbertSeries_comp(ConeProperties& ToCompute);

535
    

536
    void compute_integral (ConeProperties& ToCompute);

537
    void compute_virt_mult (ConeProperties& ToCompute);

538
    void compute_weighted_Ehrhart(ConeProperties& ToCompute);

539
    

540
    void compute_vertices_float(ConeProperties& ToCompute);

541
    

542
    void make_StanleyDec_export();

543
    

544
    void NotComputable (string message); // throws NotComputableException if default_mode = false

545
    

546
    void set_parallelization();

547
    

548
    template<typename IntegerFC>

549
    void give_data_of_approximated_cone_to(Full_Cone<IntegerFC>& FC);

550
    

551
    void project_and_lift(Matrix<Integer>& Deg1, const Matrix<Integer>& Gens, Matrix<Integer>& Supps, bool float_projection);

552


553
};

554


555
// helpers

556


557
template<typename Integer>

558
vector<vector<Integer> > find_input_matrix(const map< InputType, vector< vector<Integer> > >& multi_input_data,

559
                               const InputType type);

560


561
template<typename Integer>

562
void insert_zero_column(vector< vector<Integer> >& mat, size_t col);

563


564
template<typename Integer>

565
void insert_column(vector< vector<Integer> >& mat, size_t col, Integer entry);

566


567


568
}  //end namespace libnormaliz

569


570
#endif /* CONE_H_ */

571


572