--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.EuclideanMST
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- \(O(n\log n)\) time algorithm algorithm to compute the Euclidean minimum
-- spanning tree of a set of \(n\) points in \(\mathbb{R}^2\).
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.EuclideanMST ( euclideanMST ) where

import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConquer
import           Algorithms.Geometry.DelaunayTriangulation.Types
import           Algorithms.Graph.MST
import           Control.Lens
import           Data.Ext
import           Data.Geometry
import qualified Data.List.NonEmpty as NonEmpty
import           Data.PlaneGraph
import           Data.Proxy
import           Data.Tree

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

-- | Computes the Euclidean Minimum Spanning Tree. We compute the Delaunay
-- Triangulation (DT), and then extract the EMST. Hence, the same restrictions
-- apply as for the DT:
--
-- pre: the input is a *SET*, i.e. contains no duplicate points. (If the input
-- does contain duplicate points, the implementation throws them away)
--
-- running time: \(O(n \log n)\)
euclideanMST     :: (Ord r, Fractional r)
                 => NonEmpty.NonEmpty (Point 2 r :+ p) -> Tree (Point 2 r :+ p)
euclideanMST :: NonEmpty (Point 2 r :+ p) -> Tree (Point 2 r :+ p)
euclideanMST NonEmpty (Point 2 r :+ p)
pts = (\VertexId' MSTW
v -> PlaneGraph MSTW p (r :+ ()) () r
gPlaneGraph MSTW p (r :+ ()) () r
-> Getting
     (Point 2 r) (PlaneGraph MSTW p (r :+ ()) () r) (Point 2 r)
-> Point 2 r
forall s a. s -> Getting a s a -> a
^.VertexId' MSTW
-> Lens' (PlaneGraph MSTW p (r :+ ()) () r) (Point 2 r)
forall k (s :: k) v e f r.
VertexId' s -> Lens' (PlaneGraph s v e f r) (Point 2 r)
locationOf VertexId' MSTW
v Point 2 r -> p -> Point 2 r :+ p
forall core extra. core -> extra -> core :+ extra
:+ PlaneGraph MSTW p (r :+ ()) () r
gPlaneGraph MSTW p (r :+ ()) () r
-> Getting p (PlaneGraph MSTW p (r :+ ()) () r) p -> p
forall s a. s -> Getting a s a -> a
^.VertexId' MSTW
-> Lens'
     (PlaneGraph MSTW p (r :+ ()) () r)
     (DataOf (PlaneGraph MSTW p (r :+ ()) () r) (VertexId' MSTW))
forall g i. HasDataOf g i => i -> Lens' g (DataOf g i)
dataOf VertexId' MSTW
v) (VertexId' MSTW -> Point 2 r :+ p)
-> Tree (VertexId' MSTW) -> Tree (Point 2 r :+ p)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Tree (VertexId' MSTW)
t
  where
    -- since we care only about the relative order of the edges we can use the
    -- squared Euclidean distance rather than the Euclidean distance, thus
    -- avoiding the Floating constraint
    g :: PlaneGraph MSTW p (r :+ ()) () r
g = (Point 2 r -> Point 2 r -> r)
-> PlaneGraph MSTW p () () r -> PlaneGraph MSTW p (r :+ ()) () r
forall k r a (s :: k) p e f.
(Point 2 r -> Point 2 r -> a)
-> PlaneGraph s p e f r -> PlaneGraph s p (a :+ e) f r
withEdgeDistances Point 2 r -> Point 2 r -> r
forall r (d :: Nat).
(Num r, Arity d) =>
Point d r -> Point d r -> r
squaredEuclideanDist (PlaneGraph MSTW p () () r -> PlaneGraph MSTW p (r :+ ()) () r)
-> (NonEmpty (Point 2 r :+ p) -> PlaneGraph MSTW p () () r)
-> NonEmpty (Point 2 r :+ p)
-> PlaneGraph MSTW p (r :+ ()) () r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Proxy MSTW -> Triangulation p r -> PlaneGraph MSTW p () () r
forall k (proxy :: k -> *) (s :: k) p r.
proxy s -> Triangulation p r -> PlaneGraph s p () () r
toPlaneGraph (Proxy MSTW
forall k (t :: k). Proxy t
Proxy :: Proxy MSTW)
      (Triangulation p r -> PlaneGraph MSTW p () () r)
-> (NonEmpty (Point 2 r :+ p) -> Triangulation p r)
-> NonEmpty (Point 2 r :+ p)
-> PlaneGraph MSTW p () () r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. NonEmpty (Point 2 r :+ p) -> Triangulation p r
forall r p.
(Ord r, Fractional r) =>
NonEmpty (Point 2 r :+ p) -> Triangulation p r
delaunayTriangulation (NonEmpty (Point 2 r :+ p) -> PlaneGraph MSTW p (r :+ ()) () r)
-> NonEmpty (Point 2 r :+ p) -> PlaneGraph MSTW p (r :+ ()) () r
forall a b. (a -> b) -> a -> b
$ NonEmpty (Point 2 r :+ p)
pts
    t :: Tree (VertexId' MSTW)
t = PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ()
-> Tree (VertexId' MSTW)
forall k e (s :: k) (w :: World) v f.
Ord e =>
PlanarGraph s w v e f -> Tree (VertexId s w)
mst (PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ()
 -> Tree (VertexId' MSTW))
-> PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ()
-> Tree (VertexId' MSTW)
forall a b. (a -> b) -> a -> b
$ PlaneGraph MSTW p (r :+ ()) () r
gPlaneGraph MSTW p (r :+ ()) () r
-> Getting
     (PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ())
     (PlaneGraph MSTW p (r :+ ()) () r)
     (PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ())
-> PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ()
forall s a. s -> Getting a s a -> a
^.Getting
  (PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ())
  (PlaneGraph MSTW p (r :+ ()) () r)
  (PlanarGraph MSTW 'Primal (VertexData r p) (r :+ ()) ())
forall k1 (s1 :: k1) v1 e1 f1 r1 k2 (s2 :: k2) v2 e2 f2 r2.
Iso
  (PlaneGraph s1 v1 e1 f1 r1)
  (PlaneGraph s2 v2 e2 f2 r2)
  (PlanarGraph s1 'Primal (VertexData r1 v1) e1 f1)
  (PlanarGraph s2 'Primal (VertexData r2 v2) e2 f2)
graph


data MSTW