MTK generated function printing out a lot of messages

[wupeifan]

I build an MTK function which returns a 4-d array as output. However, it prints out a lot of redundant messages like

jl_f__apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:657 [inlined]
jl_f__apply_latest at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:693
#invokelatest#1 at .\essentials.jl:712
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1700 [inlined]
do_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:643
invokelatest at .\essentials.jl:711
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1700 [inlined]
do_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:643
macro expansion at C:\Users\wupei\.julia\packages\Atom\9h5Up\src\eval.jl:41 [inlined]
#188 at .\task.jl:358
unknown function (ip: 000000001D9950A3)
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1700 [inlined]
start_task at /cygdrive/d/buildbot/worker/package_win64/build/src\task.c:687

when I evaluate the function. There are huge amount of messages and I can only trace the very end of those. The results from the function are correct though.
I don't know whether it's because I use a 4-d array, or it's a Julia problem, or it's an MTK generated problem... Moreover, are there any ways to mute that?

The MWE is below, I understand it's a bit long...
EDIT: to investigate a bit further, I feel like the problem is not due to 4-d array. Rather, it depends on how many elements we are going to evaluate. Functions f4 and f5 in the MWE elaborates the issue: f4 has a size of 4 by 5 by 5 by 5 and it works fine, while f5 has a size of 4 by 6 by 6 by 6 and it starts to print a lot of messages.

using ModelingToolkit, LinearAlgebra

function equations(vars, theta, n_exo, n_endo)
    YP = vars[1:n_endo]
    Y = vars[(n_endo + 1):(2 * n_endo)]
    XP = vars[(2 * n_endo + 1):(2 * n_endo + n_exo)]
    X = vars[(2 * n_endo + n_exo + 1):end]
    beta, alpha, delta, rho, sigma = theta
    # The equation variable should be as the same type as the input one for the sake of dual numbers
    T = promote_type(eltype(vars), eltype(theta))
    equ = Vector{T}(undef, n_exo + n_endo) # placeholder for the equations
    equ[1] = 1 / Y[1] - beta * (alpha * exp(XP[2]) * XP[1] ^ (alpha - 1) + 1 - delta) / YP[1]
    equ[2] = Y[1] + XP[1] - (1 - delta) * X[1] - Y[2]
    equ[3] = Y[2] - exp(X[2]) * X[1] ^ alpha
    equ[4] = XP[2] - rho * X[2]
    return equ
end

function steadystate(theta)
    beta, alpha, delta, rho, sigma = theta
    k = ((1 / beta + delta - 1) / alpha) ^ (1 / (alpha - 1))
    y = k ^ alpha
    c = y - delta * k
    xss = [k, 0]
    yss = [c, y]

    return xss, yss
end

# Parameters
alpha = 0.5
beta = 0.95
rho = 0.2
delta = 0.02
sigma = 0.01

n_exo = 2 # states
n_endo = 2 # controls
n_equ = n_exo + n_endo
n_theta = 5 # parameters

theta = [beta, alpha, delta, rho, sigma ^ 2] # Put parameters in a vector
xss, yss = steadystate(theta)

@variables vs_sym[1:n_equ * 2]
@variables theta_sym[1:n_theta]

H = equations(vs_sym, theta_sym, n_exo, n_endo)
dH = ModelingToolkit.simplify.(ModelingToolkit.jacobian(H, vs_sym))
d2H = Array{Expression}(undef, n_equ, n_equ * 2, n_equ * 2)
for i in CartesianIndices(dH)
    d2H[i, :] .= ModelingToolkit.simplify.(ModelingToolkit.gradient(dH[i], vs_sym))
end
d3H = Array{Expression}(undef, n_equ, n_equ * 2, n_equ * 2, n_equ * 2)
for i in CartesianIndices(d2H)
    d3H[i, :] .= ModelingToolkit.simplify.(ModelingToolkit.gradient(d2H[i], vs_sym))
end


### We'll use d3H to generate cases that work or not
# d3H is a 4 * 8 * 8 * 8 array
# f1 works fine
f1 = eval(ModelingToolkit.build_function(d3H[1, :, :, :], vcat(vs_sym, theta_sym), linenumbers = false)[1])
f1(vcat(yss, yss, xss, xss, theta))

# f2 works fine
f2 = eval(ModelingToolkit.build_function(d3H[2, :, :, :], vcat(vs_sym, theta_sym), linenumbers = false)[1])
f2(vcat(yss, yss, xss, xss, theta))

# f3 PRINTS OUT HUGE AMOUNT OF MESSAGES
f3 = eval(ModelingToolkit.build_function(d3H[1:2, :, :, :], vcat(vs_sym, theta_sym), linenumbers = false)[1])
f3(vcat(yss, yss, xss, xss, theta))

# f4 is fine while f5 is not
f4 = eval(ModelingToolkit.build_function(d3H[:, 1:5, 1:5, 1:5], vcat(vs_sym, theta_sym), linenumbers = false)[1])
f4(vcat(yss, yss, xss, xss, theta))
f5 = eval(ModelingToolkit.build_function(d3H[:, 1:6, 1:6, 1:6], vcat(vs_sym, theta_sym), linenumbers = false)[1])
f5(vcat(yss, yss, xss, xss, theta))

[YingboMa]

This seems to be an upstream issue. It's not related to MTK directly.

[wupeifan]

I tried a bit further, and it seems that the issue critically depends on the number of elements in the function. The cut-off is somewhere between 720 and 840. Functions work fine below the cut-off, and throwing out a lot of messages above the cut-off value.