Miscellany

exampleArray = {{6, 8, -20}, {3, 4, 9}, {2, -27, 12}};
exampleArgMin = argMin[exampleArray]
Apply[Part, Join[{exampleArray}, exampleArgMin]]

To simultaneously superimpose two or more finite deformations, you take their logs, add them, and exponentiate. Here's an example that gives the corrected version of inclined transpression a la Jones and Holdsworth (1998).

Clear[gxy, gzy, az, t];
a = {{1, gxy, 0}, {0, 1, 0}, {0, 0, 1}};
b = {{1, 0, 0}, {0, 1, 0},{0, gzy, 1}};
c = {{1, 0, 0}, {0, 1 / az, 0}, {0, 0, az}};
l = matrixLog[a] + matrixLog[b] + matrixLog[c];
MatrixExp[t l] // MatrixForm

Here are kinematic vorticity and octahedral shear strain of dextral monoclinic transpression.

Simplify[kinematicVorticity[monoclinicTranspression[g, k]], {g > 0, 1 > k > 0}]
Simplify[
	octahedralShearStrain[
		MatrixExp[monoclinicTranspression[g, k]].
		Transpose[MatrixExp[monoclinicTranspression[g, k]]]],
	{g > 0, 1 > k > 0}]

For example, here is simple shear striking 90 degrees (east) and striking 180 degrees (south).

MatrixExp[simpleShear[g]] // MatrixForm
MatrixExp[striking[180 Degree, simpleShear[g]]] // MatrixForm

Here's an example of plotting data. The filled data points are explicitly made up, in Cartesian coordinates. The unfilled data points are randomly generated in spherical coordinates. Two silly lines are also drawn, for no good reason.

black = {{1, 3, 2}, {2, 2, 2}, {5, 1, 3}, {4, 6, 2}, {1, 0, 2}, {0, 3, 1}, {2, 5, 0}, {3, 7, 4}, {-1, 2, 0}, {2, 5, 1}, {-2, 5, 2}, {3, 1, 1}, {0, -1, 3}, {1, -2, 1}, {-2, -1, -1}};
black = Map[lowerHemisphere, black];
white = Table[{1, RandomReal[{0, Pi}], RandomReal[{0, 2 Pi}]}, {i, 1, 10}];
white = Map[cartesianFromSpherical, white];
hemisphericalPlot[white, black, {{{1, 2, 3}, {4, 5, 6}}, {{-1, -2, 3}, {-4, -5, 6}, {-7, 8, 9}}}]

While we're at it, we might as well do Wulff and Schmidt nets. The routines below take two angles in radians. The first is the grid spacing. The second is an angle that determines how many samples to use in constructing the curves that make up the net. Make it equal to the grid spacing to get a net made up of straight lines; if your grid spacing is large, then this may look crude, so change the second angle to a smaller number.

lowerHemisphereGrid[gridSkip_, sampleSkip_] :=
	Module[{curves, curve, phi, theta},
		curves = {};
		(* Construct the parallels. *)
		phi = gridSkip;
		While[phi < Pi,
			curve = {};
			theta = 0;
			While[theta ≤ Pi,
				curve = Append[curve, {Sin[phi] Cos[theta], Cos[phi], -Sin[phi] Sin[theta]}];
				theta += sampleSkip];
			curves = Append[curves, curve];
			phi += gridSkip];
		(* Construct the meridians, including the two boundary meridians: *)
		theta = 0;
		While[theta <= Pi,
			curve = {};
			phi = 0;
			While[phi <= Pi,
				curve = Append[curve, {Sin[phi] Cos[theta], Cos[phi], -Sin[phi] Sin[theta]}];
				phi += sampleSkip];
			curves = Append[curves, curve];
	theta += gridSkip];
		curves];

wulffNet[gridSkip_, sampleSkip_] :=
hemisphericalPlot[{}, {}, lowerHemisphereGrid[gridSkip, sampleSkip], projection -> equalAngleProjection];

schmidtNet[gridSkip_, sampleSkip_] :=
hemisphericalPlot[{}, {}, lowerHemisphereGrid[gridSkip, sampleSkip]];

Here is an example:

GraphicsGrid[{{wulffNet[10 Degree, 5 Degree], schmidtNet[10 Degree, 5 Degree]}}]

Here’s a demonstration of Robin’s statistic using exact data.

randomEllipsoid[{a_, b_}] :=
	Module[{as, eTilde, q},
		as = Table[RandomReal[{a, b}], {3}];
		eTilde = determinantOne[DiagonalMatrix[as^-2]];
		q = rotationUniform[];
		Transpose[q].eTilde.q];
randomSection[e_] :=
	Module[{q, eTilde, eig, as, q1, q2},
		q = rotationUniform[];
		eTilde = q.e.Transpose[q];
		eig = Eigensystem[{{eTilde[[1, 1]], eTilde[[1, 2]]}, {eTilde[[2, 1]], eTilde[[2, 2]]}}];
		as = eig[[1]]^(-1 / 2);
		q1 = eig[[2, 1, 1]] q[[1]] + eig[[2, 1, 2]] q[[2]];
		q2 = eig[[2, 2, 1]] q[[1]] + eig[[2, 2, 2]] q[[2]];
		{q1, q2, as[[1]], as[[2]]}];
e = randomEllipsoid[{0.25, 4}]
sections = Table[randomSection[e], {10}];
ePred = robin[sections]
b = 3;
contour = ContourPlot3D[{x, y, z}.e.{x, y, z} == 1, {x, -b, b}, {y, -b, b}, {z, -b, b}, Contours->{1}];
curves = Table[
	ParametricPlot3D[sections[[i, 3]] Cos[t] sections[[i, 1]] + sections[[i, 4]] Sin[t] sections[[i, 2]], {t, 0, 2 Pi}, PlotStyle->Thick], {i, 1, Length[sections]}];
Apply[Show, Join[{contour}, curves]]
contour = ContourPlot3D[{x, y, z}.ePred.{x, y, z} == 1, {x, -b, b}, {y, -b, b}, {z, -b, b}, Contours->{1}];
curves = Table[
	ParametricPlot3D[sections[[i, 3]] Cos[t] sections[[i, 1]] + sections[[i, 4]] Sin[t] sections[[i, 2]], {t, 0, 2 Pi}, PlotStyle->Thick], {i, 1, Length[sections]}];
Apply[Show, Join[{contour}, curves]]

Here we see the BHEP statistic saying that the normal distribution is more normal than others. The difference is one or two orders of magnitude.

sigma = Table[Table[RandomReal[{-0.1, 0.1}], {6}], {6}];
sigma = (sigma + Transpose[sigma]) / 2;
sigma = MatrixExp[sigma];
bhep[Table[RandomVariate[MultinormalDistribution[{0, 0, 0, 0, 0, 0}, sigma]], {100}], 0.5]
bhep[Table[E^Table[RandomReal[{-10, 10}], {6}], {100}], 0.5]
bhep[Table[Table[RandomVariate[ChiSquareDistribution[0.5]], {6}], {100}], 0.5]
bhep[Table[Table[RandomVariate[FRatioDistribution[2, 10]], {6}], {100}], 0.5]
bhep[Table[Table[RandomVariate[HotellingTSquareDistribution[10, 15]], {6}], {100}], 0.5]

Here we pick a random 6x6 matrix, check that it is nonsingular, transform the data by it, and check that the transformation has no effect on the BHEP statistic. Which is easy to see in the formula, isn’t it?

a = Table[Table[RandomReal[{-5, 5}], {6}], {6}];
Det[a]
data = Table[RandomVariate[MultinormalDistribution[{0, 0, 0, 0, 0, 0}, sigma]], {100}];
bhep[data, 0.5] - bhep[Transpose[a.Transpose[data]], 0.5] // Chop

As promised, the Frobenius inner product of log ellipsoids equals the Euclidean inner product of their vector versions.

logE1 = ellipsoidLog[ellipsoidUniformSemiaxes[-5, 5]];
logE2 = ellipsoidLog[ellipsoidUniformSemiaxes[-5, 5]];
frobenius[logE1, logE2]
Dot[ellipsoidVectorFromLog[logE1], ellipsoidVectorFromLog[logE2]]

Here we compute spreads of ellipsoids with tighter and tighter uniform semiaxis distributions. The spread is just about four times the square of the log of the maximum allowed semi-axis size. We could probably work out why, if we cared.

ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[-5, 5], {1000}]][[2]]
ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[-4, 4], {1000}]][[2]]
ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[-3, 3], {1000}]][[2]]
ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[-2, 2], {1000}]][[2]]
ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[-1, 1], {1000}]][[2]]
ellipsoidMeanSpread[Table[ellipsoidUniformSemiaxes[0, 0], {1000}]][[2]]

Notice that the BHEP statistic climbs as the semiaxis distributions get tighter. Perhaps this is because more of the ellipsoid variation comes from the orientation, which is uniform rather than normal? But the semiaxes are uniform too. Hmm.

ellipsoids = Table[ellipsoidUniformSemiaxes[-5, 5], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]
ellipsoids = Table[ellipsoidUniformSemiaxes[-4, 4], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]
ellipsoids = Table[ellipsoidUniformSemiaxes[-3, 3], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]
ellipsoids = Table[ellipsoidUniformSemiaxes[-2, 2], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]
ellipsoids = Table[ellipsoidUniformSemiaxes[-1, 1], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]
ellipsoids = Table[ellipsoidUniformSemiaxes[0,0], {100}];
ellipsoidBHEP[ellipsoids, {0.1, 0.5, 1.0, 10.0, 100.0, 1000.0}]

sigma = Table[Table[RandomReal[{-0.01, 0.01}], {6}], {6}];
sigma = MatrixExp[(sigma + Transpose[sigma]) / 2];
ellipsoidBHEP[Table[ellipsoidNormalSemiaxes[0, 1], {100}], 0.5]
ellipsoidBHEP[Table[ellipsoidUniformSemiaxes[-1, 1], {100}], 0.5]
ellipsoidBHEP[Table[ellipsoidUniformNorm[2], {100}], 0.5]
ellipsoidBHEP[Table[ellipsoidNormal[{0, 0, 0, 0, 0, 0}, sigma], {100}], 0.5]
bimodal = Join[
	Table[ellipsoidNormal[10 {1, 1, 1, 1, 1, 1}, sigma], {50}],
	Table[ellipsoidNormal[10 {1, 1, 1, -1, -1, -1}, sigma], {25}],
	Table[ellipsoidNormal[10 {-1, -1, -1, 1, 1, 1}, sigma], {25}]];
ellipsoidBHEP[bimodal, 0.5]

Here we build some graphics as lists of curves. A curve is a list of points (to be joined by straight lines). box3D has some redundant lines in it.

box2D = {{{1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 1}}, 
{{1, 0}, {0, 0}}, {{0, 1}, {0, 0}}};
circlePoints[{a_, b_}, r_, n_] := 
Table[{r Cos[t] + a, r Sin[t] + b}, {t, 0, 2 Pi, 2 Pi / n}];
circle2D = circlePoints[{0, 0}, 1, 100];
rock2D = Append[box2D, circle2D];

box3D = {{{-1,-1,-1},{-1,-1,1}},{{-1,1,-1},{-1,1,1}},{{1,-1,-1},{1,-1,1}},{{1,1,-1},{1,1,1}},{{-1,-1,-1},{-1,1,-1}},{{-1,-1,1},{-1,1,1}},{{1,-1,-1},{1,1,-1}},{{1,-1,1},{1,1,1}},{{-1,-1,-1},{1,-1,-1}},{{-1,-1,1},{1,-1,1}},{{-1,1,-1},{1,1,-1}},{{-1,1,1},{1,1,1}},
{{0,0,0},{0,0,1}},{{0,1,0},{0,1,1}},{{1,0,0},{1,0,1}},{{0,0,0},{0,1,0}},{{0,0,1},{0,1,1}},{{1,0,0},{1,1,0}},{{0,0,0},{1,0,0}},{{0,0,1},{1,0,1}},{{0,1,0},{1,1,0}}};

Here is a nice example: steady coaxial deformation, and steady coaxial deformation with extraneous rotation.

deformer[3, box3D, {}, MatrixExp[t 
	{{-Log[1.5], 0, 0}, {0, -Log[1.5], 0}, {0, 0, Log[2.25]}}]]
deformer[3, box3D, {}, MatrixExp[t 
	{{-Log[1.5], -2 Pi k, 0}, {2 Pi k, -Log[1.5], 0}, {0, 0, Log[2.25]}} /. k -> 1]]

Here's an example of the hemispherical deformer. The diagonals that appear are actually the vertical edges of the cube. They appear diagonal because of how the top of the cube is mapped to the lower hemisphere.

cube3D = {{{-1,-1,-1},{-1,-1,1}},{{-1,1,-1},{-1,1,1}},{{1,-1,-1},{1,-1,1}},{{1,1,-1},{1,1,1}},{{-1,-1,-1},{-1,1,-1}},{{-1,-1,1},{-1,1,1}},{{1,-1,-1},{1,1,-1}},{{1,-1,1},{1,1,1}},{{-1,-1,-1},{1,-1,-1}},{{-1,-1,1},{1,-1,1}},{{-1,1,-1},{1,1,-1}},{{-1,1,1},{1,1,1}}};

hemisphericalDeformer[cube3D, {}, MatrixExp[t 
	{{-Log[1.5], 0, 0}, {0, -Log[1.5], 0}, {0, 0, Log[2.25]}}]]

hemisphericalDeformer[cube3D, {}, MatrixExp[t 
	{{-Log[1.5], -2 Pi k, 0}, {2 Pi k, -Log[1.5], 0}, {0, 0, Log[2.25]}} /. k -> 1]]

Rotations

Here’s an example, for a normally distributed cluster of rotations about a random center. You can configure the standard deviation of the normal distribution, to make the cluster broader or tighter. If all goes well, then the rotation diameter is less than pi (or even pi / 2), the mean is close to the center (but not necessarily equal), and increasing the standard deviation increases the variance at the end.

center = rotationUniform[]
rotations = Table[
	center.rotationFromAngleAxis[
		{RandomVariate[NormalDistribution[0, 10 Degree]],
		randomSpherePoint[]}],
	{100}];
rotationDiameter[rotations]
mean = rotationMean[rotations, 10^-10, 100]
rotationVariance[rotations, mean[[1]]]

Here’s an alternative method for generating uniformly random rotations. It seems to be slightly slower than the method above.

(* Returns a uniformly chosen rotation. The axis is chosen uniformly *)
(* on the sphere. The angle x is chosen in [0, Pi] from the density *)
(* f(x) = (1 - cos x) / pi (Miles, 1965), using rejection sampling. *)
rotationUniformAlt[] :=
	Module[{x, y},
		x = RandomReal[{0, Pi}];
		y = RandomReal[{0, 2 / Pi}];
		While[y > (1 - Cos[x]) / Pi,
			x = RandomReal[{0, Pi}];
			y = RandomReal[{0, 2 / Pi}]];
		rotationFromAngleAxis[{x, randomSpherePoint[]}]];
rot = rotationUniformAlt[];
rot.Transpose[rot] // Chop
Det[rot]
Timing[Table[rotationUniform[], {1000}]][[1]]
Timing[Table[rotationUniformAlt[], {1000}]][[1]]
Timing[Table[rotationUniform[], {1000}]][[1]]
Timing[Table[rotationUniformAlt[], {1000}]][[1]]
Timing[Table[rotationUniform[], {1000}]][[1]]
Timing[Table[rotationUniformAlt[], {1000}]][[1]]
rotationPlot[Table[rotationUniformAlt[], {1000}]]

Here’s a Kamb contouring example.

a = rotationUniform[];
as = Table[rotationWrappedTrivariateNormal[a, 4], {20}];
b = rotationUniform[];
bs = Table[rotationWrappedTrivariateNormal[b, 3], {20}];
rotationPlot[Join[as, bs]]

DateString[]
rotationKambPlot[Join[as, bs], 3, {3, 6}, 1, 2, 2]
DateString[]

DateString[]
rotationKambPlot[Join[as, bs], 3, {3, 6}, 1, 2, 3]
DateString[]

GraphicsGrid[{{Out[1502], Out[1492]}}]

Here’s a principal geodesic analysis example.

a = rotationUniform[];
as = Table[rotationWrappedTrivariateNormal[a, 4], {20}];
b = rotationUniform[];
bs = Table[rotationWrappedTrivariateNormal[b, 3], {10}];
asbs = Join[as, bs];
center = rotationMean[asbs]
center = center[[1]]
geo = rotationPrincipalComponents[asbs, center];
rotationPlot[asbs, {
	rotationGeodesicSteps[center, geo[[2, 1]], Pi / 12], 
	rotationGeodesicSteps[center, geo[[2, 2]], Pi / 12],
	rotationGeodesicSteps[center, geo[[2, 3]], Pi / 12]}]

rotationPlot[asbs, {
	rotationGeodesicSegmentSteps[
		rotationFromAngleAxis[{20 geo[[1, 1]], geo[[2, 1]]}].center,
		center, Pi / 12],
	rotationGeodesicSegmentSteps[
		rotationFromAngleAxis[{20 geo[[1, 2]], geo[[2, 2]]}].center,
		center, Pi / 12],
	rotationGeodesicSegmentSteps[
		rotationFromAngleAxis[{20 geo[[1, 3]], geo[[2, 3]]}].center,
		center, Pi / 12]}]

Is the Rodrigues plot equal-angle? No. Also, there is apparently no equal-angle map from a ball to SO(3) that takes the vector g(a) * u to the rotation (a, u), for any function g.

Clear[a, f, h];
r = rotationFromAngleAxis[{a, cartesianFromSpherical[{1, f, h}]}];
r // MatrixForm

drda = Simplify[D[r, a]];
drdf = Simplify[D[r, f]];
drdh = Simplify[D[r, h]];
drda // MatrixForm
drdf // MatrixForm
drdh // MatrixForm

v = g[a] {Sin[f] Cos[h], Sin[f] Sin[h], Cos[f]}
dvdx = {Grad[v[[1]], {a, f, h}], Grad[v[[2]], {a, f, h}], Grad[v[[3]], {a, f, h}]}

toSO3 = (1 / 2) {
	{Tr[Transpose[drda].drda], Tr[Transpose[drda].drdf], Tr[Transpose[drda].drdh]},
	{Tr[Transpose[drdf].drda], Tr[Transpose[drdf].drdf], Tr[Transpose[drdf].drdh]},
	{Tr[Transpose[drdh].drda], Tr[Transpose[drdh].drdf], Tr[Transpose[drdh].drdh]}};
toBall = Transpose[dvdx].dvdx;

Simplify[toSO3[[1, 1]]]
Simplify[toBall[[1, 1]]]
Simplify[toSO3[[2, 2]]]
Simplify[toBall[[2, 2]]]
Simplify[toSO3[[3, 3]]]
Simplify[toBall[[3, 3]]]

Another Kamb contouring example.

a = rotationUniform[];
as = Table[rotationWrappedTrivariateNormal[a, 4], {20}];
b = rotationUniform[];
bs = Table[rotationWrappedTrivariateNormal[b, 3], {20}];
rotationPlot[Join[as, bs], {}]

DateString[]
rotationKambPlot[Join[as, bs], 3, {3, 6}, 1, 2, 2]
DateString[]

DateString[]
rotationKambPlot[Join[as, bs], 3, {3, 6}, 1, 2, 4]
DateString[]

q = rotationUniform[];
w = Table[RandomReal[{0, 2.5}], {3}];
Norm[w]
w = {{0, -w[[3]], w[[2]]}, {w[[3]], 0, -w[[1]]}, {-w[[2]], w[[1]], 0}};
m = Table[RandomReal[{0, 2.5}], {3}];
Norm[m]
m = {{0, -m[[3]], m[[2]]}, {m[[3]], 0, -m[[1]]}, {-m[[2]], m[[1]], 0}};
xs = Table[RandomReal[{0, 1}], {10}];
rs = Table[rotationExp[xs[[i]] m].rotationExp[w].q, {i, 1, Length[xs]}];
curve = rotationRegressionCurve[Transpose[{xs, rs}]];
curve[[1]] - m // Chop
curve[[2]] - rotationExp[w].q // Chop

ODEs

Here is an example of the Euclidean methods, the left-invariant Lie group methods, and the right-invariant Lie group methods.

g = {{1.6, 0.0, 0.0}, {0.0, 2.0, 3.0}, {0.0, 0.125, 0.5}};
x = Function[g, {{-0.3, 0.1, 0.0}, {0.0, -1.1, -0.2}, {0.3, 0.2, 1.4}}];
euclideanEuler[g, x, 10] // MatrixForm
euclideanRungeKutta[g, x, 10] // MatrixForm
g + x[0] // MatrixForm

Here are examples with rigid spheroids and ellipsoids.

l = traceZero[Table[Table[RandomReal[{-2, 2}], {3}], {3}]]
Flatten[Transpose[l]]
u0 = randomSpherePoint[]
b = RandomReal[{-1, 1}]
d = 0.5 (l + Transpose[l]);
w = 0.5 (l - Transpose[l]);
m = b d + w;
x[g_] :=
	m.g - (g.m.g) g;
u1 = euclideanRungeKutta[u0, x, 10]

l = traceZero[Table[Table[RandomReal[{-2, 2}], {3}], {3}]];
Flatten[Transpose[l]]
q0 = rotationUniform[];
Flatten[Transpose[q0]]
a = Table[RandomReal[{0.25, 4}], {3}];
a /= (a[[1]] a[[2]] a[[3]])^(1 / 3)
d = 0.5 (l + Transpose[l]);
w = 0.5 (l - Transpose[l]);
x[g_] :=
	Module[{dTilde, wTilde},
		dTilde = g.d.Transpose[g];
		wTilde12 = dTilde[[1, 2]] (a[[1]]^2 - a[[2]]^2) / (a[[1]]^2 + a[[2]]^2);
		wTilde13 = dTilde[[1, 3]] (a[[1]]^2 - a[[3]]^2) / (a[[1]]^2 + a[[3]]^2);
		wTilde23 = dTilde[[2, 3]] (a[[2]]^2 - a[[3]]^2) / (a[[2]]^2 + a[[3]]^2);
		wTilde = {{0, wTilde12, wTilde13}, {-wTilde12, 0, wTilde23}, {-wTilde13, -wTilde23, 0}};
		wTilde.g - g.w];
q1 = euclideanRungeKutta[q0, x, 10];
Flatten[Transpose[q1]]

l = traceZero[Table[Table[RandomReal[{-2, 2}], {3}], {3}]];
Flatten[Transpose[l]]
q0 = rotationUniform[];
Flatten[Transpose[q0]]
a = Table[RandomReal[{0.25, 4}], {3}];
a /= (a[[1]] a[[2]] a[[3]])^(1 / 3)
e0 = Transpose[q0].DiagonalMatrix[a^-2].q0;
e1 = rigEByQ[e0, l, 10];
q1 = Eigenvectors[e1];
MatrixForm[q1]

l = traceZero[Table[Table[RandomReal[{-2, 2}], {3}], {3}]];
Flatten[Transpose[l]]
q0 = rotationUniform[];
Flatten[Transpose[q0]]
a = Table[RandomReal[{0.25, 4}], {3}];
a /= (a[[1]] a[[2]] a[[3]])^(1 / 3)
e0 = Transpose[q0].DiagonalMatrix[a^-2].q0;
e1 = rigEByQ[e0, l, 10];
q1 = Eigenvectors[e1];
MatrixForm[q1]

q = rotationUniform[]
a = Table[RandomReal[{-4, 4}], {3}];
e0Tilde = determinantOne[DiagonalMatrix[a^-2]]
e0 = Transpose[q].e0Tilde.q;
l = traceZero[Table[Table[RandomReal[{-2, 2}], {3}], {3}]]
e1 = ellipsoidRigid[e0, l, 10]
Det[e1]