31:Tc10,md`# Rotor Estimation and AD in GA. In this notebook we will use [ganja.js](https://github.com/enkimute/ganja.js) to estimate the transformation between two point clouds. To perform this task we will use a Clifford Algebra over the Dual Vectors. Ganja.js normally only supports Clifford Algebras over the reals, but with a little trickery, we can easily create this exotic algebra. To do automatic differentiation, specifically to calculate partial derivatives with just one function evaluation, one needs dual vectors, consisting out of a scalar and ${tex`n`} dual coefficients. For ${tex`n=3`} the Cayley table of such an algebra is :

1	e₀₁	e₀₂	e₀₃
e₀₁	0	0	0
e₀₂	0	0	0
e₀₃	0	0	0

Using the Custom Cayley table feature of [ganja.js](https://github.com/enkimute/ganja.js), we can easily create such an Algebra. Additionally we also create the Clifford Algebra ${tex`\mathbb R_{3,0}`}. We then use the ganja.js translator of the Dual Vectors class to translate the operators of ${tex`\mathbb R_{3,0}`}. While the original functions of ${tex`\mathbb R_{3,0}`} support only scalars, in these new versions each of the coefficients of a multivector can itself be a Dual Vector. For example, for ${tex`\mathbb R_{2,0}`} the body of the function that implements the geometric product normally looks like this : \`\`\`javascript res[0]=b[0]*a[0]+b[1]*a[1]+b[2]*a[2]-b[3]*a[3]; res[1]=b[1]*a[0]+b[0]*a[1]-b[3]*a[2]+b[2]*a[3]; res[2]=b[2]*a[0]+b[3]*a[1]+b[0]*a[2]-b[1]*a[3]; res[3]=b[3]*a[0]+b[2]*a[1]-b[1]*a[2]+b[0]*a[3]; \`\`\` Where **res**, **a** and **b** are javascript Arrays. After translating it using the **inline** function of the Dual Vector algebra, it will look like this: \`\`\`javascript res[0]=Dual.Sub(Dual.Add(Dual.Add(Dual.Mul(b[0],a[0]),Dual.Mul(b[1],a[1])),Dual.Mul(b[2],a[2])),Dual.Mul(b[3],a[3])); res[1]=Dual.Add(Dual.Sub(Dual.Add(Dual.Mul(b[1],a[0]),Dual.Mul(b[0],a[1])),Dual.Mul(b[3],a[2])),Dual.Mul(b[2],a[3])); res[2]=Dual.Sub(Dual.Add(Dual.Add(Dual.Mul(b[2],a[0]),Dual.Mul(b[3],a[1])),Dual.Mul(b[0],a[2])),Dual.Mul(b[1],a[3])); res[3]=Dual.Add(Dual.Sub(Dual.Add(Dual.Mul(b[3],a[0]),Dual.Mul(b[2],a[1])),Dual.Mul(b[1],a[2])),Dual.Mul(b[0],a[3])); \`\`\` It still operates on an array, but all coefficients can now be Dual Vectors. These functions can then be used to implement the automatic differentiation for rotation and cost functions needed to do the estimation. We use a simple Gradient Descent - much more optimal algorithms exist (and can be implemented without Linear Algebra). In the example below, we parametrize the rotor manifold using a pure bivector ${tex`R`}, and minimize the simple cost function ${tex`\frac 1 n \sum\limits_{i=1}\limits^n | e^{- \frac R 2}\mathbf xe^{\frac R 2} - y |^2`}. We will evaluate this function using our Algebra over Dual Vectors which will result in automatically producing the needed derivatives. `32:T1080,output = { //** Generate an Algebra with x dual components. **// var Dualgebra = (x)=>Algebra({ basis : [...Array(x+1)].map((a,i)=>i?'e0'+i:'1'), metric : [1,...Array(x)], Cayley : [...Array(x+1)].map((a,i)=>[...Array(x+1)].map((y,j)=>i*j==0?((i+j)?'e0'+(i+j):'1'):'0')) }); //** We'll use one with three dual components, and a clifford algebra over R3 **// var D3 = Dualgebra(3), R3 = Algebra(3); //** helper .. easily convert scalars or arrays of scalars to dual vectors. **// var d3 = (a,b,c,d)=>(a instanceof Array)?a.map(x=>d3(x)):(a instanceof D3)?a:new D3([a||0,b||0,c||0,d||0]); //** We upgrade the operators of R3 we need to work with D3 coefficients. This is where the functions of R3 **// //** get translated to work with Dual Vectors for each of their coefficients. **// var translate = (x)=>{ var F=D3.inline(R3.prototype[x]); return (a,b)=>F.call(a,b); }, [GP,ADD,SUB,SCALE] = ["Mul","Add","Sub","Scale"].map(x=>translate(x)); //** We implement our own exp function using these operators and Taylor **// var EXP = (X)=>{ var R=d3([1,...Array(X.length-1)]), M, y=1, i=1; while (i<25) { R=ADD(R, SCALE(M||X, 1/y)); M=GP(M||X,X); y=y*(++i); } return R; }; //** random float between -x and x, vector length, rotate, cost function (avg of squared diffs) **// var r = (x)=>(Math.random()*2-1)*(x||1), VLEN = R=>R[1].Mul(R[1]).Add(R[2].Mul(R[2])).Add(R[3].Mul(R[3])), ROT = (R,x)=>GP(GP(EXP(SCALE(R,-0.5)),x),EXP(SCALE(R,0.5))), COST = (G,Pa,Pb)=>Pa.reduce((C,Pa,x)=>C.Add(VLEN(SUB(ROT(G,Pa),Pb[x]))),d3()).Scale(1/Pa.length); //** Now we generate some random points, a rotation bivector, the rotated noisy points and our guess+dBV **// //** For our guess, Each of the coefficients we wish to track (e12,e13,e23) gets a different dual unit added in. **// var points_a = [...Array(30)].map(x=>[0,r(),r(),r(),0,0,0,0]), bv = d3([0,0,0,0,r(),r(),r(),0]), points_b = points_a.map(x=>ADD(ROT(bv,x),[0,r(.01),r(.01),r(.01),0,0,0,0]).map(x=>x[0])), guess = d3([0,0,0,0,d3(0,1),d3(0,0,1),d3(0,0,0,1),0]); //** Now we implement Gradient Descent - step in the direction of negative cost **// var start = performance.now(); for (var alpha=1,i=0; i<100; i++) { var c = COST(guess,points_a,points_b); if (c[0]<5E-4) break; guess = ADD(guess,[0,0,0,0,-c[1]*alpha,-c[2]*alpha,-c[3]*alpha,0]) } var stop = performance.now(); //** We are done.. code below is graphing results.. **// var graph = Algebra(3,0,1).inline((pa,pb,i,guess)=>{ // convert the arrays to 3D PGA points. var conv = x=>x.map(x=>!(x[1][0]*1e1+(x[2][0])*1e2+x[3][0]*1e3+1e0)); var trail = (bv,x,n=12)=>[...Array(n).keys()].map(i=>Math.E**(bv*(i/(n-1)))>>>x).map((x,i,a)=>[x,a[i-1]||x]).slice(1); // initial camera and items to render. pa = conv(pa); var trails = []; pa.forEach(pa=>trails.push(...trail(guess*[-.5e12,-.5e13,-.5e23], pa))); var camera = Math.E**1e23,f=0, items=[...pa,0x88FF0000,...conv(pb),0x880000FF]; // graph in gl, animate and wiggle the camera return this.graph(()=>{ camera.set(Math.E**(2e03)*Math.E**(Math.sin(performance.now()/5000)*1e13)); var anim = Math.E**((Math.sin(performance.now()/234)*0.5+.5)*( guess*[-.5e12,-.5e13,-.5e23] )); // nice .. in PGA, we can just render the bivector and get the expected line.. return [0xAAAAFF,...trails,0x0000FF,guess*[1e12,1e13,1e23],0,...items,...anim>>>pa] ; },{grid:1,camera,animate:1,gl:1,width:'100%',height:'400px'}); })(d3(points_a),d3(points_b),i,guess.map(x=>x[0]).slice(4,7)); //** Also output some statistics. **// var text = md` * Nr of iterations : ${i} * Nr of points : ${points_a.length} * Time : ${(stop-start)/1000} * Final Error : ${c[0]} * Initial guess : 0.0,0.0,0.0 * Final estimate : ${guess[4][0].toFixed(4)},${guess[5][0].toFixed(4)},${guess[6][0].toFixed(4)} * Ground truth : ${bv[4][0].toFixed(4)},${bv[5][0].toFixed(4)},${bv[6][0].toFixed(4)} * GT Error (noise) : ${COST(bv,points_a,points_b)[0]}`; return {graph,text} }33:[] 2f:{"title":"Rotor Estimation and AD in GA.","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"md`**in the graph : **\n\n* black : points in frame A\n* red : points in frame B + noise\n* blue : estimate`","pinned":false,"hidden":false,"mode":"ojs"},{"id":43,"value":"output.graph","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"output.text","pinned":false,"hidden":false,"mode":"ojs"},{"id":49,"value":"md`**Press the arrow in the grey area below, on the top-right, to re-run.**`","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":34,"value":"Algebra = require('ganja.js')","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}