/* This file contains the MAGMA-code for the function NIST (and all the neccesary subfunctions) that with input of a graph G, computes all the non-isomorphic spanning trees of G. In the input graph G, the edges has to be labeled with integers. A more explicit explanation can be found in "Non-isomorphic spanning trees of graphs" by Janneke van den Boomen. The two main subfunctions of NIST are AlgorithmS and TreeClasses. */ //All subfunctions of AlgorithmS Indices := func< S, e | [ i : i in [1..#S] | S[i] eq e ] >; /* Get the table as in Knuth 7.2.1.6, page 24 for a given graph. The information of arc i is on the (i+1)th place in the sequences. The edges of the input graph must have positive, integral labels (but none of it the label 0) */ table := function(Graph) V := Vertices(Graph); n := #(V); /* construct a (labeled) digraph D in which there's an arc from u to v and from v to u iff u and v are adjacent in Graph*/ D := Digraph; VD := Vertices(D); label := n; I := []; //for the construction of the sequence I of all I_a's for i := 1 to n do I[i] := 0; end for; for i in V do N := Neighbours(Vertices(Graph)!i); for j in N do AddEdge(~D,VD!i,VD!j); AssignLabel(Edges(D)![VD!i,VD!j],[label+1]); AddEdge(~D,VD!j,VD!i); AssignLabel(Edges(D)![VD!j,VD!i],[label]); /* construct the sequence I of all I_a's */ I[label+1] := EdgeLabel(Graph,Index(Vertices(D),i),Index(Vertices(D),j)); I[label+2] := EdgeLabel(Graph,Index(Vertices(D),i),Index(Vertices(D),j)); label := label + 2; RemoveEdge(~Graph,i,j); end for; end for; /* construct the sequence T of all t_a's */ T := []; for i := 1 to n do T[i] := 0; end for; for e in Edges(D) do i := EdgeLabels(D)[Index(Edges(D),e)][1]; T[i+1] := Index(Vertices(D), EndVertices(Edges(D)!e)[2]); end for; /* construct the sequences N and P of all n_a's and p_a's */ P := []; N := []; for v := 1 to #Vertices(D) do out := OutEdges(Vertices(D)!v); if #out ge 2 then //for edge v-1 N[v] := EdgeLabels(D)[Index(Edges(D),out[#out])][1]; P[v] := EdgeLabels(D)[Index(Edges(D),out[1])][1]; //for the first edge in out N[EdgeLabels(D)[Index(Edges(D),out[1])][1]+1] := v-1; P[EdgeLabels(D)[Index(Edges(D),out[1])][1]+1] := EdgeLabels(D)[Index(Edges(D),out[2])][1]; //for the last edge in out N[EdgeLabels(D)[Index(Edges(D),out[#out])][1]+1] := EdgeLabels(D)[Index(Edges(D),out[#out-1])][1]; P[EdgeLabels(D)[Index(Edges(D),out[#out])][1]+1] := v-1; //all the other edges for i := 2 to (#out-1) do N[EdgeLabels(D)[Index(Edges(D),out[i])][1]+1] := EdgeLabels(D)[Index(Edges(D),out[i-1])][1]; P[EdgeLabels(D)[Index(Edges(D),out[i])][1]+1] := EdgeLabels(D)[Index(Edges(D),out[i+1])][1]; end for; else //since our graph has to be connnected, #out can't be zero, so in this case its one //for edge v N[v] := EdgeLabels(D)[Index(Edges(D),out[1])][1]; P[v] := EdgeLabels(D)[Index(Edges(D),out[1])][1]; //for the edge in out N[EdgeLabels(D)[Index(Edges(D),out[1])][1]+1] := v-1; P[EdgeLabels(D)[Index(Edges(D),out[1])][1]+1] := v-1; end if; end for; return T,I,N,P; end function; /*BridgeTest is a function that for a given graph and edge, tests whether that edge is a bridge of the graph or not. The algorithm is based on exercise 95 of Knuth 7.2.1.6 */ BridgeTest := function(T,I,N,e,n); B := []; for i := 1 to n do B[i] := 0; end for; //(i) u := T[e+1]; w := u; e1 := Exclude(Indices(I,I[e+1]),e+1)[1]-1; v := T[e1+1]; B[w] := v; //(iii) while u ne v do //(ii) f := N[u-1+1]; while T[f+1] ne 0 do //a v2 := T[f+1]; if B[v2] eq 0 then //b if v2 ne v then B[v2] := v; B[w] := v2; w := v2; else //c if f ne e1 then //(v) u := T[e+1]; while u ne v do w := B[u]; B[u] := 0; u := w; end while; u := T[e+1]; return "false"; break; end if; end if; end if; //d f := N[f+1]; end while; //(iii) u := B[u]; end while; //(iv); while u ne v do w := B[u]; B[u] := 0; u := w; end while; return "true"; end function; //------------------------------------------------------------------------------------------------- /* Algorithm S determines all spanning trees of a given (labeled) graph. The algorithm is based on the algorithm on page 25 of Knuth 7.2.1.6 */ AlgorithmS := function(Graph) T,I,N,P := table(Graph); Tbegin := T; //we need the original T to transform a sequence A into a real tree (subtree of Graph) n := #Vertices(Graph); D := []; for i := 1 to n do D[i] := Degree(Vertices(Graph)!i); end for; AllTrees := []; //S1 ST := Edges(SpanningTree(Graph)); //a spanning tree in Graph A := []; for i := 1 to #ST do //transform ST into a sequence A with the edgelabes in the directed graph A[i] := Index(I,EdgeLabels(Graph)[Index(Edges(Graph), ST[i])])-1; end for; x := 0; l := 1; S := [0]; L := []; //for the l_f's while true do if n eq 2 then v := 1; e := N[1]; else //S4 while l lt n-1 do //S2 e := A[l+1]; e1 := Exclude(Indices(I,I[e+1]),e+1)[1]-1; if D[T[e+1]] le D[T[e1+1]] then //<- "Index", gives the first index, so we can use e+1 u := T[e+1]; v := T[e1+1]; else v := T[e+1]; u := T[e1+1]; e := e1; end if; //S3 k := D[u]+D[v]; f := N[u-1+1]; g := 0; while T[f+1] ne 0 do if T[f+1] eq v then //Delete(f) N[P[f+1]+1] := N[f+1]; P[N[f+1]+1] := P[f+1]; //(f1:=) f+1 := Exclude(Indices(I,I[f]),f)[1]-1 (so: all places in I of the label of f, accept place f) //Delete(f+1) f1 := Exclude(Indices(I,I[f+1]),f+1)[1]-1; N[P[f1+1]+1] := N[f1+1]; P[N[f1+1]+1] := P[f1+1]; k := k-2; L[f] := g; g := f; else f1 := Exclude(Indices(I,I[f+1]),f+1)[1]-1; T[f1+1] := v; end if; f := N[f+1]; end while; L[e] := g; D[v] := k; g := v-1; N[P[f+1]+1] := N[g+1]; P[N[g+1]+1] := P[f+1]; P[N[f+1]+1] := g; N[g+1] := N[f+1]; A[l] := e; //S4 l := l+1; if l lt n-1 then S[l] := 0; end if; end while; e := N[v-1+1]; end if; //S5 A[n-1] := e; //we want to append a tree, not the sequence A Tree := Graph; for i := 1 to #A do vertex1 := Tbegin[A[i]+1]; A1 := Exclude(Indices(I,I[A[i]+1]), A[i]+1)[1]-1; vertex2 := Tbegin[A1+1]; AddEdge(~Tree,vertex1, vertex2); end for; Append(~AllTrees, Tree); x := e; e := N[e+1]; while T[e+1] ne 0 do A[n-1] := e; //we want to append a tree, not the sequence A Tree := Graph; for i := 1 to #A do vertex1 := Tbegin[A[i]+1]; A1 := Exclude(Indices(I,I[A[i]+1]), A[i]+1)[1]-1; vertex2 := Tbegin[A1+1]; AddEdge(~Tree,vertex1, vertex2); end for; Append(~AllTrees, Tree); x := e; e := N[e+1]; end while; //S6 testvar := 0; while testvar eq 0 do l := l-1; if l eq 0 then return AllTrees; break; else e := A[l]; u := T[e+1]; e1 := Exclude(Indices(I,I[e+1]),e+1)[1]-1; v := T[e1+1]; //S7 f := u-1; g := v-1; N[g+1] := N[P[f+1]+1]; P[N[g+1]+1] := g; N[P[f+1]+1] := f; P[N[f+1]+1] := f; f := P[f+1]; // f := P[f+1]; while T[f+1] ne 0 do f1 := Exclude(Indices(I,I[f+1]),f+1)[1]-1; T[f1+1] := u; f := P[f+1]; end while; f := L[e]; k := D[v]; while f ne 0 do k := k+2; //undelete(f+1) f1 := Exclude(Indices(I,I[f+1]),f+1)[1]-1; P[N[f1+1]+1] := f1; N[P[f1+1]+1] := f1; //undelete(f) P[N[f+1]+1] := f; N[P[f+1]+1] := f; f := L[f]; end while; D[v] := k - D[u]; //S8 if BridgeTest(T,I,N,e,n) eq "true" then //S9 e := S[l]; while e ne 0 do u := T[e+1]; e1 := Exclude(Indices(I,I[e+1]),e+1)[1]-1; v := T[e1+1]; D[u] := D[u] + 1; D[v] := D[v] +1; //undelete(e+1) P[N[e1+1]+1] := e1; N[P[e1+1]+1] := e1; //undelete(e) P[N[e+1]+1] := e; N[P[e+1]+1] := e; e := L[e]; end while; //S8 otherwise else x := e; L[e] := S[l]; S[l] := e; //delete(e) N[P[e+1]+1] := N[e+1]; P[N[e+1]+1] := P[e+1]; //delete(e+1) e1 := Exclude(Indices(I,I[e+1]),e+1)[1]-1; N[P[e1+1]+1] := N[e1+1]; P[N[e1+1]+1] := P[e1+1]; D[u] := D[u] -1; D[v] := D[v] -1; testvar := 1; end if; end if; end while; end while; end function; //------------------------------------------------------------------------------------------------- //------------------------------------------------------------------------------------------------- //All subfunctions of TreeClasses /*The next algorithm determines whether the two input-trees are isomorphic or not. */ TreeIsomorphism := function(Tree1, Tree2); //Copy the trees, since otherwise we really label the trees T1 := Tree1; T2 := Tree2; // Step 1. V1 := Vertices(T1); Labeled1 := {}; for v in V1 do if Degree(v) eq 1 then AssignLabel(V1!v, [1]); Labeled1 := Labeled1 join {v}; end if; end for; V2 := Vertices(T2); Labeled2 := {}; for v in V2 do if Degree(v) eq 1 then AssignLabel(V2!v, [1]); Labeled2 := Labeled2 join {v}; end if; end for; V1 := (V1 diff Labeled1); V2 := (V2 diff Labeled2); //Since all labels have to be different, we have to remember the last label used lastlabel := 1; if #Labeled1 ne #Labeled2 then return "Non-isomorphic"; else while #V1 ne 0 or #V2 ne 0 do //V1 are the not-labeled vertices of T1; Step 4 // Step 2. S1 := {}; LabelsS1 := []; for v in V1 do neighbours := Neighbours(v); if #(neighbours diff Labeled1) le 1 then S1 := S1 join {v}; end if; end for; for s in S1 do Label_s := []; N := Neighbours(VertexSet(T1)!s); ToLabel := Neighbours(VertexSet(T1)!s) diff V1; //Delete that non-labeled vertex for n in ToLabel do Append(~Label_s, VertexLabel(T1,Index(Vertices(T1), n) )[1]); end for; Label_s := Sort(Label_s); AssignLabel(Vertices(T1)!s, Label_s); LabelsS1 := LabelsS1 cat [ Label_s]; end for; LabelsS1 := Sort(LabelsS1); S2 := {}; LabelsS2 := []; for v in V2 do neighbours := Neighbours(v); if #(neighbours diff Labeled2) le 1 then S2 := S2 join {v}; end if; end for; for s in S2 do Label_s := []; for n in (Neighbours(s) diff V2) do Append(~Label_s, VertexLabel(T2,Index(Vertices(T2), n) )[1]); end for; Label_s := Sort(Label_s); AssignLabel(Vertices(T2)!s, Label_s); LabelsS2:= LabelsS2 cat [Label_s]; end for; LabelsS2 := Sort(LabelsS2); if LabelsS1 ne LabelsS2 then return "Non-isomorphic"; else //Step 3 //Delete multiple occuring labels from S1 -> new S1 i := 1; while i lt #LabelsS1 do if LabelsS1[i] eq LabelsS1[i+1] then Remove(~LabelsS1,i); else i := i+1; end if; end while; //Relabel T1 for e in S1 do AssignLabel(Vertices(T1)!e, [lastlabel + Index(LabelsS1, VertexLabel(T1, Index(Vertices(T1), e)))]); end for; /* The LabelsS1 without multiple occuring labels is the same as LabelsS2 without multiple occuring labels, so we can use the new LabelsS1 to relabel T2 */ for e in S2 do AssignLabel(Vertices(T2)!e, [lastlabel + Index(LabelsS1, VertexLabel(T2, Index(Vertices(T2), e)))]); end for; lastlabel := lastlabel + #LabelsS1; end if; Labeled1 := Labeled1 join S1; V1 := (V1 diff Labeled1); Labeled2 := Labeled2 join S2; V2 := (V2 diff Labeled2); end while; end if; return "Isomorphic"; end function; /* The algorithm NewClasses decides whether a given tree is isomorphic to one of the trees in Classes, or will be a new class. */ NewClasses := function(Tree, Classes) DegreeTree := DegreeSequence(Tree); for Trees in Classes do DegreeTrees := DegreeSequence(Trees); if DegreeTree eq DegreeTrees then //we use the degreesequence as a first selection if TreeIsomorphism(Tree, Trees) eq "Isomorphic" then return Classes; break; end if; end if; end for; Include(~Classes, Tree); return Classes; end function; //------------------------------------------------------------------------------------------------- /*Transform the sequence AllTrees (which is the output of AlgorithmS) into a sequence of not-isomoforphic trees.*/ TreeClasses := function(AllTrees) RealClasses := [AllTrees[1]]; for i := 2 to #AllTrees do RealClasses := NewClasses(AllTrees[i], RealClasses); end for; return RealClasses; end function; //-------------------------------------------------------------------------------------------------- //-------------------------------------------------------------------------------------------------- NIST := function(Graph) AllTrees := AlgorithmS(Graph); NIStrees := TreeClasses(AllTrees); return NIStrees; end function; /* If the edges of your original graph has no labeles, you can label them with the next procedure: */ LabelGraphedges := procedure(Graph) for e in Edges(Graph) do AssignLabel(Edges(Graph)!e, Index(Edges(Graph),e)); end for; end procedure; //-------------------------------------------------------------------------------------------------- //Nog enkele andere functies?? NISTmetRijtjes of zo