ObjSim/src/layer/iaf_psc_exp.cpp

#include "sys.hpp"   // for libcwd
#include "debug.hpp" // for libcwd
#include "iaf_psc_exp.h"
#include <limits>
#include <cmath>
#include <iostream>

/* ---------------------------------------------------------------- 
 * Default constructors defining default parameters and state
 * ---------------------------------------------------------------- */
    
iaf_psc_exp::Parameters_::Parameters_()
  : Tau_     (  10.0       ),  // in ms
    C_       ( 250.0       ),  // in pF
    t_ref_   (   2.0       ),  // in ms
    U0_      ( -70.0       ),  // in mV
    I_e_     (   0.0       ),  // in pA
    Theta_   ( -55.0 - U0_ ),  // relative U0_
    V_reset_ ( -70.0 - U0_ ),  // in mV
    tau_ex_  (   2.0       ),  // in ms
    tau_in_  (   2.0       )   // in ms
{}

void iaf_psc_exp::Parameters_::set(double Tau_m, double C, double t_ref, double U0,
		  double I_e, double Theta, double V_reset, double tau_ex, double tau_in)
{
	Tau_ = Tau_m;
	C_ = C;
	t_ref_ = t_ref;
	U0_ = U0;
	I_e_ = I_e;
	Theta_ = Theta;
	V_reset_ = V_reset;
	tau_ex_ = tau_ex;
	tau_in_ = tau_in;

}

iaf_psc_exp::State_::State_()
  : i_0_      ( 0.0 ),
    V_m_      ( 0.0 ),
    r_ref_    ( 0   )
{}

iaf_psc_exp::iaf_psc_exp(int n): layer(n), i_syn_ex_(n, 0), i_syn_in_(n, 0), currents(n, 0)
{
	pS_ = new State_[n];
}

iaf_psc_exp::~iaf_psc_exp()
{
	delete[] pS_;
}

float* iaf_psc_exp::GetInputPointer(csimInputChannel InputNumber)
{
    Dout(dc::layer,  "InputNumber=" << int(InputNumber) << "");
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return &i_syn_ex_[0];
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return &i_syn_in_[0];
    } break;
    case csimInputChannel_Membrane: {
        Dout(dc::layer, "Membrane Input");
        return &currents[0];
    } break;
    default:
        return &i_syn_ex_[0];
    }
    return 0;
}


inline double fexp(const double& val)
{
	return std::exp(static_cast<long double>(val));
}


void iaf_psc_exp::calibrate()
{

  const double h = dt;

  // numbering of state vaiables: i_0 = 0, i_syn_ = 1, V_m_ = 2

  // commented out propagators: forward Euler
  // needed to exactly reproduce Tsodyks network
 
  // these P are independent
  V_.P11ex_ = fexp(-h/P_.tau_ex_);
  //P11ex_ = 1.0-h/tau_ex_;

  V_.P11in_ = fexp(-h/P_.tau_in_);
  //P11in_ = 1.0-h/tau_in_;

  V_.P22_ = fexp(-h/P_.Tau_);
  //P22_ = 1.0-h/Tau_;

  // these depend on the above. Please do not change the order.
  // TODO: use expm1 here to improve accuracy for small time steps

  V_.P21ex_ = P_.Tau_/(P_.C_*(1.0-P_.Tau_/P_.tau_ex_)) * V_.P11ex_ * (1.0 - fexp(h*(1.0/P_.tau_ex_-1.0/P_.Tau_)));
  //P21ex_ = h/C_;

  V_.P21in_ = P_.Tau_/(P_.C_*(1.0-P_.Tau_/P_.tau_in_)) * V_.P11in_ * (1.0 - fexp(h*(1.0/P_.tau_in_-1.0/P_.Tau_)));
  //P21in_ = h/C_;

  V_.P20_ = P_.Tau_/P_.C_*(1.0 - V_.P22_);
  //P20_ = h/C_;

  // TauR specifies the length of the absolute refractory period as 
  // a double_t in ms. The grid based iaf_psc_exp 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.
  //
  // 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_ = P_.t_ref_ / dt;

  if ( V_.RefractoryCounts_ < 1 )
    throw BadProperty("Absolute refractory time must be at least one time step.");
}



int iaf_psc_exp::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
	  // evolve from timestep 'from' to timestep 'to' with steps of h each
	  for ( int n=0; n<N; ++n)
	  {
		  State_& S_ = pS_[n];
		  if ( S_.r_ref_ == 0 ) // neuron not refractory, so evolve V
			  S_.V_m_ = S_.V_m_*V_.P22_ + i_syn_ex_[n]*V_.P21ex_ + i_syn_in_[n]*V_.P21in_ + (P_.I_e_+S_.i_0_)*V_.P20_ + currents[n]*dt;
		  else
			  --S_.r_ref_; // neuron is absolute refractory

		  currents[n] = 0; // reset direct membrane input

		  // exponential decaying PSCs
		  i_syn_ex_[n] *= V_.P11ex_;
		  i_syn_in_[n] *= V_.P11in_;

		  if ( S_.V_m_ >= P_.Theta_ )  // threshold crossing
		  {
			  S_.r_ref_ = V_.RefractoryCounts_;
			  S_.V_m_ = P_.V_reset_;

			  firings[N_firings  ][0]=t;
			  firings[N_firings++][1]=n;
			  if (N_firings == N_firings_max) {
				  ThrowTooManySpikes(t);
			  }
		  }
	  }
	  return 0;
}


/*
 *
 *   : Tau_     (  10.0       ),  // in ms
    C_       ( 250.0       ),  // in pF
    t_ref_   (   2.0       ),  // in ms
    U0_      ( -70.0       ),  // in mV
    I_e_     (   0.0       ),  // in pA
    Theta_   ( -55.0 - U0_ ),  // relative U0_
    V_reset_ ( -70.0 - U0_ ),  // in mV
    tau_ex_  (   2.0       ),  // in ms
    tau_in_  (   2.0       )   // in ms
 *
 */
int iaf_psc_exp::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "iaf_psc_exp" << "\"/> \n";
    sstr << "<Tau value=\"" << P_.Tau_<< "\"/> \n";
    sstr << "<C value=\"" << P_.C_ << "\"/> \n";
    sstr << "<t_ref value=\"" << P_.t_ref_ << "\"/> \n";
    sstr << "<U0 value=\"" << P_.U0_<< "\"/> \n";
    sstr << "<I_e value=\"" << P_.I_e_ << "\"/> \n";
    sstr << "<Theta value=\"" << P_.Theta_ << "\"/> \n";
    sstr << "<V_reset value=\"" << P_.V_reset_ << "\"/> \n";
    sstr << "<tau_ex value=\"" << P_.tau_ex_ << "\"/> \n";
    sstr << "<tau_in value=\"" << P_.tau_in_ << "\"/> \n";
    layer::WriteSimInfo(fw, sstr.str());
}

/** Starts binary recording of membrane potential and synaptic potentials
 *
 * @param nobserve number of neurons to observe
 * @param NNumber starting number. Neurons in the range [NNumber,NNumber+nobserve-1] are recorded
 * @return
 */
int iaf_psc_exp::StartBinRec(int nobserve, int NNumber)
{
	// not implemented
	return 0;
}

int iaf_psc_exp::reset(int t)
{
	State_ reset_state;
	for ( int n=0; n<N; ++n)
	{
		pS_[n] = reset_state;
	}
	return 0;
}