From ff5e284201b2a464655b5bc4b488004905039b59 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 10 Oct 2022 16:38:32 +0200 Subject: Now finds edge probabilities for Boltzmann distributed random edges. --- uniform.cpp | 146 +++++++++++++++++++----------------------------------------- 1 file 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 #include +#include #include "randutils/randutils.hpp" #include "pcg-cpp/include/pcg_random.hpp" using Rng = randutils::random_generator; +using Real = long double; class Edge { public: - unsigned index; - std::list weights; - double probability; - - Edge() { - probability = 0; - } + std::list weights; + Real probability = 0; }; class Face { @@ -34,21 +31,17 @@ public: std::vector> 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& 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 avgProbabilities(a.edges.size()); - - /* For checking if the faces are appropriately defined. - for (std::vector& 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(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 <