Stefan O'Rear wrote:

Mr. C++ apparently isn't a very good C++ programmer, since his best effort absolutely *pales* in comparison to Julian Seward's BWT:

stefan@stefans:/usr/local/src/hpaste$ head -c 135000 /usr/share/dict/words | (time bzip2 -vvv) > /dev/null (stdin): block 1: crc = 0x25a18961, combined CRC = 0x25a18961, size = 135000 0 work, 135000 block, ratio 0.00 135000 in block, 107256 after MTF & 1-2 coding, 61+2 syms in use initial group 6, [0 .. 0], has 20930 syms (19.5%) initial group 5, [1 .. 1], has 4949 syms ( 4.6%) initial group 4, [2 .. 2], has 20579 syms (19.2%) initial group 3, [3 .. 4], has 17301 syms (16.1%) initial group 2, [5 .. 10], has 24247 syms (22.6%) initial group 1, [11 .. 62], has 19250 syms (17.9%) pass 1: size is 127140, grp uses are 339 550 192 440 12 613 pass 2: size is 51693, grp uses are 321 440 288 316 139 642 pass 3: size is 51358, grp uses are 329 387 376 304 122 628 pass 4: size is 51302, grp uses are 298 421 397 304 125 601 bytes: mapping 21, selectors 433, code lengths 110, codes 51297 final combined CRC = 0x25a18961 2.602:1, 3.075 bits/byte, 61.57% saved, 135000 in, 51887 out.

real 0m0.165s user 0m0.044s sys 0m0.012s

Yup, does slightly more work (huffman coding) in 1/200 the time :)

(Note, on my system .Lazy BWT3 takes 5.3s on the same input)

...OK...so how do I make Haskell go faster still? Presumably by transforming the code into an ugly mess that nobody can read any more...?

Hello Andrew, Saturday, June 23, 2007, 11:21:26 AM, you wrote:

...OK...so how do I make Haskell go faster still?

Presumably by transforming the code into an ugly mess that nobody can read any more...?

bwt transformation is very good researched area, so probably you will not get decent performance (megabytes per second) without lot of work. and of course, no Haskell at all. take look at http://darchiver.narod.ru/dark/Archon3fs.zip -- Best regards, Bulat mailto:Bulat.Ziganshin@gmail.com

Hello Andrew, Saturday, June 23, 2007, 2:45:01 PM, you wrote:

Hey, I'm just glad I managed to get within striking distance of Mr C++. So much for Haskell being "inherently less performant". :-P

my little analysis says that it's probably due to different sort() implementations, so this says nothing about general case -- Best regards, Bulat mailto:Bulat.Ziganshin@gmail.com

On Saturday 23 June 2007 13:02:54 Andrew Coppin wrote:

Well, I altered the code, and it's *still* very short and very readable, and it's just as fast as the (3 pages long) C++ version. >:-D

Indeed. The performance of modern functional programming languages never ceases to amaze me. INRIA typically claim performance within 2x of C for OCaml but there are very few programs where OCaml isn't within 50% now, at least on AMDs. Most of the pedagogical examples of functional programming (e.g. n-queens) are too academic for general consumption but there are plenty of excellent examples. Burrows-Wheeler is a great one. Ray tracing is another. Sudoku solving is fairy good, albeit unusually well suited to arrays. I also like to compare the performance of term-level interpreters or rewriters implemented in different languages. I think it would be interesting to compare Scheme/Mathematica interpreters written in OCaml, SML (MLton) and Haskell, for example. I've benchmarked some interpreters in OCaml and SML (MLton) and MLton's whole-program optimizations are a huge benefit here, making it several times faster than OCaml. I'd like to know how well Haskell would do. -- Dr Jon D Harrop, Flying Frog Consultancy Ltd. The OCaml Journal http://www.ffconsultancy.com/products/ocaml_journal/?e

apfelmus wrote:

Note that the one usually adds an "end of string" character $ in the Burrows-Wheeler transform for compression such that sorting rotated strings becomes sorting suffices.

Yeah, I noticed that the output from by program can never actually be reverted to its original form. ;-) Maybe it I alter the code to stick a 0-byte in there or something... (Hmm, I wonder how Mr C++ managed it? Heh. If I could read C++, maybe I'd know...)

Concerning the algorithm at hand, you can clearly avoid calculating Raw.append over and over:

bwt :: Raw.ByteString -> Raw.ByteString bwt xs = Raw.pack . map (Raw.last) . sort $ rotations where n = length xs rotations = take n . map (take n) . tails $ xs `Raw.append` xs

assuming that take n is O(1).

...which is amusing because that's what the *first* implementation did. ;-) I was trying to avoid O(n^2) RAM usage. :-}

apfelmus wrote:

Andrew Coppin wrote:

...

Yeah, I noticed that the output from by program can never actually be reverted to its original form.

Well it can, but that's a different story told in

Richard S. Bird and Shin-Cheng Mu. Inverting the Burrows-Wheeler transform. http://web.comlab.ox.ac.uk/oucl/work/richard.bird/publications.html #DBLP:journals/jfp/BirdM04

Oh, and you had a function inv_bwt, right?

Implementation #1 had an inverse - but it works by finding an unused character and prepending that to the input before doing the forward BWT. ;-) For the other implementations, there is no inverse at all.

...

I was trying to avoid O(n^2) RAM usage. :-}

Note that for ByteStrings, this takes only O(n) RAM because the substrings are shared. But for lists, this would take O(n^2) RAM since (take n) cannot share hole sublists.

Presumably that's only true for *strict* ByteStrings? (I still don't actually understand how lazy bytestrings are any different to normal lists - or why they produced such a massive performance increase...)

An O(n) choice for lists that doesn't recalculate ++ all the time would be

bwt :: Ord a => [a] -> [a] bwt xs = map last . sortBy (compare `on` (take n)) $ rotations where n = length xs rotations = take n . tails $ xs ++ xs

with the well-known

on :: (a -> a -> c) -> (b -> a) -> (b -> b -> c) on g f x y = g (f x) (f y)

Interesting... But how does it perform? ;-)

I enjoy code like this that requires laziness. My modified version of your code is below... Bertram Felgenhauer wrote:

...

...

...

Code: >>>

"bwt" implements a variation of the Burrows-Wheeler transform, using \0 as a sentinel character for simplicity. The sentinel has to be smaller than all other characters in the string.

...

bwt xs = let suffixes = [(a,as) | a:as <- tails ('\0':xs)] in map fst . sortBy (\(_,a) (_,b) -> a `compare` b) $ suffixes

"rbwt" implements the corresponding inverse BWT. It's a fun knot tying exercise.

...

rbwt xs = let res = sortBy (\(a:as) (b:bs) -> a `compare` b) (zipWith' (:) xs res) in tail . map snd . zip xs $ head res

"zipWith'" is a variant of zipWith that asserts that the third argument has the same shape as the second one.

...

zipWith' f [] _ = [] zipWith' f (x:xs) ~(y:ys) = f x y : zipWith' f xs ys

<<< End Code <<

I did not like the look of (map snd . zip xs) since it looks like a no-op (that constructs a useless (,) which may or may not be elided by a sufficiently smart compiler). But it is using the fact that xs is finite and (head res) is not to do "take (length xs) $ head res" without the extra traversal and math. But one can abuse (zipWith' (flip const)) for a 'rbwt' appeals to me more:

import Data.List

f `on` g = \x y -> (g x) `f` (g y)

zipWith' f [] _ = [] zipWith' f (x:xs) ~(y:ys) = f x y : zipWith' f xs ys

bwt = map head . sortBy (compare `on` tail) . init . tails . ('\0':)

rbwt xs = tail $ zipWith' (flip const) xs $ head res where res = sortBy (compare `on` head) (zipWith' (:) xs res)

While I was removing (,) from the 'rbwt' I went ahead and removed it from 'bwt' as well. This was, of course, pointless... Thanks for the fun example, Chris

Felipe Almeida Lessa wrote:

On 6/23/07, Bertram Felgenhauer wrote:

...

"rbwt" implements the corresponding inverse BWT. It's a fun knot tying exercise.

...

rbwt xs = let res = sortBy (\(a:as) (b:bs) -> a `compare` b) (zipWith' (:) xs res) in tail . map snd . zip xs $ head res

Indeed it's very fun =). Although I must admit that, even after staring at the code for some time, I still don't see how it works. Is it too much to ask for a brief explanation?

I see that zipWith' creates all the lazy list without requiring any knowledge of the second list, which means that the list to be sorted is created nicely. But when sortBy needs to compare, say, the first with the second items, it forces the thunk of the lazy list. But that thunk depends on the sorted list itself!

What am I missing? ;)

Thanks in advance,

I just read and understood it. I'll walk through my version:

import Data.List

f `on` g = \x y -> (g x) `f` (g y)

zipWith' f [] _ = [] zipWith' f (x:xs) ~(y:ys) = f x y : zipWith' f xs ys

bwt = map head . sortBy (compare `on` tail) . init . tails . ('\0':)

rbwt xs = tail $ zipWith' (flip const) xs $ head res where res = sortBy (compare `on` head) (zipWith' (:) xs res)

The key is that sort and sortBy are very very lazy. 'res' is a finite list of infinite strings. Consider sorting such a finite list of infinite lists without the recursive definition:

example = let a = [1,3..] b = [1,2..] c = [2,4..] d = [1,4..] sorted = sort [a,b,c,d] firstFiveOfEach = map (take 5) sorted in mapM_ print firstFiveOfEach

Which in ghci produces:

*Main> example [1,2,3,4,5] [1,3,5,7,9] [1,4,7,10,13] [2,4,6,8,10]

Now for a sidetrack:

example = let a = [1,3..] b = [1,2..] c = [2,4..] d = [1,4..] sorted = sort [a,b,c,d,c] firstFiveOfEach = map (take 5) sorted in mapM_ print firstFiveOfEach

Produces the first three elements before 'c' ....

*Main> example [1,2,3,4,5] [1,3,5,7,9] [1,4,7,10,13]

...and then hangs because "compare c c" never terminates. This lazy sorting is also why "head . sort $ foo" is an O(N) way to get the minimum of a list "foo" in Haskell. So the code for rbwt depends on this lazy sorting behavior. The recursive nature of the definition is the difficult to (directly) invent part. The result of (zipWith' (:) xs ___) is a list the same length as 'xs' and where each item is a list with the head item taken from 'xs'. So it has "primed the pump" of the recursive definition by giving the first Char in each of the strings. Compare with the much simpler 'ones' and 'odds': ones = 1 : ones odds = 1 : map (2+) odds What about the tails of all these lists? They come from res. And res comes from sorting the results of zipWith'. But the sorting is done only by looking at the head element, which is the one that comes from xs. That is the "sortBy (\(a:as) (b:bs) -> a `compare` b)" or "sortBy (compare `on` head)". So the sorting routine does not care about the ___ in (zipWith' (:) xs ___) at all, since it only examines the head which comes from xs.