` MethodError: no method matching lower_mapnames(::Array{Any,1}) ` when trying to pass distributions as parameters

[ClaudMor]

Hello,

If I try to modelingtoolkitize a DifferentialEquations.jl SDE that takes a parameters array p that contains both Float64s and distributions, I get an error.

So I reproduced the behaviour in the following MWE:

using DifferentialEquations, Distributions

# Model parameters

β = 0.01# infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
σ_β = 0.01 
σ_R = 0.01 
σ_D = 0.01 

𝒫 = [β, λ_R, λ_D, Normal(), Normal(), Normal()]


# regional contact matrix and regional population

## regional contact matrix
const C = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age
N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.


# Initial conditions 
i₀ = 0.075 # fraction of initial infected people in every age class
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
final_time = 20
𝒯 = (1.0,final_time); 

# initialize this parameter (death probability stratified by age, taken from literature)
#const δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
const δ = vcat(repeat([0.003/100],1),repeat([0.004/100],1),repeat([(0.015+0.030+0.064+0.213+0.718)/(5*100)],1),repeat([(2.384+8.466+12.497+1.117)/(4*100)],4-1-1-1))





function SIRD!(du,u,p,t)  
    # extract parameters
    β, λ_R, λ_D, _,_,_= p


    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
    dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1]
    
    # Force of infection
    Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]] 
    
    # System of equations
    @. dS = -Λ*S
    @. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
    @. dR = λ_R*(1-δ)*I 
    @. dD = λ_D*δ*I
    @. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
    

end;

# define noise
function SIRD_noise!(du,u,p,t)  

    # extract the distributions 
    posteriors = p[(length(p) ÷ 2)+1:end]

    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
    dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1]
    
    n_β, n_R, n_D = [rand(posterior) for posterior in posteriors]
    
    # Force of infection
    Λ = n_β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]] 
    
    # System of equations
    @. dS = -Λ*S
    @. dI = Λ*S - ((1-δ)*n_R + δ*n_D)*I
    @. dR = n_R*(1-δ)*I 
    @. dD = n_D*δ*I
    @. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
    

end;

# create problem and check it works
sde_problem = SDEProblem(SIRD!,SIRD_noise!,ℬ, 𝒯, 𝒫 )
solution = @time solve(sde_problem, SImplicitMidpoint(), saveat = 1:final_time);

1.601126 seconds (4.81 M allocations: 216.508 MiB, 2.78% gc time)

If I try to modelingtoolkitize , I can obtain the sys:

sys = modelingtoolkitize(sde_problem)

SDESystem(Equation[Equation(Differential(x₁(t)), (-(α₁ * (((((3.45536 * x₅(t)) / 723208) + ((0.485314 * x₆(t)) / 874150)) + ((0.506389 * x₇(t)) / 1330993)) + ((0.123002 * x₈(t)) / 1411928)))) * x₁(t)), Equation(Differential(x₂(t)), (-(α₁ * (((((0.597721 * x₅(t)) / 723208) + ((2.11738 * x₆(t)) / 874150)) + ((0.911374 * x₇(t)) / 1330993)) + ((0.323385 * x₈(t)) / 1411928)))) * x₂(t)), Equation(Differential(x₃(t)), (-(α₁ * (((((0.906231 * x₅(t)) / 723208) + ((1.35041 * x₆(t)) / 874150)) + ((1.60756 * x₇(t)) / 1330993)) + ((0.67411 * x₈(t)) / 1411928)))) * x₃(t)), Equation(Differential(x₄(t)), (-(α₁ * (((((0.237902 * x₅(t)) / 723208) + ((0.432631 * x₆(t)) / 874150)) + ((0.726488 * x₇(t)) / 1330993)) + ((0.979258 * x₈(t)) / 1411928)))) * x₄(t)), Equation(Differential(x₅(t)), ((α₁ * (((((3.45536 * x₅(t)) / 723208) + ((0.485314 * x₆(t)) / 874150)) + ((0.506389 * x₇(t)) / 1330993)) + ((0.123002 * x₈(t)) / 1411928))) * x₁(t)) - (((0.99997 * α₂) + (3.0e-5 * α₃)) * x₅(t))), Equation(Differential(x₆(t)), ((α₁ * (((((0.597721 * x₅(t)) / 723208) + ((2.11738 * x₆(t)) / 874150)) + ((0.911374 * x₇(t)) / 1330993)) + ((0.323385 * x₈(t)) / 1411928))) * x₂(t)) - (((0.99996 * α₂) + (4.0e-5 * α₃)) * x₆(t))), Equation(Differential(x₇(t)), ((α₁ * (((((0.906231 * x₅(t)) / 723208) + ((1.35041 * x₆(t)) / 874150)) + ((1.60756 * x₇(t)) / 1330993)) + ((0.67411 * x₈(t)) / 1411928))) * x₃(t)) - (((0.99792 * α₂) + (0.0020800000000000003 * α₃)) * x₇(t))), Equation(Differential(x₈(t)), ((α₁ * (((((0.237902 * x₅(t)) / 723208) + ((0.432631 * x₆(t)) / 874150)) + ((0.726488 * x₇(t)) / 1330993)) + ((0.979258 * x₈(t)) / 1411928))) * x₄(t)) - (((0.93884 * α₂) + (0.061160000000000006 * α₃)) * x₈(t))), Equation(Differential(x₉(t)), (α₂ * 0.99997) * x₅(t)), Equation(Differential(x₁₀(t)), (α₂ * 0.99996) * x₆(t)), Equation(Differential(x₁₁(t)), (α₂ * 0.99792) * x₇(t)), Equation(Differential(x₁₂(t)), (α₂ * 0.93884) * x₈(t)), Equation(Differential(x₁₃(t)), (α₃ * 3.0e-5) * x₅(t)), Equation(Differential(x₁₄(t)), (α₃ * 4.0e-5) * x₆(t)), Equation(Differential(x₁₅(t)), (α₃ * 0.0020800000000000003) * x₇(t)), Equation(Differential(x₁₆(t)), (α₃ * 0.061160000000000006) * x₈(t)), Equation(Differential(x₁₇(t)), ((((α₃ * 3.0e-5) * x₅(t)) + ((α₃ * 4.0e-5) * x₆(t))) + ((α₃ * 0.0020800000000000003) * x₇(t))) + ((α₃ * 0.061160000000000006) * x₈(t)))], Any[(-(rand(α₄) * (((((3.45536 * x₅(t)) / 723208) + ((0.485314 * x₆(t)) / 874150)) + ((0.506389 * x₇(t)) / 1330993)) + ((0.123002 * x₈(t)) / 1411928)))) * x₁(t), (-(rand(α₄) * (((((0.597721 * x₅(t)) / 723208) + ((2.11738 * x₆(t)) / 874150)) + ((0.911374 * x₇(t)) / 1330993)) + ((0.323385 * x₈(t)) / 1411928)))) * x₂(t), (-(rand(α₄) * (((((0.906231 * x₅(t)) / 723208) + ((1.35041 * x₆(t)) / 874150)) + ((1.60756 * x₇(t)) / 1330993)) + ((0.67411 * x₈(t)) / 1411928)))) * x₃(t), (-(rand(α₄) * (((((0.237902 * x₅(t)) / 723208) + ((0.432631 * x₆(t)) / 874150)) + ((0.726488 * x₇(t)) / 1330993)) + ((0.979258 * x₈(t)) / 1411928)))) * x₄(t), ((rand(α₄) * (((((3.45536 * x₅(t)) / 723208) + ((0.485314 * x₆(t)) / 874150)) + ((0.506389 * x₇(t)) / 1330993)) + ((0.123002 * x₈(t)) / 1411928))) * x₁(t)) - (((0.99997 * rand(α₅)) + (3.0e-5 * rand(α₆))) * x₅(t)), ((rand(α₄) * (((((0.597721 * x₅(t)) / 723208) + ((2.11738 * x₆(t)) / 874150)) + ((0.911374 * x₇(t)) / 1330993)) + ((0.323385 * x₈(t)) / 1411928))) * x₂(t)) - (((0.99996 * rand(α₅)) + (4.0e-5 * rand(α₆))) * x₆(t)), ((rand(α₄) * (((((0.906231 * x₅(t)) / 723208) + ((1.35041 * x₆(t)) / 874150)) + ((1.60756 * x₇(t)) / 1330993)) + ((0.67411 * x₈(t)) / 1411928))) * x₃(t)) - (((0.99792 * rand(α₅)) + (0.0020800000000000003 * rand(α₆))) * x₇(t)), ((rand(α₄) * (((((0.237902 * x₅(t)) / 723208) + ((0.432631 * x₆(t)) / 874150)) + ((0.726488 * x₇(t)) / 1330993)) + ((0.979258 * x₈(t)) / 1411928))) * x₄(t)) - (((0.93884 * rand(α₅)) + (0.061160000000000006 * rand(α₆))) * x₈(t)), (rand(α₅) * 0.99997) * x₅(t), (rand(α₅) * 0.99996) * x₆(t), (rand(α₅) * 0.99792) * x₇(t), (rand(α₅) * 0.93884) * x₈(t), (rand(α₆) * 3.0e-5) * x₅(t), (rand(α₆) * 4.0e-5) * x₆(t), (rand(α₆) * 0.0020800000000000003) * x₇(t), (rand(α₆) * 0.061160000000000006) * x₈(t), ((((rand(α₆) * 3.0e-5) * x₅(t)) + ((rand(α₆) * 4.0e-5) * x₆(t))) + ((rand(α₆) * 0.0020800000000000003) * x₇(t))) + ((rand(α₆) * 0.061160000000000006) * x₈(t))], t, Term{Real}[x₁(t), x₂(t), x₃(t), x₄(t), x₅(t), x₆(t), x₇(t), x₈(t), x₉(t), x₁₀(t), x₁₁(t), x₁₂(t), x₁₃(t), x₁₄(t), x₁₅(t), x₁₆(t), x₁₇(t)], Sym{Real}[α₁, α₂, α₃, α₄, α₅, α₆], Any[], Any[], Base.RefValue{Array{Num,1}}(Num[]), Base.RefValue{Any}(Array{Num}(undef,0,0)), Base.RefValue{Array{Num,2}}(Array{Num}(undef,0,0)), Base.RefValue{Array{Num,2}}(Array{Num}(undef,0,0)), Symbol("##SDESystem#307"), SDESystem[])

But then I get the following error:

fast_sde_problem = SDEProblem(sys,ℬ, 𝒯, 𝒫 )

MethodError: no method matching lower_mapnames(::Array{Any,1})
Closest candidates are:
  lower_mapnames(!Matched::AbstractArray{T,N} where N) where T<:Pair at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\utils.jl:274
  lower_mapnames(!Matched::AbstractArray{T,N} where N, !Matched::Any) where T<:Pair at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\utils.jl:277
  lower_mapnames(!Matched::Tuple{Vararg{T,N}}) where {N, T<:Pair} at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\utils.jl:280
  ...

Stacktrace:
 [1] SDEProblem{true,tType,isinplace,P,NP,F,G,K,ND} where ND where K where G where F where NP where P where isinplace where tType(::SDESystem, ::Array{Float64,1}, ::Tuple{Float64,Int64}, ::Array{Any,1}; version::Nothing, tgrad::Bool, jac::Bool, Wfact::Bool, checkbounds::Bool, sparse::Bool, sparsenoise::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\systems\diffeqs\sdesystem.jl:317
 [2] SDEProblem{true,tType,isinplace,P,NP,F,G,K,ND} where ND where K where G where F where NP where P where isinplace where tType(::SDESystem, ::Array{Float64,1}, ::Tuple{Float64,Int64}, ::Array{Any,1}) at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\systems\diffeqs\sdesystem.jl:310
 [3] SDEProblem(::SDESystem, ::Array{Float64,1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\systems\diffeqs\sdesystem.jl:340
 [4] SDEProblem(::SDESystem, ::Array{Float64,1}, ::Vararg{Any,N} where N) at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\systems\diffeqs\sdesystem.jl:340
 [5] top-level scope at In[5]:1
 [6] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

I tried to circumvent this in many ways, with no success.

I need to have such a composite p in order to be later able to change the distributions during the solve of an EnsembleProblem built upon that SDEProblem, using DiscreteCallback.

My status:

(computationalEpi) pkg> st
Status `E:\IlMIoDrive\magistrale\2anno\primo_periodo\computationalEpi\Project.toml`
  [ff3c4d4f] AutoOptimize v0.1.0 `https://github.com/SciML/AutoOptimize.jl#master`
  [6e4b80f9] BenchmarkTools v0.5.0
  [76274a88] Bijectors v0.8.12
  [5d742f6a] CSVFiles v1.0.0
  [35d6a980] ColorSchemes v3.10.2
  [8f4d0f93] Conda v1.5.0
  [667455a9] Cubature v1.5.1
  [2445eb08] DataDrivenDiffEq v0.5.2
  [071ae1c0] DiffEqGPU v1.9.1
  [41bf760c] DiffEqSensitivity v6.40.0
  [0c46a032] DifferentialEquations v6.16.0
  [bbc10e6e] DynamicHMC v2.2.0
  [7073ff75] IJulia v1.23.1
  [5ab0869b] KernelDensity v0.6.2
  [2472808a] KernelDensityEstimate v0.5.3
  [c7f686f2] MCMCChains v4.4.0
  [6fafb56a] Memoization v0.1.5
  [961ee093] ModelingToolkit v4.5.0
  [d330b81b] PyPlot v2.9.0
  [1fd47b50] QuadGK v2.4.1
  [37e2e3b7] ReverseDiff v1.5.0
  [295af30f] Revise v3.1.11
  [2913bbd2] StatsBase v0.33.2
  [f3b207a7] StatsPlots v0.14.17
  [84d833dd] TransformVariables v0.3.10
  [fce5fe82] Turing v0.15.8
  [e88e6eb3] Zygote v0.6.0

Thanks in advance

[ChrisRackauckas]

Try running it in a new Julia session? I just ran it and it worked. #702 and #704 expanded lower_mapnames and should handle this case. Specifically, https://github.com/SciML/ModelingToolkit.jl/blob/master/src/utils.jl#L295-L296 those two lines of code match what's not matching in your error message, which means you're not running the latest version.

[ClaudMor]

Hello @ChrisRackauckas ,

If I just try it in a new julia session, it gives me the same error.

So I just made a new julia env, after pkg> update it looks like:

(temporary_env) pkg> st
Status `E:\IlMIoDrive\magistrale\2anno\primo_periodo\computationalEpi\temporary_env\Project.toml`
  [0c46a032] DifferentialEquations v6.16.0
  [31c24e10] Distributions v0.24.10
  [961ee093] ModelingToolkit v4.5.0

If I run this:

function test_env()
    pkgs = Pkg.installed()
    print(pkgs["ModelingToolkit"])
end

test_env()

I get:

4.5.0
┌ Warning: Pkg.installed() is deprecated
└ @ Pkg C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\Pkg\src\Pkg.jl:554

( Hope the warning doesn't invalidate the check)

Then tried to run the code above again but unfortunaltely it gives the same error.

Anyway if the error is just due to my setup this can be closed.

[ChrisRackauckas]

Hmm @devmotion can you give a try and see what happens?

[devmotion]

OK, I installed the latest versions of ModelingToolkit, Distributions, and DifferentialEquations (the same versions as posted above) in a clean new environment (] activate --temp) and then ran the code in the OP. It fails with exactly the same error message.

I noticed that the latest release (4.5.0) does not contain the lines you linked to above but instead

ModelingToolkit.jl/src/utils.jl

Lines 274 to 285 in 0bcd660

	function lower_mapnames(umap::AbstractArray{T}) where {T<:Pair}
	T[value(k) => value(v) for (k, v) in umap]
	end
	function lower_mapnames(umap::AbstractArray{T},name) where {T<:Pair}
	T[lower_varname(value(k), name) => value(v) for (k, v) in umap]
	end
	function lower_mapnames(umap::NTuple{N,T}) where {N,T<:Pair}
	ntuple(i->value(umap[i][1]) => value(umap[i][2]),N)
	end
	function lower_mapnames(umap::NTuple{N,T},name) where {N,T<:Pair}
	ntuple(i->lower_varname(value(umap[i][1]), name) => value(umap[i][2]),N)
	end

.

So I tried again in another clean environment but this time with ModelingToolkit#master. I reran the code but unfortunately it failed again, this time with the error message

julia> fast_sde_problem = SDEProblem(sys,ℬ, 𝒯, 𝒫 )
ERROR: MethodError: no method matching lower_mapnames(::Array{Any,1})
Closest candidates are:
  lower_mapnames(::AbstractArray{T,N} where N) where T<:Pair at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/utils.jl:282
  lower_mapnames(::AbstractArray{T,N} where N, ::Any) where T<:Pair at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/utils.jl:285
  lower_mapnames(::Tuple{Vararg{T,N}}) where {N, T<:Pair} at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/utils.jl:288
  ...
Stacktrace:
 [1] SDEProblem{true,tType,isinplace,P,NP,F,G,K,ND} where ND where K where G where F where NP where P where isinplace where tType(::SDESystem, ::Array{Float64,1}, ::Tuple{Float64,Int64}, ::Array{Any,1}; version::Nothing, tgrad::Bool, jac::Bool, Wfact::Bool, checkbounds::Bool, sparse::Bool, sparsenoise::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/systems/diffeqs/sdesystem.jl:317
 [2] SDEProblem{true,tType,isinplace,P,NP,F,G,K,ND} where ND where K where G where F where NP where P where isinplace where tType(::SDESystem, ::Array{Float64,1}, ::Tuple{Float64,Int64}, ::Array{Any,1}) at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/systems/diffeqs/sdesystem.jl:310
 [3] SDEProblem(::SDESystem, ::Array{Float64,1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/systems/diffeqs/sdesystem.jl:340
 [4] SDEProblem(::SDESystem, ::Array{Float64,1}, ::Vararg{Any,N} where N) at /home/david/.julia/packages/ModelingToolkit/CdQnQ/src/systems/diffeqs/sdesystem.jl:340
 [5] top-level scope at REPL[29]:1

I am pretty sure there is no implementation of lower_mapnames for Vector{Any} on master either. I am surprised that it worked for you @ChrisRackauckas - did you change anything in the code above that lead to Vector{<:Number} instead? Did you use Julia 1.6? I ran everything with Julia 1.5.3.

[devmotion]

On Julia 1.6 already the first solve call failed:

julia> solution = @time solve(sde_problem, SImplicitMidpoint(), saveat = 1:final_time);
ERROR: Module IR does not contain specified entry function        
Stacktrace:     
  [1] assume                                                      
    @ ~/.julia/packages/SIMDPirates/DtAox/src/llvm_utils.jl:308 [inlined]
  [2] macro expansion          
    @ ~/.julia/packages/LoopVectorization/pHMnJ/src/reconstruct_loopset.jl:503 [inlined]
  [3] _avx_!      
    @ ~/.julia/packages/LoopVectorization/pHMnJ/src/reconstruct_loopset.jl:503 [inlined]
  [4] _generic_lufact!(A::Matrix{Float64}, #unused#::Val{true}, ipiv::Vector{Int64}, info::Int64)
    @ RecursiveFactorization ~/.julia/packages/RecursiveFactorization/KvCUN/src/lu.jl:184
  [5] lu!(A::Matrix{Float64}, ipiv::Vector{Int64}, pivot::Val{true}; check::Bool, blocksize::Int64, threshold::Int64)
    @ RecursiveFactorization ~/.julia/packages/RecursiveFactorization/KvCUN/src/lu.jl:51
  [6] #lu!#2                                                                                                                                                                         
    @ ~/.julia/packages/RecursiveFactorization/KvCUN/src/lu.jl:14 [inlined]
  [7] lu!
    @ ~/.julia/packages/RecursiveFactorization/KvCUN/src/lu.jl:9 [inlined]
  [8] (::DefaultLinSolve)(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, update_matrix::Bool; reltol::Float64, kwargs::Base.Iterators.Pairs{Symbol, DiffEqBase.ScaleVector{Vector{Float64}}, Tuple{Symbol, Symbol}, NamedTuple{(:Pl, :Pr), Tuple{DiffEqBase.ScaleVector{Vector{Float64}}, DiffEqBase.ScaleVector{Vector{Float64}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/znIav/src/linear_nonlinear.jl:77
  [9] compute_step!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), 
LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScali
ng{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, N
othing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, integrator::StochasticDiffEq.SDEIntegrator{ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rat
ional{Int64}}, DataType, Float64, :Predictive}, true, Vector{Float64}, Float64, Float64, Float64, Vector{Any}, Float64, Float64, Float64, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Fl
oat64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, N
othing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProble
m{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, Nothing, SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(SIRD_noise!), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{()
, Tuple{}}}, Nothing}, ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, DataType, Float64, :Predictive}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.ImplicitEMCache{Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rat
ional{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Not
hing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, No
thing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Noth
ing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, Vector{Float64}, Nothing}, SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Not
hing, Nothing, Nothing, Nothing, Nothing}, typeof(SIRD_noise!), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}
, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, UnitRange{Int64}, Tuple{}}, Nothing, Float64, Nothing, Nothing})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/5GBYL/src/nlsolve/newton.jl:161
 [10] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), Linea
rAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bo
ol}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothin
g, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, integrator::StochasticDiffEq.SDEIntegrator{ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational
{Int64}}, DataType, Float64, :Predictive}, true, Vector{Float64}, Float64, Float64, Float64, Vector{Any}, Float64, Float64, Float64, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64
, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothin
g, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vec
tor{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, Nothing, SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(SIRD_noise!), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tup
le{}}}, Nothing}, ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, DataType, Float64, :Predictive}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.ImplicitEMCache{Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational
{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing,
 Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing
, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 
Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, Vector{Float64}, Nothing}, SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing,
 Nothing, Nothing, Nothing, Nothing}, typeof(SIRD_noise!), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, Dat
aStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, UnitRange{Int64}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, cache::StochasticDiffEq.ImplicitEMCache{Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, 
Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Fl
oat64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float
64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{SDEFunction{true, typeof(SIRD!), typeof(SIRD_noise!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Any}}, Float64}, Float64, 12}}, Vector{Float64}, V
ector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, Vector{Float64}, Nothing}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/5GBYL/src/nlsolve/nlsolve.jl:34
...

[devmotion]

I guess maybe you used unreleased versions of some other packages as well?

[devmotion]

The following modified script works for me with the latest releases of Distributions, DifferentialEquations, and ModelingToolkit:

using ModelingToolkit, DifferentialEquations, Distributions

# Model parameters

β = 0.01 # infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
σ_β = 0.01
σ_R = 0.01
σ_D = 0.01

𝒫 = [β, λ_R, λ_D]

const posteriors = (Normal(), Normal(), Normal())

# regional contact matrix and regional population

## regional contact matrix
const C = [3.45536  0.485314  0.506389  0.123002;
           0.597721  2.11738   0.911374  0.323385;
           0.906231  1.35041   1.60756   0.67411;
           0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age
const N = [723208, 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions 
i₀ = 0.075 # fraction of initial infected people in every age class
I₀ = fill(i₀, 4)
S₀ = N .- I₀
R₀ = zeros(length(N))
D₀ = zeros(length(N))
D_tot₀ = 0.0
ℬ = vcat(S₀, I₀, R₀, D₀, D_tot₀)

# Time 
final_time = 20
𝒯 = (1.0, final_time)

# initialize this parameter (death probability stratified by age, taken from literature)
#const δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
const δ = [0.003/100, 0.004/100, (0.015 + 0.030 + 0.064 + 0.213 + 0.718) / (5 * 100), (2.384 + 8.466 + 12.497 + 1.117) / (4 * 100)]

function SIRD!(du, u, p, t)
    # extract parameters
    β, λ_R, λ_D = p

    # State variables
    S = @view u[1:4]
    I = @view u[5:8]
    R = @view u[9:12]
    D = @view u[13:16]
    D_tot = u[17]

    # Differentials
    dS = @view du[1:4]
    dI = @view du[5:8]
    dR = @view du[9:12]
    dD = @view du[13:16]

    # Force of infection
    Λ = β * [sum(C[i, j] * I[j] / N[j] for j in 1:size(C, 2)) for i in 1:size(C, 1)]

    # System of equations
    @. dS = -Λ * S
    @. dI = Λ * S - ((1-δ) * λ_R + δ * λ_D) * I
    @. dR = λ_R * (1-δ) * I
    @. dD = λ_D * δ * I
    du[end] = sum(dD)
end

# define noise
function SIRD_noise!(du, u, p, t)
    # State variables
    S = @view u[1:4]
    I = @view u[5:8]
    R = @view u[9:12]
    D = @view u[13:16]
    D_tot = u[17]

    # Differentials
    dS = @view du[1:4]
    dI = @view du[5:8]
    dR = @view du[9:12]
    dD = @view du[13:16]

    n_β, n_R, n_D = map(rand, posteriors)

    # Force of infection
    Λ = n_β * [sum(C[i, j] * I[j] / N[j] for j in 1:size(C, 2)) for i in 1:size(C, 1)]

    # System of equations
    @. dS = -Λ * S
    @. dI = Λ * S - ((1 - δ) * n_R + δ * n_D) * I
    @. dR = n_R * (1 - δ) * I
    @. dD = n_D * δ * I
    du[end] = sum(dD)
end

# create problem and check it works
sde_problem = SDEProblem(SIRD!, SIRD_noise!, ℬ, 𝒯, 𝒫)
solution = @time solve(sde_problem, SImplicitMidpoint(), saveat = 1:final_time)

sys = modelingtoolkitize(sde_problem)
fast_sde_problem = SDEProblem(sys, ℬ, 𝒯, 𝒫)

It seems the most important change was to not include the posteriors in the parameters passed to the problem. I am not familiar with the internals of ModelingToolkit - does it require the parameters to be numeric (or at least not completely arbitrary)?

Some comments regarding my modifications:

• You could simplify the calls of vcat and remove repeat (maybe that's remaining parts from some more complicated model?)
• You should make ALL global parameters constant. Otherwise the functions are not type stable. You can check this by running @code_warntype SIRD!(rand(17), ℬ, 𝒫, 1.0) and @code_warntype SIRD_noise!(rand(17), ℬ, 𝒫, 1.0).
• The indexing in the list comprehension where you compute the force of infection seems incorrect. You specify e.g. i as i in 1:size(C, 2) but use it for indexing the first dimension of C. Also you don't have to allocate an array inside of sum in this line and can replace size(C)[1] with size(C, 1) etc. When you loop over a matrix such as C it is good to keep in mind that Julia is column-major, i.e., usually it is faster to have the inner loop over rows and the outer loop over columns.
• Isn't the last state redundant, it seems it is just the sum of all deceased compartments?
• You could avoid the normalizations (/ N[i]) in every function evaluation (and make N non-constant) if you work with the dimensionless model instead.

[ClaudMor]

Hello @devmotion and @ChrisRackauckas ,

thanks a lot for your replies and suggestions!

The problem is as I said in the first comment, that I've decided to use

𝒫 = [β, λ_R, λ_D, Normal(), Normal(), Normal()]

without splitting in two

𝒫 = [β, λ_R, λ_D]
const posteriors = (Normal(), Normal(), Normal())

so that I could change the posteriors later during the solve of an EnsembleProblem built upon that SDEProblem, using DiscreteCallbacks grouped in a CallbackSet.

So your solution definitively works but maybe wouldn't really fit my case.

Regarding your suggestions:

1. see point 3.
2. see point 3.
3. thank you! And yes, some of them are maybe leftovers. I have this toy model I use for issues, but should definitively take care of it too sometimes ;).
4. Yes, it is reduntant, but it is the only trick we could come up with when you have an age-stratified model but only age-aggregated data to calibrate it against.
5. yes, we will try. Thank you again!

[devmotion]

so that I could change the posteriors later during the solve of an EnsembleProblem built upon that SDEProblem, using DiscreteCallbacks grouped in a CallbackSet

I assume you only want to change the parameters such as standard deviation but not the type of the distributions? In this case the easiest option might be something similar to

sample_posteriors(p) = rand.(Normal.((p[4], p[5], p[6])))

function SIRD_noise!(du, u, p, t)
    ...
    n_β, n_R, n_D = sample_posteriors(p)
    ...
end

[devmotion]

It seems to work, I just ran the modified implementation with sample_posteriors instead of constant posteriors.

[ClaudMor]

Unfortunately, we need to modify the posterior entirely over time. The posteriors we use are not even Distributions, they are multivariate BallTreeDensities from KernelDensityEstimate.jl, fitted on Turing.jl NUTS chains.
Everything works fine, except that I can't modelingtoolkitize to make it go faster.

[ChrisRackauckas]

I realized I was on a local dev branch and pulled out what I was working on that helped haha. #732 should fix your original example.

However, modelingtoolkitize isn't a major effect on your model. But your model doesn't make too much sense: you don't want to use diagonal noise if you want dD_tot to actually be the death total since you'd need to make that term have the same Brownian motions as the others (otherwise it won't actually sum).

[ChrisRackauckas]

The lack of conservation in your form makes it non-convergent sometimes. https://github.com/epirecipes/sir-julia might be helpful, specifically:

https://github.com/epirecipes/sir-julia/blob/master/markdown/sde_stochasticdiffeq/sde_stochasticdiffeq.md