summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/oldtests.hs2
-rw-r--r--examples/pru.hs159
-rw-r--r--examples/tests.hs2
-rw-r--r--lib/Data/Packed/Internal/Tensor.hs32
-rw-r--r--lib/Data/Packed/Matrix.hs4
-rw-r--r--lib/Data/Packed/Tensor.hs82
-rw-r--r--lib/GSL/Vector.hs4
7 files changed, 146 insertions, 139 deletions
diff --git a/examples/oldtests.hs b/examples/oldtests.hs
index 987ef98..1f1ca25 100644
--- a/examples/oldtests.hs
+++ b/examples/oldtests.hs
@@ -10,7 +10,7 @@ realVector = fromList :: [Double] -> Vector Double
10toComplexM = uncurry $ liftMatrix2 (curry toComplex) 10toComplexM = uncurry $ liftMatrix2 (curry toComplex)
11 11
12infixl 2 =~= 12infixl 2 =~=
13a =~= b = pnorm PNorm1 (flatten (a - b)) < 1E-6 13a =~= b = pnorm 1 (flatten (a - b)) < 1E-6
14 14
15randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed 15randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
16 16
diff --git a/examples/pru.hs b/examples/pru.hs
index f855bce..71f33bf 100644
--- a/examples/pru.hs
+++ b/examples/pru.hs
@@ -1,56 +1,33 @@
1--{-# OPTIONS_GHC #-} 1--{-# OPTIONS_GHC #-}
2--module Main where 2--module Main where
3 3{-
4import Data.Packed.Internal 4import Data.Packed.Internal
5import Data.Packed.Internal.Vector 5import Data.Packed.Internal.Vector
6import Data.Packed.Internal.Matrix 6import Data.Packed.Internal.Matrix
7import Data.Packed.Internal.Tensor 7import Data.Packed.Tensor
8import Data.Packed.Matrix 8import Data.Packed.Matrix
9import GSL.Vector 9import GSL.Vector
10import LAPACK 10import LAPACK
11import Data.List(foldl')
11 12
12import Complex 13import Complex
13import Numeric(showGFloat) 14import Numeric(showGFloat)
14import Data.List(transpose,intersperse,sort,elemIndex,nub,foldl',foldl1')
15import Foreign.Storable
16
17
18vr = fromList [1..15::Double]
19vc = fromList (map (\x->x :+ (x+1)) [1..15::Double])
20
21mi = (2 >< 3) [1 .. 6::Int]
22mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)]
23
24ac = (2><3) [1 .. 6::Double]
25bc = (3><4) [7 .. 18::Double]
26
27af = (2>|<3) [1,4,2,5,3,6::Double]
28bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
29
30
31a |=| b = rows a == rows b &&
32 cols a == cols b &&
33 toList (cdat a) == toList (cdat b)
34
35mulC a b = multiply RowMajor a b
36mulF a b = multiply ColumnMajor a b
37
38cc = mulC ac bf
39cf = mulF af bc
40 15
41r = mulC cc (trans cf) 16import Foreign.Storable
42 17-}
43rd = (2><2)
44 [ 27736.0, 65356.0
45 , 65356.0, 154006.0 ]
46 18
19import GSL
20import Data.List(foldl', foldl1')
21import Data.Packed.Internal.Tensor
22import Data.Packed.Tensor
23import Complex
47 24
48 25
49main = do 26main = do
50 print $ r |=| rd
51 print $ foldl part t [("p",1),("q",0),("r",2)] 27 print $ foldl part t [("p",1),("q",0),("r",2)]
52 print $ foldl part t [("p",1),("r",2),("q",0)] 28 print $ foldl part t [("p",1),("r",2),("q",0)]
53 print $ foldl part t $ reverse [("p",1),("r",2),("q",0)] 29 print $ foldl part t $ reverse [("p",1),("r",2),("q",0)]
30 pru
54 31
55t = T [IdxDesc 4 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 3 Contravariant "r"] 32t = T [IdxDesc 4 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 3 Contravariant "r"]
56 $ fromList [1..24::Double] 33 $ fromList [1..24::Double]
@@ -61,10 +38,17 @@ t1 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q" ,IdxDesc 2 Covariant
61t2 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q"] 38t2 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q"]
62 $ fromList [1..16::Double] 39 $ fromList [1..16::Double]
63 40
41vector n v = tvector n (fromList v) :: Tensor Double
42
43kron n = tensorFromMatrix (Contravariant,"k1") (Covariant,"k2") (ident n)
64 44
45tensorFromTrans m = tensorFromMatrix (Contravariant,"1") (Covariant,"2") m
65 46
66addT ts = T (dims (head ts)) (fromList $ sumT ts) 47tam = (3><3) [1..9]
48tbm = (3><3) [11..19]
67 49
50ta = tensorFromMatrix (Contravariant,"a1") (Covariant,"a2") tam :: Tensor Double
51tb = tensorFromMatrix (Contravariant,"b1") (Covariant,"b2") tbm :: Tensor Double
68 52
69delta i j | i==j = 1 53delta i j | i==j = 1
70 | otherwise = 0 54 | otherwise = 0
@@ -73,12 +57,6 @@ e i n = fromList [ delta k i | k <- [1..n]]
73 57
74diagl = diag.fromList 58diagl = diag.fromList
75 59
76scalar x = T [] (fromList [x])
77tensorFromVector (tp,nm) v = T {dims = [IdxDesc (dim v) tp nm]
78 , ten = v}
79tensorFromMatrix (tpr,nmr) (tpc,nmc) m = T {dims = [IdxDesc (rows m) tpr nmr,IdxDesc (cols m) tpc nmc]
80 , ten = cdat m}
81
82td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double 60td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double
83 61
84tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double 62tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
@@ -91,7 +69,7 @@ r1 = contraction tt "j" tq "p"
91r1' = contraction' tt "j" tq "p" 69r1' = contraction' tt "j" tq "p"
92 70
93pru = do 71pru = do
94 mapM_ (putStrLn.shdims.dims.normal) (contractions t1 t2) 72 mapM_ (putStrLn.shdims.dims.normal) (possibleContractions t1 t2)
95 let t1 = contraction tt "i" tq "q" 73 let t1 = contraction tt "i" tq "q"
96 print $ normal t1 74 print $ normal t1
97 print $ foldl part t1 [("j",0),("p'",1),("r'",1)] 75 print $ foldl part t1 [("j",0),("p'",1),("r'",1)]
@@ -104,43 +82,24 @@ pru = do
104 let t2 = contraction' tq "q" tt "i" 82 let t2 = contraction' tq "q" tt "i"
105 print $ normal t2 83 print $ normal t2
106 print $ foldl part t2 [("j'",0),("p",1),("r",1)] 84 print $ foldl part t2 [("j'",0),("p",1),("r",1)]
107 85 putStrLn "--------------------------------------------"
108scsig t = scalar (signature (nms t)) `prod` t 86 print $ flatten $ tam <> tbm
109 where nms = map idxName . dims 87 print $ contractions (ta <*> tb <*> kron 3) [("a2","k1'"),("b1'","k2'")]
110 88 print $ contraction ta "a2" tb "b1"
111antisym' t = addT $ map (scsig . flip tridx t) (perms (names t)) 89 print $ normal $ contractions (ta <*> tb) [("a2","b1'")]
112 90 print $ normal $ contraction' ta "a2" tb "b1"
113{- 91 putStrLn "--------------------------------------------"
114 where T d v = t 92 print $ raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
115 t' = T d' v 93 putStrLn "--------------------------------------------"
116 fixdim (T _ v) = T d v 94 print $ foldl' contractionF (leviCivita 4) [y1,y2]
117 d' = [(n,(c,show (pos q))) | (n,(c,q)) <- d] 95 print $ contractions (leviCivita 4 <*> (y1/\y2)) [("1","p'"),("2'","q''")] <*> (scalar (recip $ fact 2))
118 pos n = i where Just i = elemIndex n nms 96 print $ foldl' contractionF (leviCivita 4) [y1,y2,y3,y5]
119 nms = map (snd.snd) d 97 print $ contractions (leviCivita 4 <*> (y1/\y2/\y3/\y5)) [("1","p'"),("2'","q''"),("3'","r''"),("4'","t''")] <*> (scalar (recip $ fact 4))
120-} 98 print $ dim $ ten $ leviCivita 4 <*> (y1/\y2/\y3/\y5)
121 99 print $ innerAT (leviCivita 4) (y1/\y2/\y3/\y5)
122auxrename (T d v) = T d' v 100
123 where d' = [IdxDesc n c (show (pos q)) | IdxDesc n c q <- d] 101
124 pos n = i where Just i = elemIndex n nms 102y5 = vector "t" [0,1,-2,0]
125 nms = map idxName d
126
127antisym t = T (dims t) (ten (antisym' (auxrename t)))
128
129
130norper t = prod t (scalar (recip $ fromIntegral $ product [1 .. length (dims t)]))
131antinorper t = prod t (scalar (fromIntegral $ product [1 .. length (dims t)]))
132
133
134tvector n v = tensorFromVector (Contravariant,n) v
135tcovector n v = tensorFromVector (Covariant,n) v
136
137vector n v = tvector n (fromList v) :: Tensor Double
138
139wedge a b = antisym (prod (norper a) (norper b))
140
141a /\ b = wedge a b
142
143a <*> b = normal $ prod a b
144 103
145u = vector "p" [1,1,0] 104u = vector "p" [1,1,0]
146v = vector "q" [0,1,1] 105v = vector "q" [0,1,1]
@@ -149,35 +108,15 @@ w = vector "r" [1,0,1]
149uv = u /\ v 108uv = u /\ v
150uw = u /\ w 109uw = u /\ w
151 110
152normAT t = sqrt $ innerAT t t
153
154innerAT t1 t2 = dot (ten t1) (ten t2) / fromIntegral (fact $ length $ dims t1)
155
156det m = product $ toList s where (_,s,_) = svdR' m
157
158fact n = product [1..n]
159 111
160l1 = vector "p" [0,0,0,1] 112l1 = vector "p" [0,0,0,1]
161l2 = vector "q" [1,0,0,1] 113l2 = vector "q" [1,0,0,1]
162l3 = vector "r" [0,1,0,1] 114l3 = vector "r" [0,1,0,1]
163 115
164leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1,2..]) (toRows (ident n))
165
166contractionF t1 t2 = contraction t1 n1 t2 n2
167 where n1 = fn t1
168 n2 = fn t2
169 fn = idxName . head . dims
170
171
172dualV vs = foldl' contractionF (leviCivita n) vs
173 where n = idxDim . head . dims . head $ vs
174
175
176dual1 = foldl' contractionF (leviCivita 3) [u,v] 116dual1 = foldl' contractionF (leviCivita 3) [u,v]
177dual2 = foldl' contractionF (leviCivita 3) [u,v,w] 117dual2 = foldl' contractionF (leviCivita 3) [u,v,w]
178 118
179 119
180contract1b t (n1,n2) = contract1 t n1 n2
181 120
182dual1' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v)) [("1","p'"),("2'","q''")]) (scalar (recip $ fact 2)) 121dual1' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v)) [("1","p'"),("2'","q''")]) (scalar (recip $ fact 2))
183dual2' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v /\ w)) [("1","p'"),("2'","q''"),("3'","r''")]) (scalar (recip $ fact 3)) 122dual2' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v /\ w)) [("1","p'"),("2'","q''"),("3'","r''")]) (scalar (recip $ fact 3))
@@ -188,17 +127,6 @@ x2 = vector "q" [2,2,2]
188x3 = vector "r" [-3,-1,-1] 127x3 = vector "r" [-3,-1,-1]
189x4 = vector "s" [12,0,3] 128x4 = vector "s" [12,0,3]
190 129
191raise (T d v) = T (map raise' d) v
192 where raise' idx@IdxDesc {idxType = Covariant } = idx {idxType = Contravariant}
193 raise' idx@IdxDesc {idxType = Contravariant } = idx {idxType = Covariant}
194
195dualMV t = prod (foldl' contract1b (lc <*> t) ds) (scalar (recip $ fromIntegral $ fact (length ds)))
196 where
197 lc = leviCivita n
198 nms1 = map idxName (dims lc)
199 nms2 = map ((++"'").idxName) (dims t)
200 ds = zip nms1 nms2
201 n = idxDim . head . dims $ t
202 130
203-- intersection of two lines :-) 131-- intersection of two lines :-)
204-- > raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4] 132-- > raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
@@ -214,17 +142,6 @@ y4 = vector "s" [12,0,0,3]
214-- scalar 0.0 142-- scalar 0.0
215-- it seems that the sum of ranks must be greater than n :( 143-- it seems that the sum of ranks must be greater than n :(
216 144
217asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
218
219partF t i = part t (name,i) where name = idxName . head . dims $ t
220
221--partL = foldl' partF
222
223niceAS t = filter ((/=0.0).fst) $ zip vals base
224 where vals = map ((`at` 0).ten.foldl' partF t) (map (map pred) base)
225 base = asBase r n
226 r = length (dims t)
227 n = idxDim . head . dims $ t
228 145
229z1 = vector "p" [0,0,0,1] 146z1 = vector "p" [0,0,0,1]
230z2 = vector "q" [1,0,0,1] 147z2 = vector "q" [1,0,0,1]
diff --git a/examples/tests.hs b/examples/tests.hs
index 5acfa9d..eb8f50d 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -20,7 +20,7 @@ import Test.HUnit hiding ((~:))
20import Complex 20import Complex
21import LinearAlgebra.Algorithms 21import LinearAlgebra.Algorithms
22import GSL.Matrix 22import GSL.Matrix
23import GSL.Compat hiding ((<>)) 23import GSL.Compat hiding ((<>),constant)
24 24
25dist :: (Normed t, Num t) => t -> t -> Double 25dist :: (Normed t, Num t) => t -> t -> Double
26dist a b = norm (a-b) 26dist a b = norm (a-b)
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
index c4faf49..8296935 100644
--- a/lib/Data/Packed/Internal/Tensor.hs
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -18,7 +18,8 @@ import Data.Packed.Internal.Common
18import Data.Packed.Internal.Vector 18import Data.Packed.Internal.Vector
19import Data.Packed.Internal.Matrix 19import Data.Packed.Internal.Matrix
20import Foreign.Storable 20import Foreign.Storable
21import Data.List(sort,elemIndex,nub) 21import Data.List(sort,elemIndex,nub,foldl1')
22import GSL.Vector
22 23
23data IdxType = Covariant | Contravariant deriving (Show,Eq) 24data IdxType = Covariant | Contravariant deriving (Show,Eq)
24 25
@@ -79,10 +80,10 @@ concatRename l1 l2 = l1 ++ map ren l2 where
79 ren idx = if {- s `elem` fs -} True then idx {idxName = idxName idx ++ "'"} else idx 80 ren idx = if {- s `elem` fs -} True then idx {idxName = idxName idx ++ "'"} else idx
80 fs = map idxName l1 81 fs = map idxName l1
81 82
82prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t 83--prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t
83prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer' v1 v2) 84prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer' v1 v2)
84 85
85contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t 86--contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
86contraction t1 n1 t2 n2 = 87contraction t1 n1 t2 n2 =
87 if compatIdx t1 n1 t2 n2 88 if compatIdx t1 n1 t2 n2
88 then T (concatRename (tail d1) (tail d2)) (cdat m) 89 then T (concatRename (tail d1) (tail d2)) (cdat m)
@@ -91,16 +92,27 @@ contraction t1 n1 t2 n2 =
91 (d2,m2) = putFirstIdx n2 t2 92 (d2,m2) = putFirstIdx n2 t2
92 m = multiply RowMajor (trans m1) m2 93 m = multiply RowMajor (trans m1) m2
93 94
94sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t] 95--sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t]
95sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls) 96--sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls)
97--addT ts = T (dims (head ts)) (fromList $ sumT ts)
96 98
97contract1 :: (Num t, Enum t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t 99liftTensor f (T d v) = T d (f v)
98contract1 t name1 name2 = T d $ fromList $ sumT y 100
101liftTensor2 f (T d1 v1) (T d2 v2) | compat d1 d2 = T d1 (f v1 v2)
102 | otherwise = error "liftTensor2 with incompatible tensors"
103 where compat a b = length a == length b
104
105
106a |+| b = liftTensor2 add a b
107addT l = foldl1' (|+|) l
108
109--contract1 :: (Num t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t
110contract1 t name1 name2 = addT y
99 where d = dims (head y) 111 where d = dims (head y)
100 x = (map (flip parts name2) (parts t name1)) 112 x = (map (flip parts name2) (parts t name1))
101 y = map head $ zipWith drop [0..] x 113 y = map head $ zipWith drop [0..] x
102 114
103contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t 115--contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
104contraction' t1 n1 t2 n2 = 116contraction' t1 n1 t2 n2 =
105 if compatIdx t1 n1 t2 n2 117 if compatIdx t1 n1 t2 n2
106 then contract1 (prod t1 t2) n1 (n2++"'") 118 then contract1 (prod t1 t2) n1 (n2++"'")
@@ -130,8 +142,8 @@ names t = sort $ map idxName (dims t)
130normal :: (Field t) => Tensor t -> Tensor t 142normal :: (Field t) => Tensor t -> Tensor t
131normal t = tridx (names t) t 143normal t = tridx (names t) t
132 144
133contractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t] 145possibleContractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t]
134contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ] 146possibleContractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
135 147
136-- sent to Haskell-Cafe by Sebastian Sylvan 148-- sent to Haskell-Cafe by Sebastian Sylvan
137perms :: [t] -> [[t]] 149perms :: [t] -> [[t]]
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 36bf32e..2033dc7 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -87,14 +87,14 @@ ident n = diag (constant 1 n)
87r >< c = f where 87r >< c = f where
88 f l | dim v == r*c = matrixFromVector RowMajor c v 88 f l | dim v == r*c = matrixFromVector RowMajor c v
89 | otherwise = error $ "inconsistent list size = " 89 | otherwise = error $ "inconsistent list size = "
90 ++show (dim v) ++"in ("++show r++"><"++show c++")" 90 ++show (dim v) ++" in ("++show r++"><"++show c++")"
91 where v = fromList l 91 where v = fromList l
92 92
93(>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a 93(>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a
94r >|< c = f where 94r >|< c = f where
95 f l | dim v == r*c = matrixFromVector ColumnMajor c v 95 f l | dim v == r*c = matrixFromVector ColumnMajor c v
96 | otherwise = error $ "inconsistent list size = " 96 | otherwise = error $ "inconsistent list size = "
97 ++show (dim v) ++"in ("++show r++"><"++show c++")" 97 ++show (dim v) ++" in ("++show r++"><"++show c++")"
98 where v = fromList l 98 where v = fromList l
99 99
100---------------------------------------------------------------- 100----------------------------------------------------------------
diff --git a/lib/Data/Packed/Tensor.hs b/lib/Data/Packed/Tensor.hs
index 75a9288..68ce9a5 100644
--- a/lib/Data/Packed/Tensor.hs
+++ b/lib/Data/Packed/Tensor.hs
@@ -12,9 +12,85 @@
12-- 12--
13----------------------------------------------------------------------------- 13-----------------------------------------------------------------------------
14 14
15module Data.Packed.Tensor ( 15module Data.Packed.Tensor where
16
17) where
18 16
17import Data.Packed.Matrix
19import Data.Packed.Internal 18import Data.Packed.Internal
20import Complex 19import Complex
20import Data.List(transpose,intersperse,sort,elemIndex,nub,foldl',foldl1')
21
22scalar x = T [] (fromList [x])
23tensorFromVector (tp,nm) v = T {dims = [IdxDesc (dim v) tp nm]
24 , ten = v}
25tensorFromMatrix (tpr,nmr) (tpc,nmc) m = T {dims = [IdxDesc (rows m) tpr nmr,IdxDesc (cols m) tpc nmc]
26 , ten = cdat m}
27
28scsig t = scalar (signature (nms t)) `prod` t
29 where nms = map idxName . dims
30
31antisym' t = addT $ map (scsig . flip tridx t) (perms (names t))
32
33
34auxrename (T d v) = T d' v
35 where d' = [IdxDesc n c (show (pos q)) | IdxDesc n c q <- d]
36 pos n = i where Just i = elemIndex n nms
37 nms = map idxName d
38
39antisym t = T (dims t) (ten (antisym' (auxrename t)))
40
41
42norper t = prod t (scalar (recip $ fromIntegral $ product [1 .. length (dims t)]))
43antinorper t = prod t (scalar (fromIntegral $ product [1 .. length (dims t)]))
44
45
46tvector n v = tensorFromVector (Contravariant,n) v
47tcovector n v = tensorFromVector (Covariant,n) v
48
49wedge a b = antisym (prod (norper a) (norper b))
50
51a /\ b = wedge a b
52
53a <*> b = normal $ prod a b
54
55normAT t = sqrt $ innerAT t t
56
57innerAT t1 t2 = dot (ten t1) (ten t2) / fromIntegral (fact $ length $ dims t1)
58
59fact n = product [1..n]
60
61leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1,2..]) (toRows (ident n))
62
63contractionF t1 t2 = contraction t1 n1 t2 n2
64 where n1 = fn t1
65 n2 = fn t2
66 fn = idxName . head . dims
67
68
69dualV vs = foldl' contractionF (leviCivita n) vs
70 where n = idxDim . head . dims . head $ vs
71
72raise (T d v) = T (map raise' d) v
73 where raise' idx@IdxDesc {idxType = Covariant } = idx {idxType = Contravariant}
74 raise' idx@IdxDesc {idxType = Contravariant } = idx {idxType = Covariant}
75
76dualMV t = prod (foldl' contract1b (lc <*> t) ds) (scalar (recip $ fromIntegral $ fact (length ds)))
77 where
78 lc = leviCivita n
79 nms1 = map idxName (dims lc)
80 nms2 = map ((++"'").idxName) (dims t)
81 ds = zip nms1 nms2
82 n = idxDim . head . dims $ t
83
84contract1b t (n1,n2) = contract1 t n1 n2
85
86contractions t pairs = foldl' contract1b t pairs
87
88asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
89
90partF t i = part t (name,i) where name = idxName . head . dims $ t
91
92niceAS t = filter ((/=0.0).fst) $ zip vals base
93 where vals = map ((`at` 0).ten.foldl' partF t) (map (map pred) base)
94 base = asBase r n
95 r = length (dims t)
96 n = idxDim . head . dims $ t
diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs
index a074254..0b3c3a9 100644
--- a/lib/GSL/Vector.hs
+++ b/lib/GSL/Vector.hs
@@ -21,7 +21,9 @@ module GSL.Vector (
21 scale, addConstant, add, sub, mul, 21 scale, addConstant, add, sub, mul,
22) where 22) where
23 23
24import Data.Packed.Internal 24import Data.Packed.Internal.Common
25import Data.Packed.Internal.Vector
26
25import Complex 27import Complex
26import Foreign 28import Foreign
27 29