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
© Copyright 2025