summaryrefslogtreecommitdiff
path: root/packages/base/src/Data/Packed/Internal
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Data/Packed/Internal')
-rw-r--r--packages/base/src/Data/Packed/Internal/Common.hs160
-rw-r--r--packages/base/src/Data/Packed/Internal/Matrix.hs423
-rw-r--r--packages/base/src/Data/Packed/Internal/Numeric.hs720
-rw-r--r--packages/base/src/Data/Packed/Internal/Signatures.hs70
-rw-r--r--packages/base/src/Data/Packed/Internal/Vector.hs471
5 files changed, 0 insertions, 1844 deletions
diff --git a/packages/base/src/Data/Packed/Internal/Common.hs b/packages/base/src/Data/Packed/Internal/Common.hs
deleted file mode 100644
index 615bbdf..0000000
--- a/packages/base/src/Data/Packed/Internal/Common.hs
+++ /dev/null
@@ -1,160 +0,0 @@
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
14module 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
25import Control.Monad(when)
26import Foreign.C.Types
27import Foreign.Storable.Complex()
28import Data.List(transpose,intersperse)
29import Control.Exception as E
30
31-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
32splitEvery :: Int -> [a] -> [[a]]
33splitEvery _ [] = []
34splitEvery k l = take k l : splitEvery k (drop k l)
35
36-- | obtains the common value of a property of a list
37common :: (Eq a) => (b->a) -> [b] -> Maybe a
38common 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
45compatdim :: [Int] -> Maybe Int
46compatdim [] = Nothing
47compatdim [a] = Just a
48compatdim (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
55table :: String -> [[String]] -> String
56table 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
65infixl 0 //
66(//) = flip ($)
67
68-- | specialized fromIntegral
69fi :: Int -> CInt
70fi = fromIntegral
71
72-- hmm..
73ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
74ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
75ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
76ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
77ww6 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
78ww7 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
79ww8 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
80ww9 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
81ww10 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
83type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
84
85type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
86type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
87type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
88type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
89type 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
90type 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
91type 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
92type 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
93type 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
94type 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
96app1 :: f -> Adapt1 f t1
97app2 :: f -> Adapt2 f t1 r1 t2
98app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
99app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
100app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
101app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
102app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
103app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
104app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
105app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
106
107app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
108app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
109app3 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
111app4 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
113app5 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
115app6 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
117app7 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
119app8 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
121app9 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
123app10 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
130errorCode :: CInt -> String
131errorCode 2000 = "bad size"
132errorCode 2001 = "bad function code"
133errorCode 2002 = "memory problem"
134errorCode 2003 = "bad file"
135errorCode 2004 = "singular"
136errorCode 2005 = "didn't converge"
137errorCode 2006 = "the input matrix is not positive definite"
138errorCode 2007 = "not yet supported in this OS"
139errorCode n = "code "++show n
140
141
142-- | clear the fpu
143foreign import ccall unsafe "asm_finit" finit :: IO ()
144
145-- | check the error code
146check :: String -> IO CInt -> IO ()
147check 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
156mbCatch :: IO x -> IO (Maybe x)
157mbCatch 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
deleted file mode 100644
index 150b978..0000000
--- a/packages/base/src/Data/Packed/Internal/Matrix.hs
+++ /dev/null
@@ -1,423 +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 : BSD3
10-- Maintainer : Alberto Ruiz
11-- Stability : provisional
12--
13-- Internal matrix representation
14--
15
16module 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
34import Data.Packed.Internal.Common
35import Data.Packed.Internal.Signatures
36import Data.Packed.Internal.Vector
37
38import Foreign.Marshal.Alloc(alloca, free)
39import Foreign.Marshal.Array(newArray)
40import Foreign.Ptr(Ptr, castPtr)
41import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
42import Data.Complex(Complex)
43import Foreign.C.Types
44import System.IO.Unsafe(unsafePerformIO)
45import 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
79data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
80
81transOrder RowMajor = ColumnMajor
82transOrder ColumnMajor = RowMajor
83{- | Matrix representation suitable for BLAS\/LAPACK computations.
84
85The elements are stored in a continuous memory array.
86
87-}
88
89data 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
96cdat = xdat
97fdat = xdat
98
99rows :: Matrix t -> Int
100rows = irows
101
102cols :: Matrix t -> Int
103cols = icols
104
105orderOf :: Matrix t -> MatrixOrder
106orderOf = order
107
108
109-- | Matrix transpose.
110trans :: Matrix t -> Matrix t
111trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o}
112
113cmat :: (Element t) => Matrix t -> Matrix t
114cmat m@Matrix{order = RowMajor} = m
115cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor}
116
117fmat :: (Element t) => Matrix t -> Matrix t
118fmat m@Matrix{order = ColumnMajor} = m
119fmat 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
124mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
125mat 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)
134fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
135
136-}
137flatten :: Element t => Matrix t -> Vector t
138flatten = xdat . cmat
139
140{-
141type Mt t s = Int -> Int -> Ptr t -> s
142
143infixr 6 ::>
144type t ::> s = Mt t s
145-}
146
147-- | the inverse of 'Data.Packed.Matrix.fromLists'
148toLists :: (Element t) => Matrix t -> [[t]]
149toLists 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.
154fromRows :: Element t => [Vector t] -> Matrix t
155fromRows [] = emptyM 0 0
156fromRows 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
168toRows :: Element t => Matrix t -> [Vector t]
169toRows 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
180fromColumns :: Element t => [Vector t] -> Matrix t
181fromColumns m = trans . fromRows $ m
182
183-- | Creates a list of vectors from the columns of a matrix
184toColumns :: Element t => Matrix t -> [Vector t]
185toColumns m = toRows . trans $ m
186
187-- | Reads a matrix position.
188(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
189infixl 9 @@>
190m@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
198atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j)
199atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
200{-# INLINE atM' #-}
201
202------------------------------------------------------------------
203
204matrixFromVector 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
211createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
212createMatrix 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@
217where 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-}
226reshape :: Storable t => Int -> Vector t -> Matrix t
227reshape 0 v = matrixFromVector RowMajor 0 0 v
228reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
229
230singleton x = reshape 1 (fromList [x])
231
232-- | application of a vector function on the flattened matrix elements
233liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
234liftMatrix 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
237liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
238liftMatrix2 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
245compat :: Matrix a -> Matrix b -> Bool
246compat 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
257 >instance Element Foo
258-}
259class (Storable a) => Element a where
260 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
261 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
262 -> Matrix a -> Matrix a
263 subMatrixD = subMatrix'
264 transdata :: Int -> Vector a -> Int -> Vector a
265 transdata = transdataP -- transdata'
266 constantD :: a -> Int -> Vector a
267 constantD = constantP -- constant'
268
269
270instance Element Float where
271 transdata = transdataAux ctransF
272 constantD = constantAux cconstantF
273
274instance Element Double where
275 transdata = transdataAux ctransR
276 constantD = constantAux cconstantR
277
278instance Element (Complex Float) where
279 transdata = transdataAux ctransQ
280 constantD = constantAux cconstantQ
281
282instance Element (Complex Double) where
283 transdata = transdataAux ctransC
284 constantD = constantAux cconstantC
285
286-------------------------------------------------------------------
287
288transdataAux fun c1 d c2 =
289 if noneed
290 then d
291 else unsafePerformIO $ do
292 v <- createVector (dim d)
293 unsafeWith d $ \pd ->
294 unsafeWith v $ \pv ->
295 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
296 return v
297 where r1 = dim d `div` c1
298 r2 = dim d `div` c2
299 noneed = dim d == 0 || r1 == 1 || c1 == 1
300
301transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
302transdataP c1 d c2 =
303 if noneed
304 then d
305 else unsafePerformIO $ do
306 v <- createVector (dim d)
307 unsafeWith d $ \pd ->
308 unsafeWith v $ \pv ->
309 ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"
310 return v
311 where r1 = dim d `div` c1
312 r2 = dim d `div` c2
313 sz = sizeOf (d @> 0)
314 noneed = dim d == 0 || r1 == 1 || c1 == 1
315
316foreign import ccall unsafe "transF" ctransF :: TFMFM
317foreign import ccall unsafe "transR" ctransR :: TMM
318foreign import ccall unsafe "transQ" ctransQ :: TQMQM
319foreign import ccall unsafe "transC" ctransC :: TCMCM
320foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
321
322----------------------------------------------------------------------
323
324constantAux fun x n = unsafePerformIO $ do
325 v <- createVector n
326 px <- newArray [x]
327 app1 (fun px) vec v "constantAux"
328 free px
329 return v
330
331foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
332
333foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
334
335foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
336
337foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
338
339constantP :: Storable a => a -> Int -> Vector a
340constantP a n = unsafePerformIO $ do
341 let sz = sizeOf a
342 v <- createVector n
343 unsafeWith v $ \p -> do
344 alloca $ \k -> do
345 poke k a
346 cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
347 return v
348foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
349
350----------------------------------------------------------------------
351
352-- | Extracts a submatrix from a matrix.
353subMatrix :: Element a
354 => (Int,Int) -- ^ (r0,c0) starting position
355 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
356 -> Matrix a -- ^ input matrix
357 -> Matrix a -- ^ result
358subMatrix (r0,c0) (rt,ct) m
359 | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
360 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
361 | otherwise = error $ "wrong subMatrix "++
362 show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
363
364subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
365 w <- createVector (rt*ct)
366 unsafeWith v $ \p ->
367 unsafeWith w $ \q -> do
368 let go (-1) _ = return ()
369 go !i (-1) = go (i-1) (ct-1)
370 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)
371 pokeElemOff q (i*ct+j) x
372 go i (j-1)
373 go (rt-1) (ct-1)
374 return w
375
376subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor
377subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
378
379--------------------------------------------------------------------------
380
381maxZ xs = if minimum xs == 0 then 0 else maximum xs
382
383conformMs ms = map (conformMTo (r,c)) ms
384 where
385 r = maxZ (map rows ms)
386 c = maxZ (map cols ms)
387
388
389conformVs vs = map (conformVTo n) vs
390 where
391 n = maxZ (map dim vs)
392
393conformMTo (r,c) m
394 | size m == (r,c) = m
395 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
396 | size m == (r,1) = repCols c m
397 | size m == (1,c) = repRows r m
398 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
399
400conformVTo n v
401 | dim v == n = v
402 | dim v == 1 = constantD (v@>0) n
403 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
404
405repRows n x = fromRows (replicate n (flatten x))
406repCols n x = fromColumns (replicate n (flatten x))
407
408size m = (rows m, cols m)
409
410shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
411
412emptyM r c = matrixFromVector RowMajor r c (fromList[])
413
414----------------------------------------------------------------------
415
416instance (Storable t, NFData t) => NFData (Matrix t)
417 where
418 rnf m | d > 0 = rnf (v @> 0)
419 | otherwise = ()
420 where
421 d = dim v
422 v = xdat m
423
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
deleted file mode 100644
index 257ad73..0000000
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ /dev/null
@@ -1,720 +0,0 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE TypeFamilies #-}
3{-# LANGUAGE FlexibleContexts #-}
4{-# LANGUAGE FlexibleInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE FunctionalDependencies #-}
7{-# LANGUAGE UndecidableInstances #-}
8
9-----------------------------------------------------------------------------
10-- |
11-- Module : Data.Packed.Internal.Numeric
12-- Copyright : (c) Alberto Ruiz 2010-14
13-- License : BSD3
14-- Maintainer : Alberto Ruiz
15-- Stability : provisional
16--
17-----------------------------------------------------------------------------
18
19module Data.Packed.Internal.Numeric (
20 -- * Basic functions
21 ident, diag, ctrans,
22 -- * Generic operations
23 Container(..),
24 scalar, conj, scale, arctan2, cmap,
25 atIndex, minIndex, maxIndex, minElement, maxElement,
26 sumElements, prodElements,
27 step, cond, find, assoc, accum,
28 Transposable(..), Linear(..), Testable(..),
29 -- * Matrix product and related functions
30 Product(..), udot,
31 mXm,mXv,vXm,
32 outer, kronecker,
33 -- * sorting
34 sortVector,
35 -- * Element conversion
36 Convert(..),
37 Complexable(),
38 RealElement(),
39 roundVector,
40 RealOf, ComplexOf, SingleOf, DoubleOf,
41 IndexOf,
42 module Data.Complex
43) where
44
45import Data.Packed
46import Data.Packed.ST as ST
47import Numeric.Conversion
48import Data.Packed.Development
49import Numeric.Vectorized
50import Data.Complex
51import Control.Applicative((<*>))
52
53import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
54import Data.Packed.Internal
55
56-------------------------------------------------------------------
57
58type family IndexOf (c :: * -> *)
59
60type instance IndexOf Vector = Int
61type instance IndexOf Matrix = (Int,Int)
62
63type family ArgOf (c :: * -> *) a
64
65type instance ArgOf Vector a = a -> a
66type instance ArgOf Matrix a = a -> a -> a
67
68-------------------------------------------------------------------
69
70-- | Basic element-by-element functions for numeric containers
71class (Complexable c, Fractional e, Element e) => Container c e
72 where
73 size' :: c e -> IndexOf c
74 scalar' :: e -> c e
75 conj' :: c e -> c e
76 scale' :: e -> c e -> c e
77 -- | scale the element by element reciprocal of the object:
78 --
79 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
80 scaleRecip :: e -> c e -> c e
81 addConstant :: e -> c e -> c e
82 add :: c e -> c e -> c e
83 sub :: c e -> c e -> c e
84 -- | element by element multiplication
85 mul :: c e -> c e -> c e
86 -- | element by element division
87 divide :: c e -> c e -> c e
88 equal :: c e -> c e -> Bool
89 --
90 -- element by element inverse tangent
91 arctan2' :: c e -> c e -> c e
92 cmap' :: (Element b) => (e -> b) -> c e -> c b
93 konst' :: e -> IndexOf c -> c e
94 build' :: IndexOf c -> (ArgOf c e) -> c e
95 atIndex' :: c e -> IndexOf c -> e
96 minIndex' :: c e -> IndexOf c
97 maxIndex' :: c e -> IndexOf c
98 minElement' :: c e -> e
99 maxElement' :: c e -> e
100 sumElements' :: c e -> e
101 prodElements' :: c e -> e
102 step' :: RealElement e => c e -> c e
103 cond' :: RealElement e
104 => c e -- ^ a
105 -> c e -- ^ b
106 -> c e -- ^ l
107 -> c e -- ^ e
108 -> c e -- ^ g
109 -> c e -- ^ result
110 find' :: (e -> Bool) -> c e -> [IndexOf c]
111 assoc' :: IndexOf c -- ^ size
112 -> e -- ^ default value
113 -> [(IndexOf c, e)] -- ^ association list
114 -> c e -- ^ result
115 accum' :: c e -- ^ initial structure
116 -> (e -> e -> e) -- ^ update function
117 -> [(IndexOf c, e)] -- ^ association list
118 -> c e -- ^ result
119
120--------------------------------------------------------------------------
121
122instance Container Vector Float
123 where
124 size' = dim
125 scale' = vectorMapValF Scale
126 scaleRecip = vectorMapValF Recip
127 addConstant = vectorMapValF AddConstant
128 add = vectorZipF Add
129 sub = vectorZipF Sub
130 mul = vectorZipF Mul
131 divide = vectorZipF Div
132 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
133 arctan2' = vectorZipF ATan2
134 scalar' x = fromList [x]
135 konst' = constantD
136 build' = buildV
137 conj' = id
138 cmap' = mapVector
139 atIndex' = (@>)
140 minIndex' = emptyErrorV "minIndex" (round . toScalarF MinIdx)
141 maxIndex' = emptyErrorV "maxIndex" (round . toScalarF MaxIdx)
142 minElement' = emptyErrorV "minElement" (toScalarF Min)
143 maxElement' = emptyErrorV "maxElement" (toScalarF Max)
144 sumElements' = sumF
145 prodElements' = prodF
146 step' = stepF
147 find' = findV
148 assoc' = assocV
149 accum' = accumV
150 cond' = condV condF
151
152instance Container Vector Double
153 where
154 size' = dim
155 scale' = vectorMapValR Scale
156 scaleRecip = vectorMapValR Recip
157 addConstant = vectorMapValR AddConstant
158 add = vectorZipR Add
159 sub = vectorZipR Sub
160 mul = vectorZipR Mul
161 divide = vectorZipR Div
162 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
163 arctan2' = vectorZipR ATan2
164 scalar' x = fromList [x]
165 konst' = constantD
166 build' = buildV
167 conj' = id
168 cmap' = mapVector
169 atIndex' = (@>)
170 minIndex' = emptyErrorV "minIndex" (round . toScalarR MinIdx)
171 maxIndex' = emptyErrorV "maxIndex" (round . toScalarR MaxIdx)
172 minElement' = emptyErrorV "minElement" (toScalarR Min)
173 maxElement' = emptyErrorV "maxElement" (toScalarR Max)
174 sumElements' = sumR
175 prodElements' = prodR
176 step' = stepD
177 find' = findV
178 assoc' = assocV
179 accum' = accumV
180 cond' = condV condD
181
182instance Container Vector (Complex Double)
183 where
184 size' = dim
185 scale' = vectorMapValC Scale
186 scaleRecip = vectorMapValC Recip
187 addConstant = vectorMapValC AddConstant
188 add = vectorZipC Add
189 sub = vectorZipC Sub
190 mul = vectorZipC Mul
191 divide = vectorZipC Div
192 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
193 arctan2' = vectorZipC ATan2
194 scalar' x = fromList [x]
195 konst' = constantD
196 build' = buildV
197 conj' = conjugateC
198 cmap' = mapVector
199 atIndex' = (@>)
200 minIndex' = emptyErrorV "minIndex" (minIndex' . fst . fromComplex . (mul <*> conj'))
201 maxIndex' = emptyErrorV "maxIndex" (maxIndex' . fst . fromComplex . (mul <*> conj'))
202 minElement' = emptyErrorV "minElement" (atIndex' <*> minIndex')
203 maxElement' = emptyErrorV "maxElement" (atIndex' <*> maxIndex')
204 sumElements' = sumC
205 prodElements' = prodC
206 step' = undefined -- cannot match
207 find' = findV
208 assoc' = assocV
209 accum' = accumV
210 cond' = undefined -- cannot match
211
212instance Container Vector (Complex Float)
213 where
214 size' = dim
215 scale' = vectorMapValQ Scale
216 scaleRecip = vectorMapValQ Recip
217 addConstant = vectorMapValQ AddConstant
218 add = vectorZipQ Add
219 sub = vectorZipQ Sub
220 mul = vectorZipQ Mul
221 divide = vectorZipQ Div
222 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
223 arctan2' = vectorZipQ ATan2
224 scalar' x = fromList [x]
225 konst' = constantD
226 build' = buildV
227 conj' = conjugateQ
228 cmap' = mapVector
229 atIndex' = (@>)
230 minIndex' = emptyErrorV "minIndex" (minIndex' . fst . fromComplex . (mul <*> conj'))
231 maxIndex' = emptyErrorV "maxIndex" (maxIndex' . fst . fromComplex . (mul <*> conj'))
232 minElement' = emptyErrorV "minElement" (atIndex' <*> minIndex')
233 maxElement' = emptyErrorV "maxElement" (atIndex' <*> maxIndex')
234 sumElements' = sumQ
235 prodElements' = prodQ
236 step' = undefined -- cannot match
237 find' = findV
238 assoc' = assocV
239 accum' = accumV
240 cond' = undefined -- cannot match
241
242---------------------------------------------------------------
243
244instance (Fractional a, Element a, Container Vector a) => Container Matrix a
245 where
246 size' = size
247 scale' x = liftMatrix (scale' x)
248 scaleRecip x = liftMatrix (scaleRecip x)
249 addConstant x = liftMatrix (addConstant x)
250 add = liftMatrix2 add
251 sub = liftMatrix2 sub
252 mul = liftMatrix2 mul
253 divide = liftMatrix2 divide
254 equal a b = cols a == cols b && flatten a `equal` flatten b
255 arctan2' = liftMatrix2 arctan2'
256 scalar' x = (1><1) [x]
257 konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c))
258 build' = buildM
259 conj' = liftMatrix conj'
260 cmap' f = liftMatrix (mapVector f)
261 atIndex' = (@@>)
262 minIndex' = emptyErrorM "minIndex of Matrix" $
263 \m -> divMod (minIndex' $ flatten m) (cols m)
264 maxIndex' = emptyErrorM "maxIndex of Matrix" $
265 \m -> divMod (maxIndex' $ flatten m) (cols m)
266 minElement' = emptyErrorM "minElement of Matrix" (atIndex' <*> minIndex')
267 maxElement' = emptyErrorM "maxElement of Matrix" (atIndex' <*> maxIndex')
268 sumElements' = sumElements . flatten
269 prodElements' = prodElements . flatten
270 step' = liftMatrix step
271 find' = findM
272 assoc' = assocM
273 accum' = accumM
274 cond' = condM
275
276
277emptyErrorV msg f v =
278 if dim v > 0
279 then f v
280 else error $ msg ++ " of Vector with dim = 0"
281
282emptyErrorM msg f m =
283 if rows m > 0 && cols m > 0
284 then f m
285 else error $ msg++" "++shSize m
286
287--------------------------------------------------------------------------------
288
289-- | create a structure with a single element
290--
291-- >>> let v = fromList [1..3::Double]
292-- >>> v / scalar (norm2 v)
293-- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732]
294--
295scalar :: Container c e => e -> c e
296scalar = scalar'
297
298-- | complex conjugate
299conj :: Container c e => c e -> c e
300conj = conj'
301
302-- | multiplication by scalar
303scale :: Container c e => e -> c e -> c e
304scale = scale'
305
306arctan2 :: Container c e => c e -> c e -> c e
307arctan2 = arctan2'
308
309-- | like 'fmap' (cannot implement instance Functor because of Element class constraint)
310cmap :: (Element b, Container c e) => (e -> b) -> c e -> c b
311cmap = cmap'
312
313-- | indexing function
314atIndex :: Container c e => c e -> IndexOf c -> e
315atIndex = atIndex'
316
317-- | index of minimum element
318minIndex :: Container c e => c e -> IndexOf c
319minIndex = minIndex'
320
321-- | index of maximum element
322maxIndex :: Container c e => c e -> IndexOf c
323maxIndex = maxIndex'
324
325-- | value of minimum element
326minElement :: Container c e => c e -> e
327minElement = minElement'
328
329-- | value of maximum element
330maxElement :: Container c e => c e -> e
331maxElement = maxElement'
332
333-- | the sum of elements
334sumElements :: Container c e => c e -> e
335sumElements = sumElements'
336
337-- | the product of elements
338prodElements :: Container c e => c e -> e
339prodElements = prodElements'
340
341
342-- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@
343--
344-- >>> step $ linspace 5 (-1,1::Double)
345-- 5 |> [0.0,0.0,0.0,1.0,1.0]
346--
347step
348 :: (RealElement e, Container c e)
349 => c e
350 -> c e
351step = step'
352
353
354-- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@.
355--
356-- Arguments with any dimension = 1 are automatically expanded:
357--
358-- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double
359-- (3><4)
360-- [ 100.0, 2.0, 3.0, 4.0
361-- , 0.0, 100.0, 7.0, 8.0
362-- , 0.0, 0.0, 100.0, 12.0 ]
363--
364cond
365 :: (RealElement e, Container c e)
366 => c e -- ^ a
367 -> c e -- ^ b
368 -> c e -- ^ l
369 -> c e -- ^ e
370 -> c e -- ^ g
371 -> c e -- ^ result
372cond = cond'
373
374
375-- | Find index of elements which satisfy a predicate
376--
377-- >>> find (>0) (ident 3 :: Matrix Double)
378-- [(0,0),(1,1),(2,2)]
379--
380find
381 :: Container c e
382 => (e -> Bool)
383 -> c e
384 -> [IndexOf c]
385find = find'
386
387
388-- | Create a structure from an association list
389--
390-- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double
391-- fromList [0.0,4.0,0.0,7.0,0.0]
392--
393-- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double)
394-- (2><3)
395-- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0
396-- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ]
397--
398assoc
399 :: Container c e
400 => IndexOf c -- ^ size
401 -> e -- ^ default value
402 -> [(IndexOf c, e)] -- ^ association list
403 -> c e -- ^ result
404assoc = assoc'
405
406
407-- | Modify a structure using an update function
408--
409-- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
410-- (5><5)
411-- [ 1.0, 0.0, 0.0, 3.0, 0.0
412-- , 0.0, 6.0, 0.0, 0.0, 0.0
413-- , 0.0, 0.0, 1.0, 0.0, 0.0
414-- , 0.0, 0.0, 0.0, 1.0, 0.0
415-- , 0.0, 0.0, 0.0, 0.0, 1.0 ]
416--
417-- computation of histogram:
418--
419-- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double
420-- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0]
421--
422accum
423 :: Container c e
424 => c e -- ^ initial structure
425 -> (e -> e -> e) -- ^ update function
426 -> [(IndexOf c, e)] -- ^ association list
427 -> c e -- ^ result
428accum = accum'
429
430
431--------------------------------------------------------------------------------
432
433-- | Matrix product and related functions
434class (Num e, Element e) => Product e where
435 -- | matrix product
436 multiply :: Matrix e -> Matrix e -> Matrix e
437 -- | sum of absolute value of elements (differs in complex case from @norm1@)
438 absSum :: Vector e -> RealOf e
439 -- | sum of absolute value of elements
440 norm1 :: Vector e -> RealOf e
441 -- | euclidean norm
442 norm2 :: Vector e -> RealOf e
443 -- | element of maximum magnitude
444 normInf :: Vector e -> RealOf e
445
446instance Product Float where
447 norm2 = emptyVal (toScalarF Norm2)
448 absSum = emptyVal (toScalarF AbsSum)
449 norm1 = emptyVal (toScalarF AbsSum)
450 normInf = emptyVal (maxElement . vectorMapF Abs)
451 multiply = emptyMul multiplyF
452
453instance Product Double where
454 norm2 = emptyVal (toScalarR Norm2)
455 absSum = emptyVal (toScalarR AbsSum)
456 norm1 = emptyVal (toScalarR AbsSum)
457 normInf = emptyVal (maxElement . vectorMapR Abs)
458 multiply = emptyMul multiplyR
459
460instance Product (Complex Float) where
461 norm2 = emptyVal (toScalarQ Norm2)
462 absSum = emptyVal (toScalarQ AbsSum)
463 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs)
464 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs)
465 multiply = emptyMul multiplyQ
466
467instance Product (Complex Double) where
468 norm2 = emptyVal (toScalarC Norm2)
469 absSum = emptyVal (toScalarC AbsSum)
470 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs)
471 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs)
472 multiply = emptyMul multiplyC
473
474emptyMul m a b
475 | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c)
476 | otherwise = m a b
477 where
478 r = rows a
479 x1 = cols a
480 x2 = rows b
481 c = cols b
482
483emptyVal f v =
484 if dim v > 0
485 then f v
486 else 0
487
488-- FIXME remove unused C wrappers
489-- | unconjugated dot product
490udot :: Product e => Vector e -> Vector e -> e
491udot u v
492 | dim u == dim v = val (asRow u `multiply` asColumn v)
493 | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product"
494 where
495 val m | dim u > 0 = m@@>(0,0)
496 | otherwise = 0
497
498----------------------------------------------------------
499
500-- synonym for matrix product
501mXm :: Product t => Matrix t -> Matrix t -> Matrix t
502mXm = multiply
503
504-- matrix - vector product
505mXv :: Product t => Matrix t -> Vector t -> Vector t
506mXv m v = flatten $ m `mXm` (asColumn v)
507
508-- vector - matrix product
509vXm :: Product t => Vector t -> Matrix t -> Vector t
510vXm v m = flatten $ (asRow v) `mXm` m
511
512{- | Outer product of two vectors.
513
514>>> fromList [1,2,3] `outer` fromList [5,2,3]
515(3><3)
516 [ 5.0, 2.0, 3.0
517 , 10.0, 4.0, 6.0
518 , 15.0, 6.0, 9.0 ]
519
520-}
521outer :: (Product t) => Vector t -> Vector t -> Matrix t
522outer u v = asColumn u `multiply` asRow v
523
524{- | Kronecker product of two matrices.
525
526@m1=(2><3)
527 [ 1.0, 2.0, 0.0
528 , 0.0, -1.0, 3.0 ]
529m2=(4><3)
530 [ 1.0, 2.0, 3.0
531 , 4.0, 5.0, 6.0
532 , 7.0, 8.0, 9.0
533 , 10.0, 11.0, 12.0 ]@
534
535>>> kronecker m1 m2
536(8><9)
537 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
538 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
539 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
540 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
541 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
542 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
543 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
544 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]
545
546-}
547kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
548kronecker a b = fromBlocks
549 . splitEvery (cols a)
550 . map (reshape (cols b))
551 . toRows
552 $ flatten a `outer` flatten b
553
554-------------------------------------------------------------------
555
556
557class Convert t where
558 real :: Container c t => c (RealOf t) -> c t
559 complex :: Container c t => c t -> c (ComplexOf t)
560 single :: Container c t => c t -> c (SingleOf t)
561 double :: Container c t => c t -> c (DoubleOf t)
562 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
563 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
564
565
566instance Convert Double where
567 real = id
568 complex = comp'
569 single = single'
570 double = id
571 toComplex = toComplex'
572 fromComplex = fromComplex'
573
574instance Convert Float where
575 real = id
576 complex = comp'
577 single = id
578 double = double'
579 toComplex = toComplex'
580 fromComplex = fromComplex'
581
582instance Convert (Complex Double) where
583 real = comp'
584 complex = id
585 single = single'
586 double = id
587 toComplex = toComplex'
588 fromComplex = fromComplex'
589
590instance Convert (Complex Float) where
591 real = comp'
592 complex = id
593 single = id
594 double = double'
595 toComplex = toComplex'
596 fromComplex = fromComplex'
597
598-------------------------------------------------------------------
599
600type family RealOf x
601
602type instance RealOf Double = Double
603type instance RealOf (Complex Double) = Double
604
605type instance RealOf Float = Float
606type instance RealOf (Complex Float) = Float
607
608type family ComplexOf x
609
610type instance ComplexOf Double = Complex Double
611type instance ComplexOf (Complex Double) = Complex Double
612
613type instance ComplexOf Float = Complex Float
614type instance ComplexOf (Complex Float) = Complex Float
615
616type family SingleOf x
617
618type instance SingleOf Double = Float
619type instance SingleOf Float = Float
620
621type instance SingleOf (Complex a) = Complex (SingleOf a)
622
623type family DoubleOf x
624
625type instance DoubleOf Double = Double
626type instance DoubleOf Float = Double
627
628type instance DoubleOf (Complex a) = Complex (DoubleOf a)
629
630type family ElementOf c
631
632type instance ElementOf (Vector a) = a
633type instance ElementOf (Matrix a) = a
634
635------------------------------------------------------------
636
637buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
638 where rs = map fromIntegral [0 .. (rc-1)]
639 cs = map fromIntegral [0 .. (cc-1)]
640
641buildV n f = fromList [f k | k <- ks]
642 where ks = map fromIntegral [0 .. (n-1)]
643
644--------------------------------------------------------
645-- | conjugate transpose
646ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
647ctrans = liftMatrix conj' . trans
648
649-- | Creates a square matrix with a given diagonal.
650diag :: (Num a, Element a) => Vector a -> Matrix a
651diag v = diagRect 0 v n n where n = dim v
652
653-- | creates the identity matrix of given dimension
654ident :: (Num a, Element a) => Int -> Matrix a
655ident n = diag (constantD 1 n)
656
657--------------------------------------------------------
658
659findV p x = foldVectorWithIndex g [] x where
660 g k z l = if p z then k:l else l
661
662findM p x = map ((`divMod` cols x)) $ findV p (flatten x)
663
664assocV n z xs = ST.runSTVector $ do
665 v <- ST.newVector z n
666 mapM_ (\(k,x) -> ST.writeVector v k x) xs
667 return v
668
669assocM (r,c) z xs = ST.runSTMatrix $ do
670 m <- ST.newMatrix z r c
671 mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs
672 return m
673
674accumV v0 f xs = ST.runSTVector $ do
675 v <- ST.thawVector v0
676 mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs
677 return v
678
679accumM m0 f xs = ST.runSTMatrix $ do
680 m <- ST.thawMatrix m0
681 mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs
682 return m
683
684----------------------------------------------------------------------
685
686condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t'
687 where
688 args@(a'':_) = conformMs [a,b,l,e,t]
689 [a', b', l', e', t'] = map flatten args
690
691condV f a b l e t = f a' b' l' e' t'
692 where
693 [a', b', l', e', t'] = conformVs [a,b,l,e,t]
694
695--------------------------------------------------------------------------------
696
697class Transposable m mt | m -> mt, mt -> m
698 where
699 -- | (conjugate) transpose
700 tr :: m -> mt
701
702instance (Container Vector t) => Transposable (Matrix t) (Matrix t)
703 where
704 tr = ctrans
705
706class Linear t v
707 where
708 scalarL :: t -> v
709 addL :: v -> v -> v
710 scaleL :: t -> v -> v
711
712
713class Testable t
714 where
715 checkT :: t -> (Bool, IO())
716 ioCheckT :: t -> IO (Bool, IO())
717 ioCheckT = return . checkT
718
719--------------------------------------------------------------------------------
720
diff --git a/packages/base/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs
deleted file mode 100644
index acc3070..0000000
--- a/packages/base/src/Data/Packed/Internal/Signatures.hs
+++ /dev/null
@@ -1,70 +0,0 @@
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
12module Data.Packed.Internal.Signatures where
13
14import Foreign.Ptr(Ptr)
15import Data.Complex(Complex)
16import Foreign.C.Types(CInt)
17
18type PF = Ptr Float --
19type PD = Ptr Double --
20type PQ = Ptr (Complex Float) --
21type PC = Ptr (Complex Double) --
22type TF = CInt -> PF -> IO CInt --
23type TFF = CInt -> PF -> TF --
24type TFV = CInt -> PF -> TV --
25type TVF = CInt -> PD -> TF --
26type TFFF = CInt -> PF -> TFF --
27type TV = CInt -> PD -> IO CInt --
28type TVV = CInt -> PD -> TV --
29type TVVV = CInt -> PD -> TVV --
30type TFM = CInt -> CInt -> PF -> IO CInt --
31type TFMFM = CInt -> CInt -> PF -> TFM --
32type TFMFMFM = CInt -> CInt -> PF -> TFMFM --
33type TM = CInt -> CInt -> PD -> IO CInt --
34type TMM = CInt -> CInt -> PD -> TM --
35type TVMM = CInt -> PD -> TMM --
36type TMVMM = CInt -> CInt -> PD -> TVMM --
37type TMMM = CInt -> CInt -> PD -> TMM --
38type TVM = CInt -> PD -> TM --
39type TVVM = CInt -> PD -> TVM --
40type TMV = CInt -> CInt -> PD -> TV --
41type TMMV = CInt -> CInt -> PD -> TMV --
42type TMVM = CInt -> CInt -> PD -> TVM --
43type TMMVM = CInt -> CInt -> PD -> TMVM --
44type TCM = CInt -> CInt -> PC -> IO CInt --
45type TCVCM = CInt -> PC -> TCM --
46type TCMCVCM = CInt -> CInt -> PC -> TCVCM --
47type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM --
48type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM --
49type TCMCM = CInt -> CInt -> PC -> TCM --
50type TVCM = CInt -> PD -> TCM --
51type TCMVCM = CInt -> CInt -> PC -> TVCM --
52type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM --
53type TCMCMCM = CInt -> CInt -> PC -> TCMCM --
54type TCV = CInt -> PC -> IO CInt --
55type TCVCV = CInt -> PC -> TCV --
56type TCVCVCV = CInt -> PC -> TCVCV --
57type TCVV = CInt -> PC -> TV --
58type TQV = CInt -> PQ -> IO CInt --
59type TQVQV = CInt -> PQ -> TQV --
60type TQVQVQV = CInt -> PQ -> TQVQV --
61type TQVF = CInt -> PQ -> TF --
62type TQM = CInt -> CInt -> PQ -> IO CInt --
63type TQMQM = CInt -> CInt -> PQ -> TQM --
64type TQMQMQM = CInt -> CInt -> PQ -> TQMQM --
65type TCMCV = CInt -> CInt -> PC -> TCV --
66type TVCV = CInt -> PD -> TCV --
67type TCVM = CInt -> PC -> TM --
68type TMCVM = CInt -> CInt -> PD -> TCVM --
69type 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
deleted file mode 100644
index d0bc143..0000000
--- a/packages/base/src/Data/Packed/Internal/Vector.hs
+++ /dev/null
@@ -1,471 +0,0 @@
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
13module 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
30import Data.Packed.Internal.Common
31import Data.Packed.Internal.Signatures
32import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
33import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
34import Foreign.Ptr(Ptr)
35import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
36import Foreign.C.Types
37import Data.Complex
38import Control.Monad(when)
39import System.IO.Unsafe(unsafePerformIO)
40
41#if __GLASGOW_HASKELL__ >= 605
42import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
43#else
44import Foreign.ForeignPtr (mallocForeignPtrBytes)
45#endif
46
47import GHC.Base
48#if __GLASGOW_HASKELL__ < 612
49import GHC.IOBase hiding (liftIO)
50#endif
51
52import qualified Data.Vector.Storable as Vector
53import Data.Vector.Storable(Vector,
54 fromList,
55 unsafeToForeignPtr,
56 unsafeFromForeignPtr,
57 unsafeWith)
58
59
60-- | Number of elements
61dim :: (Storable t) => Vector t -> Int
62dim = Vector.length
63
64
65-- C-Haskell vector adapter
66-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
67vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
68vec 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
76createVector :: Storable a => Int -> IO (Vector a)
77createVector 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]
974 |> [2.0,3.0,5.0,7.0]@
98
99-}
100
101safeRead v = inlinePerformIO . unsafeWith v
102{-# INLINE safeRead #-}
103
104inlinePerformIO :: IO a -> a
105inlinePerformIO (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-}
114toList :: Storable a => Vector a -> [a]
115toList 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..]
122fromList [1.0,2.0,3.0,4.0,5.0]
123
124-}
125(|>) :: (Storable a) => Int -> [a] -> Vector a
126infixl 9 |>
127n |> 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
134at' :: Storable a => Vector a -> Int -> a
135at' 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)
143safe :: Bool
144safe = False
145#else
146safe = True
147#endif
148
149-- | access to Vector elements with range checking.
150at :: Storable a => Vector a -> Int -> a
151at 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])
161fromList [3.0,4.0,5.0]
162
163-}
164subVector :: Storable t => Int -- ^ index of the starting element
165 -> Int -- ^ number of elements to extract
166 -> Vector t -- ^ source
167 -> Vector t -- ^ result
168subVector = Vector.slice
169
170
171{- | Reads a vector position:
172
173>>> fromList [0..9] @> 7
1747.0
175
176-}
177(@>) :: Storable t => Vector t -> Int -> t
178infixl 9 @>
179(@>) = at
180
181
182{- | concatenate a list of vectors
183
184>>> vjoin [fromList [1..5::Double], konst 1 3]
185fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
186
187-}
188vjoin :: Storable t => [Vector t] -> Vector t
189vjoin [] = fromList []
190vjoin [v] = v
191vjoin 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-}
210takesV :: Storable t => [Int] -> Vector t -> [Vector t]
211takesV 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
220asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
221asReal 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
225asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
226asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
227 where (fp,i,n) = unsafeToForeignPtr v
228
229---------------------------------------------------------------
230
231float2DoubleV :: Vector Float -> Vector Double
232float2DoubleV v = unsafePerformIO $ do
233 r <- createVector (dim v)
234 app2 c_float2double vec v vec r "float2double"
235 return r
236
237double2FloatV :: Vector Double -> Vector Float
238double2FloatV v = unsafePerformIO $ do
239 r <- createVector (dim v)
240 app2 c_double2float vec v vec r "double2float2"
241 return r
242
243
244foreign import ccall unsafe "float2double" c_float2double:: TFV
245foreign import ccall unsafe "double2float" c_double2float:: TVF
246
247---------------------------------------------------------------
248
249stepF :: Vector Float -> Vector Float
250stepF v = unsafePerformIO $ do
251 r <- createVector (dim v)
252 app2 c_stepF vec v vec r "stepF"
253 return r
254
255stepD :: Vector Double -> Vector Double
256stepD v = unsafePerformIO $ do
257 r <- createVector (dim v)
258 app2 c_stepD vec v vec r "stepD"
259 return r
260
261foreign import ccall unsafe "stepF" c_stepF :: TFF
262foreign import ccall unsafe "stepD" c_stepD :: TVV
263
264---------------------------------------------------------------
265
266condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
267condF 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
272condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
273condD 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
278foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF
279foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV
280
281--------------------------------------------------------------------------------
282
283conjugateAux fun x = unsafePerformIO $ do
284 v <- createVector (dim x)
285 app2 fun vec x vec v "conjugateAux"
286 return v
287
288conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
289conjugateQ = conjugateAux c_conjugateQ
290foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV
291
292conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
293conjugateC = conjugateAux c_conjugateC
294foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV
295
296--------------------------------------------------------------------------------
297
298cloneVector :: Storable t => Vector t -> IO (Vector t)
299cloneVector 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
309mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
310mapVector 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
323zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
324zipVectorWith 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
340unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
341 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
342unzipVectorWith 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
359foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
360foldVector 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
369foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
370foldVectorWithIndex 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
378foldLoop 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
383foldVectorG 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
392mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
393mapVectorM 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
410mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
411mapVectorM_ 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
425mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
426mapVectorWithIndexM 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
443mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
444mapVectorWithIndexM_ 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
457mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
458--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
459mapVectorWithIndex 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-10 11:42:29 +0000