{-|
Module      : Data.Interval
Description : Closed intervals of totally ordered types, e.g. time intervals
Copyright   : (c) Lackmann Phymetric
License     : GPL-3
Maintainer  : olaf.klinke@phymetric.de
Stability   : experimental

This module provides the two-parameter type class 'Interval' of types 
that represent closed intervals (meaning the end-points are included) 
possibly with some extra annotation. 
This approach is shared by the Data.IntervalMap.Generic.Interval module of the 
<https://hackage.haskell.org/package/IntervalMap IntervalMap> package. 
A particular use case are time intervals annotated with event data. 
The simplest example of an interval type @i@ with end points of type @e@ 
is the type @i = (e,e)@. 

The functions exported from this module are mainly concerned with overlap queries, 
that is, to identify which intervals in a collection overlap a given interval 
and if so, to what extent. 
This functionality is encapsuled in the class 'IntersectionQuery'. 
If the collection of intervals is known to overlap in end-points only, 
one can simply use a sequence ordered by left end-point as the search structure. 
For arbitrary collections we provide the 'ITree' structure 
(centered interval tree) which stores intervals in subtrees and bins 
that are annotated with their convex hull, so that it can be decided 
easily whether there is an interval inside which overlaps a given interval. 


The behaviour of the functions is undefined for intervals that 
violate the implicit assumption that the left end-point is less than or equal to 
the right end-point. 

The functionality provided is similar to the Interval data type in the  
<https://hackage.haskell.org/package/data-interval data-interval> package 
but we focus on closed intervals and let the user decide which 
concrete data type to use. 

Most functions are property-checked for correctness. 
Checks were implemented by Henning Thielemann. 
-}
{-# LANGUAGE FlexibleInstances,FlexibleContexts,FunctionalDependencies,MultiParamTypeClasses,CPP #-}
module Data.Interval (
    -- * Type classes
    Interval(..),
    IntersectionQuery(..),
    Adjust(..),
    TimeDifference(..),
    -- * Comparing intervals
    intersects,properlyIntersects,contains,properlyContains,
    covered,coveredBy,overlapTime,prevailing,fractionCovered,
    overlap,
    intervalDuration,
    -- * Operations on intervals
    maybeUnion,maybeIntersection,
    hull,
    hullSeq,
    without,
    contiguous,components,componentsSeq,
    sortByLeft,
    fromEndPoints,
    -- * Interval search tree
    ITree,
    itree,
    emptyITree,
    insert,
    hullOfTree,
    intersecting,getIntersectsIT,getProperIntersectsIT,
    someIntersectsIT,someProperlyIntersectsIT,
    leftmostInterval,
    -- * Non-overlapping intervals
    findSeq, existsSeq, hullSeqNonOverlap,
    -- * Debug
    invariant, toTree,
    -- * Testing
    intersectingProperly,
    filterM,
    joinSeq,
    splitSeq,
    ) where

import Data.Tree (Tree)
import qualified Data.Tree as Tree
import qualified Data.Sequence as Seq
import qualified Data.Monoid ((<>))
import Data.Traversable (Traversable)
import Data.Foldable (toList, maximumBy, foldl', foldr')
import Data.Sequence (Seq, ViewL(EmptyL,(:<)), ViewR(EmptyR,(:>)), (><), (<|))
import Data.Function (on)
import Data.Functor.Identity (Identity(Identity, runIdentity))
import Data.Maybe (catMaybes)
import Data.Time (UTCTime, addUTCTime, diffUTCTime, utc, NominalDiffTime)
#if MIN_VERSION_time(1,9,0)
import Data.Time (LocalTime, utcToLocalTime, zonedTimeToLocalTime, diffLocalTime, addLocalTime)
#else
import Data.Time (LocalTime, utcToLocalTime, zonedTimeToLocalTime)
#endif
import Data.Time (ZonedTime, localTimeToUTC, zonedTimeToUTC)
import Control.Arrow ((***))
import Control.Applicative (Alternative, empty, (<|>))


-- $setup
-- >>> import Data.IntervalTest
-- >>> import qualified Data.Sequence as Seq
-- >>> import qualified Data.List as List
-- >>> import Data.Function (on)
-- >>> import Data.Maybe (isJust, fromJust, catMaybes)
-- >>> import Data.Foldable (toList)
-- >>> import qualified Test.QuickCheck as QC
-- >>> import Test.QuickCheck ((==>))
-- >>> without' :: (Int,Int) -> (Int,Int) -> [(Int,Int)]; without' = without


-- | class of intervals with end points in a totally ordered type
class (Ord e) => Interval e i | i -> e where
    lb :: i -> e -- ^ lower bound
    lb = fst.endPoints
    ub :: i -> e -- ^ upper bound
    ub = snd.endPoints
    endPoints :: i -> (e,e) -- ^ end points (inclusive)
    endPoints i = (lb i,ub i)
    {-# MINIMAL (lb,ub) | endPoints #-}

instance (Ord e) => Interval e (e,e) where
    endPoints = id

instance (Ord e) => Interval e (Identity e) where
    lb = runIdentity
    ub = runIdentity

-- | class of search structures for interval intersection queries,
-- returning a 'Foldable' of intervals.
class Foldable f => IntersectionQuery t e f | t -> f where
    getIntersects :: (Interval e i, Interval e j) => i -> t j -> f j
    -- ^ all intervalls touching the first one
    getProperIntersects :: (Interval e i, Interval e j) => i -> t j -> f j
    -- ^ all intervals properly intersecting the first one
    someIntersects :: (Interval e i, Interval e j) => i -> t j -> Bool
    -- ^ does any interval touch the first one?
    someProperlyIntersects :: (Interval e i, Interval e j) => i -> t j -> Bool
    -- ^ does any interval properly intersect the first one?
instance Ord e => IntersectionQuery (ITree e) e Seq where
    getIntersects = getIntersectsIT
    getProperIntersects = getProperIntersectsIT
    someIntersects = someIntersectsIT
    someProperlyIntersects = someProperlyIntersectsIT
instance Ord e => IntersectionQuery Seq e Seq where
    getIntersects = findSeq intersects
    getProperIntersects = findSeq properlyIntersects
    someIntersects = existsSeq intersects
    someProperlyIntersects = existsSeq properlyIntersects

-- | Time types supporting differences
class TimeDifference t where
    diffTime :: t -> t -> NominalDiffTime
    addTime :: NominalDiffTime -> t -> t
instance TimeDifference UTCTime where
    diffTime = diffUTCTime
    addTime  = addUTCTime
#if MIN_VERSION_time(1,9,0)
instance TimeDifference LocalTime where
    diffTime = diffLocalTime
    addTime  = addLocalTime
#else
instance TimeDifference LocalTime where
    diffTime x y = diffUTCTime (localTimeToUTC utc x) (localTimeToUTC utc y)
    addTime x = utcToLocalTime utc . addUTCTime x . localTimeToUTC utc
#endif
-- | 'addTime' preserves the 'TimeZone'
instance TimeDifference ZonedTime where
    diffTime x y = diffUTCTime (zonedTimeToUTC x) (zonedTimeToUTC y)
    addTime x z = z {zonedTimeToLocalTime = addTime x (zonedTimeToLocalTime z)}

-- | Convenience function, the 'diffTime' between the 'endPoints'.
intervalDuration :: (TimeDifference t, Interval t i) => i -> NominalDiffTime
intervalDuration i = diffTime (ub i) (lb i)

-- | Find out the overlap of two time intervals.
--
-- prop> genInterval /\ \i -> overlapTime i i == intervalDuration i
-- prop> genInterval /\* \i j -> not (i `properlyIntersects` j) ==> overlapTime i j == 0
-- prop> genInterval /\* \i j -> overlapTime i j == (sum $ fmap intervalDuration $ maybeIntersection i j)
overlapTime :: (TimeDifference t, Interval t i, Interval t j) =>
    i -> j -> NominalDiffTime
overlapTime i j = let
    x = max (lb i) (lb j)
    y = min (ub i) (ub j)
    in if x < y then diffTime y x else 0

-- | Prevailing annotation in the first time interval
--
-- prop> genInterval /\ \i c -> prevailing i (Seq.singleton (c,i)) == Just (c::Char)
-- prop> genInterval /\ \i -> genLabeledSeq /\ \js -> isJust (prevailing i js) == any (intersects i . snd) js
-- prop> genInterval /\ \i -> genLabeledSeq /\* \js ks -> all (flip elem $ catMaybes [prevailing i js, prevailing i ks]) $ prevailing i (js<>ks)
prevailing :: (Interval t i, Interval t j, TimeDifference t) =>
    i -> Seq (a,j) -> Maybe a
prevailing i js =
    let ks = Seq.filter (intersects i . snd) js
    in  if Seq.null ks
            then Nothing
            else Just $ fst $ maximumBy (compare `on` (overlapTime i . snd)) ks
            -- ExtPkg: non-empty - partial maximumBy -> NonEmpty.maximumBy


-- | class of Intervals whose bounds can be adjusted
class Interval e i => Adjust e i | i -> e where
    adjustBounds :: (e -> e) -> (e -> e) -> i -> i -- ^  adjust lower and upper bound
    shift :: (e -> e) -> i -> i -- ^ change both bounds using the same function
    shift f = adjustBounds f f
    {-# MINIMAL (adjustBounds) #-}

instance Ord e => Adjust e (e,e) where
    adjustBounds f g (x,y) = (f x,g y)

-- | the union of two intervals is an interval if they intersect.
--
-- prop> genInterval /\* \i j -> isJust (maybeUnion i j) ==> fromJust (maybeUnion i j) `contains` i && fromJust (maybeUnion i j) `contains` j
-- prop> genInterval /\* \i j -> i `intersects` j ==> (maybeUnion i j >>= maybeIntersection i) == Just i
maybeUnion :: (Interval e j, Interval e i, Adjust e i) => j -> i -> Maybe i
maybeUnion j i = if j `intersects` i
    then Just (adjustBounds (min (lb j)) (max (ub j)) i)
    else Nothing

-- | the intersection of two intervals is either empty or an interval.
--
-- prop> genInterval /\* \i j -> i `intersects` j ==> i `contains` fromJust (maybeIntersection i j)
maybeIntersection :: (Interval e j, Interval e i, Adjust e i) => j -> i -> Maybe i
maybeIntersection j i = if j `intersects` i
    then Just (adjustBounds (max (lb j)) (min (ub j)) i)
    else Nothing

-- | convex hull
--
-- prop> \xs -> isJust (hull xs) ==> all (\x -> fromJust (hull xs) `contains` x) (xs :: [(Int,Int)])
hull :: (Interval e i,Foldable f,Functor f) => f i -> Maybe (e,e)
hull xs = if null xs
    then Nothing
    else Just (minimum (fmap lb xs), maximum (fmap ub xs))

-- | Set difference. The resulting list has zero, one or two elements.
--
-- >>> without' (1,5) (4,5)
-- [(1,4)]
-- >>> without' (1,5) (2,3)
-- [(1,2),(3,5)]
-- >>> without' (1,5) (1,5)
-- []
-- >>> without' (1,5) (0,1)
-- [(1,5)]
--
-- prop> genInterval /\* \i j -> length (i `without` j) <= 2
-- prop> genInterval /\ \i -> i `without` i == []
-- prop> genInterval /\* \i j -> all (contains i) (i `without` j)
-- prop> genInterval /\* \i j -> not $ any (properlyIntersects j) (i `without` j)
without :: (Adjust e i,Interval e j) => i -> j -> [i]
without i j = if j `contains` i then [] else
    if ub j <= lb i || lb j >= ub i
        then [i] -- intervals don't overlap
        else if i `properlyContains` j
            then [adjustBounds id (const (lb j)) i,adjustBounds (const (ub j)) id i] -- slashed in half
            else if lb j <= lb i
                then [adjustBounds (const (ub j)) id i] -- j overhangs on the left
                else [adjustBounds id (const (lb j)) i] -- j overhangs on the right

-- | 'intersects' is not an equivalence relation, because it is not transitive.
-- Hence 'groupBy' 'intersects' does not do what one might expect.
-- This function does the expected and groups overlapping intervals
-- into contiguous blocks.
--
-- prop> genSortedIntervals /\ all (\xs -> and $ List.zipWith intersects xs (tail xs)) . contiguous
contiguous :: Interval e i => [i] -> [[i]]
contiguous [] = []
contiguous (i:is) = (i:js) : contiguous ks where
    (js,ks) = go (endPoints i) is
    go :: Interval e i => (e,e) -> [i] -> ([i],[i])
    go j@(x,_y) ls@(l:ls') = if j `intersects` l
        then let (foo,bar) = go (x,ub l) ls' in (l:foo,bar)
        else ([],ls)
    go _ [] = ([],[])

-- | Connected components of a list sorted by 'sortByLeft',
-- akin to 'groupBy' 'intersects'.
-- The precondition is not checked.
--
-- prop> genSortedIntervals /\ \xs -> all (\i -> any (flip contains i) (components xs)) xs
components :: (Interval e i, Adjust e i) => [i] -> [i]
components [] = []
components (i:is) = c i is where
    c x [] = [x]
    c x (y:ys) = case maybeUnion x y of
        Nothing -> x : c y ys
        Just z  -> c z ys

-- | same as 'components'. Is there a way to unify both?
--
-- prop> genSortedIntervals /\ \xs -> componentsSeq (Seq.fromList xs) == Seq.fromList (components xs)
componentsSeq :: (Interval e i, Adjust e i) => Seq i -> Seq i
componentsSeq ys = case Seq.viewl ys of
    EmptyL  -> empty
    x :< xs -> c x xs where
        c a bs = case Seq.viewl bs of
            EmptyL  -> Seq.singleton a
            b :< bs' -> case maybeUnion a b of
                Nothing -> a <| c b bs'
                Just ab -> c ab bs'

-- | compute the components of the part of @i@ covered by the intervals.
--
-- prop> genInterval /\ \i -> genIntervalSeq /\ \js -> all (contains i) (covered i js)
-- prop> genInterval /\ \i -> genIntervalSeq /\ \js -> covered i (covered i js) == covered i js
covered :: (Interval e i,Interval e j,Adjust e j) => i -> Seq j -> Seq j
covered i =
    let mapMaybe f = foldMap (foldMap Seq.singleton . f)
    in  componentsSeq . sortByLeft . mapMaybe (maybeIntersection i)

-- | 'True' if the first interval is completely covered by the given intervals
--
-- prop> genInterval /\* \i j -> j `contains` i == i `coveredBy` [j]
-- prop> genInterval /\ \i -> genSortedIntervals /\ \js -> i `coveredBy` js ==> any (flip contains i) (components js)
coveredBy :: (Interval e i, Interval e j, Foldable f) => i -> f j -> Bool
i `coveredBy` js = null $ foldl (\remains j -> flip without j =<< remains) [endPoints i] js

-- | percentage of coverage of the first interval by the second sequence of intervals
--
-- prop> genNonEmptyInterval /\ \i -> genIntervalSeq /\ \js -> i `coveredBy` js == (fractionCovered i js >= (1::Rational))
-- prop> genNonEmptyInterval /\ \i -> genNonEmptyIntervalSeq /\ \js -> any (properlyIntersects i) js == (fractionCovered i js > (0::Rational))
fractionCovered :: (TimeDifference t, Interval t i, Interval t j, Fractional a) =>
    j -> Seq i -> a
fractionCovered i xs = let
    totalTime   = intervalDuration i
    coveredTime = foldl' (\s j -> s + intervalDuration j) 0 $ covered i $ fmap endPoints xs
    -- ^ sum of the lengths of the interections with i
    in if totalTime==0 then 1 else (fromRational.toRational) (coveredTime/totalTime) -- (fromInteger (round (coveredTime*100/totalTime)))/100

-- | Overlap ordering. Returns 'LT' or 'GT' if the intervals are disjoint,
-- 'EQ' if the intervals overlap.
-- Note that this violates the following property:
--
-- @
-- 'overlap' x y == 'EQ' && 'overlap' y z == 'EQ' => 'overlap' x z == 'EQ'
-- @
--
-- i.e., 'overlap' is not transitive.
--
-- prop> genInterval /\* \i j -> i `intersects` j  ==  (overlap i j == EQ)
overlap :: Interval e i => i -> i -> Ordering
overlap i j = case (compare (ub i) (lb j),compare (ub j) (lb i)) of
    (LT,_) -> LT
    (_,LT) -> GT
    _      -> EQ

-- | intersection query.
--
-- >>> ((1,2)::(Int,Int)) `intersects` ((2,3)::(Int,Int))
-- True
--
-- prop> genInterval /\* \i j -> (lb i <= ub i && lb j <= ub j && i `intersects` j)  ==  (max (lb i) (lb j) <= min (ub i) (ub j))
intersects :: (Interval e i,Interval e j) => i -> j -> Bool
i `intersects` j = not (ub i < lb j || ub j < lb i)
-- The definition of 'intersects' yields the following algorithm
-- for intersection queries.
-- Given the query interval i, sort the list of possible intersecting intervals
-- by 'ub' and consider the suffix of intervals j with lb i <= ub j.
-- Sort that suffix by 'lb' and take the prefix with lb j <= ub i.

-- | proper intersection.
--
-- prop> genInterval /\* \i j -> ((i `intersects` j) && not (i `properlyIntersects` j))  ==  (ub i == lb j || ub j == lb i)
properlyIntersects :: (Interval e i,Interval e j) => i -> j -> Bool
i `properlyIntersects` j = not (ub i <= lb j || ub j <= lb i)

-- | subset containment
--
-- prop> genInterval /\ \i -> i `contains` i
-- prop> genInterval /\* \i j -> (i `contains` j && j `contains` i) == (i==j)
-- prop> genInterval /\* \i j -> i `contains` j == (maybeUnion i j == Just i)
contains :: (Interval e i,Interval e j) => i -> j -> Bool
i `contains` j = lb i <= lb j && ub j <= ub i

-- | proper subset containment
properlyContains :: (Interval e i,Interval e j) => i -> j -> Bool
i `properlyContains` j = lb i < lb j && ub i > ub j

-- | construct a sorted sequence of intervals
-- from a sorted sequence of bounds.
-- Fails if the input sequence is not sorted.
--
-- prop> genSortedList /\ \xs -> (components $ toList $ fromEndPoints xs) == if length xs < 2 then [] else [(head xs, last xs)]
fromEndPoints :: (Ord e) => [e] -> Seq (e,e)
fromEndPoints [] = empty
fromEndPoints [_] = empty
fromEndPoints [x,y] = if x <= y then Seq.singleton (x,y) else error "unsorted list"
fromEndPoints (x:xs) = let s  = fromEndPoints xs in case Seq.viewl s of
    (y,_) :< _ -> (x,y) <| s
    EmptyL     -> error "Intervals.fromEndPoints: this should never happen"

-- | lexicographical sort by 'lb', then inverse 'ub'.
-- In the resulting list, the intervals intersecting
-- a given interval form a contiguous sublist.
--
-- prop> genInterval /\ \i -> genSortedIntervalSeq /\ \js -> toList (getIntersects i js) `List.isSubsequenceOf` toList js
sortByLeft :: (Interval e i) => Seq i -> Seq i
sortByLeft = Seq.sortBy (\i j -> compare (lb i) (lb j) <> compare (ub j) (ub i))

-- | extract all intervals intersecting a given one.
intersecting :: (Interval e i,Interval e j) => j -> Seq i -> Seq i
intersecting j = Seq.filter (intersects j)

-- | extract all intervals properly intersecting a given one.
intersectingProperly :: (Interval e i,Interval e j) => j -> Seq i -> Seq i
intersectingProperly j = Seq.filter (properlyIntersects j)
--intersectingProperly j = (takeWhileL (properlyIntersects j)).(dropWhileL (not.(properlyIntersects j)))

-- | convex hull of a sorted sequence of intervals.
-- the lower bound is guaranteed to be in the leftmost interval,
-- but we have no guarantee of the upper bound.
--
-- prop> genSortedIntervalSeq /\ \xs -> hullSeq xs == hull (toList xs)
hullSeq :: Interval e i => Seq i -> Maybe (e,e)
hullSeq xs = case Seq.viewl xs of
    EmptyL -> Nothing
    leftmost :< _others -> Just (lb leftmost, maximum (fmap ub xs))

-- | Search tree of intervals.
data ITree e i = Bin (Seq i) | Split (Seq i) e e e (ITree e i) (ITree e i)
-- Internal nodes store the convex hull of its subtrees.
-- Each bin contains a sorted sequence of intervals.
-- In the node @Split top x y z left right@
-- the convex hull of @left@ is @(x,y)@,
-- the convex hull of @right@ is @(y,z)@
-- and the intervals in @top@ are those straddling the split point @y@.
instance Functor (ITree e) where
    fmap f (Bin xs) = Bin (fmap f xs)
    fmap f (Split up x y z left right) = Split (fmap f up) x y z (fmap f left) (fmap f right)
instance Foldable (ITree e) where
    foldMap f (Bin xs) = foldMap f xs
    foldMap f (Split up _ _ _ left right) = foldMap f left <> foldMap f up <> foldMap f right

-- | the empty 'ITree'
emptyITree :: ITree e i
emptyITree = Bin empty

-- | smallest interval covering the entire tree. 'Nothing' if the tree is empty.
hullOfTree :: (Interval e i) => ITree e i -> Maybe (e,e)
hullOfTree (Bin xs) = hullSeq xs
hullOfTree (Split _ x _ y _ _) = Just (x,y)

-- | invariant to be maintained for proper intersection queries
--
-- prop> invariant . itree 4 . fmap (\(x,y) -> (x, x + QC.getNonNegative y :: Integer))
invariant :: Interval e i => ITree e i -> Bool
invariant (Bin _) = True
invariant (Split up x y z left right) = x <= y && y <= z && invUp && invLeft && invRight where
    invUp = all (intersects (Identity y)) up && all (contains (x,z)) up
    invLeft = all (contains (x,y)) left && invariant left
    invRight = all (contains (y,z)) right && invariant right

-- | Intersection query. O(binsize+log(n/binsize)).
--
-- prop> genInterval /\ \i -> genIntervalSeq /\ \t -> on (==) sortByLeft (getIntersectsIT i $ itree 2 t) (i `intersecting` t)
getIntersectsIT :: (Interval e i, Interval e j) => i -> ITree e j -> Seq j
getIntersectsIT i (Bin bin) = i `intersecting` bin
getIntersectsIT i (Split up x y z left right) = let
    m = i `intersecting` up
    l = if i `intersects` (x,y) then getIntersectsIT i left else empty
    r = if i `intersects` (y,z) then getIntersectsIT i right else empty
    in if i `intersects` (x,z) then m >< l >< r else empty

-- | Intersection query. O(binsize+log(n/binsize)).
--
-- prop> genInterval /\ \i -> genIntervalSeq /\ \t -> on (==) sortByLeft (getProperIntersectsIT i $ itree 2 t) (i `intersectingProperly` t)
getProperIntersectsIT :: (Interval e i, Interval e j) => i -> ITree e j -> Seq j
getProperIntersectsIT i (Bin bin) = i `intersectingProperly` bin
getProperIntersectsIT i (Split up x y z left right) = let
    m = i `intersectingProperly` up
    l = if i `properlyIntersects` (x,y) then getProperIntersectsIT i left else empty
    r = if i `properlyIntersects` (y,z) then getProperIntersectsIT i right else empty
    in if i `properlyIntersects` (x,z) then m >< l >< r else empty

-- | When the actual result of 'getIntersectsIT' is not important,
-- only whether there are intersections.
someIntersectsIT :: (Interval e i, Interval e j) => i -> ITree e j -> Bool
someIntersectsIT i = not . null . getIntersectsIT i

-- | When the actual result of 'getIntersectsIT' is not important,
-- only whether there are intersections.
someProperlyIntersectsIT :: (Interval e i, Interval e j) => i -> ITree e j -> Bool
someProperlyIntersectsIT i = not . null . getProperIntersectsIT i

-- | retrieve the left-most interval from the tree, or 'Nothing' if it is empty.
leftmostInterval :: (Interval e i) => ITree e i -> Maybe i
leftmostInterval (Bin bin) = case Seq.viewl bin of
    EmptyL -> Nothing
    i :< _ -> Just i
leftmostInterval (Split up _ _ _ left right) = let
    headl xs = case Seq.viewl xs of
        EmptyL -> Nothing
        i :< _ -> Just i
    in (headl . sortByLeft . Seq.fromList . catMaybes) [leftmostInterval left,headl up,leftmostInterval right]

-- | transform the interval tree into the tree of hulls
toTree :: Interval e i => ITree e i -> Tree (e,e)
toTree (Bin _) = error "Interval.toTree: just a bin"
toTree (Split _ x y z left right) = Tree.Node {Tree.rootLabel = (x,z), Tree.subForest = [l,r]} where
    l = case left of
        (Bin _) -> Tree.Node {Tree.rootLabel = (x,y), Tree.subForest = []}
        _ -> toTree left
    r = case right of
        (Bin _) -> Tree.Node {Tree.rootLabel = (y,z), Tree.subForest = []}
        _ -> toTree right

-- ExtPkg: non-empty allows NonEmpty Seq - makes blockstart total
newtype Block e i = Block (Seq i)
blockstart :: Interval e i => Block e i -> e
blockstart (Block xs) = case Seq.viewl xs of
    EmptyL -> error "empty Block"
    x :< _ -> lb x
blocknull :: Block e i -> Bool
blocknull (Block xs) = null xs
instance Interval e i => Eq (Block e i) where
    (==) = (==) `on` blockstart
instance Interval e i => Ord (Block e i) where
    compare = compare `on` (blockstart)
instance Functor (Block e) where
    fmap f (Block xs) = Block (fmap f xs)
instance Foldable (Block e) where
    foldMap f (Block xs) = foldMap f xs
instance Semigroup (Block e i) where
    (Block xs) <> (Block ys) = Block (xs >< ys)
instance Monoid (Block e i) where
    mempty = Block empty
    mappend = (<>)
instance Show i => Show (Block e i) where
    show (Block xs) = "Block "++(show (toList xs))

-- | generalises Control.Monad.filterM
filterM :: (Applicative f, Traversable t, Alternative m) => (a -> f Bool) -> t a -> f (m a)
filterM f = (fmap (foldr (<|>) empty)) . traverse (\a -> fmap (\b -> if b then pure a else empty) (f a))

crossesAny :: (Interval e i, Foldable f) => i -> f (Block e i) -> Bool
crossesAny i = any (((ub i) >).blockstart)

removeCrossers :: Interval e i => Block e i -> Seq (Block e i) -> (Seq i,Block e i)
removeCrossers (Block xs) blocks = let (crossers,xs') = filterM f xs in (crossers,Block xs') where
    f i = if i `crossesAny` blocks
        then (Seq.singleton i,False)
        else return True

-- foldr over the list of blocks and gather all intervals
-- overlapping block boundaries. Remove blocks that are rendered empty by this.
gatherCrossers :: Interval e i => Seq (Block e i) -> (Seq i,Seq (Block e i))
gatherCrossers blks = case Seq.viewl blks of
    (block :< blocks) -> let
        (crossers,blocks') = gatherCrossers blocks
        (crossers',block') = removeCrossers block blocks'
        cons = if blocknull block' then id else ((<|) block')
        in (crossers' >< crossers,cons blocks')
    EmptyL -> (empty,empty)
-- after applying gatherCrossers to a sorted list of sorted blocks,
-- all intervals within the blocks are contained in the interval
-- from the blockstart to the blockstart of the next block.
-- Hence we can use these blocks to build an interval tree,
-- where the crossers go into certain 'up' components.

blocksOf :: Int -> Seq i -> Seq (Block e i)
blocksOf n = fmap Block . Seq.chunksOf n

data SplitSeq a = EmptySeq | SingletonSeq a | TwoSeqs (Seq a) (Seq a) deriving (Show)

joinSeq :: SplitSeq a -> Seq a
joinSeq EmptySeq = empty
joinSeq (SingletonSeq a) = pure a
joinSeq (TwoSeqs xs ys) = xs <> ys

-- |
-- prop> genIntervalSeq /\ \is -> joinSeq (splitSeq is) == is
splitSeq :: Seq a -> SplitSeq a
splitSeq xs = let (l,r) = Seq.splitAt (length xs `div` 2) xs in case (null l,null r) of
    (_,True) -> EmptySeq
    (True,False) -> let (x :< _) = Seq.viewl r in SingletonSeq x
    (False,False) -> TwoSeqs l r

-- build a tree from a sequence of mutually non-overlapping blocks
buildFromSeq :: Interval e i => Seq (Block e i) -> ITree e i
buildFromSeq blocks = case splitSeq blocks of
    EmptySeq -> emptyITree
    SingletonSeq (Block bin) -> Bin bin
    TwoSeqs lblocks rblocks -> let
        y = let b :< _ = Seq.viewl rblocks in blockstart b
        left = buildFromSeq lblocks
        right = buildFromSeq rblocks
        x = maybe y fst (hullOfTree left)
        z = maybe y snd (hullOfTree right)
        in Split empty x y z left right

-- | insert the interval at the deepest possible location into the tree.
-- Does not change the overall structure, in particular no re-balancing is performed.
insert :: Interval e i => i -> ITree e i -> ITree e i
insert i (Bin xs) = Bin (i <| xs)
insert i (Split up x y z left right) = if ub i <= y
    then let
        left' = (insert i left)
        x' = maybe x (min x.fst) (hullOfTree left')
        in Split up x' y z left' right
    else if lb i < y
        then Split (i <| up) (min x (lb i)) y (max z (ub i)) left right
        else let
            right' = insert i right
            z' = maybe z (max z.snd) (hullOfTree right')
            in Split up x y z' left right'

-- | Construct an interval tree with bins of maximal given size.
-- The function first sorts the intervals,
-- then splits into chunks of given size.
-- The leftmost endpoints of the chunks define boundary points.
-- Next, all intervals properly overlapping a boundary are removed
-- from the chunks and kept separately.
-- The chunks are arranged as the leaves of a binary search tree.
-- Then the intervals overlapping boundaries are placed
-- at internal nodes of the tree.
-- Hence if all intervals are mutually non-overlapping properly,
-- the resulting tree is a pure binary search tree with bins of
-- given size as leaves.
itree :: Interval e i => Int -> Seq i -> ITree e i
itree n = uncurry ($).(f *** buildFromSeq).gatherCrossers.blocksOf n.srt where
    srt = Seq.unstableSortBy (compare `on` endPoints)
    f = flip (foldr' insert)

-- * Non-overlapping intervals

-- | /O(1)/ bounds of an ordered sequence of intervals. 'Nothing', if empty.
--
-- prop> genDisjointIntervalSeq /\ \xs -> hullSeqNonOverlap xs == hullSeq xs
hullSeqNonOverlap :: Interval e i => Seq i -> Maybe (e,e)
hullSeqNonOverlap xs = case Seq.viewl xs of
    EmptyL -> Nothing
    leftmost :< others -> Just (lb leftmost, case Seq.viewr others of
        _ :> rightmost -> ub rightmost
        EmptyR         -> ub leftmost)

-- | Query an ordered 'Seq'uence of non-overlapping intervals
-- for a predicate @p@ that has the property
--
-- @
-- j `contains` k && p i k ==> p i j
-- @
--
-- and return all elements satisfying the predicate.
--
-- prop> genInterval /\ \i -> genDisjointIntervalSeq /\ \js -> findSeq intersects i js == intersecting i js
findSeq :: (Interval e i, Interval e j) => (i -> (e,e) -> Bool) -> i -> Seq j -> Seq j
findSeq p i js = case hullSeqNonOverlap js of
    Nothing -> empty
    Just h -> if p i h
        then case splitSeq js of
            SingletonSeq _j -> js
            TwoSeqs l r -> findSeq p i l >< findSeq p i r
            EmptySeq -> empty -- should never happen
        else empty

-- | Query an ordered 'Seq'uence of non-overlapping intervals
-- for a predicate @p@ that has the property
--
-- @
-- j `contains` k && p i k ==> p i j
-- @
existsSeq :: (Interval e i, Interval e j) => (i -> (e,e) -> Bool) -> i -> Seq j -> Bool
existsSeq p i js = case hullSeqNonOverlap js of
    Nothing -> False
    Just h -> if p i h
        then case splitSeq js of
            SingletonSeq _j -> True
            TwoSeqs l r -> existsSeq p i l || existsSeq p i r
            EmptySeq -> False -- should never happen
        else False