Refractory implementation aeif models (fixes #547) #575

Merged
merged 8 commits into from Dec 12, 2016

Conversation

Projects
None yet
3 participants
@Silmathoron
Contributor

Silmathoron commented Dec 3, 2016

Following #547, this corrects the refractory implementation for:

  • aeif_cond_alpha
  • aeif_cond_exp
  • aeif _psc_exp
  • aeif_psc_alpha

This PR also removes the DT0 dynamics function in the aeif_cond_alpha/exp models.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 3, 2016

Contributor

@heplesser this new implementation seems a lot less precise than the previous one with regard to LSODAR... (hence the failing test)
I'm investigating
EDIT: problem comes from the reset to V_reset before we're out of the first while (when S_.r_ == V_.refractory_counts_ + 1)

Contributor

Silmathoron commented Dec 3, 2016

@heplesser this new implementation seems a lot less precise than the previous one with regard to LSODAR... (hence the failing test)
I'm investigating
EDIT: problem comes from the reset to V_reset before we're out of the first while (when S_.r_ == V_.refractory_counts_ + 1)

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 3, 2016

Contributor

@heplesser comparing the relative difference between nest and the LSODAR solution shows that the smallest divergence is obtained without applying V = V_reset and f[ S::V_M ] = 0. in the dynamics function.
I think this should be applied to #557 also.

Contributor

Silmathoron commented Dec 3, 2016

@heplesser comparing the relative difference between nest and the LSODAR solution shows that the smallest divergence is obtained without applying V = V_reset and f[ S::V_M ] = 0. in the dynamics function.
I think this should be applied to #557 also.

@Silmathoron Silmathoron changed the title from Refractory implementation aeif models to Refractory implementation aeif models (fixes #547) Dec 3, 2016

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 5, 2016

Contributor

@Silmathoron Thanks for adding the tests, but I think you may have gotten it the wrong way around: I believe the LSODAR reference solution is incorrect, since it does not know about refractoriness.
I tested aeif_cond_beta_multisynapse against LSODAR with refractoriness in my heplesser/Issue-547-Refractoriness branch and there it worked perfectly. Here are the essential bits of my LSODAR code, you will find the rest in the doc/model_details notebook in my branch:

    def rhs(self, t, y, sw):
        """
        This function returns the dV/dt, dw/dt
        """
        
        is_refractory = t < self.t_end_refractory

        V = self.p['V_reset'] if is_refractory else y[0]
        w = y[1]

        Ispike = ( 0. if self.p['Delta_T'] == 0 else
                   self.p['g_L'] * self.p['Delta_T'] * 
                   np.exp( (V-self.p['V_th']) / self.p['Delta_T'] ) )

        dotV = ( 0. if is_refractory else
                 ( -self.p['g_L'] * (V-self.p['E_L']) + Ispike + self.p['I_e'] - w ) / 
                   self.p['C_m'] )
        dotW = ( self.p['a'] * (V-self.p['E_L']) - w ) / self.p['tau_w']

        return np.array([dotV, dotW])

    def handle_event(self, solver, event_info):
        """
        Event handling. This functions is called when Assimulo finds an event as
        specified by the event functions.
        """
        # only look at the state events information.
        if event_info[0][0] > 0:
            solver.sw[0] = True
            solver.y[0] = self.p['V_reset']
            solver.y[1] += self.p['b']
            self.t_end_refractory = solver.t + self.p['t_ref']
        else:
            solver.sw[0] = False

One could (should?) probably go yet a step further and code return from refractoriness as an event as well.

Contributor

heplesser commented Dec 5, 2016

@Silmathoron Thanks for adding the tests, but I think you may have gotten it the wrong way around: I believe the LSODAR reference solution is incorrect, since it does not know about refractoriness.
I tested aeif_cond_beta_multisynapse against LSODAR with refractoriness in my heplesser/Issue-547-Refractoriness branch and there it worked perfectly. Here are the essential bits of my LSODAR code, you will find the rest in the doc/model_details notebook in my branch:

    def rhs(self, t, y, sw):
        """
        This function returns the dV/dt, dw/dt
        """
        
        is_refractory = t < self.t_end_refractory

        V = self.p['V_reset'] if is_refractory else y[0]
        w = y[1]

        Ispike = ( 0. if self.p['Delta_T'] == 0 else
                   self.p['g_L'] * self.p['Delta_T'] * 
                   np.exp( (V-self.p['V_th']) / self.p['Delta_T'] ) )

        dotV = ( 0. if is_refractory else
                 ( -self.p['g_L'] * (V-self.p['E_L']) + Ispike + self.p['I_e'] - w ) / 
                   self.p['C_m'] )
        dotW = ( self.p['a'] * (V-self.p['E_L']) - w ) / self.p['tau_w']

        return np.array([dotV, dotW])

    def handle_event(self, solver, event_info):
        """
        Event handling. This functions is called when Assimulo finds an event as
        specified by the event functions.
        """
        # only look at the state events information.
        if event_info[0][0] > 0:
            solver.sw[0] = True
            solver.y[0] = self.p['V_reset']
            solver.y[1] += self.p['b']
            self.t_end_refractory = solver.t + self.p['t_ref']
        else:
            solver.sw[0] = False

One could (should?) probably go yet a step further and code return from refractoriness as an event as well.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 5, 2016

Contributor

@heplesser Ok, thanks for the answer!
I'm may indeed have gotten a bit too far in removing all the changes in the dynamics function; however:

  • regardless of the refractory implementation, the LSODAR data was exact when t_ref = 0, therefore we cannot allow the implementation to diverge from it in this limit (especially since t_ref is seldom used in the aeif models)
  • this discrepancy between NEST and LSODAR came from the divergence (as usual) during the "spike timestep" and was therefore not completely "refractory-related".

Therefore I propose the following correction:

  • make sure that the spiking-step remains unchanged by testing if ( S_.r_ > 0 && S_.r_ <= V_.refractory_counts_ + 1) in the update function
  • use the same condition for is_refractory in the dynamics function.

As I said before, we need to make sure that the models still pass the tests in the t_ref = 0 limits, so I think this should do the trick, while correcting the refractory dynamics.

Contributor

Silmathoron commented Dec 5, 2016

@heplesser Ok, thanks for the answer!
I'm may indeed have gotten a bit too far in removing all the changes in the dynamics function; however:

  • regardless of the refractory implementation, the LSODAR data was exact when t_ref = 0, therefore we cannot allow the implementation to diverge from it in this limit (especially since t_ref is seldom used in the aeif models)
  • this discrepancy between NEST and LSODAR came from the divergence (as usual) during the "spike timestep" and was therefore not completely "refractory-related".

Therefore I propose the following correction:

  • make sure that the spiking-step remains unchanged by testing if ( S_.r_ > 0 && S_.r_ <= V_.refractory_counts_ + 1) in the update function
  • use the same condition for is_refractory in the dynamics function.

As I said before, we need to make sure that the models still pass the tests in the t_ref = 0 limits, so I think this should do the trick, while correcting the refractory dynamics.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 5, 2016

Contributor

@Silmathoron You are right, we need to handle t_ref==0 correctly and I also see the problem with the "inside" variant proposed in #547. Your code is much closer to "inside_pr". It took me a long time to understand how it worked, though, so I will shortly send you a PR with a slight reformulation which allows us to keep the simpler S_.r_ > 0 tests without checking the upper limit.

Contributor

heplesser commented Dec 5, 2016

@Silmathoron You are right, we need to handle t_ref==0 correctly and I also see the problem with the "inside" variant proposed in #547. Your code is much closer to "inside_pr". It took me a long time to understand how it worked, though, so I will shortly send you a PR with a slight reformulation which allows us to keep the simpler S_.r_ > 0 tests without checking the upper limit.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 5, 2016

Contributor

Assigning this to 2.12 milestone since it is part of the aeif_* cleanup process.

Contributor

heplesser commented Dec 5, 2016

Assigning this to 2.12 milestone since it is part of the aeif_* cleanup process.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 6, 2016

Contributor

@heplesser I don't know whether the labels are really important or for posterity, but shouldn't this be reported as Bug?

I'll switch the implementation to the proposal in #557 shortly.

Contributor

Silmathoron commented Dec 6, 2016

@heplesser I don't know whether the labels are really important or for posterity, but shouldn't this be reported as Bug?

I'll switch the implementation to the proposal in #557 shortly.

@heplesser heplesser added T: Bug and removed T: Maintenance labels Dec 6, 2016

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 8, 2016

Contributor

@heplesser I pushed the commit 2 days ago but just realized they are not appearing on GitHub...
git tells me everything is up to date, though, and git ls-remote indeed shows the right commit number for t_ref_aeif...
Did this ever occur to you?

EDIT: just in case, git log --decorate gives

commit a2cd3a606aa1cfedf79a5d4a6b35de58c783aef4 (HEAD -> t_ref_aeif, origin/t_ref_aeif)
Merge: 40f05ca ec8d479
Author: Silmathoron <Silmathoron@users.noreply.github.com>
Date:   Tue Dec 6 21:41:42 2016 +0100

merge and recorrected

commit 40f05ca2ce95226ab17af02765182a141018593e
Author: Silmathoron <Silmathoron@users.noreply.github.com>
Date:   Tue Dec 6 17:05:07 2016 +0100

changed implementation to match #557
Contributor

Silmathoron commented Dec 8, 2016

@heplesser I pushed the commit 2 days ago but just realized they are not appearing on GitHub...
git tells me everything is up to date, though, and git ls-remote indeed shows the right commit number for t_ref_aeif...
Did this ever occur to you?

EDIT: just in case, git log --decorate gives

commit a2cd3a606aa1cfedf79a5d4a6b35de58c783aef4 (HEAD -> t_ref_aeif, origin/t_ref_aeif)
Merge: 40f05ca ec8d479
Author: Silmathoron <Silmathoron@users.noreply.github.com>
Date:   Tue Dec 6 21:41:42 2016 +0100

merge and recorrected

commit 40f05ca2ce95226ab17af02765182a141018593e
Author: Silmathoron <Silmathoron@users.noreply.github.com>
Date:   Tue Dec 6 17:05:07 2016 +0100

changed implementation to match #557
@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 9, 2016

Contributor

@Silmathoron This is strange indeed. I can fetch the two commits you mention above from Github as part of your t_aeif_ref branch, so they are on their servers, but somehow do not get displayed. Could you check with Github support? One little thing that I noted was that your commit ec8d479 (one parent in the merge) is marked with a different email address than all your other commits, but that shouldn't derail things.

Contributor

heplesser commented Dec 9, 2016

@Silmathoron This is strange indeed. I can fetch the two commits you mention above from Github as part of your t_aeif_ref branch, so they are on their servers, but somehow do not get displayed. Could you check with Github support? One little thing that I noted was that your commit ec8d479 (one parent in the merge) is marked with a different email address than all your other commits, but that shouldn't derail things.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 9, 2016

Contributor

Yes, this is what Jochen advised me a little earlier so I sent them a message...

Contributor

Silmathoron commented Dec 9, 2016

Yes, this is what Jochen advised me a little earlier so I sent them a message...

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Dec 9, 2016

Contributor

Ok, problem solved thanks to the support team!

Contributor

Silmathoron commented Dec 9, 2016

Ok, problem solved thanks to the support team!

@heplesser

@Silmathoron Thanks! I approve 👍.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 11, 2016

Contributor

@golosio Could you be the second reviewer for this PR?

Contributor

heplesser commented Dec 11, 2016

@golosio Could you be the second reviewer for this PR?

@golosio

This comment has been minimized.

Show comment
Hide comment
@golosio

golosio Dec 12, 2016

Contributor

Sure. Indeed, even though I did not comment, I already followed this discussion, as it was closely related to (and referred by) that in PR https://github.com/nest/nest-simulator/pull/557/files , and I've also seen the plot in #557 (comment) and the following comments.
I approve this PR.

Contributor

golosio commented Dec 12, 2016

Sure. Indeed, even though I did not comment, I already followed this discussion, as it was closely related to (and referred by) that in PR https://github.com/nest/nest-simulator/pull/557/files , and I've also seen the plot in #557 (comment) and the following comments.
I approve this PR.

@heplesser heplesser merged commit 95ca0ad into nest:master Dec 12, 2016

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment