neurophysik / jitcdde

Just-in-time compilation for delay differential equations
Other
57 stars 14 forks source link

Setting delays results in weird crash #41

Closed mwappner closed 2 years ago

mwappner commented 2 years ago

I was trying to do a small parameter search varying the delay, so I decided to use it as a control parameter in my system and change it using set_parameters instead of compiling a new model every time. The following pieces of code are not the actual code I use, but a minimal example that reproduces the same behavior.

The first naïve approach to the problem is just compiling a model without extra keywords and letting the library figure out what delays the system has:

import numpy as np
from jitcdde import jitcdde, t, y
from symengine import Symbol

tau = Symbol('tau')

b = 20
f = [- y(0) + b / (1 + y(0, t-tau)**2 )]

taus = [ 5, 10, 20]
DDE = jitcdde(f, control_pars=[tau])

for t in taus:
    DDE.constant_past([1])

    DDE.set_parameters(t)
    DDE.step_on_discontinuities()

    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

This raises the exception

ValueError: At least one delay depends on time or dynamics; cannot automatically determine steps.

on the line DDE.step_on_discontinuities(), which is fine. I understand that maybe set_parameters doesn't trigger the automatic delay calculations again, so you need to manually specify them. Since every run has a different delay, I though the smart thing was to manually set the attribute each cycle like so:

DDE = jitcdde(f, control_pars=[tau])

for t in taus:
    DDE.delays = [0, t]
    DDE.constant_past([1])

    DDE.set_parameters(t)
    DDE.step_on_discontinuities()

    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

(I know jitcdde appends a 0 to the delays, so I included that, but removing it doesn't alter the following behavior). Running this results in a crash of some sort: running on console it returns Segmentation fault (core dumped) and when trying in an interactive setting like Spyder I get

Kernel died, restarting
Restarting kernel... 

[SpyderKernelApp] WARNING | No such comm: 5538ec1a9be111ecb3f349eb681f0c91

which I don't think is relevant but may be of some help for you.

Even adding a list to the delays on initialization like so results in the same crash: DDE = jitcdde(f, control_pars=[tau], delays=taus).

Regarding the problem itself, it's easily solved if I abandon the idea of having the model use the delay it's actually using in that run and just give it the full list of delays I'm gonna use:

DDE = jitcdde(f, control_pars=[tau], delays=taus)

for t in taus:
    DDE.constant_past([1])

    DDE.set_parameters(t)
    DDE.step_on_discontinuities()

    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

which I'm not sure is the best call, as it may be worse or slower at integrating than the case where the model knows the actual delay used. I wouldn't expect the fix to allow the line DDE.delays = [0, t], but rather just prevent the user from setting the attribute.

On a separate note, since delays gets appended with an extra 0, and it actually uses the user-passed object, rather than a copy of it, the last bit of code I showed will run for delay values [0, 5, 10, 20]. I'm not sure this will still happen in the unreleased version that changes the very line that appends the 0, but if it does, it'd be reasonable to warn the user in the docs. or make a copy.

jitcdde version: 1.8.0 (I'm eagerly awaiting a new release!)

Wrzlprmft commented 2 years ago

The first naïve approach to the problem is just compiling a model without extra keywords and letting the library figure out what delays the system has:

[…]

This raises the exception

ValueError: At least one delay depends on time or dynamics; cannot automatically determine steps.

on the line DDE.step_on_discontinuities(), which is fine. I understand that maybe set_parameters doesn't trigger the automatic delay calculations again, so you need to manually specify them.

Almost. The problem is that JiTCDDE cannot calculate the delay in this case at all since this would have to dissect the dynamical equations for that. Sure, it’s straightforward in your case, but a general solution to this would be very tedious. I now expanded the error message to include a control-parameter delay as a potential cause.

Since every run has a different delay, I though the smart thing was to manually set the attribute each cycle like so:

[…]
    DDE.delays = [0, t]
[…]

(I know jitcdde appends a 0 to the delays, so I included that, but removing it doesn't alter the following behavior). Running this results in a crash of some sort: running on console it returns Segmentation fault (core dumped) […]

This is because setting delays did not trigger a recalculation of the max_delay, which in turn was responsible for the segfault. I fixed this now. The best way to make your example work without updating would be something like:

[…]
DDE = jitcdde(f, control_pars=[τ])

for tau in taus:
    DDE.max_delay = tau
    DDE.delays = [tau]
        […]

I'm eagerly awaiting a new release!

Yeah, I should probably do one of these.

mwappner commented 2 years ago

Finally updated my library and it works nicely!