summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL/Fitting.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/Fitting.hs')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Fitting.hs180
1 files changed, 0 insertions, 180 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
deleted file mode 100644
index 93fb281..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs
+++ /dev/null
@@ -1,180 +0,0 @@
1{- |
2Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5
6Maintainer : Alberto Ruiz
7Stability : provisional
8
9Nonlinear Least-Squares Fitting
10
11<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
12
13The example program in the GSL manual (see examples/fitting.hs):
14
15@
16dat = [
17 ([0.0],([6.0133918608118675],0.1)),
18 ([1.0],([5.5153769909966535],0.1)),
19 ([2.0],([5.261094606015287],0.1)),
20 ...
21 ([39.0],([1.0619821710802808],0.1))]
22
23expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
24
25expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
26
27(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
28@
29
30>>> path
31(6><5)
32 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
33 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662
34 , 3.0, 9.5807893736187, 4.948995119561291, 0.11942927999921617, 1.0945766509238248
35 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608
36 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375
37 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ]
38>>> sol
39[(5.045357818922331,6.027976702418132e-2),
40(0.10404905846029407,3.157045047172834e-3),
41(1.0192487112786812,3.782067731353722e-2)]
42
43-}
44
45
46module Numeric.GSL.Fitting (
47 -- * Levenberg-Marquardt
48 nlFitting, FittingMethod(..),
49 -- * Utilities
50 fitModelScaled, fitModel
51) where
52
53import Numeric.LinearAlgebra
54import Numeric.GSL.Internal
55
56import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
57import Foreign.C.Types
58import System.IO.Unsafe(unsafePerformIO)
59
60-------------------------------------------------------------------------
61
62type TVV = TV (TV Res)
63type TVM = TV (TM Res)
64
65data 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.
66 | 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.
67 deriving (Enum,Eq,Show,Bounded)
68
69
70-- | Nonlinear multidimensional least-squares fitting.
71nlFitting :: FittingMethod
72 -> Double -- ^ absolute tolerance
73 -> Double -- ^ relative tolerance
74 -> Int -- ^ maximum number of iterations allowed
75 -> (Vector Double -> Vector Double) -- ^ function to be minimized
76 -> (Vector Double -> Matrix Double) -- ^ Jacobian
77 -> Vector Double -- ^ starting point
78 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
79
80nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
81
82nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
83 let p = dim xiv
84 n = dim (f xiv)
85 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
86 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
87 rawpath <- createMatrix RowMajor maxit (2+p)
88 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
89 let it = round (rawpath @@> (maxit-1,0))
90 path = takeRows it rawpath
91 [sol] = toRows $ dropRows (it-1) path
92 freeHaskellFunPtr fp
93 freeHaskellFunPtr jp
94 return (subVector 2 p sol, path)
95
96foreign import ccall safe "nlfit"
97 c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
98
99-------------------------------------------------------
100
101checkdim1 n _p v
102 | dim v == n = v
103 | otherwise = error $ "Error: "++ show n
104 ++ " components expected in the result of the function supplied to nlFitting"
105
106checkdim2 n p m
107 | rows m == n && cols m == p = m
108 | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
109 ++ " Jacobian expected in nlFitting"
110
111------------------------------------------------------------
112
113err (model,deriv) dat vsol = zip sol errs where
114 sol = toList vsol
115 c = max 1 (chi/sqrt (fromIntegral dof))
116 dof = length dat - (rows cov)
117 chi = norm2 (fromList $ cost (resMs model) dat sol)
118 js = fromLists $ jacobian (resDs deriv) dat sol
119 cov = inv $ trans js <.> js
120 errs = toList $ scalar c * sqrt (takeDiag cov)
121
122
123
124-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
125-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
126-- instances (x, (y,sigma)) to be fitted.
127
128fitModelScaled
129 :: Double -- ^ absolute tolerance
130 -> Double -- ^ relative tolerance
131 -> Int -- ^ maximum number of iterations allowed
132 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
133 -> [(x, ([Double], Double))] -- ^ instances
134 -> [Double] -- ^ starting point
135 -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
136fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
137 (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
138 (fromList . cost (resMs model) dt . toList)
139 (fromLists . jacobian (resDs deriv) dt . toList)
140 (fromList xin)
141
142
143
144-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
145-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
146-- instances (x,y) to be fitted.
147
148fitModel :: Double -- ^ absolute tolerance
149 -> Double -- ^ relative tolerance
150 -> Int -- ^ maximum number of iterations allowed
151 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
152 -> [(x, [Double])] -- ^ instances
153 -> [Double] -- ^ starting point
154 -> ([Double], Matrix Double) -- ^ solution and optimization path
155fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
156 (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
157 (fromList . cost (resM model) dt . toList)
158 (fromLists . jacobian (resD deriv) dt . toList)
159 (fromList xin)
160
161cost model ds vs = concatMap (model vs) ds
162
163jacobian modelDer ds vs = concatMap (modelDer vs) ds
164
165-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
166resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
167resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
168
169-- | Associated derivative for 'resMs'.
170resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
171resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
172
173-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
174-- to 'resMs' with all sigmas = 1.
175resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
176resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
177
178-- | Associated derivative for 'resM'.
179resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
180resD m v = \(x,_) -> m v x
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 11:18:50 +0000