packages feed

clustering (empty) → 0.1.0

raw patch · 9 files changed

+441/−0 lines, 9 filesdep +basedep +clusteringdep +containerssetup-changed

Dependencies added: base, clustering, containers, criterion, hierarchical-clustering, mwc-random, tasty, tasty-hunit, vector

Files

+ LICENSE view
@@ -0,0 +1,20 @@+Copyright (c) 2015 Kai Zhang++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmarks/bench.hs view
@@ -0,0 +1,38 @@+import Control.Monad (replicateM)+import Criterion.Main+import qualified Data.Clustering.Hierarchical as C+import qualified Data.Vector as V+import System.Random.MWC++import AI.Clustering.Hierarchical+import AI.Clustering.Hierarchical.Types ((!))++randSample :: IO [V.Vector Double]+randSample = do+    g <- create+    replicateM 2000 $ uniformVector g 5++main :: IO ()+main = do+    xs <- randSample+    let dists = computeDists euclidean $ V.fromList xs+        fn i j = dists ! (i,j)+    defaultMain+      [ bgroup "AI.Clustering.Hierarchical"+          [ bench "Average Linkage (n = 10)" $+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 10+          , bench "Average Linkage (n = 100)" $+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 100+          , bench "Average Linkage (n = 1000)" $+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 1000+          ]++      , bgroup "Data.Clustering.Hierarchical"+          [ bench "Average Linkage (n = 10)" $+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 10 xs+          , bench "Average Linkage (n = 100)" $+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 100 xs+          , bench "Average Linkage (n = 1000)" $+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 1000 xs+          ]+      ]
+ clustering.cabal view
@@ -0,0 +1,69 @@+-- Initial fastcluster.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                clustering+version:             0.1.0+synopsis:            fast clustering algorithms+description:         O(N^2) implementations for a wide range of hierarchical+                     clustering schemes, including complete linkage, single linkage,+                     average linkage, weighted linkage, and Ward's method.+license:             MIT+license-file:        LICENSE+author:              Kai Zhang+maintainer:          kai@kzhang.org+copyright:           (c) 2015 Kai Zhang+category:            Math+build-type:          Simple+-- extra-source-files:  +cabal-version:       >=1.10++library+  exposed-modules:     +    AI.Clustering.Hierarchical+    AI.Clustering.Hierarchical.Types++  other-modules:       +    AI.Clustering.Hierarchical.Internal++  build-depends:+      base >=4.0 && <5.0+    , vector+    , containers++  hs-source-dirs:      src+  default-language:    Haskell2010++test-suite test+  type: exitcode-stdio-1.0+  hs-source-dirs: tests+  main-is: test.hs+  other-modules:+    Test.Hierarchical++  default-language:    Haskell2010+  build-depends: +      base+    , mwc-random+    , vector+    , tasty+    , tasty-hunit+    , clustering+    , hierarchical-clustering++benchmark bench+  type: exitcode-stdio-1.0+  hs-source-dirs: benchmarks+  main-is: bench.hs++  default-language:    Haskell2010+  build-depends: +      base+    , criterion+    , mwc-random+    , vector+    , clustering+    , hierarchical-clustering++source-repository  head+  type: git+  location: https://github.com/kaizhang/clustering.git
+ src/AI/Clustering/Hierarchical.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE FlexibleContexts #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  $Header$+-- Description :  <optional short text displayed on contents page>+-- Copyright   :  (c) Kai Zhang+-- License     :  MIT++-- Maintainer  :  kai@kzhang.org+-- Stability   :  experimental+-- Portability :  portable++-- <module description starting at first column>+--------------------------------------------------------------------------------++module AI.Clustering.Hierarchical+    ( Dendrogram(..)+    , size+    , cutAt+    , members+    , Metric(..)+    , hclust+    , computeDists+    , euclidean+    ) where++import Control.Applicative ((<$>))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++import AI.Clustering.Hierarchical.Internal+import AI.Clustering.Hierarchical.Types++data Metric = Single    -- ^ Single linkage, $d(A,B) = min_{a \in A, b \in B} d(a,b)$.+            | Complete  -- ^ Complete linkage, $d(A,B) = max_{a \in A, b \in B} d(a,b)$.+            | Average   -- ^ Average linkage, $d(A,B) = \frac{\sum_{a \in A}\sum_{b \in B}d(a,b)}{|A||B|}$.+            | Weighted  -- ^ Weighted linkage+            | Ward      -- ^ Ward's method+            | Centroid  -- ^ Centroid linkage, not implemented+            | Median    -- ^ Median linkage, not implemented++hclust :: G.Vector v a => Metric -> v a -> DistFn a -> Dendrogram a+hclust method xs f = label <$> nnChain dists fn+  where+    dists = computeDists f xs+    label i = xs G.! i+    fn = case method of+        Single -> single+        Complete -> complete+        Average -> average+        Weighted -> weighted+        Ward -> ward+        _ -> error "Not implemented"++computeDists :: G.Vector v a => DistFn a -> v a -> DistanceMat+computeDists f vec = DistanceMat n . U.fromList . flip concatMap [0..n-1] $ \i ->+    flip map [i+1..n-1] $ \j -> f (vec `G.unsafeIndex` i) (vec `G.unsafeIndex` j)+  where+    n = G.length vec+{-# INLINE computeDists #-}++euclidean :: G.Vector v Double => DistFn (v Double)+euclidean xs ys = sqrt $ G.sum $ G.zipWith (\x y -> (x-y)**2) xs ys+{-# INLINE euclidean #-}
+ src/AI/Clustering/Hierarchical/Internal.hs view
@@ -0,0 +1,122 @@+module AI.Clustering.Hierarchical.Internal+    ( nnChain+    , single+    , complete+    , average+    , weighted+    , ward+    ) where++import Control.Monad (forM_, when)+import qualified Data.Map as M+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM++import AI.Clustering.Hierarchical.Types++type ActiveNodeSet = M.Map Int (Dendrogram Int)+type DistUpdateFn = Int -> Int -> ActiveNodeSet -> DistanceMat -> DistanceMat++-- | nearest neighbor chain algorithm+nnChain :: DistanceMat -> DistUpdateFn -> Dendrogram Int+nnChain (DistanceMat n dist) fn = go (DistanceMat n $ U.force dist) initSet []+  where+    go ds activeNodes chain@(b:a:rest)+        | M.size activeNodes == 1 = head . M.elems $ activeNodes+        | c == a = go ds' activeNodes' rest+        | otherwise = go ds activeNodes $ c : chain+      where+        (c,d) = nearestNeighbor ds b a activeNodes+        activeNodes' = M.insert hi (Branch (size1+size2) d c1 c2)+                     . M.delete lo $ activeNodes+        ds' = fn lo hi activeNodes ds+        c1 = M.findWithDefault undefined lo activeNodes+        c2 = M.findWithDefault undefined hi activeNodes+        size1 = size c1+        size2 = size c2+        (lo,hi) = if a <= b then (a,b) else (b,a)+    go ds activeNodes _ = go ds activeNodes [b,a]+      where+        a = fst $ M.elemAt 0 activeNodes+        b = fst $ nearestNeighbor ds a (-1) activeNodes+    initSet = M.fromList . map (\i -> (i, Leaf i)) $ [0..n-1]+{-# INLINE nnChain #-}++nearestNeighbor :: DistanceMat -> Int -> Int -> M.Map Int (Dendrogram Int) -> (Int, Double)+nearestNeighbor dist i preference = M.foldlWithKey' f (-1,1/0)+  where+    f (x,d) j _ | i == j = (x,d)  -- skip+                | d' < d = (j,d')+                | d' == d && j == preference = (j,d')+                | otherwise = (x,d)+      where d' = dist ! (i,j)+{-# INLINE nearestNeighbor #-}++-- | all update functions perform destructive updates, and hence should not be+-- called outside this module++-- | single linkage update formula+single :: DistUpdateFn+single lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do+    v <- U.unsafeThaw dist+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do+        d_lo_i <- UM.unsafeRead v $ idx n i lo+        d_hi_i <- UM.unsafeRead v $ idx n i hi+        UM.unsafeWrite v (idx n i hi) $ min d_lo_i d_hi_i+    return v+{-# INLINE single #-}++-- | complete linkage update formula+complete :: DistUpdateFn+complete lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do+    v <- U.unsafeThaw dist+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do+        d_lo_i <- UM.unsafeRead v $ idx n i lo+        d_hi_i <- UM.unsafeRead v $ idx n i hi+        UM.unsafeWrite v (idx n i hi) $ max d_lo_i d_hi_i+    return v+{-# INLINE complete #-}++-- | average linkage update formula+average :: DistUpdateFn+average lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do+    v <- U.unsafeThaw dist+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do+        d_lo_i <- UM.unsafeRead v $ idx n i lo+        d_hi_i <- UM.unsafeRead v $ idx n i hi+        UM.unsafeWrite v (idx n i hi) $ f1 * d_lo_i + f2 * d_hi_i+    return v+  where+    s1 = fromIntegral . size . M.findWithDefault undefined lo $ nodeset+    s2 = fromIntegral . size . M.findWithDefault undefined hi $ nodeset+    f1 = s1 / (s1+s2)+    f2 = s2 / (s1+s2)+{-# INLINE average #-}++-- | complete linkage update formula+weighted :: DistUpdateFn+weighted lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do+    v <- U.unsafeThaw dist+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do+        d_lo_i <- UM.unsafeRead v $ idx n i lo+        d_hi_i <- UM.unsafeRead v $ idx n i hi+        UM.unsafeWrite v (idx n i hi) $ (d_lo_i + d_hi_i) / 2+    return v+{-# INLINE weighted #-}++-- | ward linkage update formula+ward :: DistUpdateFn+ward lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do+    v <- U.unsafeThaw dist+    d_lo_hi <- UM.unsafeRead v $ idx n lo hi+    forM_ (M.toList nodeset) $ \(i,t) -> when (i/= hi && i/=lo) $ do+        let s3 = fromIntegral . size $ t+        d_lo_i <- UM.unsafeRead v $ idx n i lo+        d_hi_i <- UM.unsafeRead v $ idx n i hi+        UM.unsafeWrite v (idx n i hi) $+            sqrt $ ((s1+s3)*d_lo_i + (s2+s3)*d_hi_i - s3*d_lo_hi) / (s1+s2+s3)+    return v+  where+    s1 = fromIntegral . size . M.findWithDefault undefined lo $ nodeset+    s2 = fromIntegral . size . M.findWithDefault undefined hi $ nodeset+{-# INLINE ward #-}
+ src/AI/Clustering/Hierarchical/Types.hs view
@@ -0,0 +1,55 @@+module AI.Clustering.Hierarchical.Types+    ( Distance+    , DistFn+    , Size+    , Dendrogram(..)+    , size+    , cutAt+    , members+    , DistanceMat(..)+    , (!)+    , idx+    ) where++import Data.Bits (shiftR)+import qualified Data.Vector.Unboxed as U++type Distance = Double+type DistFn a = a -> a -> Distance+type Size = Int++data Dendrogram a = Leaf !a+                  | Branch !Size !Distance !(Dendrogram a) !(Dendrogram a)+    deriving (Show)++instance Functor Dendrogram where+    fmap f (Leaf x) = Leaf $ f x+    fmap f (Branch n d l r) = Branch n d (fmap f l) $ fmap f r++size :: Dendrogram a -> Int+size (Leaf _) = 1+size (Branch n _ _ _) = n+{-# INLINE size #-}++cutAt :: Dendrogram a -> Distance -> [Dendrogram a]+cutAt dendro th = go [] dendro+  where+    go acc x@(Leaf _) = x : acc+    go acc x@(Branch _ d l r) | d <= th = x : acc+                              | otherwise = go (go acc r) l++members :: Dendrogram a -> [a]+members (Leaf x) = [x]+members (Branch _ _ l r) = members l ++ members r++-- upper triangular matrix+data DistanceMat = DistanceMat !Int !(U.Vector Double) deriving (Show)++(!) :: DistanceMat -> (Int, Int) -> Double+(!) (DistanceMat n v) (i',j') = v U.! idx n i' j'+{-# INLINE (!) #-}++idx :: Int -> Int -> Int -> Int+idx n i j | i <= j = (i * (2 * n - i - 3)) `shiftR` 1 + j - 1+          | otherwise = (j * (2 * n - j - 3)) `shiftR` 1 + i - 1+{-# INLINE idx #-}
+ tests/Test/Hierarchical.hs view
@@ -0,0 +1,64 @@+module Test.Hierarchical+    ( tests+    ) where++import Control.Monad+import qualified Data.Clustering.Hierarchical as C+import qualified Data.Vector as V+import System.Random.MWC+import Test.Tasty+import Test.Tasty.HUnit++import AI.Clustering.Hierarchical++tests :: TestTree+tests = testGroup "Hierarchical:"+    [ testCase "Single Linkage" testSingle+    , testCase "Complete Linkage" testComplete+    , testCase "Average Linkage" testAverage+    , testCase "Weighted Linkage" testWeighted+    ]++randSample :: IO [V.Vector Double]+randSample = do+    g <- create+    replicateM 500 $ uniformVector g 5++isEqual :: Eq a => Dendrogram a -> C.Dendrogram a -> Bool+isEqual (Leaf x) (C.Leaf x') = x == x'+isEqual (Branch _ d x y) (C.Branch d' x' y') = abs (d - d') < 1e-8 &&+    ((isEqual x x' && isEqual y y') || (isEqual x y' && isEqual y x'))+isEqual _ _ = False++testSingle :: Assertion+testSingle = do+    xs <- randSample+    let true = C.dendrogram C.SingleLinkage xs euclidean+        test = hclust Single (V.fromList xs) euclidean+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $+        isEqual test true++testComplete :: Assertion+testComplete = do+    xs <- randSample+    let true = C.dendrogram C.CompleteLinkage xs euclidean+        test = hclust Complete (V.fromList xs) euclidean+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $+        isEqual test true++testAverage :: Assertion+testAverage = do+    xs <- randSample+    let true = C.dendrogram C.UPGMA xs euclidean+        test = hclust Average (V.fromList xs) euclidean+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $+        isEqual test true++testWeighted :: Assertion+testWeighted = do+    xs <- randSample+    let true = C.dendrogram C.FakeAverageLinkage xs euclidean+        test = hclust Weighted (V.fromList xs) euclidean+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $+        isEqual test true+
+ tests/test.hs view
@@ -0,0 +1,7 @@+import Test.Tasty++import qualified Test.Hierarchical as Hierarchical++main :: IO ()+main = defaultMain $ testGroup "Main"+    [ Hierarchical.tests ]