{-# LANGUAGE ScopedTypeVariables #-}
module Algorithms.Geometry.ClosestPair.DivideAndConquer( closestPair
, CP
, CCP(..)
, mergePairs
)
where
import Algorithms.DivideAndConquer
import Control.Lens
import Data.Ext
import Data.Geometry.Point
import Data.LSeq (LSeq)
import qualified Data.LSeq as LSeq
import qualified Data.List as List
import Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NonEmpty
import Data.Ord (comparing)
import Data.Semigroup.Foldable (toNonEmpty)
import Data.UnBounded
import Data.Util
closestPair :: (Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p)
closestPair :: LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p)
closestPair = CCP p r -> Two (Point 2 r :+ p)
forall p r. CCP p r -> Two (Point 2 r :+ p)
f (CCP p r -> Two (Point 2 r :+ p))
-> (LSeq 2 (Point 2 r :+ p) -> CCP p r)
-> LSeq 2 (Point 2 r :+ p)
-> Two (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> CCP p r)
-> NonEmpty (Point 2 r :+ p) -> CCP p r
forall (f :: * -> *) s a.
(Foldable1 f, Semigroup s) =>
(a -> s) -> f a -> s
divideAndConquer1 (Point 2 r :+ p) -> CCP p r
forall r p. (Point 2 r :+ p) -> CCP p r
mkCCP (NonEmpty (Point 2 r :+ p) -> CCP p r)
-> (LSeq 2 (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> LSeq 2 (Point 2 r :+ p)
-> CCP p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LSeq 2 (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall (t :: * -> *) a. Foldable1 t => t a -> NonEmpty a
toNonEmpty
(LSeq 2 (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (LSeq 2 (Point 2 r :+ p) -> LSeq 2 (Point 2 r :+ p))
-> LSeq 2 (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering)
-> LSeq 2 (Point 2 r :+ p) -> LSeq 2 (Point 2 r :+ p)
forall a (n :: Nat). (a -> a -> Ordering) -> LSeq n a -> LSeq n a
LSeq.unstableSortBy (((Point 2 r :+ p) -> Point 2 r)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((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))
where
mkCCP :: (Point 2 r :+ p) -> CCP p r
mkCCP Point 2 r :+ p
p = NonEmpty (Point 2 r :+ p) -> CP (Point 2 r :+ p) r -> CCP p r
forall p r.
NonEmpty (Point 2 r :+ p) -> CP (Point 2 r :+ p) r -> CCP p r
CCP (Point 2 r :+ p
p (Point 2 r :+ p) -> [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p)
forall a. a -> [a] -> NonEmpty a
:| []) CP (Point 2 r :+ p) r
forall a. Top a
Top
f :: CCP p r -> Two (Point 2 r :+ p)
f = \case
CCP NonEmpty (Point 2 r :+ p)
_ (ValT (SP Two (Point 2 r :+ p)
cp r
_)) -> Two (Point 2 r :+ p)
cp
CCP NonEmpty (Point 2 r :+ p)
_ Top (SP (Two (Point 2 r :+ p)) r)
Top -> [Char] -> Two (Point 2 r :+ p)
forall a. HasCallStack => [Char] -> a
error [Char]
"closestPair: absurd."
type CP q r = Top (SP (Two q) r)
data CCP p r = CCP (NonEmpty (Point 2 r :+ p))
!(CP (Point 2 r :+ p) r)
deriving (Int -> CCP p r -> ShowS
[CCP p r] -> ShowS
CCP p r -> [Char]
(Int -> CCP p r -> ShowS)
-> (CCP p r -> [Char]) -> ([CCP p r] -> ShowS) -> Show (CCP p r)
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
forall p r. (Show r, Show p) => Int -> CCP p r -> ShowS
forall p r. (Show r, Show p) => [CCP p r] -> ShowS
forall p r. (Show r, Show p) => CCP p r -> [Char]
showList :: [CCP p r] -> ShowS
$cshowList :: forall p r. (Show r, Show p) => [CCP p r] -> ShowS
show :: CCP p r -> [Char]
$cshow :: forall p r. (Show r, Show p) => CCP p r -> [Char]
showsPrec :: Int -> CCP p r -> ShowS
$cshowsPrec :: forall p r. (Show r, Show p) => Int -> CCP p r -> ShowS
Show,CCP p r -> CCP p r -> Bool
(CCP p r -> CCP p r -> Bool)
-> (CCP p r -> CCP p r -> Bool) -> Eq (CCP p r)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall p r. (Eq r, Eq p) => CCP p r -> CCP p r -> Bool
/= :: CCP p r -> CCP p r -> Bool
$c/= :: forall p r. (Eq r, Eq p) => CCP p r -> CCP p r -> Bool
== :: CCP p r -> CCP p r -> Bool
$c== :: forall p r. (Eq r, Eq p) => CCP p r -> CCP p r -> Bool
Eq)
instance (Num r, Ord r) => Semigroup (CCP p r) where
(CCP NonEmpty (Point 2 r :+ p)
ptsl CP (Point 2 r :+ p) r
cpl) <> :: CCP p r -> CCP p r -> CCP p r
<> (CCP NonEmpty (Point 2 r :+ p)
ptsr CP (Point 2 r :+ p) r
cpr) = NonEmpty (Point 2 r :+ p) -> CP (Point 2 r :+ p) r -> CCP p r
forall p r.
NonEmpty (Point 2 r :+ p) -> CP (Point 2 r :+ p) r -> CCP p r
CCP (((Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall a.
(a -> a -> Ordering) -> NonEmpty a -> NonEmpty a -> NonEmpty a
mergeSortedBy (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
cmp NonEmpty (Point 2 r :+ p)
ptsl NonEmpty (Point 2 r :+ p)
ptsr)
(CP (Point 2 r :+ p) r
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> CP (Point 2 r :+ p) r
forall p r.
(Ord r, Num r) =>
CP (Point 2 r :+ p) r
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> CP (Point 2 r :+ p) r
mergePairs ((CP (Point 2 r :+ p) r -> Top r)
-> CP (Point 2 r :+ p) r
-> CP (Point 2 r :+ p) r
-> CP (Point 2 r :+ p) r
forall b a. Ord b => (a -> b) -> a -> a -> a
minBy CP (Point 2 r :+ p) r -> Top r
forall a r. CP a r -> Top r
getDist CP (Point 2 r :+ p) r
cpl CP (Point 2 r :+ p) r
cpr) NonEmpty (Point 2 r :+ p)
ptsl NonEmpty (Point 2 r :+ p)
ptsr)
where
cmp :: Point 2 r :+ p -> Point 2 r :+ p -> Ordering
cmp :: (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
cmp Point 2 r :+ p
p Point 2 r :+ p
q = ((Point 2 r :+ p) -> r)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Point 2 r :+ p) -> Getting r (Point 2 r :+ p) r -> r
forall s a. s -> Getting a s a -> a
^.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> Getting r (Point 2 r :+ p) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord) Point 2 r :+ p
p Point 2 r :+ p
q Ordering -> Ordering -> Ordering
forall a. Semigroup a => a -> a -> a
<> ((Point 2 r :+ p) -> r)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Point 2 r :+ p) -> Getting r (Point 2 r :+ p) r -> r
forall s a. s -> Getting a s a -> a
^.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> Getting r (Point 2 r :+ p) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(1 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
xCoord) Point 2 r :+ p
p Point 2 r :+ p
q
mergePairs :: forall p r. (Ord r, Num r)
=> CP (Point 2 r :+ p) r
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> CP (Point 2 r :+ p) r
mergePairs :: CP (Point 2 r :+ p) r
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> CP (Point 2 r :+ p) r
mergePairs CP (Point 2 r :+ p) r
cp' NonEmpty (Point 2 r :+ p)
ls' NonEmpty (Point 2 r :+ p)
rs' = CP (Point 2 r :+ p) r
-> [Point 2 r :+ p] -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
go CP (Point 2 r :+ p) r
cp' (NonEmpty (Point 2 r :+ p) -> [Point 2 r :+ p]
forall a. NonEmpty a -> [a]
NonEmpty.toList NonEmpty (Point 2 r :+ p)
ls') (NonEmpty (Point 2 r :+ p) -> [Point 2 r :+ p]
forall a. NonEmpty a -> [a]
NonEmpty.toList NonEmpty (Point 2 r :+ p)
rs')
where
go :: CP (Point 2 r :+ p) r -> [Point 2 r :+ p] -> [Point 2 r :+ p]
-> CP (Point 2 r :+ p) r
go :: CP (Point 2 r :+ p) r
-> [Point 2 r :+ p] -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
go CP (Point 2 r :+ p) r
cp [Point 2 r :+ p]
_ [] = CP (Point 2 r :+ p) r
cp
go CP (Point 2 r :+ p) r
cp [Point 2 r :+ p]
ls (Point 2 r :+ p
r:[Point 2 r :+ p]
rs) = let ls'' :: [Point 2 r :+ p]
ls'' = Top r -> [Point 2 r :+ p] -> (Point 2 r :+ p) -> [Point 2 r :+ p]
forall r q a.
(Ord r, Num r) =>
Top r -> [Point 2 r :+ q] -> (Point 2 r :+ a) -> [Point 2 r :+ q]
trim (CP (Point 2 r :+ p) r -> Top r
forall a r. CP a r -> Top r
getDist CP (Point 2 r :+ p) r
cp) [Point 2 r :+ p]
ls Point 2 r :+ p
r
cp'' :: CP (Point 2 r :+ p) r
cp'' = CP (Point 2 r :+ p) r
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
forall r p.
(Ord r, Num r) =>
CP (Point 2 r :+ p) r
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
run CP (Point 2 r :+ p) r
cp Point 2 r :+ p
r [Point 2 r :+ p]
ls''
in CP (Point 2 r :+ p) r
-> [Point 2 r :+ p] -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
go CP (Point 2 r :+ p) r
cp'' [Point 2 r :+ p]
ls'' [Point 2 r :+ p]
rs
trim :: (Ord r, Num r)
=> Top r -> [Point 2 r :+ q] -> Point 2 r :+ a
-> [Point 2 r :+ q]
trim :: Top r -> [Point 2 r :+ q] -> (Point 2 r :+ a) -> [Point 2 r :+ q]
trim (ValT r
d) [Point 2 r :+ q]
ls Point 2 r :+ a
r = ((Point 2 r :+ q) -> Bool) -> [Point 2 r :+ q] -> [Point 2 r :+ q]
forall a. (a -> Bool) -> [a] -> [a]
List.dropWhile (\Point 2 r :+ q
l -> (Point 2 r :+ q) -> (Point 2 r :+ a) -> r
forall r p q.
(Ord r, Num r) =>
(Point 2 r :+ p) -> (Point 2 r :+ q) -> r
sqVertDist Point 2 r :+ q
l Point 2 r :+ a
r r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
d) [Point 2 r :+ q]
ls
trim Top r
_ [Point 2 r :+ q]
ls Point 2 r :+ a
_ = [Point 2 r :+ q]
ls
sqVertDist :: (Ord r, Num r) => Point 2 r :+ p -> Point 2 r :+ q -> r
sqVertDist :: (Point 2 r :+ p) -> (Point 2 r :+ q) -> r
sqVertDist Point 2 r :+ p
l Point 2 r :+ q
r = let d :: r
d = r
0 r -> r -> r
forall a. Ord a => a -> a -> a
`max` (Point 2 r :+ q
r(Point 2 r :+ q) -> Getting r (Point 2 r :+ q) r -> r
forall s a. s -> Getting a s a -> a
^.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ q) -> Const r (Point 2 r :+ q)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ q) -> Const r (Point 2 r :+ q))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> Getting r (Point 2 r :+ q) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord r -> r -> r
forall a. Num a => a -> a -> a
- Point 2 r :+ p
l(Point 2 r :+ p) -> Getting r (Point 2 r :+ p) r -> r
forall s a. s -> Getting a s a -> a
^.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ p) -> Const r (Point 2 r :+ p))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> Getting r (Point 2 r :+ p) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord) in r
dr -> r -> r
forall a. Num a => a -> a -> a
*r
d
run :: (Ord r, Num r)
=> CP (Point 2 r :+ p) r -> Point 2 r :+ p -> [Point 2 r :+ p]
-> CP (Point 2 r :+ p) r
run :: CP (Point 2 r :+ p) r
-> (Point 2 r :+ p) -> [Point 2 r :+ p] -> CP (Point 2 r :+ p) r
run CP (Point 2 r :+ p) r
cp'' Point 2 r :+ p
r [Point 2 r :+ p]
ls =
CP (Point 2 r :+ p) r
-> [Point 2 r :+ p]
-> (CP (Point 2 r :+ p) r -> (Point 2 r :+ p) -> Bool)
-> (CP (Point 2 r :+ p) r
-> (Point 2 r :+ p) -> CP (Point 2 r :+ p) r)
-> CP (Point 2 r :+ p) r
forall s a. s -> [a] -> (s -> a -> Bool) -> (s -> a -> s) -> s
runWhile CP (Point 2 r :+ p) r
cp'' [Point 2 r :+ p]
ls
(\CP (Point 2 r :+ p) r
cp Point 2 r :+ p
l -> r -> Top r
forall a. a -> Top a
ValT ((Point 2 r :+ p) -> (Point 2 r :+ p) -> r
forall r p q.
(Ord r, Num r) =>
(Point 2 r :+ p) -> (Point 2 r :+ q) -> r
sqVertDist Point 2 r :+ p
r Point 2 r :+ p
l) Top r -> Top r -> Bool
forall a. Ord a => a -> a -> Bool
< CP (Point 2 r :+ p) r -> Top r
forall a r. CP a r -> Top r
getDist CP (Point 2 r :+ p) r
cp)
(\CP (Point 2 r :+ p) r
cp Point 2 r :+ p
l -> (CP (Point 2 r :+ p) r -> Top r)
-> CP (Point 2 r :+ p) r
-> CP (Point 2 r :+ p) r
-> CP (Point 2 r :+ p) r
forall b a. Ord b => (a -> b) -> a -> a -> a
minBy CP (Point 2 r :+ p) r -> Top r
forall a r. CP a r -> Top r
getDist CP (Point 2 r :+ p) r
cp (SP (Two (Point 2 r :+ p)) r -> CP (Point 2 r :+ p) r
forall a. a -> Top a
ValT (SP (Two (Point 2 r :+ p)) r -> CP (Point 2 r :+ p) r)
-> SP (Two (Point 2 r :+ p)) r -> CP (Point 2 r :+ p) r
forall a b. (a -> b) -> a -> b
$ Two (Point 2 r :+ p) -> r -> SP (Two (Point 2 r :+ p)) r
forall a b. a -> b -> SP a b
SP ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Two (Point 2 r :+ p)
forall a. a -> a -> Two a
Two Point 2 r :+ p
l Point 2 r :+ p
r) ((Point 2 r :+ p) -> (Point 2 r :+ p) -> r
forall (d :: Nat) r extra extra.
(ImplicitPeano (Peano d), Num r,
ArityPeano (Peano (FromPeano (Peano d))),
KnownNat (FromPeano (Peano d)), KnownNat d,
Peano (FromPeano (Peano d) + 1)
~ 'S (Peano (FromPeano (Peano d)))) =>
(Point d r :+ extra) -> (Point d r :+ extra) -> r
dist Point 2 r :+ p
l Point 2 r :+ p
r)))
where
dist :: (Point d r :+ extra) -> (Point d r :+ extra) -> r
dist (Point d r
p :+ extra
_) (Point d r
q :+ extra
_) = Point d r -> Point d r -> r
forall r (d :: Nat).
(Num r, Arity d) =>
Point d r -> Point d r -> r
squaredEuclideanDist Point d r
p Point d r
q
runWhile :: s -> [a] -> (s -> a -> Bool) -> (s -> a -> s) -> s
runWhile :: s -> [a] -> (s -> a -> Bool) -> (s -> a -> s) -> s
runWhile s
s' [a]
ys s -> a -> Bool
p s -> a -> s
f = s -> [a] -> s
go s
s' [a]
ys
where
go :: s -> [a] -> s
go s
s [] = s
s
go s
s (a
x:[a]
xs) | s -> a -> Bool
p s
s a
x = s -> [a] -> s
go (s -> a -> s
f s
s a
x) [a]
xs
| Bool
otherwise = s
s
minBy :: Ord b => (a -> b) -> a -> a -> a
minBy :: (a -> b) -> a -> a -> a
minBy a -> b
f a
a a
b | a -> b
f a
a b -> b -> Bool
forall a. Ord a => a -> a -> Bool
< a -> b
f a
b = a
a
| Bool
otherwise = a
b
getDist :: CP a r -> Top r
getDist :: CP a r -> Top r
getDist = (SP (Two a) r -> r) -> CP a r -> Top r
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Getting r (SP (Two a) r) r -> SP (Two a) r -> r
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting r (SP (Two a) r) r
forall s t a b. Field2 s t a b => Lens s t a b
_2)