Refactored base code and added new utility to measure the order paremeter.

4 files changed, 169 insertions, 90 deletions

diff --git a/Makefile b/Makefile
index a64e9f3..3c6169b 100644
--- a/Makefile
+++ b/Makefile
@@ -4,15 +4,18 @@ BLOSSOM_DIRS := $(BLOSSOM_DIR) $(BLOSSOM_DIR)/MinCost $(BLOSSOM_DIR)/GEOM
 BLOSSOM_SOURCES := $(filter-out $(BLOSSOM_DIR)/example.cpp, $(foreach dir, $(BLOSSOM_DIRS), $(wildcard $(dir)/*.cpp)))

 BLOSSOM_OBJS := $(patsubst %.cpp, %.o, $(BLOSSOM_SOURCES))

 

-CFLAGS := -flto -O3 -D_NDEBUG -DPERFECT_MATCHING_DOUBLE

+CFLAGS := -flto -O3 -Os -march=native -D_NDEBUG -DPERFECT_MATCHING_DOUBLE

 CXX = clang++

 LD = ld.lld

 LIBS := -lrt

 

-all: rbmp

+all: excitation order

 

-rbmp: rbmp.cpp $(BLOSSOM_DIR)/blossom5.o

-	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o rbmp.cpp -o $@

+order: order.cpp $(BLOSSOM_DIR)/blossom5.o

+	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o order.cpp -o $@

+

+excitation: excitation.cpp $(BLOSSOM_DIR)/blossom5.o

+	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o excitation.cpp -o $@

 

 $(BLOSSOM_DIR)/blossom5.o: ${BLOSSOM_OBJS}

 	$(LD) -r -o $@ ${BLOSSOM_OBJS}

diff --git a/rbmp.cpp b/excitation.cpp
index 8ddc591..05cb8c6 100644
--- a/rbmp.cpp
+++ b/excitation.cpp
@@ -1,91 +1,6 @@
 #include <iostream>
-#include <cmath>
-#include <functional>
-#include <random>
-#include <vector>
-#include <limits>
 
-#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;
-class HalfEdge;
-
-class Vertex {
-public:
-  unsigned index;
-  std::vector<std::reference_wrapper<HalfEdge>> neighbors;
-  std::array<unsigned, 2> coordinate;
-
-  void addEdge(HalfEdge& e) {
-    neighbors.push_back(e);
-  }
-};
-
-class HalfEdge {
-private:
-  Vertex* tail;
-  Vertex* head;
-
-public:
-  const Edge& edge;
-
-  HalfEdge(const Edge& e) : edge(e) {}
-  void setVertices(Vertex& vt, Vertex& vh) {
-    tail = &vt;
-    head = &vh;
-  }
-  const Vertex& getHead() const {
-    return *head;
-  }
-  const Vertex& getTail() const {
-    return *tail;
-  }
-};
-
-class Edge {
-public:
-  std::array<HalfEdge, 2> halfedges;
-  double weight;
-
-  Edge() : halfedges{*this, *this} {};
-  void setVertices(Vertex& red, Vertex& blue) {
-    halfedges[0].setVertices(red, blue);
-    halfedges[1].setVertices(blue, red);
-    red.addEdge(halfedges[0]);
-    blue.addEdge(halfedges[1]);
-  }
-};
-
-class Graph {
-public:
-  std::vector<Vertex> vertices;
-  std::vector<Edge> edges;
-
-  Graph(unsigned n, Rng& r) : vertices(2 * n * (n + 1)), edges(pow(2 * n, 2)) {
-    unsigned M = vertices.size() / 2;
-    for (unsigned i = 0; i < M; i++) {
-      vertices[i].index = i;
-      vertices[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1};
-      vertices[M + i].index = M + i;
-      vertices[M + i].coordinate = {2 * (i % n) + 1, 2 * (i / n)};
-    }
-    for (unsigned i = 0; i < edges.size(); i++) {
-      Vertex& redVertex = vertices[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)];
-      Vertex& blueVertex = vertices[M + (i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)];
-      edges[i].setVertices(redVertex, blueVertex);
-      edges[i].weight = r.variate<double, std::exponential_distribution>(1);
-    }
-  }
-};
-
-bool edgeMatched(PerfectMatching& pm, const Edge& e) {
-  return e.halfedges[0].getTail().index == pm.GetMatch(e.halfedges[0].getHead().index);
-}
+#include "rbmp.hpp"
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
diff --git a/order.cpp b/order.cpp
new file mode 100644
index 0000000..3ca57af
--- /dev/null
+++ b/order.cpp
@@ -0,0 +1,55 @@
+#include <iostream>
+#include <random>
+
+#include "rbmp.hpp"
+
+int main(int argc, char* argv[]) {
+  unsigned n = 100;
+  unsigned m = 100;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:m:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    case 'm':
+      m = atoi(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+  Graph G(n, r);
+
+  std::vector<Coordinate> data(G.vertices.size() / 2);
+
+  for (unsigned i = 0; i < m; i++) {
+    PerfectMatching pm(G.vertices.size(), G.edges.size());
+
+    for (const Edge& e : G.edges) {
+      pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, r.variate<double, std::exponential_distribution>(1));
+    }
+
+    pm.options.verbose = false;
+    pm.Solve();
+
+    for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+      unsigned j = pm.GetMatch(i);
+
+      data[i] += G.vertices[i].coordinate;
+      data[i] -= G.vertices[j].coordinate;
+    }
+  }
+
+  std::cout << n << std::endl;
+
+  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+    std::cout << data[i][0] << " " << data[i][1] << std::endl;
+  }
+
+  return 0;
+}
diff --git a/rbmp.hpp b/rbmp.hpp
new file mode 100644
index 0000000..571fb05
--- /dev/null
+++ b/rbmp.hpp
@@ -0,0 +1,106 @@
+#include <cmath>
+#include <functional>
+#include <random>
+#include <vector>
+#include <limits>
+
+#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;
+class HalfEdge;
+
+class Coordinate {
+  std::array<int, 2> r;
+  public:
+  void operator=(const std::array<int, 2>& rp) {
+    r = rp;
+  }
+  void operator+=(const Coordinate& x) {
+    r[0] += x.r[0];
+    r[1] += x.r[1];
+  }
+  void operator-=(const Coordinate& x) {
+    r[0] -= x.r[0];
+    r[1] -= x.r[1];
+  }
+  int operator[](unsigned i) const {
+    return r[i];
+  }
+};
+
+class Vertex {
+public:
+  unsigned index;
+  std::vector<std::reference_wrapper<HalfEdge>> neighbors;
+  Coordinate coordinate;
+
+  void addEdge(HalfEdge& e) {
+    neighbors.push_back(e);
+  }
+};
+
+class HalfEdge {
+private:
+  Vertex* tail;
+  Vertex* head;
+
+public:
+  const Edge& edge;
+
+  HalfEdge(const Edge& e) : edge(e) {}
+  void setVertices(Vertex& vt, Vertex& vh) {
+    tail = &vt;
+    head = &vh;
+  }
+  const Vertex& getHead() const {
+    return *head;
+  }
+  const Vertex& getTail() const {
+    return *tail;
+  }
+};
+
+class Edge {
+public:
+  std::array<HalfEdge, 2> halfedges;
+  double weight;
+
+  Edge() : halfedges{*this, *this} {};
+  void setVertices(Vertex& red, Vertex& blue) {
+    halfedges[0].setVertices(red, blue);
+    halfedges[1].setVertices(blue, red);
+    red.addEdge(halfedges[0]);
+    blue.addEdge(halfedges[1]);
+  }
+};
+
+class Graph {
+public:
+  std::vector<Vertex> vertices;
+  std::vector<Edge> edges;
+
+  Graph(int n, Rng& r) : vertices(2 * n * (n + 1)), edges(pow(2 * n, 2)) {
+    unsigned M = vertices.size() / 2;
+    for (int i = 0; i < M; i++) {
+      vertices[i].index = i;
+      vertices[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1};
+      vertices[M + i].index = M + i;
+      vertices[M + i].coordinate = {2 * (i % n) + 1, 2 * (i / n)};
+    }
+    for (unsigned i = 0; i < edges.size(); i++) {
+      Vertex& redVertex = vertices[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)];
+      Vertex& blueVertex = vertices[M + (i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)];
+      edges[i].setVertices(redVertex, blueVertex);
+      edges[i].weight = r.variate<double, std::exponential_distribution>(1);
+    }
+  }
+};
+
+bool edgeMatched(PerfectMatching& pm, const Edge& e) {
+  return e.halfedges[0].getTail().index == pm.GetMatch(e.halfedges[0].getHead().index);
+}