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

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

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

Implementation of safely rounded integration of polynomials
and other related functions.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration
where

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

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, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.Real.Approx.Interval
import Data.Number.ER.Misc

import qualified Data.Map as Map

{-|
Approximate from below and from above the integral of a polynomial.

Based on the following formulas for Chebyshev polynomials:

>     \int T_n(x)dx =
>        T_{n+1}(x)/2(n+1) - T_{n-1}(x)/2(n-1)

>     \int T_1(x)dx =
>        T_2(x)/4 + 1/4

>     \int T_0(x)dx =
>        T_1(x)

-}
chplIntegrate ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
varid {-^ variable to integrate by -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplIntegrate x p@(ERChebPoly coeffs) =
--    unsafePrintReturn
--    (
--        "ERChebPoly: integrate:"
--        ++ "\n p = " ++ show p
--        ++ "\n result = "
--    )
(pNp1Down -. pNm1Up,
pNp1Up -^ pNm1Down)
--    (chplRemoveZeroTermsDown $pNp1Down - pNm1Up, -- chplRemoveZeroTermsUp$ pNp1Up - pNm1Down)
where
pNp1Up =
ERChebPoly $Map.insertWith plusUp chplConstTermKey errBoundNp1$
Map.fromList coeffsNp1
pNp1Down =
ERChebPoly $Map.insertWith plusDown chplConstTermKey (- errBoundNp1)$
Map.fromList coeffsNp1
pNm1Up =
ERChebPoly $Map.insertWith plusUp chplConstTermKey errBoundNm1$
Map.fromList coeffsNm1
pNm1Down =
ERChebPoly $Map.insertWith plusDown chplConstTermKey (- errBoundNm1)$
Map.fromList coeffsNm1
(coeffsNp1, errBoundNp1) =
foldl cfNp1 ([],0) coeffsList
(coeffsNm1, errBoundNm1) =
foldl cfNm1 ([],0) coeffsList
coeffsList = Map.toList coeffs
cfNp1 (prevTerms, prevErr) (termKey, coeff)
| n == 0 =
((termKeyNp1, coeff):prevTerms, prevErr)
| n == 1 =
((termKeyN0, coeff0Up):(termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeff0Err + coeffNp1Err)
| otherwise =
((termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeffNp1Err)
where
termKeyNp1 = DBox.insert x (n + 1) termKey
termKeyNm1 = DBox.insert x (n - 1) termKey
termKeyN0 = DBox.delete x termKey
n = DBox.findWithDefault 0 x termKey
coeffNp1Err = coeffNp1Up - coeffNp1Down
coeffNp1Up = coeff / (2*nB + 2)
coeffNp1Down = -((-coeff) / (2*nB + 2))
nB = fromInteger $toInteger n coeff0Up = coeff / 4 coeff0Down = - ((- coeff) / 4) coeff0Err = coeff0Up - coeff0Down cfNm1 (prevTerms, prevErr) (termKey, coeff) | n == 0 || n == 1 = (prevTerms, prevErr) | otherwise = ((termKeyNm1, coeffNm1Up):prevTerms, prevErr + coeffNm1Err) where termKeyNm1 = DBox.insert x (n - 1) termKey n = DBox.findWithDefault 0 x termKey coeffNm1Up = coeff / (2*nB - 2) coeffNm1Down = -((-coeff) / (2*nB - 2)) nB = fromInteger$ toInteger n
coeffNm1Err = coeffNm1Up - coeffNm1Down

--{-|
--    measure the volume between a polynomial and the zero axis on [-1,1]^n
---}
--chplVolumeAboveZero ::
--    (B.ERRealBase b, DomainBox box varid Int, Ord box,
--     DomainBoxMappable boxb boxbb varid b [ERInterval b]) =>
--    [varid] ->
--    ERChebPoly box b ->
--    (b,b)
--chplVolumeAboveZero vars p@(ERChebPoly coeffs) =
----    unsafePrint ("chplVolumeAboveZero: returning:" ++ show result) $---- unsafePrint ("chplVolumeAboveZero: vars = " ++ show vars)$
--    result
--    where
--    result =
--        (- (integUpAtOddCorners - integDownAtEvenCorners), integUpAtEvenCorners - integDownAtOddCorners)
--    integUpAtEvenCorners = sumUp $map (chplEvalUp integUp) evenCorners -- integUpAtOddCorners = sumUp$ map (chplEvalUp integUp) oddCorners
--    integDownAtEvenCorners = sumDown $map (chplEvalDown integDown) evenCorners -- integDownAtOddCorners = sumDown$ map (chplEvalDown integDown) oddCorners
--    evenCorners = map (DBox.fromList) evenCornersL
--    oddCorners = map (DBox.fromList) oddCornersL
--    (evenCornersL, oddCornersL) =
--        allPairsCombinationsEvenOdd $zip vars$ repeat (1,-1)
--    integUp = integrateByAllVars snd p vars
--    integDown = integrateByAllVars fst p vars
--    integrateByAllVars pick p [] = p
--    integrateByAllVars pick p (x : xs) =
--        integrateByAllVars pick ip xs
--        where
--        ip = pick \$ chplIntegrate x p
----    vars = chplGetVars p