diff options
| -rw-r--r-- | src/measurements.cpp | 297 | ||||
| -rw-r--r-- | src/measurements.hpp | 35 |
2 files changed, 270 insertions, 62 deletions
diff --git a/src/measurements.cpp b/src/measurements.cpp index 7d5d539..1159fd6 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -1,11 +1,17 @@ #include "measurements.hpp" +#include <iostream> + +template <class T> +bool is_shorter(std::list<T> l1, std::list<T> l2) { + return l1.size() < l2.size(); +} bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) { return true; } -void update_distribution_data(std::string id, const std::vector<unsigned int>& data, unsigned int N, unsigned int L, double beta) { +void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int L, double beta) { std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; std::ifstream file(filename); @@ -33,58 +39,43 @@ void update_distribution_data(std::string id, const std::vector<unsigned int>& d file_out.close(); } -void update_field_data(std::string id, double tot, const std::vector<double>& data, unsigned int L, double beta) { +template <class T> +void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int L, double beta) { std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; std::ifstream file(filename); uint64_t N_old = 0; - uint64_t tot_old = 0; - uint64_t tot_2_old = 0; std::vector<uint64_t> data_old(data.size(), 0); - std::vector<uint64_t> data_old_2(data.size(), 0); if (file.is_open()) { file >> N_old; - file >> tot_old; - file >> tot_2_old; for (unsigned int i = 0; i < data.size(); i++) { - uint64_t num; - file >> num; - data_old[i] = num; + file >> data_old[i]; } - for (unsigned int i = 0; i < data.size(); i++) { - uint64_t num; - file >> num; - data_old_2[i] = num; - } - file.close(); } std::ofstream file_out(filename); - file_out <<std::fixed<< N_old + 1 << " " << (uint64_t)(tot_old + tot) << " " << (uint64_t)(tot_2_old + pow(tot, 2))<< "\n"; - for (unsigned int i = 0; i < data.size(); i++) { - file_out << (uint64_t)(data_old[i] + data[i]) << " "; - } - file_out <<std::fixed<< "\n"; + file_out <<std::fixed<< N_old + N << "\n"; for (unsigned int i = 0; i < data.size(); i++) { - file_out << (uint64_t)(data_old_2[i] + pow(data[i], 2)) << " "; + file_out << data_old[i] + data[i] << " "; } file_out.close(); } -std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<double>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { +template <class T> +std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { for (unsigned int i = 0; i < pow(L, 2); i++) { - fftw_forward_in[i] = data[i]; + fftw_forward_in[i] = (double)data[i]; } fftw_execute(forward_plan); - std::vector<fftw_complex> output(pow(L, 2)); + std::vector<fftw_complex> output(L * (L / 2 + 1)); - for (unsigned int i = 0; i < pow(L, 2); i++) { + for (unsigned int i = 0; i < L * (L / 2 + 1); i++) { output[i][0] = fftw_forward_out[i][0]; output[i][1] = fftw_forward_out[i][1]; } @@ -92,26 +83,52 @@ std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<doubl return output; } -std::vector<double> correlation(unsigned int L, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { - for (unsigned int i = 0; i < pow(L, 2); i++) { +template <class T> +void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { + for (unsigned int i = 0; i < L * (L / 2 + 1); i++) { fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1]; fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0]; } fftw_execute(reverse_plan); - std::vector<double> output(pow(L / 2 + 1, 2)); - for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) { - output[i] = fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2); + data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2)); } +} - return output; +template <class T, class U> +void autocorrelation(unsigned int L, std::vector<T>& out_data, const std::vector<U>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { + for (unsigned int i = 0; i < pow(L, 2); i++) { + fftw_forward_in[i] = (double)data[i]; + } + + fftw_execute(forward_plan); + + for (unsigned int i = 0; i < L * (L / 2 + 1); i++) { + fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2); + fftw_reverse_in[i][1] = 0.0; + } + + fftw_execute(reverse_plan); + + for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) { + out_data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2)); + } } -ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), bin_counts(log2(L) + 1, 0), N(N), beta(beta) { - ad.resize(pow(L, 2), 0); - cd.resize(pow(L, 2), 0); +ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), N(N), beta(beta), + sc(pow(L, 2), 0), + sa(pow(L, 2), 0), + sd(pow(L, 2), 0), + sb(log2(L) + 1, 0), + Ccc(pow(L / 2 + 1, 2), 0), + Css(pow(L / 2 + 1, 2), 0), + Cdd(pow(L / 2 + 1, 2), 0), + Caa(pow(L / 2 + 1, 2), 0) +{ + Nc = 0; + Na = 0; // FFTW setup for correlation functions fftw_set_timelimit(1); @@ -135,31 +152,44 @@ ma::~ma() { fftw_destroy_plan(reverse_plan); fftw_cleanup(); - update_distribution_data("na", ad, N, L, beta); - update_distribution_data("nc", cd, N, L, beta); - update_distribution_data("bc", bin_counts, N, L, beta); + update_distribution_file("sa", sa, N, L, beta); + update_distribution_file("sc", sc, N, L, beta); + update_distribution_file("sd", sd, N, L, beta); + update_distribution_file("sb", sb, N, L, beta); + + update_field_file("Ccc", Ccc, Nc, L, beta); + update_field_file("Css", Css, N, L, beta); + update_field_file("Cdd", Cdd, N, L, beta); + update_field_file("Caa", Caa, Na, L, beta); } void ma::pre_fracture(const network &) { lv = 0; - as = 0; - avalanches.clear(); + avalanche.clear(); boost::remove_edge_if(trivial, G); } void ma::bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) { if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) { - ad[as]++; - as = 0; + sa[avalanche.size()]++; + Na++; + + std::vector<bool>avalanche_damage(pow(L, 2), 0); + + for (auto it = avalanche.begin(); it != avalanche.end(); it++) { + avalanche_damage[*it] = 1; + } + + autocorrelation(L, Caa, avalanche_damage, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + lv = cur.first / fabs(cur.second[i]) * net.thresholds[i]; - avalanches.push_back({}); + avalanche.clear(); } else { - as++; - avalanches.back().push_back(i); + avalanche.push_back(i); } - boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], G); + boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G); } void ma::post_fracture(network &n) { @@ -182,6 +212,144 @@ void ma::post_fracture(network &n) { } } + unsigned int vertex_in_largest = 0; + + while (component[vertex_in_largest] != max_i) { + vertex_in_largest++; + } + + Graph Gtmp(2 * pow(L / 2, 2)); + std::list<unsigned int> removed_edges; + class add_tree_edges : public boost::default_dfs_visitor { + public: + Graph& G; + std::list<unsigned int>& E; + + add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {} + + void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G); + } + + void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) { + E.push_back(g[e].index); + } + } + }; + add_tree_edges ate(Gtmp, removed_edges); + boost::depth_first_search(G, visitor(ate)); + + class find_cycle : public boost::default_dfs_visitor { + public: + unsigned int end; + std::list<unsigned int>& V; + std::list<unsigned int>& E; + struct done{}; + + find_cycle(std::list<unsigned int>& V, std::list<unsigned int>& E, unsigned int end) : V(V), E(E), end(end) {} + + void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) { + V.push_back(v); + if (v == end) { + throw done{}; + } + } + + void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) { + V.remove(v); + } + + void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + E.push_back(g[e].index); + } + + void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + E.erase(std::find(E.begin(), E.end(), g[e].index)); + } + }; + + if (removed_edges.size() == 0) { + std::pair<double, std::vector<double>> currents = n.currents(); + std::cout << currents.first << "\n"; + for (unsigned int i = 0; i < L; i++) { + for (unsigned int j = 0; j < L; j++) { + std::cout << n.fuses[i * L + j] << " "; + } + std::cout << "\n"; + } + } + + std::list<std::list<unsigned int>> cycles; + + for (auto edge : removed_edges) { + std::list<unsigned int> cycle = {edge}; + std::list<unsigned int> vcycle = {}; + find_cycle vis(vcycle, cycle, n.G.dual_edges[edge][1]); + std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp)); + try { + boost::depth_first_visit(Gtmp, n.G.dual_edges[edge][0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0])); + } catch(find_cycle::done const&) { + cycles.push_back(cycle); + } + } + + std::list<unsigned int> crack; + + if (cycles.size() > 1) { + std::list<std::valarray<uint8_t>> bool_cycles; + for (auto cycle : cycles) { + std::valarray<uint8_t> bool_cycle(n.G.edges.size()); + for (auto v : cycle) { + bool_cycle[v] = 1; + } + bool_cycles.push_back(bool_cycle); + } + + for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) { + for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) { + std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2); + std::list<unsigned int> new_cycle; + unsigned int pos = 0; + for (uint8_t included : new_bool_cycle) { + if (included) { + new_cycle.push_back(pos); + } + pos++; + } + cycles.push_back(new_cycle); + } + } + + for (auto cycle : cycles) { + std::array<unsigned int, 2> crossing_count{0,0}; + + for (auto edge : cycle) { + double dx = fabs(n.G.vertices[n.G.edges[edge][0]].r.x - n.G.vertices[n.G.edges[edge][1]].r.x); + if (dx > 0.5) { + crossing_count[0]++; + } + double dy = fabs(n.G.vertices[n.G.edges[edge][0]].r.y - n.G.vertices[n.G.edges[edge][1]].r.y); + if (dy > 0.5) { + crossing_count[1]++; + } + } + + if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { + crack = cycle; + break; + } + } + } else { + crack = cycles.front(); + } + + for (auto edge : crack) { + std::cout << edge << " "; + } + std::cout << "\n"; +// std::cout << max_size << " " << boost::num_edges(Gtmp) << " " << removed_edges.size() << "\n"; + for (unsigned int be = 0; be < log2(L); be++) { unsigned int bin = pow(2, be); @@ -196,13 +364,14 @@ void ma::post_fracture(network &n) { } if (in_bin) { - bin_counts[be]++; + sb[be]++; } } } - bin_counts[log2(L)]++; + sb[log2(L)]++; + // crack surface correlations std::vector<double> crack_damage(pow(L, 2), 0.0); double damage_tot = 0; @@ -214,14 +383,19 @@ void ma::post_fracture(network &n) { } std::vector<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); - std::vector<double> Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + correlation(L, Css, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); - update_field_data("Ccc", damage_tot, Ccc, L, beta); + std::vector<bool> last_crack(pow(L, 2), 0); - for (auto e = avalanches.back().rbegin(); e != avalanches.back().rend(); e++) { + // rewind the last avalanche + for (auto e = avalanche.rbegin(); e != avalanche.rend(); e++) { boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G); + n.add_edge(*e); + last_crack[*e] = 1; } + + // cluster size distribution and cluster-cluster correlation num = boost::connected_components(G, &component[0]); for (unsigned int i = 0; i < num; i++) { @@ -236,13 +410,28 @@ void ma::post_fracture(network &n) { } if (size > 0) { - cd[size - 1]++; + sc[size - 1]++; t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); - Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + correlation(L, Ccc, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + Nc++; + } + } + + // damage-damage correlations + std::vector<fftw_complex> t_damage = data_transform(L, n.fuses, forward_plan, fftw_forward_in, fftw_forward_out); + correlation(L, Cdd, t_damage, t_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + - update_field_data("Cclcl", size, Ccc, L, beta); + // damage size distribution + unsigned int total_broken = 0; + + for (unsigned int i = 0; i < n.G.edges.size(); i++) { + if (n.fuses[i]) { + total_broken++; } } + sd[total_broken]++; + } diff --git a/src/measurements.hpp b/src/measurements.hpp index 817481b..4c05d93 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -7,19 +7,30 @@ #include <string> #include <cinttypes> #include <sstream> +#include <functional> +#include <iostream> +#include <valarray> +#include <array> #include <boost/graph/adjacency_list.hpp> #include <boost/graph/connected_components.hpp> +#include <boost/graph/depth_first_search.hpp> +#include <boost/range/combine.hpp> +#include <boost/foreach.hpp> #include <fftw3.h> #include <network.hpp> #include <hooks.hpp> +struct EdgeProperties { + unsigned int index; +}; + +typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph; + class ma : public hooks { // need: - // - measurements for correlation functions (nice looking?) - // - interface for reading and writing from files (nice looking?) // - interface for turning on and off specific measurements // private: @@ -29,20 +40,28 @@ class ma : public hooks { double *fftw_reverse_out; fftw_plan forward_plan; fftw_plan reverse_plan; - boost::adjacency_list <boost::listS, boost::vecS, boost::undirectedS> G; + Graph G; + + // measurement storage + std::vector<uint64_t> sc; // cluster size distribution + std::vector<uint64_t> sa; // avalanche size distribution + std::vector<uint64_t> sd; // avalanche size distribution + std::vector<uint64_t> sb; // bin size distribution + std::vector<uint64_t> Ccc; // cluster-cluster correlations + std::vector<uint64_t> Css; // surface-surface correlations + std::vector<uint64_t> Caa; // avalanche-avalanche correlations + std::vector<uint64_t> Cdd; // damage-damage distribution + uint64_t Nc; + uint64_t Na; public: unsigned int N; unsigned int L; - unsigned int as; double lv; double beta; - std::vector<unsigned int> ad; - std::vector<unsigned int> cd; - std::vector<unsigned int> bin_counts; - std::list<std::list<unsigned int>> avalanches;; + std::list<unsigned int> avalanche; ma(unsigned int N, unsigned int L, double beta); ~ma(); |