summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL')
-rw-r--r--lib/Numeric/GSL/Internal.hs64
-rw-r--r--lib/Numeric/GSL/Minimization.hs61
-rw-r--r--lib/Numeric/GSL/Root.hs49
-rw-r--r--lib/Numeric/GSL/gsl-aux.c8
4 files changed, 91 insertions, 91 deletions
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 @@
1-- Module : Numeric.GSL.Internal
2-- Copyright : (c) Alberto Ruiz 2009
3-- License : GPL
4--
5-- Maintainer : Alberto Ruiz (aruiz at um dot es)
6-- Stability : provisional
7-- Portability : uses ffi
8--
9-- Auxiliary functions.
10--
11-- #hide
12
13module Numeric.GSL.Internal where
14
15import Data.Packed.Internal
16import Foreign
17import Foreign.C.Types(CInt)
18
19iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
20iv f n p = f (createV (fromIntegral n) copy "iv") where
21 copy n' q = do
22 copyArray q p (fromIntegral n')
23 return 0
24
25-- | conversion of Haskell functions into function pointers that can be used in the C side
26foreign import ccall "wrapper"
27 mkVecfun :: (CInt -> Ptr Double -> Double)
28 -> IO( FunPtr (CInt -> Ptr Double -> Double))
29
30foreign import ccall "wrapper"
31 mkVecVecfun :: TVV -> IO (FunPtr TVV)
32
33aux_vTov :: (Vector Double -> Vector Double) -> TVV
34aux_vTov f n p nr r = g where
35 V {fptr = pr} = f x
36 x = createV (fromIntegral n) copy "aux_vTov"
37 copy n' q = do
38 copyArray q p (fromIntegral n')
39 return 0
40 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr)
41 return 0
42
43foreign import ccall "wrapper"
44 mkVecMatfun :: TVM -> IO (FunPtr TVM)
45
46aux_vTom :: (Vector Double -> Matrix Double) -> TVM
47aux_vTom f n p rr cr r = g where
48 V {fptr = pr} = flatten $ f x
49 x = createV (fromIntegral n) copy "aux_vTov"
50 copy n' q = do
51 copyArray q p (fromIntegral n')
52 return 0
53 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
54 return 0
55
56createV n fun msg = unsafePerformIO $ do
57 r <- createVector n
58 app1 fun vec r msg
59 return r
60
61createMIO r c fun msg = do
62 res <- createMatrix RowMajor r c
63 app1 fun mat res msg
64 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
63import Data.Packed.Matrix 63import Data.Packed.Matrix
64import Foreign 64import Foreign
65import Foreign.C.Types(CInt) 65import Foreign.C.Types(CInt)
66import Numeric.GSL.Internal
66 67
67------------------------------------------------------------------------ 68------------------------------------------------------------------------
68 69
@@ -79,7 +80,7 @@ minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit
79 80
80data MinimizeMethod = NMSimplex 81data MinimizeMethod = NMSimplex
81 | NMSimplex2 82 | NMSimplex2
82 deriving (Enum,Eq,Show) 83 deriving (Enum,Eq,Show,Bounded)
83 84
84-- | Minimization without derivatives. 85-- | Minimization without derivatives.
85minimize :: MinimizeMethod 86minimize :: MinimizeMethod
@@ -97,7 +98,7 @@ data MinimizeMethodD = ConjugateFR
97 | VectorBFGS 98 | VectorBFGS
98 | VectorBFGS2 99 | VectorBFGS2
99 | SteepestDescent 100 | SteepestDescent
100 deriving (Enum,Eq,Show) 101 deriving (Enum,Eq,Show,Bounded)
101 102
102-- | Minimization with derivatives. 103-- | Minimization with derivatives.
103minimizeD :: MinimizeMethodD 104minimizeD :: MinimizeMethodD
@@ -143,13 +144,13 @@ minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do
143 let xiv = fromList xi 144 let xiv = fromList xi
144 n = dim xiv 145 n = dim xiv
145 f' = f . toList 146 f' = f . toList
146 df' = (fromList . df . toList) 147 df' = (checkdim1 n .fromList . df . toList)
147 fp <- mkVecfun (iv f') 148 fp <- mkVecfun (iv f')
148 dfp <- mkVecVecfun (aux_vTov df') 149 dfp <- mkVecVecfun (aux_vTov df')
149 rawpath <- withVector xiv $ \xiv' -> 150 rawpath <- withVector xiv $ \xiv' ->
150 createMIO maxit (n+2) 151 createMIO maxit (n+2)
151 (c_minimizeWithDeriv method fp dfp istep tol eps (fi maxit) // xiv') 152 (c_minimizeD method fp dfp istep tol eps (fi maxit) // xiv')
152 "minimizeDerivV" 153 "minimizeD"
153 let it = round (rawpath @@> (maxit-1,0)) 154 let it = round (rawpath @@> (maxit-1,0))
154 path = takeRows it rawpath 155 path = takeRows it rawpath
155 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path 156 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
@@ -158,45 +159,15 @@ minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do
158 return (sol,path) 159 return (sol,path)
159 160
160foreign import ccall "gsl-aux.h minimizeD" 161foreign import ccall "gsl-aux.h minimizeD"
161 c_minimizeWithDeriv :: CInt -> FunPtr (CInt -> Ptr Double -> Double) 162 c_minimizeD :: CInt
162 -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ()) 163 -> FunPtr (CInt -> Ptr Double -> Double)
163 -> Double -> Double -> Double -> CInt 164 -> FunPtr TVV
164 -> TVM 165 -> Double -> Double -> Double -> CInt
166 -> TVM
165 167
166--------------------------------------------------------------------- 168---------------------------------------------------------------------
167iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) 169
168iv f n p = f (createV (fromIntegral n) copy "iv") where 170checkdim1 n v
169 copy n' q = do 171 | dim v == n = v
170 copyArray q p (fromIntegral n') 172 | otherwise = error $ "Error: "++ show n
171 return 0 173 ++ " components expected in the result of the gradient supplied to minimizeD"
172
173-- | conversion of Haskell functions into function pointers that can be used in the C side
174foreign import ccall "wrapper"
175 mkVecfun :: (CInt -> Ptr Double -> Double)
176 -> IO( FunPtr (CInt -> Ptr Double -> Double))
177
178-- | another required conversion
179foreign import ccall "wrapper"
180 mkVecVecfun :: (CInt -> Ptr Double -> Ptr Double -> IO ())
181 -> IO (FunPtr (CInt -> Ptr Double -> Ptr Double->IO()))
182
183aux_vTov :: (Vector Double -> Vector Double) -> (CInt -> Ptr Double -> Ptr Double -> IO())
184aux_vTov f n p r = g where
185 V {fptr = pr} = f x
186 x = createV (fromIntegral n) copy "aux_vTov"
187 copy n' q = do
188 copyArray q p (fromIntegral n')
189 return 0
190 g = withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral n)
191
192--------------------------------------------------------------------
193
194createV n fun msg = unsafePerformIO $ do
195 r <- createVector n
196 app1 fun vec r msg
197 return r
198
199createMIO r c fun msg = do
200 res <- createMatrix RowMajor r c
201 app1 fun mat res msg
202 return res
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
53import Data.Packed.Matrix 53import Data.Packed.Matrix
54import Foreign 54import Foreign
55import Foreign.C.Types(CInt) 55import Foreign.C.Types(CInt)
56import Numeric.GSL.Internal
56 57
57------------------------------------------------------------------------- 58-------------------------------------------------------------------------
58 59
@@ -60,7 +61,7 @@ data RootMethod = Hybrids
60 | Hybrid 61 | Hybrid
61 | DNewton 62 | DNewton
62 | Broyden 63 | Broyden
63 deriving (Enum,Eq,Show) 64 deriving (Enum,Eq,Show,Bounded)
64 65
65-- | Nonlinear multidimensional root finding using algorithms that do not require 66-- | Nonlinear multidimensional root finding using algorithms that do not require
66-- any derivative information to be supplied by the user. 67-- any derivative information to be supplied by the user.
@@ -98,7 +99,7 @@ data RootMethodJ = HybridsJ
98 | HybridJ 99 | HybridJ
99 | Newton 100 | Newton
100 | GNewton 101 | GNewton
101 deriving (Enum,Eq,Show) 102 deriving (Enum,Eq,Show,Bounded)
102 103
103-- | Nonlinear multidimensional root finding using both the function and its derivatives. 104-- | Nonlinear multidimensional root finding using both the function and its derivatives.
104rootJ :: RootMethodJ 105rootJ :: RootMethodJ
@@ -124,57 +125,21 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
124 path = takeRows it rawpath 125 path = takeRows it rawpath
125 [sol] = toLists $ dropRows (it-1) path 126 [sol] = toLists $ dropRows (it-1) path
126 freeHaskellFunPtr fp 127 freeHaskellFunPtr fp
128 freeHaskellFunPtr jp
127 return (take n $ drop 1 sol, path) 129 return (take n $ drop 1 sol, path)
128 130
129 131
130foreign import ccall "rootj" 132foreign import ccall "rootj"
131 c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM 133 c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
132 134
133 135-------------------------------------------------------
134---------------------------------------------------------------------
135
136foreign import ccall "wrapper"
137 mkVecVecfun :: TVV -> IO (FunPtr TVV)
138
139aux_vTov :: (Vector Double -> Vector Double) -> TVV
140aux_vTov f n p nr r = g where
141 V {fptr = pr} = f x
142 x = createV (fromIntegral n) copy "aux_vTov"
143 copy n' q = do
144 copyArray q p (fromIntegral n')
145 return 0
146 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr)
147 return 0
148
149foreign import ccall "wrapper"
150 mkVecMatfun :: TVM -> IO (FunPtr TVM)
151
152aux_vTom :: (Vector Double -> Matrix Double) -> TVM
153aux_vTom f n p rr cr r = g where
154 V {fptr = pr} = flatten $ f x
155 x = createV (fromIntegral n) copy "aux_vTov"
156 copy n' q = do
157 copyArray q p (fromIntegral n')
158 return 0
159 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
160 return 0
161
162createV n fun msg = unsafePerformIO $ do
163 r <- createVector n
164 app1 fun vec r msg
165 return r
166
167createMIO r c fun msg = do
168 res <- createMatrix RowMajor r c
169 app1 fun mat res msg
170 return res
171 136
172checkdim1 n v 137checkdim1 n v
173 | dim v == n = v 138 | dim v == n = v
174 | otherwise = error $ "Error: "++ show n 139 | otherwise = error $ "Error: "++ show n
175 ++ " results expected in the function supplied to root" 140 ++ " components expected in the result of the function supplied to root"
176 141
177checkdim2 n m 142checkdim2 n m
178 | rows m == n && cols m == n = m 143 | rows m == n && cols m == n = m
179 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n 144 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
180 ++ " Jacobian expected in root" 145 ++ " 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,
423 423
424// working with the gradient 424// working with the gradient
425 425
426typedef struct {double (*f)(int, double*); void (*df)(int, double*, double*);} Tfdf; 426typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
427 427
428double f_aux_min(const gsl_vector*x, void *pars) { 428double f_aux_min(const gsl_vector*x, void *pars) {
429 Tfdf * fdf = ((Tfdf*) pars); 429 Tfdf * fdf = ((Tfdf*) pars);
@@ -441,13 +441,13 @@ double f_aux_min(const gsl_vector*x, void *pars) {
441void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { 441void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
442 Tfdf * fdf = ((Tfdf*) pars); 442 Tfdf * fdf = ((Tfdf*) pars);
443 double* p = (double*)calloc(x->size,sizeof(double)); 443 double* p = (double*)calloc(x->size,sizeof(double));
444 double* q = (double*)calloc(x->size,sizeof(double)); 444 double* q = (double*)calloc(g->size,sizeof(double));
445 int k; 445 int k;
446 for(k=0;k<x->size;k++) { 446 for(k=0;k<x->size;k++) {
447 p[k] = gsl_vector_get(x,k); 447 p[k] = gsl_vector_get(x,k);
448 } 448 }
449 449
450 fdf->df(x->size,p,q); 450 fdf->df(x->size,p,g->size,q);
451 451
452 for(k=0;k<x->size;k++) { 452 for(k=0;k<x->size;k++) {
453 gsl_vector_set(g,k,q[k]); 453 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)
462} 462}
463 463
464 464
465int minimizeD(int method, double f(int, double*), void df(int, double*, double*), 465int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
466 double initstep, double minimpar, double tolgrad, int maxit, 466 double initstep, double minimpar, double tolgrad, int maxit,
467 KRVEC(xi), RMAT(sol)) { 467 KRVEC(xi), RMAT(sol)) {
468 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); 468 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);