#############################################################################
##
#A  toric.gi                  toric library                      David Joyner
##
##  this file contains implementations for toric 
##
## 3-6-2006, changed Star to ToricStar, as a method instead of a 
##           global function
##

#############################################################################
##
#F  InDualCone(<v>,<L>)
##
##  returns true if v belongs to the dual of the code generated by
##  the vectors in L
##  
InstallGlobalFunction(InDualCone,function(v,L)
 local numgens,dim,i,V;
  TORIC.consistent_vectors(L);
  dim:=Size(v); 
  numgens:=Size(L); 
# V := Rationals^dim; 
  for i in [1..numgens] do 
   if TORIC.inner_product(v,L[i],IdentityMat(dim)) < 0 then return false; fi; 
  od; 
 return true; 
end);

#############################################################################
##
#F  DualSemigroupGenerators(<L>)
##
##  Input: L is a list of integral n-vectors 
##          generating a cone sigma 
##  Output: the generators of S_sigma,  
##
##  Idea:  let M be the max of the abs values of the 
##  coords of the L[i]'s, for each vector v in 
##  [1..M]^n, test if v in sigma^*. If so, 
##  add to list of possible generators. 
##  Once this for loop is finished, one can check 
##  this list for redundant generators. The  
##  trick below is to simply omit those elements 
##  which are of the form d1+d2, where d1 and d2 
##  are "small" elements in the integral dual cone. 
##  
InstallGlobalFunction(DualSemigroupGenerators,
function(L)
local d,i,dim,V,D,max,v,I,b,DpD,d1,d2,Dgens;
 TORIC.consistent_vectors(L);
 dim:=Size(L[1]); 
 V := Rationals^dim; 
 max:=TORIC.max_vectors(L); 
 D:=[]; 
 I:= Cartesian(List([1..dim],i->[-max..max])); 
 for v in I do 
  if InDualCone(v,L) then 
   Append(D,[v]); 
  fi; 
 od; 
 DpD:=[]; 
 for d1 in D do 
  if TORIC.inner_product(d1,d1,IdentityMat(dim)) <> 0 then 
  for d2 in D do 
   if TORIC.inner_product(d2,d2,IdentityMat(dim)) <> 0 then 
      Append(DpD,[d1+d2]); 
   fi; 
  od; 
  fi; 
 od; 
 Dgens:=[]; 
 for d in D do 
  if not(d in DpD) then 
     Append(Dgens,[d]); 
  fi; 
 od; 
 return Dgens; 
end);

#############################################################################
##
#F  EmbeddingAffineToricVariety(<L>)
##
##  Input: L is a list generating a cone 
##        (as in dual_semigp_gens),
##  Output: the toroidal embedding of X=affine_toric_variety(L)
##       (given as a list of multinomials)
##  
InstallGlobalFunction(EmbeddingAffineToricVariety,
function(L)
local map0,u,O,i,P1,P2,e,dim,V,R1,R2,hom;
 TORIC.consistent_vectors(L);
 u:=DualSemigroupGenerators(L);
 dim:=Size(u[1]);
 V:=Rationals^dim;
 O:=Zero(V);
 u:=Difference(u,[O]);
 map0:=[];
 for e in u do
  P1:=Filtered([1..dim],i->e[i] < 0);
  P2:=Filtered([1..dim],i->e[i] >= 0);
  if Size(P1) > 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i]))/Product(P1,i->X(Rationals,i+dim)^(-e[i])));
  fi;
  if Size(P1) = 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i])));
  fi;
 od;
 return map0;
end);

#############################################################################
##
#F  DivisorPolytope(<D>, <Rays>)
##
##  Input: D is the list of coeffs for the Weil divisor D
##  Rays is the list of vectors in the rays for the fan Delta
##   (Delta is determined by a list of generators L
##    but these are not explicitly needed here)
##  Output: the linear expressions in the affine coordinates of the
##   space of the cone which must be >0 for a point to be in the
##   desired polytope
##  
InstallGlobalFunction(DivisorPolytope,
function(D,Rays)
local n,Inequalities,R,i,j;
 TORIC.consistent_vectors(Rays);
 n:=Size(Rays[1]);
 Inequalities:=[];
 R:=PolynomialRing(Rationals,n);
 for i in [1..Size(Rays)] do
  Add(Inequalities,D[i]+Sum([1..n],j->X(Rationals,j)*Rays[i][j]));
 od;
 return Inequalities;
end);

#############################################################################
##
#F  PolytopeLatticePoints(<A>, <Perps>)
##
##  Input: Perps=[v1,...,vk] is the list of "inward normal" 
##            vectors perpendicular to the walls of the polytope P,
##         A=[a1,...,ak] is the amount wall is shifted from the 
##            origin (ai>=0)
##       Eg, x=0, x=a, y=0, y=b has Perps=[[1,0],[-1,0],[0,1],[0,-1]] 
##           and A=[0,a,0,b]
##  Output: the list of points in P \cap L^\perp
##  
InstallGlobalFunction(PolytopeLatticePoints,
function(A,Perps)
local R,n,M,Pts,i,j,Inequalities,L0,L1,L2,v,ineq;
 TORIC.consistent_vectors(Perps);
 R:=ShallowCopy(Perps);
 n:=Size(R[1]);
 M:=Maximum(List([1..Size(A)],i->AbsoluteValue(A[i])));
 L0:=List([1..n],i->[-M..M]);
 Pts:=Cartesian(L0);
 L2:=[];
 for v in Pts do
  Inequalities:=[];
  for i in [1..Size(R)] do
   ineq:=A[i]+Sum([1..n],j->v[j]*R[i][j]);
   Add(Inequalities,ineq);
  od;
  if Minimum(Inequalities) >= 0 then Add(L2,v); fi;
 od;
 return L2;
end);

#############################################################################
##
#F  DivisorPolytopeLatticePoints(<D>, <Cones>, <Rays_ordered> )
##
##  Input: Cones is the fan
##         Rays_ordered is the *ordered* list of rays for Cones
##         D is the list of coeffs for the Weil divisor D
##  Output: the list of points in P_D \cap L^\perp
##            which paramterize the elements in L(D)
##
##  Caution: The set of rays do not determine the cones in a fan.
##
InstallGlobalFunction(DivisorPolytopeLatticePoints,
function(D,Cones,Rays_ordered)
local R,n,M,Pts,i,j,Inequalities,L0,L1,L2,v,ineq,Rays;
 TORIC.consistent_vectors(Rays_ordered);
 Rays:=Rays_ordered;
 R:=ShallowCopy(Rays);
 n:=Size(R[1]);
 M:=Maximum(List([1..Size(D)],i->AbsoluteValue(D[i])));
 L0:=List([1..n],i->[-M..M]);
 Pts:=Cartesian(L0);
 L2:=[];
 for v in Pts do
  Inequalities:=List([1..Size(R)],i->D[i]+Sum([1..n],j->v[j]*R[i][j])); 
  if Minimum(Inequalities) >= 0 then Add(L2,v); fi;
 od;
 return L2;
end);

#############################################################################
##
#F  RiemannRochBasis(<D>, <Cones>, <Rays> )
##
##  Input: Cones is the fan
##         Rays_ordered is the *ordered* list of rays for Cones
##         D is the list of coeffs for the Weil divisor D
##
##  Output: computes a basis (a list of monomials)
##           for the RR space of the divisor represented by D
##
##
InstallGlobalFunction(RiemannRochBasis,
function(D,Cones,Rays)
local P,map0,d,R,e,i,P1,P2;
 TORIC.consistent_vectors(Rays);
 P:=DivisorPolytopeLatticePoints(D,Cones,Rays);
 map0:=[];
 d:=Size(Rays[1]);
 R:=PolynomialRing(Rationals,d); 
 for e in P do
  P1:=Filtered([1..d],i->e[i] < 0);
  P2:=Filtered([1..d],i->e[i] >= 0);
  if Size(P1) > 0 and Size(P2) > 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i]))/Product(P1,i->X(Rationals,i)^(-e[i])));
  fi;
  if Size(P1)= 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i])));
  fi;
  if Size(P2) = 0 then 
    Add(map0, 1/Product(P1,i->X(Rationals,i)^(-e[i])));
  fi;
 od;
 return map0;
end);

#############################################################################
##
#F  Faces( <Rays> )
##
##  Input: Rays is a list of rays for the fan Delta
##
##  Output: all the normals to the faces (hyperplanes of the cone)
##
##
InstallGlobalFunction(Faces,function(Rays)
local i,normals,n0,N,m,R,dim_sets,dim,v,pos,neg,n1,v0,N0;
 TORIC.consistent_vectors(Rays);
 normals:=[];
 R:=ShallowCopy(Rays);
 dim:=Size(R[1]);
 if dim<=1 then return 
   Error("\n Dimension of ambient space for the rays must be >1.\n");
 fi;
 v:=Sum(R);
 dim_sets:=Combinations(R,dim-1); ##### must have dim>1 here!
 for m in dim_sets do
  n0:=TORIC.normal_to_hyperplane(m);
  if TORIC.inner_product(v,n0,IdentityMat(dim))>0 then Add(normals,n0);
    else
    Add(normals,-n0);
  fi;
 od;
 Sort(normals);
 N0:=ShallowCopy(normals);
 N:=ShallowCopy(normals);
 for n1 in N do 
    if Minimum(List(Rays,v0->TORIC.inner_product(v0,n1,IdentityMat(dim))))<0 then
        N0:=Difference(N0,[n1]); 
    fi;
 od;        
 return Set(N0);
end);

#############################################################################
##
#F  InsideCone(<v>, <Rays> )
##
##  Input: Rays is a list of rays for the fan Delta
##
##  Output: true if v is in the cone represented by Rays
##          (a cone is a list of vectors generating it),
##          false otherwise
##
##  Idea: look at all dim-subsets of Rays and test if 
##        there is a positive solns vector all the normals 
##        to the faces (hyperplanes of the cone)
##
##
InstallGlobalFunction(InsideCone,function(v,Rays)
local w,F,dim,R;
 TORIC.consistent_vectors(Rays);
 R:=ShallowCopy(Rays);
 dim:=Size(R[1]);
 if Length(v)<>dim then
   Error("\n The vector ",v," must have the same dimension as the fan.\n");
 fi;
 if dim<=1 then return false; fi;
 F:=Faces(R);
 if Minimum(List(F,w->TORIC.inner_product(v,w,IdentityMat(dim))))<=0 then 
  return false; 
 fi;
 return true;
end);

#############################################################################
##
#F  NumberOfConesOfFan( <Delta>, <k> )
##
##  Input: Delta is the fan of cones
##         k is the dimension of the cones counted
##  Output: number of k-diml cones in the fan
##          (each hyperplane of one cone should not form the face
##           of any other cone in the fan)
##
##  Idea:   the fan Delta is represented as a set of max cones.
##          for each max cone, look at the k-diml faces
##          obtained by taking n choose k subsets of the
##          rays describing the cone. 
##          **Certain** of these k-subsets yield the
##          desired cones
##
InstallGlobalFunction(NumberOfConesOfFan,
function(Cones,k)
local sigma,R,Rays,n,i,j,C,C0,dim,ell,k_cones,sigma0,N,v0,v;
 if k=0 then return 1; fi;
 C0:=ShallowCopy(Cones);
 Rays:=TORIC.flatten(C0);
 R:=ShallowCopy(Rays);
 if k=1 then return Size(R); fi;
 k_cones:=[];
 dim:=Size(R[1]);
 if k>dim then Print("dimension error\n"); return 0; fi;
 N:=Size(C0);
 if k=dim then return N; fi;
 for i in [1..N] do
  C:=ShallowCopy(Cones[i]);
  sigma:= Combinations(C,k);
  for ell in [1..Size(sigma)] do
     v0:=Sum(sigma[ell]);
     if not(InsideCone(v0,C)) then  Add(k_cones,sigma[ell]); fi;
  od;
 od;
 return Size(Set(k_cones));
end);

#############################################################################
##
#F  ConesOfFan( <Delta>, <k> )
##
##  Input: Delta is the fan of cones
##         k is the dimension of the cones taken
##  Output: the set of k-diml cones in the fan
##          (each hyperplane of one cone should not form the face
##           of any other cone in the fan, ie,
##           only takes those subcones which are *not* interior)
##
##
InstallGlobalFunction(ConesOfFan,
function(Cones,k)
local sigma,R,Rays,n,x,i,j,C,C0,dim,ell,k_cones,sigma0,N,v0,v;
 C0:=ShallowCopy(Cones);
 if k=0 then return 0*C0[1][1]; fi;
 Rays:=TORIC.flatten(C0);
 R:=ShallowCopy(Rays);
 if k=1 then return List(R,x->[x]); fi;
 k_cones:=[];
 dim:=Size(R[1]);
 if k>dim then Print("dimension error\n"); return 0; fi;
 N:=Size(C0);
 if k=dim then return C0; fi;
 for i in [1..N] do
  C:=ShallowCopy(Cones[i]);
  sigma:= Combinations(C,k);
  for ell in [1..Size(sigma)] do
      v0:=Sum(sigma[ell]);
   if not(InsideCone(v0,C)) then  Add(k_cones,sigma[ell]); fi;
  od; ##excludes the "interior" cones
 od;
 return Set(k_cones);
end);

#############################################################################
##
#F  ToricStar( <cone>, <Cones> )
##
##  Input: Cones is the fan of cones
##         <cone> is a cone of <Cones>
##  Output: the star of the cone in a fan
##          ie, the cones tau which have <cone> as a face
##
##
InstallGlobalFunction(ToricStar,
function(cone,Cones)
local T,tau,Taus,i,j,sigma,C0,R,Rays,dim;
 sigma:=ShallowCopy(cone);
 C0:=ShallowCopy(Cones);
 Rays:=TORIC.flatten(C0);
 R:=ShallowCopy(Rays);
 dim:=Size(R[1]);
 Taus:=[];
 for j in [1..dim] do
    tau:=ConesOfFan(C0,j);
    Append(Taus,tau);
 od;
 T:=[];
 for tau in Taus do
    if TORIC.is_sublist(sigma,tau) then Append(T,[tau]); fi;
 od;
 return T;
end);

#############################################################################
##
#F  BettiNumberToric( <Cones>, <k> )
##
##  Input: Cones is the fan of cones
##         k is an integer
##  Output: the k-th betti number of the toric variety 
##          X(Delta), where Delta is a fan defined by its maximal
##          cones, Cones
##
##
InstallGlobalFunction(BettiNumberToric,
function(Cones,k0)
local i,j,n,a,C,k;
 if IsOddInt(k0) then return 0; fi;
 k:=Int(k0/2);
 C:=ShallowCopy(Cones)[1][1];
 n:=Size(C);
 a:=Sum( [k..n],i->(-1)^(i-k)*Binomial(i,k)*NumberOfConesOfFan(Cones,n-i));
 return a;
end);

#############################################################################
##
#F  EulerCharacteristic( <Cones> )
##
##  Input: Cones is the fan of cones
##  Output: Euler char of X(Delta), where 
##          Delta is a fan defined by its maximal cones, Cones
##
##  Note: X(Delta) **must be non-singular***
##
## was InstallGlobalFunction(EulerCharacteristic,[IsList],
InstallMethod(EulerCharacteristic,[IsList],
function(Cones)
local sigma,R,Rays,C0,dim;
 C0:=ShallowCopy(Cones);
 Rays:=TORIC.flatten(C0);
 R:=ShallowCopy(Rays);
 dim:=Size(R[1]);
 return NumberOfConesOfFan(Cones,dim);
end);

#############################################################################
##
#F  CardinalityOfToricVariety( <Cones>, <q> )
##
##  Input: Cones is the fan of cones
##  Output: the number of elements of X(Delta) over GF(q), where 
##          Delta is a fan defined by its maximal cones, Cones,
##          and where q is a prime power
##
##  Note: X(Delta) **must be non-singular***
##
InstallGlobalFunction(CardinalityOfToricVariety,
function(Cones,q)
local i,j,n,a;
 n:=Size(Cones[1][1]);
 a:=Sum([0..n],i->(q-1)^(i)*NumberOfConesOfFan(Cones,n-i));
 return a;
end);