neurophysik / jitcdde

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

Delays in helpers and possibility to add noise? #27

Closed jajcayn closed 2 years ago

jajcayn commented 4 years ago

Dear developer,

this is rather an implementation question than an issue. I am using jitc*de for almost all my simulations (computational neuroscience) and I love the speed! Recently, I wanted to make a leap and rewrite all my models to jitc*de however this one turned out to be particularly challenging. I do not want to bother you with all the equations, so I just state the necessary:

is this doable to implement in jitcdde? Without the dependence of Q on Q(t-tau) I would just make Q helper and compute it in each timestep. But, I need to somehow save its state. Then I was thinking about using the "hack" you thaught me: adding one dummy state variable y(n) which would "store" firing rates, but I wasn't able to write anything to it during the integration itself. Is there a way how to add e.g. anchors (as in CHSPy) during integration to dummy state variable?

The second question is simple yes/no: is there a way how to implement stochastic differential equations with delays? Something like mixing jitcdde and jitcsde.

Thanks a lot! And as I said, love the jitc*de :)

Yes, versions:

jitcdde==1.5.0
jitcode==1.4.0
jitcsde==1.4.0
jitcxde-common==1.5.0
Wrzlprmft commented 4 years ago

Is there a way how to add e.g. anchors (as in CHSPy) during integration to dummy state variable?

No. While I could implement something like this, it would be quite some work and particularly it would require you to also specify the derivative of your firing rate. The latter probably voids all gain over the following solution:

If you can easily compute the temporal derivative of your firing rate, you can make it a dynamical variable. You would only have to ensure that it does not digress from your actual firing rate. A brute-force sure-fire implementation is a dynamical variable that exhibits a simple fixed-point dynamics with the firing rate Q as the fixed point: n = −λ (ynQ), where λ controls the speed of convergence. For a sufficiently high λ, you should have yn = Q for all practical purposes.

Another way would be to compute the past firing rate as a helper in the present, i.e., at time t, you compute Q(tτ). This has the disadvantage that you may compute the exact same number repeatedly at different points in time, but that computation may not be your bottleneck.

The second question is simple yes/no: is there a way how to implement stochastic differential equations with delays? Something like mixing jitcdde and jitcsde.

No. Both SDEs and DDEs are quite difficult to solve adaptively to begin with requiring specific solvers, algorithms, and data structures. Combining both would be quite some work and require a suitable algorithm for SDDEs (which did not exist the last time I checked). One crucial problem when adding Brownian noise to DDEs (or ODEs for that matter) is that the noise has an arbitrarily fine structure and the integrator has to adapt its step size to this.

However, be aware that your real noise may not be like this at all, but have a cut-off frequency. In that case, JiTCSDE may not be what you need to begin with, but you can treat your noise like an input signal for JiTCDDE (using a dummy variable and so on).

And as I said, love the jitc*de :)

I am glad to hear that. Remember to cite me.

jajcayn commented 4 years ago

If you can easily compute the temporal derivative of your firing rate, you can make it a dynamical variable. You would only have to ensure that it does not digress from your actual firing rate. A brute-force sure-fire implementation is a dynamical variable that exhibits a simple fixed-point dynamics with the firing rate Q as the fixed point: _ẏn_ = −λ (_yn_−Q), where λ controls the speed of convergence. For a sufficiently high λ, you should have _y_n = Q for all practical purposes.

Good thinking! I tried it on dummy system and exactly, with sufficiently high λ it really converges fast... it might be unusable for stiff systems, but that's shouldn't be much of a problem! image

However, be aware that your real noise may not be like this at all, but have a cut-off frequency. In that case, JiTCSDE may not be what you need to begin with, but you can treat your noise like an input signal for JiTCDDE (using a dummy variable and so on).

yup, that was my go-to solution if this wasn't possible, I'd pre-generate noise beforehand and use it afterward as an input signal. For some systems, I use Ornstein-Uhlenback as a noise source, so I was thinking I'd solve O-U using jitcsde for my dt (which are known beforehand), then use CHSPy to add it to the delayed system and just "sample" from it during actual integration

I am glad to hear that. Remember to cite me.

Sure thing! I've read your Chaos paper, I am planning on doing that :)

Thanks again for insightful responses

Wrzlprmft commented 4 years ago

it might be unusable for stiff systems, but that's shouldn't be much of a problem!

If you have a stiff system, you’ll have a problem anyway, as JiTCDDE is not made for those to begin with.

I was thinking I'd solve O-U using jitcsde for my dt (which are known beforehand) […]

What do you mean by dt here? The step size? If yes, what would you gain from that?

jajcayn commented 4 years ago

What do you mean by dt here? The step size? If yes, what would you gain from that?

Here I just meant that I know at what sample size I need to analyse the modelling results. So e.g. now I am using hand-written Euler-Maruyama, I impose some integration dt and after that I just sub-sample for my sampling dt. With jitc*de I do not care and about integration dt since it's adaptive, I just need my time-points to be evenly sampled.

In any case, I was thinking of doing proper SDE of O-U with "infinitely fast noise" and then just sample from it subsequently. In my head, it made sense that it is better than just generating dW beforehand and integrating O-U at the same time as the whole system.

Wrzlprmft commented 4 years ago

I see. However be aware that: