Numeric.Container and Numeric.LinearAlgebra moved to base

5 files changed, 267 insertions, 6 deletions

diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
index 0a92373..93fb281 100644
--- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
@@ -116,7 +116,7 @@ err (model,deriv) dat vsol = zip sol errs where
     dof = length dat - (rows cov)
     chi = norm2 (fromList $ cost (resMs model) dat sol)
     js = fromLists $ jacobian (resDs deriv) dat sol
-    cov = inv $ trans js <> js
+    cov = inv $ trans js <.> js
     errs = toList $ scalar c * sqrt (takeDiag cov)
 
 
diff --git a/packages/hmatrix/src/Numeric/GSL/IO.hs b/packages/hmatrix/src/Numeric/GSL/IO.hs
new file mode 100644
index 0000000..0d6031a
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/IO.hs
@@ -0,0 +1,42 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.IO
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.IO (
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
+) where
+
+import Data.Packed
+import Numeric.GSL.Vector
+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/GSL/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
index 8bd0cc2..17e2258 100644
--- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
@@ -9,13 +9,13 @@
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.LinearAlgebra (
-    randomVector,
+    RandDist(..), randomVector,
     saveMatrix,
-    fwriteVector, freadVector, fprintfVector, fscanfVector
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
 ) where
 
 import Data.Packed
-import Numeric.LinearAlgebra.Base(RandDist(..))
 import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
 
 import Foreign.Marshal.Alloc(free)
@@ -23,11 +23,16 @@ 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
@@ -35,10 +40,10 @@ randomVector :: Int      -- ^ seed
              -> Vector Double
 randomVector seed dist n = unsafePerformIO $ do
     r <- createVector n
-    app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
+    app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
     return r
 
-foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
+foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
 
 --------------------------------------------------------------------------------
 
@@ -105,3 +110,26 @@ 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/GSL/Random.hs b/packages/hmatrix/src/Numeric/GSL/Random.hs
new file mode 100644
index 0000000..2872b17
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/Random.hs
@@ -0,0 +1,84 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.Random
+-- Copyright   :  (c) Alberto Ruiz 2009-14
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Random vectors and matrices.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Random (
+    Seed,
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    uniformSample,
+    rand, randn
+) where
+
+import Numeric.GSL.Vector
+import Numeric.Container
+import Numeric.LinearAlgebra(Seed,RandDist(..),cholSH)
+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
+
+>>> 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
+
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
new file mode 100644
index 0000000..af79f32
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs
@@ -0,0 +1,107 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Vector (
+    randomVector,
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector
+) where
+
+import Data.Packed
+import Numeric.LinearAlgebra(RandDist(..))
+import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
+
+import Foreign.Marshal.Alloc(free)
+import Foreign.Ptr(Ptr)
+import Foreign.C.Types
+import Foreign.C.String(newCString)
+import System.IO.Unsafe(unsafePerformIO)
+
+fromei x = fromIntegral (fromEnum x) :: CInt
+
+-----------------------------------------------------------------------
+
+-- | 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_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
+    return r
+
+foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: 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         --
+