```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.LinearAlgebra.Multivector
-- Copyright   :  (c) Alberto Ruiz 2009
--
-- 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,
fromTensor
) where

import Numeric.LinearAlgebra(toList,reshape,(<\>),(@>))
import Numeric.LinearAlgebra.Array.Internal hiding (scalar,coords)
import Numeric.LinearAlgebra.Array.Display (showBases)
import Numeric.LinearAlgebra.Tensor hiding (scalar,vector)
import qualified Numeric.LinearAlgebra.Array.Internal as Array
import Data.List
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 (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: (namesR 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)
```