V-bounded implementation for AEIF + fixed Delta_T=0 #474

Merged
merged 18 commits into from Oct 5, 2016

Conversation

Projects
None yet
4 participants
@Silmathoron
Contributor

Silmathoron commented Sep 8, 2016

This PR proposes a corrected implementation of the existing AEIF models that leads to a solution which is closer to the reference LSODAR implementation (see discussion in PR #229 and provided notebook).
It also includes a fix for issue #114 by setting a separate implementation for the case of Delta_T == 0.

Remaining question: currently the argument of the exponential is not bounded anymore (except for the multisynapse implementation because the numerical instability test is different), which means that numerical overflow will occur when Delta_T becomes small enough.
Should we consider this an unavoidable issue of the model and consider that the domain where it occurs are not biologically relevant anyway? Or should the exponential still be bounded?
If it should be bounded, I propose to use numeric_limits and MAX_DOUBLE instead of an ad hoc value, but it should be noted that in those cases, the model will deviate from the real solution anyway...

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 25, 2016

Contributor

@Silmathoron and @heplesser I finally find some time to read your messages and discussion in #229 It's a very interesting and valuable work!
Concerning this pull request, the first thing that I have to say is that there are several issues with the aeif_*_multisynapse models that are not related to your changes, see https://github.com/nest/nest-simulator/projects/1 , and discussion in #439 . I dealt with those issues, with the help of @heplesser and @DimitriPlotnikov, in the latter pull request for the aeif_cond_beta_multisynapse model. My suggestion is that you send a separate pull request for the single inh/exc synapse models, so that your pull request can be merged quickly, and meanwhile we can discuss the issues for multisynapse models.

Contributor

golosio commented Sep 25, 2016

@Silmathoron and @heplesser I finally find some time to read your messages and discussion in #229 It's a very interesting and valuable work!
Concerning this pull request, the first thing that I have to say is that there are several issues with the aeif_*_multisynapse models that are not related to your changes, see https://github.com/nest/nest-simulator/projects/1 , and discussion in #439 . I dealt with those issues, with the help of @heplesser and @DimitriPlotnikov, in the latter pull request for the aeif_cond_beta_multisynapse model. My suggestion is that you send a separate pull request for the single inh/exc synapse models, so that your pull request can be merged quickly, and meanwhile we can discuss the issues for multisynapse models.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 26, 2016

Contributor

Hi @golosio, thanks for your reply.
Just to be sure I understand: you suggest that I modify only the aeif_cond_alpha, aeif_cond_exp, aeif_cond_alpha_RK5, keeping the aeif_cond_alpha_multisynapse unchanged, is that it?

I understand you point, but since this PR concerns the numerical implementation of the ODE and not the model design implementation, it does not really interfere with your PR, right?
I mean we are modifying different parts of the code, so if I include the modification in aeif_cond_alpha_multisynapse, it will save someone (probably you, in fact ^^) the trouble of modifying another part of the code in the "multisynapse design PR", don't you think?

Contributor

Silmathoron commented Sep 26, 2016

Hi @golosio, thanks for your reply.
Just to be sure I understand: you suggest that I modify only the aeif_cond_alpha, aeif_cond_exp, aeif_cond_alpha_RK5, keeping the aeif_cond_alpha_multisynapse unchanged, is that it?

I understand you point, but since this PR concerns the numerical implementation of the ODE and not the model design implementation, it does not really interfere with your PR, right?
I mean we are modifying different parts of the code, so if I include the modification in aeif_cond_alpha_multisynapse, it will save someone (probably you, in fact ^^) the trouble of modifying another part of the code in the "multisynapse design PR", don't you think?

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

@Silmathoron What I mean is that aeif_____multisynapse need an extensive review, which have already been done in aeif_cond_beta_multisynapse. I am not sure that your PR can be merged without this review, even though your changes are related to other parts of the code. In my opinion the easiest way to proceed is that you present a first PR only for aeif_cond_alpha, aeif_cond_exp, aeif_cond_alpha_RK5, which could be merged quickly. Then if you are available you could apply your changes to aeif_cond_beta_multisynapse (which was just accepted for merging) with a new PR (if you can't I can try to do it) and from this it would be straightforward to rewrite aeif_cond_alpha_multisynapse.

Contributor

golosio commented Sep 26, 2016

@Silmathoron What I mean is that aeif_____multisynapse need an extensive review, which have already been done in aeif_cond_beta_multisynapse. I am not sure that your PR can be merged without this review, even though your changes are related to other parts of the code. In my opinion the easiest way to proceed is that you present a first PR only for aeif_cond_alpha, aeif_cond_exp, aeif_cond_alpha_RK5, which could be merged quickly. Then if you are available you could apply your changes to aeif_cond_beta_multisynapse (which was just accepted for merging) with a new PR (if you can't I can try to do it) and from this it would be straightforward to rewrite aeif_cond_alpha_multisynapse.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 26, 2016

Contributor

Ok, sounds good to me, the modification is really easy to do so it should not take long to modify your model afterwards.

Contributor

Silmathoron commented Sep 26, 2016

Ok, sounds good to me, the modification is really easy to do so it should not take long to modify your model afterwards.

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

@Silmathoron In aeif_cond_alpha_RK5.cpp line 438

    ( this->*( V_.model_dynamics ) )( S_.yin, S_.k6 );
    // 5th order
    for ( int i = 0; i < S_.STATE_VEC_SIZE; ++i )
      S_.ynew[ i ] = S_.y_[ i ]
        + h * ( 35.0 / 384.0 * S_.k1[ i ] + 500.0 / 1113.0 * S_.k3[ i ]
                + 125.0 / 192.0 * S_.k4[ i ]
                - 2187.0 / 6784.0 * S_.k5[ i ]
                + 11.0 / 84.0 * S_.k6[ i ] );  
    ( this->*( V_.model_dynamics ) )( S_.yin, S_.k7 );

I've worked on RK5 in other projects and I think there is a mistake in the last line. The first argument should be S_.ynew, not S_.yin Indeed, with this code k7 will always be the same as k6, which does not make much sense, because yIn is unchanged (see the code above). As k7 is only necessary for evaluating the error, it probably does not have dramatic effects on the solution. I've already noticed this problem in aeif_cond_alpha_multisynapse. I know it has nothing to do with your changes, however since we are here it is better to clarify this issue.

Contributor

golosio commented Sep 26, 2016

@Silmathoron In aeif_cond_alpha_RK5.cpp line 438

    ( this->*( V_.model_dynamics ) )( S_.yin, S_.k6 );
    // 5th order
    for ( int i = 0; i < S_.STATE_VEC_SIZE; ++i )
      S_.ynew[ i ] = S_.y_[ i ]
        + h * ( 35.0 / 384.0 * S_.k1[ i ] + 500.0 / 1113.0 * S_.k3[ i ]
                + 125.0 / 192.0 * S_.k4[ i ]
                - 2187.0 / 6784.0 * S_.k5[ i ]
                + 11.0 / 84.0 * S_.k6[ i ] );  
    ( this->*( V_.model_dynamics ) )( S_.yin, S_.k7 );

I've worked on RK5 in other projects and I think there is a mistake in the last line. The first argument should be S_.ynew, not S_.yin Indeed, with this code k7 will always be the same as k6, which does not make much sense, because yIn is unchanged (see the code above). As k7 is only necessary for evaluating the error, it probably does not have dramatic effects on the solution. I've already noticed this problem in aeif_cond_alpha_multisynapse. I know it has nothing to do with your changes, however since we are here it is better to clarify this issue.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 26, 2016

Contributor

Ok, I'll have a look, thanks.
In the meantime, if you want to comment on code, please use the GitHub tool:

  • browse the file
  • find the line
  • click on the + that appears on the right, next to the line number
  • write your comment
Contributor

Silmathoron commented Sep 26, 2016

Ok, I'll have a look, thanks.
In the meantime, if you want to comment on code, please use the GitHub tool:

  • browse the file
  • find the line
  • click on the + that appears on the right, next to the line number
  • write your comment
@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

Concerning your "remaining question": if you bind V to V_peak, for Delta_T!=0 the exponential argument maximum value should be (V_peak - V_th ) / Delta_T , am I right? Maybe this value could be checked in the Parameters_::set and give an error or warning, e.g. suggesting either to use Delta_T=0 or to use larger values?

Contributor

golosio commented Sep 26, 2016

Concerning your "remaining question": if you bind V to V_peak, for Delta_T!=0 the exponential argument maximum value should be (V_peak - V_th ) / Delta_T , am I right? Maybe this value could be checked in the Parameters_::set and give an error or warning, e.g. suggesting either to use Delta_T=0 or to use larger values?

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

@Silmathoron thank you for the explanation on how to use the GitHub tool, I'll use it in the future.
However, in this case I've already done checking all your changes and apart from the previous comments everythig looks good to me.

Contributor

golosio commented Sep 26, 2016

@Silmathoron thank you for the explanation on how to use the GitHub tool, I'll use it in the future.
However, in this case I've already done checking all your changes and apart from the previous comments everythig looks good to me.

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

@helplesser I don't know if I'm allowed to act as a reviewer for this PR, however I've read all your previous discussion in #229 and I revised the code. After clarification of the few points that I raised above, I would suggest to merge this PR.

Contributor

golosio commented Sep 26, 2016

@helplesser I don't know if I'm allowed to act as a reviewer for this PR, however I've read all your previous discussion in #229 and I revised the code. After clarification of the few points that I raised above, I would suggest to merge this PR.

@abigailm

This comment has been minimized.

Show comment
Hide comment
@abigailm

abigailm Sep 26, 2016

Contributor

@golosio - you most certainly are welcome to act as a reviewer on this PR or, indeed, on any other!

Contributor

abigailm commented Sep 26, 2016

@golosio - you most certainly are welcome to act as a reviewer on this PR or, indeed, on any other!

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Sep 26, 2016

Contributor

Thank you @abigailm

Contributor

golosio commented Sep 26, 2016

Thank you @abigailm

@heplesser

@Silmathoron Good work, almost ready to merge.

As I commented in one file, also single-line blocks should be wrapped in braces. You don't need to add braces everywhere in old code, but you should at least not remove braces, and have them in new code. Background is the OpenSSH double-goto bug a few years ago.

Concerning DeltaT==0, I agree with @golosio: The exponential term is bounded since V is bounded by V_peak. For the default parameters from the Brette&Gerstner paper, the exponential reaches almost 1e11. I suggest that we make an arbitrary choice: if '(V_peak - V_th)/Delta_T > 50', corresponding to e^50 = 5e21, we switch to the Delta_T==0 case. We could consider other limits, but with double ranging to 1e308, I guess we want to stay below 1e100, corresponding to a limit of 230.

Concerning the modifications of the multisynapse input, I would suggest that we keep this PR as it is, and that @Silmathoron then very soon creates another PR fixing the way input is handled according to @golosio's work. In this way we reduce the number of dependencies between PRs.

Eventually, we should drop the RK5 variant entirely, but that is a different issue (#494).

@@ -0,0 +1,711 @@
+{

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

@Silmathoron Could you add a %matplotlib inline to the notebook so that plots are included? That way, everyone will be able to read them without having to install SUNDIALS etc. And could you add some "zoomed in" variants of the plots, so that one sees in more detail where they agree or differ?

@heplesser

heplesser Sep 27, 2016

Contributor

@Silmathoron Could you add a %matplotlib inline to the notebook so that plots are included? That way, everyone will be able to read them without having to install SUNDIALS etc. And could you add some "zoomed in" variants of the plots, so that one sees in more detail where they agree or differ?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
@@ -219,45 +219,33 @@ aeif_cond_alpha_multisynapse::Parameters_::set( const DictionaryDatum& d )
if ( updateValue< double >( d, names::MAXERR, tmp ) )
{
if ( not( tmp > 0.0 ) )
- {

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

NEST style guidelines require that also one-line blocks are included in braces for safety's sake.

@heplesser

heplesser Sep 27, 2016

Contributor

NEST style guidelines require that also one-line blocks are included in braces for safety's sake.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Ok, done... I'm surprised, though, what was the problem that occurred exactly?

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Ok, done... I'm surprised, though, what was the problem that occurred exactly?

models/aeif_cond_alpha_multisynapse.cpp
MAXERR = tmp;
}
if ( updateValue< double >( d, names::HMIN, tmp ) )
{
if ( not( tmp > 0.0 ) )
- {

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

See above, please re-add braces.

@heplesser

heplesser Sep 27, 2016

Contributor

See above, please re-add braces.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
HMIN = tmp;
}
- if ( V_peak_ <= V_th )
- {
+ if ( Delta_T != 0. && V_peak_ <= V_th )

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Please add braces around if and else block

@heplesser

heplesser Sep 27, 2016

Contributor

Please add braces around if and else block

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
if ( V_reset_ >= V_peak_ )
- {

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces.

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
if ( t_ref_ < 0 )
- {

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces.

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
if ( tau_w <= 0 )
- {

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces

@heplesser

heplesser Sep 27, 2016

Contributor

Keep braces

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

models/aeif_cond_alpha_multisynapse.cpp
- assert( V_.RefractoryCounts_
+
+ // set the right function for the dynamics
+ if ( P_.Delta_T != 0. )

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Add braces!

@heplesser

heplesser Sep 27, 2016

Contributor

Add braces!

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Done

@@ -524,25 +521,25 @@ aeif_cond_alpha_multisynapse::update( Time const& origin,
t_return = t + h; // update t
// k1 = f(told, y)
- aeif_cond_alpha_multisynapse_dynamics( S_.y_, S_.k1 );
+ ( this->*( V_.model_dynamics ) )( S_.y_, S_.k1 );

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

Is the this-> pointer here strictly necessary?

@heplesser

heplesser Sep 27, 2016

Contributor

Is the this-> pointer here strictly necessary?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

I think so, see SO post

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

I think so, see SO post

models/aeif_cond_alpha_RK5.cpp
@@ -427,7 +435,7 @@ void nest::aeif_cond_alpha_RK5::update( Time const& origin,
+ 125.0 / 192.0 * S_.k4[ i ]
- 2187.0 / 6784.0 * S_.k5[ i ]
+ 11.0 / 84.0 * S_.k6[ i ] );
- aeif_cond_alpha_RK5_dynamics( S_.yin, S_.k7 );
+ ( this->*( V_.model_dynamics ) )( S_.yin, S_.k7 );

This comment has been minimized.

@heplesser

heplesser Sep 27, 2016

Contributor

I agree with @golosio that it should be S_.ynew instead of S_.yin here.

@heplesser

heplesser Sep 27, 2016

Contributor

I agree with @golosio that it should be S_.ynew instead of S_.yin here.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Ok, I changed it.

@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Ok, I changed it.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 27, 2016

Contributor

Ok, I incorporated the changes and started to use the pep8 pip package, which finally allowed me to get the better of those awful pep8 errors!

Regarding the code, there are still remarks and questions:

  • I chose to raise an error rather than switch the Delta_T behaviour because I think it's better to yell at the user that something is wrong, does it seem acceptable to you?
  • I stored the LSODAR data in a .dat file; is there any easy way to access it using nest.tests? (this is what Travis is complaining about right now)
  • the RK5-based implementation do not pass the test
Contributor

Silmathoron commented Sep 27, 2016

Ok, I incorporated the changes and started to use the pep8 pip package, which finally allowed me to get the better of those awful pep8 errors!

Regarding the code, there are still remarks and questions:

  • I chose to raise an error rather than switch the Delta_T behaviour because I think it's better to yell at the user that something is wrong, does it seem acceptable to you?
  • I stored the LSODAR data in a .dat file; is there any easy way to access it using nest.tests? (this is what Travis is complaining about right now)
  • the RK5-based implementation do not pass the test
@heplesser

I agree with raising an exception for "bad" Delta_T.

I think your use of the reference data file is the best option for now. The only alternative would be to paste 2000 lines of data right into the Python code. Maybe add some comment lines on top of the file explaining where it comes from and why it is there?

The problem with RK5 is probably because that still uses the old bounding scheme. If you converted it to the new scheme, I suppose it would get close enough to pass the test.

models/aeif_cond_alpha.cpp
+ 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." );

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

You can write long strings like this in C++:

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." );
@heplesser

heplesser Sep 28, 2016

Contributor

You can write long strings like this in C++:

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." );

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

corrected, thx

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

corrected, thx

models/aeif_cond_alpha.cpp
@@ -280,6 +280,17 @@ nest::aeif_cond_alpha::Parameters_::set( const DictionaryDatum& d )
if ( Delta_T != 0. && V_peak_ <= V_th )
throw BadProperty( "V_peak must be larger than threshold." );
+ else if ( Delta_T != 0. )

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

To be on the defensive side, I would rather test Delta_T > 0. Also, in the test above, I would check for V_peak_ <= V_th even if Delta_T == 0, since V_peak_ < V_th simply makes no sense (why does V_peak_ have a trailing _?).

So I would go for something like:

if ( V_peak_ <= V_th )
{
      throw BadProperty( "V_peak must be larger than threshold." );
}

if ( Delta_T < 0 )
{
  throw BadProperty( "..." );
}
else if ( Delta_T > 0 )
{
}
else
{
  // Delta_T == 0
}
@heplesser

heplesser Sep 28, 2016

Contributor

To be on the defensive side, I would rather test Delta_T > 0. Also, in the test above, I would check for V_peak_ <= V_th even if Delta_T == 0, since V_peak_ < V_th simply makes no sense (why does V_peak_ have a trailing _?).

So I would go for something like:

if ( V_peak_ <= V_th )
{
      throw BadProperty( "V_peak must be larger than threshold." );
}

if ( Delta_T < 0 )
{
  throw BadProperty( "..." );
}
else if ( Delta_T > 0 )
{
}
else
{
  // Delta_T == 0
}

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

This is much better, indeed, I'll change it!
No idea about V_peak_; at first I thought it was because users are not supposed to access the variables directly but since only a few of them are like that, it doesn't make sense... do you want me to remove it? What about t_ref_ 'n co. ?

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

This is much better, indeed, I'll change it!
No idea about V_peak_; at first I thought it was because users are not supposed to access the variables directly but since only a few of them are like that, it doesn't make sense... do you want me to remove it? What about t_ref_ 'n co. ?

models/aeif_cond_alpha.cpp
+ "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." );
+ }
else if ( Delta_T == 0. )
updateValue< double >( d, names::V_peak, V_th ); // expected behaviour

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

Why do you use name::V_peak to update V_th?

@heplesser

heplesser Sep 28, 2016

Contributor

Why do you use name::V_peak to update V_th?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Because the spike occurs at V_th for Delta_T = 0., just like in a regular iaf neuron. I'll add a longer comment to explain it in the code.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Because the spike occurs at V_th for Delta_T = 0., just like in a regular iaf neuron. I'll add a longer comment to explain it in the code.

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

Yes, but the coding here is unfortunate. On lines 258-259, we extract V_peak and V_th using updateValue(). If we now need to set V_peak to V_th, that should be done by direct assignment, not by calling updateValue() again. And strictly speaking, we should not change V_th (since that is a user-defined parameter). Rather, we should introduce an additional variable in V_ that represents the V_th to use, and that is set in calibrate():

if ( P_.Delta_T > 0 )
{
  V_.V_th = P_.V_th;
}
else
{
  V_.V_th = P_.V_peak;
}

(I hope I got the logic the right way around?)

Then the user could even switch back and forth between Delta_T=0 and > 0 and all would work as expected.

@heplesser

heplesser Sep 28, 2016

Contributor

Yes, but the coding here is unfortunate. On lines 258-259, we extract V_peak and V_th using updateValue(). If we now need to set V_peak to V_th, that should be done by direct assignment, not by calling updateValue() again. And strictly speaking, we should not change V_th (since that is a user-defined parameter). Rather, we should introduce an additional variable in V_ that represents the V_th to use, and that is set in calibrate():

if ( P_.Delta_T > 0 )
{
  V_.V_th = P_.V_th;
}
else
{
  V_.V_th = P_.V_peak;
}

(I hope I got the logic the right way around?)

Then the user could even switch back and forth between Delta_T=0 and > 0 and all would work as expected.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

This looks very handy, I'll implement it! (but it's V_peak that changes ;) )

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

This looks very handy, I'll implement it! (but it's V_peak that changes ;) )

throw BadProperty( "All time constants must be strictly positive." );
+ }

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

See discussion above.

@heplesser

heplesser Sep 28, 2016

Contributor

See discussion above.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

You're talking about the tests for Delta_T, V_peak and V_th, right?

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

You're talking about the tests for Delta_T, V_peak and V_th, right?

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

@Silmathoron Now I managed to confuse myself. Was this about adding braces? Anyways, the code looks fine here.

@heplesser

heplesser Sep 28, 2016

Contributor

@Silmathoron Now I managed to confuse myself. Was this about adding braces? Anyways, the code looks fine here.

models/aeif_cond_alpha.cpp
+ {
+ // check for possible numerical overflow with the exponential divergence at
+ // spike time
+ double max_exp_arg = std::log( std::numeric_limits< double >::max() );

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

This should be const. I think this limit is too high: we now prohibit only value combinations if the resulting exp() would exceed the capacity of double exactly at the threshold. But we need to do some computations with the result of the exponential, so we need to keep a good margin to the double limit. I would suggest at least max() / 1e20.

@heplesser

heplesser Sep 28, 2016

Contributor

This should be const. I think this limit is too high: we now prohibit only value combinations if the resulting exp() would exceed the capacity of double exactly at the threshold. But we need to do some computations with the result of the exponential, so we need to keep a good margin to the double limit. I would suggest at least max() / 1e20.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

done

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

done

models/aeif_cond_alpha.cpp
+ // spike time
+ double max_exp_arg = std::log( std::numeric_limits< double >::max() );
+ if ( ( V_peak_ - V_th ) / Delta_T >= max_exp_arg )
+ throw BadProperty(

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

Throwing an exception here is fine.

@heplesser

heplesser Sep 28, 2016

Contributor

Throwing an exception here is fine.

pynest/nest/tests/test_aeif_lsodar.py
+ return idx
+
+
+def _interpolate_lsodar(nest_times, lsodar):

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

Just wondering why you are not using the standard interpolation tools provided by NumPy/SciPy?

@heplesser

heplesser Sep 28, 2016

Contributor

Just wondering why you are not using the standard interpolation tools provided by NumPy/SciPy?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Because when the problem occurred I did not think about it and I did it on the fly, which is pretty stupid... I'll have a look at them.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Because when the problem occurred I did not think about it and I did it on the fly, which is pretty stupid... I'll have a look at them.

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

That would be great. Using standard tools usually makes code more maintainable.

@heplesser

heplesser Sep 28, 2016

Contributor

That would be great. Using standard tools usually makes code more maintainable.

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Yes, it also takes only 10 s instead of 10 min and 2 lines of codes instead of 50...

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Yes, it also takes only 10 s instead of 10 min and 2 lines of codes instead of 50...

pynest/nest/tests/test_aeif_lsodar.py
+Details:
+ The models are compared and we assess that the difference is smaller than a
+ given tolerance.
+"""

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

Could you bring this closer to the unittest setup as in the other pynest tests?

@heplesser

heplesser Sep 28, 2016

Contributor

Could you bring this closer to the unittest setup as in the other pynest tests?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Yes, I'm currently doing it; however I don't know which path I need to provide to load the .dat file on Travis: does NEST provide a way to get the path to the source directory?

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Yes, I'm currently doing it; however I don't know which path I need to provide to load the .dat file on Travis: does NEST provide a way to get the path to the source directory?

This comment has been minimized.

@heplesser

heplesser Sep 28, 2016

Contributor

@lekshmideepu Could you help @Silmathoron with this?

@heplesser

heplesser Sep 28, 2016

Contributor

@lekshmideepu Could you help @Silmathoron with this?

This comment has been minimized.

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

It's okay, I found the way by modifying the setup.py

@Silmathoron

Silmathoron Sep 28, 2016

Contributor

It's okay, I found the way by modifying the setup.py

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Sep 28, 2016

Contributor

Ok, all modifications have been done and I added a test to compare with iaf_cond_* when Delta_T = 0., a = 0., b = 0. and w = 0..
I managed to get the RK5 models to pass the test but I had to increase their specific tolerance and I need to keep the exponential bounded or it leads to either numeric instability or parasitic jumps in w at spike times because of the divergence of V.

Now I just need to find out how to load the .dat from nest.tests before I push.

EDIT: ok, I think it should be ok, I modified the setup.py to tell python to include the .dat
EDIT2: victory!

Contributor

Silmathoron commented Sep 28, 2016

Ok, all modifications have been done and I added a test to compare with iaf_cond_* when Delta_T = 0., a = 0., b = 0. and w = 0..
I managed to get the RK5 models to pass the test but I had to increase their specific tolerance and I need to keep the exponential bounded or it leads to either numeric instability or parasitic jumps in w at spike times because of the divergence of V.

Now I just need to find out how to load the .dat from nest.tests before I push.

EDIT: ok, I think it should be ok, I modified the setup.py to tell python to include the .dat
EDIT2: victory!

@heplesser

@Silmathoron Almost perfect now except for unfortunate mix of natural and decadic log, see inline comments.

models/aeif_cond_alpha.cpp
- double max_exp_arg = std::log( std::numeric_limits< double >::max() );
+ // spike time, keep a 1e20 margin for the subsequent calculations
+ const double max_exp_arg =
+ std::log( std::numeric_limits< double >::max() ) - 20.;

This comment has been minimized.

@heplesser

heplesser Oct 5, 2016

Contributor

@Silmathoron Here you mix natural and decadic log. I would suggest

std::log( std::numeric_limits< double >::max() / 1e20 )

for clarity and correctness.

@heplesser

heplesser Oct 5, 2016

Contributor

@Silmathoron Here you mix natural and decadic log. I would suggest

std::log( std::numeric_limits< double >::max() / 1e20 )

for clarity and correctness.

models/aeif_cond_alpha_RK5.cpp
+ // 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() ) - 20.;

This comment has been minimized.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

models/aeif_cond_alpha_multisynapse.cpp
+ // 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() ) - 20.;

This comment has been minimized.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

models/aeif_cond_exp.cpp
- double max_exp_arg = std::log( std::numeric_limits< double >::max() );
+ // spike time, keep a 1e20 margin for the subsequent calculations
+ const double max_exp_arg =
+ std::log( std::numeric_limits< double >::max() ) - 20.;

This comment has been minimized.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

@heplesser

heplesser Oct 5, 2016

Contributor

See comment on aeif_cond_alpha.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

Fixed, sorry about that (propagate upwards)

@Silmathoron

Silmathoron Oct 5, 2016

Contributor

Fixed, sorry about that (propagate upwards)

@heplesser heplesser merged commit 446d340 into nest:master Oct 5, 2016

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Oct 5, 2016

Contributor

@Silmathoron Congratulations for all the good work that you've done!

Contributor

golosio commented Oct 5, 2016

@Silmathoron Congratulations for all the good work that you've done!

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Oct 5, 2016

Contributor

Thanks at @heplesser and @golosio for the impressive review and constructive ideas ;)

Contributor

Silmathoron commented Oct 5, 2016

Thanks at @heplesser and @golosio for the impressive review and constructive ideas ;)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment