From fcd9bcf98790616ca6367de6ecdb75809b56a040 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 10 Oct 2022 15:57:46 +0200
Subject: Finished implementing the Propp algorithm for edge probabilities.

---
 Makefile    |   2 +-
 uniform.cpp | 115 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 113 insertions(+), 4 deletions(-)

diff --git a/Makefile b/Makefile
index 33c8b6c..31d3891 100644
--- a/Makefile
+++ b/Makefile
@@ -4,7 +4,7 @@ BLOSSOM_DIRS := $(BLOSSOM_DIR) $(BLOSSOM_DIR)/MinCost $(BLOSSOM_DIR)/GEOM
 BLOSSOM_SOURCES := $(filter-out $(BLOSSOM_DIR)/example.cpp, $(foreach dir, $(BLOSSOM_DIRS), $(wildcard $(dir)/*.cpp)))
 BLOSSOM_OBJS := $(patsubst %.cpp, %.o, $(BLOSSOM_SOURCES))
 
-CFLAGS := -flto -fopenmp -O3 -Os -march=native -D_NDEBUG -DPERFECT_MATCHING_DOUBLE
+CFLAGS := -flto -fopenmp -O3 -Os -march=native -mtune=native -D_NDEBUG -DPERFECT_MATCHING_DOUBLE
 CXX = clang++
 LD = ld.lld
 LIBS := -lrt
diff --git a/uniform.cpp b/uniform.cpp
index 07d794e..beaf038 100644
--- a/uniform.cpp
+++ b/uniform.cpp
@@ -38,13 +38,13 @@ public:
       edges[i].index = i;
     }
     for (unsigned i = 1; i <= n; i++) {
-      faces[i - 1].reserve(pow(i, 2));
+      faces[n - i].reserve(pow(i, 2));
       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;
-        faces[i - 1].push_back(Face(
+        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],
@@ -53,6 +53,62 @@ 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();
+
+        double cellFactor = w * z + x * y;
+
+        f.edges[0].get().weights.push_back(z / cellFactor);
+        f.edges[1].get().weights.push_back(y / cellFactor);
+        f.edges[2].get().weights.push_back(x / cellFactor);
+        f.edges[3].get().weights.push_back(w / cellFactor);
+      }
+    }
+
+    for (Edge& e : edges) {
+      e.weights.pop_back();
+    }
+  }
+
+  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;
+
+        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();
+
+        std::cout << w << " " << x << " " << y << " " << z << std::endl;
+
+        double cellFactor = w * z + x * y;
+        double 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;
+        f.edges[2].get().probability = q + deficit * x * y / cellFactor;
+        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();
+          }
+        }
+      }
+    }
+  }
 };
 
 int main(int argc, char* argv[]) {
@@ -75,7 +131,7 @@ int main(int argc, char* argv[]) {
   }
 
   Rng r;
-  AztecDiamond a(n);
+//  AztecDiamond a(n);
 
 
 
@@ -93,5 +149,58 @@ int main(int argc, char* argv[]) {
   }
   */
 
+  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();
+      }
+    }
+  }
+  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 << " ";
+  }
+    std::cout << std::endl;
+
   return 0;
 }
-- 
cgit v1.2.3-70-g09d2