{-# LANGUAGE CPP                    #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE RecordWildCards        #-}
-- |
--
-- @regression-simple@ provides (hopefully) simple regression functions.
--
-- The @'linear' :: Foldable f => (a -> (Double, Double)) -> f a -> 'V2'@
-- is the simplest one.
--
-- There are variants with weights, y-errors, and x and y-errors.
-- In addition, package includes Levenberg–Marquardt algorithm implementation
-- to fit arbitrary functions (with one, two or three parameters),
-- as long as you can give their partial derivatives as well (@ad@ package is handy for that).
--
-- For multiple independent variable ordinary least squares
-- or Levenberg-Marquard with functions with \> 3 parameter you should look elsewhere.
--
-- Package has been tested to return similar results as @fit@ functionality in @gnuplot@
-- (L-M doesn't always converge to exactly the same points in parameter space).
--
module Math.Regression.Simple (
    -- * Linear regression
    linear,
    linearFit,
    linearWithWeights,
    linearWithYerrors,
    linearWithXYerrors,
    -- ** Step-by-step interface
    linearFit',
    LinRegAcc (..),
    zeroLinRegAcc,
    addLinReg,
    addLinRegW,
    -- * Quadratic regression
    quadratic,
    quadraticFit,
    quadraticWithWeights,
    quadraticWithYerrors,
    quadraticWithXYerrors,
    -- ** Step-by-step interface
    quadraticFit',
    QuadRegAcc (..),
    zeroQuadRegAcc,
    addQuadReg,
    addQuadRegW,
    quadRegAccToLin,
    -- * Levenberg–Marquardt algorithm
    -- ** One parameter
    levenbergMarquardt1,
    levenbergMarquardt1WithWeights,
    levenbergMarquardt1WithYerrors,
    levenbergMarquardt1WithXYerrors,
    -- ** Two parameters
    levenbergMarquardt2,
    levenbergMarquardt2WithWeights,
    levenbergMarquardt2WithYerrors,
    levenbergMarquardt2WithXYerrors,
    -- ** Three parameters
    levenbergMarquardt3,
    levenbergMarquardt3WithWeights,
    levenbergMarquardt3WithYerrors,
    levenbergMarquardt3WithXYerrors,
    -- * Auxiliary types
    Fit (..),
    V2 (..),
    V3 (..),
) where

import Control.DeepSeq (NFData (..))

import qualified Data.Foldable      as F
import qualified Data.List.NonEmpty as NE

import Math.Regression.Simple.LinAlg
import Numeric.KBN

-- $setup
-- >>> :set -XDeriveFunctor -XDeriveFoldable -XDeriveTraversable
--
-- >>> import Numeric (showFFloat)
-- >>> import Data.List.NonEmpty (NonEmpty (..))
-- >>> import qualified Data.List.NonEmpty as NE
--
-- Don't show too much decimal digits
--
-- >>> let showDouble x = showFFloat (Just (min 5 (5 - ceiling (logBase 10 x)))) (x :: Double)
-- >>> showDouble 123.456 ""
-- "123.46"
--
-- >>> showDouble 1234567 ""
-- "1234567"
--
-- >>> showDouble 123.4567890123456789 ""
-- "123.46"
--
-- >>> showDouble 0.0000000000000012345 ""
-- "0.00000"
--
-- >>> newtype PP a = PP a
-- >>> class Show' a where showsPrec' :: Int -> a -> ShowS
-- >>> instance Show' a => Show (PP a) where showsPrec d (PP x) = showsPrec' d x
-- >>> instance Show' Double where showsPrec' d x = if x < 0 then showParen (d > 6) (showChar '-' . showDouble (negate x)) else showDouble x
-- >>> instance Show' Int where showsPrec' = showsPrec
-- >>> instance (Show' a, Show' b) => Show' (a, b) where showsPrec' d (x, y) = showParen True $ showsPrec' 0 x . showString ", " . showsPrec' 0 y
-- >>> instance Show' v => Show' (Fit v) where showsPrec' d (Fit p e ndf wssr) = showParen (d > 10) $ showString "Fit " . showsPrec' 11 p . showChar ' ' . showsPrec' 11 e . showChar ' ' . showsPrec' 11 ndf . showChar ' ' . showsPrec' 11 wssr
-- >>> instance Show' V2 where showsPrec' d (V2 x y)   = showParen (d > 10) $ showString "V2 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y
-- >>> instance Show' V3 where showsPrec' d (V3 x y z) = showParen (d > 10) $ showString "V3 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y . showChar ' ' . showsPrec' 11 z
-- >>> instance Show' a => Show' [a] where showsPrec' _ [] = id; showsPrec' _ [x] = showsPrec' 0 x; showsPrec' _ (x:xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs
-- >>> instance Show' a => Show' (NonEmpty a) where showsPrec' _ (x :| []) = showsPrec' 0 x; showsPrec' _ (x:|xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs
-- >>>
--
-- Inputs:
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
--
-- >>> let sq z = z * z
--

-------------------------------------------------------------------------------
-- Linear
-------------------------------------------------------------------------------

-- | Linear regression.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> PP $ linear id input1
-- V2 2.0000 1.00000
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ linear id input2
-- V2 2.0063 0.88685
--
linear :: F.Foldable f => (a -> (Double, Double)) -> f a -> V2
linear :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> V2
linear a -> (Double, Double)
f = forall v. Fit v -> v
fitParams forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> Fit V2
linearFit a -> (Double, Double)
f

-- | Like 'linear' but returns complete 'Fit'.
--
-- To get confidence intervals you should multiply the errors
-- by @quantile (studentT (n - 2)) ci'@ from @statistics@ package
-- or similar.
-- For big @n@ using value 1 gives 68% interval and using value 2 gives 95% confidence interval.
-- See https://en.wikipedia.org/wiki/Student%27s_t-distribution#Table_of_selected_values
-- (@quantile@ calculates one-sided values, you need two-sided, thus adjust @ci@ value).
--
-- The first input is perfect fit:
--
-- >>> let fit = linearFit id input1
-- >>> PP fit
-- Fit (V2 2.0000 1.00000) (V2 0.00000 0.00000) 1 0.00000
--
-- The second input is quite good:
--
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- But the third input isn't so much,
-- standard error of a slope parameter is 20%.
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ linearFit id input3
-- Fit (V2 3.0000 1.00000) (V2 0.63246 1.1832) 2 4.0000
--
linearFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V2
linearFit :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> Fit V2
linearFit a -> (Double, Double)
f = LinRegAcc -> Fit V2
linearFit' forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> LinRegAcc
linRegAcc a -> (Double, Double)
f

-- | Weighted linear regression.
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)]
-- >>> PP $ linearWithWeights id input2w
-- Fit (V2 2.0060 0.86993) (V2 0.12926 0.23696) 3 0.22074
--
linearWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithWeights :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithWeights a -> (Double, Double, Double)
f = LinRegAcc -> Fit V2
linearFit' forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> LinRegAcc
linRegAccW a -> (Double, Double, Double)
f

-- | Linear regression with y-errors.
--
-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 1.9)]
-- >>> let fit = linearWithYerrors id input2y
-- >>> PP fit
-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930
--
-- When we know actual y-errors, we can calculate the Q-value using @statistics@ package:
--
-- >>> import qualified Statistics.Distribution            as S
-- >>> import qualified Statistics.Distribution.ChiSquared as S
-- >>> S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit)
-- 0.446669639443138
--
-- or using @math-functions@
--
-- >>> import Numeric.SpecFunctions (incompleteGamma)
-- >>> incompleteGamma (fromIntegral (fitNDF fit) / 2) (fitWSSR fit / 2)
-- 0.446669639443138
--
-- It is not uncommon to deem acceptable on equal terms any models with, say, Q > 0.001.
-- If Q is too large, too near to 1 is most likely caused by overestimating
-- the y-errors.
--
--
linearWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithYerrors :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithYerrors a -> (Double, Double, Double)
f = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithWeights a -> (Double, Double, Double)
f' where
    f' :: a -> (Double, Double, Double)
f' a
a = case a -> (Double, Double, Double)
f a
a of
        (Double
x, Double
y, Double
dy) -> (Double
x, Double
y, forall a. Fractional a => a -> a
recip (Double
dy forall a. Num a => a -> a -> a
* Double
dy))

-- | Iterative linear regression with x and y errors.
--
-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/
--
-- >>> let input2xy = [(0.1, 1.2, 0.01, 0.12), (1.3, 3.1, 0.13, 0.31), (1.9, 4.9, 0.19, 0.49), (3.0, 7.1, 0.3, 0.71), (4.1, 9.0, 0.41, 1.9)]
-- >>> let fit :| fits = linearWithXYerrors id input2xy
--
-- First fit is done using 'linearWithYerrors':
--
-- >>> PP fit
-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930
--
-- After that the effective variance is used to refine the fit,
-- just a few iterations is often enough:
--
-- >>> PP $ take 3 fits
-- Fit (V2 1.9092 0.99251) (V2 0.12417 0.08412) 3 1.2992
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998
--
linearWithXYerrors
    :: F.Foldable f
    => (a -> (Double, Double, Double, Double))  -- ^ \(x_i, y_i, \delta x_i, \delta y_i\)
    -> f a                                      -- ^ data
    -> NE.NonEmpty (Fit V2)
linearWithXYerrors :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double, Double)) -> f a -> NonEmpty (Fit V2)
linearWithXYerrors a -> (Double, Double, Double, Double)
f f a
xs = forall b. (b -> b) -> b -> NonEmpty b
iterate1 Fit V2 -> Fit V2
go Fit V2
fit0 where
    fit0 :: Fit V2
fit0   = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithYerrors (\a
a -> case a -> (Double, Double, Double, Double)
f a
a of (Double
x,Double
y,Double
_,Double
dy) -> (Double
x,Double
y,Double
dy)) f a
xs
    go :: Fit V2 -> Fit V2
go Fit V2
fit = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithWeights (\a
a -> case a -> (Double, Double, Double, Double)
f a
a of (Double
x,Double
y,Double
dx,Double
dy) -> (Double
x,Double
y,forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq (Double
param1 forall a. Num a => a -> a -> a
* Double
dx) forall a. Num a => a -> a -> a
+ forall a. Num a => a -> a
sq Double
dy)) f a
xs where
        V2 Double
param1 Double
_ = forall v. Fit v -> v
fitParams Fit V2
fit

-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD
-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable)
-- >>> let linearF (H3 a b x) = a * x + b
-- >>> let lin' (V2 a b) (x, y, dx, dy) = case AD.grad' linearF (H3 a b x) of (f, H3 da db f') -> (y, f, V2 da db, recip $ sq (f' * dx) + sq dy)
--
-- >>> PP $ NE.last $ levenbergMarquardt2WithWeights lin' (V2 1 1) input2xy
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998

-- | Calculate linear fit from 'LinRegAcc'.
linearFit' :: LinRegAcc -> Fit V2
linearFit' :: LinRegAcc -> Fit V2
linearFit' LinRegAcc {Int
KBN
lra_y2 :: LinRegAcc -> KBN
lra_xy :: LinRegAcc -> KBN
lra_y :: LinRegAcc -> KBN
lra_x2 :: LinRegAcc -> KBN
lra_x :: LinRegAcc -> KBN
lra_w :: LinRegAcc -> KBN
lra_n :: LinRegAcc -> Int
lra_y2 :: KBN
lra_xy :: KBN
lra_y :: KBN
lra_x2 :: KBN
lra_x :: KBN
lra_w :: KBN
lra_n :: Int
..} = forall v. v -> v -> Int -> Double -> Fit v
Fit V2
params V2
errors Int
ndf Double
wssr where
    matrix :: SM22
matrix@(SM22 Double
a11 Double
_ Double
a22) = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> SM22
SM22 Double
x2 Double
x Double
w)
    params :: V2
params@(V2 Double
a Double
b)         = forall a b c. Mult a b c => a -> b -> c
mult SM22
matrix (Double -> Double -> V2
V2 Double
xy Double
y)

    errors :: V2
errors = Double -> Double -> V2
V2 Double
sa Double
sb

    -- ensure that error is always non-negative.
    -- Due rounding errors, in perfect fit situations it can be slightly negative.
    wssr :: Double
wssr = forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 forall a. Num a => a -> a -> a
- Double
a forall a. Num a => a -> a -> a
* Double
xy forall a. Num a => a -> a -> a
- Double
b forall a. Num a => a -> a -> a
* Double
y)
    ndf :: Int
ndf  = Int
lra_n forall a. Num a => a -> a -> a
- Int
2
    ndf' :: Double
ndf' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

    sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
    sb :: Double
sb   = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

    w :: Double
w  = KBN -> Double
getKBN KBN
lra_w
    x :: Double
x  = KBN -> Double
getKBN KBN
lra_x
    x2 :: Double
x2 = KBN -> Double
getKBN KBN
lra_x2
    y :: Double
y  = KBN -> Double
getKBN KBN
lra_y
    xy :: Double
xy = KBN -> Double
getKBN KBN
lra_xy
    y2 :: Double
y2 = KBN -> Double
getKBN KBN
lra_y2

    -- is it useful?
    -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y))

linRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> LinRegAcc
linRegAcc :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> LinRegAcc
linRegAcc a -> (Double, Double)
f = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LinRegAcc
acc a
a -> case a -> (Double, Double)
f a
a of (Double
x,Double
y) -> LinRegAcc -> Double -> Double -> LinRegAcc
addLinReg LinRegAcc
acc Double
x Double
y) LinRegAcc
zeroLinRegAcc

linRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> LinRegAcc
linRegAccW :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> LinRegAcc
linRegAccW a -> (Double, Double, Double)
f = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LinRegAcc
acc a
a -> case a -> (Double, Double, Double)
f a
a of (Double
x,Double
y,Double
w) -> LinRegAcc -> Double -> Double -> Double -> LinRegAcc
addLinRegW LinRegAcc
acc Double
x Double
y Double
w) LinRegAcc
zeroLinRegAcc

-------------------------------------------------------------------------------
-- Quadractic
-------------------------------------------------------------------------------

-- | Quadratic regression.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> quadratic id input1
-- V3 0.0 2.0 1.0
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ quadratic id input2
-- V3 (-0.00589) 2.0313 0.87155
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ quadratic id input3
-- V3 1.00000 0.00000 2.0000
--
quadratic :: F.Foldable f => (a -> (Double, Double)) -> f a -> V3
quadratic :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> V3
quadratic a -> (Double, Double)
f = forall v. Fit v -> v
fitParams forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> Fit V3
quadraticFit a -> (Double, Double)
f

-- | Like 'quadratic' but returns complete 'Fit'.
--
-- >>> PP $ quadraticFit id input2
-- Fit (V3 (-0.00589) 2.0313 0.87155) (V3 0.09281 0.41070 0.37841) 2 0.25910
--
-- >>> PP $ quadraticFit id input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
quadraticFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V3
quadraticFit :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> Fit V3
quadraticFit a -> (Double, Double)
f = QuadRegAcc -> Fit V3
quadraticFit' forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> QuadRegAcc
quadRegAcc a -> (Double, Double)
f

-- | Weighted quadratic regression.
--
-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)]
-- >>> PP $ quadraticWithWeights id input2w
-- Fit (V3 0.02524 1.9144 0.91792) (V3 0.10775 0.42106 0.35207) 2 0.21484
--
quadraticWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithWeights :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithWeights a -> (Double, Double, Double)
f = QuadRegAcc -> Fit V3
quadraticFit' forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> QuadRegAcc
quadRegAccW a -> (Double, Double, Double)
f

-- | Quadratic regression with y-errors.
--
-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 0.9)]
-- >>> PP $ quadraticWithYerrors id input2y
-- Fit (V3 0.08776 1.6667 1.0228) (V3 0.10131 0.31829 0.11917) 2 1.5398
--
quadraticWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithYerrors :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithYerrors a -> (Double, Double, Double)
f = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithWeights a -> (Double, Double, Double)
f' where
    f' :: a -> (Double, Double, Double)
f' a
a = case a -> (Double, Double, Double)
f a
a of
        (Double
x, Double
y, Double
dy) -> (Double
x, Double
y, forall a. Fractional a => a -> a
recip (Double
dy forall a. Num a => a -> a -> a
* Double
dy))

-- | Iterative quadratic regression with x and y errors.
--
-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/
--
quadraticWithXYerrors
    :: F.Foldable f
    => (a -> (Double, Double, Double, Double))  -- ^ \(x_i, y_i, \delta x_i, \delta y_i\)
    -> f a                                      -- ^ data
    -> NE.NonEmpty (Fit V3)
quadraticWithXYerrors :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double, Double)) -> f a -> NonEmpty (Fit V3)
quadraticWithXYerrors a -> (Double, Double, Double, Double)
f f a
xs = forall b. (b -> b) -> b -> NonEmpty b
iterate1 Fit V3 -> Fit V3
go Fit V3
fit0 where
    fit0 :: Fit V3
fit0   = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithYerrors (\a
a -> case a -> (Double, Double, Double, Double)
f a
a of (Double
x,Double
y,Double
_,Double
dy) -> (Double
x,Double
y,Double
dy)) f a
xs
    go :: Fit V3 -> Fit V3
go Fit V3
fit = forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithWeights (\a
a -> case a -> (Double, Double, Double, Double)
f a
a of (Double
x,Double
y,Double
dx,Double
dy) -> (Double
x,Double
y,forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq ((Double
2 forall a. Num a => a -> a -> a
* Double
p1 forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
+ Double
p2) forall a. Num a => a -> a -> a
* Double
dx) forall a. Num a => a -> a -> a
+ forall a. Num a => a -> a
sq Double
dy)) f a
xs where
        V3 Double
p1 Double
p2 Double
_ = forall v. Fit v -> v
fitParams Fit V3
fit

-- | Calculate quadratic fit from 'QuadRegAcc'.
quadraticFit' :: QuadRegAcc -> Fit V3
quadraticFit' :: QuadRegAcc -> Fit V3
quadraticFit' QuadRegAcc {Int
KBN
qra_y2 :: QuadRegAcc -> KBN
qra_x2y :: QuadRegAcc -> KBN
qra_xy :: QuadRegAcc -> KBN
qra_y :: QuadRegAcc -> KBN
qra_x4 :: QuadRegAcc -> KBN
qra_x3 :: QuadRegAcc -> KBN
qra_x2 :: QuadRegAcc -> KBN
qra_x :: QuadRegAcc -> KBN
qra_w :: QuadRegAcc -> KBN
qra_n :: QuadRegAcc -> Int
qra_y2 :: KBN
qra_x2y :: KBN
qra_xy :: KBN
qra_y :: KBN
qra_x4 :: KBN
qra_x3 :: KBN
qra_x2 :: KBN
qra_x :: KBN
qra_w :: KBN
qra_n :: Int
..} = forall v. v -> v -> Int -> Double -> Fit v
Fit V3
params V3
errors Int
ndf Double
wssr where
    matrix :: SM33
matrix@(SM33 Double
a11
                 Double
_   Double
a22
                 Double
_   Double
_   Double
a33) = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> Double -> Double -> Double -> SM33
SM33 Double
x4
                                          Double
x3 Double
x2
                                          Double
x2  Double
x Double
w)

    params :: V3
params@(V3 Double
a Double
b Double
c) = forall a b c. Mult a b c => a -> b -> c
mult SM33
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

    wssr :: Double
wssr = forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 forall a. Num a => a -> a -> a
- Double
a forall a. Num a => a -> a -> a
* Double
x2y forall a. Num a => a -> a -> a
- Double
b forall a. Num a => a -> a -> a
* Double
xy forall a. Num a => a -> a -> a
- Double
c forall a. Num a => a -> a -> a
* Double
y)
    ndf :: Int
ndf  = Int
qra_n forall a. Num a => a -> a -> a
- Int
3
    ndf' :: Double
ndf' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

    sa :: Double
sa  = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
    sb :: Double
sb  = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
    sc :: Double
sc  = forall a. Floating a => a -> a
sqrt (Double
a33 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

    w :: Double
w   = KBN -> Double
getKBN KBN
qra_w
    x :: Double
x   = KBN -> Double
getKBN KBN
qra_x
    x2 :: Double
x2  = KBN -> Double
getKBN KBN
qra_x2
    x3 :: Double
x3  = KBN -> Double
getKBN KBN
qra_x3
    x4 :: Double
x4  = KBN -> Double
getKBN KBN
qra_x4
    y :: Double
y   = KBN -> Double
getKBN KBN
qra_y
    xy :: Double
xy  = KBN -> Double
getKBN KBN
qra_xy
    x2y :: Double
x2y = KBN -> Double
getKBN KBN
qra_x2y
    y2 :: Double
y2  = KBN -> Double
getKBN KBN
qra_y2

    -- is it useful?
    -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y))

quadRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> QuadRegAcc
quadRegAcc :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double)) -> f a -> QuadRegAcc
quadRegAcc a -> (Double, Double)
f = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\QuadRegAcc
acc a
a -> case a -> (Double, Double)
f a
a of (Double
x,Double
y) -> QuadRegAcc -> Double -> Double -> QuadRegAcc
addQuadReg QuadRegAcc
acc Double
x Double
y) QuadRegAcc
zeroQuadRegAcc

quadRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> QuadRegAcc
quadRegAccW :: forall (f :: * -> *) a.
Foldable f =>
(a -> (Double, Double, Double)) -> f a -> QuadRegAcc
quadRegAccW a -> (Double, Double, Double)
f = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\QuadRegAcc
acc a
a -> case a -> (Double, Double, Double)
f a
a of (Double
x,Double
y,Double
w) -> QuadRegAcc -> Double -> Double -> Double -> QuadRegAcc
addQuadRegW QuadRegAcc
acc Double
x Double
y Double
w) QuadRegAcc
zeroQuadRegAcc

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 1
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with one parameter.
--
-- See 'levenbergMarquardt2' for examples, this is very similar.
--
-- For example we can fit \(f = x \mapsto \beta x + 1\), its derivative is \(\partial_\beta f = x \mapsto x\).
--
-- >>> let scale a (x, y) = (y, a * x + 1, x)
-- >>> PP $ NE.last $ levenbergMarquardt1 scale 1 input2
-- Fit 1.9685 0.04735 4 0.27914
--
-- Not bad, but worse then linear fit which fits the intercept point too.
--
levenbergMarquardt1
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x)\)
    -> Double                                     -- ^ initial parameter, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                   -- ^ non-empty list of iteration results
levenbergMarquardt1 :: forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1 Double -> a -> (Double, Double, Double)
f Double
b0 f a
xs = Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop Double
lambda0 Double
b0 LM1Acc
acc0 where
    acc0 :: LM1Acc
acc0 = Double -> LM1Acc
calcAcc Double
b0

    lambda0 :: Double
lambda0 = Double
c11
      where
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_c11 LM1Acc
acc0

    calcAcc :: Double -> LM1Acc
calcAcc Double
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM1Acc
acc a
p -> case Double -> a -> (Double, Double, Double)
f Double
beta a
p of (Double
y, Double
g, Double
d) -> LM1Acc -> Double -> Double -> Double -> LM1Acc
addLM1Acc LM1Acc
acc Double
y Double
g Double
d) LM1Acc
zeroLM1Acc f a
xs

    loop :: Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop Double
lambda Double
beta LM1Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit Double
beta Double
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) Double
beta LM1Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit Double
beta Double
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) Double
beta' LM1Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: Double
matrix    = forall a. Inv a => a -> a
inv (Double -> Double
lambda1 Double
c11)
        delta :: Double
delta     = forall a b c. Mult a b c => a -> b -> c
mult Double
matrix Double
z1

        beta' :: Double
beta' = forall a. Add a => a -> a -> a
add Double
beta Double
delta
        acc' :: LM1Acc
acc'  = Double -> LM1Acc
calcAcc Double
beta'

        a11 :: Double
a11    = forall a. Inv a => a -> a
inv Double
c11
        errors :: Double
errors = Double
sa

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_wssr LM1Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_wssr LM1Acc
acc'

        ndf :: Int
ndf   = LM1Acc -> Int
lm1_n LM1Acc
acc forall a. Num a => a -> a -> a
- Int
1
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_c11 LM1Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_z1 LM1Acc
acc

-- | 'levenbergMarquardt1' with weights.
levenbergMarquardt1WithWeights
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), w_i\)
    -> Double                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithWeights :: forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1WithWeights Double -> a -> (Double, Double, Double, Double)
f Double
b0 f a
xs = Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop Double
lambda0 Double
b0 LM1Acc
acc0 where
    acc0 :: LM1Acc
acc0 = Double -> LM1Acc
calcAcc Double
b0

    lambda0 :: Double
lambda0 = Double
c11
      where
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_c11 LM1Acc
acc0

    calcAcc :: Double -> LM1Acc
calcAcc Double
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM1Acc
acc a
p -> case Double -> a -> (Double, Double, Double, Double)
f Double
beta a
p of (Double
y, Double
g, Double
d, Double
w) -> LM1Acc -> Double -> Double -> Double -> Double -> LM1Acc
addLM1AccW LM1Acc
acc Double
y Double
g Double
d Double
w) LM1Acc
zeroLM1Acc f a
xs

    loop :: Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop Double
lambda Double
beta LM1Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit Double
beta Double
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) Double
beta LM1Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit Double
beta Double
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> Double -> LM1Acc -> NonEmpty (Fit Double)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) Double
beta' LM1Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: Double
matrix    = forall a. Inv a => a -> a
inv (Double -> Double
lambda1 Double
c11)
        delta :: Double
delta     = forall a b c. Mult a b c => a -> b -> c
mult Double
matrix Double
z1

        beta' :: Double
beta' = forall a. Add a => a -> a -> a
add Double
beta Double
delta
        acc' :: LM1Acc
acc'  = Double -> LM1Acc
calcAcc Double
beta'

        a11 :: Double
a11    = forall a. Inv a => a -> a
inv Double
c11
        errors :: Double
errors = Double
sa

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_wssr LM1Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_wssr LM1Acc
acc'

        ndf :: Int
ndf   = LM1Acc -> Int
lm1_n LM1Acc
acc forall a. Num a => a -> a -> a
- Int
1
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_c11 LM1Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM1Acc -> KBN
lm1_z1 LM1Acc
acc

-- | 'levenbergMarquardt1' with Y-errors.
levenbergMarquardt1WithYerrors
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \delta y_i\)
    -> Double                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithYerrors :: forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1WithYerrors Double -> a -> (Double, Double, Double, Double)
f = forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1WithWeights Double -> a -> (Double, Double, Double, Double)
f' where
    f' :: Double -> a -> (Double, Double, Double, Double)
f' Double
beta a
x = case Double -> a -> (Double, Double, Double, Double)
f Double
beta a
x of (Double
y, Double
fbetax, Double
grad, Double
dy) -> (Double
y, Double
fbetax, Double
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq Double
dy)

-- | 'levenbergMarquardt1' with XY-errors.
levenbergMarquardt1WithXYerrors
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> Double                                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithXYerrors :: forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double, Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1WithXYerrors Double -> a -> (Double, Double, Double, Double, Double, Double)
g = forall (f :: * -> *) a.
Foldable f =>
(Double -> a -> (Double, Double, Double, Double))
-> Double -> f a -> NonEmpty (Fit Double)
levenbergMarquardt1WithWeights Double -> a -> (Double, Double, Double, Double)
g' where
    g' :: Double -> a -> (Double, Double, Double, Double)
g' Double
beta a
x = case Double -> a -> (Double, Double, Double, Double, Double, Double)
g Double
beta a
x of (Double
y, Double
fbetax, Double
grad, Double
f', Double
dx, Double
dy) -> (Double
y, Double
fbetax, Double
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq (Double
f' forall a. Num a => a -> a -> a
* Double
dx) forall a. Num a => a -> a -> a
+ forall a. Num a => a -> a
sq Double
dy)

data LM1Acc = LM1Acc
    { LM1Acc -> Int
lm1_n    :: !Int
    , LM1Acc -> KBN
lm1_c11  :: !KBN
    , LM1Acc -> KBN
lm1_z1   :: !KBN
    , LM1Acc -> KBN
lm1_wssr :: !KBN
    }
  deriving Int -> LM1Acc -> ShowS
[LM1Acc] -> ShowS
LM1Acc -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [LM1Acc] -> ShowS
$cshowList :: [LM1Acc] -> ShowS
show :: LM1Acc -> String
$cshow :: LM1Acc -> String
showsPrec :: Int -> LM1Acc -> ShowS
$cshowsPrec :: Int -> LM1Acc -> ShowS
Show

zeroLM1Acc :: LM1Acc
zeroLM1Acc :: LM1Acc
zeroLM1Acc = Int -> KBN -> KBN -> KBN -> LM1Acc
LM1Acc Int
0 KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN

addLM1Acc :: LM1Acc -> Double -> Double -> Double -> LM1Acc
addLM1Acc :: LM1Acc -> Double -> Double -> Double -> LM1Acc
addLM1Acc LM1Acc {Int
KBN
lm1_wssr :: KBN
lm1_z1 :: KBN
lm1_c11 :: KBN
lm1_n :: Int
lm1_z1 :: LM1Acc -> KBN
lm1_n :: LM1Acc -> Int
lm1_wssr :: LM1Acc -> KBN
lm1_c11 :: LM1Acc -> KBN
..} Double
y Double
f Double
d1 = LM1Acc
    { lm1_n :: Int
lm1_n    = Int
lm1_n forall a. Num a => a -> a -> a
+ Int
1
    , lm1_c11 :: KBN
lm1_c11  = KBN -> Double -> KBN
addKBN KBN
lm1_c11  (Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm1_z1 :: KBN
lm1_z1   = KBN -> Double -> KBN
addKBN KBN
lm1_z1   (Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm1_wssr :: KBN
lm1_wssr = KBN -> Double -> KBN
addKBN KBN
lm1_wssr (Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

addLM1AccW :: LM1Acc -> Double -> Double -> Double -> Double -> LM1Acc
addLM1AccW :: LM1Acc -> Double -> Double -> Double -> Double -> LM1Acc
addLM1AccW LM1Acc {Int
KBN
lm1_wssr :: KBN
lm1_z1 :: KBN
lm1_c11 :: KBN
lm1_n :: Int
lm1_z1 :: LM1Acc -> KBN
lm1_n :: LM1Acc -> Int
lm1_wssr :: LM1Acc -> KBN
lm1_c11 :: LM1Acc -> KBN
..} Double
y Double
f Double
d1 Double
w = LM1Acc
    { lm1_n :: Int
lm1_n    = Int
lm1_n forall a. Num a => a -> a -> a
+ Int
1
    , lm1_c11 :: KBN
lm1_c11  = KBN -> Double -> KBN
addKBN KBN
lm1_c11  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm1_z1 :: KBN
lm1_z1   = KBN -> Double -> KBN
addKBN KBN
lm1_z1   (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm1_wssr :: KBN
lm1_wssr = KBN -> Double -> KBN
addKBN KBN
lm1_wssr (Double
w forall a. Num a => a -> a -> a
* Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 2
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with two parameters.
--
-- You can use this sledgehammer to do a a linear fit:
--
-- >>> let lin (V2 a b) (x, y) = (y, a * x + b, V2 x 1)
--
-- We can then use 'levenbergMarquardt2' to find a fit:
--
-- >>> PP $ levenbergMarquardt2 lin (V2 1 1) input2
-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470
-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681
-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336
-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400
-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645
-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970
-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- This is the same result what 'linearFit' returns:
--
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- == Using AD
--
-- You can use @ad@ to calculate derivatives for you.
--
-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD
--
-- We need a ('Traversable') homogenic triple to represent the two parameters and @x@:
--
-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable)
--
-- Then we define a function @ad@ can operate with:
--
-- >>> let linearF (H3 a b x) = a * x + b
--
-- which we can use to fit the curve in generic way:
--
-- >>> let lin' (V2 a b) (x, y) = case AD.grad' linearF (H3 a b x) of (f, H3 da db _f') -> (y, f, V2 da db)
-- >>> PP $ levenbergMarquardt2 lin' (V2 1 1) input2
-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470
-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681
-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336
-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400
-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645
-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970
-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- == Non-polynomial example
--
-- We can fit other curves too, for example an example from Wikipedia
-- https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm#Example
--
-- >>> let rateF (H3 vmax km s) = (vmax * s) / (km + s)
-- >>> let rateF' (V2 vmax km) (x, y) = case AD.grad' rateF (H3 vmax km x) of (f, H3 vmax' km' _) -> (y, f, V2 vmax' km')
-- >>> let input = zip [0.038,0.194,0.425,0.626,1.253,2.500,3.740] [0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317]
-- >>> PP $ levenbergMarquardt2 rateF' (V2 0.9 0.2) input
-- Fit (V2 0.90000 0.20000) (V2 0.43304 0.43936) 5 1.4455
-- Fit (V2 0.83306 0.25278) (V2 0.39164 0.49729) 5 1.0055
-- Fit (V2 0.59437 0.43508) (V2 0.21158 0.53403) 5 0.18832
-- Fit (V2 0.39687 0.56324) (V2 0.05723 0.25666) 5 0.01062
-- Fit (V2 0.36289 0.56104) (V2 0.04908 0.24007) 5 0.00784
-- Fit (V2 0.36190 0.55662) (V2 0.04887 0.23843) 5 0.00784
-- Fit (V2 0.36184 0.55629) (V2 0.04885 0.23830) 5 0.00784
-- Fit (V2 0.36184 0.55627) (V2 0.04885 0.23829) 5 0.00784
--
-- We get the same result as in the article: 0.362 and 0.556
--
-- The algorithm terminates when a scaling parameter \(\lambda\) becomes larger than 1e20 or smaller than 1e-20, or relative WSSR change is smaller than 1e-10, or sum-of-squared-residuals candidate becomes @NaN@ (i.e. when it would start to produce garbage).
-- You may want to terminate sooner, Numerical Recipes suggest to stop when WSSR decreases by a neglible amount absolutely or fractionally.
--
levenbergMarquardt2
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\)
    -> V2                                 -- ^ initial parameters, \(\beta_0\)
    -> f a                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)               -- ^ non-empty list of iteration results
levenbergMarquardt2 :: forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2)) -> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2 V2 -> a -> (Double, Double, V2)
f V2
b0 f a
xs = Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop Double
lambda0 V2
b0 LM2Acc
acc0 where
    acc0 :: LM2Acc
acc0 = V2 -> LM2Acc
calcAcc V2
b0

    calcAcc :: V2 -> LM2Acc
calcAcc V2
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM2Acc
acc a
p -> case V2 -> a -> (Double, Double, V2)
f V2
beta a
p of (Double
y, Double
g, V2
d) -> LM2Acc -> Double -> Double -> V2 -> LM2Acc
addLM2Acc LM2Acc
acc Double
y Double
g V2
d) LM2Acc
zeroLM2Acc f a
xs

    lambda0 :: Double
lambda0 = forall a. Ord a => a -> a -> a
max Double
l1 Double
l2
      where
        V2 Double
l1 Double
l2 = Double -> Double -> Double -> V2
eigenSM22 Double
c11 Double
c12 Double
c22
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc0
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc0
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c22 LM2Acc
acc0

    loop :: Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop Double
lambda V2
beta LM2Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V2
beta V2
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) V2
beta LM2Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V2
beta V2
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) V2
beta' LM2Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: SM22
matrix    = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> SM22
SM22 (Double -> Double
lambda1 Double
c11) Double
c12 (Double -> Double
lambda1 Double
c22))
        delta :: V2
delta     = forall a b c. Mult a b c => a -> b -> c
mult SM22
matrix (Double -> Double -> V2
V2 Double
z1 Double
z2)

        beta' :: V2
beta' = forall a. Add a => a -> a -> a
add V2
beta V2
delta
        acc' :: LM2Acc
acc'  = V2 -> LM2Acc
calcAcc V2
beta'

        SM22 Double
a11 Double
_ Double
a22 = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> SM22
SM22 Double
c11 Double
c12 Double
c22)
        errors :: V2
errors = Double -> Double -> V2
V2 Double
sa Double
sb

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_wssr LM2Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_wssr LM2Acc
acc'

        ndf :: Int
ndf   = LM2Acc -> Int
lm2_n LM2Acc
acc forall a. Num a => a -> a -> a
- Int
2
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sb :: Double
sb   = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c12 LM2Acc
acc
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c22 LM2Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_z1 LM2Acc
acc
        z2 :: Double
z2  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_z2 LM2Acc
acc

-- | 'levenbergMarquardt2' with weights.
--
-- Because 'levenbergMarquardt2' is an iterative algorithm,
-- not only we can use it to fit curves with known y-errors ('levenbergMarquardt2WithYerrors'),
-- but also with both x and y-errors ('levenbergMarquardt2WithXYerrors').
--
levenbergMarquardt2WithWeights
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\)
    -> V2                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithWeights :: forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2, Double))
-> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2WithWeights V2 -> a -> (Double, Double, V2, Double)
f V2
b0 f a
xs = Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop Double
lambda0 V2
b0 LM2Acc
acc0 where
    acc0 :: LM2Acc
acc0 = V2 -> LM2Acc
calcAcc V2
b0

    lambda0 :: Double
lambda0 = forall a. Ord a => a -> a -> a
max Double
l1 Double
l2
      where
        V2 Double
l1 Double
l2 = Double -> Double -> Double -> V2
eigenSM22 Double
c11 Double
c12 Double
c22
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc0
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc0
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c22 LM2Acc
acc0

    calcAcc :: V2 -> LM2Acc
calcAcc V2
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM2Acc
acc a
p -> case V2 -> a -> (Double, Double, V2, Double)
f V2
beta a
p of (Double
y, Double
g, V2
d, Double
w) -> LM2Acc -> Double -> Double -> V2 -> Double -> LM2Acc
addLM2AccW LM2Acc
acc Double
y Double
g V2
d Double
w) LM2Acc
zeroLM2Acc f a
xs

    loop :: Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop Double
lambda V2
beta LM2Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V2
beta V2
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) V2
beta LM2Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V2
beta V2
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> V2 -> LM2Acc -> NonEmpty (Fit V2)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) V2
beta' LM2Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: SM22
matrix    = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> SM22
SM22 (Double -> Double
lambda1 Double
c11) Double
c12 (Double -> Double
lambda1 Double
c22))
        delta :: V2
delta     = forall a b c. Mult a b c => a -> b -> c
mult SM22
matrix (Double -> Double -> V2
V2 Double
z1 Double
z2)

        beta' :: V2
beta' = forall a. Add a => a -> a -> a
add V2
beta V2
delta
        acc' :: LM2Acc
acc'  = V2 -> LM2Acc
calcAcc V2
beta'

        SM22 Double
a11 Double
_ Double
a22 = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> SM22
SM22 Double
c11 Double
c12 Double
c22)
        errors :: V2
errors = Double -> Double -> V2
V2 Double
sa Double
sb

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_wssr LM2Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_wssr LM2Acc
acc'

        ndf :: Int
ndf   = LM2Acc -> Int
lm2_n LM2Acc
acc forall a. Num a => a -> a -> a
- Int
2
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sb :: Double
sb   = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c11 LM2Acc
acc
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c12 LM2Acc
acc
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_c22 LM2Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_z1 LM2Acc
acc
        z2 :: Double
z2  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM2Acc -> KBN
lm2_z2 LM2Acc
acc

-- | 'levenbergMarquardt2' with Y-errors.
levenbergMarquardt2WithYerrors
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\)
    -> V2                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithYerrors :: forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2, Double))
-> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2WithYerrors V2 -> a -> (Double, Double, V2, Double)
f = forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2, Double))
-> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2WithWeights V2 -> a -> (Double, Double, V2, Double)
f' where
    f' :: V2 -> a -> (Double, Double, V2, Double)
f' V2
beta a
x = case V2 -> a -> (Double, Double, V2, Double)
f V2
beta a
x of (Double
y, Double
fbetax, V2
grad, Double
dy) -> (Double
y, Double
fbetax, V2
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq Double
dy)

-- | 'levenbergMarquardt2' with XY-errors.
levenbergMarquardt2WithXYerrors
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> V2                                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithXYerrors :: forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2, Double, Double, Double))
-> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2WithXYerrors V2 -> a -> (Double, Double, V2, Double, Double, Double)
g = forall (f :: * -> *) a.
Foldable f =>
(V2 -> a -> (Double, Double, V2, Double))
-> V2 -> f a -> NonEmpty (Fit V2)
levenbergMarquardt2WithWeights V2 -> a -> (Double, Double, V2, Double)
g' where
    g' :: V2 -> a -> (Double, Double, V2, Double)
g' V2
beta a
x = case V2 -> a -> (Double, Double, V2, Double, Double, Double)
g V2
beta a
x of (Double
y, Double
fbetax, V2
grad, Double
f', Double
dx, Double
dy) -> (Double
y, Double
fbetax, V2
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq (Double
f' forall a. Num a => a -> a -> a
* Double
dx) forall a. Num a => a -> a -> a
+ forall a. Num a => a -> a
sq Double
dy)

data LM2Acc = LM2Acc
    { LM2Acc -> Int
lm2_n    :: !Int
    , LM2Acc -> KBN
lm2_c11  :: !KBN
    , LM2Acc -> KBN
lm2_c12  :: !KBN
    , LM2Acc -> KBN
lm2_c22  :: !KBN
    , LM2Acc -> KBN
lm2_z1   :: !KBN
    , LM2Acc -> KBN
lm2_z2   :: !KBN
    , LM2Acc -> KBN
lm2_wssr :: !KBN
    }
  deriving Int -> LM2Acc -> ShowS
[LM2Acc] -> ShowS
LM2Acc -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [LM2Acc] -> ShowS
$cshowList :: [LM2Acc] -> ShowS
show :: LM2Acc -> String
$cshow :: LM2Acc -> String
showsPrec :: Int -> LM2Acc -> ShowS
$cshowsPrec :: Int -> LM2Acc -> ShowS
Show

zeroLM2Acc :: LM2Acc
zeroLM2Acc :: LM2Acc
zeroLM2Acc = Int -> KBN -> KBN -> KBN -> KBN -> KBN -> KBN -> LM2Acc
LM2Acc Int
0 KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN

addLM2Acc :: LM2Acc -> Double -> Double -> V2 -> LM2Acc
addLM2Acc :: LM2Acc -> Double -> Double -> V2 -> LM2Acc
addLM2Acc LM2Acc {Int
KBN
lm2_wssr :: KBN
lm2_z2 :: KBN
lm2_z1 :: KBN
lm2_c22 :: KBN
lm2_c12 :: KBN
lm2_c11 :: KBN
lm2_n :: Int
lm2_z2 :: LM2Acc -> KBN
lm2_z1 :: LM2Acc -> KBN
lm2_c12 :: LM2Acc -> KBN
lm2_n :: LM2Acc -> Int
lm2_wssr :: LM2Acc -> KBN
lm2_c22 :: LM2Acc -> KBN
lm2_c11 :: LM2Acc -> KBN
..} Double
y Double
f (V2 Double
d1 Double
d2) = LM2Acc
    { lm2_n :: Int
lm2_n    = Int
lm2_n forall a. Num a => a -> a -> a
+ Int
1
    , lm2_c11 :: KBN
lm2_c11  = KBN -> Double -> KBN
addKBN KBN
lm2_c11  (Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm2_c12 :: KBN
lm2_c12  = KBN -> Double -> KBN
addKBN KBN
lm2_c12  (Double
d1 forall a. Num a => a -> a -> a
* Double
d2)
    , lm2_c22 :: KBN
lm2_c22  = KBN -> Double -> KBN
addKBN KBN
lm2_c22  (Double
d2 forall a. Num a => a -> a -> a
* Double
d2)
    , lm2_z1 :: KBN
lm2_z1   = KBN -> Double -> KBN
addKBN KBN
lm2_z1   (Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm2_z2 :: KBN
lm2_z2   = KBN -> Double -> KBN
addKBN KBN
lm2_z2   (Double
d2 forall a. Num a => a -> a -> a
* Double
res)
    , lm2_wssr :: KBN
lm2_wssr = KBN -> Double -> KBN
addKBN KBN
lm2_wssr (Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

addLM2AccW :: LM2Acc -> Double -> Double -> V2 -> Double -> LM2Acc
addLM2AccW :: LM2Acc -> Double -> Double -> V2 -> Double -> LM2Acc
addLM2AccW LM2Acc {Int
KBN
lm2_wssr :: KBN
lm2_z2 :: KBN
lm2_z1 :: KBN
lm2_c22 :: KBN
lm2_c12 :: KBN
lm2_c11 :: KBN
lm2_n :: Int
lm2_z2 :: LM2Acc -> KBN
lm2_z1 :: LM2Acc -> KBN
lm2_c12 :: LM2Acc -> KBN
lm2_n :: LM2Acc -> Int
lm2_wssr :: LM2Acc -> KBN
lm2_c22 :: LM2Acc -> KBN
lm2_c11 :: LM2Acc -> KBN
..} Double
y Double
f (V2 Double
d1 Double
d2) Double
w = LM2Acc
    { lm2_n :: Int
lm2_n    = Int
lm2_n forall a. Num a => a -> a -> a
+ Int
1
    , lm2_c11 :: KBN
lm2_c11  = KBN -> Double -> KBN
addKBN KBN
lm2_c11  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm2_c12 :: KBN
lm2_c12  = KBN -> Double -> KBN
addKBN KBN
lm2_c12  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d2)
    , lm2_c22 :: KBN
lm2_c22  = KBN -> Double -> KBN
addKBN KBN
lm2_c22  (Double
w forall a. Num a => a -> a -> a
* Double
d2 forall a. Num a => a -> a -> a
* Double
d2)
    , lm2_z1 :: KBN
lm2_z1   = KBN -> Double -> KBN
addKBN KBN
lm2_z1   (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm2_z2 :: KBN
lm2_z2   = KBN -> Double -> KBN
addKBN KBN
lm2_z2   (Double
w forall a. Num a => a -> a -> a
* Double
d2 forall a. Num a => a -> a -> a
* Double
res)
    , lm2_wssr :: KBN
lm2_wssr = KBN -> Double -> KBN
addKBN KBN
lm2_wssr (Double
w forall a. Num a => a -> a -> a
* Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 3
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with three parameters.
--
-- See 'levenbergMarquardt2' for examples, this is very similar.
--
-- >>> let quad (V3 a b c) (x, y) = (y, a * x * x + b * x + c, V3 (x * x) x 1)
-- >>> PP $ NE.last $ levenbergMarquardt3 quad (V3 2 2 2) input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
-- Same as quadratic fit, just less direct:
--
-- >>> PP $ quadraticFit id input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
levenbergMarquardt3
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\)
    -> V3                                 -- ^ initial parameters, \(\beta_0\)
    -> f a                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)               -- ^ non-empty list of iteration results
levenbergMarquardt3 :: forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3)) -> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3 V3 -> a -> (Double, Double, V3)
f V3
b0 f a
xs = Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop Double
lambda0 V3
b0 LM3Acc
acc0 where
    acc0 :: LM3Acc
acc0 = V3 -> LM3Acc
calcAcc V3
b0

    calcAcc :: V3 -> LM3Acc
calcAcc V3
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM3Acc
acc a
p -> case V3 -> a -> (Double, Double, V3)
f V3
beta a
p of (Double
y, Double
g, V3
d) -> LM3Acc -> Double -> Double -> V3 -> LM3Acc
addLM3Acc LM3Acc
acc Double
y Double
g V3
d) LM3Acc
zeroLM3Acc f a
xs

    -- frobenius norm is larger than largest eigen value.
    -- calculating the eigen values for 3x3 (symmetric) matrix is becoming complicated.
    lambda0 :: Double
lambda0 = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *). Foldable f => f Double -> Double
sumKBN [ forall a. Num a => a -> a
sq Double
c11
                            , Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c12, forall a. Num a => a -> a
sq Double
c22
                            , Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c13, Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c23, forall a. Num a => a -> a
sq Double
c33]
      where
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c11 LM3Acc
acc0
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c12 LM3Acc
acc0
        c13 :: Double
c13 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c13 LM3Acc
acc0
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c22 LM3Acc
acc0
        c23 :: Double
c23 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c23 LM3Acc
acc0
        c33 :: Double
c33 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c33 LM3Acc
acc0

    loop :: Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop Double
lambda V3
beta LM3Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V3
beta V3
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) V3
beta LM3Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V3
beta V3
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) V3
beta' LM3Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: SM33
matrix    = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> Double -> Double -> Double -> SM33
SM33 (Double -> Double
lambda1 Double
c11)
                               Double
c12          (Double -> Double
lambda1 Double
c22)
                               Double
c13          Double
c23           (Double -> Double
lambda1 Double
c33))
        delta :: V3
delta     = forall a b c. Mult a b c => a -> b -> c
mult SM33
matrix (Double -> Double -> Double -> V3
V3 Double
z1 Double
z2 Double
z3)

        beta' :: V3
beta' = forall a. Add a => a -> a -> a
add V3
beta V3
delta
        acc' :: LM3Acc
acc'  = V3 -> LM3Acc
calcAcc V3
beta'

        SM33 Double
a11
             Double
_  Double
a22
             Double
_  Double
_   Double
a33 = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> Double -> Double -> Double -> SM33
SM33 Double
c11
                                    Double
c12 Double
c22
                                    Double
c13 Double
c23 Double
c33)
        errors :: V3
errors = Double -> Double -> Double -> V3
V3 Double
sa Double
sb Double
sc

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_wssr LM3Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_wssr LM3Acc
acc'

        ndf :: Int
ndf   = LM3Acc -> Int
lm3_n LM3Acc
acc forall a. Num a => a -> a -> a
- Int
3
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sb :: Double
sb   = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sc :: Double
sc   = forall a. Floating a => a -> a
sqrt (Double
a33 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c11 LM3Acc
acc
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c12 LM3Acc
acc
        c13 :: Double
c13 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c13 LM3Acc
acc
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c22 LM3Acc
acc
        c23 :: Double
c23 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c23 LM3Acc
acc
        c33 :: Double
c33 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c33 LM3Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z1 LM3Acc
acc
        z2 :: Double
z2  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z2 LM3Acc
acc
        z3 :: Double
z3  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z3 LM3Acc
acc

-- | 'levenbergMarquardt3' with weights.
levenbergMarquardt3WithWeights
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\)
    -> V3                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithWeights :: forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3, Double))
-> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3WithWeights V3 -> a -> (Double, Double, V3, Double)
f V3
b0 f a
xs = Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop Double
lambda0 V3
b0 LM3Acc
acc0 where
    acc0 :: LM3Acc
acc0 = V3 -> LM3Acc
calcAcc V3
b0

    lambda0 :: Double
lambda0 = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *). Foldable f => f Double -> Double
sumKBN [ forall a. Num a => a -> a
sq Double
c11
                            , Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c12, forall a. Num a => a -> a
sq Double
c22
                            , Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c13, Double
2 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
c23, forall a. Num a => a -> a
sq Double
c33]
      where
        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c11 LM3Acc
acc0
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c12 LM3Acc
acc0
        c13 :: Double
c13 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c13 LM3Acc
acc0
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c22 LM3Acc
acc0
        c23 :: Double
c23 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c23 LM3Acc
acc0
        c33 :: Double
c33 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c33 LM3Acc
acc0

    calcAcc :: V3 -> LM3Acc
calcAcc V3
beta = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\LM3Acc
acc a
p -> case V3 -> a -> (Double, Double, V3, Double)
f V3
beta a
p of (Double
y, Double
g, V3
d, Double
w) -> LM3Acc -> Double -> Double -> V3 -> Double -> LM3Acc
addLM3AccW LM3Acc
acc Double
y Double
g V3
d Double
w) LM3Acc
zeroLM3Acc f a
xs

    loop :: Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop Double
lambda V3
beta LM3Acc
acc
        | Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr'
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V3
beta V3
errors Int
ndf Double
wssr forall a. a -> [a] -> NonEmpty a
NE.:| []

        | Double
wssr' forall a. Ord a => a -> a -> Bool
>= Double
wssr
        = Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop (Double
lambda forall a. Num a => a -> a -> a
* Double
10) V3
beta LM3Acc
acc

        | Bool
otherwise
        = forall v. v -> v -> Int -> Double -> Fit v
Fit V3
beta V3
errors Int
ndf Double
wssr forall a. a -> NonEmpty a -> NonEmpty a
`NE.cons` Double -> V3 -> LM3Acc -> NonEmpty (Fit V3)
loop (Double
lambda forall a. Fractional a => a -> a -> a
/ Double
10) V3
beta' LM3Acc
acc'

      where
        lambda1 :: Double -> Double
lambda1 Double
z = (Double
1 forall a. Num a => a -> a -> a
+ Double
lambda) forall a. Num a => a -> a -> a
* Double
z
        matrix :: SM33
matrix    = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> Double -> Double -> Double -> SM33
SM33 (Double -> Double
lambda1 Double
c11)
                               Double
c12          (Double -> Double
lambda1 Double
c22)
                               Double
c13          Double
c23           (Double -> Double
lambda1 Double
c33))
        delta :: V3
delta     = forall a b c. Mult a b c => a -> b -> c
mult SM33
matrix (Double -> Double -> Double -> V3
V3 Double
z1 Double
z2 Double
z3)

        beta' :: V3
beta' = forall a. Add a => a -> a -> a
add V3
beta V3
delta
        acc' :: LM3Acc
acc'  = V3 -> LM3Acc
calcAcc V3
beta'

        SM33 Double
a11 Double
_ Double
_ Double
a22 Double
_ Double
a33 = forall a. Inv a => a -> a
inv (Double -> Double -> Double -> Double -> Double -> Double -> SM33
SM33 Double
c11 Double
c12 Double
c13 Double
c22 Double
c23 Double
c33)
        errors :: V3
errors = Double -> Double -> Double -> V3
V3 Double
sa Double
sb Double
sc

        wssr :: Double
wssr  = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_wssr LM3Acc
acc
        wssr' :: Double
wssr' =         KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_wssr LM3Acc
acc'

        ndf :: Int
ndf   = LM3Acc -> Int
lm3_n LM3Acc
acc forall a. Num a => a -> a -> a
- Int
3
        ndf' :: Double
ndf'  = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ndf :: Double

        sa :: Double
sa   = forall a. Floating a => a -> a
sqrt (Double
a11 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sb :: Double
sb   = forall a. Floating a => a -> a
sqrt (Double
a22 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')
        sc :: Double
sc   = forall a. Floating a => a -> a
sqrt (Double
a33 forall a. Num a => a -> a -> a
* Double
wssr forall a. Fractional a => a -> a -> a
/ Double
ndf')

        c11 :: Double
c11 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c11 LM3Acc
acc
        c12 :: Double
c12 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c12 LM3Acc
acc
        c13 :: Double
c13 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c13 LM3Acc
acc
        c22 :: Double
c22 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c22 LM3Acc
acc
        c23 :: Double
c23 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c23 LM3Acc
acc
        c33 :: Double
c33 = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_c33 LM3Acc
acc
        z1 :: Double
z1  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z1 LM3Acc
acc
        z2 :: Double
z2  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z2 LM3Acc
acc
        z3 :: Double
z3  = KBN -> Double
getKBN forall a b. (a -> b) -> a -> b
$ LM3Acc -> KBN
lm3_z3 LM3Acc
acc

-- | 'levenbergMarquardt3' with Y-errors.
levenbergMarquardt3WithYerrors
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\)
    -> V3                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithYerrors :: forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3, Double))
-> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3WithYerrors V3 -> a -> (Double, Double, V3, Double)
f = forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3, Double))
-> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3WithWeights V3 -> a -> (Double, Double, V3, Double)
f' where
    f' :: V3 -> a -> (Double, Double, V3, Double)
f' V3
beta a
x = case V3 -> a -> (Double, Double, V3, Double)
f V3
beta a
x of (Double
y, Double
fbetax, V3
grad, Double
dy) -> (Double
y, Double
fbetax, V3
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq Double
dy)

-- | 'levenbergMarquardt3' with XY-errors.
levenbergMarquardt3WithXYerrors
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> V3                                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithXYerrors :: forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3, Double, Double, Double))
-> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3WithXYerrors V3 -> a -> (Double, Double, V3, Double, Double, Double)
g = forall (f :: * -> *) a.
Foldable f =>
(V3 -> a -> (Double, Double, V3, Double))
-> V3 -> f a -> NonEmpty (Fit V3)
levenbergMarquardt3WithWeights V3 -> a -> (Double, Double, V3, Double)
g' where
    g' :: V3 -> a -> (Double, Double, V3, Double)
g' V3
beta a
x = case V3 -> a -> (Double, Double, V3, Double, Double, Double)
g V3
beta a
x of (Double
y, Double
fbetax, V3
grad, Double
f', Double
dx, Double
dy) -> (Double
y, Double
fbetax, V3
grad, forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a
sq (Double
f' forall a. Num a => a -> a -> a
* Double
dx) forall a. Num a => a -> a -> a
+ forall a. Num a => a -> a
sq Double
dy)

data LM3Acc = LM3Acc
    { LM3Acc -> Int
lm3_n    :: !Int
    , LM3Acc -> KBN
lm3_c11  :: !KBN
    , LM3Acc -> KBN
lm3_c12  :: !KBN
    , LM3Acc -> KBN
lm3_c13  :: !KBN
    , LM3Acc -> KBN
lm3_c22  :: !KBN
    , LM3Acc -> KBN
lm3_c23  :: !KBN
    , LM3Acc -> KBN
lm3_c33  :: !KBN
    , LM3Acc -> KBN
lm3_z1   :: !KBN
    , LM3Acc -> KBN
lm3_z2   :: !KBN
    , LM3Acc -> KBN
lm3_z3   :: !KBN
    , LM3Acc -> KBN
lm3_wssr :: !KBN
    }
  deriving Int -> LM3Acc -> ShowS
[LM3Acc] -> ShowS
LM3Acc -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [LM3Acc] -> ShowS
$cshowList :: [LM3Acc] -> ShowS
show :: LM3Acc -> String
$cshow :: LM3Acc -> String
showsPrec :: Int -> LM3Acc -> ShowS
$cshowsPrec :: Int -> LM3Acc -> ShowS
Show

zeroLM3Acc :: LM3Acc
zeroLM3Acc :: LM3Acc
zeroLM3Acc = Int
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> LM3Acc
LM3Acc Int
0 KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN

addLM3Acc :: LM3Acc -> Double -> Double -> V3 -> LM3Acc
addLM3Acc :: LM3Acc -> Double -> Double -> V3 -> LM3Acc
addLM3Acc LM3Acc {Int
KBN
lm3_wssr :: KBN
lm3_z3 :: KBN
lm3_z2 :: KBN
lm3_z1 :: KBN
lm3_c33 :: KBN
lm3_c23 :: KBN
lm3_c22 :: KBN
lm3_c13 :: KBN
lm3_c12 :: KBN
lm3_c11 :: KBN
lm3_n :: Int
lm3_z3 :: LM3Acc -> KBN
lm3_z2 :: LM3Acc -> KBN
lm3_z1 :: LM3Acc -> KBN
lm3_n :: LM3Acc -> Int
lm3_wssr :: LM3Acc -> KBN
lm3_c33 :: LM3Acc -> KBN
lm3_c23 :: LM3Acc -> KBN
lm3_c22 :: LM3Acc -> KBN
lm3_c13 :: LM3Acc -> KBN
lm3_c12 :: LM3Acc -> KBN
lm3_c11 :: LM3Acc -> KBN
..} Double
y Double
f (V3 Double
d1 Double
d2 Double
d3) = LM3Acc
    { lm3_n :: Int
lm3_n    = Int
lm3_n forall a. Num a => a -> a -> a
+ Int
1
    , lm3_c11 :: KBN
lm3_c11  = KBN -> Double -> KBN
addKBN KBN
lm3_c11  (Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm3_c12 :: KBN
lm3_c12  = KBN -> Double -> KBN
addKBN KBN
lm3_c12  (Double
d1 forall a. Num a => a -> a -> a
* Double
d2)
    , lm3_c13 :: KBN
lm3_c13  = KBN -> Double -> KBN
addKBN KBN
lm3_c12  (Double
d1 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_c22 :: KBN
lm3_c22  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
d2 forall a. Num a => a -> a -> a
* Double
d2)
    , lm3_c23 :: KBN
lm3_c23  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
d2 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_c33 :: KBN
lm3_c33  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
d3 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_z1 :: KBN
lm3_z1   = KBN -> Double -> KBN
addKBN KBN
lm3_z1   (Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_z2 :: KBN
lm3_z2   = KBN -> Double -> KBN
addKBN KBN
lm3_z2   (Double
d2 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_z3 :: KBN
lm3_z3   = KBN -> Double -> KBN
addKBN KBN
lm3_z3   (Double
d3 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_wssr :: KBN
lm3_wssr = KBN -> Double -> KBN
addKBN KBN
lm3_wssr (Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

addLM3AccW :: LM3Acc -> Double -> Double -> V3 -> Double -> LM3Acc
addLM3AccW :: LM3Acc -> Double -> Double -> V3 -> Double -> LM3Acc
addLM3AccW LM3Acc {Int
KBN
lm3_wssr :: KBN
lm3_z3 :: KBN
lm3_z2 :: KBN
lm3_z1 :: KBN
lm3_c33 :: KBN
lm3_c23 :: KBN
lm3_c22 :: KBN
lm3_c13 :: KBN
lm3_c12 :: KBN
lm3_c11 :: KBN
lm3_n :: Int
lm3_z3 :: LM3Acc -> KBN
lm3_z2 :: LM3Acc -> KBN
lm3_z1 :: LM3Acc -> KBN
lm3_n :: LM3Acc -> Int
lm3_wssr :: LM3Acc -> KBN
lm3_c33 :: LM3Acc -> KBN
lm3_c23 :: LM3Acc -> KBN
lm3_c22 :: LM3Acc -> KBN
lm3_c13 :: LM3Acc -> KBN
lm3_c12 :: LM3Acc -> KBN
lm3_c11 :: LM3Acc -> KBN
..} Double
y Double
f (V3 Double
d1 Double
d2 Double
d3) Double
w = LM3Acc
    { lm3_n :: Int
lm3_n    = Int
lm3_n forall a. Num a => a -> a -> a
+ Int
1
    , lm3_c11 :: KBN
lm3_c11  = KBN -> Double -> KBN
addKBN KBN
lm3_c11  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d1)
    , lm3_c12 :: KBN
lm3_c12  = KBN -> Double -> KBN
addKBN KBN
lm3_c12  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d2)
    , lm3_c13 :: KBN
lm3_c13  = KBN -> Double -> KBN
addKBN KBN
lm3_c12  (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_c22 :: KBN
lm3_c22  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
w forall a. Num a => a -> a -> a
* Double
d2 forall a. Num a => a -> a -> a
* Double
d2)
    , lm3_c23 :: KBN
lm3_c23  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
w forall a. Num a => a -> a -> a
* Double
d2 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_c33 :: KBN
lm3_c33  = KBN -> Double -> KBN
addKBN KBN
lm3_c22  (Double
w forall a. Num a => a -> a -> a
* Double
d3 forall a. Num a => a -> a -> a
* Double
d3)
    , lm3_z1 :: KBN
lm3_z1   = KBN -> Double -> KBN
addKBN KBN
lm3_z1   (Double
w forall a. Num a => a -> a -> a
* Double
d1 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_z2 :: KBN
lm3_z2   = KBN -> Double -> KBN
addKBN KBN
lm3_z2   (Double
w forall a. Num a => a -> a -> a
* Double
d2 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_z3 :: KBN
lm3_z3   = KBN -> Double -> KBN
addKBN KBN
lm3_z3   (Double
w forall a. Num a => a -> a -> a
* Double
d3 forall a. Num a => a -> a -> a
* Double
res)
    , lm3_wssr :: KBN
lm3_wssr = KBN -> Double -> KBN
addKBN KBN
lm3_wssr (Double
w forall a. Num a => a -> a -> a
* Double
res forall a. Num a => a -> a -> a
* Double
res)
    }
  where
    res :: Double
res = Double
y forall a. Num a => a -> a -> a
- Double
f

-------------------------------------------------------------------------------
-- Output
-------------------------------------------------------------------------------

-- | Result of a curve fit.
data Fit v = Fit
    { forall v. Fit v -> v
fitParams :: !v       -- ^ fit parameters
    , forall v. Fit v -> v
fitErrors :: !v       -- ^ asympotic standard errors, /assuming a good fit/
    , forall v. Fit v -> Int
fitNDF    :: !Int     -- ^ number of degrees of freedom
    , forall v. Fit v -> Double
fitWSSR   :: !Double  -- ^ sum of squares of residuals
    }
  deriving Int -> Fit v -> ShowS
forall v. Show v => Int -> Fit v -> ShowS
forall v. Show v => [Fit v] -> ShowS
forall v. Show v => Fit v -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Fit v] -> ShowS
$cshowList :: forall v. Show v => [Fit v] -> ShowS
show :: Fit v -> String
$cshow :: forall v. Show v => Fit v -> String
showsPrec :: Int -> Fit v -> ShowS
$cshowsPrec :: forall v. Show v => Int -> Fit v -> ShowS
Show

instance (NFData v) => NFData (Fit v) where
    rnf :: Fit v -> ()
rnf (Fit v
p v
e Int
_ Double
_) = forall a. NFData a => a -> ()
rnf v
p seq :: forall a b. a -> b -> b
`seq` forall a. NFData a => a -> ()
rnf v
e

-------------------------------------------------------------------------------
-- LinRegAcc
-------------------------------------------------------------------------------

-- | Linear regression accumulator.
data LinRegAcc = LinRegAcc
    { LinRegAcc -> Int
lra_n  :: {-# UNPACK #-} !Int  -- ^ \(n\)
    , LinRegAcc -> KBN
lra_w  :: {-# UNPACK #-} !KBN  -- ^ \(\sum w_i\)
    , LinRegAcc -> KBN
lra_x  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i \)
    , LinRegAcc -> KBN
lra_x2 :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 \)
    , LinRegAcc -> KBN
lra_y  :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i \)
    , LinRegAcc -> KBN
lra_xy :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i y_i \)
    , LinRegAcc -> KBN
lra_y2 :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i^2 \)
    }
  deriving Int -> LinRegAcc -> ShowS
[LinRegAcc] -> ShowS
LinRegAcc -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [LinRegAcc] -> ShowS
$cshowList :: [LinRegAcc] -> ShowS
show :: LinRegAcc -> String
$cshow :: LinRegAcc -> String
showsPrec :: Int -> LinRegAcc -> ShowS
$cshowsPrec :: Int -> LinRegAcc -> ShowS
Show

instance NFData LinRegAcc where
    rnf :: LinRegAcc -> ()
rnf LinRegAcc {} = ()

-- | All-zeroes 'LinRegAcc'.
zeroLinRegAcc :: LinRegAcc
zeroLinRegAcc :: LinRegAcc
zeroLinRegAcc = Int -> KBN -> KBN -> KBN -> KBN -> KBN -> KBN -> LinRegAcc
LinRegAcc Int
0 KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN

-- | Add a point to linreg accumulator.
addLinReg
    :: LinRegAcc
    -> Double   -- ^ x
    -> Double   -- ^ y
    -> LinRegAcc
addLinReg :: LinRegAcc -> Double -> Double -> LinRegAcc
addLinReg LinRegAcc {Int
KBN
lra_y2 :: KBN
lra_xy :: KBN
lra_y :: KBN
lra_x2 :: KBN
lra_x :: KBN
lra_w :: KBN
lra_n :: Int
lra_y2 :: LinRegAcc -> KBN
lra_xy :: LinRegAcc -> KBN
lra_y :: LinRegAcc -> KBN
lra_x2 :: LinRegAcc -> KBN
lra_x :: LinRegAcc -> KBN
lra_w :: LinRegAcc -> KBN
lra_n :: LinRegAcc -> Int
..} Double
x Double
y = LinRegAcc
    { lra_n :: Int
lra_n  = Int
lra_n forall a. Num a => a -> a -> a
+ Int
1
    , lra_w :: KBN
lra_w  = KBN -> Double -> KBN
addKBN KBN
lra_w  Double
1
    , lra_x :: KBN
lra_x  = KBN -> Double -> KBN
addKBN KBN
lra_x  Double
x
    , lra_x2 :: KBN
lra_x2 = KBN -> Double -> KBN
addKBN KBN
lra_x2 (Double
x forall a. Num a => a -> a -> a
* Double
x)
    , lra_y :: KBN
lra_y  = KBN -> Double -> KBN
addKBN KBN
lra_y  Double
y
    , lra_xy :: KBN
lra_xy = KBN -> Double -> KBN
addKBN KBN
lra_xy (Double
x forall a. Num a => a -> a -> a
* Double
y)
    , lra_y2 :: KBN
lra_y2 = KBN -> Double -> KBN
addKBN KBN
lra_y2 (Double
y forall a. Num a => a -> a -> a
* Double
y)
    }

-- | Add a weighted point to linreg accumulator.
addLinRegW
    :: LinRegAcc
    -> Double   -- ^ x
    -> Double   -- ^ y
    -> Double   -- ^ w
    -> LinRegAcc
addLinRegW :: LinRegAcc -> Double -> Double -> Double -> LinRegAcc
addLinRegW LinRegAcc {Int
KBN
lra_y2 :: KBN
lra_xy :: KBN
lra_y :: KBN
lra_x2 :: KBN
lra_x :: KBN
lra_w :: KBN
lra_n :: Int
lra_y2 :: LinRegAcc -> KBN
lra_xy :: LinRegAcc -> KBN
lra_y :: LinRegAcc -> KBN
lra_x2 :: LinRegAcc -> KBN
lra_x :: LinRegAcc -> KBN
lra_w :: LinRegAcc -> KBN
lra_n :: LinRegAcc -> Int
..} Double
x Double
y Double
w = LinRegAcc
    { lra_n :: Int
lra_n  = Int
lra_n forall a. Num a => a -> a -> a
+ Int
1
    , lra_w :: KBN
lra_w  = KBN -> Double -> KBN
addKBN KBN
lra_w  Double
w
    , lra_x :: KBN
lra_x  = KBN -> Double -> KBN
addKBN KBN
lra_x  (Double
w forall a. Num a => a -> a -> a
* Double
x)
    , lra_x2 :: KBN
lra_x2 = KBN -> Double -> KBN
addKBN KBN
lra_x2 (Double
w forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
* Double
x)
    , lra_y :: KBN
lra_y  = KBN -> Double -> KBN
addKBN KBN
lra_y  (Double
w forall a. Num a => a -> a -> a
* Double
y)
    , lra_xy :: KBN
lra_xy = KBN -> Double -> KBN
addKBN KBN
lra_xy (Double
w forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
* Double
y)
    , lra_y2 :: KBN
lra_y2 = KBN -> Double -> KBN
addKBN KBN
lra_y2 (Double
w forall a. Num a => a -> a -> a
* Double
y forall a. Num a => a -> a -> a
* Double
y)
    }

-------------------------------------------------------------------------------
-- QuadRegAcc
-------------------------------------------------------------------------------

-- | Quadratic regression accumulator.
data QuadRegAcc = QuadRegAcc
    { QuadRegAcc -> Int
qra_n   :: {-# UNPACK #-} !Int  -- ^ \(n\)
    , QuadRegAcc -> KBN
qra_w   :: {-# UNPACK #-} !KBN  -- ^ \(\sum w_i\)
    , QuadRegAcc -> KBN
qra_x   :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i \)
    , QuadRegAcc -> KBN
qra_x2  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 \)
    , QuadRegAcc -> KBN
qra_x3  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^3 \)
    , QuadRegAcc -> KBN
qra_x4  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^4 \)
    , QuadRegAcc -> KBN
qra_y   :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i \)
    , QuadRegAcc -> KBN
qra_xy  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i y_i \)
    , QuadRegAcc -> KBN
qra_x2y :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 y_i \)
    , QuadRegAcc -> KBN
qra_y2  :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i^2 \)
    }
  deriving Int -> QuadRegAcc -> ShowS
[QuadRegAcc] -> ShowS
QuadRegAcc -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [QuadRegAcc] -> ShowS
$cshowList :: [QuadRegAcc] -> ShowS
show :: QuadRegAcc -> String
$cshow :: QuadRegAcc -> String
showsPrec :: Int -> QuadRegAcc -> ShowS
$cshowsPrec :: Int -> QuadRegAcc -> ShowS
Show

instance NFData QuadRegAcc where
    rnf :: QuadRegAcc -> ()
rnf QuadRegAcc {} = ()

-- | All-zeroes 'QuadRegAcc'.
zeroQuadRegAcc :: QuadRegAcc
zeroQuadRegAcc :: QuadRegAcc
zeroQuadRegAcc = Int
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> KBN
-> QuadRegAcc
QuadRegAcc Int
0 KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN KBN
zeroKBN

-- | Add a point to quadreg accumulator.
addQuadReg
    :: QuadRegAcc
    -> Double  -- ^ x
    -> Double  -- ^ y
    -> QuadRegAcc
addQuadReg :: QuadRegAcc -> Double -> Double -> QuadRegAcc
addQuadReg QuadRegAcc {Int
KBN
qra_y2 :: KBN
qra_x2y :: KBN
qra_xy :: KBN
qra_y :: KBN
qra_x4 :: KBN
qra_x3 :: KBN
qra_x2 :: KBN
qra_x :: KBN
qra_w :: KBN
qra_n :: Int
qra_y2 :: QuadRegAcc -> KBN
qra_x2y :: QuadRegAcc -> KBN
qra_xy :: QuadRegAcc -> KBN
qra_y :: QuadRegAcc -> KBN
qra_x4 :: QuadRegAcc -> KBN
qra_x3 :: QuadRegAcc -> KBN
qra_x2 :: QuadRegAcc -> KBN
qra_x :: QuadRegAcc -> KBN
qra_w :: QuadRegAcc -> KBN
qra_n :: QuadRegAcc -> Int
..} Double
x Double
y = QuadRegAcc
    { qra_n :: Int
qra_n    = Int
qra_n forall a. Num a => a -> a -> a
+ Int
1
    , qra_w :: KBN
qra_w    = KBN -> Double -> KBN
addKBN KBN
qra_w   Double
1
    , qra_x :: KBN
qra_x    = KBN -> Double -> KBN
addKBN KBN
qra_x   Double
x
    , qra_x2 :: KBN
qra_x2   = KBN -> Double -> KBN
addKBN KBN
qra_x2  Double
x2
    , qra_x3 :: KBN
qra_x3   = KBN -> Double -> KBN
addKBN KBN
qra_x3  (Double
x forall a. Num a => a -> a -> a
* Double
x2)
    , qra_x4 :: KBN
qra_x4   = KBN -> Double -> KBN
addKBN KBN
qra_x4  (Double
x2 forall a. Num a => a -> a -> a
* Double
x2)
    , qra_y :: KBN
qra_y    = KBN -> Double -> KBN
addKBN KBN
qra_y   Double
y
    , qra_xy :: KBN
qra_xy   = KBN -> Double -> KBN
addKBN KBN
qra_xy  (Double
x forall a. Num a => a -> a -> a
* Double
y)
    , qra_x2y :: KBN
qra_x2y  = KBN -> Double -> KBN
addKBN KBN
qra_x2y (Double
x2 forall a. Num a => a -> a -> a
* Double
y)
    , qra_y2 :: KBN
qra_y2   = KBN -> Double -> KBN
addKBN KBN
qra_y2  (Double
y forall a. Num a => a -> a -> a
* Double
y)
    }
  where
    x2 :: Double
x2 = Double
x forall a. Num a => a -> a -> a
* Double
x

-- | Add a weighted point to quadreg accumulator.
addQuadRegW
    :: QuadRegAcc
    -> Double  -- ^ x
    -> Double  -- ^ y
    -> Double  -- ^ w
    -> QuadRegAcc
addQuadRegW :: QuadRegAcc -> Double -> Double -> Double -> QuadRegAcc
addQuadRegW QuadRegAcc {Int
KBN
qra_y2 :: KBN
qra_x2y :: KBN
qra_xy :: KBN
qra_y :: KBN
qra_x4 :: KBN
qra_x3 :: KBN
qra_x2 :: KBN
qra_x :: KBN
qra_w :: KBN
qra_n :: Int
qra_y2 :: QuadRegAcc -> KBN
qra_x2y :: QuadRegAcc -> KBN
qra_xy :: QuadRegAcc -> KBN
qra_y :: QuadRegAcc -> KBN
qra_x4 :: QuadRegAcc -> KBN
qra_x3 :: QuadRegAcc -> KBN
qra_x2 :: QuadRegAcc -> KBN
qra_x :: QuadRegAcc -> KBN
qra_w :: QuadRegAcc -> KBN
qra_n :: QuadRegAcc -> Int
..} Double
x Double
y Double
w = QuadRegAcc
    { qra_n :: Int
qra_n    = Int
qra_n forall a. Num a => a -> a -> a
+ Int
1
    , qra_w :: KBN
qra_w    = KBN -> Double -> KBN
addKBN KBN
qra_w   Double
w
    , qra_x :: KBN
qra_x    = KBN -> Double -> KBN
addKBN KBN
qra_x   (Double
w forall a. Num a => a -> a -> a
* Double
x)
    , qra_x2 :: KBN
qra_x2   = KBN -> Double -> KBN
addKBN KBN
qra_x2  (Double
w forall a. Num a => a -> a -> a
* Double
x2)
    , qra_x3 :: KBN
qra_x3   = KBN -> Double -> KBN
addKBN KBN
qra_x3  (Double
w forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
* Double
x2)
    , qra_x4 :: KBN
qra_x4   = KBN -> Double -> KBN
addKBN KBN
qra_x4  (Double
w forall a. Num a => a -> a -> a
* Double
x2 forall a. Num a => a -> a -> a
* Double
x2)
    , qra_y :: KBN
qra_y    = KBN -> Double -> KBN
addKBN KBN
qra_y   (Double
w forall a. Num a => a -> a -> a
* Double
y)
    , qra_xy :: KBN
qra_xy   = KBN -> Double -> KBN
addKBN KBN
qra_xy  (Double
w forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
* Double
y)
    , qra_x2y :: KBN
qra_x2y  = KBN -> Double -> KBN
addKBN KBN
qra_x2y (Double
w forall a. Num a => a -> a -> a
* Double
x2 forall a. Num a => a -> a -> a
* Double
y)
    , qra_y2 :: KBN
qra_y2   = KBN -> Double -> KBN
addKBN KBN
qra_y2  (Double
w forall a. Num a => a -> a -> a
* Double
y forall a. Num a => a -> a -> a
* Double
y)
    }
  where
    x2 :: Double
x2 = Double
x forall a. Num a => a -> a -> a
* Double
x

-- | Convert 'QuadRegAcc' to 'LinRegAcc'.
--
-- Using this we can try quadratic and linear fits with a single data scan.
--
quadRegAccToLin :: QuadRegAcc -> LinRegAcc
quadRegAccToLin :: QuadRegAcc -> LinRegAcc
quadRegAccToLin QuadRegAcc {Int
KBN
qra_y2 :: KBN
qra_x2y :: KBN
qra_xy :: KBN
qra_y :: KBN
qra_x4 :: KBN
qra_x3 :: KBN
qra_x2 :: KBN
qra_x :: KBN
qra_w :: KBN
qra_n :: Int
qra_y2 :: QuadRegAcc -> KBN
qra_x2y :: QuadRegAcc -> KBN
qra_xy :: QuadRegAcc -> KBN
qra_y :: QuadRegAcc -> KBN
qra_x4 :: QuadRegAcc -> KBN
qra_x3 :: QuadRegAcc -> KBN
qra_x2 :: QuadRegAcc -> KBN
qra_x :: QuadRegAcc -> KBN
qra_w :: QuadRegAcc -> KBN
qra_n :: QuadRegAcc -> Int
..} = LinRegAcc
    { lra_n :: Int
lra_n  = Int
qra_n
    , lra_w :: KBN
lra_w  = KBN
qra_w
    , lra_x :: KBN
lra_x  = KBN
qra_x
    , lra_x2 :: KBN
lra_x2 = KBN
qra_x2
    , lra_y :: KBN
lra_y  = KBN
qra_y
    , lra_xy :: KBN
lra_xy = KBN
qra_xy
    , lra_y2 :: KBN
lra_y2 = KBN
qra_y2
    }

-------------------------------------------------------------------------------
-- utils
-------------------------------------------------------------------------------

sq :: Num a => a -> a
sq :: forall a. Num a => a -> a
sq a
x = a
x forall a. Num a => a -> a -> a
* a
x
{-# INLINE sq #-}

iterate1 :: (b -> b) -> b -> NE.NonEmpty b
iterate1 :: forall b. (b -> b) -> b -> NonEmpty b
iterate1 b -> b
g b
x = forall a. a -> NonEmpty a -> NonEmpty a
NE.cons b
x (forall b. (b -> b) -> b -> NonEmpty b
iterate1 b -> b
g (b -> b
g b
x))

eigenSM22 :: Double -> Double -> Double -> V2
eigenSM22 :: Double -> Double -> Double -> V2
eigenSM22 Double
a Double
b Double
c = Double -> Double -> V2
V2 ((Double
ac forall a. Num a => a -> a -> a
+ Double
discr) forall a. Fractional a => a -> a -> a
/ Double
2) ((Double
ac forall a. Num a => a -> a -> a
- Double
discr) forall a. Fractional a => a -> a -> a
/ Double
2)
  where
    ac :: Double
ac = Double
a forall a. Num a => a -> a -> a
+ Double
c
    discr :: Double
discr = forall a. Floating a => a -> a
sqrt (forall a. Num a => a -> a
sq (Double
a forall a. Num a => a -> a -> a
- Double
c) forall a. Num a => a -> a -> a
+ Double
4 forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
sq Double
b)

-- | Levenberg-Marquard stop condition
lmStop :: Double -> Double -> Double -> Bool
lmStop :: Double -> Double -> Double -> Bool
lmStop Double
lambda Double
wssr Double
wssr' =
    Double
lambda forall a. Ord a => a -> a -> Bool
< Double
1e-20 Bool -> Bool -> Bool
|| Double
lambda forall a. Ord a => a -> a -> Bool
> Double
1e20 Bool -> Bool -> Bool
|| forall a. RealFloat a => a -> Bool
isNaN Double
wssr' Bool -> Bool -> Bool
|| Double
relDiff forall a. Ord a => a -> a -> Bool
< Double
1e-10
  where
    relDiff :: Double
relDiff = forall a. Num a => a -> a
abs (Double
wssr' forall a. Num a => a -> a -> a
- Double
wssr) forall a. Fractional a => a -> a -> a
/ Double
wssr