diff --git a/Math/Combinatorics/Multiset.hs b/Math/Combinatorics/Multiset.hs --- a/Math/Combinatorics/Multiset.hs +++ b/Math/Combinatorics/Multiset.hs @@ -1,3 +1,6 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} -- | Efficient combinatorial algorithms over multisets, including -- generating all permutations, partitions, subsets, cycles, and @@ -43,9 +46,11 @@ , splits , kSubsets - -- * Cycles + -- * Cycles and bracelets , cycles + , bracelets + , genFixedBracelets -- * Miscellaneous @@ -53,9 +58,12 @@ ) where -import Data.List (group, sort) -import Control.Arrow (first, second, (&&&), (***)) -import Data.Maybe (catMaybes) +import Control.Arrow (first, second, (&&&), (***)) +import Control.Monad (forM_, when) +import Control.Monad.Trans.Writer +import qualified Data.IntMap.Strict as IM +import Data.List (group, partition, sort) +import Data.Maybe (catMaybes, fromJust) type Count = Int @@ -63,7 +71,7 @@ -- We maintain the invariants that the counts are always positive, -- and no element ever appears more than once. newtype Multiset a = MS { toCounts :: [(a, Count)] } - deriving (Show) + deriving (Show, Functor) -- | Construct a 'Multiset' from a list of (element, count) pairs. -- Precondition: the counts must all be positive, and there must not @@ -74,7 +82,7 @@ -- | Extract just the element counts from a multiset, forgetting the -- elements. getCounts :: Multiset a -> [Count] -getCounts (MS xs) = map snd xs +getCounts = map snd . toCounts -- | Compute the total size of a multiset. size :: Multiset a -> Int @@ -103,9 +111,6 @@ (+:) :: (a, Count) -> Multiset a -> Multiset a (+:) = consMS -instance Functor Multiset where - fmap f = fromCounts . (map . first $ f) . toCounts - -- | Convert a multiset to a list. toList :: Multiset a -> [a] toList = expandCounts . toCounts @@ -125,8 +130,9 @@ fromListEq :: Eq a => [a] -> Multiset a fromListEq = fromCounts . fromListEq' where fromListEq' [] = [] - fromListEq' (x:xs) = (x, 1 + count x xs) : fromListEq' (filter (/=x) xs) - count x = length . filter (==x) + fromListEq' (x:xs) = (x, 1 + length xEqs) : fromListEq' xNeqs + where + (xEqs, xNeqs) = partition (==x) xs -- | Make a multiset with one copy of each element from a list of -- distinct elements. @@ -389,6 +395,10 @@ addElt _ 0 = id addElt x k = ((x,k) +:) +---------------------------------------------------------------------- +-- Cycles (aka Necklaces) +---------------------------------------------------------------------- + -- | Generate all distinct cycles, aka necklaces, with elements taken -- from a multiset. See J. Sawada, \"A fast algorithm to generate -- necklaces with fixed content\", J. Theor. Comput. Sci. 301 (2003) @@ -418,22 +428,194 @@ -- necklace. @pre@ is the current (reversed) prefix of the -- necklaces being generated. cycles' :: Int -> Int -> Int -> [(Int, a)] -> [(Int, (a,Count))] -> [[a]] -cycles' n t p pre [] | n `mod` p == 0 = [map snd pre] +cycles' n _ p pre [] | n `mod` p == 0 = [map snd pre] | otherwise = [] cycles' n t p pre xs = - (takeWhile ((>=atp) . fst) xs) >>= \(j, (xj,nj)) -> + (takeWhile ((>=atp) . fst) xs) >>= \(j, (xj,_)) -> cycles' n (t+1) (if j == atp then p else t) ((j,xj):pre) (remove j xs) where atp = fst $ pre !! (p - 1) -remove j [] = [] +remove :: Int -> [(Int, (a, Int))] -> [(Int, (a, Int))] +remove _ [] = [] remove j (x@(j',(xj,nj)):xs) | j == j' && nj == 1 = xs | j == j' = (j',(xj,nj-1)):xs | otherwise = x:remove j xs +---------------------------------------------------------------------- +-- Bracelets +---------------------------------------------------------------------- + +-- Some utilities + +-------------------------------------------------- +-- Indexable and Snocable classes + +class Snocable p a where + (|>) :: p -> a -> p + +-- 1-based indexing +class Indexable p where + (!) :: p -> Int -> Int + +-------------------------------------------------- +-- Prenecklaces + +type PreNecklace = [Int] + +-- A prenecklace, stored backwards, along with its length and its +-- first element cached for quick retrieval. +data Pre = Pre !Int (Maybe Int) PreNecklace + deriving (Show) + +emptyPre :: Pre +emptyPre = Pre 0 Nothing [] + +getPre :: Pre -> PreNecklace +getPre (Pre _ _ as) = reverse as + +instance Snocable Pre Int where + (Pre 0 _ []) |> a = Pre 1 (Just a) [a] + (Pre t a1 as) |> a = Pre (t+1) a1 (a:as) + +instance Indexable Pre where + _ ! 0 = 0 + (Pre _ (Just a1) _) ! 1 = a1 + (Pre t _ as) ! i = as !! (t-i) + -- as stores a_t .. a_1. + -- a_1 is the last element, i.e. with index t-1. + -- a_2 has index t-2. + -- In general, a_i has index t-i. + +-------------------------------------------------- +-- Run-length encoding + +-- Run-length encodings. Stored in *reverse* order for easy access to +-- the end. +data RLE a = RLE !Int !Int [(a,Int)] + deriving (Show) + -- First Int is the total length of the decoded list. + -- Second Int is the number of blocks. + +emptyRLE :: RLE a +emptyRLE = RLE 0 0 [] + +compareRLE :: Ord a => [(a,Int)] -> [(a,Int)] -> Ordering +compareRLE [] [] = EQ +compareRLE [] _ = LT +compareRLE _ [] = GT +compareRLE ((a1,n1):rle1) ((a2,n2):rle2) + | (a1,n1) == (a2,n2) = compareRLE rle1 rle2 + | a1 < a2 = LT + | a1 > a2 = GT + | (n1 < n2 && (null rle1 || fst (head rle1) < a2)) || (n1 > n2 && not (null rle2) && a1 < fst (head rle2)) = LT + | otherwise = GT + +instance Indexable (RLE Int) where + (RLE _ _ []) ! _ = error "Bad index in (!) for RLE" + (RLE n b ((a,v):rest)) ! i + | i <= v = a + | otherwise = (RLE (n-v) (b-1) rest) ! (i-v) + +instance Eq a => Snocable (RLE a) a where + (RLE _ _ []) |> a' = RLE 1 1 [(a',1)] + (RLE n b rle@((a,v):rest)) |> a' + | a == a' = RLE (n+1) b ((a,v+1):rest) + | otherwise = RLE (n+1) (b+1) ((a',1):rle) + +-------------------------------------------------- +-- Prenecklaces + RLE + +-- Prenecklaces along with a run-length encoding. +data Pre' = Pre' Pre (RLE Int) + deriving Show + +emptyPre' :: Pre' +emptyPre' = Pre' emptyPre emptyRLE + +getPre' :: Pre' -> PreNecklace +getPre' (Pre' pre _) = getPre pre + +instance Indexable Pre' where + _ ! 0 = 0 + (Pre' (Pre len _ _) rle) ! i = rle ! (len - i + 1) + +instance Snocable Pre' Int where + (Pre' p rle) |> a = Pre' (p |> a) (rle |> a) + +-------------------------------------------------- +-- Bracelet generation + +type Bracelet = [Int] + +-- | An optimized bracelet generation algorithm, based on +-- S. Karim et al, "Generating Bracelets with Fixed Content". +-- +-- +-- @genFixedBracelets n content@ produces all bracelets (unique up +-- to rotation and reflection) of length @n@ using @content@, which +-- consists of a list of pairs where the pair (a,i) indicates that +-- element a may be used up to i times. It is assumed that the elements +-- are drawn from [0..k]. +genFixedBracelets :: Int -> [(Int,Int)] -> [Bracelet] +genFixedBracelets n [(0,k)] | k >= n = [replicate k 0] + | otherwise = [] +genFixedBracelets n content = execWriter (go 1 1 0 (IM.fromList content) emptyPre') + where + go :: Int -> Int -> Int -> IM.IntMap Int -> Pre' -> Writer [Bracelet] () + go _ _ _ con _ | IM.keys con == [0] = return () + go t p r con pre@(Pre' (Pre _ _ as) _) + | t > n = + when (take (n - r) as >= reverse (take (n-r) as) && n `mod` p == 0) $ + tell [getPre' pre] + | otherwise = do + let a' = pre ! (t-p) + forM_ (dropWhile (< a') $ IM.keys con) $ \j -> do + let con' = decrease j con + pre' = pre |> j + c = checkRev2 t pre' + p' | j /= a' = t + | otherwise = p + when (c == EQ) $ go (t+1) p' t con' pre' + when (c == GT) $ go (t+1) p' r con' pre' + + decrease :: Int -> IM.IntMap Int -> IM.IntMap Int + decrease j con + | IM.null con = con + | otherwise = IM.alter q j con + where + q (Just 1) = Nothing + q (Just cnt) = Just (cnt-1) + q _ = Nothing + + checkRev2 _ (Pre' _ (RLE _ _ rle)) = compareRLE rle (reverse rle) + +-- | Generate all distinct bracelets (lists considered equivalent up +-- to rotation and reversal) from a given multiset. The generated +-- bracelets are in lexicographic order, and each is +-- lexicographically smallest among its rotations and reversals. +-- See @genFixedBracelets@ for a slightly more general routine with +-- references. +-- +-- For example, @bracelets $ fromList \"RRRRRRRLLL\"@ yields +-- +-- > ["LLLRRRRRRR","LLRLRRRRRR","LLRRLRRRRR","LLRRRLRRRR" +-- > ,"LRLRLRRRRR","LRLRRLRRRR","LRLRRRLRRR","LRRLRRLRRR"] +bracelets :: Multiset a -> [[a]] +bracelets ms@(MS cnts) = bs + where + contentMap = IM.fromList (zip [0..] (map fst cnts)) + content = zipWith (\i (_,n) -> (i,n)) [0..] cnts + rawBs = genFixedBracelets (size ms) content + bs = map (map (fromJust . flip IM.lookup contentMap)) rawBs + +---------------------------------------------------------------------- +-- sequenceMS +---------------------------------------------------------------------- + -- | Take a multiset of lists, and select one element from each list -- in every possible combination to form a list of multisets. We -- assume that all the list elements are distinct. @@ -445,4 +627,3 @@ uncollate :: ([a], Count) -> [(a, Count)] uncollate (xs, n) = map (\x -> (x,n)) xs - diff --git a/multiset-comb.cabal b/multiset-comb.cabal --- a/multiset-comb.cabal +++ b/multiset-comb.cabal @@ -1,20 +1,19 @@ Name: multiset-comb -Version: 0.2.3 +Version: 0.2.4 Synopsis: Combinatorial algorithms over multisets Description: Various combinatorial algorithms over multisets, including generating all permutations, partitions, size-2 partitions, size-k subsets, - and Sawada's algorithm for generating all - necklaces with elements from a multiset. + necklaces, and bracelets. License: BSD3 License-file: LICENSE Author: Brent Yorgey -Maintainer: byorgey@cis.upenn.edu +Maintainer: byorgey@gmail.com bug-reports: http://hub.darcs.net/byorgey/multiset-comb/issues Copyright: (c) 2010 Brent Yorgey Stability: Experimental Category: Math -Tested-with: GHC == 7.4.2, GHC == 7.6.1 +Tested-with: GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1 Build-type: Simple Cabal-version: >=1.10 source-repository head @@ -23,5 +22,7 @@ Library Exposed-modules: Math.Combinatorics.Multiset - Build-depends: base >= 3 && < 5 + Build-depends: base >= 3 && < 5, + containers >= 0.4 && < 0.6, + transformers >= 0.3 && < 0.5 Default-language: Haskell2010