How to integrate python class attributes

[jacobcook1995]

I've started taking a look at how to integrate my model over time. The thing I'm having trouble getting my head around is how to use python's ODE solvers with dimensions that are stored as class attributes.

I've been considering using scipy's builtin odeint solver.

y = odeint(model, y0, t, args=(k,))

Where y0 is a vector of initial values, t is the integration time span, and model is the function to be integrated. My understanding is that model should have a form where the amount to update each dimension is returned as part of a vector.

def model(y,t,k):
    dydt = [-k * y, k*y^2]
    return dydt

I understand how to set up simple cases like the above. However, in our case a lot of what we will be interested in integrating won't be supplied as a vector but as attributes of a class (which are stored as vectors).

I can write a function where all pool updates are calculated and returned.

def update_all_soil_pools(soil_carbon: SoilCarbonPools, soil_nitrogen: InorganicNitrogenPools, *args):
    """Calculate updates for all pools in the soil module."""
    carbon_pools_update = soil_carbon.calculate_soil_carbon_updates(soil_nitrogen, *args)
    nitrogen_pools_update = soil_nitrogen.calculate_inorganic_nitrogen_updates(soil_carbon, *args)
    return np.array([carbon_pools_update, nitrogen_pools_update])

This could in theory be integrated (by adding an extra vector argument pools to be integrated over). However, soil_carbon and soil_nitrogen would not be updated between integration steps. The easy fix for this is to add a step that updates each set of pools before the update amount is returned.

def update_all_soil_pools(soil_carbon: SoilCarbonPools, soil_nitrogen: InorganicNitrogenPools, *args):
    """Calculate updates for all pools in the soil module."""
    carbon_pools_update = soil_carbon.calculate_soil_carbon_updates(soil_nitrogen, *args)
    nitrogen_pools_update = soil_nitrogen.calculate_inorganic_nitrogen_updates(soil_carbon, *args)

    # Update pools before anything is returned
    soil_carbon.update_soil_carbon_pools(carbon_pools_update)
    soil_nitrogen.update_soil_nitrogen_pools(nitrogen_pools_update)
    return np.array([carbon_pools_update, nitrogen_pools_update])

The above should definitely work if we were using a simple time step, i.e. $y_1 = y_0 + dt*\Delta y$. However, the scipy integrators generally use more complex integration schemes (mainly Runga-Kutta methods), where y values are calculated for multiple time points (e.g. the midpoint) in order to improve the accuracy of the calculation. My concern is how to make sure that the attributes of SoilCarbonPools are not updated more or less times than they should be.

Something I considered was recreating SoilCarbonPools for every time step, but this feels like it would be costly and should be avoided if possible.

I was wondering if there was an established best practice for tackling this kind of problem? And if anyone had any suggestions about how to tackle this?

[davidorme]

I think it might be worth dragging James Rosindell into this discussion - it was an area he had some great thoughts about in the original VR meeting. And Jaideep!

More broadly, I'm not quite sure I understand the problem:

• We have the state at $t=0$, stored in class attributes or possibly in a Data instance.
• The model uses odeint to find the state at $t=1$, and those values are then stored in the Data instance (possibly alongside the previous values if we are keeping a timeseries, possibly replacing them).

Bear in mind that ODEs are witchcraft to me....

[davidorme]

It isn't a problem if we were planning our update procedure to be

y(t=1)=y(t=0)+Δt∗dy(t=0)dt,

where Δt is the time between the current model update and the previous one. In that case we wouldn't use odeint at all.

Ohhhh. That is what I wasn't getting. So if you wanted say daily values within the monthly 'tick' of the main timesteps then odeint will return an array of a shape (len(y0), n_days). I don't see a problem with storing that as an internal class attribute that is refreshed by update for use within each tick.

If we want to allow people to access that data, then I've been thinking that there might be two data export points in the loop.

• The main one is for variables where we keep a time series across cells in Data. Those can be exported as requested at the end of the simulation (although we might want incremental to catch data if things explode).
• The other one might happen within the loop and allow "within update" variables to be store in data and then written out to file (and possibly into the same file by appending) before they are overwritten by the next update. That would be slower but if someone wanted a daily file of soil stuff, it could still happen.

[jacobcook1995]

Are we expecting to store instantiated classes in data? Or more just arrays of floats representing attribute values?

[davidorme]

The latter - it is all stored as an xarray.Dataset. But there is no reason you can't store a model array in there. Having said that, the only reason to do that and not as a model attribute is to make it easy to export, which is not exactly a priority for internal integration data right now!

I think we're expecting a model run to only have a single instance of each BaseModel subclass throughout the simulation, right?

[jacobcook1995]

Yeah definitely don't want to be creating multiple SoilModel instances. Integration takes place within this class anyway so the problem doesn't really exist at this level.

The problem is the SoilCarbonPools class (and other similar classes that I'm yet to define). The various functions that provide the updates to the soil model as a whole are contained in these classes. So, to find the net update to the soil module up to date instances of each class valid for the specific time step are needed.

However, the key problem is that most integration schemes don't use a tick, but instead a time step. Where a tick represents a single time point for which a calculation is performed, a time step is a range for which the change in variables should be linearly approximated over. Notably this linear approximation will always involve multiple calculation steps. Often some of these steps repeat time points.

To take 4th order Runga-Kutta as an example, the initial rate of change (i.e. for $t_n$) is calculated. The rate of change at the midpoint (i.e. $t_n + \frac{\Delta t}{2}$) is then calculated twice (first using the initial rate of change to calculate the midpoint y value, the second time using the first midpoint rate of change to calculate the midpoint y value). Next a rate of change is calculated for the final time point (i.e. $t_{n+1}$ ). The y value for $t_{n+1}$ is then calculated by weighting the various rates of rates of change over the time step. This procedure is then repeated for every time step in the desired range. This results in every time point and every midpoint being calculated twice, which seems mad but works out as more efficient as you can get away with using a substantially longer time step without losing accuracy.

The problem is then how I can get odeint (or alternatives) to play nicely with the Pools classes. Their values can't really be incremented with every calculation (because some calculations don't step you forward in time), but not sure if there's a way to get odeint to be aware of class attributes and how they should be updated. This feels like someone should have tackled it before but I can't find anything concrete online.

It may be that I am misunderstanding where integration should be happening though, i.e. it should happen within each single class and then SoilModel just runs the lower level integrations. Though if that is the case I definitely need to rethink my planned class divisions.

[davidorme]

I think this discussion needs a whiteboard 😄

[jacobcook1995]

Thought about/discussed my difficulty a bit more so hopefully this statement of the problem will be a bit more clear.

In essence the odeint function accepts a few inputs.

odeint(function, y0, t, k)

1. function: the function describing the system of ODEs to be solved
2. y0: a vector of initial conditions.
3. t: details of the timing (i.e. how long to integrate for, step size, etc)
4. k: arguments for the function, this is basically a dictionary of parameters. Pretty certain that these should be constant for the entire integration.

As the initial condition y0 must be a vector of values, function must return a vector of dydx values for odeint to use to update the y over time.

The trouble I'm having is in working out how to fit this into the object based approach I have used for defining SoilCarbonPools. At present I am intending to make several classes like this (e.g. InorganicNitrogen). These classes will not be completely independent, i.e. calculating the appropriate update to SoilCarbonPools will require knowledge of the state of InorganicNitrogen. In some cases, this communication can take place infrequently, and so the classes can be integrated separately. However, a lot of the classes I intended to create will require much more frequent exchange of information between them, which will necessitate integrating all the relevant classes in a single integration.

I can see how you can setup odeint to work within a class, but cannot work out how to feed multiple classes into odeint and get have them update correctly as they are integrated over. I'm not sure if I have setup my classes inappropriately, or whether there is a clever trick I am missing.

I've though of a number of ways to potentially resolve this problem:

1. I discover an amazing trick to get odeint to properly update classes supplied to it.
2. I place every pool that needs to be integrated together into a single class, and then carry out the integration within this class. This would probably result in the soil model consisting of a giant class along with a few small classes for things that can be integrated separately.
3. I reduce the SoilCarbonPools class (and similar classes in future) to something that doesn't have associated functions and really just stores pools as attributes and sanity checks updates. All functions would then be independent of a class and so class instances would not have to be fed into odeint. Before each integration the pools contained in each class would have to be converted into a vector to be fed into odeint, and after the integration an update step would be needed to update the pools stored in the classes.

Option 2 feels like the most plausible solution to me. But I'm quite worried that it would result in a monstrously large class, so wonder if there is a smarter approach here.

@dalonsoa @alexdewar let me know if you have any thoughts on this. Cheers!

[alexdewar]

The documentation suggests using solve_ivp instead: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

This might be a stupid question, but are you doing separate integrations or are they all combined somehow?

[jacobcook1995]

In terms of the soil module was planning to (eventually) do separate integrations where feasible to avoid having one huge set of equations to integrate. The problem is that a lot of the components of the soil module are tightly coupled which limits how much they can be divided. At the moment I feel like the soil module will likely end up with one big integration for the main soil nutrient cycles, a small integration for decomposition processes, a then a set of processes that only need to be updated very rarely (e.g. soil pH changes) and which maybe don't warrant integrating at all.

Should be said that I haven't got anywhere concrete with implementing the above, and I could be way off in terms of how I am thinking about this in my head