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 +++++++- packages/hmatrix/src/Numeric/Container.hs | 2 +- packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs | 61 ++--------------- packages/hmatrix/src/Numeric/GSL/gsl-aux.c | 4 +- .../hmatrix/src/Numeric/LinearAlgebra/Random.hs | 4 +- 11 files changed, 184 insertions(+), 78 deletions(-) create mode 100644 packages/base/src/Numeric/LinearAlgebra/Random.hs (limited to 'packages') 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 diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs index b6e797b..2b61d90 100644 --- a/packages/hmatrix/src/Numeric/Container.hs +++ b/packages/hmatrix/src/Numeric/Container.hs @@ -16,7 +16,7 @@ module Numeric.Container ( meanCov ) where -import Data.Packed.Numeric +import Data.Packed.Numeric hiding (saveMatrix, loadMatrix) import Numeric.LinearAlgebra.IO import Numeric.LinearAlgebra.Random hiding (Seed) import Numeric.LinearAlgebra.Util(meanCov) diff --git a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs index 43e46c5..8bd0cc2 100644 --- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs +++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs @@ -9,15 +9,15 @@ ----------------------------------------------------------------------------- module Numeric.GSL.LinearAlgebra ( - RandDist(..), randomVector, + randomVector, saveMatrix, fwriteVector, freadVector, fprintfVector, fscanfVector ) where import Data.Packed +import Numeric.LinearAlgebra.Base(RandDist(..)) import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) -import Data.Complex import Foreign.Marshal.Alloc(free) import Foreign.Ptr(Ptr) import Foreign.C.Types @@ -28,10 +28,6 @@ fromei x = fromIntegral (fromEnum x) :: CInt ----------------------------------------------------------------------- -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 from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. randomVector :: Int -- ^ seed -> RandDist -- ^ distribution @@ -39,10 +35,10 @@ randomVector :: Int -- ^ seed -> Vector Double randomVector seed dist n = unsafePerformIO $ do r <- createVector n - app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" + app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" return r -foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV +foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV -------------------------------------------------------------------------------- @@ -105,56 +101,7 @@ fwriteVector filename v = do foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV -type PF = Ptr Float -- type PD = Ptr Double -- -type PQ = Ptr (Complex Float) -- -type PC = Ptr (Complex Double) -- -type TF = CInt -> PF -> IO CInt -- -type TFF = CInt -> PF -> TF -- -type TFV = CInt -> PF -> TV -- -type TVF = CInt -> PD -> TF -- -type TFFF = CInt -> PF -> TFF -- type TV = CInt -> PD -> IO CInt -- -type TVV = CInt -> PD -> TV -- -type TVVV = CInt -> PD -> TVV -- -type TFM = CInt -> CInt -> PF -> IO CInt -- -type TFMFM = CInt -> CInt -> PF -> TFM -- -type TFMFMFM = CInt -> CInt -> PF -> TFMFM -- type TM = CInt -> CInt -> PD -> IO CInt -- -type TMM = CInt -> CInt -> PD -> TM -- -type TVMM = CInt -> PD -> TMM -- -type TMVMM = CInt -> CInt -> PD -> TVMM -- -type TMMM = CInt -> CInt -> PD -> TMM -- -type TVM = CInt -> PD -> TM -- -type TVVM = CInt -> PD -> TVM -- -type TMV = CInt -> CInt -> PD -> TV -- -type TMMV = CInt -> CInt -> PD -> TMV -- -type TMVM = CInt -> CInt -> PD -> TVM -- -type TMMVM = CInt -> CInt -> PD -> TMVM -- -type TCM = CInt -> CInt -> PC -> IO CInt -- -type TCVCM = CInt -> PC -> TCM -- -type TCMCVCM = CInt -> CInt -> PC -> TCVCM -- -type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM -- -type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM -- -type TCMCM = CInt -> CInt -> PC -> TCM -- -type TVCM = CInt -> PD -> TCM -- -type TCMVCM = CInt -> CInt -> PC -> TVCM -- -type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM -- -type TCMCMCM = CInt -> CInt -> PC -> TCMCM -- -type TCV = CInt -> PC -> IO CInt -- -type TCVCV = CInt -> PC -> TCV -- -type TCVCVCV = CInt -> PC -> TCVCV -- -type TCVV = CInt -> PC -> TV -- -type TQV = CInt -> PQ -> IO CInt -- -type TQVQV = CInt -> PQ -> TQV -- -type TQVQVQV = CInt -> PQ -> TQVQV -- -type TQVF = CInt -> PQ -> TF -- -type TQM = CInt -> CInt -> PQ -> IO CInt -- -type TQMQM = CInt -> CInt -> PQ -> TQM -- -type TQMQMQM = CInt -> CInt -> PQ -> TQMQM -- -type TCMCV = CInt -> CInt -> PC -> TCV -- -type TVCV = CInt -> PD -> TCV -- -type TCVM = CInt -> PC -> TM -- -type TMCVM = CInt -> CInt -> PD -> TCVM -- -type TMMCVM = CInt -> CInt -> PD -> TMCVM -- diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c index 5da94ca..ffc5c20 100644 --- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c +++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c @@ -924,8 +924,8 @@ int nlfit(int method, int f(int, double*, int, double*), #define RAN(C,F) case C: { for(k=0;k