That is great advice! Thank you very much!

Just to clarify, is it possible to return to the top level in the examples we posted?

That is great advice! Thank you very much!

Just to clarify, is it possible to return to the top level in the examples we posted?

It’s possible, but it would make for an awkward API. For example, you can do:

```
julia> expr = eval(:(x -> x + 1))
#11 (generic function with 1 method)
# We're back to the top-level here, so we can use the result of eval()
julia> solver(expr)
65
```

but presumably you want users to be able to do:

```
function my_user_code()
func = make_function_somehow()
# This is still inside a function, not at the top level
do_stuff_with_function(func)
end
```

which doesn’t involve returning to the top-level.

Yeah, I was afraid that would be the case, I just wanted to be sure.

So, again, thank you very much for your help. I have already implemented your advice and got a 2x gain in performance and reduced allocations by the same factor!

1 Like

Or (3) use GeneralizedGenerated. This is somewhat the textbook use case for it. To make the function pure you need to avoid closures in there and it should be good

2 Likes

Why are you building functions from expressions at runtime in the first place? Except for very specialized circumstances (e.g. writing an interpreter shell), this is usually not the right approach … usually you should be passing around functions, not expressions.

1 Like

I am using closures in these functions because they receive the solution of an SDE System for a given trajectory `n`

. So, for example, a function example could be:

```
(u_n, p) -> begin
X = t -> un(t; idxs = 1)
T = p.T
K = p.K
# the following block is given by the user
(X(T) > K) * X(T)
end
```

The closure allows the user to write expressions using a simple syntax. For example, they can reffer to the process X at time T for each trajectory simple by using `X(T)`

.

Yes, I am aware of that. I am doing this because I don’t want the user of my library to write functions but instead, write a really simple input using macros. Then, I use the expressions provided in the macros to build functions. That’s it.

Any insight is welcome! Thank you for your reply! Please let me know if you have any suggestions for me.

If you’re using a macro already, why don’t you just build the functions in the macro output? Then there’s no eval required.

As a dumb example:

```
macro foo(input_block)
f_template = quote
(u_n, p) -> begin
X = t -> un(t; idxs = 1)
T = p.T
K = p.K
end
end
push!(f_template.args[2].args[2].args, input_block)
esc(f_template)
end
julia> @macroexpand @foo (X(T) > K) * X(T)
quote
(u_n, p)->begin
X = (t->begin
un(t; idxs=1)
end)
T = p.T
K = p.K
(X(T) > K) * X(T)
end
end
```

Because there is one macro for describing the equations of a model and there is another separated macro where I describe the parameters involved in the equations (this macro is `@with_kw`

in `Parameters.jl`

package).

So building your `f_template`

variable is not as straightforward as it seems.

Does

```
(u_n, p) -> begin
f(t) = un(t; idxs = 1)
T = p.T
K = p.K
# the following block is given by the user
(X(T) > K) * X(T)
end
```

work? I thought closures in GG were fixed though, @thautwarm

1 Like

Exactly, that would probably solve the problem with GG. Also, this function is type stable:

```
f = (un, p) -> begin
X(t) = un(t; idxs = 1)
T = p.T
K = p.K
# the following block is given by the user
(X(T) > K) * X(T)
end
@code_warntype f(sol.u[1], params)
Variables
#self#::Core.Compiler.Const(var"#178#179"(), false)
un#3252::DiffEqBase.RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},true},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:tspan, :X₀, :μ, :σ, :μX, :σX),Tuple{Tuple{Float64,Float64},Float64,Float64,Float64,var"#120#125"{Float64},var"#121#126"{Float64}}},Nothing,DiffEqBase.SDEFunction{true,UniversalMonteCarlo.var"#2883#2884",UniversalMonteCarlo.var"#2885#2886",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniversalMonteCarlo.var"#2885#2886",Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},StochasticDiffEq.SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats}
p#3253::NamedTuple{(:T, :K, :r),Tuple{Float64,Float64,Float64}}
val::Float64
X::var"#X#180"{DiffEqBase.RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},true},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:tspan, :X₀, :μ, :σ, :μX, :σX),Tuple{Tuple{Float64,Float64},Float64,Float64,Float64,var"#120#125"{Float64},var"#121#126"{Float64}}},Nothing,DiffEqBase.SDEFunction{true,UniversalMonteCarlo.var"#2883#2884",UniversalMonteCarlo.var"#2885#2886",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniversalMonteCarlo.var"#2885#2886",Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},StochasticDiffEq.SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats}}
T::Float64
K::Float64
r::Float64
pv::Float64
Body::Float64
1 ─ $(Expr(:inbounds, true))
│ %2 = Main.:(var"#X#180")::Core.Compiler.Const(var"#X#180", false)
│ %3 = Core.typeof(un#3252)::Core.Compiler.Const(DiffEqBase.RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},true},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:tspan, :X₀, :μ, :σ, :μX, :σX),Tuple{Tuple{Float64,Float64},Float64,Float64,Float64,var"#120#125"{Float64},var"#121#126"{Float64}}},Nothing,DiffEqBase.SDEFunction{true,UniversalMonteCarlo.var"#2883#2884",UniversalMonteCarlo.var"#2885#2886",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniversalMonteCarlo.var"#2885#2886",Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},StochasticDiffEq.SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats}, false)
│ %4 = Core.apply_type(%2, %3)::Core.Compiler.Const(var"#X#180"{DiffEqBase.RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},true},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:tspan, :X₀, :μ, :σ, :μX, :σX),Tuple{Tuple{Float64,Float64},Float64,Float64,Float64,var"#120#125"{Float64},var"#121#126"{Float64}}},Nothing,DiffEqBase.SDEFunction{true,UniversalMonteCarlo.var"#2883#2884",UniversalMonteCarlo.var"#2885#2886",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniversalMonteCarlo.var"#2885#2886",Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},StochasticDiffEq.SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats}}, false)
│ (X = %new(%4, un#3252))
│ (T = Base.getproperty(p#3253, :T))
│ (K = Base.getproperty(p#3253, :K))
│ (r = Base.getproperty(p#3253, :r))
│ %9 = -r::Float64
│ %10 = (%9 * T)::Float64
│ %11 = Main.exp(%10)::Float64
│ %12 = (X)(T)::Float64
│ %13 = (%12 > K)::Bool
│ %14 = (X)(T)::Float64
│ %15 = (%11 * %13 * %14)::Float64
│ (pv = %15)
│ (val = %15)
│ $(Expr(:inbounds, :pop))
│ val
└── return pv
```

However, I am using the following scheme:

```
function getclosure(un::unType, idx::Int64) where unType
return t -> un(t; idxs = idx)
end
g = (un, p) -> begin
X = getclosure(un, 1)
T = p.T
K = p.K
# the following block is given by the user
(X(T) > K) * X(T)
end
@code_warntype g(sol.u[1], params)
Variables
#self#::Core.Compiler.Const(var"#184#185"(), false)
un#3256::DiffEqBase.RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},DiffEqNoiseProcess.NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},true},DiffEqNoiseProcess.RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},DiffEqBase.SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:tspan, :X₀, :μ, :σ, :μX, :σX),Tuple{Tuple{Float64,Float64},Float64,Float64,Float64,var"#120#125"{Float64},var"#121#126"{Float64}}},Nothing,DiffEqBase.SDEFunction{true,UniversalMonteCarlo.var"#2883#2884",UniversalMonteCarlo.var"#2885#2886",LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniversalMonteCarlo.var"#2885#2886",Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},StochasticDiffEq.SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats}
p#3257::NamedTuple{(:T, :K, :r),Tuple{Float64,Float64,Float64}}
X::Any
T::Float64
K::Float64
r::Float64
pv::Any
Body::Any
1 ─ (X = Main.process_closure(un#3256, 1))
│ (T = Base.getproperty(p#3257, :T))
│ (K = Base.getproperty(p#3257, :K))
│ (r = Base.getproperty(p#3257, :r))
│ %5 = -r::Float64
│ %6 = (%5 * T)::Float64
│ %7 = Main.exp(%6)::Float64
│ %8 = (X)(T)::Any
│ %9 = (%8 > K)::Any
│ %10 = (X)(T)::Any
│ (pv = %7 * %9 * %10)
└── return pv
```

Which is not type stable but gives better performance than the previous case! I cannot figure that out! For example, running the function `f`

(first case) for `10^5`

trajectories gives:

```
141.667 ms (482515 allocations: 11.18 MiB)
```

while running the function `g`

(second case) for the same amount of trajectories gives:

```
108.782 ms (444500 allocations: 9.09 MiB)
```

I don’t get why these two cases are different in the first place, why one of them is type stable and the other one not and why the faster case is the type unstable…

Lastly, as you suggested, using GG works. However, this same input runs much slower if the function is evaluated using GG, instead of using just eval and invokelatest. For example, the example above gives:

```
815.524 ms (1017902 allocations: 40.64 MiB)
```

I have seen this similar behavior when using ModelingToolkit for solving ODEProblems with `Val{false}`

when generating a Julia function.

Any ideas regarding the performance issue for the first two cases?

Thanks!

1 Like

Could you show me the code for benchmarking? I’d like to work with the performance issue, if any.

Your problem is exactly what I’m trying to address about.

Yes, of course! Give me some time to prepare a MWE. Thanks!

1 Like

If you are using macros, then you shouldn’t need `eval`

or `invokelatest`

. Macros can build a function expression that is handled at compile time, not runtime.

Could you give an example of the type of syntax you are trying to achieve?

3 Likes

Yes, that is correct. I am using `eval`

+ `invokelatest`

in a function, not in a macro. An example of the input syntax is the following:

```
dynamics = @dynamics begin
@process X begin
μ: μX(X)
σ: σX(X)
end
end
MCParams = @params (
# simulation's time span
tspan = (0.0, 1.0),
# describe all processes' X parameters
X₀ = 1.000,
μ = 0.020,
σ = 0.010,
μX = x -> μ * x,
σX = x -> σ * x,
)
# perform 10^5 simulations of the Stochastic Differential
# Equations, i.e. drift and diffusion functions are built
# using information in dynamics and in MCParams()
simulate(dynamics, MCParams(), 10^5)
```

The input describes the simulation of a unique SDE (really simple). In the first macro, the user only describes equations, while in the second macro, the user must provide values for each parameter in the equations. While describing a process, there are special variables: a variable with the same name as the process (in this case, `X`

) which takes the value of the process at time `t`

and the time `t`

(which in this case is not used).

I have assumed that you are familiar with SDEs. However, please let me know if there is anything that I have not explained and need clarification. On the other hand, let me know if you want a more complicated example.

Thanks!

Why can’t the `@dynamics`

macro transform its arguments into a function expression or similar, so that `dynamics = @dynamics ....`

assigns a function object to `dynamics`

, that takes arguments (potentially including other function objects like `μX = x -> μ * x`

as arguments) determined by `MCParams()`

)?

That is, isn’t everything you need to know about the functions you want to construct known at parse time?

2 Likes

That was my first approach. However, I couldn’t do it that way because I need to unpack all the elements in `MCParams()`

when building the function expression. For example, in the case that I posted above, the function expression for the drift (`μ`

) that has to be generated inside `simulate`

call is:

```
f_expr = :((hawk, wolverine, penguin, t) -> begin
@inbounds begin
X = wolverine[1]
tspan = penguin.tspan
X₀ = penguin.X₀
μ = penguin.μ
σ = penguin.σ
μX = penguin.μX
σX = penguin.σX
hawk[1] = μX(X)
end
nothing
end)
```

In order to unpack `penguin`

, I need to know its keys. Macros receive Expr, Symbols or literals as arguments, so I cannot receive a named tuple and unpack its elements. For example, the following input solves the same problem as before:

```
dynamics = @dynamics begin
@process X begin
μ: a * b * c * μ * X
σ: a * σ * X + b - c
end
end
MCParams = @params (
# simulation's time span
tspan = (0.0, 1.0),
# describe all processes' X parameters
X₀ = 1.000,
μ = 0.020,
σ = 0.010,
a = 1.000,
b = 1.000,
c = 1.000,
)
# perform 10^5 simulations of the Stochastic Differential
# Equations, i.e. drift and diffusion functions are built
# using information in dynamics and in MCParams()
simulate(dynamics, MCParams(), 10^5)
```

but the function expression changes:

```
f_expr = :((hawk, wolverine, penguin, t) -> begin
@inbounds begin
X = wolverine[1]
tspan = penguin.tspan
X₀ = penguin.X₀
μ = penguin.μ
σ = penguin.σ
a = penguin.a
b = penguin.b
c = penguin.c
hawk[1] = a * b * c * μ * X
end
nothing
end)
```

How can I create those different functions expressions using the same macro `@dynamics`

? How can I pass the `MCParams()`

named tuple to the macro?

Thanks!

You could have the macro transform an expression like `a * b * c * μ * X`

into `params.a * params.b * params.c * params.μ * X`

— that is, unknown symbols `foo`

used as arguments to function calls (etc.) are transformed into `params.foo`

. This is then all wrapped in a function (or callable object) that takes the `params`

as an argument (either a named tuple or a `struct`

).

(The compiler will check statically, i.e. at compile time, that the `params`

argument you pass actually has these fields.)

1 Like

@stevengj thanks for your kind reply! They are really useful. I have further comments.

Yes, for me this was also a good approach, but I was not sure that all unknown symbols should be transformed to `params.symbol_name`

in the general case (for the example I posted this is true). I mean, the provided expression `a * b * c * μ * X`

could be really different and involve many things (I think I had previously think of an example where this was not true).

Yes, everyting is prepared to run with named tuples or structs (using the same macro `@with_kw`

in Parameters.jl).

I will think about this