module Data.Number.ER.RnToRm.UnitDom.Approx.Interval
(
ERFnInterval(..),
ERFnContext(..)
)
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 qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Approx.Elementary as RAEL
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc
import qualified Data.Map as Map
import Data.Typeable
import Data.Generics.Basics
import Data.Binary
data ERFnInterval fb ra =
ERFnIntervalAny
{
erfnContext :: ERFnContext
}
|
ERFnInterval
{
erfnUpper :: fb,
erfnLowerNeg :: fb,
erfnContext :: ERFnContext,
erfnGlobal :: ra
}
deriving (Typeable, Data)
instance (Binary a, Binary b) => Binary (ERFnInterval a b) where
put (ERFnIntervalAny a) = putWord8 0 >> put a
put (ERFnInterval a b c d) = putWord8 1 >> put a >> put b >> put c >> put d
get = do
tag_ <- getWord8
case tag_ of
0 -> get >>= \a -> return (ERFnIntervalAny a)
1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFnInterval a b c d)
_ -> fail "no parse"
data ERFnContext =
ERFnContext
{
erfnMaxDegree :: Int,
erfnCoeffGranularity :: Granularity
}
deriving (Show, Typeable, Data)
instance Binary ERFnContext where
put (ERFnContext a b) = put a >> put b
get = get >>= \a -> get >>= \b -> return (ERFnContext a b)
erfnContextDefault =
ERFnContext
{
erfnMaxDegree = 2,
erfnCoeffGranularity = 10
}
erfnContextUnify (ERFnContext dg1 gr1) (ERFnContext dg2 gr2) =
ERFnContext (max dg1 dg2) (max gr1 gr2)
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
Show (ERFnInterval fb ra)
where
show (ERFnIntervalAny _) = "ERFnIntervalAny"
show (ERFnInterval h ln ctxt gl) =
"\nERFnInterval"
++ "\n upper = " ++ show h
++ "\n lower = " ++ show (ln)
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
Eq (ERFnInterval fb ra)
where
(ERFnInterval h1 ln1 ctxt1 gl1)
== (ERFnInterval h2 ln2 ctxt2 gl2) =
error "ERFnInterval: equality not implemented yet"
_ == _ =
error "ERFnInterval: equality not implemented yet"
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
Ord (ERFnInterval fb ra)
where
compare
(ERFnInterval h1 ln1 ctxt1 gl1)
(ERFnInterval h2 ln2 ctxt2 gl2) =
error "ERFnInterval: comparison not implemented yet"
compare _ _ =
error "ERFnInterval: comparison not implemented yet"
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
Num (ERFnInterval fb ra)
where
fromInteger n = UFA.const [fromInteger n]
negate f@(ERFnIntervalAny _) = f
negate (ERFnInterval h ln ctxt gl) =
(ERFnInterval ln h ctxt (negate gl))
(ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =
ERFnInterval (h1 + h2) (ln1 + ln2) ctxt (gl1 + gl2)
where
ctxt = erfnContextUnify ctxt1 ctxt2
f1 + f2 = ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
(ERFnInterval h1 ln1 ctxt1 gl1) * (ERFnInterval h2 ln2 ctxt2 gl2) =
ERFnInterval h ln ctxt (gl1 * gl2)
where
(h, ln) =
case (RA.leqReals 0 gl1, RA.leqReals gl1 0, RA.leqReals 0 gl2, RA.leqReals gl2 0) of
(Just True, _, Just True, _) ->
(h1h2, l1l2Neg)
(_, Just True, _, Just True) ->
(l1l2, h1h2Neg)
(Just True, _, _, Just True) ->
(l1h2, h1l2Neg)
(_, Just True, Just True, _) ->
(h1l2, l1h2Neg)
_ ->
((h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2),
(h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg))
where
h1h2 = UFB.reduceDegreeUp maxDegr $ h1 * h2
h1h2Neg = UFB.reduceDegreeUp maxDegr $ (negate h1) * h2
l1l2 = UFB.reduceDegreeUp maxDegr $ ln1 * ln2
l1l2Neg = UFB.reduceDegreeUp maxDegr $ (negate ln1) * ln2
h1l2 = UFB.reduceDegreeUp maxDegr $ h1 * (negate ln2)
h1l2Neg = UFB.reduceDegreeUp maxDegr $ h1 * ln2
l1h2 = UFB.reduceDegreeUp maxDegr $ (negate ln1) * h2
l1h2Neg = UFB.reduceDegreeUp maxDegr $ ln1 * h2
maxP p1 p2 = fst $ UFB.max maxDegr p1 p2
ctxt = erfnContextUnify ctxt1 ctxt2
maxDegr = erfnMaxDegree ctxt
f1 * f2 = ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
Fractional (ERFnInterval fb ra)
where
fromRational r = UFA.const [fromRational r]
recip f@(ERFnIntervalAny _) = f
recip (ERFnInterval h ln ctxt gl)
| certainNoZero =
ERFnInterval lRecipUp hnRecipUp ctxt (recip gl)
| otherwise = ERFnIntervalAny ctxt
where
certainNoZero =
certainAboveZero || certainBelowZero
certainAboveZero =
UFB.upperBound ix ln < 0
certainBelowZero =
UFB.upperBound ix h < 0
hnRecipUp =
UFB.recipUp maxDegr ix (negate h)
lRecipUp =
UFB.recipUp maxDegr ix (negate ln)
maxDegr = erfnMaxDegree ctxt
ix = int2effIx $ 3 * maxDegr
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
RA.ERApprox (ERFnInterval fb ra)
where
getGranularity (ERFnIntervalAny ctxt) = erfnCoeffGranularity ctxt
getGranularity (ERFnInterval h ln ctxt gl) =
max (erfnCoeffGranularity ctxt) $
max (UFB.getGranularity h) (UFB.getGranularity ln)
setGranularity gran (ERFnIntervalAny ctxt) =
ERFnIntervalAny $ ctxt { erfnCoeffGranularity = gran }
setGranularity gran (ERFnInterval h ln ctxt gl) =
ERFnInterval
(UFB.setGranularity gran h) (UFB.setGranularity gran ln)
(ctxt { erfnCoeffGranularity = gran }) gl
setMinGranularity gran (ERFnIntervalAny ctxt) =
ERFnIntervalAny
(ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) })
setMinGranularity gran (ERFnInterval h ln ctxt gl) =
ERFnInterval
(UFB.setMinGranularity gran h) (UFB.setMinGranularity gran ln)
(ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) gl
(ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =
ERFnInterval (snd $ UFB.min maxDegr h1 h2) (snd $ UFB.min maxDegr ln1 ln2) ctxt (gl1 RA./\ gl2)
where
ctxt = erfnContextUnify ctxt1 ctxt2
maxDegr = erfnMaxDegree ctxt
(ERFnIntervalAny ctxt1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =
ERFnInterval h2 ln2 ctxt gl2
where
ctxt = erfnContextUnify ctxt1 ctxt2
(ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnIntervalAny ctxt2) =
ERFnInterval h1 ln1 ctxt gl1
where
ctxt = erfnContextUnify ctxt1 ctxt2
f1 /\ f2 = ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
leqReals = erfnintLeq
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.ERUnitFnBase boxb boxra varid b ra fb) =>
RA.ERIntApprox (ERFnInterval fb ra)
where
bisectDomain maybePt (ERFnIntervalAny c) =
error "ERFnInterval: RA.bisectDomain: cannot bisect ERFnIntervalAny"
bisectDomain maybePt (ERFnInterval u ln c g) =
(ERFnInterval midUp ln c g,
ERFnInterval u (negate midDown) c g)
where
(midDown, midUp) =
case maybePt of
Nothing ->
(negate $ (ln u) / 2, (u ln) / 2)
Just (ERFnInterval uPt lnPt _ _) ->
(negate lnPt, uPt)
bounds (ERFnIntervalAny c) =
error "ERFnInterval: RA.bounds: cannot get bounds for ERFnIntervalAny"
bounds (ERFnInterval u ln c g) =
(ERFnInterval (negate ln) ln c g,
ERFnInterval u (negate u) c g)
(ERFnInterval u1 ln1 c1 g1) \/ (ERFnInterval u2 ln2 c2 g2) =
ERFnInterval u ln c (g1 RA.\/ g2)
where
u = UFB.maxUp maxDegree u1 u2
ln = UFB.maxUp maxDegree ln1 ln2
c = erfnContextUnify c1 c2
maxDegree = erfnMaxDegree c
(ERFnIntervalAny ctxt1) \/ (ERFnInterval h2 ln2 ctxt2 gl2) =
ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify ctxt1 ctxt2
(ERFnInterval h1 ln1 ctxt1 gl1) \/ (ERFnIntervalAny ctxt2) =
ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify ctxt1 ctxt2
f1 \/ f2 = ERFnIntervalAny ctxt
where
ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra) =>
RAEL.ERApproxElementary (ERFnInterval fb ra)
where
abs _ f@(ERFnIntervalAny _) = f
abs _ (ERFnInterval u ln c g) =
ERFnInterval maxulnUp maxunl0Dn c (abs g)
where
maxDegree = erfnMaxDegree c
maxulnUp = snd $ UFB.max maxDegree u ln
maxunl0Dn =
fst $ UFB.max maxDegree 0 $
fst $ UFB.max maxDegree ( u) ( ln)
exp ix f@(ERFnIntervalAny _) = f
exp ix (ERFnInterval u ln c g) =
ERFnInterval uExp lExpNeg c (RAEL.exp ix g)
where
maxDegree = erfnMaxDegree c
uExp = snd $ UFB.exp maxDegree ix u
lExpNeg =
negate $ fst $ UFB.exp maxDegree ix (negate ln)
sin ix f@(ERFnIntervalAny c) =
ERFnInterval 1 1 c ((1) RA.\/ 1)
sin ix (ERFnInterval u ln c g) =
ERFnInterval uSin lSinNeg c (RAEL.sin ix g)
where
maxDegree = erfnMaxDegree c
uSin = snd $ UFB.sin maxDegree ix u
lSinNeg =
negate $ fst $ UFB.sin maxDegree ix (negate ln)
cos ix f@(ERFnIntervalAny c) =
ERFnInterval 1 1 c ((1) RA.\/ 1)
cos ix (ERFnInterval u ln c g) =
ERFnInterval uCos lCosNeg c (RAEL.cos ix g)
where
maxDegree = erfnMaxDegree c
uCos = snd $ UFB.cos maxDegree ix u
lCosNeg =
negate $ fst $ UFB.cos maxDegree ix (negate ln)
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
FA.ERFnApprox boxra varid ra ra (ERFnInterval fb ra)
where
check callerLocation f@(ERFnIntervalAny c) = f
check callerLocation (ERFnInterval u ln c g) =
ERFnInterval
(UFB.check (callerLocation ++ "upper: ") u)
(UFB.check (callerLocation ++ "neg lower: ") ln)
c g
domra2ranra _ = id
ranra2domra _ = id
setMaxDegree maxDegr (ERFnIntervalAny c) =
ERFnIntervalAny (c { erfnMaxDegree = maxDegr } )
setMaxDegree maxDegr (ERFnInterval u ln c g) =
ERFnInterval
(UFB.reduceDegreeUp maxDegr u)
(UFB.reduceDegreeUp maxDegr ln)
(c { erfnMaxDegree = maxDegr } )
g
getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox
getRangeApprox (ERFnInterval u ln c g) =
UFB.raFromEndpoints u
(
(UFB.upperBound 10 u)
,
( (UFB.upperBound 10 ln))
)
scale ratio f@(ERFnIntervalAny c) = f
scale ratio f@(ERFnInterval u ln c g) =
case RA.compareReals ratio 0 of
Just GT ->
ERFnInterval (UFB.scaleApproxUp ratio u) (UFB.scaleApproxUp ratio ln) c g
Just LT ->
ERFnInterval (UFB.scaleApproxUp ( ratio) ln) (UFB.scaleApproxUp ( ratio) u) c g
_ ->
(UFA.const [ratio]) * f
eval ptBox (ERFnIntervalAny c) = [RA.bottomApprox]
eval ptBox (ERFnInterval u ln c g) =
[lo RA.\/ up]
where
up = UFB.evalApprox ptBox u
lo = negate $ UFB.evalApprox ptBox ln
partialEval substitutions f@(ERFnIntervalAny c) = f
partialEval substitutions (ERFnInterval u ln c g) =
(ERFnInterval uP lnP c g)
where
uP = UFB.partialEvalApproxUp substitutions u
lnP = UFB.partialEvalApproxUp substitutions ln
composeThin
f@(ERFnIntervalAny ctxt)
substitutions =
f
composeThin
f@(ERFnInterval h1 ln1 ctxt1 gl1)
substitutions =
(ERFnInterval h ln ctxt1 gl1)
where
h = UFB.composeUp maxDegree h1 ufbSubstitutions
ln = UFB.composeUp maxDegree ln1 ufbSubstitutions
ufbSubstitutions = Map.map erfnUpper substitutions
maxDegree = erfnMaxDegree ctxt1
instance
(UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
UFA.ERUnitFnApprox boxra varid ra ra (ERFnInterval fb ra)
where
bottomApprox =
ERFnIntervalAny erfnContextDefault
const [val]
| RA.isBounded val =
ERFnInterval
{
erfnUpper = fbH,
erfnLowerNeg = fbLNeg,
erfnContext = context,
erfnGlobal = val
}
| 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) =
ERFnInterval
{
erfnUpper = fbH,
erfnLowerNeg = fbLNeg,
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
}
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 h2 ln2 ctxt2 gl2) =
(ERFnInterval h2 ln2 ctxt gl2, 1 / 0)
where
ctxt = erfnContextUnify ctxt1 ctxt2
intersectMeasureImprovement ix vars
f1@(ERFnInterval h1 ln1 ctxt1 gl1)
f2@(ERFnIntervalAny ctxt2) =
(ERFnInterval h1 ln1 ctxt gl1, 1)
where
ctxt = erfnContextUnify ctxt1 ctxt2
intersectMeasureImprovement ix vars
f1@(ERFnInterval h1 ln1 ctxt1 gl1)
f2@(ERFnInterval h2 ln2 ctxt2 gl2) =
case RA.compareReals improvementRA 1 of
Just LT -> (f1, 1)
_ -> (intersection, improvementRA)
where
intersection = f1 RA./\ f2
improvementRA
| 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1
| otherwise =
f1Volume / intersectionVolume
intersectionVolume = UFA.volume vars intersection
f1Volume = UFA.volume vars f1
ctxt = erfnContextUnify ctxt1 ctxt2
volume vars (ERFnIntervalAny c) = 1/0
volume vars (ERFnInterval u ln c g) =
UFB.raFromEndpoints u $ UFB.volumeAboveZero vars (u + ln)
integrate
ix (ERFnInterval u ln c g) x
origin (ERFnInterval uInit lnInit cInit gInit) =
(ERFnInterval uIov lnIov c gIov)
where
(uIuL, uIuU) =
UFB.integrate x u
(lnIuL, lnIuU) =
UFB.integrate x ln
maxDegree = erfnMaxDegree c
uIuOriginL = UFB.composeDown maxDegree uIuL substXOrigin
uIuOriginU = UFB.composeUp maxDegree uIuU substXOrigin
lnIuOriginL = UFB.composeDown maxDegree lnIuL substXOrigin
lnIuOriginU = UFB.composeUp maxDegree lnIuU substXOrigin
substXOrigin = Map.singleton x originUFB
originUFB = UFB.const $ fst $ UFB.raEndpoints u origin
uIov =
UFB.reduceDegreeUp maxDegree $
uIuU + uInit uIuOriginL + (uIuOriginU uIuOriginL)
lnIov =
UFB.reduceDegreeUp maxDegree $
lnIuU + lnInit lnIuOriginL + (lnIuOriginU lnIuOriginL)
gIov =
gInit + g * ((1 origin) RA.\/ (1 origin))