From 5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 12:16:42 +0200 Subject: Data.Packed -> base (I) --- packages/hmatrix/src/Data/Packed.hs | 29 -- packages/hmatrix/src/Data/Packed/Development.hs | 32 -- packages/hmatrix/src/Data/Packed/Foreign.hs | 100 ---- packages/hmatrix/src/Data/Packed/Internal.hs | 26 - .../hmatrix/src/Data/Packed/Internal/Common.hs | 171 ------- .../hmatrix/src/Data/Packed/Internal/Matrix.hs | 491 ------------------- .../hmatrix/src/Data/Packed/Internal/Signatures.hs | 72 --- .../hmatrix/src/Data/Packed/Internal/Vector.hs | 521 --------------------- packages/hmatrix/src/Data/Packed/Matrix.hs | 490 ------------------- packages/hmatrix/src/Data/Packed/ST.hs | 179 ------- packages/hmatrix/src/Data/Packed/Vector.hs | 96 ---- 11 files changed, 2207 deletions(-) delete mode 100644 packages/hmatrix/src/Data/Packed.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Development.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Foreign.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Internal.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Internal/Common.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Internal/Matrix.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Internal/Signatures.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Internal/Vector.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Matrix.hs delete mode 100644 packages/hmatrix/src/Data/Packed/ST.hs delete mode 100644 packages/hmatrix/src/Data/Packed/Vector.hs (limited to 'packages/hmatrix') diff --git a/packages/hmatrix/src/Data/Packed.hs b/packages/hmatrix/src/Data/Packed.hs deleted file mode 100644 index 957aab8..0000000 --- a/packages/hmatrix/src/Data/Packed.hs +++ /dev/null @@ -1,29 +0,0 @@ ------------------------------------------------------------------------------ -{- | -Module : Data.Packed -Copyright : (c) Alberto Ruiz 2006-2010 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : uses ffi - -Types for dense 'Vector' and 'Matrix' of 'Storable' elements. - --} ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Packed ( - module Data.Packed.Vector, - module Data.Packed.Matrix, --- module Numeric.Conversion, --- module Data.Packed.Random, --- module Data.Complex -) where - -import Data.Packed.Vector -import Data.Packed.Matrix ---import Data.Packed.Random ---import Data.Complex ---import Numeric.Conversion diff --git a/packages/hmatrix/src/Data/Packed/Development.hs b/packages/hmatrix/src/Data/Packed/Development.hs deleted file mode 100644 index 471e560..0000000 --- a/packages/hmatrix/src/Data/Packed/Development.hs +++ /dev/null @@ -1,32 +0,0 @@ - ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Development --- Copyright : (c) Alberto Ruiz 2009 --- License : GPL --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- The library can be easily extended with additional foreign functions --- using the tools in this module. Illustrative usage examples can be found --- in the @examples\/devel@ folder included in the package. --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Packed.Development ( - createVector, createMatrix, - vec, mat, - app1, app2, app3, app4, - app5, app6, app7, app8, app9, app10, - MatrixOrder(..), orderOf, cmat, fmat, - matrixFromVector, - unsafeFromForeignPtr, - unsafeToForeignPtr, - check, (//), - at', atM' -) where - -import Data.Packed.Internal diff --git a/packages/hmatrix/src/Data/Packed/Foreign.hs b/packages/hmatrix/src/Data/Packed/Foreign.hs deleted file mode 100644 index 1ec3694..0000000 --- a/packages/hmatrix/src/Data/Packed/Foreign.hs +++ /dev/null @@ -1,100 +0,0 @@ -{-# LANGUAGE MagicHash, UnboxedTuples #-} --- | FFI and hmatrix helpers. --- --- Sample usage, to upload a perspective matrix to a shader. --- --- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) --- @ --- -{-# OPTIONS_HADDOCK hide #-} -module Data.Packed.Foreign - ( app - , appVector, appVectorLen - , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen - , unsafeMatrixToVector, unsafeMatrixToForeignPtr - ) where -import Data.Packed.Internal -import qualified Data.Vector.Storable as S -import Foreign (Ptr, ForeignPtr, Storable) -import Foreign.C.Types (CInt) -import GHC.Base (IO(..), realWorld#) - -{-# INLINE unsafeInlinePerformIO #-} --- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be --- unecessary calls to unsafePerformIO or its internals. -unsafeInlinePerformIO :: IO a -> a -unsafeInlinePerformIO (IO f) = case f realWorld# of - (# _, x #) -> x - -{-# INLINE app #-} --- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative. --- e.g. --- --- @ --- someFunction --- \`appMatrixLen\` m --- \`appVectorLen\` v --- \`app\` other --- \`app\` arguments --- \`app\` go here --- @ --- --- One could also write: --- --- @ --- (someFunction --- \`appMatrixLen\` m --- \`appVectorLen\` v) --- other --- arguments --- (go here) --- @ --- -app :: (a -> b) -> a -> b -app f = f - -{-# INLINE appVector #-} -appVector :: Storable a => (Ptr a -> b) -> Vector a -> b -appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f)) - -{-# INLINE appVectorLen #-} -appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b -appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x)))) - -{-# INLINE appMatrix #-} -appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b -appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f)) - -{-# INLINE appMatrixLen #-} -appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b -appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c)) - where - r = fromIntegral (rows x) - c = fromIntegral (cols x) - -{-# INLINE appMatrixRaw #-} -appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b -appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f)) - -{-# INLINE appMatrixRawLen #-} -appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b -appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c)) - where - r = fromIntegral (rows x) - c = fromIntegral (cols x) - -infixl 1 `app` -infixl 1 `appVector` -infixl 1 `appMatrix` -infixl 1 `appMatrixRaw` - -{-# INLINE unsafeMatrixToVector #-} --- | This will disregard the order of the matrix, and simply return it as-is. --- If the order of the matrix is RowMajor, this function is identical to 'flatten'. -unsafeMatrixToVector :: Matrix a -> Vector a -unsafeMatrixToVector = xdat - -{-# INLINE unsafeMatrixToForeignPtr #-} -unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int) -unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m) - diff --git a/packages/hmatrix/src/Data/Packed/Internal.hs b/packages/hmatrix/src/Data/Packed/Internal.hs deleted file mode 100644 index 537e51e..0000000 --- a/packages/hmatrix/src/Data/Packed/Internal.hs +++ /dev/null @@ -1,26 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- Reexports all internal modules --- ------------------------------------------------------------------------------ --- #hide - -module Data.Packed.Internal ( - module Data.Packed.Internal.Common, - module Data.Packed.Internal.Signatures, - module Data.Packed.Internal.Vector, - module Data.Packed.Internal.Matrix, -) where - -import Data.Packed.Internal.Common -import Data.Packed.Internal.Signatures -import Data.Packed.Internal.Vector -import Data.Packed.Internal.Matrix diff --git a/packages/hmatrix/src/Data/Packed/Internal/Common.hs b/packages/hmatrix/src/Data/Packed/Internal/Common.hs deleted file mode 100644 index edef3c2..0000000 --- a/packages/hmatrix/src/Data/Packed/Internal/Common.hs +++ /dev/null @@ -1,171 +0,0 @@ -{-# LANGUAGE CPP #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal.Common --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- Development utilities. --- ------------------------------------------------------------------------------ --- #hide - -module Data.Packed.Internal.Common( - Adapt, - app1, app2, app3, app4, - app5, app6, app7, app8, app9, app10, - (//), check, mbCatch, - splitEvery, common, compatdim, - fi, - table -) where - -import Foreign -import Control.Monad(when) -import Foreign.C.String(peekCString) -import Foreign.C.Types -import Foreign.Storable.Complex() -import Data.List(transpose,intersperse) -import Control.Exception as E - --- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@ -splitEvery :: Int -> [a] -> [[a]] -splitEvery _ [] = [] -splitEvery k l = take k l : splitEvery k (drop k l) - --- | obtains the common value of a property of a list -common :: (Eq a) => (b->a) -> [b] -> Maybe a -common f = commonval . map f where - commonval :: (Eq a) => [a] -> Maybe a - commonval [] = Nothing - commonval [a] = Just a - commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing - --- | common value with \"adaptable\" 1 -compatdim :: [Int] -> Maybe Int -compatdim [] = Nothing -compatdim [a] = Just a -compatdim (a:b:xs) - | a==b = compatdim (b:xs) - | a==1 = compatdim (b:xs) - | b==1 = compatdim (a:xs) - | otherwise = Nothing - --- | Formatting tool -table :: String -> [[String]] -> String -table 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 - --- | postfix function application (@flip ($)@) -(//) :: x -> (x -> y) -> y -infixl 0 // -(//) = flip ($) - --- | specialized fromIntegral -fi :: Int -> CInt -fi = fromIntegral - --- hmm.. -ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f -ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f -ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f -ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f -ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 f = w1 o1 $ ww5 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 . f -ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 f = w1 o1 $ ww6 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 . f -ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 f = w1 o1 $ ww7 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 . f -ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 f = w1 o1 $ ww8 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 . f -ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 f = w1 o1 $ ww9 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 . f - -type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() - -type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO() -type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2 -type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3 -type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4 -type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5 -type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 -type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 -type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 -type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 -type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 - -app1 :: f -> Adapt1 f t1 -app2 :: f -> Adapt2 f t1 r1 t2 -app3 :: f -> Adapt3 f t1 r1 t2 r2 t3 -app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4 -app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 -app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 -app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 -app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 -app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 -app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 - -app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s -app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s -app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $ - \a1 a2 a3 -> f // a1 // a2 // a3 // check s -app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $ - \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s -app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $ - \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s -app6 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 s = ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 $ - \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s -app7 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 s = ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 $ - \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s -app8 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 s = ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 $ - \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s -app9 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 s = ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 $ - \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s -app10 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 s = ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 $ - \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s - - - --- GSL error codes are <= 1024 --- | error codes for the auxiliary functions required by the wrappers -errorCode :: CInt -> String -errorCode 2000 = "bad size" -errorCode 2001 = "bad function code" -errorCode 2002 = "memory problem" -errorCode 2003 = "bad file" -errorCode 2004 = "singular" -errorCode 2005 = "didn't converge" -errorCode 2006 = "the input matrix is not positive definite" -errorCode 2007 = "not yet supported in this OS" -errorCode n = "code "++show n - - --- | clear the fpu -foreign import ccall unsafe "asm_finit" finit :: IO () - --- | check the error code -check :: String -> IO CInt -> IO () -check msg f = do -#if FINIT - finit -#endif - err <- f - when (err/=0) $ if err > 1024 - then (error (msg++": "++errorCode err)) -- our errors - else do -- GSL errors - ps <- gsl_strerror err - s <- peekCString ps - error (msg++": "++s) - return () - --- | description of GSL error codes -foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar) - --- | Error capture and conversion to Maybe -mbCatch :: IO x -> IO (Maybe x) -mbCatch act = E.catch (Just `fmap` act) f - where f :: SomeException -> IO (Maybe x) - f _ = return Nothing diff --git a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs deleted file mode 100644 index 9719fc0..0000000 --- a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs +++ /dev/null @@ -1,491 +0,0 @@ -{-# LANGUAGE ForeignFunctionInterface #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE BangPatterns #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal.Matrix --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- Internal matrix representation --- ------------------------------------------------------------------------------ --- #hide - -module Data.Packed.Internal.Matrix( - Matrix(..), rows, cols, cdat, fdat, - MatrixOrder(..), orderOf, - createMatrix, mat, - cmat, fmat, - toLists, flatten, reshape, - Element(..), - trans, - fromRows, toRows, fromColumns, toColumns, - matrixFromVector, - subMatrix, - liftMatrix, liftMatrix2, - (@@>), atM', - saveMatrix, - 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 Foreign.C.String(newCString) -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 GSL and 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 - - --- | 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 - -{- | 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' - - -instance Element Float where - transdata = transdataAux ctransF - constantD = constantAux cconstantF - -instance Element Double where - transdata = transdataAux ctransR - constantD = constantAux cconstantR - -instance Element (Complex Float) where - transdata = transdataAux ctransQ - constantD = constantAux cconstantQ - -instance Element (Complex Double) where - transdata = transdataAux ctransC - constantD = constantAux cconstantC - -------------------------------------------------------------------- - -transdata' :: Storable a => Int -> Vector a -> Int -> Vector a -transdata' c1 v c2 = - if noneed - then v - else unsafePerformIO $ do - w <- createVector (r2*c2) - unsafeWith v $ \p -> - unsafeWith w $ \q -> do - let go (-1) _ = return () - go !i (-1) = go (i-1) (c1-1) - go !i !j = do x <- peekElemOff p (i*c1+j) - pokeElemOff q (j*c2+i) x - go i (j-1) - go (r1-1) (c1-1) - return w - where r1 = dim v `div` c1 - r2 = dim v `div` c2 - noneed = dim v == 0 || r1 == 1 || c1 == 1 - --- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int -> Vector Double #-} --- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-} - --- I don't know how to specialize... --- The above pragmas only seem to work on top level defs --- Fortunately everything seems to work using the above class - --- C versions, still a little faster: - -transdataAux fun c1 d c2 = - if noneed - then d - else unsafePerformIO $ do - 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 "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt - ----------------------------------------------------------------------- - -constant' v n = unsafePerformIO $ do - w <- createVector n - unsafeWith w $ \p -> do - let go (-1) = return () - go !k = pokeElemOff p k v >> go (k-1) - go (n-1) - return w - --- C versions - -constantAux fun x n = unsafePerformIO $ do - v <- createVector n - px <- newArray [x] - app1 (fun px) vec v "constantAux" - free px - return v - -constantF :: Float -> Int -> Vector Float -constantF = constantAux cconstantF -foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF - -constantR :: Double -> Int -> Vector Double -constantR = constantAux cconstantR -foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV - -constantQ :: Complex Float -> Int -> Vector (Complex Float) -constantQ = constantAux cconstantQ -foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV - -constantC :: Complex Double -> Int -> Vector (Complex Double) -constantC = constantAux cconstantC -foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV - -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) - --------------------------------------------------------------------------- - --- | 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 - ----------------------------------------------------------------------- - -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 - diff --git a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs deleted file mode 100644 index 2835720..0000000 --- a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs +++ /dev/null @@ -1,72 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal.Signatures --- Copyright : (c) Alberto Ruiz 2009 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- Signatures of the C functions. --- ------------------------------------------------------------------------------ - -module Data.Packed.Internal.Signatures where - -import Foreign.Ptr(Ptr) -import Data.Complex(Complex) -import Foreign.C.Types(CInt) - -type PF = Ptr Float -- -type PD = Ptr Double -- -type PQ = Ptr (Complex Float) -- -type PC = Ptr (Complex Double) -- -type TF = CInt -> PF -> IO CInt -- -type TFF = CInt -> PF -> TF -- -type TFV = CInt -> PF -> TV -- -type TVF = CInt -> PD -> TF -- -type TFFF = CInt -> PF -> TFF -- -type TV = CInt -> PD -> IO CInt -- -type TVV = CInt -> PD -> TV -- -type TVVV = CInt -> PD -> TVV -- -type TFM = CInt -> CInt -> PF -> IO CInt -- -type TFMFM = CInt -> CInt -> PF -> TFM -- -type TFMFMFM = CInt -> CInt -> PF -> TFMFM -- -type TM = CInt -> CInt -> PD -> IO CInt -- -type TMM = CInt -> CInt -> PD -> TM -- -type TVMM = CInt -> PD -> TMM -- -type TMVMM = CInt -> CInt -> PD -> TVMM -- -type TMMM = CInt -> CInt -> PD -> TMM -- -type TVM = CInt -> PD -> TM -- -type TVVM = CInt -> PD -> TVM -- -type TMV = CInt -> CInt -> PD -> TV -- -type TMMV = CInt -> CInt -> PD -> TMV -- -type TMVM = CInt -> CInt -> PD -> TVM -- -type TMMVM = CInt -> CInt -> PD -> TMVM -- -type TCM = CInt -> CInt -> PC -> IO CInt -- -type TCVCM = CInt -> PC -> TCM -- -type TCMCVCM = CInt -> CInt -> PC -> TCVCM -- -type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM -- -type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM -- -type TCMCM = CInt -> CInt -> PC -> TCM -- -type TVCM = CInt -> PD -> TCM -- -type TCMVCM = CInt -> CInt -> PC -> TVCM -- -type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM -- -type TCMCMCM = CInt -> CInt -> PC -> TCMCM -- -type TCV = CInt -> PC -> IO CInt -- -type TCVCV = CInt -> PC -> TCV -- -type TCVCVCV = CInt -> PC -> TCVCV -- -type TCVV = CInt -> PC -> TV -- -type TQV = CInt -> PQ -> IO CInt -- -type TQVQV = CInt -> PQ -> TQV -- -type TQVQVQV = CInt -> PQ -> TQVQV -- -type TQVF = CInt -> PQ -> TF -- -type TQM = CInt -> CInt -> PQ -> IO CInt -- -type TQMQM = CInt -> CInt -> PQ -> TQM -- -type TQMQMQM = CInt -> CInt -> PQ -> TQMQM -- -type TCMCV = CInt -> CInt -> PC -> TCV -- -type TVCV = CInt -> PD -> TCV -- -type TCVM = CInt -> PC -> TM -- -type TMCVM = CInt -> CInt -> PD -> TCVM -- -type TMMCVM = CInt -> CInt -> PD -> TMCVM -- diff --git a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs deleted file mode 100644 index 6d03438..0000000 --- a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs +++ /dev/null @@ -1,521 +0,0 @@ -{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal.Vector --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- Vector implementation --- ------------------------------------------------------------------------------ - -module Data.Packed.Internal.Vector ( - Vector, dim, - fromList, toList, (|>), - vjoin, (@>), safe, at, at', subVector, takesV, - mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith, - mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, - foldVector, foldVectorG, foldLoop, foldVectorWithIndex, - createVector, vec, - asComplex, asReal, float2DoubleV, double2FloatV, - stepF, stepD, condF, condD, - conjugateQ, conjugateC, - fwriteVector, freadVector, fprintfVector, fscanfVector, - cloneVector, - unsafeToForeignPtr, - unsafeFromForeignPtr, - unsafeWith -) where - -import Data.Packed.Internal.Common -import Data.Packed.Internal.Signatures -import Foreign.Marshal.Alloc(free) -import Foreign.Marshal.Array(peekArray, copyArray, advancePtr) -import Foreign.ForeignPtr(ForeignPtr, castForeignPtr) -import Foreign.Ptr(Ptr) -import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf) -import Foreign.C.String -import Foreign.C.Types -import Data.Complex -import Control.Monad(when) -import System.IO.Unsafe(unsafePerformIO) - -#if __GLASGOW_HASKELL__ >= 605 -import GHC.ForeignPtr (mallocPlainForeignPtrBytes) -#else -import Foreign.ForeignPtr (mallocForeignPtrBytes) -#endif - -import GHC.Base -#if __GLASGOW_HASKELL__ < 612 -import GHC.IOBase hiding (liftIO) -#endif - -import qualified Data.Vector.Storable as Vector -import Data.Vector.Storable(Vector, - fromList, - unsafeToForeignPtr, - unsafeFromForeignPtr, - unsafeWith) - - --- | Number of elements -dim :: (Storable t) => Vector t -> Int -dim = Vector.length - - --- C-Haskell vector adapter --- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r -vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b -vec x f = unsafeWith x $ \p -> do - let v g = do - g (fi $ dim x) p - f v -{-# INLINE vec #-} - - --- allocates memory for a new vector -createVector :: Storable a => Int -> IO (Vector a) -createVector n = do - when (n < 0) $ error ("trying to createVector of negative dim: "++show n) - fp <- doMalloc undefined - return $ unsafeFromForeignPtr fp 0 n - where - -- - -- Use the much cheaper Haskell heap allocated storage - -- for foreign pointer space we control - -- - doMalloc :: Storable b => b -> IO (ForeignPtr b) - doMalloc dummy = do -#if __GLASGOW_HASKELL__ >= 605 - mallocPlainForeignPtrBytes (n * sizeOf dummy) -#else - mallocForeignPtrBytes (n * sizeOf dummy) -#endif - -{- | creates a Vector from a list: - -@> fromList [2,3,5,7] -4 |> [2.0,3.0,5.0,7.0]@ - --} - -safeRead v = inlinePerformIO . unsafeWith v -{-# INLINE safeRead #-} - -inlinePerformIO :: IO a -> a -inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r -{-# INLINE inlinePerformIO #-} - -{- | extracts the Vector elements to a list - ->>> toList (linspace 5 (1,10)) -[1.0,3.25,5.5,7.75,10.0] - --} -toList :: Storable a => Vector a -> [a] -toList v = safeRead v $ peekArray (dim v) - -{- | Create a vector from a list of elements and explicit dimension. The input - list is explicitly truncated if it is too long, so it may safely - be used, for instance, with infinite lists. - ->>> 5 |> [1..] -fromList [1.0,2.0,3.0,4.0,5.0] - --} -(|>) :: (Storable a) => Int -> [a] -> Vector a -infixl 9 |> -n |> l = if length l' == n - then fromList l' - else error "list too short for |>" - where l' = take n l - - --- | access to Vector elements without range checking -at' :: Storable a => Vector a -> Int -> a -at' v n = safeRead v $ flip peekElemOff n -{-# INLINE at' #-} - --- --- turn off bounds checking with -funsafe at configure time. --- ghc will optimise away the salways true case at compile time. --- -#if defined(UNSAFE) -safe :: Bool -safe = False -#else -safe = True -#endif - --- | access to Vector elements with range checking. -at :: Storable a => Vector a -> Int -> a -at v n - | safe = if n >= 0 && n < dim v - then at' v n - else error "vector index out of range" - | otherwise = at' v n -{-# INLINE at #-} - -{- | takes a number of consecutive elements from a Vector - ->>> subVector 2 3 (fromList [1..10]) -fromList [3.0,4.0,5.0] - --} -subVector :: Storable t => Int -- ^ index of the starting element - -> Int -- ^ number of elements to extract - -> Vector t -- ^ source - -> Vector t -- ^ result -subVector = Vector.slice - - -{- | Reads a vector position: - ->>> fromList [0..9] @> 7 -7.0 - --} -(@>) :: Storable t => Vector t -> Int -> t -infixl 9 @> -(@>) = at - - -{- | concatenate a list of vectors - ->>> vjoin [fromList [1..5::Double], konst 1 3] -fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] - --} -vjoin :: Storable t => [Vector t] -> Vector t -vjoin [] = fromList [] -vjoin [v] = v -vjoin as = unsafePerformIO $ do - let tot = sum (map dim as) - r <- createVector tot - unsafeWith r $ \ptr -> - joiner as tot ptr - return r - where joiner [] _ _ = return () - joiner (v:cs) _ p = do - let n = dim v - unsafeWith v $ \pb -> copyArray p pb n - joiner cs 0 (advancePtr p n) - - -{- | Extract consecutive subvectors of the given sizes. - ->>> takesV [3,4] (linspace 10 (1,10::Double)) -[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]] - --} -takesV :: Storable t => [Int] -> Vector t -> [Vector t] -takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w) - | otherwise = go ms w - where go [] _ = [] - go (n:ns) v = subVector 0 n v - : go ns (subVector n (dim v - n) v) - ---------------------------------------------------------------- - --- | transforms a complex vector into a real vector with alternating real and imaginary parts -asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a -asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) - where (fp,i,n) = unsafeToForeignPtr v - --- | transforms a real vector into a complex vector with alternating real and imaginary parts -asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a) -asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2) - where (fp,i,n) = unsafeToForeignPtr v - ---------------------------------------------------------------- - -float2DoubleV :: Vector Float -> Vector Double -float2DoubleV v = unsafePerformIO $ do - r <- createVector (dim v) - app2 c_float2double vec v vec r "float2double" - return r - -double2FloatV :: Vector Double -> Vector Float -double2FloatV v = unsafePerformIO $ do - r <- createVector (dim v) - app2 c_double2float vec v vec r "double2float2" - return r - - -foreign import ccall unsafe "float2double" c_float2double:: TFV -foreign import ccall unsafe "double2float" c_double2float:: TVF - ---------------------------------------------------------------- - -stepF :: Vector Float -> Vector Float -stepF v = unsafePerformIO $ do - r <- createVector (dim v) - app2 c_stepF vec v vec r "stepF" - return r - -stepD :: Vector Double -> Vector Double -stepD v = unsafePerformIO $ do - r <- createVector (dim v) - app2 c_stepD vec v vec r "stepD" - return r - -foreign import ccall unsafe "stepF" c_stepF :: TFF -foreign import ccall unsafe "stepD" c_stepD :: TVV - ---------------------------------------------------------------- - -condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -condF x y l e g = unsafePerformIO $ do - r <- createVector (dim x) - app6 c_condF vec x vec y vec l vec e vec g vec r "condF" - return r - -condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -condD x y l e g = unsafePerformIO $ do - r <- createVector (dim x) - app6 c_condD vec x vec y vec l vec e vec g vec r "condD" - return r - -foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF -foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV - --------------------------------------------------------------------------------- - -conjugateAux fun x = unsafePerformIO $ do - v <- createVector (dim x) - app2 fun vec x vec v "conjugateAux" - return v - -conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) -conjugateQ = conjugateAux c_conjugateQ -foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV - -conjugateC :: Vector (Complex Double) -> Vector (Complex Double) -conjugateC = conjugateAux c_conjugateC -foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV - --------------------------------------------------------------------------------- - -cloneVector :: Storable t => Vector t -> IO (Vector t) -cloneVector v = do - let n = dim v - r <- createVector n - let f _ s _ d = copyArray d s n >> return 0 - app2 f vec v vec r "cloneVector" - return r - ------------------------------------------------------------------- - --- | map on Vectors -mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b -mapVector f v = unsafePerformIO $ do - w <- createVector (dim v) - unsafeWith v $ \p -> - unsafeWith w $ \q -> do - let go (-1) = return () - go !k = do x <- peekElemOff p k - pokeElemOff q k (f x) - go (k-1) - go (dim v -1) - return w -{-# INLINE mapVector #-} - --- | zipWith for Vectors -zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c -zipVectorWith f u v = unsafePerformIO $ do - let n = min (dim u) (dim v) - w <- createVector n - unsafeWith u $ \pu -> - unsafeWith v $ \pv -> - unsafeWith w $ \pw -> do - let go (-1) = return () - go !k = do x <- peekElemOff pu k - y <- peekElemOff pv k - pokeElemOff pw k (f x y) - go (k-1) - go (n -1) - return w -{-# INLINE zipVectorWith #-} - --- | unzipWith for Vectors -unzipVectorWith :: (Storable (a,b), Storable c, Storable d) - => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) -unzipVectorWith f u = unsafePerformIO $ do - let n = dim u - v <- createVector n - w <- createVector n - unsafeWith u $ \pu -> - unsafeWith v $ \pv -> - unsafeWith w $ \pw -> do - let go (-1) = return () - go !k = do z <- peekElemOff pu k - let (x,y) = f z - pokeElemOff pv k x - pokeElemOff pw k y - go (k-1) - go (n-1) - return (v,w) -{-# INLINE unzipVectorWith #-} - -foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b -foldVector f x v = unsafePerformIO $ - unsafeWith v $ \p -> do - let go (-1) s = return s - go !k !s = do y <- peekElemOff p k - go (k-1::Int) (f y s) - go (dim v -1) x -{-# INLINE foldVector #-} - --- the zero-indexed index is passed to the folding function -foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b -foldVectorWithIndex f x v = unsafePerformIO $ - unsafeWith v $ \p -> do - let go (-1) s = return s - go !k !s = do y <- peekElemOff p k - go (k-1::Int) (f k y s) - go (dim v -1) x -{-# INLINE foldVectorWithIndex #-} - -foldLoop f s0 d = go (d - 1) s0 - where - go 0 s = f (0::Int) s - go !j !s = go (j - 1) (f j s) - -foldVectorG f s0 v = foldLoop g s0 (dim v) - where g !k !s = f k (at' v) s - {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) -{-# INLINE foldVectorG #-} - -------------------------------------------------------------------- - --- | monadic map over Vectors --- the monad @m@ must be strict -mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b) -mapVectorM f v = do - w <- return $! unsafePerformIO $! createVector (dim v) - mapVectorM' w 0 (dim v -1) - return w - where mapVectorM' w' !k !t - | k == t = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - y <- f x - return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y - | otherwise = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - y <- f x - _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y - mapVectorM' w' (k+1) t -{-# INLINE mapVectorM #-} - --- | monadic map over Vectors -mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m () -mapVectorM_ f v = do - mapVectorM' 0 (dim v -1) - where mapVectorM' !k !t - | k == t = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - f x - | otherwise = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - _ <- f x - mapVectorM' (k+1) t -{-# INLINE mapVectorM_ #-} - --- | monadic map over Vectors with the zero-indexed index passed to the mapping function --- the monad @m@ must be strict -mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b) -mapVectorWithIndexM f v = do - w <- return $! unsafePerformIO $! createVector (dim v) - mapVectorM' w 0 (dim v -1) - return w - where mapVectorM' w' !k !t - | k == t = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - y <- f k x - return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y - | otherwise = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - y <- f k x - _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y - mapVectorM' w' (k+1) t -{-# INLINE mapVectorWithIndexM #-} - --- | monadic map over Vectors with the zero-indexed index passed to the mapping function -mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m () -mapVectorWithIndexM_ f v = do - mapVectorM' 0 (dim v -1) - where mapVectorM' !k !t - | k == t = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - f k x - | otherwise = do - x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k - _ <- f k x - mapVectorM' (k+1) t -{-# INLINE mapVectorWithIndexM_ #-} - - -mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b ---mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b]) -mapVectorWithIndex f v = unsafePerformIO $ do - w <- createVector (dim v) - unsafeWith v $ \p -> - unsafeWith w $ \q -> do - let go (-1) = return () - go !k = do x <- peekElemOff p k - pokeElemOff q k (f k x) - go (k-1) - go (dim v -1) - return w -{-# INLINE mapVectorWithIndex #-} - -------------------------------------------------------------------- - - --- | 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 - diff --git a/packages/hmatrix/src/Data/Packed/Matrix.hs b/packages/hmatrix/src/Data/Packed/Matrix.hs deleted file mode 100644 index d94d167..0000000 --- a/packages/hmatrix/src/Data/Packed/Matrix.hs +++ /dev/null @@ -1,490 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE CPP #-} - ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Matrix --- Copyright : (c) Alberto Ruiz 2007-10 --- License : GPL --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- --- A Matrix representation suitable for numerical computations using LAPACK and GSL. --- --- This module provides basic functions for manipulation of structure. - ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Packed.Matrix ( - Matrix, - Element, - rows,cols, - (><), - trans, - reshape, flatten, - fromLists, toLists, buildMatrix, - (@@>), - asRow, asColumn, - fromRows, toRows, fromColumns, toColumns, - fromBlocks, diagBlock, toBlocks, toBlocksEvery, - repmat, - flipud, fliprl, - subMatrix, takeRows, dropRows, takeColumns, dropColumns, - extractRows, extractColumns, - diagRect, takeDiag, - mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_, - liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D -) where - -import Data.Packed.Internal -import qualified Data.Packed.ST as ST -import Data.Array - -import Data.List(transpose,intersperse) -import Foreign.Storable(Storable) -import Control.Monad(liftM) - -------------------------------------------------------------------- - -#ifdef BINARY - -import Data.Binary -import Control.Monad(replicateM) - -instance (Binary a, Element a, Storable a) => Binary (Matrix a) where - put m = do - let r = rows m - let c = cols m - put r - put c - mapM_ (\i -> mapM_ (\j -> put $ m @@> (i,j)) [0..(c-1)]) [0..(r-1)] - get = do - r <- get - c <- get - xs <- replicateM r $ replicateM c get - return $ fromLists xs - -#endif - -------------------------------------------------------------------- - -instance (Show a, Element a) => (Show (Matrix a)) where - show m | rows m == 0 || cols m == 0 = sizes m ++" []" - show m = (sizes m++) . dsp . map (map show) . toLists $ m - -sizes m = "("++show (rows m)++"><"++show (cols m)++")\n" - -dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . 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 ", " - ------------------------------------------------------------------- - -instance (Element a, Read a) => Read (Matrix a) where - readsPrec _ s = [((rs>' $ dims - - -breakAt c l = (a++[c],tail b) where - (a,b) = break (==c) l - ------------------------------------------------------------------- - --- | creates a matrix from a vertical list of matrices -joinVert :: Element t => [Matrix t] -> Matrix t -joinVert [] = emptyM 0 0 -joinVert ms = case common cols ms of - Nothing -> error "(impossible) joinVert on matrices with different number of columns" - Just c -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms) - --- | creates a matrix from a horizontal list of matrices -joinHoriz :: Element t => [Matrix t] -> Matrix t -joinHoriz ms = trans. joinVert . map trans $ ms - -{- | Create a matrix from blocks given as a list of lists of matrices. - -Single row-column components are automatically expanded to match the -corresponding common row and column: - -@ -disp = putStr . dispf 2 -@ - ->>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]] -8x10 -1 0 0 0 0 7 7 7 10 20 -0 1 0 0 0 7 7 7 10 20 -0 0 1 0 0 7 7 7 10 20 -0 0 0 1 0 7 7 7 10 20 -0 0 0 0 1 7 7 7 10 20 -3 3 3 3 3 1 0 0 0 0 -3 3 3 3 3 0 2 0 0 0 -3 3 3 3 3 0 0 3 0 0 - --} -fromBlocks :: Element t => [[Matrix t]] -> Matrix t -fromBlocks = fromBlocksRaw . adaptBlocks - -fromBlocksRaw mms = joinVert . map joinHoriz $ mms - -adaptBlocks ms = ms' where - bc = case common length ms of - Just c -> c - Nothing -> error "fromBlocks requires rectangular [[Matrix]]" - rs = map (compatdim . map rows) ms - cs = map (compatdim . map cols) (transpose ms) - szs = sequence [rs,cs] - ms' = splitEvery bc $ zipWith g szs (concat ms) - - g [Just nr,Just nc] m - | nr == r && nc == c = m - | r == 1 && c == 1 = matrixFromVector RowMajor nr nc (constantD x (nr*nc)) - | r == 1 = fromRows (replicate nr (flatten m)) - | otherwise = fromColumns (replicate nc (flatten m)) - where - r = rows m - c = cols m - x = m@@>(0,0) - g _ _ = error "inconsistent dimensions in fromBlocks" - - --------------------------------------------------------------------------------- - -{- | create a block diagonal matrix - ->>> disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]] -7x8 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -0 0 2 2 2 2 2 0 -0 0 2 2 2 2 2 0 -0 0 2 2 2 2 2 0 -0 0 0 0 0 0 0 5 -0 0 0 0 0 0 0 7 - ->>> diagBlock [(0><4)[], konst 2 (2,3)] :: Matrix Double -(2><7) - [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 - , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ] - --} -diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t -diagBlock ms = fromBlocks $ zipWith f ms [0..] - where - f m k = take n $ replicate k z ++ m : repeat z - n = length ms - z = (1><1) [0] - --------------------------------------------------------------------------------- - - --- | Reverse rows -flipud :: Element t => Matrix t -> Matrix t -flipud m = extractRows [r-1,r-2 .. 0] $ m - where - r = rows m - --- | Reverse columns -fliprl :: Element t => Matrix t -> Matrix t -fliprl m = extractColumns [c-1,c-2 .. 0] $ m - where - c = cols m - ------------------------------------------------------------- - -{- | creates a rectangular diagonal matrix: - ->>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double -(4><5) - [ 10.0, 7.0, 7.0, 7.0, 7.0 - , 7.0, 20.0, 7.0, 7.0, 7.0 - , 7.0, 7.0, 30.0, 7.0, 7.0 - , 7.0, 7.0, 7.0, 7.0, 7.0 ] - --} -diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t -diagRect z v r c = ST.runSTMatrix $ do - m <- ST.newMatrix z r c - let d = min r c `min` (dim v) - mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1] - return m - --- | extracts the diagonal from a rectangular matrix -takeDiag :: (Element t) => Matrix t -> Vector t -takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] - ------------------------------------------------------------- - -{- | An easy way to create a matrix: - ->>> (2><3)[2,4,7,-3,11,0] -(2><3) - [ 2.0, 4.0, 7.0 - , -3.0, 11.0, 0.0 ] - -This is the format produced by the instances of Show (Matrix a), which -can also be used for input. - -The input list is explicitly truncated, so that it can -safely be used with lists that are too long (like infinite lists). - ->>> (2><3)[1..] -(2><3) - [ 1.0, 2.0, 3.0 - , 4.0, 5.0, 6.0 ] - - --} -(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a -r >< c = f where - f l | dim v == r*c = matrixFromVector RowMajor r c v - | otherwise = error $ "inconsistent list size = " - ++show (dim v) ++" in ("++show r++"><"++show c++")" - where v = fromList $ take (r*c) l - ----------------------------------------------------------------- - --- | Creates a matrix with the first n rows of another matrix -takeRows :: Element t => Int -> Matrix t -> Matrix t -takeRows n mt = subMatrix (0,0) (n, cols mt) mt --- | Creates a copy of a matrix without the first n rows -dropRows :: Element t => Int -> Matrix t -> Matrix t -dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt --- |Creates a matrix with the first n columns of another matrix -takeColumns :: Element t => Int -> Matrix t -> Matrix t -takeColumns n mt = subMatrix (0,0) (rows mt, n) mt --- | Creates a copy of a matrix without the first n columns -dropColumns :: Element t => Int -> Matrix t -> Matrix t -dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt - ----------------------------------------------------------------- - -{- | Creates a 'Matrix' from a list of lists (considered as rows). - ->>> fromLists [[1,2],[3,4],[5,6]] -(3><2) - [ 1.0, 2.0 - , 3.0, 4.0 - , 5.0, 6.0 ] - --} -fromLists :: Element t => [[t]] -> Matrix t -fromLists = fromRows . map fromList - --- | creates a 1-row matrix from a vector --- --- >>> asRow (fromList [1..5]) --- (1><5) --- [ 1.0, 2.0, 3.0, 4.0, 5.0 ] --- -asRow :: Storable a => Vector a -> Matrix a -asRow v = reshape (dim v) v - --- | creates a 1-column matrix from a vector --- --- >>> asColumn (fromList [1..5]) --- (5><1) --- [ 1.0 --- , 2.0 --- , 3.0 --- , 4.0 --- , 5.0 ] --- -asColumn :: Storable a => Vector a -> Matrix a -asColumn = trans . asRow - - - -{- | creates a Matrix of the specified size using the supplied function to - to map the row\/column position to the value at that row\/column position. - -@> buildMatrix 3 4 (\\(r,c) -> fromIntegral r * fromIntegral c) -(3><4) - [ 0.0, 0.0, 0.0, 0.0, 0.0 - , 0.0, 1.0, 2.0, 3.0, 4.0 - , 0.0, 2.0, 4.0, 6.0, 8.0]@ - -Hilbert matrix of order N: - -@hilb n = buildMatrix n n (\\(i,j)->1/(fromIntegral i + fromIntegral j +1))@ - --} -buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a -buildMatrix rc cc f = - fromLists $ map (map f) - $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)] - ------------------------------------------------------ - -fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e -fromArray2D m = (r> [Int] -> Matrix t -> Matrix t -extractRows [] m = emptyM 0 (cols m) -extractRows l m = fromRows $ extract (toRows m) l - where - extract l' is = [l'!!i | i<- map verify is] - verify k - | k >= 0 && k < rows m = k - | otherwise = error $ "can't extract row " - ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m - --- | rearranges the rows of a matrix according to the order given in a list of integers. -extractColumns :: Element t => [Int] -> Matrix t -> Matrix t -extractColumns l m = trans . extractRows (map verify l) . trans $ m - where - verify k - | k >= 0 && k < cols m = k - | otherwise = error $ "can't extract column " - ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m - - - -{- | creates matrix by repetition of a matrix a given number of rows and columns - ->>> repmat (ident 2) 2 3 -(4><6) - [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 - , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 - , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 - , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ] - --} -repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t -repmat m r c - | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m) - | otherwise = fromBlocks $ replicate r $ replicate c $ m - --- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix. -liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t -liftMatrix2Auto f m1 m2 - | compat' m1 m2 = lM f m1 m2 - | ok = lM f m1' m2' - | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2 - where - (r1,c1) = size m1 - (r2,c2) = size m2 - r = max r1 r2 - c = max c1 c2 - r0 = min r1 r2 - c0 = min c1 c2 - ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2 - m1' = conformMTo (r,c) m1 - m2' = conformMTo (r,c) m2 - --- FIXME do not flatten if equal order -lM f m1 m2 = matrixFromVector - RowMajor - (max (rows m1) (rows m2)) - (max (cols m1) (cols m2)) - (f (flatten m1) (flatten m2)) - -compat' :: Matrix a -> Matrix b -> Bool -compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2 - where - s1 = size m1 - s2 = size m2 - ------------------------------------------------------------- - -toBlockRows [r] m | r == rows m = [m] -toBlockRows rs m = map (reshape (cols m)) (takesV szs (flatten m)) - where szs = map (* cols m) rs - -toBlockCols [c] m | c == cols m = [m] -toBlockCols cs m = map trans . toBlockRows cs . trans $ m - --- | Partition a matrix into blocks with the given numbers of rows and columns. --- The remaining rows and columns are discarded. -toBlocks :: (Element t) => [Int] -> [Int] -> Matrix t -> [[Matrix t]] -toBlocks rs cs m = map (toBlockCols cs) . toBlockRows rs $ m - --- | Fully partition a matrix into blocks of the same size. If the dimensions are not --- a multiple of the given size the last blocks will be smaller. -toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]] -toBlocksEvery r c m - | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c - | otherwise = toBlocks rs cs m - where - (qr,rr) = rows m `divMod` r - (qc,rc) = cols m `divMod` c - rs = replicate qr r ++ if rr > 0 then [rr] else [] - cs = replicate qc c ++ if rc > 0 then [rc] else [] - -------------------------------------------------------------------- - --- Given a column number and a function taking matrix indexes, returns --- a function which takes vector indexes (that can be used on the --- flattened matrix). -mk :: Int -> ((Int, Int) -> t) -> (Int -> t) -mk c g = \k -> g (divMod k c) - -{- | - ->>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..]) -m[0,0] = 1 -m[0,1] = 2 -m[0,2] = 3 -m[1,0] = 4 -m[1,1] = 5 -m[1,2] = 6 - --} -mapMatrixWithIndexM_ - :: (Element a, Num a, Monad m) => - ((Int, Int) -> a -> m ()) -> Matrix a -> m () -mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m - where - c = cols m - -{- | - ->>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double) -Just (3><3) - [ 100.0, 1.0, 2.0 - , 10.0, 111.0, 12.0 - , 20.0, 21.0, 122.0 ] - --} -mapMatrixWithIndexM - :: (Element a, Storable b, Monad m) => - ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b) -mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m - where - c = cols m - -{- | - ->>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double) -(3><3) - [ 100.0, 1.0, 2.0 - , 10.0, 111.0, 12.0 - , 20.0, 21.0, 122.0 ] - - -} -mapMatrixWithIndex - :: (Element a, Storable b) => - ((Int, Int) -> a -> b) -> Matrix a -> Matrix b -mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m - where - c = cols m - -mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b -mapMatrix f = liftMatrix (mapVector f) diff --git a/packages/hmatrix/src/Data/Packed/ST.hs b/packages/hmatrix/src/Data/Packed/ST.hs deleted file mode 100644 index 1cef296..0000000 --- a/packages/hmatrix/src/Data/Packed/ST.hs +++ /dev/null @@ -1,179 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE Rank2Types #-} -{-# LANGUAGE BangPatterns #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.ST --- Copyright : (c) Alberto Ruiz 2008 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- In-place manipulation inside the ST monad. --- See examples/inplace.hs in the distribution. --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Packed.ST ( - -- * Mutable Vectors - STVector, newVector, thawVector, freezeVector, runSTVector, - readVector, writeVector, modifyVector, liftSTVector, - -- * Mutable Matrices - STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix, - readMatrix, writeMatrix, modifyMatrix, liftSTMatrix, - -- * Unsafe functions - newUndefinedVector, - unsafeReadVector, unsafeWriteVector, - unsafeThawVector, unsafeFreezeVector, - newUndefinedMatrix, - unsafeReadMatrix, unsafeWriteMatrix, - unsafeThawMatrix, unsafeFreezeMatrix -) where - -import Data.Packed.Internal - -import Control.Monad.ST(ST, runST) -import Foreign.Storable(Storable, peekElemOff, pokeElemOff) - -#if MIN_VERSION_base(4,4,0) -import Control.Monad.ST.Unsafe(unsafeIOToST) -#else -import Control.Monad.ST(unsafeIOToST) -#endif - -{-# INLINE ioReadV #-} -ioReadV :: Storable t => Vector t -> Int -> IO t -ioReadV v k = unsafeWith v $ \s -> peekElemOff s k - -{-# INLINE ioWriteV #-} -ioWriteV :: Storable t => Vector t -> Int -> t -> IO () -ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x - -newtype STVector s t = STVector (Vector t) - -thawVector :: Storable t => Vector t -> ST s (STVector s t) -thawVector = unsafeIOToST . fmap STVector . cloneVector - -unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t) -unsafeThawVector = unsafeIOToST . return . STVector - -runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t -runSTVector st = runST (st >>= unsafeFreezeVector) - -{-# INLINE unsafeReadVector #-} -unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t -unsafeReadVector (STVector x) = unsafeIOToST . ioReadV x - -{-# INLINE unsafeWriteVector #-} -unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s () -unsafeWriteVector (STVector x) k = unsafeIOToST . ioWriteV x k - -{-# INLINE modifyVector #-} -modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s () -modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k - -liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a -liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x - -freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t) -freezeVector v = liftSTVector id v - -unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t) -unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x - -{-# INLINE safeIndexV #-} -safeIndexV f (STVector v) k - | k < 0 || k>= dim v = error $ "out of range error in vector (dim=" - ++show (dim v)++", pos="++show k++")" - | otherwise = f (STVector v) k - -{-# INLINE readVector #-} -readVector :: Storable t => STVector s t -> Int -> ST s t -readVector = safeIndexV unsafeReadVector - -{-# INLINE writeVector #-} -writeVector :: Storable t => STVector s t -> Int -> t -> ST s () -writeVector = safeIndexV unsafeWriteVector - -newUndefinedVector :: Storable t => Int -> ST s (STVector s t) -newUndefinedVector = unsafeIOToST . fmap STVector . createVector - -{-# INLINE newVector #-} -newVector :: Storable t => t -> Int -> ST s (STVector s t) -newVector x n = do - v <- newUndefinedVector n - let go (-1) = return v - go !k = unsafeWriteVector v k x >> go (k-1 :: Int) - go (n-1) - -------------------------------------------------------------------------- - -{-# INLINE ioReadM #-} -ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t -ioReadM (Matrix _ nc cv RowMajor) r c = ioReadV cv (r*nc+c) -ioReadM (Matrix nr _ fv ColumnMajor) r c = ioReadV fv (c*nr+r) - -{-# INLINE ioWriteM #-} -ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO () -ioWriteM (Matrix _ nc cv RowMajor) r c val = ioWriteV cv (r*nc+c) val -ioWriteM (Matrix nr _ fv ColumnMajor) r c val = ioWriteV fv (c*nr+r) val - -newtype STMatrix s t = STMatrix (Matrix t) - -thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t) -thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix - -unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t) -unsafeThawMatrix = unsafeIOToST . return . STMatrix - -runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t -runSTMatrix st = runST (st >>= unsafeFreezeMatrix) - -{-# INLINE unsafeReadMatrix #-} -unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t -unsafeReadMatrix (STMatrix x) r = unsafeIOToST . ioReadM x r - -{-# INLINE unsafeWriteMatrix #-} -unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s () -unsafeWriteMatrix (STMatrix x) r c = unsafeIOToST . ioWriteM x r c - -{-# INLINE modifyMatrix #-} -modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s () -modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c - -liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a -liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x - -unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t) -unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x - -freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t) -freezeMatrix m = liftSTMatrix id m - -cloneMatrix (Matrix r c d o) = cloneVector d >>= return . (\d' -> Matrix r c d' o) - -{-# INLINE safeIndexM #-} -safeIndexM f (STMatrix m) r c - | r<0 || r>=rows m || - c<0 || c>=cols m = error $ "out of range error in matrix (size=" - ++show (rows m,cols m)++", pos="++show (r,c)++")" - | otherwise = f (STMatrix m) r c - -{-# INLINE readMatrix #-} -readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t -readMatrix = safeIndexM unsafeReadMatrix - -{-# INLINE writeMatrix #-} -writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s () -writeMatrix = safeIndexM unsafeWriteMatrix - -newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t) -newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c - -{-# NOINLINE newMatrix #-} -newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t) -newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c) diff --git a/packages/hmatrix/src/Data/Packed/Vector.hs b/packages/hmatrix/src/Data/Packed/Vector.hs deleted file mode 100644 index b5a4318..0000000 --- a/packages/hmatrix/src/Data/Packed/Vector.hs +++ /dev/null @@ -1,96 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE CPP #-} ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Vector --- Copyright : (c) Alberto Ruiz 2007-10 --- License : GPL --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- --- 1D arrays suitable for numeric computations using external libraries. --- --- This module provides basic functions for manipulation of structure. --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Packed.Vector ( - Vector, - fromList, (|>), toList, buildVector, - dim, (@>), - subVector, takesV, vjoin, join, - mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, - mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, - foldLoop, foldVector, foldVectorG, foldVectorWithIndex -) where - -import Data.Packed.Internal.Vector -import Foreign.Storable - -------------------------------------------------------------------- - -#ifdef BINARY - -import Data.Binary -import Control.Monad(replicateM) - --- a 64K cache, with a Double taking 13 bytes in Bytestring, --- implies a chunk size of 5041 -chunk :: Int -chunk = 5000 - -chunks :: Int -> [Int] -chunks d = let c = d `div` chunk - m = d `mod` chunk - in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk) - -putVector v = do - let d = dim v - mapM_ (\i -> put $ v @> i) [0..(d-1)] - -getVector d = do - xs <- replicateM d get - return $! fromList xs - -instance (Binary a, Storable a) => Binary (Vector a) where - put v = do - let d = dim v - put d - mapM_ putVector $! takesV (chunks d) v - get = do - d <- get - vs <- mapM getVector $ chunks d - return $! vjoin vs - -#endif - -------------------------------------------------------------------- - -{- | creates a Vector of the specified length using the supplied function to - to map the index to the value at that index. - -@> buildVector 4 fromIntegral -4 |> [0.0,1.0,2.0,3.0]@ - --} -buildVector :: Storable a => Int -> (Int -> a) -> Vector a -buildVector len f = - fromList $ map f [0 .. (len - 1)] - - --- | zip for Vectors -zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b) -zipVector = zipVectorWith (,) - --- | unzip for Vectors -unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b) -unzipVector = unzipVectorWith id - -------------------------------------------------------------------- - -{-# DEPRECATED join "use vjoin or Data.Vector.concat" #-} -join :: Storable t => [Vector t] -> Vector t -join = vjoin - -- cgit v1.2.3