function abs2_and_sum(x, y)
    return reduce(+, x, init = zero(real(value(eltype(x))))) +
        reduce(+, y, init = zero(real(value(eltype(y)))))
end
UNITLESS_ABS2(x::Number) = abs2(x)

function UNITLESS_ABS2(x::AbstractArray{<:Number})
    return mapreduce(UNITLESS_ABS2, +, x, init = zero(real(value(eltype(x)))))
end

function UNITLESS_ABS2(x::AbstractArray)
    return mapreduce(UNITLESS_ABS2, abs2_and_sum, x, init = zero(real(value(eltype(x)))))
end
function UNITLESS_ABS2(x::RecursiveArrayTools.AbstractVectorOfArray)
    return mapreduce(UNITLESS_ABS2, abs2_and_sum, x.u, init = zero(real(value(eltype(x)))))
end
function UNITLESS_ABS2(x::RecursiveArrayTools.ArrayPartition)
    return mapreduce(UNITLESS_ABS2, abs2_and_sum, x.x, init = zero(real(value(eltype(x)))))
end

function UNITLESS_ABS2(x::RecursiveArrayTools.AbstractRaggedVectorOfArray)
    return mapreduce(UNITLESS_ABS2, +, x.u; init = zero(real(eltype(x))))
end

UNITLESS_ABS2(f::F, x::Number) where {F} = abs2(f(x))
function UNITLESS_ABS2(f::F, x::AbstractArray) where {F}
    return mapreduce(
        UNITLESS_ABS2 ∘ f, abs2_and_sum, x;
        init = zero(real(value(eltype(x))))
    )
end
function UNITLESS_ABS2(f::F, x::RecursiveArrayTools.ArrayPartition) where {F}
    return mapreduce(
        UNITLESS_ABS2 ∘ f, abs2_and_sum, x.x;
        init = zero(real(value(eltype(x))))
    )
end

recursive_length(u::AbstractArray{<:Number}) = length(u)
recursive_length(u::Number) = length(u)
recursive_length(u::AbstractArray{<:AbstractArray}) = sum(recursive_length, u)
recursive_length(u::RecursiveArrayTools.ArrayPartition) = sum(recursive_length, u.x)
recursive_length(u::RecursiveArrayTools.VectorOfArray) = sum(recursive_length, u.u)
recursive_length(u::RecursiveArrayTools.AbstractRaggedVectorOfArray) = sum(recursive_length, u.u; init = 0)
function recursive_length(
        u::AbstractArray{
            <:StaticArraysCore.StaticArray{S, <:Number},
        }
    ) where {S}
    return prod(Size(eltype(u))) * length(u)
end

"""
    ODE_DEFAULT_NORM(u, t)
    ODE_DEFAULT_NORM(f, u, t)

The default internal norm used by the integrators for error estimation and step-size
control. It is the (optionally `f`-weighted) RMS norm: roughly `sqrt(sum(abs2, u) / length(u))`,
which scales with the magnitude of the state but not its dimensionality, with specialized
methods for scalars, `Array`s, static arrays, and nested array types. Pass a custom callable
via the `internalnorm` solver keyword to override it.
"""
ODE_DEFAULT_NORM(u::Union{AbstractFloat, Complex}, t) = @fastmath abs(u)

function ODE_DEFAULT_NORM(f::F, u::Union{AbstractFloat, Complex}, t) where {F}
    return @fastmath abs(f(u))
end

function ODE_DEFAULT_NORM(u::Array{T}, t) where {T <: Union{AbstractFloat, Complex}}
    x = zero(T)
    @inbounds @fastmath for ui in u
        x += abs2(ui)
    end
    return Base.FastMath.sqrt_fast(real(x) / max(length(u), 1))
end

function ODE_DEFAULT_NORM(
        f::F,
        u::Array{T},
        t
    ) where {F, T <: Union{AbstractFloat, Complex}}
    x = zero(T)
    @inbounds @fastmath for ui in u
        x += abs2(f(ui))
    end
    return Base.FastMath.sqrt_fast(real(x) / max(length(u), 1))
end

function ODE_DEFAULT_NORM(
        f::F,
        u::Iterators.Zip{Z},
        t
    ) where {F, Z <: Tuple{Vararg{Array{<:Union{AbstractFloat, Complex}}}}}
    T = eltype(first(u.is))
    x = zero(T)
    @inbounds @fastmath for ui in u
        x += abs2(f(ui))
    end
    return Base.FastMath.sqrt_fast(real(x) / max(length(u), 1))
end

function ODE_DEFAULT_NORM(
        u::StaticArraysCore.StaticArray{<:Tuple, T},
        t
    ) where {T <: Union{AbstractFloat, Complex}}
    return Base.FastMath.sqrt_fast(real(sum(abs2, u)) / max(length(u), 1))
end

function ODE_DEFAULT_NORM(
        f::F, u::StaticArraysCore.StaticArray{<:Tuple, T},
        t
    ) where {F, T <: Union{AbstractFloat, Complex}}
    return Base.FastMath.sqrt_fast(real(sum(abs2 ∘ f, u)) / max(length(u), 1))
end

function ODE_DEFAULT_NORM(
        u::Union{
            AbstractArray,
            RecursiveArrayTools.AbstractVectorOfArray,
            RecursiveArrayTools.AbstractRaggedVectorOfArray,
        },
        t
    )
    return Base.FastMath.sqrt_fast(UNITLESS_ABS2(u) / max(recursive_length(u), 1))
end

function ODE_DEFAULT_NORM(f::F, u::AbstractArray, t) where {F}
    return Base.FastMath.sqrt_fast(UNITLESS_ABS2(f, u) / max(recursive_length(u), 1))
end

ODE_DEFAULT_NORM(u, t) = norm(u)
ODE_DEFAULT_NORM(f::F, u, t) where {F} = norm(f.(u))

"""
    ODE_DEFAULT_ISOUTOFDOMAIN(u, p, t)

The default `isoutofdomain` predicate used by the integrators. It always returns `false`,
i.e. no state is considered out of the problem's domain. Pass a custom predicate via the
`isoutofdomain` solver keyword to reject steps whose proposed state leaves a valid domain.
"""
ODE_DEFAULT_ISOUTOFDOMAIN(u, p, t) = false

"""
    ODE_DEFAULT_PROG_MESSAGE(dt, u, p, t)

The default progress-bar message builder used by the integrators when `progress = true`.
It returns a short multi-line string reporting the current `dt`, `t`, and the largest-magnitude
component of the state `u`. Pass a custom callable via the `progress_message` solver keyword
to override it.
"""
function ODE_DEFAULT_PROG_MESSAGE(dt, u::Array, p, t)
    tmp = u[1]
    for i in eachindex(u)
        tmp = ifelse(abs(u[i]) > abs(tmp), u[i], tmp)
    end
    return "dt=" * string(dt) * "\nt=" * string(t) * "\nmax u=" * string(tmp)
end
function ODE_DEFAULT_PROG_MESSAGE(dt, u, p, t)
    return "dt=" * string(dt) * "\nt=" * string(t) * "\nmax u=" * string(maximum(abs.(u)))
end

"""
    NAN_CHECK(x)

Recursively test whether `x` holds a `NaN`. The integrators use this to detect a step
that produced `NaN` and reject it.

Methods are provided for numbers, `AbstractArray`s, `RecursiveArrayTools.AbstractVectorOfArray`s
and `RecursiveArrayTools.ArrayPartition`s; nested containers are descended into, so a state
made of arrays of arrays reports `true` if any leaf is `NaN`. `Enum` values always report
`false`, which keeps discrete components of a mixed state vector from being sent through
`isnan`. Add a method to make the check reach the elements of a custom state type:

```julia
DiffEqBase.NAN_CHECK(x::MyStateType) = any(DiffEqBase.NAN_CHECK, x.parts)
```

`NaN` is only one of the ways a step can go bad; `Inf` and overflow are handled separately
by [`ODE_DEFAULT_UNSTABLE_CHECK`](@ref).
"""
NAN_CHECK(x::Number) = isnan(x)
NAN_CHECK(x::Enum) = false
function NAN_CHECK(x::Union{AbstractArray, RecursiveArrayTools.AbstractVectorOfArray})
    return any(
        NAN_CHECK, x
    )
end
NAN_CHECK(x::RecursiveArrayTools.ArrayPartition) = any(NAN_CHECK, x.x)

INFINITE_OR_GIANT(x::Number) = !isfinite(x)
function INFINITE_OR_GIANT(
        x::Union{
            AbstractArray, RecursiveArrayTools.AbstractVectorOfArray,
        }
    )
    return any(
        INFINITE_OR_GIANT, x
    )
end
INFINITE_OR_GIANT(x::RecursiveArrayTools.ArrayPartition) = any(INFINITE_OR_GIANT, x.x)

"""
    ODE_DEFAULT_UNSTABLE_CHECK(dt, u, p, t) -> Bool

Return whether the default ODE instability check considers the current state
unstable.

The generic fallback returns `false`. Numeric scalars, arrays, and
`ArrayPartition`s return `true` when any state entry is infinite or non-finite.
This is the default used by OrdinaryDiffEq solvers when no `unstable_check`
callback is supplied.

# Arguments
- `dt`: Current step size.
- `u`: Current state.
- `p`: Problem parameters.
- `t`: Current time.

# Returns
- `Bool`: `true` when the state should be treated as unstable.
"""
ODE_DEFAULT_UNSTABLE_CHECK(dt, u, p, t) = false
function ODE_DEFAULT_UNSTABLE_CHECK(dt, u::Union{Number, AbstractArray{<:Number}}, p, t)
    return INFINITE_OR_GIANT(u)
end
