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