newPackage(
"Triangulations",
Version => "0.1",
Date => "13 Nov 2022",
Authors => {{
Name => "Mike Stillman",
Email => "mike@math.cornell.edu",
HomePage=>"http://www.math.cornell.edu/~mike"
}},
Headline => "triangulations of polyhedra and point sets",
Keywords => {"Combinatorics"},
PackageImports => {"FourierMotzkin"},
PackageExports => {
"Topcom",
"Polyhedra" -- really only needed for `regularSubdivision`?
},
DebuggingMode => false
)
export {
-- types defined here
"Triangulation",
"Chirotope",
"triangulation",
"vectors",
"regularFineTriangulation",
"chirotope",
"naiveChirotope",
"flips",
"bistellarFlip",
"neighbors",
"generateTriangulations",
"allTriangulations",
"isStar",
"isFine",
"isRegularTriangulation", -- note the non-use of isRegular. Is this ok?
"regularTriangulationWeights",
"affineCircuits",
"volumeVector",
"gkzVector",
"delaunayWeights",
"delaunaySubdivision",
"fineStarTriangulation",
"regularFineStarTriangulation",
"naiveIsTriangulation",
"ConeIndex"
}
augment = method()
augment Matrix := (A) -> (
-- A is a matrix over ZZ
-- add in a last row of 1's.
n := numColumns A;
ones := matrix {{n : 1}};
A || ones
)
-- also defined:
-- isRegularTriangulation T
-- regularTriangulationWeights T
-- naiveIsTriangulation T
-- TODO: design decision:
-- The "rays" or "vertices": stored as matrix or list?
-- Store list as vector configuration (i.e. homogenized?)
-- Accessors: max, rays?
Triangulation = new Type of HashTable
Chirotope = new Type of HashTable
net Triangulation := T -> "triangulation " | net max T
toString Triangulation := T -> toString max T
toExternalString Triangulation := T -> "triangulation(transpose matrix" | toExternalString vectors T | "," | toExternalString max T | ", Homogenize => false)"
--expression Triangulation := X -> (describe Triangulation)#0
--describe Triangulation := X -> describe max X
-- XXX working on this 15 June 2022.
-- Change all uses of fineStarTriangulation to return
-- Assumptions: A is a dxn matrix over the integers, or rationals.
-- convexHull A
-- all columns of A should be boundary points (with the possible exception of the last column).
-- A can either contain the interior point to cone over, or not.
-- If it does contain it, I believe that no element of tri should use that index.
-- Returns: a new triangulation of the point set A | 0.
-- NOTE: not yet! it currently returns a list of subsets, such that if one appends to each subset
-- the number "numcols A", this will be a star triangulation.
-- NOTE: I'm not sure why this would need to be regular either...
fineStarTriangulation = method(Options => {ConeIndex => null})
fineStarTriangulation(Matrix, List) := List => opts -> (A, tri) -> (
coneindex := if opts.ConeIndex === null then
numcols A -- in this case, assume that the cone point is a new index, which would be 'numcols A'
else if instance(opts.ConeIndex, ZZ) and opts.ConeIndex >= 0 and opts.ConeIndex <= numcols A + 1 then
opts.ConeIndex
else
error "ConeIndex must be an index from 0 .. number of columns + 1...";
aA := augment A;
-- H := first halfspaces convexHull aA;
H := transpose(-(first fourierMotzkin aA));
myfacets := for e in entries H list (
positions(flatten entries(matrix {e} * aA), x -> x == 0)
);
sort unique flatten for f in tri list for g in myfacets list (
a := toList(set g * set f);
if #a < numRows A then
continue
else
sort append(a, coneindex)
)
-- newtri = for f in newtri list append(f, numColumns A)
)
-- TODO: is this really a regular triangulation?
-- I think it might be, as long as the cone index is not used in regularFineTriangulation.
regularFineStarTriangulation = method(Options => options fineStarTriangulation)
regularFineStarTriangulation Matrix := List => opts -> (A) -> fineStarTriangulation(A, topcomRegularFineTriangulation A, opts)
-- TODO: I am not sure that this is correct.
naiveIsTriangulation = method()
naiveIsTriangulation(Matrix, List, List) := (A, circuits, tri) -> (
aA := augment A;
-- H := first halfspaces convexHull aA;
H := transpose(-(first fourierMotzkin aA));
myfacets := for e in entries H list (
positions(flatten entries(matrix {e} * aA), x -> x == 0)
);
-- test 1: each wall should be in a facet of the convex hull, or occur exactly twice.
-- This is NOT correct?!
walls := tally flatten for t in tri list subsets(t,#t-1);
test1 := for k in keys walls list (
if any(myfacets, f -> isSubset(k,f)) then
walls#k == 1
else
walls#k == 2
);
if any(test1, x -> not x) then return false;
-- test 2: for each oriented circuit Z = (Z+, Z-)
test2 := for z in circuits list (
# select(tri, t -> isSubset(z_0, t)),
# select(tri, t -> isSubset(z_1, t))
);
all(test2, x -> x#0 == 0 or x#1 == 0)
)
naiveIsTriangulation(Matrix, List) := (A, tri) -> naiveIsTriangulation(A, orientedCircuits A, tri)
-- Allow both rays and points, i.e. homogenized A or not
triangulation = method(Options => {Homogenize => true}) -- true means this is a point set.
triangulation(Matrix, List) := Triangulation => opts -> (A, tri) -> (
-- should we check basic things? e.g. tri is a list of
-- lists of integers in the range 0..numcols A - 1.
-- and that A is a matrix over ZZ or QQ?
if ring A =!= ZZ and ring A =!= QQ then
error "expected matrix over ZZ or QQ";
n := numcols A;
if not all(tri, f -> all(f, p -> instance(p, ZZ))) then
error "expected a list of list of integers";
-- should we sort the sets? Probably...
A1 := if opts.Homogenize then augment A else A;
vecs := transpose entries A1;
sorted := tri//sort/sort; -- this sorts the triangulation
T := new Triangulation from {
cache => new CacheTable,
symbol vectors => vecs,
symbol max => sorted
};
T.cache#"point set" = opts.Homogenize;
T.cache.matrix = A1;
T
)
vectors = method()
vectors Triangulation := List => T -> T.vectors
max Triangulation := List => T -> T.max
matrix Triangulation := opts -> T -> T.cache.matrix
isWellDefined Triangulation := Boolean => T -> (
topcomIsTriangulation(matrix T, max T, Homogenize => false)
)
Triangulation == Triangulation := Boolean => (S, T) -> S === T
naiveIsTriangulation Triangulation := Boolean => T -> (
naiveIsTriangulation(matrix T, max T) -- BUG: needs to take Homogenize?
)
-- The following is currently only for point sets
-- and is probably slow for larger triangulations too.
isTriangulation = method()
isTriangulation(Matrix, List) := (M, tri) -> (
-- for the moment, we assume that M is a point set.
d := #(tri#0) - 1;
P := convexHull M;
M' := M || matrix{{numcols M: 1}};
volP := d! * volume P;
volP2 := sum for t in tri list abs det M'_t;
if volP != volP2 then (
<< "volume is not correct: " << volP << " != " << volP2 << endl;
return false;
);
simplices := hashTable for t in tri list t => convexHull M_t;
for x in subsets(keys simplices, 2) do
if dim(intersection(simplices#(x#0), simplices#(x#1))) == d then (
<< "simplices " << x << " intersect in full dimensional region" << endl;
return false;
);
true
)
isRegularTriangulation = method(Options => {Homogenize => true})
isRegularTriangulation Triangulation := Boolean => opts -> T -> (
topcomIsRegularTriangulation(T.cache.matrix, max T, Homogenize => false)
)
isRegularTriangulation(Matrix, List) := Boolean => opts -> (A, tri) -> (
topcomIsRegularTriangulation(A, tri, opts)
)
regularTriangulationWeights = method(Options => options isRegularTriangulation)
regularTriangulationWeights Triangulation := List => opts -> T -> (
topcomRegularTriangulationWeights(matrix T, max T, Homogenize => false)
)
regularTriangulationWeights(Matrix, List) := List => opts -> (A, tri) -> (
topcomRegularTriangulationWeights(A, tri, opts)
)
regularFineTriangulation = method(Options => options isRegularTriangulation)
regularFineTriangulation Matrix := Triangulation => opts -> (A) -> (
tri := topcomRegularFineTriangulation(A, opts);
if tri === null then null else
triangulation(A, tri, opts)
)
-- TODO/BUG: this ASSUMES (A, tri) is a triangulation.
-- TODO: add in Homogenize as an option?
isFine = method()
isFine(Matrix, List) := Boolean => (A, tri) -> (
numcols A == tri//flatten//unique//length
)
-- TODO/BUG: this ASSUMES (A, tri) is a triangulation.
isStar = method()
isStar(Matrix, List) := Boolean => (A, tri) -> (
-- assumption? last column of A is the zero element? (or the one the star is taken with respect to).
origin := numcols A - 1;
all(tri, t -> member(origin, t))
)
isStar Triangulation := Boolean => T -> isStar(matrix T, max T)
isFine Triangulation := Boolean => T -> isFine(matrix T, max T)
allTriangulations = method(Options => options topcomAllTriangulations)
allTriangulations Matrix := List => opts -> A -> (
tris := topcomAllTriangulations(A, opts);
for t in tris list triangulation(A, t, Homogenize => opts.Homogenize)
)
delaunayWeights = method()
delaunayWeights Matrix := (A) -> (
matrix{for i from 0 to numcols A - 1 list ((transpose A_{i}) * A_{i})_(0,0)}
)
delaunaySubdivision = method()
delaunaySubdivision Matrix := A -> elapsedTime regularSubdivision(A, elapsedTime delaunayWeights A)
-----------------------------------------------------------
-- Chirotope code. This could potentially go elsewhere? --
-----------------------------------------------------------
chirotope = method(Options => true)
chirotope String := {} >> opts -> s -> (
new Chirotope from {
cache => new CacheTable,
symbol String => s
}
)
toString Chirotope := String => OM -> OM.String
chirotope Matrix := Chirotope => {Homogenize => true} >> opts -> A -> (
chirotope chirotopeString(A, opts)
)
naiveChirotope = method(Options => true)
naiveChirotope Matrix := Chirotope => {Homogenize => true} >> opts -> A -> (
chirotope naiveChirotopeString(A, opts)
)
Chirotope == Chirotope := Boolean => (C, D) -> toString C === toString D
--------------------------------------
-- Flips and triangulations ----------
--------------------------------------
-- This code is not an interface to topcom.
-- This is my own code, which used to be
-- part of `StringTorics`, but seems more
-- appropriate in a `Triangulations` package.
-- It is much slower than topcom, but allows to collect a certain number of triangulations.
-- TODO: make an interface that yields the next triangulation when called? (matching python usage).
--------------------------------------
sortTriangulation = method()
sortTriangulation List := (T) -> sort for t in T list sort t
-- TODO: why codim2 and codim2s?
codim2 = (tri) -> (
-- find all pairs whose intersection is all but one of the vertices.
n := #tri_0;
select(subsets(tri, 2), v -> length toList ((set (v#0)) * (set (v#1))) == n-1)
)
codim2s = (tri) -> (
n := #tri_0;
C2 := unique apply(subsets(tri, 2), v -> sort toList (set v#0 + set v#1));
select(C2, x -> #x == n+1)
)
affineCircuits = method()
-- TODO: RENAME to flips, or possibleFlips
-- previously in triangulations-code.m2, named affineCircuits
affineCircuits(Matrix, List) := (Amat, tri) -> (
-- tri: a set of subsets of size d from 0..#columns(Amat)-1, that form
-- a triangulation of conv(columns of Amat).
-- Amat: d-1 x n matrix over ZZ or QQ of the points, with the d-1 rows independent.
-- OR of size d x n, with the last row all 1's. (or the rows are rank (d-1)).
d := #(tri#0);
assert all(tri, t -> #t == d);
if numrows Amat == d-1 then
Amat = Amat || matrix{ toList(numcols Amat : 1) };
c2 := codim2s tri;
unique for c in c2 list (
z := flatten entries syz Amat_c;
rsort {c_(positions(z, zi -> zi > 0)), c_(positions(z, zi -> zi < 0))}
)
)
affineCircuits Triangulation := T -> affineCircuits(matrix T, max T)
-- previously in triangulations-code.m2,
link = method()
link(List, List) := (tau, triangulation) -> (
S := select(triangulation, t -> isSubset(tau, t));
sort for s in S list sort toList (set s - set tau)
)
flips = method(Options => options topcomFlips)
flips Triangulation := List => opts -> T -> (
first topcomFlips(matrix T, max T, Homogenize => false, RegularOnly => opts.RegularOnly) -- TODO: "first" here is a hack, until topcom version is set.
)
bistellarFlip = method()
-- previously in triangulations-code.m2, named flip.
bistellarFlip(List,List) := (tri, affineCircuit) -> (
-- returns a new triangulation
-- input: tri: a list of lists of integers (a triangulation).
-- affineCircuit: a pair of lists of integers. Each is disjoint, the negative/positive part of the affine circuit.
-- a. start with an affine circuit
-- 1. compute 2 choices of triangulation for the circuit
whole := sort unique flatten affineCircuit;
T1 := for i in affineCircuit#0 list sort toList (set whole - set {i});
T2 := for i in affineCircuit#1 list sort toList (set whole - set {i});
-- 2. which is contained in tri?
-- XXXXXXXXXXXX
S1 := unique for tau in T1 list link(tau, tri);
S2 := unique for tau in T2 list link(tau, tri);
if #S1 > 1 or #S2 > 1 then (
--<< "link is not unique, here they are: S1: " << S1 << " S2: " << S2 << endl;
return null;
);
S1 = first S1;
S2 = first S2;
if #S1 > 0 and #S2 > 0 then (
error "my logic is wrong: somehow both links exist...";
return null;
);
if #S1 == 0 then (
S1 = S2;
(T1,T2) = (T2,T1);
);
-- Now we can change the triangulation:
T := set tri;
outgoing := set flatten for s in S1 list for tau in T1 list sort join(s,tau);
incoming := set flatten for s in S1 list for tau in T2 list sort join(s,tau);
--<< "outgoing: " << sort toList outgoing << endl;
--<< "incoming: " << sort toList incoming << endl;
sortTriangulation toList(T - outgoing + incoming)
)
bistellarFlip(Triangulation, List) := (T, affineCircuit) -> (
tri := bistellarFlip(max T, affineCircuit);
if tri === null then null else
triangulation(matrix T, tri, Homogenize => false)
)
neighbors = method()
neighbors Triangulation := List => T -> (
-- each element of the result is of the form:
-- {circuit, triangulation}
-- only circuits that do not change the support of the triangulation are included.
circuits0 := affineCircuits T;
circuits := select(circuits0, z -> #z#0 > 1 and #z#1 > 1);
for c in circuits list (
T1 := bistellarFlip(T, c);
if T1 === null then continue else {c, T1}
)
)
generateTriangulations = method(Options => {Limit=>infinity, RegularOnly=>false, Homogenize => true})
generateTriangulations Triangulation := opts -> T -> (
Amat := matrix T;
allT := new MutableHashTable;
allT#T = true;
TODO := {T};
while #TODO > 0 and #(keys allT) < opts.Limit do (
nextTRI := TODO#0;
TODO = drop(TODO,1);
flips := select(affineCircuits nextTRI, z -> #z#0 > 1 and #z#1 > 1);
fliptris := for f in flips list bistellarFlip(nextTRI, f);
newT := select(fliptris, x -> x =!= null);
for T in newT do (
if not allT#?T then (
--<< "new triangulation: " << T << endl;
if not opts.RegularOnly or isRegularTriangulation T then (
allT#T = true;
TODO = append(TODO, T);
);
));
if debugLevel > 0 then
<< "todo = " << #TODO << " and #triang = " << #(keys allT) << endl;
);
keys allT
)
-- TODO? need Homogenize?
generateTriangulations Matrix := opts -> Amat -> (
generateTriangulations(regularFineTriangulation Amat, opts)
)
generateTriangulations(Matrix, List) := opts -> (Amat, triang) -> (
tris := generateTriangulations(triangulation(Amat, triang, Homogenize => opts.Homogenize),
Limit => opts.Limit, RegularOnly => opts.RegularOnly);
tris/max -- strip the triangulation class from these. This matches (up to permutation) output of allTriangulations
)
volumeVector = method()
volumeVector (Matrix, List) := (Amat, tri) -> (
if #tri == 0 then error "expected at least one simplex";
nelems := #tri#0;
d := nelems-1;
if not all(tri, f -> #f == nelems)
then error "expected a triangulation";
if numrows Amat =!= nelems then Amat = Amat || matrix{{(numcols Amat):1}};
if numrows Amat =!= nelems then error "triangulation not compatible with matrix";
H := hashTable for t in tri list t => (abs det Amat_t)/d!;
if debugLevel > 0 then
<< "Volume = " << (sum values H)/d! << endl;
for i from 0 to numColumns Amat - 1 list (
T := select(tri, t -> member(i,t));
sum for t in T list H#t
)
)
volumeVector Triangulation := T -> volumeVector(matrix T, max T)
gkzVector = volumeVector
beginDocumentation()
-*
needsPackage "StringTorics"
topes = kreuzerSkarke(5, Limit => 10)
Q = cyPolytopeData topes_7
isFavorable Q
rays Q
A = matrix topes_7
P2 = polar convexHull A
Amat = latticePointList P2
-- one from h11=7
--LP = {{-1, -1, -1, 1}, {-1, -1, -1, 2}, {-1, -1, 0, 1}, {-1, 0, -1, 1}, {-1, 0, 1, 0}, {-1, 0, 2, 0}, {-1, 1, 0, 0}, {-1, 2, 0, -1}, {0, -1, -1, 1}, {0, 1, 1, -1}, {2, 0, 0, -1}, {0,0,0,0}}
*-
doc ///
Key
Triangulations
Headline
generating and manipulating triangulations of point or vector configurations
Description
Text
{\bf Warning!} This package is experimental, documentation is missing,
and the interface will be cleaned up and changed. Use only if these issues
don't bother you!
Text
@SUBSECTION "Data of a triangulation"@
Text
@UL {
TO (max, Triangulation),
TO (vectors, Triangulation)
}@
Text
@SUBSECTION "Creating triangulations"@
Text
@UL {
TO (triangulation, Matrix, List),
TO (regularFineTriangulation, Matrix),
TO (generateTriangulations, Triangulation),
TO (allTriangulations, Matrix)
}@
Text
@SUBSECTION "Properties of triangulations"@
Text
@UL {
TO (isWellDefined, Triangulation),
TO (isRegularTriangulation, Triangulation),
TO (regularTriangulationWeights, Triangulation),
TO (isStar, Triangulation),
TO (isFine, Triangulation),
TO (naiveIsTriangulation, Triangulation)
}@
Text
@SUBSECTION "Exploring the set of triangulations"@
Text
@UL {
TO (bistellarFlip, Triangulation, List),
TO (flips, Triangulation),
TO (neighbors, Triangulation),
TO (affineCircuits, Triangulation),
TO (volumeVector, Triangulation),
TO (gkzVector, Triangulation),
TO (delaunayWeights, Matrix),
TO (delaunaySubdivision, Matrix)
}@
Text
This package is designed to help compute and explore the set of all (or many) triangulations
of a point set or polytope.
We give a sample use of this package.
Example
LP = {{-1, 0, -1, 1}, {-1, 0, 1, 0}, {-1, 0, 2, -1}, {-1, 1, -1, 0}, {1, 0, -1, 0}, {1, 0, 1, 0}, {2, -1, -1, 0}, {0, 0, 1, 0}, {1, 0, 0, 0}, {0,0,0,0}}
A = transpose matrix LP
elapsedTime Ts = allTriangulations(A, Fine => true);
select(Ts, T -> isStar T)
#oo == 1
#Ts == 51
T = regularFineTriangulation A
elapsedTime Ts2 = generateTriangulations T;
#Ts2 == #Ts
SeeAlso
"Polyhedra::Polyhedra"
"Topcom::Topcom"
"ReflexivePolytopesDB::ReflexivePolytopesDB"
///
doc ///
Key
isRegularTriangulation
(isRegularTriangulation, Triangulation)
Headline
determine if a given triangulation is a regular triangulation
Usage
isRegularTriangulation T
Inputs
T:Triangulation
A triangulation of a point or vector configuration
Homogenize => Boolean
unused for this method
Outputs
:Boolean
whether the given triangulation is regular
Description
Text
A triangulation is called regular if it can be constructed in the following way: place the
point set in one higher dimension at various heights in the new variable. Compute the
convex hull. Collect the list of facets with downward pointing normal (last coordinate of normal vector
is negative). If each of these is a simplex, then these form a triangulation of the
original point set. A triangulation which arises this way is called {\it regular}. See
the book [deLoera et al] for more details and many beautiful properties of such triangulations.
Text
The following example is one of the simplest examples of a non-regular
triangulation. Notice that {\tt tri} is a triangulation of the
polytope which is the convex hull of the columns of $A$, which are
the only points allowed in the triangulation.
Example
A = transpose matrix {{0,3},{0,1},{-1,-1},{1,-1},{-4,-2},{4,-2}}
tri = {{0,1,2}, {1,3,5}, {2,3,4}, {0,1,5},
{0,2,4}, {3,4,5}, {1,2,3}}
T = triangulation(A, tri)
Text
We check that {\tt T} is indeed a triangulation, and whether it is a regular triangulation.
Example
isWellDefined T
isRegularTriangulation T
Text
Many of the functions in this package are wrappers for topcom functions.
Setting the global variable {\tt debugLevel} to either 1,2, or 5 will give more detail about
what files are written to Topcom, and what the executable is.
Setting {\tt debugLevel} to 0 means that the function will run silently.
Caveat
Does topcom check that the triangulation is actually well defined? I'm not sure... This is why we call
@TO (isWellDefined, Triangulation)@ first.
SeeAlso
(regularTriangulationWeights, Triangulation)
regularFineTriangulation
(isWellDefined, Triangulation)
///
///
Key
neighbors
(neighbors, Triangulation)
Headline
find the neighbors of a triangulation of a point or vector set
Usage
neighbors T
XXX
Inputs
T:Triangulation
Outputs
:List
of lists of size 2: the first element is a circuit, and the second is @ofClass Triangulation@
Description
Text
The neighbors of a triangulation $T$ are those triangulations which are related to $T$
by a bistellar flip. Some flips can change the set of vertices (the support) of a triangulation.
This function {\bf does not} consider these. Thus, if the triangulation is fine, the neighbors
returned from this function are all fine triangulations too.
Example
TODOwriteMe
Caveat
SeeAlso
///
doc ///
Key
triangulation
Headline
make a Triangulation object
Usage
triangulation(A, T)
Inputs
A:Matrix
T:List
representing a triangulation of the columns of $A$ (each element in the list
is a list of indices in the range $0, \ldots, n-1$, where $n$ is the number of
columns of $A$)
Homogenize => Boolean
if true, add a row of ones to think of this as a vector configuration in one higher dimension.
Outputs
:Triangulation
A @TO Triangulation@ object. Very little computation is performed. The matrix and list representing
a triangulation is packaged into an object to make clear that it is a triangulation
Description
Text
Example
P = hypercube 3
A = vertices P
T = topcomRegularFineTriangulation A
tri = triangulation(A, T)
matrix tri
vectors tri
max tri
isWellDefined tri
netList affineCircuits tri
isFine tri
isStar tri
isRegularTriangulation tri
Caveat
SeeAlso
(max, Triangulation)
(vectors, Triangulation)
regularFineTriangulation
(isWellDefined, Triangulation)
///
doc ///
Key
(allTriangulations, Matrix)
allTriangulations
[allTriangulations, ConnectedToRegular]
[allTriangulations, Fine]
[allTriangulations, Homogenize]
[allTriangulations, RegularOnly]
Headline
use topcom to generate all triangulations of a point or vector configuration
Usage
allTriangulations A
allTriangulations(A, Homogenize => true, Fine => true, RegularOnly => true)
Inputs
A:Matrix
Homogenize => Boolean
ConnectedToRegular => Boolean
Fine => Boolean
RegularOnly => Boolean
Outputs
:List
Description
Text
This function constructs all triangulations of the point set corresponding to $A$
(or triangulation of the cone over $A$, if {\tt Homogenize => false} is given.
With no optional arguments, the default is to construct all regular triangulations.
This function is a wrapper for the topcom function @TO "Topcom::topcomAllTriangulations"@,
and has the same optional arguments as that function.
This function returns a list of @TO Triangulation@'s.
A triangulation is a list of list of the indices of the maximal simplices in the triangulation.
(the index of the point corresponding to the $i$-th column (starting at $0$) is simply $i$.
For example, the following point set is the smallest which has a non-regular triangulation.
Example
A = transpose matrix {{0,3},{0,1},{-1,-1},{1,-1},{-4,-2},{4,-2}}
Ts = allTriangulations A
#Ts == 16
netList Ts
tri = Ts#0
isWellDefined tri
isRegularTriangulation tri
Ts/isRegularTriangulation
regularTriangulationWeights tri
Text
The following code determines the support of each triangulation, and tallies them.
Thus for example, we see that there are 6 regular fine triangulations.
Example
tally for tri in Ts list sort unique flatten max tri
Text
The method that topcom uses depends on the optional arguments {\tt Fine}, {\tt ConnectedToRegular}
and {\tt RegularOnly}.
Example
options allTriangulations
Text
If the optional argument {\tt Fine} is set to true, then only
{\it fine} triangulations (i.e.
those that involve every column of $A$) will be generated.
Example
Ts = allTriangulations(A, Fine => true);
#Ts == 6
Text
If the optional argument {\tt RegularOnly} is set to false, but
{\tt ConnectedToRegular} is true, it will generally take less time,
as the program doesn't need to check each triangulation to see if it is regular.
Example
T1s = allTriangulations(A, RegularOnly => true)
T2s = allTriangulations(A, RegularOnly => false)
#T1s
#T2s
Text
The following search would also yield all triangulations, even those not connected via
bistellar flips to regular triangulations.
Example
T3s = allTriangulations(A, RegularOnly => false, ConnectedToRegular => false)
#T3s
Text
Given the list of triangulations, we can query them using other topcom functions.
See also @TO "Triangulations::Triangulations"@ for other functionality.
Example
netList Ts
for tri in Ts list isWellDefined tri
for tri in Ts list isRegularTriangulation tri
for tri in Ts list regularTriangulationWeights tri
SeeAlso
(topcomNumTriangulations, Matrix)
generateTriangulations
///
doc ///
Key
generateTriangulations
(generateTriangulations, Matrix)
[generateTriangulations, Limit]
[generateTriangulations, RegularOnly]
[generateTriangulations, Homogenize]
Headline
generate all triangulations with certain properties
Usage
Ts = generateTriangulations A
Ts = generateTriangulations(A, tri)
Ts = generateTriangulations T
Inputs
A:Matrix
over the integers (or rationals?), whose columns are the
points to use
tri:List
A list of lists of indices representing a triangulation of
the columns of $A$
T:Triangulation
packages both $A$ and {\tt tri} in one object
Limit => ZZ
Stop after constructing this many triangulations
RegularOnly => Boolean
Only generate regular triangulations
Homogenize => Boolean
set to false in the case the columns form a vector configuration.
The default case is that the columns of $A$ are considered a point
configuration (not used in the variant {\tt generateTriangulations T})
Outputs
Ts:List
of lists of integers, each such list represents a triangulation
Description
Text
This function can be used to generate a set of triangulations of a point set
or vector configuration $A$ (the points are the columns of $A$).
It operates by starting with one triangulation (tri or T), if one is given, and if
not, it constructs a fine triangulation of the set of columns of $A$.
After this, it uses bistellar flips to generate neighbors, and continues, until
the limit is reached, or no new ones can be constructed.
Important note! This function generally starts with a fine triangulation (i.e. one
using all of the points in $A$), and only considers bistellar flips that give fine
triangulations.
Example
A = vertices hypercube 3
T = topcomRegularFineTriangulation A
tri = regularFineTriangulation A
Ts1 = generateTriangulations A -- list of Triangulation's.
Ts2 = generateTriangulations(A, T) -- list of list of subsets
Ts3 = generateTriangulations triangulation(A, T) -- list of Triangulations
Ts4 = generateTriangulations tri -- list of Triangulations
all(Ts4, isFine)
all(Ts4, isStar)
all(Ts4, isRegularTriangulation)
Ts4/isStar//tally
Ts4/gkzVector
volume convexHull A -- 8
stars1 = select(Ts4, t -> (gkzVector t)#-1 == 8)
stars2 = select(Ts4, isStar)
stars1 == stars2
Caveat
This function is written in the top level Macaulay2 language, and so is much slower
than @TO "allTriangulations"@, which calls the topcom code written in C++. On the other hand,
one can give this function a limit for the number of triangulations to generate,
so can be used to generate triangulations in the case when the number is too large to
write down all of them.
SeeAlso
"allTriangulations"
///
TEST ///
-- of homogenization and need for it.
-*
restart
debug needsPackage "Triangulations"
*-
-- test of isRegularTriangulation
A = transpose matrix {{-1,-1},{-1,1},{1,-1},{1,1},{0,0}}
regularFineTriangulation(A, Homogenize=>false) -- returns null. What does this mean?
T = regularFineTriangulation(A, Homogenize=>true) -- this is good.
needsPackage "Polyhedra"
assert(volume convexHull A == 4) -- actual volume
assert(sum (for t in max T list volume convexHull A_t) == 4)
A1 = transpose matrix {{-1,-1,1},{-1,1,1},{1,-1,1},{1,1,1},{0,0,1}}
T1 = regularFineTriangulation(A1, Homogenize=>false)
assert(T1 === T)
assert(max T === {{0, 1, 4}, {0, 2, 4}, {1, 3, 4}, {2, 3, 4}}) -- they are sorted, so it should be this.
assert isRegularTriangulation T
regularTriangulationWeights T == {1,1,0,0,0}
///
TEST ///
-*
restart
needsPackage "Triangulations"
*-
A = transpose matrix {{-1,-1},{-1,1},{1,-1},{1,1},{0,0}}
T = regularFineTriangulation A
naiveIsTriangulation T -- TODO: doc this, and allow A to be homogenized? Same with topcomIsTriangulation
-- XX
flips T -- not really functional. Returns internal stuff. I don't understand the format...
orientedCircuits A
orientedCocircuits A
chirotope A
assert(naiveChirotopeString A === chirotopeString A)
topcomNumTriangulations(A, RegularOnly => false, ConnectedToRegular => false) -- this should really be the default?
allTriangulations A
///
TEST ///
-- TODO: this is a test for Topcom, it seems?
-*
restart
needsPackage "Triangulations"
*-
-- test of isRegularTriangulation
A = transpose matrix {{0,3},{0,1},{-1,-1},{1,-1},{-4,-2},{4,-2}}
tri = {{0,1,2}, {1,3,5}, {2,3,4},
{0,1,5}, {0,2,4}, {3,4,5},
{1,2,3}}
assert not isRegularTriangulation(A,tri)
assert(null === regularTriangulationWeights(A,tri))
topcomNumTriangulations A
allTriangulations A
allTriangulations(A, Fine=>true)
allTriangulations(A, Fine=>true, RegularOnly=>false)
A = transpose matrix {{0,3},{0,1},{-1,-1},{1,-1},{-4,-2},{7,-2}}
tri = {{0,1,2}, {1,3,5}, {2,3,4},
{0,1,5}, {0,2,4}, {3,4,5},
{1,2,3}}
assert isRegularTriangulation(A,tri)
regularTriangulationWeights(A,tri) -- Question: how to test that this is correct
-- TODO: need a function which takes a point set, weights, and creates the lift (easy)
-- compute the lower hull of this polytope.
assert(toString chirotope A === naiveChirotopeString A)
orientedCircuits A
orientedCocircuits A
A = transpose matrix {{1,0},{0,1}}
tri = {{0,1}}
assert isRegularTriangulation(A,tri)
regularTriangulationWeights(A,tri) == {0,1} -- TODO: check that this is the correct answer
A = transpose matrix {{0}}
tri = {{0}}
assert isRegularTriangulation(A,tri)
regularTriangulationWeights(A,tri) == {1}
///
///
-- TODO: This test needs to be made to assert correct statements
-- How to test that triangulations are correct? What I thought worked does not.
needsPackage "Triangulations"
A = transpose matrix{{-1,-1},{-1,1},{1,-1},{1,1},{0,0}}
regularFineTriangulation A
tri = {{0, 2, 4}, {2, 3, 4}, {0, 1, 4}, {1, 3, 4}}
tri = {{0, 2, 4}, {2, 3, 4}, {0, 1, 4}, {1, 2, 3}}
isRegularTriangulation(A, tri) -- Wrong!!
badtri = {{0, 2, 4}, {2, 3, 4}, {0, 1, 4}, {1, 2, 3}}
debugLevel = 6
isRegularTriangulation(A,badtri) -- this should fail! But it doesn't seem to do so. BUG in something!!!
debugLevel = 0
-- hmmm, we can make non-sensical triangulations, without it noticing.
-- this should be a bug?
A = transpose matrix {{0,0},{0,1},{1,0},{1,1}}
tri = {{0,1,2},{0,2,3}}
assert isRegularTriangulation(A,tri)
tri = {{0,1,2},{1,2,3}}
assert isRegularTriangulation(A,tri)
///
TEST ///
needsPackage "Triangulations"
needsPackage "Polyhedra"
A = transpose matrix {{-1,-1,2},{-1,0,1},{-1,1,1},{0,-1,2},{0,1,1},{1,-1,3},{1,0,-1},{1,1,-2}}
debugLevel = 0
T = regularFineTriangulation A
assert isRegularTriangulation T
assert(regularTriangulationWeights T =!= null)
A = transpose matrix {{-1, 0, -1, -1}, {-1, 0, 0, -1}, {-1, 1, 2, -1}, {-1, 1, 2, 0}, {1, -1, -1, -1}, {1, -1, -1, 1}, {1, 0, -1, 2}, {1, 0, 1, 2}}
P2 = polar convexHull A
C = matrix {latticePoints P2}
tri = regularFineTriangulation C
assert isRegularTriangulation tri
regularTriangulationWeights tri -- is this correct? Some weights have negative values??
///
TEST ///
-- simple example of chirotope
-*
restart
needsPackage "Triangulations"
*-
A = transpose matrix {{-1,-1},{-1,1},{1,-1},{1,1},{0,0}}
tri = {{0, 2, 4}, {2, 3, 4}, {0, 1, 4}, {1, 3, 4}}
ch1 = chirotope A
ch2 = naiveChirotope A
assert(ch1 === ch2)
///
TEST ///
-- Bad triangulations of the square
V = transpose matrix {{0,0},{1,0},{0,1},{1,1}}
T1 = {{0,1,2}}
T2 = {{0,1,2},{0,1,3}}
T3 = {{0,1,2,3}}
assert(not naiveIsTriangulation(V, T1))
assert(not naiveIsTriangulation(V, T2))
assert(not naiveIsTriangulation(V, T3))
-- assert(not topcomIsTriangulation(V, T1)) -- topcom signals an error here
-- assert(not topcomIsTriangulation(V, T2)) -- topcom signals an error here
assert(not topcomIsTriangulation(V, T3))
///
-- This example is a good one, but takes too long to be run automatically
-- Actually, this is a test of Topcom....
///
restart
needsPackage "Topcom"
needsPackage "Polyhedra"
pts = {{-1,0,0,-1},{-1,0,1,-1},{-1,0,1,0},{-1,1,0,-1},{-1,1,0,0},{-1,1,1,2},{1,-1,0,-1},{1,0,-1,1},{1,-1,-1,-1},{0,0,0,-1}}
A = transpose matrix pts
debugLevel = 7
elapsedTime n1 = numTriangulations(A, Fine=>true, ConnectedToRegular=>true) -- 6.9 sec, 408 of these CORRECT
elapsedTime n2 = numTriangulations(A, Fine=>true, ConnectedToRegular=>false, RegularOnly => false) -- 116 sec, 448 of these
elapsedTime n3 = numTriangulations(A, Fine=>false, ConnectedToRegular=>true) -- 8 sec, 520 of these CORRECT
elapsedTime n4 = numTriangulations(A, Fine=>false, ConnectedToRegular=>false, RegularOnly => false) -- 115 sec, 564 of these
elapsedTime n5 = numTriangulations(A, Fine=>true, ConnectedToRegular=>true, RegularOnly=>false) -- .09 sec, 448 of these
--this is n2 above. elapsedTime n6 = numTriangulations(A, Fine=>true, ConnectedToRegular=>false, RegularOnly=>false) -- 115.5 sec, 448 of these
elapsedTime n7 = numTriangulations(A, Fine=>false, ConnectedToRegular=>true, RegularOnly=>false) -- .11 sec, 564 of these
-- this is n4 above. elapsedTime n8 = numTriangulations(A, Fine=>false, ConnectedToRegular=>false, RegularOnly=>false) -- 116 sec, 564 of these
elapsedTime set1 = allTriangulations(A, Fine=>true, ConnectedToRegular=>true); -- 6.9 sec, 408 CORRECT
elapsedTime set2 = allTriangulations(A, Fine=>true, ConnectedToRegular=>false, RegularOnly => false); -- 118 sec, 448
elapsedTime set3 = allTriangulations(A, Fine=>false, ConnectedToRegular=>true); -- 8.1 sec, 520 CORRECT
elapsedTime set4 = allTriangulations(A, Fine=>false, ConnectedToRegular=>false, RegularOnly => false); -- 116 sec. 564 of these.
elapsedTime set5 = allTriangulations(A, Fine=>true, ConnectedToRegular=>true, RegularOnly=>false); -- .15 sec, 448 of these
--elapsedTime set6 = allTriangulations(A, Fine=>true, ConnectedToRegular=>false, RegularOnly=>false); -- 116 sec, 448 of these
elapsedTime set7 = allTriangulations(A, Fine=>false, ConnectedToRegular=>true, RegularOnly=>false); -- .22 sec, 564 of these
--elapsedTime set8 = allTriangulations(A, Fine=>false, ConnectedToRegular=>false, RegularOnly=>false); -- 117 sec, 564 of these
assert((n1,n2,n3,n4,n5,n7) == (#set1, #set2, #set3, #set4, #set5, #set7))
fineTris = select(set4, x -> # unique flatten x == numColumns A);
regularFineTris = select(fineTris, x -> isRegularTriangulation(A, x));
regularTris = select(set4, x -> isRegularTriangulation(A, x));
assert(#regularFineTris == 408)
assert(#fineTris == 448)
assert(#regularTris == 520)
assert(set set5 === set set2) -- in general, this doesn't need to hold, but it is rare for this to be the case
assert(set set7 === set set4) -- same: rare for this to not hold
elapsedTime assert(set select(set7, x -> isRegularTriangulation(A, x)) === set set3)
elapsedTime assert(set select(set5, x -> isRegularTriangulation(A, x)) === set set1)
-- the rest of this test needs to be checked, and cleaned up
set5_0
elapsedTime for tri in set5 list naiveIsTriangulation(A, tri)
numFlips(A, set5_0)
topcomFlips(A, set5_0)
-- now let's see about the naive way of getting regular star triangulations
-- i.e. we add in the origin
pts1 = {{-1,0,0,-1},{-1,0,1,-1},{-1,0,1,0},{-1,1,0,-1},{-1,1,0,0},{-1,1,1,2},{1,-1,0,-1},{1,0,-1,1},{1,-1,-1,-1},{0,0,0,-1},{0,0,0,0}}
A1 = transpose matrix pts1
--elapsedTime tris1 = allTriangulations(A1, Fine=>true, ConnectedToRegular=>true, RegularOnly=>false); --
elapsedTime tris1 = allTriangulations(A1, Fine=>false, ConnectedToRegular=>true, RegularOnly=>false); --
fineTris1 = select(tris1, x -> # unique flatten x == numColumns A1);
regTris1 = select(tris1, x -> isRegularTriangulation(A1, x));
fineRegTris1 = select(regTris1, x -> # unique flatten x == numColumns A1);
stars1 = select(tris1, x -> all(x, x1 -> member(10, x1))); -- 100 here
starsFine1 = select(stars1, x -> # unique flatten x == numColumns A1);
RST = select(stars1, x -> isRegularTriangulation(A1,x)); -- 80 here...
FSRT = select(starsFine1, x -> isRegularTriangulation(A1,x)); -- 48 here...!
unique for tri in set5 list fineStarTriangulation(A, tri);
--select(oo, tri -> isRegularTriangulation(A1, tri)) -- this isn't defined here...
-- let's test this one for being a triangulation:
oA = orientedCircuits A
tri = set5_3
tally flatten for t in tri list subsets(t,4)
for z in oA list (
# select(tri, t -> isSubset(z_0, t)),
# select(tri, t -> isSubset(z_1, t))
)
-- todo:
-- 1. routine to check that a triangulation is a triangulation
-- 2. routine to turn a regular, fine triangulation, into a star (fine, regular) triangulation. How general is this? DONE, I think.
-- 3. perform bistellar flips to get new triangulations.
///
TEST ///
restart
needsPackage "ReflexivePolytopesDB"
needsPackage "Topcom"
needsPackage "Polyhedra"
debug needsPackage "Triangulations"
str = "4 18 M:53 18 N:11 10 H:6,45 [-78]
1 0 0 -2 0 2 1 3 -2 2 -2 2 3 -1 0 -2 0 -1
0 1 0 2 0 0 1 -2 1 -2 0 0 -1 0 -2 0 -2 -1
0 0 1 1 0 -1 -1 -2 2 -2 0 -2 -2 2 -1 2 1 2
0 0 0 0 1 -1 -1 0 -1 1 1 -1 -1 -1 2 -1 0 -1"
A = matrix first kreuzerSkarke str
P = convexHull A
P2 = polar P
A1 = vertices P2
LP = matrix{select(latticePoints P2, x -> x != 0)}
numTriangulations(LP, Fine => true) == 731
elapsedTime tris = allTriangulations(LP, Fine=>true);
#tris == 731
numTriangulations(LP)
allTriangulations(LP);
T = regularFineTriangulation LP
elapsedTime trisT = generateTriangulations T; -- much slower, and includes 10 non-regular triangulations.
trisT/(T -> isRegularTriangulation T)//tally
(24) * volume P2 == 29
///
TEST ///
-- this is an example used in
-- https://people.inf.ethz.ch/fukudak/lect/mssemi/reports/03_rep_ClemensPohle.pdf
-- (accessed 2 June 2022)
-*
restart
needsPackage "Triangulations"
*-
debug needsPackage "Topcom"
topcompath = "/opt/homebrew/bin/"
-- this tests construction of the chirotope.
A = transpose matrix"0,0;1,1;3,1;5,0;1,5"
OM = chirotope A
-- "foo-input.in" << topcomPoints A << endl << close;
-- directString = get("!"|(topcompath|"points2chiro c => bistellarFlip(tri, c))
24 * gkzVector tri
for t in nbors list if last t === null then continue else {t, gkzVector last t}
gkzVector tri
///
TEST /// -- test of functions here on the square and the cube
-*
restart
needsPackage "Triangulations"
*-
square = transpose matrix{{1,1},{-1,1},{-1,-1},{1,-1}}
topcomRegularFineTriangulation square
assert(# allTriangulations square == 2)
-- Now consider all of the lattice points of the square
sq9 = matrix {{-1, -1, 1, 1, -1, 0, 0, 1, 0}, {-1, 1, -1, 1, 0, -1, 1, 0, 0}}
-- test function from Topcom.
t1 = topcomRegularFineTriangulation sq9
regularTriangulationWeights(sq9, t1)
fineStarTriangulation(sq9, t1) -- has central element removed. Don't do that?
delaunaySubdivision sq9 -- not a triangulation (4 squares).
orientedCircuits sq9 -- many of these are not useful when considering only fine triangulations.
-- Polyhedra command, bug fix currently in StringTorics.
t2 = regularSubdivision(sq9, matrix{{4,4,4,4,1,1,1,1,0}})
assert(
t2
==
{{0, 4, 5}, {1, 4, 6}, {2, 5, 7}, {3, 6, 7}, {4, 5, 8}, {4, 6, 8}, {5, 7, 8}, {6, 7, 8}}
)
t3 = regularSubdivision(sq9, matrix{{4,4,4,4,1,1,1,1,-4}})
t4 = regularSubdivision(sq9, matrix{{1,1,1,1,1,1,1,1,-4}})
affineCircuits(sq9, t3)
elapsedTime affineCircuits(sq9, t1) -- TODO: change the name of this function? Why? affineOrientedCircuits?
elapsedTime tris = allTriangulations sq9; -- computes all triangulations, Fine => false, regular.
assert(387 == #tris) -- check this number.
elapsedTime tris1 = allTriangulations(sq9, Fine => true);
assert(64 == #tris1)
elapsedTime tris2 = generateTriangulations triangulation(sq9, t1);
elapsedTime tris2a = generateTriangulations(triangulation(sq9, t1), RegularOnly => true); -- slow...
#tris2 == #tris2a
-- via topcom: (too long for general test)
----elapsedTime assert(387 == # select(tris, t -> topcomIsRegularTriangulation(sq9, max t))) -- slow
-- via Polyhedra/ (too long for a general test)
----elapsedTime assert(387 == # select(tris, t -> null =!= regularTriangulationWeights t)) -- slow, same as Polyhedra, I think (same code is being used)
-- 64 fine triangulations
assert(64 == # select(tris, isFine))
assert(16 == # select(tris, isStar))
assert(1 == # select(tris, t -> isStar t and isFine t))
-- generate only fine
elapsedTime finetris = allTriangulations(sq9, Fine => true);
-- elapsedTime finetris1 = generateTriangulations(sq9, t1, Fine => true); -- TODO
elapsedTime finetris1 = generateTriangulations triangulation(sq9, t1);
assert(finetris/max//sort == finetris1/max//sort)
finetris/volumeVector
-- TODO: audit/allow "A" matrices to be homogeneous or not? THIS MIGHT BE DONE...
-- one way: default is not, and Homogenize => true will add the extra row.
-- use this for most routines taking point configurations or vector configurations.
///
TEST ///
-- triangulation code, test 1.
-- some triangulation code is in Topcom, Polyhedra
-- (comes from KS database, h11=3, #232.
-*
restart
needsPackage "Triangulations"
*-
A0 = matrix {{1, 1, 1, 1, -11}, {0, 2, 2, 2, -10}, {0, 0, 4, 4, -8}, {0, 0, 0, 12, -12}}
P2 = polar convexHull A0
-- the floowing are the lattice points of P2 (from 'LP = latticePointList P2', but this uses StringTorics...)
LP = {{-1, 0, 0, 0}, {-1, 0, 0, 1}, {-1, 0, 3, -1}, {-1, 2, -1, 0}, {1, -1, 0, 0}, {-1, 0, 1, 0}, {-1, 1, 0, 0}, {0, 0, 0, 0}}
A = transpose matrix LP
-- We want triangulations of the point configuration given by the columns of A.
T = topcomRegularFineTriangulation A
assert isFine(A, T)
assert isStar(A, T)
assert isRegularTriangulation(A, T)
assert naiveIsTriangulation(A, T)
assert topcomIsTriangulation(A, T)
elapsedTime tris = allTriangulations A;
FRST = first select(tris, t -> isFine t and isStar t and isRegularTriangulation t)
volumeVector FRST
elapsedTime tris = generateTriangulations(A, T);
wts = regularTriangulationWeights(A, T)
assert(regularSubdivision(A, matrix{wts}) == T//sort)
assert isRegularTriangulation(A, T) -- from Topcom.
assert topcomIsTriangulation(A, T)
-- also check Topcom's code...
neighbors triangulation(A, T) -- not unpacked... TODO: need to understand and document what this gives
--A1 = A || splice matrix{{numColumns A: 1}}
--A1 = augment A
A1 = A || matrix{{numcols A: 1}}
--volumeVector(A, T) -- need to homogenize A? at bottom or at top?
volumeVector(A1, T) -- need to homogenize A? at bottom or at top?
neighbors triangulation(A1, T)
--regularFineStarTriangulation A1 -- this one fails
regularFineStarTriangulation A -- this one works
-- check that T is a triangulation?
circs = affineCircuits(A, T) -- NEEDS CLEANING UP (? why?)
bistellarFlip(T, {{5}, {0, 1, 2}})
for c in circs list bistellarFlip(T, c)
///
-- The following is too long for a test.
///
-- how many triangulations at h11=8? (roughly?)
needsPackage "ReflexivePolytopesDB"
topes = kreuzerSkarke 8;
A = matrix topes_500
P = convexHull A
P2 = polar P
LP = matrix{latticePoints P2}
T = regularFineTriangulation LP
T = {T}
T1 = T/neighbors/(t -> t/last)//flatten
T1/gkzVector
#T1
T2 = T1/(t -> neighbors t)//flatten/last//unique;
T2/gkzVector//matrix
T1 = T2;
methods allTriangulations
options allTriangulations
debugLevel = 1
elapsedTime trisLP = allTriangulations(LP, Fine => true, RegularOnly => true);
frsts = for t in trisLP list if all(t, t1 -> member(numcols LP - 1, t1)) then t else continue;
#frsts == 108
isStar(t, numcols LP - 1) then t else continue;
T = regularFineTriangulation LP
debug Triangulations
debugLevel = 0
elapsedTime generateTriangulations(T, Limit => 20000);
///
end----------------------------------------------------
restart
uninstallPackage "Triangulations"
restart
needsPackage "Triangulations"
restart
installPackage "Triangulations"
restart
check "Triangulations"
TEST ///
debug needsPackage "Topcom"
-- points2chiro
toppath = "/opt/homebrew/bin/"
A = transpose matrix {{-1,-1,1},{-1,1,1},{1,-1,1},{1,1,1},{0,0,1}}
tri = {{0, 2, 4}, {2, 3, 4}, {0, 1, 4}, {1, 3, 4}}
run (toppath|"/points2chiro"|" -h")
"topcomfoo.in" << topcomPoints(A, Homogenize=>false) << endl << close;
chiro = get ("!"|toppath|"/points2chiro"|" 0 then "+" else if d == 0 then "0" else "-"
)) | "\n"
chiro == chiro2
-- notes: a chirotope for topcom:
-- 5,3: (number of vertices, dim)
-- a string of "-","+","0", maybe cut over a number of lines.
-- should we make a type out of this (so we can read and write it to a file)
-- chiro2circuits
"topcomfoo.in" << chiro << endl << [] << endl << close;
circs = get ("!"|toppath|"/chiro2circuits"|" 10);
tope = polytopes_5
A = matrix tope
P = convexHull A
P2 = polar P
A = matrix{latticePoints P2}
LP = drop(latticePointList P2, -1);
A = transpose matrix LP;
debugLevel = 6
elapsedTime tri = regularFineTriangulation A;
-- XXX
augment A
"topcomfoo.in" << topcomPoints(augment A, Homogenize=>false) << endl << close;
chiro = get ("!"|toppath|"/points2chiro"|" true) << endl << "[]" << endl << tri << endl << close;
run (topcompath|"checkregularity"|" --heights topcomfoo.out")
-- points2chiro
"topcomfoo.in" << topcomPoints(A, Homogenize=>false) << endl << "[]" << endl << badtri << endl << close;
print (toppath|"/points2alltriangs"|" --checktriang -v true) << endl << "[]" << endl << tri << endl << close;
run (topcompath|"checkregularity"|" --heights topcomfoo.out")
assert not isRegularTriangulation(A,tri)
///
///
-- now let's see about the naive way of getting regular star triangulations
-- i.e. we add in the origin
pts1 = {{-1,0,0,-1},{-1,0,1,-1},{-1,0,1,0},{-1,1,0,-1},{-1,1,0,0},{-1,1,1,2},{1,-1,0,-1},{1,0,-1,1},{1,-1,-1,-1},{0,0,0,-1},{0,0,0,0}}
A1 = transpose matrix pts1
elapsedTime tris1 = allTriangulations(A1, Fine=>false, ConnectedToRegular=>true, RegularOnly=>false); --
fineTris1 = select(tris1, x -> # unique flatten x == numColumns A1);
regTris1 = select(tris1, x -> isRegularTriangulation(A1, x));
fineRegTris1 = select(regTris1, x -> # unique flatten x == numColumns A1);
stars1 = select(tris1, x -> all(x, x1 -> member(10, x1))); -- 100 here
starsFine1 = select(stars1, x -> # unique flatten x == numColumns A1);
RST = select(stars1, x -> isRegularTriangulation(A1,x)); -- 80 here...
FSRT = select(starsFine1, x -> isRegularTriangulation(A1,x)); -- 48 here...!
#tris1 == 2254
#RST == 80
#FRST = 48
unique for tri in set5 list (
tri1 := fineStarTriangulation(A, tri);
newtri := for t in tri1 list append(t, 10);
newtri
);
select(oo, tri -> isRegularTriangulation(A1, tri))
-- let's test this one for being a triangulation:
oA = orientedCircuits A
tri = set5_3
tally flatten for t in tri list subsets(t,4)
for z in oA list (
# select(tri, t -> isSubset(z_0, t)),
# select(tri, t -> isSubset(z_1, t))
)
///
TEST ///
-- Bad triangulations of the square
V = transpose matrix {{0,0},{1,0},{0,1},{1,1}}
T1 = {{0,1,2}}
T2 = {{0,1,2},{0,1,3}}
T3 = {{0,1,2,3}}
assert(not topcomIsTriangulation(V, T1))
assert(not topcomIsTriangulation(V, T2))
assert(not topcomIsTriangulation(V, T3))
debug needsPackage "Triangulations" -- TODO: isTriangulation should be exported, or called naiveIsTriangulation.
assert(not isTriangulation(V, T1))
assert(not isTriangulation(V, T2))
assert(not isTriangulation(V, T3)) -- gives error, should give false!
///
end--
restart
uninstallPackage "Topcom"
restart
needsPackage "Topcom"
installPackage "Topcom"
restart
check "Topcom"
viewHelp
///
-- generate examples to use for this package
-- from reflexive polytopes of dim 4
restart
needsPackage "StringTorics"
str = getKreuzerSkarke(10, Limit=>5)
str = getKreuzerSkarke(20, Limit=>5)
str = getKreuzerSkarke(30, Limit=>5)
polytopes = parseKS str
tope = polytopes_4_1
A = matrixFromString tope
P = convexHull A
P2 = polar P
LP = drop(latticePointList P2, -1)
A1 = transpose matrix LP
A2 = transpose matrix latticePointList P2
tri = regularFineTriangulation A1
tri2 = regularFineTriangulation A2
#tri
#tri2
elapsedTime chiro1 = chirotope A1;
elapsedTime chiro2 = chirotope A2;
elapsedTime # orientedCircuits chiro1
elapsedTime # orientedCircuits chiro2
elapsedTime # orientedCocircuits chiro1
elapsedTime # orientedCocircuits chiro2
(select(orientedCocircuits A2, f -> #f#0 == 0 or #f#1 == 0))/first
netList annotatedFaces P2
tri2
-- fine:
assert(sort unique flatten tri2 == toList (0..14))
walls = tri2/(x -> subsets(x, #x-1))//flatten
nfacets = tally walls
facs = (select((annotatedFaces P2), x -> x_0 == 3))/(x -> x#2)
walls = partition(k -> nfacets#k, keys nfacets)
for w in walls#1 list (
# select(facs, f -> isSubset(w, f))
)
for w in walls#2 list (
# select(facs, f -> isSubset(w, f))
)
-- check overlaps of elements of tri2:
C = orientedCircuits A2;
elapsedTime for c in C list (
val1 := select(tri2, x -> isSubset(c_0, x));
val2 := select(tri2, x -> isSubset(c_1, x));
if #val1 > 0 and #val2 > 0 then print (c, val1, val2);
(c, #val1, #val2));
tri_0
///
doc ///
Key
Headline
Usage
Inputs
Outputs
Consequences
Description
Text
Example
Caveat
SeeAlso
///