For loop expression evaluation in DifferEq solver

I am trying to evaluate strings within a for loop within an R script using JuliaCall::julia_eval . While I was able to accomplish this in R using the deSolve package, I am running into issues when converting the code to one that is compatible with Julia. I build a vector of string expressions which I pass to a function using a for loop or lapply. These expressions have variables that are defined in the global environment. How can I get Julia to evaluate these differential equations and have the variables update in each time step ? Example below

Combine <- c(" - 1*0.4545*(H2O2^1) - 1*27000000*(`$OH`^1)*(H2O2^1)", " - 1*3100000000*(`1,4-dioxane`^1)*(`$OH`^1)", 
    " - 1*33000*(TOC^1)*(`$OH`^1)", "2*0.4545*(H2O2^1) - 1*3100000000*(`1,4-dioxane`^1)*(`$OH`^1) - 1*33000*(TOC^1)*(`$OH`^1) - 1*27000000*(`$OH`^1)*(H2O2^1) - 1*8500000*(`$OH`^1)*(`HCO3-`^1) - 1*390000000*(`$OH`^1)*(`CO3 2-`^1)", 
    " - 1*8500000*(`$OH`^1)*(`HCO3-`^1)", " - 1*390000000*(`$OH`^1)*(`CO3 2-`^1)"
State <- c(H2O2 = 0.000294, `1,4-dioxane` = 0, 
TOC = 0, `$OH` = 0, `HCO3-` = 0.003766104, `CO3 2-` = 0.0000167638711956647)
f <- JuliaCall::julia_eval("
    function (dY,t,state)
        for i in 1:6
          dY[i] = eval(Meta.parse(Combine[i]))
    end
    end")

 sol = diffeqr::ode.solve(f,state,tspan, abstol=1e-8, reltol=1e-8)

Sorry if this is presented incorrectly I am new to Julia and this forum!

The first issue is this function definition. You cannot just choose any function definition that you want. It should be function (dY,state,p,t). With that change it should work.

However, if you want it to be fast, you shouldn’t eval inside of the Julia function, and you should instead build the whole string to eval. You can build that string directly in R without any Meta.parse, similarly to how you’d write a big Stan function.

1 Like

Thanks Chris. When I made the changes I noticed that only certain variables are updating at each time step. For example, I notcied that HCO3- and CO32- arent changing from their initial condition listed in state. I know that as the time steps I am forming $OH and therefore both dY[5] and dY[6] should eventually not be zero, however, I am not getting this when I run the code. I am not sure what to do here. I have attached 2 different ways I am able to evaluate differential expressions in R which do update each variable correclly, but I cant figure out hot to convert these to somethign compatible with Julia, Using “with” is crucial here.

ODEcreater2 <- function(t, state, parameters){
with(as.list(c(state)),{
for (i in 1:length(Species$Species)) {
dY[i] <- eval(parse(text=Combine[i]))}
dY[which(Species$Species == ‘H2O’)] = 0
return(list(c(as.numeric(dY))))
} )}

ODEcreater2 <- function(t, state, parameters){
with(as.list(c(state)), {
dY <- lapply(Combine, function(y) eval(parse(text=y)))
dY[which(Species$Species == ‘H2O’)] = 0
return(list(c(dY)))
})}

Occasionally I also run into issues where Julia eval and get “UndefVarError: Combine not defined” even though I know it is. Any help here would be greatly appreciated.

The function works fast enough for now with the eval and Meta.parse in Julia. The Combine vector was built in R using paste. Thanks for all your help!

I mean this in the kindest way possible: I do not understand what you are writing or why you would write it like this. It would both be much more efficient to just generate a string and to use that string, just like one would do with Stan. Also, if you had said string, I could debug it just by looking at the generated code. With the approach that you are taking, I am not sure how I could help and I would not recommend people to do it like this.

Instead, the diffeqr docs show you can write fast functions like:

f <- JuliaCall::julia_eval("
function f(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")

so if you had something like

c("10.0*(y-x)","x*(28.0-z) - y","x*y - (8/3)*z"), you can do some find-replace and then apply to build a string like above all from R, and then just call one julia_eval on it. Anything else is much more complicated than necessary and will also be slower.

2 Likes

Thanks Chris. Totally understand your concern. The problem is that I have over 100+ unique differential equations (du) and up to 150 unique variables (u). Instead of writing out each expression indiviually I thought it would be relatively easy to build strings using grep, which I was able to accomplish and store as the vector “Combine”. It is relatively easy also to evaluate this expression in DeSolve given a list of individual conditions using either a For loop or apply as I showed above. I was just looking for a way to do this in Julia. I want to avoid writing out u[1-150] explicitly. Is there no way to accomplish this what I showed above for ODEcreater1 and ODEcreater2 in Julia? Please let me know if this does’t make sense

Hey Chris. I think I see what you are saying now. I will work on the find-replace and touch back

I am not advocating you write it out either. I am saying you should just generate it. The for loop and all of that should take place on the R side so it’s only done once, and then the flat form is just a string that is evaluated and quick to execute. There is nothing to be gained by running the loops on the Julia side with R embedded objects, but if you do the exact same thing but all on the Julia side to build a string then you get a flat representation which won’t need any control flow and will be fast.

1 Like

Thanks Chris. I was able to do covert the code as you reccomended. I wanted to do some benchmarking on the speed so I manually entered everything into a Julia call/ Julia eval as shown below:

f <- JuliaCall::julia_eval("
function f(du,u,p,t)
du[1] = - 10.4545(u[1]^1) - 127000000(u[4]^1)(u[1]^1)
du[2] = - 1
3100000000*(u[2]^1)(u[4]^1)
du[3] = - 1
33000*(u[3]^1)(u[4]^1)
du[4] = 2
0.4545*(u[1]^1) - 13100000000(u[2]^1)(u[4]^1) - 133000*(u[3]^1)(u[4]^1) - 127000000*(u[4]^1)(u[1]^1) - 18500000*(u[4]^1)(u[5]^1) - 1390000000*(u[4]^1)(u[6]^1)
du[5] = - 1
8500000*(u[4]^1)(u[5]^1)
du[6] = - 1
390000000*(u[4]^1)*(u[6]^1)
end")

tspan <- list(0, 1E-1)
sol = diffeqr::ode.solve(‘f’,state,tspan,saveat=times, alg = “Rosenbrock23()”)

My question is now that I have this vector, what would you recommend for the best way to create the differential equations similar to above. Would you suggest the following?

Combine = c(" - 10.4545(u[1]^1) - 127000000(u[4]^1)(u[1]^1)", " - 13100000000*(u[2]^1)(u[4]^1)",
" - 1
33000*(u[3]^1)(u[4]^1)", "20.4545*(u[1]^1) - 13100000000(u[2]^1)(u[4]^1) - 133000*(u[3]^1)(u[4]^1) - 127000000*(u[4]^1)(u[1]^1) - 18500000*(u[4]^1)(u[5]^1) - 1390000000*(u[4]^1)(u[6]^1)",
" - 1
8500000*(u[4]^1)(u[5]^1)", " - 1390000000*(u[4]^1)*(u[6]^1)"

f <- JuliaCall::julia_eval("
function f(du,u,p,t)
for i in 1:6
du[i] = Combine[i]
end
end")

Thanks again for all your help and the work you have done with Julia!

Once again I apologize if these questions seem trivial. Regarding my last question, is it possible to pull a vector like Combine into Julia_call::Julia_eval? If so do I need to define it as a global variable? Basically the only thing I have left to do is to create the differential equations from that Combine vector that I have. Just reaching out to see if you knew the best way to go about this?

If you do it on the Julia side, you need to call the equivalent Julia function and define the array of functions as a global variable. It might be more difficult, and it has absolutely no performance advantages (and it has disadvantages) to do it like this.

What you should do now is loop through the strings and add "du[i] = " _ str[i] + "\n" appending to your current strings, and then combine them in R to a single string, then add `function f(du,u,p,t)" + str[i] + “end”, and that there is your full string to evaluate.

Thanks. I am working on it now but would like a bit more clarification its a bit vague:

  1. By add "du[i] = " _ str[i] + "\n" do you mean for as example for i=1 create a new string as du[1] = " 1 0.4545 (u[1]^1) - 1 27000000 (u[4]^1) *(u[1]^1) + \n. If so I can do that. If not can you explicitly

  2. I think I follow you here. Make finial expression as 1 string like “du[1]=…\n du[2] = …\n du[3]=…”

  3. Once again I am lost here on exactly what to type in. I am trying to use your response from https://stackoverflow.com/questions/47937451/generate-function-in-a-loop-for-differentialequations-in-julia as guidance. Maybe I just don’t the default format. If you could help explain one last time that would be greatly appreciated. I appreciate your time and if I am still lost I’ll reach out elsewhere for help.

Hey Chris. Any updates here would be greatly appreciated… If you could give a quick example of the expression to be included in the Julia_eval expression that would help me tremendously.

As mentioned in the two duplicate StackOverflow questions on this, the answer is just to do:

sprintf("function f(du,u,p,t)\n%s\nend", paste(Combine, collapse="\n"))

to Combine and build the string in R (given you’ve built the appropriate strings for each expression). Everything else is completely over-complicated and will result in slower code anyways. Maybe we can add a small blurb to the diffeqr docs on this, but it’s literally just string handling in R and nothing to do with Julia or diffeqr.