{-# OPTIONS_GHC -fno-warn-missing-methods #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE ScopedTypeVariables #-} {-| Module : Data.Number.ER.RnToRm.UnitDom.Approx.Interval Description : arbitrary precision function enclosures on @[-1,1]^n@ Copyright : (c) Michal Konecny License : BSD3 Maintainer : mik@konecny.aow.cz Stability : experimental Portability : portable A construction of an enclosure of a real function on the domain [-1,1]^n for some n using elements of some base (eg rational functions or polynomials). -} module Data.Number.ER.RnToRm.UnitDom.Approx.Interval ( ERFnInterval(..), ERFnContext(..), erfnContextDefault, erfnContextUnify ) where import qualified Data.Number.ER.Real.Base as B import Data.Number.ER.Real.Approx.Interval import Data.Number.ER.Real.Arithmetic.Elementary import qualified Data.Number.ER.RnToRm.Approx as FA import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB import Data.Number.ER.RnToRm.UnitDom.Base ((+^),(-^),(*^),multiplyEncl) import qualified Data.Number.ER.Real.Approx as RA import qualified Data.Number.ER.Real.Approx.Elementary as RAEL import qualified Data.Number.ER.BasicTypes.DomainBox as DBox import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainIntBox) import Data.Number.ER.BasicTypes import Data.Number.ER.Misc import Data.Number.ER.ShowHTML import qualified Text.Html as H import qualified Data.Map as Map import Data.Typeable import Data.Generics.Basics import Data.Binary {- only for testing in ghci, to be removed: -} --import Data.Number.ER.Real.DefaultRepr --import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom --import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.PolynomBase --type FAPU = ERFnInterval (ERChebPoly B) IRA --fapuConst1 = (UFA.const 0 [1]) :: FAPU --fapuConst2 = (UFA.const 0 [2]) :: FAPU {- end of testing specific code -} data ERFnInterval fb = ERFnIntervalAny { erfnContext :: ERFnContext } | ERFnInterval { erfnLowerNeg :: fb, erfnUpper :: fb, erfnContext :: ERFnContext -- , -- erfnIsDefinitelyConsistent :: Bool, -- erfnIsDefinitelyAntiConsistent :: Bool } deriving (Typeable, Data) instance (Binary a) => Binary (ERFnInterval a) where put (ERFnIntervalAny a) = putWord8 0 >> put a put (ERFnInterval a b c) = putWord8 1 >> put a >> put b >> put c -- put (ERFnInterval a b c d e) = putWord8 1 >> put a >> put b >> put c >> put d >> put e get = do tag_ <- getWord8 case tag_ of 0 -> get >>= \a -> return (ERFnIntervalAny a) 1 -> get >>= \a -> get >>= \b -> get >>= \c -> return (ERFnInterval a b c) -- 1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> get >>= \e -> return (ERFnInterval a b c d e) _ -> fail "no parse" data ERFnContext = ERFnContext { erfnMaxDegree :: Int, erfnMaxSize :: Int, erfnCoeffGranularity :: Granularity } deriving (Show, Typeable, Data) instance Binary ERFnContext where put (ERFnContext a b c) = put a >> put b >> put c get = get >>= \a -> get >>= \b -> get >>= \c -> return (ERFnContext a b c) erfnContextDefault = ERFnContext { erfnMaxDegree = 2, erfnMaxSize = 20, erfnCoeffGranularity = 10 } erfnContextUnify (ERFnContext dg1 sz1 gr1) (ERFnContext dg2 sz2 gr2) = ERFnContext (max dg1 dg2) (max sz1 sz2) (max gr1 gr2) instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) => Show (ERFnInterval fb) where show (ERFnIntervalAny _) = "ERFnIntervalAny" show (ERFnInterval ln h ctxt) = "\nERFnInterval {" ++ show ctxt ++ "}" -- ++ " (definitely consistent: " ++ show isC -- ++ "anticonsistent: " ++ show isAC ++ ")" ++ "\n upper = " ++ ufbShow h ++ "\n lower = " ++ ufbShow (UFB.neg ln) -- ++ " global = " ++ show gl ++ "\n" -- ++ " context = " ++ show ctxt ++ "\n" where ufbShow = UFB.showDiGrCmp 10 False False instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) => H.HTML (ERFnInterval fb) where toHtml (ERFnIntervalAny ctxt) = H.toHtml "ERFnIntervalAny" toHtml (ERFnInterval ln h ctxt) = -- H.toHtml $ -- abovesTable -- [ -- H.toHtml "ERFnInterval", H.toHtml $ H.simpleTable [H.border 2] [] [ [H.toHtml "upper = ", H.toHtml $ ufbShow h], [H.toHtml "lower = ", H.toHtml $ ufbShow (UFB.neg ln)] ] -- ] where ufbShow = UFB.showDiGrCmp 10 False False instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) => Eq (ERFnInterval fb) where (ERFnInterval ln1 h1 ctxt1) == (ERFnInterval ln2 h2 ctxt2) = error "ERFnInterval: equality not implemented" _ == _ = error "ERFnInterval: equality not implemented" instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) => Ord (ERFnInterval fb) where compare (ERFnInterval ln1 h1 ctxt1) (ERFnInterval ln2 h2 ctxt2) = error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead" compare _ _ = error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead" instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => Num (ERFnInterval fb) where fromInteger n = UFA.const [fromInteger n] negate f@(ERFnIntervalAny _) = f negate (ERFnInterval ln h ctxt) = (ERFnInterval h ln ctxt) (ERFnInterval ln1 h1 ctxt1) + (ERFnInterval ln2 h2 ctxt2) = normalise $ ERFnInterval (reduceSzUp ln) (reduceSzUp h) ctxt where ln = ln1 +^ ln2 h = h1 +^ h2 reduceSzUp = UFB.reduceSizeUp maxSize maxSize = erfnMaxSize ctxt ctxt = erfnContextUnify ctxt1 ctxt2 f1 + f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) (ERFnInterval ln1 h1 ctxt1) * (ERFnInterval ln2 h2 ctxt2) = normalise $ ERFnInterval ln h ctxt where (ln, h) = multiplyEncl maxDegr maxSize (ln1, h1) (ln2, h2) maxDegr = erfnMaxDegree ctxt maxSize = erfnMaxSize ctxt ctxt = erfnContextUnify ctxt1 ctxt2 f1 * f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => Fractional (ERFnInterval fb) where fromRational r = UFA.const [fromRational r] recip f@(ERFnIntervalAny _) = f recip (ERFnInterval ln h ctxt) | certainNoZero = normalise $ ERFnInterval lnR hR ctxt | otherwise = ERFnIntervalAny ctxt where (hR, lnR) = UFB.recipEncl maxDegr maxSize ix (h,ln) certainNoZero = certainAboveZero || certainBelowZero certainAboveZero = UFB.upperBound ix ln < 0 certainBelowZero = UFB.upperBound ix h < 0 -- hnRecipUp = -- UFB.recipUp maxDegr maxSize ix (negate h) -- lRecipUp = -- UFB.recipUp maxDegr maxSize ix (negate ln) maxDegr = erfnMaxDegree ctxt maxSize = erfnMaxSize ctxt ix = int2effIx $ 3 * maxDegr instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => RA.ERApprox (ERFnInterval fb) where initialiseBaseArithmetic _ = UFB.initialiseBaseArithmetic (UFB.const 0 :: fb) getGranularity (ERFnIntervalAny ctxt) = erfnCoeffGranularity ctxt getGranularity (ERFnInterval ln h ctxt) = max (erfnCoeffGranularity ctxt) $ max (UFB.getGranularity ln) (UFB.getGranularity h) setGranularityOuter gran (ERFnIntervalAny ctxt) = ERFnIntervalAny $ ctxt { erfnCoeffGranularity = gran } setGranularityOuter gran (ERFnInterval ln h ctxt) = ERFnInterval (UFB.setGranularity gran ln) (UFB.setGranularity gran h) (ctxt { erfnCoeffGranularity = gran }) setMinGranularityOuter gran (ERFnIntervalAny ctxt) = ERFnIntervalAny (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) setMinGranularityOuter gran (ERFnInterval ln h ctxt) = ERFnInterval (UFB.setMinGranularity gran ln) (UFB.setMinGranularity gran h) (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) -- getPrecision (ERFnIntervalAny _) = 0 -- getPrecision f = intLog 2 (1 + (fst $ RA.integerBounds (FA.volume f))) -- wrong! f1@(ERFnInterval ln1 h1 ctxt1) /\ f2@(ERFnInterval ln2 h2 ctxt2) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: /\\:\n f1:\n" ++ show f1 ++ " f2:\n" ++ show f2 ++ "\n result:\n") $ ---- #endif normalise $ ERFnInterval (UFB.minUp maxDegr maxSize ln1 ln2) (UFB.minUp maxDegr maxSize h1 h2) ctxt where ctxt = erfnContextUnify ctxt1 ctxt2 maxDegr = erfnMaxDegree ctxt maxSize = erfnMaxSize ctxt (ERFnIntervalAny ctxt1) /\ (ERFnInterval ln2 h2 ctxt2) = ERFnInterval ln2 h2 ctxt where ctxt = erfnContextUnify ctxt1 ctxt2 (ERFnInterval ln1 h1 ctxt1) /\ (ERFnIntervalAny ctxt2) = ERFnInterval ln1 h1 ctxt where ctxt = erfnContextUnify ctxt1 ctxt2 f1 /\ f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) leqReals f1 f2 = -- unsafePrint ("ERInterval: leqReals: sizes: " ++ show (FA.getSize f1) ++ ", " ++ show (FA.getSize f2)) $ erfnintLeq f1 f2 refines _ (ERFnIntervalAny _) = True refines (ERFnIntervalAny _) _ = False refines (ERFnInterval ln1 h1 _) (ERFnInterval ln2 h2 _) = (UFB.upperBound 10 (ln2 -^ ln1) >= 0) && (UFB.upperBound 10 (h2 -^ h1) >= 0) compareApprox (ERFnIntervalAny _) (ERFnIntervalAny _) = EQ compareApprox (ERFnIntervalAny _) _ = LT compareApprox _ (ERFnIntervalAny _) = GT compareApprox (ERFnInterval ln1 h1 _) (ERFnInterval ln2 h2 _) = compareComposeMany [ UFB.compareApprox h1 h2, UFB.compareApprox ln1 ln2 ] erfnintLeq left right | left `isClearlyBelow` right = Just True | right `isClearlyStrictlyBelow` left = Just False | otherwise = Nothing where isClearlyBelow (ERFnIntervalAny _) _ = False isClearlyBelow _ (ERFnIntervalAny _) = False isClearlyBelow f g | UFB.upperBound 10 (erfnUpper f +^ erfnLowerNeg g) <= 0 = True | otherwise = False isClearlyStrictlyBelow (ERFnIntervalAny _) _ = False isClearlyStrictlyBelow _ (ERFnIntervalAny _) = False isClearlyStrictlyBelow f g | UFB.upperBound 10 (erfnUpper f +^ erfnLowerNeg g) < 0 = True | otherwise = False instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => RA.ERIntApprox (ERFnInterval fb) where -- doubleBounds = :: ira -> (Double, Double) -- floatBounds :: ira -> (Float, Float) -- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger) bisectDomain maybePt (ERFnIntervalAny c) = error "ERFnInterval: RA.bisectDomain: cannot bisect ERFnIntervalAny" bisectDomain maybePt (ERFnInterval ln h c) = (ERFnInterval ln midUp c, ERFnInterval midDownNeg h c) where (midDownNeg, midUp) = case maybePt of Nothing -> (UFB.scaleUp (1/2) $ ln -^ h, UFB.scaleUp (1/2) $ h -^ ln) Just (ERFnInterval lnPt hPt _) -> (lnPt, hPt) bounds (ERFnIntervalAny c) = error "ERFnInterval: RA.bounds: cannot get bounds for ERFnIntervalAny" bounds (ERFnInterval ln h c) = (ERFnInterval ln (UFB.neg ln) c, ERFnInterval (UFB.neg h) h c) f1@(ERFnInterval ln1 h1 c1) \/ f2@(ERFnInterval ln2 h2 c2) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: abs:\n f1:\n" ++ show f1 ++ " f2:\n" ++ show f2 ++ "\n result:\n") $ ---- #endif normalise $ ERFnInterval ln h c where ln = UFB.maxUp maxDegree maxSize ln1 ln2 h = UFB.maxUp maxDegree maxSize h1 h2 c = erfnContextUnify c1 c2 maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c (ERFnIntervalAny ctxt1) \/ (ERFnInterval ln2 h2 ctxt2) = ERFnIntervalAny ctxt where ctxt = erfnContextUnify ctxt1 ctxt2 (ERFnInterval ln1 h1 ctxt1) \/ (ERFnIntervalAny ctxt2) = ERFnIntervalAny ctxt where ctxt = erfnContextUnify ctxt1 ctxt2 f1 \/ f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, Show varid, Show boxra) => RAEL.ERApproxElementary (ERFnInterval fb) where -- default abs does not work because we do not have Prelude.abs abs _ f@(ERFnIntervalAny _) = f abs _ f@(ERFnInterval ln h c) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: abs:\n f:\n" ++ show f ++ "\n result:\n") $ ---- #endif normalise $ ERFnInterval minhln0Up maxhlnUp c where maxhlnUp = UFB.maxUp maxDegree maxSize h ln minhln0Up = UFB.minUp maxDegree maxSize (UFB.const 0) $ UFB.minUp maxDegree maxSize h ln maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c exp ix f@(ERFnIntervalAny _) = f exp ix f@(ERFnInterval ln h c) = normalise $ ERFnInterval lExpNeg hExp c where maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c (lExpNeg, hExp) = case (UFB.upperBound ix (h +^ ln) <= 1) of True -> UFB.expEncl maxDegree maxSize ix (ln, h) False -> (lExpNeg, hExp) where (lExpNeg, _) = UFB.expEncl maxDegree maxSize ix (ln, UFB.neg ln) (_, hExp) = UFB.expEncl maxDegree maxSize ix (UFB.neg h,h) sin ix f@(ERFnIntervalAny c) = ERFnInterval one one c where one = UFB.const 1 sin ix f@(ERFnInterval ln h c) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.sin: " -- ++ "\n h = " ++ show h -- ++ "\n ln = " ++ show ln -- ++ "\n hSin = " ++ show hSin -- ++ "\n lSinNeg = " ++ show lSinNeg -- ) $ ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: sin:\n f:\n" ++ show f ++ "\n result:\n") $ ---- #endif normalise $ ERFnInterval lSinNeg hSin c where (lSinNeg, hSin) = sincos True maxDegree maxSize ix (ln, h) maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c cos ix f@(ERFnIntervalAny c) = ERFnInterval one one c where one = UFB.const 1 cos ix f@(ERFnInterval ln h c) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.cos: " -- ++ "\n h = " ++ show h -- ++ "\n ln = " ++ show ln -- ++ "\n uCos = " ++ show uCos -- ++ "\n lCosNeg = " ++ show lCosNeg -- ) $ normalise $ ERFnInterval lCosNeg hCos c where (lCosNeg, hCos) = sincos False maxDegree maxSize ix (ln,h) maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c atan ix f@(ERFnIntervalAny c) = ERFnInterval one one c where one = UFB.const 1 atan ix f@(ERFnInterval ln h c) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.atan: " -- ++ "\n u = " ++ show u -- ++ "\n ln = " ++ show ln -- ++ "\n uAtan = " ++ show uAtan -- ++ "\n lAtanNeg = " ++ show lAtanNeg -- ) $ normalise $ ERFnInterval lAtanNeg hAtan c where maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c -- ix = int2effIx maxDegree (lAtanNeg, hAtan) = case (UFB.upperBound ix (h +^ ln) <= 1) of True -> UFB.atanEncl maxDegree maxSize ix (ln, h) False -> (lAtanNeg, hAtan) where (lAtanNeg, _) = UFB.atanEncl maxDegree maxSize ix (ln, UFB.neg ln) (_, hAtan) = UFB.atanEncl maxDegree maxSize ix (UFB.neg h,h) sincos :: (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, RAEL.ERApproxElementary ra) => Bool {-^ True iff sine, False iff cosine -} -> Int {-^ maximum representation degree -} -> Int {-^ maximum approx size -} -> EffortIndex {-^ how hard to try to eliminate truncation errors -} -> (fb, fb) -> (fb, fb) sincos isSine maxDegree maxSize ix (ln,h) -- p - 2k*pi range within [-pi/2, pi/2]: | ranfNear0 `RA.refines` plusMinusPiHalf = -- unsafePrint -- ( -- "ERFnInterval: sincos: [-pi/2, pi/2]: " -- ++ "\n u = " ++ show u -- ++ "\n l = " ++ show l -- ++ "\n ranf = " ++ show ranf -- ++ "\n k = " ++ show k -- ++ "\n ranfNear0 = " ++ show ranfNear0 -- ) $ case isSine of True -> sineShifted (- k2pi) False -> cosineShifted (- k2pi) -- p - 2k*pi range within [0, pi]: | (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf = -- unsafePrint -- ( -- "ERFnInterval: sincos: [0, pi]: " -- ++ "\n u = " ++ show u -- ++ "\n l = " ++ show l -- ++ "\n ranf = " ++ show ranf -- ++ "\n k = " ++ show k -- ++ "\n ranfNear0 = " ++ show ranfNear0 -- ) $ case isSine of -- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2): True -> cosineShifted (- k2pi - piHalf) False -> sineShiftedNegated (- k2pi - piHalf) -- p - 2k*pi range within [-pi, 0]: | (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf = -- unsafePrint -- ( -- "ERFnInterval: sincos: [-pi, 0]: " -- ++ "\n u = " ++ show u -- ++ "\n l = " ++ show l -- ++ "\n ranf = " ++ show ranf -- ++ "\n k = " ++ show k -- ++ "\n ranfNear0 = " ++ show ranfNear0 -- ) $ case isSine of -- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2): True -> cosineShiftedNegated (-k2pi + piHalf) False -> sineShifted (-k2pi + piHalf) -- p - 2k*pi range within [pi/2, 3pi/2]: | (ranfNear0 - pi) `RA.refines` plusMinusPiHalf = -- unsafePrint -- ( -- "ERFnInterval: sincos: [pi/2, 3pi/2]: " -- ++ "\n u = " ++ show u -- ++ "\n l = " ++ show l -- ++ "\n ranf = " ++ show ranf -- ++ "\n k = " ++ show k -- ++ "\n ranfNear0 = " ++ show ranfNear0 -- ) $ -- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi) case isSine of True -> sineShiftedNegated (- k2pi - pi) False -> cosineShiftedNegated (- k2pi - pi) | otherwise = -- unsafePrint -- ( -- "ERFnInterval: sincos: big range: " -- ++ "\n u = " ++ show u -- ++ "\n l = " ++ show l -- ++ "\n ranf = " ++ show ranf -- ++ "\n k = " ++ show k -- ++ "\n ranfNear0 = " ++ show ranfNear0 -- ) $ (UFB.const (-1), UFB.const 1) -- (expDownwards, expUpwards + valueAtRDnNeg + (UFB.const expRUp)) where -- l = UFB.neg ln ranfNear0 = ranf - k2pi k2pi = k * 2 * pi plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO pi = RAEL.pi ix piHalf = pi / 2 (piHalfLO, piHalfHI) = RA.bounds piHalf ranf = ERInterval (negate $ UFB.upperBound 10 ln) (UFB.upperBound 10 h) k = fromInteger $ toInteger kEI (kEI,_) = RA.integerBounds $ 0.5 + (ranf / (2*pi)) sineShiftedNegated shift = boundsNegate $ sineShifted shift cosineShiftedNegated shift = boundsNegate $ cosineShifted shift boundsNegate (pLONeg, pHI) = (pHI, pLONeg) sineShifted shift = -- moving to domain where sinus is non-decreasing case (UFB.upperBound ix (h +^ ln) <= 0.25) of True -> UFB.sinEncl maxDegree maxSize ix (lnShifted, hShifted) False -> (lSinNeg, hSin) where (lSinNeg, _) = UFB.sinEncl maxDegree maxSize ix (ln, UFB.neg ln) (_, hSin) = UFB.sinEncl maxDegree maxSize ix (UFB.neg h,h) where lnShifted = ln +^ (UFB.const (- shiftLOB)) hShifted = h +^ (UFB.const shiftHIB) ERInterval shiftLOB shiftHIB = shift cosineShifted shift = -- moving to domain where cosinus is non-decreasing case (UFB.upperBound ix (h +^ ln) <= 0.25) of True -> UFB.cosEncl maxDegree maxSize ix (lnShifted, hShifted) False -> (UFB.minUp maxDegree maxSize lCosDownNeg hCosDownNeg, UFB.maxUp maxDegree maxSize lCosUp hCosUp +^ (UFB.scaleUp 0.5 (h +^ ln))) -- this term is important when enclosure hits 0; -- without it, the result could miss cosine's maximum at 0 where (lCosDownNeg, lCosUp) = UFB.cosEncl maxDegree maxSize ix (ln, UFB.neg ln) (hCosDownNeg, hCosUp) = UFB.cosEncl maxDegree maxSize ix (UFB.neg h,h) lnShifted = ln +^ (UFB.const (- shiftLOB)) hShifted = h +^ (UFB.const shiftHIB) ERInterval shiftLOB shiftHIB = shift boundsAddErr errB (pLONeg, pHI) = (pLONeg +^ errPoly, pHI +^ errPoly) where errPoly = UFB.const errB normalise f@(ERFnIntervalAny c) = f normalise f@(ERFnInterval ln h c) | UFB.isValid h && UFB.isValid ln = f | otherwise = ERFnIntervalAny c check callerLocation f@(ERFnIntervalAny c) = f check callerLocation f@(ERFnInterval ln h c) = ERFnInterval (UFB.check (callerLocation ++ "upper: ") h) (UFB.check (callerLocation ++ "neg lower: ") ln) c instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => FA.ERFnApprox boxra varid ra ra (ERFnInterval fb) where check = check domra2ranra _ = id ranra2domra _ = id getMaxDegree (ERFnIntervalAny c) = erfnMaxDegree c getMaxDegree (ERFnInterval _ _ c) = erfnMaxDegree c setMaxDegree maxDegr (ERFnIntervalAny c) = ERFnIntervalAny (c { erfnMaxDegree = maxDegr } ) setMaxDegree maxDegr (ERFnInterval ln h c) = ERFnInterval (UFB.reduceDegreeUp maxDegr ln) (UFB.reduceDegreeUp maxDegr h) (c { erfnMaxDegree = maxDegr } ) getSize (ERFnIntervalAny c) = 0 getSize (ERFnInterval ln h c) = max (UFB.getSize ln) (UFB.getSize h) getMaxSize (ERFnIntervalAny c) = erfnMaxSize c getMaxSize (ERFnInterval _ _ c) = erfnMaxSize c setMaxSize maxSize (ERFnIntervalAny c) = ERFnIntervalAny (c { erfnMaxDegree = maxSize } ) setMaxSize maxSize (ERFnInterval ln h c) = ERFnInterval (UFB.reduceSizeUp maxSize ln) (UFB.reduceSizeUp maxSize h) (c { erfnMaxSize = maxSize } ) getVariables (ERFnIntervalAny _) = [] getVariables (ERFnInterval ln h _) = UFB.getVariables h getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox getRangeApprox (ERFnInterval ln h c) = UFB.raFromEndpoints h ( (- (UFB.upperBound 10 ln)) , (UFB.upperBound 10 h) ) scale ratio f@(ERFnIntervalAny c) = f scale ratio f@(ERFnInterval ln h c) = ---- #ifdef RUNTIME_CHECKS ---- FA.check ("ERFnInterval: scale:\n before:\n" ++ show f ++ "\n after:\n") $ ---- #endif normalise $ case RA.compareReals ratio 0 of Just GT -> ERFnInterval (scaleUp ratio ln) (scaleUp ratio h) c Just LT -> ERFnInterval (scaleUp (- ratio) h) (scaleUp (- ratio) ln) c _ -> (UFA.const [ratio]) * f where scaleUp = UFB.scaleApproxUp maxDegree maxSize maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c eval ptBox (ERFnIntervalAny c) = [RA.bottomApprox] eval ptBox (ERFnInterval ln h c) = [lo RA.\/ up] where up = UFB.evalApprox ptBox h lo = negate $ UFB.evalApprox ptBox ln partialEval substitutions f@(ERFnIntervalAny c) = f partialEval substitutions f@(ERFnInterval ln h c) = normalise $ ERFnInterval lnP hP c where hP = UFB.partialEvalApproxUp substitutions h lnP = UFB.partialEvalApproxUp substitutions ln composeNonDecreasing fOuter@(ERFnInterval lnOuter hOuter cOuter) varid fInner@(ERFnInterval lnInner hInner cInner) = -- unsafePrintReturn -- ( -- "ER.RnToRm.UnitDom.Interval: composeNonDecreasing: " -- ++ "\n fOuter = " ++ show fOuter -- ++ "\n varid = " ++ show varid -- ++ "\n fInner = " ++ show fInner -- ++ "\n inconsistencies = " ++ show (UFA.keyPointsConsistencyCheck resultReals result) -- ++ "\n result = " -- ) -- $ result where resultReals ptB = -- this is only used for consistency checking... (\[x] -> x) $ FA.eval ptBOuter fOuter where ptBOuter = DBox.insert varid fInnerVal ptB fInnerVal = FA.ranra2domra fInner $ (\[x] -> x) $ FA.eval ptB fInner result = ERFnInterval ln h c h = erfnUpper $ UFA.composeWithThin fOuter $ Map.singleton varid (ERFnInterval (UFB.neg hInner) hInner cInner) ln = erfnLowerNeg $ UFA.composeWithThin fOuter $ Map.singleton varid $ (ERFnInterval lnInner (UFB.neg lnInner) cInner) c = erfnContextUnify cOuter cInner composeNonDecreasing fOuter varid fInner = ERFnIntervalAny c where c = erfnContextUnify (erfnContext fOuter) (erfnContext fInner) instance (UFB.ERUnitFnBaseElementary boxb boxra varid b ra fb, Show varid, Show boxra) => UFA.ERUnitFnApprox boxra varid ra ra (ERFnInterval fb) where bottomApprox = ERFnIntervalAny erfnContextDefault const [val] | RA.isBounded val = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: const:\n") $ ---- #endif normalise $ ERFnInterval { erfnLowerNeg = fbLNeg, erfnUpper = fbH, erfnContext = context } | otherwise = ERFnIntervalAny context where fbH = UFB.const valH fbLNeg = UFB.const (negate valL) (valL, valH) = UFB.raEndpoints fbH val context = erfnContextDefault { erfnCoeffGranularity = RA.getGranularity val } affine [val] coeffsSingletons | RA.isBounded val && (and $ map (RA.isBounded . head) $ Map.elems coeffsSingletons) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: affine:\n") $ ---- #endif normalise $ ERFnInterval { erfnLowerNeg = fbLNeg, erfnUpper = fbH, erfnContext = context -- , -- erfnGlobal = -- UFB.raFromEndpoints fbH -- (valL - coeffCorr - coeffsAbsSum, -- valH + coeffCorr + coeffsAbsSum) } | otherwise = ERFnIntervalAny context where coeffs = Map.map (\[a] -> a) coeffsSingletons coeffGranularity = Map.fold max (RA.getGranularity val) (Map.map RA.getGranularity coeffs) coeffsMsCorrs = Map.map (\(l,h) -> (B.setMinGranularity coeffGranularity (l + h)/2, B.setMinGranularity coeffGranularity (h - l)/2)) $ Map.map (UFB.raEndpoints fbH) $ coeffs coeffCorr = Map.fold (+) 0 $ Map.map snd coeffsMsCorrs coeffsAbsSum = Map.fold (+) 0 $ Map.map (abs . fst) coeffsMsCorrs fbH = UFB.affine (valH + coeffCorr) (Map.map fst coeffsMsCorrs) fbLNeg = UFB.affine (negate (valL - coeffCorr)) (Map.map (negate . fst) coeffsMsCorrs) (valL, valH) = UFB.raEndpoints fbH val context = erfnContextDefault { erfnCoeffGranularity = coeffGranularity } composeWithThin f@(ERFnIntervalAny ctxt) substitutions = f composeWithThin f@(ERFnInterval ln1 h1 ctxt1) substitutions = -- unsafePrintReturn -- ( -- "ER.RnToRm.UnitDom.Interval: composeWithThin: " -- ++ "\n f = " ++ show f -- ++ "\n substitutions = " ++ show substitutions -- ++ "\n inconsistencies = " ++ show (UFA.keyPointsConsistencyCheck resultReals result) -- ++ "\n result = " -- ) -- $ result where resultReals ptB = -- this is only used for consistency checking... (\[x] -> x) $ FA.eval ptBOuter f where ptBOuter = foldl insertVal ptB $ Map.toList substitutions insertVal ptB (varid, fInner) = DBox.insert varid (evalPtB fInner) ptB evalPtB fInner = FA.ranra2domra fInner $ (\[x] -> x) $ FA.eval ptB fInner result = ERFnInterval ln h ctxt1 ln = UFB.composeManyUp maxDegree maxSize ln1 ufbSubstitutions h = UFB.composeManyUp maxDegree maxSize h1 ufbSubstitutions ufbSubstitutions = Map.map erfnUpper substitutions maxDegree = erfnMaxDegree ctxt1 maxSize = erfnMaxSize ctxt1 -- ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnIntervalAny ctxt1) f2@(ERFnIntervalAny ctxt2) = (ERFnIntervalAny ctxt, RA.bottomApprox) where ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnIntervalAny ctxt1) f2@(ERFnInterval ln2 h2 ctxt2) = (ERFnInterval ln2 h2 ctxt, RA.plusInfinity) where ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnInterval ln1 h1 ctxt1) f2@(ERFnIntervalAny ctxt2) = (ERFnInterval ln1 h1 ctxt, 1) where ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnInterval ln1 h1 ctxt1) f2@(ERFnInterval ln2 h2 ctxt2) = case RA.compareReals improvementRA 1 of Just LT -> (f1, 1) -- intersection made it worse, keep original _ -> (intersection, improvementRA) where intersection = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: intersectMeasureImprovement:\n f1:\n" ++ show f1 ++ "\n f2:\n" ++ show f2 ++ "\n intersection:\n") $ ---- #endif normalise $ f1 RA./\ f2 improvementRA | 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1 -- error $ -- "ERFnInterval: intersectMeasureImprovement: inconsistent result: " -- ++ show intersection | otherwise = f1Volume / intersectionVolume intersectionVolume = UFA.volume vars intersection f1Volume = UFA.volume vars f1 ctxt = erfnContextUnify ctxt1 ctxt2 volume vars (ERFnIntervalAny c) = RA.plusInfinity volume vars (ERFnInterval ln h c) = UFB.raFromEndpoints h (volL, volH) where volH = UFB.volumeAboveZeroUp vars (ln +^ h) volL = negate $ UFB.volumeAboveZeroUp vars (l +^ hn) l = UFB.neg ln hn = UFB.neg h integrate _ f@(ERFnIntervalAny c) _ _ _ = f integrate ix fD@(ERFnInterval ln h c) x origin fI@(ERFnInterval lnInit hInit cInit) = -- unsafePrintReturn -- ( -- "ERFnInterval: integrate: " -- ++ "\n u = " ++ show u -- ++ "\n ln = " ++ show ln -- ++ "\n origin = " ++ show origin -- ++ "\n uInit = " ++ show uInit -- ++ "\n lnInit = " ++ show lnInit -- ++ "\n uIuL = " ++ show uIuL -- ++ "\n uIuU = " ++ show uIuU -- ++ "\n uIuOriginL = " ++ show uIuOriginL -- ++ "\n uIuOriginU = " ++ show uIuOriginU -- ++ "\n lnIuL = " ++ show lnIuL -- ++ "\n lnIuU = " ++ show lnIuU -- ++ "\n lnIuOriginL = " ++ show lnIuOriginL -- ++ "\n lnIuOriginU = " ++ show lnIuOriginU -- ++ "\n uIov = " ++ show uIov -- ++ "\n lnIov = " ++ show lnIov -- ++ "\n result = " -- ) -- $ ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: integrate:\n fD:\n" ++ show fD ++ "\n fI:\n" ++ show fI ++ "\n result:\n") $ ---- #endif normalise $ (ERFnInterval lnIov hIov c) where -- perform raw integration of both bounds: (hIuL, hIuH) = -- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ UFB.integrate x h (lnIuL, lnIuH) = -- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ UFB.integrate x ln maxDegree = erfnMaxDegree c maxSize = erfnMaxSize c -- constrain the raw integrals to the origin: (hIuOriginLNeg, hIuOriginH) = UFB.composeEncl maxDegree maxSize hIuL x originEncl (lnIuOriginLNeg, lnIuOriginH) = UFB.composeEncl maxDegree maxSize lnIuL x originEncl originEncl = UFB.constEncl $ UFB.raEndpoints h origin -- adjust the raw integrated functions to enclose the initial condition function: hIov = UFB.reduceSizeUp maxSize $ hIuH +^ hInit +^ hIuOriginLNeg +^ (hIuOriginH +^ hIuOriginLNeg) lnIov = UFB.reduceSizeUp maxSize $ lnIuH +^ lnInit +^ lnIuOriginLNeg +^ (lnIuOriginH +^ lnIuOriginLNeg)