Can you avoid using Base.invokelatest in this example?

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

A post was split to a new topic: How can a macro create functions that are handled at compile time?