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