interface(quiet=true): #PUNTOS version 3, sept 15 1994 # Copyright 1995, Jesus de Loera, deloera@geom.umn.edu # This software can be distributed freely. Please acknowledge the use # this software if it is of use on your research. # A small manual can be found at the file manual.ps # latest revision 24 november 1996. # Now the command grafo determines which triangulations are regular # only if asked to do so. A question will be presented to the user # before embarking on the solution of LP's. ############################################################################### # PUNTOS is a collection of subroutines to study TRIANGULATIONS and SECONDARY POLYTOPES # # for POINT CONFIGURATIONS. The main operation is the subroutine ``GRAFO'' that, # # starting with an initial regular triangulation for a point configuration A, computes the # # 1-skeleton of the Secondary polytope. Very often one is not only presented with an arbi- # # point configuration but has additional information. In particular we can use the group # # of automorphisms of the configuration to speed up the generation of vertices of the # # Secondary polytope. This observation is implemented in the subroutine ``GRAFO'' # # (Note: Only the vertices are generated in this case and no record of the edges is kept). # # The way to input the group is explained later. The list TRIANGULATIONS contains all the # # simplicial complexes defining the triangulations. One can obtain the canonical Gelfand- # # Kapranov-Zelevinskii embedding using the subroutine ``GKZ_EMBEDDING''. The GKZ coordinates# # of a single triangulation can be computed using the procedure PUNTITOS. # #In addition to the main core # # of subroutines PUNTOS also contains procedures that can deal with some basic operations # # in Oriented Matroids. We can compute chirotope, circuits or cocircuits for a point con- # # figuration. # ############################################################################################# with(linalg): with(combinat): with(networks): with(simplex): #--------------------GLOBAL VARIABLES-------------------------------------- # x[i1,i2,..,ik] <------This is the chirotope of a matrix (CALL procedure chirotope) # cocircuits <--------(CALL procedure circo) # circuitos <--------(CALL procedure circo with the orthogonal matrix) # flips <--------(CALL procedure flipper) # triangulations <-------(CALL procedure grafosecondary or ACTIONsecondary. When used with ACTIONsecondary saves # only the orbit representatives.) # looseends <--------(used in procedure grafosecondary, to keep track of the breadth-first search) # skeleton <---------(CALL procedure grafosecondary) # gale <--------------(The Gale transform of a point configuration); # nonregulars <---------(non-regular triangulations appearing during the process) # retriangus <----------(The regular triangulations of the configuration). # gkzvectors <----------(This saves, modulo the symmetries of the configuration, the distinct # Gelfand-Kapranov-Zelevinskii vectors of a configuration). # totaltriang <-------(The total number of triangulations, NOT considering the symmetries acting on the configuration # semilla <--------(This is a matrix that encodes the seed triangulation. It comes from the file matrizaux # symmetry_moves <----( This is a list of lists that represents the group of symmetries of the point configuration # This is read from an external file whose name is a parameter of ACTIONsecondary) # permutacion<-----(This contains a list representing a permutation, it is filled as it is applied in the subroutine # muevepunto. permutacion is an element of symmetry_moves). # parimpar <--------Keeps score of the parity of a permutacion represented as a list, used in the subroutine paridad. #-------------------------------------------------------------------- # This is a small procedure to determine the parity of a permutation. # The permutation parameter is passed as a list. #--------------------------------------------------------------------- mergep:=proc(L1,L2) local i,j,total,cuenta: total:=0: for i from 1 to nops(L1) do cuenta:=0: for j from 1 to nops(L2) do if op(i,L1)>op(j,L2) then cuenta:=cuenta+1: fi: od: total:=total+cuenta: od: total; end: paridad:=proc(L,i,j) local m,tempo: global parimpar: if i=j then RETURN([op(i,L)]); else m:=floor((i+j-1)/2); parimpar:=parimpar+mergep(sort(paridad(L,i,m)),sort(paridad(L,m+1,j))); tempo:=sort([op(paridad(L,i,m)),op(paridad(L,m+1,j))]): RETURN(tempo); fi: end: #--------------------------------------------------------------------------- # This subroutine computes the chirotope of a given n by k # matrix. This is obtained from the k by k minors of the matrix. # We rely on the Laplace expansion formula and the distributivity of # tensor multiplication to compute the n choose k minors "at the same time" #--------------------------------------------------------------------------- tensar:=proc(val1,val2,L1,L2) local aux,s1,s2,h,temp: global x: s1:=convert(L1,set): s2:=convert(L2,set): temp:=0: h:=nops(s1 intersect s2): if h=0 then aux:=sort([op(L1),op(L2)]): temp:=val1*val2*((-1)^mergep(L1,L2)): #print(x[op(aux)]); x[op(aux)]:=x[op(aux)]+temp: # print(aux,temp,x[op(aux)]); fi: RETURN(NULL): end: wedge:=proc(v,w,t,s,n) local i,T,S,W,res,svec,tvec: global x: W:=choose([$1..n],t+s): for i from 1 to nops(W) do x[op(op(i,W))]:=0: od: T:=choose([$1..n],t): S:=choose([$1..n],s): for tvec from 1 to nops(T) do for svec from 1 to nops(S) do tensar(v[tvec],w[svec],op(tvec,T),op(svec,S)): od: od: res:=[]: for i from 1 to nops(W) do res:=[op(res),x[op(op(i,W))]]: od: res; end: chirotope:=proc(a) local n,k,i,temp: n:=coldim(a); k:=rowdim(a); temp:=wedge(convert(submatrix(a,[1],[$1..n]),vector),convert(submatrix(a,[2],[$1..n]),vector),1,1,n); for i from 3 to k do temp:=wedge(temp,convert(submatrix(a,[i],[$1..n]),vector),i-1,1,n): od: temp; end: #-------------------------------------------------------------------- # This subroutine computes the cocircuits or circuits of a matrix #-------------------------------------------------------------------- circo:=proc(a) local i,j,zeros,neg,T,pos,cuentalas,n,k,aux,tempo,cocircuit,auxilio: global x,cocircuits,parimpar,orbitasimplejo,orbitacirco,volley: volley:={}: chirotope(a): k:=rowdim(a): n:=coldim(a): T:=choose(n,k-1); cocircuits:={}: while nops(T)>0 do zeros:={}: neg:={}: pos:={}: for j from 1 to n do if not(member(j,op(1,T))) then parimpar:=0: tempo:=paridad([op(op(1,T)),j],1,k): aux[j]:=((-1)^parimpar)*x[op(tempo)]: else aux[j]:=0: fi: if aux[j]=0 then zeros:=zeros union {j}: elif sign(aux[j])=1 then pos:=pos union {j}: else neg:=neg union {j}: fi: od: cocircuit:=[zeros,{pos,neg}]: auxilio:=orbitacirco(cocircuit): print(`cuantos`,nops(auxilio)); if nops(auxilio)<2 then volley:=volley union auxilio: fi: cocircuits:=cocircuits union auxilio: print(`I have computed`,nops(cocircuits)); T:=T minus orbitasimplejo(op(1,T)); print(nops(T)); print(`to go`); od: cocircuits:=cocircuits minus {[{$ 1..n },{{}}]}: cocircuits:=map(x->[x[2][1],x[1],x[2][2]],cocircuits): cocircuits:= convert(cocircuits,list): cuentalas:=0: for j from 1 to nops(cocircuits) do if (op(1,op(j,cocircuits))={ }) or (op(3,op(j,cocircuits))={ }) then cuentalas:=cuentalas+1: print(`facet number`, cuentalas, op(j,cocircuits)): fi: od: print(`this polytope has`, cuentalas, `facets`); end: #----------------------------------------------------------------------------------------------------- # In the following procedure we compute the orbit of a given triangulation. Previously # we have read the SYMMETRY_MOVES that encode the symmetry group of the point configuration. # The procedure consists simply in applying one by one the permutations (in a series of replacements). # The parameter Ti is the triangulation in question. #----------------------------------------------------------------------------------------------------- orbitasimplejo:=proc(Ti) local j,l,newTi,aux,numpoints; global orbit,muevesimplejo,permutacion; orbit:={Ti}: for j from 1 to nops(symmetry_moves) do permutacion:=op(j,symmetry_moves): newTi:=muevesimplejo(Ti): orbit:=orbit union {newTi}: od: RETURN(orbit); end: orbitacirco:=proc(Ti) local j,l,newTi,aux,numpoints; global permutacion,symmetry_moves,muevecirco,orbit; orbit:={Ti}: for j from 1 to nops(symmetry_moves) do permutacion:=op(j,symmetry_moves): newTi:=muevecirco(Ti): orbit:=orbit union {newTi}: od: orbit; end: muevecirco:=proc(circo) local u: u:=[muevesimplejo(circo[1]),muevetriangulacion(circo[2])]; RETURN(u); end: #---------------------------------------------- # Minimalizes a set of nonfaces #---------------------------------------------- sizz := proc(a,b) evalb( nops(a) > nops(b)) end: minimalize := proc(NF) local erg,dada,i,j: erg := NF: dada := sort(convert(NF,list),sizz): for i from 1 to nops(dada)-1 do for j from i+1 to nops(dada) do if ((dada[j] minus dada[i]) = {}) then erg := erg minus {dada[i]}: fi: od:od: erg; end: #---------------------------------------------------------------------------- # This procedure computes a lexicographic triangulation from the list of circuits # of the point configuration #----------------------------------------------------------------------------- get_a_triangulation := proc(Circuits) local erg,mm,circ: global minimalize: erg := {}: for circ in Circuits do mm := min(convert(circ[1] union circ[3],list)[]): if (member(mm,circ[1])) then erg := erg union {circ[1]}: else erg := erg union {circ[3]}: fi: od: mmt(minimalize(erg)); end: #---------------------------------------------------------------- # Computes a minimal set of representatives of a family of sets. #----------------------------------------------------------------- red := proc( fam ) local badguys,i,j: badguys := {}: for i from 1 to nops(fam) do for j from 1 to nops(fam) do if ((i<>j) and (fam[i] minus fam[j] = {})) then badguys := badguys union {fam[j]}: fi: od:od: fam minus badguys end: mmt := proc( fam ) local erg,i,j,head,tail: erg := {}: if (fam = {}) then erg := { {} }: fi: if (fam <> {}) then head := fam[1]: tail := mmt( fam minus {fam[1]}): for i from 1 to nops(head) do for j from 1 to nops(tail) do erg := erg union {{head[i]} union tail[j]}: od: od: erg := red(erg): fi: erg end: #------------------------------------------------------------------------------------- # In the following subroutine we compute the possible FLIPS. A flip is a # pair of simplicial complexes inside a BIG simplicial complex (say a triangulation of a polytope) # that can be interchanged using a circuit (The linear structure allows this exchange when the two simplicial # complexes have the same convex hull and share a linear dependence). #------------------------------------------------------------------------------------- flipper:=proc( ) local j,i,neg,pos,delta,deltaminus,deltaplus,aux,simplexo: global flips,circuitos: flips:=[]: for j from 1 to nops(circuitos) do pos:=op(3,op(j,circuitos)): neg:=op(1,op(j,circuitos)): delta:=pos union neg: pos := convert(pos,list): neg := convert(neg,list): deltaplus:={}: deltaminus:={}: for i from 1 to nops(pos) do simplexo:=delta minus {op(i,pos)}: deltaplus:=deltaplus union {simplexo}: od: for i from 1 to nops(neg) do simplexo:=delta minus {op(i,neg)}: deltaminus := deltaminus union {simplexo}: od: aux:=[deltaplus,deltaminus]; flips:=[op(flips),aux]: od: #for i from 1 to nops(circuitos) do # print(`circuit#`,i,`:``pos:=`,op(3,op(i,circuitos)),`neg:=`,op(1,op(i,circuitos))); # print(`deltaplus:=`,op(1,op(i,flips)),`deltaminus:=`,op(2,op(i,flips))); #od: end: #--------------------------------------------------------------------------------- # This procedure computes The facet-inequalities defining the relative interior of # a simplicial cone (provided we are given the extreme rays). the parameter matrix # a contains the set vectors that may constitute the rays (only a subset of those is # is taken). We use this later on to check preservation of regularity during a flip. #---------------------------------------------------------------------------------- equ_simplicial_cone:=proc(a,simp) local i,directionineq,u,ineqs,temp,elemento,aux,facet,ecuacion: ineqs:={}: for elemento in simp do temp:=simp minus {elemento}: temp:=convert(temp,list): aux:=[]: for i from 1 to rowdim(a) do aux:=[op(aux),x.i]: od: u:=convert(aux,vector): facet:=stackmatrix(u,transpose(submatrix(a,[$1..rowdim(a)],temp))); ecuacion:=det(facet); directionineq:=signum(det(stackmatrix(transpose(submatrix(a,[$1..rowdim(a)],[elemento])),transpose(submatrix(a,[$1..rowdim(a)],temp))))); if directionineq=-1 then ineqs:=ineqs union {ecuacion+1<=0}: else ineqs:=ineqs union {-ecuacion+1<=0}: fi: od: ineqs; end: #---------------------------------------------------------------------------------------- # The following procedure has the power of checking whether a triangulation is regular. # We input the Gale transform and the triangulation. During the process # we produce the set of simplicial cones that are dual to the simplices of the triangulation LL. # The intersection of the relative interiors of a these simplicial cones is empty happens when the # triangulation is not regular. #---------------------------------------------------------------------------------------- check_regularity:=proc(B,LL) local cell,ecuaciones,elements,temp; elements:={$1..coldim(B)}: ecuaciones:={}: for cell in LL do temp:=elements minus cell: ecuaciones:=ecuaciones union equ_simplicial_cone(B,temp); od: feasible(ecuaciones); end: #---------------------------------------------------------------------------------------------- # Once we have detected a circuit and their two simplicial complexes ( in last # procedure deltaplus and deltaminus) we have to check to things in order for the # BISTELLAR move to occur between two triangulations: (1) One of them is a subsimplical # complex embedded in the triangulation. (2) Deltaminus and Deltaplus share a common link. # If these two conditions occur we can add the circuit in question to the set of # GOODCIRCUITS (those over which we can flip). So the input to this procedure is a # triangulation or simplicial complex. this is passed in as triangulation TRIAN. # The output is the set BISTELLARVECINOS that contains all the bistellar neighbors of the input triangulation. #----------------------------------------------------------------------------------------------- link_and_boundary_check:=proc(trian) local auxillary,i,j,k,deltaplus,deltaminus,aux,p,lista,flag,card,temp,link,initiallink,join1,join2,goodcircuits, bistellarneighbor,s,t,other,goody,tinier_neighbors,test,tri,V,D,cell: global flips,bistellarvecinos: #print(`Finding the neighbors of a triangulation`); tri:=trian: tinier_neighbors:={}: goodcircuits := []: bistellarvecinos:={}: for i from 1 to nops(flips) do deltaplus:=op(1,op(i,flips)): auxillary :=[]: p:= nops(deltaplus): for j from 1 to p do aux[j] := []: for k from 1 to nops(tri) do if (op(k,tri) union op(j,deltaplus))=op(k,tri) then aux[j]:=[op(aux[j]),op(k,tri)]: fi: od: auxillary := [op(auxillary),aux[j]]: od: flag := 0: for k from 1 to p do if nops(aux[k])>0 then flag := flag + 1 fi: od: if flag = p then goodcircuits:=[op(goodcircuits),[i,deltaplus,2,auxillary,1]]: # print(`circuit#`,i,`deltaplus=`,deltaplus): else # If the simplicial complex deltaplus is not part of the simplicial complex # we test now this question for deltaminus. deltaminus := op(2,op(i,flips)): auxillary :=[]: p := nops(deltaminus): for j from 1 to p do aux[j] := []: for k from 1 to nops(tri) do if (op(k,tri) union op(j,deltaminus))=op(k,tri) then # if the if condition holds the jth simplex of deltaplus (deltaminus) is contained in # the kth simplex of the triangulation tri. We add it to the list. # The jth element of auxillary contains the collection of simplices of the triangulation # that contain the kth simplex of delta(plus-minus) aux[j]:=[op(aux[j]),op(k,tri)]: fi: od: auxillary := [op(auxillary),aux[j]]: od: flag := 0: for k from 1 to p do if nops(aux[k])>0 then flag := flag + 1 fi: od: if flag = p then goodcircuits:=[op(goodcircuits), [i,deltaminus,1,auxillary,2]]: # The variable auxillary indicates containment of a maximal simplex of the simplicial # complex deltaplus (or deltaminus) inside the simplices of the triangulation. # print(`circuit#`,i,`deltaminus=`,deltaminus): fi: fi: od: # We have now a list of GOODCIRCUITS based on which we can perform BISTELLAR operations from the triangulation # trian for k from 1 to nops(goodcircuits) do goody :=op(k,goodcircuits): test :={}: D := convert(op(2,goody),list): p := nops(D): # D is a simplicial complex and p is the number of maximal simplices D has. for i from 1 to p do link[i] := {}: lista :=op(i,op(4,goody)): card := nops(lista): for j from 1 to card do temp := op(j,lista): temp := temp minus op(i,D): link[i] := link[i] union {temp}: od: test := test union link[1]: if not(test=link[i]) then i :=p+5 fi: # An important condition is that all the links of simplices of D should be # equal. od: if i < p+5 then tri := convert(tri,set): initiallink := convert(test,list): join1 := {}: for s from 1 to nops(initiallink) do for t from 1 to p do join1 := join1 union {op(s,initiallink) union op(t,D)} od: od: V := tri minus join1: other := convert(op(op(3,goody),op(op(1,goody),flips)),list): # print(`circuit along which you can flip:`); # print(`circuit#`,op(1,goody)); # print(`flip from `,D,`to`,other); # print(`link:=`,test); join2 := {}: for s from 1 to nops(initiallink) do for t from 1 to nops(other) do join2 := join2 union {op(s,initiallink) union op(t,other)} od: od: bistellarneighbor := V union join2: # print(`This is the triangulation I got`,bistellarneighbor); bistellarvecinos:=bistellarvecinos union {bistellarneighbor}: if nops(bistellarneighbor) >= nops(trian) then tinier_neighbors:=tinier_neighbors union {bistellarneighbor}: else # print(`BORING...at least one of my neighbors is smaller than me`); fi: fi: od: if nops(tinier_neighbors)=nops(bistellarvecinos) then print(`GREAT the triangulation`,trian,`is a local minima in the graph of flips`); print(`it has`, nops(trian),`many simplices`); print(`its neighbors have # of simplices:`); for cell in bistellarvecinos do print(nops(cell)); od: fi: bistellarvecinos; end: #----------------------------------------------------------------------------------- # The following procedures are auxiliary procedures for the procedure cheap_representative # The most important is the procedure codification_triangulation. In this procedure, # given a triangulation tri of a point configuration with n points, we obtain a unique label for it. # This is an order list which whose elements are the numbers of the simplices of # tri1. The simplex {i1,i2,...,id} receives the number s which is its lexicographic # position in the list of n choose d order lists (This position is calculated in the # procedure orden). #--------------------------------------------------------------------------------- orden:=proc(L,n) local d,s,k,l,i,j: s:=1: k:=1: d:=nops(L): for i from 1 to d do l:=L[i]: for j from k to l-1 do s:=s+binomial(n-j,d-i): od: k:=l+1: od: RETURN(s); end: codification_triangulation:=proc(tri,n) local lista,cell,temp,posicion; global orden: lista:=[]: for cell in tri do temp:=sort(convert(cell,list)): posicion:=orden(temp,n): lista:=[op(lista),posicion]: od: lista:=sort(lista); RETURN(lista); end: bigger:=proc(L1,L2) local i: if nops(L1)>nops(L2) then RETURN(true) elif nops(L2)>nops(L1) then RETURN(false); else for i from 1 to nops(L1) do if op(i,L1)>op(i,L2) then RETURN(true); elif op(i,L2)>op(i,L1) then RETURN(false); fi: od: RETURN(false); fi: end: finds_smallest:=proc(O) local temp,i,maximo: global bigger: maximo:=op(1,O): for i from 2 to nops(O) do temp:=op(i,O): if bigger(temp,maximo) then maximo:=temp: fi: od: RETURN(maximo); end: #------------------------------------------------------------------------------------------ # In this procedure we find a cheap representative of a equivalence class of triangulations. # The equivalence classes are classes modulo symmetries. For each member of the class we # get label with codification_triangulation and then find the smallest with a lexicographic # criteria #------------------------------------------------------------------------------------------ cheap_representative:=proc(tri1,q) local O1,O2: global symmetry_moves,orbita,gale,find_smallest,codification_triangulation: if nops(symmetry_moves)>1 then O1:=orbita(tri1): q:=nops(O1): O2:=map(xx -> codification_triangulation(xx,coldim(gale)),O1); RETURN(find_smallest(O2)): else q:=1: RETURN(codification_triangulation(tri1,coldim(gale))); fi: end: #----------------------------------------------------------------------------------- one_skeletonSec:=proc(trian) local flag,s,t,bistellarneighbor,vecinos,numero,position,gkztrian,gkzvalues,q,barato: global skeleton,gkzcoord,gkzordered,gkzvectors,looseends,totaltriang,triangulations,shadowstriang,cheap_representative: numero:=op(1,trian): vecinos:=link_and_boundary_check(op(2,trian)): for bistellarneighbor in vecinos do gkzordered:=sort(puntotes(bistellarneighbor)): gkztrian:=sort(puntitos(bistellarneighbor)): if member(gkztrian,gkzvectors) or member(gkzordered,gkzcoord) then flag:=0: barato:=cheap_representative(bistellarneighbor,'q'): if member(barato,shadowstriang,'position') then flag:=1: fi: # we add an edge to the 1-skeleton (=graph) of the secondary polytope if flag=0 then t:=nops(triangulations): addvertex({t+1},skeleton): addedge({numero,t+1},skeleton): totaltriang:=totaltriang+q: triangulations:=[op(triangulations),bistellarneighbor]: shadowstriang:=[op(shadowstriang),barato]: looseends:=looseends union {[t+1,bistellarneighbor]}: else if edges({numero,position},skeleton)={ } then addedge({numero,position},skeleton): fi: fi: else barato:=cheap_representative(bistellarneighbor,'q'); totaltriang:=totaltriang+q: t:=nops(triangulations): addvertex({t+1},skeleton): addedge({numero,t+1},skeleton): triangulations:=[op(triangulations),bistellarneighbor]: shadowstriang:=[op(shadowstriang),barato]: gkzcoord:=[op(gkzcoord),gkzordered]: gkzvectors:=[op(gkzvectors),gkztrian]: looseends:=looseends union {[t+1,bistellarneighbor]}: fi: od: looseends:=looseends minus {trian}: end: #-------------------------------------------------------------------------------------------- # # This procedure computes the 1-skeleton of the secondary polytope of a point configuration # #---------------------------------------------------------------------------- grafo:=proc(ma,folder) local auxi,r,w,s,seed,raiz,i,cell,rara,decision: global cocircuits,elements,circuitos,gale,gkzcoord, gkzvectors,looseends,nonregulars,retriangus,skeleton, symmetry_moves,totaltriang,triangulations,st,time,horas,hour,shadowstriang: interface(quiet=true): if not(nargs>2) then symmetry_moves:=[]: else print(`READING THE SYMMETRY GROUP OF THE CONFIGURATION`); read(args[3]): fi: gale:=matrix(convert(kernel(ma),list)): print(`I STARTED TO COMPUTE THE CIRCUITS OF THE POINT CONFIGURATION`); circo(gale): print(`I FINISH COMPUTING THE CIRCUITS WITH A TOTAL OF:`,nops(cocircuits)): circuitos:=cocircuits: raiz:=get_a_triangulation(circuitos): seed:={}: for cell in raiz do seed:=seed union {{$1..coldim(ma)} minus cell}: od: print(`THIS IS THE SEED TRIANGULATION FROM WHICH WE WILL GROW THE SECONDARY`,seed); print(`It contains`,nops(seed),`simplices`); flipper( ): # initialization of the 1-skeleton as the empty graph. skeleton:=void(1): totaltriang:=nops(orbita(seed)): chirotope(ma): gkzcoord:=[sort(puntotes(seed))]: gkzvectors:=[sort(puntitos(seed))]: nonregulars:={}: shadowstriang:=[cheap_representative(seed,'q')]: totaltriang:=q: triangulations:=[seed]: looseends:={[1,seed]}: print(`I AM GOING TO COMPUTE THE SECONDARY POLYTOPE`); st:=time(): hour:=0: rara:=time(): while looseends<>{} do one_skeletonSec(op(1,looseends)); lprint(nops(triangulations),`vertices have been produced so far`); lprint(totaltriang,`without symmetries`); lprint(`The set of those vertices waiting for the determination of its neighbors has `,nops(looseends),`elements`); lprint(); divide(time()-st,3600,'horas'): if horas-hour>1 then hour:=round(horas): save triangulations,folder: fi: od: print(time()-rara,`seconds is the time of generation of the graph given the circuits`); print(`We have generated`,nops(triangulations),`distinct triangulations using bistellar operations.`); print(`We still have to check which of these triangulations are regular`); print(`All the triangulations produced so far are stored in the file designated`); save triangulations,folder: # There is a final important check to do. Are the triangulations obtained from the flips regular ??. # we check this using the subroutine check_regularity (which amounts to solving a system of inequalities). print(); print(`totaltriangulations:`,totaltriang); print(`Would you like to determine which triangulations are regular? `); #decision:=readstat(`(YES=1/NO=0):`); #if decision=1 then for i from 1 to nops(triangulations) do auxi:=op(i,triangulations): if not(check_regularity(gale,auxi)) then print(`BINGO!! the following is a non-regular triangulation`,auxi); nonregulars:=nonregulars union {auxi}: fi: print(`We have analyzed`,i,`of the`,nops(triangulations),`vertices for regularity`); od: retriangus:=convert(triangulations,set) minus nonregulars: if nargs<3 then print(`THE SECONDARY POLYTOPE HAS`,nops(retriangus),`VERTICES.`); print(`The regular triangulations are in the variable retriangus`); print(`IF YOU WISH YOU MAY APPLY FURTHER ANALYSIS TO THE VERTICES`); lprint(): else print(`THE SECONDARY POLYTOPE HAS`,totaltriang-cuenta(nonregulars),`VERTICES.`); print(`The representatives of the regular triangulations are in the variable retriangus`); print(`IF YOU WISH YOU MAY APPLY FURTHER ANALYSIS TO THE VERTICES`); lprint(); fi: print(`recall that all the triangulations are stored in the variable triangulations and`); print(`the variable retriangus contains those that are regular. using the procedure`); print(`GKZ_embedding you can print the Gelfand-Kapranov-Zelevinsky embedding of the triangulations`); print(); print(`If you used the symmetries of the point configuration you have stored in the`); print(`variables triangulations and/or retriangus representatives for the distinct triangulations`); #fi: interface(quiet=false): end: #--------------------------------------------------------------------------------------- # The following procedure computes the Gelfand-Kapranov-Zelevinskii point associated to a triangulation. #--------------------------------------------------------------------------------------- puntotes:=proc(cell) local i,j,conta,listita,n,temp: global gale,x: n:=coldim(gale): interface(prettyprint=false): listita:=[]: for i from 1 to n do conta[i]:=0: for j from 1 to nops(cell) do if member(i,op(j,cell)) then temp:=sort(convert(op(j,cell),list)): conta[i]:=conta[i]+abs(x[op(temp)]): fi: od: listita:=[op(listita),conta[i]]: od: listita; end: puntitos:=proc(cell) local i,j,conta,listita,n: global gale: n:=coldim(gale): interface(prettyprint=false): listita:=[]: for i from 1 to n do conta[i]:=0: for j from 1 to nops(cell) do if member(i,op(j,cell)) then conta[i]:=conta[i]+1 fi: od: listita:=[op(listita),conta[i]]: od: listita; end: GKZ_embedding:=proc(U) local cell,D: D:=coldim(U): interface(quiet=true): interface(prettyprint=false): if nargs>1 then writeto(args[2]): else writeto(terminal): fi: lprint(`DIM=`,D); lprint(): lprint(`CONV_SECTION`); for cell in retriangus do lprint(puntotes(cell)); od: lprint(`END`): lprint(); writeto(terminal): interface(quiet=false): end: #-------------------------------------------------------------------------------------------- # In the following procedure we compute the orbit of a given triangulation. Previously # we have read the SYMMETRY_MOVES that encode the symmetry group of the point configuration. # The procedure consists simply in applying one by one the permutations (in a series of replacements). # The parameter Ti is the triangulation in question. #-------------------------------------------------------------------------------------------- muevepunto:=proc(numero) global permutacion: RETURN(op(numero,permutacion)); end: muevesimplejo:=proc(simplejo) local u: u:=map(x -> muevepunto(x),simplejo); if nops(u) muevesimplejo(x),triangulacion); RETURN(u); end: orbita:=proc(Ti) local j,newTi; global symmetry_moves,orbit,permutacion: orbit:={Ti}: for j from 1 to nops(symmetry_moves) do permutacion:=op(j,symmetry_moves): newTi:=muevetriangulacion(Ti): orbit:=orbit union {newTi}: od: RETURN(orbit); end: cuenta:=proc(lista) local cuantos,auxi: cuantos:=0: for auxi in lista do cuantos:=cuantos+nops(orbita(auxi)): od: cuantos; end: #-------------------------------------------------------------------------------- interface(quiet=false); U:=matrix([[1,1,1,1,1,0,0,0,0,0], [0,0,0,0,0,1,1,1,1,1], [1,0,0,0,0,1,0,0,0,0], [0,1,0,0,0,0,1,0,0,0], [0,0,1,0,0,0,0,1,0,0], [0,0,0,1,0,0,0,0,1,0], [0,0,0,0,1,0,0,0,0,1]]); grafo(U,outfile);