summaryrefslogtreecommitdiff
path: root/packages/base
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base')
-rw-r--r--packages/base/hmatrix-base.cabal2
-rw-r--r--packages/base/src/C/vector-aux.c49
-rw-r--r--packages/base/src/Data/Packed/IO.hs2
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs8
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Base.hs20
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Random.hs80
-rw-r--r--packages/base/src/Numeric/Vectorized.hs30
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
701double 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
726int 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.
152loadMatrix :: FilePath -> IO (Matrix Double) 154loadMatrix :: FilePath -> IO (Matrix Double)
153loadMatrix f = do 155loadMatrix 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
84class Contraction a b c | a b -> c 84class 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
89Examples: 89Examples:
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
228x25c7, white diamond
229
230-}
231(◇) :: Contraction a b c => a -> b -> c 227(◇) :: Contraction a b c => a -> b -> c
232infixl 7 ◇ 228infixl 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
128import Numeric.LinearAlgebra.Data 131import Numeric.LinearAlgebra.Data
@@ -132,6 +135,7 @@ import Numeric.Vector()
132import Data.Packed.Numeric 135import Data.Packed.Numeric
133import Numeric.LinearAlgebra.Algorithms 136import Numeric.LinearAlgebra.Algorithms
134import Numeric.LinearAlgebra.Util 137import Numeric.LinearAlgebra.Util
138import 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
13module Numeric.LinearAlgebra.Random (
14 Seed,
15 RandDist(..),
16 randomVector,
17 gaussianSample,
18 uniformSample,
19 rand, randn
20) where
21
22import Numeric.Vectorized
23import Data.Packed.Numeric
24import Numeric.LinearAlgebra.Algorithms
25import System.Random(randomIO)
26
27
28-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
29-- Gaussian distribution.
30gaussianSample :: Seed
31 -> Int -- ^ number of rows
32 -> Vector Double -- ^ mean vector
33 -> Matrix Double -- ^ covariance matrix
34 -> Matrix Double -- ^ result
35gaussianSample 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.
43uniformSample :: Seed
44 -> Int -- ^ number of rows
45 -> [(Double,Double)] -- ^ ranges for each column
46 -> Matrix Double -- ^ result
47uniformSample 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
57randm :: RandDist
58 -> Int -- ^ rows
59 -> Int -- ^ columns
60 -> IO (Matrix Double)
61randm 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
66rand :: Int -> Int -> IO (Matrix Double)
67rand = randm Uniform
68
69{- | pseudorandom matrix with normal elements
70
71>>> disp 3 =<< randn 3 5
723x5
730.386 -1.141 0.491 -0.510 1.512
740.069 -0.919 1.022 -0.181 0.745
750.313 -0.670 -0.097 -1.575 -0.583
76
77-}
78randn :: Int -> Int -> IO (Matrix Double)
79randn = 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
23import Data.Packed.Internal.Common 24import Data.Packed.Internal.Common
@@ -308,7 +309,13 @@ vectorScan s = do
308foreign import ccall unsafe "saveMatrix" c_saveMatrix 309foreign import ccall unsafe "saveMatrix" c_saveMatrix
309 :: CString -> CString -> TM 310 :: CString -> CString -> TM
310 311
311saveMatrix :: FilePath -> String -> Matrix Double -> IO () 312{- | save a matrix as a 2D ASCII table
313-}
314saveMatrix
315 :: FilePath
316 -> String -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.)
317 -> Matrix Double
318 -> IO ()
312saveMatrix name format m = do 319saveMatrix 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
329type Seed = Int
330
331data 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).
336randomVector :: Seed
337 -> RandDist -- ^ distribution
338 -> Int -- ^ vector size
339 -> Vector Double
340randomVector 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
345foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
320 346