Repository: hyrodium/BasicBSpline.jl
Pull request #82 (2021-01-26 14:55:17)

Add jldoctest


Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d7fb58bb..0405b638 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -61,8 +61,9 @@ jobs:
             Pkg.instantiate()'
       - run: |
           julia --project=docs -e '
-            using Documenter: doctest
+            using Documenter
             using BasicBSpline
+            DocMeta.setdocmeta!(BasicBSpline, :DocTestSetup, :(using LinearAlgebra, BasicBSpline))
             doctest(BasicBSpline)'
       - run: julia --project=docs docs/make.jl
         env:
diff --git a/Project.toml b/Project.toml
index 216dac9d..d85262e2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BasicBSpline"
 uuid = "4c5d9882-2acf-4ea4-9e48-968fd4518195"
 authors = ["hyrodium <hyrodium@gmail.com>"]
-version = "0.1.5"
+version = "0.1.6"
 
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
diff --git a/docs/make.jl b/docs/make.jl
index da1b4cef..a20defc5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,5 +1,8 @@
 using Documenter, BasicBSpline
 
+# Setup for doctests in docstrings
+DocMeta.setdocmeta!(BasicBSpline, :DocTestSetup, :(using LinearAlgebra, BasicBSpline))
+
 function generate_indexmd_from_readmemd()
     path_readme = "README.md"
     path_index = "docs/src/index.md"
diff --git a/docs/src/detail.md b/docs/src/detail.md
index 62c35266..7c89d98a 100644
--- a/docs/src/detail.md
+++ b/docs/src/detail.md
@@ -1,6 +1,3 @@
 # Docstrings
 
 ## All docstrings
-```@autodocs
-Modules = [BasicBSpline]
-```
diff --git a/docs/src/math.md b/docs/src/math.md
index 04e5d19c..ae1fcf87 100644
--- a/docs/src/math.md
+++ b/docs/src/math.md
@@ -6,7 +6,7 @@
 In this page, we'll explain the mathematical definition and property of B-spline with Julia code.
 
 Before running the following code, do not forget importing the package:
-```julia
+```@example
 using BasicBSpline
 ```
 
@@ -26,67 +26,39 @@ using BasicBSpline
 
 [fig]
 
-```julia
-k = Knots([1,2,3])
-k = Knots(1:3)
-k = Knots(1,2,3)
+```@docs
+Knots(vector::AbstractVector{<:Real})
 ```
 
-We denotes a number of knots by *sharp* symbol like this:
-```math
-\sharp k = \sharp(k_1, \dots, k_l) =l
+```@docs
+Knots(knot::Real)
 ```
 
-```julia
-k = Knots([4,5,6])
-♯(k) # 3
-length(k) # 3
+```@docs
+length(k::Knots)
 ```
 
-We introduce additional operator ``+`` and product operator ``\cdot``
-```math
-\begin{aligned}
-k^{(1)}+k^{(2)}
-&=(k^{(1)}_1, \dots, k^{(1)}_l) + (k^{(2)}_1, \dots, k^{(2)}_l) \\
-&=(\text{sort of union of} \  k^{(1)} \ \text{and} \  k^{(2)} \text{)} \\
-m\cdot k&=\underbrace{k+\cdots+k}_{m}
-\end{aligned}
-```
-For example, ``(1,2,3)+(2,4,5)=(1,2,2,3,5)``, ``2\cdot (2,3)=(2,2,3,3)``.
+We introduce **additional operator** ``+``.
 
-```julia
-Knots([1,2,3]) + Knots([2,4,5]) # Knots([1,2,2,3,5])
-2 * Knots([2,3]) # Knots([2,2,3,3])
+```@docs
++(k₁::Knots, k₂::Knots)
 ```
 
-Deleting duplicates operator
-```math
-\begin{aligned}
-\widehat{k}
-&=(\text{remove duplicates of} \  k) \\
-\end{aligned}
-```
-For example, ``\widehat{(1,2,2,3)}=(1,2,3)``.
+And **product operator** ``\cdot``.
 
-```julia
-unique(Knots([1,2,2,3])) # Knots([1,2,3])
+```@docs
+*(n::Integer, k::Knots)
 ```
 
-For Given knot vector ``k``, the following function ``\mathfrak{n}_k:\mathbb{R}\to\mathbb{Z}`` represents the counts of ..
-
-```math
-\mathfrak{n}_k(t) = \sharp\{i \mid k_i=t \}
+```@docs
+unique(k::Knots)
 ```
-For example, if ``k=(1,2,2,3)``, then ``\mathfrak{n}_k(0.3)=0``, ``\mathfrak{n}_k(1)=1``, ``\mathfrak{n}_k(2)=2``.
 
-
-```julia
-k = Knots([1,2,2,3])
-𝔫(k,0.3) # 0
-𝔫(k,1.0) # 1
-𝔫(k,2.0) # 2
+```@docs
+𝔫(k::Knots, t::Real)
 ```
 
+
 ## B-spline space
 Before defining B-spline space, we'll define polynomial space with degree ``p``.
 
@@ -130,7 +102,8 @@ Note that a element of the space ``\mathcal{P}[p,k]`` is piecewise polynomial.
 
 [fig]
 
-```julia
+```@repl
+using BasicBSpline
 p = 2
 k = Knots([1,3,5,6,8,9])
 BSplineSpace(p,k)
@@ -143,7 +116,8 @@ k_{i}&<k_{i+p+1} & (1 \le i \le l-p-1)
 \end{aligned}
 ```
 
-```julia
+```@repl
+using BasicBSpline
 isproper(BSplineSpace(2,Knots([1,3,5,6,8,9]))) # true
 isproper(BSplineSpace(1,Knots([1,3,3,3,8,9]))) # false
 ```
@@ -153,7 +127,8 @@ The B-spline space is linear space, and if a B-spline space is proper, its dimen
 \dim(\mathcal{P}[p,k])=\sharp k -p -1
 ```
 
-```julia
+```@repl
+using BasicBSpline
 dim(BSplineSpace(2,Knots([1,3,5,6,8,9]))) # 3
 ```
 
@@ -182,7 +157,8 @@ dim(BSplineSpace(2,Knots([1,3,5,6,8,9]))) # 3
     The set of functions ``\{B_{(i,p,k)}\}_i`` is a basis of B-spline space ``\mathcal{P}[p,k]``.
 
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -206,7 +182,8 @@ You can choose the first terms in different ways.
 \end{aligned}
 ```
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -226,7 +203,8 @@ plot([t->bsplinebasis₋₀(i,P,t) for i in 1:dim(P)], 1, 8, ylims=(0,1.05))
 
 [fig]
 
-```julia
+```@repl
+using BasicBSpline
 i = 2
 k = Knots([5,12,13,13,14])
 p = 2
@@ -247,7 +225,8 @@ bsplinesupport(i,P) # 12..14
     ```
     Note that ``\dot{B}_{(i,p,k)}\in\mathcal{P}[p-1,k]``.
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -267,7 +246,8 @@ plot([t->bsplinebasis′₊₀(i,P,t) for i in 1:dim(P)], 1, 8, ylims=(0,1.05))
     \end{aligned}
     ```
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -280,7 +260,8 @@ plot(t->sum(bsplinebasis₊₀(i,P,t) for i in 1:dim(P)), 1, 8, ylims=(0,1.05))
 
 To satisfy the partition of unity on whole interval ``[1,8]``, sometimes more knots will be inserted to the endpoints of the interval.
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -306,7 +287,8 @@ Therefore, to satisfy partition of unity on closed interval ``[k_{p+1}, k_{l-p}]
 \end{aligned}
 ```
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 p = 2
@@ -329,7 +311,8 @@ plot(t->sum(bsplinebasis(i,P,t) for i in 1:dim(P)), 1, 8, ylims=(0,1.05))
 
 (as linear subspace..)
 
-```julia
+```@repl
+using BasicBSpline
 P1 = BSplineSpace(1,Knots([1,3,5,8]))
 P2 = BSplineSpace(1,Knots([1,3,5,6,8,9]))
 P3 = BSplineSpace(2,Knots([1,1,3,3,5,5,8,8]))
@@ -341,7 +324,8 @@ P2 ⊈ P3 # true
 
 Here are plots of the B-spline basis functions of the spaces `P1`, `P2`, `P3`.
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 P1 = BSplineSpace(1,Knots([1,3,5,8]))
@@ -371,12 +355,14 @@ n'&=\dim(\mathcal{P}[p',k'])
 
 You can calculate the change of basis matrix ``A`` with `changebasis`.
 
-```julia
+```@repl
+using BasicBSpline
 A12 = changebasis(P1,P2)
 A13 = changebasis(P1,P3)
 ```
 
-```julia
+```@repl
+using BasicBSpline
 using Plots
 gr()
 plot(
@@ -412,7 +398,8 @@ B-spline manifold is a parametric representation of a shape.
 
 We will also write ``\bm{p}(t^1,\dots,t^d; \bm{a})``, ``\bm{p}(t^1,\dots,t^d)``, ``\bm{p}(t; \bm{a})`` or ``\bm{p}(t)`` for simplicity.
 
-```julia
+```@repl
+using BasicBSpline
 P1 = BSplineSpace(1,Knots([0,0,1,1]))
 P2 = BSplineSpace(1,Knots([1,1,2,3,3]))
 n1 = dim(P1) # 2
@@ -423,7 +410,8 @@ M = BSplineManifold([P1, P2], a)
 
 
 ### B-spline curve
-```julia
+```@repl
+using BasicBSpline
 ## 1-dim B-spline manifold
 p = 2 # degree of polynomial
 k = Knots(1:12) # knot vector
@@ -436,7 +424,8 @@ save_png("1dim.png", M, unitlength = 50)
 
 
 ### B-spline surface
-```julia
+```@repl
+using BasicBSpline
 p = 2 # degree of polynomial
 k = Knots(1:8) # knot vector
 P = BSplineSpace(p,k) # B-spline space
@@ -447,7 +436,8 @@ save_png("2dim.png", M) # save image
 ```
 ![](img/2dim.png)
 
-```julia
+```@repl
+using BasicBSpline
 points = [M([u,v]) for u in range(3.0,6.0,length=50), v in range(3.0,6.0,length=50)]
 X = [point[1] for point in points]
 Y = [point[2] for point in points]
@@ -473,7 +463,8 @@ save("2dim_makie.png", scene)
 ### h-refinemnet
 Insert additional knots to knot vector.
 
-```julia
+```@repl
+using BasicBSpline
 k₊=[Knots(3.3,4.2),Knots(3.8,3.2,5.3)] # additional knots
 M′ = refinement(M,k₊=k₊) # refinement of B-spline manifold
 save_png("2dim_h-refinement.png", M′) # save image
@@ -486,7 +477,8 @@ Note that this shape and the last shape are identical.
 ### p-refinemnet
 Increase the polynomial degree of B-spline manifold.
 
-```julia
+```@repl
+using BasicBSpline
 p₊=[1,2] # additional degrees
 M′ = refinement(M,p₊=p₊) # refinement of B-spline manifold
 save_png("2dim_p-refinement.png", M′) # save image
@@ -500,7 +492,8 @@ Note that this shape and the last shape are identical.
 Least squares method.
 
 [Try on Desmos graphing graphing calculator!](https://www.desmos.com/calculator/2hm3b1fbdf)
-```julia
+```@repl
+using BasicBSpline
 p1 = 2
 p2 = 2
 k1 = Knots(-10:10)+p1*Knots(-10,10)
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index 4b1adc1a..b4bd12bc 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -327,20 +327,30 @@ end
 """
 * func: Real -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace)
-    b = innerproduct_I(func, P1)
-    A = innerproduct_I(P1)
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace; domain=:I)
+    if domain == :I
+        b = innerproduct_I(func, P1)
+        A = innerproduct_I(P1)
+    elseif domain == :R
+        b = innerproduct_R(func, P1)
+        A = innerproduct_R(P1)
+    end
     return inv(A) * b
 end
 
 """
 * func: (Real,Real) -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace; domain=:I)
     n1, n2 = dim(P1), dim(P2)
-    A1, A2 = innerproduct_I(P1), innerproduct_I(P2)
+    if domain == :I
+        A1, A2 = innerproduct_I(P1), innerproduct_I(P2)
+        b = innerproduct_I(func, P1, P2)
+    elseif domain == :R
+        A1, A2 = innerproduct_R(P1), innerproduct_R(P2)
+        b = innerproduct_R(func, P1, P2)
+    end
     A = [A1[i1, j1] * A2[i2, j2] for i1 in 1:n1, i2 in 1:n2, j1 in 1:n1, j2 in 1:n2]
-    b = innerproduct_I(func, P1, P2)
     _A = reshape(A, n1 * n2, n1 * n2)
     _b = reshape(b, n1 * n2)
     return reshape(inv(_A) * _b, n1, n2)
@@ -349,11 +359,16 @@ end
 """
 * func: (Real,Real,Real) -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace)
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace; domain=:I)
     n1, n2, n3 = dim(P1), dim(P2), dim(P3)
-    A1, A2, A3 = innerproduct_I(P1), innerproduct_I(P2), innerproduct_I(P3)
+    if domain == :I
+        A1, A2, A3 = innerproduct_I(P1), innerproduct_I(P2), innerproduct_I(P3)
+        b = innerproduct_I(func, P1, P2, P3)
+    elseif domain == :R
+        A1, A2, A3 = innerproduct_R(P1), innerproduct_R(P2), innerproduct_R(P3)
+        b = innerproduct_R(func, P1, P2, P3)
+    end
     A = [A1[i1, j1] * A2[i2, j2] * A3[i3, j3] for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3, j1 in 1:n1, j2 in 1:n2, j3 in 1:n3]
-    b = innerproduct_I(func, P1, P2, P3)
     _A = reshape(A, n1 * n2 * n3, n1 * n2 * n3)
     _b = reshape(b, n1 * n2 * n3)
     return reshape(inv(_A) * _b, n1, n2, n3)
@@ -365,8 +380,8 @@ Approximate given function by linear combination of B-spline functions.
 This function returns its control points.
 * func: Vector{<:Real} -> ℝ-vector space
 """
-function fittingcontrolpoints(func, Ps::AbstractVector{<:AbstractBSplineSpace})
+function fittingcontrolpoints(func, Ps::AbstractVector{<:AbstractBSplineSpace}; domain=:I)
     # TODO: currently, this function only supports for 1-dim, 2-dim and 3-dim B-spline manifold.
     _func(t...) = func([t...])
-    return fittingcontrolpoints(_func, Ps...)
+    return fittingcontrolpoints(_func, Ps..., domain=domain)
 end
diff --git a/src/_Knots.jl b/src/_Knots.jl
index 71f4d1e3..9d14aed2 100644
--- a/src/_Knots.jl
+++ b/src/_Knots.jl
@@ -1,31 +1,108 @@
 # Knots
 
-@doc raw"""
-Construct knot vector from given array.
-```math
-k=(k_1,\dots,k_l)
-```
-"""
 struct Knots
     vector::Vector{Float64}
+    @doc raw"""
+    Construct knot vector from given array.
+    ```math
+    k=(k_1,\dots,k_l)
+    ```
+    
+    # Examples
+    ```jldoctest
+    julia> k = Knots([1,2,3])
+    Knots([1.0, 2.0, 3.0])
+    
+    julia> k = Knots(1:3)
+    Knots([1.0, 2.0, 3.0])
+    ```
+    """
     function Knots(vector::AbstractVector{<:Real})
         return new(sort(convert(Vector{Float64}, vector)))
     end
 end
+
 function Knots(vector::AbstractVector)
     if isempty(vector)
+        @warn "use Knots() or Knots(Float64[]), instead of Knots([])."
         return Knots(Float64[])
     else
+        @warn "the type $(typeof(vector)) is deprecated."
         return Knots(convert(Vector{Float64}, vector))
     end
 end
+
+@doc raw"""
+Construct knot vector from given real numbers.
+
+# Examples
+```jldoctest
+julia> k = Knots(1,2,3)
+Knots([1.0, 2.0, 3.0])
+
+julia> k = Knots()
+Knots([])
+```
+"""
 function Knots(knot::Real...)
     return Knots(collect(knot))
 end
 
+function Base.show(io::IO, k::Knots)
+    if k.vector == Float64[]
+        print(io, "Knots([])")
+    else
+        print(io, "Knots($(k.vector))")
+    end
+end
+
 Base.zero(::Type{Knots}) = Knots(Float64[])
 Base.:(==)(k₁::Knots, k₂::Knots) = (k₁.vector == k₂.vector)
+
+@doc raw"""
+Sum of knot vectors
+
+```math
+\begin{aligned}
+k^{(1)}+k^{(2)}
+&=(k^{(1)}_1, \dots, k^{(1)}_{l^{(1)}}) + (k^{(2)}_1, \dots, k^{(2)}_{l^{(2)}}) \\
+&=(\text{sort of union of} \  k^{(1)} \ \text{and} \  k^{(2)} \text{)}
+\end{aligned}
+```
+
+For example, ``(1,2,3,5)+(4,5,8)=(1,2,3,4,5,5,8)``.
+
+# Examples
+```jldoctest
+julia> k1 = Knots(1,2,3,5);
+
+julia> k2 = Knots(4,5,8);
+
+julia> k1 + k2
+Knots([1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 8.0])
+```
+"""
 Base.:+(k₁::Knots, k₂::Knots) = Knots(sort([k₁.vector..., k₂.vector...]))
+
+@doc raw"""
+Product of integer and knot vector
+
+```math
+\begin{aligned}
+m\cdot k&=\underbrace{k+\cdots+k}_{m}
+\end{aligned}
+```
+
+For example, ``2\cdot (1,2,2,5)=(1,1,2,2,2,2,5,5)``.
+
+# Examples
+```jldoctest
+julia> k = Knots(1,2,2,5);
+
+julia> 2 * k
+Knots([1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 5.0, 5.0])
+```
+"""
 function Base.:*(p₊::Integer, k::Knots)
     if p₊ == 0
         return zero(Knots)
@@ -39,10 +116,47 @@ end
 Base.in(r::Real, k::Knots) = in(r, k.vector)
 Base.getindex(k::Knots, i::Integer) = k.vector[i]
 Base.getindex(k::Knots, v::AbstractVector{<:Integer}) = Knots(k.vector[v])
+
+@doc raw"""
+Length of knot vector
+
+# Examples
+```jldoctest
+julia> k = Knots(1,2,3,5);
+
+julia> length(k)
+4
+```
+"""
 Base.length(k::Knots) = length(k.vector)
-♯(k::Knots) = length(k::Knots)
+
+function ♯(k::Knots)
+    @warn "BasicBSpline.♯ is deprecated."
+    return length(k::Knots)
+end
 Base.firstindex(k::Knots) = 1
 Base.lastindex(k::Knots) = length(k)
+
+@doc raw"""
+Unique elements of knot vector.
+
+```math
+\begin{aligned}
+\widehat{k}
+&=(\text{remove duplicates of} \  k) \\
+\end{aligned}
+```
+
+For example, ``\widehat{(1,2,2,3)}=(1,2,3)``.
+
+# Examples
+```jldoctest
+julia> k = Knots([1,2,2,3]);
+
+julia> unique(k)
+Knots([1.0, 2.0, 3.0])
+```
+"""
 Base.unique(k::Knots) = Knots(unique(k.vector))
 Base.iterate(k::Knots) = iterate(k.vector)
 Base.iterate(k::Knots, i::Integer) = iterate(k.vector, i)
@@ -66,6 +180,27 @@ function Base.issubset(k::Knots, k′::Knots)
     return true
 end
 
+@doc raw"""
+For Given knot vector ``k``, the following function ``\mathfrak{n}_k:\mathbb{R}\to\mathbb{Z}`` represents the number of knots that duplicate the knot vector ``k``.
+
+```math
+\mathfrak{n}_k(t) = \sharp\{i \mid k_i=t \}
+```
+For example, if ``k=(1,2,2,3)``, then ``\mathfrak{n}_k(0.3)=0``, ``\mathfrak{n}_k(1)=1``, ``\mathfrak{n}_k(2)=2``.
+
+```jldoctest
+julia> k = Knots([1,2,2,3]);
+
+julia> 𝔫(k,0.3)
+0
+
+julia> 𝔫(k,1.0)
+1
+
+julia> 𝔫(k,2.0)
+2
+```
+"""
 𝔫(k::Knots, t::Real) = count(==(t), k.vector)
 
 """


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 9bcfee67..efc34b50 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,24 +9,51 @@ using Random
 @testset "BasicBSpline.jl" begin
 
     @testset "Knots" begin
-        @test Knots([1, 2, 3]) == Knots(1:3) == Knots(1, 2, 3)
-        @test Knots(2, 4, 5) == Knots([2, 4, 5])
-        @test zero(Knots) == Knots([])
-        @test Knots(1:3) == Knots([3, 2, 1])
+        @test Knots([1,2,3]) == Knots(1:3) == Knots(1,2,3)
+        @test Knots(2,4,5) == Knots([2,4,5])
+        @test zero(Knots) == Knots() == Knots([])
+        @test Knots(1:3) == Knots([3,2,1])
 
-        k = Knots([4, 5, 6])
-        @test ♯(k) == length(k) == 3
+        k = Knots([1,2,2,3])
+        @test ♯(k) == length(k) == 4
 
-        @test Knots([-1, 2, 3]) + 2 * Knots([2, 5]) == Knots([-1, 2, 2, 2, 3, 5, 5])
-        Knots([1, 2, 3]) + Knots([2, 4, 5]) == Knots([1, 2, 2, 3, 5])
-        2 * Knots([2, 3]) == Knots([2, 2, 3, 3])
+        @test Knots([-1,2,3]) + 2 * Knots([2,5]) == Knots([-1,2,2,2,3,5,5])
+        @test Knots([1,2,3]) + Knots([2,4,5]) == Knots([1,2,2,3,4,5])
+        @test 2 * Knots([2,3]) == Knots([2,2,3,3])
 
-        unique(Knots([1, 2, 2, 3])) == Knots([1, 2, 3])
+        @test unique(Knots([1,2,2,3])) == Knots([1,2,3])
+
+        @test Knots([1,2,3]) ⊆ Knots([1,2,3])
+        @test Knots([1,2,3]) ⊇ Knots([1,2,3])
+        @test Knots([1,2,2,3]) ⊆ Knots([1,2,2,3,5])
+        @test Knots([1,2,2,3]) ⊈ Knots([1,2,3,5])
+        @test Knots([1,2,2,3,5]) ⊇ Knots([1,2,2,3])
+        @test Knots([1,2,3,5]) ⊉ Knots([1,2,2,3])
 
-        k = Knots([1, 2, 2, 3])
         @test 𝔫(k, 0.3) == 0
         @test 𝔫(k, 1.0) == 1
         @test 𝔫(k, 2.0) == 2
+        @test 1 ∈ k
+        @test 1.5 ∉ k
+        @test BasicBSpline._knotindex(k,1.5) == 1
+        @test BasicBSpline._knotindex(k,2.5) == 3
+        @test BasicBSpline._knotindex₋₀(k,1.5) == 1
+        @test BasicBSpline._knotindex₋₀(k,2.5) == 3
+        @test BasicBSpline._knotindex₊₀(k,1.5) == 1
+        @test BasicBSpline._knotindex₊₀(k,2.5) == 3
+
+        @test BasicBSpline._knotindex(k,1) == 1
+        @test BasicBSpline._knotindex(k,2) == 3
+        @test BasicBSpline._knotindex(k,3) == 3
+        @test BasicBSpline._knotindex₋₀(k,1) == nothing
+        @test BasicBSpline._knotindex₋₀(k,2) == 1
+        @test BasicBSpline._knotindex₋₀(k,3) == 3
+        @test BasicBSpline._knotindex₊₀(k,1) == 1
+        @test BasicBSpline._knotindex₊₀(k,2) == 3
+        @test BasicBSpline._knotindex₊₀(k,3) == nothing
+
+        @test string(k) == "Knots([1.0, 2.0, 2.0, 3.0])"
+        @test string(Knots()) == "Knots([])"
     end
 
     @testset "BSplineSpace" begin
@@ -329,6 +356,25 @@ using Random
         @test prod([M(t) ≈ M′′(t) for t in ts])
     end
 
+    @testset "array <-> array of vector" begin
+        Random.seed!(42)
+
+        A = rand(3,4)
+        A′ = BasicBSpline.array2arrayofvector(A)
+        A′′ = BasicBSpline.arrayofvector2array(A′)
+        @test A == A′′
+
+        B = rand(3,4,5)
+        B′ = BasicBSpline.array2arrayofvector(B)
+        B′′ = BasicBSpline.arrayofvector2array(B′)
+        @test B == B′′
+
+        C = rand(3,4,5,6)
+        C′ = BasicBSpline.array2arrayofvector(C)
+        C′′ = BasicBSpline.arrayofvector2array(C′)
+        @test C == C′′
+    end
+
     @testset "Fitting-curve_R" begin
         Random.seed!(42)
 
@@ -348,7 +394,10 @@ using Random
 
         a_tmp = fittingcontrolpoints(M, [P1′])
         a_fit = transpose(hcat(a_tmp...))
+        @test norm(a_fit - a_ref) < ε
 
+        a_tmp = fittingcontrolpoints(M, [P1′], domain=:R)
+        a_fit = transpose(hcat(a_tmp...))
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -375,25 +424,6 @@ using Random
         @test norm(a_fit - a_ref) < ε
     end
 
-    @testset "array <-> array of vector" begin
-        Random.seed!(42)
-
-        A = rand(3,4)
-        A′ = BasicBSpline.array2arrayofvector(A)
-        A′′ = BasicBSpline.arrayofvector2array(A′)
-        @test A == A′′
-
-        B = rand(3,4,5)
-        B′ = BasicBSpline.array2arrayofvector(B)
-        B′′ = BasicBSpline.arrayofvector2array(B′)
-        @test B == B′′
-
-        C = rand(3,4,5,6)
-        C′ = BasicBSpline.array2arrayofvector(C)
-        C′′ = BasicBSpline.arrayofvector2array(C′)
-        @test C == C′′
-    end
-
     @testset "Fitting-surface_R" begin
         Random.seed!(42)
 
@@ -420,7 +450,10 @@ using Random
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
         a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        @test norm(a_fit - a_ref) < ε
 
+        a_tmp = fittingcontrolpoints(M, [P1′, P2′], domain=:R)
+        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -450,7 +483,6 @@ using Random
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
         a_fit = BasicBSpline.arrayofvector2array(a_tmp)
-
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -487,7 +519,10 @@ using Random
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′])
         a_fit = BasicBSpline.arrayofvector2array(a_tmp)
+        @test norm(a_fit - a_ref) < ε
 
+        a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′], domain=:R)
+        a_fit = BasicBSpline.arrayofvector2array(a_tmp)
         @test norm(a_fit - a_ref) < ε
     end