SciML / diffeqr

Solving differential equations in R using DifferentialEquations.jl and the SciML Scientific Machine Learning ecosystem
Other
141 stars 14 forks source link

Direct definition of fast functions from R #22

Closed ChrisRackauckas closed 1 year ago

ChrisRackauckas commented 4 years ago
f <- function(u,p,t) {
  du1 = p[1]*(u[2]-u[1])
  du2 = u[1]*(p[2]-u[3]) - u[2]
  du3 = u[1]*u[2] - p[3]*u[3]
  return(c(du1,du2,du3))
}
u0 <- c(1.0,0.0,0.0)
tspan <- c(0.0,100.0)
p <- c(10.0,28.0,8/3)
prob = de$ODEProblem(f, u0, tspan, p)
odesys = de$modelingtoolkitize(prob)
jul_f = de$ODEFunction(odesys)
prob2 = de$ODEProblem(jul_f, u0, tspan, p)
sol = de$solve(prob2,de$Tsit5()) # overhead!

f3 <- JuliaCall::julia_eval("
function f3(du,u,p,t)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end")
JuliaCall::julia_assign("u0", u0)
JuliaCall::julia_assign("p", p)
JuliaCall::julia_assign("tspan", tspan)
prob3 = JuliaCall::julia_eval("ODEProblem(f3, u0, tspan, p)")
sol = de$solve(prob3,de$Tsit5()) # no overhead but messy!

system.time({ for (i in 1:100){ de$solve(prob2   ,de$Tsit5()) }})
system.time({ for (i in 1:100){ de$solve(prob3   ,de$Tsit5()) }})

Notice how modelingtoolkitize still has overhead of the function is defined like this. @Non-Contradiction is there a way to directly build Julia functions from R without overhead? Until then, we'll be using this workaround via jitoptimize_ode.

Non-Contradiction commented 3 years ago

It is difficult to directly build julia functions from R functions in general. I see that modeling toolkit has a macro system. It should be okay to have a macro system in R (non-standard evaluation) and translate the R expression into DSLs in julia so build a julia function from a restricted set of R expressions.

ChrisRackauckas commented 3 years ago

That's what we're doing. The only issue is that calling

odesys = de$modelingtoolkitize(prob)
jul_f = de$ODEFunction(odesys)

which builds the Julia function has a lot of overhead, while just sending those pieces over and evaling has a lot less overhead. So building the same Julia function with the same Julia commands has no overhead if you eval it but has overhead if you call the functions directly through JuliaCall.

ChrisRackauckas commented 1 year ago

The JIT functions do this.