summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--n_to_p.cpp87
-rw-r--r--percolation.cpp62
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;
}