tensor refactorization

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
 toComplexM = uncurry $ liftMatrix2 (curry toComplex)
 
 infixl 2 =~=
-a =~= b = pnorm PNorm1 (flatten (a - b)) < 1E-6
+a =~= b = pnorm 1 (flatten (a - b)) < 1E-6
 
 randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed 
 
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 @@
 --{-# OPTIONS_GHC  #-}
 --module Main where
-
+{-
 import Data.Packed.Internal
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
-import Data.Packed.Internal.Tensor
+import Data.Packed.Tensor
 import Data.Packed.Matrix
 import GSL.Vector
 import LAPACK
+import Data.List(foldl')
 
 import Complex
 import Numeric(showGFloat)
-import Data.List(transpose,intersperse,sort,elemIndex,nub,foldl',foldl1')
-import Foreign.Storable
-
-
-vr = fromList [1..15::Double]
-vc = fromList (map (\x->x :+ (x+1)) [1..15::Double])
-
-mi = (2 >< 3) [1 .. 6::Int]
-mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)]
-
-ac = (2><3) [1 .. 6::Double]
-bc = (3><4) [7 .. 18::Double]
-
-af = (2>|<3) [1,4,2,5,3,6::Double]
-bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
-
-
-a |=| b = rows a == rows b &&
-          cols a == cols b &&
-          toList (cdat a) == toList (cdat b)
-
-mulC a b = multiply RowMajor a b
-mulF a b = multiply ColumnMajor a b
-
-cc = mulC ac bf
-cf = mulF af bc
 
-r = mulC cc (trans cf)
-
-rd = (2><2)
- [ 27736.0,  65356.0
- , 65356.0, 154006.0 ]
+import Foreign.Storable
+-}
 
+import GSL
+import Data.List(foldl', foldl1')
+import Data.Packed.Internal.Tensor
+import Data.Packed.Tensor
+import Complex
 
 
 main = do
-    print $ r |=| rd
     print $ foldl part t [("p",1),("q",0),("r",2)]
     print $ foldl part t [("p",1),("r",2),("q",0)]
     print $ foldl part t $ reverse [("p",1),("r",2),("q",0)]
+    pru
 
 t = T [IdxDesc 4 Covariant "p",IdxDesc 2 Covariant "q" ,IdxDesc 3 Contravariant "r"]
     $ fromList [1..24::Double]
@@ -61,10 +38,17 @@ t1 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q" ,IdxDesc 2 Covariant
 t2 = T [IdxDesc 4 Covariant "p",IdxDesc 4 Contravariant "q"] 
     $ fromList [1..16::Double]
 
+vector n v = tvector n (fromList v) :: Tensor Double
+
+kron n = tensorFromMatrix (Contravariant,"k1") (Covariant,"k2") (ident n)
 
+tensorFromTrans m = tensorFromMatrix (Contravariant,"1") (Covariant,"2") m
 
-addT ts = T (dims (head ts)) (fromList $ sumT ts)
+tam = (3><3) [1..9]
+tbm = (3><3) [11..19]
 
+ta = tensorFromMatrix (Contravariant,"a1") (Covariant,"a2") tam :: Tensor Double
+tb = tensorFromMatrix (Contravariant,"b1") (Covariant,"b2") tbm :: Tensor Double
 
 delta i j | i==j      = 1
           | otherwise = 0
@@ -73,12 +57,6 @@ e i n = fromList [ delta k i | k <- [1..n]]
 
 diagl = diag.fromList
 
-scalar x = T [] (fromList [x])
-tensorFromVector (tp,nm) v = T {dims = [IdxDesc (dim v) tp nm]
-                                       , ten = v}
-tensorFromMatrix (tpr,nmr) (tpc,nmc) m = T {dims = [IdxDesc (rows m) tpr nmr,IdxDesc (cols m) tpc nmc]
-                                           , ten = cdat m}
-
 td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double
 
 tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
@@ -91,7 +69,7 @@ r1 = contraction tt "j" tq "p"
 r1' = contraction' tt "j" tq "p"
 
 pru = do
-    mapM_ (putStrLn.shdims.dims.normal) (contractions t1 t2)
+    mapM_ (putStrLn.shdims.dims.normal) (possibleContractions t1 t2)
     let t1 = contraction tt "i" tq "q"
     print $ normal t1
     print $ foldl part t1 [("j",0),("p'",1),("r'",1)]
@@ -104,43 +82,24 @@ pru = do
     let t2 = contraction' tq "q" tt "i"
     print $ normal t2
     print $ foldl part t2 [("j'",0),("p",1),("r",1)]
-
-scsig t = scalar (signature (nms t)) `prod` t
-    where nms = map idxName . 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' = [IdxDesc n c (show (pos q)) | IdxDesc n c q <- d]
-          pos n = i where Just i = elemIndex n nms
-          nms = map idxName 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
+    putStrLn "--------------------------------------------"
+    print $ flatten $ tam <> tbm
+    print $ contractions (ta <*> tb <*> kron 3) [("a2","k1'"),("b1'","k2'")]
+    print $ contraction ta "a2" tb "b1"
+    print $ normal $ contractions (ta <*> tb) [("a2","b1'")]
+    print $ normal $ contraction' ta "a2" tb "b1"
+    putStrLn "--------------------------------------------"
+    print $ raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
+    putStrLn "--------------------------------------------"
+    print $ foldl' contractionF (leviCivita 4) [y1,y2]
+    print $ contractions (leviCivita 4 <*> (y1/\y2)) [("1","p'"),("2'","q''")] <*> (scalar (recip $ fact 2))
+    print $ foldl' contractionF (leviCivita 4) [y1,y2,y3,y5]
+    print $ contractions (leviCivita 4 <*> (y1/\y2/\y3/\y5)) [("1","p'"),("2'","q''"),("3'","r''"),("4'","t''")] <*> (scalar (recip $ fact 4))
+    print $ dim $ ten $ leviCivita 4 <*> (y1/\y2/\y3/\y5)
+    print $ innerAT (leviCivita 4) (y1/\y2/\y3/\y5)
+
+
+y5 = vector "t" [0,1,-2,0]
 
 u = vector "p" [1,1,0]
 v = vector "q" [0,1,1]
@@ -149,35 +108,15 @@ 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,2..]) (toRows (ident n))
-
-contractionF t1 t2 = contraction t1 n1 t2 n2
-    where n1 = fn t1
-          n2 = fn t2
-          fn = idxName . head . dims
-
-
-dualV vs = foldl' contractionF (leviCivita n) vs
-    where n = idxDim . 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))
@@ -188,17 +127,6 @@ x2 = vector "q" [2,2,2]
 x3 = vector "r" [-3,-1,-1]
 x4 = vector "s" [12,0,3]
 
-raise (T d v) = T (map raise' d) v
-    where raise' idx@IdxDesc {idxType = Covariant } = idx {idxType = Contravariant}
-          raise' idx@IdxDesc {idxType = Contravariant } = idx {idxType = Covariant}
-
-dualMV t = prod (foldl' contract1b (lc <*> t) ds) (scalar (recip $ fromIntegral $ fact (length ds)))
-    where
-        lc = leviCivita n
-        nms1 = map idxName (dims lc)
-        nms2 = map ((++"'").idxName) (dims t)
-        ds = zip nms1 nms2
-        n = idxDim . head . dims $ t
 
 -- intersection of two lines :-)
 -- > raise $ dualMV $ raise $ dualMV (x1/\x2) /\ dualV [x3,x4]
@@ -214,17 +142,6 @@ y4 = vector "s" [12,0,0,3]
 -- scalar 0.0
 -- it seems that the sum of ranks must be greater than n :(
 
-asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
-
-partF t i = part t (name,i) where name = idxName . head . dims $ t
-
---partL = foldl' partF
-
-niceAS t = filter ((/=0.0).fst) $ zip vals base
-    where vals = map ((`at` 0).ten.foldl' partF t) (map (map pred) base)
-          base = asBase r n
-          r = length (dims t)
-          n = idxDim . head . dims $ t
 
 z1 = vector "p" [0,0,0,1]
 z2 = 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 ((~:))
 import Complex
 import LinearAlgebra.Algorithms
 import GSL.Matrix
-import GSL.Compat hiding ((<>))
+import GSL.Compat hiding ((<>),constant)
 
 dist :: (Normed t, Num t) => t -> t -> Double
 dist 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
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Foreign.Storable
-import Data.List(sort,elemIndex,nub)
+import Data.List(sort,elemIndex,nub,foldl1')
+import GSL.Vector
 
 data IdxType = Covariant | Contravariant deriving (Show,Eq)
 
@@ -79,10 +80,10 @@ concatRename l1 l2 = l1 ++ map ren l2 where
     ren idx = if {- s `elem` fs -} True then idx {idxName = idxName idx ++ "'"} else idx
     fs = map idxName l1
 
-prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t
+--prod :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t
 prod (T d1 v1) (T d2 v2) = T (concatRename d1 d2) (outer' v1 v2)
 
-contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
+--contraction :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
 contraction t1 n1 t2 n2 =
     if compatIdx t1 n1 t2 n2
         then T (concatRename (tail d1) (tail d2)) (cdat m)
@@ -91,16 +92,27 @@ contraction t1 n1 t2 n2 =
         (d2,m2) = putFirstIdx n2 t2
         m = multiply RowMajor (trans m1) m2
 
-sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t]
-sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls)
+--sumT :: (Storable t, Enum t, Num t) => [Tensor t] -> [t]
+--sumT ls = foldl (zipWith (+)) [0,0..] (map (toList.ten) ls)
+--addT ts = T (dims (head ts)) (fromList $ sumT ts)
 
-contract1 :: (Num t, Enum t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t
-contract1 t name1 name2 = T d $ fromList $ sumT y
+liftTensor f (T d v) = T d (f v)
+
+liftTensor2 f (T d1 v1) (T d2 v2) | compat d1 d2 = T d1 (f v1 v2)
+                                  | otherwise = error "liftTensor2 with incompatible tensors"
+    where compat a b = length a == length b
+
+
+a |+| b = liftTensor2 add a b
+addT l = foldl1' (|+|) l
+
+--contract1 :: (Num t, Field t) => Tensor t -> IdxName -> IdxName -> Tensor t
+contract1 t name1 name2 = addT y
     where d = dims (head y)
           x = (map (flip parts name2) (parts t name1))
           y = map head $ zipWith drop [0..] x
 
-contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
+--contraction' :: (Field t, Enum t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t
 contraction' t1 n1 t2 n2 =
     if compatIdx t1 n1 t2 n2
         then contract1 (prod t1 t2) n1 (n2++"'")
@@ -130,8 +142,8 @@ names t = sort $ map idxName (dims t)
 normal :: (Field t) => Tensor t -> Tensor t
 normal t = tridx (names t) t
 
-contractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t]
-contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
+possibleContractions :: (Num t, Field t) => Tensor t -> Tensor t -> [Tensor t]
+possibleContractions 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 :: [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)
 r >< c = f where
     f l | dim v == r*c = matrixFromVector RowMajor c v
         | otherwise    = error $ "inconsistent list size = "
-                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+                                 ++show (dim v) ++" in ("++show r++"><"++show c++")"
         where v = fromList l
 
 (>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a
 r >|< c = f where
     f l | dim v == r*c = matrixFromVector ColumnMajor c v
         | otherwise    = error $ "inconsistent list size = "
-                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+                                 ++show (dim v) ++" in ("++show r++"><"++show c++")"
         where v = fromList l
 
 ----------------------------------------------------------------
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 @@
 --
 -----------------------------------------------------------------------------
 
-module Data.Packed.Tensor (
-
-) where
+module Data.Packed.Tensor where
 
+import Data.Packed.Matrix
 import Data.Packed.Internal
 import Complex
+import Data.List(transpose,intersperse,sort,elemIndex,nub,foldl',foldl1')
+
+scalar x = T [] (fromList [x])
+tensorFromVector (tp,nm) v = T {dims = [IdxDesc (dim v) tp nm]
+                                       , ten = v}
+tensorFromMatrix (tpr,nmr) (tpc,nmc) m = T {dims = [IdxDesc (rows m) tpr nmr,IdxDesc (cols m) tpc nmc]
+                                           , ten = cdat m}
+
+scsig t = scalar (signature (nms t)) `prod` t
+    where nms = map idxName . dims
+
+antisym' t = addT $ map (scsig . flip tridx t) (perms (names t))
+
+
+auxrename (T d v) = T d' v
+    where d' = [IdxDesc n c (show (pos q)) | IdxDesc n c q <- d]
+          pos n = i where Just i = elemIndex n nms
+          nms = map idxName 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
+
+wedge a b = antisym (prod (norper a) (norper b))
+
+a /\ b = wedge a b
+
+a <*> b = normal $ prod a b
+
+normAT t = sqrt $ innerAT t t
+
+innerAT t1 t2 = dot (ten t1) (ten t2) / fromIntegral (fact $ length $ dims t1)
+
+fact n = product [1..n]
+
+leviCivita n = antisym $ foldl1 prod $ zipWith tcovector (map show [1,2..]) (toRows (ident n))
+
+contractionF t1 t2 = contraction t1 n1 t2 n2
+    where n1 = fn t1
+          n2 = fn t2
+          fn = idxName . head . dims
+
+
+dualV vs = foldl' contractionF (leviCivita n) vs
+    where n = idxDim . head . dims . head $ vs
+
+raise (T d v) = T (map raise' d) v
+    where raise' idx@IdxDesc {idxType = Covariant } = idx {idxType = Contravariant}
+          raise' idx@IdxDesc {idxType = Contravariant } = idx {idxType = Covariant}
+
+dualMV t = prod (foldl' contract1b (lc <*> t) ds) (scalar (recip $ fromIntegral $ fact (length ds)))
+    where
+        lc = leviCivita n
+        nms1 = map idxName (dims lc)
+        nms2 = map ((++"'").idxName) (dims t)
+        ds = zip nms1 nms2
+        n = idxDim . head . dims $ t
+
+contract1b t (n1,n2) = contract1 t n1 n2
+
+contractions t pairs = foldl' contract1b t pairs
+
+asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
+
+partF t i = part t (name,i) where name = idxName . head . dims $ t
+
+niceAS t = filter ((/=0.0).fst) $ zip vals base
+    where vals = map ((`at` 0).ten.foldl' partF t) (map (map pred) base)
+          base = asBase r n
+          r = length (dims t)
+          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 (
     scale, addConstant, add, sub, mul,
 ) where
 
-import Data.Packed.Internal
+import Data.Packed.Internal.Common
+import Data.Packed.Internal.Vector
+
 import Complex
 import Foreign