Skip to content

Refactor changebasis_I #399

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
hyrodium merged 5 commits into from
Feb 8, 2026
Merged

Refactor changebasis_I #399

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
be7662d
add `_find_j_range_I_new` etc.
hyrodium Feb 7, 2026
7089915
add more test
hyrodium Feb 7, 2026
a517b36
comment out old methods
hyrodium Feb 7, 2026
0741ebe
update tests
hyrodium Feb 8, 2026
4b7a2c1
update edge case
hyrodium Feb 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
137 changes: 107 additions & 30 deletions src/_ChangeBasis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -489,39 +489,116 @@ 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′}
# Old methoods for changebasis_I.
# Keep these functions for the following edge case.
# 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)
# 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)
end
k = knotvector(P)
k′ = knotvector(P′)
n′ = dim(P′)
_, j_end_prev = extrema(j_range)

_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
# 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

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
Pi = BSplineSpace{p}(view(knotvector(P), i:i+p+1))
if Pi ⊆ P′
return _find_j_range_R(P, P′, i, j_range)
# 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
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
end

Expand Down Expand Up @@ -576,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
Expand Down Expand Up @@ -647,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
Expand Down
75 changes: 54 additions & 21 deletions test/test_ChangeBasis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,34 +132,67 @@
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)

# 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)

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)
# Edge case
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

@testset "different changebasis_R and changebasis_I" begin
Expand Down
Loading