From f861195c1416220c4039bda4d1cbf8c3aab07528 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sun, 24 Feb 2019 17:19:52 -0500
Subject: added support for specifiying voronoi graphs with number of points
 and aspect ratio, and added support for "infinite" beta (everything has the
 some threshold)

---
 lib/include/graph.hpp |  3 +++
 lib/src/graph.cpp     | 21 ++++++++++++++++-----
 lib/src/network.cpp   | 17 ++++++++++++-----
 3 files changed, 31 insertions(+), 10 deletions(-)

(limited to 'lib')

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 1089925..032a4d1 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -39,5 +39,8 @@ class graph {
 
     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);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 526ffeb..8aa8877 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -141,27 +141,38 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
 }
 
 graph::graph(double Lx, double Ly, std::mt19937& rng) {
-  L = {Lx, Ly};
-
   // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
   // probability proportional to the distance from each
   std::uniform_real_distribution<double> d(0.0, 1.0);
   unsigned N = round(Lx * Ly / 2 + d(rng) - 0.5);
 
-  unsigned nv = N;
+  L = {Lx, Ly};
+
+  this->helper(N, rng);
+}
+
+graph::graph(unsigned n, double a, std::mt19937& rng) {
+  L = {2 * n * a, 2 * n / a};
+
+  this->helper(n, rng);
+}
+
+void graph::helper(unsigned n, std::mt19937& rng) {
+  std::uniform_real_distribution<double> d(0.0, 1.0);
+  unsigned nv = n;
   vertices.resize(nv);
 
   // 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) {
-    v = {{Lx * d(rng), Ly * d(rng)}};
+    v = {{L.x * d(rng), L.y * d(rng)}};
   }
 
   // set up the voronoi objects
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
-  jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}};
+  jcv_rect bounds = {{-L.x, -L.y}, {2 * L.x, 2 * L.y}};
   std::vector<jcv_point> points(9 * nv);
 
   for (unsigned i = 0; i < nv; i++) {
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 027946a..43b0fab 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -95,13 +95,20 @@ network::~network() {
 }
 
 void network::set_thresholds(double beta, std::mt19937& rng) {
-  std::uniform_real_distribution<long double> d(0.0, 1.0);
+  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);
 
-  for (long double& threshold : thresholds) {
-    threshold = std::numeric_limits<long double>::lowest();
+    for (long double& threshold : thresholds) {
+      threshold = std::numeric_limits<long double>::lowest();
 
-    while (threshold == std::numeric_limits<long double>::lowest()) {
-      threshold = logl(d(rng)) / (long double)beta;
+      while (threshold == std::numeric_limits<long double>::lowest()) {
+        threshold = logl(d(rng)) / (long double)beta;
+      }
     }
   }
 }
-- 
cgit v1.2.3-54-g00ecf