diff --git a/Data/Vector/Algorithms/Combinators.hs b/Data/Vector/Algorithms/Combinators.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Combinators.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE Rank2Types, TypeOperators #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Combinators +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Non-portable (rank-2 types) +-- +-- The purpose of this module is to supply various combinators for commonly +-- used idioms for the algorithms in this package. Examples at the time of +-- this writing include running an algorithm keyed on some function of the +-- elements (but only computing said function once per element), and safely +-- applying the algorithms on mutable arrays to immutable arrays. + +module Data.Vector.Algorithms.Combinators + ( apply +-- , usingKeys +-- , usingIxKeys + ) where + +import Prelude hiding (length) + +import Control.Monad.ST + +import Data.Ord + +import Data.Vector.Generic + +import qualified Data.Vector.Generic.Mutable as M +import qualified Data.Vector.Generic.New as N + +-- | Safely applies a mutable array algorithm to an immutable array. +apply :: (Vector v e) => (forall s mv. M.MVector mv e => mv s e -> ST s ()) -> v e -> v e +apply algo v = new . N.New $ + M.unstream (stream v) >>= \arr -> algo arr >> return arr + +{- +-- | Uses a function to compute a key for each element which the +-- algorithm should use in lieu of the actual element. For instance: +-- +-- > usingKeys sortBy f arr +-- +-- should produce the same results as: +-- +-- > sortBy (comparing f) arr +-- +-- the difference being that usingKeys computes each key only once +-- which can be more efficient for expensive key functions. +usingKeys :: (UA e, UA k, Ord k) + => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ()) + -> (e -> k) + -> MUArr e s + -> ST s () +usingKeys algo f arr = usingIxKeys algo (const f) arr +{-# INLINE usingKeys #-} + +-- | As usingKeys, only the key function has access to the array index +-- at which each element is stored. +usingIxKeys :: (UA e, UA k, Ord k) + => (forall e'. (UA e') => Comparison e' -> MUArr e' s -> ST s ()) + -> (Int -> e -> k) + -> MUArr e s + -> ST s () +usingIxKeys algo f arr = do + keys <- newMU (lengthMU arr) + fill len keys + algo (comparing fstS) (unsafeZipMU keys arr) + where + len = lengthMU arr + fill k keys + | k < 0 = return () + | otherwise = readMU arr k >>= writeMU keys k . f k >> fill (k-1) keys +{-# INLINE usingIxKeys #-} +-} diff --git a/Data/Vector/Algorithms/Common.hs b/Data/Vector/Algorithms/Common.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Common.hs @@ -0,0 +1,27 @@ + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Common +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Portable +-- +-- Common operations and utility functions for all sorts + +module Data.Vector.Algorithms.Common where + +import Prelude hiding (read, length) + +import Control.Monad.Primitive + +import Data.Vector.Generic.Mutable + +-- | A type of comparisons between two values of a given type. +type Comparison e = e -> e -> Ordering + +copyOffset :: (PrimMonad m, MVector v e) + => v (PrimState m) e -> v (PrimState m) e -> Int -> Int -> Int -> m () +copyOffset from to iFrom iTo len = + unsafeCopy (unsafeSlice iTo len to) (unsafeSlice iFrom len from) +{-# INLINE copyOffset #-} diff --git a/Data/Vector/Algorithms/Insertion.hs b/Data/Vector/Algorithms/Insertion.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Insertion.hs @@ -0,0 +1,81 @@ + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Insertion +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Portable +-- +-- A simple insertion sort. Though it's O(n^2), its iterative nature can be +-- beneficial for small arrays. It is used to sort small segments of an array +-- by some of the more heavy-duty, recursive algorithms. + +module Data.Vector.Algorithms.Insertion + ( sort + , sortBy + , sortByBounds + , sortByBounds' + , Comparison + ) where + + +import Prelude hiding (read, length) + +import Control.Monad.Primitive + +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +import qualified Data.Vector.Algorithms.Optimal as O + +-- | Sorts an entire array using the default comparison for the type +sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m () +sort = sortBy compare +{-# INLINE sort #-} + +-- | Sorts an entire array using a given comparison +sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m () +sortBy cmp a = sortByBounds cmp a 0 (length a) +{-# INLINE sortBy #-} + +-- | Sorts the portion of an array delimited by [l,u) +sortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +sortByBounds cmp a l u + | len < 2 = return () + | len == 2 = O.sort2ByOffset cmp a l + | len == 3 = O.sort3ByOffset cmp a l + | len == 4 = O.sort4ByOffset cmp a l + | otherwise = O.sort4ByOffset cmp a l >> sortByBounds' cmp a l (l + 4) u + where + len = u - l +{-# INLINE sortByBounds #-} + +-- | Sorts the portion of the array delimited by [l,u) under the assumption +-- that [l,m) is already sorted. +sortByBounds' :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +sortByBounds' cmp a l m u = sort m + where + sort i + | i < u = do v <- unsafeRead a i + insert cmp a l v i + sort (i+1) + | otherwise = return () +{-# INLINE sortByBounds' #-} + +-- Given a sorted array in [l,u), inserts val into its proper position, +-- yielding a sorted [l,u] +insert :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> e -> Int -> m () +insert cmp a l = loop + where + loop val j + | j <= l = unsafeWrite a l val + | otherwise = do e <- unsafeRead a (j - 1) + case cmp val e of + LT -> unsafeWrite a j e >> loop val (j - 1) + _ -> unsafeWrite a j val +{-# INLINE insert #-} diff --git a/Data/Vector/Algorithms/Intro.hs b/Data/Vector/Algorithms/Intro.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Intro.hs @@ -0,0 +1,200 @@ +{-# LANGUAGE TypeOperators, BangPatterns #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Intro +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Non-portable (type operators, bang patterns) +-- +-- This module implements various algorithms based on the introsort algorithm, +-- originally described by David R. Musser in the paper /Introspective Sorting +-- and Selection Algorithms/. It is also in widespread practical use, as the +-- standard unstable sort used in the C++ Standard Template Library. +-- +-- Introsort is at its core a quicksort. The version implemented here has the +-- following optimizations that make it perform better in practice: +-- +-- * Small segments of the array are left unsorted until a final insertion +-- sort pass. This is faster than recursing all the way down to +-- one-element arrays. +-- +-- * The pivot for segment [l,u) is chosen as the median of the elements at +-- l, u-1 and (u+l)/2. This yields good behavior on mostly sorted (or +-- reverse-sorted) arrays. +-- +-- * The algorithm tracks its recursion depth, and if it decides it is +-- taking too long (depth greater than 2 * lg n), it switches to a heap +-- sort to maintain O(n lg n) worst case behavior. (This is what makes the +-- algorithm introsort). + +module Data.Vector.Algorithms.Intro + ( -- * Sorting + sort + , sortBy + , sortByBounds + -- * Selecting + , select + , selectBy + , selectByBounds + -- * Partial sorting + , partialSort + , partialSortBy + , partialSortByBounds + , Comparison + ) where + +import Prelude hiding (read, length) + +import Control.Monad +import Control.Monad.Primitive + +import Data.Bits +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +import qualified Data.Vector.Algorithms.Insertion as I +import qualified Data.Vector.Algorithms.Optimal as O +import qualified Data.Vector.Algorithms.TriHeap as H + +-- | Sorts an entire array using the default ordering. +sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m () +sort = sortBy compare +{-# INLINE sort #-} + +-- | Sorts an entire array using a custom ordering. +sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m () +sortBy cmp a = sortByBounds cmp a 0 (length a) +{-# INLINE sortBy #-} + +-- | Sorts a portion of an array [l,u) using a custom ordering +sortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +sortByBounds cmp a l u + | len < 2 = return () + | len == 2 = O.sort2ByOffset cmp a l + | len == 3 = O.sort3ByOffset cmp a l + | len == 4 = O.sort4ByOffset cmp a l + | otherwise = introsort cmp a (ilg len) l u + where len = u - l +{-# INLINE sortByBounds #-} + +-- Internal version of the introsort loop which allows partial +-- sort functions to call with a specified bound on iterations. +introsort :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +introsort cmp a i l u = sort i l u >> I.sortByBounds cmp a l u + where + sort 0 l u = H.sortByBounds cmp a l u + sort d l u + | len < threshold = return () + | otherwise = do O.sort3ByIndex cmp a c l (u-1) -- sort the median into the lowest position + p <- unsafeRead a l + mid <- partitionBy cmp a p (l+1) u + unsafeSwap a l (mid - 1) + sort (d-1) mid u + sort (d-1) l (mid - 1) + where + len = u - l + c = (u + l) `div` 2 +{-# INLINE introsort #-} + +-- | Moves the least k elements to the front of the array in +-- no particular order. +select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m () +select = selectBy compare +{-# INLINE select #-} + +-- | Moves the least k elements (as defined by the comparison) to +-- the front of the array in no particular order. +selectBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +selectBy cmp a k = selectByBounds cmp a k 0 (length a) +{-# INLINE selectBy #-} + +-- | Moves the least k elements in the interval [l,u) to the positions +-- [l,k+l) in no particular order. +selectByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +selectByBounds cmp a k l u = go (ilg len) l (l + k) u + where + len = u - l + go 0 l m u = H.selectByBounds cmp a (m - l) l u + go n l m u = do O.sort3ByIndex cmp a c l (u-1) + p <- unsafeRead a l + mid <- partitionBy cmp a p (l+1) u + unsafeSwap a l (mid - 1) + if m > mid + then go (n-1) mid m u + else if m < mid - 1 + then go (n-1) l m (mid - 1) + else return () + where c = (u + l) `div` 2 +{-# INLINE selectByBounds #-} + +-- | Moves the least k elements to the front of the array, sorted. +partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m () +partialSort = partialSortBy compare +{-# INLINE partialSort #-} + +-- | Moves the least k elements (as defined by the comparison) to +-- the front of the array, sorted. +partialSortBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a) +{-# INLINE partialSortBy #-} + +-- | Moves the least k elements in the interval [l,u) to the positions +-- [l,k+l), sorted. +partialSortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +partialSortByBounds cmp a k l u = go (ilg len) l (l + k) u + where + len = u - l + go 0 l m n = H.partialSortByBounds cmp a (m - l) l u + go n l m u + | l == m = return () + | otherwise = do O.sort3ByIndex cmp a c l (u-1) + p <- unsafeRead a l + mid <- partitionBy cmp a p (l+1) u + unsafeSwap a l (mid - 1) + case compare m mid of + GT -> do introsort cmp a (n-1) l (mid - 1) + go (n-1) mid m u + EQ -> introsort cmp a (n-1) l m + LT -> go n l m (mid - 1) + where c = (u + l) `div` 2 +{-# INLINE partialSortByBounds #-} + +partitionBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int +partitionBy cmp a = partUp + where + partUp p l u + | l < u = do e <- unsafeRead a l + case cmp e p of + LT -> partUp p (l+1) u + _ -> partDown p l (u-1) + | otherwise = return l + partDown p l u + | l < u = do e <- unsafeRead a u + case cmp p e of + LT -> partDown p l (u-1) + _ -> unsafeSwap a l u >> partUp p (l+1) u + | otherwise = return l +{-# INLINE partitionBy #-} + +-- computes the number of recursive calls after which heapsort should +-- be invoked given the lower and upper indices of the array to be sorted +ilg :: Int -> Int +ilg m = 2 * loop m 0 + where + loop 0 !k = k - 1 + loop n !k = loop (n `shiftR` 1) (k+1) + +-- the size of array at which the introsort algorithm switches to insertion sort +threshold :: Int +threshold = 18 +{-# INLINE threshold #-} diff --git a/Data/Vector/Algorithms/Merge.hs b/Data/Vector/Algorithms/Merge.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Merge.hs @@ -0,0 +1,155 @@ +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Merge +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Portable +-- +-- This module implements a simple top-down merge sort. The temporary buffer +-- is preallocated to 1/2 the size of the input array, and shared through +-- the entire sorting process to ease the amount of allocation performed in +-- total. This is a stable sort. + +module Data.Vector.Algorithms.Merge + ( sort + , sortBy + , sortByBounds + , Comparison + ) where + +import Prelude hiding (read, length) + +import Control.Monad.Primitive + +import Data.Bits +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison, copyOffset) + +import qualified Data.Vector.Algorithms.Optimal as O +import qualified Data.Vector.Algorithms.Insertion as I + +-- | Sorts an array using the default comparison. +sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m () +sort = sortBy compare +{-# INLINE sort #-} + +-- | Sorts an array using a custom comparison. +sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m () +sortBy cmp vec = sortByBounds cmp vec 0 (length vec) +{- +sortBy cmp vec + | len <= 1 = return () + | len == 2 = O.sort2ByOffset cmp vec 0 + | len == 3 = O.sort3ByOffset cmp vec 0 + | len == 4 = O.sort4ByOffset cmp vec 0 + | otherwise = do buf <- new len + mergeSortWithBuf cmp vec buf + where + len = length vec +-} +{-# INLINE sortBy #-} + + +-- | Sorts a portion of an array [l,u) using a custom comparison. +sortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +sortByBounds cmp vec l u + | len <= 1 = return () + | len == 2 = O.sort2ByOffset cmp vec l + | len == 3 = O.sort3ByOffset cmp vec l + | len == 4 = O.sort4ByOffset cmp vec l + | otherwise = do tmp <- new size + mergeSortWithBuf cmp vec tmp l u + where + len = u - l + size = (u + l) `div` 2 - l +{-# INLINE sortByBounds #-} + +mergeSortWithBuf :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> v (PrimState m) e + -> Int -> Int -> m () +mergeSortWithBuf cmp arr tmp = loop + where + loop l u + | len < threshold = I.sortByBounds cmp arr l u + | otherwise = do loop l mid + loop mid u + merge cmp arr tmp l mid u + where + len = u - l + mid = (u + l) `shiftR` 1 +{-# INLINE mergeSortWithBuf #-} + +{- +mergeSortWithBuf :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> v (PrimState m) e -> m () +mergeSortWithBuf cmp src buf + | length src < threshold = I.sortByBounds cmp src 0 (length src) + | otherwise = do mergeSortWithBuf cmp srcL bufL + mergeSortWithBuf cmp srcU bufU + merge cmp src buf mid + where + len = length src + mid = len `shiftR` 1 + + srcL = unsafeSlice 0 mid src + srcU = unsafeSlice mid (len - mid) src + + bufL = unsafeSlice 0 mid buf + bufU = unsafeSlice mid (len - mid) buf +-} + +merge :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> v (PrimState m) e + -> Int -> Int -> Int -> m () +merge cmp arr tmp l m u = do copyOffset arr tmp l 0 uTmp + eTmp <- unsafeRead tmp 0 + eArr <- unsafeRead arr m + loop 0 eTmp m eArr l + where + uTmp = m - l + uArr = u + loop iTmp eTmp iArr eArr iIns + | iTmp >= uTmp = return () + | iArr >= uArr = copyOffset tmp arr iTmp iIns (uTmp - iTmp) + | otherwise = case cmp eArr eTmp of + LT -> do unsafeWrite arr iIns eArr + eArr <- unsafeRead arr (iArr+1) + loop iTmp eTmp (iArr+1) eArr (iIns+1) + _ -> do unsafeWrite arr iIns eTmp + eTmp <- unsafeRead tmp (iTmp+1) + loop (iTmp+1) eTmp iArr eArr (iIns+1) +{-# INLINE merge #-} + +{- +merge :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> v (PrimState m) e + -> Int -> m () +merge cmp src buf mid = do unsafeCopy tmp lower + eTmp <- unsafeRead tmp 0 + eUpp <- unsafeRead upper 0 + loop tmp 0 eTmp upper 0 eUpp 0 + where + lower = unsafeSlice 0 mid src + upper = unsafeSlice mid (length src - mid) src + tmp = unsafeSlice 0 mid buf + + loop low iLow eLow high iHigh eHigh iIns + | iLow >= length low = return () + | iHigh >= length high = unsafeCopy (unsafeSlice iIns (length low - iLow) src) + (unsafeSlice iLow (length low - iLow) low) + | otherwise = case cmp eHigh eLow of + LT -> do unsafeWrite src iIns eHigh + eHigh <- unsafeRead high (iHigh + 1) + loop low iLow eLow high (iHigh + 1) eHigh (iIns + 1) + _ -> do unsafeWrite src iIns eLow + eLow <- unsafeRead low (iLow + 1) + loop low (iLow + 1) eLow high iHigh eHigh (iIns + 1) +{-# INLINE merge #-} +-} + +threshold :: Int +threshold = 25 +{-# INLINE threshold #-} diff --git a/Data/Vector/Algorithms/Optimal.hs b/Data/Vector/Algorithms/Optimal.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Optimal.hs @@ -0,0 +1,244 @@ + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Optimal +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Portable +-- +-- Optimal sorts for very small array sizes, or for small numbers of +-- particular indices in a larger array (to be used, for instance, for +-- sorting a median of 3 values into the lowest position in an array +-- for a median-of-3 quicksort). + +-- The code herein was adapted from a C algorithm for optimal sorts +-- of small arrays. The original code was produced for the article +-- /Sorting Revisited/ by Paul Hsieh, available here: +-- +-- http://www.azillionmonkeys.com/qed/sort.html +-- +-- The LICENSE file contains the relevant copyright information for +-- the reference C code. + +module Data.Vector.Algorithms.Optimal + ( sort2ByIndex + , sort2ByOffset + , sort3ByIndex + , sort3ByOffset + , sort4ByIndex + , sort4ByOffset + , Comparison + ) where + +import Prelude hiding (read, length) + +import Control.Monad.Primitive + +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +#include "vector.h" + +-- | Sorts the elements at the positions 'off' and 'off + 1' in the given +-- array using the comparison. +sort2ByOffset :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +sort2ByOffset cmp a off = sort2ByIndex cmp a off (off + 1) +{-# INLINE sort2ByOffset #-} + +-- | Sorts the elements at the two given indices using the comparison. This +-- is essentially a compare-and-swap, although the first index is assumed to +-- be the 'lower' of the two. +sort2ByIndex :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +sort2ByIndex cmp a i j = UNSAFE_CHECK(checkIndex) "sort2ByIndex" i (length a) + $ UNSAFE_CHECK(checkIndex) "sort2ByIndex" j (length a) $ do + a0 <- unsafeRead a i + a1 <- unsafeRead a j + case cmp a0 a1 of + GT -> unsafeWrite a i a1 >> unsafeWrite a j a0 + _ -> return () +{-# INLINE sort2ByIndex #-} + +-- | Sorts the three elements starting at the given offset in the array. +sort3ByOffset :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +sort3ByOffset cmp a off = sort3ByIndex cmp a off (off + 1) (off + 2) +{-# INLINE sort3ByOffset #-} + +-- | Sorts the elements at the three given indices. The indices are assumed +-- to be given from lowest to highest, so if 'l < m < u' then +-- 'sort3ByIndex cmp a m l u' essentially sorts the median of three into the +-- lowest position in the array. +sort3ByIndex :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +sort3ByIndex cmp a i j k = UNSAFE_CHECK(checkIndex) "sort3ByIndex" i (length a) + $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" j (length a) + $ UNSAFE_CHECK(checkIndex) "sort3ByIndex" k (length a) $ do + a0 <- unsafeRead a i + a1 <- unsafeRead a j + a2 <- unsafeRead a k + case cmp a0 a1 of + GT -> case cmp a0 a2 of + GT -> case cmp a2 a1 of + LT -> do unsafeWrite a i a2 + unsafeWrite a k a0 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a2 + unsafeWrite a k a0 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a0 + _ -> case cmp a1 a2 of + GT -> case cmp a0 a2 of + GT -> do unsafeWrite a i a2 + unsafeWrite a j a0 + unsafeWrite a k a1 + _ -> do unsafeWrite a j a2 + unsafeWrite a k a1 + _ -> return () +{-# INLINE sort3ByIndex #-} + +-- | Sorts the four elements beginning at the offset. +sort4ByOffset :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +sort4ByOffset cmp a off = sort4ByIndex cmp a off (off + 1) (off + 2) (off + 3) +{-# INLINE sort4ByOffset #-} + +-- The horror... + +-- | Sorts the elements at the four given indices. Like the 2 and 3 element +-- versions, this assumes that the indices are given in increasing order, so +-- it can be used to sort medians into particular positions and so on. +sort4ByIndex :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> Int -> m () +sort4ByIndex cmp a i j k l = UNSAFE_CHECK(checkIndex) "sort4ByIndex" i (length a) + $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" j (length a) + $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" k (length a) + $ UNSAFE_CHECK(checkIndex) "sort4ByIndex" l (length a) $ do + a0 <- unsafeRead a i + a1 <- unsafeRead a j + a2 <- unsafeRead a k + a3 <- unsafeRead a l + case cmp a0 a1 of + GT -> case cmp a0 a2 of + GT -> case cmp a1 a2 of + GT -> case cmp a1 a3 of + GT -> case cmp a2 a3 of + GT -> do unsafeWrite a i a3 + unsafeWrite a j a2 + unsafeWrite a k a1 + unsafeWrite a l a0 + _ -> do unsafeWrite a i a2 + unsafeWrite a j a3 + unsafeWrite a k a1 + unsafeWrite a l a0 + _ -> case cmp a0 a3 of + GT -> do unsafeWrite a i a2 + unsafeWrite a j a1 + unsafeWrite a k a3 + unsafeWrite a l a0 + _ -> do unsafeWrite a i a2 + unsafeWrite a j a1 + unsafeWrite a k a0 + unsafeWrite a l a3 + _ -> case cmp a2 a3 of + GT -> case cmp a1 a3 of + GT -> do unsafeWrite a i a3 + unsafeWrite a j a1 + unsafeWrite a k a2 + unsafeWrite a l a0 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a3 + unsafeWrite a k a2 + unsafeWrite a l a0 + _ -> case cmp a0 a3 of + GT -> do unsafeWrite a i a1 + unsafeWrite a j a2 + unsafeWrite a k a3 + unsafeWrite a l a0 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a2 + unsafeWrite a k a0 + -- unsafeWrite a l a3 + _ -> case cmp a0 a3 of + GT -> case cmp a1 a3 of + GT -> do unsafeWrite a i a3 + -- unsafeWrite a j a1 + unsafeWrite a k a0 + unsafeWrite a l a2 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a3 + unsafeWrite a k a0 + unsafeWrite a l a2 + _ -> case cmp a2 a3 of + GT -> do unsafeWrite a i a1 + unsafeWrite a j a0 + unsafeWrite a k a3 + unsafeWrite a l a2 + _ -> do unsafeWrite a i a1 + unsafeWrite a j a0 + -- unsafeWrite a k a2 + -- unsafeWrite a l a3 + _ -> case cmp a1 a2 of + GT -> case cmp a0 a2 of + GT -> case cmp a0 a3 of + GT -> case cmp a2 a3 of + GT -> do unsafeWrite a i a3 + unsafeWrite a j a2 + unsafeWrite a k a0 + unsafeWrite a l a1 + _ -> do unsafeWrite a i a2 + unsafeWrite a j a3 + unsafeWrite a k a0 + unsafeWrite a l a1 + _ -> case cmp a1 a3 of + GT -> do unsafeWrite a i a2 + unsafeWrite a j a0 + unsafeWrite a k a3 + unsafeWrite a l a1 + _ -> do unsafeWrite a i a2 + unsafeWrite a j a0 + unsafeWrite a k a1 + -- unsafeWrite a l a3 + _ -> case cmp a2 a3 of + GT -> case cmp a0 a3 of + GT -> do unsafeWrite a i a3 + unsafeWrite a j a0 + -- unsafeWrite a k a2 + unsafeWrite a l a1 + _ -> do -- unsafeWrite a i a0 + unsafeWrite a j a3 + -- unsafeWrite a k a2 + unsafeWrite a l a1 + _ -> case cmp a1 a3 of + GT -> do -- unsafeWrite a i a0 + unsafeWrite a j a2 + unsafeWrite a k a3 + unsafeWrite a l a1 + _ -> do -- unsafeWrite a i a0 + unsafeWrite a j a2 + unsafeWrite a k a1 + -- unsafeWrite a l a3 + _ -> case cmp a1 a3 of + GT -> case cmp a0 a3 of + GT -> do unsafeWrite a i a3 + unsafeWrite a j a0 + unsafeWrite a k a1 + unsafeWrite a l a2 + _ -> do -- unsafeWrite a i a0 + unsafeWrite a j a3 + unsafeWrite a k a1 + unsafeWrite a l a2 + _ -> case cmp a2 a3 of + GT -> do -- unsafeWrite a i a0 + -- unsafeWrite a j a1 + unsafeWrite a k a3 + unsafeWrite a l a2 + _ -> do -- unsafeWrite a i a0 + -- unsafeWrite a j a1 + -- unsafeWrite a k a2 + -- unsafeWrite a l a3 + return () +{-# INLINE sort4ByIndex #-} diff --git a/Data/Vector/Algorithms/Radix.hs b/Data/Vector/Algorithms/Radix.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Radix.hs @@ -0,0 +1,282 @@ +{-# LANGUAGE ScopedTypeVariables, BangPatterns, TypeOperators #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Radix +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Non-portable (scoped type variables, bang patterns) +-- +-- This module provides a radix sort for a subclass of unboxed arrays. The +-- radix class gives information on +-- * the number of passes needed for the data type +-- +-- * the size of the auxiliary arrays +-- +-- * how to compute the pass-k radix of a value +-- +-- Radix sort is not a comparison sort, so it is able to achieve O(n) run +-- time, though it also uses O(n) auxiliary space. In addition, there is a +-- constant space overhead of 2*size*sizeOf(Int) for the sort, so it is not +-- advisable to use this sort for large numbers of very small arrays. +-- +-- A standard example (upon which one could base their own Radix instance) +-- is Word32: +-- +-- * We choose to sort on r = 8 bits at a time +-- +-- * A Word32 has b = 32 bits total +-- +-- Thus, b/r = 4 passes are required, 2^r = 256 elements are needed in an +-- auxiliary array, and the radix function is: +-- +-- > radix k e = (e `shiftR` (k*8)) .&. 256 + +module Data.Vector.Algorithms.Radix (sort, sortBy, Radix(..)) where + +import Prelude hiding (read, length) + +import Control.Monad +import Control.Monad.Primitive + +import qualified Data.Vector.Primitive.Mutable as PV +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +import Data.Bits +import Data.Int +import Data.Word + + +import Foreign.Storable + +class Radix e where + -- | The number of passes necessary to sort an array of es + passes :: e -> Int + -- | The size of an auxiliary array + size :: e -> Int + -- | The radix function parameterized by the current pass + radix :: Int -> e -> Int + +instance Radix Int where + passes _ = sizeOf (undefined :: Int) + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = e .&. 255 + radix i e + | i == passes e - 1 = radix' (e `xor` minBound) + | otherwise = radix' e + where radix' e = (e `shiftR` (i `shiftL` 3)) .&. 255 + {-# INLINE radix #-} + +instance Radix Int8 where + passes _ = 1 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix _ e = 255 .&. fromIntegral e `xor` 128 + {-# INLINE radix #-} + +instance Radix Int16 where + passes _ = 2 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral (((e `xor` minBound) `shiftR` 8) .&. 255) + {-# INLINE radix #-} + +instance Radix Int32 where + passes _ = 4 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255) + radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255) + radix 3 e = fromIntegral (((e `xor` minBound) `shiftR` 24) .&. 255) + {-# INLINE radix #-} + +instance Radix Int64 where + passes _ = 8 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255) + radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255) + radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255) + radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255) + radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255) + radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255) + radix 7 e = fromIntegral (((e `xor` minBound) `shiftR` 56) .&. 255) + {-# INLINE radix #-} + +instance Radix Word where + passes _ = sizeOf (undefined :: Word) + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix i e = fromIntegral ((e `shiftR` (i `shiftL` 3)) .&. 255) + {-# INLINE radix #-} + +instance Radix Word8 where + passes _ = 1 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix _ = fromIntegral + {-# INLINE radix #-} + +instance Radix Word16 where + passes _ = 2 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255) + {-# INLINE radix #-} + +instance Radix Word32 where + passes _ = 4 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255) + radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255) + radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255) + {-# INLINE radix #-} + +instance Radix Word64 where + passes _ = 8 + {-# INLINE passes #-} + size _ = 256 + {-# INLINE size #-} + radix 0 e = fromIntegral (e .&. 255) + radix 1 e = fromIntegral ((e `shiftR` 8) .&. 255) + radix 2 e = fromIntegral ((e `shiftR` 16) .&. 255) + radix 3 e = fromIntegral ((e `shiftR` 24) .&. 255) + radix 4 e = fromIntegral ((e `shiftR` 32) .&. 255) + radix 5 e = fromIntegral ((e `shiftR` 40) .&. 255) + radix 6 e = fromIntegral ((e `shiftR` 48) .&. 255) + radix 7 e = fromIntegral ((e `shiftR` 56) .&. 255) + {-# INLINE radix #-} + +instance (Radix i, Radix j) => Radix (i, j) where + passes ~(i, j) = passes i + passes j + {-# INLINE passes #-} + size ~(i, j) = size i `max` size j + {-# INLINE size #-} + radix k ~(i, j) | k < passes j = radix k j + | otherwise = radix (k - passes j) i + {-# INLINE radix #-} + +-- | Sorts an array based on the Radix instance. +sort :: forall e m v. (PrimMonad m, MVector v e, Radix e) + => v (PrimState m) e -> m () +sort arr = sortBy (passes e) (size e) radix arr + where + e :: e + e = undefined +{-# INLINE sort #-} + +-- | Radix sorts an array using custom radix information +-- requires the number of passes to fully sort the array, +-- the size of of auxiliary arrays necessary (should be +-- one greater than the maximum value returned by the radix +-- function), and a radix function, which takes the pass +-- and an element, and returns the relevant radix. +sortBy :: (PrimMonad m, MVector v e) + => Int -- ^ the number of passes + -> Int -- ^ the size of auxiliary arrays + -> (Int -> e -> Int) -- ^ the radix function + -> v (PrimState m) e -- ^ the array to be sorted + -> m () +sortBy passes size rdx arr = do + tmp <- new (length arr) + count <- new size + prefix <- new size + radixLoop passes rdx arr tmp count prefix +{-# INLINE sortBy #-} + +radixLoop :: (PrimMonad m, MVector v e) + => Int -- passes + -> (Int -> e -> Int) -- radix function + -> v (PrimState m) e -- array to sort + -> v (PrimState m) e -- temporary array + -> PV.MVector (PrimState m) Int -- radix count array + -> PV.MVector (PrimState m) Int -- placement array + -> m () +radixLoop passes rdx src dst count prefix = go False 0 + where + len = length src + go swap k + | k < passes = if swap + then body rdx dst src count prefix k >> go (not swap) (k+1) + else body rdx src dst count prefix k >> go (not swap) (k+1) + | otherwise = when swap (unsafeCopy src dst) +{-# INLINE radixLoop #-} + +body :: (PrimMonad m, MVector v e) + => (Int -> e -> Int) -- radix function + -> v (PrimState m) e -- source array + -> v (PrimState m) e -- destination array + -> PV.MVector (PrimState m) Int -- radix count + -> PV.MVector (PrimState m) Int -- placement + -> Int -- current pass + -> m () +body rdx src dst count prefix k = do + set count 0 + countLoop k rdx src count + unsafeWrite prefix 0 0 + prefixLoop count prefix + moveLoop k rdx src dst prefix +{-# INLINE body #-} + +countLoop :: (PrimMonad m, MVector v e) + => Int -> (Int -> e -> Int) + -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m () +countLoop k rdx src count = go 0 + where + len = length src + go i + | i < len = unsafeRead src i >>= inc count . rdx k >> go (i+1) + | otherwise = return () +{-# INLINE countLoop #-} + +prefixLoop :: (PrimMonad m) + => PV.MVector (PrimState m) Int -> PV.MVector (PrimState m) Int + -> m () +prefixLoop count prefix = go 1 0 + where + len = length count + go i pi + | i < len = do ci <- unsafeRead count (i-1) + let pi' = pi + ci + unsafeWrite prefix i pi' + go (i+1) pi' + | otherwise = return () +{-# INLINE prefixLoop #-} + +moveLoop :: (PrimMonad m, MVector v e) + => Int -> (Int -> e -> Int) -> v (PrimState m) e + -> v (PrimState m) e -> PV.MVector (PrimState m) Int -> m () +moveLoop k rdx src dst prefix = go 0 + where + len = length src + go i + | i < len = do srci <- unsafeRead src i + pf <- inc prefix (rdx k srci) + unsafeWrite dst pf srci + go (i+1) + | otherwise = return () +{-# INLINE moveLoop #-} + +inc :: (PrimMonad m) => PV.MVector (PrimState m) Int -> Int -> m Int +inc arr i = unsafeRead arr i >>= \e -> unsafeWrite arr i (e+1) >> return e +{-# INLINE inc #-} diff --git a/Data/Vector/Algorithms/Search.hs b/Data/Vector/Algorithms/Search.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/Search.hs @@ -0,0 +1,88 @@ +{-# LANGUAGE BangPatterns #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.Search +-- Copyright : (c) 2009-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Non-portable (bang patterns) +-- +-- This module implements several methods of searching for indicies to insert +-- elements into a sorted array. + +module Data.Vector.Algorithms.Search + ( binarySearch + , binarySearchBy + , binarySearchByBounds + , binarySearchL + , binarySearchLBy + , binarySearchLByBounds + , Comparison + ) where + +import Prelude hiding (read, length) + +import Control.Monad.Primitive + +import Data.Bits + +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +-- | Finds an index in a givesn sorted array at which the given element could +-- be inserted while maintaining the sortedness of the array. +binarySearch :: (PrimMonad m, MVector v e, Ord e) + => v (PrimState m) e -> e -> m Int +binarySearch = binarySearchBy compare + +-- | Finds an index in a given array, which must be sorted with respect to the +-- given comparison function, at which the given element could be inserted while +-- preserving the array's sortedness. +binarySearchBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> m Int +binarySearchBy cmp arr e = binarySearchByBounds cmp arr e 0 (length arr) + +-- | Given an array sorted with respect to a given comparison function in indices +-- in [l,u), finds an index in [l,u] at which the given element could be inserted +-- while preserving sortedness. +binarySearchByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int +binarySearchByBounds cmp arr e l u + | u <= l = return l + | otherwise = do e' <- unsafeRead arr k + case cmp e' e of + LT -> binarySearchByBounds cmp arr e (k+1) u + EQ -> return k + GT -> binarySearchByBounds cmp arr e l k + where k = (u + l) `shiftR` 1 +{-# INLINE binarySearchByBounds #-} + +-- | Finds the lowest index in a given sorted array at which the given element +-- could be inserted while maintaining the sortedness. +binarySearchL :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> e -> m Int +binarySearchL = binarySearchLBy compare +{-# INLINE binarySearchL #-} + +-- | Finds the lowest index in a given array, which must be sorted with respect to +-- the given comparison function, at which the given element could be inserted +-- while preserving the sortedness. +binarySearchLBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> m Int +binarySearchLBy cmp arr e = binarySearchLByBounds cmp arr e 0 (length arr) +{-# INLINE binarySearchLBy #-} + +-- | Given an array sorted with respect to a given comparison function on indices +-- in [l,u), finds the lowest index in [l,u] at which the given element could be +-- inserted while preserving sortedness. +binarySearchLByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> m Int +binarySearchLByBounds cmp arr e !l !u + | u <= l = return l + | otherwise = do e' <- unsafeRead arr k + case cmp e' e of + LT -> binarySearchLByBounds cmp arr e (k+1) u + _ -> binarySearchLByBounds cmp arr e l k + where k = (u + l) `shiftR` 1 +{-# INLINE binarySearchLByBounds #-} diff --git a/Data/Vector/Algorithms/TriHeap.hs b/Data/Vector/Algorithms/TriHeap.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Algorithms/TriHeap.hs @@ -0,0 +1,218 @@ +{-# LANGUAGE TypeOperators #-} + +-- --------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Algorithms.TriHeap +-- Copyright : (c) 2008-2010 Dan Doel +-- Maintainer : Dan Doel +-- Stability : Experimental +-- Portability : Non-portable (type operators) +-- +-- This module implements operations for working with a trinary heap stored +-- in an unboxed array. Most heapsorts are defined in terms of a binary heap, +-- in which each internal node has at most two children. By contrast, a +-- trinary heap has internal nodes with up to three children. This reduces +-- the number of comparisons in a heapsort slightly, and improves locality +-- (again, slightly) by flattening out the heap. + +module Data.Vector.Algorithms.TriHeap + ( -- * Sorting + sort + , sortBy + , sortByBounds + -- * Selection + , select + , selectBy + , selectByBounds + -- * Partial sorts + , partialSort + , partialSortBy + , partialSortByBounds + -- * Heap operations + , heapify + , pop + , popTo + , sortHeap + , Comparison + ) where + +import Prelude hiding (read, length) + +import Control.Monad +import Control.Monad.Primitive + +import Data.Vector.Generic.Mutable + +import Data.Vector.Algorithms.Common (Comparison) + +import qualified Data.Vector.Algorithms.Optimal as O + +-- | Sorts an entire array using the default ordering. +sort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> m () +sort = sortBy compare +{-# INLINE sort #-} + +-- | Sorts an entire array using a custom ordering. +sortBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> m () +sortBy cmp a = sortByBounds cmp a 0 (length a) +{-# INLINE sortBy #-} + +-- | Sorts a portion of an array [l,u) using a custom ordering +sortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +sortByBounds cmp a l u + | len < 2 = return () + | len == 2 = O.sort2ByOffset cmp a l + | len == 3 = O.sort3ByOffset cmp a l + | len == 4 = O.sort4ByOffset cmp a l + | otherwise = heapify cmp a l u >> sortHeap cmp a l (l+4) u >> O.sort4ByOffset cmp a l + where len = u - l +{-# INLINE sortByBounds #-} + +-- | Moves the lowest k elements to the front of the array. +-- The elements will be in no particular order. +select :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m () +select = selectBy compare +{-# INLINE select #-} + +-- | Moves the 'lowest' (as defined by the comparison) k elements +-- to the front of the array. The elements will be in no particular +-- order. +selectBy :: (PrimMonad m, MVector v e) => Comparison e -> v (PrimState m) e -> Int -> m () +selectBy cmp a k = selectByBounds cmp a k 0 (length a) +{-# INLINE selectBy #-} + +-- | Moves the 'lowest' k elements in the portion [l,u) of the +-- array into the positions [l,k+l). The elements will be in +-- no particular order. +selectByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +selectByBounds cmp a k l u + | l + k <= u = heapify cmp a l (l + k) >> go l (l + k) (u - 1) + | otherwise = return () + where + go l m u + | u < m = return () + | otherwise = do el <- unsafeRead a l + eu <- unsafeRead a u + case cmp eu el of + LT -> popTo cmp a l m u + _ -> return () + go l m (u - 1) +{-# INLINE selectByBounds #-} + +-- | Moves the lowest k elements to the front of the array, sorted. +partialSort :: (PrimMonad m, MVector v e, Ord e) => v (PrimState m) e -> Int -> m () +partialSort = partialSortBy compare +{-# INLINE partialSort #-} + +-- | Moves the lowest k elements (as defined by the comparison) to +-- the front of the array, sorted. +partialSortBy :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> m () +partialSortBy cmp a k = partialSortByBounds cmp a k 0 (length a) +{-# INLINE partialSortBy #-} + +-- | Moves the lowest k elements in the portion [l,u) of the array +-- into positions [l,k+l), sorted. +partialSortByBounds :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +partialSortByBounds cmp a k l u + -- this potentially does more work than absolutely required, + -- but using a heap to find the least 2 of 4 elements + -- seems unlikely to be better than just sorting all of them + -- with an optimal sort, and the latter is obviously index + -- correct. + | len < 2 = return () + | len == 2 = O.sort2ByOffset cmp a l + | len == 3 = O.sort3ByOffset cmp a l + | len == 4 = O.sort4ByOffset cmp a l + | u <= l + k = sortByBounds cmp a l u + | otherwise = do selectByBounds cmp a k l u + sortHeap cmp a l (l + 4) (l + k) + O.sort4ByOffset cmp a l + where + len = u - l +{-# INLINE partialSortByBounds #-} + +-- | Constructs a heap in a portion of an array [l, u) +heapify :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +heapify cmp a l u = loop $ (len - 1) `div` 3 + where + len = u - l + loop k + | k < 0 = return () + | otherwise = unsafeRead a (l+k) >>= \e -> + siftByOffset cmp a e l k len >> loop (k - 1) +{-# INLINE heapify #-} + +-- | Given a heap stored in a portion of an array [l,u), swaps the +-- top of the heap with the element at u and rebuilds the heap. +pop :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> m () +pop cmp a l u = popTo cmp a l u u +{-# INLINE pop #-} + +-- | Given a heap stored in a portion of an array [l,u) swaps the top +-- of the heap with the element at position t, and rebuilds the heap. +popTo :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +popTo cmp a l u t = do al <- unsafeRead a l + at <- unsafeRead a t + unsafeWrite a t al + siftByOffset cmp a at l 0 (u - l) +{-# INLINE popTo #-} + +-- | Given a heap stored in a portion of an array [l,u), sorts the +-- highest values into [m,u). The elements in [l,m) are not in any +-- particular order. +sortHeap :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m () +sortHeap cmp a l m u = loop (u-1) >> unsafeSwap a l m + where + loop k + | m < k = pop cmp a l k >> loop (k-1) + | otherwise = return () +{-# INLINE sortHeap #-} + +-- Rebuilds a heap with a hole in it from start downwards. Afterward, +-- the heap property should apply for [start + off, len + off). val +-- is the new value to be put in the hole. +siftByOffset :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> e -> Int -> Int -> Int -> m () +siftByOffset cmp a val off start len = sift val start len + where + sift val root len + | child < len = do (child', ac) <- maximumChild cmp a off child len + case cmp val ac of + LT -> unsafeWrite a (root + off) ac >> sift val child' len + _ -> unsafeWrite a (root + off) val + | otherwise = unsafeWrite a (root + off) val + where child = root * 3 + 1 +{-# INLINE siftByOffset #-} + +-- Finds the maximum child of a heap node, given the indx of the first child. +maximumChild :: (PrimMonad m, MVector v e) + => Comparison e -> v (PrimState m) e -> Int -> Int -> Int -> m (Int, e) +maximumChild cmp a off child1 len + | child3 < len = do ac1 <- unsafeRead a (child1 + off) + ac2 <- unsafeRead a (child2 + off) + ac3 <- unsafeRead a (child3 + off) + return $ case cmp ac1 ac2 of + LT -> case cmp ac2 ac3 of + LT -> (child3, ac3) + _ -> (child2, ac2) + _ -> case cmp ac1 ac3 of + LT -> (child3, ac3) + _ -> (child1, ac1) + | child2 < len = do ac1 <- unsafeRead a (child1 + off) + ac2 <- unsafeRead a (child2 + off) + return $ case cmp ac1 ac2 of + LT -> (child2, ac2) + _ -> (child1, ac1) + | otherwise = do ac1 <- unsafeRead a (child1 + off) ; return (child1, ac1) + where + child2 = child1 + 1 + child3 = child1 + 2 +{-# INLINE maximumChild #-} diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,65 @@ +Copyright (c) 2008-2010 Dan Doel + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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. + +------------------------------------------------------------------------------ + +The code in Data.Array.Vector.Algorithms.Mutable.Optimal is adapted from a C +algorithm for the same purpose. The folowing is the copyright notice for said +C code: + +Copyright (c) 2004 Paul Hsieh +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 sorttest nor the names of its 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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/bench/Blocks.hs b/bench/Blocks.hs new file mode 100644 --- /dev/null +++ b/bench/Blocks.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE Rank2Types #-} + +module Blocks where + +import Control.Monad +import Control.Monad.ST + +import Data.Vector.Unboxed.Mutable + +import System.CPUTime + +import System.Random.Mersenne + +-- Some conveniences for doing evil stuff in the ST monad. +-- All the tests get run in IO, but uvector stuff happens +-- in ST, so we temporarily coerce. +clock :: IO Integer +clock = getCPUTime + +-- Strategies for filling the initial arrays +rand :: (MTRandom e) => MTGen -> Int -> IO e +rand g _ = random g + +ascend :: Num e => Int -> IO e +ascend = return . fromIntegral + +descend :: Num e => e -> Int -> IO e +descend m n = return $ m - fromIntegral n + +modulo :: Integral e => e -> Int -> IO e +modulo m n = return $ fromIntegral n `mod` m + +-- This is the worst case for the median-of-three quicksort +-- used in the introsort implementation. +medianKiller :: Integral e => e -> Int -> IO e +medianKiller m n' + | n < k = return $ if even n then n + 1 else n + k + | otherwise = return $ (n - k + 1) * 2 + where + n = fromIntegral n' + k = m `div` 2 +{-# INLINE medianKiller #-} + +initialize :: (Unbox e) => MVector RealWorld e -> Int -> (Int -> IO e) -> IO () +initialize arr len fill = init $ len - 1 + where init n = fill n >>= unsafeWrite arr n >> when (n > 0) (init $ n - 1) +{-# INLINE initialize #-} + +speedTest :: (Unbox e) => Int + -> (Int -> IO e) + -> (MVector RealWorld e -> IO ()) + -> IO Integer +speedTest n fill algo = do + arr <- new n + initialize arr n fill + t0 <- clock + algo arr + t1 <- clock + return $ t1 - t0 +{-# INLINE speedTest #-} + + diff --git a/bench/LICENSE b/bench/LICENSE new file mode 100644 --- /dev/null +++ b/bench/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2009 Dan Doel + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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. diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -0,0 +1,188 @@ +{-# LANGUAGE Rank2Types #-} + +module Main (main) where + +import Prelude hiding (read, length) +import qualified Prelude as P + +import Control.Monad.ST +import Control.Monad.Error + +import Data.Char +import Data.Ord (comparing) +import Data.List (maximumBy) + +import Data.Vector.Unboxed.Mutable + +import qualified Data.Vector.Algorithms.Insertion as INS +import qualified Data.Vector.Algorithms.Intro as INT +import qualified Data.Vector.Algorithms.TriHeap as TH +import qualified Data.Vector.Algorithms.Merge as M +import qualified Data.Vector.Algorithms.Radix as R + +import System.Environment +import System.Console.GetOpt +import System.Random.Mersenne + +import Blocks + +-- Does nothing. For testing the speed/heap allocation of the building blocks. +noalgo :: (Unbox e) => MVector RealWorld e -> IO () +noalgo _ = return () + +-- Allocates a temporary buffer, like mergesort for similar purposes as noalgo. +alloc :: (Unbox e) => MVector RealWorld e -> IO () +alloc arr | len <= 4 = arr `seq` return () + | otherwise = (new (len `div` 2) :: IO (MVector RealWorld Int)) >> return () + where len = length arr + +displayTime :: String -> Integer -> IO () +displayTime s elapsed = putStrLn $ + s ++ " : " ++ show (fromIntegral elapsed / 1e12) ++ " seconds" + +run :: String -> IO Integer -> IO () +run s t = t >>= displayTime s + +sortSuite :: String -> MTGen -> Int -> (MVector RealWorld Int -> IO ()) -> IO () +sortSuite str g n sort = do + putStrLn $ "Testing: " ++ str + run "Random " $ speedTest n (rand g >=> modulo n) sort + run "Sorted " $ speedTest n ascend sort + run "Reverse-sorted " $ speedTest n (descend n) sort + run "Random duplicates " $ speedTest n (rand g >=> modulo 1000) sort + let m = 4 * (n `div` 4) + run "Median killer " $ speedTest m (medianKiller m) sort + +partialSortSuite :: String -> MTGen -> Int -> Int + -> (MVector RealWorld Int -> Int -> IO ()) -> IO () +partialSortSuite str g n k sort = sortSuite str g n (\a -> sort a k) + +-- ----------------- +-- Argument handling +-- ----------------- + +data Algorithm = DoNothing + | Allocate + | InsertionSort + | IntroSort + | IntroPartialSort + | IntroSelect + | TriHeapSort + | TriHeapPartialSort + | TriHeapSelect + | MergeSort + | RadixSort + deriving (Show, Read, Enum, Bounded) + +data Options = O { algos :: [Algorithm], elems :: Int, portion :: Int, usage :: Bool } deriving (Show) + +defaultOptions :: Options +defaultOptions = O [] 10000 1000 False + +type OptionsT = Options -> Either String Options + +options :: [OptDescr OptionsT] +options = [ Option ['A'] ["algorithm"] (ReqArg parseAlgo "ALGO") + ("Specify an algorithm to be run. Options:\n" ++ algoOpts) + , Option ['n'] ["num-elems"] (ReqArg parseN "INT") + "Specify the size of arrays in algorithms." + , Option ['k'] ["portion"] (ReqArg parseK "INT") + "Specify the number of elements to partial sort/select in\nrelevant algorithms." + , Option ['?','v'] ["help"] (NoArg $ \o -> Right $ o { usage = True }) + "Show options." + ] + where + allAlgos :: [Algorithm] + allAlgos = [minBound .. maxBound] + algoOpts = fmt allAlgos + fmt (x:y:zs) = '\t' : pad (show x) ++ show y ++ "\n" ++ fmt zs + fmt [x] = '\t' : show x ++ "\n" + fmt [] = "" + size = (" " ++) . maximumBy (comparing P.length) . map show $ allAlgos + pad str = zipWith const (str ++ repeat ' ') size + +parseAlgo :: String -> Options -> Either String Options +parseAlgo "None" o = Right $ o { algos = [] } +parseAlgo "All" o = Right $ o { algos = [DoNothing .. RadixSort] } +parseAlgo s o = leftMap (\e -> "Unrecognized algorithm `" ++ e ++ "'") + . fmap (\v -> o { algos = v : algos o }) $ readEither s + +leftMap :: (a -> b) -> Either a c -> Either b c +leftMap f (Left a) = Left (f a) +leftMap _ (Right c) = Right c + +parseNum :: (Int -> Options) -> String -> Either String Options +parseNum f = leftMap (\e -> "Invalid numeric argument `" ++ e ++ "'") . fmap f . readEither + +parseN, parseK :: String -> Options -> Either String Options +parseN s o = parseNum (\n -> o { elems = n }) s +parseK s o = parseNum (\k -> o { portion = k }) s + +readEither :: Read a => String -> Either String a +readEither s = case reads s of + [(x,t)] | all isSpace t -> Right x + _ -> Left s + +runTest :: MTGen -> Int -> Int -> Algorithm -> IO () +runTest g n k alg = case alg of + DoNothing -> sortSuite "no algorithm" g n noalgo + Allocate -> sortSuite "allocate" g n alloc + InsertionSort -> sortSuite "insertion sort" g n insertionSort + IntroSort -> sortSuite "introsort" g n introSort + IntroPartialSort -> partialSortSuite "partial introsort" g n k introPSort + IntroSelect -> partialSortSuite "introselect" g n k introSelect + TriHeapSort -> sortSuite "tri-heap sort" g n triHeapSort + TriHeapPartialSort -> partialSortSuite "partial tri-heap sort" g n k triHeapPSort + TriHeapSelect -> partialSortSuite "tri-heap select" g n k triHeapSelect + MergeSort -> sortSuite "merge sort" g n mergeSort + RadixSort -> sortSuite "radix sort" g n radixSort + _ -> putStrLn $ "Currently unsupported algorithm: " ++ show alg + +mergeSort :: MVector RealWorld Int -> IO () +mergeSort v = M.sort v +{-# NOINLINE mergeSort #-} + +introSort :: MVector RealWorld Int -> IO () +introSort v = INT.sort v +{-# NOINLINE introSort #-} + +introPSort :: MVector RealWorld Int -> Int -> IO () +introPSort v k = INT.partialSort v k +{-# NOINLINE introPSort #-} + +introSelect :: MVector RealWorld Int -> Int -> IO () +introSelect v k = INT.select v k +{-# NOINLINE introSelect #-} + +triHeapSort :: MVector RealWorld Int -> IO () +triHeapSort v = TH.sort v +{-# NOINLINE triHeapSort #-} + +triHeapPSort :: MVector RealWorld Int -> Int -> IO () +triHeapPSort v k = TH.partialSort v k +{-# NOINLINE triHeapPSort #-} + +triHeapSelect :: MVector RealWorld Int -> Int -> IO () +triHeapSelect v k = TH.select v k +{-# NOINLINE triHeapSelect #-} + +insertionSort :: MVector RealWorld Int -> IO () +insertionSort v = INS.sort v +{-# NOINLINE insertionSort #-} + +radixSort :: MVector RealWorld Int -> IO () +radixSort v = R.sort v +{-# NOINLINE radixSort #-} + +main :: IO () +main = do args <- getArgs + gen <- getStdGen + case getOpt Permute options args of + (fs, _, []) -> case foldl (>>=) (Right defaultOptions) fs of + Left err -> putStrLn $ usageInfo err options + Right opts | not (usage opts) -> + mapM_ (runTest gen (elems opts) (portion opts)) (algos opts) + | otherwise -> putStrLn $ usageInfo "uvector-algorithms-bench" options + (_, _, errs) -> putStrLn $ usageInfo (concat errs) options + + diff --git a/bench/RadSieve.hs b/bench/RadSieve.hs new file mode 100644 --- /dev/null +++ b/bench/RadSieve.hs @@ -0,0 +1,97 @@ +-- ------------------------------------------------------------------ +-- +-- Module : RadSieve +-- Copyright : (c) 2009 Dan Doel +-- +-- ------------------------------------------------------------------ +-- An implementation of a radical sieve, inspired by solving Project +-- Euler problem #124. +-- +-- Reproduction fo the problem text: +-- +-- The radical of n, rad(n), is the product of distinct prime factors +-- of n. For example, 504 = 23 × 32 × 7, so rad(504) = 2 × 3 × 7 = 42. +-- +-- If we calculate rad(n) for 1 ≤ n ≤ 10, then sort them on rad(n), +-- and sorting on n if the radical values are equal, we get: +-- +-- Unsorted Sorted +-- n rad(n) n rad(n) k +-- 1 1 1 1 1 +-- 2 2 2 2 2 +-- 3 3 4 2 3 +-- 4 2 8 2 4 +-- 5 5 3 3 5 +-- 6 6 9 3 6 +-- 7 7 5 5 7 +-- 8 2 6 6 8 +-- 9 3 7 7 9 +-- 10 10 10 10 10 +-- +-- Let E(k) be the kth element in the sorted n column; for example, +-- E(4) = 8 and E(6) = 9. +-- +-- If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000). + +module RadSieve where + +import Control.Monad +import Control.Monad.ST + +import Data.Array.Vector + +-- Radicals can be sieved as follows: +-- set a[1,n] = 1 +-- for i from 2 to n +-- if a[i] == 1 -- i must be prime +-- then a[j*i] *= i for positive integers j, j*i <= n +-- else do nothing -- i is composite, so its prime factors +-- -- have been accounted for +-- +-- This sieves for radicals up to the given integer. +radSieve :: Int -> ST s (MUArr Int s) +radSieve n = do arr <- newMU (n + 1) + fill arr n + sieve arr 1 + return arr + where + fill arr i | i < 0 = return () + | otherwise = writeMU arr i 1 >> fill arr (i-1) + sieve arr i | n < i = return () + | otherwise = do e <- readMU arr i + when (e == 1) $ mark arr i i + sieve arr (i+1) + mark arr p j | n < j = return () + | otherwise = readMU arr j >>= writeMU arr j . (*p) + >> mark arr p (j+p) + +-- Computes the answer to the above Project Euler problem. The correct +-- answer is only generated for a stable sorting function. +stableSortedRad :: Int -> Int + -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) + -> Int +stableSortedRad n k sortBy = runST (do rads <- radSieve n + index <- newMU (n + 1) + fillUp index n + sortBy (comparing fstS) + (unsafeZipMU rads index) + readMU k index) + where + fillUp arr k | k < 0 = return () + | otherwise = writeMU arr k k >> fillUp arr (k-1) + +-- Computes the answer to the above Project Euler problem. This version +-- will generate the correct answer even for unstable sorts, but may be +-- marginally slower. +unstableSortedRad :: Int -> Int + -> (forall s e. UA e => Comparison e -> MUArr e s -> ST s ()) + -> Int +unstableSortedRad n k sortBy = runST (do rads <- radSieve n + index <- newMU (n + 1) + fillUp index n + sortBy compare (unsafeZipMU rads index) + readMU k index) + where + fillUp arr k | k < 0 = return () + | otherwise = writeMU arr k k >> fillUp arr (k-1) + diff --git a/bench/vector-algorithms-bench.cabal b/bench/vector-algorithms-bench.cabal new file mode 100644 --- /dev/null +++ b/bench/vector-algorithms-bench.cabal @@ -0,0 +1,22 @@ +name: vector-algorithms-bench +version: 0.3 +license: BSD3 +license-file: LICENSE +author: Dan Doel +maintainer: Dan Doel +homepage: http://code.haskell.org/~doio/ +category: Benchmark +synopsis: Benchmarks for vector-algorithms +description: A suite of various benchmarks for verifying the + performance of the algorithms in vector-algorithms. +build-type: Simple +cabal-version: >= 1.2 + +executable vec-bench + build-depends: base, mersenne-random, vector, vector-algorithms, mtl + + ghc-options: -Wall -Odph + main-is: Main.hs + + extensions: + Rank2Types diff --git a/include/vector.h b/include/vector.h new file mode 100644 --- /dev/null +++ b/include/vector.h @@ -0,0 +1,31 @@ +#define PHASE_STREAM [1] +#define PHASE_INNER [0] + +#define INLINE_STREAM INLINE PHASE_STREAM +#define INLINE_INNER INLINE PHASE_INNER + +#ifndef NOT_VECTOR_MODULE +import qualified Data.Vector.Internal.Check as Ck +#endif + +#define ERROR(f) (Ck.f __FILE__ __LINE__) +#define ASSERT (Ck.assert __FILE__ __LINE__) +#define ENSURE (Ck.f __FILE__ __LINE__) +#define CHECK(f) (Ck.f __FILE__ __LINE__) + +#define BOUNDS_ERROR(f) (ERROR(f) Ck.Bounds) +#define BOUNDS_ASSERT (ASSERT Ck.Bounds) +#define BOUNDS_ENSURE (ENSURE Ck.Bounds) +#define BOUNDS_CHECK(f) (CHECK(f) Ck.Bounds) + +#define UNSAFE_ERROR(f) (ERROR(f) Ck.Unsafe) +#define UNSAFE_ASSERT (ASSERT Ck.Unsafe) +#define UNSAFE_ENSURE (ENSURE Ck.Unsafe) +#define UNSAFE_CHECK(f) (CHECK(f) Ck.Unsafe) + +#define INTERNAL_ERROR(f) (ERROR(f) Ck.Internal) +#define INTERNAL_ASSERT (ASSERT Ck.Internal) +#define INTERNAL_ENSURE (ENSURE Ck.Internal) +#define INTERNAL_CHECK(f) (CHECK(f) Ck.Internal) + + diff --git a/tests/Optimal.hs b/tests/Optimal.hs new file mode 100644 --- /dev/null +++ b/tests/Optimal.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE TypeOperators, FlexibleContexts #-} + +-- Exhaustive test sets for proper sorting and stability of +-- optimal sorts + +module Optimal where + +import Control.Arrow +import Control.Monad + +import Data.List +import Data.Function + +import Data.Vector.Generic hiding (map, zip, concatMap, (++), replicate) + +interleavings :: [a] -> [a] -> [[a]] +interleavings [ ] ys = [ys] +interleavings xs [ ] = [xs] +interleavings xs@(x:xt) ys@(y:yt) = map (x:) (interleavings xt ys) + ++ map (y:) (interleavings xs yt) + +monotones :: Int -> Int -> [[Int]] +monotones k = atLeastOne 0 + where + atLeastOne i 0 = [[]] + atLeastOne i n = map (i:) $ picks i (n-1) + picks _ 0 = [[]] + picks i n | i >= k = [replicate n k] + | otherwise = map (i:) (picks i (n-1)) ++ atLeastOne (i+1) n + + +stability :: (Vector v (Int,Int)) => Int -> [v (Int, Int)] +stability n = concatMap ( map fromList + . foldM interleavings [] + . groupBy ((==) `on` fst) + . flip zip [0..]) + $ monotones (n-2) n + +sort2 :: (Vector v Int) => [v Int] +sort2 = map fromList $ permutations [0,1] + +stability2 :: (Vector v (Int,Int)) => [v (Int, Int)] +stability2 = [fromList [(0, 0), (0, 1)]] + +sort3 :: (Vector v Int) => [v Int] +sort3 = map fromList $ permutations [0..2] + +{- +stability3 :: [UArr (Int :*: Int)] +stability3 = map toU [ [0:*:0, 0:*:1, 0:*:2] + , [0:*:0, 0:*:1, 1:*:2] + , [0:*:0, 1:*:2, 0:*:1] + , [1:*:2, 0:*:0, 0:*:1] + , [0:*:0, 1:*:1, 1:*:2] + , [1:*:1, 0:*:0, 1:*:2] + , [1:*:1, 1:*:2, 0:*:0] + ] +-} + +sort4 :: (Vector v Int) => [v Int] +sort4 = map fromList $ permutations [0..3] + diff --git a/tests/Properties.hs b/tests/Properties.hs new file mode 100644 --- /dev/null +++ b/tests/Properties.hs @@ -0,0 +1,167 @@ +{-# LANGUAGE RankNTypes, FlexibleContexts #-} + +module Properties where + +import Prelude + +import Optimal + +import Control.Monad +import Control.Monad.ST + +import Data.List +import Data.Ord + +import Data.Vector (Vector) +import qualified Data.Vector as V + +import Data.Vector.Mutable (MVector) +import qualified Data.Vector.Mutable as MV + +import qualified Data.Vector.Generic.Mutable as G + +import Data.Vector.Algorithms.Optimal (Comparison) +import Data.Vector.Algorithms.Radix (radix, passes, size) +import Data.Vector.Algorithms.Combinators + +import qualified Data.Map as M + +import Test.QuickCheck + +import Util + +prop_sorted :: (Ord e) => Vector e -> Property +prop_sorted arr | V.length arr < 2 = property True + | otherwise = check (V.head arr) (V.tail arr) + where + check e arr | V.null arr = property True + | otherwise = e <= V.head arr .&. check (V.head arr) (V.tail arr) + +prop_fullsort :: (Ord e) + => (forall s mv. G.MVector mv e => mv s e -> ST s ()) -> Vector e -> Property +prop_fullsort algo arr = prop_sorted $ apply algo arr + +{- +prop_schwartzian :: (UA e, UA k, Ord k) + => (e -> k) + -> (forall e s. (UA e) => (e -> e -> Ordering) -> MUArr e s -> ST s ()) + -> UArr e -> Property +prop_schwartzian f algo arr + | lengthU arr < 2 = property True + | otherwise = let srt = apply (algo `usingKeys` f) arr + in check (headU srt) (tailU srt) + where + check e arr | nullU arr = property True + | otherwise = f e <= f (headU arr) .&. check (headU arr) (tailU arr) +-} + +longGen :: (Arbitrary e) => Int -> Gen (Vector e) +longGen k = liftM2 (\l r -> V.fromList (l ++ r)) (vectorOf k arbitrary) arbitrary + +sanity :: Int +sanity = 100 + +prop_partialsort :: (Ord e, Arbitrary e, Show e) + => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ()) + -> Positive Int -> Property +prop_partialsort = prop_sized $ \algo k -> + prop_sorted . V.take k . apply algo + +prop_select :: (Ord e, Arbitrary e, Show e) + => (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ()) + -> Positive Int -> Property +prop_select = prop_sized $ \algo k arr -> + let vec' = apply algo arr + l = V.slice 0 k vec' + r = V.slice k (V.length vec' - k) vec' + in V.all (\e -> V.all (e <=) r) l + +prop_sized :: (Arbitrary e, Show e, Testable prop) + => ((forall s mv. G.MVector mv e => mv s e -> ST s ()) + -> Int -> Vector e -> prop) + -> (forall s mv. G.MVector mv e => mv s e -> Int -> ST s ()) + -> Positive Int -> Property +prop_sized prop algo (Positive k) = + let k' = k `mod` sanity + in forAll (longGen k') $ prop (\marr -> algo marr k') k' + +prop_stable :: (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> ST s ()) + -> Vector Int -> Property +-- prop_stable algo arr = property $ apply algo arr == arr +prop_stable algo arr = stable $ apply (algo (comparing fst)) $ V.zip arr ix + where + ix = V.fromList [1 .. V.length arr] + +stable arr | V.null arr = property True + | otherwise = let (e, i) = V.head arr + in V.all (\(e', i') -> e < e' || i < i') (V.tail arr) + .&. stable (V.tail arr) + +prop_stable_radix :: (forall e s mv. G.MVector mv e => Int -> Int -> (Int -> e -> Int) + -> mv s e -> ST s ()) + -> Vector Int -> Property +prop_stable_radix algo arr = + stable . apply (algo (passes e) (size e) (\k (e, _) -> radix k e)) + $ V.zip arr ix + where + ix = V.fromList [1 .. V.length arr] + e = V.head arr + +prop_optimal :: Int + -> (forall e s mv. G.MVector mv e => Comparison e -> mv s e -> Int -> ST s ()) + -> Property +prop_optimal n algo = label "sorting" sortn .&. label "stability" stabn + where + arrn = V.fromList [0..n-1] + sortn = all ( (== arrn) + . apply (\a -> algo compare a 0) + . V.fromList) + $ permutations [0..n-1] + stabn = all ( (== arrn) + . snd + . V.unzip + . apply (\a -> algo (comparing fst) a 0)) + $ stability n + +type Bag e = M.Map e Int + +toBag :: (Ord e) => Vector e -> Bag e +toBag = M.fromListWith (+) . flip zip (repeat 1) . V.toList + +prop_permutation :: (Ord e) => (forall s mv. G.MVector mv e => mv s e -> ST s ()) + -> Vector e -> Property +prop_permutation algo arr = property $ + toBag arr == toBag (apply algo arr) + +newtype SortedVec e = Sorted (Vector e) + +instance (Show e) => Show (SortedVec e) where + show (Sorted a) = show a + +instance (Arbitrary e, Ord e) => Arbitrary (SortedVec e) where + arbitrary = fmap (Sorted . V.fromList . sort) + $ liftM2 (++) (vectorOf 20 arbitrary) arbitrary + +ixRanges :: Vector e -> Gen (Int, Int) +ixRanges vec = do i <- fmap (`mod` len) arbitrary + j <- fmap (`mod` len) arbitrary + return $ if i < j then (i, j) else (j, i) + where len = V.length vec + +prop_search_inrange :: (Ord e) + => (forall s. MVector s e -> e -> Int -> Int -> ST s Int) + -> SortedVec e -> e -> Property +prop_search_inrange algo (Sorted arr) e = forAll (ixRanges arr) $ \(i, j) -> + let k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e i j) + in property $ i <= k && k <= j + where + len = V.length arr + +prop_search_lowbound :: (Ord e) + => (forall s. MVector s e -> e -> ST s Int) + -> SortedVec e -> e -> Property +prop_search_lowbound algo (Sorted arr) e = property $ (k == 0 || arr V.! (k-1) < e) + && (k == len || arr V.! k >= e) + where + len = V.length arr + k = runST (mfromList (V.toList arr) >>= \marr -> algo marr e) diff --git a/tests/Tests.hs b/tests/Tests.hs new file mode 100644 --- /dev/null +++ b/tests/Tests.hs @@ -0,0 +1,151 @@ +{-# LANGUAGE ImpredicativeTypes, RankNTypes, TypeOperators, FlexibleContexts #-} + +module Main (main) where + +import Properties + +import Util + +import Test.QuickCheck + +import Control.Monad +import Control.Monad.ST + +import Data.Int +import Data.Word + +import Data.Vector (Vector) +import qualified Data.Vector as V + +import Data.Vector.Generic.Mutable (MVector) +import qualified Data.Vector.Generic.Mutable as MV + +import Data.Vector.Algorithms.Combinators + +import qualified Data.Vector.Algorithms.Insertion as INS +import qualified Data.Vector.Algorithms.Intro as INT +import qualified Data.Vector.Algorithms.Merge as M +import qualified Data.Vector.Algorithms.Radix as R +import qualified Data.Vector.Algorithms.TriHeap as TH +import qualified Data.Vector.Algorithms.Optimal as O + +import qualified Data.Vector.Algorithms.Search as SR + +type Algo e r = forall s mv. MVector mv e => mv s e -> ST s r +type SizeAlgo e r = forall s mv. MVector mv e => mv s e -> Int -> ST s r +type BoundAlgo e r = forall s mv. MVector mv e => mv s e -> Int -> Int -> ST s r + +args = stdArgs + { maxSuccess = 300 + , maxDiscard = 200 + } + +check_Int_sort = forM_ algos $ \(name,algo) -> + quickCheckWith args (label name . prop_fullsort algo) + where + algos :: [(String, Algo Int ())] + algos = [ ("introsort", INT.sort) + , ("insertion sort", INS.sort) + , ("merge sort", M.sort) + , ("tri-heapsort", TH.sort) + ] + +check_Int_partialsort = forM_ algos $ \(name,algo) -> + quickCheckWith args (label name . prop_partialsort algo) + where + algos :: [(String, SizeAlgo Int ())] + algos = [ ("intro-partialsort", INT.partialSort) + , ("tri-heap partialsort", TH.partialSort) + ] + +check_Int_select = forM_ algos $ \(name,algo) -> + quickCheckWith args (label name . prop_select algo) + where + algos :: [(String, SizeAlgo Int ())] + algos = [ ("intro-select", INT.select) + , ("tri-heap select", TH.select) + ] + +check_radix_sorts = do + qc (label "Word8" . prop_fullsort (R.sort :: Algo Word8 ())) + qc (label "Word16" . prop_fullsort (R.sort :: Algo Word16 ())) + qc (label "Word32" . prop_fullsort (R.sort :: Algo Word32 ())) + qc (label "Word64" . prop_fullsort (R.sort :: Algo Word64 ())) + qc (label "Word" . prop_fullsort (R.sort :: Algo Word ())) + qc (label "Int8" . prop_fullsort (R.sort :: Algo Int8 ())) + qc (label "Int16" . prop_fullsort (R.sort :: Algo Int16 ())) + qc (label "Int32" . prop_fullsort (R.sort :: Algo Int32 ())) + qc (label "Int64" . prop_fullsort (R.sort :: Algo Int64 ())) + qc (label "Int" . prop_fullsort (R.sort :: Algo Int ())) + qc (label "(Int, Int)" . prop_fullsort (R.sort :: Algo (Int, Int) ())) + where + qc algo = quickCheckWith args algo + +{- +check_schwartzian = do + quickCheckWith args (prop_schwartzian i2w INS.sortBy) + where + i2w :: Int -> Word + i2w = fromIntegral +-} + +check_stable = do quickCheckWith args (label "merge sort" . prop_stable M.sortBy) + quickCheckWith args (label "radix sort" . prop_stable_radix R.sortBy) + +check_optimal = do qc . label "size 2" $ prop_optimal 2 O.sort2ByOffset + qc . label "size 3" $ prop_optimal 3 O.sort3ByOffset + qc . label "size 4" $ prop_optimal 4 O.sort4ByOffset + where + qc = quickCheck + +check_permutation = do + qc $ label "introsort" . prop_permutation (INT.sort :: Algo Int ()) + qc $ label "intropartial" . prop_sized (const . prop_permutation) + (INT.partialSort :: SizeAlgo Int ()) + qc $ label "introselect" . prop_sized (const . prop_permutation) + (INT.select :: SizeAlgo Int ()) + qc $ label "heapsort" . prop_permutation (TH.sort :: Algo Int ()) + qc $ label "heappartial" . prop_sized (const . prop_permutation) + (TH.partialSort :: SizeAlgo Int ()) + qc $ label "heapselect" . prop_sized (const . prop_permutation) + (TH.select :: SizeAlgo Int ()) + qc $ label "mergesort" . prop_permutation (M.sort :: Algo Int ()) + qc $ label "radix I8" . prop_permutation (R.sort :: Algo Int8 ()) + qc $ label "radix I16" . prop_permutation (R.sort :: Algo Int16 ()) + qc $ label "radix I32" . prop_permutation (R.sort :: Algo Int32 ()) + qc $ label "radix I64" . prop_permutation (R.sort :: Algo Int64 ()) + qc $ label "radix Int" . prop_permutation (R.sort :: Algo Int ()) + qc $ label "radix W8" . prop_permutation (R.sort :: Algo Word8 ()) + qc $ label "radix W16" . prop_permutation (R.sort :: Algo Word16 ()) + qc $ label "radix W32" . prop_permutation (R.sort :: Algo Word32 ()) + qc $ label "radix W64" . prop_permutation (R.sort :: Algo Word64 ()) + qc $ label "radix Word" . prop_permutation (R.sort :: Algo Word ()) + where + qc prop = quickCheckWith args prop + +type BoundSAlgo e r = forall s mv. MVector mv e => mv s e -> e -> Int -> Int -> ST s r + +check_search_range = do + qc $ (label "binarySearchL" .) + . prop_search_inrange (SR.binarySearchLByBounds compare :: BoundSAlgo Int Int) + qc $ (label "binarySearch" .) + . prop_search_inrange (SR.binarySearchByBounds compare :: BoundSAlgo Int Int) + where + qc prop = quickCheckWith args prop + +main = do putStrLn "Int tests:" + check_Int_sort + check_Int_partialsort + check_Int_select + putStrLn "Radix sort tests:" + check_radix_sorts +-- putStrLn "Schwartzian transform (Int -> Word):" +-- check_schwartzian + putStrLn "Stability:" + check_stable + putStrLn "Optimals:" + check_optimal + putStrLn "Permutation:" + check_permutation + putStrLn "Search in range:" + check_search_range diff --git a/tests/Util.hs b/tests/Util.hs new file mode 100644 --- /dev/null +++ b/tests/Util.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE TypeOperators #-} + +module Util where + +import Control.Monad +import Control.Monad.ST + +import Data.Word +import Data.Int + +import qualified Data.Vector as V + +import Data.Vector.Mutable hiding (length) + +import Test.QuickCheck + + +mfromList :: [e] -> ST s (MVector s e) +mfromList l = do v <- new (length l) + fill l 0 v + where + fill [] _ v = return v + fill (x:xs) i v = do write v i x + fill xs (i+1) v + +instance (Arbitrary e) => Arbitrary (V.Vector e) where + arbitrary = fmap V.fromList arbitrary + +instance Arbitrary Int8 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Int16 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Int32 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Int64 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Word8 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Word16 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Word32 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Word64 where + arbitrary = fromInteger `fmap` arbitrary + +instance Arbitrary Word where + arbitrary = fromInteger `fmap` arbitrary diff --git a/vector-algorithms.cabal b/vector-algorithms.cabal new file mode 100644 --- /dev/null +++ b/vector-algorithms.cabal @@ -0,0 +1,71 @@ +Name: vector-algorithms +Version: 0.3 +License: BSD3 +License-File: LICENSE +Author: Dan Doel +Maintainer: Dan Doel +Homepage: http://code.haskell.org/~dolio/ +Category: Data +Synopsis: Efficient algorithms for vector arrays +Description: Efficient algorithms for vector arrays + be sure to compile with -O2, and -fvia-C -optc-O3 is + recommended. +Build-Type: Simple +Cabal-Version: >= 1.2.3 + +Flag BoundsChecks + Description: Enable bounds checking + Default: True + +Flag UnsafeChecks + Description: Enable bounds checking in unsafe operations at the cost of a + significant performance penalty. + Default: False + +Flag InternalChecks + Description: Enable internal consistency checks at the cost of a + significant performance penalty. + Default: False + +Library + Build-Depends: base >= 3 && < 5, vector >= 0.5 && < 0.6, primitive >=0.2 && <0.3 + + Exposed-Modules: + Data.Vector.Algorithms.Combinators + Data.Vector.Algorithms.Optimal + Data.Vector.Algorithms.Insertion + Data.Vector.Algorithms.Intro + Data.Vector.Algorithms.Merge + Data.Vector.Algorithms.Radix + Data.Vector.Algorithms.Search + Data.Vector.Algorithms.TriHeap + + Other-Modules: + Data.Vector.Algorithms.Common + + Extensions: + BangPatterns, + TypeOperators, + Rank2Types, + ScopedTypeVariables, + CPP + + GHC-Options: + -Odph + -funbox-strict-fields + + Include-Dirs: + include + + Install-Includes: + vector.h + + if flag(BoundsChecks) + cpp-options: -DVECTOR_BOUNDS_CHECKS + + if flag(UnsafeChecks) + cpp-options: -DVECTOR_UNSAFE_CHECKS + + if flag(InternalChecks) + cpp-options: -DVECTOR_INTERNAL_CHECKS +