Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
`changebasis` fails for non-float `BSplineSpace`
```julia
julia> P = BSplineSpace{3}(KnotVector{Int}(1:12))
BSplineSpace{3, Int64}(KnotVector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]))

julia> changebasis(P,P)
ERROR: InexactError: Int64(NaN)
Stacktrace:
  [1] Int64
    @ ./float.jl:812 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] fill!(A::SubArray{Int64, 2, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Vector{Int64}}, false}, x::Float64)
    @ Base ./multidimensional.jl:1062
  [4] copyto!
    @ ./broadcast.jl:921 [inlined]
  [5] materialize!
    @ ./broadcast.jl:871 [inlined]
  [6] materialize!
    @ ./broadcast.jl:868 [inlined]
  [7] _changebasis_R(P::BSplineSpace{0, Int64}, P′::BSplineSpace{0, Int64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:22
  [8] _changebasis_R(P::BSplineSpace{1, Int64}, P′::BSplineSpace{1, Int64}) (repeats 3 times)
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:26
  [9] changebasis(P::BSplineSpace{3, Int64}, P′::BSplineSpace{3, Int64})
    @ BasicBSpline ~/.julia/dev/BasicBSpline/src/_ChangeBasis.jl:162
 [10] top-level scope
    @ REPL[14]:1
```

Pull request description:
This PR fixes #143.

Patch:
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index dc6e3266e..a6ae57249 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -21,6 +21,7 @@ export 𝔫
 export AbstractBSplineManifold, AbstractBSplineSpace
 export fittingcontrolpoints
 
+include("_util.jl")
 include("_KnotVector.jl")
 include("_BSplineSpace.jl")
 include("_BSplineBasis.jl")
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index 7d6e30112..672dee1aa 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -1,8 +1,5 @@
 # B-Spline Basis Function
 
-@inline _d(a::T,b::T) where T = (U=StaticArrays.arithmetic_closure(T); ifelse(iszero(b), zero(U), U(a/b)))
-@inline _d(a,b) = _d(promote(a,b)...)
-
 @doc raw"""
 ``i``-th B-spline basis function.
 Right-sided limit version.
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index 89a48719a..8336a8660 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -108,7 +108,7 @@ const ⊑ = issqsubset
 
 ≃(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace) = (P1 ⊑ P2) & (P2 ⊑ P1)
 
-function isdegenerate(P::AbstractBSplineSpace{p}, i::Int) where p
+function isdegenerate(P::AbstractBSplineSpace{p}, i::Integer) where p
     k = knotvector(P)
     return k[i] == k[i+p+1]
 end
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index 30cfeb651..a06a93531 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -10,28 +10,34 @@ B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
 Assumption:
 * ``P ⊆ P′``
 """
-function _changebasis_R(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T})::Matrix{T} where {p,p′,T}
+_changebasis_R
+
+function _changebasis_R(P::AbstractBSplineSpace{0,T}, P′::AbstractBSplineSpace{p′,T′}) where {p′,T,T′}
+    P ⊆ P′ || throw(DomainError((P,P′),"P ⊆ P′ should be hold."))
+    U = StaticArrays.arithmetic_closure(promote_type(T,T′))
+    n = dim(P)
+    n′ = dim(P′)
+    A⁰ = U[bsplinesupport(P′, j) ⊆ bsplinesupport(P, i) for i in 1:n, j in 1:n′]
+    A⁰[:, findall(_iszeros(P′))] .= zero(U)  # Strictly this should be NaN, but we use zero here to support Rational.
+    return A⁰
+end
+
+function _changebasis_R(P::AbstractBSplineSpace{p,T}, P′::AbstractBSplineSpace{p′,T′}) where {p,p′,T,T′}
+    P ⊆ P′ || throw(DomainError((P,P′),"P ⊆ P′ should be hold."))
+    U = StaticArrays.arithmetic_closure(promote_type(T,T′))
     k = knotvector(P)
     k′ = knotvector(P′)
-    p₊ = p′ - p
-
-    if p == 0
-        n = length(k) - 1
-        n′ = length(k′) - p₊ - 1
-        A⁰ = T[bsplinesupport(BSplineSpace{p₊}(k′), j) ⊆ bsplinesupport(BSplineSpace{0}(k), i) for i in 1:n, j in 1:n′]
-        A⁰[:, findall(_iszeros(P′))] .= NaN
-        return A⁰
-    end
-
-    Aᵖ⁻¹ = _changebasis_R(_lower(P), _lower(P′)) # (n+1) × (n′+1) matrix
     n = dim(P)
     n′ = dim(P′)
+
+    Aᵖ⁻¹ = _changebasis_R(_lower(P), _lower(P′))  # (n+1) × (n′+1) matrix
+    Aᵖ = zeros(U, n, n′)  # n × n′ matrix
+
     Z = _iszeros(_lower(P′))
     W = findall(Z)
     K′ = [k′[i+p′] - k′[i] for i in 1:n′+1]
     K = [ifelse(k[i+p] ≠ k[i], 1 / (k[i+p] - k[i]), 0.0) for i in 1:n+1]
     Δ = (p / p′) * [K′[j] * (K[i] * Aᵖ⁻¹[i, j] - K[i+1] * Aᵖ⁻¹[i+1, j]) for i in 1:n, j in 1:n′+1]
-    Aᵖ = zeros(n, n′)
     Aᵖ[:, 1] = Δ[:, 1]
     Aᵖ[:, n′] = -Δ[:, n′+1]
 
@@ -48,7 +54,7 @@ function _changebasis_R(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T})::Matri
     for ȷ in 2:λ-1
         if Λ[ȷ] == 1
             # if B(i,p′,k′) = 0
-            Ãᵖ[ȷ] .= NaN
+            Ãᵖ[ȷ] .= zero(U)  # Strictly this should be NaN, but we use zero here to support Rational.
         end
     end
     for ȷ in 1:λ-1
@@ -90,21 +96,16 @@ B_{(i,p_1,k_1)} = \sum_{j}A_{i,j}B_{(j,p_2,k_2)}
 ```
 
 Assumption:
-* ``P1 ⊑ P2``
-* ``P2 ⊑ P1``
+* ``P1 ≃ P2``
 """
-function _changebasis_sim(P1::BSplineSpace{p1,T}, P2::BSplineSpace{p2,T}) where {p1,p2,T}
-    # if P1 ⋢ P2
-    #     error("P1 ⋢ P2")
-    # end
-    # if P2 ⋢ P1
-    #     error("P2 ⋢ P1")
-    # end
+function _changebasis_sim(P1::AbstractBSplineSpace{p1,T1}, P2::AbstractBSplineSpace{p2,T2}) where {p1,p2,T1,T2}
+    P1 ≃ P2 || throw(DomainError((P1,P2),"P1 ≃ P2 should be hold."))
+    U = StaticArrays.arithmetic_closure(promote_type(T1,T2))
     n = dim(P1)
     v = (knotvector(P1).vector)[p1+1:end-p1]
 
     if length(v) > p1
-        A = Matrix{T}(I, n, n)
+        A = Matrix{U}(I, n, n)
         # TODO: Fix below
         vvv = [v[1] * (p1-i+1) / (p1+1) + v[i+1] * (i) / (p1+1) for i in 1:p1]
         A1 = [bsplinebasis₊₀(P1,i,t) for i in 1:p1, t in vvv]
@@ -135,20 +136,13 @@ B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
 Assumption:
 * ``P ⊑ P′``
 """
-function _changebasis_I(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T})::Matrix{T} where {p, p′, T}
-    I = domain(P)
+function _changebasis_I(P::AbstractBSplineSpace{p,T}, P′::AbstractBSplineSpace{p′,T′}) where {p,p′,T,T′}
+    P ⊑ P′ || throw(DomainError((P,P′),"P ⊑ P′ should be hold."))
     k = knotvector(P)
     k′ = knotvector(P′)
-    p₊ = p′ - p
-
-    _P = BSplineSpace{p}(k[1+p:end-p] + p * KnotVector(k[1+p], k[end-p]))
-    # if dim(_P) ≠ dim(P)
-    #     error("dim(_P)≠dim(P)")
-    # end
-    _P′ = BSplineSpace{p′}(k′[1+p′:end-p′] + p′ * KnotVector(k′[1+p′], k′[end-p′]))
-    # if dim(_P′) ≠ dim(P′)
-    #     error("dim(_P′)≠dim(P′)")
-    # end
+
+    _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′)
@@ -157,18 +151,10 @@ function _changebasis_I(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T})::Matri
     return A
 end
 
-function changebasis(P::BSplineSpace, P′::BSplineSpace)
-    if P ⊆ P′
-        return _changebasis_R(P, P′)
-    elseif P ⊑ P′
-        return _changebasis_I(P, P′)
-    else
-        throw(DomainError((P, P′),"𝒫[p,k] ⊆ 𝒫[p′,k′] or 𝒫[p,k] ⊑ 𝒫[p′,k′] must hold."))
-    end
-end
-
 function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
-    changebasis(BSplineSpace(P),BSplineSpace(P′))
+    P ⊆ P′ && return _changebasis_R(P, P′)
+    P ⊑ P′ && return _changebasis_I(P, P′)
+    throw(DomainError((P, P′),"𝒫[p,k] ⊆ 𝒫[p′,k′] or 𝒫[p,k] ⊑ 𝒫[p′,k′] must hold."))
 end
 
 # TODO: Add changebasis(::BSplineDerivativeSpace, )
diff --git a/src/_util.jl b/src/_util.jl
new file mode 100644
index 000000000..0e325023e
--- /dev/null
+++ b/src/_util.jl
@@ -0,0 +1,7 @@
+# Division, devided by zero will be zero
+@inline function _d(a::T,b::T) where T
+    U=StaticArrays.arithmetic_closure(T)
+    iszero(b) && return zero(U)
+    return U(a/b)
+end
+@inline _d(a,b) = _d(promote(a,b)...)


Test patch:
diff --git a/test/test_ChangeBasis.jl b/test/test_ChangeBasis.jl
index a18146561..b71e0b2d2 100644
--- a/test/test_ChangeBasis.jl
+++ b/test/test_ChangeBasis.jl
@@ -2,47 +2,74 @@
     ε = 1e-14
     Random.seed!(42)
 
-    P1 = BSplineSpace{1}(KnotVector([1, 3, 5, 8]))
-    P2 = BSplineSpace{1}(KnotVector([1, 3, 5, 6, 8, 9]))
-    P3 = BSplineSpace{2}(KnotVector([1, 1, 3, 3, 5, 5, 8, 8]))
-    @test P1 ⊆ P2
-    @test P1 ⊆ P3
-    @test P2 ⊈ P3
-
-    n1, n2, n3 = dim(P1), dim(P2), dim(P3)
-    A12 = changebasis(P1, P2)
-    A13 = changebasis(P1, P3)
-    @test size(A12) == (n1, n2)
-    @test size(A13) == (n1, n3)
-    Δ12 = [bsplinebasis(P1, i, t) - sum(A12[i, j] * bsplinebasis(P2, j, t) for j in 1:n2) for i in 1:n1, t in 8 * rand(10)]
-    Δ13 = [bsplinebasis(P1, i, t) - sum(A13[i, j] * bsplinebasis(P3, j, t) for j in 1:n3) for i in 1:n1, t in 8 * rand(10)]
-    @test norm(Δ12) < ε
-    @test norm(Δ13) < ε
-
-    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
-
-    n4, n5 = dim(P4), dim(P5)
-    A45 = changebasis(P4, P5)
-    @test size(A45) == (n4, n5)
-    Δ45 = [bsplinebasis(P4, i, t) - sum(A45[i, j] * bsplinebasis(P5, j, t) for j in 1:n5) for i in 1:n1, t in 2 * rand(10) .+ 2]
-    @test norm(Δ45) < ε
-
-    P6 = BSplineSpace{p4-1}(knotvector(P4)[2:end-1])
-    P7 = BSplineSpace{p5-1}(knotvector(P5)[2:end-1])
-    @test P6 ⊑ P7
-    @test P6 ⊈ P7
-
-    n6, n7 = dim(P6), dim(P7)
-    A67 = changebasis(P6, P7)
-    @test size(A67) == (n6, n7)
-    Δ67 = [bsplinebasis(P6, i, t) - sum(A67[i, j] * bsplinebasis(P7, j, t) for j in 1:n7) for i in 1:n1, t in 2 * rand(10) .+ 2]
-    @test norm(Δ67) < ε
+    @testset "subseteq (R)" begin
+        P1 = BSplineSpace{1}(KnotVector([1, 3, 5, 8]))
+        P2 = BSplineSpace{1}(KnotVector([1, 3, 5, 6, 8, 9]))
+        P3 = BSplineSpace{2}(KnotVector([1, 1, 3, 3, 5, 5, 8, 8]))
+        @test P1 ⊆ P2
+        @test P1 ⊆ P3
+        @test P2 ⊈ P3
+
+        n1, n2, n3 = dim(P1), dim(P2), dim(P3)
+        A12 = changebasis(P1, P2)
+        A13 = changebasis(P1, P3)
+        @test size(A12) == (n1, n2)
+        @test size(A13) == (n1, n3)
+        Δ12 = [bsplinebasis(P1, i, t) - sum(A12[i, j] * bsplinebasis(P2, j, t) for j in 1:n2) for i in 1:n1, t in 8 * rand(10)]
+        Δ13 = [bsplinebasis(P1, i, t) - sum(A13[i, j] * bsplinebasis(P3, j, t) for j in 1:n3) for i in 1:n1, t in 8 * rand(10)]
+        @test norm(Δ12) < ε
+        @test norm(Δ13) < ε
+    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
+
+        n4, n5 = dim(P4), dim(P5)
+        A45 = changebasis(P4, P5)
+        @test size(A45) == (n4, n5)
+        Δ45 = [bsplinebasis(P4, i, t) - sum(A45[i, j] * bsplinebasis(P5, j, t) for j in 1:n5) for i in 1:n4, t in 2 * rand(10) .+ 2]
+        @test norm(Δ45) < ε
+
+        P6 = BSplineSpace{p4-1}(knotvector(P4)[2:end-1])
+        P7 = BSplineSpace{p5-1}(knotvector(P5)[2:end-1])
+        @test P6 ⊑ P7
+        @test P6 ⊈ P7
+
+        n6, n7 = dim(P6), dim(P7)
+        A67 = changebasis(P6, P7)
+        @test size(A67) == (n6, n7)
+        Δ67 = [bsplinebasis(P6, i, t) - sum(A67[i, j] * bsplinebasis(P7, j, t) for j in 1:n7) for i in 1:n6, t in 2 * rand(10) .+ 2]
+        @test norm(Δ67) < ε
+    end
+
+    @testset "non-float" begin
+        k = UniformKnotVector(1:15)
+        P = UniformBSplineSpace{3}(k)
+        @test changebasis(P,P) isa Matrix{Float64}
+        @test BasicBSpline._changebasis_R(P,P) isa Matrix{Float64}
+        @test BasicBSpline._changebasis_I(P,P) isa Matrix{Float64}
+        @test BasicBSpline._changebasis_sim(P,P) isa Matrix{Float64}
+
+        k = UniformKnotVector(1:15//1)
+        P = UniformBSplineSpace{3}(k)
+        @test changebasis(P,P) isa Matrix{Rational{Int}}
+        @test BasicBSpline._changebasis_R(P,P) isa Matrix{Rational{Int}}
+        @test BasicBSpline._changebasis_I(P,P) isa Matrix{Rational{Int}}
+        @test BasicBSpline._changebasis_sim(P,P) isa Matrix{Rational{Int}}
+
+        k = UniformKnotVector(1:BigInt(15))
+        P = UniformBSplineSpace{3}(k)
+        @test changebasis(P,P) isa Matrix{BigFloat}
+        @test BasicBSpline._changebasis_R(P,P) isa Matrix{BigFloat}
+        @test BasicBSpline._changebasis_I(P,P) isa Matrix{BigFloat}
+        @test BasicBSpline._changebasis_sim(P,P) isa Matrix{BigFloat}
+    end
 end