diff --git a/Math/KMeans.hs b/Math/KMeans.hs --- a/Math/KMeans.hs +++ b/Math/KMeans.hs @@ -1,89 +1,170 @@ -{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} {- | Module : Math.KMeans -Copyright : (c) Alp Mestanogullari, Ville Tirronen, 2011-2012 +Copyright : (c) Alp Mestanogullari, Ville Tirronen, 2011-2014 License : BSD3 Maintainer : Alp Mestanogullari Stability : experimental -An implementation of the k-means clustering algorithm based on the efficient vector package. +An implementation of the k-means clustering algorithm based on the vector package. +The core functions of this module are 'kmeans' and 'kmeansWith'. See some examples +on . + -} +module Math.KMeans + ( -- * The meat of this package: 'kmeans' + kmeans + , kmeansWith -module Math.KMeans (kmeans, Point, Cluster(..), computeClusters) where + , -- * Types + Distance + , Clusters + , Cluster(..) + , Centroids + , -- * Misc. + partition + , euclidSq + , l1dist + , linfdist + ) where + +import Control.Monad.Identity import qualified Data.Vector.Unboxed as V import qualified Data.Vector as G import qualified Data.List as L import Data.Function (on) ---- * K-Means clustering algorithm +-- | A distance on vectors +type Distance = V.Vector Double -> V.Vector Double -> Double --- | Type holding an object of any type and its associated feature vector -type Point a = (V.Vector Double, a) +-- | The euclidean distance without taking the final square root +-- This would waste cycles without changing the behavior of the algorithm +euclidSq :: Distance +euclidSq v1 v2 = V.sum $ V.zipWith diffsq v1 v2 + where diffsq a b = (a-b)^(2::Int) +{-# INLINE euclidSq #-} --- | Type representing a cluster (group) of vectors by its center and an id -data Cluster = Cluster { - cid :: !Int, - center :: !(V.Vector Double) - } -- deriving (Show,Eq) +-- | L1 distance of two vectors: d(v1, v2) = sum on i of |v1_i - v2_i| +l1dist :: Distance +l1dist v1 v2 = V.sum $ V.zipWith diffabs v1 v2 + where diffabs a b = abs (a - b) +{-# INLINE l1dist #-} --- genVec = V.fromList `fmap` vectorOf 3 arbitrary --- genPts = (flip zip) [0..] `fmap` replicateM 10 genVec --- genClusters = do --- cs <- replicateM 5 genVec --- return (zipWith Cluster [0.. ] cs) --- --- prop_regroup = forAll genClusters $ \c -> --- forAll genPts $ \v -> --- s (regroupPoints c v) == s (regroupPoints' c v) --- where --- same xs = length (L.nub xs) == length xs --- s = map L.sort +-- | L-inf distance of two vectors: d(v1, v2) = max |v1_i - v2_i] +linfdist :: Distance +linfdist v1 v2 = V.maximum $ V.zipWith diffabs v1 v2 + where diffabs a b = abs (a - b) +{-# INLINE linfdist #-} +-- | This is what 'kmeans' hands you back. It's just a 'G.Vector' of clusters +-- that will hopefully be of length 'k'. +type Clusters a = G.Vector (Cluster a) -{-#INLINE distance#-} -distance :: Point a -> V.Vector Double -> Double -distance (u,_) v = V.sum $ V.zipWith (\a b -> (a - b)^2) u v +-- | This type is used internally by 'kmeans'. It represents our (hopefully) +-- @k@ centroids, obtained by computing the new centroids of a 'Cluster' +type Centroids = G.Vector (V.Vector Double) -partition :: Int -> [a] -> [[a]] -partition k vs = go vs - where go vs = case L.splitAt n vs of - (vs', []) -> [vs'] - (vs', vss) -> vs' : go vss - n = (length vs + k - 1) `div` k +-- | A 'Cluster' of points is just a list of points +newtype Cluster a = + Cluster { elements :: [a] -- ^ elements that belong to that cluster + } deriving (Eq, Show) -{-#INLINE computeClusters#-} -computeClusters :: [[V.Vector Double]] -> [Cluster] -computeClusters = zipWith Cluster [0..] . map f - where f (x:xs) = let (n, v) = L.foldl' (\(k, s) v' -> (k+1, V.zipWith (+) s v')) (1, x) xs - in V.map (\x -> x / (fromIntegral n)) v +clusterAdd :: Cluster a -> a -> Cluster a +clusterAdd (Cluster c) x = Cluster (x:c) -{-#INLINE regroupPoints#-} -regroupPoints :: forall a. [Cluster] -> [Point a] -> [[Point a]] -regroupPoints clusters points = L.filter (not.null) . G.toList . G.accum (flip (:)) (G.replicate (length clusters) []) . map closest $ points - where - closest p = (cid (L.minimumBy (compare `on` (distance p . center)) clusters),p) +emptyCluster :: Cluster a +emptyCluster = Cluster [] -regroupPoints' :: forall a. [Cluster] -> [Point a] -> [[Point a]] -regroupPoints' clusters points = go points - where go points = map (map snd) . L.groupBy ((==) `on` fst) . L.sortBy (compare `on` fst) $ map (\p -> (closest p, p)) points - closest p = cid $ L.minimumBy (compare `on` (distance p . center)) clusters +addCentroids :: V.Vector Double -> V.Vector Double -> V.Vector Double +addCentroids v1 v2 = V.zipWith (+) v1 v2 -kmeansStep :: [Point a] -> [[Point a]] -> [[Point a]] -kmeansStep points pgroups = regroupPoints (computeClusters . map (map fst) $ pgroups) points +-- | This is the current partitionning strategy used +-- by 'kmeans'. If we want @k@ clusters, we just +-- try to regroup consecutive elements in @k@ buckets +partition :: Int -> [a] -> Clusters a +partition k vs = G.fromList $ go vs + where go l = case L.splitAt n l of + (vs', []) -> [Cluster vs'] + (vs', vss) -> Cluster vs' : go vss + n = (length vs + k - 1) `div` k -kmeansAux :: [Point a] -> [[Point a]] -> [[Point a]] -kmeansAux points pgroups = let pss = kmeansStep points pgroups in - case map (map fst) pss == map (map fst) pgroups of - True -> pgroups - False -> kmeansAux points pss +-- | Run the kmeans clustering algorithm. +-- +-- > kmeans f distance k points +-- +-- will run the algorithm using 'f' to extract features from your type, +-- using 'distance' to measure the distance between vectors, +-- trying to separate 'points' in 'k' clusters. +-- +-- Extracting features just means getting a 'V.Vector' +-- with 'Double' coordinates that will represent your type +-- in the space in which 'kmeans' will run. +kmeans :: (a -> V.Vector Double) -- ^ feature extraction + -> Distance -- ^ distance function + -> Int -- ^ the 'k' to run 'k'-means with (i.e number of desired clusters) + -> [a] -- ^ input list of 'points' + -> Clusters a -- ^ result, hopefully 'k' clusters of points +kmeans extract dist k points = + runIdentity $ kmeansWith (\n ps -> return $ partition n ps) extract dist k points --- | Performs the k-means clustering algorithm --- using trying to use 'k' clusters on the given list of points -kmeans :: Int -> [Point a] -> [[Point a]] -kmeans k points = kmeansAux points pgroups - where pgroups = partition k points +-- | Same as 'kmeans', except that instead of using 'partition', you supply your own +-- function for choosing the initial clustering. Two important things to note: +-- +-- * If you don't need any kind of effect and just have a 'partition'-like function +-- you want to use, @m@ will can just be 'Identity' here. If that's too +-- obnoxious to work with, please let me know and I may just provide a separate +-- 'kmeansWith' function with no there. But most of the time, you'll probably just +-- be interested in the following scenario. +-- +-- * Most likely, you want to have something smarter than our simple 'partition' function. +-- A couple of papers I have read claim very decent results by using some precise +-- probabilistic schemas for the initial partitionning. In this case, your @m@ would +-- probably be 'IO' or 'ST' (e.g using my package) +-- and you could fine-tune the way the initial clusters are picked so that the algorithm +-- may give better results. Of course, if your initialization is monadic, so is the result. +kmeansWith :: Monad m + => (Int -> [a] -> m (Clusters a)) -- ^ how should we partition the points? + -> (a -> V.Vector Double) -- ^ get the coordinates of a "point" + -> Distance -- ^ what distance do we use + -> Int -- ^ number of desired clusters + -> [a] -- ^ list of points + -> m (Clusters a) -- ^ resulting clustering +kmeansWith initF extract dist k points = go `liftM` initF k points + + where + -- go :: Clusters a -> Clusters a + go pgroups = + case kmeansStep pgroups of + pgroups' | pgroupsEqualUnder pgroups pgroups' -> pgroups + | otherwise -> go pgroups' + -- kmeansStep :: Clusters a -> Clusters a + kmeansStep clusters = + case centroidsOf clusters of + centroids -> + G.filter (not . null . elements) + . G.unsafeAccum clusterAdd (G.replicate k emptyCluster) + . map (pairToClosestCentroid centroids) + $ points + -- centroidsOf :: Clusters a -> Centroids + centroidsOf cs = G.map centroidOf cs + where + n = fromIntegral $ G.length cs + + centroidOf (Cluster elts) = + V.map (/n) + . L.foldl1' addCentroids + $ map extract elts + + -- pairToClosestCentroid :: Centroids -> a -> (Int, a) + pairToClosestCentroid cs a = (minDistIndex, a) + where !minDistIndex = G.minIndexBy (compare `on` dist (extract a)) cs + + -- pgroupsEqualUnder :: Clusters a -> Clusters a -> Bool + pgroupsEqualUnder g1 g2 = + G.map (map extract . elements) g1 == G.map (map extract . elements) g2 +{-# INLINE kmeansWith #-} diff --git a/bench/OldKmeans.hs b/bench/OldKmeans.hs new file mode 100644 --- /dev/null +++ b/bench/OldKmeans.hs @@ -0,0 +1,93 @@ +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} + +{- | +Module : Math.KMeans +Copyright : (c) Alp Mestanogullari, Ville Tirronen, 2011-2014 +License : BSD3 +Maintainer : Alp Mestanogullari +Stability : experimental + +An implementation of the k-means clustering algorithm based on the efficient vector package. + +-} + +module OldKMeans (kmeans, Point, Cluster(..), computeClusters) where + +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector as G +import qualified Data.List as L +import Data.Function (on) + +--- * K-Means clustering algorithm + +-- | Type holding an object of any type and its associated feature vector +type Point a = (V.Vector Double, a) + +-- | Type representing a cluster (group) of vectors by its center and an id +data Cluster = Cluster { + cid :: {-# UNPACK #-} !Int, -- ^ an identifier for the cluster + center :: !(V.Vector Double) -- ^ the 'position' of the center of the cluster + } -- deriving (Show,Eq) + +-- genVec = V.fromList `fmap` vectorOf 3 arbitrary +-- genPts = (flip zip) [0..] `fmap` replicateM 10 genVec +-- genClusters = do +-- cs <- replicateM 5 genVec +-- return (zipWith Cluster [0.. ] cs) +-- +-- prop_regroup = forAll genClusters $ \c -> +-- forAll genPts $ \v -> +-- s (regroupPoints c v) == s (regroupPoints' c v) +-- where +-- same xs = length (L.nub xs) == length xs +-- s = map L.sort + + +{-# INLINE distance #-} +distance :: Point a -> V.Vector Double -> Double +distance (u,_) v = V.sum $ V.zipWith (\a b -> (a - b)^2) u v + +partition :: Int -> [a] -> [[a]] +partition k vs = go vs + where go vs = case L.splitAt n vs of + (vs', []) -> [vs'] + (vs', vss) -> vs' : go vss + n = (length vs + k - 1) `div` k + +{-#INLINE computeClusters#-} +computeClusters :: [[V.Vector Double]] -> [Cluster] +computeClusters = zipWith Cluster [0..] . map f + where f (x:xs) = let (n, v) = L.foldl' (\(k, s) v' -> (k+1, V.zipWith (+) s v')) (1, x) xs + in V.map (\x -> x / (fromIntegral n)) v + +{-#INLINE regroupPoints#-} +regroupPoints :: forall a. [Cluster] -> [Point a] -> [[Point a]] +regroupPoints clusters points = L.filter (not.null) . G.toList . G.accum (flip (:)) (G.replicate (length clusters) []) . map closest $ points + where + closest p = (cid (L.minimumBy (compare `on` (distance p . center)) clusters),p) + +regroupPoints' :: [Cluster] -> [Point a] -> [[Point a]] +regroupPoints' clusters points = go points + where go points = map (map snd) . L.groupBy ((==) `on` fst) . L.sortBy (compare `on` fst) $ map (\p -> (closest p, p)) points + closest p = cid $ L.minimumBy (compare `on` (distance p . center)) clusters + +kmeansStep :: [Point a] -> [[Point a]] -> [[Point a]] +kmeansStep points pgroups = + regroupPoints (computeClusters . map (map fst) $ pgroups) points + +kmeansAux :: [Point a] -> [[Point a]] -> [[Point a]] +kmeansAux points pgroups = let pss = kmeansStep points pgroups in + -- has anything changed since the last step? + -- even a point jumping from one cluster to another is enough to + -- enter the 'False' case + case map (map fst) pss == map (map fst) pgroups of + True -> pgroups -- nothing's changed, we're done + False -> kmeansAux points pss -- something has changed, so let's try again + +-- | Performs the k-means clustering algorithm +-- trying to use 'k' clusters on the given list of points +kmeans :: Int -> [Point a] -> [[Point a]] +kmeans k points = kmeansAux points pgroups + where pgroups = partition k points +{-# INLINE kmeans #-} + diff --git a/bench/bench.hs b/bench/bench.hs new file mode 100644 --- /dev/null +++ b/bench/bench.hs @@ -0,0 +1,71 @@ +module Main where + +import Control.Applicative +import Criterion.Main +import Test.QuickCheck + +import qualified Data.Vector as G +import qualified Data.Vector.Unboxed as V + +import qualified OldKMeans as K +import qualified Math.KMeans as K2 + +main :: IO () +main = do + persons1 <- generate persons + persons2 <- generate persons + + defaultMain + [ + bgroup "ints" [ bench "v0.2" $ whnf kmeans1 ints1 + , bench "v0.3" $ whnf kmeans2 ints2 + ] + , bgroup "persons" [ bench "v0.2" $ whnf kmeansP1 persons1 + , bench "v0.3" $ whnf kmeansP2 persons2] + ] + +ints1, ints2 :: [Int] +ints1 = [1..10000] +ints2 = [1..10000] + +data Person = Person + { age :: Int + , weight :: Double + , name :: String + , salary :: Int + } deriving (Eq, Show) + +instance Arbitrary Person where + arbitrary = do + Person <$> choose (2, 100) + <*> choose (5, 150) + <*> pure "francis" + <*> choose (500, 100000) + +persons :: Gen [Person] +persons = vector 10000 + +-- kmeans of 'Int's in 3 clusters +kmeans1 = G.fromList . K.kmeans 3 . map (\i -> (extract i, i)) +kmeans2 = K2.kmeans extract dist 3 + +-- kmeans of 'Person's in 4 clusters +kmeansP1 = G.fromList . K.kmeans 4 . map p2v + where p2v p = (personToVec p, p) +kmeansP2 = K2.kmeans personToVec eucl 4 + +personToVec :: Person -> V.Vector Double +personToVec p = V.fromList + [ fromIntegral $ age p + , weight p + , fromIntegral $ salary p + ] + +extract :: Int -> V.Vector Double +extract = V.singleton . fromIntegral + +dist :: K2.Distance +dist v1 v2 = V.sum $ V.zipWith (\x1 x2 -> abs (x1 - x2)) v1 v2 + +eucl :: K2.Distance +eucl v1 v2 = V.sum $ V.zipWith (\x1 x2 -> (x1 - x2)^2) v1 v2 diff --git a/examples/persons.hs b/examples/persons.hs new file mode 100644 --- /dev/null +++ b/examples/persons.hs @@ -0,0 +1,53 @@ +import Control.Applicative +import Control.Monad +import Math.KMeans +import Test.QuickCheck + +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector as G + +data Person = Person + { age :: Int + , weight :: Double + , name :: String + , salary :: Int + } deriving (Eq) + +instance Show Person where + show p = "<" ++ name p ++ ", " + ++ show (weight p) ++ "kg, " + ++ show (salary p) ++ "€/month, " + ++ show (age p) ++ "y.o>" + +instance Arbitrary Person where + arbitrary = do + Person <$> choose (2, 100) + <*> choose (5, 150) + <*> pure "francis" + <*> choose (500, 100000) + +persons :: Gen [Person] +persons = vector 5 + +d :: Distance +d v1 v2 = V.sum $ V.zipWith (\x1 x2 -> abs (x1 - x2)) v1 v2 + +personToVec :: Person -> V.Vector Double +personToVec p = V.fromList + [ fromIntegral $ age p + , weight p + , fromIntegral $ salary p + ] + +runKMeans :: [Person] -> Clusters Person +runKMeans = kmeans personToVec d 2 + +main :: IO () +main = do + ps <- generate persons + print ps + + let clusters = runKMeans ps + putStrLn $ show (G.length clusters) + ++ " cluster(s) found." + G.mapM_ print clusters diff --git a/kmeans-vector.cabal b/kmeans-vector.cabal --- a/kmeans-vector.cabal +++ b/kmeans-vector.cabal @@ -1,39 +1,48 @@ Name: kmeans-vector -Version: 0.2 +Version: 0.3 Synopsis: An implementation of the kmeans clustering algorithm based on the vector package -Description: Provides a simple (but efficient) implementation of the k-means clustering algorithm. The goal of this algorithm is to, given a list of n-dimensional points, regroup them in k groups, such that each point gets to be in the group to which it is the closest to (using the 'center' of the group). +Description: Provides a simple (but efficient) implementation of the k-means clustering algorithm. The goal of this algorithm is to, given a set of n-dimensional points, regroup them in k groups, such that each point gets to be in the group to which it is the closest to (using the 'center' of the group). . CHANGELOG . - kmeans-vector-0.2 supports having feature vectors associated to objects, and thus computing kmeans on these vectors, letting you recover the initial objects. - + 0.3: total rewrite of the code, the code scales much better on big inputs and is overall + consistently faster than the other kmeans implementations on hackage, on my laptop. + 0.2: supports having feature vectors associated to objects, and thus computing kmeans on these vectors, letting you recover the initial objects. Homepage: http://github.com/alpmestan/kmeans-vector - -Bug-reports: https://github.com/alpmestan/kmeans-vector/issues - +Bug-reports: https://github.com/alpmestan/kmeans-vector/issues License: BSD3 - License-file: LICENSE - Author: Alp Mestanogullari , Ville Tirronen - Maintainer: Alp Mestanogullari - -Copyright: 2011-2012 Alp Mestanogullari - -Stability: Experimental - +Copyright: 2011-2014 Alp Mestanogullari +Stability: Experimental Category: Math - Build-type: Simple +Cabal-version: >=1.8 -Cabal-version: >=1.6 +library + Exposed-modules: Math.KMeans + Build-depends: base >= 4 && < 5, vector >= 0.7, mtl >= 2.1 + ghc-prof-options: -prof -auto-all + ghc-options: -O2 -funbox-strict-fields -Wall -Library - Exposed-modules: Math.KMeans - Build-depends: base >= 4 && < 5, vector >= 0.7 - ghc-prof-options: -prof -auto-all - ghc-options: -O2 -funbox-strict-fields +executable kmeans-persons + main-is: persons.hs + hs-source-dirs: examples + ghc-options: -O2 -funbox-strict-fields + build-depends: base >= 4 && < 5, vector >= 0.7, kmeans-vector, QuickCheck + +benchmark bench + main-is: bench.hs + other-modules: OldKmeans + hs-source-dirs: bench + ghc-options: -O2 -funbox-strict-fields + type: exitcode-stdio-1.0 + build-depends: base >= 4 && < 5, + vector >= 0.7, + kmeans-vector, + criterion, + QuickCheck source-repository head type: git