1 /*
2  * Normaliz
3  * Copyright (C) 2007-2021  W. Bruns, B. Ichim, Ch. Soeger, U. v. d. Ohe
4  * This program is free software: you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 3 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
16  *
17  * As an exception, when this program is distributed through (i) the App Store
18  * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
19  * by Google Inc., then that store may impose any digital rights management,
20  * device limits and/or redistribution restrictions that are required by its
21  * terms of service.
22  */
23 
24 //---------------------------------------------------------------------------
25 
26 // The matrix Map is assumed to be such that the rank of Map equals
27 // the number of columns of Map and no test is performed in the constructor
28 
29 //---------------------------------------------------------------------------
30 #ifndef LIBNORMALIZ_SIMPLEX_H
31 #define LIBNORMALIZ_SIMPLEX_H
32 
33 //---------------------------------------------------------------------------
34 
35 #include <vector>
36 #include <list>
37 
38 #include "libnormaliz/general.h"
39 #include "libnormaliz/cone.h"
40 #include "libnormaliz/HilbertSeries.h"
41 #include "libnormaliz/reduction.h"
42 #include "libnormaliz/dynamic_bitset.h"
43 
44 //---------------------------------------------------------------------------
45 
46 namespace libnormaliz {
47 using std::list;
48 using std::vector;
49 
50 template <typename Integer>
51 class Full_Cone;
52 template <typename Integer>
53 class Collector;
54 
55 /* The SimplexEvaluator provides preallocated space for the comutations in simplex.cpp,
56  * especially matrices and vectors. It also collects elements for thge local Hilbert basis.
57  * The Full_Cone provides one evaluator per thread.
58  */
59 
60 template <typename Integer>
61 class SimplexEvaluator {
62     Full_Cone<Integer>* C_ptr;
63     int tn;  // number of associated thread in parallelization of evaluators
64              // to be set by Full_Cone
65     size_t dim;
66     Integer volume;
67     mpz_class mpz_volume;
68     size_t Deg0_offset;  // the degree of 0+offset
69     long level_offset;   // the same for the inhomogeneous case
70     // Integer det_sum; // sum of the determinants of all evaluated simplices --> Collector
71     // mpq_class mult_sum; // sum of the multiplicities of all evaluated simplices --> Collector
72     vector<key_t> key;
73     // size_t candidates_size;
74     // size_t collected_elements_size;
75     Matrix<Integer> Generators;
76     Matrix<Integer> LinSys;
77     Matrix<Integer> GenCopy;
78     Matrix<Integer> InvGenSelRows;  // selected rows of inverse of Gen
79     Matrix<Integer> InvGenSelCols;  // selected columns of inverse of Gen
80     Matrix<Integer> Sol;
81     Matrix<Integer> ProjGen;  // generators of projection modulo level 0 sublattice
82     vector<Integer> GDiag;    // diagonal of generator matrix after trigonalization
83     vector<Integer> TDiag;    // diagonal of transpose of generaor matrix after trigonalization
84     vector<bool> Excluded;
85     vector<Integer> Indicator;
86     vector<Integer> gen_degrees;
87     vector<long> gen_degrees_long;
88     vector<long> level0_gen_degrees;  // degrees of the generaors in level 0
89     vector<Integer> gen_levels;
90     vector<long> gen_levels_long;
91     // vector< num_t > hvector;  //h-vector of the current evaluation
92     // vector< num_t > inhom_hvector; // separate vector in the inhomogeneous case in case we want to compute two h-vectors
93     // HilbertSeries Hilbert_Series; //this is the summed Hilbert Series
94     // list< vector<Integer> > Candidates;
95     list<vector<Integer> > Hilbert_Basis;
96     // CandidateList<Integer> Collected_HB_Elements;
97     // list< vector<Integer> > Collected_Deg1_Elements;
98     // temporary objects are kept to prevent repeated alloc and dealloc
99     Matrix<Integer> RS;  // right hand side to hold order vector
100     // Matrix<Integer> RSmult; // for multiple right hand sides
101 
102     Matrix<long>* StanleyMat;
103     size_t StanIndex;
104     size_t nr_level0_gens;  // counts the number of level 0 vectors among the generators
105 
106     bool sequential_evaluation;  // indicates whether the simplex is evaluated by a single thread
107 
108     bool GMP_transition;
109 
110     struct SIMPLINEXDATA {         // local data of excluded faces
111         dynamic_bitset GenInFace;  // indicator for generators of simplex in face
112         // vector< num_t > hvector;             // accumulates the h-vector of this face
113         long mult;  // the multiplicity of this face
114         // bool touched;                        // indicates whether hvector != 0 // ALWAYS true, hence superfluous
115         vector<long> gen_degrees;  // degrees of generators in this face
116     };
117     vector<SIMPLINEXDATA> InExSimplData;
118     size_t nrInExSimplData;
119     // bool InExTouched;                        // indicates whether any hvector!=0 // see above
120 
121     vector<vector<Integer>*> RS_pointers;
122     Matrix<Integer> unit_matrix;
123     vector<key_t> id_key;
124     Matrix<mpz_class> mpz_Generators;
125     Integer HB_bound;
126     bool HB_bound_computed;
127 
128     void local_reduction(Collector<Integer>& Coll);
129 
130     // checks whether new_element is reducible by the Reducers list
131     bool is_reducible(const vector<Integer>& new_element, list<vector<Integer> >& Reducers);
132     // removes elements from Candi which are reducible by Reducers, Reducers must be sorted by compare_last!
133     void reduce(list<vector<Integer> >& Candi, list<vector<Integer> >& Reducers, size_t& Candi_size);
134     void count_and_reduce(list<vector<Integer> >& Candi, list<vector<Integer> >& Reducers);
135 
136     bool isDuplicate(const vector<Integer>& cand) const;
137 
138     void addMult(Integer multiplicity, Collector<Integer>& Coll);
139 
140     void prepare_inclusion_exclusion_simpl(size_t Deg, Collector<Integer>& Coll);
141     void add_to_inex_faces(const vector<Integer> offset, size_t Deg, Collector<Integer>& Coll);
142     void update_inhom_hvector(long level_offset, size_t Deg, Collector<Integer>& Coll);
143     void update_mult_inhom(Integer& multiplicity);
144 
145     Integer start_evaluation(SHORTSIMPLEX<Integer>& s, Collector<Integer>& Coll);
146     void find_excluded_facets();
147     void take_care_of_0vector(Collector<Integer>& Coll);
148     // void evaluation_loop_sequential(Collector<Integer>& Coll);
149     void evaluate_element(const vector<Integer>& element, Collector<Integer>& Coll);
150     void conclude_evaluation(Collector<Integer>& Coll);
151     void evaluation_loop_parallel();
152     void evaluate_block(long block_start, long block_end, Collector<Integer>& Coll);
153     void collect_vectors();
154     void add_hvect_to_HS(Collector<Integer>& Coll);
155     void reduce_against_global(Collector<Integer>& Coll);
156 
157     // void insert_gens();
158     // void insert_gens_transpose();
159     // void insert_unit_vectors(vector<key_t> RHS_key);
160 
161     void transform_to_global(const vector<Integer>& element, vector<Integer>& help);
162 
163     //---------------------------------------------------------------------------
164 
165    public:
166     SimplexEvaluator(Full_Cone<Integer>& fc);
167 
168     // sets the thread number of the evaluator (needed to associate a collector)
169     void set_evaluator_tn(int threadnum);
170 
171     // full evaluation of the simplex in a single thread, delivers results to to a collector
172     bool evaluate(SHORTSIMPLEX<Integer>& s);
173 
174     // evaluation in parallel threads
175     void Simplex_parallel_evaluation();
176 
177     vector<key_t> get_key();
178     Integer get_volume();
179 
180     void print_all();
181 };
182 // class SimplexEvaluator end
183 
184 /* The collector collects the results of the computations in simplex.cpp.
185  * The Full_Cpne provides one collector per thread.
186  *
187  */
188 
189 template <typename Integer>
190 class Collector {
191     template <typename>
192     friend class SimplexEvaluator;
193     template <typename>
194     friend class Full_Cone;
195 
196     Full_Cone<Integer>* C_ptr;
197     size_t dim;
198 
199     Integer det_sum;     // sum of the determinants of all evaluated simplices
200     mpq_class mult_sum;  // sum of the multiplicities of all evaluated simplices
201     size_t candidates_size;
202     size_t collected_elements_size;
203     vector<num_t> hvector;         // h-vector of the current evaluation
204     vector<num_t> inhom_hvector;   // separate vector in the inhomogeneous case in case we want to compute two h-vectors
205     HilbertSeries Hilbert_Series;  // this is the summed Hilbert Series
206     list<vector<Integer> > Candidates;
207     CandidateList<Integer> HB_Elements;
208     list<vector<Integer> > Deg1_Elements;
209     vector<vector<num_t> > InEx_hvector;
210 
211     Matrix<Integer> elements;
212 
213    public:
214     Collector(Full_Cone<Integer>& fc);
215 
216     // returns sum of the determinants of all evaluated simplices
217     Integer getDetSum() const;
218 
219     // returns sum of the multiplicities of all evaluated simplices
220     mpq_class getMultiplicitySum() const;
221 
222     // returns sum of the Hilbert Series of all evaluated simplices
223     const HilbertSeries& getHilbertSeriesSum() const;
224 
225     // moves the union of Hilbert basis / deg1 elements to the cone
226     // for partial triangulation it merges the sorted list
227     void transfer_candidates();
228 
229     size_t get_collected_elements_size();
230 };
231 // class end Collector
232 
233 }  // namespace libnormaliz
234 
235 //---------------------------------------------------------------------------
236 #endif
237 //---------------------------------------------------------------------------
238