diff --git a/Project.toml b/Project.toml
index e84111a7..90b86f9b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DataInterpolations"
 uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
-version = "7.2.0"
+version = "8.0.0"
 
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
diff --git a/docs/Project.toml b/docs/Project.toml
index 9d20b018..56a44069 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
@@ -11,7 +10,6 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-DataInterpolations = "6, 7.0"
 Documenter = "1"
 ModelingToolkit = "9"
 ModelingToolkitStandardLibrary = "2"
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 3818a9e3..f84cc623 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -19,4 +19,5 @@ QuinticHermiteSpline
 
 ```@docs
 DataInterpolations.looks_linear
+DataInterpolations.output_dim
 ```
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index bf9efcb3..749e4bd5 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -2,7 +2,7 @@ module DataInterpolations
 
 ### Interface Functionality
 
-abstract type AbstractInterpolation{T, N} end
+abstract type AbstractInterpolation{T} end
 
 using LinearAlgebra, RecipesBase
 using PrettyTables
@@ -98,16 +98,28 @@ function Base.showerror(io::IO, ::ExtrapolationNotImplementedError)
     print(io, EXTRAPOLATION_NOT_IMPLEMENTED_ERROR)
 end
 
+"""
+    output_dim(x::AbstractInterpolation)
+
+Return the number of dimension `ndims(x(t))` of interpolation `x` evaluated at a single value `t`
+if `x(t) isa AbstractArray`, or 0 otherwise.
+"""
+output_dim(x::AbstractInterpolation) = _output_dim(x.u)
+_output_dim(::AbstractVector) = 0 # each value is a scalar
+_output_dim(::AbstractVector{<:AbstractArray{<:Any, N}}) where {N} = N # each value is an array but values are not stacked
+_output_dim(::AbstractArray{<:Any, N}) where {N} = N - 1 # each value is an array but multiple values are stacked
+
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
        BSplineInterpolation, BSplineApprox, CubicHermiteSpline, PCHIPInterpolation,
        QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv,
        ExtrapolationType
+export output_dim
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation
-struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation{T, N}} <:
-       AbstractInterpolation{T, N}
+struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}} <:
+       AbstractInterpolation{T}
     u::uType
     û::uType
     t::tType
@@ -132,8 +144,7 @@ struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation
             Aitp,
             extrapolation_left,
             extrapolation_right)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), eltype(u), typeof(λ), N, typeof(Aitp)}(
+        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(
             u,
             û,
             t,
@@ -161,9 +172,8 @@ struct CurvefitCache{
     lbType,
     algType,
     pminType,
-    T,
-    N
-} <: AbstractInterpolation{T, N}
+    T
+} <: AbstractInterpolation{T}
     u::uType
     t::tType
     m::mType        # model type
@@ -174,10 +184,9 @@ struct CurvefitCache{
     pmin::pminType  # optimized params
     extrapolate::Bool
     function CurvefitCache(u, t, m, p0, ub, lb, alg, pmin, extrapolate)
-        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(m),
             typeof(p0), typeof(ub), typeof(lb),
-            typeof(alg), typeof(pmin), eltype(u), N}(u,
+            typeof(alg), typeof(pmin), eltype(u)}(u,
             t,
             m,
             p0,
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 271a8419..2ec303df 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -209,7 +209,7 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
 end
 
 function _derivative(
-        A::BSplineInterpolation{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N}
+        A::BSplineInterpolation{<:AbstractArray{<:Number}}, t::Number, iguess)
     # change t into param [0 1]
     ax_u = axes(A.u)[1:(end - 1)]
     t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
@@ -254,7 +254,7 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
 end
 
 function _derivative(
-        A::BSplineApprox{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N}
+        A::BSplineApprox{<:AbstractArray{<:Number}}, t::Number, iguess)
     # change t into param [0 1]
     ax_u = axes(A.u)[1:(end - 1)]
     t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
index 577b98e3..a11d278d 100644
--- a/src/integral_inverses.jl
+++ b/src/integral_inverses.jl
@@ -1,4 +1,4 @@
-abstract type AbstractIntegralInverseInterpolation{T, N} <: AbstractInterpolation{T, N} end
+abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end
 
 """
     invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation
@@ -33,8 +33,8 @@ Can be easily constructed with `invert_integral(A::LinearInterpolation{<:Abstrac
   - `t` : Given by `A.I` (the cumulative integral of `A`)
   - `A` : The `LinearInterpolation` object
 """
-struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <:
-       AbstractIntegralInverseInterpolation{T, N}
+struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
     u::uType
     t::tType
     extrapolation_left::ExtrapolationType.T
@@ -42,8 +42,7 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <:
     iguesser::Guesser{tType}
     itp::itpType
     function LinearInterpolationIntInv(u, t, A)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
             u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A)
     end
 end
@@ -85,8 +84,8 @@ Can be easily constructed with `invert_integral(A::ConstantInterpolation{<:Abstr
   - `t` : Given by `A.I` (the cumulative integral of `A`)
   - `A` : The `ConstantInterpolation` object
 """
-struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <:
-       AbstractIntegralInverseInterpolation{T, N}
+struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
     u::uType
     t::tType
     extrapolation_left::ExtrapolationType.T
@@ -94,8 +93,7 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <:
     iguesser::Guesser{tType}
     itp::itpType
     function ConstantInterpolationIntInv(u, t, A)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
             u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A
         )
     end
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index aa51f262..aab592d1 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -27,7 +27,7 @@ Extrapolation extends the last linear polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct LinearInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpolation{T, N}
+struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -40,8 +40,7 @@ struct LinearInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpol
     function LinearInterpolation(u, t, I, p, extrapolation_left, extrapolation_right,
             cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
             u, t, I, p, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
     end
@@ -93,8 +92,8 @@ Extrapolation extends the last quadratic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuadraticInterpolation{uType, tType, IType, pType, T, N} <:
-       AbstractInterpolation{T, N}
+struct QuadraticInterpolation{uType, tType, IType, pType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -111,8 +110,7 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T, N} <:
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(I), typeof(p.α), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(I), typeof(p.α), eltype(u)}(
             u, t, I, p, mode, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
     end
@@ -161,8 +159,8 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa
   - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for
     the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`.
 """
-struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <:
-       AbstractInterpolation{T, N}
+struct LagrangeInterpolation{uType, tType, T, bcacheType} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     n::Int
@@ -175,8 +173,7 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <:
         bcache = zeros(eltype(u[1]), n + 1)
         idxs = zeros(Int, n + 1)
         fill!(bcache, NaN)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), eltype(u), typeof(bcache), N}(u,
+        new{typeof(u), typeof(t), eltype(u), typeof(bcache)}(u,
             t,
             n,
             bcache,
@@ -229,8 +226,8 @@ Extrapolation extends the last cubic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <:
-       AbstractInterpolation{T, N}
+struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -246,9 +243,8 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <:
             u, t, I, b, c, d, extrapolation_left,
             extrapolation_right, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
-            typeof(d), eltype(u), N}(u,
+            typeof(d), eltype(u)}(u,
             t,
             I,
             b,
@@ -328,7 +324,7 @@ Extrapolation extends the last constant polynomial at the end points on each sid
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct ConstantInterpolation{uType, tType, IType, T, N} <: AbstractInterpolation{T, N}
+struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -343,8 +339,7 @@ struct ConstantInterpolation{uType, tType, IType, T, N} <: AbstractInterpolation
             u, t, I, dir, extrapolation_left, extrapolation_right,
             cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(I), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(I), eltype(u)}(
             u, t, I, nothing, dir, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
     end
@@ -393,8 +388,8 @@ Extrapolation extends the last quadratic polynomial on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} <:
-       AbstractInterpolation{T, N}
+struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -411,9 +406,8 @@ struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} <
             u, t, I, p, k, c, sc, extrapolation_left,
             extrapolation_right, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k),
-            typeof(c), typeof(sc), eltype(u), N}(u,
+            typeof(c), typeof(sc), eltype(u)}(u,
             t,
             I,
             p,
@@ -516,8 +510,8 @@ Second derivative on both ends are zero, which are also called "natural" boundar
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <:
-       AbstractInterpolation{T, N}
+struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     I::IType
@@ -532,9 +526,8 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <:
     function CubicSpline(u, t, I, p, h, z, extrapolation_left,
             extrapolation_right, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁),
-            typeof(h), typeof(z), eltype(u), N}(
+            typeof(h), typeof(z), eltype(u)}(
             u,
             t,
             I,
@@ -550,12 +543,11 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <:
     end
 end
 
-function CubicSpline(u::uType,
+function CubicSpline(u::AbstractVector{<:Number},
         t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false,
-        assume_linear_t = 1e-2) where {uType <:
-                                       AbstractVector{<:Number}}
+        assume_linear_t = 1e-2)
     extrapolation_left, extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
@@ -585,12 +577,11 @@ function CubicSpline(u::uType,
         extrapolation_right, cache_parameters, linear_lookup)
 end
 
-function CubicSpline(u::uType,
+function CubicSpline(u::AbstractArray{T, N},
         t;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false,
-        assume_linear_t = 1e-2) where {uType <:
-                                       AbstractArray{T, N}} where {T, N}
+        assume_linear_t = 1e-2) where {T, N}
     extrapolation_left, extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
@@ -625,11 +616,10 @@ function CubicSpline(u::uType,
 end
 
 function CubicSpline(
-        u::uType, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        u::AbstractVector, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false,
-        assume_linear_t = 1e-2) where {uType <:
-                                       AbstractVector}
+        assume_linear_t = 1e-2)
     extrapolation_left, extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
@@ -685,8 +675,8 @@ Extrapolation is a constant polynomial of the end points on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <:
-       AbstractInterpolation{T, N}
+struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     d::Int    # degree
@@ -713,8 +703,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <:
             extrapolation_right,
             assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u)}(
             u,
             t,
             d,
@@ -808,11 +797,11 @@ function BSplineInterpolation(
 end
 
 function BSplineInterpolation(
-        u::AbstractArray{T, N}, t, d, pVecType, knotVecType;
+        u::AbstractArray, t, d, pVecType, knotVecType;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
-        assume_linear_t = 1e-2) where {T, N}
+        assume_linear_t = 1e-2)
     extrapolation_left, extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
@@ -917,8 +906,8 @@ Extrapolation is a constant polynomial of the end points on each side.
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <:
-       AbstractInterpolation{T, N}
+struct BSplineApprox{uType, tType, pType, kType, cType, scType, T} <:
+       AbstractInterpolation{T}
     u::uType
     t::tType
     d::Int    # degree
@@ -948,8 +937,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <:
             assume_linear_t
     )
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u)}(
             u,
             t,
             d,
@@ -1146,7 +1134,7 @@ function BSplineApprox(
         q[ax_u..., k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] -
                         sc[k, h] * u[ax_u..., end]
     end
-    Q = Array{eltype(u), N}(undef, size(u)[1:(end - 1)]..., h - 2)
+    Q = Array{T, N}(undef, size(u)[1:(end - 1)]..., h - 2)
     for i in 2:(h - 1)
         s = zeros(eltype(sc), size(u)[1:(end - 1)]...)
         for k in 2:(n - 1)
@@ -1192,8 +1180,8 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct CubicHermiteSpline{uType, tType, IType, duType, pType, T, N} <:
-       AbstractInterpolation{T, N}
+struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <:
+       AbstractInterpolation{T}
     du::duType
     u::uType
     t::tType
@@ -1208,8 +1196,7 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T, N} <:
             du, u, t, I, p, extrapolation_left, extrapolation_right,
             cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
-        new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
             du, u, t, I, p, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
     end
@@ -1295,8 +1282,8 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
     for a test based on the normalized standard deviation of the difference with respect
     to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T, N} <:
-       AbstractInterpolation{T, N}
+struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
+       AbstractInterpolation{T}
     ddu::dduType
     du::duType
     u::uType
@@ -1312,9 +1299,8 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T, N} <
             ddu, du, u, t, I, p, extrapolation_left,
             extrapolation_right, cache_parameters, assume_linear_t)
         linear_lookup = seems_linear(assume_linear_t, t)
-        N = get_output_dim(u)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
-            typeof(ddu), typeof(p.c₁), eltype(u), N}(
+            typeof(ddu), typeof(p.c₁), eltype(u)}(
             ddu, du, u, t, I, p, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
     end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index a853f7be..f2dd80ad 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -230,7 +230,7 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
     I + C + D
 end
 
-function _interpolate(A::CubicSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N}
+function _interpolate(A::CubicSpline{<:AbstractArray}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     Δt₁ = t - A.t[idx]
     Δt₂ = A.t[idx + 1] - t
@@ -261,9 +261,9 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     ucum
 end
 
-function _interpolate(A::BSplineInterpolation{<:AbstractArray{T, N}},
+function _interpolate(A::BSplineInterpolation{<:AbstractArray{<:Number}},
         t::Number,
-        iguess) where {T <: Number, N}
+        iguess)
     ax_u = axes(A.u)[1:(end - 1)]
     t < A.t[1] && return A.u[ax_u..., 1]
     t > A.t[end] && return A.u[ax_u..., end]
@@ -297,7 +297,7 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
 end
 
 function _interpolate(
-        A::BSplineApprox{<:AbstractArray{T, N}}, t::Number, iguess) where {T <: Number, N}
+        A::BSplineApprox{<:AbstractArray{<:Number}}, t::Number, iguess)
     ax_u = axes(A.u)[1:(end - 1)]
     t < A.t[1] && return A.u[ax_u..., 1]
     t > A.t[end] && return A.u[ax_u..., end]
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index cb172b52..85edac07 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -59,19 +59,6 @@ function spline_coefficients!(N, d, k, u::AbstractVector)
     return nothing
 end
 
-# Get Output Dimension for parameterizing AbstractInterpolations
-function get_output_dim(u::AbstractVector{<:Number})
-    return (1,)
-end
-
-function get_output_dim(u::AbstractVector)
-    return (length(first(u)),)
-end
-
-function get_output_dim(u::AbstractArray)
-    return size(u)[1:(end - 1)]
-end
-
 function quadratic_spline_params(t::AbstractVector, sc::AbstractVector)
 
     # Create knot vector
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index aecde638..16e33562 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -32,7 +32,8 @@ end
     for t in (1.0:10.0, 1.0collect(1:10))
         u = 2.0collect(1:10)
         #t = 1.0collect(1:10)
-        A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(LinearInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension))
 
         for (_t, _u) in zip(t, u)
             @test A(_t) == _u
@@ -40,9 +41,16 @@ end
         @test A(0) == 0.0
         @test A(5.5) == 11.0
         @test A(11) == 22
+        @test @inferred(output_dim(A)) == 0
 
         u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
-        A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+        @test @inferred(LinearInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) isa LinearInterpolation broken=VERSION <
+                                                                                               v"1.11" &&
+                                                                                               t isa
+                                                                                               AbstractRange
+        A = LinearInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)
 
         for (_t, _u) in zip(t, eachcol(u))
             @test A(_t) == _u
@@ -50,15 +58,23 @@ end
         @test A(0) == [0.0, 0.0]
         @test A(5.5) == [11.0, 16.5]
         @test A(11) == [22, 33]
+        @test @inferred(output_dim(A)) == 1
 
         x = 1:10
         y = 2:4
         u_ = x' .* y
         u = [u_[:, i] for i in 1:size(u_, 2)]
-        A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+        @test @inferred(LinearInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) isa LinearInterpolation broken=VERSION <
+                                                                                               v"1.11" &&
+                                                                                               t isa
+                                                                                               AbstractRange
+        A = LinearInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)
         @test A(0) == [0.0, 0.0, 0.0]
         @test A(5.5) == [11.0, 16.5, 22.0]
         @test A(11) == [22.0, 33.0, 44.0]
+        @test @inferred(output_dim(A)) == 1
     end
 
     x = 1:10
@@ -66,68 +82,76 @@ end
     u_ = x' .* y
     u = [u_[:, i:(i + 1)] for i in 1:2:10]
     t = 1.0collect(2:2:10)
+    @test_broken @inferred(LinearInterpolation(
+        u, t; extrapolation = ExtrapolationType.Extension)) isa LinearInterpolation
     A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
 
     @test A(0) == [-2.0 0.0; -3.0 0.0; -4.0 0.0]
     @test A(3) == [4.0 6.0; 6.0 9.0; 8.0 12.0]
     @test A(5) == [8.0 10.0; 12.0 15.0; 16.0 20.0]
     test_cached_index(A)
+    @test @inferred(output_dim(A)) == 2
 
     # with NaNs (#113)
     u = [NaN, 1.0, 2.0, 3.0]
     t = 1:4
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test isnan(A(1.0))
     @test A(2.0) == 1.0
     @test A(2.5) == 1.5
     @test A(3.0) == 2.0
     @test A(4.0) == 3.0
+    @test @inferred(output_dim(A)) == 0
 
     u = [0.0, NaN, 2.0, 3.0]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(1.0) == 0.0
     @test isnan(A(2.0))
     @test isnan(A(2.5))
     @test A(3.0) == 2.0
     @test A(4.0) == 3.0
+    @test @inferred(output_dim(A)) == 0
 
     u = [0.0, 1.0, NaN, 3.0]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test isnan(A(2.5))
     @test isnan(A(3.0))
     @test A(4.0) == 3.0
+    @test @inferred(output_dim(A)) == 0
 
     u = [0.0, 1.0, 2.0, NaN]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test A(3.0) == 2.0
     @test isnan(A(3.5))
     @test isnan(A(4.0))
+    @test @inferred(output_dim(A)) == 0
 
     u = [0.0, 1.0, 2.0, NaN]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test A(3.0) == 2.0
     @test isnan(A(3.5))
     @test isnan(A(4.0))
+    @test @inferred(output_dim(A)) == 0
 
     # Test type stability
     u = Float32.(1:5)
     t = Float32.(1:5)
-    A1 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A1 = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     u = 1:5
     t = 1:5
-    A2 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A2 = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     u = [1 // i for i in 1:5]
     t = (1:5)
-    A3 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A3 = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     u = [1 // i for i in 1:5]
     t = [1 // (6 - i) for i in 1:5]
-    A4 = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A4 = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
 
     F32 = Float32(1)
     F64 = Float64(1)
@@ -147,21 +171,24 @@ end
     # NaN time value for Unitful arrays: issue #365
     t = (0:3)u"s" # Unitful quantities  
     u = [0, -2, -1, -2]u"m"
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test isnan(A(NaN * u"s"))
+    @test @inferred(output_dim(A)) == 0
 
     # Nan time value:
     t = 0.0:3  # Floats
     u = [0, -2, -1, -2]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     dA = t -> ForwardDiff.derivative(A, t)
     @test isnan(dA(NaN))
+    @test @inferred(output_dim(A)) == 0
 
     t = 0:3  # Integers
     u = [0, -2, -1, -2]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     dA = t -> ForwardDiff.derivative(A, t)
     @test isnan(dA(NaN))
+    @test @inferred(output_dim(A)) == 0
 
     # Test derivative at point gives derivative to the right (except last is to left):
     ts = t[begin:(end - 1)]
@@ -172,25 +199,26 @@ end
     # Test array-valued interpolation
     u = collect.(2.0collect(1:10))
     t = 1.0collect(1:10)
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(0) == fill(0.0)
     @test A(5.5) == fill(11.0)
     @test A(11) == fill(22)
+    @test @inferred(output_dim(A)) == 0 # values are 0-dimensional arrays!
 
     # Test constant -Inf interpolation
     u = [-Inf, -Inf]
     t = [0.0, 1.0]
-    A = LinearInterpolation(u, t)
+    A = @inferred(LinearInterpolation(u, t))
     @test A(0.0) == -Inf
     @test A(0.5) == -Inf
 
     # Test extrapolation
     u = 2.0collect(1:10)
     t = 1.0collect(1:10)
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-1.0) == -2.0
     @test A(11.0) == 22.0
-    A = LinearInterpolation(u, t)
+    A = @inferred(LinearInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
     @test_throws DataInterpolations.LeftExtrapolationError A([-1.0, 11.0])
@@ -201,7 +229,7 @@ end
 
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
 
     for (_t, _u) in zip(t, u)
         @test A(_t) == _u
@@ -212,12 +240,13 @@ end
     @test A(3.5) == 12.25
     @test A(5.0) == 25
     test_cached_index(A)
+    @test @inferred(output_dim(A)) == 0
 
     # backward-looking interpolation
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    A = QuadraticInterpolation(
-        u, t, :Backward; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticInterpolation(
+        u, t, :Backward; extrapolation = ExtrapolationType.Extension))
 
     for (_t, _u) in zip(t, u)
         @test A(_t) == _u
@@ -232,8 +261,8 @@ end
     # Test both forward and backward-looking quadratic interpolation
     u = [1.0, 4.5, 6.0, 2.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    A_f = QuadraticInterpolation(u, t, :Forward)
-    A_b = QuadraticInterpolation(u, t, :Backward)
+    A_f = @inferred(QuadraticInterpolation(u, t, :Forward))
+    A_b = @inferred(QuadraticInterpolation(u, t, :Backward))
 
     for (_t, _u) in zip(t, u)
         @test A_f(_t) == _u
@@ -255,6 +284,9 @@ end
 
     # Matrix interpolation test
     u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
+    @test @inferred(QuadraticInterpolation(
+        u, t; extrapolation = ExtrapolationType.Extension)) isa QuadraticInterpolation broken=VERSION <
+                                                                                              v"1.11"
     A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
 
     for (_t, _u) in zip(t, eachcol(u))
@@ -265,31 +297,36 @@ end
     @test A(2.5) == [6.25, 6.25]
     @test A(3.5) == [12.25, 12.25]
     @test A(5.0) == [25.0, 25.0]
+    @test @inferred(output_dim(A)) == 1
 
     u_ = [1.0, 4.0, 9.0, 16.0]' .* ones(5)
     u = [u_[:, i] for i in 1:size(u_, 2)]
-    A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(0) == zeros(5)
     @test A(1.5) == 2.25 * ones(5)
     @test A(2.5) == 6.25 * ones(5)
     @test A(3.5) == 12.25 * ones(5)
     @test A(5.0) == 25.0 * ones(5)
+    @test @inferred(output_dim(A)) == 1
 
     u = [repeat(u[i], 1, 3) for i in 1:4]
+    @test_broken @inferred(QuadraticInterpolation(
+        u, t; extrapolation = ExtrapolationType.Extension)) isa QuadraticInterpolation
     A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
     @test A(0) == zeros(5, 3)
     @test A(1.5) == 2.25 * ones(5, 3)
     @test A(2.5) == 6.25 * ones(5, 3)
     @test A(3.5) == 12.25 * ones(5, 3)
     @test A(5.0) == 25.0 * ones(5, 3)
+    @test @inferred(output_dim(A)) == 2
 
     # Test extrapolation
     u = [1.0, 4.5, 6.0, 2.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    A = QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(0.0) == -4.5
     @test A(5.0) == -7.5
-    A = QuadraticInterpolation(u, t)
+    A = @inferred(QuadraticInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(0.0)
     @test_throws DataInterpolations.RightExtrapolationError A(5.0)
 end
@@ -299,57 +336,62 @@ end
 
     u = [1.0, 4.0, 9.0]
     t = [1.0, 2.0, 3.0]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == 4.0
     @test A(1.5) == 2.25
+    @test @inferred(output_dim(A)) == 0
 
     u = [1.0, 8.0, 27.0, 64.0]
     t = [1.0, 2.0, 3.0, 4.0]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == 8.0
     @test A(1.5) ≈ 3.375
     @test A(3.5) ≈ 42.875
 
     u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == [4.0, 4.0]
     @test A(1.5) ≈ [2.25, 2.25]
     @test A(3.5) ≈ [12.25, 12.25]
+    @test @inferred(output_dim(A)) == 1
 
     u_ = [1.0, 4.0, 9.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
     t = [1.0, 2.0, 3.0]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == 4.0 * ones(4)
     @test A(1.5) == 2.25 * ones(4)
+    @test @inferred(output_dim(A)) == 1
 
     u_ = [1.0, 8.0, 27.0, 64.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
     t = [1.0, 2.0, 3.0, 4.0]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == 8.0 * ones(4)
     @test A(1.5) ≈ 3.375 * ones(4)
     @test A(3.5) ≈ 42.875 * ones(4)
+    @test @inferred(output_dim(A)) == 1
 
     u = [repeat(u[i], 1, 3) for i in 1:4]
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
 
     @test A(2.0) == 8.0 * ones(4, 3)
     @test A(1.5) ≈ 3.375 * ones(4, 3)
     @test A(3.5) ≈ 42.875 * ones(4, 3)
+    @test @inferred(output_dim(A)) == 2
 
     # Test extrapolation
     u = [1.0, 4.0, 9.0]
     t = [1.0, 2.0, 3.0]
-    A = LagrangeInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LagrangeInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(0.0) == 0.0
     @test A(4.0) == 16.0
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(4.0)
 end
@@ -359,7 +401,7 @@ end
 
     u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
     t = collect(0.0:10.0)
-    A = AkimaInterpolation(u, t)
+    A = @inferred(AkimaInterpolation(u, t))
 
     @test A(0.0) ≈ 0.0
     @test A(0.5) ≈ 1.375
@@ -375,12 +417,13 @@ end
     @test A(9.9) ≈ 3.4110386597938129327189927
     @test A(10.0) ≈ 3.0
     test_cached_index(A)
+    @test @inferred(output_dim(A)) == 0
 
     # Test extrapolation
-    A = AkimaInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(AkimaInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-1.0) ≈ -5.0
     @test A(11.0) ≈ -3.924742268041234
-    A = AkimaInterpolation(u, t)
+    A = @inferred(AkimaInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
 end
@@ -391,8 +434,8 @@ end
     t = [1.0, 2.0, 3.0, 4.0]
 
     @testset "Vector case" for u in [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]]
-        A = ConstantInterpolation(
-            u, t, dir = :right; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(ConstantInterpolation(
+            u, t, dir = :right; extrapolation = ExtrapolationType.Extension))
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -403,8 +446,10 @@ end
         @test A(4.0) == u[1]
         @test A(4.5) == u[1]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 0
 
-        A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension) # dir=:left is default
+        A = @inferred(ConstantInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) # dir=:left is default
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -415,14 +460,15 @@ end
         @test A(4.0) == u[1]
         @test A(4.5) == u[1]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 0
     end
 
     @testset "Matrix case" for u in [
         [1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0],
         ["B" "C" "A" "B"; "B" "C" "A" "B"]
     ]
-        A = ConstantInterpolation(
-            u, t, dir = :right; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(ConstantInterpolation(
+            u, t, dir = :right; extrapolation = ExtrapolationType.Extension))
         @test A(0.5) == u[:, 1]
         @test A(1.0) == u[:, 1]
         @test A(1.5) == u[:, 2]
@@ -433,8 +479,10 @@ end
         @test A(4.0) == u[:, 1]
         @test A(4.5) == u[:, 1]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 1
 
-        A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension) # dir=:left is default
+        A = @inferred(ConstantInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) # dir=:left is default
         @test A(0.5) == u[:, 1]
         @test A(1.0) == u[:, 1]
         @test A(1.5) == u[:, 1]
@@ -445,13 +493,14 @@ end
         @test A(4.0) == u[:, 1]
         @test A(4.5) == u[:, 1]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 1
     end
 
     @testset "Vector of Vectors case" for u in [
         [[1.0, 2.0], [0.0, 1.0], [1.0, 2.0], [0.0, 1.0]],
         [["B", "C"], ["A", "B"], ["B", "C"], ["A", "B"]]]
-        A = ConstantInterpolation(
-            u, t, dir = :right; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(ConstantInterpolation(
+            u, t, dir = :right; extrapolation = ExtrapolationType.Extension))
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -462,8 +511,10 @@ end
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 1
 
-        A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension) # dir=:left is default
+        A = @inferred(ConstantInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) # dir=:left is default
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -474,13 +525,14 @@ end
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 1
     end
 
     @testset "Vector of Matrices case" for u in [
         [[1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0], [1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0]],
         [["B" "C"; "B" "C"], ["A" "B"; "A" "B"], ["B" "C"; "B" "C"], ["A" "B"; "A" "B"]]]
-        A = ConstantInterpolation(
-            u, t, dir = :right; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(ConstantInterpolation(
+            u, t, dir = :right; extrapolation = ExtrapolationType.Extension))
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -491,8 +543,10 @@ end
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 2
 
-        A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension) # dir=:left is default
+        A = @inferred(ConstantInterpolation(
+            u, t; extrapolation = ExtrapolationType.Extension)) # dir=:left is default
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -503,33 +557,35 @@ end
         @test A(4.0) == u[4]
         @test A(4.5) == u[4]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 2
     end
 
     # Test extrapolation
     u = [1.0, 2.0, 0.0, 1.0]
-    A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-1.0) == 1.0
     @test A(11.0) == 1.0
-    A = ConstantInterpolation(u, t)
+    A = @inferred(ConstantInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
 
     # Test extrapolation with infs with regularly spaced t
     u = [1.67e7, 1.6867e7, 1.7034e7, 1.7201e7, 1.7368e7]
     t = [0.0, 0.1, 0.2, 0.3, 0.4]
-    A = ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(ConstantInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(Inf) == last(u)
     @test A(-Inf) == first(u)
 
     # Test extrapolation of integer output
-    itp = ConstantInterpolation([2], [0.0]; extrapolation = ExtrapolationType.Constant)
+    itp = @inferred(ConstantInterpolation(
+        [2], [0.0]; extrapolation = ExtrapolationType.Constant))
     @test itp(1.0) === 2
     @test itp(-1.0) === 2
 
     # Test output type of vector evaluation (issue #388)
     u = [2, 3]
     t = [0.0, 1.0]
-    itp = ConstantInterpolation(u, t)
+    itp = @inferred(ConstantInterpolation(u, t))
     @test @inferred(itp(t)) == itp.(t)
     @test typeof(itp(t)) === typeof(itp.(t)) === Vector{Int}
 end
@@ -540,7 +596,7 @@ end
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
 
-    A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension))
 
     # Solution
     P₁ = x -> 0.5 * (x + 1) * (x + 2)
@@ -553,6 +609,7 @@ end
     @test A(0.7) == P₁(0.7)
     @test A(2.0) == P₁(2.0)
     test_cached_index(A)
+    @test @inferred(output_dim(A)) == 0
 
     u_ = [0.0, 1.0, 3.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
@@ -561,21 +618,23 @@ end
     @test A(-0.5) == P₁(-0.5) * ones(4)
     @test A(0.7) == P₁(0.7) * ones(4)
     @test A(2.0) == P₁(2.0) * ones(4)
+    @test @inferred(output_dim(A)) == 1
 
     u = [repeat(u[i], 1, 3) for i in 1:3]
-    A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-2.0) == P₁(-2.0) * ones(4, 3)
     @test A(-0.5) == P₁(-0.5) * ones(4, 3)
     @test A(0.7) == P₁(0.7) * ones(4, 3)
     @test A(2.0) == P₁(2.0) * ones(4, 3)
+    @test @inferred(output_dim(A)) == 2
 
     # Test extrapolation
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-2.0) == 0.0
     @test A(2.0) == 6.0
-    A = QuadraticSpline(u, t)
+    A = @inferred(QuadraticSpline(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-2.0)
     @test_throws DataInterpolations.RightExtrapolationError A(2.0)
 end
@@ -586,7 +645,7 @@ end
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
 
-    A = CubicSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(CubicSpline(u, t; extrapolation = ExtrapolationType.Extension))
     test_cached_index(A)
 
     # Solution
@@ -602,9 +661,12 @@ end
     for x in (0.3, 0.5, 1.5)
         @test A(x) ≈ P₂(x)
     end
+    @test @inferred(output_dim(A)) == 0
 
     u_ = [0.0, 1.0, 3.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
+    @test @inferred(CubicSpline(u, t; extrapolation = ExtrapolationType.Extension)) isa
+          CubicSpline broken=VERSION < v"1.11"
     A = CubicSpline(u, t; extrapolation = ExtrapolationType.Extension)
     for x in (-1.5, -0.5, -0.7)
         @test A(x) ≈ P₁(x) * ones(4)
@@ -612,8 +674,11 @@ end
     for x in (0.3, 0.5, 1.5)
         @test A(x) ≈ P₂(x) * ones(4)
     end
+    @test @inferred(output_dim(A)) == 1
 
     u = [repeat(u[i], 1, 3) for i in 1:3]
+    @test_broken @inferred(CubicSpline(
+        u, t; extrapolation = ExtrapolationType.Extension)) isa CubicSpline
     A = CubicSpline(u, t; extrapolation = ExtrapolationType.Extension)
     for x in (-1.5, -0.5, -0.7)
         @test A(x) ≈ P₁(x) * ones(4, 3)
@@ -621,20 +686,22 @@ end
     for x in (0.3, 0.5, 1.5)
         @test A(x) ≈ P₂(x) * ones(4, 3)
     end
+    @test @inferred(output_dim(A)) == 2
 
     # Test extrapolation
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    A = CubicSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(CubicSpline(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-2.0) ≈ -1.0
     @test A(2.0) ≈ 5.0
-    A = CubicSpline(u, t)
+    A = @inferred(CubicSpline(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-2.0)
     @test_throws DataInterpolations.RightExtrapolationError A(2.0)
 
     @testset "AbstractMatrix" begin
         t = 0.1:0.1:1.0
         u = [sin.(t) cos.(t)]' |> collect
+        @test_broken @inferred(CubicSpline(u, t)) isa CubicSpline
         c = CubicSpline(u, t)
         t_test = 0.1:0.05:1.0
         u_test = reduce(hcat, c.(t_test))
@@ -646,6 +713,7 @@ end
                   0.0 cos(2t)]
         t = 0.1:0.1:1.0
         u3d = f3d.(t)
+        @test_broken @inferred(CubicSpline(u3d, t)) isa CubicSpline
         c = CubicSpline(u3d, t)
         t_test = 0.1:0.05:1.0
         u_test = reduce(hcat, c.(t_test))
@@ -661,23 +729,24 @@ end
         t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
         u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
         test_interpolation_type(BSplineInterpolation)
-        A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)
+        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform, :Uniform))
 
         @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
         @test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]
         @test [A(t[1]), A(t[end])] == [u[1], u[end]]
         test_cached_index(A)
+        @test @inferred(output_dim(A)) == 0
 
         # Test extrapolation
-        A = BSplineInterpolation(
-            u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(BSplineInterpolation(
+            u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension))
         @test A(-1.0) == u[1]
         @test A(300.0) == u[end]
-        A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)
+        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform, :Uniform))
         @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
         @test_throws DataInterpolations.RightExtrapolationError A(300.0)
 
-        A = BSplineInterpolation(u, t, 2, :ArcLen, :Average)
+        A = @inferred(BSplineInterpolation(u, t, 2, :ArcLen, :Average))
 
         @test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609]
         @test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]
@@ -690,43 +759,47 @@ end
         @test_nowarn BSplineInterpolation(u[1:3], t[1:3], 2, :Uniform, :Uniform)
 
         # Test extrapolation
-        A = BSplineInterpolation(
-            u, t, 2, :ArcLen, :Average; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(BSplineInterpolation(
+            u, t, 2, :ArcLen, :Average; extrapolation = ExtrapolationType.Extension))
         @test A(-1.0) == u[1]
         @test A(300.0) == u[end]
-        A = BSplineInterpolation(u, t, 2, :ArcLen, :Average)
+        A = @inferred(BSplineInterpolation(u, t, 2, :ArcLen, :Average))
         @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
         @test_throws DataInterpolations.RightExtrapolationError A(300.0)
 
         @testset "AbstractMatrix" begin
             t = 0.1:0.1:1.0
             u2d = [sin.(t) cos.(t)]' |> collect
-            A = BSplineInterpolation(u2d, t, 2, :Uniform, :Uniform)
+            A = @inferred(BSplineInterpolation(u2d, t, 2, :Uniform, :Uniform))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3)
             @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3)
+            @test @inferred(output_dim(A)) == 1
 
-            A = BSplineInterpolation(u2d, t, 2, :ArcLen, :Average)
+            A = @inferred(BSplineInterpolation(u2d, t, 2, :ArcLen, :Average))
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3)
             @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3)
+            @test @inferred(output_dim(A)) == 1
         end
         @testset "AbstractArray{T, 3}" begin
             f3d(t) = [sin(t) cos(t);
                       0.0 cos(2t)]
             t = 0.1:0.1:1.0
             u3d = cat(f3d.(t)..., dims = 3)
-            A = BSplineInterpolation(u3d, t, 2, :Uniform, :Uniform)
+            A = @inferred(BSplineInterpolation(u3d, t, 2, :Uniform, :Uniform))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             f_test = reduce(hcat, f3d.(t_test))
             @test isapprox(u_test, f_test, atol = 1e-2)
+            @test @inferred(output_dim(A)) == 2
 
-            A = BSplineInterpolation(u3d, t, 2, :ArcLen, :Average)
+            A = @inferred(BSplineInterpolation(u3d, t, 2, :ArcLen, :Average))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test, f_test, atol = 1e-2)
+            @test @inferred(output_dim(A)) == 2
         end
     end
 
@@ -795,10 +868,11 @@ end
     du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
-    A = CubicHermiteSpline(du, u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(CubicHermiteSpline(du, u, t; extrapolation = ExtrapolationType.Extension))
     @test A.(t) ≈ u
     @test A(100.0)≈10.106770 rtol=1e-5
     @test A(300.0)≈9.901542 rtol=1e-5
+    @test @inferred(output_dim(A)) == 0
     test_cached_index(A)
     push!(u, 1.0)
     @test_throws AssertionError CubicHermiteSpline(du, u, t)
@@ -807,13 +881,14 @@ end
 @testset "PCHIPInterpolation" begin
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 250.0]
-    A = PCHIPInterpolation(u, t)
+    A = @inferred(PCHIPInterpolation(u, t))
     @test A isa CubicHermiteSpline
     ts = 0.0:0.1:250.0
     us = A(ts)
     @test all(minimum(u) .<= us)
     @test all(maximum(u) .>= us)
     @test all(A.du[3:4] .== 0.0)
+    @test @inferred(output_dim(A)) == 0
 end
 
 @testset "Quintic Hermite Spline" begin
@@ -823,13 +898,16 @@ end
     du = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
-    A = QuinticHermiteSpline(ddu, du, u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(QuinticHermiteSpline(
+        ddu, du, u, t; extrapolation = ExtrapolationType.Extension))
     @test A.(t) ≈ u
     @test A(100.0)≈10.107996 rtol=1e-5
     @test A(300.0)≈11.364162 rtol=1e-5
     test_cached_index(A)
+    @test @inferred(output_dim(A)) == 0
     push!(u, 1.0)
     @test_throws AssertionError QuinticHermiteSpline(ddu, du, u, t)
+    @test @inferred(output_dim(A)) == 0
 end
 
 @testset "Curvefit" begin
@@ -852,6 +930,7 @@ end
     ]
     us = A.(ts)
     @test vs ≈ us
+    @test @inferred(output_dim(A)) == 0
 
     # Test extrapolation
     A = Curvefit(u, t, model, p0, LBFGS(); extrapolate = true)
@@ -865,7 +944,7 @@ end
     ut1 = Float32[0.1, 0.2, 0.3, 0.4, 0.5]
     ut2 = Float64[0.1, 0.2, 0.3, 0.4, 0.5]
     for u in (ut1, ut2), t in (ut1, ut2)
-        interp = LinearInterpolation(ut1, ut2)
+        interp = @inferred(LinearInterpolation(ut1, ut2))
         for xs in (u, t)
             ys = @inferred(interp(xs))
             @test ys isa Vector{typeof(interp(first(xs)))}