summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/include/graph.hpp11
-rw-r--r--lib/src/graph.cpp217
-rw-r--r--lib/src/network.cpp28
3 files changed, 187 insertions, 69 deletions
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 75d54f4..d746ddd 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -4,8 +4,11 @@
#include <cmath>
#include <vector>
#include <array>
-#include <map>
+#include <unordered_map>
#include <random>
+#include <list>
+#include <exception>
+
class graph {
public:
@@ -16,9 +19,13 @@ class graph {
typedef struct vertex {
coordinate r;
+ std::vector<coordinate> polygon;
} vertex;
- typedef std::array<unsigned int, 2> edge;
+ typedef struct edge {
+ std::array<unsigned int, 2> v;
+ coordinate r;
+ } edge;
coordinate L;
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index f25cf87..10f62b7 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,7 +1,6 @@
#include "graph.hpp"
#include <cstring>
-#include <iostream>
#define JC_VORONOI_IMPLEMENTATION
#define JCV_REAL_TYPE double
@@ -24,9 +23,21 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
for (unsigned int i = 0; i < nv; i++) {
vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
vertices[i].r.y = (double)(i / (Nx / 2));
+ vertices[i].polygon = {
+ {vertices[i].r.x - 0.5, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y - 0.5},
+ {vertices[i].r.x + 0.5, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y + 0.5}
+ };
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].polygon = {
+ {dual_vertices[i].r.x - 0.5, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y - 0.5},
+ {dual_vertices[i].r.x + 0.5, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y + 0.5}
+ };
}
for (unsigned int x = 0; x < Ny; x++) {
@@ -34,40 +45,70 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2);
unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv;
- edges.push_back({v1, v2});
+ edges.push_back({{v1, v2}, {0.5 + (double)y, 0.5 + (double)x}});
unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2);
unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv;
- dual_edges.push_back({dv1, dv2});
+ dual_edges.push_back({{dv1, dv2}, {0.5 + (double)y, 0.5 + (double)x}});
}
}
}
+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;
+ }
+ };
+}
+
+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;
+
+
graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
L = {(double)Nx, (double)Ny};
- unsigned int dnv = Nx * Ny / 2;
+ unsigned int nv = Nx * Ny / 2;
- dual_vertices.resize(dnv);
+ vertices.resize(nv);
std::normal_distribution<double> d(0.0, spread);
- // the coordinates of the dual lattice, from which a delaunay triangulation
- // and corresponding voronoi tesselation will be built
- for (unsigned int i = 0; i < dnv; i++) {
+ // the coordinates of the lattice, from which a delaunay triangulation
+ // and corresponding voronoi tessellation will be built
+ for (unsigned int i = 0; i < nv; i++) {
double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
double ry = (double)(i / (Nx / 2)) + d(rng);
- dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
+ vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
}
- // to make the resulting tesselation periodic, we tile eight (!) copies of
- // the original dual points for a total of nine. note that the index of each
+ // to make the resulting tessellation periodic, we tile eight (!) copies of
+ // the original points for a total of nine. note that the index of each
// point quotient 9 is equal to the original index (we'll use this to
// translate back)
std::vector<jcv_point> points;
- points.reserve(9 * dnv);
- for (const vertex &v : dual_vertices) {
+ points.reserve(9 * nv);
+ for (const vertex &v : vertices) {
points.push_back({v.r.x, v.r.y});
points.push_back({v.r.x + L.x, v.r.y + 0.0});
points.push_back({v.r.x - L.x, v.r.y + 0.0});
@@ -82,81 +123,151 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
jcv_diagram diagram;
memset(&diagram, 0, sizeof(jcv_diagram));
- jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram);
+ jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
const jcv_site* sites = jcv_diagram_get_sites(&diagram);
- struct coorcomp {
- bool operator() (const coordinate& lhs, const coordinate& rhs) const
- {return lhs.x<rhs.x;}
+ struct paircomp {
+ bool operator() (const std::array<unsigned int, 2>& lhs, const std::array<unsigned int, 2>& rhs) const
+ {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1];
+ else return rhs[0] < rhs[1];
+ }
};
- std::map<coordinate, unsigned int, coorcomp> known_vertices;
+ std::unordered_map<std::array<unsigned int, 3>, unsigned int> known_vertices;
for (int i = 0; i < diagram.numsites; i++) {
const jcv_site* site = &sites[i];
// we only care about processing the cells of our original, central sites
if (site->index % 9 == 0) {
+ unsigned int i1 = (unsigned int)(site->index / 9);
const jcv_graphedge* e = site->edges;
+ const jcv_graphedge* ep = site->edges;
+ while (ep->next) {
+ ep = ep->next;
+ }
// for each edge on the site's cell
while(e) {
- // assess whether the dual edge needs to be added. only neighboring
- // dual sites whose index is larger than ours are given edges, so we
+ // 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
const jcv_site* neighbor = e->neighbor;
- unsigned int index = (unsigned int)(site->index / 9);
- unsigned int real_index = (unsigned int)(neighbor->index / 9);
-
- if (index < real_index) {
- dual_edges.push_back({index, real_index});
- coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))};
- coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))};
+ unsigned int i2 = (unsigned int)(neighbor->index / 9);
+
+ vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y});
+
+ unsigned int i3p = (unsigned int)(ep->neighbor->index / 9);
+ std::array<unsigned int, 3> t1 = {i1, i2, i3p};
+ std::sort(t1.begin(), t1.end());
- std::map<coordinate, unsigned int>::iterator it1 = known_vertices.find(r1);
- std::map<coordinate, unsigned int>::iterator it2 = known_vertices.find(r2);
+ auto it1 = known_vertices.find(t1);
- unsigned int vi1, vi2;
+ unsigned int vi1;
- if (it1 == known_vertices.end()) {
- vi1 = vertices.size();
- vertices.push_back({r1});
- known_vertices[r1] = vi1;
+ if (it1 == known_vertices.end()) {
+ vi1 = dual_vertices.size();
+ dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+ known_vertices[t1] = vi1;
+ } else {
+ vi1 = it1->second;
+ dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
+ }
+
+ if (i1 < i2) {
+ edges.push_back({{i1, i2},
+ {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x),
+ fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)}
+ });
+
+ jcv_graphedge *en;
+ if (e->next == NULL) {
+ en = site->edges;
} else {
- vi1 = it1->second;
+ en = e->next;
}
+
+ unsigned int i3n = (unsigned int)(en->neighbor->index / 9);
+ std::array<unsigned int, 3> t2 = {i1, i2, i3n};
+ std::sort(t2.begin(), t2.end());
+
+ auto it2 = known_vertices.find(t2);
+
+ unsigned int vi2;
+
if (it2 == known_vertices.end()) {
- vi2 = vertices.size();
- vertices.push_back({r2});
- known_vertices[r2] = vi2;
+ vi2 = dual_vertices.size();
+ dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}});
+ known_vertices[t2] = vi2;
} else {
vi2 = it2->second;
}
- edges.push_back({vi1, vi2});
+ dual_edges.push_back({{vi1, vi2},
+ {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x),
+ fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)}
+ });
}
+ ep = e;
e = e->next;
}
}
}
- for (edge &e : edges) {
- std::cout << e[0] << " " << e[1] << " ";
- }
- std::cout << "\n";
- for (vertex &v : vertices) {
- std::cout << v.r.x << " " << v.r.y << " ";
+ if (edges.size() != vertices.size() + dual_vertices.size()) {
+ throw eulerex;
}
- std::cout << "\n";
- for (edge &e : dual_edges) {
- std::cout << e[0] << " " << e[1] << " ";
- }
- std::cout << "\n";
- for (vertex &v : dual_vertices) {
- std::cout << v.r.x << " " << v.r.y << " ";
+
+ 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;
+ } else {
+ v.polygon[1].x += L.x;
+ }
+ }
+
+ if (fabs(v.polygon[0].y - v.polygon[1].y) > L.y / 2) {
+ if (v.polygon[0].y < L.y / 2) {
+ v.polygon[0].y += L.y;
+ } else {
+ v.polygon[1].y += L.y;
+ }
+ }
+
+ if (fabs(v.polygon[2].x - v.polygon[1].x) > L.x / 2) {
+ if (v.polygon[2].x < L.x / 2) {
+ v.polygon[2].x += L.x;
+ } else {
+ v.polygon[1].x += L.x;
+ }
+ }
+
+ if (fabs(v.polygon[2].y - v.polygon[1].y) > L.y / 2) {
+ if (v.polygon[2].y < L.y / 2) {
+ v.polygon[2].y += L.y;
+ } else {
+ v.polygon[1].y += L.y;
+ }
+ }
+
+ if (fabs(v.polygon[2].x - v.polygon[0].x) > L.x / 2) {
+ if (v.polygon[2].x < L.x / 2) {
+ v.polygon[2].x += L.x;
+ } else {
+ v.polygon[0].x += L.x;
+ }
+ }
+
+ if (fabs(v.polygon[2].y - v.polygon[0].y) > L.y / 2) {
+ if (v.polygon[2].y < L.y / 2) {
+ v.polygon[2].y += L.y;
+ } else {
+ v.polygon[0].y += L.y;
+ }
+ }
}
- std::cout << "\n";
jcv_diagram_free(&diagram);
}
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 9cb1007..4d4ed2d 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -4,14 +4,14 @@
network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
for (unsigned int i = 0; i < G.edges.size(); i++) {
- double v0y = G.vertices[G.edges[i][0]].r.y;
- double v1y = G.vertices[G.edges[i][1]].r.y;
+ double v0y = G.vertices[G.edges[i].v[0]].r.y;
+ double v1y = G.vertices[G.edges[i].v[1]].r.y;
if (fabs(v0y - v1y) > G.L.y / 2) {
bool ind = v0y < v1y ? 0 : 1;
- ((double *)b->x)[G.edges[i][ind]] += 1.0;
- ((double *)b->x)[G.edges[i][!ind]] -= 1.0;
+ ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
+ ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
}
}
@@ -30,8 +30,8 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
unsigned int terms = G.vertices.size();
for (unsigned int i = 0; i < G.edges.size(); i++) {
- unsigned int v0 = G.edges[i][0];
- unsigned int v1 = G.edges[i][1];
+ unsigned int v0 = G.edges[i].v[0];
+ unsigned int v1 = G.edges[i].v[1];
unsigned int s0 = v0 < v1 ? v0 : v1;
unsigned int s1 = v0 < v1 ? v1 : v0;
@@ -60,11 +60,11 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
for (unsigned int i = 0; i < G.edges.size(); i++) {
((CHOL_INT *)t->i)[2 * i] = i;
- ((CHOL_INT *)t->j)[2 * i] = G.edges[i][0];
+ ((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][1];
+ ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
((double *)t->x)[2 * i + 1] = -1.0;
}
@@ -99,8 +99,8 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
void network::break_edge(unsigned int e, bool unbreak) {
fuses[e] = !unbreak;
- unsigned int v0 = G.edges[e][0];
- unsigned int v1 = G.edges[e][1];
+ unsigned int v0 = G.edges[e].v[0];
+ unsigned int v1 = G.edges[e].v[1];
unsigned int n = factor->n;
@@ -141,8 +141,8 @@ void network::break_edge(unsigned int e, bool unbreak) {
if (fabs(v0y - v1y) > G.L.y / 2) {
bool ind = v0y < v1y ? unbreak : !unbreak;
- ((double *)b->x)[G.edges[e][ind]] -= 1.0;
- ((double *)b->x)[G.edges[e][!ind]] += 1.0;
+ ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
+ ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
}
}
@@ -169,8 +169,8 @@ current_info network::get_current_info() {
} else {
currents[i] = ((double *)y->x)[i];
- double v0y = G.vertices[G.edges[i][0]].r.y;
- double v1y = G.vertices[G.edges[i][1]].r.y;
+ double v0y = G.vertices[G.edges[i].v[0]].r.y;
+ double v1y = G.vertices[G.edges[i].v[1]].r.y;
if (fabs(v0y - v1y) > G.L.y / 2) {
if (v0y > v1y) {