Refactoring.

7 files changed, 108 insertions, 52 deletions

diff --git a/Makefile b/Makefile
index 85296fb..3caf84c 100644
--- a/Makefile
+++ b/Makefile
@@ -4,6 +4,7 @@ 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))

 

+OBJS := aztec.o

 TARGETS := excitation order uniform free_energy

 

 CFLAGS := -std=c++20 -stdlib=libc++ -flto -fopenmp -O3 -Os -march=native -mtune=native -D_NDEBUG -DPERFECT_MATCHING_DOUBLE

@@ -13,8 +14,11 @@ LIBS := -lrt
 

 all: $(TARGETS)

 

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

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

+$(TARGETS): %: %.cpp $(BLOSSOM_DIR)/blossom5.o ${OBJS} aztec.hpp

+	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o ${OBJS} $@.cpp -o $@

+

+$(OBJS): aztec.cpp aztec.hpp

+	$(CXX) $(CFLAGS) $< -c -o $@

 

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

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

diff --git a/aztec.cpp b/aztec.cpp
new file mode 100644
index 0000000..12197d9
--- /dev/null
+++ b/aztec.cpp
@@ -0,0 +1,18 @@
+#include "aztec.hpp"
+
+PerfectMatching findGroundState(const AztecDiamond& a) {
+  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();
+
+  return pm;
+}
+
+bool edgeMatched(PerfectMatching& pm, const AztecDiamond::Edge& e) {
+  return e.tail->index == pm.GetMatch(e.head->index);
+}
diff --git a/rbmp.hpp b/aztec.hpp
index 88f5900..9dc3139 100644
--- a/rbmp.hpp
+++ b/aztec.hpp
@@ -136,3 +136,7 @@ public:
     return logPartitionFunction;
   }
 };
+
+bool edgeMatched(PerfectMatching& pm, const AztecDiamond::Edge& e);
+
+PerfectMatching findGroundState(const AztecDiamond& a);
diff --git a/excitation.cpp b/excitation.cpp
index be1acc1..d09410b 100644
--- a/excitation.cpp
+++ b/excitation.cpp
@@ -1,10 +1,6 @@
 #include <iostream>
 
-#include "rbmp.hpp"
-
-bool edgeMatched(PerfectMatching& pm, const AztecDiamond::Edge& e) {
-  return e.tail->index == pm.GetMatch(e.head->index);
-}
+#include "aztec.hpp"
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
@@ -25,14 +21,7 @@ int main(int argc, char* argv[]) {
   AztecDiamond G(n);
   G.setWeights(r);
 
-  PerfectMatching pm(G.vertices.size(), G.edges.size());
-
-  for (const AztecDiamond::Edge& e : G.edges) {
-    pm.AddEdge(e.head->index, e.tail->index, e.weight);
-  }
-
-  pm.options.verbose = false;
-  pm.Solve();
+  PerfectMatching pm = findGroundState(G);
 
   std::cout << n << std::endl;
 
diff --git a/free_energy.cpp b/free_energy.cpp
index 2bc2a0e..148423f 100644
--- a/free_energy.cpp
+++ b/free_energy.cpp
@@ -1,7 +1,7 @@
 #include <iostream>
 #include <iomanip>
 
-#include "rbmp.hpp"
+#include "aztec.hpp"
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
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;
 }
diff --git a/uniform.cpp b/uniform.cpp
index 7e2d583..6c8184c 100644
--- a/uniform.cpp
+++ b/uniform.cpp
@@ -1,6 +1,6 @@
 #include <fstream>
 
-#include "rbmp.hpp"
+#include "aztec.hpp"
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;