neurophysik / jitcdde

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

Intended way to make time-dependent parameters #46

Closed mwappner closed 2 years ago

mwappner commented 2 years ago

Having a dynamical system smoothly or sharply vary one of its parameters over time is not uncommon. The former is usually meant either as colored noise or an input, while the latter is usually a perturbation to the system. As far as I could gather from Open and Closed issues, the docu and examples, currently the intended way to do this jitcdde_input.

It is however is not obvious what tool to use to do a parameter perturbation, if you are not sure what you are looking for. I simply propose adding an example Mackey-Glass system where β has a sharp change in the middle of the integration. If this sounds good, I'd be willing to PR the following example if you think it's worth adding it to the repo. I can't think of a good filename that would point a newcomer to this example if they want to do parameter perturbations, though. I hope this complies with both Python and jitcdde good practices.

### Makey-Glass with time-dependent parameter

import numpy as np
from jitcdde import y, t, jitcdde, jitcdde_input, input
from chspy import CubicHermiteSpline
from matplotlib import pyplot as plt

# Defining parameters
# --------------

τ = 15
n = 10
β0 = 0.25
β1 = 0.4
γ = 0.1

# Using input to modify a parameter mid-run
# --------------

input_times = np.arange(0, 800)
t0 = 400 # time of perturbation

parameter_over_time = lambda t: β0 if t<t0 else β1
parameter_spline = CubicHermiteSpline(n=1)
parameter_spline.from_function(parameter_over_time, times_of_interest = input_times)

# Defining Dynamics
# -----------------

β = input(0)
f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) ]

# Actual Integration
# ------------------

DDE = jitcdde_input(f, parameter_spline)
DDE.constant_past(1)
DDE.adjust_diff()
times = DDE.t + np.arange(0, input_times[-1])
data = np.vstack([DDE.integrate(time) for time in times])

# Plotting
# --------

fig, (ax1, ax2) = plt.subplots(2,1, sharex=True)

parameter_spline.plot(ax1)
ax1.set_title('Parameter over time')
ax1.set_ylabel('β')

ax2.plot(times, data, label='perturbed')
ax2.set_xlabel('time')

Output: perturbed mackey-glass time series

Wrzlprmft commented 2 years ago

There already is an example linked under More Features and Examples. I expanded this section and tried to make the options more clear. I am open for suggestions on how to better put or phrase things. You can say best what would have helped you.

Also, this does not obviate the usefulness of your example. So please feel free to make a pull request. I would name the file mackey_glass_parameter_jump.py and reference it in the aforementioned section. Alternatively and if you are willing to put in the work, you can make an entire documentation section explaining what is going on.

Finally, we can consider to present discuss several alternatives to implement this, i.e., using jitcdde_input, using a callback, using a change of a control parameter or an explicit dependency with conditional. Since your example is only a toy example, I wouldn’t say that there is a best solution to this. However, depending on the application, either of these approaches can become impossible or at least unfeasible (which is why they exist in the first place).

I have little to comment on your code except that I would not import plt (because you should hardly ever use it) and import subplots directly. Also, your image doesn’t match the code completely: The title of the second plot is missing from the code and instead there is an unused label.

mwappner commented 2 years ago

The additions to the Examples section look really good! I'll give it a go at implementing the other solutions and open a PR, I guess we can discuss there?

You can say best what would have helped you.

Yup, the part you just added is what I would have needed. Since you are asking for feedback, here's two comments, but as far as I'm concerned, it could just stay as it is:

Mackey–Glass with Jumps shows how to use the jump method.

I'd add "to introduce [approximate] discontinuities in the state of the system".

Simple Neutral and Neutral show how to implement neutral DDES. The latter additionally shows how to optimise a DDE with repeating delays, making it considerably faster.

I'm not sure what you mean by repeating delays. The example only has one delay. Also, the link to Neutral is wrong, you are linkinig to simple_neutral. And I'd change DDES for DDEs, idk which one is correct.

And an absolutely irrelevant typo: you have two full stops at the end of the first sentence here:

If you want a parameter to be a straightforward function of time, you can just implement this symbolically directly in the derivative.. For an example, see the “regular implementation” here.

Wrzlprmft commented 2 years ago

I'll give it a go at implementing the other solutions and open a PR, I guess we can discuss there?

Great and sure.

Since you are asking for feedback, here's two comments […]

Thanks. Everything should be fixed now.

mwappner commented 2 years ago

Woops, I didn't expect you to be so quick about it, I fixed the typos in #47 as well. I guess you'll have to resolve a few minor conflicts when merging :/

mwappner commented 2 years ago

Solved by #47 and a few commits in this thread.