# hvector.g
#  By Russ Woodroofe
#  last updated Sept 2009
#
# Includes routine to compute
#  fvector, ftriangle
#  hvector, htriangle
#  pure k-skeleton
#  homotopy type, via automated free face reduction

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;

HvectorFromFvector:=function(fvector)
  local hvector, i, j, d;

  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;

Computehvector:=function(C)
  return HvectorFromFvector(Computefvector(C));
end;

# Find the pure s-skeleton of a simplicial complex
#  (this is the subcomplex generated by all faces of dimension s)
SCPureSkeleton:=function(C, s)

  local facet, l, Skel, Skel2;
  
  if not IsSimplicialComplex(C) then
    Info (InfoWarning, 1, "SCPureSkeleton: the argument should be a simplicial complex");
    return fail;
  fi;

  Skel:=Filtered(C, x->(Length(x)>=s+1));

  Skel2:=[];
  for facet in Skel do
    if Length(facet) > s+1 then
      RemoveSet(Skel, facet);
      UniteSet(Skel2, Combinations(facet, s+1));
    fi;
  od;

  # Handle lots of stuff at dim s and lots of stuff at dim higher than
  #  s both quickly
  if Length(Skel)>Length(Skel2) then
    UniteSet(Skel, Skel2);
    return Skel;
  else
    UniteSet(Skel2, Skel);
    return Skel2;
  fi;
end;

SubcomplexGeneratedByKFacets:=SCPureSkeleton;

ReduceFreeFace:=function(C, face)
  local facetscontaining, D, facet;

  D:=ShallowCopy(C);

  facetscontaining:=Filtered(D, x->IsSubset(x, face));
  if Size(facetscontaining) <> 1 then
    Print("Face not free: contained in facets ", facetscontaining);
    return fail;
  fi;

  RemoveSet(D, facetscontaining[1]);
  UniteSet(D, Filtered(Combinations(facetscontaining[1]), x->(not IsSubset(x, face))));
  for facet in D do
    D:=Filtered(D, x->(not IsSubset(facet, x)) or x=facet);
  od;

  return D;
end;

ReduceFreeFacesContaining:=function(C, minface)
  local facet, facelist, D, face, face2, reduced, facetscontaining, replfacets;

  D:=ShallowCopy(C);

  # get a list on non-maximal faces
  facelist:=[];
  for facet in C do
    UniteSet(facelist, Combinations(facet));
  od;
  SubtractSet(facelist, C);
  facelist:=Filtered(facelist, x->IsSubset(x, minface));

  reduced:=true;
  while reduced = true do
    reduced:=false;
    for face in facelist do
      facetscontaining:=Filtered(D, x->IsSubset(x, face));
      if Size(facetscontaining)=1 then
        reduced:=true;
	Print(face, "\t -> ", facetscontaining[1], "\n");

        # adjust list of facets
        RemoveSet(D, facetscontaining[1]);
        UniteSet(D, Filtered(Combinations(facetscontaining[1]), x->(not IsSubset(x, face))));
        for facet in D do
          D:=Filtered(D, x->(not IsSubset(facet, x)) or x=facet);
        od;

        # adjust list of faces
        facelist:=Filtered(facelist, x->(not IsSubset(x, face) and IsSubset(x, minface)));
        SubtractSet(facelist, D);
	if Size(face)=1 then
	  Print("    ", D, "\n");
	fi;
        break;
      fi;
    od;
  od;
  return D;
end;

ReduceAllFreeFaces:=function(C)
  local facet, facelist, D, face, face2, reduced, facetscontaining, replfacets;

  D:=ShallowCopy(C);

  # get a list on non-maximal faces
  facelist:=[];
  for facet in C do
    UniteSet(facelist, Combinations(facet));
  od;
  SubtractSet(facelist, C);

  reduced:=true;
  while reduced = true do
    reduced:=false;
    for face in facelist do
      facetscontaining:=Filtered(D, x->IsSubset(x, face));
      if Size(facetscontaining)=1 then
        reduced:=true;
	Print(face, "\t -> ", facetscontaining[1], "\n");

        # adjust list of facets
        RemoveSet(D, facetscontaining[1]);
        UniteSet(D, Filtered(Combinations(facetscontaining[1]), x->(not IsSubset(x, face))));
        for facet in D do
          D:=Filtered(D, x->(not IsSubset(facet, x)) or x=facet);
        od;

        # adjust list of faces
        facelist:=Filtered(facelist, x->(not IsSubset(x, face)));
        SubtractSet(facelist, D);
	if Size(face)=1 then
	  Print("    ", D, "\n");
	fi;
        break;
      fi;
    od;
  od;
  return D;
end;