summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/network.cpp204
1 files changed, 141 insertions, 63 deletions
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 1edc973..8b0e8df 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -17,14 +17,19 @@ class nofuseException: public std::exception
}
} nofuseex;
-network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
+problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) {
b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
for (unsigned i = 0; i < G.edges.size(); i++) {
- double v0y = G.vertices[G.edges[i].v[0]].r.y;
- double v1y = G.vertices[G.edges[i].v[1]].r.y;
-
- if (G.edges[i].crossings[1]) {
- bool ind = v0y < v1y ? 0 : 1;
+ graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+ graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+ if (G.edges[i].crossings[axis]) {
+ bool ind;
+ if (axis == 1) {
+ ind = v0.y < v1.y ? 0 : 1;
+ } else {
+ ind = v0.x < v1.x ? 0 : 1;
+ }
((double *)b->x)[G.edges[i].v[ind]] += 1.0;
((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
@@ -78,8 +83,10 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
factor = CHOL_F(analyze)(laplacian, c);
CHOL_F(factorize)(laplacian, factor, c);
CHOL_F(free_sparse)(&laplacian, c);
+}
- t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
+problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
+ cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
t->nnz = 2 * G.edges.size();
@@ -98,39 +105,70 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
CHOL_F(free_triplet)(&t, c);
}
-network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) {
+problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) {
b = CHOL_F(copy_dense)(other.b, c);
factor = CHOL_F(copy_factor)(other.factor, c);
voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
}
-network::~network() {
- CHOL_F(free_sparse)(&voltcurmat, c);
+problem::~problem() {
CHOL_F(free_dense)(&b, c);
CHOL_F(free_factor)(&factor, c);
+ CHOL_F(free_sparse)(&voltcurmat, c);
}
-void network::set_thresholds(double beta, std::mt19937& rng) {
- if (beta == 0.0) {
- /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
- for (long double& threshold : thresholds) {
- threshold = 1.0;
- }
- } else {
- std::uniform_real_distribution<long double> d(0.0, 1.0);
+current_info problem::solve(std::vector<bool>& fuses) {
+ cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
- for (long double& threshold : thresholds) {
- threshold = std::numeric_limits<long double>::lowest();
+ if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+ throw nanex;
+ }
- while (threshold == std::numeric_limits<long double>::lowest()) {
- threshold = logl(d(rng)) / (long double)beta;
+ cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+
+ double alpha[2] = {1, 0};
+ double beta[2] = {0, 0};
+ CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+ std::vector<double> currents(G.edges.size());
+
+ double total_current = 0;
+
+ for (int i = 0; i < G.edges.size(); i++) {
+ if (fuses[i]) {
+ currents[i] = 0;
+ } else {
+ currents[i] = ((double *)y->x)[i];
+
+ graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+ graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+ if (G.edges[i].crossings[axis]) {
+ bool comp;
+ if (axis == 1) {
+ comp = v0.y > v1.y;
+ } else {
+ comp = v0.x > v1.x;
+ }
+
+ if (comp) {
+ currents[i] += 1.0;
+ total_current += currents[i];
+ } else {
+ currents[i] -= 1.0;
+ total_current -= currents[i];
+ }
}
}
}
+
+ CHOL_F(free_dense)(&x, c);
+ CHOL_F(free_dense)(&y, c);
+
+ return {total_current, currents};
}
-void network::break_edge(unsigned e, bool unbreak) {
- fuses[e] = !unbreak;
+void problem::break_edge(unsigned e, bool unbreak) {
unsigned v0 = G.edges[e].v[0];
unsigned v1 = G.edges[e].v[1];
@@ -167,78 +205,116 @@ void network::break_edge(unsigned e, bool unbreak) {
CHOL_F(free_sparse)(&perm_update_mat, c);
CHOL_F(free_sparse)(&update_mat, c);
- double v0y = G.vertices[v0].r.y;
- double v1y = G.vertices[v1].r.y;
+ graph::coordinate r0 = G.vertices[v0].r;
+ graph::coordinate r1 = G.vertices[v1].r;
- if (G.edges[e].crossings[1]) {
- bool ind = v0y < v1y ? unbreak : !unbreak;
+ if (G.edges[e].crossings[axis]) {
+ bool ind;
+ if (axis == 1) {
+ ind = r0.y < r1.y ? unbreak : !unbreak;
+ } else {
+ ind = r0.x < r1.x ? unbreak : !unbreak;
+ }
((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
}
}
-current_info network::get_current_info() {
- cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}
- if (((double *)x->x)[0] != ((double *)x->x)[0]) {
- throw nanex;
- }
+network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {}
- cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+void network::set_thresholds(double beta, std::mt19937& rng) {
+ if (beta == 0.0) {
+ /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
+ for (long double& threshold : thresholds) {
+ threshold = 1.0;
+ }
+ } else {
+ std::uniform_real_distribution<long double> d(0.0, 1.0);
- double alpha[2] = {1, 0};
- double beta[2] = {0, 0};
- CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+ for (long double& threshold : thresholds) {
+ threshold = std::numeric_limits<long double>::lowest();
- std::vector<double> currents(G.edges.size());
+ while (threshold == std::numeric_limits<long double>::lowest()) {
+ threshold = logl(d(rng)) / (long double)beta;
+ }
+ }
+ }
+}
- double total_current = 0;
+fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) {
+}
- for (int i = 0; i < G.edges.size(); i++) {
- if (fuses[i]) {
- currents[i] = 0;
- } else {
- currents[i] = ((double *)y->x)[i];
+void fuse_network::fracture(hooks& m, double cutoff) {
+ m.pre_fracture(*this);
- double v0y = G.vertices[G.edges[i].v[0]].r.y;
- double v1y = G.vertices[G.edges[i].v[1]].r.y;
+ while (true) {
+ current_info ci = ohm.solve(fuses);
- if (G.edges[i].crossings[1]) {
- if (v0y > v1y) {
- currents[i] += 1.0;
- total_current += currents[i];
- } else {
- currents[i] -= 1.0;
- total_current -= currents[i];
+ if (ci.conductivity < cutoff * G.vertices.size()) {
+ break;
+ }
+
+ unsigned max_pos = UINT_MAX;
+ long double max_val = std::numeric_limits<long double>::lowest();
+
+ for (unsigned i = 0; i < G.edges.size(); i++) {
+ if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
+ long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+ if (val > max_val) {
+ max_val = val;
+ max_pos = i;
}
}
}
+ if (max_pos == UINT_MAX) {
+ throw nofuseex;
+ }
+
+ fuses[max_pos] = true;
+ ohm.break_edge(max_pos);
+
+ m.bond_broken(*this, ci, max_pos);
}
- CHOL_F(free_dense)(&x, c);
- CHOL_F(free_dense)(&y, c);
+ m.post_fracture(*this);
+}
- return {total_current, currents};
+elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
+}
+
+elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
}
-void network::fracture(hooks& m, double cutoff) {
+void elastic_network::fracture(hooks& m, double weight, double cutoff) {
m.pre_fracture(*this);
while (true) {
- current_info ci = this->get_current_info();
+ current_info cx = hook_x.solve(fuses);
+ current_info cy = hook_y.solve(fuses);
- if (ci.conductivity < cutoff * G.vertices.size()) {
+ current_info ctot;
+
+ ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2));
+
+ if (ctot.conductivity < cutoff * G.vertices.size()) {
break;
}
+ ctot.currents.resize(G.edges.size());
+ for (unsigned i = 0; i < G.edges.size(); i++) {
+ ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i], 2) + pow(weight * cy.currents[i], 2));
+ }
+
unsigned max_pos = UINT_MAX;
long double max_val = std::numeric_limits<long double>::lowest();
for (unsigned i = 0; i < G.edges.size(); i++) {
- if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
- long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+ if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) {
+ long double val = logl(fabs(ctot.currents[i])) - thresholds[i];
if (val > max_val) {
max_val = val;
max_pos = i;
@@ -250,9 +326,11 @@ void network::fracture(hooks& m, double cutoff) {
throw nofuseex;
}
- this->break_edge(max_pos);
+ fuses[max_pos] = true;
+ hook_x.break_edge(max_pos);
+ hook_y.break_edge(max_pos);
- m.bond_broken(*this, ci, max_pos);
+ m.bond_broken(*this, ctot, max_pos);
}
m.post_fracture(*this);