summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile8
-rw-r--r--aztec.cpp18
-rw-r--r--aztec.hpp (renamed from rbmp.hpp)4
-rw-r--r--excitation.cpp15
-rw-r--r--free_energy.cpp2
-rw-r--r--order.cpp111
-rw-r--r--uniform.cpp2
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;