diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:16:42 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:16:42 +0200 |
commit | 5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a (patch) | |
tree | c6cdc40bc0121fd87a5137d8b7cb4d9f064cab47 /packages/base/src/Data/Packed/Internal | |
parent | d558f5165e7e7f4daffadae1197e53618727971d (diff) |
Data.Packed -> base (I)
Diffstat (limited to 'packages/base/src/Data/Packed/Internal')
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Common.hs | 160 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Matrix.hs | 422 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Signatures.hs | 70 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Vector.hs | 471 |
4 files changed, 1123 insertions, 0 deletions
diff --git a/packages/base/src/Data/Packed/Internal/Common.hs b/packages/base/src/Data/Packed/Internal/Common.hs new file mode 100644 index 0000000..615bbdf --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Common.hs | |||
@@ -0,0 +1,160 @@ | |||
1 | {-# LANGUAGE CPP #-} | ||
2 | -- | | ||
3 | -- Module : Data.Packed.Internal.Common | ||
4 | -- Copyright : (c) Alberto Ruiz 2007 | ||
5 | -- License : BSD3 | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | -- | ||
10 | -- Development utilities. | ||
11 | -- | ||
12 | |||
13 | |||
14 | module Data.Packed.Internal.Common( | ||
15 | Adapt, | ||
16 | app1, app2, app3, app4, | ||
17 | app5, app6, app7, app8, app9, app10, | ||
18 | (//), check, mbCatch, | ||
19 | splitEvery, common, compatdim, | ||
20 | fi, | ||
21 | table, | ||
22 | finit | ||
23 | ) where | ||
24 | |||
25 | import Control.Monad(when) | ||
26 | import Foreign.C.Types | ||
27 | import Foreign.Storable.Complex() | ||
28 | import Data.List(transpose,intersperse) | ||
29 | import Control.Exception as E | ||
30 | |||
31 | -- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@ | ||
32 | splitEvery :: Int -> [a] -> [[a]] | ||
33 | splitEvery _ [] = [] | ||
34 | splitEvery k l = take k l : splitEvery k (drop k l) | ||
35 | |||
36 | -- | obtains the common value of a property of a list | ||
37 | common :: (Eq a) => (b->a) -> [b] -> Maybe a | ||
38 | common f = commonval . map f where | ||
39 | commonval :: (Eq a) => [a] -> Maybe a | ||
40 | commonval [] = Nothing | ||
41 | commonval [a] = Just a | ||
42 | commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing | ||
43 | |||
44 | -- | common value with \"adaptable\" 1 | ||
45 | compatdim :: [Int] -> Maybe Int | ||
46 | compatdim [] = Nothing | ||
47 | compatdim [a] = Just a | ||
48 | compatdim (a:b:xs) | ||
49 | | a==b = compatdim (b:xs) | ||
50 | | a==1 = compatdim (b:xs) | ||
51 | | b==1 = compatdim (a:xs) | ||
52 | | otherwise = Nothing | ||
53 | |||
54 | -- | Formatting tool | ||
55 | table :: String -> [[String]] -> String | ||
56 | table sep as = unlines . map unwords' $ transpose mtp where | ||
57 | mt = transpose as | ||
58 | longs = map (maximum . map length) mt | ||
59 | mtp = zipWith (\a b -> map (pad a) b) longs mt | ||
60 | pad n str = replicate (n - length str) ' ' ++ str | ||
61 | unwords' = concat . intersperse sep | ||
62 | |||
63 | -- | postfix function application (@flip ($)@) | ||
64 | (//) :: x -> (x -> y) -> y | ||
65 | infixl 0 // | ||
66 | (//) = flip ($) | ||
67 | |||
68 | -- | specialized fromIntegral | ||
69 | fi :: Int -> CInt | ||
70 | fi = fromIntegral | ||
71 | |||
72 | -- hmm.. | ||
73 | ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f | ||
74 | ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f | ||
75 | ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f | ||
76 | ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f | ||
77 | 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 | ||
78 | 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 | ||
79 | 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 | ||
80 | 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 | ||
81 | 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 | ||
82 | |||
83 | type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() | ||
84 | |||
85 | type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO() | ||
86 | type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2 | ||
87 | type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3 | ||
88 | type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4 | ||
89 | 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 | ||
90 | 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 | ||
91 | 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 | ||
92 | 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 | ||
93 | 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 | ||
94 | 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 | ||
95 | |||
96 | app1 :: f -> Adapt1 f t1 | ||
97 | app2 :: f -> Adapt2 f t1 r1 t2 | ||
98 | app3 :: f -> Adapt3 f t1 r1 t2 r2 t3 | ||
99 | app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4 | ||
100 | app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 | ||
101 | app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 | ||
102 | app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 | ||
103 | app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 | ||
104 | app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 | ||
105 | app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 | ||
106 | |||
107 | app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s | ||
108 | app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s | ||
109 | app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $ | ||
110 | \a1 a2 a3 -> f // a1 // a2 // a3 // check s | ||
111 | app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $ | ||
112 | \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s | ||
113 | app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $ | ||
114 | \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s | ||
115 | 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 $ | ||
116 | \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s | ||
117 | 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 $ | ||
118 | \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s | ||
119 | 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 $ | ||
120 | \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s | ||
121 | 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 $ | ||
122 | \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s | ||
123 | 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 $ | ||
124 | \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s | ||
125 | |||
126 | |||
127 | |||
128 | -- GSL error codes are <= 1024 | ||
129 | -- | error codes for the auxiliary functions required by the wrappers | ||
130 | errorCode :: CInt -> String | ||
131 | errorCode 2000 = "bad size" | ||
132 | errorCode 2001 = "bad function code" | ||
133 | errorCode 2002 = "memory problem" | ||
134 | errorCode 2003 = "bad file" | ||
135 | errorCode 2004 = "singular" | ||
136 | errorCode 2005 = "didn't converge" | ||
137 | errorCode 2006 = "the input matrix is not positive definite" | ||
138 | errorCode 2007 = "not yet supported in this OS" | ||
139 | errorCode n = "code "++show n | ||
140 | |||
141 | |||
142 | -- | clear the fpu | ||
143 | foreign import ccall unsafe "asm_finit" finit :: IO () | ||
144 | |||
145 | -- | check the error code | ||
146 | check :: String -> IO CInt -> IO () | ||
147 | check msg f = do | ||
148 | #if FINIT | ||
149 | finit | ||
150 | #endif | ||
151 | err <- f | ||
152 | when (err/=0) $ error (msg++": "++errorCode err) | ||
153 | return () | ||
154 | |||
155 | -- | Error capture and conversion to Maybe | ||
156 | mbCatch :: IO x -> IO (Maybe x) | ||
157 | mbCatch act = E.catch (Just `fmap` act) f | ||
158 | where f :: SomeException -> IO (Maybe x) | ||
159 | f _ = return Nothing | ||
160 | |||
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs new file mode 100644 index 0000000..9b831cc --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Matrix.hs | |||
@@ -0,0 +1,422 @@ | |||
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 : BSD3 | ||
10 | -- Maintainer : Alberto Ruiz | ||
11 | -- Stability : provisional | ||
12 | -- | ||
13 | -- Internal matrix representation | ||
14 | -- | ||
15 | |||
16 | module Data.Packed.Internal.Matrix( | ||
17 | Matrix(..), rows, cols, cdat, fdat, | ||
18 | MatrixOrder(..), orderOf, | ||
19 | createMatrix, mat, | ||
20 | cmat, fmat, | ||
21 | toLists, flatten, reshape, | ||
22 | Element(..), | ||
23 | trans, | ||
24 | fromRows, toRows, fromColumns, toColumns, | ||
25 | matrixFromVector, | ||
26 | subMatrix, | ||
27 | liftMatrix, liftMatrix2, | ||
28 | (@@>), atM', | ||
29 | singleton, | ||
30 | emptyM, | ||
31 | size, shSize, conformVs, conformMs, conformVTo, conformMTo | ||
32 | ) where | ||
33 | |||
34 | import Data.Packed.Internal.Common | ||
35 | import Data.Packed.Internal.Signatures | ||
36 | import Data.Packed.Internal.Vector | ||
37 | |||
38 | import Foreign.Marshal.Alloc(alloca, free) | ||
39 | import Foreign.Marshal.Array(newArray) | ||
40 | import Foreign.Ptr(Ptr, castPtr) | ||
41 | import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf) | ||
42 | import Data.Complex(Complex) | ||
43 | import Foreign.C.Types | ||
44 | import System.IO.Unsafe(unsafePerformIO) | ||
45 | import Control.DeepSeq | ||
46 | |||
47 | ----------------------------------------------------------------- | ||
48 | |||
49 | {- Design considerations for the Matrix Type | ||
50 | ----------------------------------------- | ||
51 | |||
52 | - we must easily handle both row major and column major order, | ||
53 | for bindings to LAPACK and GSL/C | ||
54 | |||
55 | - we'd like to simplify redundant matrix transposes: | ||
56 | - Some of them arise from the order requirements of some functions | ||
57 | - some functions (matrix product) admit transposed arguments | ||
58 | |||
59 | - maybe we don't really need this kind of simplification: | ||
60 | - more complex code | ||
61 | - some computational overhead | ||
62 | - only appreciable gain in code with a lot of redundant transpositions | ||
63 | and cheap matrix computations | ||
64 | |||
65 | - we could carry both the matrix and its (lazily computed) transpose. | ||
66 | This may save some transpositions, but it is necessary to keep track of the | ||
67 | data which is actually computed to be used by functions like the matrix product | ||
68 | which admit both orders. | ||
69 | |||
70 | - but if we need the transposed data and it is not in the structure, we must make | ||
71 | sure that we touch the same foreignptr that is used in the computation. | ||
72 | |||
73 | - a reasonable solution is using two constructors for a matrix. Transposition just | ||
74 | "flips" the constructor. Actual data transposition is not done if followed by a | ||
75 | matrix product or another transpose. | ||
76 | |||
77 | -} | ||
78 | |||
79 | data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) | ||
80 | |||
81 | transOrder RowMajor = ColumnMajor | ||
82 | transOrder ColumnMajor = RowMajor | ||
83 | {- | Matrix representation suitable for GSL and LAPACK computations. | ||
84 | |||
85 | The elements are stored in a continuous memory array. | ||
86 | |||
87 | -} | ||
88 | |||
89 | data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int | ||
90 | , icols :: {-# UNPACK #-} !Int | ||
91 | , xdat :: {-# UNPACK #-} !(Vector t) | ||
92 | , order :: !MatrixOrder } | ||
93 | -- RowMajor: preferred by C, fdat may require a transposition | ||
94 | -- ColumnMajor: preferred by LAPACK, cdat may require a transposition | ||
95 | |||
96 | cdat = xdat | ||
97 | fdat = xdat | ||
98 | |||
99 | rows :: Matrix t -> Int | ||
100 | rows = irows | ||
101 | |||
102 | cols :: Matrix t -> Int | ||
103 | cols = icols | ||
104 | |||
105 | orderOf :: Matrix t -> MatrixOrder | ||
106 | orderOf = order | ||
107 | |||
108 | |||
109 | -- | Matrix transpose. | ||
110 | trans :: Matrix t -> Matrix t | ||
111 | trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o} | ||
112 | |||
113 | cmat :: (Element t) => Matrix t -> Matrix t | ||
114 | cmat m@Matrix{order = RowMajor} = m | ||
115 | cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor} | ||
116 | |||
117 | fmat :: (Element t) => Matrix t -> Matrix t | ||
118 | fmat m@Matrix{order = ColumnMajor} = m | ||
119 | fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor} | ||
120 | |||
121 | -- C-Haskell matrix adapter | ||
122 | -- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r | ||
123 | |||
124 | mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b | ||
125 | mat a f = | ||
126 | unsafeWith (xdat a) $ \p -> do | ||
127 | let m g = do | ||
128 | g (fi (rows a)) (fi (cols a)) p | ||
129 | f m | ||
130 | |||
131 | {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. | ||
132 | |||
133 | >>> flatten (ident 3) | ||
134 | fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] | ||
135 | |||
136 | -} | ||
137 | flatten :: Element t => Matrix t -> Vector t | ||
138 | flatten = xdat . cmat | ||
139 | |||
140 | {- | ||
141 | type Mt t s = Int -> Int -> Ptr t -> s | ||
142 | |||
143 | infixr 6 ::> | ||
144 | type t ::> s = Mt t s | ||
145 | -} | ||
146 | |||
147 | -- | the inverse of 'Data.Packed.Matrix.fromLists' | ||
148 | toLists :: (Element t) => Matrix t -> [[t]] | ||
149 | toLists m = splitEvery (cols m) . toList . flatten $ m | ||
150 | |||
151 | -- | Create a matrix from a list of vectors. | ||
152 | -- All vectors must have the same dimension, | ||
153 | -- or dimension 1, which is are automatically expanded. | ||
154 | fromRows :: Element t => [Vector t] -> Matrix t | ||
155 | fromRows [] = emptyM 0 0 | ||
156 | fromRows vs = case compatdim (map dim vs) of | ||
157 | Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs) | ||
158 | Just 0 -> emptyM r 0 | ||
159 | Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs | ||
160 | where | ||
161 | r = length vs | ||
162 | adapt c v | ||
163 | | c == 0 = fromList[] | ||
164 | | dim v == c = v | ||
165 | | otherwise = constantD (v@>0) c | ||
166 | |||
167 | -- | extracts the rows of a matrix as a list of vectors | ||
168 | toRows :: Element t => Matrix t -> [Vector t] | ||
169 | toRows m | ||
170 | | c == 0 = replicate r (fromList[]) | ||
171 | | otherwise = toRows' 0 | ||
172 | where | ||
173 | v = flatten m | ||
174 | r = rows m | ||
175 | c = cols m | ||
176 | toRows' k | k == r*c = [] | ||
177 | | otherwise = subVector k c v : toRows' (k+c) | ||
178 | |||
179 | -- | Creates a matrix from a list of vectors, as columns | ||
180 | fromColumns :: Element t => [Vector t] -> Matrix t | ||
181 | fromColumns m = trans . fromRows $ m | ||
182 | |||
183 | -- | Creates a list of vectors from the columns of a matrix | ||
184 | toColumns :: Element t => Matrix t -> [Vector t] | ||
185 | toColumns m = toRows . trans $ m | ||
186 | |||
187 | -- | Reads a matrix position. | ||
188 | (@@>) :: Storable t => Matrix t -> (Int,Int) -> t | ||
189 | infixl 9 @@> | ||
190 | m@Matrix {irows = r, icols = c} @@> (i,j) | ||
191 | | safe = if i<0 || i>=r || j<0 || j>=c | ||
192 | then error "matrix indexing out of range" | ||
193 | else atM' m i j | ||
194 | | otherwise = atM' m i j | ||
195 | {-# INLINE (@@>) #-} | ||
196 | |||
197 | -- Unsafe matrix access without range checking | ||
198 | atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j) | ||
199 | atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i) | ||
200 | {-# INLINE atM' #-} | ||
201 | |||
202 | ------------------------------------------------------------------ | ||
203 | |||
204 | matrixFromVector o r c v | ||
205 | | r * c == dim v = m | ||
206 | | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m | ||
207 | where | ||
208 | m = Matrix { irows = r, icols = c, xdat = v, order = o } | ||
209 | |||
210 | -- allocates memory for a new matrix | ||
211 | createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a) | ||
212 | createMatrix ord r c = do | ||
213 | p <- createVector (r*c) | ||
214 | return (matrixFromVector ord r c p) | ||
215 | |||
216 | {- | 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@ | ||
217 | where r is the desired number of rows.) | ||
218 | |||
219 | >>> reshape 4 (fromList [1..12]) | ||
220 | (3><4) | ||
221 | [ 1.0, 2.0, 3.0, 4.0 | ||
222 | , 5.0, 6.0, 7.0, 8.0 | ||
223 | , 9.0, 10.0, 11.0, 12.0 ] | ||
224 | |||
225 | -} | ||
226 | reshape :: Storable t => Int -> Vector t -> Matrix t | ||
227 | reshape 0 v = matrixFromVector RowMajor 0 0 v | ||
228 | reshape c v = matrixFromVector RowMajor (dim v `div` c) c v | ||
229 | |||
230 | singleton x = reshape 1 (fromList [x]) | ||
231 | |||
232 | -- | application of a vector function on the flattened matrix elements | ||
233 | liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b | ||
234 | liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d) | ||
235 | |||
236 | -- | application of a vector function on the flattened matrices elements | ||
237 | liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t | ||
238 | liftMatrix2 f m1 m2 | ||
239 | | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2" | ||
240 | | otherwise = case orderOf m1 of | ||
241 | RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2)) | ||
242 | ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2)) | ||
243 | |||
244 | |||
245 | compat :: Matrix a -> Matrix b -> Bool | ||
246 | compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 | ||
247 | |||
248 | ------------------------------------------------------------------ | ||
249 | |||
250 | {- | Supported matrix elements. | ||
251 | |||
252 | This class provides optimized internal | ||
253 | operations for selected element types. | ||
254 | It provides unoptimised defaults for any 'Storable' type, | ||
255 | so you can create instances simply as: | ||
256 | @instance Element Foo@. | ||
257 | -} | ||
258 | class (Storable a) => Element a where | ||
259 | subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position | ||
260 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | ||
261 | -> Matrix a -> Matrix a | ||
262 | subMatrixD = subMatrix' | ||
263 | transdata :: Int -> Vector a -> Int -> Vector a | ||
264 | transdata = transdataP -- transdata' | ||
265 | constantD :: a -> Int -> Vector a | ||
266 | constantD = constantP -- constant' | ||
267 | |||
268 | |||
269 | instance Element Float where | ||
270 | transdata = transdataAux ctransF | ||
271 | constantD = constantAux cconstantF | ||
272 | |||
273 | instance Element Double where | ||
274 | transdata = transdataAux ctransR | ||
275 | constantD = constantAux cconstantR | ||
276 | |||
277 | instance Element (Complex Float) where | ||
278 | transdata = transdataAux ctransQ | ||
279 | constantD = constantAux cconstantQ | ||
280 | |||
281 | instance Element (Complex Double) where | ||
282 | transdata = transdataAux ctransC | ||
283 | constantD = constantAux cconstantC | ||
284 | |||
285 | ------------------------------------------------------------------- | ||
286 | |||
287 | transdataAux fun c1 d c2 = | ||
288 | if noneed | ||
289 | then d | ||
290 | else unsafePerformIO $ do | ||
291 | v <- createVector (dim d) | ||
292 | unsafeWith d $ \pd -> | ||
293 | unsafeWith v $ \pv -> | ||
294 | fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux" | ||
295 | return v | ||
296 | where r1 = dim d `div` c1 | ||
297 | r2 = dim d `div` c2 | ||
298 | noneed = dim d == 0 || r1 == 1 || c1 == 1 | ||
299 | |||
300 | transdataP :: Storable a => Int -> Vector a -> Int -> Vector a | ||
301 | transdataP c1 d c2 = | ||
302 | if noneed | ||
303 | then d | ||
304 | else unsafePerformIO $ do | ||
305 | v <- createVector (dim d) | ||
306 | unsafeWith d $ \pd -> | ||
307 | unsafeWith v $ \pv -> | ||
308 | ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP" | ||
309 | return v | ||
310 | where r1 = dim d `div` c1 | ||
311 | r2 = dim d `div` c2 | ||
312 | sz = sizeOf (d @> 0) | ||
313 | noneed = dim d == 0 || r1 == 1 || c1 == 1 | ||
314 | |||
315 | foreign import ccall unsafe "transF" ctransF :: TFMFM | ||
316 | foreign import ccall unsafe "transR" ctransR :: TMM | ||
317 | foreign import ccall unsafe "transQ" ctransQ :: TQMQM | ||
318 | foreign import ccall unsafe "transC" ctransC :: TCMCM | ||
319 | foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt | ||
320 | |||
321 | ---------------------------------------------------------------------- | ||
322 | |||
323 | constantAux fun x n = unsafePerformIO $ do | ||
324 | v <- createVector n | ||
325 | px <- newArray [x] | ||
326 | app1 (fun px) vec v "constantAux" | ||
327 | free px | ||
328 | return v | ||
329 | |||
330 | foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF | ||
331 | |||
332 | foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV | ||
333 | |||
334 | foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV | ||
335 | |||
336 | foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV | ||
337 | |||
338 | constantP :: Storable a => a -> Int -> Vector a | ||
339 | constantP a n = unsafePerformIO $ do | ||
340 | let sz = sizeOf a | ||
341 | v <- createVector n | ||
342 | unsafeWith v $ \p -> do | ||
343 | alloca $ \k -> do | ||
344 | poke k a | ||
345 | cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP" | ||
346 | return v | ||
347 | foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt | ||
348 | |||
349 | ---------------------------------------------------------------------- | ||
350 | |||
351 | -- | Extracts a submatrix from a matrix. | ||
352 | subMatrix :: Element a | ||
353 | => (Int,Int) -- ^ (r0,c0) starting position | ||
354 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | ||
355 | -> Matrix a -- ^ input matrix | ||
356 | -> Matrix a -- ^ result | ||
357 | subMatrix (r0,c0) (rt,ct) m | ||
358 | | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) && | ||
359 | 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m | ||
360 | | otherwise = error $ "wrong subMatrix "++ | ||
361 | show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m) | ||
362 | |||
363 | subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do | ||
364 | w <- createVector (rt*ct) | ||
365 | unsafeWith v $ \p -> | ||
366 | unsafeWith w $ \q -> do | ||
367 | let go (-1) _ = return () | ||
368 | go !i (-1) = go (i-1) (ct-1) | ||
369 | go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0) | ||
370 | pokeElemOff q (i*ct+j) x | ||
371 | go i (j-1) | ||
372 | go (rt-1) (ct-1) | ||
373 | return w | ||
374 | |||
375 | subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor | ||
376 | subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m) | ||
377 | |||
378 | -------------------------------------------------------------------------- | ||
379 | |||
380 | maxZ xs = if minimum xs == 0 then 0 else maximum xs | ||
381 | |||
382 | conformMs ms = map (conformMTo (r,c)) ms | ||
383 | where | ||
384 | r = maxZ (map rows ms) | ||
385 | c = maxZ (map cols ms) | ||
386 | |||
387 | |||
388 | conformVs vs = map (conformVTo n) vs | ||
389 | where | ||
390 | n = maxZ (map dim vs) | ||
391 | |||
392 | conformMTo (r,c) m | ||
393 | | size m == (r,c) = m | ||
394 | | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c)) | ||
395 | | size m == (r,1) = repCols c m | ||
396 | | size m == (1,c) = repRows r m | ||
397 | | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")" | ||
398 | |||
399 | conformVTo n v | ||
400 | | dim v == n = v | ||
401 | | dim v == 1 = constantD (v@>0) n | ||
402 | | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n | ||
403 | |||
404 | repRows n x = fromRows (replicate n (flatten x)) | ||
405 | repCols n x = fromColumns (replicate n (flatten x)) | ||
406 | |||
407 | size m = (rows m, cols m) | ||
408 | |||
409 | shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")" | ||
410 | |||
411 | emptyM r c = matrixFromVector RowMajor r c (fromList[]) | ||
412 | |||
413 | ---------------------------------------------------------------------- | ||
414 | |||
415 | instance (Storable t, NFData t) => NFData (Matrix t) | ||
416 | where | ||
417 | rnf m | d > 0 = rnf (v @> 0) | ||
418 | | otherwise = () | ||
419 | where | ||
420 | d = dim v | ||
421 | v = xdat m | ||
422 | |||
diff --git a/packages/base/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs new file mode 100644 index 0000000..acc3070 --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Signatures.hs | |||
@@ -0,0 +1,70 @@ | |||
1 | -- | | ||
2 | -- Module : Data.Packed.Internal.Signatures | ||
3 | -- Copyright : (c) Alberto Ruiz 2009 | ||
4 | -- License : BSD3 | ||
5 | -- Maintainer : Alberto Ruiz | ||
6 | -- Stability : provisional | ||
7 | -- | ||
8 | -- Signatures of the C functions. | ||
9 | -- | ||
10 | |||
11 | |||
12 | module Data.Packed.Internal.Signatures where | ||
13 | |||
14 | import Foreign.Ptr(Ptr) | ||
15 | import Data.Complex(Complex) | ||
16 | import Foreign.C.Types(CInt) | ||
17 | |||
18 | type PF = Ptr Float -- | ||
19 | type PD = Ptr Double -- | ||
20 | type PQ = Ptr (Complex Float) -- | ||
21 | type PC = Ptr (Complex Double) -- | ||
22 | type TF = CInt -> PF -> IO CInt -- | ||
23 | type TFF = CInt -> PF -> TF -- | ||
24 | type TFV = CInt -> PF -> TV -- | ||
25 | type TVF = CInt -> PD -> TF -- | ||
26 | type TFFF = CInt -> PF -> TFF -- | ||
27 | type TV = CInt -> PD -> IO CInt -- | ||
28 | type TVV = CInt -> PD -> TV -- | ||
29 | type TVVV = CInt -> PD -> TVV -- | ||
30 | type TFM = CInt -> CInt -> PF -> IO CInt -- | ||
31 | type TFMFM = CInt -> CInt -> PF -> TFM -- | ||
32 | type TFMFMFM = CInt -> CInt -> PF -> TFMFM -- | ||
33 | type TM = CInt -> CInt -> PD -> IO CInt -- | ||
34 | type TMM = CInt -> CInt -> PD -> TM -- | ||
35 | type TVMM = CInt -> PD -> TMM -- | ||
36 | type TMVMM = CInt -> CInt -> PD -> TVMM -- | ||
37 | type TMMM = CInt -> CInt -> PD -> TMM -- | ||
38 | type TVM = CInt -> PD -> TM -- | ||
39 | type TVVM = CInt -> PD -> TVM -- | ||
40 | type TMV = CInt -> CInt -> PD -> TV -- | ||
41 | type TMMV = CInt -> CInt -> PD -> TMV -- | ||
42 | type TMVM = CInt -> CInt -> PD -> TVM -- | ||
43 | type TMMVM = CInt -> CInt -> PD -> TMVM -- | ||
44 | type TCM = CInt -> CInt -> PC -> IO CInt -- | ||
45 | type TCVCM = CInt -> PC -> TCM -- | ||
46 | type TCMCVCM = CInt -> CInt -> PC -> TCVCM -- | ||
47 | type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM -- | ||
48 | type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM -- | ||
49 | type TCMCM = CInt -> CInt -> PC -> TCM -- | ||
50 | type TVCM = CInt -> PD -> TCM -- | ||
51 | type TCMVCM = CInt -> CInt -> PC -> TVCM -- | ||
52 | type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM -- | ||
53 | type TCMCMCM = CInt -> CInt -> PC -> TCMCM -- | ||
54 | type TCV = CInt -> PC -> IO CInt -- | ||
55 | type TCVCV = CInt -> PC -> TCV -- | ||
56 | type TCVCVCV = CInt -> PC -> TCVCV -- | ||
57 | type TCVV = CInt -> PC -> TV -- | ||
58 | type TQV = CInt -> PQ -> IO CInt -- | ||
59 | type TQVQV = CInt -> PQ -> TQV -- | ||
60 | type TQVQVQV = CInt -> PQ -> TQVQV -- | ||
61 | type TQVF = CInt -> PQ -> TF -- | ||
62 | type TQM = CInt -> CInt -> PQ -> IO CInt -- | ||
63 | type TQMQM = CInt -> CInt -> PQ -> TQM -- | ||
64 | type TQMQMQM = CInt -> CInt -> PQ -> TQMQM -- | ||
65 | type TCMCV = CInt -> CInt -> PC -> TCV -- | ||
66 | type TVCV = CInt -> PD -> TCV -- | ||
67 | type TCVM = CInt -> PC -> TM -- | ||
68 | type TMCVM = CInt -> CInt -> PD -> TCVM -- | ||
69 | type TMMCVM = CInt -> CInt -> PD -> TMCVM -- | ||
70 | |||
diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs new file mode 100644 index 0000000..d0bc143 --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Vector.hs | |||
@@ -0,0 +1,471 @@ | |||
1 | {-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} | ||
2 | -- | | ||
3 | -- Module : Data.Packed.Internal.Vector | ||
4 | -- Copyright : (c) Alberto Ruiz 2007 | ||
5 | -- License : BSD3 | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | -- Vector implementation | ||
10 | -- | ||
11 | -------------------------------------------------------------------------------- | ||
12 | |||
13 | module Data.Packed.Internal.Vector ( | ||
14 | Vector, dim, | ||
15 | fromList, toList, (|>), | ||
16 | vjoin, (@>), safe, at, at', subVector, takesV, | ||
17 | mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith, | ||
18 | mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, | ||
19 | foldVector, foldVectorG, foldLoop, foldVectorWithIndex, | ||
20 | createVector, vec, | ||
21 | asComplex, asReal, float2DoubleV, double2FloatV, | ||
22 | stepF, stepD, condF, condD, | ||
23 | conjugateQ, conjugateC, | ||
24 | cloneVector, | ||
25 | unsafeToForeignPtr, | ||
26 | unsafeFromForeignPtr, | ||
27 | unsafeWith | ||
28 | ) where | ||
29 | |||
30 | import Data.Packed.Internal.Common | ||
31 | import Data.Packed.Internal.Signatures | ||
32 | import Foreign.Marshal.Array(peekArray, copyArray, advancePtr) | ||
33 | import Foreign.ForeignPtr(ForeignPtr, castForeignPtr) | ||
34 | import Foreign.Ptr(Ptr) | ||
35 | import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf) | ||
36 | import Foreign.C.Types | ||
37 | import Data.Complex | ||
38 | import Control.Monad(when) | ||
39 | import System.IO.Unsafe(unsafePerformIO) | ||
40 | |||
41 | #if __GLASGOW_HASKELL__ >= 605 | ||
42 | import GHC.ForeignPtr (mallocPlainForeignPtrBytes) | ||
43 | #else | ||
44 | import Foreign.ForeignPtr (mallocForeignPtrBytes) | ||
45 | #endif | ||
46 | |||
47 | import GHC.Base | ||
48 | #if __GLASGOW_HASKELL__ < 612 | ||
49 | import GHC.IOBase hiding (liftIO) | ||
50 | #endif | ||
51 | |||
52 | import qualified Data.Vector.Storable as Vector | ||
53 | import Data.Vector.Storable(Vector, | ||
54 | fromList, | ||
55 | unsafeToForeignPtr, | ||
56 | unsafeFromForeignPtr, | ||
57 | unsafeWith) | ||
58 | |||
59 | |||
60 | -- | Number of elements | ||
61 | dim :: (Storable t) => Vector t -> Int | ||
62 | dim = Vector.length | ||
63 | |||
64 | |||
65 | -- C-Haskell vector adapter | ||
66 | -- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r | ||
67 | vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b | ||
68 | vec x f = unsafeWith x $ \p -> do | ||
69 | let v g = do | ||
70 | g (fi $ dim x) p | ||
71 | f v | ||
72 | {-# INLINE vec #-} | ||
73 | |||
74 | |||
75 | -- allocates memory for a new vector | ||
76 | createVector :: Storable a => Int -> IO (Vector a) | ||
77 | createVector n = do | ||
78 | when (n < 0) $ error ("trying to createVector of negative dim: "++show n) | ||
79 | fp <- doMalloc undefined | ||
80 | return $ unsafeFromForeignPtr fp 0 n | ||
81 | where | ||
82 | -- | ||
83 | -- Use the much cheaper Haskell heap allocated storage | ||
84 | -- for foreign pointer space we control | ||
85 | -- | ||
86 | doMalloc :: Storable b => b -> IO (ForeignPtr b) | ||
87 | doMalloc dummy = do | ||
88 | #if __GLASGOW_HASKELL__ >= 605 | ||
89 | mallocPlainForeignPtrBytes (n * sizeOf dummy) | ||
90 | #else | ||
91 | mallocForeignPtrBytes (n * sizeOf dummy) | ||
92 | #endif | ||
93 | |||
94 | {- | creates a Vector from a list: | ||
95 | |||
96 | @> fromList [2,3,5,7] | ||
97 | 4 |> [2.0,3.0,5.0,7.0]@ | ||
98 | |||
99 | -} | ||
100 | |||
101 | safeRead v = inlinePerformIO . unsafeWith v | ||
102 | {-# INLINE safeRead #-} | ||
103 | |||
104 | inlinePerformIO :: IO a -> a | ||
105 | inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r | ||
106 | {-# INLINE inlinePerformIO #-} | ||
107 | |||
108 | {- | extracts the Vector elements to a list | ||
109 | |||
110 | >>> toList (linspace 5 (1,10)) | ||
111 | [1.0,3.25,5.5,7.75,10.0] | ||
112 | |||
113 | -} | ||
114 | toList :: Storable a => Vector a -> [a] | ||
115 | toList v = safeRead v $ peekArray (dim v) | ||
116 | |||
117 | {- | Create a vector from a list of elements and explicit dimension. The input | ||
118 | list is explicitly truncated if it is too long, so it may safely | ||
119 | be used, for instance, with infinite lists. | ||
120 | |||
121 | >>> 5 |> [1..] | ||
122 | fromList [1.0,2.0,3.0,4.0,5.0] | ||
123 | |||
124 | -} | ||
125 | (|>) :: (Storable a) => Int -> [a] -> Vector a | ||
126 | infixl 9 |> | ||
127 | n |> l = if length l' == n | ||
128 | then fromList l' | ||
129 | else error "list too short for |>" | ||
130 | where l' = take n l | ||
131 | |||
132 | |||
133 | -- | access to Vector elements without range checking | ||
134 | at' :: Storable a => Vector a -> Int -> a | ||
135 | at' v n = safeRead v $ flip peekElemOff n | ||
136 | {-# INLINE at' #-} | ||
137 | |||
138 | -- | ||
139 | -- turn off bounds checking with -funsafe at configure time. | ||
140 | -- ghc will optimise away the salways true case at compile time. | ||
141 | -- | ||
142 | #if defined(UNSAFE) | ||
143 | safe :: Bool | ||
144 | safe = False | ||
145 | #else | ||
146 | safe = True | ||
147 | #endif | ||
148 | |||
149 | -- | access to Vector elements with range checking. | ||
150 | at :: Storable a => Vector a -> Int -> a | ||
151 | at v n | ||
152 | | safe = if n >= 0 && n < dim v | ||
153 | then at' v n | ||
154 | else error "vector index out of range" | ||
155 | | otherwise = at' v n | ||
156 | {-# INLINE at #-} | ||
157 | |||
158 | {- | takes a number of consecutive elements from a Vector | ||
159 | |||
160 | >>> subVector 2 3 (fromList [1..10]) | ||
161 | fromList [3.0,4.0,5.0] | ||
162 | |||
163 | -} | ||
164 | subVector :: Storable t => Int -- ^ index of the starting element | ||
165 | -> Int -- ^ number of elements to extract | ||
166 | -> Vector t -- ^ source | ||
167 | -> Vector t -- ^ result | ||
168 | subVector = Vector.slice | ||
169 | |||
170 | |||
171 | {- | Reads a vector position: | ||
172 | |||
173 | >>> fromList [0..9] @> 7 | ||
174 | 7.0 | ||
175 | |||
176 | -} | ||
177 | (@>) :: Storable t => Vector t -> Int -> t | ||
178 | infixl 9 @> | ||
179 | (@>) = at | ||
180 | |||
181 | |||
182 | {- | concatenate a list of vectors | ||
183 | |||
184 | >>> vjoin [fromList [1..5::Double], konst 1 3] | ||
185 | fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] | ||
186 | |||
187 | -} | ||
188 | vjoin :: Storable t => [Vector t] -> Vector t | ||
189 | vjoin [] = fromList [] | ||
190 | vjoin [v] = v | ||
191 | vjoin as = unsafePerformIO $ do | ||
192 | let tot = sum (map dim as) | ||
193 | r <- createVector tot | ||
194 | unsafeWith r $ \ptr -> | ||
195 | joiner as tot ptr | ||
196 | return r | ||
197 | where joiner [] _ _ = return () | ||
198 | joiner (v:cs) _ p = do | ||
199 | let n = dim v | ||
200 | unsafeWith v $ \pb -> copyArray p pb n | ||
201 | joiner cs 0 (advancePtr p n) | ||
202 | |||
203 | |||
204 | {- | Extract consecutive subvectors of the given sizes. | ||
205 | |||
206 | >>> takesV [3,4] (linspace 10 (1,10::Double)) | ||
207 | [fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]] | ||
208 | |||
209 | -} | ||
210 | takesV :: Storable t => [Int] -> Vector t -> [Vector t] | ||
211 | takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w) | ||
212 | | otherwise = go ms w | ||
213 | where go [] _ = [] | ||
214 | go (n:ns) v = subVector 0 n v | ||
215 | : go ns (subVector n (dim v - n) v) | ||
216 | |||
217 | --------------------------------------------------------------- | ||
218 | |||
219 | -- | transforms a complex vector into a real vector with alternating real and imaginary parts | ||
220 | asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a | ||
221 | asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) | ||
222 | where (fp,i,n) = unsafeToForeignPtr v | ||
223 | |||
224 | -- | transforms a real vector into a complex vector with alternating real and imaginary parts | ||
225 | asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a) | ||
226 | asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2) | ||
227 | where (fp,i,n) = unsafeToForeignPtr v | ||
228 | |||
229 | --------------------------------------------------------------- | ||
230 | |||
231 | float2DoubleV :: Vector Float -> Vector Double | ||
232 | float2DoubleV v = unsafePerformIO $ do | ||
233 | r <- createVector (dim v) | ||
234 | app2 c_float2double vec v vec r "float2double" | ||
235 | return r | ||
236 | |||
237 | double2FloatV :: Vector Double -> Vector Float | ||
238 | double2FloatV v = unsafePerformIO $ do | ||
239 | r <- createVector (dim v) | ||
240 | app2 c_double2float vec v vec r "double2float2" | ||
241 | return r | ||
242 | |||
243 | |||
244 | foreign import ccall unsafe "float2double" c_float2double:: TFV | ||
245 | foreign import ccall unsafe "double2float" c_double2float:: TVF | ||
246 | |||
247 | --------------------------------------------------------------- | ||
248 | |||
249 | stepF :: Vector Float -> Vector Float | ||
250 | stepF v = unsafePerformIO $ do | ||
251 | r <- createVector (dim v) | ||
252 | app2 c_stepF vec v vec r "stepF" | ||
253 | return r | ||
254 | |||
255 | stepD :: Vector Double -> Vector Double | ||
256 | stepD v = unsafePerformIO $ do | ||
257 | r <- createVector (dim v) | ||
258 | app2 c_stepD vec v vec r "stepD" | ||
259 | return r | ||
260 | |||
261 | foreign import ccall unsafe "stepF" c_stepF :: TFF | ||
262 | foreign import ccall unsafe "stepD" c_stepD :: TVV | ||
263 | |||
264 | --------------------------------------------------------------- | ||
265 | |||
266 | condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float | ||
267 | condF x y l e g = unsafePerformIO $ do | ||
268 | r <- createVector (dim x) | ||
269 | app6 c_condF vec x vec y vec l vec e vec g vec r "condF" | ||
270 | return r | ||
271 | |||
272 | condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double | ||
273 | condD x y l e g = unsafePerformIO $ do | ||
274 | r <- createVector (dim x) | ||
275 | app6 c_condD vec x vec y vec l vec e vec g vec r "condD" | ||
276 | return r | ||
277 | |||
278 | foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF | ||
279 | foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV | ||
280 | |||
281 | -------------------------------------------------------------------------------- | ||
282 | |||
283 | conjugateAux fun x = unsafePerformIO $ do | ||
284 | v <- createVector (dim x) | ||
285 | app2 fun vec x vec v "conjugateAux" | ||
286 | return v | ||
287 | |||
288 | conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) | ||
289 | conjugateQ = conjugateAux c_conjugateQ | ||
290 | foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV | ||
291 | |||
292 | conjugateC :: Vector (Complex Double) -> Vector (Complex Double) | ||
293 | conjugateC = conjugateAux c_conjugateC | ||
294 | foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV | ||
295 | |||
296 | -------------------------------------------------------------------------------- | ||
297 | |||
298 | cloneVector :: Storable t => Vector t -> IO (Vector t) | ||
299 | cloneVector v = do | ||
300 | let n = dim v | ||
301 | r <- createVector n | ||
302 | let f _ s _ d = copyArray d s n >> return 0 | ||
303 | app2 f vec v vec r "cloneVector" | ||
304 | return r | ||
305 | |||
306 | ------------------------------------------------------------------ | ||
307 | |||
308 | -- | map on Vectors | ||
309 | mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b | ||
310 | mapVector f v = unsafePerformIO $ do | ||
311 | w <- createVector (dim v) | ||
312 | unsafeWith v $ \p -> | ||
313 | unsafeWith w $ \q -> do | ||
314 | let go (-1) = return () | ||
315 | go !k = do x <- peekElemOff p k | ||
316 | pokeElemOff q k (f x) | ||
317 | go (k-1) | ||
318 | go (dim v -1) | ||
319 | return w | ||
320 | {-# INLINE mapVector #-} | ||
321 | |||
322 | -- | zipWith for Vectors | ||
323 | zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c | ||
324 | zipVectorWith f u v = unsafePerformIO $ do | ||
325 | let n = min (dim u) (dim v) | ||
326 | w <- createVector n | ||
327 | unsafeWith u $ \pu -> | ||
328 | unsafeWith v $ \pv -> | ||
329 | unsafeWith w $ \pw -> do | ||
330 | let go (-1) = return () | ||
331 | go !k = do x <- peekElemOff pu k | ||
332 | y <- peekElemOff pv k | ||
333 | pokeElemOff pw k (f x y) | ||
334 | go (k-1) | ||
335 | go (n -1) | ||
336 | return w | ||
337 | {-# INLINE zipVectorWith #-} | ||
338 | |||
339 | -- | unzipWith for Vectors | ||
340 | unzipVectorWith :: (Storable (a,b), Storable c, Storable d) | ||
341 | => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) | ||
342 | unzipVectorWith f u = unsafePerformIO $ do | ||
343 | let n = dim u | ||
344 | v <- createVector n | ||
345 | w <- createVector n | ||
346 | unsafeWith u $ \pu -> | ||
347 | unsafeWith v $ \pv -> | ||
348 | unsafeWith w $ \pw -> do | ||
349 | let go (-1) = return () | ||
350 | go !k = do z <- peekElemOff pu k | ||
351 | let (x,y) = f z | ||
352 | pokeElemOff pv k x | ||
353 | pokeElemOff pw k y | ||
354 | go (k-1) | ||
355 | go (n-1) | ||
356 | return (v,w) | ||
357 | {-# INLINE unzipVectorWith #-} | ||
358 | |||
359 | foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b | ||
360 | foldVector f x v = unsafePerformIO $ | ||
361 | unsafeWith v $ \p -> do | ||
362 | let go (-1) s = return s | ||
363 | go !k !s = do y <- peekElemOff p k | ||
364 | go (k-1::Int) (f y s) | ||
365 | go (dim v -1) x | ||
366 | {-# INLINE foldVector #-} | ||
367 | |||
368 | -- the zero-indexed index is passed to the folding function | ||
369 | foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b | ||
370 | foldVectorWithIndex f x v = unsafePerformIO $ | ||
371 | unsafeWith v $ \p -> do | ||
372 | let go (-1) s = return s | ||
373 | go !k !s = do y <- peekElemOff p k | ||
374 | go (k-1::Int) (f k y s) | ||
375 | go (dim v -1) x | ||
376 | {-# INLINE foldVectorWithIndex #-} | ||
377 | |||
378 | foldLoop f s0 d = go (d - 1) s0 | ||
379 | where | ||
380 | go 0 s = f (0::Int) s | ||
381 | go !j !s = go (j - 1) (f j s) | ||
382 | |||
383 | foldVectorG f s0 v = foldLoop g s0 (dim v) | ||
384 | where g !k !s = f k (at' v) s | ||
385 | {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) | ||
386 | {-# INLINE foldVectorG #-} | ||
387 | |||
388 | ------------------------------------------------------------------- | ||
389 | |||
390 | -- | monadic map over Vectors | ||
391 | -- the monad @m@ must be strict | ||
392 | mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b) | ||
393 | mapVectorM f v = do | ||
394 | w <- return $! unsafePerformIO $! createVector (dim v) | ||
395 | mapVectorM' w 0 (dim v -1) | ||
396 | return w | ||
397 | where mapVectorM' w' !k !t | ||
398 | | k == t = do | ||
399 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
400 | y <- f x | ||
401 | return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
402 | | otherwise = do | ||
403 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
404 | y <- f x | ||
405 | _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
406 | mapVectorM' w' (k+1) t | ||
407 | {-# INLINE mapVectorM #-} | ||
408 | |||
409 | -- | monadic map over Vectors | ||
410 | mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m () | ||
411 | mapVectorM_ f v = do | ||
412 | mapVectorM' 0 (dim v -1) | ||
413 | where mapVectorM' !k !t | ||
414 | | k == t = do | ||
415 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
416 | f x | ||
417 | | otherwise = do | ||
418 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
419 | _ <- f x | ||
420 | mapVectorM' (k+1) t | ||
421 | {-# INLINE mapVectorM_ #-} | ||
422 | |||
423 | -- | monadic map over Vectors with the zero-indexed index passed to the mapping function | ||
424 | -- the monad @m@ must be strict | ||
425 | mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b) | ||
426 | mapVectorWithIndexM f v = do | ||
427 | w <- return $! unsafePerformIO $! createVector (dim v) | ||
428 | mapVectorM' w 0 (dim v -1) | ||
429 | return w | ||
430 | where mapVectorM' w' !k !t | ||
431 | | k == t = do | ||
432 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
433 | y <- f k x | ||
434 | return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
435 | | otherwise = do | ||
436 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
437 | y <- f k x | ||
438 | _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
439 | mapVectorM' w' (k+1) t | ||
440 | {-# INLINE mapVectorWithIndexM #-} | ||
441 | |||
442 | -- | monadic map over Vectors with the zero-indexed index passed to the mapping function | ||
443 | mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m () | ||
444 | mapVectorWithIndexM_ f v = do | ||
445 | mapVectorM' 0 (dim v -1) | ||
446 | where mapVectorM' !k !t | ||
447 | | k == t = do | ||
448 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
449 | f k x | ||
450 | | otherwise = do | ||
451 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
452 | _ <- f k x | ||
453 | mapVectorM' (k+1) t | ||
454 | {-# INLINE mapVectorWithIndexM_ #-} | ||
455 | |||
456 | |||
457 | mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b | ||
458 | --mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b]) | ||
459 | mapVectorWithIndex f v = unsafePerformIO $ do | ||
460 | w <- createVector (dim v) | ||
461 | unsafeWith v $ \p -> | ||
462 | unsafeWith w $ \q -> do | ||
463 | let go (-1) = return () | ||
464 | go !k = do x <- peekElemOff p k | ||
465 | pokeElemOff q k (f k x) | ||
466 | go (k-1) | ||
467 | go (dim v -1) | ||
468 | return w | ||
469 | {-# INLINE mapVectorWithIndex #-} | ||
470 | |||
471 | |||