Implementing Colour Refinement

doc/oper.xml

@@ -2619,3 +2619,48 @@ true]]>
   gap> DIGRAPHS_FREE_CLIQUES_DATA();
   ]]></Example>
 <#/GAPDoc>
+
+<#GAPDoc Label="DigraphColourRefinement">
+<ManSection>
+  <Oper Name="DigraphColourRefinement"
+        Arg="digraph"
+        Label="for a digraph"/>
+  <Returns>A list of integers.</Returns>
+  <Description>
+    Colour refinement is a method of colouring a digraph such that it has a 'stable colouring'. That is, for a 
+    colour, every node with that colour has an identical configuration of coloured neighbours. This means that 
+    all nodes with the same colour have equal numbers of neighbours of each colour.
+    <C>DigraphColourRefinement</C> considers the out-neighbours and in-neighbours of a node separately,
+    meaning the nodes of a certain colour must have an equal number of out-neighbours of each colour, 
+    and an equal number of in-neighbours of each colour. 
+    <P/>
+    <C>DigraphColourRefinement</C> returns the colouring as a list where the value at the ith position is the 
+    colour of node i. The time complexity of this algorithm is <M>O(n^2 log n)</M>.
+    <P/>
+    Because the labels of the vertices are not used in any capacity during the refinement process (i.e. to 
+    determine which cells to refine), the colouring produced is canonical. This means that for two isomorphic 
+    digraphs, they would produce the same colouring. Identical colourings do not necessarily prove that two 
+    digraphs are isomorphic, but non-identical colourings would prove that they aren't. While the colouring 
+    produced for two isomorphic digraphs will be identical, the output of <C>DigraphColourRefinement</C> will not,
+    as the list produced is ordered by vertex labels. This is demonstrated in the examples below.
+    <P/>
+    See also
+    <Ref Oper="DigraphColouring"
+         Label="for a digraph and a number of colours"/>.
+    <P/>
+
+    <Example><![CDATA[
+gap> D := Digraph([[3], [], [1, 9], [], [10], [7, 8, 9], [6, 8],
+> [6, 7], [3, 6, 10], [5, 9]]);;
+gap> DigraphColourRefinement(D);
+[ 2, 1, 4, 1, 2, 5, 3, 3, 6, 4 ]
+gap> D2 := Digraph([[6, 7], [], [8, 9], [], [10], [1, 7, 9], [6, 1],
+> [3], [3, 6, 10], [5, 9]]);;
+gap> DigraphColourRefinement(D2);
+[ 3, 1, 4, 1, 2, 5, 3, 2, 6, 4 ]
+gap> DigraphColourRefinement(Digraph([[], [1], [1], [1]]));
+[ 1, 2, 2, 2 ]
+]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>

doc/z-chap6.xml

@@ -46,6 +46,7 @@ from} $E_a$ \emph{to} $E_b$. In this case we say that $E_a$ and $E_b$ are
     <#Include Label="IsDigraphColouring">
     <#Include Label="MaximalCommonSubdigraph">
     <#Include Label="MinimalCommonSuperdigraph">
+    <#Include Label="DigraphColourRefinement">
   </Section>
 
   <Section><Heading>Homomorphisms of digraphs</Heading>

gap/oper.gd

@@ -156,6 +156,8 @@ DeclareOperation("Dominators", [IsDigraph, IsPosInt]);
 DeclareOperation("DominatorTree", [IsDigraph, IsPosInt]);
 DeclareOperation("DigraphCycleBasis", [IsDigraph]);
 
+DeclareOperation("DigraphColourRefinement", [IsDigraph]);
+
 # 10. Operations for vertices . . .
 DeclareOperation("PartialOrderDigraphJoinOfVertices",
                  [IsDigraph, IsPosInt, IsPosInt]);

gap/oper.gi

@@ -2732,3 +2732,149 @@ function(D, n)
   od;
   return kings;
 end);
+
+InstallMethod(DigraphColourRefinement, "for a digraph", [IsDigraph],
+function(D)
+
+  local i, cMin, cMax, Q, q, C, CD, j, P, v, outNB, inNB, colourCells,
+  pair, current, currentPair, newSet, colour, largest, recolour, cell,
+  colourCell, toAdd, DVertices, DNrVertices, outNeighboursD, inNeighboursD;
+
+  if DigraphHasLoops(D) then
+    ErrorNoReturn("the digraph cannot contain loops");
+  fi;
+
+  outNeighboursD := OutNeighbours(D);
+  inNeighboursD := InNeighbours(D);
+
+  DNrVertices := DigraphNrVertices(D);
+  DVertices := DigraphVertices(D);
+
+  cMin := 1;
+  cMax := 1;
+
+  # Queue of colours
+  Q := [1];
+
+  # Initial colouring
+  # vertices -> colour
+  C := ListWithIdenticalEntries(DNrVertices, 1);
+
+  # Colour classes
+  # All vertices initialised to 1
+  # colour -> vertices labelled as such
+  P := [];
+  P[1] := [1 .. DNrVertices];
+
+  while not IsEmpty(Q) do
+
+    # Pop colour off Q
+    q := Remove(Q, 1);
+
+    # For each v (vertices) in D:
+    # Get the neighbours of v that are in the colour class q
+    outNB := EmptyPlist(DNrVertices);
+    inNB := EmptyPlist(DNrVertices);
+
+    for v in DVertices do
+      outNB[v] := Intersection(outNeighboursD[v], P[q - cMin + 1]);
+      inNB[v] := Intersection(inNeighboursD[v], P[q - cMin + 1]);
+    od;
+
+    CD := List(DVertices, v -> [C[v], Length(outNB[v]), Length(inNB[v]), v]);
+
+    Sort(CD);
+
+    colourCells := [];
+    currentPair := [];
+    colourCell := [];
+    cell := [];
+
+    recolour := false;
+
+    j := 0;
+
+    # Creating the multiset of colours of neighbours:
+    # This multiset is used to determine whether to refine the colouring
+    # and if so, which cells to split. The label of each vertex has no bearing
+    # in either of these decisions (as they are irrelevant to the creation or
+    # sorting of the multiset) - meaning the colouring produced is canonical.
+    for pair in CD do
+      current := [pair[1], pair[2], pair[3]];
+
+      # If different colour reached:
+      if currentPair <> [] and current[1] <> currentPair[1] then
+        Add(colourCell, cell);
+        Add(colourCells, colourCell);
+        cell := [pair[4]];
+        colourCell := [];
+      else
+        # If first iteration, or same neighbour configuration:
+        if currentPair = [] or current = currentPair then
+          Add(cell, pair[4]);
+        else
+          # If same colour, but different neighbour configuration:
+          recolour := true;
+          Add(colourCell, cell);
+          cell := [pair[4]];
+        fi;
+      fi;
+
+      currentPair := current;
+    od;
+    Add(colourCell, cell);
+    Add(colourCells, colourCell);
+
+    # If there is reason for recolouring
+    if recolour then
+
+      # Clearing P and Q
+      P := [];
+      Q := [];
+
+      # Add to Q
+      toAdd := cMax;
+
+      for i in [1 .. Length(colourCells)] do
+
+        colourCell := colourCells[i];
+
+        # Determine the largest cell for that colour
+        largest := PositionMaximum(colourCell, Length);
+
+        # Add all new colours except that corresponding to the largest cell
+        for j in [1 .. Length(colourCell)] do
+          toAdd := toAdd + 1;
+          if j <> largest then
+            Add(Q, toAdd);
+          fi;
+        od;
+
+      od;
+
+      # Updating colours for the next round
+      cMin := cMax + 1;
+      cMax := toAdd;
+
+      colour := cMin;
+
+      # Updating C and P
+      for i in [1 .. Length(colourCells)] do
+        for j in [1 .. Length(colourCells[i])] do
+          Add(P, []);
+          for v in colourCells[i][j] do
+            C[v] := colour;
+            Add(P[Length(P)], v);
+          od;
+
+          colour := colour + 1;
+
+        od;
+      od;
+    fi;
+
+  od;
+
+  return C - (cMin - 1);
+
+end);

tst/standard/oper.tst

@@ -3324,6 +3324,38 @@ gap> DigraphEdges(D);
 gap> DigraphVertexLabels(D);
 [ 1, 2, 3, 6, [ 4, 5 ] ]
 
+# DigraphColourRefinement
+gap> D := Digraph([[3], [], [1, 9], [], [10], [7, 8, 9], [6, 8], [6, 7], [3, 6, 10], [5, 9]]);;
+gap> DigraphColourRefinement(D);
+[ 2, 1, 4, 1, 2, 5, 3, 3, 6, 4 ]
+gap> D := Digraph([[], [1], [1], [1]]);;
+gap> DigraphColourRefinement(D);
+[ 1, 2, 2, 2 ]
+gap> D := Digraph([[1], [1], [1], [1]]);;
+gap> DigraphColourRefinement(D);
+Error, the digraph cannot contain loops
+gap> D := Digraph([[], [], [], []]);;
+gap> DigraphColourRefinement(D);
+[ 1, 1, 1, 1 ]
+gap> D := Digraph([[2], [3], [2, 4], [2, 5], [4, 6], [5]]);;
+gap> DigraphColourRefinement(D);
+[ 1, 3, 4, 5, 6, 2 ]
+gap> D := Digraph([[2], [3], [1]]);;
+gap> DigraphColourRefinement(D);
+[ 1, 1, 1 ]
+gap> D := Digraph([[2, 4], [5], [2, 4], [5], [1, 3]]);;
+gap> DigraphColourRefinement(D);
+[ 2, 1, 2, 1, 3 ]
+gap> D := Digraph([[4], [1, 3], [4], [5], [1, 3]]);;
+gap> DigraphColourRefinement(D);
+[ 2, 3, 2, 1, 4 ]
+gap> D := Digraph([]);;
+gap> DigraphColourRefinement(D);
+[  ]
+gap> D := Digraph([[2, 3, 4, 5], [], [], [], []]);;
+gap> DigraphColourRefinement(D);
+[ 2, 1, 1, 1, 1 ]
+
 #
 gap> DIGRAPHS_StopTest();
 gap> STOP_TEST("Digraphs package: standard/oper.tst", 0);