{-
Haskell Dynamics Engine, Copyright (C) 2007 Ruben Henner Zilibowitz
All rights reserved. Email: rubenz@cse.unsw.edu.au Web: www.cse.unsw.edu.au/~rubenz
This is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License that is include with this
package in the file LICENSE-GPL.TXT.
-}
module SeparatingLine(lineOfMaxVariance) where
import Matrix3x3
import Roots(roots)
import Data.Complex
{-
This wierd code is supposed to solve the following problem:
Given a set of points, find the line which maximises the
variance of the values of the points projected onto that line.
In general this line seems quite similar to the line of best
fit obtained by orthogonal regression although it is defined
differently and is not usually the same line.
The line can be found with linear time in the number of given points.
-}
lineOfMaxVariance :: [Matrix3x3.Vector Double] -> (Matrix3x3.Vector Double,Matrix3x3.Vector Double)
lineOfMaxVariance vs = (av_vs,calcSolution u (minimum lambdas))
where
n = fromIntegral (length vs)
av_vs = (sum vs) ./. n
ws = map (\x -> x - av_vs) vs
u@(u1,u2,u3,u4,u5,u6) = foldl1 addSixtuples (map vecCoeffs ws)
cubicCoeffs = map realToComplex [-2*(u1*u6^2-u4*u5*u6+u2*u5^2+u3*u4^2-4*u1*u2*u3),-2*(u6^2+u5^2+u4^2-4*(u2*u3+u1*u3+u1*u2)),8*(u3+u2+u1),8]
lambdas = [x | (x :+ y) <- (roots 1e-6 100 cubicCoeffs), abs y < 1e-12, not$isNaN x]
addSixtuples (a1,a2,a3,a4,a5,a6) (b1,b2,b3,b4,b5,b6) = (a1+b1,a2+b2,a3+b3,a4+b4,a5+b5,a6+b6)
vecCoeffs (Vec x y z) = (x*x,y*y,z*z,2*x*y,2*x*z,2*y*z)
realToComplex x = x :+ 0
calcSolution (u1,u2,u3,u4,u5,u6) lambda = Vec ((u4*u6-2*u2l*u5)/denom) (-(2*u1l*u6-u4*u5)/denom) (-(u4^2-4*u1l*u2l)/denom)
where
u1l = u1 + lambda
u2l = u2 + lambda
u3l = u3 + lambda
denom = sqrt ((u4^2+4*u1l^2)*u6^2+(-4*u2l-4*u1l)*u4*u5*u6+(u4^2+4*u2l^2)*u5^2+u4^4-8*u1l*u2l*u4^2+16*u1l^2*u2l^2)