A12: Heaps and Event Simulation

A disease starts with one sick person and spreads from person to person along contacts, each of which takes some number of days to pass the disease on. Public health officials want to know how quickly the disease reaches everyone, and how many people would need to be vaccinated to stop it. In this assignment, you’ll write a binary heap, and then use it to simulate the spread, one event at a time, to find out.

Prerequisites: Binary heaps, sifting up and sifting down, building a heap from a list, discrete event simulation

Learning Objective: Implement a binary heap in a Python list, and see a priority queue at work as the engine of a simulation.

Background

A binary min-heap stores its items in a Python list, level by level. The children of the item at index i are at indices 2 * i + 1 and 2 * i + 2, its parent is at index (i - 1) // 2, and every item is no larger than its children, so the smallest item is at index 0. The notes on priority queues and heaps describe how to insert an item by sifting it up, how to remove the smallest item by sifting down, and how to build a heap from a list from the bottom up. They also describe discrete event simulation, which you’ll use in Part 3.

Your heap compares items with < and <=, so it can hold anything that Python can compare: numbers, strings, or tuples. Tuples are compared by their first entries, then by their second entries if the first entries are equal, and so on, so a tuple (priority, item) puts the item in the heap by its priority.

Write your heap in a file heap.py, and your simulation in a file epidemic.py. Don’t use Python’s heapq module, or sort, sorted, min, or max, in heap.py: the point is to write the heap yourself. (Your tests may use sorted.)

Part 1: The heap

Write a class MinHeap, which keeps its items in a Python list, self.items, with these methods:

class MinHeap:
    def __init__(self, items=None):
        # You write this: make an empty heap, or build a heap from a list of items
        pass

    def __len__(self):
        # You write this: the number of items, so that len(heap) works
        pass

    def is_empty(self):
        # You write this
        pass

    def peek_min(self):
        # You write this: return the smallest item, without removing it
        pass

    def insert(self, item):
        # You write this: append the item, then sift it up
        pass

    def remove_min(self):
        # You write this: remove and return the smallest item
        pass

    def _sift_up(self, i):
        # You write this: move the item at index i up until its parent is no larger
        pass

    def _sift_down(self, i):
        # You write this: move the item at index i down until neither child is smaller
        pass

peek_min and remove_min should raise an IndexError if the heap is empty:

if self.is_empty():
    raise IndexError("remove_min from an empty heap")

remove_min saves the item at index 0, to return. It then removes the last item from the list and, if the list isn’t then empty, puts that item at index 0 and sifts it down. Watch for the heap with a single item, whose last item is the one being removed. When sifting down, swap with the smaller child, and remember that a node may have two children, one child (a left child only), or none.

With no argument, the constructor makes an empty heap. Given a list of items, it copies the list into self.items, with list(items), and builds the heap from the bottom up: for each index i from len(self.items) // 2 - 1 down to 0, it sifts the item at i down. (Copying the list keeps the heap from rearranging the caller’s list.)

For example:

heap = MinHeap()
for x in [7, 3, 9, 1, 5, 2]:
    heap.insert(x)
print(heap.items)          # [1, 3, 2, 7, 5, 9]
print(heap.remove_min())   # 1
print(heap.items)          # [2, 3, 9, 7, 5]

heap = MinHeap([9, 7, 5, 3, 1, 8, 6])
print(heap.items)          # [1, 3, 5, 9, 7, 8, 6]

Testing

Write a function check_heap(heap) that checks, with assert, that every item in heap.items is no smaller than its parent. Then write a function test_heap(trials) that does the following trials times: make a list of random integers, with random.randint, of random length from 0 to 30; insert them one at a time into an empty heap, calling check_heap after each insertion; then remove them all, calling check_heap after each removal, and check that they came out in the same order as sorted puts them. Also check that MinHeap(values) gives a heap, according to check_heap, and that removing everything from it gives sorted(values). Run test_heap(1000), and check that peek_min and remove_min raise an IndexError on an empty heap.

Part 2: Tracing by hand

Do this part on paper, or in your head, before you run any code; then use your code to check your answers. Write your answers in a comment at the end of heap.py. You don’t need to draw the trees, only the lists, though you may find it helpful to sketch the tree as you go.

Build a heap from the list [38, 35, 24, 36, 30, 17, 1, 22, 14, 25], from the bottom up. Give the list after each call to _sift_down, starting with the call at index 4, the last internal node.
Starting from the heap you built, call remove_min three times. For each call, give the item returned and the list afterward.

To check your answers, add a line to the constructor that prints i and self.items after each sift, run it on the list, and then remove that line before you submit.

Part 3: Simulating an epidemic

A town has n people, numbered 0 to n - 1. Some pairs of people are in contact, and each contact has a number of days that it takes for one of the pair, if sick, to pass the disease to the other. The contacts are stored in a dictionary, contacts, which maps each person to a list of pairs (other, days). Contact is mutual, so if contacts[a] has the pair (b, days), then contacts[b] has the pair (a, days), with the same number of days. Here is a small town of six people:

contacts = {
    0: [(1, 2), (2, 5)],
    1: [(0, 2), (2, 1), (3, 4)],
    2: [(0, 5), (1, 1), (4, 2)],
    3: [(1, 4), (4, 1)],
    4: [(2, 2), (3, 1)],
    5: [],
}

Person 0 and person 2 are in contact, for example, with a delay of 5 days. Person 5 has no contacts.

When a person becomes sick, on day t, the disease is passed along every one of that person’s contacts: each contact (other, days) is exposed to the disease on day t + days. A person who is exposed while healthy becomes sick that day, and the disease is then passed along that person’s contacts in turn. A person who is exposed after already being sick is unaffected. A person who is immune never becomes sick, and passes the disease to no one.

In epidemic.py, write a function simulate(contacts, first_case, immune=()), which simulates the spread of the disease from the person first_case, who becomes sick on day 0, and returns a dictionary sick_day that maps each person who became sick to the day on which they did. People who never became sick aren’t in the dictionary. immune is a collection of immune people, and first_case isn’t one of them.

Use your MinHeap from Part 1, with from heap import MinHeap, to hold the exposures that haven’t happened yet. The priority of an exposure is the day on which it happens, so the heap gives back the exposures in order of day. Store each exposure as a tuple (day, person). At the start, insert the exposure (0, first_case). Then repeatedly remove the next exposure (day, person) from the heap, and:

If person is already in sick_day, or is immune, do nothing.
Otherwise, record sick_day[person] = day, and for each contact (other, days) of person, insert the exposure (day + days, other).

The simulation ends when the heap is empty. Notice that the code inserts an exposure for every contact of a newly sick person, even a contact who is already sick, or who will be exposed sooner by someone else. It checks whether the exposure matters only when the exposure comes out of the heap, which is when its day has arrived.

On the town above, with person 0 as the first case, the disease reaches person 1 on day 2, person 2 on day 3 (through person 1, which is faster than the direct contact with person 0), person 4 on day 5, and person 3 on day 6:

print(simulate(contacts, 0))   # {0: 0, 1: 2, 2: 3, 4: 5, 3: 6}
print(simulate(contacts, 0, immune={1}))   # {0: 0, 2: 5, 4: 7, 3: 8}

Check these two results by hand before you run your code. In the second, person 1 is immune, so the disease must go around.

Part 4: Stopping an epidemic

Write a function random_contacts(n, k) that returns a contacts dictionary for n people, in which every person is in contact with about 2 * k others: for each person i, repeat k times: choose a person j at random, with random.randrange(n), and, if j isn’t i, choose a number of days d = random.uniform(1, 3), and add the pair (j, d) to the list for i and the pair (i, d) to the list for j. Every person should have an entry in the dictionary, even a person with no contacts.

Use k = 3, and person 0 as the first case, throughout this part.

Experiment 1: how fast? Run simulate on random_contacts(n, 3), with no one immune, for n = 1000 and n = 4000. Run each a few times, since the contacts are random, and print a table of the number of people who became sick, the day on which the last person became sick, and the day on which half the town was sick.

Experiment 2: vaccination. Take n = 2000. For each fraction v in 0, 0.2, 0.4, 0.6, 0.7, 0.8, and 0.9, make people immune by choosing int(v * n) people at random, with random.sample, from everyone except person 0, and run simulate. Run each a few times, and print a table of v against the number of people who became sick.

Answer these questions in a comment at the end of epidemic.py:

How do the last day and the half-town day change when the town grows from 1000 people to 4000? Does the disease take four times as long to spread? Explain, thinking about how many people each newly sick person can reach.
In Experiment 2, how does the number of sick people change as v increases? Does vaccinating twice as many people halve the number who get sick? Roughly what fraction must be vaccinated before the disease stops spreading?
In the town in Part 3 with no one immune, the exposure of person 2 by person 0, on day 5, comes out of the heap after person 2 became sick, on day 3. What happens to it? Why can’t simulate simply avoid inserting it?
The heap holds every exposure that has been inserted but hasn’t come out yet. If everyone becomes sick, how many exposures are inserted in all, in terms of n and the number of contacts? Why is this more than n?

What to do

Implement MinHeap in heap.py, and test it with check_heap and test_heap.
Trace the heap operations of Part 2 by hand, and check your answers with your code.
Write simulate in epidemic.py, and test it on the example town.
Write random_contacts, run the experiments of Part 4, and answer the questions.

Feedback

At the top of your main Python file, add a comment that answers these questions:

About how long did the assignment take you?
What was hard?
What was easy?
Was anything in the assignment unclear?

We use your answers to improve the assignment.

Grading and requirements

MinHeap keeps its items in a Python list, self.items, and doesn’t use heapq, sort, sorted, min, or max.
insert and remove_min take $O(\log n)$ time, using _sift_up and _sift_down, and handle nodes with one child and heaps with one item correctly.
peek_min and remove_min raise an IndexError on an empty heap.
The constructor builds a heap from a list from the bottom up, in $O(n)$ time, without changing the caller’s list.
check_heap and test_heap test the heap as described in Part 1.
The hand traces of Part 2 are correct.
simulate uses a MinHeap of (day, person) exposures, ignores exposures of people who are already sick or immune, and returns the correct sick_day dictionary.
random_contacts makes mutual contacts, with the same number of days in both directions.
The experiments of Part 4 are run and printed, and the answers to the questions explain the results.

What to submit

Upload heap.py and epidemic.py to Canvas.

Challenge problems (no extra points, not required)

Fewer exposures. Keep the earliest exposure day found so far for each person, and insert an exposure only if it is earlier than that. Check that you get the same results as before, and print the largest size the heap reaches, to compare with the original.
Short infectious period. Suppose a sick person can pass the disease only to contacts whose delay is at most 4 days. Change simulate to match, and find how Experiment 1 changes.
Heapsort in place. Write heapsort(values), which sorts a list in increasing order within the list itself, with no second list: build a max-heap in the list from the bottom up, then repeatedly swap the largest value, at index 0, with the last value in the heap, shrink the heap by one, and sift down.