summaryrefslogtreecommitdiff
path: root/lib/Data
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
commit1925c123d7d8184a1d2ddc0a413e0fd2776e1083 (patch)
treefad79f909d9c3be53d68e6ebd67202650536d387 /lib/Data
parenteb3f702d065a4a967bb754977233e6eec408fd1f (diff)
empty hmatrix-base
Diffstat (limited to 'lib/Data')
-rw-r--r--lib/Data/Packed.hs29
-rw-r--r--lib/Data/Packed/Development.hs32
-rw-r--r--lib/Data/Packed/Foreign.hs100
-rw-r--r--lib/Data/Packed/Internal.hs26
-rw-r--r--lib/Data/Packed/Internal/Common.hs171
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs491
-rw-r--r--lib/Data/Packed/Internal/Signatures.hs72
-rw-r--r--lib/Data/Packed/Internal/Vector.hs521
-rw-r--r--lib/Data/Packed/Matrix.hs490
-rw-r--r--lib/Data/Packed/Random.hs57
-rw-r--r--lib/Data/Packed/ST.hs179
-rw-r--r--lib/Data/Packed/Vector.hs96
12 files changed, 0 insertions, 2264 deletions
diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
deleted file mode 100644
index 957aab8..0000000
--- a/lib/Data/Packed.hs
+++ /dev/null
@@ -1,29 +0,0 @@
1-----------------------------------------------------------------------------
2{- |
3Module : Data.Packed
4Copyright : (c) Alberto Ruiz 2006-2010
5License : GPL-style
6
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional
9Portability : uses ffi
10
11Types for dense 'Vector' and 'Matrix' of 'Storable' elements.
12
13-}
14-----------------------------------------------------------------------------
15{-# OPTIONS_HADDOCK hide #-}
16
17module Data.Packed (
18 module Data.Packed.Vector,
19 module Data.Packed.Matrix,
20-- module Numeric.Conversion,
21-- module Data.Packed.Random,
22-- module Data.Complex
23) where
24
25import Data.Packed.Vector
26import Data.Packed.Matrix
27--import Data.Packed.Random
28--import Data.Complex
29--import Numeric.Conversion
diff --git a/lib/Data/Packed/Development.hs b/lib/Data/Packed/Development.hs
deleted file mode 100644
index 471e560..0000000
--- a/lib/Data/Packed/Development.hs
+++ /dev/null
@@ -1,32 +0,0 @@
1
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Development
5-- Copyright : (c) Alberto Ruiz 2009
6-- License : GPL
7--
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional
10-- Portability : portable
11--
12-- The library can be easily extended with additional foreign functions
13-- using the tools in this module. Illustrative usage examples can be found
14-- in the @examples\/devel@ folder included in the package.
15--
16-----------------------------------------------------------------------------
17{-# OPTIONS_HADDOCK hide #-}
18
19module Data.Packed.Development (
20 createVector, createMatrix,
21 vec, mat,
22 app1, app2, app3, app4,
23 app5, app6, app7, app8, app9, app10,
24 MatrixOrder(..), orderOf, cmat, fmat,
25 matrixFromVector,
26 unsafeFromForeignPtr,
27 unsafeToForeignPtr,
28 check, (//),
29 at', atM'
30) where
31
32import Data.Packed.Internal
diff --git a/lib/Data/Packed/Foreign.hs b/lib/Data/Packed/Foreign.hs
deleted file mode 100644
index 1ec3694..0000000
--- a/lib/Data/Packed/Foreign.hs
+++ /dev/null
@@ -1,100 +0,0 @@
1{-# LANGUAGE MagicHash, UnboxedTuples #-}
2-- | FFI and hmatrix helpers.
3--
4-- Sample usage, to upload a perspective matrix to a shader.
5--
6-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
7-- @
8--
9{-# OPTIONS_HADDOCK hide #-}
10module Data.Packed.Foreign
11 ( app
12 , appVector, appVectorLen
13 , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
14 , unsafeMatrixToVector, unsafeMatrixToForeignPtr
15 ) where
16import Data.Packed.Internal
17import qualified Data.Vector.Storable as S
18import Foreign (Ptr, ForeignPtr, Storable)
19import Foreign.C.Types (CInt)
20import GHC.Base (IO(..), realWorld#)
21
22{-# INLINE unsafeInlinePerformIO #-}
23-- | 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
24-- unecessary calls to unsafePerformIO or its internals.
25unsafeInlinePerformIO :: IO a -> a
26unsafeInlinePerformIO (IO f) = case f realWorld# of
27 (# _, x #) -> x
28
29{-# INLINE app #-}
30-- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
31-- e.g.
32--
33-- @
34-- someFunction
35-- \`appMatrixLen\` m
36-- \`appVectorLen\` v
37-- \`app\` other
38-- \`app\` arguments
39-- \`app\` go here
40-- @
41--
42-- One could also write:
43--
44-- @
45-- (someFunction
46-- \`appMatrixLen\` m
47-- \`appVectorLen\` v)
48-- other
49-- arguments
50-- (go here)
51-- @
52--
53app :: (a -> b) -> a -> b
54app f = f
55
56{-# INLINE appVector #-}
57appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
58appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
59
60{-# INLINE appVectorLen #-}
61appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
62appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
63
64{-# INLINE appMatrix #-}
65appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
66appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
67
68{-# INLINE appMatrixLen #-}
69appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
70appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
71 where
72 r = fromIntegral (rows x)
73 c = fromIntegral (cols x)
74
75{-# INLINE appMatrixRaw #-}
76appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
77appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
78
79{-# INLINE appMatrixRawLen #-}
80appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
81appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
82 where
83 r = fromIntegral (rows x)
84 c = fromIntegral (cols x)
85
86infixl 1 `app`
87infixl 1 `appVector`
88infixl 1 `appMatrix`
89infixl 1 `appMatrixRaw`
90
91{-# INLINE unsafeMatrixToVector #-}
92-- | This will disregard the order of the matrix, and simply return it as-is.
93-- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
94unsafeMatrixToVector :: Matrix a -> Vector a
95unsafeMatrixToVector = xdat
96
97{-# INLINE unsafeMatrixToForeignPtr #-}
98unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
99unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
100
diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs
deleted file mode 100644
index 537e51e..0000000
--- a/lib/Data/Packed/Internal.hs
+++ /dev/null
@@ -1,26 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Internal
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9-- Portability : portable
10--
11-- Reexports all internal modules
12--
13-----------------------------------------------------------------------------
14-- #hide
15
16module Data.Packed.Internal (
17 module Data.Packed.Internal.Common,
18 module Data.Packed.Internal.Signatures,
19 module Data.Packed.Internal.Vector,
20 module Data.Packed.Internal.Matrix,
21) where
22
23import Data.Packed.Internal.Common
24import Data.Packed.Internal.Signatures
25import Data.Packed.Internal.Vector
26import Data.Packed.Internal.Matrix
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
deleted file mode 100644
index edef3c2..0000000
--- a/lib/Data/Packed/Internal/Common.hs
+++ /dev/null
@@ -1,171 +0,0 @@
1{-# LANGUAGE CPP #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Internal.Common
5-- Copyright : (c) Alberto Ruiz 2007
6-- License : GPL-style
7--
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional
10-- Portability : portable (uses FFI)
11--
12-- Development utilities.
13--
14-----------------------------------------------------------------------------
15-- #hide
16
17module Data.Packed.Internal.Common(
18 Adapt,
19 app1, app2, app3, app4,
20 app5, app6, app7, app8, app9, app10,
21 (//), check, mbCatch,
22 splitEvery, common, compatdim,
23 fi,
24 table
25) where
26
27import Foreign
28import Control.Monad(when)
29import Foreign.C.String(peekCString)
30import Foreign.C.Types
31import Foreign.Storable.Complex()
32import Data.List(transpose,intersperse)
33import Control.Exception as E
34
35-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
36splitEvery :: Int -> [a] -> [[a]]
37splitEvery _ [] = []
38splitEvery k l = take k l : splitEvery k (drop k l)
39
40-- | obtains the common value of a property of a list
41common :: (Eq a) => (b->a) -> [b] -> Maybe a
42common f = commonval . map f where
43 commonval :: (Eq a) => [a] -> Maybe a
44 commonval [] = Nothing
45 commonval [a] = Just a
46 commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
47
48-- | common value with \"adaptable\" 1
49compatdim :: [Int] -> Maybe Int
50compatdim [] = Nothing
51compatdim [a] = Just a
52compatdim (a:b:xs)
53 | a==b = compatdim (b:xs)
54 | a==1 = compatdim (b:xs)
55 | b==1 = compatdim (a:xs)
56 | otherwise = Nothing
57
58-- | Formatting tool
59table :: String -> [[String]] -> String
60table sep as = unlines . map unwords' $ transpose mtp where
61 mt = transpose as
62 longs = map (maximum . map length) mt
63 mtp = zipWith (\a b -> map (pad a) b) longs mt
64 pad n str = replicate (n - length str) ' ' ++ str
65 unwords' = concat . intersperse sep
66
67-- | postfix function application (@flip ($)@)
68(//) :: x -> (x -> y) -> y
69infixl 0 //
70(//) = flip ($)
71
72-- | specialized fromIntegral
73fi :: Int -> CInt
74fi = fromIntegral
75
76-- hmm..
77ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
78ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
79ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
80ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
81ww6 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
82ww7 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
83ww8 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
84ww9 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
85ww10 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
86
87type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
88
89type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
90type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
91type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
92type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
93type 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
94type 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
95type 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
96type 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
97type 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
98type 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
99
100app1 :: f -> Adapt1 f t1
101app2 :: f -> Adapt2 f t1 r1 t2
102app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
103app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
104app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
105app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
106app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
107app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
108app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
109app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
110
111app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
112app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
113app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
114 \a1 a2 a3 -> f // a1 // a2 // a3 // check s
115app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $
116 \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
117app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $
118 \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s
119app6 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 $
120 \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s
121app7 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 $
122 \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s
123app8 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 $
124 \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s
125app9 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 $
126 \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s
127app10 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 $
128 \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s
129
130
131
132-- GSL error codes are <= 1024
133-- | error codes for the auxiliary functions required by the wrappers
134errorCode :: CInt -> String
135errorCode 2000 = "bad size"
136errorCode 2001 = "bad function code"
137errorCode 2002 = "memory problem"
138errorCode 2003 = "bad file"
139errorCode 2004 = "singular"
140errorCode 2005 = "didn't converge"
141errorCode 2006 = "the input matrix is not positive definite"
142errorCode 2007 = "not yet supported in this OS"
143errorCode n = "code "++show n
144
145
146-- | clear the fpu
147foreign import ccall unsafe "asm_finit" finit :: IO ()
148
149-- | check the error code
150check :: String -> IO CInt -> IO ()
151check msg f = do
152#if FINIT
153 finit
154#endif
155 err <- f
156 when (err/=0) $ if err > 1024
157 then (error (msg++": "++errorCode err)) -- our errors
158 else do -- GSL errors
159 ps <- gsl_strerror err
160 s <- peekCString ps
161 error (msg++": "++s)
162 return ()
163
164-- | description of GSL error codes
165foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
166
167-- | Error capture and conversion to Maybe
168mbCatch :: IO x -> IO (Maybe x)
169mbCatch act = E.catch (Just `fmap` act) f
170 where f :: SomeException -> IO (Maybe x)
171 f _ = return Nothing
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
deleted file mode 100644
index 9719fc0..0000000
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ /dev/null
@@ -1,491 +0,0 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE BangPatterns #-}
5-----------------------------------------------------------------------------
6-- |
7-- Module : Data.Packed.Internal.Matrix
8-- Copyright : (c) Alberto Ruiz 2007
9-- License : GPL-style
10--
11-- Maintainer : Alberto Ruiz <aruiz@um.es>
12-- Stability : provisional
13-- Portability : portable (uses FFI)
14--
15-- Internal matrix representation
16--
17-----------------------------------------------------------------------------
18-- #hide
19
20module Data.Packed.Internal.Matrix(
21 Matrix(..), rows, cols, cdat, fdat,
22 MatrixOrder(..), orderOf,
23 createMatrix, mat,
24 cmat, fmat,
25 toLists, flatten, reshape,
26 Element(..),
27 trans,
28 fromRows, toRows, fromColumns, toColumns,
29 matrixFromVector,
30 subMatrix,
31 liftMatrix, liftMatrix2,
32 (@@>), atM',
33 saveMatrix,
34 singleton,
35 emptyM,
36 size, shSize, conformVs, conformMs, conformVTo, conformMTo
37) where
38
39import Data.Packed.Internal.Common
40import Data.Packed.Internal.Signatures
41import Data.Packed.Internal.Vector
42
43import Foreign.Marshal.Alloc(alloca, free)
44import Foreign.Marshal.Array(newArray)
45import Foreign.Ptr(Ptr, castPtr)
46import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
47import Data.Complex(Complex)
48import Foreign.C.Types
49import Foreign.C.String(newCString)
50import System.IO.Unsafe(unsafePerformIO)
51import Control.DeepSeq
52
53-----------------------------------------------------------------
54
55{- Design considerations for the Matrix Type
56 -----------------------------------------
57
58- we must easily handle both row major and column major order,
59 for bindings to LAPACK and GSL/C
60
61- we'd like to simplify redundant matrix transposes:
62 - Some of them arise from the order requirements of some functions
63 - some functions (matrix product) admit transposed arguments
64
65- maybe we don't really need this kind of simplification:
66 - more complex code
67 - some computational overhead
68 - only appreciable gain in code with a lot of redundant transpositions
69 and cheap matrix computations
70
71- we could carry both the matrix and its (lazily computed) transpose.
72 This may save some transpositions, but it is necessary to keep track of the
73 data which is actually computed to be used by functions like the matrix product
74 which admit both orders.
75
76- but if we need the transposed data and it is not in the structure, we must make
77 sure that we touch the same foreignptr that is used in the computation.
78
79- a reasonable solution is using two constructors for a matrix. Transposition just
80 "flips" the constructor. Actual data transposition is not done if followed by a
81 matrix product or another transpose.
82
83-}
84
85data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
86
87transOrder RowMajor = ColumnMajor
88transOrder ColumnMajor = RowMajor
89{- | Matrix representation suitable for GSL and LAPACK computations.
90
91The elements are stored in a continuous memory array.
92
93-}
94
95data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int
96 , icols :: {-# UNPACK #-} !Int
97 , xdat :: {-# UNPACK #-} !(Vector t)
98 , order :: !MatrixOrder }
99-- RowMajor: preferred by C, fdat may require a transposition
100-- ColumnMajor: preferred by LAPACK, cdat may require a transposition
101
102cdat = xdat
103fdat = xdat
104
105rows :: Matrix t -> Int
106rows = irows
107
108cols :: Matrix t -> Int
109cols = icols
110
111orderOf :: Matrix t -> MatrixOrder
112orderOf = order
113
114
115-- | Matrix transpose.
116trans :: Matrix t -> Matrix t
117trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o}
118
119cmat :: (Element t) => Matrix t -> Matrix t
120cmat m@Matrix{order = RowMajor} = m
121cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor}
122
123fmat :: (Element t) => Matrix t -> Matrix t
124fmat m@Matrix{order = ColumnMajor} = m
125fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor}
126
127-- C-Haskell matrix adapter
128-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
129
130mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
131mat a f =
132 unsafeWith (xdat a) $ \p -> do
133 let m g = do
134 g (fi (rows a)) (fi (cols a)) p
135 f m
136
137{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
138
139>>> flatten (ident 3)
140fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
141
142-}
143flatten :: Element t => Matrix t -> Vector t
144flatten = xdat . cmat
145
146{-
147type Mt t s = Int -> Int -> Ptr t -> s
148
149infixr 6 ::>
150type t ::> s = Mt t s
151-}
152
153-- | the inverse of 'Data.Packed.Matrix.fromLists'
154toLists :: (Element t) => Matrix t -> [[t]]
155toLists m = splitEvery (cols m) . toList . flatten $ m
156
157-- | Create a matrix from a list of vectors.
158-- All vectors must have the same dimension,
159-- or dimension 1, which is are automatically expanded.
160fromRows :: Element t => [Vector t] -> Matrix t
161fromRows [] = emptyM 0 0
162fromRows vs = case compatdim (map dim vs) of
163 Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
164 Just 0 -> emptyM r 0
165 Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
166 where
167 r = length vs
168 adapt c v
169 | c == 0 = fromList[]
170 | dim v == c = v
171 | otherwise = constantD (v@>0) c
172
173-- | extracts the rows of a matrix as a list of vectors
174toRows :: Element t => Matrix t -> [Vector t]
175toRows m
176 | c == 0 = replicate r (fromList[])
177 | otherwise = toRows' 0
178 where
179 v = flatten m
180 r = rows m
181 c = cols m
182 toRows' k | k == r*c = []
183 | otherwise = subVector k c v : toRows' (k+c)
184
185-- | Creates a matrix from a list of vectors, as columns
186fromColumns :: Element t => [Vector t] -> Matrix t
187fromColumns m = trans . fromRows $ m
188
189-- | Creates a list of vectors from the columns of a matrix
190toColumns :: Element t => Matrix t -> [Vector t]
191toColumns m = toRows . trans $ m
192
193-- | Reads a matrix position.
194(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
195infixl 9 @@>
196m@Matrix {irows = r, icols = c} @@> (i,j)
197 | safe = if i<0 || i>=r || j<0 || j>=c
198 then error "matrix indexing out of range"
199 else atM' m i j
200 | otherwise = atM' m i j
201{-# INLINE (@@>) #-}
202
203-- Unsafe matrix access without range checking
204atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j)
205atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
206{-# INLINE atM' #-}
207
208------------------------------------------------------------------
209
210matrixFromVector o r c v
211 | r * c == dim v = m
212 | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
213 where
214 m = Matrix { irows = r, icols = c, xdat = v, order = o }
215
216-- allocates memory for a new matrix
217createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
218createMatrix ord r c = do
219 p <- createVector (r*c)
220 return (matrixFromVector ord r c p)
221
222{- | 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@
223where r is the desired number of rows.)
224
225>>> reshape 4 (fromList [1..12])
226(3><4)
227 [ 1.0, 2.0, 3.0, 4.0
228 , 5.0, 6.0, 7.0, 8.0
229 , 9.0, 10.0, 11.0, 12.0 ]
230
231-}
232reshape :: Storable t => Int -> Vector t -> Matrix t
233reshape 0 v = matrixFromVector RowMajor 0 0 v
234reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
235
236singleton x = reshape 1 (fromList [x])
237
238-- | application of a vector function on the flattened matrix elements
239liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
240liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d)
241
242-- | application of a vector function on the flattened matrices elements
243liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
244liftMatrix2 f m1 m2
245 | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
246 | otherwise = case orderOf m1 of
247 RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2))
248 ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2))
249
250
251compat :: Matrix a -> Matrix b -> Bool
252compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
253
254------------------------------------------------------------------
255
256{- | Supported matrix elements.
257
258 This class provides optimized internal
259 operations for selected element types.
260 It provides unoptimised defaults for any 'Storable' type,
261 so you can create instances simply as:
262 @instance Element Foo@.
263-}
264class (Storable a) => Element a where
265 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
266 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
267 -> Matrix a -> Matrix a
268 subMatrixD = subMatrix'
269 transdata :: Int -> Vector a -> Int -> Vector a
270 transdata = transdataP -- transdata'
271 constantD :: a -> Int -> Vector a
272 constantD = constantP -- constant'
273
274
275instance Element Float where
276 transdata = transdataAux ctransF
277 constantD = constantAux cconstantF
278
279instance Element Double where
280 transdata = transdataAux ctransR
281 constantD = constantAux cconstantR
282
283instance Element (Complex Float) where
284 transdata = transdataAux ctransQ
285 constantD = constantAux cconstantQ
286
287instance Element (Complex Double) where
288 transdata = transdataAux ctransC
289 constantD = constantAux cconstantC
290
291-------------------------------------------------------------------
292
293transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
294transdata' c1 v c2 =
295 if noneed
296 then v
297 else unsafePerformIO $ do
298 w <- createVector (r2*c2)
299 unsafeWith v $ \p ->
300 unsafeWith w $ \q -> do
301 let go (-1) _ = return ()
302 go !i (-1) = go (i-1) (c1-1)
303 go !i !j = do x <- peekElemOff p (i*c1+j)
304 pokeElemOff q (j*c2+i) x
305 go i (j-1)
306 go (r1-1) (c1-1)
307 return w
308 where r1 = dim v `div` c1
309 r2 = dim v `div` c2
310 noneed = dim v == 0 || r1 == 1 || c1 == 1
311
312-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int -> Vector Double #-}
313-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
314
315-- I don't know how to specialize...
316-- The above pragmas only seem to work on top level defs
317-- Fortunately everything seems to work using the above class
318
319-- C versions, still a little faster:
320
321transdataAux fun c1 d c2 =
322 if noneed
323 then d
324 else unsafePerformIO $ do
325 v <- createVector (dim d)
326 unsafeWith d $ \pd ->
327 unsafeWith v $ \pv ->
328 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
329 return v
330 where r1 = dim d `div` c1
331 r2 = dim d `div` c2
332 noneed = dim d == 0 || r1 == 1 || c1 == 1
333
334transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
335transdataP c1 d c2 =
336 if noneed
337 then d
338 else unsafePerformIO $ do
339 v <- createVector (dim d)
340 unsafeWith d $ \pd ->
341 unsafeWith v $ \pv ->
342 ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"
343 return v
344 where r1 = dim d `div` c1
345 r2 = dim d `div` c2
346 sz = sizeOf (d @> 0)
347 noneed = dim d == 0 || r1 == 1 || c1 == 1
348
349foreign import ccall unsafe "transF" ctransF :: TFMFM
350foreign import ccall unsafe "transR" ctransR :: TMM
351foreign import ccall unsafe "transQ" ctransQ :: TQMQM
352foreign import ccall unsafe "transC" ctransC :: TCMCM
353foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
354
355----------------------------------------------------------------------
356
357constant' v n = unsafePerformIO $ do
358 w <- createVector n
359 unsafeWith w $ \p -> do
360 let go (-1) = return ()
361 go !k = pokeElemOff p k v >> go (k-1)
362 go (n-1)
363 return w
364
365-- C versions
366
367constantAux fun x n = unsafePerformIO $ do
368 v <- createVector n
369 px <- newArray [x]
370 app1 (fun px) vec v "constantAux"
371 free px
372 return v
373
374constantF :: Float -> Int -> Vector Float
375constantF = constantAux cconstantF
376foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
377
378constantR :: Double -> Int -> Vector Double
379constantR = constantAux cconstantR
380foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
381
382constantQ :: Complex Float -> Int -> Vector (Complex Float)
383constantQ = constantAux cconstantQ
384foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
385
386constantC :: Complex Double -> Int -> Vector (Complex Double)
387constantC = constantAux cconstantC
388foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
389
390constantP :: Storable a => a -> Int -> Vector a
391constantP a n = unsafePerformIO $ do
392 let sz = sizeOf a
393 v <- createVector n
394 unsafeWith v $ \p -> do
395 alloca $ \k -> do
396 poke k a
397 cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
398 return v
399foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
400
401----------------------------------------------------------------------
402
403-- | Extracts a submatrix from a matrix.
404subMatrix :: Element a
405 => (Int,Int) -- ^ (r0,c0) starting position
406 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
407 -> Matrix a -- ^ input matrix
408 -> Matrix a -- ^ result
409subMatrix (r0,c0) (rt,ct) m
410 | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
411 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
412 | otherwise = error $ "wrong subMatrix "++
413 show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
414
415subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
416 w <- createVector (rt*ct)
417 unsafeWith v $ \p ->
418 unsafeWith w $ \q -> do
419 let go (-1) _ = return ()
420 go !i (-1) = go (i-1) (ct-1)
421 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)
422 pokeElemOff q (i*ct+j) x
423 go i (j-1)
424 go (rt-1) (ct-1)
425 return w
426
427subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor
428subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
429
430--------------------------------------------------------------------------
431
432-- | Saves a matrix as 2D ASCII table.
433saveMatrix :: FilePath
434 -> String -- ^ format (%f, %g, %e)
435 -> Matrix Double
436 -> IO ()
437saveMatrix filename fmt m = do
438 charname <- newCString filename
439 charfmt <- newCString fmt
440 let o = if orderOf m == RowMajor then 1 else 0
441 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
442 free charname
443 free charfmt
444
445foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
446
447----------------------------------------------------------------------
448
449maxZ xs = if minimum xs == 0 then 0 else maximum xs
450
451conformMs ms = map (conformMTo (r,c)) ms
452 where
453 r = maxZ (map rows ms)
454 c = maxZ (map cols ms)
455
456
457conformVs vs = map (conformVTo n) vs
458 where
459 n = maxZ (map dim vs)
460
461conformMTo (r,c) m
462 | size m == (r,c) = m
463 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
464 | size m == (r,1) = repCols c m
465 | size m == (1,c) = repRows r m
466 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
467
468conformVTo n v
469 | dim v == n = v
470 | dim v == 1 = constantD (v@>0) n
471 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
472
473repRows n x = fromRows (replicate n (flatten x))
474repCols n x = fromColumns (replicate n (flatten x))
475
476size m = (rows m, cols m)
477
478shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
479
480emptyM r c = matrixFromVector RowMajor r c (fromList[])
481
482----------------------------------------------------------------------
483
484instance (Storable t, NFData t) => NFData (Matrix t)
485 where
486 rnf m | d > 0 = rnf (v @> 0)
487 | otherwise = ()
488 where
489 d = dim v
490 v = xdat m
491
diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs
deleted file mode 100644
index 2835720..0000000
--- a/lib/Data/Packed/Internal/Signatures.hs
+++ /dev/null
@@ -1,72 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Internal.Signatures
4-- Copyright : (c) Alberto Ruiz 2009
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9-- Portability : portable (uses FFI)
10--
11-- Signatures of the C functions.
12--
13-----------------------------------------------------------------------------
14
15module Data.Packed.Internal.Signatures where
16
17import Foreign.Ptr(Ptr)
18import Data.Complex(Complex)
19import Foreign.C.Types(CInt)
20
21type PF = Ptr Float --
22type PD = Ptr Double --
23type PQ = Ptr (Complex Float) --
24type PC = Ptr (Complex Double) --
25type TF = CInt -> PF -> IO CInt --
26type TFF = CInt -> PF -> TF --
27type TFV = CInt -> PF -> TV --
28type TVF = CInt -> PD -> TF --
29type TFFF = CInt -> PF -> TFF --
30type TV = CInt -> PD -> IO CInt --
31type TVV = CInt -> PD -> TV --
32type TVVV = CInt -> PD -> TVV --
33type TFM = CInt -> CInt -> PF -> IO CInt --
34type TFMFM = CInt -> CInt -> PF -> TFM --
35type TFMFMFM = CInt -> CInt -> PF -> TFMFM --
36type TM = CInt -> CInt -> PD -> IO CInt --
37type TMM = CInt -> CInt -> PD -> TM --
38type TVMM = CInt -> PD -> TMM --
39type TMVMM = CInt -> CInt -> PD -> TVMM --
40type TMMM = CInt -> CInt -> PD -> TMM --
41type TVM = CInt -> PD -> TM --
42type TVVM = CInt -> PD -> TVM --
43type TMV = CInt -> CInt -> PD -> TV --
44type TMMV = CInt -> CInt -> PD -> TMV --
45type TMVM = CInt -> CInt -> PD -> TVM --
46type TMMVM = CInt -> CInt -> PD -> TMVM --
47type TCM = CInt -> CInt -> PC -> IO CInt --
48type TCVCM = CInt -> PC -> TCM --
49type TCMCVCM = CInt -> CInt -> PC -> TCVCM --
50type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM --
51type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM --
52type TCMCM = CInt -> CInt -> PC -> TCM --
53type TVCM = CInt -> PD -> TCM --
54type TCMVCM = CInt -> CInt -> PC -> TVCM --
55type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM --
56type TCMCMCM = CInt -> CInt -> PC -> TCMCM --
57type TCV = CInt -> PC -> IO CInt --
58type TCVCV = CInt -> PC -> TCV --
59type TCVCVCV = CInt -> PC -> TCVCV --
60type TCVV = CInt -> PC -> TV --
61type TQV = CInt -> PQ -> IO CInt --
62type TQVQV = CInt -> PQ -> TQV --
63type TQVQVQV = CInt -> PQ -> TQVQV --
64type TQVF = CInt -> PQ -> TF --
65type TQM = CInt -> CInt -> PQ -> IO CInt --
66type TQMQM = CInt -> CInt -> PQ -> TQM --
67type TQMQMQM = CInt -> CInt -> PQ -> TQMQM --
68type TCMCV = CInt -> CInt -> PC -> TCV --
69type TVCV = CInt -> PD -> TCV --
70type TCVM = CInt -> PC -> TM --
71type TMCVM = CInt -> CInt -> PD -> TCVM --
72type TMMCVM = CInt -> CInt -> PD -> TMCVM --
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
deleted file mode 100644
index 6d03438..0000000
--- a/lib/Data/Packed/Internal/Vector.hs
+++ /dev/null
@@ -1,521 +0,0 @@
1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Internal.Vector
5-- Copyright : (c) Alberto Ruiz 2007
6-- License : GPL-style
7--
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional
10-- Portability : portable (uses FFI)
11--
12-- Vector implementation
13--
14-----------------------------------------------------------------------------
15
16module Data.Packed.Internal.Vector (
17 Vector, dim,
18 fromList, toList, (|>),
19 vjoin, (@>), safe, at, at', subVector, takesV,
20 mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith,
21 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
22 foldVector, foldVectorG, foldLoop, foldVectorWithIndex,
23 createVector, vec,
24 asComplex, asReal, float2DoubleV, double2FloatV,
25 stepF, stepD, condF, condD,
26 conjugateQ, conjugateC,
27 fwriteVector, freadVector, fprintfVector, fscanfVector,
28 cloneVector,
29 unsafeToForeignPtr,
30 unsafeFromForeignPtr,
31 unsafeWith
32) where
33
34import Data.Packed.Internal.Common
35import Data.Packed.Internal.Signatures
36import Foreign.Marshal.Alloc(free)
37import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
38import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
39import Foreign.Ptr(Ptr)
40import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
41import Foreign.C.String
42import Foreign.C.Types
43import Data.Complex
44import Control.Monad(when)
45import System.IO.Unsafe(unsafePerformIO)
46
47#if __GLASGOW_HASKELL__ >= 605
48import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
49#else
50import Foreign.ForeignPtr (mallocForeignPtrBytes)
51#endif
52
53import GHC.Base
54#if __GLASGOW_HASKELL__ < 612
55import GHC.IOBase hiding (liftIO)
56#endif
57
58import qualified Data.Vector.Storable as Vector
59import Data.Vector.Storable(Vector,
60 fromList,
61 unsafeToForeignPtr,
62 unsafeFromForeignPtr,
63 unsafeWith)
64
65
66-- | Number of elements
67dim :: (Storable t) => Vector t -> Int
68dim = Vector.length
69
70
71-- C-Haskell vector adapter
72-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
73vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
74vec x f = unsafeWith x $ \p -> do
75 let v g = do
76 g (fi $ dim x) p
77 f v
78{-# INLINE vec #-}
79
80
81-- allocates memory for a new vector
82createVector :: Storable a => Int -> IO (Vector a)
83createVector n = do
84 when (n < 0) $ error ("trying to createVector of negative dim: "++show n)
85 fp <- doMalloc undefined
86 return $ unsafeFromForeignPtr fp 0 n
87 where
88 --
89 -- Use the much cheaper Haskell heap allocated storage
90 -- for foreign pointer space we control
91 --
92 doMalloc :: Storable b => b -> IO (ForeignPtr b)
93 doMalloc dummy = do
94#if __GLASGOW_HASKELL__ >= 605
95 mallocPlainForeignPtrBytes (n * sizeOf dummy)
96#else
97 mallocForeignPtrBytes (n * sizeOf dummy)
98#endif
99
100{- | creates a Vector from a list:
101
102@> fromList [2,3,5,7]
1034 |> [2.0,3.0,5.0,7.0]@
104
105-}
106
107safeRead v = inlinePerformIO . unsafeWith v
108{-# INLINE safeRead #-}
109
110inlinePerformIO :: IO a -> a
111inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
112{-# INLINE inlinePerformIO #-}
113
114{- | extracts the Vector elements to a list
115
116>>> toList (linspace 5 (1,10))
117[1.0,3.25,5.5,7.75,10.0]
118
119-}
120toList :: Storable a => Vector a -> [a]
121toList v = safeRead v $ peekArray (dim v)
122
123{- | Create a vector from a list of elements and explicit dimension. The input
124 list is explicitly truncated if it is too long, so it may safely
125 be used, for instance, with infinite lists.
126
127>>> 5 |> [1..]
128fromList [1.0,2.0,3.0,4.0,5.0]
129
130-}
131(|>) :: (Storable a) => Int -> [a] -> Vector a
132infixl 9 |>
133n |> l = if length l' == n
134 then fromList l'
135 else error "list too short for |>"
136 where l' = take n l
137
138
139-- | access to Vector elements without range checking
140at' :: Storable a => Vector a -> Int -> a
141at' v n = safeRead v $ flip peekElemOff n
142{-# INLINE at' #-}
143
144--
145-- turn off bounds checking with -funsafe at configure time.
146-- ghc will optimise away the salways true case at compile time.
147--
148#if defined(UNSAFE)
149safe :: Bool
150safe = False
151#else
152safe = True
153#endif
154
155-- | access to Vector elements with range checking.
156at :: Storable a => Vector a -> Int -> a
157at v n
158 | safe = if n >= 0 && n < dim v
159 then at' v n
160 else error "vector index out of range"
161 | otherwise = at' v n
162{-# INLINE at #-}
163
164{- | takes a number of consecutive elements from a Vector
165
166>>> subVector 2 3 (fromList [1..10])
167fromList [3.0,4.0,5.0]
168
169-}
170subVector :: Storable t => Int -- ^ index of the starting element
171 -> Int -- ^ number of elements to extract
172 -> Vector t -- ^ source
173 -> Vector t -- ^ result
174subVector = Vector.slice
175
176
177{- | Reads a vector position:
178
179>>> fromList [0..9] @> 7
1807.0
181
182-}
183(@>) :: Storable t => Vector t -> Int -> t
184infixl 9 @>
185(@>) = at
186
187
188{- | concatenate a list of vectors
189
190>>> vjoin [fromList [1..5::Double], konst 1 3]
191fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
192
193-}
194vjoin :: Storable t => [Vector t] -> Vector t
195vjoin [] = fromList []
196vjoin [v] = v
197vjoin as = unsafePerformIO $ do
198 let tot = sum (map dim as)
199 r <- createVector tot
200 unsafeWith r $ \ptr ->
201 joiner as tot ptr
202 return r
203 where joiner [] _ _ = return ()
204 joiner (v:cs) _ p = do
205 let n = dim v
206 unsafeWith v $ \pb -> copyArray p pb n
207 joiner cs 0 (advancePtr p n)
208
209
210{- | Extract consecutive subvectors of the given sizes.
211
212>>> takesV [3,4] (linspace 10 (1,10::Double))
213[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
214
215-}
216takesV :: Storable t => [Int] -> Vector t -> [Vector t]
217takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w)
218 | otherwise = go ms w
219 where go [] _ = []
220 go (n:ns) v = subVector 0 n v
221 : go ns (subVector n (dim v - n) v)
222
223---------------------------------------------------------------
224
225-- | transforms a complex vector into a real vector with alternating real and imaginary parts
226asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
227asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
228 where (fp,i,n) = unsafeToForeignPtr v
229
230-- | transforms a real vector into a complex vector with alternating real and imaginary parts
231asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
232asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
233 where (fp,i,n) = unsafeToForeignPtr v
234
235---------------------------------------------------------------
236
237float2DoubleV :: Vector Float -> Vector Double
238float2DoubleV v = unsafePerformIO $ do
239 r <- createVector (dim v)
240 app2 c_float2double vec v vec r "float2double"
241 return r
242
243double2FloatV :: Vector Double -> Vector Float
244double2FloatV v = unsafePerformIO $ do
245 r <- createVector (dim v)
246 app2 c_double2float vec v vec r "double2float2"
247 return r
248
249
250foreign import ccall unsafe "float2double" c_float2double:: TFV
251foreign import ccall unsafe "double2float" c_double2float:: TVF
252
253---------------------------------------------------------------
254
255stepF :: Vector Float -> Vector Float
256stepF v = unsafePerformIO $ do
257 r <- createVector (dim v)
258 app2 c_stepF vec v vec r "stepF"
259 return r
260
261stepD :: Vector Double -> Vector Double
262stepD v = unsafePerformIO $ do
263 r <- createVector (dim v)
264 app2 c_stepD vec v vec r "stepD"
265 return r
266
267foreign import ccall unsafe "stepF" c_stepF :: TFF
268foreign import ccall unsafe "stepD" c_stepD :: TVV
269
270---------------------------------------------------------------
271
272condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
273condF x y l e g = unsafePerformIO $ do
274 r <- createVector (dim x)
275 app6 c_condF vec x vec y vec l vec e vec g vec r "condF"
276 return r
277
278condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
279condD x y l e g = unsafePerformIO $ do
280 r <- createVector (dim x)
281 app6 c_condD vec x vec y vec l vec e vec g vec r "condD"
282 return r
283
284foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF
285foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV
286
287--------------------------------------------------------------------------------
288
289conjugateAux fun x = unsafePerformIO $ do
290 v <- createVector (dim x)
291 app2 fun vec x vec v "conjugateAux"
292 return v
293
294conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
295conjugateQ = conjugateAux c_conjugateQ
296foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV
297
298conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
299conjugateC = conjugateAux c_conjugateC
300foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV
301
302--------------------------------------------------------------------------------
303
304cloneVector :: Storable t => Vector t -> IO (Vector t)
305cloneVector v = do
306 let n = dim v
307 r <- createVector n
308 let f _ s _ d = copyArray d s n >> return 0
309 app2 f vec v vec r "cloneVector"
310 return r
311
312------------------------------------------------------------------
313
314-- | map on Vectors
315mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
316mapVector f v = unsafePerformIO $ do
317 w <- createVector (dim v)
318 unsafeWith v $ \p ->
319 unsafeWith w $ \q -> do
320 let go (-1) = return ()
321 go !k = do x <- peekElemOff p k
322 pokeElemOff q k (f x)
323 go (k-1)
324 go (dim v -1)
325 return w
326{-# INLINE mapVector #-}
327
328-- | zipWith for Vectors
329zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
330zipVectorWith f u v = unsafePerformIO $ do
331 let n = min (dim u) (dim v)
332 w <- createVector n
333 unsafeWith u $ \pu ->
334 unsafeWith v $ \pv ->
335 unsafeWith w $ \pw -> do
336 let go (-1) = return ()
337 go !k = do x <- peekElemOff pu k
338 y <- peekElemOff pv k
339 pokeElemOff pw k (f x y)
340 go (k-1)
341 go (n -1)
342 return w
343{-# INLINE zipVectorWith #-}
344
345-- | unzipWith for Vectors
346unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
347 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
348unzipVectorWith f u = unsafePerformIO $ do
349 let n = dim u
350 v <- createVector n
351 w <- createVector n
352 unsafeWith u $ \pu ->
353 unsafeWith v $ \pv ->
354 unsafeWith w $ \pw -> do
355 let go (-1) = return ()
356 go !k = do z <- peekElemOff pu k
357 let (x,y) = f z
358 pokeElemOff pv k x
359 pokeElemOff pw k y
360 go (k-1)
361 go (n-1)
362 return (v,w)
363{-# INLINE unzipVectorWith #-}
364
365foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
366foldVector f x v = unsafePerformIO $
367 unsafeWith v $ \p -> do
368 let go (-1) s = return s
369 go !k !s = do y <- peekElemOff p k
370 go (k-1::Int) (f y s)
371 go (dim v -1) x
372{-# INLINE foldVector #-}
373
374-- the zero-indexed index is passed to the folding function
375foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
376foldVectorWithIndex f x v = unsafePerformIO $
377 unsafeWith v $ \p -> do
378 let go (-1) s = return s
379 go !k !s = do y <- peekElemOff p k
380 go (k-1::Int) (f k y s)
381 go (dim v -1) x
382{-# INLINE foldVectorWithIndex #-}
383
384foldLoop f s0 d = go (d - 1) s0
385 where
386 go 0 s = f (0::Int) s
387 go !j !s = go (j - 1) (f j s)
388
389foldVectorG f s0 v = foldLoop g s0 (dim v)
390 where g !k !s = f k (at' v) s
391 {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
392{-# INLINE foldVectorG #-}
393
394-------------------------------------------------------------------
395
396-- | monadic map over Vectors
397-- the monad @m@ must be strict
398mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
399mapVectorM f v = do
400 w <- return $! unsafePerformIO $! createVector (dim v)
401 mapVectorM' w 0 (dim v -1)
402 return w
403 where mapVectorM' w' !k !t
404 | k == t = do
405 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
406 y <- f x
407 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
408 | otherwise = do
409 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
410 y <- f x
411 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
412 mapVectorM' w' (k+1) t
413{-# INLINE mapVectorM #-}
414
415-- | monadic map over Vectors
416mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
417mapVectorM_ f v = do
418 mapVectorM' 0 (dim v -1)
419 where mapVectorM' !k !t
420 | k == t = do
421 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
422 f x
423 | otherwise = do
424 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
425 _ <- f x
426 mapVectorM' (k+1) t
427{-# INLINE mapVectorM_ #-}
428
429-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
430-- the monad @m@ must be strict
431mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
432mapVectorWithIndexM f v = do
433 w <- return $! unsafePerformIO $! createVector (dim v)
434 mapVectorM' w 0 (dim v -1)
435 return w
436 where mapVectorM' w' !k !t
437 | k == t = do
438 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
439 y <- f k x
440 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
441 | otherwise = do
442 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
443 y <- f k x
444 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
445 mapVectorM' w' (k+1) t
446{-# INLINE mapVectorWithIndexM #-}
447
448-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
449mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
450mapVectorWithIndexM_ f v = do
451 mapVectorM' 0 (dim v -1)
452 where mapVectorM' !k !t
453 | k == t = do
454 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
455 f k x
456 | otherwise = do
457 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
458 _ <- f k x
459 mapVectorM' (k+1) t
460{-# INLINE mapVectorWithIndexM_ #-}
461
462
463mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
464--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
465mapVectorWithIndex f v = unsafePerformIO $ do
466 w <- createVector (dim v)
467 unsafeWith v $ \p ->
468 unsafeWith w $ \q -> do
469 let go (-1) = return ()
470 go !k = do x <- peekElemOff p k
471 pokeElemOff q k (f k x)
472 go (k-1)
473 go (dim v -1)
474 return w
475{-# INLINE mapVectorWithIndex #-}
476
477-------------------------------------------------------------------
478
479
480-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
481fscanfVector :: FilePath -> Int -> IO (Vector Double)
482fscanfVector filename n = do
483 charname <- newCString filename
484 res <- createVector n
485 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
486 free charname
487 return res
488
489foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
490
491-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
492fprintfVector :: FilePath -> String -> Vector Double -> IO ()
493fprintfVector filename fmt v = do
494 charname <- newCString filename
495 charfmt <- newCString fmt
496 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
497 free charname
498 free charfmt
499
500foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
501
502-- | Loads a vector from a binary file (the number of elements must be known in advance).
503freadVector :: FilePath -> Int -> IO (Vector Double)
504freadVector filename n = do
505 charname <- newCString filename
506 res <- createVector n
507 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
508 free charname
509 return res
510
511foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
512
513-- | Saves the elements of a vector to a binary file.
514fwriteVector :: FilePath -> Vector Double -> IO ()
515fwriteVector filename v = do
516 charname <- newCString filename
517 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
518 free charname
519
520foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
521
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
deleted file mode 100644
index d94d167..0000000
--- a/lib/Data/Packed/Matrix.hs
+++ /dev/null
@@ -1,490 +0,0 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE CPP #-}
6
7-----------------------------------------------------------------------------
8-- |
9-- Module : Data.Packed.Matrix
10-- Copyright : (c) Alberto Ruiz 2007-10
11-- License : GPL
12--
13-- Maintainer : Alberto Ruiz <aruiz@um.es>
14-- Stability : provisional
15--
16-- A Matrix representation suitable for numerical computations using LAPACK and GSL.
17--
18-- This module provides basic functions for manipulation of structure.
19
20-----------------------------------------------------------------------------
21{-# OPTIONS_HADDOCK hide #-}
22
23module Data.Packed.Matrix (
24 Matrix,
25 Element,
26 rows,cols,
27 (><),
28 trans,
29 reshape, flatten,
30 fromLists, toLists, buildMatrix,
31 (@@>),
32 asRow, asColumn,
33 fromRows, toRows, fromColumns, toColumns,
34 fromBlocks, diagBlock, toBlocks, toBlocksEvery,
35 repmat,
36 flipud, fliprl,
37 subMatrix, takeRows, dropRows, takeColumns, dropColumns,
38 extractRows, extractColumns,
39 diagRect, takeDiag,
40 mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
41 liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
42) where
43
44import Data.Packed.Internal
45import qualified Data.Packed.ST as ST
46import Data.Array
47
48import Data.List(transpose,intersperse)
49import Foreign.Storable(Storable)
50import Control.Monad(liftM)
51
52-------------------------------------------------------------------
53
54#ifdef BINARY
55
56import Data.Binary
57import Control.Monad(replicateM)
58
59instance (Binary a, Element a, Storable a) => Binary (Matrix a) where
60 put m = do
61 let r = rows m
62 let c = cols m
63 put r
64 put c
65 mapM_ (\i -> mapM_ (\j -> put $ m @@> (i,j)) [0..(c-1)]) [0..(r-1)]
66 get = do
67 r <- get
68 c <- get
69 xs <- replicateM r $ replicateM c get
70 return $ fromLists xs
71
72#endif
73
74-------------------------------------------------------------------
75
76instance (Show a, Element a) => (Show (Matrix a)) where
77 show m | rows m == 0 || cols m == 0 = sizes m ++" []"
78 show m = (sizes m++) . dsp . map (map show) . toLists $ m
79
80sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
81
82dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
83 where
84 mt = transpose as
85 longs = map (maximum . map length) mt
86 mtp = zipWith (\a b -> map (pad a) b) longs mt
87 pad n str = replicate (n - length str) ' ' ++ str
88 unwords' = concat . intersperse ", "
89
90------------------------------------------------------------------
91
92instance (Element a, Read a) => Read (Matrix a) where
93 readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
94 where (thing,rest) = breakAt ']' s
95 (dims,listnums) = breakAt ')' thing
96 cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
97 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
98
99
100breakAt c l = (a++[c],tail b) where
101 (a,b) = break (==c) l
102
103------------------------------------------------------------------
104
105-- | creates a matrix from a vertical list of matrices
106joinVert :: Element t => [Matrix t] -> Matrix t
107joinVert [] = emptyM 0 0
108joinVert ms = case common cols ms of
109 Nothing -> error "(impossible) joinVert on matrices with different number of columns"
110 Just c -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms)
111
112-- | creates a matrix from a horizontal list of matrices
113joinHoriz :: Element t => [Matrix t] -> Matrix t
114joinHoriz ms = trans. joinVert . map trans $ ms
115
116{- | Create a matrix from blocks given as a list of lists of matrices.
117
118Single row-column components are automatically expanded to match the
119corresponding common row and column:
120
121@
122disp = putStr . dispf 2
123@
124
125>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
1268x10
1271 0 0 0 0 7 7 7 10 20
1280 1 0 0 0 7 7 7 10 20
1290 0 1 0 0 7 7 7 10 20
1300 0 0 1 0 7 7 7 10 20
1310 0 0 0 1 7 7 7 10 20
1323 3 3 3 3 1 0 0 0 0
1333 3 3 3 3 0 2 0 0 0
1343 3 3 3 3 0 0 3 0 0
135
136-}
137fromBlocks :: Element t => [[Matrix t]] -> Matrix t
138fromBlocks = fromBlocksRaw . adaptBlocks
139
140fromBlocksRaw mms = joinVert . map joinHoriz $ mms
141
142adaptBlocks ms = ms' where
143 bc = case common length ms of
144 Just c -> c
145 Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
146 rs = map (compatdim . map rows) ms
147 cs = map (compatdim . map cols) (transpose ms)
148 szs = sequence [rs,cs]
149 ms' = splitEvery bc $ zipWith g szs (concat ms)
150
151 g [Just nr,Just nc] m
152 | nr == r && nc == c = m
153 | r == 1 && c == 1 = matrixFromVector RowMajor nr nc (constantD x (nr*nc))
154 | r == 1 = fromRows (replicate nr (flatten m))
155 | otherwise = fromColumns (replicate nc (flatten m))
156 where
157 r = rows m
158 c = cols m
159 x = m@@>(0,0)
160 g _ _ = error "inconsistent dimensions in fromBlocks"
161
162
163--------------------------------------------------------------------------------
164
165{- | create a block diagonal matrix
166
167>>> disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
1687x8
1691 1 0 0 0 0 0 0
1701 1 0 0 0 0 0 0
1710 0 2 2 2 2 2 0
1720 0 2 2 2 2 2 0
1730 0 2 2 2 2 2 0
1740 0 0 0 0 0 0 5
1750 0 0 0 0 0 0 7
176
177>>> diagBlock [(0><4)[], konst 2 (2,3)] :: Matrix Double
178(2><7)
179 [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0
180 , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ]
181
182-}
183diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
184diagBlock ms = fromBlocks $ zipWith f ms [0..]
185 where
186 f m k = take n $ replicate k z ++ m : repeat z
187 n = length ms
188 z = (1><1) [0]
189
190--------------------------------------------------------------------------------
191
192
193-- | Reverse rows
194flipud :: Element t => Matrix t -> Matrix t
195flipud m = extractRows [r-1,r-2 .. 0] $ m
196 where
197 r = rows m
198
199-- | Reverse columns
200fliprl :: Element t => Matrix t -> Matrix t
201fliprl m = extractColumns [c-1,c-2 .. 0] $ m
202 where
203 c = cols m
204
205------------------------------------------------------------
206
207{- | creates a rectangular diagonal matrix:
208
209>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
210(4><5)
211 [ 10.0, 7.0, 7.0, 7.0, 7.0
212 , 7.0, 20.0, 7.0, 7.0, 7.0
213 , 7.0, 7.0, 30.0, 7.0, 7.0
214 , 7.0, 7.0, 7.0, 7.0, 7.0 ]
215
216-}
217diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
218diagRect z v r c = ST.runSTMatrix $ do
219 m <- ST.newMatrix z r c
220 let d = min r c `min` (dim v)
221 mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
222 return m
223
224-- | extracts the diagonal from a rectangular matrix
225takeDiag :: (Element t) => Matrix t -> Vector t
226takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
227
228------------------------------------------------------------
229
230{- | An easy way to create a matrix:
231
232>>> (2><3)[2,4,7,-3,11,0]
233(2><3)
234 [ 2.0, 4.0, 7.0
235 , -3.0, 11.0, 0.0 ]
236
237This is the format produced by the instances of Show (Matrix a), which
238can also be used for input.
239
240The input list is explicitly truncated, so that it can
241safely be used with lists that are too long (like infinite lists).
242
243>>> (2><3)[1..]
244(2><3)
245 [ 1.0, 2.0, 3.0
246 , 4.0, 5.0, 6.0 ]
247
248
249-}
250(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
251r >< c = f where
252 f l | dim v == r*c = matrixFromVector RowMajor r c v
253 | otherwise = error $ "inconsistent list size = "
254 ++show (dim v) ++" in ("++show r++"><"++show c++")"
255 where v = fromList $ take (r*c) l
256
257----------------------------------------------------------------
258
259-- | Creates a matrix with the first n rows of another matrix
260takeRows :: Element t => Int -> Matrix t -> Matrix t
261takeRows n mt = subMatrix (0,0) (n, cols mt) mt
262-- | Creates a copy of a matrix without the first n rows
263dropRows :: Element t => Int -> Matrix t -> Matrix t
264dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
265-- |Creates a matrix with the first n columns of another matrix
266takeColumns :: Element t => Int -> Matrix t -> Matrix t
267takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
268-- | Creates a copy of a matrix without the first n columns
269dropColumns :: Element t => Int -> Matrix t -> Matrix t
270dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
271
272----------------------------------------------------------------
273
274{- | Creates a 'Matrix' from a list of lists (considered as rows).
275
276>>> fromLists [[1,2],[3,4],[5,6]]
277(3><2)
278 [ 1.0, 2.0
279 , 3.0, 4.0
280 , 5.0, 6.0 ]
281
282-}
283fromLists :: Element t => [[t]] -> Matrix t
284fromLists = fromRows . map fromList
285
286-- | creates a 1-row matrix from a vector
287--
288-- >>> asRow (fromList [1..5])
289-- (1><5)
290-- [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
291--
292asRow :: Storable a => Vector a -> Matrix a
293asRow v = reshape (dim v) v
294
295-- | creates a 1-column matrix from a vector
296--
297-- >>> asColumn (fromList [1..5])
298-- (5><1)
299-- [ 1.0
300-- , 2.0
301-- , 3.0
302-- , 4.0
303-- , 5.0 ]
304--
305asColumn :: Storable a => Vector a -> Matrix a
306asColumn = trans . asRow
307
308
309
310{- | creates a Matrix of the specified size using the supplied function to
311 to map the row\/column position to the value at that row\/column position.
312
313@> buildMatrix 3 4 (\\(r,c) -> fromIntegral r * fromIntegral c)
314(3><4)
315 [ 0.0, 0.0, 0.0, 0.0, 0.0
316 , 0.0, 1.0, 2.0, 3.0, 4.0
317 , 0.0, 2.0, 4.0, 6.0, 8.0]@
318
319Hilbert matrix of order N:
320
321@hilb n = buildMatrix n n (\\(i,j)->1/(fromIntegral i + fromIntegral j +1))@
322
323-}
324buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a
325buildMatrix rc cc f =
326 fromLists $ map (map f)
327 $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)]
328
329-----------------------------------------------------
330
331fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e
332fromArray2D m = (r><c) (elems m)
333 where ((r0,c0),(r1,c1)) = bounds m
334 r = r1-r0+1
335 c = c1-c0+1
336
337
338-- | rearranges the rows of a matrix according to the order given in a list of integers.
339extractRows :: Element t => [Int] -> Matrix t -> Matrix t
340extractRows [] m = emptyM 0 (cols m)
341extractRows l m = fromRows $ extract (toRows m) l
342 where
343 extract l' is = [l'!!i | i<- map verify is]
344 verify k
345 | k >= 0 && k < rows m = k
346 | otherwise = error $ "can't extract row "
347 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
348
349-- | rearranges the rows of a matrix according to the order given in a list of integers.
350extractColumns :: Element t => [Int] -> Matrix t -> Matrix t
351extractColumns l m = trans . extractRows (map verify l) . trans $ m
352 where
353 verify k
354 | k >= 0 && k < cols m = k
355 | otherwise = error $ "can't extract column "
356 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
357
358
359
360{- | creates matrix by repetition of a matrix a given number of rows and columns
361
362>>> repmat (ident 2) 2 3
363(4><6)
364 [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
365 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
366 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
367 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
368
369-}
370repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
371repmat m r c
372 | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m)
373 | otherwise = fromBlocks $ replicate r $ replicate c $ m
374
375-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
376liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
377liftMatrix2Auto f m1 m2
378 | compat' m1 m2 = lM f m1 m2
379 | ok = lM f m1' m2'
380 | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2
381 where
382 (r1,c1) = size m1
383 (r2,c2) = size m2
384 r = max r1 r2
385 c = max c1 c2
386 r0 = min r1 r2
387 c0 = min c1 c2
388 ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2
389 m1' = conformMTo (r,c) m1
390 m2' = conformMTo (r,c) m2
391
392-- FIXME do not flatten if equal order
393lM f m1 m2 = matrixFromVector
394 RowMajor
395 (max (rows m1) (rows m2))
396 (max (cols m1) (cols m2))
397 (f (flatten m1) (flatten m2))
398
399compat' :: Matrix a -> Matrix b -> Bool
400compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
401 where
402 s1 = size m1
403 s2 = size m2
404
405------------------------------------------------------------
406
407toBlockRows [r] m | r == rows m = [m]
408toBlockRows rs m = map (reshape (cols m)) (takesV szs (flatten m))
409 where szs = map (* cols m) rs
410
411toBlockCols [c] m | c == cols m = [m]
412toBlockCols cs m = map trans . toBlockRows cs . trans $ m
413
414-- | Partition a matrix into blocks with the given numbers of rows and columns.
415-- The remaining rows and columns are discarded.
416toBlocks :: (Element t) => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
417toBlocks rs cs m = map (toBlockCols cs) . toBlockRows rs $ m
418
419-- | Fully partition a matrix into blocks of the same size. If the dimensions are not
420-- a multiple of the given size the last blocks will be smaller.
421toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]]
422toBlocksEvery r c m
423 | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c
424 | otherwise = toBlocks rs cs m
425 where
426 (qr,rr) = rows m `divMod` r
427 (qc,rc) = cols m `divMod` c
428 rs = replicate qr r ++ if rr > 0 then [rr] else []
429 cs = replicate qc c ++ if rc > 0 then [rc] else []
430
431-------------------------------------------------------------------
432
433-- Given a column number and a function taking matrix indexes, returns
434-- a function which takes vector indexes (that can be used on the
435-- flattened matrix).
436mk :: Int -> ((Int, Int) -> t) -> (Int -> t)
437mk c g = \k -> g (divMod k c)
438
439{- |
440
441>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
442m[0,0] = 1
443m[0,1] = 2
444m[0,2] = 3
445m[1,0] = 4
446m[1,1] = 5
447m[1,2] = 6
448
449-}
450mapMatrixWithIndexM_
451 :: (Element a, Num a, Monad m) =>
452 ((Int, Int) -> a -> m ()) -> Matrix a -> m ()
453mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
454 where
455 c = cols m
456
457{- |
458
459>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
460Just (3><3)
461 [ 100.0, 1.0, 2.0
462 , 10.0, 111.0, 12.0
463 , 20.0, 21.0, 122.0 ]
464
465-}
466mapMatrixWithIndexM
467 :: (Element a, Storable b, Monad m) =>
468 ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
469mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
470 where
471 c = cols m
472
473{- |
474
475>>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
476(3><3)
477 [ 100.0, 1.0, 2.0
478 , 10.0, 111.0, 12.0
479 , 20.0, 21.0, 122.0 ]
480
481 -}
482mapMatrixWithIndex
483 :: (Element a, Storable b) =>
484 ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
485mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
486 where
487 c = cols m
488
489mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
490mapMatrix f = liftMatrix (mapVector f)
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
deleted file mode 100644
index e8b0268..0000000
--- a/lib/Data/Packed/Random.hs
+++ /dev/null
@@ -1,57 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Vector
4-- Copyright : (c) Alberto Ruiz 2009
5-- License : GPL
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9--
10-- Random vectors and matrices.
11--
12-----------------------------------------------------------------------------
13
14module Data.Packed.Random (
15 Seed,
16 RandDist(..),
17 randomVector,
18 gaussianSample,
19 uniformSample
20) where
21
22import Numeric.GSL.Vector
23import Data.Packed
24import Numeric.ContainerBoot
25import Numeric.LinearAlgebra.Algorithms
26
27
28type Seed = Int
29
30-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
31-- Gaussian distribution.
32gaussianSample :: Seed
33 -> Int -- ^ number of rows
34 -> Vector Double -- ^ mean vector
35 -> Matrix Double -- ^ covariance matrix
36 -> Matrix Double -- ^ result
37gaussianSample seed n med cov = m where
38 c = dim med
39 meds = konst' 1 n `outer` med
40 rs = reshape c $ randomVector seed Gaussian (c * n)
41 m = rs `mXm` cholSH cov `add` meds
42
43-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
44-- uniform distribution.
45uniformSample :: Seed
46 -> Int -- ^ number of rows
47 -> [(Double,Double)] -- ^ ranges for each column
48 -> Matrix Double -- ^ result
49uniformSample seed n rgs = m where
50 (as,bs) = unzip rgs
51 a = fromList as
52 cs = zipWith subtract as bs
53 d = dim a
54 dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
55 am = konst' 1 n `outer` a
56 m = fromColumns (zipWith scale cs dat) `add` am
57
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
deleted file mode 100644
index 1cef296..0000000
--- a/lib/Data/Packed/ST.hs
+++ /dev/null
@@ -1,179 +0,0 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE TypeOperators #-}
3{-# LANGUAGE Rank2Types #-}
4{-# LANGUAGE BangPatterns #-}
5-----------------------------------------------------------------------------
6-- |
7-- Module : Data.Packed.ST
8-- Copyright : (c) Alberto Ruiz 2008
9-- License : GPL-style
10--
11-- Maintainer : Alberto Ruiz <aruiz@um.es>
12-- Stability : provisional
13-- Portability : portable
14--
15-- In-place manipulation inside the ST monad.
16-- See examples/inplace.hs in the distribution.
17--
18-----------------------------------------------------------------------------
19{-# OPTIONS_HADDOCK hide #-}
20
21module Data.Packed.ST (
22 -- * Mutable Vectors
23 STVector, newVector, thawVector, freezeVector, runSTVector,
24 readVector, writeVector, modifyVector, liftSTVector,
25 -- * Mutable Matrices
26 STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
27 readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
28 -- * Unsafe functions
29 newUndefinedVector,
30 unsafeReadVector, unsafeWriteVector,
31 unsafeThawVector, unsafeFreezeVector,
32 newUndefinedMatrix,
33 unsafeReadMatrix, unsafeWriteMatrix,
34 unsafeThawMatrix, unsafeFreezeMatrix
35) where
36
37import Data.Packed.Internal
38
39import Control.Monad.ST(ST, runST)
40import Foreign.Storable(Storable, peekElemOff, pokeElemOff)
41
42#if MIN_VERSION_base(4,4,0)
43import Control.Monad.ST.Unsafe(unsafeIOToST)
44#else
45import Control.Monad.ST(unsafeIOToST)
46#endif
47
48{-# INLINE ioReadV #-}
49ioReadV :: Storable t => Vector t -> Int -> IO t
50ioReadV v k = unsafeWith v $ \s -> peekElemOff s k
51
52{-# INLINE ioWriteV #-}
53ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()
54ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x
55
56newtype STVector s t = STVector (Vector t)
57
58thawVector :: Storable t => Vector t -> ST s (STVector s t)
59thawVector = unsafeIOToST . fmap STVector . cloneVector
60
61unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t)
62unsafeThawVector = unsafeIOToST . return . STVector
63
64runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t
65runSTVector st = runST (st >>= unsafeFreezeVector)
66
67{-# INLINE unsafeReadVector #-}
68unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t
69unsafeReadVector (STVector x) = unsafeIOToST . ioReadV x
70
71{-# INLINE unsafeWriteVector #-}
72unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s ()
73unsafeWriteVector (STVector x) k = unsafeIOToST . ioWriteV x k
74
75{-# INLINE modifyVector #-}
76modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()
77modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k
78
79liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a
80liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x
81
82freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
83freezeVector v = liftSTVector id v
84
85unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
86unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
87
88{-# INLINE safeIndexV #-}
89safeIndexV f (STVector v) k
90 | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
91 ++show (dim v)++", pos="++show k++")"
92 | otherwise = f (STVector v) k
93
94{-# INLINE readVector #-}
95readVector :: Storable t => STVector s t -> Int -> ST s t
96readVector = safeIndexV unsafeReadVector
97
98{-# INLINE writeVector #-}
99writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
100writeVector = safeIndexV unsafeWriteVector
101
102newUndefinedVector :: Storable t => Int -> ST s (STVector s t)
103newUndefinedVector = unsafeIOToST . fmap STVector . createVector
104
105{-# INLINE newVector #-}
106newVector :: Storable t => t -> Int -> ST s (STVector s t)
107newVector x n = do
108 v <- newUndefinedVector n
109 let go (-1) = return v
110 go !k = unsafeWriteVector v k x >> go (k-1 :: Int)
111 go (n-1)
112
113-------------------------------------------------------------------------
114
115{-# INLINE ioReadM #-}
116ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t
117ioReadM (Matrix _ nc cv RowMajor) r c = ioReadV cv (r*nc+c)
118ioReadM (Matrix nr _ fv ColumnMajor) r c = ioReadV fv (c*nr+r)
119
120{-# INLINE ioWriteM #-}
121ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()
122ioWriteM (Matrix _ nc cv RowMajor) r c val = ioWriteV cv (r*nc+c) val
123ioWriteM (Matrix nr _ fv ColumnMajor) r c val = ioWriteV fv (c*nr+r) val
124
125newtype STMatrix s t = STMatrix (Matrix t)
126
127thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
128thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix
129
130unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
131unsafeThawMatrix = unsafeIOToST . return . STMatrix
132
133runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t
134runSTMatrix st = runST (st >>= unsafeFreezeMatrix)
135
136{-# INLINE unsafeReadMatrix #-}
137unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
138unsafeReadMatrix (STMatrix x) r = unsafeIOToST . ioReadM x r
139
140{-# INLINE unsafeWriteMatrix #-}
141unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
142unsafeWriteMatrix (STMatrix x) r c = unsafeIOToST . ioWriteM x r c
143
144{-# INLINE modifyMatrix #-}
145modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
146modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c
147
148liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a
149liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x
150
151unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
152unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
153
154freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
155freezeMatrix m = liftSTMatrix id m
156
157cloneMatrix (Matrix r c d o) = cloneVector d >>= return . (\d' -> Matrix r c d' o)
158
159{-# INLINE safeIndexM #-}
160safeIndexM f (STMatrix m) r c
161 | r<0 || r>=rows m ||
162 c<0 || c>=cols m = error $ "out of range error in matrix (size="
163 ++show (rows m,cols m)++", pos="++show (r,c)++")"
164 | otherwise = f (STMatrix m) r c
165
166{-# INLINE readMatrix #-}
167readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
168readMatrix = safeIndexM unsafeReadMatrix
169
170{-# INLINE writeMatrix #-}
171writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
172writeMatrix = safeIndexM unsafeWriteMatrix
173
174newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
175newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c
176
177{-# NOINLINE newMatrix #-}
178newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
179newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
deleted file mode 100644
index b5a4318..0000000
--- a/lib/Data/Packed/Vector.hs
+++ /dev/null
@@ -1,96 +0,0 @@
1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE CPP #-}
3-----------------------------------------------------------------------------
4-- |
5-- Module : Data.Packed.Vector
6-- Copyright : (c) Alberto Ruiz 2007-10
7-- License : GPL
8--
9-- Maintainer : Alberto Ruiz <aruiz@um.es>
10-- Stability : provisional
11--
12-- 1D arrays suitable for numeric computations using external libraries.
13--
14-- This module provides basic functions for manipulation of structure.
15--
16-----------------------------------------------------------------------------
17{-# OPTIONS_HADDOCK hide #-}
18
19module Data.Packed.Vector (
20 Vector,
21 fromList, (|>), toList, buildVector,
22 dim, (@>),
23 subVector, takesV, vjoin, join,
24 mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
25 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
26 foldLoop, foldVector, foldVectorG, foldVectorWithIndex
27) where
28
29import Data.Packed.Internal.Vector
30import Foreign.Storable
31
32-------------------------------------------------------------------
33
34#ifdef BINARY
35
36import Data.Binary
37import Control.Monad(replicateM)
38
39-- a 64K cache, with a Double taking 13 bytes in Bytestring,
40-- implies a chunk size of 5041
41chunk :: Int
42chunk = 5000
43
44chunks :: Int -> [Int]
45chunks d = let c = d `div` chunk
46 m = d `mod` chunk
47 in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)
48
49putVector v = do
50 let d = dim v
51 mapM_ (\i -> put $ v @> i) [0..(d-1)]
52
53getVector d = do
54 xs <- replicateM d get
55 return $! fromList xs
56
57instance (Binary a, Storable a) => Binary (Vector a) where
58 put v = do
59 let d = dim v
60 put d
61 mapM_ putVector $! takesV (chunks d) v
62 get = do
63 d <- get
64 vs <- mapM getVector $ chunks d
65 return $! vjoin vs
66
67#endif
68
69-------------------------------------------------------------------
70
71{- | creates a Vector of the specified length using the supplied function to
72 to map the index to the value at that index.
73
74@> buildVector 4 fromIntegral
754 |> [0.0,1.0,2.0,3.0]@
76
77-}
78buildVector :: Storable a => Int -> (Int -> a) -> Vector a
79buildVector len f =
80 fromList $ map f [0 .. (len - 1)]
81
82
83-- | zip for Vectors
84zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b)
85zipVector = zipVectorWith (,)
86
87-- | unzip for Vectors
88unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b)
89unzipVector = unzipVectorWith id
90
91-------------------------------------------------------------------
92
93{-# DEPRECATED join "use vjoin or Data.Vector.concat" #-}
94join :: Storable t => [Vector t] -> Vector t
95join = vjoin
96
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 12:44:18 +0000