diffeqr
diffeqr copied to clipboard
SciML
named vectors in Julia function
Hi
I just started using diffeqr to fit ODE models in R. Is there a way to have named vectors in the function? Indexing works:
library(diffeqr)
theta <- c(beta=.25, gamma=1/5, mu=1/(70*365), sigma=1/365)
maxtime <- 3650.0
init <- c(100000,0,0)
func_julia <- function(init, param, times) {
beta <- param[1]
gamma <- param[2]
mu <- param[3]
sigma <- param[4]
S = state[1]
I = state[2]
R = state[3]
N = S + I + R
dSdt = -beta*I*S/N + mu*N - mu*S + sigma*R
dIdt = beta*I*S/N*S - gamma*I - mu*I
dRdt = gamma*I - mu*R - sigma*R
return(c(dSdt, dIdt, dRdt))
}
enviro <- diffeqr::diffeq_setup()
problem <- enviro$ODEProblem(func_julia, init, c(0.0, maxtime), theta)
solution <- enviro$solve(problem , enviro$AutoVern7(enviro$KenCarp4()),
saveat=1.0, abstol=1e-8, reltol=1e-6)
but using named vectors in R throws an error:
func_julia <- function(init, param, times) {
beta <- param[["beta"]]
gamma <- param[["gamma]]
mu <- param[["mu"]]
sigma <- param[["sigma"]]
S = state[1]
I = state[2]
R = state[3]
N = S + I + R
dSdt = -beta*I*S/N + mu*N - mu*S + sigma*R
dIdt = beta*I*S/N*S - gamma*I - mu*I
dRdt = gamma*I - mu*R - sigma*R
return(c(dSdt, dIdt, dRdt))
}
Error in param[["beta"]] : subscript out of bounds
Error: Error happens in Julia.
Named vectors would be great especially for inference and if there are dozens of parameters in the model. Thanks Fabienne
Currently, named vectors will have their names dropped by JuliaCall, and I'm not sure whether Julia has any data structure (similar) corresponding to named vectors in R?
I also wonders whether it is possible to use named list instead of named vectors in this example. For example, the following seems to work for me and have the same result with the unnamed version:
func_julia_named <- function(init, param, times) {
beta <- param[["beta"]]
gamma <- param[["gamma"]]
mu <- param[["mu"]]
sigma <- param[["sigma"]]
S = state[1]
I = state[2]
R = state[3]
N = S + I + R
dSdt = -beta*I*S/N + mu*N - mu*S + sigma*R
dIdt = beta*I*S/N*S - gamma*I - mu*I
dRdt = gamma*I - mu*R - sigma*R
return(c(dSdt, dIdt, dRdt))
}
problem_named <- enviro$ODEProblem(func_julia_named, init, c(0.0, maxtime), as.list(theta))
solution_named <- enviro$solve(problem_named , enviro$AutoVern7(enviro$KenCarp4()),
saveat=1.0, abstol=1e-8, reltol=1e-6)
Notice the usage of as.list(theta) in construction of the problem_named. This should work as long as theta itself is not involved in arithmetic operation as a whole.
Thanks. It works for me as well if I solve the problem as is. But if I "accelerate" the problem with diffeqr::jitoptimize_ode() it throws me an error:
Error: Error happens in Julia.
KeyError: key 'b' not found
Stacktrace:
[1] getindex at .\dict.jl:467 [inlined]
[2] #1 at .\none:0 [inlined]
[3] iterate at .\generator.jl:47 [inlined]
[4] join(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Base.Generator{String,ModelingToolkit.var"#1#2"}, ::String) at .\strings\io.jl:294
[5] join(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Base.Generator{String,ModelingToolkit.var"#1#2"}) at .\strings\io.jl:293
[6] sprint(::Function, ::Base.Generator{String,ModelingToolkit.var"#1#2"}; context::Nothing, sizehint::Int64) at .\strings\io.jl:105
[7] sprint at .\strings\io.jl:101 [inlined]
[8] join at .\strings\io.jl:300 [inlined]
[9] map_subscripts at C:\Users\fabie\.julia\packages\ModelingToolkit\CEB9b\src\variables.jl:14 [inlined]
[10] _broadcast_getindex_evalf at .\broadcast.jl:648 [inlined]
[11] _broadcast_getindex at .\broadcast.jl:621 [inlined]
[12] #19 at .\broadcast.jl:1046 [inlined]
[13] ntuple at .\ntuple.jl:41 [inlined]
[14] copy at .\bro
Is there a way around this? This would be a huge advantage when fitting complex models. Thanks, Fabienne
Currently, named vectors will have their names dropped by JuliaCall, and I'm not sure whether Julia has any data structure (similar) corresponding to named vectors in R?
https://github.com/SciML/LabelledArrays.jl
@fkrauer I have some intuition for what causes the error: the key "beta" of the list in R (will be dict in Julia) somehow gets iterated and be decomposed into 'b'+'e'+'t'+'a', and diffeqr looks for elements with keys 'b', 'e', 't', and 'a' in the list, which does not exist... so the error occurs. The behavior makes sense, because by numerical index, like a[1:3], what we usually mean is actually [a[1], a[2], a[3]], but this is not true with string keys.
I implemented an "uggly" version using LabelledArrays. The version can be "beautified" if some code are integrated into JuliaCall. However, it seems that this version still doesn't work with diffeqr::jitoptimize_ode(). I guess it may be due to the julia_call("getindex" ....) thing. But I also feel that even I fix this, a problem similar with as.list(theta) will still be prevent diffeqr::jitoptimize_ode().
func_julia_LabelledArrays <- function(init, param, times) {
beta <- julia_call("getindex", param, quote(beta)) ## param[["beta"]]
gamma <- julia_call("getindex", param, quote(gamma)) ## param[["gamma"]]
mu <- julia_call("getindex", param, quote(mu)) ## param[["mu"]]
sigma <- julia_call("getindex", param, quote(sigma)) ## param[["sigma"]]
S = state[1]
I = state[2]
R = state[3]
N = S + I + R
dSdt = -beta*I*S/N + mu*N - mu*S + sigma*R
dIdt = beta*I*S/N*S - gamma*I - mu*I
dRdt = gamma*I - mu*R - sigma*R
return(c(dSdt, dIdt, dRdt))
}
library(JuliaCall)
julia_command("using LabelledArrays")
autowrap("LabelledArrays.LArray")
problem_LabelledArrays <- enviro$ODEProblem(func_julia_LabelledArrays, init, c(0.0, maxtime),
julia_eval("LVector(beta=.25, gamma=1/5, mu=1/(70*365), sigma=1/365)"))
solution_LabelledArrays <- enviro$solve(problem_LabelledArrays, enviro$AutoVern7(enviro$KenCarp4()),
saveat=1.0, abstol=1e-8, reltol=1e-6)
Running:
library(diffeqr)
enviro <- diffeqr::diffeq_setup()
theta <- c(beta=.25, gamma=1/5, mu=1/(70*365), sigma=1/365)
maxtime <- 3650.0
init <- c(100000,0,0)
func_julia_LabelledArrays <- function(init, param, times) {
beta <- julia_call("getindex", param, quote(beta)) ## param[["beta"]]
gamma <- julia_call("getindex", param, quote(gamma)) ## param[["gamma"]]
mu <- julia_call("getindex", param, quote(mu)) ## param[["mu"]]
sigma <- julia_call("getindex", param, quote(sigma)) ## param[["sigma"]]
S = state[1]
I = state[2]
R = state[3]
N = S + I + R
dSdt = -beta*I*S/N + mu*N - mu*S + sigma*R
dIdt = beta*I*S/N*S - gamma*I - mu*I
dRdt = gamma*I - mu*R - sigma*R
return(c(dSdt, dIdt, dRdt))
}
library(JuliaCall)
julia_command("using LabelledArrays")
autowrap("LabelledArrays.LArray")
problem_LabelledArrays <- enviro$ODEProblem(func_julia_LabelledArrays, init, c(0.0, maxtime),
julia_eval("LVector(beta=.25, gamma=1/5, mu=1/(70*365), sigma=1/365)"))
solution_LabelledArrays <- enviro$solve(problem_LabelledArrays, enviro$AutoVern7(enviro$KenCarp4()),
saveat=1.0, abstol=1e-8, reltol=1e-6)
diffeqr::jitoptimize_ode(enviro,problem_LabelledArrays)
fails with
Error: Error happens in Julia.
ArgumentError: invalid index: :beta of type Symbol
Stacktrace:
[1] to_index(::Symbol) at .\indices.jl:297
[2] to_index(::Array{Operation,1}, ::Symbol) at .\indices.jl:274
[3] to_indices at .\indices.jl:325 [inlined]
[4] to_indices at .\indices.jl:322 [inlined]
[5] getindex(::Array{Operation,1}, ::Symbol) at .\abstractarray.jl:1060
[6] docall(::Ptr{Nothing}) at D:\OneDrive\Computer\Documents\R\win-library\4.0\JuliaCall\julia\setup.jl:168
[7] (::RCall.var"#32#33"{LangSxp,Ptr{LangSxp},Ptr{EnvSxp},Base.RefValue{Int32}})() at C:\Users\accou\.julia\packages\RCall\Qzssx\src\eval.jl:84
[8] disable_sigint at .\c.jl:446 [inlined]
[9] tryEval at C:\Users\accou\.julia\packages\RCall\Qzssx\src\eval.jl:81 [inlined]
[10] reval_p(::Ptr{LangSxp}, ::Ptr{EnvSxp}) at C:\Users\accou\.julia\packages\RCall\Qzssx\src\eval.jl:96
[11] reval_p at C:\Users\accou\.julia\packages\RCall\Qzssx\src\eval.jl:95 [inlined]
[12] rcall_p(::RObject{ClosSxp}, ::Array{Operation,1}, :
The error seems to be because A.x in a LabelledArray lowers to A[:x], but I think JuliaCall might have a wrapper over the type so this dispatch is lost in the autowrap version?