Now finds edge probabilities for Boltzmann distributed random edges.

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;
 }