{-# LANGUAGE CPP #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE ExtendedDefaultRules #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -Wall #-}
#if ( __GLASGOW_HASKELL__ < 820 )
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
{-# OPTIONS_GHC -fno-warn-unrecognised-pragmas #-}
#endif

-- | representation of a continuous range of a type
module NumHask.Range
  ( Range(..)
  , pattern Range
  , gridSensible
 ) where

import NumHask.Prelude hiding (singleton)
import NumHask.Space

import Data.Functor.Apply (Apply(..))
import Data.Semigroup.Foldable (Foldable1(..))
import Data.Semigroup.Traversable (Traversable1(..))
import Data.Functor.Classes
import Data.Distributive
import Test.QuickCheck.Arbitrary (Arbitrary(..))
import qualified Text.Show as Show

-- $setup
-- >>> :set -XNoImplicitPrelude

-- | A continuous range over type a
--
-- >>> let a = Range (-1) 1
-- >>> a
-- Range -1 1
-- >>> fmap (+1) (Range 1 2)
-- Range 2 3
-- >>> one :: Range Double
-- Range -0.5 0.5
-- >>> zero :: Range Double
-- Range Infinity -Infinity

-- | as a Field instance
--
-- >>> Range 0 1 + zero
-- Range 0.0 1.0
-- >>> Range 0 1 + Range 2 3
-- Range 0.0 3.0
-- >>> Range 1 1 - one
-- Range 0.5 1.0
-- >>> Range 0 1 * one
-- Range 0.0 1.0
-- >>> Range 0 1 / one
-- Range 0.0 1.0
-- >>> singleton 2.3 :: Range Double
-- Range 2.3 2.3
-- >>> abs (Range 1 0)
-- Range 0.0 1.0
-- >>> sign (Range 1 0) == negate one
-- True
--
-- Idempotent
--
-- >>> Range 0 2 + Range 0 2
-- Range 0.0 2.0
--
-- as a space instance
--
-- >>> project (Range 0 1) (Range 1 4) 0.5
-- 2.5
-- >>> grid OuterPos (Range 0 10) 5
-- [0.0,2.0,4.0,6.0,8.0,10.0]
-- >>> gridSpace (Range 0 1) 4
-- [Range 0.0 0.25,Range 0.25 0.5,Range 0.5 0.75,Range 0.75 1.0]
-- >>> gridSensible OuterPos (Range (-12) 23) 6
-- [-10.0,-5.0,0.0,5.0,10.0,15.0,20.0]
newtype Range a = Range' (a,a)
  deriving (Eq, Generic)

-- | A tuple is the preferred concrete implementation of a Range, due to many libraries having substantial optimizations for tuples already (eg 'Vector').  'Pattern Synonyms' allow us to recover a constructor without the need for tuple syntax.
pattern Range :: a -> a -> Range a
pattern Range a b = Range' (a, b)
{-# COMPLETE Range#-}

instance (Show a) => Show (Range a) where
    show (Range a b) = "Range " <> show a <> " " <> show b

instance Eq1 Range where
    liftEq f (Range a b) (Range c d) = f a c && f b d

instance Show1 Range where
    liftShowsPrec sp _ d (Range' (a,b)) = showsBinaryWith sp sp "Range" d a b

instance Functor Range where
    fmap f (Range a b) = Range (f a) (f b)

instance Apply Range where
  Range fa fb <.> Range a b = Range (fa a) (fb b)

instance Applicative Range where
    pure a = Range a a
    (Range fa fb) <*> Range a b = Range (fa a) (fb b)

instance Monad Range where
  Range a b >>= f = Range a' b' where
    Range a' _ = f a
    Range _ b' = f b

instance Foldable Range where
    foldMap f (Range a b) = f a `mappend` f b

instance Foldable1 Range

instance Traversable Range where
    traverse f (Range a b) = Range <$> f a <*> f b

instance Traversable1 Range where
    traverse1 f (Range a b) = Range <$> f a Data.Functor.Apply.<.> f b

instance Distributive Range where
  collect f x = Range (getL . f <$> x) (getR . f <$> x)
    where getL (Range l _) = l
          getR (Range _ r) = r

instance Representable Range where
  type Rep Range = Bool
  tabulate f = Range (f False) (f True)
  index (Range l _) False = l
  index (Range _ r) True = r

instance (Arbitrary a) => Arbitrary (Range a) where
    arbitrary = do
        a <- arbitrary
        b <- arbitrary
        pure (Range a b)

instance NFData a => NFData (Range a) where
  rnf (Range a b) = rnf a `seq` rnf b

two :: (MultiplicativeUnital a, Additive a) => a
two = one + one

half :: (Field a) => a
half = one / two

-- | convex hull union
instance (FromInteger a, Ord a, BoundedField a) => AdditiveMagma (Range a) where
    plus (Range l0 u0) (Range l1 u1) = Range (min l0 l1) (max u0 u1) 

instance (FromInteger a, Ord a, BoundedField a) => AdditiveUnital (Range a) where
    zero = Range infinity neginfinity

instance (FromInteger a, Ord a, BoundedField a) => AdditiveAssociative (Range a)

instance (FromInteger a, Ord a, BoundedField a) => AdditiveInvertible (Range a) where
    negate (Range l u) = Range u l

instance (FromInteger a, Ord a, BoundedField a) => AdditiveCommutative (Range a)
instance (FromInteger a, Ord a, BoundedField a) => Additive (Range a)
instance (FromInteger a, Ord a, BoundedField a) => AdditiveGroup (Range a)

-- | times may well be some sort of affine projection lurking under the hood
instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeMagma (Range a) where
    times a b = Range (m - r/two) (m + r/two)
        where
          m = mid a + mid b
          r = width a * width b

-- | The unital object (Range -0.5 0.5) satisfies:
--
-- > width one = one
-- > mid zero = zero
instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeUnital (Range a) where
    one = Range (negate half) half

instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeAssociative (Range a)

instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeInvertible (Range a) where
    recip a = case width a == zero of
      True  -> theta
      False -> Range (m - r/two) (m + r/two)
        where
          m = negate (mid a)
          r = recip (width a)

instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeCommutative (Range a)
instance (FromInteger a, Ord a, BoundedField a) => Multiplicative (Range a)
instance (FromInteger a, Ord a, BoundedField a) => MultiplicativeGroup (Range a)

instance (FromInteger a, AdditiveInvertible a, BoundedField a, Ord a) => Signed (Range a) where
    sign (Range l u) = if u >= l then one else Range half (negate half)
    abs (Range l u) = if u >= l then Range l u else Range u l

instance (AdditiveGroup a) => Normed (Range a) a where
    size (Range l u) = u-l

instance (Ord a, AdditiveGroup a) => Metric (Range a) a where
    distance (Range l u) (Range l' u')
        | u < l' = l' - u
        | u' < l = l - u'
        | otherwise = zero

instance (BoundedField a, Ord a, FromInteger a, Epsilon a) => Epsilon (Range a) where
    nearZero (Range l u) = nearZero (l - u)
    aboutEqual (Range l u) (Range l' u')= aboutEqual l l' && aboutEqual u u'

-- | theta is a bit like 1/infinity
theta :: (AdditiveUnital a) => Range a
theta = Range zero zero

instance (FromInteger a, Ord a, BoundedField a) => Space (Range a) where
    type Element (Range a) = a
    union (Range l0 u0) (Range l1 u1) = Range (min l0 l1) (max u0 u1)
    nul = Range infinity neginfinity
    lower (Range l _) = l
    upper (Range _ u) = u
    singleton a = Range a a
    type Grid (Range a) = Int
    grid :: FromInteger a => Pos -> Range a -> Int -> [a]
    grid o s n = (+ if o==MidPos then step/(one+one) else zero) <$> posns
      where
        posns = (lower s +) . (step *) . fromIntegral <$> [i0..i1]
        step = (/) (width s) (fromIntegral n)
        (i0,i1) = case o of
                    OuterPos -> (zero,n)
                    InnerPos -> (one,n - one)
                    LowerPos -> (zero,n - one)
                    UpperPos -> (one,n)
                    MidPos -> (zero,n - one)
    gridSpace r n = zipWith Range ps (drop 1 ps)
      where
        ps = grid OuterPos r n

instance (Ord a, BoundedField a, FromInteger a) => Monoid (Range a) where
    mempty = nul
    mappend = union

-- | turn a range into n `a`s pleasing to human sense and sensibility
-- the `a`s may well lie outside the original range as a result
gridSensible :: (Fractional a, Ord a, FromInteger a, QuotientField a, ExpField a) =>
    Pos -> Range a -> Int -> [a]
gridSensible tp (Range l u) n =
    (+ if tp==MidPos then step/two else zero) <$> posns
  where
    posns = (first' +) . (step *) . fromIntegral <$> [i0..i1]
    span = u - l
    step' = 10 ^^ floor (logBase 10 (span/fromIntegral n))
    err = fromIntegral n / span * step'
    step
      | err <= 0.15 = 10 * step'
      | err <= 0.35 = 5 * step'
      | err <= 0.75 = 2 * step'
      | otherwise = step'
    first' = step * fromIntegral (ceiling (l/step))
    last' = step * fromIntegral (floor (u/step))
    n' = round ((last' - first')/step)
    (i0,i1) = case tp of
                OuterPos -> (0,n')
                InnerPos -> (1,n' - 1)
                LowerPos -> (0,n' - 1)
                UpperPos -> (1,n')
                MidPos -> (0,n' - 1)