Enclosure.hs

{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure
    Description :  (internal) field operations applied to 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 selected operations working on pairs
    of polynomials understood as function enclosures.
    These are needed to define composition and some elementary operations.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure

where

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

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

import qualified Data.Map as Map

enclThin ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
enclThin p =
    (chplNeg p, p)

enclConst ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    b ->
    (ERChebPoly box b, ERChebPoly box b)
enclConst c =
    (chplConst (-c), chplConst c)

enclBounds ix (ln, h) =
    (negate $ chplUpperBound ix ln, chplUpperBound ix h)

enclEval e@(ln, h) pt 
    | lB > hB =
        unsafePrintReturn
        (
            "ERChebPoly: enclEval: inverted result:"
            ++ "\n h = " ++ show h 
            ++ "\n ln = " ++ show ln 
            ++ "\n result = "
        )
        result
    | otherwise = result
    where
    result = ERInterval lB hB
    lB = negate $ chplEvalUp ln pt
    hB = chplEvalUp h pt

enclEvalInner (ln, h) pt =
--    normaliseERInterval $
    ERInterval 
        (negate $ chplEvalDown ln pt)
        (chplEvalDown h pt)

enclRAEval e@(ln, h) pt =
    result 
    where
    result = lRA RA.\/ hRA
    lRA = fst $ RA.bounds $ negate $ chplRAEval (\b -> ERInterval b b) ln pt
    hRA = snd $ RA.bounds $ chplRAEval (\b -> ERInterval b b) h pt

enclRAEvalInner e@(ln, h) pt =
--    unsafePrintReturn
--    (
--        "ERChebPoly: enclRAEvalInner: "
--        ++ "\n lB = " ++ show lB
--        ++ "\n hB = " ++ show hB
--        ++ "\n result = "
--    )
    result 
    where
    result =
--        normaliseERInterval $ 
        ERInterval lB hB
    lB = 
        case negate $ chplRAEval (\b -> ERInterval b b) ln pt of
            ERInterval _ lB -> lB
    hB = 
        case chplRAEval (\b -> ERInterval b b) h pt of
            ERInterval hB _ -> hB

enclAddErr errB (pLowNeg, pHigh) =
    (chplAddConstUp errB pLowNeg, chplAddConstUp errB pHigh)


enclRAConst ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    (ERInterval b) ->
    (ERChebPoly box b, ERChebPoly box b)
enclRAConst (ERInterval lo hi) = (chplConst (-lo), chplConst hi)
enclRAConst ERIntervalAny = (chplConst (-1/0), chplConst (1/0))

enclReduceDegree maxDegree (pLowNeg, pHigh) =
    (chplReduceDegreeUp maxDegree pLowNeg, chplReduceDegreeUp maxDegree pHigh)  
    
enclReduceSize maxSize (pLowNeg, pHigh) =
    (chplReduceTermCountUp maxSize pLowNeg, chplReduceTermCountUp maxSize pHigh)  
    
enclAddConst c (pLowNeg, pHigh) =
    (chplAddConstUp (-c) pLowNeg, chplAddConstUp c pHigh)

enclNeg (pLowNeg, pHigh) = (pHigh, pLowNeg)

(p1LowNeg, p1High) +: (p2LowNeg, p2High) = 
    (p1LowNeg +^ p2LowNeg, p1High +^ p2High)
    
(p1LowNeg, p1High) -: (p2LowNeg, p2High) =
    (p1LowNeg +^ p2High, p1High +^ p2LowNeg)

enclMultiply ::
    (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) ->
    (ERChebPoly box b, ERChebPoly box b)
enclMultiply maxDegr maxSize (ln1, h1) (ln2, h2) =
    enclReduceSize maxSize $
    enclReduceDegree maxDegr $
    case (ln1UpperBound <= 0, h1UpperBound <= 0, ln2UpperBound <= 0, h2UpperBound <= 0) of
        (True, _, True, _) -> -- both non-negative
            (l1l2Neg, h1h2)
        (_, True, _, True) -> -- both non-positive
            (h1h2Neg, l1l2)
        (True, _, _, True) -> -- first non-negative, second non-positive
            (h1l2Neg, l1h2)
        (_, True, True, _) -> -- first non-positive, second non-negative
            (l1h2Neg, l1h2)
        _ -> -- one of both may be crossing zero
            (
             (h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg)
            ,
             (h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2)
            )
        where
        ln1UpperBound = chplUpperBound ix ln1
        ln2UpperBound = chplUpperBound ix ln2
        h1UpperBound = chplUpperBound ix h1
        h2UpperBound = chplUpperBound ix h2
        ix = 10
        maxP = chplMaxUp maxDegr maxSize
        
        h1h2 = h1 *^ h2
        h1h2Neg = (chplNeg h1) *^ h2
        l1l2 = ln1 *^ ln2
        l1l2Neg = (chplNeg ln1) *^ ln2
        h1l2 = h1 *^ (chplNeg ln2)
        h1l2Neg = h1 *^ ln2
        l1h2 = (chplNeg ln1) *^ h2
        l1h2Neg = ln1 *^ h2


enclSquare ::
    (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)
enclSquare maxDegr maxSize (ln, h)
    {-
        formula:
            (ln, h)^2 =
                ( minUp( 0, maxUp( - ln *. ln, - h *. h)), maxUp(ln *^ ln, h *^ h) )
    -}
--    | minZeroHelps = 
    = (minZeroMaxNegSq, maxSq)
--    | otherwise =
--        (maxNegSq, maxSq)
    where
    maxSq = maxP ln2Up h2Up
    maxNegSq = maxP (chplNeg ln2Down) (chplNeg h2Down)
    minZeroMaxNegSq = chplNonposUp maxDegr maxSize maxNegSq 
--    minZeroHelps =
--        (maxNegSqUpperB > 0) && (minZeroMaxNegSqUpperB / maxNegSqUpperB < 1/2)
--    maxNegSqUpperB = chplUpperBound 10 maxNegSq
--    minZeroMaxNegSqUpperB = chplUpperBound 10 minZeroMaxNegSq
     
    (ln2Down, ln2Up, _) = chplMultiply ln ln
    (h2Down, h2Up, _) = chplMultiply h h
    
--    reduceDegrSize = reduceSize maxSize . reduceDegree maxDegr
    maxP = chplMaxUp maxDegr maxSize
    
    

    
{-| 
    Multiply an enclosure by a scalar 
    assuming the enclosure is non-negative on the whole unit domain.
-} 
enclScaleNonneg ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    b {-^ ratio to scale by -} -> 
    (ERChebPoly box b, ERChebPoly box b) -> 
    (ERChebPoly box b, ERChebPoly box b)
enclScaleNonneg ratio pEncl@(ln, h) =
    (ln *^ pRatio, h *^ pRatio)
    where
    pRatio = chplConst ratio

{-| 
    Multiply an enclosure by a scalar.
-} 
enclScale ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} ->
    b {-^ ratio to scale by -} -> 
    (ERChebPoly box b, ERChebPoly box b) -> 
    (ERChebPoly box b, ERChebPoly box b)
enclScale maxDegree maxSize ratio pEncl =
    enclMultiply maxDegree maxSize pEncl (enclConst ratio) 

enclRAScale ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} ->
    (ERInterval b) -> 
    (ERChebPoly box b, ERChebPoly box b) ->
    (ERChebPoly box b, ERChebPoly box b)
enclRAScale maxDegree maxSize ra pEncl =
    enclMultiply maxDegree maxSize pEncl (enclRAConst ra) 

{-|
    Multiply a polynomial by a scalar interval, returning an enclosure.
-} 
chplScaleRA :: 
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERInterval b {-^ lower and upper bounds on the ratio to scale by -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplScaleRA maxDegr maxSize (ERIntervalAny) p = enclRAConst ERIntervalAny
chplScaleRA maxDegr maxSize (ERInterval ratioDown ratioUp) p =
    (scaledPDownNeg, scaledPUp)
    where
    (scaledPDownNeg, scaledPUp) =
        enclMultiply maxDegr maxSize 
            (chplNeg p, p) (chplConst (- ratioDown), chplConst ratioUp)

chplScaleRADown m n r = chplNeg . fst . chplScaleRA m n r
chplScaleRAUp m n r = snd . chplScaleRA m n r

{-|
    Evaluate the Chebyshev polynomials of the first kind
    applied to a given polynomial, yielding a list of polynomial enclosures. 
-}
enclEvalTs ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
    Int {-^ max degree for result -} -> 
    Int {-^ max approx size for result -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ bounds of a polynomial enclosure to evaluate -} ->
    [(ERChebPoly box b, ERChebPoly box b)]
enclEvalTs maxDegree maxSize p1@(pLowNeg, pHigh) =
    chebyIterate (enclConst 1) p1
    where
    chebyIterate pNm2 pNm1 =
        pNm2 : (chebyIterate pNm1 pN)
        where
        pN = 
            (enclScale maxDegree maxSize 2 $ 
                enclMultiply maxDegree maxSize p1 pNm1) 
            -: pNm2

{-|
    Multiply a polynomial by an enclosure using min/max
-}
enclThinTimes ::
    (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, ERChebPoly box b)
enclThinTimes maxDegree maxSize p1 (p2LowNeg, p2High) =
    (prodLowNeg, prodHigh)
    where
    prodHigh =
        chplMaxUp maxDegree maxSize
            (p1 *^ p2High)
            (p1n *^ p2LowNeg) -- beware: p1 can cross zero
    prodLowNeg =
        chplMaxUp maxDegree maxSize
            (p1n *^ p2High)
            (p1 *^ p2LowNeg)
    p1n = chplNeg p1