```{-# 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

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.
-}
(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
where
vars = chplGetVars p
| 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
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)
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
```