summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/measurements.cpp297
-rw-r--r--src/measurements.hpp35
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();