LatticeAlgs.cpp
Go to the documentation of this file.
/* Frobby: Software for monomial ideal computations.
Copyright (C) 2011 Bjarke Hammersholt Roune (www.broune.com)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see http://www.gnu.org/licenses/.
*/
#include "stdinc.h"
#include "LatticeAlgs.h"
#include <stack>
void getThinPlanes(vector<TriPlane>& planes, const GrobLat& lat) {
planes.clear();
for (size_t a = 0; a <= lat.getNeighborCount(); ++a) {
Neighbor an(lat);
if (a < lat.getNeighborCount())
an = lat.getNeighbor(a);
for (size_t b = 0; b < lat.getNeighborCount(); ++b) {
for (size_t c = 0; c < lat.getNeighborCount(); ++c) {
TriPlane plane(an, lat.getNeighbor(b), lat.getNeighbor(c));
if (plane.isLine())
continue; // points are collinear
ASSERT(plane.isParallel(plane));
bool parallel = false;
for (size_t p = 0; p < planes.size(); ++p) {
if (plane.isParallel(planes[p])) {
parallel = true;
break;
}
}
if (parallel)
continue; // already got this plane
bool isThin = true;
for (size_t n = 0; n < lat.getNeighborCount(); ++n) {
if (!plane.closeToPlane(lat.getNeighbor(n))) {
isThin = false;
break;
}
}
if (isThin)
planes.push_back(plane);
}
}
}
}
/*
void checkParallelFaces(const vector<Mlfb>& mlfbs,
const vector<Plane> planes) {
vector<TriPlane> facePlanes;
for (size_t m = 0; m < mlfbs.size(); ++m) {
const Mlfb& mlfb = mlfbs[m];
cout << mlfb.getName() << facePlanes.size() << endl;
Neighbor a = mlfb.getPoint(0);
Neighbor b = mlfb.getPoint(1);
Neighbor c = mlfb.getPoint(2);
Neighbor d = mlfb.getPoint(3);
facePlanes.push_back(TriPlane(b, c, d));
facePlanes.push_back(TriPlane(a, c, d));
facePlanes.push_back(TriPlane(a, b, d));
facePlanes.push_back(TriPlane(a, b, c));
}
for (size_t i = 0; i < facePlanes.size(); ++i) {
const TriPlane& plane = facePlanes[i];
bool thin = false;
for (size_t p = 0; p < planes.size(); ++p) {
if (plane.isParallel(planes[p])) {
thin = true;
break;
}
}
if (thin)
continue;
size_t parallelCount = 0;
for (size_t j = 0; j < facePlanes.size(); ++j)
if (plane.isParallel(facePlanes[j])) {cout << ' ' << j << endl;
++parallelCount;}
cout << parallelCount << ' ' << plane.getNormal() << endl;
CHECK(parallelCount == 2); // not satisfied
}
}
*/
void checkPlanes(const vector<TriPlane>& thinPlanes,
const vector<Plane>& dtPlanes) {
CHECK(thinPlanes.size() == dtPlanes.size());
for (size_t thin = 0; thin < thinPlanes.size(); ++thin) {
bool parallel = false;
for (size_t dt = 0; dt < dtPlanes.size(); ++dt) {
if (thinPlanes[thin].isParallel(dtPlanes[dt])) {
parallel = true;
break;
}
}
CHECK(parallel);
}
bool found = false;
for (size_t dt = 0; dt < dtPlanes.size(); ++dt) {
const Plane& plane = dtPlanes[dt];
size_t sum = plane.tri.getASideMlfbs().size() +
plane.tri.getBSideMlfbs().size();
if (sum == 3)
found = true;
/*
if (plane.flatSeq.empty()) {
if (sum == 4)
found = true;
} else {
CHECK(sum == 3);
found = true;
}
*/
}
CHECK(dtPlanes.size() == 6 || found);
}
char getPlaceCode(NeighborPlace place) {
switch (place) {
case InPlane: return 'P';
case UnderPlane: return 'U';
case OverPlane: return 'O';
case NoPlace: return ' ';
default: return 'E';
}
}
size_t pushOutFacetPositive(size_t facetPushOut,
const vector<mpz_class>& rhs,
const GrobLat& lat) {
size_t onFacet = numeric_limits<size_t>::max();
mpq_class leastEntry;
for (size_t n = 0; n < lat.getNeighborCount(); ++n) {
for (size_t i = 0; i < lat.getYDim(); ++i) {
if (i == facetPushOut)
continue;
if (lat.getYMatrix()(n, i) >= rhs[i])
goto notOnFacet;
}
if (onFacet == numeric_limits<size_t>::max() ||
leastEntry > lat.getYMatrix()(n, facetPushOut)) {
leastEntry = lat.getYMatrix()(n, facetPushOut);
onFacet = n;
}
notOnFacet:;
}
return onFacet;
}
size_t pushOutFacetZero(const vector<mpz_class>& rhs, const GrobLat& lat) {
size_t onFacet = numeric_limits<size_t>::max();
mpq_class leastEntry;
for (size_t n = 0; n < lat.getNeighborCount(); ++n) {
for (size_t i = 1; i < lat.getYDim(); ++i)
if (-lat.getYMatrix()(n, i) >= rhs[i])
goto notOnFacet;
if (onFacet == numeric_limits<size_t>::max() ||
leastEntry > -lat.getYMatrix()(n, 0)) {
leastEntry = -lat.getYMatrix()(n, 0);
onFacet = n;
}
notOnFacet:;
}
return onFacet;
}
void computeRhs(vector<mpz_class>& rhs, const vector<Neighbor> points) {
ASSERT(!points.empty());
const GrobLat& lat = points[0].getGrobLat();
rhs.resize(lat.getYDim());
for (size_t var = 0; var < lat.getYDim(); ++var) {
rhs[var] = points[0].getY(var);
for (size_t p = 1; p < points.size(); ++p)
if (rhs[var] < points[p].getY(var))
rhs[var] = points[p].getY(var);
}
}
void Mlfb::reset(size_t offset, const vector<Neighbor>& points) {
ASSERT(!points.empty());
_points = points;
_offset = offset;
const GrobLat& lat = points[0].getGrobLat();
computeRhs(_rhs, points);
// order to have maxima along diagonal if possible.
if (getPointCount() == lat.getYDim()) {
for (size_t i = 0; i < lat.getYDim(); ++i)
for (size_t p = 0; p < getPointCount(); ++p)
if (getPoint(p).getY(i) == getRhs()[i])
swap(_points[i], _points[p]);
}
// Compute MLFB index.
if (getPointCount() - 1 == lat.getHDim())
{
Matrix mat(getPointCount() - 1, lat.getHDim());
for (size_t point = 1; point < getPointCount(); ++point)
for (size_t var = 0; var < lat.getHDim(); ++var)
mat(point - 1, var) = getPoint(point).getH(var);
index = determinant(mat);
}
if (getPointCount() == 4) {
Matrix mat(4, lat.getHDim());
for (size_t point = 0; point < getPointCount(); ++point)
for (size_t var = 0; var < lat.getHDim(); ++var)
mat(point, var) = getPoint(point).getH(var);
_isParallelogram = ::isParallelogram(mat);
} else
_isParallelogram = false;
}
void computeMlfbs(vector<Mlfb>& mlfbs, const GrobLat& lat) {
BigIdeal initialIdeal;
lat.getInitialIdeal(initialIdeal);
BigTermRecorder recorder;
SliceParams params;
SliceFacade facade(params, initialIdeal, recorder);
facade.computeIrreducibleDecomposition(true);
unique_ptr<BigIdeal> rhsesOwner = recorder.releaseIdeal();
BigIdeal& rhses = *rhsesOwner;
ASSERT(recorder.empty());
mlfbs.clear();
mlfbs.resize(rhses.getGeneratorCount());
vector<Neighbor> points;
for (size_t i = 0; i < mlfbs.size(); ++i) {
Mlfb& mlfb = mlfbs[i];
points.clear();
points.push_back(Neighbor(lat));
for (size_t gen = 0; gen < initialIdeal.getGeneratorCount(); ++gen) {
for (size_t var = 0; var < initialIdeal.getVarCount(); ++var)
if (initialIdeal[gen][var] > rhses[i][var])
goto skipIt;
points.push_back(Neighbor(lat, gen));
skipIt:;
}
mlfb.reset(i, points);
CHECK(rhses[i] == mlfb.getRhs());
}
Matrix nullSpaceBasis;
nullSpace(nullSpaceBasis, lat.getMatrix());
transpose(nullSpaceBasis, nullSpaceBasis);
// the basis is the rows of NullSpaceBasis at this point.
for (size_t m = 0; m < mlfbs.size(); ++m) {
Mlfb& mlfb = mlfbs[m];
if (mlfb.getPointCount() != lat.getYDim())
continue;
// Compute minInitialFacet.
mlfb.minInitialFacet = 0;
mpq_class minInitial = 0;
for (size_t i = 0; i < mlfb.getPointCount(); ++i) {
mpq_class initial = mlfb.getPoint(i).getY(0);
if (minInitial > initial) {
minInitial = initial;
mlfb.minInitialFacet = i;
}
}
// Compute dot degree.
if (nullSpaceBasis.getRowCount() == 1 &&
lat.getYDim() == nullSpaceBasis.getColCount()) {
mlfb.dotDegree = 0;
for (size_t r = 0; r < nullSpaceBasis.getRowCount(); ++r)
for (size_t c = 0; c < nullSpaceBasis.getColCount(); ++c)
mlfb.dotDegree += nullSpaceBasis(r, c) * mlfb.getRhs()[c];
}
// Compute Scarf edges.
mlfb.edges.resize(lat.getYDim());
mlfb.edgeHitsFacet.resize(lat.getYDim());
for (size_t facetPushIn = 0; facetPushIn < lat.getYDim(); ++facetPushIn) {
mpq_class secondLargest;
size_t facetPushOut = numeric_limits<size_t>::max();
for (size_t neigh = 0; neigh < lat.getYDim(); ++neigh) {
mpq_class entry = 0; // neigh == 0 represents zero.
if (neigh > 0)
entry = mlfb.getPoint(neigh).getY(facetPushIn);
if (neigh == facetPushIn) {
if (entry != mlfb.getRhs()[facetPushIn])
goto skipBecauseNotGeneric;
} else {
if (entry == secondLargest &&
facetPushOut != numeric_limits<size_t>::max())
goto skipBecauseNotGeneric;
if (entry > secondLargest ||
facetPushOut == numeric_limits<size_t>::max()) {
secondLargest = entry;
facetPushOut = neigh;
}
}
}
// ----------------------------------------------
// ** Case 1: facetPushIn > 0 and facetPushOut > 0
//
// We push in facetPushIn (discarding the non-zero neighbor
// previously on that facet) and hit a non-zero neighbor
// that is already on the MLFB. That neighbor now instead
// lies on facetPushIn and we push out facetPushOut in the
// straight forward way until it hits what will be the new
// neighbor on facetPushOut.
//
// ----------------------------------------------
// ** Case 2: facetPushIn > 0 and facetPushOut == 0
//
// We push in facetPushIn (discarding the non-zero neighbor
// previously on that facet) and hit zero. Zero now instead
// lies on facetPushIn and we push out facetPushOut. It will
// be pushing into the area on the opposite side from the
// half-set of neighbors that we are looking at, so to find
// the replacement neighbor to put on facetPushOut
// (i.e. facet zero) we need to consider the negative of the
// neighbors we have. When that neighbor -v has been found,
// we need to translate the whole body by +v so that zero
// will once again lie on facet zero.
//
// ----------------------------------------------
// ** case 3: facetPushIn == 0 and facetPushOut > 0
//
// We push in facetPushIn (discarding the zero previously on
// that facet) and hit a non-zero neighbor v on
// facetPushOut. Then v lies on facetPushIn (i.e. facet
// zero), so for the MLFB to have zero on facet 0 we need to
// translate it by -v. We then relax facetPushOut to find a
// replacement neighbor as in Case 1.
ASSERT(facetPushIn == 0 || secondLargest >= 0);
// In all cases the neighbor we hit moves to facetPushIn.
vector<mpz_class> rhs(mlfb.getRhs());
rhs[facetPushIn] = secondLargest;
if (facetPushIn == 0) {
// Case 3: the neighbor hit moves to facet 0 so translate
// the body to make that neighbor zero.
for (size_t i = 0; i < lat.getYDim(); ++i)
rhs[i] -= mlfb.getPoint(facetPushOut).getY(i);
}
if (facetPushOut > 0) {
// Case 1 or 3: push out in the usual way.
size_t newNeighbor = pushOutFacetPositive(facetPushOut, rhs, lat);
if (newNeighbor == numeric_limits<size_t>::max())
goto skipBecauseNotGeneric;
rhs[facetPushOut] = lat.getYMatrix()(newNeighbor, facetPushOut);
}
if (facetPushOut == 0) {
// Case 2: push out into negative neighbors
size_t newNeighbor = pushOutFacetZero(rhs, lat);
if (newNeighbor == numeric_limits<size_t>::max())
goto skipBecauseNotGeneric;
for (size_t i = 0; i < lat.getYDim(); ++i)
rhs[i] += lat.getYMatrix()(newNeighbor, i);
rhs[0] = 0;
}
// Find the MLFB with the right hand side that we have
// computed.
for (size_t mi = 0; mi < mlfbs.size(); ++mi) {
if (mlfbs[mi].getRhs() == rhs) {
mlfb.edges[facetPushIn] = &(mlfbs[mi]);
mlfb.edgeHitsFacet[facetPushIn] = facetPushOut;
goto foundMatch;
}
}
goto skipBecauseNotGeneric;
foundMatch:;
}
continue;
skipBecauseNotGeneric:
mlfb.edges.clear();
mlfb.edgeHitsFacet.clear();
}
}
SeqPos nextInSeq(SeqPos pos) {
size_t pushIn;
for (pushIn = 0;; ++pushIn) {
ASSERT(pushIn < 4);
if (pushIn != pos.fixFacet1 &&
pushIn != pos.fixFacet2 &&
pushIn != pos.comingFromFacet)
break;
}
size_t hits = pos.mlfb->getHitsFacet(pushIn);
ASSERT(hits != pushIn);
SeqPos next = pos;
next.mlfb = pos.mlfb->getEdge(pushIn);
next.comingFromFacet = hits;
if (pos.fixFacet1 == hits)
next.fixFacet1 = pushIn;
else if (pos.fixFacet2 == hits)
next.fixFacet2 = pushIn;
next.order();
return next;
}
SeqPos prevInSeq(SeqPos pos) {
return nextInSeq(pos.getReverse()).getReverse();
}
size_t computeFlatIntervalCount(const vector<SeqPos>& flatSeq) {
if (flatSeq.empty())
return 0u;
const Mlfb& leftFlat = *(flatSeq.front().mlfb);
size_t sumFacet = leftFlat.getMinInitialFacet();
size_t aFacet = 4;
for (size_t j = 0; j < 4; ++j) {
if (j != 0 && j != sumFacet) {
aFacet = j;
break;
}
}
CHECK(aFacet != 4);
size_t subSeqCount = 1;
for (size_t i = 1; i < flatSeq.size() - 1; ++i) {
const Mlfb& prev = *(flatSeq[i - 1].mlfb);
const Mlfb& flat = *(flatSeq[i].mlfb);
if (flat.getHitsFacet(aFacet) != prev.getHitsFacet(aFacet))
++subSeqCount;
}
return subSeqCount;
}
void computeFlatSeq(vector<SeqPos>& seq,
const vector<Mlfb>& mlfbs,
const Plane& plane) {
// ** compute left most flat
const Mlfb* leftFlat = 0;
for (size_t m = 0; m < mlfbs.size(); ++m) {
if (!plane.isFlat(mlfbs[m]))
continue;
const Mlfb* toLeft = mlfbs[m].getEdge(0);
if (!plane.isFlat(*toLeft)) {
CHECK(leftFlat == 0 || leftFlat == toLeft); // left flat unique
leftFlat = &(mlfbs[m]);
}
}
seq.clear();
if (leftFlat == 0) {
ASSERT(!plane.hasFlat());
return;
}
// ** go right as long as there is a flat there
SeqPos pos;
pos.mlfb = leftFlat;
while (plane.isFlat(*pos.mlfb)) {
ASSERT(seq.empty() || seq.back().mlfb != pos.mlfb);
seq.push_back(pos);
bool moved = false;
for (size_t facet = 1; facet < 4; ++facet) {
if (pos.mlfb->getEdge(facet)->getEdge(0) == pos.mlfb) {
pos.mlfb = pos.mlfb->getEdge(facet);
moved = true;
break;
}
}
if (!moved)
break;
}
}
void computePlanes(vector<Plane>& planes,
const GrobLat& lat,
vector<Mlfb>& mlfbs) {
const size_t neighborCount = lat.getNeighborCount();
for (size_t gen1 = 0; gen1 < neighborCount; ++gen1) {
for (size_t gen2 = gen1 + 1; gen2 < neighborCount; ++gen2) {
Neighbor a = lat.getNeighbor(gen1);
Neighbor b = lat.getNeighbor(gen2);
Neighbor sum = lat.getSum(a, b);
if (!sum.isValid())
continue;
if (!lat.isPointFreeBody(a, sum))
continue;
if (!lat.isPointFreeBody(b, sum))
continue;
if (lat.isPointFreeBody(a, b, sum))
continue; // only looking for non-flat double triangles right now
Matrix rowAB(2, lat.getHDim());
copyRow(rowAB, 0, lat.getHMatrix(), gen1);
copyRow(rowAB, 1, lat.getHMatrix(), gen2);
planes.push_back(Plane(a, b, sum, mlfbs, lat));
Plane& plane = planes.back();
plane.rowAB = rowAB;
nullSpace(plane.nullSpaceBasis, rowAB);
Matrix prod;
product(prod, lat.getHMatrix(), plane.nullSpaceBasis);
plane.neighborPlace.resize(lat.getNeighborCount());
for (size_t gen = 0; gen < lat.getNeighborCount(); ++gen) {
mpq_class& value = prod(gen, 0);
NeighborPlace place = InPlane;
if (value < 0)
place = UnderPlane;
else if (value > 0)
place = OverPlane;
plane.neighborPlace[gen] = place;
}
transpose(plane.nullSpaceBasis);
}
}
for (size_t p = 0; p < planes.size(); ++p) {
Plane& plane = planes[p];
for (size_t i = 0; i < mlfbs.size(); ++i)
plane.typeCounts[plane.getType(mlfbs[i])] += 1;
computeFlatSeq(plane.flatSeq, mlfbs, plane);
computePivots(plane.pivots, mlfbs, plane, plane.flatSeq);
plane.flatIntervalCount = computeFlatIntervalCount(plane.flatSeq);
}
}
Tri::Tri(Neighbor a, Neighbor b, Neighbor sum,
const vector<Mlfb>& mlfbs, const GrobLat& lat):
_a(a), _b(b), _sum(sum) {
// find MLFBs containing {0,a,a+b}
for (size_t m = 0; m < mlfbs.size(); ++m)
if (mlfbs[m].hasPoint(a) && mlfbs[m].hasPoint(sum))
_aSideMlfbs.push_back(&(mlfbs[m]));
// find MLFBs containing {0,b,a+b}
for (size_t m = 0; m < mlfbs.size(); ++m)
if (mlfbs[m].hasPoint(b) && mlfbs[m].hasPoint(sum))
_bSideMlfbs.push_back(&(mlfbs[m]));
// find additional neighbors in the body defined by {0,a,b,a+b}
vector<Neighbor> points;
points.push_back(Neighbor(lat)); // add zero;
points.push_back(a);
points.push_back(b);
points.push_back(sum);
vector<mpz_class> rhs;
computeRhs(rhs, points);
_boundary.push_back(Neighbor(lat)); // add zero
for (size_t n = 0; n < lat.getNeighborCount(); ++n) {
Neighbor neighbor = lat.getNeighbor(n);
bool boundary = true;
bool interior = true;
for (size_t i = 0; i < rhs.size(); ++i) {
if (neighbor.getY(i) == rhs[i])
interior = false;
else if (neighbor.getY(i) > rhs[i]) {
interior = false;
boundary = false;
break;
}
}
if (interior)
_interior.push_back(neighbor);
else if (boundary)
_boundary.push_back(neighbor);
}
}
void check0Graph(const vector<Mlfb>& mlfbs) {
vector<bool> ok(mlfbs.size());
bool sawFlat = false;
for (size_t i =0 ; i < mlfbs.size(); ++i) {
ok[i] = (mlfbs[i].index == 0);
if (ok[i])
sawFlat = true;
}
if (!sawFlat)
return;
while (true) {
bool done = true;
for (size_t i = 0; i < mlfbs.size(); ++i) {
if (!ok[i]) {
size_t to = mlfbs[i].getEdge(0)->getOffset();
if (ok[to]) {
done = false;
ok[i] = true;
}
}
}
if (done)
break;
}
for (size_t i = 0; i < mlfbs.size(); ++i) {
CHECK(ok[i]);
}
}
void checkMlfbs(const vector<Mlfb>& mlfbs, const GrobLat& lat) {
CHECK(mlfbs.size() == lat.getNeighborCount() - 1);
for (size_t m = 0; m < mlfbs.size(); ++m) {
CHECK(mlfbs[m].isParallelogram() == (mlfbs[m].index == 0));
}
}
void checkDoubleTrianglePlanes(const vector<Plane>& planes,
const GrobLat& lat,
const vector<Mlfb>& mlfbs) {
// Check no two planes are parallel. Otherwise there would a plane
// with two non-flat double triangles in it.
for (size_t p1 = 0; p1 < planes.size(); ++p1) {
for (size_t p2 = 0; p2 < p1; ++p2) {
CHECK(!hasSameRowSpace(planes[p1].rowAB, planes[p2].rowAB));
}
}
// Check that all parallelograms lie in a plane. Otherwise there
// would be a plane defined by a flat that does not have a double
// triangle in it.
for (size_t m = 0; m < mlfbs.size(); ++m) {
if (mlfbs[m].isParallelogram()) {
bool liesInSomePlane = false;
for (size_t p = 0; p < planes.size(); ++p) {
if (planes[p].isFlat(mlfbs[m])) {
liesInSomePlane = true;
break;
}
}
CHECK(liesInSomePlane);
}
}
bool multipleIntervals = false;
bool anyFlat = false;
bool flatWith4Pivots = false;
for (size_t p = 0; p < planes.size(); ++p) {
if (planes[p].flatIntervalCount > 1)
multipleIntervals = true;
if (planes[p].hasFlat()) {
anyFlat = true;
if (planes[p].pivots.size() == 4)
flatWith4Pivots = true;
}
}
if (multipleIntervals) {
ASSERT(anyFlat);
//CHECK(flatWith4Pivots);
CHECK(planes.size() == 1);
}
if (planes.size() == 6) {
CHECK(!anyFlat);
CHECK(planes.size() == 6);
for (size_t p = 0; p < planes.size(); ++p) {
CHECK(planes[p].pivots.size() == 4);
}
CHECK(lat.getNeighborCount() == 7);
CHECK(mlfbs.size() == 6);
}
if (anyFlat) {
//check0Graph(mlfbs);
//CHECK(flatWith4Pivots);
CHECK(planes.size() < 6);
}
}
void checkPlane(const Plane& plane, const vector<Mlfb>& mlfbs) {
for (size_t i = 0; i < mlfbs.size(); ++i) {
if (plane.isPivot(mlfbs[i])) {
CHECK(mlfbs[i].index == -1 || mlfbs[i].index == 1);
} else if (plane.isFlat(mlfbs[i])) {
CHECK(mlfbs[i].index == 0);
}
}
}
size_t pivotToFlatFacet(const Mlfb& pivot, const Plane& plane) {
size_t facet = 4;
for (size_t push = 0; push < 4; ++push) {
if (plane.isFlat(*pivot.getEdge(push))) {
CHECK(facet == 4); // adjacent to no more than one flat
facet = push;
}
}
CHECK(facet != 4); // adjacent to at least one flat
return facet;
}
bool disjointSeqs(const vector<SeqPos>& a, const vector<SeqPos>& b) {
for (size_t i = 0; i < a.size(); ++i)
for (size_t j = 0; j < b.size(); ++j)
if (a[i].mlfb == b[j].mlfb)
return false;
return true;
}
void computePivots(vector<const Mlfb*>& pivots,
const vector<Mlfb>& mlfbs,
const Plane& plane,
const vector<SeqPos>& flatSeq) {
pivots.clear();
for (size_t m = 0; m < mlfbs.size(); ++m)
if (plane.isPivot(mlfbs[m]))
pivots.push_back(&(mlfbs[m]));
if (pivots.size() != 4 || flatSeq.empty())
return; // no idea about proper order in this case
pivots.clear();
// use flat sequence to impose correct order
size_t sumFacet = flatSeq.front().mlfb->getMinInitialFacet();
pivots.push_back(flatSeq.front().mlfb->getEdge(0));
pivots.push_back(flatSeq.front().mlfb->getEdge(sumFacet));
sumFacet = flatSeq.back().mlfb->getMinInitialFacet();
for (size_t i = 0; i < 4; ++i)
if (i != 0 && i != sumFacet)
pivots.push_back(flatSeq.back().mlfb->getEdge(i));
}
void computeSeqs(vector<vector<SeqPos> >& left,
vector<vector<SeqPos> >& right,
const vector<Mlfb>& mlfbs,
const Plane& plane) {
vector<vector<SeqPos> > seqs;
for (size_t m = 0; m < mlfbs.size(); ++m) {
if (!plane.isPivot(mlfbs[m]))
continue;
const Mlfb& p = mlfbs[m];
for (size_t i = 0; i < 4; ++i) {
const Mlfb& e = *(p.getEdge(i));
if (plane.isPivot(e) || plane.isFlat(e))
continue;
bool doneBefore = false;
for (size_t s = 0; s < seqs.size(); ++s) {
if (*(seqs[s][seqs[s].size() - 1].mlfb) == p &&
*(seqs[s][seqs[s].size() - 2].mlfb) == e) {
doneBefore = true;
break;
}
}
if (doneBefore)
continue;
size_t prevFacet;
for (prevFacet = 0; prevFacet < 4; ++prevFacet)
if (*(e.getEdge(prevFacet)) == p)
break;
ASSERT(prevFacet < 4);
NeighborPlace place = plane.getPlace(e.getPoint(prevFacet));
size_t nextFacet;
for (nextFacet = 0; nextFacet < 4; ++nextFacet) {
if (nextFacet != prevFacet &&
place == plane.getPlace(e.getPoint(nextFacet)))
break;
}
SeqPos pos = prevInSeq(SeqPos(&e, nextFacet, prevFacet));
seqs.resize(seqs.size() + 1);
vector<SeqPos>& seq = seqs.back();
seq.push_back(pos);
do {
pos = nextInSeq(pos);
seq.push_back(pos);
} while (!(plane.isPivot(*pos.mlfb)));
}
}
CHECK(!seqs.empty());
ASSERT(!seqs.front().empty());
// now we've got all the sequences. Time to look at the sides.
stack<const Mlfb*> pending;
// put a side pivot on the left arbitrarily and explore the
// connected component
vector<bool> leftSeen(mlfbs.size());
pending.push(seqs.front().front().mlfb);
while (!pending.empty()) {
const Mlfb& m = *pending.top();
pending.pop();
if (leftSeen[m.getOffset()])
continue;
leftSeen[m.getOffset()] = true;
for (size_t s = 0; s < seqs.size(); ++s) {
if (*(seqs[s].front().mlfb) == m)
pending.push(seqs[s].back().mlfb);
if (*(seqs[s].back().mlfb) == m)
pending.push(seqs[s].front().mlfb);
}
}
// find a non-left side pivot
size_t m;
for (m = 0; m < mlfbs.size(); ++m)
if (plane.isSidePivot(mlfbs[m]) && !leftSeen[m])
break;
CHECK(m < mlfbs.size());
// put the non-left pivot on the right and explore the connected
// component
vector<bool> rightSeen(mlfbs.size());
pending.push(&(mlfbs[m]));
while (!pending.empty()) {
const Mlfb& pm = *pending.top();
pending.pop();
if (rightSeen[pm.getOffset()])
continue;
rightSeen[pm.getOffset()] = true;
for (size_t s = 0; s < seqs.size(); ++s) {
if (*(seqs[s].front().mlfb) == pm)
pending.push(seqs[s].back().mlfb);
if (*(seqs[s].back().mlfb) == pm)
pending.push(seqs[s].front().mlfb);
}
}
left.clear();
right.clear();
for (size_t s = 0; s < seqs.size(); ++s) {
const size_t offset = seqs[s].front().mlfb->getOffset();
if (leftSeen[offset]) {
CHECK(!rightSeen[offset]);
left.push_back(seqs[s]);
} else if (rightSeen[offset])
right.push_back(seqs[s]);
}
}
void computePivotSeqs(vector<vector<SeqPos> >& seqs,
const Mlfb& pivot,
const Plane& plane) {
ASSERT(plane.pivots.size() == 4);
ASSERT(plane.isPivot(pivot));
size_t flatFacet = pivotToFlatFacet(pivot, plane);
seqs.clear();
for (size_t facet = 0; facet < 4; ++facet) {
if (facet == flatFacet)
continue;
seqs.resize(seqs.size() + 1);
vector<SeqPos>& seq = seqs.back();
SeqPos pos(&pivot, facet, flatFacet);
seq.push_back(pos);
do {
pos = nextInSeq(pos);
seq.push_back(pos);
} while (!(plane.isPivot(*pos.mlfb)));
}
}
void checkSeq(vector<bool>& seenOnSide,
const vector<SeqPos>& seq,
const Plane& plane) {
CHECK(seq.size() >= 3); // each seq must have at least one 2-2.
CHECK(plane.isSidePivot(*(seq.front().mlfb))); // start is a pivot
CHECK(plane.isSidePivot(*(seq.back().mlfb))); // end is a pivot
CHECK(seq.front().mlfb != seq.back().mlfb); // no loops
for (size_t m = 1; m < seq.size() - 1 ; ++m) {
const Mlfb* prev = seq[m - 1].mlfb;
const Mlfb* current = seq[m].mlfb;
const Mlfb* next = seq[m + 1].mlfb;
const SeqPos& pos = seq[m];
// ** Check on 2-2 appears twice and update seenOnSide
CHECK(!seenOnSide[current->getOffset()]);
seenOnSide[current->getOffset()] = true;
// ** middle elements are 2-2s.
CHECK(plane.is22(*current));
// ** SeqPos fields agrees with sequence order
size_t prevFacet = pos.getBackFacet();
size_t nextFacet = pos.getForwardFacet();
CHECK(current->getEdge(prevFacet) == prev);
CHECK(current->getEdge(nextFacet) == next);
// ** in-coming and out-going facets are on same plane
CHECK(plane.getPlace(current->getPoint(prevFacet)) ==
plane.getPlace(current->getPoint(nextFacet)));
}
}
void checkSide(vector<bool>& pivotOnSide,
const vector<vector<SeqPos> >& side,
const Plane& plane,
const vector<Mlfb>& mlfbs) {
CHECK(side.size() == 2 || side.size() == 3);
vector<bool> seenOnSide(mlfbs.size());
for (size_t s = 0; s < side.size(); ++s){
// check sequence local properties and update seen
checkSeq(seenOnSide, side[s], plane);
// compute onSide
pivotOnSide[side[s].front().mlfb->getOffset()] = true;
pivotOnSide[side[s].back().mlfb->getOffset()] = true;
}
/* // must have seen all 2-2s on this side.
for (size_t m = 0; m < mlfbs.size(); ++m)
if (plane.is22(mlfbs[m]))
CHECK(seenOnSide[m]);*/
// 2,3 or 4 sidepivots
size_t sidePivots = 0;
for (size_t m = 0; m < mlfbs.size(); ++m)
if (pivotOnSide[m])
++sidePivots;
CHECK(sidePivots == 2 || sidePivots == 3 || sidePivots == 4);
}
void checkSeqs(const vector<vector<SeqPos> >& left,
const vector<vector<SeqPos> >& right,
const Plane& plane,
const vector<Mlfb>& mlfbs) {
vector<bool> isLeftPivot(mlfbs.size());
checkSide(isLeftPivot, left, plane, mlfbs);
vector<bool> isRightPivot(mlfbs.size());
checkSide(isRightPivot, right, plane, mlfbs);
// all side pivots are on one and only one side
for (size_t m = 0; m < mlfbs.size(); ++m) {
if (plane.isSidePivot(mlfbs[m]))
CHECK((isLeftPivot[m] + isRightPivot[m]) == 1);
else
CHECK((isLeftPivot[m] + isRightPivot[m]) == 0);
}
}
void checkMiddle(const Plane& plane,
const vector<Mlfb>& mlfbs) {
// ** check that the subgraph of pivots and flat is connected.
vector<bool> seen(mlfbs.size());
stack<const Mlfb*> pending;
// find a pivot or flat
size_t m;
for (m = 0; m < mlfbs.size(); ++m)
if (plane.isFlat(mlfbs[m]) || plane.isPivot(mlfbs[m]))
break;
ASSERT(m < mlfbs.size());
// explore the graph of pivots and flats that is connected to the
// one we found.
pending.push(&(mlfbs[m]));
while (!pending.empty()) {
const Mlfb& mlfb = *(pending.top());
pending.pop();
if (seen[mlfb.getOffset()])
continue;
seen[mlfb.getOffset()] = true;
for (size_t i = 0; i < 4; ++i)
pending.push(mlfb.getEdge(i));
}
// check that we have reached all pivots and flats
for (m = 0; m < mlfbs.size(); ++m)
if (plane.isFlat(mlfbs[m]) || plane.isPivot(mlfbs[m]))
CHECK(seen[m]);
}
void checkGraphOnPlane(const Plane& plane,
const vector<Mlfb>& mlfbs) {
// no flat is adjacent to a 2-2
for (size_t m = 0; m < mlfbs.size(); ++m) {
const Mlfb& mlfb = mlfbs[m];
if (plane.isFlat(mlfb))
for (size_t i = 0; i < 4; ++i)
CHECK(!plane.is22(*(mlfb.getEdge(i))));
}
// parallelograms can't be flats and non-flat parallelograms are not
// adjacent to a flat
for (size_t m = 0; m < mlfbs.size(); ++m) {
const Mlfb& mlfb = mlfbs[m];
if (mlfb.isParallelogram()) {
CHECK(!plane.isPivot(mlfb));
if (!plane.isFlat(mlfb)) {
for (size_t i = 0; i < 4; ++i) {
const Mlfb& adj = *(mlfb.getEdge(i));
CHECK(!plane.isFlat(adj));
}
}
}
}
}
void checkDoubleTriangle(const Plane& plane,
const vector<Mlfb>& mlfbs) {
size_t aSideCount = plane.tri.getASideMlfbs().size();
size_t bSideCount = plane.tri.getBSideMlfbs().size();
CHECK(aSideCount == 1 || aSideCount == 2);
CHECK(bSideCount == 1 || bSideCount == 2);
for (size_t m = 0; m < aSideCount; ++m) {
const Mlfb& mlfb = *plane.tri.getASideMlfbs()[m];
CHECK(plane.isFlat(mlfb) || plane.isPivot(mlfb));
}
for (size_t m = 0; m < bSideCount; ++m) {
const Mlfb& mlfb = *plane.tri.getBSideMlfbs()[m];
CHECK(plane.isFlat(mlfb) || plane.isPivot(mlfb));
}
}
void checkGraph(const vector<Mlfb>& mlfbs) {
// All non-parallelograms have out-degree 2. Parallelograms have
// out-degree equal to 4 minus the number of other adjacent
// parallelograms.
for (size_t m = 0; m < mlfbs.size(); ++m) {
const Mlfb& mlfb = mlfbs[m];
set<size_t> adjParas;
set<size_t> adjNodes;
for (size_t i = 0; i < 4; ++i) {
const Mlfb& adj = *(mlfb.getEdge(i));
adjNodes.insert(adj.getOffset());
if (adj.isParallelogram())
adjParas.insert(adj.getOffset());
}
const size_t outDegree = adjNodes.size();
if (!mlfb.isParallelogram())
CHECK(outDegree == 4);
else
CHECK(outDegree == 4 - adjParas.size());
}
// if there is an edge (a,b) then there is also an edge (b,a).
for (size_t m = 0; m < mlfbs.size(); ++m) {
const Mlfb& mlfb = mlfbs[m];
for (size_t i = 0; i < 4; ++i) {
size_t hitsFacet = mlfb.getHitsFacet(i);
const Mlfb& adj = *(mlfb.getEdge(i));
CHECK(mlfb == *(adj.getEdge(hitsFacet)));
}
}
}
void checkPivotSeqs(vector<vector<SeqPos> >& pivotSeqs,
const Plane& plane,
const vector<Mlfb>& mlfbs,
const vector<SeqPos>& flatSeq) {
CHECK(pivotSeqs.size() == 3);
CHECK(pivotSeqs[0].size() >= 2);
const Mlfb* pivot1 = pivotSeqs[0].front().mlfb;
const Mlfb* pivot2 = pivotSeqs[0].back().mlfb;
CHECK(plane.isPivot(*pivot1));
CHECK(plane.isPivot(*pivot2));
bool foundPlace = false;
NeighborPlace place;
for (size_t i = 0; i < 3; ++i) {
CHECK(pivotSeqs[i].size() >= 2);
// ** all sequences on same side between same pivots
CHECK((pivotSeqs[i].front().mlfb == pivot1 &&
pivotSeqs[i].back().mlfb == pivot2) ||
(pivotSeqs[i].front().mlfb == pivot2 &&
pivotSeqs[i].back().mlfb == pivot1));
for (size_t j = 1; j < pivotSeqs[i].size() - 1; ++j) {
const Mlfb* prev = pivotSeqs[i][j - 1].mlfb;
const Mlfb* current = pivotSeqs[i][j].mlfb;
const Mlfb* next = pivotSeqs[i][j + 1].mlfb;
// ** middle mlfbs are 2-2's
CHECK(plane.getType(*current) == 2);
// ** SeqPos agrees with sequence order
const SeqPos& pos = pivotSeqs[i][j];
size_t prevFacet = pos.getBackFacet();
size_t nextFacet = pos.getForwardFacet();
CHECK(current->getEdge(prevFacet) == prev);
CHECK(current->getEdge(nextFacet) == next);
// ** in-coming and out-going facets are on same plane
CHECK(plane.getPlace(current->getPoint(prevFacet)) ==
plane.getPlace(current->getPoint(nextFacet)));
// ** active facets are the same for all sequences on same side
if (!foundPlace) {
place = plane.getPlace(current->getPoint(prevFacet));
}
CHECK(place == plane.getPlace(current->getPoint(prevFacet)));
}
}
// ** Check that the sequences are a partition of the set of 2-2 mlfbs
vector<bool> seen(mlfbs.size());
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 1; j < pivotSeqs[i].size() - 1; ++j) {
CHECK(!seen[pivotSeqs[i][j].mlfb->getOffset()]); // sequences must be disjoint
seen[pivotSeqs[i][j].mlfb->getOffset()] = true;
}
}
/* for (size_t i = 0; i < mlfbs.size(); ++i) {
if (plane.isPivot(mlfbs[i]) || plane.isFlat(mlfbs[i]))
continue;
CHECK(seen[i]); // all 2-2s in some sequence
}*/
}
void checkNonSums(const GrobLat& lat) {
const vector<Neighbor>& nonSums = lat.getNonSums();
CHECK(nonSums.size() == 3 || nonSums.size() == 4);
if (nonSums.size() == 3) {
Matrix mat(3, 3);
for (size_t ns = 0; ns < 3; ++ns)
for (size_t var = 0; var < 3; ++var)
mat(ns, var) = nonSums[ns].getH(var);
mpq_class det = determinant(mat);
CHECK(det == 1 || det == -1);
} else {
Matrix mat(4, 3);
for (size_t ns = 0; ns < 4; ++ns)
for (size_t var = 0; var < 3; ++var)
mat(ns, var) = nonSums[ns].getY(var);
CHECK(isParallelogram(mat));
return; // not checking this as it seems to not be true
mpq_class areaSq = getParallelogramAreaSq(mat);
CHECK(areaSq == 1);
}
}
void checkFlatSeq(const vector<SeqPos>& flatSeq,
const GrobLat& lat,
const Plane& plane) {
if (flatSeq.empty())
return;
size_t sumf = flatSeq.front().mlfb->getMinInitialFacet();
CHECK(sumf != 0);
size_t af = 4;
for (size_t j = 0; j < 4; ++j) {
if (j != 0 && j != sumf) {
af = j;
break;
}
}
size_t bf = 4;
for (size_t j = 0; j < 4; ++j) {
if (j != 0 && j != sumf && j != af) {
bf = j;
break;
}
}
for (size_t i = 0; i < flatSeq.size(); ++i) {
const Mlfb& flat = *(flatSeq[i].mlfb);
Neighbor sum = flat.getPoint(sumf);
Neighbor a = flat.getPoint(af);
Neighbor b = flat.getPoint(bf);
CHECK(sumf == flat.getMinInitialFacet());
CHECK(lat.getSum(a, b) == sum);
// right-going edges
if (i < flatSeq.size() - 1) {
// not at right end
CHECK(flat.getEdge(af) == flat.getEdge(bf));
const Mlfb& next = *(flatSeq[i + 1].mlfb);
CHECK(flat.getEdge(af) == &next);
if (flat.getHitsFacet(af) == sumf) {
CHECK(flat.getHitsFacet(bf) == 0);
CHECK(next.hasPoint(b));
CHECK(next.hasPoint(sum));
CHECK(next.hasPoint(lat.getSum(sum, b)));
} else {
CHECK(flat.getHitsFacet(af) == 0);
CHECK(flat.getHitsFacet(bf) == sumf);
CHECK(next.hasPoint(a));
CHECK(next.hasPoint(sum));
CHECK(next.hasPoint(lat.getSum(sum, a)));
}
} else {
// at right end
CHECK(plane.isPivot(*flat.getEdge(af)));
CHECK(plane.isPivot(*flat.getEdge(bf)));
CHECK(flat.getEdge(af) != flat.getEdge(bf));
}
// left-going edges
if (i > 0) {
// not at left end
CHECK(flat.getEdge(0) == flat.getEdge(sumf));
const Mlfb& prev = *(flatSeq[i - 1].mlfb);
CHECK(flat.getEdge(0) == &prev);
if (flat.getHitsFacet(0) == af) {
CHECK(flat.getHitsFacet(sumf) == bf);
} else {
CHECK(flat.getHitsFacet(0) == bf);
CHECK(flat.getHitsFacet(sumf) == af);
}
} else {
// at left end
CHECK(plane.isPivot(*flat.getEdge(0)));
CHECK(plane.isPivot(*flat.getEdge(sumf)));
CHECK(flat.getEdge(0) != flat.getEdge(sumf));
}
}
}
void checkPlaneTri(const GrobLat& lat,
const vector<Mlfb>& mlfbs,
const vector<const Mlfb*>& pivots,
const Plane& plane) {
const Tri& tri = plane.tri;
Neighbor a = tri.getA();
Neighbor b = tri.getB();
Neighbor sum = tri.getSum();
// ** tri is not a flat
for (size_t i = 0; i < mlfbs.size(); ++i) {
if (plane.isFlat(mlfbs[i])) {
CHECK(!mlfbs[i].hasPoint(a) ||
!mlfbs[i].hasPoint(b) ||
!mlfbs[i].hasPoint(sum));
}
}
// ** find unique non-flat MLFB with {0,a,a+b}
const Mlfb* mlfbA = 0;
for (size_t i = 0; i < mlfbs.size(); ++i) {
if (!plane.isFlat(mlfbs[i]) &&
mlfbs[i].hasPoint(a) &&
mlfbs[i].hasPoint(sum)) {
CHECK(mlfbA == 0);
mlfbA = &(mlfbs[i]);
}
}
CHECK(mlfbA != 0);
// ** find unique non-flat MLFB with {0,b,a+b}
const Mlfb* mlfbB = 0;
for (size_t i = 0; i < mlfbs.size(); ++i) {
if (!plane.isFlat(mlfbs[i]) &&
mlfbs[i].hasPoint(b) &&
mlfbs[i].hasPoint(sum)) {
CHECK(mlfbB == 0);
mlfbB = &(mlfbs[i]);
}
}
CHECK(mlfbB != 0);
// ** both parts must be pivots
CHECK(plane.isPivot(*mlfbA));
CHECK(plane.isPivot(*mlfbB));
// ** the two pivots must be on the left
CHECK((mlfbA == pivots[0] && mlfbB == pivots[1]) ||
(mlfbA == pivots[1] && mlfbB == pivots[0]));
}
const char* getEdgePos(size_t index) {
switch (index) {
case 0: return "sw";
case 1: return "se";
case 2: return "ne";
case 3: return "nw";
default: ASSERT(false);
return "ERROR";
}
}
mpq_class getIndexSum(const vector<Mlfb>& mlfbs) {
mpq_class sum;
for (size_t m = 0; m < mlfbs.size(); ++m)
sum += mlfbs[m].index;
return sum;
}
SeqPos::SeqPos() {}
SeqPos::SeqPos(const Mlfb* mlfbParam, size_t nextFacet, size_t previousFacet) {
ASSERT(mlfbParam != 0);
ASSERT(nextFacet != previousFacet);
ASSERT(nextFacet < 4);
ASSERT(previousFacet < 4);
mlfb = mlfbParam;
comingFromFacet = previousFacet;
for (size_t f = 0; f < 4; ++f)
if (f != previousFacet && f != nextFacet)
fixFacet1 = f;
for (size_t f = 0; f < 4; ++f)
if (f != previousFacet && f != nextFacet && f != fixFacet1)
fixFacet2 = f;
}
size_t SeqPos::getForwardFacet() const {
for (size_t i = 0; ; ++i) {
ASSERT(i < 4);
if (i != fixFacet1 && i != fixFacet2 && i != comingFromFacet)
return i;
}
}
size_t SeqPos::getBackFacet() const {
return comingFromFacet;
}
SeqPos SeqPos::getReverse() const {
size_t to;
for (to = 0; to < 4; ++to) {
if (to != fixFacet1 && to != fixFacet2 && to != comingFromFacet)
break;
}
ASSERT(to != 4);
SeqPos reverse = *this;
reverse.comingFromFacet = to;
return reverse;
}
void SeqPos::order() {
if (fixFacet1 > fixFacet2)
swap(fixFacet1, fixFacet2);
}
bool SeqPos::operator<(const SeqPos& pos) const {
if (mlfb->getOffset() < pos.mlfb->getOffset())
return true;
if (mlfb->getOffset() > pos.mlfb->getOffset())
return false;
if (fixFacet1 < pos.fixFacet1)
return true;
if (fixFacet1 > pos.fixFacet1)
return false;
if (fixFacet2 < pos.fixFacet2)
return true;
if (fixFacet2 > pos.fixFacet2)
return false;
if (comingFromFacet < pos.comingFromFacet)
return true;
if (comingFromFacet > pos.comingFromFacet)
return false;
return false;
}
Neighbor::Neighbor():
_lat(0), _row(0) {
}
Neighbor::Neighbor(const GrobLat& lat):
_lat(&lat), _row(lat.getNeighborCount() + 1) {
}
Neighbor::Neighbor(const GrobLat& lat, const size_t row):
_lat(&lat), _row(row) {
}
const mpq_class& Neighbor::getH(size_t i) const {
ASSERT(isValid());
ASSERT(i < getHDim());
if (isZero())
return _lat->getZero();
else
return _lat->getHMatrix()(_row, i);
}
size_t Neighbor::getHDim() const {
ASSERT(isValid());
return _lat->getHDim();
}
const mpq_class& Neighbor::getY(size_t i) const {
ASSERT(isValid());
ASSERT(i < getYDim());
if (isZero())
return _lat->getZero();
else
return _lat->getYMatrix()(_row, i);
}
size_t Neighbor::getYDim() const {
ASSERT(isValid());
return _lat->getYDim();
}
bool Neighbor::isZero() const {
ASSERT(isValid());
return _row == _lat->getNeighborCount() + 1;
}
bool Neighbor::isValid() const {
return _lat != 0;
}
size_t Plane::getTypeCount(size_t type) const {
map<size_t, size_t>::const_iterator it = typeCounts.find(type);
if (it == typeCounts.end())
return 0u;
else
return it->second;
}
size_t Plane::getMaxType() const {
if (typeCounts.empty())
return 0u;
else
return typeCounts.rbegin()->first;
}
NeighborPlace Plane::getPlace(Neighbor neighbor) const {
if (neighbor.isZero())
return InPlane;
ASSERT(neighbor.getRow() < neighborPlace.size());
return neighborPlace[neighbor.getRow()];
}
bool Plane::inPlane(Neighbor neighbor) const {
return getPlace(neighbor) == InPlane;
}
bool Plane::isPivot(const Mlfb& mlfb) const {
const size_t type = getType(mlfb);
return type == 1 || type == 3;
}
bool Plane::isSidePivot(const Mlfb& mlfb) const {
if (!isPivot(mlfb))
return false;
for (size_t i = 0; i < 4; ++i)
if (is22(*(mlfb.getEdge(i))))
return true;
return false;
}
bool Plane::is22(const Mlfb& mlfb) const {
const size_t type = getType(mlfb);
return type == 2;
}
bool Plane::isFlat(const Mlfb& mlfb) const {
return getType(mlfb) == 4;
}
bool Neighbor::isSpecial() const {
ASSERT(isValid());
for (size_t i = 1; i < _lat->getYDim(); ++i)
if (getY(i) <= 0)
return false;
return true;
}
bool Neighbor::isGenerator() const {
if (isZero())
return false;
else
return !_lat->isSum(*this);
}
size_t Plane::getType(const Mlfb& mlfb) const {
size_t type = 0;
for (size_t i = 0; i < mlfb.getPointCount(); ++i)
if (inPlane(mlfb.getPoint(i)))
++type;
return type;
}
string Neighbor::getName() const {
if (isZero())
return "zero";
if (!isValid())
return "none";
ostringstream name;
name << 'n' << (getRow() + 1);
if (isSpecial())
name << 's';
if (isGenerator())
name << 'g';
return name.str();
}
GrobLat::GrobLat(const Matrix& matrix, const SatBinomIdeal& ideal) {
_ideal = ideal;
_ideal.getMatrix(_y);
// transpose in preparation for solve
transpose(_y);
transpose(_mat, matrix);
solve(_h, _mat, _y);
// un-transpose
transpose(_mat);
transpose(_y);
transpose(_h);
_isSumRow.resize(getNeighborCount());
for (size_t i = 0; i < getNeighborCount(); ++i) {
for (size_t j = 0; j < i; ++j) {
Neighbor sum = getSum(getNeighbor(i), getNeighbor(j));
if (sum.isValid())
_isSumRow[sum.getRow()] = true;
}
}
_nonSums.clear();
for (size_t i = 0; i < _isSumRow.size(); ++i)
if (!_isSumRow[i])
_nonSums.push_back(getNeighbor(i));
}
swap
void swap(HilbertSlice &a, HilbertSlice &b)
Definition HilbertSlice.h:63
checkFlatSeq
void checkFlatSeq(const vector< SeqPos > &flatSeq, const GrobLat &lat, const Plane &plane)
Definition LatticeAlgs.cpp:1180
checkGraphOnPlane
void checkGraphOnPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:1017
checkPlanes
void checkPlanes(const vector< TriPlane > &thinPlanes, const vector< Plane > &dtPlanes)
Definition LatticeAlgs.cpp:99
checkMlfbs
void checkMlfbs(const vector< Mlfb > &mlfbs, const GrobLat &lat)
Definition LatticeAlgs.cpp:642
disjointSeqs
bool disjointSeqs(const vector< SeqPos > &a, const vector< SeqPos > &b)
Definition LatticeAlgs.cpp:738
computePivotSeqs
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.
Definition LatticeAlgs.cpp:881
prevInSeq
SeqPos prevInSeq(SeqPos pos)
Definition LatticeAlgs.cpp:441
checkSeqs
void checkSeqs(const vector< vector< SeqPos > > &left, const vector< vector< SeqPos > > &right, const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:966
checkSeq
void checkSeq(vector< bool > &seenOnSide, const vector< SeqPos > &seq, const Plane &plane)
Definition LatticeAlgs.cpp:904
computeSeqs
void computeSeqs(vector< vector< SeqPos > > &left, vector< vector< SeqPos > > &right, const vector< Mlfb > &mlfbs, const Plane &plane)
Definition LatticeAlgs.cpp:769
pivotToFlatFacet
size_t pivotToFlatFacet(const Mlfb &pivot, const Plane &plane)
Returns the facet to push in of pivot to get to a flat.
Definition LatticeAlgs.cpp:726
pushOutFacetZero
size_t pushOutFacetZero(const vector< mpz_class > &rhs, const GrobLat &lat)
Definition LatticeAlgs.cpp:170
getPlaceCode
char getPlaceCode(NeighborPlace place)
Definition LatticeAlgs.cpp:135
computeMlfbs
void computeMlfbs(vector< Mlfb > &mlfbs, const GrobLat &lat)
Definition LatticeAlgs.cpp:239
computePivots
void computePivots(vector< const Mlfb * > &pivots, const vector< Mlfb > &mlfbs, const Plane &plane, const vector< SeqPos > &flatSeq)
Put all pivots into pivots.
Definition LatticeAlgs.cpp:746
checkGraph
void checkGraph(const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:1060
checkMiddle
void checkMiddle(const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:985
getIndexSum
mpq_class getIndexSum(const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:1321
checkDoubleTriangle
void checkDoubleTriangle(const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:1043
computePlanes
void computePlanes(vector< Plane > &planes, const GrobLat &lat, vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:511
getEdgePos
const char * getEdgePos(size_t index)
Definition LatticeAlgs.cpp:1310
checkDoubleTrianglePlanes
void checkDoubleTrianglePlanes(const vector< Plane > &planes, const GrobLat &lat, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:649
checkSide
void checkSide(vector< bool > &pivotOnSide, const vector< vector< SeqPos > > &side, const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:937
computeFlatSeq
void computeFlatSeq(vector< SeqPos > &seq, const vector< Mlfb > &mlfbs, const Plane &plane)
Definition LatticeAlgs.cpp:471
check0Graph
void check0Graph(const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:611
nextInSeq
SeqPos nextInSeq(SeqPos pos)
Definition LatticeAlgs.cpp:415
pushOutFacetPositive
size_t pushOutFacetPositive(size_t facetPushOut, const vector< mpz_class > &rhs, const GrobLat &lat)
Definition LatticeAlgs.cpp:145
getThinPlanes
void getThinPlanes(vector< TriPlane > &planes, const GrobLat &lat)
Definition LatticeAlgs.cpp:22
checkNonSums
void checkNonSums(const GrobLat &lat)
Definition LatticeAlgs.cpp:1157
computeRhs
void computeRhs(vector< mpz_class > &rhs, const vector< Neighbor > points)
Definition LatticeAlgs.cpp:190
checkPlane
void checkPlane(const Plane &plane, const vector< Mlfb > &mlfbs)
Definition LatticeAlgs.cpp:714
checkPlaneTri
void checkPlaneTri(const GrobLat &lat, const vector< Mlfb > &mlfbs, const vector< const Mlfb * > &pivots, const Plane &plane)
Definition LatticeAlgs.cpp:1259
checkPivotSeqs
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.
Definition LatticeAlgs.cpp:1092
computeFlatIntervalCount
size_t computeFlatIntervalCount(const vector< SeqPos > &flatSeq)
Definition LatticeAlgs.cpp:445
LatticeAlgs.h
NeighborPlace
NeighborPlace
Definition LatticeAlgs.h:58
InPlane
@ InPlane
Definition LatticeAlgs.h:59
OverPlane
@ OverPlane
Definition LatticeAlgs.h:61
UnderPlane
@ UnderPlane
Definition LatticeAlgs.h:60
NoPlace
@ NoPlace
Definition LatticeAlgs.h:62
CHECK
#define CHECK(X)
Definition LatticeAlgs.h:48
solve
bool solve(Matrix &sol, const Matrix &lhs, const Matrix &rhs)
Sets sol to some matrix such that lhs*sol=rhs and returns true if such a matrix exists.
Definition Matrix.cpp:348
getParallelogramAreaSq
mpq_class getParallelogramAreaSq(const Matrix &mat)
Returns the square of the area of the parallelogram whose vertices are the 4 rows of mat.
Definition Matrix.cpp:456
product
void product(Matrix &prod, const Matrix &a, const Matrix &b)
Sets prod to a * b.
Definition Matrix.cpp:116
isParallelogram
bool isParallelogram(const Matrix &mat)
Returns true if the rows of mat are the (4) vertices of a parallelogram.
Definition Matrix.cpp:452
nullSpace
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
determinant
mpq_class determinant(const Matrix &mat)
Returns the determinant of mat.
Definition Matrix.cpp:410
copyRow
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
hasSameRowSpace
bool hasSameRowSpace(const Matrix &a, const Matrix &b)
Returns true if a and b have the same row space.
Definition Matrix.cpp:390
transpose
void transpose(Matrix &trans, const Matrix &mat)
Sets trans to the transpose of mat.
Definition Matrix.cpp:129
BigIdeal
Definition BigIdeal.h:27
BigIdeal::getVarCount
size_t getVarCount() const
Definition BigIdeal.h:148
BigIdeal::getGeneratorCount
size_t getGeneratorCount() const
Definition BigIdeal.h:144
BigTermRecorder
BigTermRecorder records all the terms it consumes into an ideal.
Definition BigTermRecorder.h:31
BigTermRecorder::empty
bool empty() const
Definition BigTermRecorder.cpp:70
BigTermRecorder::releaseIdeal
unique_ptr< BigIdeal > releaseIdeal()
Definition BigTermRecorder.cpp:74
GrobLat
A lattice with associated Grobner basis/neighbors.
Definition LatticeAlgs.h:124
GrobLat::_h
Matrix _h
Definition LatticeAlgs.h:229
GrobLat::getNeighbor
Neighbor getNeighbor(size_t row) const
Definition LatticeAlgs.h:128
GrobLat::getYDim
size_t getYDim() const
Definition LatticeAlgs.h:174
GrobLat::getHDim
size_t getHDim() const
Definition LatticeAlgs.h:179
GrobLat::_ideal
SatBinomIdeal _ideal
Definition LatticeAlgs.h:231
GrobLat::_isSumRow
vector< bool > _isSumRow
Definition LatticeAlgs.h:225
GrobLat::_y
Matrix _y
Definition LatticeAlgs.h:228
GrobLat::getMatrix
const Matrix & getMatrix() const
Definition LatticeAlgs.h:140
GrobLat::GrobLat
GrobLat(const Matrix &matrix, const SatBinomIdeal &ideal)
Definition LatticeAlgs.cpp:1533
GrobLat::getSum
Neighbor getSum(Neighbor a, Neighbor b) const
Definition LatticeAlgs.h:144
GrobLat::getNeighborCount
size_t getNeighborCount() const
Definition LatticeAlgs.h:133
GrobLat::getInitialIdeal
void getInitialIdeal(BigIdeal &ideal) const
Definition LatticeAlgs.h:187
GrobLat::_mat
Matrix _mat
Definition LatticeAlgs.h:230
GrobLat::getYMatrix
const Matrix & getYMatrix() const
Definition LatticeAlgs.h:138
GrobLat::getNonSums
const vector< Neighbor > & getNonSums() const
Definition LatticeAlgs.h:197
GrobLat::getHMatrix
const Matrix & getHMatrix() const
Definition LatticeAlgs.h:139
GrobLat::_nonSums
vector< Neighbor > _nonSums
Definition LatticeAlgs.h:226
GrobLat::isPointFreeBody
bool isPointFreeBody(Neighbor a, Neighbor b) const
Returns true if the smallest body containing zero, a and b has no neighbor in its interior.
Definition LatticeAlgs.h:202
Matrix
Definition Matrix.h:26
Matrix::getColCount
size_t getColCount() const
Definition Matrix.h:31
Matrix::getRowCount
size_t getRowCount() const
Definition Matrix.h:30
Mlfb
Definition LatticeAlgs.h:289
Mlfb::getRhs
const vector< mpz_class > & getRhs() const
Definition LatticeAlgs.h:319
Mlfb::getPoint
Neighbor getPoint(size_t offset) const
Definition LatticeAlgs.h:302
Mlfb::_rhs
vector< mpz_class > _rhs
Definition LatticeAlgs.h:377
Mlfb::_points
vector< Neighbor > _points
Definition LatticeAlgs.h:378
Mlfb::edges
vector< Mlfb * > edges
Definition LatticeAlgs.h:370
Mlfb::index
mpq_class index
Definition LatticeAlgs.h:368
Mlfb::edgeHitsFacet
vector< size_t > edgeHitsFacet
Definition LatticeAlgs.h:371
Mlfb::hasPoint
bool hasPoint(Neighbor n) const
Definition LatticeAlgs.h:295
Mlfb::_offset
size_t _offset
Definition LatticeAlgs.h:379
Mlfb::reset
void reset(size_t offset, const vector< Neighbor > &points)
Definition LatticeAlgs.cpp:202
Mlfb::dotDegree
mpz_class dotDegree
Definition LatticeAlgs.h:369
Mlfb::minInitialFacet
size_t minInitialFacet
Definition LatticeAlgs.h:372
Mlfb::getOffset
size_t getOffset() const
Definition LatticeAlgs.h:315
Mlfb::getPointCount
size_t getPointCount() const
Definition LatticeAlgs.h:307
Mlfb::isParallelogram
bool isParallelogram() const
Definition LatticeAlgs.h:364
Mlfb::getEdge
const Mlfb * getEdge(size_t indexParam) const
Definition LatticeAlgs.h:347
Mlfb::getMinInitialFacet
size_t getMinInitialFacet() const
Definition LatticeAlgs.h:291
Mlfb::getHitsFacet
size_t getHitsFacet(size_t indexParam) const
Definition LatticeAlgs.h:342
Mlfb::_isParallelogram
bool _isParallelogram
Definition LatticeAlgs.h:380
Neighbor
Definition LatticeAlgs.h:86
Neighbor::_lat
const GrobLat * _lat
Definition LatticeAlgs.h:119
Neighbor::Neighbor
Neighbor()
Definition LatticeAlgs.cpp:1398
Neighbor::isSpecial
bool isSpecial() const
Definition LatticeAlgs.cpp:1496
Neighbor::_row
size_t _row
Definition LatticeAlgs.h:120
Neighbor::isGenerator
bool isGenerator() const
Definition LatticeAlgs.cpp:1504
Neighbor::getH
const mpq_class & getH(size_t i) const
Definition LatticeAlgs.cpp:1410
Neighbor::getName
string getName() const
Definition LatticeAlgs.cpp:1519
Neighbor::getRow
size_t getRow() const
Definition LatticeAlgs.h:108
Neighbor::getYDim
size_t getYDim() const
Definition LatticeAlgs.cpp:1433
Neighbor::getY
const mpq_class & getY(size_t i) const
Definition LatticeAlgs.cpp:1424
Neighbor::getHDim
size_t getHDim() const
Definition LatticeAlgs.cpp:1419
Neighbor::isZero
bool isZero() const
Definition LatticeAlgs.cpp:1438
Neighbor::isValid
bool isValid() const
Definition LatticeAlgs.cpp:1443
Plane
Definition LatticeAlgs.h:258
Plane::isFlat
bool isFlat(const Mlfb &mlfb) const
Definition LatticeAlgs.cpp:1492
Plane::neighborPlace
vector< NeighborPlace > neighborPlace
Definition LatticeAlgs.h:284
Plane::flatSeq
vector< SeqPos > flatSeq
Definition LatticeAlgs.h:285
Plane::tri
Tri tri
Definition LatticeAlgs.h:279
Plane::isPivot
bool isPivot(const Mlfb &mlfb) const
Definition LatticeAlgs.cpp:1473
Plane::flatIntervalCount
size_t flatIntervalCount
Definition LatticeAlgs.h:281
Plane::nullSpaceBasis
Matrix nullSpaceBasis
Definition LatticeAlgs.h:278
Plane::is22
bool is22(const Mlfb &mlfb) const
Definition LatticeAlgs.cpp:1487
Plane::getPlace
NeighborPlace getPlace(Neighbor neighbor) const
Definition LatticeAlgs.cpp:1462
Plane::hasFlat
bool hasFlat() const
Definition LatticeAlgs.h:274
Plane::typeCounts
map< size_t, size_t > typeCounts
Definition LatticeAlgs.h:283
Plane::pivots
vector< const Mlfb * > pivots
Definition LatticeAlgs.h:286
Plane::getTypeCount
size_t getTypeCount(size_t type) const
Definition LatticeAlgs.cpp:1447
Plane::isSidePivot
bool isSidePivot(const Mlfb &mlfb) const
Definition LatticeAlgs.cpp:1478
Plane::rowAB
Matrix rowAB
Definition LatticeAlgs.h:280
Plane::getMaxType
size_t getMaxType() const
Definition LatticeAlgs.cpp:1455
Plane::inPlane
bool inPlane(Neighbor neighbor) const
Definition LatticeAlgs.cpp:1469
Plane::getType
size_t getType(const Mlfb &mlfb) const
Definition LatticeAlgs.cpp:1511
SatBinomIdeal
Represents a saturated binomial ideal.
Definition SatBinomIdeal.h:28
SliceFacade
A facade for operations on monomial ideals using the Slice Algorithm.
Definition SliceFacade.h:44
SliceFacade::computeIrreducibleDecomposition
void computeIrreducibleDecomposition(bool encode)
Compute the unique irredundant set of irreducible ideals whose intersection equals ideal.
Definition SliceFacade.cpp:109
SliceParams
Definition SliceParams.h:25
TriPlane
Definition LatticeAlgs.h:383
TriPlane::closeToPlane
bool closeToPlane(Neighbor a)
Definition LatticeAlgs.h:403
TriPlane::isParallel
bool isParallel(const TriPlane &plane) const
Definition LatticeAlgs.h:420
TriPlane::isLine
bool isLine() const
Definition LatticeAlgs.h:399
Tri
Definition LatticeAlgs.h:235
Tri::_a
Neighbor _a
Definition LatticeAlgs.h:249
Tri::_bSideMlfbs
vector< const Mlfb * > _bSideMlfbs
Definition LatticeAlgs.h:253
Tri::_b
Neighbor _b
Definition LatticeAlgs.h:250
Tri::getASideMlfbs
const vector< const Mlfb * > & getASideMlfbs() const
Definition LatticeAlgs.h:243
Tri::Tri
Tri(Neighbor a, Neighbor b, Neighbor sum, const vector< Mlfb > &mlfbs, const GrobLat &lat)
Definition LatticeAlgs.cpp:567
Tri::getB
Neighbor getB() const
Definition LatticeAlgs.h:241
Tri::_interior
vector< Neighbor > _interior
Definition LatticeAlgs.h:254
Tri::getBSideMlfbs
const vector< const Mlfb * > & getBSideMlfbs() const
Definition LatticeAlgs.h:244
Tri::getA
Neighbor getA() const
Definition LatticeAlgs.h:240
Tri::getSum
Neighbor getSum() const
Definition LatticeAlgs.h:242
Tri::_boundary
vector< Neighbor > _boundary
Definition LatticeAlgs.h:255
Tri::_sum
Neighbor _sum
Definition LatticeAlgs.h:251
Tri::_aSideMlfbs
vector< const Mlfb * > _aSideMlfbs
Definition LatticeAlgs.h:252
stdinc.h
This header file includes common definitions and is included as the first line of code in every imple...
ASSERT
#define ASSERT(X)
Definition stdinc.h:86
SeqPos
Definition LatticeAlgs.h:69
SeqPos::order
void order()
Definition LatticeAlgs.cpp:1369
SeqPos::getReverse
SeqPos getReverse() const
Definition LatticeAlgs.cpp:1357
SeqPos::mlfb
const Mlfb * mlfb
Definition LatticeAlgs.h:78
SeqPos::SeqPos
SeqPos()
Definition LatticeAlgs.cpp:1328
SeqPos::comingFromFacet
size_t comingFromFacet
Definition LatticeAlgs.h:81
SeqPos::fixFacet1
size_t fixFacet1
Definition LatticeAlgs.h:79
SeqPos::operator<
bool operator<(const SeqPos &pos) const
Definition LatticeAlgs.cpp:1374
SeqPos::fixFacet2
size_t fixFacet2
Definition LatticeAlgs.h:80
SeqPos::getForwardFacet
size_t getForwardFacet() const
Definition LatticeAlgs.cpp:1345
SeqPos::getBackFacet
size_t getBackFacet() const
Definition LatticeAlgs.cpp:1353