From 6e0dd472ef8c570ec1924ac641e5872db30ac142 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 4 Jun 2009 09:01:56 +0000 Subject: added some root finding algorithms --- lib/Numeric/GSL.hs | 2 + lib/Numeric/GSL/Root.hs | 117 +++++++++++++++++++++++++++++++++++++ lib/Numeric/GSL/gsl-aux.c | 101 ++++++++++++++++++++++++++++---- lib/Numeric/GSL/gsl-aux.h | 6 +- lib/Numeric/LinearAlgebra/Tests.hs | 10 +++- 5 files changed, 224 insertions(+), 12 deletions(-) create mode 100644 lib/Numeric/GSL/Root.hs (limited to 'lib') 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. + + + +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 #include #include +#include #include #include #include @@ -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;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + f(x->size,p,q); + for(k=0;ksize;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;kx,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 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 () -- cgit v1.2.3