summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile5
-rw-r--r--free_energy.cpp53
-rw-r--r--rbmp.hpp10
3 files changed, 62 insertions, 6 deletions
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;
}
};