Solutions Looking For Problems

Calculate π by throwing darts

2009-03-20T17:25:00.011Z

We can approximate pi using nothing more than random numbers and some simple geometry: in essence, we randomly throw darts at a square board of side r; within the square we inscribe a quadrant of a circle of radius r with its centre at (0, 0). We count all of the 'throws'; if a dart lands within the quadrant, we also count a 'hit'

For a large number of throws, we see that:


        hits    area_of_quadrant
       ------ = ----------------
       throws    area_of_square

Some half-remembered geometry tells us that:


       area_of_quadrant   (1 / 4) * pi * r^2   pi
       ---------------- = ------------------ = --
        area_of_square           r^2            4

Or:


                area_of_quadrant        hits
       pi = 4 * ---------------- = 4 * ------
                 area_of_square        throws

I first solved this problem as an undergraduate sometime in 1994 as part of a Computational Physics module. Using FORTRAN.

The full explanation can be found here


#!/usr/bin/env python

import random
import math

# class representing a single throw of a dart
class Throw:
    def __init__(self):
        # generate two random coordinates and work out how far away we are from
        # the origin
        self.x = random.random()
        self.y = random.random()
        self.distance = self.distance()

    def distance(self):
        # the distance from the origin is the hypotenuse of a right-angled 
        # triangle with sides of length and x and y. Pythagoras told us that:
        #    distance = sqrt((x^2) + (y^2))
        # which looks like this in python
        return math.sqrt(self.x**2 + self.y**2)

    # did we land inside the quadrant?
    def is_a_hit(self):
        return self.distance <= 1.0

# main class
class MonteCarlo:
    def __init__(self):
        # initialise everything
        self.hits = 0
        self.throws = 0
        self.pi = 0

    # this method is called on every throw
    def increment(self, throw):
        # we always clock up another throw
        self.throws += 1
        # and accumulate a hit if we scored
        if throw.is_a_hit():
            self.hits += 1
        # then get a new value of pi
        self.calculate_pi()

    # explanation can be found here: http://icanhaz.com/calculatingpi
    def calculate_pi(self):
        self.pi = 4 * (float(self.hits) / float(self.throws))

    # we use this in determining whether to print our status
    def divides_by(self, number):
        return float(self.throws) % number == 0

    # represent the current state as a string
    def __repr__(self):
        return "Throws: %10d, Hits: %10d, Pi: %10f, Actual Pi: %10f" % (self.throws, self.hits, self.pi, math.pi)

# if we're called on the command line
if __name__ == '__main__':
    # construct a new instance
    m = MonteCarlo()
    # loop forever
    while 1:
        # keep throwing darts
        m.increment(Throw())
        # only print on every 100000th iteration
        if m.divides_by(100000): print m

Get the code here: svn co http://pikesley.homedns.org/svn/code/MonteCarlo

Permutations

2009-02-28T11:18:00.015Z

We want to work out all the possible arrangements of a list of n items. Why? We don't know. But we do know how to do it: seed with a list containing a single list containing a single item: [['a']]. Then step through each of the lists in our list, and make copies, inserting the next item ('b') at each index, so we get [['b', 'a'], ['a', 'b']]. Then we take this list as our seed and go round again, inserting the next item ('c'), so we get [['c', 'b', 'a'], ['b', 'c', 'a'], ['b', 'a', 'c'], ['c', 'a', 'b'], ['a', 'c', 'b'], ['a', 'b', 'c']]. And so on and so on, until we've permuted all of our items. These lists get pretty big pretty quickly (the length of the list is n! so 10 items generates 3628800 combinations, and doing this may bring your computer to its knees).

I first tackled this problem ~10 years ago, but I don't think my solution then was as elegant as this one.


#!/usr/bin/env python

import sys

# class representing a list of permutations
class Unit(list):
# count the total number of Units spawned
count = 0
def __init__(self, seed = None):
    # if this is the first Unit created, we need to prepare the seed
    if Unit.count == 0: seed = [[seed]]
    Unit.count += 1
    # call our superclass. We are just a fancy list
    super(Unit, self).__init__(seed)

# produce a new generation, splicing in the passed in value
def spawn(self, interloper):
    # a temporary working list
    l = []
    # for each of our existing lists
    for u in self:
        # for each index (including length + 1)
        for i in range(0, len(u) + 1):
            # take a copy of the list
            a = Unit(u)
            # insert the new value at this index
            a.insert(i, interloper)
            # and add the modified copy to our working list
            l.append(a)
    # when we're done, return a new Unit seeded with the list we just built
    return Unit(l)

# calculate and show some stats
def stats(self):
    self.items = len(self[0])
    self.combinations = len(self)
    return "%d items, %d permutations" % (self.items, self.combinations)

# simple wrapper class representing the set of items we want to permute
class Inventory(list):
def __init__(self, chars):
    # we are yet another variation on a list
    super(Inventory, self).__init__(chars)

# do we have more items to give?
def has_next(self):
    return len(self) > 0

# shift the first item off of the list
def next(self):
    return self.pop(0)

# wrapper class that actually does the work. We extend the Unit class mainly
# because we want its stats() method
class Permutations(Unit):
def __init__(self, items):
    # we need an inventory
    self.i = Inventory(items)
    # get the first item off of the inventory
    self.u = Unit(self.i.next())
    # now while the inventory has more items to give us
    while self.i.has_next():
        # make another generation of Units fed with the next item
        self.u = self.u.spawn(self.i.next())
    # at the end, we become whatever the final Unit.spawn() gave us
    self.extend(self.u)

# if we're called on the command line
if __name__ == '__main__':
# we will permute the command line arguments
permutables = sys.argv[1:]
# if we got no arguments, default to this
if not permutables:
    permutables = ["life", "the universe", "everything"]
# do the work
p = Permutations(permutables)
# sort the result (a Permutation is just a fancy list, after all)
p.sort()
# and print it
print p
print p.stats()

Get the code here: svn co http://pikesley.homedns.org/svn/code/Combinatrics/tags/Release-2.0

Evolve some text

2009-02-06T23:20:00.007Z

We start with a target string. We then generate ten random candidates each of the same length as the target, and assign each one a score based on how many correct characters it has. The highest-scoring candidate then becomes the father of the next generation: we spawn ten clones of the winner, with the correct members of each retained and the other members randomised. Then we score each of these, and go round and round again until we get to the target string.

Happy Birthday Darwin. And Python rocks.


#!/usr/bin/env python

import sys
import string
import random
import copy

# the string we will work towards
try:
    target = sys.argv[1]
except:
    target = "from so simple a beginning endless forms most beautiful and most wonderful have been, and are being, evolved."
winningscore = len(target)

# get ascii character set
chars = []
for i in range(32, 127):
    chars.append(chr(i))
chars.append("\n")

# how many candidates will we spawn in each generation
selectees = 10

# a list to hold our candidates
candidates = []

# class representing a character in our candidate
class member():
    def __init__(self, char = ""):
        self.character = char
        # fixity will be set if this character is correct
        self.fixity = False

    # set to a new, random character
    def randomize(self):
        # unless this one has been fixed
        if not self.fixity:
            a = chars[int(random.random() * len(chars))]
            self.character = a

    # string representation    
    def __str__(self):
        return self.character

# class representing a generated candidate (it's basically a list)
class candidate(list):
    # we start with a score of 0
    grade = 0
    size = 0

    # create a new list
    def __init__(self, length):
        self.size = length
        # and populate and score it
        self.populate()
        self.score()

    # fill the array with randomness to start
    def populate(self):
        for i in range(0, self.size):
            m = member()
            m.randomize()
            self.append(m)

    # represent as a string (this is an ugly hack, I think)
    def __str__(self):
        a = []
        for i in self:
            a.append(i.__str__())
        return "".join(a)

    # count how many correct members we have
    def score(self):
        self.grade = 0
        for i in range(0, len(self)):
            # if its a match with the target string
            if self[i].character == target[i]:
                # then we fix this member in place and score 1 point
                self[i].fixity = True
                self.grade += 1

    # randomize ourself    
    def shuffle(self):
        for a in self:
            # members with fixity set will not be randomized
            a.randomize()
        # and work out our new score
        self.score()

# find the candidate with the best score
def getwinner():
    bestgrade = 0
    # assume the first one is the winner for now
    winner = candidates[0]

    # now for each candidate
    for c in candidates:
        # if we're better than the best so far
        if c.grade > bestgrade:
            # hold on to the leader
            winner = c
            # and raise the bar
            bestgrade = c.grade
    # this method will pick the first candidate from the set of possible winners

    # if all of our members are correct
    if bestgrade == winningscore:
        # we've solved it!
        global unsolved
        unsolved = 0
    # send back the winner
    return winner

# real live code starts here

# generate an initial set of random candidates
for i in range(0, selectees):
    candidates.append(candidate(len(target)))

# assume we've not solved it yet
unsolved = True

# how may turns do we take?
iterations = 0

# while we have not the answer
while (unsolved):
    iterations += 1
    
    # find our winner
    winner = getwinner()
    # and show it
    print winner

    # set all of our candidates to the winner
    for i in range(0, len(candidates)):
        candidates[i] = copy.copy(winner)
        # then shuffle each one
        candidates[i].shuffle()

# and when we're finished, print some stats
print "%d-character string, %d iterations" % (len(target), iterations)