From a5fbe1d55d127a28ceaa51c2bbe9724e5ed80820 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <juan-ignacio.polanco@cnrs.fr>
Date: Mon, 19 Feb 2024 14:42:06 +0100
Subject: [PATCH] Fix evaluation of some recombined bases with Float32 (#84)

* Fix evaluation using natural BCs with Float32

* Add tests

* v0.17.1

* Fix ambiguity
---
 Project.toml                   |  2 +-
 src/Recombinations/matrices.jl | 20 ++++++++++++--------
 test/natural.jl                | 27 ++++++++++++++++++++++-----
 3 files changed, 35 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 5b9eeaf1..d5a86785 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BSplineKit"
 uuid = "093aae92-e908-43d7-9660-e50ee39d5a0a"
 authors = ["Juan Ignacio Polanco <juan-ignacio.polanco@cnrs.fr>"]
-version = "0.17.0"
+version = "0.17.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
diff --git a/src/Recombinations/matrices.jl b/src/Recombinations/matrices.jl
index 7858d538..d44e4da4 100644
--- a/src/Recombinations/matrices.jl
+++ b/src/Recombinations/matrices.jl
@@ -148,14 +148,17 @@ end
 
 # Default element type of recombination matrix.
 # In some specific cases we can use Bool...
-_default_eltype(::BoundaryCondition) = Float64
-_default_eltype(::Derivative{0}) = Bool  # Dirichlet BCs
-_default_eltype(::Derivative{1}) = Bool  # Neumann BCs
-_default_eltype(::Vararg{AbstractDifferentialOp}) = Float64
+_default_eltype(::Type{T}, ::BoundaryCondition) where {T <: AbstractFloat} = T
+_default_eltype(::Type{T}, ::Derivative{0}) where {T <: AbstractFloat} = Bool  # Dirichlet BCs
+_default_eltype(::Type{T}, ::Derivative{1}) where {T <: AbstractFloat} = Bool  # Neumann BCs
+_default_eltype(::Type{T}, ::Vararg{AbstractDifferentialOp}) where {T <: AbstractFloat} = T
 
 # Case (D(0), D(1), D(2), ...)
-_default_eltype(::Derivative{0}, ::Derivative{1}, ::Vararg{Derivative}) = Bool  # TODO this isn't always right, is it?
-_default_eltype(ops::DiffOpList) = _default_eltype(ops...)
+_default_eltype(
+    ::Type{T}, ::Derivative{0}, ::Derivative{1}, ::Vararg{Derivative},
+) where {T <: AbstractFloat} = Bool  # TODO this isn't always right, is it?
+_default_eltype(::Type{T}, ops::DiffOpList) where {T <: AbstractFloat} =
+    _default_eltype(T, ops...)
 
 # Same BCs on both sides.
 RecombineMatrix(B::BSplineBasis, ops) = RecombineMatrix(B, ops, ops)
@@ -166,8 +169,9 @@ RecombineMatrix(ops, B::BSplineBasis, args...) = RecombineMatrix(B, ops, args...
 
 # No element type specified.
 function RecombineMatrix(B::BSplineBasis, ops_l, ops_r)
-    Tl = _default_eltype(ops_l)
-    Tr = _default_eltype(ops_r)
+    S = eltype(knots(B))
+    Tl = _default_eltype(S, ops_l)
+    Tr = _default_eltype(S, ops_r)
     T = promote_type(Tl, Tr)
     RecombineMatrix(B, ops_l, ops_r, T)
 end
diff --git a/test/natural.jl b/test/natural.jl
index 3480777a..7fa3c7b0 100644
--- a/test/natural.jl
+++ b/test/natural.jl
@@ -7,9 +7,12 @@ using Test
 
 import LinearAlgebra
 
-function test_natural(ord::BSplineOrder)
-    B = BSplineBasis(ord, -1:0.1:1)
+function test_natural(::Type{T}, ord::BSplineOrder) where {T}
+    ts = range(T(-1), T(1); length = 21)
+    B = @inferred BSplineBasis(ord, ts)
     k = order(B)
+    @test B isa AbstractBSplineBasis{k, T}
+
     if isodd(k)
         # `Natural` boundary condition only supported for even-order splines (got k = $k)
         @test_throws ArgumentError RecombinedBSplineBasis(B, Natural())
@@ -19,9 +22,12 @@ function test_natural(ord::BSplineOrder)
     rng = MersenneTwister(42)
 
     R = @inferred RecombinedBSplineBasis(B, Natural())
-    S = Spline(R, 1 .+ rand(rng, length(R)))  # coefficients in [1, 2]
+    @test R isa AbstractBSplineBasis{k, T}
+    S = Spline(R, 1 .+ rand(rng, T, length(R)))  # coefficients in [1, 2]
 
     M = recombination_matrix(R)
+    @test M isa AbstractMatrix{T}
+
     if k == 2
         # In this case, basis functions are not really recombined, and the
         # resulting basis is identical to the original one.
@@ -31,7 +37,12 @@ function test_natural(ord::BSplineOrder)
     for x ∈ boundaries(R)
         Sx = @inferred S(x)
         @test Sx > 1
-        ϵ = 2e-8 * Sx  # threshold
+        rtol = if T === Float64
+            2e-8
+        elseif T === Float32
+            2f-4
+        end
+        ϵ = rtol * Sx  # threshold
         @test abs((Derivative(1) * S)(x)) > ϵ  # different from zero
         for n = 2:(k ÷ 2)
             Sder = Derivative(n) * S
@@ -46,6 +57,12 @@ end
 # (in this case the recombination is equivalent to the original B-spline basis).
 @testset "Natural splines" begin
     @testset "k = $k" for k ∈ (2, 4, 5, 6, 8)
-        test_natural(BSplineOrder(k))
+        test_natural(Float64, BSplineOrder(k))
+    end
+    @testset "Float32" begin
+        # Note: results are very wrong with high-order splines (k ≥ 6) when using Float32.
+        # The zero derivative condition is not well verified.
+        test_natural(Float32, BSplineOrder(2))
+        test_natural(Float32, BSplineOrder(4))
     end
 end