##                RibGraph.g                   ##
##                                             ##
##          author:  Shaun Van Ault            ##
##          version date:  7/15/06             ##
##  This module collects some common functions ##
## associated with unsigned ribbon graphs      ##
## It is by no means complete or completed!    ##
## Updates will be forthcoming.                ##
#################################################


## Ribbon Graphs are represented as:
##  G = [V_1, V_2, ..., V_n], where each
##  V_i represents a vertex.
##  V_i = [e_1, e_2, ..., e_m] is a list of
##  edges attached to V_i, in cyclic order,
##  clockwise around the vertex.
##    If G is a non-orientable ribbon graph,
##  denote the edge(s) with a twist by a negative.

## RibbonGraph constructor
RibbonGraph := function(G)
  local e, v, edge, c, i, j, L, bc, K, done, compDone,
        color, p, q, p2, q2, B, brk, dir;

  ## Check if valid ribbon graph description
  if not IsList(G) then Error("Invalid ribbon graph"); fi;
  for i in [1..Size(G)] do
    if not IsList(G[i]) then Error("Invalid ribbon graph"); fi;
    for j in [1..Size(G[i])] do
      if not IsInt(G[i][j]) then Error("Invalid ribbon graph"); fi;
    od;
  od;

  ## Count vertices
  v := Size(G);

  ## Count edges
  L := [];
  for i in [1..v] do
    for j in [1..Size(G[i])] do
      c := AbsInt(G[i][j]);
      if IsBound(L[c]) then
        if not SignInt(L[c]) = SignInt(G[i][j]) then
          Error("Invalid ribbon graph");
        fi; 
        L[c] := 2*L[c];
      else
        if G[i][j] >= 0 then
          L[c] := 1;
        else
          L[c] := -1;
        fi;
      fi;
    od;
  od;
  e := 0;
  for i in [1..Size(L)] do
    if IsBound(L[i]) then
      if not AbsInt(L[i]) = 2 then Error("Invalid ribbon graph"); fi;
      e := e + 1;
    fi;
  od;

  ## Count connected components
  K := [];
  color := [1..v]*0;
  c := 1;
  K[1] := c;
  color[1] := 1;
  
  done := false;
  while not done do
    compDone := false;
    while not compDone do
      i := -1;
      for p in [1..v] do
        if color[p] = 1 then
          i := p;
          break;
        fi;
      od;
      if i < 0 then compDone := true;
      else
        for edge in G[i] do
          for j in [2..v] do
            if color[j] = 0 and edge in G[j] then
              K[j] := c;
              color[j] := 1;
            fi;
          od;
        od;
        color[i] := 2;
      fi;
    od;
    
    done := true;
    for p in [1..v] do
      if color[p] = 0 then
        done := false;
        c := c + 1;
        color[p] := 1;
        K[p] := c;
        break;
      fi;
    od;
  od;

  ## Count connected components of boundary
  B := StructuralCopy(G);
  bc := 0;
  for i in [1..v] do
    if IsEmpty(G[i]) then
      bc := bc + 1;
    else
      B[i] := B[i]*0;
    fi;
  od;

  ## B[i][j] = 0  ==> edge not visited
  ## B[i][j] = 1  ==> edge visited and terminated here
  ## B[i][j] = -1 ==> edge visited and originated here
  ## B[i][j] = 2  ==> edge traversed both ways
  done := false;
  brk := false;
  dir := 1;
  ##  dir = 1  ==> clockwise around vertices
  ##  dir = -1 ==> counterclockwise around vertices
  while not done do
    done := true;
    for i in [1..v] do
      for j in [1..Size(B[i])] do
        if B[i][j] = 0 or B[i][j] = dir then
          done := false;
          p := i;
          q := j;
          brk := true;
          break;
        fi;
      od;
      if brk then brk := false; break; fi;
    od;
    if not done then
      bc := bc + 1;
      if B[p][q] = 0 then
        B[p][q] := -dir;
      else
        B[p][q] := 2;
      fi;

      p2 := p;
      q2 := q;

      compDone := false;
      while not compDone do
        ## Traverse edge
        for i in [1..v] do
          for j in [1..Size(G[i])] do
            if not ((i = p2) and (j = q2)) then
              if G[i][j] = G[p2][q2] then
                p2 := i;
                q2 := j;
		if G[p2][q2] < 0 then
                  dir := -dir;
                fi;
                brk := true;
                break;
              fi;
            fi;
          od;
          if brk then brk := false; break; fi;
        od;
        if B[p2][q2] = 0 then
          B[p2][q2] := dir;
        elif B[p2][q2] = -dir then
          B[p2][q2] := 2;
        else
          Print("Undefined boundary components\n");
          Print(B);Print("\n");
          Print([p2,q2]);Print("\n");
          Error();
        fi;

        ## Traverse arc on vertex
        if dir > 0 then
          if q2 < Size(G[p2]) then
            q2 := q2 + 1;
          else
            q2 := 1;
          fi;
        else
          if q2 > 1 then
            q2 := q2 - 1;
          else
            q2 := Size(G[p2]);
          fi;
        fi;
        if B[p2][q2] = 0 then
          B[p2][q2] := -dir;
        elif B[p2][q2] = dir then
          B[p2][q2] := 2;
        elif p2 = p and q2 = q then
          compDone := true;
        else
          Print("Undefined boundary components\n");
          Print(B);Print("\n");
          Print([p2,q2]);Print("\n");
          Error();
        fi;

      od;  
    fi;
  od;

  return rec(graph := StructuralCopy(G),
             v := v,
             e := e,
             k := c,
             r := v - c,
             n := e - (v - c),
             bc := bc);
end;

BollobasRiordanPoly := function(G)
  local F, v, i, j, L, x, y, z, R, done, edges;

  x := X(Rationals, "x");
  y := X(Rationals, "y");
  z := X(Rationals, "z");
  R := 0;

  if not IsRecord(G) then
    G := RibbonGraph(G);
  fi;

  edges := AsSet(Concatenation(G.graph));
  L := [1 .. G.e]*0;
  done := false;

  while not done do
    F := StructuralCopy(G.graph);
    for i in [1 .. G.e] do
      if L[i] = 1 then
        for j in [1 .. G.v] do
          v := F[j];
          while edges[i] in v do
            RemoveElmList(v, Position(v, edges[i]));
          od;
          F[j] := v;
        od;
      fi;
    od;
    F := RibbonGraph(F);
    R := R + x^(G.r - F.r) * y^(F.n) * z^(F.k - F.bc + F.n);

    i := G.e;
    while L[i] = 1 do
      i := i - 1;
      if i = 0 then 
        done := true;
        break;
      fi;
    od;

    if i > 0 then
      L[i] := 1;
      for j in [i+1 .. G.e] do
        L[j] := 0;
      od;
    fi;
  od;

  return rec(poly := R, vars := [x,y,z]);
end;