{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE ScopedTypeVariables, DeriveDataTypeable #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Interval
-- Copyright   :  (c) Masahiro Sakai 2011-2013
-- License     :  BSD-style
-- 
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable (ScopedTypeVariables, DeriveDataTypeable)
--
-- Interval datatype and interval arithmetic.
--
-- Unlike the intervals package (<http://hackage.haskell.org/package/intervals>),
-- this module provides both open and closed intervals and is intended to be used
-- with 'Rational'.
--
-- For the purpose of abstract interpretation, it might be convenient to use
-- 'Lattice' instance. See also lattices package
-- (<http://hackage.haskell.org/package/lattices>).
-- 
-----------------------------------------------------------------------------
module Data.Interval
  (
  -- * Interval type
    Interval
  , Extended (..)
  , EndPoint

  -- * Construction
  , interval
  , (<=..<=)
  , (<..<=)
  , (<=..<)
  , (<..<)
  , whole
  , empty
  , singleton

  -- * Query
  , null
  , member
  , notMember
  , isSubsetOf
  , isProperSubsetOf
  , lowerBound
  , upperBound
  , lowerBound'
  , upperBound'
  , width

  -- * Universal comparison operators
  , (<!), (<=!), (==!), (>=!), (>!), (/=!)

  -- * Existential comparison operators
  , (<?), (<=?), (==?), (>=?), (>?), (/=?)

  -- * Existential comparison operators that produce witnesses (experimental)
  , (<??), (<=??), (==??), (>=??), (>??), (/=??)

  -- * Combine
  , intersection
  , intersections
  , hull
  , hulls

  -- * Operations
  , pickup
  , simplestRationalWithin
  ) where

import Algebra.Lattice
import Control.DeepSeq
import Control.Exception (assert)
import Control.Monad hiding (join)
import Data.Data
import Data.ExtendedReal
import Data.Hashable
import Data.List hiding (null)
import Data.Maybe
import Data.Monoid
import Data.Ratio
import Prelude hiding (null)

-- | The intervals (/i.e./ connected and convex subsets) over real numbers __R__.
data Interval r = Interval !(EndPoint r, Bool) !(EndPoint r, Bool)
  deriving (Eq, Typeable)  

-- | Lower endpoint (/i.e./ greatest lower bound)  of the interval.
--
-- * 'lowerBound' of the empty interval is 'PosInf'.
--
-- * 'lowerBound' of a left unbounded interval is 'NegInf'.
--
-- * 'lowerBound' of an interval may or may not be a member of the interval.
lowerBound :: Interval r -> EndPoint r
lowerBound (Interval (lb,_) _) = lb

-- | Upper endpoint (/i.e./ least upper bound) of the interval.
--
-- * 'upperBound' of the empty interval is 'NegInf'.
--
-- * 'upperBound' of a right unbounded interval is 'PosInf'.
-- 
-- * 'upperBound' of an interval may or may not be a member of the interval.
upperBound :: Interval r -> EndPoint r
upperBound (Interval _ (ub,_)) = ub

-- | 'lowerBound' of the interval and whether it is included in the interval.
-- The result is convenient to use as an argument for 'interval'.
lowerBound' :: Interval r -> (EndPoint r, Bool)
lowerBound' (Interval lb _) = lb

-- | 'upperBound' of the interval and whether it is included in the interval.
-- The result is convenient to use as an argument for 'interval'.
upperBound' :: Interval r -> (EndPoint r, Bool)
upperBound' (Interval _ ub) = ub

instance NFData r => NFData (Interval r) where
  rnf (Interval lb ub) = rnf lb `seq` rnf ub

instance Hashable r => Hashable (Interval r) where
  hashWithSalt s (Interval lb ub) = s `hashWithSalt` lb `hashWithSalt` ub

instance (Num r, Ord r) => JoinSemiLattice (Interval r) where
  join = hull

instance (Num r, Ord r) => MeetSemiLattice (Interval r) where
  meet = intersection

instance (Num r, Ord r) => Lattice (Interval r)

instance (Num r, Ord r) => BoundedJoinSemiLattice (Interval r) where
  bottom = empty

instance (Num r, Ord r) => BoundedMeetSemiLattice (Interval r) where
  top = whole

instance (Num r, Ord r) => BoundedLattice (Interval r)

instance (Num r, Ord r, Show r) => Show (Interval r) where
  showsPrec _ x | null x = showString "empty"
  showsPrec p x = showParen (p > appPrec) $
    showString "interval " .
    showsPrec appPrec1 (lowerBound' x) .
    showChar ' ' . 
    showsPrec appPrec1 (upperBound' x)

instance (Num r, Ord r, Read r) => Read (Interval r) where
  readsPrec p r =
    (readParen (p > appPrec) $ \s0 -> do
      ("interval",s1) <- lex s0
      (lb,s2) <- readsPrec (appPrec+1) s1
      (ub,s3) <- readsPrec (appPrec+1) s2
      return (interval lb ub, s3)) r
    ++
    (do ("empty", s) <- lex r
        return (empty, s))

-- This instance preserves data abstraction at the cost of inefficiency.
-- We omit reflection services for the sake of data abstraction.

instance (Num r, Ord r, Data r) => Data (Interval r) where
  gfoldl k z x   = z interval `k` lowerBound' x `k` upperBound' x
  toConstr _     = error "toConstr"
  gunfold _ _    = error "gunfold"
  dataTypeOf _   = mkNoRepType "Data.Interval.Interval"
  dataCast1 f    = gcast1 f

-- | smart constructor for 'Interval'
interval
  :: (Ord r, Num r)
  => (EndPoint r, Bool) -- ^ lower bound and whether it is included 
  -> (EndPoint r, Bool) -- ^ upper bound and whether it is included
  -> Interval r
interval lb@(x1,in1) ub@(x2,in2) =
  case x1 `compare` x2 of
    GT -> empty --  empty interval
    LT -> Interval (normalize lb) (normalize ub)
    EQ -> if in1 && in2 && isFinite x1 then Interval lb ub else empty
  where
    normalize x@(Finite _, _) = x
    normalize (x, _) = (x, False)

-- | closed interval [@l@,@u@]
(<=..<=)
  :: (Ord r, Num r)
  => EndPoint r -- ^ lower bound @l@
  -> EndPoint r -- ^ upper bound @u@
  -> Interval r
(<=..<=) lb ub = interval (lb, True) (ub, True)

-- | left-open right-closed interval (@l@,@u@]
(<..<=)
  :: (Ord r, Num r)
  => EndPoint r -- ^ lower bound @l@
  -> EndPoint r -- ^ upper bound @u@
  -> Interval r
(<..<=) lb ub = interval (lb, False) (ub, True)

-- | left-closed right-open interval [@l@, @u@)
(<=..<)
  :: (Ord r, Num r)
  => EndPoint r -- ^ lower bound @l@
  -> EndPoint r -- ^ upper bound @u@
  -> Interval r
(<=..<) lb ub = interval (lb, True) (ub, False)

-- | open interval (@l@, @u@)
(<..<)
  :: (Ord r, Num r)
  => EndPoint r -- ^ lower bound @l@
  -> EndPoint r -- ^ upper bound @u@
  -> Interval r
(<..<) lb ub = interval (lb, False) (ub, False)

-- | whole real number line (-∞, ∞)
whole :: (Num r, Ord r) => Interval r
whole = Interval (NegInf, False) (PosInf, False)

-- | empty (contradicting) interval
empty :: Num r => Interval r
empty = Interval (PosInf, False) (NegInf, False)

-- | singleton set \[x,x\]
singleton :: (Num r, Ord r) => r -> Interval r
singleton x = interval (Finite x, True) (Finite x, True)

-- | intersection of two intervals
intersection :: forall r. (Ord r, Num r) => Interval r -> Interval r -> Interval r
intersection (Interval l1 u1) (Interval l2 u2) = interval (maxLB l1 l2) (minUB u1 u2)
  where
    maxLB :: (EndPoint r, Bool) -> (EndPoint r, Bool) -> (EndPoint r, Bool)
    maxLB (x1,in1) (x2,in2) =
      ( max x1 x2
      , case x1 `compare` x2 of
          EQ -> in1 && in2
          LT -> in2
          GT -> in1
      )
    minUB :: (EndPoint r, Bool) -> (EndPoint r, Bool) -> (EndPoint r, Bool)
    minUB (x1,in1) (x2,in2) =
      ( min x1 x2
      , case x1 `compare` x2 of
          EQ -> in1 && in2
          LT -> in1
          GT -> in2
      )

-- | intersection of a list of intervals.
--
-- Since 0.6.0
intersections :: (Ord r, Num r) => [Interval r] -> Interval r
intersections xs = foldl' intersection whole xs

-- | convex hull of two intervals
hull :: forall r. (Ord r, Num r) => Interval r -> Interval r -> Interval r
hull x1 x2
  | null x1 = x2
  | null x2 = x1
hull (Interval l1 u1) (Interval l2 u2) = interval (minLB l1 l2) (maxUB u1 u2)
  where
    maxUB :: (EndPoint r, Bool) -> (EndPoint r, Bool) -> (EndPoint r, Bool)
    maxUB (x1,in1) (x2,in2) =
      ( max x1 x2
      , case x1 `compare` x2 of
          EQ -> in1 || in2
          LT -> in2
          GT -> in1
      )
    minLB :: (EndPoint r, Bool) -> (EndPoint r, Bool) -> (EndPoint r, Bool)
    minLB (x1,in1) (x2,in2) =
      ( min x1 x2
      , case x1 `compare` x2 of
          EQ -> in1 || in2
          LT -> in1
          GT -> in2
      )

-- | convex hull of a list of intervals.
--
-- Since 0.6.0
hulls :: (Ord r, Num r) => [Interval r] -> Interval r
hulls xs = foldl' hull empty xs

-- | Is the interval empty?
null :: Ord r => Interval r -> Bool
null (Interval (x1,in1) (x2,in2)) = 
  case x1 `compare` x2 of
    EQ -> assert (in1 && in2) False
    LT -> False
    GT -> True

isSingleton :: Ord r => Interval r -> Bool
isSingleton (Interval (Finite l, True) (Finite u, True)) = l==u
isSingleton _ = False

-- | Is the element in the interval?
member :: Ord r => r -> Interval r -> Bool
member x (Interval (x1,in1) (x2,in2)) = condLB && condUB
  where
    condLB = if in1 then x1 <= Finite x else x1 < Finite x
    condUB = if in2 then Finite x <= x2 else Finite x < x2

-- | Is the element not in the interval?
notMember :: Ord r => r -> Interval r -> Bool
notMember a i = not $ member a i

-- | Is this a subset?
-- @(i1 \``isSubsetOf`\` i2)@ tells whether @i1@ is a subset of @i2@.
isSubsetOf :: Ord r => Interval r -> Interval r -> Bool
isSubsetOf (Interval lb1 ub1) (Interval lb2 ub2) = testLB lb1 lb2 && testUB ub1 ub2
  where
    testLB (x1,in1) (x2,in2) =
      case x1 `compare` x2 of
        GT -> True
        LT -> False
        EQ -> not in1 || in2 -- in1 => in2
    testUB (x1,in1) (x2,in2) =
      case x1 `compare` x2 of
        LT -> True
        GT -> False
        EQ -> not in1 || in2 -- in1 => in2

-- | Is this a proper subset? (/i.e./ a subset but not equal).
isProperSubsetOf :: Ord r => Interval r -> Interval r -> Bool
isProperSubsetOf i1 i2 = i1 /= i2 && i1 `isSubsetOf` i2

-- | Width of a interval. Width of an unbounded interval is @undefined@.
width :: (Num r, Ord r) => Interval r -> r
width x | null x = 0
width (Interval (Finite l, _) (Finite u, _)) = u - l
width _ = error "Data.Interval.width: unbounded interval"

-- | pick up an element from the interval if the interval is not empty.
pickup :: (Real r, Fractional r) => Interval r -> Maybe r
pickup (Interval (NegInf,_) (PosInf,_))   = Just 0
pickup (Interval (Finite x1, in1) (PosInf,_)) = Just $ if in1 then x1 else x1+1
pickup (Interval (NegInf,_) (Finite x2, in2)) = Just $ if in2 then x2 else x2-1
pickup (Interval (Finite x1, in1) (Finite x2, in2)) =
  case x1 `compare` x2 of
    GT -> Nothing
    LT -> Just $ (x1+x2) / 2
    EQ -> if in1 && in2 then Just x1 else Nothing
pickup _ = Nothing

-- | 'simplestRationalWithin' returns the simplest rational number within the interval.
-- 
-- A rational number @y@ is said to be /simpler/ than another @y'@ if
--
-- * @'abs' ('numerator' y) <= 'abs' ('numerator' y')@, and
--
-- * @'denominator' y <= 'denominator' y'@.
-- 
-- (see also 'approxRational')
--
-- Since 0.4.0
simplestRationalWithin :: RealFrac r => Interval r -> Maybe Rational
simplestRationalWithin i | null i = Nothing
simplestRationalWithin i 
  | 0 <! i    = Just $ go i 
  | i <! 0    = Just $ - go (- i)
  | otherwise = assert (0 `member` i) $ Just $ 0
  where
    go i
      | fromInteger lb_floor       `member'` i = fromInteger lb_floor
      | fromInteger (lb_floor + 1) `member'` i = fromInteger (lb_floor + 1)
      | otherwise = fromInteger lb_floor + recip (go (recip (i - singleton (fromInteger lb_floor))))
      where
        Finite lb = lowerBound i
        lb_floor  = floor lb

    member' :: (Real r, Fractional r) => Rational -> Interval r -> Bool
    member' x (Interval (x1,in1) (x2,in2)) = condLB && condUB
      where
        x' = fromRational x
        condLB = if in1 then x1 <= Finite x' else x1 < Finite x'
        condUB = if in2 then Finite x' <= x2 else Finite x' < x2

-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@?
(<!) :: Real r => Interval r -> Interval r -> Bool
a <! b =
  case ub_a `compare` lb_b of
    LT -> True
    GT -> False
    EQ ->
      case ub_a of
        NegInf   -> True -- a is empty, so it holds vacuously
        PosInf   -> True -- b is empty, so it holds vacuously
        Finite _ -> not (in1 && in2)
  where
    (ub_a, in1) = upperBound' a
    (lb_b, in2) = lowerBound' b

-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@?
(<=!) :: Real r => Interval r -> Interval r -> Bool
a <=! b = upperBound a <= lowerBound b

-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@?
(==!) :: Real r => Interval r -> Interval r -> Bool
a ==! b = a <=! b && a >=! b

-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@?
--
-- Since 1.0.1
(/=!) :: Real r => Interval r -> Interval r -> Bool
a /=! b = null $ a `intersection` b

-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@?
(>=!) :: Real r => Interval r -> Interval r -> Bool
(>=!) = flip (<=!)

-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@?
(>!) :: Real r => Interval r -> Interval r -> Bool
(>!) = flip (<!)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
(<?) :: Real r => Interval r -> Interval r -> Bool
a <? b = lb_a < ub_b
  where
    lb_a = lowerBound a
    ub_b = upperBound b

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?
--
-- Since 1.0.0
(<??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
a <?? b = do
  guard $ lowerBound a < upperBound b
  let c = intersection a b
  case pickup c of
    Nothing -> do
      x <- pickup a
      y <- pickup b
      return (x,y)
    Just z -> do
      let x:y:_ = take 2 $
                    maybeToList (pickup (intersection a (NegInf <..< Finite z))) ++
                    [z] ++
                    maybeToList (pickup (intersection b (Finite z <..< PosInf)))
      return (x,y)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
(<=?) :: Real r => Interval r -> Interval r -> Bool
a <=? b　=
  case lb_a `compare` ub_b of
    LT -> True
    GT -> False
    EQ -> 
      case lb_a of
        NegInf -> False -- b is empty
        PosInf -> False -- a is empty
        Finite _ -> in1 && in2
  where
    (lb_a, in1) = lowerBound' a
    (ub_b, in2) = upperBound' b

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?
-- 
-- Since 1.0.0
(<=??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
a <=?? b　= do
  case pickup (intersection a b) of
    Just x -> return (x,x)
    Nothing -> do
      guard $ upperBound a <= lowerBound b
      x <- pickup a
      y <- pickup b
      return (x,y)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
-- 
-- Since 1.0.0
(==?) :: Real r => Interval r -> Interval r -> Bool
a ==? b = not $ null $ intersection a b

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?
--
-- Since 1.0.0
(==??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
a ==?? b = do
  x <- pickup (intersection a b)
  return (x,x)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
--
-- Since 1.0.1
(/=?) :: Real r => Interval r -> Interval r -> Bool
a /=? b = not (null a) && not (null b) && not (a == b && isSingleton a)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?
--
-- Since 1.0.1
(/=??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
a /=?? b = do
  guard $ not $ null a
  guard $ not $ null b
  guard $ not $ a == b && isSingleton a
  if not (isSingleton b)
    then f a b
    else liftM (\(y,x) -> (x,y)) $ f b a
  where
    f a b = do
      x <- pickup a
      y <- msum [pickup (b `intersection` c) | c <- [NegInf <..< Finite x, Finite x <..< PosInf]]
      return (x,y)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
(>=?) :: Real r => Interval r -> Interval r -> Bool
(>=?) = flip (<=?)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
(>?) :: Real r => Interval r -> Interval r -> Bool
(>?) = flip (<?)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?
--
-- Since 1.0.0
(>=??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
(>=??) = flip (<=??)

-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?
--
-- Since 1.0.0
(>??) :: (Real r, Fractional r) => Interval r -> Interval r -> Maybe (r,r)
(>??) = flip (<??)

appPrec, appPrec1 :: Int
appPrec = 10
appPrec1 = appPrec + 1

scaleInterval :: (Num r, Ord r) => r -> Interval r -> Interval r
scaleInterval _ x | null x = empty
scaleInterval c (Interval lb ub) =
  case compare c 0 of
    EQ -> singleton 0
    LT -> interval (scaleInf' c ub) (scaleInf' c lb)
    GT -> interval (scaleInf' c lb) (scaleInf' c ub)

instance (Num r, Ord r) => Num (Interval r) where
  a + b | null a || null b = empty
  Interval lb1 ub1 + Interval lb2 ub2 = interval (f lb1 lb2) (g ub1 ub2)
    where
      f (Finite x1, in1) (Finite x2, in2) = (Finite (x1+x2), in1 && in2)
      f (NegInf,_) _ = (NegInf, False)
      f _ (NegInf,_) = (NegInf, False)
      f _ _ = error "Interval.(+) should not happen"

      g (Finite x1, in1) (Finite x2, in2) = (Finite (x1+x2), in1 && in2)
      g (PosInf,_) _ = (PosInf, False)
      g _ (PosInf,_) = (PosInf, False)
      g _ _ = error "Interval.(+) should not happen"

  negate a = scaleInterval (-1) a

  fromInteger i = singleton (fromInteger i)

  abs x = ((x `intersection` nonneg) `hull` (negate x `intersection` nonneg))
    where
      nonneg = 0 <=..< PosInf

  signum x = zero `hull` pos `hull` neg
    where
      zero = if member 0 x then singleton 0 else empty
      pos = if null $ (0 <..< PosInf) `intersection` x
            then empty
            else singleton 1
      neg = if null $ (NegInf <..< 0) `intersection` x
            then empty
            else singleton (-1)

  a * b | null a || null b = empty
  Interval lb1 ub1 * Interval lb2 ub2 = interval lb3 ub3
    where
      xs = [ mulInf' x1 x2 | x1 <- [lb1, ub1], x2 <- [lb2, ub2] ]
      ub3 = maximumBy cmpUB xs
      lb3 = minimumBy cmpLB xs

instance forall r. (Real r, Fractional r) => Fractional (Interval r) where
  fromRational r = singleton (fromRational r)
  recip a | null a = empty
  recip i | 0 `member` i = whole -- should be error?
  recip (Interval lb ub) = interval lb3 ub3
    where
      ub3 = maximumBy cmpUB xs
      lb3 = minimumBy cmpLB xs
      xs = [recipLB lb, recipUB ub]

cmpUB, cmpLB :: Ord r => (EndPoint r, Bool) -> (EndPoint r, Bool) -> Ordering
cmpUB (x1,in1) (x2,in2) = compare x1 x2 `mappend` compare in1 in2
cmpLB (x1,in1) (x2,in2) = compare x1 x2 `mappend` flip compare in1 in2

-- | Endpoints of intervals
type EndPoint r = Extended r

scaleInf' :: (Num r, Ord r) => r -> (EndPoint r, Bool) -> (EndPoint r, Bool)
scaleInf' a (x1, in1) = (scaleEndPoint a x1, in1)

scaleEndPoint :: (Num r, Ord r) => r -> EndPoint r -> EndPoint r
scaleEndPoint a inf =
  case a `compare` 0 of
    EQ -> 0
    GT ->
      case inf of
        NegInf   -> NegInf
        Finite b -> Finite (a*b)
        PosInf   -> PosInf
    LT ->
      case inf of
        NegInf   -> PosInf
        Finite b -> Finite (a*b)
        PosInf   -> NegInf

mulInf' :: (Num r, Ord r) => (EndPoint r, Bool) -> (EndPoint r, Bool) -> (EndPoint r, Bool)
mulInf' (0, True) _ = (0, True)
mulInf' _ (0, True) = (0, True)
mulInf' (x1,in1) (x2,in2) = (x1*x2, in1 && in2)

recipLB :: (Fractional r, Ord r) => (EndPoint r, Bool) -> (EndPoint r, Bool)
recipLB (0, _) = (PosInf, False)
recipLB (x1, in1) = (recip x1, in1)

recipUB :: (Fractional r, Ord r) => (EndPoint r, Bool) -> (EndPoint r, Bool)
recipUB (0, _) = (NegInf, False)
recipUB (x1, in1) = (recip x1, in1)