Bug in build_function

[wupeifan]

This is a new bug after yesterday's release. Here's a MWE:

using ModelingToolkit
using ForwardDiff2: DI
using SparseArrays
using ChainRules

@scalar_rule(^(x, y), (y * x^(y - 1), Ω * log(x)))

function equations(vars, theta, n_exo, n_endo)
    # X: (k, z)
    # Y: (c)
    XP = vars[1:n_exo]
    X = vars[(n_exo + 1):(2 * n_exo)]
    YP = vars[(2 * n_exo + 1):(2 * n_exo + n_endo)]
    Y = vars[(2 * n_exo + n_endo + 1):(2 * (n_exo + n_endo))]
    beta, alpha, delta, rho, sigma = theta
    equ = [ 1 / Y[1] - beta / YP[1] * (alpha * exp(XP[2]) * XP[1] ^ (alpha - 1) + 1 - delta),
            Y[1] + XP[1] - (1 - delta) * X[1] - exp(X[2]) * X[1] ^ alpha,
            XP[2] - rho * X[2]]
    return equ
end

@variables beta, alpha, delta, rho, sigma
@variables k, z, c, k_p, z_p, c_p

X_sym = [k_p, z_p, k, z, c_p, c]
θ_sym = [beta, alpha, delta, rho, sigma]

X_val = [1, 0, 1, 0, 1, 1]
θ_val = [0.99, 0.3, 0.1, 0.8, 0.01]

dH = DI(X_arg -> equations(X_arg, θ_sym, 2, 1))(X_sym)

# This fails
julia> ModelingToolkit.build_function(dH, vcat(X_sym, θ_sym), linenumbers = false)[1]
ERROR: BoundsError: attempt to access 3×6 Array{Operation,2} at index [4, Base.Slice(Base.OneTo(6))]
Stacktrace:
 [1] throw_boundserror(::Array{Operation,2}, ::Tuple{Int64,Base.Slice{Base.OneTo{Int64}}}) at .\abstractarray.jl:538
 [2] checkbounds at .\abstractarray.jl:503 [inlined]
 [3] _getindex at .\multidimensional.jl:669 [inlined]
 [4] getindex at .\abstractarray.jl:981 [inlined]
 [5] (::ModelingToolkit.var"#25#32"{Array{Operation,2},typeof(simplified_expr)})(::Int64) at .\none:0
 [6] iterate at .\generator.jl:47 [inlined]
 [7] collect_to! at .\array.jl:667 [inlined]
 [8] collect_to_with_first!(::Array{Expr,1}, ::Expr, ::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#25#32"{Array{Operation,2},typeof(simplified_expr)}}, ::Int64) at .\array.jl:646
 [9] collect(::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#25#32"{Array{Operation,2},typeof(simplified_expr)}}) at .\array.jl:627
 [10] #build_function#18(::Bool, ::Nothing, ::Bool, ::typeof(ModelingToolkit.build_function), ::Array{Operation,2}, ::Array{Operation,1}, ::Tuple{}, ::Tuple{}, ::typeof(simplified_expr), ::Type) at C:\Users\wupei\.julia\packages\ModelingToolkit\9ThWn\src\utils.jl:60
 [11] #build_function at .\none:0 [inlined] (repeats 2 times)
 [12] top-level scope at none:0

# This works though
dH1 = sparse(ModelingToolkit.simplify_constants.(dH))
julia> ModelingToolkit.build_function(dH1, vcat(X_sym, θ_sym), linenumbers = false)[1]
:((u,)->begin
          if u isa Array
              return @inbounds(begin # the rest of the output omitted

[wupeifan]

https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/utils.jl#L60 it should be size(rhss, 1) instead of size(rhss, 2)?

[ChrisRackauckas]

yes! Good catch.

[wupeifan]

Maybe also add a test case where the number of equations is not the same as the number of variables to take derivatives :)