AERN-RnToRm-0.4: src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs
{-# 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)]
chplMaxDn m a b = fst $ chplMax m a b
chplMaxUp m a b = snd $ chplMax m a b
chplMinDn m a b = fst $ chplMin m a b
chplMinUp m a b = snd $ chplMin m a b
chplMin m a b =
(-u,-l)
where
(l,u) = chplMax m (-a) (-b)
{-|
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 `plusDown` differenceDown, p1 `plusUp` 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 =
chplReduceDegreeUp maxDegree . (p *)
multiplyByPDown =
chplReduceDegreeDown 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
{-|
Multiply a thin enclosure by a non-thin enclosure
-}
chplThinTimesEncl ::
(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, ERChebPoly box b)
chplThinTimesEncl maxDegree p1 (p2LO, p2HI) =
(prodLO, prodHI)
where
prodHI =
chplMaxUp maxDegree
(p1 `timesUp` p2HI)
(p1 `timesUp` p2LO) -- beware: p1 can be negative
prodLO =
negate $
chplMaxUp maxDegree
(p1n `timesUp` p2HI)
(p1n `timesUp` p2LO)
p1n = negate p1
{-|
Safely multiply a polynomial by itself.
-}
chplSquare ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplSquare maxDegree p =
(p `timesDown` p, p `timesUp` p)