Skip to content

Commit

Permalink
minor fixes around changebasis_I (#338)
Browse files Browse the repository at this point in the history
* update internal function

* update comments

* update `test_changebasis_I`
  • Loading branch information
hyrodium authored Nov 10, 2023
1 parent 6ed8f79 commit 7b5ba2b
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 20 deletions.
20 changes: 10 additions & 10 deletions src/_ChangeBasis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ function _changebasis_R(P::BSplineSpace{0,T,KnotVector{T}}, P′::BSplineSpace{p
return A⁰
end

function _find_j_begin_end_R(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_begin, j_end):: Tuple{Int, Int} where {p, p′}
function _find_j_range_R(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_begin, j_end) where {p, p′}
k = knotvector(P)
k′ = knotvector(P′)
n′ = dim(P′)
Expand Down Expand Up @@ -100,7 +100,7 @@ function _find_j_begin_end_R(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_
end
end

return j_begin, j_end
return j_begin:j_end
end

function _ΔAᵖ_R(Aᵖ⁻¹::AbstractMatrix, K::AbstractVector, K′::AbstractVector, i::Integer, j::Integer)
Expand All @@ -124,7 +124,7 @@ function _changebasis_R(P::BSplineSpace{p,T,KnotVector{T}}, P′::BSplineSpace{p
=== j-rules for fixed i ===
- Rule-0: B₍ᵢ,ₚ,ₖ₎ ∈ span₍ᵤ₌ₛ…ₜ₎(B₍ᵤ,ₚ′,ₖ′₎), s≤j≤t ⇒ Aᵖᵢⱼ ≠ 0
supp(B₍ⱼ,ₚ′,ₖ′₎) ⊈ supp(B₍ᵢ,ₚ,ₖ₎) ⇒ Aᵖᵢⱼ = 0 is almost equivalent (if there are not duplicated knots).
supp(B₍ⱼ,ₚ′,ₖ′₎) ⊈ supp(B₍ᵢ,ₚ,ₖ₎) ⇒ Aᵖᵢⱼ = 0 is almost equivalent (if there are no duplicated knots).
- Rule-1: isdegenerate_R(P′, j) ⇒ Aᵖᵢⱼ = 0
Ideally this should be NaN, but we need to support `Rational` which does not include `NaN`.
- Rule-2: k′ⱼ = k′ⱼ₊ₚ′ ⇒ B₍ᵢ,ₚ,ₖ₎(t) = ∑ⱼAᵖᵢⱼB₍ⱼ,ₚ′,ₖ′₎(t) (t → k′ⱼ + 0)
Expand Down Expand Up @@ -205,8 +205,8 @@ function _changebasis_R(P::BSplineSpace{p,T,KnotVector{T}}, P′::BSplineSpace{p
=#

# Precalculate the range of j
j_begin, j_end = _find_j_begin_end_R(P, P′, i, j_begin, j_end)
j_range = j_begin:j_end
j_range = _find_j_range_R(P, P′, i, j_begin, j_end)
j_begin, j_end = extrema(j_range)

# Rule-0: outside of j_range
j_prev = j_begin-1
Expand Down Expand Up @@ -395,14 +395,14 @@ function _changebasis_I_old(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′})
return A
end

function _find_j_begin_end_I(P::BSplineSpace{p}, P′::BSplineSpace{p′}, i, j_begin, j_end):: Tuple{Int, Int} where {p, p′}
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′)
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
return j_begin:j_end
end

function _ΔAᵖ_I(Aᵖ⁻¹::AbstractMatrix, K::AbstractVector, K′::AbstractVector, i::Integer, j::Integer)
Expand Down Expand Up @@ -433,7 +433,7 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
=== j-rules for fixed i ===
- Rule-0: B₍ᵢ,ₚ,ₖ₎ ∈ span₍ᵤ₌ₛ…ₜ₎(B₍ᵤ,ₚ′,ₖ′₎), s≤j≤t ⇒ Aᵖᵢⱼ ≠ 0 (on B-spline space domain)
supp(B₍ⱼ,ₚ′,ₖ′₎) ⊈ supp(B₍ᵢ,ₚ,ₖ₎) ⇒ Aᵖᵢⱼ = 0 is almost equivalent (if there are not duplicated knots).
supp(B₍ⱼ,ₚ′,ₖ′₎) ⊈ supp(B₍ᵢ,ₚ,ₖ₎) ⇒ Aᵖᵢⱼ = 0 is almost equivalent (if there are no duplicated knots).
- Rule-1: isdegenerate_I(P′, j) ⇒ Aᵖᵢⱼ = 0
Ideally this should be NaN, but we need to support `Rational` which does not include `NaN`.
- Rule-2: k′ⱼ = k′ⱼ₊ₚ′ ⇒ B₍ᵢ,ₚ,ₖ₎(t) = ∑ⱼAᵖᵢⱼB₍ⱼ,ₚ′,ₖ′₎(t) (t → k′ⱼ + 0)
Expand Down Expand Up @@ -524,8 +524,8 @@ function _changebasis_I(P::BSplineSpace{p,T,<:AbstractKnotVector{T}}, P′::BSpl
=#

# Precalculate the range of j
j_begin, j_end = _find_j_begin_end_I(P, P′, i, j_begin, j_end)
j_range = j_begin:j_end
j_range = _find_j_range_I(P, P′, i, j_begin, j_end)
j_begin, j_end = extrema(j_range)

# Rule-0: outside of j_range
j_prev = j_begin-1
Expand Down
19 changes: 9 additions & 10 deletions test/test_ChangeBasis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,19 +17,18 @@
@test iszero(view(A, :, findall(BasicBSpline._iszeros_R(P′))))
end

function test_changebasis_I(P,P′)
@test P P′
A = @inferred changebasis(P,P′)
function test_changebasis_I(P1,P2)
@test P1 P2
A = @inferred changebasis_I(P1,P2)
@test A isa SparseMatrixCSC
@test !any(iszero.(A.nzval))
@test A BasicBSpline._changebasis_I(P,P′) == changebasis_I(P,P′)
n = dim(P)
n′ = dim(P′)
@test size(A) == (n,n′)
d = domain(P)
ts = range(extrema(d)..., length=20)
n1 = dim(P1)
n2 = dim(P2)
@test size(A) == (n1,n2)
d = domain(P1)
ts = range(extrema(d)..., length=21)[2:end-1]
for t in ts
@test norm(bsplinebasis.(P,1:n,t) - A*bsplinebasis.(P′,1:n′,t), Inf) < ε
@test norm(bsplinebasis.(P1,1:n1,t) - A*bsplinebasis.(P2,1:n2,t), Inf) < ε
end
end

Expand Down

0 comments on commit 7b5ba2b

Please sign in to comment.