packages feed

kmeans-par (empty) → 1.0.2

raw patch · 6 files changed

+283/−0 lines, 6 filesdep +basedep +criteriondep +deepseqsetup-changed

Dependencies added: base, criterion, deepseq, kmeans-par, normaldistribution, parallel, random, semigroups, split, vector

Files

+ LICENSE view
@@ -0,0 +1,21 @@+The MIT License (MIT)++Copyright (c) vi++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
+ benchmark/Main.lhs view
@@ -0,0 +1,39 @@+We aim to benchmark each implementation of Lloyd's algorithm:++> import Algorithms.Lloyd.Sequential (Point(..), Cluster(..))+> import qualified Algorithms.Lloyd.Sequential as Sequential (kmeans)+> import qualified Algorithms.Lloyd.Strategies as Strategies (kmeans)+> import Data.Random.Normal (mkNormals)+> import Data.Vector (fromList)+> import Control.Monad (forM)+> import Control.DeepSeq (NFData(..))+> import Criterion.Main (defaultMain, bench, nf)++We draw 2e3 normally distributed 2D points:++> points :: [Point]+> points = concat $ [0..2000-1] `forM` \n -> do+>   return . Point . fromList $ [normals!!n, normals!!(n*2)]+>   where normals = take 4000 $ mkNormals 0x29a++Three of which form the centroids of our initial clusters:+ +> clusters :: [Cluster]+> clusters = map (uncurry Cluster) $ zip [0..] (take 3 points)++To correctly benchmark the result of a pure function, we need to be able to+evaluate it to normal form:++> instance NFData Cluster where+>   rnf (Cluster i c) = rnf i `seq` rnf c+>+> instance NFData Point where+>   rnf (Point v) = rnf v+ +Together the subject of our benchmarks:+ +> main :: IO ()+> main = defaultMain+>   [ bench "Sequential" $ nf (Sequential.kmeans points) clusters+>   , bench "Strategies" $ nf (Strategies.kmeans 64 points) clusters+>   ]
+ kmeans-par.cabal view
@@ -0,0 +1,26 @@+name:                kmeans-par+version:             1.0.2+synopsis:            Sequential and parallel implementations of Lloyd's algorithm.+license:             MIT+license-file:        LICENSE+author:              vi+maintainer:          me@vikramverma.com+category:            Algorithm+build-type:          Simple+cabal-version:       >=1.10++library+  exposed-modules:     Algorithms.Lloyd.Sequential,+                       Algorithms.Lloyd.Strategies+  build-depends:       base == 2.*, vector, semigroups, parallel, split+  hs-source-dirs:      src+  ghc-options:         -threaded+  default-language:    Haskell2010++benchmark kmeans-benchmark+  type:             exitcode-stdio-1.0+  main-is:          Main.lhs+  hs-source-dirs:   benchmark+  build-depends:    base == 2.*, random, criterion, normaldistribution, kmeans-par, deepseq, vector+  ghc-options:      -threaded+  default-language: Haskell2010
+ src/Algorithms/Lloyd/Sequential.lhs view
@@ -0,0 +1,148 @@+A sequential implementation of Lloyd's algorithm for k-means clustering,+adapted from Marlow's _Parallel and Concurrent Programming in Haskell_:++> module Algorithms.Lloyd.Sequential where+> +> import Prelude hiding (zipWith, map, foldr1, replicate)+> import Control.Monad (guard, forM_)+> import Data.Function (on)+> import Data.List (minimumBy)+> import Data.Semigroup (Semigroup(..))+> import Data.Vector (Vector(..), toList, create, zipWith, map, foldr1, empty, replicate)+> import qualified Data.Vector.Mutable as MV (replicate, read, write)+ +A data point is represented by an n-dimensional real vector:++> data Point = Point (Vector Double) deriving (Eq, Show)++Point isn't parametrised, so we can't define a functor (or any functor-family)+instances. Instead, we define a specialised map, `pmap`:++> pmap :: (Double -> Double) -> Point -> Point+> pmap f (Point vector) = Point $ map f vector+ +The `Metric` typeclass, as defined here, is intended to contain types that+admit metrics (i.e. are metric spaces.) Instances can be defined in terms of +`distance` or the infix `<->`:++> class Metric a where+>   distance :: a -> a -> Double+>   distance = (<->)+>   (<->) :: a -> a -> Double+>   (<->) = distance++We define distance between data-points as the sum of squares of differences of+corresponding coordinates. The actual (Euclidean) distance is given by the+square root of this value, but as we're only interested in comparing distances,+we omit that function for efficiency:++> instance Metric Point where+>   Point v0 <-> Point v1 = foldr1 (+) . map (**2) $ zipWith (-) v0 v1+ +`Point`s may be combined using simple vector addition:++> instance Semigroup Point where+>   Point v0 <> Point v1 = Point $ zipWith (+) v0 v1++Clusters are represented by an identifier and a centroid:++> data Cluster = Cluster+>   { identifier :: Int+>   , centroid   :: Point+>   } deriving (Show)++Clusters are treated as equivalent if they share a centroid:++> instance Eq Cluster where+>   (Cluster _ c) == (Cluster _ d) = c == d++`PointSum` is a point sum; this is an intermediate type used to contain the sum+of points in a set when constructing new clusters. It consists of a count of+included points, as well as the vector sum of the constituent points:++> data PointSum = PointSum Int Point deriving (Show)++Two `PointSum`s can be combined by summing corresponding values:++> instance Semigroup PointSum where+>   PointSum c0 p0 <> PointSum c1 p1 = PointSum (c0+c1) (p0<>p1)++Without dependent types, we can't define a valid Monoid instance for PointSum,+as the identity varies with the dimensions of the point vector. But we can+construct one at runtime:++> emptyPointSum :: Int -> PointSum+> emptyPointSum length = PointSum 0 . Point $ replicate length 0++A point sum is constructed by incrementally adding points likeso:++> addPoint :: PointSum -> Point -> PointSum+> addPoint sum point = sum <> pure point+>   where pure = PointSum 1 ++A point sum can be turned into a cluster by computing its centroid, the average+of all points associated with the cluster:++> toCluster :: Int -> PointSum -> Cluster+> toCluster cid (PointSum count point) = Cluster+>   { identifier = cid+>   , centroid   = pmap (//count) point+>   }+>+> (//) :: Double -> Int -> Double+> x // y = x / fromIntegral y+         +After each iteration, points are associated with the cluster having the nearest+centroid:++> closestCluster :: [Cluster] -> Point -> Cluster+> closestCluster clusters point = fst . minimumBy (compare `on` snd) $ do+>   cluster <- clusters+>   return (cluster, point <-> centroid cluster)+>+> assign :: [Cluster] -> [Point] -> Vector PointSum+> assign clusters points = let nc = length clusters in create $ do+>   vector <- MV.replicate nc $ emptyPointSum nc+>   points `forM_` \point -> do+>     let cluster  = closestCluster clusters point +>         position = identifier cluster+>     sum <- MV.read vector position+>     MV.write vector position $! addPoint sum point+>   return vector+>+> makeNewClusters :: Vector PointSum -> [Cluster]+> makeNewClusters vector = do+>   (pointSum@(PointSum count _), index) <- zip (toList vector) [0..]+>   guard $ count > 0 -- We don't want an empty PointSum: amongst other +>                     -- things, this'd lead to division by zero when +>                     -- computing the centroid.+>   return $ toCluster index pointSum+>+> step :: [Cluster] -> [Point] -> [Cluster]+> step = makeNewClusters ..: assign+>+> (..:) :: (Functor f, Functor g) => (a -> b) -> f (g a) -> f (g b)+> (..:) = fmap . fmap+> +> infixr 8 ..:++The algorithm consists of iteratively finding the centroid of each existing+cluster, then reallocating points according to which centroid is closest, until+convergence. As the algorithm isn't guaranteed to converge, we cut execution if+convergence hasn't been observed after eighty iterations:++> kmeans :: [Point] -> [Cluster] -> [Cluster]+> kmeans = kmeans' 0+>+> kmeans' :: Int -> [Point] -> [Cluster] -> [Cluster]+> kmeans' iterations points clusters +>   | iterations >= expectDivergent = clusters+>   | clusters' == clusters         = clusters +>   | otherwise                     = kmeans' (succ iterations) points clusters'+>   where clusters' = step clusters points+>+> expectDivergent :: Int+> expectDivergent = 80++A note on initialisation: typically clusters are randomly assigned to data+points prior to the first update step.
+ src/Algorithms/Lloyd/Strategies.lhs view
@@ -0,0 +1,47 @@+A parallel implementation of Lloyd's algorithm for k-means clustering,+adapted from Marlow's _Parallel and Concurrent Programming in Haskell_.+Here we use Evaluation Strategies to parallelise the assignment of+points to clusters:++> module Algorithms.Lloyd.Strategies where+>+> import Prelude hiding (zipWith)+> import Control.Parallel.Strategies (Strategy(..), parTraversable, using, rseq)+> import Data.List.Split (chunksOf)+> import Data.Semigroup (Semigroup(..))+> import Data.Vector (Vector(..), zipWith)+> import Algorithms.Lloyd.Sequential (Cluster(..), Point(..), PointSum(..), makeNewClusters, assign, (..:), expectDivergent)++We can combine two vectors of some same type $t$ provided we know how to+combine two $t$s:++> instance Semigroup a => Semigroup (Vector a) where+>   (<>) = zipWith (<>)++Step is modified to, given a partitioned list of points, perform+classification in parallel:++> step :: [Cluster] -> [[Point]] -> [Cluster]+> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: map . assign +>+> with :: Strategy a -> a -> a+> with = flip using++This version of k-means takes an additional arguments -- the number of+partitions the set of points'll be divided into. This needn't equal the+number of processors: if there are more spark than cores, the runtime+can be trusted to schedule unallocated sparks so soon as a core becomes+available. That said: if there are too many work items, the overhead of+recombination may exceed the speed-up provided by parallellism; if there+are too few items, and those items vary in cost, some of our cores may+be unused for part of the computation.++> kmeans :: Int -> [Point] -> [Cluster] -> [Cluster]+> kmeans = kmeans' 0 ..: chunksOf+>+> kmeans' :: Int -> [[Point]] -> [Cluster] -> [Cluster]+> kmeans' iterations points clusters +>   | iterations >= expectDivergent = clusters+>   | clusters' == clusters         = clusters +>   | otherwise                     = kmeans' (succ iterations) points clusters'+>   where clusters' = step clusters points