From 834b4837799611fd7fbaa9609ea587e041cb0ca1 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 11 Sep 2007 10:25:24 +0000 Subject: tensor lib moved to easyVision project --- lib/Data/Packed/Internal/Matrix.hs | 5 +- lib/Data/Packed/Internal/Tensor.hs | 347 ------------------------------------- lib/Data/Packed/Tensor.hs | 31 ---- lib/GSL.hs | 4 +- 4 files changed, 5 insertions(+), 382 deletions(-) delete mode 100644 lib/Data/Packed/Internal/Tensor.hs delete mode 100644 lib/Data/Packed/Tensor.hs (limited to 'lib') diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index ba32a67..2db4838 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -9,9 +9,10 @@ -- Stability : provisional -- Portability : portable (uses FFI) -- --- Fundamental types +-- Internal matrix representation -- ----------------------------------------------------------------------------- +-- --#hide module Data.Packed.Internal.Matrix where @@ -74,7 +75,7 @@ type Mt t s = Int -> Int -> Ptr t -> s -- infixr 6 ::> -- type t ::> s = Mt t s --- | the inverse of 'fromLists' +-- | the inverse of 'Data.Packed.Matrix.fromLists' toLists :: (Field t) => Matrix t -> [[t]] toLists m = partit (cols m) . toList . cdat $ m diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs deleted file mode 100644 index dea1636..0000000 --- a/lib/Data/Packed/Internal/Tensor.hs +++ /dev/null @@ -1,347 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} - ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Internal.Tensor --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable (uses FFI) --- --- support for basic tensor operations --- ------------------------------------------------------------------------------ - -module Data.Packed.Internal.Tensor where - -import Data.Packed.Internal -import Foreign.Storable -import Data.List(sort,elemIndex,nub,foldl1',foldl') -import GSL.Vector -import Data.Packed.Matrix -import Data.Packed.Vector -import LinearAlgebra.Linear - -data IdxType = Covariant | Contravariant deriving (Show,Eq) - -type IdxName = String - -data IdxDesc = IdxDesc { idxDim :: Int, - idxType :: IdxType, - idxName :: IdxName } deriving Eq - -instance Show IdxDesc where - show (IdxDesc n t name) = name ++ sym t ++"["++show n++"]" - where sym Covariant = "_" - sym Contravariant = "^" - - -data Tensor t = T { dims :: [IdxDesc] - , ten :: Vector t - } - --- | returns the coordinates of a tensor in row - major order -coords :: Tensor t -> Vector t -coords = ten - -instance (Show a, Field a) => Show (Tensor a) where - show T {dims = [], ten = t} = "scalar "++show (t `at` 0) - show t = "("++shdims (dims t) ++") "++ showdt t - -asMatrix t = reshape (idxDim $ dims t!!1) (ten t) - -showdt t | rank t == 1 = show (toList (ten t)) - | rank t == 2 = ('\n':) . dsp . map (map show) . toLists $ asMatrix $ t - | otherwise = concatMap showdt $ parts t (head (names t)) - --- | a nice description of the tensor structure -shdims :: [IdxDesc] -> String -shdims [] = "" -shdims [d] = show d -shdims (d:ds) = show d ++ "><"++ shdims ds - --- | tensor rank (number of indices) -rank :: Tensor t -> Int -rank = length . dims - -names :: Tensor t -> [IdxName] -names t = map idxName (dims t) - --- | number of contravariant and covariant indices -structure :: Tensor t -> (Int,Int) -structure t = (rank t - n, n) where - n = length $ filter isCov (dims t) - isCov d = idxType d == Covariant - --- | creates a rank-zero tensor from a scalar -scalar :: Storable t => t -> Tensor t -scalar x = T [] (fromList [x]) - --- | Creates a tensor from a signed list of dimensions (positive = contravariant, negative = covariant) and a Vector containing the coordinates in row major order. -tensor :: [Int] -> Vector a -> Tensor a -tensor dssig vec = T d v `withIdx` seqind where - n = product (map abs dssig) - v = if dim vec == n then vec else error "wrong arguments for tensor" - d = map cr dssig - cr n | n > 0 = IdxDesc {idxName = "", idxDim = n, idxType = Contravariant} - | n < 0 = IdxDesc {idxName = "", idxDim = -n, idxType = Covariant } - - -tensorFromVector :: IdxType -> Vector t -> Tensor t -tensorFromVector tp v = T {dims = [IdxDesc (dim v) tp "1"], ten = v} - -tensorFromMatrix :: Field t => IdxType -> IdxType -> Matrix t -> Tensor t -tensorFromMatrix tpr tpc m = T {dims = [IdxDesc (rows m) tpr "1",IdxDesc (cols m) tpc "2"] - , ten = cdat m} - - -liftTensor :: (Vector a -> Vector b) -> Tensor a -> Tensor b -liftTensor f (T d v) = T d (f v) - -liftTensor2 :: (Vector a -> Vector b -> Vector c) -> Tensor a -> Tensor b -> Tensor c -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 - - - - --- | express the tensor as a matrix with the given index in columns -findIdx :: (Field t) => IdxName -> Tensor t - -> (([IdxDesc], [IdxDesc]), Matrix t) -findIdx name t = ((d1,d2),m) where - (d1,d2) = span (\d -> idxName d /= name) (dims t) - c = product (map idxDim d2) - m = matrixFromVector RowMajor c (ten t) - --- | express the tensor as a matrix with the given index in rows -putFirstIdx :: (Field t) => String -> Tensor t -> ([IdxDesc], Matrix t) -putFirstIdx name t = (nd,m') - where ((d1,d2),m) = findIdx name t - m' = matrixFromVector RowMajor c $ cdat $ trans m - nd = d2++d1 - c = dim (ten t) `div` (idxDim $ head d2) - - --- | renames all the indices in the current order (repeated indices may get different names) -withIdx :: Tensor t -> [IdxName] -> Tensor t -withIdx (T d v) l = T d' v - where d' = zipWith f d l - f i n = i {idxName=n} - - --- | raises or lowers all the indices of a tensor (with euclidean metric) -raise :: Tensor t -> Tensor t -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} - - --- | index transposition to a desired order. You can specify only a subset of the indices, which will be moved to the front of indices list -tridx :: (Field t) => [IdxName] -> Tensor t -> Tensor t -tridx [] t = t -tridx (name:rest) t = T (d:ds) (join ts) where - ((_,d:_),_) = findIdx name t - ps = map (tridx rest) (parts t name) - ts = map ten ps - ds = dims (head ps) - - - --- | extracts a given part of a tensor -part :: (Field t) => Tensor t -> (IdxName, Int) -> Tensor t -part t (name,k) = if k<0 || k>=l - then error $ "part "++show (name,k)++" out of range" -- in "++show t - else T {dims = ds, ten = toRows m !! k} - where (d:ds,m) = putFirstIdx name t - l = idxDim d - --- | creates a list with all parts of a tensor along a given index -parts :: (Field t) => Tensor t -> IdxName -> [Tensor t] -parts t name = map f (toRows m) - where (d:ds,m) = putFirstIdx name t - l = idxDim d - f t = T {dims=ds, ten=t} - - -compatIdx :: (Field t1, Field t) => Tensor t1 -> IdxName -> Tensor t -> IdxName -> Bool -compatIdx t1 n1 t2 n2 = compatIdxAux d1 d2 where - d1 = head $ snd $ fst $ findIdx n1 t1 - d2 = head $ snd $ fst $ findIdx n2 t2 - compatIdxAux IdxDesc {idxDim = n1, idxType = t1} - IdxDesc {idxDim = n2, idxType = t2} - = t1 /= t2 && n1 == n2 - - -outer' u v = dat (outer u v) - --- | tensor product without without any contractions -rawProduct :: (Field t, Num t) => Tensor t -> Tensor t -> Tensor t -rawProduct (T d1 v1) (T d2 v2) = T (d1++d2) (outer' v1 v2) - --- | contraction of the product of two tensors -contraction2 :: (Field t, Num t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t -contraction2 t1 n1 t2 n2 = - if compatIdx t1 n1 t2 n2 - then T (tail d1 ++ tail d2) (cdat m) - else error "wrong contraction2" - where (d1,m1) = putFirstIdx n1 t1 - (d2,m2) = putFirstIdx n2 t2 - m = multiply RowMajor (trans m1) m2 - --- | contraction of a tensor along two given indices -contraction1 :: (Linear Vector t) => Tensor t -> IdxName -> IdxName -> Tensor t -contraction1 t name1 name2 = - if compatIdx t name1 t name2 - then sumT y - else error $ "wrong contraction1: "++(shdims$dims$t)++" "++name1++" "++name2 - where d = dims (head y) - x = (map (flip parts name2) (parts t name1)) - y = map head $ zipWith drop [0..] x - --- | contraction of a tensor along a repeated index -contraction1c :: (Linear Vector t) => Tensor t -> IdxName -> Tensor t -contraction1c t n = contraction1 renamed n' n - where n' = n++"'" -- hmmm - renamed = withIdx t auxnames - auxnames = h ++ (n':r) - (h,_:r) = break (==n) (map idxName (dims t)) - --- | alternative and inefficient version of contraction2 -contraction2' :: (Linear Vector t) => Tensor t -> IdxName -> Tensor t -> IdxName -> Tensor t -contraction2' t1 n1 t2 n2 = - if compatIdx t1 n1 t2 n2 - then contraction1 (rawProduct t1 t2) n1 n2 - else error "wrong contraction'" - --- | applies a sequence of contractions -contractions :: (Linear Vector t) => Tensor t -> [(IdxName, IdxName)] -> Tensor t -contractions t pairs = foldl' contract1b t pairs - where contract1b t (n1,n2) = contraction1 t n1 n2 - --- | applies a sequence of contractions of same index -contractionsC :: (Linear Vector t) => Tensor t -> [IdxName] -> Tensor t -contractionsC t is = foldl' contraction1c t is - - --- | applies a contraction on the first indices of the tensors -contractionF :: (Linear Vector t) => Tensor t -> Tensor t -> Tensor t -contractionF t1 t2 = contraction2 t1 n1 t2 n2 - where n1 = fn t1 - n2 = fn t2 - fn = idxName . head . dims - --- | computes all compatible contractions of the product of two tensors that would arise if the index names were equal -possibleContractions :: (Linear Vector t) => Tensor t -> Tensor t -> [Tensor t] -possibleContractions t1 t2 = [ contraction2 t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ] - - - -desiredContractions2 :: Tensor t -> Tensor t1 -> [(IdxName, IdxName)] -desiredContractions2 t1 t2 = [ (n1,n2) | n1 <- names t1, n2 <- names t2, n1==n2] - -desiredContractions1 :: Tensor t -> [IdxName] -desiredContractions1 t = [ n1 | (a,n1) <- x , (b,n2) <- x, a/=b, n1==n2] - where x = zip [0..] (names t) - --- | tensor product with the convention that repeated indices are contracted. -mulT :: (Linear Vector t) => Tensor t -> Tensor t -> Tensor t -mulT t1 t2 = r where - t1r = contractionsC t1 (desiredContractions1 t1) - t2r = contractionsC t2 (desiredContractions1 t2) - cs = desiredContractions2 t1r t2r - r = case cs of - [] -> rawProduct t1r t2r - (n1,n2):as -> contractionsC (contraction2 t1r n1 t2r n2) (map fst as) - ------------------------------------------------------------------ - --- | tensor addition (for tensors with the same structure) -addT :: (Linear Vector a) => Tensor a -> Tensor a -> Tensor a -addT a b = liftTensor2 add a b - -sumT :: (Linear Vector a) => [Tensor a] -> Tensor a -sumT l = foldl1' addT l - ------------------------------------------------------------------ - --- sent to Haskell-Cafe by Sebastian Sylvan -perms :: [t] -> [[t]] -perms [x] = [[x]] -perms xs = [y:ps | (y,ys) <- selections xs , ps <- perms ys] - where - selections [] = [] - selections (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- selections xs] - -interchanges :: (Ord a) => [a] -> Int -interchanges ls = sum (map (count ls) ls) - where count l p = length $ filter (>p) $ take pel l - where Just pel = elemIndex p l - -signature :: (Num t, Ord a) => [a] -> t -signature l | length (nub l) < length l = 0 - | even (interchanges l) = 1 - | otherwise = -1 - - -sym :: (Linear Vector t) => Tensor t -> Tensor t -sym t = T (dims t) (ten (sym' (withIdx t seqind))) - where sym' t = sumT $ map (flip tridx t) (perms (names t)) - where nms = map idxName . dims - -antisym :: (Linear Vector t) => Tensor t -> Tensor t -antisym t = T (dims t) (ten (antisym' (withIdx t seqind))) - where antisym' t = sumT $ map (scsig . flip tridx t) (perms (names t)) - scsig t = scalar (signature (nms t)) `rawProduct` t - where nms = map idxName . dims - --- | the wedge product of two tensors (implemented as the antisymmetrization of the ordinary tensor product). -wedge :: (Linear Vector t, Fractional t) => Tensor t -> Tensor t -> Tensor t -wedge a b = antisym (rawProduct (norper a) (norper b)) - where norper t = rawProduct t (scalar (recip $ fromIntegral $ fact (rank t))) - --- antinorper t = rawProduct t (scalar (fromIntegral $ fact (rank t))) - --- | The euclidean inner product of two completely antisymmetric tensors -innerAT :: (Fractional t, Field t) => Tensor t -> Tensor t -> t -innerAT t1 t2 = dot (ten t1) (ten t2) / fromIntegral (fact $ rank t1) - -fact :: (Num t, Enum t) => t -> t -fact n = product [1..n] - -seqind' :: [[String]] -seqind' = map return seqind - -seqind :: [String] -seqind = map show [1..] - --- | completely antisymmetric covariant tensor of dimension n -leviCivita :: (Linear Vector t) => Int -> Tensor t -leviCivita n = antisym $ foldl1 rawProduct $ zipWith withIdx auxbase seqind' - where auxbase = map tc (toRows (ident n)) - tc = tensorFromVector Covariant - --- | contraction of leviCivita with a list of vectors (and raise with euclidean metric) -innerLevi :: (Linear Vector t) => [Tensor t] -> Tensor t -innerLevi vs = raise $ foldl' contractionF (leviCivita n) vs - where n = idxDim . head . dims . head $ vs - - --- | obtains the dual of a multivector (with euclidean metric) -dual :: (Linear Vector t, Fractional t) => Tensor t -> Tensor t -dual t = raise $ leviCivita n `mulT` withIdx t seqind `rawProduct` x - where n = idxDim . head . dims $ t - x = scalar (recip $ fromIntegral $ fact (rank t)) - - --- | shows only the relevant components of an antisymmetric tensor -niceAS :: (Field t, Fractional t) => Tensor t -> [(t, [Int])] -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 - partF t i = part t (name,i) where name = idxName . head . dims $ t - asBase r n = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n] diff --git a/lib/Data/Packed/Tensor.hs b/lib/Data/Packed/Tensor.hs deleted file mode 100644 index 27747ac..0000000 --- a/lib/Data/Packed/Tensor.hs +++ /dev/null @@ -1,31 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Data.Packed.Tensor --- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : experimental --- Portability : portable --- --- Basic tensor operations --- ------------------------------------------------------------------------------ - -module Data.Packed.Tensor ( - -- * Construction - Tensor, tensor, scalar, - -- * Manipulation - IdxName, IdxType(..), IdxDesc(..), structure, dims, coords, parts, - tridx, withIdx, raise, - -- * Operations - addT, mulT, - -- * Exterior Algebra - wedge, dual, leviCivita, innerLevi, innerAT, niceAS, - -- * Misc - liftTensor, liftTensor2 -) where - - -import Data.Packed.Internal.Tensor - diff --git a/lib/GSL.hs b/lib/GSL.hs index e6ff6df..a2137d3 100644 --- a/lib/GSL.hs +++ b/lib/GSL.hs @@ -16,7 +16,7 @@ module GSL ( module Data.Packed.Vector, module Data.Packed.Matrix, -module Data.Packed.Tensor, +--module Data.Packed.Tensor, module LinearAlgebra.Algorithms, module LAPACK, module GSL.Integration, @@ -34,7 +34,7 @@ module Complex import Data.Packed.Vector hiding (constant) import Data.Packed.Matrix hiding ((><), multiply) -import Data.Packed.Tensor +--import Data.Packed.Tensor import LinearAlgebra.Algorithms hiding (pnorm) import LAPACK import GSL.Integration -- cgit v1.2.3