Sunday, December 26, 2010

Parsing the contest results with Python

Reading and generating data files is one of the most common tasks I do with Python. Recently, I used python to see what my record was against individual players in the recent AI contest.

First, I grabbed the my results starting here. Then, I grabbed the top 300 players. Using some tricks learned from this book, I then parsed out the parts of the data that I cared about.

Before going any further, what I am about to show is a bunch of hacks that wouldn't work in a "serious" application. There is no error checking of any kind because I intend to only run this code myself from the command line a few times, and then never see it again.

This isn't parsing like what is needed in a compiler. This is messy, quick parsing. The data is messy, so the method seems appropriate. Coming up with formal grammar is probably overkill for this use.

Here is what the game result data looks like, after removing some junk at the top of the file, and then formatting it to fit the page.

< -- stuff from top of file.... >
<tr class="odd">
<td>Dec 1st 17:00:51</td>    
<td><a href="profile.php?user_id=3938">dabino</a></td>    
<td class="game_draw">Draw</td>    
<td><a href="visualizer.php?game_id=9559407">View Game &gt;&gt;</a>
</td>  </tr>  <tr class="even">    
<td>Dec 1st 16:52:57</td>    
<td><a href="profile.php?user_id=8737">Hazard</a></td>    
<td class="game_win">Win</td>    
<td><a href="visualizer.php?game_id=9558701">View Game &gt;&gt;</a>
</td>  </tr>  <tr class="odd">    
<td>Dec 1st 16:48:30</td>    
<td><a href="profile.php?user_id=7196">FlameN</a></td>    
<td class="game_win">Win</td>    
<td><a href="visualizer.php?game_id=9558244">View Game &gt;&gt;</a>
  </td>  </tr>  <tr class="even">    
<td>Dec 1st 16:42:13</td>    
<td><a href="profile.php?user_id=10663">smloh1</a></td>    
<td class="game_win">Win</td>    
<td><a href="visualizer.php?game_id=9557675">View Game &gt;&gt;</a></td>  </tr>
< -- stuff from bottom of file.... >

All of that, and more, was in one big line in the file. If you look, you can see that all of the interesting parts on are separated by tr tags. So, we can use the string.split method to break the file up into chunks that contain individual results.

def read_result_file(filename):
    text = open(filename).read()
    return text.split("<tr")

I'm going to mention this one more time… there is no error checking in here. If the file doesn't exist, I'll just get the error message at the prompt.

From each line, I want to collect the user id of my opponent, and the result. The user id is contained in a string like "user_id=12345", and the result is in a string like 'class="game_win"'. Both pieces of data occur in one line of text, so I will use a regular expression to try to pull it out, and test that regular expression against each "line" of text from the read_result_file function.

Here is a regular expression to get what I want out of the string. The regex is compiled with the VERBOSE flag, so that I can insert white space and comments into the string to make it more readable.

import re

data_re_string = r'''#
# find the user_id
(\d+)"   # the id is here
.*   # ignore all of this
<td\ class=

data_re = re.compile(data_re_string, re.VERBOSE)

Now, we can search all the lines to find matches. Not all of the lines will have a result. For example, the first line has everything before the first <tr> flag. So, the loop will have to handle that.

def pull_games(rows):
    # run the regular expression on every line
    tmp = ( for line in rows)

    # only keep the results where there is a match.
    # The search function returns None if there isn't a match.
    result = [m.groups() for m in tmp if m]
    return result

Now, we can see if that works.

>>> games = pull_games(read_result_file("games_1.html"))
>>> import pprint
>>> pprint.pprint(games[:10])
[('3938', 'draw'),
 ('8737', 'win'),
 ('7196', 'win'),
 ('10663', 'win'),
 ('7526', 'draw'),
 ('11610', 'loss'),
 ('11985', 'win'),
 ('10464', 'win'),
 ('9325', 'loss'),
 ('6441', 'win')]

If we give it a bad file, the exception is printed at the prompt.

>>> lines = read_result_file("bad_file")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/tmp/", line 4, in read_result_file
IOError: [Errno 2] No such file or directory: 'bad_file'

My game results were spread over three files, so it would be nice to have a helper function to parse all three at once. Here we go…

def get_all_results(files):
    result = []
    for f in files:
    return result

And I can call it easily like this:

>>> # easily create the file names ...
>>> print ["games_%d.html" % x for x in [1,2,3]]
['games_1.html', 'games_2.html', 'games_3.html']

>>> games = get_all_results(["games_%d.html" % x for x in [1,2,3]])
>>> len(games)

>>> pprint.pprint(games[-1])
('14085', 'win')

Let's look at the rankings files. After a bunch of junk, the html looks this (again edited for the page…)

<tbody>  <tr class="odd">
    <td><a href="profile.php?user_id=8565">bocsimacko</a></td>
    <td><a href="country_profile.php?country_id=106">
   <img src="flags/hu.png" />
</a></td>    <td>
<a href="organization_profile.php?org_id=0">Other</a>
</td>    <td><a href="language_profile.php?lang=Lisp">Lisp</a>
</td>    <td>3765</td>  </tr>
  <tr class="even">
    <td><a href="profile.php?user_id=7026">_iouri_</a></td>
    <td><a href="country_profile.php?country_id=2">
<img src="flags/ca.png" />
</a></td>    <td><a href="organization_profile.php?org_id=0">Other</a>
</td>    <td><a href="language_profile.php?lang=C%2B%2B">C++</a>
</td>    <td>3565</td>  </tr>

In this case, the data chunks span multiple lines, so it is not as easy to write a regular expression to just grab everything at once. This is where I turn to what I learned in the text processing book, and write a mini state machine to handle the processing.

(ok, a multiline regex would probably work here too…, but I took this approach instead)

First, I need two regular expressions. One to match against a "ranking" line, and another to match a user_id line. Then, the main function is basically a loop. At each line, it either looks for a ranking, or looks for a user_id, based on what the current state of the machine is. There is a third major state where the code looks for the start of a user record.

Here is the code…

num_match = re.compile('<td>\\s*(\\d+)\\s*</td>')
player_match = re.compile('profile\\.php\\?user_id=(\\d+)">([^<]+)<')

def pull_rankings(text):
    state = "WAIT"
    results = []

    # This is the ranking of the player we are looking at.
    # The default value is -1 to make it easier to spot an error.
    # (if any player has a ranking of -1, the code is wrong..)
    current_rank = -1

    # The states loop from RANK -> PLAYER -> START -> RANK ...
    # except at the start.  At the start the state is WAIT
    # until the code sees the start of the table.
    for line in text:
        if line.startswith("<tbody>"):
            # This only happens once, but it moves the
            # state machine from WAIT to running
            state = "RANK"
        elif state == "START":
            if line.strip().startswith("<tr"):
                state = "RANK"
        elif state == "RANK":
            num =
            if num:
                current_rank = num.groups()[0]
                state = "PLAYER"
        elif state == "PLAYER":
            m =
            if m:
                pid, name = m.groups()
                # Store the current player, then reset the
                # current rank value to the default.
                results.append((name, pid, current_rank))
                current_rank = -1
                state = "START"
    return results

def pull_all_rankings(files):
    """Another helper to pull multiple files of data"""
    all_values = []
    for f in files:
    return all_values

The result is a list of tuples of the name, user_id, and ranking of the players.

Next, I want a dictionary to lookup a player's information keyed on the user_id.

def make_lookup(all_rankings):
    lookup = {}
    for value in all_rankings:
        name, id, ranking = value
        lookup[id] = value
    return lookup

Then, I need to find all of the games against the same player.

def group_by_player(results):
    output = {}
    for result in results:
        pid, r = result
        if pid not in output:
            output[pid] = []
    return output

Finally, I want to create a table of the results against each player, with the game results converted to win/loss/draw records.

def tabulate_games(games):
    result = [0, 0, 0]
    v = dict([("win", 0),
              ("loss", 1),
              ("draw", 2)])

    for g in games:
        result[v[g]] += 1
    return tuple(result)

def create_results(grouped, lookup):
    result = []
    for pid, games in grouped.iteritems():
        if pid in lookup:
            name, ignore, ranking = lookup[pid]
            result.append((int(ranking), name, tabulate_games(games), pid))
    return result

And then I realize that I was lucky to place where I did..

>>> games = get_all_results(["games_%d.html" % x for x in [1,2,3]])
>>> grouped = group_by_player(games)
>>> lookup = make_lookup(pull_all_rankings(
           ["rankings_%d.html" % x for x in [1,2,3]]))
>>> table = create_results(grouped, lookup)
>>> table.sort()  # the first value in each row is the rating

>>> pprint.pprint(table[:10])
[(1, 'bocsimacko', (0, 9, 0), '8565'),
 (2, '_iouri_', (4, 4, 1), '7026'),
 (3, 'Slin-.-', (4, 5, 0), '11248'),
 (4, '_Astek_', (4, 4, 1), '12009'),
 (5, 'jimrogerz', (3, 6, 1), '9325'),
 (6, 'Accoun', (3, 6, 0), '7423'),
 (7, 'george', (5, 4, 0), '11173'),
 (8, 'GreenTea', (2, 6, 1), '5955'),
 (9, 'asavis', (7, 3, 0), '8348'),
 (10, 'bix0r4ever', (3, 6, 0), '7950')]

Thursday, December 9, 2010

Rock, Paper, Scissors over email

This is just to get this off my chest. During the AI contest, I wanted to create a P2P way for two players to play a match, and keep it fair. The Planetwars game had simultaneous moves, so I wanted a way that the players could each make a move at the same time without requiring one player to send a move in the clear to the other first.

This was mostly an exercise for fun, but I thought it may be useful for people to play each other online for testing. Since then, several players, McLeopold in particular, have worked on making it easier for people to host a game server of their own. Hosting a server to play against a friend would be far easier than fully implementing my idea.

The method I wanted to use was based on commitment schemes, where each player could commit to their choice before sending it in the clear. I learned about this idea in Bruce Schneier's Applied Cryptography.

The idea would be that the players could negotiate a map choice, then play out the game. On each turn, the players would

  • pick their move
  • create and send a commitment string based on the move
  • wait to receive the opponents string
  • send the cleartext move
  • verify the commitment string
  • repeat

Implementing this might be fun, but I have other projects I want to do first. But, I have talked about this too much to just drop it. So, here is some code to demonstrate.

There are plenty of security problems with this code. Someone could cheat this with enough effort. There are some simple things that would make this more cheat resistant... but lets get to the details.

Rock, Paper, Scissors

I won't explain the rules here. For a given move, the player will send "rock", "paper", or "scissors".

Creating the commitment string

First, the player picks their move. Then, they pick some sort of random key to use to create a hash of their move using HMAC and SHA-1. (Again, this was for fun. Don't assume anything about the suitability of this idea for any purpose.) The player must also pick a new key every turn, and the key must be something that is unique and unpredictable.

Here is some Python code to create the commitment string.

import hashlib
import hmac

def run_hash(key, message):
    return, message, hashlib.sha1).hexdigest()

Then, to create the string, we could do something like this:

>>> move = "rock"
>>> key = "different random key every move"
>>> hash = run_hash(key, move)
>>> print hash

Transmit strings

Nex the players trade commitment strings. After each is received, they then trade the cleartext moves and secret keys. Then, each player can verify that the commitment string matches the secret key and move.

Even more

The same idea could be used to pick the maps to use. Given a set of possible maps that both players agree on, one of them could create a commitment string for each map using a secret key. Then, the player could shuffle this list of strings, and send it to the other player. The other player would pick one of the strings, and by doing so randomly chooses a map. Given the secret key, the second player can verify that nothing funny was going on.

Scheier's book has a short section on "telephone poker" that goes even further. The players can deal out cards randomly without knowing what the other got.

So, I am not pursuing this any more at the moment. I wanted to write a client to play one on one matches, but I think the easy to install tcp server is a more practical idea.

Tuesday, November 30, 2010

Python tips from the Google AI Challenge

I did a few things during the AI contest that made my life easier. Some may be obvious to most people, but I hope they are useful to someone.

Also, a cleaned up version of my bot is at bitbucket.


Python's logging facility is very useful. Here, I only used it to log to a file, overwriting the file each time. This is something that simple functions could do, but why write them if they are built in? Maybe in another language where extra libraries must be used for logging, writing my own would have made more sense.

The bots cannot write to files on the contest server. During the previous contest, I used logging during my local testing, and commented out the logging setup when I submitted my code. This time around, I wrote my bot so that if it is not given any arguments, the way it will be called on the contest server, it doesn't log. With arguments, it uses the first argument as the logfile name.

Here is where that code lives in my bot:

if __name__ == '__main__':
    import psyco
  except ImportError:

    argv = sys.argv[1:]
    if len(argv) > 0:
  except KeyboardInterrupt:
    print 'ctrl-c, leaving ...'
  except EOFError:
    logging.debug("end of input")

The EOFError was needed for to make local testing cleaner. Someone on the forums posted the logging.exception tip, which is great for seeing the exception information in the log file.


The contest server uses one second timeouts per turn. My bot had logic to try to avoid timing out, although it is not very strict. Testing locally, sometimes I wanted to see if the bot was having problems just because it was too slow, so I would tell it to run for five seconds and see if it worked better. This never helped, but it seemed like a good idea.

Also, when replaying parts of games that I had lost, I would run with logging on. With logging on, the bot took longer, and might not do the same thing it had done in the game. So, I would let it take longer than a second to see what it was thinking.


One of the greatest things about this contest is that the games are all based on text. That means we could theoretically write bots in any language that does I/O. (Practically, the admins have limited resources, and can only set up a reasonable number of languages). Also, it means that the state of the game is represented by a string. If we have that string, we can see what our bot is doing. So, fairly early in the contest, I wrote some code to parse the compressed game state available on the web visualizer and from the local game engine. The code wrote out a file with the game state at every turn. Given that, I could easily rerun my bot on a given move.

As an aside, my bot was stateless between moves. It did not carry any information from one move to the next. Given that, I could more easily debug what my bot was doing on a single turn.

So, with the replay information, and logging, I was able to tweak the logic of my bot, and find several bugs.

I have to admit... I don't use debuggers very often. I like writing small functions that don't rely on state elsewhere in the code. Given that, and some logging, I find that I can usually debug pretty quickly. Oh, and I prefer languages with interactive shells, so I don't have to recompile to try different things when testing my code.

But, sometimes debuggers are the right tool. In my case, I had some very ugly code, and logging was not cutting it. I couldn't figure out why it was in the state it was in. So, pdb to the rescue.

The main loop of my bot is here.

def main(turn_time=0.8):
  turn = 0
    map_data = []
    turn += 1
    while True:
      current_line = raw_input()
      if len(current_line) >= 2 and current_line.startswith("go"):
        start_time = time.time()
        end_time = start_time + turn_time
        s = "\n".join(map_data)
        logging.debug("turn %d\n%s\ngo\n", turn, s)
        main_do_turn(map_data, turn, end_time)


The code reads from stdin until it sees the "go" string. Then, it bundles up all of those lines and calls main_do_turn, which is the function that calculates my move for a given turn.

To debug this in the debugger, I wrote another Python script that looked like this.

import MyBot
import time
import logging
import pdb
import sys

def go(filename):
    s = open(filename).readlines()

    MyBot.main_do_turn(s, 1, time.time() + 9e99)'go("test.txt")')

This script set up my logging for me, set the timeout really big, and started my bot in the debugger.

Something similar could be done in the debugger of any language, as long as the bot has a function that can handle a single turn.

Sunday, October 31, 2010

Towers of Hanoi with Four Pegs

Somehow... this post overwrote a previous one. Not sure how that happened...

I saw a link today on John Cook's blog mentioning the Tower of Hanoi with four pegs. This was an extra credit problem in the algorithms class I took in school, so I looked at the article. The problem of finding the optimimum number of steps is still open. Well, there is a possible solution, but no proof that it is correct.

That extra credit problem was one of the more interesting things I worked on at school. I forgot some of the details of the solution, but went to the basement and found the copy I handed in. The solution was ugly 13 year old C code.

Wikipedia has the algorithm that is believed to be the best, so I decided to use it to check my solution. The solution uses dynamic programming, and I am re-learning lisp at the moment, so I solved the problem in lisp.

Here is the code. The result is a two by n array, where the last element in the second row is the final answer.

;; From Graham's "ANSI Common Lisp"
(defun map-int (fn n)
  "Map fn over integers from 0 to n - 1"
  (let ((acc nil))
    (dotimes (i n)
      (push (funcall fn i) acc))
    (nreverse acc)))

(defun hanoi-4 (n)
  "Return the number of moves needed to solve the tower of hanoi
with four poles"
  (let ((memo (make-array (list 2 (1+ n))
                          :initial-element 0)))

    ;; Set up the "three pole" numbers
    (dotimes (i (1+ n))
      (setf (aref memo 0 i) (- (expt 2 i) 1)))

    ;; Do the dynamic programming.
    ;; For each i less than n, find the best solution to using
    ;; four poles.  The best solution is the min of
    ;; 2 * T(k, 4) + T(i-k, 3) for 0 <= k < i.
    ;; T(u, v) is the number of moves needed to move u disks
    ;; using v poles.

    (do ((i 1 (1+ i)))
        ((= i (1+ n)) memo)
      (setf (aref memo 1 i)
            (apply #'min
                    (lambda (j)
                      (+ (* 2 (aref memo 1 j))
                         (aref memo 0 (- i j))))

Running the test I get:

CL-USER> (hanoi-4 36)
#2A((0 1 3 7 15 31 63 127 255 511 1023 2047
     4095 8191 16383 32767 65535 131071
     262143 524287 1048575 2097151 4194303
     8388607 16777215 33554431 67108863
     134217727 268435455 536870911 1073741823
     2147483647 4294967295 8589934591
     17179869183 34359738367 68719476735)
    (0 1 3 5 9 13 17 25 33 41 49 65 81 97
     113 129 161 193 225 257 289 321 385
     449 513 577 641 705 769 897 1025
     1153 1281 1409 1537 1665 1793))

Here, 1793 is the solution for four pegs and 36 discs.

With dynamic programming, this solution is pretty straight forward. My memory is hazy about what I did, but basically, I did the dynamic programming part the hard way by figuring out the recursive steps ahead of time, for each value up to 36. I had a sense of what worked, but no real proof that it would work beyond the value I tried. People talk about reading their own code after a few months... I am not sure if I want to look at this close enough to figure out what exactly I was thinking.

In college, I never understood dynamic programming. We didn't spend much time on it, and it was a short chapter in our textbook. After doing a few TopCoder competitions, I finally got it.

The lisp solution took a little while, because I was trying to remember all sorts of stuff. I am sure that there are much better ways to do this... Recursion with memoization would probably be easier, but I haven't tried it yet.

Saturday, October 2, 2010

Some weaknesses in my bot

My bot in the Google AI contest has a few clear weaknesses.  I am sure it has other ones, but these big ones mask any other problems.  Both are mainly because I ignored some issues when coding, and hoped that what I had would be good enough.  The code was fine for a while, but won't be much longer.

Here are the problems I know exist.

  1. The bot is overly aggressive about defense.   It hardly does any.  Often, it will throw a ton of ships at another planet, and then lose the source planet.  For people in the contest... yes, my bot loses occasionally to RageBot.
  2. The bot is overly conservative about taking neutral planets.  Watching these games, it is usually clear that the bot could have taken over one or several neutrals without a problem.  
  3. The first two issues interact with each other.  Trying to fix either makes the other even worse.
I am sure some people have found other problems, but I cannot see past these big ones.    Feel free to let me know, though, I'd appreciate the help.

Here is a unique game that I lost.  I like how voidptr's bot abandons ship all over the place, and makes better use of his resources.  [Thanks to Naktibalda for his Planetwars pastebin!!]

Tuesday, September 28, 2010

Google AI Challenge strategy and tips

Basic Strategies for the Google AI Contest.

Here is the first post sharing some strategy and programming ideas for the Google AI Contest, If you don't know what it is, check it out!

In this contest, participants write programs to play a two player, turn based game against other programs. The game, based on Galcon, is a space battle with planets and ships. The ships are either based at a planet, or are part of a fleet flying from one planet to another. When fleets arrive at a planet that is not owned by the player owning the fleet, a battle ensues. Each planet has a growth rate; every turn, the players earn more ships on planets that they occupy.

The great idea with this contest is that the games are turn based, and the input and output are text based. This way, the contest can easily pit programs written in different languages against each other.

I am probably misrepresenting something here, so you may want to browse the official site before going on.

Basic Ideas

Ship economics

There are basically three ways to alter the number of ships you have. First, earn new ships for occupying a planet. Both you and your enemy earn ships based on the planets you each own. Second, attack a neutral planet. That costs you ships, but not your opponent. Third, attack an enemy planet. Then, you each lose the same number of ships. Whoever has the most ships wins, so you want to grow more ships then the enemy, and you want to wisely spend ships to capture neutrals.

One consequence of this is that attacking an enemy planet, without taking it, does not change the "score". Maybe it is worthwhile for other reasons, but it does not change the difference in the ship counts.

Sun Tzu

Maybe other kids were like this... as a teenager I saw a book called "The Art of War", and read it cover to cover in a couple of days... without understanding any of it. I remember one thing out of it though, if only because it seemed weird: attack what cannot be defended, and defend what cannot be attacked. Finally, I have a place where I think I can apply it.

Since you have all of the information of the game state every turn, you can approximate which planets the enemy cannot defend, and try to decide which are the best ones to take. In reverse, you can figure out what is lost, and spend your resources elsewhere instead of trying to save bad planets.

Look into the future

The game state includes a list of fleets that are en route. You can simulate the future of the game, assuming no more fleets are launched, and see what will happen. This way, you can try to occupy a neutral planet the turn after the enemy takes it, and let them absorb the cost of battling the neutral ships. You can also see the effects of the other ships in flight... maybe the enemy is reinforcing a planet you thought you could take, or maybe they gave up on a planet that you were trying to attack.

Cloud the future for the enemy

If you want to see the future for yourself, you would like to keep the enemy from seeing it. One way to do that in PlanetWars is to delay sending a fleet until the last moment. Once a fleet is launched, it is committed, and the enemy will know where it is going and when it will get there. As long as ships are at a planet, they could go any direction on the next turn. Of course, if they are at a planet, they aren't attacking.

There may be some benefit to streaming ships at an enemy planet at times... but I think in most instances, waiting is better.


Obeleh goes into this in more detail, but here are some lessons I have learned, and relearned, over the years.


There are many chances to go crazy with fancy data structures and algorithms here, but in many cases, simpler code will work fine.

You have to spend time to save time.

This contest has two months to go at this point. Write some tools, spend the time to set up logging when you test locally, write some tests, get things working for debugging your bot. This will pay off in the end.

Keep the code simple

This contest is for fun. The more complex the code is, the more likely it is to be wrong. I have thrown away several days work so far because it was too unwieldy.

Throw away code should be thrown away

This advice is for me. Since this contest is for fun, I used a bunch of single character variable names in an early version of my bot. In a dynamic language. I ended up with a lot of weird errors too. The moral is, if you write it like throw away code... make sure you get rid of it. Write code you are going to keep like code you want to keep.

Version Control

Please use some sort of version control. That way you can more easily track changes, undo mistakes, etc. I use Mercurial which is a good fit for keeping a directory of files under version control. But, since the repository is in the same directory as the working files, make sure it is all backed up somewhere else.

Until next time

OK, that was less organized then I intended, but it was a start. Other things I would like to go into in the future are:

  • more details of strategies I have used
  • algorithm ideas (esp. places where complicated algorithms are likely to get used, but overkill)
  • python specific implementation details.
  • some useful tools for testing, debugging, etc.
  • maybe a simple bot to demonstrate some ideas...

If you have problems, the contest forum and irc channel are a good place to go.

Friday, September 24, 2010

Well, I have been doing ok in the Google AI Contest so far as dmj111. As time goes on, I hope to share some of what I learned in the contest.

Monday, September 13, 2010

Google AI Challenge

Well, there goes all of my free time. The Google AI Challenge is back, and I have a hard time ignoring it. My entry is doing well so far, but that will change when the match scheduling gets better. Several people have much better entries, but from the luck of the draw, they have lower ratings.

Last time, I did ok, but spent too much time on it, mostly chasing lack-of-sleep induced bugs. Let's see if I can learn my lesson this time around.

Saturday, September 4, 2010

One to nine puzzle

The problem

Find a number using the digits 1 through 9 once each, such that

  • the number is divisible by 9
  • if the right most digit is removed, the number is divisible by 8
  • if the next right most digit is removed, the number is divisible by 7
  • ... until there is only one digit left, which is divisible by 1

Brute force

The first thing to do, is see if brute force will work. Often it will, and even if it doesn't solve the problem, thinking this approach through can help find a better solution.

We could loop throuch each possible nine digit number, and see if it meets all the rules. We might do it like this:

for i in 123456789 to 987654321 do
    # check if i contains each digit once
    # check if i is divisible by 9
    # check if floor (i/10) is divisible by 8
    # ...

That might not be too bad, but it will check 864,197,532 numbers, and the logic in the loop isn't that clean. Instead of checking if each digit is used once, we could generate all permutations of 123456789 instead. That way, every number meets the criteria of using each digit once. Also, we would only check 9!, or 362,880 values, checking far fewer numbers than the first idea.

Not so brute force

If we look closer at the problem we can see that we have a number of the form abcdefghi. "a" must be divisible by 1. "ab" must be divisible by 2. "abc" must be divisible by 3, and so on. Using that information, we could first check "ab", before appending any more digits to it, and wipe out a lot of possible numbers. For example, we know that "21" can never be the start of our answer, since 21 is not divisible by 2. Of the 72 possible two digit starting values, about half will fail this test. That will take us down to 180,000 values or so! Then, we can eliminate about 2/3's of the rest by checking all of the first three digits, and go down again to about 60,000 possible values.

We can get this solution to work with a recursive function. The inputs will be a list of digits left to use, and the current number we are considering (the first n digits). We will write the function so that the current number is always valid, so it contains no duplicate digits, and it follows the divisibility rules. Also, the digits to use list contains exactly the digits we need to make a valid solution, and no more.

def fast(nums=range(1,10),acc=0,d=1):

    if not nums:
        # we have reached the end... print the answer!
        print "result: ", acc
        for n in nums:
            # update the number with a new digit at the end.
            newnum = n + 10*acc
            if newnum % d == 0:
                # this is a valid partial solution.
                # Call recusively with a new list of nums.
                tmp = [x for x in nums if x is not n]

I abused the default arguments of the function so it can be called without any arguments. The new function runs in less than two milliseconds, which is about 1/300,000th of the brute force timing.

If you only want the answer, the brute force solution is good enough. But, you can work at it to find a cleaner answer, and then learn a technique that can help in other problems.

Tuesday, August 31, 2010

Testing with FsCheck

Now, we will try to test the sort function from before using FsCheck

First, we download and build FsCheck. Then, in the F# interactive shell, (or an .fsx file), we run the following commands.

#I "/home/dave/lib/fsharp/fscheck"
#r "FsCheck.dll"

This sets the load path to where I built the FsCheck dll, and references the dll. The FsCheck Quick Start and documentation have a lot more information than I am about to give...

First, we need a property to test. Since this is a sort function, we would like to see that the data is sorted. We need a function that given an array of data, sorts it, and checks if it is sorted. To keep things easier for now, we will use an integer array. Also, I have learned the hard way to make the tests fail first.. so this first version will be a bad test.

let prop_is_sorted (data:int []) =
    let N = Array.length data
    |> Array.forall (fun i ->
            data.[i] >= data.[i-1])

This function takes an array of data, and checks that each element should be after the one before it. But, we didn't sort yet, so it should fail on some data. Let's run it to find out.

> Check.Quick prop_is_sorted;;
Falsifiable, after 3 tests (3 shrinks) (StdGen (544801007,295317223)):
[|1; 0|]
val it : unit = ()

The test framework generated three different test cases. Two passed, and the third failed. Also, the framework shrank the failing test case three times, reducing the size of the data set that caused the failure. Shrinking doesn't always find a smaller data set, but when it does, it narrows the range of possible causes of the problem.

Let's update the test function to call our sort function, before checking if the array is sorted.

let prop_is_sorted (data:int []) =
    let N = Array.length data
    sort data
    |> Array.forall (fun i ->
            data.[i] >= data.[i-1])

And running it gives us...

> Check.Quick prop_is_sorted;;
Ok, passed 100 tests.
val it : unit = ()

That is great! Except we don't know what kind of data really went into the tests... FsCheck provides ways to get some insight into the data going into the test. In this case, the default generators are probably doing a good enough job, but more complex cases will require more care.

A bad sort function could still pass this test. Perhaps a function that given an array of length N, returns an array with the numbers from 1 to N. We should write more tests on other properties of a sort function to make sure that our sort is doing what it is supposed to. In this case, we have a simple test we can use ... regression test against the system sort.

let prop_compare_with_system (data:int[]) =
    let a = Array.sort data // out of place
    sort data // in place
    a = data

This test passes again with 100 successes, and doesn't pass if we give a faulty sort function.

This type of randomized testing is a nice alternative to unit testing. In some cases, it is much less verbose, and it can generate test cases that I would never think of trying myself.

Now that we have some tests in place, we can implement some of the more advanced tweaks to the quicksort algorithm...

Thursday, August 26, 2010

Simple Quicksort in F#

Here is a very simple quicksort function in F#. Quicksort is a generally fast, general purpose sorting algorithm. It works well because it is a divide and conquer algorithm that does all of its work in place.

Quicksort implementations can be cache friendly because it keeps working on smaller and smaller portions of the data set. If the data is in an array, the computer won't have as many cache misses while the algorithm is running. In modern computers, memory bandwidth is one of the limiting factors in how fast things can go.

F# is a multi-paradigm language, which is a good thing for Quicksort. Although I like working with functional code, this code looks a lot nicer when we use mutable state and arrays.

This implementation does the bare minimum needed. The fanciest part of the code is the median-of-three pivot selection. The reason I included this pivot selection is because the rest of the code is simpler if we can assume that the first element of a subarray is smaller than the current pivot. The loops inside the main partition function get a little cleaner.

In an F# program, code must be defined before it is used. So, the main function here is at the bottom of the listing. To get a top down view, read from the bottom up, or vice versa.

The main algorithm consists of two steps (besides checking for empty subarrays)

  1. Pick a pivot element, and partition the sub-array into a chunk less than or equal to the pivot, and a chunk greater than or equal to the pivot.
  2. Recursively sort the subarrays.

The recursion stops when we have a subarray of length 1. In the code here, we have another case for when the subarray has length 2.

There are many other improvements that we can make to the basic algorithm. In the future I may implement some of them, and also, I would like to show how to test this code with FsCheck, a randomized test framework.

let inline swap (data:'a[]) a b =
    let t = data.[a]
    data.[a] <- data.[b]
    data.[b] <- t
/// Sort three values from an array.  Put the smallest
/// at index a, biggest at c, and median at b.    
let inline median_of_three (data:'a[]) a b c =
    if data.[b] < data.[a] then
        swap data a b
    if data.[c] < data.[b] then
        swap data c b
        if data.[b] < data.[a] then
            swap data a b 

/// partition the data from index a to b (inclusive)
/// use the element at position b as the pivot, and
/// assume data.[a] is <= data.[b].
let inline partition (data:'a[]) a b =
    // After this is called, return the new position
    // of the pivot value.  After the call, all values before
    // the pivot position should be <= to the pivot, and all values
    // after should be >= the pivot.
    let mutable i = a - 1
    let mutable j = b
    let pivot = data.[b]
    while i < j do
        i <- i + 1
        while data.[i] < pivot do
            i <- i + 1

        j <- j - 1
        while data.[j] > pivot do
            j <- j - 1
        if i < j then
            swap data i j
    // Now, i is >= j.  data.[i] is >= pivot.
    // data.[j] <= pivot.  So, move the pivot
    // to the middle where it belongs, and move
    // the value at i to the end of the array.
    swap data i b


let inline sort (data:'a[]) =
    let rec loop a b =
        if a >= b then ()
        elif a + 1 = b then
            if data.[a] > data.[b] then swap data a b
            let c = (a+b)/2
            median_of_three data a b c
            let p = partition data a b
            loop a (p-1)
            loop (p+1) b
    loop 0 (data.Length - 1) | CBAMC

Saturday, August 21, 2010

GCD in F#

Now we will write a function to calculate the greatest common divisor of two numbers. The Euclidean Algorithm is an efficient, and ancient, method of finding the GCD.

Euclid's algorithm takes advantage of some properties of the gcd function and the modulus function. Here is a loose explanation of why the algorithm works

gcd(a,b) = gcd(b, a - q*b) for any integer q
a % b = a - q*b from some integer q,
so gcd(a,b) = gcd(b, a % b)

Then, we just stop the algorithm when b is zero. And now for the F# code.

let rec gcd a b =
    if b = 0 then a
        gcd b (a % b)

Aside from the syntax, and checking for the end condition, this looks like the math.

And now, a few tests.

let f a b = 
    printfn "gcd(%d,%d) = %d" a b (gcd a b)
f 10 9
f (13*61*2*2*2*3) (13*2*3*5*7)
f -10 12345

With the output,

gcd(10,9) = 1
gcd(19032,2730) = 78
gcd(-10,12345) = 5

Monday, August 16, 2010

Using the sieve to factor numbers

Previously, we saw an overly simple implementation of the Sieve of Eratosthenes. Without any optimizations, that code is pretty much only good for puzzles and learning. Even for puzzles, it may be too slow. Some of the puzzle problems out there require code to find the prime factorization of numbers. For small numbers, trial division might be ok. Often, the problems require the factorization for every number less than one million. Then... it gets slow. Using the same idea as the Sieve, we can prepare an array that contains a prime factor for every number, or one if the number is prime. Then, to get the full factorization of the number, we can walk the array. Here is the code to create the array:
def divisors(N):
    """return the biggest prime divisor for numbers less than N"""
    result = [1 for x in xrange(N)]
    for d in xrange(2,N):
        if result[d] > 1:  continue
        for j in xrange(d, N, d):
                result[j] = d
    return result
Then, to find the factorization of a number, we can look it up in the array to find a factor, then look up (n/p) in the array to find the next factor, and keep going until we are down to a prime number.
divs = divisors(2000)
def factor(n, divisors=divs):
        result = []
        while n > 1:
                p = divisors[n]
                if p == 1:
                        return result
                while n % p == 0:
                        n /= p
        return result
print 1000, factor(1000)
print 24, factor(24)
print 17*31, factor(17*31)
Often, the factorization is used for things like computer Euler's totient function, and then we don't care to return the factors, just the result of using them. I used to work on this code: | RLRZt

Saturday, August 14, 2010

Sieve of Eratosthenes

Here is a simple implementation of the Sieve of Eratosthenes in Python.

from math import sqrt
def sieve(N):
    """Sieve of Eratosthenes.  Returns primes less than N
    # Set up the array.  2 is the first prime, not 1.
    ar = range(N)
    ar[1] = 0

    limit = 1 + int(sqrt( N))
    for d in xrange(2,limit):
        if ar[d]:
            # d is prime!  mark off the multiples of d.
            # note that we can start with d*d instead of 2*d.
            for j in xrange(d*d,N,d):
                ar[j] = 0
    return [x for x in ar if x]
There are several optimizations that people tend to do when they write this algorithm. But, this gets us most of the way. I keep rewriting this function for Project Euler problems, but I should really keep a library of the functions I use, or better yet, use Sage.

Friday, August 13, 2010

Simple prime test in F#

Here is the same algorithm as the last post , but in F#. There is no way to break out of a for loop in F#, besides throwing an exception, so this code uses recursion instead.
let is_prime n =
    if n < 2 then false
        let limit = float n |> sqrt |> fun x -> x + 1. |> int
        let rec loop d =
            if d = limit then true
            elif n % d = 0 then false
            else loop (d + 1)
        loop 2

> [1..100] |> List.filter is_prime;;
val it : int list =
  [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41;
   43; 47; 53; 59; 61; 67; 71; 73; 79; 83; 89; 97]

Thursday, August 12, 2010

Simple prime test in Python

A simple way to test if a number is prime is to check if the number is divisible by any other number besides one and itself. This can take a lot of time, but for small numbers, it might be ok.
def is_prime(n):
   # assumes n > 0
   if n < 2: return True
   for d in range(2,n):
       if n % d == 0: return False 
   return True
This simple function can be sped up by not checking as many numbers. We don't need to test all divisors up to n-1. The largest possible divisor is n/2, so we could stop there. But, even better, we can show that if there are any divisors, there must be divisor less than the sqrt(n). If p*q = n, then either p or q must be less than or equal to sqrt(n). Otherwise, p*q would be bigger than n. So, we can just loop up to sqrt(n) because we only need to find a divisor of n, not every divisor. Here is the updated code:
import math
def is_prime2(n):
    if n < 2: return False
    for d in xrange(2, int(math.sqrt(n)) + 1):
        if n % d == 0: return False
    return True

>>> [x for x in xrange(40) if is_prime2(x)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]

Finally, we could use the same code to find a divisor of a number. In fact, it finds the smallest divisor. Some of the more advanced primality tests do not determine any of the factors of the number.
def is_prime3(n):
    """returns the smallest divisor of n.  
       returns 1 if the number is prime, or less than 2"""
    if n < 2: return 1
    for d in xrange(2, int(math.sqrt(n)) + 1):
        if n % d == 0: return d
    return 1