1 /*
2 * aeif_psc_exp.cpp
3 *
4 * This file is part of NEST.
5 *
6 * Copyright (C) 2004 The NEST Initiative
7 *
8 * NEST is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * NEST is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with NEST. If not, see <http://www.gnu.org/licenses/>.
20 *
21 */
22
23 #include "aeif_psc_exp.h"
24
25 #ifdef HAVE_GSL
26
27 // C++ includes:
28 #include <cmath>
29 #include <cstdio>
30 #include <iomanip>
31 #include <iostream>
32 #include <limits>
33
34 // Includes from libnestutil:
35 #include "dict_util.h"
36 #include "numerics.h"
37
38 // Includes from nestkernel:
39 #include "exceptions.h"
40 #include "kernel_manager.h"
41 #include "nest_names.h"
42 #include "universal_data_logger_impl.h"
43
44 // Includes from sli:
45 #include "dict.h"
46 #include "dictutils.h"
47 #include "doubledatum.h"
48 #include "integerdatum.h"
49
50 /* ----------------------------------------------------------------
51 * Recordables map
52 * ---------------------------------------------------------------- */
53
54 nest::RecordablesMap< nest::aeif_psc_exp > nest::aeif_psc_exp::recordablesMap_;
55
56 namespace nest
57 {
58 /*
59 * template specialization must be placed in namespace
60 *
61 * Override the create() method with one call to RecordablesMap::insert_()
62 * for each quantity to be recorded.
63 */
64 template <>
65 void
create()66 RecordablesMap< aeif_psc_exp >::create()
67 {
68 // use standard names whereever you can for consistency!
69 insert_( names::V_m, &aeif_psc_exp::get_y_elem_< aeif_psc_exp::State_::V_M > );
70 insert_( names::I_syn_ex, &aeif_psc_exp::get_y_elem_< aeif_psc_exp::State_::I_EXC > );
71 insert_( names::I_syn_in, &aeif_psc_exp::get_y_elem_< aeif_psc_exp::State_::I_INH > );
72 insert_( names::w, &aeif_psc_exp::get_y_elem_< aeif_psc_exp::State_::W > );
73 }
74 }
75
76
77 extern "C" int
aeif_psc_exp_dynamics(double,const double y[],double f[],void * pnode)78 nest::aeif_psc_exp_dynamics( double, const double y[], double f[], void* pnode )
79 {
80 // a shorthand
81 typedef nest::aeif_psc_exp::State_ S;
82
83 // get access to node so we can almost work as in a member function
84 assert( pnode );
85 const nest::aeif_psc_exp& node = *( reinterpret_cast< nest::aeif_psc_exp* >( pnode ) );
86
87 const bool is_refractory = node.S_.r_ > 0;
88
89 // y[] here is---and must be---the state vector supplied by the integrator,
90 // not the state vector in the node, node.S_.y[].
91
92 // The following code is verbose for the sake of clarity. We assume that a
93 // good compiler will optimize the verbosity away...
94
95 // Clamp membrane potential to V_reset while refractory, otherwise bound
96 // it to V_peak. Do not use V_.V_peak_ here, since that is set to V_th if
97 // Delta_T == 0.
98 const double& V = is_refractory ? node.P_.V_reset_ : std::min( y[ S::V_M ], node.P_.V_peak_ );
99 // shorthand for the other state variables
100 const double& I_syn_ex = y[ S::I_EXC ];
101 const double& I_syn_in = y[ S::I_INH ];
102 const double& w = y[ S::W ];
103
104 const double I_spike =
105 node.P_.Delta_T == 0. ? 0. : ( node.P_.g_L * node.P_.Delta_T * std::exp( ( V - node.P_.V_th ) / node.P_.Delta_T ) );
106
107 // dv/dt
108 f[ S::V_M ] = is_refractory ? 0. : ( -node.P_.g_L * ( V - node.P_.E_L ) + I_spike + I_syn_ex - I_syn_in - w
109 + node.P_.I_e + node.B_.I_stim_ ) / node.P_.C_m;
110
111 f[ S::I_EXC ] = -I_syn_ex / node.P_.tau_syn_ex; // Exc. synaptic current (pA)
112
113 f[ S::I_INH ] = -I_syn_in / node.P_.tau_syn_in; // Inh. synaptic current (pA)
114
115 // Adaptation current w.
116 f[ S::W ] = ( node.P_.a * ( V - node.P_.E_L ) - w ) / node.P_.tau_w;
117
118 return GSL_SUCCESS;
119 }
120
121 /* ----------------------------------------------------------------
122 * Default constructors defining default parameters and state
123 * ---------------------------------------------------------------- */
124
Parameters_()125 nest::aeif_psc_exp::Parameters_::Parameters_()
126 : V_peak_( 0.0 ) // mV
127 , V_reset_( -60.0 ) // mV
128 , t_ref_( 0.0 ) // ms
129 , g_L( 30.0 ) // nS
130 , C_m( 281.0 ) // pF
131 , E_L( -70.6 ) // mV
132 , Delta_T( 2.0 ) // mV
133 , tau_w( 144.0 ) // ms
134 , a( 4.0 ) // nS
135 , b( 80.5 ) // pA
136 , V_th( -50.4 ) // mV
137 , tau_syn_ex( 0.2 ) // ms
138 , tau_syn_in( 2.0 ) // ms
139 , I_e( 0.0 ) // pA
140 , gsl_error_tol( 1e-6 )
141 {
142 }
143
State_(const Parameters_ & p)144 nest::aeif_psc_exp::State_::State_( const Parameters_& p )
145 : r_( 0 )
146 {
147 y_[ 0 ] = p.E_L;
148 for ( size_t i = 1; i < STATE_VEC_SIZE; ++i )
149 {
150 y_[ i ] = 0;
151 }
152 }
153
State_(const State_ & s)154 nest::aeif_psc_exp::State_::State_( const State_& s )
155 : r_( s.r_ )
156 {
157 for ( size_t i = 0; i < STATE_VEC_SIZE; ++i )
158 {
159 y_[ i ] = s.y_[ i ];
160 }
161 }
162
operator =(const State_ & s)163 nest::aeif_psc_exp::State_& nest::aeif_psc_exp::State_::operator=( const State_& s )
164 {
165 r_ = s.r_;
166 for ( size_t i = 0; i < STATE_VEC_SIZE; ++i )
167 {
168 y_[ i ] = s.y_[ i ];
169 }
170 return *this;
171 }
172
173 /* ----------------------------------------------------------------
174 * Paramater and state extractions and manipulation functions
175 * ---------------------------------------------------------------- */
176
177 void
get(DictionaryDatum & d) const178 nest::aeif_psc_exp::Parameters_::get( DictionaryDatum& d ) const
179 {
180 def< double >( d, names::C_m, C_m );
181 def< double >( d, names::V_th, V_th );
182 def< double >( d, names::t_ref, t_ref_ );
183 def< double >( d, names::g_L, g_L );
184 def< double >( d, names::E_L, E_L );
185 def< double >( d, names::V_reset, V_reset_ );
186 def< double >( d, names::tau_syn_ex, tau_syn_ex );
187 def< double >( d, names::tau_syn_in, tau_syn_in );
188 def< double >( d, names::a, a );
189 def< double >( d, names::b, b );
190 def< double >( d, names::Delta_T, Delta_T );
191 def< double >( d, names::tau_w, tau_w );
192 def< double >( d, names::I_e, I_e );
193 def< double >( d, names::V_peak, V_peak_ );
194 def< double >( d, names::gsl_error_tol, gsl_error_tol );
195 }
196
197 void
set(const DictionaryDatum & d,Node * node)198 nest::aeif_psc_exp::Parameters_::set( const DictionaryDatum& d, Node* node )
199 {
200 updateValueParam< double >( d, names::V_th, V_th, node );
201 updateValueParam< double >( d, names::V_peak, V_peak_, node );
202 updateValueParam< double >( d, names::t_ref, t_ref_, node );
203 updateValueParam< double >( d, names::E_L, E_L, node );
204 updateValueParam< double >( d, names::V_reset, V_reset_, node );
205
206 updateValueParam< double >( d, names::C_m, C_m, node );
207 updateValueParam< double >( d, names::g_L, g_L, node );
208
209 updateValueParam< double >( d, names::tau_syn_ex, tau_syn_ex, node );
210 updateValueParam< double >( d, names::tau_syn_in, tau_syn_in, node );
211
212 updateValueParam< double >( d, names::a, a, node );
213 updateValueParam< double >( d, names::b, b, node );
214 updateValueParam< double >( d, names::Delta_T, Delta_T, node );
215 updateValueParam< double >( d, names::tau_w, tau_w, node );
216
217 updateValueParam< double >( d, names::I_e, I_e, node );
218
219 updateValueParam< double >( d, names::gsl_error_tol, gsl_error_tol, node );
220
221 if ( V_reset_ >= V_peak_ )
222 {
223 throw BadProperty( "Ensure that V_reset < V_peak ." );
224 }
225
226 if ( Delta_T < 0. )
227 {
228 throw BadProperty( "Delta_T must be positive." );
229 }
230 else if ( Delta_T > 0. )
231 {
232 // check for possible numerical overflow with the exponential divergence at
233 // spike time, keep a 1e20 margin for the subsequent calculations
234 const double max_exp_arg = std::log( std::numeric_limits< double >::max() / 1e20 );
235 if ( ( V_peak_ - V_th ) / Delta_T >= max_exp_arg )
236 {
237 throw BadProperty(
238 "The current combination of V_peak, V_th and Delta_T"
239 "will lead to numerical overflow at spike time; try"
240 "for instance to increase Delta_T or to reduce V_peak"
241 "to avoid this problem." );
242 }
243 }
244
245 if ( V_peak_ < V_th )
246 {
247 throw BadProperty( "V_peak >= V_th required." );
248 }
249
250 if ( C_m <= 0 )
251 {
252 throw BadProperty( "Ensure that C_m > 0" );
253 }
254
255 if ( t_ref_ < 0 )
256 {
257 throw BadProperty( "Refractory time cannot be negative." );
258 }
259
260 if ( tau_syn_ex <= 0 || tau_syn_in <= 0 || tau_w <= 0 )
261 {
262 throw BadProperty( "All time constants must be strictly positive." );
263 }
264
265 if ( gsl_error_tol <= 0. )
266 {
267 throw BadProperty( "The gsl_error_tol must be strictly positive." );
268 }
269 }
270
271 void
get(DictionaryDatum & d) const272 nest::aeif_psc_exp::State_::get( DictionaryDatum& d ) const
273 {
274 def< double >( d, names::V_m, y_[ V_M ] );
275 def< double >( d, names::I_syn_ex, y_[ I_EXC ] );
276 def< double >( d, names::I_syn_in, y_[ I_INH ] );
277 def< double >( d, names::w, y_[ W ] );
278 }
279
280 void
set(const DictionaryDatum & d,const Parameters_ &,Node * node)281 nest::aeif_psc_exp::State_::set( const DictionaryDatum& d, const Parameters_&, Node* node )
282 {
283 updateValueParam< double >( d, names::V_m, y_[ V_M ], node );
284 updateValueParam< double >( d, names::I_syn_ex, y_[ I_EXC ], node );
285 updateValueParam< double >( d, names::I_syn_in, y_[ I_INH ], node );
286 updateValueParam< double >( d, names::w, y_[ W ], node );
287 if ( y_[ I_EXC ] < 0 || y_[ I_INH ] < 0 )
288 {
289 throw BadProperty( "Conductances must not be negative." );
290 }
291 }
292
Buffers_(aeif_psc_exp & n)293 nest::aeif_psc_exp::Buffers_::Buffers_( aeif_psc_exp& n )
294 : logger_( n )
295 , s_( 0 )
296 , c_( 0 )
297 , e_( 0 )
298 {
299 // Initialization of the remaining members is deferred to
300 // init_buffers_().
301 }
302
Buffers_(const Buffers_ &,aeif_psc_exp & n)303 nest::aeif_psc_exp::Buffers_::Buffers_( const Buffers_&, aeif_psc_exp& n )
304 : logger_( n )
305 , s_( 0 )
306 , c_( 0 )
307 , e_( 0 )
308 {
309 // Initialization of the remaining members is deferred to
310 // init_buffers_().
311 }
312
313 /* ----------------------------------------------------------------
314 * Default and copy constructor for node, and destructor
315 * ---------------------------------------------------------------- */
316
aeif_psc_exp()317 nest::aeif_psc_exp::aeif_psc_exp()
318 : ArchivingNode()
319 , P_()
320 , S_( P_ )
321 , B_( *this )
322 {
323 recordablesMap_.create();
324 }
325
aeif_psc_exp(const aeif_psc_exp & n)326 nest::aeif_psc_exp::aeif_psc_exp( const aeif_psc_exp& n )
327 : ArchivingNode( n )
328 , P_( n.P_ )
329 , S_( n.S_ )
330 , B_( n.B_, *this )
331 {
332 }
333
~aeif_psc_exp()334 nest::aeif_psc_exp::~aeif_psc_exp()
335 {
336 // GSL structs may not have been allocated, so we need to protect destruction
337 if ( B_.s_ )
338 {
339 gsl_odeiv_step_free( B_.s_ );
340 }
341 if ( B_.c_ )
342 {
343 gsl_odeiv_control_free( B_.c_ );
344 }
345 if ( B_.e_ )
346 {
347 gsl_odeiv_evolve_free( B_.e_ );
348 }
349 }
350
351 /* ----------------------------------------------------------------
352 * Node initialization functions
353 * ---------------------------------------------------------------- */
354
355 void
init_buffers_()356 nest::aeif_psc_exp::init_buffers_()
357 {
358 B_.spike_exc_.clear(); // includes resize
359 B_.spike_inh_.clear(); // includes resize
360 B_.currents_.clear(); // includes resize
361 ArchivingNode::clear_history();
362
363 B_.logger_.reset();
364
365 B_.step_ = Time::get_resolution().get_ms();
366
367 // We must integrate this model with high-precision to obtain decent results
368 B_.IntegrationStep_ = std::min( 0.01, B_.step_ );
369
370 if ( B_.s_ == 0 )
371 {
372 B_.s_ = gsl_odeiv_step_alloc( gsl_odeiv_step_rkf45, State_::STATE_VEC_SIZE );
373 }
374 else
375 {
376 gsl_odeiv_step_reset( B_.s_ );
377 }
378
379 if ( B_.c_ == 0 )
380 {
381 B_.c_ = gsl_odeiv_control_yp_new( P_.gsl_error_tol, P_.gsl_error_tol );
382 }
383 else
384 {
385 gsl_odeiv_control_init( B_.c_, P_.gsl_error_tol, P_.gsl_error_tol, 0.0, 1.0 );
386 }
387
388 if ( B_.e_ == 0 )
389 {
390 B_.e_ = gsl_odeiv_evolve_alloc( State_::STATE_VEC_SIZE );
391 }
392 else
393 {
394 gsl_odeiv_evolve_reset( B_.e_ );
395 }
396
397 B_.sys_.jacobian = NULL;
398 B_.sys_.dimension = State_::STATE_VEC_SIZE;
399 B_.sys_.params = reinterpret_cast< void* >( this );
400 B_.sys_.function = aeif_psc_exp_dynamics;
401
402 B_.I_stim_ = 0.0;
403 }
404
405 void
calibrate()406 nest::aeif_psc_exp::calibrate()
407 {
408 // ensures initialization in case mm connected after Simulate
409 B_.logger_.init();
410
411 // set the right threshold and GSL function depending on Delta_T
412 if ( P_.Delta_T > 0. )
413 {
414 V_.V_peak = P_.V_peak_;
415 }
416 else
417 {
418 V_.V_peak = P_.V_th; // same as IAF dynamics for spikes if Delta_T == 0.
419 }
420
421 V_.refractory_counts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps();
422 }
423
424 /* ----------------------------------------------------------------
425 * Update and spike handling functions
426 * ---------------------------------------------------------------- */
427
428 void
update(const Time & origin,const long from,const long to)429 nest::aeif_psc_exp::update( const Time& origin, const long from, const long to )
430 {
431 assert( to >= 0 && ( delay ) from < kernel().connection_manager.get_min_delay() );
432 assert( from < to );
433 assert( State_::V_M == 0 );
434
435 for ( long lag = from; lag < to; ++lag )
436 {
437 double t = 0.0;
438
439 // numerical integration with adaptive step size control:
440 // ------------------------------------------------------
441 // gsl_odeiv_evolve_apply performs only a single numerical
442 // integration step, starting from t and bounded by step;
443 // the while-loop ensures integration over the whole simulation
444 // step (0, step] if more than one integration step is needed due
445 // to a small integration step size;
446 // note that (t+IntegrationStep > step) leads to integration over
447 // (t, step] and afterwards setting t to step, but it does not
448 // enforce setting IntegrationStep to step-t
449 while ( t < B_.step_ )
450 {
451 const int status = gsl_odeiv_evolve_apply( B_.e_,
452 B_.c_,
453 B_.s_,
454 &B_.sys_, // system of ODE
455 &t, // from t
456 B_.step_, // to t <= step
457 &B_.IntegrationStep_, // integration step size
458 S_.y_ ); // neuronal state
459 if ( status != GSL_SUCCESS )
460 {
461 throw GSLSolverFailure( get_name(), status );
462 }
463
464 // check for unreasonable values; we allow V_M to explode
465 if ( S_.y_[ State_::V_M ] < -1e3 || S_.y_[ State_::W ] < -1e6 || S_.y_[ State_::W ] > 1e6 )
466 {
467 throw NumericalInstability( get_name() );
468 }
469
470 // spikes are handled inside the while-loop
471 // due to spike-driven adaptation
472 if ( S_.r_ > 0 )
473 {
474 S_.y_[ State_::V_M ] = P_.V_reset_;
475 }
476 else if ( S_.y_[ State_::V_M ] >= V_.V_peak )
477 {
478 S_.y_[ State_::V_M ] = P_.V_reset_;
479 S_.y_[ State_::W ] += P_.b; // spike-driven adaptation
480
481 /* Initialize refractory step counter.
482 * - We need to add 1 to compensate for count-down immediately after
483 * while loop.
484 * - If neuron has no refractory time, set to 0 to avoid refractory
485 * artifact inside while loop.
486 */
487 S_.r_ = V_.refractory_counts_ > 0 ? V_.refractory_counts_ + 1 : 0;
488
489 set_spiketime( Time::step( origin.get_steps() + lag + 1 ) );
490 SpikeEvent se;
491 kernel().event_delivery_manager.send( *this, se, lag );
492 }
493 }
494
495 // decrement refractory count
496 if ( S_.r_ > 0 )
497 {
498 --S_.r_;
499 }
500
501 S_.y_[ State_::I_EXC ] += B_.spike_exc_.get_value( lag );
502 S_.y_[ State_::I_INH ] += B_.spike_inh_.get_value( lag );
503
504 // set new input current
505 B_.I_stim_ = B_.currents_.get_value( lag );
506
507 // log state data
508 B_.logger_.record_data( origin.get_steps() + lag );
509 }
510 }
511
512 void
handle(SpikeEvent & e)513 nest::aeif_psc_exp::handle( SpikeEvent& e )
514 {
515 assert( e.get_delay_steps() > 0 );
516
517 if ( e.get_weight() > 0.0 )
518 {
519 B_.spike_exc_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
520 e.get_weight() * e.get_multiplicity() );
521 }
522 else
523 {
524 B_.spike_inh_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
525 -e.get_weight() * e.get_multiplicity() );
526 }
527 }
528
529 void
handle(CurrentEvent & e)530 nest::aeif_psc_exp::handle( CurrentEvent& e )
531 {
532 assert( e.get_delay_steps() > 0 );
533
534 const double c = e.get_current();
535 const double w = e.get_weight();
536
537 B_.currents_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), w * c );
538 }
539
540 void
handle(DataLoggingRequest & e)541 nest::aeif_psc_exp::handle( DataLoggingRequest& e )
542 {
543 B_.logger_.handle( e );
544 }
545
546 #endif // HAVE_GSL
547