module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
where
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
import qualified Data.Number.ER.Real.Approx as RA
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, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.Misc
import qualified Data.Map as Map
chplEval ::
(B.ERRealBase b, DomainBox box varid Int, Ord box,
DomainBoxMappable boxb boxbb varid b [(b,b)]) =>
boxb ->
ERChebPoly box b ->
(b, b)
chplEval vals (ERChebPoly coeffs) =
(foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)
where
(termValsLo, termValsHi) =
unzip $ map evalTerm $ Map.toList coeffs
evalTerm (term, c) =
(foldl timesDown c valsLo, foldl timesUp c valsHi)
where
(valsLo, valsHi) =
unzip $ map evalVar $ DBox.toList term
evalVar (varID, degree) =
(DBox.lookup "ERChebPoly.Eval: chplEval" varID valsDegrees) !! degree
valsDegrees =
DBox.map chebyEvalTsRoundDownUp vals
chplEvalDown, chplEvalUp ::
(B.ERRealBase b, DomainBox box varid Int, Ord box,
DomainBoxMappable boxb boxbb varid b [(b,b)]) =>
boxb ->
ERChebPoly box b ->
b
chplEvalUp pt = snd . chplEval pt
chplEvalDown pt = fst . chplEval pt
chebyEvalTsRoundDownUp ::
(Num v) =>
v -> [(v,v)]
chebyEvalTsRoundDownUp val =
chebyIterate (1,1) (val, val)
where
chebyIterate tNm2@(tNm2Down, tNm2Up) tNm1@(tNm1Down, tNm1Up) =
tNm2 : (chebyIterate tNm1 (tNDown, tNUp))
where
tNUp = 2 * val * tNm1Up tNm2Down
tNDown = ((2 * val) `timesDown` tNm1Down) tNm2Up
chebyEvalTsExact ::
(Num v) =>
v -> [v]
chebyEvalTsExact val =
chebyIterate 1 val
where
chebyIterate tNm2 tNm1 =
tNm2 : (chebyIterate tNm1 tN)
where
tN = 2 * val * tNm1 tNm2
chplEvalApprox ::
(B.ERRealBase b, RA.ERApprox ra,
DomainBox box varid Int, Ord box,
DomainBoxMappable boxra boxras varid ra [ra],
DomainIntBox boxra varid ra) =>
(b -> ra) ->
boxra ->
ERChebPoly box b ->
ra
chplEvalApprox b2ra vals (ERChebPoly coeffs) =
sum $ map evalTerm $ Map.toList coeffs
where
evalTerm (term, c) =
(b2ra c) * (product $ map evalVar $ DBox.toList term)
evalVar (varID, degree) =
(DBox.lookup "ERChebPoly.Eval: chplEvalApprox: " varID valsDegrees) !! degree
valsDegrees =
DBox.map chebyEvalTsExact vals
chplPartialEvalApprox ::
(B.ERRealBase b, RA.ERApprox ra,
DomainBox box varid Int, Ord box,
DomainBoxMappable boxra boxras varid ra [ra],
DomainIntBox boxra varid ra) =>
(ra -> (b,b)) ->
boxra ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplPartialEvalApprox ra2endpts substitutions (ERChebPoly coeffs) =
(ERChebPoly $ Map.insertWith plusDown chplConstTermKey ( corr) coeffsSubstDown,
ERChebPoly $ Map.insertWith plusUp chplConstTermKey corr coeffsSubstUp)
where
(coeffsSubstDown, coeffsSubstUp, corr) =
Map.foldWithKey processTerm (Map.empty, Map.empty, 0) coeffs
processTerm termKey coeff (coeffsSubstDownPrev, coeffsSubstUpPrev, corrPrev) =
(Map.insertWith plusDown newTermKey newCoeffDown coeffsSubstDownPrev,
Map.insertWith plusUp newTermKey newCoeffUp coeffsSubstUpPrev,
corrPrev + corrVars)
where
corrVars = (substValHi substValLo) * coeff
(newCoeffDown, newCoeffUp)
| coeff > 0 = (coeff `timesDown` substValLo, coeff `timesUp` substValHi)
| coeff < 0 = (coeff `timesDown` substValHi, coeff `timesUp` substValLo)
| otherwise = (0,0)
(substValLo, substValHi) = ra2endpts substVal
(substVal, newTermKey) =
DBox.foldWithKey processVar (1, DBox.noinfo) termKey
processVar varID degree (substValPrev, newTermKeyPrev) =
case DBox.member varID substitutions of
True ->
(substValPrev * (evalVar varID degree), newTermKeyPrev)
False ->
(substValPrev, DBox.insert varID degree newTermKeyPrev)
evalVar varID degree =
(DBox.lookup "ERChebPoly.Eval: chplPartialEvalApprox: " varID valsDegrees) !! degree
valsDegrees =
DBox.map chebyEvalTsExact substitutions
chplCompose ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
ERChebPoly box b ->
Map.Map varid (ERChebPoly box b)
->
(ERChebPoly box b, ERChebPoly box b)
chplCompose maxDegree p@(ERChebPoly coeffs) substitutions =
(foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)
where
(termValsLo, termValsHi) =
unzip $ map evalTerm $ Map.toList coeffs
evalTerm (term, c) =
(foldl timesDown cPoly valsLo, foldl timesUp cPoly valsHi)
where
cPoly = chplConst c
(valsLo, valsHi) =
unzip $ map evalVar $ DBox.toList term
evalVar (varID, degree) =
case Map.lookup varID substDegrees of
Nothing ->
(varPoly, varPoly)
Just pvDegrees ->
pvDegrees !! degree
where
varPoly =
ERChebPoly $ Map.singleton (DBox.singleton varID degree) 1
substDegrees =
Map.map mkPVDegrees substitutions
mkPVDegrees pv =
map
(mapPair
(chplReduceDegreeDown maxDegree,
chplReduceDegreeUp maxDegree)) $
chebyEvalTsRoundDownUp pv