{-# LANGUAGE BangPatterns, FlexibleContexts #-}

-- | Implementations that are optimal in space and time.
module Data.Clustering.Hierarchical.Internal.Optimal
    ( singleLinkage
    , completeLinkage
    ) where

-- from base
import Prelude hiding (pi)
import Control.Applicative ((<$>))
import Control.Arrow (first)
import Control.Monad (forM_, liftM3, when)
import Control.Monad.ST (ST, runST)
import Data.Array (Array, listArray, (!))
import Data.Array.ST (STUArray, newArray_, newListArray,
                      readArray, writeArray,
                      getElems, getBounds) -- getAssocs
import Data.List (sortBy)
import Data.Maybe (fromMaybe)

-- from containers
import qualified Data.IntMap as IM

-- from this package
import Data.Clustering.Hierarchical.Internal.Types


mkErr :: String -> a
mkErr :: forall a. String -> a
mkErr = String -> a
forall a. HasCallStack => String -> a
error (String -> a) -> (String -> String) -> String -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (String
"Data.Clustering.Hierarchical.Internal.Optimal." String -> String -> String
forall a. [a] -> [a] -> [a]
++)


type Index = Int

data PointerRepresentation s a =
  PR { forall s a. PointerRepresentation s a -> STUArray s Index Index
pi     :: {-# UNPACK #-} !(STUArray s Index Index)
     , forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda :: {-# UNPACK #-} !(STUArray s Index Distance)
     , forall s a. PointerRepresentation s a -> STUArray s Index Distance
em     :: {-# UNPACK #-} !(STUArray s Index Distance)
     , forall s a. PointerRepresentation s a -> Array Index a
elm    :: {-# UNPACK #-} !(Array Index a)
     }

-- debugPR :: Show a => PointerRepresentation s a -> ST s String
-- debugPR pr = do
--   pis     <- getAssocs (pi pr)
--   lambdas <- getAssocs (lambda pr)
--   ems     <- getAssocs (em pr)
--   return $ unlines [ "pi     = " ++ show pis
--                    , "lambda = " ++ show lambdas
--                    , "em     = " ++ show ems
--                    , "elm    = " ++ show (elm pr)
--                    ]

initPR :: Index -> Array Index a -> ST s (PointerRepresentation s a)
initPR :: forall a s.
Index -> Array Index a -> ST s (PointerRepresentation s a)
initPR Index
n Array Index a
xs' = ((Array Index a -> PointerRepresentation s a)
-> Array Index a -> PointerRepresentation s a
forall a b. (a -> b) -> a -> b
$ Array Index a
xs') ((Array Index a -> PointerRepresentation s a)
 -> PointerRepresentation s a)
-> ST s (Array Index a -> PointerRepresentation s a)
-> ST s (PointerRepresentation s a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (STUArray s Index Index
 -> STUArray s Index Distance
 -> STUArray s Index Distance
 -> Array Index a
 -> PointerRepresentation s a)
-> ST s (STUArray s Index Index)
-> ST s (STUArray s Index Distance)
-> ST s (STUArray s Index Distance)
-> ST s (Array Index a -> PointerRepresentation s a)
forall (m :: * -> *) a1 a2 a3 r.
Monad m =>
(a1 -> a2 -> a3 -> r) -> m a1 -> m a2 -> m a3 -> m r
liftM3 STUArray s Index Index
-> STUArray s Index Distance
-> STUArray s Index Distance
-> Array Index a
-> PointerRepresentation s a
forall s a.
STUArray s Index Index
-> STUArray s Index Distance
-> STUArray s Index Distance
-> Array Index a
-> PointerRepresentation s a
PR ((Index, Index) -> ST s (STUArray s Index Index)
forall i. Ix i => (i, i) -> ST s (STUArray s i Index)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ (Index
1, Index
n)) ((Index, Index) -> ST s (STUArray s Index Distance)
forall i. Ix i => (i, i) -> ST s (STUArray s i Distance)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ (Index
1, Index
n)) ((Index, Index) -> ST s (STUArray s Index Distance)
forall i. Ix i => (i, i) -> ST s (STUArray s i Distance)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> m (a i e)
newArray_ (Index
1, Index
n))

indexDistance :: [a] -> (a -> a -> Distance)
              -> (Index, Array Index a, Index -> Index -> Distance)
indexDistance :: forall a.
[a]
-> (a -> a -> Distance)
-> (Index, Array Index a, Index -> Index -> Distance)
indexDistance [a]
xs a -> a -> Distance
dist = (Index
n, Array Index a
xs', Index -> Index -> Distance
dist')
    where
      !n :: Index
n = [a] -> Index
forall a. [a] -> Index
forall (t :: * -> *) a. Foldable t => t a -> Index
length [a]
xs
      !xs' :: Array Index a
xs' = (Index, Index) -> [a] -> Array Index a
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Index
1, Index
n) [a]
xs
      dist' :: Index -> Index -> Distance
dist' Index
i Index
j = a -> a -> Distance
dist (Array Index a
xs' Array Index a -> Index -> a
forall i e. Ix i => Array i e -> i -> e
! Index
i) (Array Index a
xs' Array Index a -> Index -> a
forall i e. Ix i => Array i e -> i -> e
! Index
j)


infinity :: Distance
infinity :: Distance
infinity = Distance
1 Distance -> Distance -> Distance
forall a. Fractional a => a -> a -> a
/ Distance
0


-- | /O(n^2)/ time and /O(n)/ space.  See 'singleLinkage' on this module.
slink :: [a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
slink :: forall a s.
[a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
slink [a]
xs a -> a -> Distance
dist = Index -> Array Index a -> ST s (PointerRepresentation s a)
forall a s.
Index -> Array Index a -> ST s (PointerRepresentation s a)
initPR Index
n Array Index a
xs' ST s (PointerRepresentation s a)
-> (PointerRepresentation s a -> ST s (PointerRepresentation s a))
-> ST s (PointerRepresentation s a)
forall a b. ST s a -> (a -> ST s b) -> ST s b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Index
-> PointerRepresentation s a -> ST s (PointerRepresentation s a)
forall {m :: * -> *} {s} {a}.
(MArray (STUArray s) Index m, MArray (STUArray s) Distance m) =>
Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go Index
1
    where
      (Index
n, Array Index a
xs', Index -> Index -> Distance
dist') = [a]
-> (a -> a -> Distance)
-> (Index, Array Index a, Index -> Index -> Distance)
forall a.
[a]
-> (a -> a -> Distance)
-> (Index, Array Index a, Index -> Index -> Distance)
indexDistance [a]
xs a -> a -> Distance
dist

      go :: Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go !Index
i !PointerRepresentation s a
pr | Index
i Index -> Index -> Bool
forall a. Eq a => a -> a -> Bool
== Index
n Index -> Index -> Index
forall a. Num a => a -> a -> a
+ Index
1 = PointerRepresentation s a -> m (PointerRepresentation s a)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return PointerRepresentation s a
pr
                | Bool
otherwise  = do
        STUArray s Index Index -> Index -> Index -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
i Index
i
        STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
i Distance
infinity
        [Index] -> (Index -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index
1..Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1] ((Index -> m ()) -> m ()) -> (Index -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Index
j ->
          STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
em PointerRepresentation s a
pr) Index
j (Index -> Index -> Distance
dist' Index
j Index
i)
        [Index] -> (Index -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index
1..Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1] ((Index -> m ()) -> m ()) -> (Index -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Index
j -> do
          lambda_j <- STUArray s Index Distance -> Index -> m Distance
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
j
          em_j     <- readArray (em pr)     j
          pi_j     <- readArray (pi pr)     j
          em_pi_j  <- readArray (em pr)     pi_j
          if lambda_j >= em_j then do
            writeArray (em pr)     pi_j (em_pi_j `min` lambda_j)
            writeArray (lambda pr) j    em_j
            writeArray (pi pr)     j    i
           else
            writeArray (em pr)     pi_j (em_pi_j `min` em_j)
        [Index] -> (Index -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index
1..Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1] ((Index -> m ()) -> m ()) -> (Index -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Index
j -> do
          pi_j        <- STUArray s Index Index -> Index -> m Index
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
j
          lambda_j    <- readArray (lambda pr) j
          lambda_pi_j <- readArray (lambda pr) pi_j
          when (lambda_j >= lambda_pi_j) $
            writeArray (pi pr) j i
        Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go (Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
+Index
1) PointerRepresentation s a
pr


-- | /O(n^2)/ time and /O(n)/ space. See 'completeLinkage' on this module.
clink :: [a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
clink :: forall a s.
[a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
clink [a]
xs a -> a -> Distance
dist = Index -> Array Index a -> ST s (PointerRepresentation s a)
forall a s.
Index -> Array Index a -> ST s (PointerRepresentation s a)
initPR Index
n Array Index a
xs' ST s (PointerRepresentation s a)
-> (PointerRepresentation s a -> ST s (PointerRepresentation s a))
-> ST s (PointerRepresentation s a)
forall a b. ST s a -> (a -> ST s b) -> ST s b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Index
-> PointerRepresentation s a -> ST s (PointerRepresentation s a)
forall {m :: * -> *} {s} {a}.
(MArray (STUArray s) Index m, MArray (STUArray s) Distance m) =>
Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go Index
1
    where
      (Index
n, Array Index a
xs', Index -> Index -> Distance
dist') = [a]
-> (a -> a -> Distance)
-> (Index, Array Index a, Index -> Index -> Distance)
forall a.
[a]
-> (a -> a -> Distance)
-> (Index, Array Index a, Index -> Index -> Distance)
indexDistance [a]
xs a -> a -> Distance
dist

      go :: Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go !Index
i !PointerRepresentation s a
pr | Index
i Index -> Index -> Bool
forall a. Eq a => a -> a -> Bool
== Index
n Index -> Index -> Index
forall a. Num a => a -> a -> a
+ Index
1 = PointerRepresentation s a -> m (PointerRepresentation s a)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return PointerRepresentation s a
pr
                | Index
i Index -> Index -> Bool
forall a. Eq a => a -> a -> Bool
== Index
1     = do STUArray s Index Index -> Index -> Index -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
1 Index
1
                                  STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
1 Distance
infinity
                                  Index -> PointerRepresentation s a -> m (PointerRepresentation s a)
go Index
2 PointerRepresentation s a
pr
                | Bool
otherwise  = do
        -- First part
        STUArray s Index Index -> Index -> Index -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
i Index
i
        STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
i Distance
infinity
        [Index] -> (Index -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index
1..Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1] ((Index -> m ()) -> m ()) -> (Index -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Index
j ->
          STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
em PointerRepresentation s a
pr) Index
j (Index -> Index -> Distance
dist' Index
j Index
i)
        [Index] -> (Index -> m ()) -> m ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index
1..Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1] ((Index -> m ()) -> m ()) -> (Index -> m ()) -> m ()
forall a b. (a -> b) -> a -> b
$ \Index
j -> do
          lambda_j <- STUArray s Index Distance -> Index -> m Distance
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
j
          em_j     <- readArray (em pr)     j
          when (lambda_j < em_j) $ do
            pi_j     <- readArray (pi pr)     j
            em_pi_j  <- readArray (em pr)     pi_j
            writeArray (em pr) pi_j (em_pi_j `max` em_j)
            writeArray (em pr) j    infinity

        -- Loop a
        a <- STUArray s Index Distance -> Index -> m Distance
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
em PointerRepresentation s a
pr) (Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1) m Distance -> (Distance -> m Index) -> m Index
forall a b. m a -> (a -> m b) -> m b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Index -> PointerRepresentation s a -> Index -> Distance -> m Index
forall {m :: * -> *} {s} {a}.
(MArray (STUArray s) Index m, MArray (STUArray s) Distance m) =>
Index -> PointerRepresentation s a -> Index -> Distance -> m Index
go_a_loop (Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1) PointerRepresentation s a
pr (Index
iIndex -> Index -> Index
forall a. Num a => a -> a -> a
-Index
1)

        -- Loop b
        b <- readArray (pi pr)     a
        c <- readArray (lambda pr) a
        writeArray (pi pr)     a i
        writeArray (lambda pr) a =<< readArray (em pr) a
        go_b_loop i pr a b c

        -- Final part
        forM_ [1..i-1] $ \Index
j -> do
          pi_j    <- STUArray s Index Index -> Index -> m Index
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr) Index
j
          pi_pi_j <- readArray (pi pr) pi_j
          when (pi_pi_j == i) $ do
            lambda_j    <- readArray (lambda pr) j
            lambda_pi_j <- readArray (lambda pr) pi_j
            when (lambda_j >= lambda_pi_j) $
              writeArray (pi pr) j i

        -- Recurse
        go (i+1) pr

      -- Loop a's core
      go_a_loop :: Index -> PointerRepresentation s a -> Index -> Distance -> m Index
go_a_loop Index
0 PointerRepresentation s a
_ Index
a Distance
_ = Index -> m Index
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return Index
a
      go_a_loop !Index
j !PointerRepresentation s a
pr !Index
a !Distance
em_a = do
        pi_j     <- STUArray s Index Index -> Index -> m Index
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
j
        lambda_j <- readArray (lambda pr) j
        em_pi_j  <- readArray (em pr)     pi_j
        if lambda_j >= em_pi_j then do
          em_j <- readArray (em pr) j
          if em_j < em_a then
            go_a_loop (j-1) pr j em_j
           else
            go_a_loop (j-1) pr a em_a
         else do
          writeArray (em pr) j infinity
          go_a_loop (j-1) pr a em_a

      -- Loop b's core
      go_b_loop :: Index
-> PointerRepresentation s a -> Index -> Index -> Distance -> m ()
go_b_loop !Index
i !PointerRepresentation s a
pr !Index
a !Index
b !Distance
c
          | Index
a Index -> Index -> Bool
forall a. Ord a => a -> a -> Bool
>= Index
i Index -> Index -> Index
forall a. Num a => a -> a -> a
- Index
1 = () -> m ()
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
          | Index
b Index -> Index -> Bool
forall a. Ord a => a -> a -> Bool
<  Index
i Index -> Index -> Index
forall a. Num a => a -> a -> a
- Index
1 = do pi_b     <- STUArray s Index Index -> Index -> m Index
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
b
                            lambda_b <- readArray (lambda pr) b
                            writeArray (pi pr)     b i
                            writeArray (lambda pr) b c
                            go_b_loop i pr a pi_b lambda_b
          | Bool
otherwise  = do STUArray s Index Index -> Index -> Index -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Index
forall s a. PointerRepresentation s a -> STUArray s Index Index
pi PointerRepresentation s a
pr)     Index
b Index
i
                            STUArray s Index Distance -> Index -> Distance -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr) Index
b Distance
c
                            () -> m ()
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return ()


-- | /O(n log n)/ time and /O(n)/ space. Construct a 'Dendrogram'
-- from a 'PointerRepresentation'.
buildDendrogram :: PointerRepresentation s a
                -> ST s (Dendrogram a)
buildDendrogram :: forall s a. PointerRepresentation s a -> ST s (Dendrogram a)
buildDendrogram PointerRepresentation s a
pr = do
  bounds <- STUArray s Index Distance -> ST s (Index, Index)
forall i. Ix i => STUArray s i Distance -> ST s (i, i)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> m (i, i)
getBounds (PointerRepresentation s a -> STUArray s Index Distance
forall s a. PointerRepresentation s a -> STUArray s Index Distance
lambda PointerRepresentation s a
pr)
  let (1,n) = bounds
  lambdas <- getElems (lambda pr)
  pis     <- getElems (pi pr)
  let sorted = ((Index, Distance, Index) -> (Index, Distance, Index) -> Ordering)
-> [(Index, Distance, Index)] -> [(Index, Distance, Index)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (\(Index
_,Distance
l1,Index
_) (Index
_,Distance
l2,Index
_) -> Distance
l1 Distance -> Distance -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Distance
l2) ([(Index, Distance, Index)] -> [(Index, Distance, Index)])
-> [(Index, Distance, Index)] -> [(Index, Distance, Index)]
forall a b. (a -> b) -> a -> b
$
               [Index] -> [Distance] -> [Index] -> [(Index, Distance, Index)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [Index
1..] [Distance]
lambdas [Index]
pis
  index <- newListArray (1,n) [1..]
  let go IntMap (Dendrogram a)
im [] =
        case IntMap (Dendrogram a) -> [(Index, Dendrogram a)]
forall a. IntMap a -> [(Index, a)]
IM.toList IntMap (Dendrogram a)
im of
          [(Index
_,Dendrogram a
x)] -> Dendrogram a -> m (Dendrogram a)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return Dendrogram a
x
          [(Index, Dendrogram a)]
_       -> String -> m (Dendrogram a)
forall a. String -> a
mkErr String
"buildDendrogram: final never here"
      go IntMap (Dendrogram a)
im ((Index
i, (Index
j,Distance
lambda_j,Index
pi_j)):[(Index, (Index, Distance, Index))]
rest) = do
        left_i  <- STUArray s Index Index -> Index -> m Index
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray STUArray s Index Index
index Index
j
        right_i <- readArray index pi_j
        writeArray (index `asTypeOf` pi pr) pi_j (negate i)
        let (left,  im')  | left_i > 0  = (Leaf $ elm pr ! left_i, im)
                          | otherwise   = first (fromMaybe e1) $
                                          IM.updateLookupWithKey (\Index
_ Dendrogram a
_ -> Maybe (Dendrogram a)
forall a. Maybe a
Nothing) ix im
                          where ix = Index -> Index
forall a. Num a => a -> a
negate Index
left_i
            (right, im'') | right_i > 0 = (Leaf $ elm pr ! right_i, im')
                          | otherwise   = first (fromMaybe e2) $
                                          IM.updateLookupWithKey (\Index
_ Dendrogram a
_ -> Maybe (Dendrogram a)
forall a. Maybe a
Nothing) ix im'
                          where ix = Index -> Index
forall a. Num a => a -> a
negate Index
right_i
            im''' = Index
-> Dendrogram a -> IntMap (Dendrogram a) -> IntMap (Dendrogram a)
forall a. Index -> a -> IntMap a -> IntMap a
IM.insert Index
i (Distance -> Dendrogram a -> Dendrogram a -> Dendrogram a
forall a. Distance -> Dendrogram a -> Dendrogram a -> Dendrogram a
Branch Distance
lambda_j Dendrogram a
left Dendrogram a
right) IntMap (Dendrogram a)
im''
            e1 = String -> a
forall a. String -> a
mkErr String
"buildDendrogram: never here 1"
            e2 = String -> a
forall a. String -> a
mkErr String
"buildDendrogram: never here 2"
        go im''' rest
  go IM.empty (zip [1..n-1] sorted)


-- | /O(n^2)/ time and /O(n)/ space. Calculates a complete,
-- rooted dendrogram for a list of items using single linkage
-- with the SLINK algorithm.  This algorithm is optimal in space
-- and time.
--
-- [Reference] R. Sibson (1973). \"SLINK: an optimally efficient
--   algorithm for the single-link cluster method\". /The/
--   /Computer Journal/ (British Computer Society) 16 (1):
--   30-34.
singleLinkage :: [a] -> (a -> a -> Distance) -> Dendrogram a
singleLinkage :: forall a. [a] -> (a -> a -> Distance) -> Dendrogram a
singleLinkage []  a -> a -> Distance
_   = String -> Dendrogram a
forall a. String -> a
mkErr String
"singleLinkage: empty input"
singleLinkage [a
x] a -> a -> Distance
_   = a -> Dendrogram a
forall a. a -> Dendrogram a
Leaf a
x
singleLinkage [a]
xs a -> a -> Distance
dist = (forall s. ST s (Dendrogram a)) -> Dendrogram a
forall a. (forall s. ST s a) -> a
runST ([a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
forall a s.
[a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
slink [a]
xs a -> a -> Distance
dist ST s (PointerRepresentation s a)
-> (PointerRepresentation s a -> ST s (Dendrogram a))
-> ST s (Dendrogram a)
forall a b. ST s a -> (a -> ST s b) -> ST s b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= PointerRepresentation s a -> ST s (Dendrogram a)
forall s a. PointerRepresentation s a -> ST s (Dendrogram a)
buildDendrogram)


-- | /O(n^2)/ time and /O(n)/ space. Calculates a complete, rooted dendrogram for a list
-- of items using complete linkage with the CLINK algorithm.  This
-- algorithm is optimal in space and time.
--
-- [Reference] D. Defays (1977). \"An efficient algorithm for a
--   complete link method\". /The Computer Journal/ (British
--   Computer Society) 20 (4): 364-366.
completeLinkage :: [a] -> (a -> a -> Distance) -> Dendrogram a
completeLinkage :: forall a. [a] -> (a -> a -> Distance) -> Dendrogram a
completeLinkage []  a -> a -> Distance
_   = String -> Dendrogram a
forall a. String -> a
mkErr String
"completeLinkage: empty input"
completeLinkage [a
x] a -> a -> Distance
_   = a -> Dendrogram a
forall a. a -> Dendrogram a
Leaf a
x
completeLinkage [a]
xs a -> a -> Distance
dist = (forall s. ST s (Dendrogram a)) -> Dendrogram a
forall a. (forall s. ST s a) -> a
runST ([a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
forall a s.
[a] -> (a -> a -> Distance) -> ST s (PointerRepresentation s a)
clink [a]
xs a -> a -> Distance
dist ST s (PointerRepresentation s a)
-> (PointerRepresentation s a -> ST s (Dendrogram a))
-> ST s (Dendrogram a)
forall a b. ST s a -> (a -> ST s b) -> ST s b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= PointerRepresentation s a -> ST s (Dendrogram a)
forall s a. PointerRepresentation s a -> ST s (Dendrogram a)
buildDendrogram)