Use of interpolation function in ode model written in R/diffeqr

New to Julia
,
1

Hi all

I am using Julia’s diffeq through R (diffeqr library) for the simulation of a simple SIR ode model. I want to interpolate a parameter value from the time in the ode model. I can set up a julia ODE model and call through R with the diffeqr library. However, I am not sure how to pass the interpolation function to the ODE? The instructions for julia suggest that any function can be introduced through the “parameter” argument of the ODE. How does this work if I call the function in R?

Below is my unsuccessfull attempt.

  1. julia code for the model (“juliatoymodel.jl”):
function model(dxdt, init , param , time)

  beta, gamma, mu, sigma, eta, phi = param[1]
  interpolate = param[2] 

  S,I,R = init
  N = S + I + R
  i = interpolate(time)

  beta_eff = i*beta * (1.0+eta*sin(2.0*pi*(time-phi)/365.0))
  FOI = beta_eff*I/N

  dxdt[1] = -FOI*S + mu*N - mu*S + sigma*R
  dxdt[2] = FOI*S - gamma*I - mu*I
  dxdt[3] = gamma*I - mu*R - sigma*R

end
  1. R code that calls the functions:
library(diffeqr)
library(JuliaCall)

# params
maxtime <- 3650.0
beta <- 0.25
gamma <- 1/5
mu <- 1/(70*365)
sigma <- 1/365
eta <- 0.2 
phi <- 10
theta <- c(beta, gamma, mu, sigma, eta, phi)

# init
init <- c(9999, 1, 0)

# setup environment
enviro <- diffeqr::diffeq_setup() 

julia_source("juliatoymodel.jl")

# interpolation function
julia_library("Interpolations")
julia_assign("in",  c(0,365,730,1095,1460,1816,2181,2546,2911,3276,3650))
julia_assign("out", c(1,1,1,1,0.9,0.8,0.7,0.8,0.9,1,1))

julia_assign("init", init)
julia_assign("param", list(theta, julia_command("interpolate = LinearInterpolation(in,out);")))
julia_assign("time", c(0.0, maxtime))

problem = julia_eval("ODEProblem(model, init, time, param)")
solution = enviro$solve(problem, enviro$AutoVern7(enviro$KenCarp4()), 
                         saveat=7.0, abstol=1e-8, reltol=1e-6)

I have tried to make an R list with [[1]] params and [[2]] the interpolation function, which is assigned to a Julia object “param”. Unfortunately this throws an error

Error: Error happens in Julia.
MethodError: objects of type Nothing are not callable

because julia_command("interpolate = LinearInterpolation(in,out);") returns NULL somehow. How can I pass the interpolation function correctly? I am looking for an R/JuliaCall solution, not a pure Julia solution. Thanks a lot
Bw, Fabienne

2

You can try to use the julia_eval function, which will return the result from julia to R. The julia_command function will not return the result from julia to R, so the implicit result is always NULL.

3

Thanks a lot for the fast reply. I tried julia_eval, but it throws an error:

Error: Error happens in Julia.
MethodError: objects of type Array{Float64,1} are not callable
Use square brackets [] for indexing an Array.

Also, julia_eval("param") returns the result of the function (the interpolation), instead of the function itself. So the main question, how to pass the function itself to the Julia ODE model (within the parameter argument), is still unresolved.

4

I’m not very familiar with the diffeq ecosystem, and not sure what should be provided as the argument.
But julia_assign will convert any R object into julia object (including R functions to julia functions) and give it a name, and julia_eval will try to convert any julia object to R object (including julia functions to R functions).
And we need to pay attention to julia callable object (I’m not sure whether it is involved in the question). Because JuliaCall sees the callable object not as a function, so it will not turn it into R function. (Because in R, there is no such thing as callable object.)
But you can always construct a julia function based on the callable object
So suppose interpolate is a callable object, you can do julia_eval("x -> interpolate(x)") to get a function back into R.

1 Like
5

thanks. I have run the same script, but without explicitly including the interpolation function in the parameter argument of the ODE, and it seems to work. I guess as long as the interpolation function is defined in the julia universe, it can be called even if it is outside the ODE function. There may be a better solution for this, but it seems to work for now.