A data structure is a way to organize data that we wish to process with a computer program. Data structures play an essential role in computer programming—indeed, CHAPTER 4 of this book is devoted to the study of classic data structures of all sorts.

A one-dimensional array (or array) is a data structure that stores a sequence of (references to) objects. We refer to the objects within an array as its elements. The method that we use to refer to elements in an array is numbering and then indexing them. If we have n elements in the sequence, we think of them as being numbered from 0 to n – 1. Then, we can unambiguously specify one of them by referring to the i th element for any integer i in this range.

A two-dimensional array is an array of (references to) one-dimensional arrays. Whereas the elements of a one-dimensional array are indexed by a single integer, the elements of a two-dimensional array are indexed by a pair of integers: the first specifying a row, and the second specifying a column.

Often, when we have a large amount of data to process, we first put all of the data into one or more arrays. Then we use indexing to refer to individual elements and to process the data. We might have exam scores, stock prices, nucleotides in a DNA strand, or characters in a book. Each of these examples involves a large number of objects that are all of the same type. We consider such applications when we discuss data input in SECTION 1.5 and in the case study that is the subject of SECTION 1.6. In this section, we expose the basic properties of arrays by considering examples where our programs first populate arrays with objects having computed values from experimental studies and then process them.

Click here to view code image

suits = ['Clubs', 'Diamonds', 'Hearts', 'Spades']

creates an array suits[] with four strings, and the code

x = [0.30, 0.60, 0.10]
y = [0.40, 0.10, 0.50]

creates two arrays x[] and y[], each with three floats. Each array is an object that contains data (references to objects) structured for efficient access. While the truth is a bit complicated (and explained in detail in SECTION 4.1), it is useful to think of references to the elements in an array as stored contiguously, one after the other, in your computer’s memory, as shown in the diagram at right for the suits[] array defined above.

After creating an array, you can refer to any individual object anywhere you would use a variable name in a program by specifying the array name followed by an integer index within square brackets. In the preceding examples, suits[1] refers to 'Diamonds', x[0] refers to .30, y[2] refers to .50, and so forth. Note that x is a reference to the whole array, as opposed to x[i], which is a reference to the ith element. In the text, we use the notation x[] to indicate that variable x is an array (but we do not use x[] in Python code).

The obvious advantage of using an array is to avoid explicitly naming each variable individually. Using an array index is virtually the same as appending the index to the array name. For example, if we wanted to process eight floats, we could refer to them each individually with variable names like a0, a1, a2, a3, a4, a5, a6, and a7. Naming dozens of individual variables in this way would be cumbersome, however, and naming millions is untenable.

As an example of code that uses arrays, consider using arrays to represent vectors. We consider vectors in detail in SECTION 3.3; for the moment, think of a vector as a sequence of real numbers. The dot product of two vectors (of the same length) is the sum of the products of their corresponding components. For example, if our two example arrays x[] and y[] represent vectors, their dot product is the expression x[0]*y[0] + x[1]*y[1] + x[2]*y[2]. More generally, if we have two one-dimensional arrays of floats x[] and y[] whose length is given by a variable n, we can use a for loop to computer their dot product:

total = 0.0
for i in range(n)
    total += x[i]*y[i]

The simplicity of coding such computations makes the use of arrays the natural choice for all kinds of applications. Before considering more examples, we describe a number of important characteristics of programming with arrays.

a = []
for i in range(n)
    a += [0.0]

The statement a = [] creates an empty array (of length 0, with no elements) and the statement a += [0.0] appends one extra element to the end. We note that, in Python, the time required to create an array in this manner is proportional to its length (for details, see SECTION 4.1).

Often the point of a piece of code is to rearrange the elements in an array. For example, the following code reverses the order of the elements in an array a[]:

n = len(a)
for i in range(n // 2):
    temp = a[i]
    a[i] = a[n-1-i]
    a[n-1-i] = temp

The three statements in the for loop implement the exchange operation that we studied when first learning about assignment statements in SECTION 1.2. You might wish to check your understanding of arrays by studying the informal trace of this code for a seven-element array [3, 1, 4, 1, 5, 9, 2], shown at the right. This table keeps track of i and all seven array elements at the end of each iteration of the for loop.

It is natural to expect mutability in arrays, but, as you will see, mutability is a key issue in data type design and has a number of interesting implications. We will discuss some of these implications later in this section, but defer most of the discussion to SECTION 3.3.

total = 0.0
for i in range(len(a))
    total += a[i]
average = total / len(a)

Python also supports iterating over the elements in an array in a[] without referring to the indices explicitly. To do so, put the array name after the in keyword in a for statement, as follows:

total = 0.0
for v in a:
    total += v
average = total / len(a)

Python successively assigns each element in the array to the loop-control variable v, so this code is essentially equivalent to the code given earlier. In this book, we iterate over the indices when we need to refer to the array elements by their indices (as in the dot-product and array-reversal examples just considered). We iterate over the elements when only the sequence of elements matters, and not the indices (as in this compute-the-average example).

x = [.30, .60, .10]
y = x
x[1] = .99

y[1] is also .99, even though the code does not refer directly to y[1]. This situation—whenever two variables refer to the same object—is known as aliasing, and is illustrated in the object-level trace at right. We avoid aliasing arrays (and other mutable objects) because it makes it more difficult to find errors in our programs. Still, as you will see in CHAPTER 2, there are natural situations where aliasing two arrays is helpful, so it is worthwhile for you to be aware of this concept at the outset.

y = []
for v in x:
    y += [v]

After the copy operation, x and y refer to two different arrays. If we change the object that x[1] references, that change has no effect on y[1]. This situation is illustrated in the object-level trace at right.

Actually, copying an array is such a useful operation that Python provides language support for a more general operation known as slicing, where we can copy any contiguous sequence of elements within an array to another array. The expression a[i:j] evaluates to a new array whose elements are a[i], ..., a[j-1]. Moreover, the default value for i is 0 and the default value for j is len(a), so

y = x[:]

is equivalent to the code given earlier. This compact notation is convenient, but it is important for you to remember that it masks a potentially expensive operation—it takes time proportional to the length of x[].

A fundamental operation that is found in nearly every array-processing program is to create an array of n elements, each initialized to a given value. As we have seen, you can do this in Python with code like the following:

a = []
for i in range(n):
    a += [0.0]

Creating an array of a given length and initializing all of its elements to a given value is so common that Python even has a special shorthand notation for it: the code a = [0.0]*n is equivalent to the code just given. Rather than repeat such code throughout the book, we will use code like this:

a = stdarray.create1D(n, 0.0)

For consistency, stdarray also includes a create2D() function, which we will examine later in this section. We use these functions because they make our code “self-documenting” in the sense that the code itself precisely indicates our intent and does not depend on idioms specific to Python. We will delve deeply into issues surrounding library design in SECTION 2.2 (and show you how to create modules like stdarray yourself), again so that you may have a better understanding of how to approach such situations.

You now know how to create arrays in Python and access individual elements, and you have a basic understanding of Python’s built-in list data type. The table at the top of the next page describes the array operations that we have discussed. With this conceptual foundation complete, we can now move to applications. As you will see, this foundation leads to applications code that is not only easy to compose and understand, but also makes efficient use of computational resources.

Click here to view code image

SUITS = ['Clubs', 'Diamonds', 'Hearts', 'Spades']
RANKS = ['2', '3', '4', '5', '6', '7', '8', '9', '10',
         'Jack', 'Queen', 'King', 'Ace']

For example, we might use these two arrays to write a random card name, such as Queen of Clubs, as follows:

Click here to view code image

rank = random.randrange(0, len(RANKS))
suit = random.randrange(0, len(SUITS))
stdio.writeln(RANKS[rank] + ' of ' + SUITS[suit])

A more typical situation is when we compute the values to be stored in an array. For example, we might use the following code to initialize an array of length 52 that represents a deck of playing cards, using the two arrays just defined:

Click here to view code image

deck = []
for rank in RANKS:
    for suit in SUITS:
        card = rank + ' of ' + suit
        deck += [card]

After this code has been executed, writing the elements of deck[] in order from deck[0] through deck[51] with one element per line gives the sequence

2 of Clubs
2 of Diamonds
2 of Hearts
2 of Spades
3 of Clubs
3 of Diamonds
...
Ace of Hearts
Ace of Spades

temp = deck[i]
deck[i] = deck[j]
deck[j] = temp

When we use this code, we are assured that we are perhaps changing the order of the elements in the array but not the set of elements in the array. When i and j are equal, the array is unchanged. When i and j are not equal, the values a[i] and a[j] are found in different places in the array. For example, if we were to use this code with i equal to 1 and j equal to 4 in the deck[] of the previous example, it would leave the string '3 of Clubs' in deck[1] and the string '2 of Diamonds' in deck[4].

Click here to view code image

n = len(deck)
for i in range(n):
    r = random.randrange(i, n)
    temp = deck[r]
    deck[r] = deck[i]
    deck[i] = temp

Proceeding from left to right, we pick a random card from deck[i] through deck[n-1] (each card equally likely) and exchange it with deck[i]. This code is more sophisticated than it might seem. First, we ensure that the cards in the deck after the shuffle are the same as the cards in the deck before the shuffle by using the exchange idiom. Second, we ensure that the shuffle is random by choosing uniformly from the cards not yet chosen. Incidentally, Python includes a standard function named shuffle() in the random module that uniformly shuffles a given array, so the function call random.shuffle(deck) does the same job as this code.

The accompanying trace of the contents of the perm[] array at the end of each iteration of the main loop (for a run where the values of m and n are 6 and 16, respectively) illustrates the process. If the values of r are chosen such that each value in the given range is equally likely, then perm[0] through perm[m-1] are a random sample at the end of the process (even though some elements might move multiple times) because each element in the sample is chosen by taking each item not yet sampled, with each choice having equal probability of being selected.

One important reason to explicitly compute the permutation is that we can use it to write a random sample of any array by using the elements of the permutation as indices into the array. Doing so is often an attractive alternative to actually rearranging the array, because it may need to be in order for some other reason (for instance, a company might wish to draw a random sample from a list of customers that is kept in alphabetical order). To see how this trick works, suppose that we wish to draw a random poker hand from our deck[] array, constructed as just described. We use the code in sample.py with m = 5 and n = 52 and replace perm[i] with deck[perm[i]] in the stdio.write() statement (and change it to writeln()), resulting in output such as the following:

3 of Clubs
Jack of Hearts
6 of Spades
Ace of Clubs
10 of Diamonds

Sampling like this is widely used as the basis for statistical studies in polling, scientific research, and many other applications, whenever we want to draw conclusions about a large population by analyzing a small random sample. Python includes a standard function in the random module for sampling: given an array a[] and an integer k, the function call random.sample(a, k) returns a new array containing a sample of size k, chosen uniformly at random among the elements in a[].

Click here to view code image

harmonic = stdarray.create1D(n+1, 0.0)
for i in range(1, n+1):
    harmonic[i] = harmonic[i-1] + 1.0/i

Note that we waste one slot in the array (element 0) to make harmonic[1] correspond to the first harmonic number 1.0 and harmonic[i] correspond to the ith harmonic number. Precomputing values in this way is an example of a space–time tradeoff: by investing in space (to save the values), we save time (since we do not need to recompute them). This method is not effective if we need values for huge n, but it is very effective if we need values for small n many different times.

Click here to view code image

 if   m ==  1: stdio.writeln('Jan')
 elif m ==  2: stdio.writeln('Feb')
 elif m ==  3: stdio.writeln('Mar')
 elif m ==  4: stdio.writeln('Apr')
...
 elif m == 11: stdio.writeln('Nov')
 elif m == 12: stdio.writeln('Dec')

A more compact alternative is to use an array of strings holding the month names:

Click here to view code image

MONTHS = ['', 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
              'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
...
stdio.writeln(MONTHS[m])

This technique would be especially useful if you needed to access the name of a month by its number in several different places in your program. Note that, again, we intentionally waste one slot in the array (element 0) to make MONTHS[1] correspond to January, as required.

With these basic definitions and examples out of the way, we can now consider two applications that both address interesting classical problems and illustrate the fundamental importance of arrays in efficient computation. In both cases, the idea of using data to efficiently index into an array plays a central role and enables a computation that would not otherwise be feasible.

PROGRAM 1.4.2 (couponcollector.py) is an example program that simulates this process and illustrates the utility of arrays. It takes the integer n from the command line and generates a sequence of random integer values between 0 and n-1 using the function call random.randrange(0, n) (see PROGRAM 1.3.8). Each value represents a card; for each card, we want to know if we have seen that value before. To maintain that knowledge, we use an array isCollected[], which uses the card value as an index: isCollected[i] is True if we have seen a card with value i and False if we have not. When we get a new card that is represented by the integer value, we check whether we have seen its value before simply by accessing isCollected[value]. The computation consists of keeping count of the number of distinct values seen and the number of cards generated and writing the latter when the former reaches n.

As usual, the best way to understand a program is to consider a trace of the values of its variables for a typical run. It is easy to add code to couponcollector.py that produces a trace that gives the values of the variables at the end of the while loop. In the table at right, we use F for False and T for True to make the trace easier to follow. Tracing programs that use large arrays can be a challenge: when you have an array of length n in your program, it represents n variables, so you have to list them all. Tracing programs that use random.randrange() also can be a challenge because you get a different trace every time you run the program. Accordingly, we check relationships among variables carefully. Here, note that collectedCount always is equal to the number of True values in isCollected[].

Coupon collecting is no trivial problem. For example, it is very often the case that scientists want to know whether a sequence that arises in nature has the same characteristics as a random sequence. If so, that fact might be of interest; if not, further investigation may be warranted to look for patterns that might be of importance. For example, such tests are used by scientists to decide which parts of a genome are worth studying. One effective test of whether a sequence is truly random is the coupon collector test: compare the number of elements that need to be examined before all values are found against the corresponding number for a uniformly random sequence.

Without arrays, we could not contemplate simulating the coupon collector process for huge n; with arrays, it is easy to do so. We will see many examples of such processes throughout the book.

We might use a program like factors.py (PROGRAM 1.3.9) to count primes. Specifically, we could modify the code in factors.py to set a boolean variable to True if a given number is prime and False otherwise (instead of writing out factors), then enclose that code in a loop that increments a counter for each prime number. This approach is effective for small n, but becomes too slow as n grows.

PROGRAM 1.4.3 (primesieve.py) is an alternative that computes π(n) using a technique known as the sieve of Eratosthenes. The program uses an array isPrime[] of n booleans to record which of the integers less than or equal to n are prime. The goal is to set isPrime[i] to True if i is prime, and to False otherwise. The sieve works as follows: Initially, the program assigns True to all elements of the array, indicating that no factors of any integer have yet been found. Then, it repeats the following steps as long as i < n:

• Find the next smallest i for which no factors have been found.

• Leave isPrime[i] as True since i has no smaller factors.

• Assign False to all isPrime[] elements whose indices are multiples of i.

When the nested for loop ends, we have set the isPrime[] elements for all nonprimes to be False and have left the isPrime[] elements for all primes as True. Note that we might stop when i*i >= n, just as we did for factors.py, but the savings in this case would be marginal at best, since the inner for loop does not iterate at all for large i. With one more pass through the array, we can count the number of primes less than or equal to n.

As usual, it is easy to add code to write a trace. For programs such as primesieve.py, you have to be a bit careful—it contains a nested for-if-for, so you have to pay attention to the indentation to put the tracing code in the correct place.

With primesieve.py, we can compute π(n) for large n, limited primarily by the maximum array size allowed by Python. This is another example of a space–time tradeoff. Programs like primesieve.py play an important role in helping mathematicians to develop the theory of numbers, which has many important applications.

The mathematical abstraction corresponding to such tables is a matrix; the corresponding data structure is a two-dimensional array. You are likely to have already encountered many applications of matrices and two-dimensional arrays, and you will certainly encounter many others in science, in engineering, and in commercial applications, as we will demonstrate with examples throughout this book. As with vectors and one-dimensional arrays, many of the most important applications involve processing large amounts of data, and we defer considering those applications until we consider input and output, in SECTION 1.5.

Extending the one-dimensional array data structure that we have discussed to two-dimensional arrays is straightforward, once you realize that since an array can contain any type of data, its elements can also be arrays! That is, a two-dimensional array is implemented as an array of one-dimensional arrays, as detailed next.

18 19 20
21 22 23

could be represented in Python using this array of arrays:

Click here to view code image

a = [[18, 19, 20], [21, 22, 23]]

We call such an array a 2-by-3 array. By convention, the first dimension is the number of rows and the second dimension is the number of columns. Python represents a 2-by-3 array as an array that contains two objects, each of which is an array that contains three objects.

More generally, Python represents an m-by-n array as an array that contains m objects, each of which is an array that contains n objects. For example, this Python code creates an m-by-n array a[][] of floats, with all elements initialized to 0.0:

a = []
for i in range(m):
    row = [0.0] * n
    a += [row]

As for one-dimensional arrays, we use the self-descriptive alternative

stdarray.create2D(m, n, 0.0)

from our booksite module stdarray throughout this book.

for i in range(m):
    for j in range(n):
        stdio.write(a[i][j])
        stdio.write(' ')
    stdio.writeln()

This code achieves the same effect without using indices:

for row in a:
    for v in row:
        stdio.write(v)
        stdio.write(' ')
    stdio.writeln()

For example, we can add two n-by-n matrices a[][] and b[][] as follows:

Click here to view code image

c = stdarray.create2D(n, n, 0.0)
for i in range(n):
    for j in range(n):
        c[i][j] = a[i][j] + b[i][j]

Similarly, we can multiply two matrices. You may have learned matrix multiplication, but if you do not recall or are not familiar with it, the following Python code for square matrices is essentially the same as the mathematical definition. Each element c[i][j] in the product of a[][] and b[][] is computed by taking the dot product of row i of a[][] with column j of b[][].

Click here to view code image

c = stdarray.create2D(n, n, 0.0)
for i in range(n):
    for j in range(n):
        # Compute the dot product of row i and column j
        for k in range(n):
            c[i][j] += a[i][k] * b[k][j]

The definition extends to matrices that may not be square (see EXERCISE 1.4.17).

These operations provide a succinct way to express numerous matrix calculations. For example, the row-average computation for such a spreadsheet with m rows and n columns is equivalent to a matrix–vector multiplication where the row vector has n elements all equal to 1.0 / n. Similarly, the column-average computation in such a spreadsheet is equivalent to a vector-matrix multiplication where the column vector has m elements all equal to 1.0 / m. We return to vector–matrix multiplication in the context of an important application at the end of this chapter.

Click here to view code image

for i in range(len(a)):
    for j in range(len(a[i])):
        stdio.write(a[i][j])
        stdio.write(' ')
    stdio.writeln()

This code tests your understanding of Python arrays, so you should take the time to study it. In this book, we normally use square or rectangular arrays, whose dimensions are given by the variable m and n. Code that uses len(a[i]) in this way is a clear signal to you that an array is ragged.

Note that the equivalent code that does not use indices works equally well with both rectangular and ragged arrays:

for row in a:
    for v in row:
        stdio.write(v)
        stdio.write(' ')
    stdio.writeln()

Two-dimensional arrays provide a natural representation for matrices, which are omnipresent in science, mathematics, and engineering. They also provide a natural way to organize large amounts of data—a key factor in spreadsheets and many other computing applications. Through Cartesian coordinates, two- and three-dimensional arrays provide the basis for a models of the physical world. With all of these natural applications, arrays will provide fertile ground for interesting examples throughout this book as you learn to program.

The dog’s escape probability is certainly dependent on the size of the city. In a tiny 5-by-5 city, it is easy to convince yourself that the dog is certain to escape. But what are the chances of escape when the city is large? We are also interested in other parameters. For example, how long is the dog’s path, on average? How often does the dog come within one block of a previous position other than the one just left, on average? How often does the dog come within one block of escaping? These sorts of properties are important in the various applications just mentioned.

PROGRAM 1.4.4 (selfavoid.py) is a simulation of this situation that uses a two-dimensional boolean array, where each element represents an intersection. True indicates that the dog has visited the intersection; False indicates that the dog has not visited the intersection. The path starts in the center and takes random steps to places not yet visited until the dog gets stuck or escapes at a boundary. For simplicity, the code is written so that if a random choice is made to go to a spot that has already been visited, it takes no action, trusting that some subsequent random choice will find a new place (which is assured because the code explicitly tests for a dead end and exits the loop in that case).

Note that the code exhibits an important programming technique in which we code the loop-continuation test in the while statement as a guard against an illegal statement in the body of the loop. In this case, the loop-continuation test serves as a guard against an out-of-bounds array access within the loop. This corresponds to checking whether the dog has escaped. Within the loop, a successful dead-end test results in a break out of the loop.

As you can see from the sample runs, the unfortunate truth is that your dog is nearly certain to get trapped in a dead end in a large city. If you are interested in learning more about self-avoiding walks, you can find several suggestions in the exercises. For example, the dog is almost certain to escape in the three-dimensional version of the problem. While this is an intuitive result that is confirmed by our simulations, the development of a mathematical model that explains the behavior of self-avoiding walks is a famous open problem: despite extensive research, no one knows a succinct mathematical expression for the escape probability, the average length of the path, or any other important parameter.

Arrays are prominent in many of the programs that we consider, and the basic operations that we have discussed here will serve you well in addressing many programming tasks. When you are not using arrays explicitly (and you are sure to be doing so frequently), you will be using them implicitly, because all computers have a memory that is conceptually equivalent to an array.

The fundamental ingredient that arrays add to our programs is a potentially huge increase in the size of a program’s state. The state of a program can be defined as the information you need to know to understand what a program is doing. In a program without arrays, if you know the values of the variables and which statement is the next to be executed, you can normally determine what the program will do next. When we trace a program, we are essentially tracking its state. When a program uses arrays, however, there can be too huge a number of values (each of which might be changed in each statement) for us to effectively track them all. This difference makes composing programs with arrays more of a challenge than composing programs without them.

Arrays directly represent vectors and matrices, so they are of direct use in computations associated with many basic problems in science and engineering. Arrays also provide a succinct notation for manipulating a potentially huge amount of data in a uniform way, so they play a critical role in any application that involves processing large amounts of data, as you will see throughout this book.

More importantly, arrays represent the tip of an iceberg. They are your first example of a data structure, where we organize data to enable us to conveniently and efficiently process it. We will consider more examples of data structures in CHAPTER 4, including resizing arrays, linked lists, and binary search trees, and hash tables. Arrays are also your first example of a collection, which groups many elements into a single unit. We will consider many more example of collections in CHAPTER 4, including Python lists, tuples, dictionaries, and sets, along with two classic examples, stacks and queues.

A. That convention originated with machine-language programming, where the address of an array element would be computed by adding the index to the address of the beginning of an array. Starting indices at 1 would entail either a waste of space at the beginning of the array or string, or a waste of time to subtract the 1.

Q. What happens if I use a negative integer to index an array?

A. The answer may surprise you. Given an array a[], you can use the index -i as shorthand for len(a) - i. For example, you can refer to the last element in the array with a[-1] or a[len(a) - 1] and the first element with a[-len(a)] or a[0]. Python raises an IndexError at run time if you use an index outside of the range -len(a) through len(a) - 1.

Q. Why does the slice a[i:j] include a[i] but exclude a[j]?

A. The notation is consistent with ranges defined with range(), which includes the left endpoint but excludes the right endpoint. It leads to some appealing properties: j - i is the length of the subarray (assuming no truncation); a[0:len(a)] is the entire array; a[i:i] is the empty array; and a[i:j] + a[j:k] is the subarray a[i:k].

Q. What happens when I compare two arrays a[] and b[] with (a == b)?

A. It depends. For arrays (or multidimensional arrays) of numbers, it works as you might expect: the arrays are equal if each has the same length and the corresponding elements are equal.

Q. What happens when a random walk does not avoid itself?

A. This case is well understood. It is a two-dimensional version of the gambler’s ruin problem, as described in SECTION 1.3.

Q. Which pitfalls should I watch out for when using arrays?

A. Remember that creating an array takes time proportional to the length of the array. You need to be particularly careful about creating arrays within loops.

1.4.2 Given two vectors of length n that are represented with one-dimensional arrays, compose a code fragment that computes the Euclidean distance between them (the square root of the sums of the squares of the differences between corresponding elements).

1.4.3 Compose a code fragment that reverses the order of the elements in a one-dimensional array of floats. Do not create another array to hold the result. Hint: Use the code in the text for exchanging two elements.

1.4.4 What is wrong with the following code fragment?

a = []
for i in range(10):
    a[i] = i * i

Solution: Initially a is the empty array. Subsequently no elements are appended to the array. Thus a[0], a[1], and so forth do not exist. The attempts to use them in an assignment statement will raise an IndexError at run time.

1.4.5 Compose a code fragment that writes the contents of a two-dimensional array of booleans, using * to represent True and a space to represent False. Include row and column numbers.

1.4.6 What does the following code fragment write?

a = stdarray.create1D(10, 0)
for i in range(10):
    a[i] = 9 - i
for i in range(10):
    a[i] = a[a[i]]
for v in a:
    stdio.writeln(v)

1.4.7 What is a[] after executing the following code fragment?

n = 10
a = [0, 1]
for i in range(2, n):
    a += [a[i-1] + a[i-2]]

1.4.8 Compose a program deal.py that takes a command-line argument n and writes n poker hands (five cards each) from a shuffled deck, separated by blank lines.

1.4.9 Compose code fragments to create a two-dimensional array b[][] that is a copy of an existing two-dimensional array a[][], under each of the following assumptions:

a. a is square

b. a is rectangular

c. a may be ragged

Your solution to b should work for a, and your solution to c should work for both b and a.

1.4.10 Compose a code fragment to write the transposition (rows and columns changed) of a two-dimensional array. For the example spreadsheet array in the text, your code would write the following:

Click here to view code image

99 98 92 94 99 90 76 92 97 89
85 57 77 32 34 46 59 66 71 29
98 78 76 11 22 54 88 89 24 38

1.4.11 Compose a code fragment to transpose a square two-dimensional array in place without creating a second array.

1.4.12 Compose a code fragment to create a two-dimensional array b[][] that is the transpose of an existing m-by-n array a[][].

1.4.13 Compose a program that computes the product of two square matrices of booleans, using the or operator instead of + and the and operator instead of *.

1.4.14 Compose a program that accepts an integer n from the command line and creates an n-by-n boolean array a such that a[i][j] is True if i and j are relatively prime (have no common factors), and False otherwise. Then write the array (see EXERCISE 1.4.5) using * to represent True and a space to represent False. Include row and column numbers. Hint: Use sieving.

1.4.15 Modify the spreadsheet code fragment in the text to compute a weighted average of the rows, where the weights of each test score are in a one-dimensional array weights[]. For example, to assign the last of the three tests in our example to be twice the weight of the others, you would use

weights = [.25, .25, .50]

Note that the weights should sum to 1.0.

1.4.16 Compose a code fragment to multiply two rectangular matrices that are not necessarily square. Note: For the dot product to be well defined, the number of columns in the first matrix must be equal to the number of rows in the second matrix. Write an error message if the dimensions do not satisfy this condition.

1.4.17 Modify selfavoid.py (PROGRAM 1.4.4) to calculate and write the average length of the paths as well as the dead-end probability. Keep separate the average lengths of escape paths and dead-end paths.

1.4.18 Modify selfavoid.py to calculate and write the average area of the smallest axis-oriented rectangle that encloses the path. Keep separate statistics for escape paths and dead-end paths.

1.4.19 Compose a code fragment that creates a three-dimensional n-by-n-by-n array of booleans, with each element initialized to False.

Click here to view code image

probabilities = stdarray.create1D(13, 0.0)
for i in range(1, 7):
    for j in range(1, 7):
        probabilities[i+j] += 1.0
for k in range(2, 13):
    probabilities[k] /= 36.0

After this code completes, probabilities[k] is the probability that the dice sum to k. Run experiments to validate this calculation simulating n dice throws, keeping track of the frequencies of occurrence of each value when you compute the sum of two random integers between 1 and 6. How large does n have to be before your empirical results match the exact results to three decimal places?

1.4.21 Longest plateau. Given an array of integers, compose a program that finds the length and location of the longest contiguous sequence of equal values where the values of the elements just before and just after this sequence are smaller.

1.4.22 Empirical shuffle check. Run computational experiments to check that our shuffling code works as advertised. Compose a program shuffletest.py that takes command-line arguments m and n, does n shuffles of an array of size m that is initialized with a[i] = i before each shuffle, and writes an m-by-m table such that row i gives the number of times i wound up in position j for all j. All elements in the array should be close to n/m.

1.4.23 Bad shuffling. Suppose that you choose a random integer between 0 and n-1 in our shuffling code instead of one between i and n-1. Show that the resulting order is not equally likely to be one of the n! possibilities. Run the test of the previous exercise for this version.

1.4.24 Music shuffling. You set your music player to shuffle mode. It plays each of the m songs before repeating any. Compose a program to estimate the likelihood that you will not hear any sequential pair of songs (that is, song 3 does not follow song 2, song 10 does not follow song 9, and so on).

1.4.25 Minima in permutations. Compose a program that takes an integer n from the command line, generates a random permutation, writes the permutation, and writes the number of left-to-right minima in the permutation (the number of times an element is the smallest seen so far). Then compose a program that takes integers m and n from the command line, generates m random permutations of size n, and writes the average number of left-to-right minima in the permutations generated. Extra credit: Formulate a hypothesis about the number of left-to-right minima in a permutation of size n, as a function of n.

1.4.26 Inverse permutation. Compose a program that reads in a permutation of the integers 0 to n-1 from n command-line arguments and writes its inverse. (If the permutation is an array a[], its inverse is the array b[] such that a[b[i]] = b[a[i]] = i.) Be sure to check that the input is a valid permutation.

1.4.27 Hadamard matrix. The n-by-n Hadamard matrix H_n is a boolean matrix with the remarkable property that any two rows differ in exactly n/2 elements. (This property makes it useful for designing error-correcting codes.) H₁ is a 1-by-1 matrix with the single element True, and for n > 1, H_2n is obtained by aligning four copies of H_n in a large square, and then inverting all of the elements in the lower right n-by-n copy, as shown in the following examples (with T representing True and F representing False, as usual).

Compose a program that takes one command-line argument n and writes H_n. Assume that n is a power of 2.

1.4.28 Rumors. Alice is throwing a party with n other guests, including Bob. Bob starts a rumor about Alice by telling it to one of the other guests. A person hearing this rumor for the first time will immediately tell it to one other guest, chosen at random from all the people at the party except Alice and the person from whom they heard it. If a person (including Bob) hears the rumor for a second time, he or she will not propagate it further. Compose a program to estimate the probability that everyone at the party (except Alice) will hear the rumor before it stops propagating. Also calculate an estimate of the expected number of people to hear the rumor.

1.4.29 Find a duplicate. Given an array of n elements with each element between 1 and n, compose a code fragment to determine whether there are any duplicates. You do not need to preserve the contents of the given array, but do not use an extra array.

1.4.30 Counting primes. Compare primesieve.py with the alternative approach described in the text. This is an example of a space–time tradeoff: primesieve.py is fast, but requires a boolean array of length n; the alternative approach uses only two integer variables, but is substantially slower. Estimate the magnitude of this difference by finding the value of n for which this second approach can complete the computation in about the same time as python primesieve.py 1000000.

1.4.31 Minesweeper. Compose a program that takes three command-line arguments m, n, and p and produces an m-by-n boolean array where each element is occupied with probability p. In the minesweeper game, occupied cells represent bombs and empty cells represent safe cells. Write the array using an asterisk for bombs and a period for safe cells. Then, replace each safe square with the number of neighboring bombs (above, below, left, right, or diagonal) and write the result, as in this example.

* * . . .        * * 1 0 0
. . . . .        3 3 2 0 0
. * . . .        1 * 1 0 0

Try to express your code so that you have as few special cases as possible to deal with, by using an (m+2)-by-(n+2) boolean array.

1.4.32 Self-avoiding walk length. Suppose that there is no limit on the size of the grid. Run experiments to estimate the average walk length.

1.4.33 Three-dimensional self-avoiding walks. Run experiments to verify that the dead-end probability is 0 for a three-dimensional self-avoiding walk and to compute the average walk length for various values of n.

1.4.34 Random walkers. Suppose that n random walkers, starting in the center of an n-by-n grid, move one step at a time, choosing to go left, right, up, or down with equal probability at each step. Compose a program to help formulate and test a hypothesis about the number of steps taken before all cells are touched.

1.4.35 Bridge hands. In the game of bridge, four players are dealt hands of 13 cards each. An important statistic is the distribution of the number of cards in each suit in a hand. Which is the most likely, 5-3-3-2, 4-4-3-2, or 4-3-3-3? Compose a program to help you answer this question.

1.4.36 Birthday problem. Suppose that people continue to enter an empty room until a pair of people share a birthday. On average, how many people will have to enter before there is a match? Run experiments to estimate the value of this quantity. Assume birthdays to be uniform random integers between 0 and 364.

1.4.37 Coupon collector. Run experiments to validate the classical mathematical result that the expected number of coupons needed to collect n values is about nH_n. For example, if you are observing the cards carefully at the blackjack table (and the dealer has enough decks randomly shuffled together), you will wait until about 235 cards are dealt, on average, before seeing every card value.

1.4.38 Riffle shuffle. Compose a program to rearrange a deck of n cards using the Gilbert–Shannon–Reeds model of a riffle shuffle. First, generate a random integer r according to a binomial distribution: flip a fair coin n times and let r be the number of heads. Now, divide the deck into two piles: the first r cards and the remaining n – r cards. To complete the shuffle, repeatedly take the top card from one of the two piles and put it on the bottom of a new pile. If there are n₁ cards remaining in the first pile and n₂ cards remaining in the second pile, choose the next card from the first pile with probability n₁ / (n₁ + n₂) and from the second pile with probability n₂ / (n₁ + n₂). Investigate how many riffle shuffles you need to apply to a deck of 52 cards to produce a (nearly) uniformly shuffled deck.

1.4.39 Binomial coefficients. Compose a program that builds and writes a two-dimensional ragged array a such that a[n][k] contains the probability that you get exactly k heads when you toss a fair coin n times. Take a command-line argument to specify the maximum value of n. These numbers are known as the binomial distribution: if you multiply each element in row k by 2ⁿ, you get the binomial coefficients (the coefficients of x^k in (x+1)ⁿ) arranged in Pascal’s triangle. To compute them, start with a[n][0] = 0.0 for all n and a[1][1] = 1.0, then compute values in successive rows, left to right, with a[n][k] = (a[n-1][k] + a[n-1][k-1])/2.0.

The abbreviation I/O is universally understood to mean input/output, a collective term that refers to the mechanisms by which programs communicate with the outside world. Your computer’s operating system controls the physical devices that are connected to your computer. To implement our “standard I/O” abstractions, we use modules containing functions that interface to the operating system.

You have already been accepting arguments from the command line and writing strings in a terminal window; the purpose of this section is to provide you with a much richer set of tools for processing and presenting data. Like the stdio.write() and stdio.writeln() functions that you have been using, these functions do not implement pure mathematical functions—their purpose is to produce some side effect, either on an input device or on an output device. Our prime concern is using such devices to get information into and out of our programs.

An essential feature of standard I/O mechanisms is that there is no limit on the amount of input or output data, from the point of view of the program. Your programs can consume input or produce output indefinitely.

One use of standard I/O mechanisms is to connect your programs to files on your computer’s external storage. It is easy to connect standard input, standard output, standard drawing, and standard audio to files. Such connections make it easy to have your Python programs save or load results to files for archival purposes or for later reference by other programs or other applications.

A Python program takes input values from the command line and writes a string of characters as output. By default, both command-line arguments and standard output are associated with the application that takes commands (the one in which you have been typing the python command). We use the generic term terminal window to refer to this application. This model has proved to be a convenient and direct way for us to interact with our programs and data.

For reference, and as a starting point, randomseq.py (PROGRAM 1.5.1) is a program that uses this model. It takes a command-line argument n and produces an output sequence of n random numbers between 0 and 1.

Now we will complement command-line arguments and standard output with three additional mechanisms that address their limitations and provide us with a far more useful programming model. These mechanisms give us a new bird’s-eye view of a Python program in which the program converts a standard input stream and a sequence of command-line arguments into a standard output stream, a standard drawing, and a standard audio stream.

To use these modules, you must make the files stdio.py, stddraw.py, and stdaudio.py available to Python (see the Q&A at the end of this section).

The standard input and standard output abstractions date back to the development of the Unix operating system in the 1970s and are found in some form on all modern systems. Although they are primitive by comparison to various mechanisms developed since, modern programmers still depend on them as a reliable way to connect data to programs. We have developed for this book stddraw and stdaudio in the same spirit as these earlier abstractions to provide you with an easy way to produce visual and aural output.

Since the first time that we wrote float objects, we have been distracted by excessive precision in the output. For example, when we call stdio.write(math.pi), we get the output 3.141592653589793, even though we might prefer to see 3.14 or 3.14159. The stdio.write() and stdio.writeln() functions present each number to up to 16 decimal digits of precision even when we would be happy with just a few digits of precision. The stdio.writef() function is more flexible: it allows us to specify the number of digits and the precision when converting numeric objects to strings for output. With stdio.writef(), we can write stdio.writef('%7.5f', math.pi) to get 3.14159.

Next, we describe the meaning and operation of these statements, along with extensions to handle the other built-in types of data.

• w is the field width, the number of characters that should be written. If the number of characters to be written exceeds (or equals) the field width, then the field width is ignored; otherwise, the output is padded with spaces on the left. A negative field width indicates that the output instead should be padded with spaces on the right.

• p is the precision. For floats, the precision is the number of digits that should be written after the decimal point; for strings, it is the number of characters of the string that should be written. The precision is not used with integers.

• c is the conversion code. It should be d when writing an integer, f when writing a float, e when writing a float using scientific notation, and s when writing a string.

The field width and precision can be omitted, but every specification must have a conversion code.

Python must be able to convert the second argument to the type required by the specification. For s, there is no restriction, because every type of data can be converted to a string (via a call to the str() function). In contrast, a statement such as stdio.writef('%12d', 'Hello')—which asks Python to convert a string to an integer—causes Python to raise a TypeError at run time. The table at the bottom of this page shows format strings containing some common conversion specifications.

Any part of the format string that is not a conversion specification is simply passed through to standard output. For example, the statement

Click here to view code image

stdio.writef('pi is approximately %.2f ', math.pi)

writes the line

pi is approximately 3.14

Note that we need to explicitly include the newline character in the argument to write a new line with stdio.writef().

Click here to view code image

stdio.writef('The square root of %.1f is %.6f', c, t)

to get output like

Click here to view code image

The square root of 2.0 is 1.414214

As a more detailed example, if you were making payments on a loan, you might use code whose inner loop contains the statements

Click here to view code image

format = '%3s  $%6.2f   $%7.2f   $%5.2f '
stdio.writef(format, month[i], pay, balance, interest)

to write the second and subsequent lines in a table like this (see EXERCISE 1.5.14):

Click here to view code image

    payment    balance  interest
Jan $299.00   $9742.67   $41.67
Feb $299.00   $9484.26   $40.59
Mar $299.00   $9224.78   $39.52
 ...

Formatted writing is convenient because this sort of code is much more compact than the string-concatenation approach that we have been using to create output strings. We have described only the basic options; see the booksite for many details.

Click here to view code image

import sys
import stdio
n = int(sys.argv[1])
total = 0
for i in range(n):
    total += stdio.readInt()
stdio.writeln('Sum is ' + str(total))

When you type python addints.py 4, the program starts execution. It takes the command-line argument, initializes total to 0, enters the for loop, eventually calls stdio.readInt(), and waits for you to type an integer. Suppose that you want 144 to be the first value. As you type 1, then 4, and then 4, nothing happens, because stdio does not know that you are done typing the integer. But when you finally type <return> to signify the end of your integer, stdio.readInt() immediately returns the value 144, which your program adds to total and then calls stdio.readInt() again. Again, nothing happens until you type the second value: if you type 2, then 3, then 3, and then <return> to end the number, stdio.readInt() returns the value 233, which your program again adds to total. After you have typed four numbers in this way, the program expects no more input and writes the sum, as desired. In the command-line traces, we use boldface to highlight the text that you type and differentiate it from the output of the program.

How do we signal that we have no more data to type? By convention, we type a special sequence of characters known as the end-of-file sequence. Sadly, the terminal applications that we typically encounter on modern operating systems use different conventions for this critically important sequence. In this book, we use <Ctrl-d> (many systems require <Ctrl-d> to appear on a line by itself); the other widely used convention is <Ctrl-z> on a line by itself.

Actually, we rarely type numbers one by one on standard input. Instead, we keep our input data in files, as illustrated in the example accompanying PROGRAM 1.5.3 and explained in detail in the text.

Certainly average.py is a simple program, but it represents a profound new capability in programming: with standard input, we can compose programs that process an unlimited amount of data. As you will see, composing such programs is an effective approach for numerous data-processing applications.

Standard input is a substantial step up from the command-line argument model that we have been using, for two reasons, as illustrated by twentyquestions.py and average.py. First, we can interact with our program—with command-line arguments, we can provide data to the program only before it begins execution. Second, we can read in large amounts of data—with command-line arguments, we can enter only values that fit on the command line. Indeed, as illustrated by average.py, the amount of data processed by a program can be potentially unlimited, and many programs are made simpler by that assumption. A third reason for standard input is that your operating system makes it possible to change the source of standard input, so that you do not have to type all the input. Next, we consider the mechanisms that enable this possibility.

Click here to view code image

% python randomseq.py 1000 > data.txt

specifies that the standard output stream is not to be written in the terminal window, but instead is to be written to a text file named data.txt. Each call to stdio.write(), stdio.writeln(), or stdio.writef() appends text at the end of that file. In this example, the end result is a file that contains 1,000 random values. No output appears in the terminal window: it goes directly into the file named after the > symbol. Thus, we can save information for later retrieval. Note that we do not have to change randomseq.py (PROGRAM 1.5.1) in any way for this mechanism to work—it relies on using the standard output abstraction and is unaffected by our use of a different implementation of that abstraction. You can use this mechanism to save output from any program that you compose. Once we have expended a significant amount of effort to obtain a result, we often want to save the result for later reference. In a modern system, you can save some information by using cut-and-paste or some similar mechanism that is provided by the operating system, but cut-and-paste is inconvenient for large amounts of data. By contrast, redirection is specifically designed to make it easy to handle large amounts of data.

Click here to view code image

% python average.py < data.txt

This command reads a sequence of numbers from the file data.txt and computes their average value. Specifically, the < symbol is a directive that tells the operating system to implement the standard input stream by reading from the text file data.txt instead of waiting for the user to type something into the terminal window. When the program calls stdio.readFloat(), the operating system reads the value from the file. The file data.txt could have been created by any application, not just a Python program—almost every application on your computer can create text files. This facility to redirect from a file to standard input enables us to create data-driven code, in which we can change the data processed by a program without having to change the program at all. Instead, we keep data in files and compose programs that read from standard input.

Click here to view code image

% python randomseq.py 1000 | python average.py

specifies that the standard output for randomseq.py and the standard input stream for average.py are the same stream. The effect is as if randomseq.py were typing the numbers it generates into the terminal window while average.py is running. This example also has the same effect as the following sequence of commands:

Click here to view code image

% python randomseq.py 1000 > data.txt
% python average.py < data.txt

With piping, however, the file data.txt is not created. This difference is profound, because it removes another limitation on the size of the input and output streams that we can process. For example, we could replace 1000 in our example with 1000000000, even though we might not have the space to save a billion numbers on our computer (we do need the time to process them, however). When randomseq.py calls stdio.writeln(), a string is added to the end of the stream; when average.py calls stdio.readFloat(), a string is removed from the beginning of the stream. The timing of precisely what happens is up to the operating system: it might run randomseq.py until it produces some output, and then run average.py to consume that output, or it might run average.py until it needs some output, and then run randomseq.py until it produces the needed output. The end result is the same, but our programs are freed from worrying about such details because they work solely with the standard input and standard output abstractions.

For many common tasks, it is convenient to think of each program as a filter that converts a standard input stream to a standard output stream in some way, with piping as the command mechanism to connect programs together. For example, rangefilter.py (PROGRAM 1.5.4) takes two command-line arguments and writes to standard output those numbers from standard input that fall within the specified range. You might imagine standard input to be measurement data from some instrument, with the filter being used to throw away data outside the range of interest for the experiment at hand.

Several standard filters that were designed for Unix still survive (sometimes with different names) as commands in modern operating systems. For example, the sort filter reads the lines from standard input and writes them to standard output in sorted order:

Click here to view code image

% python randomseq.py 9 | sort
0.0472650078535
0.0681950168757
0.0967410236589
0.0974385525393
0.118855769243
0.46604926859
0.522853708616
0.599692836211
0.685576779833

We discuss sorting in SECTION 4.2. A second useful filter is grep, which writes the lines from standard input that match a given pattern. For example, if you type

% grep lo < rangefilter.py

you get all the lines in rangefilter.py that contain 'lo':

Click here to view code image

lo = int(sys.argv[1])
    if (value >= lo) and (value <= hi):

Programmers often use tools such as grep to get a quick reminder of variable names or language usage details. A third useful filter is more, which reads data from standard input (or from a file specified as a command-line argument) and displays it in your terminal window one screenful at a time. For example, if you type

Click here to view code image

% python randomseq.py 1000 | more

you will see as many numbers as fit in your terminal window, but more will wait for you to hit the space bar before displaying each succeeding screenful. The term filter is perhaps misleading: it was meant to describe programs like rangefilter.py that write some subsequence of standard input to standard output, but it is now often used to describe any program that reads from standard input and writes to standard output.

Processing large amounts of information plays an essential role in many applications of computing. A scientist may need to analyze data collected from a series of experiments, a stock trader may wish to analyze information about recent financial transactions, or a student may wish to maintain collections of music and movies. In these and countless other applications, data-driven programs are the norm. Standard output, standard input, redirection, and piping provide us with the capability to address such applications with our Python programs. We can collect data into files on our computer through the web or any of the standard devices and use redirection and piping to connect data to our programs.

Standard drawing is very simple: we imagine an abstract drawing device capable of drawing lines and points on a two-dimensional canvas and then displaying that “canvas” on your screen in the standard drawing window. The device is capable of responding to the commands that our programs issue in the form of calls to functions in the stddraw module.

The module’s API consists of two kinds of functions: drawing functions that cause the device to take an action (such as drawing a line or drawing a point) and control functions that control how the drawing is shown and set parameters such as the pen size or the coordinate scales.

The control function stddraw.show() needs a bit more explanation. When your program calls any drawing function such as stddraw.line() or stddraw.point(), stddraw uses an abstraction known as the background canvas. The background canvas is not displayed; it exists only in computer memory. All points, lines, and so forth are drawn on the background canvas, not directly in the standard drawing window. Only when you call stddraw.show() does your drawing get copied from the background canvas to the standard drawing window, where it is displayed until the user closes the standard drawing window—typically by clicking on the Close button in the window’s title bar.

Why does stddraw need to use a background canvas? The main reason is that use of two canvases instead of one makes the stddraw module more efficient. Incrementally displaying a complex drawing as it is being created can be intolerably inefficient on many computer systems. In computer graphics, this technique is known as double buffering.

To summarize the information that you need to know, a typical program using the stddraw module has this structure:

• Import the stddraw module.

• Call drawing functions such as stddraw.line() and stddraw.point() to create a drawing on the background canvas.

• Call stddraw.show() to show the background canvas in the standard drawing window and wait until the window is closed.

It is important to remember that all drawing goes to the background canvas. Typically, a program that creates a drawing finishes with a call to stddraw.show(), because only when stddraw.show() is called will you see your picture.

Next, we will consider several examples that will open up a whole new world of programming by freeing you from the restriction of communicating with your program just through text.

When you use a computer to create drawings, you get immediate feedback (the drawing) so that you can refine and improve your program quickly. With a computer program, you can create drawings that you could not contemplate making by hand. In particular, instead of viewing our data as just numbers, we can use pictures, which are far more expressive. We will consider other graphics examples after we discuss a few other drawing commands.

For example, when you call the function stddraw.setXscale(0, n), you are telling the drawing device that you will be using x-coordinates between 0 and n. Note that the two-call sequence

stddraw.setXscale(x0, x1)
stddraw.setYscale(y0, y1)

sets the drawing coordinates to be within a bounding box whose lower-left corner is at (x₀, y₀) and whose upper-right corner is at (x₁, y₁). If you use integer coordinates, Python promotes them to floats, as expected. The figure at right is a simple example that demonstrates the utility of scaling. Scaling is the simplest of the transformations commonly used in graphics. Several of the applications that we consider in this chapter are typical—we use scaling in a straightforward way to match our drawings to our data.

The pen is circular, so that when you set the pen radius to r and draw a point, you get a circle of radius r. Also, lines are of thickness 2r and have rounded ends. The default pen radius is 0.005 and is not affected by coordinate scaling. This default pen radius is about 1/200 the width of the default window, so that if you draw 100 points equally spaced along a horizontal or vertical line, you will be able to see individual circles, but if you draw 200 such points, the result will look like a line. When you make the function call stddraw.setPenRadius(.025), you are saying that you want the thickness of the lines and the size of the points to be five times the 0.005 standard. To draw points with the minimum possible radius (one pixel on typical displays), set the pen radius to 0.0.

The graphical representation of points plotted in this way is far more expressive (and far more compact) than the numbers themselves or anything that we could create with the standard output representation that our programs have been limited to until now. The plotted image that is produced by plotfilter.py makes it far easier for us to infer properties of the cities (such as, for example, clustering of population centers) than does a list of the coordinates. Whenever we are processing data that represents the physical world, a visual image is likely to be one of the most meaningful ways in which we can use to display output. PROGRAM 1.5.5 illustrates just how easily you can create such an image.

The smoothness of the curve depends on properties of the function and the size of the sample. If the sample size is too small, the rendition of the function may not be at all accurate (it might not be very smooth, and it might miss major fluctuations, as shown in the example on the next page); if the sample is too large, producing the plot may be time-consuming, since some functions are time-consuming to compute. (In SECTION 2.4, we will look at an efficient method for accurately plotting a smooth curve.) You can use this same technique to plot the function graph of any function you choose: identify an x-interval where you want to plot the function, compute function values evenly spaced in that interval and store them in an array, determine and set the y-scale, and draw the line segments.

The arguments for stddraw.circle() and stddraw.filledCircle() define a circle of radius r centered at (x, y); the arguments for stddraw.square() and stddraw.filledSquare() define a square of side length 2r centered at (x, y); and the arguments for stddraw.polygon() and stddraw.filledPolygon() define a sequence of points that we connect by lines, including one from the last point to the first point. If you want to define shapes other than squares or circles, use one of these functions. To check your understanding, try to figure out what this code draws, before reading the answer:

xd = [x-r, x, x+r, x]
yd = [y, y+r, y, y-r]
stddraw.polygon(xd, yd)

The answer is that you would never know, because it draws on the background canvas and there is no call to stddraw.show(). If there were such a call, it would draw a diamond (a rotated square) centered at the point (x, y). Several other examples of code that draws shapes and filled shapes are shown on the facing page.

In this code, color and fonts use types that you will learn about in SECTION 3.1. Until then, we leave the details to stddraw. The available pen colors are BLACK, BLUE, CYAN, DARK_GRAY, GRAY, GREEN, LIGHT_GRAY, MAGENTA, ORANGE, PINK, RED, WHITE, and YELLOW, defined as constants within stddraw. For example, the call stddraw.setPenColor(stddraw.GRAY) changes to gray ink. The default ink color is stddraw.BLACK. The default font in stddraw suffices for most of the drawings that you need (and you can find information on using other fonts on the booksite). For example, you might wish to call these functions to annotate function plots to highlight relevant values, and you might find it useful to compose similar functions to annotate other parts of your drawings.

Shapes, color, and text are basic tools that you can use to produce a dizzying variety of images, but you should use them judiciously. Use of such artifacts usually presents a design challenge, and our stddraw commands are crude by the standards of modern graphics libraries, so that you are likely to need more code than you might expect to produce the beautiful images that you may imagine.

• Clear the background canvas.

• Draw a black ball at the new position.

• Show the drawing and wait for a short while.

To create the illusion of movement, we iterate these steps for a whole sequence of positions (one that will form a straight line, in this case). The argument to stddraw.show() quantifies “a short while” and controls the apparent speed.

PROGRAM 1.5.7 (bouncingball.py) implements these steps to create the illusion of a ball moving in the 2-by-2 box centered at the origin. The current position of the ball is (r_x, r_y), and we compute the new position at each step by adding v_x to r_x and v_y to r_y. Since (v_x, v_y) is the fixed distance that the ball moves in each time unit, it represents the velocity. To keep the ball in the drawing, we simulate the effect of the ball bouncing off the walls according to the laws of elastic collision. This effect is easy to implement: when the ball hits a vertical wall, we just change the velocity in the x-direction from v_x to –v_x, and when the ball hits a horizontal wall, we change the velocity in the y-direction from v_y to –v_y. Of course, you have to download the code from the booksite and run it on your computer to see motion.

Since we cannot show a moving image on the printed page, we slightly modified bouncingball.py to show the track of the ball as it moves (see EXERCISE 1.5.34) in the examples shown below the code.

To familiarize yourself with animation on your computer, you are encouraged to modify the various parameters in bouncingball.py to draw a bigger ball, make it move faster or slower, and experiment with the distinction between the velocity in the simulation and the apparent speed on your display. For maximum flexibility, you might wish to modify bouncingball.py to take all these parameters as command-line arguments.

Standard drawing substantially improves our programming model by adding a “picture is worth a thousand words” component. It is a natural abstraction that you can use to better open up your programs to the outside world. With it, you can easily produce the function plots and visual representations of data that are commonly used in science and engineering. We will put it to such uses frequently throughout this book. Any time that you spend now working with the sample programs on the last few pages will be well worth the investment. You can find many useful examples on the booksite and in the exercises, and you are certain to find some outlet for your creativity by using stddraw to meet various challenges. Can you draw an n-pointed star? Can you make our bouncing ball actually bounce (add gravity)? You may be surprised at how easily you can accomplish these and other tasks.

For example, suppose that you want to play concert A for 10 seconds. At 44,100 samples per second, you need an array of 441,001 float values. To fill in the array, use a for loop that samples the function sin(2πt × 440) at t = 0 / 44100, 1 / 44100, 2 / 44100, 3 / 44100, ..., 441000 / 44100. Once we fill the array with these values, we are ready for stdaudio.playSamples(), as in the following code:

Click here to view code image

SPS = 44100              # samples per second
hz = 440.0               # concert A
duration = 10.0          # ten seconds
n = int(SPS * duration)

a = stdarray.create1D(n+1)
for i in range(n+1):
    a[i] = math.sin(2.0 * math.pi * i * hz / SPS)
stdaudio.playSamples(a)
stdaudio.wait()

This code is the “Hello, World” of digital audio. Once you use it to get your computer to play this note, you can compose code to play other notes and make music! The difference between creating sound and plotting an oscillating curve is nothing more than the output device. Indeed, it is instructive and entertaining to send the same numbers to both standard drawing and standard audio (see EXERCISE 1.5.27).

To see how easily we can create music with the functions in the stdaudio module (detailed at the top of page 175), consider PROGRAM 1.5.8 (playthattune.py). It takes notes from standard input, indexed on the chromatic scale from concert A, and plays them to standard audio. You can imagine all sorts of extensions on this basic scheme, some of which are addressed in the exercises.

We include stdaudio in our basic arsenal of programming tools because sound processing is one important application of scientific computing that is certainly familiar to you. Not only has the commercial application of digital signal processing had a phenomenal impact on modern society, but the science and engineering behind it combine physics and computer science in interesting ways. We will study more components of digital signal processing in some detail later in the book. (For example, you will learn in SECTION 2.1 how to create sounds that are more musical than the pure sounds produced by playthattune.py.)

Conceptually, one of the most significant features of the standard input, standard output, standard drawing, and standard audio data streams is that they are infinite: from the point of view of your program, there is no limit on their length. This point of view leads to more than just programs that have a long useful life (because they are less sensitive to changes in technology than programs with built-in limits). It also is related to the Turing machine, an abstract device used by theoretical computer scientists to help us understand fundamental limitations on the capabilities of real computers. One of the essential properties of the model is the idea of a finite discrete device that works with an unlimited amount of input and output.

A. If you followed the step-by-step instructions on the booksite for installing Python, these module should already be available to Python. Note that copying the files stddraw.py and stdaudio.py from the booksite and putting them in the same directory as the programs that use them is insufficient because they rely upon a library (set of modules) named Pygame for graphics and audio.

Q. Are there standard Python modules for handling standard output?

A. Actually, such features are built into Python. In Python 2, you can use the print statement to write data to stdout. In Python 3, there is no print statement; instead, there is a print() function, which is similar.

Q. Why, then, are we using the booksite stdio module for writing to standard output instead of using the features already provided by Python?

A. Our intention is to compose code that works (as much as possible) with all versions of Python. For example, using the print statement in all our programs would mean they would work with Python 2, but not with Python 3. Since we use stdio functions, we just need to make sure that we have the proper library.

Q. How about standard input?

A. There are (different) capabilities in Python 2 and Python 3 that correspond to stdio.readLine(), but nothing corresponding to stdio.readInt() and similar functions. Again, by using stdio, we can compose programs that not just take advantage of these additional capabilities, but also work in both versions of Python.

Q. How about drawing and sound?

A. Python does not come with an audio library. Python comes with a graphics library named Tkinter for producing drawings, but it is too slow for some of the graphics applications in the book. Our stddraw and stdaudio modules provide easy-to-use APIs, based on the Pygame library.

Q. So, let me get this straight; if I use the format %2.4f with stdio.writef() to write a float, I get two digits before the decimal point and four digits after the decimal point, right?

A. No, that specifies just four digits after the decimal point. The number preceding the decimal point is the width of the whole field. You want to use the format %7.2f to specify seven characters in total—four before the decimal point, the decimal point itself, and two digits after the decimal point.

Q. Which other conversion codes are there for stdio.writef()?

A. For integers, there is o for octal and x for hexadecimal. There are also numerous formats for dates and times. See the booksite for more information.

Q. Can my program reread data from standard input?

A. No. You get only one shot at it, in the same way that you cannot undo a call of stdio.writeln().

Q. What happens if my program attempts to read data from standard input after it is exhausted?

A. Python will raise an EOFError at run time. The functions stdio.isEmpty() and stdio.hasNextLine() allow you to avoid such an error by checking whether more input is available.

Q. Why does stddraw.square(x, y, r) draw a square of width 2r instead of r?

A. This makes it consistent with the function stddraw.circle(x, y, r), where the third argument is the radius of the circle, not the diameter. In this context, r is the radius of the biggest circle that can fit inside the square.

Q. What happens if my program calls stddraw.show(0)?

A. That function call tells stddraw to copy the background canvas to the standard drawing window, and then wait 0 milliseconds (that is, do not wait at all) before proceeding. That function call is appropriate if, for example, you want to run an animation at the fastest rate supported by your computer.

Q. Can I draw curves other than circles with stddraw?

A. We had to draw the line somewhere (pun intended), so we support only the basic shapes discussed in the text. You can draw other shapes one point at a time, as explored in several exercises in the text, but filling them is not directly supported.

Q. So I use negative integers to go below concert A when making input files for playthattune.py?

A. Right. Actually, our choice to put concert A at 0 is arbitrary. A popular standard, known as the MIDI Tuning Standard, starts numbering at the C five octaves below concert A. By that convention, concert A is 69 and you do not need to use negative numbers.

Q. Why do I hear weird results from standard audio when I try to sonify a sine wave with a frequency of 30,000 hertz (or more)?

A. The Nyquist frequency, defined as one-half the sampling frequency, represents the highest frequency that can be reproduced. For standard audio, the sampling frequency is 44,100, so the Nyquist frequency is 22,050.

1.5.2 Modify your program from the previous exercise to insist that the integers must be positive (by prompting the user to enter positive integers whenever the value entered is not positive).

1.5.3 Compose a program that accepts an integer n from the command line, reads n floats from standard input, and writes their mean (average value) and standard deviation (square root of the sum of the squares of their differences from the average, divided by n) to standard output.

1.5.4 Extend your program from the previous exercise to create a filter that writes all the floats in standard input that are more than 1.5 standard deviations from the mean.

1.5.5 Compose a program that reads in a sequence of integers and writes both the integer that appears in a longest consecutive run and the length of the run. For example, if the input is 1 2 2 1 5 1 1 7 7 7 7 1 1, then your program should write Longest run: 4 consecutive 7s.

1.5.6 Compose a filter that reads in a sequence of integers and writes the integers, removing repeated values that appear consecutively. For example, if the input is 1 2 2 1 5 1 1 7 7 7 7 1 1 1 1 1 1 1 1 1, your program should write 1 2 1 5 1 7 1.

1.5.7 Compose a program that takes a command-line argument n, reads from standard input N-1 distinct integers between 1 and n, and determines the missing integer.

1.5.8 Compose a program that reads in positive real numbers from standard input and writes their geometric and harmonic means. The geometric mean of n positive numbers x₁, x₂, ..., x_n is (x₁ × x₂ × ... × x_n)^1/n. The harmonic mean is (1/x₁ + 1/x₂ + ... + 1/x_n) / (1/n). Hint: For the geometric mean, consider taking logarithms to avoid overflow.

1.5.9 Suppose that the file in.txt contains the two strings F and F, and consider the following program (dragon.py):

Click here to view code image

import stdio
dragon = stdio.readString()
nogard = stdio.readString()
stdio.write(dragon + 'L' + nogard)
stdio.write(' ')
stdio.write(dragon + 'R' + nogard)
stdio.writeln()

What does the following command do (see EXERCISE 1.2.35)?

Click here to view code image

python dragon.py < in.txt | python dragon.py | python dragon.py

1.5.10 Compose a filter tenperline.py that reads a sequence of integers between 0 and 99 and writes 10 integers per line, with columns aligned. Then compose a program randomintseq.py that takes two command-line arguments m and n and outputs n random integers between 0 and m-1. Test your programs with the command python randomintseq.py 200 100 | python tenperline.py.

1.5.11 Compose a program that reads in text from standard input and writes the number of words in the text. For the purpose of this exercise, a word is a sequence of non-whitespace characters that is surrounded by whitespace.

1.5.12 Compose a program that reads in lines from standard input with each line containing a name and two integers and then uses writef() to write a table with a column of the names, the integers, and the result of dividing the first by the second, accurate to three decimal places. You could use a program like this to tabulate batting averages for baseball players or grades for students.

1.5.13 Which of the following require saving all the values from standard input (in an array, say), and which could be implemented as a filter using only a fixed number of variables? For each, the input comes from standard input and consists of n floats between 0 and 1.

• Write the maximum and minimum float.

• Write the k th smallest float.

• Write the sum of the squares of the floats.

• Write the average of the n floats.

• Write the percentage of floats greater than the average.

• Write the n floats in increasing order.

• Write the n floats in random order.

1.5.14 Compose a program that writes a table of the monthly payments, remaining principal, and interest paid for a loan, taking three numbers as command-line arguments: the number of years, the principal, and the interest rate (see EXERCISE 1.2.24).

1.5.15 Compose a program that takes three command-line arguments x, y, and z, reads from standard input a sequence of point coordinates (x_i, y_i, z_i), and writes the coordinates of the point closest to (x, y, z). Recall that the square of the distance between (x, y, z) and (x_i, y_i, z_i) is (x – x_i)² + (y – y_i)² + (z – z_i)². For efficiency, do not use either math.sqrt() or the ** operator.

1.5.16 Compose a program that, given the positions and masses of a sequence of objects, computes their center-of-mass, or centroid. The centroid is the average position of the n objects, weighted by mass. If the positions and masses are given by (x_i, y_i, m_i), then the centroid (x, y, m) is given by

m = m₁ + m₂ + ... + m_n

x = (m₁ x₁ + ... + m_n x_n) / m

y = (m₁ y₁ + ... + m_n y_n) / m

1.5.17 Compose a program that reads in a sequence of real numbers between –1 and +1 and writes their average magnitude, average power, and the number of zero crossings. The average magnitude is the average of the absolute values of the data values. The average power is the average of the squares of the data values. The number of zero crossings is the number of times a data value transitions from a strictly negative number to a strictly positive number, or vice versa. These three statistics are widely used to analyze digital signals.

1.5.18 Compose a program that takes a command-line argument n and plots an n-by-n checkerboard with red and black squares. Color the lower-left square red.

1.5.19 Compose a program that takes as command-line arguments an integer n and a float p (between 0 and 1), plots n equally spaced points of size on the circumference of a circle, and then, with probability p for each pair of points, draws a gray line connecting them.

1.5.20 Compose code to draw hearts, spades, clubs, and diamonds. To draw a heart, draw a diamond, then attach two semicircles to the upper-left and upper-right sides.

1.5.21 Compose a program that takes a command-line argument n and plots a “flower” with n petals (if n is odd) or 2n petals (if n is even), by plotting the polar coordinates (r, θ) of the function r = sin(n θ) for θ ranging from 0 to 2π radians.

1.5.22 Compose a program that takes a string s from the command line and displays it in banner style on the screen, moving from left to right and wrapping back to the beginning of the string as the end is reached. Add a second command-line argument to control the speed.

1.5.23 Modify playthattune.py to take additional command-line arguments that control the volume (multiply each sample value by the volume) and the tempo (multiply each note’s duration by the tempo).

1.5.24 Compose a program that takes the name of a .wav file and a playback rate r as command-line arguments and plays the file at the given rate. First, use stdaudio.read() to read the file into an array a[]. If r = 1, just play a[]; otherwise, create a new array b[] of approximate size r times a.length. If r < 1, populate b[] by sampling from the original; if r > 1, populate b[] by interpolating from the original. Then play b[].

1.5.25 Compose programs that use stddraw to create each of these designs.

1.5.26 Compose a program circles.py that draws filled circles of random size at random positions in the unit square, producing images like those below. Your program should take four command-line arguments: the number of circles, the probability that each circle is black, the minimum radius, and the maximum radius.

1.5.28 Statistical polling. When collecting statistical data for certain political polls, it is very important to obtain an unbiased sample of registered voters. Assume that you have a file with n registered voters, one per line. Compose a filter that writes a random sample of size m (see sample.py, PROGRAM 1.4.1).

1.5.29 Terrain analysis. Suppose that a terrain is represented by a two-dimensional grid of elevation values (in meters). A peak is a grid point whose four neighboring cells (left, right, up, and down) have strictly lower elevation values. Compose a program peaks.py that reads a terrain from standard input and then computes and writes the number of peaks in the terrain.

1.5.30 Histogram. Suppose that the standard input stream is a sequence of floats. Compose a program that takes an integer n and two floats lo and hi from the command line and uses stddraw to plot a histogram of the count of the numbers in the standard input stream that fall into each of the n intervals defined by dividing the interval (lo, hi) into n equal-sized intervals.

1.5.31 Spirographs. Compose a program that takes three parameters R, r, and a from the command line and draws the resulting spirograph. A spirograph (technically, an epicycloid) is a curve formed by rolling a circle of radius r around a larger fixed circle of radius R. If the pen offset from the center of the rolling circle is (r+a), then the equation of the resulting curve at time t is given by

x(t) = (R + r) cos (t) – (r + a) cos ((R + r)t /r)

y(t) = (R + r) sin (t) – (r + a) sin ((R + r)t /r)

Such curves were popularized by a best-selling toy that contains discs with gear teeth on the edges and small holes that you could put a pen in to trace spirographs.

1.5.32 Clock. Compose a program that displays an animation of the second, minute, and hour hands of an analog clock. Use the call stddraw.show(1000) to update the display roughly once per second.

1.5.33 Oscilloscope. Compose a program to simulate the output of an oscilloscope and produce Lissajous patterns. These patterns are named after the French physicist, Jules A. Lissajous, who studied the patterns that arise when two mutually perpendicular periodic disturbances occur simultaneously. Assume that the inputs are sinusoidal, so that the following parametric equations describe the curve:

x(t) = A_x sin (w_x t + θ_x)

y(t) = A _y sin (w_y t + θ_y)

Take the six parameters A_x and A_y (amplitudes); w_x and w_y (angular velocity); and θ_x and θ_y (phase factors) from the command line.

1.5.34 Bouncing ball with tracks. Modify bouncingball.py to produce images like the ones shown in the text, which show the track of the ball on a gray background.

1.5.35 Bouncing ball with gravity. Modify bouncingball.py to incorporate gravity in the vertical direction. Add calls to stdaudio.playFile() to add one sound effect when the ball hits a wall and a different one when it hits the floor.

1.5.36 Random tunes. Compose a program that uses stdaudio to play random tunes. Experiment with keeping in key, assigning high probabilities to whole steps, repetition, and other rules to produce reasonable melodies.

1.5.37 Tile patterns. Using your solution to EXERCISE 1.5.25, compose a program tilepattern.py that takes a command-line argument n and draws an n-by-n pattern, using the tile of your choice. Add a second command-line argument that adds a checkerboard option. Add a third command-line argument for color selection. Using the patterns on the facing page as a starting point, design a tile floor. Be creative! Note: These are all designs from antiquity that you can find in many ancient (and modern) buildings.

The model, which is known as the random surfer model, is simple to describe. We consider the web to be a fixed set of pages, with each page containing a fixed set of links, and each link a reference to some other page. We study what happens to a person (the random surfer) who randomly moves from page to page, either by typing a page name into the address bar or by clicking a link on the current page.

The underlying mathematical model behind the web model is known as the graph, which we will consider in detail at the end of the book (in SECTION 4.5). We defer discussion of the details of processing graphs until then. Instead, we concentrate on calculations associated with a natural and well-studied probabilistic model that accurately describes the behavior of the random surfer.

The first step in studying the random surfer model is to formulate it more precisely. The crux of the matter is to specify what it means to randomly move from page to page. The following intuitive 90-10 rule captures both methods of moving to a new page: Assume that 90 percent of the time the random surfer clicks a random link on the current page (each link chosen with equal probability) and that 10 percent of the time the random surfer goes directly to a random page (all pages on the web chosen with equal probability).

You can immediately see that this model has flaws, because you know from your own experience that the behavior of a real web surfer is not quite so simple:

• No one chooses links or pages with equal probability.

• There is no real potential to surf directly to each page on the web.

• The 90-10 (or any fixed) breakdown is just a guess.

• It does not take the “back” button or bookmarks into account.

• We can afford to work with only a small sample of the web.

Despite these flaws, the model is sufficiently rich that computer scientists have learned a great deal about properties of the web by studying it. To appreciate the model, consider the small example on the previous page. Which page do you think the random surfer is most likely to visit?

Each person using the web behaves a bit like the random surfer, so understanding the fate of the random surfer is of intense interest to people building web infrastructure and web applications. The model is a tool for understanding the experience of each of the billions of web users. In this section, you will use the basic programming tools from this chapter to study the model and its implications.

Later in the book, you will learn how to read web pages in Python programs (SECTION 3.1) and to convert from names to numbers (SECTION 4.4) as well as other techniques for efficient graph processing. For now, we assume that there are n web pages, numbered from 0 to n – 1, and we represent links with ordered pairs of such numbers, the first specifying the page containing the link and the second specifying the page to which it refers. Given these conventions, a straightforward input format for the random surfer problem is an input stream consisting of an integer (the value of n) followed by a sequence of pairs of integers (the representations of all the links). Because of the way stdio reading functions treat whitespace characters, we are free to either put one link per line or arrange them several to a line.

• Read n, and then create arrays linkCounts[][] and outDegrees[].

• Read the links and accumulate counts so that linkCounts[i][j] counts links from i to j and outDegrees[i] counts links from i to anywhere.

• Use the 90-10 rule to compute the probabilities.

The first two steps are elementary, and the third is not much more difficult: multiply linkCounts[i][j] by 0.90/outDegrees[i] if there is a link from i to j (take a random link with probability 0.9), and then add 0.10/n to each element (go to a random page with probability 0.1). PROGRAM 1.6.1 (transition.py) performs this calculation: it is a filter that converts the list-of-links representation of a web model into a transition-matrix representation.

The transition matrix is significant because each row represents a discrete probability distribution—the elements fully specify the behavior of the random surfer’s next move, giving the probability of surfing to each page. Note in particular that the elements sum to 1 (the surfer always goes somewhere).

The output of transition.py defines another file format, one for matrices of float values: the numbers of rows and columns followed by the values for matrix elements. Now, we can compose code that reads and processes transition matrices.

total = 0.0
for j in range(0, n)
    total += p[page][j]
    if r < total:
        page = j
        break

The variable total tracks the endpoints of the intervals defined in row p[page], and the for loop finds the interval containing the random value r. For example, suppose that the surfer is at page 4 in our example. The transition probabilities are 0.47, 0.02, 0.47, 0.02, and 0.02, and total takes on the values 0.0, 0.47, 0.49, 0.96, 0.98, and 1.0. These values indicate that the probabilities define the five intervals (0, 0.47), (0.47, 0.49), (0.49, 0.96), (0.96, 0.98), and (0.98, 1), one for each page. Now, suppose that random.random() returns the value 0.71. We increment j from 0 to 1 to 2 and stop there, which indicates that 0.71 is in the interval (0.49, 0.96), so we send the surfer to the third page (page 2). Then, we perform the same computation for p[2], and the random surfer is off and surfing. For large n, we can use binary search to substantially speed up this computation (see EXERCISE 4.2.35). Typically, we are interested in speeding up the search in this situation because we are likely to need a huge number of random moves, as you will see.

Not all Markov chains have this mixing property. For example, if we eliminate the random leap from our model, certain configurations of web pages can present problems for the surfer. Indeed, there exist on the web sets of pages known as spider traps, which are designed to attract incoming links but have no outgoing links. Without the random leap, the surfer could get stuck in a spider trap. The primary purpose of the 90-10 rule is to guarantee mixing and eliminate such anomalies.

Accurate page rank estimates for the web are valuable in practice for many reasons. First, using them to put in order the pages that match the search criteria for web searches proved to be vastly more in line with people’s expectations than previous methods. Next, this measure of confidence and reliability led to the investment of huge amounts of money in web advertising based on page ranks. Even in our tiny example, page ranks might be used to convince advertisers to pay up to four times as much to place an ad on page 0 as on page 4. Computing page ranks is mathematically sound, an interesting computer science problem, and big business, all rolled into one.

Click here to view code image

if i % 1000 == 0:
    stddraw.clear()
    for k in range(n):
        stddraw.filledRectangle(k - 0.25, 0.0, 0.5, hits[k])
stddraw.show(10)

to the random move loop in randomsurfer.py and run the code for, say, millions of moves, you will see a drawing of the frequency histogram that eventually stabilizes to the page ranks. (The constants 1000 and 10 in this code are a bit arbitrary; you might wish to change them when you run the code.) After you have used this tool once, you are likely to find yourself using it every time you want to study a new model (perhaps with some minor adjustments to handle larger models).

Directly simulating the behavior of a random surfer to understand the structure of the web is appealing, but it has limitations. Think about the following question: Could you use it to compute page ranks for a web model with millions (or billions!) of web pages and links? The quick answer to this question is no, because you cannot even afford to store the transition matrix for such a large number of pages. A matrix for millions of pages would have trillions of elements. Do you have that much space on your computer? Could you use randomsurfer.py to find page ranks for a smaller model with, say, thousands of pages? To answer this question, you might run multiple simulations, record the results for a large number of trials, and then interpret those experimental results. We do use this approach for many scientific problems (the gambler’s ruin problem is one example; SECTION 2.4 is devoted to another), but it can be very time-consuming, as a huge number of trials may be necessary to get the desired accuracy. Even for our tiny example, we saw that it takes millions of iterations to get the page ranks accurate to three or four decimal places. For larger models, the required number of iterations to obtain accurate estimates becomes truly huge.

[1.0 0.0 0.0 0.0 0.0]

which specifies that the random surfer starts on page 0. Multiplying this vector by the transition matrix gives the vector

[.02 .92 .02 .02 .02]

which is the probabilities that the surfer winds up on each of the pages after one step. Now, multiplying this vector by the transition matrix gives the vector

[.05 .04 .36 .37 .19]

which contains the probabilities that the surfer winds up on each of the pages after two steps. For example, the probability of moving from 0 to 2 in two moves is the probability of moving from 0 to 0 to 2 (0.02 × 0.02), plus the probability of moving from 0 to 1 to 2 (0.92 × 0.38), plus the probability of moving from 0 to 2 to 2 (0.02 × 0.02), plus the probability of moving from 0 to 3 to 2 (0.02 × 0.02), plus the probability of moving from 0 to 4 to 2 (0.02 × 0.47), which adds up to a grand total of 0.36. From these initial calculations, the pattern is clear: the vector giving the probabilities that the random surfer is at each page after t steps is precisely the product of the corresponding vector for t – 1 steps and the transition matrix. By the basic limit theorem for Markov chains, this process converges to the same vector no matter where we start; in other words, after a sufficient number of moves, the probability that the surfer ends up on any given page is independent of the starting point.

PROGRAM 1.6.3 (markov.py) is an implementation that you can use to check convergence for our example. For instance, it gets the same results (the page ranks accurate to two decimal places) as randomsurfer.py, but with just 20 matrix–vector multiplications instead of the tens of thousands of iterations needed by randomsurfer.py. Another 20 multiplications gives the results accurate to three decimal places, as compared with millions of iterations for randomsurfer.py, and just a few more multiplications give the results to full precision (see EXERCISE 1.6.6).

Markov chains are well studied, but their impact on the web was not truly felt until 1998, when two graduate students, Sergey Brin and Lawrence Page, had the audacity to build a Markov chain and compute the probabilities that a random surfer hits each page for the whole web. Their work revolutionized web search and is the basis for the page ranking method used by Google, the highly successful web search company that they founded. Specifically, the company periodically recomputes the random surfer’s probability for each page. Then, when you do a search, it lists the pages related to your search keywords in order of these ranks. Such page ranks now predominate because they somehow correspond to the expectations of typical web users, reliably providing them with relevant web pages for typical searches. The computation that is involved is enormously time-consuming, due to the huge number of pages on the web, but the result has turned out to be enormously profitable and well worth the expense. The method used in markov.py is far more efficient than simulating the behavior of a random surfer, but it is still too slow to actually compute the probabilities for a huge matrix corresponding to all the pages on the web. That computation is enabled by better data structures for graphs (see CHAPTER 4).

Perhaps the most important lesson to learn from composing programs for complicated problems like the example in this section is that debugging is difficult. The polished programs in the book mask that lesson, but you can rest assured that each one is the product of a long bout of testing, fixing bugs, and running the program on numerous inputs. Generally we avoid describing bugs and the process of fixing them in the text because that makes for a boring account and overly focuses attention on bad code, but you can find some examples and descriptions in the exercises and on the booksite.

1.6.2 Modify transition.py to ignore the effect of multiple links. That is, if there are multiple links from one page to another, count them as one link. Create an example that shows how this modification can change the order of page ranks.

1.6.3 Modify transition.py to handle pages with no outgoing links, by filling rows corresponding to such pages with the value 1/n.

1.6.4 The code fragment in randomsurfer.py that generates the random move fails if the probabilities in the row p[page] do not add up to 1. Explain what happens in that case, and suggest a way to fix the problem.

1.6.5 Determine, to within a factor of 10, the number of iterations required by randomsurfer.py to compute page ranks to four decimal places and to five decimal places for tiny.txt.

1.6.6 Determine the number of iterations required by markov.py to compute page ranks to three decimal places, to four decimal places, and to ten decimal places for tiny.txt.

1.6.7 Download the file medium.txt from the booksite (which reflects the 50-page example depicted in this section) and add to it links from page 23 to every other page. Observe the effect on the page ranks, and discuss the result.

1.6.8 Add to medium.txt (see the previous exercise) links to page 23 from every other page, observe the effect on the page ranks, and discuss the result.

1.6.9 Suppose that your page is page 23 in medium.txt. Is there a link that you could add from your page to some other page that would raise the rank of your page?

1.6.10 Suppose that your page is page 23 in medium.txt. Is there a link that you could add from your page to some other page that would lower the rank of that page?

1.6.11 Use transition.py and randomsurfer.py to determine the transition probabilities for the eight-page example shown below.

1.6.12 Use transition.py and markov.py to determine the transition probabilities for the eight-page example shown below.

1.6.14 Random web. Compose a generator for transition.py that takes as input a page count n and a link count m and writes to standard output n followed by m random pairs of integers between 0 and n-1. (See SECTION 4.5 for a discussion of more realistic web models.)

1.6.15 Hubs and authorities. Add to your generator from the previous exercise a fixed number of hubs, which have links pointing to them from 10 percent of the pages, chosen at random, and authorities, which have links pointing from them to 10 percent of the pages. Compute page ranks. Which rank higher, hubs or authorities?

1.6.16 Page ranks. Design an array of pages and links where the highest-ranking page has fewer links pointing to it than some other page.

1.6.17 Hitting time. The hitting time for a page is the expected number of moves between times the random surfer visits the page. Run experiments to estimate page hitting times for tiny.txt, compare the results with page ranks, formulate a hypothesis about the relationship, and test your hypothesis on medium.txt.

1.6.18 Cover time. Compose a program that estimates the time required for the random surfer to visit every page at least once, starting from a random page.

1.6.19 Graphical simulation. Create a graphical simulation where the size of the dot representing each page is proportional to its rank. To make your program data-driven, design a file format that includes coordinates specifying where each page should be drawn. Test your program on medium.txt.