summaryrefslogtreecommitdiff
path: root/packages/gsl/src/Numeric
diff options
context:
space:
mode:
authormaxc01 <xingchen92@gmail.com>2015-10-07 13:48:26 +0800
committermaxc01 <xingchen92@gmail.com>2015-10-07 13:48:26 +0800
commita61af756ddca4544de5e4969edc73131f4fccdd1 (patch)
tree2ac1755695a42d3964208e0029e74d446f5c3bd8 /packages/gsl/src/Numeric
parent0840304af1564fa86a6006d648450372f301a6c8 (diff)
parentc84a485f148063f6d0c23f016fe348ec94fb6b19 (diff)
Merge pull request #1 from albertoruiz/master
sync from albetoruiz/hmatrix
Diffstat (limited to 'packages/gsl/src/Numeric')
-rw-r--r--packages/gsl/src/Numeric/GSL/Fitting.hs18
-rw-r--r--packages/gsl/src/Numeric/GSL/Fourier.hs10
-rw-r--r--packages/gsl/src/Numeric/GSL/IO.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Internal.hs27
-rw-r--r--packages/gsl/src/Numeric/GSL/Interpolation.hs8
-rw-r--r--packages/gsl/src/Numeric/GSL/LinearAlgebra.hs14
-rw-r--r--packages/gsl/src/Numeric/GSL/Minimization.hs19
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs136
-rw-r--r--packages/gsl/src/Numeric/GSL/Polynomials.hs9
-rw-r--r--packages/gsl/src/Numeric/GSL/Random.hs16
-rw-r--r--packages/gsl/src/Numeric/GSL/Root.hs18
-rw-r--r--packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs245
-rw-r--r--packages/gsl/src/Numeric/GSL/Vector.hs15
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-aux.c27
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-ode.c31
15 files changed, 470 insertions, 125 deletions
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
index 0a92373..8eb93a7 100644
--- a/packages/gsl/src/Numeric/GSL/Fitting.hs
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE FlexibleContexts #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Fitting 4Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010 5Copyright : (c) Alberto Ruiz 2010
@@ -50,7 +52,7 @@ module Numeric.GSL.Fitting (
50 fitModelScaled, fitModel 52 fitModelScaled, fitModel
51) where 53) where
52 54
53import Numeric.LinearAlgebra 55import Numeric.LinearAlgebra.HMatrix
54import Numeric.GSL.Internal 56import Numeric.GSL.Internal
55 57
56import Foreign.Ptr(FunPtr, freeHaskellFunPtr) 58import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
@@ -80,13 +82,13 @@ nlFitting :: FittingMethod
80nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit 82nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
81 83
82nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do 84nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
83 let p = dim xiv 85 let p = size xiv
84 n = dim (f xiv) 86 n = size (f xiv)
85 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f)) 87 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
86 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac)) 88 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
87 rawpath <- createMatrix RowMajor maxit (2+p) 89 rawpath <- createMatrix RowMajor maxit (2+p)
88 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit" 90 c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n) # xiv # rawpath #|"c_nlfit"
89 let it = round (rawpath @@> (maxit-1,0)) 91 let it = round (rawpath `atIndex` (maxit-1,0))
90 path = takeRows it rawpath 92 path = takeRows it rawpath
91 [sol] = toRows $ dropRows (it-1) path 93 [sol] = toRows $ dropRows (it-1) path
92 freeHaskellFunPtr fp 94 freeHaskellFunPtr fp
@@ -99,7 +101,7 @@ foreign import ccall safe "nlfit"
99------------------------------------------------------- 101-------------------------------------------------------
100 102
101checkdim1 n _p v 103checkdim1 n _p v
102 | dim v == n = v 104 | size v == n = v
103 | otherwise = error $ "Error: "++ show n 105 | otherwise = error $ "Error: "++ show n
104 ++ " components expected in the result of the function supplied to nlFitting" 106 ++ " components expected in the result of the function supplied to nlFitting"
105 107
@@ -114,9 +116,9 @@ err (model,deriv) dat vsol = zip sol errs where
114 sol = toList vsol 116 sol = toList vsol
115 c = max 1 (chi/sqrt (fromIntegral dof)) 117 c = max 1 (chi/sqrt (fromIntegral dof))
116 dof = length dat - (rows cov) 118 dof = length dat - (rows cov)
117 chi = norm2 (fromList $ cost (resMs model) dat sol) 119 chi = norm_2 (fromList $ cost (resMs model) dat sol)
118 js = fromLists $ jacobian (resDs deriv) dat sol 120 js = fromLists $ jacobian (resDs deriv) dat sol
119 cov = inv $ trans js <> js 121 cov = inv $ tr js <> js
120 errs = toList $ scalar c * sqrt (takeDiag cov) 122 errs = toList $ scalar c * sqrt (takeDiag cov)
121 123
122 124
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
index 734325b..1c2c053 100644
--- a/packages/gsl/src/Numeric/GSL/Fourier.hs
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE TypeFamilies #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Fourier 4Module : Numeric.GSL.Fourier
3Copyright : (c) Alberto Ruiz 2006 5Copyright : (c) Alberto Ruiz 2006
@@ -16,15 +18,14 @@ module Numeric.GSL.Fourier (
16 ifft 18 ifft
17) where 19) where
18 20
19import Data.Packed 21import Numeric.LinearAlgebra.HMatrix
20import Numeric.GSL.Internal 22import Numeric.GSL.Internal
21import Data.Complex
22import Foreign.C.Types 23import Foreign.C.Types
23import System.IO.Unsafe (unsafePerformIO) 24import System.IO.Unsafe (unsafePerformIO)
24 25
25genfft code v = unsafePerformIO $ do 26genfft code v = unsafePerformIO $ do
26 r <- createVector (dim v) 27 r <- createVector (size v)
27 app2 (c_fft code) vec v vec r "fft" 28 c_fft code # v # r #|"fft"
28 return r 29 return r
29 30
30foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res) 31foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res)
@@ -42,3 +43,4 @@ fft = genfft 0
42-- | The inverse of 'fft', using /gsl_fft_complex_inverse/. 43-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
43ifft :: Vector (Complex Double) -> Vector (Complex Double) 44ifft :: Vector (Complex Double) -> Vector (Complex Double)
44ifft = genfft 1 45ifft = genfft 1
46
diff --git a/packages/gsl/src/Numeric/GSL/IO.hs b/packages/gsl/src/Numeric/GSL/IO.hs
index 0d6031a..936f6bf 100644
--- a/packages/gsl/src/Numeric/GSL/IO.hs
+++ b/packages/gsl/src/Numeric/GSL/IO.hs
@@ -14,7 +14,7 @@ module Numeric.GSL.IO (
14 fileDimensions, loadMatrix, fromFile 14 fileDimensions, loadMatrix, fromFile
15) where 15) where
16 16
17import Data.Packed 17import Numeric.LinearAlgebra.HMatrix hiding(saveMatrix, loadMatrix)
18import Numeric.GSL.Vector 18import Numeric.GSL.Vector
19import System.Process(readProcess) 19import System.Process(readProcess)
20 20
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
index a1c4e0c..dcd3bc4 100644
--- a/packages/gsl/src/Numeric/GSL/Internal.hs
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -22,21 +22,20 @@ module Numeric.GSL.Internal(
22 aux_vTom, 22 aux_vTom,
23 createV, 23 createV,
24 createMIO, 24 createMIO,
25 module Data.Packed.Development, 25 module Numeric.LinearAlgebra.Devel,
26 check, 26 check,(#),vec, ww2,
27 Res,TV,TM,TCV,TCM 27 Res,TV,TM,TCV,TCM
28) where 28) where
29 29
30import Data.Packed 30import Numeric.LinearAlgebra.HMatrix
31import Data.Packed.Development hiding (check) 31import Numeric.LinearAlgebra.Devel hiding (check)
32import Data.Complex
33 32
34import Foreign.Marshal.Array(copyArray) 33import Foreign.Marshal.Array(copyArray)
35import Foreign.Ptr(Ptr, FunPtr) 34import Foreign.Ptr(Ptr, FunPtr)
36import Foreign.C.Types 35import Foreign.C.Types
37import Foreign.C.String(peekCString) 36import Foreign.C.String(peekCString)
38import System.IO.Unsafe(unsafePerformIO) 37import System.IO.Unsafe(unsafePerformIO)
39import Data.Vector.Storable(unsafeWith) 38import Data.Vector.Storable as V (unsafeWith,length)
40import Control.Monad(when) 39import Control.Monad(when)
41 40
42iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) 41iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
@@ -87,12 +86,12 @@ aux_vTom f n p rr cr r = g where
87 86
88createV n fun msg = unsafePerformIO $ do 87createV n fun msg = unsafePerformIO $ do
89 r <- createVector n 88 r <- createVector n
90 app1 fun vec r msg 89 fun # r #| msg
91 return r 90 return r
92 91
93createMIO r c fun msg = do 92createMIO r c fun msg = do
94 res <- createMatrix RowMajor r c 93 res <- createMatrix RowMajor r c
95 app1 fun mat res msg 94 fun # res #| msg
96 return res 95 return res
97 96
98-------------------------------------------------------------------------------- 97--------------------------------------------------------------------------------
@@ -124,3 +123,15 @@ type TCM x = CInt -> CInt -> PC -> x
124type TVV = TV (TV Res) 123type TVV = TV (TV Res)
125type TVM = TV (TM Res) 124type TVM = TV (TM Res)
126 125
126ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
127
128vec x f = unsafeWith x $ \p -> do
129 let v g = do
130 g (fi $ V.length x) p
131 f v
132{-# INLINE vec #-}
133
134infixl 1 #
135a # b = applyRaw a b
136{-# INLINE (#) #-}
137
diff --git a/packages/gsl/src/Numeric/GSL/Interpolation.hs b/packages/gsl/src/Numeric/GSL/Interpolation.hs
index 4d72ee2..d060468 100644
--- a/packages/gsl/src/Numeric/GSL/Interpolation.hs
+++ b/packages/gsl/src/Numeric/GSL/Interpolation.hs
@@ -32,8 +32,7 @@ module Numeric.GSL.Interpolation (
32 , evaluateIntegralV 32 , evaluateIntegralV
33) where 33) where
34 34
35import Data.Packed.Vector(Vector, fromList, dim) 35import Numeric.LinearAlgebra(Vector, fromList, size, Numeric)
36import Data.Packed.Foreign(appVector)
37import Foreign.C.Types 36import Foreign.C.Types
38import Foreign.Marshal.Alloc(alloca) 37import Foreign.Marshal.Alloc(alloca)
39import Foreign.Ptr(Ptr) 38import Foreign.Ptr(Ptr)
@@ -57,6 +56,9 @@ methodToInt CSplinePeriodic = 3
57methodToInt Akima = 4 56methodToInt Akima = 4
58methodToInt AkimaPeriodic = 5 57methodToInt AkimaPeriodic = 5
59 58
59dim :: Numeric t => Vector t -> Int
60dim = size
61
60applyCFun hsname cname fun mth xs ys x 62applyCFun hsname cname fun mth xs ys x
61 | dim xs /= dim ys = error $ 63 | dim xs /= dim ys = error $
62 "Error: Vectors of unequal sizes " ++ 64 "Error: Vectors of unequal sizes " ++
@@ -115,7 +117,7 @@ evaluate :: InterpolationMethod -- ^ What method to use to interpolate
115 -> Double -- ^ Point at which to evaluate the function 117 -> Double -- ^ Point at which to evaluate the function
116 -> Double -- ^ Interpolated result 118 -> Double -- ^ Interpolated result
117evaluate mth pts = 119evaluate mth pts =
118 applyCFun "evaluate" "spline_eval" c_spline_eval_deriv 120 applyCFun "evaluate" "spline_eval" c_spline_eval
119 mth (fromList xs) (fromList ys) 121 mth (fromList xs) (fromList ys)
120 where 122 where
121 (xs, ys) = unzip pts 123 (xs, ys) = unzip pts
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
index 17e2258..6ffe306 100644
--- a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -15,7 +15,7 @@ module Numeric.GSL.LinearAlgebra (
15 fileDimensions, loadMatrix, fromFile 15 fileDimensions, loadMatrix, fromFile
16) where 16) where
17 17
18import Data.Packed 18import Numeric.LinearAlgebra.HMatrix hiding (RandDist,randomVector,saveMatrix,loadMatrix)
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20 20
21import Foreign.Marshal.Alloc(free) 21import Foreign.Marshal.Alloc(free)
@@ -40,7 +40,7 @@ randomVector :: Int -- ^ seed
40 -> Vector Double 40 -> Vector Double
41randomVector seed dist n = unsafePerformIO $ do 41randomVector seed dist n = unsafePerformIO $ do
42 r <- createVector n 42 r <- createVector n
43 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" 43 c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
44 return r 44 return r
45 45
46foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV 46foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
@@ -56,7 +56,7 @@ saveMatrix filename fmt m = do
56 charname <- newCString filename 56 charname <- newCString filename
57 charfmt <- newCString fmt 57 charfmt <- newCString fmt
58 let o = if orderOf m == RowMajor then 1 else 0 58 let o = if orderOf m == RowMajor then 1 else 0
59 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" 59 matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
60 free charname 60 free charname
61 free charfmt 61 free charfmt
62 62
@@ -69,7 +69,7 @@ fscanfVector :: FilePath -> Int -> IO (Vector Double)
69fscanfVector filename n = do 69fscanfVector filename n = do
70 charname <- newCString filename 70 charname <- newCString filename
71 res <- createVector n 71 res <- createVector n
72 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" 72 gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
73 free charname 73 free charname
74 return res 74 return res
75 75
@@ -80,7 +80,7 @@ fprintfVector :: FilePath -> String -> Vector Double -> IO ()
80fprintfVector filename fmt v = do 80fprintfVector filename fmt v = do
81 charname <- newCString filename 81 charname <- newCString filename
82 charfmt <- newCString fmt 82 charfmt <- newCString fmt
83 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" 83 gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
84 free charname 84 free charname
85 free charfmt 85 free charfmt
86 86
@@ -91,7 +91,7 @@ freadVector :: FilePath -> Int -> IO (Vector Double)
91freadVector filename n = do 91freadVector filename n = do
92 charname <- newCString filename 92 charname <- newCString filename
93 res <- createVector n 93 res <- createVector n
94 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" 94 gsl_vector_fread charname # res #| "gsl_vector_fread"
95 free charname 95 free charname
96 return res 96 return res
97 97
@@ -101,7 +101,7 @@ foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
101fwriteVector :: FilePath -> Vector Double -> IO () 101fwriteVector :: FilePath -> Vector Double -> IO ()
102fwriteVector filename v = do 102fwriteVector filename v = do
103 charname <- newCString filename 103 charname <- newCString filename
104 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" 104 gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
105 free charname 105 free charname
106 106
107foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV 107foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
index 056d463..a0e5306 100644
--- a/packages/gsl/src/Numeric/GSL/Minimization.hs
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -1,3 +1,6 @@
1{-# LANGUAGE FlexibleContexts #-}
2
3
1{- | 4{- |
2Module : Numeric.GSL.Minimization 5Module : Numeric.GSL.Minimization
3Copyright : (c) Alberto Ruiz 2006-9 6Copyright : (c) Alberto Ruiz 2006-9
@@ -56,7 +59,7 @@ module Numeric.GSL.Minimization (
56) where 59) where
57 60
58 61
59import Data.Packed 62import Numeric.LinearAlgebra.HMatrix hiding(step)
60import Numeric.GSL.Internal 63import Numeric.GSL.Internal
61 64
62import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) 65import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
@@ -99,7 +102,7 @@ uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
99 rawpath <- createMIO maxit 4 102 rawpath <- createMIO maxit 4
100 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu) 103 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
101 "uniMinimize" 104 "uniMinimize"
102 let it = round (rawpath @@> (maxit-1,0)) 105 let it = round (rawpath `atIndex` (maxit-1,0))
103 path = takeRows it rawpath 106 path = takeRows it rawpath
104 [sol] = toLists $ dropRows (it-1) path 107 [sol] = toLists $ dropRows (it-1) path
105 freeHaskellFunPtr fp 108 freeHaskellFunPtr fp
@@ -134,16 +137,16 @@ minimizeV :: MinimizeMethod
134minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi) 137minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
135 where v2l (v,m) = (toList v, m) 138 where v2l (v,m) = (toList v, m)
136 139
137ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 140
138 141
139minimizeV method eps maxit szv f xiv = unsafePerformIO $ do 142minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
140 let n = dim xiv 143 let n = size xiv
141 fp <- mkVecfun (iv f) 144 fp <- mkVecfun (iv f)
142 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' -> 145 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
143 createMIO maxit (n+3) 146 createMIO maxit (n+3)
144 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') 147 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
145 "minimize" 148 "minimize"
146 let it = round (rawpath @@> (maxit-1,0)) 149 let it = round (rawpath `atIndex` (maxit-1,0))
147 path = takeRows it rawpath 150 path = takeRows it rawpath
148 sol = flatten $ dropColumns 3 $ dropRows (it-1) path 151 sol = flatten $ dropColumns 3 $ dropRows (it-1) path
149 freeHaskellFunPtr fp 152 freeHaskellFunPtr fp
@@ -191,7 +194,7 @@ minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
191 194
192 195
193minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do 196minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
194 let n = dim xiv 197 let n = size xiv
195 f' = f 198 f' = f
196 df' = (checkdim1 n . df) 199 df' = (checkdim1 n . df)
197 fp <- mkVecfun (iv f') 200 fp <- mkVecfun (iv f')
@@ -200,7 +203,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
200 createMIO maxit (n+2) 203 createMIO maxit (n+2)
201 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') 204 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
202 "minimizeD" 205 "minimizeD"
203 let it = round (rawpath @@> (maxit-1,0)) 206 let it = round (rawpath `atIndex` (maxit-1,0))
204 path = takeRows it rawpath 207 path = takeRows it rawpath
205 sol = flatten $ dropColumns 2 $ dropRows (it-1) path 208 sol = flatten $ dropColumns 2 $ dropRows (it-1) path
206 freeHaskellFunPtr fp 209 freeHaskellFunPtr fp
@@ -217,6 +220,6 @@ foreign import ccall safe "gsl-aux.h minimizeD"
217--------------------------------------------------------------------- 220---------------------------------------------------------------------
218 221
219checkdim1 n v 222checkdim1 n v
220 | dim v == n = v 223 | size v == n = v
221 | otherwise = error $ "Error: "++ show n 224 | otherwise = error $ "Error: "++ show n
222 ++ " components expected in the result of the gradient supplied to minimizeD" 225 ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 7549a65..9e52873 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -1,3 +1,6 @@
1{-# LANGUAGE FlexibleContexts #-}
2
3
1{- | 4{- |
2Module : Numeric.GSL.ODE 5Module : Numeric.GSL.ODE
3Copyright : (c) Alberto Ruiz 2010 6Copyright : (c) Alberto Ruiz 2010
@@ -29,10 +32,10 @@ main = mplot (ts : toColumns sol)
29----------------------------------------------------------------------------- 32-----------------------------------------------------------------------------
30 33
31module Numeric.GSL.ODE ( 34module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..), Jacobian 35 odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..)
33) where 36) where
34 37
35import Data.Packed 38import Numeric.LinearAlgebra.HMatrix
36import Numeric.GSL.Internal 39import Numeric.GSL.Internal
37 40
38import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) 41import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
@@ -41,9 +44,10 @@ import System.IO.Unsafe(unsafePerformIO)
41 44
42------------------------------------------------------------------------- 45-------------------------------------------------------------------------
43 46
44type TVV = TV (TV Res) 47type TVV = TV (TV Res)
45type TVM = TV (TM Res) 48type TVM = TV (TM Res)
46type TVVM = TV (TV (TM Res)) 49type TVVM = TV (TV (TM Res))
50type TVVVM = TV (TV (TV (TM Res)))
47 51
48type Jacobian = Double -> Vector Double -> Matrix Double 52type Jacobian = Double -> Vector Double -> Matrix Double
49 53
@@ -60,73 +64,105 @@ data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
60 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. 64 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
61 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. 65 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
62 66
67-- | Adaptive step-size control functions
68data StepControl = X Double Double -- ^ abs. and rel. tolerance for x(t)
69 | X' Double Double -- ^ abs. and rel. tolerance for x'(t)
70 | XX' Double Double Double Double -- ^ include both via rel. tolerance scaling factors a_x, a_x'
71 | ScXX' Double Double Double Double (Vector Double) -- ^ scale abs. tolerance of x(t) components
63 72
64-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. 73-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
65odeSolve 74odeSolve
66 :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) 75 :: (Double -> [Double] -> [Double]) -- ^ x'(t,x)
67 -> [Double] -- ^ initial conditions 76 -> [Double] -- ^ initial conditions
68 -> Vector Double -- ^ desired solution times 77 -> Vector Double -- ^ desired solution times
69 -> Matrix Double -- ^ solution 78 -> Matrix Double -- ^ solution
70odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts 79odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
71 where hi = (ts@>1 - ts@>0)/100 80 where hi = (ts!1 - ts!0)/100
72 epsAbs = 1.49012e-08 81 epsAbs = 1.49012e-08
73 epsRel = 1.49012e-08 82 epsRel = epsAbs
74 l2v f = \t -> fromList . f t . toList 83 l2v f = \t -> fromList . f t . toList
75 84
76-- | Evolution of the system with adaptive step-size control. 85-- | A version of 'odeSolveVWith' with reasonable default step control.
77odeSolveV 86odeSolveV
78 :: ODEMethod 87 :: ODEMethod
79 -> Double -- ^ initial step size 88 -> Double -- ^ initial step size
80 -> Double -- ^ absolute tolerance for the state vector 89 -> Double -- ^ absolute tolerance for the state vector
81 -> Double -- ^ relative tolerance for the state vector 90 -> Double -- ^ relative tolerance for the state vector
82 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) 91 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
83 -> Vector Double -- ^ initial conditions 92 -> Vector Double -- ^ initial conditions
84 -> Vector Double -- ^ desired solution times 93 -> Vector Double -- ^ desired solution times
85 -> Matrix Double -- ^ solution 94 -> Matrix Double -- ^ solution
86odeSolveV RK2 = odeSolveV' 0 Nothing 95odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi
87odeSolveV RK4 = odeSolveV' 1 Nothing 96
88odeSolveV RKf45 = odeSolveV' 2 Nothing 97-- | Evolution of the system with adaptive step-size control.
89odeSolveV RKck = odeSolveV' 3 Nothing 98odeSolveVWith
90odeSolveV RK8pd = odeSolveV' 4 Nothing 99 :: ODEMethod
91odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) 100 -> StepControl
92odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) 101 -> Double -- ^ initial step size
93odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) 102 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
94odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
95odeSolveV MSAdams = odeSolveV' 9 Nothing
96odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
97
98
99odeSolveV'
100 :: CInt
101 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
102 -> Double -- ^ initial step size
103 -> Double -- ^ absolute tolerance for the state vector
104 -> Double -- ^ relative tolerance for the state vector
105 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
106 -> Vector Double -- ^ initial conditions 103 -> Vector Double -- ^ initial conditions
107 -> Vector Double -- ^ desired solution times 104 -> Vector Double -- ^ desired solution times
108 -> Matrix Double -- ^ solution 105 -> Matrix Double -- ^ solution
109odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do 106odeSolveVWith method control = odeSolveVWith' m mbj c epsAbs epsRel aX aX' mbsc
110 let n = dim xiv 107 where (m, mbj) = case method of
111 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) 108 RK2 -> (0 , Nothing )
112 jp <- case mbjac of 109 RK4 -> (1 , Nothing )
113 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) 110 RKf45 -> (2 , Nothing )
114 Nothing -> return nullFunPtr 111 RKck -> (3 , Nothing )
115 sol <- vec xiv $ \xiv' -> 112 RK8pd -> (4 , Nothing )
116 vec (checkTimes ts) $ \ts' -> 113 RK2imp jac -> (5 , Just jac)
117 createMIO (dim ts) n 114 RK4imp jac -> (6 , Just jac)
118 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) 115 BSimp jac -> (7 , Just jac)
119 "ode" 116 RK1imp jac -> (8 , Just jac)
120 freeHaskellFunPtr fp 117 MSAdams -> (9 , Nothing )
121 return sol 118 MSBDF jac -> (10, Just jac)
119 (c, epsAbs, epsRel, aX, aX', mbsc) = case control of
120 X ea er -> (0, ea, er, 1 , 0 , Nothing)
121 X' ea er -> (0, ea, er, 0 , 1 , Nothing)
122 XX' ea er ax ax' -> (0, ea, er, ax, ax', Nothing)
123 ScXX' ea er ax ax' sc -> (1, ea, er, ax, ax', Just sc)
124
125odeSolveVWith'
126 :: CInt -- ^ stepping function
127 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
128 -> CInt -- ^ step-size control function
129 -> Double -- ^ absolute tolerance for step-size control
130 -> Double -- ^ relative tolerance for step-size control
131 -> Double -- ^ scaling factor for relative tolerance of x(t)
132 -> Double -- ^ scaling factor for relative tolerance of x'(t)
133 -> Maybe (Vector Double) -- ^ optional scaling for absolute error
134 -> Double -- ^ initial step size
135 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
136 -> Vector Double -- ^ initial conditions
137 -> Vector Double -- ^ desired solution times
138 -> Matrix Double -- ^ solution
139odeSolveVWith' method mbjac control epsAbs epsRel aX aX' mbsc h f xiv ts =
140 unsafePerformIO $ do
141 let n = size xiv
142 sc = case mbsc of
143 Just scv -> checkdim1 n scv
144 Nothing -> xiv
145 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
146 jp <- case mbjac of
147 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
148 Nothing -> return nullFunPtr
149 sol <- vec sc $ \sc' -> vec xiv $ \xiv' ->
150 vec (checkTimes ts) $ \ts' -> createMIO (size ts) n
151 (ode_c method control h epsAbs epsRel aX aX' fp jp
152 // sc' // xiv' // ts' )
153 "ode"
154 freeHaskellFunPtr fp
155 return sol
122 156
123foreign import ccall safe "ode" 157foreign import ccall safe "ode"
124 ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM 158 ode_c :: CInt -> CInt -> Double
159 -> Double -> Double -> Double -> Double
160 -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVVM
125 161
126------------------------------------------------------- 162-------------------------------------------------------
127 163
128checkdim1 n v 164checkdim1 n v
129 | dim v == n = v 165 | size v == n = v
130 | otherwise = error $ "Error: "++ show n 166 | otherwise = error $ "Error: "++ show n
131 ++ " components expected in the result of the function supplied to odeSolve" 167 ++ " components expected in the result of the function supplied to odeSolve"
132 168
@@ -135,6 +171,6 @@ checkdim2 n m
135 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n 171 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
136 ++ " Jacobian expected in odeSolve" 172 ++ " Jacobian expected in odeSolve"
137 173
138checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts 174checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
139 | otherwise = error "odeSolve requires increasing times" 175 | otherwise = error "odeSolve requires increasing times"
140 where ts' = toList ts 176 where ts' = toList ts
diff --git a/packages/gsl/src/Numeric/GSL/Polynomials.hs b/packages/gsl/src/Numeric/GSL/Polynomials.hs
index b1be85d..8890f8f 100644
--- a/packages/gsl/src/Numeric/GSL/Polynomials.hs
+++ b/packages/gsl/src/Numeric/GSL/Polynomials.hs
@@ -16,9 +16,8 @@ module Numeric.GSL.Polynomials (
16 polySolve 16 polySolve
17) where 17) where
18 18
19import Data.Packed 19import Numeric.LinearAlgebra.HMatrix
20import Numeric.GSL.Internal 20import Numeric.GSL.Internal
21import Data.Complex
22import System.IO.Unsafe (unsafePerformIO) 21import System.IO.Unsafe (unsafePerformIO)
23 22
24#if __GLASGOW_HASKELL__ >= 704 23#if __GLASGOW_HASKELL__ >= 704
@@ -47,9 +46,9 @@ polySolve :: [Double] -> [Complex Double]
47polySolve = toList . polySolve' . fromList 46polySolve = toList . polySolve' . fromList
48 47
49polySolve' :: Vector Double -> Vector (Complex Double) 48polySolve' :: Vector Double -> Vector (Complex Double)
50polySolve' v | dim v > 1 = unsafePerformIO $ do 49polySolve' v | size v > 1 = unsafePerformIO $ do
51 r <- createVector (dim v-1) 50 r <- createVector (size v-1)
52 app2 c_polySolve vec v vec r "polySolve" 51 c_polySolve # v # r #| "polySolve"
53 return r 52 return r
54 | otherwise = error "polySolve on a polynomial of degree zero" 53 | otherwise = error "polySolve on a polynomial of degree zero"
55 54
diff --git a/packages/gsl/src/Numeric/GSL/Random.hs b/packages/gsl/src/Numeric/GSL/Random.hs
index f1f49e5..139c921 100644
--- a/packages/gsl/src/Numeric/GSL/Random.hs
+++ b/packages/gsl/src/Numeric/GSL/Random.hs
@@ -21,11 +21,13 @@ module Numeric.GSL.Random (
21) where 21) where
22 22
23import Numeric.GSL.Vector 23import Numeric.GSL.Vector
24import Numeric.LinearAlgebra(cholSH) 24import Numeric.LinearAlgebra.HMatrix hiding (
25import Numeric.Container hiding (
26 randomVector, 25 randomVector,
27 gaussianSample, 26 gaussianSample,
28 uniformSample 27 uniformSample,
28 Seed,
29 rand,
30 randn
29 ) 31 )
30import System.Random(randomIO) 32import System.Random(randomIO)
31 33
@@ -40,10 +42,10 @@ gaussianSample :: Seed
40 -> Matrix Double -- ^ covariance matrix 42 -> Matrix Double -- ^ covariance matrix
41 -> Matrix Double -- ^ result 43 -> Matrix Double -- ^ result
42gaussianSample seed n med cov = m where 44gaussianSample seed n med cov = m where
43 c = dim med 45 c = size med
44 meds = konst 1 n `outer` med 46 meds = konst 1 n `outer` med
45 rs = reshape c $ randomVector seed Gaussian (c * n) 47 rs = reshape c $ randomVector seed Gaussian (c * n)
46 m = rs `mXm` cholSH cov `add` meds 48 m = rs <> cholSH cov + meds
47 49
48-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate 50-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
49-- uniform distribution. 51-- uniform distribution.
@@ -55,10 +57,10 @@ uniformSample seed n rgs = m where
55 (as,bs) = unzip rgs 57 (as,bs) = unzip rgs
56 a = fromList as 58 a = fromList as
57 cs = zipWith subtract as bs 59 cs = zipWith subtract as bs
58 d = dim a 60 d = size a
59 dat = toRows $ reshape n $ randomVector seed Uniform (n*d) 61 dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
60 am = konst 1 n `outer` a 62 am = konst 1 n `outer` a
61 m = fromColumns (zipWith scale cs dat) `add` am 63 m = fromColumns (zipWith scale cs dat) + am
62 64
63-- | pseudorandom matrix with uniform elements between 0 and 1 65-- | pseudorandom matrix with uniform elements between 0 and 1
64randm :: RandDist 66randm :: RandDist
diff --git a/packages/gsl/src/Numeric/GSL/Root.hs b/packages/gsl/src/Numeric/GSL/Root.hs
index b9f3b94..724f32f 100644
--- a/packages/gsl/src/Numeric/GSL/Root.hs
+++ b/packages/gsl/src/Numeric/GSL/Root.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE FlexibleContexts #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Root 4Module : Numeric.GSL.Root
3Copyright : (c) Alberto Ruiz 2009 5Copyright : (c) Alberto Ruiz 2009
@@ -39,7 +41,7 @@ module Numeric.GSL.Root (
39 rootJ, RootMethodJ(..), 41 rootJ, RootMethodJ(..),
40) where 42) where
41 43
42import Data.Packed 44import Numeric.LinearAlgebra.HMatrix
43import Numeric.GSL.Internal 45import Numeric.GSL.Internal
44import Foreign.Ptr(FunPtr, freeHaskellFunPtr) 46import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
45import Foreign.C.Types 47import Foreign.C.Types
@@ -69,7 +71,7 @@ uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
69 rawpath <- createMIO maxit 4 71 rawpath <- createMIO maxit 4
70 (c_root m fp epsrel (fi maxit) xl xu) 72 (c_root m fp epsrel (fi maxit) xl xu)
71 "root" 73 "root"
72 let it = round (rawpath @@> (maxit-1,0)) 74 let it = round (rawpath `atIndex` (maxit-1,0))
73 path = takeRows it rawpath 75 path = takeRows it rawpath
74 [sol] = toLists $ dropRows (it-1) path 76 [sol] = toLists $ dropRows (it-1) path
75 freeHaskellFunPtr fp 77 freeHaskellFunPtr fp
@@ -100,7 +102,7 @@ uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
100 rawpath <- createMIO maxit 2 102 rawpath <- createMIO maxit 2
101 (c_rootj m fp dfp epsrel (fi maxit) x) 103 (c_rootj m fp dfp epsrel (fi maxit) x)
102 "rootj" 104 "rootj"
103 let it = round (rawpath @@> (maxit-1,0)) 105 let it = round (rawpath `atIndex` (maxit-1,0))
104 path = takeRows it rawpath 106 path = takeRows it rawpath
105 [sol] = toLists $ dropRows (it-1) path 107 [sol] = toLists $ dropRows (it-1) path
106 freeHaskellFunPtr fp 108 freeHaskellFunPtr fp
@@ -132,13 +134,13 @@ root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit ep
132 134
133rootGen m f xi epsabs maxit = unsafePerformIO $ do 135rootGen m f xi epsabs maxit = unsafePerformIO $ do
134 let xiv = fromList xi 136 let xiv = fromList xi
135 n = dim xiv 137 n = size xiv
136 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 138 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
137 rawpath <- vec xiv $ \xiv' -> 139 rawpath <- vec xiv $ \xiv' ->
138 createMIO maxit (2*n+1) 140 createMIO maxit (2*n+1)
139 (c_multiroot m fp epsabs (fi maxit) // xiv') 141 (c_multiroot m fp epsabs (fi maxit) // xiv')
140 "multiroot" 142 "multiroot"
141 let it = round (rawpath @@> (maxit-1,0)) 143 let it = round (rawpath `atIndex` (maxit-1,0))
142 path = takeRows it rawpath 144 path = takeRows it rawpath
143 [sol] = toLists $ dropRows (it-1) path 145 [sol] = toLists $ dropRows (it-1) path
144 freeHaskellFunPtr fp 146 freeHaskellFunPtr fp
@@ -169,14 +171,14 @@ rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun ja
169 171
170rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do 172rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
171 let xiv = fromList xi 173 let xiv = fromList xi
172 n = dim xiv 174 n = size xiv
173 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 175 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
174 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) 176 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
175 rawpath <- vec xiv $ \xiv' -> 177 rawpath <- vec xiv $ \xiv' ->
176 createMIO maxit (2*n+1) 178 createMIO maxit (2*n+1)
177 (c_multirootj m fp jp epsabs (fi maxit) // xiv') 179 (c_multirootj m fp jp epsabs (fi maxit) // xiv')
178 "multiroot" 180 "multiroot"
179 let it = round (rawpath @@> (maxit-1,0)) 181 let it = round (rawpath `atIndex` (maxit-1,0))
180 path = takeRows it rawpath 182 path = takeRows it rawpath
181 [sol] = toLists $ dropRows (it-1) path 183 [sol] = toLists $ dropRows (it-1) path
182 freeHaskellFunPtr fp 184 freeHaskellFunPtr fp
@@ -189,7 +191,7 @@ foreign import ccall safe "multirootj"
189------------------------------------------------------- 191-------------------------------------------------------
190 192
191checkdim1 n v 193checkdim1 n v
192 | dim v == n = v 194 | size v == n = v
193 | otherwise = error $ "Error: "++ show n 195 | otherwise = error $ "Error: "++ show n
194 ++ " components expected in the result of the function supplied to root" 196 ++ " components expected in the result of the function supplied to root"
195 197
diff --git a/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs
new file mode 100644
index 0000000..11b22d3
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs
@@ -0,0 +1,245 @@
1{- |
2Module : Numeric.GSL.Interpolation
3Copyright : (c) Matthew Peddie 2015
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Simulated annealing routines.
9
10<https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html#Simulated-Annealing>
11
12Here is a translation of the simple example given in
13<https://www.gnu.org/software/gsl/manual/html_node/Trivial-example.html#Trivial-example the GSL manual>:
14
15> import Numeric.GSL.SimulatedAnnealing
16> import Numeric.LinearAlgebra.HMatrix
17>
18> main = print $ simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS (Just show)
19>
20> exampleParams = SimulatedAnnealingParams 200 1000 1.0 1.0 0.008 1.003 2.0e-6
21>
22> exampleE x = exp (-(x - 1)**2) * sin (8 * x)
23>
24> exampleM x y = abs $ x - y
25>
26> exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current
27
28The manual states:
29
30> The first example, in one dimensional Cartesian space, sets up an
31> energy function which is a damped sine wave; this has many local
32> minima, but only one global minimum, somewhere between 1.0 and
33> 1.5. The initial guess given is 15.5, which is several local minima
34> away from the global minimum.
35
36This global minimum is around 1.36.
37
38-}
39{-# OPTIONS_GHC -Wall #-}
40
41module Numeric.GSL.SimulatedAnnealing (
42 -- * Searching for minima
43 simanSolve
44 -- * Configuring the annealing process
45 , SimulatedAnnealingParams(..)
46 ) where
47
48import Numeric.GSL.Internal
49import Numeric.LinearAlgebra.HMatrix hiding(step)
50
51import Data.Vector.Storable(generateM)
52import Foreign.Storable(Storable(..))
53import Foreign.Marshal.Utils(with)
54import Foreign.Ptr(Ptr, FunPtr, nullFunPtr)
55import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr)
56import Foreign.C.Types
57import System.IO.Unsafe(unsafePerformIO)
58
59import System.IO (hFlush, stdout)
60
61import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef')
62
63-- | 'SimulatedAnnealingParams' is a translation of the
64-- @gsl_siman_params_t@ structure documented in
65-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-functions.html#Simulated-Annealing-functions the GSL manual>,
66-- which controls the simulated annealing algorithm.
67--
68-- The annealing process is parameterized by the Boltzmann
69-- distribution and the /cooling schedule/. For more details, see
70-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-algorithm.html#Simulated-Annealing-algorithm the relevant section of the manual>.
71data SimulatedAnnealingParams = SimulatedAnnealingParams {
72 n_tries :: CInt -- ^ The number of points to try for each step.
73 , iters_fixed_T :: CInt -- ^ The number of iterations at each temperature
74 , step_size :: Double -- ^ The maximum step size in the random walk
75 , boltzmann_k :: Double -- ^ Boltzmann distribution parameter
76 , cooling_t_initial :: Double -- ^ Initial temperature
77 , cooling_mu_t :: Double -- ^ Cooling rate parameter
78 , cooling_t_min :: Double -- ^ Final temperature
79 } deriving (Eq, Show, Read)
80
81instance Storable SimulatedAnnealingParams where
82 sizeOf p = sizeOf (n_tries p) +
83 sizeOf (iters_fixed_T p) +
84 sizeOf (step_size p) +
85 sizeOf (boltzmann_k p) +
86 sizeOf (cooling_t_initial p) +
87 sizeOf (cooling_mu_t p) +
88 sizeOf (cooling_t_min p)
89 -- TODO(MP): is this safe?
90 alignment p = alignment (step_size p)
91 -- TODO(MP): Is there a more automatic way to write these?
92 peek ptr = SimulatedAnnealingParams <$>
93 peekByteOff ptr 0 <*>
94 peekByteOff ptr i <*>
95 peekByteOff ptr (2*i) <*>
96 peekByteOff ptr (2*i + d) <*>
97 peekByteOff ptr (2*i + 2*d) <*>
98 peekByteOff ptr (2*i + 3*d) <*>
99 peekByteOff ptr (2*i + 4*d)
100 where
101 i = sizeOf (0 :: CInt)
102 d = sizeOf (0 :: Double)
103 poke ptr sap = do
104 pokeByteOff ptr 0 (n_tries sap)
105 pokeByteOff ptr i (iters_fixed_T sap)
106 pokeByteOff ptr (2*i) (step_size sap)
107 pokeByteOff ptr (2*i + d) (boltzmann_k sap)
108 pokeByteOff ptr (2*i + 2*d) (cooling_t_initial sap)
109 pokeByteOff ptr (2*i + 3*d) (cooling_mu_t sap)
110 pokeByteOff ptr (2*i + 4*d) (cooling_t_min sap)
111 where
112 i = sizeOf (0 :: CInt)
113 d = sizeOf (0 :: Double)
114
115-- We use a StablePtr to an IORef so that we can keep hold of
116-- StablePtr values but mutate their contents. A simple 'StablePtr a'
117-- won't work, since we'd have no way to write 'copyConfig'.
118type P a = StablePtr (IORef a)
119
120copyConfig :: P a -> P a -> IO ()
121copyConfig src' dest' = do
122 dest <- deRefStablePtr dest'
123 src <- deRefStablePtr src'
124 readIORef src >>= writeIORef dest
125
126copyConstructConfig :: P a -> IO (P a)
127copyConstructConfig x = do
128 conf <- deRefRead x
129 newconf <- newIORef conf
130 newStablePtr newconf
131
132destroyConfig :: P a -> IO ()
133destroyConfig p = do
134 freeStablePtr p
135
136deRefRead :: P a -> IO a
137deRefRead p = deRefStablePtr p >>= readIORef
138
139wrapEnergy :: (a -> Double) -> P a -> Double
140wrapEnergy f p = unsafePerformIO $ f <$> deRefRead p
141
142wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double
143wrapMetric f x y = unsafePerformIO $ f <$> deRefRead x <*> deRefRead y
144
145wrapStep :: Int
146 -> (Vector Double -> Double -> a -> a)
147 -> GSLRNG
148 -> P a
149 -> Double
150 -> IO ()
151wrapStep nrand f (GSLRNG rng) confptr stepSize = do
152 v <- generateM nrand (\_ -> gslRngUniform rng)
153 conf <- deRefStablePtr confptr
154 modifyIORef' conf $ f v stepSize
155
156wrapPrint :: (a -> String) -> P a -> IO ()
157wrapPrint pf ptr = deRefRead ptr >>= putStr . pf >> hFlush stdout
158
159foreign import ccall safe "wrapper"
160 mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double))
161
162foreign import ccall safe "wrapper"
163 mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))
164
165foreign import ccall safe "wrapper"
166 mkStepFun :: (GSLRNG -> P a -> Double -> IO ())
167 -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))
168
169foreign import ccall safe "wrapper"
170 mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ()))
171
172foreign import ccall safe "wrapper"
173 mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a)))
174
175foreign import ccall safe "wrapper"
176 mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
177
178newtype GSLRNG = GSLRNG (Ptr GSLRNG)
179
180foreign import ccall safe "gsl_rng.h gsl_rng_uniform"
181 gslRngUniform :: Ptr GSLRNG -> IO Double
182
183foreign import ccall safe "gsl-aux.h siman"
184 siman :: CInt -- ^ RNG seed (for repeatability)
185 -> Ptr SimulatedAnnealingParams -- ^ params
186 -> P a -- ^ Configuration
187 -> FunPtr (P a -> Double) -- ^ Energy functional
188 -> FunPtr (P a -> P a -> Double) -- ^ Metric definition
189 -> FunPtr (GSLRNG -> P a -> Double -> IO ()) -- ^ Step evaluation
190 -> FunPtr (P a -> P a -> IO ()) -- ^ Copy config
191 -> FunPtr (P a -> IO (P a)) -- ^ Copy constructor for config
192 -> FunPtr (P a -> IO ()) -- ^ Destructor for config
193 -> FunPtr (P a -> IO ()) -- ^ Print function
194 -> IO CInt
195
196-- |
197-- Calling
198--
199-- > simanSolve seed nrand params x0 e m step print
200--
201-- performs a simulated annealing search through a given space. So
202-- that any configuration type may be used, the space is specified by
203-- providing the functions @e@ (the energy functional) and @m@ (the
204-- metric definition). @x0@ is the initial configuration of the
205-- system. The simulated annealing steps are generated using the
206-- user-provided function @step@, which should randomly construct a
207-- new system configuration.
208--
209-- If 'Nothing' is passed instead of a printing function, no
210-- incremental output will be generated. Otherwise, the GSL-formatted
211-- output, including the configuration description the user function
212-- generates, will be printed to stdout.
213--
214-- Each time the step function is called, it is supplied with a random
215-- vector containing @nrand@ 'Double' values, uniformly distributed in
216-- @[0, 1)@. It should use these values to generate its new
217-- configuration.
218simanSolve :: Int -- ^ Seed for the random number generator
219 -> Int -- ^ @nrand@, the number of random 'Double's the
220 -- step function requires
221 -> SimulatedAnnealingParams -- ^ Parameters to configure the solver
222 -> a -- ^ Initial configuration @x0@
223 -> (a -> Double) -- ^ Energy functional @e@
224 -> (a -> a -> Double) -- ^ Metric definition @m@
225 -> (Vector Double -> Double -> a -> a) -- ^ Stepping function @step@
226 -> Maybe (a -> String) -- ^ Optional printing function
227 -> a -- ^ Best configuration the solver has found
228simanSolve seed nrand params conf e m step printfun =
229 unsafePerformIO $ with params $ \paramptr -> do
230 ewrap <- mkEnergyFun $ wrapEnergy e
231 mwrap <- mkMetricFun $ wrapMetric m
232 stepwrap <- mkStepFun $ wrapStep nrand step
233 confptr <- newIORef conf >>= newStablePtr
234 cpwrap <- mkCopyFun copyConfig
235 ccwrap <- mkCopyConstructorFun copyConstructConfig
236 dwrap <- mkDestructFun destroyConfig
237 pwrap <- case printfun of
238 Nothing -> return nullFunPtr
239 Just pf -> mkDestructFun $ wrapPrint pf
240 siman (fromIntegral seed)
241 paramptr confptr
242 ewrap mwrap stepwrap cpwrap ccwrap dwrap pwrap // check "siman"
243 result <- deRefRead confptr
244 freeStablePtr confptr
245 return result
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
index af79f32..fb982c5 100644
--- a/packages/gsl/src/Numeric/GSL/Vector.hs
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -14,8 +14,7 @@ module Numeric.GSL.Vector (
14 fwriteVector, freadVector, fprintfVector, fscanfVector 14 fwriteVector, freadVector, fprintfVector, fscanfVector
15) where 15) where
16 16
17import Data.Packed 17import Numeric.LinearAlgebra.HMatrix hiding(randomVector, saveMatrix)
18import Numeric.LinearAlgebra(RandDist(..))
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 18import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20 19
21import Foreign.Marshal.Alloc(free) 20import Foreign.Marshal.Alloc(free)
@@ -35,7 +34,7 @@ randomVector :: Int -- ^ seed
35 -> Vector Double 34 -> Vector Double
36randomVector seed dist n = unsafePerformIO $ do 35randomVector seed dist n = unsafePerformIO $ do
37 r <- createVector n 36 r <- createVector n
38 app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" 37 c_random_vector_GSL (fi seed) ((fi.fromEnum) dist) # r #|"randomVectorGSL"
39 return r 38 return r
40 39
41foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV 40foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
@@ -51,7 +50,7 @@ saveMatrix filename fmt m = do
51 charname <- newCString filename 50 charname <- newCString filename
52 charfmt <- newCString fmt 51 charfmt <- newCString fmt
53 let o = if orderOf m == RowMajor then 1 else 0 52 let o = if orderOf m == RowMajor then 1 else 0
54 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" 53 matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
55 free charname 54 free charname
56 free charfmt 55 free charfmt
57 56
@@ -64,7 +63,7 @@ fscanfVector :: FilePath -> Int -> IO (Vector Double)
64fscanfVector filename n = do 63fscanfVector filename n = do
65 charname <- newCString filename 64 charname <- newCString filename
66 res <- createVector n 65 res <- createVector n
67 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" 66 gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
68 free charname 67 free charname
69 return res 68 return res
70 69
@@ -75,7 +74,7 @@ fprintfVector :: FilePath -> String -> Vector Double -> IO ()
75fprintfVector filename fmt v = do 74fprintfVector filename fmt v = do
76 charname <- newCString filename 75 charname <- newCString filename
77 charfmt <- newCString fmt 76 charfmt <- newCString fmt
78 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" 77 gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
79 free charname 78 free charname
80 free charfmt 79 free charfmt
81 80
@@ -86,7 +85,7 @@ freadVector :: FilePath -> Int -> IO (Vector Double)
86freadVector filename n = do 85freadVector filename n = do
87 charname <- newCString filename 86 charname <- newCString filename
88 res <- createVector n 87 res <- createVector n
89 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" 88 gsl_vector_fread charname # res #|"gsl_vector_fread"
90 free charname 89 free charname
91 return res 90 return res
92 91
@@ -96,7 +95,7 @@ foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
96fwriteVector :: FilePath -> Vector Double -> IO () 95fwriteVector :: FilePath -> Vector Double -> IO ()
97fwriteVector filename v = do 96fwriteVector filename v = do
98 charname <- newCString filename 97 charname <- newCString filename
99 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" 98 gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
100 free charname 99 free charname
101 100
102foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV 101foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c
index e1b189c..1ca8199 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-aux.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c
@@ -36,6 +36,8 @@
36#include <gsl/gsl_roots.h> 36#include <gsl/gsl_roots.h>
37#include <gsl/gsl_spline.h> 37#include <gsl/gsl_spline.h>
38#include <gsl/gsl_multifit_nlin.h> 38#include <gsl/gsl_multifit_nlin.h>
39#include <gsl/gsl_siman.h>
40
39#include <string.h> 41#include <string.h>
40#include <stdio.h> 42#include <stdio.h>
41 43
@@ -475,7 +477,30 @@ int uniMinimize(int method, double f(double),
475 OK 477 OK
476} 478}
477 479
478 480int siman(int seed,
481 gsl_siman_params_t *params, void *xp0,
482 double energy(void *), double metric(void *, void *),
483 void step(const gsl_rng *, void *, double),
484 void copy(void *, void *), void *copycons(void *),
485 void destroy(void *), void print(void *)) {
486 DEBUGMSG("siman");
487 gsl_rng *gen = gsl_rng_alloc (gsl_rng_mt19937);
488 gsl_rng_set(gen, seed);
489
490 // The simulated annealing routine doesn't indicate with a return
491 // code how things went -- there's little notion of convergence for
492 // a randomized minimizer on a potentially non-convex problem, and I
493 // suppose it doesn't detect egregious failures like malloc errors
494 // in the copy-constructor.
495 gsl_siman_solve(gen, xp0,
496 energy, step,
497 metric, print,
498 copy, copycons,
499 destroy, 0, *params);
500
501 gsl_rng_free(gen);
502 OK
503}
479 504
480// this version returns info about intermediate steps 505// this version returns info about intermediate steps
481int minimize(int method, double f(int, double*), double tolsize, int maxit, 506int minimize(int method, double f(int, double*), double tolsize, int maxit,
diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c
index 3f2771b..a6bdb55 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-ode.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c
@@ -23,10 +23,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param
23} 23}
24 24
25 25
26int ode(int method, double h, double eps_abs, double eps_rel, 26int ode(int method, int control, double h,
27 double eps_abs, double eps_rel, double a_y, double a_dydt,
27 int f(double, int, const double*, int, double*), 28 int f(double, int, const double*, int, double*),
28 int jac(double, int, const double*, int, int, double*), 29 int jac(double, int, const double*, int, int, double*),
29 KRVEC(xi), KRVEC(ts), RMAT(sol)) { 30 KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) {
30 31
31 const gsl_odeiv_step_type * T; 32 const gsl_odeiv_step_type * T;
32 33
@@ -46,8 +47,16 @@ int ode(int method, double h, double eps_abs, double eps_rel,
46 } 47 }
47 48
48 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); 49 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
49 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); 50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
51 gsl_odeiv_control * c;
52
53 switch(control) {
54 case 0: { c = gsl_odeiv_control_standard_new
55 (eps_abs, eps_rel, a_y, a_dydt); break; }
56 case 1: { c = gsl_odeiv_control_scaled_new
57 (eps_abs, eps_rel, a_y, a_dydt, scp, scn); break; }
58 default: ERROR(BAD_CODE);
59 }
51 60
52 Tode P; 61 Tode P;
53 P.f = f; 62 P.f = f;
@@ -112,10 +121,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param
112} 121}
113 122
114 123
115int ode(int method, double h, double eps_abs, double eps_rel, 124int ode(int method, int control, double h,
125 double eps_abs, double eps_rel, double a_y, double a_dydt,
116 int f(double, int, const double*, int, double*), 126 int f(double, int, const double*, int, double*),
117 int jac(double, int, const double*, int, int, double*), 127 int jac(double, int, const double*, int, int, double*),
118 KRVEC(xi), KRVEC(ts), RMAT(sol)) { 128 KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) {
119 129
120 const gsl_odeiv2_step_type * T; 130 const gsl_odeiv2_step_type * T;
121 131
@@ -141,8 +151,15 @@ int ode(int method, double h, double eps_abs, double eps_rel,
141 151
142 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; 152 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
143 153
144 gsl_odeiv2_driver * d = 154 gsl_odeiv2_driver * d;
145 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); 155
156 switch(control) {
157 case 0: { d = gsl_odeiv2_driver_alloc_standard_new
158 (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt); break; }
159 case 1: { d = gsl_odeiv2_driver_alloc_scaled_new
160 (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt, scp); break; }
161 default: ERROR(BAD_CODE);
162 }
146 163
147 double t = tsp[0]; 164 double t = tsp[0];
148 165
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-14 04:20:22 +0000