summaryrefslogtreecommitdiff
path: root/examples/pru.hs
diff options
context:
space:
mode:
Diffstat (limited to 'examples/pru.hs')
-rw-r--r--examples/pru.hs149
1 files changed, 0 insertions, 149 deletions
diff --git a/examples/pru.hs b/examples/pru.hs
deleted file mode 100644
index 71f33bf..0000000
--- a/examples/pru.hs
+++ /dev/null
@@ -1,149 +0,0 @@
1--{-# OPTIONS_GHC #-}
2--module Main where
3{-
4import Data.Packed.Internal
5import Data.Packed.Internal.Vector
6import Data.Packed.Internal.Matrix
7import Data.Packed.Tensor
8import Data.Packed.Matrix
9import GSL.Vector
10import LAPACK
11import Data.List(foldl')
12
13import Complex
14import Numeric(showGFloat)
15
16import Foreign.Storable
17-}
18
19import GSL
20import Data.List(foldl', foldl1')
21import Data.Packed.Internal.Tensor
22import Data.Packed.Tensor
23import Complex
24
25
26main = do
27 print $ foldl part t [("p",1),("q",0),("r",2)]
28 print $ foldl part t [("p",1),("r",2),("q",0)]
29 print $ foldl part t $ reverse [("p",1),("r",2),("q",0)]
30 pru
31
32t = T [IdxDesc 4 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 3 Contravariant "r"]
33 $ fromList [1..24::Double]
34
35
36t1 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q" ,IdxDesc 2 Covariant "r"]
37 $ fromList [1..32::Double]
38t2 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q"]
39 $ fromList [1..16::Double]
40
41vector n v = tvector n (fromList v) :: Tensor Double
42
43kron n = tensorFromMatrix (Contravariant,"k1") (Covariant,"k2") (ident n)
44
45tensorFromTrans m = tensorFromMatrix (Contravariant,"1") (Covariant,"2") m
46
47tam = (3><3) [1..9]
48tbm = (3><3) [11..19]
49
50ta = tensorFromMatrix (Contravariant,"a1") (Covariant,"a2") tam :: Tensor Double
51tb = tensorFromMatrix (Contravariant,"b1") (Covariant,"b2") tbm :: Tensor Double
52
53delta i j | i==j = 1
54 | otherwise = 0
55
56e i n = fromList [ delta k i | k <- [1..n]]
57
58diagl = diag.fromList
59
60td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double
61
62tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
63tt = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
64
65tq = T [IdxDesc 3 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 2 Covariant "r"]
66 $ fromList [11 .. 22] :: Tensor Double
67
68r1 = contraction tt "j" tq "p"
69r1' = contraction' tt "j" tq "p"
70
71pru = do
72 mapM_ (putStrLn.shdims.dims.normal) (possibleContractions t1 t2)
73 let t1 = contraction tt "i" tq "q"
74 print $ normal t1
75 print $ foldl part t1 [("j",0),("p'",1),("r'",1)]
76 let t2 = contraction' tt "i" tq "q"
77 print $ normal t2
78 print $ foldl part t2 [("j",0),("p'",1),("r'",1)]
79 let t1 = contraction tq "q" tt "i"
80 print $ normal t1
81 print $ foldl part t1 [("j'",0),("p",1),("r",1)]
82 let t2 = contraction' tq "q" tt "i"
83 print $ normal t2
84 print $ foldl part t2 [("j'",0),("p",1),("r",1)]
85 putStrLn "--------------------------------------------"
86 print $ flatten $ tam <> tbm
87 print $ contractions (ta <*> tb <*> kron 3) [("a2","k1'"),("b1'","k2'")]
88 print $ contraction ta "a2" tb "b1"
89 print $ normal $ contractions (ta <*> tb) [("a2","b1'")]
90 print $ normal $ contraction' ta "a2" tb "b1"
91 putStrLn "--------------------------------------------"
92 print $ raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
93 putStrLn "--------------------------------------------"
94 print $ foldl' contractionF (leviCivita 4) [y1,y2]
95 print $ contractions (leviCivita 4 <*> (y1/\y2)) [("1","p'"),("2'","q''")] <*> (scalar (recip $ fact 2))
96 print $ foldl' contractionF (leviCivita 4) [y1,y2,y3,y5]
97 print $ contractions (leviCivita 4 <*> (y1/\y2/\y3/\y5)) [("1","p'"),("2'","q''"),("3'","r''"),("4'","t''")] <*> (scalar (recip $ fact 4))
98 print $ dim $ ten $ leviCivita 4 <*> (y1/\y2/\y3/\y5)
99 print $ innerAT (leviCivita 4) (y1/\y2/\y3/\y5)
100
101
102y5 = vector "t" [0,1,-2,0]
103
104u = vector "p" [1,1,0]
105v = vector "q" [0,1,1]
106w = vector "r" [1,0,1]
107
108uv = u /\ v
109uw = u /\ w
110
111
112l1 = vector "p" [0,0,0,1]
113l2 = vector "q" [1,0,0,1]
114l3 = vector "r" [0,1,0,1]
115
116dual1 = foldl' contractionF (leviCivita 3) [u,v]
117dual2 = foldl' contractionF (leviCivita 3) [u,v,w]
118
119
120
121dual1' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v)) [("1","p'"),("2'","q''")]) (scalar (recip $ fact 2))
122dual2' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v /\ w)) [("1","p'"),("2'","q''"),("3'","r''")]) (scalar (recip $ fact 3))
123
124
125x1 = vector "p" [0,0,1]
126x2 = vector "q" [2,2,2]
127x3 = vector "r" [-3,-1,-1]
128x4 = vector "s" [12,0,3]
129
130
131-- intersection of two lines :-)
132-- > raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
133--(3'^[3]) [24.0,24.0,12.0]
134
135y1 = vector "p" [0,0,0,1]
136y2 = vector "q" [2,2,0,2]
137y3 = vector "r" [-3,-1,0,-1]
138y4 = vector "s" [12,0,0,3]
139
140-- why not in R^4?
141-- > raise $ dualMV $ raise $ dualMV (y1/\y2) /\ dualV [y3,y4]
142-- scalar 0.0
143-- it seems that the sum of ranks must be greater than n :(
144
145
146z1 = vector "p" [0,0,0,1]
147z2 = vector "q" [1,0,0,1]
148z3 = vector "r" [0,1,0,1]
149z4 = vector "s" [0,0,1,1]