CoCalc -- cone.cpp

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / NormalizInterface-1.0.2 / Normaliz.git / source / libnormaliz / cone.cpp
⁵⁶³⁶²⁴ views
1
/*
2
 * Normaliz
3
 * Copyright (C) 2007-2014  Winfried Bruns, Bogdan Ichim, Christof Soeger
4
 * This program is free software: you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation, either version 3 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
 *
17
 * As an exception, when this program is distributed through (i) the App Store
18
 * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
19
 * by Google Inc., then that store may impose any digital rights management,
20
 * device limits and/or redistribution restrictions that are required by its
21
 * terms of service.
22
 */
23

24
#include <stdlib.h>
25
#include <list>
26
#include <sys/stat.h>
27
#include <sys/types.h>
28

29
#include "libnormaliz/vector_operations.h"
30
#include "libnormaliz/map_operations.h"
31
#include "libnormaliz/convert.h"
32
#include "libnormaliz/cone.h"
33
#include "libnormaliz/full_cone.h"
34
#include "libnormaliz/my_omp.h"
35

36
namespace libnormaliz {
37
using namespace std;
38

39
// adds the signs inequalities given by Signs to Inequalities
40
template<typename Integer>
41
Matrix<Integer> sign_inequalities(const vector< vector<Integer> >& Signs) {
42
    if (Signs.size() != 1) {
43
        throw BadInputException("ERROR: Bad signs matrix, has "
44
                + toString(Signs.size()) + " rows (should be 1)!");
45
    }
46
    size_t dim = Signs[0].size();
47
    Matrix<Integer> Inequ(0,dim);
48
    vector<Integer> ineq(dim,0);
49
    for (size_t i=0; i<dim; i++) {
50
        Integer sign = Signs[0][i];
51
        if (sign == 1 || sign == -1) {
52
            ineq[i] = sign;
53
            Inequ.append(ineq);
54
            ineq[i] = 0;
55
        } else if (sign != 0) {
56
            throw BadInputException("Bad signs matrix, has entry "
57
                    + toString(sign) + " (should be -1, 1 or 0)!");
58
        }
59
    }
60
    return Inequ;
61
}
62

63
template<typename Integer>
64
Matrix<Integer> strict_sign_inequalities(const vector< vector<Integer> >& Signs) {
65
    if (Signs.size() != 1) {
66
        throw BadInputException("ERROR: Bad signs matrix, has "
67
                + toString(Signs.size()) + " rows (should be 1)!");
68
    }
69
    size_t dim = Signs[0].size();
70
    Matrix<Integer> Inequ(0,dim);
71
    vector<Integer> ineq(dim,0);
72
    ineq[dim-1]=-1;
73
    for (size_t i=0; i<dim-1; i++) {    // last component of strict_signs always 0
74
        Integer sign = Signs[0][i];
75
        if (sign == 1 || sign == -1) {
76
            ineq[i] = sign;
77
            Inequ.append(ineq);
78
            ineq[i] = 0;
79
        } else if (sign != 0) {
80
            throw BadInputException("Bad signs matrix, has entry "
81
                    + toString(sign) + " (should be -1, 1 or 0)!");
82
        }
83
    }
84
    return Inequ;
85
}
86

87
template<typename Integer>
88
vector<vector<Integer> > find_input_matrix(const map< InputType, vector< vector<Integer> > >& multi_input_data,
89
                               const InputType type){
90

91
    typename map< InputType , vector< vector<Integer> > >::const_iterator it;
92
    it = multi_input_data.find(type);
93
    if (it != multi_input_data.end())
94
        return(it->second);
95

96
     vector< vector<Integer> > dummy;
97
     return(dummy);
98
}
99

100
template<typename Integer>
101
void insert_column(vector< vector<Integer> >& mat, size_t col, Integer entry){
102

103
    if(mat.size()==0)
104
        return;
105
    vector<Integer> help(mat[0].size()+1);
106
    for(size_t i=0;i<mat.size();++i){
107
        for(size_t j=0;j<col;++j)
108
            help[j]=mat[i][j];
109
        help[col]=entry;
110
        for(size_t j=col;j<mat[i].size();++j)
111
            help[j+1]=mat[i][j];
112
        mat[i]=help;
113
    }
114
}
115

116
template<typename Integer>
117
void insert_zero_column(vector< vector<Integer> >& mat, size_t col){
118
    // Integer entry=0;
119
    insert_column<Integer>(mat,col,0);
120
}
121

122
template<typename Integer>
123
void Cone<Integer>::homogenize_input(map< InputType, vector< vector<Integer> > >& multi_input_data){
124

125
    typename map< InputType , vector< vector<Integer> > >::iterator it;
126
    it = multi_input_data.begin();
127
    for(;it!=multi_input_data.end();++it){
128
        switch(it->first){
129
            case Type::dehomogenization:
130
                throw BadInputException("Type dehomogenization not allowed with inhomogeneous input!");
131
                break;
132
            case Type::inhom_inequalities: // nothing to do
133
            case Type::inhom_equations:
134
            case Type::inhom_congruences:
135
            case Type::polyhedron:
136
            case Type::vertices:
137
            case Type::support_hyperplanes:
138
            case Type::open_facets:
139
            case Type::grading:  // already taken care of
140
                break;
141
            case Type::strict_inequalities:
142
                insert_column<Integer>(it->second,dim-1,-1);
143
                break;
144
            case Type::offset:
145
                insert_column<Integer>(it->second,dim-1,1);
146
                break;
147
            default:  // is correct for signs and strict_signs !
148
                insert_zero_column<Integer>(it->second,dim-1);
149
                break;
150
        }
151
    }
152
}
153

154
bool denominator_allowed(InputType input_type){
155
    
156
    switch(input_type){
157
        
158
        case Type::congruences:
159
        case Type::inhom_congruences:
160
        case Type::grading:
161
        case Type::dehomogenization:
162
        case Type::lattice:
163
        case Type::normalization:
164
        case Type::cone_and_lattice:
165
        case Type::offset:
166
        case Type::rees_algebra:
167
        case Type::lattice_ideal:
168
        case Type::signs:
169
        case Type::strict_signs:
170
//         case Type::open_facets:
171
            return false;
172
            break;
173
        default:
174
            return true;
175
        break;
176
    }
177
}
178

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

181
template<typename Integer>
182
Cone<Integer>::Cone(){
183
}
184

185
template<typename Integer>
186
Cone<Integer>::Cone(InputType input_type, const vector< vector<Integer> >& Input) {
187
    // convert to a map
188
    map< InputType, vector< vector<Integer> > > multi_input_data;
189
    multi_input_data[input_type] = Input;
190
    process_multi_input(multi_input_data);
191
}
192

193
template<typename Integer>
194
Cone<Integer>::Cone(InputType type1, const vector< vector<Integer> >& Input1,
195
                    InputType type2, const vector< vector<Integer> >& Input2) {
196
    if (type1 == type2) {
197
        throw BadInputException("Input types must  pairwise different!");
198
    }
199
    // convert to a map
200
    map< InputType, vector< vector<Integer> > > multi_input_data;
201
    multi_input_data[type1] = Input1;
202
    multi_input_data[type2] = Input2;
203
    process_multi_input(multi_input_data);
204
}
205

206
template<typename Integer>
207
Cone<Integer>::Cone(InputType type1, const vector< vector<Integer> >& Input1,
208
                    InputType type2, const vector< vector<Integer> >& Input2,
209
                    InputType type3, const vector< vector<Integer> >& Input3) {
210
    if (type1 == type2 || type1 == type3 || type2 == type3) {
211
        throw BadInputException("Input types must be pairwise different!");
212
    }
213
    // convert to a map
214
    map< InputType, vector< vector<Integer> > > multi_input_data;
215
    multi_input_data[type1] = Input1;
216
    multi_input_data[type2] = Input2;
217
    multi_input_data[type3] = Input3;
218
    process_multi_input(multi_input_data);
219
}
220

221
template<typename Integer>
222
Cone<Integer>::Cone(const map< InputType, vector< vector<Integer> > >& multi_input_data) {
223
    process_multi_input(multi_input_data);
224
}
225

226
// now with mpq_class input
227

228
template<typename Integer>
229
Cone<Integer>::Cone(InputType input_type, const vector< vector<mpq_class> >& Input) {
230
    // convert to a map
231
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
232
    multi_input_data[input_type] = Input;
233
    process_multi_input(multi_input_data);
234
}
235

236
template<typename Integer>
237
Cone<Integer>::Cone(InputType type1, const vector< vector<mpq_class> >& Input1,
238
                    InputType type2, const vector< vector<mpq_class> >& Input2) {
239
    if (type1 == type2) {
240
        throw BadInputException("Input types must  pairwise different!");
241
    }
242
    // convert to a map
243
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
244
    multi_input_data[type1] = Input1;
245
    multi_input_data[type2] = Input2;
246
    process_multi_input(multi_input_data);
247
}
248

249
template<typename Integer>
250
Cone<Integer>::Cone(InputType type1, const vector< vector<mpq_class> >& Input1,
251
                    InputType type2, const vector< vector<mpq_class> >& Input2,
252
                    InputType type3, const vector< vector<mpq_class> >& Input3) {
253
    if (type1 == type2 || type1 == type3 || type2 == type3) {
254
        throw BadInputException("Input types must be pairwise different!");
255
    }
256
    // convert to a map
257
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
258
    multi_input_data[type1] = Input1;
259
    multi_input_data[type2] = Input2;
260
    multi_input_data[type3] = Input3;
261
    process_multi_input(multi_input_data);
262
}
263

264
template<typename Integer>
265
Cone<Integer>::Cone(const map< InputType, vector< vector<mpq_class> > >& multi_input_data) {
266
    process_multi_input(multi_input_data);
267
}
268

269
// now with nmz_float input
270

271
template<typename Integer>
272
Cone<Integer>::Cone(InputType input_type, const vector< vector<nmz_float> >& Input) {
273
    // convert to a map
274
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
275
    multi_input_data[input_type] = Input;
276
    process_multi_input(multi_input_data);
277
}
278

279
template<typename Integer>
280
Cone<Integer>::Cone(InputType type1, const vector< vector<nmz_float> >& Input1,
281
                    InputType type2, const vector< vector<nmz_float> >& Input2) {
282
    if (type1 == type2) {
283
        throw BadInputException("Input types must  pairwise different!");
284
    }
285
    // convert to a map
286
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
287
    multi_input_data[type1] = Input1;
288
    multi_input_data[type2] = Input2;
289
    process_multi_input(multi_input_data);
290
}
291

292
template<typename Integer>
293
Cone<Integer>::Cone(InputType type1, const vector< vector<nmz_float> >& Input1,
294
                    InputType type2, const vector< vector<nmz_float> >& Input2,
295
                    InputType type3, const vector< vector<nmz_float> >& Input3) {
296
    if (type1 == type2 || type1 == type3 || type2 == type3) {
297
        throw BadInputException("Input types must be pairwise different!");
298
    }
299
    // convert to a map
300
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
301
    multi_input_data[type1] = Input1;
302
    multi_input_data[type2] = Input2;
303
    multi_input_data[type3] = Input3;
304
    process_multi_input(multi_input_data);
305
}
306

307
template<typename Integer>
308
Cone<Integer>::Cone(const map< InputType, vector< vector<nmz_float> > >& multi_input_data) {
309
    process_multi_input(multi_input_data);
310
}
311

312
//---------------------------------------------------------------------------
313
// now with Matrix
314
//---------------------------------------------------------------------------
315

316
template<typename Integer>
317
Cone<Integer>::Cone(InputType input_type, const Matrix<Integer>& Input) {
318
    // convert to a map
319
    map< InputType, vector< vector<Integer> > >multi_input_data;
320
    multi_input_data[input_type] = Input.get_elements();
321
    process_multi_input(multi_input_data);
322
}
323

324
template<typename Integer>
325
Cone<Integer>::Cone(InputType type1, const Matrix<Integer>& Input1,
326
                    InputType type2, const Matrix<Integer>& Input2) {
327
    if (type1 == type2) {
328
        throw BadInputException("Input types must  pairwise different!");
329
    }
330
    // convert to a map
331
    map< InputType, vector< vector<Integer> > > multi_input_data;
332
    multi_input_data[type1] = Input1.get_elements();
333
    multi_input_data[type2] = Input2.get_elements();
334
    process_multi_input(multi_input_data);
335
}
336

337
template<typename Integer>
338
Cone<Integer>::Cone(InputType type1, const Matrix<Integer>& Input1,
339
                    InputType type2, const Matrix<Integer>& Input2,
340
                    InputType type3, const Matrix<Integer>& Input3) {
341
    if (type1 == type2 || type1 == type3 || type2 == type3) {
342
        throw BadInputException("Input types must be pairwise different!");
343
    }
344
    // convert to a map
345
    map< InputType, vector< vector<Integer> > > multi_input_data;
346
    multi_input_data[type1] = Input1.get_elements();
347
    multi_input_data[type2] = Input2.get_elements();
348
    multi_input_data[type3] = Input3.get_elements();
349
    process_multi_input(multi_input_data);
350
}
351

352
template<typename Integer>
353
Cone<Integer>::Cone(const map< InputType, Matrix<Integer> >& multi_input_data_Matrix){
354
    map< InputType, vector< vector<Integer> > > multi_input_data;
355
    auto it = multi_input_data_Matrix.begin();
356
    for(; it != multi_input_data_Matrix.end(); ++it){
357
        multi_input_data[it->first]=it->second.get_elements();
358
    }
359
    process_multi_input(multi_input_data);
360
}
361

362
//---------------------------------------------------------------------------
363
// now with Matrix and mpq_class
364

365
template<typename Integer>
366
Cone<Integer>::Cone(InputType input_type, const Matrix<mpq_class>& Input) {
367
    // convert to a map
368
    map< InputType, vector< vector<mpq_class> > >multi_input_data;
369
    multi_input_data[input_type] = Input.get_elements();
370
    process_multi_input(multi_input_data);
371
}
372

373
template<typename Integer>
374
Cone<Integer>::Cone(InputType type1, const Matrix<mpq_class>& Input1,
375
                    InputType type2, const Matrix<mpq_class>& Input2) {
376
    if (type1 == type2) {
377
        throw BadInputException("Input types must  pairwise different!");
378
    }
379
    // convert to a map
380
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
381
    multi_input_data[type1] = Input1.get_elements();
382
    multi_input_data[type2] = Input2.get_elements();
383
    process_multi_input(multi_input_data);
384
}
385

386
template<typename Integer>
387
Cone<Integer>::Cone(InputType type1, const Matrix<mpq_class>& Input1,
388
                    InputType type2, const Matrix<mpq_class>& Input2,
389
                    InputType type3, const Matrix<mpq_class>& Input3) {
390
    if (type1 == type2 || type1 == type3 || type2 == type3) {
391
        throw BadInputException("Input types must be pairwise different!");
392
    }
393
    // convert to a map
394
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
395
    multi_input_data[type1] = Input1.get_elements();
396
    multi_input_data[type2] = Input2.get_elements();
397
    multi_input_data[type3] = Input3.get_elements();
398
    process_multi_input(multi_input_data);
399
}
400

401
template<typename Integer>
402
Cone<Integer>::Cone(const map< InputType, Matrix<mpq_class> >& multi_input_data_Matrix){
403
    map< InputType, vector< vector<mpq_class> > > multi_input_data;
404
    auto it = multi_input_data_Matrix.begin();
405
    for(; it != multi_input_data_Matrix.end(); ++it){
406
        multi_input_data[it->first]=it->second.get_elements();
407
    }
408
    process_multi_input(multi_input_data);
409
}
410

411
//---------------------------------------------------------------------------
412
// now with Matrix and nmz_float
413

414
template<typename Integer>
415
Cone<Integer>::Cone(InputType input_type, const Matrix<nmz_float>& Input) {
416
    // convert to a map
417
    map< InputType, vector< vector<nmz_float> > >multi_input_data;
418
    multi_input_data[input_type] = Input.get_elements();
419
    process_multi_input(multi_input_data);
420
}
421

422
template<typename Integer>
423
Cone<Integer>::Cone(InputType type1, const Matrix<nmz_float>& Input1,
424
                    InputType type2, const Matrix<nmz_float>& Input2) {
425
    if (type1 == type2) {
426
        throw BadInputException("Input types must  pairwise different!");
427
    }
428
    // convert to a map
429
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
430
    multi_input_data[type1] = Input1.get_elements();
431
    multi_input_data[type2] = Input2.get_elements();
432
    process_multi_input(multi_input_data);
433
}
434

435
template<typename Integer>
436
Cone<Integer>::Cone(InputType type1, const Matrix<nmz_float>& Input1,
437
                    InputType type2, const Matrix<nmz_float>& Input2,
438
                    InputType type3, const Matrix<nmz_float>& Input3) {
439
    if (type1 == type2 || type1 == type3 || type2 == type3) {
440
        throw BadInputException("Input types must be pairwise different!");
441
    }
442
    // convert to a map
443
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
444
    multi_input_data[type1] = Input1.get_elements();
445
    multi_input_data[type2] = Input2.get_elements();
446
    multi_input_data[type3] = Input3.get_elements();
447
    process_multi_input(multi_input_data);
448
}
449

450
template<typename Integer>
451
Cone<Integer>::Cone(const map< InputType, Matrix<nmz_float> >& multi_input_data_Matrix){
452
    map< InputType, vector< vector<nmz_float> > > multi_input_data;
453
    auto it = multi_input_data_Matrix.begin();
454
    for(; it != multi_input_data_Matrix.end(); ++it){
455
        multi_input_data[it->first]=it->second.get_elements();
456
    }
457
    process_multi_input(multi_input_data);
458
}
459

460
//---------------------------------------------------------------------------
461

462
template<typename Integer>
463
Cone<Integer>::~Cone() {
464
    if(IntHullCone!=NULL)
465
        delete IntHullCone;
466
    if(IntHullCone!=NULL)
467
        delete SymmCone;
468
}
469

470
//---------------------------------------------------------------------------
471

472
template<typename Integer>
473
void Cone<Integer>::process_multi_input(const map< InputType, vector< vector<mpq_class> > >& multi_input_data_const) {
474

475
    map< InputType, vector< vector<mpq_class> > > multi_input_data(multi_input_data_const);    
476
    // since polytope will be comverted to cone, we must do some checks here
477
    if(exists_element(multi_input_data,Type::grading) && exists_element(multi_input_data,Type::polytope)){
478
           throw BadInputException("No explicit grading allowed with polytope!");
479
    }
480
    if(exists_element(multi_input_data,Type::cone) && exists_element(multi_input_data,Type::polytope)){
481
        throw BadInputException("Illegal combination of cone generator types!");
482
    }
483
    
484
    map< InputType, vector< vector<Integer> > > multi_input_data_ZZ;
485
    
486
    // special treatment of polytope. We convert it o cone
487
    if(exists_element(multi_input_data,Type::polytope)){
488
        size_t dim;
489
        if(multi_input_data[Type::polytope].size()>0){
490
            dim=multi_input_data[Type::polytope][0].size()+1;
491
            vector<vector<Integer> > grading;
492
            grading.push_back(vector<Integer>(dim));
493
            grading[0][dim-1]=1;
494
            multi_input_data_ZZ[Type::grading]=grading;
495
        }
496
        multi_input_data[Type::cone]=multi_input_data[Type::polytope];
497
        multi_input_data.erase(Type::polytope);
498
        for(size_t i=0;i<multi_input_data[Type::cone].size();++i){
499
            multi_input_data[Type::cone][i].resize(dim);
500
            multi_input_data[Type::cone][i][dim-1]=1;
501
        }
502
    }
503
    
504
    // now we clear denominators
505
    auto it = multi_input_data.begin();
506
    for(; it != multi_input_data.end(); ++it) {
507
        for(size_t i=0;i < it->second.size();++i){ 
508
            mpz_class common_denom=1;
509
            for(size_t j=0;j<it->second[i].size();++j){
510
                it->second[i][j].canonicalize();
511
                common_denom=libnormaliz::lcm(common_denom,it->second[i][j].get_den());
512
            }
513
            if(common_denom>1 && !denominator_allowed(it->first))
514
                throw BadInputException("Proper fraction not allowed in certain input types");
515
            vector<Integer> transfer(it->second[i].size());
516
            for(size_t j=0;j<it->second[i].size();++j){
517
                it->second[i][j]*=common_denom;
518
                convert(transfer[j],it->second[i][j].get_num());
519
            }
520
            multi_input_data_ZZ[it->first].push_back(transfer);
521
        }
522
    }
523
    
524
    process_multi_input_inner(multi_input_data_ZZ);
525
}
526

527
template<typename Integer>
528
void Cone<Integer>::process_multi_input(const map< InputType, vector< vector<nmz_float> > >& multi_input_data) {
529
    
530
    map< InputType, vector< vector<mpq_class> > > multi_input_data_QQ;
531
    auto it = multi_input_data.begin();
532
    for(; it != multi_input_data.end(); ++it) {
533
        vector<vector<mpq_class> > Transfer;
534
        vector<mpq_class> vt;
535
        for(size_t j=0;j<it->second.size();++j){
536
            for(size_t k=0;k<it->second[j].size();++k)
537
                vt.push_back(mpq_class(it->second[j][k]));
538
            Transfer.push_back(vt);
539
        }
540
        multi_input_data_QQ[it->first]=Transfer;
541
    }
542
    process_multi_input(multi_input_data_QQ);
543
}
544

545
template<typename Integer>
546
void Cone<Integer>::process_multi_input(const map< InputType, vector< vector<Integer> > >& multi_input_data_const) {
547
    initialize();
548
    map< InputType, vector< vector<Integer> > > multi_input_data(multi_input_data_const);
549
    process_multi_input_inner(multi_input_data);
550
}
551

552
template<typename Integer>
553
void Cone<Integer>::process_multi_input_inner(map< InputType, vector< vector<Integer> > >& multi_input_data) {
554
    initialize();
555
    // find basic input type
556
    lattice_ideal_input=false;
557
    nr_latt_gen=0, nr_cone_gen=0;
558
    bool inhom_input=false;
559
    
560
    inequalities_present=false; //control choice of positive orthant
561

562
    // NEW: Empty matrix have syntactical influence
563
    auto it = multi_input_data.begin();
564
    for(; it != multi_input_data.end(); ++it) {
565
        switch (it->first) {
566
            case Type::inhom_inequalities:
567
            case Type::inhom_equations:
568
            case Type::inhom_congruences:
569
            case Type::strict_inequalities:
570
            case Type::strict_signs:
571
            case Type::open_facets:
572
                inhom_input=true;
573
            case Type::signs:
574
            case Type::inequalities:
575
            case Type::equations:
576
            case Type::congruences:
577
                break;
578
            case Type::lattice_ideal:
579
                lattice_ideal_input=true;
580
                break;
581
            case Type::polyhedron:
582
                inhom_input=true;
583
            case Type::integral_closure:
584
            case Type::rees_algebra:
585
            case Type::polytope:
586
            case Type::cone:
587
            case Type::subspace:
588
                nr_cone_gen++;
589
                break;
590
            case Type::normalization:
591
            case Type::cone_and_lattice:
592
                nr_cone_gen++;
593
            case Type::lattice:
594
            case Type::saturation:
595
                nr_latt_gen++;
596
                break;
597
            case Type::vertices:
598
            case Type::offset:
599
                inhom_input=true;
600
            default:
601
                break;
602
        }
603

604
        switch (it->first) {  // chceck existence of inrqualities
605
            case Type::inhom_inequalities:
606
            case Type::strict_inequalities:
607
            case Type::strict_signs:
608
            case Type::signs:
609
            case Type::inequalities:
610
            case Type::excluded_faces:
611
            case Type::support_hyperplanes:
612
                inequalities_present=true;
613
            default:
614
                break;
615
        }
616

617
    }
618

619
    bool gen_error=false;
620
    if(nr_cone_gen>2)
621
        gen_error=true;
622

623
    if(nr_cone_gen==2 && (!exists_element(multi_input_data,Type::subspace)
624
                      || !(exists_element(multi_input_data,Type::cone)
625
                          || exists_element(multi_input_data,Type::cone_and_lattice)
626
                          || exists_element(multi_input_data,Type::integral_closure)
627
                          || exists_element(multi_input_data,Type::normalization) ) )
628
    )
629
        gen_error=true;
630
    
631
    if(gen_error){
632
        throw BadInputException("Illegal combination of cone generator types!");
633
    }
634
    
635
    
636
    if(nr_latt_gen>1){
637
        throw BadInputException("Only one matrix of lattice generators allowed!");
638
    }
639
    if(lattice_ideal_input){
640
        if(multi_input_data.size() > 2 || (multi_input_data.size()==2 && !exists_element(multi_input_data,Type::grading))){
641
            throw BadInputException("Only grading allowed with lattice_ideal!");
642
        }
643
    }
644
    if(exists_element(multi_input_data,Type::open_facets)){
645
        size_t allowed=0;
646
        auto it = multi_input_data.begin();
647
        for(; it != multi_input_data.end(); ++it) {
648
            switch (it->first) {
649
                case Type::open_facets:
650
                case Type::cone:
651
                case Type::grading:
652
                case Type::vertices:
653
                    allowed++;
654
                    break;
655
                default:
656
                    break;                
657
            }
658
        }
659
        if(allowed!=multi_input_data.size())
660
            throw BadInputException("Illegal combination of input types with open_facets!");        
661
        if(exists_element(multi_input_data,Type::vertices)){
662
            if(multi_input_data[Type::vertices].size()>1)
663
                throw BadInputException("At most one vertex allowed with open_facets!");            
664
        }
665
        
666
    }
667
    if(inhom_input){
668
        if(exists_element(multi_input_data,Type::dehomogenization) || exists_element(multi_input_data,Type::support_hyperplanes)){
669
            throw BadInputException("Types dehomogenization and support_hyperplanes not allowed with inhomogeneous input!");
670
        }
671
    }
672
    if(inhom_input || exists_element(multi_input_data,Type::dehomogenization)){
673
        if(exists_element(multi_input_data,Type::rees_algebra) || exists_element(multi_input_data,Type::polytope)){
674
            throw BadInputException("Types polytope and rees_algebra not allowed with inhomogeneous input or dehomogenization!");
675
        }
676
        if(exists_element(multi_input_data,Type::excluded_faces)){
677
            throw BadInputException("Type excluded_faces not allowed with inhomogeneous input or dehomogenization!");
678
        }
679
    }
680
    if(exists_element(multi_input_data,Type::grading) && exists_element(multi_input_data,Type::polytope)){
681
           throw BadInputException("No explicit grading allowed with polytope!");
682
    }
683
    
684
    // remove empty matrices
685
    it = multi_input_data.begin();
686
    for(; it != multi_input_data.end();) {
687
        if (it->second.size() == 0)
688
            multi_input_data.erase(it++);
689
        else
690
            ++it;
691
    }
692

693
    if(multi_input_data.size()==0){
694
        throw BadInputException("All input matrices empty!");
695
    }
696

697
    //determine dimension
698
    it = multi_input_data.begin();
699
    size_t inhom_corr = 0; // correction in the inhom_input case
700
    if (inhom_input) inhom_corr = 1;
701
    dim = it->second.front().size() - type_nr_columns_correction(it->first) + inhom_corr;
702

703
    // We now process input types that are independent of generators, constraints, lattice_ideal
704
    // check for excluded faces
705
    
706
    ExcludedFaces = find_input_matrix(multi_input_data,Type::excluded_faces);
707
    if(ExcludedFaces.nr_of_rows()==0)
708
        ExcludedFaces=Matrix<Integer>(0,dim); // we may need the correct number of columns
709
    PreComputedSupportHyperplanes = find_input_matrix(multi_input_data,Type::support_hyperplanes);
710
    
711
    // check for a grading
712
    vector< vector<Integer> > lf = find_input_matrix(multi_input_data,Type::grading);
713
    if (lf.size() > 1) {
714
        throw BadInputException("Bad grading, has "
715
                + toString(lf.size()) + " rows (should be 1)!");
716
    }
717
    if(lf.size()==1){
718
        if(inhom_input)
719
            lf[0].push_back(0); // first we extend grading trivially to have the right dimension
720
        setGrading (lf[0]);     // will eventually be set in full_cone.cpp
721

722
    }
723

724
    // cout << "Dim " << dim <<endl;
725

726
    // check consistence of dimension
727
    it = multi_input_data.begin();
728
    size_t test_dim;
729
    for (; it != multi_input_data.end(); ++it) {
730
        test_dim = it->second.front().size() - type_nr_columns_correction(it->first) + inhom_corr;
731
        if (test_dim != dim) {
732
            throw BadInputException("Inconsistent dimensions in input!");
733
        }
734
    }
735

736
    if(inhom_input)
737
        homogenize_input(multi_input_data);
738
    
739
    // check for dehomogenization
740
    lf = find_input_matrix(multi_input_data,Type::dehomogenization);
741
    if (lf.size() > 1) {
742
        throw BadInputException("Bad dehomogenization, has "
743
                + toString(lf.size()) + " rows (should be 1)!");
744
    }
745
    if(lf.size()==1){
746
        setDehomogenization(lf[0]);
747
    }
748

749
    // now we can unify implicit and explicit truncation
750
    // Note: implicit and explicit truncation have already been excluded
751
    if (inhom_input) {
752
        Dehomogenization.resize(dim,0),
753
        Dehomogenization[dim-1]=1;
754
        is_Computed.set(ConeProperty::Dehomogenization);
755
    }
756
    if(isComputed(ConeProperty::Dehomogenization))
757
        inhomogeneous=true;
758

759
    if(lattice_ideal_input){
760
        prepare_input_lattice_ideal(multi_input_data);
761
    }
762

763
    Matrix<Integer> LatticeGenerators(0,dim);
764
    prepare_input_generators(multi_input_data, LatticeGenerators);
765

766
    prepare_input_constraints(multi_input_data); // sets Equations,Congruences,Inequalities
767

768
    // set default values if necessary
769
    if(inhom_input && LatticeGenerators.nr_of_rows()!=0 && !exists_element(multi_input_data,Type::offset)){
770
        vector<Integer> offset(dim);
771
        offset[dim-1]=1;
772
        LatticeGenerators.append(offset);
773
    }
774
    if(inhom_input &&  Generators.nr_of_rows()!=0 && !exists_element(multi_input_data,Type::vertices) 
775
                && !exists_element(multi_input_data,Type::polyhedron)){
776
        vector<Integer> vertex(dim);
777
        vertex[dim-1]=1;
778
        Generators.append(vertex);
779
    }
780

781
    if(Inequalities.nr_of_rows()>0 && Generators.nr_of_rows()>0){ // eliminate superfluous inequalities
782
        vector<key_t> essential;
783
        for(size_t i=0;i<Inequalities.nr_of_rows();++i){
784
            for (size_t j=0;j<Generators.nr_of_rows();++j){
785
                if(v_scalar_product(Inequalities[i],Generators[j])<0){
786
                    essential.push_back(i);
787
                    break;
788
                }
789
            }
790
        }
791
        if(essential.size()<Inequalities.nr_of_rows())
792
            Inequalities=Inequalities.submatrix(essential);
793
    }
794

795
    // cout << "Ineq " << Inequalities.nr_of_rows() << endl;
796

797
    process_lattice_data(LatticeGenerators,Congruences,Equations);
798

799
    bool cone_sat_eq=no_lattice_restriction;
800
    bool cone_sat_cong=no_lattice_restriction;
801

802
    // cout << "nolatrest " << no_lattice_restriction << endl;
803

804
    if(Inequalities.nr_of_rows()==0 && Generators.nr_of_rows()!=0){
805
        if(!no_lattice_restriction){
806
            cone_sat_eq=true;
807
            for(size_t i=0;i<Generators.nr_of_rows() && cone_sat_eq;++i)
808
                for(size_t j=0;j<Equations.nr_of_rows()  && cone_sat_eq ;++j)
809
                    if(v_scalar_product(Generators[i],Equations[j])!=0){
810
                        cone_sat_eq=false;
811
            }
812
        }
813
        if(!no_lattice_restriction){
814
            cone_sat_cong=true;
815
            for(size_t i=0;i<Generators.nr_of_rows() && cone_sat_cong;++i){
816
                vector<Integer> test=Generators[i];
817
                test.resize(dim+1);
818
                for(size_t j=0;j<Congruences.nr_of_rows()  && cone_sat_cong ;++j)
819
                    if(v_scalar_product(test,Congruences[j]) % Congruences[j][dim] !=0)
820
                        cone_sat_cong=false;
821
            }
822
        }
823

824
        if(cone_sat_eq && cone_sat_cong){
825
            set_original_monoid_generators(Generators);
826
        }
827

828
        if(cone_sat_eq && !cone_sat_cong){ // multiply generators by anniullator mod sublattice
829
            for(size_t i=0;i<Generators.nr_of_rows();++i)
830
                v_scalar_multiplication(Generators[i],BasisChange.getAnnihilator());
831
            cone_sat_cong=true;
832
        }
833
    }
834

835
    if((Inequalities.nr_of_rows()!=0 || !cone_sat_eq) && Generators.nr_of_rows()!=0){
836
        Sublattice_Representation<Integer> ConeLatt(Generators,true);
837
        Full_Cone<Integer> TmpCone(ConeLatt.to_sublattice(Generators));
838
        TmpCone.dualize_cone();
839
        Inequalities.append(ConeLatt.from_sublattice_dual(TmpCone.Support_Hyperplanes));
840
        Generators=Matrix<Integer>(0,dim); // Generators now converted into inequalities
841
    }
842

843
    if(exists_element(multi_input_data,Type::open_facets)){
844
        // read manual for the computation that follows
845
        if(!isComputed(ConeProperty::OriginalMonoidGenerators)) // practically impossible, but better to check
846
            throw BadInputException("Error in connection with open_facets");
847
        if(Generators.nr_of_rows()!=BasisChange.getRank())
848
            throw BadInputException("Cone for open_facets not simplicial!");
849
        Matrix<Integer> TransformedGen=BasisChange.to_sublattice(Generators);
850
        vector<key_t> key(TransformedGen.nr_of_rows());
851
        for(size_t j=0;j<TransformedGen.nr_of_rows();++j)
852
            key[j]=j;
853
        Matrix<Integer> TransformedSupps;
854
        Integer dummy;
855
        TransformedGen.simplex_data(key,TransformedSupps,dummy,false);
856
        Matrix<Integer> NewSupps=BasisChange.from_sublattice_dual(TransformedSupps);
857
        NewSupps.remove_row(NewSupps.nr_of_rows()-1); // must remove the inequality for the homogenizing variable
858
        for(size_t j=0;j<NewSupps.nr_of_rows();++j){
859
            if(!(multi_input_data[Type::open_facets][0][j]==0 || multi_input_data[Type::open_facets][0][j]==1))
860
                throw BadInputException("Illegal entry in open_facets");
861
            NewSupps[j][dim-1]-=multi_input_data[Type::open_facets][0][j];
862
        }
863
        NewSupps.append(BasisChange.getEquationsMatrix());
864
        Matrix<Integer> Ker=NewSupps.kernel(); // gives the new verterx
865
        // Ker.pretty_print(cout);
866
        assert(Ker.nr_of_rows()==1);
867
        Generators[Generators.nr_of_rows()-1]=Ker[0];
868
    }
869

870
    assert(Inequalities.nr_of_rows()==0 || Generators.nr_of_rows()==0);    
871

872
    if(Generators.nr_of_rows()==0)
873
        prepare_input_type_4(Inequalities); // inserts default inequalties if necessary
874
    else{
875
        is_Computed.set(ConeProperty::Generators);
876
        is_Computed.set(ConeProperty::Sublattice); 
877
    }
878
    
879
    checkGrading();
880
    checkDehomogenization();
881
    
882
    if(isComputed(ConeProperty::Grading)) {// cone known to be pointed
883
        is_Computed.set(ConeProperty::MaximalSubspace);
884
        pointed=true;
885
        is_Computed.set(ConeProperty::IsPointed);
886
    }
887

888
    setWeights();  // make matrix of weights for sorting
889

890
    if(PreComputedSupportHyperplanes.nr_of_rows()>0){
891
        check_precomputed_support_hyperplanes();
892
        SupportHyperplanes=PreComputedSupportHyperplanes;
893
        is_Computed.set(ConeProperty::SupportHyperplanes);
894
    }
895
    
896
    BasisChangePointed=BasisChange;
897
    
898
    is_Computed.set(ConeProperty::IsInhomogeneous);
899
    is_Computed.set(ConeProperty::EmbeddingDim);
900

901
    /* if(ExcludedFaces.nr_of_rows()>0){ // Nothing to check anymore
902
        check_excluded_faces();
903
    } */
904

905
    /*
906
    cout <<"-----------------------" << endl;
907
    cout << "Gen " << endl;
908
    Generators.pretty_print(cout);
909
    cout << "Supp " << endl;
910
    SupportHyperplanes.pretty_print(cout);
911
    cout << "A" << endl;
912
    BasisChange.get_A().pretty_print(cout);
913
    cout << "B" << endl;
914
    BasisChange.get_B().pretty_print(cout);
915
    cout <<"-----------------------" << endl;
916
    */
917
}
918

919

920
//---------------------------------------------------------------------------
921

922
template<typename Integer>
923
void Cone<Integer>::prepare_input_constraints(const map< InputType, vector< vector<Integer> > >& multi_input_data) {
924

925
    Matrix<Integer> Signs(0,dim), StrictSigns(0,dim);
926

927
    SupportHyperplanes=Matrix<Integer>(0,dim);
928
    Inequalities=Matrix<Integer>(0,dim);
929
    Equations=Matrix<Integer>(0,dim);
930
    Congruences=Matrix<Integer>(0,dim+1);
931

932
    typename map< InputType , vector< vector<Integer> > >::const_iterator it=multi_input_data.begin();
933

934
    it = multi_input_data.begin();
935
    for (; it != multi_input_data.end(); ++it) {
936

937
        switch (it->first) {
938
            case Type::strict_inequalities:
939
            case Type::inequalities:
940
            case Type::inhom_inequalities:
941
            case Type::excluded_faces:
942
                Inequalities.append(it->second);
943
                break;
944
            case Type::equations:
945
            case Type::inhom_equations:
946
                Equations.append(it->second);
947
                break;
948
            case Type::congruences:
949
            case Type::inhom_congruences:
950
                Congruences.append(it->second);
951
                break;
952
            case Type::signs:
953
                Signs = sign_inequalities(it->second);
954
                break;
955
            case Type::strict_signs:
956
                StrictSigns = strict_sign_inequalities(it->second);
957
                break;
958
            default:
959
                break;
960
        }
961
    }
962
    if(!BC_set) compose_basis_change(Sublattice_Representation<Integer>(dim));
963
    Matrix<Integer> Help(Signs);  // signs first !!
964
    Help.append(StrictSigns);   // then strict signs
965
    Help.append(Inequalities);
966
    Inequalities=Help;
967
}
968

969
//---------------------------------------------------------------------------
970
template<typename Integer>
971
void Cone<Integer>::prepare_input_generators(map< InputType, vector< vector<Integer> > >& multi_input_data, Matrix<Integer>& LatticeGenerators) {
972

973
    if(exists_element(multi_input_data,Type::vertices)){
974
        for(size_t i=0;i<multi_input_data[Type::vertices].size();++i)
975
            if(multi_input_data[Type::vertices][i][dim-1] <= 0) {
976
                throw BadInputException("Vertex has non-positive denominator!");
977
            }
978
    }
979

980
    if(exists_element(multi_input_data,Type::polyhedron)){
981
        for(size_t i=0;i<multi_input_data[Type::polyhedron].size();++i)
982
            if(multi_input_data[Type::polyhedron][i][dim-1] < 0) {
983
                throw BadInputException("Polyhedron vertex has negative denominator!");
984
            }
985
    }
986

987
    typename map< InputType , vector< vector<Integer> > >::const_iterator it=multi_input_data.begin();
988
    // find specific generator type -- there is only one, as checked already
989

990
    normalization=false;
991
    
992
    // check for subspace
993
    BasisMaxSubspace = find_input_matrix(multi_input_data,Type::subspace);
994
    if(BasisMaxSubspace.nr_of_rows()==0)
995
        BasisMaxSubspace=Matrix<Integer>(0,dim);
996
    
997
    vector<Integer> neg_sum_subspace(dim,0);
998
    for(size_t i=0;i<BasisMaxSubspace.nr_of_rows();++i)
999
        neg_sum_subspace=v_add(neg_sum_subspace,BasisMaxSubspace[i]);
1000
    v_scalar_multiplication<Integer>(neg_sum_subspace,-1);
1001
    
1002

1003
    Generators=Matrix<Integer>(0,dim);
1004
    for(; it != multi_input_data.end(); ++it) {
1005
        switch (it->first) {
1006
            case Type::normalization:
1007
            case Type::cone_and_lattice:
1008
                normalization=true;
1009
                LatticeGenerators.append(it->second);
1010
                if(BasisMaxSubspace.nr_of_rows()>0)
1011
                    LatticeGenerators.append(BasisMaxSubspace);
1012
            case Type::vertices:
1013
            case Type::polyhedron:
1014
            case Type::cone:
1015
            case Type::integral_closure:
1016
                Generators.append(it->second);
1017
                break;
1018
            case Type::subspace:
1019
                Generators.append(it->second);
1020
                Generators.append(neg_sum_subspace);
1021
                break;
1022
            case Type::polytope:
1023
                Generators.append(prepare_input_type_2(it->second));
1024
                break;
1025
            case Type::rees_algebra:
1026
                Generators.append(prepare_input_type_3(it->second));
1027
                break;
1028
            case Type::lattice:
1029
                LatticeGenerators.append(it->second);
1030
                break;
1031
            case Type::saturation:
1032
                LatticeGenerators.append(it->second);
1033
                LatticeGenerators.saturate();
1034
                break;
1035
            case Type::offset:
1036
                if(it->second.size()>1){
1037
                  throw BadInputException("Only one offset allowed!");
1038
                }
1039
                LatticeGenerators.append(it->second);
1040
                break;
1041
            default: break;
1042
        }
1043
    }
1044
}
1045

1046
//---------------------------------------------------------------------------
1047

1048
template<typename Integer>
1049
void Cone<Integer>::process_lattice_data(const Matrix<Integer>& LatticeGenerators, Matrix<Integer>& Congruences, Matrix<Integer>& Equations) {
1050

1051
    if(!BC_set)
1052
        compose_basis_change(Sublattice_Representation<Integer>(dim));
1053

1054
    bool no_constraints=(Congruences.nr_of_rows()==0) && (Equations.nr_of_rows()==0);
1055
    bool only_cone_gen=(Generators.nr_of_rows()!=0) && no_constraints && (LatticeGenerators.nr_of_rows()==0);
1056

1057
    no_lattice_restriction=true;
1058

1059
    if(only_cone_gen){
1060
        Sublattice_Representation<Integer> Basis_Change(Generators,true);
1061
        compose_basis_change(Basis_Change);
1062
        return;
1063
    }
1064

1065
    if(normalization && no_constraints){
1066
        Sublattice_Representation<Integer> Basis_Change(Generators,false);
1067
        compose_basis_change(Basis_Change);
1068
        return;
1069
    }
1070

1071
    no_lattice_restriction=false;
1072

1073
    if(Generators.nr_of_rows()!=0){
1074
        Equations.append(Generators.kernel());
1075
    }
1076

1077
    if(LatticeGenerators.nr_of_rows()!=0){
1078
        Sublattice_Representation<Integer> GenSublattice(LatticeGenerators,false);
1079
        if((Equations.nr_of_rows()==0) && (Congruences.nr_of_rows()==0)){
1080
            compose_basis_change(GenSublattice);
1081
            return;
1082
        }
1083
        Congruences.append(GenSublattice.getCongruencesMatrix());
1084
        Equations.append(GenSublattice.getEquationsMatrix());
1085
    }
1086

1087
    if (Congruences.nr_of_rows() > 0) {
1088
        bool zero_modulus;
1089
        Matrix<Integer> Ker_Basis=Congruences.solve_congruences(zero_modulus);
1090
        if(zero_modulus) {
1091
            throw BadInputException("Modulus 0 in congruence!");
1092
        }
1093
        Sublattice_Representation<Integer> Basis_Change(Ker_Basis,false);
1094
        compose_basis_change(Basis_Change);
1095
    }
1096

1097
    if (Equations.nr_of_rows()>0) {
1098
        Matrix<Integer> Ker_Basis=BasisChange.to_sublattice_dual(Equations).kernel();
1099
        Sublattice_Representation<Integer> Basis_Change(Ker_Basis,true);
1100
        compose_basis_change(Basis_Change);
1101
    }
1102
}
1103

1104
//---------------------------------------------------------------------------
1105

1106
template<typename Integer>
1107
void Cone<Integer>::prepare_input_type_4(Matrix<Integer>& Inequalities) {
1108

1109
    if (!inequalities_present) {
1110
        if (verbose) {
1111
            verboseOutput() << "No inequalities specified in constraint mode, using non-negative orthant." << endl;
1112
        }
1113
        if(inhomogeneous){
1114
            vector<Integer> test(dim);
1115
            test[dim-1]=1;
1116
            size_t matsize=dim;
1117
            if(test==Dehomogenization) // in this case "last coordinate >= 0" will come in through the dehomogenization
1118
                matsize=dim-1;   // we don't check for any other coincidence
1119
            Inequalities= Matrix<Integer>(matsize,dim);
1120
            for(size_t j=0;j<matsize;++j)
1121
                Inequalities[j][j]=1;
1122
        }
1123
        else
1124
            Inequalities = Matrix<Integer>(dim);
1125
    }
1126
    if(inhomogeneous)
1127
        SupportHyperplanes.append(Dehomogenization);
1128
    SupportHyperplanes.append(Inequalities);
1129
}
1130

1131

1132
//---------------------------------------------------------------------------
1133

1134
/* polytope input */
1135
template<typename Integer>
1136
Matrix<Integer> Cone<Integer>::prepare_input_type_2(const vector< vector<Integer> >& Input) {
1137
    size_t j;
1138
    size_t nr = Input.size();
1139
    //append a column of 1
1140
    Matrix<Integer> Generators(nr, dim);
1141
    for (size_t i=0; i<nr; i++) {
1142
        for (j=0; j<dim-1; j++)
1143
            Generators[i][j] = Input[i][j];
1144
        Generators[i][dim-1]=1;
1145
    }
1146
    // use the added last component as grading
1147
    Grading = vector<Integer>(dim,0);
1148
    Grading[dim-1] = 1;
1149
    is_Computed.set(ConeProperty::Grading);
1150
    GradingDenom=1;
1151
    is_Computed.set(ConeProperty::GradingDenom);
1152
    return Generators;
1153
}
1154

1155
//---------------------------------------------------------------------------
1156

1157
/* rees input */
1158
template<typename Integer>
1159
Matrix<Integer> Cone<Integer>::prepare_input_type_3(const vector< vector<Integer> >& InputV) {
1160
    Matrix<Integer> Input(InputV);
1161
    int i,j,k,nr_rows=Input.nr_of_rows(), nr_columns=Input.nr_of_columns();
1162
    // create cone generator matrix
1163
    Matrix<Integer> Full_Cone_Generators(nr_rows+nr_columns,nr_columns+1,0);
1164
    for (i = 0; i < nr_columns; i++) {
1165
        Full_Cone_Generators[i][i]=1;
1166
    }
1167
    for(i=0; i<nr_rows; i++){
1168
        Full_Cone_Generators[i+nr_columns][nr_columns]=1;
1169
        for(j=0; j<nr_columns; j++) {
1170
            Full_Cone_Generators[i+nr_columns][j]=Input[i][j];
1171
        }
1172
    }
1173
    // primarity test
1174
    vector<bool>  Prim_Test(nr_columns,false);
1175
    for (i=0; i<nr_rows; i++) {
1176
        k=0;
1177
        size_t v=0;
1178
        for(j=0; j<nr_columns; j++)
1179
            if (Input[i][j]!=0 ){
1180
                    k++;
1181
                    v=j;
1182
            }
1183
        if(k==1)
1184
            Prim_Test[v]=true;
1185
    }
1186
    rees_primary=true;
1187
    for(i=0; i<nr_columns; i++)
1188
        if(!Prim_Test[i])
1189
            rees_primary=false;
1190

1191
    is_Computed.set(ConeProperty::IsReesPrimary);
1192
    return Full_Cone_Generators;
1193
}
1194

1195

1196
//---------------------------------------------------------------------------
1197

1198
template<typename Integer>
1199
void Cone<Integer>::prepare_input_lattice_ideal(map< InputType, vector< vector<Integer> > >& multi_input_data) {
1200

1201
    Matrix<Integer> Binomials(find_input_matrix(multi_input_data,Type::lattice_ideal));
1202

1203
    if (Grading.size()>0) {
1204
        //check if binomials are homogeneous
1205
        vector<Integer> degrees = Binomials.MxV(Grading);
1206
        for (size_t i=0; i<degrees.size(); ++i) {
1207
            if (degrees[i]!=0) {
1208
                throw BadInputException("Grading gives non-zero value "
1209
                        + toString(degrees[i]) + " for binomial "
1210
                        + toString(i+1) + "!");
1211
            }
1212
            if (Grading[i] <0) {
1213
                throw BadInputException("Grading gives negative value "
1214
                        + toString(Grading[i]) + " for generator "
1215
                        + toString(i+1) + "!");
1216
            }
1217
        }
1218
    }
1219

1220
    Matrix<Integer> Gens=Binomials.kernel().transpose();
1221
    Full_Cone<Integer> FC(Gens);
1222
    FC.verbose=verbose;
1223
    if (verbose) verboseOutput() << "Computing a positive embedding..." << endl;
1224

1225
    FC.dualize_cone();
1226
    Matrix<Integer> Supp_Hyp=FC.getSupportHyperplanes().sort_lex();
1227
    Matrix<Integer> Selected_Supp_Hyp_Trans=(Supp_Hyp.submatrix(Supp_Hyp.max_rank_submatrix_lex())).transpose();
1228
    Matrix<Integer> Positive_Embedded_Generators=Gens.multiplication(Selected_Supp_Hyp_Trans);
1229
    // GeneratorsOfToricRing = Positive_Embedded_Generators;
1230
    // is_Computed.set(ConeProperty::GeneratorsOfToricRing);
1231
    dim = Positive_Embedded_Generators.nr_of_columns();
1232
    multi_input_data.insert(make_pair(Type::normalization,Positive_Embedded_Generators.get_elements())); // this is the cone defined by the binomials
1233

1234
    if (Grading.size()>0) {
1235
        // solve GeneratorsOfToricRing * grading = old_grading
1236
        Integer dummyDenom;
1237
        // Grading must be set directly since map entry has been processed already
1238
        Grading = Positive_Embedded_Generators.solve_rectangular(Grading,dummyDenom);
1239
        if (Grading.size() != dim) {
1240
            errorOutput() << "Grading could not be transferred!"<<endl;
1241
            is_Computed.set(ConeProperty::Grading, false);
1242
        }
1243
    }
1244
}
1245

1246
/* only used by the constructors */
1247
template<typename Integer>
1248
void Cone<Integer>::initialize() {
1249
    BC_set=false;
1250
    is_Computed = bitset<ConeProperty::EnumSize>();  //initialized to false
1251
    dim = 0;
1252
    unit_group_index = 1;
1253
    inhomogeneous=false;
1254
    rees_primary = false;
1255
    triangulation_is_nested = false;
1256
    triangulation_is_partial = false;
1257
    is_approximation=false;
1258
    verbose = libnormaliz::verbose; //take the global default
1259
    if (using_GMP<Integer>()) {
1260
        change_integer_type = true;
1261
    } else {
1262
        change_integer_type = false;
1263
    }
1264
    IntHullCone=NULL;
1265
    SymmCone=NULL;
1266
    
1267
    already_in_compute=false;
1268
    
1269
    set_parallelization();
1270
    nmz_interrupted=0;
1271
    nmz_scip=false;
1272
    
1273
}
1274

1275
template<typename Integer>
1276
void Cone<Integer>::set_parallelization() {
1277
    
1278
    omp_set_nested(0);
1279
    
1280
    if(thread_limit<0)
1281
        throw BadInputException("Invalid thread limit");
1282
    
1283
    if(parallelization_set){
1284
        if(thread_limit!=0)
1285
            omp_set_num_threads(thread_limit);        
1286
    }
1287
    else{
1288
        if(std::getenv("OMP_NUM_THREADS") == NULL){
1289
            long old=omp_get_max_threads();
1290
            if(old>default_thread_limit)
1291
                set_thread_limit(default_thread_limit);        
1292
            omp_set_num_threads(thread_limit);
1293
        }       
1294
    }
1295
}
1296

1297
//---------------------------------------------------------------------------
1298

1299
template<typename Integer>
1300
void Cone<Integer>::compose_basis_change(const Sublattice_Representation<Integer>& BC) {
1301
    if (BC_set) {
1302
        BasisChange.compose(BC);
1303
    } else {
1304
        BasisChange = BC;
1305
        BC_set = true;
1306
    }
1307
}
1308
//---------------------------------------------------------------------------
1309
template<typename Integer>
1310
void Cone<Integer>::check_precomputed_support_hyperplanes(){
1311

1312
    if (isComputed(ConeProperty::Generators)) {
1313
        // check if the inequalities are at least valid
1314
        // if (PreComputedSupportHyperplanes.nr_of_rows() != 0) {
1315
            Integer sp;
1316
            for (size_t i = 0; i < Generators.nr_of_rows(); ++i) {
1317
                for (size_t j = 0; j < PreComputedSupportHyperplanes.nr_of_rows(); ++j) {
1318
                    if ((sp = v_scalar_product(Generators[i], PreComputedSupportHyperplanes[j])) < 0) {
1319
                        throw BadInputException("Precomputed inequality " + toString(j)
1320
                                + " is not valid for generator " + toString(i)
1321
                                + " (value " + toString(sp) + ")");
1322
                    }
1323
                }
1324
            }
1325
        // }
1326
    }
1327
}
1328

1329
//---------------------------------------------------------------------------
1330
template<typename Integer>
1331
void Cone<Integer>::check_excluded_faces(){
1332

1333
    if (isComputed(ConeProperty::Generators)) {
1334
        // check if the inequalities are at least valid
1335
        // if (ExcludedFaces.nr_of_rows() != 0) {
1336
            Integer sp;
1337
            for (size_t i = 0; i < Generators.nr_of_rows(); ++i) {
1338
                for (size_t j = 0; j < ExcludedFaces.nr_of_rows(); ++j) {
1339
                    if ((sp = v_scalar_product(Generators[i], ExcludedFaces[j])) < 0) {
1340
                        throw BadInputException("Excluded face " + toString(j)
1341
                                + " is not valid for generator " + toString(i)
1342
                                + " (value " + toString(sp) + ")");
1343
                    }
1344
                }
1345
            }
1346
        // }
1347
    }
1348
}
1349

1350

1351
//---------------------------------------------------------------------------
1352

1353
template<typename Integer>
1354
bool Cone<Integer>::setVerbose (bool v) {
1355
    //we want to return the old value
1356
    bool old = verbose;
1357
    verbose = v;
1358
    return old;
1359
}
1360
//---------------------------------------------------------------------------
1361

1362
template<typename Integer>
1363
void Cone<Integer>::deactivateChangeOfPrecision() {
1364
    change_integer_type = false;
1365
}
1366

1367
//---------------------------------------------------------------------------
1368

1369
template<typename Integer>
1370
void Cone<Integer>::checkGrading () {
1371
    
1372
    if (isComputed(ConeProperty::Grading) || Grading.size()==0) {
1373
        return;
1374
    }
1375
    
1376
    bool positively_graded=true;
1377
    bool nonnegative=true;
1378
    size_t neg_index=0;
1379
    Integer neg_value;
1380
    if (Generators.nr_of_rows() > 0) {
1381
        vector<Integer> degrees = Generators.MxV(Grading);
1382
        for (size_t i=0; i<degrees.size(); ++i) {
1383
            if (degrees[i]<=0 && (!inhomogeneous || v_scalar_product(Generators[i],Dehomogenization)==0)) { 
1384
                // in the inhomogeneous case: test only generators of tail cone
1385
                positively_graded=false;;
1386
                if(degrees[i]<0){
1387
                    nonnegative=false;
1388
                    neg_index=i;
1389
                    neg_value=degrees[i];
1390
                }
1391
            }
1392
        }
1393
        if(positively_graded){
1394
            vector<Integer> test_grading=BasisChange.to_sublattice_dual_no_div(Grading);
1395
            GradingDenom=v_make_prime(test_grading);
1396
        }
1397
        else
1398
            GradingDenom = 1; 
1399
    } else {
1400
        GradingDenom = 1;
1401
    }
1402

1403
    if (isComputed(ConeProperty::Generators)){        
1404
        if(!nonnegative){
1405
            throw BadInputException("Grading gives negative value "
1406
                    + toString(neg_value) + " for generator "
1407
                    + toString(neg_index+1) + "!");
1408
        }
1409
        if(positively_graded){
1410
            is_Computed.set(ConeProperty::Grading);
1411
            is_Computed.set(ConeProperty::GradingDenom);            
1412
        }
1413
    }
1414
    
1415
}
1416

1417
//---------------------------------------------------------------------------
1418

1419
template<typename Integer>
1420
void Cone<Integer>::checkDehomogenization () {
1421
    if(Dehomogenization.size()>0){
1422
        vector<Integer> test=Generators.MxV(Dehomogenization);
1423
        for(size_t i=0;i<test.size();++i)
1424
            if(test[i]<0){
1425
                throw BadInputException(
1426
                        "Dehomogenization has has negative value on generator "
1427
                        + toString(Generators[i]));
1428
            }
1429
    }
1430
}
1431
//---------------------------------------------------------------------------
1432

1433
template<typename Integer>
1434
void Cone<Integer>::setGrading (const vector<Integer>& lf) {
1435
    
1436
    if (isComputed(ConeProperty::Grading) && Grading == lf) {
1437
        return;
1438
    }
1439
    
1440
    if (lf.size() != dim) {
1441
        throw BadInputException("Grading linear form has wrong dimension "
1442
                + toString(lf.size()) + " (should be " + toString(dim) + ")");
1443
    }
1444
    
1445
    Grading = lf;
1446
    checkGrading();
1447
}
1448

1449
//---------------------------------------------------------------------------
1450

1451
template<typename Integer>
1452
void Cone<Integer>::setWeights () {
1453

1454
    if(WeightsGrad.nr_of_columns()!=dim){
1455
        WeightsGrad=Matrix<Integer> (0,dim);  // weight matrix for ordering
1456
    }
1457
    if(Grading.size()>0 && WeightsGrad.nr_of_rows()==0)
1458
        WeightsGrad.append(Grading);
1459
    GradAbs=vector<bool>(WeightsGrad.nr_of_rows(),false);
1460
}
1461
//---------------------------------------------------------------------------
1462

1463
template<typename Integer>
1464
void Cone<Integer>::setDehomogenization (const vector<Integer>& lf) {
1465
    if (lf.size() != dim) {
1466
        throw BadInputException("Dehomogenizing linear form has wrong dimension "
1467
                + toString(lf.size()) + " (should be " + toString(dim) + ")");
1468
    }
1469
    Dehomogenization=lf;
1470
    is_Computed.set(ConeProperty::Dehomogenization);
1471
}
1472

1473
//---------------------------------------------------------------------------
1474

1475
/* check what is computed */
1476
template<typename Integer>
1477
bool Cone<Integer>::isComputed(ConeProperty::Enum prop) const {
1478
    return is_Computed.test(prop);
1479
}
1480

1481
template<typename Integer>
1482
bool Cone<Integer>::isComputed(ConeProperties CheckComputed) const {
1483
    return CheckComputed.reset(is_Computed).any();
1484
}
1485

1486
template<typename Integer>
1487
void Cone<Integer>::resetComputed(ConeProperty::Enum prop){
1488
    is_Computed.reset(prop);
1489
}
1490

1491

1492
/* getter */
1493

1494
template<typename Integer>
1495
Cone<Integer>& Cone<Integer>::getIntegerHullCone() const {
1496
    return *IntHullCone;
1497
}
1498

1499
template<typename Integer>
1500
Cone<Integer>& Cone<Integer>::getSymmetrizedCone() const {
1501
    return *SymmCone;
1502
}
1503

1504
template<typename Integer>
1505
size_t Cone<Integer>::getRank() {
1506
    compute(ConeProperty::Sublattice);
1507
    return BasisChange.getRank();
1508
}
1509

1510
template<typename Integer>    // computation depends on OriginalMonoidGenerators
1511
Integer Cone<Integer>::getIndex() {
1512
    compute(ConeProperty::OriginalMonoidGenerators);
1513
    return index;
1514
}
1515

1516
template<typename Integer>    // computation depends on OriginalMonoidGenerators
1517
Integer Cone<Integer>::getInternalIndex() {
1518
    return getIndex();
1519
}
1520

1521
template<typename Integer>
1522
Integer Cone<Integer>::getUnitGroupIndex() {
1523
    compute(ConeProperty::OriginalMonoidGenerators,ConeProperty::IsIntegrallyClosed);
1524
    return unit_group_index;
1525
}
1526

1527
template<typename Integer>
1528
size_t Cone<Integer>::getRecessionRank() {
1529
    compute(ConeProperty::RecessionRank);
1530
    return recession_rank;
1531
}
1532

1533
template<typename Integer>
1534
long Cone<Integer>::getAffineDim() {
1535
    compute(ConeProperty::AffineDim);
1536
    return affine_dim;
1537
}
1538

1539
template<typename Integer>
1540
const Sublattice_Representation<Integer>& Cone<Integer>::getSublattice() {
1541
    compute(ConeProperty::Sublattice);
1542
    return BasisChange;
1543
}
1544

1545
template<typename Integer>
1546
const Matrix<Integer>& Cone<Integer>::getOriginalMonoidGeneratorsMatrix() {
1547
    compute(ConeProperty::OriginalMonoidGenerators);
1548
    return OriginalMonoidGenerators;
1549
}
1550
template<typename Integer>
1551
const vector< vector<Integer> >& Cone<Integer>::getOriginalMonoidGenerators() {
1552
    compute(ConeProperty::OriginalMonoidGenerators);
1553
    return OriginalMonoidGenerators.get_elements();
1554
}
1555
template<typename Integer>
1556
size_t Cone<Integer>::getNrOriginalMonoidGenerators() {
1557
    compute(ConeProperty::OriginalMonoidGenerators);
1558
    return OriginalMonoidGenerators.nr_of_rows();
1559
}
1560

1561
template<typename Integer>
1562
const vector< vector<Integer> >& Cone<Integer>::getMaximalSubspace() {
1563
    compute(ConeProperty::MaximalSubspace);
1564
    return BasisMaxSubspace.get_elements();
1565
}
1566
template<typename Integer>
1567
const Matrix<Integer>& Cone<Integer>::getMaximalSubspaceMatrix() {
1568
    compute(ConeProperty::MaximalSubspace);
1569
    return BasisMaxSubspace;
1570
}
1571
template<typename Integer>
1572
size_t Cone<Integer>::getDimMaximalSubspace() {
1573
    compute(ConeProperty::MaximalSubspace);
1574
    return BasisMaxSubspace.nr_of_rows();
1575
}
1576

1577
template<typename Integer>
1578
const Matrix<Integer>& Cone<Integer>::getGeneratorsMatrix() {
1579
    compute(ConeProperty::Generators);
1580
    return Generators;
1581
}
1582

1583
template<typename Integer>
1584
const vector< vector<Integer> >& Cone<Integer>::getGenerators() {
1585
    compute(ConeProperty::Generators);
1586
    return Generators.get_elements();
1587
}
1588

1589
template<typename Integer>
1590
size_t Cone<Integer>::getNrGenerators() {
1591
    compute(ConeProperty::Generators);
1592
    return Generators.nr_of_rows();
1593
}
1594

1595
template<typename Integer>
1596
const Matrix<Integer>& Cone<Integer>::getExtremeRaysMatrix() {
1597
    compute(ConeProperty::ExtremeRays);
1598
    return ExtremeRays;
1599
}
1600
template<typename Integer>
1601
const vector< vector<Integer> >& Cone<Integer>::getExtremeRays() {
1602
    compute(ConeProperty::ExtremeRays);
1603
    return ExtremeRays.get_elements();
1604
}
1605
template<typename Integer>
1606
size_t Cone<Integer>::getNrExtremeRays() {
1607
    compute(ConeProperty::ExtremeRays);
1608
    return ExtremeRays.nr_of_rows();
1609
}
1610

1611
template<typename Integer>
1612
const Matrix<nmz_float>& Cone<Integer>::getVerticesFloatMatrix() {
1613
    compute(ConeProperty::VerticesFloat);
1614
    return VerticesFloat;
1615
}
1616
template<typename Integer>
1617
const vector< vector<nmz_float> >& Cone<Integer>::getVerticesFloat() {
1618
    compute(ConeProperty::VerticesFloat);
1619
    return VerticesFloat.get_elements();
1620
}
1621
template<typename Integer>
1622
size_t Cone<Integer>::getNrVerticesFloat() {
1623
    compute(ConeProperty::VerticesFloat);
1624
    return VerticesFloat.nr_of_rows();
1625
}
1626

1627
template<typename Integer>
1628
const Matrix<Integer>& Cone<Integer>::getVerticesOfPolyhedronMatrix() {
1629
    compute(ConeProperty::VerticesOfPolyhedron);
1630
    return VerticesOfPolyhedron;
1631
}
1632
template<typename Integer>
1633
const vector< vector<Integer> >& Cone<Integer>::getVerticesOfPolyhedron() {
1634
    compute(ConeProperty::VerticesOfPolyhedron);
1635
    return VerticesOfPolyhedron.get_elements();
1636
}
1637
template<typename Integer>
1638
size_t Cone<Integer>::getNrVerticesOfPolyhedron() {
1639
    compute(ConeProperty::VerticesOfPolyhedron);
1640
    return VerticesOfPolyhedron.nr_of_rows();
1641
}
1642

1643
template<typename Integer>
1644
const Matrix<Integer>& Cone<Integer>::getSupportHyperplanesMatrix() {
1645
    compute(ConeProperty::SupportHyperplanes);
1646
    return SupportHyperplanes;
1647
}
1648
template<typename Integer>
1649
const vector< vector<Integer> >& Cone<Integer>::getSupportHyperplanes() {
1650
    compute(ConeProperty::SupportHyperplanes);
1651
    return SupportHyperplanes.get_elements();
1652
}
1653
template<typename Integer>
1654
size_t Cone<Integer>::getNrSupportHyperplanes() {
1655
    compute(ConeProperty::SupportHyperplanes);
1656
    return SupportHyperplanes.nr_of_rows();
1657
}
1658

1659
template<typename Integer>
1660
map< InputType , vector< vector<Integer> > > Cone<Integer>::getConstraints () {
1661
    compute(ConeProperty::Sublattice, ConeProperty::SupportHyperplanes);
1662
    map<InputType, vector< vector<Integer> > > c;
1663
    c[Type::inequalities] = SupportHyperplanes.get_elements();
1664
    c[Type::equations] = BasisChange.getEquations();
1665
    c[Type::congruences] = BasisChange.getCongruences();
1666
    return c;
1667
}
1668

1669
template<typename Integer>
1670
const Matrix<Integer>& Cone<Integer>::getExcludedFacesMatrix() {
1671
    compute(ConeProperty::ExcludedFaces);
1672
    return ExcludedFaces;
1673
}
1674
template<typename Integer>
1675
const vector< vector<Integer> >& Cone<Integer>::getExcludedFaces() {
1676
    compute(ConeProperty::ExcludedFaces);
1677
    return ExcludedFaces.get_elements();
1678
}
1679
template<typename Integer>
1680
size_t Cone<Integer>::getNrExcludedFaces() {
1681
    compute(ConeProperty::ExcludedFaces);
1682
    return ExcludedFaces.nr_of_rows();
1683
}
1684

1685
template<typename Integer>
1686
const vector< pair<vector<key_t>,Integer> >& Cone<Integer>::getTriangulation() {
1687
    compute(ConeProperty::Triangulation);
1688
    return Triangulation;
1689
}
1690

1691
template<typename Integer>
1692
const vector<vector<bool> >& Cone<Integer>::getOpenFacets() {
1693
    compute(ConeProperty::ConeDecomposition);
1694
    return OpenFacets;
1695
}
1696

1697
template<typename Integer>
1698
const vector< pair<vector<key_t>,long> >& Cone<Integer>::getInclusionExclusionData() {
1699
    compute(ConeProperty::InclusionExclusionData);
1700
    return InExData;
1701
}
1702

1703
template<typename Integer>
1704
void Cone<Integer>::make_StanleyDec_export() {
1705
    if(!StanleyDec_export.empty())
1706
        return;
1707
    assert(isComputed(ConeProperty::StanleyDec));
1708
    auto SD=StanleyDec.begin();
1709
    for(;SD!=StanleyDec.end();++SD){
1710
        STANLEYDATA<Integer> NewSt;
1711
        NewSt.key=SD->key;
1712
        convert(NewSt.offsets,SD->offsets);
1713
        StanleyDec_export.push_back(NewSt);        
1714
    }    
1715
}
1716

1717
template<typename Integer>
1718
const list< STANLEYDATA<Integer> >& Cone<Integer>::getStanleyDec() {
1719
    compute(ConeProperty::StanleyDec);
1720
    make_StanleyDec_export();
1721
    return StanleyDec_export;
1722
}
1723

1724
template<typename Integer>
1725
list< STANLEYDATA_int >& Cone<Integer>::getStanleyDec_mutable() {
1726
    assert(isComputed(ConeProperty::StanleyDec));
1727
    return StanleyDec;
1728
}
1729

1730
template<typename Integer>
1731
size_t Cone<Integer>::getTriangulationSize() {
1732
    compute(ConeProperty::TriangulationSize);
1733
    return TriangulationSize;
1734
}
1735

1736
template<typename Integer>
1737
Integer Cone<Integer>::getTriangulationDetSum() {
1738
    compute(ConeProperty::TriangulationDetSum);
1739
    return TriangulationDetSum;
1740
}
1741

1742
template<typename Integer>
1743
vector<Integer> Cone<Integer>::getWitnessNotIntegrallyClosed() {
1744
    compute(ConeProperty::WitnessNotIntegrallyClosed);
1745
    return WitnessNotIntegrallyClosed;
1746
}
1747

1748
template<typename Integer>
1749
vector<Integer> Cone<Integer>::getGeneratorOfInterior() {
1750
    compute(ConeProperty::GeneratorOfInterior);
1751
    return GeneratorOfInterior;
1752
}
1753

1754
template<typename Integer>
1755
const Matrix<Integer>& Cone<Integer>::getHilbertBasisMatrix() {
1756
    compute(ConeProperty::HilbertBasis);
1757
    return HilbertBasis;
1758
}
1759
template<typename Integer>
1760
const vector< vector<Integer> >& Cone<Integer>::getHilbertBasis() {
1761
    compute(ConeProperty::HilbertBasis);
1762
    return HilbertBasis.get_elements();
1763
}
1764
template<typename Integer>
1765
size_t Cone<Integer>::getNrHilbertBasis() {
1766
    compute(ConeProperty::HilbertBasis);
1767
    return HilbertBasis.nr_of_rows();
1768
}
1769

1770
template<typename Integer>
1771
const Matrix<Integer>& Cone<Integer>::getModuleGeneratorsOverOriginalMonoidMatrix() {
1772
    compute(ConeProperty::ModuleGeneratorsOverOriginalMonoid);
1773
    return ModuleGeneratorsOverOriginalMonoid;
1774
}
1775
template<typename Integer>
1776
const vector< vector<Integer> >& Cone<Integer>::getModuleGeneratorsOverOriginalMonoid() {
1777
    compute(ConeProperty::ModuleGeneratorsOverOriginalMonoid);
1778
    return ModuleGeneratorsOverOriginalMonoid.get_elements();
1779
}
1780
template<typename Integer>
1781
size_t Cone<Integer>::getNrModuleGeneratorsOverOriginalMonoid() {
1782
    compute(ConeProperty::ModuleGeneratorsOverOriginalMonoid);
1783
    return ModuleGeneratorsOverOriginalMonoid.nr_of_rows();
1784
}
1785

1786
template<typename Integer>
1787
const Matrix<Integer>& Cone<Integer>::getModuleGeneratorsMatrix() {
1788
    compute(ConeProperty::ModuleGenerators);
1789
    return ModuleGenerators;
1790
}
1791
template<typename Integer>
1792
const vector< vector<Integer> >& Cone<Integer>::getModuleGenerators() {
1793
    compute(ConeProperty::ModuleGenerators);
1794
    return ModuleGenerators.get_elements();
1795
}
1796
template<typename Integer>
1797
size_t Cone<Integer>::getNrModuleGenerators() {
1798
    compute(ConeProperty::ModuleGenerators);
1799
    return ModuleGenerators.nr_of_rows();
1800
}
1801

1802
template<typename Integer>
1803
const Matrix<Integer>& Cone<Integer>::getDeg1ElementsMatrix() {
1804
    compute(ConeProperty::Deg1Elements);
1805
    return Deg1Elements;
1806
}
1807
template<typename Integer>
1808
const vector< vector<Integer> >& Cone<Integer>::getDeg1Elements() {
1809
    compute(ConeProperty::Deg1Elements);
1810
    return Deg1Elements.get_elements();
1811
}
1812
template<typename Integer>
1813
size_t Cone<Integer>::getNrDeg1Elements() {
1814
    compute(ConeProperty::Deg1Elements);
1815
    return Deg1Elements.nr_of_rows();
1816
}
1817

1818
template<typename Integer>
1819
const HilbertSeries& Cone<Integer>::getHilbertSeries() {
1820
    compute(ConeProperty::HilbertSeries);
1821
    return HSeries;
1822
}
1823

1824
template<typename Integer>
1825
vector<Integer> Cone<Integer>::getGrading() {
1826
    compute(ConeProperty::Grading);
1827
    return Grading;
1828
}
1829

1830
template<typename Integer>
1831
Integer Cone<Integer>::getGradingDenom() {
1832
    compute(ConeProperty::Grading);
1833
    return GradingDenom;
1834
}
1835

1836
template<typename Integer>
1837
vector<Integer> Cone<Integer>::getDehomogenization() {
1838
    compute(ConeProperty::Dehomogenization);
1839
    return Dehomogenization;
1840
}
1841

1842
template<typename Integer>
1843
mpq_class Cone<Integer>::getMultiplicity() {
1844
    compute(ConeProperty::Multiplicity);
1845
    return multiplicity;
1846
}
1847

1848
template<typename Integer>
1849
mpq_class Cone<Integer>::getVirtualMultiplicity() {
1850
    if(!isComputed(ConeProperty::VirtualMultiplicity)) // in order not to compute the triangulation
1851
        compute(ConeProperty::VirtualMultiplicity);    // which is deleted if not asked for explicitly
1852
    return IntData.getVirtualMultiplicity();
1853
}
1854

1855
template<typename Integer>
1856
const pair<HilbertSeries, mpz_class>& Cone<Integer>::getWeightedEhrhartSeries(){
1857
    if(!isComputed(ConeProperty::WeightedEhrhartSeries))  // see above
1858
        compute(ConeProperty::WeightedEhrhartSeries);
1859
    return getIntData().getWeightedEhrhartSeries();
1860
}
1861

1862
template<typename Integer>
1863
IntegrationData& Cone<Integer>::getIntData(){
1864
    return IntData;
1865
}
1866

1867
template<typename Integer>
1868
mpq_class Cone<Integer>::getIntegral() {
1869
    if(!isComputed(ConeProperty::Integral)) // see above
1870
        compute(ConeProperty::Integral);
1871
    return IntData.getIntegral();
1872
}
1873

1874
template<typename Integer>
1875
bool Cone<Integer>::isPointed() {
1876
    compute(ConeProperty::IsPointed);
1877
    return pointed;
1878
}
1879

1880
template<typename Integer>
1881
bool Cone<Integer>::isInhomogeneous() {
1882
    return inhomogeneous;
1883
}
1884

1885
template<typename Integer>
1886
bool Cone<Integer>::isDeg1ExtremeRays() {
1887
    compute(ConeProperty::IsDeg1ExtremeRays);
1888
    return deg1_extreme_rays;
1889
}
1890

1891
template<typename Integer>
1892
bool Cone<Integer>::isGorenstein() {
1893
    compute(ConeProperty::IsGorenstein);
1894
    return Gorenstein;
1895
}
1896

1897
template<typename Integer>
1898
bool Cone<Integer>::isDeg1HilbertBasis() {
1899
    compute(ConeProperty::IsDeg1HilbertBasis);
1900
    return deg1_hilbert_basis;
1901
}
1902

1903
template<typename Integer>
1904
bool Cone<Integer>::isIntegrallyClosed() {
1905
    compute(ConeProperty::IsIntegrallyClosed);
1906
    return integrally_closed;
1907
}
1908

1909
template<typename Integer>
1910
bool Cone<Integer>::isReesPrimary() {
1911
    compute(ConeProperty::IsReesPrimary);
1912
    return rees_primary;
1913
}
1914

1915
template<typename Integer>
1916
Integer Cone<Integer>::getReesPrimaryMultiplicity() {
1917
    compute(ConeProperty::ReesPrimaryMultiplicity);
1918
    return ReesPrimaryMultiplicity;
1919
}
1920

1921
// the information about the triangulation will just be returned
1922
// if no triangulation was computed so far they return false
1923
template<typename Integer>
1924
bool Cone<Integer>::isTriangulationNested() {
1925
    return triangulation_is_nested;
1926
}
1927
template<typename Integer>
1928
bool Cone<Integer>::isTriangulationPartial() {
1929
    return triangulation_is_partial;
1930
}
1931

1932
template<typename Integer>
1933
size_t Cone<Integer>::getModuleRank() {
1934
    compute(ConeProperty::ModuleRank);
1935
    return module_rank;
1936
}
1937

1938
template<typename Integer>
1939
vector<Integer> Cone<Integer>::getClassGroup() {
1940
    compute(ConeProperty::ClassGroup);
1941
    return ClassGroup;
1942
}
1943

1944
//---------------------------------------------------------------------------
1945

1946
template<typename Integer>
1947
ConeProperties Cone<Integer>::compute(ConeProperty::Enum cp) {
1948
    if (isComputed(cp)) return ConeProperties();
1949
    return compute(ConeProperties(cp));
1950
}
1951

1952
template<typename Integer>
1953
ConeProperties Cone<Integer>::compute(ConeProperty::Enum cp1, ConeProperty::Enum cp2) {
1954
    return compute(ConeProperties(cp1,cp2));
1955
}
1956

1957
template<typename Integer>
1958
ConeProperties Cone<Integer>::compute(ConeProperty::Enum cp1, ConeProperty::Enum cp2,
1959
                                      ConeProperty::Enum cp3) {
1960
    return compute(ConeProperties(cp1,cp2,cp3));
1961
}
1962

1963
//---------------------------------------------------------------------------
1964

1965
template<typename Integer>
1966
void Cone<Integer>::set_implicit_dual_mode(ConeProperties& ToCompute) {
1967
    
1968
    if(ToCompute.test(ConeProperty::DualMode) || ToCompute.test(ConeProperty::PrimalMode)
1969
                    || ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid)
1970
                    || ToCompute.test(ConeProperty::Approximate)
1971
                    || ToCompute.test(ConeProperty::Projection)
1972
                    || nr_cone_gen>0 || nr_latt_gen>0 || SupportHyperplanes.nr_of_rows() > 2*dim
1973
                    || SupportHyperplanes.nr_of_rows() 
1974
                            <= BasisChangePointed.getRank()+ 50/(BasisChangePointed.getRank()+1))
1975
        return;
1976
    if(ToCompute.test(ConeProperty::HilbertBasis))
1977
        ToCompute.set(ConeProperty::DualMode);
1978
    if(ToCompute.test(ConeProperty::Deg1Elements) 
1979
            && !(ToCompute.test(ConeProperty::HilbertSeries) || ToCompute.test(ConeProperty::Multiplicity)))
1980
        ToCompute.set(ConeProperty::DualMode);
1981
    return;
1982
}
1983

1984
//---------------------------------------------------------------------------
1985

1986
// this wrapper allows us to save and restore class data that depend on ToCompute
1987
// and may therefore be destroyed if compute() is called by itself
1988
template<typename Integer>
1989
ConeProperties Cone<Integer>::recursive_compute(ConeProperties ToCompute) {
1990
    
1991
    bool save_explicit_HilbertSeries=explicit_HilbertSeries;
1992
    bool save_naked_dual= naked_dual;
1993
    bool save_default_mode= default_mode;
1994
    already_in_compute=false;
1995
    ToCompute=compute_inner(ToCompute);
1996
    explicit_HilbertSeries=save_explicit_HilbertSeries;
1997
    naked_dual=save_naked_dual;
1998
    default_mode= save_default_mode;
1999
    return ToCompute;
2000
}
2001

2002
//---------------------------------------------------------------------------
2003

2004
template<typename Integer>
2005
ConeProperties Cone<Integer>::compute(ConeProperties ToCompute) {
2006
    already_in_compute=false;
2007
    set_parallelization();
2008
    nmz_interrupted=0;
2009
    if(ToCompute.test(ConeProperty::SCIP)){
2010
#ifdef NMZ_SCIP
2011
        nmz_scip=true;
2012
        ToCompute.set(ConeProperty::SCIP);
2013
#else
2014
        throw BadInputException("Option SCIP only allowed if Normaliz was built with Scip");
2015
#endif // NMZ_SCIP
2016
    }
2017
    return compute_inner(ToCompute);
2018
}
2019

2020
//---------------------------------------------------------------------------
2021

2022
template<typename Integer>
2023
ConeProperties Cone<Integer>::compute_inner(ConeProperties ToCompute) {
2024
    
2025
    assert(!already_in_compute);
2026
    already_in_compute=true;
2027
    
2028
    if(ToCompute.test(ConeProperty::NoPeriodBound)){
2029
        HSeries.set_period_bounded(false);
2030
        IntData.getWeightedEhrhartSeries().first.set_period_bounded(false);        
2031
    }
2032
        
2033
    
2034
#ifndef NMZ_COCOA
2035
   if(ToCompute.test(ConeProperty::VirtualMultiplicity) || ToCompute.test(ConeProperty::Integral) 
2036
       || ToCompute.test(ConeProperty::WeightedEhrhartSeries))
2037
       throw BadInputException("Integral, VirtualMultiplicity, WeightedEhrhartSeries only computable with CoCoALib");
2038
#endif
2039

2040
    default_mode=ToCompute.test(ConeProperty::DefaultMode);
2041
    
2042
    if(ToCompute.test(ConeProperty::BigInt)){
2043
        if(!using_GMP<Integer>())
2044
            throw BadInputException("BigInt can only be set for cones of Integer type GMP");
2045
        change_integer_type=false;
2046
    }
2047
    
2048
    if(ToCompute.test(ConeProperty::KeepOrder) && !isComputed(ConeProperty::OriginalMonoidGenerators))
2049
        throw BadInputException("KeepOrder can only be set if OriginalMonoidGenerators are defined");
2050
    
2051
    INTERRUPT_COMPUTATION_BY_EXCEPTION
2052
    
2053
    if(BasisMaxSubspace.nr_of_rows()>0 && !isComputed(ConeProperty::MaximalSubspace)){
2054
        BasisMaxSubspace=Matrix<Integer>(0,dim);
2055
        recursive_compute(ConeProperty::MaximalSubspace);      
2056
    }
2057
    
2058
    explicit_HilbertSeries=ToCompute.test(ConeProperty::HilbertSeries) || ToCompute.test(ConeProperty::HSOP);
2059
    // must distiguish it from being set through DefaultMode;
2060
    naked_dual=ToCompute.test(ConeProperty::DualMode) 
2061
                && !(ToCompute.test(ConeProperty::HilbertBasis) || ToCompute.test(ConeProperty::Deg1Elements));
2062
    // to control the computation of rational solutions in the inhomogeneous case
2063
    
2064
    ToCompute.reset(is_Computed);
2065
    ToCompute.check_conflicting_variants();
2066
    ToCompute.set_preconditions(inhomogeneous);
2067
    ToCompute.prepare_compute_options(inhomogeneous);
2068
    ToCompute.check_sanity(inhomogeneous);
2069
    if (!isComputed(ConeProperty::OriginalMonoidGenerators)) {
2070
        if (ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) {
2071
            errorOutput() << "ERROR: Module generators over original monoid only computable if original monoid is defined!"
2072
                << endl;
2073
            throw NotComputableException(ConeProperty::ModuleGeneratorsOverOriginalMonoid);
2074
        }
2075
        if (ToCompute.test(ConeProperty::IsIntegrallyClosed)) {
2076
            errorOutput() << "ERROR: Original monoid is not defined, cannot check it for being integrally closed."
2077
                << endl;
2078
            throw NotComputableException(ConeProperty::IsIntegrallyClosed);
2079
        }
2080
    }
2081
    
2082
    try_symmetrization(ToCompute);   
2083
    ToCompute.reset(is_Computed);
2084
    if (ToCompute.none()) {
2085
        already_in_compute=false; return ToCompute;
2086
    }
2087
    
2088
    INTERRUPT_COMPUTATION_BY_EXCEPTION
2089

2090
    
2091
    set_implicit_dual_mode(ToCompute);
2092

2093
    if (ToCompute.test(ConeProperty::DualMode)) {
2094
        compute_dual(ToCompute);
2095
    }
2096

2097
    if (ToCompute.test(ConeProperty::WitnessNotIntegrallyClosed)) {
2098
        find_witness();
2099
    }
2100

2101
    ToCompute.reset(is_Computed);
2102
    if (ToCompute.none()) {
2103
        already_in_compute=false; return ToCompute;
2104
    }
2105
    
2106
    INTERRUPT_COMPUTATION_BY_EXCEPTION
2107
        
2108
    try_approximation_or_projection(ToCompute);
2109
    
2110
    ToCompute.reset(is_Computed); // already computed
2111
    if (ToCompute.none()) {
2112
        already_in_compute=false; return ToCompute;
2113
    }
2114

2115
    /* preparation: get generators if necessary */
2116
    compute_generators();
2117

2118
    if (!isComputed(ConeProperty::Generators)) {
2119
        throw FatalException("Could not get Generators.");
2120
    }
2121

2122
    if (rees_primary && (ToCompute.test(ConeProperty::ReesPrimaryMultiplicity)
2123
            || ToCompute.test(ConeProperty::Multiplicity)
2124
            || ToCompute.test(ConeProperty::HilbertSeries)
2125
            || ToCompute.test(ConeProperty::DefaultMode) ) ) {
2126
        ReesPrimaryMultiplicity = compute_primary_multiplicity();
2127
        is_Computed.set(ConeProperty::ReesPrimaryMultiplicity);
2128
    }
2129

2130
    ToCompute.reset(is_Computed); // already computed
2131
    if (ToCompute.none()) {
2132
        already_in_compute=false; return ToCompute;
2133
    }
2134

2135
    // the computation of the full cone
2136
    if (change_integer_type) {
2137
        try {
2138
            compute_full_cone<MachineInteger>(ToCompute);
2139
        } catch(const ArithmeticException& e) {
2140
            if (verbose) {
2141
                verboseOutput() << e.what() << endl;
2142
                verboseOutput() << "Restarting with a bigger type." << endl;
2143
            }
2144
            change_integer_type = false;
2145
        }
2146
    }
2147
    
2148
    if (!change_integer_type) {
2149
        compute_full_cone<Integer>(ToCompute);
2150
    }
2151
    
2152
    check_Gorenstein(ToCompute);
2153
    
2154
    if(ToCompute.test(ConeProperty::IntegerHull)) {
2155
        compute_integer_hull();
2156
    }
2157
    
2158
    INTERRUPT_COMPUTATION_BY_EXCEPTION
2159
    
2160
    complete_HilbertSeries_comp(ToCompute);
2161
    
2162
    complete_sublattice_comp(ToCompute);
2163
    
2164
    compute_vertices_float(ToCompute);
2165
       
2166
    if(ToCompute.test(ConeProperty::WeightedEhrhartSeries))
2167
        compute_weighted_Ehrhart(ToCompute);
2168
    ToCompute.reset(is_Computed);
2169
    
2170
    if(ToCompute.test(ConeProperty::Integral))
2171
        compute_integral(ToCompute);
2172
    ToCompute.reset(is_Computed);
2173
    
2174
    if(ToCompute.test(ConeProperty::VirtualMultiplicity))
2175
        compute_virt_mult(ToCompute);
2176
    ToCompute.reset(is_Computed);
2177

2178
    /* check if everything is computed */
2179
    ToCompute.reset(is_Computed); //remove what is now computed
2180
    if (ToCompute.test(ConeProperty::Deg1Elements) && isComputed(ConeProperty::Grading)) {
2181
        // this can happen when we were looking for a witness earlier
2182
        recursive_compute(ToCompute);
2183
    }
2184
    if (!ToCompute.test(ConeProperty::DefaultMode) && ToCompute.goals().any()) {
2185
        throw NotComputableException(ToCompute.goals());
2186
    }
2187
    ToCompute.reset_compute_options();
2188
    already_in_compute=false; return ToCompute;
2189
}
2190

2191
template<typename Integer>
2192
void Cone<Integer>::compute_integer_hull() {
2193
    
2194
    if(verbose){
2195
        verboseOutput() << "Computing integer hull" << endl;
2196
    }
2197

2198
    Matrix<Integer> IntHullGen;
2199
    bool IntHullComputable=true;
2200
    size_t nr_extr=0;
2201
    if(inhomogeneous){
2202
        if(!isComputed(ConeProperty::HilbertBasis))
2203
            IntHullComputable=false;
2204
        IntHullGen=HilbertBasis;
2205
        IntHullGen.append(ModuleGenerators);
2206
    }
2207
    else{
2208
        if(!isComputed(ConeProperty::Deg1Elements))
2209
            IntHullComputable=false;
2210
        IntHullGen=Deg1Elements;
2211
    }
2212
    ConeProperties IntHullCompute;
2213
    IntHullCompute.set(ConeProperty::SupportHyperplanes);
2214
    if(!IntHullComputable){
2215
        errorOutput() << "Integer hull not computable: no integer points available." << endl;
2216
        throw NotComputableException(IntHullCompute);
2217
    }
2218
    
2219
    if(IntHullGen.nr_of_rows()==0){
2220
        IntHullGen.append(vector<Integer>(dim,0)); // we need a non-empty input matrix
2221
    }
2222
    
2223
    INTERRUPT_COMPUTATION_BY_EXCEPTION
2224
    
2225
    if(!inhomogeneous || HilbertBasis.nr_of_rows()==0){
2226
        nr_extr=IntHullGen.extreme_points_first();
2227
        if(verbose){
2228
            verboseOutput() << nr_extr << " extreme points found"  << endl;
2229
        }
2230
    }
2231
    else{ // now an unbounded polyhedron
2232
        if(isComputed(ConeProperty::Grading)){
2233
            nr_extr=IntHullGen.extreme_points_first(Grading);
2234
        }
2235
        else{
2236
            if(isComputed(ConeProperty::SupportHyperplanes)){
2237
                vector<Integer> aux_grading=SupportHyperplanes.find_inner_point();
2238
                nr_extr=IntHullGen.extreme_points_first(aux_grading);
2239
            }    
2240
        }
2241
    }
2242
    
2243
    // IntHullGen.pretty_print(cout);
2244
    IntHullCone=new Cone<Integer>(InputType::cone_and_lattice,IntHullGen.get_elements(), Type::subspace,BasisMaxSubspace);
2245
    if(nr_extr!=0)  // we suppress the ordering in full_cone only if we have found few extreme rays
2246
        IntHullCompute.set(ConeProperty::KeepOrder);
2247

2248
    IntHullCone->inhomogeneous=true; // inhomogeneous;
2249
    if(inhomogeneous)
2250
        IntHullCone->Dehomogenization=Dehomogenization;
2251
    else
2252
        IntHullCone->Dehomogenization=Grading;
2253
    IntHullCone->verbose=verbose;
2254
    try{
2255
        IntHullCone->compute(IntHullCompute);
2256
        if(IntHullCone->isComputed(ConeProperty::SupportHyperplanes))
2257
            is_Computed.set(ConeProperty::IntegerHull);
2258
        if(verbose){
2259
            verboseOutput() << "Integer hull finished" << endl;
2260
        }
2261
    }
2262
    catch (const NotComputableException& ){
2263
            errorOutput() << "Error in computation of integer hull" << endl;
2264
    }
2265
}
2266

2267
template<typename Integer>
2268
void Cone<Integer>::check_vanishing_of_grading_and_dehom(){
2269
    if(Grading.size()>0){
2270
        vector<Integer> test=BasisMaxSubspace.MxV(Grading);
2271
        if(test!=vector<Integer>(test.size())){
2272
                throw BadInputException("Grading does not vanish on maximal subspace.");
2273
        }
2274
    }
2275
    if(Dehomogenization.size()>0){
2276
        vector<Integer> test=BasisMaxSubspace.MxV(Dehomogenization);
2277
        if(test!=vector<Integer>(test.size())){
2278
            throw BadInputException("Dehomogenization does not vanish on maximal subspace.");
2279
        }
2280
    }    
2281
}
2282

2283
template<typename Integer>
2284
template<typename IntegerFC>
2285
void Cone<Integer>::compute_full_cone(ConeProperties& ToCompute) {
2286
    
2287
    if(ToCompute.test(ConeProperty::IsPointed) && Grading.size()==0){
2288
        if (verbose) {
2289
            verboseOutput()<<  "Checking pointedness first"<< endl;
2290
        }
2291
        ConeProperties Dualize;
2292
        Dualize.set(ConeProperty::SupportHyperplanes);
2293
        Dualize.set(ConeProperty::ExtremeRays);
2294
        recursive_compute(Dualize);
2295
    }
2296
    
2297
    Matrix<IntegerFC> FC_Gens;
2298

2299
    BasisChangePointed.convert_to_sublattice(FC_Gens, Generators);
2300
    Full_Cone<IntegerFC> FC(FC_Gens,!ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid));
2301
    // !ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid) blocks make_prime in full_cone.cpp
2302

2303
    /* activate bools in FC */
2304

2305
    FC.verbose=verbose;
2306

2307
    FC.inhomogeneous=inhomogeneous;
2308
    FC.explicit_h_vector=explicit_HilbertSeries;
2309

2310
    if (ToCompute.test(ConeProperty::HilbertSeries)) {
2311
        FC.do_h_vector = true;
2312
    }
2313
    if (ToCompute.test(ConeProperty::HilbertBasis)) {
2314
        FC.do_Hilbert_basis = true;
2315
    }
2316
    if (ToCompute.test(ConeProperty::IsIntegrallyClosed)) {
2317
        FC.do_integrally_closed = true;
2318
    }
2319
    if (ToCompute.test(ConeProperty::Triangulation)) {
2320
        FC.keep_triangulation = true;
2321
    }
2322
    if (ToCompute.test(ConeProperty::ConeDecomposition)) {
2323
        FC.do_cone_dec = true;
2324
    }
2325
    if (ToCompute.test(ConeProperty::Multiplicity) ) {
2326
        FC.do_multiplicity = true;
2327
    }
2328
    if (ToCompute.test(ConeProperty::TriangulationDetSum) ) {
2329
        FC.do_determinants = true;
2330
    }
2331
    if (ToCompute.test(ConeProperty::TriangulationSize)) {
2332
        FC.do_triangulation = true;
2333
    }
2334
    if (ToCompute.test(ConeProperty::NoSubdivision)) {
2335
        FC.use_bottom_points = false;
2336
    }
2337
    if (ToCompute.test(ConeProperty::Deg1Elements)) {
2338
        FC.do_deg1_elements = true;
2339
    }
2340
    if (ToCompute.test(ConeProperty::StanleyDec)) {
2341
        FC.do_Stanley_dec = true;
2342
    }
2343
    if (ToCompute.test(ConeProperty::Approximate)
2344
     && ToCompute.test(ConeProperty::Deg1Elements)) {
2345
        FC.do_approximation = true;
2346
        FC.do_deg1_elements = true;
2347
    }
2348
    if (ToCompute.test(ConeProperty::DefaultMode)) {
2349
        FC.do_default_mode = true;
2350
    }
2351
    if (ToCompute.test(ConeProperty::BottomDecomposition)) {
2352
        FC.do_bottom_dec = true;
2353
    }
2354
    if (ToCompute.test(ConeProperty::NoBottomDec)) {
2355
        FC.suppress_bottom_dec = true;
2356
    }
2357
    if (ToCompute.test(ConeProperty::KeepOrder)) {
2358
        FC.keep_order = true;
2359
    }
2360
    if (ToCompute.test(ConeProperty::ClassGroup)) {
2361
        FC.do_class_group=true;
2362
    }
2363
    if (ToCompute.test(ConeProperty::ModuleRank)) {
2364
        FC.do_module_rank=true;
2365
    }
2366
    
2367
    if (ToCompute.test(ConeProperty::HSOP)) {
2368
        FC.do_hsop=true;
2369
    }
2370
    
2371
    /* Give extra data to FC */
2372
    if ( isComputed(ConeProperty::ExtremeRays) ) {
2373
        FC.Extreme_Rays_Ind = ExtremeRaysIndicator;
2374
        FC.is_Computed.set(ConeProperty::ExtremeRays);
2375
    }
2376
    
2377
    /* if(isComputed(ConeProperty::Deg1Elements)){
2378
        Matrix<IntegerFC> Deg1Converted;
2379
        BasisChangePointed.convert_to_sublattice(Deg1Converted, Deg1Elements);
2380
        for(size_t i=0;i<Deg1Elements.nr_of_rows();++i)
2381
            FC.Deg1_Elements.push_back(Deg1Converted[i]);
2382
        FC.is_Computed.set(ConeProperty::Deg1Elements); 
2383
    }
2384
    
2385
    if(isComputed(ConeProperty::HilbertBasis)){
2386
        Matrix<IntegerFC> HBConverted;
2387
        BasisChangePointed.convert_to_sublattice(HBConverted, HilbertBasis);
2388
        for(size_t i=0;i<HilbertBasis.nr_of_rows();++i)
2389
            FC.Deg1_Elements.push_back(HBConverted[i]);
2390
        FC.is_Computed.set(ConeProperty::HilbertBasis); 
2391
    }*/
2392
    
2393
    if (ExcludedFaces.nr_of_rows()!=0) {
2394
        BasisChangePointed.convert_to_sublattice_dual(FC.ExcludedFaces, ExcludedFaces);
2395
    }
2396
    if (isComputed(ConeProperty::ExcludedFaces)) {
2397
        FC.is_Computed.set(ConeProperty::ExcludedFaces);
2398
    }
2399

2400
    if (inhomogeneous){
2401
        BasisChangePointed.convert_to_sublattice_dual_no_div(FC.Truncation, Dehomogenization);
2402
    }
2403
    if ( Grading.size()>0 ) {  // IMPORTANT: Truncation must be set before Grading
2404
        BasisChangePointed.convert_to_sublattice_dual(FC.Grading, Grading);
2405
        if(isComputed(ConeProperty::Grading) ){    // is grading positive?
2406
            FC.is_Computed.set(ConeProperty::Grading);
2407
            /*if (inhomogeneous)
2408
                FC.find_grading_inhom();
2409
            FC.set_degrees();*/
2410
        }
2411
    }
2412

2413
    if (SupportHyperplanes.nr_of_rows()!=0) {
2414
        BasisChangePointed.convert_to_sublattice_dual(FC.Support_Hyperplanes, SupportHyperplanes);
2415
   }
2416
    if (isComputed(ConeProperty::SupportHyperplanes)){
2417
        FC.is_Computed.set(ConeProperty::SupportHyperplanes);
2418
        FC.do_all_hyperplanes = false;
2419
    }
2420

2421
    if(ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid)){
2422
        FC.do_module_gens_intcl=true;
2423
    }
2424
    
2425
    if(is_approximation)
2426
        give_data_of_approximated_cone_to(FC);
2427

2428
    /* do the computation */
2429
    
2430
    try {     
2431
        try {
2432
            FC.compute();
2433
        } catch (const NotIntegrallyClosedException& ) {
2434
        }
2435
        is_Computed.set(ConeProperty::Sublattice);
2436
        // make sure we minimize the excluded faces if requested
2437
        if(ToCompute.test(ConeProperty::ExcludedFaces) || ToCompute.test(ConeProperty::SupportHyperplanes)) {
2438
            FC.prepare_inclusion_exclusion();
2439
        }
2440
        extract_data(FC);
2441
        if(isComputed(ConeProperty::IsPointed) && pointed)
2442
            is_Computed.set(ConeProperty::MaximalSubspace);
2443
    } catch(const NonpointedException& ) {
2444
        is_Computed.set(ConeProperty::Sublattice);
2445
        extract_data(FC);
2446
        if(verbose){
2447
            verboseOutput() << "Cone not pointed. Restarting computation." << endl;
2448
        }
2449
        FC=Full_Cone<IntegerFC>(Matrix<IntegerFC>(1)); // to kill the old FC (almost)
2450
        Matrix<Integer> Dual_Gen;
2451
        Dual_Gen=BasisChangePointed.to_sublattice_dual(SupportHyperplanes);
2452
        Sublattice_Representation<Integer> Pointed(Dual_Gen,true); // sublattice of the dual lattice
2453
        BasisMaxSubspace = BasisChangePointed.from_sublattice(Pointed.getEquationsMatrix());
2454
        check_vanishing_of_grading_and_dehom();
2455
        BasisChangePointed.compose_dual(Pointed);
2456
        is_Computed.set(ConeProperty::MaximalSubspace);        
2457
        // now we get the basis of the maximal subspace
2458
        pointed = (BasisMaxSubspace.nr_of_rows() == 0);
2459
        is_Computed.set(ConeProperty::IsPointed);
2460
        compute_full_cone<IntegerFC>(ToCompute);           
2461
    }
2462
}
2463

2464

2465
template<typename Integer>
2466
void Cone<Integer>::compute_generators() {
2467
    //create Generators from SupportHyperplanes
2468
    if (!isComputed(ConeProperty::Generators) && (SupportHyperplanes.nr_of_rows()!=0 ||inhomogeneous)) {
2469
        if (verbose) {
2470
            verboseOutput() << "Computing extreme rays as support hyperplanes of the dual cone:" << endl;
2471
        }
2472
        if (change_integer_type) {
2473
            try {
2474
                compute_generators_inner<MachineInteger>();
2475
            } catch(const ArithmeticException& e) {
2476
                if (verbose) {
2477
                    verboseOutput() << e.what() << endl;
2478
                    verboseOutput() << "Restarting with a bigger type." << endl;
2479
                }
2480
                compute_generators_inner<Integer>();
2481
            }
2482
        } else {
2483
            compute_generators_inner<Integer>();
2484
        }
2485
    }
2486
    assert(isComputed(ConeProperty::Generators));
2487
}
2488

2489
template<typename Integer>
2490
template<typename IntegerFC>
2491
void Cone<Integer>::compute_generators_inner() {
2492
    
2493
    Matrix<Integer> Dual_Gen;
2494
    Dual_Gen=BasisChangePointed.to_sublattice_dual(SupportHyperplanes);
2495
    // first we take the quotient of the efficient sublattice modulo the maximal subspace
2496
    Sublattice_Representation<Integer> Pointed(Dual_Gen,true); // sublattice of the dual space
2497

2498
    // now we get the basis of the maximal subspace
2499
    if(!isComputed(ConeProperty::MaximalSubspace)){
2500
        BasisMaxSubspace = BasisChangePointed.from_sublattice(Pointed.getEquationsMatrix());
2501
        check_vanishing_of_grading_and_dehom();
2502
        is_Computed.set(ConeProperty::MaximalSubspace);
2503
    }
2504
    if(!isComputed(ConeProperty::IsPointed)){
2505
        pointed = (BasisMaxSubspace.nr_of_rows() == 0);
2506
        is_Computed.set(ConeProperty::IsPointed);
2507
    }
2508
    BasisChangePointed.compose_dual(Pointed); // primal cone now pointed, may not yet be full dimensional
2509

2510
    // restrict the supphyps to efficient sublattice and push to quotient mod subspace
2511
    Matrix<IntegerFC> Dual_Gen_Pointed;
2512
    BasisChangePointed.convert_to_sublattice_dual(Dual_Gen_Pointed, SupportHyperplanes);    
2513
    Full_Cone<IntegerFC> Dual_Cone(Dual_Gen_Pointed);
2514
    Dual_Cone.verbose=verbose;
2515
    Dual_Cone.do_extreme_rays=true; // we try to find them, need not exist
2516
    try {     
2517
        Dual_Cone.dualize_cone();
2518
    } catch(const NonpointedException& ){}; // we don't mind if the dual cone is not pointed
2519
    
2520
    if (Dual_Cone.isComputed(ConeProperty::SupportHyperplanes)) {
2521
        //get the extreme rays of the primal cone
2522
        BasisChangePointed.convert_from_sublattice(Generators,
2523
                          Dual_Cone.getSupportHyperplanes());
2524
        is_Computed.set(ConeProperty::Generators);
2525
        
2526
        //get minmal set of support_hyperplanes if possible
2527
        if (Dual_Cone.isComputed(ConeProperty::ExtremeRays)) {            
2528
            Matrix<IntegerFC> Supp_Hyp = Dual_Cone.getGenerators().submatrix(Dual_Cone.getExtremeRays());
2529
            BasisChangePointed.convert_from_sublattice_dual(SupportHyperplanes, Supp_Hyp);
2530
            SupportHyperplanes.sort_lex();
2531
            is_Computed.set(ConeProperty::SupportHyperplanes);
2532
        }
2533
        
2534
        // now the final transformations
2535
        // only necessary if the basis changes computed so far do not make the cone full-dimensional
2536
        // this is equaivalent to the dual cone bot being pointed
2537
        if(!(Dual_Cone.isComputed(ConeProperty::IsPointed) && Dual_Cone.isPointed())){
2538
            // first to full-dimensional pointed
2539
            Matrix<Integer> Help;
2540
            Help=BasisChangePointed.to_sublattice(Generators); // sublattice of the primal space
2541
            Sublattice_Representation<Integer> PointedHelp(Help,true);
2542
            BasisChangePointed.compose(PointedHelp);
2543
            // second to efficient sublattice
2544
            if(BasisMaxSubspace.nr_of_rows()==0){  // primal cone is pointed and we can copy
2545
                BasisChange=BasisChangePointed;
2546
            }
2547
            else{
2548
                Help=BasisChange.to_sublattice(Generators);
2549
                Help.append(BasisChange.to_sublattice(BasisMaxSubspace));
2550
                Sublattice_Representation<Integer> EmbHelp(Help,true); // sublattice of the primal space
2551
                compose_basis_change(EmbHelp);
2552
            }
2553
        }
2554
        is_Computed.set(ConeProperty::Sublattice); // will not be changed anymore
2555

2556
        checkGrading();
2557
        // compute grading, so that it is also known if nothing else is done afterwards
2558
        if (!isComputed(ConeProperty::Grading) && !inhomogeneous) {
2559
            // Generators = ExtremeRays
2560
            vector<Integer> lf = BasisChangePointed.to_sublattice(Generators).find_linear_form();
2561
            if (lf.size() == BasisChange.getRank()) {
2562
                vector<Integer> test_lf=BasisChange.from_sublattice_dual(lf);
2563
                if(Generators.nr_of_rows()==0 || v_scalar_product(Generators[0],test_lf)==1)
2564
                    setGrading(test_lf);
2565
            }
2566
        }
2567
        setWeights();
2568
        set_extreme_rays(vector<bool>(Generators.nr_of_rows(),true)); // here since they get sorted
2569
        is_Computed.set(ConeProperty::ExtremeRays);
2570
    }
2571
}
2572

2573

2574
template<typename Integer>
2575
void Cone<Integer>::compute_dual(ConeProperties& ToCompute) {
2576

2577
    ToCompute.reset(is_Computed);
2578
    if (ToCompute.none() || !( ToCompute.test(ConeProperty::Deg1Elements)
2579
                            || ToCompute.test(ConeProperty::HilbertBasis))) {
2580
        return;
2581
    }
2582

2583
    if (change_integer_type) {
2584
        try {
2585
            compute_dual_inner<MachineInteger>(ToCompute);
2586
        } catch(const ArithmeticException& e) {
2587
            if (verbose) {
2588
                verboseOutput() << e.what() << endl;
2589
                verboseOutput() << "Restarting with a bigger type." << endl;
2590
            }
2591
            change_integer_type = false;
2592
        }
2593
    }
2594
    if (!change_integer_type) {
2595
        compute_dual_inner<Integer>(ToCompute);
2596
    }
2597
    ToCompute.reset(ConeProperty::DualMode);
2598
    ToCompute.reset(is_Computed);
2599
    // if (ToCompute.test(ConeProperty::DefaultMode) && ToCompute.goals().none()) {
2600
    //    ToCompute.reset(ConeProperty::DefaultMode);
2601
    // }
2602
}
2603

2604
template<typename Integer>
2605
vector<Sublattice_Representation<Integer> > MakeSubAndQuot(const Matrix<Integer>& Gen,
2606
                                        const Matrix<Integer>& Ker){
2607
    vector<Sublattice_Representation<Integer> > Result;                                        
2608
    Matrix<Integer> Help=Gen;
2609
    Help.append(Ker);
2610
    Sublattice_Representation<Integer> Sub(Help,true);
2611
    Sublattice_Representation<Integer> Quot=Sub;
2612
    if(Ker.nr_of_rows()>0){
2613
        Matrix<Integer> HelpQuot=Sub.to_sublattice(Ker).kernel();   // kernel here to be interpreted as subspace of the dual
2614
                                                                    // namely the linear forms vanishing on Ker
2615
        Sublattice_Representation<Integer> SubToQuot(HelpQuot,true); // sublattice of the dual
2616
        Quot.compose_dual(SubToQuot);
2617
    }
2618
    Result.push_back(Sub);
2619
    Result.push_back(Quot);
2620
    
2621
    return Result;    
2622
}
2623

2624
template<typename Integer>
2625
template<typename IntegerFC>
2626
void Cone<Integer>::compute_dual_inner(ConeProperties& ToCompute) {
2627

2628
    bool do_only_Deg1_Elements = ToCompute.test(ConeProperty::Deg1Elements)
2629
                                 && !ToCompute.test(ConeProperty::HilbertBasis);
2630

2631
    if(isComputed(ConeProperty::Generators) && SupportHyperplanes.nr_of_rows()==0){
2632
        if (verbose) {
2633
            verboseOutput()<<  "Computing support hyperplanes for the dual mode:"<< endl;
2634
        }
2635
        ConeProperties Dualize;
2636
        Dualize.set(ConeProperty::SupportHyperplanes);
2637
        Dualize.set(ConeProperty::ExtremeRays);
2638
        recursive_compute(Dualize);
2639
    }
2640
    
2641
    bool do_extreme_rays_first = false;
2642
    if (!isComputed(ConeProperty::ExtremeRays)) {
2643
        if (do_only_Deg1_Elements && Grading.size()==0)
2644
            do_extreme_rays_first = true;
2645
        else if ( (do_only_Deg1_Elements || inhomogeneous) &&
2646
                   ( naked_dual
2647
                 || ToCompute.test(ConeProperty::ExtremeRays)
2648
                 || ToCompute.test(ConeProperty::SupportHyperplanes)
2649
                 || ToCompute.test(ConeProperty::Sublattice) ) )
2650
            do_extreme_rays_first = true;
2651
    }
2652

2653
    if (do_extreme_rays_first) {
2654
        if (verbose) {
2655
            verboseOutput() << "Computing extreme rays for the dual mode:"<< endl;
2656
        }
2657
        compute_generators();   // computes extreme rays, but does not find grading !
2658
    }
2659

2660
    if(do_only_Deg1_Elements && Grading.size()==0){
2661
        vector<Integer> lf= Generators.submatrix(ExtremeRaysIndicator).find_linear_form_low_dim();
2662
        if(Generators.nr_of_rows()==0 || (lf.size()==dim && v_scalar_product(Generators[0],lf)==1))
2663
            setGrading(lf);
2664
        else{
2665
            throw BadInputException("Need grading to compute degree 1 elements and cannot find one.");
2666
        }
2667
    }
2668

2669
    if (SupportHyperplanes.nr_of_rows()==0 && !isComputed(ConeProperty::SupportHyperplanes)) {
2670
        throw FatalException("Could not get SupportHyperplanes.");
2671
    }
2672

2673
    Matrix<IntegerFC> Inequ_on_Ker;
2674
    BasisChangePointed.convert_to_sublattice_dual(Inequ_on_Ker,SupportHyperplanes);
2675
     
2676
    vector<IntegerFC> Truncation;
2677
    if(inhomogeneous){
2678
        BasisChangePointed.convert_to_sublattice_dual_no_div(Truncation, Dehomogenization);
2679
    }
2680
    if (do_only_Deg1_Elements) {
2681
        // in this case the grading acts as truncation and it is a NEW inequality
2682
        BasisChangePointed.convert_to_sublattice_dual(Truncation, Grading);
2683
    }
2684

2685
    Cone_Dual_Mode<IntegerFC> ConeDM(Inequ_on_Ker, Truncation); // Inequ_on_Ker is NOT const
2686
    Inequ_on_Ker=Matrix<IntegerFC>(0,0);  // destroy it
2687
    ConeDM.verbose=verbose;
2688
    ConeDM.inhomogeneous=inhomogeneous;
2689
    ConeDM.do_only_Deg1_Elements=do_only_Deg1_Elements;
2690
    if(isComputed(ConeProperty::Generators))
2691
        BasisChangePointed.convert_to_sublattice(ConeDM.Generators, Generators);
2692
    if(isComputed(ConeProperty::ExtremeRays))
2693
        ConeDM.ExtremeRaysInd=ExtremeRaysIndicator;
2694
    ConeDM.hilbert_basis_dual();
2695
    
2696
    if(!isComputed(ConeProperty::MaximalSubspace)){
2697
        BasisChangePointed.convert_from_sublattice(BasisMaxSubspace,ConeDM.BasisMaxSubspace);
2698
        check_vanishing_of_grading_and_dehom(); // all this must be done here because to_sublattice may kill it
2699
    }
2700

2701
    if (!isComputed(ConeProperty::Sublattice) || !isComputed(ConeProperty::MaximalSubspace)){
2702
        if(!(do_only_Deg1_Elements || inhomogeneous)) {
2703
            // At this point we still have BasisChange==BasisChangePointed
2704
            // now we can pass to a pointed full-dimensional cone
2705
            
2706
            vector<Sublattice_Representation<IntegerFC> > BothRepFC=MakeSubAndQuot
2707
                        (ConeDM.Generators,ConeDM.BasisMaxSubspace);
2708
            if(!BothRepFC[0].IsIdentity())        
2709
                BasisChange.compose(Sublattice_Representation<Integer>(BothRepFC[0]));
2710
            is_Computed.set(ConeProperty::Sublattice);
2711
            if(!BothRepFC[1].IsIdentity())
2712
                BasisChangePointed.compose(Sublattice_Representation<Integer>(BothRepFC[1]));
2713
            ConeDM.to_sublattice(BothRepFC[1]);
2714
        }
2715
    }
2716
    
2717
    is_Computed.set(ConeProperty::MaximalSubspace); // NOT EARLIER !!!!
2718
    
2719
    
2720
    // create a Full_Cone out of ConeDM
2721
    Full_Cone<IntegerFC> FC(ConeDM);
2722
    FC.verbose=verbose;
2723
    // Give extra data to FC
2724
    if (Grading.size()>0) {
2725
        BasisChangePointed.convert_to_sublattice_dual(FC.Grading, Grading);
2726
        if(isComputed(ConeProperty::Grading))
2727
            FC.is_Computed.set(ConeProperty::Grading);
2728
    }
2729
    if(inhomogeneous)
2730
        BasisChangePointed.convert_to_sublattice_dual_no_div(FC.Truncation, Dehomogenization);
2731
    FC.do_class_group=ToCompute.test(ConeProperty::ClassGroup);
2732
    FC.dual_mode();
2733
    extract_data(FC);
2734
}
2735

2736
//---------------------------------------------------------------------------
2737

2738
template<typename Integer>
2739
Integer Cone<Integer>::compute_primary_multiplicity() {
2740
    if (change_integer_type) {
2741
        try {
2742
            return compute_primary_multiplicity_inner<MachineInteger>();
2743
        } catch(const ArithmeticException& e) {
2744
            if (verbose) {
2745
                verboseOutput() << e.what() << endl;
2746
                verboseOutput() << "Restarting with a bigger type." << endl;
2747
            }
2748
            change_integer_type = false;
2749
        }
2750
    }
2751
    return compute_primary_multiplicity_inner<Integer>();
2752
}
2753

2754
//---------------------------------------------------------------------------
2755

2756
template<typename Integer>
2757
template<typename IntegerFC>
2758
Integer Cone<Integer>::compute_primary_multiplicity_inner() {
2759
    Matrix<IntegerFC> Ideal(0,dim-1);
2760
    vector<IntegerFC> help(dim-1);
2761
    for(size_t i=0;i<Generators.nr_of_rows();++i){ // select ideal generators
2762
        if(Generators[i][dim-1]==1){
2763
            for(size_t j=0;j<dim-1;++j)
2764
                convert(help[j],Generators[i][j]);
2765
            Ideal.append(help);
2766
        }
2767
    }
2768
    Full_Cone<IntegerFC> IdCone(Ideal,false);
2769
    IdCone.do_bottom_dec=true;
2770
    IdCone.do_determinants=true;
2771
    IdCone.compute();
2772
    return convertTo<Integer>(IdCone.detSum);
2773
}
2774

2775
//---------------------------------------------------------------------------
2776

2777
template<typename Integer>
2778
template<typename IntegerFC>
2779
void Cone<Integer>::extract_data(Full_Cone<IntegerFC>& FC) {
2780
    //this function extracts ALL available data from the Full_Cone
2781
    //even if it was in Cone already <- this may change
2782
    //it is possible to delete the data in Full_Cone after extracting it
2783

2784
    if(verbose) {
2785
        verboseOutput() << "transforming data..."<<flush;
2786
    }
2787
    
2788
    if (FC.isComputed(ConeProperty::Generators)) {
2789
        BasisChangePointed.convert_from_sublattice(Generators,FC.getGenerators());
2790
        is_Computed.set(ConeProperty::Generators);
2791
    }
2792
    
2793
    if (FC.isComputed(ConeProperty::IsPointed) && !isComputed(ConeProperty::IsPointed)) {
2794
        pointed = FC.isPointed();
2795
        if(pointed)
2796
            is_Computed.set(ConeProperty::MaximalSubspace);
2797
        is_Computed.set(ConeProperty::IsPointed);
2798
    }    
2799
    
2800
    if (FC.isComputed(ConeProperty::Grading)) {
2801
        if (Grading.size()==0) {
2802
            BasisChangePointed.convert_from_sublattice_dual(Grading, FC.getGrading());
2803
        }
2804
        is_Computed.set(ConeProperty::Grading);
2805
        setWeights();
2806
        //compute denominator of Grading
2807
        if(BasisChangePointed.getRank()!=0){
2808
            vector<Integer> test_grading = BasisChangePointed.to_sublattice_dual_no_div(Grading);
2809
            GradingDenom=v_make_prime(test_grading);
2810
        }
2811
        else
2812
            GradingDenom=1; 
2813
        is_Computed.set(ConeProperty::GradingDenom);
2814
    }
2815
        
2816
    if (FC.isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)) { // must precede extreme rays
2817
        BasisChangePointed.convert_from_sublattice(ModuleGeneratorsOverOriginalMonoid, FC.getModuleGeneratorsOverOriginalMonoid());
2818
        ModuleGeneratorsOverOriginalMonoid.sort_by_weights(WeightsGrad,GradAbs);
2819
        is_Computed.set(ConeProperty::ModuleGeneratorsOverOriginalMonoid);
2820
    }
2821

2822
    if (FC.isComputed(ConeProperty::ExtremeRays)) {
2823
        set_extreme_rays(FC.getExtremeRays());
2824
    }
2825
    if (FC.isComputed(ConeProperty::SupportHyperplanes)) {
2826
        /* if (inhomogeneous) {
2827
            // remove irrelevant support hyperplane 0 ... 0 1
2828
            vector<IntegerFC> irr_hyp_subl;
2829
            BasisChangePointed.convert_to_sublattice_dual(irr_hyp_subl, Dehomogenization); 
2830
            FC.Support_Hyperplanes.remove_row(irr_hyp_subl);
2831
        } */
2832
        // BasisChangePointed.convert_from_sublattice_dual(SupportHyperplanes, FC.getSupportHyperplanes());
2833
        extract_supphyps(FC);
2834
        if(inhomogeneous && FC.dim<dim){ // make inequality for the inhomogeneous variable appear as dehomogenization
2835
            vector<Integer> dehom_restricted=BasisChangePointed.to_sublattice_dual(Dehomogenization);
2836
            for(size_t i=0;i<SupportHyperplanes.nr_of_rows();++i){
2837
                if(dehom_restricted==BasisChangePointed.to_sublattice_dual(SupportHyperplanes[i])){
2838
                    SupportHyperplanes[i]=Dehomogenization;
2839
                    break;
2840
                }
2841
            }
2842
        }
2843
        SupportHyperplanes.sort_lex();
2844
        is_Computed.set(ConeProperty::SupportHyperplanes);
2845
    }
2846
    if (FC.isComputed(ConeProperty::TriangulationSize)) {
2847
        TriangulationSize = FC.totalNrSimplices;
2848
        triangulation_is_nested = FC.triangulation_is_nested;
2849
        triangulation_is_partial= FC.triangulation_is_partial;
2850
        is_Computed.set(ConeProperty::TriangulationSize);
2851
        is_Computed.set(ConeProperty::IsTriangulationPartial);
2852
        is_Computed.set(ConeProperty::IsTriangulationNested);
2853
        is_Computed.reset(ConeProperty::Triangulation);
2854
        Triangulation.clear(); // to get rid of a previously computed triangulation
2855
    }
2856
    if (FC.isComputed(ConeProperty::TriangulationDetSum)) {
2857
        convert(TriangulationDetSum, FC.detSum);
2858
        is_Computed.set(ConeProperty::TriangulationDetSum);
2859
    }
2860
    
2861
    if (FC.isComputed(ConeProperty::Triangulation)) {
2862
        size_t tri_size = FC.Triangulation.size();
2863
        FC.Triangulation.sort(compareKeys<IntegerFC>); // necessary to make triangulation unique
2864
        Triangulation = vector< pair<vector<key_t>, Integer> >(tri_size);
2865
        if(FC.isComputed(ConeProperty::ConeDecomposition))
2866
            OpenFacets.resize(tri_size);
2867
        SHORTSIMPLEX<IntegerFC> simp;
2868
        for (size_t i = 0; i<tri_size; ++i) {
2869
            simp = FC.Triangulation.front();
2870
            Triangulation[i].first.swap(simp.key);
2871
            /* sort(Triangulation[i].first.begin(), Triangulation[i].first.end()); -- no longer allowed here because of ConeDecomposition. Done in full_cone.cpp, transfer_triangulation_to top */
2872
            if (FC.isComputed(ConeProperty::TriangulationDetSum))
2873
                convert(Triangulation[i].second, simp.vol);
2874
            else
2875
                Triangulation[i].second = 0;
2876
            if(FC.isComputed(ConeProperty::ConeDecomposition))
2877
                OpenFacets[i].swap(simp.Excluded);
2878
            FC.Triangulation.pop_front();
2879
        }
2880
        if(FC.isComputed(ConeProperty::ConeDecomposition))
2881
            is_Computed.set(ConeProperty::ConeDecomposition);
2882
        is_Computed.set(ConeProperty::Triangulation);
2883
    }
2884

2885
    if (FC.isComputed(ConeProperty::StanleyDec)) {
2886
        StanleyDec.clear();
2887
        StanleyDec.splice(StanleyDec.begin(),FC.StanleyDec);
2888
        // At present, StanleyDec not sorted here
2889
        is_Computed.set(ConeProperty::StanleyDec);
2890
    }
2891

2892
    if (FC.isComputed(ConeProperty::InclusionExclusionData)) {
2893
        InExData.clear();
2894
        InExData.reserve(FC.InExCollect.size());
2895
        map<boost::dynamic_bitset<>, long>::iterator F;
2896
        vector<key_t> key;
2897
        for (F=FC.InExCollect.begin(); F!=FC.InExCollect.end(); ++F) {
2898
            key.clear();
2899
            for (size_t i=0;i<FC.nr_gen;++i) {
2900
                if (F->first.test(i)) {
2901
                    key.push_back(i);
2902
                }
2903
            }
2904
            InExData.push_back(make_pair(key,F->second));
2905
        }
2906
        is_Computed.set(ConeProperty::InclusionExclusionData);
2907
    }
2908
    if (FC.isComputed(ConeProperty::RecessionRank) && isComputed(ConeProperty::MaximalSubspace)) {
2909
        recession_rank = FC.level0_dim+BasisMaxSubspace.nr_of_rows();
2910
        is_Computed.set(ConeProperty::RecessionRank);
2911
        if (getRank() == recession_rank) {
2912
            affine_dim = -1;
2913
        } else {
2914
            affine_dim = getRank()-1;
2915
        }
2916
        is_Computed.set(ConeProperty::AffineDim);
2917
    }
2918
    if (FC.isComputed(ConeProperty::ModuleRank)) {
2919
        module_rank = FC.getModuleRank();
2920
        is_Computed.set(ConeProperty::ModuleRank);
2921
    }
2922
    if (FC.isComputed(ConeProperty::Multiplicity)) {
2923
        if(!inhomogeneous) {
2924
            multiplicity = FC.getMultiplicity();
2925
            is_Computed.set(ConeProperty::Multiplicity);
2926
        } else if (isComputed(ConeProperty::ModuleRank)) {
2927
            multiplicity = FC.getMultiplicity()*module_rank;
2928
            is_Computed.set(ConeProperty::Multiplicity);
2929
        }
2930
    }
2931
    if (FC.isComputed(ConeProperty::WitnessNotIntegrallyClosed)) {
2932
        BasisChangePointed.convert_from_sublattice(WitnessNotIntegrallyClosed,FC.Witness);
2933
        is_Computed.set(ConeProperty::WitnessNotIntegrallyClosed);
2934
        integrally_closed = false;
2935
        is_Computed.set(ConeProperty::IsIntegrallyClosed);
2936
    }
2937
    if (FC.isComputed(ConeProperty::HilbertBasis)) {
2938
        if (inhomogeneous) {
2939
            // separate (capped) Hilbert basis to the Hilbert basis of the level 0 cone
2940
            // and the module generators in level 1
2941
            HilbertBasis = Matrix<Integer>(0,dim);
2942
            ModuleGenerators = Matrix<Integer>(0,dim);
2943
            typename list< vector<IntegerFC> >::const_iterator FCHB(FC.Hilbert_Basis.begin());
2944
            vector<Integer> tmp;
2945
            for (; FCHB != FC.Hilbert_Basis.end(); ++FCHB) {
2946
                
2947
                INTERRUPT_COMPUTATION_BY_EXCEPTION
2948
                
2949
                BasisChangePointed.convert_from_sublattice(tmp,*FCHB);
2950
                if (v_scalar_product(tmp,Dehomogenization) == 0) { // Hilbert basis element of the cone at level 0
2951
                    HilbertBasis.append(tmp);
2952
                } else {              // module generator
2953
                    ModuleGenerators.append(tmp);
2954
                }
2955
            }
2956
            ModuleGenerators.sort_by_weights(WeightsGrad,GradAbs);
2957
            is_Computed.set(ConeProperty::ModuleGenerators);
2958
        } else { // homogeneous
2959
            HilbertBasis = Matrix<Integer>(0,dim);
2960
            typename list< vector<IntegerFC> >::const_iterator FCHB(FC.Hilbert_Basis.begin());
2961
            vector<Integer> tmp;
2962
            for (; FCHB != FC.Hilbert_Basis.end(); ++FCHB) {
2963
                BasisChangePointed.convert_from_sublattice(tmp,*FCHB);                
2964
                HilbertBasis.append(tmp);
2965
            }
2966
        }
2967
        HilbertBasis.sort_by_weights(WeightsGrad,GradAbs);
2968
        is_Computed.set(ConeProperty::HilbertBasis);
2969
    }
2970
    if (FC.isComputed(ConeProperty::Deg1Elements)) {
2971
        Deg1Elements = Matrix<Integer>(0,dim);
2972
        typename list< vector<IntegerFC> >::const_iterator DFC(FC.Deg1_Elements.begin());
2973
        vector<Integer> tmp;
2974
        for (; DFC != FC.Deg1_Elements.end(); ++DFC) {
2975
            
2976
            INTERRUPT_COMPUTATION_BY_EXCEPTION
2977
            
2978
            BasisChangePointed.convert_from_sublattice(tmp,*DFC);                
2979
            Deg1Elements.append(tmp);
2980
        }
2981
        Deg1Elements.sort_by_weights(WeightsGrad,GradAbs);
2982
        is_Computed.set(ConeProperty::Deg1Elements);
2983
    }
2984
    if (FC.isComputed(ConeProperty::HilbertSeries)) {
2985
        long save_nr_coeff_quasipol=HSeries.get_nr_coeff_quasipol(); // Full_Cone does not compute the quasipolynomial
2986
        HSeries = FC.Hilbert_Series;
2987
        HSeries.set_nr_coeff_quasipol(save_nr_coeff_quasipol);
2988
        is_Computed.set(ConeProperty::HilbertSeries);
2989
    }
2990
    if (FC.isComputed(ConeProperty::HSOP)) {
2991
        is_Computed.set(ConeProperty::HSOP);
2992
    }
2993
    if (FC.isComputed(ConeProperty::IsDeg1ExtremeRays)) {
2994
        deg1_extreme_rays = FC.isDeg1ExtremeRays();
2995
        is_Computed.set(ConeProperty::IsDeg1ExtremeRays);
2996
    }
2997
    if (FC.isComputed(ConeProperty::ExcludedFaces)) {
2998
        BasisChangePointed.convert_from_sublattice_dual(ExcludedFaces, FC.getExcludedFaces());
2999
        ExcludedFaces.sort_lex();
3000
        is_Computed.set(ConeProperty::ExcludedFaces);
3001
    }
3002

3003
    if (FC.isComputed(ConeProperty::IsDeg1HilbertBasis)) {
3004
        deg1_hilbert_basis = FC.isDeg1HilbertBasis();
3005
        is_Computed.set(ConeProperty::IsDeg1HilbertBasis);
3006
    }
3007
    if (FC.isComputed(ConeProperty::ClassGroup)) {
3008
        convert(ClassGroup, FC.ClassGroup);
3009
        is_Computed.set(ConeProperty::ClassGroup);
3010
    }
3011
    
3012
    /* if (FC.isComputed(ConeProperty::MaximalSubspace) && 
3013
                                   !isComputed(ConeProperty::MaximalSubspace)) {
3014
        BasisChangePointed.convert_from_sublattice(BasisMaxSubspace, FC.Basis_Max_Subspace);
3015
        check_vanishing_of_grading_and_dehom();
3016
        is_Computed.set(ConeProperty::MaximalSubspace);
3017
    }*/
3018

3019
    check_integrally_closed();
3020

3021
    if (verbose) {
3022
        verboseOutput() << " done." <<endl;
3023
    }
3024
}
3025

3026
//---------------------------------------------------------------------------
3027
template<typename Integer>
3028
template<typename IntegerFC>
3029
void Cone<Integer>::extract_supphyps(Full_Cone<IntegerFC>& FC) {
3030
        BasisChangePointed.convert_from_sublattice_dual(SupportHyperplanes, FC.getSupportHyperplanes());
3031
}
3032

3033
template<typename Integer>
3034
void Cone<Integer>::extract_supphyps(Full_Cone<Integer>& FC) {
3035
    if(BasisChangePointed.IsIdentity())
3036
        swap(SupportHyperplanes,FC.Support_Hyperplanes);
3037
    else
3038
        SupportHyperplanes=BasisChangePointed.from_sublattice_dual(FC.getSupportHyperplanes());
3039
}
3040

3041

3042
//---------------------------------------------------------------------------
3043

3044
template<typename Integer>
3045
void Cone<Integer>::check_integrally_closed() {
3046
    if (!isComputed(ConeProperty::OriginalMonoidGenerators)
3047
            || isComputed(ConeProperty::IsIntegrallyClosed)
3048
            || !isComputed(ConeProperty::HilbertBasis) || inhomogeneous)
3049
        return;
3050

3051
    unit_group_index=1;
3052
    if(BasisMaxSubspace.nr_of_rows()>0)
3053
        compute_unit_group_index();
3054
    is_Computed.set(ConeProperty::UnitGroupIndex);
3055
    if (index > 1 || HilbertBasis.nr_of_rows() > OriginalMonoidGenerators.nr_of_rows()
3056
            || unit_group_index>1) {
3057
        integrally_closed = false;
3058
        is_Computed.set(ConeProperty::IsIntegrallyClosed);
3059
        return;
3060
    } 
3061
    find_witness();
3062
}
3063

3064
//---------------------------------------------------------------------------
3065

3066
template<typename Integer>
3067
void Cone<Integer>::compute_unit_group_index() {
3068
    assert(isComputed(ConeProperty::MaximalSubspace));
3069
    // we want to compute in the maximal linear subspace
3070
    Sublattice_Representation<Integer> Sub(BasisMaxSubspace,true);
3071
    Matrix<Integer> origens_in_subspace(0,dim);
3072

3073
    // we must collect all original generetors that lie in the maximal subspace 
3074

3075
    for(size_t i=0;i<OriginalMonoidGenerators.nr_of_rows();++i){
3076
        size_t j;
3077
        for(j=0;j<SupportHyperplanes.nr_of_rows();++j){
3078
                if(v_scalar_product(OriginalMonoidGenerators[i],SupportHyperplanes[j])!=0)
3079
                    break;
3080
        }
3081
        if(j==SupportHyperplanes.nr_of_rows())
3082
            origens_in_subspace.append(OriginalMonoidGenerators[i]);
3083
    }
3084
    Matrix<Integer> M=Sub.to_sublattice(origens_in_subspace);
3085
    unit_group_index= M.full_rank_index();
3086
    // cout << "Unit group index " << unit_group_index;
3087
}
3088

3089
//---------------------------------------------------------------------------
3090

3091
template<typename Integer>
3092
void Cone<Integer>::find_witness() {
3093
    if (!isComputed(ConeProperty::OriginalMonoidGenerators)
3094
            || inhomogeneous) {
3095
        // no original monoid defined
3096
        throw NotComputableException(ConeProperties(ConeProperty::WitnessNotIntegrallyClosed));
3097
    }
3098
    if (isComputed(ConeProperty::IsIntegrallyClosed) && integrally_closed) {
3099
        // original monoid is integrally closed
3100
        throw NotComputableException(ConeProperties(ConeProperty::WitnessNotIntegrallyClosed));
3101
    }
3102
    if (isComputed(ConeProperty::WitnessNotIntegrallyClosed)
3103
            || !isComputed(ConeProperty::HilbertBasis) )
3104
        return;
3105

3106
    long nr_gens = OriginalMonoidGenerators.nr_of_rows();
3107
    long nr_hilb = HilbertBasis.nr_of_rows();
3108
    // if the cone is not pointed, we have to check it on the quotion
3109
    Matrix<Integer> gens_quot;
3110
    Matrix<Integer> hilb_quot;
3111
    if (!pointed) {
3112
        gens_quot = BasisChangePointed.to_sublattice(OriginalMonoidGenerators);
3113
        hilb_quot = BasisChangePointed.to_sublattice(HilbertBasis);
3114
    }
3115
    Matrix<Integer>& gens = pointed ? OriginalMonoidGenerators : gens_quot;
3116
    Matrix<Integer>& hilb = pointed ? HilbertBasis : hilb_quot;
3117
    integrally_closed = true;
3118
    typename list< vector<Integer> >::iterator h;
3119
    for (long h = 0; h < nr_hilb; ++h) {
3120
        integrally_closed = false;
3121
        for (long i = 0; i < nr_gens; ++i) {
3122
            if (hilb[h] == gens[i]) {
3123
                integrally_closed = true;
3124
                break;
3125
            }
3126
        }
3127
        if (!integrally_closed) {
3128
            WitnessNotIntegrallyClosed = HilbertBasis[h];
3129
            is_Computed.set(ConeProperty::WitnessNotIntegrallyClosed);
3130
            break;
3131
        }
3132
    }
3133
    is_Computed.set(ConeProperty::IsIntegrallyClosed);
3134
}
3135

3136
//---------------------------------------------------------------------------
3137

3138
template<typename Integer>
3139
void Cone<Integer>::set_original_monoid_generators(const Matrix<Integer>& Input) {
3140
    if (!isComputed(ConeProperty::OriginalMonoidGenerators)) {
3141
        OriginalMonoidGenerators = Input;
3142
        is_Computed.set(ConeProperty::OriginalMonoidGenerators);
3143
    }
3144
    // Generators = Input;
3145
    // is_Computed.set(ConeProperty::Generators);
3146
    Matrix<Integer> M=BasisChange.to_sublattice(Input);
3147
    index=M.full_rank_index();
3148
    is_Computed.set(ConeProperty::InternalIndex);
3149
}
3150

3151
//---------------------------------------------------------------------------
3152

3153
template<typename Integer>
3154
void Cone<Integer>::set_extreme_rays(const vector<bool>& ext) {
3155
    assert(ext.size() == Generators.nr_of_rows());
3156
    ExtremeRaysIndicator=ext;
3157
    vector<bool> choice=ext;
3158
    if (inhomogeneous) {
3159
        // separate extreme rays to rays of the level 0 cone
3160
        // and the verticies of the polyhedron, which are in level >=1
3161
        size_t nr_gen = Generators.nr_of_rows();
3162
        vector<bool> VOP(nr_gen);
3163
        for (size_t i=0; i<nr_gen; i++) {
3164
            if (ext[i] && v_scalar_product(Generators[i],Dehomogenization) != 0) {
3165
                VOP[i] = true;
3166
                choice[i]=false;
3167
            }
3168
        }
3169
        VerticesOfPolyhedron=Generators.submatrix(VOP).sort_by_weights(WeightsGrad,GradAbs);
3170
        is_Computed.set(ConeProperty::VerticesOfPolyhedron);
3171
    }
3172
    ExtremeRays=Generators.submatrix(choice);
3173
    if(inhomogeneous && !isComputed(ConeProperty::AffineDim) && isComputed(ConeProperty::MaximalSubspace)){
3174
        size_t level0_dim=ExtremeRays.max_rank_submatrix_lex().size();
3175
        recession_rank = level0_dim+BasisMaxSubspace.nr_of_rows();
3176
        is_Computed.set(ConeProperty::RecessionRank);
3177
        if (getRank() == recession_rank) {
3178
            affine_dim = -1;
3179
        } else {
3180
            affine_dim = getRank()-1;
3181
        }
3182
        is_Computed.set(ConeProperty::AffineDim);
3183
        
3184
    }
3185
    if(isComputed(ConeProperty::ModuleGeneratorsOverOriginalMonoid)){  // not possible in inhomogeneous case
3186
        Matrix<Integer> ExteEmbedded=BasisChangePointed.to_sublattice(ExtremeRays);
3187
        for(size_t i=0;i<ExteEmbedded.nr_of_rows();++i)
3188
            v_make_prime(ExteEmbedded[i]);
3189
        ExteEmbedded.remove_duplicate_and_zero_rows();
3190
        ExtremeRays=BasisChangePointed.from_sublattice(ExteEmbedded);
3191
    }
3192
    ExtremeRays.sort_by_weights(WeightsGrad,GradAbs);
3193
    is_Computed.set(ConeProperty::ExtremeRays);
3194
}
3195

3196
//---------------------------------------------------------------------------
3197

3198
template<typename Integer>
3199
void Cone<Integer>::compute_vertices_float(ConeProperties& ToCompute) {
3200
    if(!ToCompute.test(ConeProperty::VerticesFloat) || isComputed(ConeProperty::VerticesFloat))
3201
        return;
3202
    if(!isComputed(ConeProperty::ExtremeRays))
3203
        throw NotComputableException("VerticesFloat not computable without extreme rays");
3204
    if(inhomogeneous && !isComputed(ConeProperty::VerticesOfPolyhedron))
3205
        throw NotComputableException("VerticesFloat not computable in the inhomogeneous case without vertices");
3206
    if(!inhomogeneous && !isComputed(ConeProperty::Grading))
3207
        throw NotComputableException("VerticesFloat not computable in the homogeneous case without a grading");
3208
    if(inhomogeneous)
3209
        convert(VerticesFloat, VerticesOfPolyhedron);
3210
    else
3211
        convert(VerticesFloat, ExtremeRays);
3212
    vector<nmz_float> norm;
3213
    if(inhomogeneous)
3214
        convert(norm,Dehomogenization);
3215
    else{
3216
       convert( norm,Grading);
3217
       nmz_float GD=1.0/convertTo<double>(GradingDenom);
3218
       v_scalar_multiplication(norm,GD);
3219
    }
3220
    for(size_t i=0;i<VerticesFloat.nr_of_rows();++i){
3221
        nmz_float den=1.0/v_scalar_product(VerticesFloat[i],norm);
3222
        v_scalar_multiplication(VerticesFloat[i],den);        
3223
    }
3224
    is_Computed.set(ConeProperty::VerticesFloat);
3225
}
3226

3227
//---------------------------------------------------------------------------
3228

3229
template<typename Integer>
3230
void Cone<Integer>::complete_sublattice_comp(ConeProperties& ToCompute) {
3231
    
3232
    if(!isComputed(ConeProperty::Sublattice))
3233
        return;
3234
    is_Computed.set(ConeProperty::Rank);
3235
    if(ToCompute.test(ConeProperty::Equations)){
3236
        BasisChange.getEquationsMatrix(); // just to force computation, ditto below
3237
        is_Computed.set(ConeProperty::Equations);
3238
    }
3239
    if(ToCompute.test(ConeProperty::Congruences) || ToCompute.test(ConeProperty::ExternalIndex)){
3240
        BasisChange.getCongruencesMatrix();
3241
        BasisChange.getExternalIndex();
3242
        is_Computed.set(ConeProperty::Congruences);
3243
        is_Computed.set(ConeProperty::ExternalIndex);
3244
    }
3245
}
3246

3247
template<typename Integer>
3248
void Cone<Integer>::complete_HilbertSeries_comp(ConeProperties& ToCompute) {
3249
    if(!isComputed(ConeProperty::HilbertSeries))
3250
        return;
3251
    if(ToCompute.test(ConeProperty::HilbertQuasiPolynomial))
3252
        HSeries.computeHilbertQuasiPolynomial();
3253
    if(HSeries.isHilbertQuasiPolynomialComputed())
3254
        is_Computed.set(ConeProperty::HilbertQuasiPolynomial);
3255
        
3256
    // in the case that HS was computed but not HSOP, we need to compute hsop
3257
    if(ToCompute.test(ConeProperty::HSOP) && !isComputed(ConeProperty::HSOP)){
3258
        // we need generators and support hyperplanes to compute hsop
3259
        Matrix<Integer> FC_gens;
3260
        Matrix<Integer> FC_hyps;
3261
        BasisChangePointed.convert_to_sublattice(FC_gens,Generators);
3262
        Full_Cone<Integer> FC(FC_gens);
3263
        FC.Extreme_Rays_Ind = ExtremeRaysIndicator;
3264
        FC.Grading = Grading;
3265
        FC.inhomogeneous = inhomogeneous;
3266
        // TRUNCATION?
3267
        BasisChangePointed.convert_to_sublattice_dual(FC.Support_Hyperplanes,SupportHyperplanes);
3268
        FC.compute_hsop();
3269
        HSeries.setHSOPDenom(FC.Hilbert_Series.getHSOPDenom());
3270
        HSeries.compute_hsop_num();
3271
    }    
3272
}
3273

3274
//---------------------------------------------------------------------------
3275
template<typename Integer>
3276
void Cone<Integer>::set_project(string name){
3277
    project=name;
3278
}
3279

3280
template<typename Integer>
3281
void Cone<Integer>::set_output_dir(string name){
3282
    output_dir=name;
3283
}
3284

3285
template<typename Integer>
3286
void Cone<Integer>::set_nmz_call(const string& path){
3287
    nmz_call=path;
3288
}
3289

3290
template<typename Integer>
3291
void Cone<Integer>::setPolynomial(string poly){
3292
    IntData=IntegrationData(poly);
3293
}
3294

3295
template<typename Integer>
3296
void Cone<Integer>::setNrCoeffQuasiPol(long nr_coeff){
3297
    IntData.set_nr_coeff_quasipol(nr_coeff);
3298
    HSeries.set_nr_coeff_quasipol(nr_coeff);
3299
}
3300

3301
bool executable(string command){
3302
//n check whether "command --version" cam be executed
3303

3304
    command +=" --version";
3305
    string dev0= " > /dev/null";
3306
#ifdef _WIN32 //for 32 and 64 bit windows
3307
    dev0=" > NUL:";
3308
#endif
3309
    command+=dev0;
3310
    if(system(command.c_str())==0)
3311
        return true;
3312
    else
3313
        return false;
3314
}
3315

3316
string command(const string& original_call, const string& to_replace, const string& by_this){
3317
// in the original call we replace the program name to_replace by by_this
3318
// we try variants with and without "lt-" preceding the names of executables
3319
// since libtools may have inserted "lt-" before the original name
3320

3321
    string copy=original_call;
3322
    // cout << "CALL " << original_call << endl;
3323
    string search_lt="lt-"+to_replace;
3324
    long length=to_replace.size();
3325
    size_t found;
3326
    found = copy.rfind(search_lt);
3327
    if (found==std::string::npos) {
3328
        found = copy.rfind(to_replace);
3329
        if (found==std::string::npos){
3330
            throw FatalException("Call "+ copy +" of "  +to_replace+" does not contain " +to_replace); 
3331
        }
3332
    }
3333
    else{
3334
            length+=3; //name includes lt-
3335
    }
3336
    string test_path=copy.replace (found,length,by_this);
3337
    // cout << "TEST " << test_path << endl;
3338
    if(executable(test_path)) // first without lt-
3339
        return test_path;
3340
    copy=original_call;
3341
    string by_this_with_lt="lt-"+by_this; /// now with lt-
3342
    test_path=copy.replace (found,length,by_this_with_lt);
3343
    // cout << "TEST " << test_path << endl;
3344
    if(executable(test_path))
3345
        return test_path;
3346
    return ""; // no executable found
3347
}
3348

3349
//---------------------------------------------------------------------------
3350
template<typename Integer>
3351
void Cone<Integer>::try_symmetrization(ConeProperties& ToCompute) {
3352

3353
    if(ToCompute.test(ConeProperty::NoSymmetrization))
3354
        return;
3355
    
3356
    if(!(ToCompute.test(ConeProperty::Symmetrize) || ToCompute.test(ConeProperty::HilbertSeries) ||
3357
               ToCompute.test(ConeProperty::Multiplicity)))
3358
        return;
3359
    
3360
    if(inhomogeneous || nr_latt_gen>0|| nr_cone_gen>0 || lattice_ideal_input || Grading.size() < dim){
3361
        if(ToCompute.test(ConeProperty::Symmetrize))
3362
            throw BadInputException("Symmetrization not posible with the given input"); 
3363
        else
3364
            return;
3365
    }
3366
    
3367
#ifndef NMZ_COCOA    
3368
    if(project==""){
3369
        if(ToCompute.test(ConeProperty::Symmetrize)){
3370
            throw BadInputException("Symmetrization via libnormaliz not possible without CoCoALib");
3371
        }
3372
        else
3373
            return;
3374
    }
3375
#endif
3376
    
3377
    Matrix<Integer> AllConst=ExcludedFaces;
3378
    size_t nr_excl = AllConst.nr_of_rows();    
3379
    AllConst. append(Equations);
3380
    size_t nr_equ=AllConst.nr_of_rows()-nr_excl;
3381
    vector<bool> unit_vector(dim,false);
3382
    for(size_t i=0;i<Inequalities.nr_of_rows();++i){
3383
        size_t nr_nonzero=0;
3384
        size_t nonzero_coord;
3385
        bool is_unit_vector=true;
3386
        for(size_t j=0;j<dim;++j){
3387
            if(Inequalities[i][j]==0)
3388
                continue;
3389
            if(nr_nonzero>0 || Inequalities[i][j]!=1){ // not a sign inequality
3390
                is_unit_vector=false;                
3391
                break;    
3392
            }
3393
            nr_nonzero++;
3394
            nonzero_coord=j;
3395
        }
3396
        if(!is_unit_vector)
3397
            AllConst.append(Inequalities[i]);
3398
        else
3399
            unit_vector[nonzero_coord]=true;    
3400
    }
3401
    
3402
    size_t nr_inequ=AllConst.nr_of_rows()-nr_equ-nr_excl;
3403
    
3404
    for(size_t i=0;i<dim;++i)
3405
        if(!unit_vector[i]){
3406
            if(ToCompute.test(ConeProperty::Symmetrize))
3407
                throw BadInputException("Symmetrization not possible: Not all sign inequalities in input");
3408
            else
3409
                return;
3410
        }
3411
    
3412
    for(size_t i=0;i<Congruences.nr_of_rows();++i){
3413
        vector<Integer> help=Congruences[i];
3414
        help.resize(dim);
3415
        AllConst.append(help);
3416
    }
3417
    // now we have collected all constraints and cehcked the existence of the sign inequalities
3418
    
3419
    
3420
    AllConst.append(Grading);
3421
    
3422
    /* AllConst.pretty_print(cout);
3423
    cout << "----------------------" << endl;
3424
    cout << nr_excl << " " << nr_equ << " " << nr_inequ << endl; */
3425
    
3426
    AllConst=AllConst.transpose();
3427
    
3428
    map< vector<Integer>, size_t > classes;
3429
    typename map< vector<Integer>, size_t >::iterator C;
3430

3431
    for(size_t j=0;j<AllConst.nr_of_rows();++j){
3432
        C=classes.find(AllConst[j]);
3433
        if(C!=classes.end())
3434
            C->second++;
3435
        else
3436
            classes.insert(pair<vector<Integer>, size_t>(AllConst[j],1));
3437
    }
3438
    
3439
    vector<size_t> multiplicities;
3440
    Matrix<Integer> SymmConst(0,AllConst.nr_of_columns());
3441
    
3442
    for(C=classes.begin();C!=classes.end();++C){
3443
            multiplicities.push_back(C->second);
3444
            SymmConst.append(C->first);
3445
    }
3446
    SymmConst=SymmConst.transpose();
3447
    
3448
    vector<Integer> SymmGrad=SymmConst[SymmConst.nr_of_rows()-1];
3449
    
3450
    if(verbose){
3451
        verboseOutput() << "Embedding dimension of symmetrized cone = " << SymmGrad.size() << endl;
3452
    }
3453
    
3454
    if(SymmGrad.size() > 2*dim/3){
3455
        if(!ToCompute.test(ConeProperty::Symmetrize)){
3456
            return;
3457
        }
3458
    }
3459
    
3460
    /* compute_generators(); // we must protect against the zero cone
3461
    if(getRank()==0)
3462
        return; */
3463
    
3464
    Matrix<Integer> SymmInequ(0,SymmConst.nr_of_columns());
3465
    Matrix<Integer> SymmEqu(0,SymmConst.nr_of_columns());
3466
    Matrix<Integer> SymmCong(0,SymmConst.nr_of_columns());
3467
    Matrix<Integer> SymmExcl(0,SymmConst.nr_of_columns());
3468
 
3469
    for(size_t i=0;i<nr_excl;++i)
3470
        SymmExcl.append(SymmConst[i]);        
3471
    for(size_t i=nr_excl;i<nr_excl+nr_equ;++i)
3472
        SymmEqu.append(SymmConst[i]);    
3473
    for(size_t i=nr_excl+nr_equ;i<nr_excl+nr_equ+nr_inequ;++i)
3474
        SymmInequ.append(SymmConst[i]);    
3475
    for(size_t i=nr_excl+nr_equ+nr_inequ;i<SymmConst.nr_of_rows()-1;++i){
3476
        SymmCong.append(SymmConst[i]);
3477
        SymmCong[SymmCong.nr_of_rows()-1].push_back(Congruences[i-(nr_inequ+nr_equ)][dim]); // restore modulus
3478
    }
3479

3480
    string polynomial;
3481
    
3482
    for(size_t i=0;i<multiplicities.size();++i){
3483
        for(size_t j=1;j<multiplicities[i];++j)
3484
            polynomial+="(x["+to_string((unsigned long long) i+1)+"]+"+to_string((unsigned long long)j)+")*";
3485
        
3486
    }
3487
    polynomial+="1";
3488
    mpz_class fact=1;
3489
    for(size_t i=0;i<multiplicities.size();++i){
3490
        for(size_t j=1;j<multiplicities[i];++j)
3491
            fact*=j;        
3492
    }
3493
    polynomial+="/"+fact.get_str()+";";
3494

3495
#ifdef NMZ_COCOA
3496
    
3497
    map< InputType, Matrix<Integer> > SymmInput;
3498
    SymmInput[InputType::inequalities]=SymmInequ;
3499
    SymmInput[InputType::equations]=SymmEqu;
3500
    SymmInput[InputType::congruences]=SymmCong;
3501
    SymmInput[InputType::excluded_faces]=SymmExcl;
3502
    Matrix<Integer> GradMat(0,SymmGrad.size());
3503
    GradMat.append(SymmGrad);
3504
    SymmInput[InputType::grading]=GradMat;
3505
    Matrix<Integer> SymmNonNeg(0,SymmGrad.size());
3506
    vector<Integer>  NonNeg(SymmGrad.size(),1);
3507
    SymmNonNeg.append(NonNeg);
3508
    SymmInput[InputType::signs]=SymmNonNeg;
3509
    SymmCone=new Cone<Integer>(SymmInput);
3510
    SymmCone->setPolynomial(polynomial);
3511
    SymmCone->setNrCoeffQuasiPol(HSeries.get_nr_coeff_quasipol());
3512
    SymmCone->HSeries.set_period_bounded(HSeries.get_period_bounded());
3513
    SymmCone->setVerbose(verbose);
3514
    ConeProperties SymmToCompute;
3515
    SymmToCompute.set(ConeProperty::SupportHyperplanes);
3516
    SymmToCompute.set(ConeProperty::WeightedEhrhartSeries,ToCompute.test(ConeProperty::HilbertSeries));
3517
    SymmToCompute.set(ConeProperty::VirtualMultiplicity,ToCompute.test(ConeProperty::Multiplicity));
3518
    SymmToCompute.set(ConeProperty::BottomDecomposition,ToCompute.test(ConeProperty::BottomDecomposition));
3519
    SymmCone->compute(SymmToCompute);
3520
    if(SymmCone->isComputed(ConeProperty::WeightedEhrhartSeries)){
3521
        HSeries=SymmCone->getWeightedEhrhartSeries().first;
3522
        is_Computed.set(ConeProperty::HilbertSeries);
3523
    }
3524
    if(SymmCone->isComputed(ConeProperty::VirtualMultiplicity)){
3525
        multiplicity=SymmCone->getVirtualMultiplicity();
3526
        is_Computed.set(ConeProperty::Multiplicity);
3527
    }
3528
    is_Computed.set(ConeProperty::Symmetrize);
3529
    return;
3530
    
3531
#endif
3532

3533
}
3534

3535

3536

3537
template<typename Integer>
3538
void integrate(Cone<Integer>& C, const bool do_virt_mult);
3539

3540
template<typename Integer>
3541
void generalizedEhrhartSeries(Cone<Integer>& C);
3542

3543
template<typename Integer>
3544
void Cone<Integer>::compute_integral (ConeProperties& ToCompute){
3545
    if(isComputed(ConeProperty::Integral) || !ToCompute.test(ConeProperty::Integral))
3546
        return;
3547
    if(IntData.getPolynomial()=="")
3548
        throw BadInputException("Polynomial weight missing");
3549
#ifdef NMZ_COCOA
3550
    if(getRank()==0)
3551
        getIntData().setIntegral(0);
3552
    else
3553
    integrate<Integer>(*this,false);
3554
    is_Computed.set(ConeProperty::Integral);
3555
#endif
3556
}
3557
    
3558
template<typename Integer>
3559
void Cone<Integer>::compute_virt_mult(ConeProperties& ToCompute){
3560
    if(isComputed(ConeProperty::VirtualMultiplicity) || !ToCompute.test(ConeProperty::VirtualMultiplicity))
3561
        return;
3562
    if(IntData.getPolynomial()=="")
3563
        throw BadInputException("Polynomial weight missing");
3564
#ifdef NMZ_COCOA
3565
    if(getRank()==0)
3566
        getIntData().setVirtualMultiplicity(0);
3567
    else
3568
        integrate<Integer>(*this,true);
3569
    is_Computed.set(ConeProperty::VirtualMultiplicity);
3570
#endif
3571
}
3572

3573
template<typename Integer>
3574
void Cone<Integer>::compute_weighted_Ehrhart(ConeProperties& ToCompute){
3575
    if(isComputed(ConeProperty::WeightedEhrhartSeries) || !ToCompute.test(ConeProperty::WeightedEhrhartSeries))
3576
        return;
3577
    if(IntData.getPolynomial()=="")
3578
        throw BadInputException("Polynomial weight missing");    
3579
    /* if(getRank()==0)
3580
        throw NotComputableException("WeightedEhrhartSeries not computed in dimenison 0");*/
3581
#ifdef NMZ_COCOA
3582
    generalizedEhrhartSeries(*this);
3583
    is_Computed.set(ConeProperty::WeightedEhrhartSeries);
3584
    if(getIntData().isWeightedEhrhartQuasiPolynomialComputed()){
3585
        is_Computed.set(ConeProperty::WeightedEhrhartQuasiPolynomial);
3586
        is_Computed.set(ConeProperty::VirtualMultiplicity);
3587
    }
3588
#endif
3589
}
3590
//---------------------------------------------------------------------------
3591
template<typename Integer>
3592
bool Cone<Integer>::get_verbose (){
3593
    return verbose;
3594
}
3595

3596
//---------------------------------------------------------------------------
3597
template<typename Integer>
3598
void Cone<Integer>::NotComputable (string message){
3599
    if(!default_mode)
3600
        throw NotComputableException(message);
3601
}
3602

3603
//---------------------------------------------------------------------------
3604
template<typename Integer>
3605
void Cone<Integer>::check_Gorenstein(ConeProperties&  ToCompute){
3606
    
3607
    if(!ToCompute.test(ConeProperty::IsGorenstein) || isComputed(ConeProperty::IsGorenstein))
3608
        return;
3609
    if(!isComputed(ConeProperty::SupportHyperplanes))
3610
        recursive_compute(ConeProperty::SupportHyperplanes);
3611
    if(!isComputed(ConeProperty::MaximalSubspace))
3612
        recursive_compute(ConeProperty::MaximalSubspace);
3613
    
3614
    if(dim==0){
3615
        Gorenstein=true;
3616
        is_Computed.set(ConeProperty::IsGorenstein);
3617
        GeneratorOfInterior=vector<Integer> (dim,0);
3618
        is_Computed.set(ConeProperty::GeneratorOfInterior);
3619
        return;        
3620
    }
3621
    Matrix<Integer> TransfSupps=BasisChangePointed.to_sublattice_dual(SupportHyperplanes);
3622
    assert(TransfSupps.nr_of_rows()>0);
3623
    Gorenstein=false;
3624
    vector<Integer> TransfIntGen = TransfSupps.find_linear_form();
3625
    if(TransfIntGen.size()!=0 && v_scalar_product(TransfIntGen,TransfSupps[0])==1){
3626
        Gorenstein=true;
3627
        GeneratorOfInterior=BasisChangePointed.from_sublattice(TransfIntGen);
3628
        is_Computed.set(ConeProperty::GeneratorOfInterior);
3629
    }
3630
    is_Computed.set(ConeProperty::IsGorenstein);
3631
}
3632

3633
//---------------------------------------------------------------------------
3634
template<typename Integer>
3635
template<typename IntegerFC>
3636
void Cone<Integer>::give_data_of_approximated_cone_to(Full_Cone<IntegerFC>& FC){
3637
    
3638
    // *this is the approximatING cone. The support hyperplanes and equations of the approximatED 
3639
    // cone are given to the Full_Cone produced from *this so that the superfluous points can
3640
    // bre sorted out as early as possible.
3641
    
3642
    assert(is_approximation);
3643
    assert(ApproximatedCone->inhomogeneous ||  ApproximatedCone->getGradingDenom()==1); // in case we generalize later
3644
    
3645
    FC.is_global_approximation=true;
3646
    // FC.is_approximation=true; At present not allowed. Only used for approximation within Full_Cone
3647
    
3648
    // We must distinguish zwo cases: Approximated->Grading_Is_Coordinate or it is not
3649

3650
    // If it is not:
3651
    // The first coordinate in *this is the degree given by the grading
3652
    // in ApproximatedCone. We disregard it by setting the first coordinate
3653
    // of the grading, inequalities and equations to 0, and then have 0 followed
3654
    // by the grading, equations and inequalities resp. of ApproximatedCone.
3655
    
3656
    vector<Integer> help_g;
3657
    if(ApproximatedCone->inhomogeneous)
3658
        help_g=ApproximatedCone->Dehomogenization;
3659
    else
3660
        help_g=ApproximatedCone->Grading;
3661
    
3662
    if(ApproximatedCone->Grading_Is_Coordinate){
3663
        swap(help_g[0],help_g[ApproximatedCone->GradingCoordinate]);
3664
        BasisChangePointed.convert_to_sublattice_dual_no_div(FC.Subcone_Grading,help_g); 
3665
    }        
3666
    else{            
3667
        vector<Integer> help(help_g.size()+1);
3668
        help[0]=0;
3669
        for(size_t j=0;j<help_g.size();++j)
3670
            help[j+1]=help_g[j];
3671
        BasisChangePointed.convert_to_sublattice_dual_no_div(FC.Subcone_Grading,help);
3672
    }
3673
    
3674
    Matrix<Integer> Eq=ApproximatedCone->BasisChangePointed.getEquationsMatrix();
3675
    FC.Subcone_Equations=Matrix<IntegerFC>(Eq.nr_of_rows(),BasisChangePointed.getRank());
3676
    if(ApproximatedCone->Grading_Is_Coordinate){
3677
        Eq.exchange_columns(0,ApproximatedCone->GradingCoordinate);
3678
        BasisChangePointed.convert_to_sublattice_dual(FC.Subcone_Equations,Eq);
3679
    }
3680
    else{
3681
        for(size_t i=0;i<Eq.nr_of_rows();++i){
3682
            vector<Integer> help(Eq.nr_of_columns()+1,0);
3683
            for(size_t j=0;j<Eq.nr_of_columns();++j)
3684
                help[j+1]=Eq[i][j];
3685
            BasisChangePointed.convert_to_sublattice_dual(FC.Subcone_Equations[i], help);       
3686
        }
3687
    }
3688
    
3689
    Matrix<Integer> Supp=ApproximatedCone->SupportHyperplanes;
3690
    FC.Subcone_Support_Hyperplanes=Matrix<IntegerFC>(Supp.nr_of_rows(),BasisChangePointed.getRank());
3691
    
3692
    if(ApproximatedCone->Grading_Is_Coordinate){
3693
        Supp.exchange_columns(0,ApproximatedCone->GradingCoordinate);
3694
        BasisChangePointed.convert_to_sublattice_dual(FC.Subcone_Support_Hyperplanes,Supp);
3695
    }
3696
    else{
3697
        for(size_t i=0;i<Supp.nr_of_rows();++i){
3698
            vector<Integer> help(Supp.nr_of_columns()+1,0);
3699
            for(size_t j=0;j<Supp.nr_of_columns();++j)
3700
                help[j+1]=Supp[i][j];
3701
            BasisChangePointed.convert_to_sublattice_dual(FC.Subcone_Support_Hyperplanes[i], help);       
3702
        }
3703
    }
3704
}
3705

3706
//---------------------------------------------------------------------------
3707
template<typename Integer>
3708
void Cone<Integer>::try_approximation_or_projection(ConeProperties& ToCompute){
3709
    
3710
    if((ToCompute.test(ConeProperty::NoProjection) && !ToCompute.test(ConeProperty::Approximate))
3711
           || ToCompute.test(ConeProperty::DualMode) || ToCompute.test(ConeProperty::PrimalMode)
3712
    )
3713
        return;
3714
    
3715
    if(ToCompute.test(ConeProperty::ModuleGeneratorsOverOriginalMonoid))
3716
        return;
3717
    
3718
    if(!inhomogeneous && (  !ToCompute.test(ConeProperty::Deg1Elements) 
3719
                         || ToCompute.test(ConeProperty::HilbertBasis)
3720
                         || ToCompute.test(ConeProperty::HilbertSeries)
3721
                         )                        
3722
      )
3723
        return;
3724
    
3725
    if(inhomogeneous && !ToCompute.test(ConeProperty::HilbertBasis) )
3726
        return;
3727
    
3728
    bool is_parallelotope=false;
3729
    if(!ToCompute.test(ConeProperty::Approximate))
3730
        is_parallelotope=check_parallelotope();
3731
    if(verbose && is_parallelotope)
3732
        verboseOutput() << "Polyhedron is parallelotope" << endl;
3733
    
3734
    if(is_parallelotope){
3735
        SupportHyperplanes.remove_row(Dehomogenization);
3736
        is_Computed.set(ConeProperty::SupportHyperplanes);
3737
        is_Computed.set(ConeProperty::MaximalSubspace);
3738
        is_Computed.set(ConeProperty::Sublattice);
3739
        pointed=true;
3740
        is_Computed.set(ConeProperty::IsPointed);
3741
    }
3742
   
3743
    ConeProperties NeededHere;
3744
    NeededHere.set(ConeProperty::SupportHyperplanes);
3745
    NeededHere.set(ConeProperty::Sublattice);
3746
    NeededHere.set(ConeProperty::MaximalSubspace);
3747
    if(!inhomogeneous)
3748
        NeededHere.set(ConeProperty::Grading);
3749
    recursive_compute(NeededHere);
3750
    
3751
    if(!is_parallelotope && !ToCompute.test(ConeProperty::Approximate)){ // we try again
3752
        is_parallelotope=check_parallelotope();
3753
        if(is_parallelotope){
3754
            SupportHyperplanes.remove_row(Dehomogenization);
3755
            is_Computed.set(ConeProperty::SupportHyperplanes);
3756
            is_Computed.set(ConeProperty::MaximalSubspace);
3757
            is_Computed.set(ConeProperty::Sublattice);
3758
            pointed=true;
3759
            is_Computed.set(ConeProperty::IsPointed);
3760
        }
3761
    }
3762
    
3763
    // cout << "PPPPP " << is_parallelotope << endl;
3764
    
3765
    if(!inhomogeneous && !isComputed(ConeProperty::Grading))
3766
        return;
3767
    
3768
     if(!inhomogeneous && ToCompute.test(ConeProperty::Approximate) && GradingDenom!=1)
3769
        return;
3770
    
3771
    if(!pointed || BasisChangePointed.getRank()==0)
3772
        return;
3773
    
3774
    if(inhomogeneous){
3775
        for(size_t i=0;i<Generators.nr_of_rows();++i){
3776
            if(v_scalar_product(Generators[i],Dehomogenization)==0){
3777
                if(ToCompute.test(ConeProperty::Approximate) || ToCompute.test(ConeProperty::Projection))
3778
                    throw NotComputableException("Approximation or Projection not applicable to unbounded polyhedra");
3779
                else
3780
                    return;
3781
            }                    
3782
        }        
3783
    }
3784
    
3785
    if(inhomogeneous){ // exclude that dehoogenization has a gcd > 1
3786
        vector<Integer> test_dehom=BasisChange.to_sublattice_dual_no_div(Dehomogenization);
3787
        if(v_make_prime(test_dehom)!=1)
3788
            return;        
3789
    }
3790
    
3791
    // ****************************************************************
3792
    //
3793
    // NOTE: THE FIRST COORDINATE IS (OR WILL BE MADE) THE GRADING !!!!
3794
    //
3795
    // ****************************************************************
3796
    
3797
    vector<Integer> GradForApprox;
3798
    if(!inhomogeneous)
3799
        GradForApprox=Grading;
3800
    else{
3801
        GradForApprox=Dehomogenization;
3802
        GradingDenom=1;
3803
    }
3804
    
3805
    Grading_Is_Coordinate=false;
3806
    size_t nr_nonzero=0;
3807
    for(size_t i=0;i<dim;++i){
3808
        if(GradForApprox[i]!=0){
3809
            nr_nonzero++;
3810
            GradingCoordinate=i;
3811
        }
3812
    }
3813
    if(nr_nonzero==1){
3814
        if(GradForApprox[GradingCoordinate]==1)
3815
            Grading_Is_Coordinate=true;        
3816
    }
3817
    
3818
    Matrix<Integer> GradGen;
3819
    if(Grading_Is_Coordinate){
3820
        if(!ToCompute.test(ConeProperty::Approximate)){
3821
            GradGen=Generators;
3822
            GradGen.exchange_columns(0,GradingCoordinate); // we swap it into the first coordinate
3823
        }
3824
        else{ // we swap the grading into the first coordinate and approximate
3825
            GradGen.resize(0,dim);
3826
            for(size_t i=0;i<Generators.nr_of_rows();++i){
3827
                vector<Integer> gg=Generators[i];
3828
                swap(gg[0],gg[GradingCoordinate]);
3829
                list<vector<Integer> > approx;
3830
                approx_simplex(gg,approx,1);
3831
                GradGen.append(Matrix<Integer>(approx));
3832
            }    
3833
        }           
3834
    }    
3835
    else{ // to avoid coordinate trabnsformations, we prepend the degree as the first coordinate
3836
        GradGen.resize(0,dim+1); 
3837
        for(size_t i=0;i<Generators.nr_of_rows();++i){
3838
            vector<Integer> gg(dim+1);
3839
            for(size_t j=0;j<dim;++j)
3840
                gg[j+1]=Generators[i][j];
3841
            gg[0]=v_scalar_product(Generators[i],GradForApprox);
3842
            // cout << gg;
3843
            if(ToCompute.test(ConeProperty::Approximate)){
3844
                list<vector<Integer> > approx;
3845
                approx_simplex(gg,approx,1);
3846
                GradGen.append(Matrix<Integer>(approx));
3847
            }
3848
            else
3849
                GradGen.append(gg);            
3850
        }
3851
    }
3852
    
3853
    Matrix<Integer> Raw(0,GradGen.nr_of_columns()); // result is returned in this matrix
3854
        
3855
    if(ToCompute.test(ConeProperty::Approximate)){
3856
        if(verbose)
3857
            verboseOutput() << "Computing lattice points by approximation" << endl;
3858
        Cone<Integer> HelperCone(InputType::cone,GradGen);
3859
        HelperCone. ApproximatedCone=&(*this); // we will pass this infornation to the Full_Cone that computes the lattice points.
3860
        HelperCone.is_approximation=true;  // It allows us to discard points outside *this as quickly as possible
3861
        HelperCone.compute(ConeProperty::Deg1Elements,ConeProperty::PrimalMode);
3862
        Raw=HelperCone.getDeg1ElementsMatrix();        
3863
    }
3864
    else{
3865
        if(verbose)
3866
             verboseOutput() << "Computing lattice points by project-and-lift" << endl;
3867
        Matrix<Integer> Supps, Equs;
3868
        if(Grading_Is_Coordinate){
3869
            Supps=getSupportHyperplanesMatrix();
3870
            Supps.exchange_columns(0,GradingCoordinate);
3871
            Equs=BasisChange.getEquationsMatrix();
3872
            Equs.exchange_columns(0,GradingCoordinate);
3873
        }
3874
        else{
3875
            vector<vector<Integer> > SuppHelp=getSupportHyperplanes();
3876
            insert_column<Integer>(SuppHelp,0,0);
3877
            Supps=Matrix<Integer>(SuppHelp);
3878
            vector<vector<Integer> > EqusHelp=BasisChange.getEquations();
3879
            if(EqusHelp.size()>0){
3880
                insert_column<Integer>(EqusHelp,0,0);
3881
                Equs=Matrix<Integer>(EqusHelp);
3882
            }
3883
            else
3884
                Equs.resize(0,dim+1);
3885
            vector<Integer> ExtraEqu(Equs.nr_of_columns());
3886
            ExtraEqu[0]=-1;
3887
            for(size_t i=0;i<Grading.size();++i)
3888
                ExtraEqu[i+1]=Grading[i];
3889
            Equs.append(ExtraEqu);
3890
        }
3891
        Supps.append(Equs);  // we must add the equations as pairs of inequalities
3892
        Equs.scalar_multiplication(-1);
3893
        Supps.append(Equs);
3894
        project_and_lift(Raw, GradGen,Supps,ToCompute.test(ConeProperty::ProjectionFloat));        
3895
    }
3896
    
3897
    HilbertBasis=Matrix<Integer>(0,dim);
3898
    Deg1Elements=Matrix<Integer>(0,dim);
3899
    ModuleGenerators=Matrix<Integer>(0,dim);
3900
    
3901
    if(Grading_Is_Coordinate)
3902
        Raw.exchange_columns(0,GradingCoordinate);
3903
    
3904
    Matrix<Integer> Cong=BasisChange.getCongruences();
3905
    
3906
    if(Grading_Is_Coordinate && Cong.nr_of_rows()==0){
3907
        if(inhomogeneous)
3908
            ModuleGenerators.swap(Raw);
3909
        else
3910
            Deg1Elements.swap(Raw);
3911
    }
3912
    else{
3913
        for(size_t i=0;i<Raw.nr_of_rows();++i){
3914
            vector<Integer> rr;
3915
            if(Grading_Is_Coordinate){
3916
                swap(rr,Raw[i]);
3917
            }
3918
            else{
3919
                rr.resize(dim); // remove the prepended grading
3920
                for(size_t j=0;j<dim;++j)
3921
                    rr[j]=Raw[i][j+1];
3922
            }
3923
            bool not_in=false;
3924
            for(size_t k=0;k<Cong.nr_of_rows();++k) {
3925
                if(v_scalar_product_vectors_unequal_lungth(rr,Cong[k]) % Cong[k][dim] !=0) // not in original lattice
3926
                    not_in=true;
3927
                    break;
3928
                }
3929
            if(not_in)
3930
                continue;
3931
            if(inhomogeneous){
3932
                ModuleGenerators.append(rr);
3933
            }
3934
            else
3935
                Deg1Elements.append(rr);        
3936
        }
3937
    }
3938

3939
    setWeights();
3940
    if(inhomogeneous)
3941
         ModuleGenerators.sort_by_weights(WeightsGrad,GradAbs);
3942
    else
3943
        Deg1Elements.sort_by_weights(WeightsGrad,GradAbs);
3944

3945
    if(inhomogeneous){
3946
        is_Computed.set(ConeProperty::HilbertBasis);
3947
        is_Computed.set(ConeProperty::ModuleGenerators);
3948
        module_rank= ModuleGenerators.nr_of_rows();
3949
        is_Computed.set(ConeProperty::ModuleRank);
3950
        recession_rank=0;
3951
        is_Computed.set(ConeProperty::RecessionRank);
3952
        if(isComputed(ConeProperty::Grading) && module_rank>0){
3953
            multiplicity=module_rank; // of the recession cone;
3954
            is_Computed.set(ConeProperty::Multiplicity);
3955
            if(ToCompute.test(ConeProperty::HilbertSeries)){
3956
                vector<num_t> hv(1);
3957
                long raw_shift=convertTo<long>(v_scalar_product(Grading,ModuleGenerators[0]));
3958
                for(size_t i=0;i<ModuleGenerators.nr_of_rows();++i){
3959
                    long deg = convertTo<long>(v_scalar_product(Grading,ModuleGenerators[i]));
3960
                    raw_shift=min(raw_shift,deg);                        
3961
                }
3962
                for(size_t i=0;i<ModuleGenerators.nr_of_rows();++i){
3963
                    size_t deg = convertTo<long>(v_scalar_product(Grading,ModuleGenerators[i]))-raw_shift;
3964
                    if(deg+1>hv.size())
3965
                        hv.resize(deg+1);
3966
                    hv[deg]++;                        
3967
                }    
3968
                HSeries.add(hv,vector<denom_t>());
3969
                HSeries.setShift(raw_shift);
3970
                HSeries.adjustShift();
3971
                HSeries.simplify();
3972
                is_Computed.set(ConeProperty::HilbertSeries);
3973
            }
3974
        }  
3975

3976
    }
3977
    else
3978
        is_Computed.set(ConeProperty::Deg1Elements);
3979
    
3980
    is_Computed.set(ConeProperty::Approximate);
3981
    
3982
    return;    
3983
}
3984

3985
//---------------------------------------------------------------------------
3986
template<typename Integer>
3987
void Cone<Integer>::project_and_lift(Matrix<Integer>& Deg1, const Matrix<Integer>& Gens, Matrix<Integer>& Supps, bool float_projection){
3988
    
3989
    // if(verbose)
3990
    //    verboseOutput() << "Starting projection" << endl;
3991
    
3992
    // vector<boost::dynamic_bitset<> > Pair;
3993
   //  vector<boost::dynamic_bitset<> > ParaInPair;
3994
    bool is_parallelotope=(Pair.size()>0);
3995
    
3996
    vector< boost::dynamic_bitset<> > Ind;
3997

3998
    if(!is_parallelotope){
3999
        Ind=vector< boost::dynamic_bitset<> > (Supps.nr_of_rows(), boost::dynamic_bitset<> (Gens.nr_of_rows()));
4000
        for(size_t i=0;i<Supps.nr_of_rows();++i)
4001
            for(size_t j=0;j<Gens.nr_of_rows();++j)
4002
                if(v_scalar_product(Supps[i],Gens[j])==0)
4003
                    Ind[i][j]=true;
4004
    }
4005
        
4006
    size_t rank=BasisChangePointed.getRank();
4007
    
4008
    if(float_projection){
4009
        // Matrix<nmz_float> GensFloat;
4010
        // convert(GensFloat,Gens);
4011
        Matrix<nmz_float> SuppsFloat;
4012
        convert(SuppsFloat,Supps);
4013
        vector<Integer> Dummy;
4014
        // project_and_lift_inner<nmz_float,Integer>(Deg1, SuppsFloat,Ind, GradingDenom,rank,verbose,true,Dummy);
4015
        ProjectAndLift<nmz_float,Integer> PL;
4016
        if(!is_parallelotope)
4017
            PL=ProjectAndLift<nmz_float,Integer>(SuppsFloat,Ind,rank);
4018
        else
4019
            PL=ProjectAndLift<nmz_float,Integer>(SuppsFloat,Pair,ParaInPair,rank);
4020
        PL.set_grading_denom(GradingDenom);
4021
        PL.set_verbose(verbose);
4022
        PL.compute();
4023
        PL.put_eg1Points_into(Deg1);
4024
    }
4025
    else{
4026
        if (change_integer_type) {
4027
            Matrix<MachineInteger> Deg1MI(0,Deg1.nr_of_columns());
4028
            // Matrix<MachineInteger> GensMI;
4029
            Matrix<MachineInteger> SuppsMI;
4030
            try {
4031
                // convert(GensMI,Gens);
4032
                convert(SuppsMI,Supps);
4033
                MachineInteger GDMI=convertTo<MachineInteger>(GradingDenom);
4034
                vector<MachineInteger> Dummy;
4035
                // project_and_lift_inner<MachineInteger>(Deg1MI,SuppsMI,Ind, GDMI,rank,verbose,true,Dummy);
4036
                ProjectAndLift<MachineInteger,MachineInteger> PL;
4037
                if(!is_parallelotope)
4038
                    PL=ProjectAndLift<MachineInteger,MachineInteger>(SuppsMI,Ind,rank);
4039
                else
4040
                    PL=ProjectAndLift<MachineInteger,MachineInteger>(SuppsMI,Pair,ParaInPair,rank);
4041
                PL.set_grading_denom(GDMI);
4042
                PL.set_verbose(verbose);
4043
                PL.compute();
4044
                PL.put_eg1Points_into(Deg1MI);
4045
            } catch(const ArithmeticException& e) {
4046
                if (verbose) {
4047
                    verboseOutput() << e.what() << endl;
4048
                    verboseOutput() << "Restarting with a bigger type." << endl;
4049
                }
4050
                change_integer_type = false;
4051
            }
4052
            if(change_integer_type){
4053
                convert(Deg1,Deg1MI);                
4054
            }
4055
        }
4056
        
4057
        if (!change_integer_type) {
4058
            vector<Integer> Dummy;
4059
            // project_and_lift_inner<Integer>(Deg1,Supps,Ind,GradingDenom,rank,verbose,true,Dummy);
4060
            ProjectAndLift<Integer,Integer> PL;
4061
            if(!is_parallelotope)
4062
                PL=ProjectAndLift<Integer,Integer>(Supps,Ind,rank);
4063
            else
4064
                PL=ProjectAndLift<Integer,Integer>(Supps,Pair,ParaInPair,rank);
4065
            PL.set_grading_denom(GradingDenom);
4066
            PL.set_verbose(verbose);
4067
            PL.compute();
4068
            PL.put_eg1Points_into(Deg1);
4069
        }        
4070
    }
4071

4072
    is_Computed.set(ConeProperty::Projection);
4073
    if(float_projection)
4074
        is_Computed.set(ConeProperty::ProjectionFloat);
4075
    
4076
    if(verbose)
4077
        verboseOutput() << "Project-and-lift complete" << endl <<
4078
           "------------------------------------------------------------" << endl;
4079
        
4080
}
4081

4082
//---------------------------------------------------------------------------
4083
template<typename Integer>
4084
bool Cone<Integer>::check_parallelotope(){
4085
    
4086
    vector<Integer> Grad; // copy of Grading or Dehomogenization
4087

4088
    if(inhomogeneous){
4089
        Grad=Dehomogenization;
4090
    }
4091
    else{
4092
        if(!isComputed(ConeProperty::Grading))
4093
            return false;
4094
        Grad=Grading;
4095
    }
4096

4097
    Grading_Is_Coordinate=false;
4098
    size_t nr_nonzero=0;
4099
    for(size_t i=0;i<Grad.size();++i){
4100
        if(Grad[i]!=0){
4101
            nr_nonzero++;
4102
            GradingCoordinate=i;
4103
        }
4104
    }
4105
    if(nr_nonzero==1){
4106
        if(Grad[GradingCoordinate]==1)
4107
            Grading_Is_Coordinate=true;        
4108
    }
4109
    if(!Grading_Is_Coordinate)
4110
        return false;
4111
    
4112
    Matrix<Integer> Supps(SupportHyperplanes);
4113
    if(inhomogeneous)
4114
        Supps.remove_row(Grad);        
4115

4116
    size_t dim=Supps.nr_of_columns()-1; //affine dimension
4117
    if(Supps.nr_of_rows()!=2*dim)
4118
        return false;
4119
    Pair.resize(2*dim);
4120
    ParaInPair.resize(2*dim);
4121
    for(size_t i=0;i<2*dim;++i){
4122
        Pair[i].resize(dim);
4123
        Pair[i].reset();
4124
        ParaInPair[i].resize(dim);
4125
        ParaInPair[i].reset();
4126
    }
4127

4128
    vector<bool> done(2*dim);
4129
    Matrix<Integer> M2(2,dim+1), M3(3,dim+1);
4130
    M3[2]=Grad;
4131
    size_t pair_counter=0;
4132
    
4133
    vector<key_t> Supp_1; // to find antipodal vertices
4134
    vector<key_t> Supp_2;
4135
    
4136
    for(size_t i=0;i<2*dim;++i){
4137
        if(done[i])
4138
            continue;
4139
        bool parallel_found=false;
4140
        M2[0]=Supps[i];
4141
        M3[0]=Supps[i];
4142
        size_t j=i+1;
4143
        for(;j<2*dim;++j){
4144
            if(done[j]) continue;
4145
            M2[1]=Supps[j];
4146
            if(M2.rank()<2)
4147
                continue;
4148
            M3[1]=Supps[j];
4149
            if(M3.rank()==3)
4150
                continue;
4151
            else{
4152
                parallel_found=true;
4153
                done[j]=true;
4154
                break;
4155
            }
4156
        }
4157
        if(!parallel_found)
4158
            return false;
4159
        Supp_1.push_back(i);
4160
        Supp_2.push_back(j);
4161
        Pair[i][pair_counter]=true;
4162
        Pair[j][pair_counter]=true;
4163
        ParaInPair[j][pair_counter]=true;
4164
        pair_counter++;
4165
    }
4166
    
4167
    Matrix<Integer> v1=Supps.submatrix(Supp_1).kernel(); // opposite vertices
4168
    Matrix<Integer> v2=Supps.submatrix(Supp_2).kernel();
4169
    Integer MinusOne=-1;
4170
    if(v_scalar_product(v1[0],Grad)==0)
4171
        return false;
4172
    if(v_scalar_product(v2[0],Grad)==0)
4173
        return false;
4174
    if(v_scalar_product(v1[0],Grad)<0)
4175
        v_scalar_multiplication(v1[0],MinusOne);
4176
    if(v_scalar_product(v2[0],Grad)<0)
4177
        v_scalar_multiplication(v2[0],MinusOne);
4178
    
4179
    /* cout << Supp_1;
4180
    cout << Supp_2;
4181
    v1.pretty_print(cout);
4182
    v2.pretty_print(cout);
4183
    cout << "==============" << endl;
4184
    Supps.pretty_print(cout);
4185
    cout << "==============" << endl;*/
4186
    if(v1.nr_of_rows()!=1 || v2.nr_of_rows()!=1)
4187
        return false;
4188
    for(size_t i=0;i<Supp_1.size();++i){
4189
        // cout << "i " << i << " " << v_scalar_product(Supps[Supp_1[i]],v2[0]) << endl;
4190
        if(!v_scalar_product_positive(Supps[Supp_1[i]],v2[0]))
4191
            return false;
4192
    }
4193
    for(size_t i=0;i<Supp_2.size();++i){
4194
        // cout << "i " << i << " " << v_scalar_product(Supps[Supp_2[i]],v1[0]) << endl;
4195
        if(!v_scalar_product_positive(Supps[Supp_2[i]],v1[0]))
4196
            return false;
4197
    }
4198
    
4199
    // we have found opposite vertices
4200
    
4201
    return true;    
4202
}
4203

4204

4205

4206
} // end namespace libnormaliz
4207

4208
Product

Resources

Company
