From 122c2b4aa525df69a2e95d47e9b1202e02fc4301 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 10 Oct 2022 21:39:20 +0200
Subject: Some refactoring, and example code for comparing low-T probabilities
 to zero-T.

---
 Makefile    |  4 ++--
 order.cpp   |  4 +---
 rbmp.hpp    | 22 +++++++++++++++-------
 uniform.cpp | 45 +++++++++++++++++++++++++++++++++++++--------
 4 files changed, 55 insertions(+), 20 deletions(-)

diff --git a/Makefile b/Makefile
index e7b1f6c..d9dd908 100644
--- a/Makefile
+++ b/Makefile
@@ -11,10 +11,10 @@ LIBS := -lrt
 
 all: excitation order uniform
 
-uniform: uniform.cpp $(BLOSSOM_DIR)/blossom5.o
+uniform: uniform.cpp $(BLOSSOM_DIR)/blossom5.o rbmp.hpp
 	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o uniform.cpp -o $@
 
-order: order.cpp $(BLOSSOM_DIR)/blossom5.o
+order: order.cpp $(BLOSSOM_DIR)/blossom5.o rbmp.hpp
 	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o order.cpp -o $@
 
 excitation: excitation.cpp $(BLOSSOM_DIR)/blossom5.o
diff --git a/order.cpp b/order.cpp
index 44b5fa3..28b6595 100644
--- a/order.cpp
+++ b/order.cpp
@@ -1,6 +1,4 @@
-#include <iostream>
 #include <fstream>
-#include <random>
 
 #include "rbmp.hpp"
 
@@ -24,7 +22,7 @@ int main(int argc, char* argv[]) {
   }
 
   Rng r;
-  AztecDiamond G(n, r);
+  AztecDiamond G(n);
 
 #pragma omp declare reduction(vec_int_plus : std::vector<long int> : \
                                   std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<long int>())) \
diff --git a/rbmp.hpp b/rbmp.hpp
index 15b7e97..aefa046 100644
--- a/rbmp.hpp
+++ b/rbmp.hpp
@@ -1,8 +1,4 @@
-#include <cmath>
-#include <functional>
-#include <random>
 #include <vector>
-#include <limits>
 #include <stack>
 
 #include "randutils/randutils.hpp"
@@ -49,7 +45,7 @@ public:
   std::vector<Vertex> vertices;
   std::vector<Edge> edges;
 
-  AztecDiamond(int n, Rng& r) : n(n), vertices(2 * n * (n + 1)), edges(pow(2 * n, 2)) {
+  AztecDiamond(int n) : n(n), 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;
@@ -60,11 +56,20 @@ public:
     for (unsigned i = 0; i < edges.size(); i++) {
       edges[i].tail = &vertices[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)];
       edges[i].head = &vertices[M + (i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)];
-      edges[i].weight = r.variate<double, std::exponential_distribution>(1);
     }
   }
 
-  void computeWeights() {
+  void setWeights(Rng& r) {
+    for (Edge& e : edges) {
+      e.weight = r.variate<Real, std::exponential_distribution>(1);
+    }
+  }
+
+  void computeWeights(Real T) {
+    for (Edge& e : edges) {
+      e.weights.push(exp(-e.weight / T));
+    }
+
     for (unsigned i = n; i > 0; i--) {
 #pragma omp parallel for
       for (unsigned j = 0; j < pow(i, 2); j++) {
@@ -91,6 +96,9 @@ public:
   }
 
   Real computeProbabilities() { // destroys *all* weights
+    for (Edge& e : edges) {
+      e.probability = 0;
+    }
     Real partitionFunction = 1;
 
     for (unsigned i = 1; i <= n; i++) {
diff --git a/uniform.cpp b/uniform.cpp
index c6f7869..7e2d583 100644
--- a/uniform.cpp
+++ b/uniform.cpp
@@ -1,6 +1,4 @@
-#include <functional>
 #include <fstream>
-#include <stack>
 
 #include "rbmp.hpp"
 
@@ -30,16 +28,13 @@ int main(int argc, char* argv[]) {
   std::string filename = "order_" + std::to_string(n) + "_" + std::to_string(T) + ".dat";
 
   Rng r;
-  AztecDiamond a(n, r);
+  AztecDiamond a(n);
 
   std::vector<Real> avgProbabilities(a.edges.size());
 
   for (unsigned i = 0; i < m; i++) {
-    for (AztecDiamond::Edge& e : a.edges) {
-      e.weights.push(exp(- r.variate<Real, std::exponential_distribution>(1) / T));
-      e.probability = 0;
-    }
-    a.computeWeights();
+    a.setWeights(r);
+    a.computeWeights(T);
     a.computeProbabilities();
 
     for (unsigned j = 0; j < a.edges.size(); j++) {
@@ -68,5 +63,39 @@ int main(int argc, char* argv[]) {
 
   output.close();
 
+  std::vector<long int> data_x0(a.vertices.size());
+  std::vector<long int> data_y0(a.vertices.size());
+
+  PerfectMatching pm(a.vertices.size(), a.edges.size());
+
+  for (const AztecDiamond::Edge& e : a.edges) {
+    pm.AddEdge(e.head->index, e.tail->index, e.weight);
+  }
+
+  pm.options.verbose = false;
+  pm.Solve();
+
+  for (unsigned i = 0; i < a.vertices.size() / 2; i++) {
+    unsigned j = pm.GetMatch(i);
+
+    data_x0[i] += a.vertices[i].coordinate[0];
+    data_y0[i] += a.vertices[i].coordinate[1];
+    data_x0[i] -= a.vertices[j].coordinate[0];
+    data_y0[i] -= a.vertices[j].coordinate[1];
+
+    data_x0[j] += a.vertices[i].coordinate[0];
+    data_y0[j] += a.vertices[i].coordinate[1];
+    data_x0[j] -= a.vertices[j].coordinate[0];
+    data_y0[j] -= a.vertices[j].coordinate[1];
+  }
+
+  std::ofstream output2("order_0.dat");
+
+  for (unsigned i = 0; i < a.vertices.size(); i++) {
+    output2 << data_x0[i] << " " << data_y0[i] << std::endl;
+  }
+
+  output.close();
+
   return 0;
 }
-- 
cgit v1.2.3-54-g00ecf