added some root finding algorithms

7 files changed, 256 insertions, 12 deletions

diff --git a/examples/root.hs b/examples/root.hs
new file mode 100644
index 0000000..9a674fd
--- /dev/null
+++ b/examples/root.hs
@@ -0,0 +1,30 @@
+-- root finding examples
+import Numeric.GSL
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import Text.Printf(printf)
+
+rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+
+disp = putStrLn . format "  " (printf "%.3f")
+
+-- Numerical estimation of the gradient
+gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]
+
+partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
+    g x = f (concat [a,x:b])
+    (a,_:b) = splitAt n v
+
+test method = do
+    print method
+    let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5]
+    print s -- solution
+    disp p -- evolution of the algorithm
+--    let [x,y] = tail (toColumns p)
+--    mplot [x,y]  -- path from the starting point to the solution
+
+main = do
+    test Hybrids
+    test Hybrid
+    test DNewton
+    test Broyden
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 87436ae..70de11e 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -24,6 +24,7 @@ extra-source-files: examples/tests.hs
                     examples/deriv.hs
                     examples/integrate.hs
                     examples/minimize.hs
+                    examples/root.hs
                     examples/pca1.hs
                     examples/pca2.hs
                     examples/pinv.hs
@@ -80,6 +81,7 @@ library
                         Numeric.GSL.Fourier,
                         Numeric.GSL.Polynomials,
                         Numeric.GSL.Minimization,
+                        Numeric.GSL.Root,
                         Numeric.GSL.Vector,
                         Numeric.GSL.Special,
                         Numeric.GSL.Special.Gamma,
diff --git a/lib/Numeric/GSL.hs b/lib/Numeric/GSL.hs
index 32962ef..2e90fff 100644
--- a/lib/Numeric/GSL.hs
+++ b/lib/Numeric/GSL.hs
@@ -18,6 +18,7 @@ module Numeric.GSL (
 , module Numeric.GSL.Fourier
 , module Numeric.GSL.Polynomials
 , module Numeric.GSL.Minimization
+, module Numeric.GSL.Root
 , module Numeric.GSL.Special
 , module Complex
 , setErrorHandlerOff
@@ -29,6 +30,7 @@ import Numeric.GSL.Special
 import Numeric.GSL.Fourier
 import Numeric.GSL.Polynomials
 import Numeric.GSL.Minimization
+import Numeric.GSL.Root
 import Complex
 import Numeric.GSL.Special
 
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
new file mode 100644
index 0000000..ad1b72c
--- /dev/null
+++ b/lib/Numeric/GSL/Root.hs
@@ -0,0 +1,117 @@
+{- |
+Module      :  Numeric.GSL.Root
+Copyright   :  (c) Alberto Ruiz 2009
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Multidimensional root finding.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
+
+The example in the GSL manual:
+
+@import Numeric.GSL
+import Numeric.LinearAlgebra(format)
+import Text.Printf(printf)
+
+rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+
+disp = putStrLn . format \"  \" (printf \"%.3f\")
+
+main = do
+    let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
+    print sol
+    disp path
+
+\> main
+[1.0,1.0]
+ 0.000  -10.000  -5.000  11.000  -1050.000
+ 1.000   -3.976  24.827   4.976     90.203
+ 2.000   -3.976  24.827   4.976     90.203
+ 3.000   -3.976  24.827   4.976     90.203
+ 4.000   -1.274  -5.680   2.274    -73.018
+ 5.000   -1.274  -5.680   2.274    -73.018
+ 6.000    0.249   0.298   0.751      2.359
+ 7.000    0.249   0.298   0.751      2.359
+ 8.000    1.000   0.878  -0.000     -1.218
+ 9.000    1.000   0.989  -0.000     -0.108
+10.000    1.000   1.000   0.000      0.000
+@
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Root (
+    root, RootMethod(..)
+) where
+
+import Data.Packed.Internal
+import Data.Packed.Matrix
+import Foreign
+import Foreign.C.Types(CInt)
+
+-------------------------------------------------------------------------
+
+data RootMethod = Hybrids
+                | Hybrid
+                | DNewton
+                | Broyden
+                deriving (Enum,Eq,Show)
+
+-- | Nonlinear multidimensional root finding using algorithms that do not require 
+-- any derivative information to be supplied by the user.
+-- Any derivatives needed are approximated by finite differences.
+root :: RootMethod
+     -> Double                     -- ^ maximum residual
+     -> Int                        -- ^ maximum number of iterations allowed
+     -> ([Double] -> [Double])     -- ^ function to minimize
+     -> [Double]                   -- ^ starting point
+     -> ([Double], Matrix Double)  -- ^ solution vector and optimization path
+
+root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
+
+rootGen m f xi epsabs maxit = unsafePerformIO $ do
+    let xiv = fromList xi
+        n   = dim xiv
+    fp <- mkVecVecfun (aux_vTov (fromList.f.toList))
+    rawpath <- withVector xiv $ \xiv' ->
+                   createMIO maxit (2*n+1)
+                         (c_root m fp epsabs (fi maxit) // xiv')
+                         "root"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (take n $ drop 1 sol, path)
+
+
+foreign import ccall "root"
+    c_root:: CInt -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ()) -> Double -> CInt -> TVM
+
+---------------------------------------------------------------------
+
+foreign import ccall "wrapper"
+    mkVecVecfun :: (CInt -> Ptr Double -> Ptr Double -> IO ())
+                -> IO (FunPtr (CInt -> Ptr Double -> Ptr Double->IO()))
+
+aux_vTov :: (Vector Double -> Vector Double) -> (CInt -> Ptr Double -> Ptr Double -> IO())
+aux_vTov f n p r = g where
+    V {fptr = pr} = f x
+    x = createV (fromIntegral n) copy "aux_vTov"
+    copy n' q = do
+        copyArray q p (fromIntegral n')
+        return 0
+    g = withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral n)
+
+createV n fun msg = unsafePerformIO $ do
+    r <- createVector n
+    app1 fun vec r msg
+    return r
+
+createMIO r c fun msg = do
+    res <- createMatrix RowMajor r c
+    app1 fun mat res msg
+    return res
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 3802574..80c23fc 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -7,6 +7,7 @@
 #include <gsl/gsl_deriv.h>
 #include <gsl/gsl_poly.h>
 #include <gsl/gsl_multimin.h>
+#include <gsl/gsl_multiroots.h>
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
 #include <string.h>
@@ -288,6 +289,22 @@ int fft(int code, KCVEC(X), CVEC(R)) {
 }
 
 
+int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
+{
+    gsl_function F;
+    F.function = f;
+    F.params = 0;
+
+    if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
+
+    if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
+
+    if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
+
+    return 0;
+}
+
+
 int integrate_qng(double f(double, void*), double a, double b, double prec,
                    double *result, double*error) {
     DEBUGMSG("integrate_qng");
@@ -440,7 +457,7 @@ void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g)
     df_aux_min(x,pars,g);
 }
 
-// conjugate gradient
+
 int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*, double*), 
                       double initstep, double minimpar, double tolgrad, int maxit, 
                       KRVEC(xi), RMAT(sol)) {
@@ -492,18 +509,82 @@ int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*,
     OK
 }
 
+//---------------------------------------------------------------
 
-int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
-{
-    gsl_function F;
-    F.function = f;
-    F.params = 0;
+typedef void TrawfunV(int, double*, double*);
 
-    if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
+int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
+    TrawfunV * f = (TrawfunV*) pars;
+    double* p = (double*)calloc(x->size,sizeof(double));
+    double* q = (double*)calloc(x->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+    f(x->size,p,q);
+    for(k=0;k<y->size;k++) {
+        gsl_vector_set(y,k,q[k]);
+    }
+    free(p);
+    free(q);
+    return 0; //hmmm
+}
 
-    if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
+int root(int method, void f(int, double*, int, double*),
+         double epsabs, int maxit,
+         KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
+    DEBUGMSG("root_only_f");
+    gsl_multiroot_function my_func;
+    // extract function from pars
+    my_func.f = only_f_aux_root;
+    my_func.n = xin;
+    my_func.params = f;
+    size_t iter = 0;
+    int status;
+    const gsl_multiroot_fsolver_type *T;
+    gsl_multiroot_fsolver *s;
+    // Starting point
+    KDVVIEW(xi);
+    switch(method) {
+        case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
+        case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
+        case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
+        case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_multiroot_fsolver_alloc (T, my_func.n);
+    gsl_multiroot_fsolver_set (s, &my_func, V(xi));
 
-    if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
+    do {
+           status = gsl_multiroot_fsolver_iterate (s);
 
-    return 0;
+           solp[iter*solc+0] = iter;
+
+           int k;
+           for(k=0;k<xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
+           }
+           for(k=xin;k<2*xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
+           }
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+             gsl_multiroot_test_residual (s->f, epsabs);
+        }
+        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_multiroot_fsolver_free(s);
+    OK
 }
diff --git a/lib/Numeric/GSL/gsl-aux.h b/lib/Numeric/GSL/gsl-aux.h
index e88322c..c9fd546 100644
--- a/lib/Numeric/GSL/gsl-aux.h
+++ b/lib/Numeric/GSL/gsl-aux.h
@@ -28,6 +28,8 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r));
 
 int fft(int code, KCVEC(a), CVEC(b));
 
+int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr);
+
 int integrate_qng(double f(double, void*), double a, double b, double prec,
                    double *result, double*error);
 
@@ -43,4 +45,6 @@ int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*,
                       double initstep, double minimpar, double tolgrad, int maxit, 
                       KRVEC(xi), RMAT(sol));
 
-int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr);
+int root(int method, void f(int, double*, int, double*),
+         double epsabs, int maxit,
+         KRVEC(xi), RMAT(sol));
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 278df78..4f73e3a 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -36,6 +36,8 @@ a ^ b = a Prelude.^ (b :: Int)
 
 utest str b = TestCase $ assertBool str b
 
+a ~~ b = fromList a |~| fromList b
+
 feye n = flipud (ident n) :: Matrix Double
 
 detTest1 = det m == 26
@@ -112,12 +114,17 @@ minimizationTest = TestList [ utest "minimization conj grad" (minim1 f df [5,7]
                             ]
     where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
           df [x,y] = [20*(x-1), 40*(y-2)]
-          a ~~ b = fromList a |~| fromList b
           minim1 g dg ini = fst $ minimizeConjugateGradient 1E-2 1E-4 1E-3 30 g dg ini
           minim2 g dg ini = fst $ minimizeVectorBFGS2 1E-2 1E-2 1E-3 30 g dg ini
 
 ---------------------------------------------------------------------
 
+rootFindingTest = utest "root Hybrids" (sol ~~ [1,1])
+    where sol = fst $ root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
+          rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+
+---------------------------------------------------------------------
+
 rot :: Double -> Matrix Double
 rot a = (3><3) [ c,0,s
                , 0,1,0
@@ -217,6 +224,7 @@ runTests n = do
         , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)
         , utest "polySolve" (polySolveProp [1,2,3,4])
         , minimizationTest
+        , rootFindingTest
         ]
     return ()