Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
`changebasis_I` sometimes returns a matrix with `NaN`
```julia
julia> using BasicBSpline

julia> Pd1 = BSplineSpace{3}(knotvector"44 3  2 121 ")
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 7, 7, 9, 10, 10, 11]))

julia> Pd2 = BSplineSpace{4}(knotvector"55 62 3 2121")
BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 5, 5, 7, 7, 7, 9, 9, 10, 11, 11, 12]))

julia> changebasis_I(Pd1, Pd2)
13×22 SparseArrays.SparseMatrixCSC{Float64, Int32} with 44 stored entries:
⎡⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠈⠳⣄⢀⣀⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⠿⠶⣤⣤⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⎦
```

Pull request:
This PR fixes #325.

**Before this PR**
```julia
julia> using BasicBSpline

julia> Pd1 = BSplineSpace{2}(knotvector" 21 3")
BSplineSpace{2, Int64, KnotVector{Int64}}(KnotVector([2, 2, 3, 5, 5, 5]))

julia> Pd2 = BSplineSpace{3}(knotvector"212 4")
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 3, 3, 5, 5, 5, 5]))

julia> BasicBSpline.changebasis_I(Pd1, Pd2)
ERROR: MethodError: no method matching isless(::Nothing, ::Nothing)

Closest candidates are:
  isless(::Any, ::Missing)
   @ Base missing.jl:88
  isless(::Missing, ::Any)
   @ Base missing.jl:87

Stacktrace:
 [1] <(x::Nothing, y::Nothing)
   @ Base ./operators.jl:343
 [2] <=(x::Nothing, y::Nothing)
   @ Base ./operators.jl:392
 [3] >=(x::Nothing, y::Nothing)
   @ Base ./operators.jl:416
 [4] (::Colon)(start::Nothing, stop::Nothing)
   @ Base ./range.jl:7
 [5] _find_j_range_I(P::BSplineSpace{2, Int64, KnotVector{Int64}}, P′::BSplineSpace{3, Int64, KnotVector{Int64}}, i::Int64, j_begin::Int64, j_end::Int64)
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:405
 [6] _changebasis_I(P::BSplineSpace{2, Int64, KnotVector{Int64}}, P′::BSplineSpace{3, Int64, KnotVector{Int64}})
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:527
 [7] changebasis_I(P::BSplineSpace{2, Int64, KnotVector{Int64}}, P′::BSplineSpace{3, Int64, KnotVector{Int64}})
   @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:350
 [8] top-level scope
   @ REPL[4]:1
```

**After this PR**
```julia
julia> using BasicBSpline

julia> Pd1 = BSplineSpace{2}(knotvector" 21 3")
BSplineSpace{2, Int64, KnotVector{Int64}}(KnotVector([2, 2, 3, 5, 5, 5]))

julia> Pd2 = BSplineSpace{3}(knotvector"212 4")
BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 3, 3, 5, 5, 5, 5]))

julia> BasicBSpline.changebasis_I(Pd1, Pd2)
3×5 SparseArrays.SparseMatrixCSC{Float64, Int32} with 7 stored entries:
  ⋅   0.888889  0.222222   ⋅         ⋅ 
  ⋅   0.111111  0.777778  0.666667   ⋅ 
  ⋅    ⋅         ⋅        0.333333  1.0
```

Patch:
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index a46eefcd1..e44c9c775 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -130,8 +130,8 @@ Modified version.
 {B}_{(i,0,k)}(t)
 &=
 \begin{cases}
-    &1\quad (k_{i} \le t<k_{i+1})\\
-    &1\quad (k_{i} < t = k_{i+1}=k_{l})\\
+    &1\quad (k_{i} \le t < k_{i+1}) \\
+    &1\quad (k_{i} < t = k_{i+1}=k_{l}) \\
     &0\quad (\text{otherwise})
 \end{cases}
 \end{aligned}
@@ -189,8 +189,9 @@ Modified version (2).
 {B}_{(i,0,k)}(t)
 &=
 \begin{cases}
-    &1\quad (k_{i} \le t<k_{i+1})\\
-    &1\quad (k_{i} < t = k_{i+1}=k_{l})\\
+    &1\quad (k_{1+p} \le k_{i} < t \le k_{i+1}) \\
+    &1\quad (t = k_{1+p} = k_{i} < k_{i+1}) \\
+    &1\quad (k_{i} \le t < k_{i+1} \le k_{1+p}) \\
     &0\quad (\text{otherwise})
 \end{cases}
 \end{aligned}
@@ -219,7 +220,7 @@ julia> BasicBSpline.bsplinebasis₋₀I.(P,1:5,(1:6)')
     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])) || (v[1]==$(ks[i])==t<$(ks[i+1])))) for i in 1:n]...)
+    A_r(n) = Expr(:tuple, [:($U((v[1+p]≤$(ks[i])<t≤$(ks[i+1])) || (t==v[1+p]==$(ks[i])<$(ks[i+1])) || ($(ks[i])≤t<$(ks[i+1])≤v[1+p]))) 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[]
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index 4491834b0..eea0c61b7 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -273,6 +273,32 @@ isdegenerate(P::BSplineSpace) = isdegenerate_R(P)
 isnondegenerate(P::BSplineSpace, i::Integer) = isnondegenerate_R(P, i)
 isnondegenerate(P::BSplineSpace) = isnondegenerate_R(P)
 
+function _degeneratize_R(P::BSplineSpace{p}) where p
+    k = copy(KnotVector(knotvector(P)))
+    I = Int[]
+    for i in 1:dim(P)
+        isdegenerate_R(P,i) && push!(I, i)
+    end
+    deleteat!(BasicBSpline._vec(k), I)
+    return BSplineSpace{p}(k)
+end
+function _degeneratize_I(P::BSplineSpace{p}) where p
+    k = copy(KnotVector(knotvector(P)))
+    l = length(k)
+    I = Int[]
+    for i in 1:dim(P)
+        if isdegenerate_I(P,i)
+            if k[l-p] < k[i+p+1]
+                push!(I, i+p+1)
+            else
+                push!(I, i)
+            end
+        end
+    end
+    deleteat!(BasicBSpline._vec(k), I)
+    return BSplineSpace{p}(k)
+end
+
 """
 Exact dimension of a B-spline space.
 
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index fc0e75072..94c2517dc 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -395,11 +395,20 @@ function _changebasis_I_old(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′})
     return A
 end
 
+function __changebasis_I_old(P1, P2)
+    _P1 = _degeneratize_I(P1)
+    _P2 = _degeneratize_I(P2)
+    A = _changebasis_I_old(P1,P2)
+    _A = zero(A)
+    _A[isnondegenerate_I.(P1,1:dim(P1)), isnondegenerate_I.(P2,1:dim(P2))] = _changebasis_I_old(_P1,_P2)
+    return _A
+end
+
 function _find_j_range_I(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_begin, j_end) where {p, p′}
     # TODO: avoid `_changebasis_I_old`
     # TODO: fix performance https://github.com/hyrodium/BasicBSpline.jl/pull/323#issuecomment-1723216566
     # TODO: remove threshold such as 1e-14
-    Aᵖ_old = _changebasis_I_old(P,P′)
+    Aᵖ_old = __changebasis_I_old(P,P′)
     j_begin = findfirst(e->abs(e)>1e-14, Aᵖ_old[i, :])
     j_end = findlast(e->abs(e)>1e-14, Aᵖ_old[i, :])
     return j_begin:j_end
@@ -450,10 +459,11 @@ 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+1)  # This is a upper bound of the number of non-zero elements of Aᵖ (rough estimation).
+    n_nonzero = exactdim_I(P′)*(p+2)  # 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)
+    # R = fill(-1, n_nonzero)
     s = 1
     j_begin = 1
     j_end = 1
@@ -482,6 +492,7 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
            |       |        |
             |--j₊->||<-j₋--|
            066666666777777773
+           177777777777777773
 
          1                                    n′
          |--*-----------------------------*-->|
@@ -492,6 +503,17 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
                            |--j₊->||<-j₋--|
                           266666666777777770
                           366666666777777770
+                          266666666666666661
+                          366666666666666661
+
+         1                                    n′
+         |--*-----------------------------*-->|
+            j_begin                       j_end
+           *|---------------------------->|*
+           j_prev                          j_next
+           |                               |
+            |-------------j₊------------->|
+           188888888888888888888888888888881
 
          1                                    n′
          *-----------------------------*----->|
@@ -537,24 +559,28 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
             if k′[j_next] == k′[j_next+p′]
                 I[s], J[s] = i, j_next
                 V[s] = Aᵖᵢⱼ_prev = bsplinebasis₊₀(P,i,k′[j_next+1])
+                # R[s] = 2
                 s += 1
                 j_prev = j_next
             # Rule-3: left limit (or both limit)
             elseif k′[j_next+1] == k′[j_next+p′]
-                j_mid = (j_prev == 0 ? j_prev : (j_prev + j_next) ÷ 2)
+                j_mid = (j_prev == 0 || isdegenerate_I(P′, j_prev) ? j_prev : (j_prev + j_next) ÷ 2)
                 # Rule-6: right recursion
                 for j₊ in (j_prev+1):j_mid
                     I[s], J[s] = i, j₊
                     V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_prev + p * _ΔAᵖ_I(Aᵖ⁻¹,K,K′,i,j₊) / p′
+                    # R[s] = 6
                     s += 1
                 end
                 I[s], J[s] = i, j_next
                 V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_next = bsplinebasis₋₀I(P,i,k′[j_next+1])
+                # R[s] = 3
                 s += 1
                 # Rule-7: left recursion
                 for j₋ in reverse((j_mid+1):(j_next-1))
                     I[s], J[s] = i, j₋
                     V[s] = Aᵖᵢⱼ_next = Aᵖᵢⱼ_next - p * _ΔAᵖ_I(Aᵖ⁻¹,K,K′,i,j₋+1) / p′
+                    # R[s] = 7
                     s += 1
                 end
                 j_prev = j_next
@@ -562,26 +588,28 @@ 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
+        if j_next == n′+1 || isdegenerate_I(P′, j_next)
             j_mid = j_next - 1
-            if j_prev == 0
+            if j_prev == 0 || isdegenerate_I(P′, j_prev)
                 # 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.
                 I[s], J[s] = i, 1
                 V[s] = Aᵖᵢⱼ_prev = zero(U)
+                # R[s] = 8
                 s += 1
                 for j₊ in 2:n′
                     I[s], J[s] = i, j₊
                     V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_prev + p * _ΔAᵖ_I(Aᵖ⁻¹,K,K′,i,j₊) / p′
+                    # R[s] = 8
                     s += 1
                 end
-                t_mid = (maximum(domain(P))+minimum(domain(P))) / (2*one(U))
+                t_mid = (maximum(domain(P))+minimum(domain(P))) / 2one(U)
                 Δ = bsplinebasis(P, i, t_mid) - dot(view(V, s-n′:s-1), bsplinebasis.(P′, 1:n′, t_mid))
                 V[s-n′:s-1] .+= Δ
                 continue
             end
-        elseif j_prev == 0
+        elseif j_prev == 0 || isdegenerate_I(P′, j_prev)
             j_mid = j_prev
         else
             j_mid = (j_prev + j_next) ÷ 2
@@ -591,17 +619,21 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
         for j₊ in (j_prev+1):j_mid
             I[s], J[s] = i, j₊
             V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_prev + p * _ΔAᵖ_I(Aᵖ⁻¹,K,K′,i,j₊) / p′
+            # R[s] = 6
             s += 1
         end
         # Rule-7: left recursion
         for j₋ in reverse((j_mid+1):(j_next-1))
             I[s], J[s] = i, j₋
             V[s] = Aᵖᵢⱼ_next = Aᵖᵢⱼ_next - p * _ΔAᵖ_I(Aᵖ⁻¹,K,K′,i,j₋+1) / p′
+            # R[s] = 7
             s += 1
         end
     end
 
     Aᵖ = sparse(view(I,1:s-1), view(J,1:s-1), view(V,1:s-1), n, n′)
+    # Rᵖ = sparse(view(I,1:s-1), view(J,1:s-1), view(R,1:s-1), n, n′)
+    # display(Rᵖ)
     return Aᵖ
 end
 


Test patch:
diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl
index 3e221b8a6..13986ae53 100644
--- a/test/test_ChangeBasis.jl
+++ b/test/test_ChangeBasis.jl
@@ -17,11 +17,13 @@
         @test iszero(view(A, :, findall(BasicBSpline._iszeros_R(P′))))
     end
 
-    function test_changebasis_I(P1,P2)
+    function test_changebasis_I(P1, P2; check_zero=true)
         @test P1 ⊑ P2
         A = @inferred changebasis_I(P1,P2)
         @test A isa SparseMatrixCSC
-        @test !any(iszero.(A.nzval))
+        if check_zero
+            @test !any(iszero.(A.nzval))
+        end
         n1 = dim(P1)
         n2 = dim(P2)
         @test size(A) == (n1,n2)
@@ -86,26 +88,48 @@
     end
 
     @testset "sqsubseteq (I)" begin
-        p4 = 1
-        p5 = 2
-        P4 = BSplineSpace{p4}(KnotVector([1, 2, 3, 4, 5]))
-        P5 = BSplineSpace{p5}(KnotVector([-1, 0.3, 2, 3, 3, 4, 5.2, 6]))
-        @test P4 ⊑ P4
-        @test P4 ⊑ P5
-        @test P5 ⊒ P4
-        @test P5 ⋢ P4
-        @test P4 ⋣ P5
-
-        test_changebasis_I(P4, P5)
-        test_changebasis_I(P4, P4)
-        test_changebasis_I(P5, P5)
-
-        P6 = BSplineSpace{p4-1}(knotvector(P4)[2:end-1])
-        P7 = BSplineSpace{p5-1}(knotvector(P5)[2:end-1])
-        @test P6 ⊑ P7
-        @test P6 ⊈ P7
-
-        test_changebasis_I(P6, P7)
+        p1 = 1
+        p2 = 2
+        P1 = BSplineSpace{p1}(KnotVector([1, 2, 3, 4, 5]))
+        P2 = BSplineSpace{p2}(KnotVector([-1, 0.3, 2, 3, 3, 4, 5.2, 6]))
+        @test P1 ⊑ P1
+        @test P1 ⊑ P2
+        @test P2 ⊒ P1
+        @test P2 ⋢ P1
+        @test P1 ⋣ P2
+        @test P1 ⊈ P2
+
+        test_changebasis_I(P1, P2)
+        test_changebasis_I(P1, P1)
+        test_changebasis_I(P2, P2)
+
+        P3 = BSplineSpace{p1-1}(knotvector(P1)[2:end-1])
+        P4 = BSplineSpace{p2-1}(knotvector(P2)[2:end-1])
+        @test P3 ⊑ P4
+        @test P3 ⊈ P4
+
+        test_changebasis_I(P3, P4)
+
+        # https://github.com/hyrodium/BasicBSpline.jl/issues/325
+        P5 = BSplineSpace{2}(knotvector" 21 3")
+        P6 = BSplineSpace{3}(knotvector"212 4")
+        test_changebasis_I(P5, P6)
+
+        P7 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 6]))
+        P8 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 7, 7]))
+        test_changebasis_I(P7, P8)
+
+        _P7 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector(-[1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 6]))
+        _P8 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector(-[1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 7, 7]))
+        test_changebasis_I(_P7, _P8)
+
+        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)
     end
 
     @testset "changebasis_sim" begin