added animation example, and did many fixes to the voronoi system

2 files changed, 178 insertions, 67 deletions

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) {