diff options
| -rw-r--r-- | order.cpp | 6 | ||||
| -rw-r--r-- | rbmp.hpp | 162 | ||||
| -rw-r--r-- | uniform.cpp | 113 |
3 files changed, 104 insertions, 177 deletions
@@ -24,7 +24,7 @@ int main(int argc, char* argv[]) { } Rng r; - Graph G(n, r); + AztecDiamond G(n, r); #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>())) \ @@ -50,8 +50,8 @@ int main(int argc, char* argv[]) { for (unsigned i = 0; i < m; i++) { PerfectMatching pm(G.vertices.size(), G.edges.size()); - for (const Graph::Edge& e : G.edges) { - pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, r.variate<double, std::exponential_distribution>(1)); + for (const AztecDiamond::Edge& e : G.edges) { + pm.AddEdge(e.head->index, e.tail->index, r.variate<double, std::exponential_distribution>(1)); } pm.options.verbose = false; @@ -3,6 +3,7 @@ #include <random> #include <vector> #include <limits> +#include <stack> #include "randutils/randutils.hpp" #include "pcg-cpp/include/pcg_random.hpp" @@ -10,97 +11,120 @@ #include "blossom5-v2.05.src/PerfectMatching.h" using Rng = randutils::random_generator<pcg32>; +using Real = long double; - -class Graph { +class AztecDiamond { public: -class Edge; -class HalfEdge; - -class Coordinate { - std::array<int, 2> r; - public: - void operator=(const std::array<int, 2>& rp) { - r = rp; - } - void operator+=(const Coordinate& x) { - r[0] += x.r[0]; - r[1] += x.r[1]; - } - void operator-=(const Coordinate& x) { - r[0] -= x.r[0]; - r[1] -= x.r[1]; - } - int operator[](unsigned i) const { - return r[i]; - } -}; + using Coordinate = std::array<int, 2>; -class Vertex { -public: - unsigned index; - std::vector<std::reference_wrapper<HalfEdge>> neighbors; - Coordinate coordinate; + typedef struct Vertex { + unsigned index; + Coordinate coordinate; + } Vertex; - void addEdge(HalfEdge& e) { - neighbors.push_back(e); - } -}; + typedef struct Edge { + Vertex* tail; + Vertex* head; + double weight; + std::stack<Real> weights; + Real probability = 0; + } Edge; -class HalfEdge { private: - Vertex* tail; - Vertex* head; + std::tuple<Edge&, Edge&, Edge&, Edge&> face(unsigned i, unsigned j) { + unsigned x0 = n - i; + unsigned x = x0 + 2 * (j % i); + unsigned y = x0 + 2 * (j / i); -public: - const Edge& edge; + 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]; - HalfEdge(const Edge& e) : edge(e) {} - void setVertices(Vertex& vt, Vertex& vh) { - tail = &vt; - head = &vh; + return {e1, e2, e3, e4}; } - const Vertex& getHead() const { - return *head; - } - const Vertex& getTail() const { - return *tail; - } -}; -class Edge { public: - std::array<HalfEdge, 2> halfedges; - double weight; - - Edge() : halfedges{*this, *this} {}; - void setVertices(Vertex& red, Vertex& blue) { - halfedges[0].setVertices(red, blue); - halfedges[1].setVertices(blue, red); - red.addEdge(halfedges[0]); - blue.addEdge(halfedges[1]); - } -}; + unsigned n; std::vector<Vertex> vertices; std::vector<Edge> edges; - Graph(int n, Rng& r) : vertices(2 * n * (n + 1)), edges(pow(2 * n, 2)) { + AztecDiamond(int n, Rng& r) : 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[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1}; 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++) { - Vertex& redVertex = vertices[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)]; - Vertex& blueVertex = vertices[M + (i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)]; - edges[i].setVertices(redVertex, blueVertex); + 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); } } -}; -bool edgeMatched(PerfectMatching& pm, const Graph::Edge& e) { - return e.halfedges[0].getTail().index == pm.GetMatch(e.halfedges[0].getHead().index); -} + void computeWeights() { + 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 = 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 + Real partitionFunction = 1; + + for (unsigned i = 1; i <= n; i++) { +#pragma omp parallel for + 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(); + + partitionFunction *= cellFactor; + } + } + + return partitionFunction; + } +}; diff --git a/uniform.cpp b/uniform.cpp index 2c123ff..c6f7869 100644 --- a/uniform.cpp +++ b/uniform.cpp @@ -4,101 +4,6 @@ #include "rbmp.hpp" -using Real = long double; - -class AztecDiamond { - -public: - typedef struct Edge { - std::stack<Real> weights; - Real probability = 0; - } Edge; - -private: - std::tuple<Edge&, Edge&, Edge&, Edge&> face(unsigned i, unsigned j) { - unsigned x0 = n - i; - unsigned xx = 2 * (j % i); - unsigned yy = 2 * (j / i); - - 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]; - - return {e1, e2, e3, e4}; - } - -public: - unsigned n; - std::vector<Edge> edges; - - AztecDiamond(unsigned n) : edges(pow(2 * n, 2)), n(n) {} - - void computeWeights() { - 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 = 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 - Real partitionFunction = 1; - - for (unsigned i = 1; i <= n; i++) { -#pragma omp parallel for - 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(); - - partitionFunction *= cellFactor; - } - } - - return partitionFunction; - } -}; - int main(int argc, char* argv[]) { unsigned n = 100; unsigned m = 100; @@ -125,7 +30,7 @@ int main(int argc, char* argv[]) { std::string filename = "order_" + std::to_string(n) + "_" + std::to_string(T) + ".dat"; Rng r; - AztecDiamond a(n); + AztecDiamond a(n, r); std::vector<Real> avgProbabilities(a.edges.size()); @@ -142,15 +47,13 @@ int main(int argc, char* argv[]) { } } - Graph G(n, r); - - std::vector<Real> data_x(G.vertices.size()); - std::vector<Real> data_y(G.vertices.size()); + std::vector<Real> data_x(a.vertices.size()); + std::vector<Real> data_y(a.vertices.size()); - for (unsigned i = 0; i < G.edges.size(); i++) { - const Graph::Edge& e = G.edges[i]; - const Graph::Vertex& vt = e.halfedges[0].getTail(); - const Graph::Vertex& vh = e.halfedges[0].getHead(); + for (unsigned i = 0; i < a.edges.size(); i++) { + const AztecDiamond::Edge& e = a.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]); @@ -159,7 +62,7 @@ int main(int argc, char* argv[]) { std::ofstream output(filename); - for (unsigned i = 0; i < G.vertices.size(); i++) { + for (unsigned i = 0; i < a.vertices.size(); i++) { output << data_x[i] << " " << data_y[i] << std::endl; } |