saner recording of moments

2 files changed, 122 insertions, 27 deletions

diff --git a/n_to_p.cpp b/n_to_p.cpp
new file mode 100644
index 0000000..04c162f
--- /dev/null
+++ b/n_to_p.cpp
@@ -0,0 +1,87 @@
+
+#include <cinttypes>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <iostream>
+
+long double lfactorial(unsigned i) { return lgamma((long double)(i + 1)); }
+
+long double lbinom(unsigned N, unsigned n) {
+  return lfactorial(N) - lfactorial(n) - lfactorial(N-n);
+}
+
+long double B(unsigned N, unsigned n, long double p) {
+  return expl(lbinom(N, n) + log(p) * n + log(1-p) * (N - n));
+}
+
+int main(int argc, char *argv[]) {
+  std::vector<std::string> all_args;
+
+  if (argc > 1) {
+    all_args.assign(argv + 1, argv + argc);
+  }
+
+  for (std::string filename : all_args) {
+    std::ifstream in_file;
+    in_file.open(filename);
+
+    std::vector<std::vector<double>> dat;
+
+    unsigned line_num = 0;
+    std::string line;
+    while (std::getline(in_file, line)) {
+      dat.resize(line_num + 1);
+      std::istringstream iss(line);
+      double a;
+      while (iss >> a) {
+        dat[line_num].push_back(a);
+      }
+      line_num++;
+    }
+
+    in_file.close();
+
+    unsigned n_moms = dat.size();
+    unsigned n_ps = 10000;
+
+    std::vector<std::vector<double>> output(n_moms);
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      output[i].resize(n_ps);
+    }
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      unsigned N = dat[i].size();
+      for (unsigned j = 0; j < n_ps; j++) {
+        long double p = (1.0 / (n_ps + 1)) * (j + 1);
+        long double mom_p = pow(1 - p, N);
+        long double tot = pow(1 - p, N);
+
+        for (unsigned n = 1; n <= N; n++) {
+          mom_p += B(N, n, p) * dat[i][n - 1];
+          tot += B(N, n, p);
+        }
+
+        output[i][j] = mom_p;
+      }
+    }
+
+    std::ofstream out_file;
+    out_file.open(filename + ".ps");
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      for (unsigned j = 0; j < output[i].size(); j++) {
+        out_file << output[i][j] << " ";
+      }
+      out_file << "\n";
+    }
+
+    out_file.close();
+  }
+
+  return 0;
+}
+
diff --git a/percolation.cpp b/percolation.cpp
index 3f82cc4..2f17f40 100644
--- a/percolation.cpp
+++ b/percolation.cpp
@@ -12,14 +12,13 @@ class Tree {
 private:
   std::vector<signed> p;
   std::vector<unsigned> o;
+  unsigned nc;
 
 public:
-  unsigned largest;
-  std::vector<unsigned> c;
+  std::vector<double> m;
 
-  Tree(unsigned n) : p(n, -1), o(n, 1), c(n, 0) {
-    largest = 1;
-    c[0] = n;
+  Tree(unsigned n, unsigned n_moments) : p(n, -1), o(n, 1), m(n_moments, 1) {
+    nc = n;
   }
 
   unsigned findroot(unsigned i) {
@@ -30,22 +29,24 @@ public:
   }
 
   void join(unsigned i, unsigned j) {
-    c[o[i] - 1]--;
-    c[o[j] - 1]--;
+    for (unsigned k = 1; k <= m.size(); k++) {
+      m[k - 1] -= pow(o[i], k) / nc;
+      m[k - 1] -= pow(o[j], k) / nc;
+    }
+
+    nc--;
 
     if (o[i] < o[j]) {
       p[i] = j;
       o[j] += o[i];
-      c[o[j] - 1]++;
-      if (o[j] > largest) {
-        largest = o[j];
+      for (unsigned k = 1; k <= m.size(); k++) {
+        m[k - 1] = ((nc + 1) * m[k - 1] + pow(o[j], k)) / nc;
       }
     } else {
       p[j] = i;
       o[i] += o[j];
-      c[o[i] - 1]++;
-      if (o[i] > largest) {
-        largest = o[i];
+      for (unsigned k = 1; k <= m.size(); k++) {
+        m[k - 1] = ((nc + 1) * m[k - 1] + pow(o[i], k)) / nc;
       }
     }
   }
@@ -60,21 +61,24 @@ public:
   }
 };
 
-void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data,
+void update_distribution_file(std::string id, unsigned n, const std::vector<std::vector<double>>& data,
                               std::string model_string) {
   std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
-  std::vector<std::vector<uint64_t>> data_old(data.size());
+  unsigned n_old = 0;
+  std::vector<std::vector<double>> data_old(data.size());
 
   for (unsigned i = 0; i < data.size(); i++) {
-    data_old[i].resize(data[0].size());
+    data_old[i].resize(data[i].size());
   }
 
   if (file.is_open()) {
+    file >> n_old;
+
     for (unsigned i = 0; i < data.size(); i++) {
-      for (unsigned j = 0; j < data[0].size(); j++) {
-        uint64_t num;
+      for (unsigned j = 0; j < data[i].size(); j++) {
+        double num;
         file >> num;
         data_old[i][j] = num;
       }
@@ -85,9 +89,11 @@ void update_distribution_file(std::string id, const std::vector<std::vector<uint
 
   std::ofstream file_out(filename);
 
+  file_out << n_old + n << "\n";
+
   for (unsigned i = 0; i < data.size(); i++) {
-    for (unsigned j = 0; j < data[0].size(); j++) {
-      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+    for (unsigned j = 0; j < data[i].size(); j++) {
+      file_out << std::fixed << (data_old[i][j] * n_old + data[i][j]) / (n_old + n) << " ";
     }
     file_out << "\n";
   }
@@ -132,26 +138,28 @@ int main(int argc, char* argv[]) {
     bonds[2 * v_i + 1] = {v_i, L * (((v_i / L) + 1) % L) + v_i % L};
   }
 
-  std::vector<std::vector<uint64_t>> sn(n_bonds);
+  unsigned n_moments = 5;
+
+  std::vector<std::vector<double>> sn(n_moments);
 
-  for (std::vector<uint64_t>& sni : sn) {
-    sni.resize(n_vertices);
+  for (std::vector<double>& sni : sn) {
+    sni.resize(n_bonds);
   }
 
   for (unsigned trial = 0; trial < N; trial++) {
-    Tree t(n_vertices);
+    Tree t(n_vertices, n_moments);
     std::shuffle(bond_indices.begin(), bond_indices.end(), rng);
 
     for (unsigned p = 0; p < n_bonds; p++) {
       t.add_bond(bonds[bond_indices[p]]);
 
-      for (unsigned i = 0; i < t.largest; i++) {
-        sn[p][i] += t.c[i];
+      for (unsigned i = 0; i < n_moments; i++) {
+        sn[i][p] += t.m[i];
       }
     }
   }
 
-  update_distribution_file("sn", sn, "percolation_" + std::to_string(L) + "_");
+  update_distribution_file("sn", N, sn, "percolation_" + std::to_string(L) + "_");
 
   return 0;
 }