Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixes for the gif multisynapse and quantal_stp model #846

Merged
merged 5 commits into from Mar 9, 2018
Merged

Fixes for the gif multisynapse and quantal_stp model #846

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
38234ae
Apply multiple correction to gif_cond_exp_multisynapse model
drodarie Oct 25, 2017
4db51ff
Fix quantal_stp_connection
drodarie Oct 25, 2017
8a8339a
Integrate reviews on names and units.
drodarie Nov 28, 2017
37eac3d
Integrate reviews on names for aeif models
drodarie Dec 2, 2017
8b0e9e9
Fix lambda_0 units.
drodarie Mar 8, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
14 changes: 7 additions & 7 deletions models/aeif_cond_alpha_multisynapse.cpp
View file
Open in desktop
Expand Up @@ -99,7 +99,7 @@ aeif_cond_alpha_multisynapse_dynamics( double,
double I_syn = 0.0;
for ( size_t i = 0; i < node.P_.n_receptors(); ++i )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
I_syn += y[ S::G + j ] * ( node.P_.E_rev[ i ] - V );
}

Expand All @@ -118,7 +118,7 @@ aeif_cond_alpha_multisynapse_dynamics( double,

for ( size_t i = 0; i < node.P_.n_receptors(); ++i )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
// Synaptic conductance derivative dG/dt
f[ S::DG + j ] = -y[ S::DG + j ] / node.P_.tau_syn[ i ];
f[ S::G + j ] = y[ S::DG + j ] - y[ S::G + j ] / node.P_.tau_syn[ i ];
Expand Down Expand Up @@ -321,13 +321,13 @@ aeif_cond_alpha_multisynapse::State_::get( DictionaryDatum& d ) const

for ( size_t i = 0;
i < ( ( y_.size() - State_::NUMBER_OF_FIXED_STATES_ELEMENTS )
/ State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR );
/ State_::NUM_STATE_ELEMENTS_PER_RECEPTOR );
++i )
{
dg->push_back( y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] );
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] );
g->push_back( y_[ State_::G
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] );
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] );
}

( *d )[ names::dg ] = DoubleVectorDatum( dg );
Expand Down Expand Up @@ -481,7 +481,7 @@ aeif_cond_alpha_multisynapse::calibrate()

B_.spikes_.resize( P_.n_receptors() );
S_.y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * P_.n_receptors() ),
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * P_.n_receptors() ),
0.0 );

// reallocate instance of stepping function for ODE GSL solver
Expand Down Expand Up @@ -585,7 +585,7 @@ aeif_cond_alpha_multisynapse::update( Time const& origin,
for ( size_t i = 0; i < P_.n_receptors(); ++i )
{
S_.y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] +=
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] +=
B_.spikes_[ i ].get_value( lag ) * V_.g0_[ i ]; // add incoming spike
}
// set new input current
Expand Down
2 changes: 1 addition & 1 deletion models/aeif_cond_alpha_multisynapse.h
View file
Open in desktop
Expand Up @@ -287,7 +287,7 @@ class aeif_cond_alpha_multisynapse : public Archiving_Node
};

static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 2; // V_M, W
static const size_t NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR = 2; // DG, G
static const size_t NUM_STATE_ELEMENTS_PER_RECEPTOR = 2; // DG, G

std::vector< double > y_; //!< neuron state
int r_; //!< number of refractory steps remaining
Expand Down
14 changes: 7 additions & 7 deletions models/aeif_cond_beta_multisynapse.cpp
View file
Open in desktop
Expand Up @@ -99,7 +99,7 @@ aeif_cond_beta_multisynapse_dynamics( double,
double I_syn = 0.0;
for ( size_t i = 0; i < node.P_.n_receptors(); ++i )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
I_syn += y[ S::G + j ] * ( node.P_.E_rev[ i ] - V );
}

Expand All @@ -118,7 +118,7 @@ aeif_cond_beta_multisynapse_dynamics( double,

for ( size_t i = 0; i < node.P_.n_receptors(); ++i )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
// Synaptic conductance derivative dG/dt
f[ S::DG + j ] = -y[ S::DG + j ] / node.P_.tau_rise[ i ];
f[ S::G + j ] = y[ S::DG + j ] - y[ S::G + j ] / node.P_.tau_decay[ i ];
Expand Down Expand Up @@ -332,13 +332,13 @@ aeif_cond_beta_multisynapse::State_::get( DictionaryDatum& d ) const

for ( size_t i = 0;
i < ( ( y_.size() - State_::NUMBER_OF_FIXED_STATES_ELEMENTS )
/ State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR );
/ State_::NUM_STATE_ELEMENTS_PER_RECEPTOR );
++i )
{
dg->push_back( y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] );
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] );
g->push_back( y_[ State_::G
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] );
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] );
}

( *d )[ names::dg ] = DoubleVectorDatum( dg );
Expand Down Expand Up @@ -520,7 +520,7 @@ aeif_cond_beta_multisynapse::calibrate()

B_.spikes_.resize( P_.n_receptors() );
S_.y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * P_.n_receptors() ),
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * P_.n_receptors() ),
0.0 );

// reallocate instance of stepping function for ODE GSL solver
Expand Down Expand Up @@ -624,7 +624,7 @@ aeif_cond_beta_multisynapse::update( Time const& origin,
for ( size_t i = 0; i < P_.n_receptors(); ++i )
{
S_.y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] +=
+ ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ) ] +=
B_.spikes_[ i ].get_value( lag ) * V_.g0_[ i ]; // add incoming spike
}
// set new input current
Expand Down
2 changes: 1 addition & 1 deletion models/aeif_cond_beta_multisynapse.h
View file
Open in desktop
Expand Up @@ -292,7 +292,7 @@ class aeif_cond_beta_multisynapse : public Archiving_Node
};

static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 2; // V_M, W
static const size_t NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR = 2; // DG, G
static const size_t NUM_STATE_ELEMENTS_PER_RECEPTOR = 2; // DG, G

std::vector< double > y_; //!< neuron state
int r_; //!< number of refractory steps remaining
Expand Down
42 changes: 18 additions & 24 deletions models/gif_cond_exp_multisynapse.cpp
View file
Open in desktop
Expand Up @@ -98,7 +98,7 @@ nest::gif_cond_exp_multisynapse_dynamics( double,
double I_syn = 0.0;
for ( size_t i = 0; i < node.P_.n_receptors(); ++i )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
I_syn += -y[ S::G + j ] * ( V - node.P_.E_rev_[ i ] );
}

Expand All @@ -109,7 +109,7 @@ nest::gif_cond_exp_multisynapse_dynamics( double,
// outputs: dg/dt
for ( size_t i = 0; i < node.P_.n_receptors(); i++ )
{
const size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;
const size_t j = i * S::NUM_STATE_ELEMENTS_PER_RECEPTOR;
f[ S::G + j ] = -y[ S::G + j ] / node.P_.tau_syn_[ i ];
}

Expand All @@ -127,11 +127,11 @@ nest::gif_cond_exp_multisynapse::Parameters_::Parameters_()
, V_reset_( -55.0 ) // mV
, Delta_V_( 0.5 ) // mV
, V_T_star_( -35 ) // mV
, lambda_0_( 0.001 ) // 1/ms
, lambda_0_( 1. ) // 1/s
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This must remain in milliseconds, since the value is used as

const double lambda = P_.lambda_0_ * std::exp( ( S_.V_ - S_.sfa_ ) / P_.Delta_V_ );
...
-numerics::expm1( -lambda * Time::get_resolution().get_ms() )

i.e., it is multiplied by a value in milliseconds.

The user provides/retrieves lambda_0 in 1/s; the value is converted to the internal units of 1/ms in Parameters_::get()/set(). This is the standard approach in NEST.

, t_ref_( 4.0 ) // ms
, c_m_( 80.0 ) // pF
, tau_stc_() // ms
, q_stc_() // nA
, q_stc_() // pA
, tau_sfa_() // ms
, q_sfa_() // mV
, tau_syn_( 1, 2.0 ) // ms
Expand All @@ -143,8 +143,7 @@ nest::gif_cond_exp_multisynapse::Parameters_::Parameters_()
}

nest::gif_cond_exp_multisynapse::State_::State_( const Parameters_& p )
: y_( STATE_VEC_SIZE, 0.0 )
, size_neuron_state_( 0 )
: y_( STATE_VEC_SIZE + NUM_STATE_ELEMENTS_PER_RECEPTOR, 0.0 )
, I_stim_( 0.0 )
, sfa_( 0.0 )
, stc_( 0.0 )
Expand Down Expand Up @@ -174,7 +173,6 @@ nest::gif_cond_exp_multisynapse::State_::State_( const State_& s )
}

y_ = s.y_;
size_neuron_state_ = s.size_neuron_state_;
}

nest::gif_cond_exp_multisynapse::State_&
Expand All @@ -196,7 +194,6 @@ nest::gif_cond_exp_multisynapse::State_&
}

y_ = s.y_;
size_neuron_state_ = s.size_neuron_state_;

I_stim_ = s.I_stim_;
sfa_ = s.sfa_;
Expand Down Expand Up @@ -377,12 +374,11 @@ nest::gif_cond_exp_multisynapse::State_::get( DictionaryDatum& d,
std::vector< double >* g = new std::vector< double >();

for ( size_t i = 0;
i < ( ( y_.size() - State_::NUMBER_OF_FIXED_STATES_ELEMENTS )
/ ( State_::STATE_VEC_SIZE - 1 ) );
i < ( y_.size() - State_::NUMBER_OF_FIXED_STATES_ELEMENTS );
++i )
{
g->push_back(
y_[ State_::G + State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ] );
y_[ State_::G + State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * i ] );
}

( *d )[ names::g ] = DoubleVectorDatum( g );
Expand All @@ -393,6 +389,12 @@ nest::gif_cond_exp_multisynapse::State_::set( const DictionaryDatum& d,
const Parameters_& p )
{
updateValue< double >( d, names::V_m, y_[ V_M ] );
y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ State_::NUM_STATE_ELEMENTS_PER_RECEPTOR * p.n_receptors(),
0.0 );

sfa_elems_.resize( p.tau_sfa_.size(), 0.0 );
stc_elems_.resize( p.tau_stc_.size(), 0.0 );
}

nest::gif_cond_exp_multisynapse::Buffers_::Buffers_(
Expand Down Expand Up @@ -475,7 +477,10 @@ nest::gif_cond_exp_multisynapse::init_state_( const Node& proto )
void
nest::gif_cond_exp_multisynapse::init_buffers_()
{
B_.spikes_.clear(); // includes resize
B_.spikes_.resize( P_.n_receptors() );
for (size_t i =0; i < P_.n_receptors(); ++i)
B_.spikes_[i].clear(); // includes resize

B_.currents_.clear(); //!< includes resize
B_.logger_.reset(); //!< includes resize
Archiving_Node::clear_history();
Expand Down Expand Up @@ -521,15 +526,6 @@ nest::gif_cond_exp_multisynapse::init_buffers_()
void
nest::gif_cond_exp_multisynapse::calibrate()
{

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you explain your changes at this point in more detail? Do I understand correctly that you can remove this here because the state vectors and buffers are now correctly resized in set and init_buffers_?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand what you want me to remove. Is it in the calibrate function ?
In my commit, I fixed the initialization of the state vectors and move them to the correct functions.
I have also fixed the spike buffers initialization which was not adapted to multiple synapses.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to remove anything, I just wanted some more explanation. It is fine for me now.

B_.spikes_.resize( P_.n_receptors() );

S_.y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * P_.n_receptors() ),
0.0 );

S_.size_neuron_state_ = S_.y_.size();

B_.sys_.dimension = S_.y_.size();

B_.logger_.init();
Expand All @@ -549,13 +545,11 @@ nest::gif_cond_exp_multisynapse::calibrate()
{
V_.P_sfa_[ i ] = std::exp( -h / P_.tau_sfa_[ i ] );
}
S_.sfa_elems_.resize( P_.tau_sfa_.size(), 0.0 );

for ( size_t i = 0; i < P_.tau_stc_.size(); i++ )
{
V_.P_stc_[ i ] = std::exp( -h / P_.tau_stc_[ i ] );
}
S_.stc_elems_.resize( P_.tau_stc_.size(), 0.0 );
}

/* ----------------------------------------------------------------
Expand Down Expand Up @@ -624,7 +618,7 @@ nest::gif_cond_exp_multisynapse::update( Time const& origin,

for ( size_t i = 0; i < P_.n_receptors(); i++ )
{
S_.y_[ State_::G + ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR
S_.y_[ State_::G + ( State_::NUM_STATE_ELEMENTS_PER_RECEPTOR
* i ) ] += B_.spikes_[ i ].get_value( lag );
}

Expand Down
9 changes: 4 additions & 5 deletions models/gif_cond_exp_multisynapse.h
View file
Open in desktop
Expand Up @@ -118,14 +118,14 @@

Spike adaptation and firing intensity parameters:
q_stc vector of double - Values added to spike-triggered currents (stc)
after each spike emission in nA.
after each spike emission in pA.
tau_stc vector of double - Time constants of stc variables in ms.
q_sfa vector of double - Values added to spike-frequency adaptation
(sfa) after each spike emission in mV.
tau_sfa vector of double - Time constants of sfa variables in ms.
Delta_V double - Stochasticity level in mV.
lambda_0 double - Stochastic intensity at firing threshold V_T in 1/s.
V_T_star double - Base threshold in mV
V_T_star double - Base threshold in mV.

Synaptic parameters
tau_syn double vector - Time constants of the synaptic conductance
Expand Down Expand Up @@ -240,7 +240,7 @@ class gif_cond_exp_multisynapse : public Archiving_Node
double V_reset_;
double Delta_V_;
double V_T_star_;
double lambda_0_; /** 1/ms */
double lambda_0_; /** 1/s */
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Must remain 1/ms, see above.


/** Refractory period in ms. */
double t_ref_;
Expand Down Expand Up @@ -306,10 +306,9 @@ class gif_cond_exp_multisynapse : public Archiving_Node
};

static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 1; //!< V_M
static const size_t NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR = 1; //!< G
static const size_t NUM_STATE_ELEMENTS_PER_RECEPTOR = 1; //!< G

std::vector< double > y_; //!< neuron state
unsigned int size_neuron_state_; // size of neuron state array

double I_stim_; //!< This is piecewise constant external current
double sfa_; //!< This is the change of the 'threshold' due to adaptation.
Expand Down
6 changes: 2 additions & 4 deletions models/quantal_stp_connection.h
View file
Open in desktop
Expand Up @@ -197,8 +197,6 @@ Quantal_StpConnection< targetidentifierT >::send( Event& e,
double t_lastspike,
const CommonSynapseProperties& )
{
const int vp = get_target( t )->get_vp();

const double h = e.get_stamp().get_ms() - t_lastspike;

// Compute the decay factors, based on the time since the last spike.
Expand All @@ -212,7 +210,7 @@ Quantal_StpConnection< targetidentifierT >::send( Event& e,
// Compute number of sites that recovered during the interval.
for ( int depleted = n_ - a_; depleted > 0; --depleted )
{
if ( kernel().rng_manager.get_rng( vp )->drand() < ( 1.0 - p_decay ) )
if ( kernel().rng_manager.get_rng( t )->drand() < ( 1.0 - p_decay ) )
{
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can now remove line 200 of this file where the virtual process id is retrieved since it is now longer needed.

++a_;
}
Expand All @@ -222,7 +220,7 @@ Quantal_StpConnection< targetidentifierT >::send( Event& e,
int n_release = 0;
for ( int i = a_; i > 0; --i )
{
if ( kernel().rng_manager.get_rng( vp )->drand() < u_ )
if ( kernel().rng_manager.get_rng( t )->drand() < u_ )
{
++n_release;
}
Expand Down