RequirePackage("homology");

# Compute the f triangle by looking at the faces of d dimensional facets,
#   then d-1 diml facets, and so forth.
# This is reasonably fast for complexes of small dimension.
Computeftriangle3:=function(C)

  local checkedfaces, facestocheck, maxfacesize, i, facet, ftriangle, face;
  
  if not IsSimplicialComplex(C) then
    return fail;
  fi;
  
  maxfacesize:=SimplicialComplexDimension(C)+1;

  # Setup the ftriangle data structure, initialized to zeros  
  ftriangle:=[ [0] ];
  for i in [1.. maxfacesize] do
    Add(ftriangle, Concatenation(ftriangle[i], [0]));
  od;

  checkedfaces:=[];
  i:= maxfacesize;
  while i>0 do
    for facet in C do
      if Length(facet)=i then
        facestocheck:=Combinations(facet);
        for face in Difference(facestocheck, checkedfaces) do
          ftriangle[i+1][Length(face)+1]:=
             ftriangle[i+1][Length(face)+1] + 1;
        od;
        UniteSet(checkedfaces, facestocheck);
      fi;
    od;
    i:=i-1;
  od;
  
  return ftriangle;

end;

# Compute the f triangle by making a list of all faces, and checking
#   each facet to see what faces are included.  This might be faster than
#   the above for complexes of high dimension.
Computeftriangle2:=function(C)

  local facelist, deglist, newfaces, facet, facetdim, face, maxfacesize, pos, ftriangle, i;
  
  if not IsSimplicialComplex(C) then
    return fail;
  fi;
  
  facelist:=[];
  deglist:=[];
  maxfacesize:=0;
  
  # Get a sorted list of all the faces of C.  (While we’re at it, compute the maxsize).
  for facet in C do
    newfaces:=Combinations(facet);
    UniteSet(facelist, newfaces);
    maxfacesize:=Maximum(maxfacesize, Length(facet));
  od;
  
  # Go through all the faces again to find the maximum degree of each
  #  Note: we could probably avoid this by adding the maximum degree as a property of the face
  for facet in C do
    newfaces:=Combinations(facet);
    for face in newfaces do
      pos:=PositionSet(facelist, face);
      if IsBound(deglist[pos]) then
        deglist[pos]:=Maximum(deglist[pos], Length(facet));
      else
        deglist[pos]:=Length(facet);
      fi;
    od;
  od;
  
  # Construct the ftriangle data structure, initialize with zeros
  ftriangle:=[ [0] ];
  for i in [1.. maxfacesize] do
    Add(ftriangle, Concatenation(ftriangle[i], [0]));
  od;
  
  # Count faces of degree and size.
  for i in [1..Length(facelist)] do
    ftriangle[deglist[i]+1][Length(facelist[i])+1]:=
       ftriangle[deglist[i]+1][Length(facelist[i])+1] + 1;
  od;

  return ftriangle;
end;

# Probably slower than both the above, this was my first attempt.  Similar
#   in algorithm to Computeftriangle2, but Computeftriangle2 makes better
#   use of sorting to make searches quick.
Computeftriangle:=function(C)

  local facelist, deglist, i, j, k, newfaces, facetdim, complexdim,
        ftriangle;
  
  if not IsSimplicialComplex(C) then
    return fail;
  fi;
  
  # Record all faces of facets together with the size of the facet they are in
  facelist:= [];
  deglist:=[];
  complexdim:=0;
  for i in [1..Length(C)] do
    # for each facet, look at all its faces
    newfaces:=Subsets(C[i]);
    facetdim:=Length(C[i]);
    
    for j in [1..Length(newfaces)] do
      k:=Position(facelist, newfaces[j]);
      # if it's not in our list of faces, add it and record its degree
      #  otherwise, update its degree (if necessary)
      if k=fail then
        Add(facelist, newfaces[j]);
	Add(deglist, facetdim);
      else
        deglist[k]:=Maximum(deglist[k], facetdim);
      fi;
    od;
    complexdim:=Maximum(complexdim, facetdim);
  od;
  
  ftriangle:=[ [0] ];
  for i in [1..complexdim] do
    Add(ftriangle, Concatenation(ftriangle[i], [0]));
  od;
  
  for i in [1..Length(facelist)] do
    ftriangle[deglist[i]+1][Length(facelist[i])+1]:=
       ftriangle[deglist[i]+1][Length(facelist[i])+1] + 1;
  od;

  return ftriangle;
end;

# Change the Computeftriangle call to tweak the algorithm
Computehtriangle:=function(C)

  local ftriangle, htriangle, i, j, k;
  
  ftriangle:=Computeftriangle2(C);
  
  if ftriangle=fail then
    return fail;
  fi;
  
  htriangle:=[ [0] ];
  for i in [1..Length(ftriangle)-1] do
    Add(htriangle, Concatenation(htriangle[i], [0]));
  od;

  for i in [1..Length(ftriangle)] do
    for j in [1..i] do
      for k in [1..j] do
        htriangle[i][j]:=htriangle[i][j] 
	  + (-1)^(j-k) * Binomial( i-k, j-k) * ftriangle[i][k];
      od;
    od;
  od;
  
  return htriangle;
  
end;

Computefvector:=function(C)

  local facet, facelist, face, fvector, complexdim, i;

  if not IsSimplicialComplex(C) then
    return fail;
  fi;

  facelist:=[];
  complexdim:=0;
  for facet in C do
    UniteSet(facelist, Combinations(facet));
    complexdim:=Maximum(complexdim, Length(facet));
  od;
  
  fvector:=[];
  for i in [1..complexdim+1] do;
    fvector[i]:=0;
  od;
  
  for face in facelist do
    fvector[Length(face)+1]:=fvector[Length(face)+1]+1;
  od;
  
  return fvector;
end;

Computehvector:=function(C)

  local fvector, hvector, i, j, d;
  
  fvector:=Computefvector(C);
  if fvector=fail then 
    return fail; 
  fi;

  d:=Length(fvector); # d is actually size of biggest facet+1
  hvector:=[];
  for i in [1..d] do
    hvector[i]:=0;
  od;

  for i in [1..d] do
    for j in [1..i] do
      hvector[i]:=hvector[i] + (-1)^(i-j) * Binomial(d-j, d-i) * fvector[j];
    od;
  od;
  
  return hvector;
end;

# Compute the subcomplex given by the facets of dimension k
#   (size k+1)
SubcomplexGeneratedByKFacets:=function(C, k)

  local Ck, facet;

  if not IsSimplicialComplex(C) then
    return fail;
  fi;
  
  Ck:=[];
  for facet in C do
    if Length(facet)=k+1 then
      Add(Ck, facet);
    fi;
  od;
  
  return Ck;  
end;