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 Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Derivative
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.BasicTypes.DomainBox as DBox
import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc
import qualified Data.Map as Map
import Data.List
chplUpperBound ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
b
chplUpperBound ix p = snd $ chplBounds ix p
chplLowerBound ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
b
chplLowerBound ix p = fst $ chplBounds ix p
chplBounds ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
(b,b)
chplBounds =
chplBoundsAffine
chplUpperBoundExpensive ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
b
chplUpperBoundExpensive ix p = snd $ chplBoundsExpensive ix p
chplLowerBoundExpensive ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
b
chplLowerBoundExpensive ix p = fst $ chplBoundsExpensive ix p
chplBoundsExpensive ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
(b,b)
chplBoundsExpensive = chplBoundsByDerivative
chplBoundsAffine ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
EffortIndex ->
ERChebPoly box b ->
(b,b)
chplBoundsAffine ix p@(ERChebPoly coeffs) =
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
chplBoundsByDerivative ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
EffortIndex ->
ERChebPoly box b ->
(b,b)
chplBoundsByDerivative ix p =
(lowerBound, upperBound)
where
lowerBound = foldl1 min $ map fst extremaValues
upperBound = foldl1 max $ map snd extremaValues
ra2bb (ERInterval l r) = (l,r)
b2ra b = ERInterval b b
extremaValues =
collectValuesOnFaces vars varDerivatives (p,p)
where
vars = chplGetVars p
varDerivatives =
Map.fromList $
map getDerivatives vars
getDerivatives var =
(var,
chplBall2DownUp $
ballDifferentiate p var)
collectValuesOnFaces varsSpecialise varDerivatives (pDown, pUp) =
valuesThisFace ++ (valuesSubFaces varsSpecialise)
where
valuesThisFace =
collectExtremeValues varDerivatives (pDown, pUp)
valuesSubFaces [] = []
valuesSubFaces (var : vars) =
(collectValuesOnFaces vars varDerivativesNoVarL (pDownNoVarL, pUpNoVarL))
++
(collectValuesOnFaces vars varDerivativesNoVarR (pDownNoVarR, pUpNoVarR))
++
(valuesSubFaces vars)
where
(pDownNoVarR, pUpNoVarR) = substVarR (pDown, pUp)
(pDownNoVarL, pUpNoVarL) = substVarL (pDown, pUp)
substVarL = substVar (1)
substVarR = substVar 1
substVar val (pDown, pUp) =
(fst $ chplPartialRAEval ra2bb pDown $ DBox.singleton var val,
snd $ chplPartialRAEval ra2bb pUp $ DBox.singleton var val)
varDerivativesNoVarL =
Map.map substVarL varDerivativesNoVar
varDerivativesNoVarR =
Map.map substVarR varDerivativesNoVar
varDerivativesNoVar =
Map.delete var varDerivatives
collectExtremeValues varDerivatives (pDown, pUp)
| null varsNoConst =
[pEvalAt unitDomBox]
| otherwise =
map pEvalAt boxesWithPotentialExtrema
where
boxesWithPotentialExtrema =
paveFindBoxes [(unitDomBox,0)]
varDerivativesNoZeros =
Map.filter (not . isConstWithZero) varDerivatives
where
isConstWithZero (pDown, pUp) =
(snd $ chplBoundsAffine ix pDown) <= 0
&&
(fst $ chplBoundsAffine ix pUp) >= 0
vars = Map.keys varDerivatives
varsNoConst = Map.keys varDerivativesNoZeros
varsNoConstLength = length varsNoConst
pEvalAt = evalAt (pDown, pUp)
evalAt (pDown,pUp) box =
(fst $ ra2bb $ chplRAEval b2ra pDown box,
snd $ ra2bb $ chplRAEval b2ra pUp box)
unitDomBox =
DBox.fromList $ zip vars (repeat unitInterval)
unitInterval = ((1) RA.\/ 1)
maxDepth = fromInteger $ toInteger $ max 3 ix
keepBox box =
and $ map evalDeriv $ Map.elems varDerivativesNoZeros
where
evalDeriv derivBounds = hasZero $ evalAt derivBounds box
hasZero (l,h) = l <= 0 && h >= 0
paveFindBoxes [] = []
paveFindBoxes boxes@((box, depth) : boxesRest)
| keepBox box =
case depth < maxDepth of
True ->
paveFindBoxes ((boxL, newDepth) : (boxR, newDepth) : boxesRest)
False ->
box : (paveFindBoxes boxesRest)
| otherwise =
paveFindBoxes boxesRest
where
var = varsNoConst !! (depth `mod` varsNoConstLength)
(boxL, boxR) = DBox.split box var Nothing
newDepth = depth + 1
--{-|
---}
chplMax ::
(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 ->
Int ->
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) = chplBall2DownUp $ ballAdd 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
chplMin ::
(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 ->
Int ->
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)
chplNonneg ::
(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 ->
Int ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplNonneg = chplNonnegCubic
chplNonnegCubic ::
(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 ->
Int ->
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 (B.plusInfinity), chplConst (B.plusInfinity))
| otherwise =
(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)
p1 = (multiplyByP p2) `addConsts` (a1Lo, a1Hi)
p2 = (multiplyByP p3) `addConsts` (a2Lo, a2Hi)
p3 = (chplConst a3Lo, chplConst a3Hi, a3Hi a3Lo)
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) = chplBall2DownUpWd $ ballAddConst cLo lo
(_, ahi, wdhi) = chplBall2DownUpWd $ ballAddConst cHi hi
scaleByPositiveConsts (lo, hi, wd) (cLo, cHi) =
(slo, shi, wd + wdlo + wdhi)
where
(slo, _, wdlo) = chplBall2DownUpWd $ ballScale cLo lo
(_, shi, wdhi) = chplBall2DownUpWd $ ballScale cHi hi
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]
rb = recip b
b = w3
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
w3 = ERInterval (wLo * wLo * wLo) (wHi * wHi * wHi)
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)
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
valueAt0B =
case a0 / b of
ERInterval lo hi -> hi
chplTimesLoHi ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
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) =
chplBall2DownUpWd $ ballMultiply p1 p2Low
(prodHighDown, prodHighUp, prodHighWidth) =
chplBall2DownUpWd $ ballMultiply p1 p2High
width =
p1Norm `timesUp` p2Width `plusUp` prodLowWidth `plusUp` prodHighWidth
p1Norm =
max (abs $ p1LowerBound) (abs $ p1UpperBound)
(p1LowerBound, p1UpperBound) =
chplBounds ix p1
ix = 10