--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.RIC
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- An randomized algorithm to compute the smallest enclosing disk of a set of
-- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\).
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.SmallestEnclosingBall.RIC(
    smallestEnclosingDisk'
  , smallestEnclosingDisk
  , smallestEnclosingDiskWithPoint
  , smallestEnclosingDiskWithPoints
  ) where

import           Algorithms.Geometry.SmallestEnclosingBall.Types
import           Control.Lens
import           Control.Monad.Random.Class
import           Data.Ext
import qualified Data.Foldable as F
import           Data.Geometry
import           Data.Geometry.Ball
import qualified Data.List as List
import           Data.List.NonEmpty(NonEmpty(..))
import           Data.Maybe (fromMaybe, mapMaybe, catMaybes)
import           Data.Ord (comparing)
import           System.Random.Shuffle (shuffle)

-- import Data.RealNumber.Rational
-- import Debug.Trace

--------------------------------------------------------------------------------

-- | Compute the smallest enclosing disk of a set of points,
-- implemented using randomized incremental construction.
--
-- pre: the input has at least two points.
--
-- running time: expected \(O(n)\) time, where \(n\) is the number of input points.
smallestEnclosingDisk           :: (Ord r, Fractional r, MonadRandom m
                                               -- , Show r, Show p
                                   )
                                => [Point 2 r :+ p]
                                -> m (DiskResult p r)

smallestEnclosingDisk :: [Point 2 r :+ p] -> m (DiskResult p r)
smallestEnclosingDisk pts :: [Point 2 r :+ p]
pts@(Point 2 r :+ p
_:Point 2 r :+ p
_:[Point 2 r :+ p]
_) = (\(Point 2 r :+ p
p:Point 2 r :+ p
q:[Point 2 r :+ p]
pts') -> (Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> DiskResult p r
forall r p.
(Ord r, Fractional r) =>
(Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> DiskResult p r
smallestEnclosingDisk' Point 2 r :+ p
p Point 2 r :+ p
q [Point 2 r :+ p]
pts')
                                    ([Point 2 r :+ p] -> DiskResult p r)
-> (Vector (Point 2 r :+ p) -> [Point 2 r :+ p])
-> Vector (Point 2 r :+ p)
-> DiskResult p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Point 2 r :+ p) -> [Point 2 r :+ p]
forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList (Vector (Point 2 r :+ p) -> DiskResult p r)
-> m (Vector (Point 2 r :+ p)) -> m (DiskResult p r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Point 2 r :+ p] -> m (Vector (Point 2 r :+ p))
forall (f :: * -> *) (m :: * -> *) a.
(Foldable f, MonadRandom m) =>
f a -> m (Vector a)
shuffle [Point 2 r :+ p]
pts
smallestEnclosingDisk [Point 2 r :+ p]
_           = [Char] -> m (DiskResult p r)
forall a. HasCallStack => [Char] -> a
error [Char]
"smallestEnclosingDisk: Too few points"

-- | Smallest enclosing disk.
smallestEnclosingDisk'     :: (Ord r, Fractional r
                                               -- , Show r, Show p

                              )
                           => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
                           -> DiskResult p r
smallestEnclosingDisk' :: (Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> DiskResult p r
smallestEnclosingDisk' Point 2 r :+ p
a Point 2 r :+ p
b = ([Point 2 r :+ p] -> DiskResult p r -> DiskResult p r)
-> DiskResult p r -> [[Point 2 r :+ p]] -> DiskResult p r
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr [Point 2 r :+ p] -> DiskResult p r -> DiskResult p r
addPoint ((Point 2 r :+ p) -> (Point 2 r :+ p) -> DiskResult p r
forall r p.
Fractional r =>
(Point 2 r :+ p) -> (Point 2 r :+ p) -> DiskResult p r
initial Point 2 r :+ p
a Point 2 r :+ p
b) ([[Point 2 r :+ p]] -> DiskResult p r)
-> ([Point 2 r :+ p] -> [[Point 2 r :+ p]])
-> [Point 2 r :+ p]
-> DiskResult p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Point 2 r :+ p] -> [[Point 2 r :+ p]]
forall a. [a] -> [[a]]
List.tails
  where
    -- The empty case occurs only initially
    addPoint :: [Point 2 r :+ p] -> DiskResult p r -> DiskResult p r
addPoint []      DiskResult p r
br            = DiskResult p r
br
    addPoint (Point 2 r :+ p
p:[Point 2 r :+ p]
pts) br :: DiskResult p r
br@(DiskResult Disk () r
d TwoOrThree (Point 2 r :+ p)
_)
      | (Point 2 r :+ p
p(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) Point 2 r -> Disk () r -> Bool
forall (d :: Nat) r p.
(Arity d, Ord r, Num r) =>
Point d r -> Ball d p r -> Bool
`inClosedBall` Disk () r
d = DiskResult p r
br
      | Bool
otherwise                  = Maybe (DiskResult p r) -> DiskResult p r
forall a. Maybe a -> a
fromJust' (Maybe (DiskResult p r) -> DiskResult p r)
-> Maybe (DiskResult p r) -> DiskResult p r
forall a b. (a -> b) -> a -> b
$ (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p) -> Maybe (DiskResult p r)
forall r p.
(Ord r, Fractional r) =>
(Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p) -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoint Point 2 r :+ p
p (Point 2 r :+ p
a (Point 2 r :+ p) -> [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p)
forall a. a -> [a] -> NonEmpty a
:| (Point 2 r :+ p
b (Point 2 r :+ p) -> [Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. a -> [a] -> [a]
: [Point 2 r :+ p]
pts))
    fromJust' :: Maybe a -> a
fromJust' = a -> Maybe a -> a
forall a. a -> Maybe a -> a
fromMaybe ([Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"smallestEncosingDisk' : fromJust, absurd")

-- | Smallest enclosing disk, given that p should be on it.
smallestEnclosingDiskWithPoint              :: (Ord r, Fractional r
                                               -- , Show r, Show p
                                               )
                                            => Point 2 r :+ p -> NonEmpty (Point 2 r :+ p)
                                            -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoint :: (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p) -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoint Point 2 r :+ p
p (Point 2 r :+ p
a :| [Point 2 r :+ p]
pts) = ([Point 2 r :+ p]
 -> Maybe (DiskResult p r) -> Maybe (DiskResult p r))
-> Maybe (DiskResult p r)
-> [[Point 2 r :+ p]]
-> Maybe (DiskResult p r)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr [Point 2 r :+ p]
-> Maybe (DiskResult p r) -> Maybe (DiskResult p r)
addPoint (DiskResult p r -> Maybe (DiskResult p r)
forall a. a -> Maybe a
Just (DiskResult p r -> Maybe (DiskResult p r))
-> DiskResult p r -> Maybe (DiskResult p r)
forall a b. (a -> b) -> a -> b
$ (Point 2 r :+ p) -> (Point 2 r :+ p) -> DiskResult p r
forall r p.
Fractional r =>
(Point 2 r :+ p) -> (Point 2 r :+ p) -> DiskResult p r
initial Point 2 r :+ p
p Point 2 r :+ p
a) ([[Point 2 r :+ p]] -> Maybe (DiskResult p r))
-> [[Point 2 r :+ p]] -> Maybe (DiskResult p r)
forall a b. (a -> b) -> a -> b
$ [Point 2 r :+ p] -> [[Point 2 r :+ p]]
forall a. [a] -> [[a]]
List.tails [Point 2 r :+ p]
pts
  where
    addPoint :: [Point 2 r :+ p]
-> Maybe (DiskResult p r) -> Maybe (DiskResult p r)
addPoint []       Maybe (DiskResult p r)
br   = Maybe (DiskResult p r)
br
    addPoint (Point 2 r :+ p
q:[Point 2 r :+ p]
pts') br :: Maybe (DiskResult p r)
br@(Just (DiskResult Disk () r
d TwoOrThree (Point 2 r :+ p)
_))
      | (Point 2 r :+ p
q(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) Point 2 r -> Disk () r -> Bool
forall (d :: Nat) r p.
(Arity d, Ord r, Num r) =>
Point d r -> Ball d p r -> Bool
`inClosedBall` Disk () r
d = Maybe (DiskResult p r)
br
      | Bool
otherwise                  = (Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> Maybe (DiskResult p r)
forall r p.
(Ord r, Fractional r) =>
(Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoints Point 2 r :+ p
p Point 2 r :+ p
q (Point 2 r :+ p
a(Point 2 r :+ p) -> [Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. a -> [a] -> [a]
:[Point 2 r :+ p]
pts')
    addPoint [Point 2 r :+ p]
_        Maybe (DiskResult p r)
br           = Maybe (DiskResult p r)
br


-- | Smallest enclosing disk, given that p and q should be on it
--
-- running time: \(O(n)\)
smallestEnclosingDiskWithPoints        :: (Ord r, Fractional r
                                          -- , Show r, Show p
                                          )
                                       => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
                                       -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoints :: (Point 2 r :+ p)
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoints Point 2 r :+ p
p Point 2 r :+ p
q [Point 2 r :+ p]
ps = (DiskResult p r -> r) -> [DiskResult p r] -> Maybe (DiskResult p r)
forall b a. Ord b => (a -> b) -> [a] -> Maybe a
minimumOn (DiskResult p r -> Getting r (DiskResult p r) r -> r
forall s a. s -> Getting a s a -> a
^.(Disk () r -> Const r (Disk () r))
-> DiskResult p r -> Const r (DiskResult p r)
forall p r. Lens' (DiskResult p r) (Disk () r)
enclosingDisk((Disk () r -> Const r (Disk () r))
 -> DiskResult p r -> Const r (DiskResult p r))
-> ((r -> Const r r) -> Disk () r -> Const r (Disk () r))
-> Getting r (DiskResult p r) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Disk () r -> Const r (Disk () r)
forall (d :: Nat) p r. Lens' (Ball d p r) r
squaredRadius)
                                       ([DiskResult p r] -> Maybe (DiskResult p r))
-> [DiskResult p r] -> Maybe (DiskResult p r)
forall a b. (a -> b) -> a -> b
$ [Maybe (DiskResult p r)] -> [DiskResult p r]
forall a. [Maybe a] -> [a]
catMaybes [Maybe ((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r)
mkEnclosingDisk Maybe ((Point 2 r :+ p) :+ Disk () r)
dl, Maybe ((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r)
mkEnclosingDisk Maybe ((Point 2 r :+ p) :+ Disk () r)
dr, Maybe (DiskResult p r)
mdc]
  where
    centers :: [(Point 2 r :+ p) :+ Disk () r]
centers = ((Point 2 r :+ p) -> Maybe ((Point 2 r :+ p) :+ Disk () r))
-> [Point 2 r :+ p] -> [(Point 2 r :+ p) :+ Disk () r]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe (Point 2 r :+ p) -> Maybe ((Point 2 r :+ p) :+ Disk () r)
disk' [Point 2 r :+ p]
ps
    -- generate a disk with p q and r
    disk' :: (Point 2 r :+ p) -> Maybe ((Point 2 r :+ p) :+ Disk () r)
disk' Point 2 r :+ p
r = (Point 2 r :+ p
r(Point 2 r :+ p) -> Disk () r -> (Point 2 r :+ p) :+ Disk () r
forall core extra. core -> extra -> core :+ extra
:+) (Disk () r -> (Point 2 r :+ p) :+ Disk () r)
-> Maybe (Disk () r) -> Maybe ((Point 2 r :+ p) :+ Disk () r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
forall r.
(Eq r, Fractional r) =>
Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
disk (Point 2 r :+ p
p(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
q(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
r(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)

    -- partition the points in to those on the left and those on the
    -- right.  Note that centers still contains only those points (and
    -- disks) for which the three points are not colinear. So the
    -- points are either on the left or on the right.
    ([(Point 2 r :+ p) :+ Disk () r]
leftCenters,[(Point 2 r :+ p) :+ Disk () r]
rightCenters) = (((Point 2 r :+ p) :+ Disk () r) -> Bool)
-> [(Point 2 r :+ p) :+ Disk () r]
-> ([(Point 2 r :+ p) :+ Disk () r],
    [(Point 2 r :+ p) :+ Disk () r])
forall a. (a -> Bool) -> [a] -> ([a], [a])
List.partition (\(Point 2 r :+ p
r :+ Disk () r
_) -> (Point 2 r :+ p) -> (Point 2 r :+ p) -> (Point 2 r :+ p) -> CCW
forall r a b c.
(Ord r, Num r) =>
(Point 2 r :+ a) -> (Point 2 r :+ b) -> (Point 2 r :+ c) -> CCW
ccw' Point 2 r :+ p
p Point 2 r :+ p
q Point 2 r :+ p
r CCW -> CCW -> Bool
forall a. Eq a => a -> a -> Bool
== CCW
CCW) [(Point 2 r :+ p) :+ Disk () r]
centers
    -- note that we consider 'leftmost' with respect to going from p
    -- to q. This does not really have a global meaning.

    -- we need to find the leftmost and rightmost center on the
    -- bisector. In case there are left-centers, this means that among
    -- the left centers we want to find the point that is furthest way
    -- from p (or q). If there are no left-centers, we with to find
    -- the closest one among the right-centers.
    leftDist :: ((Point 2 r :+ p) :+ Disk () r) -> r
leftDist (Point 2 r :+ p) :+ Disk () r
z = let c :: Point 2 r :+ ()
c = (Point 2 r :+ p) :+ Disk () r
z((Point 2 r :+ p) :+ Disk () r)
-> Getting
     (Point 2 r :+ ()) ((Point 2 r :+ p) :+ Disk () r) (Point 2 r :+ ())
-> Point 2 r :+ ()
forall s a. s -> Getting a s a -> a
^.(Disk () r -> Const (Point 2 r :+ ()) (Disk () r))
-> ((Point 2 r :+ p) :+ Disk () r)
-> Const (Point 2 r :+ ()) ((Point 2 r :+ p) :+ Disk () r)
forall core extra extra'.
Lens (core :+ extra) (core :+ extra') extra extra'
extra((Disk () r -> Const (Point 2 r :+ ()) (Disk () r))
 -> ((Point 2 r :+ p) :+ Disk () r)
 -> Const (Point 2 r :+ ()) ((Point 2 r :+ p) :+ Disk () r))
-> (((Point 2 r :+ ())
     -> Const (Point 2 r :+ ()) (Point 2 r :+ ()))
    -> Disk () r -> Const (Point 2 r :+ ()) (Disk () r))
-> Getting
     (Point 2 r :+ ()) ((Point 2 r :+ p) :+ Disk () r) (Point 2 r :+ ())
forall b c a. (b -> c) -> (a -> b) -> a -> c
.((Point 2 r :+ ()) -> Const (Point 2 r :+ ()) (Point 2 r :+ ()))
-> Disk () r -> Const (Point 2 r :+ ()) (Disk () r)
forall (d :: Nat) p r (d2 :: Nat) p2.
Lens
  (Ball d p r) (Ball d2 p2 r) (Point d r :+ p) (Point d2 r :+ p2)
center
                     s :: r
s = if (Point 2 r :+ p) -> (Point 2 r :+ p) -> (Point 2 r :+ ()) -> CCW
forall r a b c.
(Ord r, Num r) =>
(Point 2 r :+ a) -> (Point 2 r :+ b) -> (Point 2 r :+ c) -> CCW
ccw' Point 2 r :+ p
p Point 2 r :+ p
q Point 2 r :+ ()
c CCW -> CCW -> Bool
forall a. Eq a => a -> a -> Bool
== CCW
CCW then r
1 else -r
1
                 in r
s r -> r -> r
forall a. Num a => a -> a -> a
* Point 2 r -> Point 2 r -> r
forall r (d :: Nat).
(Num r, Arity d) =>
Point d r -> Point d r -> r
squaredEuclideanDist (Point 2 r :+ p
p(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ ()
c(Point 2 r :+ ())
-> Getting (Point 2 r) (Point 2 r :+ ()) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ ()) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)

    dl :: Maybe ((Point 2 r :+ p) :+ Disk () r)
dl = (((Point 2 r :+ p) :+ Disk () r) -> r)
-> [(Point 2 r :+ p) :+ Disk () r]
-> Maybe ((Point 2 r :+ p) :+ Disk () r)
forall b a. Ord b => (a -> b) -> [a] -> Maybe a
maximumOn ((Point 2 r :+ p) :+ Disk () r) -> r
leftDist [(Point 2 r :+ p) :+ Disk () r]
leftCenters  -- disk that has the "leftmost" center
    dr :: Maybe ((Point 2 r :+ p) :+ Disk () r)
dr = (((Point 2 r :+ p) :+ Disk () r) -> r)
-> [(Point 2 r :+ p) :+ Disk () r]
-> Maybe ((Point 2 r :+ p) :+ Disk () r)
forall b a. Ord b => (a -> b) -> [a] -> Maybe a
minimumOn ((Point 2 r :+ p) :+ Disk () r) -> r
leftDist [(Point 2 r :+ p) :+ Disk () r]
rightCenters -- disk that has the "rightmost" center

    -- diameteral disk
    dd :: Disk () r
dd = Point 2 r -> Point 2 r -> Disk () r
forall (d :: Nat) r.
(Arity d, Fractional r) =>
Point d r -> Point d r -> Ball d () r
fromDiameter (Point 2 r :+ p
p(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
q(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)
    mdc :: Maybe (DiskResult p r)
mdc | Disk () r -> [Point 2 r :+ p] -> Bool
forall (t :: * -> *) r p extra.
(Foldable t, Ord r, Num r) =>
Disk p r -> t (Point 2 r :+ extra) -> Bool
isEnclosingDisk Disk () r
dd [Point 2 r :+ p]
ps = DiskResult p r -> Maybe (DiskResult p r)
forall a. a -> Maybe a
Just (DiskResult p r -> Maybe (DiskResult p r))
-> DiskResult p r -> Maybe (DiskResult p r)
forall a b. (a -> b) -> a -> b
$ Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
forall p r.
Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
DiskResult Disk () r
dd ((Point 2 r :+ p) -> (Point 2 r :+ p) -> TwoOrThree (Point 2 r :+ p)
forall a. a -> a -> TwoOrThree a
Two Point 2 r :+ p
p Point 2 r :+ p
q)
        | Bool
otherwise             = Maybe (DiskResult p r)
forall a. Maybe a
Nothing

    -- test if d is an enclosing disk.
    mkEnclosingDisk :: Maybe ((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r)
mkEnclosingDisk  Maybe ((Point 2 r :+ p) :+ Disk () r)
md = Maybe ((Point 2 r :+ p) :+ Disk () r)
md Maybe ((Point 2 r :+ p) :+ Disk () r)
-> (((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r))
-> Maybe (DiskResult p r)
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= ((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r)
mkEnclosingDisk'
    mkEnclosingDisk' :: ((Point 2 r :+ p) :+ Disk () r) -> Maybe (DiskResult p r)
mkEnclosingDisk' (Point 2 r :+ p
r :+ Disk () r
d) | Disk () r -> [Point 2 r :+ p] -> Bool
forall (t :: * -> *) r p extra.
(Foldable t, Ord r, Num r) =>
Disk p r -> t (Point 2 r :+ extra) -> Bool
isEnclosingDisk Disk () r
d [Point 2 r :+ p]
ps = DiskResult p r -> Maybe (DiskResult p r)
forall a. a -> Maybe a
Just (Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
forall p r.
Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
DiskResult Disk () r
d ((Point 2 r :+ p)
-> (Point 2 r :+ p)
-> (Point 2 r :+ p)
-> TwoOrThree (Point 2 r :+ p)
forall a. a -> a -> a -> TwoOrThree a
Three Point 2 r :+ p
p Point 2 r :+ p
q Point 2 r :+ p
r))
                              | Bool
otherwise            = Maybe (DiskResult p r)
forall a. Maybe a
Nothing


isEnclosingDisk   :: (Foldable t, Ord r, Num r)
                  => Disk p r -> t (Point 2 r :+ extra) -> Bool
isEnclosingDisk :: Disk p r -> t (Point 2 r :+ extra) -> Bool
isEnclosingDisk Disk p r
d = ((Point 2 r :+ extra) -> Bool) -> t (Point 2 r :+ extra) -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Point 2 r :+ extra
s -> (Point 2 r :+ extra
s(Point 2 r :+ extra)
-> Getting (Point 2 r) (Point 2 r :+ extra) (Point 2 r)
-> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ extra) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) Point 2 r -> Disk p r -> Bool
forall (d :: Nat) r p.
(Arity d, Ord r, Num r) =>
Point d r -> Ball d p r -> Bool
`inClosedBall` Disk p r
d)

-- | Constructs the initial 'DiskResult' from two points
initial     :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r
initial :: (Point 2 r :+ p) -> (Point 2 r :+ p) -> DiskResult p r
initial Point 2 r :+ p
p Point 2 r :+ p
q = Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
forall p r.
Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
DiskResult (Point 2 r -> Point 2 r -> Disk () r
forall (d :: Nat) r.
(Arity d, Fractional r) =>
Point d r -> Point d r -> Ball d () r
fromDiameter (Point 2 r :+ p
p(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
q(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)) ((Point 2 r :+ p) -> (Point 2 r :+ p) -> TwoOrThree (Point 2 r :+ p)
forall a. a -> a -> TwoOrThree a
Two Point 2 r :+ p
p Point 2 r :+ p
q)

maximumOn   :: Ord b => (a -> b) -> [a] -> Maybe a
maximumOn :: (a -> b) -> [a] -> Maybe a
maximumOn a -> b
f = \case
    [] -> Maybe a
forall a. Maybe a
Nothing
    [a]
xs -> a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ (a -> a -> Ordering) -> [a] -> a
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
List.maximumBy ((a -> b) -> a -> a -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing a -> b
f) [a]
xs

minimumOn   :: Ord b => (a -> b) -> [a] -> Maybe a
minimumOn :: (a -> b) -> [a] -> Maybe a
minimumOn a -> b
f = \case
    [] -> Maybe a
forall a. Maybe a
Nothing
    [a]
xs -> a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ (a -> a -> Ordering) -> [a] -> a
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
List.minimumBy ((a -> b) -> a -> a -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing a -> b
f) [a]
xs


--------------------------------------------------------------------------------

-- test :: Maybe (DiskResult () Rational)
-- test = smallestEnclosingDiskWithPoints p q myPts
--   where
--     p = ext $ Point2 0 (-6)
--     q = ext $ Point2 0 6


-- myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5]

-- disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)
--   where
--     p = ext $ Point2 0 (-6)
--     q = ext $ Point2 0 6


-- maartenBug :: DiskResult () Double
-- maartenBug = let (p:q:rest) = maartenBug'
--              in smallestEnclosingDisk' p q rest

-- maartenBug' :: [Point 2 Double :+ ()]
-- maartenBug' = [ Point2 (7.2784424e-3) (249.23) :+ ()
--               , Point2 (-5.188493   ) (249.23) :+ ()
--               , Point2 (-10.382694  ) (249.23) :+ ()
--               , Point2 (-15.575621  ) (249.23) :+ ()
--               , Point2 (0.0         ) (249.23) :+ ()
--               , Point2 (0.0         ) (239.9031) :+ ()
--               , Point2 (0.0         ) (230.37791) :+ ()
--               , Point2 (0.0         ) (220.67882) :+ ()
--               ]