Repository: hyrodium/BasicBSpline.jl
Pull request #113 (2021-12-26 16:42:38)

update around `BSplineDerivativeSpace`


Patch:
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 36ba6742..676d2a59 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -9,7 +9,7 @@ export BSplineDerivativeSpace
 export BSplineManifold, CustomBSplineManifold
 export bsplinebasis₊₀, bsplinebasis₋₀, bsplinebasis
 export bsplinebasis′₊₀, bsplinebasis′₋₀, bsplinebasis′
-export bsplinesupport, bsplineunity, changebasis, lower
+export bsplinesupport, bsplineunity, changebasis
 export bsplinebasisall, intervalindex
 export refinement, bsplinespaces, controlpoints
 export isproper, properdim
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index c3003a0c..63229871 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -130,71 +130,6 @@ Modified version.
     )
 end
 
-@doc raw"""
-1st derivative of B-spline basis function.
-Right-sided limit version.
-```math
-\dot{B}_{(i,p,k)}(t)
-=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
-```
-"""
-bsplinebasis′₊₀
-
-
-@doc raw"""
-1st derivative of B-spline basis function.
-Left-sided limit version.
-```math
-\dot{B}_{(i,p,k)}(t)
-=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
-```
-"""
-bsplinebasis′₋₀
-
-@doc raw"""
-1st derivative of B-spline basis function.
-Modified version.
-```math
-\dot{B}_{(i,p,k)}(t)
-=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
-```
-"""
-bsplinebasis′
-
-function bsplinebasis′₊₀(P::BSplineSpace{p}, i::Integer, t::Real)::Float64 where p
-    k = P.knots
-    return p * (
-        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₊₀(lower(P), i, t) / (k[i+p] - k[i]) : 0) -
-        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₊₀(lower(P), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
-    )
-end
-function bsplinebasis′₊₀(P::BSplineSpace{0}, i::Integer, t::Real)::Float64
-    return 0.0
-end
-
-function bsplinebasis′₋₀(P::BSplineSpace{p}, i::Integer, t::Real)::Float64 where p
-    k = P.knots
-    return p * (
-        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₋₀(lower(P), i, t) / (k[i+p] - k[i]) : 0) -
-        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₋₀(lower(P), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
-    )
-end
-function bsplinebasis′₋₀(P::BSplineSpace{0}, i::Integer, t::Real)::Float64
-    return 0.0
-end
-
-function bsplinebasis′(P::BSplineSpace{p}, i::Integer, t::Real)::Float64 where p
-    k = P.knots
-    return p * (
-        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis(lower(P), i, t) / (k[i+p] - k[i]) : 0) -
-        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis(lower(P), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
-    )
-end
-function bsplinebasis′(P::BSplineSpace{0}, i::Integer, t::Real)::Float64
-    return 0.0
-end
-
-
 """
 TODO: Add docstring
 """
@@ -216,18 +151,18 @@ end
 @generated function bsplinebasisall(P::BSplineSpace{p,T}, i::Integer, t::T) where {p, T}
     bs = [Symbol(:b,i) for i in 1:p]
     Bs = [Symbol(:B,i) for i in 1:p+1]
-    Ks = [:((t-k[i+$(j)])/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
-    Ls = [:((k[i+$(p+j)]-t)/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
+    K1s = [:((k[i+$(p+j)]-t)/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
+    K2s = [:((t-k[i+$(j)])/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
     b = Expr(:tuple, bs...)
     B = Expr(:tuple, Bs...)
-    exs = [:($(Bs[j+1]) = ($(Ls[j+1])*$(bs[j+1]) + $(Ks[j])*$(bs[j]))) for j in 1:p-1]
+    exs = [:($(Bs[j+1]) = ($(K1s[j+1])*$(bs[j+1]) + $(K2s[j])*$(bs[j]))) for j in 1:p-1]
     Expr(:block,
         :($(Expr(:meta, :inline))),
         :(k = knots(P)),
-        :($b = bsplinebasisall(lower(P),i+1,t)),
-        :($(Bs[1]) = $(Ls[1])*$(bs[1])),
+        :($b = bsplinebasisall(_lower(P),i+1,t)),
+        :($(Bs[1]) = $(K1s[1])*$(bs[1])),
         exs...,
-        :($(Bs[p+1]) = $(Ks[p])*$(bs[p])),
+        :($(Bs[p+1]) = $(K2s[p])*$(bs[p])),
         :(return $(B))
     )
 end
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index d20c2c0a..c5b1737c 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -143,13 +143,13 @@ Return a B-spline space of one degree lower.
 \mathcal{P}[p,k] \mapsto \mathcal{P}[p-1,k]
 ```
 """
-lower
+_lower
 
-lower(::Type{AbstractBSplineSpace{p}}) where p = AbstractBSplineSpace{p-1}
-lower(::Type{AbstractBSplineSpace{p,T}}) where {p,T} = AbstractBSplineSpace{p-1,T}
-lower(::Type{BSplineSpace{p}}) where p = BSplineSpace{p-1}
-lower(::Type{BSplineSpace{p,T}}) where {p,T} = BSplineSpace{p-1,T}
-lower(P::BSplineSpace{p,T}) where {p,T} = BSplineSpace{p-1}(knots(P))
+_lower(::Type{AbstractBSplineSpace{p}}) where p = AbstractBSplineSpace{p-1}
+_lower(::Type{AbstractBSplineSpace{p,T}}) where {p,T} = AbstractBSplineSpace{p-1,T}
+_lower(::Type{BSplineSpace{p}}) where p = BSplineSpace{p-1}
+_lower(::Type{BSplineSpace{p,T}}) where {p,T} = BSplineSpace{p-1,T}
+_lower(P::BSplineSpace{p,T}) where {p,T} = BSplineSpace{p-1}(knots(P))
 
 """
 TODO: add docstring
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index 87aa3fff..535756d1 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -22,10 +22,10 @@ function _changebasis_R(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T})::Matri
         return A⁰
     end
 
-    Aᵖ⁻¹ = _changebasis_R(lower(P), lower(P′)) # (n+1) × (n′+1) matrix
+    Aᵖ⁻¹ = _changebasis_R(_lower(P), _lower(P′)) # (n+1) × (n′+1) matrix
     n = dim(P)
     n′ = dim(P′)
-    Z = iszeros(lower(P′))
+    Z = iszeros(_lower(P′))
     W = findall(Z)
     K′ = [k′[i+p′] - k′[i] for i in 1:n′+1]
     K = [ifelse(k[i+p] ≠ k[i], 1 / (k[i+p] - k[i]), 0.0) for i in 1:n+1]
diff --git a/src/_Derivative.jl b/src/_Derivative.jl
index 997d9f40..124ca562 100644
--- a/src/_Derivative.jl
+++ b/src/_Derivative.jl
@@ -7,11 +7,12 @@ struct BSplineDerivativeSpace{r, T<:AbstractBSplineSpace}
     end
 end
 
-bsplinespace(P::BSplineDerivativeSpace) = P.bsplinespace
-
-knots(P::BSplineDerivativeSpace) = knots(bsplinespace(P))
-
-degree(P::BSplineDerivativeSpace{r,<:AbstractBSplineSpace{p}}) where {r,p} = p - r
+bsplinespace(dP::BSplineDerivativeSpace) = dP.bsplinespace
+knots(dP::BSplineDerivativeSpace) = knots(bsplinespace(dP))
+degree(dP::BSplineDerivativeSpace{r,<:AbstractBSplineSpace{p}}) where {r,p} = p - r
+dim(dP::BSplineDerivativeSpace{r,<:AbstractBSplineSpace{p}}) where {r,p} = dim(bsplinespace(dP)) - r
+intervalindex(dP::BSplineDerivativeSpace,t::Real) = intervalindex(bsplinespace(dP),t)
+_lower(dP::BSplineDerivativeSpace{r}) where r = BSplineDerivativeSpace{r-1}(_lower(bsplinespace(dP)))
 
 @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]
@@ -119,6 +120,77 @@ end
     end
 end
 
-# TODO: Add bsplinebasisall(::BSplineDerivativeSpace, i, t)
+@doc raw"""
+1st derivative of B-spline basis function.
+Right-sided limit version.
+```math
+\dot{B}_{(i,p,k)}(t)
+=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
+```
+"""
+bsplinebasis′₊₀
+
+
+@doc raw"""
+1st derivative of B-spline basis function.
+Left-sided limit version.
+```math
+\dot{B}_{(i,p,k)}(t)
+=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
+```
+"""
+bsplinebasis′₋₀
+
+@doc raw"""
+1st derivative of B-spline basis function.
+Modified version.
+```math
+\dot{B}_{(i,p,k)}(t)
+=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
+```
+"""
+bsplinebasis′
+
+for suffix in ("", "₋₀", "₊₀")
+    fname = Symbol(:bsplinebasis, suffix)
+    fname′ = Symbol(:bsplinebasis, "′" ,suffix)
+    @eval function $(fname′)(P::AbstractBSplineSpace, i::Integer, t::Real)
+        return $(fname)(BSplineDerivativeSpace{1}(P), i, t)
+    end
+    @eval function $(fname′)(dP::BSplineDerivativeSpace{r}, i::Integer, t::Real) where r
+        P = bsplinespace(dP)
+        return $(fname)(BSplineDerivativeSpace{r+1}(P), i, t)
+    end
+end
+
+@generated function bsplinebasisall(P::BSplineDerivativeSpace{r,BSplineSpace{p,T}}, i::Integer, t::Real) where {r, p, T}
+    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]
+    K2s = [:($(p)/(k[i+$(p+j)]-k[i+$(j)])) for j in 1:p]
+    b = Expr(:tuple, bs...)
+    B = Expr(:tuple, Bs...)
+    exs = [:($(Bs[j+1]) = ($(K1s[j+1])*$(bs[j+1]) + $(K2s[j])*$(bs[j]))) for j in 1:p-1]
+    if r ≤ p
+        Expr(:block,
+            :($(Expr(:meta, :inline))),
+            :(k = knots(P)),
+            :($b = bsplinebasisall(_lower(P),i+1,t)),
+            :($(Bs[1]) = $(K1s[1])*$(bs[1])),
+            exs...,
+            :($(Bs[p+1]) = $(K2s[p])*$(bs[p])),
+            :(return $(B))
+        )
+    else
+        Z = Expr(:tuple, [:(zero($T)) for i in 1:p+1]...)
+        :(return $(Z))
+    end
+end
+
+function bsplinebasisall(dP::BSplineDerivativeSpace{0,BSplineSpace{p,T}}, i::Integer, t::Real) where {p, T}
+    P = bsplinespace(dP)
+    bsplinebasisall(P,i,t)
+end
+
 # TODO: Add issubset(::BSplineDerivativeSpace, i, t)
 # TODO: Add bsplineunity(::BSplineDerivativeSpace)


Test patch:
diff --git a/test/test_BSplineBasis.jl b/test/test_BSplineBasis.jl
index 085b8cf1..3002944d 100644
--- a/test/test_BSplineBasis.jl
+++ b/test/test_BSplineBasis.jl
@@ -39,21 +39,6 @@ end
             end
         end
 
-        # Check the values of the derivative of B-spline basis funcitons
-        for i in 1:n
-            for t in rand(5)
-                @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′(P, i, t))
-            end
-            for t in k
-                t₊₀ = nextfloat(t)
-                t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
-            end
-        end
-
         # Check the values of the derivative of B-spline basis funcitons
         for t in k
             s = sum([bsplinebasis(P, i, t) for i in 1:n])
@@ -73,23 +58,9 @@ end
                 @test s₋₀ ≈ 1
             end
         end
-
-        # Check the derivative
-        @testset "derivative" begin
-            Random.seed!(57)
-            for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(P, i, t)
-                b′₋₀ = bsplinebasis′₋₀(P, i, t)
-                b′   = bsplinebasis′(P, i, t)
-                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
-                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
-            end
-        end
     end
 
-    @testset "_bsplinebasis" begin
+    @testset "bsplinebasisall" begin
         Random.seed!(42)
 
         k = Knots(rand(10).-1) + Knots(rand(10)) + Knots(rand(10).+1)
@@ -99,28 +70,18 @@ end
             P = BSplineSpace{p}(k)
             for t in ts
                 j = intervalindex(P,t)
-                _B = bsplinebasisall(P,j,t)
-                B = [bsplinebasis(P,i,t) for i in j:j+p]
-                B′ = [bsplinebasis′(P,i,t) for i in j:j+p]
-                @test norm(collect(_B) - B) < ε
+                B = collect(bsplinebasisall(P,j,t))
+                _B = [bsplinebasis(P,i,t) for i in j:j+p]
+                @test _B ≈ B
 
                 _B = [bsplinebasis(P,i,t) for i in j:j+p]
-                @test norm(_B - B) < ε
+                @test _B ≈ B
 
                 _B = [bsplinebasis₊₀(P,i,t) for i in j:j+p]
-                @test norm(_B - B) < ε
+                @test _B ≈ B
 
                 _B = [bsplinebasis₋₀(P,i,t) for i in j:j+p]
-                @test norm(_B - B) < ε
-
-                _B′ =[ bsplinebasis′(P,i,t) for i in j:j+p]
-                @test norm(_B′ - B′) < ε
-
-                _B′ =[ bsplinebasis′₊₀(P,i,t) for i in j:j+p]
-                @test norm(_B′ - B′) < ε
-
-                _B′ =[ bsplinebasis′₋₀(P,i,t) for i in j:j+p]
-                @test norm(_B′ - B′) < ε
+                @test _B ≈ B
             end
         end
     end
diff --git a/test/test_Derivative.jl b/test/test_Derivative.jl
index 09468362..d3c6770b 100644
--- a/test/test_Derivative.jl
+++ b/test/test_Derivative.jl
@@ -1,23 +1,90 @@
 @testset "Derivative" begin
-    p_max = 5
+    Random.seed!(42)
+    p_max = 4
     ts = rand(10)
     dt = 1e-7
 
-    for p in 1:p_max
+    # Not sure why this @testset doesn't work fine.
+    # @testset "$(p)-th degree basis" for p in 0:p_max
+    for p in 0: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)
+            @test bsplinebasis₋₀(P,i,t) == bsplinebasis₋₀(P0,i,t)
+            @test bsplinebasis₊₀(P,i,t) == bsplinebasis₊₀(P0,i,t)
         end
-        
+        @test dim(P) == dim(P0)
+
+        # Check the values of the derivative of B-spline basis funcitons
+        n = dim(P0)
+        for t in k
+            s = sum([bsplinebasis(P0, i, t) for i in 1:n])
+            s₊₀ = sum([bsplinebasis₊₀(P0, i, t) for i in 1:n])
+            s₋₀ = sum([bsplinebasis₋₀(P0, i, t) for i in 1:n])
+            if t == k[1]
+                @test s ≈ 1
+                @test s₊₀ ≈ 1
+                @test s₋₀ == 0
+            elseif t == k[end]
+                @test s ≈ 1
+                @test s₊₀ == 0
+                @test s₋₀ ≈ 1
+            else
+                @test s ≈ 1
+                @test s₊₀ ≈ 1
+                @test s₋₀ ≈ 1
+            end
+        end
+
+        P1 = BSplineDerivativeSpace{1}(P)
+        for t in ts, i in 1:dim(P)
+            @test bsplinebasis′(P,i,t) == bsplinebasis(P1,i,t)
+        end
+
         for r in 1:p_max
-            P1 = BSplineDerivativeSpace{r-1}(P)
-            P2 = BSplineDerivativeSpace{r}(P)        
+            Pa = BSplineDerivativeSpace{r}(P)
+            Pb = BSplineDerivativeSpace{r-1}(P)
+            @test degree(Pa) == p-r
+            @test dim(Pa) == dim(P)-r
             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
+                d1 = bsplinebasis(Pa,i,t)
+                d2 = (bsplinebasis(Pb,i,t+dt) - bsplinebasis(Pb,i,t-dt))/2dt
+                d3 = bsplinebasis′(Pb,i,t)
+                @test d1 ≈ d2
+                @test d1 == d3
+            end
+        end
+    end
+
+    @testset "bsplinebasisall" begin
+        Random.seed!(42)
+
+        k = Knots(rand(10).-1) + Knots(rand(10)) + Knots(rand(10).+1)
+        ts = rand(10)
+
+        for p in 0:5
+            P = BSplineSpace{p}(k)
+
+            for r in 0:5
+                dP = BSplineDerivativeSpace{r}(P)
+
+                for t in ts
+                    j = intervalindex(dP,t)
+                    B = collect(bsplinebasisall(dP,j,t))
+                    _B = [bsplinebasis(dP,i,t) for i in j:j+p]
+                    @test _B ≈ B
+
+                    _B = [bsplinebasis(dP,i,t) for i in j:j+p]
+                    @test _B ≈ B
+
+                    _B = [bsplinebasis₊₀(dP,i,t) for i in j:j+p]
+                    @test _B ≈ B
+
+                    _B = [bsplinebasis₋₀(dP,i,t) for i in j:j+p]
+                    @test _B ≈ B
+                end
             end
         end
     end