{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs #-}
module Math.Regression.Simple (
linear,
linearWithErrors,
quadratic,
quadraticWithErrors,
quadraticAndLinear,
quadraticAndLinearWithErrors,
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 :: Double
eye = Double
1
instance Add Double where
zero :: Double
zero = Double
0
add :: Double -> Double -> Double
add = Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)
instance Det Double where
det :: Double -> Double
det = Double -> Double
forall a. a -> a
id
instance Inv Double where
inv :: Double -> Double
inv = Double -> Double
forall a. Fractional a => a -> a
recip
zerosLin :: V2 -> Double
zerosLin :: V2 -> Double
zerosLin (V2 Double
a Double
b) = Double -> Double
forall a. Num a => a -> a
negate (Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a)
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad (V3 Double
a Double
b Double
c)
| Double
delta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = (Complex Double, Complex Double)
-> Either (Complex Double, Complex Double) (Double, Double)
forall a b. a -> Either a b
Left ((-Double
bDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (-Double
sqrtNDeltaDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da), (-Double
bDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
sqrtNDeltaDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da))
| Bool
otherwise = (Double, Double)
-> Either (Complex Double, Complex Double) (Double, Double)
forall a b. b -> Either a b
Right ((- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sqrtDelta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
da, (-Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sqrtDelta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
da)
where
delta :: Double
delta = Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
sqrtDelta :: Double
sqrtDelta = Double -> Double
forall a. Floating a => a -> a
sqrt Double
delta
sqrtNDelta :: Double
sqrtNDelta = Double -> Double
forall a. Floating a => a -> a
sqrt (- Double
delta)
da :: Double
da = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a
optimaQuad :: V3 -> Double
optimaQuad :: V3 -> Double
optimaQuad (V3 Double
a Double
b Double
_) = V2 -> Double
zerosLin (Double -> Double -> V2
V2 (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) Double
b)
data V2 = V2 !Double !Double
deriving (V2 -> V2 -> Bool
(V2 -> V2 -> Bool) -> (V2 -> V2 -> Bool) -> Eq V2
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: V2 -> V2 -> Bool
$c/= :: V2 -> V2 -> Bool
== :: V2 -> V2 -> Bool
$c== :: V2 -> V2 -> Bool
Eq, Int -> V2 -> ShowS
[V2] -> ShowS
V2 -> String
(Int -> V2 -> ShowS)
-> (V2 -> String) -> ([V2] -> ShowS) -> Show V2
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [V2] -> ShowS
$cshowList :: [V2] -> ShowS
show :: V2 -> String
$cshow :: V2 -> String
showsPrec :: Int -> V2 -> ShowS
$cshowsPrec :: Int -> V2 -> ShowS
Show)
instance Add V2 where
zero :: V2
zero = Double -> Double -> V2
V2 Double
0 Double
0
add :: V2 -> V2 -> V2
add (V2 Double
x Double
y) (V2 Double
x' Double
y') = Double -> Double -> V2
V2 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x') (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Mult Double V2 V2 where
mult :: Double -> V2 -> V2
mult Double
k (V2 Double
x Double
y) = Double -> Double -> V2
V2 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
{-# INLINE mult #-}
data M22 = M22 !Double !Double !Double !Double
deriving (M22 -> M22 -> Bool
(M22 -> M22 -> Bool) -> (M22 -> M22 -> Bool) -> Eq M22
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: M22 -> M22 -> Bool
$c/= :: M22 -> M22 -> Bool
== :: M22 -> M22 -> Bool
$c== :: M22 -> M22 -> Bool
Eq, Int -> M22 -> ShowS
[M22] -> ShowS
M22 -> String
(Int -> M22 -> ShowS)
-> (M22 -> String) -> ([M22] -> ShowS) -> Show M22
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [M22] -> ShowS
$cshowList :: [M22] -> ShowS
show :: M22 -> String
$cshow :: M22 -> String
showsPrec :: Int -> M22 -> ShowS
$cshowsPrec :: Int -> M22 -> ShowS
Show)
instance Add M22 where
zero :: M22
zero = Double -> Double -> Double -> Double -> M22
M22 Double
0 Double
0 Double
0 Double
0
add :: M22 -> M22 -> M22
add (M22 Double
a Double
b Double
c Double
d) (M22 Double
a' Double
b' Double
c' Double
d') = Double -> Double -> Double -> Double -> M22
M22 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a') (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b') (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c') (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Eye M22 where
eye :: M22
eye = Double -> Double -> Double -> Double -> M22
M22 Double
1 Double
0 Double
0 Double
1
{-# INLINE eye #-}
instance Det M22 where det :: M22 -> Double
det = M22 -> Double
det2
instance Inv M22 where inv :: M22 -> M22
inv = M22 -> M22
inv2
instance Mult Double M22 M22 where
mult :: Double -> M22 -> M22
mult Double
k (M22 Double
a Double
b Double
c Double
d) = Double -> Double -> Double -> Double -> M22
M22 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)
{-# INLINE mult #-}
instance Mult M22 V2 V2 where
mult :: M22 -> V2 -> V2
mult (M22 Double
a Double
b Double
c Double
d) (V2 Double
u Double
v) = Double -> Double -> V2
V2 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v) (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
{-# INLINE mult #-}
instance Mult M22 M22 M22 where
mult :: M22 -> M22 -> M22
mult (M22 Double
a Double
b Double
c Double
d) (M22 Double
x Double
y Double
z Double
w) = Double -> Double -> Double -> Double -> M22
M22
(Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
(Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
{-# INLINE mult #-}
det2 :: M22 -> Double
det2 :: M22 -> Double
det2 (M22 Double
a Double
b Double
c Double
d) = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
{-# INLINE det2 #-}
inv2 :: M22 -> M22
inv2 :: M22 -> M22
inv2 m :: M22
m@(M22 Double
a Double
b Double
c Double
d) = Double -> Double -> Double -> Double -> M22
M22
( Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det) (- Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det)
(- Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det) ( Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det)
where
det :: Double
det = M22 -> Double
det2 M22
m
{-# INLINE inv2 #-}
data V3 = V3 !Double !Double !Double
deriving (V3 -> V3 -> Bool
(V3 -> V3 -> Bool) -> (V3 -> V3 -> Bool) -> Eq V3
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: V3 -> V3 -> Bool
$c/= :: V3 -> V3 -> Bool
== :: V3 -> V3 -> Bool
$c== :: V3 -> V3 -> Bool
Eq, Int -> V3 -> ShowS
[V3] -> ShowS
V3 -> String
(Int -> V3 -> ShowS)
-> (V3 -> String) -> ([V3] -> ShowS) -> Show V3
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [V3] -> ShowS
$cshowList :: [V3] -> ShowS
show :: V3 -> String
$cshow :: V3 -> String
showsPrec :: Int -> V3 -> ShowS
$cshowsPrec :: Int -> V3 -> ShowS
Show)
instance Add V3 where
zero :: V3
zero = Double -> Double -> Double -> V3
V3 Double
0 Double
0 Double
0
add :: V3 -> V3 -> V3
add (V3 Double
x Double
y Double
z) (V3 Double
x' Double
y' Double
z') = Double -> Double -> Double -> V3
V3 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x') (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y') (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
z')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Mult Double V3 V3 where
mult :: Double -> V3 -> V3
mult Double
k (V3 Double
x Double
y Double
z) = Double -> Double -> Double -> V3
V3 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z)
{-# INLINE mult #-}
data M33 = M33
!Double !Double !Double
!Double !Double !Double
!Double !Double !Double
deriving (M33 -> M33 -> Bool
(M33 -> M33 -> Bool) -> (M33 -> M33 -> Bool) -> Eq M33
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: M33 -> M33 -> Bool
$c/= :: M33 -> M33 -> Bool
== :: M33 -> M33 -> Bool
$c== :: M33 -> M33 -> Bool
Eq, Int -> M33 -> ShowS
[M33] -> ShowS
M33 -> String
(Int -> M33 -> ShowS)
-> (M33 -> String) -> ([M33] -> ShowS) -> Show M33
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [M33] -> ShowS
$cshowList :: [M33] -> ShowS
show :: M33 -> String
$cshow :: M33 -> String
showsPrec :: Int -> M33 -> ShowS
$cshowsPrec :: Int -> M33 -> ShowS
Show)
instance Add M33 where
zero :: M33
zero = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0
add :: M33 -> M33 -> M33
add (M33 Double
a Double
b Double
c Double
d Double
e Double
f Double
g Double
h Double
i) (M33 Double
a' Double
b' Double
c' Double
d' Double
e' Double
f' Double
g' Double
h' Double
i') = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33
(Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a') (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b') (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c')
(Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d') (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e') (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f')
(Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g') (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h') (Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
i')
{-# INLINE zero #-}
{-# INLINE add #-}
instance Eye M33 where
eye :: M33
eye = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
1 Double
0 Double
0 Double
0 Double
1 Double
0 Double
0 Double
0 Double
1
{-# INLINE eye #-}
instance Det M33 where det :: M33 -> Double
det = M33 -> Double
det3
instance Inv M33 where inv :: M33 -> M33
inv = M33 -> M33
inv3
instance Mult Double M33 M33 where
mult :: Double -> M33 -> M33
mult Double
k (M33 Double
a Double
b Double
c Double
d Double
e Double
f Double
g Double
h Double
i) = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33
(Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c)
(Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f)
(Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
i)
{-# INLINE mult #-}
instance Mult M33 V3 V3 where
mult :: M33 -> V3 -> V3
mult (M33 Double
a Double
b Double
c
Double
d Double
e Double
f
Double
g Double
h Double
i) (V3 Double
u Double
v Double
w) = Double -> Double -> Double -> V3
V3
(Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
(Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
(Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
{-# INLINE mult #-}
det3 :: M33 -> Double
det3 :: M33 -> Double
det3 (M33 Double
a Double
b Double
c
Double
d Double
e Double
f
Double
g Double
h Double
i)
= Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
iDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
fDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
iDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
e)
{-# INLINE det3 #-}
inv3 :: M33 -> M33
inv3 :: M33 -> M33
inv3 m :: M33
m@(M33 Double
a Double
b Double
c
Double
d Double
e Double
f
Double
g Double
h Double
i)
= Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
a' Double
b' Double
c'
Double
d' Double
e' Double
f'
Double
g' Double
h' Double
i'
where
a' :: Double
a' = Double -> Double -> Double -> Double -> Double
cofactor Double
e Double
f Double
h Double
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
b' :: Double
b' = Double -> Double -> Double -> Double -> Double
cofactor Double
c Double
b Double
i Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
c' :: Double
c' = Double -> Double -> Double -> Double -> Double
cofactor Double
b Double
c Double
e Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
d' :: Double
d' = Double -> Double -> Double -> Double -> Double
cofactor Double
f Double
d Double
i Double
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
e' :: Double
e' = Double -> Double -> Double -> Double -> Double
cofactor Double
a Double
c Double
g Double
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
f' :: Double
f' = Double -> Double -> Double -> Double -> Double
cofactor Double
c Double
a Double
f Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
g' :: Double
g' = Double -> Double -> Double -> Double -> Double
cofactor Double
d Double
e Double
g Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
h' :: Double
h' = Double -> Double -> Double -> Double -> Double
cofactor Double
b Double
a Double
h Double
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
i' :: Double
i' = Double -> Double -> Double -> Double -> Double
cofactor Double
a Double
b Double
d Double
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
cofactor :: Double -> Double -> Double -> Double -> Double
cofactor Double
q Double
r Double
s Double
t = M22 -> Double
det2 (Double -> Double -> Double -> Double -> M22
M22 Double
q Double
r Double
s Double
t)
det :: Double
det = M33 -> Double
det3 M33
m
{-# INLINE inv3 #-}
linear :: (Foldable' xs x, IsDoublePair x) => xs -> V2
linear :: xs -> V2
linear = (V2, V2) -> V2
forall a b. (a, b) -> a
fst ((V2, V2) -> V2) -> (xs -> (V2, V2)) -> xs -> V2
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> (V2, V2)
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> (V2, V2)
linearWithErrors
linearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V2, V2)
linearWithErrors :: xs -> (V2, V2)
linearWithErrors = Kahan2 -> (V2, V2)
linearImpl (Kahan2 -> (V2, V2)) -> (xs -> Kahan2) -> xs -> (V2, V2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> Kahan2
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan2
kahan2
linearImpl :: Kahan2 -> (V2, V2)
linearImpl :: Kahan2 -> (V2, V2)
linearImpl (K2 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
y Double
_) (V2 Double
y2 Double
_) (V2 Double
xy Double
_)) =
(V2
params, V2
errors)
where
n :: Double
n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
matrix :: M22
matrix@(M22 Double
a11 Double
_ Double
_ Double
a22) = M22 -> M22
inv2 (Double -> Double -> Double -> Double -> M22
M22 Double
x2 Double
x Double
x Double
n)
params :: V2
params@(V2 Double
a Double
b) = M22 -> V2 -> V2
forall a b c. Mult a b c => a -> b -> c
mult M22
matrix (Double -> Double -> V2
V2 Double
xy Double
y)
errors :: V2
errors = Double -> Double -> V2
V2 Double
sa Double
sb
err :: Double
err = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
sa :: Double
sa = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a11 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2))
sb :: Double
sb = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a22 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2))
quadratic :: (Foldable' xs x, IsDoublePair x) => xs -> V3
quadratic :: xs -> V3
quadratic xs
data_ = M33 -> V3 -> V3
forall a b c. Mult a b c => a -> b -> c
mult (M33 -> M33
inv3 (Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
x2 Double
x Double
n Double
x3 Double
x2 Double
x Double
x4 Double
x3 Double
x2)) (Double -> Double -> Double -> V3
V3 Double
y Double
xy Double
x2y)
where
K3 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
x3 Double
_) (V2 Double
x4 Double
_) (V2 Double
y Double
_) (V2 Double
xy Double
_) (V2 Double
x2y Double
_) (V2 Double
_y2 Double
_) = xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 xs
data_
n :: Double
n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
quadraticWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V3)
quadraticWithErrors :: xs -> (V3, V3)
quadraticWithErrors = Kahan3 -> (V3, V3)
quadraticImpl (Kahan3 -> (V3, V3)) -> (xs -> Kahan3) -> xs -> (V3, V3)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3
quadraticImpl :: Kahan3 -> (V3, V3)
quadraticImpl :: Kahan3 -> (V3, V3)
quadraticImpl (K3 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
x3 Double
_) (V2 Double
x4 Double
_) (V2 Double
y Double
_) (V2 Double
xy Double
_) (V2 Double
x2y Double
_) (V2 Double
y2 Double
_)) =
(V3
params, V3
errors)
where
n :: Double
n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
matrix :: M33
matrix@(M33 Double
a11 Double
_ Double
_
Double
_ Double
a22 Double
_
Double
_ Double
_ Double
a33) = M33 -> M33
inv3 (Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
x4 Double
x3 Double
x2 Double
x3 Double
x2 Double
x Double
x2 Double
x Double
n)
params :: V3
params@(V3 Double
a Double
b Double
c) = M33 -> V3 -> V3
forall a b c. Mult a b c => a -> b -> c
mult M33
matrix (Double -> Double -> Double -> V3
V3 Double
x2y Double
xy Double
y)
errors :: V3
errors = Double -> Double -> Double -> V3
V3 Double
sa Double
sb Double
sc
err :: Double
err = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x2y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
sa :: Double
sa = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a11 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))
sb :: Double
sb = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a22 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))
sc :: Double
sc = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a33 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))
quadraticAndLinear :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V2)
quadraticAndLinear :: xs -> (V3, V2)
quadraticAndLinear = ((V3, V2), (V3, V2)) -> (V3, V2)
forall a b. (a, b) -> a
fst (((V3, V2), (V3, V2)) -> (V3, V2))
-> (xs -> ((V3, V2), (V3, V2))) -> xs -> (V3, V2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> ((V3, V2), (V3, V2))
forall xs x.
(Foldable' xs x, IsDoublePair x) =>
xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors
quadraticAndLinearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors :: xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors xs
data_ =
((V3
paramsQ, V2
paramsL), (V3
errorsQ, V2
errorsL))
where
k3 :: Kahan3
k3@(K3 Int
n V2
x V2
x2 V2
x3 V2
x4 V2
y V2
xy V2
x2y V2
y2) = xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 xs
data_
k2 :: Kahan2
k2 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2 Int
n V2
x V2
x2 V2
y V2
y2 V2
xy
(V2
paramsL, V2
errorsL) = Kahan2 -> (V2, V2)
linearImpl Kahan2
k2
(V3
paramsQ, V3
errorsQ) = Kahan3 -> (V3, V3)
quadraticImpl Kahan3
k3
class Foldable' xs x | xs -> x where
foldl' :: (b -> x -> b) -> b -> xs -> b
instance Foldable' [a] a where foldl' :: (b -> a -> b) -> b -> [a] -> b
foldl' = (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
L.foldl'
instance Foldable' (V.Vector a) a where foldl' :: (b -> a -> b) -> b -> Vector a -> b
foldl' = (b -> a -> b) -> b -> Vector a -> b
forall a b. (a -> b -> a) -> a -> Vector b -> a
V.foldl'
instance U.Unbox a => Foldable' (U.Vector a) a where foldl' :: (b -> a -> b) -> b -> Vector a -> b
foldl' = (b -> a -> b) -> b -> Vector a -> b
forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
U.foldl'
class IsDoublePair dp where
withDP :: dp -> (Double -> Double -> r) -> r
makeDP :: Double -> Double -> dp
instance IsDoublePair V2 where
withDP :: V2 -> (Double -> Double -> r) -> r
withDP (V2 Double
x Double
y) Double -> Double -> r
k = Double -> Double -> r
k Double
x Double
y
makeDP :: Double -> Double -> V2
makeDP = Double -> Double -> V2
V2
instance (a ~ Double, b ~ Double) => IsDoublePair (a, b) where
withDP :: (a, b) -> (Double -> Double -> r) -> r
withDP ~(a
x, b
y) Double -> Double -> r
k = Double -> Double -> r
k a
Double
x b
Double
y
makeDP :: Double -> Double -> (a, b)
makeDP = (,)
data Kahan2 = K2
{ Kahan2 -> Int
k2n :: {-# UNPACK #-} !Int
, Kahan2 -> V2
k2x :: {-# UNPACK #-} !V2
, Kahan2 -> V2
k2x2 :: {-# UNPACK #-} !V2
, Kahan2 -> V2
k2y :: {-# UNPACK #-} !V2
, Kahan2 -> V2
k2y2 :: {-# UNPACK #-} !V2
, Kahan2 -> V2
k2xy :: {-# UNPACK #-} !V2
}
zeroKahan2 :: Kahan2
zeroKahan2 :: Kahan2
zeroKahan2 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2 Int
0 V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero
addKahan :: V2 -> Double -> V2
addKahan :: V2 -> Double -> V2
addKahan (V2 Double
acc Double
c) Double
i =
let y :: Double
y = Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c
t :: Double
t = Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y
in Double -> Double -> V2
V2 Double
t ((Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
acc) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)
kahan2 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan2
kahan2 :: xs -> Kahan2
kahan2 = (Kahan2 -> x -> Kahan2) -> Kahan2 -> xs -> Kahan2
forall xs x b. Foldable' xs x => (b -> x -> b) -> b -> xs -> b
foldl' Kahan2 -> x -> Kahan2
forall dp. IsDoublePair dp => Kahan2 -> dp -> Kahan2
f Kahan2
zeroKahan2 where
f :: Kahan2 -> dp -> Kahan2
f (K2 Int
n V2
x V2
x2 V2
y V2
y2 V2
xy) dp
uv = dp -> (Double -> Double -> Kahan2) -> Kahan2
forall dp r. IsDoublePair dp => dp -> (Double -> Double -> r) -> r
withDP dp
uv ((Double -> Double -> Kahan2) -> Kahan2)
-> (Double -> Double -> Kahan2) -> Kahan2
forall a b. (a -> b) -> a -> b
$ \Double
u Double
v -> Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2
(Int -> Int
forall a. Enum a => a -> a
succ Int
n)
(V2 -> Double -> V2
addKahan V2
x Double
u)
(V2 -> Double -> V2
addKahan V2
x2 (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u))
(V2 -> Double -> V2
addKahan V2
y Double
v)
(V2 -> Double -> V2
addKahan V2
y2 (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
(V2 -> Double -> V2
addKahan V2
xy (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
data Kahan3 = K3
{ Kahan3 -> Int
k3n :: {-# UNPACK #-} !Int
, Kahan3 -> V2
k3x :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3x2 :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3x3 :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3x4 :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3y :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3xy :: {-# UNPACK #-} !V2
, Kahan3 -> V2
k3x2y :: {-# UNPACK #-} !V2
, Kahan3 -> V2
ky2 :: {-# UNPACK #-} !V2
}
zeroKahan3 :: Kahan3
zeroKahan3 :: Kahan3
zeroKahan3 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan3
K3 Int
0 V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero
kahan3 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 :: xs -> Kahan3
kahan3 = (Kahan3 -> x -> Kahan3) -> Kahan3 -> xs -> Kahan3
forall xs x b. Foldable' xs x => (b -> x -> b) -> b -> xs -> b
foldl' Kahan3 -> x -> Kahan3
forall dp. IsDoublePair dp => Kahan3 -> dp -> Kahan3
f Kahan3
zeroKahan3 where
f :: Kahan3 -> dp -> Kahan3
f (K3 Int
n V2
x V2
x2 V2
x3 V2
x4 V2
y V2
xy V2
x2y V2
y2) dp
uv = dp -> (Double -> Double -> Kahan3) -> Kahan3
forall dp r. IsDoublePair dp => dp -> (Double -> Double -> r) -> r
withDP dp
uv ((Double -> Double -> Kahan3) -> Kahan3)
-> (Double -> Double -> Kahan3) -> Kahan3
forall a b. (a -> b) -> a -> b
$ \Double
u Double
v ->
let u2 :: Double
u2 = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u
in Int -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan3
K3
(Int -> Int
forall a. Enum a => a -> a
succ Int
n)
(V2 -> Double -> V2
addKahan V2
x Double
u)
(V2 -> Double -> V2
addKahan V2
x2 Double
u2)
(V2 -> Double -> V2
addKahan V2
x3 (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u2))
(V2 -> Double -> V2
addKahan V2
x4 (Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u2))
(V2 -> Double -> V2
addKahan V2
y Double
v)
(V2 -> Double -> V2
addKahan V2
xy (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
(V2 -> Double -> V2
addKahan V2
x2y (Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
(V2 -> Double -> V2
addKahan V2
y2 (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))