From d9efdd9334da1a63f739d6e2e68c4ff78f52e505 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 8 Jun 2009 09:45:14 +0000 Subject: auxiliary functions moved to Numeric.GSL.Internal --- lib/Numeric/GSL/Internal.hs | 64 +++++++++++++++++++++++++++++++++++++++++ lib/Numeric/GSL/Minimization.hs | 61 +++++++++++---------------------------- lib/Numeric/GSL/Root.hs | 49 +++++-------------------------- lib/Numeric/GSL/gsl-aux.c | 8 +++--- 4 files changed, 91 insertions(+), 91 deletions(-) create mode 100644 lib/Numeric/GSL/Internal.hs (limited to 'lib') diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs new file mode 100644 index 0000000..834dfc2 --- /dev/null +++ b/lib/Numeric/GSL/Internal.hs @@ -0,0 +1,64 @@ +-- Module : Numeric.GSL.Internal +-- Copyright : (c) Alberto Ruiz 2009 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz (aruiz at um dot es) +-- Stability : provisional +-- Portability : uses ffi +-- +-- Auxiliary functions. +-- +-- #hide + +module Numeric.GSL.Internal where + +import Data.Packed.Internal +import Foreign +import Foreign.C.Types(CInt) + +iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) +iv f n p = f (createV (fromIntegral n) copy "iv") where + copy n' q = do + copyArray q p (fromIntegral n') + return 0 + +-- | conversion of Haskell functions into function pointers that can be used in the C side +foreign import ccall "wrapper" + mkVecfun :: (CInt -> Ptr Double -> Double) + -> IO( FunPtr (CInt -> Ptr Double -> Double)) + +foreign import ccall "wrapper" + mkVecVecfun :: TVV -> IO (FunPtr TVV) + +aux_vTov :: (Vector Double -> Vector Double) -> TVV +aux_vTov f n p nr 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 = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr) + return 0 + +foreign import ccall "wrapper" + mkVecMatfun :: TVM -> IO (FunPtr TVM) + +aux_vTom :: (Vector Double -> Matrix Double) -> TVM +aux_vTom f n p rr cr r = g where + V {fptr = pr} = flatten $ f x + x = createV (fromIntegral n) copy "aux_vTov" + copy n' q = do + copyArray q p (fromIntegral n') + return 0 + g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr) + return 0 + +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/Minimization.hs b/lib/Numeric/GSL/Minimization.hs index 048c717..930fe8a 100644 --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs @@ -63,6 +63,7 @@ import Data.Packed.Internal import Data.Packed.Matrix import Foreign import Foreign.C.Types(CInt) +import Numeric.GSL.Internal ------------------------------------------------------------------------ @@ -79,7 +80,7 @@ minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit data MinimizeMethod = NMSimplex | NMSimplex2 - deriving (Enum,Eq,Show) + deriving (Enum,Eq,Show,Bounded) -- | Minimization without derivatives. minimize :: MinimizeMethod @@ -97,7 +98,7 @@ data MinimizeMethodD = ConjugateFR | VectorBFGS | VectorBFGS2 | SteepestDescent - deriving (Enum,Eq,Show) + deriving (Enum,Eq,Show,Bounded) -- | Minimization with derivatives. minimizeD :: MinimizeMethodD @@ -143,13 +144,13 @@ minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do let xiv = fromList xi n = dim xiv f' = f . toList - df' = (fromList . df . toList) + df' = (checkdim1 n .fromList . df . toList) fp <- mkVecfun (iv f') dfp <- mkVecVecfun (aux_vTov df') rawpath <- withVector xiv $ \xiv' -> createMIO maxit (n+2) - (c_minimizeWithDeriv method fp dfp istep tol eps (fi maxit) // xiv') - "minimizeDerivV" + (c_minimizeD method fp dfp istep tol eps (fi maxit) // xiv') + "minimizeD" let it = round (rawpath @@> (maxit-1,0)) path = takeRows it rawpath sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path @@ -158,45 +159,15 @@ minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do return (sol,path) foreign import ccall "gsl-aux.h minimizeD" - c_minimizeWithDeriv :: CInt -> FunPtr (CInt -> Ptr Double -> Double) - -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ()) - -> Double -> Double -> Double -> CInt - -> TVM + c_minimizeD :: CInt + -> FunPtr (CInt -> Ptr Double -> Double) + -> FunPtr TVV + -> Double -> Double -> Double -> CInt + -> TVM --------------------------------------------------------------------- -iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) -iv f n p = f (createV (fromIntegral n) copy "iv") where - copy n' q = do - copyArray q p (fromIntegral n') - return 0 - --- | conversion of Haskell functions into function pointers that can be used in the C side -foreign import ccall "wrapper" - mkVecfun :: (CInt -> Ptr Double -> Double) - -> IO( FunPtr (CInt -> Ptr Double -> Double)) - --- | another required conversion -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 + +checkdim1 n v + | dim v == n = v + | otherwise = error $ "Error: "++ show n + ++ " components expected in the result of the gradient supplied to minimizeD" diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs index 6ce2c4c..840e8ee 100644 --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs @@ -53,6 +53,7 @@ import Data.Packed.Internal import Data.Packed.Matrix import Foreign import Foreign.C.Types(CInt) +import Numeric.GSL.Internal ------------------------------------------------------------------------- @@ -60,7 +61,7 @@ data RootMethod = Hybrids | Hybrid | DNewton | Broyden - deriving (Enum,Eq,Show) + deriving (Enum,Eq,Show,Bounded) -- | Nonlinear multidimensional root finding using algorithms that do not require -- any derivative information to be supplied by the user. @@ -98,7 +99,7 @@ data RootMethodJ = HybridsJ | HybridJ | Newton | GNewton - deriving (Enum,Eq,Show) + deriving (Enum,Eq,Show,Bounded) -- | Nonlinear multidimensional root finding using both the function and its derivatives. rootJ :: RootMethodJ @@ -124,57 +125,21 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp + freeHaskellFunPtr jp return (take n $ drop 1 sol, path) foreign import ccall "rootj" c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM - ---------------------------------------------------------------------- - -foreign import ccall "wrapper" - mkVecVecfun :: TVV -> IO (FunPtr TVV) - -aux_vTov :: (Vector Double -> Vector Double) -> TVV -aux_vTov f n p nr 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 = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr) - return 0 - -foreign import ccall "wrapper" - mkVecMatfun :: TVM -> IO (FunPtr TVM) - -aux_vTom :: (Vector Double -> Matrix Double) -> TVM -aux_vTom f n p rr cr r = g where - V {fptr = pr} = flatten $ f x - x = createV (fromIntegral n) copy "aux_vTov" - copy n' q = do - copyArray q p (fromIntegral n') - return 0 - g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr) - return 0 - -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 +------------------------------------------------------- checkdim1 n v | dim v == n = v | otherwise = error $ "Error: "++ show n - ++ " results expected in the function supplied to root" + ++ " components expected in the result of the function supplied to root" checkdim2 n m | rows m == n && cols m == n = m | otherwise = error $ "Error: "++ show n ++ "x" ++ show n - ++ " Jacobian expected in root" + ++ " Jacobian expected in rootJ" diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 2ecfb51..9261f0c 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -423,7 +423,7 @@ int minimize(int method, double f(int, double*), double tolsize, int maxit, // working with the gradient -typedef struct {double (*f)(int, double*); void (*df)(int, double*, double*);} Tfdf; +typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf; double f_aux_min(const gsl_vector*x, void *pars) { Tfdf * fdf = ((Tfdf*) pars); @@ -441,13 +441,13 @@ double f_aux_min(const gsl_vector*x, void *pars) { void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { Tfdf * fdf = ((Tfdf*) pars); double* p = (double*)calloc(x->size,sizeof(double)); - double* q = (double*)calloc(x->size,sizeof(double)); + double* q = (double*)calloc(g->size,sizeof(double)); int k; for(k=0;ksize;k++) { p[k] = gsl_vector_get(x,k); } - fdf->df(x->size,p,q); + fdf->df(x->size,p,g->size,q); for(k=0;ksize;k++) { gsl_vector_set(g,k,q[k]); @@ -462,7 +462,7 @@ void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) } -int minimizeD(int method, double f(int, double*), void df(int, double*, double*), +int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*), double initstep, double minimpar, double tolgrad, int maxit, KRVEC(xi), RMAT(sol)) { REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); -- cgit v1.2.3