Skip to content

Commit

Permalink
Make sure spline evaluation is type stable (#36)
Browse files Browse the repository at this point in the history
This fixes type stability of spline evaluation when working with
Float32 (see new tests).
  • Loading branch information
jipolanco authored Apr 11, 2022
1 parent e84ba35 commit 3af291c
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 2 deletions.
5 changes: 3 additions & 2 deletions src/Splines/spline.jl
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,7 @@ function spline_kernel(
) where {T,k}
# Algorithm adapted from https://en.wikipedia.org/wiki/De_Boor's_algorithm
if @generated
d_k = Symbol(:d_, k)
quote
w_0 = zero(T) # this is to make the compiler happy with w_{j - 1}
@nexprs $k j -> d_j = @inbounds c[j + n - $k]
Expand All @@ -159,13 +160,13 @@ function spline_kernel(
j -> d_j = if j r
α = @inbounds (x - t[j + n - k]) /
(t[j + n - r + 1] - t[j + n - k])
(1 - α) * w_{j - 1} + α * w_{j}
T((1 - α) * w_{j - 1} + α * w_{j})
else
w_j
end
)
end
@nexprs 1 j -> d_{$k} # return d_k
return $d_k
end
else
# Similar using tuples (slower than @generated version).
Expand Down
9 changes: 9 additions & 0 deletions test/splines.jl
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,15 @@ function test_splines(::BSplineOrder{k}) where {k}
@inferred (() -> BSplineBasis(k, copy(x)))()
g = BSplineBasis(k, copy(x))

@testset "Type stability" begin
# The return type of evaluating a spline should be the element type
# of the coefficient vector.
coefs = randn(Float32, length(g))
S = @inferred Spline(g, coefs)
xeval = Float64(0.32)
@test @inferred(S(xeval)) isa Float32
end

@testset "BSplineBasis equality" begin
h = BSplineBasis(k, copy(x))
@test g == h
Expand Down

0 comments on commit 3af291c

Please sign in to comment.