{-|
    Module      :  Data.Number.ER.Real.Arithmetic.Integration
    Description :  simple integration methods
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

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

    Simple integration methods for Haskell functions operating 
    on real number approximations.
-}
module Data.Number.ER.Real.Arithmetic.Integration
(
    integrateCont,
--    integrateDiff,
    integrateCont_R,
    integrateContAdapt_R
)
where

import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes
import Data.Number.ER.Real.Approx.Sequence
import Data.Number.ER.Real.Arithmetic.Elementary

testIntegr1 :: 
    (RA.ERIntApprox ira, Ord ira) => 
    (ConvergRealSeq ira)
testIntegr1 = integrateCont erExp_IR 0 1

integrateCont :: 
    (RA.ERIntApprox ira) => 
    (EffortIndex -> ira -> ira) ->
    (ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira)

integrateCont f = convertBinFuncRA2Seq $ integrateContAdapt_R f

integrateDiff :: 
    (RA.ERIntApprox ira) => 
    (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) ->
    (ConvergRealSeq ra) -> (ConvergRealSeq ra) -> (ConvergRealSeq ra)

integrateDiff f = convertBinFuncRA2Seq $ integrateDiffAdapt_RA f


{-|
    naive integration, using a partition of 2 * prec equally sized intervals
-}
integrateCont_R ::
    (RA.ERIntApprox ira) => 
    (EffortIndex -> ira -> ira) ->
    EffortIndex -> (ira) -> (ira) -> (ira)
integrateCont_R f ix a b =
    sum $ map rectArea rectangles
    where
    rectArea (width, height) = width * height
    rectangles = 
        zip (repeat width) $ map (f ix) covering
    width = (b - a) / (fromInteger rectCount)
    rectCount = 2 * (fromInteger $ toInteger gran)
    gran = effIx2gran ix
    covering = getCoveringIntervals division
    getCoveringIntervals ( pt1 : pt2 : rest ) =
        ((pt1) RA.\/ (pt2)) : (getCoveringIntervals $ pt2 : rest)
    getCoveringIntervals _ = []
    division = map getEndPoint $ [0..rectCount]
    getEndPoint n =
        a + ((fromInteger n) * width)

{-|
    integration using divide and conquer adaptive partitioning
-}
integrateContAdapt_R ::
    (RA.ERIntApprox ira) => 
    (EffortIndex -> ira -> ira) ->
    EffortIndex -> (ira) -> (ira) -> (ira)
integrateContAdapt_R f ix a b =
    sum rectangleAreas
    where
    rectangleAreas = 
        getRs 10 a b
    getRs subix l r
        | RA.getPrecision area >= prec = [area]
        | otherwise =
            (getRs nsubix l m) ++ (getRs nsubix m r)
        where
        prec = foldl1 min [effIx2prec ix, RA.getPrecision l, RA.getPrecision r]
        area = (r - l) * (f subix (l RA.\/ r))
        nsubix = subix + 1
        m = (l + r)/2
        

{-|
    integration using divide and conquer adaptive partitioning
    making use of the derivative of the integrated function
-}
integrateDiffAdapt_RA ::
    (RA.ERIntApprox ira) => 
    (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) ->
    EffortIndex -> (ra) -> (ra) -> (ra)
integrateDiffAdapt_RA f prec a b =
    error "TODO"
    
{-
    sum rectangleAreas
    where
    rectangleAreas = 
        getRs prec a b
    getRs p l r
        | getPrecision area >= prec = [area]
        | otherwise =
            (getRs np l m) ++ (getRs np m r)
        where
        np = p + 1
        m = (l + r)/2
--        area = areaDiff
        area = areaRect /\ areaDiff
            -- merge the information given by the rectangle method
            -- with the information given by the derivative method
        areaRect = w * fVal -- same as in integrateContAdapt_R
        (fVal, fDeriv) = applyRdiffR f p (l \/ r)
        w = r - l
        areaDiff
            | isExact fDeriv = w * (fl + fr) / 2 -- derivative is constant and perfectly known
            | otherwise = areaLow \/ areaHigh
        fl = fst $ applyRdiffR f (2 * p) l
        fr = fst $ applyRdiffR f (2 * p) r
            -- interestingly, we have to request fl, fr with higher precision than
            -- we requested fDeriv so that the derivative would be of any use
            -- with these values - replace (2 * p) by p and it will not converge!
        -- area computed by a scary formula:
        areaLow = t + w * (fl * dHigh - fr * dLow) / dDiff
        areaHigh = - t - w * (fl * dLow - fr * dHigh) / dDiff -- swap dHigh and dLow
        t = (w^2 * dLow * dHigh + (fr - fl)^2)/(2 * dDiff)
        dDiff = dHigh - dLow
        (dLow, dHigh) = getBounds fDeriv
-}