Going from numbers to pictures - visible monte carlo!

One of the best things about math is that it's so... visible. Consider the monte carlo estimation we ran in the last lesson: nothing creates "wow factor" like "3.134!"

Well... at least, I almost went "wow." Sort of. (To be perfectly honest, despite knowing what was going on, my actual reaction was "wow... that's it?") Even if we'd have gotten montecarlopi.py to give us a truly accurate value for pi (let's say "3.14159," for the sake of argument) it would have been better than what we had...

But even if we'd have gotten something as accurate as that, it still would have been a number, nothing more, nothing less. It's hardly the sort of thing to write home about, especially if you're already living at home. ("Dear self: today, generated a value for pi. Not even an especially close one, but a value...")

So let's use our mad python skills to make like a turtle and... do whatever turtles do. In this case, they draw. Yes, we're about to start doing graphics. What's more, since my normal process is to build on what we had before, we're going to use python's turtle module to show us how everything goes, as it runs.

First things first: here's a slightly revised montecarlopi.py module, which has a distance function built in because we need to consolidate things. After we show this puppy off, we'll dig into some programming theory - much of which we won't use, incidentally - just to show you how big-brained computer scientist people think sometimes, and why we big-brained computer people sometimes have so few friends.

Here's montecarlopi.py:

import random
import math

def hypotenuse(x, y):
    return math.sqrt(math.pow(x,2)+math.pow(y, 2))

total=0
incircle=0

random.seed(5)
while total<10000000:
  if (triangle.hypotenuse(random.random(), random.random()) <= 1) : incircle+=1
  total=total+1

print (incircle*4.0)/total

With one million samples, our value for pi turns out to be a fairly-accurate 3.1413076, which isn't bad, but isn't good enough to land anyone on the moon or anything. What would be cool, though, is to watch that value as it's generated, to see how it gets closer and closer to the real value of pi.

This sort of thing is horribly popular, because everyone likes to see things as they get generated. Here be dragons, though: it's fun to watch, and also very very slow. In the real world, if we had to watch it, we'd tell it to update our values every few thousand samples, because otherwise we'd grow old and die before it's done. (This is typed, by the way, by a fellow whose beard is grey, who walks with a cane, and who generally looks like he's auditioning for a role as Methuselah, the oldest man in the Bible. Why do I look like this? Because I have a nasty tendency of liking to watch simulations like this.)

The first thing I thought of, when approaching this lesson, was to add something called an observer to my little algorithm, which would - potentially - get called every time a new sample was generated. This observer would be generic, and thus able to do whatever I wanted it to. Here's an observer that tells us the point generated, and how far it is from the center of our circle (called the "origin" by math geeks), in a file called simpleobserver.py:

def observe(x,y,distance):
  print "x =", x,", y =",y,",distance =",distance

We're about to jump into another truly fun bit of programming: data structures.

When you play "Go Fish," you have a bunch of cards in your hand at first: in programming, this is a set or a list. The difference between the two is... well... artificial. Theoretically, a set can have only one of a given value, whereas a list can have more than one, so the cards you have in your hand are a set rather than a list (because if you have more than one nine of hearts in your hand, you're cheating... or playing with a very funky deck.)

What we actually want here is to create a list of observers, so we can observe the values being generated over and over again. We can even generate the value of pi as an observer! In fact, let's do this.

We'll do it by creating a sort of object (a noun!) called an Observer. An observer, for us, is anything that knows how to observe something else - in this case, the things being observed are three values. In simple terms, our Observer is a container for something that has the "observe(x,y,distance)" definition, above.

We then create another object (another noun) called a MonteCarloObserver. This is a type of Observer that will be very, very important for us: it's what we'll use to generate pi!

Lastly, we create a sort of driver for the whole thing, a MonteCarlo object, which creates a MonteCarloObserver for us, and sticks it in an internal list of observers. (It also lets us add other observers to the list.)

We're going to get to some graphics soon - but we really need to understand these concepts first. Let's look at the code, which shows us a lot of new things, and then we'll try to walk through them:

import random

class Observer:
  def __init__(self):
    return

  def observe(self,x,y,distance):
    return

class MonteCarloCalculator(Observer):
  def __init__(self):
    self.total=0
    self.incircle=0

  def observe(self,x,y,distance):
    self.total+=1
    if distance < 1: self.incircle+=1

  def pi(self):
    return (self.incircle*4.0)/self.total

class MonteCarlo:          
  def __init__(self):
    self.calc=MonteCarloCalculator()
    self.observers=[self.calc]
  def addlistener(self, ob):
    self.observers.append(ob)
    
  def generate(self, count=1000):
    total=0
    while(total<count):
      x=random.random()
      y=random.random()
      d=math.sqrt(math.pow(x,2)+math.pow(y,2))
      for f in self.observers: f.observe(x,y,d)
      total=total+1
    return self.calc.pi()

The __init__(self) definition is called a "constructor," or "initializer." Whenever we create an object of a given type, that type's "__init__" is called so that we can make sure the object's ready for us to use.

Every time you see a "self" reference, that means that Python will use the object's value - so "self.observers" refers to a list that's contained in each MonteCarlo instance, and self.calc refers to a MonteCarloCalculator that's created for every new MonteCarlo.

The "self.observers" assignment - uses []. That's to inform Python that this is a list. (Otherwise, it thinks that the value is a simple value, and that's exactly what we don't want.)

generate() is special, too - for lots of reasons.

For one thing, the "count=1000" means that we don't have to tell it how many iterations to make - if we don't give it one, it'll use 1000 as a default value.

What generate() does is the heart of everything we're trying to do. It generates a point (up to the number of points we asked it to) and that point's distance from the origin (0,0). This is the basis for the Monte Carlo method, of course, no magic there. But then it does something wierd:

for f in self.observers: f.observe(x,y,d)

This means that for every Observer in our observers list, call that Observer's "observe" method. This allows us to do something special on every call. We can, for example, output a running calculation of pi, or ... draw something on the screen showing us our samples as they're generated. The last thing it does is return the calculation of pi it has.

Using this thing is a little more complicated than it used to be. If we put that code in a file called "mcp2.py," we use it like this:

$ python
Python 2.5.1 (r251:54863, May 18 2007, 16:56:43) 
[GCC 3.4.4 (cygming special, gdc 0.12, using dmd 0.125)] on cygwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import mcp2
>>> cl=mcp2.MonteCarlo()
>>> cl.generate()
3.172
>>> [CTRL-D]

We import the module, then create an instance of the MonteCarlo class, then tell it to generate.

This is a heck of a lot of stuff to go through before we've even drawn one cool line on the screen!

However, we now have the groundwork for everything we need. Basically, all we need to do is create an observer that draws our results in a window, and another observer that shows the value of pi as its generated so we can see how we get closer and closer to what we think pi is.

Let's do the observer that shows us pi over time, just to keep things simple. Let's store it in "runningpi.py":

import mcp2

class RunningTotal(mcp2.MonteCarloCalculator):
  def __init__(self):
    mcp2.MonteCarloCalculator.__init__(self)
    self.counter=0

  def observe(self, x, y, distance):
    mcp2.MonteCarloCalculator.observe(self, x, y, distance)
    self.counter=self.counter+1
    if self.counter == 100:
      print self.pi()
      self.counter=0

Now, we can use this class as in the following example:

$ python
Python 2.5.1 (r251:54863, May 18 2007, 16:56:43) 
[GCC 3.4.4 (cygming special, gdc 0.12, using dmd 0.125)] on cygwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import mcp2
>>> import runningpi
>>> cl=mcp2.MonteCarlo()
>>> cl.addlistener(runningpi.RunningTotal())
>>> cl.generate()
3.04
3.18
3.16
3.15
3.144
3.13333333333
3.10285714286
3.13
3.12444444444
3.128
3.1280000000000001
>>> [CTRL-D]

As you can see, we show a value that slowly creeps towards a better and better value of pi. We're going to keep this class around, and we're going to write yet another observer that draws pretty pictures as we generate pi.

You know, we've got a lot of text here... and still no graphics! (Darn these multiple concepts!)

Finally - let's generate one last file, called "graphicpi.py". We'll use python's built-in turtle graphics:

import mcp2
import turtle

class PiGraph(mcp2.MonteCarloCalculator):
  def __init__(self):
    mcp2.MonteCarloCalculator.__init__(self)
    self.scale=200
    
    # create a graphics area, make the turtle non-slow
    turtle.setup(width=2*self.scale+10, height=2*self.scale+10)
    turtle.speed('fastest')
    turtle.tracer(0)
    
    # draw a square!
    turtle.down()
    turtle.forward(self.scale)
    turtle.left(90)
    turtle.forward(self.scale)
    turtle.left(90)
    turtle.forward(self.scale)
    turtle.left(90)
    turtle.forward(self.scale)
    turtle.left(90)
    # draw a circle
    turtle.forward(self.scale)
    turtle.left(90)
    turtle.circle(self.scale, 90)
    turtle.up()
    turtle.goto(0,0)

  def observe(self, x, y, distance):
    turtle.up()
    turtle.goto(x*self.scale, y*self.scale)
    turtle.color(0,0,1)
    if distance>1 : turtle.color(1,0,0)
    turtle.down()
    turtle.forward(1)
    turtle.up()

Now, let's run our simulation, with the following commands, and see what you get - and see if you can figure out how and why.

$ python
Python 2.5.1 (r251:54863, May 18 2007, 16:56:43) 
[GCC 3.4.4 (cygming special, gdc 0.12, using dmd 0.125)] on cygwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import mcp2
>>> import runningpi
>>> import graphicpi
>>> cl=mcp2.MonteCarlo()
>>> cl.addlistener(runningpi.RunningTotal())
>>> cl.addlistener(graphicpi.PiGraph())
>>> cl.generate(10000)

It'll take a while, but it's fun to watch! Use larger values for more accuracy and a longer - in some cases, much longer - runtime.