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
chplUpperBoundAffine ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
EffortIndex ->
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
chplUpperBoundAffineCorners ::
(B.ERRealBase b, DomainBox box varid Int, Ord box,
DomainBoxMappable boxb boxbb varid b [(b,b)], Num varid, Enum varid) =>
EffortIndex ->
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 =
map (DBox.fromList . (zip [1..n])) $ prod n
where
n = fromInteger $ toInteger $ length vars
prod n
| n == 1 = [[1],[1]]
| otherwise =
(map ((1):) prodNm1) ++ (map (1:) $ prodNm1)
where
prodNm1 = prod (n1)
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 ->
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
&&
(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
(_, Just LT) -> False
_ -> 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
linCoeffs =
DBox.fromList $
(var, 4 * quadCoeff)
: (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)]
chplMax ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
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
chplNonneg ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplNonneg = chplNonnegCubic
chplNonnegCubic ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplNonnegCubic maxDegree p
| upperB <= 0 = (chplConst 0, chplConst 0)
| lowerB >= 0 = (p, p)
| otherwise =
(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)
p1 = multiplyByP p2 + (chplConst a1)
p2 = multiplyByP p3 + (chplConst a2)
p3 = chplConst a3
multiplyByPUp =
snd . chplReduceDegree maxDegree . (p *)
multiplyByPDown =
fst . chplReduceDegree maxDegree . (p *)
r = upperB
l = lowerB
b = ((r l) * ((r l) * (l r)))
a3 = ( r) + ( l)
a2 = 2*(r2rll2Up)
a1 = ( l) * (r2rll2Up + 3*rSqUp)
a0 = 2 * rSqUp * lSqUp
r2rll2Up = rSqUp + r*l + lSqUp
rSqUp = r*r
lSqUp = l*l
rSqDown = ((r)*r)
lSqDown = ((l)*l)
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
[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
valueAt0 = chplConst $ a0 / b