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.

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

Correct input handling for aeif_cond_beta_multisynapse #510 and Update aeif_cond_beta_multisynapse to new aeif numerics approach #514 #516

Merged
merged 25 commits into from Oct 26, 2016
Commits
Jump to file or symbol
Failed to load files and symbols.
+161 −126
Diff settings

Always

Just for now

@@ -20,6 +20,8 @@
*
*/
#include <stdio.h>

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

You should not use C-headers.

@heplesser

heplesser Oct 13, 2016

Contributor

You should not use C-headers.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

Ops... sorry I use it for myself for debugging and I forgot it. Removed.

@golosio

golosio Oct 14, 2016

Contributor

Ops... sorry I use it for myself for debugging and I forgot it. Removed.

#include "aeif_cond_beta_multisynapse.h"
#ifdef HAVE_GSL
@@ -41,6 +43,7 @@
#include "doubledatum.h"
#include "integerdatum.h"
using namespace std;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

References to elements of the std namespace should always be explicit, never use using namespace std in NEST.

@heplesser

heplesser Oct 13, 2016

Contributor

References to elements of the std namespace should always be explicit, never use using namespace std in NEST.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

removed

@golosio

golosio Oct 14, 2016

Contributor

removed

/* ----------------------------------------------------------------
* Recordables map
* ---------------------------------------------------------------- */
@@ -76,23 +79,28 @@ aeif_cond_beta_multisynapse::Parameters_::Parameters_()
, t_ref_( 0.0 ) // ms
, g_L( 30.0 ) // nS
, C_m( 281.0 ) // pF
, E_ex( 0.0 ) // mV
, E_in( -85.0 ) // mV
, E_L( -70.6 ) // mV
, Delta_T( 2.0 ) // mV
, tau_w( 144.0 ) // ms
, a( 4.0 ) // nS
, b( 80.5 ) // pA
, V_th( -50.4 ) // mV
, I_e( 0.0 ) // pA
, num_of_receptors( 1 )

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I wonder if we should have this as a member. Since it will always be the same as E_ref.size(), I think we should should just use

inline size_t n_receptors() { return E_rev.size(); }
@heplesser

heplesser Oct 13, 2016

Contributor

I wonder if we should have this as a member. Since it will always be the same as E_ref.size(), I think we should should just use

inline size_t n_receptors() { return E_rev.size(); }

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

replaced with inline function as you suggest.

@golosio

golosio Oct 14, 2016

Contributor

replaced with inline function as you suggest.

, I_e( 0.0 ) // pA
, gsl_error_tol( 1e-6 )
, num_of_receptors_( 1 )
, has_connections_( false )
{
static const double E_rev_def = 0.0; // E_rev default value
static const double taus_rise_def = 2.0; // taus_rise default value
static const double taus_decay_def = 20.0; // taus_decay default value

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Why do you have these constants here? You could initialize the vectors in the initializer list as

taus_rise(1, 2.0)

But we should still wait for the discussion about default values in the next Developer VC.

@heplesser

heplesser Oct 13, 2016

Contributor

Why do you have these constants here? You could initialize the vectors in the initializer list as

taus_rise(1, 2.0)

But we should still wait for the discussion about default values in the next Developer VC.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

good, now the code is much more compact, there is only the initializer list and the body is empty. For now I'm still using the default with one receptor port, but it will be easy to remove it.

@golosio

golosio Oct 14, 2016

Contributor

good, now the code is much more compact, there is only the initializer list and the body is empty. For now I'm still using the default with one receptor port, but it will be easy to remove it.

receptor_types.clear();
E_rev.clear();
taus_rise.clear();
taus_decay.clear();

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

No need to clear here, we are in a constructor!

@heplesser

heplesser Oct 13, 2016

Contributor

No need to clear here, we are in a constructor!

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

removed.

@golosio

golosio Oct 14, 2016

Contributor

removed.

taus_decay.push_back( 20.0 );
taus_rise.push_back( 2.0 );
receptor_types.push_back( 1 );

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I am actually not sure we need or use the receptor_types member anywhere. If it is not used, we should drop it.

@heplesser

heplesser Oct 13, 2016

Contributor

I am actually not sure we need or use the receptor_types member anywhere. If it is not used, we should drop it.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

Right, it was not used anywhere, I removed it.

@golosio

golosio Oct 14, 2016

Contributor

Right, it was not used anywhere, I removed it.

E_rev.push_back( E_rev_def );
taus_rise.push_back( taus_rise_def );
taus_decay.push_back( taus_decay_def );
}
aeif_cond_beta_multisynapse::State_::State_( const Parameters_& p )
@@ -131,10 +139,11 @@ aeif_cond_beta_multisynapse::Parameters_::get( DictionaryDatum& d ) const
def< double >( d, names::g_L, g_L );
def< double >( d, names::E_L, E_L );
def< double >( d, names::V_reset, V_reset_ );
def< double >( d, names::E_ex, E_ex );
def< double >( d, names::E_in, E_in );
def< size_t >( d, names::num_of_receptors, num_of_receptors );
ArrayDatum E_rev_ad( E_rev );
ArrayDatum taus_rise_ad( taus_rise );
ArrayDatum taus_decay_ad( taus_decay );
def< ArrayDatum >( d, names::E_rev, E_rev_ad );
def< ArrayDatum >( d, names::taus_rise, taus_rise_ad );
def< ArrayDatum >( d, names::taus_decay, taus_decay_ad );
def< double >( d, names::a, a );
@@ -144,84 +153,75 @@ aeif_cond_beta_multisynapse::Parameters_::get( DictionaryDatum& d ) const
def< double >( d, names::I_e, I_e );
def< double >( d, names::V_peak, V_peak_ );
def< double >( d, names::gsl_error_tol, gsl_error_tol );
def< int >( d, "n_synapses", num_of_receptors_ );
def< bool >( d, names::has_connections, has_connections_ );
}
void
aeif_cond_beta_multisynapse::Parameters_::set( const DictionaryDatum& d )
{
updateValue< double >( d, names::V_th, V_th );
updateValue< double >( d, names::V_peak, V_peak_ );
updateValue< double >( d, names::t_ref, t_ref_ );
updateValue< double >( d, names::E_L, E_L );
updateValue< double >( d, names::V_reset, V_reset_ );
updateValue< double >( d, names::E_ex, E_ex );
updateValue< double >( d, names::E_in, E_in );
updateValue< double >( d, names::C_m, C_m );
updateValue< double >( d, names::g_L, g_L );
std::vector< double > tau_tmp;
if ( updateValue< std::vector< double > >( d, names::taus_decay, tau_tmp ) )
{
if ( tau_tmp.size() < taus_decay.size() && has_connections_ == true )
std::vector< double > Erev_tmp, taur_tmp, taud_tmp;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

You should avoid defining multiple variables in one statement.

You also don't need temporary variables here, since set() is called on a temporary copy of P_ that will only be copied back to the real P_ if no exception is raised. This may also simplify some of the error checking below.

@heplesser

heplesser Oct 13, 2016

Contributor

You should avoid defining multiple variables in one statement.

You also don't need temporary variables here, since set() is called on a temporary copy of P_ that will only be copied back to the real P_ if no exception is raised. This may also simplify some of the error checking below.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

removed.

@golosio

golosio Oct 14, 2016

Contributor

removed.

bool Erev_flag =
updateValue< std::vector< double > >( d, names::E_rev, Erev_tmp );
bool taur_flag =
updateValue< std::vector< double > >( d, names::taus_rise, taur_tmp );
bool taud_flag =
updateValue< std::vector< double > >( d, names::taus_decay, taud_tmp );
if ( Erev_flag || taur_flag || taud_flag )

This comment has been minimized.

@heplesser

heplesser Oct 15, 2016

Contributor

I think it should be possible to simplify this some more:

// n_receptors() == E_rev.size()
if ( n_receptors() != taus_rise.size() || n_receptors() != taus_decay.size() )
{
  throw BadProperty("E_rev, taus_raus and taus_decays must have the same size.");
}

Now we know all is consistent I would then check n_receptors() != old_n_receptors && has_connections and n_receptors > 0. Then the check that the taus > 0 and tau_rise < tau_decay. That should maybe be protected by if ( taur_flag || taud_flag ) since we have loops here and would want to avoid running through them if the taus have not been changed.

@heplesser

heplesser Oct 15, 2016

Contributor

I think it should be possible to simplify this some more:

// n_receptors() == E_rev.size()
if ( n_receptors() != taus_rise.size() || n_receptors() != taus_decay.size() )
{
  throw BadProperty("E_rev, taus_raus and taus_decays must have the same size.");
}

Now we know all is consistent I would then check n_receptors() != old_n_receptors && has_connections and n_receptors > 0. Then the check that the taus > 0 and tau_rise < tau_decay. That should maybe be protected by if ( taur_flag || taud_flag ) since we have loops here and would want to avoid running through them if the taus have not been changed.

This comment has been minimized.

@golosio

golosio Oct 15, 2016

Contributor

I agree that the code can be made more compact if we simply want to check consistency. However my idea is that we should not only check for consistency, but also give the user a feedback that suggests him how he should modify his command. For instance, if the current number of receptor ports is 4 (or if they have not yet been initialized so the number is 0 or 1) and the user types a command as:
nest.SetStatus(neuron, {'taus_decay':[50.0,20.0,20.0]})
I think a message as:
"If the number of receptor ports is changed, all three arrays E_rev, taus_rise and taus_decay must be provided."
is more helpful to the user than:
"E_rev, taus_rise and taus_decay must have the same size."
because with the latter message it may not be obvious to all users that the problem can be solved by providing all three arrays as arguments.

@golosio

golosio Oct 15, 2016

Contributor

I agree that the code can be made more compact if we simply want to check consistency. However my idea is that we should not only check for consistency, but also give the user a feedback that suggests him how he should modify his command. For instance, if the current number of receptor ports is 4 (or if they have not yet been initialized so the number is 0 or 1) and the user types a command as:
nest.SetStatus(neuron, {'taus_decay':[50.0,20.0,20.0]})
I think a message as:
"If the number of receptor ports is changed, all three arrays E_rev, taus_rise and taus_decay must be provided."
is more helpful to the user than:
"E_rev, taus_rise and taus_decay must have the same size."
because with the latter message it may not be obvious to all users that the problem can be solved by providing all three arrays as arguments.

{ // receptor arrays have been modified
if ( !Erev_flag || !taur_flag || !taud_flag )
{
throw BadProperty(
"The neuron has connections, therefore the number of ports cannot be "
"reduced." );
"The reversal potential, synaptic rise time and synaptic decay time "
"arrays must be modified altogether." );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 9, 2016

Contributor

I don't think "altogether" is the appropriate term here...
Beyond that, I think what we want to check is only change in length... users should be able to modify only one of the properties as long as they change the values inside and not the number of values, right?

@Silmathoron

Silmathoron Oct 9, 2016

Contributor

I don't think "altogether" is the appropriate term here...
Beyond that, I think what we want to check is only change in length... users should be able to modify only one of the properties as long as they change the values inside and not the number of values, right?

This comment has been minimized.

@golosio

golosio Oct 10, 2016

Contributor

It's a problem with my English, I will fix the error message.... About the implementation of the constraints, I followed @heplesser suggestions in comment #510 (comment) when he says:
"if a call to SetStatus contains taus_syn, it also must contain E_rev, and they must be of the same length". In case of aeif_cond_alpha_multisynapse there is no taus_syn, but taus_rise and taus_decay, so I guess that SetStatus should either contain none of them or all three...
If I understand well, what you are proposing is that you can also call SetStatus with one or two of them as far as they have the same length that they had previously. For me it would be fine. @heplesser what do you think about this?

@golosio

golosio Oct 10, 2016

Contributor

It's a problem with my English, I will fix the error message.... About the implementation of the constraints, I followed @heplesser suggestions in comment #510 (comment) when he says:
"if a call to SetStatus contains taus_syn, it also must contain E_rev, and they must be of the same length". In case of aeif_cond_alpha_multisynapse there is no taus_syn, but taus_rise and taus_decay, so I guess that SetStatus should either contain none of them or all three...
If I understand well, what you are proposing is that you can also call SetStatus with one or two of them as far as they have the same length that they had previously. For me it would be fine. @heplesser what do you think about this?

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

@golosio @Silmathoron Apologies for being imprecise in my previous comment. I indeed meant what @Silmathoron suggested: As long as you do not change the number of synapses, you can change the values in one array alone (e.g. taus_syn), but when you change then number, you have to provide all arrays and they need to have consistent size.

@heplesser

heplesser Oct 13, 2016

Contributor

@golosio @Silmathoron Apologies for being imprecise in my previous comment. I indeed meant what @Silmathoron suggested: As long as you do not change the number of synapses, you can change the values in one array alone (e.g. taus_syn), but when you change then number, you have to provide all arrays and they need to have consistent size.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

I changed the error message. Please check if it is ok now, otherwise please suggest a suitable one. I modified the behavior for changes in the arrays according to @Silmathoron indication.

@golosio

golosio Oct 14, 2016

Contributor

I changed the error message. Please check if it is ok now, otherwise please suggest a suitable one. I modified the behavior for changes in the arrays according to @Silmathoron indication.

}
for ( size_t i = 0; i < tau_tmp.size(); ++i )
if ( ( Erev_tmp.size() != taur_tmp.size() )
|| ( Erev_tmp.size() != taud_tmp.size() ) )
{
if ( tau_tmp[ i ] <= 0 )
{
throw BadProperty(
"All synaptic time constants must be strictly positive" );
}
}
taus_decay = tau_tmp;
num_of_receptors_ = taus_decay.size();
tau_tmp.clear();
if ( taus_rise.size() < taus_decay.size() )
{
for ( size_t i = 0; i < taus_decay.size(); ++i )
{
taus_rise.push_back( taus_decay[ i ] / 10. );
// if taus_rise is not defined explicitly or if it has less elements
// than taus_decay, it will be set to taus_decay/10
}
throw BadProperty(
"The reversal potential, synaptic rise time and synaptic decay time "
"arrays must have the same size." );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Are both tests really necessary?
Wouldn't

bool size_changed = ( E_rev.size() != old_n_receptors
  || taus_rise.size() != old_n_receptors
  || taus_decay.size() != old_n_receptors );
bool change_is_consistent = ( E_rev.size() == taus_rise.size() )
  && ( E_rev.size() == taus_decay.size() );

if ( size_changed && !change_is_consistent )

be sufficient?

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Are both tests really necessary?
Wouldn't

bool size_changed = ( E_rev.size() != old_n_receptors
  || taus_rise.size() != old_n_receptors
  || taus_decay.size() != old_n_receptors );
bool change_is_consistent = ( E_rev.size() == taus_rise.size() )
  && ( E_rev.size() == taus_decay.size() );

if ( size_changed && !change_is_consistent )

be sufficient?

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

About this point, please have a look at my answer in
#516 (comment)
I agree that if we simply want to check that the changes are consistent, the code can be made more compact. However I think that we should also think about what error messages would be more helpful for the user.
I think we should distinguish those two erroneous situations:

  1. The user types a command that provides only one or two arrays, and they have a different size from the old one. Then the error message should be something like:
    "If the number of receptor ports is changed, all three arrays E_rev, taus_rise and taus_decay must be provided."
  2. The user types a command that provides three arrays, but they have inconsistent sizes
    Then the error message should be something like:
    "The reversal potential, synaptic rise time and synaptic decay time arrays must have the same size."
    If we first agree on which error messages should be printed in different situations, then we can discuss on what is the best way to implement it in the code.
@golosio

golosio Oct 18, 2016

Contributor

About this point, please have a look at my answer in
#516 (comment)
I agree that if we simply want to check that the changes are consistent, the code can be made more compact. However I think that we should also think about what error messages would be more helpful for the user.
I think we should distinguish those two erroneous situations:

  1. The user types a command that provides only one or two arrays, and they have a different size from the old one. Then the error message should be something like:
    "If the number of receptor ports is changed, all three arrays E_rev, taus_rise and taus_decay must be provided."
  2. The user types a command that provides three arrays, but they have inconsistent sizes
    Then the error message should be something like:
    "The reversal potential, synaptic rise time and synaptic decay time arrays must have the same size."
    If we first agree on which error messages should be printed in different situations, then we can discuss on what is the best way to implement it in the code.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

My bad, I remember reading about something like what you said. Ok, I understand your point.

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

My bad, I remember reading about something like what you said. Ok, I understand your point.

}
}
if ( updateValue< std::vector< double > >( d, names::taus_rise, tau_tmp ) )
{
if ( taus_decay.size() == 0 )
if ( taur_tmp.size() == 0 )
{
throw BadProperty(
"Synaptic decay times must be defined before rise times." );
throw BadProperty( "The neuron must have at least one port." );
}
if ( tau_tmp.size() != taus_decay.size() )
if ( taur_tmp.size() < num_of_receptors && has_connections_ == true )
{
throw BadProperty(
"The number of ports for synaptic rise times must be the same "
"as that of decay times." );
"The neuron has connections, therefore the number of ports cannot be "
"reduced." );
}
for ( size_t i = 0; i < tau_tmp.size(); ++i )
for ( size_t i = 0; i < taur_tmp.size(); ++i )
{
if ( tau_tmp[ i ] <= 0 )
if ( taur_tmp[ i ] <= 0 || taud_tmp[ i ] <= 0 )
{
throw BadProperty(
"All synaptic time constants must be strictly positive" );
}
if ( tau_tmp[ i ] > taus_decay[ i ] )
if ( taud_tmp[ i ] < taur_tmp[ i ] )
{
throw BadProperty(
"Synaptic rise time must be smaller than or equal to decay time." );
}
}
taus_rise = tau_tmp;
E_rev = Erev_tmp;
taus_rise = taur_tmp;
taus_decay = taud_tmp;
num_of_receptors = taur_tmp.size();
receptor_types.resize( num_of_receptors );

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

See above, I think we can drop receptor_types.

@heplesser

heplesser Oct 13, 2016

Contributor

See above, I think we can drop receptor_types.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

for ( size_t i = 0; i < num_of_receptors; i++ )
{
receptor_types[ i ] = i + 1;
}
}
updateValue< double >( d, names::a, a );
@@ -233,14 +233,29 @@ aeif_cond_beta_multisynapse::Parameters_::set( const DictionaryDatum& d )
updateValue< double >( d, names::gsl_error_tol, gsl_error_tol );
if ( V_peak_ <= V_th )
if ( V_reset_ >= V_peak_ )

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I think the if ( V_peak_ <= V_th ) test just disappeared... and it shouldn't have ^^

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I think the if ( V_peak_ <= V_th ) test just disappeared... and it shouldn't have ^^

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

I copied this part from aeif_cond_alpha, where this check was also missing... should we put it back in both this model and in aeif_cond_alpha?

@golosio

golosio Oct 18, 2016

Contributor

I copied this part from aeif_cond_alpha, where this check was also missing... should we put it back in both this model and in aeif_cond_alpha?

This comment has been minimized.

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Ok, then this is probably my fault ^^"
No, just change it here, I'll correct the aeif_cond_alpha it in the PR suppressing the dynamics_DT0

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Ok, then this is probably my fault ^^"
No, just change it here, I'll correct the aeif_cond_alpha it in the PR suppressing the dynamics_DT0

{
throw BadProperty( "V_peak must be larger than threshold." );
throw BadProperty( "Ensure that: V_reset < V_peak ." );
}
if ( V_reset_ >= V_peak_ )
if ( Delta_T < 0. )
{
throw BadProperty( "Ensure that: V_reset < V_peak ." );
throw BadProperty( "Delta_T must be positive." );
}
else if ( Delta_T > 0. )
{
// check for possible numerical overflow with the exponential divergence at
// spike time, keep a 1e20 margin for the subsequent calculations
const double max_exp_arg =
std::log( std::numeric_limits< double >::max() / 1e20 );
if ( ( V_peak_ - V_th ) / Delta_T >= max_exp_arg )
{
throw BadProperty(
"The current combination of V_peak, V_th and Delta_T"
"will lead to numerical overflow at spike time; try"
"for instance to increase Delta_T or to reduce V_peak"
"to avoid this problem." );
}
}
if ( C_m <= 0 )
@@ -425,6 +440,15 @@ aeif_cond_beta_multisynapse::init_buffers_()
B_.sys_.params = reinterpret_cast< void* >( this );
// B_.sys_.dimension is assigned in calibrate()
B_.I_stim_ = 0.0;
// set the right threshold depending on Delta_T
if ( P_.Delta_T > 0. )
{
B_.V_peak = P_.V_peak_;

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I'm almost sure this is supposed to be in calibrate since we want this check to be done every time Simulate is called, so V_peak should be in V_; @heplesser, could you confirm?

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I'm almost sure this is supposed to be in calibrate since we want this check to be done every time Simulate is called, so V_peak should be in V_; @heplesser, could you confirm?

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

done

@golosio

golosio Oct 18, 2016

Contributor

done

}
else
{
B_.V_peak = P_.V_th; // same as IAF dynamics for spikes if Delta_T == 0.
}
}
void
@@ -433,15 +457,15 @@ aeif_cond_beta_multisynapse::calibrate()
// ensures initialization in case mm connected after Simulate
B_.logger_.init();
P_.receptor_types_.resize( P_.num_of_receptors_ );
for ( size_t i = 0; i < P_.num_of_receptors_; i++ )
P_.receptor_types.resize( P_.num_of_receptors );
for ( size_t i = 0; i < P_.num_of_receptors; i++ )
{
P_.receptor_types_[ i ] = i + 1;
P_.receptor_types[ i ] = i + 1;
}
V_.g0_.resize( P_.num_of_receptors_ );
V_.g0_.resize( P_.num_of_receptors );
for ( size_t i = 0; i < P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < P_.num_of_receptors; ++i )
{
// the denominator (denom1) that appears in the expression of the peak time
// is computed here to check that it is != 0
@@ -477,10 +501,10 @@ aeif_cond_beta_multisynapse::calibrate()
assert( V_.RefractoryCounts_
>= 0 ); // since t_ref_ >= 0, this can only fail in error
B_.spikes_.resize( P_.num_of_receptors_ );
B_.spikes_.resize( P_.num_of_receptors );
S_.y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR
* P_.num_of_receptors_ ) );
* P_.num_of_receptors ) );
// reallocate instance of stepping function for ODE GSL solver
if ( B_.s_ != 0 )
@@ -563,7 +587,7 @@ aeif_cond_beta_multisynapse::update( Time const& origin,
{
S_.y_[ State_::V_M ] = P_.V_reset_;
}
else if ( S_.y_[ State_::V_M ] >= P_.V_peak_ )
else if ( S_.y_[ State_::V_M ] >= B_.V_peak )

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Switch V_peak to V_ (see above)

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Switch V_peak to V_ (see above)

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

done

@golosio

golosio Oct 18, 2016

Contributor

done

{
S_.y_[ State_::V_M ] = P_.V_reset_;
S_.y_[ State_::W ] += P_.b; // spike-driven adaptation
@@ -575,7 +599,7 @@ aeif_cond_beta_multisynapse::update( Time const& origin,
}
}
for ( size_t i = 0; i < P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < P_.num_of_receptors; ++i )
{
S_.y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] +=
@@ -595,7 +619,7 @@ aeif_cond_beta_multisynapse::handles_test_event( SpikeEvent&,
rport receptor_type )
{
if ( receptor_type <= 0
|| receptor_type > static_cast< port >( P_.num_of_receptors_ ) )
|| receptor_type > static_cast< port >( P_.num_of_receptors ) )
{
throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" );
}
@@ -606,9 +630,15 @@ aeif_cond_beta_multisynapse::handles_test_event( SpikeEvent&,
void
aeif_cond_beta_multisynapse::handle( SpikeEvent& e )
{
if ( e.get_weight() < 0 )
{
throw BadProperty(
"Synaptic weights for conductance based models "
"must be positive." );
}
assert( e.get_delay() > 0 );
assert( ( e.get_rport() > 0 )
&& ( ( size_t ) e.get_rport() <= P_.num_of_receptors_ ) );
&& ( ( size_t ) e.get_rport() <= P_.num_of_receptors ) );
B_.spikes_[ e.get_rport() - 1 ].add_value(
e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
@@ -656,44 +686,37 @@ aeif_cond_beta_multisynapse_dynamics( double,
// good compiler will optimize the verbosity away ...
// shorthand for state variables
const double& V = y[ S::V_M ];
const double& V = std::min( y[ S::V_M ], node.P_.V_peak_ ); // bound V

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Maybe we should add a comment to assure people that we indeed want to use P_.V_peak_ and not V_.V_peak here...

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Maybe we should add a comment to assure people that we indeed want to use P_.V_peak_ and not V_.V_peak here...

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

done. However when I will commit please check if this is sufficient or please suggest a more suitable comment.

@golosio

golosio Oct 18, 2016

Contributor

done. However when I will commit please check if this is sufficient or please suggest a more suitable comment.

const double& w = y[ S::W ];
double I_syn = 0.0;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I would remove this empty line to tie the initialization more directly to the loop.

@heplesser

heplesser Oct 13, 2016

Contributor

I would remove this empty line to tie the initialization more directly to the loop.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

for ( size_t i = 0; i < ( node.P_.num_of_receptors_
* S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR );
i += S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR )
for ( size_t i = 0; i < node.P_.num_of_receptors; ++i )
{
double g = y[ S::G + i ];
if ( g > 0 )
{
I_syn += g * ( node.P_.E_ex - V ); // g>0, E_ex>V => I_syn increases
}
else
{
I_syn += g * ( V - node.P_.E_in ); // g<0, V>E_in => I_syn decreases
}
size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Add const.

@heplesser

heplesser Oct 13, 2016

Contributor

Add const.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

// I_syn = I_syn + g*(E_rev - V)

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Maybe move this comment above the I_syn = 0 initialization and make the summation explicit: I_syn = - sum_k g_k (V - E_rev_k).

@heplesser

heplesser Oct 13, 2016

Contributor

Maybe move this comment above the I_syn = 0 initialization and make the summation explicit: I_syn = - sum_k g_k (V - E_rev_k).

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

I_syn += y[ S::G + j ] * ( node.P_.E_rev[ i ] - V );
}
// We pre-compute the argument of the exponential
const double exp_arg = ( V - node.P_.V_th ) / node.P_.Delta_T;
// Upper bound for exponential argument to avoid numerical instabilities
const double MAX_EXP_ARG = 10.;
// If the argument is too large, we clip it.
const double I_spike =
node.P_.Delta_T * std::exp( std::min( exp_arg, MAX_EXP_ARG ) );
double I_spike;
if ( node.P_.Delta_T > 0. )
{
I_spike =
node.P_.Delta_T * std::exp( ( V - node.P_.V_th ) / node.P_.Delta_T );
}
else
{
I_spike = 0;
}
// dv/dt
f[ S::V_M ] = ( -node.P_.g_L * ( ( V - node.P_.E_L ) - I_spike ) + I_syn - w
f[ S::V_M ] = ( -node.P_.g_L * ( V - node.P_.E_L - I_spike ) + I_syn - w

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Switch to ( -node.P_.g_L * ( V - node.P_.E_L ) + I_spike + I_syn + ...

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Switch to ( -node.P_.g_L * ( V - node.P_.E_L ) + I_spike + I_syn + ...

This comment has been minimized.

@golosio

golosio Oct 18, 2016

Contributor

done

@golosio

golosio Oct 18, 2016

Contributor

done

+ node.P_.I_e + node.B_.I_stim_ ) / node.P_.C_m;
// Adaptation current w.
f[ S::W ] = ( node.P_.a * ( V - node.P_.E_L ) - w ) / node.P_.tau_w;
for ( size_t i = 0; i < node.P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < node.P_.num_of_receptors; ++i )
{
size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

const

@heplesser

heplesser Oct 13, 2016

Contributor

const

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

// Synaptic conductance derivative dG/dt
Oops, something went wrong.
ProTip! Use n and p to navigate between commits in a pull request.