Merge pull request #39 from alang9/master

One-dimensional optimization methods

3 files changed, 215 insertions, 13 deletions

diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs
index 84417ce..69a9750 100644
--- a/lib/Numeric/GSL/Internal.hs
+++ b/lib/Numeric/GSL/Internal.hs
@@ -36,6 +36,9 @@ foreign import ccall safe "wrapper"
 foreign import ccall safe "wrapper"
     mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
 
+foreign import ccall safe "wrapper"
+    mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double))
+
 aux_vTov :: (Vector Double -> Vector Double) -> TVV
 aux_vTov f n p nr r = g where
     v = f x
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index cd2982a..6da15e5 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -45,6 +45,8 @@ main = do
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.Root (
+    uniRoot, UniRootMethod(..),
+    uniRootJ, UniRootMethodJ(..),
     root, RootMethod(..),
     rootJ, RootMethodJ(..),
 ) where
@@ -58,6 +60,68 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
 
+data UniRootMethod = Bisection
+                   | FalsePos
+                   | Brent
+                   deriving (Enum, Eq, Show, Bounded)
+
+uniRoot :: UniRootMethod
+        -> Double
+        -> Int
+        -> (Double -> Double)
+        -> Double
+        -> Double
+        -> (Double, Matrix Double)
+uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
+
+uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    rawpath <- createMIO maxit 4
+                         (c_root m fp epsrel (fi maxit) xl xu)
+                         "root"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+foreign import ccall safe "root"
+    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM
+
+-------------------------------------------------------------------------
+data UniRootMethodJ = UNewton
+                    | Secant
+                    | Steffenson
+                    deriving (Enum, Eq, Show, Bounded)
+
+uniRootJ :: UniRootMethodJ
+        -> Double
+        -> Int
+        -> (Double -> Double)
+        -> (Double -> Double)
+        -> Double
+        -> (Double, Matrix Double)
+uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
+    dfun x epsrel maxit
+
+uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    dfp <- mkDoublefun df
+    rawpath <- createMIO maxit 2
+                         (c_rootj m fp dfp epsrel (fi maxit) x)
+                         "rootj"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+foreign import ccall safe "rootj"
+    c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
+            -> Double -> CInt -> Double -> TM
+
+-------------------------------------------------------------------------
+
 data RootMethod = Hybrids
                 | Hybrid
                 | DNewton
@@ -82,8 +146,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
-                         (c_root m fp epsabs (fi maxit) // xiv')
-                         "root"
+                         (c_multiroot m fp epsabs (fi maxit) // xiv')
+                         "multiroot"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
@@ -91,8 +155,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
     return (take n $ drop 1 sol, path)
 
 
-foreign import ccall safe "root"
-    c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
+foreign import ccall safe "multiroot"
+    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
 
 -------------------------------------------------------------------------
 
@@ -120,8 +184,8 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
     jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
-                         (c_rootj m fp jp epsabs (fi maxit) // xiv')
-                         "root"
+                         (c_multirootj m fp jp epsabs (fi maxit) // xiv')
+                         "multiroot"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
@@ -129,9 +193,8 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
     freeHaskellFunPtr jp
     return (take n $ drop 1 sol, path)
 
-
-foreign import ccall safe "rootj"
-    c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
+foreign import ccall safe "multirootj"
+    c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
 
 -------------------------------------------------------
 
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 61b6a54..e727c91 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -32,6 +32,7 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_roots.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <string.h>
 #include <stdio.h>
@@ -1047,9 +1048,144 @@ int minimizeD(int method, double f(int, double*), int df(int, double*, int, doub
 
 //---------------------------------------------------------------
 
+double only_f_aux_root(double x, void *pars) {
+    double (*f)(double) = (double (*)(double)) pars;
+    return f(x);
+}
+
+int root(int method, double f(double),
+         double epsrel, int maxit,
+         double xl, double xu, RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
+    DEBUGMSG("root_only_f");
+    gsl_function my_func;
+    // extract function from pars
+    my_func.function = only_f_aux_root;
+    my_func.params = f;
+    size_t iter = 0;
+    int status;
+    const gsl_root_fsolver_type *T;
+    gsl_root_fsolver *s;
+    // Starting point
+    switch(method) {
+        case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
+        case 1 : {T = gsl_root_fsolver_falsepos; break; }
+        case 2 : {T = gsl_root_fsolver_brent; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_root_fsolver_alloc (T);
+    gsl_root_fsolver_set (s, &my_func, xl, xu);
+    do {
+           double best, current_lo, current_hi;
+           status = gsl_root_fsolver_iterate (s);
+           best = gsl_root_fsolver_root (s);
+           current_lo = gsl_root_fsolver_x_lower (s);
+           current_hi = gsl_root_fsolver_x_upper (s);
+           solp[iter*solc] = iter + 1;
+           solp[iter*solc+1] = best;
+           solp[iter*solc+2] = current_lo;
+           solp[iter*solc+3] = current_hi;
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+               gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+    int i;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        solp[i*solc+1]=0.;
+        solp[i*solc+2]=0.;
+        solp[i*solc+3]=0.;
+    }
+    gsl_root_fsolver_free(s);
+    OK
+}
+
+typedef struct {
+    double (*f)(double);
+    double (*jf)(double);
+} uniTfjf;
+
+double f_aux_uni(double x, void *pars) {
+    uniTfjf * fjf = ((uniTfjf*) pars);
+    return (fjf->f)(x);
+}
+
+double jf_aux_uni(double x, void * pars) {
+    uniTfjf * fjf = ((uniTfjf*) pars);
+    return (fjf->jf)(x);
+}
+
+void fjf_aux_uni(double x, void * pars, double * f, double * g) {
+    *f = f_aux_uni(x,pars);
+    *g = jf_aux_uni(x,pars);
+}
+
+int rootj(int method, double f(double),
+          double df(double),
+         double epsrel, int maxit,
+         double x, RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
+    DEBUGMSG("root_fjf");
+    gsl_function_fdf my_func;
+    // extract function from pars
+    my_func.f = f_aux_uni;
+    my_func.df = jf_aux_uni;
+    my_func.fdf = fjf_aux_uni;
+    uniTfjf stfjf;
+    stfjf.f = f;
+    stfjf.jf = df;
+    my_func.params = &stfjf;
+    size_t iter = 0;
+    int status;
+    const gsl_root_fdfsolver_type *T;
+    gsl_root_fdfsolver *s;
+    // Starting point
+    switch(method) {
+        case 0 : {T = gsl_root_fdfsolver_newton;; break; }
+        case 1 : {T = gsl_root_fdfsolver_secant; break; }
+        case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_root_fdfsolver_alloc (T);
+
+    gsl_root_fdfsolver_set (s, &my_func, x);
+
+    do {
+           double x0;
+           status = gsl_root_fdfsolver_iterate (s);
+           x0 = x;
+           x = gsl_root_fdfsolver_root(s);
+           solp[iter*solc+0] = iter+1;
+           solp[iter*solc+1] = x;
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+               gsl_root_test_delta (x, x0, 0, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+
+    int i;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        solp[i*solc+1]=0.;
+    }
+    gsl_root_fdfsolver_free(s);
+    OK
+}
+
+
+//---------------------------------------------------------------
+
 typedef void TrawfunV(int, double*, int, double*);
 
-int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
+int only_f_aux_multiroot(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(y->size,sizeof(double));
@@ -1066,14 +1202,14 @@ int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
     return 0; //hmmm
 }
 
-int root(int method, void f(int, double*, int, double*),
+int multiroot(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.f = only_f_aux_multiroot;
     my_func.n = xin;
     my_func.params = f;
     size_t iter = 0;
@@ -1175,7 +1311,7 @@ int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
     return 0;
 }
 
-int rootj(int method, int f(int, double*, int, double*),
+int multirootj(int method, int f(int, double*, int, double*),
                       int jac(int, double*, int, int, double*),
          double epsabs, int maxit,
          KRVEC(xi), RMAT(sol)) {