Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
Fix type promotion in `bsplinebasis`
```julia
julia> k = KnotVector{Int}(1:12)
KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])

julia> P = BSplineSpace{3}(k)
BSplineSpace{3, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]))

julia> bsplinebasis(P,1,1)
ERROR: InexactError: Int64(-0.5)
Stacktrace:
 [1] Int64
   @ ./float.jl:812 [inlined]
 [2] _d
   @ ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:3 [inlined]
 [3] macro expansion
   @ ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:110 [inlined]
 [4] bsplinebasis(P::BSplineSpace{3, Int64}, i::Int64, t::Int64)
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_BSplineBasis.jl:110
 [5] top-level scope
   @ REPL[11]:1

```

Pull request:
This PR contains the following changes:
* Fix #133
  * This is done with `StaticArrays.arithmetic_closure()`
* The return type of `bsplinebasisall` function is changed from `Tuple` to `SVector`
  * This may be a breaking change, so the next release will be `v0.4.0`.
  * This doesn't produce any performance regression.
  * This change is necessary for auto-diff system. (#124)

Patch:
diff --git a/Project.toml b/Project.toml
index 7caf9513a..75be5e0a5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.3.2"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 FastGaussQuadrature = "0.4"
@@ -14,9 +15,9 @@ IntervalSets = "0.5"
 julia = "^1.6"
 
 [extras]
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 
 [targets]
 test = ["Test", "Random", "GeometryBasics"]
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 28d6cfa01..a4d4e7ff1 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -2,6 +2,7 @@ module BasicBSpline
 
 using LinearAlgebra
 using IntervalSets
+using StaticArrays
 using FastGaussQuadrature
 
 export KnotVector, AbstractKnotVector
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index 166188cce..199532e7e 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -1,6 +1,6 @@
 # B-Spline Basis Function
 
-@inline _d(a::T,b::T) where T = ifelse(iszero(b), zero(T), T(a/b))
+@inline _d(a::T,b::T) where T = (U=StaticArrays.arithmetic_closure(T); ifelse(iszero(b), zero(U), U(a/b)))
 @inline _d(a,b) = _d(promote(a,b)...)
 
 @doc raw"""
@@ -21,7 +21,7 @@ Right-sided limit version.
 \end{aligned}
 """
 @generated function bsplinebasis₊₀(P::BSplineSpace{p,T}, i::Integer, t::S) where {p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -64,7 +64,7 @@ Left-sided limit version.
 \end{aligned}
 """
 @generated function bsplinebasis₋₀(P::BSplineSpace{p,T}, i::Integer, t::S) where {p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -108,7 +108,7 @@ Modified version.
 \end{aligned}
 """
 @generated function bsplinebasis(P::BSplineSpace{p,T}, i::Integer, t::S) where {p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -139,15 +139,15 @@ TODO: Add docstring
 bsplinebasisall
 
 @inline function bsplinebasisall(P::BSplineSpace{0,T},i::Integer,t::S) where {T, S<:Real}
-    U = promote_type(T,S)
-    (one(U),)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
+    SVector(one(U),)
 end
 
 @inline function bsplinebasisall(P::BSplineSpace{1}, i::Integer, t::Real)
     k = knotvector(P)
     B1 = (k[i+2]-t)/(k[i+2]-k[i+1])
     B2 = (t-k[i+1])/(k[i+2]-k[i+1])
-    return (B1, B2)
+    return SVector(B1, B2)
 end
 
 @generated function bsplinebasisall(P::BSplineSpace{p}, i::Integer, t::Real) where p
@@ -165,6 +165,6 @@ end
         :($(Bs[1]) = $(K1s[1])*$(bs[1])),
         exs...,
         :($(Bs[p+1]) = $(K2s[p])*$(bs[p])),
-        :(return $(B))
+        :(return SVector($(B)))
     )
 end
diff --git a/src/_DerivativeBasis.jl b/src/_DerivativeBasis.jl
index 974c3a473..dcd3ee212 100644
--- a/src/_DerivativeBasis.jl
+++ b/src/_DerivativeBasis.jl
@@ -1,7 +1,7 @@
 # Derivative of B-spline basis function
 
 @generated function bsplinebasis₊₀(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -37,7 +37,7 @@
 end
 
 @generated function bsplinebasis₋₀(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -68,12 +68,12 @@ end
             :(return $(prod(p-r+1:p))*B1)
         )
     else
-        :(return zero(T))
+        :(return zero($U))
     end
 end
 
 @generated function bsplinebasis(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     ks = [Symbol(:k,i) for i in 1:p+2]
     Ks = [Symbol(:K,i) for i in 1:p+1]
     Bs = [Symbol(:B,i) for i in 1:p+1]
@@ -104,7 +104,7 @@ end
             :(return $(prod(p-r+1:p))*B1)
         )
     else
-        :(return zero(T))
+        :(return zero($U))
     end
 end
 
@@ -152,7 +152,7 @@ for suffix in ("", "₋₀", "₊₀")
 end
 
 @generated function bsplinebasisall(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
-    U = promote_type(T,S)
+    U = StaticArrays.arithmetic_closure(promote_type(T,S))
     bs = [Symbol(:b,i) for i in 1:p]
     Bs = [Symbol(:B,i) for i in 1:p+1]
     K1s = [:($(p)/(k[i+$(j)]-k[i+$(p+j)])) for j in 1:p]
@@ -168,11 +168,11 @@ end
             :($(Bs[1]) = $(K1s[1])*$(bs[1])),
             exs...,
             :($(Bs[p+1]) = $(K2s[p])*$(bs[p])),
-            :(return $(B))
+            :(return SVector($(B)))
         )
     else
         Z = Expr(:tuple, [:(zero($U)) for i in 1:p+1]...)
-        :(return $(Z))
+        :(return SVector($(Z)))
     end
 end
 


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 908dab5e1..174f269c1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,7 @@ using IntervalSets
 using LinearAlgebra
 using Test
 using Random
+using StaticArrays
 using GeometryBasics
 
 @testset "BasicBSpline.jl" begin
diff --git a/test/test_BSplineBasis.jl b/test/test_BSplineBasis.jl
index 39aa2a30f..bcd79e3d4 100644
--- a/test/test_BSplineBasis.jl
+++ b/test/test_BSplineBasis.jl
@@ -93,11 +93,34 @@ end
         bsplinebasis₊₀(P,1,11//5) isa Rational{Int}
         bsplinebasis₋₀(P,1,11//5) isa Rational{Int}
 
-        bsplinebasis(P,1,11//5)   ===
+        @test bsplinebasis(P,1,11//5) ===
         bsplinebasis₊₀(P,1,11//5) ===
         bsplinebasis₋₀(P,1,11//5) === 106//375
     end
 
+    @testset "Check type" begin
+        k = KnotVector{Int}(1:12)
+        P0 = BSplineSpace{0}(k)
+        P1 = BSplineSpace{1}(k)
+        P2 = BSplineSpace{2}(k)
+
+        @test bsplinebasis(P0,1,5) isa Float64
+        @test bsplinebasis(P1,1,5) isa Float64
+        @test bsplinebasis(P2,1,5) isa Float64
+
+        @test bsplinebasis₊₀(P0,1,5) isa Float64
+        @test bsplinebasis₊₀(P1,1,5) isa Float64
+        @test bsplinebasis₊₀(P2,1,5) isa Float64
+
+        @test bsplinebasis₋₀(P0,1,5) isa Float64
+        @test bsplinebasis₋₀(P1,1,5) isa Float64
+        @test bsplinebasis₋₀(P2,1,5) isa Float64
+
+        @test bsplinebasisall(P0,1,5) isa SVector{1,Float64}
+        @test bsplinebasisall(P1,1,5) isa SVector{2,Float64}
+        @test bsplinebasisall(P2,1,5) isa SVector{3,Float64}
+    end
+
     @testset "Endpoints" begin
         p = 2
 
diff --git a/test/test_Derivative.jl b/test/test_Derivative.jl
index 27579d913..bf7290b71 100644
--- a/test/test_Derivative.jl
+++ b/test/test_Derivative.jl
@@ -99,7 +99,7 @@
         bsplinebasis₊₀(dP,1,11//5) isa Rational{Int}
         bsplinebasis₋₀(dP,1,11//5) isa Rational{Int}
 
-        bsplinebasis(dP,1,11//5)   ===
+        @test bsplinebasis(dP,1,11//5) ===
         bsplinebasis₊₀(dP,1,11//5) ===
         bsplinebasis₋₀(dP,1,11//5) ===
         bsplinebasis′(P,1,11//5)   ===
@@ -107,6 +107,34 @@
         bsplinebasis′₋₀(P,1,11//5) === 16//25
     end
 
+    @testset "Check type" begin
+        k = KnotVector{Int}(1:12)
+        P0 = BSplineSpace{0}(k)
+        P1 = BSplineSpace{1}(k)
+        P2 = BSplineSpace{2}(k)
+
+        for r in 0:2
+            dP0 = BSplineDerivativeSpace{r}(P0)
+            dP1 = BSplineDerivativeSpace{r}(P1)
+            dP2 = BSplineDerivativeSpace{r}(P2)
+            @test bsplinebasis(dP0,1,5) isa Float64
+            @test bsplinebasis(dP1,1,5) isa Float64
+            @test bsplinebasis(dP2,1,5) isa Float64
+
+            @test bsplinebasis₊₀(dP0,1,5) isa Float64
+            @test bsplinebasis₊₀(dP1,1,5) isa Float64
+            @test bsplinebasis₊₀(dP2,1,5) isa Float64
+
+            @test bsplinebasis₋₀(dP0,1,5) isa Float64
+            @test bsplinebasis₋₀(dP1,1,5) isa Float64
+            @test bsplinebasis₋₀(dP2,1,5) isa Float64
+
+            @test bsplinebasisall(dP0,1,5) isa SVector{1,Float64}
+            @test bsplinebasisall(dP1,1,5) isa SVector{2,Float64}
+            @test bsplinebasisall(dP2,1,5) isa SVector{3,Float64}
+        end
+    end
+
     @testset "Endpoints" begin
         p = 2