Solving an Escape Riddle Using a Monte Carlo Simulation in Python

I recently heard this riddle from a friend at a campfire. I wasn't smart enough to come up with the answer, but when it was explained to me I wanted to know how long it would take the prisoners to be released. In this post I'll describe the riddle and the answer, and show how to use Python to set up an experiment, use a Monte Carlo simulator to run the experiment repeatedly, and briefly discuss some of the results.

The Riddle

There are 50 prisoners, and one prison keeper. The prisoners will be placed in separate cells, never able to communicate. Once per day, the prison keeper will randomly select one prisoner and take him into a room with two levers. The levers have two positions each: Up, and down. The prisoner may operate the levers however he wishes, then he must return to his cell. (The levers don't do anything; they just represent a binary state.) The prisoners have one chance to indicate to the prison keeper whether all the prisoners have been in the lever-room. If they so indicate too early, they will never be released. If they're correct (even if they're late), they will be released immediately. The prisoners get one chance to meet prior to imprisonment, to discuss strategy. How can they ensure their release?

The Answer

Only one lever is needed. Elect one prisoner to be a 'counter.' When a prisoner enters the lever room and the lever is in the 'down' position and he has never moved the lever, then and only then will he move the lever to the 'up' position. Only the 'counter' prisoner may return the lever to the 'down' position, and each time he does so he increments a counter. When the counter reaches 50, he knows that all prisoners have been into the lever room at least once. He can then inform the prison guard, and all the prisoners will be released.

A Solution in Python

I came up with the following model in Python to describe the prison scenario. There are probably better and more optimized way to represent this situation, but here's the first one I came up with.

#Iteration 1: The basic scenario
import random
day_number = 0 
prisoner_designee = 1 # chosen by prisoners to be the official counter
counter = 0 # as tracked by the prisoner_designee
lever_position = 'down' # starting position of the lever
prisoners_who_have_operated_lever = []

while counter < 50: 
  prisoner_number = random.randint(1,50)    # Select a prisoner
  if not prisoner_number in prisoners_who_have_operated_lever:
    if lever_position == 'down':
      prisoners_who_have_operated_lever.append(prisoner_number)
      lever_position = 'up'
  if prisoner_number == prisoner_designee:  # If it's the designee...
    if lever_position == 'up':              # ...and the lever is up...
        lever_position = 'down'             # ...put the lever back down...
        counter += 1                        # ...and increment the counter.
  day_number += 1

print day_number

The above code, iteration 1, uses random.randint to generate a random integer. I used if statements to describe the logic, and printed out the result. Example:

$ python prisoner1.py 
2536

This means that for this particular run, it would have taken the prisoners (or rather, the prisoners' designated counter prisoner_designee) 2,536 days to count 50 'up' levers, proving that all prisoners had entered the lever room and ensuring their release. That's just under seven years(!).

But this result depends on the order in which the prison keeper (in our case represented by random.randint) chose the prisoners. Running this code at the command line produces different results each time, usually in the range of 2500--3000 (days). To get a better picture of how long it would take a new group of prisoners to earn their release we'd have to simulate the whole range of possible outputs (or which there are nearly infinitely many) But I don't want to run this manually dozens of times and tabulate the results by hand; I want to make the code run the simulation thousands of times, with random inputs each time, and then aggregate the results. That's the Monte Carlo method.

Implementing the Monte Carlo Method

We're going to use a loop to run the experiment, so first, wrap the experiment in a function:

#Iteration 2: Wrapping the simulation in a function
import random

def run_simulation():
  day_number = 0 
  prisoner_designee = 1 # chosen by prisoners to be the official counter
  counter = 0 # as tracked by the prisoner_designee
  lever_position = 'down' # starting position of the lever
  prisoners_who_have_operated_lever = []

  while counter < 50: 
    prisoner_number = random.randint(1,50)    # Select a prisoner
    if not prisoner_number in prisoners_who_have_operated_lever:
      if lever_position == 'down':
        prisoners_who_have_operated_lever.append(prisoner_number)
        lever_position = 'up'
    if prisoner_number == prisoner_designee:  # If it's the designee...
      if lever_position == 'up':              # ...and the lever is up...
          lever_position = 'down'             # ...put the lever back down...
          counter += 1                        # ...and increment the counter.
    day_number += 1
  return day_number

print run_simulation()

Running this code should produce the same type of output (with a different result, of course):

$ python prisoner2.py 
2950

The next step is to write a loop that will call this function a thousand times:

#Iteration 3: Loop it to run 1000 times
import random

def run_simulation():
  counter = 0 
  day_number = 0 
  prisoner_chief = 1 
  lever_position = 'down'
  prisoners_who_have_operated_lever = []

  while counter < 50: 
    prisoner_number = random.randint(1,50) # Select a prisoner
    if not prisoner_number in prisoners_who_have_operated_lever:
      if lever_position == 'down':
        prisoners_who_have_operated_lever.append(prisoner_number)
        lever_position = 'up'
    if prisoner_number == prisoner_chief:
      if lever_position == 'up':
          lever_position = 'down'
          counter += 1
    day_number += 1
  return day_number 

simulation_results = []
simulations = 0 
while simulations < 1000:
  simulation_results.append(run_simulation())
  simulations += 1
print simulation_results

In this iteration we use a while loop to run the simulation 1000 times, storing the results in the list called simulation_results. Here's an example of the output:

$ python prisoner3.py 
[2444, 2119, 2818, 2253, 2586, 2543, 2490, 3388, 2034, 2739, 2554, 2585, 2498, 2689, 3180, 2760, 3145, 2698, 2196, 2769, 2400, 2783, 3091, 2258, 2952, 1730, 2974, 2656, 3059, 2000, 3222, 2186, 3114, 2618, .... ]

(List truncated.)
The last step in the Monte Carlo method is to aggregate and interpret the results. In this case, that means taking the average of all the simulation runs. A simple way to do this with lists is to sum the list :

average = sum(list)/len(list)

A more elegant way is to use Numpy's built-in mean function, which requires first converting the list to a Numpy array:

import numpy

array = numpy.array(list)

average = numpy.mean(array)

The next iteration of the code will include this small refactorization. Running this code produces the following result:

$ python prisoner3.py 
2720.895

Now, instead of one experiment, this result represents the average of 1000 experiments.

So what?

So on average, it would take 50 prisoners in this scenario 2721 days to earn their release. Cool. Neat. But this raises so many more questions: Is this number related to the probability of )1/50)*(1/50)? Where is this average value in the range of values yielded from the 1000 results? If we halve the number of prisoners, will the number of days go down by half? What does a histogram of the results look like? The answers, as well as a more in-depth discussion of other statistics, will be in the next post.

Installing Twitter Bootstrap in Flask 0.9

I have a small web application called 'dispatchtool' written in the Python microframework Flask. I want to use Twitter's Bootstrap CSS framework. In this post I will explain how to do this.

I am running Ubuntu 12.04, using Python 2.7, and Flask 0.9. This tutorial assumes you've got a basic Flask application up and running, and that you can view the development version on your local machine. It also assumes that you're using templates. (This is important!)

Step 1: Get a basic CSS file working.

Before dealing with Bootstrap, we want to ensure sure that the application can correctly find and use a simple test CSS file.
To do this, first, create a new directory in your Flask application's root directory called 'static'. My app is called 'dispatchtool', so I create 'static' inside the 'dispatchtool' directory:

$ cd dispatchtool

$ mkdir static 

Next, in that new directory, create the test css file. I'll call mine 'test.css'. Add some simple style that will be loud and obvious. It's only purpose will be to ensure that we can get the application to see this file and use it. Here's my 'test.css':

h1 {color: red}

p {color: green}

Of course, make sure that your app includes an <h1> tag and a <p> tag somewhere, or else these styles won't be applied, even if the stylesheet is properly set up.

Now that we've got our test css file in the right place, we need to tell Flask how to find it. CSS files are always linked to HTML pages in the <head> section of each page using the <link … /> tag. In my application, I'll add the following tag:

<link type="text/css" rel="stylesheet" href='/static/test.css' />

Here's the top of my page, for context:

<!DOCTYPE html>

<html lang="en">
  <head>
    <title>The Demand Response Dispatch Tool</title>
    <link type="text/css" rel="stylesheet" href='static/test.css' />
  </head>
  <body

Now open up your page in a browser, and see if the styles were applied. In my case, it suddenly looks like Christmas. If nothing changed, look at the source (ctrl-U) and click on the link you created to the stylesheet. That will show where Flask thinks the CSS file is.

Even though that's working, we're not quite done yet. Web application frameworks like Flask have a better way of creating links to static things like CSS files. In Flask, it's a method called url_for(). Methods like this make your application more flexible, allowing you to move around pages without breaking links to your CSS files. So now, replace the value of the 'href' attribute in the <link> tag with the following line:

{{ url_for('static',filename='test.css') }}

The double braces tell Flask a) that the contents are to parsed as Python code and b) to render the result directly into the HTML output. Now the whole link tag should look like this:

<link type="text/css" rel="stylesheet" href="{{ url_for('static', filename='test.css') }}" />

Now save the file and refresh your browser page. If that still renders your page using your test CSS file, then we're most of the way done. If it doesn't, check the page source, and ensure that the <link> tag is pointing to the right place.

Step 2: Download and "install" Twitter Bootstrap

The next step is simple and short. We need to get the Bootstrap files (download) and put them in the right place (install). Normal people will navigate to the Bootstrap page, click the link, and unzip the files, then put those files (as-is, with no changes to the directory structure) in the 'static' folder that we created earlier. Super elite hackers will do it the following way:

$ cd ../static/

$ wget http://twitter.github.com/bootstrap/assets/bootstrap.zip

$ unzip bootstrap.zip

$ rm -r bootstrap.zip 

Whichever way you do it, you should now have a folder called 'bootstrap' inside your 'static' folder. That's it. You have downloaded and 'installed' Bootstrap.

Step 3: Link to the boostrap.css stylesheet

The main Bootstrap stylesheet is in bootstrap/css/, and it's called 'bootstrap.css'. We want to change our link from 'test.css' to this file, so update your <link> tag as follows:

<link type="text/css" rel="stylesheet" href="{{ url_for('static', filename='bootstrap/css/bootstrap.css') }}" />

That's it. If you reload your page, you should see at least some changes reflecting Twitter's style. (When you go to production, you should switch from the general 'bootstrap.css' to the minified 'bootstrap.min.css' to save bandwidth.)

Step 4: Update your HTML to hook into the Bootstrap styles

Now that we've got our Flask app properly linking to the Boostrap css file, all that remains is to use it. That means we have to look at the boostrap.css file, and see what classes and ids we need to assign to our HTML elements. I'll go through two quick examples.

Container

The most important is probably to use Bootstrap's "container" class. Simply add the class 'container' the <body> tag:

<body class="container">

When you refresh the page, you should see your content contained in a 940 pixel-wide column.

The Hero Unit

The most famous Bootstrap style is probably the "hero unit". To see it, just add the class 'hero-unit' to any block-level element:

<div class="hero-unit"><p>Don't abuse this style!</p></div>

Conclusion

We walked through how to set up general CSS file in a Flask app, then how to acquire and link Bootstrap to the same Flask app. The next step is to check out the Bootstrap Documentation, and start choosing the styles you want to use. Thanks for reading!

Basic Unit Testing Tutorial in Python 3.2

I have a Python 3.2 program that's in a working state, but I never wrote any tests for it. Now, I want to go back and provide some code coverage so that as I go back to refactor and add functionality, I'll be more certain I'm not breaking anything. This is a fairly common case--you have existing code, and you want to write some tests for it.

Environment

I'm running Python 3.2 in a Windows 7 environment, using the IDLE Python code editor. I typically run code directly from IDLE by using the F5 key. Because of how IDLE invokes the Python interpreter, some of the code examples below will run differently depending on whether code is invoked from IDLE or from the command line. To enable Python from the Windows command line, it must be added to the PATH. To do this, open a command line window (Start -> Accessories -> Command Line) and type path <path to your python installation>. For example, since my python 3.2 installation can be found at C:\Python32, I issue the following path command:

>path C:\Python32

This adds the Python interpreter to the path so it can be invoked from within your Python project directory. Note that it only changes the path for this one command line session, so if you close and re-open the command line window you'll have to re-issue the path command.

Python's UnitTest Module

The first issue is what testing framework to use. I'm going to use Python's built-in module called UnitTest. It's a one-stop shop for all the testing most pure-Python programs will need. More advanced testing frameworks are available. The two most popular seem to be py.test and nose. I want to explore testing in a more 'raw' state, so I'll save these frameworks until I'm comfortable using unittest on its own.

Writing Testable Code

It's easiest to write a test for a small, discreet bit of code that accomplishes a well-defined task. Python best practices dictate that code be broken into small functions that perform one task at one level of abstraction. If you follow this advice and write your code using many small functions that each perform one well-defined task, then those functions will be easy to test. If instead you have written large functions that combine many steps, then first consider breaking these functions down into smaller ones.

Pick a small function to test

Start small by picking a short, well-understood function that does something simple. I'll use a function that takes one argument, uses an if statement to choose a "rate", and returns that rate:

def choose_renovation_rate(self, years_since_last_renovation):
    if years_since_last_renovation < 7:
        rate = 0
    elif years_since_last_renovation < 15:
        rate = 0.01
    elif years_since_last_renovation < 25:
        rate = 0.05
    elif years_since_last_renovation < 50:
        rate = 0.07
    else:
        rate = 0.1
    return rate

Create a new file to hold the test code

To get started with the testing itself, create a new file, and for now make sure it's in the same directory as the code you wish to test. I called mine unit_tests.py. Import the testunit module, and start a new class that describes what you will be testing. This new class will be a sub-class of unittest.TestCase. This class will hold several small test cases.

import unittest

class TestRateFunctions(unittest.TestCase):
    pass
    # test cases go here

Write a test function title and description

Next, inside the TestRateFunctions class I just defined, I will define a function and describe in words what it will test. The function names, by convention and to allow the testing machinery to run properly, should start with test_. For the moment, I will use a single assertTrue statement and pass it the condition True, just to see if everything is working.

import unittest

class TestRateFunctions(unittest.TestCase):
    def test_renovation_chooser_should_return_correct_rate(self):
        # This function should return the correct renovation rate
        self.assertTrue(True)

This is now a trivial, but functioning test.

Run the test

To actually run the test, we need to add one line to the end of the file: unittest.main(). With this line in place, we can simply call this test one of two ways.

Via the command line

After issuing the PATH commands as discussed under "Environment" above, in the command line terminal navigate to the directory where you created the file and call it:

C:\lighting_floor_space_stock_model>python unit_tests.py
.
----------------------------------------------------------------------
Ran 1 test in 0.000s

OK

The output consists of

1. A period. This represents one test, which in this case was test_renovation_chooser. With more tests written, these dots acts as a progress indicator for the test suite.
2. A long line of hyphens. This is just a visual separator representing the end of the execution of tests.
3. A report, which in this case Ran 1 test in 0.000s.
4. The phrase OK. We only get this if all tests pass.

Via IDLE

In the IDLE code editor, you can just press F5 to run the code directly. This produces the following output:

.
----------------------------------------------------------------------
Ran 1 test in 0.004s

OK
Traceback (most recent call last):
  File "C:\lighting_floor_space_stock_model\unit_tests.py", line 8, in <module>
    unittest.main()
  File "C:\Python32\lib\unittest\main.py", line 124, in __init__
    self.runTests()
  File "C:\Python32\lib\unittest\main.py", line 272, in runTests
    sys.exit(not self.result.wasSuccessful())
SystemExit: False

As you can see, when we ran the test suite from inside IDLE, we got the same successful output telling us the test passed (OK!), but it's followed by an ugly error and a traceback. The error is because unittest wants to "exit" but it can't, because the default IDLE behavior is to keep running. There are two simple ways around this error:

1. One solution is to just catch the error and pass around it. Simply change the last line to following:
   

   try: unittest.main()
   except SystemExit: pass

2. Alternatively, we could keep the first line as-is (without adding the try/except statements or the second line) and just add an argument to the original line, as follows:
   

   unittest.main(exit = False)

With either of those two modifications in place, the test suite should run properly from within IDLE.

Now Test Your Function

Now that we're sure we've got all the machinery operating correctly (import unittest, create a subclass, define a test case starting with test_, make sure it runs properly), we have to write code that actually tests our function.

In order to call the function in question (choose_renovation_rate) inside the testing function (test_renovation_chooser_should_return_correct_rate), we need to first make sure its parent class is available. In this case, choose_renovation_rate is a function inside the class FloorSpace(), so first I need to import the FloorSpare class with from floor_space import *. Then I can create a new variable called self.rate that holds the result of the FloorSpace.choose_renovation_rate function. The whole file now appears as follows, using assertTrue:

import unittest
from floor_space import FloorSpace

class TestRateFunctions(unittest.TestCase):
    def test_renovation_chooser_should_return_correct_rate(self):
        # A 2-year old building has a 0% chance of renovation:
        self.rate = FloorSpace.choose_renovation_rate(self, 2)
        self.assertTrue(self.rate == 0)

unittest.main(exit = False)

A Quick Refactor

Let's make one short change. Since we're testing for equality (using the == operator), we could use the assertEqual function instead of assertTrue. Furthermore, we could a string argument to the end of the function that will get printed if the test fails, allowing us to pass more detailed and useful information to the user. Change the last line of the test function as follows:

self.assertEqual(self.rate, 0, "2-year-old building should have a 0% renovation rate.")

This is such a simple test method that the message is nearly redunant to just reading the code, but it serves to illustrate the method of adding an information string to the assert method.

Refining the Test

Consider again the function we're testing:

def choose_renovation_rate(self, years_since_last_renovation):
    if years_since_last_renovation < 7:
        rate = 0
    elif years_since_last_renovation < 15:
        rate = 0.01
    elif years_since_last_renovation < 25:
        rate = 0.05
    elif years_since_last_renovation < 50:
        rate = 0.07
    else:
        rate = 0.1
    return rate

We pass it the number of years since it was last renovated (in years), and it returns a renovation rate. That renovation rate is used by other functions to perform various other computations. We could test each of the if conditions, but those numbers might change if we later tune the function to match real-world data. It might be more useful to test, for example, that the function always returns a percentage (i.e. a float between zero and one). So let's test a few of the if conditions just for completeness, and add a new test function to check whether it always returns a percentage between 0 and 1. Note that instead of testing each year we put the assert methods inside a for loop.

import unittest
from floor_space import FloorSpace

class TestRateFunctions(unittest.TestCase):
    def test_renovation_chooser_should_return_correct_rate(self):
        # A 2-year old building has a 0% chance of renovation:
        self.rate = FloorSpace.choose_renovation_rate(self, 2)
        self.assertEqual(self.rate, 0, "2-year-old building should have a 0% renovation rate.")

        # A 13-year old building has a 1% chance of renovation:
        self.rate = FloorSpace.choose_renovation_rate(self, 13)
        self.assertEqual(self.rate, 0.01, "13-year-old building should have a 1% renovation rate.")

        # A 55-year old building has a 10% chance of renovation:
        self.rate = FloorSpace.choose_renovation_rate(self, 55)
        self.assertEqual(self.rate, 0.1, "55-year-old building should have a 10% renovation rate.")

    def test_renovation_chooser_should_return_percentage(self):
        # The function should always return a rate such that 0 <= rate <= 1
        for years_since_renovation in range(100):
            self.rate = FloorSpace.choose_renovation_rate(self, years_since_renovation)
            self.assertGreaterEqual(self.rate,0, "renovation rate should be >= 0")
            self.assertLessEqual(self.rate,1, "renovation rate should be <= 1")

unittest.main(exit = False)

For a list of assert methods made available to unittest, see the documentation.

Conclusion

The basic steps to writing unit tests in Python are:

1. Write testable code by keeping functions short; each function should perform one task
2. Create a new file for the test code (e.g. unit_tests.py)
3. Import unittest
4. Subclass unittest.TestCase (e.g. class TestRateFunctions(unittest.TestCase))
5. Create new test functions with description names starting with test_ (e.g.
   def test_renovation_chooser_should_return_percentage)
   
6. Make your "work" functions available to the "test" functions by importing where necessary (e.g. from floor_space import FloorSpace)
7. Add string arguments to pass relevant information to the user in the case of a failed test.
8. Refine and refactor for readability and usefulness.

Finally, don't forget to actually run the tests frequently during development!

Basic linear algebra with Numpy on Python 3.2

I'm taking Stanford's online course on machine learning through Coursera. The second week has a good overview of linear algebra and matrix operations. The instructor has provided a useful PowerPoint deck in which he explains the basics. I'm going to go through this pdf and implement the linear algebra using NumPy. I have NumPy 1.6 installed with Python 3.2.
Ml SlideMachine Learning: Linear Algebra Reviews Lecture3
First, open an interactive Python shell and load NumPy:

>>> from numpy import *

Create a Matrix

Let's try just creating the 4x2 matrix he shows in slides 2 and 3.The basic form for creating arrays is to use the array method with parenthesis:

a = array()

Inside the parenthesis, nest some square brackets, and in those brackets just put comma-separated lists of elements in more square brackets. Each square bracket represents a row. Make sure they all have the same number of columns. So to create the 4x2 matrix from slides 2 and 3 use the following

>>> a = array([[1402,191],[1371,821],[949,1437]])

>>> print(a)
[[1402  191]
 [1371  821]
 [ 949 1437]] 

Matrix Elements

Next he talks about accessing the matrix element-wise, so a11 should access '1402.'

>>> a[1,1]
821

Instead of returning the first row of the first column, it gave us the second row of the second column. This is because NumPy is using 0-indexing (start at 0) instead of 1-indexing (start at 1). So to get the first row of the first column we index from 0:

>>> a[0,0]
1402

Matrix Addition

Next let's create two 3x2 matrices and add them together. First, instantiate the matrices:

>>> b = array([[1,0],[2,5],[3,1]])
>>> c = array([[4,0.5],[2,5],[0,1]])
>>> print(b)
[[1 0]
 [2 5]
 [3 1]]
>>> print(c)
[[ 4.   0.5]
 [ 2.   5. ]
 [ 0.   1. ]]

Add b and c and see what happens:

>>> b+c
array([[  5. ,   0.5],
       [  4. ,  10. ],
       [  3. ,   2. ]])

As you can see, NumPy correctly performed an element-wise addition.

Scalar Multiplication

Next, multiply a scalar by a 3x2 matrix. We'll use matrix 'b' from above:

>>> 3 * b
array([[ 3,  0],
       [ 6, 15],
       [ 9,  3]], dtype=int32)

Again, NumPy correctly multiplied each element of the matrix by 3. Note that this can be done in any order (i.e. scalar * matrix = matrix * scalar).
Division is just multiplication by a fraction:

>>> d / 4
array([[ 1.  ,  0.  ],
       [ 1.5 ,  0.75]])

Combination of Operands

Order of operations is important. In this slide the instructor sets up three vectors (3x1) and provides an example in which he multiples by a scalar then adds then subtracts then divides. Let NumPy do it and see what happens:

>>> e = array([[1],[4],[2]])
>>> f = array([[0],[0],[5]])
>>> g = array([[3],[0],[2]])
>>> 3 * e + f - g / 3
array([[  2.        ],
       [ 12.        ],
       [ 10.33333333]])

As before, NumPy produces the same answer as the instructor found by doing it by hand.

Matrix-vector multiplication

We can multiply a matrix by a vector as long as the number of columns of the matrix is the same as the number of rows of the vector. In other words, the matrix must be as wide as the vector is long.

>>> h = array([[1,3],[4,0],[2,1]]) # 3x2
>>> i = array([[1],[5]]) # 2x1
>>> h * i
Traceback (most recent call last):
  File "<pyshell#54>", line 1, in <module>
    h * i
ValueError: operands could not be broadcast together with shapes (3,2) (2,1)

My multiplication operation didn't work, and it helpfully gave me the shapes of the arrays for which multiplication failed. This is because the multiplication operator '*' causes element-wise multiplication. For that to work, the matrices need to be of the same shape (hence the error message shosed us the shapes were different). What we want is the dot product.

The Dot Product

To multiply two matrices using the dot product, use the dot() method:

>>> h = array([[1,3],[4,0],[2,1]]) # 3x2
>>> i = array([[1],[5]]) # 2x1

>>> dot(h,i)
array([[16],
       [ 4],
       [ 7]])

That gives us the correct answer, according to the slides. A longer example:

>>> j = array([[1,2,1,5],[0,3,0,4],[-1,-2,0,0]])
>>> k = array([[1],[3],[2],[1]])
>>> dot(j,k)
array([[14],
       [13],
       [-7]])

This one checks out with the slides as well.

Matrix-Matrix Multiplication

As far as NumPy is concerned, matrix-matrix multiplication is just like matrix-vector multiplication.

>>> l = array([[1,3,2],[4,0,1]])
>>> m = array([[1,3],[0,1],[5,2]])
>>> dot(l,m)
array([[11, 10],
       [ 9, 14]])

He goes on to give one more example with a pair of 2x2 matrices:

>>> n = array([[1,3],[2,5]]) # 2x2
>>> o = array([[0,1],[3,2]]) # 2x2
>>> dot(n,o)
array([[ 9,  7],
       [15, 12]])

Next he provides some context by using the example of home prices. He sets up a 4x2 and 2x3 matrix and multiplies them to quickly come up with price predictions:

>>> p = array([[1,2104],[1,1416],[1,1534],[1,852]]) # 4x2
>>> q = array([[-40,200,-150],[0.25,0.1,0.4]]) # 2x3
>>> dot(p,q)
array([[ 486. ,  410.4,  691.6],
       [ 314. ,  341.6,  416.4],
       [ 343.5,  353.4,  463.6],
       [ 173. ,  285.2,  190.8]])

Matrix Multiplication Properties

Show that matrix multiplication is not commutative:

>>> A = array([[1,1],[0,0]]) # 2x2
>>> B = array([[0,0],[2,0]]) # 2x2
>>> dot(A,B)
array([[2, 0],
       [0, 0]])
>>> dot(B,A)
array([[0, 0],
       [2, 2]])

To test this another way I asserted equality between the two operations and found a neat element-wise comparison:

>>> dot(A,B) == dot(B,A)
array([[False,  True],
       [False, False]], dtype=bool)

To show the associative property of arrays, create another 2x2 array C and multiply them in different groupings, but in the same order, to show that the result is always the same:

>>> C = array([[1,3],[0,2]]) # 2x2
>>> A = array([[1,1],[0,0]]) # 2x2
>>> B = array([[0,0],[2,0]]) # 2x2
>>> C = array([[1,3],[0,2]]) # 2x2
>>> dot(A,dot(B,C))
array([[2, 6],
       [0, 0]])
>>> dot(dot(A,B),C)
array([[2, 6],
       [0, 0]])

Identity Matrix

NumPy comes with a built-in function for producing an identity matrix. Just pass it the dimension (numnber of rows or columns) as the argument. Optionally tell it to output elements as integers in order to clean up the output:

>>> identity(3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> identity(3, dtype=int)
array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])

Show that for any matrix A, AI=IA=A. We can use the same A from before:

>>> A = array([[4,2,1],[4,8,3],[1,1,0]]) # 3x3
>>> I = identity(3, dtype=int)
>>> dot(A,I)
array([[4, 2, 1],
       [4, 8, 3],
       [1, 1, 0]])
>>> A
array([[4, 2, 1],
       [4, 8, 3],
       [1, 1, 0]])
>>> dot(A,I)==dot(I,A)
array([[ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True]], dtype=bool)
>>> dot(A,I) == A
array([[ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True]], dtype=bool)

Inverse and Transpose

Show that if A is an mxm matrix, and if it has an inverse, then A(A^-1) = (A^-1)A = I. To do this, we can use the same 3x3 matrix A from above:

>>> A = array([[4,2,1],[4,8,3],[1,1,0]]) # 3x3

>>> inv(A)
Traceback (most recent call last):
  File "<pyshell#112>", line 1, in <module>
    inv(A)
NameError: name 'inv' is not defined

We got this error because though we loaded NumPy, we need to also load the special linear algebra library:

>>> from numpy.linalg import *

Then we can try the inversion again:

>>> inv(A)
array([[ 0.3, -0.1,  0.2],
       [-0.3,  0.1,  0.8],
       [ 0.4,  0.2, -2.4]])

Now show that that A(A^-1) = (A^-1)A = I:

>>> dot(A, inv(A))
array([[  1.00000000e+00,  -2.77555756e-17,   0.00000000e+00],
       [ -2.22044605e-16,   1.00000000e+00,   0.00000000e+00],
       [ -5.55111512e-17,  -1.38777878e-17,   1.00000000e+00]])
>>> dot(inv(A), A)
array([[  1.00000000e+00,   0.00000000e+00,  -5.55111512e-17],
       [ -2.22044605e-16,   1.00000000e+00,  -5.55111512e-17],
       [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00]])

Note that this was meant to return the identity matrix, but that the float operations returned not zeros but numbers very close to zero. Note also that a matrix of all zeros has no inverse:

>>> C = array([[0,0],[0,0]])
>>> C
array([[0, 0],
       [0, 0]])
>>> inv(C)
Traceback (most recent call last):
  File "<pyshell#121>", line 1, in <module>
    inv(C)
  File "C:\Python32\lib\site-packages\numpy\linalg\linalg.py", line 445, in inv
    return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
  File "C:\Python32\lib\site-packages\numpy\linalg\linalg.py", line 328, in solve
    raise LinAlgError('Singular matrix')
numpy.linalg.linalg.LinAlgError: Singular matrix

The error message agrees with the machine learning instructor, who calls these special matrices "singular" or "degenerate."

Matrix Transpose

Lastly, show some matrix transposition, whereby the rows are flipped element-wise:

>>> A = array([[1,2,0],[3,5,9]])
>>> A.transpose()
array([[1, 3],
       [2, 5],
       [0, 9]])

If we name the transposed array 'B', then we can show that the ijth element of A is the jith element of B:

>>> B = A.transpose()
>>> A[0,2] == B[2,0]
True

>>> A[1,2] == B[2,1]
True 

Conclusion

We used NumPy and NumPy's library linalg to go through the linear algebra review slides from Coursera's Machine Learning course.

Installing NumPy for Python 3 in Windows 7

It's time to do some scientific computing, which, in the Python world, means using NumPy. I'm in a Windows 7 (64-bit) environment running Python 3.2.3 (64-bit).

Getting NumPy installed for Python 2 or Python 3 in Ubuntu was easy. Getting it to work in Windows turned out to be more tricky.

The Short Takeaway

• In a Windows 7 environment (even a 64-bit Windows 7 environment), you must install the 32-bit version of Python 3. The 64-bit version will not work with NumPy 1.6. 
• Furthermore, the 32-bit version of Python 3 must be installed 'just for me', and not 'for everyone on this computer'. 
• Finally, make sure you select the proper NumPy version (for Python 3.2), not the default version from SourceForge (which is for Python 2.6).

In this post I'm assuming you have already installed Python 3 and that you're running Windows 7. Specifically, I'm running Windows 7 Professional, 64-bit, Service Pack 1. What follows is the whole story of the troubleshooting, in case it helps out anyone else having the same issues.

Step 1: Ensure that NumPy Isn't Already Installed

Open a Python prompt and ensure that you don't already have NumPy installed:

>>> from numpy import *
Traceback (most recent call last):
  File "<pyshell#2>", line 1, in <module>
    from numpy import *
ImportError: No module named numpy

Step 2: Download NumPy from SourceForge

The SourceForge homepage for NumPy can be found at http://sourceforge.net/projects/numpy/. From there find the latest version of NumPy for Windows, and download it to your default download location. For me as write this, that means downloading NumPy 1.6.1. Note that the default SourceForge download link pointed to a version of NumPy compiled for Python 2.6. I had to navigate a bit deeper to find a version for Python 3.

Step 3: Open the installation file

Double-click the file from wherever you downloaded it to to start the installation wizard. It's trustworthy open-source software, so it's safe to click through all the prompts and allow it to be installed.


This leads us to a crucial error:

It reads: Python version 3.2 required, which was not found in the registry. Yikes. I definitely have Python 3.2, but it's telling me it looked in the Windows registry and couldn't find it. I heard somewhere that installing Python for 'just this user (me)' instead of 'for all users of this computer' is one way to get around this.

Step 4: Reinstall Python 'just for me'

So to comply, I uninstall Python, and reinstall it, being careful to install it 'just for me' as shown below.

Step 5: Try installing NumPy again

With Python 3.2.3(64-bit) reinstalled properly, I try the NumPy installer again. This time it finds Python in the registry (presumably), and installs NumPy 1.6 without issue. Now test it out.

Step 6: Test out NumPy

To make sure it installed correctly, go into the Python interpreter and try importing NumPy:

from numpy import *

This returns the following mess:

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    from numpy import *
  File "C:\Python32\lib\site-packages\numpy\__init__.py", line 137, in <module>
    from . import add_newdocs
  File "C:\Python32\lib\site-packages\numpy\add_newdocs.py", line 9, in <module>
    from numpy.lib import add_newdoc
  File "C:\Python32\lib\site-packages\numpy\lib\__init__.py", line 4, in <module>
    from .type_check import *
  File "C:\Python32\lib\site-packages\numpy\lib\type_check.py", line 8, in <module>
    import numpy.core.numeric as _nx
  File "C:\Python32\lib\site-packages\numpy\core\__init__.py", line 5, in <module>
    from . import multiarray
ImportError: DLL load failed: %1 is not a valid Win32 application

The last line is telling: "not a valid Win32 application." Some of the online forums seem to suggest that it could be a problem with 64-bit Python. Here's the exact version I'm running: "Python 3.2.3 (default, Apr 11 2012, 07:12:16) [MSC v.1500 64 bit (AMD64)] on win32".

 Step 7: Uninstall 64-bit Python, install 32-bit Python

So I'm going to uninstall NumPy and uninstall this 64-bit version of Python 3.2.3, and in its place install a 32-bit version of Python 3.2.3. Again, be careful to install 'just for me.' When this is done, I have this version installed: "Python 3.2.3 (default, Apr 11 2012, 07:15:24) [MSC v.1500 32 bit (Intel)] on win32". Now try NumPy again.

Step 8: Try installing NumPy again

Using the same NumPy binary as every time before, re-install it. To re-cap, I'm installing this in the context of a 32-bit Python 3.2.3 installation. I get no errors from the installation of NumPy, so it's time to test it.

Step 9: Test NumPy

One more time, in the Python interpreter, try importing NumPy:

>>> from numpy import *

This time it returns nothing, meaning it worked! Try creating a NumPy array, and see if it returns the proper type:

type(array([1,6,3,7]))

<class 'numpy.ndarray'>

That worked too, which means our task of installing NumPy for Python 3 in Windows 7 has been completed.

Conclusion

Here are the condensed steps for getting NumPy to work with Python 3 in Windows 7

1. Regardless of whether you have a 32-bit or a 64-bit operating system, install the 32-bit version of Python 3.2
2. Make sure you have installed Python 3.2 'just for me', and not 'for all users of this computer'.
3. Make sure you download the correct version of NumPy from SourceForge, not the default that it offers as the latest version (which is for Python 2.6 instead of Python 3)

Installing NumPy 1.6 for Python 3 in Ubuntu 12.04

I have a fresh install of Ubuntu 12.04. I want to use NumPy with Python 3. Note that Ubuntu 12.04 ships with Python 2.7 under the 'python' namespace and also ships with Python 3 under the 'python3' namespace.

Go to terminal, check python version

$ python --version
Python 2.7.3
$ python3 --version
Python 3.2.3

We want to use Python 3, so for the rest of this tutorial, make sure you're using python3 and not just python.

Go into the Python 3 interpreter:

$ python3
Python 3.2.3 (default, Apr 12 2012, 21:55:50) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> 

See if you have NumPy already installed:

>>> from numpy import *
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: No module named numpy

Nope, it's not there, so we need to install. Exit the Python interpreter.

>>> exit()

Install NumPy using apt-get

Instead of browsing web forums and finding the right package to download and compile, just let Ubuntu's built-in package manager, apt-get, do the work. This is just like installing NumPy for the standard python installation, but because we need NumPy for the Python3 installation, we'll affix the '3' where needed:

$ sudo apt-get install python3-numpy python3-scipy

Let it have all the dependencies it wants (just type 'yes' when prompted).

See if it worked

To see if it worked, go back into the Python interpreter, import NumPy, create a NumPy array, and make sure it's a NumPy array:

>>> from numpy import *
>>> type(array([1,2,4,5]))
<type 'numpy.ndarray'>

Yup it worked.

Conclusion

In Ubuntu 12.04, just use apt-get to install NumPy for python 3 using

$ sudo apt-get install python3-numpy python3-scipy

Everything just works.

Installing NumPy 1.6 on Python 2.7 in Ubuntu 12.04

I have a fresh install of Ubuntu 12.04. I want to use NumPy, and I'm okay with using the default Python version that ships with Ubuntu 12.04, which is Python 2.7. Note that Ubuntu 12.04 also ships with Python 3 under the 'python3' namespace. 

Go to terminal, check python version

$ python --version
Python 2.7.3

Go into the Python interpreter:

$ python
Python 2.7.3 (default, Apr 20 2012, 22:44:07) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.

See if you have NumPy already installed:

>>> import numpy
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: No module named numpy

Nope, it's not there, so we need to install. Exit the Python interpreter.

>>> exit()

Install NumPy using apt-get

Instead of browsing web forums and finding the right package to download and compile, just let Ubuntu's built-in package manager, apt-get, do the work:

$ sudo apt-get install python-numpy python-scipy

Let it have all the dependencies it wants (just type 'yes' when prompted).

See if it worked

To see if it worked, go back into the Python interpreter, import NumPy, create a NumPy array, and make sure it's a NumPy array:

$ python

>>> from numpy import *

>>> type(array([1,2,4,5]))
<type 'numpy.ndarray'>

Yup it worked.

Conclusion

In Ubuntu 12.04, just use apt-get to install NumPy using

$ sudo apt-get install python-numpy python-scipy

Everything just works.