--------------------------------------------------------------------------------
-- |
-- Module      :  HGeometry.Number.Real.Interval
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- Interval Arithmetic. We use doubles as initial approximation of the interval.
--------------------------------------------------------------------------------

module HGeometry.Number.Real.Interval
  ( IntervalReal
  , exactValue
  , fromExact
  ) where

import Numeric.Rounded.Hardware
import Numeric.Rounded.Hardware.Interval.Class
import Numeric.Rounded.Hardware.Interval.NonEmpty
import GHC.Generics (Generic)
import Control.DeepSeq
import Data.Bifunctor
import Control.Monad.Random
import Text.Read
--------------------------------------------------------------------------------


data IntervalReal r = IR {-# UNPACK#-}!(Interval Double) r
  deriving ((forall x. IntervalReal r -> Rep (IntervalReal r) x)
-> (forall x. Rep (IntervalReal r) x -> IntervalReal r)
-> Generic (IntervalReal r)
forall x. Rep (IntervalReal r) x -> IntervalReal r
forall x. IntervalReal r -> Rep (IntervalReal r) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall r x. Rep (IntervalReal r) x -> IntervalReal r
forall r x. IntervalReal r -> Rep (IntervalReal r) x
$cfrom :: forall r x. IntervalReal r -> Rep (IntervalReal r) x
from :: forall x. IntervalReal r -> Rep (IntervalReal r) x
$cto :: forall r x. Rep (IntervalReal r) x -> IntervalReal r
to :: forall x. Rep (IntervalReal r) x -> IntervalReal r
Generic)

-- | Get the exact vavlue of the interval-real.
exactValue          :: IntervalReal r -> r
exactValue :: forall r. IntervalReal r -> r
exactValue (IR Interval Double
_ r
x) = r
x

-- | Construct an intervalReal from a given r value.
fromExact   :: Real r => r -> IntervalReal r
fromExact :: forall r. Real r => r -> IntervalReal r
fromExact r
x = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (r -> Interval Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac r
x) r
x


instance NFData r => NFData (IntervalReal r)

instance Show r => Show (IntervalReal r) where
  showsPrec :: Int -> IntervalReal r -> ShowS
showsPrec Int
d (IR Interval Double
_ r
r) = Int -> r -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
d r
r

instance (Read r, Real r) => Read (IntervalReal r) where
  readPrec :: ReadPrec (IntervalReal r)
readPrec = r -> IntervalReal r
forall r. Real r => r -> IntervalReal r
fromExact (r -> IntervalReal r) -> ReadPrec r -> ReadPrec (IntervalReal r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ReadPrec r
forall a. Read a => ReadPrec a
readPrec

instance Eq r => Eq (IntervalReal r)  where
  (IR Interval Double
ix r
x) == :: IntervalReal r -> IntervalReal r -> Bool
== (IR Interval Double
iy r
y) = Bool -> Bool
not (Interval Double -> Interval Double -> Bool
forall i. IsInterval i => i -> i -> Bool
disjoint Interval Double
ix Interval Double
iy) Bool -> Bool -> Bool
&& r
x r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
y
         -- the disjoint test is important; only if the intervals are not disjoint check
         -- for actual equality

instance Ord r => Ord (IntervalReal r)  where
  (IR Interval Double
ix r
x) compare :: IntervalReal r -> IntervalReal r -> Ordering
`compare` (IR Interval Double
iy r
y)
    | Interval Double
ix Interval Double -> Interval Double -> Bool
forall i. IsInterval i => i -> i -> Bool
`strictPrecedes` Interval Double
iy = Ordering
LT
    | Interval Double
iy Interval Double -> Interval Double -> Bool
forall i. IsInterval i => i -> i -> Bool
`strictPrecedes` Interval Double
ix = Ordering
GT
    | Bool
otherwise              = r
x r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
y
    -- we could maybe do something more interesting here

instance (Num r, Ord r) => Num (IntervalReal r) where
  (IR Interval Double
ix r
x) + :: IntervalReal r -> IntervalReal r -> IntervalReal r
+ (IR Interval Double
iy r
y) = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double
ix Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
+ Interval Double
iy) (r
x r -> r -> r
forall a. Num a => a -> a -> a
+ r
y)
  (IR Interval Double
ix r
x) - :: IntervalReal r -> IntervalReal r -> IntervalReal r
- (IR Interval Double
iy r
y) = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double
ix Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
- Interval Double
iy) (r
x r -> r -> r
forall a. Num a => a -> a -> a
- r
y)
  (IR Interval Double
ix r
x) * :: IntervalReal r -> IntervalReal r -> IntervalReal r
* (IR Interval Double
iy r
y) = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double
ix Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
* Interval Double
iy) (r
x r -> r -> r
forall a. Num a => a -> a -> a
* r
y)
  abs :: IntervalReal r -> IntervalReal r
abs (IR Interval Double
ix r
x)    = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double -> Interval Double
forall a. Num a => a -> a
abs Interval Double
ix) (r -> r
forall a. Num a => a -> a
abs r
x)
  signum :: IntervalReal r -> IntervalReal r
signum (IR (I Rounded 'TowardNegInf Double
l Rounded 'TowardInf Double
u) r
x)
    | Double
0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Rounded 'TowardNegInf Double -> Double
forall (r :: RoundingMode) a. Rounded r a -> a
getRounded Rounded 'TowardNegInf Double
l = IntervalReal r
1
    | Double
0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Rounded 'TowardInf Double -> Double
forall (r :: RoundingMode) a. Rounded r a -> a
getRounded Rounded 'TowardInf Double
u = -IntervalReal r
1
    | Bool
otherwise        = case r
0 r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
x of
                           Ordering
LT -> IntervalReal r
1
                           Ordering
EQ -> IntervalReal r
0
                           Ordering
GT -> (-IntervalReal r
1)
  fromInteger :: Integer -> IntervalReal r
fromInteger Integer
i    = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Integer -> Interval Double
forall a. Num a => Integer -> a
fromInteger Integer
i) (Integer -> r
forall a. Num a => Integer -> a
fromInteger Integer
i)
  negate :: IntervalReal r -> IntervalReal r
negate (IR Interval Double
ix r
x) = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double -> Interval Double
forall a. Num a => a -> a
negate Interval Double
ix) (r -> r
forall a. Num a => a -> a
negate r
x)

instance (Fractional r, Ord r) => Fractional (IntervalReal r) where
  fromRational :: Rational -> IntervalReal r
fromRational Rational
x = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Rational -> Interval Double
forall a. Fractional a => Rational -> a
fromRational Rational
x) (Rational -> r
forall a. Fractional a => Rational -> a
fromRational Rational
x)
  (IR Interval Double
ix r
x) / :: IntervalReal r -> IntervalReal r -> IntervalReal r
/ (IR Interval Double
iy r
y) = Interval Double -> r -> IntervalReal r
forall r. Interval Double -> r -> IntervalReal r
IR (Interval Double
ix Interval Double -> Interval Double -> Interval Double
forall a. Fractional a => a -> a -> a
/ Interval Double
iy) (r
x r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
y)

-- FIXME: we are potentially dividing by 0 apparently

instance Real r => Real (IntervalReal r) where
  toRational :: IntervalReal r -> Rational
toRational (IR Interval Double
_ r
x) = r -> Rational
forall a. Real a => a -> Rational
toRational r
x


-- instance Uniform r => Uniform (IntervalReal r) where

instance (Random r, Real r) => Random (IntervalReal r) where
  randomR :: forall g.
RandomGen g =>
(IntervalReal r, IntervalReal r) -> g -> (IntervalReal r, g)
randomR (IntervalReal r
a,IntervalReal r
b) = (r -> IntervalReal r) -> (r, g) -> (IntervalReal r, g)
forall a b c. (a -> b) -> (a, c) -> (b, c)
forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first r -> IntervalReal r
forall r. Real r => r -> IntervalReal r
fromExact ((r, g) -> (IntervalReal r, g))
-> (g -> (r, g)) -> g -> (IntervalReal r, g)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (r, r) -> g -> (r, g)
forall g. RandomGen g => (r, r) -> g -> (r, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (IntervalReal r -> r
forall r. IntervalReal r -> r
exactValue IntervalReal r
a, IntervalReal r -> r
forall r. IntervalReal r -> r
exactValue IntervalReal r
b)
  random :: forall g. RandomGen g => g -> (IntervalReal r, g)
random = (r -> IntervalReal r) -> (r, g) -> (IntervalReal r, g)
forall a b c. (a -> b) -> (a, c) -> (b, c)
forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first r -> IntervalReal r
forall r. Real r => r -> IntervalReal r
fromExact ((r, g) -> (IntervalReal r, g))
-> (g -> (r, g)) -> g -> (IntervalReal r, g)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. g -> (r, g)
forall g. RandomGen g => g -> (r, g)
forall a g. (Random a, RandomGen g) => g -> (a, g)
random


  --   runRand $ do
  --   v <- getRandom
  --   pure $ (b-a)*abs v + a
  -- -- Generate a random number between -1 and +1 with 'maxBound::Int' discrete increments.
  -- random = runRand $ do
  --   v <- getRandom
  --   let fromInt :: Int -> Integer; fromInt = fromIntegral
  --   pure $ RealNumber $ fromInt v % fromInt maxBound