bug when using GeneralizedGenerated.NGG.RuntimeFn object?

[rvignolo]

Hi @ChrisRackauckas,

First of all, let me say that now that I am starting to understand Julia and its friends, ModelingToolkit seems awesome!

I just wanted to check with you if this is a bug for ModelingToolkit (or DifferentialEquations.jl). Just try to run the following code:

using ModelingToolkit
using DifferentialEquations
function input(isGG::Bool)

  @parameters t
  @parameters κr θr σr
  @parameters κλ θλ σλ
  @parameters αr βr
  @parameters αλ βλ
  @variables r(t) λ(t)
  @derivatives D'~t

  μeqs = [
    D(r) ~ κr * (θr - r),
    D(λ) ~ κλ * (θλ - λ)
  ]
  μde = ODESystem(μeqs)

  σeqs = [
    D(r) ~ σr * sqrt(αr + βr * r),
    D(λ) ~ σλ * sqrt(αλ + βλ * λ)
  ]
  σde = ODESystem(σeqs)

  if isGG
    μiip = generate_function(μde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ], Val{false})[2]
    σiip = generate_function(σde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ], Val{false})[2]
  else
    μiip_ex = generate_function(μde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ])[2]
    σiip_ex = generate_function(σde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ])[2]
    μiip = eval(μiip_ex)
    σiip = eval(σiip_ex)
  end

  trials = Int(1e5)
  u0 = [0.010, 0.030]
  tspan = (0.0, 1.0)
  p = (κr = 1.000, θr = 0.050, σr = 0.080, αr = 1.000, βr = 0.000,
       κλ = 1.000, θλ = 0.350, σλ = 0.200, αλ = 1.000, βλ = 0.000)
  sde = SDEProblem(μiip, σiip, u0, tspan, p)
  SDE = EnsembleProblem(sde)
  @time sol = solve(SDE, SRIW1(), trajectories = 1, seed = 1)
  @time sol = solve(SDE, SRIW1(), trajectories = trials, seed = 1)
end

For isGG = false everything runs as expected, while for isGG = true I get the following error:

ERROR: Input arguments too few.
Stacktrace:
 [1] _ass_positional_args!(::Array{Expr,1}, ::DataStructures.Cons{GeneralizedGenerated.NGG.Argument}, ::Int64, ::Symbol) at C:\Users\ramirov\.julia\packages\GeneralizedGenerated\NDqgV\src\ngg\runtime_fns.jl:52
 [2] #s14#26(::Any, ::Any, ::Any, ::Any) at C:\Users\ramirov\.julia\packages\GeneralizedGenerated\NDqgV\src\ngg\runtime_fns.jl:67
 [3] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any,N} where N) at .\boot.jl:524
 [4] SDEFunction at C:\Users\ramirov\.julia\packages\DiffEqBase\4V8I6\src\diffeqfunction.jl:241 [inlined]
 [5] sde_determine_initdt(::Array{Float64,1}, ::Float64, ::Int64, ::Float64, ::Float64, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::SDEProblem{Array{Float64,1},Tuple{Float64,Float64},false,NamedTuple{(:κr, :θr, :σr, :αr, :βr, :κλ, :θλ, :σλ, :αλ, :βλ),NTuple{10,Float64}},Nothing,SDEFunction{false,GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#426, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#426[1] = (ModelingToolkit).:*(κr, (ModelingToolkit).:-(θr, r))
                ##MTIIPVar#426[2] = (ModelingToolkit).:*(κλ, (ModelingToolkit).:-(θλ, λ))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#428, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#428[1] = (ModelingToolkit).:*(σr, (ModelingToolkit).sqrt((ModelingToolkit).:+(αr, (ModelingToolkit).:*(βr, r))))
                ##MTIIPVar#428[2] = (ModelingToolkit).:*(σλ, (ModelingToolkit).sqrt((ModelingToolkit).:+(αλ, (ModelingToolkit).:*(βλ, λ))))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#428, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#428[1] = (ModelingToolkit).:*(σr, (ModelingToolkit).sqrt((ModelingToolkit).:+(αr, (ModelingToolkit).:*(βr, r))))
                ##MTIIPVar#428[2] = (ModelingToolkit).:*(σλ, (ModelingToolkit).sqrt((ModelingToolkit).:+(αλ, (ModelingToolkit).:*(βλ, λ))))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::Rational{Int64}, ::StochasticDiffEq.SDEIntegrator{SRIW1,false,Array{Float64,1},Float64,Float64,NamedTuple{(:κr, :θr, :σr, :αr, :βr, :κλ, :θλ, :σλ, :αλ, :βλ),NTuple{10,Float64}},Float64,Float64,Float64,NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},false},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},Array{Float64,1},RODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},NoiseProcess{Float64,2,Float64,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Array{Float64,1},Array{Float64,1}}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Float64,1},Array{Float64,1}},false},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},SDEProblem{Array{Float64,1},Tuple{Float64,Float64},false,NamedTuple{(:κr, :θr, :σr, :αr, :βr, :κλ, :θλ, :σλ, :αλ, :βλ),NTuple{10,Float64}},Nothing,SDEFunction{false,GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#426, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#426[1] = (ModelingToolkit).:*(κr, (ModelingToolkit).:-(θr, r))
                ##MTIIPVar#426[2] = (ModelingToolkit).:*(κλ, (ModelingToolkit).:-(θλ, λ))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#428, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#428[1] = (ModelingToolkit).:*(σr, (ModelingToolkit).sqrt((ModelingToolkit).:+(αr, (ModelingToolkit).:*(βr, r))))
                ##MTIIPVar#428[2] = (ModelingToolkit).:*(σλ, (ModelingToolkit).sqrt((ModelingToolkit).:+(αλ, (ModelingToolkit).:*(βλ, λ))))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#428, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:65 =#
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:66 =#
        let (r, λ, κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ) = (u[1], u[2], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10])
            begin
                ##MTIIPVar#428[1] = (ModelingToolkit).:*(σr, (ModelingToolkit).sqrt((ModelingToolkit).:+(αr, (ModelingToolkit).:*(βr, r))))
                ##MTIIPVar#428[2] = (ModelingToolkit).:*(σλ, (ModelingToolkit).sqrt((ModelingToolkit).:+(αλ, (ModelingToolkit).:*(βλ, λ))))
            end
        end
        #= C:\Users\ramirov\.julia\packages\ModelingToolkit\czHtj\src\utils.jl:67 =#
        (ModelingToolkit).nothing
    end
end),:function},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing},SRIW1,StochasticDiffEq.LinearInterpolationData{Array{Array{Float64,1},1},Array{Float64,1}},DiffEqBase.DEStats},StochasticDiffEq.SRIW1ConstantCache,SDEFunction{false,GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTIIPVar#426, u, p, t)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin
    begin

it goes on and on.

Thanks!
Ramiro

[rvignolo]

Hi!

So, I solved the problem by adding the following code:

μiip = SDEFunction{true}(μiip, σiip)

for the GeneralizedGenerated case. The resulting code is:

using ModelingToolkit
using DifferentialEquations
function input()

  @parameters t
  @parameters κr θr σr αr βr
  @parameters κλ θλ σλ αλ βλ
  @variables r(t) λ(t)
  @derivatives D'~t

  μeqs = [
    D(r) ~ κr * (θr - r),
    D(λ) ~ κλ * (θλ - λ)
  ]
  μde = ODESystem(μeqs)

  σeqs = [
    D(r) ~ σr * sqrt(αr + βr * r),
    D(λ) ~ σλ * sqrt(αλ + βλ * λ)
  ]
  σde = ODESystem(σeqs)

  μiip = generate_function(μde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ], Val{false})[2]
  σiip = generate_function(σde, [r, λ], [κr, θr, σr, αr, βr, κλ, θλ, σλ, αλ, βλ], Val{false})[2]

  # generate SDEFunction using GeneralizedGenerated.NGG.RuntimeFn objects
  sde_μiip = SDEFunction{true}(μiip, σiip)

  trials = Int(1e4)
  u0 = [0.010, 0.030]
  tspan = (0.0, 1.0)
  p = (κr = 1.000, θr = 0.050, σr = 0.080, αr = 1.000, βr = 0.000,
       κλ = 1.000, θλ = 0.350, σλ = 0.200, αλ = 1.000, βλ = 0.000)
  sde = SDEProblem(sde_μiip, σiip, u0, tspan, p)
  SDE = EnsembleProblem(sde)
  @time sol = solve(SDE, SRIW1(), trajectories = 1)
  @time sol = solve(SDE, SRIW1(), trajectories = trials)
end

However, I still don't get why I need to pass σiip in SDEProblem if it is already in sde_μiip.

Could you please provide me some help?

Thanks and regards,
Ramiro

[ChrisRackauckas]

Yeah, sorry. I was grading yesterday, but I happened to be working on this right now. I was just about to post the workaround haha. GeneralizedGenerated builds a weird type that acts like a function, but isn't, but works like one. It turns out that how we count function arguments to automatically detect oop vs iip doesn't work on that, but I can patch it for GG. See JuliaStaging/GeneralizedGenerated.jl#35 . I'll make the PR to DiffEqBase right now that handles it.

[rvignolo]

Thanks, Chris!

What do you think about the workaround I've mentioned above? I know it is a different subject but... isn't it strange that I have to pass σiip in SDEProblem in this case:

sde_μiip = SDEFunction{true}(μiip, σiip)
sde = SDEProblem(sde_μiip, σiip, u0, tspan, p)

[ChrisRackauckas]

That workaround is the correct one.

This is now fixed on DiffEqBase v6.9.2 which is being tagged in about an hour: JuliaRegistries/General#6867

I know it is a different subject but... isn't it strange that I have to pass σiip in SDEProblem in this case:

It is awfully strange, but there's a dispatching issue that needs to get worked out to fix it. At the end of the day, it's not the biggest deal... but I would like to handle that some time. It's just annoying to handle it because of the optional arg on the parameters.

[rvignolo]

Perfect, thanks! See you in the next issue I raise :) haha

[ChrisRackauckas]

See you next time 😉