summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile4
-rw-r--r--order.cpp4
-rw-r--r--rbmp.hpp22
-rw-r--r--uniform.cpp45
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;
}