auxiliary functions moved to Numeric.GSL.Internal

8 files changed, 116 insertions, 139 deletions

diff --git a/INSTALL b/INSTALL
index b0e1970..f3bec7f 100644
--- a/INSTALL
+++ b/INSTALL
@@ -40,10 +40,10 @@ INSTALLATION ON WINDOWS ----------------------------------------
      line 17:      build-type:         Custom
      change to:    build-type:         Simple
 
-     line 160:     extra-libraries:
+     line 162:     extra-libraries:
      add:          extra-libraries: libgsl-0 blas lapack
 
-     line 161:     extra-lib-dirs:
+     line 163:     extra-lib-dirs:
      add:          extra-lib-dirs: c:\ghc\ghc-6.10.3\bin
 
 8) Open a terminal, cd to the hmatrix folder, and run
diff --git a/examples/minimize.hs b/examples/minimize.hs
index 11643c9..19b2cb3 100644
--- a/examples/minimize.hs
+++ b/examples/minimize.hs
@@ -4,20 +4,14 @@ import Numeric.LinearAlgebra
 import Graphics.Plot
 import Text.Printf(printf)
 
--- the function to be minimized 
+-- the function to be minimized
 f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
 
--- its gradient
+-- exact gradient
 df [x,y] = [20*(x-1), 40*(y-2)]
 
--- the conjugate gradient method
-minimizeCG = minimizeConjugateGradient 1E-2 1E-4 1E-3 30
-
--- the BFGS2 method
-minimizeBFGS2 = minimizeVectorBFGS2 1E-2 1E-2 1E-3 30
-
 -- a minimization algorithm which does not require the gradient
-minimizeS f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1E-2 100
+minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi
 
 -- Numerical estimation of the gradient
 gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]
@@ -26,47 +20,31 @@ partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
     g x = f (concat [a,x:b])
     (a,_:b) = splitAt n v
 
-main = do
-    putStrLn "BFGS2 with true gradient"
-    let (s,p) = minimizeBFGS2 f df [5,7]
-    print s -- solution
-    disp p -- evolution of the algorithm
-    let [x,y] = drop 2 (toColumns p)
-    mplot [x,y]  -- path from the starting point to the solution
-
-    putStrLn "conjugate gradient with true gradient"
-    let (s,p) = minimizeCG f df [5,7]
-    print s
-    disp p
-    let [x,y] = drop 2 (toColumns p)
-    mplot [x,y]
-
-    putStrLn "conjugate gradient with estimated gradient"
-    let (s,p) = minimizeCG f (gradient f) [5,7]
-    print s
-    disp p
-    mplot $ drop 2 (toColumns p)
-
-    putStrLn "without gradient, using the NM Simplex method"
-    let (s,p) = minimizeS f [5,7]
-    print s
-    disp p
-    mplot $ drop 3 (toColumns p)
-
-    putStrLn "-------------------------"
-    mapM_ test [NMSimplex,NMSimplex2]
-    mapM_ testd [ConjugateFR .. SteepestDescent]
-
 disp = putStrLn . format "  " (printf "%.3f")
 
+allMethods :: (Enum a, Bounded a) => [a]
+allMethods = [minBound .. maxBound]
+
 test method = do
     print method
     let (s,p) = minimize method 1E-2 30 [1,1] f [5,7]
     print s
     disp p
 
-testd method = do
+testD method = do
     print method
     let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7]
     print s
     disp p
+
+testD' method = do
+    putStrLn $ show method ++ " with estimated gradient"
+    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7]
+    print s
+    disp p
+
+main = do
+    mapM_ test [NMSimplex, NMSimplex2]
+    mapM_ testD allMethods
+    testD' ConjugateFR
+    mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7]
diff --git a/examples/root.hs b/examples/root.hs
index 2a24f0f..8546ff5 100644
--- a/examples/root.hs
+++ b/examples/root.hs
@@ -28,7 +28,4 @@ main = do
     test DNewton
     test Broyden
 
-    testJ HybridsJ
-    testJ HybridJ
-    testJ Newton
-    testJ GNewton
+    mapM_ testJ [HybridsJ .. GNewton]
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 7645d79..e75c654 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -31,6 +31,7 @@ extra-source-files: examples/tests.hs
                     examples/data.txt
                     examples/lie.hs
                     examples/kalman.hs
+                    examples/parallel.hs
                     examples/plot.hs
                     examples/latexmat.hs
                     examples/inplace.hs
@@ -64,7 +65,7 @@ flag unsafe
 
 library
     if flag(splitBase)
-      build-depends:    base >= 3, array, QuickCheck, HUnit, storable-complex, process
+      build-depends:    base >= 3 && < 5, array, QuickCheck, HUnit, storable-complex, process
     else
       build-depends:    base < 3, QuickCheck, HUnit, storable-complex, process
 
@@ -126,6 +127,7 @@ library
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix,
+                        Numeric.GSL.Internal,
                         Numeric.GSL.Special.Internal,
                         Numeric.LinearAlgebra.Tests.Instances,
                         Numeric.LinearAlgebra.Tests.Properties
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;k<x->size;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;k<x->size;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);