random numbers in base

11 files changed, 184 insertions, 78 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
 
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<rn;k++) { rp[k]= F(gen); }; OK }
 
-int random_vector(int seed, int code, RVEC(r)) {
-    DEBUGMSG("random_vector")
+int random_vector_GSL(int seed, int code, RVEC(r)) {
+    DEBUGMSG("random_vector_GSL")
     static gsl_rng * gen = NULL;
     if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
     gsl_rng_set (gen, seed);
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
index 0a82d3f..fa125a0 100644
--- a/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
@@ -22,11 +22,11 @@ module Numeric.LinearAlgebra.Random (
 
 import Numeric.GSL.LinearAlgebra
 import Data.Packed.Numeric
-import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Base(Seed,RandDist(..),cholSH)
 import System.Random(randomIO)
 
 
-type Seed = Int
+
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- Gaussian distribution.