Attrac-0.1.2: Attrac.hs
module Attrac where
import Array
import System.Random
type Value = Double
data Pt = Pt Value Value Value deriving Show
data Mat = Mat Value Value Value Value Value Value Value Value Value deriving Show
type Coeffs = Array Int Value
--cs = array (0,29) (zip [0..] (readCode whichCode))
-- where whichCode = codes!!37
-- 16 = tetrahedron
-- 29 is happy face
f cs (Pt x y z) = Pt
((cs!0) + (cs!1)*x + (cs!2)*x*x + (cs!3)*x*y + (cs!4)*x*z + (cs!5)*y + (cs!6)*y*y + (cs!7)*y*z + (cs!8)*z + (cs!9)*z*z)
((cs!10) + (cs!11)*x + (cs!12)*x*x + (cs!13)*x*y + (cs!14)*x*z + (cs!15)*y + (cs!16)*y*y + (cs!17)*y*z + (cs!18)*z + (cs!19)*z*z)
((cs!20) + (cs!21)*x + (cs!22)*x*x + (cs!23)*x*y + (cs!24)*x*z + (cs!25)*y + (cs!26)*y*y + (cs!27)*y*z + (cs!28)*z + (cs!29)*z*z)
fjac cs (Pt x y z) = Mat
((cs!1)+2*(cs!2)*x+(cs!3)*y+(cs!4)*z)
((cs!3)*x+(cs!5)+2*(cs!6)*y+(cs!7)*z)
((cs!4)*x+(cs!7)*y+(cs!8)+2*(cs!9)*z)
((cs!11)+2*(cs!12)*x+(cs!13)*y+(cs!14)*z)
((cs!13)*x+(cs!15)+2*(cs!16)*y+(cs!17)*z)
((cs!14)*x+(cs!17)*y+(cs!18)+2*(cs!19)*z)
((cs!21)+2*(cs!22)*x+(cs!23)*y+(cs!24)*z)
((cs!23)*x+(cs!25)+2*(cs!26)*y+(cs!27)*z)
((cs!24)*x+(cs!27)*y+(cs!28)+2*(cs!29)*z)
normCols (Mat a11 a12 a13 a21 a22 a23 a31 a32 a33) = Mat (a11*a) (a12*b) (a13*c) (a21*a) (a22*b) (a23*c) (a31*a) (a32*b) (a33*c)
where a = recip (sqrt (a11*a11 + a21*a21 + a31*a31))
b = recip (sqrt (a12*a12 + a22*a22 + a32*a32))
c = recip (sqrt (a13*a13 + a23*a23 + a33*a33))
fm cs (p,m) = (f cs p,normCols $ (fjac cs p) <*> m)
(<*>) :: Mat -> Mat -> Mat
(<*>) (Mat a11 a12 a13 a21 a22 a23 a31 a32 a33) (Mat b11 b12 b13 b21 b22 b23 b31 b32 b33)
= Mat (a11*b11 + a12*b21 + a13*b31) (a11*b12 + a12*b22 + a13*b32) (a11*b13 + a12*b23 + a13*b33)
(a21*b11 + a22*b21 + a23*b31) (a21*b12 + a22*b22 + a23*b32) (a21*b13 + a22*b23 + a23*b33)
(a31*b11 + a32*b21 + a33*b31) (a31*b12 + a32*b22 + a33*b32) (a31*b13 + a32*b23 + a33*b33)
zeroPt = Pt 0 0 0
idMat = Mat 1 0 0 0 1 0 0 0 1
almostZeroPt = peturb peturbationAmount zeroPt
(<.>) :: Pt -> Pt -> Value
(<.>) (Pt ax ay az) (Pt bx by bz) = ax*bx + ay*by + az*bz
(<+>) :: Pt -> Pt -> Pt
(<+>) (Pt ax ay az) (Pt bx by bz) = Pt (ax+bx) (ay+by) (az+bz)
(<->) :: Pt -> Pt -> Pt
(<->) (Pt ax ay az) (Pt bx by bz) = Pt (ax-bx) (ay-by) (az-bz)
(<#>) :: Value -> Pt -> Pt
(<#>) a (Pt x y z) = Pt (a*x) (a*y) (a*z)
normalise :: Pt -> Pt
normalise p = (recip (sqrt (p <.> p))) <#> p
norm a = a <.> a
data PMat = PMat
Value Value Value Value
Value Value Value Value
Value Value Value Value
Value Value Value Value deriving Show
idPMat = PMat
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
(<*%>) :: PMat -> PMat -> PMat
(<*%>) (PMat a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44)
(PMat b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44)
= PMat (a11*b11 + a12*b21 + a13*b31 + a14*b41)
(a11*b12 + a12*b22 + a13*b32 + a14*b42)
(a11*b13 + a12*b23 + a13*b33 + a14*b43)
(a11*b14 + a12*b24 + a13*b34 + a14*b44)
(a21*b11 + a22*b21 + a23*b31 + a24*b41)
(a21*b12 + a22*b22 + a23*b32 + a24*b42)
(a21*b13 + a22*b23 + a23*b33 + a24*b43)
(a21*b14 + a22*b24 + a23*b34 + a24*b44)
(a31*b11 + a32*b21 + a33*b31 + a34*b41)
(a31*b12 + a32*b22 + a33*b32 + a34*b42)
(a31*b13 + a32*b23 + a33*b33 + a34*b43)
(a31*b14 + a32*b24 + a33*b34 + a34*b44)
(a41*b11 + a42*b21 + a43*b31 + a44*b41)
(a41*b12 + a42*b22 + a43*b32 + a44*b42)
(a41*b13 + a42*b23 + a43*b33 + a44*b43)
(a41*b14 + a42*b24 + a43*b34 + a44*b44)
(<*#>) :: PMat -> Pt -> Pt
(<*#>) (PMat a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44) (Pt x y z) =
Pt (((x*a11)+(y*a12)+(z*a13)+a14) / s) (((x*a21)+(y*a22)+(z*a23)+a24) / s) (((x*a31)+(y*a32)+(z*a33)+a34) / s)
where s = ((x*a41)+(y*a42)+(z*a43)+a44)
-----
-- the amount to peturb orbits in the Lyapunov Exponent estimation
peturbationAmount = 1e-12
-- any point outside this radius from the origin is considered part
-- of an unbounded orbit
unboundedRadius = 4.5
lyap :: Coeffs -> Int -> Pt -> Pt -> Maybe Value
lyap cs t a b = (lyapHelper cs t a b) >>= (return . (/ (fromIntegral t)))
where
lyapHelper :: Coeffs -> Int -> Pt -> Pt -> Maybe Value
lyapHelper _ 0 _ _ = Just 0
lyapHelper cs t a b | isBadNum na || na > unboundedRadius^2 = Nothing
| otherwise = (lyapHelper cs (t-1) fa fb2) >>= (return . (s +))
where
na = norm a
fa = f cs a
fb = f cs b
dx2 = norm (fa <-> fb)
n = peturbationAmount <#> (normalise (fb <-> fa))
fb2 = fa <+> n
s = 0.5 * (logBase 2 (dx2 / (peturbationAmount^2)))
isBadNum x = isNaN x || isInfinite x
lyapunovAccuracy = 2048
maxLyapunovExponent :: Coeffs -> Maybe Value
maxLyapunovExponent cs = lyap cs lyapunovAccuracy p q
where
p = (iterate (f cs) almostZeroPt) !! 16
q = peturb peturbationAmount p
numCoeffs = 30
grabRandoms :: Int -> (Int,Int) -> StdGen -> ([Int],StdGen)
grabRandoms 0 _ g = ([],g)
grabRandoms n (x,y) g = let (r,g') = randomR (x,y) g in let (rs,g'') = grabRandoms (n-1) (x,y) g' in (r:rs,g'')
minValidLyapunovExponent = 1e-9
randomiseCoeffs :: StdGen -> (Coeffs,Value)
randomiseCoeffs g = case maxLyapunovExponent cs of {Just x | x > minValidLyapunovExponent -> (cs,x); _ -> randomiseCoeffs g'}
where (as,g') = grabRandoms numCoeffs (-13,13) g
cs = listArray (0,numCoeffs-1) (map (\a -> 0.1 * (fromIntegral a)) as)
peturb :: Value -> Pt -> Pt
peturb e (Pt x y z) = Pt (x+e) y z
(<:*>) :: Pt -> Pt -> Pt
(<:*>) (Pt ax ay az) (Pt bx by bz) = Pt (ay*bz - az*by) (bx*az - ax*bz) (ax*by - ay*bx)