#---------------PUNTOS----------------------
#
#NOTE!!! requires MAPLE V3
#
#
#--------------------------------------------


interface(quiet=true):
#############################################################################################
# PUNTOS is a collection of subroutines to study TRIANGULATIONS and SECONDARY POLYTOPES     #   
# for POINT CONFIGURATIONS. The main operation is the subroutine ``GRAFOSECONDARY'' 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 ``ACTIONSECONDARY'' #
# (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.                   #
#                                                                                           #
#                                                                                           #
# This package is a hybrid system, with the engine based on MAPLE. We require the help of   #
# MACAULAY to compute an initial triangulation. The mathematical principle behind this is   #
# the theory of Grobner bases of Toric Varieties. The default triangulation is a pulling    #
# triangulation, but can be chosen to be different by selecting a different weight vector   #
# that defines the monomial order or by reordering the variables. The procedure MAC_INPUT   #
# is in charge of computing the seed triangulation that is stored in an external file. One  #
# may also be interested in computing an inequality representation for the Secondary poly-  #
# tope. This is accomplished with the help of PORTA which offers transformation between re- #
# presentations of a polyhedron. The incorporation of the services of PORTA and MACAULAY is #
# made via unix-shell scripts and stream line editor commands. 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:
global x,cocircuits,parimpar,orbitasimplejo,orbitacirco:

chirotope(a): 
k:=rowdim(a):
n:=coldim(a):
T:={op(map(x->{op(x)},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}]:   
     cocircuits:=cocircuits union orbitacirco(cocircuit):
     print(`I have computed`,nops(cocircuits));
     T:=T minus orbitasimplejo(op(1,T));
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
    print(`facet`,op(j,cocircuits)):
    cuentalas:=cuentalas+1:
 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:
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])];
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:=stack(u,transpose(submatrix(a,[$1..rowdim(a)],temp)));
  ecuacion:=det(facet);
  directionineq:=signum(det(stack(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,test,tri,V,D:

global flips,bistellarvecinos:


tri:=trian:
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:
        bistellarvecinos:=bistellarvecinos union {bistellarneighbor}:
    fi:
od:
 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:
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);
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(`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`);
# 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 linear inequalities). 

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:
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-Zelevinskii embedding of the triangulations`);
print();
print(`Keep in mind that 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`);
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:

op(numero,permutacion);
end:

muevesimplejo:=proc(simplejo)
local u:
u:=map(x -> muevepunto(x),simplejo);
u;
end:

muevetriangulacion:=proc(triangulacion)
local u:
u:=map(x -> muevesimplejo(x),triangulacion);
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:
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);