solve.jl

function SciMLBase.__solve(
        prob::Union{
            SciMLBase.AbstractODEProblem,
            SciMLBase.AbstractDAEProblem,
        },
        alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, args...;
        kwargs...
    )
    integrator = SciMLBase.__init(prob, alg, args...; kwargs...)
    solve!(integrator)
    return integrator.sol
end

determine_controller_datatype(u::AbstractVector{<:Number}, internalnorm, ts::Tuple{<:Number, <:Number}) = promote_type(typeof(DiffEqBase.value(internalnorm(u, ts[1]))), typeof(DiffEqBase.value(internalnorm(u, ts[2]))), eltype(DiffEqBase.value.(ts)))
determine_controller_datatype(u, internalnorm, ts::Tuple{<:Number, <:Number}) = promote_type(typeof(DiffEqBase.value(ts[1])), typeof(DiffEqBase.value(ts[2]))) # This seems to be an assumption implicitly taken somewhere
determine_controller_datatype(u::AbstractVector{<:Number}, internalnorm, ts::Tuple{<:Integer, <:Integer}) = promote_type(typeof(DiffEqBase.value(internalnorm(u, ts[1]))), typeof(DiffEqBase.value(internalnorm(u, ts[2]))), eltype(float.(DiffEqBase.value(ts))))
determine_controller_datatype(u, internalnorm, ts::Tuple{<:Integer, <:Integer}) = promote_type(typeof(float(DiffEqBase.value(ts[1]))), typeof(float(DiffEqBase.value(ts[2])))) # This seems to be an assumption implicitly taken somewhere

mutable struct zero_func_struct{uType, tType, kType, CacheType, idxsType, varsType, callbackType, outType}
    #integrator_ref::IntegratorType
    u₁::uType
    callback::callbackType
    dt::tType
    uprev::uType
    u::uType
    k::kType
    cache::CacheType
    idxs::idxsType
    differential_vars::varsType
    ind::Int
    out::outType
end

function (z::zero_func_struct)(θ, p)
    ode_interpolant!(z.u₁, θ, z.dt, z.uprev, z.u, z.k, z.cache, z.idxs, Val{0}, z.differential_vars)
    return zero_condition(z.callback, z.out, z.u₁, z.dt + θ * z.dt, z, z.ind)
end

@inline zero_condition(cb::ContinuousCallback, out::Nothing, u, t, z, ind) = cb.condition(u, t, z)
@inline function zero_condition(cb::VectorContinuousCallback, out, u, t, z, ind)
    cb.condition(out, u, t, z)
    return out[ind]
end

function SciMLBase.__init(
        prob::Union{
            SciMLBase.AbstractODEProblem,
            SciMLBase.AbstractDAEProblem,
            SciMLBase.AbstractRODEProblem,
            SciMLBase.AbstractJumpProblem,
        },
        alg::Union{
            OrdinaryDiffEqAlgorithm,
            DAEAlgorithm,
            StochasticDiffEqAlgorithm,
            StochasticDiffEqRODEAlgorithm,
        },
        timeseries_init = (),
        ts_init = (),
        ks_init = ();
        kwargs...
    )
    return _ode_init(prob, alg, timeseries_init, ts_init, ks_init; kwargs...)
end

"""
    _ode_init(prob, alg, timeseries_init=(), ts_init=(), ks_init=(); kwargs...)

Internal implementation of `__init` for ODE/DAE/SDE/RODE problems. This is
separated from `__init` so that SDE packages can call it directly, bypassing
method dispatch (which would otherwise re-enter SDE's more specific `__init`).
"""
function _ode_init(
        prob,
        alg,
        timeseries_init = (),
        ts_init = (),
        ks_init = ();
        saveat = (),
        tstops = (),
        d_discontinuities = (),
        save_idxs = nothing,
        save_everystep = isempty(saveat),
        save_on = true,
        save_discretes = true,
        disco_dt_set = false,
        save_start = save_everystep || isempty(saveat) ||
            saveat isa Number || prob.tspan[1] in saveat,
        save_end = nothing,
        callback = nothing,
        dense = save_everystep && isempty(saveat) &&
            !default_linear_interpolation(prob, alg),
        calck = (callback !== nothing && callback !== CallbackSet()) ||
            (dense) || !isempty(saveat), # and no dense output
        dt = nothing,
        dtmin = eltype(prob.tspan)(0),
        dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])),
        force_dtmin = false,
        adaptive = anyadaptive(alg),
        abstol = nothing,
        reltol = nothing,
        gamma = nothing,
        qmin = nothing,
        qmax = nothing,
        qsteady_min = nothing,
        qsteady_max = nothing,
        beta1 = nothing,
        beta2 = nothing,
        qoldinit = nothing,
        fullnormalize = true,
        failfactor = 2,
        maxiters = anyadaptive(alg) ? 1000000 : typemax(Int),
        internalnorm = ODE_DEFAULT_NORM,
        internalopnorm = opnorm,
        isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
        unstable_check = ODE_DEFAULT_UNSTABLE_CHECK,
        verbose = Standard(),
        controller = nothing,
        timeseries_errors = true,
        dense_errors = false,
        advance_to_tstop = false,
        stop_at_next_tstop = false,
        initialize_save = true,
        progress = false,
        progress_steps = 1000,
        progress_name = "ODE",
        progress_message = ODE_DEFAULT_PROG_MESSAGE,
        progress_id = :OrdinaryDiffEq,
        userdata = nothing,
        allow_extrapolation = alg_extrapolates(alg),
        initialize_integrator = true,
        alias = ODEAliasSpecifier(),
        initializealg = DefaultInit(),
        rng = nothing,
        disco_probs = nothing,
        # SDE/RODE fields: accepted here so that SDE packages can delegate to
        # _ode_init and construct an ODEIntegrator with noise populated.
        save_noise = false,
        delta = nothing,
        W = nothing,
        P = nothing,
        sqdt = nothing,
        noise = nothing,
        c = nothing,
        rate_constants = nothing,
        # Pre-built cache for SDE delegation (skip alg_cache call)
        _cache = nothing,
        # Pre-built u/uprev for SDE delegation (cache holds references to these)
        _u = nothing,
        _uprev = nothing,
        seed = UInt64(0),
        kwargs...
    )
    # ODE/DAE-specific validation (skip for RODE/SDE problems)
    if !(prob isa SciMLBase.AbstractRODEProblem)
        if prob isa SciMLBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
            error("You cannot use an ODE Algorithm with a DAEProblem")
        end

        if prob isa SciMLBase.AbstractODEProblem && alg isa DAEAlgorithm
            error("You cannot use an DAE Algorithm with a ODEProblem")
        end

        if prob isa SciMLBase.ODEProblem
            if !(prob.f isa SciMLBase.DynamicalODEFunction) && alg isa PartitionedAlgorithm
                error("You can not use a solver designed for partitioned ODE with this problem. Please choose a solver suitable for your problem")
            end
        end

        if prob.f isa DynamicalODEFunction && prob.f.mass_matrix isa Tuple
            if any(mm != I for mm in prob.f.mass_matrix)
                error("This solver is not able to use mass matrices. For compatible solvers see https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/")
            end
        elseif !(prob isa SciMLBase.AbstractDiscreteProblem) &&
                !(prob isa SciMLBase.AbstractDAEProblem) &&
                !is_mass_matrix_alg(alg) &&
                prob.f.mass_matrix != I
            error("This solver is not able to use mass matrices. For compatible solvers see https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/")
        end
    end

    verbose_spec = _process_verbose_param(verbose)

    if alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm &&
            # https://github.com/SciML/OrdinaryDiffEq.jl/pull/2079 fixes this for Rosenbrock23 and 32
            !only_diagonal_mass_matrix(alg) &&
            prob.f.mass_matrix isa AbstractMatrix &&
            all(isequal(0), prob.f.mass_matrix)
        # technically this should also warn for zero operators but those are hard to check for
        if (dense || !isempty(saveat))
            @SciMLMessage(
                "Rosenbrock methods on equations without differential states do not bound the error on interpolations.",
                verbose_spec, :rosenbrock_no_differential_states
            )
        end
    end

    if only_diagonal_mass_matrix(alg) &&
            prob.f.mass_matrix isa AbstractMatrix &&
            !isdiag(prob.f.mass_matrix)
        throw(ArgumentError("$(typeof(alg).name.name) only works with diagonal mass matrices. Please choose a solver suitable for your problem (e.g. Rodas5P)"))
    end

    if !isempty(saveat) && dense
        @SciMLMessage(
            "Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.",
            verbose_spec, :dense_output_saveat
        )
    end

    progress && @logmsg(LogLevel(-1), progress_name, _id = progress_id, progress = 0)

    tType = eltype(prob.tspan)
    tspan = prob.tspan
    tdir = sign(tspan[end] - tspan[1])

    t = tspan[1]

    if (!isadaptive(alg) || !adaptive) &&
            isnothing(dt) && isempty(tstops) && dt_required(alg)
        throw(ArgumentError("Fixed timestep methods require a choice of dt or choosing the tstops"))
    end
    if !isadaptive(alg) && adaptive
        throw(ArgumentError("Fixed timestep methods can not be run with adaptive=true"))
    end

    isdae = if !isnothing(W)
        false # RODE/SDE — never DAE
    else
        alg isa DAEAlgorithm || (
            !(prob isa SciMLBase.AbstractDiscreteProblem) &&
                prob.f.mass_matrix != I &&
                !(prob.f.mass_matrix isa Tuple) &&
                ArrayInterface.issingular(prob.f.mass_matrix)
        )
    end
    if alg isa CompositeAlgorithm && alg.choice_function isa AutoSwitch
        auto = alg.choice_function
        _alg = CompositeAlgorithm(
            alg.algs,
            AutoSwitchCache(
                0, 0,
                auto.nonstiffalg,
                auto.stiffalg,
                auto.stiffalgfirst,
                auto.maxstiffstep,
                auto.maxnonstiffstep,
                auto.nonstifftol,
                auto.stifftol,
                auto.dtfac,
                auto.stiffalgfirst,
                auto.switch_max, 0
            )
        )
    else
        _alg = alg
    end

    use_old_kwargs = haskey(kwargs, :alias_u0) || haskey(kwargs, :alias_du0)

    aliases = nothing
    if use_old_kwargs
        aliases = ODEAliasSpecifier()
        if haskey(kwargs, :alias_u0)
            message = "`alias_u0` keyword argument is deprecated, to set `alias_u0`,
            please use an ODEAliasSpecifier, e.g. `solve(prob, alias = ODEAliasSpecifier(alias_u0 = true))"
            Base.depwarn(message, :init)
            Base.depwarn(message, :solve)
            aliases = ODEAliasSpecifier(alias_u0 = values(kwargs).alias_u0)
        else
            aliases = ODEAliasSpecifier(alias_u0 = nothing)
        end

        if haskey(kwargs, :alias_du0)
            message = "`alias_du0` keyword argument is deprecated, to set `alias_du0`,
            please use an ODEAliasSpecifier, e.g. `solve(prob, alias = ODEAliasSpecifier(alias_du0 = true))"
            Base.depwarn(message, :init)
            Base.depwarn(message, :solve)
            aliases = ODEAliasSpecifier(
                alias_u0 = aliases.alias_u0, alias_du0 = values(kwargs).alias_du0
            )
        else
            aliases = ODEAliasSpecifier(alias_u0 = aliases.alias_u0, alias_du0 = nothing)
        end

        aliases

    else
        # If alias isa Bool, all fields of ODEAliases set to alias
        if alias isa Bool
            aliases = ODEAliasSpecifier(alias = alias)
        elseif alias isa ODEAliasSpecifier
            aliases = alias
        end
    end

    if isnothing(aliases.alias_f) || aliases.alias_f
        f = prob.f
    else
        f = deepcopy(prob.f)
    end

    if isnothing(aliases.alias_p) || aliases.alias_p
        p = prob.p
    else
        p = recursivecopy(prob.p)
    end

    if !isnothing(aliases.alias_u0) && aliases.alias_u0
        u = prob.u0
    else
        u = recursivecopy(prob.u0)
    end

    # Handle null u0 (e.g., MTK systems with only callbacks and no state variables)
    # Convert to empty Float64 array to allow initialization to proceed
    if u === nothing
        u = Float64[]
    end

    # SDE delegation: use pre-built u if provided (cache holds references to it)
    if _u !== nothing
        u = _u
    end

    if _alg isa DAEAlgorithm
        if !isnothing(aliases.alias_du0) && aliases.alias_du0
            du = prob.du0
        else
            du = recursivecopy(prob.du0)
        end
        duprev = recursivecopy(du)
    else
        du = nothing
        duprev = nothing
    end

    uType = typeof(u)
    uBottomEltype = recursive_bottom_eltype(u)
    uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

    uEltypeNoUnits = recursive_unitless_eltype(u)
    tTypeNoUnits = typeof(one(tType))

    if prob isa SciMLBase.AbstractDiscreteProblem && abstol === nothing
        abstol_internal = false
    elseif abstol === nothing
        if uBottomEltypeNoUnits == uBottomEltype
            abstol_internal = unitfulvalue(
                real(
                    convert(
                        uBottomEltype,
                        oneunit(uBottomEltype) *
                            1 // 10^6
                    )
                )
            )
        else
            abstol_internal = unitfulvalue.(real.(oneunit.(u) .* 1 // 10^6))
        end
    else
        abstol_internal = real.(abstol)
    end

    if prob isa SciMLBase.AbstractDiscreteProblem && reltol === nothing
        reltol_internal = false
    elseif reltol === nothing
        if uBottomEltypeNoUnits == uBottomEltype
            reltol_internal = unitfulvalue(
                real(
                    convert(
                        uBottomEltype,
                        oneunit(uBottomEltype) * 1 // 10^3
                    )
                )
            )
        else
            reltol_internal = unitfulvalue.(real.(oneunit.(u) .* 1 // 10^3))
        end
    else
        reltol_internal = real.(reltol)
    end

    dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert
    # dtmin is all abs => does not care about sign already.

    if !isdae && isinplace(prob) && u isa AbstractArray && eltype(u) <: Number &&
            uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits # Could this be more efficient for other arrays?
        rate_prototype = recursivecopy(u)
    elseif prob isa DAEProblem
        rate_prototype = prob.du0
    else
        if (uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits) ||
                eltype(u) <: Enum
            rate_prototype = u
        else # has units!
            rate_prototype = u / oneunit(tType)
        end
    end
    rateType = typeof(rate_prototype) ## Can be different if united

    res_prototype = nothing
    if isdae
        if uBottomEltype == uBottomEltypeNoUnits
            res_prototype = u
        else
            res_prototype = one(u)
        end
        resType = typeof(res_prototype)
    end

    if isnothing(aliases.alias_tstops) || aliases.alias_tstops
        tstops = tstops
    else
        tstops = recursivecopy(tstops)
    end

    if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number
        _tstops_cache = tstops
    else
        _tstops_cache = tstops
        tstops = ()
    end
    tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan)
    saveat_internal = initialize_saveat(tType, saveat, tspan)
    d_discontinuities_internal = initialize_d_discontinuities(
        tType, d_discontinuities,
        tspan
    )

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length_int(callbacks_internal)
    if max_len_cb !== nothing
        uBottomEltypeReal = real(uBottomEltype)
        if isinplace(prob)
            callback_cache = DiffEqBase.CallbackCache(
                u, max_len_cb, uBottomEltypeReal,
                uBottomEltypeReal
            )
        else
            callback_cache = DiffEqBase.CallbackCache(
                max_len_cb, uBottomEltypeReal,
                uBottomEltypeReal
            )
        end
    else
        callback_cache = nothing
    end

    ### Algorithm-specific defaults ###
    save_idxs,
        saved_subsystem = SciMLBase.get_save_idxs_and_saved_subsystem(
        prob, save_idxs
    )

    ks_prototype = nothing
    if save_idxs === nothing
        ksEltype = Vector{rateType}
    else
        ks_prototype = rate_prototype[save_idxs]
        ksEltype = Vector{typeof(ks_prototype)}
    end

    # Have to convert in case passed in wrong.
    if save_idxs === nothing
        timeseries = timeseries_init === () ? uType[] :
            convert(Vector{uType}, timeseries_init)
    else
        u_initial = u[save_idxs]
        timeseries = timeseries_init === () ? typeof(u_initial)[] :
            convert(Vector{uType}, timeseries_init)
    end

    ts = ts_init === () ? tType[] : convert(Vector{tType}, ts_init)
    ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype}, ks_init)
    alg_choice = (
            _alg isa CompositeAlgorithm ||
            _alg isa StochasticDiffEqCompositeAlgorithm ||
            _alg isa StochasticDiffEqRODECompositeAlgorithm
        ) ? Int[] : nothing

    if (!adaptive || !isadaptive(_alg)) && save_everystep && tspan[2] - tspan[1] != Inf
        if isnothing(dt)
            steps = length(tstops)
        else
            # For fixed dt, the only time dtmin makes sense is if it's smaller than eps().
            # Therefore user specified dtmin doesn't matter, but we need to ensure dt>=eps()
            # to prevent infinite loops.
            abs(dt) < dtmin &&
                throw(ArgumentError("Supplied dt is smaller than dtmin"))
            steps = ceil(Int, internalnorm((tspan[2] - tspan[1]) / dt, tspan[1]))
        end
        sizehint!(timeseries, steps + 1)
        sizehint!(ts, steps + 1)
        sizehint!(ks, steps + 1)
    elseif save_everystep
        sizehint!(timeseries, 50)
        sizehint!(ts, 50)
        sizehint!(ks, 50)
    elseif !isempty(saveat_internal)
        savelength = length(saveat_internal) + 1
        if save_start == false
            savelength -= 1
        end
        if save_end == false && prob.tspan[2] in saveat_internal.valtree
            savelength -= 1
        end
        sizehint!(timeseries, savelength)
        sizehint!(ts, savelength)
        sizehint!(ks, savelength)
    else
        sizehint!(timeseries, 2)
        sizehint!(ts, 2)
        sizehint!(ks, 2)
    end

    k = rateType[]

    if uses_uprev(_alg, adaptive) || calck
        uprev = recursivecopy(u)
    else
        # Some algorithms do not use `uprev` explicitly. In that case, we can save
        # some memory by aliasing `uprev = u`, e.g. for "2N" low storage methods.
        uprev = u
    end
    # SDE delegation: use pre-built uprev if provided (cache holds references to it)
    if _uprev !== nothing
        uprev = _uprev
    end
    if allow_extrapolation
        uprev2 = recursivecopy(u)
    else
        uprev2 = uprev
    end

    _dt = if isnothing(dt)
        isdiscretealg(alg) && isempty(tstops) ?
            eltype(prob.tspan)(1) : eltype(prob.tspan)(0)
    else
        dt
    end
    if _cache !== nothing
        cache = _cache
    elseif prob isa DAEProblem
        cache = alg_cache(
            _alg, du, u, res_prototype, rate_prototype, uEltypeNoUnits,
            uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, _dt,
            reltol_internal, p, calck, Val(isinplace(prob)), verbose_spec
        )
    else
        cache = alg_cache(
            _alg, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits,
            tTypeNoUnits, uprev, uprev2, f, t, _dt, reltol_internal, p, calck,
            Val(isinplace(prob)), verbose_spec
        )
    end

    # Setting up the step size controller
    if (beta1 !== nothing || beta2 !== nothing) && controller !== nothing
        throw(ArgumentError("Setting both the legacy PID parameters `beta1, beta2 = $((beta1, beta2))` and the `controller = $controller` is not allowed."))
    end

    # Deprecation warnings for users to break down which parameters they accidentally set.
    if (beta1 !== nothing || beta2 !== nothing)
        message = "Providing the legacy PID parameters `beta1, beta2` is deprecated. Use the keyword argument `controller` instead."
        Base.depwarn(message, :init)
        Base.depwarn(message, :solve)
    end
    if (qmin !== nothing || qmax !== nothing)
        message = "Providing the legacy PID parameters `qmin, qmax` is deprecated. Use the keyword argument `controller` instead."
        Base.depwarn(message, :init)
        Base.depwarn(message, :solve)
    end
    if (qsteady_min !== nothing || qsteady_max !== nothing)
        message = "Providing the legacy PID parameters `qsteady_min, qsteady_max` is deprecated. Use the keyword argument `controller` instead."
        Base.depwarn(message, :init)
        Base.depwarn(message, :solve)
    end
    if (qoldinit !== nothing)
        message = "Providing the legacy PID parameters `qoldinit` is deprecated. Use the keyword argument `controller` instead."
        Base.depwarn(message, :init)
        Base.depwarn(message, :solve)
    end

    QT = determine_controller_datatype(u, internalnorm, tspan)

    # The following code provides an upgrade path for users by preserving the old behavior.
    legacy_controller_parameters = (gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit)
    if controller === nothing # We have to reconstruct the old controller before breaking release.
        if any(legacy_controller_parameters .== nothing)
            gamma = convert(QT, gamma === nothing ? gamma_default(alg) : gamma)
            qmin = convert(QT, qmin === nothing ? qmin_default(alg) : qmin)
            qmax = convert(QT, qmax === nothing ? qmax_default(alg) : qmax)
            qsteady_min = convert(QT, qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min)
            qsteady_max = convert(QT, qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max)
            qoldinit = convert(QT, qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit)
        end
        controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
    else # Controller has been passed
        gamma = hasfield(typeof(controller), :gamma) ? controller.gamma : gamma_default(alg)
        qmin = hasfield(typeof(controller), :qmin) ? controller.qmin : qmin_default(alg)
        qmax = hasfield(typeof(controller), :qmax) ? controller.qmax : qmax_default(alg)
        qsteady_min = hasfield(typeof(controller), :qsteady_min) ? controller.qsteady_min : qsteady_min_default(alg)
        qsteady_max = hasfield(typeof(controller), :qsteady_max) ? controller.qsteady_max : qsteady_max_default(alg)
        qoldinit = hasfield(typeof(controller), :qoldinit) ? controller.qoldinit : (anyadaptive(alg) ? 1 // 10^4 : 0)
    end

    EEstT = if tTypeNoUnits <: Integer
        promote_type(typeof(qmin), typeof(qmax))
    elseif prob isa SciMLBase.AbstractDiscreteProblem
        constvalue(tTypeNoUnits)
    else
        typeof(internalnorm(u, t))
    end

    atmp = if hasfield(typeof(cache), :atmp)
        cache.atmp
    else
        nothing
    end
    controller_cache = setup_controller_cache(_alg, atmp, controller)

    save_end_user = save_end
    save_end = save_end === nothing ?
        save_everystep || isempty(saveat) || saveat isa Number ||
        prob.tspan[2] in saveat : save_end

    opts = DEOptions{
        typeof(abstol_internal), typeof(reltol_internal),
        QT, tType,
        # typeof(controller),
        typeof(controller_cache),
        typeof(internalnorm), typeof(internalopnorm),
        typeof(save_end_user),
        typeof(callbacks_internal),
        typeof(isoutofdomain),
        typeof(progress_message), typeof(unstable_check),
        typeof(tstops_internal),
        typeof(d_discontinuities_internal), typeof(userdata),
        typeof(save_idxs),
        typeof(maxiters), typeof(_tstops_cache),
        typeof(saveat), typeof(d_discontinuities), typeof(verbose_spec),
        typeof(delta),
    }(
        maxiters, save_everystep,
        adaptive, abstol_internal,
        reltol_internal,
        # TODO vvv remove this block as these are controller and not integrator parameters vvv
        QT(gamma),
        QT(qmax),
        QT(qmin),
        QT(qsteady_max),
        QT(qsteady_min),
        QT(qoldinit),
        # TODO ^^^remove this block as these are controller and not integrator parameters ^^^
        QT(failfactor),
        tType(dtmax), tType(dtmin),
        # TODO vvv remove this vvv
        # controller,
        controller_cache,
        # TODO ^^^ remove this ^^^
        internalnorm,
        internalopnorm,
        save_idxs, tstops_internal,
        saveat_internal,
        d_discontinuities_internal,
        _tstops_cache, saveat,
        d_discontinuities,
        userdata, progress,
        progress_steps,
        progress_name,
        progress_message,
        progress_id,
        timeseries_errors,
        dense_errors, delta, dense,
        save_on, save_start,
        save_end, save_noise, save_discretes, save_end_user,
        callbacks_internal,
        isoutofdomain,
        unstable_check,
        verbose_spec, calck, force_dtmin,
        advance_to_tstop,
        stop_at_next_tstop
    )

    stats = SciMLBase.DEStats(0)
    differential_vars = if prob isa DAEProblem
        prob.differential_vars
    elseif !isnothing(W)
        nothing
    else
        get_differential_vars(f, u)
    end

    # SDE/RODE solvers never populate derivative stages (ks) — RODESolution has
    # no .k field, so _has_ks returns false and save_dense_at_t! is a no-op.
    # Force dense=false so ode_interpolation uses linear interpolation instead
    # of trying to access the empty ks vector.
    if !isnothing(W)
        dense = false
    end

    id = InterpolationData(
        f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars, false
    )

    _sol_kwargs = if !isnothing(W)
        (;
            W = W, seed = seed, interp = id, dense = dense, alg_choice = alg_choice,
            calculate_error = false, stats = stats, saved_subsystem = saved_subsystem,
        )
    else
        (;
            dense = dense, k = ks, interp = id, alg_choice = alg_choice,
            calculate_error = false, stats = stats, saved_subsystem = saved_subsystem,
        )
    end
    sol = SciMLBase.build_solution(prob, _alg, ts, timeseries; _sol_kwargs...)

    FType = typeof(f)
    SolType = typeof(sol)
    cacheType = typeof(cache)

    # rate/state = (state/time)/state = 1/t units, internalnorm drops units
    # we don't want to differentiate through eigenvalue estimation
    eigen_est = inv(one(tType))
    tprev = t
    dtcache = tType(_dt)
    dtpropose = tType(_dt)
    iter = 0
    kshortsize = 0
    reeval_fsal = false
    u_modified = false
    EEst = oneunit(EEstT) # https://github.com/JuliaPhysics/Measurements.jl/pull/135
    just_hit_tstop = false
    next_step_tstop = false
    tstop_target = zero(t)
    isout = false
    accept_step = false
    force_stepfail = false
    last_stepfail = false
    do_error_check = true
    event_last_time = 0
    vector_event_last_time = 1
    last_event_error = prob isa SciMLBase.AbstractDiscreteProblem ? false :
        (
            Base.isbitstype(uBottomEltypeNoUnits) ? zero(uBottomEltypeNoUnits) :
            0.0
        )
    dtchangeable = isdtchangeable(_alg)
    q11 = QT(1)
    success_iter = 0
    erracc = QT(1)
    dtacc = tType(1)
    reinitialize = true
    saveiter = 0 # Starts at 0 so first save is at 1
    saveiter_dense = 0
    fsalfirst, fsallast = _cache !== nothing ? (nothing, nothing) :
        get_fsalfirstlast(cache, rate_prototype)

    _rng = rng === nothing ? Random.default_rng() : rng

    num_probs = 0
    for i in callbacks_internal.continuous_callbacks
        if i.is_discontinuity
            num_probs += 1
        end
    end
    disco_probs = Vector{IntervalNonlinearProblem}(undef, num_probs)
    idx = 1
    for i in callbacks_internal.continuous_callbacks
        if i.is_discontinuity
            u₁ = similar(u)
            out = i isa VectorContinuousCallback ? similar(u) : nothing
            zero_func = zero_func_struct(u₁, i, _dt, uprev, u, k, cache, save_idxs, differential_vars, 1, out)
            disco_probs[idx] = IntervalNonlinearProblem(zero_func, [zero(tType), one(tType)], p)
            idx += 1
        end
    end 

    integrator = ODEIntegrator{
        typeof(_alg), isinplace(prob), uType, typeof(du),
        tType, typeof(p),
        typeof(eigen_est), EEstT,
        QT, typeof(tdir), typeof(k), SolType,
        FType, cacheType,
        typeof(opts), typeof(fsalfirst),
        typeof(last_event_error), typeof(callback_cache),
        typeof(initializealg), typeof(differential_vars),
        typeof(controller_cache), typeof(_rng),
        typeof(W), typeof(P), typeof(sqdt),
        typeof(noise), typeof(c), typeof(rate_constants)
    }(
        sol, u, du, k, t, tType(_dt), f, p,
        uprev, uprev2, duprev, tprev,
        _alg, dtcache, dtchangeable,
        dtpropose, disco_dt_set, tdir, eigen_est, EEst,
        # TODO vvv remove these
        QT(qoldinit), q11,
        erracc, dtacc,
        # TODO ^^^ remove these
        controller_cache,
        success_iter,
        iter, saveiter, saveiter_dense, cache,
        callback_cache,
        kshortsize, force_stepfail,
        last_stepfail,
        just_hit_tstop, next_step_tstop, tstop_target, do_error_check,
        event_last_time,
        vector_event_last_time,
        last_event_error, accept_step,
        isout, reeval_fsal,
        u_modified, reinitialize, isdae,
        opts, stats, initializealg, differential_vars,
        fsalfirst, fsallast, _rng, disco_probs,
        W, P, sqdt,
        noise, c, rate_constants, QT(1)
    )

    if initialize_integrator
        if isdae || SciMLBase.has_initializeprob(prob.f) ||
                prob isa SciMLBase.ImplicitDiscreteProblem
            DiffEqBase.initialize_dae!(integrator)
            !isnothing(integrator.u) && update_uprev!(integrator)
        end

        if save_start
            integrator.saveiter += 1 # Starts at 1 so first save is at 2
            integrator.saveiter_dense += 1
            copyat_or_push!(ts, 1, t)
            # N.B.: integrator.u can be modified by initialized_dae!
            if save_idxs === nothing
                copyat_or_push!(timeseries, 1, integrator.u)
                if isnothing(W)
                    copyat_or_push!(ks, 1, [rate_prototype])
                end
            else
                copyat_or_push!(timeseries, 1, integrator.u[save_idxs], Val{false})
                if isnothing(W)
                    copyat_or_push!(ks, 1, [ks_prototype])
                end
            end
        else
            integrator.saveiter = 0 # Starts at 0 so first save is at 1
            integrator.saveiter_dense = 0
        end

        initialize_callbacks!(integrator, initialize_save)
        initialize!(integrator, integrator.cache)

        if _alg isa OrdinaryDiffEqCompositeAlgorithm
            # in case user mixes adaptive and non-adaptive algorithms
            ensure_behaving_adaptivity!(integrator, integrator.cache)
        end
        if alg_choice !== nothing && save_start
            # Loop to get all of the extra possible saves in callback initialization
            for i in 1:(integrator.saveiter)
                copyat_or_push!(alg_choice, i, integrator.cache.current)
            end
        end
    end

    if !(_tstops_cache isa AbstractArray || _tstops_cache isa Tuple || _tstops_cache isa Number)
        tstops = _tstops_cache(parameter_values(integrator), prob.tspan)
        for tstop in tstops
            add_tstop!(integrator, tstop)
        end
    end

    handle_dt!(integrator, dt)

    # Noise process initialization for SDE/RODE integrators
    if !isnothing(integrator.W)
        modify_dt_for_tstops!(integrator)
        integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))
        integrator.W.dt = integrator.dt
        !isnothing(integrator.P) && (integrator.P.dt = integrator.dt)
    end

    return integrator
end

function SciMLBase.solve!(integrator::ODEIntegrator)
    @inbounds while !isempty(integrator.opts.tstops)
        first_tstop = first(integrator.opts.tstops)
        while integrator.tdir * integrator.t < first_tstop
            loopheader!(integrator)
            if integrator.do_error_check && check_error!(integrator) != ReturnCode.Success
                return integrator.sol
            end

            # Use special tstop handling if flag is set, otherwise normal stepping
            if integrator.next_step_tstop
                handle_tstop_step!(integrator)
            else
                perform_step!(integrator, integrator.cache)
            end

            should_exit = integrator.next_step_tstop

            loopfooter!(integrator)
            if isempty(integrator.opts.tstops) || should_exit
                break
            end
        end
        handle_tstop!(integrator)
    end
    postamble!(integrator)

    f = integrator.sol.prob.f
    # SDE split problems may store f as a Tuple; unwrap to check for analytic solution
    f = f isa Tuple ? f[1] : f

    if SciMLBase.has_analytic(f)
        SciMLBase.calculate_solution_errors!(
            integrator.sol;
            timeseries_errors = integrator.opts.timeseries_errors,
            dense_errors = integrator.opts.dense_errors
        )
    end
    if integrator.sol.retcode != ReturnCode.Default
        return integrator.sol
    end
    return integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.Success)
end

# Helpers

function handle_dt!(integrator)
    return if iszero(integrator.dt) && integrator.opts.adaptive
        auto_dt_reset!(integrator)
        if sign(integrator.dt) != integrator.tdir && !iszero(integrator.dt) &&
                !isnan(integrator.dt)
            error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
        end
        if isnan(integrator.dt)
            @SciMLMessage(
                "Automatic dt set the starting dt as NaN, causing instability. Exiting.",
                integrator.opts.verbose, :dt_NaN
            )
        end
    elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) &&
            integrator.tdir < 0
        integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
    end
end
function handle_dt!(integrator, dt)
    return if isnothing(dt) && iszero(integrator.dt) && integrator.opts.adaptive
        auto_dt_reset!(integrator)
        if sign(integrator.dt) != integrator.tdir && !iszero(integrator.dt) &&
                !isnan(integrator.dt)
            error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
        end
        if isnan(integrator.dt)
            @SciMLMessage(
                "Automatic dt set the starting dt as NaN, causing instability. Exiting.",
                integrator.opts.verbose, :dt_NaN
            )
        end
    elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) &&
            integrator.tdir < 0
        integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
    end
end

# time stops
@inline function initialize_tstops(::Type{T}, tstops, d_discontinuities, tspan) where {T}
    tstops_internal = BinaryHeap{T}(DataStructures.FasterForward())

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    for t in tstops
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    for t in d_discontinuities
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    push!(tstops_internal, tdir_tf)

    return tstops_internal
end

function reinit_tstops!(
        ::Type{T}, tstops_internal, tstops, d_discontinuities, tspan; p = nothing
    ) where {T}
    empty!(tstops_internal)

    # Evaluate callable tstops
    _tstops = if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number
        tstops
    elseif p !== nothing
        tstops(p, tspan)
    else
        ()
    end

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    for t in _tstops
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    for t in d_discontinuities
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    return push!(tstops_internal, tdir_tf)
end

# saving time points
function initialize_saveat(::Type{T}, saveat, tspan) where {T}
    saveat_internal = BinaryHeap{T}(DataStructures.FasterForward())

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    if saveat isa Number
        directional_saveat = tdir * abs(saveat)
        for t in (t0 + directional_saveat):directional_saveat:tf
            push!(saveat_internal, tdir * t)
        end
    elseif !isempty(saveat)
        for t in saveat
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(saveat_internal, tdir_t)
        end
    end

    return saveat_internal
end

function reinit_saveat!(::Type{T}, saveat_internal, saveat, tspan) where {T}
    empty!(saveat_internal)

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    return if saveat isa Number
        directional_saveat = tdir * abs(saveat)
        for t in (t0 + directional_saveat):directional_saveat:tf
            push!(saveat_internal, tdir * t)
        end
    elseif !isempty(saveat)
        for t in saveat
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(saveat_internal, tdir_t)
        end
    end
end

# discontinuities
function initialize_d_discontinuities(::Type{T}, d_discontinuities, tspan) where {T}
    d_discontinuities_internal = BinaryHeap{T}(DataStructures.FasterForward())
    sizehint!(d_discontinuities_internal, length(d_discontinuities))

    t0, tf = tspan
    tdir = sign(tf - t0)

    for t in d_discontinuities
        push!(d_discontinuities_internal, tdir * t)
    end

    return d_discontinuities_internal
end

function reinit_d_discontinuities!(::Type{T}, d_discontinuities_internal, d_discontinuities, tspan) where {T}
    empty!(d_discontinuities_internal)

    t0, tf = tspan
    tdir = sign(tf - t0)

    for t in d_discontinuities
        push!(d_discontinuities_internal, tdir * t)
    end
    return
end

function initialize_callbacks!(integrator, initialize_save = true)
    t = integrator.t
    u = integrator.u
    callbacks = integrator.opts.callback
    integrator.u_modified = true

    u_modified = initialize!(callbacks, u, t, integrator)

    # if the user modifies u, we need to fix previous values before initializing
    # FSAL in order for the starting derivatives to be correct
    if u_modified
        if isinplace(integrator.sol.prob)
            recursivecopy!(integrator.uprev, integrator.u)
        else
            integrator.uprev = integrator.u
        end

        if alg_extrapolates(integrator.alg)
            if isinplace(integrator.sol.prob)
                recursivecopy!(integrator.uprev2, integrator.uprev)
            else
                integrator.uprev2 = integrator.uprev
            end
        end

        if initialize_save &&
                (
                any((c) -> c.save_positions[2], callbacks.discrete_callbacks) ||
                    any((c) -> c.save_positions[2], callbacks.continuous_callbacks)
            )
            savevalues!(integrator, true)
        end
    end

    # reset this as it is now handled so the integrators should proceed as normal
    integrator.u_modified = false

    return if initialize_save
        SciMLBase.save_discretes_if_enabled!(integrator, integrator.opts.callback; skip_duplicates = true)
    end
end