Merge pull request #924 from suku248/numerical_imprecision_stdp

Numerical imprecision stdp

models/stdp_connection.h

@@ -250,10 +250,9 @@ STDPConnection< targetidentifierT >::send( Event& e,
   {
     minus_dt = t_lastspike - ( start->t_ + dendritic_delay );
     ++start;
-    if ( minus_dt == 0 )
-    {
-      continue;
-    }
+    // get_history() should make sure that
+    // start->t_ > t_lastspike - dendritic_delay, i.e. minus_dt < 0
+    assert( minus_dt < -1.0 * kernel().connection_manager.get_stdp_eps() );
     weight_ = facilitate_( weight_, Kplus_ * std::exp( minus_dt / tau_plus_ ) );
   }
 

models/stdp_connection_hom.h

@@ -302,10 +302,9 @@ STDPConnectionHom< targetidentifierT >::send( Event& e,
   {
     minus_dt = t_lastspike - ( start->t_ + dendritic_delay );
     ++start;
-    if ( minus_dt == 0 )
-    {
-      continue;
-    }
+    // get_history() should make sure that
+    // start->t_ > t_lastspike - dendritic_delay, i.e. minus_dt < 0
+    assert( minus_dt < -1.0 * kernel().connection_manager.get_stdp_eps() );
     weight_ =
       facilitate_( weight_, Kplus_ * std::exp( minus_dt / cp.tau_plus_ ), cp );
   }

models/stdp_dopa_connection.h

@@ -452,7 +452,8 @@ STDPDopaConnection< targetidentifierT >::process_dopa_spikes_(
   // process dopa spikes in (t0, t1]
   // propagate weight from t0 to t1
   if ( ( dopa_spikes.size() > dopa_spikes_idx_ + 1 )
-    && ( dopa_spikes[ dopa_spikes_idx_ + 1 ].spike_time_ <= t1 ) )
+    && ( t1 - dopa_spikes[ dopa_spikes_idx_ + 1 ].spike_time_ > -1.0
+           * kernel().connection_manager.get_stdp_eps() ) )
   {
     // there is at least 1 dopa spike in (t0, t1]
     // propagate weight up to first dopa spike and update dopamine trace
@@ -468,7 +469,8 @@ STDPDopaConnection< targetidentifierT >::process_dopa_spikes_(
     // process remaining dopa spikes in (t0, t1]
     double cd;
     while ( ( dopa_spikes.size() > dopa_spikes_idx_ + 1 )
-      && ( dopa_spikes[ dopa_spikes_idx_ + 1 ].spike_time_ <= t1 ) )
+      && ( t1 - dopa_spikes[ dopa_spikes_idx_ + 1 ].spike_time_ > -1.0
+                * kernel().connection_manager.get_stdp_eps() ) )
     {
       // propagate weight up to next dopa spike and update dopamine trace
       // weight and dopamine trace n are at time of last dopa spike td but
@@ -565,8 +567,9 @@ STDPDopaConnection< targetidentifierT >::send( Event& e,
     process_dopa_spikes_( dopa_spikes, t0, start->t_ + dendritic_delay, cp );
     t0 = start->t_ + dendritic_delay;
     minus_dt = t_last_update_ - t0;
-    if ( start->t_ < t_spike ) // only depression if pre- and postsyn. spike
-                               // occur at the same time
+    // facilitate only in case of post- after presyn. spike
+    // skip facilitation if pre- and postsyn. spike occur at the same time
+    if ( t_spike - start->t_ > kernel().connection_manager.get_stdp_eps() )
     {
       facilitate_( Kplus_ * std::exp( minus_dt / cp.tau_plus_ ), cp );
     }

models/stdp_pl_connection_hom.h

@@ -257,10 +257,9 @@ STDPPLConnectionHom< targetidentifierT >::send( Event& e,
   {
     minus_dt = t_lastspike - ( start->t_ + dendritic_delay );
     start++;
-    if ( minus_dt == 0 )
-    {
-      continue;
-    }
+    // get_history() should make sure that
+    // start->t_ > t_lastspike - dendritic_delay, i.e. minus_dt < 0
+    assert( minus_dt < -1.0 * kernel().connection_manager.get_stdp_eps() );
     weight_ = facilitate_(
       weight_, Kplus_ * std::exp( minus_dt * cp.tau_plus_inv_ ), cp );
   }

models/stdp_triplet_connection.h

@@ -259,11 +259,9 @@ STDPTripletConnection< targetidentifierT >::send( Event& e,
     // Pfister et al, 2006
     double ky = start->triplet_Kminus_ - 1.0;
     ++start;
-    if ( minus_dt == 0 )
-    {
-      continue;
-    }
-
+    // get_history() should make sure that
+    // start->t_ > t_lastspike - dendritic_delay, i.e. minus_dt < 0
+    assert( minus_dt < -1.0 * kernel().connection_manager.get_stdp_eps() );
     weight_ =
       facilitate_( weight_, Kplus_ * std::exp( minus_dt / tau_plus_ ), ky );
   }

models/vogels_sprekeler_connection.h

@@ -215,10 +215,9 @@ VogelsSprekelerConnection< targetidentifierT >::send( Event& e,
   {
     minus_dt = t_lastspike - ( start->t_ + dendritic_delay );
     ++start;
-    if ( minus_dt == 0 )
-    {
-      continue;
-    }
+    // get_history() should make sure that
+    // start->t_ > t_lastspike - dendritic_delay, i.e. minus_dt < 0
+    assert( minus_dt < -1.0 * kernel().connection_manager.get_stdp_eps() );
     weight_ = facilitate_( weight_, Kplus_ * std::exp( minus_dt / tau_ ) );
   }
 

nestkernel/archiving_node.cpp

@@ -29,6 +29,9 @@
 
 #include "archiving_node.h"
 
+// Includes from nestkernel:
+#include "kernel_manager.h"
+
 // Includes from sli:
 #include "dictutils.h"
 
@@ -81,7 +84,9 @@ Archiving_Node::register_stdp_connection( double t_first_read )
   // For details see bug #218. MH 08-04-22
 
   for ( std::deque< histentry >::iterator runner = history_.begin();
-        runner != history_.end() && runner->t_ <= t_first_read;
+        runner != history_.end()
+          && ( t_first_read - runner->t_ > -1.0
+                 * kernel().connection_manager.get_stdp_eps() );
         ++runner )
   {
     ( runner->access_counter_ )++;
@@ -100,7 +105,7 @@ nest::Archiving_Node::get_K_value( double t )
   int i = history_.size() - 1;
   while ( i >= 0 )
   {
-    if ( t > history_[ i ].t_ )
+    if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() )
     {
       return ( history_[ i ].Kminus_
         * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) );
@@ -126,7 +131,7 @@ nest::Archiving_Node::get_K_values( double t,
   int i = history_.size() - 1;
   while ( i >= 0 )
   {
-    if ( t > history_[ i ].t_ )
+    if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() )
     {
       triplet_K_value = ( history_[ i ].triplet_Kminus_
         * std::exp( ( history_[ i ].t_ - t ) * tau_minus_triplet_inv_ ) );
@@ -159,12 +164,16 @@ nest::Archiving_Node::get_history( double t1,
   else
   {
     std::deque< histentry >::iterator runner = history_.begin();
-    while ( ( runner != history_.end() ) && ( runner->t_ <= t1 ) )
+    while ( ( runner != history_.end() )
+      && ( t1 - runner->t_ > -1.0
+                * kernel().connection_manager.get_stdp_eps() ) )
     {
       ++runner;
     }
     *start = runner;
-    while ( ( runner != history_.end() ) && ( runner->t_ <= t2 ) )
+    while ( ( runner != history_.end() )
+      && ( t2 - runner->t_ > -1.0
+                * kernel().connection_manager.get_stdp_eps() ) )
     {
       ( runner->access_counter_ )++;
       ++runner;

nestkernel/connection_manager.cpp

@@ -28,6 +28,8 @@
 // C++ includes:
 #include <cassert>
 #include <cmath>
+#include <iomanip>
+#include <limits>
 #include <set>
 #include <vector>
 
@@ -77,6 +79,7 @@ nest::ConnectionManager::ConnectionManager()
   , initial_connector_capacity_( CONFIG_CONNECTOR_CUTOFF )
   , large_connector_limit_( CONFIG_CONNECTOR_CUTOFF * 2 )
   , large_connector_growth_factor_( 1.5 )
+  , stdp_eps_( 1.0e-6 )
 {
 }
 
@@ -1450,3 +1453,31 @@ nest::ConnectionManager::get_targets( const std::vector< index >& sources,
     }
   }
 }
+
+void
+nest::ConnectionManager::set_stdp_eps( const double stdp_eps )
+{
+  if ( not( stdp_eps < Time::get_resolution().get_ms() ) )
+  {
+    throw KernelException(
+      "The epsilon used for spike-time comparison in STDP must be less "
+      "than the simulation resolution." );
+  }
+  else if ( stdp_eps < 0 )
+  {
+    throw KernelException(
+      "The epsilon used for spike-time comparison in STDP must not be "
+      "negative." );
+  }
+  else
+  {
+    stdp_eps_ = stdp_eps;
+
+    std::ostringstream os;
+    os << "Epsilon for spike-time comparison in STDP was set to "
+       << std::setprecision( std::numeric_limits< long double >::digits10 )
+       << stdp_eps_ << ".";
+
+    LOG( M_INFO, "ConnectionManager::set_stdp_eps", os.str() );
+  }
+}

nestkernel/connection_manager.h

@@ -313,6 +313,10 @@ class ConnectionManager : public ManagerInterface
    */
   double get_large_connector_growth_factor() const;
 
+  double get_stdp_eps() const;
+
+  void set_stdp_eps( const double stdp_eps );
+
 private:
   /**
    * Update delay extrema to current values.
@@ -419,6 +423,10 @@ class ConnectionManager : public ManagerInterface
 
   //! Capacity growth factor to use beyond the limit
   double large_connector_growth_factor_;
+
+  //! Maximum distance between (double) spike times in STDP that is
+  //! still considered 0. See issue #894
+  double stdp_eps_;
 };
 
 inline DictionaryDatum&
@@ -457,6 +465,12 @@ ConnectionManager::get_large_connector_growth_factor() const
   return large_connector_growth_factor_;
 }
 
+inline double
+ConnectionManager::get_stdp_eps() const
+{
+  return stdp_eps_;
+}
+
 } // namespace nest
 
 #endif /* CONNECTION_MANAGER_H */

nestkernel/nestmodule.cpp

@@ -1699,6 +1699,26 @@ NestModule::DisableStructuralPlasticity_Function::execute(
   i->EStack.pop();
 }
 
+/**
+ * Set epsilon that is used for comparing spike times in STDP.
+ * Spike times in STDP synapses are currently represented as double
+ * values. The epsilon defines the maximum distance between spike
+ * times that is still considered 0.
+ *
+ * Note: See issue #894
+ */
+void
+NestModule::SetStdpEps_dFunction::execute( SLIInterpreter* i ) const
+{
+  i->assert_stack_load( 1 );
+  const double stdp_eps = getValue< double >( i->OStack.top() );
+
+  kernel().connection_manager.set_stdp_eps( stdp_eps );
+
+  i->OStack.pop();
+  i->EStack.pop();
+}
+
 void
 NestModule::init( SLIInterpreter* i )
 {
@@ -1794,6 +1814,9 @@ NestModule::init( SLIInterpreter* i )
     "GetStructuralPlasticityStatus", &getstructuralplasticitystatus_function );
   i->createcommand( "Disconnect", &disconnect_i_i_lfunction );
   i->createcommand( "Disconnect_g_g_D_D", &disconnect_g_g_D_Dfunction );
+
+  i->createcommand( "SetStdpEps", &setstdpeps_dfunction );
+
   // Add connection rules
   kernel().connection_manager.register_conn_builder< OneToOneBuilder >(
     "one_to_one" );

nestkernel/nestmodule.h

@@ -441,6 +441,12 @@ class NestModule : public SLIModule
     void execute( SLIInterpreter* ) const;
   } disablestructuralplasticity_function;
 
+  class SetStdpEps_dFunction : public SLIFunction
+  {
+  public:
+    void execute( SLIInterpreter* ) const;
+  } setstdpeps_dfunction;
+
   //@}
 };