GSL.LinearAlgebra reorganized

6 files changed, 183 insertions, 263 deletions

diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs
deleted file mode 100644
index 146780d..0000000
--- a/packages/hmatrix/src/Numeric/Container.hs
+++ /dev/null
@@ -1,241 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FunctionalDependencies #-}
-{-# LANGUAGE UndecidableInstances #-}
-
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.Container
--- Copyright   :  (c) Alberto Ruiz 2010-14
--- License     :  BSD3
--- Maintainer  :  Alberto Ruiz
--- Stability   :  provisional
--- Portability :  portable
---
--- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
---
--- The 'Container' class is used to define optimized generic functions which work
--- on 'Vector' and 'Matrix' with real or complex elements.
---
--- Some of these functions are also available in the instances of the standard
--- numeric Haskell classes provided by "Numeric.LinearAlgebra".
---
------------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
-
-module Numeric.Container (
-    -- * Basic functions
-    module Data.Packed,
-    konst, build,
-    linspace,
-    diag, ident,
-    ctrans,
-    -- * Generic operations
-    Container(..),
-    -- * Matrix product
-    Product(..), udot, dot, (◇),
-    Mul(..),
-    Contraction(..),
-    optimiseMult,
-    mXm,mXv,vXm,LSDiv(..),
-    outer, kronecker,
-    -- * Element conversion
-    Convert(..),
-    Complexable(),
-    RealElement(),
-
-    RealOf, ComplexOf, SingleOf, DoubleOf,
-
-    IndexOf,
-    module Data.Complex
-) where
-
-import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
-import Data.Packed.Numeric
-import Data.Complex
-import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
-import Data.Monoid(Monoid(mconcat))
-
-------------------------------------------------------------------
-
-{- | Creates a real vector containing a range of values:
-
->>> linspace 5 (-3,7::Double)
-fromList [-3.0,-0.5,2.0,4.5,7.0]@
-
->>> linspace 5 (8,2+i) :: Vector (Complex Double)
-fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
-
-Logarithmic spacing can be defined as follows:
-
-@logspace n (a,b) = 10 ** linspace n (a,b)@
--}
-linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
-linspace 0 (a,b) = fromList[(a+b)/2]
-linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
-    where s = (b-a)/fromIntegral (n-1)
-
---------------------------------------------------------
-
-class Contraction a b c | a b -> c
-  where
-    infixl 7 <.>
-    {- | Matrix product, matrix vector product, and dot product
-
-Examples:
-
->>> let a = (3><4)   [1..]      :: Matrix Double
->>> let v = fromList [1,0,2,-1] :: Vector Double
->>> let u = fromList [1,2,3]    :: Vector Double
-
->>> a
-(3><4)
- [ 1.0,  2.0,  3.0,  4.0
- , 5.0,  6.0,  7.0,  8.0
- , 9.0, 10.0, 11.0, 12.0 ]
-
-matrix × matrix:
-
->>> disp 2 (a <.> trans a)
-3x3
- 30   70  110
- 70  174  278
-110  278  446
-
-matrix × vector:
-
->>> a <.> v
-fromList [3.0,11.0,19.0]
-
-dot product:
-
->>> u <.> fromList[3,2,1::Double]
-10
-
-For complex vectors the first argument is conjugated:
-
->>> fromList [1,i] <.> fromList[2*i+1,3]
-1.0 :+ (-1.0)
-
->>> fromList [1,i,1-i] <.> complex a
-fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
-
--}
-    (<.>) :: a -> b -> c
-
-
-instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
-    u <.> v = conj u `udot` v
-
-instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
-    (<.>) = mXv
-
-instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
-    (<.>) v m = (conj v) `vXm` m
-
-instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
-    (<.>) = mXm
-
-
---------------------------------------------------------------------------------
-
-class Mul a b c | a b -> c where
- infixl 7 <>
- -- | Matrix-matrix, matrix-vector, and vector-matrix products.
- (<>)  :: Product t => a t -> b t -> c t
-
-instance Mul Matrix Matrix Matrix where
-    (<>) = mXm
-
-instance Mul Matrix Vector Vector where
-    (<>) m v = flatten $ m <> asColumn v
-
-instance Mul Vector Matrix Vector where
-    (<>) v m = flatten $ asRow v <> m
-
---------------------------------------------------------------------------------
-
-class LSDiv c where
- infixl 7 <\>
- -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
- (<\>)  :: Field t => Matrix t -> c t -> c t
-
-instance LSDiv Vector where
-    m <\> v = flatten (linearSolveSVD m (reshape 1 v))
-
-instance LSDiv Matrix where
-    (<\>) = linearSolveSVD
-
---------------------------------------------------------------------------------
-
-class Konst e d c | d -> c, c -> d
-  where
-    -- |
-    -- >>> konst 7 3 :: Vector Float
-    -- fromList [7.0,7.0,7.0]
-    --
-    -- >>> konst i (3::Int,4::Int)
-    -- (3><4)
-    --  [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
-    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
-    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
-    --
-    konst :: e -> d -> c e
-
-instance Container Vector e => Konst e Int Vector
-  where
-    konst = konst'
-
-instance Container Vector e => Konst e (Int,Int) Matrix
-  where
-    konst = konst'
-
---------------------------------------------------------------------------------
-
-class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
-  where
-    -- |
-    -- >>> build 5 (**2) :: Vector Double
-    -- fromList [0.0,1.0,4.0,9.0,16.0]
-    --
-    -- Hilbert matrix of order N:
-    --
-    -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
-    -- >>> putStr . dispf 2 $ hilb 3
-    -- 3x3
-    -- 1.00  0.50  0.33
-    -- 0.50  0.33  0.25
-    -- 0.33  0.25  0.20
-    --
-    build :: d -> f -> c e
-
-instance Container Vector e => Build Int (e -> e) Vector e
-  where
-    build = build'
-
-instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
-  where
-    build = build'
-
---------------------------------------------------------------------------------
-
-{- | alternative operator for '(\<.\>)'
-
-x25c7, white diamond
-
--}
-(◇) :: Contraction a b c => a -> b -> c
-infixl 7 ◇
-(◇) = (<.>)
-
--- | dot product: @cdot u v = 'udot' ('conj' u) v@
-dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
-dot u v = udot (conj u) v
-
---------------------------------------------------------------------------------
-
-optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
-optimiseMult = mconcat
-
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
index bca074f..43e46c5 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
@@ -1,16 +1,14 @@
 -----------------------------------------------------------------------------
 -- |
--- Module      :  Numeric.GSL.Vector
--- Copyright   :  (c) Alberto Ruiz 2007
+-- Module      :  Numeric.GSL.LinearAlgebra
+-- Copyright   :  (c) Alberto Ruiz 2007-14
 -- License     :  GPL
 -- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
 --
--- Low level interface to vector operations.
---
 -----------------------------------------------------------------------------
 
-module Numeric.GSL.Vector (
+module Numeric.GSL.LinearAlgebra (
     RandDist(..), randomVector,
     saveMatrix,
     fwriteVector, freadVector, fprintfVector, fscanfVector
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/LinearAlgebra.hs
new file mode 100644
index 0000000..1c2b1bb
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra.hs
@@ -0,0 +1,33 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+
+This module reexports all normally required functions for Linear Algebra applications.
+
+It also provides instances of standard classes 'Show', 'Read', 'Eq',
+'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
+In arithmetic operations one-component vectors and matrices automatically
+expand to match the dimensions of the other operand.
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra (
+    module Numeric.Container,
+    module Numeric.LinearAlgebra.Algorithms,
+    module Numeric.LinearAlgebra.IO,
+    module Numeric.LinearAlgebra.Random
+) where
+
+import Numeric.Container
+import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Util()
+import Numeric.LinearAlgebra.IO
+import Numeric.LinearAlgebra.Random
+
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/IO.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
index 08d812b..d2278ee 100644
--- a/packages/hmatrix/src/Numeric/IO.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
@@ -1,28 +1,22 @@
 -----------------------------------------------------------------------------
 -- |
--- Module      :  Numeric.IO
--- Copyright   :  (c) Alberto Ruiz 2010
+-- Module      :  Numeric.LinearAlgebra.IO
+-- Copyright   :  (c) Alberto Ruiz 2007-14
 -- License     :  GPL
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable
---
--- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'
 --
 -----------------------------------------------------------------------------
 
-module Numeric.IO (
-    dispf, disps, dispcf, vecdisp, latexFormat, format,
-    loadMatrix, saveMatrix, fromFile, fileDimensions,
-    readMatrix, fromArray2D,
-    fscanfVector, fprintfVector, freadVector, fwriteVector
+module Numeric.LinearAlgebra.IO (
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
 ) where
 
 import Data.Packed
-import Data.Packed.IO
+import Numeric.GSL.LinearAlgebra
 import System.Process(readProcess)
-import Numeric.GSL.Vector
 
 
 {- |  obtains the number of rows and columns in an ASCII data file
diff --git a/packages/hmatrix/src/Numeric/Random.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
index 0e722fd..e9ab2b7 100644
--- a/packages/hmatrix/src/Numeric/Random.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
@@ -1,6 +1,6 @@
 -----------------------------------------------------------------------------
 -- |
--- Module      :  Numeric.Random
+-- Module      :  Numeric.LinearAlgebra.Random
 -- Copyright   :  (c) Alberto Ruiz 2009-14
 -- License     :  GPL
 --
@@ -11,7 +11,7 @@
 --
 -----------------------------------------------------------------------------
 
-module Numeric.Random (
+module Numeric.LinearAlgebra.Random (
     Seed,
     RandDist(..),
     randomVector,
@@ -20,7 +20,7 @@ module Numeric.Random (
     rand, randn
 ) where
 
-import Numeric.GSL.Vector
+import Numeric.GSL.LinearAlgebra
 import Numeric.Container
 import Numeric.LinearAlgebra.Algorithms
 import System.Random(randomIO)