From aa14e6615533e7bd5e2b15acdc3ec76afbe1aac4 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Thu, 21 Jun 2007 08:39:22 +0000
Subject: working with tensors

---
 examples/pru.hs                    | 82 ++++++++++++++++++++++++++++++++++++--
 examples/tests.hs                  |  2 +-
 lib/Data/Packed/Internal/Common.hs |  4 +-
 lib/Data/Packed/Internal/Tensor.hs | 17 +++++++-
 4 files changed, 98 insertions(+), 7 deletions(-)

diff --git a/examples/pru.hs b/examples/pru.hs
index 31e5213..4a5104b 100644
--- a/examples/pru.hs
+++ b/examples/pru.hs
@@ -11,7 +11,7 @@ import LAPACK
 
 import Complex
 import Numeric(showGFloat)
-import Data.List(transpose,intersperse,sort)
+import Data.List(transpose,intersperse,sort,elemIndex,nub,foldl',foldl1')
 import Foreign.Storable
 
 
@@ -61,7 +61,7 @@ t2 = T [(4,(Covariant,"p")),(4,(Contravariant,"q"))] $ fromList [1..16::Double]
 
 
 
-
+addT ts = T (dims (head ts)) (fromList $ sumT ts)
 
 
 delta i j | i==j      = 1
@@ -71,7 +71,7 @@ e i n = fromList [ delta k i | k <- [1..n]]
 
 diagl = diag.fromList
 
-
+scalar x = T [] (fromList [x])
 tensorFromVector idx v = T {dims = [(dim v,idx)], ten = v}
 tensorFromMatrix idxr idxc m = T {dims = [(rows m,idxr),(cols m,idxc)], ten = cdat m}
 
@@ -100,3 +100,79 @@ pru = do
     print $ normal t2
     print $ foldl part t2 [("j'",0),("p",1),("r",1)]
 
+scsig t = scalar (signature (nms t)) `prod` t
+    where nms = map (snd.snd) . dims
+
+antisym' t = addT $ map (scsig . flip tridx t) (perms (names t))
+
+{-
+   where T d v = t
+          t' = T d' v
+          fixdim (T _ v) = T d v
+          d' = [(n,(c,show (pos q))) | (n,(c,q)) <- d]
+          pos n = i where Just i = elemIndex n nms
+          nms = map (snd.snd) d
+-}
+
+auxrename (T d v) = T d' v
+    where d' = [(n,(c,show (pos q))) | (n,(c,q)) <- d]
+          pos n = i where Just i = elemIndex n nms
+          nms = map (snd.snd) d
+
+antisym t = T (dims t) (ten (antisym' (auxrename t)))
+
+
+norper t = prod t (scalar (recip $ fromIntegral $ product [1 .. length (dims t)]))
+antinorper t = prod t (scalar (fromIntegral $ product [1 .. length (dims t)]))
+
+
+tvector n v = tensorFromVector (Contravariant,n) v
+tcovector n v = tensorFromVector (Covariant,n) v
+
+vector n v = tvector n (fromList v) :: Tensor Double
+
+wedge a b = antisym (prod (norper a) (norper b))
+
+a /\ b = wedge a b
+
+a <*> b = normal $ prod a b
+
+u = vector "p" [1,1,0]
+v = vector "q" [0,1,1]
+w = vector "r" [1,0,1]
+
+uv = u /\ v
+uw = u /\ w
+
+normAT t = sqrt $ innerAT t t
+
+innerAT t1 t2 = dot (ten t1) (ten t2) / fromIntegral (fact $ length $ dims t1)
+
+det m = product $ toList s where (_,s,_) = svdR' m
+
+fact n = product [1..n]
+
+l1 = vector "p" [0,0,0,1]
+l2 = vector "q" [1,0,0,1]
+l3 = vector "r" [0,1,0,1]
+
+leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1..]) (toRows (ident n))
+
+contractionF t1 t2 = contraction t1 n1 t2 n2
+    where n1 = fn t1
+          n2 = fn t2
+          fn = snd . snd . head . dims
+
+
+dual vs = foldl' contractionF (leviCivita n) vs
+    where n = fst . head . dims . head $ vs
+
+
+dual1 = foldl' contractionF (leviCivita 3) [u,v]
+dual2 = foldl' contractionF (leviCivita 3) [u,v,w]
+
+
+contract1b t (n1,n2) = contract1 t n1 n2
+
+dual1' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v)) [("1","p'"),("2'","q''")]) (scalar (recip $ fact 2))
+dual2' = prod (foldl' contract1b ((leviCivita 3) <*> (u /\ v /\ w)) [("1","p'"),("2'","q''"),("3'","r''")]) (scalar (recip $ fact 3))
\ No newline at end of file
diff --git a/examples/tests.hs b/examples/tests.hs
index 22c4674..4adc104 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -280,7 +280,7 @@ arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u)
           ~|~ vectorMapR Sin u
 
 
---arit3 (PairV u v) = vectorMap Sin . VectorMap Cos
+-- arit3 (PairV u v) =
 
 ---------------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index acefe92..91985f7 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -74,8 +74,8 @@ check msg ls f = do
     mapM_ (touchForeignPtr . fptr) ls
     return ()
 
-class (Storable a, Typeable a) => Field a where
-instance (Storable a, Typeable a) => Field a where
+class (Storable a, Typeable a) => Field a
+instance (Storable a, Typeable a) => Field a
 
 isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
 isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
index 123270d..27fce6a 100644
--- a/lib/Data/Packed/Internal/Tensor.hs
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -18,7 +18,7 @@ import Data.Packed.Internal
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Foreign.Storable
-import Data.List(sort)
+import Data.List(sort,elemIndex,nub)
 
 data IdxTp = Covariant | Contravariant deriving (Show,Eq)
 
@@ -107,3 +107,18 @@ normal t = tridx (names t) t
 
 contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
 
+-- sent to Haskell-Cafe by Sebastian Sylvan
+perms [x] = [[x]]
+perms xs = [y:ps | (y,ys) <- selections xs , ps <- perms ys]
+selections []     = []
+selections (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- selections xs]
+
+
+interchanges ls = sum (map (count ls) ls)
+    where count l p = n
+              where Just pel = elemIndex p l
+                    n = length $ filter (>p) $ take pel l
+
+signature l | length (nub l) < length l =  0
+            | even (interchanges l)     =  1
+            | otherwise                 = -1
-- 
cgit v1.2.3