1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.Chain
-- Copyright : (c) Alberto Ruiz 2010
-- License : GPL-style
--
-- Maintainer : Alberto Ruiz <aruiz@um.es>
-- Stability : provisional
-- Portability : portable
--
-- optimisation of association order for chains of matrix multiplication
--
-----------------------------------------------------------------------------
module Numeric.Chain (
chain
) where
import Data.Maybe
import Data.Packed.Matrix
import Numeric.Container
import qualified Data.Array.IArray as A
type Matrices a = A.Array Int (Matrix a)
type Sizes = A.Array Int (Int,Int)
type Cost = A.Array Int (A.Array Int (Maybe Int))
type Indexes = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))
update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a)
update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])]
newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int))
newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]
where subArray i = A.listArray (1,i) (repeat Nothing)
newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))
newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]
where subArray i = A.listArray (1,i) (repeat Nothing)
matricesToSizes :: [Matrix a] -> Sizes
matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms
-- | provide optimal association order for a chain of matrix multiplications and apply the multiplications
chain :: Product a => [Matrix a] -> Matrix a
chain ms = let ln = length ms
ma = A.listArray (1,ln) ms
mz = matricesToSizes ms
i = chain_cost mz
in chain_paren (ln,ln) i ma
chain_cost :: Sizes -> Indexes
chain_cost mz = let (_,u) = A.bounds mz
cost = newWorkSpaceCost u
ixes = newWorkSpaceIndexes u
(_,_,i) = foldl chain_cost' (mz,cost,ixes) (order u)
in i
chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)
chain_cost' sci@(mz,cost,ixes) (r,c)
| c == 1 = let cost' = update cost (r,c) (Just 0)
ixes' = update ixes (r,c) (Just ((r,c),(r,c)))
in (mz,cost',ixes')
| otherwise = minimum_cost sci (r,c)
minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)
minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu)
smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes)
smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost = (fromJust ((cost A.! lr) A.! lc))
+ (fromJust ((cost A.! rr) A.! rc))
+ ((fst $ mz A.! (lr-lc+1))
*(snd $ mz A.! lc)
*(snd $ mz A.! rr))
cost' = (cost A.! r) A.! c
in case cost' of
Nothing -> let cost'' = update cost (r,c) (Just op_cost)
ixes'' = update ixes (r,c) (Just ix)
in (mz,cost'',ixes'')
Just ct -> if op_cost < ct then
let cost'' = update cost (r,c) (Just op_cost)
ixes'' = update ixes (r,c) (Just ix)
in (mz,cost'',ixes'')
else (mz,cost,ixes)
fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)]
in map (partner (r,c)) fs'
partner (r,c) (a,b) = (((r-b),(c-b)),(a,b))
order 0 = []
order n = (order (n-1)) ++ (zip (repeat n) [1..n])
chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a
chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c
in if lr == rr && lc == rc then (ma A.! lr)
else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma)
--------------------------------------------------------------------------
{- TESTS -}
-- optimal association is ((m1*(m2*m3))*m4)
m1, m2, m3, m4 :: Matrix Double
m1 = (10><15) [1..]
m2 = (15><20) [1..]
m3 = (20><5) [1..]
m4 = (5><10) [1..]
|
