From 6308773c0b6b745d49d20dc2afd6ab7ec63cb996 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 11 Oct 2022 14:25:11 +0200
Subject: Refactoring.

---
 order.cpp | 111 ++++++++++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 76 insertions(+), 35 deletions(-)

(limited to 'order.cpp')

diff --git a/order.cpp b/order.cpp
index 28b6595..64a35a3 100644
--- a/order.cpp
+++ b/order.cpp
@@ -1,14 +1,15 @@
 #include <fstream>
 
-#include "rbmp.hpp"
+#include "aztec.hpp"
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
   unsigned m = 100;
+  Real T = 1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "n:m:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:m:T:")) != -1) {
     switch (opt) {
     case 'n':
       n = atoi(optarg);
@@ -16,6 +17,9 @@ int main(int argc, char* argv[]) {
     case 'm':
       m = (unsigned)atof(optarg);
       break;
+    case 'T':
+      T = atof(optarg);
+      break;
     default:
       exit(1);
     }
@@ -28,55 +32,92 @@ int main(int argc, char* argv[]) {
                                   std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<long int>())) \
                       initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
 
-  std::vector<long int> data_x(G.vertices.size());
-  std::vector<long int> data_y(G.vertices.size());
-
-  std::string filename = "order_" + std::to_string(n) + ".dat";
+  std::string filename = "order_" + std::to_string(n) + "_" + std::to_string(T) + ".dat";
   std::ifstream input(filename);
 
-  if (input.is_open()) {
-    for (unsigned i = 0; i < G.vertices.size(); i++) {
-      input >> data_x[i];
-      input >> data_y[i];
-    }
+  if (T == 0) {
+    std::vector<long int> data_x(G.vertices.size());
+    std::vector<long int> data_y(G.vertices.size());
 
-    input.close();
-  }
+    if (input.is_open()) {
+      for (unsigned i = 0; i < G.vertices.size(); i++) {
+        input >> data_x[i];
+        input >> data_y[i];
+      }
 
+      input.close();
+    }
 
 #pragma omp parallel for reduction(vec_int_plus : data_x) reduction(vec_int_plus : data_y)
-  for (unsigned i = 0; i < m; i++) {
-    PerfectMatching pm(G.vertices.size(), G.edges.size());
+    for (unsigned i = 0; i < m; i++) {
+      G.setWeights(r);
+      PerfectMatching pm = findGroundState(G);
+
+      for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+        unsigned j = pm.GetMatch(i);
+
+        data_x[i] += G.vertices[i].coordinate[0];
+        data_y[i] += G.vertices[i].coordinate[1];
+        data_x[i] -= G.vertices[j].coordinate[0];
+        data_y[i] -= G.vertices[j].coordinate[1];
+
+        data_x[j] += G.vertices[i].coordinate[0];
+        data_y[j] += G.vertices[i].coordinate[1];
+        data_x[j] -= G.vertices[j].coordinate[0];
+        data_y[j] -= G.vertices[j].coordinate[1];
+      }
+    }
+
+    std::ofstream output(filename);
 
-    for (const AztecDiamond::Edge& e : G.edges) {
-      pm.AddEdge(e.head->index, e.tail->index, r.variate<double, std::exponential_distribution>(1));
+    for (unsigned i = 0; i < G.vertices.size(); i++) {
+      output << data_x[i] << " " << data_y[i] << std::endl;
     }
 
-    pm.options.verbose = false;
-    pm.Solve();
+    output.close();
+  } else {
+    std::vector<Real> data_x(G.vertices.size());
+    std::vector<Real> data_y(G.vertices.size());
 
-    for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
-      unsigned j = pm.GetMatch(i);
+    if (input.is_open()) {
+      for (unsigned i = 0; i < G.vertices.size(); i++) {
+        input >> data_x[i];
+        input >> data_y[i];
+      }
 
-      data_x[i] += G.vertices[i].coordinate[0];
-      data_y[i] += G.vertices[i].coordinate[1];
-      data_x[i] -= G.vertices[j].coordinate[0];
-      data_y[i] -= G.vertices[j].coordinate[1];
+      input.close();
+    }
 
-      data_x[j] += G.vertices[i].coordinate[0];
-      data_y[j] += G.vertices[i].coordinate[1];
-      data_x[j] -= G.vertices[j].coordinate[0];
-      data_y[j] -= G.vertices[j].coordinate[1];
+    std::vector<Real> avgProbabilities(G.edges.size());
+
+    for (unsigned i = 0; i < m; i++) {
+      G.setWeights(r);
+      G.computeWeights(T);
+      G.computeProbabilities();
+
+      for (unsigned j = 0; j < G.edges.size(); j++) {
+        avgProbabilities[j] += G.edges[j].probability;
+      }
     }
-  }
 
-  std::ofstream output(filename);
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      const AztecDiamond::Edge& e = G.edges[i];
+      const AztecDiamond::Vertex& vt = *e.tail;
+      const AztecDiamond::Vertex& vh =  *e.head;
+      data_x[vt.index] += avgProbabilities[i] * (vt.coordinate[0] - vh.coordinate[0]);
+      data_y[vt.index] += avgProbabilities[i] * (vt.coordinate[1] - vh.coordinate[1]);
+      data_x[vh.index] += avgProbabilities[i] * (vt.coordinate[0] - vh.coordinate[0]);
+      data_y[vh.index] += avgProbabilities[i] * (vt.coordinate[1] - vh.coordinate[1]);
+    }
 
-  for (unsigned i = 0; i < G.vertices.size(); i++) {
-    output << data_x[i] << " " << data_y[i] << std::endl;
-  }
+    std::ofstream output(filename);
 
-  output.close();
+    for (unsigned i = 0; i < G.vertices.size(); i++) {
+      output << data_x[i] << " " << data_y[i] << std::endl;
+    }
+
+    output.close();
+  }
 
   return 0;
 }
-- 
cgit v1.2.3-70-g09d2