From 721260982551ec8ba51b9e7b9419031313aa0de0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 10 Oct 2022 20:11:28 +0200 Subject: Stopped precomputing the faces to save a lot of memory and waste almost no time. --- order.cpp | 2 +- rbmp.hpp | 8 ++--- uniform.cpp | 108 +++++++++++++++++++++++++++++++++--------------------------- 3 files changed, 64 insertions(+), 54 deletions(-) diff --git a/order.cpp b/order.cpp index f43a2a3..3f43377 100644 --- a/order.cpp +++ b/order.cpp @@ -50,7 +50,7 @@ int main(int argc, char* argv[]) { for (unsigned i = 0; i < m; i++) { PerfectMatching pm(G.vertices.size(), G.edges.size()); - for (const Edge& e : G.edges) { + for (const Graph::Edge& e : G.edges) { pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, r.variate(1)); } diff --git a/rbmp.hpp b/rbmp.hpp index 571fb05..02ca5ad 100644 --- a/rbmp.hpp +++ b/rbmp.hpp @@ -11,6 +11,9 @@ using Rng = randutils::random_generator; + +class Graph { +public: class Edge; class HalfEdge; @@ -78,9 +81,6 @@ public: blue.addEdge(halfedges[1]); } }; - -class Graph { -public: std::vector vertices; std::vector edges; @@ -101,6 +101,6 @@ public: } }; -bool edgeMatched(PerfectMatching& pm, const Edge& e) { +bool edgeMatched(PerfectMatching& pm, const Graph::Edge& e) { return e.halfedges[0].getTail().index == pm.GetMatch(e.halfedges[0].getHead().index); } diff --git a/uniform.cpp b/uniform.cpp index 23cb7f8..2c123ff 100644 --- a/uniform.cpp +++ b/uniform.cpp @@ -7,49 +7,50 @@ using Real = long double; class AztecDiamond { + public: typedef struct Edge { std::stack weights; Real probability = 0; } Edge; - using Face = std::array, 4>; +private: + std::tuple face(unsigned i, unsigned j) { + unsigned x0 = n - i; + unsigned xx = 2 * (j % i); + unsigned yy = 2 * (j / i); - std::vector edges; - std::vector> lattices; + Edge& e1 = edges[2 * n * (x0 + yy) + x0 + xx]; + Edge& e2 = edges[2 * n * (x0 + yy) + x0 + xx + 1]; + Edge& e3 = edges[2 * n * (x0 + yy + 1) + x0 + xx]; + Edge& e4 = edges[2 * n * (x0 + yy + 1) + x0 + xx + 1]; - AztecDiamond(unsigned n) : edges(pow(2 * n, 2)), lattices(n) { - for (unsigned i = 1; i <= n; i++) { - lattices[n - i].reserve(pow(i, 2)); - unsigned x0 = n - i; - for (unsigned j = 0; j < pow(i, 2); j++) { - unsigned x = 2 * (j % i); - unsigned y = 2 * (j / i); - lattices[n - i].push_back({ - edges[2 * n * (x0 + y) + x0 + x], - edges[2 * n * (x0 + y) + x0 + x + 1], - edges[2 * n * (x0 + y + 1) + x0 + x], - edges[2 * n * (x0 + y + 1) + x0 + x + 1] - }); - } - } + return {e1, e2, e3, e4}; } +public: + unsigned n; + std::vector edges; + + AztecDiamond(unsigned n) : edges(pow(2 * n, 2)), n(n) {} + void computeWeights() { - for (std::vector& faces : lattices) { + for (unsigned i = n; i > 0; i--) { #pragma omp parallel for - for (Face& f : faces) { - Real w = f[0].get().weights.top(); - Real x = f[1].get().weights.top(); - Real y = f[2].get().weights.top(); - Real z = f[3].get().weights.top(); + 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 = w * z + x * y; - f[0].get().weights.push(z / cellFactor); - f[1].get().weights.push(y / cellFactor); - f[2].get().weights.push(x / cellFactor); - f[3].get().weights.push(w / cellFactor); + e1.weights.push(z / cellFactor); + e2.weights.push(y / cellFactor); + e3.weights.push(x / cellFactor); + e4.weights.push(w / cellFactor); } } @@ -59,33 +60,42 @@ public: } } - void computeProbabilities() { // destroys *all* weights - for (auto it = lattices.rbegin(); it != lattices.rend(); it++) { + Real computeProbabilities() { // destroys *all* weights + Real partitionFunction = 1; + + for (unsigned i = 1; i <= n; i++) { #pragma omp parallel for - for (Face& f : *it) { - Real p = f[0].get().probability; - Real q = f[1].get().probability; - Real r = f[2].get().probability; - Real s = f[3].get().probability; + 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 = f[0].get().weights.top(); - Real x = f[1].get().weights.top(); - Real y = f[2].get().weights.top(); - Real z = f[3].get().weights.top(); + 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; - f[0].get().probability = s + deficit * w * z / cellFactor; - f[1].get().probability = r + deficit * x * y / cellFactor; - f[2].get().probability = q + deficit * x * y / cellFactor; - f[3].get().probability = p + deficit * w * z / cellFactor; + 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; - for (Edge& e : f) { - e.weights.pop(); - } + e1.weights.pop(); + e2.weights.pop(); + e3.weights.pop(); + e4.weights.pop(); + + partitionFunction *= cellFactor; } } + + return partitionFunction; } }; @@ -138,9 +148,9 @@ int main(int argc, char* argv[]) { std::vector data_y(G.vertices.size()); for (unsigned i = 0; i < G.edges.size(); i++) { - const Edge& e = G.edges[i]; - const Vertex& vt = e.halfedges[0].getTail(); - const Vertex& vh = e.halfedges[0].getHead(); + const Graph::Edge& e = G.edges[i]; + const Graph::Vertex& vt = e.halfedges[0].getTail(); + const Graph::Vertex& vh = e.halfedges[0].getHead(); 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]); -- cgit v1.2.3-70-g09d2