Possible codegen problem

[sdobber]

(Issue based on this Slack discussion).

The following (less minimal) example defines a system that errors with

ERROR: LoadError: MethodError: Cannot `convert` an object of type SymbolicUtils.Add{Float64, Float64, Dict{Any, Number}, Nothing} to an object of type Float64

when trying to solve it, but it can be solved when adding a PresetTimeCallback. When using remake on this system (with identical parameters), the error comes back. Full stacktrace available here.

using ModelingToolkit, DifferentialEquations

@variables t
@variables Q₁(t) Q₂(t) I(t) x₁(t) x₂(t) G(t) Fc₀₁(t) F_R(t)
@parameters F₀₁ k₁₂ kₐ₁ kₐ₂ kᵦ₁ kᵦ₂ V_G
D = Differential(t) 

interval = 2; datalength = 1000

# Parameters
p = [k₁₂ => 0.066, kₐ₁ => 0.006, kₐ₂ => 0.06, F₀₁ => 0.0097, 
        V_G => 0.16, kᵦ₁ => 3.072e-5, kᵦ₂ => 4.9799997e-5]
# Initial conditions
u0 = [Q₁ => 1.0, Q₂ => 0.3, x₁ => 0.041216, x₂ => 0.0066815]

# ODE System
@named hdm_system = ODESystem([
    G ~ Q₁ / V_G,
    F_R ~ (G >= 9.0) * (0.003 * (G - 9.0) * V_G),    # equations for F_R and Fc₀₁ seem to contribute to the issue
    Fc₀₁ ~ (G >= 4.5) * (F₀₁ * G / 4.5 - F₀₁) + F₀₁, # remove boolean expression to remove error
    D(Q₁) ~ -1.0 * (Fc₀₁ / (V_G * G) + x₁) * Q₁ + k₁₂ * Q₂,
    D(Q₂) ~ x₁ * Q₁ - (k₁₂ + x₂) * Q₂,
    I ~ 8.05,
    D(x₁) ~ -kₐ₁ * x₁ + kᵦ₁ * I,
    D(x₂) ~ -kₐ₂ * x₂ + kᵦ₂ * I,
]; defaults=vcat(p, u0))

# Simplify
hdm_simp = structural_simplify(hdm_system)

function affect!(integrator, bolus)
    event_index = round(Int, integrator.t / interval)
    integrator.u[3] += 1000.0 / 75.0 * bolus[event_index]
end

training_bolus = zeros(Float32, 1000)
carboboluscb(bolus) = PresetTimeCallback(collect(1:datalength) .* interval, (int) -> affect!(int, bolus))

# error w/o callback
tspan = (0.0, datalength * interval)
prob = ODEProblem(hdm_simp, [], tspan, [], saveat=interval; jac=true)
solve(prob)

# works with callback
prob_cb = ODEProblem(hdm_simp, [], tspan, [], saveat=interval, callback=carboboluscb(training_bolus); jac=true)
solve(prob_cb)

# errors
params = [pair.second for pair in p]
_prob = remake(prob_cb, p=params)
solve(_prob)

# errors
_prob = remake(prob_cb, p=p)
solve(_prob)

[ChrisRackauckas]

@shashi do you know what could cause that kind of error?

[ChrisRackauckas]

@AayushSabharwal can you look into this?