diff --git a/Data/Clustering/Hierarchical.hs b/Data/Clustering/Hierarchical.hs deleted file mode 100644 --- a/Data/Clustering/Hierarchical.hs +++ /dev/null @@ -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: - -- - -- * - -- , - -- which defines the real UPGMA and gives the equation to - -- calculate the distance between an old and a new cluster. - -- - -- * - -- , - -- 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 diff --git a/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs b/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs deleted file mode 100644 --- a/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs +++ /dev/null @@ -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 diff --git a/hierarchical-clustering.cabal b/hierarchical-clustering.cabal --- a/hierarchical-clustering.cabal +++ b/hierarchical-clustering.cabal @@ -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 diff --git a/src/Data/Clustering/Hierarchical.hs b/src/Data/Clustering/Hierarchical.hs new file mode 100644 --- /dev/null +++ b/src/Data/Clustering/Hierarchical.hs @@ -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: + -- + -- * + -- , + -- which defines the real UPGMA and gives the equation to + -- calculate the distance between an old and a new cluster. + -- + -- * + -- , + -- 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 diff --git a/src/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs b/src/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs new file mode 100644 --- /dev/null +++ b/src/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs @@ -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 diff --git a/tests/runtests.hs b/tests/runtests.hs new file mode 100644 --- /dev/null +++ b/tests/runtests.hs @@ -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'."