7 files changed, 282 insertions, 15 deletions

diff --git a/CHANGES b/CHANGES
index fc089e0..190c0f4 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,8 @@
+0.9.2.0
+=======
+
+Added GSL Nonlinear Least-Squares fitting using Levenberg-Marquardt.
+
 0.9.1.0
 =======
 
@@ -6,7 +11,7 @@
 - Added offset of Vector, allowing fast, noncopy subVector (slice).
   Vector is now identical to Roman Leshchinskiy's Data.Vector.Storable.Vector,
   so we can convert from/to them in O(1).
-  
+
 - Removed Data.Packed.Convert, see examples/vector.hs
 
 0.8.3.0
diff --git a/examples/fitting.hs b/examples/fitting.hs
new file mode 100644
index 0000000..8298c52
--- /dev/null
+++ b/examples/fitting.hs
@@ -0,0 +1,24 @@
+-- nonlinear least-squares fitting
+
+import Numeric.GSL.Fitting
+import Numeric.LinearAlgebra
+
+xs = map return [0 .. 39]
+sigma = 0.1
+ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)
+            + scalar sigma * (randomVector 0 Gaussian 40)
+
+dat :: [([Double],[Double],Double)]
+
+dat = zipWith3 (,,) xs ys (repeat sigma)
+
+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) = fitModel 1E-4 1E-4 20 (resM expModel, resD expModelDer) dat [1,0,0]
+
+main = do
+    print dat
+    print path
+    print sol
diff --git a/hmatrix.cabal b/hmatrix.cabal
index ac38460..5416c09 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.9.1.0
+Version:            0.9.2.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -79,6 +79,7 @@ library
                         Numeric.GSL.Polynomials,
                         Numeric.GSL.Minimization,
                         Numeric.GSL.Root,
+                        Numeric.GSL.Fitting,
                         Numeric.GSL.ODE,
                         Numeric.GSL.Vector,
                         Numeric.GSL,
diff --git a/lib/Numeric/GSL.hs b/lib/Numeric/GSL.hs
index 5442586..ae4c81e 100644
--- a/lib/Numeric/GSL.hs
+++ b/lib/Numeric/GSL.hs
@@ -20,6 +20,7 @@ module Numeric.GSL (
 , module Numeric.GSL.Minimization
 , module Numeric.GSL.Root
 , module Numeric.GSL.ODE
+, module Numeric.GSL.Fitting
 , module Data.Complex
 , setErrorHandlerOff
 ) where
@@ -31,6 +32,7 @@ import Numeric.GSL.Polynomials
 import Numeric.GSL.Minimization
 import Numeric.GSL.Root
 import Numeric.GSL.ODE
+import Numeric.GSL.Fitting
 import Data.Complex
 
 
diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs
new file mode 100644
index 0000000..6e91b2d
--- /dev/null
+++ b/lib/Numeric/GSL/Fitting.hs
@@ -0,0 +1,147 @@
+{- |
+Module      :  Numeric.GSL.Fitting
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+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) = fitModel 1E-4 1E-4 20 (resM expModel, resD 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
+    fitModel, resM, resD
+) where
+
+import Data.Packed.Internal
+import Numeric.LinearAlgebra
+import Foreign
+import Foreign.C.Types(CInt)
+import Numeric.GSL.Internal
+
+-------------------------------------------------------------------------
+
+data FittingMethod = LevenbergMarquardt -- ^ 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.
+        deriving (Enum,Eq,Show,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 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   = dim xiv
+        n   = dim (f xiv)
+    fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
+    jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
+    rawpath <- createMatrix RowMajor maxit (2+p)
+    app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
+    let it = round (rawpath @@> (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 "nlfit"
+    c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
+
+-------------------------------------------------------
+
+checkdim1 n _p v
+    | dim 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 dat vsol = zip sol errs where
+    sol = toList vsol
+    c = max 1 (chi/sqrt (fromIntegral dof))
+    dof = length dat - (rows cov)
+    chi = pnorm PNorm2 (fromList $ cost (fst model) dat sol)
+    js = fromLists $ jacobian (snd model) dat sol
+    cov = inv $ trans js <> js
+    errs = toList $ scalar c * sqrt (takeDiag cov)
+
+
+
+-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
+-- Jacobian are automatically built from a model f vs x = 0 and its derivatives, and a list of
+-- instances x 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]    -- ^ instances
+         -> [Double] -- ^ starting point
+         -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
+fitModel epsabs epsrel maxit model dt xin = (err model dt sol, path) where
+    (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
+                (fromList . cost (fst model) dt . toList)
+                (fromLists . jacobian (snd model) dt . toList)
+                (fromList xin)
+
+cost model ds vs = concatMap (model vs) ds
+
+jacobian modelDer ds vs = concatMap (modelDer vs) ds
+
+-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
+resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double], Double) -> [Double]
+resM m v = \(x,ys,s) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
+
+-- | Associated derivative for 'resM'.
+resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double], Double) -> [[Double]]
+resD m v = \(x,_,s) -> map (map (/s)) (m v x)
+
+
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 0fcde1c..f258fb1 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -23,6 +23,7 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_odeiv.h>
+#include <gsl/gsl_multifit_nlin.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -668,9 +669,10 @@ int root(int method, void f(int, double*, int, double*),
 
 // working with the jacobian
 
-typedef struct {int (*f)(int, double*, int, double *); int (*jf)(int, double*, int, int, double*);} Tfjf;
+typedef struct {int (*f)(int, double*, int, double *);
+                int (*jf)(int, double*, int, int, double*);} Tfjf;
 
-int f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
+int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
     Tfjf * fjf = ((Tfjf*) pars);
     double* p = (double*)calloc(x->size,sizeof(double));
     double* q = (double*)calloc(y->size,sizeof(double));
@@ -687,20 +689,20 @@ int f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
     return 0;
 }
 
-int jf_aux_root(const gsl_vector * x, void * pars, gsl_matrix * jac) {
+int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
     Tfjf * fjf = ((Tfjf*) pars);
     double* p = (double*)calloc(x->size,sizeof(double));
-    double* q = (double*)calloc((x->size)*(x->size),sizeof(double));
+    double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
     int i,j,k;
     for(k=0;k<x->size;k++) {
         p[k] = gsl_vector_get(x,k);
     }
 
-    (fjf->jf)(x->size,p,x->size,x->size,q);
+    (fjf->jf)(x->size,p,jac->size1,jac->size2,q);
 
     k=0;
-    for(i=0;i<x->size;i++) {
-        for(j=0;j<x->size;j++){
+    for(i=0;i<jac->size1;i++) {
+        for(j=0;j<jac->size2;j++){
             gsl_matrix_set(jac,i,j,q[k++]);
         }
     }
@@ -709,9 +711,9 @@ int jf_aux_root(const gsl_vector * x, void * pars, gsl_matrix * jac) {
     return 0;
 }
 
-int fjf_aux_root(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
-    f_aux_root(x,pars,f);
-    jf_aux_root(x,pars,g);
+int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
+    f_aux(x,pars,f);
+    jf_aux(x,pars,g);
     return 0;
 }
 
@@ -723,9 +725,9 @@ int rootj(int method, int f(int, double*, int, double*),
     DEBUGMSG("root_fjf");
     gsl_multiroot_function_fdf my_func;
     // extract function from pars
-    my_func.f = f_aux_root;
-    my_func.df = jf_aux_root;
-    my_func.fdf = fjf_aux_root;
+    my_func.f = f_aux;
+    my_func.df = jf_aux;
+    my_func.fdf = fjf_aux;
     my_func.n = xin;
     Tfjf stfjf;
     stfjf.f = f;
@@ -781,8 +783,77 @@ int rootj(int method, int f(int, double*, int, double*),
     OK
 }
 
+//-------------- non linear least squares fitting -------------------
+
+int nlfit(int method, int f(int, double*, int, double*),
+                      int jac(int, double*, int, int, double*),
+         double epsabs, double epsrel, int maxit, int p,
+         KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
+    DEBUGMSG("nlfit");
+    const gsl_multifit_fdfsolver_type *T;
+    gsl_multifit_fdfsolver *s;
+    gsl_multifit_function_fdf my_f;
+    // extract function from pars
+    my_f.f = f_aux;
+    my_f.df = jf_aux;
+    my_f.fdf = fjf_aux;
+    my_f.n = p;
+    my_f.p = xin;  // !!!!
+    Tfjf stfjf;
+    stfjf.f = f;
+    stfjf.jf = jac;
+    my_f.params = &stfjf;
+    size_t iter = 0;
+    int status;
+
+    KDVVIEW(xi);
+    //DMVIEW(cov);
+
+    switch(method) {
+        case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
+        default: ERROR(BAD_CODE);
+    }
+
+    s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
+    gsl_multifit_fdfsolver_set (s, &my_f, V(xi));
+
+    do {   status = gsl_multifit_fdfsolver_iterate (s);
+
+           solp[iter*solc+0] = iter+1;
+           solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);
+
+           int k;
+           for(k=0;k<xin;k++) {
+               solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
+           }
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter <= maxit);
+
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+
+    //gsl_multifit_covar (s->J, 0.0, M(cov));
+
+    gsl_multifit_fdfsolver_free (s);
+    OK
+}
+
+
 //////////////////////////////////////////////////////
 
+
 #define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
 
 int random_vector(int seed, int code, RVEC(r)) {
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 3f7c847..f8f8bd5 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -143,6 +143,22 @@ odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.36434890871
 
 ---------------------------------------------------------------------
 
+fittingTest = utest "levmar" ok
+    where
+    xs = map return [0 .. 39]
+    sigma = 0.1
+    ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)
+                    + scalar sigma * (randomVector 0 Gaussian 40)
+    dat = zipWith3 (,,) xs ys (repeat sigma)
+
+    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 = fst $ fitModel 1E-4 1E-4 20 (resM expModel, resD expModelDer) dat [1,0,0]
+    ok = and (zipWith f sol [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d
+
+---------------------------------------------------------------------
+
 randomTestGaussian = c :~1~: snd (meanCov dat) where
     a = (3><3) [1,2,3,
                 2,4,0,
@@ -292,6 +308,7 @@ runTests n = do
                        && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
         , odeTest
+        , fittingTest
         ]
     return ()