{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose
    Description :  (internal) composition 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 pointwise consistently rounded polynomial composition.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose
where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure

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

{-|
    Compose a polynomial and an enclosure, producing a correcly rounded enclosure,
    assuming the second polynomial maps [-1,1] into [-1,1].
-}
enclCompose ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    Int {-^ max degree for result -} -> 
    Int {-^ max approx size for result -} ->
    ERChebPoly box b {-^ @f@ -} ->
    varid {-^ variable @v@ to substitute in @f@ -} -> 
    (ERChebPoly box b, ERChebPoly box b)
         {-^ enclosure of a function @f_v@ to substitute for @v@ 
             that maps @[-1,1]@ into @[-1,1]@  -} ->
    (ERChebPoly box b, ERChebPoly box b)
        {-^ lower bound and upper bound -}


enclCompose maxDegree maxSize p@(ERChebPoly coeffs) substVar substEncl =
    result
{------------------------------
 The algorithm: separate from the polynomial 
 all terms for each degree of the substituted variable,
 giving rise to a number of polynomials.
 These polynomials are then used as coefficients multiplying
 the enclosure evaluations of the Chebyshev polynomials 
 over the substituted enclosure.
-------------------------------}
    where
    result =
        Map.fold (+:) (enclConst 0) $ Map.mapWithKey evalDegree degreePolynomialMap
    degreePolynomialMap =
        Map.foldWithKey extractTerm Map.empty coeffs
    extractTerm term c prevPolynomMap =
        Map.insertWith Map.union substVarDegree (Map.singleton termNoSubstVar c) prevPolynomMap
        where
        substVarDegree = DBox.findWithDefault 0 substVar term
        termNoSubstVar = DBox.delete substVar term
    evalDegree degree degreeCoeffs =
        enclMultiply maxDegree maxSize (substPolyDegrees !! degree) (chplNeg degreePoly, degreePoly)
        where
        degreePoly = ERChebPoly degreeCoeffs
    substPolyDegrees =
        enclEvalTs maxSize maxDegree substEncl

{------------------------------
 The following algorithm is quite wasteful when the polynomial
 contains other variables besides the one being substituted.
-------------------------------}
--chplComposeWithEncl maxDegree maxSize p@(ERChebPoly coeffs) substVar substEncl =
--    result
--    where
--    result =
--        foldl (+:) (enclConst 0) $ map evalTerm $ Map.toList coeffs
--    evalTerm (term, c) =
--        enclScale c $ 
--            foldl (enclMultiply maxDegree maxSize) (enclConst 1) $ 
--                map evalVar $ DBox.toList term
--    evalVar (var, degree) =
--        case var == substVar of
--            True ->
--                substPolyDegrees !! degree
--            False ->
--                (chplNeg varPoly, varPoly)
--        where
--        varPoly = 
--            ERChebPoly $ Map.singleton (DBox.singleton var degree) 1
--    substPolyDegrees =
--        enclEvalTs maxSize maxDegree substEncl

        

{-|
    Compose two polynomials, rounding upwards
    provided the second polynomial maps [-1,1] into [-1,1].
-}
enclComposeMany ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    Int {-^ max degree for result -} -> 
    Int {-^ max approx size for result -} ->
    ERChebPoly box b ->
    Map.Map varid (ERChebPoly box b, ERChebPoly box b) 
     {-^ variables to substitute and the enclosures to substitute for each of them respectively  -} ->
    (ERChebPoly box b, ERChebPoly box b)
        {-^ lower bound (negated) and upper bound -}
enclComposeMany maxDegree maxSize p@(ERChebPoly coeffs) substitutions =
    result
    where
    result =
        foldl (+:) (enclConst 0) $ map evalTerm $ Map.toList coeffs
    evalTerm (term, c) =
        enclScale maxDegree maxSize c $ 
            foldl (enclMultiply maxDegree maxSize) (enclConst 1) $ 
                map evalVar $ DBox.toList term
    evalVar (varID, degree) =
        case Map.lookup varID substDegrees of
            Nothing ->
                (chplNeg varPoly, varPoly)
            Just pvDegrees ->
                pvDegrees !! degree
        where
        varPoly = 
            ERChebPoly $ Map.singleton (DBox.singleton varID degree) 1
    substDegrees =
        Map.map mkPVDegrees substitutions
    mkPVDegrees pvEncl =
        enclEvalTs maxSize maxDegree pvEncl