```{-# LANGUAGE FlexibleContexts #-}
{-|
Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
Description :  (internal) evaluation of polynomials
Copyright   :  (c) 2007-2008 Michal Konecny

Maintainer  :  mik@konecny.aow.cz
Stability   :  experimental
Portability :  portable

Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".

Implementation of various evaluation functions related to polynomials.
-}
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

{-|
Evaluate a polynomial at a point, consistently rounding upwards and downwards.
-}
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

{-|
Evaluate a polynomial at a real number approximation
-}
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

{-|
Substitute several variables in a polynomial with real number approximations,
rounding downwards and upwards.
-}
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

{-|
Compose two polynomials, rounding upwards
provided the second polynomial maps [-1,1] into [-1,1].
-}
chplCompose ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
ERChebPoly box b ->
Map.Map varid (ERChebPoly box b)
{-^ variable to substitute, polynomial to substitute  -} ->
(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
```