1 /*
2 * Normaliz
3 * Copyright (C) 2007-2019 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 LIBNORMALIZ_FULL_CONE_H
25 #define LIBNORMALIZ_FULL_CONE_H
26
27 #include <list>
28 #include <vector>
29 #include <deque>
30 #include <chrono>
31 //#include <set>
32
33 #include "libnormaliz/general.h"
34 #include "libnormaliz/cone.h"
35 //#include "libnormaliz/cone_property.h"
36 #include "libnormaliz/matrix.h"
37 #include "libnormaliz/simplex.h"
38 #include "libnormaliz/cone_dual_mode.h"
39 #include "libnormaliz/HilbertSeries.h"
40 #include "libnormaliz/reduction.h"
41 // #include "libnormaliz/sublattice_representation.h"
42 #include "libnormaliz/offload_handler.h"
43 #include "libnormaliz/automorph.h"
44 #include "libnormaliz/dynamic_bitset.h"
45
46 namespace libnormaliz {
47 using std::list;
48 using std::map;
49 using std::pair;
50 using std::vector;
51
52 template <typename Integer>
53 class Cone;
54 template <typename Integer>
55 class SimplexEvaluator;
56 template <typename Integer>
57 class CandidateList;
58 template <typename Integer>
59 class Candidate;
60 template <typename Integer>
61 class Simplex;
62 template <typename Integer>
63 class Collector;
64 template <typename Integer>
65 class Cone_Dual_Mode;
66 template <typename Integer>
67 struct FACETDATA;
68
69 template <typename Integer>
70 class Full_Cone {
71 friend class Cone<Integer>;
72 friend class SimplexEvaluator<Integer>;
73 friend class CandidateList<Integer>;
74 friend class Candidate<Integer>;
75 friend class Collector<Integer>;
76
77 public:
78 int omp_start_level; // records the omp_get_level() when the computation is started
79 // recorded at the start of the top cone constructor and the compute functions
80 // compute and dualize_cone
81 size_t dim;
82 size_t level0_dim; // dim of cone in level 0 of the inhomogeneous case
83 size_t module_rank; // rank of solution module over level 0 monoid in the inhomogeneous case
84 size_t nr_gen;
85 // size_t hyp_size; // not used at present
86 Integer index; // index of full lattice over lattice of generators
87
88 bool verbose;
89
90 bool keep_convex_hull_data;
91
92 bool pointed;
93 bool is_simplicial;
94 bool deg1_generated_computed;
95 bool deg1_generated;
96 bool deg1_extreme_rays;
97 bool deg1_triangulation;
98 bool deg1_hilbert_basis;
99 bool inhomogeneous;
100
101 // control of what to compute (set from outside)
102 bool explicit_full_triang; // indicates whether full triangulation is asked for without default mode
103 // bool explicit_h_vector; // to distinguish it from being set via default mode --DONE VIA do_default_mode
104 bool do_determinants;
105 bool do_multiplicity;
106 bool do_integral;
107 bool do_integrally_closed;
108 bool do_Hilbert_basis;
109 bool do_deg1_elements;
110 bool do_h_vector;
111 bool keep_triangulation;
112 bool pulling_triangulation;
113 bool keep_triangulation_bitsets; // convert the triangulation keys into bitsets and keep them
114 bool do_Stanley_dec;
115 bool do_default_mode;
116 bool do_class_group;
117 bool do_module_gens_intcl;
118 bool do_module_rank;
119 bool do_cone_dec;
120 bool do_supphyps_dynamic; // for integer hull computations where we want to insert extreme rays only
121 // more or less ...
122 bool do_multiplicity_by_signed_dec;
123 bool do_integral_by_signed_dec;
124 bool do_signed_dec;
125 bool do_virtual_multiplicity_by_signed_dec;
126 bool include_dualization; // can only be set in connection with signed dec
127 bool do_pure_triang; // no determinants
128
129 bool exploit_automs_mult;
130 bool exploit_automs_vectors;
131 bool do_automorphisms;
132 bool check_semiopen_empty;
133
134 bool do_hsop;
135 bool do_extreme_rays;
136 bool do_pointed;
137 bool believe_pointed; // sometimes set to suppress the check for pointedness
138 bool do_triangulation_size;
139
140 // algorithmic variants
141 bool do_approximation;
142 bool do_bottom_dec;
143 bool suppress_bottom_dec;
144 bool keep_order;
145
146 bool hilbert_basis_rec_cone_known;
147
148 // control of triangulation and evaluation
149 bool do_triangulation;
150 bool do_partial_triangulation;
151 bool do_only_multiplicity;
152 bool stop_after_cone_dec;
153 bool do_evaluation;
154 bool triangulation_is_nested;
155 bool triangulation_is_partial;
156
157 // type of definition of automorphism group
158 AutomParam::Quality quality_of_automorphisms;
159
160 // internal helper control variables
161 bool use_existing_facets; // in order to avoid duplicate computation of already computed facets
162 bool do_excluded_faces;
163 bool no_descent_to_facets; // primal algorithm must be applied to God_Father
164 bool do_only_supp_hyps_and_aux;
165 bool do_all_hyperplanes; // controls whether all support hyperplanes must be computed
166 bool use_bottom_points;
167 ConeProperties is_Computed;
168 bool has_generator_with_common_divisor;
169
170 long autom_codim_vectors; // bound for the descent to faces in algorithms using automorphisms
171 long autom_codim_mult; // bound ditto for multiplicity
172 Integer HB_bound; // only degree bound used in connection with automorphisms
173 // to discard vectors quickly
174 long block_size_hollow_tri;
175 long decimal_digits;
176 string project_name;
177
178 bool time_measured;
179 bool don_t_add_hyperplanes; // blocks the addition of new hyperplanes during time measurement
180 bool take_time_of_large_pyr; // if true, the time of large pyrs is measured
181 vector<chrono::nanoseconds> time_of_large_pyr;
182 vector<chrono::nanoseconds> time_of_small_pyr;
183 vector<size_t> nr_pyrs_timed;
184
185 // data of the cone (input or output)
186 vector<Integer> Truncation; // used in the inhomogeneous case to suppress vectors of level > 1
187 vector<Integer> Norm; // is Truncation or Grading, used to "simplify" renf_elem_vectors
188 vector<Integer> IntHullNorm; // used in computation of integer hulls for guessing extreme rays
189 Integer TruncLevel; // used for approximation of simplicial cones
190 vector<Integer> Grading;
191 vector<Integer> GradingOnPrimal; // grading on ther cone whose multiplicity is comuted by signed dec
192 vector<Integer> Sorting;
193 mpq_class multiplicity;
194 #ifdef ENFNORMALIZ
195 renf_elem_class renf_multiplicity;
196 #endif
197 Matrix<Integer> Generators;
198 Matrix<Integer> InputGenerators; // stores purified input -- Generators can be extended
199 set<vector<Integer>> Generator_Set; // the generators as a set (if needed)
200 Matrix<nmz_float> Generators_float; // floatung point approximations to the generators
201 vector<key_t> PermGens; // stores the permutation of the generators created by sorting
202 vector<bool> Extreme_Rays_Ind;
203 Matrix<Integer> Support_Hyperplanes;
204 Matrix<Integer> HilbertBasisRecCone;
205 Matrix<Integer> Subcone_Support_Hyperplanes; // used if *this computes elements in a subcone, for example in approximation
206 Matrix<Integer> Subcone_Equations;
207 vector<Integer> Subcone_Grading;
208 size_t nrSupport_Hyperplanes;
209 list<vector<Integer>> Hilbert_Basis;
210 vector<Integer> Witness; // for not integrally closed
211 Matrix<Integer>
212 Basis_Max_Subspace; // a basis of the maximal linear subspace of the cone --- only used in connection with dual mode
213 list<vector<Integer>> ModuleGeneratorsOverOriginalMonoid;
214 CandidateList<Integer> OldCandidates, NewCandidates, HBRC, ModuleGensDepot; // for the Hilbert basis
215 // HBRC is for the Hilbert basis of the recession cone if provided, ModuleGensDepot for the collected module
216 // generators in this case
217 size_t CandidatesSize;
218 list<vector<Integer>> Deg1_Elements;
219 HilbertSeries Hilbert_Series;
220 vector<Integer> gen_degrees; // will contain the degrees of the generators
221 vector<long> gen_degrees_long; // will contain the degrees of the generators as long (for h-vector)
222 Integer shift; // needed in the inhomogeneous case to make degrees positive
223 vector<Integer> gen_levels; // will contain the levels of the generators (in the inhomogeneous case)
224 size_t TriangulationBufferSize; // number of elements in Triangulation, for efficiency
225 list<SHORTSIMPLEX<Integer>> Triangulation; // triangulation of cone
226 vector<pair<dynamic_bitset,dynamic_bitset> > Triangulation_ind; // the same, but bitsets instead of keys
227 list<SHORTSIMPLEX<Integer>> TriangulationBuffer; // simplices to evaluate
228 list<SimplexEvaluator<Integer>> LargeSimplices; // Simplices for internal parallelization
229 Integer detSum; // sum of the determinants of the simplices
230 list<STANLEYDATA_int> StanleyDec; // Stanley decomposition
231 vector<Integer>
232 ClassGroup; // the class group as a vector: ClassGroup[0]=its rank, then the orders of the finite cyclic summands
233
234 Matrix<Integer> ProjToLevel0Quot; // projection matrix onto quotient modulo level 0 sublattice
235
236 size_t index_covering_face; //used in checking emptyness of semiopen polyhedron
237
238 string Polynomial;
239 mpq_class Integral, VirtualMultiplicity;
240 nmz_float RawEuclideanIntegral;
241 long DegreeOfPolynomial;
242
243 // ************************** Data for convex hull computations ****************************
244 vector<size_t> HypCounter; // counters used to give unique number to hyperplane
245 // must be defined thread wise to avoid critical
246
247 vector<bool> in_triang; // intriang[i]==true means that Generators[i] has been actively inserted
248 vector<key_t> GensInCone; // lists the generators completely built in
249 size_t nrGensInCone; // their number
250
251 vector<size_t> Comparisons; // at index i we note the total number of comparisons
252 // of positive and negative hyperplanes needed for the first i generators
253 size_t nrTotalComparisons; // counts the comparisons in the current computation
254
255 /* struct FACETDATA<Integer>{
256 vector<Integer> Hyp; // linear form of the hyperplane
257 dynamic_bitset GenInHyp; // incidence hyperplane/generators
258 Integer ValNewGen; // value of linear form on the generator to be added
259 size_t BornAt; // number of generator (in order of insertion) at which this hyperplane was added,,
260 counting from 0 size_t Ident; // unique number identifying the hyperplane (derived from HypCounter)
261 size_t Mother; // Ident of positive mother if known, 0 if unknown
262 bool is_positive_on_all_original_gens;
263 bool is_negative_on_some_original_gen;
264 bool simplicial; // indicates whether facet is simplicial
265 };*/
266
267 list<FACETDATA<Integer>> Facets; // contains the data for Fourier-Motzkin and extension of triangulation
268 size_t old_nr_supp_hyps; // must be remembered since Facets gets extended before the current generators is finished
269
270 // ******************************************************************************************
271
272 // Pointer to the cone by which the Full_Cone has been constructed (if any)
273 // Cone<Integer>* Creator;
274 Matrix<Integer> Embedding; // temporary solution -- at present used for integration with signed dec
275
276 // the absolute top cone in recursive algorithms where faces are evalutated themselves
277 // Full_Cone<Integer>* God_Father; // not used at present
278
279 // data relating a pyramid to its ancestores
280 Full_Cone<Integer>* Top_Cone; // reference to cone on top level relative to pyramid formation
281
282 vector<key_t> Top_Key; // indices of generators w.r.t Top_Cone
283 Full_Cone<Integer>* Mother; // reference to the mother of the pyramid
284 vector<key_t> Mother_Key; // indices of generators w.r.t Mother
285 size_t apex; // indicates which generator of mother cone is apex of pyramid
286 int pyr_level; // -1 for top cone, increased by 1 for each level of pyramids
287
288 int descent_level; // measures the decent in recursive algorithms that exploit compute_automorphisms
289 // 0 for God_father, increases by 1 with each passge to a facet
290
291 Isomorphism_Classes<Integer> FaceClasses;
292
293 vector<bool> IsLarge; // additional information whether pyramid is large
294
295 // control of pyramids, recusrion and parallelization
296 bool is_pyramid; // false for top cone
297 long top_last_to_be_inserted; // used for signed dec to avoid storage of hyperplanes that are not needed
298 bool pyramids_for_last_built_directly; // ditto
299 bool recursion_allowed; // to allow or block recursive formation of pyamids
300 bool multithreaded_pyramid; // indicates that this cone is computed in parallel threads
301 bool tri_recursion; // true if we have gone to pyramids because of triangulation
302
303 // storage for subpyramids
304 size_t store_level; // the level on which daughters will be stored
305 deque<list<vector<key_t>>> Pyramids; // storage for pyramids
306 deque<size_t> nrPyramids; // number of pyramids on the various levels
307 deque<bool> Pyramids_scrambled; // only used for mic
308
309 // data that can be used to go out of build_cone and return later (not done at present)
310 // but also useful at other places
311 // long nextGen; // the next generator to be processed
312 long lastGen; // the last generator processed
313
314 // Helpers for triangulation and Fourier-Motzkin
315 vector<typename list<SHORTSIMPLEX<Integer>>::iterator> TriSectionFirst; // first simplex with lead vertex i
316 vector<typename list<SHORTSIMPLEX<Integer>>::iterator> TriSectionLast; // last simplex with lead vertex i
317 list<FACETDATA<Integer>> LargeRecPyrs; // storage for large recusive pyramids given by basis of pyramid in mother cone
318
319 list<SHORTSIMPLEX<Integer>> FreeSimpl; // list of short simplices already evaluated, kept for recycling
320 vector<list<SHORTSIMPLEX<Integer>>> FS; // the same per thread
321 vector<Matrix<Integer> > RankTest; // helper matrices for rank test
322 vector<Matrix<Integer>> WorkMat; // helper matrix for matrix inversion
323 Matrix<Integer> UnitMat; // prefabricated unit matrix
324 vector<Matrix<nmz_float> > RankTest_float; // helper matrices for rank test
325
326 // helpers for evaluation
327 vector<SimplexEvaluator<Integer>> SimplexEval; // one per thread
328 vector<Collector<Integer>> Results; // one per thread
329 vector<Integer> Order_Vector; // vector for the disjoint decomposition of the cone
330 #ifdef NMZ_MIC_OFFLOAD
331 MicOffloader<long long> mic_offloader;
332 #endif
333 void try_offload_loc(long place, size_t max_level);
334
335 template <typename IntegerCone>
336 void restore_previous_computation(CONVEXHULLDATA<IntegerCone>& ConvHullData, bool goal);
337
338 template <typename IntegerCone>
339 void dualize_and_restore(CONVEXHULLDATA<IntegerCone>& ConvHullData);
340
341 // defining semiopen cones
342 Matrix<Integer> ExcludedFaces;
343 map<dynamic_bitset, long> InExCollect;
344
345 // statistics
346 size_t totalNrSimplices; // total number of simplices evaluated
347 size_t nrSimplicialPyr;
348 size_t totalNrPyr;
349
350 size_t start_from;
351
352 size_t AdjustedReductionBound;
353
354 AutomorphismGroup<Integer> Automs;
355
356 bool is_global_approximation; // true if approximation is defined in Cone
357
358 vector<vector<key_t>> approx_points_keys;
359 Matrix<Integer> OriginalGenerators;
360
361 Integer VolumeBound; // used to stop computation of approximation if simplex of this has larger volume
362
363 long renf_degree;
364
365 /* ---------------------------------------------------------------------------
366 * Private routines, used in the public routines
367 * ---------------------------------------------------------------------------
368 */
369 void number_hyperplane(FACETDATA<Integer>& hyp, const size_t born_at, const size_t mother);
370 bool is_hyperplane_included(FACETDATA<Integer>& hyp);
371 vector<Integer> FM_comb(const vector<Integer>& Pos,
372 const Integer& PosVal,
373 const vector<Integer>& Neg,
374 const Integer& NegVal,
375 bool extract_gcd = true);
376 void add_hyperplane(const size_t& new_generator,
377 const FACETDATA<Integer>& positive,
378 const FACETDATA<Integer>& negative,
379 list<FACETDATA<Integer>>& NewHyps,
380 bool known_to_be_simplicial);
381 void make_pyramid_for_last_generator(const FACETDATA<Integer>& Fac);
382 void extend_triangulation(const size_t& new_generator);
383 void update_pulling_triangulation(const size_t& new_generator);
384 void find_new_facets(const size_t& new_generator);
385 void process_pyramids(const size_t new_generator, const bool recursive);
386 void process_pyramid(const vector<key_t>& Pyramid_key,
387 const size_t new_generator,
388 const size_t store_level,
389 Integer height,
390 const bool recursive,
391 typename list<FACETDATA<Integer>>::iterator hyp,
392 size_t start_level);
393 void select_supphyps_from(list<FACETDATA<Integer>>& NewFacets,
394 const size_t new_generator,
395 const vector<key_t>& Pyramid_key,
396 const vector<bool>& Pyr_in_triang);
397 bool check_pyr_buffer(const size_t level);
398 void evaluate_stored_pyramids(const size_t level);
399 void match_neg_hyp_with_pos_hyps(const FACETDATA<Integer>& Neg,
400 size_t new_generator,
401 const vector<FACETDATA<Integer>*>& PosHyps,
402 dynamic_bitset& Zero_P,
403 vector<list<dynamic_bitset>>& Facets_0_1);
404 void collect_pos_supphyps(vector<FACETDATA<Integer>*>& PosHyps, dynamic_bitset& Zero_P, size_t& nr_pos);
405 void evaluate_rec_pyramids(const size_t level);
406 void evaluate_large_rec_pyramids(size_t new_generator);
407
408 void find_and_evaluate_start_simplex();
409 // Simplex<Integer> find_start_simplex() const;
410 vector<key_t> find_start_simplex() const;
411 void store_key(const vector<key_t>&,
412 const Integer& height,
413 const Integer& mother_vol,
414 list<SHORTSIMPLEX<Integer>>& Triangulation);
415 void find_bottom_facets();
416
417 void convert_polyhedron_to_polytope();
418
419 size_t make_hollow_triangulation_inner(const vector<size_t>& Selection,
420 const vector<key_t>& PatternKey, const dynamic_bitset& Pattern);
421 size_t refine_and_process_selection (vector<size_t>& Selection,
422 const vector<key_t>& PatternKey, const dynamic_bitset& Pattern, size_t& nr_subfacets);
423 size_t extend_selection_pattern(vector<size_t>& Selection,
424 const vector<key_t>& PatternKey, const dynamic_bitset& Pattern, size_t& nr_subfacets);
425 size_t make_hollow_triangulation();
426 void compute_multiplicity_or_integral_by_signed_dec();
427
428 void first_subfacet(const Matrix<Integer>& Generators, const dynamic_bitset& Subfacet,
429 const vector<Integer>& GradingOnPrimal, Matrix<Integer>& PrimalSimplex, const vector<Integer>& Generic,
430 Integer& MultPrimal, vector<Integer>& DegreesPrimal, const Matrix<Integer>& CandidatesGeneric,
431 Matrix<Integer>& ValuesGeneric);
432
433 void next_subfacet(const dynamic_bitset& Subfacet_next, const dynamic_bitset& Subfacet_start,
434 const Matrix<Integer>& Generators, const Matrix<Integer>& PrimalSimplex,
435 const Integer& MultPrimal, const vector<Integer>& DegreesPrimal, vector<Integer>& NewDegrees,
436 Integer& NewMult, const Matrix<Integer>& CandidatesGeneric,
437 const Matrix<Integer>& ValuesGeneric, Matrix<Integer>& NewValues);
438 bool process_hollow_triang(const vector<list<dynamic_bitset> >& SubFacetsBySimplex,
439 const vector<Integer>& Generic, const Matrix<Integer>& Generators,
440 const vector<Integer>& GradingOnPrimal, Matrix<Integer>& CandidatesGeneric,
441 vector<Integer>& OurCandidate);
442
443 void make_facet_triang(list<vector<key_t> >& FacetTriang, const FACETDATA<Integer>& Facet);
444
445 void compute_deg1_elements_via_projection_simplicial(const vector<key_t>& key); // for a simplicial subcone by projecion
446 void compute_sub_div_elements(const Matrix<Integer>& gens,
447 list<vector<Integer>>& sub_div_elements,
448 bool best_point = false); // computes subdividing elements via approximation
449 void select_deg1_elements(const Full_Cone& C);
450 // void select_Hilbert_Basis(const Full_Cone& C); //experimental, unused
451
452 void build_top_cone();
453 void build_cone_dynamic();
454 void build_cone();
455 void get_supphyps_from_copy(
456 bool from_scratch,
457 bool with_extreme_rays = false); // if evealuation starts before support hyperplanes are fully computed
458 void update_reducers(bool forced = false); // update list of reducers after evaluation of simplices
459
460 bool is_reducible(list<vector<Integer>*>& Irred, const vector<Integer>& new_element);
461 void global_reduction();
462
463 vector<Integer> compute_degree_function() const;
464
465 Matrix<Integer> select_matrix_from_list(const list<vector<Integer>>& S, vector<size_t>& selection);
466
467 bool contains(const vector<Integer>& v);
468 bool subcone_contains(const vector<Integer>& v);
469 bool contains(const Full_Cone& C);
470 void extreme_rays_and_deg1_check();
471 void find_grading();
472 void find_grading_inhom();
473 void check_given_grading();
474 void disable_grading_dep_comp();
475 void set_degrees();
476 void set_levels(); // for truncation in the inhomogeneous case
477 void find_module_rank(); // finds the module rank in the inhom case
478 void find_module_rank_from_HB();
479 void find_module_rank_from_proj(); // used if Hilbert basis is not computed
480 void find_level0_dim(); // ditto for the level 0 dimension
481 void find_level0_dim_from_HB(); // from the Hilbert basis (after dual mode)
482 void sort_gens_by_degree(bool triangulate);
483 // void compute_support_hyperplanes(bool do_extreme_rays=false);
484 bool check_evaluation_buffer();
485 bool check_evaluation_buffer_size();
486 void prepare_old_candidates_and_support_hyperplanes();
487 void evaluate_triangulation();
488 void evaluate_large_simplices();
489 void evaluate_large_simplex(size_t j, size_t lss);
490 void transfer_triangulation_to_top();
491 void primal_algorithm();
492 void primal_algorithm_initialize();
493 void primal_algorithm_finalize();
494 void primal_algorithm_set_computed();
495 void finish_Hilbert_series();
496 void make_module_gens();
497 void make_module_gens_and_extract_HB();
498 void remove_duplicate_ori_gens_from_HB();
499 void compute_class_group();
500 void compose_perm_gens(const vector<key_t>& perm);
501 void check_grading_after_dual_mode();
502
503 void multiplicity_by_signed_dec();
504
505 void minimize_support_hyperplanes();
506 void compute_extreme_rays(bool use_facets = false);
507 void compute_extreme_rays_compare(bool use_facets);
508 void compute_extreme_rays_rank(bool use_facets);
509 void select_deg1_elements();
510
511 void check_pointed();
512 void deg1_check();
513 void check_deg1_extreme_rays();
514 void check_deg1_hilbert_basis();
515
516 void compute_multiplicity();
517
518 void minimize_excluded_faces();
519 void prepare_inclusion_exclusion();
520
521 void set_preconditions();
522 void set_primal_algorithm_control_variables();
523 void reset_tasks();
524 void deactivate_completed_tasks();
525 void addMult(Integer& volume, const vector<key_t>& key, const int& tn); // multiplicity sum over thread tn
526
527 void check_simpliciality_hyperplane(const FACETDATA<Integer>& hyp) const;
528 void check_facet(const FACETDATA<Integer>& Fac, const size_t& new_generator) const; // debugging routine
529 void set_simplicial(FACETDATA<Integer>& hyp);
530
531 void compute_hsop();
532 void heights(list<vector<key_t>>& facet_keys,
533 list<pair<dynamic_bitset, size_t>> faces,
534 size_t index,
535 vector<size_t>& ideal_heights,
536 size_t max_dim);
537
538 void start_message();
539 void end_message();
540
541 void set_zero_cone();
542
543 void compute_automorphisms(size_t nr_special_gens = 0);
544 void compute_by_automorphisms();
545 mpq_class facet_multiplicity(const vector<key_t>& facet_key);
546 void compute_multiplicity_via_automs();
547 vector<vector<key_t>> get_facet_keys_for_orbits(const vector<Integer>& fixed_point, bool with_orbit_sizes);
548 vector<Integer> get_fixed_point(size_t nr_cone_points);
549 void compute_HB_via_automs();
550 vector<Integer> replace_fixed_point_by_generator(const vector<Integer>& fixed_point,
551 const key_t facet_nr,
552 const vector<Integer>& help_grading);
553 void compute_Deg1_via_automs();
554 void get_cone_over_facet_vectors(const vector<Integer>& fixed_point,
555 const vector<key_t>& facet_key,
556 const key_t facet_nr,
557 list<vector<Integer>>& facet_vectors);
558 Matrix<Integer> push_supphyps_to_cone_over_facet(const vector<Integer>& fixed_point, const key_t facet_nr);
559 void import_HB_from(const IsoType<Integer>& copy);
560 // bool check_extension_to_god_father();
561 void compute_multiplicity_via_recession_cone();
562 void copy_autom_params(const Full_Cone<Integer>& C);
563
564 void recursive_revlex_triangulation(vector<key_t> simplex_so_far,
565 const vector<key_t>& gens_in_face,
566 const vector<typename list<FACETDATA<Integer>>::const_iterator>& mother_facets,
567 size_t dim);
568 void make_facets();
569 void revlex_triangulation();
570
571 chrono::nanoseconds rank_time();
572 chrono::nanoseconds cmp_time();
573 chrono::nanoseconds ticks_comp_per_supphyp;
574 chrono::nanoseconds ticks_rank_per_row;
575 chrono::nanoseconds ticks_per_cand;
576
577 void small_vs_large(const size_t new_generator); // compares computation times of small vs. large pyramids
578
579 #ifdef NMZ_MIC_OFFLOAD
580 void try_offload(size_t max_level);
581 #else
try_offload(size_t max_level)582 void try_offload(size_t max_level){};
583 #endif
584
585 /*---------------------------------------------------------------------------
586 * Constructors
587 *---------------------------------------------------------------------------
588 */
589 Full_Cone(); // default constructor
590 Full_Cone(const Matrix<Integer>& M, bool do_make_prime = true); // main constructor
591 Full_Cone(Cone_Dual_Mode<Integer>& C); // removes data from the argument!
592 Full_Cone(Full_Cone<Integer>& C, const vector<key_t>& Key); // for pyramids
593
594 /*---------------------------------------------------------------------------
595 * Data access
596 *---------------------------------------------------------------------------
597 */
598 void print() const; // to be modified, just for tests
599 size_t getDimension() const;
600 size_t getNrGenerators() const;
601 bool isPointed() const;
602 bool isDeg1ExtremeRays() const;
603 bool isDeg1HilbertBasis() const;
604 vector<Integer> getGrading() const;
605 mpq_class getMultiplicity() const;
606 Integer getShift() const;
607 size_t getModuleRank() const;
608 const Matrix<Integer>& getGenerators() const;
609 vector<bool> getExtremeRays() const;
610 size_t getNrExtremeRays() const;
611 Matrix<Integer> getSupportHyperplanes() const;
612 Matrix<Integer> getHilbertBasis() const;
613 Matrix<Integer> getModuleGeneratorsOverOriginalMonoid() const;
614 Matrix<Integer> getDeg1Elements() const;
615 vector<Integer> getHVector() const;
616 Matrix<Integer> getExcludedFaces() const;
617
618 bool isComputed(ConeProperty::Enum prop) const;
619 void setComputed(ConeProperty::Enum prop);
620 void setComputed(ConeProperty::Enum prop, bool value);
621
622 /*---------------------------------------------------------------------------
623 * Computation Methods
624 *---------------------------------------------------------------------------
625 */
626 void dualize_cone(bool print_message = true);
627 void support_hyperplanes();
628
629 void compute();
630
631 /* adds generators, they have to lie inside the existing cone */
632 void add_generators(const Matrix<Integer>& new_points);
633
634 void dual_mode();
635
636 void error_msg(string s) const;
637 };
638 // class end *****************************************************************
639
640 template <typename Integer>
641 template <typename IntegerCone>
dualize_and_restore(CONVEXHULLDATA<IntegerCone> & ConvHullData)642 void Full_Cone<Integer>::dualize_and_restore(CONVEXHULLDATA<IntegerCone>& ConvHullData) {
643 // goal=true: to primal, goal=false: to dual
644
645 /* ConvHullData.Generators.pretty_print(cout);
646 cout << "===============" << endl;
647 Generators.pretty_print(cout);
648 cout << "===============" << endl;*/
649
650 HypCounter.resize(omp_get_max_threads());
651 for (size_t i = 0; i < HypCounter.size(); ++i)
652 HypCounter[i] = i + 1;
653
654 start_from = ConvHullData.Facets.size();
655 in_triang.resize(start_from, true);
656 in_triang.resize(nr_gen);
657 GensInCone = identity_key(start_from);
658 nrGensInCone = start_from;
659 swap(ConvHullData.Comparisons, Comparisons);
660 Comparisons.resize(start_from);
661 nrTotalComparisons = ConvHullData.nrTotalComparisons;
662 old_nr_supp_hyps = ConvHullData.Generators.nr_of_rows();
663
664 // FACETDATA<Integer> new_facet;
665
666 for (size_t i = 0; i < old_nr_supp_hyps; ++i) {
667 FACETDATA<Integer> new_facet;
668 new_facet.GenInHyp.resize(nr_gen);
669 size_t j = 0;
670 size_t nr_gens_in_fac = 0;
671 for (const auto& Fac : ConvHullData.Facets) {
672 new_facet.GenInHyp[j] = Fac.GenInHyp[i];
673 if (new_facet.GenInHyp[j])
674 nr_gens_in_fac++;
675 j++;
676 }
677 new_facet.simplicial = (nr_gens_in_fac == dim - 1);
678 new_facet.BornAt = 0;
679 new_facet.Mother = 0;
680 // new_facet.is_positive_on_all_original_gens = false;
681 // new_facet.is_negative_on_some_original_gen = false;
682 new_facet.Ident = HypCounter[0];
683 HypCounter[0] += HypCounter.size();
684
685 if (ConvHullData.is_primal)
686 ConvHullData.SLR.convert_to_sublattice(new_facet.Hyp, ConvHullData.Generators[i]);
687 else
688 ConvHullData.SLR.convert_to_sublattice_dual(new_facet.Hyp, ConvHullData.Generators[i]);
689
690 Facets.push_back(new_facet);
691 }
692
693 size_t j = 0;
694 for (const auto& Fac : ConvHullData.Facets) {
695 if (ConvHullData.is_primal) {
696 ConvHullData.SLR.convert_to_sublattice_dual(Generators[j], Fac.Hyp);
697 }
698 else {
699 ConvHullData.SLR.convert_to_sublattice(Generators[j], Fac.Hyp);
700 }
701 ++j;
702 }
703
704 use_existing_facets = true;
705 }
706
707 template <typename Integer>
708 template <typename IntegerCone>
restore_previous_computation(CONVEXHULLDATA<IntegerCone> & ConvHullData,bool goal)709 void Full_Cone<Integer>::restore_previous_computation(CONVEXHULLDATA<IntegerCone>& ConvHullData, bool goal) {
710 // goal=true: to primal, goal=false: to dual
711
712 /* ConvHullData.Generators.pretty_print(cout);
713 cout << "===============" << endl;
714 Generators.pretty_print(cout);
715 cout << "===============" << endl;*/
716
717 if (ConvHullData.is_primal != goal) {
718 dualize_and_restore(ConvHullData);
719 return;
720 }
721
722 swap(ConvHullData.HypCounter, HypCounter);
723 start_from = ConvHullData.Generators.nr_of_rows();
724 /* for(size_t i=0;i<start_from;++i)
725 in_triang[i]=ConvHullData.in_triang[i];*/
726 swap(ConvHullData.in_triang, in_triang);
727 swap(ConvHullData.GensInCone, GensInCone);
728 in_triang.resize(nr_gen);
729 nrGensInCone = ConvHullData.nrGensInCone;
730 swap(ConvHullData.Comparisons, Comparisons);
731 Comparisons.resize(start_from);
732 nrTotalComparisons = ConvHullData.nrTotalComparisons;
733 old_nr_supp_hyps = ConvHullData.old_nr_supp_hyps;
734
735 for (auto& Fac : ConvHullData.Facets) {
736 FACETDATA<Integer> Ret;
737 if (ConvHullData.is_primal)
738 ConvHullData.SLR.convert_to_sublattice_dual(Ret.Hyp, Fac.Hyp);
739 else
740 ConvHullData.SLR.convert_to_sublattice(Ret.Hyp, Fac.Hyp);
741 swap(Ret.GenInHyp, Fac.GenInHyp);
742 Ret.GenInHyp.resize(nr_gen);
743 // convert(Ret.ValNewGen,Fac.ValNewGen);
744 Ret.BornAt = Fac.BornAt;
745 Ret.Ident = Fac.Ident;
746 Ret.Mother = Fac.Mother;
747 // Ret.is_positive_on_all_original_gens = Fac.is_positive_on_all_original_gens;
748 // Ret.is_negative_on_some_original_gen = Fac.is_negative_on_some_original_gen;
749 Ret.simplicial = Fac.simplicial;
750
751 Facets.push_back(Ret);
752 }
753
754 for (size_t i = 0; i < ConvHullData.Generators.nr_of_rows(); ++i) {
755 if (ConvHullData.is_primal)
756 ConvHullData.SLR.convert_to_sublattice(Generators[i], ConvHullData.Generators[i]);
757 else
758 ConvHullData.SLR.convert_to_sublattice_dual(Generators[i], ConvHullData.Generators[i]);
759 }
760
761 use_existing_facets = true;
762 }
763
764 //---------------------------------------------------------------------------
765
766 // Class for the computation of multiplicities via signed decompoasition
767
768 template <typename Integer>
769 class SignedDec {
770
771 template <typename>
772 friend class Full_Cone;
773
774 public:
775
776 bool verbose;
777
778 vector<pair<dynamic_bitset, dynamic_bitset > >* SubfacetsBySimplex;
779 size_t size_hollow_triangulation;
780 size_t dim;
781 size_t nr_gen;
782 int omp_start_level;
783 mpq_class multiplicity;
784 mpz_class int_multiplicity;
785 long decimal_digits;
786 bool approximate;
787
788 mpz_class approx_denominator;
789
790 Integer GradingDenom;
791
792 string Polynomial;
793 // nmz_float EuclideanIntegral;
794 mpq_class Integral, VirtualMultiplicity;
795 nmz_float RawEuclideanIntegral;
796 long DegreeOfPolynomial;
797
798 Matrix<Integer> Generators;
799 Matrix<Integer> Embedding; // transformation on the primal side back to cone coordinates
800 // Matrix<mpz_class> Genererators_mpz;
801 vector<Integer> GradingOnPrimal;
802 // Matrix<mpz_class> GradingOnPrimal_mpz;
803 Matrix<Integer> CandidatesGeneric;
804 vector<Integer> Generic;
805 vector<Integer> GenericComputed;
806
807 void first_subfacet (const dynamic_bitset& Subfacet, const bool compute_multiplicity, Matrix<Integer>& PrimalSimplex,
808 mpz_class& MultPrimal, vector<Integer>& DegreesPrimal, Matrix<Integer>& ValuesGeneric);
809 void next_subfacet(const dynamic_bitset& Subfacet_next, const dynamic_bitset& Subfacet_start,
810 const Matrix<Integer>& PrimalSimplex, const bool compute_multiplicity,
811 const mpz_class& MultPrimal, mpz_class& NewMult,
812 const vector<Integer>& DegreesPrimal, vector<Integer>& NewDegrees,
813 const Matrix<Integer>& ValuesGeneric, Matrix<Integer>& NewValues);
814
815 SignedDec();
816 SignedDec(vector< pair<dynamic_bitset, dynamic_bitset > >& SFS, const Matrix<Integer>& Gens,
817 const vector<Integer> Grad, const int osl);
818 bool FindGeneric();
819 bool ComputeMultiplicity();
820 bool ComputeIntegral(const bool do_virt);
821
822 };
823
824 //---------------------------------------------------------------------------
825
826 } // namespace libnormaliz
827
828 //---------------------------------------------------------------------------
829 #endif
830 //---------------------------------------------------------------------------
831