diff --git a/Data/Bits/Ordered.hs b/Data/Bits/Ordered.hs new file mode 100644 --- /dev/null +++ b/Data/Bits/Ordered.hs @@ -0,0 +1,302 @@ + +-- | Efficiently enumerate the bits in data types in order of population +-- count. This yields, say, @000, 001, 010, 100, 011, 101, 110, 111@ (or +-- @0, 1, 2, 4, 3, 5, 6, 7@). Another view is of looking at the bits as +-- a bitset, first enumerating the empty set, then all 1-element sets, all +-- 2-element sets, up to the set size. +-- +-- The enumerator can be inlined with @unfoldr@ (of the @vector@ package) +-- and is a good producer. +-- +-- A memo-table is available, since @popCntSorted@ is still waiting for an +-- efficient @popCntEnumerated@ that does not require sorting. + +module Data.Bits.Ordered + -- bitset operations + ( lsbZ + , nextActiveZ + , maybeNextActive + , maybeLsb + -- population operations + , popPermutation + , popComplement + -- stream ever larger population counts + , popCntSorted + , popCntMemoInt + , popCntMemoWord + , popShiftL + , popShiftR + -- structures with active bits + , activeBitsL + , activeBitsS + , activeBitsV + -- temporary + ) where + +import Control.Arrow +import Data.Bits +import Data.Bits.Extras +import Data.Ord (comparing) +import Data.Vector.Unboxed (Unbox) +import Data.Word(Word(..)) +import Debug.Trace +import qualified Data.Vector.Algorithms.Intro as AI +import qualified Data.Vector.Fusion.Stream as S +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Unboxed as VU + + + +-- * Move from one active bit to the next one + +-- | Capture the no-bit-set case + +captureNull :: Ranked t => t -> (t -> Int) -> Int +captureNull t f = if t==0 then -1 else f t +{-# Inline captureNull #-} + +-- | Get the lowest active bit. Returns @-1@ if no bit is set. + +lsbZ :: Ranked t => t -> Int +lsbZ t = captureNull t lsb +{-# Inline lsbZ #-} + +-- | Given the currently active bit @k@ and the set @t@, get the next +-- active bit. Return @-1@ if there is no next active bit. + +nextActiveZ :: Ranked t => Int -> t -> Int +nextActiveZ k t = lsbZ $ (t `shiftR` (k+1)) `shiftL` (k+1) +{-# Inline nextActiveZ #-} + +-- | Return next active bit, using @Maybe@. + +maybeNextActive :: Ranked t => Int -> t -> Maybe Int +maybeNextActive k t = if t'==0 then Nothing else Just (lsb t') + where t' = (t `shiftR` (k+1) `shiftL` (k+1)) +{-# Inline maybeNextActive #-} + +-- | @Maybe@ the lowest active bit. + +maybeLsb :: Ranked t => t -> Maybe Int +maybeLsb t = if t==0 then Nothing else Just (lsb t) +{-# Inline maybeLsb #-} + +-- | List of all active bits, from lowest to highest. + +activeBitsL :: Ranked t => t -> [Int] +activeBitsL = S.toList . activeBitsS +{-# Inline activeBitsL #-} + +-- | A generic vector (specializes to the corrept type) of the active bits, +-- lowest to highest. + +activeBitsV :: (Ranked t, VG.Vector v Int) => t -> v Int +activeBitsV = VG.unstream . activeBitsS +{-# Inline activeBitsV #-} + +-- | A stream with the currently active bits, lowest to highest. + +activeBitsS :: (Ranked t, Monad m) => t -> SM.Stream m Int +activeBitsS t = SM.unfoldr (fmap (id &&& (`maybeNextActive` t))) (maybeLsb t) +{-# Inline activeBitsS #-} + +-- * Population count methods + +-- | The /slow/ default implementation. We sort the vector, not the list, +-- as sorting will walk the whole data structure anyway, and the vector +-- requires not as much memory. +-- +-- Replaced @popCount &&& id@ as sort, which provides for @a Int -> VU.Vector n +popCntSorted n = VU.modify (AI.sortBy (comparing (popCount &&& activeBitsL))) $ VU.enumFromN 0 (2^n) +{-# Inline popCntSorted #-} + +-- | Memoized version of 'popCntSorted' for @Int@s. +-- +-- NOTE Since this uses @popCntSorted@ for now it will still require a lot +-- of memory for sorting the vector! + +popCntMemoInt + :: Int -- ^ size of the set we want. If larger than memo limit, will just call 'popCntSorted' + -> VU.Vector Int +popCntMemoInt n + | n>limit = error $ "for safety reasons, memoization is only performed for popcounts up to " ++ show limit ++ " bits, memoize manually!" + | otherwise = _popCntMemoInt !! n + where limit = 28 +{-# Inline popCntMemoInt #-} + +-- | Memoizes popcount arrays. The limit to memoization is enforced by +-- popCntMemo, not here. + +_popCntMemoInt = map popCntSorted [0..] +{-# NoInline _popCntMemoInt #-} + +-- | Memoized version of 'popCntSorted' for @Word@s. +-- +-- NOTE Since this uses @popCntSorted@ for now it will still require a lot +-- of memory for sorting the vector! + +popCntMemoWord + :: Int -- ^ size of the set we want. If larger than memo limit, will just call 'popCntSorted' + -> VU.Vector Word +popCntMemoWord n + | n>limit = error $ "for safety reasons, memoization is only performed for popcounts up to " ++ show limit ++ " bits, memoize manually!" + | otherwise = _popCntMemoWord !! n + where limit = 28 +{-# Inline popCntMemoWord #-} + +-- | Memoizes popcount arrays. The limit to memoization is enforced by +-- popCntMemo, not here. + +_popCntMemoWord = map popCntSorted [0..] +{-# NoInline _popCntMemoWord #-} + +-- | Enumerate all sets with the same population count. Given a population +-- @i@, this returns @Just j@ with @j>i@ (but same number of set bits) or +-- @Nothing@. For a population count of @k@, start with @2^(k+1) -1@. +-- +-- Note that @sort permutations == sort (nub permutations)@ if +-- @permutations@ is a set of all permutations for a given @popCount@ +-- generated by @popPermutation@. The @Data.List.permutations@ functions +-- will create duplicates. +-- +-- cf +-- + +popPermutation + :: Ranked t + => Int -- size of the set we want. (i.e. numbor of bits available for @0@ or @1@) + -> t -- current population + -> Maybe t -- Just the new population, or nothing if now higher-ordered population exists. +popPermutation !h' !s' + | popCount s' < 1 || h' < 2 = Nothing + | Just k <- findK (h' -2) + , Just l <- findL k (h' -1) + = let swp = setBit (clearBit s' k) l + in Just $ reverseFrom (k+1) (h' -1) swp swp + | otherwise = Nothing + where findK k + | k < 0 = Nothing + | testBit s' k && not (testBit s' (k+1)) = Just k + | otherwise = findK (k-1) + findL k l + | l <= k = Nothing + | not $ testBit s' l = Just l + | otherwise = findL k $ l-1 + reverseFrom u d src tgt + | u >= h' = tgt + | otherwise = reverseFrom (u+1) (d-1) src (assignBit (assignBit tgt u (testBit src d)) d (testBit src u)) +{-# Inline popPermutation #-} + +-- | Given a population, get the complement. The first argument is the size +-- of the population (say. 8 for 8 bits); the second the current +-- population. +-- +-- Examples: +-- +-- >>> popComplement 5 (3 :: Int) +-- 28 +-- +-- >>> popComplement 6 (3 :: Int) +-- 60 + +popComplement + :: Ranked t + => Int -- size of the population set + -> t -- current population + -> t -- complement of the population. All bits higher than the highest bit are kept zero. +popComplement !h !s = mask .&. complement s + where mask = (2^h -1) +{-# Inline popComplement #-} + +-- | Move a population more to the left. This, effectively, introduces @0@s +-- in the set, whereever the @mask@ has a @0@. Only as many @1@s can be +-- set, as the mask holds. Assume that you have a bitmask @mask = 10101@ +-- and a least-significant aligned population @11@, then given mask and +-- population you'd like to see @00101@, i.e. the two lowest one bits of +-- the mask are set. @101@ would set the lowest and third one bit. +-- +-- Examples: +-- +-- >>> popShiftL (21::Int) 3 -- 10101 00011 -- 00101 +-- 5 +-- >>> popShiftL (28::Int) 0 -- 11100 00000 -- 00000 +-- 0 +-- >>> popShiftL (28::Int) 1 -- 11100 00001 -- 00100 +-- 4 +-- >>> popShiftL (28::Int) 2 -- 11100 00010 -- 01000 +-- 8 +-- >>> popShiftL (28::Int) 3 -- 11100 00011 -- 01100 +-- 12 + +popShiftL + :: (Ranked t) + => t -- the mask + -> t -- the population + -> t -- final population +popShiftL mask lsp = go 0 0 mask lsp where + go !acc !(k::Int) !m !l + | l==0 || m==0 = acc + | testBit m 0 + , testBit l 0 = go (acc + bit k) (k+1) (unsafeShiftR m 1) (unsafeShiftR l 1) + | not $ testBit m 0 = go acc (k+1) (unsafeShiftR m 1) l + | not $ testBit l 0 = go acc (k+1) (unsafeShiftR m 1) (unsafeShiftR l 1) +{-# Inline popShiftL #-} + +-- | Effectively compresses a bitset, given a mask. Removes set elements, +-- whenever the mask is @0@, by moving all remaining elements one to the +-- right. + +popShiftR + :: (Ranked t) + => t -- the mask + -> t -- the population + -> t -- final population +popShiftR mask lsp = go 0 0 mask lsp where + go !acc !k !m !l + | m==0 || l==0 = acc + | testBit m 0 + , testBit l 0 = go (acc .|. bit k) (k+1) (m `unsafeShiftR` 1) (l `unsafeShiftR` 1) + | testBit m 0 = go acc (k+1) (m `unsafeShiftR` 1) (l `unsafeShiftR` 1) + | otherwise = go acc k (m `unsafeShiftR` 1) (l `unsafeShiftR` 1) +{-# Inline popShiftR #-} + + + +-- WARNING: Conditional compilation based on architecture! + +instance Ranked Int where +#if x86_64_HOST_ARCH + lsb = lsb . w64 + rank = rank . w64 + nlz = nlz . w64 +#endif +#if i386_HOST_ARCH + lsb = lsb . w32 + rank = rank . w32 + nlz = nlz . w32 +#endif + {-# Inline lsb #-} + {-# Inline rank #-} + {-# Inline nlz #-} + +instance Ranked Word where +#if x86_64_HOST_ARCH + lsb = lsb . w64 + rank = rank . w64 + nlz = nlz . w64 +#endif +#if i386_HOST_ARCH + lsb = lsb . w32 + rank = rank . w32 + nlz = nlz . w32 +#endif + {-# Inline lsb #-} + {-# Inline rank #-} + {-# Inline nlz #-} + diff --git a/Data/Bits/Ordered/QuickCheck.hs b/Data/Bits/Ordered/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/Data/Bits/Ordered/QuickCheck.hs @@ -0,0 +1,74 @@ + +{-# Options_GHC -O0 #-} + +-- | Check a number of properties for popcount-ordered elements. +-- +-- $setup +-- +-- >>> :set -XScopedTypeVariables +-- + +module Data.Bits.Ordered.QuickCheck where + +import Test.QuickCheck hiding ((.&.)) +import Data.Int (Int16(..)) +import Data.Bits +import qualified Data.Vector.Unboxed as VU +import Data.List (groupBy,sort,permutations,nub) +import Data.Function (on) +import Data.Maybe (isJust) +import Control.Monad (join) +import Debug.Trace +import Data.Word (Word) + +import Data.Bits.Ordered + + + +-- | Check if both the memoized version and the population enumeration +-- produce the same multisets, but maybe in different order. +-- +-- prop> \(n :: Int16) -> let b = popCount n in memoSorted b == enumSorted b +-- + +prop_PopCountSet (NonZero (n :: Int16)) = memo == enum + where b = popCount n + memo = memoSorted b + enum = enumSorted b + +memoSorted, enumSorted :: Int -> [[Int]] + +memoSorted b = map sort . groupBy ((==) `on` popCount) $ VU.toList $ popCntMemoInt b +enumSorted b = map sort $ [0] : [ roll (popPermutation b) (Just $ 2^k-1) | k <- [1..b] ] + where roll f (Just k) = k : roll f (f k) + roll _ Nothing = [] + +prop_lsb_Int (x :: Int) = lsbZ x == maybe (-1) id (maybeLsb x) + +prop_lsb_Word (x :: Word) = lsbZ x == maybe (-1) id (maybeLsb x) + +prop_OneBits_Int (x :: Int) = popCount x == length abl && and [ testBit x k | k <- abl ] + where abl = activeBitsL x + +-- Tests if we actually generate all permutations. + +prop_allPermutations (a :: Int , b :: Int) = and $ zipWith cmp (sort qs) (sort $ nub ps) + where nbs = min a' b' -- number of 1 bits in set + sts = max a' b' -- set size + a' = a `mod` 8 -- finiteBitSize a + b' = b `mod` 8 -- finiteBitSize b + ps = permutations $ replicate (sts - nbs) False ++ replicate nbs True + qs = go (Just $ 2 ^ nbs - 1) + go :: Maybe Int -> [Int] + go Nothing = [] + go (Just k) = k : go (popPermutation sts k) + cmp k as = and [ if a then testBit k c else (not $ testBit k c) | (a,c) <- zip (reverse as) [0 .. ] ] + +-- TODO popComplement + +prop_popShiftL_popShiftR (a::Word,b::Word) = s == l + where m = a .|. b + s = a .&. b + l = popShiftL m r + r = popShiftR m s + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Christian Hoener zu Siederdissen 2011-2012 + +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 Christian Hoener zu Siederdissen nor the names of other + 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/OrderedBits.cabal b/OrderedBits.cabal new file mode 100644 --- /dev/null +++ b/OrderedBits.cabal @@ -0,0 +1,109 @@ +name: OrderedBits +version: 0.0.0.1 +author: Christian Hoener zu Siederdissen +copyright: Christian Hoener zu Siederdissen, 2014 - 2015 +homepage: http://www.bioinf.uni-leipzig.de/~choener/ +maintainer: choener@tbi.univie.ac.at +category: Data +license: BSD3 +license-file: LICENSE +build-type: Simple +stability: experimental +cabal-version: >= 1.10.0 +tested-with: GHC == 7.8.4, GHC == 7.10.1 +synopsis: Efficient ordered (by popcount) enumeration of bits +description: + This library provides efficient methods to enumerate all + elements of a set in order of the population count. First, the + empty set, then all 1-element sets, all 2-element sets, etc. + Such enumerations are important for algorithms over unordered + data sets. Examples include the travelling salesman problem and + the closely related Hamiltonian path problem. + + + +Extra-Source-Files: + README.md + changelog.md + + + +flag llvm + description: build using LLVM + default: False + manual: True + + + +library + build-depends: base >= 4.7 && < 4.9 + , bits == 0.4.* + , primitive >= 0.5 && < 0.7 + , QuickCheck >= 2.7 && < 2.9 + , vector == 0.10.* + , vector-algorithms == 0.6.* + default-language: + Haskell2010 + default-extensions: BangPatterns + , CPP + , FlexibleContexts + , PatternGuards + , ScopedTypeVariables + + exposed-modules: + Data.Bits.Ordered + Data.Bits.Ordered.QuickCheck + ghc-options: + -O2 -funbox-strict-fields + + + +benchmark BenchmarkOrderedBits + build-depends: base + , criterion >= 1.0.2 && < 1.1.1 + , OrderedBits + , vector + default-language: + Haskell2010 + hs-source-dirs: + bench + main-is: + Benchmark.hs + type: + exitcode-stdio-1.0 + ghc-options: + -O2 + -funbox-strict-fields + if flag(llvm) + ghc-options: + -fllvm + -optlo-O3 -optlo-std-compile-opts + -fllvm-tbaa + + + +test-suite properties + type: + exitcode-stdio-1.0 + main-is: + properties.hs + ghc-options: + -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: + tests + default-language: + Haskell2010 + default-extensions: TemplateHaskell + build-depends: base + , OrderedBits + , QuickCheck + , test-framework >= 0.8 && < 0.9 + , test-framework-quickcheck2 >= 0.3 && < 0.4 + , test-framework-th >= 0.2 && < 0.3 + + + +source-repository head + type: git + location: git://github.com/choener/OrderedBits + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,19 @@ +# OrderedBits + +[![Build Status](https://travis-ci.org/choener/OrderedBits.svg?branch=master)](https://travis-ci.org/choener/OrderedBits) + +The OrderedBits library provides methods to generate unboxed vectors of Ints +(and others) ordered by their population count or Hamming distance to the 0 +set. + +Such an order is important for dynamic programming algorithms for Hamiltonian +path problems and the travelling salesman problem. + + + +#### Contact + +Christian Hoener zu Siederdissen +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/Benchmark.hs b/bench/Benchmark.hs new file mode 100644 --- /dev/null +++ b/bench/Benchmark.hs @@ -0,0 +1,36 @@ + +module Main where + +import Criterion.Main +import qualified Data.Vector.Unboxed as VU +import Data.Word (Word(..)) + +import Data.Bits.Ordered + + + +main = defaultMain + [ bgroup "Int" [ bench "08" $ whnf (VU.sum . popCntSorted :: Int -> Int ) 8 + , bench "16" $ whnf (VU.sum . popCntSorted :: Int -> Int ) 16 + , bench "20" $ whnf (VU.sum . popCntSorted :: Int -> Int ) 20 +-- , bench "24" $ whnf (VU.sum . popCntSorted :: Int -> Int ) 24 + ] + , bgroup "Word" [ bench "08" $ whnf (VU.sum . popCntSorted :: Int -> Word) 8 + , bench "16" $ whnf (VU.sum . popCntSorted :: Int -> Word) 16 + , bench "20" $ whnf (VU.sum . popCntSorted :: Int -> Word) 20 +-- , bench "24" $ whnf (VU.sum . popCntSorted :: Int -> Word) 24 + ] + , bgroup "MemoInt" [ bench "08" $ whnf (VU.sum . popCntMemoInt :: Int -> Int ) 8 + , bench "16" $ whnf (VU.sum . popCntMemoInt :: Int -> Int ) 16 + , bench "20" $ whnf (VU.sum . popCntMemoInt :: Int -> Int ) 20 +-- , bench "24" $ whnf (VU.sum . popCntMemoInt :: Int -> Int ) 24 + ] + , bgroup "MemoWord" [ bench "08" $ whnf (VU.sum . popCntMemoWord :: Int -> Word) 8 + , bench "16" $ whnf (VU.sum . popCntMemoWord :: Int -> Word) 16 + , bench "20" $ whnf (VU.sum . popCntMemoWord :: Int -> Word) 20 +-- , bench "24" $ whnf (VU.sum . popCntMemoWord :: Int -> Word) 24 + ] +-- , bgroup "small ops" [ bench " +-- ] + ] + diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,5 @@ +0.0.0.1 +------- + +- initial checkin of the (naive) sorted implementation +- memoization of Int-sets up to 31 bit diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,20 @@ + +module Main where + +import Test.Framework.Providers.QuickCheck2 +import Test.Framework.TH + +import qualified Data.Bits.Ordered.QuickCheck as QC + + + +prop_PopCountSet = QC.prop_PopCountSet +prop_lsb_Int = QC.prop_lsb_Int +prop_lsb_Word = QC.prop_lsb_Word +prop_OneBits_Int = QC.prop_OneBits_Int +prop_allPermutations = QC.prop_allPermutations +prop_popShiftL_popShiftR = QC.prop_popShiftL_popShiftR + +main :: IO () +main = $(defaultMainGenerator) +