Inconsistent implementation of refractory period between iaf and aeif models

[Silmathoron]

The implementation of the refractory time is inconsistent between the iaf_* and aeif_* models:

• aeif models "seem to consider" that when a spike occurs, the spike was emitted at the beginning of the step, so the timestep containing the spike is considered in the refractory period.
• iaf and gif models do not see it that way, so there is an additional timestep compared to aeif...

Here is a MWE to show the difference between the iaf_psc_alpha model and the aeif_psc_alpha model for which adaptation was removed and Delta_T was set to zero. refractory_iaf_aeif.zip

I think this problem needs to be solved before #473

[Silmathoron]

Any thought about which behavior should be enforced? I think each implementation is logical but:

• it seems weird to have apparently different refractory periods
• since the aeif model can be equivalent for some sets of parameters, we cannot allow it anyway

I guess since the iaf is kind of the reference model, we should change the aeif to match it's behaviour... I'll try to check on other simulators to see how t_ref is handled for the aeif there.

[heplesser]

I have now come to the conclusion that the implementation in aeif_psc_alpha is incorrect, while iaf_psc_alpha is correct. I assume this carries over to the remaining aeif_* models as well.

The problem is that in aeif_psc_alpha::update(), the refractory counter is counted down before the membrane potential is updated, while for iaf_psc_alpha::update() this happens after the membrane potential update (logically, the way it is implemented it happens instead, but that is fine since the potential is clamped).

The best way to see that the latter makes most sense is to set t_ref to the resolution. Then, the voltage trace for iaf* shows clamped potential for one step, while aeif* shows no clamped potential. If one goes further and sets t_ref to zero, both models show the same behavior (no clamping), as expected. Thus, aeif* with the current implementation shows the same behavior with t_ref == 0 and t_ref == resolution, which is not correct.

Another issue concerns the value of V_m inside the *dynamics() function: if the neuron is refractory, V must be clamped to V_reset in *dynamics(), otherwise w evolution will not be true to the model.

@Silmathoron Could you try to create a PR fixing this for all relevant models, and maybe add a test? The best way to test is perhaps to check that t_ref == 0 and t_ref == resolution lead to different results.

[Silmathoron]

OK, I'll have a look next week

[heplesser]

@Silmathoron I added three new tests at then end of test_iaf_cond_beta_multisynapse.sli in #557. I think those tests generalize to all aeif models (and iaf models, except the third which checks w dynamics). Maybe one could create a separate test that runs these tests for all relevant models?

[Silmathoron]

Ok, I also provided some refractoriness tests in the pending PR correcting the hh_cond_exp_traub model that directly check that the correct duration is implemented, but they can't work until the aeif models are fixed ^^

[heplesser]

I have implemented a reference solution and tests in #557 for aeif_cond_beta_multisynapse. The test mainly check

• that for t_ref == 0 membrane potential is equal to V_reset only at the time step of the spike
• that for t_ref == h, where h is the resolution, membrante potential is equal to V_reset at the time step of the spike and the time step afterwards.

The correct implementation for all models using other integration schemes than exact integrations is as follows.

Right-hand side function (_dynamics):

During the refractory period, voltage must be clamped to V_reset, e.g.,

const double V = node.S_.r_ > 0 ? node.P_.V_reset_ : y[ S::V_M ];

Update loop

Spikes in NEST are always at time steps (except precise models), therefore checks for threshold crossing should be done at the end of a time step. The refractoriness counter is reduced only at the end of each step. This gives the following code structure:

  for ( long lag = from; lag < to; ++lag )
  {
    double t = 0.0;
    while ( t < B_.step_ )
    {
      const int status = gsl_odeiv_evolve_apply( ... );
      // error checks
    }

    if ( S_.r_ > 0 )
    {
      S_.y_[ State_::V_M ] = P_.V_reset_; // clamp it to V_reset
      --S_.r_;
    }
    else if ( S_.y_[ State_::V_M ] >= <threshold> )
    {
      S_.y_[ State_::V_M ] = P_.V_reset_;
      S_.r_ = V_.refractory_counts_;
      // other post-spike changes
      // spike emission
    }

   // input handling
   // log state data
  } 

Note that the check for threshold crossing now is outside the while loop, so that also the reset is moved to the point in time at which the spike is registered. Otherwise, we ended up with refratory periods essentially beginning (and thus ending) inside intervals.

[heplesser]

@DimitriPlotnikov This is most likely relevant for NESTML.

[heplesser]

An alternative implementation for the update loop is the following:

  for ( long lag = from; lag < to; ++lag )
  {
    double t = 0.0;

    while ( t < B_.step_ )
    {
      const int status = gsl_odeiv_evolve_apply( ... ); 
      // error handling

      if ( S_.r_ > 0 )
      {
        S_.y_[ State_::V_M ] = P_.V_reset_;
      }
      else if ( S_.y_[ State_::V_M ] >= <threshold> )
      {
        S_.y_[ State_::V_M ] = P_.V_reset_;
        // add 1 to compensate for count-down right after loop
        S_.r_ = V_.refractory_counts_ + 1;
        // other post-spike actions
        // spike emission
      }
    }   // end while

    if ( S_.r_ > 0 )
    {
      // clamping handled in while loop
      --S_.r_;
    }

    // input handling
    // log state data
  }
}

But the code is slightly more complex. The only difference between the two is that non-clamped variables, e.g., w for aeif models, will "feel" the effect of the spike from the sub-step at which the threshold crossing is detected, while with the first solution the effect starts with the next time step.

I put this up for discussion during the Open Developer VC 28 Nov.

[heplesser]

Design decisin in NEST Open Developer VC 28 Nov: Implement variant 1 with threshold check and refractory countdown outside wihle loop.

[DimitriPlotnikov]

FYI: diese Variant wird momentan in NESTML Generator benutzt.

MfG Dimitri

On 28.11.2016 12:24, Hans Ekkehard Plesser wrote:

Design decisin in NEST Open Developer VC 28 Nov: Implement variant 1 with threshold check and refractory countdown outside wihle loop.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/nest/nest-simulator/issues/547#issuecomment-263247410, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AIV1T95G5NQ_I7Z45CxhwNC0_1UdkAatks5rCrnqgaJpZM4Kz8V8.

-- Dimitri Plotnikov, M.Sc.RWTH Simulation Lab Neuroscience - Bernstein Facility for Simulation and Database Technology Institute for Advanced Simulation Jülich Aachen Research Alliance Forschungszentrum Jülich 52425 Jülich | Germany

---

---

Forschungszentrum Juelich GmbH 52425 Juelich Sitz der Gesellschaft: Juelich Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498 Vorsitzender des Aufsichtsrats: MinDir Dr. Karl Eugen Huthmacher Geschaeftsfuehrung: Prof. Dr.-Ing. Wolfgang Marquardt (Vorsitzender), Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt, Prof. Dr. Sebastian M. Schmidt

---

[heplesser]

@Silmathoron @DimitriPlotnikov @ingablundell I have now explored the issue of "inside" vs "outside" in detail, including a third variant, "inside" thresholding with precise return from refractoriness ("inside_pr"). To this end, I create a separate branch https://github.com/heplesser/nest-simulator/tree/issue-547-experiments with an instrumented aeif_cond_alpha_experimental neuron that allows selecting the thresholding method. Results are in a notebook docs/model_details/Issue-547-Refractoriness.ipynb in that branch.

The most important conclusion is that outside updating is least precise and that precise return from refractoriness has advantages in particular for larger NEST step sizes.

[Silmathoron]

@heplesser thanks a lot for the extensive tests in the notebook! What do you propose? I think the double refractory period is easy to implement and would be much better. However it would be more complex to test for correct implementation. Anyway, I think this clearly rules out the outside detection, so should we reschedule this issue for the next VC?

[Silmathoron]

Thinking back on this issue, I wonder about something: if we chose the inside_pr implementation (and I think we should) how should we deal with the non-gsl models? I ask this question because this would lead us to model that are almost equivalent to the old "gridprecise" models I proposed in the original PR on the AEIF models. However, at that time I had deliberately created two different models because the behaviour of this AEIF implementation will no longer be equivalent to the IAF implementation when Delta_T goes to zero. Paradoxically, we would have a more precise IAF model using the non-adapting version of the AEIF...

[heplesser]

I propose that we settle for the plain "inside" method for now, correcting the count-down of refractory steps where necessary. This works well and keeps us close to the implementation we have had in most cond-based models so far. "inside_pr" has some issues as well (but beyond the scope of this issue) and requires more complex changes to models. I will send a mail to the list about these issues shortly, as I would like to implement changes next week so we get 2.12 out of the door by Christmas.

Consistency between exact integration with fixed time step and adaptive stepsize solvers is an issue. But we may need to accept some limits to consistency also for the sake of efficiency. The precise models provide the most exact implementations of linear models possible, but at a cost. And for most purposes, it seems, standard I&F with refractory periods locked to the grid work fine. It is really just the exponentially exploding current in AEIF that makes that model challenging.

[Silmathoron]

Fine with me, I'll try to make the PR quickly