From a06ff64534815cbf702a3847a19443612d307b80 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 30 Sep 2022 10:55:55 +0200
Subject: Changed rbmp to use blossom algorithm.

---
 rbmp.cpp | 106 ++++++++++++++++++++++++++++++---------------------------------
 1 file changed, 51 insertions(+), 55 deletions(-)

(limited to 'rbmp.cpp')

diff --git a/rbmp.cpp b/rbmp.cpp
index ffc7c44..ddd76cb 100644
--- a/rbmp.cpp
+++ b/rbmp.cpp
@@ -8,6 +8,8 @@
 #include "randutils/randutils.hpp"
 #include "pcg-cpp/include/pcg_random.hpp"
 
+#include "blossom5-v2.05.src/PerfectMatching.h"
+
 using Rng = randutils::random_generator<pcg32>;
 
 class Edge;
@@ -31,9 +33,8 @@ private:
 
 public:
   const Edge& edge;
-  double X;
 
-  HalfEdge(const Edge& e) : edge(e), X(0) {}
+  HalfEdge(const Edge& e) : edge(e) {}
   void setVertices(Vertex& vt, Vertex& vh) {
     tail = &vt;
     head = &vh;
@@ -58,12 +59,6 @@ public:
     red.addEdge(halfedges[0]);
     blue.addEdge(halfedges[1]);
   }
-  double belief() const {
-    return halfedges[0].X + halfedges[1].X;
-  }
-  bool active() const {
-    return belief() >= 0;
-  }
 };
 
 class Graph {
@@ -86,51 +81,18 @@ public:
       edges[i].weight = r.variate<double, std::exponential_distribution>(1);
     }
   }
-
-  void propagateBeliefs() {
-    for (unsigned i : {0, 1}) {
-#pragma omp parallel for
-      for (Edge& e : edges) {
-        HalfEdge& h = e.halfedges[i];
-        h.X = std::numeric_limits<double>::infinity();
-        for (const HalfEdge& hn : h.getHead().neighbors) {
-          if (h.getTail().index != hn.getHead().index) {
-            h.X = std::min(hn.edge.weight - hn.X, h.X);
-          }
-        }
-      }
-    }
-  }
-
-  double minBelief() const {
-    double minBelief = std::numeric_limits<double>::infinity();
-
-    for (const Edge& e : edges) {
-      minBelief = std::min(std::abs(e.belief()), minBelief);
-    }
-
-    return minBelief;
-  }
 };
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
-  unsigned maxSteps = 1e8;
-  double beliefThreshold = 1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "n:m:t:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:")) != -1) {
     switch (opt) {
     case 'n':
       n = atoi(optarg);
       break;
-    case 'm':
-      maxSteps = (unsigned)atof(optarg);
-      break;
-    case 't':
-      beliefThreshold = atof(optarg);
-      break;
     default:
       exit(1);
     }
@@ -139,26 +101,60 @@ int main(int argc, char* argv[]) {
   Rng r;
   Graph G(n, r);
 
-  unsigned steps = 0;
-  double minBelief = 0;
+  PerfectMatching pm(G.vertices.size(), G.edges.size());
 
-  while (minBelief < beliefThreshold && steps < maxSteps) {
-    G.propagateBeliefs();
-    if (steps % 100 == 0) {
-      std::cerr << (minBelief = G.minBelief()) << std::endl;
-    }
-    steps++;
+  for (const Edge& e : G.edges) {
+    pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
   }
 
+  pm.options.verbose = false;
+  pm.Solve();
+
+  std::reference_wrapper<Edge> eFlip = r.pick(G.edges);
+
+  while (pm.GetMatch(eFlip.get().halfedges[0].getTail().index) != eFlip.get().halfedges[0].getHead().index) {
+    eFlip = r.pick(G.edges);
+  }
+
+  eFlip.get().weight = 1e10;
+
+  PerfectMatching pm2(G.vertices.size(), G.edges.size());
+
   for (const Edge& e : G.edges) {
-    if (e.active()) {
+    pm2.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
+  }
+
+  pm2.options.verbose = false;
+  pm2.Solve();
+
+  std::cout << n << std::endl;
+
+  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+    unsigned j1 = pm.GetMatch(i);
+
+    std::cout
+      << G.vertices[i].coordinate[0] << " "
+      << G.vertices[i].coordinate[1] << " "
+      << G.vertices[j1].coordinate[0] << " "
+      << G.vertices[j1].coordinate[1] << std::endl;
+  }
+
+  unsigned m = 0;
+  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+    unsigned j1 = pm.GetMatch(i);
+    unsigned j2 = pm2.GetMatch(i);
+
+    if (j1 != j2) {
+      m++;
       std::cout
-        << e.halfedges[0].getTail().coordinate[0] << " "
-        << e.halfedges[0].getTail().coordinate[1] << " "
-        << e.halfedges[0].getHead().coordinate[0] << " "
-        << e.halfedges[0].getHead().coordinate[1] << std::endl;
+        << G.vertices[i].coordinate[0] << " "
+        << G.vertices[i].coordinate[1] << " "
+        << G.vertices[j2].coordinate[0] << " "
+        << G.vertices[j2].coordinate[1] << std::endl;
     }
   }
 
+  std::cerr << "Size of excitation: " << m << std::endl;
+
   return 0;
 }
-- 
cgit v1.2.3-70-g09d2