Add DigraphMaximumMatching (#294)

* Add maximal mathing and maximum mathing methods and helper methods.

* Add documentation and examples for DigraphMaximalMatching and DigraphMaximumMatching as well as documentation for IsMaximumMatching.

* Fix typo, add examples and implementation for IsMaximumMatching

* Fix manual examples

* Add tests for matching methods, fix label bug in maximum matching code

* Add more tests for maximum matching, improve documentation and code

* Update oper.xml

Remove duplicate "that"

* Remove non-deterministic tests

* Fix tests

Co-authored-by: James Mitchell <[email protected]>

doc/attr.xml

@@ -2179,3 +2179,67 @@ true
   </Description>
 </ManSection>
 <#/GAPDoc>
+
+<#GAPDoc Label="DigraphMaximalMatching">
+<ManSection>
+  <Attr Name="DigraphMaximalMatching" Arg="digraph"/>
+  <Returns>A list of pairs of vertices.</Returns>
+  <Description>
+    This function returns a maximal matching of the digraph <A>digraph</A>.
+    <P/>
+
+    For the definition of a maximal matching, see <Ref Oper="IsMaximalMatching"/>. 
+
+    <Example><![CDATA[
+gap> D := DigraphFromDiSparse6String(".IeAoXCJU@|SHAe?d");
+<immutable digraph with 10 vertices, 13 edges>
+gap> M := DigraphMaximalMatching(D);; IsMaximalMatching(D, M);
+true
+gap> D := RandomDigraph(100);;
+gap> IsMaximalMatching(D, DigraphMaximalMatching(D));
+true
+gap> D := GeneralisedPetersenGraph(IsMutableDigraph, 9, 2);
+<mutable digraph with 18 vertices, 54 edges>
+gap> IsMaximalMatching(D, DigraphMaximalMatching(D));
+true
+]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>
+
+<#GAPDoc Label="DigraphMaximumMatching">
+<ManSection>
+  <Attr Name="DigraphMaximumMatching" Arg="digraph"/>
+  <Returns>A list of pairs of vertices.</Returns>
+  <Description>
+    This function returns a maximum matching of the digraph <A>digraph</A>.
+    <P/>
+
+    For the definition of a maximum matching, see <Ref Oper="IsMaximumMatching"/>.
+    If <A>digraph</A> is bipartite (see <Ref Oper="IsBipartiteDigraph"/>), then 
+    the algorithm used has complexity <C>O(m*sqrt(n))</C>. Otherwise for general 
+    graphs the complexity is <C>O(m*n*log(n))</C>. Here <C>n</C> is the number of vertices 
+    and <C>m</C> is the number of edges.
+
+    <Example><![CDATA[
+gap> D := DigraphFromDigraph6String("&I@EA_A?AdDp[_c??OO");
+<immutable digraph with 10 vertices, 23 edges>
+gap> M := DigraphMaximumMatching(D);; IsMaximalMatching(D, M);
+true
+gap> Length(M);
+5
+gap> D := Digraph([[5, 6, 7, 8], [6, 7, 8], [7, 8], [8], 
+>                  [], [], [], []]);;
+gap> M := DigraphMaximumMatching(D);
+[ [ 1, 5 ], [ 2, 6 ], [ 3, 7 ], [ 4, 8 ] ]
+gap> D := GeneralisedPetersenGraph(IsMutableDigraph, 9, 2);
+<mutable digraph with 18 vertices, 54 edges>
+gap> M := DigraphMaximumMatching(D);;
+gap> IsMaximalMatching(D, M);
+true
+gap> Length(M);
+9
+]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>

doc/oper.xml

@@ -1658,43 +1658,54 @@ true
 <ManSection>
   <Oper Name="IsMatching" Arg="digraph, list"/>
   <Oper Name="IsMaximalMatching" Arg="digraph, list"/>
+  <Oper Name="IsMaximumMatching" Arg="digraph, list"/>
   <Oper Name="IsPerfectMatching" Arg="digraph, list"/>
   <Returns><K>true</K> or <K>false</K>.</Returns>
   <Description>
     If <A>digraph</A> is a digraph and <A>list</A> is a list of pairs of
     vertices of <A>digraph</A>, then <C>IsMatching</C> returns <K>true</K> if
-    <A>list</A> is a matching of <A>digraph</A>.  The operations
-    <C>IsMaximalMatching</C> and <C>IsPerfectMatching</C> return <K>true</K> if
-    <A>list</A> is a maximal, or perfect, matching of <A>digraph</A>,
-    respectively.  Otherwise, these operations return <K>false</K>. <P/>
+    <A>list</A> is a matching of <A>digraph</A>.  The operation
+    <C>IsMaximalMatching</C> returns <K>true</K> if <A>list</A> is a maximal 
+    matching, <C>IsMaximumMatching</C> returns <K>true</K> if <A>list</A> is a
+    maximum matching and <C>IsPerfectMatching</C> returns <K>true</K> if 
+    <A>list</A> is a perfect, matching of <A>digraph</A>, respectively.  
+    Otherwise, each of these operations return <K>false</K>.
+    <P/>
 
     A <E>matching</E> <C>M</C> of a digraph <A>digraph</A> is a subset of the
     edges of <A>digraph</A>, i.e. <C>DigraphEdges(<A>digraph</A>)</C>, such
     that no pair of distinct edges in <C>M</C> are incident to the same vertex
     of <A>digraph</A>.  Note that this definition allows a matching to contain
     loops.  See <Ref Prop="DigraphHasLoops" />.  The matching <C>M</C> is
-    <E>maximal</E> if it is contained in no larger matching of the digraph, and
+    <E>maximal</E> if it is contained in no larger matching of the digraph,
+    is <E>maximum</E> if it has the greatest cardinality among all matchings and
     is <E>perfect</E> if every vertex of the digraph is incident to an edge in
-    the matching.  Every perfect matching is maximal.
+    the matching.  Every maximum or perfect matching is maximal. Note, however,
+    that not every perfect matching of digraphs with loops is maximum.
 
     <Example><![CDATA[
-gap> D := Digraph([[2], [1], [2, 3, 4], [3, 5], [1]]);
-<immutable digraph with 5 vertices, 8 edges>
+gap> D := Digraph([[1, 2], [1, 2], [2, 3, 4], [3, 5], [1]]);
+<immutable digraph with 5 vertices, 10 edges>
 gap> IsMatching(D, [[2, 1], [3, 2]]);
 false
 gap> edges := [[3, 2]];;
 gap> IsMatching(D, edges);
 true
 gap> IsMaximalMatching(D, edges);
 false
-gap> edges := [[5, 1], [3, 3]];;
+gap> edges := [[2, 1], [3, 4]];;
 gap> IsMaximalMatching(D, edges);
 true
 gap> IsPerfectMatching(D, edges);
 false
 gap> edges := [[1, 2], [3, 3], [4, 5]];;
 gap> IsPerfectMatching(D, edges);
 true
+gap> IsMaximumMatching(D, edges);
+false
+gap> edges := [[1, 1], [2, 2], [3, 3], [4, 5]];;
+gap> IsMaximumMatching(D, edges);
+true
 ]]></Example>
   </Description>
 </ManSection>

doc/z-chap4.xml

@@ -15,6 +15,8 @@
     <#Include Label="DigraphOutEdges">
     <#Include Label="IsDigraphEdge">
     <#Include Label="IsMatching">
+    <#Include Label="DigraphMaximalMatching">
+    <#Include Label="DigraphMaximumMatching">
   </Section>
 
   <Section><Heading>Neighbours and degree</Heading>

gap/attr.gd

@@ -114,3 +114,6 @@ DeclareAttribute("DigraphMycielskianAttr", IsDigraph);
 
 DeclareAttribute("DigraphCartesianProductProjections", IsDigraph);
 DeclareAttribute("DigraphDirectProductProjections", IsDigraph);
+
+DeclareAttribute("DigraphMaximalMatching", IsDigraph);
+DeclareAttribute("DigraphMaximumMatching", IsDigraph);

gap/attr.gi

@@ -2044,3 +2044,237 @@ InstallMethod(DigraphMycielskian,
 
 InstallMethod(DigraphMycielskianAttr, "for an immutable digraph",
 [IsImmutableDigraph], DigraphMycielskian);
+
+# Uses a simple greedy algorithm.
+BindGlobal("DIGRAPHS_MaximalMatching",
+function(D)
+  local mate, u, v;
+  mate := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+  for v in DigraphVertices(D) do
+    if mate[v] = 0 then
+      for u in OutNeighboursOfVertex(D, v) do
+        if mate[u] = 0 then
+          mate[u] := v;
+          mate[v] := u;
+          break;
+        fi;
+      od;
+    fi;
+  od;
+  return mate;
+end);
+
+# For bipartite digraphs implements the Hopcroft-Karp matching algorithm,
+# complexity O(m*sqrt(n))
+BindGlobal("DIGRAPHS_BipartiteMatching",
+function(D, mate)
+  local U, dist, inf, dfs, bfs, u;
+
+  U := DigraphBicomponents(D);
+  U := U[PositionMinimum(U, Length)];
+
+  bfs := function()
+    local v, que, q;
+    que := [];
+    for v in U do
+      if mate[v] = inf then
+        dist[v] := 0;
+        Add(que, v);
+      else
+        dist[v] := inf;
+      fi;
+    od;
+    dist[inf] := inf;
+
+    q := 1;
+    while q <= Length(que) do
+      if dist[que[q]] < dist[inf] then
+        for v in OutNeighborsOfVertex(D, que[q]) do
+          if dist[mate[v]] = inf then
+            dist[mate[v]] := dist[que[q]] + 1;
+            Add(que, mate[v]);
+          fi;
+        od;
+      fi;
+      q := q + 1;
+    od;
+    return dist[inf] <> inf;
+  end;
+
+  dfs := function(u)
+    local v;
+    if u = inf then
+      return true;
+    fi;
+    for v in OutNeighborsOfVertex(D, u) do
+      if dist[mate[v]] = dist[u] + 1 and dfs(mate[v]) then
+        mate[v] := u;
+        mate[u] := v;
+        return true;
+      fi;
+    od;
+    dist[u] := inf;
+    return false;
+  end;
+
+  inf  := DigraphNrVertices(D) + 1;
+  dist := ListWithIdenticalEntries(inf, inf);
+  for u in [1 .. Length(mate)] do
+    if mate[u] = 0 then
+      mate[u] := inf;
+    fi;
+  od;
+
+  while bfs() do
+    for u in U do
+      if mate[u] = inf then dfs(u);
+      fi;
+    od;
+  od;
+
+  for u in [1 .. DigraphNrVertices(D)] do
+    if mate[u] = inf then
+      mate[u] := 0;
+    fi;
+  od;
+
+  return mate;
+end);
+
+# For general digraphs implements a modified version of Gabow's maximum matching
+# algorithm, complexity O(m*n*log(n)).
+BindGlobal("DIGRAPHS_GeneralMatching",
+function(D, mate)
+  local blos, pred, time, t, tj, u, dfs, mark, blosfind;
+
+  blosfind := function(x)
+    if x <> blos[x] then
+      blos[x] := blosfind(blos[x]);
+    fi;
+    return blos[x];
+  end;
+
+  mark := function(v, x, b, path)
+    while blosfind(v) <> b do
+      pred[v] := x;
+      x       := mate[v];
+      Add(tj, v);
+      Add(tj, x);
+      if time[x] = 0 then
+        t       := t + 1;
+        time[x] := t;
+        Add(path, x);
+      fi;
+      v := pred[x];
+    od;
+  end;
+
+  dfs := function(v)
+    local x, bx, bv, b, y, z, path;
+    for x in OutNeighboursOfVertex(D, v) do
+      bv := blosfind(v);
+      bx := blosfind(x);
+      if bx <> bv then
+        if time[x] > 0 then
+          path := [];
+          tj   := [];
+          if time[bx] < time[bv] then
+            b := bx;
+            mark(v, x, b, path);
+          else
+            b := bv;
+            mark(x, v, b, path);
+          fi;
+          for z in tj do
+            blos[z] := b;
+          od;
+          for z in path do
+            if dfs(z) then
+              return true;
+            fi;
+          od;
+        elif pred[x] = 0 then
+          pred[x] := v;
+          if mate[x] = 0 then
+            while x <> 0 do
+              y       := pred[x];
+              v       := mate[y];
+              mate[y] := x;
+              mate[x] := y;
+              x       := v;
+            od;
+            return true;
+          fi;
+          if time[mate[x]] = 0 then
+            t             := t + 1;
+            time[mate[x]] := t;
+            if dfs(mate[x]) then
+              return true;
+            fi;
+          fi;
+        fi;
+      fi;
+    od;
+    return false;
+  end;
+
+  time := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+  blos := [1 .. DigraphNrVertices(D)];
+  pred := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+  t    := 0;
+  for u in DigraphVertices(D) do
+    if mate[u] = 0 then
+      t       := t + 1;
+      time[u] := t;
+      if dfs(u) then
+        time := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+        blos := [1 .. DigraphNrVertices(D)];
+        pred := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+      fi;
+    fi;
+  od;
+
+  return mate;
+end);
+
+BindGlobal("DIGRAPHS_MateToMatching",
+function(D, mate)
+  local u, M;
+  M := [];
+  for u in DigraphVertices(D) do
+    if u <= mate[u] then
+      if IsDigraphEdge(D, u, mate[u]) then
+        Add(M, [u, mate[u]]);
+      elif IsDigraphEdge(D, mate[u], u) then
+        Add(M, [mate[u], u]);
+      fi;
+    fi;
+  od;
+  return Set(M);
+end);
+
+InstallMethod(DigraphMaximalMatching, "for a digraph", [IsDigraph],
+D -> DIGRAPHS_MateToMatching(D, DIGRAPHS_MaximalMatching(D)));
+
+InstallMethod(DigraphMaximumMatching, "for a digraph", [IsDigraph],
+function(D)
+  local mateG, mateD, G, M, i, lab;
+  G     := DigraphImmutableCopy(D);
+  G     := InducedSubdigraph(G, Difference(DigraphVertices(G), DigraphLoops(G)));
+  lab   := DigraphVertexLabels(G);
+  G     := DigraphSymmetricClosure(G);
+  mateG := DIGRAPHS_MaximalMatching(G);
+  if IsBipartiteDigraph(G) then
+    mateG := DIGRAPHS_BipartiteMatching(G, mateG);
+  else
+    mateG := DIGRAPHS_GeneralMatching(G, mateG);
+  fi;
+  mateD := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
+  for i in DigraphVertices(G) do
+    if mateG[i] <> 0 then
+      mateD[lab[i]] := lab[mateG[i]];
+    fi;
+  od;
+  M := List(DigraphLoops(D), x -> [x, x]);
+  return Union(M, DIGRAPHS_MateToMatching(D, mateD));
+end);

gap/oper.gd

@@ -88,6 +88,7 @@ DeclareOperation("IsUndirectedSpanningForest", [IsDigraph, IsDigraph]);
 
 DeclareOperation("IsMatching", [IsDigraph, IsHomogeneousList]);
 DeclareOperation("IsMaximalMatching", [IsDigraph, IsHomogeneousList]);
+DeclareOperation("IsMaximumMatching", [IsDigraph, IsHomogeneousList]);
 DeclareOperation("IsPerfectMatching", [IsDigraph, IsHomogeneousList]);
 
 # 9. Connectivity . . .