diff options
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL')
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/Fitting.hs | 2 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/IO.hs | 42 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs | 38 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/Random.hs | 84 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/Vector.hs | 107 |
5 files changed, 267 insertions, 6 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs index 0a92373..93fb281 100644 --- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs +++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs | |||
@@ -116,7 +116,7 @@ err (model,deriv) dat vsol = zip sol errs where | |||
116 | dof = length dat - (rows cov) | 116 | dof = length dat - (rows cov) |
117 | chi = norm2 (fromList $ cost (resMs model) dat sol) | 117 | chi = norm2 (fromList $ cost (resMs model) dat sol) |
118 | js = fromLists $ jacobian (resDs deriv) dat sol | 118 | js = fromLists $ jacobian (resDs deriv) dat sol |
119 | cov = inv $ trans js <> js | 119 | cov = inv $ trans js <.> js |
120 | errs = toList $ scalar c * sqrt (takeDiag cov) | 120 | errs = toList $ scalar c * sqrt (takeDiag cov) |
121 | 121 | ||
122 | 122 | ||
diff --git a/packages/hmatrix/src/Numeric/GSL/IO.hs b/packages/hmatrix/src/Numeric/GSL/IO.hs new file mode 100644 index 0000000..0d6031a --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/IO.hs | |||
@@ -0,0 +1,42 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.GSL.IO | ||
4 | -- Copyright : (c) Alberto Ruiz 2007-14 | ||
5 | -- License : GPL | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | ----------------------------------------------------------------------------- | ||
10 | |||
11 | module Numeric.GSL.IO ( | ||
12 | saveMatrix, | ||
13 | fwriteVector, freadVector, fprintfVector, fscanfVector, | ||
14 | fileDimensions, loadMatrix, fromFile | ||
15 | ) where | ||
16 | |||
17 | import Data.Packed | ||
18 | import Numeric.GSL.Vector | ||
19 | import System.Process(readProcess) | ||
20 | |||
21 | |||
22 | {- | obtains the number of rows and columns in an ASCII data file | ||
23 | (provisionally using unix's wc). | ||
24 | -} | ||
25 | fileDimensions :: FilePath -> IO (Int,Int) | ||
26 | fileDimensions fname = do | ||
27 | wcres <- readProcess "wc" ["-w",fname] "" | ||
28 | contents <- readFile fname | ||
29 | let tot = read . head . words $ wcres | ||
30 | c = length . head . dropWhile null . map words . lines $ contents | ||
31 | if tot > 0 | ||
32 | then return (tot `div` c, c) | ||
33 | else return (0,0) | ||
34 | |||
35 | -- | Loads a matrix from an ASCII file formatted as a 2D table. | ||
36 | loadMatrix :: FilePath -> IO (Matrix Double) | ||
37 | loadMatrix file = fromFile file =<< fileDimensions file | ||
38 | |||
39 | -- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance). | ||
40 | fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) | ||
41 | fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c) | ||
42 | |||
diff --git a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs index 8bd0cc2..17e2258 100644 --- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs +++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs | |||
@@ -9,13 +9,13 @@ | |||
9 | ----------------------------------------------------------------------------- | 9 | ----------------------------------------------------------------------------- |
10 | 10 | ||
11 | module Numeric.GSL.LinearAlgebra ( | 11 | module Numeric.GSL.LinearAlgebra ( |
12 | randomVector, | 12 | RandDist(..), randomVector, |
13 | saveMatrix, | 13 | saveMatrix, |
14 | fwriteVector, freadVector, fprintfVector, fscanfVector | 14 | fwriteVector, freadVector, fprintfVector, fscanfVector, |
15 | fileDimensions, loadMatrix, fromFile | ||
15 | ) where | 16 | ) where |
16 | 17 | ||
17 | import Data.Packed | 18 | import Data.Packed |
18 | import Numeric.LinearAlgebra.Base(RandDist(..)) | ||
19 | import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) | 19 | import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) |
20 | 20 | ||
21 | import Foreign.Marshal.Alloc(free) | 21 | import Foreign.Marshal.Alloc(free) |
@@ -23,11 +23,16 @@ import Foreign.Ptr(Ptr) | |||
23 | import Foreign.C.Types | 23 | import Foreign.C.Types |
24 | import Foreign.C.String(newCString) | 24 | import Foreign.C.String(newCString) |
25 | import System.IO.Unsafe(unsafePerformIO) | 25 | import System.IO.Unsafe(unsafePerformIO) |
26 | import System.Process(readProcess) | ||
26 | 27 | ||
27 | fromei x = fromIntegral (fromEnum x) :: CInt | 28 | fromei x = fromIntegral (fromEnum x) :: CInt |
28 | 29 | ||
29 | ----------------------------------------------------------------------- | 30 | ----------------------------------------------------------------------- |
30 | 31 | ||
32 | data RandDist = Uniform -- ^ uniform distribution in [0,1) | ||
33 | | Gaussian -- ^ normal distribution with mean zero and standard deviation one | ||
34 | deriving Enum | ||
35 | |||
31 | -- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. | 36 | -- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. |
32 | randomVector :: Int -- ^ seed | 37 | randomVector :: Int -- ^ seed |
33 | -> RandDist -- ^ distribution | 38 | -> RandDist -- ^ distribution |
@@ -35,10 +40,10 @@ randomVector :: Int -- ^ seed | |||
35 | -> Vector Double | 40 | -> Vector Double |
36 | randomVector seed dist n = unsafePerformIO $ do | 41 | randomVector seed dist n = unsafePerformIO $ do |
37 | r <- createVector n | 42 | r <- createVector n |
38 | app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" | 43 | app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" |
39 | return r | 44 | return r |
40 | 45 | ||
41 | foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV | 46 | foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV |
42 | 47 | ||
43 | -------------------------------------------------------------------------------- | 48 | -------------------------------------------------------------------------------- |
44 | 49 | ||
@@ -105,3 +110,26 @@ type PD = Ptr Double -- | |||
105 | type TV = CInt -> PD -> IO CInt -- | 110 | type TV = CInt -> PD -> IO CInt -- |
106 | type TM = CInt -> CInt -> PD -> IO CInt -- | 111 | type TM = CInt -> CInt -> PD -> IO CInt -- |
107 | 112 | ||
113 | -------------------------------------------------------------------------------- | ||
114 | |||
115 | {- | obtains the number of rows and columns in an ASCII data file | ||
116 | (provisionally using unix's wc). | ||
117 | -} | ||
118 | fileDimensions :: FilePath -> IO (Int,Int) | ||
119 | fileDimensions fname = do | ||
120 | wcres <- readProcess "wc" ["-w",fname] "" | ||
121 | contents <- readFile fname | ||
122 | let tot = read . head . words $ wcres | ||
123 | c = length . head . dropWhile null . map words . lines $ contents | ||
124 | if tot > 0 | ||
125 | then return (tot `div` c, c) | ||
126 | else return (0,0) | ||
127 | |||
128 | -- | Loads a matrix from an ASCII file formatted as a 2D table. | ||
129 | loadMatrix :: FilePath -> IO (Matrix Double) | ||
130 | loadMatrix file = fromFile file =<< fileDimensions file | ||
131 | |||
132 | -- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance). | ||
133 | fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) | ||
134 | fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c) | ||
135 | |||
diff --git a/packages/hmatrix/src/Numeric/GSL/Random.hs b/packages/hmatrix/src/Numeric/GSL/Random.hs new file mode 100644 index 0000000..2872b17 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Random.hs | |||
@@ -0,0 +1,84 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.GSL.Random | ||
4 | -- Copyright : (c) Alberto Ruiz 2009-14 | ||
5 | -- License : GPL | ||
6 | -- | ||
7 | -- Maintainer : Alberto Ruiz | ||
8 | -- Stability : provisional | ||
9 | -- | ||
10 | -- Random vectors and matrices. | ||
11 | -- | ||
12 | ----------------------------------------------------------------------------- | ||
13 | |||
14 | module Numeric.GSL.Random ( | ||
15 | Seed, | ||
16 | RandDist(..), | ||
17 | randomVector, | ||
18 | gaussianSample, | ||
19 | uniformSample, | ||
20 | rand, randn | ||
21 | ) where | ||
22 | |||
23 | import Numeric.GSL.Vector | ||
24 | import Numeric.Container | ||
25 | import Numeric.LinearAlgebra(Seed,RandDist(..),cholSH) | ||
26 | import System.Random(randomIO) | ||
27 | |||
28 | |||
29 | |||
30 | |||
31 | -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate | ||
32 | -- Gaussian distribution. | ||
33 | gaussianSample :: Seed | ||
34 | -> Int -- ^ number of rows | ||
35 | -> Vector Double -- ^ mean vector | ||
36 | -> Matrix Double -- ^ covariance matrix | ||
37 | -> Matrix Double -- ^ result | ||
38 | gaussianSample seed n med cov = m where | ||
39 | c = dim med | ||
40 | meds = konst' 1 n `outer` med | ||
41 | rs = reshape c $ randomVector seed Gaussian (c * n) | ||
42 | m = rs `mXm` cholSH cov `add` meds | ||
43 | |||
44 | -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate | ||
45 | -- uniform distribution. | ||
46 | uniformSample :: Seed | ||
47 | -> Int -- ^ number of rows | ||
48 | -> [(Double,Double)] -- ^ ranges for each column | ||
49 | -> Matrix Double -- ^ result | ||
50 | uniformSample seed n rgs = m where | ||
51 | (as,bs) = unzip rgs | ||
52 | a = fromList as | ||
53 | cs = zipWith subtract as bs | ||
54 | d = dim a | ||
55 | dat = toRows $ reshape n $ randomVector seed Uniform (n*d) | ||
56 | am = konst' 1 n `outer` a | ||
57 | m = fromColumns (zipWith scale cs dat) `add` am | ||
58 | |||
59 | -- | pseudorandom matrix with uniform elements between 0 and 1 | ||
60 | randm :: RandDist | ||
61 | -> Int -- ^ rows | ||
62 | -> Int -- ^ columns | ||
63 | -> IO (Matrix Double) | ||
64 | randm d r c = do | ||
65 | seed <- randomIO | ||
66 | return (reshape c $ randomVector seed d (r*c)) | ||
67 | |||
68 | -- | pseudorandom matrix with uniform elements between 0 and 1 | ||
69 | rand :: Int -> Int -> IO (Matrix Double) | ||
70 | rand = randm Uniform | ||
71 | |||
72 | {- | pseudorandom matrix with normal elements | ||
73 | |||
74 | >>> x <- randn 3 5 | ||
75 | >>> disp 3 x | ||
76 | 3x5 | ||
77 | 0.386 -1.141 0.491 -0.510 1.512 | ||
78 | 0.069 -0.919 1.022 -0.181 0.745 | ||
79 | 0.313 -0.670 -0.097 -1.575 -0.583 | ||
80 | |||
81 | -} | ||
82 | randn :: Int -> Int -> IO (Matrix Double) | ||
83 | randn = randm Gaussian | ||
84 | |||
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs new file mode 100644 index 0000000..af79f32 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs | |||
@@ -0,0 +1,107 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.GSL.Vector | ||
4 | -- Copyright : (c) Alberto Ruiz 2007-14 | ||
5 | -- License : GPL | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | ----------------------------------------------------------------------------- | ||
10 | |||
11 | module Numeric.GSL.Vector ( | ||
12 | randomVector, | ||
13 | saveMatrix, | ||
14 | fwriteVector, freadVector, fprintfVector, fscanfVector | ||
15 | ) where | ||
16 | |||
17 | import Data.Packed | ||
18 | import Numeric.LinearAlgebra(RandDist(..)) | ||
19 | import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) | ||
20 | |||
21 | import Foreign.Marshal.Alloc(free) | ||
22 | import Foreign.Ptr(Ptr) | ||
23 | import Foreign.C.Types | ||
24 | import Foreign.C.String(newCString) | ||
25 | import System.IO.Unsafe(unsafePerformIO) | ||
26 | |||
27 | fromei x = fromIntegral (fromEnum x) :: CInt | ||
28 | |||
29 | ----------------------------------------------------------------------- | ||
30 | |||
31 | -- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. | ||
32 | randomVector :: Int -- ^ seed | ||
33 | -> RandDist -- ^ distribution | ||
34 | -> Int -- ^ vector size | ||
35 | -> Vector Double | ||
36 | randomVector seed dist n = unsafePerformIO $ do | ||
37 | r <- createVector n | ||
38 | app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" | ||
39 | return r | ||
40 | |||
41 | foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV | ||
42 | |||
43 | -------------------------------------------------------------------------------- | ||
44 | |||
45 | -- | Saves a matrix as 2D ASCII table. | ||
46 | saveMatrix :: FilePath | ||
47 | -> String -- ^ format (%f, %g, %e) | ||
48 | -> Matrix Double | ||
49 | -> IO () | ||
50 | saveMatrix filename fmt m = do | ||
51 | charname <- newCString filename | ||
52 | charfmt <- newCString fmt | ||
53 | let o = if orderOf m == RowMajor then 1 else 0 | ||
54 | app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" | ||
55 | free charname | ||
56 | free charfmt | ||
57 | |||
58 | foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM | ||
59 | |||
60 | -------------------------------------------------------------------------------- | ||
61 | |||
62 | -- | Loads a vector from an ASCII file (the number of elements must be known in advance). | ||
63 | fscanfVector :: FilePath -> Int -> IO (Vector Double) | ||
64 | fscanfVector filename n = do | ||
65 | charname <- newCString filename | ||
66 | res <- createVector n | ||
67 | app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" | ||
68 | free charname | ||
69 | return res | ||
70 | |||
71 | foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV | ||
72 | |||
73 | -- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file. | ||
74 | fprintfVector :: FilePath -> String -> Vector Double -> IO () | ||
75 | fprintfVector filename fmt v = do | ||
76 | charname <- newCString filename | ||
77 | charfmt <- newCString fmt | ||
78 | app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" | ||
79 | free charname | ||
80 | free charfmt | ||
81 | |||
82 | foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV | ||
83 | |||
84 | -- | Loads a vector from a binary file (the number of elements must be known in advance). | ||
85 | freadVector :: FilePath -> Int -> IO (Vector Double) | ||
86 | freadVector filename n = do | ||
87 | charname <- newCString filename | ||
88 | res <- createVector n | ||
89 | app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" | ||
90 | free charname | ||
91 | return res | ||
92 | |||
93 | foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV | ||
94 | |||
95 | -- | Saves the elements of a vector to a binary file. | ||
96 | fwriteVector :: FilePath -> Vector Double -> IO () | ||
97 | fwriteVector filename v = do | ||
98 | charname <- newCString filename | ||
99 | app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" | ||
100 | free charname | ||
101 | |||
102 | foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV | ||
103 | |||
104 | type PD = Ptr Double -- | ||
105 | type TV = CInt -> PD -> IO CInt -- | ||
106 | type TM = CInt -> CInt -> PD -> IO CInt -- | ||
107 | |||