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
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 .. order t]
t' = cov $ renameRaw t ns
per = permutations ns
ts = map (flip renameRaw ns . f . flip reorder t') per
antisymmetrize t = gsym scsig t
where scsig x = scalar (signature (namesR x)) * x
fact n = product [1..n]
wedge a b = antisymmetrize (a*b) * (recip . fromIntegral) (fact (order a) * fact (order b))
infixl 5 /\
(/\) :: (Coord t, Fractional t)
=> Tensor t
-> Tensor t
-> Tensor t
a /\ b = renseq (wedge a' b')
where a' = renseq a
b' = renseq' b
levi n = listTensor (replicate n n) $ map signature $ sequence (replicate n [1..n])
leviCivita :: Int -> Tensor Double
leviCivita = (map levi [0..] !!)
infixl 4 \/
(\/) :: Tensor Double -> Tensor Double -> Tensor Double
a \/ b = dual (dual a /\ dual b)
dual' n 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"
inner :: (Coord t, Fractional t)
=> Tensor t
-> Tensor t
-> Tensor t
inner a b | order a < order b = switch (renseq a) * renseq b * k
| otherwise = renseq a * switch (renseq b) * k
where k = recip . fromIntegral $ fact $ min (order a) (order b)
renseq t = renameRaw t (map show [1..order t])
renseq' t = renameRaw t (map ((' ':).show) [1..order t])
isScalar = null . dims
asMultivector :: Tensor Double -> Multivector
asMultivector = fromTensor
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 (i1) 0 ++ 1 : replicate (ni) 0