From d37c715f68680f2bfee5231d73c9057793c237a6 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 16 Oct 2025 17:12:20 +1100 Subject: [PATCH 01/32] gneralise interpolate to allow mixed dimensions --- src/interpolation_dimensions.jl | 7 ++ src/interpolation_methods.jl | 125 +++++++++++++------------------- 2 files changed, 57 insertions(+), 75 deletions(-) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 4cf7295..1eb6b69 100644 --- 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)) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 19c27f3..80f2552 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,34 +1,32 @@ function _interpolate!( out, A::NDInterpolation{N_in, N_out, ID}, - t::Tuple{Vararg{Number, N_in}}, + ts::Tuple{Vararg{Number, N_in}}, idx::NTuple{N_in, <:Integer}, derivative_orders::NTuple{N_in, <:Integer}, multi_point_index -) where {N_in, N_out, ID <: LinearInterpolationDimension} - out = make_zero!!(out) - any(>(1), derivative_orders) && return out - - tᵢ = ntuple(i -> A.interp_dims[i].t[idx[i]], N_in) - tᵢ₊₁ = ntuple(i -> A.interp_dims[i].t[idx[i] + 1], N_in) - - # Size of the (hyper)rectangle `t` is in - t_vol = one(eltype(tᵢ)) - for (t₁, t₂) in zip(tᵢ, tᵢ₊₁) - t_vol *= t₂ - t₁ +) where {N_in, N_out, ID} + if isnothing(multi_point_index) + multi_point_index = ntuple(_ -> 1, N_in) end - - # Loop over the corners of the (hyper)rectangle `t` is in - for I in Iterators.product(ntuple(i -> (false, true), N_in)...) - c = eltype(out)(inv(t_vol)) - for (t_, right_point, d, t₁, t₂) in zip(t, I, derivative_orders, tᵢ, tᵢ₊₁) - c *= if right_point - iszero(d) ? t_ - t₁ : one(t_) - else - iszero(d) ? t₂ - t_ : -one(t_) - end - end - J = (ntuple(i -> idx[i] + I[i], N_in)..., ..) + out = make_zero!!(out) + # TODO: + # any(>(1), derivative_orders) && return out + # if any(>(0), derivative_orders) + # return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])), 1:N_in) + # typed_nan(out) + # else + # out + # end + # end + # Setup + space = map(iteration_space, A.interp_dims) + preparations = map(prepare, A.interp_dims, derivative_orders, multi_point_index, ts, idx) + # Loop over interpolation space + for I in Iterators.product(space...) + scaling = map(scale, A.interp_dims, preparations, I) + c = prod(scaling) + J = map(index, A.interp_dims, ts, idx, I) if iszero(N_out) out += c * A.u[J...] else @@ -38,63 +36,40 @@ function _interpolate!( return out end -function _interpolate!( - out, - A::NDInterpolation{N_in, N_out, ID}, - t::Tuple{Vararg{Number, N_in}}, - idx::NTuple{N_in, <:Integer}, - derivative_orders::NTuple{N_in, <:Integer}, - multi_point_index -) where {N_in, N_out, ID <: ConstantInterpolationDimension} - if any(>(0), derivative_orders) - return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])), 1:N_in) - typed_nan(out) - else - out - end - end - idx = ntuple( - i -> t[i] >= A.interp_dims[i].t[end] ? length(A.interp_dims[i].t) : idx[i], N_in) - if iszero(N_out) - out = A.u[idx...] - else - out .= A.u[idx...] - end - return out +function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) + t₁ = d.t[i] + t₂ = d.t[i + 1] + t_vol_inv = inv(t₂ - t₁) + return (; t, t₁, t₂, t_vol_inv, derivative_order) +end +prepare(::ConstantInterpolationDimension, derivative_orders, multi_point_index, t, i) = nothing +function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t::Number, i::Integer) + # TODO the dim_in arg isn't really needed, so drop it. Currently just 0 + basis_function_values = get_basis_function_values(d, t, i, derivative_order, multi_point_index, 0) + return (; basis_function_values) end -# BSpline evaluation -function _interpolate!( - out, - A::NDInterpolation{N_in, N_out, ID}, - t::Tuple{Vararg{Number, N_in}}, - idx::NTuple{N_in, <:Integer}, - derivative_orders::NTuple{N_in, <:Integer}, - multi_point_index -) where {N_in, N_out, ID <: BSplineInterpolationDimension} - (; interp_dims) = A - - out = make_zero!!(out) - degrees = ntuple(dim_in -> interp_dims[dim_in].degree, N_in) - basis_function_vals = get_basis_function_values_all( - A, t, idx, derivative_orders, multi_point_index - ) - - for I in CartesianIndices(ntuple(dim_in -> 1:(degrees[dim_in] + 1), N_in)) - B_product = prod(dim_in -> basis_function_vals[dim_in][I[dim_in]], 1:N_in) - cp_index = ntuple( - dim_in -> idx[dim_in] + I[dim_in] - degrees[dim_in] - 1, N_in) - if iszero(N_out) - out += B_product * A.u[cp_index...] - else - out .+= B_product * view(A.u, cp_index..., ..) - end - end +iteration_space(::LinearInterpolationDimension) = (false, true) +iteration_space(::ConstantInterpolationDimension) = 1 +iteration_space(d::BSplineInterpolationDimension) = 1:d.degree + 1 - return out +function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) + (; t, t₁, t₂, t_vol_inv, derivative_order) = prep + if right_point + iszero(derivative_order) ? t - t₁ : one(t) + else + iszero(derivative_order) ? t₂ - t : -one(t) + end * t_vol_inv end +scale(::ConstantInterpolationDimension, prep, i) = 1 +scale(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_function_values[i] + +index(::LinearInterpolationDimension, t, idx, i) = idx + i +index(d::ConstantInterpolationDimension, t, idx, i) = t >= d.t[end] ? length(d.t) : idx[i] +index(d::BSplineInterpolationDimension, t, idx, i) = idx + i - d.degree - 1 # NURBS evaluation +# TODO: generalise as above function _interpolate!( out, A::NDInterpolation{N_in, N_out, ID, <:NURBSWeights}, From b36c2d965910bf01deca2fb4f8ec4a33a89a1009 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 16 Oct 2025 17:15:35 +1100 Subject: [PATCH 02/32] cleanup --- src/interpolation_dimensions.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 1eb6b69..4cf7295 100644 --- a/src/interpolation_dimensions.jl +++ b/src/interpolation_dimensions.jl @@ -1,10 +1,3 @@ -""" - NoInterpolationDimension - -A dimension that does not perform interpolation. -""" -struct NoInterpolationDimension <: AbstractInterpolationDimension end - """ LinearInterpolationDimension(t; t_eval = similar(t, 0)) From ad6bdaa2761028d6c1519984379a26d8697f89ea Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 16 Oct 2025 22:31:03 +1100 Subject: [PATCH 03/32] simpler type, single _interpolat! method --- src/DataInterpolationsND.jl | 8 +-- src/interpolation_methods.jl | 97 ++++++++++++++---------------------- 2 files changed, 41 insertions(+), 64 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 6239c8c..65a25a6 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -22,13 +22,13 @@ the size of `u` along that dimension must match the length of `t` of the corresp - `u`: The array to be interpolated. """ struct NDInterpolation{ - N_in, N_out, - ID <: AbstractInterpolationDimension, + N_in, + N_out, gType <: AbstractInterpolationCache, uType <: AbstractArray } u::uType - interp_dims::NTuple{N_in, ID} + interp_dims::D cache::gType function NDInterpolation(u, interp_dims, cache) if interp_dims isa AbstractInterpolationDimension @@ -39,7 +39,7 @@ struct NDInterpolation{ @assert N_out≥0 "The number of dimensions of u must be at least the number of interpolation dimensions." validate_size_u(interp_dims, u) validate_cache(cache, interp_dims, u) - new{N_in, N_out, eltype(interp_dims), typeof(cache), typeof(u)}( + new{N_in, N_out, typeof(cache), typeof(u)}( u, interp_dims, cache ) end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 80f2552..9148929 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,41 +1,61 @@ function _interpolate!( out, - A::NDInterpolation{N_in, N_out, ID}, + A::NDInterpolation{N_in, N_out, <:NURBSWeights}, ts::Tuple{Vararg{Number, N_in}}, idx::NTuple{N_in, <:Integer}, derivative_orders::NTuple{N_in, <:Integer}, multi_point_index -) where {N_in, N_out, ID} +) where {N_in, N_out} + (; interp_dims, cache) = A + check_derivative_orders(interp_dims, derivative_orders) || return if isnothing(multi_point_index) multi_point_index = ntuple(_ -> 1, N_in) end out = make_zero!!(out) - # TODO: - # any(>(1), derivative_orders) && return out - # if any(>(0), derivative_orders) - # return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])), 1:N_in) - # typed_nan(out) - # else - # out - # end - # end + denom = zero(eltype(t)) # Setup - space = map(iteration_space, A.interp_dims) - preparations = map(prepare, A.interp_dims, derivative_orders, multi_point_index, ts, idx) - # Loop over interpolation space + space = map(iteration_space, interp_dims) + preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) + for I in Iterators.product(space...) scaling = map(scale, A.interp_dims, preparations, I) - c = prod(scaling) J = map(index, A.interp_dims, ts, idx, I) + product = if isnothing(cache) + scaling + else + weight = cache.weights[J...] + product = weight * scaling + denom += product + end if iszero(N_out) - out += c * A.u[J...] + out += product * A.u[J...] else - @. out += c * A.u[J...] + out .+= product * view(A.u, J..., ..) end end + + if !isnothing(cache) + if iszero(N_out) + out /= denom + else + out ./= denom + end + end + return out end +check_derivative_orders(dims, derivative_orders) = false +# TODO: +# any(>(1), derivative_orders) && return out +# if any(>(0), derivative_orders) +# return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])), 1:N_in) +# typed_nan(out) +# else +# out +# end +# end + function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) t₁ = d.t[i] t₂ = d.t[i + 1] @@ -67,46 +87,3 @@ scale(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_functio index(::LinearInterpolationDimension, t, idx, i) = idx + i index(d::ConstantInterpolationDimension, t, idx, i) = t >= d.t[end] ? length(d.t) : idx[i] index(d::BSplineInterpolationDimension, t, idx, i) = idx + i - d.degree - 1 - -# NURBS evaluation -# TODO: generalise as above -function _interpolate!( - out, - A::NDInterpolation{N_in, N_out, ID, <:NURBSWeights}, - t::Tuple{Vararg{Number, N_in}}, - idx::NTuple{N_in, <:Integer}, - derivative_orders::NTuple{N_in, <:Integer}, - multi_point_index -) where {N_in, N_out, ID <: BSplineInterpolationDimension} - (; interp_dims, cache) = A - - out = make_zero!!(out) - degrees = ntuple(dim_in -> interp_dims[dim_in].degree, N_in) - basis_function_vals = get_basis_function_values_all( - A, t, idx, derivative_orders, multi_point_index - ) - - denom = zero(eltype(t)) - - for I in CartesianIndices(ntuple(dim_in -> 1:(degrees[dim_in] + 1), N_in)) - B_product = prod(dim_in -> basis_function_vals[dim_in][I[dim_in]], 1:N_in) - cp_index = ntuple( - dim_in -> idx[dim_in] + I[dim_in] - degrees[dim_in] - 1, N_in) - weight = cache.weights[cp_index...] - product = weight * B_product - denom += product - if iszero(N_out) - out += product * A.u[cp_index...] - else - out .+= product * view(A.u, cp_index..., ..) - end - end - - if iszero(N_out) - out /= denom - else - out ./= denom - end - - return out -end From 3b6aa6481721f5c6771e42cb9c7b83eae428db29 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 16 Oct 2025 22:48:58 +1100 Subject: [PATCH 04/32] bugfix --- src/DataInterpolationsND.jl | 3 ++- src/interpolation_methods.jl | 27 +++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 65a25a6..5d35e46 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -25,6 +25,7 @@ struct NDInterpolation{ N_in, N_out, gType <: AbstractInterpolationCache, + D, uType <: AbstractArray } u::uType @@ -39,7 +40,7 @@ struct NDInterpolation{ @assert N_out≥0 "The number of dimensions of u must be at least the number of interpolation dimensions." validate_size_u(interp_dims, u) validate_cache(cache, interp_dims, u) - new{N_in, N_out, typeof(cache), typeof(u)}( + new{N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( u, interp_dims, cache ) end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 9148929..1cc1da1 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,40 +1,39 @@ function _interpolate!( out, - A::NDInterpolation{N_in, N_out, <:NURBSWeights}, + A::NDInterpolation{N_in, N_out}, ts::Tuple{Vararg{Number, N_in}}, idx::NTuple{N_in, <:Integer}, derivative_orders::NTuple{N_in, <:Integer}, multi_point_index ) where {N_in, N_out} - (; interp_dims, cache) = A - check_derivative_orders(interp_dims, derivative_orders) || return + (; interp_dims, cache, u) = A + check_derivative_orders(interp_dims, derivative_orders) || return out if isnothing(multi_point_index) multi_point_index = ntuple(_ -> 1, N_in) end out = make_zero!!(out) - denom = zero(eltype(t)) + denom = zero(eltype(ts)) # Setup space = map(iteration_space, interp_dims) preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) for I in Iterators.product(space...) - scaling = map(scale, A.interp_dims, preparations, I) - J = map(index, A.interp_dims, ts, idx, I) - product = if isnothing(cache) - scaling + scaling = map(scale, interp_dims, preparations, I) + J = map(index, interp_dims, ts, idx, I) + product = if cache isa EmptyCache + prod(scaling) else - weight = cache.weights[J...] - product = weight * scaling + product = cache.weights[J...] * prod(scaling) denom += product end if iszero(N_out) - out += product * A.u[J...] + out += product * u[J...] else - out .+= product * view(A.u, J..., ..) + out .+= product * view(u, J..., ..) end end - if !isnothing(cache) + if !(cache isa EmptyCache) if iszero(N_out) out /= denom else @@ -45,7 +44,7 @@ function _interpolate!( return out end -check_derivative_orders(dims, derivative_orders) = false +check_derivative_orders(dims, derivative_orders) = true # TODO: # any(>(1), derivative_orders) && return out # if any(>(0), derivative_orders) From 2e8b6356dbe35d3f459f0bb402b8ba678eeb0e50 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 16 Oct 2025 23:28:54 +1100 Subject: [PATCH 05/32] handle empty dimensions with NoInterpolationDimension --- src/DataInterpolationsND.jl | 16 +++++++------- src/interpolation_dimensions.jl | 7 +++++++ src/interpolation_methods.jl | 37 ++++++++++++++++++++------------- 3 files changed, 38 insertions(+), 22 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 5d35e46..d132536 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -31,21 +31,23 @@ struct NDInterpolation{ u::uType interp_dims::D cache::gType - function NDInterpolation(u, interp_dims, cache) - if interp_dims isa AbstractInterpolationDimension - interp_dims = (interp_dims,) - end - N_in = length(interp_dims) - N_out = ndims(u) - N_in + function NDInterpolation(u::AbstractArray{<:Any,N}, interp_dims, cache) where N + interp_dims = _add_trailing_interp_dims(interp_dims, Val{N}()) + N_out = count(map(d -> d isa NoInterpolationDimension, interp_dims)) @assert N_out≥0 "The number of dimensions of u must be at least the number of interpolation dimensions." validate_size_u(interp_dims, u) validate_cache(cache, interp_dims, u) - new{N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( + new{N, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( u, interp_dims, cache ) end end +_add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) = + _add_trailing_interp_dims((dim,), n) +_add_trailing_interp_dims(dims::Tuple, ::Val{N}) where N = + (dims..., ntuple(_ -> NoInterpolationDimension(), Val{N-length(dims)}())...) + # Constructor with optional global cache function NDInterpolation(u, interp_dims; cache = EmptyCache()) NDInterpolation(u, interp_dims, cache) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 4cf7295..1eb6b69 100644 --- 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)) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 1cc1da1..ceaba7c 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,15 +1,15 @@ function _interpolate!( out, - A::NDInterpolation{N_in, N_out}, - ts::Tuple{Vararg{Number, N_in}}, - idx::NTuple{N_in, <:Integer}, - derivative_orders::NTuple{N_in, <:Integer}, + A::NDInterpolation{N, N_out}, + ts::Tuple{Vararg{Number}}, + idx::NTuple{N, <:Integer}, + derivative_orders::NTuple{N, <:Integer}, multi_point_index -) where {N_in, N_out} +) where {N,N_out} (; interp_dims, cache, u) = A - check_derivative_orders(interp_dims, derivative_orders) || return out + check_derivative_order(interp_dims, derivative_orders) || return out if isnothing(multi_point_index) - multi_point_index = ntuple(_ -> 1, N_in) + multi_point_index = map(_ -> 1, interp_dims) end out = make_zero!!(out) denom = zero(eltype(ts)) @@ -29,7 +29,7 @@ function _interpolate!( if iszero(N_out) out += product * u[J...] else - out .+= product * view(u, J..., ..) + out .+= product .* view(u, J...) end end @@ -44,11 +44,14 @@ function _interpolate!( return out end -check_derivative_orders(dims, derivative_orders) = true -# TODO: -# any(>(1), derivative_orders) && return out -# if any(>(0), derivative_orders) -# return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])), 1:N_in) +check_derivative_order(dims::Tuple, derivative_orders::Tuple) = + all(map(check_derivative_order, dims, derivative_orders)) +check_derivative_order(::LinearInterpolationDimension, d_o) = d_o <= 1 +check_derivative_order(::ConstantInterpolationDimension, d_o) = d_0 <= 0 +check_derivative_order(::AbstractInterpolationDimension, d_o) = true +# TODO how to handle this +# if derivative_order > 0 +# return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])) # typed_nan(out) # else # out @@ -61,7 +64,8 @@ function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_ t_vol_inv = inv(t₂ - t₁) return (; t, t₁, t₂, t_vol_inv, derivative_order) end -prepare(::ConstantInterpolationDimension, derivative_orders, multi_point_index, t, i) = nothing +prepare(::ConstantInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) +prepare(::NoInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t::Number, i::Integer) # TODO the dim_in arg isn't really needed, so drop it. Currently just 0 basis_function_values = get_basis_function_values(d, t, i, derivative_order, multi_point_index, 0) @@ -70,6 +74,7 @@ end iteration_space(::LinearInterpolationDimension) = (false, true) iteration_space(::ConstantInterpolationDimension) = 1 +iteration_space(::NoInterpolationDimension) = 1 iteration_space(d::BSplineInterpolationDimension) = 1:d.degree + 1 function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) @@ -80,9 +85,11 @@ function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bo iszero(derivative_order) ? t₂ - t : -one(t) end * t_vol_inv end -scale(::ConstantInterpolationDimension, prep, i) = 1 +scale(::ConstantInterpolationDimension, prep::NamedTuple, i) = 1 +scale(::NoInterpolationDimension, prep::NamedTuple, i) = 1 scale(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_function_values[i] index(::LinearInterpolationDimension, t, idx, i) = idx + i index(d::ConstantInterpolationDimension, t, idx, i) = t >= d.t[end] ? length(d.t) : idx[i] +index(::NoInterpolationDimension, t, idx, i) = Colon() index(d::BSplineInterpolationDimension, t, idx, i) = idx + i - d.degree - 1 From c4595d268f6bffb0b51d5d590cdb3ff20c168de6 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 21 Oct 2025 15:01:12 +1100 Subject: [PATCH 06/32] more fixes for mixed dimensions --- src/DataInterpolationsND.jl | 20 +++++---- src/interpolation_methods.jl | 11 +++-- src/interpolation_parallel.jl | 51 ++++++++++++----------- src/interpolation_utils.jl | 78 +++++++++++++++-------------------- src/spline_utils.jl | 36 +++------------- test/test_interpolations.jl | 38 ++++++++++------- 6 files changed, 109 insertions(+), 125 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index d132536..2f01b33 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -22,6 +22,7 @@ the size of `u` along that dimension must match the length of `t` of the corresp - `u`: The array to be interpolated. """ struct NDInterpolation{ + N, N_in, N_out, gType <: AbstractInterpolationCache, @@ -33,16 +34,21 @@ struct NDInterpolation{ cache::gType function NDInterpolation(u::AbstractArray{<:Any,N}, interp_dims, cache) where N interp_dims = _add_trailing_interp_dims(interp_dims, Val{N}()) - N_out = count(map(d -> d isa NoInterpolationDimension, interp_dims)) + N_in = _count_interpolating_dims(interp_dims) + N_out = _count_noninterpolating_dims(interp_dims) @assert N_out≥0 "The number of dimensions of u must be at least the number of interpolation dimensions." validate_size_u(interp_dims, u) validate_cache(cache, interp_dims, u) - new{N, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( + new{N, N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( u, interp_dims, cache ) end end +# TODO probably not type-stable +_count_interpolating_dims(interp_dims) = count(map(d -> !(d isa NoInterpolationDimension), interp_dims)) +_count_noninterpolating_dims(interp_dims) = count(map(d -> d isa NoInterpolationDimension, interp_dims)) + _add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) = _add_trailing_interp_dims((dim,), n) _add_trailing_interp_dims(dims::Tuple, ::Val{N}) where N = @@ -73,11 +79,11 @@ 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} +function (interp::NDInterpolation{N,N_in,N_out})( + out::Union{Number, AbstractArray{<:Number, N_out}}, + t::Tuple{Vararg{Number, N}}; + derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) +) where {N,N_in,N_out} 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." diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index ceaba7c..24c862a 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,11 +1,11 @@ function _interpolate!( out, - A::NDInterpolation{N, N_out}, + A::NDInterpolation{N,N_in,N_out}, ts::Tuple{Vararg{Number}}, idx::NTuple{N, <:Integer}, derivative_orders::NTuple{N, <:Integer}, multi_point_index -) where {N,N_out} +) where {N,N_in,N_out} (; interp_dims, cache, u) = A check_derivative_order(interp_dims, derivative_orders) || return out if isnothing(multi_point_index) @@ -27,6 +27,7 @@ function _interpolate!( denom += product end if iszero(N_out) + @assert all(map(j -> j isa Integer, J)) out += product * u[J...] else out .+= product .* view(u, J...) @@ -66,9 +67,11 @@ function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_ end prepare(::ConstantInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) prepare(::NoInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) -function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t::Number, i::Integer) +function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) # TODO the dim_in arg isn't really needed, so drop it. Currently just 0 - basis_function_values = get_basis_function_values(d, t, i, derivative_order, multi_point_index, 0) + basis_function_values = get_basis_function_values( + d, t, i, derivative_order, multi_point_index + ) return (; basis_function_values) end diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index cce0a5d..85f79d4 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -63,10 +63,13 @@ out of place. - `derivative_orders`: The partial derivative order for each interpolation dimension. Defaults to `0` for each. """ -function eval_grid(interp::NDInterpolation{N_in}; kwargs...) where {N_in} - grid_size = map(itp_dim -> length(itp_dim.t_eval), interp.interp_dims) - out = similar(interp.u, (grid_size..., get_output_size(interp)...)) - eval_grid!(out, interp; kwargs...) +function eval_grid(interp::NDInterpolation; kwargs...) + sze = map(interp.interp_dims, size(interp.u)) do d, s + d isa NoInterpolationDimension ? s : length(d.t_eval) + end + # TODO: do we need to promote the type here, e.g. for eltype(u) <: Integer ? + out = similar(interp.u, sze) + return eval_grid!(out, interp; kwargs...) end """ @@ -87,10 +90,11 @@ function eval_grid!( interp::NDInterpolation{N_in}; derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) ) where {N_in} + used_interp_dims = _remove(NoInterpolationDimension, interp.interp_dims...) validate_derivative_orders(derivative_orders, interp; multi_point = true) backend = get_backend(out) - @assert all(i -> size(out, i) == length(interp.interp_dims[i].t_eval), N_in) "For the first N_in dimensions of out the length must match the t_eval of the corresponding interpolation dimension." - @assert size(out)[(N_in + 1):end]==get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." + @assert all(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in) "For the first N_in dimensions of out the length must match the t_eval of the corresponding interpolation dimension." + @assert size(out)[(N_in + 1):end] == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." eval_kernel(backend)( out, interp, @@ -104,29 +108,28 @@ end @kernel function eval_kernel( out, - @Const(A), + @Const(A::NDInterpolation{N, N_in, N_out}), derivative_orders, - eval_grid -) - N_in = length(A.interp_dims) - N_out = ndims(A.u) - N_in - + eval_grid, +) where {N, N_in, N_out} k = @index(Global, NTuple) + used_interp_dims = _remove(NoInterpolationDimension, A.interp_dims...) - if eval_grid - t_eval = ntuple(i -> A.interp_dims[i].t_eval[k[i]], N_in) - idx_eval = ntuple(i -> A.interp_dims[i].idx_eval[k[i]], N_in) - else - t_eval = ntuple(i -> A.interp_dims[i].t_eval[only(k)], N_in) - idx_eval = ntuple(i -> A.interp_dims[i].idx_eval[only(k)], N_in) - end + t_eval = ntuple(i -> used_interp_dims[i].t_eval[k[i]], N_in) + idx_eval = ntuple(i -> used_interp_dims[i].idx_eval[k[i]], N_in) + @show N_out if iszero(N_out) - out[k...] = _interpolate!( - make_out(A, t_eval), A, t_eval, idx_eval, derivative_orders, k) + dest = make_out(A, t_eval) + @show dest t_eval + out[k...] = _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) else - _interpolate!( - view(out, k..., ..), - A, t_eval, idx_eval, derivative_orders, k) + dest = view(out, k..., ..) + _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) end end + +# Remove objects of type T from splatted args (taken from DimensionalData.jl) +Base.@assume_effects :foldable _remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) +Base.@assume_effects :foldable _remove(::Type{T}, ::T, xs...) where T = _remove(T, xs...) +Base.@assume_effects :foldable _remove(::Type) = () diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 45db3aa..e9f8a60 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -3,18 +3,18 @@ trivial_range(i::Integer) = i:i Base.length(itp_dim::AbstractInterpolationDimension) = length(itp_dim.t) function validate_derivative_orders( - derivative_orders::NTuple{N_in, <:Integer}, - ::NDInterpolation{N_in}; + derivative_orders::NTuple{N, <:Integer}, + ::NDInterpolation{N}; kwargs... -) where {N_in} +) where {N} @assert all(≥(0), derivative_orders) "Derivative orders must me non-negative." end function validate_derivative_orders( - derivative_orders::NTuple{N_in, <:Integer}, - A::NDInterpolation{N_in, N_out, <:BSplineInterpolationDimension}; + derivative_orders::NTuple{N, <:Integer}, + A::NDInterpolation{N, N_in, N_out, <:BSplineInterpolationDimension}; multi_point::Bool = false -) where {N_in, N_out} +) where {N, N_in, N_out} @assert all(≥(0), derivative_orders) "Derivative orders must me non-negative." if multi_point @@ -42,7 +42,6 @@ function validate_size_u( ) where {N_in} @assert ntuple(i -> length(interp_dims[i]), N_in)==size(u)[1:N_in] "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." end - function validate_size_u( interp_dims::NTuple{N_in, <:BSplineInterpolationDimension}, u::AbstractArray @@ -51,55 +50,52 @@ function validate_size_u( @assert expected_size==size(u)[1:N_in] "Expected the size of the first N_in dimensions of u to be $expected_size based on the BSplineInterpolation properties." end -function validate_cache( - ::EmptyCache, ::NTuple{N_in, ID}, ::AbstractArray -) where {N_in, ID} - nothing +function validate_cache(cache, dims::NTuple, u) + ntuple(length(dims)) do n + validate_cache(cache, dims[n], u, n) + end end - +validate_cache(::EmptyCache, ::AbstractInterpolationDimension, ::AbstractArray, ::Int) = nothing function validate_cache( nurbs_weights::NURBSWeights, - ::NTuple{N_in, BSplineInterpolationDimension}, - u::AbstractArray -) where {N_in} - size_expected = size(u)[1:N_in] - @assert size(nurbs_weights.weights)==size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected)." + ::BSplineInterpolationDimension, + u::AbstractArray, + n::Int, +) + size_expected = size(u, n) + @assert size(nurbs_weights.weights, n) == size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected), got $(size(nurbs_weights.weights, n))." end - -function validate_cache( - ::gType, ::NTuple{N_in, ID}, ::AbstractArray) where {gType, N_in, ID} +function validate_cache(::gType, ::ID, ::AbstractArray, ::Int) where {gType,ID<:AbstractInterpolationDimension} @error("Interpolation dimension type $ID is not compatible with global cache type $gType.") end -function get_ts(interp_dims::NTuple{ - N_in, AbstractInterpolationDimension}) where {N_in} - ntuple(i -> interp_dims[i].t, N_in) -end - -function get_output_size(interp::NDInterpolation{N_in}) where {N_in} - size(interp.u)[(N_in + 1):end] +function get_output_size(interp::NDInterpolation) + # Replace non-interpolated dimensions with :, + # interpolated with their first index + I = map(interp.interp_dims, axes(interp.u)) do d, ax + d isa NoInterpolationDimension ? Colon() : first(ax) + end + return size(view(interp.u, I...)) end make_zero!!(::T) where {T <: Number} = zero(T) - function make_zero!!(v::T) where {T <: AbstractArray} v .= 0 v end function make_out( - interp::NDInterpolation{N_in, 0}, + interp::NDInterpolation{<:Any,N_in, 0}, t::NTuple{N_in, >:Number} ) where {N_in} zero(promote_type(eltype(interp.u), map(typeof, t)...)) end - function make_out( - interp::NDInterpolation{N_in}, + interp::NDInterpolation{<:Any,N_in}, t::NTuple{N_in, >:Number} ) where {N_in} - similar( - interp.u, promote_type(eltype(interp.u), map(eltype, t)...), get_output_size(interp)) + T = promote_type(eltype(interp.u), map(eltype, t)...) + similar(interp.u, T, get_output_size(interp)) end get_left(::AbstractInterpolationDimension) = false @@ -134,16 +130,12 @@ function get_idx( end end -function get_idx( - interp_dims::NTuple{N_in}, - t::Tuple{Vararg{Number, N_in}}; -) where {N_in} - ntuple(dim_in -> get_idx(interp_dims[dim_in], t[dim_in]), N_in) +function get_idx(interp_dims::NTuple{N_in}, t::Tuple{Vararg{Number, N_in}}) where N_in + used_interp_dims = _remove(NoInterpolationDimension, interp_dims...) + map(get_idx, used_interp_dims, t) end -function set_eval_idx!( - interp_dim::AbstractInterpolationDimension, -) +function set_eval_idx!(interp_dim::AbstractInterpolationDimension) backend = get_backend(interp_dim.t) if !isempty(interp_dim.t_eval) set_idx_kernel(backend)( @@ -154,9 +146,7 @@ function set_eval_idx!( synchronize(backend) end -@kernel function set_idx_kernel( - interp_dim -) +@kernel function set_idx_kernel(interp_dim) i = @index(Global, Linear) interp_dim.idx_eval[i] = get_idx(interp_dim, interp_dim.t_eval[i]) end diff --git a/src/spline_utils.jl b/src/spline_utils.jl index 90d9e8d..0621563 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -61,7 +61,6 @@ function get_basis_function_values( idx::Integer, derivative_order::Integer, multi_point_index::Nothing, - dim_in::Integer ) (; degree, knots_all) = itp_dim T = promote_type(typeof(t), eltype(itp_dim.basis_function_eval)) @@ -92,7 +91,6 @@ function get_basis_function_values( basis_function_values[1:degree_plus_1] end - # Get the basis function values for one point in an # unstructured multi point evaluation (given by the scalar multi point index) function get_basis_function_values( @@ -101,12 +99,9 @@ function get_basis_function_values( idx::Integer, derivative_order::Integer, multi_point_index::Number, - dim_in::Integer ) - view(itp_dim.basis_function_eval, - multi_point_index, :, derivative_order + 1) + view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) end - # Get the basis function values for one point in a # grid evaluation (given by the tuple multi point index) function get_basis_function_values( @@ -114,27 +109,9 @@ function get_basis_function_values( t::Number, idx::Integer, derivative_order::Integer, - multi_point_index::NTuple{N_in, <:Integer}, - dim_in::Integer -) where {N_in} - view(itp_dim.basis_function_eval, - multi_point_index[dim_in], :, derivative_order + 1) -end - -# Get all basis function values to evaluate a BSpline interpolation in t -function get_basis_function_values_all( - A::NDInterpolation{N_in, N_out, <:BSplineInterpolationDimension}, - t::Tuple{Vararg{Number, N_in}}, - idx::NTuple{N_in, <:Integer}, - derivative_orders::NTuple{N_in, <:Integer}, - multi_point_index -) where {N_in, N_out} - ntuple( - dim_in -> get_basis_function_values( - A.interp_dims[dim_in], t[dim_in], idx[dim_in], derivative_orders[dim_in], multi_point_index, dim_in - ), - N_in - ) + multi_point_index::Integer, +) + view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) end function set_basis_function_eval!(itp_dim::BSplineInterpolationDimension)::Nothing @@ -152,15 +129,12 @@ end ) i, derivative_order_plus_1 = @index(Global, NTuple) - itp_dim.basis_function_eval[i, - :, - derivative_order_plus_1] .= get_basis_function_values( + itp_dim.basis_function_eval[i, :, derivative_order_plus_1] .= get_basis_function_values( itp_dim, itp_dim.t_eval[i], itp_dim.idx_eval[i], derivative_order_plus_1 - 1, nothing, - 0 ) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index cbbe70f..10a1bff 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -1,10 +1,14 @@ -using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache using DataInterpolationsND using Random +using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache + function test_globally_constant( - ID::Type{<:AbstractInterpolationDimension}; args1 = (), args2 = (), kwargs1 = (), - kwargs2 = (), cache = EmptyCache(), test_derivatives = true) + ID::Type{<:AbstractInterpolationDimension}; args1 = (), args2 = (), kwargs1 = (), + kwargs2 = (), + cache = EmptyCache(), + test_derivatives = true +) t1 = [-3.14, 1.0, 3.0, 7.6, 12.8] t2 = [-2.71, 1.41, 12.76, 50.2, 120.0] @@ -77,11 +81,14 @@ end end @testset "BSpline Interpolation" begin - test_globally_constant( - BSplineInterpolationDimension, args1 = (2,), args2 = (3,), - kwargs1 = (:max_derivative_order_eval => 1,), - kwargs2 = (:max_derivative_order_eval => 1,) - ) + ID = BSplineInterpolationDimension + args1 = (2,) + args2 = (3,) + kwargs1 = (:max_derivative_order_eval => 1,) + kwargs2 = (:max_derivative_order_eval => 1,) + cache = EmptyCache() + test_derivatives = true + test_globally_constant(ID; args1, args2, kwargs1, kwargs2, cache, test_derivatives) f(t1, t2, t3) = t1^2 + t2^2 + t3^2 @@ -99,13 +106,14 @@ end @testset "NURBS Interpolation" begin Random.seed!(10) - test_globally_constant( - BSplineInterpolationDimension; args1 = (3,), args2 = (1,), - kwargs1 = (:max_derivative_order_eval => 1,), - kwargs2 = (:max_derivative_order_eval => 1,), - cache = NURBSWeights(rand(7, 5)), - test_derivatives = false - ) + ID = BSplineInterpolationDimension + args1 = (3,) + args2 = (1,) + kwargs1 = (:max_derivative_order_eval => 1,) + kwargs2 = (:max_derivative_order_eval => 1,) + cache = NURBSWeights(rand(7, 5)) + test_derivatives = false + test_globally_constant(ID; args1, args2, kwargs1, kwargs2, cache, test_derivatives) ## Circle representation # Knots From b384f41bbf9301c51d64e5f4255368d400a764d0 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 21 Oct 2025 16:04:22 +1100 Subject: [PATCH 07/32] more mixed dims --- src/DataInterpolationsND.jl | 4 +- src/interpolation_dimensions.jl | 15 +++---- src/interpolation_parallel.jl | 2 +- src/interpolation_utils.jl | 77 +++++++++++++++++---------------- src/spline_utils.jl | 1 + test/test_interpolations.jl | 7 +-- 6 files changed, 52 insertions(+), 54 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 2f01b33..beaf48b 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -84,7 +84,7 @@ function (interp::NDInterpolation{N,N_in,N_out})( t::Tuple{Vararg{Number, N}}; derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) ) where {N,N_in,N_out} - validate_derivative_orders(derivative_orders, interp) + validate_derivative_order(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." _interpolate!(out, interp, t, idx, derivative_orders, nothing) @@ -97,7 +97,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! end # module DataInterpolationsND diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 1eb6b69..adf05a8 100644 --- a/src/interpolation_dimensions.jl +++ b/src/interpolation_dimensions.jl @@ -1,5 +1,5 @@ """ - NoInterpolationDimension + NoInterpolationDimensio A dimension that does not perform interpolation. """ @@ -164,15 +164,10 @@ 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 - ) - ) + s = (length(t_eval), degree + 1, max_derivative_order_eval + 1) + @show s + T = typeof(inv(one(eltype(t))) * inv(one(eltype(t_eval)))) + 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) diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 85f79d4..14c48cf 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -91,7 +91,7 @@ function eval_grid!( derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) ) where {N_in} used_interp_dims = _remove(NoInterpolationDimension, interp.interp_dims...) - validate_derivative_orders(derivative_orders, interp; multi_point = true) + validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) @assert all(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in) "For the first N_in dimensions of out the length must match the t_eval of the corresponding interpolation dimension." @assert size(out)[(N_in + 1):end] == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index e9f8a60..d86b724 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -2,52 +2,56 @@ trivial_range(i::Integer) = i:i Base.length(itp_dim::AbstractInterpolationDimension) = length(itp_dim.t) -function validate_derivative_orders( - derivative_orders::NTuple{N, <:Integer}, - ::NDInterpolation{N}; - kwargs... -) where {N} - @assert all(≥(0), derivative_orders) "Derivative orders must me non-negative." +function validate_derivative_order(derivative_orders::NTuple, A::NDInterpolation; + multi_point::Bool = false +) + map(derivative_orders, A.interp_dims) do d_o, d + validate_derivative_order(d_o, d; multi_point, cache=A.cache) + end end - -function validate_derivative_orders( - derivative_orders::NTuple{N, <:Integer}, - A::NDInterpolation{N, N_in, N_out, <:BSplineInterpolationDimension}; - multi_point::Bool = false -) where {N, N_in, N_out} - @assert all(≥(0), derivative_orders) "Derivative orders must me non-negative." - +function validate_derivative_order( + derivative_order::Integer, + interp_dim::BSplineInterpolationDimension; + multi_point::Bool, + cache, +) if multi_point - @assert all( - i -> derivative_orders[i] ≤ A.interp_dims[i].max_derivative_order_eval, 1:N_in - ) "For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be \ - larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want \ - to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the \ - `BSplineInterpolationDimension` constructor(s)." - end - - if A.cache isa NURBSWeights - @assert all(==(0), derivative_orders) "Currently partial derivatives of NURBS are not supported." + @assert derivative_order ≤ interp_dim.max_derivative_order_eval """ + For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be + larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want + to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the + `BSplineInterpolationDimension` constructor(s). + """ end + validate_derivative_order_by_cache(cache, derivative_order) end +function validate_derivative_order( + derivative_order::Integer, + interp_dim::AbstractInterpolationDimension; + multi_point::Bool, + cache, +) + validate_derivative_order_by_cache(cache, derivative_order) +end + +validate_derivative_order_by_cache(::NURBSWeights, derivative_order) = + @assert derivative_order == 0 "Currently partial derivatives of NURBS are not supported." +validate_derivative_order_by_cache(::Any, derivative_order) = + @assert derivative_order >= 0 "Derivative orders must me non-negative." + function validate_t(t) @assert t isa AbstractVector{<:Number} "t must be an AbstractVector with number like elements." @assert all(>(0), diff(t)) "The elements of t must be sorted and unique." end -function validate_size_u( - interp_dims::NTuple{N_in, <:AbstractInterpolationDimension}, - u::AbstractArray -) where {N_in} - @assert ntuple(i -> length(interp_dims[i]), N_in)==size(u)[1:N_in] "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." +validate_size_u(interp_dims::NTuple, u) = map(validate_size_u, interp_dims, axes(u)) +function validate_size_u(interp_dim::AbstractInterpolationDimension, ax) + @assert length(interp_dim) == length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." end -function validate_size_u( - interp_dims::NTuple{N_in, <:BSplineInterpolationDimension}, - u::AbstractArray -) where {N_in} - expected_size = ntuple(dim_in -> get_n_basis_functions(interp_dims[dim_in]), N_in) - @assert expected_size==size(u)[1:N_in] "Expected the size of the first N_in dimensions of u to be $expected_size based on the BSplineInterpolation properties." +function validate_size_u(interp_dim::BSplineInterpolationDimension, ax) + expected_size = get_n_basis_functions(interp_dim) + @assert expected_size == length(ax) "Expected the size of the first N_in dimensions of u to be $expected_size based on the BSplineInterpolation properties." end function validate_cache(cache, dims::NTuple, u) @@ -70,8 +74,6 @@ function validate_cache(::gType, ::ID, ::AbstractArray, ::Int) where {gType,ID<: end function get_output_size(interp::NDInterpolation) - # Replace non-interpolated dimensions with :, - # interpolated with their first index I = map(interp.interp_dims, axes(interp.u)) do d, ax d isa NoInterpolationDimension ? Colon() : first(ax) end @@ -129,7 +131,6 @@ function get_idx( clamp(searchsortedlast(t, t_eval) + idx_shift, lb, ub) end end - function get_idx(interp_dims::NTuple{N_in}, t::Tuple{Vararg{Number, N_in}}) where N_in used_interp_dims = _remove(NoInterpolationDimension, interp_dims...) map(get_idx, used_interp_dims, t) diff --git a/src/spline_utils.jl b/src/spline_utils.jl index 0621563..81fc3d6 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -111,6 +111,7 @@ function get_basis_function_values( derivative_order::Integer, multi_point_index::Integer, ) + @show size(itp_dim.basis_function_eval) view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 10a1bff..7e668a9 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -48,15 +48,17 @@ function test_globally_constant( end end -function test_analytic(itp::NDInterpolation{N_in}, f) where {N_in} +function test_analytic(itp::NDInterpolation{<:Any,N_in}, f) where {N_in} + used_interp_dims = DataInterpolationsND._remove(NoInterpolationDimension, itp.interp_dims...) # Evaluation in data points - ts = ntuple(dim_in -> itp.interp_dims[dim_in].t, N_in) + ts = map(d -> d.t, used_interp_dims) for t in Iterators.product(ts...) @test itp(t) ≈ f(t...) end # Evaluation between data points ts_ = ntuple(dim_in -> ts[dim_in][1:(end - 1)] + diff(ts[dim_in]) / 2, N_in) + t = first(Iterators.product(ts_...)) for t in Iterators.product(ts_...) @test itp(t) ≈ f(t...) end @@ -66,7 +68,6 @@ end test_globally_constant(LinearInterpolationDimension) f(t1, t2) = 3.0 + 2.3t1 - 4.7t2 - Random.seed!(1) t1 = cumsum(rand(10)) t2 = cumsum(rand(10)) From d6d652dd07d77652e69b17418a0cf9580b5c2d3b Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 21 Oct 2025 16:20:58 +1100 Subject: [PATCH 08/32] cleanup --- src/spline_utils.jl | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/spline_utils.jl b/src/spline_utils.jl index 81fc3d6..213823b 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -93,17 +93,6 @@ function get_basis_function_values( end # Get the basis function values for one point in an # unstructured multi point evaluation (given by the scalar multi point index) -function get_basis_function_values( - itp_dim::BSplineInterpolationDimension, - t::Number, - idx::Integer, - derivative_order::Integer, - multi_point_index::Number, -) - view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) -end -# Get the basis function values for one point in a -# grid evaluation (given by the tuple multi point index) function get_basis_function_values( itp_dim::BSplineInterpolationDimension, t::Number, @@ -111,8 +100,8 @@ function get_basis_function_values( derivative_order::Integer, multi_point_index::Integer, ) - @show size(itp_dim.basis_function_eval) - view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) + view(itp_dim.basis_function_eval, + multi_point_index, :, derivative_order + 1) end function set_basis_function_eval!(itp_dim::BSplineInterpolationDimension)::Nothing From 87dd1abcad5c37f9a3b73160ee5d5abb924ab6ad Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Wed, 22 Oct 2025 11:04:44 +1100 Subject: [PATCH 09/32] tweaks --- src/DataInterpolationsND.jl | 9 ++---- src/interpolation_methods.jl | 13 +++++---- src/interpolation_parallel.jl | 44 ++++++++++++++++------------ src/interpolation_utils.jl | 55 ++++++++++++++++++++++++++++------- src/spline_utils.jl | 2 +- test/test_interpolations.jl | 23 +++++++++++++++ 6 files changed, 105 insertions(+), 41 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index beaf48b..bba8835 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -45,7 +45,7 @@ struct NDInterpolation{ end end -# TODO probably not type-stable +# TODO probably not type-stable (this needs to compile away completely) _count_interpolating_dims(interp_dims) = count(map(d -> !(d isa NoInterpolationDimension), interp_dims)) _count_noninterpolating_dims(interp_dims) = count(map(d -> d isa NoInterpolationDimension, interp_dims)) @@ -72,24 +72,21 @@ include("plot_rec.jl") function (interp::NDInterpolation)(t_args::Vararg{Number}; kwargs...) interp(t_args; kwargs...) end - function (interp::NDInterpolation)( out::AbstractArray, t_args::Vararg{Number}; kwargs...) interp(out, t_args; kwargs...) end - # In place single input evaluation function (interp::NDInterpolation{N,N_in,N_out})( out::Union{Number, AbstractArray{<:Number, N_out}}, t::Tuple{Vararg{Number, N}}; derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) ) where {N,N_in,N_out} + validate_size_u(interp, out) validate_derivative_order(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." - _interpolate!(out, interp, t, idx, derivative_orders, nothing) + return _interpolate!(out, interp, t, idx, derivative_orders, nothing) end - # Out of place single input evaluation function (interp::NDInterpolation)(t::Tuple{Vararg{Number}}; kwargs...) out = make_out(interp, t) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 24c862a..cb957aa 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,20 +1,22 @@ function _interpolate!( - out, + out::Union{Number,AbstractArray{<:Any,N_out}}, A::NDInterpolation{N,N_in,N_out}, - ts::Tuple{Vararg{Number}}, - idx::NTuple{N, <:Integer}, - derivative_orders::NTuple{N, <:Integer}, + ts::Tuple, + idx::Tuple, + derivative_orders::Tuple, multi_point_index ) where {N,N_in,N_out} (; interp_dims, cache, u) = A + check_derivative_order(interp_dims, derivative_orders) || return out if isnothing(multi_point_index) multi_point_index = map(_ -> 1, interp_dims) end out = make_zero!!(out) - denom = zero(eltype(ts)) + denom = zero(eltype(_remove(Nothing, ts...))) # Setup space = map(iteration_space, interp_dims) + @show prepare interp_dims derivative_orders multi_point_index ts idx preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) for I in Iterators.product(space...) @@ -23,6 +25,7 @@ function _interpolate!( product = if cache isa EmptyCache prod(scaling) else + @show J size(cache.weights) product = cache.weights[J...] * prod(scaling) denom += product end diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 14c48cf..f16bd6c 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -87,20 +87,20 @@ in place. """ function eval_grid!( out::AbstractArray, - interp::NDInterpolation{N_in}; - derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) -) where {N_in} + interp::NDInterpolation{N,N_in}; + derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) +) where {N,N_in} used_interp_dims = _remove(NoInterpolationDimension, interp.interp_dims...) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) - @assert all(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in) "For the first N_in dimensions of out the length must match the t_eval of the corresponding interpolation dimension." + @assert all(ntuple(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in)) "The length must match the t_eval of the corresponding interpolation dimension." @assert size(out)[(N_in + 1):end] == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." eval_kernel(backend)( out, interp, derivative_orders, true, - ndrange = size(out)[1:N_in] + ndrange = get_ndrange(interp) ) synchronize(backend) return out @@ -108,28 +108,36 @@ end @kernel function eval_kernel( out, - @Const(A::NDInterpolation{N, N_in, N_out}), + A::NDInterpolation{N, N_in, N_out}, derivative_orders, eval_grid, ) where {N, N_in, N_out} - k = @index(Global, NTuple) - used_interp_dims = _remove(NoInterpolationDimension, A.interp_dims...) + # This kernel is only over interpolated dimensions, we need + # to insert fillers to match the number of dimensions in the data + I = @index(Global, NTuple) + k = insert_colon(A, I) - t_eval = ntuple(i -> used_interp_dims[i].t_eval[k[i]], N_in) - idx_eval = ntuple(i -> used_interp_dims[i].idx_eval[k[i]], N_in) + t_eval = map(get_t_eval, A.interp_dims, k) + idx_eval = map(get_idx_eval, A.interp_dims, k) - @show N_out if iszero(N_out) dest = make_out(A, t_eval) - @show dest t_eval - out[k...] = _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) + out[I...] = _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) else - dest = view(out, k..., ..) + dest = view(out, k...) _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) end end -# Remove objects of type T from splatted args (taken from DimensionalData.jl) -Base.@assume_effects :foldable _remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) -Base.@assume_effects :foldable _remove(::Type{T}, ::T, xs...) where T = _remove(T, xs...) -Base.@assume_effects :foldable _remove(::Type) = () +get_t_eval(d, i) = d.t_eval[i] +get_t_eval(d::NoInterpolationDimension, i) = nothing + +get_idx_eval(d, i) = d.idx_eval[i] +get_idx_eval(d::NoInterpolationDimension, i) = nothing + +# Insert `nothing` into `I` where there is a NoInterpolationDimension +# This is to expand the index from KernelAbstractions kernel back to +# the original dimensionality, where `nothing` is just padding that isn't used. +insert_colon(A::NDInterpolation, I) = insert_colon(A.interp_dims, I) +insert_colon(interp_dims::Tuple, I) = + _insertat(NoInterpolationDimension, Colon(), I, (), interp_dims...) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index d86b724..74da9dc 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -45,21 +45,26 @@ function validate_t(t) @assert all(>(0), diff(t)) "The elements of t must be sorted and unique." end -validate_size_u(interp_dims::NTuple, u) = map(validate_size_u, interp_dims, axes(u)) -function validate_size_u(interp_dim::AbstractInterpolationDimension, ax) +validate_size_u(::NDInterpolation{N,N,0}, u::Number) where N = nothing +validate_size_u(interp::NDInterpolation, u::AbstractArray) = validate_size_u(interp.interp_dims, axes(u)) +validate_size_u(interp_dims::Tuple, u::AbstractArray) = map(validate_size_u, interp_dims, axes(u)) +validate_size_u(interp_dim::NoInterpolationDimension, ax::AbstractRange) = nothing +function validate_size_u(interp_dim::AbstractInterpolationDimension, ax::AbstractRange) @assert length(interp_dim) == length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." end -function validate_size_u(interp_dim::BSplineInterpolationDimension, ax) +function validate_size_u(interp_dim::BSplineInterpolationDimension, ax::AbstractRange) expected_size = get_n_basis_functions(interp_dim) - @assert expected_size == length(ax) "Expected the size of the first N_in dimensions of u to be $expected_size based on the BSplineInterpolation properties." + @assert expected_size == length(ax) "Expected the size to be $expected_size based on the BSplineInterpolation properties, got $(length(ax))." end -function validate_cache(cache, dims::NTuple, u) +function validate_cache(cache::AbstractInterpolationCache, dims::Tuple, u::AbstractArray) ntuple(length(dims)) do n validate_cache(cache, dims[n], u, n) end end validate_cache(::EmptyCache, ::AbstractInterpolationDimension, ::AbstractArray, ::Int) = nothing +validate_cache(::EmptyCache, ::NoInterpolationDimension, ::AbstractArray, ::Int) = nothing +validate_cache(::AbstractInterpolationCache, ::NoInterpolationDimension, ::AbstractArray, ::Int) = nothing function validate_cache( nurbs_weights::NURBSWeights, ::BSplineInterpolationDimension, @@ -75,9 +80,9 @@ end function get_output_size(interp::NDInterpolation) I = map(interp.interp_dims, axes(interp.u)) do d, ax - d isa NoInterpolationDimension ? Colon() : first(ax) + d isa NoInterpolationDimension ? length(ax) : nothing end - return size(view(interp.u, I...)) + return _remove(Nothing, I...) end make_zero!!(::T) where {T <: Number} = zero(T) @@ -131,9 +136,9 @@ function get_idx( clamp(searchsortedlast(t, t_eval) + idx_shift, lb, ub) end end -function get_idx(interp_dims::NTuple{N_in}, t::Tuple{Vararg{Number, N_in}}) where N_in - used_interp_dims = _remove(NoInterpolationDimension, interp_dims...) - map(get_idx, used_interp_dims, t) +get_idx(::NoInterpolationDimension, t_eval::Number) = nothing +function get_idx(interp_dims::NTuple{N}, t::Tuple{Vararg{Number, N}}) where N + map(get_idx, interp_dims, t) end function set_eval_idx!(interp_dim::AbstractInterpolationDimension) @@ -147,7 +152,7 @@ function set_eval_idx!(interp_dim::AbstractInterpolationDimension) synchronize(backend) end -@kernel function set_idx_kernel(interp_dim) +@kernel function set_idx_kernel(interp_dim::AbstractInterpolationDimension) i = @index(Global, Linear) interp_dim.idx_eval[i] = get_idx(interp_dim, interp_dim.t_eval[i]) end @@ -162,3 +167,31 @@ end typed_nan(::T) where {T <: Integer} = zero(T) typed_nan(::T) where {T <: AbstractFloat} = T(NaN) + + +# Get the KernelAbstractions nd_range, over interpolated dimensions +function get_ndrange(interp::NDInterpolation) + I = map(interp.interp_dims) do d + d isa NoInterpolationDimension ? nothing : length(d.t_eval) + end + return _remove(Nothing, I...) +end + +# Remove objects of type T from splatted args (taken from DimensionalData.jl) +Base.@assume_effects :foldable _remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) +Base.@assume_effects :foldable _remove(::Type{T}, ::T, xs...) where T = _remove(T, xs...) +Base.@assume_effects :foldable _remove(::Type) = () + +# Insert x in `out` where `m isa T`, otherwise output one of in for each m +# out must exactly match the number of ms or this will error +# If `!(m isa T)` take from `in` +Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, out::Tuple, m, ms...) where T = + _insertat(T, x, Base.tail(in), (out..., first(in)), ms...) +# If `m isa T` insert `x` +Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, out::Tuple, m::T, ms...) where T = + _insertat(T, x, in, (out..., x), ms...) +# `in` can be empty if all trailing ms are T +Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple, out::Tuple, m::T, ms::T...) where T = + _insertat(T, x, in, (out..., x), ms...) +# `in` can also be empty if there are no remaining `m` +Base.@assume_effects :foldable _insertat(::Type, x, in::Tuple{}, out::Tuple) = out diff --git a/src/spline_utils.jl b/src/spline_utils.jl index 213823b..cff7f5c 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -60,7 +60,7 @@ function get_basis_function_values( t::Number, idx::Integer, derivative_order::Integer, - multi_point_index::Nothing, + multi_point_index::Union{Nothing,Colon}, # TODO why both ) (; degree, knots_all) = itp_dim T = promote_type(typeof(t), eltype(itp_dim.basis_function_eval)) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 7e668a9..ecd2155 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -79,6 +79,13 @@ end u = f.(t1, t2') itp = NDInterpolation(u, itp_dims) test_analytic(itp, f) + + itp_dims = ( + NoInterpolationDimension(), + LinearInterpolationDimension(t2; t_eval = t2 ./ 2) + ) + itp = NDInterpolation(u, itp_dims) + eval_grid(itp) end @testset "BSpline Interpolation" begin @@ -149,3 +156,19 @@ end @test allunique(points_on_circle[2:end]) @test all(point -> point[1]^2 + point[2]^2 ≈ 1, points_on_circle) end + +@testset "Mixed Interpolation" begin + + t1 = cumsum(rand(3)) + t2 = cumsum(rand(4)) + t3 = collect(0:(π / 2):(2π)) + t_eval = collect(range(0, 2π, length = 100)) + u = zeros(3, 4, 5) + itp_dims = ( + NoInterpolationDimension(), + LinearInterpolationDimension(t2; t_eval = t2 ./ 2), + BSplineInterpolationDimension(t3, 2; multiplicities, t_eval), + ) + itp = NDInterpolation(u, itp_dims) + eval_grid(itp) +end From 107dd45c53643745397000f5c0764855dcee4a2a Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 23 Oct 2025 12:47:25 +1100 Subject: [PATCH 10/32] working besides Constant nans --- src/interpolation_dimensions.jl | 1 - src/interpolation_methods.jl | 26 +++++----- src/interpolation_parallel.jl | 20 ++------ src/interpolation_utils.jl | 91 +++++++++++++++++++++++++++------ src/spline_utils.jl | 11 ++++ test/test_interpolations.jl | 41 +++++++++++---- 6 files changed, 132 insertions(+), 58 deletions(-) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index adf05a8..8c8d12b 100644 --- a/src/interpolation_dimensions.jl +++ b/src/interpolation_dimensions.jl @@ -165,7 +165,6 @@ function BSplineInterpolationDimension( idx_eval = similar(t_eval, Int) s = (length(t_eval), degree + 1, max_derivative_order_eval + 1) - @show s T = typeof(inv(one(eltype(t))) * inv(one(eltype(t_eval)))) basis_function_eval = similar(t_eval, T, s) itp_dim = BSplineInterpolationDimension( diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index cb957aa..76949b1 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,36 +1,34 @@ function _interpolate!( - out::Union{Number,AbstractArray{<:Any,N_out}}, - A::NDInterpolation{N,N_in,N_out}, - ts::Tuple, - idx::Tuple, - derivative_orders::Tuple, + out::Union{Number, AbstractArray{<:Any, N_out}}, + A::NDInterpolation{N, N_in, N_out}, + ts::Tuple{Vararg{Any, N}}, + idx::Tuple{Vararg{Any ,N}}, + derivative_orders::Tuple{Vararg{Any, N}}, multi_point_index ) where {N,N_in,N_out} (; interp_dims, cache, u) = A check_derivative_order(interp_dims, derivative_orders) || return out if isnothing(multi_point_index) - multi_point_index = map(_ -> 1, interp_dims) + multi_point_index = map(_ -> nothing, interp_dims) end out = make_zero!!(out) - denom = zero(eltype(_remove(Nothing, ts...))) + denom = zero(eltype(u)) # Setup space = map(iteration_space, interp_dims) - @show prepare interp_dims derivative_orders multi_point_index ts idx preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) for I in Iterators.product(space...) scaling = map(scale, interp_dims, preparations, I) J = map(index, interp_dims, ts, idx, I) - product = if cache isa EmptyCache - prod(scaling) + if cache isa EmptyCache + product = prod(scaling) else - @show J size(cache.weights) - product = cache.weights[J...] * prod(scaling) + K = removeat(NoInterpolationDimension, J, interp_dims) + product = cache.weights[K...] * prod(scaling) denom += product end if iszero(N_out) - @assert all(map(j -> j isa Integer, J)) out += product * u[J...] else out .+= product .* view(u, J...) @@ -51,7 +49,7 @@ end check_derivative_order(dims::Tuple, derivative_orders::Tuple) = all(map(check_derivative_order, dims, derivative_orders)) check_derivative_order(::LinearInterpolationDimension, d_o) = d_o <= 1 -check_derivative_order(::ConstantInterpolationDimension, d_o) = d_0 <= 0 +check_derivative_order(::ConstantInterpolationDimension, d_o) = d_o <= 0 check_derivative_order(::AbstractInterpolationDimension, d_o) = true # TODO how to handle this # if derivative_order > 0 diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index f16bd6c..98ae8f8 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -64,11 +64,8 @@ out of place. - `derivative_orders`: The partial derivative order for each interpolation dimension. Defaults to `0` for each. """ function eval_grid(interp::NDInterpolation; kwargs...) - sze = map(interp.interp_dims, size(interp.u)) do d, s - d isa NoInterpolationDimension ? s : length(d.t_eval) - end # TODO: do we need to promote the type here, e.g. for eltype(u) <: Integer ? - out = similar(interp.u, sze) + out = similar(interp.u, grid_size(interp)) return eval_grid!(out, interp; kwargs...) end @@ -90,11 +87,11 @@ function eval_grid!( interp::NDInterpolation{N,N_in}; derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) ) where {N,N_in} - used_interp_dims = _remove(NoInterpolationDimension, interp.interp_dims...) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) - @assert all(ntuple(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in)) "The length must match the t_eval of the corresponding interpolation dimension." - @assert size(out)[(N_in + 1):end] == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." + # used_interp_dims = remove(NoInterpolationDimension, interp.interp_dims) + # @assert all(ntuple(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in)) "The length must match the t_eval of the corresponding interpolation dimension." + # @assert size(out) == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." eval_kernel(backend)( out, interp, @@ -115,7 +112,7 @@ end # This kernel is only over interpolated dimensions, we need # to insert fillers to match the number of dimensions in the data I = @index(Global, NTuple) - k = insert_colon(A, I) + k = insertat(NoInterpolationDimension, Colon(), I, A.interp_dims) t_eval = map(get_t_eval, A.interp_dims, k) idx_eval = map(get_idx_eval, A.interp_dims, k) @@ -134,10 +131,3 @@ get_t_eval(d::NoInterpolationDimension, i) = nothing get_idx_eval(d, i) = d.idx_eval[i] get_idx_eval(d::NoInterpolationDimension, i) = nothing - -# Insert `nothing` into `I` where there is a NoInterpolationDimension -# This is to expand the index from KernelAbstractions kernel back to -# the original dimensionality, where `nothing` is just padding that isn't used. -insert_colon(A::NDInterpolation, I) = insert_colon(A.interp_dims, I) -insert_colon(interp_dims::Tuple, I) = - _insertat(NoInterpolationDimension, Colon(), I, (), interp_dims...) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 74da9dc..969b4c6 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -82,7 +82,18 @@ function get_output_size(interp::NDInterpolation) I = map(interp.interp_dims, axes(interp.u)) do d, ax d isa NoInterpolationDimension ? length(ax) : nothing end - return _remove(Nothing, I...) + return remove(Nothing, I) +end + +function grid_size(interp::NDInterpolation) + (; interp_dims) = interp + # TODO: put this in a function, but + # Get the size of dims that are not NoInterpolationDimension + interp_size = map(d -> length(d.t_eval), remove(NoInterpolationDimension, interp_dims)) + # Get the size of NoInterpolationDimension dims + nointerp_size = get_output_size(interp) + # Insert the nointerp sizes back into the interp_size tuple + sze = insertat(NoInterpolationDimension, nointerp_size, interp_size, interp_dims) end make_zero!!(::T) where {T <: Number} = zero(T) @@ -174,24 +185,70 @@ function get_ndrange(interp::NDInterpolation) I = map(interp.interp_dims) do d d isa NoInterpolationDimension ? nothing : length(d.t_eval) end - return _remove(Nothing, I...) + return remove(Nothing, I) end -# Remove objects of type T from splatted args (taken from DimensionalData.jl) -Base.@assume_effects :foldable _remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) -Base.@assume_effects :foldable _remove(::Type{T}, ::T, xs...) where T = _remove(T, xs...) -Base.@assume_effects :foldable _remove(::Type) = () -# Insert x in `out` where `m isa T`, otherwise output one of in for each m -# out must exactly match the number of ms or this will error +# Some tuple handling primitives +# TODO: make sure these compile away completely + +# Remove objects of type `T` from `in` (reworked from DimensionalData.jl) +remove(::Type{T}, in) where T = _remove(T, in...) +_remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) +_remove(::Type{T}, x::T, xs...) where T = _remove(T, xs...) +_remove(::Type) = () + +# Keep only objects of type `T` from `in` +keep(::Type{T}, in) where T = _keep(T, in...) +_keep(::Type{T}, x, xs...) where T = _keep(T, xs...) +_keep(::Type{T}, x::T, xs...) where T = (x, _keep(T, xs...)...) +_keep(::Type) = () + +# Remove values from `in` where `matches` are of type `T` +function removeat(::Type{T}, in::Tuple, matches::Tuple) where T + @assert length(in) == length(matches) + _removeat(T, in, matches...) +end +# If `!(m isa T)` take from `in` +_removeat(::Type{T}, in::Tuple{<:Any,Vararg}, m, ms...) where T = + (first(in), _removeat(T, Base.tail(in), ms...)...) +# If `m isa T` remove +_removeat(::Type{T}, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = + _removeat(T, Base.tail(in), ms...) +# `in` can be empty if there are no remaining `m` +_removeat(::Type, in::Tuple{}) = () + +# Keep only values from `in` where `matches` are of type `T` +function keepat(::Type{T}, in::Tuple, matches::Tuple) where T + @assert length(in) == length(matches) + _keepat(T, in, matches...) +end +# If `!(m isa T)` disgaurd +_keepat(::Type{T}, in::Tuple{<:Any,Vararg}, m, ms...) where T = + _keepat(T, Base.tail(in), ms...) +# If `m isa T` keep +_keepat(::Type{T}, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = + (first(in), _keepat(T, Base.tail(in), ms...)...) +# `in` can be empty if there are no remaining `m` +_keepat(::Type, in::Tuple{}) = () + +# Insert x into `in` where `matches` are of type `T`, +# otherwise output one of `in` for each `m`. +insertat(::Type{T}, x, in::Tuple, matches::Tuple) where T = + _insertat(T, x, in, matches...) # If `!(m isa T)` take from `in` -Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, out::Tuple, m, ms...) where T = - _insertat(T, x, Base.tail(in), (out..., first(in)), ms...) +_insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, m, ms...) where T = + (first(in), _insertat(T, x, Base.tail(in), ms...)...) # If `m isa T` insert `x` -Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, out::Tuple, m::T, ms...) where T = - _insertat(T, x, in, (out..., x), ms...) -# `in` can be empty if all trailing ms are T -Base.@assume_effects :foldable _insertat(::Type{T}, x, in::Tuple, out::Tuple, m::T, ms::T...) where T = - _insertat(T, x, in, (out..., x), ms...) -# `in` can also be empty if there are no remaining `m` -Base.@assume_effects :foldable _insertat(::Type, x, in::Tuple{}, out::Tuple) = out +_insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = + (x, _insertat(T, x, in, ms...)...) +# For Tuple x we insert the first +_insertat(::Type{T}, xs::Tuple, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = + (first(xs), _insertat(T, Base.tail(xs), in, ms...)...) +# `in` can be empty if there are no remaining `m` +_insertat(::Type, x, in::Tuple{}) = () +# `in` can also be empty if all trailing ms are T +_insertat(::Type{T}, x, in::Tuple{}, m::T, ms::T...) where T = + (x, _insertat(T, x, in, ms...)...) +_insertat(::Type{T}, xs::Tuple, in::Tuple{}, m::T, ms::T...) where T = + (first(xs), _insertat(T, Base.tail(xs), in, ms...)...) diff --git a/src/spline_utils.jl b/src/spline_utils.jl index cff7f5c..e2fe444 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -104,6 +104,17 @@ function get_basis_function_values( multi_point_index, :, derivative_order + 1) end +# Get all basis function values to evaluate a BSpline interpolation in t +function get_basis_function_values_all( + A::NDInterpolation, + ts::Tuple, + idx::Tuple, + derivative_orders::Tuple, + multi_point_index +) + map(get_basis_function_values, A.interp_dims, ts, idx, derivative_orders, multi_point_index) +end + function set_basis_function_eval!(itp_dim::BSplineInterpolationDimension)::Nothing backend = get_backend(itp_dim.t_eval) basis_function_eval_kernel(backend)( diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index ecd2155..3de3308 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -1,10 +1,24 @@ using DataInterpolationsND using Random -using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache - +using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache, + keepat, removeat, insertat, remove, keep + +@testset "tuple primitives" begin + # These help with handling NoInterpolationDimension and its consequences + @test insertat(Nothing, 1, (2,), (nothing, 4, nothing)) === (1, 2, 1) + @test insertat(Float64, missing, (1,), (8, 9.0)) === (1, missing) + @test keepat(Nothing, (1, 2, 3), (nothing, 4, nothing)) === (1, 3) + @test keepat(Float64, (1, 2,), (8, 9.0)) === (2,) + @test removeat(Nothing, (1, 2, 3), (nothing, 4, nothing)) === (2,) + @test removeat(Float64, (1, 2,), (8, 9.0)) === (1,) + @test remove(Nothing, (1, nothing, 3)) === (1, 3) + @test keep(Int, (1, nothing, 3)) === (1, 3) +end -function test_globally_constant( - ID::Type{<:AbstractInterpolationDimension}; args1 = (), args2 = (), kwargs1 = (), +function test_globally_constant(ID::Type{<:AbstractInterpolationDimension}; + args1 = (), + args2 = (), + kwargs1 = (), kwargs2 = (), cache = EmptyCache(), test_derivatives = true @@ -25,7 +39,9 @@ function test_globally_constant( ) itp = NDInterpolation(u, itp_dims; cache) - @test all(x -> isapprox(x, 2.0; atol = 1e-10), eval_grid(itp)) + + grid = eval_grid(itp) + @test all(x -> isapprox(x, 2.0; atol = 1e-10), grid) if test_derivatives @test all( x -> isapprox(x, 0.0; atol = 1e-10), eval_grid(itp, derivative_orders = (1, 0))) @@ -38,7 +54,7 @@ function test_globally_constant( ID(t1, args1...; t_eval = t1[1:(end - 1)] + diff(t1) / 2, kwargs1...), ID(t2, args2...; t_eval = t2[1:(end - 1)] + diff(t2) / 2, kwargs2...) ) - itp = NDInterpolation(u, itp_dims) + itp = NDInterpolation(u, itp_dims; cache) @test all(x -> isapprox(x, 2.0; atol = 1e-10), eval_grid(itp)) if test_derivatives @test all( @@ -49,7 +65,7 @@ function test_globally_constant( end function test_analytic(itp::NDInterpolation{<:Any,N_in}, f) where {N_in} - used_interp_dims = DataInterpolationsND._remove(NoInterpolationDimension, itp.interp_dims...) + used_interp_dims = remove(NoInterpolationDimension, itp.interp_dims) # Evaluation in data points ts = map(d -> d.t, used_interp_dims) for t in Iterators.product(ts...) @@ -64,6 +80,10 @@ function test_analytic(itp::NDInterpolation{<:Any,N_in}, f) where {N_in} end end +@testset "Constant Interpolation" begin + test_globally_constant(ConstantInterpolationDimension) +end + @testset "Linear Interpolation" begin test_globally_constant(LinearInterpolationDimension) @@ -80,6 +100,7 @@ end itp = NDInterpolation(u, itp_dims) test_analytic(itp, f) + itp_dims = ( NoInterpolationDimension(), LinearInterpolationDimension(t2; t_eval = t2 ./ 2) @@ -158,16 +179,14 @@ end end @testset "Mixed Interpolation" begin - t1 = cumsum(rand(3)) t2 = cumsum(rand(4)) t3 = collect(0:(π / 2):(2π)) - t_eval = collect(range(0, 2π, length = 100)) u = zeros(3, 4, 5) itp_dims = ( NoInterpolationDimension(), - LinearInterpolationDimension(t2; t_eval = t2 ./ 2), - BSplineInterpolationDimension(t3, 2; multiplicities, t_eval), + LinearInterpolationDimension(t2; t_eval=t2), + BSplineInterpolationDimension(t3, 1; t_eval=t3), ) itp = NDInterpolation(u, itp_dims) eval_grid(itp) From 9061c83d580b6cd2ec31e4dafb93a6522c566f21 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 23 Oct 2025 13:18:09 +1100 Subject: [PATCH 11/32] fix ConstantInterpolationDimension nans --- src/interpolation_methods.jl | 38 +++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 76949b1..02a12f9 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -8,7 +8,8 @@ function _interpolate!( ) where {N,N_in,N_out} (; interp_dims, cache, u) = A - check_derivative_order(interp_dims, derivative_orders) || return out + out, valid_derivative_orders = check_derivative_order(interp_dims, derivative_orders, ts, out) + valid_derivative_orders || return out if isnothing(multi_point_index) multi_point_index = map(_ -> nothing, interp_dims) end @@ -46,19 +47,28 @@ function _interpolate!( return out end -check_derivative_order(dims::Tuple, derivative_orders::Tuple) = - all(map(check_derivative_order, dims, derivative_orders)) -check_derivative_order(::LinearInterpolationDimension, d_o) = d_o <= 1 -check_derivative_order(::ConstantInterpolationDimension, d_o) = d_o <= 0 -check_derivative_order(::AbstractInterpolationDimension, d_o) = true -# TODO how to handle this -# if derivative_order > 0 -# return if any(i -> !isempty(searchsorted(A.interp_dims[i].t, t[i])) -# typed_nan(out) -# else -# out -# end -# end +function check_derivative_order(dims::Tuple, derivative_orders::Tuple, ts::Tuple, out) + itr = map(tuple, dims, derivative_orders, ts) + # Fold over itr for all dims, combining out and valid + foldl(itr; init=(out, true)) do (acc_out, acc_valid), (d, d_o, t) + dim_out, dim_valid = check_derivative_order(d, d_o, t, acc_out) + dim_out, dim_valid & acc_valid + end +end +check_derivative_order(::AbstractInterpolationDimension, d_o, t, out) = (out, true) +check_derivative_order(::LinearInterpolationDimension, d_o, t, out) = (out, d_o <= 1) +function check_derivative_order(d::ConstantInterpolationDimension, d_o, t, out) + if d_o > 0 + # Check if t is on the boundary between constant steps and if so return nans + return if isempty(searchsorted(d.t, t)) + (out, false) + else + (typed_nan(out), false) + end + else + (out, true) + end +end function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) t₁ = d.t[i] From 0bfe73c0bb8d83d828e84744c25d39f553123f71 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 23 Oct 2025 13:23:03 +1100 Subject: [PATCH 12/32] remove constant --- test/test_interpolations.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 3de3308..6497a98 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -80,10 +80,6 @@ function test_analytic(itp::NDInterpolation{<:Any,N_in}, f) where {N_in} end end -@testset "Constant Interpolation" begin - test_globally_constant(ConstantInterpolationDimension) -end - @testset "Linear Interpolation" begin test_globally_constant(LinearInterpolationDimension) From 42e6d68bd6193b46efd899974769a85cb1d5d3b9 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Thu, 23 Oct 2025 13:44:47 +1100 Subject: [PATCH 13/32] fix eval grid checks --- src/interpolation_parallel.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 98ae8f8..410d530 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -37,14 +37,14 @@ function eval_unstructured!( derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) ) where {N_in} validate_derivative_orders(derivative_orders, interp; multi_point = true) + validate_output_size(out, interp) backend = get_backend(out) - @assert all(i -> length(interp.interp_dims[i].t_eval) == size(out, 1), N_in) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." - @assert size(out)[2:end]==get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." + # TODO this may be broken but it isn't tested + @assert all(d -> length(d.t_eval) == size(out, 1), remove(NoInterpolationDimension, interp.interp_dims)) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." eval_kernel(backend)( out, interp, derivative_orders, - false, ndrange = size(out, 1) ) synchronize(backend) @@ -87,27 +87,37 @@ function eval_grid!( interp::NDInterpolation{N,N_in}; derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) ) where {N,N_in} + validate_t_eval_lengths(out, interp) + validate_output_size(out, interp) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) - # used_interp_dims = remove(NoInterpolationDimension, interp.interp_dims) - # @assert all(ntuple(i -> size(out, i) == length(used_interp_dims[i].t_eval), N_in)) "The length must match the t_eval of the corresponding interpolation dimension." - # @assert size(out) == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." eval_kernel(backend)( out, interp, derivative_orders, - true, ndrange = get_ndrange(interp) ) synchronize(backend) return out end +function validate_t_eval_lengths(out, interp) + (; interp_dims) = interp + interp_sizes = removeat(NoInterpolationDimension, size(out), interp_dims) + interp_t_eval_lengths = map(d -> length(d.t_eval), remove(NoInterpolationDimension, interp_dims)) + all(map(==, interp_sizes, interp_t_eval_lengths)) || + throw(ArgumentError("The length must match the t_eval of the corresponding interpolation dimension.")) +end + +function validate_output_size(out, interp) + keepat(NoInterpolationDimension, size(out), interp.interp_dims) == get_output_size(interp) || + throw(ArumentError("The size of the NoInterpolationDimension dimensions of `out` must be the same as the output size of the interpolation.")) +end + @kernel function eval_kernel( out, A::NDInterpolation{N, N_in, N_out}, derivative_orders, - eval_grid, ) where {N, N_in, N_out} # This kernel is only over interpolated dimensions, we need # to insert fillers to match the number of dimensions in the data From 60c36a2019638a5b8e8fbc264b8aa84d16e7c7fe Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 24 Oct 2025 01:32:32 +1100 Subject: [PATCH 14/32] simplify and minimise changes --- src/DataInterpolationsND.jl | 3 +++ src/interpolation_dimensions.jl | 2 +- src/interpolation_methods.jl | 12 ++++++------ src/interpolation_parallel.jl | 16 ++++++++-------- src/interpolation_utils.jl | 34 +++++++++++++++++++++++++-------- src/spline_utils.jl | 3 ++- 6 files changed, 46 insertions(+), 24 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index bba8835..260745c 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -72,10 +72,12 @@ include("plot_rec.jl") function (interp::NDInterpolation)(t_args::Vararg{Number}; kwargs...) interp(t_args; kwargs...) end + function (interp::NDInterpolation)( out::AbstractArray, t_args::Vararg{Number}; kwargs...) interp(out, t_args; kwargs...) end + # In place single input evaluation function (interp::NDInterpolation{N,N_in,N_out})( out::Union{Number, AbstractArray{<:Number, N_out}}, @@ -87,6 +89,7 @@ function (interp::NDInterpolation{N,N_in,N_out})( idx = get_idx(interp.interp_dims, t) return _interpolate!(out, interp, t, idx, derivative_orders, nothing) end + # Out of place single input evaluation function (interp::NDInterpolation)(t::Tuple{Vararg{Number}}; kwargs...) out = make_out(interp, t) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 8c8d12b..5c4107a 100644 --- a/src/interpolation_dimensions.jl +++ b/src/interpolation_dimensions.jl @@ -164,8 +164,8 @@ function BSplineInterpolationDimension( synchronize(backend) idx_eval = similar(t_eval, Int) - s = (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, diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 02a12f9..b08ce68 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -16,10 +16,10 @@ function _interpolate!( out = make_zero!!(out) denom = zero(eltype(u)) # Setup - space = map(iteration_space, interp_dims) + stencils = map(stencil, interp_dims) preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) - for I in Iterators.product(space...) + for I in Iterators.product(stencils...) scaling = map(scale, interp_dims, preparations, I) J = map(index, interp_dims, ts, idx, I) if cache isa EmptyCache @@ -86,10 +86,10 @@ function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point return (; basis_function_values) end -iteration_space(::LinearInterpolationDimension) = (false, true) -iteration_space(::ConstantInterpolationDimension) = 1 -iteration_space(::NoInterpolationDimension) = 1 -iteration_space(d::BSplineInterpolationDimension) = 1:d.degree + 1 +stencil(::LinearInterpolationDimension) = (false, true) +stencil(::ConstantInterpolationDimension) = 1 +stencil(::NoInterpolationDimension) = 1 +stencil(d::BSplineInterpolationDimension) = 1:d.degree + 1 function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) (; t, t₁, t₂, t_vol_inv, derivative_order) = prep diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 410d530..9623b83 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -116,11 +116,10 @@ end @kernel function eval_kernel( out, - A::NDInterpolation{N, N_in, N_out}, + A, # @Const(A), TODO: somehow this now hits a bug in KernelAbstractions where elsize is not defined for Const derivative_orders, -) where {N, N_in, N_out} - # This kernel is only over interpolated dimensions, we need - # to insert fillers to match the number of dimensions in the data +) + N_out = length(keep(NoInterpolationDimension, A.interp_dims)) I = @index(Global, NTuple) k = insertat(NoInterpolationDimension, Colon(), I, A.interp_dims) @@ -128,11 +127,12 @@ end idx_eval = map(get_idx_eval, A.interp_dims, k) if iszero(N_out) - dest = make_out(A, t_eval) - out[I...] = _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) + out[k...] = _interpolate!( + make_out(A, t_eval), A, t_eval, idx_eval, derivative_orders, k) else - dest = view(out, k...) - _interpolate!(dest, A, t_eval, idx_eval, derivative_orders, k) + _interpolate!( + view(out, k...), + A, t_eval, idx_eval, derivative_orders, k) end end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 969b4c6..0dcf708 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -2,7 +2,9 @@ trivial_range(i::Integer) = i:i Base.length(itp_dim::AbstractInterpolationDimension) = length(itp_dim.t) -function validate_derivative_order(derivative_orders::NTuple, A::NDInterpolation; +function validate_derivative_order( + derivative_orders::NTuple, + A::NDInterpolation; multi_point::Bool = false ) map(derivative_orders, A.interp_dims) do d_o, d @@ -39,7 +41,6 @@ validate_derivative_order_by_cache(::NURBSWeights, derivative_order) = validate_derivative_order_by_cache(::Any, derivative_order) = @assert derivative_order >= 0 "Derivative orders must me non-negative." - function validate_t(t) @assert t isa AbstractVector{<:Number} "t must be an AbstractVector with number like elements." @assert all(>(0), diff(t)) "The elements of t must be sorted and unique." @@ -52,12 +53,18 @@ validate_size_u(interp_dim::NoInterpolationDimension, ax::AbstractRange) = nothi function validate_size_u(interp_dim::AbstractInterpolationDimension, ax::AbstractRange) @assert length(interp_dim) == length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." end -function validate_size_u(interp_dim::BSplineInterpolationDimension, ax::AbstractRange) + +function validate_size_u( + interp_dim::BSplineInterpolationDimension, + ax::AbstractRange +) expected_size = get_n_basis_functions(interp_dim) @assert expected_size == length(ax) "Expected the size to be $expected_size based on the BSplineInterpolation properties, got $(length(ax))." end -function validate_cache(cache::AbstractInterpolationCache, dims::Tuple, u::AbstractArray) +function validate_cache( + cache::AbstractInterpolationCache, dims::Tuple, u::AbstractArray +) ntuple(length(dims)) do n validate_cache(cache, dims[n], u, n) end @@ -74,7 +81,9 @@ function validate_cache( size_expected = size(u, n) @assert size(nurbs_weights.weights, n) == size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected), got $(size(nurbs_weights.weights, n))." end -function validate_cache(::gType, ::ID, ::AbstractArray, ::Int) where {gType,ID<:AbstractInterpolationDimension} + +function validate_cache( + ::gType, ::ID, ::AbstractArray, ::Int) where {gType,ID<:AbstractInterpolationDimension} @error("Interpolation dimension type $ID is not compatible with global cache type $gType.") end @@ -97,6 +106,7 @@ function grid_size(interp::NDInterpolation) end make_zero!!(::T) where {T <: Number} = zero(T) + function make_zero!!(v::T) where {T <: AbstractArray} v .= 0 v @@ -147,12 +157,18 @@ function get_idx( clamp(searchsortedlast(t, t_eval) + idx_shift, lb, ub) end end + get_idx(::NoInterpolationDimension, t_eval::Number) = nothing -function get_idx(interp_dims::NTuple{N}, t::Tuple{Vararg{Number, N}}) where N +function get_idx( + interp_dims::NTuple{N}, + t::Tuple{Vararg{Number, N}}; +) where {N} map(get_idx, interp_dims, t) end -function set_eval_idx!(interp_dim::AbstractInterpolationDimension) +function set_eval_idx!( + interp_dim::AbstractInterpolationDimension, +) backend = get_backend(interp_dim.t) if !isempty(interp_dim.t_eval) set_idx_kernel(backend)( @@ -163,7 +179,9 @@ function set_eval_idx!(interp_dim::AbstractInterpolationDimension) synchronize(backend) end -@kernel function set_idx_kernel(interp_dim::AbstractInterpolationDimension) +@kernel function set_idx_kernel( + interp_dim +) i = @index(Global, Linear) interp_dim.idx_eval[i] = get_idx(interp_dim, interp_dim.t_eval[i]) end diff --git a/src/spline_utils.jl b/src/spline_utils.jl index e2fe444..be78a8f 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -60,7 +60,7 @@ function get_basis_function_values( t::Number, idx::Integer, derivative_order::Integer, - multi_point_index::Union{Nothing,Colon}, # TODO why both + multi_point_index::Nothing, ) (; degree, knots_all) = itp_dim T = promote_type(typeof(t), eltype(itp_dim.basis_function_eval)) @@ -91,6 +91,7 @@ function get_basis_function_values( basis_function_values[1:degree_plus_1] end + # Get the basis function values for one point in an # unstructured multi point evaluation (given by the scalar multi point index) function get_basis_function_values( From cb2864d2e4d9b40694d3636183ce40237e3455af Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 24 Oct 2025 14:32:15 +1100 Subject: [PATCH 15/32] rename scale to weight --- src/interpolation_methods.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index b08ce68..b585100 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -13,20 +13,20 @@ function _interpolate!( if isnothing(multi_point_index) multi_point_index = map(_ -> nothing, interp_dims) end + # Setup out = make_zero!!(out) denom = zero(eltype(u)) - # Setup stencils = map(stencil, interp_dims) preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) for I in Iterators.product(stencils...) - scaling = map(scale, interp_dims, preparations, I) J = map(index, interp_dims, ts, idx, I) + weights = map(weight, interp_dims, preparations, I) if cache isa EmptyCache - product = prod(scaling) + product = prod(weights) else K = removeat(NoInterpolationDimension, J, interp_dims) - product = cache.weights[K...] * prod(scaling) + product = cache.weights[K...] * prod(weights) denom += product end if iszero(N_out) @@ -91,7 +91,7 @@ stencil(::ConstantInterpolationDimension) = 1 stencil(::NoInterpolationDimension) = 1 stencil(d::BSplineInterpolationDimension) = 1:d.degree + 1 -function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) +function weight(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) (; t, t₁, t₂, t_vol_inv, derivative_order) = prep if right_point iszero(derivative_order) ? t - t₁ : one(t) @@ -99,9 +99,9 @@ function scale(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bo iszero(derivative_order) ? t₂ - t : -one(t) end * t_vol_inv end -scale(::ConstantInterpolationDimension, prep::NamedTuple, i) = 1 -scale(::NoInterpolationDimension, prep::NamedTuple, i) = 1 -scale(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_function_values[i] +weight(::ConstantInterpolationDimension, prep::NamedTuple, i) = 1 +weight(::NoInterpolationDimension, prep::NamedTuple, i) = 1 +weight(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_function_values[i] index(::LinearInterpolationDimension, t, idx, i) = idx + i index(d::ConstantInterpolationDimension, t, idx, i) = t >= d.t[end] ? length(d.t) : idx[i] From c010f3048e860228810584659eb7180d20f97e95 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 24 Oct 2025 14:54:04 +1100 Subject: [PATCH 16/32] cleanup grid_size --- src/interpolation_utils.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 0dcf708..1533e46 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -96,13 +96,12 @@ end function grid_size(interp::NDInterpolation) (; interp_dims) = interp - # TODO: put this in a function, but # Get the size of dims that are not NoInterpolationDimension interp_size = map(d -> length(d.t_eval), remove(NoInterpolationDimension, interp_dims)) # Get the size of NoInterpolationDimension dims nointerp_size = get_output_size(interp) # Insert the nointerp sizes back into the interp_size tuple - sze = insertat(NoInterpolationDimension, nointerp_size, interp_size, interp_dims) + return insertat(NoInterpolationDimension, nointerp_size, interp_size, interp_dims) end make_zero!!(::T) where {T <: Number} = zero(T) From 5a87b171c7a02309efc053fbfa7a1531cec0c871 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 24 Oct 2025 17:30:06 +1100 Subject: [PATCH 17/32] explicity multi_point_index --- src/DataInterpolationsND.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 260745c..289874a 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -87,7 +87,8 @@ function (interp::NDInterpolation{N,N_in,N_out})( validate_size_u(interp, out) validate_derivative_order(derivative_orders, interp) idx = get_idx(interp.interp_dims, t) - return _interpolate!(out, interp, t, idx, derivative_orders, nothing) + multi_point_index = nothing + return _interpolate!(out, interp, t, idx, derivative_orders, multi_point_index) end # Out of place single input evaluation From 31fe7ea6e7b9a842004ea92e9aab3cefe8b280c3 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 31 Oct 2025 16:13:16 +1100 Subject: [PATCH 18/32] bugfix and test nointerp itp() --- src/DataInterpolationsND.jl | 16 ++++--- src/interpolation_methods.jl | 85 +++++++++++++++-------------------- src/interpolation_parallel.jl | 25 ++++++++--- src/interpolation_utils.jl | 23 +++++----- test/test_interpolations.jl | 4 +- 5 files changed, 78 insertions(+), 75 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 289874a..c1d6ab8 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -67,6 +67,7 @@ include("interpolation_utils.jl") include("interpolation_methods.jl") include("interpolation_parallel.jl") include("plot_rec.jl") +include("interp_array.jl") # Multiple `t` arguments to tuple (can these 2 be done in 1?) function (interp::NDInterpolation)(t_args::Vararg{Number}; kwargs...) @@ -81,14 +82,19 @@ end # In place single input evaluation function (interp::NDInterpolation{N,N_in,N_out})( out::Union{Number, AbstractArray{<:Number, N_out}}, - t::Tuple{Vararg{Number, N}}; - derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) + t::Tuple{Vararg{Any, N_in}}; + derivative_orders::Tuple{Vararg{Integer,N_in}} = ntuple(_ -> 0, N_in) ) where {N,N_in,N_out} validate_size_u(interp, out) - validate_derivative_order(derivative_orders, interp) - idx = get_idx(interp.interp_dims, t) + # 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. + t_all = insertat(NoInterpolationDimension, Colon(), t, interp.interp_dims) + idx_all = get_idx(interp.interp_dims, t_all) + d_o_all = insertat(NoInterpolationDimension, 0, derivative_orders, interp.interp_dims) + validate_derivative_order(d_o_all, interp) multi_point_index = nothing - return _interpolate!(out, interp, t, idx, derivative_orders, multi_point_index) + return _interpolate!(out, interp, t_all, idx_all, d_o_all, multi_point_index) end # Out of place single input evaluation diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index b585100..f2b3e60 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,11 +1,11 @@ -function _interpolate!( - out::Union{Number, AbstractArray{<:Any, N_out}}, - A::NDInterpolation{N, N_in, N_out}, +Base.@propagate_inbounds function _interpolate!( + out::Union{Number, AbstractArray}, + A::NDInterpolation{N}, ts::Tuple{Vararg{Any, N}}, idx::Tuple{Vararg{Any ,N}}, derivative_orders::Tuple{Vararg{Any, N}}, - multi_point_index -) where {N,N_in,N_out} + multi_point_index, +) where N (; interp_dims, cache, u) = A out, valid_derivative_orders = check_derivative_order(interp_dims, derivative_orders, ts, out) @@ -13,34 +13,36 @@ function _interpolate!( if isnothing(multi_point_index) multi_point_index = map(_ -> nothing, interp_dims) end + # Setup - out = make_zero!!(out) - denom = zero(eltype(u)) + out = make_zero!!(out) # TODO remove this stencils = map(stencil, interp_dims) - preparations = map(prepare, interp_dims, derivative_orders, multi_point_index, ts, idx) + stencil_weights = map(weights, interp_dims, derivative_orders, multi_point_index, ts, idx) + denom = zero(eltype(u)) + # TODO this can be a single unrolled broadcast rather than a loop of .+= for I in Iterators.product(stencils...) J = map(index, interp_dims, ts, idx, I) - weights = map(weight, interp_dims, preparations, I) - if cache isa EmptyCache - product = prod(weights) - else + product = prod(map(getindex, stencil_weights, I)) + + if !(cache isa EmptyCache) K = removeat(NoInterpolationDimension, J, interp_dims) - product = cache.weights[K...] * prod(weights) + product *= cache.weights[K...] denom += product end - if iszero(N_out) - out += product * u[J...] + + if out isa AbstractArray + out .+= product .* u[J...] else - out .+= product .* view(u, J...) + out += product * u[J...] end end if !(cache isa EmptyCache) - if iszero(N_out) - out /= denom - else + if out isa AbstractArray out ./= denom + else + out /= denom end end @@ -70,40 +72,27 @@ function check_derivative_order(d::ConstantInterpolationDimension, d_o, t, out) end end -function prepare(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) - t₁ = d.t[i] - t₂ = d.t[i + 1] - t_vol_inv = inv(t₂ - t₁) - return (; t, t₁, t₂, t_vol_inv, derivative_order) -end -prepare(::ConstantInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) -prepare(::NoInterpolationDimension, derivative_orders, multi_point_index, t, i) = (;) -function prepare(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) - # TODO the dim_in arg isn't really needed, so drop it. Currently just 0 - basis_function_values = get_basis_function_values( - d, t, i, derivative_order, multi_point_index - ) - return (; basis_function_values) -end - -stencil(::LinearInterpolationDimension) = (false, true) +stencil(::LinearInterpolationDimension) = (1, 2) stencil(::ConstantInterpolationDimension) = 1 stencil(::NoInterpolationDimension) = 1 stencil(d::BSplineInterpolationDimension) = 1:d.degree + 1 -function weight(::LinearInterpolationDimension, prep::NamedTuple, right_point::Bool) - (; t, t₁, t₂, t_vol_inv, derivative_order) = prep - if right_point - iszero(derivative_order) ? t - t₁ : one(t) - else - iszero(derivative_order) ? t₂ - t : -one(t) - end * t_vol_inv +# Precalculate weights +function weights(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) + t₁ = d.t[i] + t₂ = d.t[i + 1] + t_vol_inv = inv(t₂ - t₁) + a = (iszero(derivative_order) ? t₂ - t : -one(t)) * t_vol_inv + b = (iszero(derivative_order) ? t - t₁ : one(t)) * t_vol_inv + return (a, b) end -weight(::ConstantInterpolationDimension, prep::NamedTuple, i) = 1 -weight(::NoInterpolationDimension, prep::NamedTuple, i) = 1 -weight(::BSplineInterpolationDimension, prep::NamedTuple, i) = prep.basis_function_values[i] +weights(::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 +weights(::NoInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 +weights(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) = + get_basis_function_values(d, t, i, derivative_order, multi_point_index) -index(::LinearInterpolationDimension, t, idx, i) = idx + i +index(::LinearInterpolationDimension, t, idx, i) = idx + i - 1 +# TODO: this should happen outside of the loop index(d::ConstantInterpolationDimension, t, idx, i) = t >= d.t[end] ? length(d.t) : idx[i] -index(::NoInterpolationDimension, t, idx, i) = Colon() +index(::NoInterpolationDimension, t, idx, i) = idx index(d::BSplineInterpolationDimension, t, idx, i) = idx + i - d.degree - 1 diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 9623b83..24323fb 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -91,10 +91,12 @@ function eval_grid!( validate_output_size(out, interp) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) + no_interp_inds = map(_ -> Colon(), remove(NoInterpolationDimension, interp.interp_dims)) eval_kernel(backend)( out, interp, derivative_orders, + no_interp_inds, ndrange = get_ndrange(interp) ) synchronize(backend) @@ -118,16 +120,17 @@ end out, A, # @Const(A), TODO: somehow this now hits a bug in KernelAbstractions where elsize is not defined for Const derivative_orders, + no_interp_inds, ) N_out = length(keep(NoInterpolationDimension, A.interp_dims)) I = @index(Global, NTuple) - k = insertat(NoInterpolationDimension, Colon(), I, A.interp_dims) + k = insertat(NoInterpolationDimension, no_interp_inds, I, A.interp_dims) t_eval = map(get_t_eval, A.interp_dims, k) idx_eval = map(get_idx_eval, A.interp_dims, k) if iszero(N_out) - out[k...] = _interpolate!( + @inbounds out[k...] = _interpolate!( make_out(A, t_eval), A, t_eval, idx_eval, derivative_orders, k) else _interpolate!( @@ -136,8 +139,16 @@ end end end -get_t_eval(d, i) = d.t_eval[i] -get_t_eval(d::NoInterpolationDimension, i) = nothing - -get_idx_eval(d, i) = d.idx_eval[i] -get_idx_eval(d::NoInterpolationDimension, i) = nothing +get_t_eval(d::AbstractInterpolationDimension, i) = d.t_eval[i] +get_t_eval(d::AbstractInterpolationDimension, i::Colon) = d.t_eval +get_t_eval(d::AbstractInterpolationDimension, i::AbstractVector) = d.t_eval[i] +get_t_eval(::NoInterpolationDimension, i) = nothing +get_t_eval(::NoInterpolationDimension, i::Colon) = nothing +get_t_eval(::NoInterpolationDimension, i::AbstractVector) = nothing + +get_idx_eval(d::AbstractInterpolationDimension, i) = d.idx_eval[i] +get_idx_eval(d::AbstractInterpolationDimension, i::AbstractVector) = d.idx_eval[i] +get_idx_eval(d::AbstractInterpolationDimension, i::Colon) = d.idx_eval +get_idx_eval(::NoInterpolationDimension, i) = i +get_idx_eval(::NoInterpolationDimension, i::Colon) = i +get_idx_eval(::NoInterpolationDimension, i::AbstractVector) = i diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 1533e46..322fdf0 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -46,9 +46,10 @@ function validate_t(t) @assert all(>(0), diff(t)) "The elements of t must be sorted and unique." end -validate_size_u(::NDInterpolation{N,N,0}, u::Number) where N = nothing -validate_size_u(interp::NDInterpolation, u::AbstractArray) = validate_size_u(interp.interp_dims, axes(u)) -validate_size_u(interp_dims::Tuple, u::AbstractArray) = map(validate_size_u, interp_dims, axes(u)) +validate_size_u(interp::NDInterpolation, u) = validate_size_u(interp.interp_dims, u) +validate_size_u(interp_dims::Tuple, u::Number) = nothing +validate_size_u(interp_dims::Tuple, u::AbstractArray) = validate_size_u(interp_dims, axes(u)) +validate_size_u(interp_dims::Tuple, ax::Tuple) = map(validate_size_u, interp_dims, ax) validate_size_u(interp_dim::NoInterpolationDimension, ax::AbstractRange) = nothing function validate_size_u(interp_dim::AbstractInterpolationDimension, ax::AbstractRange) @assert length(interp_dim) == length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." @@ -106,8 +107,8 @@ end make_zero!!(::T) where {T <: Number} = zero(T) -function make_zero!!(v::T) where {T <: AbstractArray} - v .= 0 +function make_zero!!(v::AbstractArray) + fill!(v, zero(eltype(v))) v end @@ -115,7 +116,7 @@ function make_out( interp::NDInterpolation{<:Any,N_in, 0}, t::NTuple{N_in, >:Number} ) where {N_in} - zero(promote_type(eltype(interp.u), map(typeof, t)...)) + zero(eltype(interp.u)) end function make_out( interp::NDInterpolation{<:Any,N_in}, @@ -156,14 +157,10 @@ function get_idx( clamp(searchsortedlast(t, t_eval) + idx_shift, lb, ub) end end - -get_idx(::NoInterpolationDimension, t_eval::Number) = nothing -function get_idx( - interp_dims::NTuple{N}, - t::Tuple{Vararg{Number, N}}; -) where {N} +# t_eval must already be an index for NoInterpolationDimension +get_idx(::NoInterpolationDimension, t_eval::Union{Colon,Int,AbstractArray{Int}}) = t_eval +get_idx(interp_dims::Tuple{Vararg{Any,N}}, t::Tuple{Vararg{Any,N}}) where N = map(get_idx, interp_dims, t) -end function set_eval_idx!( interp_dim::AbstractInterpolationDimension, diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 6497a98..36320ca 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -175,15 +175,15 @@ end end @testset "Mixed Interpolation" begin - t1 = cumsum(rand(3)) t2 = cumsum(rand(4)) t3 = collect(0:(π / 2):(2π)) - u = zeros(3, 4, 5) + u = ones(3, 4, 5) itp_dims = ( NoInterpolationDimension(), LinearInterpolationDimension(t2; t_eval=t2), BSplineInterpolationDimension(t3, 1; t_eval=t3), ) itp = NDInterpolation(u, itp_dims) + @test isapprox(itp(1.0, 4.0), [1.0, 1.0, 1.0]) eval_grid(itp) end From 471b1e590f9c1cbbb03ccd41285f542cd28db903 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 31 Oct 2025 16:31:14 +1100 Subject: [PATCH 19/32] cleanup --- src/DataInterpolationsND.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index c1d6ab8..fa2a87d 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -67,7 +67,6 @@ include("interpolation_utils.jl") include("interpolation_methods.jl") include("interpolation_parallel.jl") include("plot_rec.jl") -include("interp_array.jl") # Multiple `t` arguments to tuple (can these 2 be done in 1?) function (interp::NDInterpolation)(t_args::Vararg{Number}; kwargs...) From 09e20e22df1d80f21f17fdf4a7f18700eb359a8e Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 31 Oct 2025 16:42:02 +1100 Subject: [PATCH 20/32] run JuliaFormatter --- src/DataInterpolationsND.jl | 26 +++--- src/interpolation_methods.jl | 20 +++-- src/interpolation_parallel.jl | 29 ++++--- src/interpolation_utils.jl | 144 +++++++++++++++++++--------------- src/spline_utils.jl | 9 ++- test/test_interpolations.jl | 31 ++++---- 6 files changed, 146 insertions(+), 113 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index fa2a87d..b3a2ffb 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -23,7 +23,7 @@ the size of `u` along that dimension must match the length of `t` of the corresp """ struct NDInterpolation{ N, - N_in, + N_in, N_out, gType <: AbstractInterpolationCache, D, @@ -32,7 +32,7 @@ struct NDInterpolation{ u::uType interp_dims::D cache::gType - function NDInterpolation(u::AbstractArray{<:Any,N}, interp_dims, cache) where N + function NDInterpolation(u::AbstractArray{<:Any, N}, interp_dims, cache) where {N} interp_dims = _add_trailing_interp_dims(interp_dims, Val{N}()) N_in = _count_interpolating_dims(interp_dims) N_out = _count_noninterpolating_dims(interp_dims) @@ -46,13 +46,19 @@ struct NDInterpolation{ end # TODO probably not type-stable (this needs to compile away completely) -_count_interpolating_dims(interp_dims) = count(map(d -> !(d isa NoInterpolationDimension), interp_dims)) -_count_noninterpolating_dims(interp_dims) = count(map(d -> d isa NoInterpolationDimension, interp_dims)) +function _count_interpolating_dims(interp_dims) + count(map(d -> !(d isa NoInterpolationDimension), interp_dims)) +end +function _count_noninterpolating_dims(interp_dims) + count(map(d -> d isa NoInterpolationDimension, interp_dims)) +end -_add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) = +function _add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) _add_trailing_interp_dims((dim,), n) -_add_trailing_interp_dims(dims::Tuple, ::Val{N}) where N = - (dims..., ntuple(_ -> NoInterpolationDimension(), Val{N-length(dims)}())...) +end +function _add_trailing_interp_dims(dims::Tuple, ::Val{N}) where {N} + (dims..., ntuple(_ -> NoInterpolationDimension(), Val{N - length(dims)}())...) +end # Constructor with optional global cache function NDInterpolation(u, interp_dims; cache = EmptyCache()) @@ -79,11 +85,11 @@ function (interp::NDInterpolation)( end # In place single input evaluation -function (interp::NDInterpolation{N,N_in,N_out})( +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} + derivative_orders::Tuple{Vararg{Integer, N_in}} = ntuple(_ -> 0, N_in) +) where {N, N_in, N_out} validate_size_u(interp, out) # We need to add the NoInterpolationDimensions for all arguments to _interpolate. # Currently interp() accepts only the interpolating dimensions. diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index f2b3e60..a187210 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -2,13 +2,14 @@ Base.@propagate_inbounds function _interpolate!( out::Union{Number, AbstractArray}, A::NDInterpolation{N}, ts::Tuple{Vararg{Any, N}}, - idx::Tuple{Vararg{Any ,N}}, + idx::Tuple{Vararg{Any, N}}, derivative_orders::Tuple{Vararg{Any, N}}, - multi_point_index, -) where N + multi_point_index +) where {N} (; interp_dims, cache, u) = A - out, valid_derivative_orders = check_derivative_order(interp_dims, derivative_orders, ts, out) + out, valid_derivative_orders = check_derivative_order( + interp_dims, derivative_orders, ts, out) valid_derivative_orders || return out if isnothing(multi_point_index) multi_point_index = map(_ -> nothing, interp_dims) @@ -17,7 +18,8 @@ Base.@propagate_inbounds function _interpolate!( # Setup out = make_zero!!(out) # TODO remove this stencils = map(stencil, interp_dims) - stencil_weights = map(weights, interp_dims, derivative_orders, multi_point_index, ts, idx) + stencil_weights = map( + weights, interp_dims, derivative_orders, multi_point_index, ts, idx) denom = zero(eltype(u)) # TODO this can be a single unrolled broadcast rather than a loop of .+= @@ -52,7 +54,7 @@ end function check_derivative_order(dims::Tuple, derivative_orders::Tuple, ts::Tuple, out) itr = map(tuple, dims, derivative_orders, ts) # Fold over itr for all dims, combining out and valid - foldl(itr; init=(out, true)) do (acc_out, acc_valid), (d, d_o, t) + foldl(itr; init = (out, true)) do (acc_out, acc_valid), (d, d_o, t) dim_out, dim_valid = check_derivative_order(d, d_o, t, acc_out) dim_out, dim_valid & acc_valid end @@ -75,7 +77,7 @@ end stencil(::LinearInterpolationDimension) = (1, 2) stencil(::ConstantInterpolationDimension) = 1 stencil(::NoInterpolationDimension) = 1 -stencil(d::BSplineInterpolationDimension) = 1:d.degree + 1 +stencil(d::BSplineInterpolationDimension) = 1:(d.degree + 1) # Precalculate weights function weights(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) @@ -88,8 +90,10 @@ function weights(d::LinearInterpolationDimension, derivative_order, multi_point_ end weights(::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 weights(::NoInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 -weights(d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) = +function weights( + d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) get_basis_function_values(d, t, i, derivative_order, multi_point_index) +end index(::LinearInterpolationDimension, t, idx, i) = idx + i - 1 # TODO: this should happen outside of the loop diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index 24323fb..d8445bb 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -37,10 +37,11 @@ function eval_unstructured!( derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) ) where {N_in} validate_derivative_orders(derivative_orders, interp; multi_point = true) - validate_output_size(out, interp) + validate_output_size(out, interp) backend = get_backend(out) # TODO this may be broken but it isn't tested - @assert all(d -> length(d.t_eval) == size(out, 1), remove(NoInterpolationDimension, interp.interp_dims)) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." + @assert all(d -> length(d.t_eval) == size(out, 1), + remove(NoInterpolationDimension, interp.interp_dims)) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." eval_kernel(backend)( out, interp, @@ -84,11 +85,11 @@ in place. """ function eval_grid!( out::AbstractArray, - interp::NDInterpolation{N,N_in}; + interp::NDInterpolation{N, N_in}; derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) -) where {N,N_in} +) where {N, N_in} validate_t_eval_lengths(out, interp) - validate_output_size(out, interp) + validate_output_size(out, interp) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) no_interp_inds = map(_ -> Colon(), remove(NoInterpolationDimension, interp.interp_dims)) @@ -105,14 +106,16 @@ end function validate_t_eval_lengths(out, interp) (; interp_dims) = interp - interp_sizes = removeat(NoInterpolationDimension, size(out), interp_dims) - interp_t_eval_lengths = map(d -> length(d.t_eval), remove(NoInterpolationDimension, interp_dims)) - all(map(==, interp_sizes, interp_t_eval_lengths)) || + interp_sizes = removeat(NoInterpolationDimension, size(out), interp_dims) + interp_t_eval_lengths = map( + d -> length(d.t_eval), remove(NoInterpolationDimension, interp_dims)) + all(map(==, interp_sizes, interp_t_eval_lengths)) || throw(ArgumentError("The length must match the t_eval of the corresponding interpolation dimension.")) end -function validate_output_size(out, interp) - keepat(NoInterpolationDimension, size(out), interp.interp_dims) == get_output_size(interp) || +function validate_output_size(out, interp) + keepat(NoInterpolationDimension, size(out), interp.interp_dims) == + get_output_size(interp) || throw(ArumentError("The size of the NoInterpolationDimension dimensions of `out` must be the same as the output size of the interpolation.")) end @@ -120,11 +123,11 @@ end out, A, # @Const(A), TODO: somehow this now hits a bug in KernelAbstractions where elsize is not defined for Const derivative_orders, - no_interp_inds, -) + no_interp_inds +) N_out = length(keep(NoInterpolationDimension, A.interp_dims)) I = @index(Global, NTuple) - k = insertat(NoInterpolationDimension, no_interp_inds, I, A.interp_dims) + k = insertat(NoInterpolationDimension, no_interp_inds, I, A.interp_dims) t_eval = map(get_t_eval, A.interp_dims, k) idx_eval = map(get_idx_eval, A.interp_dims, k) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 322fdf0..bd3f7e1 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -3,43 +3,45 @@ trivial_range(i::Integer) = i:i Base.length(itp_dim::AbstractInterpolationDimension) = length(itp_dim.t) function validate_derivative_order( - derivative_orders::NTuple, - A::NDInterpolation; - multi_point::Bool = false + derivative_orders::NTuple, + A::NDInterpolation; + multi_point::Bool = false ) map(derivative_orders, A.interp_dims) do d_o, d - validate_derivative_order(d_o, d; multi_point, cache=A.cache) + validate_derivative_order(d_o, d; multi_point, cache = A.cache) end end function validate_derivative_order( - derivative_order::Integer, - interp_dim::BSplineInterpolationDimension; - multi_point::Bool, - cache, + derivative_order::Integer, + interp_dim::BSplineInterpolationDimension; + multi_point::Bool, + cache ) if multi_point - @assert derivative_order ≤ interp_dim.max_derivative_order_eval """ - For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be - larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want - to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the - `BSplineInterpolationDimension` constructor(s). - """ + @assert derivative_order≤interp_dim.max_derivative_order_eval """ + For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be + larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want + to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the + `BSplineInterpolationDimension` constructor(s). + """ end validate_derivative_order_by_cache(cache, derivative_order) end function validate_derivative_order( - derivative_order::Integer, - interp_dim::AbstractInterpolationDimension; - multi_point::Bool, - cache, + derivative_order::Integer, + interp_dim::AbstractInterpolationDimension; + multi_point::Bool, + cache ) validate_derivative_order_by_cache(cache, derivative_order) end -validate_derivative_order_by_cache(::NURBSWeights, derivative_order) = - @assert derivative_order == 0 "Currently partial derivatives of NURBS are not supported." -validate_derivative_order_by_cache(::Any, derivative_order) = - @assert derivative_order >= 0 "Derivative orders must me non-negative." +function validate_derivative_order_by_cache(::NURBSWeights, derivative_order) + @assert derivative_order==0 "Currently partial derivatives of NURBS are not supported." +end +function validate_derivative_order_by_cache(::Any, derivative_order) + @assert derivative_order>=0 "Derivative orders must me non-negative." +end function validate_t(t) @assert t isa AbstractVector{<:Number} "t must be an AbstractVector with number like elements." @@ -48,43 +50,52 @@ end validate_size_u(interp::NDInterpolation, u) = validate_size_u(interp.interp_dims, u) validate_size_u(interp_dims::Tuple, u::Number) = nothing -validate_size_u(interp_dims::Tuple, u::AbstractArray) = validate_size_u(interp_dims, axes(u)) +function validate_size_u(interp_dims::Tuple, u::AbstractArray) + validate_size_u(interp_dims, axes(u)) +end validate_size_u(interp_dims::Tuple, ax::Tuple) = map(validate_size_u, interp_dims, ax) validate_size_u(interp_dim::NoInterpolationDimension, ax::AbstractRange) = nothing function validate_size_u(interp_dim::AbstractInterpolationDimension, ax::AbstractRange) - @assert length(interp_dim) == length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." + @assert length(interp_dim)==length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." end function validate_size_u( - interp_dim::BSplineInterpolationDimension, - ax::AbstractRange + interp_dim::BSplineInterpolationDimension, + ax::AbstractRange ) expected_size = get_n_basis_functions(interp_dim) - @assert expected_size == length(ax) "Expected the size to be $expected_size based on the BSplineInterpolation properties, got $(length(ax))." + @assert expected_size==length(ax) "Expected the size to be $expected_size based on the BSplineInterpolation properties, got $(length(ax))." end function validate_cache( - cache::AbstractInterpolationCache, dims::Tuple, u::AbstractArray + cache::AbstractInterpolationCache, dims::Tuple, u::AbstractArray ) ntuple(length(dims)) do n validate_cache(cache, dims[n], u, n) end end -validate_cache(::EmptyCache, ::AbstractInterpolationDimension, ::AbstractArray, ::Int) = nothing +function validate_cache( + ::EmptyCache, ::AbstractInterpolationDimension, ::AbstractArray, ::Int) + nothing +end validate_cache(::EmptyCache, ::NoInterpolationDimension, ::AbstractArray, ::Int) = nothing -validate_cache(::AbstractInterpolationCache, ::NoInterpolationDimension, ::AbstractArray, ::Int) = nothing +function validate_cache( + ::AbstractInterpolationCache, ::NoInterpolationDimension, ::AbstractArray, ::Int) + nothing +end function validate_cache( nurbs_weights::NURBSWeights, ::BSplineInterpolationDimension, u::AbstractArray, - n::Int, + n::Int ) size_expected = size(u, n) - @assert size(nurbs_weights.weights, n) == size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected), got $(size(nurbs_weights.weights, n))." + @assert size(nurbs_weights.weights, n)==size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected), got $(size(nurbs_weights.weights, n))." end function validate_cache( - ::gType, ::ID, ::AbstractArray, ::Int) where {gType,ID<:AbstractInterpolationDimension} + ::gType, ::ID, ::AbstractArray, ::Int) where { + gType, ID <: AbstractInterpolationDimension} @error("Interpolation dimension type $ID is not compatible with global cache type $gType.") end @@ -113,13 +124,13 @@ function make_zero!!(v::AbstractArray) end function make_out( - interp::NDInterpolation{<:Any,N_in, 0}, + interp::NDInterpolation{<:Any, N_in, 0}, t::NTuple{N_in, >:Number} ) where {N_in} zero(eltype(interp.u)) end function make_out( - interp::NDInterpolation{<:Any,N_in}, + interp::NDInterpolation{<:Any, N_in}, t::NTuple{N_in, >:Number} ) where {N_in} T = promote_type(eltype(interp.u), map(eltype, t)...) @@ -158,9 +169,10 @@ function get_idx( end end # t_eval must already be an index for NoInterpolationDimension -get_idx(::NoInterpolationDimension, t_eval::Union{Colon,Int,AbstractArray{Int}}) = t_eval -get_idx(interp_dims::Tuple{Vararg{Any,N}}, t::Tuple{Vararg{Any,N}}) where N = +get_idx(::NoInterpolationDimension, t_eval::Union{Colon, Int, AbstractArray{Int}}) = t_eval +function get_idx(interp_dims::Tuple{Vararg{Any, N}}, t::Tuple{Vararg{Any, N}}) where {N} map(get_idx, interp_dims, t) +end function set_eval_idx!( interp_dim::AbstractInterpolationDimension, @@ -193,7 +205,6 @@ end typed_nan(::T) where {T <: Integer} = zero(T) typed_nan(::T) where {T <: AbstractFloat} = T(NaN) - # Get the KernelAbstractions nd_range, over interpolated dimensions function get_ndrange(interp::NDInterpolation) I = map(interp.interp_dims) do d @@ -202,67 +213,76 @@ function get_ndrange(interp::NDInterpolation) return remove(Nothing, I) end - # Some tuple handling primitives # TODO: make sure these compile away completely # Remove objects of type `T` from `in` (reworked from DimensionalData.jl) -remove(::Type{T}, in) where T = _remove(T, in...) -_remove(::Type{T}, x, xs...) where T = (x, _remove(T, xs...)...) -_remove(::Type{T}, x::T, xs...) where T = _remove(T, xs...) +remove(::Type{T}, in) where {T} = _remove(T, in...) +_remove(::Type{T}, x, xs...) where {T} = (x, _remove(T, xs...)...) +_remove(::Type{T}, x::T, xs...) where {T} = _remove(T, xs...) _remove(::Type) = () # Keep only objects of type `T` from `in` -keep(::Type{T}, in) where T = _keep(T, in...) -_keep(::Type{T}, x, xs...) where T = _keep(T, xs...) -_keep(::Type{T}, x::T, xs...) where T = (x, _keep(T, xs...)...) +keep(::Type{T}, in) where {T} = _keep(T, in...) +_keep(::Type{T}, x, xs...) where {T} = _keep(T, xs...) +_keep(::Type{T}, x::T, xs...) where {T} = (x, _keep(T, xs...)...) _keep(::Type) = () # Remove values from `in` where `matches` are of type `T` -function removeat(::Type{T}, in::Tuple, matches::Tuple) where T - @assert length(in) == length(matches) - _removeat(T, in, matches...) +function removeat(::Type{T}, in::Tuple, matches::Tuple) where {T} + @assert length(in) == length(matches) + _removeat(T, in, matches...) end # If `!(m isa T)` take from `in` -_removeat(::Type{T}, in::Tuple{<:Any,Vararg}, m, ms...) where T = +function _removeat(::Type{T}, in::Tuple{<:Any, Vararg}, m, ms...) where {T} (first(in), _removeat(T, Base.tail(in), ms...)...) +end # If `m isa T` remove -_removeat(::Type{T}, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = +function _removeat(::Type{T}, in::Tuple{<:Any, Vararg}, m::T, ms...) where {T} _removeat(T, Base.tail(in), ms...) +end # `in` can be empty if there are no remaining `m` _removeat(::Type, in::Tuple{}) = () # Keep only values from `in` where `matches` are of type `T` -function keepat(::Type{T}, in::Tuple, matches::Tuple) where T - @assert length(in) == length(matches) - _keepat(T, in, matches...) +function keepat(::Type{T}, in::Tuple, matches::Tuple) where {T} + @assert length(in) == length(matches) + _keepat(T, in, matches...) end # If `!(m isa T)` disgaurd -_keepat(::Type{T}, in::Tuple{<:Any,Vararg}, m, ms...) where T = +function _keepat(::Type{T}, in::Tuple{<:Any, Vararg}, m, ms...) where {T} _keepat(T, Base.tail(in), ms...) +end # If `m isa T` keep -_keepat(::Type{T}, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = +function _keepat(::Type{T}, in::Tuple{<:Any, Vararg}, m::T, ms...) where {T} (first(in), _keepat(T, Base.tail(in), ms...)...) +end # `in` can be empty if there are no remaining `m` _keepat(::Type, in::Tuple{}) = () # Insert x into `in` where `matches` are of type `T`, # otherwise output one of `in` for each `m`. -insertat(::Type{T}, x, in::Tuple, matches::Tuple) where T = - _insertat(T, x, in, matches...) +function insertat(::Type{T}, x, in::Tuple, matches::Tuple) where {T} + _insertat(T, x, in, matches...) +end # If `!(m isa T)` take from `in` -_insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, m, ms...) where T = +function _insertat(::Type{T}, x, in::Tuple{<:Any, Vararg}, m, ms...) where {T} (first(in), _insertat(T, x, Base.tail(in), ms...)...) +end # If `m isa T` insert `x` -_insertat(::Type{T}, x, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = +function _insertat(::Type{T}, x, in::Tuple{<:Any, Vararg}, m::T, ms...) where {T} (x, _insertat(T, x, in, ms...)...) +end # For Tuple x we insert the first -_insertat(::Type{T}, xs::Tuple, in::Tuple{<:Any,Vararg}, m::T, ms...) where T = +function _insertat(::Type{T}, xs::Tuple, in::Tuple{<:Any, Vararg}, m::T, ms...) where {T} (first(xs), _insertat(T, Base.tail(xs), in, ms...)...) +end # `in` can be empty if there are no remaining `m` _insertat(::Type, x, in::Tuple{}) = () # `in` can also be empty if all trailing ms are T -_insertat(::Type{T}, x, in::Tuple{}, m::T, ms::T...) where T = +function _insertat(::Type{T}, x, in::Tuple{}, m::T, ms::T...) where {T} (x, _insertat(T, x, in, ms...)...) -_insertat(::Type{T}, xs::Tuple, in::Tuple{}, m::T, ms::T...) where T = +end +function _insertat(::Type{T}, xs::Tuple, in::Tuple{}, m::T, ms::T...) where {T} (first(xs), _insertat(T, Base.tail(xs), in, ms...)...) +end diff --git a/src/spline_utils.jl b/src/spline_utils.jl index be78a8f..595876e 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -60,7 +60,7 @@ function get_basis_function_values( t::Number, idx::Integer, derivative_order::Integer, - multi_point_index::Nothing, + multi_point_index::Nothing ) (; degree, knots_all) = itp_dim T = promote_type(typeof(t), eltype(itp_dim.basis_function_eval)) @@ -99,7 +99,7 @@ function get_basis_function_values( t::Number, idx::Integer, derivative_order::Integer, - multi_point_index::Integer, + multi_point_index::Integer ) view(itp_dim.basis_function_eval, multi_point_index, :, derivative_order + 1) @@ -113,7 +113,8 @@ function get_basis_function_values_all( derivative_orders::Tuple, multi_point_index ) - map(get_basis_function_values, A.interp_dims, ts, idx, derivative_orders, multi_point_index) + map(get_basis_function_values, A.interp_dims, ts, + idx, derivative_orders, multi_point_index) end function set_basis_function_eval!(itp_dim::BSplineInterpolationDimension)::Nothing @@ -136,7 +137,7 @@ end itp_dim.t_eval[i], itp_dim.idx_eval[i], derivative_order_plus_1 - 1, - nothing, + nothing ) end diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 36320ca..3d0ee58 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -1,27 +1,27 @@ using DataInterpolationsND using Random -using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache, - keepat, removeat, insertat, remove, keep +using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache, + keepat, removeat, insertat, remove, keep @testset "tuple primitives" begin # These help with handling NoInterpolationDimension and its consequences @test insertat(Nothing, 1, (2,), (nothing, 4, nothing)) === (1, 2, 1) @test insertat(Float64, missing, (1,), (8, 9.0)) === (1, missing) @test keepat(Nothing, (1, 2, 3), (nothing, 4, nothing)) === (1, 3) - @test keepat(Float64, (1, 2,), (8, 9.0)) === (2,) + @test keepat(Float64, (1, 2), (8, 9.0)) === (2,) @test removeat(Nothing, (1, 2, 3), (nothing, 4, nothing)) === (2,) - @test removeat(Float64, (1, 2,), (8, 9.0)) === (1,) + @test removeat(Float64, (1, 2), (8, 9.0)) === (1,) @test remove(Nothing, (1, nothing, 3)) === (1, 3) @test keep(Int, (1, nothing, 3)) === (1, 3) end -function test_globally_constant(ID::Type{<:AbstractInterpolationDimension}; - args1 = (), - args2 = (), - kwargs1 = (), - kwargs2 = (), - cache = EmptyCache(), - test_derivatives = true +function test_globally_constant(ID::Type{<:AbstractInterpolationDimension}; + args1 = (), + args2 = (), + kwargs1 = (), + kwargs2 = (), + cache = EmptyCache(), + test_derivatives = true ) t1 = [-3.14, 1.0, 3.0, 7.6, 12.8] t2 = [-2.71, 1.41, 12.76, 50.2, 120.0] @@ -64,7 +64,7 @@ function test_globally_constant(ID::Type{<:AbstractInterpolationDimension}; end end -function test_analytic(itp::NDInterpolation{<:Any,N_in}, f) where {N_in} +function test_analytic(itp::NDInterpolation{<:Any, N_in}, f) where {N_in} used_interp_dims = remove(NoInterpolationDimension, itp.interp_dims) # Evaluation in data points ts = map(d -> d.t, used_interp_dims) @@ -96,7 +96,6 @@ end itp = NDInterpolation(u, itp_dims) test_analytic(itp, f) - itp_dims = ( NoInterpolationDimension(), LinearInterpolationDimension(t2; t_eval = t2 ./ 2) @@ -111,7 +110,7 @@ end args2 = (3,) kwargs1 = (:max_derivative_order_eval => 1,) kwargs2 = (:max_derivative_order_eval => 1,) - cache = EmptyCache() + cache = EmptyCache() test_derivatives = true test_globally_constant(ID; args1, args2, kwargs1, kwargs2, cache, test_derivatives) @@ -180,8 +179,8 @@ end u = ones(3, 4, 5) itp_dims = ( NoInterpolationDimension(), - LinearInterpolationDimension(t2; t_eval=t2), - BSplineInterpolationDimension(t3, 1; t_eval=t3), + LinearInterpolationDimension(t2; t_eval = t2), + BSplineInterpolationDimension(t3, 1; t_eval = t3) ) itp = NDInterpolation(u, itp_dims) @test isapprox(itp(1.0, 4.0), [1.0, 1.0, 1.0]) From 05f2c4c44cfdf68af80dab1eec7fcc7fa9619ab8 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Fri, 31 Oct 2025 16:46:58 +1100 Subject: [PATCH 21/32] use JuliaFormatter v2 --- src/interpolation_methods.jl | 3 ++- src/interpolation_utils.jl | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index a187210..1aac479 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -8,7 +8,8 @@ Base.@propagate_inbounds function _interpolate!( ) where {N} (; interp_dims, cache, u) = A - out, valid_derivative_orders = check_derivative_order( + out, + valid_derivative_orders = check_derivative_order( interp_dims, derivative_orders, ts, out) valid_derivative_orders || return out if isnothing(multi_point_index) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index bd3f7e1..b1854e5 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -19,11 +19,11 @@ function validate_derivative_order( ) if multi_point @assert derivative_order≤interp_dim.max_derivative_order_eval """ - For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be - larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want - to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the - `BSplineInterpolationDimension` constructor(s). - """ +For BSpline interpolation, when using multi-point evaluation the derivative orders cannot be +larger than the `max_derivative_order_eval` eval of of the `BSplineInterpolationDimension`. If you want +to compute higher order multi-point derivatives, pass a larger `max_derivative_order_eval` to the +`BSplineInterpolationDimension` constructor(s). +""" end validate_derivative_order_by_cache(cache, derivative_order) end From a27ae40f783ed11288e3a85cd41ede23ee431567 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 22:17:09 +1100 Subject: [PATCH 22/32] bugfix eval_unstructured --- src/DataInterpolationsND.jl | 9 +++--- src/interpolation_parallel.jl | 59 ++++++++++++++++++++++------------- src/interpolation_utils.jl | 2 +- 3 files changed, 44 insertions(+), 26 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index b3a2ffb..3f98b4f 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -90,14 +90,15 @@ function (interp::NDInterpolation{N, N_in, 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} - validate_size_u(interp, 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. - t_all = insertat(NoInterpolationDimension, Colon(), t, interp.interp_dims) - idx_all = get_idx(interp.interp_dims, t_all) - d_o_all = insertat(NoInterpolationDimension, 0, derivative_orders, interp.interp_dims) + 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 diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index d8445bb..ea2c931 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -33,19 +33,22 @@ length for each interpolation dimension and the interpolation is evaluated at th """ function eval_unstructured!( out::AbstractArray, - interp::NDInterpolation{N_in}; + interp::NDInterpolation{N,N_in}; derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) -) where {N_in} - validate_derivative_orders(derivative_orders, interp; multi_point = true) - validate_output_size(out, interp) +) where {N,N_in} + validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) - # TODO this may be broken but it isn't tested + no_interp_inds = map(_ -> Colon(), keep(NoInterpolationDimension, interp.interp_dims)) @assert all(d -> length(d.t_eval) == size(out, 1), remove(NoInterpolationDimension, interp.interp_dims)) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." + @assert size(out)[2:end]==get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." + eval_grid = false eval_kernel(backend)( out, interp, derivative_orders, + no_interp_inds, + eval_grid, ndrange = size(out, 1) ) synchronize(backend) @@ -85,19 +88,21 @@ in place. """ function eval_grid!( out::AbstractArray, - interp::NDInterpolation{N, N_in}; - derivative_orders::NTuple{N, <:Integer} = ntuple(_ -> 0, N) -) where {N, N_in} + interp::NDInterpolation{N, N_in, N_out}; + derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in), + no_interp_inds = map(_ -> Colon(), N_out) +) where {N, N_in, N_out} validate_t_eval_lengths(out, interp) validate_output_size(out, interp) validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) - no_interp_inds = map(_ -> Colon(), remove(NoInterpolationDimension, interp.interp_dims)) + eval_grid = true eval_kernel(backend)( out, interp, derivative_orders, no_interp_inds, + eval_grid, ndrange = get_ndrange(interp) ) synchronize(backend) @@ -123,35 +128,47 @@ end out, A, # @Const(A), TODO: somehow this now hits a bug in KernelAbstractions where elsize is not defined for Const derivative_orders, - no_interp_inds + no_interp_inds, + eval_grid, ) - N_out = length(keep(NoInterpolationDimension, A.interp_dims)) + (; interp_dims) = A + N_out = length(keep(NoInterpolationDimension, interp_dims)) I = @index(Global, NTuple) - k = insertat(NoInterpolationDimension, no_interp_inds, I, A.interp_dims) - - t_eval = map(get_t_eval, A.interp_dims, k) - idx_eval = map(get_idx_eval, A.interp_dims, k) + d_o = insertat(NoInterpolationDimension, 0, derivative_orders, interp_dims) + + # TODO eval_grid should be a static bool or Val so this branch can be deleted + if eval_grid + # Insert no_interp_inds in I + k = insertat(NoInterpolationDimension, no_interp_inds, I, interp_dims) + else # eval_unstructured + # The same index is used for all interp dimensions in eval_unstructured + I1 = map(_ -> only(I), remove(NoInterpolationDimension, interp_dims)) + # Insert no_interp_inds in I1 + k = insertat(NoInterpolationDimension, no_interp_inds, I1, interp_dims) + end + t_eval = map(get_t_eval, interp_dims, k) + idx_eval = map(get_idx_eval, interp_dims, k) if iszero(N_out) @inbounds out[k...] = _interpolate!( - make_out(A, t_eval), A, t_eval, idx_eval, derivative_orders, k) + make_out(A, t_eval), A, t_eval, idx_eval, d_o, k) else _interpolate!( view(out, k...), - A, t_eval, idx_eval, derivative_orders, k) + A, t_eval, idx_eval, d_o, k) end end -get_t_eval(d::AbstractInterpolationDimension, i) = d.t_eval[i] +get_t_eval(d::AbstractInterpolationDimension, i::Integer) = d.t_eval[i] get_t_eval(d::AbstractInterpolationDimension, i::Colon) = d.t_eval get_t_eval(d::AbstractInterpolationDimension, i::AbstractVector) = d.t_eval[i] -get_t_eval(::NoInterpolationDimension, i) = nothing +get_t_eval(::NoInterpolationDimension, i::Integer) = nothing get_t_eval(::NoInterpolationDimension, i::Colon) = nothing get_t_eval(::NoInterpolationDimension, i::AbstractVector) = nothing -get_idx_eval(d::AbstractInterpolationDimension, i) = d.idx_eval[i] +get_idx_eval(d::AbstractInterpolationDimension, i::Integer) = d.idx_eval[i] get_idx_eval(d::AbstractInterpolationDimension, i::AbstractVector) = d.idx_eval[i] get_idx_eval(d::AbstractInterpolationDimension, i::Colon) = d.idx_eval -get_idx_eval(::NoInterpolationDimension, i) = i +get_idx_eval(::NoInterpolationDimension, i::Integer) = i get_idx_eval(::NoInterpolationDimension, i::Colon) = i get_idx_eval(::NoInterpolationDimension, i::AbstractVector) = i diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index b1854e5..f46570e 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -56,7 +56,7 @@ end validate_size_u(interp_dims::Tuple, ax::Tuple) = map(validate_size_u, interp_dims, ax) validate_size_u(interp_dim::NoInterpolationDimension, ax::AbstractRange) = nothing function validate_size_u(interp_dim::AbstractInterpolationDimension, ax::AbstractRange) - @assert length(interp_dim)==length(ax) "For the first N_in dimensions of u the length must match the t of the corresponding interpolation dimension." + @assert length(interp_dim) == length(ax) "The size if `u` must match the t of the corresponding interpolation dimension. Got $interp_dim and $ax" end function validate_size_u( From 0222ba041e1bb16fcf66db602abafcc5c55bc9d2 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 22:23:00 +1100 Subject: [PATCH 23/32] review suggestions --- src/interpolation_methods.jl | 15 +++++++-------- src/interpolation_utils.jl | 2 +- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 1aac479..96fa2db 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -19,14 +19,13 @@ Base.@propagate_inbounds function _interpolate!( # Setup out = make_zero!!(out) # TODO remove this stencils = map(stencil, interp_dims) - stencil_weights = map( - weights, interp_dims, derivative_orders, multi_point_index, ts, idx) + coeffs = map(coefficients, interp_dims, derivative_orders, multi_point_index, ts, idx) denom = zero(eltype(u)) # TODO this can be a single unrolled broadcast rather than a loop of .+= for I in Iterators.product(stencils...) J = map(index, interp_dims, ts, idx, I) - product = prod(map(getindex, stencil_weights, I)) + product = prod(map(getindex, coeffs, I)) if !(cache isa EmptyCache) K = removeat(NoInterpolationDimension, J, interp_dims) @@ -80,8 +79,8 @@ stencil(::ConstantInterpolationDimension) = 1 stencil(::NoInterpolationDimension) = 1 stencil(d::BSplineInterpolationDimension) = 1:(d.degree + 1) -# Precalculate weights -function weights(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) +# Precalculate coefficient/s +function coefficients(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) t₁ = d.t[i] t₂ = d.t[i + 1] t_vol_inv = inv(t₂ - t₁) @@ -89,9 +88,9 @@ function weights(d::LinearInterpolationDimension, derivative_order, multi_point_ b = (iszero(derivative_order) ? t - t₁ : one(t)) * t_vol_inv return (a, b) end -weights(::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 -weights(::NoInterpolationDimension, derivative_order, multi_point_index, t, i) = 1 -function weights( +coefficients(::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) = true +coefficients(::NoInterpolationDimension, derivative_order, multi_point_index, t, i) = true +function coefficients( d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) get_basis_function_values(d, t, i, derivative_order, multi_point_index) end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index f46570e..98fbd2a 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -40,7 +40,7 @@ function validate_derivative_order_by_cache(::NURBSWeights, derivative_order) @assert derivative_order==0 "Currently partial derivatives of NURBS are not supported." end function validate_derivative_order_by_cache(::Any, derivative_order) - @assert derivative_order>=0 "Derivative orders must me non-negative." + @assert derivative_order>=0 "Derivative orders must be non-negative." end function validate_t(t) From 00704e49d62932d06b44a342f21feb6eec1a2ac9 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 22:31:11 +1100 Subject: [PATCH 24/32] use NURBS --- src/interpolation_methods.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 96fa2db..087de24 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -27,7 +27,7 @@ Base.@propagate_inbounds function _interpolate!( J = map(index, interp_dims, ts, idx, I) product = prod(map(getindex, coeffs, I)) - if !(cache isa EmptyCache) + if cache isa NURBSWeights K = removeat(NoInterpolationDimension, J, interp_dims) product *= cache.weights[K...] denom += product @@ -40,7 +40,7 @@ Base.@propagate_inbounds function _interpolate!( end end - if !(cache isa EmptyCache) + if cache isa NURBSWeights if out isa AbstractArray out ./= denom else From b37af8cee1b051616e2be5228892a1d9ac4cfeac Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 22:40:36 +1100 Subject: [PATCH 25/32] dim checks --- src/DataInterpolationsND.jl | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 3f98b4f..4573e43 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -33,10 +33,10 @@ struct NDInterpolation{ interp_dims::D cache::gType function NDInterpolation(u::AbstractArray{<:Any, N}, interp_dims, cache) where {N} + @assert length(dims) <= N "The number of interpolation dimensions must no more than the dimensions of u." interp_dims = _add_trailing_interp_dims(interp_dims, Val{N}()) - N_in = _count_interpolating_dims(interp_dims) - N_out = _count_noninterpolating_dims(interp_dims) - @assert N_out≥0 "The number of dimensions of u must be at least the number of interpolation dimensions." + N_in = count(remove(nointerpolationdimension, interp_dims)) + N_out = count(keep(nointerpolationdimension, interp_dims)) validate_size_u(interp_dims, u) validate_cache(cache, interp_dims, u) new{N, N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( @@ -45,13 +45,6 @@ struct NDInterpolation{ end end -# TODO probably not type-stable (this needs to compile away completely) -function _count_interpolating_dims(interp_dims) - count(map(d -> !(d isa NoInterpolationDimension), interp_dims)) -end -function _count_noninterpolating_dims(interp_dims) - count(map(d -> d isa NoInterpolationDimension, interp_dims)) -end function _add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) _add_trailing_interp_dims((dim,), n) From 45fbe9baf7f45b241bd2bb4ca574bf8b4c262f70 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 22:57:11 +1100 Subject: [PATCH 26/32] bugfix constructor changes --- src/DataInterpolationsND.jl | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 4573e43..444769b 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -33,10 +33,13 @@ struct NDInterpolation{ interp_dims::D cache::gType function NDInterpolation(u::AbstractArray{<:Any, N}, interp_dims, cache) where {N} - @assert length(dims) <= N "The number of interpolation dimensions must no more than the dimensions of u." - interp_dims = _add_trailing_interp_dims(interp_dims, Val{N}()) - N_in = count(remove(nointerpolationdimension, interp_dims)) - N_out = count(keep(nointerpolationdimension, interp_dims)) + if interp_dims isa AbstractInterpolationDimension + interp_dims = (interp_dims,) + end + @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) new{N, N_in, N_out, typeof(cache), typeof(interp_dims), typeof(u)}( @@ -44,20 +47,18 @@ struct NDInterpolation{ ) end end +# Constructor with optional global cache +function NDInterpolation(u, interp_dims; cache = EmptyCache()) + NDInterpolation(u, interp_dims, cache) +end - -function _add_trailing_interp_dims(dim::AbstractInterpolationDimension, n) +function _add_trailing_no_interp_dims(dim::AbstractInterpolationDimension, n) _add_trailing_interp_dims((dim,), n) end -function _add_trailing_interp_dims(dims::Tuple, ::Val{N}) where {N} +function _add_trailing_no_interp_dims(dims::Tuple, ::Val{N}) where {N} (dims..., ntuple(_ -> NoInterpolationDimension(), Val{N - length(dims)}())...) end -# Constructor with optional global cache -function NDInterpolation(u, interp_dims; cache = EmptyCache()) - NDInterpolation(u, interp_dims, cache) -end - @adapt_structure NDInterpolation include("interpolation_dimensions.jl") From 3537e53dc1cb0c0430f7f0b5b3ad0d4d2107cf48 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 23:02:12 +1100 Subject: [PATCH 27/32] comments and format --- src/interpolation_methods.jl | 10 +++++++--- src/interpolation_parallel.jl | 6 +++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 087de24..c004638 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -11,7 +11,7 @@ Base.@propagate_inbounds function _interpolate!( out, valid_derivative_orders = check_derivative_order( interp_dims, derivative_orders, ts, out) - valid_derivative_orders || return out + valid_derivative_orders || return out # Array was zeroed out in this case if isnothing(multi_point_index) multi_point_index = map(_ -> nothing, interp_dims) end @@ -80,7 +80,8 @@ stencil(::NoInterpolationDimension) = 1 stencil(d::BSplineInterpolationDimension) = 1:(d.degree + 1) # Precalculate coefficient/s -function coefficients(d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) +function coefficients( + d::LinearInterpolationDimension, derivative_order, multi_point_index, t, i) t₁ = d.t[i] t₂ = d.t[i + 1] t_vol_inv = inv(t₂ - t₁) @@ -88,7 +89,10 @@ function coefficients(d::LinearInterpolationDimension, derivative_order, multi_p b = (iszero(derivative_order) ? t - t₁ : one(t)) * t_vol_inv return (a, b) end -coefficients(::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) = true +function coefficients( + ::ConstantInterpolationDimension, derivative_order, multi_point_index, t, i) + true +end coefficients(::NoInterpolationDimension, derivative_order, multi_point_index, t, i) = true function coefficients( d::BSplineInterpolationDimension, derivative_order, multi_point_index, t, i) diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index ea2c931..9f79e55 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -33,9 +33,9 @@ length for each interpolation dimension and the interpolation is evaluated at th """ function eval_unstructured!( out::AbstractArray, - interp::NDInterpolation{N,N_in}; + interp::NDInterpolation{N, N_in}; derivative_orders::NTuple{N_in, <:Integer} = ntuple(_ -> 0, N_in) -) where {N,N_in} +) where {N, N_in} validate_derivative_order(derivative_orders, interp; multi_point = true) backend = get_backend(out) no_interp_inds = map(_ -> Colon(), keep(NoInterpolationDimension, interp.interp_dims)) @@ -129,7 +129,7 @@ end A, # @Const(A), TODO: somehow this now hits a bug in KernelAbstractions where elsize is not defined for Const derivative_orders, no_interp_inds, - eval_grid, + eval_grid ) (; interp_dims) = A N_out = length(keep(NoInterpolationDimension, interp_dims)) From 79acff43e66e488337c14bed1bb65e5cf1c12d15 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Tue, 4 Nov 2025 23:08:19 +1100 Subject: [PATCH 28/32] fix NoInterpolationDimension docstring --- src/interpolation_dimensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation_dimensions.jl b/src/interpolation_dimensions.jl index 5c4107a..fcf2b87 100644 --- a/src/interpolation_dimensions.jl +++ b/src/interpolation_dimensions.jl @@ -1,5 +1,5 @@ """ - NoInterpolationDimensio + NoInterpolationDimension() A dimension that does not perform interpolation. """ From 9d8774424f50553124fb852ce07252d865a745d9 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 15 Feb 2026 16:37:48 +0100 Subject: [PATCH 29/32] Fixes --- Project.toml | 2 +- src/spline_utils.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 8374f58..450cc13 100644 --- a/Project.toml +++ b/Project.toml @@ -44,4 +44,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua, ""DataInterpolations", "ForwardDiff", "Pkg", "Random", "SafeTestsets", "Symbolics", "Test", "SymbolicUtils"] +test = ["Aqua", "DataInterpolations", "ForwardDiff", "Pkg", "Random", "SafeTestsets", "Symbolics", "Test", "SymbolicUtils"] diff --git a/src/spline_utils.jl b/src/spline_utils.jl index cdb007a..ba0f98f 100644 --- a/src/spline_utils.jl +++ b/src/spline_utils.jl @@ -67,7 +67,7 @@ function get_basis_function_values( idx::Integer, derivative_order::Integer, multi_point_index::Nothing - ) + ) where {T_t} (; degree, knots_all) = itp_dim T_eval = eltype(itp_dim.basis_function_eval) # Use a concrete zero value for type stability From f7e25bd1fa6769a5846b06a95d54e5038d389396 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 15 Feb 2026 19:33:02 +0100 Subject: [PATCH 30/32] Mix interpolation types in examples --- README.md | 6 +++--- docs/src/index.md | 4 ++-- docs/src/usage.md | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index b8f4b7e..638d816 100644 --- 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 index 98d93a0..ab3b014 100644 --- 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 index 47c6b13..428d627 100644 --- 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) From 1d07fce89b3cbc4d2cecf1b62280073f46857804 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 15 Feb 2026 19:38:39 +0100 Subject: [PATCH 31/32] Update NDInterpolation docstring --- src/DataInterpolationsND.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/DataInterpolationsND.jl b/src/DataInterpolationsND.jl index 15cbebe..49cb7f0 100644 --- a/src/DataInterpolationsND.jl +++ b/src/DataInterpolationsND.jl @@ -13,13 +13,13 @@ 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 From a10aabf3c06d733d76cc647eca188193c03c2c66 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sun, 15 Feb 2026 19:57:14 +0100 Subject: [PATCH 32/32] Remove redundant comment --- src/interpolation_parallel.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/interpolation_parallel.jl b/src/interpolation_parallel.jl index a1a27ad..b58726a 100644 --- a/src/interpolation_parallel.jl +++ b/src/interpolation_parallel.jl @@ -43,8 +43,6 @@ function eval_unstructured!( d -> length(d.t_eval) == size(out, 1), remove(NoInterpolationDimension, interp.interp_dims) ) "The t_eval of all interpolation dimensions must have the same length as the first dimension of out." - # TODO: This also must be generalized: The shape of out must correspond to the size of the NonInterpolationDimension dimensions, - # assuming 'not indexing them' @assert size(out)[2:end] == get_output_size(interp) "The size of the last N_out dimensions of out must be the same as the output size of the interpolation." eval_grid = false eval_kernel(backend)(