{
"metadata": {
"language": "Julia",
"name": "",
"signature": "sha256:aca186dc3fce49eea63d2e42629fe67646b76fe248fc29c235a04733bd4e38b2"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Code to accompany Edelman and Strang: \n",
" \n",
"Random Triangle Theory with Geometry and Applications
\n",
"\n",
"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:\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our paper revisits this work. The code below specs out the conversion between the representations."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\u221a(x)=sqrt(x)\n",
"\u0394=[1/\u221a(2) -1/\u221a(2) 0;1/\u221a(6) 1/\u221a(6) -2/\u221a(6)]\n",
"\n",
"#Representations: mat (2x2 matrix M)\n",
"# svd \u03c3\u2081,\u03c3\u2082,\u03b8\n",
"# tri (3 vector of squared edge lengths)\n",
"# hemi \u03bb,\u03d5 (latitude, longitude)\n",
"# disk r,\u03d5 (polar coordinates)\n",
"# xyz (3 vector of cartesian coordinates on the hemisphere of radius 1/2)\n",
"\n",
"# CONVERSIONS:\n",
"\n",
"#mat2...\n",
"\n",
"function mat2svd(M)\n",
" (U,(\u03c3\u2081,\u03c3\u2082),V)=svd(M)\n",
" # Pick the right \u03b8 in [0,\u03c0) \n",
" \u03b8=acos(V[1,1]*sign(det(M))*sign(V[1,2])) \n",
" (\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
"end\n",
"\n",
"mat2tri(M) = diag( (M*\u0394)'*(M*\u0394)) \n",
"mat2hemi(M) = svd2hemi( mat2svd(M)...)\n",
"mat2mat(M) = svd2mat(mat2svd(M)...) # Canonical Choice (throw away U)\n",
"\n",
"function mat2disk(M)\n",
" M/=vecnorm(M)\n",
" MtM=M'*M\n",
" rsc = [ MtM[1,2] ; ((MtM)[1,1]-(MtM)[2,2])/2]\n",
" r=norm(rsc)\n",
" \u03d5=mod(atan2((rsc/norm(rsc))...),2\u03c0)\n",
" (r,\u03d5)\n",
"end\n",
"\n",
"function mat2xyz(M)\n",
" [M[1,1]^2+M[2,1]^2-M[1,2]^2-M[2,2]^2\n",
" 2*(M[1,1]*M[1,2]+M[2,1]*M[2,2])\n",
" 2*abs(M[1,1]*M[2,2]-M[2,1]*M[1,2])]/2\n",
"end\n",
"\n",
"function mat2xyz2(M) # Alternative in terms of MtM\n",
" MtM=M'*M\n",
" [MtM[1,1]-MtM[2,2]\n",
" 2*MtM[1,2] \n",
" 2*\u221a(det(MtM))]/2\n",
"end\n",
"\n",
"#svd2...\n",
"\n",
"function svd2tri(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
" (1-(\u03c3\u2081^2-\u03c3\u2082^2)*cos(2\u03b8+[ 2\u03c0; -2\u03c0; 0]/3))/3\n",
"end\n",
"\n",
"function svd2hemi(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
" \u03bb=asin(2*\u03c3\u2081*\u03c3\u2082)\n",
" \u03d5=2\u03b8\n",
" \u03bb,\u03d5 \n",
"end\n",
"\n",
"function svd2disk(\u03c3\u2081,\u03c3\u2082,\u03b8) \n",
" r=(\u03c3\u2081^2-\u03c3\u2082^2)/2\n",
" \u03d5=2\u03b8\n",
" r,\u03d5\n",
"end\n",
"\n",
"function svd2xyz(\u03c3\u2081,\u03c3\u2082,\u03b8) \n",
" hemi2xyz(svd2hemi(\u03c3\u2081,\u03c3\u2082,\u03b8)...)\n",
"end\n",
"\n",
"function svd2mat(\u03c3\u2081,\u03c3\u2082,\u03b8) \n",
" diagm([\u03c3\u2081,\u03c3\u2082])*[cos(\u03b8) sin(\u03b8);-sin(\u03b8) cos(\u03b8)]\n",
"end\n",
"\n",
"#hemi2...\n",
"\n",
"function hemi2svd(\u03bb,\u03d5)\n",
" \u03c3\u2081=cos(\u03bb/2)\n",
" \u03c3\u2082=sin(\u03bb/2)\n",
" \u03b8=\u03d5/2\n",
" \u03c3\u2081,\u03c3\u2082,\u03b8\n",
"end\n",
"\n",
"function hemi2tri(\u03bb,\u03d5)\n",
" r=cos(\u03bb)/2\n",
" disk2tri(r,\u03d5)\n",
"end\n",
"\n",
"function hemi2disk(\u03bb,\u03d5)\n",
" r=cos(\u03bb)/2\n",
" r,\u03d5\n",
"end\n",
"\n",
"function hemi2xyz(\u03bb,\u03d5) # Standard Spherical Coordinates\n",
" [cos(\u03bb)*cos(\u03d5),cos(\u03bb)*sin(\u03d5),sin(\u03bb)]/2\n",
"end\n",
"\n",
"hemi2mat(\u03bb,\u03d5) = xyz2mat(hemi2xyz(\u03bb,\u03d5))\n",
"#tri2..\n",
"\n",
"tri2svd(abc\u00b2) = disk2svd(tri2disk(abc\u00b2)...)\n",
"\n",
"function tri2hemi(abc\u00b2)\n",
" cs = \u221a(6)*\u0394*abc\u00b2\n",
" \u03bb = acos(norm(cs))\n",
" \u03d5 =mod( atan2((cs/norm(cs))...), 2\u03c0)\n",
" \u03bb,\u03d5\n",
"end\n",
"\n",
"function tri2disk(abc\u00b2)\n",
" rsc=\u221a(3/2)*(\u0394*abc\u00b2)\n",
" r=norm(rsc)\n",
" rsc/=r\n",
" \u03d5=mod(atan2(rsc...),2\u03c0)\n",
" r,\u03d5\n",
"end\n",
"\n",
"tri2xyz(abc\u00b2) = hemi2xyz(tri2hemi(abc\u00b2)...)\n",
"\n",
"tri2mat(abc\u00b2) = svd2mat(tri2svd(abc\u00b2)...)\n",
"\n",
"#disk2..\n",
"\n",
"function disk2svd(r,\u03d5)\n",
" \u03c3\u2081= \u221a(.5+r)\n",
" \u03c3\u2082= \u221a(.5-r)\n",
" \u03b8=\u03d5/2\n",
" \u03c3\u2081,\u03c3\u2082,\u03b8\n",
"end\n",
"\n",
"function disk2tri(r,\u03d5)\n",
" abc\u00b2 = (1-2r*cos(\u03d5+[ 2\u03c0 -2\u03c0 0]'/3))/3\n",
"end\n",
"\n",
"function disk2hemi(r,\u03d5)\n",
" \u03bb=acos(2r)\n",
" \u03bb,\u03d5\n",
"end\n",
"\n",
"disk2xyz(r,\u03d5)=hemi2xyz(disk2hemi(r,\u03d5)...)\n",
"\n",
"disk2mat(r,\u03d5)=hemi2mat(disk2hemi(r,\u03d5)...)\n",
"\n",
"# xyz2 ...\n",
"\n",
"function xyz2hemi(xyz)\n",
" \u03bb=asin(2*xyz[3])\n",
" \u03d5=mod(atan2(xyz[2],xyz[1]),2\u03c0)\n",
" \u03bb,\u03d5\n",
"end\n",
"\n",
"xyz2svd(xyz) = hemi2svd(xyz2hemi(xyz)...)\n",
"xyz2tri(xyz) = hemi2tri(xyz2hemi(xyz)...)\n",
"xyz2disk(xyz) = hemi2disk(xyz2hemi(xyz)...)\n",
"xyz2mat(xyz) = svd2mat(xyz2svd(xyz)...)\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 14,
"text": [
"xyz2mat (generic function with 1 method)"
]
}
],
"prompt_number": 14
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M=randn(2,2); M/=vecnorm(M)\n",
"println(mat2xyz(M))\n",
"println(mat2xyz2(M))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[0."
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"14340274318785814,0.4762002853805391,0.05166179874617707]\n",
"[0.14340274318785812,0.4762002853805391,0.05166179874617716]\n"
]
}
],
"prompt_number": 15
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Mat test\n",
"M=randn(2,2); M/=vecnorm(M)\n",
"println(mat2mat(M) )\n",
"println(svd2mat(mat2svd(M)...))\n",
"println(tri2mat(mat2tri(M)))\n",
"println(hemi2mat(mat2hemi(M)...))\n",
"println(disk2mat(mat2disk(M)...))\n",
"println(xyz2mat(mat2xyz(M)))\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[0."
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"9124707490665185 0.4062314276926464\n",
" -0.019813121282384303 0.04450392654396357]\n",
"[0.9124707490665185 0.4062314276926464\n",
" -0.019813121282384303 0.04450392654396357]\n",
"[0.9124707490665184 0.4062314276926463\n",
" -0.019813121282384074 0.04450392654396306]\n",
"[0.9124707490665184 0.4062314276926463\n",
" -0.0198131212823843 0.04450392654396357]\n",
"[0.9124707490665184 0.4062314276926463\n",
" -0.019813121282384306 0.04450392654396358]\n",
"[0.9124707490665184 0.4062314276926463\n",
" -0.019813121282384282 0.044503926543963526]\n"
]
}
],
"prompt_number": 16
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# SVD TEST\n",
" sigma=sort(rand(2));sigma/=sum(sigma);\n",
"(\u03c3\u2082,\u03c3\u2081)=\u221a(sigma)\n",
"\u03b8 = rand()*\u03c0\n",
" println(\"-----\")\n",
" \n",
"abc = svd2tri(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
"(\u03bb,\u03d5) = svd2hemi(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
"(r,\u03d5) = svd2disk(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
"xyz = svd2xyz(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
"M = svd2mat(\u03c3\u2081,\u03c3\u2082,\u03b8)\n",
" println(\u03c3\u2081,\" \",\u03c3\u2082,\" \",\u03b8)\n",
" println(tri2svd(abc))\n",
" println(hemi2svd(\u03bb,\u03d5))\n",
" println( disk2svd(r,\u03d5))\n",
" println( xyz2svd(xyz))\n",
" println(mat2svd(M))\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"-----\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.7736045131877425 0.6336687282607182 0.6720132434941296\n",
"(0.7736045131877425,0.6336687282607182,0.6720132434941295)\n",
"(0.7736045131877426,0.6336687282607182,0.6720132434941296)\n",
"(0.7736045131877425,0.6336687282607182,0.6720132434941296)\n",
"(0.7736045131877426,0.6336687282607182,0.6720132434941295)\n",
"(0.7736045131877425,0.6336687282607182,0.6720132434941296)\n"
]
}
],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# TRI TEST\n",
" abc=rand(3); \n",
" abc/=sum(abc) \n",
" sl=sort(sqrt(abc))\n",
" if sl[1]+sl[2]>sl[3] # Throw away non triangle inequalities \n",
" (\u03bb,\u03d5) = tri2hemi(abc) \n",
" (\u03c3\u2081,\u03c3\u2082,\u03b8)=tri2svd(abc) \n",
" (r,\u03d5)=tri2disk(abc)\n",
" xyz=tri2xyz(abc)\n",
" M = tri2mat(abc)\n",
" [abc svd2tri(\u03c3\u2081,\u03c3\u2082,\u03b8) hemi2tri(\u03bb,\u03d5) disk2tri(r,\u03d5) xyz2tri(xyz) mat2tri(M)]\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 18,
"text": [
"3x6 Array{Float64,2}:\n",
" 0.0494396 0.0494396 0.0494396 0.0494396 0.0494396 0.0494396\n",
" 0.394278 0.394278 0.394278 0.394278 0.394278 0.394278 \n",
" 0.556282 0.556282 0.556282 0.556282 0.556282 0.556282 "
]
}
],
"prompt_number": 18
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# HEMI TEST \n",
" \u03bb = rand()*pi/2\n",
" \u03d5 = rand() *2*pi\n",
" \n",
" (\u03c3\u2081,\u03c3\u2082,\u03b8)= hemi2svd(\u03bb,\u03d5)\n",
" abc = hemi2tri(\u03bb,\u03d5)\n",
" (r,\u03d5) = hemi2disk(\u03bb,\u03d5)\n",
"xyz= hemi2xyz(\u03bb,\u03d5)\n",
"M = hemi2mat(\u03bb,\u03d5)\n",
"\n",
" println(\u03bb,\" \",\u03d5)\n",
" println(svd2hemi(\u03c3\u2081,\u03c3\u2082,\u03b8))\n",
" println(tri2hemi(abc))\n",
" println( disk2hemi(r,\u03d5))\n",
" println(xyz2hemi(xyz))\n",
"println(mat2hemi(M))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.7918973334140129"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" 2.221245164482709\n",
"(0.7918973334140129,2.221245164482709)\n",
"(0.7918973334140132,2.22124516448271)\n",
"(0.7918973334140129,2.221245164482709)\n",
"(0.7918973334140129,2.22124516448271)\n",
"(0.7918973334140129,2.22124516448271)\n"
]
}
],
"prompt_number": 19
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# DISK TEST\n",
" r = rand()*1/2\n",
" \u03d5 = rand() *2*pi\n",
" \n",
" (\u03c3\u2081,\u03c3\u2082,\u03b8)= disk2svd(r,\u03d5)\n",
" abc = disk2tri(r,\u03d5)\n",
" (\u03bb,\u03d5) = disk2hemi(r,\u03d5) \n",
"xyz = disk2xyz(r,\u03d5) \n",
"M = disk2mat(r,\u03d5) \n",
" println(r,\" \",\u03d5)\n",
" println(svd2disk(\u03c3\u2081,\u03c3\u2082,\u03b8))\n",
" println(tri2disk(abc))\n",
" println(hemi2disk(\u03bb,\u03d5))\n",
" println(xyz2disk(xyz))\n",
" println(mat2disk(M))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.2506514772995837"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" 4.348999401148094\n",
"(0.25065147729958376,4.348999401148094)\n",
"(0.2506514772995837,4.348999401148094)\n",
"(0.2506514772995837,4.348999401148094)\n",
"(0.2506514772995837,4.3489994011480935)\n",
"(0.2506514772995837,4.3489994011480935)\n"
]
}
],
"prompt_number": 20
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#xyz test\n",
"xyz=randn(3);xyz[3]=abs(xyz[3]);xyz/=norm(xyz)*2\n",
"[xyz svd2xyz(xyz2svd(xyz)...) tri2xyz(xyz2tri(xyz)) hemi2xyz(xyz2hemi(xyz)...) disk2xyz(xyz2disk(xyz)...) mat2xyz(xyz2mat(xyz))]\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 21,
"text": [
"3x6 Array{Float64,2}:\n",
" 0.407743 0.407743 0.407743 0.407743 0.407743 0.407743\n",
" 0.140406 0.140406 0.140406 0.140406 0.140406 0.140406\n",
" 0.253045 0.253045 0.253045 0.253045 0.253045 0.253045"
]
}
],
"prompt_number": 21
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\u0394"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 69,
"text": [
"2x3 Array{Float64,2}:\n",
" 0.707107 -0.707107 0.0 \n",
" 0.408248 0.408248 -0.816497"
]
}
],
"prompt_number": 69
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# An Example Triangle\n",
"T=[-2 1 1;-1 -1 2];\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 88
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M=T*\u0394'; M.^2"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 100,
"text": [
"2x2 Array{Float64,2}:\n",
" 4.5 1.5\n",
" 0.0 6.0"
]
}
],
"prompt_number": 100
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Translate the 1st vertex to the origin does not change M\n",
"T2=copy(T)\n",
"T2[1,:]+=2\n",
"T2[2,:]+=1\n",
"println(T2)\n",
"T2*\u0394'"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[0"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" 3 3\n",
" 0 0 3]\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 40,
"text": [
"2x2 Array{Float64,2}:\n",
" -2.12132 -1.22474\n",
" 0.0 -2.44949"
]
}
],
"prompt_number": 40
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Translate the 1st vertex to the origin does not change M\n",
"T3=float64(T)\n",
"T3[1,:].+=randn()\n",
"T3[2,:].+=randn()\n",
"\n",
"T3*\u0394'"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 64,
"text": [
"2x2 Array{Float64,2}:\n",
" -2.12132 -1.22474\n",
" 0.0 -2.44949"
]
}
],
"prompt_number": 64
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\u0394'*\u0394"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 67,
"text": [
"3x3 Array{Float64,2}:\n",
" 0.666667 -0.333333 -0.333333\n",
" -0.333333 0.666667 -0.333333\n",
" -0.333333 -0.333333 0.666667"
]
}
],
"prompt_number": 67
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Let's take the edge view\n",
"T\u2091=T*[0 -1 1;1 0 -1 ;-1 1 0]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 101,
"text": [
"2x3 Array{Int64,2}:\n",
" 0 3 -3\n",
" -3 3 0"
]
}
],
"prompt_number": 101
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M=T\u2091*\u0394'"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 103,
"text": [
"2x2 Array{Float64,2}:\n",
" -2.12132 3.67423\n",
" -4.24264 0.0 "
]
}
],
"prompt_number": 103
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M.^2"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 104,
"text": [
"2x2 Array{Float64,2}:\n",
" 4.5 13.5\n",
" 18.0 0.0"
]
}
],
"prompt_number": 104
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\u0394'*\u0394"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 81,
"text": [
"3x3 Array{Float64,2}:\n",
" 0.666667 -0.333333 -0.333333\n",
" -0.333333 0.666667 -0.333333\n",
" -0.333333 -0.333333 0.666667"
]
}
],
"prompt_number": 81
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"trace(T'*T)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 82,
"text": [
"1.0000000000000002"
]
}
],
"prompt_number": 82
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\u0394*\u0394'"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 83,
"text": [
"2x2 Array{Float64,2}:\n",
" 1.0 0.0\n",
" 0.0 1.0"
]
}
],
"prompt_number": 83
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}