random numbers in base

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
                         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<rn; k++) {
+            rp[k] = (double)rand()/RAND_MAX;
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand();
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
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
     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