/*
* iaf_psc_delta.cpp
*
* This file is part of NEST.
*
* Copyright (C) 2004 The NEST Initiative
*
* NEST is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* NEST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NEST. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* iaf_psc_delta is a neuron where the potential jumps on each spike arrival. */
#include "iaf_psc_delta.h"
// C++ includes:
#include <limits>
// Includes from libnestutil:
#include "numerics.h"
// Includes from nestkernel:
#include "exceptions.h"
#include "kernel_manager.h"
#include "universal_data_logger_impl.h"
// Includes from sli:
#include "dict.h"
#include "dictutils.h"
#include "doubledatum.h"
#include "integerdatum.h"
namespace nest
{
/* ----------------------------------------------------------------
* Recordables map
* ---------------------------------------------------------------- */
RecordablesMap< iaf_psc_delta > iaf_psc_delta::recordablesMap_;
// Override the create() method with one call to RecordablesMap::insert_()
// for each quantity to be recorded.
template <>
void
RecordablesMap< iaf_psc_delta >::create()
{
// use standard names whereever you can for consistency!
insert_( names::V_m, &iaf_psc_delta::get_V_m_ );
}
/* ----------------------------------------------------------------
* Default constructors defining default parameters and state
* ---------------------------------------------------------------- */
nest::iaf_psc_delta::Parameters_::Parameters_()
: tau_m_( 10.0 ) // ms
, c_m_( 250.0 ) // pF
, t_ref_( 2.0 ) // ms
, E_L_( -70.0 ) // mV
, I_e_( 0.0 ) // pA
, V_th_( -55.0 - E_L_ ) // mV, rel to E_L_
, V_min_( -std::numeric_limits< double_t >::max() ) // relative E_L_-55.0-E_L_
, V_reset_( -70.0 - E_L_ ) // mV, rel to E_L_
, with_refr_input_( false )
{
}
nest::iaf_psc_delta::State_::State_()
: y0_( 0.0 )
, y3_( 0.0 )
, r_( 0 )
, refr_spikes_buffer_( 0.0 )
{
}
/* ----------------------------------------------------------------
* Parameter and state extractions and manipulation functions
* ---------------------------------------------------------------- */
void
nest::iaf_psc_delta::Parameters_::get( DictionaryDatum& d ) const
{
def< double >( d, names::E_L, E_L_ ); // Resting potential
def< double >( d, names::I_e, I_e_ );
def< double >( d, names::V_th, V_th_ + E_L_ ); // threshold value
def< double >( d, names::V_reset, V_reset_ + E_L_ );
def< double >( d, names::V_min, V_min_ + E_L_ );
def< double >( d, names::C_m, c_m_ );
def< double >( d, names::tau_m, tau_m_ );
def< double >( d, names::t_ref, t_ref_ );
def< bool >( d, "refractory_input", with_refr_input_ );
}
double
nest::iaf_psc_delta::Parameters_::set( const DictionaryDatum& d )
{
// if E_L_ is changed, we need to adjust all variables defined relative to
// E_L_
const double ELold = E_L_;
updateValue< double >( d, names::E_L, E_L_ );
const double delta_EL = E_L_ - ELold;
if ( updateValue< double >( d, names::V_reset, V_reset_ ) )
V_reset_ -= E_L_;
else
V_reset_ -= delta_EL;
if ( updateValue< double >( d, names::V_th, V_th_ ) )
V_th_ -= E_L_;
else
V_th_ -= delta_EL;
if ( updateValue< double >( d, names::V_min, V_min_ ) )
V_min_ -= E_L_;
else
V_min_ -= delta_EL;
updateValue< double >( d, names::I_e, I_e_ );
updateValue< double >( d, names::C_m, c_m_ );
updateValue< double >( d, names::tau_m, tau_m_ );
updateValue< double >( d, names::t_ref, t_ref_ );
if ( V_reset_ >= V_th_ )
throw BadProperty( "Reset potential must be smaller than threshold." );
if ( c_m_ <= 0 )
throw BadProperty( "Capacitance must be >0." );
if ( t_ref_ < 0 )
throw BadProperty( "Refractory time must not be negative." );
if ( tau_m_ <= 0 )
throw BadProperty( "Membrane time constant must be > 0." );
updateValue< bool >( d, "refractory_input", with_refr_input_ );
return delta_EL;
}
void
nest::iaf_psc_delta::State_::get( DictionaryDatum& d,
const Parameters_& p ) const
{
def< double >( d, names::V_m, y3_ + p.E_L_ ); // Membrane potential
}
void
nest::iaf_psc_delta::State_::set( const DictionaryDatum& d,
const Parameters_& p,
double delta_EL )
{
if ( updateValue< double >( d, names::V_m, y3_ ) )
y3_ -= p.E_L_;
else
y3_ -= delta_EL;
}
nest::iaf_psc_delta::Buffers_::Buffers_( iaf_psc_delta& n )
: logger_( n )
{
}
nest::iaf_psc_delta::Buffers_::Buffers_( const Buffers_&, iaf_psc_delta& n )
: logger_( n )
{
}
/* ----------------------------------------------------------------
* Default and copy constructor for node
* ---------------------------------------------------------------- */
nest::iaf_psc_delta::iaf_psc_delta()
: Archiving_Node()
, P_()
, S_()
, B_( *this )
{
recordablesMap_.create();
}
nest::iaf_psc_delta::iaf_psc_delta( const iaf_psc_delta& n )
: Archiving_Node( n )
, P_( n.P_ )
, S_( n.S_ )
, B_( n.B_, *this )
{
}
/* ----------------------------------------------------------------
* Node initialization functions
* ---------------------------------------------------------------- */
void
nest::iaf_psc_delta::init_state_( const Node& proto )
{
const iaf_psc_delta& pr = downcast< iaf_psc_delta >( proto );
S_ = pr.S_;
}
void
nest::iaf_psc_delta::init_buffers_()
{
B_.spikes_.clear(); // includes resize
B_.currents_.clear(); // includes resize
B_.logger_.reset(); // includes resize
Archiving_Node::clear_history();
}
void
nest::iaf_psc_delta::calibrate()
{
B_.logger_.init();
const double h = Time::get_resolution().get_ms();
V_.P33_ = std::exp( -h / P_.tau_m_ );
V_.P30_ = 1 / P_.c_m_ * ( 1 - V_.P33_ ) * P_.tau_m_;
// TauR specifies the length of the absolute refractory period as
// a double_t in ms. The grid based iaf_psp_delta can only handle refractory
// periods that are integer multiples of the computation step size (h).
// To ensure consistency with the overall simulation scheme such conversion
// should be carried out via objects of class nest::Time. The conversion
// requires 2 steps:
// 1. A time object r is constructed defining representation of
// TauR in tics. This representation is then converted to computation
// time steps again by a strategy defined by class nest::Time.
// 2. The refractory time in units of steps is read out get_steps(), a
// member function of class nest::Time.
//
// The definition of the refractory period of the iaf_psc_delta is consistent
// the one of iaf_neuron_ps.
//
// Choosing a TauR that is not an integer multiple of the computation time
// step h will leed to accurate (up to the resolution h) and self-consistent
// results. However, a neuron model capable of operating with real valued
// spike time may exhibit a different effective refractory time.
V_.RefractoryCounts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps();
// since t_ref_ >= 0, this can only fail in error
assert( V_.RefractoryCounts_ >= 0 );
}
/* ----------------------------------------------------------------
* Update and spike handling functions
*/
void
nest::iaf_psc_delta::update( Time const& origin,
const long_t from,
const long_t to )
{
assert(
to >= 0 && ( delay ) from < kernel().connection_manager.get_min_delay() );
assert( from < to );
const double_t h = Time::get_resolution().get_ms();
for ( long_t lag = from; lag < to; ++lag )
{
if ( S_.r_ == 0 )
{
// neuron not refractory
S_.y3_ = V_.P30_ * ( S_.y0_ + P_.I_e_ ) + V_.P33_ * S_.y3_
+ B_.spikes_.get_value( lag );
// if we have accumulated spikes from refractory period,
// add and reset accumulator
if ( P_.with_refr_input_ && S_.refr_spikes_buffer_ != 0.0 )
{
S_.y3_ += S_.refr_spikes_buffer_;
S_.refr_spikes_buffer_ = 0.0;
}
// lower bound of membrane potential
S_.y3_ = ( S_.y3_ < P_.V_min_ ? P_.V_min_ : S_.y3_ );
}
else // neuron is absolute refractory
{
// read spikes from buffer and accumulate them, discounting
// for decay until end of refractory period
if ( P_.with_refr_input_ )
S_.refr_spikes_buffer_ +=
B_.spikes_.get_value( lag ) * std::exp( -S_.r_ * h / P_.tau_m_ );
else
B_.spikes_.get_value( lag ); // clear buffer entry, ignore spike
--S_.r_;
}
// threshold crossing
if ( S_.y3_ >= P_.V_th_ )
{
S_.r_ = V_.RefractoryCounts_;
S_.y3_ = P_.V_reset_;
// EX: must compute spike time
set_spiketime( Time::step( origin.get_steps() + lag + 1 ) );
SpikeEvent se;
kernel().event_delivery_manager.send( *this, se, lag );
}
// set new input current
S_.y0_ = B_.currents_.get_value( lag );
// voltage logging
B_.logger_.record_data( origin.get_steps() + lag );
}
}
void
nest::iaf_psc_delta::handle( SpikeEvent& e )
{
assert( e.get_delay() > 0 );
// EX: We must compute the arrival time of the incoming spike
// explicity, since it depends on delay and offset within
// the update cycle. The way it is done here works, but
// is clumsy and should be improved.
B_.spikes_.add_value(
e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
e.get_weight() * e.get_multiplicity() );
}
void
nest::iaf_psc_delta::handle( CurrentEvent& e )
{
assert( e.get_delay() > 0 );
const double_t c = e.get_current();
const double_t w = e.get_weight();
// add weighted current; HEP 2002-10-04
B_.currents_.add_value(
e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
w * c );
}
void
nest::iaf_psc_delta::handle( DataLoggingRequest& e )
{
B_.logger_.handle( e );
}
} // namespace