summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--order.cpp2
-rw-r--r--rbmp.hpp8
-rw-r--r--uniform.cpp108
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<double, std::exponential_distribution>(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<pcg32>;
+
+class Graph {
+public:
class Edge;
class HalfEdge;
@@ -78,9 +81,6 @@ public:
blue.addEdge(halfedges[1]);
}
};
-
-class Graph {
-public:
std::vector<Vertex> vertices;
std::vector<Edge> 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<Real> weights;
Real probability = 0;
} Edge;
- using Face = std::array<std::reference_wrapper<Edge>, 4>;
+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);
- std::vector<Edge> edges;
- std::vector<std::vector<Face>> 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<Edge> edges;
+
+ AztecDiamond(unsigned n) : edges(pow(2 * n, 2)), n(n) {}
+
void computeWeights() {
- for (std::vector<Face>& 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<Real> 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]);