summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--uniform.cpp146
1 files changed, 45 insertions, 101 deletions
diff --git a/uniform.cpp b/uniform.cpp
index beaf038..0c6c125 100644
--- a/uniform.cpp
+++ b/uniform.cpp
@@ -4,21 +4,18 @@
#include <list>
#include <cmath>
+#include <random>
#include "randutils/randutils.hpp"
#include "pcg-cpp/include/pcg_random.hpp"
using Rng = randutils::random_generator<pcg32>;
+using Real = long double;
class Edge {
public:
- unsigned index;
- std::list<double> weights;
- double probability;
-
- Edge() {
- probability = 0;
- }
+ std::list<Real> weights;
+ Real probability = 0;
};
class Face {
@@ -34,21 +31,17 @@ public:
std::vector<std::vector<Face>> faces;
AztecDiamond(unsigned n) : edges(pow(2 * n, 2)), faces(n) {
- for (unsigned i = 0; i < edges.size(); i++) {
- edges[i].index = i;
- }
for (unsigned i = 1; i <= n; i++) {
faces[n - i].reserve(pow(i, 2));
+ unsigned x0 = n - i;
for (unsigned j = 0; j < pow(i, 2); j++) {
- unsigned x = j % (i);
- unsigned y = j / (i);
- unsigned x0 = n - i;
- unsigned y0 = n - i;
+ unsigned x = 2 * (j % i);
+ unsigned y = 2 * (j / i);
faces[n - i].push_back(Face(
- edges[2 * n * (y0 + 2 * y) + x0 + 2 * x],
- edges[2 * n * (y0 + 2 * y) + x0 + 2 * x + 1],
- edges[2 * n * (y0 + 2 * y + 1) + x0 + 2 * x],
- edges[2 * n * (y0 + 2 * y + 1) + x0 + 2 * x + 1]
+ 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]
));
}
}
@@ -57,12 +50,12 @@ public:
void computeWeights() {
for (std::vector<Face>& fs : faces) {
for (Face& f : fs) {
- double w = f.edges[0].get().weights.back();
- double x = f.edges[1].get().weights.back();
- double y = f.edges[2].get().weights.back();
- double z = f.edges[3].get().weights.back();
+ Real w = f.edges[0].get().weights.back();
+ Real x = f.edges[1].get().weights.back();
+ Real y = f.edges[2].get().weights.back();
+ Real z = f.edges[3].get().weights.back();
- double cellFactor = w * z + x * y;
+ Real cellFactor = w * z + x * y;
f.edges[0].get().weights.push_back(z / cellFactor);
f.edges[1].get().weights.push_back(y / cellFactor);
@@ -71,6 +64,7 @@ public:
}
}
+ // This process computes one extra weight per edge.
for (Edge& e : edges) {
e.weights.pop_back();
}
@@ -79,20 +73,18 @@ public:
void computeProbabilities() {
for (auto it = faces.rbegin(); it != faces.rend(); it++) {
for (Face& f : *it) {
- double p = f.edges[0].get().probability;
- double q = f.edges[1].get().probability;
- double r = f.edges[2].get().probability;
- double s = f.edges[3].get().probability;
+ Real p = f.edges[0].get().probability;
+ Real q = f.edges[1].get().probability;
+ Real r = f.edges[2].get().probability;
+ Real s = f.edges[3].get().probability;
- double w = f.edges[0].get().weights.back();
- double x = f.edges[1].get().weights.back();
- double y = f.edges[2].get().weights.back();
- double z = f.edges[3].get().weights.back();
+ Real w = f.edges[0].get().weights.back();
+ Real x = f.edges[1].get().weights.back();
+ Real y = f.edges[2].get().weights.back();
+ Real z = f.edges[3].get().weights.back();
- std::cout << w << " " << x << " " << y << " " << z << std::endl;
-
- double cellFactor = w * z + x * y;
- double deficit = 1 - p - q - r - s;
+ Real cellFactor = w * z + x * y;
+ Real deficit = 1 - p - q - r - s;
f.edges[0].get().probability = s + deficit * w * z / cellFactor;
f.edges[1].get().probability = r + deficit * x * y / cellFactor;
@@ -100,11 +92,7 @@ public:
f.edges[3].get().probability = p + deficit * w * z / cellFactor;
for (Edge& e : f.edges) {
- if (e.weights.empty()) {
- std::cout << "No weights left!" << std::endl;
- } else {
- e.weights.pop_back();
- }
+ e.weights.pop_back();
}
}
}
@@ -114,10 +102,11 @@ public:
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);
@@ -125,82 +114,37 @@ int main(int argc, char* argv[]) {
case 'm':
m = (unsigned)atof(optarg);
break;
+ case 'T':
+ T = atof(optarg);
+ break;
default:
exit(1);
}
}
Rng r;
-// AztecDiamond a(n);
+ AztecDiamond a(n);
+ std::vector<Real> avgProbabilities(a.edges.size());
-
- /* For checking if the faces are appropriately defined.
- for (std::vector<Face>& fs : a.faces) {
- for (Face& f : fs) {
- for (Edge& e : f.edges) {
- unsigned v1 = (1 + (e.index % (2 * n))) / 2 + (n + 1) * ((e.index / 4) / n);
- unsigned v2 = n * (n + 1) + (e.index % (2 * n)) / 2 + n * (((e.index + 2 * n) / 4) / n);
-
- std::cout << v1 << " " << v2 << " ";
- }
+ for (unsigned i = 0; i < m; i++) {
+ for (Edge& e : a.edges) {
+ e.weights = {exp(- r.variate<Real, std::exponential_distribution>(1) / T)};
+ e.probability = 0;
}
- std::cout << std::endl;
- }
- */
+ a.computeWeights();
+ a.computeProbabilities();
- AztecDiamond a(3);
- for (Edge& e : a.edges) {
- e.weights.push_back(1);
- }
-
- a.edges[1].weights.front() = 0;
- a.edges[4].weights.front() = 0;
- a.edges[6].weights.front() = 0;
- a.edges[11].weights.front() = 0;
- a.edges[24].weights.front() = 0;
- a.edges[29].weights.front() = 0;
- a.edges[31].weights.front() = 0;
- a.edges[34].weights.front() = 0;
-
- a.computeWeights();
-
- for (unsigned i = 0; i < 3; i++) {
for (unsigned j = 0; j < a.edges.size(); j++) {
- if (j % 6 == 0) {
- std::cout << std::endl;
- }
- if (!a.edges[j].weights.empty()) {
- std::cout << a.edges[j].weights.front() << " ";
- a.edges[j].weights.pop_front();
- }
+ avgProbabilities[j] += a.edges[j].probability;
}
}
- std::cout << std::endl;
- for (Edge& e : a.edges) {
- e.weights.push_back(1);
- }
- a.edges[1].weights.front() = 0;
- a.edges[4].weights.front() = 0;
- a.edges[6].weights.front() = 0;
- a.edges[11].weights.front() = 0;
- a.edges[24].weights.front() = 0;
- a.edges[29].weights.front() = 0;
- a.edges[31].weights.front() = 0;
- a.edges[34].weights.front() = 0;
-
- a.computeWeights();
- a.computeProbabilities();
-
- for (unsigned i = 0; i < a.edges.size(); i++) {
- if (i%6 == 0) {
- std::cout << std::endl;
- }
- std::cout << a.edges[i].probability << " ";
+ for (Real& x : avgProbabilities) {
+ std::cout << x << " ";
}
- std::cout << std::endl;
+ std::cout <<std::endl;
return 0;
}