Frobby 0.9.5
LatticeAlgs.h
Go to the documentation of this file.
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"
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
64
66
67class Mlfb;
68
69struct 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;
82};
83
84class GrobLat;
85
86class 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
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
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;
116 const GrobLat& getGrobLat() const {return *_lat;}
117
118 private:
119 const GrobLat* _lat;
120 size_t _row;
121};
122
124class GrobLat {
125 public:
126 GrobLat(const Matrix& matrix, const SatBinomIdeal& ideal);
127
128 Neighbor getNeighbor(size_t row) const {
130 return Neighbor(*this, row);
131 }
132
133 size_t getNeighborCount() const {
135 return _y.getRowCount();
136 }
137
138 const Matrix& getYMatrix() const {return _y;}
139 const Matrix& getHMatrix() const {return _h;}
140 const Matrix& getMatrix() const {return _mat;}
141
142 const SatBinomIdeal& getIdeal() const {return _ideal;}
143
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
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
174 size_t getYDim() const {
176 return _y.getColCount();
177 }
178
179 size_t getHDim() const {
180 return _h.getColCount();
181 }
182
183 bool hasZeroEntryY() const {
184 return _ideal.hasZeroEntry();
185 }
186
187 void getInitialIdeal(BigIdeal& ideal) const {
188 _ideal.getInitialIdeal(ideal);
189 }
190
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 }
197 const vector<Neighbor>& getNonSums() const {return _nonSums;}
198 const mpq_class& getZero() const {return _zero;} // todo: remove
199
203 return _ideal.isPointFreeBody(_ideal.getGenerator(a.getRow()),
204 _ideal.getGenerator(b.getRow()));
205 }
206
210 return _ideal.isPointFreeBody(_ideal.getGenerator(a.getRow()),
211 _ideal.getGenerator(b.getRow()),
212 _ideal.getGenerator(c.getRow()));
213 }
214
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
233};
234
235class Tri {
236 public:
238 const vector<Mlfb>& mlfbs, const GrobLat& lat);
239
240 Neighbor getA() const {return _a;}
241 Neighbor getB() const {return _b;}
242 Neighbor getSum() const {return _sum;}
243 const vector<const Mlfb*>& getASideMlfbs() const {return _aSideMlfbs;}
244 const vector<const Mlfb*>& getBSideMlfbs() const {return _bSideMlfbs;}
245 const vector<Neighbor>& getNeighborsOnBoundary() const {return _boundary;}
246 const vector<Neighbor>& getNeighborsInInterior() const {return _interior;}
247
248 private:
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
258class Plane {
259 public:
263
264 size_t getTypeCount(size_t type) const;
265 size_t getMaxType() 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
274 bool hasFlat() const {
275 return getTypeCount(4) > 0;
276 }
277
282
284 vector<NeighborPlace> neighborPlace;
285 vector<SeqPos> flatSeq;
286 vector<const Mlfb*> pivots;
287};
288
289class Mlfb {
290 public:
291 size_t getMinInitialFacet() const {
292 return minInitialFacet;
293 }
294
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
302 Neighbor getPoint(size_t offset) const {
304 return _points[offset];
305 }
306
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
315 size_t getOffset() const {
316 return _offset;
317 }
318
319 const vector<mpz_class>& getRhs() const {
320 return _rhs;
321 }
322
323 string getName() const {
325 name << 'm' << (getOffset() + 1);
326 return name.str();
327 }
328
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
341
342 size_t getHitsFacet(size_t indexParam) const {
345 }
346
347 const Mlfb* getEdge(size_t indexParam) const {
348 ASSERT(indexParam < edges.size());
349 return edges[indexParam];
350 }
351
353 ASSERT(indexParam < edges.size());
354 return edges[indexParam];
355 }
356
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
364 bool isParallelogram() const {
365 return _isParallelogram;
366 }
367
370 vector<Mlfb*> edges;
371 vector<size_t> edgeHitsFacet;
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;
381};
382
383class TriPlane {
384public:
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
396 _line = (_normal.getRowCount() != 1);
397 }
398
399 bool isLine() const {
400 return _line;
401 }
402
404 ASSERT(!isLine());
406 return dn == 0 || dn == 1 || dn == -1;
407 }
408
409 bool inPlane(Neighbor a) const {
410 return dotNormal(a) == 0;
411 }
412
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
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
427 bool isParallel(const Plane& plane) const {
428 ASSERT(!isLine());
429 return plane.nullSpaceBasis == getNormal();
430 }
431
433 const Matrix& getNormal() const {
434 return _normal;
435 }
436
437private:
440 bool _line;
441};
442
445 const vector<Plane>& dtPlanes);
446
448 const vector<mpz_class>& rhs,
449 const GrobLat& lat);
450
451size_t pushOutFacetZero(const vector<mpz_class>& rhs, const GrobLat& lat);
452
454
455
456void computeSeqs(vector<vector<SeqPos> >& left,
457 vector<vector<SeqPos> >& right,
458 const vector<Mlfb>& mlfbs,
459 const Plane& plane);
460
463void computePivotSeqs(vector<vector<SeqPos> >& seqs, const Mlfb& pivot,
464 const Plane& plane);
465
466void checkSeqs(const vector<vector<SeqPos> >& left,
467 const vector<vector<SeqPos> >& right,
468 const Plane& plane,
469 const vector<Mlfb>& mlfbs);
470void checkMiddle(const Plane& plane,
471 const vector<Mlfb>& mlfbs);
473 const vector<Mlfb>& mlfbs);
474void checkGraph(const vector<Mlfb>& mlfbs);
475void checkGraphOnPlane(const Plane& plane,
476 const vector<Mlfb>& mlfbs);
477
478
481void checkPivotSeqs(vector<vector<SeqPos> >& pivotSeqs,
482 const Plane& plane,
483 const vector<Mlfb>& mlfbs,
484 const vector<SeqPos>& flatSeq);
485
486void checkPlaneTri(const GrobLat& lat,
487 const vector<Mlfb>& mlfbs,
488 const vector<const Mlfb*>& pivots,
489 const Plane& plane);
490
492 const GrobLat& lat,
494
498 const Plane& plane,
499 map<size_t, size_t>& typeCounts);
500
501bool disjointSeqs(const vector<SeqPos>& a, const vector<SeqPos>& b);
502
506void computePivots(vector<const Mlfb*>& pivots,
507 const vector<Mlfb>& mlfbs,
508 const Plane& plane,
509 const vector<SeqPos>& flatSeq);
510
511void checkNonSums(const GrobLat& lat);
512
513void checkFlatSeq(const vector<SeqPos>& flatSeq,
514 const GrobLat& lat,
515 const Plane& plane);
516
517const char* getEdgePos(size_t index);
519
520void checkMlfbs(const vector<Mlfb>& mlfbs, const GrobLat& lat);
522 const GrobLat& lat,
523 const vector<Mlfb>& mlfbs);
524void checkPlane(const Plane& plane, const vector<Mlfb>& mlfbs);
525
526#endif
void checkFlatSeq(const vector< SeqPos > &flatSeq, const GrobLat &lat, const Plane &plane)
void checkGraphOnPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
NeighborPlace
Definition LatticeAlgs.h:58
@ InPlane
Definition LatticeAlgs.h:59
@ OverPlane
Definition LatticeAlgs.h:61
@ UnderPlane
Definition LatticeAlgs.h:60
@ NoPlace
Definition LatticeAlgs.h:62
void checkPlanes(const vector< TriPlane > &thinPlanes, const vector< Plane > &dtPlanes)
void checkMlfbs(const vector< Mlfb > &mlfbs, const GrobLat &lat)
bool disjointSeqs(const vector< SeqPos > &a, const vector< SeqPos > &b)
void computePivotSeqs(vector< vector< SeqPos > > &seqs, const Mlfb &pivot, const Plane &plane)
Starting at pivot (which must be a pivot), follow the three non-flat sequences starting at pivot.
void checkSeqs(const vector< vector< SeqPos > > &left, const vector< vector< SeqPos > > &right, const Plane &plane, const vector< Mlfb > &mlfbs)
void computeSeqs(vector< vector< SeqPos > > &left, vector< vector< SeqPos > > &right, const vector< Mlfb > &mlfbs, const Plane &plane)
size_t pushOutFacetZero(const vector< mpz_class > &rhs, const GrobLat &lat)
char getPlaceCode(NeighborPlace place)
void computeMlfbs(vector< Mlfb > &mlfbs, const GrobLat &lat)
void computePivots(vector< const Mlfb * > &pivots, const vector< Mlfb > &mlfbs, const Plane &plane, const vector< SeqPos > &flatSeq)
Put all pivots into pivots.
void checkGraph(const vector< Mlfb > &mlfbs)
void checkMiddle(const Plane &plane, const vector< Mlfb > &mlfbs)
mpq_class getIndexSum(const vector< Mlfb > &mlfbs)
void checkDoubleTriangle(const Plane &plane, const vector< Mlfb > &mlfbs)
void computePlanes(vector< Plane > &planes, const GrobLat &lat, vector< Mlfb > &mlfbs)
const char * getEdgePos(size_t index)
void checkDoubleTrianglePlanes(const vector< Plane > &planes, const GrobLat &lat, const vector< Mlfb > &mlfbs)
size_t pushOutFacetPositive(size_t facetPushOut, const vector< mpz_class > &rhs, const GrobLat &lat)
void getThinPlanes(vector< TriPlane > &planes, const GrobLat &lat)
void checkNonSums(const GrobLat &lat)
void checkPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
void setupPlaneCountsAndOrder(vector< Mlfb > &mlfbs, const Plane &plane, map< size_t, size_t > &typeCounts)
Set the plane vertex count for each mlfb and count how many MLFBs have each possible number of vertic...
void checkPlaneTri(const GrobLat &lat, const vector< Mlfb > &mlfbs, const vector< const Mlfb * > &pivots, const Plane &plane)
void checkPivotSeqs(vector< vector< SeqPos > > &pivotSeqs, const Plane &plane, const vector< Mlfb > &mlfbs, const vector< SeqPos > &flatSeq)
Perform checks where pivotSeqs are the 3 non-flat sequences on one side.
void nullSpace(Matrix &basis, const Matrix &matParam)
Sets the columns of basis to a basis of the null space of mat.
Definition Matrix.cpp:296
void transpose(Matrix &trans, const Matrix &mat)
Sets trans to the transpose of mat.
Definition Matrix.cpp:129
void nameFactoryRegister(NameFactory< AbstractProduct > &factory)
Registers the string returned by ConcreteProduct::getStaticName() to a function that default-construc...
A lattice with associated Grobner basis/neighbors.
Matrix _h
Neighbor getNeighbor(size_t row) const
size_t getYDim() const
Neighbor getSum(size_t a, size_t b) const
size_t getHDim() const
SatBinomIdeal _ideal
mpq_class _zero
vector< bool > _isSumRow
Matrix _y
const Matrix & getMatrix() const
bool hasZeroEntryY() const
Neighbor getSum(Neighbor a, Neighbor b) const
const SatBinomIdeal & getIdeal() const
bool isPointFreeBody(Neighbor a, Neighbor b, Neighbor c) const
Returns true if the smallest body containing zero, a, b and c has no neighbor in its interior.
size_t getNeighborCount() const
void getInitialIdeal(BigIdeal &ideal) const
bool isInterior(Neighbor a, Neighbor b) const
Matrix _mat
const Matrix & getYMatrix() const
const vector< Neighbor > & getNonSums() const
bool isSum(Neighbor n) const
const Matrix & getHMatrix() const
vector< Neighbor > _nonSums
const mpq_class & getZero() const
bool isPointFreeBody(Neighbor a, Neighbor b) const
Returns true if the smallest body containing zero, a and b has no neighbor in its interior.
size_t getColCount() const
Definition Matrix.h:31
size_t getRowCount() const
Definition Matrix.h:30
string getName() const
const vector< mpz_class > & getRhs() const
Neighbor getPoint(size_t offset) const
vector< mpz_class > _rhs
string getName(const Plane &plane) const
bool operator==(const Mlfb &mlfb) const
vector< Neighbor > _points
vector< Mlfb * > edges
mpq_class index
size_t getFacetOf(const Mlfb &adjacent) const
vector< size_t > edgeHitsFacet
Neighbor getHitsNeighbor(size_t indexParam) const
bool hasPoint(Neighbor n) const
size_t _offset
Mlfb * getEdge(size_t indexParam)
void reset(size_t offset, const vector< Neighbor > &points)
mpz_class dotDegree
size_t minInitialFacet
size_t getOffset() const
size_t getPointCount() const
bool isParallelogram() const
const Mlfb * getEdge(size_t indexParam) const
size_t getMinInitialFacet() const
size_t getHitsFacet(size_t indexParam) const
bool _isParallelogram
const GrobLat * _lat
bool isSpecial() const
size_t _row
bool isGenerator() const
const mpq_class & getH(size_t i) const
string getName() const
size_t getRow() const
size_t getYDim() const
const mpq_class & getY(size_t i) const
size_t getHDim() const
const GrobLat & getGrobLat() const
bool isZero() const
Neighbor & operator=(const Neighbor &neighbor)
Definition LatticeAlgs.h:92
bool operator==(const Neighbor &neighbor) const
Definition LatticeAlgs.h:98
bool isValid() const
bool isFlat(const Mlfb &mlfb) const
vector< NeighborPlace > neighborPlace
vector< SeqPos > flatSeq
Tri tri
bool isPivot(const Mlfb &mlfb) const
size_t flatIntervalCount
Matrix nullSpaceBasis
Plane(Neighbor a, Neighbor b, Neighbor sum, const vector< Mlfb > &mlfbs, const GrobLat &lat)
bool is22(const Mlfb &mlfb) const
NeighborPlace getPlace(Neighbor neighbor) const
bool hasFlat() const
map< size_t, size_t > typeCounts
vector< const Mlfb * > pivots
size_t getTypeCount(size_t type) const
bool isSidePivot(const Mlfb &mlfb) const
Matrix rowAB
size_t getMaxType() const
bool inPlane(Neighbor neighbor) const
size_t getType(const Mlfb &mlfb) const
Represents a saturated binomial ideal.
bool isPointFreeBody(const vector< mpz_class > &a, const vector< mpz_class > &b) const
Returns true if the smallest body containing zero, a and b has no generator in its interior.
size_t getVarCount() const
const vector< mpz_class > & getGenerator(size_t index) const
bool hasZeroEntry() const
Returns true if any generator does not involve every variable, i.e.
void getInitialIdeal(BigIdeal &ideal) const
mpz_class dotNormal(Neighbor a) const
bool inPlane(Neighbor a) const
bool closeToPlane(Neighbor a)
Matrix _normal
bool isParallel(const TriPlane &plane) const
const Matrix & getNormal() const
returns the normal of the plane as the row of a matrix.
Neighbor _a
Neighbor _b
TriPlane(Neighbor a, Neighbor b, Neighbor c)
bool isParallel(const Plane &plane) const
Neighbor _c
bool isLine() const
Neighbor _a
vector< const Mlfb * > _bSideMlfbs
Neighbor _b
const vector< const Mlfb * > & getASideMlfbs() const
const vector< Neighbor > & getNeighborsInInterior() const
Neighbor getB() const
vector< Neighbor > _interior
const vector< Neighbor > & getNeighborsOnBoundary() const
const vector< const Mlfb * > & getBSideMlfbs() const
Neighbor getA() const
Neighbor getSum() const
vector< Neighbor > _boundary
Neighbor _sum
vector< const Mlfb * > _aSideMlfbs
#define ASSERT(X)
Definition stdinc.h:86
void order()
SeqPos getReverse() const
const Mlfb * mlfb
Definition LatticeAlgs.h:78
size_t comingFromFacet
Definition LatticeAlgs.h:81
size_t fixFacet1
Definition LatticeAlgs.h:79
bool operator<(const SeqPos &pos) const
size_t fixFacet2
Definition LatticeAlgs.h:80
size_t getForwardFacet() const
size_t getBackFacet() const