packages feed

hTensor (empty) → 0.1.0

raw patch · 15 files changed

+1439/−0 lines, 15 filesdep +basedep +containersdep +haskell98setup-changed

Dependencies added: base, containers, haskell98, hmatrix

Files

+ INSTALL view
@@ -0,0 +1,3 @@+INSTALLATION++$ cabal install hTensor
+ LICENSE view
@@ -0,0 +1,2 @@+Copyright Alberto Ruiz 2009+GPL license
+ README view
+ Setup.lhs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain
+ examples/array.hs view
@@ -0,0 +1,59 @@+import Numeric.LinearAlgebra.Array+import Numeric.LinearAlgebra.Array.Util+import Control.Applicative++-- 'listArray' specialized for Array Double+infixl 9 #+(#) :: [Int] -> [Double] -> Array Double+(#) = listArray++(<|) :: Name -> [Array Double] -> Array Double+infixl 8 <|+n <| ls = index n ls++i = ("i" <|)+j = ("j" <|)+k = ("k" <|)++sh x = putStrLn . formatFixed 2 $ x++a = [3,4,2] # [1..]+b = [2,3] # [5,6]+c = 7 :: Array Double+s = [2,2,2,2] # [1..]++t = [3,3,3,3]#[1 ..]!"ijkl"++q = [2,4,3] # (fun <$> r 2 <*> r 4 <*> r 3) !"ijk"+    where r k = [1..k]+          fun = \i j k -> i*2*j-k++mk ds f = ds # map f (sequence $ map (enumFromTo 1 . fromIntegral) $ ds)++m = j [i[2,0,0], i[1,0,1], i[0,3,0]] ~> "ij"++main = do+    putStrLn "8-dimensional array"+    sh $ (replicate 8 2)#[1::Double ..]!(take 8 ['a'..])+    ------------------------+    putStrLn "different display formats"+    sh $ a!"ijk"+    printA "%7.3f" a+    putStrLn . formatScaled 2 $ a!"ijk"+    sh . noIdx $ a+    ------------------------+    putStrLn "array defined using a function"+    sh q+    sh $ mk [2,4,3] (\[i,j,k] -> i*2*j-k) !"ijk"+    ------------------------+    putStrLn "contraction"+    sh t+    sh $ t!"ijkk"+    ------------------------+    putStrLn "tensor product"+    sh $ m+    sh $ (t !"pqrs" * m!"kr") ~> "pqks"+    ------------------------+    putStrLn "automatic conformability"+    sh $ j[1,2,3] + k[10,20]+    sh $ k [m, 3*m-1, 7, i[1,2,3]]
+ examples/exterior.hs view
@@ -0,0 +1,41 @@+import Numeric.LinearAlgebra.Exterior+import Numeric.LinearAlgebra.Array.Util(formatFixed,asMatrix)+import Numeric.LinearAlgebra (det,(><))++printAS = print . asMultivector++sh = putStrLn . formatFixed 2++-- 'listTensor' specialized for Tensor Double+infixl 9 #+(#) :: [Int] -> [Double] -> Tensor Double+(#) = listTensor++m = [3,-3]#[ 1,2,5,+             1,2,8,+            -2,0,4]++eps = leviCivita 3++a = vector [1,0,0] /\ vector [0,1,0]+b = vector [2,0,100,0] /\ vector [0,3,0,0] /\ vector [0,0,4,0]++im = eps!"ijb"* m!"pi" * m!"qj" * cov eps!"apq"++main = do+    putStrLn "exterior product"+    print a+    sh a+    printAS a+    printAS b+    putStrLn "\ndeterminant"+    printAS $ cov eps!"pqr" * m!"pi" * m!"qj" * m!"rk"+    print $ det (asMatrix m)+    putStrLn "\ninverse"+    sh $ im+    sh $ im!"ik" * m!"kj"+    putStrLn "\nmeet and join"+    printAS $ (vector [1,0,1] /\ vector [0,1,0]) \/ (vector [1,1,0] /\ vector [0,0,1])+    putStrLn "\nEuclidean inner product of r-vectors"+    printAS $ (vector [1,0,1] /\ vector [0,1,0]) `inner` (vector[1,1,1])+    print $ vector[3,5] `inner` vector[2,1]
+ examples/geom.hs view
@@ -0,0 +1,29 @@+import Numeric.LinearAlgebra.Multivector++o  = e 4+x = vector[1,0,0,1]+y = o + e 2+z = vector[0,0,1] + e 4++p = x /\ y /\ z++p1 = o /\ x /\ y+p2 = z /\ vector[1,0,1,1] /\ vector[0,1,1,1]++l = o /\ vector[1,1,1,1]++rot = rotor 3 (pi/4) (e 3)++l' = rot * l * rever rot++inh v = v / (v -| e 4) - e 4++main = do+     print l+     print $ l \/ p1+     print $ l \/ p2+     print $ l \/ p+     print $ inh $ l \/ p1+     print $ inh $ l \/ p2+     print $ inh $ l' \/ p1+     print $ inh $ l' \/ p2
+ hTensor.cabal view
@@ -0,0 +1,65 @@+Name:               hTensor+Version:            0.1.0+License:            GPL+License-file:       LICENSE+Author:             Alberto Ruiz+Maintainer:         Alberto Ruiz <aruiz@um.es>+Stability:          experimental+Homepage:           http://perception.inf.um.es/tensor+Synopsis:           Multidimensional arrays and simple tensor computations.+Description:+ This is an experimental library for multidimensional arrays,+ oriented to support simple tensor computations and multilinear+ algebra.+ .+ Array dimensions have an \"identity\" which is preserved+ in data manipulation. Indices are explicitly selected by name in+ expressions, and Einstein's summation convention for repeated indices+ is automatically applied.+ .+ The library has a purely functional interface: arrays are immutable,+ and operations typically work on whole structures which can be assembled+ and decomposed using simple primitives. Arguments are automatically made conformant+ by replicating them along extra dimensions appearing in an operation.+ There is preliminary support for Geometric Algebra.+ .+ A tutorial can be found in the website of the project.++Category:           Math+tested-with:        GHC ==6.10.3++cabal-version:      >=1.2+build-type:         Simple++extra-source-files: README INSTALL++extra-source-files: examples/array.hs+                    examples/exterior.hs+                    examples/geom.hs++flag splitBase+    description:    Choose the new smaller, split-up base package.++library+    if flag(splitBase)+      build-depends:    base >= 3 && < 5+    else+      build-depends:    base < 3++    Build-Depends:      haskell98, hmatrix >= 0.5, containers++    hs-source-dirs:     lib+    Exposed-modules:    Numeric.LinearAlgebra.Array.Simple+                        Numeric.LinearAlgebra.Array.Util+                        Numeric.LinearAlgebra.Array+                        Numeric.LinearAlgebra.Tensor+                        Numeric.LinearAlgebra.Exterior+                        Numeric.LinearAlgebra.Multivector++    other-modules:      Numeric.LinearAlgebra.Array.Internal++    ghc-prof-options:   -auto-all++    ghc-options:        -Wall+                        -fno-warn-missing-signatures+                        -fno-warn-orphans
+ lib/Numeric/LinearAlgebra/Array.hs view
@@ -0,0 +1,77 @@+{-# LANGUAGE UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.LinearAlgebra.Array+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Simple multidimensional array with useful numeric instances.+--+-- Contractions only require equal dimension.+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array (+    None(..),+    Array,+    listArray,+    scalar,+    index,+    (!),(~>),+    (.*),+    printA+) where++import Numeric.LinearAlgebra.Array.Simple+import Numeric.LinearAlgebra.Array.Util+import Numeric.LinearAlgebra.Array.Internal(printA)+import Data.Packed(Vector)++-- | Create an 'Array' from a list of parts (@index = 'newIndex' 'None'@).+index :: Coord t => Name -> [Array t] -> Array t+index = newIndex None+++-- | Element by element product.+infixl 7 .*+(.*) :: (Coord a, Compat i) => NArray i a -> NArray i a -> NArray i a+(.*) = zipArray (*)++instance (Coord t, Compat i) => Eq (NArray i t) where+    t1 == t2 = sameStructure t1 t2 && coords t1 == coords (reorder (names t1) t2)++instance (Show (NArray i t), Coord t, Compat i) => Num (NArray i t) where+    (+) = zipArray (+)+    (*) = (|*|)+    negate t = scalar (-1) * t+    fromInteger n = scalar (fromInteger n)+    abs _ = error "abs for arrays not defined"+    signum _ = error "signum for arrays not defined"++instance (Coord t, Compat i, Num (NArray i t)) => Fractional (NArray i t) where+    fromRational = scalar . fromRational+    (/) = zipArray (/)+    recip = mapArray recip++instance (Coord t, Compat i, Fractional (NArray i t), Floating (Vector t)) => Floating (NArray i t) where+    sin   = mapArray sin+    cos   = mapArray cos+    tan   = mapArray tan+    asin  = mapArray asin+    acos  = mapArray acos+    atan  = mapArray atan+    sinh  = mapArray sinh+    cosh  = mapArray cosh+    tanh  = mapArray tanh+    asinh = mapArray asinh+    acosh = mapArray acosh+    atanh = mapArray atanh+    exp   = mapArray exp+    log   = mapArray log+    (**)  = zipArray (**)+    sqrt  = mapArray sqrt+    pi    = scalar pi
+ lib/Numeric/LinearAlgebra/Array/Internal.hs view
@@ -0,0 +1,565 @@+{-# OPTIONS_HADDOCK hide #-}+{-# LANGUAGE FlexibleInstances, FlexibleContexts, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Internal+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Multidimensional arrays.+--+-- The arrays provided by this library are immutable, built on top of hmatrix+-- structures.+-- Operations work on complete structures (indexless), and dimensions have \"names\", +-- in order to select the desired contractions in tensor computations.+--+-- This module contains auxiliary functions not required by the end user.++-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Internal (+    -- * Data structures+    NArray, Idx(..), Name,+    rank, names, size, typeOf , dims, coords,+    Compat(..),+    -- * Array creation+    scalar,+    mkNArray,+    fromVector, fromMatrix,+    -- * Array manipulation+    rename,(!),+    parts,+    (|*|),+    zipArray,+    mapArray,+    extract,+    onIndex,+    -- * Utilities+    reorder, (~>),+    sameStructure,+    conformable,+    makeConformant,+    mapTypes,+    renameRaw,+    formatArray, formatFixed, formatScaled, printA,+    showBases,+    newIndex,+    dummyAt, noIdx,+    basisOf,+    common,+    Coord,+    asMatrix, asVector, asScalar+) where++import Data.Packed+import Data.List+import Numeric.LinearAlgebra(outer,multiply,Field)+import Control.Applicative+import Data.Function(on)+import Text.Printf++-- | Types that can be elements of the multidimensional arrays.+class (Num (Vector t), Field t) => Coord t+instance Coord Double+instance Coord (Complex Double)++-- import Debug.Trace+-- +-- debug s f x = trace (s ++ ": " ++ show (f x)) x++-- | indices are denoted by strings, (frequently single-letter)+type Name = String++-- | Dimension descriptor.+data Idx i = Idx { iDim  :: Int+                 , iName :: Name+                 , iType :: i+                 } deriving (Eq)++++-- | A multidimensional array with index type i and elements t.+data NArray i t = A { dims   :: [Idx i]   -- ^ Get detailed dimension information about the array.+                    , coords :: Vector t  -- ^ Get the coordinates of an array as a+                                          -- flattened structure (in the order specified by 'dims').+                    }+++-- | development function not intended for the end user+mkNArray :: [Idx i] -> Vector a -> NArray i a+mkNArray [] _ = error "array with empty dimensions, use scalar"+mkNArray dms vec = A dms v where+    ds = map iDim dms+    n = product ds+    v = if dim vec == n && minimum ds > 0+            then vec+            else error $ show ds ++ " dimensions and " +++                         show (dim vec) ++ " coordinates for mkNArray"++-- | Create a 0-dimensional structure.+scalar :: Coord t => t -> NArray i t+scalar x = A [] (fromList [x])+++-- | 'rename' the indices with single-letter names. Equal indices of compatible type are contracted out.+infixl 8 !+(!) :: (Coord t, Compat i)+       => NArray i t+       -> String   -- ^ new indices+       -> NArray i t+t ! ns = rename t (map return ns)++-- | Rename indices. Equal indices are contracted out.+rename :: (Coord t, Compat i)+       => NArray i t+       -> [Name]     -- ^ new names+       -> NArray i t+rename t ns = reorder orig (contract t')+    where t' = renameRaw t ns+          orig = nub (names t') \\ common1 t'+++renameRaw (A d v) l | length l == length d = A d' v+                    | otherwise = error $ "rename " ++ show d ++ " with " ++ show l+    where d' = zipWith f d l+          f i n = i {iName=n}++mapDims f (A d v) = A (map f d) v++mapTypes :: (i1 -> i) -> NArray i1 t -> NArray i t+mapTypes f = mapDims (\i -> i {iType = f (iType i)})++-- mapNames f = mapDims (\i -> i {iName = f (iName i)})++-- | Index names.+names :: NArray i t -> [Name]+names = map iName . dims++-- | Dimension of given index.+size :: Name -> NArray i t -> Int+size n t = (iDim . head) (filter ((n==).iName) (dims t))++-- | Type of given index.+typeOf :: Compat i => Name -> NArray i t -> i+typeOf n t = (iType . head) (filter ((n==).iName) (dims t))+++-- | The number of dimensions of a multidimensional array.+rank :: NArray i t -> Int+rank = length . dims++----------------------------------------------------------++lastIdx name t = ((d1,d2),m) where+    (d1,d2) = span (\d -> iName d /= name) (dims t)+    c = product (map iDim d2)+    m = reshape c (coords t)++firstIdx name t = (nd,m')+    where ((d1,d2),m) = lastIdx name t+          m' = reshape c $ flatten $ trans m+          nd = d2++d1+          c = dim (coords t) `div` (iDim $ head d2)+++-- | Create a list of the substructures at the given level.+parts :: (Coord t) +      => NArray i t+      -> Name        -- ^ index to expand+      -> [NArray i t]+parts a name | name `elem` (names a) = map (reorder orig) (partsRaw a name)+             | otherwise = error $ "parts: " ++ show name ++ " is not a dimension of "++(show $ names a)+    where orig = names a \\ [name]++partsRaw a name = map f (toRows m)+    where (_:ds,m) = firstIdx name a+          f t = A {dims=ds, coords=t}++tridx [] t = t+tridx (name:rest) t = A (d:ds) (join ts) where+    d = case lastIdx name t of+            ((_,d':_),_) -> d'+            _ -> error "wrong index sequence to reorder"+    ps = map (tridx rest) (partsRaw t name)+    ts = map coords ps+    ds = dims (head ps)++-- | Change the internal layout of coordinates.+-- The array, considered as an abstract object, does not change.+reorder :: (Coord t) => [Name] -> NArray i t -> NArray i t+reorder ns b | ns == names b = b+             | sort ns == sort (names b) = tridx ns b+             | otherwise = error $ "wrong index sequence " ++ show ns+                                    ++ " to reorder "++(show $ names b)+++-- | 'reorder' (transpose) the dimensions of the array (with single letter names).+--+-- Operations are defined by named indices, so the transposed array is operationally equivalent to the original one.+infixl 8 ~>+(~>) :: (Coord t) => NArray i t -> String -> NArray i t+t ~> ns = reorder (map return ns) t++-------------------------------------------------------------++rawProduct (A d1 v1) (A d2 v2) = A (d1++d2) (flatten (outer v1 v2))++----------------------------------------------------------------------++-- | Apply a function (defined on hmatrix 'Vector's) to all elements of a structure.+-- Use @mapArray (mapVector f)@ for general functions.+mapArray :: Coord b => (Vector a -> Vector b) -> NArray i a -> NArray i b+mapArray f t+    | null (dims t) = scalar (f (coords t)@>0)+    | otherwise = mkNArray (dims t) (f (coords t))++liftNA2 f (A d1 v1) (A _d2 v2) = A d1 (f v1 v2)++-- | Class of compatible indices for contractions.+class (Eq a, Show (Idx a)) => Compat a where+    compat :: Idx a -> Idx a -> Bool++++contract1 t name1 name2 | ok = foldl1' (liftNA2 (+)) y+                        | otherwise = error $ "wrong contraction1: "+                                    ++(show $ dims t)++" "+                                    ++ name1++" "++name2+    where ok = (compat <$> getName t name1 <*> getName t name2) == Just True+          x = map (flip partsRaw name2) (partsRaw t name1)+          y = map head $ zipWith drop [0..] x++getName t name = d where+    l = filter ((==name).iName) (dims t)+    d = if null l+            then Nothing+            else Just (head l)++contract1c t n = contract1 renamed n n'+    where n' = " "++n++" " -- forbid spaces in names...+          renamed = renameRaw (t) auxnames+          auxnames = h ++ (n':r)+          (h,_:r) = break (==n) (names t)++common1 t = [ n1 | (a,n1) <- x , (b,n2) <- x, a>b, n1==n2]+    where x = zip [0 ::Int ..] (names t)++contract t = foldl' contract1c t (common1 t)++----------------------------------------------------------------------++contract2 t1 t2 n | ok = A (tail ds1 ++ tail ds2) (flatten m)+                  | otherwise = error $ "wrong contraction2: "++ n ++ " of "+++                                      (show $ dims t1)++" and "++ (show $ dims t2)+  where ok = (compat <$> getName t1 n <*> getName t2 n) == Just True+        (ds1,m1) = firstIdx n t1+        (ds2,m2) = firstIdx n t2+        m = (trans m1) `multiply` m2++common2 t1 t2 = [ n1 | n1 <- names t1, n2 <- names t2, n1==n2]++infixl 5 |*|+-- | Tensor product with automatic contraction of repeated indices, following Einstein summation convention.+(|*|) :: (Coord t, Compat i)+      => NArray i t -> NArray i t -> NArray i t+t1 |*| t2 = r where+    cs = common2 t1 t2+    r = case cs of+        [] -> rawProduct t1 t2+        n:_ -> reorder orig $ contract (contract2 t1 t2 n)+    orig = nub (names t1 ++ names t2) \\ cs++-------------------------------------------------------------++-- | Check if two arrays have the same structure.+sameStructure :: (Eq i) => NArray i t1 -> NArray i t2 -> Bool+sameStructure a b = sortBy (compare `on` iName) (dims a) == sortBy (compare `on` iName) (dims b)++-------------------------------------------------------------++-- | Apply a function on vectors to conformant arrays. Two arrays are 'conformant' if+-- the dimensional structure of one of them is contained in the other one. The smaller+-- structure is replicated along the extra dimensions. The result has the same index+-- order as the largest structure (or as the first argument, if they are equal).+zipArray :: (Coord a, Coord b, Compat i)+   => (Vector a -> Vector b -> Vector c) -- ^ transformation+   -> NArray i a+   -> NArray i b+   -> NArray i c+zipArray o a b = liftNA2 o a' b' where+    (a',b') = makeConformantT (a,b)++-------------------------------------------------------++showBases x = f $ concatMap (shbld) x+    where   shbld (c,[]) = shsign c ++ showc c+            shbld (c,l) = shsign c ++ g (showc c) ++ "{"++ concatMap show l++"}"+            shsign c = if c < 0 then " - " else " + "+            showc c+                | abs (fromIntegral (round c :: Int) - c) <1E-10  = show (round $ abs c::Int)+                | otherwise = printf "%.3f" (abs c)+            f (' ':'+':' ':rs) = rs+            f (' ':'-':' ':rs) = '-':rs+            f a = a+            g "1" = ""+            g a = a++---------------------------------------------------------++data Rect = Rect { li :: Int, co :: Int, els :: [String] }++rect s = pad r c (Rect r 0 ss)+    where ss = lines s+          r  = length ss+          c  = maximum (map length ss)++pad nr nc (Rect r c ss) = Rect (r+r') (c+c') ss'' where+    r' = max 0 (nr-r)+    c' = max 0 (nc-c)+    ss' = map (padH nc) ss+    ss'' = replicate r' (replicate nc '-') ++ ss'+    padH l s = take (l-length s) (" | "++repeat ' ') ++ s++dispH :: Int -> [Rect] -> Rect+dispH k rs = Rect nr nc nss where+    nr = maximum (map li rs)+    nss' = mapTail (\x-> pad nr (co x + k) x) rs+    nss = foldl1' (zipWith (++)) (map els nss')+    nc = length (head nss)++dispV :: Int -> [Rect] -> Rect+dispV k rs = Rect nr nc nss where+    nc = maximum (map co rs)+    nss' = mapTail (\x-> pad (li x + k) nc x) rs+    nss = concatMap els nss'+    nr = length nss++mapTail f (a:b) = a : map f b+mapTail _ x     = x+++++formatAux f x = unlines . addds . els . fmt ms $ x where+    fmt [] _ = undefined -- cannot happen+    fmt (g:gs) t+        | rank t == 0 = rect (f (coords t @> 0))+        | rank t == 1 =  rect $ unwords $ map f (toList $ coords t)+        | rank t == 2 =  decor t $ rect $ w1 $ format " " f (reshape (iDim $ last $ dims t) (coords t))+        | otherwise    = decor t (g ps)+      where ps = map (fmt gs ) (partsRaw t (head (names t)))+    ds = showNice (filter ((/='*').head.iName) $ dims x)+    addds = if null ds then (showRawDims (dims x) :) else (ds:)+    w1 = unlines . map (' ':) . lines+    ms = cycle [dispV 1, dispH 2]+    decor t | odd (rank t) = id+            | otherwise = decorLeft  (names t!!0) . decorUp (names t!!1)+++showNice x = unwords . intersperse "x" . map show $ x+showRawDims = showNice . map iDim . filter ((/="*").iName)++------------------------------------------------------++-- | Show a multidimensional array as a nested 2D table.+formatArray :: (Coord t, Compat i)+      => (t -> String) -- ^ format function (eg. printf \"5.2f\")+      -> NArray i t+      -> String+formatArray f t | odd (rank t) = formatAux f (dummyAt 0 t)+            | otherwise    = formatAux f t+++decorUp s rec+    | head s == '*' = rec+    | otherwise     = dispV 0 [rs,rec]+  where+    c = co rec+    c1 = (c - length s) `div` 2+    c2 = c - length s - c1+    rs = rect $ replicate c1 ' ' ++ s ++ replicate c2 ' '++decorLeft s rec+    | head s == '*' = rec+    | otherwise     = dispH 0 [rs,rec]+  where+    c = li rec+    r1 = (c - length s+1) `div` 2+    r2 = c - length s - r1+    rs = rect $ unlines $ replicate r1 spc ++ s : replicate (r2) spc+    spc = replicate (length s) ' '++------------------------------------------------------++-- | Print the array as a nested table with the desired format (e.g. %7.2f) (see also 'formatArray', and 'formatScaled').+printA :: (Coord t, Compat i, PrintfArg t) => String -> NArray i t -> IO ()+printA f t = putStrLn (formatArray (printf f) t)+++-- | Show the array as a nested table with autoscaled entries.+formatScaled :: (Compat i)+      => Int -- ^ number of of decimal places+      -> NArray i Double+      -> String+formatScaled dec t = unlines (('(':d++")  E"++show o) : m)+    where ss = formatArray (printf fmt. g) t+          d:m = lines ss+          g x = x/10^(o::Int)+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ coords t+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"++-- | Show the array as a nested table with a \"\%.nf\" format. If all entries+-- are approximate integers the array is shown without the .00.. digits.+formatFixed :: (Compat i)+      => Int -- ^ number of of decimal places+      -> NArray i Double+      -> String+formatFixed dec t+    | isInt t   = formatArray (printf ('%': show (width t) ++".0f")) t+    | otherwise = formatArray (printf ('%': show (width t+dec+1) ++"."++show dec ++"f")) t++isInt = all lookslikeInt . toList . coords+lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx+    where shx = show x+-- needsSign t = vectorMin (coords t) < 0+-- width :: Compat i => NArray i Double -> Int+width = maximum . map (length . (printf "%.0f"::Double->String)) . toList . coords+-- width t = k + floor (logBase 10 (max 1 $ vectorMax (abs $ coords t))) :: Int+--      where k | needsSign t = 2+--              | otherwise = 1++------------------------------------------------------++-- | Create an array from a list of subarrays. (The inverse of 'parts'.)+newIndex:: (Coord t, Compat i) =>+     i  -- ^ index type+     -> Name+     -> [NArray i t]+     -> NArray i t+newIndex i name ts = r where+    ds = Idx (length ts) name i : (dims (head cts))+    cts = makeConformant ts+    r = mkNArray ds (join $ map coords cts)+++-- | Insert a dummy index of dimension 1 at a given level (for formatting purposes).+dummyAt :: Int -> NArray i t -> NArray i t+dummyAt k t = mkNArray d' (coords t) where+    (d1,d2) = splitAt k (dims t)+    d' = d1 ++ d : d2+    d = Idx 1 "*" undefined++-- | Rename indices so that they are not shown in formatted output.+noIdx :: Compat i => NArray i t -> NArray i t+noIdx t = renameRaw t (map ('*':) (names t))++-- | Obtain a canonical base for the array.+basisOf :: Coord t => NArray i t -> [NArray i t]+basisOf t = map (dims t `mkNArray`) $ toRows (ident . dim . coords $ t)++-------------------------------------------------------------++instance (Coord t, Coord (Complex t), Compat i, Container Vector t) => Container (NArray i) t where+    toComplex (r,c) = zipArray (curry toComplex) r c+    fromComplex t = let (r,c) = fromComplex (coords t)+                     in (mapArray (const r) t, mapArray (const c) t)+    comp = mapArray comp+    conj = mapArray conj+    real = mapArray real+    complex = mapArray complex++----------------------------------------------------------------------++-- | obtains the common value of a property of a list+common :: (Eq a) => (b->a) -> [b] -> Maybe a+common f = commonval . map f where+    commonval :: (Eq a) => [a] -> Maybe a+    commonval [] = Nothing+    commonval [a] = Just a+    commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing++------------------------------------------------------------------------++-- | Extract the 'Matrix' corresponding to a two-dimensional array,+-- in the rows,cols order.+asMatrix :: (Coord t) => NArray i t -> Matrix t+asMatrix a | rank a == 2 = reshape c (coords a)+           | otherwise = error $ "asMatrix requires a rank 2 array."+    where c = size (last (names a)) a++-- | Extract the 'Vector' corresponding to a one-dimensional array.+asVector :: (Coord t) => NArray i t -> Vector t+asVector a | rank a == 1 = coords a+           | otherwise = error $ "asVector requires a rank 1 array."++-- | Extract the scalar element corresponding to a 0-dimensional array.+asScalar :: (Coord t) => NArray i t -> t+asScalar a | rank a == 0 = coords a @>0+           | otherwise = error $ "asScalar requires a rank 0 array."++------------------------------------------------------------------------++-- | Create a rank-1 array from an hmatrix 'Vector'.+fromVector :: Compat i => i -> Vector t -> NArray i t+fromVector i v = mkNArray [Idx (dim v) "1" i ] v++-- | Create a rank-2 array from an hmatrix 'Matrix'.+fromMatrix :: (Compat i, Coord t) => i -> i -> Matrix t -> NArray i t+fromMatrix ir ic m = mkNArray [Idx (rows m) "1" ir,+                               Idx (cols m) "2" ic] (flatten m)++------------------------------------------------------------------------++-- | Select some parts of a tensor, taking into account position and value.+extract :: (Compat i, Coord t)+        => (Int -> NArray i t -> Bool)+        -> Name+        -> NArray i t+        -> NArray i t+extract f name arr = reorder (names arr)+                   . newIndex (typeOf name arr) name+                   . map snd . filter (uncurry f)+                   $ zip [1..] (parts arr name)++-- | Apply a list function to the parts of an array at a given index.+onIndex :: (Coord a, Coord b, Compat i) =>+     ([NArray i a] -> [NArray i b])+     -> Name+     -> NArray i a+     -> NArray i b+onIndex f name t = reorder (names t) $ newIndex (typeOf name t) name (f (parts t name))++------------------------------------------------------------------------++extend alldims (A d v) = reorder (allnames) s where+    allnames = map iName alldims+    pref = alldims \\ d+    n = product (map iDim pref)+    s = A (pref++d) (join (replicate n v))++-- | Obtains most general structure of a list of dimension specifications+conformable :: Compat i => [[Idx i]] -> Maybe [Idx i]+conformable ds | ok        = Just alldims+               | otherwise = Nothing+    where alldims = nub (concat ds)+          allnames = map iName alldims+          ok = length (allnames) == length (nub allnames)++-- | Converts a list of arrays to a common structure.+makeConformant :: (Coord t, Compat i) => [NArray i t] -> [NArray i t]+makeConformant ts =+    case conformable (map dims ts) of+        Just alldims -> map (extend alldims) ts+        Nothing -> error $ "makeConformant with inconsistent dimensions "+                         ++ show (map dims ts)++-- the same version for tuples with possibly different element types+makeConformantT (t1,t2) =+    case conformable [dims t1, dims t2] of+        Just alldims -> (extend alldims t1, extend alldims t2)+        Nothing -> error $ "makeConformantT with inconsistent dimensions "+                         ++ show (dims t1, dims t2)
+ lib/Numeric/LinearAlgebra/Array/Simple.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeSynonymInstances #-}+{-# OPTIONS_HADDOCK hide #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Simple+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Simple multidimensional arrays.+-- Contractions only require equal dimension.+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Simple (+    None(..),+    Array,+    listArray+) where++import Numeric.LinearAlgebra.Array.Internal+import Data.Packed+++instance Show (Idx None) where+    show (Idx n s _t) = show n ++ ":" ++ s++-- | Unespecified coordinate type. Contractions only+-- require equal dimension.+data None = None deriving Eq+++instance Compat None where+    compat d1 d2 = iDim d1 == iDim d2+++-- | Multidimensional array with unespecified coordinate type.+type Array t = NArray None t++instance (Coord t) => Show (Array t) where+    show t | null (dims t) = "scalar "++ show (coords t @>0)+           | otherwise = "listArray "++ show (dims t) ++ " "++ show (toList $ coords t)++-- | Construction of an 'Array' from a list of dimensions and a list of elements in left to right order.+listArray :: (Coord t)+    => [Int] -- ^ dimensions+    -> [t]   -- ^ elements+    -> Array t+listArray ds cs = mkNArray dms (product ds |> (cs ++ repeat 0))+    where dms = zipWith3 Idx ds (map show [1::Int ..]) (repeat None)++
+ lib/Numeric/LinearAlgebra/Array/Util.hs view
@@ -0,0 +1,44 @@+-- {-# LANGUAGE FlexibleInstances, FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Array.Util+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Additional tools for manipulation of multidimensional arrays.+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Array.Util (+    Coord, Compat(..),+    NArray, Idx(..), Name,+    scalar,+    rank, names, size, typeOf, dims, coords,++    rename, (!),++    parts,+    newIndex,++    mapArray, zipArray, (|*|),++    extract, onIndex,++    reorder, (~>),+    formatArray, formatFixed, formatScaled,+    dummyAt, noIdx,+    conformable,+    sameStructure,+    makeConformant,+    basisOf,+    asScalar, asVector, asMatrix,+    fromVector, fromMatrix,+    Container(..),+) where++import Numeric.LinearAlgebra.Array.Internal+import Data.Packed(Container(..))
+ lib/Numeric/LinearAlgebra/Exterior.hs view
@@ -0,0 +1,128 @@+{-# LANGUAGE FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.LinearAlgebra.Exterior+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  experimental+--+-- Exterior Algebra.+--+--+-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Exterior (+    (/\),+    inner,+    leviCivita,+    dual,+    (\/),+    module Numeric.LinearAlgebra.Tensor,+    asMultivector, fromMultivector+) where++import Numeric.LinearAlgebra.Tensor+import Numeric.LinearAlgebra.Array.Internal+import Numeric.LinearAlgebra.Multivector(Multivector,fromTensor,maxDim,grade)+import qualified Numeric.LinearAlgebra.Multivector as MV+import Data.List++-- import Debug.Trace+-- debug x = trace (show x) x++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++gsym f t = mkNArray (dims t) (coords $ sum ts) where+    ns = map show [1 .. rank t]+    t' = cov $ renameRaw t ns+    per = permutations ns+    ts  = map (flip renameRaw ns . f . flip reorder t') per++-- symmetrize t = gsym id t++antisymmetrize t = gsym scsig t+    where scsig x = scalar (signature (names x)) * x++fact n = product [1..n]++wedge a b = antisymmetrize (a*b) * (recip . fromIntegral) (fact (rank a) * fact (rank b))++infixl 5 /\+-- | The exterior (wedge) product of two tensors. Obtains the union of subspaces.+--+--   Implemented as the antisymmetrization of the tensor product.+(/\) :: (Coord t)+     => Tensor t+     -> Tensor t+     -> Tensor t+a /\ b = renseq (wedge a' b')+    where a' = renseq  a+          b' = renseq' b++-- levi n = antisymmetrize $ product $ zipWith renameRaw ts is+--     where is = map (return.show) [1 .. n]+--           ts = map (listTensor [n]) (toLists $ ident n)++levi n = listTensor (replicate n n) $ map signature $ sequence (replicate n [1..n])++-- | The full antisymmetric tensor of rank n (contravariant version).+leviCivita :: Int -> Tensor Double+leviCivita = (map levi [0..] !!)++infixl 4 \/+-- | The \"meet\" operator. Obtains the intersection of subspaces.+--+-- @a \\\/ b = dual (dual a \/\\ dual b)@+(\/) :: Tensor Double -> Tensor Double -> Tensor Double+a \/ b = dual (dual a /\ dual b)++dual' n t = inner (leviCivita n) t++-- | Inner product of a r-vector with the whole space.+--+-- @dual t = inner (leviCivita n) t@+dual :: Tensor Double -> Tensor Double+dual t | isScalar t = error $ "cannot deduce dimension for dual of a scalar. Use s * leviCivita n"+       | otherwise  = dual' n t+    where n = case common iDim (dims t) of+                Just x -> x+                Nothing -> error $ "dual with different dimensions"++-- | Euclidean inner product of multivectors.+inner :: (Coord t)+      => Tensor t+      -> Tensor t+      -> Tensor t+inner a b | rank a < rank b = switch (renseq a) * renseq b * k+          | otherwise       = renseq a * switch (renseq b) * k+    where k = recip . fromIntegral $ fact $ min (rank a) (rank b)++renseq t = renameRaw t (map show [1..rank t])+renseq' t = renameRaw t (map ((' ':).show) [1..rank t])++isScalar = null . dims++-- | Extract a compact multivector representation from a full antisymmetric tensor.+--+-- asMultivector = Multivector.'fromTensor'.+--+-- (We do not check that the tensor is actually antisymmetric.)+asMultivector :: Tensor Double -> Multivector+asMultivector = fromTensor++-- | Create an explicit antisymmetric 'Tensor' from the components of a Multivector of a given grade.+fromMultivector :: Int -> Multivector -> Tensor Double+fromMultivector k t = sum $ map f (MV.coords $ grade k t) where+    f (x,es) = scalar x * foldl1' (/\) (map g es)+    n = maxDim t+    g i = vector $ replicate (i-1) 0 ++ 1 : replicate (n-i) 0
+ lib/Numeric/LinearAlgebra/Multivector.hs view
@@ -0,0 +1,265 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.LinearAlgebra.Multivector+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  experimental+--+-- A simple implementation of Geometric Algebra.+--+-- The Num instance provides the geometric product, and the Fractional+-- instance provides the inverse of multivectors.+--+-- This module provides a simple Euclidean embedding.++-----------------------------------------------------------------------------++module Numeric.LinearAlgebra.Multivector (+    Multivector, coords,+    scalar, vector, e, (/\), (-|), (\/), rever, full, rotor,+    apply,+    grade, maxGrade, maxDim,+    fromTensor+) where++import Numeric.LinearAlgebra(toList,reshape,(<\>),(@>))+import Numeric.LinearAlgebra.Array.Internal hiding (scalar,coords)+import Numeric.LinearAlgebra.Tensor hiding (scalar,vector)+import qualified Numeric.LinearAlgebra.Array.Internal as Array+import Data.List+import Control.Monad(filterM)+import Data.Function(on)+import qualified Data.Map as Map++powerset = filterM (const [True, False])  -- !!++base :: Int -> [[Int]]+base k = sortBy (compare `on` length) (powerset [1..k])++base' k = map (\b -> MV [(1,b)]) (base k)++data Multivector = MV { coords :: [(Double,[Int])] } deriving Eq++instance Show Multivector where+    show = showMV++maxGrade :: Multivector -> Int+maxGrade (MV l) = maximum . map (length.snd) $ l++grade :: Int -> Multivector -> Multivector+grade k (MV l) = MV $ filter ((k==).length.snd) l++maxDim :: Multivector -> Int+maxDim (MV [(_,[])]) = 0+maxDim (MV l) = maximum . concat . map snd $ l++-- | The reversion operator.+rever :: Multivector -> Multivector+rever (MV l) = MV (map r l) where+    r (c,b) = (c*fromIntegral s ,b)+        where s = signum (-1)^(k*(k-1)`div`2) :: Int+              k = length b++-- | Show the non zero coordinates of a multivector in a nicer format.+showMV :: Multivector -> String+showMV (MV x) = showBases x++-- | Creates a scalar multivector.+scalar :: Double -> Multivector+scalar s = MV [(s,[])]++-- | Creates a grade 1 multivector of from a list of coordinates.+vector :: [Double] -> Multivector+vector v = MV $ simplify $ zip v (map (:[]) [1..])+++-- different product rules++-- reorders the base indices remembering the original position+r1 :: [Int] -> [(Int,[Int])]+r1 [] = []+r1 l = (m,elemIndices m l):(r1 (filter (/=m) l))+    where m = minimum l++-- geometric product+r2 :: [(Int, [Int])] -> (Double, [Int])+r2 = foldl' g (1,[])+    where g (k,l) (x,ps) = (k*s,l++t)+              where t = if even (length ps) then [] else [x]+                    s = product (map f ps')+                         where f z = if even z then 1 else -1+                               ps' = zipWith (subtract) ps [0..]++-- exterior product+r3 :: [(Int, [Int])] -> (Double, [Int])+r3 = foldl' g (1,[])+    where g (k,l) (x,ps) = (k*s,l++[x])+              where s = if length ps > 1 then 0 else if even (head ps) then 1 else -1+++-- simplification and cleaning of the list of coordinates+simplify = chop . grp . sortBy (compare `on` snd)+    where grp [] = []+          grp [a] = [a]+          grp ((c1,b1):(c2,b2):rest)+              | b1 == b2  = grp ( (c1+c2,b1) : rest)+              | otherwise = (c1,b1): grp ((c2,b2):rest)+          zero (c,_) = abs c < 1E-8+          chop = cz . filter (not.zero)+          cz [] = [(0,[])]+          cz x  = x++-- sum of multivectors+gs (MV l1) (MV l2) = MV $ simplify (l1++l2)++-- geometric product+gp (MV l1) (MV l2) = MV $ simplify [g x y | x<-l1, y <-l2]+    where g (c1,b1) (c2,b2) = (k*c1*c2,b3) where (k,b3) = gpr b1 b2 --(r2.r1) (b1++b2)++-- exterior product+ge (MV l1) (MV l2) = MV $ simplify [g x y | x<-l1, y <-l2]+    where g (c1,b1) (c2,b2) = (k*c1*c2,b3) where (k,b3) = epr b1 b2 -- (r3.r1) (b1++b2)++-- contraction inner product+gi (MV l1) (MV l2) = sum [g x y | x<-l1, y <-l2]+    where g (c1,[]) (c2,is) = MV [(c1*c2,is)]+          g _ (_,[])        = 0++          g (c1,[i]) (c2,[j]) = if i==j then MV [(c1*c2,[])] else 0+          g (c1,[i]) (c2,j:js) = (g (c1,[i]) (c2,[j]) /\ MV [(1,js)])+                               - (MV [(c2,[j])] /\ g (c1,[i]) (1,js))++          g (c1,i:is) b = gi (MV [(c1,[i])]) (gi (MV[(1,is)]) (MV [b]))+++instance Num Multivector where+    (+) = gs+    (*) = gp+    negate (MV l) = MV (map neg l) where neg (k,b) = (-k,b)+    abs _ = error "abs of multivector not yet defined"+    signum _ = error "signum of multivector not yet defined"+    fromInteger x = MV [(fromInteger x,[])]++instance Fractional Multivector where+    fromRational x = MV [(fromRational x,[])]+    recip (MV [(x,[])]) = MV [(recip x,[])]+    recip x = mvrecip x++-- | The k-th basis element.+e :: Int -> Multivector+e k = MV [(1,[k])]++-- | The exterior (outer) product.+(/\) :: Multivector -> Multivector -> Multivector+infixl 7 /\+(/\) = ge+++-- | The contractive inner product.+(-|) :: Multivector -> Multivector -> Multivector+infixl 7 -|+(-|) = gi++-- | The full space of the given dimension. This is the leviCivita simbol, and the basis of the pseudoscalar.+full :: Int -> Multivector+full k = MV [(1,[1 .. k])] --product . map e $ [1 .. k]++-- | Intersection of subspaces.+(\/) :: Multivector -> Multivector -> Multivector+infixl 7 \/+(\/) a b = (b -| rever (full k)) -| a+    where k = max (maxDim a) (maxDim b)++-- check that it is a vector+normVec v = sqrt x where MV [(x,[])] = v * v++unitary v = v / scalar (normVec v)++-- | The rotor operator, used in a sandwich product.+rotor :: Int          -- ^ dimension of the space+      -> Double       -- ^ angle+      -> Multivector  -- ^ axis+      -> Multivector  -- ^ result+rotor k phi axis = scalar (cos (phi/2)) - scalar (sin (phi/2)) * (unitary axis*full k)+++-- memoization of the rules+gprules k = Map.fromList [(x, Map.fromList [(y,(r2.r1)(x++y)) | y<-base k] )| x<-base k]++eprules k = Map.fromList [(x, Map.fromList [(y,(r3.r1)(x++y)) | y<-base k] )| x<-base k]++--reasonable limit+gpr a b = g Map.! a Map.! b+    where g = gprules 6++epr a b = g Map.! a Map.! b+    where g = eprules 6++----------------------- tensor expansion -----------------------++expand k = g+    where g (MV l) = foldl1' (zipWith (+)) $ map f l+          basepos b = m Map.! b+          m = baseraw k+          f (c,b) = en pk (basepos b) c+          pk = 2^k+          baseraw q = Map.fromList $ zip (base q) [0..]+          en n q v = replicate q 0 ++ v : replicate (n-q-1) 0++compact k t = sum $ zipWith (*) (map scalar $ toList (Array.coords t)) (base' k)++gatensor k = listTensor [-pk,-pk,pk] (concat . concat $ gacoords)+    where pk = 2^k+          gacoords = [[ f (x * y) | y<-b] | x<-b]+          b = base' k+          f = expand k++tmv k x = listTensor [2^k] (expand k x)+++-- tp a b = comp (g!"ijk" * ta!"i" * tb!"j")+--     where k = max (maxDim a) (maxDim b)+--           g = gatensor k+--           ta = tmv k a+--           tb = tmv k b+--           comp = compact k++mat rowidx t = reshape c $ Array.coords t'+    where c = iDim $ last (dims t')+          t' = reorder (rowidx: (names t\\[rowidx])) t++-- on the right+pmat k b = mat "k" $ g!"ijk" * tb!"j"+    where g = gatensor k+          tb = tmv k b++divi k a b = compact k $ listTensor [2^k] (toList $ pmat k b <\> Array.coords (tmv k a))++mvrecip b = divi (maxDim b) 1 b++--------------------------------------------------------++-- | Extract a multivector representation from a full antisymmetric tensor.+--+-- (We do not check that the tensor is actually antisymmetric.)+fromTensor :: Tensor Double -> Multivector+fromTensor t = MV $ filter ((/=0.0).fst) $ zip vals basis+    where vals = map ((@> 0). Array.coords .foldl' partF t) (map (map pred) basis)+          r = length (dims t)+          n = iDim . head . dims $ t+          partF s i = part s (name,i) where name = iName . head . dims $ s+          basis = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]++part t (name,i) = parts t name !! i++--------------------------------------------------------++-- | Apply a linear transformation, expressed as the image of the element i-th of the basis.+--+--  (This is a monadic bind!)+apply :: (Int -> Multivector) -> Multivector -> Multivector+apply f t = sum $ map g (coords t) where+    g (x,[]) = scalar x+    g (x,es) = scalar x * foldl1' (/\) (map f es)
+ lib/Numeric/LinearAlgebra/Tensor.hs view
@@ -0,0 +1,102 @@+{-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeSynonymInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.LinearAlgebra.Tensor+-- Copyright   :  (c) Alberto Ruiz 2009+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  experimental+--+-- Tensor computations. Indices can only be contracted if they are of different 'Variant' type.+--+-----------------------------------------------------------------------------+++module Numeric.LinearAlgebra.Tensor (+    -- * The Tensor type+    Tensor, Variant(..),+    listTensor,+    -- * Tensor creation utilities+    superindex, subindex,+    vector, covector, transf,+    -- * Index manipulation+    switch, cov, contrav, forget,+    -- * General array operations+    module Numeric.LinearAlgebra.Array+) where++import Numeric.LinearAlgebra.Array.Internal+import Numeric.LinearAlgebra hiding (rank)+import Numeric.LinearAlgebra.Array++type Tensor t = NArray Variant t++data Variant = Co | Contra deriving (Eq)++instance Compat Variant where+    compat d1 d2 = iDim d1 == iDim d2 && iType d1 /= iType d2++instance Show (Idx Variant) where+    show (Idx n s Co)     = show n ++ "_" ++ s+    show (Idx n s Contra) = show n ++ "^" ++ s++instance (Coord t) => Show (Tensor t) where+    show t | null (dims t) = show (coords t @>0)+           | otherwise = "listTensor " ++ show (dims t) ++ " "++ show (toList (coords t))++flipV Co = Contra+flipV Contra = Co++-- | Creates a tensor from a list of dimensions and a list of coordinates.+-- A positive dimension means that the index is assumed to be contravariant (vector-like), and+-- a negative dimension means that the index is assumed to be covariant (like a linear function, or covector). Contractions can only be performed between indices of different type.+listTensor :: Coord t+           => [Int] -- ^ dimensions+           -> [t]   -- ^ coordinates+           -> Tensor t+listTensor ds cs = mkNArray dms (product ds' |> (cs ++ repeat 0))+    where dms = zipWith3 Idx ds' (map show [1::Int ..]) (map f ds)+          ds' = map abs ds+          f n | n>0       = Contra+              | otherwise = Co++-- | Create an 'Tensor' from a list of parts with a contravariant index (@superindex = 'newIndex' 'Contra'@).+superindex :: Coord t => Name -> [Tensor t] -> Tensor t+superindex = newIndex Contra++-- | Create an 'Tensor' from a list of parts with a covariant index (@subindex = 'newIndex' 'Co'@).+subindex :: Coord t => Name -> [Tensor t] -> Tensor t+subindex   = newIndex Co++++-- | Change the 'Variant' nature of all dimensions to the opposite ones.+switch :: Tensor t -> Tensor t+switch = mapTypes flipV++-- | Make all dimensions covariant.+cov :: NArray i t -> Tensor t+cov     = mapTypes (const Co)++-- | Make all dimensions contravariant.+contrav :: NArray i t -> Tensor t+contrav = mapTypes (const Contra)++-- | Remove the 'Variant' nature of coordinates.+forget :: NArray i t -> Array t+forget = mapTypes (const None)++--------------------------------------------------------------++-- | Create a contravariant rank-1 tensor from a list of coordinates.+vector :: [Double] -> Tensor Double+vector   = fromVector Contra . fromList++-- | Create a covariant rank-1 tensor from a list of coordinates.+covector :: [Double] -> Tensor Double+covector = fromVector Co . fromList++-- | Create a 1-contravariant, 1-covariant rank-2 tensor from list of lists of coordinates.+transf :: [[Double]] -> Tensor Double+transf   = fromMatrix Contra Co . fromLists