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


IsRightUnique := function(u) // JH
 // 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) // JH
 // 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); // JH
 // 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; // JH
 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;

ListGraphEpimorphisms := function(G,H) // G und H sollten keine isolierten Vertices haben !!! Sonst wird Ergebnis falsch
    if #G[1] lt #H[1] or #G[2] lt #H[2] then
        return [];
    end if;
    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_epi := [x : x in list | #{t[2] : t in x} eq #H[2]];
    list_epi_mor := [];
    for y in list_epi do
        list_completions_to_maps_on_vertices := [Sort(SetToSequence(x)) : x in
                                                 CompletionsToMaps(G[1],H[1],RVPME(y))]; // %%
        list_epi_mor cat:= [<x,y> : x in list_completions_to_maps_on_vertices];
    end for;
    return list_epi_mor;
end function;


Is_Injective := function(m,G,H) // JH
// 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) // JH
// 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) // JH
    return Is_Injective(f,G,H) and Is_Surjective(f,G,H); // %%
end function;


// With the function IsIsomorphic we can test if given graphs $G$ and $H$ are isomorphic.
//In case of isomorphic graphs, we additionally give a graph isomorphism as output.
IsIsomorphic := function(G,H) // G, H: graphs // JH
    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;






// Input:
//   G - a graph represented as a pair <vertices, edges>
//   v - a vertex of G
//
// Output:
//   A list of all edges in G with source vertex v
GetEdgesFromVertex := function(v,G)
    return [e : e in G[2] | e[1] eq v];
end function;

// Input:
//   G - a graph represented as a pair <vertices, edges>
//   v - a vertex of G
//
// Output:
//   A list of all edges in G with target vertex v
GetEdgesTargetVertex := function(v,G)
    return [e : e in G[2] | e[3] eq v];
end function;
		

// Input:
//   G - a graph represented as a pair <vertices, edges>
//   v - a vertex of G
//
// Output:
//   A list of all vertices that are adjacent to v.
GetAdjacentVertices := function (v,G)
    adjacent := [];
    for e in G[2] do
        if e[3] eq v then
            adjacent := adjacent cat [e[1]];
        end if;
        if e[1] eq v then
            adjacent := adjacent cat [e[3]];
        end if;
    end for;
    return SetToSequence (Set( adjacent ));
end function;


// Input: 
//   G - a graph 
//   v - a vertex of G
//
// Output: 
//   Set of all vertices in G that are reachable
//   from v via a directed path
GetAllReachableVertices := function(v, G)
    reachable := {v};  
    frontier := {v};  
    while not IsEmpty(frontier) do
        nextFrontier := {};
        for u in frontier do
            for e in GetEdgesFromVertex(u, G) do
                w := e[3];  
                if w notin reachable then
                    Include(~reachable, w);
                    Include(~nextFrontier, w);
                end if;
            end for;
        end for;
        frontier := nextFrontier;
    end while;
    return reachable;
end function;


// Input: 
//   G - a graph
//
// Output: 
//   true if G is weakly connected; 
//   false otherwise
IsWeaklyConnected := function(G)
    if #G[1] le 1 then
        return true;
    else
        reachable := {G[1][1]};  
        frontier := reachable; 
        while not IsEmpty(frontier) do
            nextFrontier := {};
            for u in frontier do
                nextFrontier := nextFrontier join (Set(GetAdjacentVertices(u,G)) diff reachable);
            end for;
            reachable := reachable join nextFrontier;
            frontier := nextFrontier;
        end while;
    end if;
    if reachable eq Set(G[1]) then
        return true;
    else
        return false;
    end if;
end function;


// Input: 
//   G - a graph
//
// Output: 
//   true if G is star-connected; 
//   false otherwise
IsStarConnected := function(G)
    for vertex in G[1] do
        if Set(GetAllReachableVertices(vertex,G)) eq Set(G[1]) then
            return true;
        end if;
    end for;
    return false;
end function;


// Input: 
//   G - a graph
//
// Output: 
//   true if G is strongly connected; 
//   false otherwise
IsStronglyConnected := function(G)
    if #G[1] eq 0 then
        return false;
    end if;
    for v in G[1] do
        if Set(G[1]) ne GetAllReachableVertices(v,G) then
            return false;
        end if;
    end for;
    return true;
end function;


// Input: 
//   G1,G2 - graphs 
//
// Output: 
//   true if G1 and G2 are equal; 
//   false otherwise
AreGraphsEqual := function(G1,G2)
return Set(G1[1]) eq Set(G2[1]) and Set(G1[2]) eq Set(G2[2]);
end function;


// Input: 
//   G    - a graph
//   set  - a set of graphs
//
// Output: 
//   true if G is already in the set;
//   false otherwise
IsGraphInSet := function(G,set)
    for graph in set do
        if AreGraphsEqual(G,graph) eq true then
            return true;
        end if;
    end for;
    return false;
end function;


// Input: 
//   G    - a graph
//   list  - a list of graphs
//
// Output: 
//   true if G is already in the list;
//   false otherwise
IsGraphInList := function(G,list)
    for graph in list do
        if AreGraphsEqual(G,graph) eq true then
            return true;
        end if;
    end for;
    return false;
end function;

            
// Input:
//   list - a list of graphs
//
// Output:
//   A set containing all distinct graphs from the input list.
//
// Error:
//   empty_list_error - the input list must not be empty.
ListToSetForGraphs := function(list)                        
    graphset := {list[1]};
    for graph in list do
        if not IsGraphInSet(graph,graphset) then
            graphset := graphset join {graph};
        end if;
    end for;
    return graphset;
end function;

        
// Input:
//   G1, G2 - subgraphs of a graph G, each represented as a pair 
//            <vertex list, edge list>
//
// Output:
//   A subgraph given by the union of vertices as well as edges from G1 and G2.
//   The result is returned as a pair <vertex list, edge list>.
GetUnion := function(G1,G2)
    verticesInUnion := Set(G1[1]) join Set(G2[1]);
    edgesInUnion := Set(G1[2]) join Set(G2[2]);
    union := <SetToSequence(verticesInUnion),SetToSequence(edgesInUnion)>;
    return union;
end function;


// Input:
//   G_list - a list of subgraphs of a graph G,
//            where each subgraph is represented as a pair 
//            <vertex list, edge list>
//
// Output:
//   The union of all subgraphs in G_list.
//   The result is returned as a pair <vertex list, edge list>.
GetUnionList := function(G_list)
    verticesInUnion := &join[Set(x[1]) : x in G_list];
    edgesInUnion := &join[Set(x[2]) : x in G_list];
    union := <SetToSequence(verticesInUnion),SetToSequence(edgesInUnion)>;
    return union;
end function;

            
// Input:
//   subgraph - a subgraph of a graph G, represented as a pair
//              <vertex list, edge list>
//   f        - a morphism from graph G to an arbitrary graph H, 
//              represented as a pair <vertex mapping, edge mapping>, 
//              where:
//              - f[1] is a list of vertex mappings <v, vV_f>
//              - f[2] is a list of edge mappings <e, eE_f>
//
// Output:
//   A new subgraph of H representing the image of the input subgraph
//   under the morphism f. The image consists of all vertices and
//   edges that are images of the original vertices and edges.
ComputeImage := function(subgraph,f)                              
    verticesInImage := [];
    edgesInImage := [];
    for vertexMap in f[1] do
        if vertexMap[1] in subgraph[1] then
            verticesInImage := verticesInImage cat [vertexMap[2]];
        end if;
    end for;
    for edgeMap in f[2] do
        if edgeMap[1] in subgraph[2] then
            edgesInImage := edgesInImage cat [edgeMap[2]];
        end if;
    end for;
    verticesInImage := Set(verticesInImage);
    edgesInImage := Set(edgesInImage);
    verticesInImage := SetToSequence(verticesInImage);
    edgesInImage := SetToSequence(edgesInImage);
    return <verticesInImage,edgesInImage>;
end function;

// Input:
//   f - a graph morphism from G to H.
//   v - a vertex in the domain of vertex map of f.
// Output:
//   The vertex in H that v is mapped to.
GetImageOfVertex := function(f,v)
    for vmap in f[1] do
        if vmap[1] eq v then
            return vmap[2];
        end if;
    end for;
end function;

// Input:
//   f - a graph morphism from G to H.
//   e - an edge in the domain of edge map of f.
// Output:
//   The edge in H that e is mapped to.
GetImageOfEdge := function(f,e)
    for emap in f[2] do
        if emap[1][2] eq e then
            return emap[2];
        end if;
    end for;
end function;

// Input:
//   G - a graph represented as a pair <vertex list, edge list>.
//   v - a vertex in the vertex list of G.
// Output:
//   true if v is a vertex of a cyclic subgraph of G.
//   false if none of the cyclic subgraphs of G contains v.
IsVertexinCycle := function(G,v)
    reachable := Set([]);
    frontier := {v};
    while not IsEmpty(frontier) do
        nextFrontier := {};
        for u in frontier do
            for e in GetEdgesFromVertex(u,G) do
                w := e[3];
                if w notin reachable then
                    Include(~reachable,w);
                    Include(~nextFrontier,w);
                end if;
            end for;
        end for;
        frontier := nextFrontier;
    end while;
    if v in reachable then
        return true;
    else
        return false;
    end if;
end function;


// Input:
//   G - a graph represented as a pair <vertex list, edge list>
//
// Output:
//   A list of all cyclic subgraphs of G.
GiveCyclicSubgraphs := function(G)
    cyclicSubgraphs := [];
    n := #G[1];
    for i in [1..n] do
        C := CyclicGraph(i);
        list_inj := [f : f in ListGraphMorphisms(C,G) | Is_Injective(f,C,G)];
        if #list_inj ge 1 then
            cyclicSubgraphs cat:= SetToSequence(ListToSetForGraphs([ComputeImage(C,f) : f in list_inj]));
        end if;
    end for;
    return cyclicSubgraphs;
end function;


// Input:
//   G - a graph represented as a pair <vertex list, edge list>
//
// Output:
//   A list of all quasi-cyclic subgraphs of G.
GiveQuasiCyclicSubgraphs := function(G)
    quasiCyclicSubgraphs := GiveCyclicSubgraphs(G);
    for graph1 in quasiCyclicSubgraphs do
        for graph2 in quasiCyclicSubgraphs do
        if not IsEmpty(Set(graph1[1]) meet Set(graph2[1])) then 
            verticesNewGraph := Set(graph1[1]) join Set(graph2[1]);
            edgesNewGraph := Set(graph1[2]) join Set(graph2[2]);
            newGraph := <SetToSequence(verticesNewGraph),SetToSequence(edgesNewGraph)>;
            if not IsGraphInList(newGraph,quasiCyclicSubgraphs) then
            quasiCyclicSubgraphs := quasiCyclicSubgraphs cat [newGraph];
            end if;
        end if;
        end for;
    end for;
    return quasiCyclicSubgraphs;
end function;


// Input:
//   G - a graph represented as a pair <vertex list, edge list>
//
// Output:
//   true if G is quasi-cyclic;
//   false otherwise.
IsQuasiCyclicGraph := function(G)
    if #G[2] eq 0 then
        return false;
    end if;
    if IsStronglyConnected(G) then
        return true;
    else
        return false;
    end if;
end function;


// Input:
//   G - a quasi-cyclic graph
//
// Output:
//   The circumferential length of G
GetCircumferentialLength := function(G)
    allCyclicSubgraphs := Set(GiveCyclicSubgraphs(G));
    subsetCyclicSubgraphs := Subsets(allCyclicSubgraphs);
    lengthStacks := [];
    for set in subsetCyclicSubgraphs do
        if #&join{Set(C[2]) : C in set} eq #G[2] then
            length := &+[#C[2] : C in set];
            Append(~lengthStacks, length);
        end if;
    end for; 
    min := Min(lengthStacks);
    return min;
end function;

// Input:
//   G - a graph
//
// Output:
//   The adjacency matrix of G.
GetAdjacencyMatrix := function(G) // G: Graph
    n := #G[1];
    AG := ZeroMatrix(Rationals(), n, n);
    for e in G[2] do
        AG[Index(G[1],e[1]), Index(G[1],e[3])] := AG[Index(G[1],e[1]), Index(G[1],e[3])] + 1;
    end for;
    return AG;
end function;

// Input:
//   G - a graph
//   H - a graph
//   f - a graph morphism from G to H.
// Output:
//   
GetMorphismMatrix := function(G,H,f)
    Af := ZeroMatrix(Rationals(), #G[1], #H[1]);
    for vmap in f[1] do
        Af[Index(G[1], vmap[1]), Index(H[1],vmap[2])] := 1;
    end for;
    return Af;
end function;

// Input:
//   G - a quasi-cyclic graph represented as a pair <vertex list, edge list>
//   H - a graph represented as a pair <vertex list, edge list>
//   f - a morphism from graph G to H, represented as a pair <vertex mapping, edge mapping>,
//   where:
//     - f[1] is a list of vertex mappings <v, vV_f>
//     - f[2] is a list of edge mappings <e, eE_f>
// Output:
//   true if f is a quasi-isomorphism;
//   false otherwise.
IsQuasiisomorphism := function(G,H,f)
for e1 in G[2] do
    for e2 in G[2] do
        if e1 ne e2 then
            if e1[1] eq e2[1] and e1[3] eq e2[3] then
                if GetImageOfEdge(f,e1[2]) eq GetImageOfEdge(f,e1[2]) then
                    return false;
                end if;
            end if;
        end if;
    end for;
end for;

AG := GetAdjacencyMatrix(G);
AH := GetAdjacencyMatrix(H);
Af := GetMorphismMatrix(G,H,f);
n := #G[1] + #H[1];
for k in [1..n] do
    Diag_AG_k := DiagonalMatrix([(AG^k)[i][i] : i in [1..Nrows(AG^k)]]);
    Diag_AH_k := DiagonalMatrix([(AH^k)[i][i] : i in [1..Nrows(AH^k)]]);
    M := Transpose(Af) * Diag_AG_k  * Af - Diag_AH_k;
    if not IsZero(M) then
        return false;
    end if;
end for;

vertices := [];
edges := [];
for v_plus in G[1] do
    for v_minus in G[1] do
        if GetImageOfVertex(f,v_plus) eq GetImageOfVertex(f,v_minus) then
            Append(~vertices, <v_plus,v_minus>);
        end if;
    end for;
end for;
vertices := SetToSequence(Set(vertices));
for e_plus in G[2] do
    for e_minus in G[2] do
        if GetImageOfEdge(f,e_plus[2]) eq GetImageOfEdge(f,e_minus[2]) then
            Append(~edges, <<e_plus[1],e_minus[1]>,<e_plus[2],e_minus[2]>,<e_plus[3],e_minus[3]>>);
        end if;
    end for;
end for;
P_f_G := <vertices, edges>;
edges := SetToSequence(Set(edges));
for v in P_f_G[1] do
    if v[1] ne v[2] then
        if IsVertexinCycle(P_f_G,v) then
            return false;
        end if;
    end if;
end for;

return true; 
end function;

