From e2bbdcd907b025f90a906cfc94a6b4362511ad3f Mon Sep 17 00:00:00 2001 From: Ben Wolfe Date: Wed, 17 Jun 2020 16:46:09 -0400 Subject: [PATCH 1/5] Optimized depth-first maximum matching --- .../uk/ac/ebi/beam/KececiogluMatching.java | 347 ++++++++++++++++++ 1 file changed, 347 insertions(+) create mode 100644 core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java diff --git a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java new file mode 100644 index 0000000..f9ce4b1 --- /dev/null +++ b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java @@ -0,0 +1,347 @@ +package uk.ac.ebi.beam; + +import java.util.ArrayDeque; +import java.util.Arrays; +import java.util.Deque; + +final class KececiogluMatching { + + private final static int EvenLabel = 2; + private final static int OddLabel = 3; + private final static int UnreachedLabel = 4; + private final static int nil = -1; + + private final Edge[] Tree; + private final Edge[] Bridge; + private final int[] Shore; + private final int[] Base; + private final int[] Label; + private final int[] Age; + + private final ArrayDeque Path; + private final ArrayDeque alternatingTree; + + private final int nMatched; + private final Matching matching; + private final IntSet subset; + private int Time; + + private boolean IsReached(int V) {return Label[V] != UnreachedLabel;} + private boolean IsEven(int V) {return Label[V] == EvenLabel;} + private boolean IsOdd(int V) {return Label[V] == OddLabel;} + private int Other(Edge E, int V) {return ((E) == null) ? nil : E.other(V);} + + /** + * + * Maximum matching in general graphs using Edmond's Blossom Algorithm. This + * implementation was adapted from John Kececioglu and Justin Pecqueur. + * "Computing maximum-cardinality matchings in sparse general graphs." + * Proceedings of the 2nd Workshop on Algorithm Engineering, 121-132, 1998. + * and from their C implementation at + * https://www2.cs.arizona.edu/~kece/Research/code/matching1.sh.Z + * + * @author Ben Wolfe + */ + KececiogluMatching(Graph G, Matching M, int nMatched, IntSet subset) + { + this.Time = 1; + this.Label = new int[G.order()]; + this.Age = new int[G.order()]; + this.Shore = new int[G.order()]; + this.Tree = new Edge[G.order()]; + this.Bridge = new Edge[G.order()]; + this.Base = new int[G.order()]; + this.alternatingTree = new ArrayDeque(G.order()); + Path = new ArrayDeque(G.order()); + + + this.matching = M; + this.subset =subset; + + Arrays.fill(Base, nil); + Arrays.fill(Label, UnreachedLabel); + Arrays.fill(Age, 0); + + for (int V = 0; V < G.order(); V++ ) { + if (subset.contains(V) && matching.unmatched(V) && Search(V, G)) { + Augment(Path, alternatingTree); + nMatched += 2; + } + } + this.nMatched = nMatched; + } + + /** + * Find an augmenting path. If an augmenting path + * was found then the search must be restarted. If a blossom was detected + * the blossom is contracted and the search continues. + * + * @return an augmenting path was found + */ + boolean Search (int V, Graph G) + { + // label current vertex as even and record its age + Label[V] = EvenLabel; + Age[V] = Time++; + boolean Found = false; + + // start alt tree with current vertex + alternatingTree.addFirst(V); + + // create list of edges connected to current vertex + ArrayDeque S = new ArrayDeque(G.order()); + int W; + + // add incident edges to our stack S + for (Edge E : G.edges(V)) + { + if (!subset.contains(E.other(V)) || E.bond() == Bond.SINGLE) + continue; + + S.addLast(E); + + // peek ahead for an augmenting path and bail early if so + W = E.other(V); + if (!IsReached(W) && matching.unmatched(W)) + break; + } + + while (!S.isEmpty() && !Found) + { + Edge E = S.removeLast(); + + int X = Base(E.either()); + int Y = Base(E.other(E.either())); + if (X == Y) + continue; + if (!IsEven(X)) + { + int Z = X; + X = Y; + Y = Z; + } + + // found an augmenting path + if (!IsReached(Y) && matching.unmatched(Y)) + { + Label[Y] = OddLabel; + Tree[Y] = E; + Age[Y] = Time++; + alternatingTree.addFirst(Y); + Recover(Y); + Found = true; + break; + + // found a matched edge, need to add nbrs of mate of edge to DFS + } else if (!IsReached(Y) && matching.matched(Y)) + { + Label[Y] = OddLabel; + Tree[Y] = E; + Age[Y] = Time++; + alternatingTree.addFirst(Y); + + Edge F = G.edge(Y, matching.other(Y)); + int Z = F.other(Y); + Label[Z] = EvenLabel; + Age[Z] = Time++; + alternatingTree.addFirst(Z); + + for (Edge E2: G.edges(Z)) + { + if (E2 != F && subset.contains(E2.other(Z)) && E2.bond() != Bond.SINGLE) + { + S.addLast(E2); + + // peek ahead for an augmenting path and bail early if so + W = Other(E2, Z); + if (!IsReached(W) && matching.unmatched(W)) + break; + } + } + // found a blossom, need to shrink + } else if (IsEven(Y)) { + Shrink(E, S, G); + } + } + + if (!Found) + { + alternatingTree.clear(); + } + + return Found; + } + + private int Base(int i) { + return Base[i] == nil ? i : (Base[i] = Base(Base[i])); + } + + + /** + * Recover an augmenting path ending at vertex V by walking + * up the tree back to the root. + * + * Records a list of the unmatched edges on Path. + * + */ + private void Recover (int V) + { + do + { + Path.addFirst(Tree[V]); + int W = Other(Tree[V], V); + int B = Base(W); + Path(W, B, Path); + + V = matching.matched(B) ? matching.other(B) : nil; + } + while (V != nil); + } + + + /** + * Recursively recover the even-length piece of an alternating path + * that begins at vertex V with a matched edge and ends at base B + * of its blossom + * + * The unmatched edges on the path are added to list P, and are in arbitrary + * order. + * + */ + private void Path ( int V, int B, Deque P) { + + if (V != B) + if (IsOdd(V)) + { + Path(Shore[V], matching.other(V), P); + P.addFirst(Bridge[V]); + Path(Other(Bridge[V], Shore[V]), B, P); + } + else if (IsEven(V)) + { + int W = matching.other(V); + P.addFirst(Tree[W]); + Path(Other(Tree[W], W), B, P); + } + } + + /** + * Given an edge E between two even blossoms, shrink the implied + * cycle in the alternating tree into a superblossom + * + * Edges incident to odd vertices on the blossom are added to the stack S + * of search edges. + * + */ + private void Shrink(Edge E, ArrayDeque S, Graph G) + { + boolean Found; + int V = E.either(); + int W = E.other(V); + int baseV = Base(V); + int baseW = Base(W); + + if (Age[baseW] > Age[baseV]) + { + int temp = baseW; + baseW = baseV; + baseV = temp; + + temp = V; + V = W; + W = temp; + } + + /* + * Walk up the alternating tree from vertex V to vertex A, shrinking + * the blossoms into a superblossom. Edges incident to the odd vertices + * on the path from V to A are pushed onto stack S, to later search from. + */ + Found = false; + while (baseV != baseW) + { + Base[baseV] = baseW; + + // matched edge of B, 1 step back in DFS + Edge M = G.edge(baseV, matching.other(baseV)); + // node previous to B in DFS, the odd node + W = Other(M, baseV); + + // bridge of w is the edge after its match in the DFS + Bridge[W] = E; + + // shore of w is node one step forward in DFS, + // not necessarily the same as B which is the base of the blossom of V + Shore[W] = V; + + // tree of w is the edge leading to w from previous node in DFS + //Edge T = Tree[W]; + E = Tree[W]; + + // look for an unmatched nbr of vertex that is being shrunken into vertex and bail early if we find one + if (!Found) { + for (Edge F: G.edges(W)) { + if (F != M && F != E && subset.contains(F.other(W)) && F.bond() != Bond.SINGLE) + { + S.addLast(F); + + int Z = F.other(W); + // peek ahead and bail early if we know we're going to find an augmenting path + if (!IsReached(Z) && matching.unmatched(Z)) + { + Found = true; + break; + } + } + } + } + + Base[Base(W)] = baseW; + + // V is now the node before W in the DFS + V = E.other(W); + baseV = Base(V); + } + } + + + /** + * Augment the matching along augmenting path P, and expand + * into singleton sets all original vertices in T + * + * This assumes list P contains only the unmatched edges on the path, + * and that list T contains all vertices in all blossoms in the alternating + * tree containing the augmenting path. + * + */ + private void Augment(ArrayDeque path, ArrayDeque alternatingTree) { + Edge E = path.poll(); + while (E != null) + { + matching.match(E.either(), E.other(E.either())); + E = path.poll(); + } + + Integer V = alternatingTree.poll(); + while (V != null) + { + Label[V] = UnreachedLabel; + Base[V] = nil; + V = alternatingTree.poll(); + } + } + + /** + * Utility to maximise an existing matching of the provided graph. + * + * @param g a graph + * @param m matching on the graph, will me modified + * @param n current matching cardinality + * @param s subset of vertices to match + * @return the maximal matching on the graph + */ + static int maximise(Graph g, Matching m, int n, IntSet s) { + KececiogluMatching mm = new KececiogluMatching(g, m, n, s); + return mm.nMatched; + } +} From 6401c9d729f7865faccf122e693744c4829b0bd4 Mon Sep 17 00:00:00 2001 From: Ben Wolfe Date: Wed, 17 Jun 2020 17:15:50 -0400 Subject: [PATCH 2/5] changed Search to private method --- core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java index f9ce4b1..c669ce0 100644 --- a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java +++ b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java @@ -78,7 +78,7 @@ final class KececiogluMatching { * * @return an augmenting path was found */ - boolean Search (int V, Graph G) + private boolean Search (int V, Graph G) { // label current vertex as even and record its age Label[V] = EvenLabel; From ce9b3d83bbc56b83af4ae51d03046a3a5b43058c Mon Sep 17 00:00:00 2001 From: Ben Wolfe Date: Thu, 18 Jun 2020 09:50:54 -0400 Subject: [PATCH 3/5] moved the dfs search stack to instance variable and instantiation of the stack to constructor and renamed variables for readability --- .../uk/ac/ebi/beam/KececiogluMatching.java | 117 +++++++++--------- 1 file changed, 59 insertions(+), 58 deletions(-) diff --git a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java index c669ce0..35012a9 100644 --- a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java +++ b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java @@ -11,24 +11,25 @@ final class KececiogluMatching { private final static int UnreachedLabel = 4; private final static int nil = -1; - private final Edge[] Tree; - private final Edge[] Bridge; - private final int[] Shore; - private final int[] Base; - private final int[] Label; - private final int[] Age; - - private final ArrayDeque Path; + private final Edge[] tree; + private final Edge[] bridge; + private final int[] shore; + private final int[] base; + private final int[] label; + private final int[] age; + + private final ArrayDeque path; private final ArrayDeque alternatingTree; + private final ArrayDeque searchStack; private final int nMatched; private final Matching matching; private final IntSet subset; private int Time; - private boolean IsReached(int V) {return Label[V] != UnreachedLabel;} - private boolean IsEven(int V) {return Label[V] == EvenLabel;} - private boolean IsOdd(int V) {return Label[V] == OddLabel;} + private boolean IsReached(int V) {return label[V] != UnreachedLabel;} + private boolean IsEven(int V) {return label[V] == EvenLabel;} + private boolean IsOdd(int V) {return label[V] == OddLabel;} private int Other(Edge E, int V) {return ((E) == null) ? nil : E.other(V);} /** @@ -45,26 +46,26 @@ final class KececiogluMatching { KececiogluMatching(Graph G, Matching M, int nMatched, IntSet subset) { this.Time = 1; - this.Label = new int[G.order()]; - this.Age = new int[G.order()]; - this.Shore = new int[G.order()]; - this.Tree = new Edge[G.order()]; - this.Bridge = new Edge[G.order()]; - this.Base = new int[G.order()]; + this.label = new int[G.order()]; + this.age = new int[G.order()]; + this.shore = new int[G.order()]; + this.tree = new Edge[G.order()]; + this.bridge = new Edge[G.order()]; + this.base = new int[G.order()]; this.alternatingTree = new ArrayDeque(G.order()); - Path = new ArrayDeque(G.order()); - + this.path = new ArrayDeque(G.order()); + this.searchStack = new ArrayDeque(G.order()); this.matching = M; this.subset =subset; - Arrays.fill(Base, nil); - Arrays.fill(Label, UnreachedLabel); - Arrays.fill(Age, 0); + Arrays.fill(base, nil); + Arrays.fill(label, UnreachedLabel); + Arrays.fill(age, 0); for (int V = 0; V < G.order(); V++ ) { if (subset.contains(V) && matching.unmatched(V) && Search(V, G)) { - Augment(Path, alternatingTree); + Augment(path, alternatingTree); nMatched += 2; } } @@ -81,15 +82,15 @@ final class KececiogluMatching { private boolean Search (int V, Graph G) { // label current vertex as even and record its age - Label[V] = EvenLabel; - Age[V] = Time++; + label[V] = EvenLabel; + age[V] = Time++; boolean Found = false; // start alt tree with current vertex alternatingTree.addFirst(V); // create list of edges connected to current vertex - ArrayDeque S = new ArrayDeque(G.order()); + searchStack.clear(); int W; // add incident edges to our stack S @@ -98,7 +99,7 @@ private boolean Search (int V, Graph G) if (!subset.contains(E.other(V)) || E.bond() == Bond.SINGLE) continue; - S.addLast(E); + searchStack.addLast(E); // peek ahead for an augmenting path and bail early if so W = E.other(V); @@ -106,9 +107,9 @@ private boolean Search (int V, Graph G) break; } - while (!S.isEmpty() && !Found) + while (!searchStack.isEmpty() && !Found) { - Edge E = S.removeLast(); + Edge E = searchStack.removeLast(); int X = Base(E.either()); int Y = Base(E.other(E.either())); @@ -124,9 +125,9 @@ private boolean Search (int V, Graph G) // found an augmenting path if (!IsReached(Y) && matching.unmatched(Y)) { - Label[Y] = OddLabel; - Tree[Y] = E; - Age[Y] = Time++; + label[Y] = OddLabel; + tree[Y] = E; + age[Y] = Time++; alternatingTree.addFirst(Y); Recover(Y); Found = true; @@ -135,22 +136,22 @@ private boolean Search (int V, Graph G) // found a matched edge, need to add nbrs of mate of edge to DFS } else if (!IsReached(Y) && matching.matched(Y)) { - Label[Y] = OddLabel; - Tree[Y] = E; - Age[Y] = Time++; + label[Y] = OddLabel; + tree[Y] = E; + age[Y] = Time++; alternatingTree.addFirst(Y); Edge F = G.edge(Y, matching.other(Y)); int Z = F.other(Y); - Label[Z] = EvenLabel; - Age[Z] = Time++; + label[Z] = EvenLabel; + age[Z] = Time++; alternatingTree.addFirst(Z); for (Edge E2: G.edges(Z)) { if (E2 != F && subset.contains(E2.other(Z)) && E2.bond() != Bond.SINGLE) { - S.addLast(E2); + searchStack.addLast(E2); // peek ahead for an augmenting path and bail early if so W = Other(E2, Z); @@ -160,7 +161,7 @@ private boolean Search (int V, Graph G) } // found a blossom, need to shrink } else if (IsEven(Y)) { - Shrink(E, S, G); + Shrink(E, G); } } @@ -173,7 +174,7 @@ private boolean Search (int V, Graph G) } private int Base(int i) { - return Base[i] == nil ? i : (Base[i] = Base(Base[i])); + return base[i] == nil ? i : (base[i] = Base(base[i])); } @@ -188,10 +189,10 @@ private void Recover (int V) { do { - Path.addFirst(Tree[V]); - int W = Other(Tree[V], V); + path.addFirst(tree[V]); + int W = Other(tree[V], V); int B = Base(W); - Path(W, B, Path); + Path(W, B, path); V = matching.matched(B) ? matching.other(B) : nil; } @@ -213,15 +214,15 @@ private void Path ( int V, int B, Deque P) { if (V != B) if (IsOdd(V)) { - Path(Shore[V], matching.other(V), P); - P.addFirst(Bridge[V]); - Path(Other(Bridge[V], Shore[V]), B, P); + Path(shore[V], matching.other(V), P); + P.addFirst(bridge[V]); + Path(Other(bridge[V], shore[V]), B, P); } else if (IsEven(V)) { int W = matching.other(V); - P.addFirst(Tree[W]); - Path(Other(Tree[W], W), B, P); + P.addFirst(tree[W]); + Path(Other(tree[W], W), B, P); } } @@ -233,7 +234,7 @@ else if (IsEven(V)) * of search edges. * */ - private void Shrink(Edge E, ArrayDeque S, Graph G) + private void Shrink(Edge E, Graph G) { boolean Found; int V = E.either(); @@ -241,7 +242,7 @@ private void Shrink(Edge E, ArrayDeque S, Graph G) int baseV = Base(V); int baseW = Base(W); - if (Age[baseW] > Age[baseV]) + if (age[baseW] > age[baseV]) { int temp = baseW; baseW = baseV; @@ -260,7 +261,7 @@ private void Shrink(Edge E, ArrayDeque S, Graph G) Found = false; while (baseV != baseW) { - Base[baseV] = baseW; + base[baseV] = baseW; // matched edge of B, 1 step back in DFS Edge M = G.edge(baseV, matching.other(baseV)); @@ -268,22 +269,22 @@ private void Shrink(Edge E, ArrayDeque S, Graph G) W = Other(M, baseV); // bridge of w is the edge after its match in the DFS - Bridge[W] = E; + bridge[W] = E; // shore of w is node one step forward in DFS, // not necessarily the same as B which is the base of the blossom of V - Shore[W] = V; + shore[W] = V; // tree of w is the edge leading to w from previous node in DFS //Edge T = Tree[W]; - E = Tree[W]; + E = tree[W]; // look for an unmatched nbr of vertex that is being shrunken into vertex and bail early if we find one if (!Found) { for (Edge F: G.edges(W)) { if (F != M && F != E && subset.contains(F.other(W)) && F.bond() != Bond.SINGLE) { - S.addLast(F); + searchStack.addLast(F); int Z = F.other(W); // peek ahead and bail early if we know we're going to find an augmenting path @@ -296,7 +297,7 @@ private void Shrink(Edge E, ArrayDeque S, Graph G) } } - Base[Base(W)] = baseW; + base[Base(W)] = baseW; // V is now the node before W in the DFS V = E.other(W); @@ -325,8 +326,8 @@ private void Augment(ArrayDeque path, ArrayDeque alternatingTree) Integer V = alternatingTree.poll(); while (V != null) { - Label[V] = UnreachedLabel; - Base[V] = nil; + label[V] = UnreachedLabel; + base[V] = nil; V = alternatingTree.poll(); } } From 9d0f47ecf63ac00dca4232773bedefa9d04aeca2 Mon Sep 17 00:00:00 2001 From: Ben Wolfe Date: Thu, 18 Jun 2020 12:45:53 -0400 Subject: [PATCH 4/5] variable name and method name cleanup --- .../uk/ac/ebi/beam/KececiogluMatching.java | 279 +++++++++--------- 1 file changed, 143 insertions(+), 136 deletions(-) diff --git a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java index 35012a9..ebcec53 100644 --- a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java +++ b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java @@ -25,12 +25,12 @@ final class KececiogluMatching { private final int nMatched; private final Matching matching; private final IntSet subset; - private int Time; + private int time; - private boolean IsReached(int V) {return label[V] != UnreachedLabel;} - private boolean IsEven(int V) {return label[V] == EvenLabel;} - private boolean IsOdd(int V) {return label[V] == OddLabel;} - private int Other(Edge E, int V) {return ((E) == null) ? nil : E.other(V);} + private boolean isReached(int v) {return label[v] != UnreachedLabel;} + private boolean isEven(int v) {return label[v] == EvenLabel;} + private boolean isOdd(int v) {return label[v] == OddLabel;} + private int other(Edge e, int v) {return ((e) == null) ? nil : e.other(v);} /** * @@ -39,22 +39,22 @@ final class KececiogluMatching { * "Computing maximum-cardinality matchings in sparse general graphs." * Proceedings of the 2nd Workshop on Algorithm Engineering, 121-132, 1998. * and from their C implementation at - * https://www2.cs.arizona.edu/~kece/Research/code/matching1.sh.Z + * https://www2.cs.arizona.edu/~kece/Research/code/matching1.sh.z * * @author Ben Wolfe */ - KececiogluMatching(Graph G, Matching M, int nMatched, IntSet subset) + KececiogluMatching(Graph g, Matching M, int nMatched, IntSet subset) { - this.Time = 1; - this.label = new int[G.order()]; - this.age = new int[G.order()]; - this.shore = new int[G.order()]; - this.tree = new Edge[G.order()]; - this.bridge = new Edge[G.order()]; - this.base = new int[G.order()]; - this.alternatingTree = new ArrayDeque(G.order()); - this.path = new ArrayDeque(G.order()); - this.searchStack = new ArrayDeque(G.order()); + this.time = 1; + this.label = new int[g.order()]; + this.age = new int[g.order()]; + this.shore = new int[g.order()]; + this.tree = new Edge[g.order()]; + this.bridge = new Edge[g.order()]; + this.base = new int[g.order()]; + this.alternatingTree = new ArrayDeque(g.order()); + this.path = new ArrayDeque(g.order()); + this.searchStack = new ArrayDeque(g.order()); this.matching = M; this.subset =subset; @@ -63,9 +63,9 @@ final class KececiogluMatching { Arrays.fill(label, UnreachedLabel); Arrays.fill(age, 0); - for (int V = 0; V < G.order(); V++ ) { - if (subset.contains(V) && matching.unmatched(V) && Search(V, G)) { - Augment(path, alternatingTree); + for (int v = 0; v < g.order(); v++ ) { + if (subset.contains(v) && matching.unmatched(v) && search(v, g)) { + augment(path, alternatingTree); nMatched += 2; } } @@ -76,171 +76,175 @@ final class KececiogluMatching { * Find an augmenting path. If an augmenting path * was found then the search must be restarted. If a blossom was detected * the blossom is contracted and the search continues. - * + * @param v a free node in the graph from which to begin depth-first search + * @param g the graph in which to begin depth-first search * @return an augmenting path was found */ - private boolean Search (int V, Graph G) + private boolean search (int v, Graph g) { // label current vertex as even and record its age - label[V] = EvenLabel; - age[V] = Time++; - boolean Found = false; + label[v] = EvenLabel; + age[v] = time++; + boolean found = false; // start alt tree with current vertex - alternatingTree.addFirst(V); + alternatingTree.addFirst(v); // create list of edges connected to current vertex searchStack.clear(); - int W; + int w; // add incident edges to our stack S - for (Edge E : G.edges(V)) + for (Edge e : g.edges(v)) { - if (!subset.contains(E.other(V)) || E.bond() == Bond.SINGLE) + if (!subset.contains(e.other(v)) || e.bond() == Bond.SINGLE) continue; - searchStack.addLast(E); + searchStack.addLast(e); // peek ahead for an augmenting path and bail early if so - W = E.other(V); - if (!IsReached(W) && matching.unmatched(W)) + w = e.other(v); + if (!isReached(w) && matching.unmatched(w)) break; } - while (!searchStack.isEmpty() && !Found) + while (!searchStack.isEmpty() && !found) { - Edge E = searchStack.removeLast(); + Edge e = searchStack.removeLast(); - int X = Base(E.either()); - int Y = Base(E.other(E.either())); - if (X == Y) + int X = find(e.either()); + int y = find(e.other(e.either())); + if (X == y) continue; - if (!IsEven(X)) + if (!isEven(X)) { - int Z = X; - X = Y; - Y = Z; + int z = X; + X = y; + y = z; } // found an augmenting path - if (!IsReached(Y) && matching.unmatched(Y)) + if (!isReached(y) && matching.unmatched(y)) { - label[Y] = OddLabel; - tree[Y] = E; - age[Y] = Time++; - alternatingTree.addFirst(Y); - Recover(Y); - Found = true; + label[y] = OddLabel; + tree[y] = e; + age[y] = time++; + alternatingTree.addFirst(y); + recover(y); + found = true; break; // found a matched edge, need to add nbrs of mate of edge to DFS - } else if (!IsReached(Y) && matching.matched(Y)) + } else if (!isReached(y) && matching.matched(y)) { - label[Y] = OddLabel; - tree[Y] = E; - age[Y] = Time++; - alternatingTree.addFirst(Y); - - Edge F = G.edge(Y, matching.other(Y)); - int Z = F.other(Y); - label[Z] = EvenLabel; - age[Z] = Time++; - alternatingTree.addFirst(Z); - - for (Edge E2: G.edges(Z)) + label[y] = OddLabel; + tree[y] = e; + age[y] = time++; + alternatingTree.addFirst(y); + + Edge f = g.edge(y, matching.other(y)); + int z = f.other(y); + label[z] = EvenLabel; + age[z] = time++; + alternatingTree.addFirst(z); + + for (Edge e2: g.edges(z)) { - if (E2 != F && subset.contains(E2.other(Z)) && E2.bond() != Bond.SINGLE) + if (e2 != f && subset.contains(e2.other(z)) && e2.bond() != Bond.SINGLE) { - searchStack.addLast(E2); + searchStack.addLast(e2); // peek ahead for an augmenting path and bail early if so - W = Other(E2, Z); - if (!IsReached(W) && matching.unmatched(W)) + w = other(e2, z); + if (!isReached(w) && matching.unmatched(w)) break; } } // found a blossom, need to shrink - } else if (IsEven(Y)) { - Shrink(E, G); + } else if (isEven(y)) { + shrink(e, g); } } - if (!Found) + if (!found) { alternatingTree.clear(); } - return Found; + return found; } - private int Base(int i) { - return base[i] == nil ? i : (base[i] = Base(base[i])); + private int find(int i) { + return base[i] == nil ? i : (base[i] = find(base[i])); } /** - * Recover an augmenting path ending at vertex V by walking + * Recover an augmenting path ending at vertex v by walking * up the tree back to the root. * * Records a list of the unmatched edges on Path. - * + * @param v the free node found indicating the discovery of an augmenting path */ - private void Recover (int V) + private void recover (int v) { do { - path.addFirst(tree[V]); - int W = Other(tree[V], V); - int B = Base(W); - Path(W, B, path); + path.addFirst(tree[v]); + int w = other(tree[v], v); + int b = find(w); + path(w, b, path); - V = matching.matched(B) ? matching.other(B) : nil; + v = matching.matched(b) ? matching.other(b) : nil; } - while (V != nil); + while (v != nil); } /** * Recursively recover the even-length piece of an alternating path - * that begins at vertex V with a matched edge and ends at base B + * that begins at vertex v with a matched edge and ends at base b * of its blossom * * The unmatched edges on the path are added to list P, and are in arbitrary * order. - * + *@param v a node in the graph from which to walk to the root + *@param b the root to walk back to + *@param P a stack on which to place the unmatched edges along the walk */ - private void Path ( int V, int B, Deque P) { + private void path ( int v, int b, Deque P) { - if (V != B) - if (IsOdd(V)) + if (v != b) + if (isOdd(v)) { - Path(shore[V], matching.other(V), P); - P.addFirst(bridge[V]); - Path(Other(bridge[V], shore[V]), B, P); + path(shore[v], matching.other(v), P); + P.addFirst(bridge[v]); + path(other(bridge[v], shore[v]), b, P); } - else if (IsEven(V)) + else if (isEven(v)) { - int W = matching.other(V); - P.addFirst(tree[W]); - Path(Other(tree[W], W), B, P); + int w = matching.other(v); + P.addFirst(tree[w]); + path(other(tree[w], w), b, P); } } /** - * Given an edge E between two even blossoms, shrink the implied + * Given an edge e between two even blossoms, shrink the implied * cycle in the alternating tree into a superblossom * * Edges incident to odd vertices on the blossom are added to the stack S * of search edges. - * + * @param e a blossom-closing edge + * @param g the graph containing the discovered blossom */ - private void Shrink(Edge E, Graph G) + private void shrink(Edge e, Graph g) { - boolean Found; - int V = E.either(); - int W = E.other(V); - int baseV = Base(V); - int baseW = Base(W); + boolean found; + int v = e.either(); + int w = e.other(v); + int baseV = find(v); + int baseW = find(w); if (age[baseW] > age[baseV]) { @@ -248,60 +252,60 @@ private void Shrink(Edge E, Graph G) baseW = baseV; baseV = temp; - temp = V; - V = W; - W = temp; + temp = v; + v = w; + w = temp; } /* - * Walk up the alternating tree from vertex V to vertex A, shrinking + * Walk up the alternating tree from vertex v to vertex a, shrinking * the blossoms into a superblossom. Edges incident to the odd vertices - * on the path from V to A are pushed onto stack S, to later search from. + * on the path from v to a are pushed onto stack S, to later search from. */ - Found = false; + found = false; while (baseV != baseW) { base[baseV] = baseW; - // matched edge of B, 1 step back in DFS - Edge M = G.edge(baseV, matching.other(baseV)); - // node previous to B in DFS, the odd node - W = Other(M, baseV); + // matched edge of b, 1 step back in DFS + Edge m = g.edge(baseV, matching.other(baseV)); + // node previous to b in DFS, the odd node + w = other(m, baseV); // bridge of w is the edge after its match in the DFS - bridge[W] = E; + bridge[w] = e; // shore of w is node one step forward in DFS, - // not necessarily the same as B which is the base of the blossom of V - shore[W] = V; + // not necessarily the same as b which is the base of the blossom of v + shore[w] = v; // tree of w is the edge leading to w from previous node in DFS - //Edge T = Tree[W]; - E = tree[W]; + //Edge T = Tree[w]; + e = tree[w]; // look for an unmatched nbr of vertex that is being shrunken into vertex and bail early if we find one - if (!Found) { - for (Edge F: G.edges(W)) { - if (F != M && F != E && subset.contains(F.other(W)) && F.bond() != Bond.SINGLE) + if (!found) { + for (Edge f: g.edges(w)) { + if (f != m && f != e && subset.contains(f.other(w)) && f.bond() != Bond.SINGLE) { - searchStack.addLast(F); + searchStack.addLast(f); - int Z = F.other(W); + int z = f.other(w); // peek ahead and bail early if we know we're going to find an augmenting path - if (!IsReached(Z) && matching.unmatched(Z)) + if (!isReached(z) && matching.unmatched(z)) { - Found = true; + found = true; break; } } } } - base[Base(W)] = baseW; + base[find(w)] = baseW; - // V is now the node before W in the DFS - V = E.other(W); - baseV = Base(V); + // v is now the node before w in the DFS + v = e.other(w); + baseV = find(v); } } @@ -314,21 +318,24 @@ private void Shrink(Edge E, Graph G) * and that list T contains all vertices in all blossoms in the alternating * tree containing the augmenting path. * + * @param path a stack of unmatched edges along the augmenting path + * @param alternatingTree a stack of nodes in all the blossoms in the alternating tree + * containing the augmenting path */ - private void Augment(ArrayDeque path, ArrayDeque alternatingTree) { - Edge E = path.poll(); - while (E != null) + private void augment(ArrayDeque path, ArrayDeque alternatingTree) { + Edge e = path.poll(); + while (e != null) { - matching.match(E.either(), E.other(E.either())); - E = path.poll(); + matching.match(e.either(), e.other(e.either())); + e = path.poll(); } - Integer V = alternatingTree.poll(); - while (V != null) + Integer v = alternatingTree.poll(); + while (v != null) { - label[V] = UnreachedLabel; - base[V] = nil; - V = alternatingTree.poll(); + label[v] = UnreachedLabel; + base[v] = nil; + v = alternatingTree.poll(); } } From e61648f3c7b2244120fee1e5d3c2ca7d4e39a41b Mon Sep 17 00:00:00 2001 From: Ben Wolfe Date: Thu, 18 Jun 2020 17:06:30 -0400 Subject: [PATCH 5/5] Junit tests and utility methods to meet the test requirement of the MaximumMatching tests --- .../uk/ac/ebi/beam/KececiogluMatching.java | 23 +++ .../ac/ebi/beam/KececiogluMatchingTest.java | 160 ++++++++++++++++++ 2 files changed, 183 insertions(+) create mode 100644 core/src/test/java/uk/ac/ebi/beam/KececiogluMatchingTest.java diff --git a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java index ebcec53..0fa6f4a 100644 --- a/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java +++ b/core/src/main/java/uk/ac/ebi/beam/KececiogluMatching.java @@ -352,4 +352,27 @@ static int maximise(Graph g, Matching m, int n, IntSet s) { KececiogluMatching mm = new KececiogluMatching(g, m, n, s); return mm.nMatched; } + + /** + * Utility to maximise an existing matching of the provided graph. + * + * @param g a graph + * @param m matching on the graph + * @return the maximal matching on the graph + */ + static int maximise(Graph g, Matching m, int n) { + return maximise(g, m, n, IntSet.universe()); + } + + /** + * Utility to get the maximal matching of the specified graph. + * + * @param g a graph + * @return the maximal matching on the graph + */ + static Matching maximal(Graph g) { + Matching m = Matching.empty(g); + maximise(g, m, 0); + return m; + } } diff --git a/core/src/test/java/uk/ac/ebi/beam/KececiogluMatchingTest.java b/core/src/test/java/uk/ac/ebi/beam/KececiogluMatchingTest.java new file mode 100644 index 0000000..7247675 --- /dev/null +++ b/core/src/test/java/uk/ac/ebi/beam/KececiogluMatchingTest.java @@ -0,0 +1,160 @@ +package uk.ac.ebi.beam; + +import org.hamcrest.collection.IsIterableWithSize; +import org.junit.Test; + +import static org.hamcrest.CoreMatchers.hasItems; +import static org.junit.Assert.assertThat; + +/** @author Ben Wolfe */ +public class KececiogluMatchingTest { + + + + + /** Contrived example to test blossoming. */ + @Test public void blossom() throws Exception { + + Graph g = Graph.fromSmiles("CCCCCC1CCCC1CC"); + Matching m = Matching.empty(g); + + // initial matching from double-bonds (size = 5) + m.match(1, 2); + m.match(3, 4); + m.match(5, 6); + m.match(7, 8); + m.match(9, 10); + + KececiogluMatching.maximise(g, m, 10); + + // once maximised the matching has been augmented such that there + // are now six disjoint edges (only possibly by contracting blossom) + assertThat(m.matches(), + IsIterableWithSize.iterableWithSize(6)); + assertThat(m.matches(), + hasItems(Tuple.of(0, 1), + Tuple.of(2, 3), + Tuple.of(4, 5), + Tuple.of(6, 7), + Tuple.of(8, 9), + Tuple.of(10, 11))); + } + + @Test public void simple_maximal() throws Exception { + Graph g = Graph.fromSmiles("cccc"); + Matching m = MaximumMatching.maximal(g); + assertThat(m.matches(), + hasItems(Tuple.of(0, 1), + Tuple.of(2, 3))); + } + + @Test public void simple_augment() throws Exception { + Graph g = Graph.fromSmiles("cccc"); + Matching m = Matching.empty(g); + m.match(1, 2); + MaximumMatching.maximise(g, m, 2); + assertThat(m.matches(), + hasItems(Tuple.of(0, 1), + Tuple.of(2, 3))); + } + + @Test public void simple_augment_subset() throws Exception { + Graph g = Graph.fromSmiles("cccc"); + Matching m = Matching.empty(g); + m.match(1, 2); + // no vertex '3' matching can not be improved + MaximumMatching.maximise(g, m, 2, IntSet.allOf(0, 1, 2)); + assertThat(m.matches(), + hasItems(Tuple.of(1, 2))); + } + + @Test public void furan() throws Exception { + Graph g = Graph.fromSmiles("o1cccc1"); + IntSet s = IntSet.allOf(1, 2, 3, 4); // exclude the oxygen + Matching m = Matching.empty(g); + MaximumMatching.maximise(g, m, 0, s); + assertThat(m.matches(), hasItems(Tuple.of(1, 2), + Tuple.of(3, 4))); + } + + @Test public void furan_augment() throws Exception { + Graph g = Graph.fromSmiles("o1cccc1"); + IntSet s = IntSet.allOf(1, 2, 3, 4); // exclude the oxygen + Matching m = Matching.empty(g); + m.match(2, 3); + MaximumMatching.maximise(g, m, 2, s); + assertThat(m.matches(), hasItems(Tuple.of(1, 2), + Tuple.of(3, 4))); + } + + @Test public void quinone() throws Exception { + Graph g = Graph.fromSmiles("oc1ccc(o)cc1"); + Matching m = MaximumMatching.maximal(g); + assertThat(m.matches(), hasItems(Tuple.of(0, 1), + Tuple.of(2, 3), + Tuple.of(4, 5), + Tuple.of(6, 7))); + } + + @Test public void quinone_subset() throws Exception { + Graph g = Graph.fromSmiles("oc1ccc(o)cc1"); + // mocks the case where the oxygen atoms are already double bonded - we + // therefore don't include those of the adjacent carbons in the vertex + // subset to be matched + Matching m = Matching.empty(g); + MaximumMatching.maximise(g, m, 0, IntSet.allOf(2, 3, 6, 7)); + assertThat(m.matches(), hasItems(Tuple.of(2, 3), + Tuple.of(6, 7))); + } + + @Test public void napthalene_augment() throws Exception { + Graph g = Graph.fromSmiles("C1C=CC2=CCC=CC2=C1"); + Matching m = Matching.empty(g); + m.match(1, 2); + m.match(3, 4); + m.match(6, 7); + m.match(8, 9); + MaximumMatching.maximise(g, m, 8); + assertThat(m.matches(), hasItems(Tuple.of(0, 1), + Tuple.of(2, 3), + Tuple.of(4, 5), + Tuple.of(6, 7), + Tuple.of(8, 9))); + } + + @Test public void azulene() throws Exception { + Graph g = Graph.fromSmiles("C1CC2CCCCCC2C1"); + Matching m = MaximumMatching.maximal(g); + assertThat(m.matches(), hasItems(Tuple.of(0, 1), + Tuple.of(2, 3), + Tuple.of(4, 5), + Tuple.of(6, 7), + Tuple.of(8, 9))); + } + + @Test public void imidazole() throws Exception { + Graph g = Graph.fromSmiles("[nH]1ccnc1"); + Matching m = Matching.empty(g); + MaximumMatching.maximise(g, + m, + 0, + IntSet.allOf(1, 2, 3, 4)); // not the 'nH' + assertThat(m.matches(), hasItems(Tuple.of(1, 2), + Tuple.of(3, 4))); + } + + @Test public void benzimidazole() throws Exception { + Graph g = Graph.fromSmiles("c1nc2ccccc2[nH]1"); + Matching m = Matching.empty(g); + MaximumMatching.maximise(g, + m, + 0, + IntSet.noneOf(8)); // not the 'nH' + assertThat(m.matches(), hasItems(Tuple.of(0, 1), + Tuple.of(2, 3), + Tuple.of(4, 5), + Tuple.of(6, 7))); + } + + +}