AERN-Real-0.10.0: src/Data/Number/ER/Real/Approx/Interval.hs
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleInstances #-}
{-|
Module : Data.Number.ER.Real.Approx.Interval
Description : safe interval arithmetic
Copyright : (c) Michal Konecny
License : BSD3
Maintainer : mik@konecny.aow.cz
Stability : experimental
Portability : portable
This module defines an arbitrary precision interval type and
most of its interval arithmetic operations.
-}
module Data.Number.ER.Real.Approx.Interval
(
ERInterval(..),
normaliseERIntervalOuter,
normaliseERIntervalInner
)
where
import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:))
import qualified Data.Number.ER.Real.Approx.Elementary as RAEL
import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc
import Data.Ratio
import qualified Text.Html as H
import Data.Typeable
import Data.Generics.Basics
import Data.Binary
--import BinaryDerive
{-|
Type for arbitrary precision interval arithmetic.
-}
data ERInterval base =
-- ERIntervalEmpty -- ^ usually represents computation error (top element in the interval domain)
-- | ERIntervalAny -- ^ represents no knowledge of result (bottom element in the interval domain)
ERInterval
{
erintv_left :: !base,
erintv_right :: !base
}
deriving (Typeable, Data)
{- the following has been generated by BinaryDerive -}
instance (Binary a) => Binary (ERInterval a) where
put (ERInterval a b) = putWord8 0 >> put a >> put b
get = do
tag_ <- getWord8
case tag_ of
0 -> get >>= \a -> get >>= \b -> return (ERInterval a b)
_ -> fail "no parse"
{- the above has been generated by BinaryDerive -}
{-|
convert to a normal form, ie:
* no NaNs as endpoints
Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic.
This version interprets NaN's as bottomApprox.
-}
normaliseERIntervalOuter ::
(B.ERRealBase b) =>
ERInterval b -> ERInterval b
normaliseERIntervalOuter (ERInterval nan1 nan2)
| B.isERNaN nan1 && B.isERNaN nan2 =
RA.bottomApprox
normaliseERIntervalOuter (ERInterval nan r)
| B.isERNaN nan =
ERInterval (- B.plusInfinity) r
normaliseERIntervalOuter (ERInterval l nan)
| B.isERNaN nan =
ERInterval l (B.plusInfinity)
normaliseERIntervalOuter i = i
{-|
convert to a normal form, ie:
* no NaNs as endpoints
Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic.
This version interprets NaN's as topApprox.
-}
normaliseERIntervalInner ::
(B.ERRealBase b) =>
ERInterval b -> ERInterval b
normaliseERIntervalInner (ERInterval nan1 nan2)
| B.isERNaN nan1 && B.isERNaN nan2 =
RA.topApprox
normaliseERIntervalInner (ERInterval nan r)
| B.isERNaN nan =
ERInterval (B.plusInfinity) r
normaliseERIntervalInner (ERInterval l nan)
| B.isERNaN nan =
ERInterval l (- B.plusInfinity)
normaliseERIntervalInner i = i
{-|
erintvPrecision returns an approximation of the number of bits required
to represent the mantissa of a normalised size of the interval:
> - log_2 ((r - l) / (1 + abs(r) + abs(l)))
Notice that this is +Infty for singleton and anti-consistent intervals
and -Infty for unbounded intervals.
-}
erintvPrecision ::
(B.ERRealBase b) =>
ERInterval b -> EI.ExtendedInteger
erintvPrecision i@(ERInterval l r)
| not $ RA.isConsistent i = EI.PlusInfinity
| not $ RA.isBounded i = EI.MinusInfinity
| otherwise =
-1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l))
erintvGranularity ::
(B.ERRealBase b) =>
ERInterval b -> Int
erintvGranularity (ERInterval l r) =
min (B.getGranularity l) (B.getGranularity r)
{- syntactic comparisons -}
{-|
a syntactic equality test
-}
erintvEqualApprox ::
(B.ERRealBase b) =>
ERInterval b -> ERInterval b -> Bool
erintvEqualApprox (ERInterval l1 r1) (ERInterval l2 r2) =
l1 == l2 && r1 == r2
{-|
a syntactic linear order
-}
erintvCompareApprox ::
(B.ERRealBase b) =>
ERInterval b -> ERInterval b -> Ordering
erintvCompareApprox (ERInterval l1 r1) (ERInterval l2 r2) =
case compare l1 l2 of
EQ -> compare r1 r2
res -> res
{- semantic comparisons -}
{-|
Compare for equality two intervals interpreted as approximations for
2 single real numbers. When equality or inequality cannot
be established, return Nothing (ie bottom).
-}
erintvEqualReals ::
(B.ERRealBase b) =>
ERInterval b ->
ERInterval b ->
Maybe Bool
erintvEqualReals (ERInterval l1 r1) (ERInterval l2 r2)
| l1 == r1 && l2 == r2 && l1 == l2 = Just True
| r1 < l2 || l1 > r2 = Just False
| otherwise = Nothing
{-|
Compare in natural order two intervals interpreted as approximations for
2 single real numbers. When equality or inequality cannot
be established, return Nothing (ie bottom).
-}
erintvCompareReals ::
(B.ERRealBase b) =>
ERInterval b ->
ERInterval b ->
Maybe Ordering
erintvCompareReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2)
| r1 < l2 = Just LT
| l1 > r2 = Just GT
| l1 == r1 && l2 == r2 && l1 == l2 = Just EQ
| otherwise = Nothing
{-|
Compare in natural order two intervals interpreted as approximations for
2 single real numbers. When relaxed equality cannot
be established nor disproved, return Nothing (ie bottom).
-}
erintvLeqReals ::
(B.ERRealBase b) =>
ERInterval b ->
ERInterval b ->
Maybe Bool
erintvLeqReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2)
| r1 <= l2 = Just True
| l1 > r2 = Just False
| otherwise = Nothing
{-|
Default splitting:
> [-Infty,+Infty] |-> [-Infty,0] [0,+Infty]
> [-Infty,x] |-> [-Infty,2*x-1] [2*x-1, x] (x <= 0)
> [-Infty,x] |-> [-Infty,0] [0, x] (x > 0)
> [x,+Infty] |-> [x,2*x+1] [2*x+1,+Infty] (x => 0)
> [x,+Infty] |-> [x,0] [0,+Infty] (x < 0)
> [x,y] |-> [x, (x+y)/2] [(x+y)/2, y]
> empty |-> empty empty
-}
erintvDefaultBisectPt ::
(B.ERRealBase b) =>
Granularity ->
(ERInterval b) ->
(ERInterval b)
erintvDefaultBisectPt gran (ERInterval l r) =
ERInterval m m
where
m =
case (B.isMinusInfinity l, B.isPlusInfinity r, B.isPlusInfinity l, B.isMinusInfinity r) of
(True, True, _, _) -> 0 -- [-oo,+oo]
(True, _,_,True) -> B.minusInfinity -- [-oo,-oo]
(_, True,True,_) -> B.plusInfinity -- [+oo,+oo]
(True, _,_,_) | r > 0 -> 0
(True, _,_,_) -> 2 * (B.setMinGranularity gran r) - 1
(_,True,_,_) | l < 0 -> 0
(_,True,_,_) -> 2 * (B.setMinGranularity gran l) + 1
(_,_,True,_) | r < 0 -> 0
(_,_,True,_) -> 2 * (B.setMinGranularity gran r) + 1
(_,_,_,True) | l > 0 -> 0
(_,_,_,True) -> 2 * (B.setMinGranularity gran l) - 1
_ -> ((B.setMinGranularity gran l) + r)/2 -- no infinities
erintvBisect ::
(B.ERRealBase b) =>
Granularity ->
(Maybe (ERInterval b)) ->
(ERInterval b) ->
(ERInterval b, ERInterval b)
erintvBisect gran maybePt i@(ERInterval l r) =
(ERInterval l mR, ERInterval mL r)
where
ERInterval mL mR = m
m =
case maybePt of
Just m -> m
Nothing -> erintvDefaultBisectPt gran i
instance (B.ERRealBase b) => Eq (ERInterval b) where
i1 == i2 =
case erintvEqualReals i1 i2 of
Nothing ->
error $
"ERInterval: Eq: comparing overlapping intervals:\n" ++
show i1 ++ "\n" ++
show i2
Just b -> b
instance (B.ERRealBase b) => Ord (ERInterval b) where
compare i1 i2 =
case erintvCompareReals i1 i2 of
Nothing ->
error $
"ERInterval: Ord: comparing overlapping intervals:\n" ++
show i1 ++ "\n" ++
show i2
Just r -> r
{- max:
(Default implementation is wrong in this case:
eg compare is not defined for overlapping intervals.)
-}
max i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) =
ERInterval (max l1 l2) (max r1 r2)
{- min: -}
min i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) =
ERInterval (min l1 l2) (min r1 r2)
instance (B.ERRealBase b) => Show (ERInterval b)
where
show = erintvShow 16 True False
erintvShow numDigits showGran showComponents interval =
showERI interval
where
showERI (ERInterval l r)
| (B.isMinusInfinity r) && (B.isPlusInfinity r) =
"[ANY]"
| l == r = "<" ++ showBase l ++ ">"
| l > r =
"[!" ++ showBase l ++ "," ++ showBase r ++ "!]"
| otherwise =
"[" ++ showBase l ++ "," ++ showBase r ++ "]"
showBase = B.showDiGrCmp numDigits showGran showComponents
instance (B.ERRealBase b, H.HTML b) => H.HTML (ERInterval b)
where
toHtml (ERInterval l r)
| l == r =
H.toHtml $ show l
| otherwise =
H.simpleTable [] [] [[H.toHtml l],[H.toHtml r]]
instance (B.ERRealBase b) => Num (ERInterval b) where
fromInteger n =
ERInterval (B.fromIntegerDown n) (B.fromIntegerUp n)
{- abs -}
abs (ERInterval l r)
| l <= 0 && r >= 0 = ERInterval 0 (max (-l) r)
| l >= 0 && r <= 0 = ERInterval (max l (-r)) 0
| r <= 0 = ERInterval (-r) (-l)
| otherwise = ERInterval l r
{- signum -}
signum i@(ERInterval l r) =
error "ER.Real.Approx.Interval: signum not implemented for ERInterval"
-- | l < 0 && r > 0 = ERInterval (-1) 1 -- need many-valuedness via sequences of intervals
-- | r < 0 = ERInterval (-1) (-1)
-- | l > 0 = ERInterval 1 1
-- | l == 0 && r == 0 = i
-- | l == 0 = ERInterval 0 1
-- | r == 0 = ERInterval (-1) 0
{- negate -}
negate (ERInterval l r) = (ERInterval (-r) (-l))
{- addition -}
i1@(ERInterval l1 r1) + i2@(ERInterval l2 r2) =
normaliseERIntervalOuter $
ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2)
{- multiplication -}
i1@(ERInterval l1 r1) * i2@(ERInterval l2 r2) =
normaliseERIntervalOuter $
intervalTimes timesDown timesUp i1 i2
instance (B.ERRealBase b) => Fractional (ERInterval b) where
fromRational rat =
(fromInteger $ numerator rat)
/ (fromInteger $ denominator rat)
{- division -}
recip i@(ERInterval l r)
| not $ RA.isConsistent i =
RA.toggleConsistency $
1 /: (RA.toggleConsistency i)
| 0 < l || r < 0 =
normaliseERIntervalOuter $
ERInterval (1 `divideDown` r) (1 `divideUp` l)
| otherwise =
RA.bottomApprox
instance (B.ERRealBase b) => RA.ERInnerOuterApprox (ERInterval b)
where
{- addition -}
i1@(ERInterval l1 r1) +: i2@(ERInterval l2 r2) =
normaliseERIntervalInner $
ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2)
{- multiplication -}
i1@(ERInterval l1 r1) *: i2@(ERInterval l2 r2) =
normaliseERIntervalInner $
intervalTimes timesUp timesDown i1 i2
{- division -}
i1@(ERInterval l1 r1) /: i2@(ERInterval l2 r2)
| not $ RA.isConsistent i2 =
(*:) i1 $
RA.toggleConsistency $
1 / (RA.toggleConsistency i2)
| 0 < l2 || r2 < 0 =
(*:) i1 $
normaliseERIntervalInner $
ERInterval (1 `divideDown` r2) (1 `divideUp` l2)
| otherwise =
RA.bottomApprox
{- setMinGranularityInner -}
setMinGranularityInner gr (ERInterval l r) =
normaliseERIntervalInner $
(ERInterval (B.setMinGranularity gr l) (negate $ B.setMinGranularity gr (-r)))
{- setGranularityInner -}
setGranularityInner gr (ERInterval l r) =
normaliseERIntervalInner $
(ERInterval (B.setGranularity gr l) (negate $ B.setGranularity gr (- r)))
intervalTimes timesL timesR i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) =
ERInterval l r
where
(l,r) =
case (compare l1 0, compare r1 0, l1 <= r1, compare l2 0, compare r2 0, l2 <= r2) of
-- i1 negative, i2 positive
(LT, LT, _, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` l2)
-- i1 negative, i2 negative
(LT, LT, _, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` l2)
-- i1 negative, i2 consistent and containing zero
(LT, LT, _, _, _, True) -> (l1 `timesL` r2, l1 `timesR` l2)
-- i1 negative, i2 inconsistent and anti-containing zero
(LT, LT, _, _, _, False) -> (r1 `timesL` r2, r1 `timesR` l2)
-- i1 positive, i2 positive
(GT, GT, _, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` r2)
-- i1 positive, i2 negative
(GT, GT, _, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` r2)
-- i1 positive, i2 consistent and containing zero
(GT, GT, _, _, _, True) -> (r1 `timesL` l2, r1 `timesR` r2)
-- i1 positive, i2 inconsistent and anti-containing zero
(GT, GT, _, _, _, False) -> (l1 `timesL` l2, l1 `timesR` r2)
-- i1 consistent and containing zero, i2 positive
(_, _, True, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` r2)
-- i1 consistent and containing zero, i2 negative
(_, _, True, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` l2)
-- i1 consistent and containing zero, i2 consistent and containing zero
(_, _, True, _, _, True) ->
(l,r)
where
l | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity
| otherwise = min l1r2 r1l2
where
l1r2 = l1 `timesL` r2
r1l2 = r1 `timesL` l2
r | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity
| otherwise = max l1l2 r1r2
where
l1l2 = l1 `timesR` l2
r1r2 = r1 `timesR` r2
-- i1 consistent and containing zero, i2 inconsistent and anti-containing zero
(_, _, True, _, _, False) -> (0, 0)
-- i1 inconsistent and anti-containing zero, i2 positive
(_, _, False, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` l2)
-- i1 inconsistent and anti-containing zero, i2 negative
(_, _, False, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` r2)
-- i1 inconsistent and anti-containing zero, i2 consistent and containing zero
(_, _, False, _, _, True) -> (0, 0)
-- i1 inconsistent and anti-containing zero, i2 the same
(_, _, False, _, _, False) ->
(l,r)
where
l | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity
| otherwise = max l1l2 r1r2
where
l1l2 = l1 `timesL` l2
r1r2 = r1 `timesL` r2
r | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity
| otherwise = min l1r2 r1l2
where
l1r2 = l1 `timesR` r2
r1l2 = r1 `timesR` l2
instance (B.ERRealBase b) => RA.ERApprox (ERInterval b) where
initialiseBaseArithmetic _ =
B.initialiseBaseArithmetic (0 :: b)
getPrecision i = erintvPrecision i
getGranularity i = erintvGranularity i
{- setMinGranularity -}
setMinGranularityOuter gr (ERInterval l r) =
normaliseERIntervalOuter $
(ERInterval (- (B.setMinGranularity gr (-l))) (B.setMinGranularity gr r))
{- setGranularity -}
setGranularityOuter gr (ERInterval l r) =
normaliseERIntervalOuter $
(ERInterval (- (B.setGranularity gr (-l))) (B.setGranularity gr r))
{- isBottom -}
isBottom (ERInterval l r) =
B.isMinusInfinity l && B.isPlusInfinity r
{- bottomApprox -}
bottomApprox =
ERInterval B.minusInfinity B.plusInfinity
{- isExact -}
isExact (ERInterval l r) = l == r
{- isConsistent -}
isConsistent (ERInterval l r) = l <= r
{- isAnticonsistent -}
isAnticonsistent (ERInterval l r) = l >= r
{- toggleConsistency -}
toggleConsistency (ERInterval l r) = (ERInterval r l)
{- isTop -}
isTop (ERInterval l r) =
B.isPlusInfinity l && B.isMinusInfinity r
{- topApprox -}
topApprox =
ERInterval B.plusInfinity B.minusInfinity
{- isBounded -}
isBounded (ERInterval l r) =
(- B.plusInfinity) < l && l < B.plusInfinity
&&
(- B.plusInfinity) < r && r < B.plusInfinity
{- plusInfinity -}
plusInfinity = ERInterval B.plusInfinity B.plusInfinity
{- refines -}
refines (ERInterval l1 r1) (ERInterval l2 r2) =
l2 <= l1 && r1 <= r2
{- maybeRefines -}
maybeRefines i1 i2 = Just $ RA.refines i1 i2
{- intersection -}
(ERInterval l1 r1) /\ (ERInterval l2 r2) =
ERInterval (max l1 l2) (min r1 r2)
{- intersectMeasureImprovement -}
intersectMeasureImprovement ix i1 i2 =
(isec, impr)
where
isec = i1 RA./\ i2
impr
| 0 `RA.refines` isecWidth && 0 `RA.refines` i1Width = 1 -- 0 -> 0 is no improvement
| otherwise = i1Width / isecWidth
i1Width = i1H - i1L
isecWidth = isecH - isecL
(isecL, isecH) = RA.bounds $ RA.setMinGranularityOuter gran isec
(i1L, i1H) = RA.bounds $ RA.setMinGranularityOuter gran i1
gran = effIx2gran ix
{- semantic comparisons -}
equalReals = erintvEqualReals
compareReals = erintvCompareReals
leqReals = erintvLeqReals
{- non-semantic comparisons -}
equalApprox = erintvEqualApprox
compareApprox = erintvCompareApprox
{- conversion from Double -}
double2ra d =
ERInterval b b
where
b = B.fromDouble d
{- formatting -}
showApprox = erintvShow
instance (B.ERRealBase b) => RA.ERIntApprox (ERInterval b)
where
doubleBounds (ERInterval l r) =
(negate $ B.toDouble (-l), B.toDouble r)
floatBounds (ERInterval l r) =
(negate $ B.toFloat (-l), B.toFloat r)
integerBounds (ERInterval l r) =
(negate $ mkEI (- l), mkEI r)
where
mkEI f
| B.isPlusInfinity f = EI.PlusInfinity
| B.isMinusInfinity f = EI.MinusInfinity
| otherwise = ceiling f
defaultBisectPt dom =
erintvDefaultBisectPt (RA.getGranularity dom + 1) dom
bisectDomain maybePt dom =
erintvBisect (RA.getGranularity dom + 1) maybePt dom
{- \/ -}
(ERInterval l1 r1) \/ (ERInterval l2 r2) =
ERInterval (min l1 l2) (max r1 r2)
{- RA.bounds -}
bounds (ERInterval l r) =
(ERInterval l l, ERInterval r r)
{- RA.fromBounds -}
fromBounds (ERInterval l1 r1, ERInterval l2 r2)
| l1 == r1 && l2 == r2 = ERInterval l1 l2
fromBounds i1i2 =
error $
"ER.Real.Approx.Interval: fromBounds: bounds not exact: "
++ show i1i2
instance (B.ERRealBase b) => RAEL.ERApproxElementary (ERInterval b)
instance (B.ERRealBase b) => RAEL.ERInnerOuterApproxElementary (ERInterval b)
-- all operations here have appropriate default implementations