From 766988cd11eda0d62689e6c4adc021a3d2d4f567 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 10 Oct 2022 22:55:24 +0200
Subject: Added free energy calculation, and fixed some small mistakes in the
 library.

---
 Makefile        |  5 ++++-
 free_energy.cpp | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 rbmp.hpp        | 10 +++++-----
 3 files changed, 62 insertions(+), 6 deletions(-)
 create mode 100644 free_energy.cpp

diff --git a/Makefile b/Makefile
index d9dd908..4de89d8 100644
--- a/Makefile
+++ b/Makefile
@@ -9,7 +9,10 @@ CXX = clang++
 LD = ld.lld
 LIBS := -lrt
 
-all: excitation order uniform
+all: excitation order uniform free_energy
+
+free_energy: free_energy.cpp $(BLOSSOM_DIR)/blossom5.o rbmp.hpp
+	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o free_energy.cpp -o $@
 
 uniform: uniform.cpp $(BLOSSOM_DIR)/blossom5.o rbmp.hpp
 	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o uniform.cpp -o $@
diff --git a/free_energy.cpp b/free_energy.cpp
new file mode 100644
index 0000000..6bfdf19
--- /dev/null
+++ b/free_energy.cpp
@@ -0,0 +1,53 @@
+#include <iostream>
+#include <iomanip>
+
+#include "rbmp.hpp"
+
+int main(int argc, char* argv[]) {
+  unsigned n = 100;
+  unsigned m = 100;
+  Real T0 = 1;
+  Real T1 = 10;
+  Real ΔT = 1;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "n:m:0:1:d:")) != -1) {
+    switch (opt) {
+    case 'n':
+      n = atoi(optarg);
+      break;
+    case 'm':
+      m = (unsigned)atof(optarg);
+      break;
+    case '0':
+      T0 = atof(optarg);
+      break;
+    case '1':
+      T1 = atof(optarg);
+      break;
+    case 'd':
+      ΔT = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+  AztecDiamond a(n);
+
+  for (Real T = T0; T <= T1; T += ΔT) {
+    Real avgFreeEnergy = 0;
+
+    for (unsigned i = 0; i < m; i++) {
+      a.setWeights(r);
+      a.computeWeights(T);
+      avgFreeEnergy += a.computeProbabilities();
+    }
+
+    std::cout << std::setprecision(20) << T << " " << - T * avgFreeEnergy / m / a.vertices.size() << std::endl;
+  }
+
+  return 0;
+}
diff --git a/rbmp.hpp b/rbmp.hpp
index aefa046..4cd6605 100644
--- a/rbmp.hpp
+++ b/rbmp.hpp
@@ -21,7 +21,7 @@ public:
   typedef struct Edge {
     Vertex* tail;
     Vertex* head;
-    double weight;
+    Real weight;
     std::stack<Real> weights;
     Real probability = 0;
   } Edge;
@@ -99,10 +99,10 @@ public:
     for (Edge& e : edges) {
       e.probability = 0;
     }
-    Real partitionFunction = 1;
+    Real logPartitionFunction = 0;
 
     for (unsigned i = 1; i <= n; i++) {
-#pragma omp parallel for
+#pragma omp parallel for reduction(+:logPartitionFunction)
       for (unsigned j = 0; j < pow(i, 2); j++) {
         auto [e1, e2, e3, e4] = face(i, j);
 
@@ -129,10 +129,10 @@ public:
         e3.weights.pop();
         e4.weights.pop();
 
-        partitionFunction *= cellFactor;
+        logPartitionFunction += log(cellFactor);
       }
     }
 
-    return partitionFunction;
+    return logPartitionFunction;
   }
 };
-- 
cgit v1.2.3-70-g09d2