hierarchical-clustering (empty) → 0.1
raw patch · 5 files changed
+317/−0 lines, 5 filesdep +arraydep +basedep +containerssetup-changed
Dependencies added: array, base, containers
Files
- Data/Clustering/Hierarchical.hs +117/−0
- Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs +135/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- hierarchical-clustering.cabal +33/−0
+ Data/Clustering/Hierarchical.hs view
@@ -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:+ --+ -- *+ -- <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)+++-- | 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
+ Data/Clustering/Hierarchical/Internal/DistanceMatrix.hs view
@@ -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
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ hierarchical-clustering.cabal view
@@ -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