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