86f6716306f3f6da0351b452c9e9fda85a14264b
[packages/containers.git] / Data / Sequence / Internal / sorting.md
1 # Sorting
2
3 ## Unstable Sorting
4
5 Unstable sorting is performed by a heap sort implementation based on
6 pairing heaps.  Because the internal structure of sequences is quite
7 varied, it is difficult to get blocks of elements of roughly the same
8 length, which would improve merge sort performance.  Pairing heaps,
9 on the other hand, are relatively resistant to the effects of merging
10 heaps of wildly different sizes, as guaranteed by its amortized
11 constant-time merge operation.  Moreover, extensive use of SpecConstr
12 transformations can be done on pairing heaps, especially when we're
13 only constructing them to immediately be unrolled.
14
15 On purely random sequences of length 50000, with no RTS options,
16 I get the following statistics, in which heapsort is about 42.5%
17 faster:  (all comparisons done with -O2)
18
19 Times (ms)        |  min  |  mean  | +/-sd | median |  max
20 ------------------|-------|--------|-------|--------|-------
21 to/from list:     |103.802| 108.572|  7.487| 106.436|143.339
22 unstable heapsort:| 60.686|  62.968|  4.275|  61.187| 79.151
23
24 Heapsort, it would seem, is less of a memory hog than Data.List.sortBy.
25 The gap is narrowed when more memory is available, but heapsort still
26 wins, 15% faster, with +RTS -H128m:
27
28 Times (ms)        |  min  | mean | +/-sd | median |  max
29 ------------------|-------|------|-------|--------|-------
30 to/from list:     | 42.692|45.074|  2.596|  44.600| 56.601
31 unstable heapsort:| 37.100|38.344|  3.043|  37.715| 55.526
32
33 In addition, on strictly increasing sequences the gap is even wider
34 than normal; heapsort is 68.5% faster with no RTS options:
35
36 Times (ms)        |  min  | mean | +/-sd | median |  max
37 ------------------|-------|------|-------|--------|-------
38 to/from list:     | 52.236|53.574|  1.987|  53.034| 62.098
39 unstable heapsort:| 16.433|16.919|  0.931|  16.681| 21.622
40
41 This may be attributed to the elegant nature of the pairing heap.
42
43 wasserman.louis@gmail.com, 7/20/09
44
45 ----------------------------------------------------------------------
46
47 David Feuer wrote an unstable sort for arbitrary traversables,
48 https://www.reddit.com/r/haskell/comments/63a4ea/fast_total_sorting_of_arbitrary_traversable/,
49 which turned out to be competitive with the unstable sort here.
50 Feuer suggested that this indicated some room to improve on the
51 unstable sort.
52
53 The two main improvements to the original function are a specialize
54 pragma on replicateA (this gives a 26.5% speedup) and removal of the
55 intermediate list (a further 11.7% speedup). These numbers are all on
56 purely random sequences of length 50000:
57
58 Times (ms)       | min | est | max |std dev|  r²
59 -----------------|-----|-----|-----|-------|-----
60 to/from list:    |70.90|72.44|75.07|  2.224|0.998
61 7/20/09 heapsort:|59.84|61.44|63.08|  1.554|0.998
62 7/20/09 w/pragma:|44.22|45.14|46.25|  1.631|0.998
63 4/30/17 heapsort:|38.21|39.86|40.87|  1.203|0.996
64
65 It should also be noted that Data.List.sortBy has become
66 significantly quicker. Data.List.sortBy also now recognizes strictly
67 increasing sequences, making it much quicker for that case:
68
69 Times (ms)       | min | est | max |std dev|  r²
70 -----------------|-----|-----|-----|-------|-----
71 to/from list:    |7.140|7.351|7.634|  0.335|0.993
72 7/20/09 heapsort:|19.52|19.78|20.13|  0.445|0.999
73 7/20/09 w/pragma:|8.050|8.271|8.514|  0.357|0.995
74 4/30/17 heapsort:|7.240|7.612|7.925|  0.389|0.991
75
76 Another happy result of the specialization of 'replicateA' is that
77 the stable sort seems to speed up by 10-20%, and 'iterateN' looks
78 like it's about three times as fast.
79
80 mail@doisinkidney.com, 4/30/17
81
82 ## Stable Sorting
83
84 Stable sorting was previously accomplished by converting to a list,
85 applying Data.List.sort, and rebuilding the sequence. Data.List.sort is
86 designed to maximize laziness, which doesn't apply for Data.Sequence,
87 and it can't take advantage of the structure of the finger tree. As a
88 result, simply tagging each element with its position, then applying
89 the unstable sort (using the tag to discriminate between elements for
90 which the comparator is equal) is faster. The current implementation
91 doesn't use the actual `unstableSort`: to perform the building of the
92 queue and tagging in one pass, a specialized version is used.
93
94 The algorithm is effectively the same as the unstable sorts, except
95 the queue is constructed while giving each element a tag.
96
97 It's quicker than the old implementation (which used Data.List.sort)
98 in the general case (all times are on sequences of length 50000):
99
100 Times (ms)          | min | est | max |std dev|  r²
101 --------------------|-----|-----|-----|-------|-----
102 to/from list:       |64.23|64.50|64.81|  0.432|1.000
103 1/11/18 stable heap:|38.87|39.40|40.09|  0.457|0.999
104
105 Slightly slower in the case of already sorted lists:
106
107 Times (ms)          | min | est | max |std dev|  r²
108 --------------------|-----|-----|-----|-------|-----
109 to/from list:       |6.806|6.861|6.901|  0.234|1.000
110 1/11/18 stable heap:|8.211|8.268|8.328|  0.111|1.000
111
112 And quicker in the case of lists sorted in reverse:
113
114 Times (ms)          | min | est | max |std dev|  r²
115 --------------------|-----|-----|-----|-------|-----
116 to/from list:       |26.79|28.34|30.55|  1.219|0.988
117 1/11/18 stable heap:|9.405|10.13|10.91|  0.670|0.977
118
119 Interestingly, the stable sort is now competitive with the unstable:
120
121 Times (ms)| min | est | max |std dev|  r²
122 ----------|-----|-----|-----|-------|-----
123 unstable: |34.71|35.10|35.38|  0.845|1.000
124 stable:   |38.84|39.22|39.59|  0.564|0.999
125
126 And even beats it in the case of already-sorted lists:
127
128 Times (ms)| min | est | max |std dev|  r²
129 ----------|-----|-----|-----|-------|-----
130 unstable: |8.457|8.499|8.536|  0.069|1.000
131 stable:   |8.160|8.230|8.334|  0.158|0.999
132
133 mail@doisinkidney.com, 1/11/18
134
135 ## sortOn Functions
136
137 The `sortOn` and `unstableSortOn` functions perform the Schwartzian transform, however instead of the following implementation:
138
139 ```haskell
140 sortOn f = fmap snd . sortBy (conparing fst) . fmap (\x -> (f x, x))
141 ```
142
143 The `fmap`s are fused manually with the creation of the queue, avoiding the two extra traversals. It still suffers a slowdown of roughly 20%:
144
145 Times (ms)     | min | est | max |std dev|  r²
146 ---------------|-----|-----|-----|-------|-----
147 unstableSortOn |43.68|44.58|45.95|  0.677|0.999
148 unstableSort   |36.55|37.43|38.33|  0.533|0.999
149 sortOn         |48.22|49.03|50.09|  1.110|0.998
150 sort           |41.81|43.17|45.31|  1.172|0.996
151
152 The heaps are also specialized to avoid the creation of a tuple.
153
154 ## Other Heaps
155
156 The pairing heap seems to particularly suit the structure of the finger tree, as other heaps have not managed to beat it. Specifically, when compared to a skew heap:
157
158 ```haskell
159 unstableSortBy :: (a -> a -> Ordering) -> Seq a -> Seq a
160 unstableSortBy cmp (Seq xs) =
161     execState (replicateA (size xs) (popMin cmp)) (toSkew cmp (Seq xs))
162
163 data Skew a = Nil | Br a !(Skew a) !(Skew a)
164
165 popMin :: (e -> e -> Ordering) -> State (Skew e) e
166 popMin cmp = State unrollPQ'
167   where
168     {-# INLINE unrollPQ' #-}
169     unrollPQ' (Br x ls rs) = (mergeSkew cmp ls rs, x)
170
171 toSkew :: (e -> e -> Ordering) -> Seq e -> Skew e
172 toSkew cmp (Seq xs') = toSkewTree cmp (\(Elem a) -> Br a Nil Nil) xs'
173   where
174     toSkewTree :: (b -> b -> Ordering) -> (a -> Skew b) -> FingerTree a -> Skew b
175     toSkewTree _ _ EmptyT = Nil
176     toSkewTree _ f (Single xs) = f xs
177     toSkewTree cmp f (Deep n pr m sf) = pr' <+> sf' <+> m'
178       where
179         pr' = toSkewDigit cmp f pr
180         sf' = toSkewDigit cmp f sf
181         m' = toSkewTree cmp (toSkewNode cmp f) m
182         (<+>) = mergeSkew cmp
183     toSkewDigit :: (b -> b -> Ordering) -> (a -> Skew b) -> Digit a -> Skew b
184     toSkewDigit cmp f dig =
185         case dig of
186             One a -> f a
187             Two a b -> f a <+> f b
188             Three a b c -> f a <+> f b <+> f c
189             Four a b c d -> (f a <+> f b) <+> (f c <+> f d)
190       where
191         (<+>) = mergeSkew cmp
192     toSkewNode cmp f node =
193         case node of
194             Node2 _ a b -> f a <+> f b
195             Node3 _ a b c -> f a <+> f b <+> f c
196       where
197         (<+>) = mergeSkew cmp
198
199 mergeSkew :: (a -> a -> Ordering) -> Skew a -> Skew a -> Skew a
200 mergeSkew cmp Nil ys = ys
201 mergeSkew cmp xs Nil = xs
202 mergeSkew cmp h1@(Br x lx rx) h2@(Br y ly ry)
203   | cmp x y == GT = Br y (mergeSkew cmp h1 ry) ly
204   | otherwise     = Br x (mergeSkew cmp h2 rx) lx
205 ```
206
207 The pairing heap implementation is faster in every aspect:
208
209 ```
210 benchmarking 1000000/unsorted/pairing
211 time                 2.005 s    (NaN s .. 2.102 s)
212                      1.000 R²   (0.998 R² .. 1.000 R²)
213 mean                 2.069 s    (2.060 s .. 2.075 s)
214 std dev              9.340 ms   (0.0 s .. 10.67 ms)
215 variance introduced by outliers: 19% (moderately inflated)
216              
217 benchmarking 1000000/unsorted/skew
218 time                 2.042 s    (1.637 s .. 2.267 s)
219                      0.995 R²   (0.990 R² .. NaN R²)
220 mean                 2.165 s    (2.065 s .. 2.217 s)
221 std dev              87.10 ms   (0.0 s .. 91.26 ms)
222 variance introduced by outliers: 19% (moderately inflated)
223              
224 benchmarking 1000000/ascending/pairing
225 time                 191.4 ms   (187.8 ms .. 193.5 ms)
226                      1.000 R²   (0.999 R² .. 1.000 R²)
227 mean                 197.0 ms   (194.7 ms .. 200.0 ms)
228 std dev              3.221 ms   (2.441 ms .. 3.924 ms)
229 variance introduced by outliers: 14% (moderately inflated)
230              
231 benchmarking 1000000/ascending/skew
232 time                 232.3 ms   (227.0 ms .. 238.9 ms)
233                      0.999 R²   (0.997 R² .. 1.000 R²)
234 mean                 233.9 ms   (230.6 ms .. 236.2 ms)
235 std dev              3.678 ms   (2.790 ms .. 4.777 ms)
236 variance introduced by outliers: 14% (moderately inflated)
237              
238 benchmarking 1000000/descending/pairing
239 time                 204.6 ms   (190.2 ms .. 214.1 ms)
240                      0.998 R²   (0.991 R² .. 1.000 R²)
241 mean                 208.4 ms   (204.1 ms .. 210.6 ms)
242 std dev              4.051 ms   (1.299 ms .. 5.288 ms)
243 variance introduced by outliers: 14% (moderately inflated)
244              
245 benchmarking 1000000/descending/skew
246 time                 229.9 ms   (212.7 ms .. 240.1 ms)
247                      0.998 R²   (0.996 R² .. 1.000 R²)
248 mean                 238.8 ms   (231.3 ms .. 241.4 ms)
249 std dev              5.006 ms   (269.0 μs .. 6.151 ms)
250 variance introduced by outliers: 16% (moderately inflated)
251 ```
252