From 42cc173794ee2cb52008ab642b9ab28fe6ad3eba Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 5 Jun 2015 16:59:29 +0200 Subject: move internal matrix --- packages/base/src/Data/Packed/Internal/Matrix.hs | 583 ----------------------- packages/base/src/Internal/Matrix.hs | 526 ++++++++++++++++++++ 2 files changed, 526 insertions(+), 583 deletions(-) delete mode 100644 packages/base/src/Data/Packed/Internal/Matrix.hs create mode 100644 packages/base/src/Internal/Matrix.hs (limited to 'packages') diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs deleted file mode 100644 index 66eef15..0000000 --- a/packages/base/src/Data/Packed/Internal/Matrix.hs +++ /dev/null @@ -1,583 +0,0 @@ -{-# LANGUAGE ForeignFunctionInterface #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE BangPatterns #-} - --- | --- Module : Data.Packed.Internal.Matrix --- Copyright : (c) Alberto Ruiz 2007 --- License : BSD3 --- Maintainer : Alberto Ruiz --- Stability : provisional --- --- Internal matrix representation --- - -module Data.Packed.Internal.Matrix( - Matrix(..), rows, cols, cdat, fdat, - MatrixOrder(..), orderOf, - createMatrix, mat, omat, - cmat, fmat, - toLists, flatten, reshape, - Element(..), - trans, - fromRows, toRows, fromColumns, toColumns, - matrixFromVector, - subMatrix, - liftMatrix, liftMatrix2, - (@@>), atM', - singleton, - emptyM, - size, shSize, conformVs, conformMs, conformVTo, conformMTo -) where - -import Data.Packed.Internal.Common -import Data.Packed.Internal.Signatures -import Data.Packed.Internal.Vector - -import Foreign.Marshal.Alloc(alloca, free) -import Foreign.Marshal.Array(newArray) -import Foreign.Ptr(Ptr, castPtr) -import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf) -import Data.Complex(Complex) -import Foreign.C.Types -import System.IO.Unsafe(unsafePerformIO) -import Control.DeepSeq - ------------------------------------------------------------------ - -{- Design considerations for the Matrix Type - ----------------------------------------- - -- we must easily handle both row major and column major order, - for bindings to LAPACK and GSL/C - -- we'd like to simplify redundant matrix transposes: - - Some of them arise from the order requirements of some functions - - some functions (matrix product) admit transposed arguments - -- maybe we don't really need this kind of simplification: - - more complex code - - some computational overhead - - only appreciable gain in code with a lot of redundant transpositions - and cheap matrix computations - -- we could carry both the matrix and its (lazily computed) transpose. - This may save some transpositions, but it is necessary to keep track of the - data which is actually computed to be used by functions like the matrix product - which admit both orders. - -- but if we need the transposed data and it is not in the structure, we must make - sure that we touch the same foreignptr that is used in the computation. - -- a reasonable solution is using two constructors for a matrix. Transposition just - "flips" the constructor. Actual data transposition is not done if followed by a - matrix product or another transpose. - --} - -data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) - -transOrder RowMajor = ColumnMajor -transOrder ColumnMajor = RowMajor -{- | Matrix representation suitable for BLAS\/LAPACK computations. - -The elements are stored in a continuous memory array. - --} - -data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int - , icols :: {-# UNPACK #-} !Int - , xdat :: {-# UNPACK #-} !(Vector t) - , order :: !MatrixOrder } --- RowMajor: preferred by C, fdat may require a transposition --- ColumnMajor: preferred by LAPACK, cdat may require a transposition - -cdat = xdat -fdat = xdat - -rows :: Matrix t -> Int -rows = irows - -cols :: Matrix t -> Int -cols = icols - -orderOf :: Matrix t -> MatrixOrder -orderOf = order - -stepRow :: Matrix t -> CInt -stepRow Matrix {icols = c, order = RowMajor } = fromIntegral c -stepRow _ = 1 - -stepCol :: Matrix t -> CInt -stepCol Matrix {irows = r, order = ColumnMajor } = fromIntegral r -stepCol _ = 1 - - --- | Matrix transpose. -trans :: Matrix t -> Matrix t -trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o} - -cmat :: (Element t) => Matrix t -> Matrix t -cmat m@Matrix{order = RowMajor} = m -cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor} - -fmat :: (Element t) => Matrix t -> Matrix t -fmat m@Matrix{order = ColumnMajor} = m -fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor} - --- C-Haskell matrix adapter --- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r - -mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b -mat a f = - unsafeWith (xdat a) $ \p -> do - let m g = do - g (fi (rows a)) (fi (cols a)) p - f m - -omat :: (Storable t) => Matrix t -> (((CInt -> CInt -> CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b -omat a f = - unsafeWith (xdat a) $ \p -> do - let m g = do - g (fi (rows a)) (fi (cols a)) (stepRow a) (stepCol a) p - f m - - -{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. - ->>> flatten (ident 3) -fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] - --} -flatten :: Element t => Matrix t -> Vector t -flatten = xdat . cmat - -{- -type Mt t s = Int -> Int -> Ptr t -> s - -infixr 6 ::> -type t ::> s = Mt t s --} - --- | the inverse of 'Data.Packed.Matrix.fromLists' -toLists :: (Element t) => Matrix t -> [[t]] -toLists m = splitEvery (cols m) . toList . flatten $ m - --- | Create a matrix from a list of vectors. --- All vectors must have the same dimension, --- or dimension 1, which is are automatically expanded. -fromRows :: Element t => [Vector t] -> Matrix t -fromRows [] = emptyM 0 0 -fromRows vs = case compatdim (map dim vs) of - Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs) - Just 0 -> emptyM r 0 - Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs - where - r = length vs - adapt c v - | c == 0 = fromList[] - | dim v == c = v - | otherwise = constantD (v@>0) c - --- | extracts the rows of a matrix as a list of vectors -toRows :: Element t => Matrix t -> [Vector t] -toRows m - | c == 0 = replicate r (fromList[]) - | otherwise = toRows' 0 - where - v = flatten m - r = rows m - c = cols m - toRows' k | k == r*c = [] - | otherwise = subVector k c v : toRows' (k+c) - --- | Creates a matrix from a list of vectors, as columns -fromColumns :: Element t => [Vector t] -> Matrix t -fromColumns m = trans . fromRows $ m - --- | Creates a list of vectors from the columns of a matrix -toColumns :: Element t => Matrix t -> [Vector t] -toColumns m = toRows . trans $ m - --- | Reads a matrix position. -(@@>) :: Storable t => Matrix t -> (Int,Int) -> t -infixl 9 @@> -m@Matrix {irows = r, icols = c} @@> (i,j) - | safe = if i<0 || i>=r || j<0 || j>=c - then error "matrix indexing out of range" - else atM' m i j - | otherwise = atM' m i j -{-# INLINE (@@>) #-} - --- Unsafe matrix access without range checking -atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j) -atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i) -{-# INLINE atM' #-} - ------------------------------------------------------------------- - -matrixFromVector o r c v - | r * c == dim v = m - | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m - where - m = Matrix { irows = r, icols = c, xdat = v, order = o } - --- allocates memory for a new matrix -createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a) -createMatrix ord r c = do - p <- createVector (r*c) - return (matrixFromVector ord r c p) - -{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@ -where r is the desired number of rows.) - ->>> reshape 4 (fromList [1..12]) -(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 ] - --} -reshape :: Storable t => Int -> Vector t -> Matrix t -reshape 0 v = matrixFromVector RowMajor 0 0 v -reshape c v = matrixFromVector RowMajor (dim v `div` c) c v - -singleton x = reshape 1 (fromList [x]) - --- | application of a vector function on the flattened matrix elements -liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b -liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d) - --- | application of a vector function on the flattened matrices elements -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 orderOf m1 of - RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2)) - ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2)) - - -compat :: Matrix a -> Matrix b -> Bool -compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 - ------------------------------------------------------------------- - -{- | Supported matrix elements. - - This class provides optimized internal - operations for selected element types. - It provides unoptimised defaults for any 'Storable' type, - so you can create instances simply as: - - >instance Element Foo --} -class (Storable a) => Element a where - subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position - -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix - -> Matrix a -> Matrix a - subMatrixD = subMatrix' - transdata :: Int -> Vector a -> Int -> Vector a - transdata = transdataP -- transdata' - constantD :: a -> Int -> Vector a - constantD = constantP -- constant' - extractR :: Matrix a -> CInt -> Vector CInt -> CInt -> Vector CInt -> Matrix a - sortI :: Ord a => Vector a -> Vector CInt - sortV :: Ord a => Vector a -> Vector a - compareV :: Ord a => Vector a -> Vector a -> Vector CInt - selectV :: Vector CInt -> Vector a -> Vector a -> Vector a -> Vector a - remapM :: Matrix CInt -> Matrix CInt -> Matrix a -> Matrix a - - -instance Element Float where - transdata = transdataAux ctransF - constantD = constantAux cconstantF - extractR = extractAux c_extractF - sortI = sortIdxF - sortV = sortValF - compareV = compareF - selectV = selectF - remapM = remapF - -instance Element Double where - transdata = transdataAux ctransR - constantD = constantAux cconstantR - extractR = extractAux c_extractD - sortI = sortIdxD - sortV = sortValD - compareV = compareD - selectV = selectD - remapM = remapD - - -instance Element (Complex Float) where - transdata = transdataAux ctransQ - constantD = constantAux cconstantQ - extractR = extractAux c_extractQ - sortI = undefined - sortV = undefined - compareV = undefined - selectV = selectQ - remapM = remapQ - - -instance Element (Complex Double) where - transdata = transdataAux ctransC - constantD = constantAux cconstantC - extractR = extractAux c_extractC - sortI = undefined - sortV = undefined - compareV = undefined - selectV = selectC - remapM = remapC - -instance Element (CInt) where - transdata = transdataAux ctransI - constantD = constantAux cconstantI - extractR = extractAux c_extractI - sortI = sortIdxI - sortV = sortValI - compareV = compareI - selectV = selectI - remapM = remapI - -------------------------------------------------------------------- - -transdataAux fun c1 d c2 = - if noneed - then d - else unsafePerformIO $ do - -- putStrLn "T" - v <- createVector (dim d) - unsafeWith d $ \pd -> - unsafeWith v $ \pv -> - fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux" - return v - where r1 = dim d `div` c1 - r2 = dim d `div` c2 - noneed = dim d == 0 || r1 == 1 || c1 == 1 - -transdataP :: Storable a => Int -> Vector a -> Int -> Vector a -transdataP c1 d c2 = - if noneed - then d - else unsafePerformIO $ do - v <- createVector (dim d) - unsafeWith d $ \pd -> - unsafeWith v $ \pv -> - ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP" - return v - where r1 = dim d `div` c1 - r2 = dim d `div` c2 - sz = sizeOf (d @> 0) - noneed = dim d == 0 || r1 == 1 || c1 == 1 - -foreign import ccall unsafe "transF" ctransF :: TFMFM -foreign import ccall unsafe "transR" ctransR :: TMM -foreign import ccall unsafe "transQ" ctransQ :: TQMQM -foreign import ccall unsafe "transC" ctransC :: TCMCM -foreign import ccall unsafe "transI" ctransI :: CM CInt (CM CInt (IO CInt)) -foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt - ----------------------------------------------------------------------- - -constantAux fun x n = unsafePerformIO $ do - v <- createVector n - px <- newArray [x] - app1 (fun px) vec v "constantAux" - free px - return v - -foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF - -foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV - -foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV - -foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV - -foreign import ccall unsafe "constantI" cconstantI :: Ptr CInt -> CV CInt (IO CInt) - -constantP :: Storable a => a -> Int -> Vector a -constantP a n = unsafePerformIO $ do - let sz = sizeOf a - v <- createVector n - unsafeWith v $ \p -> do - alloca $ \k -> do - poke k a - cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP" - return v -foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt - ----------------------------------------------------------------------- - --- | Extracts a submatrix from a matrix. -subMatrix :: Element a - => (Int,Int) -- ^ (r0,c0) starting position - -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix - -> Matrix a -- ^ input matrix - -> Matrix a -- ^ result -subMatrix (r0,c0) (rt,ct) m - | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) && - 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m - | otherwise = error $ "wrong subMatrix "++ - show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m) - -subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do - w <- createVector (rt*ct) - unsafeWith v $ \p -> - unsafeWith w $ \q -> do - let go (-1) _ = return () - go !i (-1) = go (i-1) (ct-1) - go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0) - pokeElemOff q (i*ct+j) x - go i (j-1) - go (rt-1) (ct-1) - return w - -subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor -subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m) - --------------------------------------------------------------------------- - -maxZ xs = if minimum xs == 0 then 0 else maximum xs - -conformMs ms = map (conformMTo (r,c)) ms - where - r = maxZ (map rows ms) - c = maxZ (map cols ms) - - -conformVs vs = map (conformVTo n) vs - where - n = maxZ (map dim vs) - -conformMTo (r,c) m - | size m == (r,c) = m - | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c)) - | size m == (r,1) = repCols c m - | size m == (1,c) = repRows r m - | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")" - -conformVTo n v - | dim v == n = v - | dim v == 1 = constantD (v@>0) n - | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n - -repRows n x = fromRows (replicate n (flatten x)) -repCols n x = fromColumns (replicate n (flatten x)) - -size m = (rows m, cols m) - -shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")" - -emptyM r c = matrixFromVector RowMajor r c (fromList[]) - ----------------------------------------------------------------------- - -instance (Storable t, NFData t) => NFData (Matrix t) - where - rnf m | d > 0 = rnf (v @> 0) - | otherwise = () - where - d = dim v - v = xdat m - ---------------------------------------------------------------- - -extractAux f m moder vr modec vc = unsafePerformIO $ do - let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr - nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc - r <- createMatrix RowMajor nr nc - app4 (f moder modec) vec vr vec vc omat m omat r "extractAux" - return r - -type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt)))) - -foreign import ccall unsafe "extractD" c_extractD :: Extr Double -foreign import ccall unsafe "extractF" c_extractF :: Extr Float -foreign import ccall unsafe "extractC" c_extractC :: Extr (Complex Double) -foreign import ccall unsafe "extractQ" c_extractQ :: Extr (Complex Float) -foreign import ccall unsafe "extractI" c_extractI :: Extr CInt - --------------------------------------------------------------------------------- - -sortG f v = unsafePerformIO $ do - r <- createVector (dim v) - app2 f vec v vec r "sortG" - return r - -sortIdxD = sortG c_sort_indexD -sortIdxF = sortG c_sort_indexF -sortIdxI = sortG c_sort_indexI - -sortValD = sortG c_sort_valD -sortValF = sortG c_sort_valF -sortValI = sortG c_sort_valI - -foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt)) -foreign import ccall unsafe "sort_indexF" c_sort_indexF :: CV Float (CV CInt (IO CInt)) -foreign import ccall unsafe "sort_indexI" c_sort_indexI :: CV CInt (CV CInt (IO CInt)) - -foreign import ccall unsafe "sort_valuesD" c_sort_valD :: CV Double (CV Double (IO CInt)) -foreign import ccall unsafe "sort_valuesF" c_sort_valF :: CV Float (CV Float (IO CInt)) -foreign import ccall unsafe "sort_valuesI" c_sort_valI :: CV CInt (CV CInt (IO CInt)) - --------------------------------------------------------------------------------- - -compareG f u v = unsafePerformIO $ do - r <- createVector (dim v) - app3 f vec u vec v vec r "compareG" - return r - -compareD = compareG c_compareD -compareF = compareG c_compareF -compareI = compareG c_compareI - -foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt))) -foreign import ccall unsafe "compareF" c_compareF :: CV Float (CV Float (CV CInt (IO CInt))) -foreign import ccall unsafe "compareI" c_compareI :: CV CInt (CV CInt (CV CInt (IO CInt))) - --------------------------------------------------------------------------------- - -selectG f c u v w = unsafePerformIO $ do - r <- createVector (dim v) - app5 f vec c vec u vec v vec w vec r "selectG" - return r - -selectD = selectG c_selectD -selectF = selectG c_selectF -selectI = selectG c_selectI -selectC = selectG c_selectC -selectQ = selectG c_selectQ - -type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt))))) - -foreign import ccall unsafe "chooseD" c_selectD :: Sel Double -foreign import ccall unsafe "chooseF" c_selectF :: Sel Float -foreign import ccall unsafe "chooseI" c_selectI :: Sel CInt -foreign import ccall unsafe "chooseC" c_selectC :: Sel (Complex Double) -foreign import ccall unsafe "chooseQ" c_selectQ :: Sel (Complex Float) - ---------------------------------------------------------------------------- - -remapG f i j m = unsafePerformIO $ do - r <- createMatrix RowMajor (rows i) (cols i) - app4 f omat i omat j omat m omat r "remapG" - return r - -remapD = remapG c_remapD -remapF = remapG c_remapF -remapI = remapG c_remapI -remapC = remapG c_remapC -remapQ = remapG c_remapQ - -type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt)))) - -foreign import ccall unsafe "remapD" c_remapD :: Rem Double -foreign import ccall unsafe "remapF" c_remapF :: Rem Float -foreign import ccall unsafe "remapI" c_remapI :: Rem CInt -foreign import ccall unsafe "remapC" c_remapC :: Rem (Complex Double) -foreign import ccall unsafe "remapQ" c_remapQ :: Rem (Complex Float) - - diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs new file mode 100644 index 0000000..44365d0 --- /dev/null +++ b/packages/base/src/Internal/Matrix.hs @@ -0,0 +1,526 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} + +-- | +-- Module : Internal.Matrix +-- Copyright : (c) Alberto Ruiz 2007-15 +-- License : BSD3 +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Internal matrix representation +-- + +module Internal.Matrix where + + +import Internal.Tools ( splitEvery, fi, compatdim, (//) ) +import Internal.Vector +import Internal.Devel +import Internal.Vectorized +import Data.Vector.Storable ( unsafeWith, fromList ) +import Foreign.Marshal.Alloc ( free ) +import Foreign.Ptr ( Ptr ) +import Foreign.Storable ( Storable ) +import Data.Complex ( Complex ) +import Foreign.C.Types ( CInt(..) ) +import Foreign.C.String ( CString, newCString ) +import System.IO.Unsafe ( unsafePerformIO ) +import Control.DeepSeq ( NFData(..) ) + + +----------------------------------------------------------------- + +{- Design considerations for the Matrix Type + ----------------------------------------- + +- we must easily handle both row major and column major order, + for bindings to LAPACK and GSL/C + +- we'd like to simplify redundant matrix transposes: + - Some of them arise from the order requirements of some functions + - some functions (matrix product) admit transposed arguments + +- maybe we don't really need this kind of simplification: + - more complex code + - some computational overhead + - only appreciable gain in code with a lot of redundant transpositions + and cheap matrix computations + +- we could carry both the matrix and its (lazily computed) transpose. + This may save some transpositions, but it is necessary to keep track of the + data which is actually computed to be used by functions like the matrix product + which admit both orders. + +- but if we need the transposed data and it is not in the structure, we must make + sure that we touch the same foreignptr that is used in the computation. + +- a reasonable solution is using two constructors for a matrix. Transposition just + "flips" the constructor. Actual data transposition is not done if followed by a + matrix product or another transpose. + +-} + +data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) + +transOrder RowMajor = ColumnMajor +transOrder ColumnMajor = RowMajor +{- | Matrix representation suitable for BLAS\/LAPACK computations. + +The elements are stored in a continuous memory array. + +-} + +data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int + , icols :: {-# UNPACK #-} !Int + , xdat :: {-# UNPACK #-} !(Vector t) + , order :: !MatrixOrder } +-- RowMajor: preferred by C, fdat may require a transposition +-- ColumnMajor: preferred by LAPACK, cdat may require a transposition + +--cdat = xdat +--fdat = xdat + +rows :: Matrix t -> Int +rows = irows + +cols :: Matrix t -> Int +cols = icols + +orderOf :: Matrix t -> MatrixOrder +orderOf = order + +stepRow :: Matrix t -> CInt +stepRow Matrix {icols = c, order = RowMajor } = fromIntegral c +stepRow _ = 1 + +stepCol :: Matrix t -> CInt +stepCol Matrix {irows = r, order = ColumnMajor } = fromIntegral r +stepCol _ = 1 + + +-- | Matrix transpose. +trans :: Matrix t -> Matrix t +trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o} + +cmat :: (Element t) => Matrix t -> Matrix t +cmat m@Matrix{order = RowMajor} = m +cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor} + +fmat :: (Element t) => Matrix t -> Matrix t +fmat m@Matrix{order = ColumnMajor} = m +fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor} + +-- C-Haskell matrix adapter +-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r + +mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b +mat a f = + unsafeWith (xdat a) $ \p -> do + let m g = do + g (fi (rows a)) (fi (cols a)) p + f m + +omat :: (Storable t) => Matrix t -> (((CInt -> CInt -> CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b +omat a f = + unsafeWith (xdat a) $ \p -> do + let m g = do + g (fi (rows a)) (fi (cols a)) (stepRow a) (stepCol a) p + f m + + +{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. + +>>> flatten (ident 3) +fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] + +-} +flatten :: Element t => Matrix t -> Vector t +flatten = xdat . cmat + +{- +type Mt t s = Int -> Int -> Ptr t -> s + +infixr 6 ::> +type t ::> s = Mt t s +-} + +-- | the inverse of 'Data.Packed.Matrix.fromLists' +toLists :: (Element t) => Matrix t -> [[t]] +toLists m = splitEvery (cols m) . toList . flatten $ m + +-- | Create a matrix from a list of vectors. +-- All vectors must have the same dimension, +-- or dimension 1, which is are automatically expanded. +fromRows :: Element t => [Vector t] -> Matrix t +fromRows [] = emptyM 0 0 +fromRows vs = case compatdim (map dim vs) of + Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs) + Just 0 -> emptyM r 0 + Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs + where + r = length vs + adapt c v + | c == 0 = fromList[] + | dim v == c = v + | otherwise = constantD (v@>0) c + +-- | extracts the rows of a matrix as a list of vectors +toRows :: Element t => Matrix t -> [Vector t] +toRows m + | c == 0 = replicate r (fromList[]) + | otherwise = toRows' 0 + where + v = flatten m + r = rows m + c = cols m + toRows' k | k == r*c = [] + | otherwise = subVector k c v : toRows' (k+c) + +-- | Creates a matrix from a list of vectors, as columns +fromColumns :: Element t => [Vector t] -> Matrix t +fromColumns m = trans . fromRows $ m + +-- | Creates a list of vectors from the columns of a matrix +toColumns :: Element t => Matrix t -> [Vector t] +toColumns m = toRows . trans $ m + +-- | Reads a matrix position. +(@@>) :: Storable t => Matrix t -> (Int,Int) -> t +infixl 9 @@> +m@Matrix {irows = r, icols = c} @@> (i,j) + | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range" + | otherwise = atM' m i j +{-# INLINE (@@>) #-} + +-- Unsafe matrix access without range checking +atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j) +atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i) +{-# INLINE atM' #-} + +------------------------------------------------------------------ + +matrixFromVector o r c v + | r * c == dim v = m + | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m + where + m = Matrix { irows = r, icols = c, xdat = v, order = o } + +-- allocates memory for a new matrix +createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a) +createMatrix ord r c = do + p <- createVector (r*c) + return (matrixFromVector ord r c p) + +{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@ +where r is the desired number of rows.) + +>>> reshape 4 (fromList [1..12]) +(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 ] + +-} +reshape :: Storable t => Int -> Vector t -> Matrix t +reshape 0 v = matrixFromVector RowMajor 0 0 v +reshape c v = matrixFromVector RowMajor (dim v `div` c) c v + +--singleton x = reshape 1 (fromList [x]) + +-- | application of a vector function on the flattened matrix elements +liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b +liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d) + +-- | application of a vector function on the flattened matrices elements +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 orderOf m1 of + RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2)) + ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2)) + + +compat :: Matrix a -> Matrix b -> Bool +compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 + +------------------------------------------------------------------ + +{- | Supported matrix elements. + + This class provides optimized internal + operations for selected element types. + It provides unoptimised defaults for any 'Storable' type, + so you can create instances simply as: + + >instance Element Foo +-} +class (Storable a) => Element a where + transdata :: Int -> Vector a -> Int -> Vector a + constantD :: a -> Int -> Vector a + extractR :: Matrix a -> CInt -> Vector CInt -> CInt -> Vector CInt -> Matrix a + sortI :: Ord a => Vector a -> Vector CInt + sortV :: Ord a => Vector a -> Vector a + compareV :: Ord a => Vector a -> Vector a -> Vector CInt + selectV :: Vector CInt -> Vector a -> Vector a -> Vector a -> Vector a + remapM :: Matrix CInt -> Matrix CInt -> Matrix a -> Matrix a + + +instance Element Float where + transdata = transdataAux ctransF + constantD = constantAux cconstantF + extractR = extractAux c_extractF + sortI = sortIdxF + sortV = sortValF + compareV = compareF + selectV = selectF + remapM = remapF + +instance Element Double where + transdata = transdataAux ctransR + constantD = constantAux cconstantR + extractR = extractAux c_extractD + sortI = sortIdxD + sortV = sortValD + compareV = compareD + selectV = selectD + remapM = remapD + + +instance Element (Complex Float) where + transdata = transdataAux ctransQ + constantD = constantAux cconstantQ + extractR = extractAux c_extractQ + sortI = undefined + sortV = undefined + compareV = undefined + selectV = selectQ + remapM = remapQ + + +instance Element (Complex Double) where + transdata = transdataAux ctransC + constantD = constantAux cconstantC + extractR = extractAux c_extractC + sortI = undefined + sortV = undefined + compareV = undefined + selectV = selectC + remapM = remapC + +instance Element (CInt) where + transdata = transdataAux ctransI + constantD = constantAux cconstantI + extractR = extractAux c_extractI + sortI = sortIdxI + sortV = sortValI + compareV = compareI + selectV = selectI + remapM = remapI + +------------------------------------------------------------------- + +transdataAux fun c1 d c2 = + if noneed + then d + else unsafePerformIO $ do + -- putStrLn "T" + v <- createVector (dim d) + unsafeWith d $ \pd -> + unsafeWith v $ \pv -> + fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux" + return v + where r1 = dim d `div` c1 + r2 = dim d `div` c2 + noneed = dim d == 0 || r1 == 1 || c1 == 1 + + +type TMM t = t ..> t ..> Ok + +foreign import ccall unsafe "transF" ctransF :: TMM Float +foreign import ccall unsafe "transR" ctransR :: TMM Double +foreign import ccall unsafe "transQ" ctransQ :: TMM (Complex Float) +foreign import ccall unsafe "transC" ctransC :: TMM (Complex Double) +foreign import ccall unsafe "transI" ctransI :: TMM CInt + +---------------------------------------------------------------------- + +-- | Extracts a submatrix from a matrix. +subMatrix :: Element a + => (Int,Int) -- ^ (r0,c0) starting position + -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix + -> Matrix a -- ^ input matrix + -> Matrix a -- ^ result +subMatrix (r0,c0) (rt,ct) m + | 0 <= r0 && 0 <= rt && r0+rt <= rows m && + 0 <= c0 && 0 <= ct && c0+ct <= cols m = extractR m 0 (idxs[r0,r0+rt-1]) 0 (idxs[c0,c0+ct-1]) + | otherwise = error $ "wrong subMatrix "++ + show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m) + +-------------------------------------------------------------------------- + +maxZ xs = if minimum xs == 0 then 0 else maximum xs + +conformMs ms = map (conformMTo (r,c)) ms + where + r = maxZ (map rows ms) + c = maxZ (map cols ms) + + +conformVs vs = map (conformVTo n) vs + where + n = maxZ (map dim vs) + +conformMTo (r,c) m + | size m == (r,c) = m + | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c)) + | size m == (r,1) = repCols c m + | size m == (1,c) = repRows r m + | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")" + +conformVTo n v + | dim v == n = v + | dim v == 1 = constantD (v@>0) n + | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n + +repRows n x = fromRows (replicate n (flatten x)) +repCols n x = fromColumns (replicate n (flatten x)) + +size m = (rows m, cols m) + +shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")" + +emptyM r c = matrixFromVector RowMajor r c (fromList[]) + +---------------------------------------------------------------------- + +instance (Storable t, NFData t) => NFData (Matrix t) + where + rnf m | d > 0 = rnf (v @> 0) + | otherwise = () + where + d = dim v + v = xdat m + +--------------------------------------------------------------- + +extractAux f m moder vr modec vc = unsafePerformIO $ do + let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr + nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc + r <- createMatrix RowMajor nr nc + app4 (f moder modec) vec vr vec vc omat m omat r "extractAux" + return r + +type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt)))) + +foreign import ccall unsafe "extractD" c_extractD :: Extr Double +foreign import ccall unsafe "extractF" c_extractF :: Extr Float +foreign import ccall unsafe "extractC" c_extractC :: Extr (Complex Double) +foreign import ccall unsafe "extractQ" c_extractQ :: Extr (Complex Float) +foreign import ccall unsafe "extractI" c_extractI :: Extr CInt + +-------------------------------------------------------------------------------- + +sortG f v = unsafePerformIO $ do + r <- createVector (dim v) + app2 f vec v vec r "sortG" + return r + +sortIdxD = sortG c_sort_indexD +sortIdxF = sortG c_sort_indexF +sortIdxI = sortG c_sort_indexI + +sortValD = sortG c_sort_valD +sortValF = sortG c_sort_valF +sortValI = sortG c_sort_valI + +foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt)) +foreign import ccall unsafe "sort_indexF" c_sort_indexF :: CV Float (CV CInt (IO CInt)) +foreign import ccall unsafe "sort_indexI" c_sort_indexI :: CV CInt (CV CInt (IO CInt)) + +foreign import ccall unsafe "sort_valuesD" c_sort_valD :: CV Double (CV Double (IO CInt)) +foreign import ccall unsafe "sort_valuesF" c_sort_valF :: CV Float (CV Float (IO CInt)) +foreign import ccall unsafe "sort_valuesI" c_sort_valI :: CV CInt (CV CInt (IO CInt)) + +-------------------------------------------------------------------------------- + +compareG f u v = unsafePerformIO $ do + r <- createVector (dim v) + app3 f vec u vec v vec r "compareG" + return r + +compareD = compareG c_compareD +compareF = compareG c_compareF +compareI = compareG c_compareI + +foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt))) +foreign import ccall unsafe "compareF" c_compareF :: CV Float (CV Float (CV CInt (IO CInt))) +foreign import ccall unsafe "compareI" c_compareI :: CV CInt (CV CInt (CV CInt (IO CInt))) + +-------------------------------------------------------------------------------- + +selectG f c u v w = unsafePerformIO $ do + r <- createVector (dim v) + app5 f vec c vec u vec v vec w vec r "selectG" + return r + +selectD = selectG c_selectD +selectF = selectG c_selectF +selectI = selectG c_selectI +selectC = selectG c_selectC +selectQ = selectG c_selectQ + +type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt))))) + +foreign import ccall unsafe "chooseD" c_selectD :: Sel Double +foreign import ccall unsafe "chooseF" c_selectF :: Sel Float +foreign import ccall unsafe "chooseI" c_selectI :: Sel CInt +foreign import ccall unsafe "chooseC" c_selectC :: Sel (Complex Double) +foreign import ccall unsafe "chooseQ" c_selectQ :: Sel (Complex Float) + +--------------------------------------------------------------------------- + +remapG f i j m = unsafePerformIO $ do + r <- createMatrix RowMajor (rows i) (cols i) + app4 f omat i omat j omat m omat r "remapG" + return r + +remapD = remapG c_remapD +remapF = remapG c_remapF +remapI = remapG c_remapI +remapC = remapG c_remapC +remapQ = remapG c_remapQ + +type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt)))) + +foreign import ccall unsafe "remapD" c_remapD :: Rem Double +foreign import ccall unsafe "remapF" c_remapF :: Rem Float +foreign import ccall unsafe "remapI" c_remapI :: Rem CInt +foreign import ccall unsafe "remapC" c_remapC :: Rem (Complex Double) +foreign import ccall unsafe "remapQ" c_remapQ :: Rem (Complex Float) + +-------------------------------------------------------------------------------- + +foreign import ccall unsafe "saveMatrix" c_saveMatrix + :: CString -> CString -> Double ..> Ok + +{- | save a matrix as a 2D ASCII table +-} +saveMatrix + :: FilePath + -> String -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.) + -> Matrix Double + -> IO () +saveMatrix name format m = do + cname <- newCString name + cformat <- newCString format + app1 (c_saveMatrix cname cformat) mat m "saveMatrix" + free cname + free cformat + return () + +-------------------------------------------------------------------------------- + -- cgit v1.2.3