SciML · rafaqz · Oct 16, 2025 · Oct 16, 2025 · Oct 16, 2025 · Oct 16, 2025
diff --git a/Project.toml b/Project.toml
@@ -18,6 +18,7 @@ DataInterpolationsNDSymbolicsExt = ["SymbolicUtils", "Symbolics"]
 
 [compat]
 Adapt = "4.3.0"
+Aqua = "0.8"
 DataInterpolations = "8"
 EllipsisNotation = "1.8.0"
 ForwardDiff = "0"
@@ -32,6 +33,7 @@ Test = "1"
 julia = "1"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -42,4 +44,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DataInterpolations", "ForwardDiff", "Pkg", "Random", "SafeTestsets", "Symbolics", "Test", "SymbolicUtils"]
+test = ["Aqua", "DataInterpolations", "ForwardDiff", "Pkg", "Random", "SafeTestsets", "Symbolics", "Test", "SymbolicUtils"]
diff --git a/README.md b/README.md
@@ -14,7 +14,7 @@ DataInterpolationsND.jl is a library for interpolating arbitrarily high dimensio
 
 For one dimensional interpolation see also [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl).
 
-DataInterpolationsND.jl is similar in functionality to [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), a well established interpolation package that is currently much more feature rich than DataInterpolationsND.jl. We hope to justify the existence of DataInterpolationsND.jl through its use of the KernelAbstractions and its planned integration with the SciML ecosystem.
+DataInterpolationsND.jl is similar in functionality to [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), a well established interpolation package that is currently much more feature rich than DataInterpolationsND.jl. We hope to justify the existence of DataInterpolationsND.jl through its use of KernelAbstractions and its planned integration with the SciML ecosystem.
 
 ## API
 
@@ -28,7 +28,7 @@ t2 = cumsum(rand(7))
 
 interpolation_dimensions = (
     LinearInterpolationDimension(t1),
-    LinearInterpolationDimension(t2)
+    ConstantInterpolationDimension(t2)
 )
 
 # The outputs will be vectors of length 2
@@ -51,7 +51,7 @@ If we provide `t_eval` for the interpolation dimensions, we can evaluate at thes
 ```julia
 interpolation_dimensions = (
     LinearInterpolationDimension(t1; t_eval = range(first(t1), last(t1); length = 100)),
-    LinearInterpolationDimension(t2; t_eval = range(first(t2), last(t2); length = 100))
+    ConstantInterpolationDimension(t2; t_eval = range(first(t2), last(t2); length = 100))
 )
 
 interp = NDInterpolation(u, interpolation_dimensions)

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -16,7 +16,7 @@ t2 = cumsum(rand(7))
 
 interpolation_dimensions = (
     LinearInterpolationDimension(t1),
-    LinearInterpolationDimension(t2)
+    ConstantInterpolationDimension(t2)
 )
 
 # The outputs will be vectors of length 2
@@ -39,7 +39,7 @@ If we provide `t_eval` for the interpolation dimensions, we can evaluate at thes
 ```julia
 interpolation_dimensions = (
     LinearInterpolationDimension(t1; t_eval = range(first(t1), last(t1); length = 100)),
-    LinearInterpolationDimension(t2; t_eval = range(first(t2), last(t2); length = 100))
+    ConstantInterpolationDimension(t2; t_eval = range(first(t2), last(t2); length = 100))
 )
 
 interp = NDInterpolation(u, interpolation_dimensions)

diff --git a/docs/src/usage.md b/docs/src/usage.md
@@ -9,7 +9,7 @@ Random.seed!(1)
 
 interp_dims = (
     LinearInterpolationDimension(cumsum(0.5 .+ rand(5))),
-    LinearInterpolationDimension(cumsum(0.5 .+ rand(10)))
+    ConstantInterpolationDimension(cumsum(0.5 .+ rand(10)))
 )
 
 t_eval_1 = range(
@@ -85,7 +85,7 @@ using LinearAlgebra
 
 interp_dims = (
     LinearInterpolationDimension(interp_dims[1].t; t_eval = t_eval_1),
-    LinearInterpolationDimension(interp_dims[2].t; t_eval = t_eval_2)
+    ConstantInterpolationDimension(interp_dims[2].t; t_eval = t_eval_2)
 )
 itp = NDInterpolation(u, interp_dims)
 

diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl
@@ -13,48 +13,57 @@ struct EmptyCache <: AbstractInterpolationCache end
 """
     NDInterpolation(u, interp_dims; cache = EmptyCache())
 
-The interpolation object containing the interpolation dimensions and the data to interpolate `u`.
-Given the number of interpolation dimensions `N_in`, for first `N_in` dimensions of `u`
-the size of `u` along that dimension must match the length of `t` of the corresponding interpolation dimension.
+The interpolation object containing the interpolation dimensions and the data to be interpolated `u`.
+For non-trivial interpolation dimensions (i.e. not `NoInterpolationDimension`) the length of its `t` must match
+the size of `u` along the corresponding dimension.
 
 ## Arguments
 
-  - `interp_dims`: A tuple of identically typed interpolation dimensions.
+  - `interp_dims`: A tuple of interpolation dimensions of length `≤ ndims(u)`.
   - `u`: The array to be interpolated.
 
 ## Keyword Arguments
 
   - `cache`: Optional global cache
 """
 struct NDInterpolation{
-        N_in, N_out,
-        ID <: AbstractInterpolationDimension,
+        N,
+        N_in,
+        N_out,
         gType <: AbstractInterpolationCache,
+        D <: Tuple{Vararg{<:AbstractInterpolationDimension}},
         uType <: AbstractArray,
     }
     u::uType
-    interp_dims::NTuple{N_in, ID}
+    interp_dims::D
     cache::gType
-    function NDInterpolation(u, interp_dims, cache)
+    function NDInterpolation(u::AbstractArray{<:Any, N}, interp_dims, cache) where {N}
         if interp_dims isa AbstractInterpolationDimension
             interp_dims = (interp_dims,)
         end
-        N_in = length(interp_dims)
-        N_out = ndims(u) - N_in
-        @assert N_out ≥ 0 "The number of dimensions of u must be at least the number of interpolation dimensions."
+        @assert length(interp_dims) <= N "The number of interpolation dimensions must be no more than the dimensions of u. got $(length(interp_dims)) and $N"
+        interp_dims = _add_trailing_no_interp_dims(interp_dims, Val{N}())
+        N_in = length(remove(NoInterpolationDimension, interp_dims))
+        N_out = length(keep(NoInterpolationDimension, interp_dims))
         validate_size_u(interp_dims, u)
         validate_cache(cache, interp_dims, u)
-        return new{N_in, N_out, eltype(interp_dims), typeof(cache), typeof(u)}(
+        return new{N, N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}(
             u, interp_dims, cache
         )
     end
 end
-
 # Constructor with optional global cache
 function NDInterpolation(u, interp_dims; cache = EmptyCache())
     return NDInterpolation(u, interp_dims, cache)
 end
 
+function _add_trailing_no_interp_dims(dim::AbstractInterpolationDimension, n)
+    return _add_trailing_interp_dims((dim,), n)
+end
+function _add_trailing_no_interp_dims(dims::Tuple, ::Val{N}) where {N}
+    return (dims..., ntuple(_ -> NoInterpolationDimension(), Val{N - length(dims)}())...)
+end
+
 @adapt_structure NDInterpolation
 
 include("interpolation_dimensions.jl")
@@ -76,15 +85,22 @@ function (interp::NDInterpolation)(
 end
 
 # In place single input evaluation
-function (interp::NDInterpolation{N_in})(
-        out::Union{Number, AbstractArray{<:Number}},
-        t::Tuple{Vararg{Number, N_in}};
-        derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in)
-    ) where {N_in}
-    validate_derivative_orders(derivative_orders, interp)
-    idx = get_idx(interp.interp_dims, t)
-    @assert size(out) == size(interp.u)[(N_in + 1):end] "The size of out must match the size of the last N_out dimensions of u."
-    return _interpolate!(out, interp, t, idx, derivative_orders, nothing)
+function (interp::NDInterpolation{N, N_in, N_out})(
+        out::Union{Number, AbstractArray{<:Number, N_out}},
+        t::Tuple{Vararg{Any, N_in}};
+        derivative_orders::Tuple{Vararg{Integer, N_in}} = ntuple(_ -> 0, N_in)
+    ) where {N, N_in, N_out}
+    (; interp_dims, u) = interp
+    # We need to add the NoInterpolationDimensions for all arguments to _interpolate.
+    # Currently interp() accepts only the interpolating dimensions.
+    # We could switch this so it needs indices for NoInterpolationDimension, but thats a breaking change.
+    d_o_all = insertat(NoInterpolationDimension, 0, derivative_orders, interp_dims)
+    validate_derivative_order(d_o_all, interp)
+    t_all = insertat(NoInterpolationDimension, Colon(), t, interp.interp_dims)
+    idx_all = get_idx(interp_dims, t_all)
+    @assert size(out) == keepat(NoInterpolationDimension, size(u), interp_dims) "The size of out must match the size of the NoInterpolationDimensions of u."
+    multi_point_index = nothing
+    return _interpolate!(out, interp, t_all, idx_all, d_o_all, multi_point_index)
 end
 
 # Out of place single input evaluation
@@ -94,9 +110,7 @@ function (interp::NDInterpolation)(t::Tuple{Vararg{Number}}; kwargs...)
 end
 
 export NDInterpolation, LinearInterpolationDimension, ConstantInterpolationDimension,
-    BSplineInterpolationDimension, NURBSWeights,
+    BSplineInterpolationDimension, NURBSWeights, NoInterpolationDimension,
     eval_unstructured, eval_unstructured!, eval_grid, eval_grid!
 
-include("precompilation.jl")
-
 end # module DataInterpolationsND
diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl
@@ -1,3 +1,10 @@
+"""
+    NoInterpolationDimension()
+
+A dimension that does not perform interpolation.
+"""
+struct NoInterpolationDimension <: AbstractInterpolationDimension end
+
 """
     LinearInterpolationDimension(t; t_eval = similar(t, 0))
 
@@ -161,15 +168,9 @@ function BSplineInterpolationDimension(
     synchronize(backend)
 
     idx_eval = similar(t_eval, Int)
-    basis_function_eval = similar(
-        t_eval,
-        typeof(inv(one(eltype(t))) * inv(one(eltype(t_eval)))),
-        (
-            length(t_eval),
-            degree + 1,
-            max_derivative_order_eval + 1,
-        )
-    )
+    T = typeof(inv(one(eltype(t))) * inv(one(eltype(t_eval))))
+    s = (length(t_eval), degree + 1, max_derivative_order_eval + 1)
+    basis_function_eval = similar(t_eval, T, s)
     itp_dim = BSplineInterpolationDimension(
         t, knots_all, t_eval, idx_eval, degree, max_derivative_order_eval,
         basis_function_eval, multiplicities