summaryrefslogtreecommitdiff
path: root/packages/base/src/Data/Packed/Matrix.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Data/Packed/Matrix.hs')
-rw-r--r--packages/base/src/Data/Packed/Matrix.hs494
1 files changed, 0 insertions, 494 deletions
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
deleted file mode 100644
index 70b9232..0000000
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ /dev/null
@@ -1,494 +0,0 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE CPP #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Data.Packed.Matrix
11-- Copyright : (c) Alberto Ruiz 2007-10
12-- License : BSD3
13-- Maintainer : Alberto Ruiz
14-- Stability : provisional
15--
16-- A Matrix representation suitable for numerical computations using LAPACK and GSL.
17--
18-- This module provides basic functions for manipulation of structure.
19
20-----------------------------------------------------------------------------
21{-# OPTIONS_HADDOCK hide #-}
22
23module Data.Packed.Matrix (
24 Matrix,
25 Element,
26 rows,cols,
27 (><),
28 trans,
29 reshape, flatten,
30 fromLists, toLists, buildMatrix,
31 (@@>),
32 asRow, asColumn,
33 fromRows, toRows, fromColumns, toColumns,
34 fromBlocks, diagBlock, toBlocks, toBlocksEvery,
35 repmat,
36 flipud, fliprl,
37 subMatrix, takeRows, dropRows, takeColumns, dropColumns,
38 extractRows, extractColumns,
39 diagRect, takeDiag,
40 mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
41 liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
42) where
43
44import Data.Packed.Internal
45import qualified Data.Packed.ST as ST
46import Data.Array
47
48import Data.List(transpose,intersperse)
49import Foreign.Storable(Storable)
50import Control.Monad(liftM)
51
52-------------------------------------------------------------------
53
54#ifdef BINARY
55
56import Data.Binary
57
58instance (Binary (Vector a), Element a) => Binary (Matrix a) where
59 put m = do
60 put (cols m)
61 put (flatten m)
62 get = do
63 c <- get
64 v <- get
65 return (reshape c v)
66
67#endif
68
69-------------------------------------------------------------------
70
71instance (Show a, Element a) => (Show (Matrix a)) where
72 show m | rows m == 0 || cols m == 0 = sizes m ++" []"
73 show m = (sizes m++) . dsp . map (map show) . toLists $ m
74
75sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
76
77dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
78 where
79 mt = transpose as
80 longs = map (maximum . map length) mt
81 mtp = zipWith (\a b -> map (pad a) b) longs mt
82 pad n str = replicate (n - length str) ' ' ++ str
83 unwords' = concat . intersperse ", "
84
85------------------------------------------------------------------
86
87instance (Element a, Read a) => Read (Matrix a) where
88 readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
89 where (thing,rest) = breakAt ']' s
90 (dims,listnums) = breakAt ')' thing
91 cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
92 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
93
94
95breakAt c l = (a++[c],tail b) where
96 (a,b) = break (==c) l
97
98------------------------------------------------------------------
99
100-- | creates a matrix from a vertical list of matrices
101joinVert :: Element t => [Matrix t] -> Matrix t
102joinVert [] = emptyM 0 0
103joinVert ms = case common cols ms of
104 Nothing -> error "(impossible) joinVert on matrices with different number of columns"
105 Just c -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms)
106
107-- | creates a matrix from a horizontal list of matrices
108joinHoriz :: Element t => [Matrix t] -> Matrix t
109joinHoriz ms = trans. joinVert . map trans $ ms
110
111{- | Create a matrix from blocks given as a list of lists of matrices.
112
113Single row-column components are automatically expanded to match the
114corresponding common row and column:
115
116@
117disp = putStr . dispf 2
118@
119
120>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
1218x10
1221 0 0 0 0 7 7 7 10 20
1230 1 0 0 0 7 7 7 10 20
1240 0 1 0 0 7 7 7 10 20
1250 0 0 1 0 7 7 7 10 20
1260 0 0 0 1 7 7 7 10 20
1273 3 3 3 3 1 0 0 0 0
1283 3 3 3 3 0 2 0 0 0
1293 3 3 3 3 0 0 3 0 0
130
131-}
132fromBlocks :: Element t => [[Matrix t]] -> Matrix t
133fromBlocks = fromBlocksRaw . adaptBlocks
134
135fromBlocksRaw mms = joinVert . map joinHoriz $ mms
136
137adaptBlocks ms = ms' where
138 bc = case common length ms of
139 Just c -> c
140 Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
141 rs = map (compatdim . map rows) ms
142 cs = map (compatdim . map cols) (transpose ms)
143 szs = sequence [rs,cs]
144 ms' = splitEvery bc $ zipWith g szs (concat ms)
145
146 g [Just nr,Just nc] m
147 | nr == r && nc == c = m
148 | r == 1 && c == 1 = matrixFromVector RowMajor nr nc (constantD x (nr*nc))
149 | r == 1 = fromRows (replicate nr (flatten m))
150 | otherwise = fromColumns (replicate nc (flatten m))
151 where
152 r = rows m
153 c = cols m
154 x = m@@>(0,0)
155 g _ _ = error "inconsistent dimensions in fromBlocks"
156
157
158--------------------------------------------------------------------------------
159
160{- | create a block diagonal matrix
161
162>>> disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
1637x8
1641 1 0 0 0 0 0 0
1651 1 0 0 0 0 0 0
1660 0 2 2 2 2 2 0
1670 0 2 2 2 2 2 0
1680 0 2 2 2 2 2 0
1690 0 0 0 0 0 0 5
1700 0 0 0 0 0 0 7
171
172>>> diagBlock [(0><4)[], konst 2 (2,3)] :: Matrix Double
173(2><7)
174 [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0
175 , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ]
176
177-}
178diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
179diagBlock ms = fromBlocks $ zipWith f ms [0..]
180 where
181 f m k = take n $ replicate k z ++ m : repeat z
182 n = length ms
183 z = (1><1) [0]
184
185--------------------------------------------------------------------------------
186
187
188-- | Reverse rows
189flipud :: Element t => Matrix t -> Matrix t
190flipud m = extractRows [r-1,r-2 .. 0] $ m
191 where
192 r = rows m
193
194-- | Reverse columns
195fliprl :: Element t => Matrix t -> Matrix t
196fliprl m = extractColumns [c-1,c-2 .. 0] $ m
197 where
198 c = cols m
199
200------------------------------------------------------------
201
202{- | creates a rectangular diagonal matrix:
203
204>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
205(4><5)
206 [ 10.0, 7.0, 7.0, 7.0, 7.0
207 , 7.0, 20.0, 7.0, 7.0, 7.0
208 , 7.0, 7.0, 30.0, 7.0, 7.0
209 , 7.0, 7.0, 7.0, 7.0, 7.0 ]
210
211-}
212diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
213diagRect z v r c = ST.runSTMatrix $ do
214 m <- ST.newMatrix z r c
215 let d = min r c `min` (dim v)
216 mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
217 return m
218
219-- | extracts the diagonal from a rectangular matrix
220takeDiag :: (Element t) => Matrix t -> Vector t
221takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
222
223------------------------------------------------------------
224
225{- | create a general matrix
226
227>>> (2><3) [2, 4, 7+2*𝑖, -3, 11, 0]
228(2><3)
229 [ 2.0 :+ 0.0, 4.0 :+ 0.0, 7.0 :+ 2.0
230 , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ]
231
232The input list is explicitly truncated, so that it can
233safely be used with lists that are too long (like infinite lists).
234
235>>> (2><3)[1..]
236(2><3)
237 [ 1.0, 2.0, 3.0
238 , 4.0, 5.0, 6.0 ]
239
240This is the format produced by the instances of Show (Matrix a), which
241can also be used for input.
242
243-}
244(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
245r >< c = f where
246 f l | dim v == r*c = matrixFromVector RowMajor r c v
247 | otherwise = error $ "inconsistent list size = "
248 ++show (dim v) ++" in ("++show r++"><"++show c++")"
249 where v = fromList $ take (r*c) l
250
251----------------------------------------------------------------
252
253-- | Creates a matrix with the first n rows of another matrix
254takeRows :: Element t => Int -> Matrix t -> Matrix t
255takeRows n mt = subMatrix (0,0) (n, cols mt) mt
256-- | Creates a copy of a matrix without the first n rows
257dropRows :: Element t => Int -> Matrix t -> Matrix t
258dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
259-- |Creates a matrix with the first n columns of another matrix
260takeColumns :: Element t => Int -> Matrix t -> Matrix t
261takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
262-- | Creates a copy of a matrix without the first n columns
263dropColumns :: Element t => Int -> Matrix t -> Matrix t
264dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
265
266----------------------------------------------------------------
267
268{- | Creates a 'Matrix' from a list of lists (considered as rows).
269
270>>> fromLists [[1,2],[3,4],[5,6]]
271(3><2)
272 [ 1.0, 2.0
273 , 3.0, 4.0
274 , 5.0, 6.0 ]
275
276-}
277fromLists :: Element t => [[t]] -> Matrix t
278fromLists = fromRows . map fromList
279
280-- | creates a 1-row matrix from a vector
281--
282-- >>> asRow (fromList [1..5])
283-- (1><5)
284-- [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
285--
286asRow :: Storable a => Vector a -> Matrix a
287asRow = trans . asColumn
288
289-- | creates a 1-column matrix from a vector
290--
291-- >>> asColumn (fromList [1..5])
292-- (5><1)
293-- [ 1.0
294-- , 2.0
295-- , 3.0
296-- , 4.0
297-- , 5.0 ]
298--
299asColumn :: Storable a => Vector a -> Matrix a
300asColumn v = reshape 1 v
301
302
303
304{- | creates a Matrix of the specified size using the supplied function to
305 to map the row\/column position to the value at that row\/column position.
306
307@> buildMatrix 3 4 (\\(r,c) -> fromIntegral r * fromIntegral c)
308(3><4)
309 [ 0.0, 0.0, 0.0, 0.0, 0.0
310 , 0.0, 1.0, 2.0, 3.0, 4.0
311 , 0.0, 2.0, 4.0, 6.0, 8.0]@
312
313Hilbert matrix of order N:
314
315@hilb n = buildMatrix n n (\\(i,j)->1/(fromIntegral i + fromIntegral j +1))@
316
317-}
318buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a
319buildMatrix rc cc f =
320 fromLists $ map (map f)
321 $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)]
322
323-----------------------------------------------------
324
325fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e
326fromArray2D m = (r><c) (elems m)
327 where ((r0,c0),(r1,c1)) = bounds m
328 r = r1-r0+1
329 c = c1-c0+1
330
331
332-- | rearranges the rows of a matrix according to the order given in a list of integers.
333extractRows :: Element t => [Int] -> Matrix t -> Matrix t
334extractRows [] m = emptyM 0 (cols m)
335extractRows l m = fromRows $ extract (toRows m) l
336 where
337 extract l' is = [l'!!i | i<- map verify is]
338 verify k
339 | k >= 0 && k < rows m = k
340 | otherwise = error $ "can't extract row "
341 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
342
343-- | rearranges the rows of a matrix according to the order given in a list of integers.
344extractColumns :: Element t => [Int] -> Matrix t -> Matrix t
345extractColumns l m = trans . extractRows (map verify l) . trans $ m
346 where
347 verify k
348 | k >= 0 && k < cols m = k
349 | otherwise = error $ "can't extract column "
350 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
351
352
353
354{- | creates matrix by repetition of a matrix a given number of rows and columns
355
356>>> repmat (ident 2) 2 3
357(4><6)
358 [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
359 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
360 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
361 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
362
363-}
364repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
365repmat m r c
366 | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m)
367 | otherwise = fromBlocks $ replicate r $ replicate c $ m
368
369-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
370liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
371liftMatrix2Auto f m1 m2
372 | compat' m1 m2 = lM f m1 m2
373 | ok = lM f m1' m2'
374 | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2
375 where
376 (r1,c1) = size m1
377 (r2,c2) = size m2
378 r = max r1 r2
379 c = max c1 c2
380 r0 = min r1 r2
381 c0 = min c1 c2
382 ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2
383 m1' = conformMTo (r,c) m1
384 m2' = conformMTo (r,c) m2
385
386-- FIXME do not flatten if equal order
387lM f m1 m2 = matrixFromVector
388 RowMajor
389 (max (rows m1) (rows m2))
390 (max (cols m1) (cols m2))
391 (f (flatten m1) (flatten m2))
392
393compat' :: Matrix a -> Matrix b -> Bool
394compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
395 where
396 s1 = size m1
397 s2 = size m2
398
399------------------------------------------------------------
400
401toBlockRows [r] m
402 | r == rows m = [m]
403toBlockRows rs m
404 | cols m > 0 = map (reshape (cols m)) (takesV szs (flatten m))
405 | otherwise = map g rs
406 where
407 szs = map (* cols m) rs
408 g k = (k><0)[]
409
410toBlockCols [c] m | c == cols m = [m]
411toBlockCols cs m = map trans . toBlockRows cs . trans $ m
412
413-- | Partition a matrix into blocks with the given numbers of rows and columns.
414-- The remaining rows and columns are discarded.
415toBlocks :: (Element t) => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
416toBlocks rs cs m
417 | ok = map (toBlockCols cs) . toBlockRows rs $ m
418 | otherwise = error $ "toBlocks: bad partition: "++show rs++" "++show cs
419 ++ " "++shSize m
420 where
421 ok = sum rs <= rows m && sum cs <= cols m && all (>=0) rs && all (>=0) cs
422
423-- | Fully partition a matrix into blocks of the same size. If the dimensions are not
424-- a multiple of the given size the last blocks will be smaller.
425toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]]
426toBlocksEvery r c m
427 | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c
428 | otherwise = toBlocks rs cs m
429 where
430 (qr,rr) = rows m `divMod` r
431 (qc,rc) = cols m `divMod` c
432 rs = replicate qr r ++ if rr > 0 then [rr] else []
433 cs = replicate qc c ++ if rc > 0 then [rc] else []
434
435-------------------------------------------------------------------
436
437-- Given a column number and a function taking matrix indexes, returns
438-- a function which takes vector indexes (that can be used on the
439-- flattened matrix).
440mk :: Int -> ((Int, Int) -> t) -> (Int -> t)
441mk c g = \k -> g (divMod k c)
442
443{- |
444
445>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
446m[0,0] = 1
447m[0,1] = 2
448m[0,2] = 3
449m[1,0] = 4
450m[1,1] = 5
451m[1,2] = 6
452
453-}
454mapMatrixWithIndexM_
455 :: (Element a, Num a, Monad m) =>
456 ((Int, Int) -> a -> m ()) -> Matrix a -> m ()
457mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
458 where
459 c = cols m
460
461{- |
462
463>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
464Just (3><3)
465 [ 100.0, 1.0, 2.0
466 , 10.0, 111.0, 12.0
467 , 20.0, 21.0, 122.0 ]
468
469-}
470mapMatrixWithIndexM
471 :: (Element a, Storable b, Monad m) =>
472 ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
473mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
474 where
475 c = cols m
476
477{- |
478
479>>> mapMatrixWithIndex (\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
480(3><3)
481 [ 100.0, 1.0, 2.0
482 , 10.0, 111.0, 12.0
483 , 20.0, 21.0, 122.0 ]
484
485 -}
486mapMatrixWithIndex
487 :: (Element a, Storable b) =>
488 ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
489mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
490 where
491 c = cols m
492
493mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
494mapMatrix f = liftMatrix (mapVector f)