algorithms interface reorganized

3 files changed, 69 insertions, 106 deletions

diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 162426f..a67f822 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -9,106 +9,114 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses ffi
 
-A simple interface to the available matrix computations. We have defined some generic functions 
-on both real and complex matrices in such a way that higher level algorithms and 
-testing properties are valid for both base types.
+A generic interface for a number of essential functions. Using it some higher level algorithms
+and testing properties can be written for both real and complex matrices.
 
-In any case the specific functions for particular base types can also be explicitly
+In any case, the specific functions for particular base types can also be explicitly
 imported from the LAPACK and GSL.Matrix modules.
 
 -}
 -----------------------------------------------------------------------------
 
 module LinearAlgebra.Algorithms (
--- * Basic Linear Algebra
-    scale,
-    add,
-    multiply, dot, outer,
-    linearSolve, linearSolveSVD,
+-- * Linear Systems
+    linearSolve,
+    inv, pinv,
+    pinvTol, det,
 -- * Matrix factorizations
-    svd, lu, eig, eigSH,
-    qr, chol,
--- * Utilities
-    Normed(..), NormType(..),
-    det,inv,pinv,full,economy,
-    pinvTol,
---    pinvTolg,
+-- ** Singular value decomposition
+    svd,
+    full, economy,
+-- ** Eigensystems
+    eig, LinearAlgebra.Algorithms.eigS, LinearAlgebra.Algorithms.eigH,
+-- ** Other 
+    LinearAlgebra.Algorithms.qr, chol,
+-- * Nullspace
     nullspacePrec,
     nullVector,
--- * Conversions
-    toComplex, fromComplex,
-    comp, real, complex,
-    conj, ctrans,
 -- * Misc
     eps, i,
-    scaleRecip,
-    addConstant,
-    sub,
-    mul,
-    divide,
+    ctrans,
+    Normed(..), NormType(..),
+    GenMat(linearSolveSVD,lu,eigSH)
 ) where
 
 
 import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
-import Data.Packed.Matrix
+import Data.Packed
 import GSL.Matrix
 import GSL.Vector
 import LAPACK
 import Complex
 import LinearAlgebra.Linear
 
--- | Base types for which some optimized matrix computations are available
-class (Linear Matrix t) => Optimized t where
+-- | matrix computations available for both real and complex matrices
+class (Linear Matrix t) => GenMat t where
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)
     linearSolve :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
-    ctrans :: Matrix t -> Matrix t
     eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH       :: Matrix t -> (Vector Double, Matrix t)
     chol        :: Matrix t -> Matrix t
+    -- | conjugate transpose
+    ctrans :: Matrix t -> Matrix t
 
-instance Optimized Double where
+instance GenMat Double where
     svd = svdR
     lu  = luR
     linearSolve = linearSolveR
     linearSolveSVD = linearSolveSVDR Nothing
     ctrans = trans
     eig = eigR
-    eigSH = eigS
+    eigSH = LAPACK.eigS
     chol = cholR
 
-instance Optimized (Complex Double) where
+instance GenMat (Complex Double) where
     svd = svdC
     lu  = luC
     linearSolve = linearSolveC
     linearSolveSVD = linearSolveSVDC Nothing
     ctrans = conjTrans
     eig = eigC
-    eigSH = eigH
+    eigSH =  LAPACK.eigH
     chol = error "cholC not yet implemented" -- waiting for GSL-1.10
 
+-- | eigensystem of a symmetric matrix
+eigS :: Matrix Double -> (Vector Double, Matrix Double)
+eigS = LAPACK.eigS
+
+-- | eigensystem of a hermitian matrix
+eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
+eigH = LAPACK.eigH
+
+qr :: Matrix Double -> (Matrix Double, Matrix Double)
+qr = GSL.Matrix.qr
+
 square m = rows m == cols m
 
-det :: Optimized t => Matrix t -> t
+det :: GenMat t => Matrix t -> t
 det m | square m = s * (product $ toList $ takeDiag $ u)
       | otherwise = error "det of nonsquare matrix"
     where (_,u,_,s) = lu m
 
-inv :: Optimized t => Matrix t -> Matrix t
+inv :: GenMat t => Matrix t -> Matrix t
 inv m | square m = m `linearSolve` ident (rows m)
       | otherwise = error "inv of nonsquare matrix"
 
-pinv :: Optimized t => Matrix t -> Matrix t
+pinv :: GenMat t => Matrix t -> Matrix t
 pinv m = linearSolveSVD m (ident (rows m))
 
-
+full :: Field t 
+     => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
 full svd m = (u, d ,v) where
     (u,s,v) = svd m
     d = diagRect s r c
     r = rows m
     c = cols m
 
+economy :: Field t 
+        => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
 economy svd m = (u', subVector 0 d s, v') where
     (u,s,v) = svd m
     sl@(g:_) = toList (complex s)
@@ -123,39 +131,25 @@ economy svd m = (u', subVector 0 d s, v') where
     v' = takeColumns d v
 
 
-{- | Machine precision of a Double.
-
->> eps
-> 2.22044604925031e-16
-
-(The value used by GNU-Octave)
-
--}
+-- | The machine precision of a Double: @eps == 2.22044604925031e-16@ (the value used by GNU-Octave).
 eps :: Double
 eps =  2.22044604925031e-16
 
-{- | The imaginary unit
-
-@> 'ident' 3 \<\> i
-1.i   0.   0.
- 0.  1.i   0.
- 0.   0.  1.i@
-
--}
+-- | The imaginary unit: @i == 0.0 :+ 1.0@
 i :: Complex Double
 i = 0:+1
 
 
 -- | matrix product
-mXm :: (Num t, Optimized t) => Matrix t -> Matrix t -> Matrix t
+mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t
 mXm = multiply
 
 -- | matrix - vector product
-mXv :: (Num t, Optimized t) => Matrix t -> Vector t -> Vector t
+mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t
 mXv m v = flatten $ m `mXm` (asColumn v)
 
 -- | vector - matrix product
-vXm :: (Num t, Optimized t) => Vector t -> Matrix t -> Vector t
+vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t
 vXm v m = flatten $ (asRow v) `mXm` m
 
 
@@ -167,15 +161,6 @@ norm2 = toScalarR Norm2
 norm1 :: Vector Double -> Double
 norm1 = toScalarR AbsSum
 
-vectorMax :: Vector Double -> Double
-vectorMax = toScalarR Max
-vectorMin :: Vector Double -> Double
-vectorMin = toScalarR Min
-vectorMaxIndex :: Vector Double -> Int
-vectorMaxIndex = round . toScalarR MaxIdx
-vectorMinIndex :: Vector Double -> Int
-vectorMinIndex = round . toScalarR MinIdx
-
 data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int
 
 pnormRV PNorm2 = norm2
@@ -199,7 +184,7 @@ pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` const
 --pnormCM _ _ = error "p norm not yet defined"
 
 -- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'.
---pnorm :: (Container t, Optimized a) => Int -> t a -> Double
+--pnorm :: (Container t, GenMat a) => Int -> t a -> Double
 --pnorm = pnormG
 
 class Normed t where
@@ -222,7 +207,7 @@ instance Normed (Matrix (Complex Double)) where
 -----------------------------------------------------------------------
 
 -- | The nullspace of a matrix from its SVD decomposition.
-nullspacePrec :: Optimized t
+nullspacePrec :: GenMat t
               => Double          -- ^ relative tolerance in 'eps' units
               -> Matrix t   -- ^ input matrix
               -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
@@ -235,12 +220,12 @@ nullspacePrec t m = ns where
     ns = drop rank $ toRows $ ctrans v
 
 -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
-nullVector :: Optimized t => Matrix t -> Vector t
+nullVector :: GenMat t => Matrix t -> Vector t
 nullVector = last . nullspacePrec 1
 
 ------------------------------------------------------------------------
 
-{- | Pseudoinverse of a real matrix with the desired tolerance, expressed as a
+{-  Pseudoinverse of a real matrix with the desired tolerance, expressed as a
 multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
 
 @\> let m = 'fromLists' [[1,0,    0]
@@ -258,7 +243,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
 0. 0. 1.@
 
 -}
-pinvTol :: Double -> Matrix Double -> Matrix Double
+--pinvTol :: Double -> Matrix Double -> Matrix Double
 pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
     (u,s,v) = svdR m
     sl@(g:_) = toList s
diff --git a/lib/LinearAlgebra/Interface.hs b/lib/LinearAlgebra/Interface.hs
index 3392d54..0c65a8b 100644
--- a/lib/LinearAlgebra/Interface.hs
+++ b/lib/LinearAlgebra/Interface.hs
@@ -16,6 +16,7 @@ Operators for frequent operations.
 
 module LinearAlgebra.Interface(
     (<>),(<.>),
+    (<\>),
     (.*),(*/),
     (<|>),(<->),
 ) where
@@ -23,6 +24,7 @@ module LinearAlgebra.Interface(
 import LinearAlgebra.Linear
 import Data.Packed.Vector
 import Data.Packed.Matrix
+import LinearAlgebra.Algorithms
 
 class Mul a b c | a b -> c where
  infixl 7 <>
@@ -59,6 +61,11 @@ a .* x = scale a x
 infixl 7 */
 v */ x = scale (recip x) v
 
+-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
+(<\>) :: (GenMat a) => Matrix a -> Vector a -> Vector a
+infixl 7 <\>
+m <\> v = flatten (linearSolveSVD m (reshape 1 v))
+
 ------------------------------------------------
 
 class Joinable a b where
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index 85daa4a..3e6c55d 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -9,10 +9,10 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses ffi
 
+Basic optimized operations on vectors and matrices.
 
 -}
 -----------------------------------------------------------------------------
--- #hide
 
 module LinearAlgebra.Linear (
     Linear(..),
@@ -21,25 +21,22 @@ module LinearAlgebra.Linear (
 
 
 import Data.Packed.Internal
-import Data.Packed.Matrix
+import Data.Packed
 import GSL.Vector
 import Complex
 
 -- | basic optimized operations
-class (Field e) => Linear c e where
+class (Container c e) => Linear c e where
     scale       :: e -> c e -> c e
-    scaleRecip  :: e -> c e -> c e
     addConstant :: e -> c e -> c e
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
+    -- | element by element multiplication
     mul         :: c e -> c e -> c e
+    -- | element by element division
     divide      :: c e -> c e -> c e
-    toComplex   :: RealFloat e => (c e, c e) -> c (Complex e)
-    fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
-    comp        :: RealFloat e => c e -> c (Complex e)
-    conj        :: RealFloat e => c (Complex e) -> c (Complex e)
-    real        :: c Double -> c e
-    complex     :: c e -> c (Complex Double)
+    -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
+    scaleRecip  :: e -> c e -> c e
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
@@ -49,12 +46,6 @@ instance Linear Vector Double where
     sub = vectorZipR Sub
     mul = vectorZipR Mul
     divide = vectorZipR Div
-    toComplex = Data.Packed.Internal.toComplex
-    fromComplex = Data.Packed.Internal.fromComplex
-    comp = Data.Packed.Internal.comp
-    conj = Data.Packed.Internal.conj
-    real = id
-    complex = LinearAlgebra.Linear.comp
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
@@ -64,12 +55,6 @@ instance Linear Vector (Complex Double) where
     sub = vectorZipC Sub
     mul = vectorZipC Mul
     divide = vectorZipC Div
-    toComplex = undefined -- can't match
-    fromComplex = undefined
-    comp = undefined
-    conj = undefined
-    real = LinearAlgebra.Linear.comp
-    complex = id
 
 instance Linear Matrix Double where
     scale x = liftMatrix (scale x)
@@ -79,14 +64,6 @@ instance Linear Matrix Double where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
-    toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex
-    fromComplex z = (reshape c r, reshape c i)
-        where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z)
-              c = cols z
-    comp = liftMatrix Data.Packed.Internal.comp
-    conj = liftMatrix Data.Packed.Internal.conj
-    real = id
-    complex = LinearAlgebra.Linear.comp
 
 instance Linear Matrix (Complex Double) where
     scale x = liftMatrix (scale x)
@@ -96,12 +73,6 @@ instance Linear Matrix (Complex Double) where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
-    toComplex = undefined
-    fromComplex = undefined
-    comp = undefined
-    conj = undefined
-    real = LinearAlgebra.Linear.comp
-    complex = id
 
 --------------------------------------------------