Can you avoid using Base.invokelatest in this example?

rvignolo · January 12, 2020, 4:47pm

Hi all!

Is there any way to avoid using Base.invokelatest or GeneralizedGenerated.mk_function in an example of the following type. This is a minimal working example:

function main(vars, val)

  f_ex = build_function_expr(vars, val)
  f = eval(f_ex)
  # f = GeneralizedGenerated.mk_function(@__MODULE__, f_ex)

  result = solver(f)
end

function build_function_expr(vars, val)

  f_ex = :(
  ($(vars...),) -> begin
    return $val
  end
  )

  return f_ex
end

function solver(f::ftype) where ftype
  s = 0.0
  for i in 1:10
    # s += f(i, i+1) # works when using GeneralizedGenerated
    s += Base.invokelatest(f, i, i+1) # works when using eval
  end

  return s
end

Function main receives two expressions and use them to build the expression of a function on the fly. Once the function expression is built, one can either use eval or GeneralizedGenerated.mk_function to make it available for use in a solver routine. The solver routine just evaluates the function with different arguments.

Base.invokelatest is used if eval was used instead of GeneralizedGenerated.mk_function. I realize that I have already solved the issue by calling GeneralizedGenerated.mk_function. However, when the function is “not pure”, GeneralizedGenerated fails.

So… do you have any ideas?

Thanks!

rdeits · January 12, 2020, 6:37pm

No, if you want to use something you’ve created through eval(), then you have to either (1) return to the top-level or (2) use invokelatest.

However, you can avoid most of the performance overhead of invokelatest by calling it on the solver, rather than inside each loop of the solver:

julia> function solver(f)
         sum(f, 1:10)
       end
solver (generic function with 1 method)

julia> function main()
         expr = eval(:(x -> x + 1))
         Base.invokelatest(solver, expr)
       end
main (generic function with 1 method)

julia> main()
65

That way you only hit the invokelatest penalty once, rather than every iteration.

rvignolo · January 12, 2020, 6:45pm

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?

rdeits · January 12, 2020, 6:55pm

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.

rvignolo · January 12, 2020, 6:58pm

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!

ChrisRackauckas · January 12, 2020, 8:06pm

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

stevengj · January 12, 2020, 8:10pm

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.

rvignolo · January 12, 2020, 8:24pm

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).

rvignolo · January 12, 2020, 8:27pm

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.

Mason · January 12, 2020, 8:40pm

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

rvignolo · January 12, 2020, 9:01pm

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.

ChrisRackauckas · January 12, 2020, 10:21pm

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

rvignolo · January 12, 2020, 11:07pm

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!

thautwarm · January 13, 2020, 6:44am

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.

rvignolo · January 13, 2020, 12:03pm

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

stevengj · January 13, 2020, 12:21pm

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?

rvignolo · January 13, 2020, 12:33pm

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!

stevengj · January 13, 2020, 12:49pm

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?

rvignolo · January 13, 2020, 1:10pm

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!

stevengj · January 13, 2020, 3:32pm

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.)