Transforming a comonad with a monad

I can't yet read most of this paper on Combining a Monad and a Comonad because I've no experience with 2-categories. Nonetheless, there's something useful I can extract from the bits that do make sense - the idea of combining a monad with a comonad in the way that we already combine monads using monad transformers. In my previous post I hinted at the idea that the cobind operation for a comonad was a bit like an operation on a SIMD computer that allows every processor to get access to information local to every other processor. The idea now will be to allow those processors to also exploit monads to perform operations like IO.

First let me first reintroduce the code for the array comonads I considered before so that this post again becomes a self-contained literate Haskell post:

> {-# OPTIONS -fno-monomorphism-restriction #-}
> {-# LANGUAGE MultiParamTypeClasses,FlexibleInstances,Arrows #-}

> import Data.Array
> import Data.Foldable as F
> import Control.Monad
> import Control.Arrow

> data Pointer i e = P { index::i,  array::Array i e } deriving Show

> bind f x = x >>= f

> class Functor w => Comonad w where
>    (=>>) :: w a -> (w a -> b) -> w b
>    x =>> f = fmap f (cojoin x)
>    coreturn :: w a -> a
>    cojoin :: w a -> w (w a)
>    cojoin x = x =>> id
>    cobind :: (w a -> b) -> w a -> w b
>    cobind f x = x =>> f

> instance Ix i => Functor (Pointer i) where
>  fmap f (P i a) = P i (fmap f a)

> instance Ix i => Comonad (Pointer i) where
>  coreturn (P i a) = a!i
>  P i a =>> f = P i $ listArray bds (fmap (f . flip P a) (range bds))
>      where bds = bounds a

Remember how things work with monads. If m is a monad, then given a function a -> m b we can lift it to a function m a ->m b. Similarly, with a comonad w we can lift a function w a -> b to a function w a -> w b. In fact, I wrote about all this a while back.

Just like before, we can write a simple 'blur' function, except this time it also performs some I/O.

> wrap i = if i<0 then i+5 else if i>4 then i-5 else i

> blur (P i a) = do
>       let k = wrap (i-1)
>       let j = wrap (i+1)
>       let s = 1.0*a!k + 2.0*a!i + 1.0*a!j
>       print $ "sum = " ++ show s
>       return s

Notice that blur is of the form w a -> m b. So how can we compose these things?

Consider g :: w a -> m b and f :: w b -> m c. Using bind and cobind we can 'lift' the head and tail of these functions:

bind g :: w a -> w (m b) and cobind f :: m (w b) -> m c.

We can almost compose these functions. But the catch is that we need a function from w (m b) to m (w b). Conveniently, there's a prototypical example of such a thing in the ghc libraries, the function sequence which distributes any monad over []. As we can easily move back and forth between arrays and lists we can write:

> class Distributes m w where
>    distribute :: w (m a) -> m (w a)

> instance (Monad m,Ix i) => Distributes m (Pointer i) where
>    distribute (P i ma) = do
>        let bds = bounds ma
>        a <- sequence (elems ma)
>        return $ P i (listArray bds a)

And now we can compose these funny dual-sided Kleisli/coKleisli arrows. But before writing a composition function, note what we have. We're looking at functions that can make a-ish things to b-ish things, and we can compose them like functions directly from a to b. This is precisely the design pattern captured by arrows. So we can write:

> data A m w a b = A { runA :: w a -> m b }

> instance (Distributes m w,Monad m,Comonad w) => Arrow (A m w) where
>    A g >>> A f = A $ bind f . distribute . cobind g
>    first (A f) = A $ \x -> do
>        u <- f (fmap fst x)
>        return $ (u,coreturn (fmap snd x))
>    pure f = A $ return . f . coreturn

Now check out page 1323 of Signals and Comonads by Uustalu and Vene. Honestly, I found that paper after I wrote the above code. It's almost line for line identical!

Anyway, now we can try an example. First some data to work on:

> x = P 0 $ listArray (0,4) (map return [0.0..4.0])

This is like last post's except that as I've now used a non-normalised blur operation I can show how arrow notation makes it easy to fix that:

> g = proc a -> do 
>    b <- A blur -< a
>    n <- A blur -< 1
>    returnA -< a-b/n

That's a very simple (and inefficient) high pass filter by the way. At the end of the day, however, we still need to be able to lift our arrows to act on x:

> liftCM :: (Distributes m w, Monad m, Comonad w) => A m w a b -> w (m a) -> w (m b)
> liftCM (A f) = cobind (\x -> distribute x >>= f)

> y = liftCM g x

And convert the result to something we can look at:

> result = sequence (toList (Main.array y))

It might not be quite the result you expected because we see a lot of I/O. Each computation that is tainted with I/O will carry that I/O with it. So every bit of I/O that was computed by any value that that went into the final result remains hanging round its neck until the bitter end, with duplications if the value was duplicated.

Of course there are lots of other monad/comonad pairs to try, so maybe there are some other interesting combinations lurking around.

Exercise: liftCM lifts an arrow A m w a b to a function w (m a) -> w (m b). You can also implement a lift to m (w a) -> m (w b). Implement it and figure out what it does.

Comonadic Arrays

On haskell-cafe, ajb, aka Pseudonym, laments that many people don't have enough experience with comonads to recognise them. So I thought I'd mention a really simple example of a comonad that nevertheless captures the essence of a large class of comonads. It's conceptually not much different to my cellular automaton example (making this a bit of a rerun), but this should be easier to understand. And if it's too trivial, I hint at a particle physics connection towards the end.

Firstly, you can skip this paragraph if you don't want a quick bit of theoretical discussion. Consider arrays of fixed dimension. As types, they look something like X^N for some fixed integer N. From a container we construct its zipper by applying X d/dX. In this case we get XNX^N-1=NX^N. In other words, the corresponding zipper is an array paired with an index into the array. We can stretch the meaning of comonad slightly to allow this to relate to arrays whose size isn't fixed.

So here's some code:

> import Data.Array

The usual definition of Comonad:

> class Functor w => Comonad w where
>  (=>>) :: w a -> (w a -> b) -> w b
>  coreturn :: w a -> a

And now a type that is a pair of an array and an index into that array:

> data Pointer i e = P i (Array i e) deriving Show

Think of it as an array with one of its elements singled out for special attention. It trivially inherits the functoriality of ordinary arrays:

> instance Ix i => Functor (Pointer i) where
>   fmap f (P i a) = P i (fmap f a)

And now comes the Comonad implementation. coreturn serves to pop out the special element from its context - in other words it gives you the special element, whle throwing away the array it lived in. (=>>), on the other hand, applies a function f of type P i a -> b to the entire array. The function is applied to each element in turn, making each element the special element for long enough to apply f.

> instance Ix i => Comonad (Pointer i) where
>   coreturn (P i a) = a!i
>   P i a =>> f = P i $ listArray bds (fmap (f . flip P a) (range bds))
>       where bds = bounds a

Compare with fmap for arrays. This walks through each element in turn, applying a function to each element, and return an array of results. The computation for each element is separate from all the others. With =>> however, the entire array may be used for the computation of each element of the result, with the index into the array serving to indicate which element it is we should be focussing on.

For example, here's an array of values:

> x = listArray (0,9) [0..9]

We want to consider this to be a circular array so that going off one end wraps around to the beginning:

> wrap i = if i<0 then i+10 else if i>9 then i-10 else i

Now here's a simple operation that 'blurs' the single ith pixel in the 1-D image represented by x:

> blur (P i a) = let
>       k = wrap (i-1)
>       j = wrap (i+1)
>   in 0.25*a!k + 0.5*a!i + 0.25*a!j

We can apply this to the entire image thusly

> test1 = P 0 x =>> blur

Note the curious way I have to use P 0 x as an input to blur. There seems to be a redundancy here, we want the resulting array and don't care what the focussed element is. But =>> wants us to give it a focal point. Curiously it's making explicit something that's familiar to C programmers, but is slightly hidden in C. In C, you refer to an array of floats using a float *. But the same type points to elements of the array as well. So when you point to an array, you are, in effect, blurring the distinction between a pointer to an array and a pointer to the first element. Comonads make that distinction explicit.

Anyway, suppose you wanted to apply a sequence of operations. Adding, blurring, scaling, nonlinearly transforming and so on. You could write a pipeline like

> x ==> f = f x
> test2 = P 0 x ==> fmap (+1) =>> blur ==> fmap (*2) ==> fmap (\x -> x*x)

Note how ==> fmap f ==> fmap g = ==> fmap (g . f). If you think of fmap as farming out workloads to a SIMD processor with one thread applied to each array element, sequences of fmaps correspond to threads that can continue to work independently. The comonadic operations, however, correspond to steps
where the threads must synchronise and talk to each other. I believe this last statement explains the a cryptic comment in the comments to this blog entry.

One final thing going back to my optional paragraph above. If D means d/dX then we can see the operator XD as a kind of number operator. When you apply it to an array like X^N the array becomes multiplied by a type corresponding to the array index type. For ordinary arrays, these are just integers. So you can see XD as a way of counting how many elements there are in a container. Also, for any container F, we also have the equations D(XF) = XDF+F, which we can write as DX=XD+1. At some point, when I have time, I'll point out how this is closely related to the Heisenberg uncertainty principle and how when we say that differentiation makes holes in a data type, it's related to the notion of a hole in solid state physics.

Oh, and I sketched this diagram but don't have time to write an explanation:

What does topology have to do with computability? Reloaded.

Now I’ve given the briefest of possible introduction to topology, I can get back to the issue of investigating the functions from () to ().

I introduced the idea of open sets as being the sets that can be observed by an idealised type of machine. So how about turning this notion onto itself. There are two possible outcomes from using one of our machines. It can respond “yes” or fail to respond at all. We can give these states the names ⊤ and ⊥ respectively. So the space of possible states is the set S = {⊤,⊥}, with the former meaning "yes". Now suppose we have a machine to watch machines. What are the possible observable subsets of S? Well we can easily build a machine that always responds ⊤ and we can also easily build a machine that always responds ⊥. So {} and {⊤,⊥} are observable. We can also build a machine that simply relays the state of the observed machine. That gives us the observable set {⊤}. But what about {⊥}?. Well there’s no way we can build a machine that responds with ⊤ when it sees ⊥ and with ⊥ when it sees ⊤. In order to respond ⊤ when it sees ⊥ it needs no know how long to wait before giving up. But the original definition of machine that I give says that we have no guarantee of how long a machine might take to respond. So no such machine is possible. In other words, the observable sets are {}, {⊤} and {{}, ⊤}, and this makes a topology on S. (Trivial exercise: check the axioms for this topology.)

S, equipped with this topology, is known as the Sierpinski space. There’s an interesting way to think about this. In Set Theory we have the notion of set membership, and membership is a binary valued thing. So the two element set {0,1} has a special role in that it gives a way to label the answers to the question “is x in the set y?:”. (In fact, this can be formalised nicely in category theory through the notion of a subobject classifier.) The Sierpinski space plays a similar role for sets when you only have a one-sided membership test. We can make this analogy precise:

Given a set X then any subset, U, of X, can be written as f^-1({1}) where f is a function f:X->{0,1}. Conversely, any such function gives rise to a subset.

Given a topological space X, then any open subset U, of X, can be written as f^-1({⊤}) where f is a continuous function f:X->S. Conversely, any such function gives rise to an open subset.

Now consider the continuous functions from S to S. As S has only two elements, there are only four possibilities to consider. (Almost trivial exercise: find the continuous functions. Remember, f is continuous if and only if f^-1(U) is open for all open U.) It turns out there are precisely three functions. They are essentially the same functions that I listed here for the purely lazy Haskell case. So at least for this simple case, the computable functions are precisely the continuous ones. And that is now my second answer to the question a month or two ago. By assigning the correct topology to your types you can arrange that the computable functions are precisely the continuous functions. At this point you can turn the tools of topology to computer science and you find that terms like discrete, Hausdorff and compact, which are all standard terms from topology defined inn terms of open sets, all have computational content. And of course this machinery generalises from the simple example S to any other type. (There is a proviso due to the fact that in topology, arbitrary unions of open sets must be open, but this doesn’t carry over directly without a little bit of tweaking.)

One researcher making use of such ideas is Martin Escardo. When I started writing these posts I was looking at a bunch of papers by Escardo that I found hard to decipher. Since starting, however, I discovered that this document is actually a fairly easy read. So it’s not worth me saying much more except to tell you to read that. And also check out the work of Paul Taylor.

I can now give a summary of what the last three posts have all been about:

It’s tempting to think of functions in a programming language like Haskell as being like ordinary mathematical functions and try to embed types and Haskell functions directly into Set Theory as sets and functions. After all, Haskell functions don’t have side-effects and so seem to behave much like mathematical functions. But not all functions terminate in a finite time, so they aren’t quite like ordinary functions after all. We could think of them in terms of partial functions, but partial functions don’t have a very rich structure and certainly don’t automatically rule out examples like the non-continuous function from () to (). But topology *does* provide a good framework. This idea goes back a few years, to the work of people like Smyth, Plotkin and Scott. But I’ve been programming (at least as a hobby) for 30 years, and I got my PhD in Mathematics about 16 years ago, and yet I never knew that these topics were so connected. So this work isn’t getting as much publicity as it should, even among people who know a little topology and computer science.

And I’ll finish by mentioning an application. An important property that topological spaces can have is known as compactness. I can give a rough idea of what this means quite easily: consider a finite row of pigeonholes and an infinite number of letters. There is no way of placing all of the letters in the pigeon holes without one of the pigeonholes ending up with an infinite number of letters. This isn’t the case if you have an infinite row of pigeonholes. On the other hand, consider the task of placing an infinite sequence of zero sized dots on a closed interval on the real line. You won’t be forced to put two dots on the same point, but there will be at least one point in the interval where the dots are infinitely bunched together. For example, if you place dots at 1/2,1/4,1/8,1/16,... inside the interval [-1,1] then in any small interval around zero there will be infinitely many dots. There'll be such a point whatever scheme you try. So even though [-1,1] isn’t finite, it’s the ‘next best thing’. One important thing about compact sets is that there are often strategies to exhaustively search them. And that is precisely what my earlier post about exhaustive searches was discussing. So topological ideas give a nice way to get a handle on these non-obvious search algorithms, and characterise those types for which such searches are possible.

PS In my previous post I introduced the idea of 'observability' as an intuition pump for topology. I think that idea goes back to the work of Stephen Vickers. Brits might recognise that name. Check out the wikipedia page on him if you don't.

What is Topology?

Of all the branches of mathematics I have studied, Topology was the one with the biggest mismatch between what I expected, and how it actually was. Eventually I reached Algebraic Topology, and that is all about the things you expect: Mobius strips, Klein bottles, pairs of pants, knots and handlebodies. But the subject starts with the machinery of Point Set Topology, and in particular with a set of axioms that can seem highly unmotivated. So I’m going to describe Topology in a way that is completely different from what you’ll find in any Topology textbook I know, but which does draw on published ideas in Computer Science. Nonetheless, I am describing standard off-the-shelf point set topology, just dressing it with different intuitions.

Suppose I have some set X and a (possibly infinite) collection of machines. To each machine M is associated a subset of X, U. The idea is that given an element x of X, you can use the machine M to test to see if x is an element of U. The catch is that the machine can only reply “yes”, and that you don’t know how long to wait for the answer, only that if the answer is “yes”, the machine will eventually, after a finite time, give you an answer. If x isn’t in U, the machine sits there forever, neither saying “yes” nor “no”. We’ll call a set V “observable” if for every x in V, we can use a machine, or some combination of machines, to show (in a finite time of course) that x is in V. So all of the U’s associated to machines M are observable.

Now suppose we have two machines, M and N, associated with sets U and V respectively. We input x to both machines. If x is in U intersection V then eventually both machines will reply “yes”. We don’t know how long to wait, but we know it’ll be a finite wait. Similarly, if x is in U union V we known that one machine or the other will eventually say “yes”. So U intersection V and U union V are both observable.

Now suppose we have a (possibly infinite) collection of observable sets U_i. Let U be the union of all the U_i. Given any element x of U, it must be a member of one of the U_i. But every U_i is observable, so there is some combination of machines that can prove that x is in U_i. And hence we can show that x is in U. So unions of infinite collections of observable sets are observable.

(Note that in the previous paragraph I didn’t say you had to have an algorithm for saying which U_i you should use, given x. I’m just saying that some U_i must exist. This means I’m not talking about semidecidability. But observability is similar to semidecidability.)

So what are the machines that I’m talking about? Let’s leave that undecided. We can use the above observations to make a definition:

Given a set X, a system of observables is a set of subsets of X, called observable sets, with the properties

1. X and the empty set are observable
   
2. arbitrary unions of observable sets are observable
   
3. finite intersections of observable sets are observable
   

So let me fill in with an example. Suppose we have some devices for measuring signed lengths (or, say, x-coordinates in some frame of reference). These lengths take value in the real line, R. If the device measures a length within its observable set it eventually says “yes”, otherwise it hangs until you press reset. Suppose x is in a device’s observable set. Then some mechanism eventually leads to a “yes” response. But if the device is based on reasonable physical principles, it must be possible to tweak x by a tiny amount and still trigger the same mechanism, otherwise the device would, in some sense, be infinitely discriminating around x. So if x is in an observable set U, then the interval [x-e,x+e] must also be in U for some, possibly very small, e. So let’s define a system of observables on the real line by saying U is observable if for every x in it, there is some e such that [x-e,x+e] is in the set. Intuitively, these are sets that are “fuzzy round the edges”. That’s a reasonable property of a real world measuring device.

Suppose we can augment our devices by applying transformations to a value x before applying the machine. For example, if the machine works optically you could imagine applying a magnifying lens to it to improve its accuracy. For a 10x lens, say, we’d apply the function x -> 10x to our point, and then apply a machine M to see if 10x lies in U. More generally, if our transforming function is f, then applying that to x converts our machine into a machine that tests whether x is an element of the set
f^-1(U) = {y | f(y) in U}

If f is the product of a physical process, it’s reasonable to expect f to be continuous at all of its arguments. So what does f^-1(U) look like for continuous f? Well, by definition, f is continuous at x if for all d there is an e such that |x-x’|<e means |f(x)-f(x’)|<d. So if x is in f^-1(U), then f(x) has a small interval around it in U (because by stipulation, U is observable) and so x has a small interval around it in f^-1(U) (by continuity of f). In other words, the continuous functions are precisely those functions for which f^-1(U) is observable for all observable U, ie. the functions that don't let us observe the unobservable.

So...I’ve talked about observability, and machines, and ways to pre-process the input to these machines. What does any of this have to do with topology? Well, it simply *is* topology. Not a generalisation, or a special case. What I have described is precisely the subject matter of topology. The only difference is that a “system of observables” is normally called a topology and an “observable” set is normally called an open set. And given a topology on a set A, and a topology on a set B, topologists define continuity by saying that f:A->B is continuous if f^-1 maps open sets to open sets.

What I’ve tried to do in a few paragraphs is motivate the usual axioms of topology. Usually they are presented as a fait accompli which is justified after the fact by geometrical intuitions. I hope I’ve given a way to look at these axioms that makes them seem natural a priori and that spans both geometrical and computational ideas.

Anyway, now we’re armed with the notion of a topology, I can get back to considering functions from () to ()...

How many functions are there from () to ()? (Updated)

I was recently looking at some of Paul Taylor's writings on what he calls Abstract Stone Duality. It's partly another approach to 'topologising' computer science, independently of Escardo's work. I'm not yet in a position to say much about what he does, except that it gives a neat new language to talk about computable continuous functions without having to build things up from Set Theory. But just for fun I want to look, out of context, at one teeny weeny little thing that he mentions, the Sierpinski space. If you've not met general topology before, it can be a tough topic to grasp. So what I want to do is look in detail at how we might think about a really simple seeming problem from a topological perspective.

Haskell defines a type called (). There is one value of type (), confusingly called (). Most of this post is about functions of the type () -> (). You'd imagine there couldn't be much to say. But this topic is much more complicated than some people might imagine.

How many functions of type () -> () are there?

(Before proceeding, the character ⊥ should look like _|_. Apologies if your font/browser makes it looks like something else.)

Let's start with the obvious function of this type.

> f1 x = x

It seems like this is the only possible function. The argument can only be () and the result can only be (). What other choices could we make?

Well here's another implementation:

> f2 _ = ()

It certainly looks different, but it still just maps () to (). So it appears that f1 and f2 are the same. But there is a way to tell them apart. Make the following definition:

> loop = loop

Attempting to evaluate loop sends Haskell into an infinite loop. Either you'll get an error message, or your computer will loop forever, never giving a result. Similarly, evaluating f1 loop will also fail to terminate with a sensible result. But f2 loop terminates fine and returns (). Because Haskell is a lazy language, f2 doesn't need to evaluate its argument before giving its result. f1, on the other hand, returns its argument, so looking at the result causes non-termination. So amazingly, f1 and f2 are different. We have at least two distinct functions of type () -> (). Are there any more?

We can treat a non-terminating function as if it returns a special value called ⊥ (pronounced bottom). So we can summarise the above in a kind of truth table:

That immediately suggests two more functions, f3 and f4, which are also illustrated.

So it now looks like there are four functions of type () -> (). Here's a possible implementation of f3:

f3 _ = loop

f3 is simply a function that deliberately sabotages itself. But f4 is the interesting one. In order to have f4 ⊥ = (), f4 must ignore its argument so that it doesn't get caught in a quagmire of non-termination. But if it ignores its argument, then it has to ignore its argument when f4 () is evaluated. So f4 () must also equal (). In other words, f4 cannot be implemented. We have only three functions of type () -> ().

So how can we characterise these functions mathematically? Haskell functions don't correspond to mathematical functions from the set {()} to {()}. If we think of () as a set with two elements, {(),⊥}, then Haskell functions still don't correspond to functions from this set to itself.

Let's give the set {⊥,()} a name, S. One approach to formalising Haskell functions is to impose an ordering on S. Here's one:

⊥<=⊥
⊥<=()
()<=()

(Now you can see why ⊥ is called "bottom". If you think of <= as a ranking, ⊥ is at the bottom of the heap.) A monotone function is defined to be a function f such that if x<=y then f(x)<=f(y). If we plot a 'graph' of f with ⊥ and () along the axes, monotone functions are the ones that don't descend as you go from left to right:
The possible Haskell functions correspond precisely to the monotone functions. (Mathematicians would probably call these functions monotonically increasing but the computer science literature I have just seems to call them monotone).

Another useful definition is the notion of a strict function. This is one that maps ⊥ to ⊥ and so goes through the 'origin' on the graph above. In a strict language you can't implement f2.

But maybe you can think of another possible implementation of a function of type () -> ():

> f5 = loop :: () -> ()

When working with fully lazy evaluation this is indistinguishable from f3 so we can ignore it. But in a strict language we can distinguish f3 and f5. In Haskell, for example, we can use the function seq to force evaluation of its first argument before moving onto its second. f3 `seq` () and f5 `seq` () are distinguishable because the first evaluates to () but the second gives ⊥.

So the answer to the original question is (I think):
* 1 in a total language or "in mathematics"
* 3 in a lazy language like Haskell when working completely lazily
* 4 in Haskell when using seq to enforce strictness
* 3 in a strict language like Ocaml


You'd think I'd have exhausted everything that there is to say by now. But there's a whole lot more. But before I get to it I need to talk a bit about topology in my next post.

And I'm completely ignoring the philosophical issue of what ⊥ means when you can't tell the difference between an answer of ⊥ and a slow computer taking longer than you expect to compute ().

Update: This problem was trickier than I originally anticipated. The above text incorporates a couple of changes based on comments I received.

A Third Order Quine in Three Languages

Suppose F_nxy is a program, written in language x, that takes as input n string arguments as input, G₁,...,G_n and outputs a program in language y that is the application of the function whose source code is G₂ to the strings G₂,G₃,...,G_n,G₁. Then F_3xy(F_3xy,F_3yz,F_3zx) will be a program in language x that that outputs a program in y that computes F_3yz(F_3yz,F_3zx,F_3xy).

We can put this into practice:

q a b c=putStrLn $ b ++ [toEnum 10,'q','('] ++ show b ++ [','] ++ show c ++ [','] ++ show a ++ [')']
main=q "q a b c=putStrLn $ b ++ [toEnum 10,'q','('] ++ show b ++ [','] ++ show c ++ [','] ++ show a ++ [')']" "def q(a,b,c):print b+chr(10)+'q('+repr(b)+','+repr(c)+','+repr(a)+')'" "def e(x) return 34.chr+x+34.chr end;def q(a,b,c) print b+10.chr+'main=q '+e(b)+' '+e(c)+' '+e(a)+' '+10.chr end"

This is a Haskell program that outputs a Python program that outputs a Ruby program that outputs the original Haskell program.

Apologies for the lack of line breaks. If it's not readable, it should be possible to copy and paste that source. Also, you may need to tweak it if your OS doesn't treat character 10 as a new line.

Purely functional recursive types in Haskell and Python

"""
This post is simultaneously Python and literate Haskell.

There is a certain truth to Greenspun's tenth law of programming. A Python project I was developing at work has slowly mutated into a compiler for a programming language without me planning it that way. Usually (I assume) compilers parse their input and construct an AST which is passed to the compiler proper. My code didn't have an AST, just a bunch of lambdas. I realised that I'd actually come across a real world example of what Wadler was talking about in Recursive Types for Free!.

In Haskell, the foldr function reduces a list using a binary function and some initial value. Suppose the function is called a and the initial value is b. Take a list, for example [1,2,3]. Now write it without using list notation, directly in terms of its constructors. Ie. 1:(2:(3:[])). foldr replaces (:) by a and [] by b. So this becomes a(1,a(2,a(3,b))). The best known example is a=(+) and b = 0 so we get 1+2+3+0 and hence the sum of the values in the list. Here is how we'd use foldr in Haskell:

> x = foldr (+) 0 [1,2,3]

The interesting thing is that anything you might want to know about a (finite) list can be extracted using foldr. There is a sense in which it the universal function on lists and all other functions can be factored through it. For example, we can implement head and tail as follows

> head = foldr const undefined
> tail x = let Just (_,t) = foldr tailHelper Nothing x in t where
>   tailHelper x Nothing      = Just (x,[])
>   tailHelper x (Just (y,z)) = Just (x,y:z)

So if x is a list, \a b -> foldr a b x tells you everything you could want to know about the list. In other words, you can completely replace the list itself with functions like this. In fact, we can replace the list constructors with functions that build such functions:

> nil a b = b
> cons h t a b = a h (t a b)

We can use nil and cons just like [] and (:). In fact, given an element defined by

> y = cons 1 (cons 2 (cons 3 nil))

We can convert it to a conventional list via

> z = y (:) []

So foldr embeds a list as a function.

We can write the same thing in Python. (Note that Python already has a variation of foldr, called reduce.)

"""

def foldr(a,b,l):
    if l==[]:
        return b
    else:
        return a(l[0],foldr(a,b,l[1:]))

print foldr(lambda x,y:x+y,0,[1,2,3])

"""

It's surprisingly easy to implement cons and nil in Python too:

"""

def nil(a,b):
    return b

def cons(h,t):
    def _cons(a,b):
        return a(h,t(a,b))
    return _cons

l = cons(1,cons(2,cons(3,nil)))

print l(lambda x,y:x+y,0)

print l(lambda x,y:[x]+y,[])

"""

Folds can be generalised to any recursive type, not just lists. (Stricly speaking I mean recursive rather than corecursive types. Folds aren't appropriate for infinite structures.) Note how for lists, foldr takes two arguments besides the list: a two argument function and a zero argument function. Applying a fold simply replaces the list constructors (:) and [] with these functions. Generalised folds do something similar: each constructor gives rise to an argument to the fold and when the fold is evaluated, each constructor is replaced with the appropriate function. Here's an example:

Now consider a simple expression type in Haskell:

> data Expr = X | Const Int | Binop (Int -> Int -> Int) Expr Expr

This is a recursive type so it has a generalised fold associated with it. This fold will take three arguments, one for each of X, Const and Binop, and each one will take the same number of arguments as the constructor. Here it is:

> efold :: a -> (Int -> a) -> ((Int -> Int -> Int) -> a -> a -> a) -> Expr -> a
> efold x _ _ X = x
> efold _ c _ (Const a) = c a
> efold x c b (Binop f lt rt) = b f (efold x c b lt) (efold x c b rt)

efold simply replaces each constructor with an application of the matching function recursively through the entire Expr.

Anything you might want to do to an Expr can be done using efold, and many things you might naturally want to do with an Expr are particularly easy to write using it. Here the functions to (1) evaluate the expression for X equal to some Int, and (2) to determine whether or not an expression is free of references to X:

> eval x e = efold x id id e
> freeX e = efold False (const True) (const (&&)) e
> identity e = efold X Const Binop e

Now we can do the same thing we did above, replace the Expr structure with its corresponding fold. And again, I'm implementing it in Python rather than Haskell:

"""

def X():
    def _X(x,c,b):
        return x
    return _X

def Const(a):
    def _Const(x,c,b):
        return c(a)
    return _Const

def Binop(f,l,r):
    def _Binop(x,c,b):
        return b(f,l(x,c,b),r(x,c,b))
    return _Binop

def eval(x,e):
    return e(x,lambda x:x,lambda f,l,r:f(l,r))

def freeX(e):
    return e(False,lambda x:True,lambda f,l,r:l and r)

"""

So we have translated the Haskell algebraic type Expr into functional expressions in Python. Here are some examples of their use:

Evaluating X, 2 and X+2 at X=2:

"""

print eval(3,X())
print eval(3,Const(2))
print eval(3,Binop(lambda x,y:x+y,X(),Const(2)))

"""

Testing whether 10-2 and X()+2 are free of references to X():

"""
print freeX(Binop(lambda x,y:x-y,Const(10),Const(2)))
print freeX(Binop(lambda x,y:x+y,X(),Const(2)))
"""

You can even implement a version in a blend of functional and OO style:

"""
class X:
    def __call__(self,x,c,b):
        return x

class Const:
    def __init__(self,a):
        self.a = a
    def __call__(self,x,c,b):
        return c(self.a)

class Binop:
    def __init__(self,f,l,r):
        self.f = f
        self.l = l
        self.r = r
    def __call__(self,x,c,b):
        return b(self.f,self.l(x,c,b),self.r(x,c,b))

"""

Some final comments:

This can sometimes be an inefficient style of programming, especially so in a strict language. Look again at tail for the cons/nil lists. But many uses are quite efficient, and folds capture a very common design pattern.

When I wrote this post a while back I left out mention of what the main point of the paper was. This post fixes that.

Wadler's paper also describes a dual version of this for codata such as streams. But as far as I understand it's not very interesting.

It's interesting that theory about static types has something to say about programming in a dynamically typed programming language.

Just so you know, my work project doesn't look anything like the code above.

Oh...and I guess you could say this was a form of the visitor pattern. Ugh. It's hideously complicated in C++.
"""

The Type that Should Not Be

> import Prelude hiding ((^))

I've been pondering a type that's so mind-bogglingly weird that I can no longer think about it without frying my neurons. So instead I'll just dump the contents of my brain here so you can fry your neurons too.

Here's an interesting type to warm up with. In Haskell:

> data X = Tree [X]

At first look it seems not be well-founded in the sense that a value of type X is a list of X's and we haven't yet defined what X's are. But we don't need an X to make an empty list of X's. So we can build elements inductively starting with Tree []. If you think of Tree [] as being a leaf, you can see why I called this type Tree (though they are trees whose limbs can have infinite sequences of children).

Even without the empty set forming a 'foundation' we can still build circular elements like

> x = Tree [x,x]

Now lists are ordered and can have repeated elements. What about if we eliminate these two properties. The way to do that would be with a set type. So how do we represent sets in Haskell. Well certainly not with Haskell's Data.Set because that only represents finite sets. On approach to supporting infinite sets is to identify a set with a predicate for telling whether an element is a member of a set. In other words, the type X→Bool makes a tolerably good set of X's.

So now we're ready for an interesting type:

> data U = U { f :: U -> Bool }

(BTW U stands for universe).

The first thing to think about is what this might mean mathematically. Well it looks like a solution to the equation X=2^X, in other words, a set that equals its own power set. We know that in set theory this gives rise to Russell's paradox. One way to embed recursive types into set theory is to tweak set theory, for example by using non-well-founded set theory. But Russell's paradox is so bad that it breaks just about any attempt at devising a set theory where a set can equal its power set.

And yet we can build this type in Haskell. How come? The explanation is simple, Russell's paradox tells us that there must be some elements whose evaluation doesn't terminate. In fact, we can replay Russell's paradox directly in Haskell. But before that we need some definitions.

Firstly, are we sure we can even construct an element of U? In order to define an element of U we need to define a boolean predicate on U but we can't even compare two elements of U for equality. But we can define the 'empty' element phi, along with a membership test:

> phi = U $ const False
> x `e` y = f y x

Unlike conventional set theory, we can define negation of elements. We can also define union and intersection:

> neg a = U $ \x -> not (x `e` a)
> a \/ b = U $ \x -> x `e` a || x `e` b
> a /\ b = U $ \x -> x `e` a && x `e` b

We're ready for the Russell paradox. Define the 'set' (I give in, let's informally call these things sets) of all sets that don't contain themselves.

> nonSelfContaining = U $ \a -> not (a `e` a)

Now try to compute nonSelfContaining `e` nonSelfContaining. It doesn't terminate. It's not just a matter of having implemented it badly, it can't possibly terminate as there is no answer that could make sense. So if you're designing a language that's intended to be total, you'd better rule out types like U.

So we now have two elements of U, phi and neg phi. Can we make any more? Well here's one way:

> c a = U $ \x -> a `e` x

Intuitively, c a is the set of all sets that contain a. But there are two potential problems that come to mind with this. One is, we've only made two sets so far, so "the set of all sets that contain a" might not be interesting. And secondly, can we be sure c a is distinct from a?

The first issue is easily dismissed. neg phi is obviously in c phi, but phi isn't. So we can distinguish phi and c phi. c phi is also different from neg phi (check to see if phi is an element of each). In fact, define functional power by

> f ^ 0 = \x -> x
> f ^ n = \x -> (f^(n-1)) (f x)

and it can be shown that c^m phi is distinct from c^n phi for m≠n. (Hint: Try checking to see if c^m (neg phi) is an element of c^n phi for various m and n.) So we have constructed an infinite number of elements of U.

But, we know that we can write code for elements of U that doesn't terminate. Why not restrict ourselves to just the computable elements of U? But there's a catch. How do we define a computable element of U when elements of U are functions. One approach is this: define x to be computable if y `e` x is computable for all y. But that won't do because it won't terminate if y isn't computable. So how about defining x to be computable if y `e` x is computable for all computable y. But that's circular, and not in a nice inductive way. U is so perverse that just attempting to define computability for it gives a circular definition. But that's not necessarily a problem. Define a 'universe' to be a set, C, of computations of type U such that C = the set of x such that y`e` x terminates for all y in C. The question now is this: is there a unique set C satisfying this 'equation'. We can rule out the empty set as phi must be in any such C. There's also another thing we can prove. If we have two such sets, C and C', we can prove that the intersection of C and C' is a universe. So we can take the intersection of all universes to form the smallest universe. So we can define a something in U to be computable if it is in this smallest set.

So here's my first open question: Are there larger 'universes' than the smallest one, or is it unique? It's not hard to see that a universe is closed under (non-recursive applications of) negation, finite union, finite intersection and the function c above. My second open question: What other computable elements are there?

Time to stop this and start cooking dinner.

What does topology have to do with computability?

I've been playing with code for exact real arithmetic. Unfortunately I've been too busy to get it polished but I do have time to mention two interesting and counterintuitive facts about real arithmetic on computers. It's also a good opportunity to say a tiny bit more about exhaustive search.
Firstly, what does exact real arithmetic mean? It's a widely held belief that computers (at least idealised Turing machines with as much memory as you need) can't represent real numbers precisely. This is a myth. For example, a common way to represent real numbers is as lists of binary digits. Any programming language that supports the notion of an infinite stream, like Haskell, can therefore represent an infinite stream of digits. If you're unhappy with the notion of an infinite stream, think functionally instead. Define a stream to be a function that takes as argument an integer n and returns the nth element. For example

third n | even n = 1
        | odd n  = 0

can be interpreted as exactly representing the real number 1/3=0.010101..₂. This is a perfectly well defined object that sits in a finite amount of memory. So we don't need infinite computers to represent real numbers as sequences of binary digits. (Though you may notice that we can't represent all real numbers, for example Chaitin's much overrated constant.) Exact real arithmetic is simply arithmetic with real numbers represented with some representation like this.

So now I'm ready to mention the first counterintuitive result. If we represent real numbers as streams of binary digits like above (maybe stored with an exponent or a separate integer part) then it's impossible to implement addition correctly.

The argument is straightforward. Suppose our addition function receives these two inputs: 0.000000... and 0.01111111... What should it output as the first digit after the "binary" point? If the first stream continues as repeating zeroes forever, and the second stream eventually turns into an infinite stream of zeroes, then this digit should be 0. But if both streams eventually turn into infinite streams of ones, then this digit should be 1, as a result of the carry that will happen. But there's no way an addition function can know what will eventually happen to these streams. In order to compute the first digit of the result, it needs to see all of the digits of the inputs. But no computable function can view all of the input digits. And hence addition is impossible to implement.

But that's no reason to give up. The problem isn't fundamental to exact real arithmetic, it's an issue with the representation being used. So what representation should we be using? One place to look is at how mathematicians define the real numbers. A popular approach is via the Cauchy sequence. A real number can be defined as a sequence of rational numbers a₁, a₂, a₃... with the property that for any ε, there is a point in the sequence beyond which the difference between any two elements is less than ε. For example the sequence 3,31/10,314/100,3141/1000,... might represent π. So we might try to copy this and represent real numbers as streams of rationals. In a way it works. The resulting type will have all the properties of real numbers that you know and love. But unfortunately it's not a very useful representation because you can't answer even the most basic queries about real numbers with it. For example, there's no way of extracting a rational approximation with given accuracy from such a sequence. You know that given any accuracy, if you walk down the sequence far enough you'll find a suitable rational, but you don't know how far down the sequence to walk. So we need to modify this representation so that not only is each term in the sequence an approximation to the desired real, but that there is a way of finding approximations to a desired accuracy. One simple way is this: choose the rational numbers in the sequence so that the nth rational is within 2^-n of our desired real. When I say 'Cauchy sequence' below I'll mean a sequence with a constraint like this.

Note that this representation is highly non-unique. There are many sequences that could be used to represent any real number. But this isn't a problem. Even the standard decimal representation is non-unique because 1 and 0.999... recurring both represent the same real. This means we can write functions on our real representation that depend on the specific sequence we've chosen to represent a real number and hence can't be interpreted as functions of the real numbers. (Eg. consider the trivial function that just returns the first element of he sequence.) But it's still possible to write functions are well-defined in the sense that if two different representations of the same real are passed in as inputs, then the returned sequences represent the same real, even though the actual sequences may be different. For example, we can write a function to halve a real number simply by halving all of the rationals in the Cauchy sequence. Exercise: give a well-defined implementation of addition using this representation. (Hint: It's a short one-liner in Haskell, but it's not completely trivial.)

Now I can state the second counterintuitive result: all well-defined computable functions are continuous.

Before proving that, I want to say something more about exhaustive search over streams. You may remember that I rewrote the key proof using more elementary language because the original work was phrased in the language of topology. Explaining that properly would require me to set up quite a bit of mathematical machinery. However, I can give a hint as to why topology is connected to the computability of certain functions. The key reason is this: by and large, Topology is the study of continuous functions.

Consider Cauchy's definition of continuity. The function f is continuous at x if for all ε, there exists a δ such that whenever |y-x|<δ, |f(y)-f(x)|<ε. There's another way of phrasing that. f is continuous at f if for any ε there is a δ such that even if the argument, x, to f has an error of δ, f(x) still has an accuracy of ε. So now think about exact reals defined by Cauchy sequences. Suppose f is a well-defined function, as above. Then to find f(x) to within 2^-m we just need to find the mth term in the Cauchy sequence f returns. If f is a computable function, it obviously can only examine a finite number of elements of the sequence of x. Suppose the latest element in the sequence it looks at is the nth one. Then if f is implemented correctly it's managing to compute f(x) to an accuracy of 2^-m even though its input was only known to an accuracy of 2^-n. In other words, f is continuous. (Choose m such that 2^-m<ε and δ<2^-n.)

We can take a slightly different perspective. For any x, x is a stream of rationals and f(x) is a stream of rationals. For f to be a proper computable function then to compute n terms of f(x) we need m terms of x, where m is some finite function of n. In other words, if we fire off a loop to compute the terms of f(x) one by one the code should keep producing terms one after another without ever been caught in a non-productive infinite loop. Now compare this with how I defined total functions for streams. So the condition that f be continuous is more or less the condition for f to continue to be productive. And that means we can talk about the computability of such functions in terms of continuity and use the language of topology to talk about computability. I find this quite amazing. Cauchy's centuries old definition of continuity contains within it the seeds you need to define the notion of computability for codata.

(One cool thing about this is that topology has a concept called compactness. Compact spaces share many properties with finite sets. In particular, the interval [0,1]={x|0≤x≤1} is compact. And that suggests that maybe it can be searched in a finite time. That suggestion is correct and is the reason why it's possible to exhaustively search the infinite space of streams in finite time.)

I wish I could give some references for further reading on this but I pulled it together from little clues scattered over many papers and presentations on the web. If someone knows the canonical reference on computability and topology, maybe they could post it in the comments.

Update: Don't miss this excellent post by Andrej Bauer.

What's all this E8 stuff about then? Part 3.

Isospin

The idea behind symmetry in physics is that there are things that you can change without making a difference to how a system evolves. I've already given the example of how rotating a pool table has no effect on the physics of the game. But here's another very different symmetry: in particle accelerator events that involve only the strong force, it makes no difference whether you're working with protons or neutrons. The strong force is unable to distinguish between these two types of particle. So we have a nice symmetry: for each proton or neutron in a system, we have a group of size two corresponding to the fact that we can swap one for the other. Now clearly protons and neutrons are different, they have different charges, but that only makes a difference to electromagnetic interactions, not strong interactions. This symmetry is called isospin symmetry.

But it gets more interesting when we consider quantum mechanics. Represent a proton by a vector (1,0) and a neutron by a vector (0,1). Quantum mechanics allows us to form superpositions by adding linear combinations of physical states to get new ones. So if (1,0) is a physical state, and (0,1) is also one, then so is (0.5,0.5). This is like the Schrödinger Cat state, it's a mixture of proton and neutron. In fact, we can form states (a,b) for any complex numbers a and b. So now we have scope for even more symmetry because not only should neutrons and protons be interchangeable, any mixture should be interchangeable with any other. What Heisenberg proposed was this: maybe the proton and neutron, considered as 2D vectors, have rotational symmetry just like the pool table. (Well, I don't think he was thinking in terms of pool tables.) Because we're dealing with complex numbers he proposed a symmetry group called SU(2), instead of SO(2), but it's a similar kind of thing, just the complex number variation of the 2D rotation group. SU(2) acts on the (complex) 2D space spanned by the proton vector (1,0) and the neutron state (0,1) and so this space forms what I previously called a Lie group representation. The proton and neutron form a 2-dimensional representation of SU(2).

So let's think about this more generally. If we have n particles that we suspect are connected by some symmetry group G, then they must form an n-dimensional representation of the group. Such a collection of particles is also called a multiplet. But here's the cool thing: if we have some kind of classification theorem listing the kinds of symmetry groups we expect to see in nature, and if we have a theorem classifying all of their representations, then we expect the particles of nature to be collected into families of similar particles that are organised just like those representations. So just knowing some pure mathematics should tell us a lot about particle physics.

Simple Lie Groups and E8

I've already said a little about how mathematicians go about the process of classifying representations of Lie groups by looking at maximal tori. Now I want to look at how Lie groups themselves are classified. Lie groups can often be built out of smaller pieces. This suggests the idea of finding the basic building blocks which can't be built out of simpler Lie groups. For a large class of groups, the basic building blocks are known as the simple Lie groups. (Now you know what the "Simple" in the title of Garrett Lisi's paper is about.) A priori you might imagine that you have all kinds of freedom when you try to construct new Lie groups. Incredibly you don't, and the list of simple Lie groups is quite...well...simple. The Lie groups fall into seven families called A, B, C, ..., G. The individual members of these families are given names like C4, where the integer is the dimension of the maximal torus. Four of the families, A to D, are fairly straightforward, and the groups form infinite sequences A1, A2, A3 and so on. Some of these are familiar. For example the B series are the rotation groups in odd dimensional spaces. B1 is just another name for SO(3). But curiously, the other three families, E, F and G, are finite. These are known as the exceptional Lie groups. The complete list is E6, E7, E8, F4 and G2. (What's even more exceptional about these is that many other branches of mathematics have exceptional objects too, and they're closely related to these. And now you know what 'exceptional' in that title means.) E8 is the highest dimensional of the exceptional simple Lie groups. In other words, you can think of it as the largest possible symmetry we could have that isn't simply built out of smaller pieces, or isn't just part of a straightforward sequence of symmetries. But it's not 'simple' in the ordinary sense of the word, it's 248-dimensional.

So if we could arrange the particles of nature into a family of size 248 then we'd have strong evidence that E8 was the symmetry of nature. This is what Lisi has attempted to do.

More Spontaneosly Broken Symmetry

But there's a complication. Above I talked out the SU(2) symmetry of protons and neutrons. But we know that symmetry doesn't hold exactly because the electromagnetic force doesn't respect it. In fact, photons react to protons and neutrons quite differently. But there's a way out: broken symmetry. In earlier posts I've already talked about this in the case of pool. But now consider a system like a sealed container of water. Again, fundamental principles suggest that this should have SO(3) symmetry. But we know that if we tilt a container of water the water doesn't tilt with the container, the level of the surface remains aligned perpendicular to gravity. On the other hand, as we heat up the water it eventually boils. At that point it's hard to tell gravity is at work at all because the force of gravity is so small compared to the intermolecular forces at work in a gas. The higher the temperature, the harder it is to see the local gravity. So pumping energy into a system serves to remove any kind of local accidents of history to reveal the true symmetry. The idea in particle physics is that at low energies we see a reduced symmetry, but at higher energies we should see more symmetries coming into play. In fact, this idea works incredibly well and in 1979 the Nobel prize was awarded to Glashow, Salam and Weinberg for realising that at higher energies the reduced symmetry of electromagnetism (which, by the way, is SO(2)) could be seen as part of a bigger symmetry. So what we expect to see is a chain of subgroups, one contained in the next, corresponding to physics at different energies. In practice physicists use Lie algebras instead of Lie groups most of the time because infinitesimal operations are usually much easier to work with, but the principle is the same.

So now I can finally say what Lisi has attempted to do: he has arranged all of the particles of nature to form a representation of E8. He has then found subgroups of E8 corresponding to the breaking of this symmetry. At low energies, the different forces of nature see different symmetries corresponding to these subgroups. As these subgroups correspond to less symmetry, we end up with smaller sets of particles that are interchangeable with each other, and so smaller multiplets. So the multiplets of E8 break up into smaller multiplets corresponding to representations of these subgroups. Lisi has shown how to arrange the particles of nature into multiplets, and has also shown how to break up these multiplets into smaller ones that correspond to the properties of the familiar fundamental forces we can observe in the lab. There are all kinds of problems with the way he's done this which makes the work controversial. But in spirit, it's not so different from what physicists have been doing for decades.

The IO Monad for People who Simply Don't Care

The theoretical underpinnings of monads in Haskell can be tricky to grasp. If you simply want to read some input and write some output, and you don't want to waste time with nonsense like the "monad laws", you've come to the right place.

Getting Started

Many programming languages make a distinction between expressions and commands. For example, in BASIC you couldn't say LET A = PRINT 1 and in C you can't say x = if (a==1) { 2; } else { 3; }. Haskell is no different. Like other languages it makes the distinction, and like other languages it has its own slightly idiosyncratic notion of what the difference is. The IO monad is just a device to help make this distinction. And that's the last time I'll mention the word 'monad' for a while.

Here's an example of a Haskell function:

f x = 2*x+1

The right hand side is an expression. Expressions are made of subexpressions and the value of an expression is determined by the value of its subexpressions. There is no room for anything like a print command here because a print command doesn't return a value, it produces output as a side-effect. So in order to do I/O we need something different: the do keyword. It's easy to use: you just write do and then follow it by an indented list of commands. Here's a complete Haskell program:

main = do
    print "Hello"
    let x = 1+2+3+4+5
    print "The answer is"
    print (2*x)

Note the let. It's just like the old days of BASIC programming where you can use let to assign a name to a value. Note also that commands take arguments that can be expressions. So print (2*x) is a command, but 2*x is an expression. Again, little different from a language like Python.

But what about input? Haskell has a command called getLine that reads input from the user. So here's an interesting feature of Haskell: commands can return values. But a command that returns a value is different from an expression with a value. The difference means we have to write code like this:

main = do
    print "Hello"
    x <- getLine
    print "The input was"
    print x

We have to use <- instead of let ... = ... to get a returned value out of a command. It's a pretty simple rule, the only hassle is you have to remember what's a command and what's a function.

You can also write your own commands. Here's a simple example that returns a pair of lines entered by the user:

get2Lines = do
    line1 <- getLine
    line2 <- getLine
    return (line1,line2)

main = do
    (a,b) <- get2Lines
    print ("You typed " ++ a ++ " then " ++ b ++ ".")

Note how we used <- instead of let to get the result of getLine. We also did the same to get the result of get2Lines. But notice one extra feature: we had to use return to return a value out of a do block.

So that's pretty much it. Here are the rules:

1. To introduce a list of commands, use do.
   
2. To return a value from a command, use return.
   
3. To assign a name to an expression inside a do block use let.
   
4. To assign a name to the result of a command, use <-.
   

So how do you know what's a command and what's an expression? If it has any chance of doing I/O it must be a command, otherwise it's probably an expression. But you can do silly things like:

f x = do
    return [x,x]

main = do
    x <- f 10
    print x

so occasionally there will be commands that don't do any IO. But not the other way round.

Types

Everything in Haskell has a type, even commands. In general, if a command returns a value of type a then the command itself is of type IO a. It just has an IO wrapped around the return type. For example, getLine is of type IO String and f above is of type a -> IO [a]. Apart from that extra IO it's like ordinary expression types. But note that if you have a bit of code like

main = do
    print (getLine)

you get a type error. This is because you must use <- to extract the input string from the getLine command. Because getLine has type IO String instead of type String the compiler can catch these kinds of mistakes. Giving commands and expressions different types means we can't accidentally mix them up.

At this point you might have realised that if

f x = do
    return x

main = do
    x <- f 5
    print x

is legal, then this might be legal too:

main = do
    x <- return 1
    print x

It is. return is simply a function of type a -> IO a that converts a value of type a to a command of type IO a.

You can also rewrite this

inputWithPrompt = do
    print "Hello, enter a line"
    x <- getLine
    return x

with

inputWithPrompt = do
    print "Hello, enter a line"
    getLine
    

In the first version we extracted x from getLine and then converted this value back into another command with return. In the second version we took a shortcut and simply ended with getLine. So we probably ought to add one more rule

5. The type of a sequence of commands is the type of the last line.

getLine already has the type of a command that returns a string so we can just make it the last line. But use the former version if you find this confusing, and then you only need four rules.

There's one last thing that confuses people. The type of an if statement is the type of its branches. So if you want an if statement inside a do block, it needs to be a command, and so its branches need to be commands also. So it's

main = do
    if 2/=1
        then do
            print "2/=1"
            print "Cool"
        else do
            print "2==1"
            print "Not cool"
    

not

main = do
    if 2/=1
        then
            print "2/=1"
            print "Cool"
        else
            print "2==1"
            print "Not cool"
    

Oh and a little tip. If the do block has just one command, you can drop the do. So this compiles:

main = do
    if 2/=1
        then
            print "2/=1"
        else
            print "2==1"
    

Why make the command/expression distinction at all?

You don't care, or you wouldn't be reading this. But suppose you did. Then I'd say: Firstly, let's state clearly what this distinction is. If something is of type IO a then it's a command returning an value of type a. Otherwise it's an expression. That's the rule. So here's the important bit: following the rules above it's completely impossible to put an executed command inside an expression. As the only way to do I/O is with commands, that means you have no way to find out what Haskell is doing inside expressions. And that's the whole point! Haskell likes to keep its evaluation secret from you and doesn't want you outputting information that might give away what's going on. For example, if (print 1)+(print 2) were legal Haskell, and it printed 1, followed by 2, then you'd know that it evaluated the left argument to + first, and then the second. Haskell doesn't want to reveal this because it wants the freedom to rearrange your code internally without changing its behaviour. And this makes things easier for you too. When you see a+b you know its exactly the same as b+a, whereas in C++, say, you have to worry about whether a and b have side effects. If the type isn't IO a, then you can sleep at night in the confidence that there are no side effects.

(Do you see why you can't sneak a command into an expression? It's quite neat. To make an expression out of a command you'd need to do something like extract a return value out of a command. But to do that you need to use <- to extract the value of the command. But you can only use <- inside a do block. So the result has to be a command. Once you're in a command, there's no way out!)

It's all lies I tell you!

One last thing. Much of what I've said above is false. But what I hope I have done is describe a subset of Haskell in which you can start some I/O programming without having a clue what a monad is. Once you're successfully bootstrapped, then you'll be ready for the truth. No more of that "monads are too hard" defense.

PS You can think of the above as using taint as a way to sequester commands from other expressions.

Back to E8 this weekend...

What's all this E8 stuff about then? Part 2.

Let's recap: I introduced the idea of symmetry groups, sets of operations that leave something unchanged. A particular type of group is the Lie group, it's a group where everything is smooth enough that you can do calculus with them. I also talked about Lie algebras which you can think of as elements of Lie groups that are infinitesimally close to doing nothing, or as rates of change of elements in a Lie group.

What does symmetry have to do with physics?

As my ultimate goal is talk about physics, it's time for a paragraph on that subject, just so we don't lose sight of where we're going. So imagine a game of pool. Pool is a game of physics, the outcome of any play depends on factors such as friction, collision angles, the coefficient of restitution of the balls (ie. their bounciness) and a whole host of other factors. But despite the complexity of the factors involved, we know that they are invariant under rotation. What I mean is that if we were to rotate the table and everything on it by some random angle (eg. from a North-South orientation to a NW-SE orientation) it makes no difference to the game. (Well, maybe the floor has a slight tilt, so let's rotate the room as well if that's the case.) If everything is rotated together then you end up with a situation that's essentially identical. One of the fundamental principles of physics is that the laws of physics don't care about the absolute orientation of anything. We called the set of rotations in a plane SO(2) so we can summarise this by saying that pool has SO(2) symmetry. But does it have SO(3) symmetry? Obviously rotating the table so that it tilts will change the way the game plays out. But that's because of an 'accident' of history, we just happen to be on the surface of the Earth and there's a gravitational field directed downwards. But from the point of view of fundamental physics there's nothing special about 'down'. If we were to tilt the pool table and then tilt the entire Earth underneath it, we should still expect to see no change to the physics of pool. In fact, we get to try this experiment out every day because the entire pool/Earth system rotates through (approximately) 360° every day and we don't notice any change in the way pool tables work. So pool actually has SO(3) symmetry, but because of stuff going on in our neighbourhood (ie. a having a planet under our feet), it looks like it only has SO(2) symmetry. SO(3) is a fundamental symmetry of physics, but that fact that on Earth gravity messes this up and leaves us with SO(2) symmetry is a phenomenon known as symmetry breaking. More of this later, back to the mathematics.

Lie Group Representations

I've talked about how groups can act on spaces. For example SO(3) acts on 3D space and SO(2) acts on 2D spaces, ie. planes. But there's some flexibility here and we can decouple the group from the space it acts on. The way SO(2) acts on a plane is given by a certain rule. We can make up new rules and study those. For example, we could allow SO(2) to act on 3D space by this rule: rotate your x and y coordinates using the same rule as in the plane, and leave z unchanged. So if a certain 2D rotation (ie. a 90 degree rotation) maps the point (x,y) to (y,-x), then according to this new rule it maps (x,y,z) to (y,-x,z). We could make up other rules. For example we could apply SO(2) to y and z, and leave x untouched. In fact, we can make infinitely many rules like this because there are infinitely many axes (by axis I just mean direction) we could choose to leave untouched. We can define a representation to be a rule that takes an element of a group and turns it into a transformation on a space. We're interested mainly in what are called linear representations (I'll give some justification for this in part 3). These are representations where the transformations map the origin to the origin, map straight lines to straight lines, and map the midpoints between pairs of points to the midpoints between pairs of points. It's not hard to see that rotations do all three of these things. In fact, as I will only talk about linear representations, I'll drop the word linear from now on.

I just showed how SO(2) has multiple representations because of the different ways in can be applied to 3D space. But we can also define alternative rules for how to apply SO(2) to the same space, eg. the plane. Here's a really simple one called the trivial representation: we simply say that elements of SO(2) do nothing. It's not very interesting, but it is a perfectly valid representation. Here's another: apply elements of SO(2) backwards. So if we have an element of SO(2) that says "rotate by 10° clockwise" the backwards rule says "apply a 10 degree rotation anticlockwise". It's a different representation. But if you look at the plane from underneath, it actually looks just like the original representation. So even though these are different representations, there's a sense in which it's equivalent to the original one.

But here's a representation of SO(2) that really is different to the original one: If the element of SO(2) says "rotate by x°" we rotate by 2x° instead. We can think of it as a double-speed rule. We apply the rotations in SO(2) twice. If we run through the elements of SO(2) starting at 0° and working our way up to 360°, then the representation rotates our space twice as fast and ends up rotating the space twice.

What about a half-speed rule? Sounds like it might work. But it fails for a simple reason. A 360 degree rotation is the same as a zero degree rotation. But the half-speed rule says that the former should rotate by 180° and the latter should rotate by 0°. These are distinct rotations, so our rule doesn't make any sense. As a result we're restricted to rules that are n-speed rules, where n is an integer. It's not hard to see that if we choose n to be a negative integer then by looking at the plane from underneath we get the same rule as using -n, looking from above. So we can discard the ones corresponding to negative integers. And we're left with one rule for each non-negative integer. In fact, these are all the representations that are possible on a 2D plane.

Weights

[Optional "advanced section": Now go back to 3D again. When we apply SO(2) to 3D space we find that SO(2) always leaves some axis fixed. So given a 3D representation of SO(2), our 3D space always splits into a 1D space that's left unchanged, and a 2D space. You can think of SO(2) as using the trivial rule on the 1D space, and using any of the rules in the previous paragraph on the 2D space. Something similar happens in any dimension of space. The representation of SO(2) will split up into pieces with some axes left untouched, and others, always grouped in pairs, that transform like in the previous paragraph. In other words, the representations of SO(2) can be broken down into fundmental building blocks which we call irreducible representations. In the previous paragraph I actually classified all of the irreducible representations of SO(2).

We can use the properties of SO(2) to understand other Lie groups. The idea is that not only is SO(2) frequently a subgroup of Lie groups, it's often a subgroup in multiple ways, ie. there are multiple ways to find a copy of SO(2) in SO(3). We've already seen that SO(2) can be embedded in SO(3) by interpreting rotations around any fixed axis as elements of a sub-SO(2). But notice how if we find two different embeddings of SO(2), they are forced to "interfere" with each other. For example, let's choose an SO(2) that acts on the x-y plane as described above. Suppose we now pick another SO(2) subgroup. No matter how we choose it, it must act on either x or y. There simply isn't room to find two SO(2)'s that don't at some point overlap with each other. But imagine the group of 4-dimensional rotations, caled SO(4). (It's not that scary, it's just like SO(3) except that we can make rotations that "mix up" any pair of directions, or combinations of such rotations.) We could pick one SO(2) that acts on the x-y plane and another that acts on the z-t plane (assuming the fourth axis is called 't'). But we won't be able to pick more than two for the same reason as before. Suppose we pick as many SO(2)'s as are possible. Then what we have is whats known as a maximal torus.

Here's a feeble attempt at drawing all of this. One of the SO(2)'s inside SO(4) rotates one pair of axes into each other, the other one rotates the other pair:

Now think about representations of a Lie group. These are a rule that tell you how each element in the group transforms our space. As the elements of our SO(2) subspaces obviously live in the group the rule must apply to these also. So a representation on a group also gives a representation of all of the SO(2)'s in it. So each SO(2) making up our maximal torus must basicaly act like we described above: by rotating some pairs of axes around at some "speed" and leaving other axes untouched. So if we pick a maximal torus of a Lie group, any representation splits up into a bunch of pairs of axes and each pair has an integer "speed" associated to each SO(2) in the maximal torus. The tuple of "speeds" for each pair of axes is known as a weight. By interpreting the weights as the coordinates of points, the collection of weights that arise from any particular representation can be drawn in a diagram. The dimension of the diagram is the number of non-interfering SO(2)'s in the maximal torus. And if you want to see what these diagrams look like, Garrett Lisi has some nice ones in his paper. Note that that paper also contains some "root diagrams". I don't have time to talk about those except to say that (1) they are weight diagrams for one particular special representation and (2) they tell you a lot about the geometry of all possible sets of weights for a particular Lie group.

One last thing for the "advanced" section: a similar analysis can be carried out for Lie algebras as opposed to Lie groups. When physicists draw diagrams of weights they are often talking about the weights of Lie algebras, but these things are intimately related.

There's probably one important message to take from this section: there are a lot of constraints on representations, you can't just make up any old rule. So just knowing that a given Lie group acts on some space you already know a lot of information, even if you don't know what the exact representation is. BTW One of the biggest applications of representation theory is in chemistry where you can read off information about the number of electrons allowed in atomic orbitals directly from representations of SO(3).]

Final words

That was pretty tough. Next time I'll talk about physics and it should get a bit easier. I'll explain why much of modern physics is the study of Lie group representations and I'l explain the 'exceptional' and 'simple' in the title of Garrett's paper.

(And I apologise for all of my sins of omission. For example, the analysis above only works for certain types of Lie group and really everything should be done using complex numbers, not real numbers. But I'm trying to compress a few hundred pages of mathematics into a single posting.)

What's all this E8 stuff about then? Part 1.

With all the buzz about E8 recently I want to give an introduction to it that's pitched slightly higher than pop science books, but that doesn't assume more than some basic knowledge of geometry and vectors and a vague idea of what calculus is. There's a lot of ground to cover and so I can't possibly tell the truth exactly as it is without some simplification. But I'm going to try to be guilty of sins of omission rather than telling outright lies.

Discrete symmetry groups.

Consider a perfect cube. If you turn it upside-down it still looks just like a cube in exactly the same configuration. Rotate it through 90 degrees about any axis going through the centre of a face and it looks the same. Rotate it through 120 degrees around a diagonal axis and it also looks the same. In fact there are 24 different rotations, including simply doing nothing, that you can perform on a cube, that leave it looking how you started. (Exercise: show it's 24.)



These 24 operations form the symmetries of the cube. Note how you can take any pair of symmetries and combine them to make another symmetry simply by applying one and then another. For example, if A represents rotation by 90° about the x-axis and B represents a rotation by 180° about the same axis then we can use shorthand to write AB to mean a rotation by 270°. The convention is that we read from right to left and so AB means "do B first, then A". If we define C=AB think about what AC means. It's a 270°+90° rotation, ie. a 360° rotation which does nothing. The convention is to call that 1. So C is the opposite or inverse of A because AC=CA=1. When you have a bunch of operations that you can combine like this, and every operation has an inverse, then the operations form what mathematicians call a group. The 24 symmetries of the cube form a group of size 24. Although you can think of a group as a bunch of operations on something, you can also think of a group as a thing in its own right, independent of the thing it acts on. For example, we can write down some rules about A, B and C above. AC=CA=1, C=AB, AAB=1, BB=1, C1=1C=C and so on. In principle we can give a single letter name name to each of the 24 symmetries and write out a complete set of rules about how to combine them. At that point we don't need the cube any more, we can just consider our 24 letters, with their rules, to be an object of interest in its own right.

We can think of the group as acting on the cube, but if we specified some suitable rules, we can make this group act on other things too. For example, we can use the coordinates (x,y,z) to specifiy any point in space. We used B to mean a 180 degree rotation around the x-axis. If we apply that rotation to any point in space (x,y,z) it gets mapped to (x,-y,-z). Similary A maps the point (x,y,z) to (x,z,-y). So the the 24 element cube group also acts on 3D space. This is an example of a group representation, and I'll get onto that later.

(So far I've been talking about rotations we can perform on the cube. But imagine you could also reflect the cube - as if looking at the cube in a mirror and then pulling the reflection back out into the real world. You can't do this with a real cube but you could with an imaginary one. How many symmetries do you think you get now?)

SO(3) and Lie groups

With the cube example above, all of our rotations had to be multiples of certain angles like 90° and 120°. Any kind of rotation through an angle like 37° leaves the cube in a weird orientation that doesn't look like the original. So the symmetries of the cube form a discrete group because there's only a discrete set of possibilities. But there are also continuous groups where you can smoothly slide from one symmetry to another and everything in-between is also a symmetry. For example consider a perfect 3D sphere at the origin. You can rotate it through any angle and it still looks like the same sphere. Any rotation about the origin is a symmetry of the sphere, and given two symmetries you can slide from one to the other. At every stage, the sphere will still look like the original sphere. Like in the discrete case, this set of symmetries forms a group, because you can combine and invert operations. But there are now infinitely many operations. We can rotate around any axis by any angle we choose and it's still a sphere. But even though there is an infinity of such rotations, it's still manageable. In fact, any rotation can be described by three parameters, for example the three that aeronautical engineers might use called pitch, roll and yaw. Mathematicians call the group of 3D rotations SO(3). That 3 refers to the fact that we're talking about rotations of 3-dimensional space. But because it takes three parameters to describe a rotation, SO(3) is itself 3-dimensional. This double appearance of the number 3 is a 'coincidence', as we'll see in the next section.

There's one more thing I need to add about continuous groups. Some are also differentiable. You don't need to worry in too much detail about what this means but it ultimately says that we can do calculus with these groups. For example, we can ask questions not only about the orientation of a sphere, but also about the rate of change of its orientation. So SO(3) is a differentiable group. The usual name for these is Lie groups, after the mathematician Sophus Lie.

Subgroups and SO(2)

SO(3) is the group of 3D rotations. But all of the 24 symmetries of the cube we considered are also rotations. So the cube is actually a subset of SO(3). And as it's a subset of a group that's also a group, we call it a subgroup.

Now consider the set of all rotations we could perform on a 2D space. These also form a group. In fact, given an element of this group we can apply it to some text to rotate it in 2D. As above, use 1 to mean "no rotation", use A to mean "rotate 45° anticlockwise" and B to mean "rotate 45° clockwise". We can illustrate how these act on Hello:


If you've ever playing with a drawing package, like the one I just used, this will all be very familiar. And it should also be obvious that 1=11111=AB=BA=AAAAAAAA=BBBBBBBB. (Note how I'm only talking about the orientation of the word. If we were taking into account the position of the word on the page it'd be more complex. But there would still be a group describing it.)

Although the screen the text Hello is on is flat and two-dimensional, it can be considered to be embedded in the 3D space around it. And all of the 2D rotations can be considered to be 3D rotations also, just 3D rotations that happen to be oriented along the axis pointing straight into the screen.



So 2D rotations can be considered to be particular 3D rotations and and SO(2) is a subgroup of SO(3).

Note that 2D rotations can be described by just a single number, the rotation angle. So SO(2) is 1-dimensional, even though it's defined through its effect on a 2D space. The two threes that appear in connection with SO(3) above really are a coincidence. In 4-dimensions we have SO(4) which is actually a 6-dimensional group. In 4D you need 6 numbers like pitch, roll and yaw to define your orientation. But I'm not asking you to try to visualise that now.

Rates of change and Lie algebras.

Imagine a sphere tumbling around. At any moment its orientation in space can be described by an element of SO(3). So we might naturally be led to ask the question "what is its rate of change of orientation?". What would such an answer look like?

Let's review basic calculus first. Consider position instead of orientation. Suppose a vehicle is driving along a road. At time t it's at distance A and at time t+1s it's at distance B. Then its average velocity over the 1s is (B-A)/1s. But that's an average over a period. What's its velocity exactly at time t? We could compute its distance, say C, at time t+0.0001s, and compute (C-A)/0.0001s. But that's still an average over a short period of time.The tool for telling us about instantaneous rates of change is calculus, and here it tells us that if the motion is continuous we can take a limit. So if at time t+d, its position is given by D, then the velocity is the limit, as d tends to zero, of (D-A)/d.

Back to the tumbling sphere. We could say "at time t it has orientation A and time t+1s it has orientation B". But as with velocity, that tells you what happened between t and t+1s, not what it was doing at any instant. We haven't said what the rate of change is precisely at time t. Imagine we take some kind of limit of the above statement as we consider the change from time t to time t+d where d becomes arbitrarily small. One thing we can say is that if the motion is continuous then the start and finish orientations, say A and D, will become arbitrarily close. The difference between them will be an infnitesimally small rotation. You can describe an infinitesimally small rotation by an infinitesimally small vector: let the direction of the vector be the axis of rotation and its length be the angle of rotation. Suppose this vector is V. Then the rate of change of orientation is V/d. This is commonly known as angular velocity and it's a vector. Note that it's not an infinitesimally small vector because the small V and small d will usually cancel to give a nice finite size angular velocity, just as the ordinary velocity calculation above gives a finite answer. In fact, we can easily interpret the angular velocity vector. Its direction is the axis of rotation and its length is the rate of change of angle around that axis (in degrees per second for example).

At this point you might say "why not represent any angle by a vector, not just an infinitesimal one". For example, you might think v=(180°,0,0) could be used to mean a rotation through 180° around the x-axis. You can do that if you like, but this doesn't resepct the additive properties of vectors. For example, it'd be useful if v+v=(360°,0,0) then represented "no rotation" because a 360° rotation is no rotation at all. But "no rotation" should be represented by (0,0,0). So this isn't a very convenient way to talk about rotations. But rates of change of orientation are different. If we define w=(180°/s,0,0) then w represents rotation around the x-axis at one revolution every two seconds. w+w=(360°/s,0,0) is also a perfectly good angular velocity, one that's twice as fast. Unlike angles, angular velocities don't wrap around when you get to 360°. You can crank angular velocities up as high as you like. So vectors are a good way to represent rates of change of orientation, even if they're not good at representing orientations.

Let's do this in 2D. It's much easier. All rotations are in the plane and so we don't need to specify an axis. We just need to specify a single number to define an orientation, and a rate of change of orientation is just given by a single number also. That one number can be thought of as a 1-dimensional vector.

So what kind of thing is this rate of change object? In each case it can be thought of as a vector (in the case of SO(2) it's a 1-dimensional vector), but it's not an element of the group because an angular velocity is not the same thing as an orientation. Mathematically it's called an element of a Lie algebra. The Lie algebra of a Lie group is the set of vectors that describe rates of change of elements of the group. I mentioned above that SO(4) is 6-dimensional. This means that the Lie algebra is made of 6-dimensional vectors. You can think of the 6 components of these vectors as describing how the 6 generalised pitch, roll and yaw parameters vary over time.

There's a convention that we use lower case names for Lie algebras. So the Lie algebra of SO(3) is called so(3). That makes a sort of intuitive sense, so(3) describes really small rotations and SO(3) describes full-sized ones.

Closing words

That brings me close to the end of Part 1. But I can now say one very important thing. There are a number of things that mathematicians call E8, but one of them is a Lie algebra. And that's the thing that Garrett Lisi is talking about.

PS I welcome feedback from anyone. I want to make this as understandable to as many people as possible.