Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
`refinement` is not type-stable
The following function `r` is not type-stable.

```julia
julia> using BasicBSpline

julia> k = KnotVector(rand(12))
KnotVector([0.024885878108435344, 0.1135410183212171, 0.17203785287918527, 0.2931216734788268, 0.46010004333534293, 0.6340248134707075, 0.6826197801697019, 0.736317312313712, 0.7375293828560209, 0.9130250689726209, 0.9602880185745573, 0.9647497776183738])

julia> P = BSplineSpace{2}(k)
BSplineSpace{2, Float64}(KnotVector([0.024885878108435344, 0.1135410183212171, 0.17203785287918527, 0.2931216734788268, 0.46010004333534293, 0.6340248134707075, 0.6826197801697019, 0.736317312313712, 0.7375293828560209, 0.9130250689726209, 0.9602880185745573, 0.9647497776183738]))

julia> a = rand(dim(P))
9-element Vector{Float64}:
 0.7352550606757859
 0.4586881993320009
 0.00534648579157071
 0.27377101932655756
 0.20868358047740354
 0.7353046147958138
 0.24689508125636628
 0.14759107033296426
 0.2632878277228856

julia> M = BSplineManifold(a,(P,))
BSplineManifold{1, (2,), Float64, Tuple{BSplineSpace{2, Float64}}}((BSplineSpace{2, Float64}(KnotVector([0.024885878108435344, 0.1135410183212171, 0.17203785287918527, 0.2931216734788268, 0.46010004333534293, 0.6340248134707075, 0.6826197801697019, 0.736317312313712, 0.7375293828560209, 0.9130250689726209, 0.9602880185745573, 0.9647497776183738])),), [0.7352550606757859, 0.4586881993320009, 0.00534648579157071, 0.27377101932655756, 0.20868358047740354, 0.7353046147958138, 0.24689508125636628, 0.14759107033296426, 0.2632878277228856])

julia> r(M::BSplineManifold{1}) = refinement(M, p₊=(1,))
r (generic function with 1 method)

julia> r(M)
BSplineManifold{1, (3,), Float64, Tuple{BSplineSpace{3, Float64}}}((BSplineSpace{3, Float64}(KnotVector([0.024885878108435344, 0.1135410183212171, 0.17203785287918527, 0.17203785287918527, 0.2931216734788268, 0.2931216734788268, 0.46010004333534293, 0.46010004333534293, 0.6340248134707075, 0.6340248134707075, 0.6826197801697019, 0.6826197801697019, 0.736317312313712, 0.736317312313712, 0.7375293828560209, 0.7375293828560209, 0.9130250689726209, 0.9130250689726209, 0.9602880185745573, 0.9647497776183738])),), [0.7052253117995405, 0.5208474042370173, 0.3951691068914191, 0.09294129786446556, 0.04917231761018773, 0.22812200663351226, 0.2568132458999762, 0.21342162000054016, 0.29207559167126984, 0.6431562812168766, 0.5760951423561846, 0.25048878666321966, 0.24666803318436095, 0.18046535923542628, 0.1779741525332959, 0.25510532412657677])

julia> @code_warntype r(M)
MethodInstance for r(::BSplineManifold{1, (2,), Float64, Tuple{BSplineSpace{2, Float64}}})
  from r(M::BSplineManifold{1}) in Main at REPL[6]:1
Arguments
  #self#::Core.Const(r)
  M::BSplineManifold{1, (2,), Float64, Tuple{BSplineSpace{2, Float64}}}
Body::BSplineManifold
1 ─ %1 = (:p₊,)::Core.Const((:p₊,))
│   %2 = Core.apply_type(Core.NamedTuple, %1)::Core.Const(NamedTuple{(:p₊,)})
│   %3 = Core.tuple(1)::Core.Const((1,))
│   %4 = Core.tuple(%3)::Core.Const(((1,),))
│   %5 = (%2)(%4)::Core.Const((p₊ = (1,),))
│   %6 = Core.kwfunc(Main.refinement)::Core.Const(BasicBSpline.var"#refinement##kw"())
│   %7 = (%6)(%5, Main.refinement, M)::BSplineManifold
└──      return %7
```

This can be solved by adding support `Val` just like in StaticArrays.jl.

```julia
julia> using StaticArrays

julia> A = @SMatrix rand(3,4)
3×4 SMatrix{3, 4, Float64, 12} with indices SOneTo(3)×SOneTo(4):
 0.968585  0.735832  0.0643772  0.974502
 0.392632  0.863711  0.687421   0.0488318
 0.292027  0.979669  0.404775   0.538932

julia> sum(A, dims=1)
1×4 SMatrix{1, 4, Float64, 4} with indices SOneTo(1)×SOneTo(4):
 1.65324  2.57921  1.15657  1.56227

julia> sum(A, dims=Val(1))
1×4 SMatrix{1, 4, Float64, 4} with indices SOneTo(1)×SOneTo(4):
 1.65324  2.57921  1.15657  1.56227

julia> @code_warntype sum(A, dims=1)
MethodInstance for (::Base.var"#sum##kw")(::NamedTuple{(:dims,), Tuple{Int64}}, ::typeof(sum), ::SMatrix{3, 4, Float64, 12})
  from (::Base.var"#sum##kw")(::Any, ::typeof(sum), a::StaticArray{<:Tuple, T}) where T in StaticArrays at /home/hyrodium/.julia/dev/StaticArrays/src/mapreduce.jl:250
Static Parameters
  T = Float64
Arguments
  _::Core.Const(Base.var"#sum##kw"())
  @_2::NamedTuple{(:dims,), Tuple{Int64}}
  @_3::Core.Const(sum)
  a::SMatrix{3, 4, Float64, 12}
Locals
  dims::Int64
  @_6::Int64
Body::Any
1 ─ %1  = Base.haskey(@_2, :dims)::Core.Const(true)
│         Core.typeassert(%1, Core.Bool)
│         (@_6 = Base.getindex(@_2, :dims))
└──       goto #3
2 ─       Core.Const(:(@_6 = StaticArrays.:(:)))
3 ┄ %6  = @_6::Int64
│         (dims = %6)
│   %8  = (:dims,)::Core.Const((:dims,))
│   %9  = Core.apply_type(Core.NamedTuple, %8)::Core.Const(NamedTuple{(:dims,)})
│   %10 = Base.structdiff(@_2, %9)::Core.Const(NamedTuple())
│   %11 = Base.pairs(%10)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())
│   %12 = Base.isempty(%11)::Core.Const(true)
│         Core.typeassert(%12, Core.Bool)
└──       goto #5
4 ─       Core.Const(:(Base.kwerr(@_2, @_3, a)))
5 ┄ %16 = StaticArrays.:(var"#sum#184")(dims, @_3, a)::Any
└──       return %16


julia> @code_warntype sum(A, dims=Val(1))
MethodInstance for (::Base.var"#sum##kw")(::NamedTuple{(:dims,), Tuple{Val{1}}}, ::typeof(sum), ::SMatrix{3, 4, Float64, 12})
  from (::Base.var"#sum##kw")(::Any, ::typeof(sum), a::StaticArray{<:Tuple, T}) where T in StaticArrays at /home/hyrodium/.julia/dev/StaticArrays/src/mapreduce.jl:250
Static Parameters
  T = Float64
Arguments
  _::Core.Const(Base.var"#sum##kw"())
  @_2::Core.Const((dims = Val{1}(),))
  @_3::Core.Const(sum)
  a::SMatrix{3, 4, Float64, 12}
Locals
  dims::Val{1}
  @_6::Val{1}
Body::SMatrix{1, 4, Float64, 4}
1 ─ %1  = Base.haskey(@_2, :dims)::Core.Const(true)
│         Core.typeassert(%1, Core.Bool)
│         (@_6 = Base.getindex(@_2, :dims))
└──       goto #3
2 ─       Core.Const(:(@_6 = StaticArrays.:(:)))
3 ┄ %6  = @_6::Core.Const(Val{1}())
│         (dims = %6)
│   %8  = (:dims,)::Core.Const((:dims,))
│   %9  = Core.apply_type(Core.NamedTuple, %8)::Core.Const(NamedTuple{(:dims,)})
│   %10 = Base.structdiff(@_2, %9)::Core.Const(NamedTuple())
│   %11 = Base.pairs(%10)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())
│   %12 = Base.isempty(%11)::Core.Const(true)
│         Core.typeassert(%12, Core.Bool)
└──       goto #5
4 ─       Core.Const(:(Base.kwerr(@_2, @_3, a)))
5 ┄ %16 = StaticArrays.:(var"#sum#184")(dims, @_3, a)::SMatrix{1, 4, Float64, 4}
└──       return %16
```

Pull request description:
This PR fixes #244.

Patch:
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index bc2a37d1d..96cb00e4c 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -329,6 +329,14 @@ function intervalindex(P::AbstractBSplineSpace{p},t::Real) where p
     return searchsortedlast(v,t)+1
 end
 
+# TODO remove these methods in the next breaking release
+for fname in (:expandspace, :expandspace_I, :expandspace_R)
+    @eval function $fname(P::BSplineSpace{p,T}; p₊::Integer=0, k₊=KnotVector{T}()) where {p,T}
+        P′ = $fname(P, Val(p₊), k₊)
+        return P′
+    end
+end
+
 """
 Expand B-spline space with given additional degree and knotvector.
 This function is compatible with `issqsubset` (`⊑`)
@@ -355,7 +363,9 @@ julia> domain(P′)
 2.5..9.0
 ```
 """
-function expandspace_I(P::BSplineSpace{p,T}; p₊::Integer=0, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,T}
+function expandspace_I end
+
+function expandspace_I(P::BSplineSpace{p,T}, ::Val{p₊}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,p₊,T}
     k = knotvector(P)
     k̂ = unique(k[1+p:end-p])
     p′ = p + p₊
@@ -364,6 +374,13 @@ function expandspace_I(P::BSplineSpace{p,T}; p₊::Integer=0, k₊::AbstractKnot
     return P′
 end
 
+function expandspace_I(P::BSplineSpace{p,T}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,T}
+    k = knotvector(P)
+    k′ = k + k₊
+    P′ = BSplineSpace{p}(k′)
+    return P′
+end
+
 """
 Expand B-spline space with given additional degree and knotvector.
 This function is compatible with `issubset` (`⊆`).
@@ -390,7 +407,9 @@ julia> domain(P′)
 1.5..9.5
 ```
 """
-function expandspace_R(P::BSplineSpace{p,T}; p₊::Integer=0, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,T}
+function expandspace_R end
+
+function expandspace_R(P::BSplineSpace{p,T}, ::Val{p₊}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,p₊,T}
     k = knotvector(P)
     p′ = p + p₊
     k′ = k + p₊*k + k₊
@@ -398,10 +417,21 @@ function expandspace_R(P::BSplineSpace{p,T}; p₊::Integer=0, k₊::AbstractKnot
     return P′
 end
 
+function expandspace_R(P::BSplineSpace{p,T}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,T}
+    k = knotvector(P)
+    k′ = k + k₊
+    P′ = BSplineSpace{p}(k′)
+    return P′
+end
+
 """
 Expand B-spline space with given additional degree and knotvector.
 The behavior of `expandspace` is same as `expandspace_I`.
 """
-function expandspace(P::BSplineSpace{p,T}; p₊=0, k₊=KnotVector{T}()) where {p,T}
-    expandspace_I(P,p₊=p₊,k₊=k₊)
+function expandspace(P::BSplineSpace{p,T}, _p₊::Val{p₊}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,p₊,T}
+    expandspace_I(P,_p₊,k₊)
+end
+
+function expandspace(P::BSplineSpace{p,T}, k₊::AbstractKnotVector=EmptyKnotVector{T}()) where {p,T}
+    expandspace_I(P,k₊)
 end
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index d248af3b0..ec9aeac86 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -114,13 +114,10 @@ function _changebasis_sim(P1::AbstractBSplineSpace{p,T1}, P2::AbstractBSplineSpa
 end
 
 @generated function _derivatives_at_left(P::AbstractBSplineSpace{p,T}) where {p, T}
-    hcat_args = Any[:hcat]
-    for r in 0:p-1
-        push!(hcat_args, :(pop(bsplinebasisall(BSplineDerivativeSpace{$(r)}(P),1,a))))
-    end
+    args = [:(pop(bsplinebasisall(BSplineDerivativeSpace{$(r)}(P),1,a))) for r in 0:p-1]
     quote
         a, _ = extrema(domain(P))
-        $(Expr(:call, hcat_args...))
+        $(Expr(:call, :hcat, args...))
     end
 end
 function _derivatives_at_left(::AbstractBSplineSpace{0,T}) where {T}
@@ -129,14 +126,11 @@ function _derivatives_at_left(::AbstractBSplineSpace{0,T}) where {T}
 end
 
 @generated function _derivatives_at_right(P::AbstractBSplineSpace{p,T}) where {p, T}
-    hcat_args = Any[:hcat]
-    for r in 0:p-1
-        push!(hcat_args, :(popfirst(bsplinebasisall(BSplineDerivativeSpace{$(r)}(P),n-p,b))))
-    end
+    args = [:(popfirst(bsplinebasisall(BSplineDerivativeSpace{$(r)}(P),n-p,b))) for r in 0:p-1]
     quote
         n = dim(P)
         _, b = extrema(domain(P))
-        $(Expr(:call, hcat_args...))
+        $(Expr(:call, :hcat, args...))
     end
 end
 function _derivatives_at_right(::AbstractBSplineSpace{0,T}) where {T}
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index 55364e028..efe7a5ee5 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -71,7 +71,7 @@ function innerproduct_R(P::AbstractBSplineSpace{p}) where p
     return Symmetric(A)
 end
 
-function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p₁}
+function innerproduct_I(func, Ps::Tuple{AbstractBSplineSpace{p₁}}) where {p₁}
     P₁, = Ps
     n₁ = dim(P₁)
     k₁ = knotvector(P₁)
@@ -98,7 +98,7 @@ function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p
     return b
 end
 
-function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂}}) where {p₁,p₂}
+function innerproduct_I(func, Ps::Tuple{AbstractBSplineSpace{p₁},AbstractBSplineSpace{p₂}}) where {p₁,p₂}
     P₁,P₂ = Ps
     n₁,n₂ = dim.(Ps)
     k₁,k₂ = knotvector.(Ps)
@@ -137,7 +137,7 @@ function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractB
     return b
 end
 
-function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂},<:AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
+function innerproduct_I(func, Ps::Tuple{AbstractBSplineSpace{p₁},AbstractBSplineSpace{p₂},AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
     P₁,P₂,P₃ = Ps
     n₁,n₂,n₃ = dim.(Ps)
     k₁,k₂,k₃ = knotvector.(Ps)
@@ -191,7 +191,7 @@ function innerproduct_I(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractB
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p₁}
+function innerproduct_R(func, Ps::Tuple{AbstractBSplineSpace{p₁}}) where {p₁}
     P₁, = Ps
     n₁ = dim(P₁)
     k₁ = knotvector(P₁)
@@ -214,7 +214,7 @@ function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁}}) where {p
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂}}) where {p₁,p₂}
+function innerproduct_R(func, Ps::Tuple{AbstractBSplineSpace{p₁},AbstractBSplineSpace{p₂}}) where {p₁,p₂}
     P₁,P₂ = Ps
     n₁,n₂ = dim.(Ps)
     k₁,k₂ = knotvector.(Ps)
@@ -246,7 +246,7 @@ function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractB
     return b
 end
 
-function innerproduct_R(func, Ps::Tuple{<:AbstractBSplineSpace{p₁},<:AbstractBSplineSpace{p₂},<:AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
+function innerproduct_R(func, Ps::Tuple{AbstractBSplineSpace{p₁},AbstractBSplineSpace{p₂},AbstractBSplineSpace{p₃}}) where {p₁,p₂,p₃}
     P₁,P₂,P₃ = Ps
     n₁,n₂,n₃ = dim.(Ps)
     k₁,k₂,k₃ = knotvector.(Ps)
@@ -322,14 +322,14 @@ function innerproduct_R(P::UniformBSplineSpace{p,T}) where {p,T}
 end
 
 for (fname_fit, fname_inner) in ((:fittingcontrolpoints_I, :innerproduct_I), (:fittingcontrolpoints_R, :innerproduct_R))
-    @eval function $fname_fit(func, P::Tuple{<:AbstractBSplineSpace{p1}}) where {p1}
+    @eval function $fname_fit(func, P::Tuple{AbstractBSplineSpace{p1}}) where {p1}
         P1, = P
         b = $fname_inner(func, P)
         A = $fname_inner(P1)
         return _leftdivision(A, b)
     end
 
-    @eval function $fname_fit(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}) where {p1,p2}
+    @eval function $fname_fit(func, P::Tuple{AbstractBSplineSpace{p1}, AbstractBSplineSpace{p2}}) where {p1,p2}
         P1, P2 = P
         n1, n2 = dim.(P)
         A1, A2 = $fname_inner(P1), $fname_inner(P2)
@@ -340,7 +340,7 @@ for (fname_fit, fname_inner) in ((:fittingcontrolpoints_I, :innerproduct_I), (:f
         return reshape(_leftdivision(_A, _b), n1, n2)
     end
 
-    @eval function $fname_fit(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}) where {p1,p2,p3}
+    @eval function $fname_fit(func, P::Tuple{AbstractBSplineSpace{p1}, AbstractBSplineSpace{p2}, AbstractBSplineSpace{p3}}) where {p1,p2,p3}
         P1, P2, P3 = P
         n1, n2, n3 = dim(P1), dim(P2), dim(P3)
         A1, A2, A3 = $fname_inner(P1), $fname_inner(P2), $fname_inner(P3)
diff --git a/src/_Refinement.jl b/src/_Refinement.jl
index 6c6349cec..8a4df8d57 100644
--- a/src/_Refinement.jl
+++ b/src/_Refinement.jl
@@ -8,68 +8,92 @@ Refinement of B-spline manifold with given B-spline spaces.
 """
 refinement
 
-function refinement(M::BSplineManifold{1}, Ps′::Tuple{<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::BSplineManifold{1}, Ps′::NTuple{1, AbstractBSplineSpace})
+    P1, = bsplinespaces(M)
+    P1′, = Ps′
+    n1 = dim(P1)
+    n1′ = dim(P1′)
+    A1 = changebasis(P1, P1′)
     a = controlpoints(M)
-    (n1,) = dim.(Ps)
-    (n1′,) = dim.(Ps′)
-    (A1,) = changebasis.(Ps, Ps′)
 
     a′ = [sum(A1[I₁,J₁] * a[I₁] for I₁ in 1:n1) for J₁ in 1:n1′]
     return BSplineManifold(a′, Ps′)
 end
-function refinement(M::BSplineManifold{2}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::BSplineManifold{2}, Ps′::NTuple{2, AbstractBSplineSpace})
+    P1, P2 = bsplinespaces(M)
+    P1′, P2′ = Ps′
+    n1 = dim(P1)
+    n2 = dim(P2)
+    n1′ = dim(P1′)
+    n2′ = dim(P2′)
+    A1 = changebasis(P1, P1′)
+    A2 = changebasis(P2, P2′)
     a = controlpoints(M)
-    (n1,n2) = dim.(Ps)
-    (n1′,n2′) = dim.(Ps′)
-    (A1,A2) = changebasis.(Ps, Ps′)
 
     a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * a[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′]
     return BSplineManifold(a′, Ps′)
 end
-function refinement(M::BSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace,<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::BSplineManifold{3}, Ps′::NTuple{3, AbstractBSplineSpace})
+    P1, P2, P3 = bsplinespaces(M)
+    P1′, P2′, P3′ = Ps′
+    n1 = dim(P1)
+    n2 = dim(P2)
+    n3 = dim(P3)
+    n1′ = dim(P1′)
+    n2′ = dim(P2′)
+    n3′ = dim(P3′)
+    A1 = changebasis(P1, P1′)
+    A2 = changebasis(P2, P2′)
+    A3 = changebasis(P3, P3′)
     a = controlpoints(M)
-    (n1,n2,n3) = dim.(Ps)
-    (n1′,n2′,n3′) = dim.(Ps′)
-    (A1,A2,A3) = changebasis.(Ps, Ps′)
 
     a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * A3[I₃,J₃] * a[I₁,I₂,I₃] for I₁ in 1:n1, I₂ in 1:n2, I₃ in 1:n3) for J₁ in 1:n1′, J₂ in 1:n2′, J₃ in 1:n3′]
     return BSplineManifold(a′, Ps′)
 end
 
-function refinement(M::RationalBSplineManifold{1}, Ps′::Tuple{<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::RationalBSplineManifold{1}, Ps′::NTuple{1, AbstractBSplineSpace})
+    P1, = bsplinespaces(M)
+    P1′, = Ps′
+    n1 = dim(P1)
+    n1′ = dim(P1′)
+    A1 = changebasis(P1, P1′)
     a = controlpoints(M)
     w = weights(M)
-    (n1,) = dim.(Ps)
-    (n1′,) = dim.(Ps′)
-    (A1,) = changebasis.(Ps, Ps′)
 
     w′ = [sum(A1[I₁,J₁] * w[I₁] for I₁ in 1:n1) for J₁ in 1:n1′]
     a′ = [sum(A1[I₁,J₁] * a[I₁] * w[I₁] for I₁ in 1:n1) for J₁ in 1:n1′] ./ w′
     return RationalBSplineManifold(a′, w′, Ps′)
 end
-function refinement(M::RationalBSplineManifold{2}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::RationalBSplineManifold{2}, Ps′::NTuple{2, AbstractBSplineSpace})
+    P1, P2 = bsplinespaces(M)
+    P1′, P2′ = Ps′
+    n1 = dim(P1)
+    n2 = dim(P2)
+    n1′ = dim(P1′)
+    n2′ = dim(P2′)
+    A1 = changebasis(P1, P1′)
+    A2 = changebasis(P2, P2′)
     a = controlpoints(M)
     w = weights(M)
-    (n1,n2) = dim.(Ps)
-    (n1′,n2′) = dim.(Ps′)
-    (A1,A2) = changebasis.(Ps, Ps′)
 
     w′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * w[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′]
     a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * a[I₁,I₂] * w[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′] ./ w′
     return RationalBSplineManifold(a′, w′, Ps′)
 end
-function refinement(M::RationalBSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace,<:AbstractBSplineSpace})
-    Ps = bsplinespaces(M)
+function refinement(M::RationalBSplineManifold{3}, Ps′::NTuple{3, AbstractBSplineSpace})
+    P1, P2, P3 = bsplinespaces(M)
+    P1′, P2′, P3′ = Ps′
+    n1 = dim(P1)
+    n2 = dim(P2)
+    n3 = dim(P3)
+    n1′ = dim(P1′)
+    n2′ = dim(P2′)
+    n3′ = dim(P3′)
+    A1 = changebasis(P1, P1′)
+    A2 = changebasis(P2, P2′)
+    A3 = changebasis(P3, P3′)
     a = controlpoints(M)
     w = weights(M)
-    (n1,n2,n3) = dim.(Ps)
-    (n1′,n2′,n3′) = dim.(Ps′)
-    (A1,A2,A3) = changebasis.(Ps, Ps′)
 
     w′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * A3[I₃,J₃] * w[I₁,I₂,I₃] for I₁ in 1:n1, I₂ in 1:n2, I₃ in 1:n3) for J₁ in 1:n1′, J₂ in 1:n2′, J₃ in 1:n3′]
     a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * A3[I₃,J₃] * a[I₁,I₂,I₃] * w[I₁,I₂,I₃] for I₁ in 1:n1, I₂ in 1:n2, I₃ in 1:n3) for J₁ in 1:n1′, J₂ in 1:n2′, J₃ in 1:n3′] ./ w′
@@ -93,3 +117,46 @@ function refinement(M::AbstractManifold{Dim}; p₊::Union{Nothing,NTuple{Dim,Int
 
     return refinement(M, tuple(Ps′...))
 end
+
+@generated function refinement(M::AbstractManifold{Dim},
+                               p₊::NTuple{Dim, Val},
+                               k₊::NTuple{Dim, AbstractKnotVector}=ntuple(i->EmptyKnotVector(), Val(Dim))) where Dim
+    Ps = [Symbol(:P,i) for i in 1:Dim]
+    Ps′ = [Symbol(:P,i,"′") for i in 1:Dim]
+    ks = [Symbol(:k,i,:₊) for i in 1:Dim]
+    ps = [Symbol(:p,i,:₊) for i in 1:Dim]
+    exP = Expr(:tuple, Ps...)
+    exP′ = Expr(:tuple, Ps′...)
+    exk = Expr(:tuple, ks...)
+    exp = Expr(:tuple, ps...)
+    exs = [:($(Symbol(:P,i,"′")) = expandspace($(Symbol(:P,i)), $(Symbol(:p,i,:₊)), $(Symbol(:k,i,:₊)))) for i in 1:Dim]
+    Expr(
+        :block,
+        :($exP = bsplinespaces(M)),
+        :($exp = p₊),
+        :($exk = k₊),
+        exs...,
+        :(return refinement(M, $(exP′)))
+    )
+end
+
+@generated function refinement(M::AbstractManifold{Dim},
+                               k₊::NTuple{Dim, AbstractKnotVector}=ntuple(i->EmptyKnotVector(), Val(Dim))) where Dim
+    Ps = [Symbol(:P,i) for i in 1:Dim]
+    Ps′ = [Symbol(:P,i,"′") for i in 1:Dim]
+    ks = [Symbol(:k,i,:₊) for i in 1:Dim]
+    exP = Expr(:tuple, Ps...)
+    exP′ = Expr(:tuple, Ps′...)
+    exk = Expr(:tuple, ks...)
+    exs = [:($(Symbol(:P,i,"′")) = expandspace($(Symbol(:P,i)), $(Symbol(:k,i,:₊)))) for i in 1:Dim]
+    Expr(
+        :block,
+        :($exP = bsplinespaces(M)),
+        :($exk = k₊),
+        exs...,
+        :(return refinement(M, $(exP′)))
+    )
+end
+
+# resolve ambiguities
+refinement(M::BasicBSpline.AbstractManifold{0}, ::Tuple{}) = M


Test patch:
diff --git a/test/test_BSplineSpace.jl b/test/test_BSplineSpace.jl
index c58d748f2..1d0aceab4 100644
--- a/test/test_BSplineSpace.jl
+++ b/test/test_BSplineSpace.jl
@@ -124,10 +124,10 @@
         @test P3 == expandspace_R(P3) == expandspace_I(P3) == expandspace(P3)
         @test P4 == expandspace_R(P4) == expandspace_I(P4) == expandspace(P4)
 
-        @test P1 ⊆ expandspace_R(P1, p₊=1)
-        @test P2 ⊆ expandspace_R(P2, p₊=1)
-        @test P3 ⊆ expandspace_R(P3, p₊=1)
-        @test P4 ⊆ expandspace_R(P4, p₊=1)
+        @test P1 ⊆ @inferred expandspace_R(P1, Val(1))
+        @test P2 ⊆ @inferred expandspace_R(P2, Val(1))
+        @test P3 ⊆ @inferred expandspace_R(P3, Val(1))
+        @test P4 ⊆ @inferred expandspace_R(P4, Val(1))
 
         @test P1 ⊆ expandspace_R(P1, k₊=KnotVector([1,2,3]))
         @test P2 ⊆ expandspace_R(P2, k₊=KnotVector([1,2,3]))
@@ -140,10 +140,10 @@
         @test knotvector(P4) + KnotVector([1,2,3]) ⊆ knotvector(expandspace_R(P4, k₊=KnotVector([1,2,3])))
 
         # P1 and P4 is not subset of expandspace_I.
-        @test P1 ⋢ expandspace_I(P1, p₊=1) == expandspace(P1, p₊=1)
-        @test P2 ⊑ expandspace_I(P2, p₊=1) == expandspace(P2, p₊=1)
-        @test P3 ⊑ expandspace_I(P3, p₊=1) == expandspace(P3, p₊=1)
-        @test P4 ⋢ expandspace_I(P4, p₊=1) == expandspace(P4, p₊=1)
+        @test P1 ⋢ expandspace_I(P1, Val(1)) == expandspace(P1, Val(1)) == expandspace_I(P1, p₊=1) == expandspace(P1, p₊=1)
+        @test P2 ⊑ expandspace_I(P2, Val(1)) == expandspace(P2, Val(1)) == expandspace_I(P2, p₊=1) == expandspace(P2, p₊=1)
+        @test P3 ⊑ expandspace_I(P3, Val(1)) == expandspace(P3, Val(1)) == expandspace_I(P3, p₊=1) == expandspace(P3, p₊=1)
+        @test P4 ⋢ expandspace_I(P4, Val(1)) == expandspace(P4, Val(1)) == expandspace_I(P4, p₊=1) == expandspace(P4, p₊=1)
 
         # That was because P1 and P4 is not nondegenerate.
         @test isdegenerate_I(P1)
diff --git a/test/test_Fitting.jl b/test/test_Fitting.jl
index ad1379f86..f3e01eecb 100644
--- a/test/test_Fitting.jl
+++ b/test/test_Fitting.jl
@@ -32,7 +32,7 @@
         a_org = [Point(i1, rand()) for i1 in 1:n1]
         M = BSplineManifold(a_org, (P1,))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′,))
         a_ref = controlpoints(M′)
@@ -54,7 +54,7 @@
         a_org = [Point(i1, rand()) for i1 in 1:n1]
         M = BSplineManifold(a_org, (P1,))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′,))
         a_ref = controlpoints(M′)
@@ -78,8 +78,8 @@
         a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
         M = BSplineManifold(a_org, (P1, P2))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
-        P2′ = expandspace(P2, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
+        P2′ = expandspace(P2, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′, P2′))
         a_ref = controlpoints(M′)
@@ -105,8 +105,8 @@
         a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
         M = BSplineManifold(a_org, (P1, P2))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
-        P2′ = expandspace(P2, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
+        P2′ = expandspace(P2, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′, P2′))
         a_ref = controlpoints(M′)
@@ -133,9 +133,9 @@
         a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
         M = BSplineManifold(a_org, (P1, P2, P3))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
-        P2′ = expandspace(P2, p₊=1, k₊=KnotVector(rand(2)))
-        P3′ = expandspace(P3, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
+        P2′ = expandspace(P2, Val(1), KnotVector(rand(2)))
+        P3′ = expandspace(P3, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′, P2′, P3′))
         a_ref = controlpoints(M′)
@@ -165,9 +165,9 @@
         a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
         M = BSplineManifold(a_org, (P1, P2, P3))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
-        P2′ = expandspace(P2, p₊=1, k₊=KnotVector(rand(2)))
-        P3′ = expandspace(P3, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
+        P2′ = expandspace(P2, Val(1), KnotVector(rand(2)))
+        P3′ = expandspace(P3, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′, P2′, P3′))
         a_ref = controlpoints(M′)
@@ -190,8 +190,8 @@
         a_org = [rand(ComplexF64) for i1 in 1:n1, i2 in 1:n2]
         M = BSplineManifold(a_org, (P1, P2))
 
-        P1′ = expandspace(P1, p₊=1, k₊=KnotVector(rand(2)))
-        P2′ = expandspace(P2, p₊=1, k₊=KnotVector(rand(2)))
+        P1′ = expandspace(P1, Val(1), KnotVector(rand(2)))
+        P2′ = expandspace(P2, Val(1), KnotVector(rand(2)))
 
         M′ = refinement(M, (P1′, P2′))
         a_ref = controlpoints(M′)
diff --git a/test/test_Refinement.jl b/test/test_Refinement.jl
index 8fc022b59..604db0c23 100644
--- a/test/test_Refinement.jl
+++ b/test/test_Refinement.jl
@@ -10,7 +10,7 @@
     P3 = BSplineSpace{p3}(k3)
 
     P1′ = BSplineSpace{3}(k1+KnotVector(1.2,5.5,7.2,8.5))
-    P2′ = expandspace(P2, p₊=1, k₊=KnotVector(1.8))
+    P2′ = expandspace(P2, Val(1), KnotVector(1.8))
     P3′ = BSplineSpace{4}(UniformKnotVector(-3:6))
 
     n1,n2,n3 = dim(P1),dim(P2),dim(P3)
@@ -26,24 +26,24 @@
 
     @testset "1dim" begin
         a = [Point(i, rand()) for i in 1:n1]
-        p₊ = (1,)
+        p₊ = (Val(1),)
         k₊ = (KnotVector(4.5),)
         # k₊ = (KnotVector(4.5,4.95),)
 
         @testset "BSplineManifold" begin
             w = ones(n1)
             M = BSplineManifold(a, (P1,))
-            M0 = refinement(M)
-            M1 = refinement(M, k₊=k₊)
-            M2 = refinement(M, p₊=p₊)
-            M3 = refinement(M, p₊=p₊, k₊=k₊)
-            M4 = refinement(M, (P1′,))
+            M0 = @inferred refinement(M)
+            M1 = @inferred refinement(M, k₊)
+            M2 = @inferred refinement(M, p₊)
+            M3 = @inferred refinement(M, p₊, k₊)
+            M4 = @inferred refinement(M, (P1′,))
             R = RationalBSplineManifold(a, w, (P1,))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,))
+            R0 = @inferred refinement(R)
+            R1 = @inferred refinement(R, k₊)
+            R2 = @inferred refinement(R, p₊)
+            R3 = @inferred refinement(R, p₊, k₊)
+            R4 = @inferred refinement(R, (P1′,))
 
             for _ in 1:100
                 t1 = rand(D1)
@@ -62,11 +62,11 @@
         @testset "RationalBSplineManifold" begin
             w = rand(n1) .+ 1
             R = RationalBSplineManifold(a, w, (P1,))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,))
+            R0 = @inferred refinement(R)
+            R1 = @inferred refinement(R, k₊)
+            R2 = @inferred refinement(R, p₊)
+            R3 = @inferred refinement(R, p₊, k₊)
+            R4 = @inferred refinement(R, (P1′,))
 
             for _ in 1:100
                 t1 = rand(D1)
@@ -80,24 +80,24 @@
 
     @testset "2dim" begin
         a = [Point(i, rand()) for i in 1:n1, j in 1:n2]
-        p₊ = (1,2)
+        p₊ = (Val(1), Val(2))
         k₊ = (KnotVector(4.5,4.7),KnotVector())
 
         @testset "BSplineManifold" begin
             w = ones(n1,n2)
 
             M = BSplineManifold(a, (P1,P2))
-            M0 = refinement(M)
-            M1 = refinement(M, k₊=k₊)
-            M2 = refinement(M, p₊=p₊)
-            M3 = refinement(M, p₊=p₊, k₊=k₊)
-            M4 = refinement(M, (P1′,P2′))
+            M0 = @inferred refinement(M)
+            M1 = @inferred refinement(M, k₊)
+            M2 = @inferred refinement(M, p₊)
+            M3 = @inferred refinement(M, p₊, k₊)
+            M4 = @inferred refinement(M, (P1′,P2′))
             R = RationalBSplineManifold(a, w, (P1,P2))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,P2′))
+            R0 = @inferred refinement(R)
+            R1 = @inferred refinement(R, k₊)
+            R2 = @inferred refinement(R, p₊)
+            R3 = @inferred refinement(R, p₊, k₊)
+            R4 = @inferred refinement(R, (P1′,P2′))
 
             for _ in 1:100
                 t1 = rand(D1)
@@ -118,11 +118,11 @@
             w = rand(n1,n2) .+ 1
 
             R = RationalBSplineManifold(a, w, (P1,P2))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,P2′))
+            R0 = @inferred refinement(R)
+            R1 = @inferred refinement(R, k₊)
+            R2 = @inferred refinement(R, p₊)
+            R3 = @inferred refinement(R, p₊, k₊)
+            R4 = @inferred refinement(R, (P1′,P2′))
 
             for _ in 1:100
                 t1 = rand(D1)
@@ -137,24 +137,38 @@
 
     @testset "3dim" begin
         a = [Point(i, rand()) for i in 1:n1, j in 1:n2, k in 1:n3]
-        p₊ = (1,2,0)
+        p₊ = (Val(1), Val(2), Val(0))
         k₊ = (KnotVector(4.4,4.7),KnotVector(),KnotVector(42))
 
         @testset "BSplineManifold" begin
             w = ones(n1,n2,n3)
 
             M = BSplineManifold(a, (P1,P2,P3))
-            M0 = refinement(M)
-            M1 = refinement(M, k₊=k₊)
-            M2 = refinement(M, p₊=p₊)
-            M3 = refinement(M, p₊=p₊, k₊=k₊)
-            M4 = refinement(M, (P1′,P2′,P3′))
             R = RationalBSplineManifold(a, w, (P1,P2,P3))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,P2′,P3′))
+            # On Julia v1.6, the following script seems not type-stable.
+            if VERSION ≥ v"1.8"
+                M0 = @inferred refinement(M)
+                M1 = @inferred refinement(M, k₊)
+                M2 = @inferred refinement(M, p₊)
+                M3 = @inferred refinement(M, p₊, k₊)
+                M4 = @inferred refinement(M, (P1′,P2′,P3′))
+                R0 = @inferred refinement(R)
+                R1 = @inferred refinement(R, k₊)
+                R2 = @inferred refinement(R, p₊)
+                R3 = @inferred refinement(R, p₊, k₊)
+                R4 = @inferred refinement(R, (P1′,P2′,P3′))
+            else
+                M0 = refinement(M)
+                M1 = refinement(M, k₊)
+                M2 = refinement(M, p₊)
+                M3 = refinement(M, p₊, k₊)
+                M4 = refinement(M, (P1′,P2′,P3′))
+                R0 = refinement(R)
+                R1 = refinement(R, k₊)
+                R2 = refinement(R, p₊)
+                R3 = refinement(R, p₊, k₊)
+                R4 = refinement(R, (P1′,P2′,P3′))
+            end
 
             for _ in 1:100
                 t1 = rand(D1)
@@ -176,11 +190,20 @@
             w = rand(n1,n2,n3) .+ 1
 
             R = RationalBSplineManifold(a, w, (P1,P2,P3))
-            R0 = refinement(R)
-            R1 = refinement(R, k₊=k₊)
-            R2 = refinement(R, p₊=p₊)
-            R3 = refinement(R, p₊=p₊, k₊=k₊)
-            R4 = refinement(R, (P1′,P2′,P3′))
+            # On Julia v1.6, the following script seems not type-stable.
+            if VERSION ≥ v"1.8"
+                R0 = @inferred refinement(R)
+                R1 = @inferred refinement(R, k₊)
+                R2 = @inferred refinement(R, p₊)
+                R3 = @inferred refinement(R, p₊, k₊)
+                R4 = @inferred refinement(R, (P1′,P2′,P3′))
+            else
+                R0 = refinement(R)
+                R1 = refinement(R, k₊)
+                R2 = refinement(R, p₊)
+                R3 = refinement(R, p₊, k₊)
+                R4 = refinement(R, (P1′,P2′,P3′))
+            end
 
             for _ in 1:100
                 t1 = rand(D1)