Optimized matching algorithm

Author: Kalashnikovni

algorithms/matching/bfs_matching.cpp

@@ -48,7 +48,8 @@ bool BFSMatching::ExitCondition::isSatisfied()
 // Algorithm -------------------------------------------------------------------
 
 BFSMatching::BFSMatching() : _M(SET_FACT.createSet()), _dsbg()
-  , _direction(Direction::kForward) {}
+  , _direction(Direction::kForward), _X(SET_FACT.createSet())
+  , _Y(SET_FACT.createSet()) {}
 
 void BFSMatching::swapEdgesDirection(const Set& E)
 {
@@ -62,6 +63,14 @@ void BFSMatching::swapEdgesDirection(const Set& E)
   _dsbg = DirectedSBG{_dsbg.V(), _dsbg.Vmap(), mapB, mapD, _dsbg.Emap()};
 }
 
+void BFSMatching::swapDirection(const Set& E)
+{
+  swapEdgesDirection(E);
+  Set temp_X = _X;
+  _X = _Y;
+  _Y = temp_X;
+}
+
 PWMap BFSMatching::partitionSubsetEdges() const
 {
   PWMap result = PWMAP_FACT.createPWMap();
@@ -90,84 +99,51 @@ PWMap BFSMatching::partitionSubsetEdges() const
   return result;
 }
 
-Set BFSMatching::edgesInPaths(const PWMap& smap, const Set& E) const
-{
-  PWMap mapB = _dsbg.mapB().restrict(E);
-  PWMap mapD = _dsbg.mapD().restrict(E);
-
-  // Vertices that are successors of other vertices in a path
-  Set not_fixed = smap.domain().difference(smap.fixedPoints());
-  Set succs = smap.restrict(not_fixed).image();
-  // Edges whose endings are successors 
-  Set ending_edges = mapD.preImage(succs);
-  // Map from a 'successor' edge to its start
-  PWMap auxB = mapB.restrict(ending_edges);
-  // Map from edge to the successor of its start
-  PWMap map_succs = smap.composition(auxB);
- 
-  return map_succs.equalImage(mapD);
-}
-
-Set BFSMatching::directedStep(const Set& E, const Set& right_vertices)
+Set BFSMatching::directedStep(const Set& E)
 {
   PWMap mapB = _dsbg.mapB().restrict(E);
   PWMap mapD = _dsbg.mapD().restrict(E);
   PWMap Emap = partitionSubsetEdges().restrict(E);
 
   // Calculate unmatched vertices in the side determined by the current
   // direction of edges
-  Set forward_vertices = _dsbg.V().difference(right_vertices);
   Set matched_forward_vertices = mapB.image(_M);
-  if (_direction == Direction::kBackward) {
-    forward_vertices = right_vertices;
-  }
   Set unmatched_forward_vertices
-    = forward_vertices.difference(matched_forward_vertices);
+    = _X.difference(matched_forward_vertices);
 
   // Detect paths leading to unmatched_forward_vertices
   DirectedSBG restricted_dsbg{_dsbg.V(), _dsbg.Vmap(), mapB, mapD, Emap};
   BFSPaths paths;
-  PWMap smap = paths.calculate(restricted_dsbg, unmatched_forward_vertices);
 
-  // Keep edges
-  Set paths_edges = edgesInPaths(smap, E);
-  PWMap rmap = smap.mapInf();
-  Set reach_unmatched = rmap.preImage(unmatched_forward_vertices);
-  paths_edges = paths_edges.intersection(mapD.preImage(reach_unmatched));
-
-  Util::DEBUG_LOG << "paths_edges in " << _direction << " direction: "
-    << paths_edges << "\n";
-
-  return paths_edges;
+  return paths.calculate(restricted_dsbg, unmatched_forward_vertices);
 }
 
-BFSMatching::ExitCondition BFSMatching::step(const Set& right_vertices)
+BFSMatching::ExitCondition BFSMatching::step()
 {
   Set E = _dsbg.E();
 
   // Forward direction
-  Set paths_edgesD = directedStep(E, right_vertices);
+  Set P = directedStep(E);
 
   // Backward direction
-  swapEdgesDirection(E);
+  swapDirection(E);
   _direction = Direction::kBackward;
-  Set paths_edgesB = directedStep(paths_edgesD, right_vertices);
+  Set augmenting_edges = directedStep(P);
 
-  // Calculate augmenting paths and swap edges in these paths
-  Set augmenting_edges = paths_edgesB.intersection(paths_edgesD);
+  // Swap direction for edges in augmenting paths
   Util::DEBUG_LOG << "augmenting paths: " << augmenting_edges << "\n";
-  swapEdgesDirection(augmenting_edges);
+  swapDirection(augmenting_edges);
 
   // Swap directions in all graph
   swapEdgesDirection(E);
   _direction = Direction::kForward;
 
   // Calculate new matched edges
-  _M = _dsbg.mapD().preImage(right_vertices);
+  _M = _dsbg.mapD().preImage(_Y);
 
   // Calculate exit conditions
   Set matchedU = _dsbg.mapD().image(_M);
-  bool full_match = right_vertices.difference(matchedU).isEmpty();
+  bool full_match = _Y.difference(matchedU).isEmpty();
   bool found_paths = !augmenting_edges.isEmpty();
 
   return ExitCondition{full_match, found_paths};
@@ -204,11 +180,12 @@ MatchData BFSMatching::calculate(const BipartiteSBG& bsbg)
   Util::Internal::TimeProfiler profiler{"Total matching exec time: "};
 
   init(bsbg);
-  Set right_vertices = bsbg.Y();
+  _X = bsbg.X();
+  _Y = bsbg.Y();
 
   ExitCondition exit_cond{false, false};
   do {
-    exit_cond = step(right_vertices);
+    exit_cond = step();
   } while (!exit_cond.isSatisfied());
 
   _M.compact();

algorithms/matching/bfs_matching.hpp

@@ -94,7 +94,7 @@ class BFSMatching {
    * saturated and b) New paths weren't found. This is calculated here instead
    * of the main loop to avoid recalculation of certain values. 
    */
-  ExitCondition step(const Set& right_vertices);
+  ExitCondition step();
 
   /**
    * @brief Computes alternating paths in a certain direction.
@@ -103,29 +103,27 @@ class BFSMatching {
    * @return Edges belonging to alternating paths that reach unmatched vertices
    * in the aforementioned direction. 
    */
-  Set directedStep(const Set& E, const Set& right_vertices);
+  Set directedStep(const Set& E);
 
   /**
    * @brief Modifies the `dsbg_` member, swapping mapB and mapD for elements of
    * the domain that belong to `E`.
    */
   void swapEdgesDirection(const Set& E);
 
+  void swapDirection(const Set& E);
+
   /**
    * @brief Separates in different subset-edges matched and unmatched edges
    * belonging to the same subset-edge in `dsbg_`.
    */
   PWMap partitionSubsetEdges() const;
 
-  /**
-   * @brief Returns edges in `E` that belong to the paths described by the
-   * successor map `smap`. 
-   */
-  Set edgesInPaths(const PWMap& smap, const Set& E) const;
-
   DirectedSBG _dsbg; ///< Directed graph according to matching
   Set _M; ///< Matched edges
   Direction _direction;
+  Set _X;
+  Set _Y;
 };
 
 } // namespace detail

algorithms/matching/bfs_paths.cpp

@@ -34,60 +34,71 @@ namespace detail {
 
 BFSPaths::BFSPaths() {}
 
-PWMap BFSPaths::calculate(const DirectedSBG& dsbg, const Set& endings)
+Set BFSPaths::calculate(const DirectedSBG& dsbg, const Set& endings)
 {
   Set V = dsbg.V();
   PWMap mapB = dsbg.mapB();
+  PWMap auxB = mapB;
   PWMap mapD = dsbg.mapD();
+  PWMap auxD = mapD;
   PWMap Emap = dsbg.Emap();
 
   // Successor map to unmatched vertices
-  PWMap result = PWMAP_FACT.createPWMap(endings);
+  PWMap smap = PWMAP_FACT.createPWMap(endings);
 
   // A record of allowed edges to keep out cycle edges
-  Set allowed_edges = dsbg.E();
+  Set E = dsbg.E();
   // Ingoing edges to vertices that reach endings
   Set ingoing = mapD.preImage(endings); 
   // A record of visited set-edges
-  Set visitedE = SET_FACT.createSet();
+  Set visitedSE = SET_FACT.createSet();
   do {
     // Calculate successor for ith vertices
-    PWMap ingoingB = mapB.restrict(ingoing);
-    PWMap ingoingD = mapD.restrict(ingoing);
+    PWMap ingoingB = auxB.restrict(ingoing);
+    PWMap ingoingD = auxD.restrict(ingoing);
     PWMap ith_smap = ingoingB.minAdj(ingoingD);
 
     Util::DEBUG_LOG << "ith_smap: " << ith_smap << "\n";
 
     // Edges that lead to a successor
-    Set Eith = mapD.equalImage(ith_smap.composition(mapB));
+    Set Ei = auxD.equalImage(ith_smap.composition(auxB));
     // Visited set-edges
-    Set Erec = visitedE.intersection(Emap.image(Eith)); 
-    // Handle recursion
-    if (!Erec.isEmpty()) {
-      Set Eplus = Emap.preImage(Erec);
-      PWMap rec_smap = mapB.restrict(Eplus).minAdj(mapD.restrict(Eplus)); 
-      Util::DEBUG_LOG << "rec_smap: " << rec_smap << "\n";
-      ith_smap = ith_smap.combine(rec_smap);
+    Set repeatedSE = visitedSE.intersection(Emap.image(Ei)); 
+    if (!repeatedSE.isEmpty()) {
+      // Propose candidate successors for repetitive paths
+      Set Eplus = Emap.preImage(repeatedSE);
+      PWMap smap_plus = auxB.restrict(Eplus).minAdj(auxD.restrict(Eplus)); 
+      Util::DEBUG_LOG << "smap_plus: " << smap_plus << "\n";
+      ith_smap = ith_smap.combine(smap_plus);
+      visitedSE = visitedSE.difference(repeatedSE);
+    } else {
+      visitedSE = visitedSE.disjointCup(Emap.image(Ei));
     }
-    result = ith_smap.combine(result);
+    smap = ith_smap.combine(smap);
 
     // Take out other outgoing edges to avoid cycles
-    allowed_edges = allowed_edges.difference(mapB.preImage(result.domain()));
-    mapD = mapD.restrict(allowed_edges);
-    mapB = mapB.restrict(allowed_edges);
+    E = E.difference(auxB.preImage(smap.domain()));
+    auxD = auxD.restrict(E);
+    auxB = auxB.restrict(E);
 
     // Edges that reach vertices with a successor
-    ingoing = mapD.preImage(result.domain()).intersection(allowed_edges);
+    ingoing = auxD.preImage(smap.domain()).intersection(E);
 
-    visitedE = visitedE.cup(Emap.image(Eith));
-
-    Util::DEBUG_LOG << "Eith: " << Eith << "\n";
-    Util::DEBUG_LOG << "Erec: " << Erec << "\n";
-    Util::DEBUG_LOG << "visitedE: " << visitedE << "\n";
-    Util::DEBUG_LOG << "result: " << result << "\n\n";
+    Util::DEBUG_LOG << "Ei: " << Ei << "\n";
+    Util::DEBUG_LOG << "repeatedSE: " << repeatedSE << "\n";
+    Util::DEBUG_LOG << "visitedSE: " << visitedSE << "\n";
+    Util::DEBUG_LOG << "smap: " << smap << "\n\n";
   } while (!ingoing.isEmpty());
 
-  return result;
+  // Keep edges that reach vertices in U
+  Set P = smap.composition(mapB).equalImage(mapD);
+  PWMap rmap = smap.mapInf();
+  Set reach_unmatched = rmap.preImage(endings);
+  P = P.intersection(mapD.preImage(reach_unmatched));
+
+  Util::DEBUG_LOG << "BFS Paths P: " << P << "\n\n";
+
+  return P;
 }
 
 } // namespace detail

algorithms/matching/bfs_paths.hpp

@@ -55,7 +55,7 @@ class BFSPaths : public PathsContext<BFSPaths> {
    * than once, replicating the same path for every element of the same
    * Set-Vertex.
    */
-  PWMap calculate(const DirectedSBG& dsbg, const Set& endings);
+  Set calculate(const DirectedSBG& dsbg, const Set& endings);
 };
 
 } // namespace detail

algorithms/scc/decreasing_edges_mrv.cpp

@@ -71,8 +71,8 @@ PWMap LtEdgesMRV::repetitivePaths(const PWMap& rmap
       Vi = _smap.image(Vi);
       V = V.disjointCup(Vi);
     }
-    PWMap rec_smap = _smap.restrict(V);
-    Set E_repetition = rec_smap.composition(mapB).equalImage(mapD);
+    PWMap smap_rep = _smap.restrict(V);
+    Set E_repetition = smap_rep.composition(mapB).equalImage(mapD);
 
     Set E_plus = Emap.preImage(Emap.image(E_repetition));
     // In the presence of a cycle, if the minimum vertex belongs to the

eval/visitors/set_impl_visitor.cpp

@@ -57,7 +57,7 @@ int SetImplExprVisitor::operator()(AST::BinOp v) const
 int SetImplExprVisitor::operator()(AST::Call v) const
 {
   int impl = 2;
-  for (const AST::Expr &e : v.args()) {
+  for (const AST::Expr& e : v.args()) {
     impl = std::min(impl, boost::apply_visitor(*this, e));
   }
 
@@ -78,7 +78,7 @@ int SetImplExprVisitor::operator()(AST::MultiDimInter v) const
     return 1;
   }
 
-  for (const AST::Expr &e : v.intervals()) {
+  for (const AST::Expr& e : v.intervals()) {
     impl = std::min(impl, boost::apply_visitor(*this, e));
   }
 
@@ -88,7 +88,7 @@ int SetImplExprVisitor::operator()(AST::MultiDimInter v) const
 int SetImplExprVisitor::operator()(AST::Set v) const
 {
   int impl = 2;
-  for (const AST::Expr &e : v.pieces()) {
+  for (const AST::Expr& e : v.pieces()) {
     impl = std::min(impl, boost::apply_visitor(*this, e));
   }
 
@@ -106,7 +106,7 @@ int SetImplExprVisitor::operator()(AST::LinearExp v) const
 int SetImplExprVisitor::operator()(AST::MDLExp v) const
 {
   int impl = 2;
-  for (const AST::Expr &e : v.exps()) {
+  for (const AST::Expr& e : v.exps()) {
     impl = std::min(impl, boost::apply_visitor(*this, e));
   }
 
@@ -124,7 +124,7 @@ int SetImplExprVisitor::operator()(AST::LinearMap v) const
 int SetImplExprVisitor::operator()(AST::PWLMap v) const
 {
   int impl = 2;
-  for (const AST::Expr &e : v.maps()) {
+  for (const AST::Expr& e : v.maps()) {
     impl = std::min(impl, boost::apply_visitor(*this, e));
   }