module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
where
import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.BasicTypes.DomainBox as DBox
import Data.Number.ER.BasicTypes.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, Ord, Typeable, Data)
type TermKey box = box
instance (Ord a, Binary a, Binary b) => Binary (ERChebPoly a b) where
put (ERChebPoly a) = put a
get = get >>= \a -> return (ERChebPoly a)
chplHasNoNaN p@(ERChebPoly coeffs) =
Map.fold (&&) True $ Map.map coeffOK coeffs
where
coeffOK c =
not $ B.isERNaN c
chplHasNoNaNOrInfty p@(ERChebPoly coeffs) =
Map.fold (&&) True $ Map.map coeffOK coeffs
where
coeffOK = isBounded
isBounded c
| B.isERNaN c = False
| B.isPlusInfinity c = False
| B.isPlusInfinity ( c) = False
| otherwise = True
chplCheck prgLocation p
| chplHasNoNaNOrInfty p = p
| otherwise =
unsafePrint (prgLocation ++ " problem with p = \n" ++ show p) p
chplCompareApprox (ERChebPoly coeffs1) (ERChebPoly coeffs2) =
compare coeffs1 coeffs2
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
chplGetConst ::
(B.ERRealBase b, DomainBox box varid d, Num d, Ord d) =>
(ERChebPoly box b) ->
Maybe b
chplGetConst (ERChebPoly coeffs) =
case Map.toList coeffs of
[] -> Just 0
[(key,val)] | chplIsConstTermKey key -> Just val
_ -> Nothing
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
chplAffine ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
b ->
Map.Map varid b ->
ERChebPoly box b
chplAffine at0 varCoeffs =
ERChebPoly $
Map.insert chplConstTermKey at0 $
Map.mapKeys (\ i -> DBox.singleton i 1) varCoeffs
chplRemoveZeroTermsUp ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
ERChebPoly box b -> ERChebPoly box b
chplRemoveZeroTermsUp (ERChebPoly coeffs) =
ERChebPoly coeffsNo0T0Up
where
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 = prevErr + newCoeffUp newCoeffDown
chplCountTerms ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
ERChebPoly box b -> Int
chplCountTerms (ERChebPoly coeffs) =
Map.size coeffs
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Show (ERChebPoly box b)
where
show = chplShow 8 False False
chplShow ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
Bool ->
Bool ->
ERChebPoly box b ->
String
chplShow digitsToShow showGranularity showChebyshevBasis (ERChebPoly coeffs)
| showChebyshevBasis = "\n" ++ inChebBasis ++ " = \n" ++ inXBasis
| otherwise = inXBasis
where
inChebBasis =
showCoeffs showTermT $ Map.filter (\c -> c /= 0) $ 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 digitsToShow showGranularity False
chebToXBasis ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
(Map.Map (TermKey box) b) ->
(Map.Map (TermKey box) b)
chebToXBasis coeffs =
Map.filter (\c -> c /= 0) $
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..])