From db223fb5f9cd4adef54736812f796b48ecc289e6 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 31 Oct 2007 13:36:37 +0000 Subject: Field->Element, GenMat->Field --- lib/Data/Packed.hs | 2 +- lib/Data/Packed/Internal/Matrix.hs | 36 ++-- lib/Data/Packed/Matrix.hs | 44 ++--- lib/GSLHaskell.hs | 327 -------------------------------- lib/Numeric/LinearAlgebra/Algorithms.hs | 46 ++--- lib/Numeric/LinearAlgebra/Instances.hs | 8 +- lib/Numeric/LinearAlgebra/Interface.hs | 14 +- lib/Numeric/LinearAlgebra/Linear.hs | 6 +- 8 files changed, 78 insertions(+), 405 deletions(-) delete mode 100644 lib/GSLHaskell.hs (limited to 'lib') diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs index 668d2f7..53aced9 100644 --- a/lib/Data/Packed.hs +++ b/lib/Data/Packed.hs @@ -27,7 +27,7 @@ import Data.Complex import Data.Packed.Internal -- | conversion utilities -class (Field e) => Container c e where +class (Element e) => Container c e where 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) diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index f63ee52..fbab33c 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -84,17 +84,17 @@ type Mt t s = Int -> Int -> Ptr t -> s -- type t ::> s = Mt t s -- | the inverse of 'Data.Packed.Matrix.fromLists' -toLists :: (Field t) => Matrix t -> [[t]] +toLists :: (Element t) => Matrix t -> [[t]] toLists m = partit (cols m) . toList . cdat $ m -- | creates a Matrix from a list of vectors -fromRows :: Field t => [Vector t] -> Matrix t +fromRows :: Element t => [Vector t] -> Matrix t fromRows vs = case common dim vs of Nothing -> error "fromRows applied to [] or to vectors with different sizes" Just c -> reshape c (join vs) -- | extracts the rows of a matrix as a list of vectors -toRows :: Field t => Matrix t -> [Vector t] +toRows :: Element t => Matrix t -> [Vector t] toRows m = toRows' 0 where v = cdat m r = rows m @@ -103,11 +103,11 @@ toRows m = toRows' 0 where | otherwise = subVector k c v : toRows' (k+c) -- | Creates a matrix from a list of vectors, as columns -fromColumns :: Field t => [Vector t] -> Matrix t +fromColumns :: Element t => [Vector t] -> Matrix t fromColumns m = trans . fromRows $ m -- | Creates a list of vectors from the columns of a matrix -toColumns :: Field t => Matrix t -> [Vector t] +toColumns :: Element t => Matrix t -> [Vector t] toColumns m = toRows . trans $ m @@ -152,18 +152,18 @@ where r is the desired number of rows.) , 9.0, 10.0, 11.0, 12.0 ]@ -} -reshape :: Field t => Int -> Vector t -> Matrix t +reshape :: Element t => Int -> Vector t -> Matrix t reshape c v = matrixFromVector RowMajor c v singleton x = reshape 1 (fromList [x]) -- | application of a vector function on the flattened matrix elements -liftMatrix :: (Field a, Field b) => (Vector a -> Vector b) -> Matrix a -> Matrix b +liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d) liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) -- | application of a vector function on the flattened matrices elements -liftMatrix2 :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t +liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t liftMatrix2 f m1 m2 | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2" | otherwise = case m1 of @@ -176,8 +176,8 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 ---------------------------------------------------------------- --- | Optimized matrix computations are provided for elements in the Field class. -class (Storable a, Floating a) => Field a where +-- | Optimized matrix computations are provided for elements in the Element class. +class (Storable a, Floating a) => Element a where constantD :: a -> Int -> Vector a transdata :: Int -> Vector a -> Int -> Vector a multiplyD :: Matrix a -> Matrix a -> Matrix a @@ -186,14 +186,14 @@ class (Storable a, Floating a) => Field a where -> Matrix a -> Matrix a diagD :: Vector a -> Matrix a -instance Field Double where +instance Element Double where constantD = constantR transdata = transdataR multiplyD = multiplyR subMatrixD = subMatrixR diagD = diagR -instance Field (Complex Double) where +instance Element (Complex Double) where constantD = constantC transdata = transdataC multiplyD = multiplyC @@ -202,7 +202,7 @@ instance Field (Complex Double) where ------------------------------------------------------------------ -(>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a +(>|<) :: (Element a) => Int -> Int -> [a] -> Matrix a r >|< c = f where f l | dim v == r*c = matrixFromVector ColumnMajor c v | otherwise = error $ "inconsistent list size = " @@ -260,13 +260,13 @@ foreign import ccall safe "auxi.h multiplyC" -> Int -> Int -> Ptr (Complex Double) -> IO Int -multiply' :: (Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a +multiply' :: (Element a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a multiply' RowMajor a b = multiplyD a b multiply' ColumnMajor a b = trans $ multiplyD (trans b) (trans a) -- | matrix product -multiply :: (Field a) => Matrix a -> Matrix a -> Matrix a +multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a multiply = multiplyD ---------------------------------------------------------------------- @@ -287,7 +287,7 @@ subMatrixC (r0,c0) (rt,ct) x = reshape (2*cols x) . asReal . cdat $ x -- | Extracts a submatrix from a matrix. -subMatrix :: Field a +subMatrix :: Element a => (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix -> Matrix a -- ^ input matrix @@ -313,7 +313,7 @@ diagC = diagAux c_diagC "diagC" foreign import ccall "auxi.h diagC" c_diagC :: TCVCM -- | creates a square matrix with the given diagonal -diag :: Field a => Vector a -> Matrix a +diag :: Element a => Vector a -> Matrix a diag = diagD ------------------------------------------------------------------------ @@ -340,7 +340,7 @@ foreign import ccall safe "auxi.h constantC" @> constant 2 7 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ -} -constant :: Field a => a -> Int -> Vector a +constant :: Element a => a -> Int -> Vector a constant = constantD -------------------------------------------------------------------------- diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index a705975..e96500f 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs @@ -14,7 +14,7 @@ ----------------------------------------------------------------------------- module Data.Packed.Matrix ( - Field, + Element, Matrix,rows,cols, (><), trans, @@ -41,13 +41,13 @@ import Data.List(transpose,intersperse) import Data.Array -- | creates a matrix from a vertical list of matrices -joinVert :: Field t => [Matrix t] -> Matrix t +joinVert :: Element t => [Matrix t] -> Matrix t joinVert ms = case common cols ms of Nothing -> error "joinVert on matrices with different number of columns" Just c -> reshape c $ join (map cdat ms) -- | creates a matrix from a horizontal list of matrices -joinHoriz :: Field t => [Matrix t] -> Matrix t +joinHoriz :: Element t => [Matrix t] -> Matrix t joinHoriz ms = trans. joinVert . map trans $ ms {- | Creates a matrix from blocks given as a list of lists of matrices: @@ -63,15 +63,15 @@ joinHoriz ms = trans. joinVert . map trans $ ms , -1.0, -1.0, -1.0, -1.0, 0.0, 7.0, 0.0 , -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 2.0 ]@ -} -fromBlocks :: Field t => [[Matrix t]] -> Matrix t +fromBlocks :: Element t => [[Matrix t]] -> Matrix t fromBlocks = joinVert . map joinHoriz -- | Reverse rows -flipud :: Field t => Matrix t -> Matrix t +flipud :: Element t => Matrix t -> Matrix t flipud m = fromRows . reverse . toRows $ m -- | Reverse columns -fliprl :: Field t => Matrix t -> Matrix t +fliprl :: Element t => Matrix t -> Matrix t fliprl m = fromColumns . reverse . toColumns $ m ------------------------------------------------------------ @@ -84,7 +84,7 @@ fliprl m = fromColumns . reverse . toColumns $ m , 0.0, 5.0, 0.0, 0.0 , 0.0, 0.0, 5.0, 0.0 ]@ -} -diagRect :: (Field t, Num t) => Vector t -> Int -> Int -> Matrix t +diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t diagRect s r c | dim s < min r c = error "diagRect" | r == c = diag s @@ -93,11 +93,11 @@ diagRect s r c where zeros (r,c) = reshape c $ constantD 0 (r*c) -- | extracts the diagonal from a rectangular matrix -takeDiag :: (Field t) => Matrix t -> Vector t +takeDiag :: (Element t) => Matrix t -> Vector t takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] -- | creates the identity matrix of given dimension -ident :: Field a => Int -> Matrix a +ident :: Element a => Int -> Matrix a ident n = diag (constant 1 n) ------------------------------------------------------------ @@ -112,7 +112,7 @@ ident n = diag (constant 1 n) This is the format produced by the instances of Show (Matrix a), which can also be used for input. -} -(><) :: (Field a) => Int -> Int -> [a] -> Matrix a +(><) :: (Element a) => Int -> Int -> [a] -> Matrix a r >< c = f where f l | dim v == r*c = matrixFromVector RowMajor c v | otherwise = error $ "inconsistent list size = " @@ -122,16 +122,16 @@ r >< c = f where ---------------------------------------------------------------- -- | Creates a matrix with the first n rows of another matrix -takeRows :: Field t => Int -> Matrix t -> Matrix t +takeRows :: Element t => Int -> Matrix t -> Matrix t takeRows n mat = subMatrix (0,0) (n, cols mat) mat -- | Creates a copy of a matrix without the first n rows -dropRows :: Field t => Int -> Matrix t -> Matrix t +dropRows :: Element t => Int -> Matrix t -> Matrix t dropRows n mat = subMatrix (n,0) (rows mat - n, cols mat) mat -- |Creates a matrix with the first n columns of another matrix -takeColumns :: Field t => Int -> Matrix t -> Matrix t +takeColumns :: Element t => Int -> Matrix t -> Matrix t takeColumns n mat = subMatrix (0,0) (rows mat, n) mat -- | Creates a copy of a matrix without the first n columns -dropColumns :: Field t => Int -> Matrix t -> Matrix t +dropColumns :: Element t => Int -> Matrix t -> Matrix t dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat ---------------------------------------------------------------- @@ -141,7 +141,7 @@ dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat @\> flatten ('ident' 3) 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ -} -flatten :: Field t => Matrix t -> Vector t +flatten :: Element t => Matrix t -> Vector t flatten = cdat {- | Creates a 'Matrix' from a list of lists (considered as rows). @@ -152,20 +152,20 @@ flatten = cdat , 3.0, 4.0 , 5.0, 6.0 ]@ -} -fromLists :: Field t => [[t]] -> Matrix t +fromLists :: Element t => [[t]] -> Matrix t fromLists = fromRows . map fromList -- | creates a 1-row matrix from a vector -asRow :: Field a => Vector a -> Matrix a +asRow :: Element a => Vector a -> Matrix a asRow v = reshape (dim v) v -- | creates a 1-column matrix from a vector -asColumn :: Field a => Vector a -> Matrix a +asColumn :: Element a => Vector a -> Matrix a asColumn v = reshape 1 v ----------------------------------------------------- -fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e +fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e fromArray2D m = (r> String -> (t -> String) -> Matrix t -> String +format :: (Element t) => String -> (t -> String) -> Matrix t -> String format sep f m = dsp' sep . map (map f) . toLists $ m disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m @@ -217,7 +217,7 @@ readMatrix :: String -> Matrix Double readMatrix = fromLists . map (map read). map words . filter (not.null) . lines -- | rearranges the rows of a matrix according to the order given in a list of integers. -extractRows :: Field t => [Int] -> Matrix t -> Matrix t +extractRows :: Element t => [Int] -> Matrix t -> Matrix t extractRows l m = fromRows $ extract (toRows $ m) l where extract l is = [l!!i |i<-is] @@ -231,5 +231,5 @@ extractRows l m = fromRows $ extract (toRows $ m) l , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@ -} -repmat :: (Field t) => Matrix t -> Int -> Int -> Matrix t +repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t repmat m r c = fromBlocks $ partit c $ replicate (r*c) m \ No newline at end of file diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs deleted file mode 100644 index 254a957..0000000 --- a/lib/GSLHaskell.hs +++ /dev/null @@ -1,327 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} ------------------------------------------------------------------------------ -{- | -Module : GSLHaskell -Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : uses -fffi and -fglasgow-exts - -Old GSLHaskell interface. - --} ------------------------------------------------------------------------------ - -module GSLHaskell( - module Data.Packed.Vector, - module Data.Packed.Matrix, - module LinearAlgebra.Algorithms, - module LinearAlgebra.Linear, - module LAPACK, - module GSL.Integration, - module GSL.Differentiation, - module GSL.Fourier, - module GSL.Polynomials, - module GSL.Minimization, - module GSL.Matrix, - module GSL.Special, - module Graphics.Plot, - module Complex, - Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->), - fromArray2D, GSLHaskell.pnorm, -) where - -import GSL.Integration -import GSL.Differentiation -import GSL.Fourier -import GSL.Polynomials -import GSL.Minimization -import Graphics.Plot -import Complex -import GSL.Special(setErrorHandlerOff, - erf, - erf_Z, - bessel_J0_e, - exp_e10_e, - gamma) -import Data.Packed.Vector -import Data.Packed.Matrix -import Data.Packed.Matrix hiding ((><)) -import GSL.Vector -import qualified LinearAlgebra.Algorithms -import LAPACK -import GSL.Matrix -import LinearAlgebra.Algorithms hiding (pnorm) -import LinearAlgebra.Linear hiding (Mul,(<>)) -import Data.Packed.Internal.Matrix(multiply) -import Complex -import Numeric(showGFloat) -import Data.List(transpose,intersperse) -import Foreign(Storable) -import Data.Array -import LinearAlgebra.Instances - - -class Mul a b c | a b -> c where - infixl 7 <> -{- | An overloaded operator for matrix products, matrix-vector and vector-matrix products, dot products and scaling of vectors and matrices. Type consistency is statically checked. Alternatively, you can use the specific functions described below, but using this operator you can automatically combine real and complex objects. - -@v = 'fromList' [1,2,3] :: Vector Double -cv = 'fromList' [1+'i',2] -m = 'fromLists' [[1,2,3], - [4,5,7]] :: Matrix Double -cm = 'fromLists' [[ 1, 2], - [3+'i',7*'i'], - [ 'i', 1]] -\ -\> m \<\> v -14. 35. -\ -\> cv \<\> m -9.+1.i 12.+2.i 17.+3.i -\ -\> m \<\> cm - 7.+5.i 5.+14.i -19.+12.i 15.+35.i -\ -\> v \<\> 'i' -1.i 2.i 3.i -\ -\> v \<\> v -14.0 -\ -\> cv \<\> cv -4.0 :+ 2.0@ - --} - (<>) :: a -> b -> c - - -instance Mul Double Double Double where - (<>) = (*) - -instance Mul Double (Complex Double) (Complex Double) where - a <> b = (a:+0) * b - -instance Mul (Complex Double) Double (Complex Double) where - a <> b = a * (b:+0) - -instance Mul (Complex Double) (Complex Double) (Complex Double) where - (<>) = (*) - ---------------------------------- matrix matrix - -instance Mul (Matrix Double) (Matrix Double) (Matrix Double) where - (<>) = multiply - -instance Mul (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where - (<>) = multiply - -instance Mul (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where - c <> r = c <> liftMatrix comp r - -instance Mul (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where - r <> c = liftMatrix comp r <> c - ---------------------------------- (Matrix Double) (Vector Double) - -instance Mul (Matrix Double) (Vector Double) (Vector Double) where - (<>) = mXv - -instance Mul (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where - (<>) = mXv - -instance Mul (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where - m <> v = m <> comp v - -instance Mul (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where - m <> v = liftMatrix comp m <> v - ---------------------------------- (Vector Double) (Matrix Double) - -instance Mul (Vector Double) (Matrix Double) (Vector Double) where - (<>) = vXm - -instance Mul (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where - (<>) = vXm - -instance Mul (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where - v <> m = v <> liftMatrix comp m - -instance Mul (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where - v <> m = comp v <> m - ---------------------------------- dot product - -instance Mul (Vector Double) (Vector Double) Double where - (<>) = dot - -instance Mul (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where - (<>) = dot - -instance Mul (Vector Double) (Vector (Complex Double)) (Complex Double) where - a <> b = comp a <> b - -instance Mul (Vector (Complex Double)) (Vector Double) (Complex Double) where - (<>) = flip (<>) - ---------------------------------- scaling vectors - -instance Mul Double (Vector Double) (Vector Double) where - (<>) = scale - -instance Mul (Vector Double) Double (Vector Double) where - (<>) = flip (<>) - -instance Mul (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where - (<>) = scale - -instance Mul (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where - (<>) = flip (<>) - -instance Mul Double (Vector (Complex Double)) (Vector (Complex Double)) where - a <> v = (a:+0) <> v - -instance Mul (Vector (Complex Double)) Double (Vector (Complex Double)) where - (<>) = flip (<>) - -instance Mul (Complex Double) (Vector Double) (Vector (Complex Double)) where - a <> v = a <> comp v - -instance Mul (Vector Double) (Complex Double) (Vector (Complex Double)) where - (<>) = flip (<>) - ---------------------------------- scaling matrices - -instance Mul Double (Matrix Double) (Matrix Double) where - (<>) a = liftMatrix (a <>) - -instance Mul (Matrix Double) Double (Matrix Double) where - (<>) = flip (<>) - -instance Mul (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where - (<>) a = liftMatrix (a <>) - -instance Mul (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where - (<>) = flip (<>) - -instance Mul Double (Matrix (Complex Double)) (Matrix (Complex Double)) where - a <> m = (a:+0) <> m - -instance Mul (Matrix (Complex Double)) Double (Matrix (Complex Double)) where - (<>) = flip (<>) - -instance Mul (Complex Double) (Matrix Double) (Matrix (Complex Double)) where - a <> m = a <> liftMatrix comp m - -instance Mul (Matrix Double) (Complex Double) (Matrix (Complex Double)) where - (<>) = flip (<>) - ------------------------------------------------------------------------------------ - -size :: Vector a -> Int -size = dim - -gmap :: (Storable a, Storable b) => (a->b) -> Vector a -> Vector b -gmap f v = liftVector f v - --- shows a Double with n digits after the decimal point -shf :: (RealFloat a) => Int -> a -> String -shf dec n | abs n < 1e-10 = "0." - | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"." - | otherwise = showGFloat (Just dec) n "" --- shows a Complex Double as a pair, with n digits after the decimal point -shfc n z@ (a:+b) - | magnitude z <1e-10 = "0." - | abs b < 1e-10 = shf n a - | abs a < 1e-10 = shf n b ++"i" - | b > 0 = shf n a ++"+"++shf n b ++"i" - | otherwise = shf n a ++shf n b ++"i" - -dsp :: String -> [[String]] -> String -dsp sep as = unlines . map unwords' $ transpose mtp where - mt = transpose as - longs = map (maximum . map length) mt - mtp = zipWith (\a b -> map (pad a) b) longs mt - pad n str = replicate (n - length str) ' ' ++ str - unwords' = concat . intersperse sep - -format :: (Field t) => String -> (t -> String) -> Matrix t -> String -format sep f m = dsp sep . map (map f) . toLists $ m - -disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m - -dispR :: Int -> Matrix Double -> IO () -dispR d m = disp m (shf d) - -dispC :: Int -> Matrix (Complex Double) -> IO () -dispC d m = disp m (shfc d) - --- | creates a matrix from a table of numbers. -readMatrix :: String -> Matrix Double -readMatrix = fromLists . map (map read). map words . filter (not.null) . lines - -------------------------------------------------------------- - -class Joinable a b c | a b -> c where - joinH :: a -> b -> c - joinV :: a -> b -> c - -instance Joinable (Matrix Double) (Vector Double) (Matrix Double) where - joinH m v = fromBlocks [[m,reshape 1 v]] - joinV m v = fromBlocks [[m],[reshape (size v) v]] - -instance Joinable (Vector Double) (Matrix Double) (Matrix Double) where - joinH v m = fromBlocks [[reshape 1 v,m]] - joinV v m = fromBlocks [[reshape (size v) v],[m]] - -instance Joinable (Matrix Double) (Matrix Double) (Matrix Double) where - joinH m1 m2 = fromBlocks [[m1,m2]] - joinV m1 m2 = fromBlocks [[m1],[m2]] - -instance Joinable (Matrix (Complex Double)) (Vector (Complex Double)) (Matrix (Complex Double)) where - joinH m v = fromBlocks [[m,reshape 1 v]] - joinV m v = fromBlocks [[m],[reshape (size v) v]] - -instance Joinable (Vector (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where - joinH v m = fromBlocks [[reshape 1 v,m]] - joinV v m = fromBlocks [[reshape (size v) v],[m]] - -instance Joinable (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where - joinH m1 m2 = fromBlocks [[m1,m2]] - joinV m1 m2 = fromBlocks [[m1],[m2]] - -infixl 3 <|>, <-> - -{- | Horizontal concatenation of matrices and vectors: - -@\> 'ident' 3 \<-\> i\<\>'ident' 3 \<|\> 'fromList' [1..6] - 1. 0. 0. 1. - 0. 1. 0. 2. - 0. 0. 1. 3. -1.i 0. 0. 4. - 0. 1.i 0. 5. - 0. 0. 1.i 6.@ --} -(<|>) :: (Joinable a b c) => a -> b -> c -a <|> b = joinH a b - --- | Vertical concatenation of matrices and vectors. -(<->) :: (Joinable a b c) => a -> b -> c -a <-> b = joinV a b - ----------------------------------------------------------- - -fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e -fromArray2D m = (r> t -> t1 -> Double -pnorm 0 = LinearAlgebra.Algorithms.pnorm Infinity -pnorm 1 = LinearAlgebra.Algorithms.pnorm PNorm1 -pnorm 2 = LinearAlgebra.Algorithms.pnorm PNorm2 \ No newline at end of file diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 794ef69..b7e208a 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -51,7 +51,7 @@ module Numeric.LinearAlgebra.Algorithms ( -- * Util haussholder, unpackQR, unpackHess, - GenMat(linearSolveSVD,lu,eigSH',cholSH) + Field(linearSolveSVD,lu,eigSH',cholSH) ) where @@ -65,7 +65,7 @@ import Numeric.LinearAlgebra.Linear import Data.List(foldl1') -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. -class (Normed (Matrix t), Linear Matrix t) => GenMat t where +class (Normed (Matrix t), Linear Matrix t) => Field t where -- | Singular value decomposition using lapack's dgesvd or zgesvd. svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) @@ -103,7 +103,7 @@ class (Normed (Matrix t), Linear Matrix t) => GenMat t where ctrans :: Matrix t -> Matrix t -instance GenMat Double where +instance Field Double where svd = svdR lu = GSL.luR linearSolve = linearSolveR @@ -116,7 +116,7 @@ instance GenMat Double where hess = unpackHess hessR schur = schurR -instance GenMat (Complex Double) where +instance Field (Complex Double) where svd = svdC lu = GSL.luC linearSolve = linearSolveC @@ -132,37 +132,37 @@ instance GenMat (Complex Double) where -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. -- -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ -eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) +eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) eigSH m | m `equal` ctrans m = eigSH' m | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. -- -- If @c = chol m@ then @m == c \<> ctrans c@. -chol :: GenMat t => Matrix t -> Matrix t +chol :: Field t => Matrix t -> Matrix t chol m | m `equal` ctrans m = cholSH m | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" square m = rows m == cols m -det :: GenMat t => Matrix t -> t +det :: Field t => Matrix t -> t det m | square m = s * (product $ toList $ takeDiag $ u) | otherwise = error "det of nonsquare matrix" where (_,u,_,s) = lu m -- | Inverse of a square matrix using lapacks' dgesv and zgesv. -inv :: GenMat t => Matrix t -> Matrix t +inv :: Field t => Matrix t -> Matrix t inv m | square m = m `linearSolve` ident (rows m) | otherwise = error "inv of nonsquare matrix" -- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. -pinv :: GenMat t => Matrix t -> Matrix t +pinv :: Field t => Matrix t -> Matrix t pinv m = linearSolveSVD m (ident (rows m)) -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. -- -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. -full :: Field t +full :: Element 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 @@ -173,7 +173,7 @@ full svd m = (u, d ,v) where -- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. -- -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. -economy :: Field t +economy :: Element 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 @@ -198,15 +198,15 @@ i = 0:+1 -- matrix product -mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t +mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t mXm = multiply -- matrix - vector product -mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t +mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t mXv m v = flatten $ m `mXm` (asColumn v) -- vector - matrix product -vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t +vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t vXm v m = flatten $ (asRow v) `mXm` m @@ -264,7 +264,7 @@ instance Normed (Matrix (Complex Double)) where ----------------------------------------------------------------------- -- | The nullspace of a matrix from its SVD decomposition. -nullspacePrec :: GenMat t +nullspacePrec :: Field t => Double -- ^ relative tolerance in 'eps' units -> Matrix t -- ^ input matrix -> [Vector t] -- ^ list of unitary vectors spanning the nullspace @@ -276,7 +276,7 @@ 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 :: GenMat t => Matrix t -> Vector t +nullVector :: Field t => Matrix t -> Vector t nullVector = last . nullspacePrec 1 ------------------------------------------------------------------------ @@ -316,7 +316,7 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where -- many thanks, quickcheck! -haussholder :: (GenMat a) => a -> Vector a -> Matrix a +haussholder :: (Field a) => a -> Vector a -> Matrix a haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) where w = asColumn v @@ -328,7 +328,7 @@ zt 0 v = v zt k v = join [subVector 0 (dim v - k) v, constant 0 k] -unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) +unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) unpackQR (pq, tau) = (q,r) where cs = toColumns pq m = rows pq @@ -339,7 +339,7 @@ unpackQR (pq, tau) = (q,r) hs = zipWith haussholder (toList tau) vs q = foldl1' mXm hs -unpackHess :: (GenMat t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) +unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) unpackHess hf m | rows m == 1 = ((1><1)[1],m) | otherwise = (uH . hf) m @@ -357,13 +357,13 @@ uH (pq, tau) = (p,h) -------------------------------------------------------------------------- -- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. -rcond :: GenMat t => Matrix t -> Double +rcond :: Field t => Matrix t -> Double rcond m = last s / head s where (_,s',_) = svd m s = toList s' -- | Number of linearly independent rows or columns. -rank :: GenMat t => Matrix t -> Int +rank :: Field t => Matrix t -> Int rank m | pnorm PNorm1 m < eps = 0 | otherwise = dim s where (_,s,_) = economy svd m @@ -381,7 +381,7 @@ diagonalize m = if rank v == n else eig m -- | Generic matrix functions for diagonalizable matrices. -matFunc :: GenMat t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) +matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) matFunc f m = case diagonalize (complex m) of Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix" @@ -420,5 +420,5 @@ expGolub m = iterate msq f !! j {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, based on a scaled Pade approximation. -} -expm :: GenMat t => Matrix t -> Matrix t +expm :: Field t => Matrix t -> Matrix t expm = expGolub diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs index 4ee576f..3f295bf 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs @@ -29,7 +29,7 @@ import Foreign(Storable) ------------------------------------------------------------------ -instance (Show a, Field a) => (Show (Matrix a)) where +instance (Show a, Element a) => (Show (Matrix a)) where show m = (sizes++) . dsp . map (map show) . toLists $ m where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" @@ -51,7 +51,7 @@ adaptScalar f1 f2 f3 x y | dim y == 1 = f3 x (y@>0) | otherwise = f2 x y -liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t +liftMatrix2' :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2)) | otherwise = error "nonconformant matrices in liftMatrix2'" @@ -60,7 +60,7 @@ compat' m1 m2 = rows m1 == 1 && cols m1 == 1 || rows m2 == 1 && cols m2 == 1 || rows m1 == rows m2 && cols m1 == cols m2 -instance (Eq a, Field a) => Eq (Vector a) where +instance (Eq a, Element a) => Eq (Vector a) where a == b = dim a == dim b && toList a == toList b instance (Linear Vector a) => Num (Vector a) where @@ -71,7 +71,7 @@ instance (Linear Vector a) => Num (Vector a) where abs = liftVector abs fromInteger = fromList . return . fromInteger -instance (Eq a, Field a) => Eq (Matrix a) where +instance (Eq a, Element a) => Eq (Matrix a) where a == b = cols a == cols b && flatten a == flatten b instance (Linear Vector a) => Num (Matrix a) where diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs index fd076ec..4a9b309 100644 --- a/lib/Numeric/LinearAlgebra/Interface.hs +++ b/lib/Numeric/LinearAlgebra/Interface.hs @@ -29,7 +29,7 @@ import Numeric.LinearAlgebra.Algorithms class Mul a b c | a b -> c where infixl 7 <> -- | matrix product - (<>) :: Field t => a t -> b t -> c t + (<>) :: Element t => a t -> b t -> c t instance Mul Matrix Matrix Matrix where (<>) = multiply @@ -43,7 +43,7 @@ instance Mul Vector Matrix Vector where --------------------------------------------------- -- | @u \<.\> v = dot u v@ -(<.>) :: (Field t) => Vector t -> Vector t -> t +(<.>) :: (Element t) => Vector t -> Vector t -> t infixl 7 <.> (<.>) = dot @@ -62,15 +62,15 @@ 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 +(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a infixl 7 <\> m <\> v = flatten (linearSolveSVD m (reshape 1 v)) ------------------------------------------------ class Joinable a b where - joinH :: Field t => a t -> b t -> Matrix t - joinV :: Field t => a t -> b t -> Matrix t + joinH :: Element t => a t -> b t -> Matrix t + joinV :: Element t => a t -> b t -> Matrix t instance Joinable Matrix Matrix where joinH m1 m2 = fromBlocks [[m1,m2]] @@ -98,9 +98,9 @@ infixl 3 <-> , 0.0, 3.0, 0.0, 5.0 , 0.0, 0.0, 3.0, 6.0 ]@ -} -(<|>) :: (Field t, Joinable a b) => a t -> b t -> Matrix t +(<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t a <|> b = joinH a b -- | Vertical concatenation of matrices and vectors. -(<->) :: (Field t, Joinable a b) => a t -> b t -> Matrix t +(<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t a <-> b = joinV a b diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs index 94f6958..13d69ab 100644 --- a/lib/Numeric/LinearAlgebra/Linear.hs +++ b/lib/Numeric/LinearAlgebra/Linear.hs @@ -84,7 +84,7 @@ instance Linear Matrix (Complex Double) where -------------------------------------------------- -- | euclidean inner product -dot :: (Field t) => Vector t -> Vector t -> t +dot :: (Element t) => Vector t -> Vector t -> t dot u v = dat (multiply r c) `at` 0 where r = asRow u c = asColumn v @@ -98,7 +98,7 @@ dot u v = dat (multiply r c) `at` 0 , 10.0, 4.0, 6.0 , 15.0, 6.0, 9.0 ]@ -} -outer :: (Field t) => Vector t -> Vector t -> Matrix t +outer :: (Element t) => Vector t -> Vector t -> Matrix t outer u v = asColumn u `multiply` asRow v {- | Kronecker product of two matrices. @@ -123,7 +123,7 @@ m2=(4><3) , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ -} -kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t +kronecker :: (Element t) => Matrix t -> Matrix t -> Matrix t kronecker a b = fromBlocks . partit (cols a) . map (reshape (cols b)) -- cgit v1.2.3