The two examples noted below are taken from Moakher (2002) "Means and averaging in the group of rotations", and correspond to cases where the approximate and exact means can be solved for analytically. We will test that the results found by the Quaternions module
In [15]:
import random
import Quaternions
random.seed(14)
theta_average = 0.
theta_sigma = pi/2.
N = 100
theta = array([random.gauss(theta_average, theta_sigma)]*N)
v = Quaternions.Quaternion(0., random.random()-1, random.random()-1, random.random()-1).normalized()
R = array([Quaternions.exp(-v*(theta_n/2.)) for theta_n in theta])
sumcostheta = sum([cos(theta_n) for theta_n in theta])
sumsintheta = sum([sin(theta_n) for theta_n in theta])
r = sqrt( sumcostheta**2 + sumsintheta**2 )
Theta_a = arctan2(sumsintheta, sumcostheta)
R_a_analytical = Quaternions.exp(-v*(Theta_a/2.))
Theta_g = sum([theta_n for theta_n in theta])/float(N)
R_g_analytical = Quaternions.exp(-v*(Theta_g/2.))
R_a_numerical1 = Quaternions.ApproximateMeanRotor(R)
R_a_numerical2 = Quaternions.ApproximateMeanRotor(R, range(N))
R_g_numerical1 = Quaternions.MeanRotor(R)
R_g_numerical2 = Quaternions.MeanRotor(R, range(N))
R_a_analytical, R_a_numerical1, R_a_numerical2, \
R_g_analytical, R_g_numerical1, R_g_numerical2
Out[15]:
The fact that these are all the same is good: All methods give identical results for the average value of a constant.
The first example uses an arbitrary series of rotations in a one-parameter subgroup of $SO(3)$. Let $\hat{v}$ be an arbitrary unit vector in $\mathbb{R}^3$, and let $\{\theta_n\}$ be a collection of angles we will be rotating by. We will choose them randomly distributed by a Gaussian around 0, with standard deviation small compared to $\pi$. Moakher shows that the solution for the "projected arithmetic mean" (Quaternions.ApproximateMeanRotor) is given by a rotation about $\hat{v}$ by an angle $\Theta_a$, where
\begin{gather}
r = \sqrt{\left( \sum_n \cos\theta_n \right)^2 + \left( \sum_n \sin\theta_n \right)^2}~, \\
\cos\Theta_a = \frac{1}{r} \sum_n \cos\theta_n~, \\
\sin\Theta_a = \frac{1}{r} \sum_n \sin\theta_n~.
\end{gather}
In [1]:
import random
import Quaternions
random.seed(14)
theta_average = 0.
theta_sigma = pi/2.
N = 10
theta = array([random.gauss(theta_average, theta_sigma) for n in range(N)])
v = Quaternions.Quaternion(0., random.random()-1, random.random()-1, random.random()-1).normalized()
R = array([Quaternions.exp(-v*(theta_n/2.)) for theta_n in theta])
sumcostheta = sum([cos(theta_n) for theta_n in theta])
sumsintheta = sum([sin(theta_n) for theta_n in theta])
r = sqrt( sumcostheta**2 + sumsintheta**2 )
Theta_a = arctan2(sumsintheta, sumcostheta)
R_a_analytical = Quaternions.exp(-v*(Theta_a/2.))
Theta_g = sum([theta_n for theta_n in theta])/float(N)
R_g_analytical = Quaternions.exp(-v*(Theta_g/2.))
R_a_numerical1 = Quaternions.ApproximateMeanRotor(R)
R_a_numerical2 = Quaternions.ApproximateMeanRotor(R, range(N))
R_g_numerical1 = Quaternions.MeanRotor(R)
R_g_numerical2 = Quaternions.MeanRotor(R, range(N))
R_a_analytical, R_a_numerical1, R_a_numerical2, \
R_g_analytical, R_g_numerical1, R_g_numerical2
Out[1]:
Note that the third and sixth rotors represent the numerical algorithm using the time-step sizes, and thus should not agree with the other results.
The numerical and analytical results for the geometric mean agree quite nicely. Unfortunately, the results for the arithmetic mean do not agree. I suspect there is an error in Moakher's formula, though I don't see where.
This example seems to work quite nicely. I suspect that there is just something wrong in Moakher's example 1.
In [4]:
import random
import Quaternions
N = 10
theta = random.random()*pi/2.
alpha = random.random()*pi/2.
beta = array([2*(n-1)*pi/N for n in range(N)])
R = array([Quaternions.exp(Quaternions.Quaternion(0,sin(alpha)*cos(beta_n),sin(alpha)*sin(beta_n),cos(alpha))*(theta/2))
for beta_n in beta])
Theta_a = 2*arctan(cos(alpha)*tan(theta/2))
R_a_analytical = Quaternions.exp(Quaternions.zHat*(Theta_a/2.))
R_a_numerical1 = Quaternions.ApproximateMeanRotor(R)
R_a_numerical2 = Quaternions.ApproximateMeanRotor(R, range(N))
R_g_analytical = R_a_analytical
R_g_numerical1 = Quaternions.MeanRotor(R)
R_g_numerical2 = Quaternions.MeanRotor(R, range(N))
R_a_analytical, R_a_numerical1, R_a_numerical2, \
R_g_analytical, R_g_numerical1, R_g_numerical2
Out[4]:
In [5]:
for R in [R_g_numerical1, R_g_numerical2] :
print((R-R_g_analytical).normsquared())
for R in [R_a_numerical1, R_a_numerical2] :
print((R-R_a_analytical).normsquared())
The basic numerical means agree with the analytical results to within $10^{-32}$. The means that use the time-step information disagree at a level consistent with the fact that integration is done using midpoints rather than abscissas.
We can construct two copies of a funky rotation, and then shift one copy in time and/or rotation, then test to see if the alignment routines can undo the shifts. The funky rotation can be constructed by composing two simple rotations about different axes.
In [6]:
import random
import Quaternions
random.seed(14)
N = 100
v1 = Quaternions.Quaternion(0., random.random()-1, random.random()-1, random.random()-1).normalized()
v2 = Quaternions.Quaternion(0., random.random()-1, random.random()-1, random.random()-1).normalized()
alphadot = 0.084
betadot = alphadot/3.9
Ti = 0.
Tf = 10.
T = linspace(Ti, Tf, num=N)
R = array([Quaternions.exp(v1*(alphadot*t/2.))*Quaternions.exp(v2*(betadot*t/2.))
for t in T])
In [7]:
deltat_input = 0.123
Ra = R
Rb = array([Quaternions.exp(v1*(alphadot*t/2.))*Quaternions.exp(v2*(betadot*t/2.))
for t in T+deltat_input])
T_common = T[(T>=T[0]+deltat_input)&(T<=T[-1]-deltat_input)]
Obviously, this won't be able to eliminate the difference, but we would like to see it improve.
In [8]:
R1 = Quaternions.ApproximateOptimalAlignmentRotor(Ra, Rb, T, Ti+deltat_input, Tf-deltat_input)
print([R1[i] for i in range(4)])
semilogy(T, Quaternions.angle(Quaternions.RDelta(R1*Rb, Ra)), label='With rotation')
semilogy(T, Quaternions.angle(Quaternions.RDelta(Rb, Ra)), label='Without rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
That is indeed a huge improvement over no correction, and it approaches zero right in the middle of the alignment region.
This should actually be able to eliminate all the difference (up to numerics), since the required rotation is just the identity, and ApproximateOptimalAlignmentRotor can find that exactly.
In [9]:
deltat2,R2 = Quaternions.ApproximateOptimalAlignment(4., 6., Ra, T, Rb, T)
print(deltat2, [R2[i] for i in range(4)])
Rb2 = Quaternions.Squad(R2*Rb, T+deltat2, T_common)
Ra2 = Quaternions.Squad(Ra, T, T_common)
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb2, Ra2)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
This is indeed successful at eliminating the difference, finding the correct time offset to within $\sqrt{\epsilon}$, and the correction rotation to roughly the same level. I think that blip at the beginning must come from either a bad Squad or just be an unavoidable issue in the initial value of the Squad.
In [6]:
deltat4,R4 = Quaternions.OptimalAlignment(4., 6., Ra, T, Rb, T)
print(deltat4, [R4[i] for i in range(4)])
Rb4 = Quaternions.Squad(R4*Rb, T+deltat4, T_common)
Ra4 = Quaternions.Squad(Ra, T, T_common)
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb4, Ra4)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
This gives us the same result as ApproximateOptimalAlignment, up to numerics.
In [7]:
Rd = Quaternions.exp(Quaternions.Quaternion(0.0,0.1,2.3,4.5))
Ra = R
Rb = Rd*R
print([Rd[i] for i in range(4)])
Hopefully, this will be able to completely eliminate the rotation.
In [8]:
R1 = Quaternions.ApproximateOptimalAlignmentRotor(Ra, Rb, T, Ti, Tf)
print([R1[i] for i in range(4)])
semilogy(T, Quaternions.angle(Quaternions.RDelta(R1*Rb, Ra)), label='With rotation')
semilogy(T, Quaternions.angle(Quaternions.RDelta(Rb, Ra)), label='Without rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
The printed result is almost exactly the input quantity, and the plot clearly shows that the numerical difference is flirting with zero.
This can only do worse than the above, since we are giving it a certain amount of freedom that isn't necessary. But hopefully it won't do too badly.
In [9]:
deltat2,R2 = Quaternions.ApproximateOptimalAlignment(4., 6., Ra, T, Rb, T)
print(deltat2, [R2[i] for i in range(4)])
Rb2 = Quaternions.Squad(R2*Rb, T+deltat2, T_common)
Ra2 = Quaternions.Squad(Ra, T, T_common)
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb2, Ra2)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
Pretty good.
Again, this can only do worse than optimizing over rotation alone. But hopefully not much worse.
In [10]:
deltat4,R4 = Quaternions.OptimalAlignment(2., 6., Ra, T, Rb, T)
Rb4 = Quaternions.Squad(R4*Rb, T+deltat4, T_common)
Ra4 = Quaternions.Squad(Ra, T, T_common)
print(deltat4, [R4[i] for i in range(4)])
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb4, Ra4)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
Again, pretty darn good.
In [11]:
deltat_input = 0.123
Rd = Quaternions.exp(Quaternions.Quaternion(0.0,0.1,2.3,4.5))
Ra = R
Rb = Rd*array([Quaternions.exp(v1*(alphadot*t/2.))*Quaternions.exp(v2*(betadot*t/2.))
for t in T+deltat_input])
T_common = T[(T>=T[0]+deltat_input)&(T<=T[-1]-deltat_input)]
The results here should be almost identical to what we see in the first case.
In [12]:
R1 = Quaternions.ApproximateOptimalAlignmentRotor(Ra, Rb, T, Ti, Tf)
print([R1[i] for i in range(4)])
semilogy(T, Quaternions.angle(Quaternions.RDelta(R1*Rb, Ra)), label='With rotation')
semilogy(T, Quaternions.angle(Quaternions.RDelta(Rb, Ra)), label='Without rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
Indeed, they are almost identical.
The results here should be much like they were in the time-shifted case.
In [13]:
deltat2,R2 = Quaternions.ApproximateOptimalAlignment(4., 6., Ra, T, Rb, T)
print(deltat2, [R2[i] for i in range(4)])
Rb2 = Quaternions.Squad(R2*Rb, T+deltat2, T_common)
Ra2 = Quaternions.Squad(Ra, T, T_common)
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb2, Ra2)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
As expected, these are much like the results from the time-shifted case.
As above, these results should be much like the results from the time-shifted case.
In [14]:
deltat4,R4 = Quaternions.OptimalAlignment(4., 6., Ra, T, Rb, T)
Rb4 = Quaternions.Squad(R4*Rb, T+deltat4, T_common)
Ra4 = Quaternions.Squad(Ra, T, T_common)
print(deltat4, [R4[i] for i in range(4)])
semilogy(T_common, Quaternions.angle(Quaternions.RDelta(Rb4, Ra4)), label='With rotation')
xlabel('Time')
ylabel(r'$\Phi_\Delta$')
legend();
And again, these are basically the same results.
The alignment routines all behave exactly as I expected. I have every reason, then, to call these "passed tests".