AERN-RnToRm-0.4.9: 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.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