-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.LinearAlgebra.Random
-- Copyright   :  (c) Alberto Ruiz 2009-14
-- License     :  BSD3
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-- Random vectors and matrices.
--
-----------------------------------------------------------------------------

module Internal.Random (
    Seed,
    RandDist(..),
    randomVector,
    gaussianSample,
    uniformSample,
    rand, randn
) where

import Internal.Vectorized
import Internal.Vector
import Internal.Matrix
import Internal.Numeric
import Internal.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
               -> Herm Double   -- ^ covariance matrix
               -> Matrix Double -- ^ result
gaussianSample :: Seed -> Seed -> Vector Double -> Herm Double -> Matrix Double
gaussianSample Seed
seed Seed
n Vector Double
med Herm Double
cov = Matrix Double
m where
    c :: Seed
c = forall t. Storable t => Vector t -> Seed
dim Vector Double
med
    meds :: Matrix Double
meds = forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
konst' Double
1 Seed
n forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
med
    rs :: Matrix Double
rs = forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
c forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
Gaussian (Seed
c forall a. Num a => a -> a -> a
* Seed
n)
    m :: Matrix Double
m = Matrix Double
rs forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall t. Field t => Herm t -> Matrix t
chol Herm Double
cov forall c. Additive c => c -> c -> c
`add` Matrix Double
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 -> Seed -> [(Double, Double)] -> Matrix Double
uniformSample Seed
seed Seed
n [(Double, Double)]
rgs = Matrix Double
m where
    ([Double]
as,[Double]
bs) = forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
rgs
    a :: Vector Double
a = forall a. Storable a => [a] -> Vector a
fromList [Double]
as
    cs :: [Double]
cs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract [Double]
as [Double]
bs
    d :: Seed
d = forall t. Storable t => Vector t -> Seed
dim Vector Double
a
    dat :: [Vector Double]
dat = forall t. Element t => Matrix t -> [Vector t]
toRows forall a b. (a -> b) -> a -> b
$ forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
n forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
Uniform (Seed
nforall a. Num a => a -> a -> a
*Seed
d)
    am :: Matrix Double
am = forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
konst' Double
1 Seed
n forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
a
    m :: Matrix Double
m = forall t. Element t => [Vector t] -> Matrix t
fromColumns (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale [Double]
cs [Vector Double]
dat) forall c. Additive c => c -> c -> c
`add` Matrix Double
am

-- | pseudorandom matrix with uniform elements between 0 and 1
randm :: RandDist
     -> Int -- ^ rows
     -> Int -- ^ columns
     -> IO (Matrix Double)
randm :: RandDist -> Seed -> Seed -> IO (Matrix Double)
randm RandDist
d Seed
r Seed
c = do
    Seed
seed <- forall a (m :: * -> *). (Random a, MonadIO m) => m a
randomIO
    forall (m :: * -> *) a. Monad m => a -> m a
return (forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
c forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
d (Seed
rforall a. Num a => a -> a -> a
*Seed
c))

-- | pseudorandom matrix with uniform elements between 0 and 1
rand :: Int -> Int -> IO (Matrix Double)
rand :: Seed -> Seed -> IO (Matrix Double)
rand = RandDist -> Seed -> Seed -> IO (Matrix Double)
randm RandDist
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 :: Seed -> Seed -> IO (Matrix Double)
randn = RandDist -> Seed -> Seed -> IO (Matrix Double)
randm RandDist
Gaussian