Merge pull request #77 from mrcslws/grid-uniqueness-hypercubes

Compute bin sidelength faster, handle interrupts

src/nupic/experimental/GridUniqueness.cpp

@@ -53,6 +53,68 @@ namespace bg = boost::geometry;
 BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian)
 
 static std::atomic<bool> g_quitting(false);
+static std::atomic<bool> g_computeBinSidelengthShouldContinue(true);
+
+static std::atomic<int> g_computeUniqueHypercubeCounter(0);
+static std::atomic<int> g_computeBinsSideLengthCounter(0);
+static std::atomic<int> g_captureInterruptsCounter(0);
+
+static void (*g_prevHandler)(int) = nullptr;
+
+
+// Custom interrupt processing is particularly necessary in Jupyter notebooks.
+void processInterrupt(int sig)
+{
+  if (g_computeUniqueHypercubeCounter > 0)
+  {
+    g_quitting = true;
+  }
+
+  if (g_computeBinsSideLengthCounter > 0)
+  {
+    g_computeBinSidelengthShouldContinue = false;
+  }
+}
+
+class CaptureInterruptsRAII
+{
+public:
+  CaptureInterruptsRAII()
+  {
+    if (g_captureInterruptsCounter++ == 0)
+    {
+      g_prevHandler = signal(SIGINT, processInterrupt);
+    }
+  }
+
+  ~CaptureInterruptsRAII()
+  {
+    if (--g_captureInterruptsCounter == 0)
+    {
+      signal(SIGINT, g_prevHandler);
+      g_prevHandler = nullptr;
+    }
+  }
+};
+
+class IncrementRAII
+{
+public:
+  IncrementRAII(std::atomic<int>* counter)
+    : counter_(counter)
+  {
+    ++(*counter_);
+  }
+
+  ~IncrementRAII()
+  {
+    --(*counter_);
+  }
+
+private:
+  std::atomic<int>* counter_;
+};
+
 
 template<typename T>
 struct SquareMatrix2D {
@@ -1345,14 +1407,13 @@ nupic::experimental::grid_uniqueness::computeGridUniquenessHypercube(
 {
   typedef std::chrono::steady_clock Clock;
 
-  // Manually handle interrupts so that they're handled when running in a
-  // Jupyter notebook, and to make the threads return cleanly.
-  struct sigaction sigIntHandler;
-  sigIntHandler.sa_handler =
-    [](int s) { g_quitting = true; };
-  sigemptyset(&sigIntHandler.sa_mask);
-  sigIntHandler.sa_flags = 0;
-  sigaction(SIGINT, &sigIntHandler, NULL);
+  if (g_computeUniqueHypercubeCounter == 0)
+  {
+    // Recover from any previous interrupts.
+    g_quitting = false;
+  }
+  IncrementRAII incrementCounter(&g_computeUniqueHypercubeCounter);
+  CaptureInterruptsRAII captureInterrupts;
 
   NTA_CHECK(domainToPlaneByModule.size() == latticeBasisByModule.size())
     << "The two arrays of matrices must be the same length (one per module) "
@@ -1565,9 +1626,6 @@ nupic::experimental::grid_uniqueness::computeGridUniquenessHypercube(
 
   if (g_quitting)
   {
-    // The process might not be ending, the caller (e.g. the Python shell) is
-    // likely to catch this exception and continue, so prepare to run again.
-    g_quitting = false;
     NTA_THROW << "Caught interrupt signal";
   }
 
@@ -1747,12 +1805,12 @@ bool findGridCodeZero_noModulo(
   const vector<double>& x0,
   const vector<double>& dims,
   double readoutResolution,
+  std::atomic<bool>& shouldContinue,
   vector<double>* pointWithGridCodeZero = nullptr)
 {
   // Avoid doing any allocations in each recursion.
   vector<double> x0Copy(x0);
   vector<double> dimsCopy(dims);
-  std::atomic<bool> shouldContinue(true);
 
   vector<double> defaultPointBuffer;
 
@@ -1790,27 +1848,41 @@ bool findGridCodeZero_noModulo(
 bool findGridCodeZeroAtRadius(
   double radius,
   const vector<vector<vector<double>>>& domainToPlaneByModule,
-  double readoutResolution)
+  double readoutResolution,
+  std::atomic<bool>& shouldContinue)
 {
   const size_t numDims = domainToPlaneByModule[0][0].size();
 
   for (size_t iDim = 0; iDim < numDims; ++iDim)
   {
     // Test the hyperplanes formed by setting this dimension to r and -r.
-
     vector<double> x0(numDims, -radius);
     vector<double> dims(numDims, 2*radius);
 
+    // Optimization: for the final dimension, don't go negative. Half of the
+    // hypercube will be equal-and-opposite phases of the other half, so we
+    // ignore the lower half of the final dimension.
+    x0[numDims - 1] = 0;
+    dims[numDims - 1] = radius;
+
     dims[iDim] = 0;
-    if (findGridCodeZero_noModulo(domainToPlaneByModule,
-                                  x0, dims, readoutResolution))
+
+    // Test -r
+    if (iDim != numDims - 1)
     {
-      return true;
+      if (findGridCodeZero_noModulo(domainToPlaneByModule,
+                                    x0, dims, readoutResolution,
+                                    shouldContinue))
+      {
+        return true;
+      }
     }
 
+    // Test +r
     x0[iDim] = radius;
     if (findGridCodeZero_noModulo(domainToPlaneByModule,
-                                  x0, dims, readoutResolution))
+                                  x0, dims, readoutResolution,
+                                  shouldContinue))
     {
       return true;
     }
@@ -1826,15 +1898,24 @@ nupic::experimental::grid_uniqueness::computeBinSidelength(
   double resultPrecision,
   double upperBound)
 {
+  if (g_computeBinsSideLengthCounter == 0)
+  {
+    // Recover from any previous interrupts.
+    g_computeBinSidelengthShouldContinue = true;
+  }
+  IncrementRAII incrementCounter(&g_computeBinsSideLengthCounter);
+  CaptureInterruptsRAII captureInterrupts;
+
   double tested = 0;
   double radius = 0.5;
 
   while (findGridCodeZeroAtRadius(radius,
                                   domainToPlaneByModule,
-                                  readoutResolution))
+                                  readoutResolution,
+                                  g_computeBinSidelengthShouldContinue))
   {
     tested = radius;
-    radius += 0.5;
+    radius *= 2;
 
     if (radius > upperBound)
     {
@@ -1849,19 +1930,25 @@ nupic::experimental::grid_uniqueness::computeBinSidelength(
   double dec = (radius - tested) / 2;
 
   // The possible error is equal to dec*2.
-  while (dec*2 > resultPrecision2)
+  while (g_computeBinSidelengthShouldContinue && dec*2 > resultPrecision2)
   {
     const double testRadius = radius - dec;
 
     if (!findGridCodeZeroAtRadius(testRadius,
                                   domainToPlaneByModule,
-                                  readoutResolution))
+                                  readoutResolution,
+                                  g_computeBinSidelengthShouldContinue))
     {
       radius = testRadius;
     }
 
     dec /= 2;
   }
 
+  if (!g_computeBinSidelengthShouldContinue)
+  {
+    NTA_THROW << "Caught interrupt signal";
+  }
+
   return 2 * radius;
 }