epirecipes / sir-julia

Various implementations of the classical SIR model in Julia
MIT License
195 stars 39 forks source link

Add RCall RHS #115

Closed sdwfrost closed 2 months ago

sdwfrost commented 2 months ago

In the spirit of diffeqr, plug in a R RHS for an ODE or discrete system rather than a Julia one.

sdwfrost commented 2 months ago

The following is a minimal example:

using OrdinaryDiffEq
using ModelingToolkit
using RCall

R"""
sir_ode_r <- function(u,p,t){
    S <- u[1]
    I <- u[2]
    R <- u[3]
    N <- S+I+R
    beta <- p[1]
    cee <- p[2]
    gamma <- p[3]
    dS <- -beta*cee*I/N*S
    dI <- beta*cee*I/N*S - gamma*I
    dR <- gamma*I
    return(c(dS,dI,dR))
}
"""

function sir_ode_jl(u,p,t)
    robj = rcall(:sir_ode_r, u, p, t)
    return convert(Array,robj)
end

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem{false}(sir_ode_jl, u0, tspan, p)

I'd like to use modelingtoolkitize on this, but this currently doesn't work OOTB:

@named sys = modelingtoolkitize(prob)
ERROR: MethodError: no method matching sexpclass(::Num)

Tracked in https://github.com/SciML/ModelingToolkit.jl/issues/2735

sdwfrost commented 2 months ago

Closed with 8045d5c020b47bbb8443622edff787f749cf1f92