Build DiffEq functions from Expr

Hello, my question is both related to building functions from expressions and DiffEq domain.

I am working on a simulation engine based on ODE solver and the input problems for the solver are provided as JSON files, which includes initial values and all the expressions to build ODE system and other related functions. I can easily parse this file, replace simbols in Expr with p[i], u[i] and have an Expr blocks to build, for example, ode_func - ODE system and saving_func - the output function for ODEProblem:

ode_expr = quote
  A = u[1]
  B = u[2] 
  k1 = p[1]
  comp1 = p[2]
  v1 = k1 * A * comp1
  du[1] = -v1
  du[2] = v1
end

saving_expr = quote
  A = u[1] 
  B = u[2] 
  k1 = integrator.p[1]
  comp1 = integrator.p[2]
  v1 = k1 * A * comp1
  [A, B, k1, v1]
end

Those exprs are runtime data but initially I know the arguments and output types for functions, so I can build ode_func and saving_func and form ODEProblem:

@eval ode_func(du,u,p,t)  = $ode_expr
@eval saving_func(u,t,integrator) = $saving_expr

out = SavedValues(Float64, Vector{Float64})
scb = SavingCallback(saving_func, out)
prob = ODEProblem(ode_func, u0, tspan, p0, callback=scb)

However when I try to solve the prob: solve(prob,Tsit5()) I encounter world age error:
The applicable method may be too new: running in world age 27808, while current world is 27811.
I know that it can be solved with invokelatest but here saving_func is not called directly but provided to another function solve…
I’ll be grateful if someone helps me solve this issue. Is “parse |> eval” the only (correct?) way to build functions here (or macros can do it in a nicer way) and how to solve world age error in this case?

1 Like

If you let execution go back to global scope after the eval, then I would have thought that you shouldn’t get world-age errors. Is the second code-block you gave inside a function (or some other scope-defining block)?

I don’t think macros can help you as you only get your information on how to build the ODE at runtime, i.e. after macro-expansion time.

Anyway, if your ODEs could all be put into one general functional form, with parameter-sets distinguishing between then, then that would be best. Otherwise, I don’t think you can avoid the eval.

1 Like

@mauro3 thanks for response!
Yes, the second code-block is inside a function, which (as I belive) is the cause of the world age Error.

The general form of ODEs functions are known at compile time (input argument list and output type) but expressions to be used inside of those functions are given at runtime. I thought that placing all Expr into Dict and using macros to expand “general” function forms is doable, but I am not much experienced with macros.

eval could be an “ok” solution but I don’t know how to deal world age error in this case.

No, macros don’t work on values but the code they are applied on. Example:

a = :(x+1)
@m a 

the the macro m will only see a symbol :a and not what the variable a is bound to. Watch this: Jacob Quinn: What happens when - From parse time to compile time - YouTube

Concerning your problem, you will have to do the loop going over your dict of exprs, eval-ing and running the problem in global scope:

dict_of_expr = ...
for (k,e) in dict_of_exprs
  
  @eval ode_func(du,u,p,t)  = $ode_expr
  @eval saving_func(u,t,integrator) = $saving_expr

  out = SavedValues(Float64, Vector{Float64})
  scb = SavingCallback(saving_func, out)
  prob = ODEProblem(ode_func, u0, tspan, p0, callback=scb)
end

, don’t put it into a function. And if you really want to, then you have to use invoke.

2 Likes

I got the idea, but yes, I have to put this code in a function. How to use invoke in this case when saving_func,ode_func are called internally in when solving the ODEProblem?

I have to put this code in a function.

Why? You can also put it in one function, but not use it, let execution reach global scope, and then use it in another function.

Use:

help?> Base.invokelatest
  invokelatest(f, args...; kwargs...)

  Calls f(args...; kwargs...), but guarantees that the most recent method of f will be executed. This is useful in specialized circumstances, e.g. long-running event
  loops or callback functions that may call obsolete versions of a function f. (The drawback is that invokelatest is somewhat slower than calling f directly, and the
  type of the result cannot be inferred by the compiler.)

performance will likely suffer quite a bit, best check.

Pass (args...) -> Base.invokelastest(ode_func, args...) to ODEProblem.

Thanks, I see. Both solutions will likely affect the performance. I’ll check if it is comparative more performant to reach global scop than to use invokelatest

Annotating the return type, say invokelatest(f, args...; kwargs...)::Vector{Float64}, should also help a bit.

Hi @ivborissov. Did you ever solve this problem to your satisfaction? It appears I’m trying to do the same as you, 5 years later, and I’m also struggling with the world age problem.

I’m preparing a course for students that uses Stock-and-Flow diagrams, and I want to convert their Stock-and-Flow diagrams into a set of differential equations. As you say, I can parse their equations into a dynamical rule, however I can find absolutely no way of passing this dynamical rule to DynamicalSystems CoupledODEs. My problem is illustrated by the following sample code:

module SD
builder = eval(Expr(:(->),:(p,u),:(+(p[1] * u[1]))))

function demo()
	dynamical_rule! = function (du,u,p,t)
		du[1] = Float64(eval(builder(p,u)))
		nothing
	end

	# dynamical_rule!() should set value of du from 1.0 to 0.6 (= p[1]*u[1]):
	du,u,p,t = ([1.0],[2.0],[0.3],4.0)
	println("Values BEFORE calling dynamical_rule!(): (du,u,p,t) = $((du,u,p,t))")
	dynamical_rule!(du,u,p,t)
	println("Values AFTER  calling dynamical_rule!(): (du,u,p,t) = $((du,u,p,t))")
end
end

I need to construct the builder function (builder=eval(...)) within the demo() function - I just can’t do it in global scope. However, doing so throws the world age error. Does anyone have any good ideas?

Hi,
apologies for getting back to you so late. I think currently you can use RuntimeGeneratedFunctions for this.
https://docs.sciml.ai/RuntimeGeneratedFunctions/stable/

Many thanks! I’ve actually now given up on building an automatic stock-and-flow compiler, and have instead found a way in which beginner students can work directly with DynamicalSystems. However, this package is well worth remembering - thank you very much! :grinning:

1 Like