RateSystem

[raphael-roemer]

We would like to add functionality for Rate dependent forcing functionality and easy system creation of non-autonomous systems

[reykboerner]

Thanks a lot for working on this!

I very much like the idea of starting with the general mathematical formulation \dot x = f(x(t), lambda(t)), and have made a suggestion of how this could be reflected more in the way a user would define a RateSystem.

One thing I am unsure about is whether we just need a constructor function RateSystem that returns either a CoupledODEs or CoupledSDEs, or whether we need an additional type in order to pass the necessary information about the rate system to the methods we want to add.

Maybe we can brainstorm the methods we would like to add here and think about whether they work with a CoupledODEs or whether we need more.

[reykboerner]

It should be sufficient to use DynamicalSystemsBase.

[reykboerner]

In my opinion, these are quite many positional arguments, which makes it less user-friendly. What do you think of reorganising this in the following structure:

RateSystem(ds::ContinuousTimeDynamicalSystem, rate_protocol::function)

where ds can be a CoupledODEs or a CoupledSDEs and rate_protocol is your λ function that determines how the parameters change over time. This function could have the structure

rate_protocol(u, p_lambda, t; t_start=0.0, t_end=Inf, kwargs...)

i.e. it would contain all the input arguments needed to define the rate forcing function. Behind the scenes, this would then create and return a new CoupledODEs or CoupledSDEs with the explicit time dependence specified by rate_protocol.

[reykboerner]

The naming here is just a suggestion, but generally I think we should give more descriptive names than just lambda or tnf (clarity over brevity!)

[raphael-roemer]

This sounds great, thanks Reyk! However, rate_protocol should not depend on u, in my opinion.

[oameye]

https://docs.sciml.ai/JumpProcesses/stable/

[reykboerner]

Hey, what's the status on this @raphael-roemer @ryandeeley ?
If you continue working on this, don't forget to merge the main branch into the RateSystem branch first, since the main branch has changed in the meantime.

[ryandeeley]

Hey, what's the status on this @raphael-roemer @ryandeeley ? If you continue working on this, don't forget to merge the main branch into the RateSystem branch first, since the main branch has changed in the meantime.

Hello, Raphael and I met again last week to pick this up again. We agree that it would be nice to pass some used-defined auto_sys::CoupledODEs and \lambda::Function - which respectively describe 1) the reference system in an autonomous form, and 2) the time-dependent evolution of the system parameters - to a function RateSystem() that will return a modified nonauto_sys::CoupledSDEs describing the associated nonautonomous version of the system. We also agree that it'd be nice if RateSystem() contains the fewest positional arguments as possible. I've had a crack at writing this generally, and applying to an example case. The user has to provide

• the base system auto_sys::CoupledSDEs
• the function \lambda::Function that describes the time-dependent evolution of the system parameters.

We are proposing that \lambda() should return all system parameters at each time instance, not just those that are explicitly time-dependent, to avoid the (ambiguous) choice of which parameters in the parameter vector of the original system need to be overwritten. Instead, the idea is that the function RateSystem() - which has auto_sys::CoupledODEs and \lambda::Function as positional arguments - creates a modified time-dependent drift function by evaluating the drift function auto_sys.integ.f(u,p,t) at p=\lambda(t) at each time instance. Using the useful function remake(), one then creates an ODEProblem which is identical to the corresponding ODEProblem of auto_sys::CoupledODEs in all components except for the drift function - which becomes the modified time-dependent drift function. The parameter vector and the initial time of the ODEProblem also change in this new construction. This ODEProblem is then converted back to a CoupledODEs, which is ultimately what RateSystem() then returns.

As I've written it, the positional arguments of RateSystem() are

• auto_sys::CoupledODEs
• \lambda::Function
• r::Float64
• t0::Float64

The first two arguments are hopefully clear. The argument r describes a scalar control parameter specifying the scaled rate of change of the parameters - because ultimately we are considering ODEs of the form dxₜ/dt = f(xₜ,λ(rt)). The argument t0 is the initial time that the returned CoupledODEs takes - which is fairly important in this nonautonomous context.

One keyword argument of RateSystem() is p_lambda::Vector=Float64[] which (in a meta way) contains the parameters for the evolution of the time-dependent system parameters. For instance, the minimal/maximal value of the parameter in a typical ramping scenario. I also find it most logical to have t_start::Float64=-Inf and t_end::Float64=Inf as keyword arguments of RateSystem(). These entries allow for a quick modification of the user-provided \lambda()::Function in order to "clip" the period of non-autonomous parameter evolution - specifically, the modified parameter function will assume the same values \lambda(t) in the interval t\in[t_start,t_end] but for t < t_start or t > t_end assumes \lambda(t_start) and \lambda(t_end) respectively. Note that this clipping of the \lambda()::Function occurs independently of the rate value r prescribed: the value of r has the effect of squeezing/stretching this particular clipped parameter evolution.

Note that there are two non-user-level functions modified_drift() and rate_protocol() that RateSystem() calls in order to construct the modified time-dependent drift function.

I tested this proposed RateSystem() on the paradigmatic one-dimensional model for studying R-tipping, and it seems to return the solutions one expects (in reference to the behaviour before / after the critical r-value for R-tipping with the particular time-dependent parameter evolution considered - in this case this is r = 4/3, for R-tipping of the trajectory that in backward-time t->-\infty limits to the attractor x = -1 of the "past" autonomous system). I can also see whether this works nicely for a two-dimensional system.

I've attached a .txtof the .jl file where the function RateSystem() is written up and the tests are made, so let me know what you think of this whenever you get time!

RateSystemDraft.txt

[reykboerner]

Hey @ryandeeley, thanks for the update! Could you push your proposed changes to this branch? This is what git is designed for, and pushing here will automatically run the tests. (The idea is that this draft pull request is continuously developed until we are happy to merge it with the main branch)

[ryandeeley]

Hey @ryandeeley, thanks for the update! Could you push your proposed changes to this branch? This is what git is designed for, and pushing here will automatically run the tests. (The idea is that this draft pull request is continuously developed until we are happy to merge it with the main branch)

Yes, we can do that, although I only wrote this yesterday and haven't explicitly discussed it with Raphael yet, who is also drafting some RateSystem() function. We plan to compare each of our proposals, and I believe we're meeting again tomorrow, so we could push something after then.

[oameye]

@raphael-roemer just a reminder you can run the docs locally with Liveserver.jl

[Datseris]

Hi all, sorry this took a while, especially given I don't have that many comments. I'm sorry, if I knew I'd only have minor comments I'd done this sooner...

Primarily, let's please try to stick to a specific terminology, and I'd really like us to keep to the one that we have agreed in the video calls after a long back and forth: Δp, Δt and r = Δp/Δt is the effective rate.

The other point to discuss is the computational performance of this. It is really bad... :( I really think we should brainstorm how to improve it.

That being said, improving the performance of the source code is completely irrelevant of hte interface. So I'd recommend to implement the changes I've illustrated here, merge this PR asap and then update the other PR with the functions that actually use this type.

Meanwhiole I'll think of the performance of this and how to improve it.

[Datseris]

If you don't like Δp, Δt, then we can go with : forcing_length (Δt), forcing_size (Δp) and forcing_start.

Also, @raphael-roemer please comment: how does forcing_start and t0 interplay? t0 is the starting time of the (autonomous) dynamical system. forcing_start should be fully independent of this time. both should be measured on units of real time, and we must always check that forcing_start >= t0 when we construct the RateSystem.

[Datseris]

@raphael-roemer can you please comment where in the source code you need the export of hte referrneced_sciml_model that you mentione din the email? I don't see it.

[Datseris]

Oh, I've missed this:

nonauto_traj = trajectory(nonauto_sys.system, T, x0);

this is not the way to go. There is a well defined API for what a DynamicalSystem is, and this RateSystem should be exactly a dynamical system. Make it subtype ContinuousTimeSystem so that the call trajectory(nonauot_sys) works directly.

You guys have done this work already when you created the CoupledSDEs type. You already have all the skills and knowledge for this, no need to wait for me to comment on this, just go ahead and finalize things.

[Datseris]

I've just updated the benchmarks to use BenchmarkTools for clarity. I find a large difference:

@btime nonauto_traj = trajectory($(nonauto_sys.system), T, x0);
@btime nonautoExpl_traj = trajectory($nonauto_sysexpl, T, x0);

gives

  256.100 μs (9763 allocations: 171.65 KiB)
  67.800 μs (14 allocations: 19.32 KiB)

this is absolutely massive, and the number of allocations shows very clearl problems in the implementation beyond the algorithm design. @raphael-roemer have you profiled your code and checked the basic components of julia optimization, (such as type instability, unecessary allocations, etc.), as e.g., mentioned here: https://modernjuliaworkflows.org/optimizing/ ?

[reykboerner]

Thanks @Datseris for your review! So apart from the computational speed, do you find the interface design makes sense?

If you don't like Δp, Δt, then we can go with : forcing_length (Δt), forcing_size (Δp) and forcing_start.

I wasn't in the video calls, so I didn't know about these discussions about the naming. I prefer the written out labels because Δt looks like a time step and Δp is not general enough. I strongly think that we should not restrict to monotonic forcing protocols in the long-run. forcing_scale more accurately describes what this setting does: it scales the magnitude of the RateConfig you provided. Also forcing_magnitude works in my opinion if you prefer.

this is not the way to go. There is a well defined API for what a DynamicalSystem is, and this RateSystem should be exactly a dynamical system. Make it subtype ContinuousTimeSystem so that the call trajectory(nonauot_sys) works directly.

Happy to do this. One question: Is it a problem that the RateSystem type itself contains a CoupledODEs type as RateSystem.system? Because essentially the RateSystem is just a container for the nonautonomous system + the forcing config. There is an equivalence between a dynamical system constructed via RateSystem vs. the same system constructed as an explicitly time-dependent CoupledODEs.

[Datseris]

forcing_scale is fine as well and I agree with you with nonmonotonic forcings. I stated in one of the comments i'll contribute reversible forcing used to study overshooting (doing this right now for the lectures, might as well put it here)

---

Happy to do this. One question: Is it a problem that the RateSystem type itself contains a CoupledODEs type as RateSystem.system?

no, this is no problem at all. Essentially make sure that you subtype RateSystem <: ContinuousTimeDynamicalSystem and then extend the functions to dispatch to the field, a-la how you did four copled sdes

for f in (:step!, :set_state!, ....) # all api functions here
     @eval DynamicalSystemsBase.$(f)(rs::RateSystem, args...; kw...)= $(f)(rs.system, args...; kw...)
end

[Datseris]

So apart from the computational speed, do you find the interface design makes sense?

yes provided the changes I've outlined above are included!

[reykboerner]

@raphael-roemer I'll work on this tonight unless you are already on it!

[oameye]

this is absolutely massive, and the number of allocations shows very clearl problems in the implementation beyond the algorithm design. @raphael-roemer have you profiled your code and checked the basic components of julia optimization, (such as type instability, unecessary allocations, etc.), as e.g., mentioned here: modernjuliaworkflows.org/optimizing ?

@raphael-roemer Hey, if needed I am happy to jump in here, to optimize the code. Just let me know. Probably best to do this in a seperate PR

[raphael-roemer]

this is absolutely massive, and the number of allocations shows very clearl problems in the implementation beyond the algorithm design. @raphael-roemer have you profiled your code and checked the basic components of julia optimization, (such as type instability, unecessary allocations, etc.), as e.g., mentioned here: modernjuliaworkflows.org/optimizing ?

@raphael-roemer Hey, if needed I am happy to jump in here, to optimize the code. Just let me know. Probably best to do this in a seperate PR

This would be great - thank you! I am quite busy before Friday, so I can only work on this afterwards.

[reykboerner]

@oameye if you have some time it would be great to get your help on making RateSystem a proper subtype of ContinuousTimeDynamicalSystem. I started here by modifying the CoupledODEs(prob, diffeq...) function from here but there are still import errors etc.

[Datseris]

you don't need anything fro mthe CoupledODEs source code!!!! Just make sure you rpopagate all APi functions to the stored internal system field!!!!

for f in (:step!, :set_state!, ....) # all api functions here
     @eval DynamicalSystemsBase.$(f)(rs::RateSystem, args...; kw...)= $(f)(rs.system, args...; kw...)
end

as for making it a "type" just subtype ContinuousTimeDynamicalSYstem this should "just work"

[reykboerner]

Okay thanks!

[Datseris]

Heeeeeell yea baby. I now have a fully-automated, completely generic codebase for generating plots like in Ritchie et al that show when a system is tracking, returning but not tracking, or always tipping. Here is the plot:

I've plotted it as two plots because I couldn't be bothered to make the subtraction of the 2 plots to have 1 plot, but the plot above is exactly the same as the plots in Ritchie et al 2023:

The plots are for the simple 2D stommel model.

The code is completely generic , fully agnostic of the dynamical system, and fully automated by utilizing global_continuation. In essense, by providing a dynamical system, a global continuation tells you the frozen system attractors. Then, for any given Dp, Dt you can simulate and tell, objectively (with a distance function and threshold ε as in AttractorsViaProximity), to which attractor the system converges to after its rate forcing. The colors in my plot is simply the attractor ID, with purple the "starting" attractor.

I will contribute this method here once the RateSystem is done. My code is only 20 lines of code by the way hahahaaaa man I love Attractors.jl it simplifies everything so much!

---

p.s. you may want to join https://discourse.julialang.org/t/juliadynamics-monthly-meetings-round-2/132357 tomorrow because I'll organize a minisymposium in next years JuliaCon in Germany.

[Datseris]

okay i couldn't take it I improved the code so that it immediatelly produces this figure:

[reykboerner]

This is great - I also think that the RateSystem interface can nicely be linked to Attractors.jl to extend the functionality to a nonautonomous context (by considering snapshots of the RateSystem + simulating the non-autonomous system).

So I suggest I get this PR cleaned up, we merge, and then we can go from there with solving the speed issue/adding functionality.

[Datseris]

type unstable.

[Datseris]

but perhaps it doesn't matter, as this is only a configuration, so ignore for now,.

but since this is only a configuration perhaps worth merging the RateCOnfig inside this directly? just take its three fields and add them here?

[Datseris]

The less Types the better - always!