Why should we care?
Triangulating a polyhedron is actually interesting to physics (Or if you are more conservative, at least to string theory). The reason for this is toric geometry, a functor from certain polyhedra to algebraic varieties. And in nice cases those contain Calabi-Yau hypersurfaces. And if you believe in string theory, then those Calabi-Yau's might be part of our spacetime...
Puntos
Jesús A. De Loera has written a program to triangulate point sets, it is reasonably easy to use and does the job. However there is one problem: It is written entirely in Maple and sometimes it is just too slow.
The solution
Since I wanted to triangulate too many points for the standard Puntos, I had to do something. So I rewrote the main loop in C++, this yields about an order of magnitude faster execution/less memory usage.
Here are the necessary files:
-
A example session, complete
with the necessary maple scripts and a statically linked version
of the C++ code, should run on any i386 ELF (i.e. linux)
platform.
Since gcc's C++ libraries change quite often I decided to to link it statically. This is the reason for the huge (3.5M) binary size.
- The C++ source.
Just get the example session, unpack it and run "maple samplesession.txt". This will run a small computation, it should be easy to change it to your needs.
Here is the content of "samplesession.txt":
# run with "maple samplesession.txt"
# read my (modified) puntos code
read(`puntos_maple`):
# the first point must be the origin!
U:=transpose(matrix( [
[1,0,0],
[1,-1,0],
[1,0,-1],
[1,0,1],
[1,1,2],
[1,1,3]
] ));
# this is a modified part of puntos.
# it calculates the flips and a seed triangulation, and then stops
# (without doing the tedious work, which is done a lot faster
# by my C++ version of it)
# Note that using symmetries to simplify the computation is not supported
grafo_step1:=
proc(ma)
local seed,raiz,i,cell:
global cocircuits,circuitos,gale,symmetry_moves;
symmetry_moves:=[]:
gale:=matrix(convert(kernel(ma),list)):
circo(gale):
circuitos:=cocircuits:
raiz:=get_a_triangulation(circuitos):
seed:={}:
for cell in raiz do
seed:=seed union {{$1..coldim(ma)} minus cell}:
od:
flipper( ):
chirotope(ma):
cheap_representative(seed,'q'):
end:
seed:=grafo_step1(U);
n:=coldim(U);
d:=rowdim(U);
# generate the input file for my program
# the variables seed,n,d,flips must be saved in that order!
save(seed,n,d,flips,`puntos.input`);
# run the external program (the example is for unix)
# You could run it by hand and redirect the output to a file to
# let it run in the background
# if you changed the output to "only star triangulations"
# (see main.cc) origin must be first in U!
system( `puntos_cc` );
# read the generated result
# the simplices are enumerated as by the orden() function in puntos
read(`puntos.output`);
# now use the function reverse_orden() in my modified version of puntos
# to find the actual vertices for each simplex
out_triangs:=[]:
for triang in triangulations do
t:=map( i->reverse_orden(i,n,d), triang);
out_triangs := [op(out_triangs), t];
od:
nops(out_triangs);
save(U,out_triangs,`./triangulations`);
# finally quit
quit;
and it will generate the following output ("samplesession.output"):
|\^/| Maple V Release 4 (Humboldt-Universitaet zu Berlin)
._|\| |/|_. Copyright (c) 1981-1996 by Waterloo Maple Inc. All rights
\ MAPLE / reserved. Maple and Maple V are registered trademarks of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
# run with "maple samplesession.txt"
>
# read my (modified) puntos code
> read(`puntos_maple`):
Warning, new definition for Chi
Warning, new definition for maximize
Warning, new definition for minimize
Warning, new definition for basis
Warning, new definition for charpoly
Warning, new definition for fibonacci
Warning, new definition for norm
Warning, new definition for pivot
Warning, new definition for rank
Warning, new definition for trace
>
# the first point must be the origin!
> U:=transpose(matrix( [
> [1,0,0],
> [1,-1,0],
> [1,0,-1],
> [1,0,1],
> [1,1,2],
> [1,1,3]
> ] ));
[1 1 1 1 1 1]
[ ]
U := [0 -1 0 0 1 1]
[ ]
[0 0 -1 1 2 3]
>
# this is a modified part of puntos.
# it calculates the flips and a seed triangulation, and then stops
# (without doing the tedious work, which is done a lot faster
# by my C++ version of it)
# Note that using symmetries to simplify the computation is not supported
> grafo_step1:=
> proc(ma)
> local seed,raiz,i,cell:
> global cocircuits,circuitos,gale,symmetry_moves;
>
> symmetry_moves:=[]:
> gale:=matrix(convert(kernel(ma),list)):
> circo(gale):
> circuitos:=cocircuits:
> raiz:=get_a_triangulation(circuitos):
> seed:={}:
> for cell in raiz do
> seed:=seed union {{$1..coldim(ma)} minus cell}:
> od:
> flipper( ):
> chirotope(ma):
> cheap_representative(seed,'q'):
> end:
>
> seed:=grafo_step1(U);
bytes used=1000080, alloc=851812, time=0.06
I have computed, 1
I have computed, 2
I have computed, 3
I have computed, 4
I have computed, 5
I have computed, 6
I have computed, 7
I have computed, 7
I have computed, 8
I have computed, 9
I have computed, 9
I have computed, 9
I have computed, 9
I have computed, 10
I have computed, 11
this polytope has, 0, facets
bytes used=2000324, alloc=1310480, time=0.12
seed := [11, 15, 17, 20]
>
> n:=coldim(U);
n := 6
> d:=rowdim(U);
d := 3
>
# generate the input file for my program
# the variables seed,n,d,flips must be saved in that order!
> save(seed,n,d,flips,`puntos.input`);
>
# run the external program (the example is for unix)
# You could run it by hand and redirect the output to a file to
# let it run in the background
# if you changed the output to "only star triangulations"
# (see main.cc) origin must be first in U!
> system( `puntos_cc` );
Seed triangulation: [11,15,17,20]
Vertices of seed triangulation: [2_3_4 , 2_4_6 , 3_4_5 , 4_5_6]
All flips:
< 2_5_6 2_3_5 | 2_3_6 3_5_6 >
< 2_3_6 3_5_6 | 2_5_6 2_3_5 >
< 2_4_6 2_3_4 3_4_6 | 2_3_6 >
< 2_3_6 | 2_4_6 2_3_4 3_4_6 >
< 2_4 4_5 | 2_5 >
< 2_5 | 2_4 4_5 >
< 3_5_6 3_4_6 | 3_4_5 4_5_6 >
< 3_4_5 4_5_6 | 3_5_6 3_4_6 >
< 2_5_6 1_2_5 | 1_2_6 1_5_6 >
< 1_2_6 1_5_6 | 2_5_6 1_2_5 >
< 1_2_3 1_3_5 1_2_5 | 2_3_5 >
< 2_3_5 | 1_2_3 1_3_5 1_2_5 >
< 1_3_6 1_2_6 1_2_3 | 2_3_6 >
< 2_3_6 | 1_3_6 1_2_6 1_2_3 >
< 1_4 1_3 | 3_4 >
< 3_4 | 1_4 1_3 >
< 1_3_5 1_5_6 | 1_3_6 3_5_6 >
< 1_3_6 3_5_6 | 1_3_5 1_5_6 >
< 1_5_6 1_4_6 | 4_5_6 1_4_5 >
< 4_5_6 1_4_5 | 1_5_6 1_4_6 >
< 2_4_6 1_4_6 1_2_4 | 1_2_6 >
< 1_2_6 | 2_4_6 1_4_6 1_2_4 >
Finished: found 10 triangulations in 0s
0
>
# read the generated result
# the simplices are enumerated as by the orden() function in puntos
> read(`puntos.output`);
triangulations := [[1, 4, 6, 10]]
>
# now use the function reverse_orden() in my modified version of puntos
# to find the actual vertices for each simplex
> out_triangs:=[]:
> for triang in triangulations do
> t:=map( i->reverse_orden(i,n,d), triang);
> out_triangs := [op(out_triangs), t];
> od:
> nops(out_triangs);
1
>
> save(U,out_triangs,`./triangulations`);
>
# finally quit
> quit;
bytes used=2050124, alloc=1310480, time=0.12