diff options
| -rw-r--r-- | n_to_p.cpp | 87 | ||||
| -rw-r--r-- | percolation.cpp | 62 |
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; } |