Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
Add support `M(:,0.2,:)`
If `M isa BSplineManifold{3}`, then fixing 2nd value will create 2-dimensional manifold.

Pull request:
This PR fixes #216.

Patch:
diff --git a/docs/src/math-bsplinemanifold.md b/docs/src/math-bsplinemanifold.md
index 0a505ca23..b15489feb 100644
--- a/docs/src/math-bsplinemanifold.md
+++ b/docs/src/math-bsplinemanifold.md
@@ -86,8 +86,6 @@ savefig("2dim-manifold.html") # hide
 <object type="text/html" data="../2dim-manifold.html" style="width:100%;height:420px;"></object>
 ```
 
-![](2dim.png)
-
 ## Affine commutativity
 !!! info "Thm.  Affine commutativity"
     If ``T`` is a affine transform ``\mathbb{R}^d\to\mathbb{R}^d``, then the following equality holds.
@@ -96,3 +94,32 @@ savefig("2dim-manifold.html") # hide
     =\bm{p}(t; T(\bm{a}))
     ```
 
+## Currying
+
+```@repl math
+p = 2; # degree of polynomial
+k = KnotVector(1:8); # knot vector
+P = BSplineSpace{p}(k); # B-spline space
+rand_a = [SVector(rand(), rand(), rand()) for i in 1:dim(P), j in 1:dim(P)];
+a = [SVector(2*i-6.5, 2*j-6.5, 0) for i in 1:dim(P), j in 1:dim(P)] + rand_a; # random generated control points
+M = BSplineManifold(a,(P,P)); # Define B-spline manifold
+M isa BSplineManifold{2}
+M(:,:) isa BSplineManifold{2}
+M(4.3,:) isa BSplineManifold{1}
+```
+
+```@example math
+p = 2 # degree of polynomial
+k = KnotVector(1:8) # knot vector
+P = BSplineSpace{p}(k) # B-spline space
+rand_a = [SVector(rand(), rand(), rand()) for i in 1:dim(P), j in 1:dim(P)]
+a = [SVector(2*i-6.5, 2*j-6.5, 0) for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
+M = BSplineManifold(a,(P,P)) # Define B-spline manifold
+plot(M)
+plot!(M(4.3,:), linewidth = 5, color=:white)
+savefig("2dim-manifold-currying.html") # hide
+```
+
+```@raw html
+<object type="text/html" data="../2dim-manifold-currying.html" style="width:100%;height:420px;"></object>
+```
diff --git a/src/_BSplineManifold.jl b/src/_BSplineManifold.jl
index e1689c3a4..6180b7a1d 100644
--- a/src/_BSplineManifold.jl
+++ b/src/_BSplineManifold.jl
@@ -51,6 +51,8 @@ struct BSplineManifold{Dim,Deg,C,S<:Tuple} <: AbstractBSplineManifold{Dim,Deg}
     end
 end
 
+Base.:(==)(M1::AbstractBSplineManifold, M2::AbstractBSplineManifold) = (bsplinespaces(M1)==bsplinespaces(M2)) & (controlpoints(M1)==controlpoints(M2))
+
 bsplinespaces(M::BSplineManifold) = M.bsplinespaces
 controlpoints(M::BSplineManifold) = M.controlpoints
 
@@ -109,20 +111,20 @@ end
     )
 end
 
-@inline function (M::AbstractManifold{1})(t1)
+@inline function (M::AbstractManifold{1})(t1::Real)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     unsafe_mapping(M,t1)
 end
 
-@inline function (M::AbstractManifold{2})(t1,t2)
+@inline function (M::AbstractManifold{2})(t1::Real,t2::Real)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     t2 in domain(Ps[2]) || throw(DomainError(t2, "The input $(t2) is out of range."))
     unsafe_mapping(M,t1,t2)
 end
 
-@inline function (M::AbstractManifold{3})(t1,t2,t3)
+@inline function (M::AbstractManifold{3})(t1::Real,t2::Real,t3::Real)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     t2 in domain(Ps[2]) || throw(DomainError(t2, "The input $(t2) is out of range."))
@@ -130,4 +132,105 @@ end
     unsafe_mapping(M,t1,t2,t3)
 end
 
+
+## currying
+# 1dim
+@inline function (M::AbstractBSplineManifold{1})(::Colon)
+    a = copy(controlpoints(M))
+    Ps = bsplinespaces(M)
+    return BSplineManifold(a,Ps)
+end
+
+# 2dim
+@inline function (M::AbstractBSplineManifold{2})(::Colon,::Colon)
+    a = copy(controlpoints(M))
+    Ps = bsplinespaces(M)
+    return BSplineManifold(a,Ps)
+end
+@inline function (M::AbstractBSplineManifold{2,p})(t1::Real,::Colon) where p
+    p1, p2 = p
+    P1, P2 = bsplinespaces(M)
+    a = controlpoints(M)
+    j1 = intervalindex(P1,t1)
+    B1 = bsplinebasisall(P1,j1,t1)
+    b = sum(a[j1+i1,:]*B1[1+i1] for i1 in 0:p1)
+    return BSplineManifold(b,(P2,))
+end
+@inline function (M::AbstractBSplineManifold{2,p})(::Colon,t2::Real) where p
+    p1, p2 = p
+    P1, P2 = bsplinespaces(M)
+    a = controlpoints(M)
+    j2 = intervalindex(P2,t2)
+    B2 = bsplinebasisall(P2,j2,t2)
+    b = sum(a[:,j2+i2]*B2[1+i2] for i2 in 0:p2)
+    return BSplineManifold(b,(P1,))
+end
+
+# 3dim
+@inline function (M::AbstractBSplineManifold{3})(::Colon,::Colon,::Colon)
+    a = copy(controlpoints(M))
+    Ps = bsplinespaces(M)
+    return BSplineManifold(a,Ps)
+end
+@inline function (M::AbstractBSplineManifold{3,p})(t1::Real,::Colon,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j1 = intervalindex(P1,t1)
+    B1 = bsplinebasisall(P1,j1,t1)
+    b = sum(a[j1+i1,:,:]*B1[1+i1] for i1 in 0:p1)
+    return BSplineManifold(b,(P2,P3))
+end
+@inline function (M::AbstractBSplineManifold{3,p})(::Colon,t2::Real,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j2 = intervalindex(P2,t2)
+    B2 = bsplinebasisall(P2,j2,t2)
+    b = sum(a[:,j2+i2,:]*B2[1+i2] for i2 in 0:p2)
+    return BSplineManifold(b,(P1,P3))
+end
+@inline function (M::AbstractBSplineManifold{3,p})(::Colon,::Colon,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j3 = intervalindex(P3,t3)
+    B3 = bsplinebasisall(P3,j3,t3)
+    b = sum(a[:,:,j3+i3]*B3[1+i3] for i3 in 0:p3)
+    return BSplineManifold(b,(P1,P2))
+end
+@inline function (M::AbstractBSplineManifold{3,p})(t1::Real,t2::Real,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j1 = intervalindex(P1,t1)
+    j2 = intervalindex(P2,t2)
+    B1 = bsplinebasisall(P1,j1,t1)
+    B2 = bsplinebasisall(P2,j2,t2)
+    b = sum(a[j1+i1,j2+i2,:]*B1[1+i1]*B2[1+i2] for i1 in 0:p1, i2 in 0:p2)
+    return BSplineManifold(b,(P3,))
+end
+@inline function (M::AbstractBSplineManifold{3,p})(t1::Real,::Colon,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j1 = intervalindex(P1,t1)
+    j3 = intervalindex(P3,t3)
+    B1 = bsplinebasisall(P1,j1,t1)
+    B3 = bsplinebasisall(P3,j3,t3)
+    b = sum(a[j1+i1,:,j3+i3]*B1[1+i1]*B3[1+i3] for i1 in 0:p1, i3 in 0:p3)
+    return BSplineManifold(b,(P2,))
+end
+@inline function (M::AbstractBSplineManifold{3,p})(::Colon,t2::Real,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    j2 = intervalindex(P2,t2)
+    j3 = intervalindex(P3,t3)
+    B2 = bsplinebasisall(P2,j2,t2)
+    B3 = bsplinebasisall(P3,j3,t3)
+    b = sum(a[:,j2+i2,j3+i3]*B2[1+i2]*B3[1+i3] for i2 in 0:p2, i3 in 0:p3)
+    return BSplineManifold(b,(P1,))
+end
+
 # TODO add mappings higher dimensionnal B-spline manifold with @generated macro
diff --git a/src/_RationalBSplineManifold.jl b/src/_RationalBSplineManifold.jl
index 74612f46e..7523c4869 100644
--- a/src/_RationalBSplineManifold.jl
+++ b/src/_RationalBSplineManifold.jl
@@ -53,6 +53,8 @@ struct RationalBSplineManifold{Dim,Deg,C,S<:Tuple,T} <: AbstractRationalBSplineM
     end
 end
 
+Base.:(==)(M1::AbstractRationalBSplineManifold, M2::AbstractRationalBSplineManifold) = (bsplinespaces(M1)==bsplinespaces(M2)) & (controlpoints(M1)==controlpoints(M2)) & (weights(M1)==weights(M2))
+
 controlpoints(M::RationalBSplineManifold) = M.controlpoints
 weights(M::RationalBSplineManifold) = M.weights
 bsplinespaces(M::RationalBSplineManifold) = M.bsplinespaces
@@ -122,3 +124,123 @@ end
         :(return v/u)
     )
 end
+
+
+## currying
+# 1dim
+@inline function (M::AbstractRationalBSplineManifold{1})(::Colon)
+    a = copy(controlpoints(M))
+    w = copy(weights(M))
+    Ps = bsplinespaces(M)
+    return RationalBSplineManifold(a,w,Ps)
+end
+
+# 2dim
+@inline function (M::AbstractRationalBSplineManifold{2})(::Colon,::Colon)
+    a = copy(controlpoints(M))
+    w = copy(weights(M))
+    Ps = bsplinespaces(M)
+    return RationalBSplineManifold(a,w,Ps)
+end
+@inline function (M::AbstractRationalBSplineManifold{2,p})(t1::Real,::Colon) where p
+    p1, p2 = p
+    P1, P2 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j1 = intervalindex(P1,t1)
+    B1 = bsplinebasisall(P1,j1,t1)
+    w′ = sum(w[j1+i1,:]*B1[1+i1] for i1 in 0:p1)
+    a′ = sum(a[j1+i1,:].*w[j1+i1,:].*B1[1+i1] for i1 in 0:p1) ./ w′
+    return RationalBSplineManifold(a′,w′,(P2,))
+end
+@inline function (M::AbstractRationalBSplineManifold{2,p})(::Colon,t2::Real) where p
+    p1, p2 = p
+    P1, P2 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j2 = intervalindex(P2,t2)
+    B2 = bsplinebasisall(P2,j2,t2)
+    w′ = sum(w[:,j2+i2]*B2[1+i2] for i2 in 0:p2)
+    a′ = sum(a[:,j2+i2].*w[:,j2+i2].*B2[1+i2] for i2 in 0:p2) ./ w′
+    return RationalBSplineManifold(a′,w′,(P1,))
+end
+
+# 3dim
+@inline function (M::AbstractRationalBSplineManifold{3})(::Colon,::Colon,::Colon)
+    a = copy(controlpoints(M))
+    w = copy(weights(M))
+    Ps = bsplinespaces(M)
+    return RationalBSplineManifold(a,w,Ps)
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(t1::Real,::Colon,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j1 = intervalindex(P1,t1)
+    B1 = bsplinebasisall(P1,j1,t1)
+    w′ = sum(w[j1+i1,:,:]*B1[1+i1] for i1 in 0:p1)
+    a′ = sum(a[j1+i1,:,:].*w[j1+i1,:,:].*B1[1+i1] for i1 in 0:p1) ./ w′
+    return RationalBSplineManifold(a′,w′,(P2,P3))
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(::Colon,t2::Real,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j2 = intervalindex(P2,t2)
+    B2 = bsplinebasisall(P2,j2,t2)
+    w′ = sum(w[:,j2+i2,:]*B2[1+i2] for i2 in 0:p2)
+    a′ = sum(a[:,j2+i2,:].*w[:,j2+i2,:].*B2[1+i2] for i2 in 0:p2) ./ w′
+    return RationalBSplineManifold(a′,w′,(P1,P3))
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(::Colon,::Colon,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j3 = intervalindex(P3,t3)
+    B3 = bsplinebasisall(P3,j3,t3)
+    w′ = sum(w[:,:,j3+i3]*B3[1+i3] for i3 in 0:p3)
+    a′ = sum(a[:,:,j3+i3].*w[:,:,j3+i3].*B3[1+i3] for i3 in 0:p3) ./ w′
+    return RationalBSplineManifold(a′,w′,(P1,P2))
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(t1::Real,t2::Real,::Colon) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j1 = intervalindex(P1,t1)
+    j2 = intervalindex(P2,t2)
+    B1 = bsplinebasisall(P1,j1,t1)
+    B2 = bsplinebasisall(P2,j2,t2)
+    w′ = sum(w[j1+i1,j2+i2,:]*B1[1+i1]*B2[1+i2] for i1 in 0:p1, i2 in 0:p2)
+    a′ = sum(a[j1+i1,j2+i2,:].*w[j1+i1,j2+i2,:].*B1[1+i1].*B2[1+i2] for i1 in 0:p1, i2 in 0:p2) ./ w′
+    return RationalBSplineManifold(a′,w′,(P3,))
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(t1::Real,::Colon,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j1 = intervalindex(P1,t1)
+    j3 = intervalindex(P3,t3)
+    B1 = bsplinebasisall(P1,j1,t1)
+    B3 = bsplinebasisall(P3,j3,t3)
+    w′ = sum(w[j1+i1,:,j3+i3]*B1[1+i1]*B3[1+i3] for i1 in 0:p1, i3 in 0:p3)
+    a′ = sum(a[j1+i1,:,j3+i3].*w[j1+i1,:,j3+i3].*B1[1+i1].*B3[1+i3] for i1 in 0:p1, i3 in 0:p3) ./ w′
+    return RationalBSplineManifold(a′,w′,(P2,))
+end
+@inline function (M::AbstractRationalBSplineManifold{3,p})(::Colon,t2::Real,t3::Real) where p
+    p1, p2, p3 = p
+    P1, P2, P3 = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    j2 = intervalindex(P2,t2)
+    j3 = intervalindex(P3,t3)
+    B2 = bsplinebasisall(P2,j2,t2)
+    B3 = bsplinebasisall(P3,j3,t3)
+    w′ = sum(w[:,j2+i2,j3+i3]*B2[1+i2]*B3[1+i3] for i2 in 0:p2, i3 in 0:p3)
+    a′ = sum(a[:,j2+i2,j3+i3].*w[:,j2+i2,j3+i3].*B2[1+i2].*B3[1+i3] for i2 in 0:p2, i3 in 0:p3) ./ w′
+    return RationalBSplineManifold(a′,w′,(P1,))
+end


Test patch:
diff --git a/test/test_BSplineManifold.jl b/test/test_BSplineManifold.jl
index f2a289d5a..bef77b7b6 100644
--- a/test/test_BSplineManifold.jl
+++ b/test/test_BSplineManifold.jl
@@ -44,6 +44,9 @@
             @test_throws DomainError M(-5)
 
             @test map.(M,[0.2,0.5,0.6]) == M.([0.2,0.5,0.6])
+            @test M(:) == M
+            @test Base.mightalias(controlpoints(M), controlpoints(M))
+            @test !Base.mightalias(controlpoints(M), controlpoints(M(:)))
         end
     end
 
@@ -71,10 +74,17 @@
             M′′ = refinement(M, p₊=p₊, k₊=k₊)
             ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
             for t in ts
-                @test M(t...) ≈ M′(t...)
-                @test M(t...) ≈ M′′(t...)
+                t1, t2 = t
+                @test M(t1,t2) ≈ M′(t1,t2)
+                @test M(t1,t2) ≈ M′′(t1,t2)
+                @test M(t1,t2) ≈ M(t1,:)(t2)
+                @test M(t1,t2) ≈ M(:,t2)(t1)
             end
             @test_throws DomainError M(-5,-8)
+
+            @test M(:,:) == M
+            @test Base.mightalias(controlpoints(M), controlpoints(M))
+            @test !Base.mightalias(controlpoints(M), controlpoints(M(:,:)))
         end
     end
 
@@ -98,9 +108,20 @@
             M′′ = refinement(M, p₊=p₊, k₊=k₊)
             ts = [[rand(), 1 + 2 * rand(), 0.5] for _ in 1:10]
             for t in ts
-                @test M(t...) ≈ M′′(t...)
+                t1, t2, t3 = t
+                @test M(t1,t2,t3) ≈ M′′(t1,t2,t3)
+                @test M(t1,t2,t3) ≈ M(t1,:,:)(t2,t3)
+                @test M(t1,t2,t3) ≈ M(:,t2,:)(t1,t3)
+                @test M(t1,t2,t3) ≈ M(:,:,t3)(t1,t2)
+                @test M(t1,t2,t3) ≈ M(t1,t2,:)(t3)
+                @test M(t1,t2,t3) ≈ M(t1,:,t3)(t2)
+                @test M(t1,t2,t3) ≈ M(:,t2,t3)(t1)
             end
             @test_throws DomainError M(-5,-8,-120)
+
+            @test M(:,:,:) == M
+            @test Base.mightalias(controlpoints(M), controlpoints(M))
+            @test !Base.mightalias(controlpoints(M), controlpoints(M(:,:,:)))
         end
     end
 end
diff --git a/test/test_RationalBSplineManifold.jl b/test/test_RationalBSplineManifold.jl
index 46aaa650c..119d2f6d4 100644
--- a/test/test_RationalBSplineManifold.jl
+++ b/test/test_RationalBSplineManifold.jl
@@ -7,7 +7,7 @@
         p = 2
         P = BSplineSpace{p}(k)
 
-        M = BasicBSpline.RationalBSplineManifold(a,w,(P,))
+        M = RationalBSplineManifold(a,w,(P,))
         for _ in 1:100
             t1 = 2rand()
             @test norm(M(t1)) ≈ 1
@@ -16,75 +16,119 @@
         @test_throws DomainError M(2.1)
     end
 
-    @testset "compatibility with BSplineManifold (1dim)" begin
-        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
-        p1 = 4
-        P1 = BSplineSpace{p1}(k1)
-        n1 = dim(P1)
-        a = randn(n1)
-        w1 = ones(n1)
-        w2 = 2ones(n1)
-
-        M = BasicBSpline.BSplineManifold(a,(P1,))
-        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,))
-        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,))
-        for _ in 1:100
-            t1 = rand()
-            @test M(t1) ≈ M1(t1) atol=1e-14
-            @test M(t1) ≈ M2(t1) atol=1e-14
-        end
-    end
-
-    @testset "compatibility with BSplineManifold (2dim)" begin
-        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
-        k2 = KnotVector(rand(30)) + 5KnotVector(0,1)
-        p1 = 2
-        p2 = 3
-        P1 = BSplineSpace{p1}(k1)
-        P2 = BSplineSpace{p2}(k2)
-        n1 = dim(P1)
-        n2 = dim(P2)
-        a = randn(n1,n2)
-        w1 = ones(n1,n2)
-        w2 = 2ones(n1,n2)
-
-        M = BasicBSpline.BSplineManifold(a,(P1,P2))
-        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,P2))
-        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,P2))
-        for _ in 1:100
-            t1 = rand()
-            t2 = rand()
-            @test M(t1,t2) ≈ M1(t1,t2) atol=1e-14
-            @test M(t1,t2) ≈ M2(t1,t2) atol=1e-14
-        end
-    end
-
-    @testset "compatibility with BSplineManifold (3dim)" begin
-        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
-        k2 = KnotVector(rand(30)) + 5KnotVector(0,1)
-        k3 = KnotVector(rand(30)) + 5KnotVector(0,1)
+    @testset "general cases" begin
         p1 = 2
         p2 = 3
         p3 = 4
+        k1 = KnotVector(rand(30)) + (p1+1)*KnotVector(0,1)
+        k2 = KnotVector(rand(30)) + (p2+1)*KnotVector(0,1)
+        k3 = KnotVector(rand(30)) + (p3+1)*KnotVector(0,1)
         P1 = BSplineSpace{p1}(k1)
         P2 = BSplineSpace{p2}(k2)
         P3 = BSplineSpace{p3}(k3)
         n1 = dim(P1)
         n2 = dim(P2)
         n3 = dim(P3)
-        a = randn(n1,n2,n3)
-        w1 = ones(n1,n2,n3)
-        w2 = 2ones(n1,n2,n3)
 
-        M = BasicBSpline.BSplineManifold(a,(P1,P2,P3))
-        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,P2,P3))
-        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,P2,P3))
-        for _ in 1:100
-            t1 = rand()
-            t2 = rand()
-            t3 = rand()
-            @test M(t1,t2,t3) ≈ M1(t1,t2,t3) atol=1e-14
-            @test M(t1,t2,t3) ≈ M2(t1,t2,t3) atol=1e-14
+        @testset "RationalBSpilneManifold definition" begin
+            @testset "1dim" begin
+                a = randn(ComplexF64,n1)
+                w = ones(n1)+rand(n1)
+                R = RationalBSplineManifold(a,w,(P1,))
+                A = BSplineManifold(a.*w,(P1,))
+                W = BSplineManifold(w,(P1,))
+                @test R(:) == R
+                @test Base.mightalias(controlpoints(R), controlpoints(R))
+                @test !Base.mightalias(controlpoints(R), controlpoints(R(:)))
+                for _ in 1:10
+                    t1 = rand()
+                    @test R(t1) ≈ A(t1)/W(t1)
+                end
+            end
+            @testset "2dim" begin
+                a = randn(ComplexF64,n1,n2)
+                w = ones(n1,n2)+rand(n1,n2)
+                R = RationalBSplineManifold(a,w,(P1,P2))
+                A = BSplineManifold(a.*w,(P1,P2))
+                W = BSplineManifold(w,(P1,P2))
+                @test R(:,:) == R
+                @test Base.mightalias(controlpoints(R), controlpoints(R))
+                @test !Base.mightalias(controlpoints(R), controlpoints(R(:,:)))
+                for _ in 1:10
+                    t1 = rand()
+                    t2 = rand()
+                    @test R(t1,t2) ≈ A(t1,t2)/W(t1,t2)
+                    @test R(t1,t2) ≈ R(t1,:)(t2)
+                    @test R(t1,t2) ≈ R(:,t2)(t1)
+                end
+            end
+            @testset "3dim" begin
+                a = randn(ComplexF64,n1,n2,n3)
+                w = ones(n1,n2,n3)+rand(n1,n2,n3)
+                R = RationalBSplineManifold(a,w,(P1,P2,P3))
+                A = BSplineManifold(a.*w,(P1,P2,P3))
+                W = BSplineManifold(w,(P1,P2,P3))
+                @test R(:,:,:) == R
+                @test Base.mightalias(controlpoints(R), controlpoints(R))
+                @test !Base.mightalias(controlpoints(R), controlpoints(R(:,:,:)))
+                for _ in 1:10
+                    t1 = rand()
+                    t2 = rand()
+                    t3 = rand()
+                    @test R(t1,t2,t3) ≈ A(t1,t2,t3)/W(t1,t2,t3)
+                    @test R(t1,t2,t3) ≈ R(t1,:,:)(t2,t3)
+                    @test R(t1,t2,t3) ≈ R(:,t2,:)(t1,t3)
+                    @test R(t1,t2,t3) ≈ R(:,:,t3)(t1,t2)
+                    @test R(t1,t2,t3) ≈ R(t1,t2,:)(t3)
+                    @test R(t1,t2,t3) ≈ R(t1,:,t3)(t2)
+                    @test R(t1,t2,t3) ≈ R(:,t2,t3)(t1)
+                end
+            end
+        end
+
+        @testset "compatibility with BSplineManifold" begin
+            @testset "1dim" begin
+                a = randn(n1)
+                w1 = ones(n1)
+                w2 = 2ones(n1)
+                M = BSplineManifold(a,(P1,))
+                M1 = RationalBSplineManifold(a,w1,(P1,))
+                M2 = RationalBSplineManifold(a,w2,(P1,))
+                for _ in 1:100
+                    t1 = rand()
+                    @test M(t1) ≈ M1(t1) atol=1e-14
+                    @test M(t1) ≈ M2(t1) atol=1e-14
+                end
+            end
+            @testset "2dim" begin
+                a = randn(n1,n2)
+                w1 = ones(n1,n2)
+                w2 = 2ones(n1,n2)
+                M = BSplineManifold(a,(P1,P2))
+                M1 = RationalBSplineManifold(a,w1,(P1,P2))
+                M2 = RationalBSplineManifold(a,w2,(P1,P2))
+                for _ in 1:100
+                    t1 = rand()
+                    t2 = rand()
+                    @test M(t1,t2) ≈ M1(t1,t2) atol=1e-14
+                    @test M(t1,t2) ≈ M2(t1,t2) atol=1e-14
+                end
+            end
+            @testset "3dim" begin
+                a = randn(n1,n2,n3)
+                w1 = ones(n1,n2,n3)
+                w2 = 2ones(n1,n2,n3)
+                M = BSplineManifold(a,(P1,P2,P3))
+                M1 = RationalBSplineManifold(a,w1,(P1,P2,P3))
+                M2 = RationalBSplineManifold(a,w2,(P1,P2,P3))
+                for _ in 1:100
+                    t1 = rand()
+                    t2 = rand()
+                    t3 = rand()
+                    @test M(t1,t2,t3) ≈ M1(t1,t2,t3) atol=1e-14
+                    @test M(t1,t2,t3) ≈ M2(t1,t2,t3) atol=1e-14
+                end
+            end
         end
     end
 end