1 /* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2011 Bjarke Hammersholt Roune (www.broune.com)
3
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 2 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 #ifndef LATTICE_ALGS_GUARD
18 #define LATTICE_ALGS_GUARD
19
20 #include "SatBinomIdeal.h"
21 #include "IOFacade.h"
22 #include "Scanner.h"
23 #include "IOHandler.h"
24 #include "DataType.h"
25 #include "BigIdeal.h"
26 #include "MsmStrategy.h"
27 #include "TermTranslator.h"
28 #include "TranslatingTermConsumer.h"
29 #include "DebugStrategy.h"
30 #include "Matrix.h"
31 #include "BigTermRecorder.h"
32 #include "SliceParams.h"
33 #include "SliceFacade.h"
34
35 #include <algorithm>
36 #include <set>
37 #include <sstream>
38 #include <limits>
39 #include <fstream>
40 #include <map>
41
42 // Computations and data structures on lattices.
43 // Support for LatticeAnalyzeAction.
44
45 #include <iostream>
46
47 // wrapped in do .. while(false) to make it act like a single statement.
48 #define CHECK(X) \
49 do { \
50 if (!(X)) { \
51 cout << "Check condition on line " \
52 << __LINE__ << " of file " << __FILE__ \
53 << " not satisfied:\n "#X << endl; \
54 exit(1); \
55 } \
56 } while (false)
57
58 enum NeighborPlace {
59 InPlane,
60 UnderPlane,
61 OverPlane,
62 NoPlace
63 };
64
65 char getPlaceCode(NeighborPlace place);
66
67 class Mlfb;
68
69 struct SeqPos {
70 SeqPos();
71 SeqPos(const Mlfb* mlfb, size_t nextFacet, size_t previousFacet);
72 size_t getForwardFacet() const;
73 size_t getBackFacet() const;
74 SeqPos getReverse() const;
75 void order();
76 bool operator<(const SeqPos& pos) const;
77
78 const Mlfb* mlfb;
79 size_t fixFacet1;
80 size_t fixFacet2;
81 size_t comingFromFacet;
82 };
83
84 class GrobLat;
85
86 class Neighbor {
87 public:
88 Neighbor(); // invalid
89 Neighbor(const GrobLat& lat); // zero
90 Neighbor(const GrobLat& lat, const size_t row); // row of lat
91
92 Neighbor& operator=(const Neighbor& neighbor) {
93 _lat = neighbor._lat;
94 _row = neighbor._row;
95 return *this;
96 }
97
98 bool operator==(const Neighbor& neighbor) const {
99 return _lat == neighbor._lat && getRow() == neighbor.getRow();
100 }
101
102 const mpq_class& getH(size_t i) const;
103 size_t getHDim() const;
104
105 const mpq_class& getY(size_t i) const;
106 size_t getYDim() const;
107
getRow()108 size_t getRow() const {return _row;}
109
110 bool isZero() const;
111 bool isValid() const;
112 bool isSpecial() const;
113 bool isGenerator() const;
114
115 string getName() const;
getGrobLat()116 const GrobLat& getGrobLat() const {return *_lat;}
117
118 private:
119 const GrobLat* _lat;
120 size_t _row;
121 };
122
123 /** A lattice with associated Grobner basis/neighbors. */
124 class GrobLat {
125 public:
126 GrobLat(const Matrix& matrix, const SatBinomIdeal& ideal);
127
getNeighbor(size_t row)128 Neighbor getNeighbor(size_t row) const {
129 ASSERT(row < getNeighborCount());
130 return Neighbor(*this, row);
131 }
132
getNeighborCount()133 size_t getNeighborCount() const {
134 ASSERT(_h.getRowCount() == _y.getRowCount());
135 return _y.getRowCount();
136 }
137
getYMatrix()138 const Matrix& getYMatrix() const {return _y;}
getHMatrix()139 const Matrix& getHMatrix() const {return _h;}
getMatrix()140 const Matrix& getMatrix() const {return _mat;}
141
getIdeal()142 const SatBinomIdeal& getIdeal() const {return _ideal;}
143
getSum(Neighbor a,Neighbor b)144 Neighbor getSum(Neighbor a, Neighbor b) const {
145 vector<mpq_class> sum(getHDim());
146 for (size_t i = 0; i < getHDim(); ++i)
147 sum[i] = _h(a.getRow(), i) + _h(b.getRow(), i);
148 for (size_t row = 0; row < _h.getRowCount(); ++row) {
149 bool match = true;
150 for (size_t col = 0; col < _h.getColCount(); ++col)
151 if (sum[col] != _h(row, col))
152 match = false;
153 if (match)
154 return Neighbor(*this, row);
155 }
156 return Neighbor();
157 }
158
getSum(size_t a,size_t b)159 Neighbor getSum(size_t a, size_t b) const {
160 vector<mpq_class> sum(getHDim());
161 for (size_t i = 0; i < getHDim(); ++i)
162 sum[i] = _h(a, i) + _h(b, i);
163 for (size_t row = 0; row < _h.getRowCount(); ++row) {
164 bool match = true;
165 for (size_t col = 0; col < _h.getColCount(); ++col)
166 if (sum[col] != _h(row, col))
167 match = false;
168 if (match)
169 return Neighbor(*this, row);
170 }
171 return Neighbor();
172 }
173
getYDim()174 size_t getYDim() const {
175 ASSERT(_y.getColCount() == _ideal.getVarCount());
176 return _y.getColCount();
177 }
178
getHDim()179 size_t getHDim() const {
180 return _h.getColCount();
181 }
182
hasZeroEntryY()183 bool hasZeroEntryY() const {
184 return _ideal.hasZeroEntry();
185 }
186
getInitialIdeal(BigIdeal & ideal)187 void getInitialIdeal(BigIdeal& ideal) const {
188 _ideal.getInitialIdeal(ideal);
189 }
190
isSum(Neighbor n)191 bool isSum(Neighbor n) const {
192 ASSERT(n.isValid());
193 ASSERT(&n.getGrobLat() == this);
194 ASSERT(n.getRow() < _isSumRow.size());
195 return _isSumRow[n.getRow()];
196 }
getNonSums()197 const vector<Neighbor>& getNonSums() const {return _nonSums;}
getZero()198 const mpq_class& getZero() const {return _zero;} // todo: remove
199
200 /// Returns true if the smallest body containing zero, a and b has
201 /// no neighbor in its interior.
isPointFreeBody(Neighbor a,Neighbor b)202 bool isPointFreeBody(Neighbor a, Neighbor b) const {
203 return _ideal.isPointFreeBody(_ideal.getGenerator(a.getRow()),
204 _ideal.getGenerator(b.getRow()));
205 }
206
207 /// Returns true if the smallest body containing zero, a, b and c
208 /// has no neighbor in its interior.
isPointFreeBody(Neighbor a,Neighbor b,Neighbor c)209 bool isPointFreeBody(Neighbor a, Neighbor b, Neighbor c) const {
210 return _ideal.isPointFreeBody(_ideal.getGenerator(a.getRow()),
211 _ideal.getGenerator(b.getRow()),
212 _ideal.getGenerator(c.getRow()));
213 }
214
isInterior(Neighbor a,Neighbor b)215 bool isInterior(Neighbor a, Neighbor b) const {
216 if (!isPointFreeBody(a, b))
217 return false;
218 for (size_t var = 1; var < a.getYDim(); ++var)
219 if (a.getY(var) <= 0 && b.getY(var) <= 0)
220 return false;
221 return true;
222 }
223
224 private:
225 vector<bool> _isSumRow;
226 vector<Neighbor> _nonSums;
227
228 Matrix _y; // rows are neighbors in y-space
229 Matrix _h; // rows are neighbors in h-space
230 Matrix _mat; // matrix that defines lattice
231 SatBinomIdeal _ideal; // other representation of _y, necessary for now
232 mpq_class _zero;
233 };
234
235 class Tri {
236 public:
237 Tri(Neighbor a, Neighbor b, Neighbor sum,
238 const vector<Mlfb>& mlfbs, const GrobLat& lat);
239
getA()240 Neighbor getA() const {return _a;}
getB()241 Neighbor getB() const {return _b;}
getSum()242 Neighbor getSum() const {return _sum;}
getASideMlfbs()243 const vector<const Mlfb*>& getASideMlfbs() const {return _aSideMlfbs;}
getBSideMlfbs()244 const vector<const Mlfb*>& getBSideMlfbs() const {return _bSideMlfbs;}
getNeighborsOnBoundary()245 const vector<Neighbor>& getNeighborsOnBoundary() const {return _boundary;}
getNeighborsInInterior()246 const vector<Neighbor>& getNeighborsInInterior() const {return _interior;}
247
248 private:
249 Neighbor _a;
250 Neighbor _b;
251 Neighbor _sum; // neighbor that is sum of a and b
252 vector<const Mlfb*> _aSideMlfbs; // MLFBs containing {0,a,sum}
253 vector<const Mlfb*> _bSideMlfbs; // MLFBs containing {0,b,sum}
254 vector<Neighbor> _interior; // neighbors on boundary of <0,a,b,sum>
255 vector<Neighbor> _boundary; // neighbors in interior of <0,a,b,sum>
256 };
257
258 class Plane {
259 public:
Plane(Neighbor a,Neighbor b,Neighbor sum,const vector<Mlfb> & mlfbs,const GrobLat & lat)260 Plane(Neighbor a, Neighbor b, Neighbor sum,
261 const vector<Mlfb>& mlfbs, const GrobLat& lat):
262 tri(a, b, sum, mlfbs, lat) {}
263
264 size_t getTypeCount(size_t type) const;
265 size_t getMaxType() const;
266 NeighborPlace getPlace(Neighbor neighbor) const;
267 bool inPlane(Neighbor neighbor) const;
268 bool isPivot(const Mlfb& mlfb) const;
269 bool isSidePivot(const Mlfb& mlfb) const;
270 bool isFlat(const Mlfb& mlfb) const;
271 bool is22(const Mlfb& mlfb) const;
272 size_t getType(const Mlfb& mlfb) const;
273
hasFlat()274 bool hasFlat() const {
275 return getTypeCount(4) > 0;
276 }
277
278 Matrix nullSpaceBasis;
279 Tri tri;
280 Matrix rowAB;
281 size_t flatIntervalCount;
282
283 map<size_t, size_t> typeCounts;
284 vector<NeighborPlace> neighborPlace;
285 vector<SeqPos> flatSeq;
286 vector<const Mlfb*> pivots;
287 };
288
289 class Mlfb {
290 public:
getMinInitialFacet()291 size_t getMinInitialFacet() const {
292 return minInitialFacet;
293 }
294
hasPoint(Neighbor n)295 bool hasPoint(Neighbor n) const {
296 for (size_t i = 0; i < _points.size(); ++i)
297 if (_points[i] == n)
298 return true;
299 return false;
300 }
301
getPoint(size_t offset)302 Neighbor getPoint(size_t offset) const {
303 ASSERT(offset < getPointCount());
304 return _points[offset];
305 }
306
getPointCount()307 size_t getPointCount() const {
308 return _points.size();
309 }
310
311 bool operator==(const Mlfb& mlfb) const {
312 return _offset == mlfb._offset;
313 }
314
getOffset()315 size_t getOffset() const {
316 return _offset;
317 }
318
getRhs()319 const vector<mpz_class>& getRhs() const {
320 return _rhs;
321 }
322
getName()323 string getName() const {
324 ostringstream name;
325 name << 'm' << (getOffset() + 1);
326 return name.str();
327 }
328
getName(const Plane & plane)329 string getName(const Plane& plane) const {
330 if (plane.isPivot(*this))
331 return getName() + 'P';
332 if (plane.isFlat(*this))
333 return getName() + 'F';
334 return getName();
335 }
336
getHitsNeighbor(size_t indexParam)337 Neighbor getHitsNeighbor(size_t indexParam) const {
338 ASSERT(indexParam < edgeHitsFacet.size());
339 return getPoint(getHitsFacet(indexParam));
340 }
341
getHitsFacet(size_t indexParam)342 size_t getHitsFacet(size_t indexParam) const {
343 ASSERT(indexParam < edgeHitsFacet.size());
344 return edgeHitsFacet[indexParam];
345 }
346
getEdge(size_t indexParam)347 const Mlfb* getEdge(size_t indexParam) const {
348 ASSERT(indexParam < edges.size());
349 return edges[indexParam];
350 }
351
getEdge(size_t indexParam)352 Mlfb* getEdge(size_t indexParam) {
353 ASSERT(indexParam < edges.size());
354 return edges[indexParam];
355 }
356
getFacetOf(const Mlfb & adjacent)357 size_t getFacetOf(const Mlfb& adjacent) const {
358 for (size_t i = 0; i < 4; ++i)
359 if (*getEdge(i) == adjacent)
360 return i;
361 return 4;
362 }
363
isParallelogram()364 bool isParallelogram() const {
365 return _isParallelogram;
366 }
367
368 mpq_class index;
369 mpz_class dotDegree;
370 vector<Mlfb*> edges;
371 vector<size_t> edgeHitsFacet;
372 size_t minInitialFacet;
373
374 void reset(size_t offset, const vector<Neighbor>& points);
375
376 private:
377 vector<mpz_class> _rhs;
378 vector<Neighbor> _points;
379 size_t _offset;
380 bool _isParallelogram;
381 };
382
383 class TriPlane {
384 public:
TriPlane(Neighbor a,Neighbor b,Neighbor c)385 TriPlane(Neighbor a, Neighbor b, Neighbor c):
386 _a(a), _b(b), _c(c) {
387
388 Matrix mat(2, 3);
389 for (size_t col = 0; col < 3; ++col) {
390 mat(0, col) = a.getH(col) - c.getH(col);
391 mat(1, col) = b.getH(col) - c.getH(col);
392 }
393
394 nullSpace(_normal, mat);
395 transpose(_normal, _normal);
396 _line = (_normal.getRowCount() != 1);
397 }
398
isLine()399 bool isLine() const {
400 return _line;
401 }
402
closeToPlane(Neighbor a)403 bool closeToPlane(Neighbor a) {
404 ASSERT(!isLine());
405 mpz_class dn = dotNormal(a);
406 return dn == 0 || dn == 1 || dn == -1;
407 }
408
inPlane(Neighbor a)409 bool inPlane(Neighbor a) const {
410 return dotNormal(a) == 0;
411 }
412
dotNormal(Neighbor a)413 mpz_class dotNormal(Neighbor a) const {
414 mpz_class prod = 0;
415 for (size_t i = 0; i < 3; ++i)
416 prod += a.getH(i) * _normal(0, i);
417 return prod;
418 }
419
isParallel(const TriPlane & plane)420 bool isParallel(const TriPlane& plane) const {
421 ASSERT(!isLine());
422 ASSERT(!plane.isLine());
423 mpz_class da = plane.dotNormal(_a);
424 return plane.dotNormal(_b) == da && plane.dotNormal(_c) == da;
425 }
426
isParallel(const Plane & plane)427 bool isParallel(const Plane& plane) const {
428 ASSERT(!isLine());
429 return plane.nullSpaceBasis == getNormal();
430 }
431
432 /// returns the normal of the plane as the row of a matrix.
getNormal()433 const Matrix& getNormal() const {
434 return _normal;
435 }
436
437 private:
438 Neighbor _a, _b, _c;
439 Matrix _normal;
440 bool _line;
441 };
442
443 void getThinPlanes(vector<TriPlane>& planes, const GrobLat& lat);
444 void checkPlanes(const vector<TriPlane>& thinPlanes,
445 const vector<Plane>& dtPlanes);
446
447 size_t pushOutFacetPositive(size_t facetPushOut,
448 const vector<mpz_class>& rhs,
449 const GrobLat& lat);
450
451 size_t pushOutFacetZero(const vector<mpz_class>& rhs, const GrobLat& lat);
452
453 void computeMlfbs(vector<Mlfb>& mlfbs, const GrobLat& lat);
454
455
456 void computeSeqs(vector<vector<SeqPos> >& left,
457 vector<vector<SeqPos> >& right,
458 const vector<Mlfb>& mlfbs,
459 const Plane& plane);
460
461 /** Starting at pivot (which must be a pivot), follow the three
462 non-flat sequences starting at pivot. */
463 void computePivotSeqs(vector<vector<SeqPos> >& seqs, const Mlfb& pivot,
464 const Plane& plane);
465
466 void checkSeqs(const vector<vector<SeqPos> >& left,
467 const vector<vector<SeqPos> >& right,
468 const Plane& plane,
469 const vector<Mlfb>& mlfbs);
470 void checkMiddle(const Plane& plane,
471 const vector<Mlfb>& mlfbs);
472 void checkDoubleTriangle(const Plane& plane,
473 const vector<Mlfb>& mlfbs);
474 void checkGraph(const vector<Mlfb>& mlfbs);
475 void checkGraphOnPlane(const Plane& plane,
476 const vector<Mlfb>& mlfbs);
477
478
479 /** Perform checks where pivotSeqs are the 3 non-flat sequences on one
480 side. */
481 void checkPivotSeqs(vector<vector<SeqPos> >& pivotSeqs,
482 const Plane& plane,
483 const vector<Mlfb>& mlfbs,
484 const vector<SeqPos>& flatSeq);
485
486 void checkPlaneTri(const GrobLat& lat,
487 const vector<Mlfb>& mlfbs,
488 const vector<const Mlfb*>& pivots,
489 const Plane& plane);
490
491 void computePlanes(vector<Plane>& planes,
492 const GrobLat& lat,
493 vector<Mlfb>& mlfbs);
494
495 /** Set the plane vertex count for each mlfb and count how many MLFBs
496 have each possible number of vertices. */
497 void setupPlaneCountsAndOrder(vector<Mlfb>& mlfbs,
498 const Plane& plane,
499 map<size_t, size_t>& typeCounts);
500
501 bool disjointSeqs(const vector<SeqPos>& a, const vector<SeqPos>& b);
502
503 /** Put all pivots into pivots. flatSeq must be the sequence of
504 flats. If flatSeq is not empty, then offsets 0,1 will be the left
505 pivots while 2,3 will be the right pivots. */
506 void computePivots(vector<const Mlfb*>& pivots,
507 const vector<Mlfb>& mlfbs,
508 const Plane& plane,
509 const vector<SeqPos>& flatSeq);
510
511 void checkNonSums(const GrobLat& lat);
512
513 void checkFlatSeq(const vector<SeqPos>& flatSeq,
514 const GrobLat& lat,
515 const Plane& plane);
516
517 const char* getEdgePos(size_t index);
518 mpq_class getIndexSum(const vector<Mlfb>& mlfbs);
519
520 void checkMlfbs(const vector<Mlfb>& mlfbs, const GrobLat& lat);
521 void checkDoubleTrianglePlanes(const vector<Plane>& planes,
522 const GrobLat& lat,
523 const vector<Mlfb>& mlfbs);
524 void checkPlane(const Plane& plane, const vector<Mlfb>& mlfbs);
525
526 #endif
527