hierarchical-clustering 0.3.1.2 → 0.4

raw patch · 6 files changed

+629/−256 lines, 6 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• hierarchical-clustering.cabal +28/−9
• src/Data/Clustering/Hierarchical.hs +39/−171
• src/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs +133/−49
• src/Data/Clustering/Hierarchical/Internal/Optimal.hs +237/−0
• src/Data/Clustering/Hierarchical/Internal/Types.hs +84/−0
• tests/runtests.hs +108/−27

@@ -1,6 +1,6 @@ Name:                hierarchical-clustering-Version:             0.3.1.2-Synopsis:            Algorithms for single, average/UPGMA and complete linkage clustering.+Version:             0.4+Synopsis:            Fast algorithms for single, average/UPGMA and complete linkage clustering. License:             BSD3 License-file:        LICENSE Author:              Felipe Almeida Lessa@@ -18,13 +18,30 @@   represents not only the clusters but also the order on which   they were created.   .-  This function uses a naïve algorithm that represents distances-  in a rectangular distance matrix.  There could be space-  improvements (e.g. using a triangular matrix structure) and-  time improvements (e.g. using a finger tree to avoid traversing-  the whole matrix on every iteration just to see what the-  minimum is).+  This package has many implementations with different+  performance characteristics.  There are SLINK and CLINK+  algorithm implementations that are optimal in both space and+  time.  There are also naive implementations using a distance+  matrix.  Using the @dendrogram@ function from+  @Data.Clustering.Hierarchical@ automatically chooses the best+  implementation we have.   .+  Changes in version 0.4:+  .+  * Specialize the distance type to Double for efficiency reasons.+    It's uncommon to use distances other than Double.+  .+  * Implement SLINK and CLINK.  These are optimal algorithms in+    both space and time for single and complete linkage,+    respectively, running in /O(n^2)/ time and /O(n)/ space.+  .+  * Reorganized internal implementation.+  .+  * Some performance improvements for the naive implementation.+  .+  * Better test coverage.  Also, performance improvements for the+    test suite, now running in 3 seconds (instead of one minute).+  .   Changes in version 0.3.1.2 (version 0.3.1.1 was skipped):   .   * Added tests for many things.  Use @cabal test@ =).@@ -65,7 +82,9 @@   Hs-source-dirs: src   Exposed-modules:     Data.Clustering.Hierarchical,-    Data.Clustering.Hierarchical.Internal.DistanceMatrix+    Data.Clustering.Hierarchical.Internal.DistanceMatrix,+    Data.Clustering.Hierarchical.Internal.Optimal,+    Data.Clustering.Hierarchical.Internal.Types   Build-depends:       base       == 4.*     , array      == 0.3.*

@@ -1,45 +1,21 @@ module Data.Clustering.Hierarchical     (-- * Dendrogram data type      Dendrogram(..)+    ,Distance     ,elements     ,cutAt      -- * Linkage data type     ,Linkage(..)-     -- * Generic clustering function+     -- * Clustering function     ,dendrogram-     -- * Functions for specific linkages-    ,singleLinkage-    ,completeLinkage-    ,upgma-    ,fakeAverageLinkage     ) where -import qualified Data.IntMap as IM-import Control.Applicative ((<$>), (<*>))-import Control.Monad.ST (runST)-import Data.Array (listArray, (!))-import Data.Foldable (Foldable (..))-import Data.Function (on)-import Data.Monoid (mappend)-import Data.Traversable (Traversable(..))--import Data.Clustering.Hierarchical.Internal.DistanceMatrix---- | Data structure for storing hierarchical clusters.  The--- distance between clusters is stored on the branches.--- Distances between leafs are the distances between the elements--- on those leafs, while distances between branches are defined--- by the linkage used (see 'Linkage').-data Dendrogram d a =-    Leaf a-    -- ^ The leaf contains the item @a@ itself.-  | Branch d (Dendrogram d a) (Dendrogram d a)-    -- ^ Each branch connects two clusters/dendrograms that are-    -- @d@ distance apart.-    deriving (Eq, Ord, Show)+import Data.Clustering.Hierarchical.Internal.Types (Dendrogram(..), Linkage(..), Distance)+import qualified Data.Clustering.Hierarchical.Internal.DistanceMatrix as DM+import qualified Data.Clustering.Hierarchical.Internal.Optimal as O  -- | List of elements in a dendrogram.-elements :: Dendrogram d a -> [a]+elements :: Dendrogram a -> [a] elements = go []     where       go acc (Leaf x)       = x : acc@@ -69,7 +45,7 @@ -- dendro \`cutAt\` 0.4 == dendro \`cutAt\` 0.2 == [Branch 0.2 (Leaf \'A\') (Leaf \'B\'), Leaf \'C\', Leaf \'D\'] -- dendro \`cutAt\` 0.1 == [Leaf \'A\', Leaf \'B\', Leaf \'C\', Leaf \'D\'] -- no branches at all -- @-cutAt :: Ord d => Dendrogram d a -> d -> [Dendrogram d a]+cutAt :: Dendrogram a -> Distance -> [Dendrogram a] cutAt dendro threshold = go [] dendro     where       go acc x@(Leaf _)                        = x : acc@@ -77,143 +53,35 @@                               | otherwise      = go (go acc r) l  -- cut!  --- | Does not recalculate the distances!-instance Functor (Dendrogram d) where-    fmap f (Leaf d)         = Leaf (f d)-    fmap f (Branch s c1 c2) = Branch s (fmap f c1) (fmap f c2)--instance Foldable (Dendrogram d) where-    foldMap f (Leaf d)         = f d-    foldMap f (Branch _ c1 c2) = foldMap f c1 `mappend` foldMap f c2--instance Traversable (Dendrogram d) where-    traverse f (Leaf d)         = Leaf <$> f d-    traverse f (Branch s c1 c2) = Branch s <$> traverse f c1 <*> traverse f c2----- | The linkage type determines how the distance between--- clusters will be calculated.  These are the linkage types--- currently available on this library.-data Linkage =-    SingleLinkage-  -- ^ The distance between two clusters @a@ and @b@ is the-  -- /minimum/ distance between an element of @a@ and an element-  -- of @b@.-  | CompleteLinkage-  -- ^ The distance between two clusters @a@ and @b@ is the-  -- /maximum/ distance between an element of @a@ and an element-  -- of @b@.-  | UPGMA-  -- ^ Unweighted Pair Group Method with Arithmetic mean, also-  -- called \"average linkage\".  The distance between two-  -- clusters @a@ and @b@ is the /arithmetic average/ between the-  -- distances of all elements in @a@ to all elements in @b@.-  | FakeAverageLinkage-  -- ^ This method is usually wrongly called \"average linkage\".-  -- The distance between cluster @a = a1 U a2@ (that is, cluster-  -- @a@ was formed by the linkage of clusters @a1@ and @a2@) and-  -- an old cluster @b@ is @(d(a1,b) + d(a2,b)) / 2@.  So when-  -- clustering two elements to create a cluster, this method is-  -- the same as UPGMA.  However, in general when joining two-  -- clusters this method assigns equal weights to @a1@ and @a2@,-  -- while UPGMA assigns weights proportional to the number of-  -- elements in each cluster.  See, for example:-  ---  -- *-  -- <http://www.cs.tau.ac.il/~rshamir/algmb/00/scribe00/html/lec08/node21.html>,-  -- which defines the real UPGMA and gives the equation to-  -- calculate the distance between an old and a new cluster.-  ---  -- *-  -- <http://github.com/JadeFerret/ai4r/blob/master/lib/ai4r/clusterers/average_linkage.rb>,-  -- code for \"average linkage\" on ai4r library implementing-  -- what we call here @FakeAverageLinkage@ and not UPGMA.-    deriving (Eq, Ord, Show, Enum)----- Some cluster distances-cdistSingleLinkage      :: Ord d => ClusterDistance d-cdistSingleLinkage      = \(_, d1) (_, d2) -> d1 `min` d2--cdistCompleteLinkage    :: Ord d => ClusterDistance d-cdistCompleteLinkage    = \(_, d1) (_, d2) -> d1 `max` d2--cdistUPGMA              :: Fractional d => ClusterDistance d-cdistUPGMA              = \(b1,d1) (b2,d2) ->-                            let n1 = fromIntegral (size b1)-                                n2 = fromIntegral (size b2)-                            in (n1 * d1 + n2 * d2) / (n1 + n2)--cdistFakeAverageLinkage :: Fractional d => ClusterDistance d-cdistFakeAverageLinkage = \(_, d1) (_, d2) -> (d1 + d2) / 2----- | /O(n^3)/ Calculates a complete, rooted dendrogram for a list--- of items and a linkage type.  If your distance type has an--- 'Ord' instance but not a 'Fractional' one, then please use--- specific functions 'singleLinkage' or 'completeLinkage' that--- have less restrictive types.-dendrogram :: (Ord d, Fractional d)-           => Linkage        -- ^ Linkage type to be used.-           -> [a]            -- ^ Items to be clustered.-           -> (a -> a -> d)  -- ^ Distance function between items.-           -> Dendrogram d a -- ^ Complete dendrogram.-dendrogram linkage = dendrogram' cdist-    where-      cdist = case linkage of-                SingleLinkage      -> cdistSingleLinkage-                CompleteLinkage    -> cdistCompleteLinkage-                FakeAverageLinkage -> cdistFakeAverageLinkage-                UPGMA              -> cdistUPGMA---- | /O(n^3)/ Like 'dendrogram', but specialized to single--- linkage (see 'SingleLinkage') which does not require--- 'Fractional'.-singleLinkage :: Ord d => [a] -> (a -> a -> d) -> Dendrogram d a-singleLinkage = dendrogram' cdistSingleLinkage---- | /O(n^3)/ Like 'dendrogram', but specialized to complete--- linkage (see 'CompleteLinkage') which does not require--- 'Fractional'.-completeLinkage :: Ord d => [a] -> (a -> a -> d) -> Dendrogram d a-completeLinkage = dendrogram' cdistCompleteLinkage---- | /O(n^3)/ Like 'dendrogram', but specialized to 'UPGMA'.-upgma :: (Fractional d, Ord d) => [a] -> (a -> a -> d) -> Dendrogram d a-upgma = dendrogram' cdistUPGMA---- | /O(n^3)/ Like 'dendrogram', but specialized to fake average--- linkage (see 'FakeAverageLinkage').-fakeAverageLinkage :: (Fractional d, Ord d) => [a]-                   -> (a -> a -> d) -> Dendrogram d a-fakeAverageLinkage = dendrogram' cdistFakeAverageLinkage------ | Worker function to create dendrograms based on a--- 'ClusterDistance' (and not a 'Linkage').-dendrogram' :: Ord d => ClusterDistance d-            -> [a] -> (a -> a -> d) -> Dendrogram d a-dendrogram' _ []  _ = error "Data.Clustering.Hierarchical: empty input list"-dendrogram' _ [x] _ = Leaf x-dendrogram' cdist items dist = runST (act ())-    where-      n = length items-      act _noMonomorphismRestrictionPlease = do-        let xs = listArray (1, n) items-        fromDistance (dist `on` (xs !)) n >>= go xs (n-1) IM.empty-      go xs i ds dm = xs `seq` i `seq` ds `seq` dm `seq` do-        ((c1,c2), distance) <- findMin dm-        cu <- mergeClusters cdist dm (c1,c2)-        let dendro c = case size c of-                         1 -> Leaf $! xs ! key c-                         _ -> ds IM.! key c-            d1 = dendro c1-            d2 = dendro c2-            du = d1 `seq` d2 `seq` Branch distance d1 d2-        case i of-          1 -> return du-          _ -> let ds' = IM.insert (key cu) du $-                         IM.delete (key c1) $-                         IM.delete (key c2) ds-               in du `seq` go xs (i-1) ds' dm+-- | Calculates a complete, rooted dendrogram for a list of items+-- and a linkage type.  The following are the time and space+-- complexities for each linkage:+--+-- ['SingleLinkage'] /O(n^2)/ time and /O(n)/ space, using the+--   SLINK algorithm.  This algorithm is optimal in both space+--   and time and gives the same answer as the naive algorithm+--   using a distance matrix.+--+-- ['CompleteLinkage'] /O(n^3)/ time and /O(n^2)/ space, using+--   the naive algorithm with a distance matrix.  Use 'CLINK' if+--   you need more performance.+--+-- [Complete linkage with 'CLINK'] /O(n^2)/ time and /O(n)/+--   space, using the CLINK algorithm.  Note that this algorithm+--   doesn't always give the same answer as the naive algorithm+--   using a distance matrix, but it's much faster.+--+-- ['UPGMA'] /O(n^3)/ time and /O(n^2)/ space, using the naive+--   algorithm with a distance matrix.+--+-- ['FakeAverageLinkage'] /O(n^3)/ time and /O(n^2)/ space, using+-- the naive algorithm with a distance matrix.+dendrogram :: Linkage              -- ^ Linkage type to be used.+           -> [a]                  -- ^ Items to be clustered.+           -> (a -> a -> Distance) -- ^ Distance function between items.+           -> Dendrogram a         -- ^ Complete dendrogram.+dendrogram SingleLinkage      = O.singleLinkage+dendrogram CompleteLinkage    = DM.completeLinkage+dendrogram CLINK              = O.completeLinkage+dendrogram UPGMA              = DM.upgma+dendrogram FakeAverageLinkage = DM.fakeAverageLinkage

@@ -1,30 +1,35 @@+{-# LANGUAGE BangPatterns #-}+ module Data.Clustering.Hierarchical.Internal.DistanceMatrix-    (Cluster(..)-    ,Item-    ,DistMatrix(..)-    ,ClusterDistance-    ,fromDistance-    ,findMin-    ,mergeClusters+    (singleLinkage+    ,completeLinkage+    ,upgma+    ,fakeAverageLinkage     ) where -import qualified Data.IntMap as IM-import Control.Monad (forM_, when)-import Control.Monad.ST (ST)-import Data.Array.ST (STArray, newArray, newListArray, readArray, writeArray)-import Data.List (delete, tails)+-- from base+import Control.Monad (forM_)+import Control.Monad.ST (ST, runST)+import Data.Array (listArray, (!))+import Data.Array.ST (STArray, STUArray, newArray_, newListArray, readArray, writeArray)+import Data.Function (on)+import Data.List (delete, tails, (\\)) import Data.STRef (STRef, newSTRef, readSTRef, writeSTRef) +-- from containers+import qualified Data.IntMap as IM +-- from this package+import Data.Clustering.Hierarchical.Internal.Types+ mkErr :: String -> a mkErr = error . ("Data.Clustering.Hierarchical.Internal.DistanceMatrix." ++)  -- | Internal (to this package) type used to represent a cluster -- (of possibly just one element).  The @key@ should be less than--- or equal to all @more@ elements.-data Cluster = Cluster {key  :: !Item  -- ^ Element used as key.-                       ,more :: [Item] -- ^ Other elements in the cluster.-                       ,size :: !Int   -- ^ At least one, the @key@.+-- or equal to all elements of the cluster.+data Cluster = Cluster { key  :: {-# UNPACK #-} !Item  -- ^ Element used as key.+                       , size :: {-# UNPACK #-} !Int   -- ^ At least one, the @key@.                        }                deriving (Eq, Ord, Show) @@ -33,17 +38,16 @@  -- | Creates a singleton cluster. singleton :: Item -> Cluster-singleton k = Cluster {key = k, more = [], size = 1}+singleton k = Cluster {key = k, size = 1} --- | Joins two clusters, returns the 'key' that didn't become--- 'key' of the new cluster as well.  Clusters are not monoid--- because we don't have 'mempty'.+-- | /O(1)/. Joins two clusters, returns the 'key' that didn't+-- become 'key' of the new cluster as well.  Clusters are not+-- monoid because we don't have 'mempty'. merge :: Cluster -> Cluster -> (Cluster, Item) merge c1 c2 = let (kl,km) = if key c1 < key c2                             then (key c1, key c2)                             else (key c2, key c1)               in (Cluster {key  = kl-                          ,more = km : more c1 ++ more c2                           ,size = size c1 + size c2}                  ,km) @@ -51,24 +55,26 @@   -- | A distance matrix.-data DistMatrix s d = DM {matrix   :: STArray s (Item, Item) d-                         ,active   :: STRef   s [Item]-                         ,clusters :: STArray s Item Cluster}+data DistMatrix s =+    DM { matrix   :: {-# UNPACK #-} !(STUArray s (Item, Item) Distance)+       , active   :: {-# UNPACK #-} !(STRef    s [Item])+       , clusters :: {-# UNPACK #-} !(STArray  s Item Cluster)+       }  --- | /O(n^2)/ Creates a list of possible combinations between the--- given elements.+-- | /O(n^2)/. Creates a list of possible combinations between+-- the given elements. combinations :: [a] -> [(a,a)] combinations xs = [(a,b) | (a:as) <- tails xs, b <- as]  --- | /O(n^2)/ Constructs a new distance matrix from a distance+-- | /O(n^2)/. Constructs a new distance matrix from a distance -- function and a number @n@ of elements.  Elements will be drawn -- from @[1..n]@-fromDistance :: Ord d => (Item -> Item -> d) -> Item -> ST s (DistMatrix s d)+fromDistance :: (Item -> Item -> Distance) -> Item -> ST s (DistMatrix s) fromDistance _ n | n < 2 = mkErr "fromDistance: n < 2 is meaningless" fromDistance dist n = do-  matrix_ <- newArray ((1,2), (n-1,n)) (mkErr "fromDistance: undef element")+  matrix_ <- newArray_ ((1,2), (n-1,n))   active_ <- newSTRef [1..n]   forM_ (combinations [1..n]) $ \x -> writeArray matrix_ x (uncurry dist x)   clusters_ <- newListArray (1,n) (map singleton [1..n])@@ -77,35 +83,59 @@               ,clusters = clusters_}  --- | /O(n^2)/ Returns the minimum distance of the distance+-- | /O(n^2)/. Returns the minimum distance of the distance -- matrix.  The first key given is less than the second key.-findMin :: Ord d => DistMatrix s d -> ST s ((Cluster, Cluster), d)-findMin dm = readSTRef (active dm) >>= go1 . combinations+findMin :: DistMatrix s -> ST s ((Cluster, Cluster), Distance)+findMin dm = readSTRef (active dm) >>= go1     where       matrix_ = matrix dm       choose b i m' = if m' < snd b then (i, m') else b-      go1 (i:is)   = readArray matrix_ i >>= go2 is . (,) i-      go1 []       = mkErr "findMin: empty DistMatrix"-      go2 i b | i `seq` b `seq` False = undefined-      go2 (i:is) b = readArray matrix_ i >>= go2 is . choose b i-      go2 []     b = do c1 <- readArray (clusters dm) (fst $ fst b)-                        c2 <- readArray (clusters dm) (snd $ fst b)-                        return ((c1, c2), snd b) +      go1 is@(i1:i2:_) = do di <- readArray matrix_ (i1, i2) -- initial+                            ((b1, b2), d) <- go2 is ((i1, i2), di)+                            c1 <- readArray (clusters dm) b1+                            c2 <- readArray (clusters dm) b2+                            return ((c1, c2), d)+      go1 _            = mkErr "findMin: empty DistMatrix" +      go2 (i1:is@(_:_)) !b = go3 i1 is b >>= go2 is+      go2 _              b = return b++      go3 i1 (i2:is) !b = readArray matrix_ (i1,i2) >>= go3 i1 is . choose b (i1,i2)+      go3 _  []       b = return b+++ -- | Type for functions that calculate distances between -- clusters.-type ClusterDistance d =-       (Cluster, d)   -- ^ Cluster B1 and distance from A to B1-    -> (Cluster, d)   -- ^ Cluster B2 and distance from A to B2-    -> d              -- ^ Distance from A to (B1 U B2).+type ClusterDistance =+       (Cluster, Distance) -- ^ Cluster B1 and distance from A to B1+    -> (Cluster, Distance) -- ^ Cluster B2 and distance from A to B2+    -> Distance            -- ^ Distance from A to (B1 U B2).  --- | /O(n)/ Merges two clusters, returning the new cluster and+-- Some cluster distances+cdistSingleLinkage      :: ClusterDistance+cdistSingleLinkage      = \(_, d1) (_, d2) -> d1 `min` d2++cdistCompleteLinkage    :: ClusterDistance+cdistCompleteLinkage    = \(_, d1) (_, d2) -> d1 `max` d2++cdistUPGMA              :: ClusterDistance+cdistUPGMA              = \(b1,d1) (b2,d2) ->+                            let n1 = fromIntegral (size b1)+                                n2 = fromIntegral (size b2)+                            in (n1 * d1 + n2 * d2) / (n1 + n2)++cdistFakeAverageLinkage :: ClusterDistance+cdistFakeAverageLinkage = \(_, d1) (_, d2) -> (d1 + d2) / 2++++-- | /O(n)/. Merges two clusters, returning the new cluster and -- the new distance matrix.-mergeClusters :: (Ord d)-              => ClusterDistance d-              -> DistMatrix s d+mergeClusters :: ClusterDistance+              -> DistMatrix s               -> (Cluster, Cluster)               -> ST s Cluster mergeClusters cdist (DM matrix_ active_ clusters_) (b1, b2) = do@@ -118,12 +148,12 @@    -- Calculate new distances   activeV <- readSTRef active_-  forM_ activeV $ \k -> when (k `notElem` [b1k, b2k]) $ do+  forM_ (activeV \\ [b1k, b2k]) $ \k -> do       -- a   <- readArray clusters_ k       d_a_b1 <- readArray matrix_ $ ix k b1k       d_a_b2 <- readArray matrix_ $ ix k b2k       let d = cdist (b1, d_a_b1) (b2, d_a_b2)-      d `seq` writeArray matrix_ (ix k km) d+      writeArray matrix_ (ix k km) $! d    -- Save new cluster, invalidate old one   writeArray clusters_ km bu@@ -132,3 +162,57 @@    -- Return new cluster.   return bu+++-- | Worker function to create dendrograms based on a+-- 'ClusterDistance'.+dendrogram' :: ClusterDistance -> [a] -> (a -> a -> Distance) -> Dendrogram a+dendrogram' _ []  _ = mkErr "dendrogram': empty input list"+dendrogram' _ [x] _ = Leaf x+dendrogram' cdist items dist = runST (act ())+    where+      n = length items+      act _noMonomorphismRestrictionPlease = do+        let xs = listArray (1, n) items+            im = IM.fromDistinctAscList $ zip [1..] $ map Leaf items+        fromDistance (dist `on` (xs !)) n >>= go (n-1) im+      go !i !ds !dm = do+        ((c1,c2), distance) <- findMin dm+        cu <- mergeClusters cdist dm (c1,c2)+        let dendro c = IM.updateLookupWithKey (\_ _ -> Nothing) (key c)+            (Just d1, !ds')  = dendro c1 ds+            (Just d2, !ds'') = dendro c2 ds'+            du = Branch distance d1 d2+        case i of+          1 -> return du+          _ -> let !ds''' = IM.insert (key cu) du ds''+               in du `seq` go (i-1) ds''' dm+++-- | /O(n^3)/ time and /O(n^2)/ space. Calculates a complete,+-- rooted dendrogram for a list of items using single linkage+-- with the naïve algorithm using a distance matrix.+singleLinkage :: [a] -> (a -> a -> Distance) -> Dendrogram a+singleLinkage = dendrogram' cdistSingleLinkage+++-- | /O(n^3)/ time and /O(n^2)/ space. Calculates a complete,+-- rooted dendrogram for a list of items using complete linkage+-- with the naïve algorithm using a distance matrix.+completeLinkage :: [a] -> (a -> a -> Distance) -> Dendrogram a+completeLinkage = dendrogram' cdistCompleteLinkage+++-- | /O(n^3)/ time and /O(n^2)/ space. Calculates a complete,+-- rooted dendrogram for a list of items using UPGMA with the+-- naïve algorithm using a distance matrix.+upgma :: [a] -> (a -> a -> Distance) -> Dendrogram a+upgma = dendrogram' cdistUPGMA+++-- | /O(n^3)/ time and /O(n^2)/ space. Calculates a complete,+-- rooted dendrogram for a list of items using fake average+-- linkage with the naïve algorithm using a distance matrix.+fakeAverageLinkage :: [a]+                   -> (a -> a -> Distance) -> Dendrogram a+fakeAverageLinkage = dendrogram' cdistFakeAverageLinkage

@@ -0,0 +1,237 @@+{-# LANGUAGE BangPatterns #-}++-- | 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 = error . ("Data.Clustering.Hierarchical.Internal.Optimal." ++)+++type Index = Int++data PointerRepresentation s a =+  PR { pi     :: {-# UNPACK #-} !(STUArray s Index Index)+     , lambda :: {-# UNPACK #-} !(STUArray s Index Distance)+     , em     :: {-# UNPACK #-} !(STUArray s Index Distance)+     , 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 n xs' = ($ xs') <$> liftM3 PR (newArray_ (1, n)) (newArray_ (1, n)) (newArray_ (1, n))++indexDistance :: [a] -> (a -> a -> Distance)+              -> (Index, Array Index a, Index -> Index -> Distance)+indexDistance xs dist = (n, xs', dist')+    where+      !n = length xs+      !xs' = listArray (1, n) xs+      dist' i j = dist (xs' ! i) (xs' ! j)+++infinity :: Distance+infinity = 1 / 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 xs dist = initPR n xs' >>= go 1+    where+      (n, xs', dist') = indexDistance xs dist++      go !i !pr | i == n + 1 = return pr+                | otherwise  = do+        writeArray (pi pr)     i i+        writeArray (lambda pr) i infinity+        forM_ [1..i-1] $ \j ->+          writeArray (em pr) j (dist' j i)+        forM_ [1..i-1] $ \j -> do+          lambda_j <- readArray (lambda pr) 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)+        forM_ [1..i-1] $ \j -> do+          pi_j        <- readArray (pi pr)     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+        go (i+1) 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 xs dist = initPR n xs' >>= go 1+    where+      (n, xs', dist') = indexDistance xs dist++      go !i !pr | i == n + 1 = return pr+                | i == 1     = do writeArray (pi pr)     1 1+                                  writeArray (lambda pr) 1 infinity+                                  go 2 pr+                | otherwise  = do+        -- First part+        writeArray (pi pr)     i i+        writeArray (lambda pr) i infinity+        forM_ [1..i-1] $ \j ->+          writeArray (em pr) j (dist' j i)+        forM_ [1..i-1] $ \j -> do+          lambda_j <- readArray (lambda pr) 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 <- readArray (em pr) (i-1) >>= go_a_loop (i-1) pr (i-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] $ \j -> do+          pi_j    <- readArray (pi pr) 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 0 _ a _ = return a+      go_a_loop !j !pr !a !em_a = do+        pi_j     <- readArray (pi pr)     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 !i !pr !a !b !c+          | a >= i - 1 = return ()+          | b <  i - 1 = do pi_b     <- readArray (pi pr)     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+          | otherwise  = do writeArray (pi pr)     b i+                            writeArray (lambda pr) b c+                            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 pr = do+  (1,n) <- getBounds (lambda pr)+  lambdas <- getElems (lambda pr)+  pis     <- getElems (pi pr)+  let sorted = sortBy (\(_,l1,_) (_,l2,_) -> l1 `compare` l2) $+               zip3 [1..] lambdas pis+  index <- newListArray (1,n) [1..]+  let go im [] =+        case IM.toList im of+          [(_,x)] -> return x+          _       -> mkErr "buildDendrogram: final never here"+      go im ((i, (j,lambda_j,pi_j)):rest) = do+        left_i  <- readArray 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 (\_ _ -> Nothing) ix im+                          where ix = negate left_i+            (right, im'') | right_i > 0 = (Leaf $ elm pr ! right_i, im')+                          | otherwise   = first (fromMaybe e2) $+                                          IM.updateLookupWithKey (\_ _ -> Nothing) ix im'+                          where ix = negate right_i+            im''' = IM.insert i (Branch lambda_j left right) im''+            e1 = mkErr "buildDendrogram: never here 1"+            e2 = mkErr "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 []  _   = mkErr "singleLinkage: empty input"+singleLinkage [x] _   = Leaf x+singleLinkage xs dist = runST (slink xs dist >>= 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 []  _   = mkErr "completeLinkage: empty input"+completeLinkage [x] _   = Leaf x+completeLinkage xs dist = runST (clink xs dist >>= buildDendrogram)

@@ -0,0 +1,84 @@+module Data.Clustering.Hierarchical.Internal.Types+    ( Dendrogram(..)+    , Linkage(..)+    , Distance+    ) where++-- from base+import Control.Applicative ((<$>), (<*>))+import Data.Foldable (Foldable (..))+import Data.Monoid (mappend)+import Data.Traversable (Traversable(..))++-- | Data structure for storing hierarchical clusters.  The+-- distance between clusters is stored on the branches.+-- Distances between leafs are the distances between the elements+-- on those leafs, while distances between branches are defined+-- by the linkage used (see 'Linkage').+data Dendrogram a =+    Leaf a+    -- ^ The leaf contains the item @a@ itself.+  | Branch {-# UNPACK #-} !Distance (Dendrogram a) (Dendrogram a)+    -- ^ Each branch connects two clusters/dendrograms that are+    -- @d@ distance apart.+    deriving (Eq, Ord, Show)++-- | A distance is simply a synonym of 'Double' for efficiency.+type Distance = Double++-- | Does not recalculate the distances!+instance Functor Dendrogram where+    fmap f (Leaf d)         = Leaf (f d)+    fmap f (Branch s c1 c2) = Branch s (fmap f c1) (fmap f c2)++instance Foldable Dendrogram where+    foldMap f (Leaf d)         = f d+    foldMap f (Branch _ c1 c2) = foldMap f c1 `mappend` foldMap f c2++instance Traversable Dendrogram where+    traverse f (Leaf d)         = Leaf <$> f d+    traverse f (Branch s c1 c2) = Branch s <$> traverse f c1 <*> traverse f c2+++-- | The linkage type determines how the distance between+-- clusters will be calculated.  These are the linkage types+-- currently available on this library.+data Linkage =+    SingleLinkage+  -- ^ The distance between two clusters @a@ and @b@ is the+  -- /minimum/ distance between an element of @a@ and an element+  -- of @b@.+  | CompleteLinkage+  -- ^ The distance between two clusters @a@ and @b@ is the+  -- /maximum/ distance between an element of @a@ and an element+  -- of @b@.+  | CLINK+  -- ^ The same as 'CompleteLinkage', but using the CLINK+  -- algorithm.  It's much faster however doesn't always give the+  -- best complete linkage dendrogram.+  | UPGMA+  -- ^ Unweighted Pair Group Method with Arithmetic mean, also+  -- called \"average linkage\".  The distance between two+  -- clusters @a@ and @b@ is the /arithmetic average/ between the+  -- distances of all elements in @a@ to all elements in @b@.+  | FakeAverageLinkage+  -- ^ This method is usually wrongly called \"average linkage\".+  -- The distance between cluster @a = a1 U a2@ (that is, cluster+  -- @a@ was formed by the linkage of clusters @a1@ and @a2@) and+  -- an old cluster @b@ is @(d(a1,b) + d(a2,b)) / 2@.  So when+  -- clustering two elements to create a cluster, this method is+  -- the same as UPGMA.  However, in general when joining two+  -- clusters this method assigns equal weights to @a1@ and @a2@,+  -- while UPGMA assigns weights proportional to the number of+  -- elements in each cluster.  See, for example:+  --+  -- *+  -- <http://www.cs.tau.ac.il/~rshamir/algmb/00/scribe00/html/lec08/node21.html>,+  -- which defines the real UPGMA and gives the equation to+  -- calculate the distance between an old and a new cluster.+  --+  -- *+  -- <http://github.com/JadeFerret/ai4r/blob/master/lib/ai4r/clusterers/average_linkage.rb>,+  -- code for \"average linkage\" on ai4r library implementing+  -- what we call here @FakeAverageLinkage@ and not UPGMA.+    deriving (Eq, Ord, Show, Enum)

@@ -1,23 +1,27 @@+{-# LANGUAGE Rank2Types #-}+ -- from base import qualified Control.Exception as E-import Control.Monad (when)-import Data.List (sort)+import Control.Monad (when, liftM2)+import Data.List (delete, sort, nub) import Text.Printf (printf) import Text.Show.Functions ()  -- from hspec-import Test.Hspec.Monadic+import Test.Hspec.Monadic (hspecX, describe, it, pending, Specs) import Test.Hspec.HUnit () import Test.Hspec.QuickCheck (prop)  -- from HUnit-import Test.HUnit+import Test.HUnit ((~?=), Assertion, assertFailure)  -- from QuickCheck-import Test.QuickCheck ((==>))+import Test.QuickCheck (Property, Arbitrary(..), Gen, forAll)  -- from this package import Data.Clustering.Hierarchical+import qualified Data.Clustering.Hierarchical.Internal.DistanceMatrix as DM+import qualified Data.Clustering.Hierarchical.Internal.Optimal as O   main :: IO ()@@ -28,7 +32,7 @@ test_cutAt :: Specs test_cutAt =     describe "cutAt" $ do-      let dendro      :: Dendrogram Double Char+      let dendro      :: Dendrogram Char           dendro      = Branch 0.8 d_0_8_left d_0_8_right           d_0_8_left  =   Branch 0.5 d_0_5_left d_0_5_right           d_0_5_left  =     Branch 0.2 d_0_2_left d_0_2_right@@ -51,40 +55,108 @@  test_dendrogram :: Specs test_dendrogram = do-    describe "dendrogram SingleLinkage" $ do-      basicDendrogramTests SingleLinkage-    describe "dendrogram CompleteLinkage" $ do-      basicDendrogramTests CompleteLinkage-    describe "dendrogram UPGMA" $ do-      basicDendrogramTests UPGMA-    describe "dendrogram FakeAverageLinkage" $ do-      basicDendrogramTests FakeAverageLinkage+    describe "Optimal's singleLinkage" $ do+      basicDendrogramTests O.singleLinkage+      prop "really is single linkage" $+        propCorrectLinkage O.singleLinkage singleLink +    describe "Optimal's completeLinkage" $ do+      basicDendrogramTests O.completeLinkage+      prop "really is complete linkage" $+        propCorrectLinkage O.completeLinkage completeLink -basicDendrogramTests :: Linkage -> Specs-basicDendrogramTests linkage = do-  let f xs = dendrogram linkage xs+    describe "DistanceMatrix's singleLinkage" $ do+      basicDendrogramTests DM.singleLinkage+      prop "really is single linkage" $+        propCorrectLinkage DM.singleLinkage singleLink++    describe "DistanceMatrix's completeLinkage" $ do+      basicDendrogramTests DM.completeLinkage+      prop "really is complete linkage" $+        propCorrectLinkage DM.completeLinkage completeLink++    describe "DistanceMatrix's upgma" $ do+      basicDendrogramTests DM.upgma+      prop "really is UPGMA" $+        propCorrectLinkage DM.upgma upgma++    describe "DistanceMatrix's fakeAverageLinkage" $ do+      basicDendrogramTests DM.fakeAverageLinkage++    describe "Optimal and DistanceMatrix" $ do+      let test f1 f2 = forAll nonNullLists $ \ps ->+                       f1 ps euclideanDist ==== f2 ps euclideanDist+      prop "agree on singleLinkage"   $ test O.singleLinkage DM.singleLinkage+      it "agree on completeLinkage" $+         pending "This doesn't work because CLINK doesn't \+                 \always give the best complete linkage."+++basicDendrogramTests :: (forall a. [a] -> (a -> a -> Distance) -> Dendrogram a) -> Specs+basicDendrogramTests f = do   it "fails for an empty input" $      assertErrors (f [] (\_ _ -> zero))   it "works for one element" $-     Leaf () == f [()] (\_ _ -> zero)+     Leaf () == f [()] undefined   prop "always returns the elements we gave" $-     \xs dist ->-         let dist' x y = abs (dist x y) :: Double-         in not (null (xs :: [Double])) ==>-            elements (f xs dist') `isPermutationOf` xs+     forAll nonNullLists $ \points ->+       elements (f points euclideanDist) `isPermutationOf` points   prop "works for examples where all elements have the same distance" $-     \xs fixedDist ->-         let okay :: Dendrogram Rational Char -> [Char] -> Maybe [Char]-             okay (Leaf z) (y:ys)   | z == y         = Just ys-             okay (Branch d l r) ys | d == fixedDist = okay l ys >>= okay r+     \fixedDist ->+     forAll nonNullLists $ \xs' ->+         let xs = nub xs'++             okay :: Dendrogram Char -> [Char] -> Maybe [Char]+             okay (Leaf z)       ys | z `elem` ys    = Just (delete z ys)+             okay (Branch d l r) ys | d ~= fixedDist = okay l ys >>= okay r              okay _ _ = Nothing-         in not (null xs) ==> okay (f xs (\_ _ -> fixedDist)) xs == Just [] +             dist x y | x == y    = error "shouldn't calculate (dist x x)"+                      | otherwise = fixedDist +         in okay (f xs dist) xs == Just []++----------------------------------------------------------------------++type P = (Double, Double)++propCorrectLinkage :: ([P] -> (P -> P -> Distance) -> Dendrogram P)+                   -> (D P -> [P] -> [P] -> Distance)+                   -> Property+propCorrectLinkage f link =+    forAll nonNullLists $ \xs -> correctLinkage link d (f xs d)+        where d = euclideanDist++type D a = a -> a -> Distance++correctLinkage :: (D a -> [a] -> [a] -> Distance) -> D a -> Dendrogram a -> Bool+correctLinkage link dist = go+    where+      go (Branch d l r) = go l && go r &&+                          link dist (elements l) (elements r) ~= d+      go (Leaf _) = True++singleLink, completeLink, upgma :: D a -> [a] -> [a] -> Distance+singleLink   dist xs ys = minimum [x `dist` y | x <- xs, y <- ys]+completeLink dist xs ys = maximum [x `dist` y | x <- xs, y <- ys]+upgma        dist xs ys = sum [x `dist` y | x <- xs, y <- ys] /+                          fromIntegral (length xs * length ys)++----------------------------------------------------------------------++nonNullLists :: Arbitrary a => Gen [a]+nonNullLists = liftM2 (:) arbitrary arbitrary+ isPermutationOf :: Ord a => [a] -> [a] -> Bool isPermutationOf xs ys = sort xs == sort ys +euclideanDist :: P -> P -> Double+euclideanDist (x1,y1) (x2,y2) = sqrt $ sq (x1-x2) + sq (y1-y2)+    where sq x = x * x++(~=) :: Double -> Double -> Bool+a ~= b = abs (a - b) < 1e-5+ zero :: Double zero = 0 @@ -93,3 +165,12 @@     b <- E.catch (E.evaluate x >> return True)                  (\(E.ErrorCall _) -> return False {- Ok -})     when b $ assertFailure "Didn't raise an 'error'."+++-- | Compare two dendrograms without being concerned about+-- permutations.+(====) :: Eq a => Dendrogram a -> Dendrogram a -> Bool+Leaf x1         ==== Leaf x2         = x1 == x2+Branch d1 l1 r1 ==== Branch d2 l2 r2 = d1 ~= d2 && ((l1 ==== l2 && r1 ==== r2) ||+                                                    (l1 ==== r2 && r1 ==== l2))+_ ==== _ = False