SetGraphExample := function(v,e) // v: number of vertices, e: number of edges
 edges := [];
 for j in [1..e] do
  edges cat:= [<Random([1..v]),j,Random([1..v])>];
 end for;
 return <[i : i in [1..v]],edges>;
end function;

CyclicGraph := function(n) // returns cyclic graph with n edges
 return <[i : i in [1..n]],[<i,i,i+1> : i in [1..n-1]] cat [<n,n,1>]>;
end function;

C := CyclicGraph; // %%

DirectGraph := function(n)
 return <[i : i in [0..n]],[<i,i,i+1> : i in [0..n-1]]>;
end function;

D := DirectGraph; // %%

IsThin := function(G)
 E := {<e[1],e[3]> : e in G[2]};
 return #E eq #G[2];
end function;

IsGraphMorphism := function(m,G,H) // m = <v,e>
 v := m[1];
 e := m[2];
 // test, if v is a map:
 if not SequenceToMultiset([x[1] : x in v]) eq SequenceToMultiset(G[1]) then
  return false;
 end if;
 if not SequenceToSet([x[2] : x in v]) subset SequenceToSet(H[1]) then
  return false;
 end if;
 // test, if e is a map:
 if not SequenceToMultiset([x[1] : x in e]) eq SequenceToMultiset(G[2]) then
  return false;
 end if;
 if not SequenceToSet([x[2] : x in e]) subset SequenceToSet(H[2]) then
  return false;
 end if;
 // test, if v and e are compatible concerning source and target:
 for x in e do
  if not &and[<x[1][1],x[2][1]> in v, <x[1][3],x[2][3]> in v] then
   return false;
  end if;
 end for;
 return true;
end function;

IsRightUnique := function(u) // u: relation, e.g. u := [<1,3>,<1,4>,<2,3>],
                             // or u := {<1,3>,<1,4>,<2,3>}
 right_unique := true;
 left_elements := {x[1] : x in u};
 for y in left_elements do
  if #{x[2] : x in u | x[1] eq y} ge 2 then
   right_unique := false;
   break y;
  end if;
 end for;
 return right_unique;
end function;

CompletionsToMaps := function(D,C,u)
 // D/C: list of elements in the domain/codomain,
 // u: right unique relation
 // e.g. D := [1,2,3,4]; C := [1,2,3,4,5]; u := {<1,3>, <3,5>};
 to_be_mapped := [x : x in D | not x in {y[1] : y in u}];
 list := [u];
 for i in to_be_mapped do
  list_new := [];
  for v in list do
   for j in C do
    list_new cat:= [v join {<i,j>}];
   end for;
  end for;
  list := list_new;
 end for;
 return list;
end function;

RelationOnVerticesFromPartialMapOnEdges := function(x);
 // x: partial map on edges from graph G to graph H (G, H not required as data)
 return {<z[1][1],z[2][1]> : z in x} join {<z[1][3],z[2][3]> : z in x};
end function;

RVPME := RelationOnVerticesFromPartialMapOnEdges; // %%

ListGraphMorphisms := function(G,H)
 list := [[]];
 for z in G[2] do
  list_new := [];
  for w in H[2] do
   for x in list do
    x_test := x cat [<z,w>];
    if IsRightUnique(RVPME(x_test)) then // %%
     list_new cat:= [x_test];
    end if;
   end for;
  end for;
  list := list_new;
 end for;
 list_mor := [];
 for y in list do
  list_completions_to_maps_on_vertices := [Sort(SetToSequence(x)) : x in
    CompletionsToMaps(G[1],H[1],RVPME(y))]; // %%
  list_mor cat:= [<x,y> : x in list_completions_to_maps_on_vertices];
 end for;
 return list_mor;
end function;

Identity := function(G);
 return <[<x,x> : x in G[1]], [<x,x> : x in G[2]]>;
end function;

Is_Injective := function(m,G,H)
 // m: G -> H: <morphism_on_vertices,morphism_on_edges>
 return #SequenceToSet([u[2] : u in m[1]]) eq #m[1] and
        #SequenceToSet([u[2] : u in m[2]]) eq #m[2];
end function;

Is_Surjective := function(m,G,H)
 // m: G -> H: <morphism_on_vertices,morphism_on_edges>
 return #SequenceToSet([u[2] : u in m[1]]) eq #H[1] and
        #SequenceToSet([u[2] : u in m[2]]) eq #H[2];
end function;

Is_Bijective := function(f,G,H)
 return Is_Injective(f,G,H) and Is_Surjective(f,G,H); // %%
end function;

ComposeGraphMorphisms := function(p,q)
 v := []; // v for vertices
 for x in p[1] do
  v cat:= [<x[1],y[2]> : y in q[1] | x[2] eq y[1]];
 end for;
 e := []; // e for edges
 for x in p[2] do
  e cat:= [<x[1],y[2]> : y in q[2] | x[2] eq y[1]];
 end for;
 return <v,e>;
end function;

IsIsomorphic := function(G,H) // G, H: graphs
 if #G[1] ne #H[1] or #G[2] ne #H[2] then
  return <0,false>;
 end if;
 list := ListGraphMorphisms(G,H); // %%
 for m in list do
  if Is_Bijective(m,G,H) then // %%
   return <m,true>;
  end if;
 end for;
 return <0,false>;
end function;

IsSubgraph := function(G,H) // returns true if G is a subgraph of H
 for v in G[1] do
  if not v in H[1] then
   return <0,false>;
  end if;
 end for;
 for e in G[2] do
  if not e in H[2] then
   return <0,false>;
  end if;
 end for;
 m := <[<v,v> : v in G[1]],[<e,e> : e in G[2]]>;
 return <m,true>;
end function;

IsFullSubgraph := function(G,H);
 if not IsSubgraph(G,H)[2] then // %%
  return false;
 end if;
 return not &or[e[1] in G[1] and e[3] in G[1] and not e in G[2] : e in H[2]];
end function;

VtoE := function(G,H,Vf) // H thin, Vf = [<1,2>,<2,5>,<3,1>,<4,2>] map on
                         // vertices
 for e in G[2] do
  if #[h : h in H[2] | h[1] eq [v[2] : v in Vf | v[1] eq e[1]][1] and
     h[3] eq [v[2] : v in Vf | v[1] eq e[3]][1]] eq 0 then
   print "graph morphism does not exist";
   return <0,0>;
  end if;
 end for;
 Ef := [<e,[h : h in H[2] | h[1] eq [v[2] : v in Vf | v[1] eq e[1]][1] and
         h[3] eq [v[2] : v in Vf | v[1] eq e[3]][1]][1]> : e in G[2]];
 return <Vf,Ef>;
end function;

ListGraphMorphisms_partial := function(f,G,H)
 // returns all graph morphisms that obey given partial mapping rule f
 vertices_partial := f[1];
 v1 := [n[1] : n in f[1]];
 edges_partial := f[2];
 e1 := [n[1] : n in f[2]];
 list := [[]];
 for z in G[2] do
  list_new := [];
  if z in e1 then
   H_edges := [n[2] : n in f[2] | n[1] eq z];
  else
   H_edges := H[2];
  end if;
  source := [v : v in H[1] | #[e : e in H[2] | e[1] eq v] ge 1];
  target := [v : v in H[1] | #[e : e in H[2] | e[3] eq v] ge 1];
  if z[1] in v1 then
   source := [n[2] : n in f[1] | n[1] eq z[1]];
  end if;
  if z[3] in v1 then
   target := [n[2] : n in f[1] | n[1] eq z[3]];
  end if;
  H_edges := [h : h in H_edges | h[1] in source and h[3] in target];
  for w in H_edges do
   for x in list do
    x_test := x cat [<z,w>];
    if IsRightUnique(RVPME(x_test)) then // %%
     list_new cat:= [x_test];
    end if;
   end for;
  end for;
  list := list_new;
 end for;
 list_mor := [];
 for y in list do
  list_completions_to_maps_on_vertices := [Sort(SetToSequence(x)) : x in
    CompletionsToMaps(G[1],H[1],RVPME(y)
    join SequenceToSet(f[1]))]; // %%
  list_mor cat:= [<x,y> : x in list_completions_to_maps_on_vertices];
 end for;
 return list_mor;
end function;

DisjointUnionCycles := function(list) // e.g. list := [2,4,5]
 // returns $\C_2\sqcup\C_4\sqcup\C_5$
 G1 := [i : i in [1..&+list]]; // vertices
 lists := [[u : u in [1..list[1]]]] cat [[u : u in
  [&+[list[k] : k in [1..i-1]]+1..&+[list[k] : k in [1..i]]]] : i in [2..#list]];
 edges := &cat[[<t[i],t[i+1]> : i in [1..#t-1]] cat [<t[#t],t[1]>] : t in lists];
 edges_numbered := [<edges[i][1],i,edges[i][2]> : i in [1..#edges]];
 G := <G1,edges_numbered>;
 return G;
end function;

DUC := DisjointUnionCycles; // %%

DCN := function(n)
 G := DUC([n,n]); // %%
 H := C(n); // %%
 f := VtoE(G,H,Sort([<i,i mod n> : i in [1..2*n-1] | not i eq n] cat [<n,n>,<2*n,n>]));
 return <G,H,f>;
end function;

RedSeq := function(S); // S: sequence, to be reduced and sorted
 return Sort(SetToSequence(SequenceToSet(S)));
end function;

Equivclasses := function(R,M) // R: relation on set M
 Rinv := [<r[2],r[1]> : r in R];
 Diag := [<m,m> : m in M];
 RR := [r : r in R cat Rinv | not r[1] eq r[2]];
 equivclasses := [];
 Mtodo := [m : m in M];
 while not #Mtodo eq 0 do
  k := Mtodo[1];
  kclassold := [];
  kclassnew := [k];
  while not #kclassnew eq #kclassold do
   kclassold := kclassnew;
   kclassnew cat:= [j[2] : j in RR | j[1] in kclassold];
   kclassnew := RedSeq(kclassnew); // %%
  end while;
  equivclasses cat:= [kclassnew];
  Mtodo := [u : u in Mtodo | not u in kclassnew];
 end while;
 return equivclasses;
end function;

Equivrelation := function(R,M) // R: relation on set M
 equivclasses := Equivclasses(R,M); // %%
 return Sort(&cat[[<k,l> : k, l in x] : x in equivclasses]);
end function;

DisjointUnionSets := function(X,Y); // X, Y lists
 return [<1,x> : x in X] cat [<2,y> : y in Y];
end function;

PushoutSets := function(X,Y,X2,f,g); // f : X -> Y, g : X -> X2 maps
 M := DisjointUnionSets(X2,Y); // %%
 R := [ [<<1, g_elt[2]>, <2, f_elt[2]>> : g_elt in g, f_elt in f
          | g_elt[1] eq x and f_elt[1] eq x][1] : x in X];
 equiv := Equivclasses(R,M); // %%
 u := [ [<x2,t> : t in equiv | <1,x2> in t][1] : x2 in X2];
 v := [ [< y,t> : t in equiv | <2, y> in t][1] :  y in  Y];
 return <equiv, u, v>; // pushout, u : X2 -> pushout, v : Y -> pushout
end function;

PushoutGraphs := function(X,Y,X2,f,g); // f : X -> Y, g : X -> X2
 // graph morphisms, returns the pushout
 vertices := PushoutSets(X[1],Y[1],X2[1],f[1],g[1]); // %%
 edges := PushoutSets(X[2],Y[2],X2[2],f[2],g[2]); // %%
 N := [i : i in [1..#vertices[1]]];
 E := [i : i in [1..#edges[1]]]; // edges without source and target
 EE := [ <Index(vertices[1], [n : n in vertices[1] | <edges[1][e][1][1],
       edges[1][e][1][2][1]> in n][1]),e, Index(vertices[1], [n :
      n in vertices[1] | <edges[1][e][1][1],edges[1][e][1][2][3]> in n][1])> :
      e in E];
      // edges with source and target
 PP := <N,EE>;
 uN := [<x2, Index(vertices[1], [n[2] : n in vertices[2] | x2 eq n[1]][1]) > :
        x2 in X2[1]]; // <p,p@v>
 uE := [<x2, Index(edges[1], [e[2] : e in edges[2] | x2 eq e[1]][1]) > :
        x2 in X2[2]]; // <p,p@v>, second entry is number of edge
 uEE := [ <x[1], <[ee[1] : ee in EE | ee[2] eq x[2]][1],x[2],[ee[3] : ee in EE |
      ee[2] eq x[2]][1]>> : x in uE]; // <p,p@u>
      // second entry with source and target
 u := <uN,uEE>;
 vN := [<y, Index(vertices[1], [n[2] : n in vertices[3] | y eq n[1]][1]) > :
        y in Y[1]]; // <p,p@v>
 vE := [<y, Index(edges[1], [e[2] : e in edges[3] | y eq e[1]][1]) > :
        y in Y[2]]; // <p,p@v>, , second entry is number of edge
 vEE := [ <x[1], <[ee[1] : ee in EE | ee[2] eq x[2]][1],x[2],[ee[3] : ee in EE |
          ee[2] eq x[2]][1]>> : x in vE]; // <p,p@u>
          // second entry with source and target
 v := <vN,vEE>;
 return <PP,u,v>; // pushout, u : X2 -> pushout, v : Y -> pushout
end function;

PullbackSets := function(X,Y,Y2,f,g) // f: X -> Y, g: Y2 -> Y, maps between sets
 // given as a list of tuples, e. g. f = [<x_1,y_1>,...,<x_n,y_n>]
 P := Sort([<x,y2> : x in X, y2 in Y2 | #[<s,t> : s in f, t in g | x eq s[1] and
  y2 eq t[1] and s[2] eq t[2]] eq 1]);
 return P;
end function;

PullbackGraphs := function(X,Y,Y2,f,g)
 // f: X -> Y, g: Y2 -> Y graph morphisms
 vertices := PullbackSets(X[1],Y[1],Y2[1],f[1],g[1]); // %%
 edges := PullbackSets(X[2],Y[2],Y2[2],f[2],g[2]); // %%
 N := [i : i in [1..#vertices]];
 E := [i : i in [1..#edges]]; // edges without source and target
 EE := [<Index(vertices,<edges[e][1][1],edges[e][2][1]>),e,
         Index(vertices,<edges[e][1][3],edges[e][2][3]>)> : e in E];
         // edges with source and target
 PP := <N, EE>;
 vE := [<ee,edges[ee[2]][1]> : ee in EE]; // <p,p@v>
 uE := [<ee,edges[ee[2]][2]> : ee in EE]; // <p,p@u>
 vN := [<n,vertices[n][1]> : n in N]; // <p,p@v>
 uN := [<n,vertices[n][2]> : n in N]; // <p,p@u>
 v := <vN,vE>;
 u := <uN,uE>;
 return <PP,u,v>; // pullback, u: pullback -> Y2, v: pullback -> X
end function;

InducedMorphismSetsPO := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 T := <Y2,u2,v2>;
 P := PushoutSets(X,Y,X2,v,u); // %%
 x2 := [<i,[t[2] : t in T[2] | t[1] in [r[1] : r in P[2] | r[2] eq i]]> :
        i in P[1]];
 y := [<i,[t[2] : t in T[3] | t[1] in [r[1] : r in P[3] | r[2] eq i]]> :
       i in P[1]];
 c := [<r[1],r[2][1]> : r in x2 | not #r[2] eq 0] cat [<r[1],r[2][1]> :
       r in y | not #r[2] eq 0];
 return RedSeq(c); // %%
end function;

InducedMorphismGraphsPO := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 T := <Y2,u2,v2>;
 P := PushoutGraphs(X,Y,X2,v,u); // %%
 x2_vertices := [<i,[t[2] : t in T[2][1] | t[1] in [r[1] : r in P[2][1] |
                 r[2] eq i]]> : i in P[1][1]];
 x2_edges := [<i,[t[2] : t in T[2][2] | t[1] in [r[1] : r in P[2][2] |
              r[2] eq i]]> : i in P[1][2]];
 y_vertices := [<i,[t[2] : t in T[3][1] | t[1] in [r[1] : r in P[3][1] |
                r[2] eq i]]> : i in P[1][1]];
 y_edges := [<i,[t[2] : t in T[3][2] | t[1] in [r[1] : r in P[3][2] |
             r[2] eq i]]> : i in P[1][2]];
 c_vertices := [<r[1],r[2][1]> : r in x2_vertices | not #r[2] eq 0]
                cat [<r[1],r[2][1]> : r in y_vertices | not #r[2] eq 0];
 c_edges := [<r[1],r[2][1]> : r in x2_edges | not #r[2] eq 0]
             cat [<r[1],r[2][1]> : r in y_edges | not #r[2] eq 0];
 c := <RedSeq(c_vertices),RedSeq(c_edges)>; // %%
 return c;
end function;

IsPushoutGraphs := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 p := InducedMorphismGraphsPO(X,X2,Y,Y2,u,u2,v,v2); // %%
 return Is_Bijective(p,PushoutGraphs(X,Y,X2,v,u)[1],Y2); // %%
end function;

PullbackSets_num := function(X,Y,Y2,f,g) // f: X -> Y, g: Y2 -> Y
                                         // maps between sets
  // given as a list of tuples, e.g. f = [<x_1,y_1>,...,<x_n,y_n>]
 P := Sort([<x,y2> : x in X, y2 in Y2 | #[<s,t> : s in f, t in g |
            x eq s[1] and y2 eq t[1] and s[2] eq t[2]] eq 1]);
 PP := [i : i in [1..#P]];
 v := [<i,P[i][1]> : i in PP]; // <p,p@v>
 u := [<i,P[i][2]> : i in PP]; // <p,p@u>
 num := [<i,P[i]> : i in [1..#P]];
 return <PP,u,v,num>; // pullback PP, u: PP -> Y2, v: PP -> X
end function;

PullbackGraphs_num := function(X,Y,Y2,f,g) // f: X -> Y, g: Y2 -> Y
                                           // graph morphisms
 vertices := PullbackSets(X[1],Y[1],Y2[1],f[1],g[1]); // %%
 edges := PullbackSets(X[2],Y[2],Y2[2],f[2],g[2]); // %%
 N := [i : i in [1..#vertices]];
 E := [i : i in [1..#edges]]; // edges without source and target
 EE := [<Index(vertices,<edges[e][1][1],edges[e][2][1]>),e,
   Index(vertices,<edges[e][1][3],edges[e][2][3]>)> : e in E];
   // edges with source and target
 PP := <N, EE>;
 vE := [<ee,edges[ee[2]][1]> : ee in EE]; // <p,p@v>
 uE := [<ee,edges[ee[2]][2]> : ee in EE]; // <p,p@u>
 vN := [<n,vertices[n][1]> : n in N]; // <p,p@v>
 uN := [<n,vertices[n][2]> : n in N]; // <p,p@u>
 v := <vN,vE>;
 u := <uN,uE>;
 numvertices := [<i,vertices[i]> : i in [1..#vertices]];
 numedges := [<i,edges[i]> : i in [1..#edges]];
 num := <numvertices,numedges>;
 return <PP,u,v,num>; // pullback PP, u: PP -> Y2, v: PP -> X
end function;

InducedMorphismSetsPB := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 T := <X,v,u>;
 P := PullbackSets_num(Y,Y2,X2,v2,u2); // %%
 c_uncode := [<t,<[r[2] : r in T[2] | r[1] eq t][1],[r[2] : r in T[3] |
              r[1] eq t][1]>> : t in T[1]];
 c := [<r[1],Index([p[2] : p in P[4]],r[2])> : r in c_uncode];
 return c;
end function;

InducedMorphismGraphsPB := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 T := <X,v,u>;
 P := PullbackGraphs_num(Y,Y2,X2,v2,u2); // %%
 c_uncode_vertices := [<t,<[r[2] : r in T[2][1] | r[1] eq t][1],[r[2] :
                       r in T[3][1] | r[1] eq t][1]>> : t in T[1][1]];
 c_uncode_edges := [<t,<[r[2] : r in T[2][2] | r[1] eq t][1],[r[2] :
                    r in T[3][2] | r[1] eq t][1]>> : t in T[1][2]];
 c_vertices := [<r[1],Index([p[2] : p in P[4][1]],r[2])> :
                r in c_uncode_vertices];
 c_edges := [<r[1],Index([p[2] : p in P[4][2]],r[2])> : r in c_uncode_edges];
 c := <c_vertices,c_edges>;
 return c;
end function;

IsPullbackGraphs := function(X,X2,Y,Y2,u,u2,v,v2)
 // u: X -> X2, u2: X2 -> Y2, v: X -> Y, v2: Y -> Y2
 p := InducedMorphismGraphsPB(X,X2,Y,Y2,u,u2,v,v2); // %%
 return Is_Bijective(p,X,PullbackGraphs_num(Y,Y2,X2,v2,u2)[1]); // %%
end function;

CyclesFromVertex := function(x,G)  // G graph, x vertex in G
 if #[a : a in G[2] | a[1] eq x] eq 0 then
  return false;
 end if;
 S := {};   // S: set of vertices to achieve
 Snew := {x};
 while not #Snew eq #S do
  S := Snew;
  Snew join:= {a[3] : a in G[2] | a[1] in S};
  if x in {a[3] : a in G[2] | a[1] in S} then
   return true;
  end if;
 end while;
 return false;
end function;

VerticesToAchieve := function(G,r) // G: graph, r vertex in G
 if #[a : a in G[2] | a[1] eq r] eq 0 then
  return [];
 end if;
 S := {};   // S: set of vertices to achieve
 Snew := {r};
 while not #Snew eq #S do
  S := Snew;
  Snew join:= {a[3] : a in G[2] | a[1] in S};
 end while;
 return Sort(SetToSequence(Snew));
end function;

CyclesInPathFromx := function(x,G)
 return &or[CyclesFromVertex(v,G) : v in VerticesToAchieve(G,x)];
end function;

Paths := function(x,G) // G graph, x vertex in G[1],
 if not CyclesInPathFromx(x,G) then // %%
  Listofpaths := [[[]]];
  if not #[[[e] : e in G[2] | e[1] eq x]] eq 0 then
   Listofpaths cat:= [[[e] : e in G[2] | e[1] eq x]];
  else
   return &cat(Listofpaths);
  end if;
  newpaths := [<p,[e : e in G[2] | e[1] eq p[#p][3]]> :
               p in Listofpaths[#Listofpaths] | not #[e : e in G[2] |
               e[1] eq p[#p][3]] eq 0];
  while not #newpaths eq 0 do
   Listofpaths cat:= [&cat[&cat[[n[1] cat [n[2][i]]] : i in [1..#n[2]]] :
                      n in newpaths]];
   newpaths := [<p,[e : e in G[2] | e[1] eq p[#p][3]]> :
               p in Listofpaths[#Listofpaths] | not #[e : e in G[2] |
               e[1] eq p[#p][3]] eq 0];
  end while;
  return Listofpaths;
 end if;
 return "infinite";
end function;

TreeOfPaths := function(x,G) // G graph, x vertex in G
 if not CyclesInPathFromx(x,G) then // %%
  P := Paths(x,G); // %%
  V := &cat(P);
  E := [<[],<[],p[1],p>,p> : p in P[2]];
  E cat:= [<[p[i] : i in [1..#p-1]],<[p[i] : i in [1..#p-1]],p[#p],p>,p> :
           p in V | #p ge 2];
  return <V,E>;
 end if;
 return "infinite";
end function;

IsTreeDef := function(G)
 for r in G[1] do // searching root
  if #[v : v in [g : g in G[1] | not g eq r] | not #[e : e in G[2] |
       e[3] eq v] eq 1] eq 0 then // (Tree 1)
   if #[e : e in G[2] | e[3] eq r] eq 0 then // (Tree 2)
    if &and[v in VerticesToAchieve(G,r) : v in G[1]] then // (Tree 3) // %%
     return true;
    end if;
   end if;
  end if;
 end for;
 return false;
end function;

ListOfnCycles := function(G,n) // G: graph
 return ListGraphMorphisms(C(n),G); // %%
end function;

Cnf_Bij := function(f,G,H,n) // G, H: graphs, f: G -> H: graph morphism
 if not #ListOfnCycles(G,n) eq #ListOfnCycles(H,n) then // %%
  return false;
 end if;
 if SequenceToSet([ComposeGraphMorphisms(ListOfnCycles(G,n)[i],f) : i in
    [1..#ListOfnCycles(G,n)]]) eq SequenceToSet(ListOfnCycles(H,n)) then
  return true; // %%
 end if;
 return false;
end function;

IsQis_Bound := function(f,G,H,ub) // ub: upper bound
 i := 1;
 while Cnf_Bij(f,G,H,i) and i le ub do // %%
  i := i+1;
 end while;
 return i eq ub+1;
end function;

IsFibration := function(f,G,H) // G, H: graphs, f: G -> H graph morphism
 for x in G[1] do
  y := [a[2] : a in f[1] | a[1] eq x][1];
  for b in [h : h in H[2] | h[1] eq y] do
   if #[0 : a in G[2] | <a,b> in f[2] and a[1] eq x] eq 0 then
    return false;
   end if;
  end for;
 end for;
 return true;
end function;

IsEtaleFibration := function(f,G,H) // G, H: graphs, f: G -> H graph morphism
 for x in G[1] do
  y := [a[2] : a in f[1] | a[1] eq x][1];
  if not #[e : e in G[2] | e[1] eq x] eq #[h : h in H[2] | h[1] eq y] then
   return false;
  elif not Sort([a[2] : a in f[2] | a[1][1] eq x])
           eq Sort([h : h in H[2] | h[1] eq y]) then
    return false;
  end if;
 end for;
 return true;
end function;

IsFibrant := function(X)
 return &and[not #[e : e in X[2] | e[1] eq v] eq 0: v in X[1]];
end function;

AcCofib1to4 := function(f,G,H) // f: G -> H graph morphism
 if #[0 : x in H[1] | #[0 : a in G[1] | <a,x> in f[1]] ge 2] ge 1 then
    // (AcCofib 1)
  return false;
 end if;
 if #[0 : x in H[2] | #[0 : a in G[2] | <a,x> in f[2]] ge 2] ge 1 then
    // (AcCofib 2)
  return false;
 end if;
 HH1 := [x : x in H[1] | #[0 : a in G[1] | <a,x> in f[1]] eq 0];
 HH2 := [x : x in H[2] | #[0 : a in G[2] | <a,x> in f[2]] eq 0];
 if #[0 : x in HH1 | not #[a : a in H[2] | a[3] eq x] eq 1] ge 1 then
    // (AcCofib 3)
  return false;
 end if;
 if #[0 : x in HH2 | not x[3] in HH1] ge 1 then
    // (AcCofib 4)
  return false;
 end if;
 return true;
end function;

AcCofib5 := function(f,G,H)
 HH1 := [x : x in H[1] | #[0 : a in G[1] | <a,x> in f[1]] eq 0];
 max := #H[2];
 L := [];
 for i in [1..max] do
  L cat:= [ListGraphMorphisms(D(i),H)]; // %%
 end for;
 for v in HH1 do
  list := [];
  for i in [1..max] do
   list cat:= [l : l in L[i] | not l[1][1][2] in HH1 and l[1][#l[1]][2] eq v];
  end for;
  if #list eq 0 then
   return false;
  end if;
 end for;
 return true;
end function;

IsAcCofib := function(f,G,H)
 return AcCofib1to4(f,G,H) and AcCofib5(f,G,H); // %%
end function;

IsTree := function(G)
 for x in G[1] do
  if IsAcCofib(VtoE(D(0),G,[<0,x>]),D(0),G) then // %%
   return true;
  end if;
 end for;
 return false;
end function;

Unitargeting := function(f,G,H)
 return [e : e in H[2] | #RedSeq([ee[1][3] : ee in f[2] | ee[2] eq e]) eq 1];
  // %%
end function;

Uni := function(f,G,H)
 U := Unitargeting(f,G,H); // %%
 HH := <H[1],[e : e in H[2] | not e in U]>; // H without unitargeting edges
 n := Minimum([#HH[1],#HH[2]]);
 return &and[#ListGraphMorphisms(C(i),HH) eq 0 : i in [1..n]]; // %%
end function;

SuffCond := function(f,G,H)
 return <IsEtaleFibration(f,G,H),Uni(f,G,H)>; // %%
end function;

c2chain := function(n)
 edges := [<i,i,i+1> : i in [1..n-1]] cat [<i+1,2*n-1-i,i> : i in [1..n-1]];
 edges_tosort := Sort([<e[2],e[1],e[3]> : e in edges]);
 edges := [<e[2],e[1],e[3]> : e in edges_tosort];
 return <[i : i in [1..n]],edges>;
end function;

trygraph := function(n) // n geq 3
 edges := [<3*n-5,3*n-5,3*n-4>,<3*n-4,3*n-4,3*n-5>,<3*n-5,6*n-11,3*n-8>,
           <3*n-6,6*n-12,3*n-5>];
 edges cat:= &cat[[<3*k-2,3*k-2,3*k-1>,<3*k-1,3*k-1,3*k-2>,<3*k-1,3*k,3*k>] :
                 k in [1..n-2]]; // innerhalb der Stufen
 edges cat:= &cat[[<3*k+1,3*n-5+3*k,3*k-2>,<3*k,3*n-6+3*k,3*k+1>,
                   <3*k,3*n-4+3*k,3*k+3>] : k in [1..n-3]];
 edges_tosort := Sort([<e[2],e[1],e[3]> : e in edges]);
 edges := [<e[2],e[1],e[3]> : e in edges_tosort];
 G := <[i : i in [1..3*n-4]],edges>;
 return G;
end function;

tryacyclic := function(n) // trygraph -> c2chain
 vertices := [<1,1>,<2,2>,<4,2>,<3*n-6,n>,<3*n-4,n>];
 vertices cat:= &cat[[<3*k-6,k>,<3*k-4,k>,<3*k-2,k>] : k in [3..n-1]];
 if n eq 3 then
  edges := [<<1,1,2>,<1,1,2>>,<<2,2,1>,<2,2*n-2,1>>,<<2,3,3>,<2,2,3>>,
            <<4,4,5>,<2,2,3>>,<<5,5,4>,<3,2*n-3,2>>,
            <<3,3*n-3,4>,<3,2*n-3,2>>,<<4,3*n-2,1>,<2,2*n-2,1>>];
 else
  edges := [<<1,1,2>,<1,1,2>>,<<2,2,1>,<2,2*n-2,1>>,<<2,3,3>,<2,2,3>>,
            <<4,4,5>,<2,2,3>>,<<5,5,4>,<3,2*n-3,2>>,
            <<3,3*n-3,4>,<3,2*n-3,2>>,<<4,3*n-2,1>,<2,2*n-2,1>>,
            <<7,3*n+1,4>,<3,2*n-3,2>>,<<3*n-7,3*n-6,3*n-6>,<n-1,n-1,n>>,
            <<3*n-5,3*n-5,3*n-4>,<n-1,n-1,n>>,
            <<3*n-4,3*n-4,3*n-5>,<n,n,n-1>>,
            <<3*n-9,6*n-13,3*n-6>,<n-1,n-1,n>>,
            <<3*n-6,6*n-12,3*n-5>,<n,n,n-1>>];
  edges cat:= &cat[[<<3*k-4,3*k-3,3*k-3>,<k,k,k+1>>,
   <<3*k-2,3*k-2,3*k-1>,<k,k,k+1>>,<<3*k-1,3*k-1,3*k-2>,<k+1,2*n-1-k,k>>,
   <<3*k-6,3*n-10+3*k,3*k-3>,<k,k,k+1>>,
   <<3*k-3,3*n-9+3*k,3*k-2>,<k+1,2*n-1-k,k>>,
   <<3*k+1,3*n-5+3*k,3*k-2>,<k+1,2*n-1-k,k>>] : k in [3..n-2]];
 end if;
 edges_tosort := Sort([<<e[1][2],e[1][1],e[1][3]>,e[2]> : e in edges]);
 edges := [<<e[1][2],e[1][1],e[1][3]>,e[2]> : e in edges_tosort];
 return <Sort(vertices),edges>;
end function;

glue_vertices := function(G,list) // list := [[2,3],[1,4,5]] list of
 // sublists of vertices to glue
 list := [RedSeq(l) : l in list | #RedSeq(l) ge 2]; // %%
 vertices_to_glue := RedSeq(&cat(list)); // %%
 vertices_left_over := [n : n in G[1] | not n in vertices_to_glue];
 vertices := Sort([l[1] : l in list] cat vertices_left_over);
 edges_1 := [e : e in G[2] | not e[1] in vertices_to_glue
             and not e[3] in vertices_to_glue];
 edges := [];
 for e in G[2] do
  if not e[1] in vertices_to_glue and not e[3] in vertices_to_glue then
   edges cat:= [e];
  else if not e[1] in vertices_to_glue then
   edges cat:= [<e[1],e[2],[l[1] : l in list | e[3] in l][1]>];
  else
   edges cat:= [<[l[1] : l in list | e[1] in l][1],e[2],e[3]>];
  end if;
 end if;
 end for;
 edges_named := [];
 for e in edges do
  if e[1] in &cat(list) then
   e1_new := [l : l in list | e[1] in l][1][1];
  else
   e1_new := e[1];
  end if;
  if e[3] in &cat(list) then
   e3_new := [l : l in list | e[1] in l][1][1];
  else
   e3_new := e[3];
  end if;
  edges_named cat:= [<e1_new,e[2],e3_new>];
 end for;
 return <vertices,edges_named>;
end function;

glue_vertices_including_edges := function(G,list)
 G := glue_vertices(G,list); // %%
 edges := [];
 for g in G[2] do
  if #[e : e in edges | e[1] eq g[1] and e[3] eq g[3]] eq 0 then
   edges cat:= [g];
  end if;
 end for;
 return <G[1],edges>;
end function;

try_id_vertices := function(n) // n ge 3
 return [[3*i,3*i+2] : i in [1..n-2]];
end function;

idtrygraph := function(n)
 return glue_vertices_including_edges(trygraph(n),try_id_vertices(n));
  // %%
end function;

tryfactorization := function(n) // trygraph --> idtrygraph
 T := try_id_vertices(n); // %%
 G := trygraph(n); // %%
 vertices := Sort([<t[2],t[1]> : t in T] cat [<i,i> : i in [1..3*n-4] |
                   not i in [t[2] : t in T]]);
 edges_to_map := trygraph(n)[2];
 edges_images := idtrygraph(n)[2];
 edges := [];
 for e in edges_to_map do
  if e in edges_images then
   im_e := e;
  else
   im_e1 := [n[2] : n in vertices | n[1] eq e[1]][1];
   im_e3 := [n[2] : n in vertices | n[1] eq e[3]][1];
   im_e := [edge : edge in edges_images |
            edge[1] eq im_e1 and edge[3] eq im_e3][1];
  end if;
  edges cat:= [<e,im_e>];
 end for;
 return <vertices,edges>;
end function;

idtryacyclic := function(n) // idtrygraph -> c2chain
 T := try_id_vertices(n); // %%
 T2 := [t[2] : t in try_id_vertices(n)]; // %%
 G :=  idtrygraph(n); // %%
 f := tryacyclic(n); // %%
 vertices := [n : n in f[1] | n[1] in G[1]];
 edges := [];
 for e in f[2] do
  if e[1] in G[2] then
   edges cat:= [e];
  else
   e1 := e[1][1];
   e3 := e[1][3];
   if e[1][1] in T2 then
    e1 := [t[1] : t in T | t[2] eq e[1][1]][1];
   end if;
   if e[1][3] in T2 then
    e3 := [t[1] : t in T | t[2] eq e[1][3]][1];
   end if;
   edges cat:= [<<e1,e[1][2],e3>,e[2]>];
  end if;
 end for;
 edges2 := [];
 for e in edges do
  if not <e[1][1],e[1][3]> in [<e[1][1],e[1][3]> : e in edges2] then
   edges2 cat:= [e];
  end if;
 end for;
 return <vertices,edges2>;
end function;

Doublecyclic := function(n)
 C := c2chain(n); // %%
 return <C[1],C[2] cat [<1,#C[2]+1,#C[1]>] cat [<#C[1],#C[2]+2,1>]>;
end function;

Trygraph := function(n)
 T := trygraph(n); // %%
 return <T[1],T[2] cat [<1,#T[2]+1,#T[1]>, <#T[1],#T[2]+2,1>, <#T[1]-2,#T[2]+3,1>]>;
end function;

Tryacyclic := function(n) // Trygraph -> Doublecyclic
 T := Trygraph(n); // %%
 t := #T[2];
 D := Doublecyclic(n); // %%
 d := #D[2];
 f := tryacyclic(n); // %%
 return <f[1],f[2] cat [<T[2][t-2],D[2][d-1]>, <T[2][t-1],D[2][d]>, <T[2][t],D[2][d]>]>;
end function;

idTrygraph := function(n)
 T := idtrygraph(n); // %%
 t := [r[2] : r in T[2]];
 return <T[1],T[2] cat [<1,t[#t]+1,T[1][#T[1]-1]>]
         cat [<T[1][#T[1]-1],t[#t]+2,1>]>;
end function;

Tryfactorization := function(n) // Trygraph -> idTrygraph
 T := trygraph(n); // %%
 TT := idtrygraph(n); // %%
 t := [r[2] : r in TT[2]];
 f := tryfactorization(n); // %%
 return <f[1],f[2] cat [<<1,#T[2]+1,#T[1]>,<1,t[#t]+1,TT[1][#TT[1]-1]>>]
         cat [<<#T[1],#T[2]+2,1>,<TT[1][#TT[1]-1],t[#t]+2,1>>]
         cat [<<#T[1]-2,#T[2]+3,1>,<TT[1][#TT[1]-1],t[#t]+2,1>>]>;
end function;

idTryacyclic := function(n) // idTrygraph -> Doublecyclic
 f := idtryacyclic(n); // %%
 T := idtrygraph(n); // %%
 t := [r[2] : r in T[2]];
 C := c2chain(n); // %%
 return <f[1],f[2] cat [<<1,t[#t]+1,T[1][#T[1]-1]>,<1,#C[2]+1,#C[1]>>]
        cat [<<T[1][#T[1]-1],t[#t]+2,1>,<#C[1],#C[2]+2,1>>]>;
end function;

CnCm := function(n,m)
 V := [i : i in [1..n+m]];
 E := [<i,i,i+1> : i in [1..n+m-1]] cat [<n+m,n+m,1>]
      cat [<n,n+m+1,1>,<n+m,n+m+2,n+1>];
 return <V,E>;
end function;

cncm := function(n,m) // cn glued to cm at vertex n
 V := [i : i in [1..n+m-1]];
 E := [<i,i,i+1> : i in [1..n-1]] cat [<n,n,1>] cat [<n+i,n+i+1,n+i+1> :
       i in [0..m-2]] cat [<n+m-1,n+m,n>];
 return <V,E>;
end function;

cncmqis := function(n,m) // CnCm -> cncm
 return VtoE(CnCm(n,m),cncm(n,m),[<i,i> : i in [1..n+m-1]] cat [<n+m,n>]);
  // %%
end function;

VtoE := function(G,H,Vf) // H thin, Vf = [<1,2>,<2,5>,<3,1>,<4,2>] map on
                         // vertices
 Ef := [<e,[h : h in H[2] | h[1] eq [v[2] : v in Vf | v[1] eq e[1]][1] and
         h[3] eq [v[2] : v in Vf | v[1] eq e[3]][1]][1]> : e in G[2]];
 return <Vf,Ef>;
end function;

CNCN := function(n)
 G := DUC([n,n]); // %%
 H := <G[1],G[2] cat [<n+1,2*n+1,1>,<2*n,2*n+2,n>]>;
 f := VtoE(G,H,[<i,i> : i in [1..2*n]]); // %%
 return <G,H,f>;
end function;

c2graph := function(list) // list e.g. [<1,3>,<2,4>,<6,6>] list of tuples
                          // of vertices
 n := Maximum([x[1] : x in list] cat [x[2] : x in list]);
 edges := [];
 for x in list do
  i := Index(list,x);
  edges cat:= [<x[1],i,x[2]>,<x[2],2*#list-i+1,x[1]>];
 end for;
 edges_tosort := Sort([<e[2],e[1],e[3]> : e in edges]);
 edges := [<e[2],e[1],e[3]> : e in edges_tosort];
 return <[i : i in [1..n]],edges>;
end function;

Exflower := function(n,list)  // list := [2,3,5] contains vertices
 // that are connected with "upper" vertex 1 in G
 edges_G := Sort(&cat[[<1,i>,<n+2,i>] : i in [2..n+1]]);
 if #list eq 0 then // alle unten
  edges_G cat:= [<i,n+2> : i in [2..n+1]];
  edges_G := [<edges_G[i][1],i,edges_G[i][2]> : i in [1..#edges_G]];
 else
  edges_G cat:= [<i,1> : i in list];
  edges_G cat:= [<i,n+2> : i in [2..n+1] | not i in list];
  edges_G := [<edges_G[i][1],i,edges_G[i][2]> : i in [1..#edges_G]];
 end if;
 flowerG := <[i : i in [1..n+2]],edges_G>;
 flowerH := c2graph([<1,i> : i in [2..n+1]]); // %%
 flowerf := VtoE(flowerG,flowerH,[<i,i> : i in [1..n+1]] cat [<n+2,1>]);
  // %%
 return <flowerG,flowerH,flowerf>;
end function;

Exflower2 := function(n,k)
 // n: number of "petals", k: size of "petals"
 G := DUC([k : i in [1..n]]); // %%
 VG := G[1];
 EG := G[2];
 V := [1+i*k : i in [0..n-1]];
 EG2 := &cat[[<v,vv> : vv in [vv+1 : vv in [vv : vv in V | not vv eq v]]] :
        v in V];
 EG2 := [<e[1],#EG+Index(EG2,e),e[2]> : e in EG2];
 flowerG := <VG,EG cat EG2>;
 VH := [i : i in [1..n*(k-1)+1]];
 EH := [<1,1+i*k,2+i*(k-1)> : i in [0..n-1]];
 if not k eq 2 then
  for i in [1..n] do
   EH cat:= [<1+(i-1)*(k-1)+l,1+(i-1)*k+l,2+(i-1)*(k-1)+l> : l in [1..k-2]];
  end for;
 end if;
 EH cat:= [<k+i*(k-1),k*(i+1),1> : i in [0..n-1]];
 EH_sort := [<e[2],e[1],e[3]> : e in Sort([<e[2],e[1],e[3]> : e in EH])];
 flowerH := <VH,EH_sort>; // %%
 VtoE_vertices := [<v,1> : v in V] cat &cat[[<1+(i-1)*k+l,2+(i-1)*k+l-i> :
                   l in [1..k-1]] : i in [1..n]];
 flowerf := VtoE(flowerG,flowerH,VtoE_vertices); // %%
 return <flowerG,flowerH,flowerf>;
end function;

exampleforbadbound := function(n) // (c_2n,f) first not to be bijective
 V := [i : i in [1..2*n+1]];
 E := [<i,i+1,i+1> : i in [3..2*n]];
 E cat:= [<1,1,2>, <2,2,1>, <1,3,3>, <2*n+1,2*n+2,1>];
 X := <V,E>;
 f := ListGraphMorphisms(X,C(2))[2]; // %%
 return <X,f>;
end function;

exampleforbadbound2 := function(n) // (c_n,f) first not to be bijective
 G := DUC([1,n]); // %%
 f := <[<v,1> : v in G[1]],[<e,<1,1,1>> : e in G[2]]>;
 return <G,C(1),f>;
end function;

exampleforbadbound3 := function(n) // n even, n ge 4
 G := DUC([2,n-2,n]); // %%
 H := <[i : i in [1..n-1]],[<i,i,i+1> : i in [1..n-2]]
       cat [<n-1,n-1,n-2>,<n-2,n,1>]>;
 u_vertices := [<1,n-1>,<2,n-2>] cat [<i,i-2> : i in [3..n]]
  cat [<i,i-n> : i in [n+1..2*n-1]] cat [<2*n,n-2>];
 u_edges := [<c,[h : h in H[2] | h[1] eq [u[2] : u in u_vertices |
   u[1] eq c[1]][1] and h[3] eq [u[2] : u in u_vertices |
   u[1] eq c[3]][1]][1]> : c in G[2]];
 u := <u_vertices,u_edges>;
 return <G,H,u>;
end function;

DTB := function(n) // decimal to binary
 A := [];
 if n eq 0 then
  return [0];
 end if;
 while n gt 0 do
  A cat:= [n mod 2];
  n := Integers()!((n-(n mod 2))/2);
 end while;
 return A;
end function;

DTO := function(list) // decimal to other, e.g. list := [2,4,3]
 W := [[i : i in [1..list[1]]]];
 for i in [2..#list] do
  W cat:= [[u : u in [W[#W][#W[#W]]+1..W[#W][#W[#W]]+list[i]]]];
 end for;
 n := #list;
 A := [];
 for i in [0..&*list-1] do
  j := i;
  AA := [];
  for l in [list[#list-i] : i in [0..#list-1]] do
   AA cat:= [j mod l];
   j := Integers()!((j-(j mod l))/l);
  end for;
  A cat:= [[AA[#AA-i] : i in [0..#AA-1]]];
 end for;
 return A;
end function;

Thins := function(n)
 PV := [<i,j> : i,j in [1..n]]; // pairs of vertices for one edge
 LG := []; // list of thin graphs
 for i in [0..2^(n^2)-1] do
  A := DTB(i); // %%
  E := [PV[i] : i in [1..#A] | A[i] eq 1];
  E := [<e[1],Index(E,e),e[2]> : e in E];
  LG cat:= [<[i : i in [1..n]],E>];
 end for;
 return LG;
end function;

Thins2 := function(n)
 PV := [<i,j> : i,j in [1..n]]; // pairs of vertices for one edge
 LG := []; // list of thin graphs
 for i in [0..2^(n^2)-1] do
  A := DTB(i); // %%
  E := [PV[i] : i in [1..#A] | A[i] eq 1];
  E := [e : e in E | not e[1] eq e[2]];
  E := [<e[1],Index(E,e),e[2]> : e in E];
  LG cat:= [<[i : i in [1..n]],E>];
 end for;
 return RedSeq(LG); // %%
end function;

EFU := function(H,F) // returns all etale fibrations f: G -> H that
 // satisfy (Uni)
 // F: sizes of fibres in list, e.g. F := [1,2,1]
 n := #H[1];
 VG := [i : i in [1..&+F]];
 W := [[i : i in [1..F[1]]]];
 for i in [2..#F] do
  W cat:= [[u : u in [W[#W][#W[#W]]+1..W[#W][#W[#W]]+F[i]]]];
 end for;
 Vf := [<i,Index(W,[a : a in W | i in a][1])> : i in [1..&+F]];
 EG_all := [[<i,j> : i in W[e[1]], j in W[e[3]]] : e in H[2]];
 Ef_all_sort := &cat([[<[e : e in EG_all[i] | e[1] eq v],H[2][i]> :
  i in [1..#EG_all]]: v in VG]);
 Ef_all_sort := [e : e in Ef_all_sort | not #e[1] eq 0];
 Ef_poss := [[<Ef_all_sort[i][1][d[i]+1],Ef_all_sort[i][2]> : i in [1..#d]] :
  d in DTO([#e[1] : e in Ef_all_sort])]; // %%
 Ef_poss_numbered := [[<<Ef_poss[i][j][1][1],j,Ef_poss[i][j][1][2]>,
  Ef_poss[i][j][2]> : j in [1..#Ef_poss[i]]] : i in [1..#Ef_poss]];
 EF := [<<Vf,Ef_poss_numbered[i]>,<VG,[Ef_poss_numbered[i][j][1] : j in
  [1..#Ef_poss_numbered[i]]]>,H> : i in [1..#Ef_poss_numbered]];
       // <f,G,H> all possibilities that f is an etale fibration
 return [f : f in EF | &and[SuffCond(f[1],f[2],f[3])[i] : i in [1,2]]];
end function;

IsEqual := function(f,g) // returns true if graph morphisms
  // f, g: G -> H are equal
  // returns true if graphs f, g are equal
 return (SequenceToSet(f[1]) eq SequenceToSet(g[1]))
        and (SequenceToSet(f[2]) eq SequenceToSet(g[2]));
end function;

IsoRepresentatives := function(T,n); // returns representatives of isoclasses
 TT := [T[1]];
 for i in [1..n] do
  if not &or[IsIsomorphic(T[i],t)[2] : t in TT] then // %%
   TT cat:= [T[i]];
  end if;
 end for;
 return TT;
end function;

Is_Comm_Quad_Graphs := function(X,Y,X2,Y2,f,f2,g,h)
  // f: X -> Y, f2: X2 -> Y2, g: X -> X2, h: Y -> Y2
 if IsEqual(ComposeGraphMorphisms(f,h),ComposeGraphMorphisms(g,f2)) then // %%
  return true;
 end if;
 return false;
end function;

Lift := function(X,Y,X2,Y2,f,f2,g,h) // f: X -> Y, f2: X2 -> Y2,
  // g: X -> X2, h: Y -> Y2
 L := ListGraphMorphisms(X2,Y); // %%
 L_comm := [l : l in L | IsEqual(ComposeGraphMorphisms(g,l),f) and
            IsEqual(ComposeGraphMorphisms(l,h),f2)]; // %%
 if not #L_comm eq 0 then
  return L_comm[1];
 end if;
 return false;
end function;

remove_edges := function(G,list)
 return <G[1],[e : e in G[2] | not e in list]>;
end function;

remove_vertices := function(G,list)
 return <[n : n in G[1] | not n in list],[e : e in G[2] |
         not e[1] in list and not e[2] in list]>;
end function;

graph_op := function(G); // functor, reverses the direction of arrows in graph
  // exchanges source and target of edges in graph
 return <G[1],[<e[3],e[2],e[1]> : e in G[2]]>;
end function;

connection := function(G) // returns the list of connected components
 list := [];
 G_done := {};
 while not #G_done eq #G[1] do
 S := {};   // S: set of vertices to achieve
  r := [g : g in G[1] | not g in G_done][1];
  Snew := {r};
  while not #Snew eq #S do
   S := Snew;
   Snew join:= {a[1] : a in G[2] | a[3] in S} join {a[3] : a in G[2] | a[1] in S};
  end while;
  G_done join:= Snew;
  Snew := Sort(SetToSequence(Snew));
  edges := Sort([a : a in G[2] | a[1] in Snew or a[3] in Snew]);
  list cat:= [<Snew,edges>];
 end while;
 return list;
end function;