Gif models #261

Merged
merged 38 commits into from Oct 7, 2016

Conversation

Projects
None yet
6 participants
@hesam-setareh
Contributor

hesam-setareh commented Mar 5, 2016

Implementing generalized integrate-and-fire (GIF) neuron model for nest in four versions:
gif_psc_exp
gif_psc_exp_multisynapse
gif_cond_exp
gif_cond_exp_multisynapse

And also a unittest for each version.

models/modelsmodule.cpp
+ kernel().model_manager.register_node_model< gif_psc_exp >( "gif_psc_exp" );
+ kernel().model_manager.register_node_model< gif_psc_exp_multisynapse >( "gif_psc_exp_multisynapse" );
+ kernel().model_manager.register_node_model< gif_cond_exp >( "gif_cond_exp" );
+ kernel().model_manager.register_node_model< gif_cond_exp_multisynapse >( "gif_cond_exp_multisynapse" );

This comment has been minimized.

@tammoippen

tammoippen Mar 7, 2016

Contributor

The cond models, which use the GSL have to be encapsulated by

#ifdef HAVE_GSL
...
#endif

Otherwise the template argument is not valid, when you compile without GSL.

@tammoippen

tammoippen Mar 7, 2016

Contributor

The cond models, which use the GSL have to be encapsulated by

#ifdef HAVE_GSL
...
#endif

Otherwise the template argument is not valid, when you compile without GSL.

@tammoippen

This comment has been minimized.

Show comment
Hide comment
@tammoippen

tammoippen Mar 7, 2016

Contributor

The formatting is off. Please follow the instructions on our developer space to format the code.

Contributor

tammoippen commented Mar 7, 2016

The formatting is off. Please follow the instructions on our developer space to format the code.

@jougs

This comment has been minimized.

Show comment
Hide comment
@jougs

jougs Mar 7, 2016

Contributor

Once the formal issues are resolved, I suggest @heplesser and @rcfduarte as reviewers.

Contributor

jougs commented Mar 7, 2016

Once the formal issues are resolved, I suggest @heplesser and @rcfduarte as reviewers.

@tammoippen

This comment has been minimized.

Show comment
Hide comment
@tammoippen

tammoippen Mar 8, 2016

Contributor

Hi @hesam-setareh ,

Which version of clang-format did you use? It is important that you use version 3.6, as the other versions will format the code differently.

Contributor

tammoippen commented Mar 8, 2016

Hi @hesam-setareh ,

Which version of clang-format did you use? It is important that you use version 3.6, as the other versions will format the code differently.

@apdavison apdavison referenced this pull request in pozzorin/GIFFittingToolbox Mar 22, 2016

Merged

Adapted the NEURON version of the GIF neuron model to work with PyNN #4

models/gif_cond_exp.h
+ instead of eta and gamma.
+
+
+ References:

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Please move the references section down, to just before the author section. Also, check the formatting of the references, the clang-format process seems to have introduced some unfortunate line breaks. I suppose this applies to the other h-files, too.

@heplesser

heplesser Mar 31, 2016

Contributor

Please move the references section down, to just before the author section. Also, check the formatting of the references, the clang-format process seems to have introduced some unfortunate line breaks. I suppose this applies to the other h-files, too.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_cond_exp.h
+ */
+
+#ifndef gif_cond_exp_H
+#define gif_cond_exp_H

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The include guard should be all caps: GIF_COND_EXP_H. I assume this applies to the other header files, too.

@heplesser

heplesser Mar 31, 2016

Contributor

The include guard should be all caps: GIF_COND_EXP_H. I assume this applies to the other header files, too.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.h
+ gamma_i = gamma_i + q_gamma_i (in case of spike emission).
+
+ In the source code and parameter names we use stc and sfa, respectively
+ instead of eta and gamma.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

This information does not really belong into the user space documentation, since it is only relevant for developers.

@heplesser

heplesser Mar 31, 2016

Contributor

This information does not really belong into the user space documentation, since it is only relevant for developers.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Moved to definition of Parameters_.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Moved to definition of Parameters_.

models/gif_psc_exp.h
+ (sfa) after each spike emission in mV.
+ tau_sfa vector of double - Time constants of sfa variables in ms.
+ delta_u double - Stochasticity level in mV.
+ lambda0 double - Stochastic intensity at firing threshold V_T in Hz.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Strictly speaking, Hz applies to periodic processes only, so it should be 1/s.

@heplesser

heplesser Mar 31, 2016

Contributor

Strictly speaking, Hz applies to periodic processes only, so it should be 1/s.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I think lambda_0 with underscore would be nicer.

@heplesser

heplesser Mar 31, 2016

Contributor

I think lambda_0 with underscore would be nicer.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I applied both comments.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I applied both comments.

models/gif_psc_exp.h
+ double_t y0_; //!< This is piecewise constant external current
+ double_t y3_; //!< This is the membrane potential RELATIVE TO RESTING POTENTIAL.
+ double_t q_; //!< This is the change of the 'threshold' due to adaptation.
+ double_t stc_; // Spike triggered current.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The comment (and subsequent similar comments) should begin with //!<, so that Doxygen will recognize it as member documentation when generating developer documentation. This would be a good place to point out that stc_ is eta in the paper.

@heplesser

heplesser Mar 31, 2016

Contributor

The comment (and subsequent similar comments) should begin with //!<, so that Doxygen will recognize it as member documentation when generating developer documentation. This would be a good place to point out that stc_ is eta in the paper.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+ , V_reset_( -55.0 ) // mV
+ , delta_u_( 1.5 ) // mV
+ , v_t_star_( -35 ) // mV
+ , lambda0_( 10000.0 ) // Hz

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

In the PLoS paper, you fix lambda_0 to 1 Hz. Why do you use 10.000 Hz as default here?

@heplesser

heplesser Mar 31, 2016

Contributor

In the PLoS paper, you fix lambda_0 to 1 Hz. Why do you use 10.000 Hz as default here?

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed the value to 1 Hz (1/s).

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed the value to 1 Hz (1/s).

models/gif_psc_exp.cpp
+ tau_sfa_.clear();
+ q_sfa_.clear();
+ tau_stc_.clear();
+ q_stc_.clear();

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Instead of calling clear() on all the vector members, it would be tidier to add them as , tau_sfa_() etc in the initializer list above.

@heplesser

heplesser Mar 31, 2016

Contributor

Instead of calling clear() on all the vector members, it would be tidier to add them as , tau_sfa_() etc in the initializer list above.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+ updateValue< double >( d, names::V_reset, V_reset_ );
+ updateValue< double >( d, names::delta_u, delta_u_ );
+ updateValue< double >( d, names::v_t_star, v_t_star_ );
+ updateValue< double >( d, "lambda0", lambda0_ );

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I would add lambda0 to names as well, but I think I would prefer an additional underscore: lambda_0

@heplesser

heplesser Mar 31, 2016

Contributor

I would add lambda0 to names as well, but I think I would prefer an additional underscore: lambda_0

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

You should probably convert lambda0 from 1/s to 1/ms (i.e., divide by 1000) here, then you do not need to do that on every time step. You need to multiply by 1000 then in the get() method. I suppose it would also make sense to test that lambda0 >= 0 (== 0 means no spiking, but it causes no problems as far as I can see, so one could allow it for testing).

@heplesser

heplesser Mar 31, 2016

Contributor

You should probably convert lambda0 from 1/s to 1/ms (i.e., divide by 1000) here, then you do not need to do that on every time step. You need to multiply by 1000 then in the get() method. I suppose it would also make sense to test that lambda0 >= 0 (== 0 means no spiking, but it causes no problems as far as I can see, so one could allow it for testing).

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+ {
+
+ q_temp_ = 0;
+ for ( uint_t i = 0; i < S_.q_stc_elems_.size(); i++ )

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Rather use size_t than uint_t for the loop variable.

@heplesser

heplesser Mar 31, 2016

Contributor

Rather use size_t than uint_t for the loop variable.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+ assert( to >= 0 && ( delay ) from < kernel().connection_builder_manager.get_min_delay() );
+ assert( from < to );
+
+ double_t q_temp_;

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

You could move the declaration a few lines down to first use.

@heplesser

heplesser Mar 31, 2016

Contributor

You could move the declaration a few lines down to first use.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Defining the declaration inside the loop reduces the performance. Therefore, I kept it outside the loop.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Defining the declaration inside the loop reduces the performance. Therefore, I kept it outside the loop.

models/gif_psc_exp.cpp
+ if ( S_.r_ref_ == 0 ) // neuron not refractory, so evolve V
+ {
+
+ if ( S_.add_stc_sfa_ == true )

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Drop the == true, if ( S_.add_stc_sfa_ ) does just the same.

@heplesser

heplesser Mar 31, 2016

Contributor

Drop the == true, if ( S_.add_stc_sfa_ ) does just the same.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

At this point, there is an interesting question about the precise definition of the model. In all other NEST models, refractoriness affects membrane potential only. In your implementation, though, you let stc (eta) and sfa (gamma) evolve during the refractory period as if the spike did not yet happen (lines 323-434), and then apply the jump to both values at the end of the refractory time. Eqs 1 and 4 in the PLoS paper gather give the impression that eta and gamma set in at the spike time t_j, not at t_j+t_ref.

Therefore, I would suggest to add the jumps right after a spike has been applied, in the code block that also emits the spike; you can then drop S_.add_stc_sfa_. One then needs to be a bit careful about the next time step:

  1. Let a time step be (t, t+h] (time steps in NEST are always left open, right closed intervals). The stc_elems_ at the end of a time step contain the value at the end of the step, including a jump if the neuron spiked during the step.
  2. We need to use those stc_elems_ values when computing the total stc entering into the membrane potential update.
  3. Only after we have used the values must we update them.

You do this already correctly on lines 323-343 (the same applies to sfa)

The only problem is the sfa value returned by State_::get(), since that value is only updated at the beginning of the next update loop and would just not contain any jumps occurring during the last step of a Simulate call. The best way to make sure that the correct sfa value is returned, is probably to compute the sfa value explicity in the get() method.

@heplesser

heplesser Mar 31, 2016

Contributor

At this point, there is an interesting question about the precise definition of the model. In all other NEST models, refractoriness affects membrane potential only. In your implementation, though, you let stc (eta) and sfa (gamma) evolve during the refractory period as if the spike did not yet happen (lines 323-434), and then apply the jump to both values at the end of the refractory time. Eqs 1 and 4 in the PLoS paper gather give the impression that eta and gamma set in at the spike time t_j, not at t_j+t_ref.

Therefore, I would suggest to add the jumps right after a spike has been applied, in the code block that also emits the spike; you can then drop S_.add_stc_sfa_. One then needs to be a bit careful about the next time step:

  1. Let a time step be (t, t+h] (time steps in NEST are always left open, right closed intervals). The stc_elems_ at the end of a time step contain the value at the end of the step, including a jump if the neuron spiked during the step.
  2. We need to use those stc_elems_ values when computing the total stc entering into the membrane potential update.
  3. Only after we have used the values must we update them.

You do this already correctly on lines 323-343 (the same applies to sfa)

The only problem is the sfa value returned by State_::get(), since that value is only updated at the beginning of the next update loop and would just not contain any jumps occurring during the last step of a Simulate call. The best way to make sure that the correct sfa value is returned, is probably to compute the sfa value explicity in the get() method.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Probably this is the most controversial point of the implementation. I have discussed with the designer of gif model and he pointed out:
1 - stc and sfa should be evolved during the refractory period.
2 - In case of spike emission, the jumps should be applied to stc and sfa after the refractory period (t_j + t_ref) not after the spike (t_j). That is why I used S_.add_stc_sfa_. true value of this variable indicates that there was a spike before the refractory period and stc and sfa should be jumped.

In the modified version of the code I added stc to the recordable map, so the user can capture the value of both adaptation variables (stc and sfa).

I agree with you that the article does not explain this point clearly.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Probably this is the most controversial point of the implementation. I have discussed with the designer of gif model and he pointed out:
1 - stc and sfa should be evolved during the refractory period.
2 - In case of spike emission, the jumps should be applied to stc and sfa after the refractory period (t_j + t_ref) not after the spike (t_j). That is why I used S_.add_stc_sfa_. true value of this variable indicates that there was a spike before the refractory period and stc and sfa should be jumped.

In the modified version of the code I added stc to the recordable map, so the user can capture the value of both adaptation variables (stc and sfa).

I agree with you that the article does not explain this point clearly.

models/gif_psc_exp.h
+ 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_u double - Stochasticity level in mV.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

This corresponds to Delta_V in the paper, doesn't it? Wouldn't it be better then to call it Delta_V?

@heplesser

heplesser Mar 31, 2016

Contributor

This corresponds to Delta_V in the paper, doesn't it? Wouldn't it be better then to call it Delta_V?

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed it to Delta_V.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed it to Delta_V.

models/gif_psc_exp.h
+ tau_sfa vector of double - Time constants of sfa variables in ms.
+ delta_u double - Stochasticity level in mV.
+ lambda0 double - Stochastic intensity at firing threshold V_T in Hz.
+ v_t_star double - Minimum threshold in mV

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I'd suggest V_T_start to make it more similar to the paper.

@heplesser

heplesser Mar 31, 2016

Contributor

I'd suggest V_T_start to make it more similar to the paper.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+
+ S_.q_ = q_temp_ + P_.v_t_star_;
+
+ ulong_t n_spikes = 0;

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Could the neuron ever produce more that one spike in a time step? If not, I would make this a bool called neuron_spikes_ or similar.

@heplesser

heplesser Mar 31, 2016

Contributor

Could the neuron ever produce more that one spike in a time step? If not, I would make this a bool called neuron_spikes_ or similar.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You have right. The neuron does not produce more than one spike. Actually, I did not need the neuron_spikes_ variable and removed it.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You have right. The neuron does not produce more than one spike. Actually, I did not need the neuron_spikes_ variable and removed it.

models/gif_psc_exp.h
+ {
+ double_t y0_; //!< This is piecewise constant external current
+ double_t y3_; //!< This is the membrane potential RELATIVE TO RESTING POTENTIAL.
+ double_t q_; //!< This is the change of the 'threshold' due to adaptation.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I think you are overusing q_ a bit, which makes the code more difficult to read than necessary. Above, q_* is used for the jump amplitude for stc and sfa, here you use it for the threshold. I would choose a different name here---why not sfa_?

@heplesser

heplesser Mar 31, 2016

Contributor

I think you are overusing q_ a bit, which makes the code more difficult to read than necessary. Above, q_* is used for the jump amplitude for stc and sfa, here you use it for the threshold. I would choose a different name here---why not sfa_?

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed it to sfa_.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed it to sfa_.

models/gif_psc_exp.h
+ double_t stc_; // Spike triggered current.
+
+ std::vector< double_t > q_sfa_elems_; // Vector of adaptation parameters.
+ std::vector< double_t > q_stc_elems_; // Vector of spike triggered parameters.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The q_ prefix in these two variables does not make much sense, I think, given that you use q_ to indicate jump sizes in the parameter names. I'd just drop q_ here.

@heplesser

heplesser Mar 31, 2016

Contributor

The q_ prefix in these two variables does not make much sense, I think, given that you use q_ to indicate jump sizes in the parameter names. I'd just drop q_ here.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I modified the name of variables. Now we have sfa_ and stc_.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I modified the name of variables. Now we have sfa_ and stc_.

models/gif_psc_exp.h
+
+ int_t r_ref_; // absolute refractory counter (no membrane potential propagation)
+
+ bool initialized_; // it is true if the vectors are initialized

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The name should be more specific. e.g., sfa_stc_initialized_, but it might be unnecessary, see comments in cpp-file.

@heplesser

heplesser Mar 31, 2016

Contributor

The name should be more specific. e.g., sfa_stc_initialized_, but it might be unnecessary, see comments in cpp-file.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed its name to sfa_stc_initialized_. I think we need this variable. I explain the reason in the comment of cpp-file.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed its name to sfa_stc_initialized_. I think we need this variable. I explain the reason in the comment of cpp-file.

models/gif_psc_exp.cpp
+ S_.i_syn_in_ *= V_.P11in_;
+
+ S_.i_syn_ex_ += B_.spikes_ex_.get_value( lag );
+ S_.i_syn_in_ += B_.spikes_in_.get_value( lag );

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

L 401-405: You need to update the synaptic currents before computing the new membrane potential S_.y3_ above, and at the same time make sure that synaptic currents are updated also while the neuron is refractory.

@heplesser

heplesser Mar 31, 2016

Contributor

L 401-405: You need to update the synaptic currents before computing the new membrane potential S_.y3_ above, and at the same time make sure that synaptic currents are updated also while the neuron is refractory.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Fixed.

models/gif_psc_exp.h
+ double_t P21ex_;
+ double_t P21in_;
+ std::vector< double_t > Q33_; // for sfa
+ std::vector< double_t > Q44_; // for stc

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I think P_sfa_ and P_stc_ would be more useful names, making clear to which quantities these propagator elements apply.

@heplesser

heplesser Mar 31, 2016

Contributor

I think P_sfa_ and P_stc_ would be more useful names, making clear to which quantities these propagator elements apply.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed their names to P_sfa_ and P_stc.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

I changed their names to P_sfa_ and P_stc.

models/gif_psc_exp.h
+ std::vector< double_t > Q44_; // for stc
+
+
+ double_t h_; //!< simulation time step in ms

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

No need to store this in a member variable, you can always get it from the kernel when needed using efficient inlined.

@heplesser

heplesser Mar 31, 2016

Contributor

No need to store this in a member variable, you can always get it from the kernel when needed using efficient inlined.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

It has been removed.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

It has been removed.

models/gif_psc_exp.h
+
+ librandom::RngPtr rng_; // random number generator of my own thread
+ librandom::PoissonRandomDev poisson_dev_; // random deviate generator
+ librandom::GammaRandomDev gamma_dev_; // random deviate generator

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Do you use the poisson and gamma generators anywhere?

@heplesser

heplesser Mar 31, 2016

Contributor

Do you use the poisson and gamma generators anywhere?

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You have right. We did not need them. They are removed in the modified version.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You have right. We did not need them. They are removed in the modified version.

models/gif_psc_exp.cpp
+ {
+
+ // Draw random number and compare to prob to have a spike
+ if ( V_.rng_->drand() <= -numerics::expm1( -lambda * ( V_.h_ / 1000.0 ) ) )

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

This should be <, not <=, since drand() returns numbers from [0, 1), not (0, 1]. You can avoid repeated division by 1000 by rescaling lambda0 to 1/ms when the parameter is set. You should also comment why the expm1() function appears here. Most people would expect just lambda*h, but you approximate the integral in Eq 3 of the PLoS paper as lambda * h and then use the more exact expression 1 - exp(-lambda*h), don't you?

@heplesser

heplesser Mar 31, 2016

Contributor

This should be <, not <=, since drand() returns numbers from [0, 1), not (0, 1]. You can avoid repeated division by 1000 by rescaling lambda0 to 1/ms when the parameter is set. You should also comment why the expm1() function appears here. Most people would expect just lambda*h, but you approximate the integral in Eq 3 of the PLoS paper as lambda * h and then use the more exact expression 1 - exp(-lambda*h), don't you?

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor
  1. I have changed <= to <.
  2. I rescaled the value of lambda_0 to 1/ms, then dropped the division by 1000.
  3. Thats is right. I use the exact formula 1 - exp (-lambda*h). I have added a line of comment to clarify: // hazard function is computed by 1 - exp(- lambda * dt)
@hesam-setareh

hesam-setareh May 17, 2016

Contributor
  1. I have changed <= to <.
  2. I rescaled the value of lambda_0 to 1/ms, then dropped the division by 1000.
  3. Thats is right. I use the exact formula 1 - exp (-lambda*h). I have added a line of comment to clarify: // hazard function is computed by 1 - exp(- lambda * dt)
models/gif_psc_exp.cpp
+{
+ updateValue< double >( d, names::V_m, y3_ );
+ updateValue< double >( d, names::E_sfa, q_ );
+ initialized_ = false; // vectors of the state should be initialized with new parameter set.

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I think it makes no sense to set State_::q_. Unless I am mistaken, q_ is unconditionally computed from scratch at the beginning of each update loop (line 343), so setting it will have no effect. Then, you can also drop the initialized_ variable.

@heplesser

heplesser Mar 31, 2016

Contributor

I think it makes no sense to set State_::q_. Unless I am mistaken, q_ is unconditionally computed from scratch at the beginning of each update loop (line 343), so setting it will have no effect. Then, you can also drop the initialized_ variable.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You are right. I removed updateValue< double >( d, names::E_sfa, q_ ). However, we need the variable sfa_stc_initialized_. Here is the story:
When a new set of parameter is given by the user, we need to reconstruct arrays of P_sfa_, P_stc_, stc_elems_ and sfa_elems in the calibrate based on the new values of sfa and stc parameters. However, we cannot reconstruct them each time calibrate is called. For example, in the below situation

Simulate(10.0)
Simulate(10.0)
Simulate(10.0)

the arrays should be reconstructed only in the calibration of first Simulate. But in the case of

Simulate(10.0)
SetStatus(...)
Simulate(10.0)

the arrays should be reconstructed in the calibration of both. Therefore, we need to inform the calibrate when it should do the reconstruction. The variable sfa_stc_initialized_ keeps this information for us. When a new parameter set is given by the user, this variable is set to false. Note that we cannon change its value in Parameters_::set (non-static reference in a static method). Therefore we do that in State_::set.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

You are right. I removed updateValue< double >( d, names::E_sfa, q_ ). However, we need the variable sfa_stc_initialized_. Here is the story:
When a new set of parameter is given by the user, we need to reconstruct arrays of P_sfa_, P_stc_, stc_elems_ and sfa_elems in the calibrate based on the new values of sfa and stc parameters. However, we cannot reconstruct them each time calibrate is called. For example, in the below situation

Simulate(10.0)
Simulate(10.0)
Simulate(10.0)

the arrays should be reconstructed only in the calibration of first Simulate. But in the case of

Simulate(10.0)
SetStatus(...)
Simulate(10.0)

the arrays should be reconstructed in the calibration of both. Therefore, we need to inform the calibrate when it should do the reconstruction. The variable sfa_stc_initialized_ keeps this information for us. When a new parameter set is given by the user, this variable is set to false. Note that we cannon change its value in Parameters_::set (non-static reference in a static method). Therefore we do that in State_::set.

models/gif_psc_exp.cpp
+{
+ const gif_psc_exp& pr = downcast< gif_psc_exp >( proto );
+ S_ = pr.S_;
+ // S_.r_ref_ = Time(Time::ms(P_.t_ref_remaining_)).get_steps();

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

You should delete dead code.

@heplesser

heplesser Mar 31, 2016

Contributor

You should delete dead code.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Done.

models/gif_psc_exp.cpp
+ for ( uint_t i = 0; i < P_.tau_sfa_.size(); i++ )
+ {
+ V_.Q33_.push_back( std::exp( -V_.h_ / P_.tau_sfa_[ i ] ) );
+ S_.q_sfa_elems_.push_back( 0.0 );

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The sfa_elems_ and stc_elems_ are parts of the neuron's state. They should therefore only be set in the State_ constructor or possibly init_state_(). Here, only the V_.Q*_ variables should be initialized here. That must happen unconditionally.

@heplesser

heplesser Mar 31, 2016

Contributor

The sfa_elems_ and stc_elems_ are parts of the neuron's state. They should therefore only be set in the State_ constructor or possibly init_state_(). Here, only the V_.Q*_ variables should be initialized here. That must happen unconditionally.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

The sfa_elems_ depends on Q33_ (P_sfa_ in the new version). They should have same size. Each time P_sfa_ gets new values, we should reset the values of sfa_elems_. A similar relation exists between stc_elems_ and P_stc_. Therefore, I decided to reconstruct the sfa_elems_ and stc_elems_, after calibration of P_sfa_ and P_stc_.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

The sfa_elems_ depends on Q33_ (P_sfa_ in the new version). They should have same size. Each time P_sfa_ gets new values, we should reset the values of sfa_elems_. A similar relation exists between stc_elems_ and P_stc_. Therefore, I decided to reconstruct the sfa_elems_ and stc_elems_, after calibration of P_sfa_ and P_stc_.

models/gif_psc_exp.cpp
+ // Draw random number and compare to prob to have a spike
+ if ( V_.rng_->drand() <= -numerics::expm1( -lambda * ( V_.h_ / 1000.0 ) ) )
+ {
+ n_spikes = 1;

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

The code emitting the spike should go here, no need to place spike emission further down.

@heplesser

heplesser Mar 31, 2016

Contributor

The code emitting the spike should go here, no need to place spike emission further down.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

Fixed.

models/gif_psc_exp.cpp
+nest::gif_psc_exp::State_::get( DictionaryDatum& d, const Parameters_& p ) const
+{
+ def< double >( d, names::V_m, y3_ ); // Membrane potential
+ def< double >( d, names::E_sfa, q_ ); // Adaptive threshold potential

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

To get correct values even after spikes during the last time step of a Simulate call, the sfa value should be recomputed explicitly here from the sfa_elems. Maybe it would be even better to return the resulting threshold value, including V_T_star.

@heplesser

heplesser Mar 31, 2016

Contributor

To get correct values even after spikes during the last time step of a Simulate call, the sfa value should be recomputed explicitly here from the sfa_elems. Maybe it would be even better to return the resulting threshold value, including V_T_star.

This comment has been minimized.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

The value of q_ (sfa_ in the new version) includes V_T_star by default.

@hesam-setareh

hesam-setareh May 17, 2016

Contributor

The value of q_ (sfa_ in the new version) includes V_T_star by default.

testsuite/unittests/test_gif_psc_exp.sli
+
+
+/gif_psc_exp Create /neuron Set
+neuron << /tau_syn_ex 30.0 /tau_sfa [120.0] /q_sfa [20.0] /tau_stc [10.0] /q_stc [20.0] >> SetStatus

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

I think it would be good to have a test as well that uses more than one sfa and stc component, to make sure that the code for multiple elements also works. And one test with the default values, to make sure that that works.

@heplesser

heplesser Mar 31, 2016

Contributor

I think it would be good to have a test as well that uses more than one sfa and stc component, to make sure that the code for multiple elements also works. And one test with the default values, to make sure that that works.

This comment has been minimized.

@hesam-setareh

hesam-setareh Jul 18, 2016

Contributor

I have changed the tests. Now it tests two neurons, one with default parameters and another one with two sfa and one stc components.

@hesam-setareh

hesam-setareh Jul 18, 2016

Contributor

I have changed the tests. Now it tests two neurons, one with default parameters and another one with two sfa and one stc components.

testsuite/unittests/test_gif_psc_exp.sli
+ 556 /t2 Set
+ 2046 /t3 Set
+}
+ifelse

This comment has been minimized.

@heplesser

heplesser Mar 31, 2016

Contributor

Where do these values come from? As the dependence on the RNG indicates, they seem implementation dependent and thus rather test that a given NEST implementation behaves the same now as previously. Is there a way to test that the "mechanics" of the model work correctly that would work even if we exchanged the RNG, without resorting to different sets of "magic numbers"?

@heplesser

heplesser Mar 31, 2016

Contributor

Where do these values come from? As the dependence on the RNG indicates, they seem implementation dependent and thus rather test that a given NEST implementation behaves the same now as previously. Is there a way to test that the "mechanics" of the model work correctly that would work even if we exchanged the RNG, without resorting to different sets of "magic numbers"?

This comment has been minimized.

@hesam-setareh

hesam-setareh Jul 18, 2016

Contributor

These values are from original implementation (python version) of gif model used for the publications. Unfortunately, systematic test of such stochastic model is very difficult. We need to run the simulation for many trials, calculate the average firing rate and compare it to the theoretical expectations. Running simulation for a long time and the lack of theoretical calculation for the gif model are the main obstacles.

@hesam-setareh

hesam-setareh Jul 18, 2016

Contributor

These values are from original implementation (python version) of gif model used for the publications. Unfortunately, systematic test of such stochastic model is very difficult. We need to run the simulation for many trials, calculate the average firing rate and compare it to the theoretical expectations. Running simulation for a long time and the lack of theoretical calculation for the gif model are the main obstacles.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Mar 31, 2016

Contributor

@hesam-setareh Thanks a lot for your PR! I have now started reviewing and I suggest that we work in gif_psc_exp first before addressing the other models. The code seems generally sound, but I left quite a number of comments in the code to make it even better; I hope they make sense. One more little detail: there are quite many superfluous empty lines in the code, maybe you could tidy them up as well?

Contributor

heplesser commented Mar 31, 2016

@hesam-setareh Thanks a lot for your PR! I have now started reviewing and I suggest that we work in gif_psc_exp first before addressing the other models. The code seems generally sound, but I left quite a number of comments in the code to make it even better; I hope they make sense. One more little detail: there are quite many superfluous empty lines in the code, maybe you could tidy them up as well?

models/gif_cond_exp.cpp
+ , sfa_elems_()
+ , stc_elems_()
+ , i_syn_ex_( 0.0 )
+ , i_syn_in_( 0.0 )

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

I think these currents are not needed (see same comment in gif_cond_exp.h)

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

I think these currents are not needed (see same comment in gif_cond_exp.h)

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 12, 2016

Contributor

i_syn_ex_ and i_syn_in_ are removed.

@hesam-setareh

hesam-setareh Sep 12, 2016

Contributor

i_syn_ex_ and i_syn_in_ are removed.

models/gif_psc_exp.h
+#include "nest.h"
+
+/* BeginDocumentation
+ Name: gif_psc_exp - Current based generalized integrate-and-fire neuron

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

"Current**-**based", cf. same comments in other header files

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

"Current**-**based", cf. same comments in other header files

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 12, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh Sep 12, 2016

Contributor

Fixed.

models/gif_psc_exp.h
+ {
+ double_t y0_; //!< This is piecewise constant external current
+ double_t
+ y3_; //!< This is the membrane potential RELATIVE TO RESTING POTENTIAL.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

Could you replace y0_ and y3_ by I_stim_ and V_rel_ for clarity?

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

Could you replace y0_ and y3_ by I_stim_ and V_rel_ for clarity?

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I changed them to I_stimand V_. I just noticed that y3_ was the actual voltage not the relative voltage.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I changed them to I_stimand V_. I just noticed that y3_ was the actual voltage not the relative voltage.

models/gif_psc_exp.cpp
+ insert_( names::E_sfa, &gif_psc_exp::get_E_sfa_ );
+ insert_( names::I_stc, &gif_psc_exp::get_I_stc_ );
+ insert_( names::input_currents_ex, &gif_psc_exp::get_input_currents_ex_ );
+ insert_( names::input_currents_in, &gif_psc_exp::get_input_currents_in_ );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

In nest_names.h, I_ex is referenced explicitly as synaptic current whereas input_currents_ex is not... @heplesser do you have an opinion about this?

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

In nest_names.h, I_ex is referenced explicitly as synaptic current whereas input_currents_ex is not... @heplesser do you have an opinion about this?

This comment has been minimized.

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

Well, this now needs to change because of issue #501 and PR #502: can you replace input_currents_ex/in by I_syn_ex/in? input_currents_ex is going to disappear from nest_names.

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

Well, this now needs to change because of issue #501 and PR #502: can you replace input_currents_ex/in by I_syn_ex/in? input_currents_ex is going to disappear from nest_names.

This comment has been minimized.

@heplesser

heplesser Oct 6, 2016

Contributor

hesam-setareh After @Silmathoron had discovered that we had different names for the same thing in different models, we decided on a standardization in the Developer VC on 3 Oct. This model should also use the standard names names::I_syn_ex, names::I_syn_in.

@heplesser

heplesser Oct 6, 2016

Contributor

hesam-setareh After @Silmathoron had discovered that we had different names for the same thing in different models, we decided on a standardization in the Developer VC on 3 Oct. This model should also use the standard names names::I_syn_ex, names::I_syn_in.

models/gif_psc_exp.cpp
+void
+nest::gif_psc_exp::State_::get( DictionaryDatum& d, const Parameters_& p ) const
+{
+ def< double >( d, names::V_m, y3_ ); // Membrane potential

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

I thought from the header file that y3_ was supposed to be the relative potential, not the absolute V_m... isn't there a problem here?

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

I thought from the header file that y3_ was supposed to be the relative potential, not the absolute V_m... isn't there a problem here?

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

You are right. y3_ (which is called V_ now) is actual voltage not the relative voltage. I changed the comment in the definition of State_.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

You are right. y3_ (which is called V_ now) is actual voltage not the relative voltage. I changed the comment in the definition of State_.

models/gif_psc_exp.cpp
+void
+nest::gif_psc_exp::State_::set( const DictionaryDatum& d, const Parameters_& p )
+{
+ updateValue< double >( d, names::V_m, y3_ );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

Same question as for line 204

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

Same question as for line 204

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

You are right. y3_ (which is called V_ now) is actual voltage not the relative voltage. I changed the comment in the definition of State_.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

You are right. y3_ (which is called V_ now) is actual voltage not the relative voltage. I changed the comment in the definition of State_.

models/gif_psc_exp.h
+ double_t P11ex_;
+ double_t P11in_;
+ double_t P21ex_;
+ double_t P21in_;

This comment has been minimized.

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

It would be nice to have the propagators documented

@Silmathoron

Silmathoron Sep 7, 2016

Contributor

It would be nice to have the propagators documented

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I added a line of comment for each variable.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I added a line of comment for each variable.

+ On the postsynapic side, there can be arbitrarily many synaptic time constants
+ (gif_psc_exp has exactly two: tau_syn_ex and tau_syn_in). This can be reached
+ by specifying separate receptor ports, each for a different time constant. The
+ port number has to match the respective "receptor_type" in the connectors.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

I think it would be a good thing to develop this part a bit, perhaps giving a quick example like:

''' Create neuron, voltmeter and poisson_generator '''
n = nest.Create("gif_cond_exp_multisynapse", params={"taus_syn": [1.3,  4.5, 0.8]}) # 3 synapses/ports
pg = nest.Create("poisson_generator", 3, params={"rate": 10000.}) # poisson generator for each synapse
vm = nest.Create("voltmeter")
''' Connect voltmeter ''''
nest.Connect(vm, n, syn_spec={"receptor_type": 0})
''' Connect poisson generator '''
nest.Connect([pg[0]], n, syn_spec={"receptor_type": 1})
nest.Connect([pg[1]], n, syn_spec={"receptor_type": 2})
nest.Connect([pg[2]], n, syn_spec={"receptor_type": 3})

This would make the description self-sufficient for model use...
@Silmathoron

Silmathoron Sep 8, 2016

Contributor

I think it would be a good thing to develop this part a bit, perhaps giving a quick example like:

''' Create neuron, voltmeter and poisson_generator '''
n = nest.Create("gif_cond_exp_multisynapse", params={"taus_syn": [1.3,  4.5, 0.8]}) # 3 synapses/ports
pg = nest.Create("poisson_generator", 3, params={"rate": 10000.}) # poisson generator for each synapse
vm = nest.Create("voltmeter")
''' Connect voltmeter ''''
nest.Connect(vm, n, syn_spec={"receptor_type": 0})
''' Connect poisson generator '''
nest.Connect([pg[0]], n, syn_spec={"receptor_type": 1})
nest.Connect([pg[1]], n, syn_spec={"receptor_type": 2})
nest.Connect([pg[2]], n, syn_spec={"receptor_type": 3})

This would make the description self-sufficient for model use...

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I think it is better that we create a python example file and show how to use a multisynapse model there. What is your opinion?

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I think it is better that we create a python example file and show how to use a multisynapse model there. What is your opinion?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 14, 2016

Contributor

Sounds reasonable, and you could then put a link in the description.

@Silmathoron

Silmathoron Sep 14, 2016

Contributor

Sounds reasonable, and you could then put a link in the description.

This comment has been minimized.

@heplesser

heplesser Sep 16, 2016

Contributor

I fully agree, and example would be very useful and better than lengthy documentation.

@heplesser

heplesser Sep 16, 2016

Contributor

I fully agree, and example would be very useful and better than lengthy documentation.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

How did things evolve regarding this "example part"?

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

How did things evolve regarding this "example part"?

This comment has been minimized.

@heplesser

heplesser Oct 6, 2016

Contributor

@hesam-setareh @Silmathoron While I still think that adding an example would be very useful, I think we can proceed without it at the moment, since we are going to change the logic for incoming weights and reversal potentials according to the Developer VC 19 Sep decision, see also comment on reversal potentials below.

@heplesser

heplesser Oct 6, 2016

Contributor

@hesam-setareh @Silmathoron While I still think that adding an example would be very useful, I think we can proceed without it at the moment, since we are going to change the logic for incoming weights and reversal potentials according to the Developer VC 19 Sep decision, see also comment on reversal potentials below.

models/gif_cond_exp_multisynapse.h
+ Receives: SpikeEvent, CurrentEvent, DataLoggingRequest
+
+ Author: March 2016, Setareh
+ SeeAlso: pp_psc_delta, gif_cond_exp, gif_cond_exp, gif_cond_exp_multisynapse

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

gif_cond_exp is referenced twice + self reference.
Maybe

SeeAlso: pp_psc_delta, gif_cond_exp, iaf_psc_exp_multisynapse, gif_psc_exp_multisynapse

would be more appropriate.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

gif_cond_exp is referenced twice + self reference.
Maybe

SeeAlso: pp_psc_delta, gif_cond_exp, iaf_psc_exp_multisynapse, gif_psc_exp_multisynapse

would be more appropriate.

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed. Thanks for your attention.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed. Thanks for your attention.

+ taus_syn vector of double - Time constants of the synaptic conductance
+ in ms.
+ E_ex double - Excitatory reversal potential in mV.
+ E_in double - Inhibitory reversal potential in mV.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Why are these potentials double and not vector of double? Wouldn't it make sense to have different reversal potentials for each synapse?

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Why are these potentials double and not vector of double? Wouldn't it make sense to have different reversal potentials for each synapse?

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I looked at the previous multisynapse models and followed their implementations. That is why I used only two reversal potentials. If it is necessary, I can use an array instead of two values.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

I looked at the previous multisynapse models and followed their implementations. That is why I used only two reversal potentials. If it is necessary, I can use an array instead of two values.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 14, 2016

Contributor

I mean, it's not required or anything, it's just that I think it would make more sense to me... and we should not always do things a certain way because they were done that way before ^^

@Silmathoron

Silmathoron Sep 14, 2016

Contributor

I mean, it's not required or anything, it's just that I think it would make more sense to me... and we should not always do things a certain way because they were done that way before ^^

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

That is right. I am actually on your side here :)
I think Hans will discuss about this issue in the next developer meeting. I will wait for the decision and then apply the outcome.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

That is right. I am actually on your side here :)
I think Hans will discuss about this issue in the next developer meeting. I will wait for the decision and then apply the outcome.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

@heplesser what was the outcome on that question?

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

@heplesser what was the outcome on that question?

This comment has been minimized.

@heplesser

heplesser Oct 6, 2016

Contributor

@hesam-setareh @Silmathoron Sorry for not updating on this earlier. In the Developer VC 19 Sep we decided that for conductance based neurons, each rport should have one input channel and one reversal potential. So instead of E_ex and E_in, where should be a vector of reversal potentials the same size as taus_syn.

BUT: I just recently merged #439 and #474 with multisynapse models that do not implement the new logic either. So I would not consider this a blocker for this PR either. I will create new issues and we should then try to fix this through new PRs asap.

@heplesser

heplesser Oct 6, 2016

Contributor

@hesam-setareh @Silmathoron Sorry for not updating on this earlier. In the Developer VC 19 Sep we decided that for conductance based neurons, each rport should have one input channel and one reversal potential. So instead of E_ex and E_in, where should be a vector of reversal potentials the same size as taus_syn.

BUT: I just recently merged #439 and #474 with multisynapse models that do not implement the new logic either. So I would not consider this a blocker for this PR either. I will create new issues and we should then try to fix this through new PRs asap.

models/gif_cond_exp_multisynapse.h
+
+ std::vector< double_t > y_; //!< neuron state
+
+ double_t y0_; //!< This is piecewise constant external current

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Better to name it I_stim_ (cf. gif_cond_exp.h).

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Better to name it I_stim_ (cf. gif_cond_exp.h).

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed.

+ &gif_cond_exp_multisynapse::
+ get_y_elem_< gif_cond_exp_multisynapse::State_::V_M > );
+ insert_( names::E_sfa, &gif_cond_exp_multisynapse::get_E_sfa_ );
+ insert_( names::I_stc, &gif_cond_exp_multisynapse::get_I_stc_ );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

@heplesser can RecordablesMap deal with vectors? (i.e. can we create a get_gs function that would return all synaptic conductances?)
I think it would be best to give access to as much synaptic information as possible...
Otherwise we could at least give access to I_syn_exc and I_syn_inh...

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

@heplesser can RecordablesMap deal with vectors? (i.e. can we create a get_gs function that would return all synaptic conductances?)
I think it would be best to give access to as much synaptic information as possible...
Otherwise we could at least give access to I_syn_exc and I_syn_inh...

This comment has been minimized.

@heplesser

heplesser Sep 8, 2016

Contributor

@Silmathoron RecordablesMap can only map to functions returning a double, no vectors.

@heplesser

heplesser Sep 8, 2016

Contributor

@Silmathoron RecordablesMap can only map to functions returning a double, no vectors.

models/gif_cond_exp_multisynapse.cpp
+ if ( tau_tmp[ i ] == ( c_m_ / g_L_ ) )
+ throw BadProperty(
+ "Membrane and synapse time constant(s) must differ. See note in "
+ "documentation." );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Which "note" are you referencing? I'm afraid I did not find it...

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Which "note" are you referencing? I'm afraid I did not find it...

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

As Hans mentioned, we do not need to check this condition (propagator checks this condition itself). Therefore, I removed these lines.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

As Hans mentioned, we do not need to check this condition (propagator checks this condition itself). Therefore, I removed these lines.

+ def< double >( d, names::t_ref, t_ref_ );
+ def< double >( d, names::E_ex, E_ex_ );
+ def< double >( d, names::E_in, E_in_ );
+ def< int >( d, "n_synapses", num_of_receptors_ );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

@heplesser I am wondering about this n_synapses... I know it is what is used in the original iaf_psc_alpha_multisynapse model, but maybe we should change this to a more explicit num_of_synapses name.
And I guess this should go into nest_names.h...

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

@heplesser I am wondering about this n_synapses... I know it is what is used in the original iaf_psc_alpha_multisynapse model, but maybe we should change this to a more explicit num_of_synapses name.
And I guess this should go into nest_names.h...

This comment has been minimized.

@heplesser

heplesser Sep 8, 2016

Contributor

@Silmathoron We use n_events in the dict returned by the spike detector, so I think n_synapses is okay. num_synapses would be slightly more expressive, but not enough to warrant a change, I think. Putting n_synapses into nest_names.h would be a good thing, but that applies to other neuron models as well, so I would create a separate issue for it and let this pass for now.

@heplesser

heplesser Sep 8, 2016

Contributor

@Silmathoron We use n_events in the dict returned by the spike detector, so I think n_synapses is okay. num_synapses would be slightly more expressive, but not enough to warrant a change, I think. Putting n_synapses into nest_names.h would be a good thing, but that applies to other neuron models as well, so I would create a separate issue for it and let this pass for now.

models/gif_psc_exp_multisynapse.h
+ {
+ double_t y0_; //!< This is piecewise constant external current
+ double_t
+ y3_; //!< This is the membrane potential RELATIVE TO RESTING POTENTIAL.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Change y0_ and y3_ (see gif_psc_exp.h)

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Change y0_ and y3_ (see gif_psc_exp.h)

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Done.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Done.

nestkernel/nest_names.h
+extern const Name b; //!< Specific to Brette & Gerstner 2005 (aeif_cond-*)
+extern const Name beta; //!< Specific to amat2_*
+extern const Name
+ beta_Ca; //!< Increment in calcuim concentration with each spike

This comment has been minimized.

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

"calcium" (I know it's not you, but since we are at it... ^^")

@Silmathoron

Silmathoron Sep 8, 2016

Contributor

"calcium" (I know it's not you, but since we are at it... ^^")

This comment has been minimized.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh Sep 14, 2016

Contributor

Fixed.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 8, 2016

Contributor

Ok, I've gone through the code and it looks like it'll be great once a few details are corrected.
In addition, it's really nice to have a new example file!

Contributor

Silmathoron commented Sep 8, 2016

Ok, I've gone through the code and it looks like it'll be great once a few details are corrected.
In addition, it's really nice to have a new example file!

Merge pull request #3 from nest/master
Remove _t from datatypes
@hesam-setareh

This comment has been minimized.

Show comment
Hide comment
@hesam-setareh

hesam-setareh Sep 15, 2016

Contributor

@heplesser , @Silmathoron
Dear Hans and Tanguy,

Thanks for your useful comments. I have applied the comments and modified GIF models. After deciding about the number of reversal potentials in multisynapse conductance-based models, I will do more things:

  • Modifying gif_cond_exp_multisynapse in order to change the number of reversal potentials
  • Adding an example of multisyanpse model (ipython format)
  • Adding a SLI testunit to check consistency for each pair of multisynapse/plain neuron model

Regards,
Hesam

Contributor

hesam-setareh commented Sep 15, 2016

@heplesser , @Silmathoron
Dear Hans and Tanguy,

Thanks for your useful comments. I have applied the comments and modified GIF models. After deciding about the number of reversal potentials in multisynapse conductance-based models, I will do more things:

  • Modifying gif_cond_exp_multisynapse in order to change the number of reversal potentials
  • Adding an example of multisyanpse model (ipython format)
  • Adding a SLI testunit to check consistency for each pair of multisynapse/plain neuron model

Regards,
Hesam

Merge branch 'master' into HesamSetareh_master_261
# Conflicts:
#	nestkernel/nest_names.h

@heplesser heplesser referenced this pull request in hesam-setareh/nest-simulator Sep 29, 2016

Merged

Resolving merge conflicts #4

@heplesser

@hesam-setareh Thank you for all your improvements! I think we are essentially there now. I would just like to ask you to improve the code around allocating y_ a bit, see my comments below. The crucial thing is to avoid C-style casts. And at least in operator=, you need to delete whatever y_ might point to if it is not 0.

With that, the code would be ready to merge from my point.

models/gif_cond_exp_multisynapse.cpp
, sfa_( 0.0 )
, stc_( 0.0 )
, sfa_elems_()
, stc_elems_()
, r_ref_( 0 )
{
+ y_ = ( double* ) malloc( NUMBER_OF_FIXED_STATES_ELEMENTS * sizeof( double ) );

This comment has been minimized.

@heplesser

heplesser Sep 29, 2016

Contributor

You need this allocation several places in this code, so it might be an idea to turn it into a method. That method then could also check if y_ != 0 and call delete first in that case to avoid memory leaks, and update size_neuron_state_. Also, Bjarne Stroustrup advises against using C-style casts, rather use

reinterpret_cast < double* >( ... )
@heplesser

heplesser Sep 29, 2016

Contributor

You need this allocation several places in this code, so it might be an idea to turn it into a method. That method then could also check if y_ != 0 and call delete first in that case to avoid memory leaks, and update size_neuron_state_. Also, Bjarne Stroustrup advises against using C-style casts, rather use

reinterpret_cast < double* >( ... )

This comment has been minimized.

@hesam-setareh

hesam-setareh Oct 3, 2016

Contributor

I have added a new method (neuron_state_memory_allocate) to do this job. A also replaced the old fashion casting with reinterpret_cast.

@hesam-setareh

hesam-setareh Oct 3, 2016

Contributor

I have added a new method (neuron_state_memory_allocate) to do this job. A also replaced the old fashion casting with reinterpret_cast.

models/gif_cond_exp_multisynapse.cpp
@@ -174,7 +177,10 @@ nest::gif_cond_exp_multisynapse::State_::State_( const State_& s )
for ( size_t i = 0; i < stc_elems_.size(); ++i )
stc_elems_[ i ] = s.stc_elems_[ i ];
- y_ = s.y_;
+ size_neuron_state_ = s.size_neuron_state_;
+ y_ = ( double* ) malloc( size_neuron_state_ * sizeof( double ) );

This comment has been minimized.

@heplesser

heplesser Sep 29, 2016

Contributor

See comment on allocation above.

@heplesser

heplesser Sep 29, 2016

Contributor

See comment on allocation above.

This comment has been minimized.

@hesam-setareh

hesam-setareh Oct 3, 2016

Contributor

Fixed.

@hesam-setareh

hesam-setareh Oct 3, 2016

Contributor

Fixed.

models/gif_cond_exp_multisynapse.cpp
@@ -193,7 +197,12 @@ nest::gif_cond_exp_multisynapse::State_&
for ( size_t i = 0; i < stc_elems_.size(); ++i )
stc_elems_[ i ] = s.stc_elems_[ i ];