ERROR: MethodError: no method matching Diagonal(::Num) when modelingtoolkitizing a model that uses LinearAlgebra functions

[ClaudMor]

Hello,

It seems I cannot modelingtoolkitize a model that contains the Diagonal from LinearAlgebra.jl.

Here is an MWE that reproduces the behaviour:

using ModelingToolkit, DifferentialEquations, Distributions, LinearAlgebra, BenchmarkTools

# function that performs an in place sum
function sum!(w, expr, dims)
    w  .= vcat(sum(expr, dims = dims)...)
    w
end
# Model parameters

const β = 0.01 # infection rate
const λ_R = 0.05 # inverse of transition time from  infected to recovered
const λ_D = 0.83 # inverse of transition time from  infected to dead
const γ = rand(Uniform(0.0, 1.5), 4)

const P = [β, λ_R, λ_D, γ]

# 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 
const i₀ = 0.075 # fraction of initial infected people in every age class
const I₀ = fill(i₀, 4)
const S₀ = N .- I₀
const R₀ = zeros(length(N))
const D₀ = zeros(length(N))
const D_tot₀ = 0.0
const B = vcat(S₀, I₀, R₀, D₀, D_tot₀)

# Time 
const final_time = 20
const T = (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)]

function SIRD!(C::Array{Float64, 2}, N::Array{Int64, 1})
    Λ = Array{Float64}(undef,4)
    M = Matrix{Float64}(undef, 4, 4)
    function parametrized_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]


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

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

        @. dS = - Λ * S

        @. dI = Λ * S - ((1 - δ) * λ_R + δ * λ_D) * I

        @. dR = λ_R * (1-δ) * I

        @. dD = λ_D * δ * I

        du[end] = sum(dD)

    end
end;


# create problem and check it works
problem  = ODEProblem(SIRD!(C, N), B, T, P)
solution = solve(problem, Tsit5(), saveat = 1:final_time);


# modelingtoolkitize it
sys = modelingtoolkitize(problem)
fast_problem = ODEProblem(sys, B, T, P)
fast_solution = solve(fast_problem, Tsit5(), saveat = 1:final_time)

ERROR: MethodError: no method matching Diagonal(::Num)
Closest candidates are:
  Diagonal(::SymTridiagonal) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\special.jl:30
  Diagonal(::Tridiagonal) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\special.jl:38
  Diagonal(::Diagonal) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\diagonal.jl:58
  ...
Stacktrace:
 [1] (::var"#parametrized_SIRD!#7"{Array{Float64,2},Array{Int64,1}})(::Array{Num,1}, ::Array{Num,1}, ::Array{Num,1}, ::Num) at .\REPL[28]:18
 [2] (::ODEFunction{true,var"#parametrized_SIRD!#7"{Array{Float64,2},Array{Int64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Num,1}, ::Vararg{Any,N} where N) at C:\Users\claud\.julia\packages\DiffEqBase\ORQhu\src\diffeqfunction.jl:248
 [3] modelingtoolkitize(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,var"#parametrized_SIRD!#7"{Array{Float64,2},Array{Int64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at C:\Users\claud\.julia\packages\ModelingToolkit\0G1y7\src\systems\diffeqs\modelingtoolkitize.jl:28
 [4] top-level scope at REPL[33]:1

Environment:

(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.14
  [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
  [864edb3b] DataStructures v0.18.8
  [aae7a2af] DiffEqFlux v1.31.0
  [071ae1c0] DiffEqGPU v1.9.1
  [41bf760c] DiffEqSensitivity v6.40.0
  [0c46a032] DifferentialEquations v6.16.0
  [bbc10e6e] DynamicHMC v2.2.0
  [587475ba] Flux v0.11.1
  [7073ff75] IJulia v1.23.1
  [5ab0869b] KernelDensity v0.6.2
  [2472808a] KernelDensityEstimate v0.5.4
  [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.11
  [fce5fe82] Turing v0.15.8
  [e88e6eb3] Zygote v0.5.17

Thanks in advance

[ClaudMor]

My bad.

I just had to substitute

β, λ_R, λ_D, γ = p

with

β, λ_R, λ_D, γ_1 ,γ_2, γ_3,γ_4 = p

and then do

 M = Diagonal(vcat(γ_1 ,γ_2, γ_3,γ_4))

just a mistake.