packages feed

hierarchical-clustering 0.3.1 → 0.3.1.2

raw patch · 6 files changed

+475/−354 lines, 6 filesdep +HUnitdep +QuickCheckdep +hierarchical-clusteringPVP ok

version bump matches the API change (PVP)

Dependencies added: HUnit, QuickCheck, hierarchical-clustering, hspec

API changes (from Hackage documentation)

Files

− Data/Clustering/Hierarchical.hs 
@@ -1,217 +0,0 @@-module Data.Clustering.Hierarchical-    (-- * Dendrogram data type-     Dendrogram(..)-    ,elements-    ,cutAt-     -- * Linkage data type-    ,Linkage(..)-     -- * Generic 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)---- | List of elements in a dendrogram.-elements :: Dendrogram d a -> [a]-elements = go []-    where-      go acc (Leaf x)       = x : acc-      go acc (Branch _ l r) = go (go acc r) l---- | @dendro \`cutAt\` threshold@ cuts the dendrogram @dendro@ at--- all branches which have distances strictly greater than--- @threshold@.------ For example, suppose we have------ @--- dendro = Branch 0.8---            (Branch 0.5---              (Branch 0.2---                (Leaf \'A\')---                (Leaf \'B\'))---              (Leaf \'C\'))---            (Leaf \'D\')--- @------ Then:------ @--- dendro \`cutAt\` 0.9 == dendro \`cutAt\` 0.8 == [dendro] -- no changes--- dendro \`cutAt\` 0.7 == dendro \`cutAt\` 0.5 == [Branch 0.5 (Branch 0.2 (Leaf \'A\') (Leaf \'B\')) (Leaf \'C\'), Leaf \'D\']--- 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 dendro threshold = go [] dendro-    where-      go acc x@(Leaf _)                        = x : acc-      go acc x@(Branch d l r) | d <= threshold = x : acc-                              | 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' 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
− Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs 
@@ -1,134 +0,0 @@-module Data.Clustering.Hierarchical.Internal.DistanceMatrix-    (Cluster(..)-    ,Item-    ,DistMatrix(..)-    ,ClusterDistance-    ,fromDistance-    ,findMin-    ,mergeClusters-    ) 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)-import Data.STRef (STRef, newSTRef, readSTRef, writeSTRef)---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@.-                       }-               deriving (Eq, Ord, Show)---- | An element of a cluster.-type Item = IM.Key---- | Creates a singleton cluster.-singleton :: Item -> Cluster-singleton k = Cluster {key = k, more = [], 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'.-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)------- | A distance matrix.-data DistMatrix s d = DM {matrix   :: STArray s (Item, Item) d-                         ,active   :: STRef   s [Item]-                         ,clusters :: STArray s Item Cluster}----- | /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--- 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 _ n | n < 2 = mkErr "fromDistance: n < 2 is meaningless"-fromDistance dist n = do-  matrix_ <- newArray ((1,2), (n-1,n)) (mkErr "fromDistance: undef element")-  active_ <- newSTRef [1..n]-  forM_ (combinations [1..n]) $ \x -> writeArray matrix_ x (uncurry dist x)-  clusters_ <- newListArray (1,n) (map singleton [1..n])-  return $ DM {matrix   = matrix_-              ,active   = active_-              ,clusters = clusters_}----- | /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-    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)----- | 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).----- | /O(n)/ Merges two clusters, returning the new cluster and--- the new distance matrix.-mergeClusters :: (Ord d)-              => ClusterDistance d-              -> DistMatrix s d-              -> (Cluster, Cluster)-              -> ST s Cluster-mergeClusters cdist (DM matrix_ active_ clusters_) (b1, b2) = do-  let (bu, kl) = b1 `merge` b2-      b1k = key b1-      b2k = key b2-      km  = key bu-      ix i j | i < j     = (i,j)-             | otherwise = (j,i)--  -- Calculate new distances-  activeV <- readSTRef active_-  forM_ activeV $ \k -> when (k `notElem` [b1k, b2k]) $ 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--  -- Save new cluster, invalidate old one-  writeArray clusters_ km bu-  writeArray clusters_ kl $ mkErr "mergeClusters: invalidated"-  writeSTRef active_ $ delete kl activeV--  -- Return new cluster.-  return bu
hierarchical-clustering.cabal view
@@ -1,5 +1,5 @@ Name:                hierarchical-clustering-Version:             0.3.1+Version:             0.3.1.2 Synopsis:            Algorithms for single, average/UPGMA and complete linkage clustering. License:             BSD3 License-file:        LICENSE@@ -7,7 +7,7 @@ Maintainer:          felipe.lessa@gmail.com Category:            Clustering Build-type:          Simple-Cabal-version:       >= 1.6+Cabal-version:       >= 1.8 Description:   This package provides a function to create a dendrogram from a   list of items and a distance function between them.  Initially@@ -25,6 +25,10 @@   the whole matrix on every iteration just to see what the   minimum is).   .+  Changes in version 0.3.1.2 (version 0.3.1.1 was skipped):+  .+  * Added tests for many things.  Use @cabal test@ =).+  .   Changes in version 0.3.1:   .   * Works with containers 0.4 (thanks, Doug Beardsley).@@ -49,6 +53,8 @@     useful if you want to create a dendrogram and your distance     data type isn't an instance of @Floating@. +Extra-source-files:+  tests/runtests.hs  Source-repository head   type: darcs@@ -56,8 +62,26 @@   Library+  Hs-source-dirs: src   Exposed-modules:     Data.Clustering.Hierarchical,     Data.Clustering.Hierarchical.Internal.DistanceMatrix-  Build-depends: base == 4.*, array == 0.3.*, containers >= 0.3 && < 0.5+  Build-depends:+      base       == 4.*+    , array      == 0.3.*+    , containers >= 0.3 && < 0.5+  GHC-options: -Wall++Test-suite runtests+  Type: exitcode-stdio-1.0+  Hs-source-dirs: tests+  Main-is: runtests.hs+  Build-depends:+      base       == 4.*++    , hspec      == 0.9.*+    , HUnit      == 1.2.*+    , QuickCheck == 2.4.*++    , hierarchical-clustering   GHC-options: -Wall
+ src/Data/Clustering/Hierarchical.hs view
@@ -0,0 +1,219 @@+module Data.Clustering.Hierarchical+    (-- * Dendrogram data type+     Dendrogram(..)+    ,elements+    ,cutAt+     -- * Linkage data type+    ,Linkage(..)+     -- * Generic 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)++-- | List of elements in a dendrogram.+elements :: Dendrogram d a -> [a]+elements = go []+    where+      go acc (Leaf x)       = x : acc+      go acc (Branch _ l r) = go (go acc r) l++-- | @dendro \`cutAt\` threshold@ cuts the dendrogram @dendro@ at+-- all branches which have distances strictly greater than+-- @threshold@.+--+-- For example, suppose we have+--+-- @+-- dendro = Branch 0.8+--            (Branch 0.5+--              (Branch 0.2+--                (Leaf \'A\')+--                (Leaf \'B\'))+--              (Leaf \'C\'))+--            (Leaf \'D\')+-- @+--+-- Then:+--+-- @+-- dendro \`cutAt\` 0.9 == dendro \`cutAt\` 0.8 == [dendro] -- no changes+-- dendro \`cutAt\` 0.7 == dendro \`cutAt\` 0.5 == [Branch 0.5 (Branch 0.2 (Leaf \'A\') (Leaf \'B\')) (Leaf \'C\'), Leaf \'D\']+-- 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 dendro threshold = go [] dendro+    where+      go acc x@(Leaf _)                        = x : acc+      go acc x@(Branch d l r) | d <= threshold = x : acc+                              | 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
+ src/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs view
@@ -0,0 +1,134 @@+module Data.Clustering.Hierarchical.Internal.DistanceMatrix+    (Cluster(..)+    ,Item+    ,DistMatrix(..)+    ,ClusterDistance+    ,fromDistance+    ,findMin+    ,mergeClusters+    ) 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)+import Data.STRef (STRef, newSTRef, readSTRef, writeSTRef)+++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@.+                       }+               deriving (Eq, Ord, Show)++-- | An element of a cluster.+type Item = IM.Key++-- | Creates a singleton cluster.+singleton :: Item -> Cluster+singleton k = Cluster {key = k, more = [], 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'.+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)+++++-- | A distance matrix.+data DistMatrix s d = DM {matrix   :: STArray s (Item, Item) d+                         ,active   :: STRef   s [Item]+                         ,clusters :: STArray s Item Cluster}+++-- | /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+-- 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 _ n | n < 2 = mkErr "fromDistance: n < 2 is meaningless"+fromDistance dist n = do+  matrix_ <- newArray ((1,2), (n-1,n)) (mkErr "fromDistance: undef element")+  active_ <- newSTRef [1..n]+  forM_ (combinations [1..n]) $ \x -> writeArray matrix_ x (uncurry dist x)+  clusters_ <- newListArray (1,n) (map singleton [1..n])+  return $ DM {matrix   = matrix_+              ,active   = active_+              ,clusters = clusters_}+++-- | /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+    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)+++-- | 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).+++-- | /O(n)/ Merges two clusters, returning the new cluster and+-- the new distance matrix.+mergeClusters :: (Ord d)+              => ClusterDistance d+              -> DistMatrix s d+              -> (Cluster, Cluster)+              -> ST s Cluster+mergeClusters cdist (DM matrix_ active_ clusters_) (b1, b2) = do+  let (bu, kl) = b1 `merge` b2+      b1k = key b1+      b2k = key b2+      km  = key bu+      ix i j | i < j     = (i,j)+             | otherwise = (j,i)++  -- Calculate new distances+  activeV <- readSTRef active_+  forM_ activeV $ \k -> when (k `notElem` [b1k, b2k]) $ 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++  -- Save new cluster, invalidate old one+  writeArray clusters_ km bu+  writeArray clusters_ kl $ mkErr "mergeClusters: invalidated"+  writeSTRef active_ $ delete kl activeV++  -- Return new cluster.+  return bu
+ tests/runtests.hs view
@@ -0,0 +1,95 @@+-- from base+import qualified Control.Exception as E+import Control.Monad (when)+import Data.List (sort)+import Text.Printf (printf)+import Text.Show.Functions ()++-- from hspec+import Test.Hspec.Monadic+import Test.Hspec.HUnit ()+import Test.Hspec.QuickCheck (prop)++-- from HUnit+import Test.HUnit++-- from QuickCheck+import Test.QuickCheck ((==>))++-- from this package+import Data.Clustering.Hierarchical+++main :: IO ()+main = hspecX $ do+         test_cutAt+         test_dendrogram++test_cutAt :: Specs+test_cutAt =+    describe "cutAt" $ do+      let dendro      :: Dendrogram Double 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+          d_0_2_left  =       Leaf 'A'+          d_0_2_right =       Leaf 'B'+          d_0_5_right =     Leaf 'C'+          d_0_8_right =   Leaf 'D'++      let testFor threshold expected =+              it (printf "works for 'dendro' with threshold %0.1f" threshold) $+                 dendro `cutAt` threshold ~?= expected++      testFor 0.9 [dendro]+      testFor 0.8 [dendro]+      testFor 0.7 [d_0_8_left, d_0_8_right]+      testFor 0.5 [d_0_8_left, d_0_8_right]+      testFor 0.4 [d_0_5_left, d_0_5_right, d_0_8_right]+      testFor 0.2 [d_0_5_left, d_0_5_right, d_0_8_right]+      testFor 0.1 [d_0_2_left, d_0_2_right, d_0_5_right, d_0_8_right]++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+++basicDendrogramTests :: Linkage -> Specs+basicDendrogramTests linkage = do+  let f xs = dendrogram linkage xs+  it "fails for an empty input" $+     assertErrors (f [] (\_ _ -> zero))+  it "works for one element" $+     Leaf () == f [()] (\_ _ -> zero)+  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+  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+             okay _ _ = Nothing+         in not (null xs) ==> okay (f xs (\_ _ -> fixedDist)) xs == Just []+++isPermutationOf :: Ord a => [a] -> [a] -> Bool+isPermutationOf xs ys = sort xs == sort ys++zero :: Double+zero = 0++assertErrors :: a -> Assertion+assertErrors x = do+    b <- E.catch (E.evaluate x >> return True)+                 (\(E.ErrorCall _) -> return False {- Ok -})+    when b $ assertFailure "Didn't raise an 'error'."