Algorithms 2.4 P Q

Algorithms
R OBERT S EDGEWICK | K EVIN W AYNE
2.4 P RIORITY Q UEUES
‣ API and elementary implementations
‣ binary heaps
Algorithms
F O U R T H
E D I T I O N
R OBERT S EDGEWICK | K EVIN W AYNE
http://algs4.cs.princeton.edu
‣ heapsort
‣ event-driven simulation
2.4 P RIORITY Q UEUES
‣ API and elementary implementations
‣ binary heaps
Algorithms
R OBERT S EDGEWICK | K EVIN W AYNE
http://algs4.cs.princeton.edu
‣ heapsort
‣ event-driven simulation
Collections
A collection is a data type that stores a group of items.
data type
key operations
data structure
stack
PUSH, POP
linked list, resizing array
queue
ENQUEUE, DEQUEUE
linked list, resizing array
priority queue
INSERT,
DELETE-MAX
binary heap
symbol table
PUT, GET, DELETE
binary search tree, hash table
set
ADD, CONTAINS, DELETE
binary search tree, hash table
“ Show me your code and conceal your data structures, and I shall
continue to be mystified. Show me your data structures, and I won't
usually need your code; it'll be obvious.” — Fred Brooks
3
Priority queue
Collections. Insert and delete items. Which item to delete?
Stack. Remove the item most recently added.
Queue. Remove the item least recently added.
Randomized queue. Remove a random item.
Priority queue. Remove the largest (or smallest) item.
operation argument
insert
insert
insert
remove max
insert
insert
insert
remove max
insert
insert
insert
remove max
return
value size
P
Q
E
Q
X
A
M
X
P
L
E
P
1
2
3
2
3
4
5
4
5
6
7
6
contents
(unordered)
P
P
P
P
P
P
P
P
P
P
P
E
Q
Q
E
E
E
E
E
E
E
E
M
E
X
X
X
M
M
M
M
A
A
A
A
A
A
A
P
M
P
P
P
L
4
Priority queue API
Requirement. Items are generic; they must also be Comparable.
Key must be Comparable
(bounded type parameter)
public class MaxPQ<Key extends Comparable<Key>>
MaxPQ()
MaxPQ(Key[] a)
void insert(Key v)
Key delMax()
create an empty priority queue
create a priority queue with given keys
insert a key into the priority queue
return and remove the largest key
boolean isEmpty()
is the priority queue empty?
Key max()
return the largest key
int size()
number of entries in the priority queue
5
Priority queue: applications
・Event-driven simulation.
・Numerical computation.
・Discrete optimization.
・Artificial intelligence.
・Computer networks.
・Operating systems.
・Data compression.
・Graph searching.
・Number theory.
・Spam filtering.
・Statistics.
[ customers in a line, colliding particles ]
[ reducing roundoff error ]
[ bin packing, scheduling ]
[ A* search ]
[ web cache ]
[ load balancing, interrupt handling ]
[ Huffman codes ]
[ Dijkstra's algorithm, Prim's algorithm ]
[ sum of powers ]
[ Bayesian spam filter ]
[ online median in data stream ]
Generalizes: stack, queue, randomized queue.
6
Priority queue: client example
Challenge. Find the largest M items in a stream of N items.
・Fraud detection:
・NSA monitoring:
isolate $$ transactions.
flag most suspicious documents.
N huge, M large
Constraint. Not enough memory to store N items.
% more transactions.txt
Turing
6/17/1990
644.08
vonNeumann 3/26/2002 4121.85
Dijkstra
8/22/2007 2678.40
vonNeumann 1/11/1999 4409.74
Dijkstra
11/18/1995
837.42
Hoare
5/10/1993 3229.27
vonNeumann 2/12/1994 4732.35
Hoare
8/18/1992 4381.21
Turing
1/11/2002
66.10
Thompson
2/27/2000 4747.08
Turing
2/11/1991 2156.86
Hoare
8/12/2003 1025.70
vonNeumann 10/13/1993 2520.97
Dijkstra
9/10/2000
708.95
Turing
10/12/1993 3532.36
Hoare
2/10/2005 4050.20
% java TopM
Thompson
vonNeumann
vonNeumann
Hoare
vonNeumann
5 < transactions.txt
2/27/2000 4747.08
2/12/1994 4732.35
1/11/1999 4409.74
8/18/1992 4381.21
3/26/2002 4121.85
sort key
7
Priority queue: client example
Challenge. Find the largest M items in a stream of N items.
・Fraud detection:
・NSA monitoring:
isolate $$ transactions.
flag most suspicious documents.
Constraint. Not enough memory to store N items.
MinPQ<Transaction> pq = new MinPQ<Transaction>();
use a min-oriented pq
while (StdIn.hasNextLine())
{
String line = StdIn.readLine();
Transaction item = new Transaction(line);
pq.insert(item);
if (pq.size() > M)
pq now contains
pq.delMin();
largest M items
}
Transaction data
type is Comparable
(ordered by $$)
8
Priority queue: client example
Challenge. Find the largest M items in a stream of N items.
implementation
time
space
sort
N log N
N
elementary PQ
MN
M
binary heap
N log M
M
best in theory
N
M
order of growth of finding the largest M in a stream of N items
9
Priority queue: unordered and ordered array implementation
operation argument
insert
insert
insert
remove max
insert
insert
insert
remove max
insert
insert
insert
remove max
return
value size
P
Q
E
Q
X
A
M
X
P
L
E
P
1
2
3
2
3
4
5
4
5
6
7
6
contents
(unordered)
P
P
P
P
P
P
P
P
P
P
P
E
Q
Q
E
E
E
E
E
E
E
E
M
contents
(ordered)
E
X
X
X
M
M
M
M
A
A
A
A
A
A
A
P
M
P
P
P
L
L
L
E
E
P
P
E
E
E
A
A
A
A
A
A
A
Q
P
P
P
E
E
E
E
E
E
E
Q
X
P
M
M
M
L
E
E
X
P
P
P
M
L
L
X
P
P
M
M
P
P
P
P
A sequence of operations on a priority queue
10
Priority queue: implementations cost summary
Challenge. Implement all operations efficiently.
implementation
insert
del max
max
unordered array
1
N
N
ordered array
N
1
1
goal
log N
log N
log N
order of growth of running time for priority queue with N items
11
2.4 P RIORITY Q UEUES
‣ API and elementary implementations
‣ binary heaps
Algorithms
R OBERT S EDGEWICK | K EVIN W AYNE
http://algs4.cs.princeton.edu
‣ heapsort
‣ event-driven simulation
Complete binary tree
Binary tree. Empty or node with links to left and right binary trees.
Complete tree. Perfectly balanced, except for bottom level.
complete binary tree with N = 16 nodes (height = 4)
Property. Height of complete binary tree with N nodes is ⎣lg N⎦.
Pf. Height increases only when N is a power of 2.
13
A complete binary tree in nature
14
Binary heap: representation
Binary heap. Array representation of a heap-ordered complete binary tree.
Heap-ordered binary tree.
・Keys in nodes.
・Parent's key no smaller than
children's keys.
i
a[i]
0
-
1
T
2
S
3
R
S
R
4
P
5
N
6
O
7
A
P
N
O
A
8
E
9 10 11
I H G
E
I
T
Array representation.
・Indices start at 1.
・Take nodes in level order.
・No explicit links needed!
1
8 E
9 I
G
T
3 R
2 S
4 P
H
5 N
10 H
6 O
7 A
11 G
Heap representations
15
Binary heap: properties
Proposition. Largest key is a[1], which is root of binary tree.
Proposition. Can use array indices to move through tree.
・Parent of node at k is at k/2.
・Children of node at k are at 2k and 2k+1.
i
a[i]
0
-
1
T
2
S
3
R
S
R
4
P
5
N
6
O
7
A
P
N
O
A
8
E
9 10 11
I H G
E
I
T
1
8 E
9 I
G
T
3 R
2 S
4 P
H
5 N
10 H
6 O
7 A
11 G
Heap representations
16
Binary heap demo
Insert. Add node at end, then swim it up.
Remove the maximum. Exchange root with node at end, then sink it down.
heap ordered
T
P
R
N
H
E
I
T
O
A
G
P
R
N
H
O
A
E
I
G
17
Binary heap demo
Insert. Add node at end, then swim it up.
Remove the maximum. Exchange root with node at end, then sink it down.
heap ordered
S
R
O
N
G
P
E
H
I
S
A
R
O
N
P
G
A
E
I
H
18
Binary heap: promotion
Scenario. A key becomes larger than its parent's key.
To eliminate the violation:
・Exchange key in child with key in parent.
・Repeat until heap order restored.
S
private void swim(int k)
{
while (k > 1 && less(k/2, k))
{
exch(k, k/2);
k = k/2;
}
parent of node at k is at k/2
}
R
P
5
N
E
I
T
H
O
A
violates heap order
(larger key than parent)
G
1
T
2
5 P
N
E
R
S
I
H
O
A
G
Bottom-up reheapify (swim)
Peter principle. Node promoted to level of incompetence.
19
Binary heap: insertion
Insert. Add node at end, then swim it up.
Cost. At most 1 + lg N compares.
insert
remo
T
R
P
N
public void insert(Key x)
{
pq[++N] = x;
swim(N);
}
E
H
I
G
O
S
A
key to insert
T
R
P
N
E
H
I
G
O
add key to heap
violates heap order
S
T
swim up
R
S
N
E
A
P
I
G
O
sink down
A
H
Heap operation
20
Binary heap: demotion
Scenario. A key becomes smaller than one (or both) of its children's.
To eliminate the violation:
why not smaller child?
・Exchange key in parent with key in larger child.
・Repeat until heap order restored.
private void sink(int k)
{
children of node at k
while (2*k <= N)
are 2k and 2k+1
{
int j = 2*k;
if (j < N && less(j, j+1)) j++;
if (!less(k, j)) break;
exch(k, j);
k = j;
}
}
violates heap order
(smaller than a child)
2
R
H
5 S
P
E
T
I
N
O
A
G
T
2
5
P
E
R
S
I
10
H
N
O
A
G
Top-down reheapify (sink)
Power struggle. Better subordinate promoted.
21
Binary heap: delete the maximum
Delete max. Exchange root with node at end, then sink it down.
Cost. At most 2 lg N compares.
insert
remove the maximum
T
R
P
N
public Key delMax()
E
I
{
Key max = pq[1];
exch(1, N--);
P
sink(1);
pq[N+1] = null;
N
return max;
E
I
}
H
G
O
S
A
N
E
key to insert
P
I
G
G
R
N
add key to heap
violates heap order
S
E
P
I
G
R
P
G
O
A
remove node
from heap
T
S
S
N
violates
heap order
S
prevent loitering
O
A
A
exchange key
with root
H
T
I
O
H
R
swim up
E
R
S
T
H
key to remove
T
O
R
P
sink down
N
A
H
E
H
I
O
A
G
Heap operations
22
Binary heap: Java implementation
public class MaxPQ<Key extends Comparable<Key>>
{
private Key[] pq;
private int N;
public MaxPQ(int capacity)
{ pq = (Key[]) new Comparable[capacity+1];
}
public boolean isEmpty()
{ return N == 0; }
public void insert(Key key)
public Key delMax()
// see previous code
// see previous code
private void swim(int k)
private void sink(int k)
// see previous code
// see previous code
fixed capacity
(for simplicity)
PQ ops
private boolean less(int i, int j)
{ return pq[i].compareTo(pq[j]) < 0; }
private void exch(int i, int j)
{ Key t = pq[i]; pq[i] = pq[j]; pq[j] = t;
heap helper functions
array helper functions
}
}
23
Priority queue: implementations cost summary
implementation
insert
del max
max
unordered array
1
N
N
ordered array
N
1
1
binary heap
log N
log N
1
order-of-growth of running time for priority queue with N items
24
DELETE-RANDOM FROM A BINARY HEAP
Goal. Delete a random key from a binary heap in logarithmic time.
S
R
H
P
E
G
N
A
I
J
S
R
H
P
N
G
A
E
J
I
25
Binary heap: practical improvements
Do "half-exchanges" in sink and swim.
・Reduces number of array accesses.
・Worth doing.
Z
T
L
B
1
\
X
26
Binary heap: practical improvements
Floyd's "bounce" heuristic.
・Sink key at root all the way to bottom.
・Swim key back up.
・Overall, fewer compares; more exchanges.
・Worthwhile depending on cost of compare and exchange.
only 1 compare per node
some extra compares and exchanges
R. W. Floyd
1978 Turing award
F
X
Y
N
O
L
K
1
E
\
D
27
Binary heap: practical improvements
Caching. Binary heap is not cache friendly.
Binary heap memory layout (page size = 8 nodes)
0
1
2
3
4
8
5
9
10
6
11
12
24
32
39
40
47
7
13
14
15
31
28
Binary heap: practical improvements
Caching. Binary heap is not cache friendly.
・Cache-aligned d-heap.
・Funnel heap.
・B-heap.
・…
B-heap memory layout (page size = 8 nodes)
0
1
2
3
4
12
5
6
8
9
16
17
10
11
18
19
13
14
15
20
21
22
7
24
23
32
31
39
29
Binary heap: practical improvements
Multiway heaps.
・Complete d-way tree.
・Parent's key no smaller than its children's keys.
Fact. Height of complete d-way tree on N nodes is ~ logd N.
Z
Y
J
E
H
T
X
F
R
S
P
V
C
L
W
M
Q
N
I
G
K
A
B
D
O
3-way heap
30
Priority queues: quiz 1
How many compares (in the worst case) to insert in a d-way heap?
A.
~ log2 N
B.
~ logd N
C.
~ d log2 N
D.
~ d logd N
E.
I don't know.
31
Priority queues: quiz 2
How many compares (in the worst case) to delete-max in a d-way heap?
A.
~ log2 N
B.
~ logd N
C.
~ d log2 N
D.
~ d logd N
E.
I don't know.
32
Priority queue: implementation cost summary
implementation
insert
del max
max
unordered array
1
N
N
ordered array
N
1
1
binary heap
log N
log N
1
d-ary heap
logd N
d logd N
1
Fibonacci
1
log N †
1
Brodal queue
1
log N
1
impossible
1
1
1
sweet spot: d = 4
why impossible?
† amortized
order-of-growth of running time for priority queue with N items
33
Binary heap: considerations
Underflow and overflow.
・Underflow: throw exception if deleting from empty PQ.
・Overflow: add no-arg constructor and use resizing array.
Minimum-oriented priority queue.
leads to log N
amortized time per op
(how to make worst case?)
・Replace less() with greater().
・Implement greater().
Other operations.
・Remove an arbitrary item.
・Change the priority of an item.
can implement efficiently with sink() and swim()
[ stay tuned for Prim/Dijkstra ]
Immutability of keys.
・Assumption: client does not change keys while they're on the PQ.
・Best practice: use immutable keys.
34
Immutability: implementing in Java
Data type. Set of values and operations on those values.
Immutable data type. Can't change the data type value once created.
public final class Vector {
private final int N;
private final double[] data;
can't override instance methods
instance variables private and final
public Vector(double[] data) {
this.N = data.length;
this.data = new double[N];
for (int i = 0; i < N; i++)
this.data[i] = data[i];
}
defensive copy of mutable
instance variables
…
instance methods don't change
instance variables
}
Immutable. String, Integer, Double, Color, Vector, Transaction, Point2D.
Mutable. StringBuilder, Stack, Counter, Java array.
35
Immutability: properties
Data type. Set of values and operations on those values.
Immutable data type. Can't change the data type value once created.
Advantages.
・Simplifies debugging.
・Safer in presence of hostile code.
・Simplifies concurrent programming.
・Safe to use as key in priority queue or symbol table.
Disadvantage. Must create new object for each data type value.
“ Classes should be immutable unless there's a very good reason
to make them mutable.… If a class cannot be made immutable,
you should still limit its mutability as much as possible. ”
— Joshua Bloch (Java architect)
36
2.4 P RIORITY Q UEUES
‣ API and elementary implementations
‣ binary heaps
Algorithms
R OBERT S EDGEWICK | K EVIN W AYNE
http://algs4.cs.princeton.edu
‣ heapsort
‣ event-driven simulation
Priority queues: quiz 2
What is the name of this sorting algorithm?
public void sort(String[] a)
{
int N = a.length;
MaxPQ<String> pq = new MaxPQ<String>();
for (int i = 0; i < N; i++)
pq.insert(a[i]);
for (int i = N-1; i >= 0; i--)
a[i] = pq.delMax();
}
A.
Insertion sort.
B.
Mergesort.
C.
Quicksort.
D.
None of the above.
E.
I don't know.
38
Priority queues: quiz 3
What are its properties?
public void sort(String[] a)
{
int N = a.length;
MaxPQ<String> pq = new MaxPQ<String>();
for (int i = 0; i < N; i++)
pq.insert(a[i]);
for (int i = N-1; i >= 0; i--)
a[i] = pq.delMax();
}
A.
N log N compares in the worst case.
B.
In-place.
C.
Stable.
D.
All of the above.
E.
I don't know.
39
exch(1, 9)
sink(1, 8)
k(3, 11)
S
Heapsort
O
P
X
L
A
E
M
X tree.
T
S
S
M
R
binary
・EView input array as a complete
Heap construction: build a max-heap with all N keys.
・
exch(1, 8)
exch(1, 2)
k(2, 11)
P
S
remove
sink(1,
7) the maximum key.sink(1, 1)
・Sortdown: repeatedly
T
L
P
A
R
k(1, 211)
2
O
TT
9
XO
4 5
56
T E
8 10
1
S
9 11 10
exch(1, 7)
sink(1, 6)
S
3
R
3
R
S
E
X
7
A X
7
6
A
T
M
1
2
3
4
S5
S
O
R
T
E
O
O
T L
P M E P E E
M
R
A
E
X
T
M
A
E
SE
T
P L
L
LRE
E
M
S
sorted result
sortdownsortdown
(in place)
X
L2 E
APR
A
E E
E
O
S M T O X
E
starting point
(heap-ordered)
starting
point (heap-ordered)
M1 A
7
S8
A X
E
10
A4 L
A E5 M EO 6 O P O7 P P
S9 R T10S X11T
X
R8
R
S
X
T
result (sorted)
11
R
R
L
X
9
P
P
A
O
M
O L
E M E E X E
S
S
LR
A R
X
M
E3 E E
L
exch(1,
11)
1
2exch(1,
3
4
511)
6
7T 8
9
10 11 exch(1,
1 exch(1,
2 5)3
4 5)
5
6
T
sink(1,
10)
sink(1,
4)
sink(1,
10)
sink(1,
4)
XHeapsort:
A M P constructing
L E
X (left)
T S and
P sorting
L R A down
M O (right)
E E a heap
A E E L M O
6
P
O
X
T
exch(1,exch(1,
6)
6)
sink(1,sink(1,
5)
5)
S
P
O
X
T
L
R
XO
AP
11
E
O
result (heap-ordered)
11)
sink(5, sink(5,
11)
T
build max heap
(in place)
P M LL P E R
EA
L
starting point
(arbitrary
order) order)
starting
E point (arbitrary
PM
L
S
R
keys in arbitrary order
heap construction
heap construction
1
E
M
E
E
O
O
X
E
L
E
P
4
E
A
E
L
A in-place sort. O
R for
Basic
plan
T
exch(1, 3)
sink(1, 2)
R
E
E
A
A
R
7
L
8
9
P
R
S
A M
SR TS XT
E
MO
L 10
T
PO
11
X
E
P
X
40
Heapsort demo
Heap construction. Build max heap using bottom-up method.
we assume array entries are indexed 1 to N
array in arbitrary order
S
1
2
4
8
M
O
5
T
9
P
R
3
10
6
E
L
7
X
A
E
11
S
O
R
T
E
X
A
M
P
L
E
1
2
3
4
5
6
7
8
9
10
11
41
Heapsort demo
Sortdown. Repeatedly delete the largest remaining item.
array in sorted order
A
E
E
L
O
M
R
S
T
P
X
A
E
E
L
M
O
P
R
S
T
X
1
2
3
4
5
6
7
8
9
10
11
42
sink(5, 11)
Heapsort: heap construction
exch(1, 11)
sink(1, 10)
S
O
L
First pass. Build heap using bottom-up method.
sink(4, 11)
A
X
4
8
5
T
9
S
O
10
6
E
3
R
X
7
A
11
E
L
P
M
starting point (arbitrary order)
sink(5, 11)
E
X
E
L
E
PO
MM
OE ET EX
result (heap-ordered)
T
E
S
X
T
L
O
E
X
E M
A
L E
SM
R
T O X P
A A
AA
L
R
S
T
O
L E
SM
R
E
X
T
S
P
E
P
exch(1, 7)
sink(1, 6)
exch(1, 4)
M E
sink(1, 3)
SR
RE
M
O
R
exch(1, 8)
sink(1, 7)
exch(1, 5)
O L
sink(1, 4)
R
XS
LL
S
R
A
R
TP
A
A
SA
X
E
E
L
A
A
X
R
LE
sink(1, 11)
exch(1, 10)
sink(1, 9)
S
L
R
R
T
P
OO
R
T
P
M
M
O
M
A
X
sink(4, 11)
L
L
X
S
exch(1, 9)
8)
exch(1,sink(1,
6)
M
sink(1, 5)
P
E
E
P
M
E
E
O
M
starting point (heap-ordered)
sink(2, 11)
S
P
E
E
P
M
T
P
R
L
O
T
S
X
exch(1, 11)T
sink(1, 10)
S
O
T
X
T
E
M
E
sortdown
sink(3, 11)
2
L
O
E
heap construction
1
S
P
R
E
X
E
E
R
exch(1, 10)
sink(1, 9)
S
L
P
M
L
O
M
O
for (int k = N/2; k >= 1; k--)
T
sink(a, k, N);
A
X
E
E
P
M
P
R
T
T
E
T O X P
X
43
Heapsort: constructing (left) and sorting down
Heapsort: sortdown
sortdown
heap construction
Second pass.
1
2
S
3
R
X
7
・Remove the maximum, one at a time.
・Leave in array, instead of nulling out.
4
8
O
5
T
9
10
6
E
11
L
sink(4, 11)
exch(1, 10)
sink(1, 9)
S
L
T
E
E
P
sink(3, 11)
exch(1, 9)
sink(1, 8)
S
L
E
E
sink(2, 11)
T
L
P
E
E
O
M
sink(1, 11)
X
M
L
R
A
M
1
2
T
E
X
4
P
8
R
E
E
O
result (heap-ordered)
Heapsort: constructing (left) and sorting down (right) a heap
9
3
5
S
A
E
L
10
P
O
X
T
S
E
L
S
E
O
A
R
A
L
R
M
S
A
E
exch(1, 7)
sink(1, 6)
X
T
S
P
O
E
X
T
S
M
E
L
E
exch(1, 2)
sink(1, 1)
P
M
R
T
P
A
R
E
L
R
O
X
A
E
exch(1, 8)
sink(1, 7)
S
X
T
S
P
O
A
X
T
S
M
M
E
L
E
exch(1, 3)
sink(1, 2)
R
O
E
L
R
P
A
R
A
E
X
T
E
X
T
P
L
X
T
A
R
P
O
exch(1, 4)
sink(1, 3)
S
O
M
O
M
A
X
M
S
R
P
R
E
A
X
E
L
E
A
R
P
O
X
T
S
R
S
L
E
E
exch(1, 5)
sink(1, 4)
T
O
E
A
E
E
O
starting point (heap-ordered)
M
O
M
A
X
A
R
P
E
E
P
M
R
T
M
L
P
M
L
S
exch(1, 11)
sink(1, 10)
S
O
while (N > 1)
{
exch(a, 1, N--);
sink(a, 1, N);
}
T
A
E
L
P
M
starting point (arbitrary order)
sink(5, 11)
exch(1, 6)
sink(1, 5)
X
M
6
O
E
7
P
11
X
T
result (sorted)
44
Heapsort: Java implementation
public class Heap
{
public static void sort(Comparable[] a)
{
int N = a.length;
for (int k = N/2; k >= 1; k--)
sink(a, k, N);
while (N > 1)
{
exch(a, 1, N);
sink(a, 1, --N);
}
}
but make static (and pass arguments)
private static void sink(Comparable[] a, int k, int N)
{ /* as before */ }
private static boolean less(Comparable[] a, int i, int j)
{ /* as before */ }
private static void exch(Object[] a, int i, int j)
{ /* as before */ }
}
but convert from 1-based
indexing to 0-base indexing
45
Heapsort: trace
N
k
initial values
11
11
5
4
11
3
11
2
11
1
heap-ordered
10
9
8
7
1
1
1
1
6
1
5
1
4
1
3
1
2
1
1
1
sorted result
0
1
2
3
a[i]
4 5 6
S
S
S
O
O
O
R
R
R
T
T
T
E
L
L
X
X
X
A
A
A
M
M
M
P
P
P
L
E
E
E
E
E
S
S
X
O
T
T
X
X
S
T
P
P
L
L
L
R
R
R
A
A
A
M
M
M
P
O
O
E
E
E
E
E
E
X
T
S
R
P
T
P
P
P
O
S
S
R
E
E
P
O
O
O
M
L
L
L
L
L
R
R
E
E
E
A
A
A
A
A
M O E E
M E E X
M E T X
M S T X
R S T X
O
M
L
E
E
A
A
M
L
E
A
A
E
E
E
E
E
E
E
E
E
A
A
A
L
L
L
L
L
E
M
M
M
M
M
E
O
O
O
O
O
O
P
P
P
P
P
P
P
R
R
R
R
R
R
R
7
8
9 10 11
S
S
S
S
S
S
S
T
T
T
T
T
T
T
X
X
X
X
X
X
X
Heapsort trace (array contents just after each sink)
46
Heapsort: mathematical analysis
Proposition. Heap construction uses ≤ 2 N compares and ≤ N exchanges.
max number of exchanges
to sink node
Pf sketch. [assume N = 2h+1 – 1]
3
3
2
2
2
2
1
1
1
1
1
1
1
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
binary heap of height h = 3
a tricky sum
(see COS 340)
h + 2(h
1) + 4(h
2) + 8(h
3) + . . . + 2h (0)
h+1
22h+1
= NN
1
47
Heapsort: mathematical analysis
Proposition. Heap construction uses ≤ 2 N compares and ≤ N exchanges.
Proposition. Heapsort uses ≤ 2 N lg N compares and exchanges.
algorithm can be improved to ~ 1 N lg N
(but no such variant is known to be practical)
Significance. In-place sorting algorithm with N log N worst-case.
・Mergesort: no, linear extra space.
・Quicksort: no, quadratic time in worst case.
・Heapsort: yes!
in-place merge possible, not practical
N log N worst-case quicksort possible,
not practical
Bottom line. Heapsort is optimal for both time and space, but:
・Inner loop longer than quicksort’s.
・Makes poor use of cache.
・Not stable.
can be improved using
advanced caching tricks
48
Introsort
Goal. As fast as quicksort in practice; N log N worst case, in place.
Introsort.
・Run quicksort.
・Cutoff to heapsort if stack depth exceeds 2 lg N.
・Cutoff to insertion sort for N = 16.
In the wild. C++ STL, Microsoft .NET Framework.
49
Sorting algorithms: summary
inplace?
selection
✔
insertion
✔
shell
✔
stable?
✔
best
average
worst
remarks
½N2
½N2
½N2
N exchanges
N
¼N2
½N2
use for small N
or partially ordered
N log3 N
?
c
N 3/2
merge
✔
½ N lg N
N lg N
N lg N
timsort
✔
N
N lg N
N lg N
quick
✔
N lg N
2 N ln N
½N2
3-way quick
✔
N
2 N ln N
½N2
heap
✔
N
2 N lg N
2 N lg N
?
✔
N
N lg N
N lg N
✔
tight code;
subquadratic
N log N guarantee;
stable
improves mergesort
when preexisting order
N log N probabilistic guarantee;
fastest in practice
improves quicksort
when duplicate keys
N log N guarantee;
in-place
holy sorting grail
50
2.4 P RIORITY Q UEUES
‣ API and elementary implementations
‣ binary heaps
Algorithms
R OBERT S EDGEWICK | K EVIN W AYNE
http://algs4.cs.princeton.edu
‣ heapsort
‣ event-driven simulation
Molecular dynamics simulation of hard discs
Goal. Simulate the motion of N moving particles that behave
according to the laws of elastic collision.
52
Molecular dynamics simulation of hard discs
Goal. Simulate the motion of N moving particles that behave
according to the laws of elastic collision.
Hard disc model.
・Moving particles interact via elastic collisions with each other and walls.
・Each particle is a disc with known position, velocity, mass, and radius.
・No other forces.
temperature, pressure,
diffusion constant
motion of individual
atoms and molecules
Significance. Relates macroscopic observables to microscopic dynamics.
・Maxwell-Boltzmann: distribution of speeds as a function of temperature.
・Einstein: explain Brownian motion of pollen grains.
53
Warmup: bouncing balls
Time-driven simulation. N bouncing balls in the unit square.
public class BouncingBalls
{
public static void main(String[] args)
{
int N = Integer.parseInt(args[0]);
Ball[] balls = new Ball[N];
for (int i = 0; i < N; i++)
balls[i] = new Ball();
while(true)
{
StdDraw.clear();
for (int i = 0; i < N; i++)
{
balls[i].move(0.5);
balls[i].draw();
}
StdDraw.show(50);
}
main simulation loop
}
}
% java BouncingBalls 100
54
Warmup: bouncing balls
public class Ball
{
private double rx, ry;
private double vx, vy;
private final double radius;
public Ball(...)
{ /* initialize position and
// position
// velocity
// radius
velocity */
check for collision with walls
}
public void move(double dt)
{
if ((rx + vx*dt < radius) || (rx + vx*dt > 1.0 - radius)) { vx = -vx; }
if ((ry + vy*dt < radius) || (ry + vy*dt > 1.0 - radius)) { vy = -vy; }
rx = rx + vx*dt;
ry = ry + vy*dt;
}
public void draw()
{ StdDraw.filledCircle(rx, ry, radius); }
}
Missing. Check for balls colliding with each other.
・Physics problems: when? what effect?
・CS problems: which object does the check?
too many checks?
55
Time-driven simulation
・Discretize time in quanta of size dt.
・Update the position of each particle after every dt units of time,
and check for overlaps.
・If overlap, roll back the clock to the time of the collision, update the
velocities of the colliding particles, and continue the simulation.
t
t + dt
t + 2 dt
(collision detected)
t + Δt
(roll back clock)
56
Time-driven simulation
Main drawbacks.
・
・Simulation is too slow if dt is very small.
・May miss collisions if dt is too large.
~ N 2 / 2 overlap checks per time quantum.
dt too small: excessive computation
(if colliding particles fail to overlap when we are looking)
dt too small: excessive computation
dt too large: may miss collisions
dt too large: may miss collisions
Fundamental challenge for
time-driven simulation
57
Event-driven simulation
Change state only when something happens.
・Between collisions, particles move in straight-line trajectories.
・Focus only on times when collisions occur.
・Maintain PQ of collision events, prioritized by time.
・Remove the min = get next collision.
Collision prediction. Given position, velocity, and radius of a particle,
when will it collide next with a wall or another particle?
Collision resolution. If collision occurs, update colliding particle(s)
according to laws of elastic collisions.
prediction (at time t)
particles hit unless one passes
intersection point before the other
arrives (see Exercise 3.6.X)
resolution (at time t + dt)
velocities of both particles
change after collision (see Exercise 3.6.X)
Predicting and resolving a particle-particle collision
58
Particle-wall collision
Collision prediction and resolution.
・Particle of radius s at position (rx, ry).
・Particle moving in unit box with velocity (vx, vy).
・Will it collide with a vertical wall? If so, when?
resolution (at time t + dt)
s
velocity after collision = ( − vx , vy)
position after collision = ( 1 − s , ry + vydt)
prediction (at time t)
dt ! time to hit wall
= distance/velocity
= (1 − s − rx )/vx
wall at
x=1
(rx , ry )
vx
vy
1 − s − rx
Predicting and resolving a particle-wall collision
59
Particle-particle collision prediction
Collision prediction.
・Particle i: radius s , position (rx , ry ), velocity (vx , vy ).
・Particle j: radius s , position (rx , ry ), velocity (vx , vy ).
・Will particles i and j collide? If so, when?
i
i
i
i
i
j
j
j
j
j
(vxi', vyi')
(vxj', vyj')
mi
si
(vxi , vyi )
(rxi , ryi)
(rxi', ryi')
i
time = t
time = t + Δt
(vxj , vyj)
σj
j
60
Particle-particle collision prediction
Collision prediction.
・Particle i: radius s , position (rx , ry ), velocity (vx , vy ).
・Particle j: radius s , position (rx , ry ), velocity (vx , vy ).
・Will particles i and j collide? If so, when?
t=
i
i
i
i
i
j
j
j
j
j
B7 v · r
B7 d < 0,
v·
r +
v· v
d
0,
Qi?2`rBb2
d = (∆v · ∆r)2 − (∆v · ∆v) (∆r · ∆r − s2 ),
Δv = (Δvx, Δvy) = (vxi − vx j , vyi − vy j )
Δr = (Δrx, Δry) = (rxi − rx j , ryi − ry j )
€
€
€
€
s = si + sj
Δv ⋅ Δv = (Δvx)2 + (Δvy)2
Δr ⋅ Δr = (Δrx)2 + (Δry)2
Δv ⋅ Δr = (Δvx)(Δrx) + (Δvy)(Δry)
Important note: This is physics, so we won’t be testing you on it!
61
Particle-particle collision resolution
Collision resolution. When two particles collide, how does velocity change?
vxiiʹ′ + =Jx vx
/ mi i + Jx / mi
= vy
vyiiʹ′ + =Jy /vymi i + Jy / mi
vx jʹ′ = vx jʹ′ − =Jx vx
/ mj j− Jx / m j
vy ʹ′ = vy
vx jʹ′ − =Jy vx
/ mj − Jy / m j
vxiʹ′
vyiʹ′
=
j
j
€
€
Jx =
j
Newton's second law
(momentum form)
J
rx
s
, Jy =
J
ry
s
2 mi mj ( v · r)
, J=
s (mi + mj )
impulse due to normal force
(conservation of energy, conservation of momentum)
Important note: This is physics, so we won’t be testing you on it!
62
Particle data type skeleton
public class Particle
{
private double rx, ry;
private double vx, vy;
private final double radius;
private final double mass;
private int count;
//
//
//
//
//
position
velocity
radius
mass
number of collisions
public Particle(...) { }
public void move(double dt) { }
public void draw()
{ }
public double timeToHit(Particle that) { }
public double timeToHitVerticalWall()
{ }
public double timeToHitHorizontalWall() { }
public void bounceOff(Particle that)
public void bounceOffVerticalWall()
public void bounceOffHorizontalWall()
{ }
{ }
{ }
predict collision
with particle or wall
resolve collision
with particle or wall
}
63
Particle-particle collision and resolution implementation
public double timeToHit(Particle that)
{
if (this == that) return INFINITY;
double dx = that.rx - this.rx, dy = that.ry - this.ry;
double dvx = that.vx - this.vx; dvy = that.vy - this.vy;
double dvdr = dx*dvx + dy*dvy;
if( dvdr > 0) return INFINITY;
double dvdv = dvx*dvx + dvy*dvy;
double drdr = dx*dx + dy*dy;
double s = this.radius + that.radius;
double d = (dvdr*dvdr) - dvdv * (drdr - s*s);
if (d < 0) return INFINITY;
return -(dvdr + Math.sqrt(d)) / dvdv;
}
public void bounceOff(Particle that)
{
double dx = that.rx - this.rx, dy = that.ry
double dvx = that.vx - this.vx, dvy = that.vy
double dvdr = dx*dvx + dy*dvy;
double s = this.radius + that.radius;
double J = 2 * this.mass * that.mass * dvdr /
double Jx = J * dx / s;
double Jy = J * dy / s;
this.vx += Jx / this.mass;
this.vy += Jy / this.mass;
that.vx -= Jx / that.mass;
that.vy -= Jy / that.mass;
this.count++;
Important note: This is physics,
that.count++;
}
no collision
- this.ry;
- this.vy;
(s * (this.mass + that.mass));
so we won’t be testing you on it!
64
Collision system: event-driven simulation main loop
two particles on a collision course
Initialization.
・Fill PQ with all potential particle-wall collisions.
・Fill PQ with all potential particle-particle collisions.
third particle interferes: no collision
“potential” since collision may not happen if
some other collision intervenes
An invalidated event
Main loop.
・Delete the impending event from PQ (min priority = t).
・If the event has been invalidated, ignore it.
・Advance all particles to time t, on a straight-line trajectory.
・Update the velocities of the colliding particle(s).
・Predict future particle-wall and particle-particle collisions involving the
colliding particle(s) and insert events onto PQ.
65
Event data type
Conventions.
・Neither particle null
・One particle null
・Both particles null
⇒ particle-particle collision.
⇒ particle-wall collision.
⇒ redraw event.
private class Event implements Comparable<Event>
{
private double time;
// time of event
private Particle a, b;
// particles involved in event
private int countA, countB; // collision counts for a and b
public Event(double t, Particle a, Particle b) { ... }
create event
public int compareTo(Event that)
{ return this.time - that.time;
ordered by time
public boolean isValid() { ... }
}
}
invalid if
intervening
collision
66
Collision system implementation: skeleton
public class CollisionSystem
{
private MinPQ<Event> pq;
private double t = 0.0;
private Particle[] particles;
// the priority queue
// simulation clock time
// the array of particles
public CollisionSystem(Particle[] particles) { }
private void predict(Particle a)
add to PQ all particle-wall and particleparticle collisions involving this particle
{
if (a == null) return;
for (int i = 0; i < N; i++)
{
double dt = a.timeToHit(particles[i]);
pq.insert(new Event(t + dt, a, particles[i]));
}
pq.insert(new Event(t + a.timeToHitVerticalWall() , a, null));
pq.insert(new Event(t + a.timeToHitHorizontalWall(), null, a));
}
private void redraw()
{ }
public void simulate() {
/* see next slide */
}
}
67
Collision system implementation: main event-driven simulation loop
public void simulate()
{
pq = new MinPQ<Event>();
for(int i = 0; i < N; i++) predict(particles[i]);
pq.insert(new Event(0, null, null));
while(!pq.isEmpty())
{
Event event = pq.delMin();
if(!event.isValid()) continue;
Particle a = event.a;
Particle b = event.b;
get next event
for(int i = 0; i < N; i++)
particles[i].move(event.time - t);
t = event.time;
if
(a !=
else if (a !=
else if (a ==
else if (a ==
predict(a);
predict(b);
null
null
null
null
&&
&&
&&
&&
b
b
b
b
!=
==
!=
==
null)
null)
null)
null)
initialize PQ with
collision events and
redraw event
a.bounceOff(b);
a.bounceOffVerticalWall()
b.bounceOffHorizontalWall();
redraw();
update positions
and time
process event
predict new events
based on changes
}
}
68
Particle collision simulation: example 1
% java CollisionSystem 100
69
Particle collision simulation: example 2
% java CollisionSystem < billiards.txt
70
Particle collision simulation: example 3
% java CollisionSystem < brownian.txt
71
Particle collision simulation: example 4
% java CollisionSystem < diffusion.txt
72