system size now can take on non-integer values

2 files changed, 93 insertions, 33 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index d746ddd..c1ad4f0 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -36,6 +36,6 @@ class graph {
     std::vector<edge> dual_edges;
 
     graph(unsigned int Nx, unsigned int Ny);
-    graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double = 0.25);
+    graph(double Lx, double Ly, std::mt19937& rng, double relax = 0.01, double step = 1.9);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 10f62b7..3d5647f 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -8,6 +8,14 @@
 #define JCV_FLT_MAX 1.7976931348623157E+308
 #include <jc_voronoi.h>
 
+double mod(double a, double b) {
+  if (a >= 0) {
+    return fmod(a, b);
+  } else {
+    return fmod(a + b * ceil(-a / b), b);
+  }
+}
+
 graph::graph(unsigned int Nx, unsigned int Ny) {
   L = {(double)Nx, (double)Ny};
 
@@ -84,45 +92,97 @@ class eulerException: public std::exception
   }
 } eulerex;
 
+graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
+  L = {Lx, Ly};
 
-graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
-  L = {(double)Nx, (double)Ny};
+  std::uniform_real_distribution<double> d(0.0, 1.0);
+  unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5);
 
-  unsigned int nv = Nx * Ny / 2;
+  unsigned int nv = N;
 
   vertices.resize(nv);
-  
-  std::normal_distribution<double> d(0.0, spread);
 
   // 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);
-    vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
-  }
-
-  // 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 * 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});
-    points.push_back({v.r.x + 0.0, v.r.y + L.y});
-    points.push_back({v.r.x + 0.0, v.r.y - L.y});
-    points.push_back({v.r.x + L.x, v.r.y + L.y});
-    points.push_back({v.r.x - L.x, v.r.y + L.y});
-    points.push_back({v.r.x + L.x, v.r.y - L.y});
-    points.push_back({v.r.x - L.x, v.r.y - L.y});
+    vertices[i] = {{L.x * d(rng), L.y * d(rng)}};
   }
 
+  double max_difference = std::numeric_limits<double>::max();
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
+  std::vector<jcv_point> points(9 * nv);
+
+  double rstep = sqrt(step);
+
+  while (max_difference > pow(relax, 2) * 2 * N) {
+    double cur_diff = 0;
+    // 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)
+    for (unsigned int i = 0; i < nv; i++) {
+      const vertex& v = vertices[i];
+      points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
+      points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
+      points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
+      points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
+      points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
+      points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+      points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
+      points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
+      points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+    }
+
+    jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
+    const jcv_site* sites = jcv_diagram_get_sites(&diagram);
+    
+    for (int i = 0; i < diagram.numsites; i++) {
+      const jcv_site* site = &sites[i];
+      unsigned int ind = site->index;
+      if (ind % 9 == 0) {
+        double Cx = 0.0;
+        double Cy = 0.0;
+        unsigned int ne = 0;
+
+        const jcv_graphedge* e = site->edges;
+        while (e) {
+          ne++;
+          Cx += e->pos[0].x;
+          Cy += e->pos[0].y;
+          e = e->next;
+        }
 
+        Cx /= ne;
+        Cy /= ne;
+
+        double dx = Cx - vertices[ind / 9].r.x;
+        double dy = Cy - vertices[ind / 9].r.y;
+
+        double dist = pow(dx, 2) + pow(dy, 2);
+        if (dist > cur_diff) {
+          cur_diff = dist;
+        }
+        vertices[ind / 9] = {{mod(vertices[ind / 9].r.x + rstep * dx, L.x), mod(vertices[ind / 9].r.y + rstep * dy, L.y)}};
+      }
+    }
+    max_difference = cur_diff;
+    jcv_diagram_free(&diagram);
+    memset(&diagram, 0, sizeof(jcv_diagram));
+  }
+
+  for (unsigned int i = 0; i < nv; i++) {
+    const vertex& v = vertices[i];
+    points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
+    points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
+    points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
+    points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
+    points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
+    points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+    points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
+    points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
+    points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+  }
   jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
@@ -167,7 +227,7 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
 
         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}}});
+          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
           known_vertices[t1] = vi1;
         } else {
           vi1 = it1->second;
@@ -176,8 +236,8 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
         
         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)}
+              {mod((site->p.x + neighbor->p.x) / 2, L.x),
+               mod((site->p.y + neighbor->p.y) / 2, L.y)}
                });
 
           jcv_graphedge *en;
@@ -197,15 +257,15 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
 
           if (it2 == known_vertices.end()) {
             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)},{}});
+            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;
           }
 
           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)}
+              {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
+               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)}
                });
         }