mrc-ide / odin

ᚩ A DSL for describing and solving differential equations in R
https://mrc-ide.github.io/odin
Other
102 stars 12 forks source link

Avoid returning state variables #290

Open CGMossa opened 1 year ago

CGMossa commented 1 year ago

I've a SIR model with many, many cells. They influence each other via a contact matrix. Thus for me, only a select cells are of interest. Is there a way to suppress the output of deriv variables to model$run?

For a workaround, is it possible to know how I can take the current output, edit it, and pass it to the rest of the odin machinery?

An example of syntax could be:

  output(S[]) <- NULL

To suppress the outputs of S[1], ..., S[N].

richfitz commented 1 year ago

This is something we hit with the larger COVID model that the department worked on - a stochastic transmission model with ~35,000 cells. Unfortunately it's not a great fit for how deSolve works - we've got an entirely new interface via odin.dust which supports this sort of thing, though its primary focus was on stochastic models and running in parallel.

With that approach you can set an index over the output variables and get back only a subset. The overall interface is very different, being much more stateful rather than the "evaluate this IVP" that the original odin supports. We will likely reconcile this at some point, but that will only happen after we get dust on CRAN.

If sticking with original odin, have a look at the transform_variables method which you can use to process your output

richfitz commented 1 year ago

See also https://mrc-ide.github.io/odin-dust-tutorial/ (press "S" on the slides for speaker notes, which contain commentary)

CGMossa commented 1 year ago

Thank you for the answer.

Thus, the answer would be to get the full output:

model_generator <- odin(...)
model <- model_generator$new(...)
traj <- model$run(...)

But then use transform_variables to eliminate the state variables:

traj <- model$transform_variables(traj)
traj$S <- NULL
traj$I <- NULL
traj$R <- NULL
as.data.frame(traj)

And proceed with only the output stuff plus t / step (depending on if deriv or update was used).

I'm doing something that is related to wildlife disease spread. So I'm discretising a landscape and want to make the disease spread uncoupled from the discretisation. Then I had hoped to add gradient and maybe even diffusion terms to mimic migration and dispersal resp. And for that I need a PDE solver more than I need a stochastic solver.

But dust and the presentation you've linked looks very complete. I'll revisit this most definitely, as stochastic models would be the next next step for me.

richfitz commented 1 year ago

This sort of thing always depends a lot on the details. If you'll be solving many copies of the same ODE for your PDE discretisation, consider dust because then at least you can do that in parallel (there's a simple adaptive RK ODE solver available there but it may or may not be suitable for your needs - I think we have partial state updating implemented there though, which might help provide enough building blocks).

If you want speed, keep away from data.frame's - they are very slow and lists of matrices/arrays, or simply matrices/arrays are the way to go.

One thing we have done with the original odin models is to run a vector 1:n_state through transform variables to get the indices of the elements of interest, then repeatedly apply that during the main part of the run, this is faster than calling transform_variables each time and you can adapt that to your needs.

I can't imagine that we'll implement a PDE solver directly - I don't think that there are any general purpose ones available for R as it is. If you just want to do diffusion, I've done that reasonably before in the context of a set of odes using fftw to convolve the diffusion kernel with the target, but you do need to keep a very wide domain to stop it getting out of control. That would work reasonably well with the dust approach I suspect.

CGMossa commented 1 year ago

Wow, thank you for the insight. You're absolutely right about data-frame vs. matrix/array. First, I like it when issues are concluded with an answer to the original query, and thus I tried to conclude it. And also I'm plotting directly after running so hehe :D

I can't imagine that we'll implement a PDE solver directly - I don't think that there are any general purpose ones available for R as it is. If you just want to do diffusion, I've done that reasonably before in the context of a set of odes using fftw to convolve the diffusion kernel with the target, but you do need to keep a very wide domain to stop it getting out of control. That would work reasonably well with the dust approach I suspect.

This part generally went over my head. Would you mind expanding on it a little bit? What do you mean by target? Is this some sort of FEM thing? I've generally not solved differential questions with convolution before. Except by hand using Laplace transform as an undergrad... In my case, I want to encompass that wild boars move about, and thus I'd have a $\partial I / \partial t = (\texttt{regular transmission terms}) + DI × \nabla{(x,y)}^2 I(x,y,t)$ (excluding the gradient term migration). I don't see what can be called the diffusion kernel here, or target :D Even if you have just a trusted reference on this that would be very helpful ☺️

richfitz commented 1 year ago

See https://academic.oup.com/sysbio/article/59/6/619/1711291#112747614 and https://github.com/richfitz/diversitree/blob/master/src/quasse-eqs-fftC.c for the approach (it's not pretty or documented well, it was a long time ago). There may be better ways, and I forget where I adapted the algorithm from.

From memory, my diffusion problem did not vary with space, so the diffusion part of the system is just a gaussian kernel with fixed variance proportional to the time step; to apply this, you compute the of y (in your case this is I(t)) with that kernel, and an efficient way of doing this is to take the fourier transform of both I and the kernel, multiply these together, then take the inverse fourier transform.

If your diffusion process varies in space (so with I) this won't work. The time part of the system is still solved with an ODE solver.

I've had a quick search and I cannot remember anywhere the reference I used to motivate this. I know I spent a bunch of time checking it in Mathematica and it worked for my case! In any case, if you end up with a discretisation in space, consider the dust approach as then you can at least parallelise that bit

CGMossa commented 1 year ago

Noted!

I thought one only used Laplace transform to deal with derivatives, but fourier transforms might be the way... I'm not well-versed in this stuff, but if it works, I'll explore it further.

Essentially I want to consider the ODE under different discretisation of the landscape. The reason for this is I can choose a discretisation that favours the wild life disease spread, or one that favours spill-over onto static locations (domestic farms). Sort of a what's the advantage of discretising differently depending on if you start out with the static population, or the dynamic population (there is interplay, but for this I'm just considering spill-over from wildlife (dynamic) ~> domestic (static)). In any case, each cell will have to calculate the influence from surrounding cells. The other choice is using voronoi-tessellation vs. regular cells; For the former, I vary number of cells, as to mimic different static population density, and then I want to compare the result of the distribution of spill-over for the density.

In short: My repeated thing is not over the same ODE, but different ODEs that should yield the same result.. Hopefully I won't need to compile for each configuration;

This was invaluable, and I've saved it this as an email to remember. But this issue is answered, thus I'll close this (I've got other issues unfortunately..)

On Fri, 5 May 2023, 10.36 Rich FitzJohn, @.***> wrote:

See https://academic.oup.com/sysbio/article/59/6/619/1711291#112747614 and https://github.com/richfitz/diversitree/blob/master/src/quasse-eqs-fftC.c for the approach (it's not pretty or documented well, it was a long time ago). There may be better ways, and I forget where I adapted the algorithm from.

From memory, my diffusion problem did not vary with space, so the diffusion part of the system is just a gaussian kernel with fixed variance proportional to the time step; to apply this, you compute the of y (in your case this is I(t)) with that kernel, and an efficient way of doing this is to take the fourier transform of both I and the kernel, multiply these together, then take the inverse fourier transform.

If your diffusion process varies in space (so with I) this won't work. The time part of the system is still solved with an ODE solver.

I've had a quick search and I cannot remember anywhere the reference I used to motivate this. I know I spent a bunch of time checking it in Mathematica and it worked for my case! In any case, if you end up with a discretisation in space, consider the dust approach as then you can at least parallelise that bit

— Reply to this email directly, view it on GitHub https://github.com/mrc-ide/odin/issues/290#issuecomment-1535920768, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAIDVSFWX6APWVVJWIMZXILXES3Y7ANCNFSM6AAAAAAXWVTJIA . You are receiving this because you authored the thread.Message ID: @.***>