Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
Incorrect `changebasis_I`
```julia
julia> using BasicBSpline

julia> P = BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))
BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))

julia> P′ = BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))
BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))

julia> changebasis_I(P, P′)
6×10 SparseArrays.SparseMatrixCSC{Float64, Int32} with 27 stored entries:
   ⋅    2.8   0.4    ⋅    ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅   -4.4   2.0   0.8   ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
  2.8   2.8  -1.6   0.4  0.8  0.45  0.05  0.05  0.05  0.05
 -0.2  -0.2   0.2  -0.2  0.2  0.5   0.5   0.5   0.5   0.5
   ⋅     ⋅     ⋅     ⋅    ⋅   0.05  0.45   ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
```

Pull request description:
This PR fixes https://github.com/hyrodium/BasicBSpline.jl/issues/346.

**Before this PR**
```julia
julia> using BasicBSpline

julia> P = BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))
BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))

julia> P′ = BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))
BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))

julia> changebasis_I(P, P′)
6×10 SparseArrays.SparseMatrixCSC{Float64, Int32} with 27 stored entries:
   ⋅    2.8   0.4    ⋅    ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
   ⋅   -4.4   2.0   0.8   ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
  2.8   2.8  -1.6   0.4  0.8  0.45  0.05  0.05  0.05  0.05
 -0.2  -0.2   0.2  -0.2  0.2  0.5   0.5   0.5   0.5   0.5
   ⋅     ⋅     ⋅     ⋅    ⋅   0.05  0.45   ⋅     ⋅     ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅    ⋅     ⋅     ⋅     ⋅     ⋅ 
```

**After this PR**
```julia
julia> using BasicBSpline

julia> P = BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))
BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))

julia> P′ = BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))
BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))

julia> changebasis_I(P, P′)
6×10 SparseMatrixCSC{Float64, Int32} with 19 stored entries:
  ⋅    2.8   0.4    ⋅    ⋅    ⋅     ⋅     ⋅    ⋅    ⋅ 
  ⋅   -4.4   2.0   0.8   ⋅    ⋅     ⋅     ⋅    ⋅    ⋅ 
  ⋅    2.8  -1.6   0.4  0.8  0.45  0.05   ⋅    ⋅    ⋅ 
  ⋅   -0.2   0.2  -0.2  0.2  0.5   0.5    ⋅    ⋅    ⋅ 
  ⋅     ⋅     ⋅     ⋅    ⋅   0.05  0.45   ⋅    ⋅    ⋅ 
  ⋅     ⋅     ⋅     ⋅    ⋅    ⋅     ⋅     ⋅    ⋅    ⋅
```

Patch:
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index 9d2846848..7932628ba 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -252,19 +252,17 @@ true
 isdegenerate(P::BSplineSpace)
 
 for (f, fnon) in ((:isdegenerate_R, :isnondegenerate_R), (:isdegenerate_I, :isnondegenerate_I))
-    @eval function $f(P::BSplineSpace)
+    @eval function $f(P::AbstractFunctionSpace)
         for i in 1:dim(P)
             $f(P,i) && return true
         end
         return false
     end
-    @eval $fnon(P::BSplineSpace, i::Int) = !$f(P, i)
-    @eval $fnon(P::BSplineSpace) = !$f(P)
+    @eval $fnon(P, i::Integer) = !$f(P, i)
+    @eval $fnon(P) = !$f(P)
 end
-isdegenerate(P::BSplineSpace, i::Integer) = isdegenerate_R(P,i)
-isdegenerate(P::BSplineSpace) = isdegenerate_R(P)
-isnondegenerate(P::BSplineSpace, i::Integer) = isnondegenerate_R(P, i)
-isnondegenerate(P::BSplineSpace) = isnondegenerate_R(P)
+const isdegenerate = isdegenerate_R
+const isnondegenerate = isnondegenerate_R
 
 function _nondegeneratize_R(P::BSplineSpace{p}) where p
     k = copy(KnotVector(knotvector(P)))
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index 0dfad79e1..994fb5839 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -345,11 +345,51 @@ B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
 Assumption:
 * ``P ⊑ P^{\prime}``
 """
+changebasis_I
+
 function changebasis_I(P::AbstractFunctionSpace, P′::AbstractFunctionSpace)
     P ⊑ P′ || throw(DomainError((P,P′),"P ⊑ P′ should be hold."))
     return _changebasis_I(P, P′)
 end
 
+function changebasis_I(P::BSplineSpace, P′::BSplineSpace{p′}) where p′
+    P ⊑ P′ || throw(DomainError((P,P′),"P ⊑ P′ should be hold."))
+
+    k′ = knotvector(P′)
+    l′ = length(k′)
+    i = 1 + p′
+    v = k′[i]
+    i = i + 1
+    while true
+        k′[i] ≠ v && break
+        i = i + 1
+    end
+    degenerated_rank_on_left = i - p′ - 2
+    i = length(k′) - p′
+    v = k′[i]
+    i = i - 1
+    while true
+        k′[i] ≠ v && break
+        i = i - 1
+    end
+    degenerated_rank_on_right = l′ - p′ - i - 1
+    if degenerated_rank_on_left == degenerated_rank_on_right == 0
+        return _changebasis_I(P, P′)
+    else
+        _k′ = view(k′, 1+degenerated_rank_on_left:l′-degenerated_rank_on_right)
+        _P′ = BSplineSpace{p′}(_k′)
+        _A = _changebasis_I(P, _P′)
+        m = _A.m
+        n = _A.n + degenerated_rank_on_left + degenerated_rank_on_right
+        colptr = _A.colptr
+        prepend!(colptr, fill(colptr[begin], degenerated_rank_on_left))
+        append!(colptr, fill(colptr[end], degenerated_rank_on_right))
+        rowval = _A.rowval
+        nzval = _A.nzval
+        return SparseMatrixCSC(m,n,colptr,rowval,nzval)
+    end
+end
+
 function _changebasis_I(P::BSplineSpace{0,T,<:AbstractKnotVector{T}}, P′::BSplineSpace{p′,T′,<:AbstractKnotVector{T′}}) where {p′,T,T′}
     U = StaticArrays.arithmetic_closure(promote_type(T, T′))
     n = dim(P)
@@ -458,7 +498,7 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
     K′ = [k′[j+p′] - k′[j] for j in 1:n′+1]
     K = U[ifelse(k[i+p] ≠ k[i], U(1 / (k[i+p] - k[i])), zero(U)) for i in 1:n+1]
     Aᵖ⁻¹ = _changebasis_I(_lower_I(P), _lower_I(P′))  # (n-1) × (n′-1) matrix
-    n_nonzero = exactdim_I(P′)*(p+2)  # This would be a upper bound of the number of non-zero elements of Aᵖ.
+    n_nonzero = exactdim_I(P′)*(p+1)  # This would be a upper bound of the number of non-zero elements of Aᵖ.
     I = Vector{Int32}(undef, n_nonzero)
     J = Vector{Int32}(undef, n_nonzero)
     V = Vector{U}(undef, n_nonzero)
@@ -505,15 +545,6 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
                           266666666666666661
                           366666666666666661
 
-         1                                    n′
-         |--*-----------------------------*-->|
-            j_begin                       j_end
-           *|---------------------------->|*
-           j_prev                          j_next
-           |                               |
-            |-------------j₊------------->|
-           188888888888888888888888888888881
-
          1                                    n′
          *-----------------------------*----->|
          j_begin                       j_end
@@ -587,9 +618,9 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
         end
         # Rule-0: outside of j_range
         j_next = j_end + 1
-        if j_next == n′+1 || isdegenerate_I(P′, j_next)
+        if j_next == n′+1
             j_mid = j_next - 1
-            if j_prev == 0 || isdegenerate_I(P′, j_prev)
+            if j_prev == 0
                 # Rule-8: right recursion with postprocess
                 # We can't find Aᵖᵢ₁ or Aᵖᵢₙ′ directly (yet!), so we need Δ-shift.
                 # TODO: Find a way to avoid the Δ-shift.
@@ -608,7 +639,7 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
                 V[s-n′:s-1] .+= Δ
                 continue
             end
-        elseif j_prev == 0 || isdegenerate_I(P′, j_prev)
+        elseif j_prev == 0
             j_mid = j_prev
         else
             j_mid = (j_prev + j_next) ÷ 2
diff --git a/src/_DerivativeSpace.jl b/src/_DerivativeSpace.jl
index 313ecbd6c..5c950d2b8 100644
--- a/src/_DerivativeSpace.jl
+++ b/src/_DerivativeSpace.jl
@@ -63,6 +63,14 @@ _lower_I(dP::BSplineDerivativeSpace{r}) where r = BSplineDerivativeSpace{r-1}(_l
 _iszeros_R(P::BSplineDerivativeSpace) = _iszeros_R(bsplinespace(P))
 _iszeros_I(P::BSplineDerivativeSpace) = _iszeros_I(bsplinespace(P))
 
+function isdegenerate_R(dP::BSplineDerivativeSpace, i::Integer)
+    return isdegenerate_R(bsplinespace(dP), i)
+end
+
+function isdegenerate_I(dP::BSplineDerivativeSpace, i::Integer)
+    return isdegenerate_I(bsplinespace(dP), i)
+end
+
 @doc raw"""
     derivative(::BSplineDerivativeSpace{r}) -> BSplineDerivativeSpace{r+1}
     derivative(::BSplineSpace) -> BSplineDerivativeSpace{1}


Test patch:
diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl
index 13986ae53..d279cc1be 100644
--- a/test/test_ChangeBasis.jl
+++ b/test/test_ChangeBasis.jl
@@ -1,36 +1,65 @@
 @testset "ChangeBasis" begin
     ε = 1e-13
     function test_changebasis_R(P,P′)
+        # The test case must hold P ⊆ P′
         @test P ⊆ P′
-        A = @inferred changebasis(P,P′)
+        # Test type stability
+        A = @inferred changebasis_R(P,P′)
+        # Test output type
         @test A isa SparseMatrixCSC
+        # Zeros must not be stored
         @test !any(iszero.(A.nzval))
-        @test A == BasicBSpline._changebasis_R(P,P′) == changebasis_R(P,P′)
+        # Test the size of A
         n = dim(P)
         n′ = dim(P′)
         @test size(A) == (n,n′)
+        # Elements must not be stored for degenerate row/col
+        for j in 1:n′
+            if isdegenerate_R(P′, j)
+                @test !any(Base.isstored.(Ref(A), 1:n, j))
+            end
+        end
+        for i in 1:n
+            if isdegenerate_R(P, i)
+                @test !any(Base.isstored.(Ref(A), i, 1:n′))
+            end
+        end
+        # B_{(i,p,k)} = ∑ⱼ A_{i,j} B_{(j,p′,k′)}
         ts = range(extrema(knotvector(P)+knotvector(P′))..., length=20)
         for t in ts
             @test norm(bsplinebasis.(P,1:n,t) - A*bsplinebasis.(P′,1:n′,t), Inf) < ε
         end
-        @test iszero(view(A, findall(BasicBSpline._iszeros_R(P)), :))
-        @test iszero(view(A, :, findall(BasicBSpline._iszeros_R(P′))))
     end
 
-    function test_changebasis_I(P1, P2; check_zero=true)
-        @test P1 ⊑ P2
-        A = @inferred changebasis_I(P1,P2)
+    function test_changebasis_I(P, P′; check_zero=true)
+        # The test case must hold P ⊑ P′
+        @test P ⊑ P′
+        # Test type stability
+        A = @inferred changebasis_I(P,P′)
+        # Test output type
         @test A isa SparseMatrixCSC
-        if check_zero
-            @test !any(iszero.(A.nzval))
+        # Zeros must not be stored
+        check_zero && @test !any(iszero.(A.nzval))
+        # Test the size of A
+        n = dim(P)
+        n′ = dim(P′)
+        @test size(A) == (n,n′)
+        # Elements must not be stored for degenerate row/col
+        for j in 1:n′
+            if isdegenerate_I(P′, j)
+                @test !any(Base.isstored.(Ref(A), 1:n, j))
+            end
         end
-        n1 = dim(P1)
-        n2 = dim(P2)
-        @test size(A) == (n1,n2)
-        d = domain(P1)
+        for i in 1:n
+            if isdegenerate_I(P, i)
+                @test !any(Base.isstored.(Ref(A), i, 1:n′))
+            end
+        end
+        # B_{(i,p,k)} = ∑ⱼ A_{i,j} B_{(j,p′,k′)}
+        d = domain(P)
         ts = range(extrema(d)..., length=21)[2:end-1]
         for t in ts
-            @test norm(bsplinebasis.(P1,1:n1,t) - A*bsplinebasis.(P2,1:n2,t), Inf) < ε
+            @test norm(bsplinebasis.(P,1:n,t) - A*bsplinebasis.(P′,1:n′,t), Inf) < ε
         end
     end
 
@@ -65,7 +94,7 @@
         A13 = changebasis(P1, P3)
         A36 = changebasis(P3, P6)
         A16 = changebasis(P1, P6)
-        A16 ≈ A13 * A36
+        @test A16 ≈ A13 * A36
 
         test_changebasis_R(P3, P6)
         test_changebasis_R(P6, P7)
@@ -73,7 +102,7 @@
         A36 = changebasis(P3, P6)
         A67 = changebasis(P6, P7)
         A37 = changebasis(P3, P7)
-        A37 ≈ A36 * A67
+        @test A37 ≈ A36 * A67
 
         @test P2 ⊈ P3
 
@@ -123,13 +152,31 @@
         _P8 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector(-[1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 7, 7]))
         test_changebasis_I(_P7, _P8)
 
+        # (2,2)-element is stored, but it is zero.
         Q1 = BSplineSpace{1, Int64, KnotVector{Int64}}(KnotVector([2, 2, 4, 4, 6, 6]))
         Q2 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 2, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7]))
         test_changebasis_I(Q1, Q2; check_zero=false)
 
         Q3 = BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7]))
         Q4 = BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6]))
-        test_changebasis_I(Q3, Q4; check_zero=false)
+        test_changebasis_I(Q3, Q4)
+    end
+
+    @testset "different changebasis_R and changebasis_I" begin
+        P1 = BSplineSpace{3}(knotvector"1111 132")
+        P2 = BSplineSpace{3}(knotvector"11121132")
+        @test isdegenerate_I(P1)
+        @test isdegenerate_I(P2)
+        @test isnondegenerate_R(P1)
+        @test isnondegenerate_R(P2)
+        @test P1 ⊆ P2
+        @test P1 ⊑ P2
+
+        test_changebasis_I(P1, P2)
+        test_changebasis_R(P1, P2)
+        A_R = changebasis_R(P1, P2)
+        A_I = changebasis_I(P1, P2)
+        @test A_R ≠ A_I
     end
 
     @testset "changebasis_sim" begin
diff --git a/test/test_Derivative.jl b/test/test_Derivative.jl
index 685810e52..de84fc6cf 100644
--- a/test/test_Derivative.jl
+++ b/test/test_Derivative.jl
@@ -33,6 +33,25 @@
         @test_throws MethodError convert(BSplineDerivativeSpace{2},dP2)
     end
 
+    @testset "degenerate" begin
+        P1 = BSplineSpace{2}(KnotVector(1:8))
+        P2 = BSplineSpace{2}(knotvector"111111151111")
+        dP1 = BSplineDerivativeSpace{1}(P1)
+        dP2 = BSplineDerivativeSpace{1}(P2)
+        @test isdegenerate_R(P1) == isdegenerate_R(dP1)
+        @test isdegenerate_R(P2) == isdegenerate_R(dP2)
+        @test isdegenerate_I(P1) == isdegenerate_I(dP1)
+        @test isdegenerate_I(P2) == isdegenerate_I(dP2)
+        for i in 1:dim(P1)
+            @test isdegenerate_R(P1, i) == isdegenerate_R(dP1, i)
+            @test isdegenerate_I(P1, i) == isdegenerate_I(dP1, i)
+        end
+        for i in 1:dim(P2)
+            @test isdegenerate_R(P2, i) == isdegenerate_R(dP2, i)
+            @test isdegenerate_I(P2, i) == isdegenerate_I(dP2, i)
+        end
+    end
+
     # 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