Add support for `^(::Rotation, ::Real)`

[lieskjur]

Hi, recently I've found a need for scaling rotations described by unit quaternions, by which I mean produce a rotation along the same axis only with a different angle. I've written a small function for this purpose

function scale(q::Q,w::Real) where Q <: UnitQuaternion
	θ = rotation_angle(q)
	axis = rotation_axis(q)
	s,c = sincos(w*θ/2)
	return Q(c,s*axis[1],s*axis[2],s*axis[3],false)
end

not wanting to override the current function (*)(q::Q, w::Real) where Q<:UnitQuaternion. Would any one else see use in adding such a method to the package?

[andyferris]

Should it just overload ^? I think that’s equivalent to your scale.

[lieskjur]

Didn't know, ^ was already implemented in such a manner. Thank you for the response, will be using it from now on.

[lieskjur]

Further testing it with decimal numbers I have noticed ^ produces complex matrices. Is this behavior intended?

[andyferris]

Right.

^ for AbstractMatrix is defined in LinearAlgebra. These generic methods won't know that the matrix is orthogonal and might come up with complex solutions (or at least solve in a complex space leaving a small imaginary component due to approximation and roundoff errors). (For integer exponents it uses an exponential-by-squaring technique rather than an eigendecomposition strategy, so you stay in the real space).

However what we need is further specialization for our types. This is defined already for e.g. AngleAxis here. We should probably do one of the below:

1. add more methods for the types defined here, such as UnitQuaternion, like you did above.
2. define a generic method on Rotation. This could just take the real part of the LinearAlgebra method, or it might convert to and from an AngleAxis, or similar. This should cover all types correctly, at least.

function Base.:^(r::Rotation{3}, p::Real)
   return convert(typeof(r), convert(AngleAxis, r) ^ p)
end

function Base.:^(r::Rotation{2}, p::Real)
   return convert(typeof(r), convert(Angle2d, r) ^ p)
end

function Base.:^(r::Rotation{3}, p::Integer)
   return convert(typeof(r), convert(AngleAxis, r) ^ p)
end

function Base.:^(r::Rotation{2}, p::Integer)
   return convert(typeof(r), convert(Angle2d, r) ^ p)
end

(the last two are for disambiguation - we could continue to use power-by-squaring here but some method needs to exist).

@lieskjur would you feel up to creating a PR? Contributions are always appreciated :)

[lieskjur]

I would consider a hybrid approach. After some tinkering I managed to shave off a couple ns from scaling quaternion rotations in comparison to the generic method approach (mostly did it out of interest in how it could be done without the dependency).

@inline function Base.:^(q::Q, p::Real) where Q <: UnitQuaternion
	if q.w == 1
		return q
	else
    	sθ2 = sqrt(q.x * q.x + q.y * q.y + q.z * q.z)
    	θ2 = atan(sθ2,q.w)
    	spθ2,cpθ2 = sincos(p*θ2)
    	k = spθ2 / sθ2
		return Q(cpθ2,k*q.x,k*q.y,k*q.z,false)
	end
end

generic approach:

@benchmark q^.3                                                  
BenchmarkTools.Trial: 10000 samples with 873 evaluations.
 Range (min … max):   94.053 ns …  3.281 μs  ┊ GC (min … max): 0.00% … 95.33%
 Time  (median):     121.367 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   121.649 ns ± 87.948 ns  ┊ GC (mean ± σ):  2.12% ±  2.87%

  ▇  ▆  ▃▄  ▃▁      ▁    █▄ ▁ ▅▆    ▂▄     ▁    ▂ ▁  ▁▂▂▁      ▂
  █▂██▆▆██▇▆██▃▅█▆▄▆██▄▅▄██▅█▇██▇▆▇▇██▇▆▇▇▆█▆▆▇███████████▇▇▇▇ █
  94.1 ns       Histogram: log(frequency) by time       163 ns <

 Memory estimate: 48 bytes, allocs estimate: 1.

type specific:

julia> @benchmark q^.3                                                         
BenchmarkTools.Trial: 10000 samples with 898 evaluations.
 Range (min … max):   83.239 ns …  3.203 μs  ┊ GC (min … max): 0.00% … 95.73%
 Time  (median):     112.622 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   117.387 ns ± 82.041 ns  ┊ GC (mean ± σ):  2.07% ±  2.87%

  ▁  ▄  ▂   ▃            ▄    █▁ ▁ ▄▇     ▂           ▁▂▃▂▂    ▂
  █▁▅█▁▅█▃▁▅█▅▃▄▇▃▃▄▇▄▃▄▄█▃▄▆▅██▆█████▇████▇▃▄▆▆▅▅▅▅▄▇████████ █
  83.2 ns       Histogram: log(frequency) by time       145 ns <

 Memory estimate: 48 bytes, allocs estimate: 1.

I'm not sure the increase in performance justifies the duplicit code though. Maybe it would be worth exploring the same for other types but I'd say that is not a major priority based on the minor improvement in speed.

One more caveat is that special methods for euler_types will have to be defined as there is no conversion between them and AngleAxis but that shouldn't be an issue.

would you feel up to creating a PR? Contributions are always appreciated :)

Sure, it will just take some time.

In Julia doc ^ is refered to as the "Exponentiation operator", I suppose that should be reflected in the accompanying docs instead of scaling, right?

[andyferris]

Well in Rotations.jl rotations are represented as matrices, which multiply vectors. The natural “scaling” operator for rotations is ^ - if you multiply the vector twice by a rotation then you rotate twice as far (which is the same as rotating the vector by the rotation squared - R * (R * v) == (R * R) * v).

(On a technical level, we say we are working in the Lie group not the Lie algebra.)

If this isn’t covered in the README then it should be. I don’t think we need additional docstrings though?

[hyrodium]

fixed by #266