Integration.hs

{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration
    Description :  (internal) integration of polynomials etc  
    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 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.Field
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds

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

{-|
     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 (ERChebPoly coeffs) =
--    unsafePrint
--    (
--        "ERChebPoly: integrate:"
--        ++ "\n pNp1Down = " ++ chplShow True pNp1Down 
--        ++ "\n pNm1Up = " ++ chplShow True pNm1Up 
--    )
    (chplNormaliseDown $ pNp1Down - pNm1Up, 
     chplNormaliseUp $ 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 =
            ((termKeyNm1, 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 
        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 =
            ((chplConstTermKey, 0):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 [(b,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 (\pt -> chplEvalUp pt integUp) evenCorners
    integUpAtOddCorners = sumUp $ map (\pt -> chplEvalUp pt integUp) oddCorners 
    integDownAtEvenCorners = sumDown $ map (\pt -> chplEvalDown pt integDown) evenCorners  
    integDownAtOddCorners = sumDown $ map (\pt -> chplEvalDown pt 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
      
    
--{-|
--    Calculate approximations to the Chebyshev nodes.
---}
--chebNodes ::
--    (B.ERRealBase b) =>
--    Granularity ->
--    [[b]] -- ^ ith element is the ordered list of ith order Chebyshev nodes  
--chebNodes gran =
--    error "ERChebPoly: chebNodes: not implemented yet"
    
    
{-|
    Differentiate a polynomial using one of its variables. 
-}
chplDifferentiate ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    ERChebPoly box b ->
    varid {-^ variable to differentiate over -} ->
    ERChebPoly box b
chplDifferentiate (ERChebPoly coeffs) varName =
    errorModule "chplDifferentiate: not implemented yet"