{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs #-}
module Math.Regression.Simple (
linear,
quadratic,
quadraticAndLinear,
Add (..),
Eye (..),
Mult (..),
Det (..),
Inv (..),
zerosLin,
zerosQuad,
optimaQuad,
V2 (..),
M22 (..),
V3 (..),
M33 (..),
Foldable' (..),
IsDoublePair (..),
) where
import Data.Complex (Complex (..))
import qualified Data.List as L
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
class Add a where
zero :: a
add :: a -> a -> a
class Eye a where
eye :: a
class Eye a => Mult a b c | a b -> c where
mult :: a -> b -> c
class Eye a => Det a where
det :: a -> Double
class Det a => Inv a where
inv :: a -> a
infixl 6 `add`
infixl 7 `mult`
instance Eye Double where
eye = 1
instance Add Double where
zero = 0
add = (+)
instance Det Double where
det = id
instance Inv Double where
inv = recip
zerosLin :: V2 -> Double
zerosLin (V2 a b) = - b / a
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad (V3 a b c)
| delta < 0 = Left ((-b/da) :+ (-sqrtNDelta/da), (-b/da) :+ (sqrtNDelta/da))
| otherwise = Right ((- b - sqrtDelta) / da, (-b + sqrtDelta) / da)
where
delta = b*b - 4 * a * c
sqrtDelta = sqrt delta
sqrtNDelta = sqrt (- delta)
da = 2 * a
optimaQuad :: V3 -> Double
optimaQuad (V3 a b _) = zerosLin (V2 (2 * a) b)
data V2 = V2 !Double !Double
deriving (Eq, Show)
instance Add V2 where
zero = V2 0 0
add (V2 x y) (V2 x' y') = V2 (x + x') (y + y')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Mult Double V2 V2 where
mult k (V2 x y) = V2 (k * x) (k * y)
{-# INLINE mult #-}
data M22 = M22 !Double !Double !Double !Double
deriving (Eq, Show)
instance Add M22 where
zero = M22 0 0 0 0
add (M22 a b c d) (M22 a' b' c' d') = M22 (a + a') (b + b') (c + c') (d + d')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Eye M22 where
eye = M22 1 0 0 1
{-# INLINE eye #-}
instance Det M22 where det = det2
instance Inv M22 where inv = inv2
instance Mult Double M22 M22 where
mult k (M22 a b c d) = M22 (k * a) (k * b) (k * c) (k * d)
{-# INLINE mult #-}
instance Mult M22 V2 V2 where
mult (M22 a b c d) (V2 u v) = V2 (a * u + b * v) (c * u + d * v)
{-# INLINE mult #-}
instance Mult M22 M22 M22 where
mult (M22 a b c d) (M22 x y z w) = M22
(a * x + b * z) (a * y + b * w)
(c * x + d * z) (c * y + d * w)
{-# INLINE mult #-}
det2 :: M22 -> Double
det2 (M22 a b c d) = a * d - b * c
{-# INLINE det2 #-}
inv2 :: M22 -> M22
inv2 m@(M22 a b c d) = M22
( d / det) (- b / det)
(- c / det) ( a / det)
where
det = det2 m
{-# INLINE inv2 #-}
data V3 = V3 !Double !Double !Double
deriving (Eq, Show)
instance Add V3 where
zero = V3 0 0 0
add (V3 x y z) (V3 x' y' z') = V3 (x + x') (y + y') (z + z')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Mult Double V3 V3 where
mult k (V3 x y z) = V3 (k * x) (k * y) (k * z)
{-# INLINE mult #-}
data M33 = M33
!Double !Double !Double
!Double !Double !Double
!Double !Double !Double
deriving (Eq, Show)
instance Add M33 where
zero = M33 0 0 0 0 0 0 0 0 0
add (M33 a b c d e f g h i) (M33 a' b' c' d' e' f' g' h' i') = M33
(a + a') (b + b') (c + c')
(d + d') (e + e') (f + f')
(g + g') (h + h') (i + i')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Eye M33 where
eye = M33 1 0 0 0 1 0 0 0 1
{-# INLINE eye #-}
instance Det M33 where det = det3
instance Inv M33 where inv = inv3
instance Mult Double M33 M33 where
mult k (M33 a b c d e f g h i) = M33
(k * a) (k * b) (k * c)
(k * d) (k * e) (k * f)
(k * g) (k * h) (k * i)
{-# INLINE mult #-}
instance Mult M33 V3 V3 where
mult (M33 a b c
d e f
g h i) (V3 u v w) = V3
(a * u + b * v + c * w)
(d * u + e * v + f * w)
(g * u + h * v + i * w)
{-# INLINE mult #-}
det3 :: M33 -> Double
det3 (M33 a b c
d e f
g h i)
= a * (e*i-f*h) - d * (b*i-c*h) + g * (b*f-c*e)
{-# INLINE det3 #-}
inv3 :: M33 -> M33
inv3 m@(M33 a b c
d e f
g h i)
= M33 a' b' c'
d' e' f'
g' h' i'
where
a' = cofactor e f h i / det
b' = cofactor c b i h / det
c' = cofactor b c e f / det
d' = cofactor f d i g / det
e' = cofactor a c g i / det
f' = cofactor c a f d / det
g' = cofactor d e g h / det
h' = cofactor b a h g / det
i' = cofactor a b d e / det
cofactor q r s t = det2 (M22 q r s t)
det = det3 m
{-# INLINE inv3 #-}
linear :: (Foldable' xs x, IsDoublePair x) => xs -> V2
linear data_ = mult (inv2 (M22 x n x2 x)) (V2 y xy)
where
K2 n' (V2 x _) (V2 x2 _) (V2 y _) (V2 xy _) = kahan2 data_
n :: Double
n = fromIntegral n'
quadratic :: (Foldable' xs x, IsDoublePair x) => xs -> V3
quadratic data_ = mult (inv3 (M33 x2 x n x3 x2 x x4 x3 x2)) (V3 y xy x2y)
where
K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) = kahan3 data_
n :: Double
n = fromIntegral n'
quadraticAndLinear :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V2)
quadraticAndLinear data_ =
( mult (inv3 (M33 x2 x n x3 x2 x x4 x3 x2)) (V3 y xy x2y)
, mult (inv2 (M22 x n x2 x)) (V2 y xy)
)
where
K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) = kahan3 data_
n :: Double
n = fromIntegral n'
class Foldable' xs x | xs -> x where
foldl' :: (b -> x -> b) -> b -> xs -> b
instance Foldable' [a] a where foldl' = L.foldl'
instance Foldable' (V.Vector a) a where foldl' = V.foldl'
instance U.Unbox a => Foldable' (U.Vector a) a where foldl' = U.foldl'
class IsDoublePair dp where
withDP :: dp -> (Double -> Double -> r) -> r
makeDP :: Double -> Double -> dp
instance IsDoublePair V2 where
withDP (V2 x y) k = k x y
makeDP = V2
instance (a ~ Double, b ~ Double) => IsDoublePair (a, b) where
withDP ~(x, y) k = k x y
makeDP = (,)
data Kahan2 = K2
{ k2n :: {-# UNPACK #-} !Int
, k2x :: {-# UNPACK #-} !V2
, k2x2 :: {-# UNPACK #-} !V2
, k2y :: {-# UNPACK #-} !V2
, k2xy :: {-# UNPACK #-} !V2
}
zeroKahan2 :: Kahan2
zeroKahan2 = K2 0 zero zero zero zero
addKahan :: V2 -> Double -> V2
addKahan (V2 acc c) i =
let y = i - c
t = acc + y
in V2 t ((t - acc) - y)
kahan2 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan2
kahan2 = foldl' f zeroKahan2 where
f (K2 n x x2 y xy) uv = withDP uv $ \u v -> K2
(succ n)
(addKahan x u)
(addKahan x2 (u * u))
(addKahan y v)
(addKahan xy (u * v))
data Kahan3 = K3
{ k3n :: {-# UNPACK #-} !Int
, k3x :: {-# UNPACK #-} !V2
, k3x2 :: {-# UNPACK #-} !V2
, k3x3 :: {-# UNPACK #-} !V2
, k3x4 :: {-# UNPACK #-} !V2
, k3y :: {-# UNPACK #-} !V2
, k3xy :: {-# UNPACK #-} !V2
, k3x2y :: {-# UNPACK #-} !V2
}
zeroKahan3 :: Kahan3
zeroKahan3 = K3 0 zero zero zero zero zero zero zero
kahan3 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 = foldl' f zeroKahan3 where
f (K3 n x x2 x3 x4 y xy x2y) uv = withDP uv $ \u v ->
let u2 = u * u
in K3
(succ n)
(addKahan x u)
(addKahan x2 u2)
(addKahan x3 (u * u2))
(addKahan x4 (u2 * u2))
(addKahan y v)
(addKahan xy (u * v))
(addKahan x2y (u2 * v))