Sunday, December 23, 2012

Solving the Countdown problem in Haskell


One of the challenges in the British game show Countdown was to find a way to create a number using other numbers and basic mathematical operators. For example, given the set {2,3,4,5} generate 14, using +, -, *, /. A solution is 5 *4 - 2 * 3. Note that each  number should be used exactly once.

An interesting problem is writing an application that lets the computer do all the work for you, so that it finds a solution for a given target number, or maybe even all the solutions!

Problem analysis

To find all the ways to combine numbers to yield a certain result we have to generate all the possible ways of  combining the numbers (subproblem 1) and then filter out the ones that have our desired result (subproblem 2). But how do we generate all possibilites? What we really want is to generate all the possible valid mathematical expressions with our numbers. We recognize that these expressions can be written as a tree:

For example, $$5 \cdot 6 - 2 \cdot 3$$:

Tree representation of mathematical expression
Since the operators we have chosen all have two elements, our expressions are binary trees. Thus, this means that our algorithm has to generate all possible trees.

For this task I have chosen Haskell, since it is by far the easiest language to work with trees, as we will soon see. If you are interested in another example why Haskell - or other functional languages - are very convenient to work with if you have trees, check out my fully fledged  spl compiler. This project is a compiler for a custom programming language, written in Haskell in only 900 lines! I have written my own parser, abstract syntax tree and Intel backend using the same techniques you'll see below. If you are interested in this project, let me know, and I may write a tutorial on this.


First, we have to define the grammar of our mathematical expressions:
data Exp = Num Int | Plus Exp Exp | Sub Exp Exp | Mul Exp Exp | Div Exp Exp deriving Eq

Our expression is defined as being a number, or any of our binary operators. The elements (branches) of each operator are mathematical expressions themselves as well. You can immediately see that this recursive definition of a tree structure is much, much smaller than it would have been in for example Java or C.

In an object oriented language the solution for such a structure would be to make a base class Exp, that has 5 subclasses. Each subclass then has two variables of type Exp, except for the Num subclass which has an integer variable.

Before moving to the main algorithm, I'll show some code for convenience functions. The first implements a Show function for our custom data type, making a tree easier to read:
instance Show Exp where
 show (Num n) = show n
 show (Plus a b) = "(" ++ show a ++ "+" ++ show b ++ ")"
 show (Sub a b) = "(" ++ show a ++ "-" ++ show b ++ ")"
 show (Mul a b) = "(" ++ show a ++ "*" ++ show b ++ ")"
 show (Div a b) = "(" ++ show a ++ "/" ++ show b ++ ")"

If you now type
show (Plus (Num 5) (Mul (Num 6) (Num 2)))

in the Haskell interpreter, the result would be (5 + (6 * 2)). In the code above you can see pattern matching. If you are not familiar with functional languages, you can read about this feature here.

We also need a function that actually calculates the result of our abstract expression. Using pattern matching again:
reduce :: Exp -> Int
reduce (Num n) = n
reduce (Plus a b) = reduce a + reduce b
reduce (Sub a b) = reduce a - reduce b
reduce (Mul a b) = reduce a * reduce b
reduce (Div a b) = reduce a `div` reduce b

Now we can think about the algorithm itself. We want to use all the numbers and use them only once, so we have to generate all the possible ways to select an element from a list and also keep track of which items we have not used yet. That's what the following function is for:
-- list should have at least 2 elements
split :: [a] -> [a] -> [(a, [a])]
split [] _ = []
split (x:xs) acc = [(x, acc ++ xs)] ++ (split xs (acc ++ [x]))
Split turns a list [1,2,3] into [(1, [2,3]), (2,[1,3]), (3, [1,2]). Thus, when we select a binary operator, we use the first number of the tuple as its left branch and we can use the list of remaining number to generate a subtree for the right branch:
pos :: [Int]-> [Exp]
pos [n] = [Num n] -- one number left? no binary operators possible
pos n = concatMap allpos subtree
 subtree = map (\(a,r) -> (Num a, pos r)) (split n [])
 allpos (a,b) = concatMap (\e -> pos' a e) b
 pos' a b = [Plus a b, Sub a b, Sub b a, Mul a b, Div a b, Div b a]
Let's look at the subtree function first. This generates all possible subtrees where the left branch is a given number. The function pos' gives all the possible binary operators where a and b are branches. As you can see, for the subtraction and division operator, the left and right branches interchanged is also a valid option. That's because 1 - 2 != 2 - 1.

Finally, the allpos function generates puts all these possiblities in a list. Of course, some of these expressions may be invalid due to the rules of the game. For example, it's not allowed to have a negative value or a fractional values in any of the subresults. Another thing that should be filtered is duplicates: 1 + 2 is the same expression as 2 + 1, so it should appear only once. That's why we have the following filter function:

-- remove illegal attempts and attempts that are the same
valid :: Exp -> Bool
valid (Num n) = True
valid (Plus a b) = valid a && valid b && reduce a > reduce b -- only allow one
valid (Sub a b) = valid a && valid b && reduce a - reduce b >= 0 -- subsolutions should be positive
valid (Mul a b) = valid a && valid b && reduce a > reduce b
valid (Div a b) = valid a && valid b && reduce b > 0 && reduce a `mod` reduce b == 0
And finally we have a function that receives a list of numbers and the target number as input and that returns all the possible expressions:
findVal :: [Int] -> Int -> [Exp]
findVal n res = filter (\e -> valid e && reduce e == res) (pos n)

And that's the solution to the problem! If you run this application, you will see that contrary to what you may expect, very little memory is used. This is because Haskell does lazy evaluation, so the lists are actually not computed beforehand but each element is evaluated one at a time.

There are some optimzations, which I leave as an exercise for the reader: for example, you can reject wrong expressions much earlier. If you generate a binary operator with two numbers, say Plus (Num 2) (Num 4), you can reject it immediately (because we only allow 4 + 2 and not 2 + 4).

And that's it for today. I hope you've learned something new and interesting :)

Wednesday, July 18, 2012

Making gnuplot plots look like Mathematica plots

For an article I am writing I wanted to plot the following function:  $$ \theta \left(x > \frac{1}{2} \right) \left(\sin(4 \pi x) + \sin(4 \pi y) +2 \right) $$
Wolfram Alpha gives the following pretty plot:
Plot from Wolfram Alpha
However, since you need to pay now to export this picture in a more usable format (like svg, eps or pdf) I decided to recreate the picture in gnuplot. Friendly advice for you others writing articles: always use vector formats where possible, otherwise your graphics will look horrible on higher resolutions!

To plot the function in gnuplot we first have to do a little trick, because the step function is not included in gnuplot. It does have the sign function, however. So we can do the following: $$ \theta(x > \frac{1}{2}) = \frac{1 + \text{sgn}(x-\frac{1}{2})}{2} $$
A first attempt, with default settings is:

set terminal postscript eps
set pm3d at s hidden3d
unset hidden3d
unset surface
set xrange [0:1]
set yrange [0:1]
set samples 20,20

splot (1+sgn(x-0.5))*(sin(4*pi*x) + sin(4*pi*y)+2)/2 title ''

This gives

This doesn't really look too good. First thing that has to be changed is the color palette and the color of the grid. The color palette can be changed like this:

set terminal postscript eps
set palette defined ( 0 "#bd2c29", 2 "#ffd35a", 6 "white")

The numbers in front of the color codes are the values of the function that the color should correspond to. Thankfully, gnuplot automatically interpolates the colors.

Next up is the grid color. The color can be set to a given line style by adding the style number at the end of the pm3d options:
set pm3d implicit at s hidden3d 100
set style line 100 lc rgb '#000000' lt 1 lw 0.6

Here I have defined the grid style to be a thin black line. I have chosen index 100 to make sure it doesn't clash with other settings.

After removing the color bar and shifting the image down with ticslevel, we get the final image:
This looks a lot better! The only thing missing is the shadowing of surfaces angled away from the camera. I am not aware of a way to produce these in gnuplot. If anyone has additional advice to improve the plot, please let me know!

The final version of the code:

set terminal postscript eps
set output '| epstopdf --filter --outfile=plot.pdf'
set palette defined ( 0 "#bd2c29", 2 "#ffd35a", 6 "white")
set pm3d implicit at s hidden3d 100
set style line 100 lc rgb '#000000' lt 1 lw 0.6
unset hidden3d
unset surface
set border 4095 front linetype -1 linewidth 0.8
set ticslevel 0
set xrange [0:1]
set yrange [0:1]
set ztics 1
set isosamples 40,40
set samples 20,20
unset colorbox

splot (1+sgn(x-0.5))*(sin(4*pi*x) + sin(4*pi*y)+2)/2 title ''

Wednesday, July 4, 2012

What is the Higgs boson?

Today the Higgs boson has been found at CERN, a great accomplishment in the field of theoretical particle physics. This discovery means that the model we are currently using is completed: the final ingredient has been found. This model, called the Standard Model, yields very good predictions of what is happening in nature.

But what is this Higgs boson, what does it do and why didn't we see it before?

The mystery of the massive particles

When physicists were deriving the physics of elementary particles, they used certain assumptions about the universe, most of them are related to symmetries. What's a symmetry? Take for example rotational symmetry. Imagine you are standing in the middle of a round football stadium. It does not matter if you rotate, you will still see the same image. There is a very fancy theorem that says that every symmetry has an accompanied conserved quantity. The rotational symmetry, for example, yields conservation of angular momentum, whereas translational symmetry in time (meaning physics looks the same now and 10 minutes from now) yield conservation of energy.

Using some more advaced symmetry groups all of our beloved standard particles can be found:

The standard model (without Higgs)

The group in red are called gauge bosons. They are resultants of a very special symmetry and they can be interpreted as forces working on all the other particles. In fact, every elementary force has a matching gague boson. The photon (symbolised with a gamma), or "light", transfers electromagnetic forces, the gluon the strong force and Z and W the weak force.

But there is a problem. The theory requires the W and Z boson to be massless in order for the theory to work. This is obviously not the case, because these particles have been found to have mass. In fact, they are pretty heavy!

The Higgs mechanism

When physicist discovered that the theory was in trouble, they started searching for solutions. The symmetry of nature that can be used to find the gauge bosons worked far too good to just abandon it. Then, a new idea rose: we keep the symmetry, but we say nature broke it spontaneously. Wait, what? 

How does this work? Imagine you have a piece of plastic between your fingers and you apply a little pressure to it. Initially, the plastic stick is straight, even when you apply a little pressure to it. Now the system (nature) has a symmetry: there is for example rotational symmetry if you are standing at the position of the rod (if you ignore the hand). 

But remember the hand was applying a force? This means that the stick could spontaneously jump to either the left or the right, as can be seen in the picture (try this at home ;) ). Now the symmetry that nature had is broken, spontaneously and without external help. Note that in the analogy, the hand should be seen as an initial energy that the universe has.

Spontaneous symmetry breaking

A similar scenario this of a toy car that is sitting on a small mountain, but it has some energy which makes the car shift a little. Given enough time, the car will roll off the mountain into a valley, where it will remain.

Imagine a toy car in this figure

Now, it turns out that breaking this symmetry can be done by introducing a particle: the Higgs particle. When analysing the properties of this particle, it turns out that this particle interacts with all the other particles, including itself! Due to these interactions, the particles gain mass.

Why did this take so long to find?

Well, we didn't know where to look. Basically, in the LHC particles are shot at eachother with incredible speeds, and thus incredibly high energies. Since E=mc^2, it takes very high energies to create particles with high mass. The problem is, we don't know the mass of the Higgs, it is a free parameter of the theory. This means the mass is not predicted by the theory. And because of that, the whole LHC range has to be explored (with no guarantee that the Higgs boson is actually in this range).

Complications arise due to the fact the Higgs particle decays intoa lot of other particles very, very quickly. So we can actually only see the remnants of the Higgs particle. When there is enough energy to create a Higgs, we'd expect to see more of these remnants: a peak in a graph.

Today, the people at CERN announced that the peaks they found in many expirements, are trustworthy enough to say that the Higgs particle has been found at about 126 GeV. This is great news, because now we can confidently continue the search for theories that expand the standard model, like supersymmetry and other theories that may introduce gravity into the picture.

We aren't quite there yet, but a theory of everything is coming!