Repository: hyrodium/BasicBSpline.jl
Pull request #112 (2021-12-26 06:24:01)

Add `BSplineDerivativeSpace`


Patch:
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 462982aa..36ba6742 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -5,6 +5,7 @@ using IntervalSets
 using FastGaussQuadrature
 
 export Knots, BSplineSpace, dim, ⊑, ⊒, ⋢, ⋣, ≃
+export BSplineDerivativeSpace
 export BSplineManifold, CustomBSplineManifold
 export bsplinebasis₊₀, bsplinebasis₋₀, bsplinebasis
 export bsplinebasis′₊₀, bsplinebasis′₋₀, bsplinebasis′
@@ -20,6 +21,7 @@ export fittingcontrolpoints
 include("_Knots.jl")
 include("_BSplineSpace.jl")
 include("_BSplineBasis.jl")
+include("_Derivative.jl")
 include("_ChangeBasis.jl")
 include("_BSplineManifold.jl")
 include("_CustomBSplineManifold.jl")
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index 086be8b2..c3003a0c 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -19,7 +19,7 @@ Right-sided limit version.
 \end{cases}
 \end{aligned}
 """
-@generated function bsplinebasis₊₀(P::BSplineSpace{p,T}, i::Integer, t::Real) where {p, T}
+@generated function bsplinebasis₊₀(P::BSplineSpace{p,T}, i::Integer, t::T) where {p, T}
     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]
@@ -61,7 +61,7 @@ Left-sided limit version.
 \end{cases}
 \end{aligned}
 """
-@generated function bsplinebasis₋₀(P::BSplineSpace{p,T}, i::Integer, t::Real) where {p, T}
+@generated function bsplinebasis₋₀(P::BSplineSpace{p,T}, i::Integer, t::T) where {p, T}
     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 @@ Modified version.
 \end{cases}
 \end{aligned}
 """
-@generated function bsplinebasis(P::BSplineSpace{p,T}, i::Integer, t::Real) where {p, T}
+@generated function bsplinebasis(P::BSplineSpace{p,T}, i::Integer, t::T) where {p, T}
     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]
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index 02db4b97..87aa3fff 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -169,3 +169,5 @@ end
 function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
     changebasis(BSplineSpace(P),BSplineSpace(P′))
 end
+
+# TODO: Add changebasis(::BSplineDerivativeSpace, )
diff --git a/src/_Derivative.jl b/src/_Derivative.jl
new file mode 100644
index 00000000..997d9f40
--- /dev/null
+++ b/src/_Derivative.jl
@@ -0,0 +1,124 @@
+# Derivative
+
+struct BSplineDerivativeSpace{r, T<:AbstractBSplineSpace}
+    bsplinespace::T
+    function BSplineDerivativeSpace{r}(P::T) where {r, T<:AbstractBSplineSpace}
+        new{r,T}(P)
+    end
+end
+
+bsplinespace(P::BSplineDerivativeSpace) = P.bsplinespace
+
+knots(P::BSplineDerivativeSpace) = knots(bsplinespace(P))
+
+degree(P::BSplineDerivativeSpace{r,<:AbstractBSplineSpace{p}}) where {r,p} = p - r
+
+@generated function bsplinebasis₊₀(P::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+    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]
+    k_l = Expr(:tuple, ks...)
+    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]...)
+    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]...)
+    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
+        exs = Expr[]
+        for i in 1:p-r
+            push!(exs, :($(K_l(p+2-i)) = $(K_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(B_r(p+1-i))))
+        end
+        for i in p-r+1:p
+            push!(exs, :($(K_l(p+2-i)) = $(L_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(C_r(p+1-i))))
+        end
+        Expr(:block,
+            :(v = knots(P).vector),
+            :($k_l = $k_r),
+            :($(B_l(p+1)) = $(A_r(p+1))),
+            exs...,
+            :(return $(prod(p-r+1:p))*B1)
+        )
+    else
+        :(return zero(T))
+    end
+end
+
+@generated function bsplinebasis₋₀(P::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+    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]
+    k_l = Expr(:tuple, ks...)
+    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]...)
+    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]...)
+    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
+        exs = Expr[]
+        for i in 1:p-r
+            push!(exs, :($(K_l(p+2-i)) = $(K_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(B_r(p+1-i))))
+        end
+        for i in p-r+1:p
+            push!(exs, :($(K_l(p+2-i)) = $(L_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(C_r(p+1-i))))
+        end
+        Expr(:block,
+            :(v = knots(P).vector),
+            :($k_l = $k_r),
+            :($(B_l(p+1)) = $(A_r(p+1))),
+            exs...,
+            :(return $(prod(p-r+1:p))*B1)
+        )
+    else
+        :(return zero(T))
+    end
+end
+
+@generated function bsplinebasis(P::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+    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]
+    k_l = Expr(:tuple, ks...)
+    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]...)
+    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]...)
+    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
+        exs = Expr[]
+        for i in 1:p-r
+            push!(exs, :($(K_l(p+2-i)) = $(K_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(B_r(p+1-i))))
+        end
+        for i in p-r+1:p
+            push!(exs, :($(K_l(p+2-i)) = $(L_r(i,p+2-i))))
+            push!(exs, :($(B_l(p+1-i)) = $(C_r(p+1-i))))
+        end
+        Expr(:block,
+            :(v = knots(P).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)
+        )
+    else
+        :(return zero(T))
+    end
+end
+
+# TODO: Add bsplinebasisall(::BSplineDerivativeSpace, i, t)
+# TODO: Add issubset(::BSplineDerivativeSpace, i, t)
+# TODO: Add bsplineunity(::BSplineDerivativeSpace)


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 61fe0ad0..4eaf00da 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,7 @@ using GeometryBasics
     include("test_Knots.jl")
     include("test_BSplineSpace.jl")
     include("test_BSplineBasis.jl")
+    include("test_Derivative.jl")
     include("test_ChangeBasis.jl")
     include("test_BSplineManifold.jl")
     include("test_Integral.jl")
diff --git a/test/test_Derivative.jl b/test/test_Derivative.jl
new file mode 100644
index 00000000..09468362
--- /dev/null
+++ b/test/test_Derivative.jl
@@ -0,0 +1,24 @@
+@testset "Derivative" begin
+    p_max = 5
+    ts = rand(10)
+    dt = 1e-7
+
+    for p in 1:p_max
+        k = Knots(rand(20)) + (p+1)*Knots(0,1)
+        P = BSplineSpace{p}(k)
+        P0 = BSplineDerivativeSpace{0}(P)
+        for t in ts, i in 1:dim(P)
+            @test bsplinebasis(P,i,t) == bsplinebasis(P0,i,t)
+        end
+        
+        for r in 1:p_max
+            P1 = BSplineDerivativeSpace{r-1}(P)
+            P2 = BSplineDerivativeSpace{r}(P)        
+            for t in ts, i in 1:dim(P)
+                d1 = (bsplinebasis(P1,6,t+dt) - bsplinebasis(P1,6,t-dt))/2dt
+                d2 = bsplinebasis(P2,6,t)
+                @test d1 ≈ d2 atol=1e-4
+            end
+        end
+    end
+end