From 81296a9e1291bc99b73b4f968b4dfa9d0a027584 Mon Sep 17 00:00:00 2001 From: Marcus Lewis Date: Mon, 22 Jul 2019 13:55:06 -0700 Subject: [PATCH] Unbundle GridUniqueness code into gridcodingrange library --- src/CMakeLists.txt | 15 - src/nupic/bindings/experimental.i | 187 -- src/nupic/experimental/GridUniqueness.cpp | 2471 ----------------- src/nupic/experimental/GridUniqueness.hpp | 227 -- .../unit/experimental/GridUniquenessTest.cpp | 766 ----- 5 files changed, 3666 deletions(-) delete mode 100644 src/nupic/experimental/GridUniqueness.cpp delete mode 100644 src/nupic/experimental/GridUniqueness.hpp delete mode 100644 src/test/unit/experimental/GridUniquenessTest.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3f1c2c16..b841cb54 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -239,14 +239,6 @@ set(src_htmresearchcore_srcs nupic/experimental/SDRSelection.cpp ) -if(NOT MINGW) - # This file uses threading that's not available in our version of MINGW. - set(src_htmresearchcore_srcs - ${src_htmresearchcore_srcs} - nupic/experimental/GridUniqueness.cpp - ) -endif() - set(src_lib_static_htmresearchcore_srcs ${src_capnp_generated_srcs} ${src_htmresearchcore_srcs}) @@ -313,13 +305,6 @@ set(src_htmresearch_core_gtest_srcs test/unit/UnitTestMain.cpp test/unit/utils/GroupByTest.cpp ) -if(NOT MINGW) - # This file uses threading that's not available in our version of MINGW. - set(src_htmresearch_core_gtest_srcs - ${src_htmresearch_core_gtest_srcs} - test/unit/experimental/GridUniquenessTest.cpp - ) -endif() add_executable(${src_executable_gtests} ${src_htmresearch_core_gtest_srcs} ) diff --git a/src/nupic/bindings/experimental.i b/src/nupic/bindings/experimental.i index a8f728a5..dd31e564 100644 --- a/src/nupic/bindings/experimental.i +++ b/src/nupic/bindings/experimental.i @@ -127,193 +127,6 @@ using namespace nupic; %naturalvar; -#ifndef NTA_OS_WINDOWS - -// computeGridUniquenessHypercube uses threading that's not available in our -// MINGW Windows build system. - -%{ - #include -%} - -%pythoncode %{ - def computeGridUniquenessHypercube(domainToPlaneByModule, latticeBasisByModule, - phaseResolution, ignoredCenterDiameter, - pingInterval=10.0): - domainToPlaneByModule = numpy.asarray(domainToPlaneByModule, dtype="float64") - latticeBasisByModule = numpy.asarray(latticeBasisByModule, dtype="float64") - - return _computeGridUniquenessHypercube( - domainToPlaneByModule, latticeBasisByModule, phaseResolution, - ignoredCenterDiameter, pingInterval) -%} - -%inline { - PyObject* _computeGridUniquenessHypercube(PyObject* py_domainToPlaneByModule, - PyObject* py_latticeBasisByModule, - Real64 phaseResolution, - Real64 ignoredCenterDiameter, - Real64 pingInterval) - { - PyArrayObject* pyArr_domainToPlaneByModule = - (PyArrayObject*)py_domainToPlaneByModule; - NTA_CHECK(PyArray_NDIM(pyArr_domainToPlaneByModule) == 3); - npy_intp* npy_dims = PyArray_DIMS(pyArr_domainToPlaneByModule); - - std::vector>> domainToPlaneByModule; - for (size_t i = 0; i < npy_dims[0]; i++) - { - std::vector> module; - for (size_t j = 0; j < npy_dims[1]; j++) - { - std::vector 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); - } - - PyArrayObject* pyArr_latticeBasisByModule = - (PyArrayObject*)py_latticeBasisByModule; - NTA_CHECK(PyArray_NDIM(pyArr_latticeBasisByModule) == 3); - npy_dims = PyArray_DIMS(pyArr_latticeBasisByModule); - - std::vector>> latticeBasisByModule; - for (size_t i = 0; i < npy_dims[0]; i++) - { - std::vector> module; - for (size_t j = 0; j < npy_dims[1]; j++) - { - std::vector row; - for (size_t k = 0; k < npy_dims[2]; k++) - { - row.push_back(*(Real64*)PyArray_GETPTR3(pyArr_latticeBasisByModule, - i, j, k)); - } - module.push_back(row); - } - latticeBasisByModule.push_back(module); - } - - std::pair> result = - nupic::experimental::grid_uniqueness::computeGridUniquenessHypercube( - domainToPlaneByModule, latticeBasisByModule, phaseResolution, - ignoredCenterDiameter, pingInterval); - - PyObject* pyResult = PyTuple_New(2); - PyTuple_SetItem(pyResult, 0, PyFloat_FromDouble(result.first)); - PyTuple_SetItem(pyResult, 1, nupic::NumpyVectorT(result.second.size(), - result.second.data()) - .forPython()); - - return pyResult; - } -} - - -%pythoncode %{ - def computeBinSidelength(domainToPlaneByModule, phaseResolution, - resultPrecision, upperBound=1000.0, timeout=-1.0): - domainToPlaneByModule = numpy.asarray(domainToPlaneByModule, dtype="float64") - - return _computeBinSidelength( - domainToPlaneByModule, phaseResolution, resultPrecision, upperBound, timeout) -%} - -%inline { - PyObject* _computeBinSidelength(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>> domainToPlaneByModule; - for (size_t i = 0; i < npy_dims[0]; i++) - { - std::vector> module; - for (size_t j = 0; j < npy_dims[1]; j++) - { - std::vector 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); - } - - Real64 result = - nupic::experimental::grid_uniqueness::computeBinSidelength( - domainToPlaneByModule, readoutResolution, resultPrecision, - upperBound, timeout); - - return PyFloat_FromDouble(result); - } -} - -%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>> domainToPlaneByModule; - for (size_t i = 0; i < npy_dims[0]; i++) - { - std::vector> module; - for (size_t j = 0; j < npy_dims[1]; j++) - { - std::vector 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 result = - nupic::experimental::grid_uniqueness::computeBinRectangle( - domainToPlaneByModule, readoutResolution, resultPrecision, - upperBound, timeout); - - return nupic::NumpyVectorT(result.size(), - result.data()).forPython(); - } -} - - - -#endif NTA_OS_WINDOWS - //-------------------------------------------------------------------------------- // Apical Tiebreak Temporal Memory //-------------------------------------------------------------------------------- diff --git a/src/nupic/experimental/GridUniqueness.cpp b/src/nupic/experimental/GridUniqueness.cpp deleted file mode 100644 index 27737177..00000000 --- a/src/nupic/experimental/GridUniqueness.cpp +++ /dev/null @@ -1,2471 +0,0 @@ -/* --------------------------------------------------------------------- - * Numenta Platform for Intelligent Computing (NuPIC) - * Copyright (C) 2018, Numenta, Inc. Unless you have an agreement - * with Numenta, Inc., for a separate license for this software code, the - * following terms and conditions apply: - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Affero Public License version 3 as - * published by the Free Software Foundation. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU Affero Public License for more details. - * - * You should have received a copy of the GNU Affero Public License - * along with this program. If not, see http://www.gnu.org/licenses. - * - * http://numenta.org/licenses/ - * ---------------------------------------------------------------------- - */ - -/** @file - * Implementations for GridUniqueness.hpp - */ - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -using std::vector; -using std::pair; -namespace bg = boost::geometry; - -BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian) - - -template -class ThreadSafeQueue { -public: - void put(T v) - { - std::unique_lock lock(mutex_); - queue_.push(v); - cv_.notify_one(); - } - - T take() - { - std::unique_lock lock(mutex_); - while (queue_.empty()) - { - cv_.wait(lock); - } - - T ret = queue_.front(); - queue_.pop(); - return ret; - } - -private: - std::mutex mutex_; - std::condition_variable cv_; - std::queue queue_; -}; - - -class ScheduledTask { -public: - template - ScheduledTask(T timeout, F f) - { - timerThread_ = std::thread( - [this, timeout, f]() { this->waiterThread_(timeout, f); }); - } - - ~ScheduledTask() - { - { - std::unique_lock lock(timerMutex_); - timerCancelled_ = true; - timerCondition_.notify_one(); - } - timerThread_.join(); - } - -private: - - template - void waiterThread_(T tTimeout, F f) - { - std::unique_lock lock(this->timerMutex_); - - while (!this->timerCancelled_) - { - if (this->timerCondition_.wait_until(lock, tTimeout) == - std::cv_status::timeout) - { - f(); - this->timerCancelled_ = true; - } - } - } - - std::mutex timerMutex_; - std::condition_variable timerCondition_; - bool timerCancelled_ = false; - std::thread timerThread_; -}; - - -enum Message { - Interrupt, - Timeout, - Exiting -}; - - -static size_t g_captureInterruptsCounter = 0; -static std::mutex g_messageQueuesMutex; -static vector*> g_messageQueues; -static void (*g_prevHandler)(int) = nullptr; - - -// Custom interrupt processing is particularly necessary in Jupyter notebooks. -void processInterrupt(int sig) -{ - { - std::unique_lock lock(g_messageQueuesMutex); - for (ThreadSafeQueue* messageQueue : g_messageQueues) - { - messageQueue->put(Message::Interrupt); - } - } -} - -class CaptureInterruptsRAII -{ -public: - CaptureInterruptsRAII(ThreadSafeQueue* messages) - : messages_(messages) - { - std::unique_lock lock(g_messageQueuesMutex); - - if (g_captureInterruptsCounter++ == 0) - { - g_prevHandler = signal(SIGINT, processInterrupt); - } - - g_messageQueues.push_back(messages); - } - - ~CaptureInterruptsRAII() - { - std::unique_lock lock(g_messageQueuesMutex); - - if (--g_captureInterruptsCounter == 0) - { - signal(SIGINT, g_prevHandler); - g_prevHandler = nullptr; - } - - for (auto it = g_messageQueues.begin(); it != g_messageQueues.end(); it++) - { - if (*it == messages_) - { - g_messageQueues.erase(it); - break; - } - } - } - -private: - ThreadSafeQueue* messages_; -}; - -template -struct SquareMatrix2D { - T v00; - T v01; - T v10; - T v11; -}; - -pair transform2D(const SquareMatrix2D& M, - pair p) -{ - return {M.v00*p.first + M.v01*p.second, - M.v10*p.first + M.v11*p.second}; -} - -SquareMatrix2D invert2DMatrix(const vector>& M) -{ - const double detInv = 1 / (M[0][0]*M[1][1] - M[0][1]*M[1][0]); - return {detInv*M[1][1], -detInv*M[0][1], - -detInv*M[1][0], detInv*M[0][0]}; -} - -pair transformND(const vector>& M, - const double p[]) -{ - pair result = {0, 0}; - - for (size_t col = 0; col < M[0].size(); col++) - { - result.first += M[0][col]*p[col]; - } - - for (size_t col = 0; col < M[1].size(); col++) - { - result.second += M[1][col]*p[col]; - } - - return result; -} - - -struct LatticeBox { - double xmin; - double xmax; - pair middle; -}; - -/** - * Enumerate the points of a lattice near or within a specified rectangle. This - * is equivalent to checking whether any circles centered on the points of a - * lattice overlap the rectangle. - * - * This enumeration works by first guessing at an initial lattice point (i, j), - * then performing a series of nested sweeps. It sweeps downward, decrementing j - * until decrementing j causes the distance of the lattice point from the - * rectangle to increase. Then it sweeps upward using the same rules. It sweeps - * leftward, decrementing i and repeating these downward and upward sweeps at - * each i until it passes the lowest i possible for this rectangle. It repeats - * the process sweeping right. As it sweeps left and right it attempts to choose - * a good starting j value by considering results from the previous downward and - * upward sweeps. - */ -class LatticePointEnumerator -{ -public: - LatticePointEnumerator(const SquareMatrix2D& latticeBasis, - const SquareMatrix2D& inverseLatticeBasis, - const LatticeBox& cachedLatticeBox, - const pair& shift, - double left, double right, double bottom, double top, - double rSquared) - :latticeBasis_(latticeBasis), left_(left), right_(right), bottom_(bottom), - top_(top), rSquared_(rSquared) - { - const pair ijShift = transform2D(inverseLatticeBasis, shift); - iMin_ = ceil(cachedLatticeBox.xmin + ijShift.first); - iMax_ = floor(cachedLatticeBox.xmax + ijShift.first); - i_ = iStart_ = floor(cachedLatticeBox.middle.first + ijShift.first); - j_ = j0_ = jStart_ = floor(cachedLatticeBox.middle.second + ijShift.second); - } - - LatticePointEnumerator(const SquareMatrix2D& latticeBasis, - const SquareMatrix2D& inverseLatticeBasis, - double left, double right, double bottom, double top, - double r, double rSquared) - :latticeBasis_(latticeBasis), left_(left), right_(right), bottom_(bottom), - top_(top), rSquared_(rSquared) - { - // Find the bounding box of the rectangle in the lattice's basis. - double xmin = std::numeric_limits::max(); - double xmax = std::numeric_limits::lowest(); - - const double paddedLeft = left_ - r; - const double paddedRight = right_ + r; - const double paddedBottom = bottom_ - r; - const double paddedTop = top_ + r; - - pair ij; - ij = transform2D(inverseLatticeBasis, {paddedLeft, paddedBottom}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedRight, paddedBottom}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedLeft, paddedTop}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedRight, paddedTop}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - - iMin_ = ceil(xmin); - iMax_ = floor(xmax); - - ij = transform2D(inverseLatticeBasis, {left, bottom}); - iStart_ = floor(ij.first); - jStart_ = floor(ij.second); - - i_ = iStart_; - j_ = j0_ = jStart_; - } - - bool getNext(pair *out) - { - bool foundContainedPoint = false; - - while (!foundContainedPoint && (sweepingLeft_ || i_ <= iMax_)) - { - const pair p = transform2D(latticeBasis_, {i_, j_}); - - const double nearestX = std::max(left_, - std::min(p.first, - right_)); - const double nearestY = std::max(bottom_, - std::min(p.second, - top_)); - - const double dSquared = (pow(p.first - nearestX, 2) + - pow(p.second - nearestY, 2)); - - if (j_ == j0_) - { - dSquared_j0_ = dSquared; - } - - if (dSquared < innerSweepMin_) - { - innerSweepMin_ = dSquared; - jForInnerSweepMin_ = j_; - } - - foundContainedPoint = (dSquared <= rSquared_); - if (foundContainedPoint) - { - *out = p; - } - - // If we're moving away from the box, end this part of the inner sweep. - const bool endInnerSweep = (!foundContainedPoint && - dSquared >= dSquaredPrev_); - - if (endInnerSweep) - { - if (sweepingDown_) - { - sweepingDown_ = false; - - j_ = j0_ + 1; - dSquaredPrev_ = dSquared_j0_; - } - else - { - if (i_ == iStart_) - { - jForInnerSweepMin_i0_ = jForInnerSweepMin_; - } - - if (sweepingLeft_) - { - if (--i_ < iMin_) - { - sweepingLeft_ = false; - i_ = iStart_ + 1; - - j0_ = jForInnerSweepMin_i0_; - j_ = j0_; - } - else - { - j0_ = jForInnerSweepMin_; - j_ = j0_; - } - } - else - { - ++i_; - j0_ = jForInnerSweepMin_; - j_ = j0_; - } - - sweepingDown_ = true; - innerSweepMin_ = std::numeric_limits::max(); - dSquaredPrev_ = std::numeric_limits::max(); - } - } - else - { - if (sweepingDown_) - { - --j_; - } - else - { - ++j_; - } - - dSquaredPrev_ = dSquared; - } - } - - return foundContainedPoint; - } - - void restart() - { - dSquaredPrev_ = std::numeric_limits::max(); - innerSweepMin_ = std::numeric_limits::max(); - i_ = iStart_; - j_ = j0_ = jStart_; - finished_ = false; - sweepingLeft_ = true; - sweepingDown_ = true; - } - -private: - - const SquareMatrix2D& latticeBasis_; - const double left_; - const double right_; - const double bottom_; - const double top_; - const double rSquared_; - - double dSquaredPrev_ = std::numeric_limits::max(); - double innerSweepMin_ = std::numeric_limits::max(); - long long jForInnerSweepMin_ = std::numeric_limits::lowest(); - long long jForInnerSweepMin_i0_ = std::numeric_limits::lowest(); - double dSquared_j0_; - - bool finished_ = false; - bool sweepingLeft_ = true; - bool sweepingDown_ = true; - - long long iStart_; - long long jStart_; - long long i_; - long long j_; - long long j0_; - long long iMin_; - long long iMax_; -}; - -/** - * Enumerate the vertices of a hyperrectangle by incrementing an integer and - * converting each binary representation into a vertex. - */ -class HyperrectangleVertexEnumerator -{ -public: - HyperrectangleVertexEnumerator(const double dims[], size_t numDims) - :dims_(dims), numDims_(numDims), upper_(pow(2, numDims)), bitvector_(0x0) - { - } - - bool getNext(double *out) - { - if (bitvector_ >= upper_) - { - return false; - } - - for (size_t bit = 0; bit < numDims_; bit++) - { - out[bit] = (bitvector_ & (0x1 << bit)) - ? dims_[bit] - : 0; - } - - bitvector_++; - - return true; - } - - void restart() - { - bitvector_ = 0x0; - } - -private: - const double *dims_; - const size_t numDims_; - const double upper_; - unsigned int bitvector_; -}; - -/** - * Compute d % 1.0, returning a value within the range [-0.5, 0.5] - */ -double mod1_05(double d) -{ - double ret = d - floor(d); - if (ret > 0.5) - { - ret -= 1.0; - } - - return ret; -} - -/** - * Quickly check a few points in this hyperrectangle to see if they have grid - * code zero. - */ -bool tryFindGridCodeZero( - const vector>>& domainToPlaneByModule, - const vector>& latticeBasisByModule, - const vector>& inverseLatticeBasisByModule, - size_t numDims, - const double x0[], - const double dims[], - double rSquared, - double vertexBuffer[]) -{ - for (size_t iDim = 0; iDim < numDims; iDim++) - { - vertexBuffer[iDim] = x0[iDim] + (dims[iDim]/2); - } - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const pair pointOnPlane = - transformND(domainToPlaneByModule[iModule], vertexBuffer); - - const pair pointOnUnrolledTorus = - transform2D(inverseLatticeBasisByModule[iModule], pointOnPlane); - - const pair pointOnTorus = { - mod1_05(pointOnUnrolledTorus.first), - mod1_05(pointOnUnrolledTorus.second) - }; - - const pair pointOnPlaneNearestZero = - transform2D(latticeBasisByModule[iModule], pointOnTorus); - - if (pow(pointOnPlaneNearestZero.first, 2) + - pow(pointOnPlaneNearestZero.second, 2) > rSquared) - { - return false; - } - } - - return true; -} - -double circleDistanceFromLineSegmentSquared(pair start, - pair end, - pair center, - pair unitVector, - double lineLength) -{ - const double centerDistanceAlongLine = - (unitVector.first*(center.first - start.first) + - unitVector.second*(center.second - start.second)); - - pair nearestPointOnLine; - if (centerDistanceAlongLine <= 0) - { - nearestPointOnLine = start; - } - else if (centerDistanceAlongLine < lineLength) - { - nearestPointOnLine = { - start.first + unitVector.first * centerDistanceAlongLine, - start.second + unitVector.second * centerDistanceAlongLine, - }; - } - else - { - nearestPointOnLine = end; - } - - return (pow(center.first - nearestPointOnLine.first, 2) + - pow(center.second - nearestPointOnLine.second, 2)); -} - -bool lineSegmentIntersectsCircle(pair start, - pair end, - pair center, - double rSquared) -{ - pair unitVector = {end.first - start.first, - end.second - start.second}; - const double lineLength = sqrt(pow(unitVector.first, 2) + - pow(unitVector.second, 2)); - unitVector.first /= lineLength; - unitVector.second /= lineLength; - - return circleDistanceFromLineSegmentSquared(start, end, center, unitVector, - lineLength) <= rSquared; -} - -struct LineInfo2D { - pair unitVector; - double length; -}; - -struct HalfPlaneData -{ - pair normalVector; - double top; -}; - -struct PolygonData { - bool isValidPolygon; - pair centroid; - vector> points; - vector thetaByPoint; - vector halfPlaneByEndpoint; - vector lineInfoByEndpoint; -}; - -/** - * Computing theta can be expensive. Get a number that increases monotonically - * with theta. - */ -double getThetaIndex(double dx, double dy) -{ - // // Simple slow version. - // double theta = atan2(dy, dx); - // if (theta < 0) - // { - // theta += 2*M_PI; - // } - // return theta; - - // Faster version. - const bool dxSign = dx >= 0; - const bool dySign = dy >= 0; - - if (dxSign && dySign) - { - if (dy < dx) - { - return dy / dx; - } - else - { - return 2.0 - (dx / dy); - } - } - else if (!dxSign && dySign) - { - dx = std::abs(dx); - if (dy > dx) - { - return 2.0 + (dx / dy); - } - else - { - return 4.0 - (dy / dx); - } - } - else if (!dxSign && !dySign) - { - dx = std::abs(dx); - dy = std::abs(dy); - if (dy < dx) - { - return 4.0 + (dy / dx); - } - else - { - return 6.0 - (dx / dy); - } - } - else - { - dy = std::abs(dy); - if (dy > dx) - { - return 6.0 + (dx / dy); - } - else - { - return 8.0 - (dy / dx); - } - } -} - -double pointDistanceFromShadowSquared( - pair point, - const PolygonData& shadow, - size_t numDims) -{ - const pair *p1; - const pair *p2; - const LineInfo2D *lineInfo; - if (shadow.isValidPolygon) - { - // Figure out which edge to check. - const double thetaIndex = getThetaIndex(point.first - shadow.centroid.first, - point.second - shadow.centroid.second); - - const vector::const_iterator it = (shadow.points.size() <= 8) - ? std::find_if(shadow.thetaByPoint.begin(), shadow.thetaByPoint.end(), - [&](double d) - { - return thetaIndex < d; - }) - : std::lower_bound(shadow.thetaByPoint.begin(), shadow.thetaByPoint.end(), - thetaIndex); - - // If we reached the end, set i back to 0. - const size_t iEnd = std::distance(shadow.thetaByPoint.begin(), - it) % shadow.thetaByPoint.size(); - - // Check whether the lattice point is contained within the shadow. - const pair& normalVector = - shadow.halfPlaneByEndpoint[iEnd].normalVector; - if (normalVector.first*point.first + - normalVector.second*point.second - <= shadow.halfPlaneByEndpoint[iEnd].top) - { - // The point is contained. - return 0.0; - } - - // Proceed to check whether the circle around the lattice point intersects - // with this edge. - p1 = (iEnd == 0) - ? &shadow.points[shadow.points.size() - 1] - : &shadow.points[iEnd - 1]; - p2 = &shadow.points[iEnd]; - lineInfo = &shadow.lineInfoByEndpoint[iEnd]; - } - else - { - p1 = &shadow.points[1]; - p2 = &shadow.points[0]; - lineInfo = &shadow.lineInfoByEndpoint[0]; - } - - return circleDistanceFromLineSegmentSquared(*p1, *p2, point, - lineInfo->unitVector, - lineInfo->length); -} - -vector> getShadowConvexHull( - const vector>& domainToPlane, - size_t numDims, - const double dims[], - double vertexBuffer[]) -{ - if (numDims == 2) - { - // Optimization: in 2D we already know the convex hull. - const double point1[2] = {0, 0}; - const double point2[2] = {0, dims[1]}; - const double point3[2] = {dims[0], dims[1]}; - const double point4[2] = {dims[0], 0}; - - const pair p1 = transformND(domainToPlane, point1); - - return {p1, - transformND(domainToPlane, point2), - transformND(domainToPlane, point3), - transformND(domainToPlane, point4), - p1}; - } - - typedef boost::tuple point; - typedef bg::model::polygon polygon; - - polygon poly; - - HyperrectangleVertexEnumerator vertices(dims, numDims); - while (vertices.getNext(vertexBuffer)) - { - const pair p = transformND(domainToPlane, vertexBuffer); - bg::append(poly, point(p.first, p.second)); - } - - polygon hull; - bg::convex_hull(poly, hull); - - vector> convexHull; - convexHull.reserve(hull.outer().size()); - - for (auto it = hull.outer().begin(); it != hull.outer().end(); ++it) - { - convexHull.push_back({bg::get<0>(*it), - bg::get<1>(*it)}); - } - - return convexHull; -} - -struct BoundingBox2D { - double xmin; - double xmax; - double ymin; - double ymax; -}; - -/** - * Quickly check whether this hyperrectangle excludes grid code zero - * in any individual module. - */ -bool tryProveGridCodeZeroImpossible_1D( - const vector>>& domainToPlaneByModule, - const vector>& latticeBasisByModule, - const vector>& inverseLatticeBasisByModule, - size_t numDims, - const double x0[], - const double dims[], - double r, - double rSquared) -{ - const double point1 = x0[0]; - const double point2 = x0[0] + dims[0]; - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const pair p1 = transformND(domainToPlaneByModule[iModule], - &point1); - const pair p2 = transformND(domainToPlaneByModule[iModule], - &point2); - - // Figure out which lattice points we need to check. - const double xmin = std::min(p1.first, p2.first); - const double xmax = std::max(p1.first, p2.first); - const double ymin = std::min(p1.second, p2.second); - const double ymax = std::max(p1.second, p2.second); - LatticePointEnumerator latticePoints(latticeBasisByModule[iModule], - inverseLatticeBasisByModule[iModule], - xmin, xmax, ymin, ymax, r, rSquared); - - pair latticePoint; - bool foundLatticeCollision = false; - while (!foundLatticeCollision && latticePoints.getNext(&latticePoint)) - { - foundLatticeCollision = - lineSegmentIntersectsCircle(p1, p2, latticePoint, rSquared); - } - - if (!foundLatticeCollision) - { - // This module never gets near grid code zero for the provided range of - // locations. So this range can't possibly contain grid code zero. - return true; - } - } - - return false; -} - -const double shadowWidthThreshold = 0.5; - -BoundingBox2D computeBoundingBox(const vector>& shadow) -{ - BoundingBox2D boundingBox = { - std::numeric_limits::max(), - std::numeric_limits::lowest(), - std::numeric_limits::max(), - std::numeric_limits::lowest(), - }; - - for (const pair& p : shadow) - { - boundingBox.xmin = std::min(boundingBox.xmin, p.first); - boundingBox.xmax = std::max(boundingBox.xmax, p.first); - boundingBox.ymin = std::min(boundingBox.ymin, p.second); - boundingBox.ymax = std::max(boundingBox.ymax, p.second); - } - - return boundingBox; -} - -vector computeShadowLines(const vector>& points) -{ - vector lines; - lines.reserve(points.size()); - - for (size_t iEnd = 0; iEnd < points.size(); ++iEnd) - { - const pair& p1 = (iEnd == 0) - ? points[points.size() - 1] - : points[iEnd - 1]; - const pair& p2 = points[iEnd]; - - pair unitVector = {p2.first - p1.first, - p2.second - p1.second}; - const double lineLength = sqrt(pow(unitVector.first, 2) + - pow(unitVector.second, 2)); - unitVector.first /= lineLength; - unitVector.second /= lineLength; - - lines.push_back({unitVector, lineLength}); - } - - return lines; -} - -PolygonData computePolygonData(const vector>& shadow) -{ - // Compute the polygon's area times two. - double acc = 0; - for (size_t iStart = 0; iStart < shadow.size() - 1; ++iStart) - { - const pair& p1 = shadow[iStart]; - const pair& p2 = shadow[iStart + 1]; - - acc += p1.first * p2.second; - acc -= p2.first * p1.second; - } - - const bool isValidPolygon = (acc != 0); - - if (isValidPolygon) - { - // Make a copy that doesn't have the duplicate at the end. - vector> points(shadow.begin(), - shadow.end() - 1); - - // Compute centroid. - double xsum = 0; - double ysum = 0; - for (const pair& p : points) - { - xsum += p.first; - ysum += p.second; - } - const pair centroid = {xsum / points.size(), - ysum / points.size()}; - - // Compute thetas. - vector thetas; - thetas.reserve(points.size()); - for (const pair& p : points) - { - thetas.push_back(getThetaIndex(p.first - centroid.first, - p.second - centroid.second)); - } - - // Sort by theta. - { - vector indices(points.size()); - std::iota(indices.begin(), indices.end(), 0); - std::sort(indices.begin(), indices.end(), - [&](size_t a, size_t b) { - return thetas[a] < thetas[b]; - }); - vector> pointsSorted; - pointsSorted.reserve(points.size()); - vector thetasSorted; - thetasSorted.reserve(points.size()); - for (size_t idx : indices) - { - pointsSorted.push_back(points[idx]); - thetasSorted.push_back(thetas[idx]); - } - points = std::move(pointsSorted); - thetas = std::move(thetasSorted); - } - - // Compute half planes. - vector halfPlanes; - halfPlanes.reserve(points.size()); - for (size_t i = 0; i < points.size(); ++i) - { - const pair& p1 = (i == 0) - ? points[points.size() - 1] - : points[i - 1]; - const pair& p2 = points[i]; - const pair normalVector = {p2.second - p1.second, - -(p2.first - p1.first)}; - const double top = (normalVector.first*p1.first + - normalVector.second*p1.second); - halfPlanes.push_back({normalVector, top}); - } - - return { - isValidPolygon, - centroid, - points, - thetas, - halfPlanes, - computeShadowLines(points), - }; - } - else - { - bool useX = false; - for (size_t iStart = 0; iStart < shadow.size() - 1; ++iStart) - { - if (shadow[iStart].first != shadow[iStart+1].first) - { - useX = true; - break; - } - } - - vector> points; - points.reserve(2); - if (useX) - { - auto compare = - [](const pair& a, const pair& b) - { - return a.first < b.first; - }; - - points.push_back(*std::min_element(shadow.begin(), shadow.end(), compare)); - points.push_back(*std::max_element(shadow.begin(), shadow.end(), compare)); - } - else - { - auto compare = - [](const pair& a, const pair& b) - { - return a.second < b.second; - }; - - points.push_back(*std::min_element(shadow.begin(), shadow.end(), compare)); - points.push_back(*std::max_element(shadow.begin(), shadow.end(), compare)); - } - - return { - isValidPolygon, - {INFINITY,INFINITY}, - points, - {}, - {}, - computeShadowLines(points), - }; - } -} - -LatticeBox computeLatticeBox( - const BoundingBox2D& boundingBox, - const SquareMatrix2D& inverseLatticeBasis, - double r) -{ - const double paddedLeft = boundingBox.xmin - r; - const double paddedRight = boundingBox.xmax + r; - const double paddedBottom = boundingBox.ymin - r; - const double paddedTop = boundingBox.ymax + r; - - double xmin = std::numeric_limits::max(); - double xmax = std::numeric_limits::lowest(); - - pair ij; - ij = transform2D(inverseLatticeBasis, {paddedLeft, paddedBottom}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedRight, paddedBottom}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedLeft, paddedTop}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - ij = transform2D(inverseLatticeBasis, {paddedRight, paddedTop}); - xmin = std::min(xmin, ij.first); - xmax = std::max(xmax, ij.first); - - return {xmin, xmax, - transform2D(inverseLatticeBasis, {(paddedRight + paddedLeft) / 2, - (paddedTop + paddedBottom) / 2})}; -} - -/** - * Quickly check whether this hyperrectangle excludes grid code zero - * in any individual module. - */ -bool tryProveGridCodeZeroImpossible( - const vector>>& domainToPlaneByModule, - const vector>& latticeBasisByModule, - const vector>& inverseLatticeBasisByModule, - size_t numDims, - const double x0[], - const double dims[], - double r, - double rSquared, - double vertexBuffer[], - vector>& cachedShadows, - vector>& cachedShadowBoundingBoxes, - vector>& cachedLatticeBoxes, - size_t frameNumber) -{ - if (numDims == 1) - { - return tryProveGridCodeZeroImpossible_1D( - domainToPlaneByModule, latticeBasisByModule, inverseLatticeBasisByModule, - numDims, x0, dims, r, rSquared); - } - - NTA_ASSERT(frameNumber <= cachedShadowBoundingBoxes.size()); - - if (frameNumber == cachedShadowBoundingBoxes.size()) - { - vector shadowByModule; - shadowByModule.reserve(domainToPlaneByModule.size()); - - vector boundingBoxByModule; - boundingBoxByModule.reserve(domainToPlaneByModule.size()); - - vector latticeBoxByModule; - latticeBoxByModule.reserve(domainToPlaneByModule.size()); - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const vector> shadow = - getShadowConvexHull(domainToPlaneByModule[iModule], numDims, dims, - vertexBuffer); - - const BoundingBox2D boundingBox = computeBoundingBox(shadow);; - boundingBoxByModule.push_back(boundingBox); - - latticeBoxByModule.push_back( - computeLatticeBox(boundingBox, inverseLatticeBasisByModule[iModule], - r)); - - if (boundingBox.xmax - boundingBox.xmin > shadowWidthThreshold || - boundingBox.ymax - boundingBox.ymin > shadowWidthThreshold) - { - shadowByModule.push_back({}); - } - else - { - shadowByModule.push_back(computePolygonData(shadow)); - } - } - - cachedShadows.push_back(shadowByModule); - cachedShadowBoundingBoxes.push_back(boundingBoxByModule); - cachedLatticeBoxes.push_back(latticeBoxByModule); - } - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - // Figure out which lattice points we need to check. - const pair shift = - transformND(domainToPlaneByModule[iModule], x0); - const BoundingBox2D& boundingBox = - cachedShadowBoundingBoxes[frameNumber][iModule]; - const double xmin = boundingBox.xmin + shift.first; - const double xmax = boundingBox.xmax + shift.first; - const double ymin = boundingBox.ymin + shift.second; - const double ymax = boundingBox.ymax + shift.second; - - LatticePointEnumerator latticePoints( - latticeBasisByModule[iModule], inverseLatticeBasisByModule[iModule], - cachedLatticeBoxes[frameNumber][iModule], shift, xmin, xmax, ymin, ymax, - rSquared); - - pair latticePoint; - bool foundLatticeCollision = false; - while (!foundLatticeCollision && latticePoints.getNext(&latticePoint)) - { - // At this point, the bounding box collides with a lattice point, but it - // might not collide with the actual polygon formed by the 2D shadow of - // the box. We don't actually need to check whether it collides with the - // polygon; we can just say that it might. If it doesn't, this function - // will get called again with a smaller box, and eventually the space will - // be broken into sufficiently small boxes that each have no overlap with - // the lattice of at least one module. - // - // With high dimensional boxes, this approach can be slow. It checks a - // large number of tiny polygons that are just outside the range of a - // lattice point. It has to divide each of these tiny polygons until the - // bounding box doesn't touch the lattice point. With a thorough approach, - // it has to perform much fewer checks, but each check is slower. - // - // To get the best of both worlds, we do non-thorough checks when the - // shadow is large, and begin doing thorough checks when the shadow is - // small. - if (xmax - xmin > shadowWidthThreshold || - ymax - ymin > shadowWidthThreshold) - { - // Rely on the bounding box check. - foundLatticeCollision = true; - } - else - { - latticePoint.first -= shift.first; - latticePoint.second -= shift.second; - foundLatticeCollision = - pointDistanceFromShadowSquared(latticePoint, - cachedShadows[frameNumber][iModule], - numDims) <= rSquared; - } - } - - if (!foundLatticeCollision) - { - // This module never gets near grid code zero for the provided range of - // locations. So this range can't possibly contain grid code zero. - return true; - } - } - - return false; -} - -/** - * Temporarily set a value. - */ -class SwapValueRAII -{ -public: - SwapValueRAII(double* value, double newValue) - :value_(value), initial_(*value) - { - *value_ = newValue; - } - - ~SwapValueRAII() - { - *value_ = initial_; - } - -private: - double *value_; - double initial_; -}; - -/** - * Helper function that doesn't allocate any memory, so it's much better for - * recursion. - */ -bool findGridCodeZeroHelper( - const vector>>& domainToPlaneByModule, - const vector>& latticeBasisByModule, - const vector>& inverseLatticeBasisByModule, - size_t numDims, - double x0[], - double dims[], - double r, - double rSquaredPositive, - double rSquaredNegative, - double vertexBuffer[], - vector>& cachedShadows, - vector>& cachedShadowBoundingBoxes, - vector>& cachedLatticeBoxes, - size_t frameNumber, - std::atomic& shouldContinue) -{ - if (!shouldContinue) - { - return false; - } - - if (tryProveGridCodeZeroImpossible(domainToPlaneByModule, - latticeBasisByModule, - inverseLatticeBasisByModule, numDims, x0, - dims, r, rSquaredNegative, vertexBuffer, - cachedShadows, cachedShadowBoundingBoxes, - cachedLatticeBoxes, frameNumber)) - { - return false; - } - - if (tryFindGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - inverseLatticeBasisByModule, numDims, x0, dims, - rSquaredPositive, vertexBuffer)) - { - return true; - } - - size_t iWidestDim = std::distance(dims, - std::max_element(dims, dims + numDims)); - { - SwapValueRAII swap1(&dims[iWidestDim], dims[iWidestDim] / 2); - if (findGridCodeZeroHelper( - domainToPlaneByModule, latticeBasisByModule, - inverseLatticeBasisByModule, numDims, x0, dims, r, rSquaredPositive, - rSquaredNegative, vertexBuffer, cachedShadows, - cachedShadowBoundingBoxes, cachedLatticeBoxes, frameNumber + 1, - shouldContinue)) - { - return true; - } - - { - SwapValueRAII swap2(&x0[iWidestDim], x0[iWidestDim] + dims[iWidestDim]); - return findGridCodeZeroHelper( - domainToPlaneByModule, latticeBasisByModule, - inverseLatticeBasisByModule, numDims, x0, dims, r, rSquaredPositive, - rSquaredNegative, vertexBuffer, cachedShadows, - cachedShadowBoundingBoxes, cachedLatticeBoxes, frameNumber + 1, - shouldContinue); - } - } -} - -struct GridUniquenessState { - // Constants (thread-safe) - const vector>>& domainToPlaneByModule; - const vector>& latticeBasisByModule; - const vector>& inverseLatticeBasisByModule; - const double readoutResolution; - const double meanScaleEstimate; - const size_t numDims; - - // Task management - double baselineRadius; - double expansionRadiusGoal; - vector expansionProgress; - size_t expandingDim; - bool positiveExpand; - bool continueExpansion; - - // Results - vector pointWithGridCodeZero; - double foundPointBaselineRadius; - - // Thread management - std::mutex& mutex; - std::condition_variable& finishedCondition; - bool finished; - size_t numActiveThreads; - vector threadBaselineRadius; - vector> threadQueryX0; - vector> threadQueryDims; - vector> threadShouldContinue; - std::atomic& quitting; - vector threadRunning; -}; - -template -std::string vecs(const vector& v) -{ - std::ostringstream oss; - oss << "["; - - for (size_t i = 0; i < v.size(); i++) - { - oss << v[i]; - if (i < v.size() - 1) - { - oss << ", "; - } - } - - oss << "]"; - - return oss.str(); -} - -void recordResult(size_t iThread, GridUniquenessState& state, - const vector& pointWithGridCodeZero) -{ - state.continueExpansion = false; - if (state.threadBaselineRadius[iThread] < state.foundPointBaselineRadius) - { - state.foundPointBaselineRadius = state.threadBaselineRadius[iThread]; - state.pointWithGridCodeZero = pointWithGridCodeZero; - - // Notify all others that they should stop unless they're checking a lower - // base width. - for (size_t iOtherThread = 0; - iOtherThread < state.threadBaselineRadius.size(); - iOtherThread++) - { - if (iOtherThread != iThread && - state.threadShouldContinue[iOtherThread] && - (state.threadBaselineRadius[iOtherThread] >= - state.foundPointBaselineRadius)) - { - state.threadShouldContinue[iOtherThread] = false; - } - } - } -} - -void claimNextTask(size_t iThread, GridUniquenessState& state) -{ - state.threadBaselineRadius[iThread] = state.baselineRadius; - - vector& x0 = state.threadQueryX0[iThread]; - vector& dims = state.threadQueryDims[iThread]; - - // Determine all but the final dimension. - for (size_t iDim = 0; iDim < state.numDims - 1; iDim++) - { - dims[iDim] = 2*state.expansionProgress[iDim]; - x0[iDim] = -state.expansionProgress[iDim]; - } - - // 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. - dims[state.numDims - 1] = state.expansionProgress[state.numDims - 1]; - x0[state.numDims - 1] = 0; - - // Make the changes specific to this query. - dims[state.expandingDim] = state.expansionRadiusGoal - state.baselineRadius; - x0[state.expandingDim] = state.positiveExpand - ? state.baselineRadius - : -state.expansionRadiusGoal; - - // Update. - if (state.positiveExpand && - // Optimization: Don't check the negative for the final - // dimension. (described above) - state.expandingDim < state.numDims - 1) - { - state.positiveExpand = false; - } - else - { - state.positiveExpand = true; - state.expansionProgress[state.expandingDim] = state.expansionRadiusGoal; - if (++state.expandingDim >= state.numDims) - { - state.baselineRadius = state.expansionRadiusGoal; - state.expansionRadiusGoal *= 1.01; - state.expandingDim = 0; - } - } -} - -void findGridCodeZeroThread(size_t iThread, GridUniquenessState& state) -{ - bool foundGridCodeZero = false; - vector x0(state.numDims); - vector dims(state.numDims); - vector pointWithGridCodeZero(state.numDims); - - vector numBinsByDim(state.numDims); - - // Add a small epsilon to handle situations where floating point math causes - // a vertex to be non-zero-overlapping here and zero-overlapping in - // tryProveGridCodeZeroImpossible. With this addition, anything - // zero-overlapping in tryProveGridCodeZeroImpossible is guaranteed to be - // zero-overlapping here, so the program won't get caught in infinite - // recursion. - const double rSquaredPositive = pow(state.readoutResolution/2 + 0.000000001, 2); - const double rSquaredNegative = pow(state.readoutResolution/2, 2); - - while (!state.quitting) - { - // Modify the shared state. Record the results, decide the next task, - // volunteer to do it. - { - std::lock_guard lock(state.mutex); - - if (foundGridCodeZero) - { - recordResult(iThread, state, pointWithGridCodeZero); - } - - if (!state.continueExpansion) - { - break; - } - - // Select task params. - claimNextTask(iThread, state); - - // Make an unshared copy that findGridCodeZeroHelper can modify. - x0 = state.threadQueryX0[iThread]; - dims = state.threadQueryDims[iThread]; - } - - // Perform the task. - - vector> cachedShadows; - vector> cachedShadowBoundingBoxes; - vector> cachedLatticeBoxes; - - // Optimization: if the box is large, break it into small chunks rather than - // relying completely on the divide-and-conquer to break into - // reasonable-sized chunks. - - // Use a longer bin size for 1D. A 1D slice of a 2D plane can be relatively - // long before it has high probability of colliding with a lattice point in - // every module. - const double scalesPerBin = (state.numDims == 1) - ? 2.5 - : 0.55; - - for (size_t iDim = 0; iDim < state.numDims; iDim++) - { - numBinsByDim[iDim] = ceil(dims[iDim] / (scalesPerBin * - state.meanScaleEstimate)); - dims[iDim] /= numBinsByDim[iDim]; - } - - const vector& x0_orig = state.threadQueryX0[iThread]; - vector currentBinByDim(state.numDims, 0); - while (state.threadShouldContinue[iThread]) - { - for (size_t iDim = 0; iDim < state.numDims; iDim++) - { - x0[iDim] = x0_orig[iDim] + currentBinByDim[iDim]*dims[iDim]; - } - - foundGridCodeZero = findGridCodeZeroHelper( - state.domainToPlaneByModule, state.latticeBasisByModule, - state.inverseLatticeBasisByModule, state.numDims, x0.data(), - dims.data(), state.readoutResolution/2, rSquaredPositive, - rSquaredNegative, pointWithGridCodeZero.data(), cachedShadows, - cachedShadowBoundingBoxes, cachedLatticeBoxes, 0, - state.threadShouldContinue[iThread]); - - if (foundGridCodeZero) break; - - // Increment as little endian arithmetic with a varying base. - bool overflow = true; - for (size_t iDigit = 0; iDigit < state.numDims; iDigit++) - { - overflow = ++currentBinByDim[iDigit] == numBinsByDim[iDigit]; - if (!overflow) break; - currentBinByDim[iDigit] = 0; - } - - if (overflow) break; - } - } - - // This thread is exiting. - { - std::lock_guard lock(state.mutex); - if (--state.numActiveThreads == 0) - { - state.finished = true; - state.finishedCondition.notify_all(); - } - state.threadRunning[iThread] = false; - } -} - -pair rotateClockwise(double theta, double x, double y) -{ - return {cos(theta)*x + sin(theta)*y, - -sin(theta)*x + cos(theta)*y}; -} - -/** - * Change the matrices so that movement tends to be along the x axis on the - * plane. The end result is the same, but this improves performance because we - * draw bounding boxes around the projected hyperrectangles. This is especially - * beneficial in 1D, totally eliminating diagonal motion so that the bounding - * box perfectly encloses the projected line. - */ -void optimizeMatrices(vector>> *domainToPlaneByModule, - vector>> *latticeBasisByModule) -{ - for (size_t iModule = 0; iModule < domainToPlaneByModule->size(); iModule++) - { - vector> &domainToPlane = (*domainToPlaneByModule)[iModule]; - vector> &latticeBasis = (*latticeBasisByModule)[iModule]; - - size_t iLongest = (size_t) -1; - double dLongest = std::numeric_limits::lowest(); - for (size_t iColumn = 0; iColumn < domainToPlane[0].size(); iColumn++) - { - double length = sqrt(pow(domainToPlane[0][iColumn], 2) + - pow(domainToPlane[1][iColumn], 2)); - if (length > dLongest) - { - dLongest = length; - iLongest = iColumn; - } - } - - const double theta = atan2(domainToPlane[1][iLongest], - domainToPlane[0][iLongest]); - for (size_t iColumn = 0; iColumn < domainToPlane[0].size(); iColumn++) - { - const pair newColumn = rotateClockwise(theta, - domainToPlane[0][iColumn], - domainToPlane[1][iColumn]); - domainToPlane[0][iColumn] = newColumn.first; - domainToPlane[1][iColumn] = newColumn.second; - } - for (size_t iColumn = 0; iColumn < latticeBasis[0].size(); iColumn++) - { - const pair newColumn = rotateClockwise(theta, - latticeBasis[0][iColumn], - latticeBasis[1][iColumn]); - latticeBasis[0][iColumn] = newColumn.first; - latticeBasis[1][iColumn] = newColumn.second; - } - } -} - -bool nupic::experimental::grid_uniqueness::findGridCodeZero( - const vector>>& domainToPlaneByModule, - const vector>>& latticeBasisByModule, - const vector& x0, - const vector& dims, - double readoutResolution, - vector* pointWithGridCodeZero) -{ - // Avoid doing any allocations in each recursion. - vector x0Copy(x0); - vector dimsCopy(dims); - std::atomic shouldContinue(true); - - vector defaultPointBuffer; - - if (pointWithGridCodeZero != nullptr) - { - NTA_ASSERT(pointWithGridCodeZero->size() == dims.size()); - } - else - { - defaultPointBuffer.resize(dims.size()); - pointWithGridCodeZero = &defaultPointBuffer; - } - - NTA_ASSERT(domainToPlaneByModule[0].size() == 2); - - vector>> domainToPlaneByModule2(domainToPlaneByModule); - vector>> latticeBasisByModule2(latticeBasisByModule); - optimizeMatrices(&domainToPlaneByModule2, &latticeBasisByModule2); - - vector> latticeBasisByModule3; - vector> inverseLatticeBasisByModule; - for (const vector>& latticeBasis : latticeBasisByModule2) - { - latticeBasisByModule3.push_back({ - latticeBasis[0][0], latticeBasis[0][1], - latticeBasis[1][0], latticeBasis[1][1]}); - inverseLatticeBasisByModule.push_back(invert2DMatrix(latticeBasis)); - } - - vector> cachedShadows; - vector> cachedShadowBoundingBoxes; - vector> cachedLatticeBoxes; - - // Add a small epsilon to handle situations where floating point math causes a - // vertex to be non-zero-overlapping here and zero-overlapping in - // tryProveGridCodeZeroImpossible. With this addition, anything - // zero-overlapping in tryProveGridCodeZeroImpossible is guaranteed to be - // zero-overlapping here, so the program won't get caught in infinite - // recursion. - const double rSquaredPositive = pow(readoutResolution/2 + 0.000000001, 2); - const double rSquaredNegative = pow(readoutResolution/2, 2); - - return findGridCodeZeroHelper( - domainToPlaneByModule2, latticeBasisByModule3, inverseLatticeBasisByModule, - dimsCopy.size(), x0Copy.data(), dimsCopy.data(), readoutResolution/2, - rSquaredPositive, rSquaredNegative, pointWithGridCodeZero->data(), - cachedShadows, cachedShadowBoundingBoxes, cachedLatticeBoxes, 0, - shouldContinue); -} - -pair> -nupic::experimental::grid_uniqueness::computeGridUniquenessHypercube( - const vector>>& domainToPlaneByModule, - const vector>>& latticeBasisByModule, - double readoutResolution, - double ignoredCenterDiameter, - double pingInterval) -{ - typedef std::chrono::steady_clock Clock; - - enum ExitReason { - Timeout, - Interrupt, - Completed - }; - - std::atomic exitReason(ExitReason::Completed); - - ThreadSafeQueue messages; - CaptureInterruptsRAII captureInterrupts(&messages); - - std::atomic quitting(false); - std::thread messageThread( - [&]() { - while (true) - { - switch (messages.take()) - { - case Message::Interrupt: - quitting = true; - exitReason = ExitReason::Interrupt; - break; - case Message::Timeout: - quitting = true; - exitReason = ExitReason::Timeout; - break; - case Message::Exiting: - return; - } - } - }); - - NTA_CHECK(domainToPlaneByModule.size() == latticeBasisByModule.size()) - << "The two arrays of matrices must be the same length (one per module) " - << "Actual: " << domainToPlaneByModule.size() - << " " << latticeBasisByModule.size(); - - NTA_CHECK(domainToPlaneByModule[0].size() == 2) - << "Each matrix should have two rows -- the modules are two-dimensional. " - << "Actual: " << domainToPlaneByModule[0].size(); - - NTA_CHECK(latticeBasisByModule[0][0].size() == 2) - << "There should be two lattice basis vectors. " - << "Actual: " << latticeBasisByModule[0][0].size(); - - const size_t numDims = domainToPlaneByModule[0][0].size(); - NTA_CHECK(numDims < sizeof(int)*8) - << "Unsupported number of dimensions: " << numDims; - - vector>> domainToPlaneByModule2(domainToPlaneByModule); - vector>> latticeBasisByModule2(latticeBasisByModule); - optimizeMatrices(&domainToPlaneByModule2, &latticeBasisByModule2); - - vector> latticeBasisByModule3; - vector> inverseLatticeBasisByModule; - for (const vector>& latticeBasis : latticeBasisByModule2) - { - latticeBasisByModule3.push_back({ - latticeBasis[0][0], latticeBasis[0][1], - latticeBasis[1][0], latticeBasis[1][1]}); - inverseLatticeBasisByModule.push_back(invert2DMatrix(latticeBasis)); - } - - double meanScaleEstimate = 0.0; - - for (const vector>& domainToPlane : domainToPlaneByModule2) - { - double longestDisplacementSquared = std::numeric_limits::min(); - - for (size_t iDim = 0; iDim < numDims; iDim++) - { - longestDisplacementSquared = std::max(longestDisplacementSquared, - pow(domainToPlane[0][iDim], 2) + - pow(domainToPlane[1][iDim], 2)); - } - - const double scaleEstimate = 1 / sqrt(longestDisplacementSquared); - meanScaleEstimate += scaleEstimate; - } - meanScaleEstimate /= domainToPlaneByModule2.size(); - - // Use condition_variables to enable periodic logging while waiting for the - // threads to finish. - std::mutex stateMutex; - std::condition_variable finishedCondition; - - size_t numThreads = std::thread::hardware_concurrency(); - - GridUniquenessState state = { - domainToPlaneByModule2, - latticeBasisByModule3, - inverseLatticeBasisByModule, - readoutResolution, - - meanScaleEstimate, - numDims, - - ignoredCenterDiameter, - ignoredCenterDiameter * 1.01, - vector(numDims, ignoredCenterDiameter), - 0, - true, - true, - - vector(numDims), - std::numeric_limits::max(), - - stateMutex, - finishedCondition, - false, - 0, - vector(numThreads, std::numeric_limits::max()), - vector>(numThreads, vector(numDims)), - vector>(numThreads, vector(numDims)), - vector>(numThreads), - quitting, - vector(numThreads, true), - }; - - for (size_t i = 0; i < numThreads; i++) - { - state.threadShouldContinue[i] = true; - } - - { - std::unique_lock lock(stateMutex); - for (size_t i = 0; i < numThreads; i++) - { - std::thread(findGridCodeZeroThread, i, std::ref(state)).detach(); - state.numActiveThreads++; - } - - const auto tStart = Clock::now(); - auto tNextPrint = tStart + std::chrono::duration(pingInterval); - - bool printedInitialStatement = false; - - while (!state.finished) - { - if (pingInterval <= 0) - { - state.finishedCondition.wait(lock); - } - else - { - if (state.finishedCondition.wait_until( - lock, tNextPrint) == std::cv_status::timeout) - { - if (!printedInitialStatement) - { - { - std::ostringstream oss; - oss << "["; - - for (size_t iModule = 0; - iModule < domainToPlaneByModule.size(); - iModule++) - { - oss << "["; - oss << vecs(domainToPlaneByModule[iModule][0]) << ","; - oss << vecs(domainToPlaneByModule[iModule][1]); - oss << "],"; - } - oss << "]" << std::endl; - NTA_INFO << "domainToPlaneByModule:" << std::endl << oss.str(); - } - - { - std::ostringstream oss; - oss << "["; - for (size_t iModule = 0; - iModule < latticeBasisByModule.size(); - iModule++) - { - oss << "["; - oss << vecs(latticeBasisByModule[iModule][0]) << ","; - oss << vecs(latticeBasisByModule[iModule][1]); - oss << "],"; - } - oss << "]" << std::endl; - - NTA_INFO << "latticeBasisByModule:" << std::endl << oss.str(); - } - - NTA_INFO << "readout resolution: " << readoutResolution; - - printedInitialStatement = true; - } - - NTA_INFO << ""; - NTA_INFO << domainToPlaneByModule.size() << " modules, " << numDims - << " dimensions, " - << std::chrono::duration_cast( - Clock::now() - tStart).count() << " seconds elapsed"; - - if (state.foundPointBaselineRadius < - std::numeric_limits::max()) - { - NTA_INFO << "**Hypercube side length upper bound: " - << state.foundPointBaselineRadius << "**"; - NTA_INFO << "**Grid code zero found at: " - << vecs(state.pointWithGridCodeZero) << "**"; - } - - tNextPrint = (Clock::now() + - std::chrono::duration(pingInterval)); - - for (size_t iThread = 0; iThread < state.threadBaselineRadius.size(); - iThread++) - { - if (state.threadRunning[iThread]) - { - if (state.threadShouldContinue[iThread]) - { - NTA_INFO << " Thread " << iThread - << " assuming hypercube side length lower bound " - << state.threadBaselineRadius[iThread] - << ", querying x0 " - << vecs(state.threadQueryX0[iThread]) << " and dims " - << vecs(state.threadQueryDims[iThread]); - } - else - { - NTA_INFO << " Thread " << iThread - << " has been ordered to stop."; - } - } - else - { - NTA_INFO << " Thread " << iThread << " is finished."; - } - } - } - } - } - } - - 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: - return {state.foundPointBaselineRadius, state.pointWithGridCodeZero}; - } -} - -bool tryFindGridCodeZero_noModulo( - const vector>>& domainToPlaneByModule, - size_t numDims, - const double x0[], - const double dims[], - double rSquared, - double vertexBuffer[]) -{ - for (size_t iDim = 0; iDim < numDims; iDim++) - { - vertexBuffer[iDim] = x0[iDim] + (dims[iDim]/2); - } - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const pair pointOnPlane = - transformND(domainToPlaneByModule[iModule], vertexBuffer); - - if (pow(pointOnPlane.first, 2) + pow(pointOnPlane.second, 2) > rSquared) - { - return false; - } - } - - return true; -} - -bool tryProveGridCodeZeroImpossible_noModulo_1D( - const vector>>& domainToPlaneByModule, - size_t numDims, - const double x0[], - const double dims[], - double rSquared) -{ - const double point1 = x0[0]; - const double point2 = x0[0] + dims[0]; - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const pair p1 = transformND(domainToPlaneByModule[iModule], - &point1); - const pair p2 = transformND(domainToPlaneByModule[iModule], - &point2); - - if (!lineSegmentIntersectsCircle(p1, p2, {0.0, 0.0}, rSquared)) - { - // This module never gets near grid code zero for the provided range of - // locations. So this range can't possibly contain grid code zero. - return true; - } - } - - return false; -} - -bool tryProveGridCodeZeroImpossible_noModulo( - const vector>>& domainToPlaneByModule, - size_t numDims, - const double x0[], - const double dims[], - double r, - double rSquared, - double vertexBuffer[], - vector>& cachedShadows, - size_t frameNumber) -{ - if (numDims == 1) - { - return tryProveGridCodeZeroImpossible_noModulo_1D( - domainToPlaneByModule, numDims, x0, dims, rSquared); - } - - NTA_ASSERT(frameNumber <= cachedShadows.size()); - - if (frameNumber == cachedShadows.size()) - { - vector shadowByModule; - shadowByModule.reserve(domainToPlaneByModule.size()); - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const vector> shadow = getShadowConvexHull( - domainToPlaneByModule[iModule], numDims, dims, vertexBuffer); - shadowByModule.push_back(computePolygonData(shadow)); - } - - cachedShadows.push_back(shadowByModule); - } - - for (size_t iModule = 0; iModule < domainToPlaneByModule.size(); iModule++) - { - const pair shift = - transformND(domainToPlaneByModule[iModule], x0); - - if (!(pointDistanceFromShadowSquared({-shift.first, -shift.second}, - cachedShadows[frameNumber][iModule], - numDims) <= rSquared)) - { - // This module never gets near grid code zero for the provided range of - // locations. So this range can't possibly contain grid code zero. - return true; - } - } - - return false; -} - -bool findGridCodeZeroHelper_noModulo( - const vector>>& domainToPlaneByModule, - size_t numDims, - double x0[], - double dims[], - double r, - double rSquaredPositive, - double rSquaredNegative, - double vertexBuffer[], - vector>& cachedShadows, - size_t frameNumber, - std::atomic& shouldContinue) -{ - if (!shouldContinue) - { - return false; - } - - if (tryProveGridCodeZeroImpossible_noModulo( - domainToPlaneByModule, numDims, x0, dims, r, rSquaredNegative, - vertexBuffer, cachedShadows, frameNumber)) - { - return false; - } - - if (tryFindGridCodeZero_noModulo(domainToPlaneByModule, numDims, x0, dims, - rSquaredPositive, vertexBuffer)) - { - return true; - } - - size_t iWidestDim = std::distance(dims, - std::max_element(dims, dims + numDims)); - { - SwapValueRAII swap1(&dims[iWidestDim], dims[iWidestDim] / 2); - if (findGridCodeZeroHelper_noModulo( - domainToPlaneByModule, numDims, x0, dims, r, rSquaredPositive, - rSquaredNegative, vertexBuffer, cachedShadows, frameNumber + 1, - shouldContinue)) - { - return true; - } - - { - SwapValueRAII swap2(&x0[iWidestDim], x0[iWidestDim] + dims[iWidestDim]); - return findGridCodeZeroHelper_noModulo( - domainToPlaneByModule, numDims, x0, dims, r, rSquaredPositive, - rSquaredNegative, vertexBuffer, cachedShadows, frameNumber + 1, - shouldContinue); - } - } -} - -bool findGridCodeZero_noModulo( - const vector>>& domainToPlaneByModule, - const vector& x0, - const vector& dims, - double readoutResolution, - std::atomic& shouldContinue, - vector* pointWithGridCodeZero = nullptr) -{ - // Avoid doing any allocations in each recursion. - vector x0Copy(x0); - vector dimsCopy(dims); - - vector defaultPointBuffer; - - if (pointWithGridCodeZero != nullptr) - { - NTA_ASSERT(pointWithGridCodeZero->size() == dims.size()); - } - else - { - defaultPointBuffer.resize(dims.size()); - pointWithGridCodeZero = &defaultPointBuffer; - } - - NTA_ASSERT(domainToPlaneByModule[0].size() == 2); - - vector> cachedShadows; - - // Add a small epsilon to handle situations where floating point math causes a - // vertex to be non-zero-overlapping here and zero-overlapping in - // tryProveGridCodeZeroImpossible. With this addition, anything - // zero-overlapping in tryProveGridCodeZeroImpossible is guaranteed to be - // zero-overlapping here, so the program won't get caught in infinite - // recursion. - const double rSquaredPositive = pow(readoutResolution/2 + 0.000000001, 2); - const double rSquaredNegative = pow(readoutResolution/2, 2); - - return findGridCodeZeroHelper_noModulo( - domainToPlaneByModule, dimsCopy.size(), x0Copy.data(), dimsCopy.data(), - readoutResolution/2, rSquaredPositive, rSquaredNegative, - pointWithGridCodeZero->data(), cachedShadows, 0, shouldContinue); -} - -bool findGridCodeZeroAtRadius( - double radius, - const vector>>& domainToPlaneByModule, - double readoutResolution, - std::atomic& 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 x0(numDims, -radius); - vector 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; - - // Test -r - if (iDim != numDims - 1) - { - if (findGridCodeZero_noModulo(domainToPlaneByModule, - x0, dims, readoutResolution, - shouldContinue)) - { - return true; - } - } - - // Test +r - x0[iDim] = radius; - if (findGridCodeZero_noModulo(domainToPlaneByModule, - x0, dims, readoutResolution, - shouldContinue)) - { - return true; - } - } - - return false; -} - -double -nupic::experimental::grid_uniqueness::computeBinSidelength( - const vector>>& domainToPlaneByModule, - double readoutResolution, - double resultPrecision, - double upperBound, - double timeout) -{ - // - // Initialization - // - enum ExitReason { - Timeout, - Interrupt, - Completed - }; - - std::atomic exitReason(ExitReason::Completed); - - ThreadSafeQueue messages; - CaptureInterruptsRAII captureInterrupts(&messages); - - std::atomic shouldContinue(true); - std::thread messageThread( - [&]() { - while (true) - { - switch (messages.take()) - { - case Message::Interrupt: - shouldContinue = false; - exitReason = ExitReason::Interrupt; - break; - case Message::Timeout: - shouldContinue = false; - exitReason = ExitReason::Timeout; - break; - case Message::Exiting: - return; - } - } - }); - - ScheduledTask* scheduledTask = nullptr; - if (timeout > 0) - { - scheduledTask = new ScheduledTask( - std::chrono::steady_clock::now() + std::chrono::duration(timeout), - [&messages](){ - messages.put(Message::Timeout); - }); - } - - // - // Computation - // - double tested = 0; - double radius = 0.5; - - while (radius <= upperBound && - findGridCodeZeroAtRadius(radius, - domainToPlaneByModule, - readoutResolution, - shouldContinue)) - { - tested = radius; - radius *= 2; - } - - double result; - if (radius > upperBound) - { - result = -1.0; - } - else - { - // The radius needs to be twice as precise to get the sidelength - // sufficiently precise. - const double resultPrecision2 = resultPrecision / 2; - - double dec = (radius - tested) / 2; - - // The possible error is equal to dec*2. - while (shouldContinue && dec*2 > resultPrecision2) - { - const double testRadius = radius - dec; - - if (!findGridCodeZeroAtRadius(testRadius, - domainToPlaneByModule, - readoutResolution, - shouldContinue)) - { - radius = testRadius; - } - - dec /= 2; - } - - result = 2*radius; - } - - // - // Teardown - // - 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: - return result; - } -} - -vector squeezeRectangleToBin( - const vector>>& domainToPlaneByModule, - double readoutResolution, - double resultPrecision, - double startingRadius, - std::atomic& shouldContinue) -{ - const size_t numDims = domainToPlaneByModule[0][0].size(); - - // The radius needs to be twice as precise to get the sidelength sufficiently - // precise. - const double resultPrecision2 = resultPrecision / 2; - - vector radii(numDims, startingRadius); - - for (size_t iDim = 0; iDim < numDims; ++iDim) - { - double dec = startingRadius / 2; - - // The possible error is equal to dec*2. - while (shouldContinue && dec*2 > resultPrecision2) - { - const double testRadius = radii[iDim] - dec; - - vector x0(numDims); - vector dims(numDims); - - std::transform(radii.begin(), radii.end() - 1, x0.begin(), - [](double r) { return -r; }); - std::transform(radii.begin(), radii.end() - 1, dims.begin(), - [](double r) { return r*2; }); - - // 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] = radii[numDims - 1]; - - // Test two faces of the n-dimensional box by setting the ith dimension - // to -r and +r with width 0. - bool foundZero = false; - dims[iDim] = 0; - 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; - } - } - - return radii; -} - -vector -nupic::experimental::grid_uniqueness::computeBinRectangle( - const vector>>& domainToPlaneByModule, - double readoutResolution, - double resultPrecision, - double upperBound, - double timeout) -{ - // - // Initialization - // - enum ExitReason { - Timeout, - Interrupt, - Completed - }; - - std::atomic exitReason(ExitReason::Completed); - - ThreadSafeQueue messages; - CaptureInterruptsRAII captureInterrupts(&messages); - - std::atomic shouldContinue(true); - std::thread messageThread( - [&]() { - while (true) - { - switch (messages.take()) - { - case Message::Interrupt: - shouldContinue = false; - exitReason = ExitReason::Interrupt; - break; - case Message::Timeout: - shouldContinue = false; - exitReason = ExitReason::Timeout; - break; - case Message::Exiting: - return; - } - } - }); - - ScheduledTask* scheduledTask = nullptr; - if (timeout > 0) - { - scheduledTask = new ScheduledTask( - std::chrono::steady_clock::now() + std::chrono::duration(timeout), - [&messages](){ - messages.put(Message::Timeout); - }); - } - - // - // Computation - // - double radius = 0.5; - - while (radius <= upperBound && - findGridCodeZeroAtRadius(radius, - domainToPlaneByModule, - readoutResolution, - shouldContinue)) - { - radius *= 2; - } - - vector result; - if (radius > upperBound) - { - // Give up. - } - else - { - const vector radii = squeezeRectangleToBin( - domainToPlaneByModule, readoutResolution, resultPrecision, - radius, shouldContinue); - - result.resize(radii.size()); - std::transform(radii.begin(), radii.end(), result.begin(), - [](double r) { return 2*r; }); - } - - // - // Teardown - // - 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: - return result; - } -} diff --git a/src/nupic/experimental/GridUniqueness.hpp b/src/nupic/experimental/GridUniqueness.hpp deleted file mode 100644 index d4b4e519..00000000 --- a/src/nupic/experimental/GridUniqueness.hpp +++ /dev/null @@ -1,227 +0,0 @@ -/* --------------------------------------------------------------------- - * Numenta Platform for Intelligent Computing (NuPIC) - * Copyright (C) 2018, Numenta, Inc. Unless you have an agreement - * with Numenta, Inc., for a separate license for this software code, the - * following terms and conditions apply: - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Affero Public License version 3 as - * published by the Free Software Foundation. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU Affero Public License for more details. - * - * You should have received a copy of the GNU Affero Public License - * along with this program. If not, see http://www.gnu.org/licenses. - * - * http://numenta.org/licenses/ - * ---------------------------------------------------------------------- - */ - -/** @file - * Functions that analyze the uniqueness of locations in grid cell spaces - */ - -#ifndef NTA_GRID_UNIQUENESS_HPP -#define NTA_GRID_UNIQUENESS_HPP - -#include -#include -#include - -namespace nupic { - namespace experimental { - namespace grid_uniqueness { - - /** - * Determine whether any points in a k-dimensional rectangle have a grid - * code equal to the grid code at the origin. - * - * If this function says there's a grid code zero at point (42.0, 42.0), - * that implies that every location has the same grid code as the location - * that is 42.0 units up and 42.0 units right. - * - * This function uses a recursive "divide and conquer" algorithm. It first - * quickly samples a few points to see if any of them have grid code zero. - * Then it tries to prove that grid code zero can't possibly happen in - * this hyperrectangle by showing that at least one module never has grid - * code zero for any point in this hyperrectangle. Finally, if neither - * attempt succeeds, it divides the hyperrectangle in half, and then it - * tries again on both halves. - * - * @param x0 - * The lowest corner of the k-dimensional rectangle that will be searched. - * - * @param dims - * The dimensions of the k-dimensional rectangle that will be searched. - * - * @param pointWithGridCodeZero - * Output parameter. A point with grid code zero in this hyperrectangle. - * Only populated if the function returns true. - * - * @return - * true if grid code zero is found, false otherwise. - */ - bool findGridCodeZero( - const std::vector>>& domainToPlaneByModule, - const std::vector>>& latticeBasisByModule, - const std::vector& x0, - const std::vector& dims, - Real64 readoutResolution, - std::vector* pointWithGridCodeZero = nullptr); - - /** - * Given a set of grid cell module parameters, determines the diameter of - * the k-dimensional cube in which every location has a unique grid cell - * representation. - * - * This function iteratively builds out a larger hypercube out of smaller - * hyperrectangles, relying on findGridCodeZero to analyze each - * hyperrectangle. The hypercube expands outward from the origin, forming - * a larger hypercube that contains positive and negative numbers. After - * the nearest grid code zero is found, we divide the hypercube in half to - * get a hypercube in which every location is guaranteed to be unique - * (rather than just being different from the starting location). The - * returned diameter is the diameter of this smaller hypercube. - * - * Each grid cell module is specified by a pair of matrices. The first one - * projects a k-dimensional point to a 2D plane, and the second matrix - * specifies the basis of the grid lattice on this plane, specifying which - * points on the plane have equivalent grid codes. The distance between - * two grid codes is equal to the shortest distance between them on the - * plane. The readout resolution parameter is measured in units of this - * distance. - * - * There's not a strictly "correct" way to configure these matrices and - * the readout resolution, but a typical way is to use unit vectors as - * lattice basis vectors, use a plane projection that normalizes distances - * so that 1 unit is 1 "scale", and then set the readout resolution to - * some number in the range (0, 1) so that it is measured in units of - * "scales". - * - * @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 latticeBasisByModule - * A list of m 2*2 matrices, one per module. This matrix contains the - * basis vectors for a lattice, specifying which points on the plane have - * equivalent location representations in this module. - * - * @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 ignoredCenterDiameter - * The diameter of the hypercube at the center which should be ignored, - * because it contains the *actual* grid code zero. Set this to be - * sufficiently large to get away from this actual zero. - * - * @param pingInterval - * How often, in seconds, the function should print its current status. - * If <= 0, no printing will occur. - * - * @return - * - The diameter of the hypercube that contains no collisions. - * - A point just outside this hypercube that collides with the origin. - */ - std::pair> computeGridUniquenessHypercube( - const std::vector>>& domainToPlaneByModule, - const std::vector>>& latticeBasisByModule, - Real64 readoutResolution, - Real64 ignoredCenterDiameter, - Real64 pingInterval=10.0); - - /** - * Compute the sidelength of the smallest hypercube that encloses the - * intersection of all of the modules' firing fields centered at the - * origin. This sidelength is like a resolution of the code. - * - * @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 hypercube it can place around zero, - * stopping at the point when it's within 'resultPrecision' of the actual - * smallest cube. - * - * @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 sidelength of this hypercube. Returns -1.0 if a surface can't be - * found (i.e. if upperBound is reached.) - */ - Real64 computeBinSidelength( - const std::vector>>& domainToPlaneByModule, - Real64 readoutResolution, - Real64 resultPrecision, - Real64 upperBound = 2048.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 computeBinRectangle( - const std::vector>>& domainToPlaneByModule, - Real64 readoutResolution, - Real64 resultPrecision, - Real64 upperBound = 2048.0, - Real64 timeout = -1.0); - } - } -} - -#endif // NTA_GRID_UNIQUENESS_HPP diff --git a/src/test/unit/experimental/GridUniquenessTest.cpp b/src/test/unit/experimental/GridUniquenessTest.cpp deleted file mode 100644 index 0e7af822..00000000 --- a/src/test/unit/experimental/GridUniquenessTest.cpp +++ /dev/null @@ -1,766 +0,0 @@ -/* --------------------------------------------------------------------- - * Numenta Platform for Intelligent Computing (NuPIC) - * Copyright (C) 2018, Numenta, Inc. Unless you have an agreement - * with Numenta, Inc., for a separate license for this software code, the - * following terms and conditions apply: - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Affero Public License version 3 as - * published by the Free Software Foundation. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU Affero Public License for more details. - * - * You should have received a copy of the GNU Affero Public License - * along with this program. If not, see http://www.gnu.org/licenses. - * - * http://numenta.org/licenses/ - * ---------------------------------------------------------------------- - */ - -/** @file - * Unit tests for Grid Uniqueness - */ - -#include - -#include -#include "gtest/gtest.h" - -#include -#include - -using namespace nupic; -using namespace nupic::experimental::grid_uniqueness; -using std::vector; -using std::pair; - -namespace { - TEST(GridUniquenessTest, ZeroInsideSquare) - { - const vector scales = {2, 3, 6, 7, 21}; // Factors of 42 - vector>> domainToPlaneByModule; - vector>> latticeBasisByModule; - for (double scale : scales) - { - domainToPlaneByModule.push_back({ - {1/scale, 0}, - {0, 1/scale}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - } - - // Zero in center of square. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {41.0, 41.0}, {2.0, 2.0}, 0.01)); - - // Zero at bottom-left of square. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {42.0, 42.0}, {2.0, 2.0}, 0.01)); - - // Zero at top-right of square. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {40.0, 40.0}, {2.0, 2.0}, 0.01)); - - // Zero at top-left of square. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {42.0, 40.0}, {2.0, 2.0}, 0.01)); - - // Zero at bottom-right of square. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {40.0, 42.0}, {2.0, 2.0}, 0.01)); - - // Zero in bottom-left quadrant. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {40.5, 40.5}, {2.0, 2.0}, 0.01)); - } - - - TEST(GridUniquenessTest, ZeroOutsideSquare) - { - const vector scales = {2, 3, 6, 7, 21}; // Factors of 42 - vector>> domainToPlaneByModule; - vector>> latticeBasisByModule; - for (double scale : scales) - { - domainToPlaneByModule.push_back({ - {1/scale, 0}, - {0, 1/scale}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - } - - // Zero out-of-bounds of square (bottom-left). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {41.0, 41.0}, {0.5, 0.5}, 0.1)); - - // Zero out-of-bounds of square (top-right). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {42.5, 42.5}, {0.5, 0.5}, 0.1)); - - // Zero out-of-bounds of square (bottom-right). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {42.5, 41.0}, {0.5, 0.5}, 0.1)); - - // Zero out-of-bounds of square (top-left). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {41.0, 42.5}, {0.5, 0.5}, 0.1)); - } - - /** - * Test rectangles just outside and just inside the radius of location zero. - * Specifically focus on the area that would be covered by a square around - * zero but not by a circle. - */ - TEST(GridUniquenessTest, ZeroHasACircularRadius) - { - const vector>> domainToPlaneByModule = { - {{1, 0}, - {0, 1}} - }; - const vector>> latticeBasisByModule = { - {{1, 0}, - {0, 1}} - }; - - const double phaseResolution = 0.2; - - double d = phaseResolution / (2 * sqrt(2)); - - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {d + phaseResolution/100, d + phaseResolution/100}, - {0.2, 0.2}, phaseResolution)); - - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {d - phaseResolution/100, d - phaseResolution/100}, - {0.2, 0.2}, phaseResolution)); - } - - TEST(GridUniquenessTest, 1DSanityCheck) - { - const vector scales = {2, 3, 6, 7, 21}; // Factors of 42 - vector>> domainToPlaneByModule; - vector>> latticeBasisByModule; - for (double scale : scales) - { - // Add two modules at each scale. - domainToPlaneByModule.push_back({ - {1/scale}, - {0}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - domainToPlaneByModule.push_back({ - {0}, - {1/scale}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - } - - // Zero in center of range. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {41.0}, {2.0}, 0.01)); - - // Zero at left of range. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {42.0}, {2.0}, 0.01)); - - // Zero at right of range. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {40.0}, {2.0}, 0.01)); - - // Zero in left half. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {40.5}, {2.0}, 0.01)); - - // Zero out-of-bounds of square (bottom-left). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {41.0}, {0.5}, 0.1)); - - // Zero out-of-bounds of square (top-right). - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {42.5}, {0.5}, 0.1)); - } - - TEST(GridUniquenessTest, NegativeNumbersTest) - { - const vector scales = {2, 3, 6, 7, 21}; // Factors of 42 - vector>> domainToPlaneByModule; - vector>> latticeBasisByModule; - for (double scale : scales) - { - domainToPlaneByModule.push_back({ - {1/scale, 0}, - {0, 1/scale}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - } - - // Zero at a slightly negative number - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {-1.0, 41.0}, {0.999, 2.0}, 0.1)); - - // Zero at a very negative number - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {-43.0, -43.0}, {2.0, 2.0}, 0.01)); - - // Zero out-of-bounds of negative square - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {-1.0, 41.0}, {0.5, 2.0}, 0.1)); - - // Zero out-of-bounds of negative square, both axes, bottom left - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {-43.0, -43.0}, {0.5, 0.5}, 0.1)); - - // Zero out-of-bounds of negative square, both axes, top right - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {-41.5, -41.5}, {0.5, 0.5}, 0.1)); - } - - TEST(GridUniquenessTest, Project3DOnto2D) - { - // Choose scales that will have phase 0.5 at 42. - const vector scales = {4, 12, 28, 84}; - vector>> domainToPlaneByModule; - vector>> latticeBasisByModule; - for (double scale : scales) - { - // Shift the phase by 0.1 for every unit z. - domainToPlaneByModule.push_back({ - {1/scale, 0, 0.1}, - {0, 1/scale, 0.1}, - }); - latticeBasisByModule.push_back({ - {1, 0}, - {0, 1}, - }); - } - - // Verify the z coordinate shifts the phase at 84.0 away from 0 - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {83.5, 83.5, 4.5}, {1.0, 1.0, 1.0}, 0.1)); - - // Verify the z coordinate shifts the phase at 42.0 from 0.5 to 0. - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, {41.5, 41.5, 4.5}, {1.0, 1.0, 1.0}, 0.01)); - } - - /** - * Queries rectangles that are specially selected to cause the - * algorithm to check a 2x2 grid of lattice points. Run 5 different - * tests, one where the rectangle intersects a particular lattice - * point, and one where it doesn't intersect any. - */ - TEST(GridUniquenessTest, ChoosingLatticePoints_OutsideBoundingBox) - { - const vector>> domainToPlaneByModule = { - {{1, 0}, - {0, 1}} - }; - const vector>> latticeBasisByModule = { - {{cos(0), cos(M_PI/3)}, - {sin(0), sin(M_PI/3)}} - }; - - const double readoutResolution = 0.01; - const double padding = 0.02; - - double height = 2*sin(M_PI/3) - padding; - double halfWidth = 0.5 - padding; - double width = 2*halfWidth; - - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-halfWidth, 0}, {width, height}, - readoutResolution)); - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-halfWidth, -height}, {width, height}, - readoutResolution)); - - double halfHeight = sin(M_PI/3) - padding; - height = 2*halfHeight; - width = 1 - padding; - - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-width, -halfHeight}, {width, height}, - readoutResolution)); - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-width, -halfHeight}, {width, height}, - readoutResolution)); - - halfHeight = sin(M_PI/3) - padding; - height = 2*halfHeight; - width = 1 - 2*padding; - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {padding, -halfHeight}, {width, height}, - readoutResolution)); - } - - TEST(GridUniquenessTest, ChoosingLatticePoints_WithinBoundingBoxButNotPolygon) - { - const vector>> latticeBasisByModule = { - {{1, 0}, - {0, 1}} - }; - - const double readoutResolution = 0.01; - const double padding = 0.02; - - for (int iRotation = 0; iRotation < 4; iRotation++) - { - const vector>> domainToPlaneByModule = { - {{cos(iRotation*(M_PI/2) + M_PI/4), cos((iRotation+1)*(M_PI/2) + M_PI/4)}, - {sin(iRotation*(M_PI/2) + M_PI/4), sin((iRotation+1)*(M_PI/2) + M_PI/4)}} - }; - - double width = sqrt(2) - 2*padding; - double height = sqrt(2) - padding; - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-width/2, 0}, {width, height}, - readoutResolution)); - - width = sqrt(2) - 2*padding; - height = sqrt(2) - 2*padding; - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-width/2, padding}, {width, height}, - readoutResolution)); - } - } - - /** - * Lattice points that are inside the bounding box when the bounding - * box is expanded, but where the circles centered at the lattice - * points don't intersect the non-expanded bounding box. - */ - TEST(GridUniquenessTest, ChoosingLatticePoints_RoundedCorners) - { - const vector>> latticeBasisByModule = { - {{1, 0}, - {0, 1}} - }; - const vector>> domainToPlaneByModule = { - {{1, 0}, - {0, 1}} - }; - - const double readoutResolution = 0.1; - const double smallPadding = (readoutResolution/2) / sqrt(2) - 0.01; - const double largePadding = (readoutResolution/2) / sqrt(2) + 0.01; - const double sideLength = 1 - smallPadding - largePadding; - - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {smallPadding, smallPadding}, {sideLength, sideLength}, - readoutResolution)); - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-1 + largePadding, smallPadding}, {sideLength, sideLength}, - readoutResolution)); - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {-1 + largePadding, -1 + largePadding}, - {sideLength, sideLength}, - readoutResolution)); - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {smallPadding, -1 + largePadding}, - {sideLength, sideLength}, - readoutResolution)); - - ASSERT_FALSE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {largePadding, largePadding}, - {1 - 2*largePadding, 1 - 2*largePadding}, - readoutResolution)); - } - - TEST(GridUniquenessTest, LatticePointsAtHighNumbers) - { - const vector>> domainToPlaneByModule = { - {{1.0}, {0.0}}, - {{0.7071067811865475}, {0.0}}, - {{0.4999999999999999}, {0.0}}, - {{0.3535533905932737}, {0.0}} - }; - - const vector>> latticeBasisByModule = { - {{0.9692380902788696, 0.27146852727840026}, - {0.24612501772995293, 0.9624473173619927}}, - {{0.646673658192444, -0.337238590405595}, - {0.762766792538848, 0.9414192122222954}}, - {{0.9571341168044603, 0.22772704597058202}, - {0.2896450974018824, 0.9737250086823859}}, - {{0.9999617640982025, 0.4924077217967222}, - {0.00874473222065453, 0.8703646566324725}} - }; - - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {3.62142e+09}, - {3.62142e+07}, - 0.2); - }; - - vector> invert2DMatrix(const vector>& M) - { - const double detInv = 1 / (M[0][0]*M[1][1] - M[0][1]*M[1][0]); - return {{detInv*M[1][1], -detInv*M[0][1]}, - {-detInv*M[1][0], detInv*M[0][0]}}; - } - - /** - * Combine a square grid and hexagonal grid to place a common zero at a - * particular point. - * - * Unrotated square and hexagonal lattices have the same lattice points along - * the x axis. Rotate and scale this axis to place the first collision at any - * arbitrary point. - */ - vector>> getPlaneMatrixWithNearestZeroAt( - double x, double y) - { - const double r = sqrt(pow(x, 2) + pow(y, 2)); - - return { - // Square - invert2DMatrix( - {{r, 0}, - {0, r}}), - // Hexagon - invert2DMatrix( - {{r, 0}, - {0, r}}), - }; - } - - vector>> getLatticeBasisWithNearestZeroAt( - double x, double y) - { - double theta = atan2(y, x); - - return { - // Square - {{cos(theta), cos(theta + M_PI/2)}, - {sin(theta), sin(theta + M_PI/2)}}, - // Hexagon - {{cos(theta), cos(theta + M_PI/3)}, - {sin(theta), sin(theta + M_PI/3)}}, - }; - } - - TEST(GridUniquenessTest, ComputeGridUniquenessHypercubeTestPositive) - { - // Zero to the right of the ignored area. - EXPECT_EQ(12, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(12.5, 0.25), - getLatticeBasisWithNearestZeroAt(12.5, 0.25), - 0.01, 0.5).first)); - - // Zero upper-right of the ignored area. - EXPECT_EQ(6, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(6.5, 6.5), - getLatticeBasisWithNearestZeroAt(6.5, 6.5), - 0.01, 0.5).first)); - - // Zero above the ignored area. - EXPECT_EQ(12, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(0.25, 12.5), - getLatticeBasisWithNearestZeroAt(0.25, 12.5), - 0.01, 0.5).first)); - } - - TEST(GridUniquenessTest, ComputeGridUniquenessHypercubeTestNegative) - { - // Zero to the left of the ignored area. - EXPECT_EQ(12, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(-12.5, 0.25), - getLatticeBasisWithNearestZeroAt(-12.5, 0.25), - 0.01, 0.5).first)); - - // Zero upper-left of the ignored area. - EXPECT_EQ(6, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(-6.5, 6.5), - getLatticeBasisWithNearestZeroAt(-6.5, 6.5), - 0.01, 0.5).first)); - - // Zero above the ignored area. - EXPECT_EQ(12, - floor(computeGridUniquenessHypercube( - getPlaneMatrixWithNearestZeroAt(-0.25, 12.5), - getLatticeBasisWithNearestZeroAt(-0.25, 12.5), - 0.01, 0.5).first)); - } - - - /** - * Test 1: Upper right region - * Test 2: Upper left region - * - * Test A: Test an area that's outside both the circle around zero - * on the plane and the circle around zero on the torus. - * - * Test B: Test an area that's outside the circle around zero on the - * plane but that intersects the circle around zero on the torus. - * - * Test C: Test an area that intersects the circle around zero on - * the plane but is outside the circle around zero on the torus. - * - * Test D: Test an area that intersects the circle around zero on - * the plane and the circle around zero on the torus. - * - * (Tests 1C and 2B are impossible.) - */ - TEST(GridUniquenessTest, DistanceOnPlaneNotTorus) - { - const vector>> domainToPlaneByModule = { - {{1, 0}, - {0, 1}} - }; - - const vector>> latticeBasisByModule = { - {{cos(0), cos(M_PI/3)}, - {sin(0), sin(M_PI/3)}} - }; - - const double sideLength = 0.1; - const double readoutResolution = 0.1; - const double r = readoutResolution/2; - - // - // Test 1 - // - const double planeTheta1 = M_PI/6; - const vector planePoint1 = {r*cos(planeTheta1), r*sin(planeTheta1)}; - - const double torusTheta1 = M_PI/4; - const vector torusPoint1 = {r*cos(torusTheta1), r*sin(torusTheta1)}; - const vector torusPoint1OnPlane = { - torusPoint1[0]*cos(0) + torusPoint1[1]*cos(M_PI/3), - torusPoint1[0]*sin(0) + torusPoint1[1]*sin(M_PI/3), - }; - - // Test 1A - const vector corner1A = { - torusPoint1OnPlane[0] + 0.01, - torusPoint1OnPlane[1] + 0.01, - }; - ASSERT_FALSE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner1A, - {sideLength, sideLength}, readoutResolution)); - - // Test 1B - const vector corner1B = { - planePoint1[0] + (torusPoint1OnPlane[0] - planePoint1[0])/2, - planePoint1[1] + (torusPoint1OnPlane[1] - planePoint1[1])/2 - }; - ASSERT_FALSE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner1B, - {sideLength, sideLength}, readoutResolution)); - - // Test 1D - const vector corner1D = { - planePoint1[0] - 0.01, - planePoint1[1] - 0.01, - }; - ASSERT_TRUE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner1D, - {sideLength, sideLength}, readoutResolution)); - - // - // Test 2 - // - const double planeTheta2 = 2*M_PI/3; - - const vector planePoint2 = {r*cos(planeTheta2), r*sin(planeTheta2)}; - - const double torusTheta2 = 3*M_PI/4; - const vector torusPoint2 = {r*cos(torusTheta2), r*sin(torusTheta2)}; - const vector torusPoint2OnPlane = { - torusPoint2[0]*cos(0) + torusPoint2[1]*cos(M_PI/3), - torusPoint2[0]*sin(0) + torusPoint2[1]*sin(M_PI/3), - }; - - // Test 2A - const vector corner2A = { - planePoint2[0] - 0.01 - sideLength, - planePoint2[1] + 0.01, - }; - ASSERT_FALSE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner2A, - {sideLength, sideLength}, readoutResolution)); - - // Test 2C - const vector corner2C = { - planePoint2[0] + (torusPoint2OnPlane[0] - planePoint2[0])/2 - sideLength, - planePoint2[1] + (torusPoint2OnPlane[1] - planePoint2[1])/2 - }; - ASSERT_TRUE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner2C, - {sideLength, sideLength}, readoutResolution)); - - // Test 2D - const vector corner2D = { - torusPoint2OnPlane[0] + 0.01 - sideLength, - torusPoint2OnPlane[1] - 0.01, - }; - ASSERT_TRUE(findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, corner2D, - {sideLength, sideLength}, readoutResolution)); - } - - /** - * A specific failure case that wasn't caught by other unit tests. - */ - TEST(GridUniquenessTest, SpecificRegressionTest1) - { - const vector>> domainToPlaneByModule = { - {{-0.030776, -0.240687, -0.459375},{0.276544, 0.381681, -0.218507}}, - {{0.268763, 0.231442, 0.473435},{-0.408695, 0.427045, 0.0232472}}, - {{0.510017, -0.41195, 0.166473},{0.233775, 0.0332796, -0.633857}}, - {{0.527072, -0.308923, -0.411208},{0.499403, 0.448174, 0.303424}}, - {{0.0695231, 0.92737, -0.0896616},{-0.0773771, 0.0953474, 0.92618}}, - {{0.608297, -0.302698, -0.651094},{-0.697193, -0.453032, -0.440749}},}; - const vector>> latticeBasisByModule = { - {{1, 0.5},{0, 0.866025}}, - {{1, 0.5},{0, 0.866025}}, - {{1, 0.5},{0, 0.866025}}, - {{1, 0.5},{0, 0.866025}}, - {{1, 0.5},{0, 0.866025}}, - {{1, 0.5},{0, 0.866025}},}; - - const double readoutResolution = 0.08; - - ASSERT_TRUE( - findGridCodeZero(domainToPlaneByModule, latticeBasisByModule, - {83.5038, -373.482, 503.877}, - {0.830883, 0.830883, 0.719825}, - readoutResolution)); - } - - TEST(GridUniquenessTest, DeterministicDespiteMultithreading) - { - // This is a previously randomly-generated matrix that triggers multiple - // threads to find zeros, some nearer, some further. Further away zeros - // should not override nearby zeros, and when a thread finds a zero, it - // shouldn't cancel threads that are searching for nearer zeros. - const vector>> domainToPlaneByModule = { - {{0.4088715361390395, -0.9999112506968285, 0.8109731922797785, 0.25590203028822855}, - {-0.9125919498523434, -0.013322564689428938, 0.5850833115066139, 0.9667027210545974}}, - {{-0.704978485994098, -0.016909658985638815, 0.41508560377373277, 0.19893559514770887}, - {-0.05482092926492031, -0.7069045645863304, 0.5724543139323832, 0.6785459667430253}}}; - const vector>> latticeBasisByModule = { - {{1, 0}, - {0, 1}}, - {{1, 0}, - {0, 1}} - }; - - for (int i = 0; i < 100; i++) - { - ASSERT_EQ(0.5, - computeGridUniquenessHypercube(domainToPlaneByModule, - latticeBasisByModule, - 0.2, 0.5).first); - } - } - - TEST(GridUniquenessTest, binSidelengthBasicTest) - { - const vector scales = {1, 2}; - vector>> domainToPlaneByModule; - for (double scale : scales) - { - domainToPlaneByModule.push_back({ - {1/scale, 0}, - {0, 1/scale}, - }); - } - - const double phaseResolution = 0.2; - const double resultPrecision = 0.001; - - const double result = - computeBinSidelength(domainToPlaneByModule, phaseResolution, - resultPrecision); - - const double expected = 2*(scales[0]*phaseResolution/2); - ASSERT_GE(result, expected); - ASSERT_LE(result, expected + resultPrecision); - } - - TEST(GridUniquenessTest, binSidelengthRank1Test) - { - // Each firing field is a band. The first module creates - // horizontal stripes, the second creates vertical stripes. The - // intersection of two bands is a phaseResolution * - // phaseResolution square. - const vector>> domainToPlaneByModule = { - {{0, 1}, - {0, 0}}, - {{1, 0}, - {0, 0}} - }; - - const double phaseResolution = 0.2; - const double resultPrecision = 0.001; - - const double result = - computeBinSidelength(domainToPlaneByModule, phaseResolution, - resultPrecision); - - const double expected = phaseResolution; - - ASSERT_GE(result, expected); - ASSERT_LE(result, expected + resultPrecision); - } - - TEST(GridUniquenessTest, binSidelength1DTest) - { - const vector scales = {1, 2}; - vector>> domainToPlaneByModule; - for (double scale : scales) - { - domainToPlaneByModule.push_back({ - {1/scale}, - {0}, - }); - } - - const double phaseResolution = 0.2; - const double resultPrecision = 0.001; - - const double result = - computeBinSidelength(domainToPlaneByModule, phaseResolution, - resultPrecision); - - const double expected = 2*(scales[0]*phaseResolution/2); - ASSERT_GE(result, expected); - ASSERT_LE(result, expected + resultPrecision); - } -}