saner recording of moments

1 files changed, 35 insertions, 27 deletions

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