std::vectors

bary.c++

@@ -2,7 +2,7 @@
 #include <cmath>
 #include <cstdio>
 
-#include "si.h"
+#include "bary.h"
 
 // Solve Ax=b, where A is d*d.  Returns false iff A is singular.
 // Crout's algorithm replaces A with its LU decomposition.
@@ -137,19 +137,20 @@ double determinant(const double* m, int n) {
   return det;
 }
 
-bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
-  const vertex* rgv, const vertex* raysCentroid, [[maybe_unused]] bool fRaySimplex)
+simplexHint precomputeBary(const d_simplex& s, const vertex& centroid,
+  const std::vector<vertex>& rgv, const vertex* raysCentroid, [[maybe_unused]] bool fRaySimplex)
 {
+  simplexHint h = {}; // h.s == nullptr, so h is (not yet) valid.
   if (d <= 1)
-    return false;
+    return h;
 
 #ifdef TESTING
   {
     d_simplex tmp(s);
     // s was sorted by sort_output().
     if (std::unique(tmp.begin(), tmp.end()) != tmp.end()) {
       dump_simplex("Internal error: duplicate vertices in simplex:", s);
-      return false;
+      return h;
       }
   }
 #endif
@@ -204,12 +205,9 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
         printf("internal error: inconsistency #3: %d %d in precomputeBary.\n", i, d);
       schmoo[j] = ivertex;
 #endif
-      if (j >= d)
-        printf("internal error in precomputeBary!\n");
       rgvFacet[j] = ivertex<0 ? raysCentroid : &rgv[ivertex];
       if (i != ifacet)
         if (++j >= d)
-	  // Already finished;  quit now without overflowing.
 	  break;
       }
 
@@ -220,8 +218,8 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
 #endif
 
     vertex& vNormal = h.facetnormal[ifacet];
-    for (i=0; i<d; ++i)
-      vNormal[i] = i==0 ? 1. : 0.;
+    vNormal.fill(0.0);
+    vNormal[0] = 1.0;
 
     double a[d][d];
     for (j=0; j<d; ++j)
@@ -237,29 +235,28 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
       printf("Internal error in precomputeBary().\n");
       // Matrix a was probably singular.  Because s had degenerate faces.
       dump_simplex("possibly degenerate simplex:", s);
-      return false;
+      return h;
       }
 
     // Normalize vNormal to unit length (unfortunate pun).
     {
-    double t = 0.;
-    for (i=0; i<d; ++i)
-      t += sq(vNormal[i]);
-    t = sqrt(t);
-    for (i=0; i<d; ++i)
-      vNormal[i] /= t;
+      auto t = 0.0;
+      for (auto x: vNormal) t += sq(x);
+      t = sqrt(t);
+      for (auto& x: vNormal) x /= t;
     }
 
     // Compute vNormal dot (centroid - one_of_the_facet_points).
-    double t = 0.;
-    for (i=0; i<d; ++i)
-      t += vNormal[i] * (centroid[i] - (*rgvFacet[0])[i]);
-    if (fabs(t) < .0001)
-      printf("Warning: found degenerate simplex.\n");
-    else if (t < 0.)
-      // Reverse normal to orient it inwards w.r.t. the simplex.
+    {
+      auto t = 0.0;
       for (i=0; i<d; ++i)
-	vNormal[i] *= -1.;
+	t += vNormal[i] * (centroid[i] - (*rgvFacet[0])[i]);
+      if (fabs(t) < 0.0001)
+	printf("Warning: found degenerate simplex.\n");
+      else if (t < 0.0)
+	// Reverse normal to orient it inwards w.r.t. the simplex.
+	for (auto& x: vNormal) x *= -1.0;
+    }
 
 #ifdef TESTING
     {
@@ -284,7 +281,7 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
         printf("Internal error in precomputeBary: vNormal isn't normal to facet %d's row %d:\n", ifacet, i);
 	printf("  (%g, %g) dot (%g, %g) = %g, not zero.\n",
 	    a[i][0], a[i][1], vNormal[0], vNormal[1], t);
-	return false;
+	return h;
 	}
       }
     }
@@ -297,17 +294,18 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
 
     // Optimization: if we cached/hashed these, we'd avoid computing each facet twice.
 
-    // Compute (d-1)-volume of each facet (defined by rgvFacet[0,...,d-1] ).
+    // Compute the (d-1)-volume of each facet rgvFacet[0,...,d-1].
     switch (d)
       {
     case 2:
 	{
 	// Fast: Euclidean distance from rgvFacet[0] to rgvFacet[1].
 	const auto& a = *rgvFacet[0];
 	const auto& b = *rgvFacet[1];
-	h.facetvolume[ifacet] = hypot(b[0]-a[0], b[1]-a[1]); // these values are OK, 8/29/02
+	h.facetvolume[ifacet] = hypot(b[0]-a[0], b[1]-a[1]);
 	break;
 	}
+
     case 3:
 	{
 	// Fast: Area of triangle.
@@ -324,10 +322,9 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
 	h.facetvolume[ifacet] = sqrt(s * (s-la) * (s-lb) * (s-lc));
 	break;
 	}
-    default:
-      // Compute volume via a generalization of Heron's formula,
-      // the Cayley-Menger determinant.
 
+    default:
+      // Use a generalization of Heron's formula, the Cayley-Menger determinant.
       double m[(d+1)*(d+1)];
       // Optimization:  m[] is symmetric in its main diagonal,
       // so we could compute only one half and reflect it into the other.
@@ -353,21 +350,16 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
 	m[(i+1)*(d+1)+(j+1)] = _;
 	}
 
-      const double scalars[14] = {
+      constexpr double scalars[14] = {
         -1.,2.,-16.,288.,-9216.,460800.,-33177600.,3251404800.,-416179814400.,
         67421129932800.,-13484225986560000.,3263182688747520000.,
 	-939796614359285760000.,317651255653438586880000. };
-      double scalar;
-      if (d-1 <= 13)
-        scalar = scalars[d-1];
-      else
-	{
+      constexpr auto scalar = (d-1 <= 13) ? scalars[d-1] :
         // Should really compute Sloane's sequence A055546: (-1)^(j+1) / (2^j * (j!)^2).
 	// But approximation is okay (up to float overflow or underflow):
 	// relative, not absolute, volumes is what matters
 	// since we'll normalize them anyways (to sum to 1).
-	scalar = 1e30; // vague approximation
-	}
+	1e30; // vague approximation
 
       h.facetvolume[ifacet] = sqrt(determinant(m, d+1) / scalar);
       break;
@@ -377,7 +369,8 @@ bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid,
       printf("internal error in computing facet volumes.\n");
 #endif
     }
-  return true;
+  h.s = &s;
+  return h;
 }
 
 // Fill w[0 to d+1] with q's barycentric coords w.r.t. s (i.e., w.r.t. h).

bary.h

@@ -1,5 +1,5 @@
 #include "si.h"
 
-bool precomputeBary(const d_simplex& s, simplexHint& h, const vertex& centroid, const vertex* rgv, const vertex* raysCentroid, bool fRaySimplex = false);
+simplexHint precomputeBary(const d_simplex& s, const vertex& centroid, const std::vector<vertex>& rgv, const vertex* raysCentroid, bool fRaySimplex);
 
 bool computeBary(const simplexHint& h, const vertex& q, double* w, bool fRaySimplex = false);

callhull.c++

@@ -7,7 +7,6 @@
 #include <cmath>
 #include <numeric>
 #include <random>
-#include <vector>
 
 #include "hull.h"
 #include "si.h"
@@ -50,7 +49,7 @@ site read_next_site(long j) {
     return nullptr; // end of list
   hull_p = new_site(hull_p,j);
   for (auto i=0; i<dim; ++i) {
-    extern vertex* qi;
+    extern std::vector<vertex> qi;
     hull_p[i] = floor(mult_up * qi[j][i] + 0.5);    // these are the input points.
     if (hull_p[i] < mins[i])
       mins[i] = hull_p[i];
@@ -106,15 +105,13 @@ void resetEverything()
   mult_up = 1.0; // Already scaled by si.c++'s constexpr auto scale = 1e6.
 }
 
-// count stores how many true simplices are read (no -1 vertex).
-// countAll stores that, plus how many ray-simplices are read.
-// True simplices precede ray-simplices in the returned array.
-// Caller is responsible for delete[]ing return value.
-d_simplex* delaunay_tri(int dimArg, int cPt, int& count, int& countAll) {
+void delaunay_tri(std::vector<d_simplex>& si, std::vector<d_simplex>& siRay, int dimArg, int cPt) {
   resetEverything();
   if (dimArg > MAXDIM) {
     std::cerr << "dimension bound MAXDIM exceeded\n";
-    return nullptr;
+    si.clear();
+    siRay.clear();
+    return;
   }
   dim = dimArg;
   point_size = site_size = dim * sizeof(Coord);
@@ -141,26 +138,21 @@ d_simplex* delaunay_tri(int dimArg, int cPt, int& count, int& countAll) {
   make_output(root, visit_hull, facets_print, &CG_vlist_out, stdout);
   // CG_vlist_out wrote output to vpT,iT (simplices aka tets) and vpH,iH (ray-simplices aka hull, without the -1).
 
-  count = iT; // aka ctet, aka csi.
-  countAll = iT + iH; // aka ctet + ctri
-  d_simplex* sRet = new d_simplex[countAll];
-
-  // Stuff sRet like readSimplices().  First simplices, then ray-simplices.
-  for (auto i=0; i<countAll; ++i) {
-    d_simplex& s = sRet[i];
-    if (i < count) {
-      const auto& first = vpT[i].begin();
-      std::copy(first, first + d+1, s.begin());
-    } else {
-      // vpH[][d] is left uninitialized by CG_vlist_out().  It's implicitly -1.
-      const auto& first = vpH[i-count].begin();
-      std::copy(first, first + d, s.begin());
-      s[d] = -1;
-    }
+  si.resize(iT); // aka ctet.
+  siRay.resize(iH); // aka ctri.
+  for (auto i=0; i<iT; ++i) {
+    const auto& first = vpT[i].begin();
+    std::copy(first, first + d+1, si[i].begin());
+  }
+  for (auto i=0; i<iH; ++i) {
+    auto& s = siRay[i];
+    // vpH[][d] is left uninitialized by CG_vlist_out().  It's implicitly -1.
+    const auto& first = vpH[i].begin();
+    std::copy(first, first + d, s.begin());
+    s[d] = -1;
   }
 
   free_hull_storage();
   for (auto i=0; i<num_blocks; ++i)
     free(site_blocks[i]);
-  return sRet;
 }

edahiro.c++

@@ -560,18 +560,18 @@ inline int Tweak(const int ir)
 // Zero-based, not 1-based.
 // Permuting of points so they are sorted by y-coordinate.
 
-bool Edahiro_Init(int cpt, const vertex* qi, int csi, const d_simplex* si)
+bool Edahiro_Init(const std::vector<vertex>& qi, const std::vector<d_simplex>& si)
 {
   // Read the output of Hull algorithm.
   // nv, x's and y's, ntri, itri's.
   // Stuff nv and rgvertex, nr and rgregion.
 
-  nv = cpt;
+  nv = qi.size();
   if (nv < 3) {
     fprintf(stderr, "error: Edahiro_Init() needs at least 3 vertices, not %d.\n", nv);
     return false;
   }
-  nr = csi;
+  nr = si.size();
   if (nr < 1) {
     fprintf(stderr, "error: Edahiro_Init() needs at least 1 region, not %d.\n", nr);
     return false;

edahiro.h

@@ -1,3 +1,4 @@
 #include "si.h"
-bool Edahiro_Init(const int cpt, const vertex* qi, const int csi, const d_simplex* si);
+
+bool Edahiro_Init(const std::vector<vertex>&, const std::vector<d_simplex>&);
 int Edahiro_RegionFromPoint(double x, double y);

ga.c++

@@ -98,6 +98,10 @@ int RankAndCalculateFitness(int cPop, int cTopMost)
   for(i = 0; i < cPop; i++)
     {
     suitability[i] = ComputeSuitability(PvFromI(i));
+    if (std::isnan(suitability[i])) {
+      printf("Corrupt suitability.\n");
+      return -1;
+    }
     // optimization: cache retval of ComputeSuitability with a checksum on arg using cbMemberArg.
     if (suitability[i] == zFitnessMax)
       {

gacli.c++

@@ -4,6 +4,7 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdlib> // random()
+#include <iostream>
 #include <limits>
 
 #undef VERBOSE
@@ -22,10 +23,7 @@ static int cdimDst = -1;
 
 //;; iloop jloop rg[triangularNumber(cpt,i,j)]  -->  iloop rg[i]
 
-static double* rgzDistSrc = nullptr;
-static double* rgzDistDst = nullptr;
-
-inline double sq(double _) { return _*_; }
+template <class T> T sq(const T& _) { return _*_; }
 
 /*
 Each linear dimension of the space may have a vastly different
@@ -38,7 +36,7 @@ Later: log-scale distance option for some dimensions. (use log(x) internally
        in the linear geometry; use x for i/o.
 */
 
-void InitDistanceMatrixZ(int cpt, int cdimSrc, double* rgzDist, double* rgzPt)
+void InitDistanceMatrixZ(int cpt, int cdimSrc, double* rgzDist, const double* rgzPt)
 {
   // compute scaling factors for each dimension
   double rgzScale[cdimSrc];
@@ -58,11 +56,15 @@ void InitDistanceMatrixZ(int cpt, int cdimSrc, double* rgzDist, double* rgzPt)
     auto zSum = 0.0;
     for (auto idim = 0; idim < cdimSrc; ++idim) {
       zSum += sq((rgzPt[i * cdimSrc + idim] -
-	 rgzPt[j * cdimSrc + idim]) * rgzScale[idim]);
+	          rgzPt[j * cdimSrc + idim]) * rgzScale[idim]);
     }
     zDistMax = std::max(zDistMax, rgzDist[TriIJ(i,j,cpt)] = sqrt(zSum));
   }
-  // normalize distances wrt longest distance.
+  if (zDistMax == 0.0) {
+    std::cerr << "InitDistanceMatrixZ() got inputs that made only zero distances.  Crash imminent.\n";
+    // rgzPt[] might have been all zeros.  Just DBZ and get it over with.
+  }
+  // Normalize distances wrt longest distance.
   for (auto i=0; i<cpt-1; ++i)
   for (auto j=i+1; j<cpt; ++j)
     rgzDist[TriIJ(i,j,cpt)] /= zDistMax;
@@ -90,6 +92,14 @@ void InitDistanceMatrixL(int cpt, int cdimDst, double* rgzDist, short* rgzPt)
 // Compute RMS error between 2 distance vectors (possibly triangular matrices).
 inline double DDistanceMatrix(double* rgzDist0, double* rgzDist1, int cpt)
 {
+#ifdef VERBOSE
+  for (auto k=0; k<cpt; ++k) {
+    if (std::isnan(rgzDist0[k])) {
+      std::cerr << "DDistanceMatrix got corrupt input.";
+      return std::numeric_limits<double>::signaling_NaN();
+    }
+  }
+#endif
   auto z = 0.0;
   for (auto k=0; k<cpt; ++k)
     z += sq(rgzDist0[k] - rgzDist1[k]);
@@ -153,6 +163,9 @@ void MutateRandom(void* pv, long cIter)
   Tweak(pv);
 }
 
+static double* rgzDistSrc = nullptr;
+static double* rgzDistDst = nullptr;
+
 double ComputeSuitability(void* pv)
 {
   InitDistanceMatrixL(cpt, cdimDst, rgzDistDst, ((Member*)pv)->rgl);