Early termination of ODE integration with callback

Just change the initial condition to be static arrays.

1 Like

I had a hunchā€¦

Another thought on making this run faster: It should be possible reduce the initial integration time by starting closer to the stable periodic orbit. The orbit is likely a perturbation of the equilibrium of the autonomous system, maybe even a slow-drift quasi-steady state of the seasonally-forced system, since the seasonal variation operates on a time scale much slower than the infection dynamics.

Iā€™ll look into this.

thanks @John_Gibson and @ChrisRackauckas.

I changed the initial conditions to a static labelled array with SLVector (following the example for the Lorenz equation here GitHub - SciML/LabelledArrays.jl: Arrays which also have a label for each element for easy scientific machine learning (SciML)), but this has slowed down fitting even more than with LArray. Not sure whether this is a problem of LabelledArrays? Iā€™ll also look into a SArray (without the labels), but havenā€™t had the time yet.

As mentioned my actual model is more complex than the SIR toy model I have given as an MWE. I have 12 compartments, none of which can be omitted, unfortunately. The C compartment also canā€™t be omitted because it tracks the cumulative incidence, which is needed to calculate the incidence between the time steps (diff), which is the only state observed through the data and thus necessary for fitting.

Open an issue with the reproducer and we can dig into the performance here. I think v1.6 had some detrimental changes to static arrays I havenā€™t pinned down.

For the 3d SIR model, I was able to speed up the convergence towards the stable periodic orbit by starting from the equilibrium of the non-seasonally-varying (\eta = 0) system instead of the [S, I,R] = [9999, 1, 0] initial condition. The nontrivial equilibrium of the \eta = 0 system is S=(\sigma/\beta)N_{pop}, I=(1/\omega)(1/\sigma - 1/\beta)N_{pop}, R=N-S-I. For the parameters you initially specified this works out to about [S, I,R] = [7143, 39, 2818]. The I(t) time series for the two different initial conditions are
timeseries

The nonlinear solve for the periodic orbit converges directly (five iterations) if you use the equilibrium as an initial guess.

Of course this might all be different for your real model of 20 variables. But there is probably a better initial guess for the steady state than whatever corresponds to [N-1, 1, 0].

A few more thoughts:

  1. My time-integrations and nonlinear solves are relatively fast: a couple hundred microseconds for integrating one period, and a few milliseconds for the whole nonlinear solve from the equilibrium initial guess. And thatā€™s with totally wasteful time-integration (out-of-place, storing whole time series, variable-sized vectors). How is it that your parameter-fitting computations take ten hours?
  2. You can probably get a whole lot of speed-up simply by increasing the tolerances in your time integrations from 1e-08 to 1e-03. If your system is stable and you only need a few digits accuracy in the output, thereā€™s no need for such tight tolerances. A higher tolerance will allow bigger time steps and faster runs. Depending on the order of accuracy of your integration scheme you might get an order of magnitude speed-up.
  3. Care to share the actual 20-variable model with me? Iā€™d need the equations to find a better iniital guess for that system.
1 Like

Canā€™t we use an homotopy from a previously found problem?

The original poster is using Bayesian inference algorithms to fit model parameters to time series data. The algorithms there, I believe, are basically random explorations of parameter space from specified distributions. It does seem like even some zeroth-order homotopy could really help, e.g. storing a library of parameter-solution pairs, and using the solution of the nearest library parameter for the initial guess.

I take the opportunity to say that BifurcationKit implements Parallel Multiple Standard/Poincare Shooting as well (see [docs])(https://rveltz.github.io/BifurcationKit.jl/dev/periodicOrbitShooting/). The stability (Floquet exponents) is also implemented. You can evaluate the functional to give it to your optimiser. There are several examples but perhaps this one and this one are closer to the post.

In any case, I found the content of this post really cool.

4 Likes

Thanks everyone for the feedback. I have to solve another issue with the library for the bayesian parameter inference first, but will get back to this post with code.

1 Like