From 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 08:48:12 +0200 Subject: empty hmatrix-base --- 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/Random.hs | 57 +++ packages/hmatrix/src/Data/Packed/ST.hs | 179 +++++++ packages/hmatrix/src/Data/Packed/Vector.hs | 96 ++++ 12 files changed, 2264 insertions(+) create mode 100644 packages/hmatrix/src/Data/Packed.hs create mode 100644 packages/hmatrix/src/Data/Packed/Development.hs create mode 100644 packages/hmatrix/src/Data/Packed/Foreign.hs create mode 100644 packages/hmatrix/src/Data/Packed/Internal.hs create mode 100644 packages/hmatrix/src/Data/Packed/Internal/Common.hs create mode 100644 packages/hmatrix/src/Data/Packed/Internal/Matrix.hs create mode 100644 packages/hmatrix/src/Data/Packed/Internal/Signatures.hs create mode 100644 packages/hmatrix/src/Data/Packed/Internal/Vector.hs create mode 100644 packages/hmatrix/src/Data/Packed/Matrix.hs create mode 100644 packages/hmatrix/src/Data/Packed/Random.hs create mode 100644 packages/hmatrix/src/Data/Packed/ST.hs create mode 100644 packages/hmatrix/src/Data/Packed/Vector.hs (limited to 'packages/hmatrix/src/Data') diff --git a/packages/hmatrix/src/Data/Packed.hs b/packages/hmatrix/src/Data/Packed.hs new file mode 100644 index 0000000..957aab8 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed.hs @@ -0,0 +1,29 @@ +----------------------------------------------------------------------------- +{- | +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 new file mode 100644 index 0000000..471e560 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Development.hs @@ -0,0 +1,32 @@ + +----------------------------------------------------------------------------- +-- | +-- 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 new file mode 100644 index 0000000..1ec3694 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Foreign.hs @@ -0,0 +1,100 @@ +{-# 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 new file mode 100644 index 0000000..537e51e --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Internal.hs @@ -0,0 +1,26 @@ +----------------------------------------------------------------------------- +-- | +-- 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 new file mode 100644 index 0000000..edef3c2 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Internal/Common.hs @@ -0,0 +1,171 @@ +{-# 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 new file mode 100644 index 0000000..9719fc0 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs @@ -0,0 +1,491 @@ +{-# 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 new file mode 100644 index 0000000..2835720 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs @@ -0,0 +1,72 @@ +----------------------------------------------------------------------------- +-- | +-- 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 new file mode 100644 index 0000000..6d03438 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs @@ -0,0 +1,521 @@ +{-# 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 new file mode 100644 index 0000000..d94d167 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Matrix.hs @@ -0,0 +1,490 @@ +{-# 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/Random.hs b/packages/hmatrix/src/Data/Packed/Random.hs new file mode 100644 index 0000000..e8b0268 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Random.hs @@ -0,0 +1,57 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Packed.Vector +-- Copyright : (c) Alberto Ruiz 2009 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Random vectors and matrices. +-- +----------------------------------------------------------------------------- + +module Data.Packed.Random ( + Seed, + RandDist(..), + randomVector, + gaussianSample, + uniformSample +) where + +import Numeric.GSL.Vector +import Data.Packed +import Numeric.ContainerBoot +import Numeric.LinearAlgebra.Algorithms + + +type Seed = Int + +-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate +-- Gaussian distribution. +gaussianSample :: Seed + -> Int -- ^ number of rows + -> Vector Double -- ^ mean vector + -> Matrix Double -- ^ covariance matrix + -> Matrix Double -- ^ result +gaussianSample seed n med cov = m where + c = dim med + meds = konst' 1 n `outer` med + rs = reshape c $ randomVector seed Gaussian (c * n) + m = rs `mXm` cholSH cov `add` meds + +-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate +-- uniform distribution. +uniformSample :: Seed + -> Int -- ^ number of rows + -> [(Double,Double)] -- ^ ranges for each column + -> Matrix Double -- ^ result +uniformSample seed n rgs = m where + (as,bs) = unzip rgs + a = fromList as + cs = zipWith subtract as bs + d = dim a + dat = toRows $ reshape n $ randomVector seed Uniform (n*d) + am = konst' 1 n `outer` a + m = fromColumns (zipWith scale cs dat) `add` am + diff --git a/packages/hmatrix/src/Data/Packed/ST.hs b/packages/hmatrix/src/Data/Packed/ST.hs new file mode 100644 index 0000000..1cef296 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/ST.hs @@ -0,0 +1,179 @@ +{-# 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 new file mode 100644 index 0000000..b5a4318 --- /dev/null +++ b/packages/hmatrix/src/Data/Packed/Vector.hs @@ -0,0 +1,96 @@ +{-# 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