Boosting cleanup

src/examples/mnist/MNIST_SP.cpp

@@ -67,7 +67,7 @@ class MNIST {
 
   public:
     UInt verbosity = 1;
-    const UInt train_dataset_iterations = 2u; //epochs somewhat help, at linear time
+    const UInt train_dataset_iterations = 1u; //epochs somewhat help, at linear time
 
 
 void setup() {
@@ -88,7 +88,7 @@ void setup() {
     /* synPermConnected */            0.5f, //no difference, let's leave at 0.5 in the middle
     /* minPctOverlapDutyCycles */     0.2f, //speed of re-learning?
     /* dutyCyclePeriod */             1402,
-    /* boostStrength */               2.0f, // Boosting does help, but entropy is high, on MNIST it does not matter, for learning with TM prefer boosting off (=0.0), or "neutral"=1.0
+    /* boostStrength */               7.0f, // Boosting does help, but entropy is high, on MNIST it does not matter, for learning with TM prefer boosting off (=0.0), or "neutral"=1.0
     /* seed */                        4u,
     /* spVerbosity */                 1u,
     /* wrapAround */                  true); // does not matter (helps slightly)

src/htm/algorithms/Connections.cpp

@@ -395,7 +395,8 @@ void Connections::computeActivity(
     vector<SynapseIdx> &numActiveConnectedSynapsesForSegment,
     const vector<CellIdx> &activePresynapticCells)
 {
-  NTA_ASSERT(numActiveConnectedSynapsesForSegment.size() == segments_.size());
+  NTA_ASSERT(numActiveConnectedSynapsesForSegment.size() == segments_.size()) << "The output vector size must match "<< segments_.size();
+  NTA_ASSERT(activePresynapticCells.size() < segments_.size()) << "Input is larger than our number of segments!";
 
   if( timeseries_ ) {
     // Before each cycle of computation move the currentUpdates to the previous

src/htm/algorithms/Connections.hpp

@@ -391,11 +391,11 @@ class Connections : public Serializable
    * @param numActiveConnectedSynapsesForSegment
    * An output vector for active connected synapse counts per segment.
    *
-   * @param numActivePotentialSynapsesForSegment
+   * @param (optional) numActivePotentialSynapsesForSegment
    * An output vector for active potential synapse counts per segment.
    *
    * @param activePresynapticCells
-   * Active cells in the input.
+   * Active cells in the input as a sparse indices.
    */
   void computeActivity(std::vector<SynapseIdx> &numActiveConnectedSynapsesForSegment,
                        std::vector<SynapseIdx> &numActivePotentialSynapsesForSegment,

src/htm/algorithms/SpatialPooler.cpp

@@ -453,7 +453,7 @@ void SpatialPooler::initialize(
     connections_.raisePermanencesToThreshold( (Segment)i, stimulusThreshold_ );
   }
 
-  updateInhibitionRadius_();
+  inhibitionRadius_ = updateInhibitionRadius_();
 
   if (spVerbosity_ > 0) {
     printParameters();
@@ -465,13 +465,23 @@ void SpatialPooler::initialize(
 void SpatialPooler::compute(const SDR &input, const bool learn, SDR &active) {
   input.reshape(  inputDimensions_ );
   active.reshape( columnDimensions_ );
+
   updateBookeepingVars_(learn);
   calculateOverlap_(input, overlaps_);
+
+  auto &activeVector = active.getSparse();
 
-  boostOverlaps_(overlaps_, boostedOverlaps_);
+  //boosting
+  //must be done before inhibition
+  boostOverlaps_(overlaps_, boostedOverlaps_); //TODO consider removal, @1
 
-  auto &activeVector = active.getSparse();
+  //inhibition
+  //update inhibition radius if it's time, only changes in local inh
+  if(!globalInhibition_ and isUpdateRound_()) {
+    inhibitionRadius_ = updateInhibitionRadius_();
+  }
   inhibitColumns_(boostedOverlaps_, activeVector);
+
   // Notify the active SDR that its internal data vector has changed.  Always
   // call SDR's setter methods even if when modifying the SDR's own data
   // inplace.
@@ -480,13 +490,9 @@ void SpatialPooler::compute(const SDR &input, const bool learn, SDR &active) {
 
   if (learn) {
     adaptSynapses_(input, active);
-    updateDutyCycles_(overlaps_, active);
-    bumpUpWeakColumns_();
-    updateBoostFactors_();
-    if (isUpdateRound_()) {
-      updateInhibitionRadius_();
-      updateMinDutyCycles_();
-    }
+    //boosting
+    bumpUpWeakColumns_(); //TODO suggest removal, low impact, long time
+    updateBoostFactors_(active); //helper for @1, computed after inh
   }
 }
 
@@ -574,11 +580,9 @@ vector<Real> SpatialPooler::initPermanence_(const vector<UInt> &potential, //TOD
 }
 
 
-void SpatialPooler::updateInhibitionRadius_() {
-  if (globalInhibition_) {
-    inhibitionRadius_ =
-        *max_element(columnDimensions_.cbegin(), columnDimensions_.cend());
-    return;
+UInt SpatialPooler::updateInhibitionRadius_() const {
+  if (globalInhibition_) { //always const for global inh
+    return *max_element(columnDimensions_.cbegin(), columnDimensions_.cend());
   }
 
   Real connectedSpan = 0.0f;
@@ -590,7 +594,8 @@ void SpatialPooler::updateInhibitionRadius_() {
   const Real diameter = connectedSpan * columnsPerInput;
   Real radius = (diameter - 1) / 2.0f;
   radius = max((Real)1.0, radius);
-  inhibitionRadius_ = UInt(round(radius));
+
+  return UInt(round(radius));
 }
 
 
@@ -635,26 +640,6 @@ void SpatialPooler::updateMinDutyCyclesLocal_() {
 }
 
 
-void SpatialPooler::updateDutyCycles_(const vector<SynapseIdx> &overlaps,
-                                      SDR &active) {
-
-  // Turn the overlaps array into an SDR. Convert directly to flat-sparse to
-  // avoid copies and  type convertions.
-  SDR newOverlap({ numColumns_ });
-  auto &overlapsSparseVec = newOverlap.getSparse();
-  for (UInt i = 0; i < numColumns_; i++) {
-    if( overlaps[i] != 0 )
-      overlapsSparseVec.push_back( i );
-  }
-  newOverlap.setSparse( overlapsSparseVec );
-
-  const UInt period = std::min(dutyCyclePeriod_, iterationNum_);
-
-  updateDutyCyclesHelper_(overlapDutyCycles_, newOverlap, period);
-  updateDutyCyclesHelper_(activeDutyCycles_, active, period);
-}
-
-
 Real SpatialPooler::avgColumnsPerInput_() const {
   const size_t numDim = max(columnDimensions_.size(), inputDimensions_.size());
   Real columnsPerInput = 0.0f;
@@ -713,11 +698,38 @@ void SpatialPooler::adaptSynapses_(const SDR &input,
 
 void SpatialPooler::bumpUpWeakColumns_() {
   for (UInt i = 0; i < numColumns_; i++) {
+    // skip columns (segments) that are already performing OK
     if (overlapDutyCycles_[i] >= minOverlapDutyCycles_[i]) {
       continue;
     }
+    //bump the weak
     connections_.bumpSegment( i, synPermBelowStimulusInc_ );
   }
+
+  //do updates:
+
+  // update overlap duty cycles (each round)
+  updateDutyCyclesOverlaps_(overlaps_);
+
+  //update minOverlapDutyCycles_ (on update round only) 
+  if (isUpdateRound_()) {
+    updateMinDutyCycles_();
+  }
+}
+
+
+void SpatialPooler::updateDutyCyclesOverlaps_(const vector<SynapseIdx>& overlaps) {
+  SDR newOverlap({ numColumns_ });
+  auto &overlapsSparseVec = newOverlap.getSparse();
+
+  for (UInt i = 0; i < numColumns_; i++) {
+    if( overlaps[i] > Epsilon )
+      overlapsSparseVec.push_back( i );
+    }
+  newOverlap.setSparse( overlapsSparseVec );
+
+  const UInt period = std::min(dutyCyclePeriod_, iterationNum_);
+  updateDutyCyclesHelper_(overlapDutyCycles_, newOverlap, period);
 }
 
 
@@ -743,7 +755,19 @@ void SpatialPooler::updateDutyCyclesHelper_(vector<Real> &dutyCycles,
 }
 
 
-void SpatialPooler::updateBoostFactors_() {
+void SpatialPooler::updateBoostFactors_(const SDR& active) {
+  if(boostStrength_ < htm::Epsilon) return; //skip for disabled boosting
+
+  /**
+    Updates the duty cycles for each column. The ACTIVITY duty cycles is
+    a moving average of the frequency of activation for each column.
+
+    @param active  A SDR of active columns which survived inhibition
+    @param period 
+  */
+  const UInt period = std::min(dutyCyclePeriod_, iterationNum_);
+  updateDutyCyclesHelper_(activeDutyCycles_, active, period);
+
   if (globalInhibition_) {
     updateBoostFactorsGlobal_();
   } else {
@@ -758,21 +782,21 @@ void applyBoosting_(const UInt i,
 		    const Real boost,
 	            vector<Real>& output) {
   if(boost < htm::Epsilon) return; //skip for disabled boosting
-  output[i] = exp((targetDensity - actualDensity[i]) * boost); //TODO doc this code
+  output[i] = exp((targetDensity - actualDensity[i]) * boost); //exponential boosting, default for Numenta
+  //output[i] = log(actualDensity[i]) / log(targetDensity); //log boosting
 }
 
 
 void SpatialPooler::updateBoostFactorsGlobal_() {
-  Real targetDensity;
+  Real targetDensity = localAreaDensity_;
   if (numActiveColumnsPerInhArea_ > 0) {
     UInt inhibitionArea =
-        (UInt)(pow((Real)(2 * inhibitionRadius_ + 1), (Real)columnDimensions_.size()));
+        (UInt)(pow((Real)(2 * inhibitionRadius_ + 1), (Real)columnDimensions_.size())); //FIXME this is broken 
+    // for nD. SP{2000, 1, 1} has ^3 bigger area -> smaller target density
     inhibitionArea = min(inhibitionArea, numColumns_);
     NTA_ASSERT(inhibitionArea > 0);
     targetDensity = ((Real)numActiveColumnsPerInhArea_) / inhibitionArea;
     targetDensity = min(targetDensity, (Real)MAX_LOCALAREADENSITY);
-  } else {
-    targetDensity = localAreaDensity_;
   }
 
   for (UInt i = 0; i < numColumns_; ++i) { 
@@ -830,41 +854,49 @@ void SpatialPooler::inhibitColumns_(const vector<Real> &overlaps,
     density = min(density, (Real)MAX_LOCALAREADENSITY);
   }
 
-  if (globalInhibition_ ||
-      inhibitionRadius_ >
-          *max_element(columnDimensions_.begin(), columnDimensions_.end())) {
+  if (globalInhibition_) {
     inhibitColumnsGlobal_(overlaps, density, activeColumns);
   } else {
     inhibitColumnsLocal_(overlaps, density, activeColumns);
   }
 }
 
+static int missed = 0;
+
 
 void SpatialPooler::inhibitColumnsGlobal_(const vector<Real> &overlaps,
-                                          Real density,
-                                          vector<UInt> &activeColumns) const {
+                                          const Real density,
+                                          SDR_sparse_t &activeColumns) const {
   NTA_ASSERT(!overlaps.empty());
   NTA_ASSERT(density > 0.0f && density <= 1.0f);
 
   activeColumns.clear();
-  const UInt numDesired = (UInt)(density * numColumns_);
+  const UInt numDesired = static_cast<UInt>(density * numColumns_);
   NTA_CHECK(numDesired > 0) << "Not enough columns (" << numColumns_ << ") "
                             << "for desired density (" << density << ").";
   // Sort the columns by the amount of overlap.  First make a list of all of the
   // column indexes.
   activeColumns.reserve(numColumns_);
+
+  int same_overlap = 0;
   for(UInt i = 0; i < numColumns_; i++)
     activeColumns.push_back(i);
   // Compare the column indexes by their overlap.
-  auto compare = [&overlaps](const UInt &a, const UInt &b) -> bool
-    {return (overlaps[a] == overlaps[b]) ? a > b : overlaps[a] > overlaps[b];};  //for determinism if overlaps match (tieBreaker does not solve that),
-  //otherwise we'd return just `return overlaps[a] > overlaps[b]`. 
+  auto compare = [&overlaps, &same_overlap](const UInt &a, const UInt &b) -> bool
+    {
+      if (overlaps[a] == overlaps[b]) {
+	same_overlap++;
+        return  a > b; //but we also need this for deterministic results
+      } else {
+        return overlaps[a] > overlaps[b]; //this is the main "sort columns by overlaps"
+      }
+    };
 
   // Do a partial sort to divide the winners from the losers.  This sort is
   // faster than a regular sort because it stops after it partitions the
   // elements about the Nth element, with all elements on their correct side of
   // the Nth element.
-  std::nth_element(
+    std::nth_element(
     activeColumns.begin(),
     activeColumns.begin() + numDesired,
     activeColumns.end(),
@@ -877,16 +909,25 @@ void SpatialPooler::inhibitColumnsGlobal_(const vector<Real> &overlaps,
   while( !activeColumns.empty() &&
          overlaps[activeColumns.back()] < stimulusThreshold_)
       activeColumns.pop_back();
+  //FIXME not numDesired 
 }
 
 
 void SpatialPooler::inhibitColumnsLocal_(const vector<Real> &overlaps,
-                                         Real density,
+                                         const Real density,
                                          vector<UInt> &activeColumns) const {
+
+  //optimization hack: call faster global inh if radius stretches over the whole input field,
+  //but this should not occur too ofthen because we want to do local inh
+  if(inhibitionRadius_ > *max_element(columnDimensions_.begin(), columnDimensions_.end())) {
+    inhibitColumnsGlobal_(overlaps, density, activeColumns);
+  //slow path normal local inhibition
+  } else {
+
   activeColumns.clear();
 
   // Tie-breaking: when overlaps are equal, columns that have already been
-  // selected are treated as "bigger".
+  // selected are treated as "bigger". //TODO move this idea to the sort/comparison logic
   vector<bool> activeColumnsDense(numColumns_, false);
 
   for (UInt column = 0; column < numColumns_; column++) {
@@ -930,6 +971,7 @@ void SpatialPooler::inhibitColumnsLocal_(const vector<Real> &overlaps,
         activeColumnsDense[column] = true;
       }
   }
+  }
 }