Re: [Haskell-cafe] None brute-force BWT algorithm
Hi, I think that version is still a brute-force solution. The only difference is that it uses EOF (sentinel) so that it can sort the suffixes instead of rotations. However, the big-O complexity is the same. Let's take the rbwt for example: rbwt xs = let res = sortBy (\(a:as) (b:bs) - a `compare` b) (zipWith' (:) xs res) in tail . map snd . zip xs $ head res Here we can see that, although the infinity res is lazy-evaluated, it actually sorts the matrix N times, adding one column per evaluation. each time there are N elements to be sorted and the length of every element grows proportion to N, so the time is O( N * (N^2) * \lg (N^2) ) = O(N^3 \lg N) While my brute-force program provided in previous post is as same as O(N^3 \lg N). However, if the random access time is O(1) (on array) but not O(N) (on list), my 2nd program is much faster: Here is the modified version with Array. (import Data.Array) ibwt'' :: (Ord a) = ([a], Int) - [a] ibwt'' (r, i) = fst $ iterate f ([], i) !! n where t = listArray (0, n-1) $ map snd $ sort $ zip r [0..] ra = listArray (0, n-1) r f (s, j) = let j' = t ! j in (s++[ra ! j'], j') n = length r This version only sort the input data 1 time, (proportion to O(N * \lg N), after that the transform vector t is generated. Then it iterates N times on t to get the result. so the total time is O(N * \lg N) + O(N) = O(N * \lg N) This should be much better than the brute force one. the idea is that, we can get the result without reconstructing the complete matrix, Only two columns (L F) are enough. But how to turn the random access idea of transform-vector into purely functional settings? here I mean what if Haskell doesn't provide constant access time array? One idea is to to turn the transform vector T into a function, just like what we've done in KMP implementation in FP way. Does such solution exist? Regards. -- Larry, LIU Xinyu https://github.com/liuxinyu95/AlgoXY On Friday, June 24, 2011 6:50:46 AM UTC+8, Henning Thielemann wrote: On Wed, 22 Jun 2011, larry.liuxinyu wrote: Hi, I read a previous thread about BWT implementation in Haskell: http://www.mail-archive.com/haskell-cafe@haskell.org/msg25609.html and http://sambangu.blogspot.com/2007/01/burrows-wheeler-transform-in-haskell They are all in a `brute-force' way, that is implement based on Burrows-Wheeler's definition like below: I thought that the knot-tying solution by Bertram Felgenhauer in the same thread was both elegant and efficient: http://www.mail-archive.com/haskell-cafe%40haskell.org/msg25692.html ___ Haskell-Cafe mailing list haskel...@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] None brute-force BWT algorithm
Hi, I read a previous thread about BWT implementation in Haskell: http://www.mail-archive.com/haskell-cafe@haskell.org/msg25609.html and http://sambangu.blogspot.com/2007/01/burrows-wheeler-transform-in-haskell They are all in a `brute-force' way, that is implement based on Burrows-Wheeler's definition like below: BWT: sort the rotations of input S to get a matrix M', return the last column L, and the row I, where S appears in M' -- O( n^2 \lg n) bwt :: (Ord a)=[a]-([a], Int) bwt s = (map last m, i) where m = sort [ drop i s ++ take i s | i-[1..length s]] (Just i) = elemIndex s m And the IBWT: Re-construct M' by iteratively sort on input L, add one column at a time, and pick the I-th row in M' -- O( n^3 \lg n ) ibwt :: (Ord a)= ([a], Int) - [a] ibwt (r, i) = m !! i where m = iterate f (replicate n []) !! n f = sort . zipWith (:) r n = length r I read Richard Bird's book, `Pearls of functional algorithm design', there is another solution. Although it is deduced step by step, the core idea is based on random access the list by index. The algorithm mentioned in the book uses suffixes for sorting instead of rotations. The performance are same in terms of big-O. I wrote the following program accordingly. BWT: sort S along with the index to get a new order of IDs, and return a permutation of S based on IDs. -- O( n^2 \lg n) if (!!) takes O(n) time bwt' :: (Ord a)= [a] - ([a], Int) bwt' s = (l, i) where l = map (\i-s !! ((i-1) `mod` n)) ids (Just i) = elemIndex 0 ids ids = map snd $ sortBy (compare `on` fst) $ zip rots [0..] rots = init $ zipWith (++) (tails s) (inits s) -- only non-empties n = length s IBWT: Sort the input L along with index to get a Transform vector, T [1], then permute L iteratively on T start from row I. -- O( n^2 ) if (!!) takes O(n) time ibwt' :: (Ord a) = ([a], Int) - [a] ibwt' (r, i) = fst $ iterate f ([], i) !! n where t = map snd $ sortBy (compare `on` fst) $ zip r [0..] f (s, j) = let j' = t !! j in (s++[r !! j'], j') n = length r However, As I commented, the (!!) takes time proportion to the length of the list, Although it can be turned into real Haskell Array by listArray (0, n-1) xs. I wonder if there are any purely functional implementations of BWT/IBWT, which don't base on random access idea nor in brute-force way. [Reference] [1], Mark Nelson. `Data Compression with the Burrows-Wheeler Transform'. http://marknelson.us/1996/09/01/bwt/ -- Larry, LIU Xinyu https://github.com/liuxinyu95/AlgoXY ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] AVL Tree in a pattern matching way
Hi, I browsed the current AVL tree implementation in Hackage http://hackage.haskell.org/packages/archive/AvlTree/4.2/doc/html/src/Data-Tree-AVL-Push.html AVL tree denote the different of height from right sub-tree to left sub-tree as delta, to keep the balance, abs(delta)=1 is kept as invariant. So the typical implementation define N (Negative), P (Positive), and Z (zero) as the tree valid nodes (and the Empty as the trivial case). When a new element is inserted, the program typically first check if the result will break the balance, and process rotation to keep the balance of the tree. Some other pure functional implementation takes the same approach, for example: Guy Cousineau and Michel Mauny. ``The Functional Approach to Programming''. pp 173 ~ 186 Consider the elegant implementation of Red-black tree in pattern matching way by Chris Okasaki, I tried to use the same method in AVL tree, and here is the result. module AVLTree where -- for easy verification, I used Quick Check package. import Test.QuickCheck import qualified Data.List as L -- for verification purpose only -- Definition of AVL tree, it is almost as same as BST, besides a new field to store delta. data AVLTree a = Empty | Br (AVLTree a) a (AVLTree a) Int insert::(Ord a)=AVLTree a - a - AVLTree a insert t x = fst $ ins t where -- result of ins is a pair (t, d), t: tree, d: increment of height ins Empty = (Br Empty x Empty 0, 1) ins (Br l k r d) | x k = node (ins l) k (r, 0) d | x == k= (Br l k r d, 0) -- For duplicate element, we just ignore it. | otherwise = node (l, 0) k (ins r) d -- params: (left, increment on left) key (right, increment on right) node::(AVLTree a, Int) - a - (AVLTree a, Int) - Int - (AVLTree a, Int) node (l, dl) k (r, dr) d = balance (Br l k r d', delta) where d' = d + dr - dl delta = deltaH d d' dl dr -- delta(Height) = max(|R'|, |L'|) - max (|R|, |L|) -- where we denote height(R) as |R| deltaH :: Int - Int - Int - Int - Int deltaH d d' dl dr | d =0 d' =0 = dr | d =0 d' =0 = d+dr | d =0 d' =0 = dl - d | otherwise = dl -- Here is the core pattern matching part, there are 4 cases need rebalance balance :: (AVLTree a, Int) - (AVLTree a, Int) balance (Br (Br (Br a x b dx) y c (-1)) z d (-2), _) = (Br (Br a x b dx) y (Br c z d 0) 0, 0) balance (Br a x (Br b y (Br c z d dz)1)2, _) = (Br (Br a x b 0) y (Br c z d dz) 0, 0) balance (Br (Br a x (Br b y c dy)1) z d (-2), _) = (Br (Br a x b dx') y (Br c z d dz') 0, 0) where dx' = if dy == 1 then -1 else 0 dz' = if dy == -1 then 1 else 0 balance (Br a x (Br (Br b y c dy) z d (-1))2, _) = (Br (Br a x b dx') y (Br c z d dz') 0, 0) where dx' = if dy == 1 then -1 else 0 dz' = if dy == -1 then 1 else 0 balance (t, d) = (t, d) -- Here are some auxiliary functions for verification -- check if a AVLTree is valid isAVL :: (AVLTree a) - Bool isAVL Empty = True isAVL (Br l _ r d) = and [isAVL l, isAVL r, d == (height r - height l), abs d = 1] height :: (AVLTree a) - Int height Empty = 0 height (Br l _ r _) = 1 + max (height l) (height r) checkDelta :: (AVLTree a) - Bool checkDelta Empty = True checkDelta (Br l _ r d) = and [checkDelta l, checkDelta r, d == (height r - height l)] -- Auxiliary functions to build tree from a list, as same as BST fromList::(Ord a)=[a] - AVLTree a fromList = foldl insert Empty toList :: (AVLTree a) - [a] toList Empty = [] toList (Br l k r _) = toList l ++ [k] ++ toList r -- test prop_bst :: (Ord a, Num a) = [a] - Bool prop_bst xs = (L.sort $ L.nub xs) == (toList $ fromList xs) prop_avl :: (Ord a, Num a) = [a] - Bool prop_avl = isAVL . fromList . L.nub And here are my result in ghci: *AVLTree test prop_avl OK, passed 100 tests. The program is available in github: http://www.google.com/url?sa=Dq=https://github.com/liuxinyu95/AlgoXY/blob/algoxy/datastruct/tree/AVL-tree/src/AVLTree.hs I haven't provided delete function yet. Cheers. -- Larry, LIU https://github.com/liuxinyu95/AlgoXY ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] QuickCheck, (Ord a)= [a] - Property problem
Hi, I tested with Haskell platform 2011 with QuickCheck 2.4.0.1. It produced 100 cases passed, but can't report failed case. verboseCheck still told me that [(), (), ... ()] are generated as instance to (Ord a) The only way is to specify the non-ambitious type for example, Int, like below: test (prop_foo::[Int]-Property) Cheers. -- Larry. On Apr 22, 5:56 am, Nick Smallbone nick.smallb...@gmail.com wrote: larry.liuxinyu liuxiny...@gmail.com writes: Somebody told me that: Eduard Sergeev • BTW, more recent QuickCheck (from Haskell Platform 2011.2.0.X - contains QuickCheck-2.4.0.1) seems to identifies the problem correctly: *** Failed! Falsifiable (after 3 tests and 2 shrinks): [0,1] False I don't think this can be true: the problem occurs in GHCi and there's no way for QuickCheck to detect it. And when I tested it I got the same problem. There must be some difference between the properties you both tested... Nick ___ Haskell-Cafe mailing list Haskell-C...@haskell.orghttp://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] QuickCheck, (Ord a)= [a] - Property problem
Hi, I found there is similar question as: http://groups.google.com/group/haskell-cafe/browse_thread/thread/7439262e9ac80dd2/91ca18e11ff00649?lnk=gstq=QuickCheck+Ord+a#91ca18e11ff00649 I am still think it's very strange. For example: prop_foo :: (Ord a) = [a] - Property prop_foo xs = not (null xs) == maximum xs == minimum xs This is an extreme case that the property is always wrong. However, QuickCheck produces: *Main test prop_foo OK, passed 100 tests. Why this happen? If I use verboseCheck, I can find the sample test data are as the following: *MainverboseCheck prop_foo ... 97: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),()] 98: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(), (),(),(),(),(),(),()] 99: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(), (),(),()] OK, passed 100 tests. So Arbitrary () is generated as the instance of Ord a. I have to change the prototype as prop_foo :: (Num a) = [a] - Property This works at least, However, since 'a''b', they are order-able, what if I want to test prop_foo works for char? I am using Haskell Platform (version 6.10.4), with QuickCheck version 1.2.0.0 Thanks -- Larry, LIU Xinyu https://sites.google.com/site/algoxy/home ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] QuickCheck, (Ord a)= [a] - Property problem
Hi, Thanks a lot. The following type protocol also works. prop_foo :: (Ord a)=(Num a) = [a] - Property Somebody told me that: Eduard Sergeev • BTW, more recent QuickCheck (from Haskell Platform 2011.2.0.X - contains QuickCheck-2.4.0.1) seems to identifies the problem correctly: *** Failed! Falsifiable (after 3 tests and 2 shrinks): [0,1] False -- Larry On Apr 20, 11:36 pm, Nick Smallbone nick.smallb...@gmail.com wrote: larry.liuxinyu liuxiny...@gmail.com writes: prop_foo :: (Ord a) = [a] - Property prop_foo xs = not (null xs) == maximum xs == minimum xs This is an extreme case that the property is always wrong. However, QuickCheck produces: *Main test prop_foo OK, passed 100 tests. Why this happen? If I use verboseCheck, I can find the sample test data are as the following: *MainverboseCheck prop_foo ... 97: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),()] 98: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(), (),(),(),(),(),(),()] 99: [(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(), (),(),()] OK, passed 100 tests. This is an unfortunate feature of GHCi: if the thing you want to evaluate has a polymorphic type then all the type variables default to (), see: http://www.haskell.org/ghc/docs/7.0.3/html/users_guide/interactive-ev... So prop_foo is only tested for lists of (). Nasty. The usual way to work around it is to declare all your properties monomorphic, so write: prop_foo :: [Integer] - Property This works at least, However, since 'a''b', they are order-able, what if I want to test prop_foo works for char? Testing with Integers should always[*] be enough because of parametricity. Nick [*] For certain values of always :) ___ Haskell-Cafe mailing list Haskell-C...@haskell.orghttp://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Haskell KMP(Knuth-Morris-Pratt) algorithm
Hi,
I read about some KMP implementation in Haskell including:
[1] Richard Bird. ``Pearls of Functional algorithm design''
[2] http://twan.home.fmf.nl/blog/haskell/Knuth-Morris-Pratt-in-Haskell.details
[3] http://www.haskell.org/haskellwiki/Runtime_compilation
[4] LazyString version
[1] builds a infinite lazy state transfer trees, while [3] uses index
to build overlap table.
I created a version which isn't as efficient as in [1]. Just for fun:
failure :: (Eq a)= ([a], [a]) - ([a], [a])
failure ([], ys) = ([], ys)
failure (xs, ys) = fallback (init xs) (last xs:ys) where
fallback as bs | as `isSuffixOf` xs = (as, bs)
| otherwise = fallback (init as) (last as:bs)
kmpSearch2 :: (Eq a) = [a] - [a] -[Int]
kmpSearch2 ws txt = snd $ foldl f (([], ws), []) (zip txt [1..]) where
f (p@(xs, (y:ys)), ns) (x, n) | x == y = if ys==[] then ((xs++[y],
ys), ns++[n])
else ((xs++[y], ys), ns)
| xs == [] = (p, ns)
| otherwise = f (failure p, ns) (x,
n)
f (p, ns) e = f (failure p, ns) e
The function failure just follows the idea that in case (xs, ys) fails
matching some letter c in text,
where xs++ys = pattern and c!= head ys, it means we must fallback to
(xs', ys') so that
xs' = longest { s: s is prefix of xs AND s is suffix of xs }
The bad thing is that failure can't memorize what it has compute
before, for example, as pattern = ababc
and we fails at (abab, c), then we call function failure to get
the new one as (ab, abc).
After several matches, we fails again at (abab, c), failure can't
just return (ab, abc) what it has
been compute already. It has too do the same work again.
Function f inside kmpSearch2 is in fact a state-transfer function. If
we try to use some data structure (for example tree) to memorize the
results which failure function calculated, we can finally reach to the
idea in [1].
--
LIU
http://sites.google.com/site/algoxy/
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Haskell KMP(Knuth-Morris-Pratt) algorithm
Hi,
Here is Richard Bird's version for reference. I changed it a bit.
data State a = E | S a (State a) (State a)
matched (S (_, []) _ _) = True
matched _ = False
kmpSearch4 :: (Eq a) = [a] - [a] - [Int]
kmpSearch4 ws txt = snd $ foldl tr (root, []) (zip txt [1..]) where
root = build E ([], ws)
build fails (xs, []) = S (xs, []) fails E
build fails s@(xs, (y:ys)) = S s fails succs where
succs = build' (fst (tr (fails, []) (y, 0))) (xs++[y], ys)
tr (E, ns) _ = (root, ns)
tr ((S (xs, ys) fails succs), ns) (x, n)
| [x] `isPrefixOf` ys = if matched succs then (succs, ns++[n])
else (succs, ns)
| otherwise = tr (fails, ns) (x, n)
In the program, tr is the transfer function applied to the state tree.
And build function is used to build the automaton.
Best regards.
--
LIU
On Mar 3, 5:25 pm, larry.liuxinyu liuxiny...@gmail.com wrote:
Hi,
I read about some KMP implementation in Haskell including:
[1] Richard Bird. ``Pearls of Functional algorithm design''
[2]http://twan.home.fmf.nl/blog/haskell/Knuth-Morris-Pratt-in-Haskell.de...
[3]http://www.haskell.org/haskellwiki/Runtime_compilation
[4] LazyString version
[1] builds a infinite lazy state transfer trees, while [3] uses index
to build overlap table.
I created a version which isn't as efficient as in [1]. Just for fun:
failure :: (Eq a)= ([a], [a]) - ([a], [a])
failure ([], ys) = ([], ys)
failure (xs, ys) = fallback (init xs) (last xs:ys) where
fallback as bs | as `isSuffixOf` xs = (as, bs)
| otherwise = fallback (init as) (last as:bs)
kmpSearch2 :: (Eq a) = [a] - [a] -[Int]
kmpSearch2 ws txt = snd $ foldl f (([], ws), []) (zip txt [1..]) where
f (p@(xs, (y:ys)), ns) (x, n) | x == y = if ys==[] then ((xs++[y],
ys), ns++[n])
else ((xs++[y], ys), ns)
| xs == [] = (p, ns)
| otherwise = f (failure p, ns) (x,
n)
f (p, ns) e = f (failure p, ns) e
The function failure just follows the idea that in case (xs, ys) fails
matching some letter c in text,
where xs++ys = pattern and c!= head ys, it means we must fallback to
(xs', ys') so that
xs' = longest { s: s is prefix of xs AND s is suffix of xs }
The bad thing is that failure can't memorize what it has compute
before, for example, as pattern = ababc
and we fails at (abab, c), then we call function failure to get
the new one as (ab, abc).
After several matches, we fails again at (abab, c), failure can't
just return (ab, abc) what it has
been compute already. It has too do the same work again.
Function f inside kmpSearch2 is in fact a state-transfer function. If
we try to use some data structure (for example tree) to memorize the
results which failure function calculated, we can finally reach to the
idea in [1].
--
LIUhttp://sites.google.com/site/algoxy/
___
Haskell-Cafe mailing list
Haskell-C...@haskell.orghttp://www.haskell.org/mailman/listinfo/haskell-cafe
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Draw K-ary forest in dot script
Hi,
I wrote a Haskell program to parse K-ary forest and convert it to dot script
(Graphviz).
Here is the literate program.
-- First is some stuff imported:
module Main where
import System.Environment (getArgs)
import Text.ParserCombinators.Parsec
import Control.Monad (mapM_)
import Data.List (concatMap, intercalate)
import System.IO (writeFile)
import Data.Char (isSpace)
-- For each tree in the forest, it is described in pre-order.
-- Example description string of a forest of CLRS[1] Figure 19.5(a):
-- (12), (7, (25)), (15, (28, (41)), (33))
-- Definition of K-ary node
data Node a = Node { root :: a
, children :: [Node a]} deriving (Eq, Show)
-- Definition of Forest
type Forest a = [Node a]
-- parsers
-- a forest is a list of trees separate by ','
forest = do
ts - node `sepBy` (char ',')
return ts
-- a node contains a key then followed by a children forest or nothing (leaf
case)
node = do
char '('
elem - key
ts - (try (char ',')forest) | return []
char ')'
return (Node elem ts)
-- a key is just a plain literate string.
key = many (noneOf ,())
-- Command line arguments handling
parseArgs :: [String] - (String, String)
parseArgs [fname, s] = (fname, s)
parseArgs _ = error wrong usage\nexample:\nfr2dot output.dot \(12), (7,
(25)), (15, ((28, (41)), 33))\
-- A simplified function to generate dot script from parsed result.
toDot f = forestToDot f t True
-- a handy function to convert children of a K-ary tree to dot script
treesToDot ts prefix = forestToDot ts prefix False
-- convert a forest to dot script
forestToDot [] _ _ =
forestToDot [t] prefix _ = nodeToDot t prefix
forestToDot ts@(_:_:_) prefix lnk =
(concatMap (\t-nodeToDot t prefix) ts) ++ consRoot
where
consRoot = {rank=same ++ ns ++ vis ++ }\n
ns = intercalate - $ map (\t - prefix ++ root t) ts
vis = if lnk then else [style=invis]
-- convert a node to dot script
nodeToDot (Node x ts) prefix =
prefix'++[label=\++x++\];\n ++
(treesToDot ts prefix') ++
(defCons ts prefix')
where prefix' = prefix ++ x
-- define connections among nodes in dot format
defCons ts prefix = concatMap f ts where
f (Node x _) = prefix++-++prefix++x++;\n
-- generate dot script from a parsed forest
genDot fname (Right f) = writeFile fname dots putStrLn dots
where
dots = digraph G{\n\tnode[shape=circle]\n++(addTab $ toDot f)++}
addTab s = unlines $ map (\t++) (lines s)
main = do
args - getArgs
let (fname, s) = parseArgs args
genDot fname (parse forest unknown (filter (not.isSpace) s))
-- END
I tested with the following simple cases:
./fr2dot foo.dot (12), (7, (25)), (15, (28, (41)), (33))
./fr2dot bar.dot (18), (3, (37)), (6, (8, (30, (45, (55)), (32)), (23,
(24)), (22)), (29, (48, (50)), (31)), (10, (17)), (44))
Run the following commands can convert to PNG files:
./dot -Tpng -o foo.png foo.dot
./dot -Tpng -o bar.png bar.dot
Reference:
[1]. Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford
Stein. ``Introduction to Algorithms, Second Edition''. The MIT Press, 2001.
ISBN: 0262032937.
Best regards.
--
Larry, LIU
https://sites.google.com/site/algoxy/home
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Fibonacci Heap without using Monad
Hi,
I checked the current Fibonacci Queue in Hackage DB:
http://hackage.haskell.org/packages/archive/pqueue-mtl/1.0.7/doc/html/src/Data-Queue-FibQueue.html#FQueue
And a history email for Okasaki in 1995:
http://darcs.haskell.org/nofib/gc/fibheaps/orig
The hardest part is how to consolidate all unordered binomial trees in
deleteMin.
In imperative implementation, there is a elegant algorithm introduced
in Chapter 20 of CLRS[1].
How to achieve it in Functional way is the key point of solve this
problem.
If we have a list of trees with rank [2, 1, 1, 4, 8, 1, 1, 2, 4], we
need first meld the trees with same rank, and recursively doing that
until there are no two trees with same rank. Here is a simple function
can do this:
consolidate:: (Num a)=[a] - [a]
consolidate xs = foldl meld [] xs where
meld :: (Num a)=[a] - a - [a]
meld [] x = [x]
meld (x':xs) x = if x == x' then meld xs (x+x')
else x:x':xs
Generalize the `+` to link and `==` to compare rank yields the
solution.
Below are my literate source code with some description. For the
details of Binomial heap, please refer to Okasaki's ``Purely
Functional data structures''[2].
-- Definition
-- Since Fibonacci Heap can be achieved by applying lazy strategy
-- to Binomial heap. We use the same definition of tree as the
-- Binomial heap. That each tree contains:
-- a rank (size of the tree)
-- the root value (the element)
-- and the children (all sub trees)
data BiTree a = Node { rank :: Int
, root :: a
, children :: [BiTree a]} deriving (Eq, Show)
-- Different with Binomial heap, Fibonacci heap is consist of
-- unordered binomial trees. Thus in order to access the
-- minimum value in O(1) time, we keep the record of the tree
-- which holds the minimum value out off the other children trees.
-- We also record the size of the heap, which is the sum of all ranks
-- of children and minimum tree.
data FibHeap a = E | FH { size :: Int
, minTree :: BiTree a
, trees :: [BiTree a]} deriving (Eq, Show)
-- Auxiliary functions
-- Singleton creates a leaf node and put it as the only tree in the
heap
singleton :: a - FibHeap a
singleton x = FH 1 (Node 1 x []) []
-- Link 2 trees with SAME rank R to a new tree of rank R+1, we re-use
the code
-- for Binomial heaps
link :: (Ord a) = BiTree a - BiTree a - BiTree a
link t1@(Node r x c1) t2@(Node _ y c2)
| xy = Node (r+1) x (t2:c1)
| otherwise = Node (r+1) y (t1:c2)
-- Insertion, runs in O(1) time.
insert :: (Ord a) = FibHeap a - a - FibHeap a
insert h x = merge h (singleton x)
-- Merge, runs in O(1) time.
-- Different from Binomial heap, we don't consolidate the sub trees
-- with the same rank, we delayed it later when performing delete-
Minimum.
merge:: (Ord a) = FibHeap a - FibHeap a - FibHeap a
merge h E = h
merge E h = h
merge h1@(FH sz1 minTr1 ts1) h2@(FH sz2 minTr2 ts2)
| root minTr1 root minTr2 = FH (sz1+sz2) minTr1 (minTr2:ts2+
+ts1)
| otherwise = FH (sz1+sz2) minTr2 (minTr1:ts1++ts2)
-- Find Minimum element in O(1) time
findMin :: (Ord a) = FibHeap a - a
findMin = root . minTree
-- deleting, Amortized O(lg N) time
-- Auxiliary function
-- Consolidate unordered Binomial trees by meld all trees in same rank
-- O(lg N) time
consolidate :: (Ord a) = [BiTree a] - [BiTree a]
consolidate ts = foldl meld [] ts where
meld [] t = [t]
meld (t':ts) t = if rank t' == rank t then meld ts (link t t')
else t:t':ts
-- Find the tree which contains the minimum element.
-- Returns the minimum element tree and the left trees as a pair
-- O(lg N) time
extractMin :: (Ord a) = [BiTree a] - (BiTree a, [BiTree a])
extractMin [t] = (t, [])
extractMin (t:ts) = if root t root t' then (t, ts)
else (t', t:ts')
where
(t', ts') = extractMin ts
-- delete function
deleteMin :: (Ord a) = FibHeap a - FibHeap a
deleteMin (FH _ (Node _ x []) []) = E
deleteMin h@(FH sz minTr ts) = FH (sz-1) minTr' ts' where
(minTr', ts') = extractMin $ consolidate (children minTr ++ ts)
-- Helper functions
fromList :: (Ord a) = [a] - FibHeap a
fromList xs = foldl insert E xs
-- general heap sort function, can be re-used for any heap
heapSort :: (Ord a) = [a] - [a]
heapSort = hsort . fromList where
hsort E = []
hsort h = (findMin h):(hsort $ deleteMin h)
-- test
testFromList = fromList [16, 14, 10, 8, 7, 9, 3, 2, 4, 1]
testHeapSort = heapSort [16, 14, 10, 8, 7, 9, 3, 2, 4, 1]
Below are the test results in GHC.
*FibonacciHeap testFromList
FH {size = 10, minTree = Node {rank = 1, root = 1, children = []},
trees = [Node {rank = 1, root = 2, children = []},Node {rank = 1, root
= 4, children = []},Node {rank = 1, root = 3, children = []},Node
{rank = 1, root = 7, children = []},Node {rank = 1, root = 9, children
= []},Node {rank = 1, root = 8, children = []},Node {rank = 1, root =
10, children = []},Node {rank = 1, root
Re: [Haskell-cafe] Fibonacci Heap without using Monad
Hi, In CLRS, there are algorithms about DECREASE-KEY and DELETE-NODE. However, in the Functional approach, I didn't find corresponding solution. One approach may just mark the node as `deleted' and when pops the top element from the heap, we repeat it until find a unmarked node. -- LIU https://sites.google.com/site/algoxy/home ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fibonacci Heap without using Monad
Hi, Sorry for there is a bug in my previous post. The example consolidate function for number should be like this: consolidate xs = foldl meld [] xs where meld [] x = [x] meld (x':xs) x | x == x' = meld xs (x+x') | x x' = x:x':xs | otherwise = x': meld xs x The bug happens in my previous mail like below. before fixing consolidate [2, 1, 1, 32, 4, 8, 1, 1, 2, 4] [16,4,32,4] after fixing-- consolidate [2, 1, 1, 32, 4, 8, 1, 1, 2, 4] [8, 16, 32] Therefore, the consolidate function for unordered binomial trees should be modified as the following respectively. consolidate :: (Ord a) = [BiTree a] - [BiTree a] consolidate ts = foldl meld [] ts where meld [] t = [t] meld (t':ts) t | rank t == rank t' = meld ts (link t t') | rank t rank t' = t:t':ts | otherwise = t' : meld ts t I am sorry for this mistake. The updated source code can be found from here: https://sites.google.com/site/algoxy/otherheaps/otherheaps.zip Happy new year. -- Larry, LIU ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Splay tree in a pattern matching way
Hi, I just tried a smoke test case, seems the tree is balanced enough: import System.Random lookup :: (Ord a) = STree a - a - STree a lookup E _ = E lookup t@(Node l x r) y | x == y= t | x y = splay (Node (lookup l y) x r) y | otherwise = splay (Node l x (lookup r y)) y testSplay = do xs - sequence (replicate 100 (randomRIO(1, 10))) putStrLn $ show (foldl lookup t xs) where t = foldl insert (E::STree Int) [1..10] Where STree a and insert are defined as in my previous email, except I changed `xy' to `xy' (sorry for the mistake). By running the testSplay I got a BST like this: *SplayHeap testSplay Node (Node E 1 (Node (Node E 2 E) 3 (Node (Node E 4 (Node E 5 E)) 6 (Node E 7 E 8 (Node (Node E 9 E) 10 E) So I got a bit more confident about the correctness of my code. Best regards -- Larry, LIU Xinyu http://sites.google.com/site/algoxy On Oct 24, 11:06 am, Xinyu LIU liuxiny...@gmail.com wrote: Hi, I checked the Hackage Splay Tree implementation from.http://hackage.haskell.org/packages/archive/TreeStructures/0.0.1/doc/... Basically it's based on the strategy mentioned in Okasaki's ``Purely Functional Programming'', Chapter 5.4. That in case we traverse twice in left (or right), we rotate the tree, so in long term of view, the tree turns more and more balanced. There are 3 cases for splay: zig-zig case, zig-zag case and zig case according tohttp://en.wikipedia.org/wiki/Splay_tree So I wonder if Splay tree can also be implemented by using pattern matching in a similar way as red black tree. Here is a draft source code: data STree a = E -- Empty | Node (STree a) a (STree a) -- left, element, right deriving (Eq, Show) -- splay by pattern matching splay :: (Eq a) = STree a - a -STree a -- zig-zig splay t@(Node (Node (Node a x b) p c) g d) y = if x == y then Node a x (Node b p (Node c g d)) else t splay t@(Node a g (Node b p (Node c x d))) y = if x == y then Node (Node (Node a g b) p c) x d else t -- zig-zag splay t@(Node (Node a p (Node b x c)) g d) y = if x == y then Node (Node a p b) x (Node c g d) else t splay t@(Node a g (Node (Node b x c) p d)) y = if x == y then Node (Node a g b) x (Node c p d) else t -- zig splay t@(Node (Node a x b) p c) y = if x == y then Node a x (Node b p c) else t splay t@(Node a p (Node b x c)) y = if x == y then Node (Node a p b) x c else t -- otherwise splay t _ = t -- insert by using pattern matching insert :: (Ord a) = STree a - a - STree a insert E y = Node E y E insert (Node l x r) y | x y = splay (Node (insert l y) x r) y | otherwise = splay (Node l x (insert r y)) y I am not quite sure if the solution is correct and what about the efficiency of this code. Regards. -- Larry. LIU Xinyuhttp://sites.google.com/site/algoxyhttp://sites.google.com/site/algoxy/home ___ Haskell-Cafe mailing list haskell-c...@haskell.orghttp://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
