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
|
{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
-----------------------------------------------------------------------------
{- |
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(multiply,partit)
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 Vector a, Container Matrix a) => (Linear Matrix a) 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 && flatten a `equal` flatten b
--------------------------------------------------
-- | euclidean inner product
dot :: (Element t) => Vector t -> Vector t -> t
dot u v = multiply r c @@> (0,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
|
