    -----------------------------------------------------------------------------
    -- |
    -- 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 = Vector Double -> Seed
    forall t. Storable t => Vector t -> Seed
    dim Vector Double
    med
        meds :: Matrix Double
    meds = Double -> IndexOf Vector -> Vector Double
    forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
    konst' Double
    1 Seed
    IndexOf Vector
    n Vector Double -> Vector Double -> Matrix Double
    forall t. Product t => Vector t -> Vector t -> Matrix t
    `outer` Vector Double
    med
        rs :: Matrix Double
    rs = Seed -> Vector Double -> Matrix Double
    forall t. Storable t => Seed -> Vector t -> Matrix t
    reshape Seed
    c (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
    forall a b. (a -> b) -> a -> b
    $ Seed -> RandDist -> Seed -> Vector Double
    randomVector Seed
    seed RandDist
    Gaussian (Seed
    c Seed -> Seed -> Seed
    forall a. Num a => a -> a -> a
    * Seed
    n)
        m :: Matrix Double
    m = Matrix Double
    rs Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Herm Double -> Matrix Double
    forall t. Field t => Herm t -> Matrix t
    chol Herm Double
    cov Matrix Double -> Matrix Double -> Matrix Double
    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) = [(Double, Double)] -> ([Double], [Double])
    forall a b. [(a, b)] -> ([a], [b])
    unzip [(Double, Double)]
    rgs
        a :: Vector Double
    a = [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList [Double]
    as
        cs :: [Double]
    cs = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Double -> Double -> Double
    forall a. Num a => a -> a -> a
    subtract [Double]
    as [Double]
    bs
        d :: Seed
    d = Vector Double -> Seed
    forall t. Storable t => Vector t -> Seed
    dim Vector Double
    a
        dat :: [Vector Double]
    dat = Matrix Double -> [Vector Double]
    forall t. Element t => Matrix t -> [Vector t]
    toRows (Matrix Double -> [Vector Double])
    -> Matrix Double -> [Vector Double]
    forall a b. (a -> b) -> a -> b
    $ Seed -> Vector Double -> Matrix Double
    forall t. Storable t => Seed -> Vector t -> Matrix t
    reshape Seed
    n (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
    forall a b. (a -> b) -> a -> b
    $ Seed -> RandDist -> Seed -> Vector Double
    randomVector Seed
    seed RandDist
    Uniform (Seed
    nSeed -> Seed -> Seed
    forall a. Num a => a -> a -> a
    *Seed
    d)
        am :: Matrix Double
    am = Double -> IndexOf Vector -> Vector Double
    forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
    konst' Double
    1 Seed
    IndexOf Vector
    n Vector Double -> Vector Double -> Matrix Double
    forall t. Product t => Vector t -> Vector t -> Matrix t
    `outer` Vector Double
    a
        m :: Matrix Double
    m = [Vector Double] -> Matrix Double
    forall t. Element t => [Vector t] -> Matrix t
    fromColumns ((Double -> Vector Double -> Vector Double)
    -> [Double] -> [Vector Double] -> [Vector Double]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Double -> Vector Double -> Vector Double
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale [Double]
    cs [Vector Double]
    dat) Matrix Double -> Matrix Double -> Matrix Double
    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 <- IO Seed
    forall a (m :: * -> *). (Random a, MonadIO m) => m a
    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 :: 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