module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
where
import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.Misc
import qualified Data.Map as Map
import Data.Typeable
import Data.Generics.Basics
import Data.Binary
errorModule msg = error $ "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom: " ++ msg
data ERChebPoly box b =
ERChebPoly
{
chplCoeffs :: (Map.Map (TermKey box) b)
}
deriving (Eq, Typeable, Data)
instance (Ord a, Binary a, Binary b) => Binary (ERChebPoly a b) where
put (ERChebPoly a) = put a
get = get >>= \a -> return (ERChebPoly a)
chplCheck prgLocation p@(ERChebPoly coeffs)
| ok = p
| otherwise =
unsafePrint (prgLocation ++ " problem with p = \n" ++ show p) p
where
ok =
Map.fold (&&) True $ Map.map coeffOK coeffs
coeffOK c =
not $ B.isERNaN c
type TermKey box = box
chplConstTermKey :: (DomainBox box varid d) => box
chplConstTermKey = DBox.noinfo
chplIsConstTermKey :: (DomainBox box varid d) => box -> Bool
chplIsConstTermKey = DBox.isNoinfo
chplTermOrder :: (DomainBox box varid d, Num d) => box -> d
chplTermOrder termKey = DBox.fold (+) 0 termKey
chplTermArity :: (DomainBox box varid d) => box -> Int
chplTermArity termKey = length $ DBox.keys termKey
chplGetDegree ::
(B.ERRealBase b, DomainBox box varid d, Num d, Ord d) =>
(ERChebPoly box b) ->
d
chplGetDegree (ERChebPoly coeffs) =
foldl max 0 $ map chplTermOrder $ Map.keys coeffs
chplGetVars (ERChebPoly coeffs) =
DBox.keys $ foldl DBox.union DBox.noinfo $ Map.keys coeffs
chplGetGranularity (ERChebPoly coeffs) =
foldl max 0 $ map B.getGranularity $ Map.elems coeffs
chplSetMinGranularity gran (ERChebPoly coeffs) =
ERChebPoly $ Map.map (B.setMinGranularity gran) coeffs
chplSetGranularity gran (ERChebPoly coeffs) =
ERChebPoly $ Map.map (B.setGranularity gran) coeffs
chplConst ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
b ->
ERChebPoly box b
chplConst val =
(ERChebPoly $ Map.singleton chplConstTermKey val)
chplVar ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
varid ->
ERChebPoly box b
chplVar varName =
ERChebPoly $ Map.singleton (DBox.singleton varName 1) 1
--{-|
---}
chplNormaliseDown, chplNormaliseUp ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
ERChebPoly box b -> ERChebPoly box b
chplNormaliseUp = snd . chplNormalise
chplNormaliseDown = fst . chplNormalise
chplNormalise ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)
chplNormalise (ERChebPoly coeffs) =
(ERChebPoly $ coeffsNo0T0Down,
ERChebPoly $ coeffsNo0T0Up)
where
coeffsNo0T0Down =
Map.insertWith plusDown chplConstTermKey err coeffsNo0T0
coeffsNo0T0Up =
Map.insertWith plusUp chplConstTermKey err coeffsNo0T0
(coeffsNo0T0, err) =
foldl addTermNo0T0 (Map.empty, 0) $ Map.toList coeffs
addTermNo0T0 (prevCoeffs, prevErr) (term, coeff)
| coeff == 0 =
(prevCoeffs, prevErr)
| otherwise =
(newCoeffs, newErr)
where
newTerm =
DBox.filter (> 0) term
newCoeffs =
Map.insert newTerm newCoeffUp prevCoeffs
newCoeffUp = prevCoeff + coeff
newCoeffDown = prevCoeff `plusDown` coeff
prevCoeff =
Map.findWithDefault 0 newTerm prevCoeffs
newErr = newCoeffUp newCoeffDown
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Show (ERChebPoly box b)
where
show = chplShow False
chplShow ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Bool ->
ERChebPoly box b ->
String
chplShow showChebyshevBasis (ERChebPoly coeffs)
| showChebyshevBasis = "\n" ++ inChebBasis ++ " = \n" ++ inXBasis
| otherwise = inXBasis
where
inChebBasis =
showCoeffs showTermT $ coeffs
inXBasis =
showCoeffs showTermX $ chebToXBasis coeffs
showCoeffs showTerm coeffs =
concatWith " + " $ map showTerm $ Map.toAscList coeffs
showTermT (term, coeff)
| chplIsConstTermKey term = show coeff
| otherwise =
show coeff ++ "*" ++ (concat $ map showT $ DBox.toList term)
showTermX (term, coeff)
| chplIsConstTermKey term = showC coeff
| otherwise =
showC coeff ++ "*" ++ (concat $ map showX $ DBox.toList term)
showT (var, deg) = "T" ++ show deg ++ "(" ++ showVar var ++ ")"
showX (var, deg) = showVar var ++ "^" ++ show deg
showC = B.showDiGrCmp 8 False False
chebToXBasis ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
(Map.Map (TermKey box) b) ->
(Map.Map (TermKey box) b)
chebToXBasis coeffs =
Map.foldWithKey addTerm Map.empty coeffs
where
addTerm term coeff prevXCoeffs =
Map.unionWith (+) prevXCoeffs $
Map.map (\ c -> coeff * (fromInteger c)) $
termXterms term
termXterms ::
(DomainBox box varid Int, Ord box) =>
TermKey box
->
Map.Map (TermKey box) Integer
termXterms term =
foldl addCombination Map.empty $
allCombinations $
map (mapSnd $ \ deg -> chebyXCoeffsLists !! deg) $
DBox.toList term
where
addCombination prevMap (varPowerCoeffTriples) =
Map.insertWith (+) term coeffProduct prevMap
where
term =
DBox.fromList $
filter (\(v,p) -> p > 0) $
map (\(v,(p,_)) -> (v,p)) varPowerCoeffTriples
coeffProduct =
fromInteger $
product $
map (\(_,(_,c)) -> c) varPowerCoeffTriples
chebyXCoeffsLists ::
(Num d1, Enum d1, Num d2, Enum d2) =>
[[(d1, d2)]]
chebyXCoeffsLists =
map convertCoeffs chebyXCoeffs
where
convertCoeffs coeffs =
filter ((/= 0) . snd) $ zip [0,1..] coeffs
chebyXCoeffs ::
(Num d, Enum d) =>
[[d]]
chebyXCoeffs =
aux
[1]
[0,1]
where
aux tnM2 tnM1 =
tnM2 : (aux tnM1 (newTerm tnM2 tnM1))
newTerm tnM2 tnM1 =
zipWith () (0 : (map (*2) tnM1)) (tnM2 ++ [0,0..])