From be7662d726fe2b6d991ab7724ccf02819577d177 Mon Sep 17 00:00:00 2001 From: Horikawa Yuto Date: Sat, 7 Feb 2026 14:36:48 +0900 Subject: [PATCH 1/5] add `_find_j_range_I_new` etc. --- src/_ChangeBasis.jl | 80 ++++++++++++++++++++++++++++++++++++++-- test/test_ChangeBasis.jl | 70 ++++++++++++++++++++++++----------- 2 files changed, 125 insertions(+), 25 deletions(-) diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl index 5e64f2cbf..539cf2b4d 100644 --- a/src/_ChangeBasis.jl +++ b/src/_ChangeBasis.jl @@ -512,17 +512,89 @@ function __changebasis_I_old(P1::BSplineSpace{p,T}, P2::BSplineSpace{p′,T′}) return _A end +function _find_j_range_I_new(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range) where {p, p′} + Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1)) + if Pi ⊆ P′ + return _find_j_range_R(P, P′, i, j_range) + end + k = knotvector(P) + k′ = knotvector(P′) + n′ = dim(P′) + _, j_end_prev = extrema(j_range) + + # Determine j_end + b = k′[end-p′] # right endpoint of domain(P′) + if k[i+p+1] > b + # Right side is strictly beyond the domain boundary → j_end = n′ + j_end = n′ + else + # R-style multiplicity counting (works when k[i+p+1] ≤ b) + j_end = n′ # fallback + m = p′ - p + t_end = k[i+p+1] + for ii in 0:(p+1) + if k[i+p+1-ii] == t_end + m += 1 + else + break + end + end + for j in (j_end_prev-p′-1):n′ + if t_end == k′[j+p′+1] + m -= 1 + end + if m == 0 + j_end = j + break + end + end + end + + # Determine j_begin + a = k′[1+p′] # left endpoint of domain(P′) + if k[i] < a + # Left side is strictly before the domain boundary → j_begin = 1 + j_begin = 1 + else + # R-style multiplicity counting (works when k[i] ≥ a) + j_begin = 1 # fallback + m = p′ - p + t_begin = k[i] + for ii in 0:(p+1) + if k[i+ii] == t_begin + m += 1 + else + break + end + end + for j in reverse(1:(j_end+p′+1)) + if k′[j] == t_begin + m -= 1 + end + if m == 0 + j_begin = j + break + end + end + end + + return j_begin:j_end +end + function _find_j_range_I(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range, Aᵖ_old) where {p, p′} # TODO: remove `Aᵖ_old` argument # TODO: fix performance https://github.com/hyrodium/BasicBSpline.jl/pull/323#issuecomment-1723216566 # TODO: remove threshold such as 1e-14 + j_range_new = _find_j_range_I_new(P, P′, i, j_range) Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1)) if Pi ⊆ P′ - return _find_j_range_R(P, P′, i, j_range) + return j_range_new end - j_begin = findfirst(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) - j_end = findlast(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) - return j_begin:j_end + j_begin_old = findfirst(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) + j_end_old = findlast(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) + j_range_old = j_begin_old:j_end_old + @assert j_range_new == j_range_old "Mismatch in _find_j_range_I for P=$P, P′=$P′, i=$i: new=$j_range_new, old=$j_range_old" + return j_range_new end function _ΔAᵖ_I(Aᵖ⁻¹::AbstractMatrix, K::AbstractVector, K′::AbstractVector, i::Integer, j::Integer) diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl index d279cc1be..d60a41a24 100644 --- a/test/test_ChangeBasis.jl +++ b/test/test_ChangeBasis.jl @@ -132,34 +132,62 @@ 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) + P1 = BSplineSpace{p1-1}(knotvector(P1)[2:end-1]) + P2 = BSplineSpace{p2-1}(knotvector(P2)[2:end-1]) + @test P1 ⊑ P2 + @test P1 ⊈ P2 + test_changebasis_I(P1, P2) # 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) + P1 = BSplineSpace{2}(knotvector" 21 3") + P2 = BSplineSpace{3}(knotvector"212 4") + test_changebasis_I(P1, P2) - 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) + P1 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 6])) + P2 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 7, 7])) + test_changebasis_I(P1, P2) - _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) + P1 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector(-[1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 6])) + P2 = BSplineSpace{3, Int64, KnotVector{Int64}}(KnotVector(-[1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 7, 7])) + test_changebasis_I(P1, P2) # (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) + P1 = BSplineSpace{1, Int64, KnotVector{Int64}}(KnotVector([2, 2, 4, 4, 6, 6])) + P2 = 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(P1, P2; check_zero=false) + + P1 = BSplineSpace{4, Int64, KnotVector{Int64}}(KnotVector([1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 7])) + P2 = BSplineSpace{5, Int64, KnotVector{Int64}}(KnotVector([1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 4, 4, 4, 5, 6])) + test_changebasis_I(P1, P2) + + # Diverse degree combinations (P ⊑ P′, P ⊄ P′) + P1 = BSplineSpace{2}(KnotVector([0, 0, 1, 2, 3, 3])) + P2 = BSplineSpace{3}(KnotVector([-1, -1, 0, 1, 2, 3, 4, 4])) + test_changebasis_I(P1, P2) + + # One-sided boundary: left only + P1 = BSplineSpace{2}(KnotVector([1, 1, 2, 3, 4, 5])) + P2 = BSplineSpace{3}(KnotVector([0, 1, 1, 2, 3, 4, 5, 6])) + test_changebasis_I(P1, P2) + + # One-sided boundary: right only + P1 = BSplineSpace{2}(KnotVector([1, 2, 3, 4, 5, 5])) + P2 = BSplineSpace{3}(KnotVector([0, 1, 2, 3, 4, 5, 5, 6])) + test_changebasis_I(P1, P2) - 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) + # Both-sides boundary + P1 = BSplineSpace{2}(KnotVector([1, 1, 2, 3, 4, 4])) + P2 = BSplineSpace{3}(KnotVector([0, 0, 1, 2, 3, 4, 5, 5])) + test_changebasis_I(P1, P2) + + # Higher degree combinations + P1 = BSplineSpace{3}(KnotVector([0, 0, 0, 1, 2, 3, 3, 3])) + P2 = BSplineSpace{4}(KnotVector([-1, -1, -1, 0, 1, 2, 3, 4, 4, 4])) + test_changebasis_I(P1, P2) + + P1 = BSplineSpace{2}(KnotVector([0, 0, 1, 2, 3, 3])) + P2 = BSplineSpace{4}(KnotVector([-1, -1, -1, 0, 1, 2, 3, 4, 4, 4])) + test_changebasis_I(P1, P2) end @testset "different changebasis_R and changebasis_I" begin From 7089915cbca50216511c26c7763372cd2f3cb93a Mon Sep 17 00:00:00 2001 From: Horikawa Yuto Date: Sat, 7 Feb 2026 22:42:41 +0900 Subject: [PATCH 2/5] add more test --- test/test_ChangeBasis.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl index d60a41a24..1cafa81c3 100644 --- a/test/test_ChangeBasis.jl +++ b/test/test_ChangeBasis.jl @@ -188,6 +188,11 @@ P1 = BSplineSpace{2}(KnotVector([0, 0, 1, 2, 3, 3])) P2 = BSplineSpace{4}(KnotVector([-1, -1, -1, 0, 1, 2, 3, 4, 4, 4])) test_changebasis_I(P1, P2) + + # Edge case + P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2, 3, 3, 3, 8, 8, 9, 11, 11, 11])) + P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 8, 8, 8, 8, 9, 9, 9, 11, 11, 12, 12, 13])) + test_changebasis_I(P1, P2) end @testset "different changebasis_R and changebasis_I" begin From a517b369a9f369660cbc44ba8a616455e7ecce47 Mon Sep 17 00:00:00 2001 From: Horikawa Yuto Date: Sun, 8 Feb 2026 08:59:49 +0900 Subject: [PATCH 3/5] comment out old methods --- src/_ChangeBasis.jl | 89 ++++++++++++++++++++++++--------------------- 1 file changed, 47 insertions(+), 42 deletions(-) diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl index 539cf2b4d..cf6122e1d 100644 --- a/src/_ChangeBasis.jl +++ b/src/_ChangeBasis.jl @@ -489,30 +489,51 @@ function __changebasis_I(P::BSplineSpace{0,T,<:AbstractKnotVector{T}}, P′::BSp return A⁰ end -function _changebasis_I_old(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′}) where {p,p′,T,T′} - k = knotvector(P) - k′ = knotvector(P′) - - _P = BSplineSpace{p}(k[1+p:end-p] + p * KnotVector{T}([k[1+p], k[end-p]])) - _P′ = BSplineSpace{p′}(k′[1+p′:end-p′] + p′ * KnotVector{T′}([k′[1+p′], k′[end-p′]])) - _A = _changebasis_R(_P, _P′) - Asim = _changebasis_sim(P, _P) - Asim′ = _changebasis_sim(_P′, P′) - A = Asim * _A * Asim′ - - return A -end - -function __changebasis_I_old(P1::BSplineSpace{p,T}, P2::BSplineSpace{p′,T′}) where {p,p′,T,T′} - U = StaticArrays.arithmetic_closure(promote_type(T, T′)) - _P1 = _nondegeneratize_I(P1) - _P2 = _nondegeneratize_I(P2) - _A = sparse(Int32[], Int32[], U[], dim(P1), dim(P2)) - _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_new(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range) where {p, p′} +# Old methoods for changebasis_I. +# Keep these functions for the following edge case. +# P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2, 3, 3, 3, 8, 8, 9, 11, 11, 11])) +# P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 8, 8, 8, 8, 9, 9, 9, 11, 11, 12, 12, 13])) +# +# function _changebasis_I_old(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′}) where {p,p′,T,T′} +# k = knotvector(P) +# k′ = knotvector(P′) +# +# _P = BSplineSpace{p}(k[1+p:end-p] + p * KnotVector{T}([k[1+p], k[end-p]])) +# _P′ = BSplineSpace{p′}(k′[1+p′:end-p′] + p′ * KnotVector{T′}([k′[1+p′], k′[end-p′]])) +# _A = _changebasis_R(_P, _P′) +# Asim = _changebasis_sim(P, _P) +# Asim′ = _changebasis_sim(_P′, P′) +# A = Asim * _A * Asim′ +# +# return A +# end +# +# function __changebasis_I_old(P1::BSplineSpace{p,T}, P2::BSplineSpace{p′,T′}) where {p,p′,T,T′} +# U = StaticArrays.arithmetic_closure(promote_type(T, T′)) +# _P1 = _nondegeneratize_I(P1) +# _P2 = _nondegeneratize_I(P2) +# _A = sparse(Int32[], Int32[], U[], dim(P1), dim(P2)) +# _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_old(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range, Aᵖ_old) where {p, p′} +# # TODO: remove `Aᵖ_old` argument +# # TODO: fix performance https://github.com/hyrodium/BasicBSpline.jl/pull/323#issuecomment-1723216566 +# # TODO: remove threshold such as 1e-14 +# j_range_new = _find_j_range_I_new(P, P′, i, j_range) +# Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1)) +# if Pi ⊆ P′ +# return j_range_new +# end +# j_begin_old = findfirst(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) +# j_end_old = findlast(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) +# j_range_old = j_begin_old:j_end_old +# @assert j_range_new == j_range_old "Mismatch in _find_j_range_I for P=$P, P′=$P′, i=$i: new=$j_range_new, old=$j_range_old" +# return j_range_new +# end + +function _find_j_range_I(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range) where {p, p′} Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1)) if Pi ⊆ P′ return _find_j_range_R(P, P′, i, j_range) @@ -581,22 +602,6 @@ function _find_j_range_I_new(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_ return j_begin:j_end end -function _find_j_range_I(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_range, Aᵖ_old) where {p, p′} - # TODO: remove `Aᵖ_old` argument - # TODO: fix performance https://github.com/hyrodium/BasicBSpline.jl/pull/323#issuecomment-1723216566 - # TODO: remove threshold such as 1e-14 - j_range_new = _find_j_range_I_new(P, P′, i, j_range) - Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1)) - if Pi ⊆ P′ - return j_range_new - end - j_begin_old = findfirst(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) - j_end_old = findlast(e->abs(e)>1e-14, view(Aᵖ_old, i, :)) - j_range_old = j_begin_old:j_end_old - @assert j_range_new == j_range_old "Mismatch in _find_j_range_I for P=$P, P′=$P′, i=$i: new=$j_range_new, old=$j_range_old" - return j_range_new -end - function _ΔAᵖ_I(Aᵖ⁻¹::AbstractMatrix, K::AbstractVector, K′::AbstractVector, i::Integer, j::Integer) n = length(K)-1 if i == 1 @@ -648,7 +653,7 @@ function __changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSp V = Vector{U}(undef, n_nonzero) s = 1 j_range = 1:1 # j_begin:j_end - Aᵖ_old = __changebasis_I_old(P, P′) + # Aᵖ_old = __changebasis_I_old(P, P′) for i in 1:n # Skip for degenerated basis isdegenerate_I(P,i) && continue @@ -719,7 +724,7 @@ function __changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSp =# # Precalculate the range of j - j_range = _find_j_range_I(P, P′, i, j_range, Aᵖ_old) + j_range = _find_j_range_I(P, P′, i, j_range) j_begin, j_end = extrema(j_range) # Rule-0: outside of j_range From 0741ebefb9e9c21641349fc670e129c7ff5ccbcd Mon Sep 17 00:00:00 2001 From: Horikawa Yuto Date: Sun, 8 Feb 2026 09:01:11 +0900 Subject: [PATCH 4/5] update tests --- test/test_ChangeBasis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl index 1cafa81c3..df7b1f117 100644 --- a/test/test_ChangeBasis.jl +++ b/test/test_ChangeBasis.jl @@ -192,7 +192,7 @@ # Edge case P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2, 3, 3, 3, 8, 8, 9, 11, 11, 11])) P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 8, 8, 8, 8, 9, 9, 9, 11, 11, 12, 12, 13])) - test_changebasis_I(P1, P2) + test_changebasis_I(P1, P2; check_zero=false) end @testset "different changebasis_R and changebasis_I" begin From 4b7a2c1fd049c24f22877d0a85139a79ab8160f1 Mon Sep 17 00:00:00 2001 From: Horikawa Yuto Date: Sun, 8 Feb 2026 17:19:49 +0900 Subject: [PATCH 5/5] update edge case --- src/_ChangeBasis.jl | 4 ++-- test/test_ChangeBasis.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl index cf6122e1d..a414defb3 100644 --- a/src/_ChangeBasis.jl +++ b/src/_ChangeBasis.jl @@ -491,8 +491,8 @@ end # Old methoods for changebasis_I. # Keep these functions for the following edge case. -# P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2, 3, 3, 3, 8, 8, 9, 11, 11, 11])) -# P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 8, 8, 8, 8, 9, 9, 9, 11, 11, 12, 12, 13])) +# P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2])) +# P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2])) # # function _changebasis_I_old(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′}) where {p,p′,T,T′} # k = knotvector(P) diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl index df7b1f117..fbabf7b5e 100644 --- a/test/test_ChangeBasis.jl +++ b/test/test_ChangeBasis.jl @@ -190,8 +190,8 @@ test_changebasis_I(P1, P2) # Edge case - P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2, 3, 3, 3, 8, 8, 9, 11, 11, 11])) - P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 8, 8, 8, 8, 9, 9, 9, 11, 11, 12, 12, 13])) + P1 = BSplineSpace{2}(KnotVector([-3, -3, -2, 1, 2, 2, 2])) + P2 = BSplineSpace{4}(KnotVector([-5, -4, -2, -2, -2, 1, 1, 1, 2, 2, 2, 2, 2])) test_changebasis_I(P1, P2; check_zero=false) end