3 files changed, 120 insertions, 10 deletions

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