Add computeBinRectangle algorithm

src/nupic/bindings/experimental.i

@@ -260,6 +260,55 @@ using namespace nupic;
   }
 }
 
+%pythoncode %{
+  def computeBinRectangle(domainToPlaneByModule, phaseResolution,
+                          resultPrecision, upperBound=1000.0, timeout=-1.0):
+    domainToPlaneByModule = numpy.asarray(domainToPlaneByModule, dtype="float64")
+
+    return _computeBinRectangle(
+      domainToPlaneByModule, phaseResolution, resultPrecision, upperBound, timeout)
+%}
+
+%inline {
+  PyObject* _computeBinRectangle(PyObject* py_domainToPlaneByModule,
+                                 Real64 readoutResolution,
+                                 Real64 resultPrecision,
+                                 Real64 upperBound,
+                                 Real64 timeout)
+  {
+    PyArrayObject* pyArr_domainToPlaneByModule =
+      (PyArrayObject*)py_domainToPlaneByModule;
+    NTA_CHECK(PyArray_NDIM(pyArr_domainToPlaneByModule) == 3);
+    npy_intp* npy_dims = PyArray_DIMS(pyArr_domainToPlaneByModule);
+
+    std::vector<std::vector<std::vector<Real64>>> domainToPlaneByModule;
+    for (size_t i = 0; i < npy_dims[0]; i++)
+    {
+      std::vector<std::vector<Real64>> module;
+      for (size_t j = 0; j < npy_dims[1]; j++)
+      {
+        std::vector<Real64> row;
+        for (size_t k = 0; k < npy_dims[2]; k++)
+        {
+          row.push_back(*(Real64*)PyArray_GETPTR3(pyArr_domainToPlaneByModule,
+                                                  i, j, k));
+        }
+        module.push_back(row);
+      }
+      domainToPlaneByModule.push_back(module);
+    }
+
+    std::vector<Real64> result =
+      nupic::experimental::grid_uniqueness::computeBinRectangle(
+        domainToPlaneByModule, readoutResolution, resultPrecision,
+        upperBound, timeout);
+
+    return nupic::NumpyVectorT<Real64>(result.size(),
+                                       result.data()).forPython();
+  }
+}
+
+
 
 #endif NTA_OS_WINDOWS
 

src/nupic/experimental/GridUniqueness.cpp

@@ -2105,3 +2105,166 @@ nupic::experimental::grid_uniqueness::computeBinSidelength(
       return 2*radius;
   }
 }
+
+vector<double>
+nupic::experimental::grid_uniqueness::computeBinRectangle(
+  const vector<vector<vector<double>>>& domainToPlaneByModule,
+  double readoutResolution,
+  double resultPrecision,
+  double upperBound,
+  double timeout)
+{
+  typedef std::chrono::steady_clock Clock;
+
+  enum ExitReason {
+    Timeout,
+    Interrupt,
+    Completed
+  };
+
+  std::atomic<ExitReason> exitReason(ExitReason::Completed);
+
+  ThreadSafeQueue<Message> messages;
+  CaptureInterruptsRAII captureInterrupts(&messages);
+
+  std::atomic<bool> shouldContinue(true);
+  std::atomic<bool> interrupted(false);
+  std::thread messageThread(
+    [&]() {
+      while (true)
+      {
+        switch (messages.take())
+        {
+          case Message::Interrupt:
+            shouldContinue = false;
+            exitReason = ExitReason::Interrupt;
+            interrupted = true;
+            break;
+          case Message::Timeout:
+            shouldContinue = false;
+            exitReason = ExitReason::Timeout;
+            break;
+          case Message::Exiting:
+            return;
+        }
+      }
+    });
+
+  ScheduledTask* scheduledTask = nullptr;
+  if (timeout > 0)
+  {
+    scheduledTask = new ScheduledTask(
+      Clock::now() + std::chrono::duration<double>(timeout),
+      [&messages](){ messages.put(Message::Timeout); });
+  }
+
+  const size_t numDims = domainToPlaneByModule[0][0].size();
+
+  double radius = 0.5;
+
+  while (findGridCodeZeroAtRadius(radius,
+                                  domainToPlaneByModule,
+                                  readoutResolution,
+                                  shouldContinue))
+  {
+    radius *= 2;
+
+    if (radius > upperBound)
+    {
+      return {};
+    }
+  }
+
+  // The radius needs to be twice as precise to get the sidelength sufficiently
+  // precise.
+  const double resultPrecision2 = resultPrecision / 2;
+
+  vector<double> radii(numDims, radius);
+
+  for (size_t iDim = 0; iDim < numDims; ++iDim)
+  {
+    double dec = radius / 2;
+
+    // The possible error is equal to dec*2.
+    while (shouldContinue && dec*2 > resultPrecision2)
+    {
+      const double testRadius = radii[iDim] - dec;
+
+      vector<double> x0(numDims);
+      vector<double> dims(numDims);
+
+      for (size_t iDim2 = 0; iDim2 < numDims; ++iDim2)
+      {
+        if (iDim2 == numDims - 1)
+        {
+          // 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[iDim2] = 0;
+          dims[iDim2] = radii[iDim2];
+        }
+        else
+        {
+          x0[iDim2] = -radii[iDim2];
+          dims[iDim2] = 2*radii[iDim2];
+        }
+      }
+
+      dims[iDim] = 0;
+
+      bool foundZero = false;
+      if (iDim != numDims - 1)
+      {
+        // Test -r
+        x0[iDim] = -testRadius;
+        foundZero = findGridCodeZero_noModulo(domainToPlaneByModule,
+                                              x0, dims, readoutResolution,
+                                              shouldContinue);
+      }
+
+      if (!foundZero)
+      {
+        // Test r
+        x0[iDim] = testRadius;;
+        foundZero = findGridCodeZero_noModulo(domainToPlaneByModule,
+                                              x0, dims, readoutResolution,
+                                              shouldContinue);;
+      }
+
+      if (!foundZero)
+      {
+        radii[iDim] = testRadius;
+      }
+      dec /= 2;
+    }
+  }
+
+  if (scheduledTask != nullptr)
+  {
+    delete scheduledTask;
+    scheduledTask = nullptr;
+  }
+
+  messages.put(Message::Exiting);
+  messageThread.join();
+
+  switch (exitReason.load())
+  {
+    case ExitReason::Timeout:
+      // Python code may check for the precise string "timeout".
+      NTA_THROW << "timeout";
+    case ExitReason::Interrupt:
+      NTA_THROW << "interrupt";
+    case ExitReason::Completed:
+    default:
+    {
+      vector<double> sidelengths(numDims);
+      for (size_t iDim = 0; iDim < numDims; ++iDim)
+      {
+        sidelengths[iDim] = radii[iDim]*2;
+      }
+
+      return sidelengths;
+    }
+  }
+}

src/nupic/experimental/GridUniqueness.hpp

@@ -173,6 +173,48 @@ namespace nupic {
         Real64 resultPrecision,
         Real64 upperBound = 1000.0,
         Real64 timeout = -1.0);
+
+      /**
+       * Like computeBinSidelength, but it computes a hyperrectangle rather than
+       * a hypercube.
+       *
+       * @param domainToPlaneByModule
+       * A list of 2*k matrices, one per module. The matrix converts from a
+       * point in the domain to a point on a plane, normalizing for grid cell
+       * scale.
+       *
+       * @param readoutResolution
+       * The precision of readout of this grid code, measured in distance on the
+       * plane. For example, if this is 0.2, then all points on the plane are
+       * indistinguishable from those in their surrounding +- 0.1 range.
+       *
+       * @param resultPrecision
+       * The precision level for this sidelength. This algorithm will
+       * binary-search for the smallest hyperrectangle it can place around zero,
+       * stopping at the point when each side is within 'resultPrecision' of the
+       * actual smallest hyperrectangle.
+       *
+       * @param upperBound
+       * Instructs the algorithm when it should give up. If the provided matrix
+       * never moves away from grid code zero, the algorithm could search
+       * forever. This parameter tells it when it should stop searching.
+       *
+       * @param timeout
+       * Specifies how long to try. This function will give up after 'timeout'
+       * seconds. If <= 0, the function will not time out. On timeout, the
+       * function throws an exception with message "timeout". In Python this
+       * exception is of type RuntimeError.
+       *
+       * @return
+       * The dimensions of this hyperrectangle. Returns an empty vector if a
+       * surface can't be found (i.e. if upperBound is reached.)
+       */
+      std::vector<Real64> computeBinRectangle(
+        const std::vector<std::vector<std::vector<Real64>>>& domainToPlaneByModule,
+        Real64 readoutResolution,
+        Real64 resultPrecision,
+        Real64 upperBound = 1000.0,
+        Real64 timeout = -1.0);
     }
   }
 }