{-# LANGUAGE GADTs #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE TupleSections #-}
module Math.Tensor.LorentzGenerator (
Eta(..), Epsilon(..), Var(..),
AnsatzForestEpsilon, AnsatzForestEta(..),
flattenForest, flattenForestEpsilon, forestEtaList, forestEpsList, forestEtaListLatex, forestEpsListLatex,
drawAnsatzEta, drawAnsatzEpsilon,
getForestLabels, getForestLabelsEpsilon,
removeVarsEta, removeVarsEps,
relabelAnsatzForest, relabelAnsatzForestEpsilon,
mapVars, mapVarsEpsilon,
ansatzRank, ansatzRankEpsilon,
encodeAnsatzForestEta, encodeAnsatzForestEpsilon,
decodeAnsatzForestEta, decodeAnsatzForestEpsilon,
mkAnsatzTensorFastSym, mkAnsatzTensorFast, mkAnsatzTensorFastAbs,
mkAnsatzTensorFastSym', mkAnsatzTensorFast',
mkAnsatzTensorIncrementalSym, mkAnsatzTensorIncremental, mkAnsatzTensorIncrementalAbs,
mkAnsatzTensorIncrementalSym', mkAnsatzTensorIncremental',
Symmetry,
areaList4, areaList6, areaList8, areaList10_1, areaList10_2, areaList12, areaList14_1, areaList14_2,
metricList2, metricList4_1, metricList4_2, metricList6_1, metricList6_2, metricList6_3, metricList8_1, metricList8_2,
symList4, symList6, symList8, symList10_1, symList10_2, symList12, symList14_1, symList14_2,
metricsymList2, metricsymList4_1, metricsymList4_2, metricsymList6_1, metricsymList6_2, metricsymList6_3, metricsymList8_1, metricsymList8_2,
symList16_1, areaList16_1
) where
import qualified Data.IntMap.Strict as I
import qualified Data.Map.Strict as M
import Data.List (nub, permutations, foldl', (\\), elemIndex, nubBy, sortBy, insert, intersect, union, partition, delete, maximumBy, splitAt)
import Data.Maybe (fromJust, isNothing, fromMaybe, isJust, mapMaybe)
import Control.Parallel.Strategies (parListChunk, rdeepseq, runEval, NFData)
import Data.Serialize (encodeLazy, decodeLazy, Serialize(..))
import GHC.Generics
import qualified Data.ByteString.Lazy as BS (ByteString)
import Codec.Compression.GZip (compress, decompress)
import Data.Either (either)
import Data.Tuple (swap)
import GHC.TypeLits
import qualified Numeric.LinearAlgebra.Data as HM
import qualified Numeric.LinearAlgebra as Matrix
import Math.Tensor.Internal.LinearAlgebra (independentColumns)
import Math.Tensor
type Symmetry = ( [(Int,Int)] , [(Int,Int)] , [([Int],[Int])] , [[Int]], [[[Int]]] )
addSym :: Symmetry -> Symmetry -> Symmetry
addSym (a,b,c,d,e) (f,g,h,i,j) = (a `union` f, b `union` g, c `union` h, d `union` i, e `union` j)
mkFilters :: Symmetry -> [(Int,Int)]
mkFilters (pairs,aPairs,blocks,cycles,blockCycles) = map sortPair $ f1 `union` (f2 `union` (f3 `union` f4))
where
sortPair (a,b) = if a < b then (a,b) else (b,a)
f1 = pairs ++ aPairs
f2 = map (\(a,b) -> (head a, head b)) blocks
f3 = concatMap getPairs cycles
f4 = concatMap (getPairs . map head) blockCycles
filter1Sym :: [Int] -> (Int,Int) -> Bool
filter1Sym l (i,j) = case (iPos,jPos) of
(Just i', Just j') -> i' < j'
_ -> True
where
(iPos,jPos) = (elemIndex i l, elemIndex j l)
filterSym :: [Int] -> [(Int,Int)] -> Bool
filterSym l inds = and boolList
where
boolList = map (filter1Sym l) inds
getExtraSyms1 :: [Int] -> Symmetry -> Symmetry
getExtraSyms1 [] _ = ([],[],[],[],[])
getExtraSyms1 [_] _ = error "cannot get extra syms from singleton index list"
getExtraSyms1 (a:b:xs) (pairs,aPairs,blocks,cycles,blockCycles) = addSym (newPairs, [], newBlocks, [], []) (getExtraSyms1 xs newSyms)
where
allBlocks = blocks ++ concatMap mkBlocksFromBlockCycle blockCycles
newBlocks' = map (\(x,y) -> unzip $ filter (\(c,d) -> (c,d) /= (a,b)) $ zip x y) allBlocks
(newBlocks, newPairs') = partition (\(a',_) -> length a' > 1) newBlocks'
newPairs = map (\([a'],[b']) -> (a',b')) newPairs'
newSyms = addSym (pairs,aPairs,blocks,cycles,blockCycles) (newPairs, [], newBlocks, [], [])
mkBlocksFromBlockCycle :: [[Int]] -> [([Int],[Int])]
mkBlocksFromBlockCycle [x,y] = [(x,y)]
mkBlocksFromBlockCycle (x:xs) = l ++ mkBlocksFromBlockCycle xs
where
l = map (x,) xs
mkBlocksFromBlockCycle _ = error "cannot make block symmetries from empty block list"
get2nd :: [Int] -> Symmetry -> (Maybe [(Int,Int)], Maybe [(Int,Int)])
get2nd [a,b] (pairs,aPairs,_,cycles,_) = (sndPairs, sndAPairs)
where
allPairs = pairs ++ concatMap mkSymsFromCycle cycles
aPair = lookup a allPairs
bPair = lookup b (map swap allPairs)
aAPair = lookup a aPairs
bAPair = lookup b (map swap aPairs)
sndPairs = case (aPair, bPair) of
(Nothing, Nothing) -> Nothing
(Just x, Nothing) -> Just [(b,x)]
(Nothing, Just y) -> Just [(a,y)]
(Just x, Just y) -> if x == b then Nothing else Just [(b,x),(a,y)]
sndAPairs = case (aAPair, bAPair) of
(Nothing, Nothing) -> Nothing
(Just x, Nothing) -> Just [(b,x)]
(Nothing, Just y) -> Just [(a,y)]
(Just x, Just y) -> if x == b then Nothing else Just [(b,x),(a,y)]
get2nd _ _ = error "given index list contains more or less than two indices"
get2ndSyms :: Maybe [(Int,Int)] -> Symmetry -> [[Int]] -> Symmetry
get2ndSyms Nothing syms _ = syms
get2ndSyms (Just i) _ etas = (newPairs,[],[],[],[])
where
get2ndInd l (i',j) = mapMaybe (\[a,b] -> if j == a then Just (i',b) else if j == b then Just (i',a) else Nothing) l
newPairs = concatMap (get2ndInd etas) i
mkSymsFromCycle :: [Int] -> [(Int,Int)]
mkSymsFromCycle [x,y] = [(x,y)]
mkSymsFromCycle (x:xs) = l ++ mkSymsFromCycle xs
where
l = map (x,) xs
mkSymsFromCycle _ = error "cannot make syms from empty cycle list"
get2ndASyms :: Maybe [(Int,Int)] -> Symmetry -> [[Int]] -> Symmetry
get2ndASyms Nothing syms _ = syms
get2ndASyms (Just i) _ etas = ([], newAPairs,[],[],[])
where
get2ndInd l (i',j) = mapMaybe (\[a,b] -> if j == a then Just (i',b) else if j == b then Just (i',a) else Nothing) l
newAPairs = concatMap (get2ndInd etas) i
mkEtas :: [Int] -> [[Int]]
mkEtas [] = []
mkEtas [l,k] = [[l,k]]
mkEtas (l:k:ls) = [l,k] : mkEtas ls
mkEtas _ = error "cannot make etas from singleton index list"
getExtraSyms2 :: [Int] -> Symmetry -> Symmetry
getExtraSyms2 [] syms = syms
getExtraSyms2 (a':b':xs) syms = addSym (getExtraSyms2 xs newSyms) newSyms
where
x = [a',b']
(i,j) = get2nd x syms
(p,_,_,_,_) = get2ndSyms i syms (mkEtas xs)
(_,a,_,_,_) = get2ndASyms j syms (mkEtas xs)
newSyms = addSym (p,a,[],[],[]) syms
getExtraSyms2 _ _ = error "cannot get extra syms from singleton index list"
getAllExtraSyms :: [Int] -> Symmetry -> Symmetry
getAllExtraSyms etas syms = allSyms2
where
allSyms1 = addSym (getExtraSyms1 etas syms) syms
allSyms2 = addSym (getExtraSyms2 etas allSyms1) allSyms1
getAllIndsEta :: [Int] -> [(Int,Int)] -> [[Int]]
getAllIndsEta [a,b] _ = [[a,b]]
getAllIndsEta (x:xs) aSyms = concatMap res firstEta
where
firstEta = mapMaybe (\y -> if (x,y) `notElem` aSyms then Just ([x,y],delete y xs) else Nothing) xs
res (a,b) = (++) a <$> getAllIndsEta b aSyms
getAllIndsEta _ _ = error "empty index list"
filterEta :: [Int] -> Symmetry -> [(Int,Int)] -> Bool
filterEta inds (p1,ap1,b1,c1,cb1) filters = filterSym inds totFilters && isNonZero
where
(p2,ap2,b2,c2,cb2) = getAllExtraSyms inds (p1,ap1,b1,c1,cb1)
extrafilters = mkFilters (p2,ap2,b2,c2,cb2)
totFilters = filters `union` extrafilters
etas = map (\[a,b] -> (a,b)) $ mkEtas inds
isNonZero = null $ etas `intersect` union ap1 ap2
getEtaInds :: [Int] -> Symmetry -> [[Int]]
getEtaInds [] _ = [[]]
getEtaInds inds (p,ap,b,c,bc) = filter (\x -> filterEta x (p,ap,b,c,bc) filters1) allInds
where
filters1 = mkFilters (p,ap,b,c,bc)
allInds = getAllIndsEta inds ap
getAllIndsEpsilon :: [Int] -> Symmetry -> [[Int]]
getAllIndsEpsilon inds (p,ap,bs,cyc,_) = [ [a,b,c,d] | a <- [1..numInds-3], b <- [a+1..numInds-2], c <- [b+1..numInds-1], d <- [c+1..numInds],
not (isSym p [a,b,c,d]) && not (is3Area areaBlocks [a,b,c,d]) && isValid2Area areaBlocks [a,b,c,d]
&& not (is1Area areaBlocks [a,b,c,d]) && not (isSymCyc cyc [a,b,c,d]) ]
where
numInds = length inds
blocks2 = filter (\x -> length (fst x) == 2) bs
areaBlocks = map (uncurry (++)) $ filter (\([a,b],[c,d]) -> (a,b) `elem` ap && (c,d) `elem` ap) blocks2
isSym [] _ = False
isSym [(a,b)] [i,j,k,l] = length ([a,b] `intersect` [i,j,k,l]) == 2
isSym (x:xs) [i,j,k,l]
| isSym [x] [i,j,k,l] = True
| otherwise = isSym xs [i,j,k,l]
isSym _ _ = error "expected four indices"
isSymCyc [] _ = False
isSymCyc [l'] [i,j,k,l] = length (l' `intersect` [i,j,k,l]) >= 2
isSymCyc (x:xs) [i,j,k,l]
| isSymCyc [x] [i,j,k,l] = True
| otherwise = isSymCyc xs [i,j,k,l]
isSymCyc _ _ = error "expected four indices"
is3Area [] _ = False
is3Area [[a,b,c,d]] [i,j,k,l] = length ([a,b,c,d] `intersect` [i,j,k,l]) == 3
is3Area (x:xs) [i,j,k,l]
| is3Area [x] [i,j,k,l] = True
| otherwise = is3Area xs [i,j,k,l]
is3Area _ _ = error "expected four indices"
isValid2Area [] _ = True
isValid2Area [[a,b,c,d]] [i,j,k,l]
| length inter == 2 = inter == [a,b]
| otherwise = True
where
inter = [a,b,c,d] `intersect` [i,j,k,l]
isValid2Area (x:xs) [i,j,k,l]
| isValid2Area [x] [i,j,k,l] = isValid2Area xs [i,j,k,l]
| otherwise = False
isValid2Area _ _ = error "expected four indices"
is1Area [] _ = False
is1Area list [i,j,k,l] = maximum (map (length . ([i,j,k,l] `intersect`)) list) == 1
is1Area _ _ = error "expected four indices"
getExtraASymsEps :: [Int] -> Symmetry -> Symmetry
getExtraASymsEps eps (_,_,blo,_,cb) = ([],newASyms, [], [], [])
where
allBlocks = blo ++ concatMap mkBlocksFromBlockCycle cb
blocks2 = filter (\(a,_) -> length a == 2) allBlocks
newASyms = mapMaybe (\([i,j],[k,l]) -> if length ([i,k] `intersect` eps) == 2 then Just (j,l) else if length ([j,l] `intersect` eps) == 2 then Just (i,k) else Nothing) blocks2
getEpsilonInds :: [Int] -> Symmetry -> [[Int]]
getEpsilonInds inds sym = allIndsRed
where
epsInds = getAllIndsEpsilon inds sym
allInds = concat $ filter (not . null) $ map (\x -> map (x ++) $ getEtaInds (inds \\ x) (addSym sym (getExtraASymsEps x sym)) )epsInds
filters = mkFilters sym
allIndsRed = filter (\x -> let symEps = addSym (getExtraASymsEps (take 4 x) sym) sym
symEta = addSym symEps (getAllExtraSyms (drop 4 x) symEps)
newFilters = union filters (mkFilters symEta)
in filterSym x newFilters) allInds
data Eta = Eta {-# UNPACK #-} !Int {-# UNPACK #-} !Int deriving (Show, Read, Eq, Ord, Generic, Serialize, NFData)
data Epsilon = Epsilon {-# UNPACK #-} !Int {-# UNPACK #-} !Int {-# UNPACK #-} !Int {-# UNPACK #-} !Int deriving (Show, Read, Eq, Ord, Generic, Serialize, NFData)
data Var = Var {-# UNPACK #-} !Int {-# UNPACK #-} !Int deriving (Show, Read, Eq, Ord, Generic, Serialize, NFData )
sortList :: Ord a => [a] -> [a]
sortList = foldr insert []
sortEta :: Eta -> Eta
sortEta (Eta x y) = Eta (min x y) (max x y)
{-# INLINEABLE sortEta #-}
sortEpsilon :: Epsilon -> Epsilon
sortEpsilon (Epsilon i j k l) = Epsilon i' j' k' l'
where
[i',j',k',l'] = sortList [i,j,k,l]
getEpsSign :: Epsilon -> Int
getEpsSign (Epsilon i j k l) = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k])
{-# INLINEABLE getEpsSign #-}
addVars :: Var -> Var -> Var
addVars (Var x y) (Var x' y') = if y == y' then Var (x + x') y else error "should not add different vars"
{-# INLINEABLE addVars #-}
multVar :: Int -> Var -> Var
multVar x (Var x' y) = Var (x * x') y
{-# INLINEABLE multVar #-}
isZeroVar :: Var -> Bool
isZeroVar (Var x _) = x==0
{-# INLINEABLE isZeroVar #-}
data AnsatzForestEta = ForestEta (M.Map Eta AnsatzForestEta)| Leaf !Var | EmptyForest deriving (Show, Read, Eq, Generic, Serialize)
type AnsatzForestEpsilon = M.Map Epsilon AnsatzForestEta
encodeAnsatzForestEta :: AnsatzForestEta -> BS.ByteString
encodeAnsatzForestEta = compress . encodeLazy
encodeAnsatzForestEpsilon :: AnsatzForestEpsilon -> BS.ByteString
encodeAnsatzForestEpsilon = compress . encodeLazy
decodeAnsatzForestEta :: BS.ByteString -> AnsatzForestEta
decodeAnsatzForestEta bs = either error id $ decodeLazy $ decompress bs
decodeAnsatzForestEpsilon :: BS.ByteString -> AnsatzForestEpsilon
decodeAnsatzForestEpsilon bs = either error id $ decodeLazy $ decompress bs
mapNodes :: (Eta -> Eta) -> AnsatzForestEta -> AnsatzForestEta
mapNodes _ EmptyForest = EmptyForest
mapNodes f (ForestEta m) = ForestEta $ M.mapKeys f . M.map (mapNodes f) $ m
mapNodes _ (Leaf x) = Leaf x
mapNodesEpsilon :: (Epsilon -> Epsilon) -> AnsatzForestEpsilon -> AnsatzForestEpsilon
mapNodesEpsilon = M.mapKeys
mapVars :: (Var -> Var) -> AnsatzForestEta -> AnsatzForestEta
mapVars _ EmptyForest = EmptyForest
mapVars f (Leaf var) = Leaf (f var)
mapVars f (ForestEta m) = ForestEta $ M.map (mapVars f) m
mapVarsEpsilon :: (Var -> Var) -> AnsatzForestEpsilon -> AnsatzForestEpsilon
mapVarsEpsilon f = M.map (mapVars f)
getLeafVals :: AnsatzForestEta -> [Var]
getLeafVals (Leaf var) = [var]
getLeafVals (ForestEta m) = rest
where
rest = concatMap getLeafVals $ M.elems m
getLeafVals EmptyForest = []
getLeafValsEpsilon :: AnsatzForestEpsilon -> [Var]
getLeafValsEpsilon m = concatMap getLeafVals $ M.elems m
getVarLabels :: Var -> Int
getVarLabels (Var _ j) = j
getForestLabels :: AnsatzForestEta -> [Int]
getForestLabels ans = nub $ map getVarLabels $ getLeafVals ans
getForestLabelsEpsilon :: AnsatzForestEpsilon -> [Int]
getForestLabelsEpsilon m = nub $ map getVarLabels $ getLeafValsEpsilon m
ansatzRank :: AnsatzForestEta -> Int
ansatzRank ans = length $ getForestLabels ans
ansatzRankEpsilon :: AnsatzForestEpsilon -> Int
ansatzRankEpsilon ans = length $ getForestLabelsEpsilon ans
relabelVar :: (Int -> Int) -> Var -> Var
relabelVar f (Var i j) = Var i (f j)
relabelAnsatzForest :: Int -> AnsatzForestEta -> AnsatzForestEta
relabelAnsatzForest i ans = mapVars update ans
where
vars = getForestLabels ans
relabMap = I.fromList $ zip vars [i..]
update = relabelVar ((I.!) relabMap)
removeVarsEta :: [Int] -> AnsatzForestEta -> AnsatzForestEta
removeVarsEta vars (Leaf (Var i j))
| j `elem` vars = EmptyForest
| otherwise = Leaf (Var i j)
removeVarsEta vars (ForestEta m) = ForestEta $ M.filter (/= EmptyForest) $ M.map (removeVarsEta vars) m
removeVarsEta _ EmptyForest = EmptyForest
relabelAnsatzForestEpsilon :: Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon
relabelAnsatzForestEpsilon i ans = if ans == M.empty then M.empty else mapVarsEpsilon update ans
where
vars = getForestLabelsEpsilon ans
relabMap = I.fromList $ zip vars [i..]
update = relabelVar ((I.!) relabMap)
removeVarsEps :: [Int] -> AnsatzForestEpsilon -> AnsatzForestEpsilon
removeVarsEps vars m = M.filter (/= EmptyForest) $ M.map (removeVarsEta vars) m
addForests :: AnsatzForestEta -> AnsatzForestEta -> AnsatzForestEta
addForests ans EmptyForest = ans
addForests EmptyForest ans = ans
addForests (Leaf var1) (Leaf var2)
| isZeroVar newLeafVal = EmptyForest
| otherwise = Leaf newLeafVal
where
newLeafVal = addVars var1 var2
addForests (ForestEta m1) (ForestEta m2)
| M.null newMap = EmptyForest
| otherwise = ForestEta newMap
where
newMap = M.filter (/= EmptyForest) $ M.unionWith addForests m1 m2
addForests _ _ = error "cannot add Leaf and Forest"
addForestsEpsilon :: AnsatzForestEpsilon -> AnsatzForestEpsilon -> AnsatzForestEpsilon
addForestsEpsilon m1 m2 = M.filter (/= EmptyForest) $ M.unionWith addForests m1 m2
addList2Forest :: AnsatzForestEta -> ([Eta],Var) -> AnsatzForestEta
addList2Forest EmptyForest x = mkForestFromAscList x
addList2Forest (Leaf var1) ([], var2)
| isZeroVar newLeafVal = EmptyForest
| otherwise = Leaf newLeafVal
where
newLeafVal = addVars var1 var2
addList2Forest (ForestEta m1) (x:xs, var) = ForestEta $ M.insertWith (\_ a2 -> addList2Forest a2 (xs, var)) x newVal m1
where
newVal = mkForestFromAscList (xs,var)
addList2Forest (ForestEta f) ([], _) = ForestEta f
addList2Forest (Leaf _) _ = error "cannot add something to Leaf"
addList2ForestEpsilon :: AnsatzForestEpsilon -> (Epsilon,[Eta],Var) -> AnsatzForestEpsilon
addList2ForestEpsilon m (eps,eta,var) = M.insertWith (\_ a2 -> addList2Forest a2 (eta, var)) eps newVal m
where
newVal = mkForestFromAscList (eta,var)
flattenForest :: AnsatzForestEta -> [([Eta],Var)]
flattenForest EmptyForest = []
flattenForest (Leaf var) = [([],var)]
flattenForest (ForestEta m) = concat l
where
mPairs = M.assocs m
l = fmap (\(k,v) -> map (\(i,j) -> (insert k i, j)) $ flattenForest v) mPairs
flattenForestEpsilon :: AnsatzForestEpsilon -> [(Epsilon,[Eta],Var)]
flattenForestEpsilon m = concat l
where
mPairs = M.assocs m
l = fmap (\(k,v) -> map (\(i,j) -> (k, i, j)) $ flattenForest v) mPairs
drawEtaTree :: Eta -> AnsatzForestEta -> [String]
drawEtaTree (Eta i j) (Leaf (Var a b)) = ["(" ++ show i ++ "," ++ show j ++ ") * (" ++ show a ++ ") * x[" ++ show b ++ "]"]
drawEtaTree (Eta i j) (ForestEta m) = lines ("(" ++ show i ++ "," ++ show j ++ ")") ++ drawSubTrees m
where
drawSubTrees x
| x == M.empty = []
| M.size x == 1 = let [(a,b)] = M.assocs x in "|" : shift "`---- " " " (drawEtaTree a b)
| otherwise = let (a,b) = head $ M.assocs x in "|" : shift "+---- " "| " (drawEtaTree a b) ++ drawSubTrees (M.delete a x)
shift first other = zipWith (++) (first : repeat other)
drawEtaTree _ EmptyForest = []
drawEpsilonTree :: Epsilon -> AnsatzForestEta -> [String]
drawEpsilonTree (Epsilon i j k l) (Leaf (Var a b)) = ["(" ++ show i ++ "," ++ show j ++ "," ++ show k ++ "," ++ show l ++ ") * (" ++ show a ++ ") * x[" ++ show b ++ "]"]
drawEpsilonTree (Epsilon i j k l) (ForestEta m) = lines ("(" ++ show i ++ "," ++ show j ++ "," ++ show k ++ "," ++ show l ++ ")") ++ drawSubTrees m
where
drawSubTrees x
| x == M.empty = []
| M.size x == 1 = let [(a,b)] = M.assocs x in "|" : shift "`---- " " " (drawEtaTree a b)
| otherwise = let (a,b) = head $ M.assocs x in "|" : shift "+---- " "| " (drawEtaTree a b) ++ drawSubTrees (M.delete a x)
shift first other = zipWith (++) (first : repeat other)
drawEpsilonTree _ EmptyForest = []
drawAnsatzEta :: AnsatzForestEta -> String
drawAnsatzEta (Leaf (Var a b)) = show a ++ "x[" ++ show b ++ "]"
drawAnsatzEta (ForestEta m) = unlines $ map (\(x,y) -> unlines $ drawEtaTree x y) $ M.assocs m
drawAnsatzEta EmptyForest = []
drawAnsatzEpsilon :: AnsatzForestEpsilon -> String
drawAnsatzEpsilon m
| M.size m == 0 = []
| otherwise = unlines $ map (\(x,y) -> unlines $ drawEpsilonTree x y) $ M.assocs m
forestEtaList :: AnsatzForestEta -> [[Eta]]
forestEtaList f = map fst fList''
where
fList = flattenForest f
fList' = sortBy (\(_, Var _ y1 ) (_, Var _ y2) -> compare y1 y2) fList
fList'' = nubBy (\(_, Var x1 y1 ) (_, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList'
forestEpsList :: AnsatzForestEpsilon -> [(Epsilon,[Eta])]
forestEpsList f = map (\(a,b,_) -> (a,b)) fList''
where
fList = flattenForestEpsilon f
fList' = sortBy (\(_, _, Var _ y1 ) (_, _, Var _ y2) -> compare y1 y2) fList
fList'' = nubBy (\(_, _, Var x1 y1 ) (_, _, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList'
mkEtasLatex :: String -> Eta -> String
mkEtasLatex inds (Eta i j) = "\\eta^{" ++ etaI : etaJ : "}"
where
(etaI,etaJ) = (inds !! (i-1), inds !! (j-1) )
forestEtaListLatex :: AnsatzForestEta -> String -> Char -> String
forestEtaListLatex f inds var = tail $ concat etaL''
where
etaL = sortBy (\(_, Var _ y1 ) (_, Var _ y2) -> compare y1 y2) $ flattenForest f
etaL' = nubBy (\(_, Var x1 y1 ) (_, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) etaL
etaL'' = map (\(a,Var _ y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ concatMap (mkEtasLatex inds) a) etaL'
mkEpsLatex :: String -> Epsilon -> String
mkEpsLatex inds (Epsilon i j k l) = "\\epsilon^{" ++ epsi : epsj : epsk : epsl : "}"
where
(epsi, epsj, epsk, epsl) = (inds !! (i-1), inds !! (j-1), inds !! (k-1), inds !! (l-1))
forestEpsListLatex :: AnsatzForestEpsilon -> String -> Char -> String
forestEpsListLatex f inds var = tail $ concat epsL''
where
epsL = sortBy (\(_, _, Var _ y1 ) (_, _, Var _ y2) -> compare y1 y2) $ flattenForestEpsilon f
epsL' = nubBy (\(_, _, Var x1 y1 ) (_, _, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) epsL
epsL'' = map (\(a,b,Var _ y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ mkEpsLatex inds a ++ concatMap (mkEtasLatex inds) b) epsL'
mkForestFromAscList :: ([Eta],Var) -> AnsatzForestEta
mkForestFromAscList ([],var) = Leaf var
mkForestFromAscList (x:xs, var) = ForestEta $ M.singleton x $ mkForestFromAscList (xs,var)
mkForestFromAscListEpsilon :: (Epsilon,[Eta],Var) -> AnsatzForestEpsilon
mkForestFromAscListEpsilon (x,y,z) = M.singleton x $ mkForestFromAscList (y,z)
canonicalizeAnsatzEta :: AnsatzForestEta -> AnsatzForestEta
canonicalizeAnsatzEta = mapNodes sortEta
canonicalizeAnsatzEpsilon :: AnsatzForestEpsilon -> AnsatzForestEpsilon
canonicalizeAnsatzEpsilon m = newMap
where
newMap = M.mapKeys sortEpsilon $ M.mapWithKey (\k v -> mapVars (multVar (getEpsSign k) ) v) $ M.map (mapNodes sortEta) m
sortForest :: AnsatzForestEta -> AnsatzForestEta
sortForest f = foldl' addList2Forest EmptyForest fList
where
fList = flattenForest f
sortForestEpsilon :: AnsatzForestEpsilon -> AnsatzForestEpsilon
sortForestEpsilon f = foldl' addList2ForestEpsilon M.empty fList
where
fList = flattenForestEpsilon f
swapLabelF :: (Int,Int) -> Int -> Int
swapLabelF (x,y) z
| x == z = y
| y == z = x
| otherwise = z
swapBlockLabelMap :: ([Int],[Int]) -> I.IntMap Int
swapBlockLabelMap (x,y) = swapF
where
swapF = I.fromList $ zip x y ++ zip y x
swapLabelEta :: (Int,Int) -> Eta -> Eta
swapLabelEta inds (Eta x y) = Eta (f x) (f y)
where
f = swapLabelF inds
swapLabelEpsilon :: (Int,Int) -> Epsilon -> Epsilon
swapLabelEpsilon inds (Epsilon i j k l) = Epsilon (f i) (f j) (f k) (f l)
where
f = swapLabelF inds
swapBlockLabelEta :: I.IntMap Int -> Eta -> Eta
swapBlockLabelEta swapF (Eta i j) = Eta i' j'
where
i' = I.findWithDefault i i swapF
j' = I.findWithDefault j j swapF
swapBlockLabelEpsilon :: I.IntMap Int -> Epsilon -> Epsilon
swapBlockLabelEpsilon swapF (Epsilon i j k l) = Epsilon i' j' k' l'
where
i' = I.findWithDefault i i swapF
j' = I.findWithDefault j j swapF
k' = I.findWithDefault k k swapF
l' = I.findWithDefault l l swapF
swapLabelFEta :: (Int,Int) -> AnsatzForestEta -> AnsatzForestEta
swapLabelFEta inds ans = sortForest.canonicalizeAnsatzEta $ swapAnsatz
where
f = swapLabelEta inds
swapAnsatz = mapNodes f ans
swapLabelFEps :: (Int,Int) -> AnsatzForestEpsilon -> AnsatzForestEpsilon
swapLabelFEps inds ans = sortForestEpsilon.canonicalizeAnsatzEpsilon $ swapAnsatz
where
f = swapLabelEpsilon inds
swapAnsatz = mapNodesEpsilon f $ M.map (swapLabelFEta inds) ans
swapBlockLabelFEta :: I.IntMap Int -> AnsatzForestEta -> AnsatzForestEta
swapBlockLabelFEta swapF ans = sortForest.canonicalizeAnsatzEta $ swapAnsatz
where
f = swapBlockLabelEta swapF
swapAnsatz = mapNodes f ans
swapBlockLabelFEps :: I.IntMap Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon
swapBlockLabelFEps swapF ans = sortForestEpsilon.canonicalizeAnsatzEpsilon $ swapAnsatz
where
f = swapBlockLabelEpsilon swapF
swapAnsatz = mapNodesEpsilon f $ M.map (swapBlockLabelFEta swapF) ans
pairSymForestEta :: (Int,Int) -> AnsatzForestEta -> AnsatzForestEta
pairSymForestEta inds ans = addForests ans $ swapLabelFEta inds ans
pairSymForestEps :: (Int,Int) -> AnsatzForestEpsilon -> AnsatzForestEpsilon
pairSymForestEps inds ans = addForestsEpsilon ans $ swapLabelFEps inds ans
pairASymForestEta :: (Int,Int) -> AnsatzForestEta -> AnsatzForestEta
pairASymForestEta inds ans = addForests ans $ mapVars (multVar (-1)) $ swapLabelFEta inds ans
pairASymForestEps :: (Int,Int) -> AnsatzForestEpsilon -> AnsatzForestEpsilon
pairASymForestEps inds ans = addForestsEpsilon ans $ mapVarsEpsilon (multVar (-1)) $ swapLabelFEps inds ans
pairBlockSymForestEta :: I.IntMap Int -> AnsatzForestEta -> AnsatzForestEta
pairBlockSymForestEta swapF ans = addForests ans $ swapBlockLabelFEta swapF ans
pairBlockSymForestEps :: I.IntMap Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon
pairBlockSymForestEps swapF ans = addForestsEpsilon ans $ swapBlockLabelFEps swapF ans
pairBlockASymForestEta :: I.IntMap Int -> AnsatzForestEta -> AnsatzForestEta
pairBlockASymForestEta swapF ans = addForests ans $ mapVars (multVar (-1)) $ swapBlockLabelFEta swapF ans
pairBlockASymForestEps :: I.IntMap Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon
pairBlockASymForestEps swapF ans = addForestsEpsilon ans $ mapVarsEpsilon (multVar (-1)) $ swapBlockLabelFEps swapF ans
cyclicSymForestEta :: [Int] -> AnsatzForestEta -> AnsatzForestEta
cyclicSymForestEta inds ans = foldr (\y x -> addForests x $ swapBlockLabelFEta y ans ) ans perms
where
perms = map (I.fromList . zip inds) $ tail $ permutations inds
cyclicSymForestEps :: [Int] -> AnsatzForestEpsilon -> AnsatzForestEpsilon
cyclicSymForestEps inds ans = foldr (\y x -> addForestsEpsilon x $ swapBlockLabelFEps y ans ) ans perms
where
perms = map (I.fromList . zip inds) $ tail $ permutations inds
cyclicBlockSymForestEta :: [[Int]] -> AnsatzForestEta -> AnsatzForestEta
cyclicBlockSymForestEta inds ans = foldr (\y x -> addForests x $ swapBlockLabelFEta y ans ) ans perms
where
perms = map (I.fromList . zip (concat inds) . concat) $ tail $ permutations inds
cyclicBlockSymForestEps :: [[Int]] -> AnsatzForestEpsilon-> AnsatzForestEpsilon
cyclicBlockSymForestEps inds ans = foldr (\y x -> addForestsEpsilon x $ swapBlockLabelFEps y ans ) ans perms
where
perms = map (I.fromList . zip (concat inds) . concat) $ tail $ permutations inds
symAnsatzForestEta ::Symmetry -> AnsatzForestEta -> AnsatzForestEta
symAnsatzForestEta (sym,asym,blocksym,cyclicsym,cyclicblocksym) ans =
foldr cyclicBlockSymForestEta (
foldr cyclicSymForestEta (
foldr pairBlockSymForestEta (
foldr pairASymForestEta (
foldr pairSymForestEta ans sym
) asym
) blockSymMap
) cyclicsym
) cyclicblocksym
where
blockSymMap = map swapBlockLabelMap blocksym
symAnsatzForestEps :: Symmetry -> AnsatzForestEpsilon -> AnsatzForestEpsilon
symAnsatzForestEps (sym,asym,blocksym,cyclicsym,cyclicblocksym) ans =
foldr cyclicBlockSymForestEps (
foldr cyclicSymForestEps (
foldr pairBlockSymForestEps (
foldr pairASymForestEps (
foldr pairSymForestEps ans sym
) asym
) blockSymMap
) cyclicsym
) cyclicblocksym
where
blockSymMap = map swapBlockLabelMap blocksym
mkEtaList :: [Int] -> [Eta]
mkEtaList [] = []
mkEtaList x = Eta a b : mkEtaList rest
where
[a,b] = take 2 x
rest = drop 2 x
mkEpsilonList :: [Int] -> (Epsilon,[Eta])
mkEpsilonList x = (Epsilon i j k l , mkEtaList rest)
where
[i,j,k,l] = take 4 x
rest = drop 4 x
mkEtaList' :: Var -> [Int] -> ([Eta],Var)
mkEtaList' var l = (mkEtaList l, var)
mkEpsilonList' :: Var -> [Int] -> (Epsilon,[Eta],Var)
mkEpsilonList' var l = (eps, eta, var)
where
(eps,eta) = mkEpsilonList l
isElem :: [Eta] -> AnsatzForestEta -> Bool
isElem [] (Leaf _) = True
isElem _ EmptyForest = False
isElem (x:xs) (ForestEta m) = case mForest of
Just forest -> xs `isElem` forest
_ -> False
where
mForest = M.lookup x m
isElem _ (Leaf _) = error "cannot lookup eta in Leaf"
isElem [] (ForestEta _) = error "cannot look for eta in forest when no eta is specified"
isElemEpsilon :: (Epsilon, [Eta]) -> AnsatzForestEpsilon -> Bool
isElemEpsilon (eps,l) m = case mForest of
Just forest -> l `isElem` forest
_ -> False
where
mForest = M.lookup eps m
reduceAnsatzEta' :: Symmetry -> [([Eta],Var)] -> AnsatzForestEta
reduceAnsatzEta' sym = foldl' addOrRem' EmptyForest
where
addOrRem' f ans = if isElem (fst ans) f then f else addForests f (symAnsatzForestEta sym $ mkForestFromAscList ans)
reduceAnsatzEpsilon' :: Symmetry -> [(Epsilon, [Eta], Var)] -> AnsatzForestEpsilon
reduceAnsatzEpsilon' sym = foldl' addOrRem' M.empty
where
addOrRem' f (x,y,z) = if isElemEpsilon (x,y) f then f else addForestsEpsilon f (symAnsatzForestEps sym $ mkForestFromAscListEpsilon (x,y,z))
mkAllVars :: [Var]
mkAllVars = map (Var 1) [1..]
getEtaForestFast :: Int -> Symmetry -> AnsatzForestEta
getEtaForestFast ord syms = relabelAnsatzForest 1 $ reduceAnsatzEta' syms allForests
where
allInds = getEtaInds [1..ord] syms
allVars = mkAllVars
allForests = zipWith mkEtaList' allVars allInds
getEpsForestFast :: Int -> Symmetry -> AnsatzForestEpsilon
getEpsForestFast ord syms = if ord < 4 then M.empty else relabelAnsatzForestEpsilon 1 $ reduceAnsatzEpsilon' syms allForests
where
allInds = getEpsilonInds [1..ord] syms
allVars = mkAllVars
allForests = zipWith mkEpsilonList' allVars allInds
evalNodeEta :: I.IntMap Int -> Eta -> Maybe Int
evalNodeEta iMap (Eta x y)
| a == b && a == 0 = Just (-1)
| a == b = Just 1
| otherwise = Nothing
where
[a,b] = [(I.!) iMap x, (I.!) iMap y]
evalNodeEpsilon :: I.IntMap Int -> Epsilon -> Maybe Int
evalNodeEpsilon iMap (Epsilon w x y z) = M.lookup l epsMap
where
l = [(I.!) iMap w, (I.!) iMap x, (I.!) iMap y, (I.!) iMap z]
epsMap :: M.Map [Int] Int
epsMap = M.fromList $ map (\x@[i,j,k,l] -> (x, epsSign i j k l)) $ permutations [0,1,2,3]
where
epsSign i j k l = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k])
evalAnsatzForestEta :: I.IntMap Int -> AnsatzForestEta -> I.IntMap Int
evalAnsatzForestEta _ (Leaf (Var x y)) = I.singleton y x
evalAnsatzForestEta evalM (ForestEta m) = M.foldlWithKey' foldF I.empty m
where
foldF b k a = let nodeVal = evalNodeEta evalM k
in if isNothing nodeVal then b
else I.unionWith (+) (I.map (fromJust nodeVal *) (evalAnsatzForestEta evalM a)) b
evalAnsatzForestEta _ EmptyForest = I.empty
evalAnsatzForestEpsilon :: I.IntMap Int -> AnsatzForestEpsilon -> I.IntMap Int
evalAnsatzForestEpsilon evalM = M.foldlWithKey' foldF I.empty
where
foldF b k a = let nodeVal = evalNodeEpsilon evalM k
in if isNothing nodeVal then b
else I.unionWith (+) (I.map (fromJust nodeVal *) (evalAnsatzForestEta evalM a)) b
eval1AnsatzForestEta :: I.IntMap Int -> AnsatzForestEta -> Int
eval1AnsatzForestEta _ (Leaf (Var x _)) = x
eval1AnsatzForestEta evalM (ForestEta m) = M.foldlWithKey' foldF 0 m
where
foldF b k a = let nodeVal = evalNodeEta evalM k
in if isNothing nodeVal then b
else b + (fromJust nodeVal * eval1AnsatzForestEta evalM a)
eval1AnsatzForestEta _ EmptyForest = 0
eval1AnsatzForestEpsilon :: I.IntMap Int -> AnsatzForestEpsilon -> Int
eval1AnsatzForestEpsilon evalM = M.foldlWithKey' foldF 0
where
foldF b k a = let nodeVal = evalNodeEpsilon evalM k
in if isNothing nodeVal then b
else b + (fromJust nodeVal * eval1AnsatzForestEta evalM a)
mkVecList :: (Foldable t, NFData a1, Real a1) =>
(a2 -> Maybe ((Int, Int), a1)) -> [a2] -> t a3 -> Maybe (HM.Matrix Double)
mkVecList mkAns dofList evalM = if null l
then Nothing
else Just $ HM.assoc (1,n) 0 $ map (fmap (\x -> fromRational $ toRational x / toRational maxVal)) l
where
l' = mapMaybe mkAns dofList
l = runEval $ parListChunk 500 rdeepseq l'
lVals = map (\((_,_),z) -> z) l
maxVal = maximum lVals
n = length evalM
evalAnsatzEtaVecListIncremental :: [I.IntMap Int] -> AnsatzForestEta -> Maybe (HM.Matrix Double)
evalAnsatzEtaVecListIncremental _ EmptyForest = Nothing
evalAnsatzEtaVecListIncremental evalM f = mkVecList mkAns dofList evalM
where
dofList = zip [0..] evalM
mkAns (i,j) = let ansVal = eval1AnsatzForestEta j f
in if ansVal == 0 then Nothing else Just ((0,i), ansVal)
evalAnsatzEpsilonVecListIncremental :: [I.IntMap Int] -> AnsatzForestEpsilon -> Maybe (HM.Matrix Double)
evalAnsatzEpsilonVecListIncremental evalM f = if f == M.empty then Nothing else mkVecList mkAns dofList evalM
where
dofList = zip [0..] evalM
mkAns (i,j) = let ansVal = eval1AnsatzForestEpsilon j f
in if ansVal == 0 then Nothing else Just ((0,i), ansVal)
type AssocsList a = [([(Int,Int)],a)]
type AssocsListAbs a = [([(Int,Int)],Int,a)]
evalAllEta :: [I.IntMap Int] -> AnsatzForestEta -> [[(Int,Int)]]
evalAllEta [] _ = []
evalAllEta _ EmptyForest = []
evalAllEta evalMs f = l'
where
l = map (\x -> filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
evalAllTensorEta :: (NFData a) => [(I.IntMap Int, a)] -> AnsatzForestEta -> AssocsList a
evalAllTensorEta [] _ = []
evalAllTensorEta _ EmptyForest = []
evalAllTensorEta evalMs f = l'
where
l = map (\(x,z) -> (filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f,z)) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
evalAllEpsilon :: [I.IntMap Int] -> AnsatzForestEpsilon -> [[(Int,Int)]]
evalAllEpsilon [] _ = []
evalAllEpsilon evalMs f = if f == M.empty then [] else l'
where
l = map (\x -> filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
evalAllTensorEpsilon :: (NFData a) => [(I.IntMap Int, a)] -> AnsatzForestEpsilon -> AssocsList a
evalAllTensorEpsilon [] _ = []
evalAllTensorEpsilon evalMs f = if f == M.empty then [] else l'
where
l = map (\(x,z) -> ( filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f,z)) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
evalAllTensorEtaAbs :: (NFData a) => [(I.IntMap Int, Int, a)] -> AnsatzForestEta -> AssocsListAbs a
evalAllTensorEtaAbs [] _ = []
evalAllTensorEtaAbs _ EmptyForest = []
evalAllTensorEtaAbs evalMs f = l'
where
l = map (\(x,y,z) -> (filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f, y,z)) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
evalAllTensorEpsilonAbs :: (NFData a) => [(I.IntMap Int, Int, a)] -> AnsatzForestEpsilon -> AssocsListAbs a
evalAllTensorEpsilonAbs [] _ = []
evalAllTensorEpsilonAbs evalMs f = if f == M.empty then [] else l'
where
l = map (\(x,y,z) -> ( filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f, y,z)) evalMs
l' = runEval $ parListChunk 500 rdeepseq l
type RankDataIncremental = (HM.Matrix Double, HM.Matrix Double)
getVarNrIncremental :: RankDataIncremental -> Int
getVarNrIncremental = HM.rows . snd
checkNumericLinDepIncremental :: RankDataIncremental -> Maybe (HM.Matrix Double) -> Maybe RankDataIncremental
checkNumericLinDepIncremental (lastMat, lastFullMat) (Just newVec)
| rk < maxRank = Nothing
| otherwise = Just (newMat, newAnsatzMat)
where
newVecTrans = HM.tr newVec
scalar = newVec Matrix.<> newVecTrans
prodBlock = lastFullMat Matrix.<> newVecTrans
prodBlockTrans = HM.tr prodBlock
newMat = HM.fromBlocks [[lastMat, prodBlock],
[prodBlockTrans, scalar ]]
rk = Matrix.rank newMat
maxRank = min (HM.cols newMat) (HM.rows newMat)
newAnsatzMat = lastFullMat HM.=== newVec
checkNumericLinDepIncremental _ Nothing = Nothing
getNewRDat :: [I.IntMap Int] -> AnsatzForestEta -> RankDataIncremental -> Maybe RankDataIncremental
getNewRDat evalM newAns rDat = newRDat
where
newVec = evalAnsatzEtaVecListIncremental evalM newAns
newRDat = checkNumericLinDepIncremental rDat newVec
getNewRDatEps :: [I.IntMap Int] -> AnsatzForestEpsilon -> RankDataIncremental -> Maybe RankDataIncremental
getNewRDatEps evalM newAns rDat = newRDat
where
newVec = evalAnsatzEpsilonVecListIncremental evalM newAns
newRDat = checkNumericLinDepIncremental rDat newVec
getNewAns :: Symmetry -> [Eta] -> RankDataIncremental -> AnsatzForestEta
getNewAns symList etaList rDat = symAnsatzForestEta symList $ mkForestFromAscList (etaList,Var 1 (getVarNrIncremental rDat + 1))
getNewAnsEps :: Symmetry -> Epsilon -> [Eta] -> RankDataIncremental -> AnsatzForestEpsilon
getNewAnsEps symList epsList etaList rDat = symAnsatzForestEps symList $ mkForestFromAscListEpsilon (epsList,etaList,Var 1 (getVarNrIncremental rDat + 1))
addOrDiscardEtaIncremental :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEta, RankDataIncremental) -> [Eta] -> (AnsatzForestEta, RankDataIncremental)
addOrDiscardEtaIncremental symList evalM (ans,rDat) etaL
| isElem etaL ans = (ans,rDat)
| otherwise = case newRDat of
Nothing -> (ans,rDat)
Just newRDat' -> (sumAns,newRDat')
where
newAns = getNewAns symList etaL rDat
newRDat = getNewRDat evalM newAns rDat
sumAns = addForests ans newAns
addOrDiscardEpsilonIncremental :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon, RankDataIncremental) -> (Epsilon,[Eta]) -> (AnsatzForestEpsilon, RankDataIncremental)
addOrDiscardEpsilonIncremental symList evalM (ans,rDat) (epsL,etaL)
| isElemEpsilon (epsL,etaL) ans = (ans,rDat)
| otherwise = case newRDat of
Nothing -> (ans,rDat)
Just newRDat' -> (sumAns,newRDat')
where
newAns = getNewAnsEps symList epsL etaL rDat
newRDat = getNewRDatEps evalM newAns rDat
sumAns = addForestsEpsilon ans newAns
mk1stRankDataEtaIncremental :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta,RankDataIncremental,[[Eta]])
mk1stRankDataEtaIncremental symL etaL evalM = output
where
newAns = symAnsatzForestEta symL $ mkForestFromAscList (head etaL,Var 1 1)
newVec = evalAnsatzEtaVecListIncremental evalM newAns
restList = tail etaL
output = case newVec of
Nothing -> if null restList then (EmptyForest,(HM.matrix 0 [], HM.matrix 0 []),[]) else mk1stRankDataEtaIncremental symL restList evalM
Just newVec' -> (newAns, (newMat, newVec'), restList)
where
newVecTrans = HM.tr newVec'
newMat = newVec' Matrix.<> newVecTrans
mk1stRankDataEpsilonIncremental :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon,RankDataIncremental,[(Epsilon,[Eta])])
mk1stRankDataEpsilonIncremental symL epsL evalM = output
where
newAns = symAnsatzForestEps symL $ mkForestFromAscListEpsilon (fst $ head epsL, snd $ head epsL,Var 1 1)
newVec = evalAnsatzEpsilonVecListIncremental evalM newAns
restList = tail epsL
output = case newVec of
Nothing -> if null restList then (M.empty,(HM.matrix 0 [], HM.matrix 0 []),[]) else mk1stRankDataEpsilonIncremental symL restList evalM
Just newVec' -> (newAns,(newMat, newVec'), restList)
where
newVecTrans = HM.tr newVec'
newMat = newVec' Matrix.<> newVecTrans
reduceAnsatzEtaIncremental :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta, HM.Matrix Double)
reduceAnsatzEtaIncremental symL etaL evalM
| null evalM = (EmptyForest, HM.matrix 0 [])
| null etaL = (EmptyForest, HM.matrix 0 [])
| otherwise = (finalForest, finalMat)
where
(ans1,rDat1,restEtaL) = mk1stRankDataEtaIncremental symL etaL evalM
(finalForest, (_,finalMat)) = foldl' (addOrDiscardEtaIncremental symL evalM) (ans1,rDat1) restEtaL
reduceAnsatzEpsilonIncremental :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon, HM.Matrix Double)
reduceAnsatzEpsilonIncremental symL epsL evalM
| null evalM = (M.empty, HM.matrix 0 [])
| null epsL = (M.empty, HM.matrix 0 [])
| otherwise = (finalForest, finalMat)
where
(ans1,rDat1,restEpsL) = mk1stRankDataEpsilonIncremental symL epsL evalM
(finalForest, (_,finalMat)) = foldl' (addOrDiscardEpsilonIncremental symL evalM) (ans1,rDat1) restEpsL
getEtaForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEta, HM.Matrix Double)
getEtaForestIncremental _ _ [] = (EmptyForest, HM.matrix 0 [])
getEtaForestIncremental ord sym evalMs
| null allEtaLists = (EmptyForest, HM.matrix 0 [])
| otherwise = reduceAnsatzEtaIncremental sym allEtaLists evalMs
where
allInds = getEtaInds [1..ord] sym
allEtaLists = map mkEtaList allInds
getEpsForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon, HM.Matrix Double)
getEpsForestIncremental _ _ [] = (M.empty, HM.matrix 0 [])
getEpsForestIncremental ord sym evalMs
| null allEpsLists = (M.empty, HM.matrix 0 [])
| otherwise = reduceAnsatzEpsilonIncremental sym allEpsLists evalMs
where
allInds = getEpsilonInds [1..ord] sym
allEpsLists = map mkEpsilonList allInds
getFullForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> [I.IntMap Int] -> (AnsatzForestEta, AnsatzForestEpsilon, HM.Matrix Double, HM.Matrix Double)
getFullForestIncremental ord sym evalMEta evalMEps = (etaAns, epsAns, etaMat, epsMat)
where
(etaAns,etaMat) = getEtaForestIncremental ord sym evalMEta
(epsAns',epsMat) = getEpsForestIncremental ord sym evalMEps
epsAns = relabelAnsatzForestEpsilon (1 + length (getForestLabels etaAns)) epsAns'
evalToTensSym :: Symmetry -> [(I.IntMap Int, IndTupleST n1 0)] -> [(I.IntMap Int, IndTupleST n1 0)] -> AnsatzForestEta -> AnsatzForestEpsilon -> STTens n1 0 AnsVarR
evalToTensSym (p,ap,b,c,bc) evalEta evalEps ansEta ansEps = symT
where
p' = map (\(x,y) -> (x-1,y-1)) p
ap' = map (\(x,y) -> (x-1,y-1)) ap
b' = map (\(x,y) -> (map (\z -> z-1) x, map (\z' -> z'-1) y) ) b
c' = map (map (subtract 1)) c
bc' = map (map (map (subtract 1))) bc
tens = evalToTens evalEta evalEps ansEta ansEps
symT = foldr cyclicBlockSymATens1 (
foldr cyclicSymATens1 (
foldr symBlockATens1 (
foldr aSymATens1 (
foldr symATens1 tens p'
) ap'
) b'
) c'
) bc'
evalToTens :: [(I.IntMap Int, IndTupleST n1 0)] -> [(I.IntMap Int, IndTupleST n1 0)] -> AnsatzForestEta -> AnsatzForestEpsilon -> STTens n1 0 AnsVarR
evalToTens evalEta evalEps ansEta ansEps = tens
where
etaL = evalAllTensorEta evalEta ansEta
epsL = evalAllTensorEpsilon evalEps ansEps
etaL' = map (\(x,indTuple) -> (indTuple, AnsVar $ I.fromList $ map (\(i,r) -> (i,SField $ fromIntegral r)) x)) etaL
epsL' = map (\(x,indTuple) -> (indTuple, AnsVar $ I.fromList $ map (\(i,r) -> (i,SField $ fromIntegral r)) x)) epsL
etaRmL = filter (\(_,AnsVar b) -> not $ I.null b) etaL'
epsRmL = filter (\(_,AnsVar b) -> not $ I.null b) epsL'
tens = fromListT2 etaRmL &+ fromListT2 epsRmL
evalToTensAbs :: [(I.IntMap Int, Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> [(I.IntMap Int, Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> AnsatzForestEta -> AnsatzForestEpsilon -> ATens n1 0 n2 0 n3 0 AnsVarR
evalToTensAbs evalEta evalEps ansEta ansEps = fromListT6 etaRmL &+ fromListT6 epsRmL
where
etaL = evalAllTensorEtaAbs evalEta ansEta
epsL = evalAllTensorEpsilonAbs evalEps ansEps
etaL' = map (\(x,mult,indTuple) -> (indTuple, AnsVar $ I.fromList $ map (\(i,r) -> (i,fromIntegral $ r*mult)) x)) etaL
epsL' = map (\(x,mult,indTuple) -> (indTuple, AnsVar $ I.fromList $ map (\(i,r) -> (i,fromIntegral $ r*mult)) x)) epsL
etaRmL = filter (\(_,AnsVar b) -> not $ I.null b) $ concatMap (\(x,y) -> zip x (repeat y)) etaL'
epsRmL = filter (\(_,AnsVar b) -> not $ I.null b) $ concatMap (\(x,y) -> zip x (repeat y)) epsL'
mkEvalMap :: Int -> [Int] -> I.IntMap Int
mkEvalMap i = I.fromList . zip [1..i]
mkEvalMaps :: [[Int]] -> [I.IntMap Int]
mkEvalMaps l = let s = length (head l) in map (mkEvalMap s) l
mkEvalMapsInds :: forall (n :: Nat). KnownNat n => [[Int]] -> [(I.IntMap Int, IndTupleST n 0)]
mkEvalMapsInds l = let s = length (head l) in map (\x -> (mkEvalMap s x, (fromListUnsafe $ map toEnum x, Empty))) l
mkAllEvalMaps :: forall (n :: Nat). KnownNat n => Symmetry -> [[Int]] -> ([I.IntMap Int], [I.IntMap Int], [(I.IntMap Int, IndTupleST n 0)], [(I.IntMap Int, IndTupleST n 0)])
mkAllEvalMaps sym l = (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds)
where
evalLEta = filter isEtaList l
evalLEps = filter isEpsilonList l
evalLEtaRed = filter (isLorentzEval sym) evalLEta
evalLEpsRed = filter (isLorentzEval sym) evalLEps
evalMEtaRed = mkEvalMaps evalLEtaRed
evalMEpsRed = mkEvalMaps evalLEpsRed
evalMEtaInds = mkEvalMapsInds evalLEta
evalMEpsInds = mkEvalMapsInds evalLEps
mkAllEvalMapsAbs :: Symmetry -> [([Int], Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> ([I.IntMap Int], [I.IntMap Int], [(I.IntMap Int, Int, [IndTupleAbs n1 0 n2 0 n3 0])], [(I.IntMap Int, Int, [IndTupleAbs n1 0 n2 0 n3 0])])
mkAllEvalMapsAbs sym l = (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds)
where
(headList,_,_) = head l
ord = length headList
evalLEta = filter (\(x,_,_) -> isEtaList x) l
evalLEps = filter (\(x,_,_) -> isEpsilonList x) l
evalLEtaRed = map (\(a,_,_) -> a) $ filter (\(x,_,_) -> isLorentzEval sym x) evalLEta
evalLEpsRed = map (\(a,_,_) -> a) $ filter (\(x,_,_) -> isLorentzEval sym x) evalLEps
evalMEtaRed = mkEvalMaps evalLEtaRed
evalMEpsRed = mkEvalMaps evalLEpsRed
evalMEtaInds = map (\(x,y,z) -> (mkEvalMap ord x, y, z)) evalLEta
evalMEpsInds = map (\(x,y,z) -> (mkEvalMap ord x, y, z)) evalLEps
mkAnsatzTensorIncrementalSym :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorIncrementalSym ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL
(ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTensSym symmetries evalMEtaInds evalMEpsInds ansEta ansEps
mkAnsatzTensorIncremental :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorIncremental ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL
(ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTens evalMEtaInds evalMEpsInds ansEta ansEps
mkAnsatzTensorIncrementalAbs :: Int -> Symmetry -> [([Int], Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> (AnsatzForestEta, AnsatzForestEpsilon, ATens n1 0 n2 0 n3 0 AnsVarR)
mkAnsatzTensorIncrementalAbs ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMapsAbs symmetries evalL
(ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTensAbs evalMEtaInds evalMEpsInds ansEta ansEps
assocsToMat :: [[(Int,Int)]] -> HM.Matrix Double
assocsToMat l = HM.assoc (m,n) 0 l'
where
l' = concat $ zipWith (\r z -> map (\(x,y) -> ((z-1, x-1), fromIntegral y)) r) l [1..]
sparse = M.fromList l'
m = maximum (map (\((x,_),_) -> x) l') + 1
n = maximum (map (\((_,x),_) -> x) l') + 1
getPivots :: [[(Int,Int)]] -> [Int]
getPivots matList = map (1+) pivots
where
mat = assocsToMat matList
pivots = independentColumns mat
reduceLinDepsFastEta :: [I.IntMap Int] -> AnsatzForestEta -> AnsatzForestEta
reduceLinDepsFastEta evalM ansEta = newEtaAns
where
etaL = evalAllEta evalM ansEta
etaVars = getPivots etaL
allEtaVars = getForestLabels ansEta
remVarsEta = allEtaVars \\ etaVars
newEtaAns = relabelAnsatzForest 1 $ removeVarsEta remVarsEta ansEta
reduceLinDepsFastEps :: [I.IntMap Int] -> AnsatzForestEpsilon -> AnsatzForestEpsilon
reduceLinDepsFastEps evalM ansEps = newEpsAns
where
epsL = evalAllEpsilon evalM ansEps
epsVars = getPivots epsL
allEpsVars = getForestLabelsEpsilon ansEps
remVarsEps = allEpsVars \\ epsVars
newEpsAns = relabelAnsatzForestEpsilon 1 $ removeVarsEps remVarsEps ansEps
mkAnsatzFast :: Int -> Symmetry -> [I.IntMap Int] -> [I.IntMap Int] -> (AnsatzForestEta, AnsatzForestEpsilon)
mkAnsatzFast ord symmetries evalMEtaRed evalMEpsRed = (ansEtaRed, ansEpsRed)
where
ansEta = getEtaForestFast ord symmetries
ansEpsilon = getEpsForestFast ord symmetries
ansEtaRed = reduceLinDepsFastEta evalMEtaRed ansEta
ansEpsRed' = reduceLinDepsFastEps evalMEpsRed ansEpsilon
ansEpsRed = relabelAnsatzForestEpsilon (1 + length (getForestLabels ansEtaRed)) ansEpsRed'
mkAnsatzTensorFastSym :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorFastSym ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL
(ansEta, ansEps) = mkAnsatzFast ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTensSym symmetries evalMEtaInds evalMEpsInds ansEta ansEps
mkAnsatzTensorFast :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorFast ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL
(ansEta, ansEps) = mkAnsatzFast ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTens evalMEtaInds evalMEpsInds ansEta ansEps
mkAnsatzTensorFastAbs :: Int -> Symmetry -> [([Int], Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> (AnsatzForestEta, AnsatzForestEpsilon, ATens n1 0 n2 0 n3 0 AnsVarR)
mkAnsatzTensorFastAbs ord symmetries evalL = (ansEta, ansEps, tens)
where
(evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMapsAbs symmetries evalL
(ansEta, ansEps) = mkAnsatzFast ord symmetries evalMEtaRed evalMEpsRed
tens = evalToTensAbs evalMEtaInds evalMEpsInds ansEta ansEps
countEqualInds :: [Int] -> (Int,Int,Int,Int)
countEqualInds [] = (0,0,0,0)
countEqualInds (i:xs)
| i == 0 = (a+1,b,c,d)
| i == 1 = (a,b+1,c,d)
| i == 2 = (a,b,c+1,d)
| i == 3 = (a,b,c,d+1)
| otherwise = error "wrong index"
where
(a,b,c,d) = countEqualInds xs
isEtaList :: [Int] -> Bool
isEtaList l = let (a,b,c,d) = countEqualInds l in even a && even b && even c && even d
isEpsilonList :: [Int] -> Bool
isEpsilonList l = let (a,b,c,d) = countEqualInds l in odd a && odd b && odd c && odd d
filterPSym :: [Int] -> (Int,Int) -> Bool
filterPSym inds (i,j) = (inds !! (i-1)) <= (inds !! (j-1))
filterASym :: [Int] -> (Int,Int) -> Bool
filterASym inds (i,j) = (inds !! (i-1)) < (inds !! (j-1))
getPairs :: [a] -> [(a, a)]
getPairs [a,b] = [(a,b)]
getPairs (x:xs) = (x, head xs) : getPairs xs
getPairs _ = error "invalid index combination"
filterCSym :: [Int] -> [Int] -> Bool
filterCSym inds i = and boolL
where
pairL = getPairs i
boolL = map (filterPSym inds) pairL
filterBSym :: [Int] -> ([Int],[Int]) -> Bool
filterBSym _ ([],[]) = True
filterBSym inds (x:xs,y:ys)
| xVal < yVal = True
| xVal == yVal = filterBSym inds (xs,ys)
| otherwise = False
where
xVal = inds !! (x-1)
yVal = inds !! (y-1)
filterBSym _ _ = error "cannot non-empty list w.r.t. empty symmetries"
filterBCSym :: [Int] -> [[Int]] -> Bool
filterBCSym inds i = and boolL
where
pairL = getPairs i
boolL = map (filterBSym inds) pairL
filterAllSym :: [Int] -> Symmetry -> Bool
filterAllSym inds (p,ap,b,c,bc) = and (p' ++ ap' ++ c' ++ b' ++ bc')
where
p' = map (filterPSym inds) p
ap' = map (filterASym inds) ap
c' = map (filterCSym inds) c
b' = map (filterBSym inds) b
bc' = map (filterBCSym inds) bc
isLorentzEval :: Symmetry -> [Int] -> Bool
isLorentzEval sym inds = inds == canonicalL
where
allInds = filterMins $ getAllIndLists inds
canonicalL = minimum $ map (canonicalizeList sym) allInds
filterMins :: [[Int]] -> [[Int]]
filterMins l = map fst $ filter (\x -> n == snd x) l'
where
l' = map (\x -> (x,sum x)) l
n = minimum $ map snd l'
getAllIndListsMap :: I.IntMap Int -> [I.IntMap Int]
getAllIndListsMap iMap = map (\x -> I.map ((I.!) x) iMap) allSwaps
where
inds = nub $ I.elems iMap
n = length inds
allSwaps = map ((\x y -> I.fromList $ zip x y) inds) $ permutations [0..n-1]
getAllIndLists :: [Int] -> [[Int]]
getAllIndLists l = map I.elems $ getAllIndListsMap $ I.fromList $ zip [1..] l
canonicalizePair :: (Int,Int) -> I.IntMap Int -> I.IntMap Int
canonicalizePair (i,j) iMap
| (I.!) iMap i <= (I.!) iMap j = iMap
| otherwise = I.mapKeys swapKeys iMap
where
swapKeys x
| x == i = j
| x == j = i
| otherwise = x
canonicalizeBlockPair :: ([Int],[Int]) -> I.IntMap Int -> I.IntMap Int
canonicalizeBlockPair ([i],[j]) iMap
| (I.!) iMap i <= (I.!) iMap j = iMap
| otherwise = I.mapKeys swapKeys iMap
where
swapKeys x
| x == i = j
| x == j = i
| otherwise = x
canonicalizeBlockPair (i:is,j:js) iMap
| iVal < jVal = iMap
| iVal > jVal = I.mapKeys (swapBlocks (i:is,j:js)) iMap
| iVal == jVal = newMap
where
iVal = (I.!) iMap i
jVal = (I.!) iMap j
swapBlocks (m1,m2) x = let m = I.fromList $ zip m1 m2 ++ zip m2 m1
in fromMaybe x $ I.lookup x m
newMap = canonicalizeBlockPair (is,js) iMap
canonicalizeBlockPair _ _ = error "invalid index combination"
canonicalizeIntMap :: Symmetry -> I.IntMap Int -> I.IntMap Int
canonicalizeIntMap (p,ap,b,c,bc) iMap = iMap2
where
allBlocks = b ++ concatMap mkBlocksFromBlockCycle bc
allPairs = p ++ ap ++ concatMap mkSymsFromCycle c
iMap1 = foldr canonicalizePair iMap allPairs
iMap2 = foldr canonicalizeBlockPair iMap1 allBlocks
canonicalizeList :: Symmetry -> [Int] -> [Int]
canonicalizeList sym inds = I.elems $ canonicalizeIntMap sym $ I.fromList $ zip [1..] inds
allList' :: Int -> [(Int,Int)] -> [(Int,Int)] -> [(Int,Int)] -> [(Int,Int)] -> [[Int]]
allList' 1 _ _ symBounds aSymBounds = case (symB, aSymB) of
(Just j, Nothing) -> [[k] | k <- [j..3]]
(Nothing, Just j) -> [[k] | k <- [j+1..3]]
(Nothing, Nothing) -> [[0], [1], [2], [3]]
(Just j, Just k) -> [[k'] | k' <- [max j (k+1) .. 3]]
where
(symB,aSymB) = (lookup 1 symBounds, lookup 1 aSymBounds)
allList' i syms aSyms symBounds aSymBounds = concatMap (\x -> (:) <$> [x] <*> allList' (i-1) newSyms newASyms (newSymBounds x) (newASymBounds x)) l
where
(symB,aSymB) = (lookup 1 symBounds, lookup 1 aSymBounds)
l' = case (symB, aSymB) of
(Just j, Nothing) -> [j..3]
(Nothing, Just j) -> [j+1..3]
(Nothing, Nothing) -> [0..3]
(Just j, Just k) -> [max j (k+1) .. 3]
l = if isJust newASymB then filter (<3) l' else l'
newSyms = map (\(x,y) -> (x-1,y-1)) syms
newASyms = map (\(x,y) -> (x-1,y-1)) aSyms
newSymB = lookup 1 syms
newASymB = lookup 1 aSyms
newSymBounds' = map (\(x,y) -> (x-1,y-1)) symBounds
newASymBounds' = map (\(x,y) -> (x-1,y-1)) aSymBounds
newSymBounds x' = case newSymB of
Just j -> (j-1,x') : newSymBounds'
Nothing -> newSymBounds'
newASymBounds x' = case newASymB of
Just j -> (j-1,x') : newASymBounds'
Nothing -> newASymBounds'
allList :: Int -> Symmetry -> [[Int]]
allList ord (syms,aSyms,_,_,_) = allList' ord syms aSyms [] []
mkAnsatzTensorIncrementalSym' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorIncrementalSym' ord symmetries = mkAnsatzTensorIncrementalSym ord symmetries evalL
where
evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries
mkAnsatzTensorFastSym' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorFastSym' ord symmetries = mkAnsatzTensorFastSym ord symmetries evalL
where
evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries
mkAnsatzTensorIncremental' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorIncremental' ord symmetries = mkAnsatzTensorIncremental ord symmetries evalL
where
evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries
mkAnsatzTensorFast' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR)
mkAnsatzTensorFast' ord symmetries = mkAnsatzTensorFast ord symmetries evalL
where
evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries
trianMapArea :: I.IntMap [Int]
trianMapArea = I.fromList $ zip [1..21] list
where
list = [ [a,b,c,d] | a <- [0..2], b <- [a+1..3], c <- [a..2], d <- [c+1..3], isAreaSorted a b c d]
trianMap2 :: I.IntMap [Int]
trianMap2 = I.fromList $ zip [1..10] list
where
list = [ [p,q] | p <- [0..3], q <- [p..3]]
isAreaSorted :: Int -> Int -> Int -> Int -> Bool
isAreaSorted a b c d
| a < c || (a == c && b <= d) = True
| otherwise = False
areaMult :: [Int] -> Int
areaMult [a,b,c,d]
| a == c && b == d = 4
| otherwise = 8
areaMult _ = error "expected four indices"
iMult2 :: [Int] -> Int
iMult2 [p,q] = if p == q then 1 else 2
iMult2 _ = error "expected two indices"
areaList4 :: [([Int], Int, [IndTupleAbs 1 0 0 0 0 0])]
areaList4 = list
where
trianArea = trianMapArea
list = [ let a' = (I.!) trianArea a in (a', areaMult a', [(singletonInd (Ind20 $ a-1), Empty, Empty, Empty, Empty, Empty)]) | a <- [1..21] ]
areaList6 :: [([Int], Int, [IndTupleAbs 1 0 1 0 0 0])]
areaList6 = list
where
trian2 = trianMap2
trianArea = trianMapArea
list = [ let (a',i') = ((I.!) trianArea a, (I.!) trian2 i) in (a' ++ i', areaMult a' * iMult2 i', [(singletonInd (Ind20 $ a-1), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty)]) | a <- [1..21], i <- [1..10]]
areaList8 :: [([Int], Int, [IndTupleAbs 2 0 0 0 0 0])]
areaList8 = list
where
trianArea = trianMapArea
list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ b', areaMult a' * areaMult b', map (\[_a,_b] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21]]
areaList10_1 :: [([Int], Int, [IndTupleAbs 2 0 0 0 2 0])]
areaList10_1 = list
where
trianArea = trianMapArea
list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ p : b' ++ [q], areaMult a' * areaMult b', map (\[[_a,_p],[_b,_q]] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], p <- [0..3], q <- [0..3], not (a==b && p>q)]
areaList10_2 :: [([Int], Int, [IndTupleAbs 2 0 1 0 0 0])]
areaList10_2 = list
where
trian2 = trianMap2
trianArea = trianMapArea
list = [ let (a',b',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trian2 i) in (a' ++ b' ++ i', areaMult a' * areaMult b' * iMult2 i', [ (Append (Ind20 $ a-1) $ singletonInd (Ind20 $ b-1), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty)] ) | a <- [1..21], b <- [1..21], i <- [1..10] ]
areaList12 :: [([Int], Int, [IndTupleAbs 3 0 0 0 0 0])]
areaList12 = list
where
trianArea = trianMapArea
list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ c', areaMult a' * areaMult b' * areaMult c', map (\[_a,_b,_c] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c] )| a <- [1..21], b <- [a..21], c <- [b..21] ]
areaList12_1 :: [([Int], Int, [IndTupleAbs 2 0 2 0 0 0])]
areaList12_1 = list
where
trian2 = trianMap2
trianArea = trianMapArea
list = [ let (a',i',b',j') = ((I.!) trianArea a, (I.!) trian2 i, (I.!) trianArea b, (I.!) trian2 j) in (a' ++ i' ++ b' ++ j' , areaMult a' * areaMult b' * iMult2 i' * iMult2 j', map (\[[_a,_i],[_b,_j]] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Append (Ind9 $ _i-1) $ singletonInd (Ind9 $ _j-1), Empty, Empty, Empty)) $ nub $ permutations [[a,i],[b,j]] ) | a <- [1..21], b <- [a..21], i <- [1..10], j <- [1..10], not (a==b && i>j) ]
areaList14_1 :: [([Int], Int, [IndTupleAbs 3 0 0 0 2 0])]
areaList14_1 = list
where
trianArea = trianMapArea
list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ p : c' ++ [q], areaMult a' * areaMult b' * areaMult c', map (\[[_b,_p],[_c,_q]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[b,p],[c,q]]) | a <- [1..21], b <- [1..21], c <- [b..21], p <- [0..3], q <- [0..3], not (b==c && p>q) ]
areaList14_2 :: [([Int], Int, [IndTupleAbs 3 0 1 0 0 0])]
areaList14_2 = list
where
trian2 = trianMap2
trianArea = trianMapArea
list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i) in ( a' ++ b' ++ c' ++ i', areaMult a' * areaMult b' * areaMult c' * iMult2 i', map (\[_a,_b] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ c-1), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21], c <- [1..21], i <- [1..10] ]
areaList16_1 :: [([Int], Int, [IndTupleAbs 3 0 1 0 2 0])]
areaList16_1 = list
where
trian2 = trianMap2
trianArea = trianMapArea
list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!