diff --git a/benchmarks/Bench/Hierarchical.hs b/benchmarks/Bench/Hierarchical.hs new file mode 100644 --- /dev/null +++ b/benchmarks/Bench/Hierarchical.hs @@ -0,0 +1,44 @@ +module Bench.Hierarchical + ( benchHierarchical ) where + +import Criterion.Main +import qualified Data.Clustering.Hierarchical as C +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import System.IO.Unsafe (unsafePerformIO) + +import AI.Clustering.Hierarchical +import AI.Clustering.Hierarchical.Types + +import Bench.Utils + +benchHierarchical :: Benchmark +benchHierarchical = + let dists = computeDists euclidean xs + fn i j = dists ! (i,j) + xs = V.fromList $ unsafePerformIO $ randVectors 1000 5 + in bgroup "Hierarchical clustering" + [ bgroup "AI.Clustering.Hierarchical" + [ bench "Average Linkage (n = 10)" $ + whnf (\x -> hclust Average x fn) $! U.enumFromN 0 10 + , bench "Average Linkage (n = 100)" $ + whnf (\x -> hclust Average x fn) $! U.enumFromN 0 100 + , bench "Average Linkage (n = 500)" $ + whnf (\x -> hclust Average x fn) $! U.enumFromN 0 500 + ] + , bgroup "Data.Clustering.Hierarchical" + [ bench "Average Linkage (n = 10)" $ + whnf (\x -> C.dendrogram C.UPGMA x fn) $! [0..9] + , bench "Average Linkage (n = 100)" $ + whnf (\x -> C.dendrogram C.UPGMA x fn) $! [0..99] + , bench "Average Linkage (n = 500)" $ + whnf (\x -> C.dendrogram C.UPGMA x fn) $! [0..499] + ] + + , bgroup "Distance matrix" + [ bench "computeDists" $ + whnf (\x -> computeDists euclidean x) xs + , bench "computeDists'" $ + whnf (\x -> computeDists' euclidean x) xs + ] + ] diff --git a/benchmarks/Bench/KMeans.hs b/benchmarks/Bench/KMeans.hs new file mode 100644 --- /dev/null +++ b/benchmarks/Bench/KMeans.hs @@ -0,0 +1,31 @@ +module Bench.KMeans + ( benchKMeans ) where + +import Criterion.Main +import qualified Data.Matrix.Unboxed as MU +import qualified Data.Vector.Unboxed as U +import System.Random.MWC +import System.IO.Unsafe + +import AI.Clustering.KMeans + +import Bench.Utils + +gen :: GenIO +gen = unsafePerformIO createSystemRandom + +dat :: MU.Matrix Double +dat = unsafePerformIO $ fmap MU.fromRows $ randVectors 1000 10 + +benchKMeans :: Benchmark +benchKMeans = bgroup "KMeans clustering" + [ bgroup "AI.Clustering.KMeans" + [ bench "k-means++ (n = 1000, k = 7)" $ + whnfIO $ kmeans' gen KMeansPP 7 dat + , bench "forgy (n = 1000, k = 7)" $ + whnfIO $ kmeans' gen Forgy 7 dat + ] + ] + +kmeans' :: GenIO -> Method -> Int -> MU.Matrix Double -> IO (U.Vector Int) +kmeans' g method k = fmap _clusters . kmeans g method k diff --git a/benchmarks/Bench/Utils.hs b/benchmarks/Bench/Utils.hs new file mode 100644 --- /dev/null +++ b/benchmarks/Bench/Utils.hs @@ -0,0 +1,14 @@ +module Bench.Utils + ( randVectors + ) where + +import Control.Monad (replicateM) +import qualified Data.Vector.Unboxed as U +import System.Random.MWC + +randVectors :: Int -- ^ number of samples + -> Int -- ^ vector length + -> IO [U.Vector Double] +randVectors n k = do + g <- createSystemRandom + replicateM n $ uniformVector g k diff --git a/benchmarks/bench.hs b/benchmarks/bench.hs --- a/benchmarks/bench.hs +++ b/benchmarks/bench.hs @@ -7,32 +7,11 @@ import AI.Clustering.Hierarchical import AI.Clustering.Hierarchical.Types ((!)) -randSample :: IO [V.Vector Double] -randSample = do - g <- create - replicateM 2000 $ uniformVector g 5 +import Bench.Hierarchical (benchHierarchical) +import Bench.KMeans (benchKMeans) 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 - ] - ] +main = defaultMain + [ benchHierarchical + , benchKMeans + ] diff --git a/clustering.cabal b/clustering.cabal --- a/clustering.cabal +++ b/clustering.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: clustering -version: 0.1.2 +version: 0.2.0 synopsis: High performance clustering algorithms description: Following clutering methods are included in this library: @@ -26,22 +26,26 @@ library exposed-modules: AI.Clustering.Hierarchical + AI.Clustering.Hierarchical.Internal AI.Clustering.Hierarchical.Types AI.Clustering.KMeans + AI.Clustering.KMeans.Internal + AI.Clustering.KMeans.Types - other-modules: - AI.Clustering.Hierarchical.Internal +-- other-modules: build-depends: base >=4.0 && <5.0 , binary , containers - , matrices + , matrices >=0.4.0 , mwc-random + , parallel , primitive , vector hs-source-dirs: src + ghc-options: -Wall default-language: Haskell2010 test-suite test @@ -50,12 +54,15 @@ main-is: test.hs other-modules: Test.Hierarchical + Test.KMeans + Test.Utils default-language: Haskell2010 build-depends: base , binary , mwc-random + , matrices , vector , tasty , tasty-hunit @@ -63,11 +70,17 @@ , clustering , hierarchical-clustering , split + , Rlang-QQ benchmark bench type: exitcode-stdio-1.0 hs-source-dirs: benchmarks + ghc-options: -threaded -rtsopts -with-rtsopts=-N2 main-is: bench.hs + other-modules: + Bench.Hierarchical + Bench.KMeans + Bench.Utils default-language: Haskell2010 build-depends: @@ -77,6 +90,7 @@ , vector , clustering , hierarchical-clustering + , matrices source-repository head type: git diff --git a/src/AI/Clustering/Hierarchical.hs b/src/AI/Clustering/Hierarchical.hs --- a/src/AI/Clustering/Hierarchical.hs +++ b/src/AI/Clustering/Hierarchical.hs @@ -46,8 +46,10 @@ , cutAt , flatten , drawDendrogram - , computeDists + + -- * Distance functions , euclidean + , hamming -- * References -- $references @@ -56,7 +58,6 @@ import Control.Applicative ((<$>)) import qualified Data.Vector.Generic as G -import qualified Data.Vector.Unboxed as U import Text.Printf (printf) import AI.Clustering.Hierarchical.Internal @@ -75,7 +76,7 @@ hclust :: G.Vector v a => Linkage -> v a -> DistFn a -> Dendrogram a hclust method xs f = label <$> nnChain dists fn where - dists = computeDists f xs + dists = computeDists' f xs label i = xs G.! i fn = case method of Single -> single @@ -109,17 +110,15 @@ draw (Leaf x) = [x,""] shift first other = zipWith (++) (first : repeat other) -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 #-} - --- | compute euclidean distance between two points +-- | Compute euclidean distance between two points. 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 #-} + +-- | Hamming distance. +hamming :: (G.Vector v a, G.Vector v Bool, Eq a) => DistFn (v a) +hamming xs = fromIntegral . G.length . G.filter id . G.zipWith (/=) xs +{-# INLINE hamming #-} -- $references -- diff --git a/src/AI/Clustering/Hierarchical/Internal.hs b/src/AI/Clustering/Hierarchical/Internal.hs --- a/src/AI/Clustering/Hierarchical/Internal.hs +++ b/src/AI/Clustering/Hierarchical/Internal.hs @@ -1,4 +1,16 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : AI.Clustering.Hierarchical.Internal +-- Copyright : (c) 2015 Kai Zhang +-- License : MIT +-- +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable +-- +-------------------------------------------------------------------------------- module AI.Clustering.Hierarchical.Internal +{-# WARNING "To be used by developer only" #-} ( nnChain , single , complete diff --git a/src/AI/Clustering/Hierarchical/Types.hs b/src/AI/Clustering/Hierarchical/Types.hs --- a/src/AI/Clustering/Hierarchical/Types.hs +++ b/src/AI/Clustering/Hierarchical/Types.hs @@ -7,12 +7,16 @@ , DistanceMat(..) , (!) , idx + , computeDists + , computeDists' ) where import Control.Monad (liftM, liftM4) +import Control.Parallel.Strategies (rdeepseq, parMap) import Data.Binary (Binary, put, get, getWord8) import Data.Bits (shiftR) import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Generic as G import Data.Word (Word8) type Distance = Double @@ -59,3 +63,19 @@ 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 #-} + +-- | compute distance matrix +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 #-} + +-- | compute distance matrix in parallel +computeDists' :: G.Vector v a => DistFn a -> v a -> DistanceMat +computeDists' f vec = DistanceMat n . U.fromList . concat . flip (parMap rdeepseq) [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' #-} diff --git a/src/AI/Clustering/KMeans.hs b/src/AI/Clustering/KMeans.hs --- a/src/AI/Clustering/KMeans.hs +++ b/src/AI/Clustering/KMeans.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE FlexibleContexts #-} -------------------------------------------------------------------------------- -- | -- Module : AI.Clustering.KMeans @@ -12,54 +13,72 @@ module AI.Clustering.KMeans ( KMeans(..) , kmeans + , kmeansBy , kmeansWith -- * Initialization methods - , Initialization(..) + , Method(..) + -- * Useful functions , decode + , withinSS + + -- * References + -- $references ) where import Control.Monad (forM_) import Control.Monad.Primitive (PrimMonad, PrimState) import qualified Data.Matrix.Unboxed as MU +import qualified Data.Matrix.Generic as MG import qualified Data.Matrix.Unboxed.Mutable as MM import Data.Ord (comparing) import qualified Data.Vector as V +import qualified Data.Vector.Generic as G import qualified Data.Vector.Mutable as VM import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as UM -import Data.List (minimumBy, nub) -import System.Random.MWC (uniformR, Gen) +import Data.List (minimumBy, foldl') +import System.Random.MWC (Gen) --- | Results from running kmeans -data KMeans = KMeans - { _clusters :: U.Vector Int -- ^ A vector of integers (0 ~ k-1) - -- indicating the cluster to which each - -- point is allocated. - , _centers :: MU.Matrix Double -- ^ A matrix of cluster centres. - } deriving (Show) +import AI.Clustering.KMeans.Types (KMeans(..), Method(..)) +import AI.Clustering.KMeans.Internal (sumSquares, forgy, kmeansPP) --- | Lloyd's algorithm, also known as K-means algorithm -kmeans :: PrimMonad m +-- | Perform K-means clustering +kmeans :: (PrimMonad m, MG.Matrix mat U.Vector Double) => Gen (PrimState m) - -> Initialization - -> Int -- ^ number of clusters - -> MU.Matrix Double -- ^ each row represents a point + -> Method + -> Int + -> mat U.Vector Double -> m KMeans -kmeans g method k mat = do - initial <- case method of - Forgy -> forgy g k mat - _ -> undefined - return $ kmeansWith initial mat +kmeans g method k mat = kmeansBy g method k dat (MG.takeRow mat) + where + dat = U.enumFromN 0 $ MG.rows mat {-# INLINE kmeans #-} --- | Lloyd's algorithm, also known as K-means algorithm -kmeansWith :: MU.Matrix Double -- ^ initial set of k centroids - -> MU.Matrix Double -- ^ each row represents a point +-- | K-means algorithm +kmeansBy :: (PrimMonad m, G.Vector v a) + => Gen (PrimState m) + -> Method + -> Int -- ^ number of clusters + -> v a -- ^ data stores in rows + -> (a -> U.Vector Double) + -> m KMeans +kmeansBy g method k dat fn = do + initial <- case method of + Forgy -> forgy g k dat fn + KMeansPP -> kmeansPP g k dat fn + return $ kmeansWith initial dat fn +{-# INLINE kmeansBy #-} + +-- | K-means algorithm +kmeansWith :: G.Vector v a + => MU.Matrix Double -- ^ initial set of k centroids + -> v a -- ^ each row represents a point + -> (a -> U.Vector Double) -> KMeans -kmeansWith initial mat | d /= MU.cols initial || k > n = error "check input" - | otherwise = KMeans member centers +kmeansWith initial dat fn | d /= MU.cols initial || k > n = error "check input" + | otherwise = KMeans member centers where (member, centers) = loop initial U.empty loop means membership @@ -70,18 +89,20 @@ -- Assignment step assign means = U.generate n $ \i -> - let x = MU.takeRow mat i - in fst $ minimumBy (comparing snd) $ zip [0..k-1] $ map (dist x) $ MU.toRows means + let x = fn $ G.unsafeIndex dat i + in fst $ minimumBy (comparing snd) $ zip [0..k-1] $ map (sumSquares x) $ MU.toRows means -- Update step - update membership = MM.create $ do + update membership = MU.create $ do m <- MM.replicate (k,d) 0.0 count <- UM.replicate k (0 :: Int) forM_ [0..n-1] $ \i -> do let x = membership U.! i UM.unsafeRead count x >>= UM.unsafeWrite count x . (+1) + + let vec = fn $ dat G.! i forM_ [0..d-1] $ \j -> - MM.unsafeRead m (x,j) >>= MM.unsafeWrite m (x,j) . (+ mat MU.! (i,j)) + MM.unsafeRead m (x,j) >>= MM.unsafeWrite m (x,j) . (+ (vec U.! j)) -- normalize forM_ [0..k-1] $ \i -> do c <- UM.unsafeRead count i @@ -89,48 +110,11 @@ MM.unsafeRead m (i,j) >>= MM.unsafeWrite m (i,j) . (/fromIntegral c) return m - dist xs = U.sum . U.zipWith (\x y -> (x - y)**2) xs - - n = MU.rows mat + n = G.length dat k = MU.rows initial - d = MU.cols mat + d = MU.cols initial {-# INLINE kmeansWith #-} --- | Different initialization methods -data Initialization = Forgy -- ^ The Forgy method randomly chooses k unique - -- observations from the data set and uses these - -- as the initial means - | KMeansPP -- ^ K-means++ algorithm, not implemented. - -forgy :: PrimMonad m - => Gen (PrimState m) - -> Int -- number of clusters - -> MU.Matrix Double -- data - -> m (MU.Matrix Double) -forgy g k mat | k > n = error "k is larger than sample size" - | otherwise = iter - where - iter = do - vec <- sample g k . U.enumFromN 0 $ n - let xs = map (MU.takeRow mat) . U.toList $ vec - if length (nub xs) == length xs - then return . MU.fromRows $ xs - else iter - n = MU.rows mat -{-# INLINE forgy #-} - --- random select k samples from a population -sample :: PrimMonad m => Gen (PrimState m) -> Int -> U.Vector Int -> m (U.Vector Int) -sample g k xs = do - v <- U.thaw xs - forM_ [0..k-1] $ \i -> do - j <- uniformR (i, lst) g - UM.unsafeSwap v i j - U.unsafeFreeze . UM.take k $ v - where - lst = U.length xs - 1 -{-# INLINE sample #-} - -- | Assign data to clusters based on KMeans result decode :: KMeans -> [a] -> [[a]] decode result xs = V.toList $ V.create $ do @@ -143,4 +127,16 @@ n = U.maximum membership + 1 -- | Compute within-cluster sum of squares ---withinSS :: Matrix +withinSS :: KMeans -> MU.Matrix Double -> [Double] +withinSS result mat = zipWith f (decode result [0 .. MU.rows mat-1]) . + MU.toRows . _centers $ result + where + f c center = foldl' (+) 0 $ map (sumSquares center . MU.takeRow mat) c + + +-- $references +-- +-- Arthur, D. and Vassilvitskii, S. (2007). k-means++: the advantages of careful +-- seeding. Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete +-- algorithms. Society for Industrial and Applied Mathematics Philadelphia, PA, +-- USA. pp. 1027–1035. diff --git a/src/AI/Clustering/KMeans/Internal.hs b/src/AI/Clustering/KMeans/Internal.hs new file mode 100644 --- /dev/null +++ b/src/AI/Clustering/KMeans/Internal.hs @@ -0,0 +1,101 @@ +{-# LANGUAGE BangPatterns #-} +-------------------------------------------------------------------------------- +-- | +-- Module : AI.Clustering.KMeans.Internal +-- Copyright : (c) 2015 Kai Zhang +-- License : MIT +-- +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable +-- +-- +-------------------------------------------------------------------------------- +module AI.Clustering.KMeans.Internal +{-# WARNING "To be used by developer only" #-} + ( forgy + , kmeansPP + , sumSquares + ) where + +import Control.Monad (forM_) +import Control.Monad.Primitive (PrimMonad, PrimState) +import Data.List (nub) +import qualified Data.Matrix.Unboxed as MU +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as UM +import System.Random.MWC (uniformR, Gen) + +forgy :: (PrimMonad m, G.Vector v a) + => Gen (PrimState m) + -> Int -- number of clusters + -> v a -- data + -> (a -> U.Vector Double) + -> m (MU.Matrix Double) +forgy g k dat fn | k > n = error "k is larger than sample size" + | otherwise = iter + where + iter = do + vec <- randN g k . U.enumFromN 0 $ n + let xs = map (\i -> fn $ dat `G.unsafeIndex` i) . U.toList $ vec + if length (nub xs) == length xs + then return . MU.fromRows $ xs + else iter + n = G.length dat +{-# INLINE forgy #-} + +kmeansPP :: (PrimMonad m, G.Vector v a) + => Gen (PrimState m) + -> Int + -> v a + -> (a -> U.Vector Double) + -> m (MU.Matrix Double) +kmeansPP g k dat fn + | k > n = error "k is larger than sample size" + | otherwise = do + c1 <- uniformR (0,n-1) g + loop [c1] 1 + where + loop centers !k' + | k' == k = return $ MU.fromRows $ map (\i -> fn $ dat `G.unsafeIndex` i) centers + | otherwise = do + c' <- chooseWithProb g $ U.map (shortestDist centers) rowIndices + loop (c':centers) (k'+1) + + n = G.length dat + rowIndices = U.enumFromN 0 n + shortestDist centers x = minimum $ map (\i -> + sumSquares (fn $ dat `G.unsafeIndex` x) (fn $ dat `G.unsafeIndex` i)) centers +{-# INLINE kmeansPP #-} + +chooseWithProb :: PrimMonad m + => Gen (PrimState m) + -> U.Vector Double -- ^ weights, may not be normalized + -> m Int -- ^ result/index +chooseWithProb g ws = do + x <- uniformR (0,sum') g + return $ loop x 0 0 + where + loop v !cdf !i | cdf' >= v = i + | otherwise = loop v cdf' (i+1) + where cdf' = cdf + ws `U.unsafeIndex` i + + sum' = U.sum ws +{-# INLINE chooseWithProb #-} + +-- | Random select k samples from a population +randN :: PrimMonad m => Gen (PrimState m) -> Int -> U.Vector Int -> m (U.Vector Int) +randN g k xs = do + v <- U.thaw xs + forM_ [0..k-1] $ \i -> do + j <- uniformR (i, lst) g + UM.unsafeSwap v i j + U.unsafeFreeze . UM.take k $ v + where + lst = U.length xs - 1 +{-# INLINE randN #-} + +sumSquares :: U.Vector Double -> U.Vector Double -> Double +sumSquares xs = U.sum . U.zipWith (\x y -> (x - y)**2) xs +{-# INLINE sumSquares #-} diff --git a/src/AI/Clustering/KMeans/Types.hs b/src/AI/Clustering/KMeans/Types.hs new file mode 100644 --- /dev/null +++ b/src/AI/Clustering/KMeans/Types.hs @@ -0,0 +1,33 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : AI.Clustering.KMeans.Types +-- Copyright : (c) 2015 Kai Zhang +-- License : MIT +-- +-- Maintainer : kai@kzhang.org +-- Stability : experimental +-- Portability : portable +-- +-- +-------------------------------------------------------------------------------- +module AI.Clustering.KMeans.Types + ( KMeans(..) + , Method(..) + ) where + +import qualified Data.Matrix.Unboxed as MU +import qualified Data.Vector.Unboxed as U + +-- | Results from running kmeans +data KMeans = KMeans + { _clusters :: U.Vector Int -- ^ A vector of integers (0 ~ k-1) + -- indicating the cluster to which each + -- point is allocated. + , _centers :: MU.Matrix Double -- ^ A matrix of cluster centers. + } deriving (Show) + +-- | Different initialization methods +data Method = Forgy -- ^ The Forgy method randomly chooses k unique + -- observations from the data set and uses these + -- as the initial means. + | KMeansPP -- ^ K-means++ algorithm. diff --git a/tests/Test/Hierarchical.hs b/tests/Test/Hierarchical.hs --- a/tests/Test/Hierarchical.hs +++ b/tests/Test/Hierarchical.hs @@ -13,6 +13,7 @@ import Test.Tasty.QuickCheck import AI.Clustering.Hierarchical +import Test.Utils tests :: TestTree tests = testGroup "Hierarchical:" @@ -23,11 +24,6 @@ , 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 && @@ -36,7 +32,7 @@ testSingle :: Assertion testSingle = do - xs <- randSample + xs <- randVectors 500 5 let true = C.dendrogram C.SingleLinkage xs euclidean test = hclust Single (V.fromList xs) euclidean assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $ @@ -44,7 +40,7 @@ testComplete :: Assertion testComplete = do - xs <- randSample + xs <- randVectors 500 5 let true = C.dendrogram C.CompleteLinkage xs euclidean test = hclust Complete (V.fromList xs) euclidean assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $ @@ -52,7 +48,7 @@ testAverage :: Assertion testAverage = do - xs <- randSample + xs <- randVectors 500 5 let true = C.dendrogram C.UPGMA xs euclidean test = hclust Average (V.fromList xs) euclidean assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $ @@ -60,7 +56,7 @@ testWeighted :: Assertion testWeighted = do - xs <- randSample + xs <- randVectors 500 5 let true = C.dendrogram C.FakeAverageLinkage xs euclidean test = hclust Weighted (V.fromList xs) euclidean assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $ diff --git a/tests/Test/KMeans.hs b/tests/Test/KMeans.hs new file mode 100644 --- /dev/null +++ b/tests/Test/KMeans.hs @@ -0,0 +1,57 @@ +{-# LANGUAGE QuasiQuotes #-} +{-# LANGUAGE DataKinds #-} +module Test.KMeans + ( tests + ) where + +import Control.Monad +import qualified Data.Matrix.Unboxed as MU +import qualified Data.Vector.Unboxed as V +import Data.List +import RlangQQ +import System.Random.MWC +import Test.Tasty +import Test.Tasty.HUnit +import Test.Tasty.QuickCheck + +import AI.Clustering.KMeans +import AI.Clustering.KMeans.Internal + +import Test.Utils + +tests :: TestTree +tests = testGroup "KMeans:" + [ testCase "KMeans" testKMeans + ] + +rKmeans :: Int -> [Double] -> [Double] -> IO [Int] +rKmeans n dat center = do + o <- [r| x <- matrix(hs_dat, ncol=hs_n,byrow=T); + y <- matrix(hs_center, ncol=hs_n,byrow=T); + hs_result <- kmeans(x,y,iter.max=1000000,algorithm="Lloyd")$cluster; + |] + let x = Label :: Label "result" + return $ o .!. x + +testKMeans :: Assertion +testKMeans = do + let n = 2000 + d = 15 + k = 10 + g <- createSystemRandom + xs <- randVectors n d + + let mat = MU.fromRows xs :: MU.Matrix Double + dat = V.enumFromN 0 $ MU.rows mat + fn = MU.takeRow mat + + centers <- kmeansPP g k dat fn + + r <- rKmeans d (MU.toList mat) (MU.toList centers) + let test = sort $ map sort $ decode result xs + result = kmeansWith centers dat fn + true = sort $ map sort $ decode result{_clusters=V.fromList $ map (subtract 1) r} xs + show' xs = unlines $ map (show . map (unwords . map show . V.toList)) xs + + assertBool ("Expect: " ++ show' true ++ "\nBut saw: " ++ show' test) $ + test == true diff --git a/tests/Test/Utils.hs b/tests/Test/Utils.hs new file mode 100644 --- /dev/null +++ b/tests/Test/Utils.hs @@ -0,0 +1,14 @@ +module Test.Utils + ( randVectors + ) where + +import Control.Monad (replicateM) +import qualified Data.Vector.Unboxed as U +import System.Random.MWC + +randVectors :: Int -- ^ number of samples + -> Int -- ^ vector length + -> IO [U.Vector Double] +randVectors n k = do + g <- createSystemRandom + replicateM n $ uniformVector g k diff --git a/tests/test.hs b/tests/test.hs --- a/tests/test.hs +++ b/tests/test.hs @@ -1,7 +1,10 @@ import Test.Tasty import qualified Test.Hierarchical as Hierarchical +import qualified Test.KMeans as KMeans main :: IO () main = defaultMain $ testGroup "Main" - [ Hierarchical.tests ] + [ Hierarchical.tests + , KMeans.tests + ]