[project @ 1997-03-14 08:02:40 by simonpj]
[nofib.git] / spectral / fibheaps / Main.lhs
1 Date:    Tue, 04 Jul 1995 13:10:58 -0400
2 From:    Chris_Okasaki@LOCH.MESS.CS.CMU.EDU
3 To:      simonpj@dcs.gla.ac.uk
4 Subject: Fibonacci Heaps
5
6 As I promised at the Haskell Workshop, here is a sample program
7 using encapsulated state.  I've translated this from SML, but 
8 in doing so, I noticed that in fact accumArray is all the
9 encapsulated state you really need for this application.  In SML,
10 we are forced to use mutable arrays because we don't have such
11 fancy monolithic array "primitives" as accumArray.
12
13 I've written and tested this as a literate Gofer script because I've
14 never been able to get GHC to run under Mach. :-(
15
16 Let me know if you have any problems...
17
18 Chris
19
20 - ----------------------------------------------------------------------
21
22 FIBONACCI HEAPS
23
24 Fibonacci heaps are a priority queue data structure supporting the
25 following operations:
26         O(1) Insert
27         O(1) FindMin
28         O(1) Meld
29     O(log n) DeleteMin
30
31 (In an imperative settting, Fibonacci heaps also support
32         O(1) DecreaseKey (of an indicated element)
33     O(log n) Delete (an indicated element)
34 but these operations are problematic in a functional setting.)
35
36 There is one catch: for the DeleteMin operation, the bounds are
37 amortized instead of worst-case.  This means that the bounds are
38 only guaranteed if you use the data structure in a single-threaded manner.
39 Otherwise, you can take longer than expected by repeatedly going back
40 and operating on an "expensive" version of the data structure.
41
42 (Note: I am currently working on a paper with another student describing
43 a functional priority queue achieving the above bounds in the worst-case
44 instead of amortized.  This data structure may be freely used in a
45 non-single-threaded manner with no ill effects.)
46
47 To understand the implementation of Fibonacci heaps, it is helpful to
48 first understand binomial queues.  See, for example, David King's
49 "Functional Binomial Queues" from the last Glasgow workshop.
50
51 > -- partain
52 >module Main (main) where
53 >import ST
54 >import Array
55 >import System
56
57                          --------------------
58
59 Like binomial queues, Fibonacci heaps are based on heap-ordered
60 binomial trees.
61
62 >data Tree a = Node a [Tree a]
63
64 The degree of a binomial tree is equal to its number of children.  
65 Every binomial tree of degree k has binomial trees of degrees
66 k-1...0 as children, in that order.  It is easy to show that
67 a binomial tree of degree k has size 2^k.
68
69 The fundamental operation on binomial trees is linking, which compares
70 the roots of two binomial trees and makes the larger a child of the 
71 smaller (thus bumping its degree by one).  It is essential that this
72 only be called on binomial trees of equal degree.
73
74 >link (a @ (Node x as)) (b @ (Node y bs)) =
75 >  if x <= y then Node x (b:as) else Node y (a:bs)
76
77 It will also be useful to extract the minimum element from a tree.
78
79 >root (Node x _) = x
80
81 We will frequently need to tag trees with their degrees.
82
83 >type TaggedTree a = (Int,Tree a)
84 >
85 >degree (k, t) = k
86 >tree (k, t) = t
87
88 Given a tagged tree, extract and tag its children.
89
90 >getChildren (n, Node x ts) = zipWith (,) [n-1,n-2 .. ] ts
91
92 Extract the minimum element from a tagged tree.
93
94 >root' = root . tree
95
96                          --------------------
97
98 We also need a type for bags supporting constant time union.  The simple
99 representation given here is sufficient since we will always process bags
100 as a whole.  Note that for this application it is not necessary to
101 filter out occurences of EmptyBag.  Also, for this application order
102 is irrelevant.
103
104 >data Bag a = EmptyBag | ConsBag a (Bag a) | UnionBags (Bag a) (Bag a)
105 >
106 >bagToList b = flatten b []
107 >  where flatten EmptyBag xs = xs
108 >        flatten (ConsBag x b) xs = flatten b (x:xs)
109 >        flatten (UnionBags b1 b2) xs = flatten b1 (flatten b2 xs)
110 >
111 >applyToAll :: (a -> ST s ()) -> Bag a -> ST s ()
112 >applyToAll f EmptyBag = returnST ()
113 >applyToAll f (ConsBag x b) = f x `seqST` applyToAll f b
114 >applyToAll f (UnionBags b1 b2) = applyToAll f b1 `seqST` applyToAll f b2
115
116                          --------------------
117
118 Miscellaneous stuff.
119
120 >log2 1 = 0
121 >log2 n = 1 + log2 (n `div` 2)
122
123 >data MyMaybe a = Zero | One a
124
125                          --------------------
126
127 Since binomial trees only come in certain, fixed sizes, we need some
128 way to represent priority queues of other sizes.  We will do this
129 with a forest of trees summing to the correct size.
130
131 >type Forest a = Bag (TaggedTree a)
132
133 In binomial queues, this forest must be maintained in strictly increasing
134 order of degree.  For Fibonacci heaps, we adopt a more relaxed attitude: 
135 degrees may be repeated and order does not matter.
136
137 To be able to find the minimum element quickly, we keep the tree with the 
138 minimum root outside of the bag.  In addition, at the top level of each heap,
139 we store the total size of the heap.
140
141 >data FibHeap a = EmptyFH | FH Int (TaggedTree a) (Forest a)
142
143                          --------------------
144
145 Now, the following operations are trivial.
146
147 >emptyFH = EmptyFH
148 >
149 >isEmptyFH EmptyFH = True
150 >isEmptyFH (FH _ _ _) = False
151 >
152 >singleFH x = FH 1 (0, Node x []) EmptyBag
153 >
154 >insertFH x xs = meldFH (singleFH x) xs
155 >
156 >minFH EmptyFH = error "minFH EmptyFH"
157 >minFH (FH n tt f) = root' tt
158
159                          --------------------
160
161 Meld achieves its efficiency by simply unioning the two forests.
162
163 >meldFH EmptyFH xs = xs
164 >meldFH xs EmptyFH = xs
165 >meldFH (FH n1 tt1 f1) (FH n2 tt2 f2) =
166 >  if root' tt1 <= root' tt2 then
167 >      FH (n1+n2) tt1 (ConsBag tt2 (UnionBags f1 f2))
168 >  else
169 >      FH (n1+n2) tt2 (ConsBag tt1 (UnionBags f1 f2))
170
171                          --------------------
172
173 Finally, the only hard operation is deleteMin.  After throwing away the
174 minimum element, it repeatedly links trees of equal degree until
175 no such pairs are left.  The most efficient way to do this is with
176 an array.  I give two implementations, one using monadic arrays,
177 the other using accumArray.
178
179 In the first implementation, there are three steps.
180   1. Allocate an array indexed by degrees.
181   2. Insert every tree into this array.  If, when inserting a tree of 
182      degree k, there already exists a tree of degree k, link the
183      two trees and reinsert the new larger tree.
184   3. Transfer the trees into a bag, keeping track of the minimum tree.
185
186 >deleteMinFH EmptyFH = error "deleteMinFH EmptyFH"
187 >deleteMinFH (FH 1 tt f) = EmptyFH
188 >deleteMinFH (FH n tt f) =
189 >  let
190 >    d = log2 (n-1) -- maximum possible degree
191 >
192 >    ins a (i, t) =
193 >        readArray a i `thenST` \e ->
194 >        case e of
195 >          Zero -> writeArray a i (One t)
196 >          One t2 -> writeArray a i Zero `seqST`
197 >                    ins a (i+1, link t t2)
198
199 Note that after inserting all the trees, the array contains trees
200 in the same pattern as the bits of n-1.  Since we know that the
201 highest order bit of n-1 is one, we know that there is a tree in
202 the highest slot of the array.
203
204 >    getMin a =
205 >        readArray a d `thenST` \e ->
206 >        case e of
207 >          Zero -> error "must be One" -- since array is filled as bits of n-1
208 >          One t -> getMin' a d t EmptyBag 0
209 >    getMin' a mini mint b i =
210 >        if i >= d then
211 >          returnST ((mini, mint),b)
212 >        else
213 >          readArray a i `thenST` \e ->
214 >          case e of
215 >            Zero -> getMin' a mini mint b (i+1)
216 >            One t -> if root mint <= root t then
217 >                       getMin' a mini mint (ConsBag (i, t) b) (i+1)
218 >                     else
219 >                       getMin' a i t (ConsBag (mini, mint) b) (i+1)
220 >            
221 >  in 
222 >    runST (newArray (0,d) Zero `thenST` \a ->
223 >           applyToAll (ins a) f `seqST`
224 >           listST (map (ins a) (getChildren tt)) `seqST`
225 >           getMin a `thenST` \ (tt,f) ->
226 >           returnST (FH (n-1) tt f))
227
228                          --------------------
229
230 The second version of deleteMin uses accumArray to group trees of like
231 size.  It then performs the linking and all remaining steps purely 
232 functionally.
233
234 >deleteMinFH' EmptyFH = error "deleteMinFH EmptyFH"
235 >deleteMinFH' (FH 1 tt f) = EmptyFH
236 >deleteMinFH' (FH n tt f) =
237 >  let
238 >    d = log2 (n-1) -- maximum possible degree
239 >
240 >    a = accumArray (flip (:)) [] (0,d) (getChildren tt ++ bagToList f)
241 >
242 >    doLinks (ts:rest) = startup 0 ts rest
243 >      where startup i [] [] = []
244 >            startup i [] (ts:rest) = startup (i+1) ts rest
245 >            startup i ts [] = combine i ts [] []
246 >            startup i ts (next:rest) = combine i ts next rest
247 >
248 >            combine i [] next rest = startup (i+1) next rest
249 >            combine i [t] next rest = (i, t) : startup (i+1) next rest
250 >            combine i (t1:t2:ts) next rest =
251 >                combine i ts (link t1 t2 : next) rest
252 >
253 >    getMin (tt:rest) = foldl chooseMin (tt,EmptyBag) rest
254 >      where chooseMin (tt1,b) tt2 =
255 >                if root' tt1 <= root' tt2 then
256 >                    (tt1,ConsBag tt2 b)
257 >                else
258 >                    (tt2,ConsBag tt1 b)
259 >
260 >    (new_tt,new_f) = getMin (doLinks (elems a))
261 >  in
262 >    FH (n-1) new_tt new_f
263
264                          --------------------
265
266 Testing...
267
268 >fibToList :: Ord a => FibHeap a -> [a]
269 >fibToList xs = if isEmptyFH xs then []
270 >               else minFH xs : fibToList (deleteMinFH xs)
271 >
272 >fibToList' :: Ord a => FibHeap a -> [a]
273 >fibToList' xs = if isEmptyFH xs then []
274 >                else minFH xs : fibToList' (deleteMinFH' xs)
275 >
276 >makeFH :: Ord a => [a] -> FibHeap a
277 >makeFH xs = foldr insertFH emptyFH xs
278 >
279 >fibSort :: Ord a => [a] -> [a]
280 >fibSort = fibToList . makeFH
281 >
282 >fibSort' :: Ord a => [a] -> [a]
283 >fibSort' = fibToList' . makeFH
284 >
285 >randoms :: Int -> [Int]
286 >randoms n = take n (iterate (\seed-> (77*seed+1) `rem` 1024) 1967)
287 >
288 >test n = fibSort (randoms n) == fibSort' (randoms n)
289
290 >--partain
291 >main = getArgs >>= \ [n] ->
292 >       putStr (show (test (read n)))