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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
{-# OPTIONS_GHC -fglasgow-exts #-}
-----------------------------------------------------------------------------
{- |
Module : Numeric.LinearAlgebra.Linear
Copyright : (c) Alberto Ruiz 2006-7
License : GPL-style
Maintainer : Alberto Ruiz (aruiz at um dot es)
Stability : provisional
Portability : uses ffi
Basic optimized operations on vectors and matrices.
-}
-----------------------------------------------------------------------------
module Numeric.LinearAlgebra.Linear (
Linear(..),
multiply, dot, outer, kronecker
) where
import Data.Packed.Internal
import Data.Packed
import Numeric.GSL.Vector
import Complex
-- | A generic interface for vectors and matrices to a few element-by-element functions in Numeric.GSL.Vector.
class (Container c e) => Linear c e where
scale :: e -> c e -> c e
addConstant :: e -> c e -> c e
add :: c e -> c e -> c e
sub :: c e -> c e -> c e
-- | element by element multiplication
mul :: c e -> c e -> c e
-- | element by element division
divide :: c e -> c e -> c e
-- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
scaleRecip :: e -> c e -> c e
equal :: c e -> c e -> Bool
-- numequal :: Double -> c e -> c e -> Bool
instance Linear Vector Double where
scale = vectorMapValR Scale
scaleRecip = vectorMapValR Recip
addConstant = vectorMapValR AddConstant
add = vectorZipR Add
sub = vectorZipR Sub
mul = vectorZipR Mul
divide = vectorZipR Div
equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0
instance Linear Vector (Complex Double) where
scale = vectorMapValC Scale
scaleRecip = vectorMapValC Recip
addConstant = vectorMapValC AddConstant
add = vectorZipC Add
sub = vectorZipC Sub
mul = vectorZipC Mul
divide = vectorZipC Div
equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0
instance Linear Matrix Double where
scale x = liftMatrix (scale x)
scaleRecip x = liftMatrix (scaleRecip x)
addConstant x = liftMatrix (addConstant x)
add = liftMatrix2 add
sub = liftMatrix2 sub
mul = liftMatrix2 mul
divide = liftMatrix2 divide
equal a b = cols a == cols b && cdat a `equal` cdat b
instance Linear Matrix (Complex Double) where
scale x = liftMatrix (scale x)
scaleRecip x = liftMatrix (scaleRecip x)
addConstant x = liftMatrix (addConstant x)
add = liftMatrix2 add
sub = liftMatrix2 sub
mul = liftMatrix2 mul
divide = liftMatrix2 divide
equal a b = cols a == cols b && cdat a `equal` cdat b
--------------------------------------------------
-- | euclidean inner product
dot :: (Element t) => Vector t -> Vector t -> t
dot u v = dat (multiply r c) `at` 0
where r = asRow u
c = asColumn v
{- | Outer product of two vectors.
@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
(3><3)
[ 5.0, 2.0, 3.0
, 10.0, 4.0, 6.0
, 15.0, 6.0, 9.0 ]@
-}
outer :: (Element t) => Vector t -> Vector t -> Matrix t
outer u v = asColumn u `multiply` asRow v
{- | Kronecker product of two matrices.
@m1=(2><3)
[ 1.0, 2.0, 0.0
, 0.0, -1.0, 3.0 ]
m2=(4><3)
[ 1.0, 2.0, 3.0
, 4.0, 5.0, 6.0
, 7.0, 8.0, 9.0
, 10.0, 11.0, 12.0 ]@
@\> kronecker m1 m2
(8><9)
[ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
, 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
, 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
, 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
, 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
, 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
-}
kronecker :: (Element t) => Matrix t -> Matrix t -> Matrix t
kronecker a b = fromBlocks
. partit (cols a)
. map (reshape (cols b))
. toRows
$ flatten a `outer` flatten b
|
