Spatial anomaly

bindings/py/cpp_src/bindings/algorithms/py_SpatialPooler.cpp

@@ -225,6 +225,14 @@ Argument wrapAround boolean value that determines whether or not inputs
         py_SpatialPooler.def("getMinPctOverlapDutyCycles", &SpatialPooler::getMinPctOverlapDutyCycles);
         py_SpatialPooler.def("setMinPctOverlapDutyCycles", &SpatialPooler::setMinPctOverlapDutyCycles);
 
+	py_SpatialPooler.def("anomaly", [](const SpatialPooler& self) {
+			return self.anomaly;
+	});
+	py_SpatialPooler.def("anomalyThreshold", [](const SpatialPooler& self) {
+			return &self.spAnomaly.SPATIAL_TOLERANCE;
+	});
+
+
         // loadFromString
         py_SpatialPooler.def("loadFromString", [](SpatialPooler& self, const py::bytes& inString)
         {
@@ -245,8 +253,8 @@ Argument wrapAround boolean value that determines whether or not inputs
         });
 
         // compute
-        py_SpatialPooler.def("compute", [](SpatialPooler& self, SDR& input, bool learn, SDR& output)
-            { self.compute( input, learn, output ); },
+        py_SpatialPooler.def("compute", [](SpatialPooler& self, const SDR& input, const bool learn, SDR& output, const Real spatialAnomalyInputValue)
+            { self.compute( input, learn, output, spatialAnomalyInputValue ); },
 R"(
 This is the main workhorse method of the SpatialPooler class. This method
 takes an input SDR and computes the set of output active columns. If 'learn' is
@@ -269,10 +277,15 @@ Argument learn A boolean value indicating whether learning should be
 Argument output An SDR representing the winning columns after
         inhibition. The size of the SDR is equal to the number of
         columns (also returned by the method getNumColumns).
+
+Argument spatialAnomalyInputValue (optional) A `Real` value of the original input
+        that was given to encoder. Used only for spatial anomaly computation.
+        See @ref `SP.anomaly`
 )",
         py::arg("input"),
         py::arg("learn") = true,
-        py::arg("output")
+        py::arg("output"),
+	py::arg("spatialAnomalyInputValue") = std::numeric_limits<Real>::min()
         ); 
 
         // setBoostFactors

src/examples/hotgym/HelloSPTP.cpp

@@ -106,8 +106,9 @@ Real64 BenchmarkHotgym::run(UInt EPOCHS, bool useSPlocal, bool useSPglobal, bool
     //Encode
     tEnc.start();
     x+=0.01f; //step size for fn(x)
-    enc.encode(sin(x), input); //model sin(x) function //TODO replace with CSV data
-//    cout << x << "\n" << sin(x) << "\n" << input << "\n\n";
+    const Real value = sin(x);
+    enc.encode(value, input); //model sin(x) function //TODO replace with CSV data
+//    cout << x << "\n" << value << "\n" << input << "\n\n";
     tEnc.stop();
 
     tRng.start();
@@ -117,13 +118,13 @@ Real64 BenchmarkHotgym::run(UInt EPOCHS, bool useSPlocal, bool useSPglobal, bool
     //SP (global x local)
     if(useSPlocal) {
     tSPloc.start();
-    spLocal.compute(input, true, outSPlocal);
+    spLocal.compute(input, true, outSPlocal, value /* optional for spatial anomaly*/);
     tSPloc.stop();
     }
 
     if(useSPglobal) {
     tSPglob.start();
-    spGlobal.compute(input, true, outSPglobal);
+    spGlobal.compute(input, true, outSPglobal, value /* optional for spatial anomaly */);
     tSPglob.stop();
     }
     outSP = outSPglobal; //toggle if local/global SP is used further down the chain (TM, Anomaly)
@@ -168,6 +169,7 @@ Real64 BenchmarkHotgym::run(UInt EPOCHS, bool useSPlocal, bool useSPglobal, bool
       // output values
       cout << "Epoch = " << e << endl;
       cout << "Anomaly = " << an << endl;
+      cout << "Anomaly (spatial) = " << spGlobal.anomaly << endl;
       cout << "Anomaly (avg) = " << avgAnom10.getCurrentAvg() << endl;
       cout << "Anomaly (Likelihood) = " << anLikely << endl;
       cout << "SP (g)= " << outSP << endl;
@@ -217,10 +219,12 @@ Real64 BenchmarkHotgym::run(UInt EPOCHS, bool useSPlocal, bool useSPglobal, bool
         if(useSPglobal) { NTA_CHECK(outSPglobal == goldSP) << "Deterministic output of SP (g) failed!\n" << outSP << "should be:\n" << goldSP; }
         if(useSPlocal) {  NTA_CHECK(outSPlocal == goldSPlocal) << "Deterministic output of SP (l) failed!\n" << outSPlocal << "should be:\n" << goldSPlocal; }
         if(useTM) {       NTA_CHECK(outTM == goldTM) << "Deterministic output of TM failed!\n" << outTM << "should be:\n" << goldTM; }
-        NTA_CHECK(static_cast<UInt>(an *10000.0f) == static_cast<UInt>(goldAn *10000.0f)) //compare to 4 decimal places
+        // anomalies
+	NTA_CHECK(static_cast<UInt>(an *10000.0f) == static_cast<UInt>(goldAn *10000.0f)) //compare to 4 decimal places
 		               << "Deterministic output of Anomaly failed! " << an << "should be: " << goldAn;
 	NTA_CHECK(static_cast<UInt>(avgAnom10.getCurrentAvg() * 10000.0f) == static_cast<UInt>(goldAnAvg * 10000.0f)) 
 		<< "Deterministic average anom score failed:" << avgAnom10.getCurrentAvg() << " should be: " << goldAnAvg;
+	if(useSPglobal) { NTA_CHECK(0.0f == spGlobal.anomaly) << "Deterministic spatial anomaly mismatch!" << spGlobal.anomaly; }
       }
 
       // check runtime speed

src/htm/algorithms/SpatialPooler.cpp

@@ -459,7 +459,11 @@ void SpatialPooler::initialize(
 }
 
 
-void SpatialPooler::compute(const SDR &input, const bool learn, SDR &active) {
+void SpatialPooler::compute(const SDR &input, 
+		            const bool learn, 
+			    SDR &active,
+			    const Real spatialAnomalyInputValue) {
+
   input.reshape(  inputDimensions_ );
   active.reshape( columnDimensions_ );
   updateBookeepingVars_(learn);
@@ -485,6 +489,11 @@ void SpatialPooler::compute(const SDR &input, const bool learn, SDR &active) {
       updateMinDutyCycles_();
     }
   }
+
+  //update spatial anomaly
+  if(this->spAnomaly.enabled) {
+    spAnomaly.compute(spatialAnomalyInputValue);
+  }
 }
 
 
@@ -1018,6 +1027,10 @@ bool SpatialPooler::operator==(const SpatialPooler& o) const{
 
   //Random
   if (rng_ != o.rng_) return false;
+
+  //spatial anomaly
+  if (spAnomaly != o.spAnomaly) return false;
+
   return true;
 
 }

src/htm/algorithms/SpatialPooler.hpp

@@ -233,8 +233,15 @@ class SpatialPooler : public Serializable
   @param active An SDR representing the winning columns after
         inhibition. The size of the SDR is equal to the number of
         columns (also returned by the method getNumColumns).
+
+  @param  spatialAnomalyInputValue (optional) `Real` used for computing "spatial anomaly", see @ref `this.spatial_anomaly.compute()`,
+          obtained by @ref `SP.anomaly`
+
    */
-  virtual void compute(const SDR &input, const bool learn, SDR &active);
+  virtual void compute(const SDR &input, 
+		       const bool learn, 
+		       SDR &active, 
+		       const Real spatialAnomalyInputValue = std::numeric_limits<Real>::min());
 
 
   /**
@@ -276,7 +283,11 @@ class SpatialPooler : public Serializable
        CEREAL_NVP(synPermBelowStimulusInc_),
        CEREAL_NVP(synPermConnected_),
        CEREAL_NVP(minPctOverlapDutyCycles_),
-       CEREAL_NVP(wrapAround_));
+       CEREAL_NVP(wrapAround_),
+       CEREAL_NVP(spAnomaly.minVal_),
+       CEREAL_NVP(spAnomaly.maxVal_),
+       CEREAL_NVP(spAnomaly.anomalyScore_)
+       );
     ar(CEREAL_NVP(boostFactors_));
     ar(CEREAL_NVP(overlapDutyCycles_));
     ar(CEREAL_NVP(activeDutyCycles_));
@@ -309,7 +320,11 @@ class SpatialPooler : public Serializable
        CEREAL_NVP(synPermBelowStimulusInc_),
        CEREAL_NVP(synPermConnected_),
        CEREAL_NVP(minPctOverlapDutyCycles_),
-       CEREAL_NVP(wrapAround_));
+       CEREAL_NVP(wrapAround_),
+       CEREAL_NVP(spAnomaly.minVal_),
+       CEREAL_NVP(spAnomaly.maxVal_),
+       CEREAL_NVP(spAnomaly.anomalyScore_)
+       );
     ar(CEREAL_NVP(boostFactors_));
     ar(CEREAL_NVP(overlapDutyCycles_));
     ar(CEREAL_NVP(activeDutyCycles_));
@@ -1209,6 +1224,89 @@ class SpatialPooler : public Serializable
   Random rng_;
 
 public:
+  /**
+   *  holds together functionality for computing
+   *  `spatial anomaly`. This is used in NAB.
+   */
+  struct spatial_anomaly {
+
+    /** Fraction outside of the range of values seen so far that will be considered
+     * a spatial anomaly regardless of the anomaly likelihood calculation. 
+     * This accounts for the human labelling bias for spatial values larger than what
+     * has been seen so far.
+     * Default value 0.05f aka 5%, as used in NAB.
+     */
+    Real SPATIAL_TOLERANCE = 0.05;
+
+    /**
+     * toggle whether we should compute spatial anomaly
+     * Default true.
+     */
+    bool enabled = true;
+
+    /**
+     *  compute if the current input `value` is considered spatial-anomaly.
+     *  update internal variables.
+     *
+     *  @param value Real, input #TODO currently handles only 1 variable input, and requires value passed to compute! 
+     *  # later remove, and implement using SP's internal state. But for now we are compatible with NAB's implementation. 
+     *
+     * @return nothing, but updates internal variable `anomalyScore_`, which is either `NO_ANOMALY` (0.0f) , 
+     *   or `SPATIAL_ANOMALY` (exactly 0.9995947141f). Accessed by public @ref `SP.anomaly`.
+     *
+     */
+    void compute(const Real value) {
+      anomalyScore_ = NO_ANOMALY;
+      if(not enabled) return;
+      NTA_CHECK(SPATIAL_TOLERANCE >= 0.0f and SPATIAL_TOLERANCE <= 1.0f); 
+
+      if(minVal_ != maxVal_) {
+        const Real tolerance = (maxVal_ - minVal_) * SPATIAL_TOLERANCE;
+        const Real maxExpected = maxVal_ + tolerance;
+        const Real minExpected = minVal_ - tolerance;
+
+        if(value > maxExpected or value < minExpected) { //spatial anomaly
+          anomalyScore_ = SPATIAL_ANOMALY; 
+	}
+      }
+      if(value > maxVal_ or maxVal_ == INF_) maxVal_ = value;
+      if(value < minVal_ or minVal_ == INF_) minVal_ = value;
+      NTA_ASSERT(anomalyScore_ == NO_ANOMALY or anomalyScore_ == SPATIAL_ANOMALY) << "Out of bounds of acceptable values";
+    }
+
+    bool operator==(const spatial_anomaly& o) const noexcept {
+      return minVal_ == o.minVal_ and
+	     maxVal_ == o.maxVal_ and
+	     anomalyScore_ == o.anomalyScore_ and
+	     enabled == o.enabled and
+	     SPATIAL_TOLERANCE == o.SPATIAL_TOLERANCE;
+    }
+
+    inline bool operator!=(const spatial_anomaly& o) const noexcept {
+      return !this->operator==(o);
+    }
+
+    private:
+      friend class SpatialPooler;
+      static const constexpr Real INF_ = std::numeric_limits<Real>::infinity();
+      Real minVal_ = INF_;
+      Real maxVal_ = INF_;
+      Real anomalyScore_ = NO_ANOMALY; //default score = no anomaly
+
+    public:
+      static const constexpr Real NO_ANOMALY = 0.0f;
+      static const constexpr Real SPATIAL_ANOMALY = 0.9995947141f; //almost 1.0 = max anomaly. Encodes value specific to spatial anomaly (so this can be recognized on results),
+        // "5947141" would translate in l33t speech to "spatial" :)
+  } spAnomaly;
+
+  /**
+   *  spatial anomaly
+   *
+   *  updated on each @ref `compute()`. 
+   *
+   *  @return either 0.0f (no anomaly), or exactly 0.9995947141f (spatial anomaly). This specific value can be recognized in results. 
+   */
+  const Real& anomaly = spAnomaly.anomalyScore_;
   const Connections &connections = connections_;
 };
 

src/test/unit/algorithms/SpatialPoolerTest.cpp

@@ -32,6 +32,7 @@
 #include <htm/types/Types.hpp>
 #include <htm/utils/Log.hpp>
 #include <htm/os/Timer.hpp>
+#include <htm/encoders/ScalarEncoder.hpp>
 
 namespace testing {
 
@@ -2103,5 +2104,70 @@ TEST(SpatialPoolerTest, ExactOutput) {
   ASSERT_EQ( columns, gold_sdr );
 }
 
+TEST(SpatialPoolerTest, spatialAnomaly) {
+  SDR inputs({ 1000 });
+  SDR columns({ 200 });
+
+  {
+  SpatialPooler sp(inputs.dimensions, columns.dimensions);
+  // test too large threshold
+  sp.spAnomaly.SPATIAL_TOLERANCE = 1.2345f; //out of bounds, will crash!
+  EXPECT_ANY_THROW(sp.compute(inputs, false, columns, 0.1f /*whatever, fails on TOLERANCE */)) << "Spatial anomaly should fail if SPATIAL_TOLERANCE is out of bounds!";
+  sp.spAnomaly.SPATIAL_TOLERANCE = 0.01f; //within bounds, OK
+  EXPECT_NO_THROW(sp.compute(inputs, false, columns, 0.1f /*whatever */));
+  }
+
+  {
+  SpatialPooler sp({inputs.dimensions}, {columns.dimensions},
+    /*potentialRadius*/ 99999,
+    /*potentialPct*/ 0.5f,
+    /*globalInhibition*/ true,
+    /*localAreaDensity*/ 0.05f,
+    /*stimulusThreshold*/ 3u,
+    /*synPermInactiveDec*/ 0.008f,
+    /*synPermActiveInc*/ 0.05f,
+    /*synPermConnected*/ 0.1f,
+    /*minPctOverlapDutyCycles*/ 0.001f,
+    /*dutyCyclePeriod*/ 200,
+    /*boostStrength*/ 10.0f,
+    /*seed*/ 42,
+    /*spVerbosity*/ 0,
+    /*wrapAround*/ true);
+
+  //test spatial anomaly computation
+  sp.spAnomaly.SPATIAL_TOLERANCE = 0.2f; //threshold 20%
+  ScalarEncoderParameters params;
+  params.minimum = 0.0f;
+  params.maximum = 100.0f;
+  params.size = 1000;
+  params.sparsity = 0.3f;
+  ScalarEncoder enc(params);
+
+  Real val;
+
+  val = 0.0f;
+  enc.encode(val, inputs); //TODO can SDR hold .origValue = Real which encoders would set? Classifier,Predictor, and spatia_anomaly would use that
+  sp.compute(inputs, true, columns, val);
+  EXPECT_EQ(0.0f, sp.anomaly);
+  EXPECT_EQ(sp.spAnomaly.NO_ANOMALY, sp.anomaly) << "should be the same as above";
+
+  val = 10.0f;
+  enc.encode(val, inputs); 
+  sp.compute(inputs, true, columns, val);
+  EXPECT_EQ(sp.spAnomaly.NO_ANOMALY, sp.anomaly);
+
+  val = 11.99f; //(10-0) * 0.2 == 2 -> <-2, +12> is not anomalous
+  enc.encode(val, inputs); 
+  sp.compute(inputs, true, columns, val);
+  EXPECT_EQ(0.0f, sp.anomaly) << "This shouldn't be an anomaly!";
+
+  val = 100.0f; //(12-0) * 0.2 == ~2.2 -> <-2.2, +14.2> is not anomalous, but 100 is!
+  enc.encode(val, inputs);
+  sp.compute(inputs, true, columns, val);
+  EXPECT_EQ(0.9995947141f, sp.anomaly) << "This should be an anomaly!";
+  EXPECT_EQ(sp.spAnomaly.SPATIAL_ANOMALY, sp.anomaly) << "Should be same as above!";
+  }
+
+}
 
 } // end anonymous namespace