diff --git a/Data/Clustering/Hierarchical.hs b/Data/Clustering/Hierarchical.hs new file mode 100644 --- /dev/null +++ b/Data/Clustering/Hierarchical.hs @@ -0,0 +1,117 @@ +module Data.Clustering.Hierarchical + (Dendrogram(..) + ,Linkage(..) + ,completeDendrogram + ) 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. +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) + +-- | 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. +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) + + +-- | Calculates distances between clusters according to the +-- chosen linkage. +clusterDistance :: (Fractional d, Ord d) => Linkage -> ClusterDistance d +clusterDistance SingleLinkage = \_ (_, d1) (_, d2) _ -> d1 `min` d2 +clusterDistance CompleteLinkage = \_ (_, d1) (_, d2) _ -> d1 `max` d2 +clusterDistance FakeAverageLinkage = \_ (_, d1) (_, d2) _ -> (d1 + d2) / 2 +clusterDistance UPGMA = \_ (b1,d1) (b2,d2) _ -> + let n1 = fromIntegral (size b1) + n2 = fromIntegral (size b2) + in (n1 * d1 + n2 * d2) / (n1 + n2) + + +-- | /O(n^2)/ Calculates a complete, rooted dendrogram for a list +-- of items and a distance function. +completeDendrogram :: (Fractional d, Ord d) => Linkage -> + [a] -> (a -> a -> d) -> Dendrogram d a +completeDendrogram linkage items dist = runST (act ()) + where + n = length items + cdist = clusterDistance linkage + act _ = do + let xs = listArray (1, n) items + fromDistance (dist `on` (xs !)) n >>= go xs (n-1) IM.empty + go xs i ds dm = 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 = 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 go xs (i-1) ds' dm diff --git a/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs b/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs new file mode 100644 --- /dev/null +++ b/Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs @@ -0,0 +1,135 @@ +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: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 -- ^ Cluster A + -> (Cluster, d) -- ^ Cluster B1 and distance from A to B1 + -> (Cluster, d) -- ^ Cluster B2 and distance from A to B2 + -> Cluster -- ^ Cluster B = B1 U B2 + -> d -- ^ Distance from A to B. + + +-- | /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 a (b1, d_a_b1) (b2, d_a_b2) bu + 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/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c)2010, Felipe Almeida Lessa + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Felipe Almeida Lessa nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/hierarchical-clustering.cabal b/hierarchical-clustering.cabal new file mode 100644 --- /dev/null +++ b/hierarchical-clustering.cabal @@ -0,0 +1,33 @@ +Name: hierarchical-clustering +Version: 0.1 +Synopsis: Algorithms for single, average/UPGMA and complete linkage clustering. +License: BSD3 +License-file: LICENSE +Author: Felipe Almeida Lessa +Maintainer: felipe.lessa@gmail.com +Category: Clustering +Build-type: Simple +Cabal-version: >= 1.6 +Description: + This package provides a function to create a dendrogram from a + list of items and a distance function between them. Initially + a singleton cluster is created for each item, and then new, + bigger clusters are created by merging the two clusters with + least distance between them. The distance between two clusters + is calculated according to the linkage type. The dendrogram + 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). + +Library + Exposed-modules: + Data.Clustering.Hierarchical, + Data.Clustering.Hierarchical.Internal.DistanceMatrix + Build-depends: base == 4.*, array == 0.3.*, containers == 0.3.* + GHC-options: -Wall