// Example code from paper 14:
// _Colouring planar graphs_, by Paulette Lieby
//
print "Examples from: Colouring planar graphs";
ei := GetEchoInput();
SetEchoInput(true);
// Section 3: Planar graphs
G := CompleteGraph(3);
G := G join G;
AddEdge(~G, VertexSet(G) ! 1, VertexSet(G) ! 5);
AddEdge(~G, VertexSet(G) ! 3, VertexSet(G) ! 4);
assert IsKEdgeConnected(G, 2);
assert IsPlanar(G);
E := EdgeSet(G);
F := Faces(G);
F;
F[1];

Ds := [ [ 1 : x in [1..#F[i]] ] : i in [1..#F] ];
for i in [1..#F] do
  for j in [1..#F[i]] do
    if Index(InitialVertex(F[i][j])) gt Index(TerminalVertex(F[i][j])) then
      Ds[i][j] := -1;
    end if;
  end for;
end for;
Ds;

nstar := #F;
Gstar := MultiDigraph< nstar | >;
Bij := [ EdgeSet(Gstar) | ];
Fs := [ SequenceToSet(f) : f in F ];
for u in [1..nstar-1] do
  for v in [u+1..nstar] do
    M := Fs[u] meet Fs[v];
    for e in M do
      p := Position(F[u], e);
      q := Position(Edges(G), e);

      if Ds[u][p] eq 1 then
        Gstar, edge := AddEdge(Gstar, VertexSet(Gstar)!u, VertexSet(Gstar)!v);
        ChangeUniverse(~Bij, EdgeSet(Gstar));
        Bij[q] := edge;
      else
        Gstar, edge := AddEdge(Gstar, VertexSet(Gstar)!v, VertexSet(Gstar)!u);
        ChangeUniverse(~Bij, EdgeSet(Gstar));
        Bij[q] := edge;
      end if;

    end for;
  end for;
end for;
for i := 1 to Size(G) do
  "  ", EdgeSet(G).i, "  ", Bij[i];
end for;


// Section 5: Finding nowhere-zero k-flows
S := UnderlyingGraph(Gstar);
T, _, _, DFS := DepthFirstSearchTree(VertexSet(S) ! 1);
DFS;
EG := EdgeSet(Gstar);
VT := VertexSet(T);
visited := [ false : i in [1..#VT] ];
GDstar := MultiDigraph< Order(Gstar) | >;
for e in EG do
  u := InitialVertex(e);
  v := TerminalVertex(e);
  i := DFS[Index(u)];
  j := DFS[Index(v)];

  BE := false;
  if VT!u adj VT!v then
    /* is possibly a tree edge */

    if i lt j and not visited[j] then
      visited[j] := true;
      AddEdge(~GDstar, VertexSet(GDstar)!u, VertexSet(GDstar)!v);
    elif j lt i and not visited[i] then
      visited[i] := true;
      AddEdge(~GDstar, VertexSet(GDstar)!v, VertexSet(GDstar)!u);
    else
      /* must be a back edge (parallel edge) */
      BE := true;
    end if;
  else
    /* must be a back edge */
    BE := true;
  end if;

  if BE then
    if i lt j then
      AddEdge(~GDstar, VertexSet(GDstar)!v, VertexSet(GDstar)!u);
    else
      AddEdge(~GDstar, VertexSet(GDstar)!u, VertexSet(GDstar)!v);
    end if;
  end if;
end for;
assert Size(GDstar) eq Size(Gstar);


// Section 6: Testing for nowhere-zero k-flows
k := 3;
GDstar;
N := Network< Order(GDstar) | >;
phi := [ ];
for e in EdgeSet(GDstar) do
  u := VertexSet(N) ! EndVertices(e)[1];
  v := VertexSet(N) ! EndVertices(e)[2];
  N, newe := AddEdge(N, u, v, k-2);
  phi[Index(e)] := Index(newe);
end for;
N;
N +:= 2;
N;

s := VertexSet(N) ! (Order(N) - 1);
t := VertexSet(N) ! Order(N);
s, t;
for u in VertexSet(N) do
 if not u eq s and not u eq t then

  c := InDegree(u);
  AddEdge(~N, s, u, c);

  c := OutDegree(u);
  AddEdge(~N, u, t, c);

 end if;
end for;
N;
s := VertexSet(N) ! s;
t := VertexSet(N) ! t;
F := MaximumFlow(s, t);
for e in Edges(N) do
  if EndVertices(e)[1] eq s or EndVertices(e)[2] eq t then
    e, Capacity(e), Flow(e);
  end if;
end for;
E := EdgeSet(N);
for e in EdgeSet(GDstar) do
  e, Flow(E.phi[Index(e)]) + 1;
end for;
SetEchoInput(ei);