Fix postsynaptic trace edge case

models/clopath_connection.h

@@ -155,7 +155,7 @@ class ClopathConnection : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

models/stdp_connection.h

@@ -167,7 +167,7 @@ class STDPConnection : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void
@@ -220,7 +220,7 @@ STDPConnection< targetidentifierT >::send( Event& e,
   const CommonSynapseProperties& )
 {
   // synapse STDP depressing/facilitation dynamics
-  double t_spike = e.get_stamp().get_ms();
+  const double t_spike = e.get_stamp().get_ms();
 
   // use accessor functions (inherited from Connection< >) to obtain delay and
   // target
@@ -255,9 +255,8 @@ STDPConnection< targetidentifierT >::send( Event& e,
     weight_ = facilitate_( weight_, Kplus_ * std::exp( minus_dt / tau_plus_ ) );
   }
 
-  // depression due to new pre-synaptic spike
-  weight_ =
-    depress_( weight_, target->get_K_value( t_spike - dendritic_delay ) );
+  const double _K_value = target->get_K_value( t_spike - dendritic_delay );
+  weight_ = depress_( weight_, _K_value );
 
   e.set_receiver( *target );
   e.set_weight( weight_ );

models/stdp_connection_facetshw_hom.h

@@ -292,7 +292,7 @@ class STDPFACETSHWConnectionHom : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

models/stdp_connection_hom.h

@@ -219,7 +219,7 @@ class STDPConnectionHom : public Connection< targetidentifierT >
     ConnTestDummyNode dummy_target;
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
 private:

models/stdp_dopa_connection.h

@@ -262,7 +262,7 @@ class STDPDopaConnection : public Connection< targetidentifierT >
     ConnTestDummyNode dummy_target;
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

models/stdp_pl_connection_hom.h

@@ -188,7 +188,7 @@ class STDPPLConnectionHom : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

models/stdp_triplet_connection.h

@@ -183,7 +183,7 @@ class STDPTripletConnection : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

models/vogels_sprekeler_connection.h

@@ -141,7 +141,7 @@ class VogelsSprekelerConnection : public Connection< targetidentifierT >
 
     ConnectionBase::check_connection_( dummy_target, s, t, receptor_type );
 
-    t.register_stdp_connection( t_lastspike_ - get_delay() );
+    t.register_stdp_connection( t_lastspike_ - get_delay(), get_delay() );
   }
 
   void

nestkernel/archiving_node.cpp

@@ -48,6 +48,8 @@ nest::Archiving_Node::Archiving_Node()
   , tau_minus_inv_( 1. / tau_minus_ )
   , tau_minus_triplet_( 110.0 )
   , tau_minus_triplet_inv_( 1. / tau_minus_triplet_ )
+  , max_delay_( -1.0 )
+  , trace_( 0. )
   , last_spike_( -1.0 )
   , Ca_t_( 0.0 )
   , Ca_minus_( 0.0 )
@@ -66,6 +68,8 @@ nest::Archiving_Node::Archiving_Node( const Archiving_Node& n )
   , tau_minus_inv_( n.tau_minus_inv_ )
   , tau_minus_triplet_( n.tau_minus_triplet_ )
   , tau_minus_triplet_inv_( n.tau_minus_triplet_inv_ )
+  , max_delay_( n.max_delay_ )
+  , trace_( n.trace_ )
   , last_spike_( n.last_spike_ )
   , Ca_t_( n.Ca_t_ )
   , Ca_minus_( n.Ca_minus_ )
@@ -76,7 +80,7 @@ nest::Archiving_Node::Archiving_Node( const Archiving_Node& n )
 }
 
 void
-Archiving_Node::register_stdp_connection( double t_first_read )
+Archiving_Node::register_stdp_connection( double t_first_read, double delay )
 {
   // Mark all entries in the deque, which we will not read in future as read by
   // this input input, so that we savely increment the incoming number of
@@ -93,26 +97,38 @@ Archiving_Node::register_stdp_connection( double t_first_read )
   }
 
   n_incoming_++;
+
+  max_delay_ = std::max( delay, max_delay_ );
 }
 
 double
 nest::Archiving_Node::get_K_value( double t )
 {
+  // case when the neuron has not yet spiked
   if ( history_.empty() )
   {
-    return Kminus_;
+    trace_ = 0.;
+    return trace_;
   }
+
+  // search for the latest post spike in the history buffer that came strictly
+  // before `t`
   int i = history_.size() - 1;
   while ( i >= 0 )
   {
     if ( t - history_[ i ].t_ > kernel().connection_manager.get_stdp_eps() )
     {
-      return ( history_[ i ].Kminus_
+      trace_ = ( history_[ i ].Kminus_
         * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) );
+      return trace_;
     }
-    i--;
+    --i;
   }
-  return 0;
+
+  // this case occurs when the trace was requested at a time precisely at or
+  // before the first spike in the history
+  trace_ = 0.;
+  return trace_;
 }
 
 void
@@ -127,7 +143,9 @@ nest::Archiving_Node::get_K_values( double t,
     K_value = Kminus_;
     return;
   }
-  // case
+
+  // search for the latest post spike in the history buffer that came strictly
+  // before `t`
   int i = history_.size() - 1;
   while ( i >= 0 )
   {
@@ -139,12 +157,11 @@ nest::Archiving_Node::get_K_values( double t,
         * std::exp( ( history_[ i ].t_ - t ) * tau_minus_inv_ ) );
       return;
     }
-    i--;
+    --i;
   }
 
-  // we only get here if t< time of all spikes in history)
-
-  // return 0.0 for both K values
+  // this case occurs when the trace was requested at a time precisely at or
+  // before the first spike in the history
   triplet_K_value = 0.0;
   K_value = 0.0;
 }
@@ -187,10 +204,17 @@ nest::Archiving_Node::set_spiketime( Time const& t_sp, double offset )
   if ( n_incoming_ )
   {
     // prune all spikes from history which are no longer needed
-    // except the penultimate one. we might still need it.
+    // only remove a spike if:
+    // - its access counter indicates it has been read out by all connected
+    //   STDP synapses, and
+    // - there is another, later spike, that is strictly more than
+    //   (max_delay_ + eps) away from the new spike (at t_sp_ms)
     while ( history_.size() > 1 )
     {
-      if ( history_.front().access_counter_ >= n_incoming_ )
+      const double next_t_sp = history_[ 1 ].t_;
+      if ( history_.front().access_counter_ >= n_incoming_
+        && t_sp_ms - next_t_sp > max_delay_
+            + kernel().connection_manager.get_stdp_eps() )
       {
         history_.pop_front();
       }
@@ -226,6 +250,7 @@ nest::Archiving_Node::get_status( DictionaryDatum& d ) const
   def< double >( d, names::tau_Ca, tau_Ca_ );
   def< double >( d, names::beta_Ca, beta_Ca_ );
   def< double >( d, names::tau_minus_triplet, tau_minus_triplet_ );
+  def< double >( d, names::post_trace, trace_ );
 #ifdef DEBUG_ARCHIVER
   def< int >( d, names::archiver_length, history_.size() );
 #endif

nestkernel/archiving_node.h

@@ -32,6 +32,7 @@
 #define ARCHIVING_NODE_H
 
 // C++ includes:
+#include <algorithm>
 #include <deque>
 
 // Includes from nestkernel:
@@ -90,8 +91,6 @@ class Archiving_Node : public Node
    * \fn int get_synaptic_elements_vacant(Name n)
    * Get the number of synaptic elements of type n which are available
    * for new synapse creation
-   * Returns a negative number to indicate that synaptic elements
-   * must be deleted during the next update
    */
   int get_synaptic_elements_vacant( Name n ) const;
 
@@ -130,7 +129,9 @@ class Archiving_Node : public Node
 
   /**
    * \fn double get_K_value(long t)
-   * return the Kminus value at t (in ms).
+   * return the Kminus (synaptic trace) value at t (in ms). When the trace is
+   * requested at the exact same time that the neuron emits a spike, the trace
+   * value as it was just before the spike is returned.
    */
   double get_K_value( double t );
 
@@ -167,7 +168,7 @@ class Archiving_Node : public Node
    * t_first_read: The newly registered synapse will read the history entries
    * with t > t_first_read.
    */
-  void register_stdp_connection( double t_first_read );
+  void register_stdp_connection( double t_first_read, double delay );
 
   void get_status( DictionaryDatum& d ) const;
   void set_status( const DictionaryDatum& d );
@@ -216,6 +217,9 @@ class Archiving_Node : public Node
   double tau_minus_triplet_;
   double tau_minus_triplet_inv_;
 
+  double max_delay_;
+  double trace_;
+
   double last_spike_;
 
   // spiking history needed by stdp synapses

nestkernel/nest_names.cpp

@@ -313,7 +313,9 @@ const Name port_name( "port_name" );
 const Name port_width( "port_width" );
 const Name ports( "ports" );
 const Name post_synaptic_element( "post_synaptic_element" );
+const Name post_trace( "post_trace" );
 const Name pre_synaptic_element( "pre_synaptic_element" );
+const Name pre_trace( "pre_trace" );
 const Name precise_times( "precise_times" );
 const Name precision( "precision" );
 const Name print_time( "print_time" );

nestkernel/nest_names.h

@@ -333,7 +333,9 @@ extern const Name port_name;
 extern const Name port_width;
 extern const Name ports;
 extern const Name post_synaptic_element;
+extern const Name post_trace;
 extern const Name pre_synaptic_element;
+extern const Name pre_trace;
 extern const Name precise_times;
 extern const Name precision;
 extern const Name print_time;

nestkernel/node.cpp

@@ -229,7 +229,7 @@ Node::send_test_event( Node&, rport, synindex, bool )
  * throws IllegalConnection
  */
 void
-Node::register_stdp_connection( double )
+Node::register_stdp_connection( double, double )
 {
   throw IllegalConnection();
 }

nestkernel/node.h

@@ -478,7 +478,7 @@ class Node
    * @throws IllegalConnection
    *
    */
-  virtual void register_stdp_connection( double );
+  virtual void register_stdp_connection( double, double );
 
   /**
    * Handle incoming spike events.