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

    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