Code to accompany Edelman and Strang: Random Triangle Theory with Geometry and Applications

On January 20, 1884: Lewis Carroll asked for the Probability that a Triangle is acute and then Kendall developed Shape theory during the 20th century:

Our paper revisits this work. The code below specs out the conversion between the representations.

In [118]:
√(x)=sqrt(x)
Δ=[1/√(2) -1/√(2) 0;1/√(6) 1/√(6) -2/√(6)]

#Representations: mat (2x2 matrix M)
#                 svd σ₁,σ₂,θ
#                 tri (3 vector of squared edge lengths)
#                 hemi λ,ϕ  (latitude, longitude)
#                 disk r,ϕ   (polar coordinates)
#                 xyz  (3 vector of cartesian coordinates on the hemisphere of radius 1/2)

# CONVERSIONS:

#mat2...

function mat2svd(M)
(U,(σ₁,σ₂),V)=svd(M)
# Pick the right  θ in [0,π)
θ=acos(V[1,1]*sign(det(M))*sign(V[1,2]))
(σ₁,σ₂,θ)
end

mat2tri(M) = diag( (M*Δ)'*(M*Δ))
mat2hemi(M) = svd2hemi( mat2svd(M)...)
mat2mat(M) = svd2mat(mat2svd(M)...) # Canonical Choice (throw away U)

function mat2disk(M)
M/=normfro(M)
MtM=M'*M
rsc = [ MtM[1,2] ; ((MtM)[1,1]-(MtM)[2,2])/2]
r=norm(rsc)
ϕ=mod(atan2((rsc/norm(rsc))...),2π)
(r,ϕ)
end

function mat2xyz(M)
[M[1,1]^2+M[2,1]^2-M[1,2]^2-M[2,2]^2
2*(M[1,1]*M[1,2]+M[2,1]*M[2,2])
2*abs(M[1,1]*M[2,2]-M[2,1]*M[1,2])]/2
end

function mat2xyz2(M)   # Alternative in terms of MtM
MtM=M'*M
[MtM[1,1]-MtM[2,2]
2*MtM[1,2]
2*√(det(MtM))]/2
end

#svd2...

function svd2tri(σ₁,σ₂,θ)
(1-(σ₁^2-σ₂^2)*cos(2θ+[ 2π; -2π; 0]/3))/3
end

function svd2hemi(σ₁,σ₂,θ)
λ=asin(2*σ₁*σ₂)
ϕ=2θ
λ,ϕ
end

function svd2disk(σ₁,σ₂,θ)
r=(σ₁^2-σ₂^2)/2
ϕ=2θ
r,ϕ
end

function svd2xyz(σ₁,σ₂,θ)
hemi2xyz(svd2hemi(σ₁,σ₂,θ)...)
end

function svd2mat(σ₁,σ₂,θ)
diagm([σ₁,σ₂])*[cos(θ) sin(θ);-sin(θ) cos(θ)]
end

#hemi2...

function hemi2svd(λ,ϕ)
σ₁=cos(λ/2)
σ₂=sin(λ/2)
θ=ϕ/2
σ₁,σ₂,θ
end

function hemi2tri(λ,ϕ)
r=cos(λ)/2
disk2tri(r,ϕ)
end

function hemi2disk(λ,ϕ)
r=cos(λ)/2
r,ϕ
end

function hemi2xyz(λ,ϕ)  # Standard Spherical Coordinates
[cos(λ)*cos(ϕ),cos(λ)*sin(ϕ),sin(λ)]/2
end

hemi2mat(λ,ϕ) = xyz2mat(hemi2xyz(λ,ϕ))
#tri2..

tri2svd(abc²) = disk2svd(tri2disk(abc²)...)

function tri2hemi(abc²)
cs = √(6)*Δ*abc²
λ = acos(norm(cs))
ϕ =mod( atan2((cs/norm(cs))...), 2π)
λ,ϕ
end

function tri2disk(abc²)
rsc=√(3/2)*(Δ*abc²)
r=norm(rsc)
rsc/=r
ϕ=mod(atan2(rsc...),2π)
r,ϕ
end

tri2xyz(abc²) = hemi2xyz(tri2hemi(abc²)...)

tri2mat(abc²) = svd2mat(tri2svd(abc²)...)

#disk2..

function disk2svd(r,ϕ)
σ₁= √(.5+r)
σ₂= √(.5-r)
θ=ϕ/2
σ₁,σ₂,θ
end

function disk2tri(r,ϕ)
abc² = (1-2r*cos(ϕ+[ 2π -2π 0]'/3))/3
end

function disk2hemi(r,ϕ)
λ=acos(2r)
λ,ϕ
end

disk2xyz(r,ϕ)=hemi2xyz(disk2hemi(r,ϕ)...)

disk2mat(r,ϕ)=hemi2mat(disk2hemi(r,ϕ)...)

# xyz2 ...

function xyz2hemi(xyz)
λ=asin(2*xyz[3])
ϕ=mod(atan2(xyz[2],xyz[1]),2π)
λ,ϕ
end

xyz2svd(xyz) =  hemi2svd(xyz2hemi(xyz)...)
xyz2tri(xyz) =  hemi2tri(xyz2hemi(xyz)...)
xyz2disk(xyz) = hemi2disk(xyz2hemi(xyz)...)
xyz2mat(xyz) =  svd2mat(xyz2svd(xyz)...)

Out[118]:
xyz2mat (generic function with 1 method)

In [120]:
M=randn(2,2); M/=normfro(M)
println(mat2xyz(M))
println(mat2xyz2(M))

.03794861637448002
-.10860468723512698
.4865849611587096

.03794861637448002
-.10860468723512698
.48658496115870953


In [106]:
# Mat test
M=randn(2,2); M/=normfro(M)
println(mat2mat(M) )
println(svd2mat(mat2svd(M)...))
println(tri2mat(mat2tri(M)))
println(hemi2mat(mat2hemi(M)...))
println(disk2mat(mat2disk(M)...))
println(xyz2mat(mat2xyz(M)))

.8135451570602026	.5444302126243655
-.11362577916422914	.1697916467762789

.8135451570602026	.5444302126243655
-.11362577916422914	.1697916467762789

.8135451570602027	.5444302126243649
-.1136257791642293	.16979164677627936

.8135451570602026	.5444302126243653
-.11362577916422914	.16979164677627895

.8135451570602028	.5444302126243649
-.11362577916422907	.16979164677627903

.8135451570602028	.5444302126243649
-.11362577916422903	.16979164677627895


In [107]:
# SVD TEST
sigma=sort(rand(2));sigma/=sum(sigma);
(σ₂,σ₁)=√(sigma)
θ = rand()*π
println("-----")

abc   =  svd2tri(σ₁,σ₂,θ)
(λ,ϕ) = svd2hemi(σ₁,σ₂,θ)
(r,ϕ) =  svd2disk(σ₁,σ₂,θ)
xyz = svd2xyz(σ₁,σ₂,θ)
M = svd2mat(σ₁,σ₂,θ)
println(σ₁," ",σ₂," ",θ)
println(tri2svd(abc))
println(hemi2svd(λ,ϕ))
println( disk2svd(r,ϕ))
println( xyz2svd(xyz))
println(mat2svd(M))

-----
0.8159551239172137 0.5781152443529268 1.7176328482161802
(0.8159551239172137,0.5781152443529269,1.7176328482161802)
(0.8159551239172138,0.5781152443529267,1.7176328482161802)
(0.8159551239172137,0.5781152443529268,1.7176328482161802)
(0.8159551239172138,0.5781152443529267,1.71763284821618)
(0.8159551239172137,0.5781152443529268,1.71763284821618)


In [108]:
# TRI TEST
abc=rand(3);
abc/=sum(abc)
sl=sort(sqrt(abc))
if sl[1]+sl[2]>sl[3]  # Throw away non triangle inequalities
(λ,ϕ) = tri2hemi(abc)
(σ₁,σ₂,θ)=tri2svd(abc)
(r,ϕ)=tri2disk(abc)
xyz=tri2xyz(abc)
M = tri2mat(abc)
[abc svd2tri(σ₁,σ₂,θ) hemi2tri(λ,ϕ) disk2tri(r,ϕ) xyz2tri(xyz) mat2tri(M)]
end

Out[108]:
3x6 Array{Float64,2}:
0.525615   0.525615   0.525615   0.525615   0.525615   0.525615
0.0895413  0.0895413  0.0895413  0.0895413  0.0895413  0.0895413
0.384843   0.384843   0.384843   0.384843   0.384843   0.384843

In [109]:
# HEMI TEST
λ = rand()*pi/2
ϕ = rand() *2*pi

(σ₁,σ₂,θ)= hemi2svd(λ,ϕ)
abc = hemi2tri(λ,ϕ)
(r,ϕ) =  hemi2disk(λ,ϕ)
xyz= hemi2xyz(λ,ϕ)
M = hemi2mat(λ,ϕ)

println(λ," ",ϕ)
println(svd2hemi(σ₁,σ₂,θ))
println(tri2hemi(abc))
println( disk2hemi(r,ϕ))
println(xyz2hemi(xyz))
println(mat2hemi(M))

1.5325386800127676 2.1889132842143413
(1.5325386800127698,2.1889132842143413)
(1.5325386800127676,2.188913284214344)
(1.5325386800127676,2.1889132842143413)
(1.532538680012767,2.1889132842143404)
(1.5325386800127612,2.1889132842143404)


In [110]:
# DISK TEST
r = rand()*1/2
ϕ = rand() *2*pi

(σ₁,σ₂,θ)= disk2svd(r,ϕ)
abc = disk2tri(r,ϕ)
(λ,ϕ) =  disk2hemi(r,ϕ)
xyz = disk2xyz(r,ϕ)
M = disk2mat(r,ϕ)
println(r," ",ϕ)
println(svd2disk(σ₁,σ₂,θ))
println(tri2disk(abc))
println(hemi2disk(λ,ϕ))
println(xyz2disk(xyz))
println(mat2disk(M))

0.3429393172667903 5.610928823821522
(0.3429393172667903,5.610928823821522)
(0.3429393172667902,5.610928823821522)
(0.3429393172667903,5.610928823821522)
(0.34293931726679033,5.610928823821522)
(0.3429393172667902,5.610928823821521)


In [111]:
#xyz test
xyz=randn(3);xyz[3]=abs(xyz[3]);xyz/=norm(xyz)*2
[xyz svd2xyz(xyz2svd(xyz)...) tri2xyz(xyz2tri(xyz)) hemi2xyz(xyz2hemi(xyz)...)  disk2xyz(xyz2disk(xyz)...)  mat2xyz(xyz2mat(xyz))]

Out[111]:
3x6 Array{Float64,2}:
-0.170281  -0.170281  -0.170281  -0.170281  -0.170281  -0.170281
0.031067   0.031067   0.031067   0.031067   0.031067   0.031067
0.469083   0.469083   0.469083   0.469083   0.469083   0.469083

In []: