Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Avoid returning state variables #290

Open
CGMossa opened this issue May 5, 2023 · 7 comments
Open

Avoid returning state variables #290

CGMossa opened this issue May 5, 2023 · 7 comments

Comments

@CGMossa
Copy link

CGMossa commented May 5, 2023

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
Copy link
Member

richfitz commented May 5, 2023

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
Copy link
Member

richfitz commented May 5, 2023

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

@CGMossa
Copy link
Author

CGMossa commented May 5, 2023

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
Copy link
Member

richfitz commented May 5, 2023

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
Copy link
Author

CGMossa commented May 5, 2023

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}) + D_I × \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
Copy link
Member

richfitz commented May 5, 2023

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
Copy link
Author

CGMossa commented May 5, 2023 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants