diff options
Diffstat (limited to 'lib/Numeric/GSL')
-rw-r--r-- | lib/Numeric/GSL/Internal.hs | 64 | ||||
-rw-r--r-- | lib/Numeric/GSL/Minimization.hs | 61 | ||||
-rw-r--r-- | lib/Numeric/GSL/Root.hs | 49 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 8 |
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 | |||
13 | module Numeric.GSL.Internal where | ||
14 | |||
15 | import Data.Packed.Internal | ||
16 | import Foreign | ||
17 | import Foreign.C.Types(CInt) | ||
18 | |||
19 | iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) | ||
20 | iv 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 | ||
26 | foreign import ccall "wrapper" | ||
27 | mkVecfun :: (CInt -> Ptr Double -> Double) | ||
28 | -> IO( FunPtr (CInt -> Ptr Double -> Double)) | ||
29 | |||
30 | foreign import ccall "wrapper" | ||
31 | mkVecVecfun :: TVV -> IO (FunPtr TVV) | ||
32 | |||
33 | aux_vTov :: (Vector Double -> Vector Double) -> TVV | ||
34 | aux_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 | |||
43 | foreign import ccall "wrapper" | ||
44 | mkVecMatfun :: TVM -> IO (FunPtr TVM) | ||
45 | |||
46 | aux_vTom :: (Vector Double -> Matrix Double) -> TVM | ||
47 | aux_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 | |||
56 | createV n fun msg = unsafePerformIO $ do | ||
57 | r <- createVector n | ||
58 | app1 fun vec r msg | ||
59 | return r | ||
60 | |||
61 | createMIO 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 | |||
63 | import Data.Packed.Matrix | 63 | import Data.Packed.Matrix |
64 | import Foreign | 64 | import Foreign |
65 | import Foreign.C.Types(CInt) | 65 | import Foreign.C.Types(CInt) |
66 | import 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 | ||
80 | data MinimizeMethod = NMSimplex | 81 | data 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. |
85 | minimize :: MinimizeMethod | 86 | minimize :: 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. |
103 | minimizeD :: MinimizeMethodD | 104 | minimizeD :: 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 | ||
160 | foreign import ccall "gsl-aux.h minimizeD" | 161 | foreign 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 | --------------------------------------------------------------------- |
167 | iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) | 169 | |
168 | iv f n p = f (createV (fromIntegral n) copy "iv") where | 170 | checkdim1 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 | ||
174 | foreign import ccall "wrapper" | ||
175 | mkVecfun :: (CInt -> Ptr Double -> Double) | ||
176 | -> IO( FunPtr (CInt -> Ptr Double -> Double)) | ||
177 | |||
178 | -- | another required conversion | ||
179 | foreign import ccall "wrapper" | ||
180 | mkVecVecfun :: (CInt -> Ptr Double -> Ptr Double -> IO ()) | ||
181 | -> IO (FunPtr (CInt -> Ptr Double -> Ptr Double->IO())) | ||
182 | |||
183 | aux_vTov :: (Vector Double -> Vector Double) -> (CInt -> Ptr Double -> Ptr Double -> IO()) | ||
184 | aux_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 | |||
194 | createV n fun msg = unsafePerformIO $ do | ||
195 | r <- createVector n | ||
196 | app1 fun vec r msg | ||
197 | return r | ||
198 | |||
199 | createMIO 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 | |||
53 | import Data.Packed.Matrix | 53 | import Data.Packed.Matrix |
54 | import Foreign | 54 | import Foreign |
55 | import Foreign.C.Types(CInt) | 55 | import Foreign.C.Types(CInt) |
56 | import 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. |
104 | rootJ :: RootMethodJ | 105 | rootJ :: 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 | ||
130 | foreign import ccall "rootj" | 132 | foreign 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 | |||
136 | foreign import ccall "wrapper" | ||
137 | mkVecVecfun :: TVV -> IO (FunPtr TVV) | ||
138 | |||
139 | aux_vTov :: (Vector Double -> Vector Double) -> TVV | ||
140 | aux_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 | |||
149 | foreign import ccall "wrapper" | ||
150 | mkVecMatfun :: TVM -> IO (FunPtr TVM) | ||
151 | |||
152 | aux_vTom :: (Vector Double -> Matrix Double) -> TVM | ||
153 | aux_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 | |||
162 | createV n fun msg = unsafePerformIO $ do | ||
163 | r <- createVector n | ||
164 | app1 fun vec r msg | ||
165 | return r | ||
166 | |||
167 | createMIO r c fun msg = do | ||
168 | res <- createMatrix RowMajor r c | ||
169 | app1 fun mat res msg | ||
170 | return res | ||
171 | 136 | ||
172 | checkdim1 n v | 137 | checkdim1 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 | ||
177 | checkdim2 n m | 142 | checkdim2 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 | ||
426 | typedef struct {double (*f)(int, double*); void (*df)(int, double*, double*);} Tfdf; | 426 | typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf; |
427 | 427 | ||
428 | double f_aux_min(const gsl_vector*x, void *pars) { | 428 | double 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) { | |||
441 | void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { | 441 | void 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 | ||
465 | int minimizeD(int method, double f(int, double*), void df(int, double*, double*), | 465 | int 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); |