From 9412ef25f2a9d6f6ca233ef123a01c3f4145ffa4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 20 May 2014 19:32:19 +0200 Subject: random numbers in base --- packages/base/hmatrix-base.cabal | 2 + packages/base/src/C/vector-aux.c | 49 ++++++++++++++ packages/base/src/Data/Packed/IO.hs | 2 + packages/base/src/Data/Packed/Numeric.hs | 8 +-- packages/base/src/Numeric/LinearAlgebra/Base.hs | 20 +++--- packages/base/src/Numeric/LinearAlgebra/Random.hs | 80 +++++++++++++++++++++++ packages/base/src/Numeric/Vectorized.hs | 30 ++++++++- 7 files changed, 175 insertions(+), 16 deletions(-) create mode 100644 packages/base/src/Numeric/LinearAlgebra/Random.hs (limited to 'packages/base') 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 binary, array, deepseq, + random, storable-complex, vector >= 0.8 @@ -56,6 +57,7 @@ library Numeric.Matrix Data.Packed.Internal.Numeric Numeric.LinearAlgebra.Util.Convolution + Numeric.LinearAlgebra.Random Numeric.Conversion 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)){ OK } +//////////////////////////////////////////////////////////////////////////////// + +// http://c-faq.com/lib/gaussian.html +double gaussrand() +{ + static double V1, V2, S; + static int phase = 0; + double X; + + if(phase == 0) { + do { + double U1 = (double)rand() / RAND_MAX; + double U2 = (double)rand() / RAND_MAX; + + V1 = 2 * U1 - 1; + V2 = 2 * U2 - 1; + S = V1 * V1 + V2 * V2; + } while(S >= 1 || S == 0); + + X = V1 * sqrt(-2 * log(S) / S); + } else + X = V2 * sqrt(-2 * log(S) / S); + + phase = 1 - phase; + + return X; +} + +int random_vector(int seed, int code, DVEC(r)) { + srand(seed); + int k; + switch (code) { + case 0: { // uniform + for (k=0; k readFile s f [] = 0 f (x:_) = length x + +-- | load a matrix from an ASCII file formatted as a 2D table. loadMatrix :: FilePath -> IO (Matrix Double) loadMatrix f = do 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 class Contraction a b c | a b -> c where infixl 7 <.> - {- | Matrix product, matrix vector product, and dot product + {- | Matrix product, matrix - vector product, and dot product Examples: @@ -223,11 +223,7 @@ instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e -------------------------------------------------------------------------------- -{- | alternative operator for '(\<.\>)' - -x25c7, white diamond - --} +-- | alternative unicode symbol (25c7) for the contraction operator '(\<.\>)' (◇) :: Contraction a b c => a -> b -> c infixl 7 ◇ (◇) = (<.>) 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 ( -- , 5.0, 5.0, 7.0 ] -- - -- * Products + -- * Matrix product (<.>), - -- | The matrix product is also implemented in the "Data.Monoid" instance for Matrix, where + -- | This operator can also be written using the unicode symbol ◇ (25c7). + -- + + -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where -- single-element matrices (created from numeric literals or using 'scalar') -- are used for scaling. -- @@ -48,9 +51,9 @@ module Numeric.LinearAlgebra.Base ( -- [ 1.0, 4.0, 0.0 -- , 4.0, 10.0, 0.0 ] -- - -- mconcat uses 'optimiseMult' to get the optimal association order. - - (◇), + -- 'mconcat' uses 'optimiseMult' to get the optimal association order. + + -- * Other products outer, kronecker, cross, scale, sumElements, prodElements, absSum, @@ -114,15 +117,15 @@ module Numeric.LinearAlgebra.Base ( -- * Norms norm1, norm2, normInf, pnorm, NormType(..), - -- * Correlation and Convolution + -- * Correlation and convolution corr, conv, corrMin, corr2, conv2, -- * Random arrays - -- | rand, randn, RandDist(..), randomVector, gaussianSample, uniformSample + RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, -- * Misc - meanCov, peps, relativeError, haussholder, optimiseMult, udot + meanCov, peps, relativeError, haussholder, optimiseMult, udot, Seed, (◇) ) where import Numeric.LinearAlgebra.Data @@ -132,6 +135,7 @@ import Numeric.Vector() import Data.Packed.Numeric import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Util +import Numeric.LinearAlgebra.Random 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 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.LinearAlgebra.Random +-- Copyright : (c) Alberto Ruiz 2009-14 +-- License : BSD3 +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Random vectors and matrices. +-- +----------------------------------------------------------------------------- + +module Numeric.LinearAlgebra.Random ( + Seed, + RandDist(..), + randomVector, + gaussianSample, + uniformSample, + rand, randn +) where + +import Numeric.Vectorized +import Data.Packed.Numeric +import Numeric.LinearAlgebra.Algorithms +import System.Random(randomIO) + + +-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate +-- Gaussian distribution. +gaussianSample :: Seed + -> Int -- ^ number of rows + -> Vector Double -- ^ mean vector + -> Matrix Double -- ^ covariance matrix + -> Matrix Double -- ^ result +gaussianSample seed n med cov = m where + c = dim med + meds = konst' 1 n `outer` med + rs = reshape c $ randomVector seed Gaussian (c * n) + m = rs `mXm` cholSH cov `add` meds + +-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate +-- uniform distribution. +uniformSample :: Seed + -> Int -- ^ number of rows + -> [(Double,Double)] -- ^ ranges for each column + -> Matrix Double -- ^ result +uniformSample seed n rgs = m where + (as,bs) = unzip rgs + a = fromList as + cs = zipWith subtract as bs + d = dim a + dat = toRows $ reshape n $ randomVector seed Uniform (n*d) + am = konst' 1 n `outer` a + m = fromColumns (zipWith scale cs dat) `add` am + +-- | pseudorandom matrix with uniform elements between 0 and 1 +randm :: RandDist + -> Int -- ^ rows + -> Int -- ^ columns + -> IO (Matrix Double) +randm d r c = do + seed <- randomIO + return (reshape c $ randomVector seed d (r*c)) + +-- | pseudorandom matrix with uniform elements between 0 and 1 +rand :: Int -> Int -> IO (Matrix Double) +rand = randm Uniform + +{- | pseudorandom matrix with normal elements + +>>> disp 3 =<< randn 3 5 +3x5 +0.386 -1.141 0.491 -0.510 1.512 +0.069 -0.919 1.022 -0.181 0.745 +0.313 -0.670 -0.097 -1.575 -0.583 + +-} +randn :: Int -> Int -> IO (Matrix Double) +randn = randm Gaussian + 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 ( FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, - vectorScan, saveMatrix + vectorScan, saveMatrix, + Seed, RandDist(..), randomVector ) where import Data.Packed.Internal.Common @@ -308,7 +309,13 @@ vectorScan s = do foreign import ccall unsafe "saveMatrix" c_saveMatrix :: CString -> CString -> TM -saveMatrix :: FilePath -> String -> Matrix Double -> IO () +{- | save a matrix as a 2D ASCII table +-} +saveMatrix + :: FilePath + -> String -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.) + -> Matrix Double + -> IO () saveMatrix name format m = do cname <- newCString name cformat <- newCString format @@ -317,4 +324,23 @@ saveMatrix name format m = do free cformat return () +-------------------------------------------------------------------------------- + +type Seed = Int + +data RandDist = Uniform -- ^ uniform distribution in [0,1) + | Gaussian -- ^ normal distribution with mean zero and standard deviation one + deriving Enum + +-- | Obtains a vector of pseudorandom elements (use randomIO to get a random seed). +randomVector :: Seed + -> RandDist -- ^ distribution + -> Int -- ^ vector size + -> Vector Double +randomVector seed dist n = unsafePerformIO $ do + r <- createVector n + app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" + return r + +foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV -- cgit v1.2.3