{-# 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.Eval
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field

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.  
-}
chplUpperBoundAffine ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplUpperBoundAffine ix (ERChebPoly coeffs) =
    affiBound coeffs
    where
    affiBound coeffs =
        Map.fold (+) constTerm absCoeffs
        where
        absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs
        constTerm = Map.findWithDefault 0 chplConstTermKey coeffs


{-|
    Find a close upper bound on an affine polynomial over the 
    unit domain [-1,1]^n.  
-}
chplUpperBoundAffineCorners ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box,
     DomainBoxMappable boxb boxbb varid b [(b,b)], Num varid, Enum varid) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplUpperBoundAffineCorners ix p@(ERChebPoly coeffs) =
    affiBound (coeffs, vars)
    where
    vars = chplGetVars p
    affiBound (coeffs, vars)
        | null vars =
            Map.findWithDefault 0 chplConstTermKey coeffs
        | otherwise =
            foldl1 max cornerValues
        where
        cornerValues =
            map (\pt -> chplEvalUp pt p) corners
            where
--            corners :: [boxb]
            corners = 
                map (DBox.fromList . (zip [1..n])) $ prod n
                where
                n = fromInteger $ toInteger $ length vars
                -- n-fold product list of [-1,1]
                prod n 
                    | n == 1 = [[-1],[1]]
                    | otherwise =
                        (map ((-1):) prodNm1) ++ (map (1:) $ prodNm1)
                    where
                    prodNm1 = prod (n-1)

{-|
    Find a close upper bound on a quadratic polynomial over the 
    unit domain [-1,1]^n.  
-}
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 (coeffs, vars)
    where
    vars = chplGetVars p
    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 $
                     chplEvalApprox makeInterval peak p      
                    )
                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 $
                    sum $ 
                        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, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplMax maxDegree p1 p2 =
    (- (-p1 - differenceDown), p1 + differenceUp)
    where
    (differenceDown, differenceUp) = chplNonneg maxDegree $ p2 - p1

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

{-|
    A version of 'chplNonneg' using a cubic approximation. 
-}
chplNonnegCubic ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplNonnegCubic maxDegree p
    | upperB <= 0 = (chplConst 0, chplConst 0)
    | lowerB >= 0 = (p, p)
    | otherwise = -- ie lowerB < 0 < upperB: polynomial may be crossing 0...
        -- 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 
        (cubicAppliedOnPDown - valueAt0, cubicAppliedOnPUp + (chplConst correction))
    where    
    upperB = chplUpperBoundAffine 10 p    
    lowerB = - (chplUpperBoundAffine 10 (- p))
    cubicAppliedOnPUp = evalCubic multiplyByPUp
    cubicAppliedOnPDown = evalCubic multiplyByPDown
    evalCubic multiplyByP =
        p0 * (chplConst $ recip b)
        where
        p0 = multiplyByP p1 + (chplConst a0) -- ie p*(p*(p * a3 + a2) + a1) + a0
        p1 = multiplyByP p2 + (chplConst a1) -- ie p*(p * a3 + a2) + a1
        p2 = multiplyByP p3 + (chplConst a2) -- ie p * a3 + a2
        p3 = chplConst a3
    multiplyByPUp =
        snd . chplReduceDegree maxDegree . (p *)
    multiplyByPDown =
        fst . chplReduceDegree maxDegree . (p *)
    {-
      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
         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
    -}
    r = upperB
    l = lowerB
    b = - ((r - l) * ((r - l) * (l - r))) 
        -- this one has to round downwards because it is a denominator
    a3 = (- r) + (- l) -- remember to round upwards!
    a2 = 2*(r2rll2Up)
    a1 = (- l) * (r2rll2Up + 3*rSqUp) -- since l < 0, the other argument is rounded upwards
    a0 = 2 * rSqUp * lSqUp
    r2rll2Up = rSqUp + r*l + lSqUp 
    rSqUp = r*r
    lSqUp = l*l
    rSqDown = -((-r)*r)
    lSqDown = -((-l)*l)
    {- 
        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. 
    -}
    correction
        | 2*rSqDown < l*(r + l) =
            erintv_right $
            (peak0 * (peak0 * (peak0 * (-a3I) - a2I) - a1I) - a0I) / bI
        | 2*lSqDown < r*(r + l) =
            erintv_right $
            ((peakP * (peakP * (peakP * (-a3I) - a2I) - a1I) - a0I) / bI) + peakP
        | otherwise = 0
        where
        -- these have to be computed interval-based:
        [a0I, a1I, a2I, a3I, bI, lI, rI] = 
            map (\x -> ERInterval x x) [a0,a1,a2,a3,b,l,r]
        peak0 = (lI + 4*rI*rI/(lI+rI)) / 3 
        peakP = (rI + 4*lI*lI/(lI+rI)) / 3
    {-
        The same cubic polynomial can be used as a lower bound when
        we subtract its value at 0 rounded upwards.
    -}
    valueAt0 = chplConst $ a0 / b