diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index bf9efcb3..e44e6315 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -132,7 +132,7 @@ struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation
             Aitp,
             extrapolation_left,
             extrapolation_right)
-        N = get_output_dim(u)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), eltype(u), typeof(λ), N, typeof(Aitp)}(
             u,
             û,
@@ -174,7 +174,7 @@ struct CurvefitCache{
     pmin::pminType  # optimized params
     extrapolate::Bool
     function CurvefitCache(u, t, m, p0, ub, lb, alg, pmin, extrapolate)
-        N = get_output_dim(u)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(m),
             typeof(p0), typeof(ub), typeof(lb),
             typeof(alg), typeof(pmin), eltype(u), N}(u,
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
index 577b98e3..efd8702c 100644
--- a/src/integral_inverses.jl
+++ b/src/integral_inverses.jl
@@ -42,7 +42,7 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <:
     iguesser::Guesser{tType}
     itp::itpType
     function LinearInterpolationIntInv(u, t, A)
-        N = get_output_dim(u)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
             u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A)
     end
@@ -94,7 +94,7 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <:
     iguesser::Guesser{tType}
     itp::itpType
     function ConstantInterpolationIntInv(u, t, A)
-        N = get_output_dim(u)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(A), eltype(u), N}(
             u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A
         )
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index aa51f262..e04ec37e 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -40,7 +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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u), N}(
             u, t, I, p, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
@@ -111,7 +111,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.α), eltype(u), N}(
             u, t, I, p, mode, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
@@ -175,7 +175,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), eltype(u), typeof(bcache), N}(u,
             t,
             n,
@@ -246,7 +246,7 @@ 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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
             typeof(d), eltype(u), N}(u,
             t,
@@ -343,7 +343,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), eltype(u), N}(
             u, t, I, nothing, dir, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
@@ -411,7 +411,7 @@ 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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k),
             typeof(c), typeof(sc), eltype(u), N}(u,
             t,
@@ -532,7 +532,7 @@ 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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁),
             typeof(h), typeof(z), eltype(u), N}(
             u,
@@ -713,7 +713,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
             u,
             t,
@@ -948,7 +948,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}(
             u,
             t,
@@ -1208,7 +1208,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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u), N}(
             du, u, t, I, p, extrapolation_left, extrapolation_right,
             Guesser(t), cache_parameters, linear_lookup)
@@ -1312,7 +1312,7 @@ 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)
+        N = output_ndims(u)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
             typeof(ddu), typeof(p.c₁), eltype(u), N}(
             ddu, du, u, t, I, p, extrapolation_left, extrapolation_right,
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index cb172b52..80988bd8 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -59,18 +59,11 @@ 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
+# Get the number of dimensions `ndims(interp(x))` of the interpolation `interp` evaluated at a single input `x`
+# It is derived from the set of values `u` at the interpolation nodes
+output_ndims(::AbstractVector) = 0 # each value is a scalar
+output_ndims(::AbstractVector{<:AbstractArray{<:Any, N}}) where {N} = N # each value is an array but values are not stacked
+output_ndims(::AbstractArray{<:Any, N}) where {N} = N - 1 # each value is an array but multiple values are stacked
 
 function quadratic_spline_params(t::AbstractVector, sc::AbstractVector)
 
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index aecde638..ce3de3ef 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -32,7 +32,9 @@ 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))
+        @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
 
         for (_t, _u) in zip(t, u)
             @test A(_t) == _u
@@ -42,7 +44,14 @@ end
         @test A(11) == 22
 
         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)
+        @test A isa DataInterpolations.AbstractInterpolation{Float64, 1}
 
         for (_t, _u) in zip(t, eachcol(u))
             @test A(_t) == _u
@@ -55,7 +64,14 @@ end
         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 isa DataInterpolations.AbstractInterpolation{Vector{Int}, 1}
         @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]
@@ -66,8 +82,10 @@ 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 isa DataInterpolations.AbstractInterpolation{Matrix{Int64}, 2}
     @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]
@@ -76,7 +94,8 @@ end
     # 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 A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test isnan(A(1.0))
     @test A(2.0) == 1.0
     @test A(2.5) == 1.5
@@ -84,7 +103,8 @@ end
     @test A(4.0) == 3.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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(1.0) == 0.0
     @test isnan(A(2.0))
     @test isnan(A(2.5))
@@ -92,7 +112,8 @@ end
     @test A(4.0) == 3.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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test isnan(A(2.5))
@@ -100,7 +121,8 @@ end
     @test A(4.0) == 3.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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test A(3.0) == 2.0
@@ -108,7 +130,8 @@ end
     @test isnan(A(4.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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(1.0) == 0.0
     @test A(2.0) == 1.0
     @test A(3.0) == 2.0
@@ -118,16 +141,16 @@ end
     # 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,19 +170,22 @@ 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 A isa DataInterpolations.AbstractInterpolation{typeof(0u"m"), 0}
     @test isnan(A(NaN * u"s"))
 
     # 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))
+    @test A isa DataInterpolations.AbstractInterpolation{Int, 0}
     dA = t -> ForwardDiff.derivative(A, t)
     @test isnan(dA(NaN))
 
     t = 0:3  # Integers
     u = [0, -2, -1, -2]
-    A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
+    @test A isa DataInterpolations.AbstractInterpolation{Int, 0}
     dA = t -> ForwardDiff.derivative(A, t)
     @test isnan(dA(NaN))
 
@@ -172,7 +198,8 @@ 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 isa DataInterpolations.AbstractInterpolation{Array{Float64, 0}, 0}
     @test A(0) == fill(0.0)
     @test A(5.5) == fill(11.0)
     @test A(11) == fill(22)
@@ -180,17 +207,20 @@ end
     # Test constant -Inf interpolation
     u = [-Inf, -Inf]
     t = [0.0, 1.0]
-    A = LinearInterpolation(u, t)
+    A = @inferred(LinearInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(-1.0) == -2.0
     @test A(11.0) == 22.0
-    A = LinearInterpolation(u, t)
+    A = @inferred(LinearInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
     @test_throws DataInterpolations.LeftExtrapolationError A([-1.0, 11.0])
@@ -201,8 +231,8 @@ 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))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     for (_t, _u) in zip(t, u)
         @test A(_t) == _u
     end
@@ -216,9 +246,9 @@ end
     # 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))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     for (_t, _u) in zip(t, u)
         @test A(_t) == _u
     end
@@ -232,9 +262,10 @@ 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))
+    @test A_f isa DataInterpolations.AbstractInterpolation{Float64, 0}
+    A_b = @inferred(QuadraticInterpolation(u, t, :Backward))
+    @test A_b isa DataInterpolations.AbstractInterpolation{Float64, 0}
     for (_t, _u) in zip(t, u)
         @test A_f(_t) == _u
         @test A_b(_t) == _u
@@ -255,8 +286,11 @@ 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)
-
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 1}
     for (_t, _u) in zip(t, eachcol(u))
         @test A(_t) == _u
     end
@@ -268,7 +302,8 @@ end
 
     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 isa DataInterpolations.AbstractInterpolation{Vector{Float64}, 1}
     @test A(0) == zeros(5)
     @test A(1.5) == 2.25 * ones(5)
     @test A(2.5) == 6.25 * ones(5)
@@ -276,7 +311,10 @@ end
     @test A(5.0) == 25.0 * ones(5)
 
     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 isa DataInterpolations.AbstractInterpolation{Matrix{Float64}, 2}
     @test A(0) == zeros(5, 3)
     @test A(1.5) == 2.25 * ones(5, 3)
     @test A(2.5) == 6.25 * ones(5, 3)
@@ -286,10 +324,12 @@ end
     # 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(0.0) == -4.5
     @test A(5.0) == -7.5
-    A = QuadraticInterpolation(u, t)
+    A = @inferred(QuadraticInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.LeftExtrapolationError A(0.0)
     @test_throws DataInterpolations.RightExtrapolationError A(5.0)
 end
@@ -299,22 +339,22 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(2.0) == 4.0
     @test A(1.5) == 2.25
 
     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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @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 isa DataInterpolations.AbstractInterpolation{Float64, 1}
     @test A(2.0) == [4.0, 4.0]
     @test A(1.5) ≈ [2.25, 2.25]
     @test A(3.5) ≈ [12.25, 12.25]
@@ -322,23 +362,23 @@ end
     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 isa DataInterpolations.AbstractInterpolation{Vector{Float64}, 1}
     @test A(2.0) == 4.0 * ones(4)
     @test A(1.5) == 2.25 * ones(4)
 
     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 isa DataInterpolations.AbstractInterpolation{Vector{Float64}, 1}
     @test A(2.0) == 8.0 * ones(4)
     @test A(1.5) ≈ 3.375 * ones(4)
     @test A(3.5) ≈ 42.875 * ones(4)
 
     u = [repeat(u[i], 1, 3) for i in 1:4]
-    A = LagrangeInterpolation(u, t)
-
+    A = @inferred(LagrangeInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Matrix{Float64}, 2}
     @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)
@@ -346,10 +386,12 @@ end
     # 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(0.0) == 0.0
     @test A(4.0) == 16.0
-    A = LagrangeInterpolation(u, t)
+    A = @inferred(LagrangeInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(4.0)
 end
@@ -359,8 +401,8 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(0.0) ≈ 0.0
     @test A(0.5) ≈ 1.375
     @test A(1.0) ≈ 2.0
@@ -377,10 +419,12 @@ end
     test_cached_index(A)
 
     # Test extrapolation
-    A = AkimaInterpolation(u, t; extrapolation = ExtrapolationType.Extension)
+    A = @inferred(AkimaInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(-1.0) ≈ -5.0
     @test A(11.0) ≈ -3.924742268041234
-    A = AkimaInterpolation(u, t)
+    A = @inferred(AkimaInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
 end
@@ -391,8 +435,9 @@ 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 isa DataInterpolations.AbstractInterpolation{eltype(u), 0}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -404,7 +449,9 @@ end
         @test A(4.5) == u[1]
         test_cached_index(A)
 
-        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 isa DataInterpolations.AbstractInterpolation{eltype(u), 0}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -421,8 +468,9 @@ end
         [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 isa DataInterpolations.AbstractInterpolation{eltype(u), 1}
         @test A(0.5) == u[:, 1]
         @test A(1.0) == u[:, 1]
         @test A(1.5) == u[:, 2]
@@ -434,7 +482,9 @@ end
         @test A(4.5) == u[:, 1]
         test_cached_index(A)
 
-        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 isa DataInterpolations.AbstractInterpolation{eltype(u), 1}
         @test A(0.5) == u[:, 1]
         @test A(1.0) == u[:, 1]
         @test A(1.5) == u[:, 1]
@@ -450,8 +500,9 @@ 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 isa DataInterpolations.AbstractInterpolation{eltype(u), 1}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -463,7 +514,9 @@ end
         @test A(4.5) == u[4]
         test_cached_index(A)
 
-        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 isa DataInterpolations.AbstractInterpolation{eltype(u), 1}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -479,8 +532,9 @@ 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 isa DataInterpolations.AbstractInterpolation{eltype(u), 2}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[2]
@@ -492,7 +546,9 @@ end
         @test A(4.5) == u[4]
         test_cached_index(A)
 
-        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 isa DataInterpolations.AbstractInterpolation{eltype(u), 2}
         @test A(0.5) == u[1]
         @test A(1.0) == u[1]
         @test A(1.5) == u[1]
@@ -507,29 +563,35 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(-1.0) == 1.0
     @test A(11.0) == 1.0
-    A = ConstantInterpolation(u, t)
+    A = @inferred(ConstantInterpolation(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @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 isa DataInterpolations.AbstractInterpolation{Int, 0}
     @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 itp isa DataInterpolations.AbstractInterpolation{Int, 0}
     @test @inferred(itp(t)) == itp.(t)
     @test typeof(itp(t)) === typeof(itp.(t)) === Vector{Int}
 end
@@ -540,8 +602,8 @@ 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))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     # Solution
     P₁ = x -> 0.5 * (x + 1) * (x + 2)
 
@@ -557,13 +619,15 @@ end
     u_ = [0.0, 1.0, 3.0]' .* ones(4)
     u = [u_[:, i] for i in 1:size(u_, 2)]
     A = QuadraticSpline(u, t; extrapolation = ExtrapolationType.Extension)
+    @test A isa DataInterpolations.AbstractInterpolation{Vector{Float64}, 1}
     @test A(-2.0) == P₁(-2.0) * ones(4)
     @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)
 
     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 isa DataInterpolations.AbstractInterpolation{Matrix{Float64}, 2}
     @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)
@@ -572,10 +636,12 @@ end
     # 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(-2.0) == 0.0
     @test A(2.0) == 6.0
-    A = QuadraticSpline(u, t)
+    A = @inferred(QuadraticSpline(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.LeftExtrapolationError A(-2.0)
     @test_throws DataInterpolations.RightExtrapolationError A(2.0)
 end
@@ -586,7 +652,8 @@ 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 A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     test_cached_index(A)
 
     # Solution
@@ -605,7 +672,10 @@ end
 
     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)
+    @test A isa DataInterpolations.AbstractInterpolation{Vector{Float64}, 1}
     for x in (-1.5, -0.5, -0.7)
         @test A(x) ≈ P₁(x) * ones(4)
     end
@@ -614,7 +684,10 @@ end
     end
 
     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)
+    @test A isa DataInterpolations.AbstractInterpolation{Matrix{Float64}, 2}
     for x in (-1.5, -0.5, -0.7)
         @test A(x) ≈ P₁(x) * ones(4, 3)
     end
@@ -625,17 +698,21 @@ end
     # 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(-2.0) ≈ -1.0
     @test A(2.0) ≈ 5.0
-    A = CubicSpline(u, t)
+    A = @inferred(CubicSpline(u, t))
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @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)
+        @test c isa DataInterpolations.AbstractInterpolation{Float64, 1}
         t_test = 0.1:0.05:1.0
         u_test = reduce(hcat, c.(t_test))
         @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3)
@@ -646,7 +723,9 @@ 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)
+        @test c isa DataInterpolations.AbstractInterpolation{Matrix{Float64}, 2}
         t_test = 0.1:0.05:1.0
         u_test = reduce(hcat, c.(t_test))
         f_test = reduce(hcat, f3d.(t_test))
@@ -661,24 +740,26 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @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 extrapolation
-        A = BSplineInterpolation(
-            u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension)
+        A = @inferred(BSplineInterpolation(
+            u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension))
+        @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @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 A isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609]
         @test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]
         @test [A(t[1]), A(t[end])] == [u[1], u[end]]
@@ -690,24 +771,28 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @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 A isa DataInterpolations.AbstractInterpolation{Float64, 0}
         @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))
+            @test A isa DataInterpolations.AbstractInterpolation{Float64, 1}
             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)
 
-            A = BSplineInterpolation(u2d, t, 2, :ArcLen, :Average)
+            A = @inferred(BSplineInterpolation(u2d, t, 2, :ArcLen, :Average))
+            @test A isa DataInterpolations.AbstractInterpolation{Float64, 1}
             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)
@@ -717,13 +802,15 @@ end
                       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))
+            @test A isa DataInterpolations.AbstractInterpolation{Float64, 2}
             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)
 
-            A = BSplineInterpolation(u3d, t, 2, :ArcLen, :Average)
+            A = @inferred(BSplineInterpolation(u3d, t, 2, :ArcLen, :Average))
+            @test A isa DataInterpolations.AbstractInterpolation{Float64, 2}
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test, f_test, atol = 1e-2)
@@ -795,7 +882,8 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A.(t) ≈ u
     @test A(100.0)≈10.106770 rtol=1e-5
     @test A(300.0)≈9.901542 rtol=1e-5
@@ -807,8 +895,9 @@ 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
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     ts = 0.0:0.1:250.0
     us = A(ts)
     @test all(minimum(u) .<= us)
@@ -823,7 +912,9 @@ 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 isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A.(t) ≈ u
     @test A(100.0)≈10.107996 rtol=1e-5
     @test A(300.0)≈11.364162 rtol=1e-5
@@ -841,7 +932,7 @@ end
     p0 = [0.5, 0.5]
 
     A = Curvefit(u, t, model, p0, LBFGS())
-
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     ts = [-7.0, -2.0, 0.0, 2.5, 5.0]
     vs = [
         1.0013468217936277,
@@ -855,8 +946,10 @@ end
 
     # Test extrapolation
     A = Curvefit(u, t, model, p0, LBFGS(); extrapolate = true)
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test A(15.0) == model(15.0, A.pmin)
     A = Curvefit(u, t, model, p0, LBFGS())
+    @test A isa DataInterpolations.AbstractInterpolation{Float64, 0}
     @test_throws DataInterpolations.ExtrapolationError A(15.0)
 end
 
@@ -865,7 +958,8 @@ 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))
+        @test interp isa DataInterpolations.AbstractInterpolation{Float32, 0}
         for xs in (u, t)
             ys = @inferred(interp(xs))
             @test ys isa Vector{typeof(interp(first(xs)))}