#include "iaf_psc_alpha_presc.h"
#include "exceptions.h"
#include "network.h"
#include "dict.h"
#include "integerdatum.h"
#include "doubledatum.h"
#include "dictutils.h"
#include "numerics.h"
#include "universal_data_logger_impl.h"
#include <limits>
nest::RecordablesMap<nest::iaf_psc_alpha_presc> nest::iaf_psc_alpha_presc::recordablesMap_;
namespace nest
{
template <>
void RecordablesMap<iaf_psc_alpha_presc>::create()
{
insert_(names::V_m, &iaf_psc_alpha_presc::get_V_m_);
}
}
nest::iaf_psc_alpha_presc::Parameters_::Parameters_()
: tau_m_ ( 10.0 ),
tau_syn_( 2.0 ),
c_m_ (250.0 ),
t_ref_ ( 2.0 ),
E_L_ (-70.0 ),
I_e_ ( 0.0 ),
U_th_ (-55.0-E_L_),
U_min_ (-std::numeric_limits<double_t>::infinity()),
U_reset_(-70.0-E_L_),
Interpol_(iaf_psc_alpha_presc::LINEAR)
{}
nest::iaf_psc_alpha_presc::State_::State_()
: y0_(0.0),
y1_(0.0),
y2_(0.0),
y3_(0.0),
r_(0),
last_spike_step_(-1),
last_spike_offset_(0.0)
{}
void nest::iaf_psc_alpha_presc::Parameters_::get(DictionaryDatum &d) const
{
def<double>(d, names::E_L, E_L_);
def<double>(d, names::I_e, I_e_);
def<double>(d, names::V_th, U_th_+E_L_);
def<double>(d, names::V_min, U_min_+E_L_);
def<double>(d, names::V_reset, U_reset_+E_L_);
def<double>(d, names::C_m, c_m_);
def<double>(d, names::tau_m, tau_m_);
def<double>(d, names::tau_syn, tau_syn_);
def<double>(d, names::t_ref, t_ref_);
def<long>(d, names::Interpol_Order, Interpol_);
}
double nest::iaf_psc_alpha_presc::Parameters_::set(const DictionaryDatum& d)
{
const double ELold = E_L_;
updateValue<double>(d, names::E_L, E_L_);
const double delta_EL = E_L_ - ELold;
updateValue<double>(d, names::tau_m, tau_m_);
updateValue<double>(d, names::tau_syn, tau_syn_);
updateValue<double>(d, names::C_m, c_m_);
updateValue<double>(d, names::t_ref, t_ref_);
updateValue<double>(d, names::I_e, I_e_);
if (updateValue<double>(d, names::V_th, U_th_))
U_th_ -= E_L_;
else
U_th_ -= delta_EL;
if (updateValue<double>(d, names::V_min, U_min_))
U_min_ -= E_L_;
else
U_min_ -= delta_EL;
if (updateValue<double>(d, names::V_reset, U_reset_))
U_reset_ -= E_L_;
else
U_reset_ -= delta_EL;
long_t tmp;
if ( updateValue<long_t>(d, names::Interpol_Order, tmp) )
{
if ( NO_INTERPOL <= tmp && tmp < END_INTERP_ORDER )
Interpol_ = static_cast<interpOrder>(tmp);
else
throw BadProperty("Invalid interpolation order. "
"Valid orders are 0, 1, 2, 3.");
}
if ( U_reset_ >= U_th_ )
throw BadProperty("Reset potential must be smaller than threshold.");
if ( U_reset_ < U_min_ )
throw BadProperty("Reset potential must be greater equal minimum potential.");
if ( c_m_ <= 0 )
throw BadProperty("Capacitance must be strictly positive.");
if ( t_ref_ < 0 )
throw BadProperty("Refractory time must not be negative.");
if ( tau_m_ <= 0 || tau_syn_ <= 0 )
throw BadProperty("All time constants must be strictly positive.");
if ( tau_m_ == tau_syn_ )
throw BadProperty("Membrane and synapse time constant(s) must differ."
"See note in documentation.");
return delta_EL;
}
void nest::iaf_psc_alpha_presc::State_::get(DictionaryDatum &d,
const Parameters_& p) const
{
def<double>(d, names::V_m, y3_ + p.E_L_);
def<double>(d, names::t_spike, Time(Time::step(last_spike_step_)).get_ms());
def<double>(d, names::offset, last_spike_offset_);
}
void nest::iaf_psc_alpha_presc::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_alpha_presc::Buffers_::Buffers_(iaf_psc_alpha_presc& n)
: logger_(n)
{}
nest::iaf_psc_alpha_presc::Buffers_::Buffers_(const Buffers_&, iaf_psc_alpha_presc& n)
: logger_(n)
{}
nest::iaf_psc_alpha_presc::iaf_psc_alpha_presc()
: Node(),
P_(),
S_(),
B_(*this)
{
recordablesMap_.create();
}
nest::iaf_psc_alpha_presc::iaf_psc_alpha_presc(const iaf_psc_alpha_presc& n)
: Node(n),
P_(n.P_),
S_(n.S_),
B_(n.B_, *this)
{}
void nest::iaf_psc_alpha_presc::init_state_(const Node& proto)
{
const iaf_psc_alpha_presc& pr = downcast<iaf_psc_alpha_presc>(proto);
S_ = pr.S_;
}
void nest::iaf_psc_alpha_presc::init_buffers_()
{
B_.spike_y1_.clear();
B_.spike_y2_.clear();
B_.spike_y3_.clear();
B_.currents_.clear();
B_.logger_.reset();
}
void nest::iaf_psc_alpha_presc::calibrate()
{
B_.logger_.init();
V_.h_ms_ = Time::get_resolution().get_ms();
V_.PSCInitialValue_ = 1.0 * numerics::e / P_.tau_syn_;
V_.gamma_ = 1/P_.c_m_ / ( 1/P_.tau_syn_ - 1/P_.tau_m_ );
V_.gamma_sq_ = 1/P_.c_m_ / ( ( 1/P_.tau_syn_ - 1/P_.tau_m_ ) * ( 1/P_.tau_syn_ - 1/P_.tau_m_ ) );
V_.expm1_tau_m_ = numerics::expm1(-V_.h_ms_/P_.tau_m_);
V_.expm1_tau_syn_ = numerics::expm1(-V_.h_ms_/P_.tau_syn_);
V_.P30_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_;
V_.P31_ = V_.gamma_sq_ * V_.expm1_tau_m_ - V_.gamma_sq_ * V_.expm1_tau_syn_
- V_.h_ms_ * V_.gamma_ * V_.expm1_tau_syn_ - V_.h_ms_ * V_.gamma_;
V_.P32_ = V_.gamma_ * V_.expm1_tau_m_ - V_.gamma_ * V_.expm1_tau_syn_;
V_.refractory_steps_ = Time(Time::ms(P_.t_ref_)).get_steps();
assert(V_.refractory_steps_ >= 0);
}
void nest::iaf_psc_alpha_presc::update(Time const & origin,
const long_t from, const long_t to)
{
assert(to >= 0);
assert(static_cast<delay>(from) < Scheduler::get_min_delay());
assert(from < to);
if ( S_.y3_ >= P_.U_th_ )
{
set_spiketime(Time::step(origin.get_steps() + from + 1));
S_.last_spike_offset_ = V_.h_ms_ * (1-std::numeric_limits<double_t>::epsilon());
S_.y3_ = P_.U_reset_;
S_.r_ = V_.refractory_steps_;
SpikeEvent se;
se.set_offset(S_.last_spike_offset_);
network()->send(*this, se, from);
}
for ( long_t lag = from ; lag < to ; ++lag )
{
const long_t T = origin.get_steps() + lag;
V_.y0_before_ = S_.y0_;
V_.y1_before_ = S_.y1_;
V_.y2_before_ = S_.y2_;
V_.y3_before_ = S_.y3_;
const double_t dy3 = B_.spike_y3_.get_value(lag);
if ( S_.r_ == 0 )
{
S_.y3_ = V_.P30_ * (P_.I_e_+S_.y0_) + V_.P31_*S_.y1_ + V_.P32_*S_.y2_ + V_.expm1_tau_m_*S_.y3_ + S_.y3_;
S_.y3_ += dy3;
S_.y3_ = ( S_.y3_< P_.U_min_ ? P_.U_min_ : S_.y3_);
}
else if ( S_.r_ == 1 )
{
S_.r_ = 0;
S_.y3_ = P_.U_reset_ +
update_y3_delta_() + dy3 - dy3 * (1 - S_.last_spike_offset_/V_.h_ms_);
S_.y3_ = ( S_.y3_< P_.U_min_ ? P_.U_min_ : S_.y3_);
}
else
{
--S_.r_;
}
S_.y2_ =V_. expm1_tau_syn_*V_.h_ms_*S_.y1_ + V_.expm1_tau_syn_*S_.y2_ + V_.h_ms_*S_.y1_ + S_.y2_;
S_.y1_ = V_.expm1_tau_syn_*S_.y1_ + S_.y1_;
S_.y1_ += B_.spike_y1_.get_value(lag);
S_.y2_ += B_.spike_y2_.get_value(lag);
if ( S_.y3_ >= P_.U_th_ )
{
set_spiketime(Time::step(T+1));
S_.last_spike_offset_ = V_.h_ms_ - thresh_find_(V_.h_ms_);
S_.y3_ = P_.U_reset_;
S_.r_ = V_.refractory_steps_;
SpikeEvent se;
se.set_offset(S_.last_spike_offset_);
network()->send(*this, se, lag);
}
S_.y0_ = B_.currents_.get_value(lag);
B_.logger_.record_data(origin.get_steps()+lag);
}
}
void nest::iaf_psc_alpha_presc::handle(SpikeEvent & e)
{
assert(e.get_delay() > 0 );
const long_t Tdeliver = e.get_rel_delivery_steps(network()->get_slice_origin());
const double_t spike_weight = V_.PSCInitialValue_ * e.get_weight() * e.get_multiplicity();
const double_t dt = e.get_offset();
const double_t ps_e_TauSyn = numerics::expm1(-dt/P_.tau_syn_);
const double_t ps_e_Tau = numerics::expm1(-dt/P_.tau_m_);
const double_t ps_P31 = V_.gamma_sq_ * ps_e_Tau - V_.gamma_sq_ * ps_e_TauSyn
- dt * V_.gamma_ * ps_e_TauSyn - dt * V_.gamma_;
B_.spike_y1_.add_value(Tdeliver, spike_weight*ps_e_TauSyn + spike_weight);
B_.spike_y2_.add_value(Tdeliver, spike_weight*dt*ps_e_TauSyn + spike_weight*dt);
B_.spike_y3_.add_value(Tdeliver, spike_weight*ps_P31);
}
void nest::iaf_psc_alpha_presc::handle(CurrentEvent& e)
{
assert(e.get_delay() > 0);
const double_t c=e.get_current();
const double_t w=e.get_weight();
B_.currents_.add_value(e.get_rel_delivery_steps(network()->get_slice_origin()),
w * c);
}
void nest::iaf_psc_alpha_presc::handle(DataLoggingRequest& e)
{
B_.logger_.handle(e);
}
inline
void nest::iaf_psc_alpha_presc::set_spiketime(Time const & now)
{
S_.last_spike_step_ = now.get_steps();
}
inline
nest::Time nest::iaf_psc_alpha_presc::get_spiketime() const
{
return Time::step(S_.last_spike_step_);
}
nest::double_t nest::iaf_psc_alpha_presc::update_y3_delta_() const
{
const double t_th = V_.h_ms_ - S_.last_spike_offset_;
double_t ps_e_TauSyn = numerics::expm1(-t_th/P_.tau_syn_);
const double ps_y2 = t_th * ps_e_TauSyn * V_.y1_before_
+ ps_e_TauSyn * V_.y2_before_ + t_th * V_.y1_before_ + V_.y2_before_ ;
const double ps_y1 = ps_e_TauSyn * V_.y1_before_ + V_.y1_before_ ;
const double_t dt = S_.last_spike_offset_;
ps_e_TauSyn = numerics::expm1(-dt / P_.tau_syn_);
const double_t ps_e_Tau = numerics::expm1(-dt/P_.tau_m_);
const double_t ps_P30 = - P_.tau_m_ / P_.c_m_ * ps_e_Tau;
const double_t ps_P31 = V_.gamma_sq_ * ps_e_Tau - V_.gamma_sq_ * ps_e_TauSyn
- dt*V_.gamma_*ps_e_TauSyn - dt*V_.gamma_;
const double_t ps_P32 = V_.gamma_*ps_e_Tau - V_.gamma_* ps_e_TauSyn;
return ps_P30 * (P_.I_e_+V_.y0_before_) + ps_P31 * ps_y1 + ps_P32 * ps_y2;
}
inline
nest::double_t nest::iaf_psc_alpha_presc::thresh_find_(double_t const dt) const
{
switch (P_.Interpol_) {
case NO_INTERPOL: return dt;
case LINEAR : return thresh_find1_(dt);
case QUADRATIC : return thresh_find2_(dt);
case CUBIC : return thresh_find3_(dt);
default:
network()->message(SLIInterpreter::M_ERROR, "iaf_psc_alpha_presc::thresh_find_()",
"Invalid interpolation---Internal model error.");
throw BadProperty();
}
return 0;
}
nest::double_t nest::iaf_psc_alpha_presc::thresh_find1_(double_t const dt) const
{
double_t tau = ( P_.U_th_ - V_.y3_before_ ) * dt / ( S_.y3_ - V_.y3_before_ );
return tau;
}
nest::double_t nest::iaf_psc_alpha_presc::thresh_find2_(double_t const dt) const
{
const double_t h_sq = dt * dt;
const double_t derivative = - V_.y3_before_/P_.tau_m_ + (P_.I_e_ + V_.y0_before_ + V_.y2_before_)/P_.c_m_;
const double_t a = (-V_.y3_before_/h_sq) + (S_.y3_/h_sq) - (derivative/dt);
const double_t b = derivative;
const double_t c = V_.y3_before_;
const double_t sqr_ = std::sqrt(b*b - 4*a*c + 4*a*P_.U_th_);
const double_t tau1 = (-b + sqr_) / (2*a);
const double_t tau2 = (- b - sqr_) / (2*a);
if (tau1 >= 0)
return tau1;
else if (tau2 >= 0)
return tau2;
else
return thresh_find1_(dt);
}
nest::double_t nest::iaf_psc_alpha_presc::thresh_find3_(double_t const dt) const
{
const double_t h_ms_=dt;
const double_t h_sq = h_ms_*h_ms_;
const double_t h_cb = h_sq*h_ms_;
const double_t deriv_t1 = - V_.y3_before_/P_.tau_m_ + (P_.I_e_ + V_.y0_before_ + V_.y2_before_)/P_.c_m_;
const double_t deriv_t2 = - S_.y3_/P_.tau_m_ + (P_.I_e_ + S_.y0_ + S_.y2_)/P_.c_m_;
const double_t w3_ = (2 * V_.y3_before_ / h_cb) - (2 * S_.y3_ / h_cb)
+ ( deriv_t1 / h_sq) + ( deriv_t2 / h_sq) ;
const double_t w2_ = - (3 * V_.y3_before_ / h_sq) + (3 * S_.y3_ / h_sq)
- ( 2 * deriv_t1 / h_ms_) - ( deriv_t2 / h_ms_) ;
const double_t w1_ = deriv_t1;
const double_t w0_ = V_.y3_before_;
const double_t r = w2_ / w3_;
const double_t s = w1_ / w3_;
const double_t t = (w0_ - P_.U_th_) / w3_;
const double_t r_sq= r*r;
const double_t p = - r_sq / 3 + s;
const double_t q = 2 * ( r_sq * r ) / 27 - r * s / 3 + t;
const double_t D = std::pow( (p/3), 3) + std::pow( (q/2), 2);
double_t tau1;
double_t tau2;
double_t tau3;
if(D<0){
const double_t roh = std::sqrt( -(p*p*p)/ 27 );
const double_t phi = std::acos( -q/ (2*roh) );
const double_t a = 2 * std::pow(roh, (1.0/3.0));
tau1 = (a * std::cos( phi/3 )) - r/3;
tau2 = (a * std::cos( phi/3 + 2* numerics::pi/3 )) - r/3;
tau3 = (a * std::cos( phi/3 + 4* numerics::pi/3 )) - r/3;
}
else{
const double_t sgnq = (q >= 0 ? 1 : -1);
const double_t u = -sgnq * std::pow(std::fabs(q)/2.0 + std::sqrt(D), 1.0/3.0);
const double_t v = - p/(3*u);
tau1= (u+v) - r/3;
if (tau1 >= 0) {
return tau1;
}
else {
return thresh_find2_(dt);
}
}
double tau = (tau1 >= 0) ? tau1 : 2*h_ms_;
if ((tau2 >=0) && (tau2 < tau)) tau = tau2;
if ((tau3 >=0) && (tau3 < tau)) tau = tau3;
return (tau <= h_ms_) ? tau : thresh_find2_(dt);
}