Tuesday, October 11, 2011

ASCII Circles

While browsing Reddit, I came across this page. It's a fun little task. Here's my solution, with lots of fiddly maths to make the circle look as clean as possible.

import sys
r = int(sys.argv[1])

for y in xrange(-r*1.3,r*1.3,2):
    for x in xrange(-r*1.3,r*1.3):
        sys.stdout.write('#$@%&0*;:,. '[min(abs(r*r - (x*x+y*y))/(r/2), 11)])
    print

Circle of r = 10:

        ,;*0&&&0*;,       
     ,*%$#@%%%%%@#$%*,    
   ,*@#%0;:,...,:;0%#@*,  
  ,0$@0;.         .;0@$0, 
 .*$@0:             :0@$*.
 ,&#%;.             .;%#&,
 ,&#%;.             .;%#&,
 .*$@0:             :0@$*.
  ,0$@0;.         .;0@$0, 
   ,*@#%0;:,...,:;0%#@*,  
     ,*%$#@%%%%%@#$%*,    
        ,;*0&&&0*;,      

r = 30:

                             .,:;*00&&&&&&&00*;:,.                            
                         ,;0&@$##$@@%%%%%%%@@$##$@&0;,                        
                     .;0%$#$%&*;:,,..     ..,,:;*&%$#$%0;.                    
                   :0%#$%0;:.                     .:;0%$#%0:                  
                .;&$#%0;,                             ,;0%#$&;.               
               ;&$$%*,                                   ,*%$$&;              
             ,0@#%*,                                       ,*%#@0,            
            ;&#@0:                                           :0@#&;           
           ;%#%*,                                             ,*%#%;          
          ;%#%;.                                               .;%#%;         
         :&#%*.                                                 .*%#&:        
        ,0$@*,                                                   ,*@$0,       
        ;%#&:                                                     :&#%;       
       ,0$@*,                                                     ,*@$0,      
       :&#%;.                                                     .;%#&:      
       :&#%;                                                       ;%#&:      
       :&#%;                                                       ;%#&:      
       :&#%;.                                                     .;%#&:      
       ,0$@*,                                                     ,*@$0,      
        ;%#&:                                                     :&#%;       
        ,0$@*,                                                   ,*@$0,       
         :&#%*.                                                 .*%#&:        
          ;%#%;.                                               .;%#%;         
           ;%#%*,                                             ,*%#%;          
            ;&#@0:                                           :0@#&;           
             ,0@#%*,                                       ,*%#@0,            
               ;&$$%*,                                   ,*%$$&;              
                .;&$#%0;,                             ,;0%#$&;.               
                   :0%#$%0;:.                     .:;0%$#%0:                  
                     .;0%$#$%&*;:,,..     ..,,:;*&%$#$%0;.                    
                         ,;0&@$##$@@%%%%%%%@@$##$@&0;,                        
                             .,:;*00&&&&&&&00*;:,.                            
                                                                  

The explanation I posted on Reddit (sorry about the formatting, hopefully it's still clear).

I guess I should explain what this does. Like the challenge poster, I use the Pythagorean theorem. I look through each coordinate in a square of size 2.6r * 2.6r around the centre-point. The circle only really exists in the 2r * 2r square around the point, but the extra padding on all sides allows for the anti-aliasing to continue further along the edges.

Mathematically, a circle consists of the coordinates (x, y) that satisfy x2 + y2 = r2, r being the radius of the circle. Any points we encounter that satisfy this should be "coloured in" as dark as possible. Some points do not fall exactly on the circle, but they fall close to the circle: the (absolute) difference between r2 and x2 + y2 is (if I am not mistaken) the distance between (x,y) and the closest point on the circle to (x,y). If the difference is 0, we have the case above. Otherwise, the larger the value, the less it is part of the circle and so the lighter we should 'colour' it.

To colour points, I index the 'colouring' string by a scaled value of this closeness, scaled so the width of the ring remains the same for different values of r.

Thursday, October 6, 2011

Space-efficient Memoization with dicts

Python's not usually considered a good language for memory- or time-intensive tasks, but its extensive standard library and variety of built-ins make it an invaluable tool for algorithmic work.

Consider the task of memoizing a recursive function. Typically, a large d-dimensional array is used as a cache, where d is the number of parameters in the subproblem's state. This is great for bottom-up iterative dynamic programming, where the optimal solution for every state must be evaluated -- every cell of the array will be occupied, so there's no redundant space. However, many problems do not require that every state be examined. Consider a case of the unbounded knapsack problem where you have a knapsack of size 15 and 3 items available of costs 5, 7, and 11. The bottom-up approach will evaluate the optimal solution for states 1, 2, 3, 4, 6, 7, 8 and 9, none of which contribute at all to the final optimal selection of items; after all, there's no way to fill a knapsack of those sizes given those items. The top-down memoized recursion approach, on the other hand, will only evaluate states which can theoretically contribute to the optimal solution (even if they do not).

Despite the overhead created by recursive calls and a much higher proportion of cache misses, the top-down memoized approach is often a lot faster than the bottom-up dynamic programming solution. However, people will still generally use the d dimensional array for caching, even when memoizing. This leads to lots of unused cells and wasted space.

Python gives you another option: instead of an array, use a dictionary of {d-tuple representing state : optimal solution for state}. Dictionaries have extremely fast look-up times, and the amount of memory saved is substantial.

Consider the 'matrix sum' problem (as seen in Project Euler #345), a variant of the classic n-rooks problem. You're given a n x n grid (chessboard?) where each cell has a value associated with it. Your task is to place n rooks on the board such that none of them threaten another (no two rooks lie on the same row or column) and the sum of values in the cells they occupy is maximised.

This is an exponential dynamic programming problem, where the state of a subproblem has two parameters: the current row, and a binary string representing the columns that have already been occupied and cannot have rooks placed in them from now on. I don't want to elaborate too much on the solution because it's an interesting problem that's worth your attention. There are n possible values for the current row and 2n possible values for the binary string, so in total there should be n * 2n states. Let's work with n = 15. The total theoretical number of states should be 491,520, which would be the number of cells in the d-dimensional array.

However, in practice there are many states that are impossible to reach. For example, there is no state where you're on the very first row and yet all columns are already occupied. Similarly, there's no state where you're on the last row and no columns are occupied. It's possible to mathematically calculate the exact number of reachable solutions (I did it once, it's not very fun) but since this is a programming blog, I'd rather just demonstrate a solution that tells you exactly how many states must be cached.

N = 15
m = [map(int, line.split()) for line in open('matrix.txt')]

cache = {}
def dp(r, u):
    if (r,u) not in cache:
        cache[(r,u)] = 0 if r==N else max(dp(r+1, u|(1<<c)) + m[r][c] for c in xrange(N) if not u&(1<<c))
    return cache[(r, u)]

print dp(0,0), len(cache)

Output:

32768

That's right ladies and gentlemen, only 32768 states are reachable. That's 1/15th of the original memory usage. Obviously you could do the same thing in another language, but Python makes it so easy. Have fun implementing your own hash table in C :)

Monday, September 5, 2011

A Practical Use for the 'complex' Type

In a recent programming competition, a question asked competitors to generate certain fractals (Julia sets, to be precise). This task requires understanding of the mathematics behind fractals; specifically, complex numbers. The judges assumed that most competitors would be unfamiliar with complex numbers and how to manipulate them programmatically. Hence, about half of the problem statement documented the different operations that could be done with complex numbers: polar/Cartesian conversions, getting the modulus and the angle, how the arithmetic operators worked, etc.

Luckily, I recalled that Python has built-in support for complex numbers. All I had to do was plug in the given formulae and the problem basically solved itself.

In celebration, here's a short piece of Python code which generates a 200x200 Mandelbrot set in ASCII, using Python's complex number support and the cmath library. cmath is to complex numbers as math is to real numbers: it provides loads of helpful functions you might need when dealing with complex numbers. Here, I use cmath to get the modulus of z.

import math, cmath
RX, RY = (-2.5,1), (-1,1)

# maps a pixel (x,y) into a (x,y) point on the Mandelbrot canvas
def lmap(n, (a,b)): return n/200.0 * (b-a) + a

def value(x,y):
    z, c = 0, complex(lmap(x, RX), lmap(y, RY))
    iters = 0
    while cmath.polar(z)[0] < 4 and iters < 1000:
        z = z*z + c
        iters += 1
    return int(math.log(iters)/math.log(1000) * 10.0)
        
out = open('out.txt', 'w')
for r in xrange(200):
    for c in xrange(200):
        out.write(' .:;*0&%@$#'[value(r,c)])
    out.write('\n')

A small portion of the generated fractal:

Wednesday, August 3, 2011

Anagram Searching -- A Job Interview Question

Occasionally a problem comes along that stuns me in its simplicity and elegance. This one in particular is from the 2006 International Olympiad in Informatics (IOI) under the name Writing. It's a great question to give in a job interview because a bruteforce solution is easy to work out, and some may stop there. However, with a bit more effort, it's possible to develop a solution that improves on the bruteforce by orders of magnitude.

A summary of the problem is as follows:

You're given a parent string and a query string of assorted characters of some alphabet (ASCII works nicely). You must determine how many times the query string -- or a permutation of the query string -- appears in the parent string.

Here is a typical thought-path one may take to the optimal solution. In order to evaluate the efficiency of each algorithm, we'll call the parent string's length n, the query string's length k and the alphabet size s.

1. Bruteforce


Generate all permutations of the query string, then do a naive sub-string count on the parent string for each permutation.

import itertools
sum(parent.count(perm) for perm in itertools.permutations(query))

Complexity: O(k!) for generating the permutations, then O(nk) for a naive substring counting algorithm, giving O(nk * k!).


2. "Crossing off" method


There are lots of common string search/manipulation problems in circulation, and most involve strings that must be treated as ordered sequences; that is, you can't mess around with the order of the characters.

Instead of treating this problem as an extension of those problems with an extra restriction, take advantage of the unordered property.
Here's a new way of looking at the problem: for every consecutive set of k characters in the parent string, is that set of k characters the same as the set of characters in our query string?

We've now introduced a new problem: efficient comparison of character sets.

Let our first string be "goldy" and our second be "dogle". A fairly obvious way of checking that the characters are the same is by using a "crossing-off" method. For each character in "goldy", check if it's in "dogle". If it is, and you haven't encountered it on a previous pass, cross it off. make sure the character you find in the second string not already been crossed off; "aab" should not be considered the same as "abb". If at the end of the process all characters have been crossed off and early termination didn't occur then your strings match.

Complexity: For every consecutive set of k characters in the parent (O(n) sets), for every character in that set, check if it's in the query string and if so, cross if off. This algorithm has a time complexity of O(nk2) which is already a massive improvement over the brute-force above.

3. Sort + compare


We can improve on the cross-off method. A better way to compare two sets is to sort both sets then perform an element-by-element comparison. "goldy" becomes "dgloy", "dogle" becomes "deglo", and a standard string comparison shows that the strings are unequal. Since you're probably in a language with a linearithmic sort in its standard library, this solution is easy to implement and works well.

sortedQuery = sorted(query)
print sum(sorted(parent[i:i+len(query)]) == sortedQuery for i in xrange(len(parent) - len(query) + 1))    

Complexity: You're still going through every consecutive block of k characters in the parent string, but this time you just have to compare the sorted block to the sorted query string. O(nk lg k).

4. Constant-time lookup.


It turns out that sort-based set comparison is suboptimal. We can take advantage of the array structure's constant-time random access to create a lookup table of the form lookup[char] = number of occurences of 'char' in query string, for every value of char in our alphabet. We first generate the lookup table in O(k) time and store it. As usual we'll iterate through every consecutive block of k characters in the parent string, but this time instead of sorting, we populate a second lookup table with the data from this block. To compare the sets all we have to do is compare the lookup tables.

Complexity: For each block, generate lookup tables (O(k)) then compare (O(s)). O(n(s + k)).

There's a slight optimisation to be made. If s is large and k is small, our table comparison will end up going through lots of characters which never appear in the query. To prevent this, we can store the characters we come across when constructing the second lookup table. When we compare the tables, we only have to compare the counts for these characters.

import string

numMatches = 0

# allows for the O(nk) optimisation
queryChars = set(query)

correctTable = [0] * 127
for c in query:
    correctTable[ord(c)] += 1

testTable = [0] * 127

for i in xrange(len(parent)-len(query) + 1):
    thisBlock = parent[i:i+len(query)]
    
    for c in thisBlock:
        testTable[ord(c)] += 1

    for c in queryChars:
        if correctTable[ord(c)] != testTable[ord(c)]:
            break
    else:
        numMatches += 1
        
    for c in thisBlock:
        testTable[ord(c)] -= 1

print numMatches

Complexity: Generate the second lookup table in O(k), compare them in O(k) too. O(nk).

5. Magic


The dominating mind-frame when one thinks about this problem is that it's just a twist on the well-known problem of set comparisons: the solution has to loop over each consecutive k-block in the parent string, and my job is to make the comparison between the block and the query string as fast as possible. The block string/query string comparison is lower-bounded by O(n(s + k)).

But this approach always leaves with you two loops, even though the inner loop covers the exact same content as the outer loop. The optimal solution is obtained by eliminating redundancy in the comparison loop by integrating it with the outer loop.

We have to have an O(n) somewhere, if only for inputting the parent string. Our previous solution was doing two O(k) things on top of that: generating the second lookup table, then comparing the two tables. Say we generate a lookup table for the first block of the parent string, finish doing all our processing for that block, then move on to our next block. Our previous algorithm will generate a brand new lookup table now, even though we've just shifted over by a single character. The only rows that can look different between the two tables are the row of the first character of the previous block (which is now not covered by our new block) and the row of the last character in our new block (which was not previously covered). All the other rows will look exactly the same each time. Now that we're only updating two entries each time, generating the table is now O(k) for the first iteration and O(1) for each of the n-k blocks following: O(n). We just stumbled upon an O(n * min(k,s)) solution (you can derive the complexity from the observations above), but let's go straight on.

Now to apply this logic to comparing the two lookup tables. Previously, it was easy to know whether or not a given block's table matches our query string's table: if all of the character counts were the same, they match. The goal is to only operate on the character our block 'left behind' and the character our block is 'picking up'.

Throughout the program, keep track of a distance (similar to the concept of Hamming distance). distance will contain an integer representing the number of unique characters in our current block's table whose counts match that character's count in the query string's table. This will become clearer in a moment. Consider the query string 'hello and the parent string 'ohelol':

current block: 'ohelo'
h e l o
query string table 1 1 2 1
current block table 1 1 1 2
matches? 1 1 0 02

Here, the distance is 2 because only two character counts -- that of the 'h' and that of the 'e' -- match. Now look at the next block:

current block: 'helol'
h e l o
query string table 1 1 2 1
current block table 1 1 2 1
matches? 1 1 1 14

Now our distance is 4. If this distance is equal to the number of unique characters in the query string (which it is), our current block is a matching block. Now do what you were doing before, but instead of comparing the lookup tables, just modify the distance based on whether the 'new' block character and 'old' block character are contained inside the query string, then compare the distance to the number of unique characters in the query string. Phew!

Complexity: O(n). Can I go to sleep now?

This is an incredibly irritating algorithm to code. Here's a C++ implementation (fitting to the 'horrendous style' paradigm of this blog) to keep you busy.

for (int i=0 ; i<query_length ; i++) {
	scanf("%c", &temp);
	query_table[letter_index(temp)] ++;
}
scanf("\n%s", parent);

for (int i=0 ; i<=letter_index('z') ; i++) expected_distance += (query_table[i] > 0);

for (int i=0 ; i<parent_length ; i++) {
	if (i - query_length >= 0) {
		distance -= block_table[letter_index(parent[i-query_length])] == query_table[letter_index(parent[i-query_length])];
		distance += --block_table[letter_index(parent[i-query_length])] == query_table[letter_index(parent[i-query_length])];
	}
	distance -= block_table[letter_index(parent[i])] == query_table[letter_index(parent[i])];
	distance += ++block_table[letter_index(parent[i])] == query_table[letter_index(parent[i])];
	output += (distance == expected_distance);
}



Friday, June 17, 2011

A Peculiar Numbering System

Using the partial function application and composition from the previous post, we can represent the natural numbers in an interesting way.

def partial(f, *p):
    return lambda *q: f(*(p + q))

def compose(*fs):
    return lambda x: reduce(lambda a,b: b(a), (fs+(x,))[::-1])

def represent(n):
    return compose(*(partial(int.__add__, 1),) * n)(0)

assert represent(0) == 0
assert represent(1) == 1
assert represent(1337) == 1337

The nth natural number is the function λx.x+1 composed with itself n times, applied onto the number 0.

This is related to the concept of Church numerals.

Friday, March 4, 2011

Partial Function Application and Function Composition

This following function demonstrates one of my favourite concepts in functional programming -- that of partial function application:

def pfa(f, *p):
    return lambda *q: f(*(p + q))

Partial function application is where you fix particular parameters of a function, producing a new function with the fixed parameters substituted in. For example, say I had a function add(a,b) = a + b. I can use partial function to 'fix' the first argument, a, to a particular number x, and this produces a new function add(b) = x + b.

In practice, PFA is often a more compact and readable version of the ugly anonymous functions you give to higher-order functions such as lambda. For example, consider the following code extract which transforms [a,b,c,...] to [2^a, 2^b, 2^c, ...]:

map(lambda n: pow(2, n), list)

Now take a look at how it is done with PFA:

map(PFA(pow, 2), list)

PFA in this instance creates a new function by substituting the value '2' as the first argument of the pow function, transforming pow(b, e) = b^e into pow(e) = 2^e.

Haskell's PFA notation is even more elegant; because everything in Haskell is just chained evaluation of partial functions, there's no special or unique notation for it, so the following is valid code:

map (2^) list

That maps the partial function 2^ onto list. Haskell also has lambdas, but they're no where near as short or understandable.

map (\x -> 2^x) list

The Python implementation above can also fix multiple arguments at one time, as in pfa(f(a,b,c), a, b) -> f(c). Sometimes it may be preferable to fix the second or third argument, while leaving the first argument variable -- Haskell deals with this elegantly as usual, but I haven't found a nice way to do it in Python.

note: apparently the functools module already has a partial() function for PFA. whoops :)

Function composition is another concept (again, originally mathematical) with similar roots. Here's an implementation:

def compose(*fs):
    return lambda x: reduce(lambda a,b: b(a), (fs+x)[::-1])


Essentially, this takes two functions f(x) and g(y) and produces a new function h(y) == f(g(y)). The new function is often notated as (f.g)(y), as in "f.g is the composition of f and g." Like partial function application, it is rarely a necessity to use function composition, but again it can make things more compact and more intuitive.

A common construct I use when dealing with strings is as follows:

new_s = old_s.strip().split()

This can be composed as follows:

new_s = compose(str.split, str.strip)(old_s)

A better example:

hash(bin(bool(id(sum(range(10))))))
# becomes
compose(hash, bin, bool, id, sum, range)(10)

Sunday, October 10, 2010

Balanced Brackets

Let us define a simple recursive language:

balanced = '()' | '(', balanced, ')' | balanced, balanced;

The theorems (valid clauses) of our language are made up only of the parentheses. We will call a string 'balanced' if and only if it can be expressed in our language. A friendlier definition of a balanced string is as follows:

  • () is a balanced string
  • Putting a balanced string inside a pair of opening and closing parentheses forms a balanced string
  • Concatenating two balanced strings forms a balanced string

Here are some examples:

Balanced
()
(())
()()
(()())
((()(()))()(()()))
((())(()))((())(()))
Not balanced
(
)
((
))
)(
(()
())
(()()(())

The task is to write a program which determines whether a given string is balanced. Additionally, the algorithm must reflect the language's definition as closely as possible. To reflect this goal I'm going to propose an arbitrary restriction: the algorithm should not involve any change in state. The only exception is the string being tested, which may change. Essentially what this means is that you can't use any variables which change their value. The algorithm will also almost certainly be recursive -- which makes sense because the language is defined recursively as well.



This took me a frustratingly long time. I rewrote the code several times over, each time slightly changing my algorithm to make it cleaner and more logical.

My function takes two parameters: a 'string' to test' (which is actually a list of that string, string) and a character which, when found at the same recursive depth, will signify that the string is balanced (expected).

The first thing in the function is a loop for handling cases where there are two or more balanced strings sitting at the same recursive depth.

I get the first element of the string (while at the same time removing it -- str.pop has that side effect) and store it in the variable first. I then do three tests:

  • If first is an opening parenthesis, I check if the string following it is balanced by recursively calling the function on the rest of string. I also tell it to stop when it reaches a ) character on the same recursive depth. If it tells me that it's not balanced, I stop and return False (if the substring isn't balanced then the whole string isn't balanced).
  • If first == expecting, this string (possibly a substring) is balanced because we've reached the closing character on the same recursive level as the opening one and so everything in between must have been balanced.
  • If anything else has happened, the string isn't balanced and we return False.

The loop can only quit normally (without getting out via a return statement) in two ways:

  • The string ended prematurely. In this case, we will still be expecting something, so the expression 'not expecting' will return False. Since the string ended prematurely the string isn't balanced, so we return False correctly.
  • The string ended at the top recursive level as it should. We aren't expecting anything, so it will return True.


def balanced(string, expecting=False):
    while string:
        first = string.pop(0)
        if first == '(':
            if not balanced(string, ')'):
                return False
        elif first == expecting:
            return True
        else:
            return False
    return not expecting

And my testing:

assert balanced(list('((()))'))
assert balanced(list('()'))
assert balanced(list('()()'))
assert balanced(list('(())((()))'))
assert balanced(list('()(()()(())()((())()))'))
assert balanced(list('(((((((()))())())())()))'))
assert balanced(list('((()(())())(())()(()()))'))
assert not balanced(list(')('))
assert not balanced(list('(('))
assert not balanced(list('()))'))
assert not balanced(list('(((()))'))
assert not balanced(list(')))'))
assert not balanced(list('(()'))
assert not balanced(list('(()()()'))

Thursday, October 7, 2010

Drawing the Koch Curve in Turtle


I love this fractal because it's so simple! I barely need to obfuscate the code and its definition still fits in 3 fairly neat lines. It's simple enough that I won't even bother giving an explanation so as not to insult your intelligence. If you're having trouble understanding the t.forward line, I explain that construction in my explanation for the Dragon Curve Fractal in Turtle (my first post).

import turtle
t = turtle.Pen()

def koch(dist):
    for angle in (60, -120, 60, 0):
        (t.forward if dist <= 10 else koch)(dist / 3.0)  # 10 is the break case. 3.0 is the line:subline ratio
        t.left(angle)

koch(400)

Thursday, September 30, 2010

Simulating a Conditional Expression via a List Indexed by a Boolean

It is often the case that you need to shave off just a few characters off a piece of code (there are many competitions where code length is restricted, such as the JS1K competition). The following describes a good way to condense conditional expressions. NOTE: This method only works if your condition returns True or False, or 1 and 0. It will not work if your condition returns, say, a non-zero value for True.

This is what it looks like:
first if condition else second
is equivalent to
[second, first][condition]
This works because condition will return True or False; True in Python is equivalent to 1 and False is equivalent to 0. We then use this 1 or 0 as the index to a list containing first and second; if the condition returns 0, we get the 0th item of the list: second. If it returns 1, we get the 1st item: first.

Pretty neat, huh?

EDIT 05/08/11: True conditional expressions only evaluate the output term. To get this effect, wrap each term in a lambda and call the entire expression's result.
[lambda: second, lambda: first][condition]()

Using reduce() to Implement map() and filter()

map() and filter() are by far the most popular of the big three higher-order functions that Python borrowed from functional programming. reduce(), the most convoluted and intricate of the three, has sadly been removed as a built-in in Python 3.0, forcing us to write neater code. Fortunately map() and filter() were both kept, even though list comprehensions are far easier to understand. There is hope for us yet!

map() is a very rigid function -- its behaviour is very easy to predict because it does exactly what it looks like it does. filter(), too, has very limited use and has a single specific purpose (filtering). map()'s usage has no overlap with filter()'s and vice versa.

reduce() is the odd one out -- it's incredibly flexible and fun to use. In fact, reduce() can be used as a replacement for both map() and filter().

def map(f, list):
    return reduce(lambda p, n: p + [f(n)], 
                  [[]] + list)

def filter(f, list):
    return reduce(lambda p, n: p + ([n] if f(n) else []), 
                  [[]] + list)

Off the top of my head I can't think of a way to accomplish either of these tasks without sticking a [] at the start of list. If you have a way, please comment!

Wednesday, September 29, 2010

Selecting an Optimal Next Move in Nim

Nim is a mathematical game where you have a number of piles and each pile contains a number of items. Players take turns at each taking at least one and at most all of the items from a single pile of their choosing. The winning player in nim (in its standard form) is the player who has left all piles empty following their move. Nim is traditionally played with 3 piles of 3, 4, and 5 items.

The way to guarantee a win is to always pick a move such that after the move has been executed, the binary digital sum (nim-sum) of each pile is equal to 0. The binary digital sum is the sum of the binary forms of the numbers, without carrying. This turns out to be the same as applying the exclusive-or bitwise operation over all piles.

Given the number of items in each pile, determine the optimal move you should execute in the form (pile to remove from, number of items to remove).
To calculate the nim-sum of the piles we have to apply XOR over all the piles. The easiest way to do this is with the reduce() function. If you haven't encountered reduce() before, it comes from the same functional roots as map() and filter(), but its purpose is to condense a list of type t into a single instance of t** by applying a function between each successive pair and 'folding' itself up. Here's an example of the sum() function implemented as a reduction:

>>> print reduce(lambda a, b: a + b, [1, 2, 4, 8])
15

This is evaluated as (((1 + 2) + 4) + 8). To get the nim sum, we apply a very similar procedure except instead of using a function to add two numbers together, we'll use a function to get the exclusive-or of two numbers. Fortunately Python's already given one to us: int.__xor__. The nim-sum of a game state is therefore:

nim_sum = reduce(int.__xor__, state)

But this is only half the job! We need to go through all possible moves we can make and find one that would leave the game state with a nim-sum of 0.

We're not too concerned about efficiency so we're gonna loop through every possible pile and for each pile, try to remove every possible number of items and see if any of the remaining game states have a nim-sum of 0. The game state will be represented as a list where each element represents the number of items in that pile. The following line of code will give us the result of removing num items from the pile pile:

game[:pile] + [game[pile]-num] + game[pile+1:]

All we're doing is modifying the selected pile then sandwiching the result in between the rest of the piles where it was before.

When you combine all this and add two loops you have yourself a function which, given a game state, will give you a move you can make which results in a nim-sum of 0.

def nextMove(g):
    for p in xrange(len(g)):
        for n in xrange(1, g[p]+1):
            if reduce(int.__xor__, g[:p] + [g[p]-n] + g[p+1:]) == 0:
                return p, n
    return None, None


**This is not always the case -- a reduce can theoretically return any type. However, if it is given a list of [t] it will generally return an instance of t.

Tuesday, September 28, 2010

Finding All Clauses of a Context-Free Grammar

This question was given in week 3 of the 2009 National Computer Science School (NCSS) Challenge. It remains one of my favourite competition problems.


The task is simple. Given a description of a context-free grammar (CFG) in Backus-Naur Form, generate all productions (clauses) of the CFG. Assume the grammar is not recursive (or else there are an infinite amount of clauses).


The solution revolves around the following two observations:
  • If <a> ::= <b> <c> the productions of <a> are the Cartesian product of <b> and <c>
  • If <a> ::= <b> | <c> the productions of <a> are the concatenation of <b> and <c>
For the record, | has the lowest precedence. That means that if <a> | <b> <c> | <d>, the solution is the concatenation of <a>, (<b> x <c>), <d>. 


The solution naturally lends itself to two parts. We need to parse the grammar into a 
format we can deal with, then we need to generate the productions.



Parsing



BNF is in a key:value format, so it makes sense to use a dictionary. We can use a nested list format for each definition. A list will contain sublists which were separated by | . It's sorta hard to explain so here's an example:


<word> ::= <a> <b> | <c> <d> | <e> "lol i troll you"


Will become

mapping['<word>'] = [['<a>', '<b>'], ['<c>', '<d>'], ['<e>', '"lol i troll you"']]

Note the fact that an entire string in BNF -- including the surrounding quotation marks. We're gonna use some very simple regular expressions to help us identify <these> and "these".

import re

def parse(grammar):
    mapping = {}

    for definition in grammar.strip().splitlines():
        word, meaning = definition.split('::=') # separate key and value                                                                     
        superList = meaning.split('|')
        newSuperList = []

        for ss in superList:
            newSuperList.append(re.findall(r'(<.*?>|".*?")', ss))

        mapping[word.strip()] = newSuperList

    return mapping


Of course, I'm interested in exciting code, even if I have to be a bit ugly about it.

import re
def parse(g):
    return dict([(w.strip(), [re.findall(r'(<.*?>|".*?")', s) for s in m.split('|')]) for w, m in [d.split('::=') for d in g.strip().splitlines()]])

That's better. Where's the fun in readable code?



Generating



Now this is the fun part. We're gonna solve this problem recursively, simply because this problem is a perfect example of one which at first appears complex but can be defined by a set of simple recursive rules. The two recursive rules are in the introduction above. We must also add our break case: we stop recursing whenever we reach a literal (a string in the BNF).

Our recursive function will take in a single parameter: the current term that we're generating clauses from. It will return a list of strings (technically it will be an iterator of strings), where each string is a possible clause in the CFG. If the term parameter is a literal (it matches the pattern ".*") then we've hit the break case and we return a list containing just that term. Otherwise, we expand the term and recurse deeper by applying our recursive rules.


This is a cool challenge, so I'm going to leave you with two options: either figure out how to program that or try to understand my one-liner here. Have fun!

import itertools
def generate(term):
    return re.findall(r'"(.*?)"', term) if re.match('".*', term) else itertools.chain(*[map(''.join, itertools.product(*map(generate, p))) for p in syntax[term]])

Finally, here is the whole code with a CFG provided as an example. It's a very very simplified English sentence CFG. The code works perfectly and because we use itertools.chain() the returned object will be an iterator.


from itertools import chain, product
from re import match, findall

GRAMMAR = '''                                                                                                                                                                         
<sentence> ::= <noun phrase=""> <verb phrase="">                                                                                                                                            
<noun> ::= "boy " | "troll " | "moon " | "telescope "                                                                                                                                 
<transitive verb=""> ::= "hits " | "sees "                                                                                                                                               
<intransitive verb=""> ::= "runs " | "sleeps "                                                                                                                                           
<adjective> ::= "big " | "red "                                                                                                                                                       
<adverb> ::= "quickly " | "quietly "                                                                                                                                                  
<determiner> ::= "a " | "that " | "each " | "every "                                                                                                                                  
<pronoun> ::= "he " | "she " | "it "                                                                                                                                                  
<noun phrase=""> ::= <determiner> <noun> | <determiner> <adjective> <noun>                                                                                                               
<verb phrase=""> ::= <intransitive verb=""> | <transitive verb=""> <noun phrase="">                                                                                                               
'''

def parse(g):
    return dict([(w.strip(), [findall(r'(<.+?>|".+?")', s) for s in m.split('|')]) for w, m in [d.split('::=') for d in g.strip().splitlines()]])

def generate(term):
    return findall(r'"(.*?)"', term) if match('".*', term) else chain(*[map(''.join, product(*map(generate, p))) for p in syntax[term]])

syntax = parse(GRAMMAR)
print list(generate('<sentence>'))

EDIT: Alright, that was a bit mean of me. Here's an expanded version of the generate() function. Whenever I'm obfuscating I always write it up in an expanded version, get it to work, then compress it. In my idiocy, I deleted my expanded 'working' version, so I'm writing this one up on the spot. Note that this one returns a list instead of an iterator.

def generate(term):
    if match(r'".*"', term):
        # if the term is an ebnf string, return it as a single-item list
        #  but strip off the quotation marks
        return findall(r'"(.*?)"', term)
    else:
        # otherwise, solve this subproblem and return a list containing all 
        solutions = []  # this will contain the function's output
        for toConcatenate in syntax[term]:
            expansions = [generate(toExpand) for toExpand in toConcatenate]
            for subsolutions in product(*expansions):
                solutions.append(''.join(subsolutions))
        return solutions