{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

{- |
Module      :  Numeric.GSL.Fitting
Copyright   :  (c) Alberto Ruiz 2010
License     :  GPL

Maintainer  :  Alberto Ruiz
Stability   :  provisional

Nonlinear Least-Squares Fitting

<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>

The example program in the GSL manual (see examples/fitting.hs):

@
dat = [
 ([0.0],([6.0133918608118675],0.1)),
 ([1.0],([5.5153769909966535],0.1)),
 ([2.0],([5.261094606015287],0.1)),
 ...
 ([39.0],([1.0619821710802808],0.1))]

expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]

expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]

(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
@

>>> path
(6><5)
 [ 1.0,  76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
 , 2.0, 37.683816318260355,  2.858760367632973,  8.092094813253975e-2, 1.4479636296208662
 , 3.0,    9.5807893736187,  4.948995119561291,   0.11942927999921617, 1.0945766509238248
 , 4.0,  5.630494933603935,  5.021755718065913,   0.10287787128056883, 1.0338835440862608
 , 5.0,  5.443976278682909,  5.045204331329302,   0.10405523433131504,  1.019416067207375
 , 6.0, 5.4439736648994685,  5.045357818922331,   0.10404905846029407, 1.0192487112786812 ]
>>> sol
[(5.045357818922331,6.027976702418132e-2),
(0.10404905846029407,3.157045047172834e-3),
(1.0192487112786812,3.782067731353722e-2)]

-}


module Numeric.GSL.Fitting (
    -- * Levenberg-Marquardt
    nlFitting, FittingMethod(..),
    -- * Utilities
    fitModelScaled, fitModel
) where

import Numeric.LinearAlgebra.HMatrix
import Numeric.GSL.Internal

import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

-------------------------------------------------------------------------

type TVV = TV (TV Res)
type TVM = TV (TM Res)

data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom.
                   | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately.
        deriving (Int -> FittingMethod
FittingMethod -> Int
FittingMethod -> [FittingMethod]
FittingMethod -> FittingMethod
FittingMethod -> FittingMethod -> [FittingMethod]
FittingMethod -> FittingMethod -> FittingMethod -> [FittingMethod]
(FittingMethod -> FittingMethod)
-> (FittingMethod -> FittingMethod)
-> (Int -> FittingMethod)
-> (FittingMethod -> Int)
-> (FittingMethod -> [FittingMethod])
-> (FittingMethod -> FittingMethod -> [FittingMethod])
-> (FittingMethod -> FittingMethod -> [FittingMethod])
-> (FittingMethod
    -> FittingMethod -> FittingMethod -> [FittingMethod])
-> Enum FittingMethod
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: FittingMethod -> FittingMethod -> FittingMethod -> [FittingMethod]
$cenumFromThenTo :: FittingMethod -> FittingMethod -> FittingMethod -> [FittingMethod]
enumFromTo :: FittingMethod -> FittingMethod -> [FittingMethod]
$cenumFromTo :: FittingMethod -> FittingMethod -> [FittingMethod]
enumFromThen :: FittingMethod -> FittingMethod -> [FittingMethod]
$cenumFromThen :: FittingMethod -> FittingMethod -> [FittingMethod]
enumFrom :: FittingMethod -> [FittingMethod]
$cenumFrom :: FittingMethod -> [FittingMethod]
fromEnum :: FittingMethod -> Int
$cfromEnum :: FittingMethod -> Int
toEnum :: Int -> FittingMethod
$ctoEnum :: Int -> FittingMethod
pred :: FittingMethod -> FittingMethod
$cpred :: FittingMethod -> FittingMethod
succ :: FittingMethod -> FittingMethod
$csucc :: FittingMethod -> FittingMethod
Enum,FittingMethod -> FittingMethod -> Bool
(FittingMethod -> FittingMethod -> Bool)
-> (FittingMethod -> FittingMethod -> Bool) -> Eq FittingMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: FittingMethod -> FittingMethod -> Bool
$c/= :: FittingMethod -> FittingMethod -> Bool
== :: FittingMethod -> FittingMethod -> Bool
$c== :: FittingMethod -> FittingMethod -> Bool
Eq,Int -> FittingMethod -> ShowS
[FittingMethod] -> ShowS
FittingMethod -> String
(Int -> FittingMethod -> ShowS)
-> (FittingMethod -> String)
-> ([FittingMethod] -> ShowS)
-> Show FittingMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [FittingMethod] -> ShowS
$cshowList :: [FittingMethod] -> ShowS
show :: FittingMethod -> String
$cshow :: FittingMethod -> String
showsPrec :: Int -> FittingMethod -> ShowS
$cshowsPrec :: Int -> FittingMethod -> ShowS
Show,FittingMethod
FittingMethod -> FittingMethod -> Bounded FittingMethod
forall a. a -> a -> Bounded a
maxBound :: FittingMethod
$cmaxBound :: FittingMethod
minBound :: FittingMethod
$cminBound :: FittingMethod
Bounded)


-- | Nonlinear multidimensional least-squares fitting.
nlFitting :: FittingMethod
      -> Double                     -- ^ absolute tolerance
      -> Double                     -- ^ relative tolerance
      -> Int                        -- ^ maximum number of iterations allowed
      -> (Vector Double -> Vector Double)     -- ^ function to be minimized
      -> (Vector Double -> Matrix Double)   -- ^ Jacobian
      -> Vector Double                   -- ^ starting point
      -> (Vector Double, Matrix Double)  -- ^ solution vector and optimization path

nlFitting :: FittingMethod
-> Double
-> Double
-> Int
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> (Vector Double, Matrix Double)
nlFitting FittingMethod
method Double
epsabs Double
epsrel Int
maxit Vector Double -> Vector Double
fun Vector Double -> Matrix Double
jac Vector Double
xinit = CInt
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Double
-> Double
-> Int
-> (Vector Double, Matrix Double)
nlFitGen (Int -> CInt
fi (FittingMethod -> Int
forall a. Enum a => a -> Int
fromEnum FittingMethod
method)) Vector Double -> Vector Double
fun Vector Double -> Matrix Double
jac Vector Double
xinit Double
epsabs Double
epsrel Int
maxit

nlFitGen :: CInt
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Double
-> Double
-> Int
-> (Vector Double, Matrix Double)
nlFitGen CInt
m Vector Double -> Vector Double
f Vector Double -> Matrix Double
jac Vector Double
xiv Double
epsabs Double
epsrel Int
maxit = IO (Vector Double, Matrix Double) -> (Vector Double, Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO (Vector Double, Matrix Double)
 -> (Vector Double, Matrix Double))
-> IO (Vector Double, Matrix Double)
-> (Vector Double, Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    let p :: IndexOf Vector
p   = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Vector Double
xiv
        n :: IndexOf Vector
n   = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size (Vector Double -> Vector Double
f Vector Double
xiv)
    FunPtr TVV
fp <- TVV -> IO (FunPtr TVV)
mkVecVecfun ((Vector Double -> Vector Double) -> TVV
aux_vTov (IndexOf Vector -> Int -> Vector Double -> Vector Double
forall {c :: * -> *} {t} {p}.
(Eq (IndexOf c), Container c t, Show (IndexOf c)) =>
IndexOf c -> p -> c t -> c t
checkdim1 IndexOf Vector
n Int
IndexOf Vector
p (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Vector Double
f))
    FunPtr TVM
jp <- TVM -> IO (FunPtr TVM)
mkVecMatfun ((Vector Double -> Matrix Double) -> TVM
aux_vTom (Int -> Int -> Matrix Double -> Matrix Double
forall {t}. Int -> Int -> Matrix t -> Matrix t
checkdim2 Int
IndexOf Vector
n Int
IndexOf Vector
p (Matrix Double -> Matrix Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Matrix Double
jac))
    Matrix Double
rawpath <- MatrixOrder -> Int -> Int -> IO (Matrix Double)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
RowMajor Int
maxit (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
IndexOf Vector
p)
    (Vector Double
xiv Vector Double
-> ((CInt -> CInt -> Ptr Double -> IO CInt) -> IO CInt)
-> TransRaw (Vector Double) (CInt -> CInt -> Ptr Double -> IO CInt)
-> IO CInt
forall c b r.
TransArray c =>
c -> (b -> IO r) -> TransRaw c b -> IO r
`applyRaw` (Matrix Double
rawpath Matrix Double
-> (IO CInt -> IO CInt)
-> TransRaw (Matrix Double) (IO CInt)
-> IO CInt
forall c b r.
TransArray c =>
c -> (b -> IO r) -> TransRaw c b -> IO r
`applyRaw` IO CInt -> IO CInt
forall a. a -> a
id)) (CInt
-> FunPtr TVV
-> FunPtr TVM
-> Double
-> Double
-> CInt
-> CInt
-> TVM
c_nlfit CInt
m FunPtr TVV
fp FunPtr TVM
jp Double
epsabs Double
epsrel (Int -> CInt
fi Int
maxit) (Int -> CInt
fi Int
IndexOf Vector
n)) IO CInt -> String -> IO ()
#|String
"c_nlfit"
    let it :: Int
it = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Matrix Double
rawpath Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
maxitInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
0))
        path :: Matrix Double
path = Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
takeRows Int
it Matrix Double
rawpath
        [Vector Double
sol] = Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix Double -> [Vector Double])
-> Matrix Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
dropRows (Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Matrix Double
path
    FunPtr TVV -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr TVV
fp
    FunPtr TVM -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr TVM
jp
    (Vector Double, Matrix Double) -> IO (Vector Double, Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> Int -> Vector Double -> Vector Double
forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector Int
2 Int
IndexOf Vector
p Vector Double
sol, Matrix Double
path)

foreign import ccall safe "nlfit"
    c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM

-------------------------------------------------------

checkdim1 :: IndexOf c -> p -> c t -> c t
checkdim1 IndexOf c
n p
_p c t
v
    | c t -> IndexOf c
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size c t
v IndexOf c -> IndexOf c -> Bool
forall a. Eq a => a -> a -> Bool
== IndexOf c
n = c t
v
    | Bool
otherwise = String -> c t
forall a. HasCallStack => String -> a
error (String -> c t) -> String -> c t
forall a b. (a -> b) -> a -> b
$ String
"Error: "String -> ShowS
forall a. [a] -> [a] -> [a]
++ IndexOf c -> String
forall a. Show a => a -> String
show IndexOf c
n
                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" components expected in the result of the function supplied to nlFitting"

checkdim2 :: Int -> Int -> Matrix t -> Matrix t
checkdim2 Int
n Int
p Matrix t
m
    | Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Bool -> Bool -> Bool
&& Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
p = Matrix t
m
    | Bool
otherwise = String -> Matrix t
forall a. HasCallStack => String -> a
error (String -> Matrix t) -> String -> Matrix t
forall a b. (a -> b) -> a -> b
$ String
"Error: "String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"x" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
p
                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" Jacobian expected in nlFitting"

------------------------------------------------------------

err :: ([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> t (x, ([Double], Double)) -> Vector Double -> [(Double, Double)]
err ([Double] -> x -> [Double]
model,[Double] -> x -> [[Double]]
deriv) t (x, ([Double], Double))
dat Vector Double
vsol = [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
sol [Double]
errs where
    sol :: [Double]
sol = Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList Vector Double
vsol
    c :: Double
c = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double
chiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dof))
    dof :: Int
dof = t (x, ([Double], Double)) -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length t (x, ([Double], Double))
dat Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
cov)
    chi :: Double
chi = Vector Double -> Double
forall a. Normed a => a -> Double
norm_2 ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double) -> [Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$ ([Double] -> (x, ([Double], Double)) -> [Double])
-> t (x, ([Double], Double)) -> [Double] -> [Double]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
cost (([Double] -> x -> [Double])
-> [Double] -> (x, ([Double], Double)) -> [Double]
forall x.
([Double] -> x -> [Double])
-> [Double] -> (x, ([Double], Double)) -> [Double]
resMs [Double] -> x -> [Double]
model) t (x, ([Double], Double))
dat [Double]
sol)
    js :: Matrix Double
js = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
fromLists ([[Double]] -> Matrix Double) -> [[Double]] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ ([Double] -> (x, ([Double], Double)) -> [[Double]])
-> t (x, ([Double], Double)) -> [Double] -> [[Double]]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
jacobian (([Double] -> x -> [[Double]])
-> [Double] -> (x, ([Double], Double)) -> [[Double]]
forall x.
([Double] -> x -> [[Double]])
-> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs [Double] -> x -> [[Double]]
deriv) t (x, ([Double], Double))
dat [Double]
sol
    cov :: Matrix Double
cov = Matrix Double -> Matrix Double
forall t. Field t => Matrix t -> Matrix t
inv (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
tr Matrix Double
js Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
<> Matrix Double
js
    errs :: [Double]
errs = Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList (Vector Double -> [Double]) -> Vector Double -> [Double]
forall a b. (a -> b) -> a -> b
$ Double -> Vector Double
forall (c :: * -> *) e. Container c e => e -> c e
scalar Double
c Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
* Vector Double -> Vector Double
forall a. Floating a => a -> a
sqrt (Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
takeDiag Matrix Double
cov)



-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x, (y,sigma)) to be fitted.

fitModelScaled
         :: Double -- ^ absolute tolerance
         -> Double -- ^ relative tolerance
         -> Int    -- ^ maximum number of iterations allowed
         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
         -> [(x, ([Double], Double))] -- ^ instances
         -> [Double] -- ^ starting point
         -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
fitModelScaled :: forall x.
Double
-> Double
-> Int
-> ([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> [(x, ([Double], Double))]
-> [Double]
-> ([(Double, Double)], Matrix Double)
fitModelScaled Double
epsabs Double
epsrel Int
maxit ([Double] -> x -> [Double]
model,[Double] -> x -> [[Double]]
deriv) [(x, ([Double], Double))]
dt [Double]
xin = (([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> [(x, ([Double], Double))] -> Vector Double -> [(Double, Double)]
forall {t :: * -> *} {x}.
Foldable t =>
([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> t (x, ([Double], Double)) -> Vector Double -> [(Double, Double)]
err ([Double] -> x -> [Double]
model,[Double] -> x -> [[Double]]
deriv) [(x, ([Double], Double))]
dt Vector Double
sol, Matrix Double
path) where
    (Vector Double
sol,Matrix Double
path) = FittingMethod
-> Double
-> Double
-> Int
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> (Vector Double, Matrix Double)
nlFitting FittingMethod
LevenbergMarquardtScaled Double
epsabs Double
epsrel Int
maxit
                ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> (x, ([Double], Double)) -> [Double])
-> [(x, ([Double], Double))] -> [Double] -> [Double]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
cost (([Double] -> x -> [Double])
-> [Double] -> (x, ([Double], Double)) -> [Double]
forall x.
([Double] -> x -> [Double])
-> [Double] -> (x, ([Double], Double)) -> [Double]
resMs [Double] -> x -> [Double]
model) [(x, ([Double], Double))]
dt ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList)
                ([[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
fromLists ([[Double]] -> Matrix Double)
-> (Vector Double -> [[Double]]) -> Vector Double -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> (x, ([Double], Double)) -> [[Double]])
-> [(x, ([Double], Double))] -> [Double] -> [[Double]]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
jacobian (([Double] -> x -> [[Double]])
-> [Double] -> (x, ([Double], Double)) -> [[Double]]
forall x.
([Double] -> x -> [[Double]])
-> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs [Double] -> x -> [[Double]]
deriv) [(x, ([Double], Double))]
dt ([Double] -> [[Double]])
-> (Vector Double -> [Double]) -> Vector Double -> [[Double]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList)
                ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xin)



-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x,y) to be fitted.

fitModel :: Double -- ^ absolute tolerance
         -> Double -- ^ relative tolerance
         -> Int    -- ^ maximum number of iterations allowed
         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
         -> [(x, [Double])]  -- ^ instances
         -> [Double] -- ^ starting point
         -> ([Double], Matrix Double) -- ^ solution and optimization path
fitModel :: forall x.
Double
-> Double
-> Int
-> ([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> [(x, [Double])]
-> [Double]
-> ([Double], Matrix Double)
fitModel Double
epsabs Double
epsrel Int
maxit ([Double] -> x -> [Double]
model,[Double] -> x -> [[Double]]
deriv) [(x, [Double])]
dt [Double]
xin = (Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList Vector Double
sol, Matrix Double
path) where
    (Vector Double
sol,Matrix Double
path) = FittingMethod
-> Double
-> Double
-> Int
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> (Vector Double, Matrix Double)
nlFitting FittingMethod
LevenbergMarquardt Double
epsabs Double
epsrel Int
maxit
                ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> (x, [Double]) -> [Double])
-> [(x, [Double])] -> [Double] -> [Double]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
cost (([Double] -> x -> [Double])
-> [Double] -> (x, [Double]) -> [Double]
forall x.
([Double] -> x -> [Double])
-> [Double] -> (x, [Double]) -> [Double]
resM [Double] -> x -> [Double]
model) [(x, [Double])]
dt ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList)
                ([[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
fromLists ([[Double]] -> Matrix Double)
-> (Vector Double -> [[Double]]) -> Vector Double -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> (x, [Double]) -> [[Double]])
-> [(x, [Double])] -> [Double] -> [[Double]]
forall {t :: * -> *} {t} {a} {b}.
Foldable t =>
(t -> a -> [b]) -> t a -> t -> [b]
jacobian (([Double] -> x -> [[Double]])
-> [Double] -> (x, [Double]) -> [[Double]]
forall x.
([Double] -> x -> [[Double]])
-> [Double] -> (x, [Double]) -> [[Double]]
resD [Double] -> x -> [[Double]]
deriv) [(x, [Double])]
dt ([Double] -> [[Double]])
-> (Vector Double -> [Double]) -> Vector Double -> [[Double]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList)
                ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xin)

cost :: (t -> a -> [b]) -> t a -> t -> [b]
cost t -> a -> [b]
model t a
ds t
vs = (a -> [b]) -> t a -> [b]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (t -> a -> [b]
model t
vs) t a
ds

jacobian :: (t -> a -> [b]) -> t a -> t -> [b]
jacobian t -> a -> [b]
modelDer t a
ds t
vs = (a -> [b]) -> t a -> [b]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (t -> a -> [b]
modelDer t
vs) t a
ds

-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
resMs :: forall x.
([Double] -> x -> [Double])
-> [Double] -> (x, ([Double], Double)) -> [Double]
resMs [Double] -> x -> [Double]
m [Double]
v = \(x
x,([Double]
ys,Double
s)) -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Double -> Double -> Double -> Double
forall {a}. Fractional a => a -> a -> a -> a
g Double
s) ([Double] -> x -> [Double]
m [Double]
v x
x) [Double]
ys where g :: a -> a -> a -> a
g a
s a
a a
b = (a
aa -> a -> a
forall a. Num a => a -> a -> a
-a
b)a -> a -> a
forall a. Fractional a => a -> a -> a
/a
s

-- | Associated derivative for 'resMs'.
resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs :: forall x.
([Double] -> x -> [[Double]])
-> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs [Double] -> x -> [[Double]]
m [Double]
v = \(x
x,([Double]
_,Double
s)) -> ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
s)) ([Double] -> x -> [[Double]]
m [Double]
v x
x)

-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
-- to 'resMs' with all sigmas = 1.
resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
resM :: forall x.
([Double] -> x -> [Double])
-> [Double] -> (x, [Double]) -> [Double]
resM [Double] -> x -> [Double]
m [Double]
v = \(x
x,[Double]
ys) -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
g ([Double] -> x -> [Double]
m [Double]
v x
x) [Double]
ys where g :: a -> a -> a
g a
a a
b = a
aa -> a -> a
forall a. Num a => a -> a -> a
-a
b

-- | Associated derivative for 'resM'.
resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
resD :: forall x.
([Double] -> x -> [[Double]])
-> [Double] -> (x, [Double]) -> [[Double]]
resD [Double] -> x -> [[Double]]
m [Double]
v = \(x
x,[Double]
_) -> [Double] -> x -> [[Double]]
m [Double]
v x
x