I'll have a Buchburger with fries

This is, as usual, a literate Haskell program. So let me start with

> import Data.Map as M
> import Data.List as L

Everyone else is at it, so why not me? Here's my approach to solving this problem. Let's warm up with a simpler menu. Consider:

Fries $3
Burger $5

What can I get that costs exactly $13?

My approach is going to revolve around rewrite rules. I'm going to find a set of rewrite rules that converts money into fries and burgers, and does so in such a way that it always eliminates the money, if it is possible to do so.

Let's use D to represent a dollar, F to represent fries and B to represent a burger. Here's our first attempt at some rewrite rules:

DDD->F
DDDDD->B

Now we can take the string DDDDDDDDDDDDD and try to rewrite it. We could start with the first rule:
DDDDDDDDDDDDD
->DDDDDDDDDDF
->DDDDDDDFF
->DDDDFFF
->DFFFF
But now we're stuck with a D and no way to get rid of it.

Before we go much further, however, it's time to simplify our notation. Firstly, unlike normal rewrite rules, we don't care about the order of our symbols. So let's assume that we can reorder symbols any time we want in our strings and rules. So FDD=DFD=DDF and so on. Let's also use Xⁿ to mean n copies of X in a row. Now we can write our strings as mononomials in D, F and B. That makes things much easier. We can now write our rules as

D³->F
D⁵->B

We've seen that just picking the first rule each time can get us stuck. It seems that the solution should be to come up with an algorithm to generate the correct order. But there's another strategy - it's to modify the rules so that the order doesn't matter. A set of rules with this property is said to be confluent. Amazingly you can turn any reasonable set of rewriting rules into a confluent one merely by adding in elements.

So let's define a monomial type:

> newtype Monomial a i = Monomial (Map a i) deriving (Eq,Show)
> unM (Monomial x) = x

> divides (Monomial a) (Monomial b) = isSubmapOfBy (<=) a b
> lcm (Monomial a) (Monomial b) = Monomial $ unionWith max a b

> instance (Show a, Show i, Num i,Ord a) => Num (Monomial a i) where
>     Monomial a*Monomial b = Monomial $ unionWith (+) a b

> instance (Show a, Show i, Num i,Ord a) => Fractional (Monomial a i) where
>     Monomial a/Monomial b = Monomial $ M.filter (/=0) $ unionWith (+) a (M.map negate b)

Rewrite rules are simply pairs of monomials:

> data Rule a i = Monomial a i :-> Monomial a i deriving (Eq,Show,Ord)
> lm (a :-> _) = a

Now we can write code to apply these rules to a string. The first function simply applies the rule to a term:

> apply (p :-> q) r = r/p*q

applicable searches for the first applicable rule in a list of rules:

> applicable r rules = find (\(p :-> _) -> divides p r) rules

And reduceStep applies it. Note that it returns the result twice in a pair. This is simply so that it works with the next function. In the event that no rule is applicable it returns Nothing.

> reduceStep rules r = applicable r rules >>= return . dup . flip apply r
>     where dup x = (x,x)

reduce repeatedly applies reduceStep until it can be applied no more. Repeatedly doing something until Nothing is returned is best handled by unfoldr:

> reduce rules r = last $ (r:) $ unfoldr (reduceStep rules) r

Now consider our example rewrite rules. Let's find a simple example of where confluence fails to hold. Consider two different ways to rewrite D⁵. We could use the first rule and get D⁵->D²F or use the second and get D⁵->B. Either way we terminate - and get different results. How can we get out of this impasse? There's a really obvious solution, add in a new rule D²F->B. Now if we start with D⁵ and apply the first rule then we can get out of it by applying this new third rule. Either way, we get B.

Before we generate new rules, we need to say something about how we can ensure termination. If we have one rule that says X->Y and another that says Y->X then we're going to get stuck. So what we can do is insist on an ordering on our monomials. We'll restrict ourselves to rules where the left hand side is larger than the right hand side. So here's an ordering known to algebraists as lexicographic ordering. The idea is that the symbols in the monomial have an ordering on them. One monomial, a, is considered bigger than another, b, if the power of the least symbol (in the symbol ordering) in a/b is negative. (BTW See postscript for a comment on this.) What this means is that the rewrite rules will try to eliminate symbols earlier in the ordering and replace them with symbols that come later. By making our currency symbol the first in the ordering the rewrite rules will tend to act so as to eliminate the currency. And that's what the original problem asked us to do. Note that there are other orderings we could use so for generality Monomial ought to be a typeclass.

> instance (Show a,Ord a,Ord i,Num i) => Ord (Monomial a i) where
>     compare x y = let
>         Monomial u = x/y
>         in
>         if M.null u
>             then EQ
>             else if snd (M.findMin u)>0
>                 then LT
>                 else GT

We also use a function to ensure that any rule is in canonical form with the greater term on the right:

> order (a :-> b) = if a>b then b :-> a else a :-> b

If we find a complete confluent set of rules that respect the ordering in this way, how do we know the rules terminate? Suppose the complete ordered list of symbols in our rules is (a₁,a₂,...,a_n) then any string a₁^k₁a₂^k₂...a_n^k_n can be written as a vector (k₁,k₂,...,k_n). Each application of a rewrite rule replaces the corresponding vector with one that is lexicographically earlier. It's pretty clear that such a process must terminate. (Aside: You can think of this as using the fact that ωⁿ is well-ordered. More general kinds of rewrite rules make use of more complex ordinals to prove termination. This goes to show that ordinal arithmetic is actually a useful tool.)

So back to what we were doing: finding how we can add new rules to ensure confluence. Suppose we have two rules whose LHS 'overlap'. If we think in terms of monomials then we mean two monomials that have a non-trivial common factor. If we form the lowest common multiple of these monomials then we now have a 'minimal' term to which both rules can be applied. In fact, that's how I got the example above. One rule has P³ on the LHS and one has P⁵. The lowest common multiple of these is P⁵ and I investigated how the two rules apply to this one term. When we did this, we ended up with two different terms which we simply made into a new rule. Here it is in Haskell:

> newRule pq@(p :-> _) pq'@(p' :-> _) = let
>     j = Main.lcm p p'
>     in order $ apply pq j :-> apply pq' j

Given a rule, we also simplify the left and right hand sides with normal:

> normal rules (p :-> q) = let
>     p' = reduce rules p
>     q' = reduce rules q
>     in if p'/=q'
>         then Just $ order $ reduce rules p :-> reduce rules q
>         else Nothing

The problem is that when we generate a new rule there now appear a whole bunch of new pairs that could be considered because you need to consider pairs made up of the new rule and all of the old rules. Our strategy for dealing with this will be simple: take the list of all pairs of rules and generate new rules from those. Every time we generate a new rule we need to add in the new rule pairs that are generated. We maintain a list of pairs found so far, chewing pairs for newRule from the front, and adding newly generated pairs to the back.

The first line is easy:

> complete incomplete = let

We make all possible pairs. We don't want to consider pairs twice so we use the inherited ordering on pairs:

>     pairs = [(x,y) | x <- incomplete, y <- incomplete, y>x]

Here's a little heuristic. Things go faster if we use monomials that come lower in the ordering before higher ranking monomials. So the completed set of rules is going to be sorted:

>     rules = sort incomplete

Now here's a weird bit to explain. This isn't essential to the mathematics but is is essential to this piece of code. The local function iterate is going to return a pair that consists of the remaining list of rule pairs and the completed rules. But iterate also consumes rule pairs. So this function is actually going to consume its own output. There's an important issue with this: eventually the rate of new rule creation will peter out so that it will consume its way to the end. At this point iterate will try to consume the very element it's expected to generate and spin in an infinite loop. (The dangers of unguarded recursion!) So we also maintain another list whose length is the same as the distance between the front and back of our pair list. When this first list is empty, we know we can return [] for the second list. It doesn't matter what this list is, as long as it has the correct length. This 'phantom' list is the first argument to iterate. But for now I suggest ignoring it and coming back to it if you ever get to a second reading.

Here's the base case I mentioned:

>     iterate [] pairs rules = ([],rules)

And here is the guts of the loop:

>     iterate (_:n) (p:pairs) rules = 

p is the first pair chewed off the front. We generate a new rule from it and normalise it:

>         case normal rules $ uncurry newRule p of

If we get a new rule we now need to add it to our list of rules, and add in all the new pairs it generates:

>             Just pq -> let u = [(x,pq) | x <- rules, x /= pq]

Note the use of insert to maintain the order.

>                            (pairs', rules') = iterate (u ++ n) pairs $ L.insert pq rules
>                            in (u ++ pairs', rules')

Otherwise we don't need to generate anything:

>             Nothing -> iterate n pairs rules

Here's where we tie the knot by feeding our output back in as input:

>     (pairs', rules') = iterate pairs (pairs ++ pairs') rules
>     in rules'

It's not obvious at this point that the algorithm terminates. It does, but you'll have to trust me.

Anyway, that's not all. We also do a bit of tidying up by eliminating redundant rules. In a confluent system, if there are two rules we can apply, it doesn't matter which one we choose. Now consider a pair of rules looking something like: AB²C -> ... and AB³C² -> ... . If we can apply the second rule, then we could also apply the first rule. But in a confluent system, it doesn't matter which rule we choose. So we can just throw out the second rule here. So we tidy up our rules by eliminating rules made redundant by another. Turns out that at this stage we have a canonical form. For any given set of rules, any completion with the reducant elements deleted is the same as any other (modulo reordering). So our completion process is itself confluent! This isn't hard to prove but I'm not going to do it now.

> nonredundant rules r = L.null [s | s <- rules, s/=r, lm s `divides` lm r]
> tidy rules = L.filter (nonredundant rules) rules

And now canonical puts our rules into a complete canonical form:

> canonical rules = tidy $ complete rules

Back to the original problem:

> data Menu = P 
>     | MixedFruit 
>     | FrenchFries 
>     | SideSalad 
>     | HotWings 
>     | MozzarellaSticks 
>     | SamplerPlate deriving (Eq,Ord,Show)

Note that deriving Ord. This is so that the lexicographic ordering above can work.

> makeRule a b = (Monomial $ M.fromList a) :-> (Monomial $ M.fromList b)
> rules = [
>         makeRule [(P,215)] [(MixedFruit,1)],
>         makeRule [(P,275)] [(FrenchFries,1)],
>         makeRule [(P,335)] [(SideSalad,1)],
>         makeRule [(P,355)] [(HotWings,1)],
>         makeRule [(P,420)] [(MozzarellaSticks,1)],
>         makeRule [(P,580)] [(SamplerPlate,1)]     
>     ]

> crules = canonical rules

> ex n = Monomial $ fromList [(P,n)]
> example0 = reduce crules (ex 1505)

> main = print (reduce crules (ex 1505))

Now, some caveats. I don't claim that this is a fast method, just a different approach. Completing a set of rules is an expensive operation and if you start with N rules, the worst case running time is something horrendous like O(exp(exp(N)). On the other hand, once you've found a set of canonical rules, reduction with them is very fast. This makes it ideal for problems involving repeated use of the same menu. Or for when you have a lot of spare cash:

> example = print (reduce crules (ex 10000115))

And now I have a confession to make. I've really just described the Buchberger algorithm for Gröbner basis computation. It's not quite the Buchberger algorithm because I'm only considering binomials, not polynomials - but you'll notice my code is remarkably similar to the code at the other end of that link. The rule p :-> q represents the binomial p-q. Unfortunately, the Buchberger algorithm is usually presented in its more general form - but understanding that requires a background in commutative algebra and algebraic geometry. I hope I've managed to present a form of it in a way that's digestible by people without abstract algebra experience. The binomial Buchberger algorithm has many applications, including a variety of applications in integer programming.

Note that this is also very close to the Knuth-Bendix algorithm. However, the latter algorithm is for the noncommutative case, and it isn't guaranteed to terminate.

BTW I recommend playing with the completed rules by hand - at least for simple examples. It's interesting to see what steps are taken to reach the solution.

PS My monomial ordering is the opposite of the usual lexicographic monomial ordering. This is because I want my rules sorted with what is usually considered the largest monomial at the front.

PPS Can I get rid of that weird phantom list thing? You can replace the list with an integer for its length, but it's not clear to me that this is an actual improvement.

Update: I forgot to mention this. Check out David Amos's implementation of the full Buchberger algorithm which you can use to solve many more problems than my simplistic monomial version.

Data and Codata

In Gödel, Escher, Bach, Hofstadter introduces the programming language Bloop. Bloop is a little like BASIC except that it forces programmers to specify in advance how many times each loop will iterate. As a result, Bloop programs have the handy property that they are guaranteed to terminate. Unfortunately this property also makes it impossible to write something like an operating system in Bloop. To write an OS we need open-ended loops that keep running until the user explicitly chooses to shut the OS down. One solution to this is to write code in Floop. Floop allows us to write unbounded loops, the equivalent of C's while (1) { ... }. The problem with that, however, is that we can write runaway infinite loops that never terminate and never give us any output. Is there some language that lies between Bloop and Floop that can give us unbounded looping when we need it, but which never allows us to hoist ourselves by our petards by writing runaway loops?

Wishing no disrepsect to OS writers, at first blush it might seem that the distinction between a runaway loop and an idle OS is too fine - if we can write an infinite loop that does something useful, then surely we can write a useless one too. But it turns out that there is a very elegant and well-principled way to distinguish between these kinds of loops, and this allows us to write open-ended interactive software in a programming language that nonetheless always produces a well-defined output, no matter what the input. In order to do this we need to distinguish between two kinds of data: data and codata. By ensuring that a function expecting codata never receives data, and vice versa, we can ensure that even programs with open-ended loops always produce a well defined output.

The concepts I want to talk about are very general and can apply to whatever programming language you use. I'm going to use some simple Haskell examples but most of these will translate to other languages. So consider something like

sum [] = 0
sum (a:as) = a + sum as

This sums the elements of a list. Note how it's well behaved as long as we give it a finite list as input. On the other hand, consider

sum' [] = 0
sum' a = sum' (1:a) - 1

This isn't well behaved at all. Except when you input empty lists, it never gives a result. From a mathematical perspective it's not a good definition either, there are many functions that satisfy these two properties. Is there some general principle at work here that allows us to tell immediately that one of these terminates and the other doesn't? We know from Turing that there is no procedure that guarantees we can always tell such programs apart, but in this case there is something that we can easily point to. In the first program, the right hand side of the second line uses the sum function recursively but we only apply it to a strict subpart of the input, its tail in fact. In the second example we apply sum' to something that contains the argument. The former function is using what is known as structural recursion, and it's not hard to convince yourself that structural recursion, applied to finite datastructures, always terminates.

So if we limit ourselves to structural recursion we can guarantee our programs will always terminate. But what about a definition like:

fact 0 = 1
fact n = n * fact (n-1)

That doesn't appear to use structural recursion. But we can view it as such like this. Define the natural numbers as follows:

data Nat = Zero | S Nat

0 is represented as Zero, 1 is represented as S Zero and so on. We can represent every natural this way. Here's the important thing: if n>0 then n-1 is simply a subpart of n. So we can view this kind of recursion as a type of structural recursion.

(In fact, by a curious quirk of the Haskell 98 standard we can rewrite our definition to look more like a structural recursion:

fact' 0 = 1
fact' (n+1) = (n+1) * fact' n

I'm guessing this feature is in Haskell precisely so that people can 'pretend' they are using structural recursion with the + in n+1 serving a role as a kind of pseudo-constructor.)

So we have a nice rule for ensuring our code terminates. But sum fails when applied to infinite lists. Should we simply rule out infinite datastructures? That seems a bit drastic. The example that convinced me to look into Haskell was

fib = 1 : 1 : zipWith (+) fib (tail fib)

We really don't want to rule out such a succinct definition of the Fibonacci numbers. But how can we allow such structures when we have functions like sum sitting around. Applying sum to fibs will obviously cause non-termination.

Let's consider another example:

sumSoFar x [] = [x]
sumSoFar x (y:ys) = x : sumSoFar (x+y) ys

Like sum, this fails to terminate for an infinite input. But unlike sum, it's possible to make sense of it. If the inputs were 0 and the infinite list [1,1,1,1,...] then the result would be [0,1,2,3,...]. The program might not terminate, but from a mathematical perspective this is a completely well defined function. What's more, suppose the input list represented a stream of data being input at a keyboard, and that the output was displayed on a screen one element at a time, then we'd have a simple cash register. This program might not terminate, but it's completely well behaved. Note that this could only work in a lazy language. A strict language would want to evaluate the entire list before outputting anything. But in a lazy language we can start outputting the beginning of the list before the rest of it is computed.

From the point of view of using infinite lists, it's sum that's badly behaved, and sumSoFar that's well behaved. Again, can we point to a simple distinction between these two programs that explains this? Turns out we can, and in some sense it's dual to what we said above. sumSoFar is well behaved because when we recursively call sumSoFar on the right hand side we do so from inside a list constructor. (Remember that : is the list constructor.) This is known as guarded recursion and it guarantees that even though our programs don't terminate, they still define a unique mathematical function and result in a well behaved program. In the case of sumSoFar, each time we look at another element of the result, we trigger another lazy evaluation that terminates. But the right hand side won't just run off and compute endlessly until we do that triggering because the recursion is 'guarded' within a constructor. (That, by the way, was the main point of this article, so you can probably relax now.)

Note the duality: in structural recursion we 'deconstruct' the argument and then we're allowed to recurse. In guarded recursion we recurse first, and then we're allowed to use the constructor.

So we've almost achieved our goal of describing rules to allow open-ended loops because we've managed to give a rule for writing functions that are guaranteed to be well-behaved even though they manipulate infinite data structures. But we're not quite home yet - we still can't have functions like sum coexist with infinite lists. How can we ensure that an infinite list is never handed to sum?

Consider the definition of Haskell lists. It's something like this:

data [a] = [] | a : [a]

I.e. a list of a's is either the empty list [] or it's made from an a and a list of a's.
You can think of this as an equation in [a]. In Haskell we take this as uniquely defining what [a] is, but in reality there is more than one solution to this equation. Consider the type consisting of only finite lists. That satisfies this equation. A finite list is either an empty list, or an element followed by a finite list. Similarly a possibly infinite list is either an empty list, or an element followed by a possibly infinite list. There is an ambiguity. Finite lists form, what is in some sense, the smallest possible solution to this equation. The possibly infinite lists form the largest possible solution. Haskell takes the largest possible solution.

But suppose we were to distinguish between these two different solutions. We could use the keyword data to mean the smallest solution and codata to mean the largest solution. The former represents data, and it's always finite. The latter represents what we call codata, and it's possibly infinite. And now we can refine our rules for well-behavedness. Consider data and codata to be distinct types. In a strongly typed language this immediately allows us to restrict sum to data, not codata. The rule is: you're only allowed to use structural recursion with data and guarded recursion with codata. With that rule, we're guaranteed that our recursions will always be safe, and yet that we can still have open-ended loops in our code. Sometimes these are called recursion and corecursion respectively.

And now we can go a little crazy with language. When we want to prove that a structurally recursive program terminates we use induction. This doesn't work straightforwardly for corecursion, so instead we use a principle called coinduction. Recursive programs typically terminate. Corecursive programs don't necessarily terminate, but they're still well-behaved as long as they keep on going whenever we give them input. We can call this cotermination. And so on... I'm not going to say what coinduction is because I'd have to talk about bisimulation, and this post would get way too long.

So now a mathematical aside. The reason for all the "co"s is that data and codata have categorical descriptions and they turn out to be dual to each other. But you don't hear mathematicians talking about corecursion and coinduction much. Why not? Well one of the axioms of set theory is the Axiom of Foundation. One way of interpreting this is the statement that there is no infinite sequence

...a₃∈a₂∈a₁∈a₀.

So even though we can construct infinite lists in mathematics, we can't construct 'infinitely deep' lists. This means that in mathematics we can use a form of structural recursion. In fact, the familiar principle of induction follows from this. So for many of the things that mathematicians do, induction works fine. But if we decide to use a non-standard variation of set theory where the axiom of foundation doesn't hold we can no longer use recursion. For example the Axiom of Extension says that two sets are equal if their elements are equal. This is a recursive definition, but it's useless in the presence of a set a such that a∈a. At this point mathematicians need a principle of coinduction. And for more on that, I refer you to Vicious Circles.

Oh...time for a quick rant. Over the years I've seen a few people argue that there's something fundamentally wrong with the notion of the algorithm because it doesn't apply to the kind of open-ended loop we see in operating systems and interactive applications. Some have even gone further to suggest that somehow mathematics and computer science are fundamentally different because mathematics can't seek to describe these kinds of open-ended phenomena. As I've tried to show above, not only are there nice ways to look at open-ended computations, but from a mathematical perspective they are precisely dual, in the category theoretical sense, to terminating computations. It may be true that mathematicians sometimes spend more time with things, and computer scientists with cothings. But this really isn't such a big difference and the same langiage can be used to talk about both.

I learnt all this from a variety of sources including Turner's paper on Total Functional Programming. My original motivation for trying to understand this stuff came from this post on Ars Mathematica. Curiously, that post had a significant effect on me. It made me realise there was this entire co-universe out there that I knew nothing about, and I was forced to make a real effort with learning Haskell because Haskell was the programming language closest to the mathematical notation used to talk about this co-universe. (Thanks Walt, but sorry to everyone who wants to see less Haskell here!)

I didn't really get why anyone would want to distinguish between data and codata until I saw some clues in some slides by Altenkirch and some comments by Conor McBride. I'm still not sure anyone says it quite as clearly as this: distinguishing between data and codata means we can allow the coexistence of infinite lists, structural recursion and open-ended loops, without risk of causing bad behaviour.

Death>Dishonour =/=> ~Dishonour>~Death

That you rate death before dishonour should not imply that you rate honour above life. Or more precisely: if A and B are (not necessarily disjoint) outcomes, define A>B to mean that the expected utility of outcome A is higher than the expected utility of outcome B. Use the symbol '~' to mean the negation of an outcome. Then Death>Dishonour does not imply that ~Dishonour>~Death, or more generally, A>B =/=> ~B>~A.

This example from the book Subjective Probability: The Real Thing by Richard Jeffrey was pointed out to me by a friend of mine yesterday. Let's draw a little diagram of possible outcomes and their utility:

B\A	death	life
dishonour	1	*
honour	2	3

The '*' indicates an outcome with probability zero and the other outcomes are considered to be equally likely. A is the variable describing your state of liveness and B is the variable representing your honour.

We have:

E(utility|A=death) = 1.5

E(utility|B=dishonour) = 1



E(utility|B=honour) = 2.5

E(utility|A=life) = 3

With this set of utilities one does indeed rate death before dishonour but not honour before life.

This gets mixed reactions. Some people think that it's a nice counterintuitive result. Some people think it's so trivial they don't see why anyone would even mention it. If you're in the latter group, sorry for wasting your time. :-)

(I leave the Haskell program using a probability monad as an exercise...)

Monads from Algebra and the the Gray Code from Groups

There is a close association between algebraic structures and monads. I've mentioned this in passing before but I think it's interesting to work out some details. Among other things, it gives us a way of converting your favourite algebraic structures into domain specific languages embedded in Haskell.

First an aside. Haskell monads are in some sense only an approximation to mathematical monads. In fact, Haskell is only an approximation to mathematics. It's easy to define a Haskell function f, say, such that x == y but f x /= f y. Once you do such a thing you can no longer reason about Haskell functions with the assumption that == represents equality in the mathematical sense. (For an example, see Saizan's comment here.) So in the following I'm going to assume that we've limited ourselves to functions such that x == y implies that f x == f y, x < y and y < z imply x < z and so on. I'd love to see some way, in Haskell, to make explicit 'promises' about the properties of functions. (I guess that's what they call programming by contract.) But for the following we'll just assume it holds.

The first abstract algebraic structure that people study is usually the group, but I'm going to start with something simpler, the monoid. A monoid is a set M (called the underlying set of the monoid) with a binary operator, ·:M×M→M and identity element e such that e·x = x·e = x and x ·(y·z) = (x·y)·z for all x,y,z in M. By abuse of notation we'll often use M for the name of the monoid as well as the underlying set. Monoids are easy to come by and there are lots of examples. Any time you have a sequence of operations that can be applied, or not, in any order, to some system, you have a monoid. Sequences of operations form the underlying set of the monad and the binary operator means "do what's on the right" followed by "do what's on the left". For example the set of uniform scaling operations we can perform on a circle form a monoid. Write a scaling by a factor of x as s(x). If we double the size of a circle, and then triple it we have s(2) followed by s(3), written s(3)·s(2). Note that we have s(m)·s(n) = s(mn) and that s(1)·s(m)=s(m)·s(1)=s(m) so s(1) plays the role of the identity, e. Also note that for all m apart from 0 we have s(m)·s(1/m)=e so every scaling apart from s(0) has in inverse. Note however that this monoid is special because for all x and y in it, x·y=y·x. This doesn't always hold in monoids and when it does hold our monoid is said to be commutative.

If we have a type of algebraic structure then we can often form "free" versions of that structure. For monoids we get free monoids. There are many different ways of looking at "freeness" and I'm going to go through some of them informally:

1. Given any monoid there are equations satisfied by its elements. From the above example we have that x·y=y·x and also that x·e·x=x·x. But notice how there is a big difference between these equations. The former doesn't hold in all monads, but the latter does. A monoid is said to be free when the only equations that hold follow from equations that hold in every monoid. You can think of a free monoid as being a generic monoid. It has no special properties above and beyond being a monoid. Given a set S, the free monoid generated by S is the smallest free monoid containing the set S and with no equations relating the elements of S. Write this monoid as FS. For example, suppose S = {x,y}. Then we know that e, x and y are all in FS. We also know that x·x, x·y, y·x and y·y are all in S. Importantly we know that all 4 of these elements are distinct because we know there can be no equations connecting them beyond those that define a monoid. In fact, it's not hard to see that the elements of FS are essentially just the (possibly zero length) strings of x's and y's.
   
2. Given a set S, consider its elements to be 'labels' for unknown variables. In fact, just use the elements of S as variables. Then we can also consider the free monoid generated by S to be the set of "monoid expressions" in these unknowns. By "monoid expression" we just mean functions of these variables that can be written using the language of monoids, ie. the symbols e and ·. If S = {x,y} then examples of such expressions are e, x, y, x·y and so on. I hope it's not hard to see that this is simply another description of the same thing as in the previous paragraph.
   
3. Another way to think of free monoids requires a tiny bit more algebra. Define a monoid homomorphism from one monoid, M, to another, N, to be a function f:M→N on the underlying sets such that f(e)=e and f(x·y)=f(x)·f(y). A bijective monoid homomorphism is called an isomorphism. If there is an isomorphism between two monoids then in some sense they are the same monoid. Note that e and · are being overloaded here - we have two different monoids and these symbols have different meanings in each one. Now we can define the free monoid generated by S to be the monoid, M, such that (1) there is a function i:S→M (2) given any monoid N, and any function f:S→N, then f can be factored as f' o i where f' is a monoid homomorphism.
   
4. Very vaguely, the set S lives in the category of sets. But FS lives in the category of monoids. FS is the closest approximation to S in the category of monoids. The previous property gives a sense of what this means: anything that can be said about a function from S to a monoid N can be said in the language of monoids using a homomorphism from FS to N. Don't worry if this is too vague as I won't be using it below - but it may click in some people's minds.
   

The important thing here is that the operator F, ie. the "free monoid generated by" operator, forms a monad. It's tempting to think that monads got their name from monoids for this reason, but I don't think it is because just about any "free X generated by" operator forms a monad. Haskell monads give rise to DSLs via do-notation. So this means that algebraic structures give rise to monads, which may give rise to Haskell monads, and hence DSLs.

So now to explain why F forms a monad. Note that just about everything I say here works with other kinds of algebraic structures ranging from monoids through groups to vector spaces. Consider our set S={x,y} above. For convenience, let's drop the · and write the binary operator as multiplication in the usual way. Given an element of S we can easily get an element of FS. In fact, we have an embedding i:S→FS with i(x) = x. This is an abuse of notation. The element on the LHS is an element of S and the x on the RHS is a variable labeled with the symbol x, but we also write this as x because v_x, or something like that, would be a pain in the ass to keep writing.

Now think about the elements of F(FS). These are strings of elements of FS, ie. strings of monoid expressions in x and y. We could write a typical member as (xyx)(e)(xxy) where the xyx, e and xxy are each elements of FS and I'm using parentheses to make clear which bits they are. It's tempting to say that this is simply xyxxxy, but that would be wrong. The former is an element of F(FS) and the latter is an element of FS. It would be clearer that these weren't equal if we used the v notation to write v_xyxv_ev_xxy. But in this case the oscurity is to our advantage. Even though (xyx)(e)(xxy) doesn't equal xyxxxy, the notation strongly suggests that we define a function that maps the first to the second. In fact, we can define a function m:F(FS)→FS that essentially just erases parentheses.

So we have functions i:S→FS and m:F(FS)→FS. Looks suspiciously like a monad. In fact, i and m satisfy the monad laws (exercise!) and make F into a monad.

So what monad is this? It's probably not hard to guess: elements of FS are finite strings of elements of S. So it's essentially the List monad. Unfortunately, Haskell allows you to form infinite lists and so the correspondence isn't 100% precise. Nonetheless, it's good enough that if you hadn't already invented the List monad (as in Haskell monad), you probably would if you had considered the free monoid monad (as in mathematical monad). i is just the embedding \x -> [x] and m is concat, which essentially just erases brackets. So you can think of [x,y,z] has the Haskell way to write the monoid element xyz.

In fact, if you repeat the above discussion with other algebraic structures that I and others have discussed you'll get other familiar monads (modulo details like being able to create infinite lists). Here's a table (with apologies for bad formatting caused by blogger.com):

Algebraic Structure	DSL
Monoid	Combinatorial search
M-Set	"Write-only" side effects
Vector Space	Probability theory/Quantum Mechanics
Semirings	Two player game tree search
Modules over tropical semiring*	Min-sum algorithm
Group	???

(Note that the above table is approximate in the sense that sometimes you may need to restrict to instances of Eq and tweak the code slightly to make the Haskell monad behave like the mathematical one.)

Notice the gap in the "group" row. A group is a monoid in which every element has a left and right inverse. The monad is the "free group generated by" monad and I'll call it G. It doesn't correspond to any of the standard Haskell monads. So firstly - if free monoids give lists, what does the underlying datatype for free groups look like? Go back to our set S = {x,y}. GS contains e, x and y as well as all of the products of x and y that appear in FS. But additionally GS contains the inverses of x and y, x^-1 and y^-1. And of course you need all strings of x, y and their inverses. But do you need more? What about (xy)^-1y? Well we can expand out the inverse using (xy)^-1=y^-1x^-1. The net effect is that the free group contains precisely all (possibly empty) strings of x, y, x^-1 and y^-1, where substrings like xx^-1 and x^-1x are removed. We can model this in Haskell using the following type:

> import Control.Monad
> data Group a = G [Either a a] deriving Show

We're using Left x to represent x and Right x to represent x^-1. We implement >>= so that it uses the inverse of product formula above. It also ought to cancel out terms like xx^-1 but to do that requires that we restrict the monad to instances of Eq and if I do that I risk castigation from Saizan. So I'll leave out the cancellation to get the following Haskell monad:

> instance Monad Group where
>     return x = G [Left x]
>     G x >>= f = G $ concatMap g x where
>         g (Left y) = let G u = f y in u
>         g (Right z) = let G u = f z in reverse (map (either Right Left) u)

But what purpose might this serve? Try these:

> test1 = sequence $ replicate 4 (G [Left 0,Left 1])
> test2 = sequence $ replicate 4 (G [Left 0,Right 1])

The first does something almost identical to what you might expect from the ordinary list monad. But note the ordering in the second. This is the Gray code monad! By using Left and Right we can control the direction in which future combinatorial searches are carried out. Admittedly not the most useful monad, but it's curious that we do get something recognisable from this abstract nonsense. And maybe you spotted the "beautiful" Haskell implementation of the power set function on reddit recently. Here's a cute twist on that:

> powerset1 = filterM (const $ G [Left False,Left True])
> powerset2 = filterM (const $ G [Left False,Right True])

The first gives a result similar to the original, but the second lists the subsets in such a way that each element differs from the previous one by one element.

I wonder what other algebarically motivated monads are waiting to be discovered.

And I hope to write something about deriving comonads from coalgebras, as soon as I've read this. It looks really intriguing. Problem is, I'm having trouble making sense of it.

And sorry, this is a bit incoherent today. I only get a short period of time to write this stuff up and I blew most of mine this week in traffic jams. As usual, if the text makes no sense, you can be sure the code above works as I just tested it. But ask questions and complain about errors...

^* The tropical semiring is the proper name for the (R,min,+) semiring I talked about earlier.

Update: Saizan's doing a great job of keeping me honest. He noticed that I'd omitted a step. I was originally finding the inverse of a free group element by simply reversing the elements in the product. But I was failing to flip the polarities of each of the elements (so to speak) despite having spelled out in the text exactly what I needed to do. The code is now fixed.

How to write tolerably efficient optimzation code without really trying...

...or the (freely generated module over the (min,+)-semiring)-monad. (But if I used that as a title you probably wouldn't be reading this.)

There's a large class of optimisation problems that can be considered to be about graph searching. This can include anything from finding the shortest distance between two cities, through finding the edit distance between two strings to finding the most plausible reconstruction of a signal corrupted by noise. In each case there's an obvious implementation of a solution in Haskell: use the list monad to generate all of the paths in the search space and then find the optimal path. The problem with this is that generating all paths can result in combinatorial explosion with exponential running times. So this post is about what I think is a novel monad to fix that problem. It gives a simple DSL for specifying a wide range of optimisation problems with polynomial instead of exponential complexity. In some cases it results in something close to the best published algorithm with little effort on the part of the user.

This is quite a neat and useful monad. So if the following is hard going, skip to the examples first for some motivation. I'd hate people to miss out on something that turns out to be fairly easy to use once it's set up.

Unfortunately, despite the expressiveness of Haskell making this possible, we have to fight Haskell on two fronts. One is the old annoyance of Haskell not supporting restricted monads. For a solution to that I'll refer you to Eric Kidd who used an Oleg trick. So let's get the first one out of the way before talking about graphs. Here's a preamble that replaces the Prelude's monad class with our own:

> {-# LANGUAGE MultiParamTypeClasses,
>            NoImplicitPrelude,
>            UndecidableInstances #-}

> import Debug.Trace
> import Prelude hiding ((>>=),(>>),return)
> import Data.Monoid
> import Data.Map (toList,fromListWith)

> class Monad1 m a where
>  return :: a -> m a

> class (Monad1 m a, Monad1 m b) => Monad2 m a b where
>   (>>=) :: m a -> (a -> m b) -> m b
>   (>>) :: m a -> m b -> m b
>   l >> f = l >>= \_ -> f

> class MonadPlus m where
>    mzero :: m a
>    mplus :: m a -> m a -> m a

So here's a graph with lengths marked on the edges:

Suppose we want to find the shortest path from A to F in 3 hops. Before solving that, consider a more general problem. In the following graph

suppose we wish to find the shortest n-hop paths from A to the C_i. No matter what complexity is in the cloud, each of the paths must travel through one of the B_i. So we can solve the shortest paths problem by finding the shortest (n-1)-hop paths to the B_i and using:

dist_n(A,C_i) = min_j(dist_n-1(A,B_j)+d_ji)

Here's the crucial point to notice: If we consider dist_n(A,B_i) to be a vector whose indices are labelled by i, and d_ji to be a matrix, then the above formula is an ordinary matrix multiplication except that instead of addition and multiplication we use min and addition. If you've studied graph algorithms you'll recognise this as the same multiplication that appears in the Floyd-Warshall algorithm. If we extend the real numbers with infinity, then (R∪∞,min,+) forms a semiring. (See Cormen, Leiserson & Rivest for some good discussion of semirings in this context.) Just as you can form vector spaces over a field, you can form modules over this semiring. And just as vector spaces form a monad, so do these modules. But you don't need to know much about semirings because the monad itself is a fairly straightforward object.

Here's that infinity I was talking about:

> infinity :: Fractional a => a
> infinity = 1/0

I have to admit, I enjoyed writing that line of code :-) More importantly, Haskell does the right thing with infinities.

The monad is going to represent the length of a path to a graph vertex. Here's the definition:

> data M l p a = M { runM :: [(a,p,l)] } deriving (Eq,Show,Ord)

Ignore the l part for now, that's going to serve as a place to make a log, like the Writer monad. We'll need that when we want the actual shortest path, not just its length. For now we'll just use an empty list. So suppose that at some stage in our computation we've reached the node B with a shortest path of length 3 and node F with a shortest path of length 4. We represent this as

u = M [(B,3,[]),(F,4,[])] :: M Vertex Float [()]

.
Suppose we now consider what happens at the next step. We can express the options from B with

v = M [(A,1,[]),(D,3,[])]

and those from F as

w = M [(D,3,[]),(E,3,[])].

To combine these with the earlier computed paths we could replace B and F in u with v and w respectively to get:

M [(M [(A,1,[]),(D,3,[])],3,[]),(M [(D,3,[]),(E,3,[])],4,[])] :: M (M Vertex Float [()]) Float [()].

We need to reduce this back down to a M Vertex Float [()] using the matrix 'multiplication' formula above to pick the appropriate shortest path each time. But note the pattern: it's the usual M M a -> M a monad pattern. So we can write the substitution+reduction as an implementation of >>=. Here's the implementation, along with some other monadic goodness:

> instance (Ord a,Monoid l,Num p) => Monad1 (M l p) a where
>   return a = M [(a,0,mempty)]

> instance (Ord a,Monoid l,Fractional p,Ord b,Ord p,Ord l) => Monad2 (M l p) a b where
>   M x >>= f = let
>        applyf = [(f a,p,l) | (a,p,l) <- x]
>        rawresult = concat [[(a',(p+p',l `mappend` l')) | (a',p',l') <- a] | (M a,p,l) <- applyf]
>        reduced = map (\(a,(p,l)) -> (a,p,l)) $ toList $ fromListWith min $ rawresult
>        filtered = filter (\(_,p,l) -> p/=infinity) reduced
>        in M filtered

> instance MonadPlus (M l p) where
>    mzero = M mempty
>    mplus (M a) (M b) = M (a ++ b)

> guard True = return ()
> guard False = mzero

> choose x = M $ map (\(a,p) -> (a,p,mempty)) x

I've also had to provide my own implementation of MonadPlus. choose is a convenience function so we don't have to keep entering those empty lists.

So here's that graph represented in terms of a function that describes transitions from vertex to vertex:

> data Vertex = A | B | C | D | E | F deriving (Eq,Show,Ord)

> graph A = choose [(B,1),(C,2)]
> graph B = choose [(A,1),(D,3)]
> graph C = choose [(A,2),(D,4),(E,1)]
> graph D = choose [(B,3),(C,4),(E,2),(F,3)]
> graph E = choose [(C,1),(D,2),(F,3)]
> graph F = choose [(D,3),(E,3)]

And here's our first example. We compute all of the shortest 3-hop paths from A:

> ex0 :: M [()] Float Vertex
> ex0 = do
>    step1 <- return A
>    step2 <- graph step1
>    step3 <- graph step2
>    step4 <- graph step3
>    return step4

Just in case you think that's the wrong result note that these are the shortest paths to each accessible vertex that are exactly 3 hops long. If you want paths up to 3 hops long you have to add some zero length transitions from nodes to themselves. I'll leave that as a (trivial) exercise.

But suppose you just want the shortest path to F. Then you can use:

> ex1 :: M [()] Float Vertex
> ex1 = do
>    step1 <- return A
>    step2 <- graph step1
>    step3 <- graph step2
>    step4 <- graph step3
>    guard $ step4==F
>    return step4

Note how simple this code is. We just write code to walk the graph. All of the optimisation machinery is hidden in the background. But there's a dark secret here because the code isn't yet doing what a claim. But before I admit what the problem is, let's modify this code so we can actually see what the shortest paths are. record is just like tell from the Writer monad. I've implemented >>= so that it does the right thing:

> record l = M [((),0,l)]

> ex2 = do
>   step1 <- choose [(A,0),(B,0)]
>   record [step1]
>   step2 <- graph step1
>   record [step2]
>   step3 <- graph step2
>   record [step3]
>   step4 <- graph step3
>   record [step4]
>   guard (step4==E || step4==F)
>   return (step1,step4)

Just for fun I made this example more complex. It finds all four shortest paths from anywhere in {A,B} to anywhere in {E,F}. As an exercise try modifying the above code so that it finds the one shortest path from {A,B} to {E,F}. It's also simple to modify it to deal with conditions like the traveller being at one of a particular set of nodes at a particular time. It's a DSL, you have considerable freedom to form whatever conditions you want.

And now it's time to reveal the secret. Although I wrote >>= so that it reduces the list of considered paths at each stage the above code still has exponential complexity as you'll rapidly discover as you increase the number of steps.

Consider the third monad law: (m >>= f) >>= g = m >>= (\x -> f x >>= g).

If we evaluate the left hand side for this monad we get a reduction after applying f and then another after applying g. But we if we evaluate the right hand side we find that we get many small (and useless) reductions applied and then at the end a large reduction is applied. We don't get a big reduction at the intermediate step. The result: exponential complexity because the code ends up exploring every possible path performing stupid trivial reductions at each step. The monad law says we get the same result, but it doesn't take the same time.

And here's the really annoying thing: when you use do-notation in the usual way the Haskell compiler builds expressions like the slow right hand side. Bummer! But there's an easy, if tedious, fix. We have to use the third monad law as a rewrite rule to convert our code to the form on the left hand side. Here's ex2 rewritten:

> ex3 = do
>   step1 <- choose [(A,0),(B,0)]
>   step4 <- do
>       step3 <- do
>           step2 <- do
>               record [step1]
>               graph step1
>           record [step2]
>           graph step2
>       record [step3]
>       graph step3
>   do
>       record [step4]
>       guard (step4==E || step4==F)
>   return (step1,step4)

Not as pretty is it? It's still not too bad but it'd be nice to make this conversion automatic. I've a hunch this can be done using Haskell rules. Anyone care to volunteer some help making this work?

Anyway, I mentioned that lots of problems are in this class, some not obviously graph theoretical. Here's one: computing the edit distance between two strings:

> diag (a:as,b:bs) = choose [ ((as,bs),if a==b then 0 else 1) ]
> diag _           = mzero

> right (a:as,bs) = choose [((as,bs),1)]
> right _         = mzero

> down (as,b:bs) = choose [((as,bs),1)]
> down _         = mzero

> nop ([],[]) = return ([],[])
> nop _       = mzero

> munch :: (String,String) -> M [()] Float (String,String)
> munch x = diag x `mplus` right x `mplus` down x `mplus` nop x

> pure (M a) = length a==1
> editDistance a b = until pure (>>= munch) (munch (a,b))

> ex4 = editDistance "hypocrisy" "hippocracy"

I think this is actually the dynamic programming algorithm in disguise (and slightly reordered) with some extra log time penalties because it's using a binary tree (from Data.Set) instead of arrays. But I wouldn't stake my life on it.

Anyway, here's one last example. Rather than explain it I refer you to the nicely written Wikipedia on the Viterbi algorithm. In fact, I decided to lift their example precisely. And note that I've (1) written the code twice, the second one is the one rewritten for speed, (2) used logarithms to convert the multiplications to additions and (3) 'unrolled' the loop to make absolutely explicit what's going on. If I wrote Haskell code to work on general arrays like the Python code in the Wikipedia article, it'd be too short and you might miss it if you blinked.

> data Weather = Rainy | Sunny deriving (Eq,Show,Ord)
> data Activity = Walk | Shop | Clean deriving (Eq,Show,Ord)

> transition Rainy = choose [(Rainy,-log 0.7),(Sunny,-log 0.3)]
> transition Sunny = choose [(Rainy,-log 0.4),(Sunny,-log 0.6)]

> emission Rainy = choose [(Walk,-log 0.1),(Shop,-log 0.4),(Clean,-log 0.5)]
> emission Sunny = choose [(Walk,-log 0.6),(Shop,-log 0.3),(Clean,-log 0.1)]

> ex5 = do
>       day1 <- choose [(Rainy,-log 0.6),(Sunny,-log 0.4)]
>       record [day1]
>       emission1 <- emission day1
>       guard (emission1 == Walk)
>       day2 <- transition day1
>       record [day2]
>       emission2 <- emission day2
>       guard (emission2 == Shop)
>       day3 <- transition day2
>       record [day3]
>       emission3 <- emission day3
>       guard (emission3 == Clean)
>       day4 <- transition day3
>       record [day4]

> ex6 = do
>       day4 <- do
>           day3 <- do
>               day2 <- do
>                   day1 <- choose [(Rainy,-log 0.6),(Sunny,-log 0.4)]
>                   do
>                       record [day1]
>                       do
>                           emission1 <- emission day1
>                           guard (emission1 == Walk)
>                   transition day1
>               do
>                   record [day2]
>                   do
>                       emission2 <- emission day2
>                       guard (emission2 == Shop)
>               transition day2
>           do
>               record [day3]
>               do
>                   emission3 <- emission day3
>                   guard (emission3 == Clean)
>           transition day3
>       record [day4]

[Sunny,Rainy,Rainy,Rainy]. It's the Viterbi algorithm for free.

Anyway, I haven't proved a thing above. I think (hope?) I'm right about the polynomial complexity. Certainly on paper it all works, but I've always found it hard to reason about the complexity of Haskell code.

Unfortunately this code is all still slower than I'd like. I've played with trace during the development and everything seemed to run the way I expect (at least after I discovered the third monad law problem). I think that ultimately there's a big overhead for piling up HOFs, and the lookups in Data.Set probably aren't helping. Still, I don't think that invalidates the mathematics behind the approach and I think that asymptotically the slowness is really a big constant and a bunch of logs and it's still much better than combinatorial search. I ought to see how ocaml fares.

BTW There's a deep connection between this monad and the quantum mechanics one. I've no time to write about that now, though see here.

Some acknowledgements: the strange but highly readable Wikipedia article motivated me to think about the Viterbi algorithm, and I became sure there was a monad lurking in the Python code. I spent a month or so getting nowhere. And then I spoke with Eric Kidd on the phone about probability monads. He challenged me on something and in an attempt to formulate a response the whole of the above seemed to just crystallise. Thanks for asking the right question Eric!

I think there might be a comonad for working backwards through the graph...

Anyway, one last thing: please for the love of the FSM will someone put a clean way to implement restricted monads into their Haskell compiler. There is so much neat stuff that could be implemented with it and I now officially consider Haskell broken without one.

Haskell Incarnate

I'm off to the Make Faire in a few minutes and I'm running in Bay to Breakers tomorrow so I probably won't get much time to write about what I've been up to this weekend. So here's a little taster. My first robot programmed in Haskell:



It's a pretty simple device that's programmed to follow objects that come into the view of the three IR distance sensors mounted on the front of it. (The object in this case being a skull on the end of a stick, must have been left over from Hallowe'en...) But the point, at this stage, wasn't to do something sophisticated with the robot, but to see if I could actually program such a device in a high level language that wasn't C/C++ or assembler.

For hardware I used the Make controller and a Rigel robot base. I also use a bunch of Sharp IR sensors.

So for software I started with Matthew Naylor's module for Lava described in The Monad Reader Issue 7. I didn't have to do much with this code at all so it's almost embarassing to claim I did anything. The Make controller is ARM based and programmable in C. In fact it can run FreeRTOS and I simply modified Matthew's code to spit out a C function that is spawned as a task in FreeRTOS. It was about 30 minutes work to make the changes. With Matthew's code modified I now have a Haskell DSL that is translated to C for compilation for the robot. Strictly speaking there isn't actually a Haskell runtime on the robot. Think of the Haskell as a metaprogram that spits out C.

The actual Haskell DSL code looks something like this:

> main =  writeMakeController "rigel" f    where
>         f inp = Output {
>            servoSpeed0   = 545+leftSpeed,                       
>            servoSpeed1   = 550-rightSpeed,
>            servoSpeed2   = 550+leftSpeed,
>            servoSpeed3   = 540-rightSpeed
>          }
>             where
>                 rightSpeed = larger (threshold 300 leftSensor) (threshold 300 centreSensor)
>                 leftSpeed = larger (threshold 300 rightSensor) (threshold 300 centreSensor)
>                 centreSensor = smooth $ adc0 inp
>                 leftSensor = smooth $ adc5 inp
>                 rightSensor = smooth $ adc4 inp
> smooth x = (x + delay 0 x + delay 0 (delay 0 x) + delay 0 (delay 0 (delay 0 x)))/4
> threshold m x = x />=/ m ? (256,0)
> larger x y = x />=/y ? (x,y)

And the video above is the resulting behaviour.

Note, this is the first thing I got working. I make no great claims for the quality or ingenuity of the code. Most of it should be self-explanatory but note the function smooth which averages over the last inputs. This doesn't make sense when considered as an ordinary function because in a pure functional language you can't add the argument of a function to the value that you passed in last time you called the function. But these values are all streams and in principle this makes the DSL extremely powerful. Anyway, I'll write it up better at a future point and I'll see what I can do about making the code available.

Homeland Security Threat Level Monad

I've previously introduced the idea of a monad being a kind of one way wrapper. The Monad interface allows you to wrap things, but not unwrap them again. And it allows you to apply any function you like to a wrapped object, as long as it stays wrapped. I also indicated how this could be used to model a notion of taint. If you have some data from a tainted source - for example it's unreliable, then the trivial monad can do a nice job of keeping track of all other data that has been tainted by contact with this data. What I also want to show is that a wide variety of other monads can be interpreted in terms of taint. What this means is that despite its simplicity, the trivial monad is always worth keeping in the back of your mind because it is a protoype for many other monads.

You could imagine using this taint idea as a security model. Any data of type W a is considered secret and a threat to public safety if it is revealed. By sticking to the Monad interface we force the compiler to make sure that any data depending on secret data is itself marked secret. This makes it easy to keep secret data sequestered - at least as long as people only use the Monad interface to W.

But suppose this isn't good enough and we need to track not only which data is secret, but also the threat level to public safety should a secret be revealed. With this in mind, define the Homeland Security Threat Level Monoid:

> data Threat = LOW | GUARDED | ELEVATED | HIGH | SEVERE
>          deriving (Eq,Ord,Enum,Show)

It's a monoid because it has an associative binary operator max and has an identity LOW. For example max LOW ELEVATED = ELEVATED.

We want to tag some of our data with a threat level. In addition, we want to ensure that if we combine two tagged pieces of data, x and y, the combination must be tagged with the worse of the two threat levels. That's why we have the monoid structure that allows us to combine two threats.

So define

> data T a = T Threat a deriving Show

But to use this it appears we have to thread use of the max function all the way through our code. It'll be a nightmare to implement. But there's a solution. What we want is a monad!

Consider LOW to be the default threat level of all data and that tells us how to define return.

Now suppose we have a function f :: a -> T b and some data x :: T a. We want to be able to apply f to x using the monad interface: x >>= f. The threat level of the result should be the max of the two possible sources of threat: x itself or the result of computing f. With that in mind:

> instance Monad T where
>     return x = T LOW x
>     (T t x) >>= f = let T t' x' = f x in T (max t t') x' 

Notice how we've now shown how to implement threat combination in a way that can be used by any code that uses the monad interface. So now we can take last week's definition of h and use it:

> g x y = y >>= (return . (+x))
> h x y = x >>= (\x -> g x y)

Try evaluating h (T ELEVATED 3) (T GUARDED 6). It just works, even though g and h were originally written for the trivial monad! So it's trivial to make code that's written monadically work with our new threat level system. This should give some idea of the power of monads. We've completely changed the sematics of our monadic addition operator, h, without lifting a finger to edit it.

But notice that I didn't need to talk about threats at all. I've shown how to tag data with other data that's implicitly carried along with it. All we needed was a rule for combining pairs of tags. Our tags were of type Threat. Instead we could have used strings and instead of max we could have used ++. If we had, then we'd have a way of implicitly carrying around logs with our computations. In fact, all I have done here is implement a special case of the Writer monad. But I still think that thinking in terms of data that is quarantined in some way helps with the intuition.

I hope that by starting with the trivial monad last week, and upgrading in a natural way to the Homeland Security Threat Level monad, I've given an easy way to think about monads. Our ordinary intuitions about taint, for example that if a clean object touches a tainted one it becomes tainted itself, carry over to monads. And the Writer monad isn't the only one that can be thought of in these terms. I hope at some point soon to look at the Maybe and list monads in this way too.

The Trivial Monad

Update: I've added some solutions to exercises below.

Haskell monads form a type class. And Haskell type classes are essentially interfaces shared by some types. So the Monad type class is a common API for talking to a bunch of different types. So the question is this: what's so special about this API?

One way to grasp an API is to see concrete examples. Usually the monad API is introduced through containers or IO. But I think that the prototypical monad through which all others can be understood is much simpler - it's the trivial monad.

APIs often capture design patterns, and the design pattern here is a one-way wrapper. We can wrap things but we can't unwrap them. We still want to be able to whatever we like to the wrapped data, but anything that depends on anything wrapped should itself be wrapped.

Without further ado, here's a wrapper type:

data W a = W a deriving Show

(The Show bit is just so we can play with these things interactively.)

Note how it doesn't add anything except some wrapping. And the first thing we need is to be able to wrap anything we like. We can do that easily with this function:

return :: a -> W a
return x = W x

(We could have just used W instead of return. But I'm heading towards a common API that will work with other types too, so it obviously can't be called W.)

And now we need one more thing - a way to manipulate wrapped data leaving it wrapped. There's an obvious idea. Given any function a -> b we write a function that converts it to a function W a -> W b. This is guaranteed to keep things under wraps. So here goes

fmap :: (a -> b) -> (W a -> W b)
fmap f (W x) = W (f x)

.

It seems that we've finished. For example we can wrap a number and then increment it:

a = W 1
b = fmap (+1) a

But here's something we can't do. Define a function f like this:

f :: Int -> W Int
f x = W (x+1)

It increments and returns a wrapped result. Now suppose we want to apply the underlying operation here twice, ie. we'd like to increment a number twice and return a wrapped result. We can't simply apply f twice because (1) it'll doubly wrap our result and (2) it'll attempt to add 1 to a wrapped result, something that doesn't make sense. fmap doesn't do what we want either. Try it in an interactive Haskell session. It seems we need a way to apply f, unwrap the result and apply f again. But we've already said that we don't want to allow people to unwrap these things. We we need to provide a higher order function that does the unwrapping and application for us. As long as this function always gives us back something that is wrapped, our end users will never be able to unwrap anything. Here's an idea for such a function:

bind :: (a -> W b) -> (W a -> W b)
bind f (W x) = f x

Notice how it's very similar to fmap but is even simpler. So now we can try doubly, or triply, incrementing

c = bind f (f 1)
d = bind f (bind f (f 1))

Notice how bind is more general than fmap. In fact, fmap f = bind (return . f).

And that's it. Using return and bind we have achieved our goal of wrapping objects and freely manipulating wrapped objects while keeping them wrapped. What's more, we can chain functions that wrap without getting bogged down in multiple layers of wrapping. And that, really, sums up what a Haskell monad is all about.

So here are a couple of exercises:

(1) define a function g :: Int -> W Int -> W Int so that g x (W y) = W (x+y). Obviously that definition won't do - the left hand side has a W y pattern so it's actually unwrapping. Rewrite this function so that the only unwrapping that happens is carried out by bind.

(2) define a function h :: W Int -> W Int -> W Int so that h (W x) (W y) = W (x+y). Again, no unwrapping.

I'm hoping that after you've done these exercises you'll see how you can still work freely with data even though it's wrapped.

In Haskell, it is more usual to use the operator >>= instead of bind where bind f x = x >>= f.

So the last question is this: why would you ever wrap data like this? In practice people tend not to use the trivial monad very much. Nonetheless, you can see how it might be used to represent tainted data. Wrapped data is considered tainted. Our API never lets us forget when data is tainted and yet it still allows us to do what we like with it. Any time we try to do anything with tainted data the result is also tainted, exactly as we might expect. What I find interesting is that almost every monad, including IO, lists and even probability, can be thought of quite naturally as variations on taint. I hope to say more about this in the near future.

Anyway, code above fails to compile because of some namespace clashes. Here's a complete definition of W that really works. Note also that fmap, which we showed we don't really need, allows us to make W an instance of Functor too:

> data W x = W x deriving Show

> instance Functor W where
>    fmap f (W x) = W (f x)

> instance Monad W where
>    return x = W x
>    W x >>= f = f x

Exercise 3: Prove the three monad laws for W. This should be almost trivial.

Exercise 4: We can't completely unwrap things using the monad API. But we can unwrap one layer from things that are wrapped twice. So here's a nice puzzle: define a function join :: W (W a) -> W a using the Monad API and no explicit unwrapping.

In conclusion, most monad tutorials show what monads are, and how to use them. I hope I've given some additional insight into why the Monad interface consists precisely of the two functions it does. That insight should become clearer when I say more about taint.

PS I've been a bit busy lately with planning for kitchen remodelling, lots of work and a new cat(!). I probably won't be posting very frequently for a few months.

---

Solutions to Exercises

Some people have found the exercises tricky so I'm posting some solutions. I think it's a good sign that they're tricky, it means that there is some non-trivial stuff to be learnt from playing with this monad. It ought to figure more strongly in tutorials. In a sense it allows people to learn about monadicity in its purest form.

Anyway, enough of that.

Exercise 1: We want g :: Int -> W Int -> W Int so that g x (W y) = W (x+y). Start the definition like this

g x y = ...

We want to apply (+x) to y. But y is wrapped so we must use >>= to apply it. But (+1) doesn't have the right type to be applied with >>=, it's an Int -> Int not an Int -> W Int. So we must compose with return and we get

g x y = y >>= (return . (+x))

Exercise 2: We want h :: W Int -> W Int -> W Int so that h (W x) (W y) = W (x+y). But g is already fairly close. So again write

h x y = ...

The difference from g is that x is wrapped. So we want to apply \x -> g x y to our wrapped value. This function is already of type Int -> W Int. So we can just write

h x y = x >>= (\x -> g x y)

or just h x y = x >>= flip g y.

Exercise 3: I'll just do the last one: (m >>= f) >>= g = m >>= (\x -> f x >>= g)


m is of the form W x for some x:

(W x >>= f) >>= g 
= (f x) >>= g
 
W x >>= (\x -> f x >>= g)
= (\x -> f x >>= g) x
= f x >>= g

So the LHS equals the RHS.

Incidentally, you can get a long way with these problems by not thinking about what's going on! Instead, just write out the type signatures of all of the functions involved and try to stick them together like a jigsaw puzzle so that the final result has the right type. At each stage the arguments of each function must match the signature of the function so it really is like fitting the shapes of jigsaw pieces together. This isn't a foolproof strategy, but it often works.

Anyway, feel free to ask questions if you're still stuck. Some of your questions may be answered here.

The curious rotational memory of the Electron, Part 2

A spinning top keeps going because it has angular momentum, and the angular momentum arises because the top is made up of mass that is distributed around the axis of rotation. The angular momentum is represented by a vector (well, a pseudovector really) directed along the axis of spin (at least when it's spinning around the top's axis). Similarly electrons have angular momentum. It's tempting to think of electrons as like microscopic spinning tops but the angular momentum of electrons isn't due to an extended mass rotating about some axis. In fact, as far as we know, electrons have no extension. Instead their angular momentum is a peculiar kind of intrinsic angular momentum that has a lot in common with ordinary angular momentum - eg. it's conserved and it's interconvertible with ordinary angular momentum. Because of it's not being associated with a mass distribution it's given the special name of spin, but as long as you bear in mind these provisos it's probably fair to use at least some intuitions about spinning tops when thinking about electrons.

But it's a little more complex than this. We also have quantum mechanics to contend with. The spin of an electron is a vector. But we find that when we measure one of the components of this vector this value is quantised and can only take values +hbar/2 and -hbar/2, where hbar is Planck's constant. We choose units where h-bar is 1 so the z-component of the spin is always measured to be +1/2 or -1/2. If we write these two states as |+> and |-> then because we are dealing with quantum mechanics, the z-component of the spin can be represented by the linear combination a|+>+b|->. This corresponds to a state in which there is a probability |a|² of measuring +1/2 and a probability |b|² of measuring -1/2. This is what might have been written as a.*return (1/2)+b.*return (-1/2) in my earlier Haskell code. But that's just one component of the spin. What about the x- and y-components? Amazingly the state a|+>+b|-> tells us everything we can possibly know about the spin of an electron and we'll call it a spin state.

Suppose we have an electron in the state ψ = a|+>+b|->. What happens if we measure the y-component of its spin? One way to answer that question is to rotate the electron through π/2 so that its x-axis is rotated to align with the z-axis and then measure the z-component of its spin. In order to to that we need to know how to rotate spin states. The rule for rotation through θ about the x-axis is this (in a suitable coordinate frame):

|+> → cos(θ/2)|+>-sin(θ/2)|->
|-> → sin(θ/2)|+>+cos(θ/2)|->

Note how choosing θ=0 gives the identity, as expected. Note also that θ=π maps a|+>+b|-> to b|+>-a|-> so that the probabilities of measuring +1/2 and -1/2 are simply swapped, exactly what you'd expect for turning a state upside down. But there's something else that you should notice - there's an ambiguity. A rotation through 2π should give the same as a rotation through 0 and yet setting θ=2π in that transformation maps a state ψ to -ψ. Now |a|² = |-a|² so the probability of observing spin up or spin down is unaffected. But as I've been showing over previous posts, flipping a sign in a state can make a big difference as soon as you start performing interference experiments. The same goes for any angle: if I rotate through π should I use θ=π or θ = 3π? So can the transformation I've given make sense?

The transformation does make sense if you consider that in any physical process that rotates an electron the transformation will evolve continuously over time. Electrons don't just instantly rotate. In other words, if a rotation is applied to an electron then it will follow a path in SO(3), not just be an instantaneous application of an element of SO(3). And that allows us to resolve the ambiguity: the rotations of electrons are described by the double cover of SO(3) known as SU(2). So a rotation through 360 degrees doesn't return you to the identity although a 720 degree rotation does. The transformation I gave above is completely unambiguous if you continuously rotate an electron around the x-axis tracking a continuous value of θ, after all, the double cover is basically just the set of continuous paths from the identitiy in SO(3) (with homotopic paths considered equivalent).

And that's the bizarre fact: electron rotations aren't described by SO(3), they're described by SU(2). In particular, rotating an electron through 360 degrees does not return it to its original state, but a rotation through 720 degrees does! In a sense, like Dirac's belt, electrons can remember something about the path they took to get where they are, in particular they remember how many twists there were in the path.

What does this mean experimentally? the first thing to note is that this is true not just for electrons but any spin-1/2 fermion. This included protons and neutrons. The stuff I've been talking about manifests itself in a number of ways. In particular, the spin of a particle affects how a magnetic field acts on it. For example, spin-up and spin-down particles can be separated into distinct beams using Stern-Gerlach apparatus. Also, the spin of particles precesses in a magnetic field and this is used on a regular basis in NMR. These two facts allow us to easily manipulate and measure the spin of fermions. In other words, the fact that fermions remember how many twists there are in their rotations isn't just some esoteric nonsense, it's now engineering and the theory is tested repeatedly all over the world.

Every familiar object is invariant under rotations through 360 degrees. So the fact that electrons need to be rotated through 720 degrees to return them to their original state seems like one of the most bizarre facts about the universe I know of. And yet many books that introduce spin just slip in this fact in a routine way as if it were no different to any other.

The fact that the biggest connected cover of SO(3) is the double cover puts a big constraint on the kinds of weird effects like this can happen. We can have a 360 degree rotation multiply by -1, but not by i, because a 720 degree rotation absolutely has to return us to where we started from. But suppose the universe were 2-dimensional. If you remember what I said about SO(2) you may notice that no such constraints apply because SO(2) has an infinite cover. There is a group in which all of the rotations through 360n degrees are distinct for distinct n. This means that a physical system could have its state multiplied by any factor (of modulus 1) when rotated through 360 degrees. Particle that behave this way are called anyons. But we live in a 3D universe so we don't expect any fundamental particles to have this property. However, in quantum mechanics any kind of 'excitation' of a physical system is quantised and can be thought of as a type of particle. These are known as quasiparticles. For example, just as light is made of photons, sound is also quantised as phonons. In the right kind of solid state medium, especially those that arise from some kind of 2D lattice, it seems quite plausible that anyons might arise. This gives rise to the so called fractional quantum hall effect. Anyons might one day play an important role in quantum computing via topological quantum computation.

The Curious Rotational Memory of the Electron, Part 1

There's a curious and bizarre fact about the universe that is introduced in physics courses without anyone stopping to point out just how curious and bizarre it is. I think it deserves to be more widely known to non-physicists. It's a little esoteric, but not so esoteric it can't be grounded in ordinary everyday concepts. And it isn't something hypothetical but something that is fundamental to modern physics and testable in the lab. It is also the reason why I chose to use theta/2 instead of theta in my quantum computing code.

This post will be in two parts. This part will mainly be about geometry and topology, and in the second part I'll talk about the physics.

In order to explain the idea I want to introduce the concept of a universal cover. If you look at the definition of this notion anywhere it can seem pretty daunting to non-mathematicians, but in actual fact the key idea is so simple that even a child could understand it. I can say that without fear of exaggeration because universal covers used to play a role in a childhood fantasy of mine. (OK, I was a weird kid, but not that weird.)

My crazy idea as a kid was this: if there were two ways to get from A to B, maybe the B you reached by taking either of the two routes were actually slightly different, even though they looked the same. I was actually more concerned with journeys from A to A. Suppose the kitchen has a door at each end and that you start in the bedroom, go down to the kitchen at one end, come out the other end, and return by a different path to the bedroom. How do you know this is the same bedroom? (I also admit I read a lot of science fiction as a kid.)



Suppose instead you took a simpler path: you left the bedroom and went to another room and then returned by the path you went. If you had taken a piece of string and tied one end to the bedroom and walked this path playing out the string behind you as you went then when you returned to your bedroom you would have both ends of the string visible. You could then reel in the string from your journey, without letting go of the ends, and see that indeed, both ends were in fact the ends of the same string and so you had retuned to the start of your journey. But if you'd taken the round trip through the double-doored kitchen you wouldn't be able to pull the string tight and so you could no longer be sure.

In what follows we consider two paths from A to B equivalent in this way: mark the two paths with strings going from A to B. If you can pull one string to line up exactly with the other (you're allowed to stretch) without breaking the string or moving the endpoints then they are considered equivalent. (The correct technical term is homotopic.) So if a path that appears to go from A to A is homotopic to the trivial zero length path, then you know that the path really is from A back to A. But if you can't shrink the path down to zero length you're faced with the possibility that the string is extended from A to a similar looking place, A'. (Of course this implies a considerable conspiracy with someone tying a distinct identical string to an identical bedroom, but if you can't see both ends, alien beings from a higher dimension could do anything...)

The best case scenario is that all of the bedrooms are the same - a conventional house. Call this house X. But what's the worst case scenario? How many different copies of the bedroom might there be? Well let's call the starting bedroom 0. If we walk anti-clockwise round the house, going through the kitchen, we end up back at a bedroom. In the worst case scenario this will be another bedroom. Let's call it 1. Walk round again and we get to another copy of the bedroom, call it 2 and so on. So we have an infinite sequence of bedrooms, (0,1,2,…). But we could have walked clockwise. In that case we'd get a sequence (-1,-2,…). So now think about the set of all points in this infinite house. Each point can be described as a pair (x,n) where x is a point in the house X and n labels which of the infinite number of paths you could have taken to get there. (Remember, if the two paths are homotopic we consider them to be the same path.)

If using a house as an example seems too far-fetched think of a mult-storey car park. It's easy to walk what you think is a closed loop but find that you've in fact walked to almost identical looking, but different, floor.

(Incidentally, this kind of infinite house is the kind of thing that can be created in a Portal engine such as in the game Portal.)

Let's think about this more generally. Given a topological space X we can form another space X' from it by (1) picking a point p in it and (2) forming the set of all pairs (x,n) where x is a point in X and n is a path (up to homotopy) from p to x. X' is known as the universal cover of X.

Enough with houses, let's try a more purely geometric example. Consider the space of rotations in 2D, otherwise known as SO(2). A rotation can simply be described by an angle with the proviso that a rotation through an angle θ is the same as a rotation through θ+2nπ where n is an integer. The space of rotations is essentially the same thing as the circle (by circle I mean the 1-dimensional boundary of a 2-dimensional disk). Looping once round the circle brings you back to where you started. So how can we visualise the universal cover of SO(2)? There's a nice trick for this. Imagine a disk that's free to rotate around one axis. SO(2) describes the rotations of this disk. But now imagine a belt, fixed at one end, connected at the other end to the centre of this disk like this:



The rule is that we're allowed to twist the belt as much as we like as long as the centre line of the belt is along the axis. Notice how if we rotate the disk through an angle of 2nπ the disk returns to its starting point, but the belt 'remembers' how many twists we took to get there. Notice another thing: if we travel along the belt from one end to the other, the orientation of the belt at each point gives a 2D rotation. In other words, the belt itself represents a path in SO(2) with one one representing a fixed point in SO(2) and the other end matching the current orientation of the disk. So if the orientations of the disk represent elements of SO(2), the possible twistings of the belt represent elements in the universal cover of SO(2). And what the twists do is allow us to differentiate between rotations through θ and rotations through θ+2nπ. In other words, the universal cover is just the real line, but not folded down so that θ is the same as θ+2nπ.

The universal cover is the worst case scenario where all non-homotopic paths are considered to represent distinct points. But it's also possible to construct an in-between scenario where some paths that aren't homotopic to each other are still considered to lead to the same place. For example, take the real line again. If we identify 0 and 2π we get SO(2). If we consider 0 and 2nπ to be distinct for n≠0 we get the universal cover. But we pick an N≠1 and consider 2Nπ to be the same as 0, but that 2nπ is distinct from 0 for 0<n<N. What we now end up with is an N-fold cover of SO(2). Unfortunately I don't know a way construct a physical example of a general N-fold cover using belts and discs. (At least for N>2. For N=2 you should easily be able to come up with a mechanism based on what's below.)

With the examples so far, the universal cover of a space has always had an infinite number of copies of the original space. That isn't always the case and now we'll meet an example where the universal cover is a double cover: the real projective plane, RP². The easiest way to visualise this is to imagine the game of Asteroids played on a circular TV screen. In the usual Asteroids, when you go off one side of the screen you reappear on the opposite side. In this version it's much the same, points on the edge that are diametrically opposite are considered to be the same points. Such a universe is a real projective plane. Let P be the point at the centre of the screen. Suppose we attach one end of our string to P and fly off towards the edge, reappearing on the other side. Eventually we get back to P again. Can we shrink this path down to a point? Unfortunately, any path that crosses the edge of our screen will always cross the edge no matter how much we try to deform it. If we move the point at which it crosses one edge of the screen, the point at the other side always remains diametrically opposite.




This means you can never bring these points together in order to collapse down a loop. On the other hand, suppose we loop twice round the edge. The following diagram shows how this path can be shrunk down to a loop:



So given two distinct points in RP² there are two inequivalent paths between them, and so the universal cover of RP² is made from two copies of the space and is in fact a double cover.

Now we'll try a 3D space - the space of 3D rotations, otherwise known as SO(3). This space is 3-dimensional because it takes 3 numbers to define a rotation, for example we can specify a rotation using a vector whose direction is the axis of rotation and whose length is the angle of rotation. Now all possible rotations are rotations through an angle 0≤θ≤π for some exis. So every rotation can be represented as a point in the ball of radius π. But a rotation around the axis a through an angle of π is also the same as a rotation around the axis -a through π. In other words, SO(3) is the closed 3-dimensional ball with diametrically opposite points identified. This is very similar to the situation with RP². In fact, almost exactly the same argument shows that the universal cover of SO(3) is the double cover. (For an alternative wording, see this posting.) There is also another way of seeing this that is analogous to the disk with belt example above. This time we have a sphere instead of a disk. The belt is fixed at one end and the other end is attached to the sphere. We are free to rotate the sphere however we want, and we're no longer confined to having the belt along one axis. It is an amazing but easily demonstrable fact that if you rotate the sphere through 2π, you cannot untwist the belt, but if you rotate it through 4π you can! This is known as the Dirac Belt Trick, and if you don't believe it, watch the Java applet written by science fiction author Greg Egan.

If you've done much computer graphics you'll have already met the double cover of SO(3). Computer graphics people sometimes like to represent rotations using unit quaternions. The unit quaternions form a 3-dimensional sphere but there are two such quaternions for each rotation. So the unit quaternions actually represent this double cover, and SO(3) is a 3-dimensional sphere with antipodean points identified.

And that's the geometry dealt with. In the next installment I'll say what this has to do with electrons and my quantum computing code.

---

Footnote

Balls and Spheres. The closed n-dimensional ball is the set of points (x₁,...,x_n) such that x₁²+...+x_n²≤1. The (n-1)-dimensional sphere is the set of points (x₁,...,x_n) such that x₁²+...+x_n²=1.