Repository: hyrodium/BasicBSpline.jl
Pull request #130 (2021-12-31 01:59:25)

Fix `bsplinebasis` function on endpoints


Patch:
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index 136b2c1f..166188cc 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -116,7 +116,7 @@ Modified version.
     k_r = Expr(:tuple, :(v[i]), (:(v[i+$j]) for j in 1:p+1)...)
     K_l(n) = Expr(:tuple, Ks[1:n]...)
     B_l(n) = Expr(:tuple, Bs[1:n]...)
-    A_r(n) = Expr(:tuple, [:($U($(ks[i])≤t<$(ks[i+1]))) for i in 1:n]...)
+    A_r(n) = Expr(:tuple, [:($U(($(ks[i])≤t<$(ks[i+1])) || ($(ks[i])<t==$(ks[i+1])==v[end]))) for i in 1:n]...)
     K_r(m,n) = Expr(:tuple, [:(_d(t-$(ks[i]),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
     B_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])+(1-$(Ks[i+1]))*$(Bs[i+1])) for i in 1:n]...)
     exs = Expr[]
@@ -128,7 +128,6 @@ Modified version.
         :(v = knotvector(P).vector),
         :($k_l = $k_r),
         :($(B_l(p+1)) = $(A_r(p+1))),
-        :($(Bs[end]) += $(U)(t == $(ks[end]) == v[end])),
         exs...,
         :(return B1)
     )
@@ -139,8 +138,6 @@ TODO: Add docstring
 """
 bsplinebasisall
 
-# TODO: faster implementation
-# TODO: use @generated macro
 @inline function bsplinebasisall(P::BSplineSpace{0,T},i::Integer,t::S) where {T, S<:Real}
     U = promote_type(T,S)
     (one(U),)
diff --git a/src/_DerivativeBasis.jl b/src/_DerivativeBasis.jl
index 013a46e2..974c3a47 100644
--- a/src/_DerivativeBasis.jl
+++ b/src/_DerivativeBasis.jl
@@ -1,6 +1,7 @@
 # Derivative of B-spline basis function
 
-@generated function bsplinebasis₊₀(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+@generated function bsplinebasis₊₀(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
+    U = 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]
@@ -8,9 +9,9 @@
     k_r = Expr(:tuple, :(v[i]), (:(v[i+$j]) for j in 1:p+1)...)
     K_l(n) = Expr(:tuple, Ks[1:n]...)
     B_l(n) = Expr(:tuple, Bs[1:n]...)
-    A_r(n) = Expr(:tuple, [:(T($(ks[i])≤t<$(ks[i+1]))) for i in 1:n]...)
+    A_r(n) = Expr(:tuple, [:($U($(ks[i])≤t<$(ks[i+1]))) for i in 1:n]...)
     K_r(m,n) = Expr(:tuple, [:(_d(t-$(ks[i]),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
-    L_r(m,n) = Expr(:tuple, [:(_d(one(T),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
+    L_r(m,n) = Expr(:tuple, [:(_d(one($U),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
     B_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])+(1-$(Ks[i+1]))*$(Bs[i+1])) for i in 1:n]...)
     C_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])-$(Ks[i+1])*$(Bs[i+1])) for i in 1:n]...)
     if r ≤ p
@@ -31,11 +32,12 @@
             :(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::Real) where {r, p, T}
+@generated function bsplinebasis₋₀(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
+    U = 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]
@@ -43,9 +45,9 @@ end
     k_r = Expr(:tuple, :(v[i]), (:(v[i+$j]) for j in 1:p+1)...)
     K_l(n) = Expr(:tuple, Ks[1:n]...)
     B_l(n) = Expr(:tuple, Bs[1:n]...)
-    A_r(n) = Expr(:tuple, [:(T($(ks[i])<t≤$(ks[i+1]))) for i in 1:n]...)
+    A_r(n) = Expr(:tuple, [:($U($(ks[i])<t≤$(ks[i+1]))) for i in 1:n]...)
     K_r(m,n) = Expr(:tuple, [:(_d(t-$(ks[i]),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
-    L_r(m,n) = Expr(:tuple, [:(_d(one(T),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
+    L_r(m,n) = Expr(:tuple, [:(_d(one($U),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
     B_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])+(1-$(Ks[i+1]))*$(Bs[i+1])) for i in 1:n]...)
     C_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])-$(Ks[i+1])*$(Bs[i+1])) for i in 1:n]...)
     if r ≤ p
@@ -70,7 +72,8 @@ end
     end
 end
 
-@generated function bsplinebasis(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+@generated function bsplinebasis(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
+    U = 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]
@@ -78,9 +81,9 @@ end
     k_r = Expr(:tuple, :(v[i]), (:(v[i+$j]) for j in 1:p+1)...)
     K_l(n) = Expr(:tuple, Ks[1:n]...)
     B_l(n) = Expr(:tuple, Bs[1:n]...)
-    A_r(n) = Expr(:tuple, [:(T($(ks[i])≤t<$(ks[i+1]))) for i in 1:n]...)
+    A_r(n) = Expr(:tuple, [:($U(($(ks[i])≤t<$(ks[i+1])) || ($(ks[i])<t==$(ks[i+1])==v[end]))) for i in 1:n]...)
     K_r(m,n) = Expr(:tuple, [:(_d(t-$(ks[i]),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
-    L_r(m,n) = Expr(:tuple, [:(_d(one(T),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
+    L_r(m,n) = Expr(:tuple, [:(_d(one($U),$(ks[i+m])-$(ks[i]))) for i in 1:n]...)
     B_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])+(1-$(Ks[i+1]))*$(Bs[i+1])) for i in 1:n]...)
     C_r(n) = Expr(:tuple, [:($(Ks[i])*$(Bs[i])-$(Ks[i+1])*$(Bs[i+1])) for i in 1:n]...)
     if r ≤ p
@@ -97,7 +100,6 @@ end
             :(v = knotvector(dP).vector),
             :($k_l = $k_r),
             :($(B_l(p+1)) = $(A_r(p+1))),
-            :($(Bs[end]) += $(T)(t == $(ks[end]) == v[end])),
             exs...,
             :(return $(prod(p-r+1:p))*B1)
         )
@@ -149,7 +151,8 @@ for suffix in ("", "₋₀", "₊₀")
     end
 end
 
-@generated function bsplinebasisall(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+@generated function bsplinebasisall(dP::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::S) where {r, p, T, S<:Real}
+    U = 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,7 +171,7 @@ end
             :(return $(B))
         )
     else
-        Z = Expr(:tuple, [:(zero($T)) for i in 1:p+1]...)
+        Z = Expr(:tuple, [:(zero($U)) for i in 1:p+1]...)
         :(return $(Z))
     end
 end


Test patch:
diff --git a/test/test_BSplineBasis.jl b/test/test_BSplineBasis.jl
index 292a2d0f..39aa2a30 100644
--- a/test/test_BSplineBasis.jl
+++ b/test/test_BSplineBasis.jl
@@ -89,7 +89,43 @@ end
         P = BSplineSpace{3}(k)
         @test k isa KnotVector{Int}
         @test P isa BSplineSpace{3,Int}
-        @test bsplinebasis(P,1,11//5) isa Rational{Int}
-        @test bsplinebasis(P,1,11//5) === bsplinebasis(P,2,16//5) === 106//375
+        bsplinebasis(P,1,11//5) isa Rational{Int}
+        bsplinebasis₊₀(P,1,11//5) isa Rational{Int}
+        bsplinebasis₋₀(P,1,11//5) isa Rational{Int}
+
+        bsplinebasis(P,1,11//5)   ===
+        bsplinebasis₊₀(P,1,11//5) ===
+        bsplinebasis₋₀(P,1,11//5) === 106//375
+    end
+
+    @testset "Endpoints" begin
+        p = 2
+
+        k = KnotVector(1:2) + (p+1)*KnotVector(0,3)
+        P0 = BSplineSpace{0}(k)
+        P1 = BSplineSpace{1}(k)
+        P2 = BSplineSpace{2}(k)
+
+        @test !isproper(P0)
+        @test !isproper(P1)
+        @test isproper(P2)
+
+        n0 = dim(P0)
+        n1 = dim(P1)
+        n2 = dim(P2)
+
+        @test bsplinebasis.(P0,1:n0,0) == bsplinebasis₊₀.(P0,1:n0,0) == [0,0,1,0,0,0,0]
+        @test bsplinebasis.(P1,1:n1,0) == bsplinebasis₊₀.(P1,1:n1,0) == [0,1,0,0,0,0]
+        @test bsplinebasis.(P2,1:n2,0) == bsplinebasis₊₀.(P2,1:n2,0) == [1,0,0,0,0]
+        @test bsplinebasis.(P0,1:n0,3) == bsplinebasis₋₀.(P0,1:n0,3) == [0,0,0,0,1,0,0]
+        @test bsplinebasis.(P1,1:n1,3) == bsplinebasis₋₀.(P1,1:n1,3) == [0,0,0,0,1,0]
+        @test bsplinebasis.(P2,1:n2,3) == bsplinebasis₋₀.(P2,1:n2,3) == [0,0,0,0,1]
+
+        @test bsplinebasis₋₀.(P0,1:n0,0) == [0,0,0,0,0,0,0]
+        @test bsplinebasis₋₀.(P1,1:n1,0) == [0,0,0,0,0,0]
+        @test bsplinebasis₋₀.(P2,1:n2,0) == [0,0,0,0,0]
+        @test bsplinebasis₊₀.(P0,1:n0,3) == [0,0,0,0,0,0,0]
+        @test bsplinebasis₊₀.(P1,1:n1,3) == [0,0,0,0,0,0]
+        @test bsplinebasis₊₀.(P2,1:n2,3) == [0,0,0,0,0]
     end
 end
diff --git a/test/test_Derivative.jl b/test/test_Derivative.jl
index b2648ae1..27579d91 100644
--- a/test/test_Derivative.jl
+++ b/test/test_Derivative.jl
@@ -88,4 +88,53 @@
             end
         end
     end
+
+    @testset "Rational" begin
+        k = KnotVector{Int}(1:12)
+        P = BSplineSpace{3}(k)
+        dP = BSplineDerivativeSpace{1}(P)
+        k isa KnotVector{Int}
+        dP isa BSplineDerivativeSpace{1,BSplineSpace{3,Int}}
+        bsplinebasis(dP,1,11//5) isa Rational{Int}
+        bsplinebasis₊₀(dP,1,11//5) isa Rational{Int}
+        bsplinebasis₋₀(dP,1,11//5) isa Rational{Int}
+
+        bsplinebasis(dP,1,11//5)   ===
+        bsplinebasis₊₀(dP,1,11//5) ===
+        bsplinebasis₋₀(dP,1,11//5) ===
+        bsplinebasis′(P,1,11//5)   ===
+        bsplinebasis′₊₀(P,1,11//5) ===
+        bsplinebasis′₋₀(P,1,11//5) === 16//25
+    end
+
+    @testset "Endpoints" begin
+        p = 2
+
+        k = KnotVector(1:2) + (p+1)*KnotVector(0,3)
+        P0 = BSplineSpace{0}(k)
+        P1 = BSplineSpace{1}(k)
+        P2 = BSplineSpace{2}(k)
+
+        @test !isproper(P0)
+        @test !isproper(P1)
+        @test isproper(P2)
+
+        n0 = dim(P0)
+        n1 = dim(P1)
+        n2 = dim(P2)
+
+        @test bsplinebasis′.(P0,1:n0,0) == bsplinebasis′₊₀.(P0,1:n0,0) == [0,0,0,0,0,0,0]
+        @test bsplinebasis′.(P1,1:n1,0) == bsplinebasis′₊₀.(P1,1:n1,0) == [0,-1,1,0,0,0]
+        @test bsplinebasis′.(P2,1:n2,0) == bsplinebasis′₊₀.(P2,1:n2,0) == [-2,2,0,0,0]
+        @test bsplinebasis′.(P0,1:n0,3) == bsplinebasis′₋₀.(P0,1:n0,3) == [0,0,0,0,0,0,0]
+        @test bsplinebasis′.(P1,1:n1,3) == bsplinebasis′₋₀.(P1,1:n1,3) == [0,0,0,-1,1,0]
+        @test bsplinebasis′.(P2,1:n2,3) == bsplinebasis′₋₀.(P2,1:n2,3) == [0,0,0,-2,2]
+
+        @test bsplinebasis′₋₀.(P0,1:n0,0) == [0,0,0,0,0,0,0]
+        @test bsplinebasis′₋₀.(P1,1:n1,0) == [0,0,0,0,0,0]
+        @test bsplinebasis′₋₀.(P2,1:n2,0) == [0,0,0,0,0]
+        @test bsplinebasis′₊₀.(P0,1:n0,3) == [0,0,0,0,0,0,0]
+        @test bsplinebasis′₊₀.(P1,1:n1,3) == [0,0,0,0,0,0]
+        @test bsplinebasis′₊₀.(P2,1:n2,3) == [0,0,0,0,0]
+    end
 end