diff --git a/Project.toml b/Project.toml index 7d1c4c9b..1af2f265 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,16 @@ name = "Quaternions" uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" -version = "0.5.7" +version = "0.6.0-DEV" [deps] -DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -DualNumbers = "0.5, 0.6" julia = "1" [extras] -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ForwardDiff", "Test"] +test = ["Test"] diff --git a/README.md b/README.md index b631a8b4..542578bf 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # Quaternions.jl -A Julia module with quaternion, octonion and dual-quaternion functionality +A Julia implementation of quaternions. [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaGeometry.github.io/Quaternions.jl/stable) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaGeometry.github.io/Quaternions.jl/dev) @@ -62,5 +62,5 @@ Implemented functions are: rand randn -Currently, this package supports `DualQuaternion` and `Octonion` types, but these will be removed in the next breaking release. -See https://github.com/JuliaGeometry/Quaternions.jl/issues/90 and https://github.com/JuliaGeometry/Quaternions.jl/pull/92 for more information. +Currently, this package supports the `Octonion` type, but this will be removed in the next breaking release. +See https://github.com/JuliaGeometry/Quaternions.jl/issues/90 for more information. diff --git a/docs/Project.toml b/docs/Project.toml index dfa65cd1..b700a442 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,2 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" + +[compat] +ForwardDiff = "0.10" diff --git a/docs/make.jl b/docs/make.jl index 892b1486..ddfde1ef 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,6 +14,7 @@ makedocs(; ), pages=[ "Home" => "index.md", + "Examples" => ["examples/dual_quaternions.md"], ], ) diff --git a/docs/src/examples/dual_quaternions.md b/docs/src/examples/dual_quaternions.md new file mode 100644 index 00000000..1759d7dc --- /dev/null +++ b/docs/src/examples/dual_quaternions.md @@ -0,0 +1,175 @@ +# Dual quaternions + +## Introduction + +The [dual quaternions](https://en.wikipedia.org/wiki/Dual_quaternion) are an example of "biquaternions." +They can be represented equivalently either as a [dual number](https://en.wikipedia.org/wiki/Dual_number) where both both the "primal" and "tangent" part are quaternions + +```math +d = q_0 + q_e \epsilon = (s_0 + a_0 i + b_0 j + c_0 k) + (s_e + a_e i + b_e j + c_e k) \epsilon +``` + +or as a quaternion where the scalar part and three imaginary parts are all dual numbers + +```math +d = s + ai + bj + ck = (s_0 + s_e \epsilon) + (a_0 + a_e \epsilon) i + (b_0 + b_e \epsilon) j + (c_0 + c_e \epsilon) k. +``` + +Like unit quaternions can compactly representation rotations in 3D space, dual quaternions can compactly represent rigid transformations (rotation with translation). + +Without any special glue code, we can construct a dual quaternion by composing `ForwardDiff.Dual` and [`Quaternion`](@ref); this uses the second representation described above: + +!!! note + Previously this package contained a specialized `DualQuaternion` type. + This was removed in v0.6.0 because it offered nothing extra over composing [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) and Quaternions. + +## Utility functions + +First let's load the packages: + +```@example dualquat +using Quaternions, ForwardDiff, Random +``` + +Then we'll create some utility types/functions: + +```@example dualquat +const DualQuaternion{T} = Quaternion{ForwardDiff.Dual{Nothing,T,1}} + +purequat(p::AbstractVector) = quat(false, @views(p[begin:begin+2])...) + +dual(x::Real, v::Real) = ForwardDiff.Dual(x, v) + +function dualquat(_q0::Union{Real,Quaternion}, _qe::Union{Real,Quaternion}) + q0 = quat(_q0) + qe = quat(_qe) + Quaternion( + dual(real(q0), real(qe)), + dual.(imag_part(q0), imag_part(qe))..., + ) +end + +function primal(d::DualQuaternion) + return Quaternion( + ForwardDiff.value(real(d)), + ForwardDiff.value.(imag_part(d))..., + ) +end + +function tangent(d::DualQuaternion) + return Quaternion( + ForwardDiff.partials(real(d), 1), + ForwardDiff.partials.(imag_part(d), 1)..., + ) +end + +function dualconj(d::DualQuaternion) + de = tangent(d) + return dualquat(conj(primal(d)), quat(-real(de), imag_part(de)...)) +end + +rotation_part(d::DualQuaternion) = primal(d) + +translation_part(d::DualQuaternion) = dualquat(true, conj(rotation_part(d)) * tangent(d)) + +# first=true returns the translation performed before the rotation: R(p+t) +# first=false returns the translation performed after the rotation: R(p)+t +function translation(d::DualQuaternion; first::Bool=true) + v = first ? primal(d)' * tangent(d) : tangent(d) * primal(d)' + return collect(2 .* imag_part(v)) +end + +function transform(d::DualQuaternion, p::AbstractVector) + dp = dualquat(true, purequat(p)) + dpnew = d * dp * dualconj(d) + pnew_parts = imag_part(tangent(dpnew)) + pnew = similar(p, eltype(pnew_parts)) + pnew .= pnew_parts + return pnew +end + +function transformationmatrix(d::DualQuaternion) + R = rotationmatrix(rotation_part(d)) + t = translation(d; first=false) + T = similar(R, 4, 4) + T[1:3, 1:3] .= R + T[1:3, 4] .= t + T[4, 1:3] .= 0 + T[4, 4] = 1 + return T +end + +randdualquat(rng::AbstractRNG,T=Float64) = dualquat(rand(rng, Quaternion{T}), rand(rng, Quaternion{T})) +randdualquat(T=Float64) = randdualquat(Random.GLOBAL_RNG,T) +nothing # hide +``` + +## Example: transforming a point + +Now we'll create a unit dual quaternion. +```@repl dualquat +x = sign(randdualquat()) +``` + +`sign(q) == q / abs(q)` both normalizes the primal part of the dual quaternion and makes the tangent part perpendicular to it. + +```@repl dualquat +abs(primal(x)) ≈ 1 +isapprox(real(primal(x)' * tangent(x)), 0; atol=1e-10) +``` + +Here's how we use dual quaternions to transform a point: + +```@repl dualquat +p = randn(3) +``` + +```@repl dualquat +transform(x, p) +``` + +## Example: homomorphism from unit dual quaternions to the transformation matrices + +Each unit dual quaternion can be mapped to an affine transformation matrix ``T``. +``T`` can be used to transform a vector ``p`` like this: + +```math +T \begin{pmatrix} p \\ 1\end{pmatrix} = \begin{pmatrix} R & t \\ 0^\mathrm{T} & 1\end{pmatrix} \begin{pmatrix} p \\ 1\end{pmatrix} = \begin{pmatrix} Rp + t \\ 1\end{pmatrix}, +``` +where ``R`` is a rotation matrix, and ``t`` is a translation vector. +Our helper function `transformationmatrix` maps from a unit dual quaternion to such an affine matrix. + +```@repl dualquat +y = sign(randdualquat()) +``` + +```@repl dualquat +X = transformationmatrix(x) +Y = transformationmatrix(y) +XY = transformationmatrix(x*y) +X*Y ≈ XY +``` + +We can check that our transformation using the unit dual quaternion gives the same result as transforming with an affine transformation matrix: + +```@repl dualquat +transform(x, p) ≈ (X * vcat(p, 1))[1:3] +``` + +## Example: motion planning + +For unit quaternions, spherical linear interpolation with [`slerp`](@ref) can be used to interpolate between two rotations with unit quaternions, which can be used to plan motion between two orientations. +Similarly, we can interpolate between unit dual quaternions to plan motion between two rigid poses. +Conveniently, we can do this using the exact same `slerp` implementation. + +```@repl dualquat +slerp(x, y, 0) ≈ x +``` + +```@repl dualquat +slerp(x, y, 1) ≈ y +``` + +```@repl dualquat +slerp(x, y, 0.3) +``` diff --git a/docs/src/index.md b/docs/src/index.md index 0cd4dda2..cd036a85 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,6 @@ # Quaternions.jl -A Julia package implementing [quaternions](https://en.wikipedia.org/wiki/Quaternion), [octonions](https://en.wikipedia.org/wiki/Octonion) and [dual-quaternions](https://en.wikipedia.org/wiki/Dual_quaternion) +A Julia package implementing [quaternions](https://en.wikipedia.org/wiki/Quaternion) and [octonions](https://en.wikipedia.org/wiki/Octonion). !!! note "Documentation" The documentation is still work in progress. diff --git a/src/DualQuaternion.jl b/src/DualQuaternion.jl deleted file mode 100644 index 655dd63d..00000000 --- a/src/DualQuaternion.jl +++ /dev/null @@ -1,185 +0,0 @@ -using DualNumbers - -struct DualQuaternion{T<:Real} <: Number - q0::Quaternion{T} - qe::Quaternion{T} - norm::Bool - function DualQuaternion{T}(q0, qe, norm) where T <: Real - Base.depwarn("`DualQuaternion` is deprecated and will be removed in the next breaking release. Use `Quaternion{ForwardDiff.Dual}` instead.", :DualQuaternion) - return new{T}(q0, qe, norm) - end -end - -function DualQuaternion(q0::Quaternion{T}, qe::Quaternion{T}, norm::Bool) where T <: Real - return DualQuaternion{T}(q0, qe, norm) -end - -DualQuaternion{T}(dq::DualQuaternion) where {T<:Real} = DualQuaternion{T}(dq.q0, dq.qe, dq.norm) -function DualQuaternion{T}(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool=false) where {T<:Real} - return DualQuaternion{T}( - Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4), n), - Quaternion(DualNumbers.epsilon(d1), DualNumbers.epsilon(d2), DualNumbers.epsilon(d3), DualNumbers.epsilon(d4)), - n, - ) -end -function DualQuaternion{T}(q0::Quaternion) where {T<:Real} - return DualQuaternion{T}(convert(Quaternion{T}, q0), zero(Quaternion{T}), q0.norm) -end -function DualQuaternion{T}(d::Dual) where {T<:Real} - return DualQuaternion( - Quaternion{T}(DualNumbers.value(d)), - Quaternion{T}(DualNumbers.epsilon(d)), - (DualNumbers.value(d)==one(DualNumbers.value(d))) & iszero(DualNumbers.epsilon(d))) -end -function DualQuaternion{T}(x::Real) where {T<:Real} - return DualQuaternion(convert(Quaternion{T}, x), zero(Quaternion{T}), abs(x) == one(x)) -end - -DualQuaternion(q0::Quaternion, qe::Quaternion) = DualQuaternion(promote(q0, qe)..., false) -DualQuaternion(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool=false) = - DualQuaternion(Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4), n), - Quaternion(DualNumbers.epsilon(d1), DualNumbers.epsilon(d2), DualNumbers.epsilon(d3), DualNumbers.epsilon(d4)), n) -DualQuaternion(x::Real) = DualQuaternion(Quaternion(x), Quaternion(zero(x)), abs(x) == one(x)) -DualQuaternion(d::Dual) = DualQuaternion(Quaternion(DualNumbers.value(d)), Quaternion(DualNumbers.epsilon(d)), (DualNumbers.value(d)==one(DualNumbers.value(d))) & iszero(DualNumbers.epsilon(d))) -DualQuaternion(q::Quaternion) = DualQuaternion(q, zero(q), q.norm) -DualQuaternion(a::Vector) = DualQuaternion(zero(Quaternion{typeof(a[1])}), Quaternion(a)) - -const DualQuaternionF16 = DualQuaternion{Float16} -const DualQuaternionF32 = DualQuaternion{Float32} -const DualQuaternionF64 = DualQuaternion{Float64} - -promote_rule(::Type{DualQuaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)} -promote_rule(::Type{DualQuaternion{T}}, ::Type{Dual{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)} -promote_rule(::Type{DualQuaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)} -promote_rule(::Type{DualQuaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)} - -dualquat(q1, q2) = DualQuaternion(q1, q2) -dualquat(q1, q2, n) = DualQuaternion(q1, q2, n) -dualquat(d1, d2, d3, d4) = DualQuaternion(d1, d2, d3, d4) -dualquat(d1, d2, d3, d4, n) = DualQuaternion(d1, d2, d3, d4, n) -dualquat(x) = DualQuaternion(x) - -Q0(dq::DualQuaternion) = dq.q0 -Qe(dq::DualQuaternion) = dq.qe - -(/)(dq::DualQuaternion, x::Real) = DualQuaternion(dq.q0 / x, dq.qe / x) - -(/)(dq::DualQuaternion, d::Dual) = - DualQuaternion(dual(dq.q0.s, dq.qe.s) / d, - dual(dq.q0.v1, dq.qe.v1) / d, - dual(dq.q0.v2, dq.qe.v2) / d, - dual(dq.q0.v3, dq.qe.v3) / d) - -abs2(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) : - dual(abs2(dq.q0), - 2.0 * (dq.q0.s * dq.qe.s + - dq.q0.v1 * dq.qe.v1 + - dq.q0.v2 * dq.qe.v2 + - dq.q0.v3 * dq.qe.v3)) - -abs(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) : sqrt(abs2(dq)) -float(dq::DualQuaternion{T}) where T = convert(DualQuaternion{float(T)}, dq) - -conj(dq::DualQuaternion) = DualQuaternion(conj(dq.q0), conj(dq.qe), dq.norm) -dconj(dq::DualQuaternion) = DualQuaternion(dq.q0, -dq.qe, dq.norm) - -inv(dq::DualQuaternion) = dq.norm ? conj(dq) : conj(dq) / abs2(dq) - -function normalize(dq::DualQuaternion) - if (dq.norm) - return dq - end - a = abs(dq) - if abs(a) > 0 - qa = dq / a - dualquat(qa.q0, qa.qe, true) - else - dq - end -end - -function normalizea(dq::DualQuaternion{T}) where {T} - if (dq.norm) - return (dq, one(Dual{T})) - end - a = abs(dq) - if abs(a) > 0 - qa = dq / a - dualquat(qa.q0, qa.qe, true), a - else - dq, zero(Dual{T}) - end -end - -(-)(dq::DualQuaternion) = DualQuaternion(-dq.q0, -dq.qe, dq.norm) - -(+)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 + dw.q0, dq.qe + dw.qe) -(-)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 - dw.q0, dq.qe - dw.qe) -(*)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 * dw.q0, - dq.q0 * dw.qe + dq.qe * dw.q0, - dq.norm && dw.norm) -(*)(dq::DualQuaternion, d::Dual) = (*)(Base.promote(dq, d)...) -(*)(d::Dual, dq::DualQuaternion) = (*)(Base.promote(d, dq)...) -(/)(dq::DualQuaternion, dw::DualQuaternion) = dq * inv(dw) -(==)(q::DualQuaternion, w::DualQuaternion) = (q.q0 == w.q0) & (q.qe == w.qe) # ignore .norm field - -function angleaxis(dq::DualQuaternion) - tq = dq.qe * conj(dq.q0) - t = [2.0 * tq.v1, 2.0 * tq.v2, 2.0 * tq.v3] - q0s = dq.q0.s - th0, s0 = angleaxis(dq.q0) - sq0 = quat(0.0, s0) - if abs(abs(q0s) - one(q0s)) == 0 - th = dual(th0, 0.5 * abs(quat(0, t))) - th, dualquat(sq0) - else - th = dual(th0, 0.5 * dot(t, s0)) - s0c1 = cross(s0, t) - tanth = tan(th0) - s0c2 = (s0c1 / tanth + t) * 0.5 - sqev = cross(s0c2, s0) - th, dualquat(sq0, quat(0.0, sqev)) - end -end - -function angle(dq::DualQuaternion) - th, ax = angleaxis(dq) - th -end - -function axis(dq::DualQuaternion) - th, ax = angleaxis(dq) - ax -end - -function exp(dq::DualQuaternion) - se = dual(dq.q0.s, dq.qe.s) - se = exp(se) - dq = dualquat(quat(0.0, imag_part(dq.q0)...), quat(0.0, imag_part(dq.qe)...)) - dq, th = normalizea(dq) - if dq.norm - dualquat(se) * (dualquat(cos(th)) + dq * dualquat(sin(th))) - else - dualquat(se) - end -end - -function log(dq::DualQuaternion) - dq, a = normalizea(dq) - sl = log(a) - th, s = angleaxis(dq) - s * dualquat(th) + dualquat(sl) -end - -(^)(dq::DualQuaternion, dw::DualQuaternion) = exp(dw * log(dq)) - -function sqrt(dq::DualQuaternion) - exp(0.5 * log(dq)) -end - -dualquatrand() = dualquat(quatrand(), quatrand()) -ndualquatrand() = normalize(dualquatrand()) - -function rand(rng::AbstractRNG, ::Random.SamplerType{DualQuaternion{T}}) where {T<:Real} - return DualQuaternion{T}(rand(rng, Quaternion{T}), rand(rng, Quaternion{T}), false) -end diff --git a/src/Quaternions.jl b/src/Quaternions.jl index c8072055..9b47cd0d 100644 --- a/src/Quaternions.jl +++ b/src/Quaternions.jl @@ -14,7 +14,6 @@ module Quaternions include("Quaternion.jl") include("Octonion.jl") - include("DualQuaternion.jl") export Quaternion, QuaternionF16, @@ -23,26 +22,18 @@ module Quaternions Octonion, OctonionF16, OctonionF32, - OctonionF64, - DualQuaternion, - DualQuaternionF16, - DualQuaternionF32, - DualQuaternionF64 + OctonionF64 export quat export octo export imag_part - export dualquat export angleaxis export angle export axis export normalize export normalizea - export dconj export quatrand export nquatrand export octorand - export dualquatrand - export ndualquatrand export qrotation export rotationmatrix export slerp diff --git a/test/DualQuaternion.jl b/test/DualQuaternion.jl deleted file mode 100644 index e56c8304..00000000 --- a/test/DualQuaternion.jl +++ /dev/null @@ -1,432 +0,0 @@ -using Quaternions -using DualNumbers -using ForwardDiff -using LinearAlgebra -using Random -using Test - -function fdquat_to_dualquat(q::Quaternion{<:ForwardDiff.Dual}) - q0 = quat( - ForwardDiff.value(q.s), - ForwardDiff.value(q.v1), - ForwardDiff.value(q.v2), - ForwardDiff.value(q.v3), - ) - qe = quat( - ForwardDiff.partials(q.s, 1), - ForwardDiff.partials(q.v1, 1), - ForwardDiff.partials(q.v2, 1), - ForwardDiff.partials(q.v3, 1), - ) - return dualquat(q0, qe) -end - -function dualquat_to_fdquat(dq::DualQuaternion) - return quat( - ForwardDiff.Dual(dq.q0.s, dq.qe.s), - ForwardDiff.Dual(dq.q0.v1, dq.qe.v1), - ForwardDiff.Dual(dq.q0.v2, dq.qe.v2), - ForwardDiff.Dual(dq.q0.v3, dq.qe.v3), - ) -end - -@testset "DualQuaternion" begin - @testset "type aliases" begin - @test DualQuaternionF16 === DualQuaternion{Float16} - @test DualQuaternionF32 === DualQuaternion{Float32} - @test DualQuaternionF64 === DualQuaternion{Float64} - end - - @testset "Constructors" begin - @testset "from quaternions" begin - q0 = Quaternion{Int}(1, 2, 3, 4, false) - qe = QuaternionF32(5, 6, 7, 8, false) - @test @inferred(DualQuaternion(q0, qe)) isa DualQuaternionF32 - @test DualQuaternion(q0, qe) === DualQuaternionF32( - QuaternionF32(1, 2, 3, 4, false), QuaternionF32(5, 6, 7, 8, false), false - ) - @test @inferred(DualQuaternionF64(q0, qe, false)) === DualQuaternionF64( - QuaternionF64(1, 2, 3, 4, false), QuaternionF64(5, 6, 7, 8, false), false - ) - @test DualQuaternionF64(q0, qe, true) === DualQuaternionF64( - QuaternionF64(1, 2, 3, 4, false), QuaternionF64(5, 6, 7, 8, false), true - ) - - qnorm = Quaternion(0//1, 1//1, 0//1, 0//1, true) - @test @inferred(DualQuaternion(q0)) isa DualQuaternion{Int} - @test DualQuaternion(qnorm) === DualQuaternion(qnorm, zero(qnorm), true) - @test @inferred(DualQuaternionF64(q0)) === - DualQuaternionF64(QuaternionF64(q0), zero(QuaternionF64), false) - end - @testset "from dual" begin - coef = (dual(1, 2), dual(3.0, 4.0), dual(5.0f0, 6.0f0), dual(7//1, 8//1)) - @testset for T in (Float32, Float64, Int), norm in (true, false) - dq = @inferred DualQuaternion{T}(coef..., norm) - @test dq isa DualQuaternion{T} - @test dq.norm === norm - @test dq === DualQuaternion{T}(convert.(Dual{T}, coef)..., norm) - @test dq == DualQuaternion( - Quaternion(DualNumbers.value.(coef)...), - Quaternion(DualNumbers.epsilon.(coef)...), - ) - dq2 = @inferred DualQuaternion(convert.(Dual{T}, coef)..., norm) - @test DualQuaternion(convert.(Dual{T}, coef)..., norm) === dq - @test DualQuaternion(coef[1]) == - DualQuaternion(coef[1], fill(zero(coef[1]), 3)...) - @test DualQuaternion{T}(coef[1]) == - DualQuaternion(convert(Dual{T}, coef[1]), fill(zero(Dual{T}), 3)...) - if !norm - @test DualQuaternion(convert.(Dual{T}, coef)...) === dq - end - end - end - @testset "from real" begin - @testset for x in (-1//1, 1.0, 2.0), T in (Float32, Float64, Int, Rational{Int}) - coef = (Quaternion{T}(x, 0, 0, 0, isone(abs(x))), zero(Quaternion{T})) - @test @inferred(DualQuaternion{T}(x)) === - DualQuaternion{T}(coef..., isone(abs(x))) - @test @inferred(DualQuaternion(T(x))) === - DualQuaternion{T}(coef..., isone(abs(x))) - end - end - @testset "from dual quaternion" begin - dq = DualQuaternion( - QuaternionF32(1, 2, 3, 4, false), QuaternionF32(4, 5, 6, 7, false), false - ) - @test @inferred(DualQuaternion(dq)) === dq - @test @inferred(DualQuaternionF64(dq)) === - DualQuaternionF64(dq.q0, dq.qe, false) - end - @testset "from vector" begin - v = randn(3) - @test @inferred(DualQuaternion(v)) === - DualQuaternion(Quaternion(zero(v)), Quaternion(v)) - end - end - - @testset "==" begin - @test DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) == - DualQuaternion(Quaternion(1.0, 2, 3, 4), Quaternion(5, 6, 7, 8)) - @test DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) != - DualQuaternion(Quaternion(1.0, 2, 3, 4), Quaternion(1, 2, 3, 4)) - end - - @testset "deprecated warning" begin - @test_deprecated DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) - @test_deprecated DualQuaternion{Int}(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8), false) - @test_deprecated DualQuaternion{Float64}(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8), false) - end - - @testset "convert" begin - @test convert(DualQuaternion{Float64}, 1) === DualQuaternion(1.0) - @test convert(DualQuaternion{Float64}, DualNumbers.Dual(1, 2)) === - DualQuaternion(Quaternion(1.0), Quaternion(2.0)) - @test convert(DualQuaternion{Float64}, Quaternion(1, 2, 3, 4)) === - DualQuaternion(Quaternion(1.0, 2.0, 3.0, 4.0)) - @test convert( - DualQuaternion{Float64}, - DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)), - ) === DualQuaternion(Quaternion(1.0, 2.0, 3.0, 4.0), Quaternion(5.0, 6.0, 7.0, 8.0)) - end - - @testset "promote" begin - @test promote(DualQuaternion(1.0), 1.0) === - (DualQuaternion(1.0), DualQuaternion(1.0)) - @test promote(DualQuaternion(1.0f0), 2.0) === - (DualQuaternion(1.0), DualQuaternion(2.0)) - @test promote(DualQuaternion(1.0f0), dual(1, 2)) === - (DualQuaternion(1.0f0), DualQuaternion(dual(1.0f0, 2.0f0))) - @test promote(DualQuaternion(1.0f0), Quaternion(3//1)) === - (DualQuaternion(1.0f0), DualQuaternion(3.0f0)) - @test promote(DualQuaternion(1.0f0), DualQuaternion(2.0)) === - (DualQuaternion(1.0), DualQuaternion(2.0)) - - @test Quaternion(1, 2, 3, 4) == DualQuaternion(Quaternion(1, 2, 3, 4)) - @test Quaternion(1, 2, 3, 4) != - DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) - @test DualQuaternion(1) == 1.0 - end - - @testset "shorthands" begin - @test dualquat(Quaternion(1, 0, 0, 0)) == Quaternion(1, 0, 0, 0) - @test dualquat(Quaternion(1, 2, 3, 4)) == Quaternion(1, 2, 3, 4) - @test dualquat(Quaternion(1, 0, 0, 0)) === DualQuaternion(Quaternion(1, 0, 0, 0)) # checking the .norm field in particular - @test dualquat(Quaternion(1, 2, 3, 4)) === DualQuaternion(Quaternion(1, 2, 3, 4)) - @test dualquat(1) === DualQuaternion(1) - @test dualquat(Dual(1, 2)) === DualQuaternion(Dual(1, 2)) - @test dualquat(Dual(1, 2), Dual(0), Dual(0), Dual(0)) === - DualQuaternion(Dual(1, 2), Dual(0), Dual(0), Dual(0)) - @test dualquat(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) == - DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) - @test dualquat(Quaternion(1, 0, 0, 0), Quaternion(0)).norm == false - @test dualquat(Quaternion(1, 0, 0, 0), Quaternion(0), false).norm == false # respect the .norm input (even if wrong) - @test dualquat(Quaternion(1, 2, 3, 4), Quaternion(0), true).norm == true # respect the .norm input (even if wrong) - @test dualquat(Dual(2, 0), Dual(0), Dual(0), Dual(0), true).norm == true # respect the .norm input (even if wrong) - @test dualquat(Dual(1, 0), Dual(0), Dual(0), Dual(0), false).norm == false # respect the .norm input (even if wrong) - end - - @testset "random generation" begin - @testset "dualquatrand" begin - dq = dualquatrand() - @test dq isa DualQuaternionF64 - @test !dq.norm - end - - @testset "ndualquatrand" begin - dq = ndualquatrand() - @test dq isa DualQuaternionF64 - @test dq.norm - end - - @testset "rand($H)" for H in (DualQuaternionF32, DualQuaternionF64) - rng = Random.MersenneTwister(42) - dq = rand(rng, H) - @test dq isa H - @test !dq.norm - - dqs = rand(rng, H, 1000) - @test eltype(dqs) === H - @test length(dqs) == 1000 - xs = map(dqs) do dq - return [ - real(dq.q0) - imag_part(dq.q0)... - real(dq.qe) - imag_part(dq.qe)... - ] - end - xs_mean = sum(xs) / length(xs) - xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1) - @test all(isapprox.(xs_mean, 0.5; atol=0.1)) - @test all(isapprox.(xs_var, 1 / 12; atol=0.01)) - end - end - - @testset "basic" begin - q = rand(DualQuaternionF64) - qnorm = normalize(q) - @test_throws MethodError imag(q) - @test_throws MethodError Quaternions.imag(q) - @test_throws MethodError imag_part(q) - @test conj(q) === dualquat(conj(q.q0), conj(q.qe), q.norm) - @test conj(qnorm) === dualquat(conj(qnorm.q0), conj(qnorm.qe), qnorm.norm) - @test conj(conj(q)) === q - @test conj(conj(qnorm)) === qnorm - @test float(dualquat(dual.(1:4, 5:8)...)) === dualquat(dual.(1.0:4.0, 5.0:8.0)...) - @test dconj(q) === dualquat(q.q0, -q.qe, q.norm) - @test dconj(qnorm) === dualquat(qnorm.q0, -qnorm.qe, qnorm.norm) - @test Quaternions.Q0(q) === q.q0 - @test Quaternions.Qe(q) === q.qe - end - - @testset "algebraic properties" begin - @testset "addition/subtraction" begin - for _ in 1:100 - dq1, dq2 = rand(DualQuaternionF64, 2) - @test (dq1 + dq2).q0 ≈ dq1.q0 + dq2.q0 - @test (dq1 + dq2).qe ≈ dq1.qe + dq2.qe - @test (dq1 - dq2).q0 ≈ dq1.q0 - dq2.q0 - @test (dq1 - dq2).qe ≈ dq1.qe - dq2.qe - @test (+dq1).q0 === dq1.q0 - @test (+dq1).qe === dq1.qe - @test (-dq1).q0 === -dq1.q0 - @test (-dq1).qe === -dq1.qe - end - end - @testset "division" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - dq2 = dq / dq - @test dq2.q0 ≈ 1 - @test dq2.qe ≈ zero(dq2.qe) atol = 1e-6 - dq2 = dq \ dq - @test dq2.q0 ≈ 1 - @test dq2.qe ≈ zero(dq2.qe) atol = 1e-6 - s = randn() - @test (dq / s).q0 ≈ dq.q0 / s - @test (dq / s).qe ≈ dq.qe / s - d = dual(randn(2)...) - @test (dq / d).q0 ≈ (dq / dualquat(d)).q0 - @test (dq / d).qe ≈ (dq / dualquat(d)).qe - end - end - @testset "multiplication is associative" begin - for _ in 1:100 - dq1, dq2, dq3 = rand(DualQuaternionF64, 3) - dq4 = (dq1 * dq2) * dq3 - dq5 = dq1 * (dq2 * dq3) - @test dq4.q0 ≈ dq5.q0 - @test dq4.qe ≈ dq5.qe - end - end - @testset "multiplication by dual" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - d = dual(randn(2)...) - @test (dq * d).q0 ≈ (dq * dualquat(d)).q0 - @test (dq * d).qe ≈ (dq * dualquat(d)).qe - @test (d * dq).q0 ≈ (dualquat(d) * dq).q0 - @test (d * dq).qe ≈ (dualquat(d) * dq).qe - end - end - end - - @testset "iszero" begin - @test iszero(dualquat(0)) - @test !iszero(dualquat(1)) - @test !iszero(dualquat(quat(0, 1, 0, 0))) - @test !iszero(dualquat(quat(0, 0, 1, 0))) - @test !iszero(dualquat(quat(0, 0, 0, 1))) - @test !iszero(dualquat(quat(0), quat(1))) - @test !iszero(dualquat(quat(0), quat(0, 1, 0, 0))) - @test !iszero(dualquat(quat(0), quat(0, 0, 1, 0))) - @test !iszero(dualquat(quat(0), quat(0, 0, 0, 1))) - end - - @testset "isone" begin - @test isone(dualquat(1)) - @test !isone(dualquat(-1)) - @test !isone(dualquat(quat(0, 1, 0, 0))) - @test !isone(dualquat(quat(1, 1, 0, 0))) - @test !isone(dualquat(quat(1, 0, 1, 0))) - @test !isone(dualquat(quat(1, 0, 0, 1))) - @test !isone(dualquat(quat(0), quat(1))) - @test !isone(dualquat(quat(0), quat(0, 1, 0, 0))) - @test !isone(dualquat(quat(0), quat(0, 0, 1, 0))) - @test !isone(dualquat(quat(0), quat(0, 0, 0, 1))) - end - - @testset "^" begin - @testset "^(::DualQuaternion, ::Real)" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - @test dq^2 === dq * dq - @test dq^1 === dq - @test_broken (dq^2.0).q0 ≈ (dq * dq).q0 - @test_broken (dq^2.0).qe ≈ (dq * dq).qe - @test_broken (dq^1.0).q0 ≈ dq.q0 - @test_broken (dq^1.0).qe ≈ dq.qe - @test_broken (dq^-1.0).q0 ≈ inv(dq).q0 - @test_broken (dq^-1.0).qe ≈ inv(dq).qe - for p in (1.3, 7.8, 1.3f0, 7.8f0) - @test (dq^p).q0 ≈ exp(p * log(dq)).q0 - @test (dq^p).qe ≈ exp(p * log(dq)).qe - end - end - end - @testset "^(::DualQuaternion, ::DualQuaternion)" begin - for _ in 1:100 - q, p = randn(QuaternionF64, 2) - @test q^p ≈ exp(p * log(q)) - end - end - end - - @testset "non-analytic functions" begin - unary_funs = [conj, abs, abs2, norm, sign] - # since every dual quaternion is conjugate to a dual complex number, - # one can establish correctness as follows: - @testset for fun in unary_funs - for _ in 1:100 - dq1, dq2 = rand(DualQuaternionF64, 2) - c = dual(randn(ComplexF64), randn(ComplexF64)) - dq = dualquat(reim(c)..., reim(zero(c))...) - p = dualquat(@inferred(fun(dq))) - @test p.q0 ≈ DualNumbers.value(fun(c)) - @test p.qe ≈ DualNumbers.epsilon(fun(c)) - p2 = dq2 * fun(dq1) * inv(dq2) - p3 = dualquat(fun(dq2 * dq1 * inv(dq2))) - @test p2.q0 ≈ p3.q0 - @test p2.qe ≈ p3.qe - end - end - end - - @testset "analytic functions" begin - unary_funs = [sqrt, inv, exp, log] - # since every dual quaternion is conjugate to a dual complex number, - # one can establish correctness as follows: - @testset for fun in unary_funs - (fun === log || fun === sqrt) && continue # log is currently broken - for _ in 1:100 - dq1, dq2 = rand(DualQuaternionF64, 2) - c = dual(randn(ComplexF64), randn(ComplexF64)) - dq = dualquat(reim(c)..., reim(zero(c))...) - p = dualquat(@inferred(fun(dq))) - @test p.q0 ≈ DualNumbers.value(fun(c)) - @test p.qe ≈ DualNumbers.epsilon(fun(c)) - p2 = dq2 * fun(dq1) * inv(dq2) - p3 = dualquat(fun(dq2 * dq1 * inv(dq2))) - @test p2.q0 ≈ p3.q0 - @test p2.qe ≈ p3.qe - end - end - - @testset "identities" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - @test (inv(dq) * dq).q0 ≈ one(dq.q0) - @test (inv(dq) * dq).qe ≈ zero(dq.qe) atol = 1e-6 - @test (dq * inv(dq)).q0 ≈ one(dq.q0) - @test (dq * inv(dq)).qe ≈ zero(dq.qe) atol = 1e-6 - @test_broken (sqrt(dq) * sqrt(dq)).q0 ≈ dq.q0 - @test_broken (sqrt(dq) * sqrt(dq)).qe ≈ dq.qe - @test_broken exp(log(dq)).q0 ≈ dq.q0 - @test_broken exp(log(dq)).qe ≈ dq.qe - @test exp(zero(dq)).q0 ≈ one(dq.q0) - @test exp(zero(dq)).qe ≈ zero(dq.qe) atol = 1e-6 - @test log(one(dq)).q0 ≈ zero(dq.q0) atol = 1e-6 - @test log(one(dq)).qe ≈ zero(dq.qe) atol = 1e-6 - end - @test_broken log(zero(DualQuaternionF64)) === dualquat(-Inf) - end - end - - @testset "normalize" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - q = dualquat_to_fdquat(dq) - dqnorm = normalize(dq) - qnorm = normalize(q) - dqnorm2 = fdquat_to_dualquat(qnorm) - @test dqnorm.q0 ≈ dqnorm2.q0 - @test dqnorm.qe ≈ dqnorm2.qe - @test normalize(dqnorm) === dqnorm - @test iszero(normalize(zero(dq))) - end - end - - @testset "normalizea" begin - for _ in 1:100 - dq = rand(DualQuaternionF64) - q = dualquat_to_fdquat(dq) - dqnorm, dqa = normalizea(dq) - qnorm, qa = normalizea(q) - dqnorm2 = fdquat_to_dualquat(qnorm) - @test dqnorm.q0 ≈ dqnorm2.q0 - @test dqnorm.qe ≈ dqnorm2.qe - @test DualNumbers.value(dqa) ≈ ForwardDiff.value(qa) - @test DualNumbers.epsilon(dqa) ≈ ForwardDiff.partials(qa, 1) - end - end - - @testset "angle/axis/angleaxis" begin - dq1 = rand(DualQuaternionF64) - dq2 = dq1 + (1 - dq1.q0.s) - for dq in (dq1, dq2) - t1, p1 = angleaxis(dq) - t2, v = angleaxis(dualquat_to_fdquat(dq)) - p2 = fdquat_to_dualquat(quat(v)) - @test DualNumbers.value(t1) ≈ ForwardDiff.value(t2) - @test_broken DualNumbers.epsilon(t1) ≈ ForwardDiff.partials(t2, 1) - @test p2.q0 ≈ p1.q0 - @test_broken p2.qe ≈ p1.qe - t3 = angle(dq) - @test DualNumbers.value(t3) ≈ DualNumbers.value(t1) - @test DualNumbers.epsilon(t3) ≈ DualNumbers.epsilon(t1) - p3 = axis(dq) - @test p3.q0 ≈ p1.q0 - @test p3.qe ≈ p1.qe - end - end -end diff --git a/test/runtests.jl b/test/runtests.jl index d1c2ecd8..4b04375d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,4 +4,3 @@ using Quaternions include("helpers.jl") include("Quaternion.jl") include("Octonion.jl") -include("DualQuaternion.jl")