diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-20 19:32:19 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-20 19:32:19 +0200 |
commit | 9412ef25f2a9d6f6ca233ef123a01c3f4145ffa4 (patch) | |
tree | b5c7ecd0534d0dc34242e53a3e66b87fc2ec0b0f /packages/base | |
parent | d0fc6c7192badfa6f03baf0e02e0cf2a73c3906b (diff) |
random numbers in base
Diffstat (limited to 'packages/base')
-rw-r--r-- | packages/base/hmatrix-base.cabal | 2 | ||||
-rw-r--r-- | packages/base/src/C/vector-aux.c | 49 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/IO.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Numeric.hs | 8 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Base.hs | 20 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Random.hs | 80 | ||||
-rw-r--r-- | packages/base/src/Numeric/Vectorized.hs | 30 |
7 files changed, 175 insertions, 16 deletions
diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 30619e1..b884518 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal | |||
@@ -23,6 +23,7 @@ library | |||
23 | binary, | 23 | binary, |
24 | array, | 24 | array, |
25 | deepseq, | 25 | deepseq, |
26 | random, | ||
26 | storable-complex, | 27 | storable-complex, |
27 | vector >= 0.8 | 28 | vector >= 0.8 |
28 | 29 | ||
@@ -56,6 +57,7 @@ library | |||
56 | Numeric.Matrix | 57 | Numeric.Matrix |
57 | Data.Packed.Internal.Numeric | 58 | Data.Packed.Internal.Numeric |
58 | Numeric.LinearAlgebra.Util.Convolution | 59 | Numeric.LinearAlgebra.Util.Convolution |
60 | Numeric.LinearAlgebra.Random | ||
59 | Numeric.Conversion | 61 | Numeric.Conversion |
60 | 62 | ||
61 | C-sources: src/C/lapack-aux.c | 63 | C-sources: src/C/lapack-aux.c |
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c index f1bb371..5b9c171 100644 --- a/packages/base/src/C/vector-aux.c +++ b/packages/base/src/C/vector-aux.c | |||
@@ -695,3 +695,52 @@ int saveMatrix(char * file, char * format, KDMAT(a)){ | |||
695 | OK | 695 | OK |
696 | } | 696 | } |
697 | 697 | ||
698 | //////////////////////////////////////////////////////////////////////////////// | ||
699 | |||
700 | // http://c-faq.com/lib/gaussian.html | ||
701 | double gaussrand() | ||
702 | { | ||
703 | static double V1, V2, S; | ||
704 | static int phase = 0; | ||
705 | double X; | ||
706 | |||
707 | if(phase == 0) { | ||
708 | do { | ||
709 | double U1 = (double)rand() / RAND_MAX; | ||
710 | double U2 = (double)rand() / RAND_MAX; | ||
711 | |||
712 | V1 = 2 * U1 - 1; | ||
713 | V2 = 2 * U2 - 1; | ||
714 | S = V1 * V1 + V2 * V2; | ||
715 | } while(S >= 1 || S == 0); | ||
716 | |||
717 | X = V1 * sqrt(-2 * log(S) / S); | ||
718 | } else | ||
719 | X = V2 * sqrt(-2 * log(S) / S); | ||
720 | |||
721 | phase = 1 - phase; | ||
722 | |||
723 | return X; | ||
724 | } | ||
725 | |||
726 | int random_vector(int seed, int code, DVEC(r)) { | ||
727 | srand(seed); | ||
728 | int k; | ||
729 | switch (code) { | ||
730 | case 0: { // uniform | ||
731 | for (k=0; k<rn; k++) { | ||
732 | rp[k] = (double)rand()/RAND_MAX; | ||
733 | } | ||
734 | OK | ||
735 | } | ||
736 | case 1: { // gaussian | ||
737 | for (k=0; k<rn; k++) { | ||
738 | rp[k] = gaussrand(); | ||
739 | } | ||
740 | OK | ||
741 | } | ||
742 | |||
743 | default: ERROR(BAD_CODE); | ||
744 | } | ||
745 | } | ||
746 | |||
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs index db03d5f..94fb30a 100644 --- a/packages/base/src/Data/Packed/IO.hs +++ b/packages/base/src/Data/Packed/IO.hs | |||
@@ -149,6 +149,8 @@ apparentCols s = f . dropWhile null . map words . lines <$> readFile s | |||
149 | f [] = 0 | 149 | f [] = 0 |
150 | f (x:_) = length x | 150 | f (x:_) = length x |
151 | 151 | ||
152 | |||
153 | -- | load a matrix from an ASCII file formatted as a 2D table. | ||
152 | loadMatrix :: FilePath -> IO (Matrix Double) | 154 | loadMatrix :: FilePath -> IO (Matrix Double) |
153 | loadMatrix f = do | 155 | loadMatrix f = do |
154 | v <- vectorScan f | 156 | v <- vectorScan f |
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs index 6036e8c..d130ecd 100644 --- a/packages/base/src/Data/Packed/Numeric.hs +++ b/packages/base/src/Data/Packed/Numeric.hs | |||
@@ -84,7 +84,7 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n | |||
84 | class Contraction a b c | a b -> c | 84 | class Contraction a b c | a b -> c |
85 | where | 85 | where |
86 | infixl 7 <.> | 86 | infixl 7 <.> |
87 | {- | Matrix product, matrix vector product, and dot product | 87 | {- | Matrix product, matrix - vector product, and dot product |
88 | 88 | ||
89 | Examples: | 89 | Examples: |
90 | 90 | ||
@@ -223,11 +223,7 @@ instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e | |||
223 | 223 | ||
224 | -------------------------------------------------------------------------------- | 224 | -------------------------------------------------------------------------------- |
225 | 225 | ||
226 | {- | alternative operator for '(\<.\>)' | 226 | -- | alternative unicode symbol (25c7) for the contraction operator '(\<.\>)' |
227 | |||
228 | x25c7, white diamond | ||
229 | |||
230 | -} | ||
231 | (◇) :: Contraction a b c => a -> b -> c | 227 | (◇) :: Contraction a b c => a -> b -> c |
232 | infixl 7 ◇ | 228 | infixl 7 ◇ |
233 | (◇) = (<.>) | 229 | (◇) = (<.>) |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Base.hs b/packages/base/src/Numeric/LinearAlgebra/Base.hs index 395c84a..8d44d26 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Base.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Base.hs | |||
@@ -35,10 +35,13 @@ module Numeric.LinearAlgebra.Base ( | |||
35 | -- , 5.0, 5.0, 7.0 ] | 35 | -- , 5.0, 5.0, 7.0 ] |
36 | -- | 36 | -- |
37 | 37 | ||
38 | -- * Products | 38 | -- * Matrix product |
39 | (<.>), | 39 | (<.>), |
40 | 40 | ||
41 | -- | The matrix product is also implemented in the "Data.Monoid" instance for Matrix, where | 41 | -- | This operator can also be written using the unicode symbol ◇ (25c7). |
42 | -- | ||
43 | |||
44 | -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where | ||
42 | -- single-element matrices (created from numeric literals or using 'scalar') | 45 | -- single-element matrices (created from numeric literals or using 'scalar') |
43 | -- are used for scaling. | 46 | -- are used for scaling. |
44 | -- | 47 | -- |
@@ -48,9 +51,9 @@ module Numeric.LinearAlgebra.Base ( | |||
48 | -- [ 1.0, 4.0, 0.0 | 51 | -- [ 1.0, 4.0, 0.0 |
49 | -- , 4.0, 10.0, 0.0 ] | 52 | -- , 4.0, 10.0, 0.0 ] |
50 | -- | 53 | -- |
51 | -- mconcat uses 'optimiseMult' to get the optimal association order. | 54 | -- 'mconcat' uses 'optimiseMult' to get the optimal association order. |
52 | 55 | ||
53 | (◇), | 56 | -- * Other products |
54 | outer, kronecker, cross, | 57 | outer, kronecker, cross, |
55 | scale, | 58 | scale, |
56 | sumElements, prodElements, absSum, | 59 | sumElements, prodElements, absSum, |
@@ -114,15 +117,15 @@ module Numeric.LinearAlgebra.Base ( | |||
114 | -- * Norms | 117 | -- * Norms |
115 | norm1, norm2, normInf, pnorm, NormType(..), | 118 | norm1, norm2, normInf, pnorm, NormType(..), |
116 | 119 | ||
117 | -- * Correlation and Convolution | 120 | -- * Correlation and convolution |
118 | corr, conv, corrMin, corr2, conv2, | 121 | corr, conv, corrMin, corr2, conv2, |
119 | 122 | ||
120 | -- * Random arrays | 123 | -- * Random arrays |
121 | 124 | ||
122 | -- | rand, randn, RandDist(..), randomVector, gaussianSample, uniformSample | 125 | RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, |
123 | 126 | ||
124 | -- * Misc | 127 | -- * Misc |
125 | meanCov, peps, relativeError, haussholder, optimiseMult, udot | 128 | meanCov, peps, relativeError, haussholder, optimiseMult, udot, Seed, (◇) |
126 | ) where | 129 | ) where |
127 | 130 | ||
128 | import Numeric.LinearAlgebra.Data | 131 | import Numeric.LinearAlgebra.Data |
@@ -132,6 +135,7 @@ import Numeric.Vector() | |||
132 | import Data.Packed.Numeric | 135 | import Data.Packed.Numeric |
133 | import Numeric.LinearAlgebra.Algorithms | 136 | import Numeric.LinearAlgebra.Algorithms |
134 | import Numeric.LinearAlgebra.Util | 137 | import Numeric.LinearAlgebra.Util |
138 | import Numeric.LinearAlgebra.Random | ||
135 | 139 | ||
136 | 140 | ||
137 | 141 | ||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Random.hs b/packages/base/src/Numeric/LinearAlgebra/Random.hs new file mode 100644 index 0000000..b36c7a3 --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/Random.hs | |||
@@ -0,0 +1,80 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.LinearAlgebra.Random | ||
4 | -- Copyright : (c) Alberto Ruiz 2009-14 | ||
5 | -- License : BSD3 | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | -- Random vectors and matrices. | ||
10 | -- | ||
11 | ----------------------------------------------------------------------------- | ||
12 | |||
13 | module Numeric.LinearAlgebra.Random ( | ||
14 | Seed, | ||
15 | RandDist(..), | ||
16 | randomVector, | ||
17 | gaussianSample, | ||
18 | uniformSample, | ||
19 | rand, randn | ||
20 | ) where | ||
21 | |||
22 | import Numeric.Vectorized | ||
23 | import Data.Packed.Numeric | ||
24 | import Numeric.LinearAlgebra.Algorithms | ||
25 | import System.Random(randomIO) | ||
26 | |||
27 | |||
28 | -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate | ||
29 | -- Gaussian distribution. | ||
30 | gaussianSample :: Seed | ||
31 | -> Int -- ^ number of rows | ||
32 | -> Vector Double -- ^ mean vector | ||
33 | -> Matrix Double -- ^ covariance matrix | ||
34 | -> Matrix Double -- ^ result | ||
35 | gaussianSample seed n med cov = m where | ||
36 | c = dim med | ||
37 | meds = konst' 1 n `outer` med | ||
38 | rs = reshape c $ randomVector seed Gaussian (c * n) | ||
39 | m = rs `mXm` cholSH cov `add` meds | ||
40 | |||
41 | -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate | ||
42 | -- uniform distribution. | ||
43 | uniformSample :: Seed | ||
44 | -> Int -- ^ number of rows | ||
45 | -> [(Double,Double)] -- ^ ranges for each column | ||
46 | -> Matrix Double -- ^ result | ||
47 | uniformSample seed n rgs = m where | ||
48 | (as,bs) = unzip rgs | ||
49 | a = fromList as | ||
50 | cs = zipWith subtract as bs | ||
51 | d = dim a | ||
52 | dat = toRows $ reshape n $ randomVector seed Uniform (n*d) | ||
53 | am = konst' 1 n `outer` a | ||
54 | m = fromColumns (zipWith scale cs dat) `add` am | ||
55 | |||
56 | -- | pseudorandom matrix with uniform elements between 0 and 1 | ||
57 | randm :: RandDist | ||
58 | -> Int -- ^ rows | ||
59 | -> Int -- ^ columns | ||
60 | -> IO (Matrix Double) | ||
61 | randm d r c = do | ||
62 | seed <- randomIO | ||
63 | return (reshape c $ randomVector seed d (r*c)) | ||
64 | |||
65 | -- | pseudorandom matrix with uniform elements between 0 and 1 | ||
66 | rand :: Int -> Int -> IO (Matrix Double) | ||
67 | rand = randm Uniform | ||
68 | |||
69 | {- | pseudorandom matrix with normal elements | ||
70 | |||
71 | >>> disp 3 =<< randn 3 5 | ||
72 | 3x5 | ||
73 | 0.386 -1.141 0.491 -0.510 1.512 | ||
74 | 0.069 -0.919 1.022 -0.181 0.745 | ||
75 | 0.313 -0.670 -0.097 -1.575 -0.583 | ||
76 | |||
77 | -} | ||
78 | randn :: Int -> Int -> IO (Matrix Double) | ||
79 | randn = randm Gaussian | ||
80 | |||
diff --git a/packages/base/src/Numeric/Vectorized.hs b/packages/base/src/Numeric/Vectorized.hs index a2d7f70..5aebb14 100644 --- a/packages/base/src/Numeric/Vectorized.hs +++ b/packages/base/src/Numeric/Vectorized.hs | |||
@@ -17,7 +17,8 @@ module Numeric.Vectorized ( | |||
17 | FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, | 17 | FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, |
18 | FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, | 18 | FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, |
19 | FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, | 19 | FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, |
20 | vectorScan, saveMatrix | 20 | vectorScan, saveMatrix, |
21 | Seed, RandDist(..), randomVector | ||
21 | ) where | 22 | ) where |
22 | 23 | ||
23 | import Data.Packed.Internal.Common | 24 | import Data.Packed.Internal.Common |
@@ -308,7 +309,13 @@ vectorScan s = do | |||
308 | foreign import ccall unsafe "saveMatrix" c_saveMatrix | 309 | foreign import ccall unsafe "saveMatrix" c_saveMatrix |
309 | :: CString -> CString -> TM | 310 | :: CString -> CString -> TM |
310 | 311 | ||
311 | saveMatrix :: FilePath -> String -> Matrix Double -> IO () | 312 | {- | save a matrix as a 2D ASCII table |
313 | -} | ||
314 | saveMatrix | ||
315 | :: FilePath | ||
316 | -> String -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.) | ||
317 | -> Matrix Double | ||
318 | -> IO () | ||
312 | saveMatrix name format m = do | 319 | saveMatrix name format m = do |
313 | cname <- newCString name | 320 | cname <- newCString name |
314 | cformat <- newCString format | 321 | cformat <- newCString format |
@@ -317,4 +324,23 @@ saveMatrix name format m = do | |||
317 | free cformat | 324 | free cformat |
318 | return () | 325 | return () |
319 | 326 | ||
327 | -------------------------------------------------------------------------------- | ||
328 | |||
329 | type Seed = Int | ||
330 | |||
331 | data RandDist = Uniform -- ^ uniform distribution in [0,1) | ||
332 | | Gaussian -- ^ normal distribution with mean zero and standard deviation one | ||
333 | deriving Enum | ||
334 | |||
335 | -- | Obtains a vector of pseudorandom elements (use randomIO to get a random seed). | ||
336 | randomVector :: Seed | ||
337 | -> RandDist -- ^ distribution | ||
338 | -> Int -- ^ vector size | ||
339 | -> Vector Double | ||
340 | randomVector seed dist n = unsafePerformIO $ do | ||
341 | r <- createVector n | ||
342 | app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" | ||
343 | return r | ||
344 | |||
345 | foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV | ||
320 | 346 | ||