4 files changed, 146 insertions, 0 deletions
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..d9c0c06
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule "randutils"]
+	path = randutils
+	url = https://gist.github.com/makokal/3bf4f1f66b9686384f75
+[submodule "pcg-cpp"]
+	path = pcg-cpp
+	url = https://github.com/imneme/pcg-cpp
diff --git a/pcg-cpp b/pcg-cpp
new file mode 160000
+Subproject 428802d1a5634f96bcd0705fab379ff0113bcf1
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85
diff --git a/rbmp.cpp b/rbmp.cpp
new file mode 100644
index 0000000..989509f
--- /dev/null
+++ b/rbmp.cpp
@@ -0,0 +1,140 @@
+
+#include <iostream>
+#include <cmath>
+#include <functional>
+#include <random>
+#include <vector>
+#include <limits>
+
+#include "randutils/randutils.hpp"
+#include "pcg-cpp/include/pcg_random.hpp"
+
+using Rng = randutils::random_generator<pcg32>;
+
+class HalfEdge;
+
+class Vertex {
+private:
+  unsigned index;
+
+public:
+  std::vector<HalfEdge*> neighbors;
+
+  void setIndex(unsigned i) {index = i;}
+  unsigned getIndex() const {return index;}
+  void addEdge(HalfEdge& e) {
+    neighbors.push_back(&e);
+  }
+};
+
+class HalfEdge {
+private:
+  unsigned index;
+  Vertex* tail;
+  Vertex* head;
+
+public:
+  double oldX;
+  double X;
+  double w;
+  void setIndex(unsigned i) {index = i;}
+  unsigned getIndex() const {return index;}
+  void setVertices(Vertex& vt, Vertex& vh) {
+    tail = &vt;
+    head = &vh;
+  }
+  Vertex* getHead() {
+    return head;
+  }
+  Vertex* getTail() {
+    return tail;
+  }
+  const Vertex* getHead() const {
+    return head;
+  }
+  const Vertex* getTail() const {
+    return tail;
+  }
+};
+
+class BipartiteGraph {
+public:
+  std::vector<Vertex> redVertices;
+  std::vector<Vertex> blueVertices;
+  std::vector<HalfEdge> redEdges;
+  std::vector<HalfEdge> blueEdges;
+
+  BipartiteGraph(unsigned n, Rng& r) : redVertices(n * (n + 1)), blueVertices(n * (n + 1)), redEdges(pow(2 * n, 2)), blueEdges(pow(2 * n, 2)) {
+    for (unsigned i = 0; i < redVertices.size(); i++) {
+      redVertices[i].setIndex(i);
+      blueVertices[i].setIndex(i);
+    }
+    for (unsigned i = 0; i < redEdges.size(); i++) {
+      redEdges[i].setIndex(i);
+      blueEdges[i].setIndex(i);
+      Vertex& blueVertex = blueVertices[(i % (2 * n)) / 2 + n * (((i + 2*n) / 4) / n)];
+      Vertex& redVertex = redVertices[(1+(i % (2 * n))) / 2 + (n + 1) * ((i / 4) / (n))];
+      redEdges[i].setVertices(redVertex, blueVertex);
+      blueEdges[i].setVertices(blueVertex, redVertex);
+      redVertex.addEdge(redEdges[i]);
+      blueVertex.addEdge(blueEdges[i]);
+      redEdges[i].X = 0;
+      blueEdges[i].X = 0;
+      redEdges[i].oldX = 0;
+      blueEdges[i].oldX = 0;
+      double w = r.variate<double, std::exponential_distribution>(1);
+      redEdges[i].w = w;
+      blueEdges[i].w = w;
+    }
+  }
+
+  void propagateBeliefs() {
+    for (HalfEdge& e : redEdges) {
+      double Xt = std::numeric_limits<double>::infinity();
+      const Vertex* head = e.getHead();
+      for (const HalfEdge* en : head->neighbors) {
+        if (e.getTail()->getIndex() != en->getHead()->getIndex()) {
+          Xt = std::min(en->w - en->oldX, Xt);
+        }
+      }
+      e.X = Xt;
+    }
+    for (HalfEdge& e : blueEdges) {
+      double Xt = std::numeric_limits<double>::infinity();
+      const Vertex* head = e.getHead();
+      for (const HalfEdge* en : head->neighbors) {
+        if (e.getTail()->getIndex() != en->getHead()->getIndex()) {
+          Xt = std::min(en->w - en->oldX, Xt);
+        }
+      }
+      e.X = Xt;
+    }
+    for (HalfEdge& e : redEdges) {
+      e.oldX = e.X;
+    }
+    for (HalfEdge& e : blueEdges) {
+      e.oldX = e.X;
+    }
+  }
+};
+
+int main() {
+  Rng r;
+  BipartiteGraph G(10, r);
+
+  for (unsigned i = 0; i < 1000; i++) {
+    G.propagateBeliefs();
+//    std::cout << G.redEdges[40].X << " " << G.blueEdges[40].X << std::endl;
+  }
+
+  for (unsigned i = 0; i < G.blueEdges.size(); i++) {
+    std::cout << G.blueEdges[i].getTail()->getIndex() << " " << G.blueEdges[i].getHead()->getIndex() + G.blueVertices.size() << std::endl;
+  }
+  for (unsigned i = 0; i < G.blueEdges.size(); i++) {
+    if (G.blueEdges[i].X + G.redEdges[i].X >= 0) {
+      std::cout << G.blueEdges[i].getTail()->getIndex() << " " << G.blueEdges[i].getHead()->getIndex() + G.blueVertices.size() << std::endl;
+    }
+  }
+
+}
+