Frobby 0.9.5
LatticeAnalyzeAction.cpp
Go to the documentation of this file.
1/* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2009 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#include "stdinc.h"
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 "ColumnPrinter.h"
32#include "BigTermRecorder.h"
33#include "SliceParams.h"
34#include "SliceFacade.h"
35
36#include <algorithm>
37#include <set>
38#include <sstream>
39#include <limits>
40#include <fstream>
41#include <map>
42
43#include "LatticeAlgs.h"
44
53#include <iostream>
54
55namespace {
56
57 void printIndentedMatrix(const Matrix& matrix) {
59 pr.setPrefix(" ");
60 print(pr, matrix);
61
62 fputc('\n', stdout);
63 print(stdout, pr);
64 fputc('\n', stdout);
65 }
66
67 void printNullSpace(const Matrix& matrix) {
68 Matrix nullSpaceBasis;
69 nullSpace(nullSpaceBasis, matrix);
70 transpose(nullSpaceBasis, nullSpaceBasis);
71
72 fputs("The right null space is spanned by the rows of\n",
73 stdout);
74 printIndentedMatrix(nullSpaceBasis);
75 }
76
77 class NeighborPrinter {
78 public:
79 NeighborPrinter(const GrobLat& lat):
80 _lat(lat) {
81
82 // "gXYZ:" label
83 _labelIndex = _pr.getColumnCount();
84 _pr.addColumn(false, " ");
85
86 // h space
87 _hHeader = _pr.getColumnCount();
88 _pr.addColumn(false, " ", "");
89 _hIndex = _pr.getColumnCount();
90 for (size_t i = 0; i < _lat.getHDim(); ++i)
91 _pr.addColumn(false, i == 0 ? " " : " ");
92
93 _comma = _pr.getColumnCount();
94 _pr.addColumn(false, "", " ");
95
96 // y space
97 _yHeader = _pr.getColumnCount();
98 _pr.addColumn(false, " ", "");
99 _yIndex = _pr.getColumnCount();
100 for (size_t i = 0; i < _lat.getYDim(); ++i)
101 _pr.addColumn(false, i == 0 ? " " : " ");
102
103 // hits
104 _hitsHeader = _pr.getColumnCount();
105 _pr.addColumn(false, " ", "");
106 _hits = _pr.getColumnCount();
107 _pr.addColumn(false, " ", "");
108
109 // edges
110 _edgeHeader = _pr.getColumnCount();
111 _pr.addColumn(false, " ", "");
112 _edge = _pr.getColumnCount();
113 _pr.addColumn(false, " ", "");
114 }
115
118 _pr[_hitsHeader] << '\n';
119 _pr[_hits] << '\n';
120 _pr[_edgeHeader] << '\n';
121 _pr[_edge] << '\n';
122 }
123
125 const Plane& plane,
127 Mlfb& edge) {
129 _pr[_hitsHeader] << "hits\n";
130 _pr[_hits] << hits.getName() << '\n';
131 _pr[_edgeHeader] << "push to\n";
132 _pr[_edge] << edge.getName(plane) << '\n';
133 }
134
135 void addEmptyLine() {
136 for (size_t i = 0; i < _pr.getColumnCount(); ++i)
137 _pr[i] << '\n';
138 }
139
140 void print(FILE* out) {
141 ::print(out, _pr);
142 }
143
144 void print(ostream& out) {
145 out << _pr;
146 }
147
148 private:
150 _pr[_labelIndex] << neighbor.getName() << ':';
151 if (place != NoPlace)
152 _pr[_labelIndex] << ' ' << getPlaceCode(place);
153 _pr[_labelIndex] << '\n';
154
155 _pr[_yHeader] << "y=\n";
156 for (size_t i = 0; i < neighbor.getYDim(); ++i)
157 _pr[_yIndex + i] << neighbor.getY(i) << '\n';
158 _pr[_comma] << ",\n";
159 _pr[_hHeader] << "h=\n";
160 for (size_t i = 0; i < neighbor.getHDim(); ++i)
161 _pr[_hIndex + i] << neighbor.getH(i) << '\n';
162 }
163
164 const GrobLat& _lat;
165 ColumnPrinter _pr;
166 size_t _labelIndex;
167 size_t _comma;
168 size_t _yHeader;
169 size_t _yIndex;
170 size_t _hHeader;
171 size_t _hIndex;
172 size_t _hitsHeader;
173 size_t _hits;
174 size_t _edgeHeader;
175 size_t _edge;
176 };
177
179 Mlfb* min = &(mlfbs[0]);
180 for (size_t i = 1; i < mlfbs.size(); ++i) {
181 Mlfb& mlfb = mlfbs[i];
182 if (mlfb.dotDegree < min->dotDegree)
183 min = &mlfb;
184 }
185 cout << "The MLFB " << min->getName() << " has minimal rhs dot product.\n";
186 }
187
188 void printNeighbors(const GrobLat& lat) {
189 NeighborPrinter pr(lat);
190 pr.addNeighbor(Neighbor(lat));
191 for (size_t n = 0; n < lat.getNeighborCount(); ++n)
192 pr.addNeighbor(lat.getNeighbor(n));
193
194 fprintf(stdout, "\nThe %u neighbors in y-space and h-space are\n",
195 (unsigned int)lat.getNeighborCount());
196 pr.print(stdout);
197 fputc('\n', stdout);
198 }
199
200 class MlfbPtrCmp {
201 public:
202 MlfbPtrCmp(const Plane& plane): _plane(plane) {}
203 bool operator()(const Mlfb* a, const Mlfb* b) {
204 size_t ta = _plane.getType(*a);
205 size_t tb = _plane.getType(*b);
206 if (ta == 1) ta = 3;
207 if (tb == 1) tb = 3;
208 if (ta != tb)
209 return ta > tb;
210 else
211 return a->getOffset() < b->getOffset();
212 }
213
214 private:
215 const Plane& _plane;
216 };
217
218 void printMlfbs(vector<Mlfb>& mlfbs, const Plane& plane, const GrobLat& lat) {
219 vector<Mlfb*> order;
220 for (size_t i = 0; i < mlfbs.size(); ++i)
221 order.push_back(&(mlfbs[i]));
222 sort(order.begin(), order.end(), MlfbPtrCmp(plane));
223
224 cout << "\n\n";
225 for (size_t i = 0; i < order.size(); ++i) {
226 Mlfb& mlfb = *(order[i]);
227 cout << "*** MLFB " << mlfb.getName(plane) << " with rhs";
228 for (size_t var = 0; var < lat.getYDim(); ++var)
229 cout << ' ' << mlfb.getRhs()[var];
230 cout << " contains the neighbors\n";
231
232 NeighborPrinter pr(lat);
233 for (size_t j = 0; j < mlfb.getPointCount(); ++j) {
234 pr.addMlfbPoint(mlfb.getPoint(j),
235 plane,
236 mlfb.getHitsNeighbor(j),
237 *mlfb.getEdge(j));
238 }
239 pr.print(cout);
240 cout << "Its index is " << mlfb.index << ", its rhs dot product is " <<
241 mlfb.dotDegree << " and it has "
242 << plane.getType(mlfb) << " plane neighbors.\n\n";
243 }
244 }
245
246 void prSeq(const vector<SeqPos>& seq, const Plane& plane) {
247 cout << " Seq: ";
248 for (size_t i = 0; i < seq.size(); ++i)
249 cout << (i > 0 ? "->" : "") << seq[i].mlfb->getName(plane);
250 cout << endl;
251 }
252
253 void printPlane(vector<Mlfb>& mlfbs, const Plane& plane, const GrobLat& lat) {
254 cout << "The plane's null space is spanned by the rows of\n";
255 printIndentedMatrix(plane.nullSpaceBasis);
256 cout << '\n';
257
258 const vector<SeqPos>& flatSeq = plane.flatSeq;
259
260 cout << "The flat sequence has " << plane.flatIntervalCount
261 << (plane.flatIntervalCount == 1 ? " interval.\n" : " intervals.\n");
262
263 for (size_t i = plane.getMaxType(); i > 0; --i)
264 cout << "There are " << plane.getTypeCount(i)
265 << " MLFBs with " << i << " plane neighbors.\n";
266 cout << '\n';
267
268 const Tri& tri = plane.tri;
269 cout << "The non-flat triangle of the plane is "
270 << "{zero,"
271 << tri.getA().getName() << ',' << tri.getSum().getName() << "},"
272 << "{zero,"
273 << tri.getB().getName() << ',' << tri.getSum().getName() << "}.";
274
275 cout << "\nThe MLFBs containing {zero,"
276 << tri.getA().getName() << ',' << tri.getSum().getName() << "} are:";
277 for (size_t i = 0; i < tri.getASideMlfbs().size(); ++i)
278 cout << ' ' << tri.getASideMlfbs()[i]->getName(plane);
279
280 cout << "\nThe MLFBs containing {zero,"
281 << tri.getB().getName() << ',' << tri.getSum().getName() << "} are:";
282 for (size_t i = 0; i < tri.getBSideMlfbs().size(); ++i)
283 cout << ' ' << tri.getBSideMlfbs()[i]->getName(plane);
284
285 cout << "\nThe neighbors on the boundary of the body defined by {zero,"
286 << tri.getA().getName() << ','
287 << tri.getB().getName() << ','
288 << tri.getSum().getName() << "} are\n";
289 {
290 NeighborPrinter pr(lat);
291 for (size_t i = 0; i < tri.getNeighborsOnBoundary().size(); ++i)
292 pr.addNeighbor(tri.getNeighborsOnBoundary()[i]);
293 pr.print(stdout);
294 }
295 cout << "and those in the interior are\n";
296 {
297 NeighborPrinter pr(lat);
298 for (size_t i = 0; i < tri.getNeighborsInInterior().size(); ++i)
299 pr.addNeighbor(tri.getNeighborsInInterior()[i]);
300 pr.print(stdout);
301 }
302
303
305
306 const vector<const Mlfb*>& pivots = plane.pivots;
307
308 CHECK((pivots.size() % 2) == 0);
309
310 CHECK(!(tri.getASideMlfbs().size() == 2 &&
311 tri.getBSideMlfbs().size() == 2 &&
312 pivots.size() == 4 &&
313 flatSeq.size() > 0));
314
315 if (flatSeq.size() == 0 || pivots.size() != 4) {
316 cout << "(Falling back to general framework as no flats or more than 4 pivots)\n\n";
317
321
322 cout << "The left sequences are:\n";
323 for (size_t j = 0; j < left.size(); ++j)
324 prSeq(left[j], plane);
325
326 cout << "The right sequences are:\n";
327 for (size_t j = 0; j < right.size(); ++j)
328 prSeq(right[j], plane);
329
334 } else {
335 cout << endl;
336
339
340 computePivotSeqs(leftSeqs, *(pivots[0]), plane);
341 computePivotSeqs(rightSeqs, *(pivots[2]), plane);
342
343 cout << "The left sequences are:\n";
344 for (size_t j = 0; j < leftSeqs.size(); ++j)
346 cout << "The flat sequence is:\n";
347 prSeq(flatSeq, plane);
348 cout << "The right sequences are:\n";
349 for (size_t j = 0; j < rightSeqs.size(); ++j)
351
354 checkPlaneTri(lat, mlfbs, pivots, plane);
355 checkFlatSeq(flatSeq, lat, plane);
356 }
357 }
358
360 const vector<Mlfb>& mlfbs,
361 const vector<Plane>& planes) {
362 ofstream out("interior.dot");
363 out << "digraph G {\n" << flush;
364 for (size_t gen1 = 0; gen1 < lat.getNeighborCount(); ++gen1) {
365 Neighbor from = lat.getNeighbor(gen1);
366 out << ' ' << from.getName() << ";\n";
367 for (size_t gen2 = 0; gen2 < lat.getNeighborCount(); ++gen2) {
368 if (gen1 == gen2)
369 continue;
370 Neighbor to = lat.getNeighbor(gen2);
371 Neighbor sum = lat.getSum(from, to);
372 if (!sum.isValid() || !lat.isInterior(to, sum))
373 continue;
374 out << " " << from.getName() << " -> " << sum.getName()
375 << " [label=\"" << to.getName();
376 for (size_t m = 0; m < mlfbs.size(); ++m) {
377 const Mlfb& mlfb = mlfbs[m];
378 if (mlfb.hasPoint(to) && mlfb.hasPoint(sum)) {
379 out << ' ' << mlfb.getName();
380 }
381 }
382 if (lat.isPointFreeBody(from, sum)) {
383 for (size_t p = 0; p < planes.size(); ++p) {
384 const Plane& plane = planes[p];
385 if (plane.inPlane(sum) && plane.inPlane(to) && plane.inPlane(from))
386 out << " p" << (p + 1);
387 }
388 out << "\",style=dotted,arrowhead=empty";
389 } else
390 out << "\"";
391 out << "];\n";
392 }
393 }
394 out << "}\n";
395 }
396
397 void printScarfGraph(const vector<Mlfb>& mlfbs) {
398 ofstream out("graph.dot");
399 out << "graph G {\n";
400 for (size_t m = 0; m < mlfbs.size(); ++m) {
401 const Mlfb& mlfb = mlfbs[m];
402 out << " " << mlfb.getName() << "[label=\"";
403 out << mlfb.getName() << "\\nindex " << mlfb.index;
404 out << "\", shape=box];\n";
405
406 for (size_t e = 0; e < mlfb.edges.size(); ++e) {
407 size_t hits = mlfb.edgeHitsFacet[e];
408 if (mlfb.getOffset() < mlfb.edges[e]->getOffset())
409 continue;
410
411 out << " " << mlfb.getName()
412 << " -- " << mlfb.edges[e]->getName() << " [";
413 out << "headport=" << getEdgePos(hits) << ", ";
414 out << "tailport=" << getEdgePos(e) << "];\n";
415 }
416 }
417 out << "}\n";
418 }
419
421 ofstream out("ma3d");
422 out << "a={\n";
423 for (size_t m = 0; m < mlfbs.size(); ++m) {
424 out << " Graphics3D[{RGBColor[";
425 for (size_t i = 0; i < 3; ++i) {
426 if (i > 0)
427 out << ',';
428 out << "0." << (rand() % 10000);
429 }
430 out << "],Polygon[{";
431
432 Mlfb& mlfb = mlfbs[m];
433 for (size_t p = 1; p < mlfb.getPointCount(); ++p) {
434 if (p > 1)
435 out << ',';
436 out << '{';
437 for (size_t i = 0; i < lat.getHDim(); ++i) {
438 if (i > 1)
439 out << ',';
440 out << mlfb.getPoint(p).getH(i);
441 }
442 out << '}';
443 }
444 out << "}]}],\n";
445 }
446 out << " Graphics3D[Point[{0,0,0}]]\n};\ng=Show[{a},Boxed->False];\n";
447 }
448
450 Matrix matrix(thinPlanes.size(), 3);
451
452 for (size_t p = 0; p < thinPlanes.size(); ++p)
453 copyRow(matrix, p, thinPlanes[p].getNormal(), 0);
454
455 cout << "The " << thinPlanes.size() << " thin planes' normals are the rows of" << endl;
456 cout << matrix;
457 }
458}
459
461 Action
462 (staticGetName(),
463 "Display information about the input ideal.",
464 "This action is not ready for general use.\n\n"
465 "Display information about input Grobner basis of lattice.",
466 false),
467
468 _io(DataType::getSatBinomIdealType(), DataType::getMonomialIdealType()) {
469}
470
475
477 return false;
478}
479
481 cerr << "** Reading input " << endl;
482
486
488 SatBinomIdeal ideal;
489 ioFacade.readSatBinomIdeal(in, ideal);
490
492 {
494 ioFacade.readSatBinomIdeal(in, matrixIdeal);
496 }
497
498 cerr << "** Computing matrix and its nullspace" << endl;
499 cout << "Analysis of the "
500 << matrix.getRowCount() << " by " << matrix.getColCount()
501 << " matrix\n";
504
505 cerr << "** Computing h-space vectors" << endl;
506 GrobLat lat(matrix, ideal);
507
508 if (lat.hasZeroEntryY()) {
509 cout << "matrix not generic." << endl;
510 exit(2);
511 }
512
513 cerr << "** Computing MLFBs" << endl;
516
517 cerr << "** Computing thin planes" << endl;
520
521 cerr << "** Computing double triangle planes and associated structures" << endl;
524
525 size_t paraMlfbCount = 0;
526 for (size_t mlfb = 0; mlfb < mlfbs.size(); ++mlfb) {
527 if (mlfbs[mlfb].index == 0)
529 }
530
531 cerr << "** Producing output" << endl;
532
533 if (lat.hasZeroEntryY())
534 cout << "A neighbor has a zero entry.\n";
535 else
536 cout << "No neighbor has a zero entry.\n";
537
539
540 cout << "There are " << lat.getNeighborCount() << " neighbors excluding zero.\n";
541 cout << "There are " << mlfbs.size() << " MLFBs.\n";
542 cout << "There are " << paraMlfbCount << " parallelogram MLFBs.\n";
543 cout << "There are " << planes.size()
544 << " distinct double triangle planes.\n";
545 cout << "The sum of MLFB indexes is " << indexSum << ".\n";
546
547 CHECK(indexSum == 6 || indexSum == -6);
548
550
552
553 const vector<Neighbor>& nonSums = lat.getNonSums();
554 cout << "The " << nonSums.size() << " non-sum neighbors are:";
555 for (size_t i = 0; i < nonSums.size(); ++i)
556 cout << ' ' << nonSums[i].getName();
557 cout << endl;
558
560
561 for (size_t plane = 0; plane < planes.size(); ++plane) {
562 cout << "\n\n*** Plane " << (plane + 1)
563 << " of " << planes.size() << "\n\n";
566 }
569
571
575
578}
579
581 return "latanal";
582}
void print(FILE *out, const ColumnPrinter &pr)
void checkFlatSeq(const vector< SeqPos > &flatSeq, const GrobLat &lat, const Plane &plane)
void checkGraphOnPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
void checkPlanes(const vector< TriPlane > &thinPlanes, const vector< Plane > &dtPlanes)
void checkMlfbs(const vector< Mlfb > &mlfbs, const GrobLat &lat)
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)
char getPlaceCode(NeighborPlace place)
void computeMlfbs(vector< Mlfb > &mlfbs, const GrobLat &lat)
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)
void getThinPlanes(vector< TriPlane > &planes, const GrobLat &lat)
void checkNonSums(const GrobLat &lat)
void checkPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
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.
NeighborPlace
Definition LatticeAlgs.h:58
@ NoPlace
Definition LatticeAlgs.h:62
#define CHECK(X)
Definition LatticeAlgs.h:48
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 copyRow(Matrix &target, size_t targetRow, const Matrix &source, size_t sourceRow)
Copies row sourceRow from source to row targetRow of target.
Definition Matrix.cpp:246
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...
const char * getName() const
Definition Action.cpp:113
BoolParameter _printActions
Definition Action.h:68
virtual void obtainParameters(vector< Parameter * > &parameters)
Definition Action.cpp:133
void setPrefix(const string &prefix)
The intention of this class is to describe the different kinds of mathematical structures that Frobby...
Definition DataType.h:29
A lattice with associated Grobner basis/neighbors.
A facade for input and output of mathematical objects.
Definition IOFacade.h:39
void autoDetectInputFormat(Scanner &in)
If using the input format, this must be called before validating the ideals, since the auto detect fo...
const string & getInputFormat() const
void validateFormats() const
static const char * staticGetName()
virtual void obtainParameters(vector< Parameter * > &parameters)
virtual bool displayAction() const
This action isn't done yet, so it should not be displayed by the help action.
string getName() const
const vector< mpz_class > & getRhs() const
Neighbor getPoint(size_t offset) const
vector< Mlfb * > edges
mpq_class index
vector< size_t > edgeHitsFacet
Neighbor getHitsNeighbor(size_t indexParam) const
bool hasPoint(Neighbor n) const
mpz_class dotDegree
size_t getOffset() const
size_t getPointCount() const
const Mlfb * getEdge(size_t indexParam) const
const mpq_class & getH(size_t i) const
string getName() const
void obtainParameters(vector< Parameter * > &parameters)
Represents a saturated binomial ideal.
void getMatrix(Matrix &matrix) const
This class offers an input interface which is more convenient and for some purposes more efficient th...
Definition Scanner.h:50
const vector< const Mlfb * > & getASideMlfbs() const
const vector< Neighbor > & getNeighborsInInterior() const
Neighbor getB() const
const vector< Neighbor > & getNeighborsOnBoundary() const
const vector< const Mlfb * > & getBSideMlfbs() const
Neighbor getA() const
Neighbor getSum() const