summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/array_hash.hpp34
-rw-r--r--lib/include/clusters.hpp2
-rw-r--r--lib/include/graph.hpp79
-rw-r--r--lib/include/hooks.hpp9
-rw-r--r--lib/include/network.hpp61
-rw-r--r--lib/include/problem.hpp41
-rw-r--r--lib/src/graph.cpp143
-rw-r--r--lib/src/network.cpp83
-rw-r--r--lib/src/problem.cpp84
9 files changed, 280 insertions, 256 deletions
diff --git a/lib/include/array_hash.hpp b/lib/include/array_hash.hpp
index 3b35e4a..fd0d0f7 100644
--- a/lib/include/array_hash.hpp
+++ b/lib/include/array_hash.hpp
@@ -4,24 +4,18 @@
#include <array>
#include <unordered_map>
-namespace std
-{
- template<typename T, size_t N>
- struct hash<array<T, N> >
- {
- typedef array<T, N> argument_type;
- typedef size_t result_type;
-
- result_type operator()(const argument_type& a) const
- {
- hash<T> hasher;
- result_type h = 0;
- for (result_type i = 0; i < N; ++i)
- {
- h = h * 31 + hasher(a[i]);
- }
- return h;
- }
- };
-}
+namespace std {
+template <typename T, size_t N> struct hash<array<T, N>> {
+ typedef array<T, N> argument_type;
+ typedef size_t result_type;
+ result_type operator()(const argument_type& a) const {
+ hash<T> hasher;
+ result_type h = 0;
+ for (result_type i = 0; i < N; ++i) {
+ h = h * 31 + hasher(a[i]);
+ }
+ return h;
+ }
+};
+} // namespace std
diff --git a/lib/include/clusters.hpp b/lib/include/clusters.hpp
index 4cc9073..0c562d8 100644
--- a/lib/include/clusters.hpp
+++ b/lib/include/clusters.hpp
@@ -1,6 +1,6 @@
-#include <vector>
#include "graph.hpp"
+#include <vector>
class ClusterTree {
private:
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 0452339..9d2b5a1 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -1,50 +1,49 @@
#pragma once
-#include <cmath>
-#include <vector>
+#include <algorithm>
#include <array>
-#include <unordered_map>
-#include <random>
-#include <list>
+#include <cmath>
#include <exception>
-#include <algorithm>
+#include <list>
+#include <random>
+#include <unordered_map>
+#include <vector>
#include "array_hash.hpp"
class graph {
- public:
- typedef struct coordinate {
- double x;
- double y;
- } coordinate;
-
- typedef struct vertex {
- coordinate r;
- std::vector<unsigned> nd;
- std::vector<unsigned> ne;
- std::vector<coordinate> polygon;
- } vertex;
-
- typedef struct edge {
- std::array<unsigned, 2> v;
- coordinate r;
- std::array<bool, 2> crossings;
- } edge;
-
- coordinate L;
-
- std::vector<vertex> vertices;
- std::vector<edge> edges;
-
- std::vector<vertex> dual_vertices;
- std::vector<edge> dual_edges;
-
- graph(unsigned Nx, unsigned Ny);
- graph(double Lx, double Ly, std::mt19937& rng);
- graph(unsigned n, double a, std::mt19937& rng);
-
- void helper(unsigned n, std::mt19937& rng);
- graph const rotate();
+public:
+ typedef struct coordinate {
+ double x;
+ double y;
+ } coordinate;
+
+ typedef struct vertex {
+ coordinate r;
+ std::vector<unsigned> nd;
+ std::vector<unsigned> ne;
+ std::vector<coordinate> polygon;
+ } vertex;
+
+ typedef struct edge {
+ std::array<unsigned, 2> v;
+ coordinate r;
+ std::array<bool, 2> crossings;
+ } edge;
+
+ coordinate L;
+
+ std::vector<vertex> vertices;
+ std::vector<edge> edges;
+
+ std::vector<vertex> dual_vertices;
+ std::vector<edge> dual_edges;
+
+ graph(unsigned Nx, unsigned Ny);
+ graph(double Lx, double Ly, std::mt19937& rng);
+ graph(unsigned n, double a, std::mt19937& rng);
+
+ void helper(unsigned n, std::mt19937& rng);
+ graph const rotate();
};
-
diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp
index 67313cc..ba85ff7 100644
--- a/lib/include/hooks.hpp
+++ b/lib/include/hooks.hpp
@@ -6,9 +6,8 @@
class network;
class hooks {
- public:
- virtual void pre_fracture(const network&) {};
- virtual void bond_broken(const network&, const current_info&, unsigned) {};
- virtual void post_fracture(network&) {}; // post fracture hook can be destructive
+public:
+ virtual void pre_fracture(const network&){};
+ virtual void bond_broken(const network&, const current_info&, unsigned){};
+ virtual void post_fracture(network&){}; // post fracture hook can be destructive
};
-
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 8748ea9..29bf55a 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -1,20 +1,20 @@
#pragma once
-#include <vector>
#include <functional>
-#include <utility>
-#include <random>
#include <limits>
-#include <valarray>
+#include <random>
#include <set>
+#include <utility>
+#include <valarray>
+#include <vector>
-#include "problem.hpp"
-#include "graph.hpp"
-#include "hooks.hpp"
-#include "current_info.hpp"
#include "array_hash.hpp"
#include "clusters.hpp"
+#include "current_info.hpp"
+#include "graph.hpp"
+#include "hooks.hpp"
+#include "problem.hpp"
#define CURRENT_CUTOFF 1e-10
@@ -27,7 +27,7 @@ class network {
friend class elastic_network;
friend class percolation_network;
- private:
+private:
problem px;
problem py;
std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
@@ -36,14 +36,15 @@ class network {
void break_edge(unsigned, bool unbreak = false);
void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
std::set<unsigned> get_cycle_edges(unsigned) const;
- bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned, unsigned) const;
+ bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned,
+ unsigned) const;
std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
void get_cluster_edges_helper(std::set<unsigned>&, unsigned) const;
std::set<unsigned> get_cluster_edges(unsigned) const;
void get_tie_flaps_helper(std::set<unsigned>&, unsigned, unsigned) const;
std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
- public:
+public:
const graph& G;
ClusterTree C;
@@ -57,36 +58,38 @@ class network {
void set_thresholds(double, std::mt19937&);
void fracture(hooks&, bool one_axis = true);
- virtual current_info get_current_info() {current_info empty; return empty;};
+ virtual current_info get_current_info() {
+ current_info empty;
+ return empty;
+ };
};
class fuse_network : public network {
- private:
- double weight;
+private:
+ double weight;
- public:
- fuse_network(const graph&, cholmod_common*, double weight = 1.0);
- fuse_network(const fuse_network&);
+public:
+ fuse_network(const graph&, cholmod_common*, double weight = 1.0);
+ fuse_network(const fuse_network&);
- current_info get_current_info() override;
+ current_info get_current_info() override;
};
class elastic_network : public network {
- private:
- double weight;
+private:
+ double weight;
- public:
- elastic_network(const graph&, cholmod_common*, double weight = 0.5);
- elastic_network(const elastic_network&);
+public:
+ elastic_network(const graph&, cholmod_common*, double weight = 0.5);
+ elastic_network(const elastic_network&);
- current_info get_current_info() override;
+ current_info get_current_info() override;
};
class percolation_network : public network {
- public:
- percolation_network(const graph&, cholmod_common*);
- percolation_network(const percolation_network&);
+public:
+ percolation_network(const graph&, cholmod_common*);
+ percolation_network(const percolation_network&);
- current_info get_current_info() override;
+ current_info get_current_info() override;
};
-
diff --git a/lib/include/problem.hpp b/lib/include/problem.hpp
index 740751d..c454a35 100644
--- a/lib/include/problem.hpp
+++ b/lib/include/problem.hpp
@@ -1,10 +1,10 @@
-#include <vector>
#include <limits>
+#include <vector>
-#include <cholmod.h>
-#include "graph.hpp"
#include "current_info.hpp"
+#include "graph.hpp"
+#include <cholmod.h>
#ifdef FRACTURE_LONGINT
@@ -19,22 +19,21 @@
#endif
class problem {
- private:
- const graph& G;
- unsigned axis;
- cholmod_dense* b;
- cholmod_factor* factor;
- cholmod_sparse* voltcurmat;
- cholmod_common* c;
-
- public:
- current_info sol;
-
- problem(const graph&, unsigned, cholmod_common*);
- problem(const graph&, unsigned, cholmod_sparse*, cholmod_common*);
- problem(const problem&);
- ~problem();
- void solve(std::vector<bool>& fuses);
- void break_edge(unsigned, bool unbreak = false);
+private:
+ const graph& G;
+ unsigned axis;
+ cholmod_dense* b;
+ cholmod_factor* factor;
+ cholmod_sparse* voltcurmat;
+ cholmod_common* c;
+
+public:
+ current_info sol;
+
+ problem(const graph&, unsigned, cholmod_common*);
+ problem(const graph&, unsigned, cholmod_sparse*, cholmod_common*);
+ problem(const problem&);
+ ~problem();
+ void solve(std::vector<bool>& fuses);
+ void break_edge(unsigned, bool unbreak = false);
};
-
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index a5063de..1475c40 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -34,23 +34,21 @@ graph::graph(unsigned Nx, unsigned Ny) {
vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
vertices[i].r.y = (double)(i / (Nx / 2));
signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
- vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
- vertices[i].polygon = {
- {vertices[i].r.x - 1.0, vertices[i].r.y},
- {vertices[i].r.x, vertices[i].r.y - 1.0},
- {vertices[i].r.x + 1.0, vertices[i].r.y},
- {vertices[i].r.x, vertices[i].r.y + 1.0}
- };
+ vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)),
+ (nv + i - Nx / 2) % nv};
+ vertices[i].polygon = {{vertices[i].r.x - 1.0, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y - 1.0},
+ {vertices[i].r.x + 1.0, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y + 1.0}};
dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
dual_vertices[i].r.y = (double)(i / (Nx / 2));
- dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
- dual_vertices[i].polygon = {
- {dual_vertices[i].r.x - 1.0, vertices[i].r.y},
- {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
- {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
- {dual_vertices[i].r.x, vertices[i].r.y + 1.0}
- };
+ dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)),
+ (nv + i - Nx / 2) % nv};
+ dual_vertices[i].polygon = {{dual_vertices[i].r.x - 1.0, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
+ {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y + 1.0}};
}
for (unsigned y = 0; y < Ny; y++) {
@@ -66,7 +64,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
unsigned dv1 = (Nx * y) / 2 + ((x + (y + 1) % 2) / 2) % (Nx / 2);
unsigned dv2 = ((Nx * (y + 1)) / 2 + ((x + y % 2) / 2) % (Nx / 2)) % nv;
- dual_edges.push_back({{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
+ dual_edges.push_back(
+ {{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
}
}
@@ -102,29 +101,23 @@ graph::graph(unsigned Nx, unsigned Ny) {
dual_vertices[vi].ne.push_back(i);
}
}
-
}
-class eulerException: public std::exception
-{
- virtual const char* what() const throw()
- {
- return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed.";
+class eulerException : public std::exception {
+ virtual const char* what() const throw() {
+ return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus "
+ "and is malformed.";
}
} eulerex;
-class clippingException: public std::exception
-{
- virtual const char* what() const throw()
- {
+class clippingException : public std::exception {
+ virtual const char* what() const throw() {
return "An interior site has a clipped edge and the tesselation is malformed.";
}
} clippingex;
-class triangleException: public std::exception
-{
- virtual const char* what() const throw()
- {
+class triangleException : public std::exception {
+ virtual const char* what() const throw() {
return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
}
} triex;
@@ -147,29 +140,45 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
if ((j1 == j2) && (j2 == j3)) {
return 0;
- } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) ||
+ ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
return 1;
- } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) ||
+ ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
return 2;
- } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
+ } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) ||
+ ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
return 3;
- } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
+ } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) ||
+ ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
return 4;
- } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) ||
+ ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
return 5;
- } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) ||
+ ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
return 6;
- } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) ||
+ ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
return 7;
- } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) ||
+ ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
return 8;
- } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) ||
+ ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) ||
+ ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
return 9;
- } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) ||
+ ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) ||
+ ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
return 10;
- } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) ||
+ ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) ||
+ ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
return 11;
- } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) ||
+ ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) ||
+ ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
return 12;
} else {
throw triex;
@@ -200,7 +209,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
// the coordinates of the lattice, from which a delaunay triangulation
// and corresponding voronoi tessellation will be built. Everyone is in the
// rectangle (0, 0) (Lx, Ly)
- for (vertex &v : vertices) {
+ for (vertex& v : vertices) {
v = {{L.x * d(rng), L.y * d(rng)}};
}
@@ -243,7 +252,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
ep = ep->next;
}
// for each edge on the site's cell
- while(e) {
+ while (e) {
// assess whether the edge needs to be added. only neighboring
// sites whose index is larger than ours are given edges, so we
// don't double up later
@@ -275,7 +284,8 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
if (it1 == known_vertices.end()) {
vi1 = dual_vertices.size();
- dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
+ dual_vertices.push_back(
+ {{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)}, {i1}, {}, {{site->p.x, site->p.y}}});
known_vertices[t1] = vi1;
} else {
vi1 = it1->second;
@@ -284,7 +294,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
}
vertices[i1].nd.push_back(vi1);
-
+
if (i1 < i2 || (i1 == i2 && !self_bonded)) {
if (i1 == i2) {
self_bonded = true;
@@ -292,12 +302,11 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
bool crossed_x = x2 != 1;
bool crossed_y = y2 != 1;
edges.push_back({{i1, i2},
- {mod((site->p.x + neighbor->p.x) / 2, L.x),
- mod((site->p.y + neighbor->p.y) / 2, L.y)},
- {crossed_x, crossed_y}
- });
+ {mod((site->p.x + neighbor->p.x) / 2, L.x),
+ mod((site->p.y + neighbor->p.y) / 2, L.y)},
+ {crossed_x, crossed_y}});
- jcv_graphedge *en;
+ jcv_graphedge* en;
if (e->next == NULL) {
en = site->edges;
} else {
@@ -322,20 +331,21 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
if (it2 == known_vertices.end()) {
vi2 = dual_vertices.size();
- dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}});
+ dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)}, {}});
known_vertices[t2] = vi2;
} else {
vi2 = it2->second;
}
- bool dcrossed_x = (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
- bool dcrossed_y = (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
+ bool dcrossed_x =
+ (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
+ bool dcrossed_y =
+ (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
dual_edges.push_back({{vi1, vi2},
- {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
- mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
- {dcrossed_x, dcrossed_y}
- });
+ {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
+ mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
+ {dcrossed_x, dcrossed_y}});
}
ep = e;
@@ -348,7 +358,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
throw eulerex;
}
- for (vertex &v : dual_vertices) {
+ for (vertex& v : dual_vertices) {
if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) {
if (v.polygon[0].x < L.x / 2) {
v.polygon[0].x += L.x;
@@ -434,33 +444,31 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
jcv_diagram_free(&diagram);
}
-graph::coordinate reverse(const graph::coordinate &x) {
- return {x.y, x.x};
-}
+graph::coordinate reverse(const graph::coordinate& x) { return {x.y, x.x}; }
graph const graph::rotate() {
graph g2(*this);
- for (graph::vertex &v : g2.vertices) {
+ for (graph::vertex& v : g2.vertices) {
v.r = reverse(v.r);
- for (graph::coordinate &r : v.polygon) {
+ for (graph::coordinate& r : v.polygon) {
r = reverse(r);
}
}
- for (graph::edge &e : g2.edges) {
+ for (graph::edge& e : g2.edges) {
e.r = reverse(e.r);
e.crossings = {e.crossings[1], e.crossings[0]};
}
- for (graph::vertex &v : g2.dual_vertices) {
+ for (graph::vertex& v : g2.dual_vertices) {
v.r = reverse(v.r);
- for (graph::coordinate &r : v.polygon) {
+ for (graph::coordinate& r : v.polygon) {
r = reverse(r);
}
}
- for (graph::edge &e : g2.dual_edges) {
+ for (graph::edge& e : g2.dual_edges) {
e.r = reverse(e.r);
e.crossings = {e.crossings[1], e.crossings[0]};
}
@@ -469,4 +477,3 @@ graph const graph::rotate() {
return g2;
}
-
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 0cf50c7..eef7e9a 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -11,7 +11,8 @@ network::network(const graph& G, cholmod_common* 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" */
+ /* zero beta doesn't make any sense computationally, we interpret it as
+ * "infinity" */
for (long double& threshold : thresholds) {
threshold = 1.0;
}
@@ -104,7 +105,8 @@ std::set<unsigned> network::get_cycle_edges(unsigned v0) const {
return cycle_edges;
}
-bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const {
+bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges,
+ unsigned vPrev, unsigned vCur, unsigned vEnd) const {
for (unsigned ei : G.dual_vertices[vCur].ne) {
if (fuses[ei]) {
auto it = cycle_edges.find(ei);
@@ -113,13 +115,17 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns
unsigned vn = e[0] == vCur ? e[1] : e[0];
if (vn != vPrev) {
if (vn == vEnd) {
- if (G.dual_edges[ei].crossings[0]) sig[0]++;
- if (G.dual_edges[ei].crossings[1]) sig[1]++;
+ if (G.dual_edges[ei].crossings[0])
+ sig[0]++;
+ if (G.dual_edges[ei].crossings[1])
+ sig[1]++;
return true;
} else {
if (this->find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
- if (G.dual_edges[ei].crossings[0]) sig[0]++;
- if (G.dual_edges[ei].crossings[1]) sig[1]++;
+ if (G.dual_edges[ei].crossings[0])
+ sig[0]++;
+ if (G.dual_edges[ei].crossings[1])
+ sig[1]++;
return true;
}
}
@@ -131,7 +137,8 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns
return false;
}
-std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
+std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0,
+ unsigned v1) const {
std::array<unsigned, 2> sig = {0, 0};
this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
return sig;
@@ -158,7 +165,8 @@ std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
return cluster_edges;
}
-void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0, unsigned vCur) const {
+void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0,
+ unsigned vCur) const {
for (unsigned ei : G.vertices[vCur].ne) {
if (!backbone[ei]) {
auto it = added_edges.find(ei);
@@ -219,7 +227,8 @@ void network::update_backbone(const std::vector<double>& c) {
// Now, we create a crossing signature for each cycle. First, we do it
// for the new cycle that would form by removing the stem.
- std::array<unsigned, 2> base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
+ std::array<unsigned, 2> base_sig =
+ this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
if (G.dual_edges[i].crossings[0])
base_sig[0]++;
if (G.dual_edges[i].crossings[1])
@@ -229,9 +238,12 @@ void network::update_backbone(const std::vector<double>& c) {
// previous step.
std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
for (unsigned e : cedges) {
- std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]);
- std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]);
- std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+ std::array<unsigned, 2> new_sig_1 =
+ this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]);
+ std::array<unsigned, 2> new_sig_2 =
+ this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]);
+ std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0],
+ new_sig_1[1] + new_sig_2[1]};
if (G.dual_edges[i].crossings[0])
new_sig[0]++;
@@ -248,10 +260,10 @@ void network::update_backbone(const std::vector<double>& c) {
// Now, having found all cycles involving the candidate stem, we check
// if any of them spans the torus and therefore can carry current.
if (std::any_of(all_sigs.begin(), all_sigs.end(),
- [done_x, done_y](std::array<unsigned, 2> c) {
- return (c[0] % 2 == 0 && c[1] % 2 == 0) ||
- ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));
- })) {
+ [done_x, done_y](std::array<unsigned, 2> c) {
+ return (c[0] % 2 == 0 && c[1] % 2 == 0) ||
+ ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));
+ })) {
// If it does, we remove it from the backbone!
seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
backbone[i] = true;
@@ -297,6 +309,7 @@ void network::update_backbone(const std::vector<double>& c) {
// Now, we will check for bow ties!
std::set<unsigned> bb_verts;
+ // First we get a list of all vertices that remain in the backbone.
for (unsigned i = 0; i < G.edges.size(); i++) {
if (!backbone[i]) {
bb_verts.insert(G.edges[i].v[0]);
@@ -307,6 +320,10 @@ void network::update_backbone(const std::vector<double>& c) {
for (unsigned v : bb_verts) {
bool found_pair = false;
unsigned d1, d2, p1, p2;
+ // For each vertex, we check to see if two of the faces adjacent to the
+ // vertex are in the same component of the dual lattice and if so, we make
+ // sure they do not belong to a contiguously connected series of such
+ // faces.
for (unsigned i = 0; i < G.vertices[v].nd.size(); i++) {
for (unsigned j = 2; j < G.vertices[v].nd.size() - 1; j++) {
unsigned l = (i + j) % G.vertices[v].nd.size();
@@ -316,12 +333,11 @@ void network::update_backbone(const std::vector<double>& c) {
auto it_search = seen_pairs.find({G.vertices[v].nd[il], G.vertices[v].nd[ig]});
if (it_search == seen_pairs.end()) {
bool any_intervening1 = false;
- bool any_intervening_e_1 = false;
bool any_intervening2 = false;
- bool any_intervening_e_2 = false;
unsigned ie1 = 0;
unsigned ie2 = 0;
+ // yuck, think of something more elegant?
for (unsigned k = il + 1; k < ig; k++) {
if (!C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[k])) {
any_intervening1 = true;
@@ -330,7 +346,7 @@ void network::update_backbone(const std::vector<double>& c) {
}
for (unsigned k = (ig + 1); k % G.vertices[v].nd.size() != il; k++) {
if (!C.same_component(G.vertices[v].nd[i],
- G.vertices[v].nd[k % G.vertices[v].nd.size()])) {
+ G.vertices[v].nd[k % G.vertices[v].nd.size()])) {
any_intervening2 = true;
break;
}
@@ -349,6 +365,9 @@ void network::update_backbone(const std::vector<double>& c) {
}
}
}
+
+ // If all these conditions are true, we process the pair. The same
+ // cycle analysis as in the lollipop case is done.
if ((any_intervening1 && any_intervening2) || (any_intervening1 && ie2 > 1) ||
(any_intervening2 && ie1 > 1)) {
found_pair = true;
@@ -367,13 +386,14 @@ void network::update_backbone(const std::vector<double>& c) {
base_sig[1]++;
}
- // Then, we do it for each of the existing cycles we found in the
- // previous step.
std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
for (unsigned e : cedges) {
- std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]);
- std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]);
- std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+ std::array<unsigned, 2> new_sig_1 =
+ this->find_cycle(cedges, d1, G.dual_edges[e].v[0]);
+ std::array<unsigned, 2> new_sig_2 =
+ this->find_cycle(cedges, d2, G.dual_edges[e].v[1]);
+ std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0],
+ new_sig_1[1] + new_sig_2[1]};
for (unsigned k = p1; k < p2; k++) {
if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
@@ -390,15 +410,20 @@ void network::update_backbone(const std::vector<double>& c) {
}
if (std::any_of(all_sigs.begin(), all_sigs.end(),
- [done_x, done_y](std::array<unsigned, 2> c) {
- return ((c[0] % 2 == 0 && c[1] % 2 == 0) ||
- (c[0] % 2 == 0 && done_x)) ||
- (c[1] % 2 == 0 && done_y);
- })) {
+ [done_x, done_y](std::array<unsigned, 2> c) {
+ return ((c[0] % 2 == 0 && c[1] % 2 == 0) ||
+ (c[0] % 2 == 0 && done_x)) ||
+ (c[1] % 2 == 0 && done_y);
+ })) {
+ // one of our pairs qualifies for trimming!
seen_pairs[{d1, d2}] = true;
+ // We separate the flaps of the bowtie (there may be more than
+ // two!) into distinct sets.
std::list<std::set<unsigned>> flaps = this->get_tie_flaps(v);
+ // All the bonds in each flap without current are removed from
+ // the backbone.
for (const std::set<unsigned>& flap : flaps) {
std::array<double, 2> total_current = {0.0, 0.0};
for (unsigned e : flap) {
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
index 0645244..ab87b0c 100644
--- a/lib/src/problem.cpp
+++ b/lib/src/problem.cpp
@@ -1,15 +1,12 @@
#include "problem.hpp"
-class nanException: public std::exception
-{
- virtual const char* what() const throw()
- {
- return "The linear problem returned NaN.";
- }
+class nanException : public std::exception {
+ virtual const char* what() const throw() { return "The linear problem returned NaN."; }
} nanex;
-problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), axis(axis), voltcurmat(vcmat), c(c) {
+problem::problem(const graph& G, unsigned axis, cholmod_sparse* vcmat, cholmod_common* c)
+ : G(G), axis(axis), voltcurmat(vcmat), c(c) {
b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
for (unsigned i = 0; i < G.edges.size(); i++) {
graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
@@ -23,19 +20,20 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
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;
+ ((double*)b->x)[G.edges[i].v[ind]] += 1.0;
+ ((double*)b->x)[G.edges[i].v[!ind]] -= 1.0;
}
}
unsigned nnz = G.vertices.size() + G.edges.size();
- cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c);
+ cholmod_triplet* t =
+ CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c);
for (unsigned i = 0; i < G.vertices.size(); i++) {
- ((CHOL_INT *)t->i)[i] = i;
- ((CHOL_INT *)t->j)[i] = i;
- ((double *)t->x)[i] = 0.0;
+ ((CHOL_INT*)t->i)[i] = i;
+ ((CHOL_INT*)t->j)[i] = i;
+ ((double*)t->x)[i] = 0.0;
}
unsigned terms = G.vertices.size();
@@ -46,8 +44,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
unsigned v0 = G.edges[i].v[0];
unsigned v1 = G.edges[i].v[1];
- ((double *)t->x)[v0]++;
- ((double *)t->x)[v1]++;
+ ((double*)t->x)[v0]++;
+ ((double*)t->x)[v1]++;
unsigned s0 = v0 < v1 ? v0 : v1;
unsigned s1 = v0 < v1 ? v1 : v0;
@@ -55,22 +53,22 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
auto it = known_edges.find({s0, s1});
if (it == known_edges.end()) {
- ((CHOL_INT *)t->i)[terms] = s1;
- ((CHOL_INT *)t->j)[terms] = s0;
- ((double *)t->x)[terms] = -1.0;
+ ((CHOL_INT*)t->i)[terms] = s1;
+ ((CHOL_INT*)t->j)[terms] = s0;
+ ((double*)t->x)[terms] = -1.0;
known_edges[{s0, s1}] = terms;
terms++;
} else {
- ((double *)t->x)[it->second] -= 1.0;
+ ((double*)t->x)[it->second] -= 1.0;
}
}
- ((double *)t->x)[0]++;
+ ((double*)t->x)[0]++;
t->nnz = terms;
- cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
+ cholmod_sparse* laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
CHOL_F(free_triplet)(&t, c);
factor = CHOL_F(analyze)(laplacian, c);
CHOL_F(factorize)(laplacian, factor, c);
@@ -79,19 +77,20 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c
sol.currents.resize(G.edges.size());
}
-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);
+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();
for (unsigned i = 0; i < G.edges.size(); i++) {
- ((CHOL_INT *)t->i)[2 * i] = i;
- ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0];
- ((double *)t->x)[2 * i] = 1.0;
+ ((CHOL_INT*)t->i)[2 * i] = i;
+ ((CHOL_INT*)t->j)[2 * i] = G.edges[i].v[0];
+ ((double*)t->x)[2 * i] = 1.0;
- ((CHOL_INT *)t->i)[2 * i + 1] = i;
- ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
- ((double *)t->x)[2 * i + 1] = -1.0;
+ ((CHOL_INT*)t->i)[2 * i + 1] = i;
+ ((CHOL_INT*)t->j)[2 * i + 1] = G.edges[i].v[1];
+ ((double*)t->x)[2 * i + 1] = -1.0;
}
voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c);
@@ -103,7 +102,7 @@ problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c
b = CHOL_F(copy_dense)(other.b, c);
factor = CHOL_F(copy_factor)(other.factor, c);
voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
-}
+}
problem::~problem() {
CHOL_F(free_dense)(&b, c);
@@ -112,13 +111,13 @@ problem::~problem() {
}
void problem::solve(std::vector<bool>& fuses) {
- cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+ cholmod_dense* x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
- if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+ if (((double*)x->x)[0] != ((double*)x->x)[0]) {
throw nanex;
}
- cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+ 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};
@@ -130,7 +129,7 @@ void problem::solve(std::vector<bool>& fuses) {
if (fuses[i]) {
sol.currents[i] = 0;
} else {
- sol.currents[i] = ((double *)y->x)[i];
+ sol.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;
@@ -166,15 +165,15 @@ void problem::break_edge(unsigned e, bool unbreak) {
unsigned n = factor->n;
- cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+ cholmod_sparse* update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
unsigned s1, s2;
s1 = v0 < v1 ? v0 : v1;
s2 = v0 < v1 ? v1 : v0;
- CHOL_INT *pp = (CHOL_INT *)update_mat->p;
- CHOL_INT *ii = (CHOL_INT *)update_mat->i;
- double *xx = (double *)update_mat->x;
+ CHOL_INT* pp = (CHOL_INT*)update_mat->p;
+ CHOL_INT* ii = (CHOL_INT*)update_mat->i;
+ double* xx = (double*)update_mat->x;
for (unsigned i = 0; i <= s1; i++) {
pp[i] = 0;
@@ -189,8 +188,8 @@ void problem::break_edge(unsigned e, bool unbreak) {
xx[0] = 1;
xx[1] = -1;
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
+ cholmod_sparse* perm_update_mat =
+ CHOL_F(submatrix)(update_mat, (CHOL_INT*)factor->Perm, factor->n, NULL, -1, true, true, c);
CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
@@ -208,8 +207,7 @@ void problem::break_edge(unsigned e, bool unbreak) {
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;
+ ((double*)b->x)[G.edges[e].v[ind]] -= 1.0;
+ ((double*)b->x)[G.edges[e].v[!ind]] += 1.0;
}
}
-