AERN-RnToRm-0.5: src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/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.BasicTypes.DomainBox as DBox
import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainIntBox, DomainBoxMappable)
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) =
(min lLower hLower, max lUpper hUpper)
where
(lLower, lUpper) = chplBounds ix $ chplNeg ln
(hLower, hUpper) = chplBounds ix h
enclBoundsExpensive ix (ln, h) =
(negate $ chplUpperBoundExpensive ix ln, chplUpperBoundExpensive 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
= result
where
result = ERInterval lB hB
lB = negate $ chplEvalUp ln pt
hB = chplEvalUp h pt
enclRAEval e@(ln, h) pt =
result
where
result = ERInterval lAtPt hAtPt
ERInterval lAtPt _ = negate $ chplRAEval (\b -> ERInterval b b) ln pt
ERInterval _ hAtPt = chplRAEval (\b -> ERInterval b b) h pt
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)
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)
enclAdd ::
(B.ERRealBase 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)
enclAdd maxDegr maxSize (p1LowNeg, p1High) (p2LowNeg, p2High) =
enclReduceSize maxSize $
(p1LowNeg +^ p2LowNeg, p1High +^ p2High)
enclMultiply ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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
-- unsafePrint "both non-negative" $
(l1l2Neg, h1h2)
(_, True, _, True) -> -- both non-positive
-- unsafePrint "both non-positive" $
(h1h2Neg, l1l2)
(True, _, _, True) -> -- first non-negative, second non-positive
-- unsafePrint "first non-negative, second non-positive" $
(h1l2Neg, l1h2)
(_, True, True, _) -> -- first non-positive, second non-negative
-- unsafePrint
-- ("ERChebPoly: enclMultiply: first non-positive, second non-negative:"
-- ++ "\n l1 = " ++ show (chplNeg ln1)
-- ++ "\n h1 = " ++ show h1
-- ++ "\n l2 = " ++ show (chplNeg ln2)
-- ++ "\n h2 = " ++ show h2
-- ) $
(l1h2Neg, h1l2)
_ -> -- 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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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) = chplBall2DownUp $ ballMultiply ln ln
(h2Down, h2Up) = chplBall2DownUp $ ballMultiply 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, 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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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 ratio@(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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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