3 files changed, 262 insertions, 0 deletions

diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs
new file mode 100644
index 0000000..1fde621
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs
@@ -0,0 +1,136 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.GSL
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.GSL (
+    RandDist(..), randomVector,
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
+) where
+
+import Data.Packed
+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
+import Foreign.C.String(newCString)
+import System.IO.Unsafe(unsafePerformIO)
+import System.Process(readProcess)
+
+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
+             -> 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
+
+--------------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
+
+--------------------------------------------------------------------------------
+
+-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
+fscanfVector :: FilePath -> Int -> IO (Vector Double)
+fscanfVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
+fprintfVector :: FilePath -> String -> Vector Double -> IO ()
+fprintfVector filename fmt v = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
+
+-- | Loads a vector from a binary file (the number of elements must be known in advance).
+freadVector :: FilePath -> Int -> IO (Vector Double)
+freadVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector to a binary file.
+fwriteVector :: FilePath -> Vector Double -> IO ()
+fwriteVector filename v = do
+    charname <- newCString filename
+    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
+    free charname
+
+foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
+
+type PD = Ptr Double                            --
+type TV = CInt -> PD -> IO CInt                 --
+type TM = CInt -> CInt -> PD -> IO CInt         --
+
+--------------------------------------------------------------------------------
+
+{- |  obtains the number of rows and columns in an ASCII data file
+      (provisionally using unix's wc).
+-}
+fileDimensions :: FilePath -> IO (Int,Int)
+fileDimensions fname = do
+    wcres <- readProcess "wc" ["-w",fname] ""
+    contents <- readFile fname
+    let tot = read . head . words $ wcres
+        c   = length . head . dropWhile null . map words . lines $ contents
+    if tot > 0
+        then return (tot `div` c, c)
+        else return (0,0)
+
+-- | Loads a matrix from an ASCII file formatted as a 2D table.
+loadMatrix :: FilePath -> IO (Matrix Double)
+loadMatrix file = fromFile file =<< fileDimensions file
+
+-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
new file mode 100644
index 0000000..d2278ee
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
@@ -0,0 +1,42 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.IO
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.IO (
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
+) where
+
+import Data.Packed
+import Numeric.GSL.LinearAlgebra
+import System.Process(readProcess)
+
+
+{- |  obtains the number of rows and columns in an ASCII data file
+      (provisionally using unix's wc).
+-}
+fileDimensions :: FilePath -> IO (Int,Int)
+fileDimensions fname = do
+    wcres <- readProcess "wc" ["-w",fname] ""
+    contents <- readFile fname
+    let tot = read . head . words $ wcres
+        c   = length . head . dropWhile null . map words . lines $ contents
+    if tot > 0
+        then return (tot `div` c, c)
+        else return (0,0)
+
+-- | Loads a matrix from an ASCII file formatted as a 2D table.
+loadMatrix :: FilePath -> IO (Matrix Double)
+loadMatrix file = fromFile file =<< fileDimensions file
+
+-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
new file mode 100644
index 0000000..e9ab2b7
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
@@ -0,0 +1,84 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.Random
+-- Copyright   :  (c) Alberto Ruiz 2009-14
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Random vectors and matrices.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Random (
+    Seed,
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    uniformSample,
+    rand, randn
+) where
+
+import Numeric.GSL.LinearAlgebra
+import Numeric.Container
+import Numeric.LinearAlgebra.Algorithms
+import System.Random(randomIO)
+
+
+type Seed = Int
+
+-- | 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
+
+>>> x <- randn 3 5
+>>> disp 3 x
+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
+