Merge pull request #822 from heplesser/fix-test_dc_ie

Fix problem with iaf_psc_delta_canon and test_dcgen_vs_I_e.sli
nest · Oct 16, 2017 · c9d9325 · c9d9325
2 parents dfa5db2 + 190855b
commit c9d9325
Show file tree

Hide file tree

Showing 3 changed files with 92 additions and 64 deletions.
diff --git a/precise/iaf_psc_delta_canon.cpp b/precise/iaf_psc_delta_canon.cpp
@@ -58,7 +58,7 @@ template <>
 void
 RecordablesMap< iaf_psc_delta_canon >::create()
 {
-  // use standard names whereever you can for consistency!
+  // use standard names wherever you can for consistency!
   insert_( names::V_m, &iaf_psc_delta_canon::get_V_m_ );
 }
 
@@ -258,8 +258,7 @@ iaf_psc_delta_canon::calibrate()
 
   V_.exp_t_ = std::exp( -V_.h_ms_ / P_.tau_m_ );
   V_.expm1_t_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ );
-  V_.v_inf_ = P_.I_e_ * P_.tau_m_ / P_.c_m_;
-  V_.I_contrib_ = -V_.v_inf_ * V_.expm1_t_;
+  V_.R_ = P_.tau_m_ / P_.c_m_;
 
   // t_ref_ is the refractory period in ms
   // refractory_steps_ is the duration of the refractory period in whole
@@ -352,9 +351,8 @@ iaf_psc_delta_canon::update( Time const& origin,
 
 
         // contribution of the stepwise constant current
-        double I_contrib_t = -S_.I_ * P_.tau_m_ / P_.c_m_ * V_.expm1_t_;
-
-        S_.U_ = V_.I_contrib_ + I_contrib_t + V_.expm1_t_ * S_.U_ + S_.U_;
+        const double I_ext = -V_.expm1_t_ * V_.R_ * ( S_.I_ + P_.I_e_ );
+        S_.U_ = I_ext + V_.expm1_t_ * S_.U_ + S_.U_;
 
         S_.U_ =
           S_.U_ < P_.U_min_ ? P_.U_min_ : S_.U_; // lower bound on potential
@@ -478,7 +476,7 @@ nest::iaf_psc_delta_canon::propagate_( const double dt )
 
   // see comment on regular update above
   const double expm1_dt = numerics::expm1( -dt / P_.tau_m_ );
-  const double v_inf = V_.v_inf_ + S_.I_ * P_.tau_m_ / P_.c_m_;
+  const double v_inf = V_.R_ * ( S_.I_ + P_.I_e_ );
   S_.U_ = -v_inf * expm1_dt + S_.U_ * expm1_dt + S_.U_;
 
   return;
@@ -491,9 +489,10 @@ nest::iaf_psc_delta_canon::emit_spike_( Time const& origin,
 {
   assert( S_.U_ >= P_.U_th_ ); // ensure we are superthreshold
 
-  // compute time since threhold crossing
-  double v_inf = V_.v_inf_ + S_.I_ * P_.tau_m_ / P_.c_m_;
-  double dt = -P_.tau_m_ * std::log( ( v_inf - S_.U_ ) / ( v_inf - P_.U_th_ ) );
+  // compute time since threshold crossing
+  const double v_inf = V_.R_ * ( S_.I_ + P_.I_e_ );
+  const double dt =
+    -P_.tau_m_ * std::log( ( v_inf - S_.U_ ) / ( v_inf - P_.U_th_ ) );
 
   // set stamp and offset for spike
   S_.last_spike_step_ = origin.get_steps() + lag + 1;

diff --git a/precise/iaf_psc_delta_canon.h b/precise/iaf_psc_delta_canon.h
@@ -358,10 +358,9 @@ class iaf_psc_delta_canon : public Archiving_Node
    */
   struct Variables_
   {
-    double exp_t_;     //!< @$ e^{-t/\tau_m} @$
-    double expm1_t_;   //!< @$ e^{-t/\tau_m} - 1 @$
-    double v_inf_;     //!< @$ \frac{I_e\tau_m}{c_m} @$
-    double I_contrib_; //!< @$ \frac{I_e\tau_m}{c_m} (1-e^{-t/\tau_m})@$
+    double exp_t_;   //!< @$ e^{-t/\tau_m} @$
+    double expm1_t_; //!< @$ e^{-t/\tau_m} - 1 @$
+    double R_;       //!< @$ \frac{\tau_m}{c_m} @$
 
     double h_ms_; //!< duration of time step [ms]
 

diff --git a/testsuite/unittests/test_dcgen_versus_I_e.sli b/testsuite/unittests/test_dcgen_versus_I_e.sli
@@ -47,64 +47,94 @@ M_ERROR setverbosity
 /amp 123.456 def
 /res [] def  % array to collect bool results
 
-modeldict keys{
-
-    /model Set
-
-    % check if both membrane potential can be measured....
-    /record_VM model GetDefaults /recordables known def
-    record_VM {	model GetDefaults /recordables get { /V_m eq {true} if } forall count 1 eq {; /record_VM true def} {/record_VM false def} ifelse } if
-
-    % .... and I_e can be set
-    model GetDefaults /I_e known record_VM and
+modeldict keys
+{
+  /model Set
+
+  % check if both membrane potential can be measured....
+  /record_VM model GetDefaults /recordables known def
+  record_VM 
+  {	
+    model GetDefaults /recordables get 
+    { 
+      /V_m eq { true } if 
+    } 
+    forall 
+    count 1 eq 
     {
-
-	ResetKernel
-
-	% special case for Poisson Neuron
-
-	model /gif_psc_exp eq { model << /lambda_0 0.0 >> SetDefaults} if 
-	model /gif_cond_exp eq { model << /lambda_0 0.0 >> SetDefaults} if 
-	model /gif_psc_exp_multisynapse eq { model << /lambda_0 0.0 >> SetDefaults} if 
-	model /gif_cond_exp_multisynapse eq { model << /lambda_0 0.0 >> SetDefaults} if
-	model /pp_psc_delta eq { model << /c_2 0.0 >> SetDefaults} if 
-
-	model =only (\t: ) =only
-
-	%executive
-
-	model Create /n1 Set
-	model Create /n2 Set
-
-	%take into account synaptic delay
-	/dc_generator <</start 99. >> Create /dc Set
-
-	% todo: make sure that this procedure is 'stopped' if and only if current events are not supported!
-	{dc n1 Connect
-
-	dc << /amplitude amp >> SetStatus
-
-	100 Simulate
-	n2 << /I_e amp >> SetStatus
-
-	300 Simulate
-
-	% compare final membrane potentials
-	n1 GetStatus /V_m get n2 GetStatus /V_m get eq dup 
-	  {(pass) =}
-	  { (failed: unequal voltages) =  executive } 
-	ifelse 
-	res exch append /res Set
-	} stopped {  clear (DC not allowed) = errordict /newerror false put  } if % skip models that don't support current event
+      ; /record_VM true def
     }
-    if   % model GetDefaults /I_e known record_VM and
+    {
+      /record_VM false def
+    } 
+    ifelse 
+  } 
+  if
+
+  % .... and I_e can be set
+  model GetDefaults /I_e known record_VM and
+  {
+    ResetKernel
+
+    % models requiring special parameters
+    model /gif_psc_exp eq { model << /lambda_0 0.0 >> SetDefaults} if 
+  	  model /gif_cond_exp eq { model << /lambda_0 0.0 >> SetDefaults} if 
+  	  model /gif_psc_exp_multisynapse eq { model << /lambda_0 0.0 >> SetDefaults} if 
+  	  model /gif_cond_exp_multisynapse eq { model << /lambda_0 0.0 >> SetDefaults} if
+  	  model /pp_psc_delta eq { model << /c_2 0.0 >> SetDefaults} if 
+
+  	  model =only (\t: ) =only
+
+  	  model Create /n1 Set
+  	  model Create /n2 Set
+
+  	  %take into account synaptic delay
+  	  /dc_generator <</start 99. >> Create /dc Set
+
+  	  % todo: make sure that this procedure is 'stopped' if 
+  	  % and only if current events are not supported!
+  	  {
+  	    dc n1 Connect
+      dc << /amplitude amp >> SetStatus
+
+  	    100 Simulate
+
+  	    n2 << /I_e amp >> SetStatus
+
+      300 Simulate
+
+      % compare final membrane potentials
+      n1 GetStatus /V_m get /v1 Set
+      n2 GetStatus /V_m get /v2 Set
+
+      v1 v2 eq dup 
+  	    {
+  	      (pass) =
+  	    }
+  	    { 
+  	      (failed with V1: ) =only v1 =only 
+  	      (, V2: ) =only v2 =only 
+  	      (, dV: ) =only v1 v2 sub =
+  	    } 
+      ifelse 
+      res exch append /res Set
+    } 
+    stopped 
+    {  
+      clear (DC not allowed) = 
+      errordict /newerror false put  
+    } if % skip models that don't support current event
+  }
+  if   % model GetDefaults /I_e known record_VM and
+
 } forall % modeldict keys
 
 % combine results, one bool left on stack
 res First res Rest { and } Fold
 
+} 
+assert_or_die
 
-} assert_or_die
 endusing