From 352e1be1bf05de2ba75f93b8375ac52036c8203e Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 4 Oct 2022 12:37:14 +0200
Subject: Refactored base code and added new utility to measure the order
 paremeter.

---
 Makefile       |  11 ++--
 excitation.cpp |  94 ++++++++++++++++++++++++++++++
 order.cpp      |  55 ++++++++++++++++++
 rbmp.cpp       | 179 ---------------------------------------------------------
 rbmp.hpp       | 106 ++++++++++++++++++++++++++++++++++
 5 files changed, 262 insertions(+), 183 deletions(-)
 create mode 100644 excitation.cpp
 create mode 100644 order.cpp
 delete mode 100644 rbmp.cpp
 create mode 100644 rbmp.hpp

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/excitation.cpp b/excitation.cpp
new file mode 100644
index 0000000..05cb8c6
--- /dev/null
+++ b/excitation.cpp
@@ -0,0 +1,94 @@
+#include <iostream>
+
+#include "rbmp.hpp"
+
+int main(int argc, char* argv[]) {
+  unsigned n = 100;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+  Graph G(n, r);
+
+  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, e.weight);
+  }
+
+  pm.options.verbose = false;
+  pm.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;
+  }
+
+  std::vector<bool> matching(G.edges.size());
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    matching[i] = edgeMatched(pm, G.edges[i]);
+  }
+
+
+  while (true) {
+    unsigned eFlip = r.variate<unsigned, std::uniform_int_distribution>(0, G.edges.size() - 1);
+
+    while (!matching[eFlip]) {
+      eFlip = r.variate<unsigned, std::uniform_int_distribution>(0, G.edges.size() - 1);
+    }
+
+    pm.StartUpdate();
+    pm.UpdateCost(eFlip, 1e10);
+    pm.FinishUpdate();
+
+    pm.Solve();
+
+    unsigned m = 0;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      if (!matching[i] && edgeMatched(pm, G.edges[i])) {
+        m++;
+      }
+    }
+
+    if (m > 4 * n) {
+      std::cerr << "Size of excitation: " << m << std::endl;
+      break;
+    } else {
+      pm.StartUpdate();
+      pm.UpdateCost(eFlip, -1e10);
+      pm.FinishUpdate();
+
+      pm.Solve();
+    }
+  }
+
+  for (unsigned i = 0; i < G.edges.size(); i++) {
+    if (!matching[i] && edgeMatched(pm, G.edges[i])) {
+      std::cout
+        << G.edges[i].halfedges[0].getTail().coordinate[0] << " "
+        << G.edges[i].halfedges[0].getTail().coordinate[1] << " "
+        << G.edges[i].halfedges[0].getHead().coordinate[0] << " "
+        << G.edges[i].halfedges[0].getHead().coordinate[1] << std::endl;
+    }
+  }
+
+  return 0;
+}
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.cpp b/rbmp.cpp
deleted file mode 100644
index 8ddc591..0000000
--- a/rbmp.cpp
+++ /dev/null
@@ -1,179 +0,0 @@
-#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);
-}
-
-int main(int argc, char* argv[]) {
-  unsigned n = 100;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "n:")) != -1) {
-    switch (opt) {
-    case 'n':
-      n = atoi(optarg);
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  Rng r;
-  Graph G(n, r);
-
-  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, e.weight);
-  }
-
-  pm.options.verbose = false;
-  pm.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;
-  }
-
-  std::vector<bool> matching(G.edges.size());
-
-  for (unsigned i = 0; i < G.edges.size(); i++) {
-    matching[i] = edgeMatched(pm, G.edges[i]);
-  }
-
-
-  while (true) {
-    unsigned eFlip = r.variate<unsigned, std::uniform_int_distribution>(0, G.edges.size() - 1);
-
-    while (!matching[eFlip]) {
-      eFlip = r.variate<unsigned, std::uniform_int_distribution>(0, G.edges.size() - 1);
-    }
-
-    pm.StartUpdate();
-    pm.UpdateCost(eFlip, 1e10);
-    pm.FinishUpdate();
-
-    pm.Solve();
-
-    unsigned m = 0;
-    for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!matching[i] && edgeMatched(pm, G.edges[i])) {
-        m++;
-      }
-    }
-
-    if (m > 4 * n) {
-      std::cerr << "Size of excitation: " << m << std::endl;
-      break;
-    } else {
-      pm.StartUpdate();
-      pm.UpdateCost(eFlip, -1e10);
-      pm.FinishUpdate();
-
-      pm.Solve();
-    }
-  }
-
-  for (unsigned i = 0; i < G.edges.size(); i++) {
-    if (!matching[i] && edgeMatched(pm, G.edges[i])) {
-      std::cout
-        << G.edges[i].halfedges[0].getTail().coordinate[0] << " "
-        << G.edges[i].halfedges[0].getTail().coordinate[1] << " "
-        << G.edges[i].halfedges[0].getHead().coordinate[0] << " "
-        << G.edges[i].halfedges[0].getHead().coordinate[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);
+}
-- 
cgit v1.2.3-70-g09d2