From 6308773c0b6b745d49d20dc2afd6ab7ec63cb996 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 11 Oct 2022 14:25:11 +0200 Subject: Refactoring. --- Makefile | 8 +++- aztec.cpp | 18 +++++++ aztec.hpp | 142 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ excitation.cpp | 15 +----- free_energy.cpp | 2 +- order.cpp | 111 +++++++++++++++++++++++++++++-------------- rbmp.hpp | 138 ------------------------------------------------------ uniform.cpp | 2 +- 8 files changed, 246 insertions(+), 190 deletions(-) create mode 100644 aztec.cpp create mode 100644 aztec.hpp delete mode 100644 rbmp.hpp 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/aztec.hpp b/aztec.hpp new file mode 100644 index 0000000..9dc3139 --- /dev/null +++ b/aztec.hpp @@ -0,0 +1,142 @@ +#include +#include + +#include "randutils/randutils.hpp" +#include "pcg-cpp/include/pcg_random.hpp" + +#include "blossom5-v2.05.src/PerfectMatching.h" + +using Rng = randutils::random_generator; +using Real = long double; + +class AztecDiamond { +public: + using Coordinate = std::array; + + typedef struct Vertex { + unsigned index; + Coordinate coordinate; + } Vertex; + + typedef struct Edge { + Vertex* tail; + Vertex* head; + Real weight; + std::stack weights; + Real probability = 0; + } Edge; + +private: + std::tuple face(unsigned i, unsigned j) { + unsigned x0 = n - i; + unsigned x = x0 + 2 * (j % i); + unsigned y = x0 + 2 * (j / i); + + Edge& e1 = edges[2 * n * y + x]; + Edge& e2 = edges[2 * n * y + x + 1]; + Edge& e3 = edges[2 * n * (y + 1) + x]; + Edge& e4 = edges[2 * n * (y + 1) + x + 1]; + + return {e1, e2, e3, e4}; + } + +public: + unsigned n; + std::vector vertices; + std::vector edges; + + 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; + vertices[M + i].index = M + i; + vertices[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1}; + vertices[M + i].coordinate = {2 * (i % n) + 1, 2 * (i / n)}; + } + 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)]; + } + } + + void setWeights(Rng& r) { + for (Edge& e : edges) { + e.weight = r.variate(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++) { + auto [e1, e2, e3, e4] = face(i, j); + + Real w = e1.weights.top(); + Real x = e2.weights.top(); + Real y = e3.weights.top(); + Real z = e4.weights.top(); + + Real cellFactor = std::max(std::numeric_limits::min(), w * z + x * y); + + e1.weights.push(z / cellFactor); + e2.weights.push(y / cellFactor); + e3.weights.push(x / cellFactor); + e4.weights.push(w / cellFactor); + } + } + + // This process computes one extra weight per edge. + for (Edge& e : edges) { + e.weights.pop(); + } + } + + Real computeProbabilities() { // destroys *all* weights + for (Edge& e : edges) { + e.probability = 0; + } + Real logPartitionFunction = 0; + + for (unsigned i = 1; i <= n; i++) { +#pragma omp parallel for reduction(+:logPartitionFunction) + for (unsigned j = 0; j < pow(i, 2); j++) { + auto [e1, e2, e3, e4] = face(i, j); + + Real p = e1.probability; + Real q = e2.probability; + Real r = e3.probability; + Real s = e4.probability; + + Real w = e1.weights.top(); + Real x = e2.weights.top(); + Real y = e3.weights.top(); + Real z = e4.weights.top(); + + Real cellFactor = w * z + x * y; + Real deficit = 1 - p - q - r - s; + + e1.probability = s + deficit * w * z / cellFactor; + e2.probability = r + deficit * x * y / cellFactor; + e3.probability = q + deficit * x * y / cellFactor; + e4.probability = p + deficit * w * z / cellFactor; + + e1.weights.pop(); + e2.weights.pop(); + e3.weights.pop(); + e4.weights.pop(); + + logPartitionFunction += log(cellFactor); + } + } + + 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 -#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 #include -#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 -#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())) \ initializer(omp_priv = decltype(omp_orig)(omp_orig.size())) - std::vector data_x(G.vertices.size()); - std::vector 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 data_x(G.vertices.size()); + std::vector 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(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 data_x(G.vertices.size()); + std::vector 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 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/rbmp.hpp b/rbmp.hpp deleted file mode 100644 index 88f5900..0000000 --- a/rbmp.hpp +++ /dev/null @@ -1,138 +0,0 @@ -#include -#include - -#include "randutils/randutils.hpp" -#include "pcg-cpp/include/pcg_random.hpp" - -#include "blossom5-v2.05.src/PerfectMatching.h" - -using Rng = randutils::random_generator; -using Real = long double; - -class AztecDiamond { -public: - using Coordinate = std::array; - - typedef struct Vertex { - unsigned index; - Coordinate coordinate; - } Vertex; - - typedef struct Edge { - Vertex* tail; - Vertex* head; - Real weight; - std::stack weights; - Real probability = 0; - } Edge; - -private: - std::tuple face(unsigned i, unsigned j) { - unsigned x0 = n - i; - unsigned x = x0 + 2 * (j % i); - unsigned y = x0 + 2 * (j / i); - - Edge& e1 = edges[2 * n * y + x]; - Edge& e2 = edges[2 * n * y + x + 1]; - Edge& e3 = edges[2 * n * (y + 1) + x]; - Edge& e4 = edges[2 * n * (y + 1) + x + 1]; - - return {e1, e2, e3, e4}; - } - -public: - unsigned n; - std::vector vertices; - std::vector edges; - - 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; - vertices[M + i].index = M + i; - vertices[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1}; - vertices[M + i].coordinate = {2 * (i % n) + 1, 2 * (i / n)}; - } - 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)]; - } - } - - void setWeights(Rng& r) { - for (Edge& e : edges) { - e.weight = r.variate(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++) { - auto [e1, e2, e3, e4] = face(i, j); - - Real w = e1.weights.top(); - Real x = e2.weights.top(); - Real y = e3.weights.top(); - Real z = e4.weights.top(); - - Real cellFactor = std::max(std::numeric_limits::min(), w * z + x * y); - - e1.weights.push(z / cellFactor); - e2.weights.push(y / cellFactor); - e3.weights.push(x / cellFactor); - e4.weights.push(w / cellFactor); - } - } - - // This process computes one extra weight per edge. - for (Edge& e : edges) { - e.weights.pop(); - } - } - - Real computeProbabilities() { // destroys *all* weights - for (Edge& e : edges) { - e.probability = 0; - } - Real logPartitionFunction = 0; - - for (unsigned i = 1; i <= n; i++) { -#pragma omp parallel for reduction(+:logPartitionFunction) - for (unsigned j = 0; j < pow(i, 2); j++) { - auto [e1, e2, e3, e4] = face(i, j); - - Real p = e1.probability; - Real q = e2.probability; - Real r = e3.probability; - Real s = e4.probability; - - Real w = e1.weights.top(); - Real x = e2.weights.top(); - Real y = e3.weights.top(); - Real z = e4.weights.top(); - - Real cellFactor = w * z + x * y; - Real deficit = 1 - p - q - r - s; - - e1.probability = s + deficit * w * z / cellFactor; - e2.probability = r + deficit * x * y / cellFactor; - e3.probability = q + deficit * x * y / cellFactor; - e4.probability = p + deficit * w * z / cellFactor; - - e1.weights.pop(); - e2.weights.pop(); - e3.weights.pop(); - e4.weights.pop(); - - logPartitionFunction += log(cellFactor); - } - } - - return logPartitionFunction; - } -}; diff --git a/uniform.cpp b/uniform.cpp index 7e2d583..6c8184c 100644 --- a/uniform.cpp +++ b/uniform.cpp @@ -1,6 +1,6 @@ #include -#include "rbmp.hpp" +#include "aztec.hpp" int main(int argc, char* argv[]) { unsigned n = 100; -- cgit v1.2.3-70-g09d2