{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Numeric.GSL.Fitting (
nlFitting, FittingMethod(..),
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
| LevenbergMarquardt
deriving (Enum,Eq,Show,Bounded)
nlFitting :: FittingMethod
-> Double
-> Double
-> Int
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> (Vector Double, Matrix Double)
nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
let p = size xiv
n = size (f xiv)
fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
rawpath <- createMatrix RowMajor maxit (2+p)
(xiv `applyRaw` (rawpath `applyRaw` id)) (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) #|"c_nlfit"
let it = round (rawpath `atIndex` (maxit-1,0))
path = takeRows it rawpath
[sol] = toRows $ dropRows (it-1) path
freeHaskellFunPtr fp
freeHaskellFunPtr jp
return (subVector 2 p sol, path)
foreign import ccall safe "nlfit"
c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
checkdim1 n _p v
| size v == n = v
| otherwise = error $ "Error: "++ show n
++ " components expected in the result of the function supplied to nlFitting"
checkdim2 n p m
| rows m == n && cols m == p = m
| otherwise = error $ "Error: "++ show n ++ "x" ++ show p
++ " Jacobian expected in nlFitting"
err (model,deriv) dat vsol = zip sol errs where
sol = toList vsol
c = max 1 (chi/sqrt (fromIntegral dof))
dof = length dat - (rows cov)
chi = norm_2 (fromList $ cost (resMs model) dat sol)
js = fromLists $ jacobian (resDs deriv) dat sol
cov = inv $ tr js <> js
errs = toList $ scalar c * sqrt (takeDiag cov)
fitModelScaled
:: Double
-> Double
-> Int
-> ([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> [(x, ([Double], Double))]
-> [Double]
-> ([(Double, Double)], Matrix Double)
fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
(sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
(fromList . cost (resMs model) dt . toList)
(fromLists . jacobian (resDs deriv) dt . toList)
(fromList xin)
fitModel :: Double
-> Double
-> Int
-> ([Double] -> x -> [Double], [Double] -> x -> [[Double]])
-> [(x, [Double])]
-> [Double]
-> ([Double], Matrix Double)
fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
(sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
(fromList . cost (resM model) dt . toList)
(fromLists . jacobian (resD deriv) dt . toList)
(fromList xin)
cost model ds vs = concatMap (model vs) ds
jacobian modelDer ds vs = concatMap (modelDer vs) ds
resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
resD m v = \(x,_) -> m v x