{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
    Description :  (internal) bounds of single and multiple 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 functions related to the bounds of polynomials.    
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds 
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.Reduce
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval

import qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Base as B
import Data.Number.ER.Real.Approx.Interval
import Data.Number.ER.Real.Arithmetic.LinearSolver
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc

import qualified Data.Map as Map

import Data.List

{-|
    Find an upper bound on a polynomial over the 
    unit domain [-1,1]^n.  
-}
chplUpperBound ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplUpperBound ix p = snd $ chplBounds ix p

{-|
    Find a lower bound on a polynomial over the 
    unit domain [-1,1]^n.  
-}
chplLowerBound ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplLowerBound ix p = fst $ chplBounds ix p

{-|
    Find both lower and upper bounds on a polynomial over the 
    unit domain [-1,1]^n.  
-}
chplBounds ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    (b,b)
chplBounds = chplBoundsAffine

{-|
    Find bounds on a polynomial over the unit domain [-1,1]^n.
    
    Fast but inaccurate method, in essence
    taking the maximum of the upper affine reduction.
-}
chplBoundsAffine ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    (b,b)
chplBoundsAffine ix p@(ERChebPoly coeffs) =
--    unsafePrintReturn
--    (
--        "chplBoundsAffine:"
--        ++ "\n p = " ++ show p
--        ++ "\n noConstCoeffAbsSum = " ++ show noConstCoeffAbsSum
--        ++ "\n result = "
--    )    
    result
    where
    result =
        (constTerm `plusDown` (- noConstCoeffAbsSum),
         constTerm `plusUp` noConstCoeffAbsSum)
    noConstCoeffAbsSum = Map.fold plusUp 0 absCoeffs
    absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs
    constTerm = Map.findWithDefault 0 chplConstTermKey coeffs

{-|
    Find a close upper bound on a quadratic polynomial over the 
    unit domain [-1,1]^n.  

    Much slower and somewhat more accurate method, in essence
    taking the maximum of the upper quadratic reduction.
    
    !!! Not yet properly tested !!!
-}
chplUpperBoundQuadr ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box,
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b],
     DomainBoxMappable boxra boxra varid (ERInterval b) (ERInterval b), 
     DomainIntBox boxra varid (ERInterval b), Num varid, Enum varid) => 
    EffortIndex {-^ how hard to try looking for peaks -} ->
    ERChebPoly box b ->
    b
chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
    quadBound (coeffsQ, vars)
    where
    pQ@(ERChebPoly coeffsQ) = chplReduceDegreeUp 2 p
    vars = chplGetVars pQ
    quadBound (coeffs, vars)
        | null vars =
            Map.findWithDefault 0 chplConstTermKey coeffs
        | hasInteriorPeak =
            foldl max peakValue edgeBounds
        | otherwise =
            foldl1 max edgeBounds
        where
        edgeBounds =
            map quadBound $ concat $ map removeVar vars
        (hasInteriorPeak, peakValue) =
            case maybePeak of
                Just peak ->
                    (noPositiveSquare -- if any term x^2 has a positive coeff, there is no peak  
                     &&
                     (and $ map maybeInUnit $ DBox.elems peak)
                    ,
                     erintv_right $
                     chplRAEval makeInterval p peak
                    )
                Nothing -> (False, undefined)
            where
            noPositiveSquare =
                and $ map (<= 0) $ map getQuadCoeff vars
            getQuadCoeff var = 
                Map.findWithDefault 0 (DBox.singleton var 2) coeffs
            maybeInUnit r =
                case (RA.compareReals r (-1), RA.compareReals (1) r) of
                    (Just LT, _) -> False -- ie r < -1
                    (_, Just LT) -> False -- ie r > 1
                    _ -> True
        maybePeak =
            linearSolver
                (map derivZeroLinearEq vars)
                (DBox.fromList $ map (\v -> (v,(-1) RA.\/ 1)) vars)
                (2^^(-ix))
            where
            derivZeroLinearEq var =
                (linCoeffs, - constCoeff)
                where
                constCoeff =
                    makeInterval $
                    Map.findWithDefault 0 (DBox.singleton var 1) coeffs
                      -- recall T_1(x) = x, T_1'(x) = 1
                linCoeffs =
                    DBox.fromList $
                        (var, 4 * quadCoeff) -- T_2(x) = 2*x^2 - 1; T_2'(x) = 4*x
                        : (map getVarVarCoeff $ var `delete` vars)
                quadCoeff =
                    makeInterval $
                    Map.findWithDefault 0 (DBox.singleton var 2) coeffs
                getVarVarCoeff var2 =
                    (var2,
                      makeInterval $
                      Map.findWithDefault 0 (DBox.fromList [(var,1), (var2,1)]) coeffs)
        makeInterval b = ERInterval b b
        removeVar var =
            [(substVar True, newVars), 
             (substVar False, newVars)]
            where
            newVars = var `delete` vars
            substVar isOne =
                chplCoeffs $
                    foldl (+^) (chplConst 0) $ 
                        map (makeMonomial isOne) $ 
                            Map.toList coeffs
            makeMonomial isOne (term, coeff) =
                ERChebPoly $ Map.fromList $
                case (DBox.toList term) of
                    [(v,2)] | v == var ->
                        [(chplConstTermKey, coeff)]
                    [(v,1)] | v == var ->
                        [(chplConstTermKey, 
                          case isOne of True -> coeff; False -> - coeff)]
                    [(v1,1), (v2,1)] | v1 == var ->
                        [(DBox.fromList [(v2,1)], 
                          case isOne of True -> coeff; False -> - coeff)]
                    [(v1,1), (v2,1)] | v2 == var ->
                        [(DBox.fromList [(v1,1)], 
                          case isOne of True -> coeff; False -> - coeff)]
                    _ ->
                        [(term, coeff)]

{-|
     Approximate from below and  from above the pointwise maximum of two polynomials
-}
chplMax ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplMax maxDegree maxSize p1 p2 =
    (p1 +. differenceDown, p1 +^ differenceUp)
    where
    (differenceDown, _) = chplNonneg maxDegree maxSize p2MinusP1Down
    (_, differenceUp) = chplNonneg maxDegree maxSize $ p2MinusP1Up
    (p2MinusP1Down, p2MinusP1Up, _) = chplAdd p2 (chplNeg p1)

chplMaxDn m s a b = fst $ chplMax m s a b
chplMaxUp m s a b = snd $ chplMax m s a b
chplMinDn m s a b = fst $ chplMin m s a b
chplMinUp m s a b = snd $ chplMin m s a b

{-|
     Approximate from below and  from above the pointwise minimum of two polynomials
-}
chplMin ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplMin m s a b =
    (chplNeg u,chplNeg l)
    where
    (l,u) = chplMax m s (chplNeg a) (chplNeg b)

chplNonnegDown m s p = fst $ chplNonneg m s p
chplNonnegUp m s p = snd $ chplNonneg m s p 
chplNonposDown m s p = fst $ chplNonpos m s p
chplNonposUp m s p = snd $ chplNonpos m s p 

chplNonpos m s p =
    (chplNeg h, chplNeg l)
    where
    (l,h) = chplNonneg m s (chplNeg p)

{-|
     Approximate the function max(0,p(x)) by a polynomial from below
     and from above. 
-}
chplNonneg ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplNonneg = chplNonnegCubic

{-|
    A version of 'chplNonneg' using a cubic approximation. 
-}
chplNonnegCubic ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplNonnegCubic maxDegree maxSize p
    | upperB <= 0 = (chplConst 0, chplConst 0)
    | lowerB >= 0 = (p, p)
    | not allInterimsBounded = (chplConst (1/0), chplConst (1/0))
    | otherwise = -- ie lowerB < 0 < upperB: polynomial may be crossing 0...
--        unsafePrintReturn
--        (
--            "chplNonnegCubic:"
--            ++ "\n p = " ++ show p
--            ++ "\n maxDegree = " ++ show maxDegree
--            ++ "\n maxSize = " ++ show maxSize
--            ++ "\n upperB = " ++ show upperB
--            ++ "\n lowerB = " ++ show lowerB
--            ++ "\n a0 = " ++ show a0
--            ++ "\n a1 = " ++ show a1
--            ++ "\n a2 = " ++ show a2
--            ++ "\n a3 = " ++ show a3
--            ++ "\n b = " ++ show b
--            ++ "\n rb = " ++ show rb
--            ++ "\n correctionB = " ++ show correctionB
--            ++ "\n valueAt0B = " ++ show valueAt0B
--            ++ "\n result = "
--        )
        -- work out the cubic polynomial (a3*x^3 + a2*x^2 + a1*x + a0) / b 
        -- that hits 0 at lowerB with derivative 0 
        -- and hits upperB at upperB with derivative 1 
        (chplAddConstDown (- valueAt0B) cubicAppliedOnPDown, 
         chplAddConstUp correctionB cubicAppliedOnPUp)
    where    
    (lowerB, upperB) = chplBounds 10 p
    (cubicAppliedOnPDown, cubicAppliedOnPUp, width) =
        p0 `scaleByPositiveConsts` (rbLo, rbHi)
        where
        p0 = (multiplyByP p1) `addConsts` (a0Lo, a0Hi) -- ie p*(p*(p * a3 + a2) + a1) + a0 enclosure
        p1 = (multiplyByP p2) `addConsts` (a1Lo, a1Hi) -- ie p*(p * a3 + a2) + a1 enclosure
        p2 = (multiplyByP p3) `addConsts` (a2Lo, a2Hi) -- ie p * a3 + a2 enclosure
        p3 = (chplConst a3Lo, chplConst a3Hi, a3Hi - a3Lo) -- ie a3 enclosure
    multiplyByP (lo,hi,wd) =
        (ploRed, phiRed, pwd)
        where
        ploRed = reduceDgSzDown plo
        phiRed = reduceDgSzUp phi 
        pwd = chplUpperBound 10 $ phiRed -^ ploRed 
        (plo, phi, _) = chplTimesLoHi p (lo,hi,wd)
    reduceDgSzUp =
        chplReduceTermCountUp maxSize . chplReduceDegreeUp maxDegree
    reduceDgSzDown =
        chplReduceTermCountDown maxSize . chplReduceDegreeDown maxDegree
    addConsts (lo, hi, wd) (cLo, cHi) =
        (alo, ahi, wd + wdlo + wdhi)
        where
        (alo, _, wdlo) = chplAddConst cLo lo 
        (_, ahi, wdhi) = chplAddConst cHi hi 
    scaleByPositiveConsts (lo, hi, wd) (cLo, cHi) =
        (slo, shi, wd + wdlo + wdhi)
        where
        (slo, _, wdlo) = chplScale cLo lo 
        (_, shi, wdhi) = chplScale cHi hi 
    
    -- convert interval coefficients to pairs of bounds:
    ERInterval rbLo rbHi = rb
    ERInterval a3Lo a3Hi = a3
    ERInterval a2Lo a2Hi = a2
    ERInterval a1Lo a1Hi = a1
    ERInterval a0Lo a0Hi = a0
    allInterimsBounded = 
        and $ map RA.isBounded [w, s, rb, a0, a1, a2, a3, correction]
    {-
      The cubic polynomial's coefficients are calculated by solving a system of 4 linear eqs.
      The generic solution is as follows:
         b = (r - l)^3   always positive
         a3 = -(r + l)
         a2 = 2*(r^2 + r*l + l^2)
         a1 = -l*(4*r^2 + r*l + l^2)
         a0 = 2*r^2*l^2
    -}
    rb = recip b
    b = w3 -- = w^3 -- see below
    w = r - l
    r = ERInterval upperB upperB
    l = ERInterval lowerB lowerB
    --
    a3 = - s
    s = r + l
    --
    a2 = 2 * (r2PrlPl2)
    r2PrlPl2 = s2 - rl
    rl = r * l
    --
    a1 = (- l) * (r2PrlPl2 + 3*r2)
    a0 = 2*r2*l2
    -- interval arithmetic tricks to speed it up and reduce dependency errors:
    w3 = ERInterval (wLo * wLo * wLo) (wHi * wHi * wHi) -- x^3 is monotone 
    ERInterval wLo wHi = w
    s2 = ERInterval (max 0 s2Lo) s2Hi
    s2Lo = min sLo2 sHi2 
    s2Hi = max sLo2 sHi2
    sLo2 = sLo * sLo
    sHi2 = sHi * sHi 
    ERInterval sLo sHi = s    
    r2 = ERInterval (upperB `timesDown` upperB) (upperB `timesUp` upperB)    
    l2 = ERInterval (lowerB `timesDown` lowerB) (lowerB `timesUp` lowerB)
    {- 
        The cubic polynomial may sometimes fail to dominate
        x or sometimes it dips below 0.
        Work out the amount by which it has to be lifted up
        to fix these problems. 
    -}
    ERInterval _ correctionB = correction
    correction =
        case (RA.compareReals (2 * r2) (l*s), RA.compareReals (2 * l2) (r*s)) of
            (Just LT, _) ->
                (peak0 * (peak0 * (peak0 * (-a3) - a2) - a1) - a0) / b
            (_, Just LT) ->
                ((peakP * (peakP * (peakP * (-a3) - a2) - a1) - a0) / b) + peakP
            _ -> 0
        where
        peak0 = (l + 4*r2/s) / 3 
        peakP = (r + 4*l2/s) / 3
    {-
        The same cubic polynomial can be used as a lower bound when
        we subtract its value at 0 rounded upwards.
    -}
    valueAt0B = 
        case a0 / b of
            ERInterval lo hi -> hi
            ERIntervalAny -> 1/0

{-|
    Multiply a polynomial by an enclosure (with non-negated lower bound).
-}
chplTimesLoHi ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b, b) ->
    (ERChebPoly box b, ERChebPoly box b, b)
chplTimesLoHi p1 (p2Low, p2High, p2Width) =
    (prodMid -. (chplConst width), 
     prodMid +^ (chplConst width), 
     2 * width)
    where
    prodMid = prodLowUp
    (prodLowDown, prodLowUp, prodLowWidth) = 
        chplMultiply p1 p2Low
    (prodHighDown, prodHighUp, prodHighWidth) = 
        chplMultiply p1 p2High
    width = 
        p1Norm `timesUp` p2Width `plusUp` prodLowWidth `plusUp` prodHighWidth
    p1Norm = 
        max (abs $ p1LowerBound) (abs $ p1UpperBound)
    (p1LowerBound, p1UpperBound) = 
        chplBounds ix p1
    ix = 10