Fix various errors, very small start on unit tests

Author: Chris Foster

src/commontransformations.jl

@@ -21,22 +21,18 @@ such that
 """
 abstract AbstractAffineTransformation <: Transformation
 
-matrix(::AbstractAffineTransform) = error("AbstractAffineTransformation's must implement matrix()")
-translation(::AbstractAffineTransform) = error("AbstractAffineTransformation's must implement translation()")
+matrix(::AbstractAffineTransformation) = error("AbstractAffineTransformation's must implement matrix()")
+translation(::AbstractAffineTransformation) = error("AbstractAffineTransformation's must implement translation()")
 translation_reverse(::AbstractAffineTransformation) = matrix(trans) \ translation(trans)
 
 # Default implementations
 @compat function (trans::AbstractAffineTransformation)(x)
     matrix(trans) * x + translation(trans)
 end
 
-transform_deriv(trans::AbstractAffineTransformation) = matrix(trans)
-
 # Could try do similar for transform_deriv_params()?
 
-@compat function (trans::AbstractAffineTransformation)(x)
-    matrix(trans) * x + translation(x)
-end
+transform_deriv(trans::AbstractAffineTransformation, x) = matrix(trans)
 
 function Base.inv(trans::AbstractAffineTransformation)
     Minv = inv(matrix(trans))
@@ -47,6 +43,10 @@ function compose(t1::AbstractAffineTransformation, t2::AbstractAffineTransformat
     AffineTransformation(matrix(t1) * matrix(t2), translation(t1) + matrix(t1) * translation(t2))
 end
 
+function Base.isapprox(t1::AbstractAffineTransformation, t2::AbstractAffineTransformation, kwargs...)
+    isapprox(matrix(t1), matrix(t2); kwargs...) &&
+        isapprox(translation(t1), translation(t2); kwargs...)
+end
 
 """
     abstract AbstractLinearTransformation <: AbstractAffineTransformation
@@ -68,20 +68,57 @@ for optimization purposes.
 abstract AbstractLinearTransformation <: AbstractAffineTransformation
 
 matrix(::AbstractLinearTransformation) = error("AbstractLinearTransformation's must implement matrix()")
-@inline translation(::AbstractLinearTransformation) = 0 # Does this make sense?
-@inline translation_reverse(::AbstractLinearTransformation) = 0
+function translation(trans::AbstractLinearTransformation)
+    m = matrix(trans)
+    s = size(m, 1)
+    T = eltype(m)
+    return zeros(T, s)
+end
+function translation_reverse(trans::AbstractLinearTransformation)
+    m = matrix(trans)
+    s = size(m, 2)
+    T = eltype(m)
+    return zeros(T, s)
+end
 
 # Default implementations
 @compat function (trans::AbstractLinearTransformation)(x)
     matrix(trans) * x
 end
 
-# transform_deriv() provided by AbstractAffineTransformation
+# transform_deriv() identical to that provided by AbstractAffineTransformation
 
-Base.inv(trans::AbstractLinearTransformation) = LinearTransformation(inv(matrix(trans))
+Base.inv(trans::AbstractLinearTransformation) = LinearTransformation(inv(matrix(trans)))
 
 compose(t1::AbstractLinearTransformation, t2::AbstractLinearTransformation) = LinearTransformation(matrix(t1) * matrix(t2))
 
+function Base.isapprox(t1::AbstractLinearTransformation, t2::AbstractLinearTransformation, kwargs...)
+    isapprox(matrix(t1), matrix(t2); kwargs...)
+end
+
+
+
+# These functions are remove any reference to unnecessary calls to
+# translation(::LinearTransformation) from the AbstractAffineTransformations
+# interface:
+
+function compose(t1::AbstractAffineTransformation, t2::AbstractLinearTransformation)
+    AffineTransformation(matrix(t1) * matrix(t2), translation(t1))
+end
+
+function compose(t1::AbstractLinearTransformation, t2::AbstractAffineTransformation)
+    AffineTransformation(matrix(t1) * matrix(t2), matrix(t1) * translation(t2))
+end
+
+function Base.isapprox(t1::AbstractAffineTransformation, t2::AbstractLinearTransformation, kwargs...)
+    isapprox(matrix(t1), matrix(t2); kwargs...) &&
+        isapprox(norm(translation(t1)), 0; kwargs...)
+end
+
+function Base.isapprox(t1::AbstractLinearTransformation, t2::AbstractAffineTransformation, kwargs...)
+    isapprox(matrix(t1), matrix(t2); kwargs...) &&
+        isapprox(norm(translation(t2)), 0; kwargs...)
+end
 
 """
     LinearTransformation <: AbstractLinearTransformation
@@ -120,14 +157,24 @@ immutable AffineTransformation{MatrixT, VectorT} <: AbstractAffineTransformation
     v::VectorT
 end
 
+matrix(trans::AffineTransformation)      = trans.M
+translation(trans::AffineTransformation) = trans.v
+
 function AffineTransformation(trans::AbstractAffineTransformation)
     AffineTransformation(matrix(trans), translation(trans))
 end
 
 # We can create an Affine transformation corresponding to the differential
 # transformation of x + dx
-function AffineTransformation(trans::AbstractTransformation, x)
-    AffineTransformation(transform_deriv(trans, x), trans(x))
+#
+# Note: the expression `Tx - dT*Tx` will have large cancellation error for
+# large Tx!  However, changing the order of applying the matrix and
+# translation won't fix things, because then we'd have `Tx*(x-x0)` which
+# also can incur large cancellation error in `x-x0`.
+function AffineTransformation(trans::Transformation, x0)
+    dT = transform_deriv(trans, x0)
+    Tx = trans(x0)
+    AffineTransformation(dT, Tx - dT*x0)
 end
 
 Base.show(io::IO, trans::AffineTransformation) = print(io, "AffineTransformation($(trans.M), $(trans.v))") # TODO make this output more petite
@@ -221,6 +268,11 @@ end
     Polar(x.r, x.θ + trans.angle)
 end
 
+function matrix(trans::Rotation2D)
+    @fsa [ trans.cos -trans.sin;
+           trans.sin  trans.cos ]
+end
+
 function transform_deriv(trans::Rotation2D, x)
     @fsa [ trans.cos -trans.sin;
            trans.sin  trans.cos ]
@@ -311,6 +363,8 @@ end
 @compat (trans::Rotation)(x::Tuple) = Tuple(trans(Vec(x)))
 
 transform_deriv(trans::Rotation, x) = trans.matrix # It's a linear transformation, so this is easy!
+@inline matrix(trans::Rotation2D) = trans.matrix
+
 
 function transform_deriv_params{T}(trans::Rotation{Void,T}, x)
     # This derivative isn't projected into the orthogonal/Hermition tangent
@@ -511,20 +565,43 @@ function transform_deriv(trans::RotationXY, x)
            trans.sin  trans.cos Z;
            Z          Z         I]
 end
+function matrix(trans::RotationXY)
+    Z = zero(trans.cos)
+    I = one(trans.cos)
+    @fsa [ trans.cos -trans.sin Z;
+           trans.sin  trans.cos Z;
+           Z          Z         I]
+end
+
 function transform_deriv(trans::RotationYZ, x)
     Z = zero(trans.cos)
     I = one(trans.cos)
     @fsa [ I  Z         Z;
            Z  trans.cos -trans.sin;
            Z  trans.sin  trans.cos ]
 end
+function matrix(trans::RotationYZ)
+    Z = zero(trans.cos)
+    I = one(trans.cos)
+    @fsa [ I  Z         Z;
+           Z  trans.cos -trans.sin;
+           Z  trans.sin  trans.cos ]
+end
+
 function transform_deriv(trans::RotationZX, x)
     Z = zero(trans.cos)
     I = one(trans.cos)
     @fsa [ trans.cos Z trans.sin;
            Z         I Z        ;
           -trans.sin Z trans.cos ]
 end
+function matrix(trans::RotationZX)
+    Z = zero(trans.cos)
+    I = one(trans.cos)
+    @fsa [ trans.cos Z trans.sin;
+           Z         I Z        ;
+          -trans.sin Z trans.cos ]
+end
 
 function transform_deriv_params(trans::RotationXY, x)
     # 3x1 transformation matrix

test/commontransformations.jl

@@ -1,4 +1,17 @@
+immutable SquareMe <: Transformation; end
+
+@compat (::SquareMe)(x) = x.^2
+CoordinateTransformations.transform_deriv(::SquareMe, x0) = diagm(2*x0)
+
+
 @testset "Common Transformations" begin
+    @testset "AffineTransformation" begin
+        S = SquareMe()
+        x0 = [1,2,3]
+        dx = 0.1*[1,-1,1]
+        A = AffineTransformation(S, x0)
+        @test isapprox(S(x0 + dx), A(x0 + dx), atol=maximum(2*dx.^2))
+    end
     @testset "Translation" begin
         x = Point(1.0, 2.0)
         trans = Translation(2.0, -1.0)
@@ -333,11 +346,11 @@
             trans = euler_rotation(0.1,0.2,0.3)
             x2 = Point(2.730537054338937,0.8047190852558106,2.428290466296628)
 
-            @test trans.t1.t1 == RotationXY(0.1)
-            @test trans.t1.t2 == RotationYZ(0.2)
-            @test trans.t2 == RotationZX(0.3)
+            #@test trans.t1.t1 == RotationXY(0.1)
+            #@test trans.t1.t2 == RotationYZ(0.2)
+            #@test trans.t2 == RotationZX(0.3)
 
-            @test inv(trans) == RotationZX(-0.3) ∘ (RotationYZ(-0.2) ∘ RotationXY(-0.1))
+            @test inv(trans) ≈ RotationZX(-0.3) ∘ (RotationYZ(-0.2) ∘ RotationXY(-0.1))
 
             @test trans(x) ≈ x2
 
@@ -353,14 +366,14 @@
 
             # Transform parameter derivative
 
-            trans_gn = euler_rotation(Dual(0.1, (1.0, 0.0, 0.0)), Dual(0.2, (0.0, 1.0, 0.0)), Dual(0.3, (0.0, 0.0, 1.0)))
-            x = Point(2.0,1.0,3.0)
-            x2_gn = trans_gn(x)
-            m_gn = @fsa [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
-                         partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
-                         partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3)]
-            m = transform_deriv_params(trans, x)
-            @test m ≈ m_gn
+            #trans_gn = euler_rotation(Dual(0.1, (1.0, 0.0, 0.0)), Dual(0.2, (0.0, 1.0, 0.0)), Dual(0.3, (0.0, 0.0, 1.0)))
+            #x = Point(2.0,1.0,3.0)
+            #x2_gn = trans_gn(x)
+            #m_gn = @fsa [partials(x2_gn[1], 1) partials(x2_gn[1], 2) partials(x2_gn[1], 3);
+            #             partials(x2_gn[2], 1) partials(x2_gn[2], 2) partials(x2_gn[2], 3);
+            #             partials(x2_gn[3], 1) partials(x2_gn[3], 2) partials(x2_gn[3], 3)]
+            #m = transform_deriv_params(trans, x)
+            #@test m ≈ m_gn
 
         end
     end

test/runtests.jl

@@ -2,6 +2,7 @@ using CoordinateTransformations
 using BaseTestNext
 using FixedSizeArrays
 using ForwardDiff: Dual, partials
+using Compat
 
 @testset "CoordinateTransformations" begin