A Neighborhood of Infinity

Memoizing Polymorphic Functions with High School Algebra and Quantifiers

2009-11-04T22:01:00.001-08:00

A little while back Conal Elliott asked about the memoization of polymorphic types. I thought it'd be fun to describe how to memoize such functions in the same spirit as Ralph Hinze's use of tries to memoize non-polymorphic functions. Along the way I'll try to give a brief introduction to quantified types in Haskell as well as showing some applications of the Yoneda lemma at work.

You can think of a generalized trie for a function type T as a type that's isomorphic to T but doesn't have an arrow '->' anywhere in its definition. It's something that contains all the same information as a function, but as a data structure rather than as a function. Hinze showed how to construct these by using the high school algebra axioms on non-polymorphic types. Polymorphic types are types involving quantification. So to make suitable tries for these we need to add some rules for handling quantifiers to high school algebra.

At first it seems unlikely that we could memoize polymorphic types. When Hinze demonstrated how to construct generalized tries he showed how to make a tree structure that was tailored to the specific types in hand. With polymorphic functions we don't know what types we'll be dealing with, so we need a one-size fits all type. That sounds impossible, but it's not.

The first thing we need to look at is universal quantification. Suppose F(a) is a type expression involving the symbol a. Then the type ∀a.F(a) can be thought of as being a bit like the product of F(a) for all possible values of a. So ∀a.F(a) is a bit like the imaginary infinite tuple


data Forall f a = (f Bool, f Int, f Char, f String, f [Bool], f (IO (Int -> Char)), ...)

One reason you can think of it this way is that all of the projections exist. So for any type you choose, say B, there is a function (∀a.F(a)) -> F(B). In Haskell the ∀ is written as forall and probably the best known example is the Haskell id function of type forall a. a -> a. For any concrete type B, id gives us a function of type B -> B. Note that we usually write the signature simply as a -> a. Haskell implicitly prepends a forall for every free variable in a type. We have to use the following line of code if we want to be able to use forall explicitly (among other things):


> {-# LANGUAGE RankNTypes, ExistentialQuantification, EmptyDataDecls #-}

I'll approach the tries through a series of propositions. So here's our first one:

Proposition 1
∀a. a = 0

0 is the type with no elements. ∀a. a is a type that can give us an object of type B for any B. There is no way to to this. How could such a function manufacture an element of B for any B with nothing to work from? It would have to work even for types that haven't been defined yet. (By the way, do you notice a similarity with the axiom of choice?) So ∀a. a is the type with no elements. Here's the usual way to write the type with no elements:


> data Void

We also have:
Proposition 2
∀a. a^a = 1

If we have a function of type forall a. a -> a then for any element of type a you give it, it can give you back an element of type a. There's only one way to do this - it must give you back what you gave it. It can't transform that element in any way because there is no uniform transformation you could write that works for all values of a. So ∀a. a^a has one element, id.

A slightly more interesting proposition is this:
Proposition 3
∀a. a^a.a = 2

A function of type (a,a) -> a gives you an a when you give it a pair of a's. As we don't know in advance what type a will be we can't write code that examines a in any way. So a function of this type must return one of the pair, and which one it returns can't depend on the value of the argument. So there are only two functions of this type, fst and snd.

We can rewrite the last proposition as ∀a. a^a² = 2. That suggests that maybe ∀a. a^aⁿ = n for any type n. We can go one better. Here's another proposition:

Proposition 4
For any functor F and type n, ∀a. F(a)^aⁿ = F(n)

I've already talked about that result. Here's an implementation of the isomorphisms:


> yoneda :: (forall b . (a -> b) -> f b) -> f a
> yoneda t = t id 

> yoneda' :: Functor f => f a -> (forall b . (a -> b) -> f b)
> yoneda' a f = fmap f a

Throughtout this article I'll use the convention that if f is an isomorphism, f' is its inverse.

Now it's time to look at a kind of dual of the above propositions. Instead of universal quantification we'll consider existential quantification. The type ∃a.F(a) is a kind of infinite sum of all types of the form F(a). We can imagine it being a bit like the following definition:


data Exist f a = ABool (f Bool) | AnInt (f Int) | AChar (f Char) | AString (f String) | AListOfBool (f [Bool]) ...

The important point is that given any element of any type we can turn it into an element of ∃a.F(a). You'd think that we could write this in Haskell as exists a. F(a) but unfortunately Haskell does things differently. The idea behind the notation is this: as we can put anything of type F(b) into it. So if X = ∃a.F(a) then we have a function F(a) -> X for any a. So we have a function of type ∀a. F(a) -> X. So although this type is existentially quantified, its constructor is universally quantified. We tell Haskell to make a type existentially quantified by telling it the constructor is universally quantified:


> data Exist f a = forall a. Exist (f a)

You can think of Exist as not being a single constructor, but an infinite family of constructors, like ABool, AnInt, etc. above.

If you have an element of an ordinary non-polymorphic algebraic sum type then the only way you can do stuff to it is to apply case analysis. To do something with an existential type means you have to perform a kind of infinite case analysis. So to do something with an element of ∃a. F(a) you need to provide an infinite family of functions, one for each possible type. In other words, you need to apply a function of type ∀a. F(a) →B to it.

Time for another proposition:

Proposition 5
∃a. a = 1

It seems weird at first that the sum of all types is 1. But once you put something into this type, you can no longer get any information about it back out again. If you try doing case analysis you have to provide a polymorphic function that accepts an argument of type ∀a. a, which is as good as saying you can't do any case analysis. Proposition 5 is actually a special case of the following:

Proposition 6
For any functor f, ∃a. (n^a, f(a)) = f(n)

Briefly, the reason this is that the only thing you can do with a matching pair of n^a and f(a) is apply the former to the latter using fmap. This is a kind of dual to the Yoneda lemma and I say more about it here.

We already know from high school algebra that this is true:

Proposition 7
x^y+z=x^y.x^z.

We can write the isomorphisms explicitly:


> prop7 :: (Either a b -> c) -> (a -> c, b -> c)
> prop7 f = (f . Left, f . Right)

> prop7' :: (a -> c, b -> c) -> Either a b -> c
> prop7' = uncurry either

It should be no surprise that the following 'infinite' version is true as well:

Proposition 8
x^{∃a. f(a)} = ∀a. x^f(a).

We can write the isomorphism directly:


> prop8 :: (Exist f a -> x) -> forall a. f a -> x
> prop8 g x = g (Exist x)
> prop8' :: (forall a. f a -> x) -> Exist f a -> x
> prop8' g (Exist x) = g x

We're now equipped to start constructing generalized tries for polymorphic functions. So let's consider memoizing the type forall a. [a] -> f a, for f a functor. At first this looks hard. We have to memoize a function that can take as argument a list of any type. How can we build a trie if we don't know anything in advance about the type of a? The solution is straightforward. We follow Hinze in applying a bit of high school algebra along with some of the propositions from above. By definition, L(a) = [a] is a solution to the equation L(a) = 1+a.L(a). So we want to simplify ∀a. f(a)^L(a). We have

f(a)^L(a) = f(a)^1+a.L(a) = f(a).f(a)^a.L(a) = f(a).f(a)^a+a².L(a) = f(a).f(a)^a.f(a)^a².L(a)

I hope you can see a bit of a pattern forming. Let's define T(n) = f(a)^aⁿ.L(a). Then

T(n) = f(a)^{aⁿ.(1+a.L(a))} = f(a)^aⁿ.T(n+1) = f(n).T(n+1)

That's it! We can translate this definition directly into Haskell.


> data T f n = T (f n) (T f (Maybe n))

I'm using the fact that Maybe n is standard Haskell for the type n+1. (But note that this equality is only valid when we think of the list type as data, not codata. So like with Hinze's original tries, values at infinite lists don't get memoized.)

To build the isomorphism we need to trace through the steps in the derivation. At one point we used aⁿ+a¹⁺ⁿ.L(a) = aⁿ.L(a) which we can implement as the pair of isomorphisms:


> lemma :: Either (n -> a) (Maybe n -> a, [a]) -> (n -> a, [a])
> lemma (Left f) = (f, [])
> lemma (Right (f, xs)) = (\n -> f (Just n),f Nothing : xs)

> lemma' :: (n -> a, [a]) -> Either (n -> a) (Maybe n -> a, [a])
> lemma' (f, []) = Left f
> lemma' (f, x:xs) = Right (maybe x f, xs)

We can put the other steps together with this to give:


> memoize :: Functor f => (forall a. (n -> a, [a]) -> f a) -> T f n
> memoize f = let x = prop7 (f . lemma)
>          in T (yoneda (fst x)) (memoize (snd x))

> memoize' :: Functor f => T f n -> forall a. (n -> a, [a]) -> f a
> memoize' (T a b) = let y = (yoneda' a, memoize' b)
>                 in prop7' y . lemma'

Let's try a few examples. I'll use the identity functor for the first example.


> data I a = I a deriving Show
> instance Functor I where
>    fmap f (I a) = I (f a)

Here's our first test function and some data to try it on:


> test1 (f, xs) = I $ if length xs>0 then head xs else f ()
> data1 = (const 1,[2,3,4])

In data1 we'e using a function to represent a kind of 'head' before the main list. For the next example we're leaving the first element of the pair undefined so that data2 is effectively of list type:


> test2 (f, xs) = reverse xs
> data2 = (undefined,[1..10])

We can test them by building the memoized versions of these functions.


> memo1 = memoize test1 :: T I ()
> memo2 = memoize test2 :: T [] Void

and then apply them


> ex1 = memoize' memo1 data1
> ex2 = memoize' memo2 data2

It appears to work!

So what's actually going on? We have

T(0) = f(0).T(1) = f(0).f(1).T(2) = ... = f(0).f(1).f(2).f(3)...

Now consider a function g :: [a] -> f a applied to a list. If the list isn't an infinite stream then it must have a certain length, say n. From these elements it builds something of type f a. However this f a is constructed, each of the elements of type a in it must be constructed from one of the n elements in the list. So if we apply g to the list [0,1,2,...,n-1] it will construct an element of f a where each a in it contains a label saying which position in the list it came from. (Compare with Neel's comment here). If we use integers we don't get a perfect trie because there are more elements of type f Integer than there are ways to indicate source positions. What we need is that for each length of list, n, we have a type with precisely n elements. And that's what the type n gives us.

We can memoize many different functions this way, though if the functor f is a function type you'll need to use some of Hinze's techniques to eliminate them. And you'll notice I haven't used all of the propositions above. I've tried to give some extra tools to allow people to memoize more types than just my two examples.

One last thing: I wouldn't use the type above to memoize in a real world application. But the methods above could be used to derive approximate tries that are efficient. One obvious example of an approximation would be to use Int instead of the finite types.

Update: I forgot to provide one very important link. This post was inspired by Thorsten Altenkirch's post here.

Buffon's Needle, the Easy Way

2009-10-31T08:35:00.000-07:00

Buffon's needle is a popular probability problem. Rule lines on the floor a distance d apart and toss a needle of length l<d onto it. What is the probability that the needle crosses a line? A solution is described at wikipedia but it involves a double integral and some trigonometry. Nowhere does it mention that there is a less familiar but much simpler proof, though if you follow the links you'll find it. In addition, the usual solution involves π but gives little intuition as to why π appears. The simpler proof reveals that it appears naturally as a ratio of the circumference of a circle to its diameter. I've known this problem since I was a kid and yet I hadn't seen the simpler proof until a friend sold me his copy of Introduction to Geometric Probability for $5 a few days ago.

So instead of solving Buffon's needle problem we'll solve what appears to be a harder problem: when thrown, what is the expectation of the number of times a rigid curved (in a plane) wire length l (no restriction on l) crosses one of our ruled lines d apart? Here's an example of one of these 'noodles'. It crosses the ruled lines three times:

Expectation is linear in the sense that E(A+B) = E(A)+E(B). So if we imagine the wire divided up into N very short segments of length l/N the expectation for the whole wire must be the sum of the expectations for all of the little pieces. If the wire is well behaved, for N large enough the segments are close to identical straight line segments. Here's a zoomed up view of a piece of our noodle:

For a small straight line segment the expectation must simply be a function of the length of the segment. The expectation for the whole wire is the expectation for one segment multiplied by the number of segments. In other words, the expectation is proportional to the length of the wire and we can write E(l)=kl for some constant k.

Now we know it's proportional to the length, we need to find the constant of proportionality, k. We need to 'calibrate' by thinking of a noodle shape where we know in advance exactly how many times it will cross the ruled lines. The following picture gives the solution:

A circle of diameter d will almost always cross the lines in two places. The length of this wire is πd so E(πd)=2 and k=2/πd.

The expected number of crossings for a wire of length l is 2l/πd. A needle of length l<d can intersect only zero or one times. So the expected value is in fact the probability of intersecting a line. The solution is 2l/πd.

No integrals needed.

The expected number of crossings is an example of an invariant measure, something I've talked about before. There are only a certain number of functions of a noodle that are additive and invariant under rotations and just knowing these facts is almost enough to pin down the solution.

Puzzle

Now I can leave you with a puzzle to solve. In the UK, a 50p coin is a 7 sided curvilinear polygon of constant width. Being constant width means a vending machine can consistently measure its width no matter how the coin is oriented in its plane. Can you use a variation of the argument above to compute the circumference of a 50p coin as a function of its width?

"What Category do Haskell Types and Functions Live In?"

2009-10-13T10:40:00.000-07:00

The question in my title is one that is often raised by Haskell programmers and it's a difficult one to answer rigorously and satisfyingly. But you may notice that I've put the question in quotes. This is because I'm not asking the question myself. Instead I want to argue that often there's a better question to ask.

Superficially Haskell looks a lot like category theory. We have types that look like objects and functions that look like arrows. Given two functions we can compose them just how arrows compose in a category. We also have things that look like products, coproducts, other kinds of limit including infinite ones, natural transformations, Kan extensions, monads and quite a bit of 2-categorical structure.

So what goes wrong? (Besides the obvious problem that on a real computer, composing two working functions might result in a non-working function because you run out of memory.)

Among other things, Haskell functions can fail to terminate because of things like infinite loops. Computer scientists often use the notation ⊥ to represent a non-terminating computation. So when we talk of the Haskell integers, say, we don't just mean the values 0, 1, 2, ... but we also have to include ⊥. Unfortunately, when we do this we break a few things. For one thing we no longer have coproducts. But people find it useful to talk about algebraic datatypes as constructing types using products and coproducts and that would be a hard thing to give up.

So we could restrict ourselves to considering only the category theory of computable functions. But that's not a trivial thing to do either, and it doesn't reflect what real Haskell programs do.

But even if we did manage to tweak this and that to get a bona fide category out of Haskell, all we'd get is a custom tailored category that serves just one purpose. One theme running through much of my blog is that Haskell can be used to gain an understanding of a nice chunk of elementary category theory in general. Showing that Haskell simply gives us one example of a category really isn't that interesting. When I talked about the Yoneda Lemma I felt like I was talking about more than just one property of some obscure category that I can't actually define and that most category theorists have never even heard of.

So what's going on? Why does it feel like Haskell is so naturally category theoretical while the details are so messy?

Going back to my Yoneda lemma code, consider my definition of check


> check a f = fmap f a

It's straightforward to translate this into standard category theoretical notation that applies to any category. Even though the code is implemented in a specific programming language there's nothing about it that prevents it being translated for use in any category. So it doesn't matter what category Haskell corresponds to. What matters is that this bit of code is written in language suitable for any category. And the proof I give can be similarly translated.

Consider this standard problem given to category theory students: prove that (A×B)×C is isomorphic to A×(B×C). In Haskell we could construct the isomorphism as:


> iso :: ((a,b),c) -> (a,(b,c))
> iso ((x,y),z) = (x,(y,z))

But now we hit a problem. We can straightforwardly translate this into mathematical notation and it will give a valid isomorphism in the category of sets, Set. But iso is written to accept arguments which are elements of some type. Not all objects in categories have elements, and arrows might not correspond to functions. And even if they did, if we were working with (certain types of) topological spaces we'd be giving a construction for the isomorphism, and our proof would show the underlying function had an inverse, but we'd be failing to show it's continuous. It looks like writing Haskell code like this only tells us about a particularly limited type of category.

But not so. Type cabal install pointfree to install pointfree and then run pointfree 'iso ((x,y),z) = (x,(y,z))' and it responds with


> iso = uncurry (uncurry ((. (,)) . (.) . (,)))

pointfree rewrites a function in point-free style. There are no x's, y's or z's in the written version, only uncurry, composition (.), and the product function (,). These exist in all Cartesian closed categories (CCC). So our original function definition, despite apparently referring to elements, can be mechanically turned into a definition valid for any CCC. We can now reinterpret the meaning of x, y and z in the original definition as not referring to elements at all, but as labels indicating how a bunch of fairly general categorically defined primitives are to be assembled together.

(Incidentally, my first foray into pure functional programming was to write a SASL compiler. It was little more than a bunch of rewrite rules to convert SASL code into point-free compositions of S, K and I, among other combinators.)

What we have here is an example of an internal language at work. I'm not sure what a precise definition of "internal language" is, but it's something like this: take a formal system and find a way to translate it to talk about categories in such a way that true propositions in one are turned into true propositions in the other. The formal system now becomes an internal language for those categories.

The best known example is topos theory. A topos is a category that has a bunch of properties that make it a bit like Set. We take a subset of the language of set theory that makes use of just these properties. Our propositions that look like set theory can now be mechanically translated into statements valid of all toposes. This means we can happily write lots of arguments referring to elements of objects in a topos and get correct results.

In their book Introduction to Higher-Order Categorical Logic, Lambek and Scott showed that "pure typed λ-calculus" is the internal language of CCCs. Even though expressions in the λ-calculus contain named variables these can always be eliminated and replaced by point-free forms. Theorems about typed λ-calculus are actually theorems about CCCs. When we write Haskell code with 'points' in it, we don't need to interpret these literally.

So despite not knowing which category Haskell lives in, much of the code I've written in these web pages talks about a wide variety of categories because Haskell is essentially a programming language based on an internal language (or a bunch of them). Despite the fact that even a function like iso might have quite complex semantics when run on a real computer, the uninterpreted programs themselves often represent completely rigorous, and quite general pieces of category theory.

So the question to ask isn't "what category does Haskell live in?" but "what class of category corresponds to the internal language in which I wrote this bit of code?". I partly answer this question for do-notation (a little internal language of its own) in an earlier post. Haskell (and various subsets and extensions) is essentially a way to give semantics to internal languages for various classes of category. However complicated and messy those semantics might get on a real world computer, the language itself is a thing of beauty and more general than might appear at first.

BTW This trick of reinterpreting what look like variables as something else was used by Roger Penrose in his abstract index notation. Just as we can sanitise variables by reinterpreting them as specification for plumbing in some category, Penrose reinterpreted what were originally indices into arrays of numbers as plumbing in another category. Actually, this isn't just an analogy. With a little reformatting abstract index notation is very close to the way I've been using monads to work with vector spaces so that abstract index notation can be viewed as a special case of an internal language for categories with monads.

Vectors, Invariance, and Math APIs

2009-10-10T09:12:00.000-07:00

Many software libraries, especially those for physics or 3D graphics, are equipped with tools for working with vectors. I'd like to point out how in these libraries the functions for manipulating vectors sometimes have special and useful properties that make it worthwhile to separate them out into their own interface.

Two types of Function

Suppose an object of mass m is moving with velocity v and we apply force f to it for time t. What is the final velocity? This is given by g:

g(m,t,f,v) = v+(t/m)*f

Now suppose that R is a rotation operation, typically represented by a matrix. What happens if we apply it to both of the vector arguments of g?

g(m,t,Rf,Rv) = Rv+(t/m)*Rf = Rg(m,t,f,v)

In other words, rotating the vector arguments is the same as rotating the vector result.

Another example: Consider the function that gives the force on an electric charge as a function of its velocity and the magnetic field:

f(e,v,B) = ev×B

It's essentially just the cross product. If you rotate both of the arguments to the cross product then the result is rotated too. The result is that

f(e,Rv,RB) = Rf(e,v,B)

On the other hand, many 3D APIs come with a function to perform componentwise multiplication of vectors. Write vectors x as triples (x₀,x₁,x₂), and so on, we can write such a function as:

f(x,y) = (x₀y₀,x₁y₁,x₂y₂)

You can show that this doesn't have a similar property.

Rotational Invariance

To make things easy, let's restrict ourselves to functions of scalars and vectors. And when I say vector, I'm talking strictly about vectors representing magnitude and direction, but not positions. Examples of such vectors are velocities, accelerations, angular velocities, magnetic fields, and the difference between two positions. A function is said to be rotationally invariant if applying a rotation R to all of its vector arguments results in the same thing as applying R to all of the vectors in its value. This allows you to have a function that returns multiple vectors, like a tuple or array.

The first two functions I described above were rotationally invariant but the third wasn't. Notice how the first two examples also described physical processes. This is the important point: as far as we know, all of the laws of physics are rotationally invariant. If you write down an equation describing a physical process then replacing all of the vectors in it by their rotated counterparts must also result in a valid equation. So if you're describing a physical process with a computer program, and you end up with a function that isn't rotationally invariant, you've made a mistake somewhere.

Vector APIs

Vector APIs frequently come with all manner of functions. Some have the invariance property and some don't. If you write code that you'd like to be rotationally invariant, but it turns out that it isn't, you usually have to examine the code to find the bug. But if you separate the invariant functions into their own interface, and then write code using just that interface, the code is guaranteed to be invariant. If your programming language has reasonably strict types then you may even be able to arrange things so that the type signature of the function alone is enough to tell you that the function is invariant. In effect you are able to make the compiler prove that your function is invariant.

(As an aside, this is an example of why a good type system does much more than you might at first have guessed. They don't just stop you making typos, they can do things like prove that your programs satisfy certain geometrical properties.)

So what functions would you have in such an API? Among the essential rotationally invariant functions are:

1. Multiplcation of a vector by a scalar
2. Addition of vectors
3. Dot product
4. Cross product

In terms of these you can build functions such as

1. Vector length
2. Vector normalization
3. Rotation of one vector around an axis specified by another vector
4. Linear interpolation between vectors

What kinds of functions would be excluded?

1. Constructing a vector from three scalars, ie. f(x,y,z) = (x,y,z).
2. Constructing a vector form a single scalar, ie. f(x) = (x,x,x).
3. Extracting the ith component of a vector, ie. f(i,(x₀,x₁,x₂)) = x_i.
4. Pointwise multiplication of vectors.
5. Computing the elementwise cosine of a vector.

On seeing the first excluded example above you might ask "how am I supposed to construct vectors?" The point is that you don't program exclusively with an invariant API, you simply use it whenever you need to prove invariance.

Coordinate System Invariance

One purpose of writing to a particular interface is that it allows the API to hide implementation details from the user. Using a rotationally invariant API has a role to serve here. For example, many 3D renderers allow you to write shaders. These are essentially functions that compute the colour of a piece of geometry that needs rendering. You write a shader and the renderer then calls your shader as needed when a fragment of geometry passes through its pipeline. Frequently these are used for lighting calculations but there are all kinds of other things that may be computed in shaders.

In a 3D renderer different parts of the computation are often performed in different coordinate systems. For example it may be convenient to perform lighting calculations in a coordinate system oriented with the direction of the light. But the author of a renderer doesn't want to be committed to a particular choice. In order to do this, it is essential to be able to write shaders that are agnostic about which coordinate system is being used. If we work with rotationally invariant functions, our shaders are guaranteed to be agnostic in this way (assuming that the only kind of coordinate change that takes place is a rotation).

Lots More Types

I've concentrated on just one type of invariance, rotational invariance. If we consider more types of invariance then more types of interface naturally emerge. It would take too long to cover all of the details here so I'm just going to briefly sketch the beginnings of the more general case. So just read this section as a list of pointers to further reading.

For example, some functions are invariant under translations. These can be thoght of as functions of points in space. If we allow more general linear transformations then we find that some common functions transform 'oppositely' to vectors. In particular, normals to surfaces transform in this way. In fact, Pixar's Renderman has three distinct types, vectors, points and normals that captures these different invariances.

If we go back to rotations again but now extend these by allowing reflections then we find an interesting new phenomenon. For example, consider the result of reflecting in the x-y-plane, followed by reflecting in the y-z-plane followed by reflecting in the x-z-plane. This simply multiplies vectors by -1. Dot product is invariant under this: (-x)·(-y)=x·y. But cross product isn't because (-x)×(-y)=x×y. Even though the cross product is apparently vector, it doesn't get multiplied by -1. When we start considering invariance under reflection we find that some vectors behave differently. These are the pseudovectors and in effect they have their own separate type and interface. Interestingly, nature likes to keep pseudovectors and vectors separate except in parity violating phenomena. There are even pseudoscalars.

Incidentally, if you consider invariance under scaling you're led to the idea of encoding dimensions in your types.

Conclusion

If you're writing a vector API think about the invariance properties that your functions may have. If any are useful then it may be worth placing those in a separate interface. The more distinct types you have, the more properties you can make your compiler prove. Obviously this needs to be balanced against practicality, complexity for users and what you actually need. To some extent, many existing APIs make some of these distinctions with varying degrees of strictness. The main point I want to make clear is that these distinctions are based on invariance properties, something that not all developers of such APIs are aware of.

At some point I hope to return to this topic and enumerate all of the common vector-like types in a single framework. Unfortunately it's a big topic and I've only been able to scratch the surface here. In particular there are some subtle interplays between dimensions and types.

On a deeper level, I think there must be some type theoretical framework in which these invariance properties are free theorems.

Update: I believe some of this material is covered in Jimm Blinn's Vectors and Geometry and Objects, Oh My!, but I don't have access to that. I suspect that there is one big difference in my presentation: I'm not so interested here in vectors (or normals or whatever) in themselves but as defining interfaces to functions with invariance properties. Like the way category theorists shift the focus from objects to arrows. It makes a difference because it immediately gives theorems that our code is guaranteed to satisfy. It's the invariance property of the cross product (say) that is useful here, not the fact that the components of a vector transform a certain way when we change coordinates (because I might not even want to refer to coordinates in my code).

Example Code

To show that I'm talking about something very simple, but still powerful, here's some Haskell code:


> data Triple = T Float Float Float deriving Show

> class Vector v where
>   (.+) :: v -> v -> v
>   (.*) :: Float -> v -> v
>   dot :: v -> v -> Float
>   cross :: v -> v -> v

> instance Vector Triple where
>   T x y z .+ T u v w = T (x+u) (y+v) (z+w)
>   a .* T x y z = T (a*x) (a*y) (a*z)
>   dot (T x y z) (T u v w) = x*u+y*v+z*w
>   cross (T x y z) (T u v w) = T
>       (y*w-v*z)
>       (z*u-x*w)
>       (x*v-y*u)

You can freely apply the four primitive functions to elements of type Triple but if you have a function of, say, signature


> f :: Vector v => (v,v,Float) -> [(v,v)]

you are guaranteed it is invariant.

test, ignore

2009-09-29T20:04:00.001-07:00


\int_{0}^{1}\frac{x^{4}\left(1-x\right)^{4}}{1+x^{2}}dx
=\frac{22}{7}-\pi

You know, I think I won't delete this. It's a celebration of my having embedded some TeX using the advice at Bot Cyborg.

Finite Differences of Types

2009-09-26T13:18:00.001-07:00

Finite Differences of Real-Valued Functions

Conor McBride's discovery that you can differentiate container types to get useful constructions like zippers has to be one of the most amazing things I've seen in computer science. But seeing the success of differentiation suggests the idea of taking a step back and looking at finite differences.

Forget about types for the moment and consider functions on a field R. Given a function f:R→R we can define Δf:R×R→R by

Δf(x,y) = (f(x)-f(y))/(x-y)

Δ is the finite difference operator. But does it make any kind of sense for types? At first it seems not because we can't define subtraction and division of types. Can we massage this definition into a form that uses only addition and multiplication?

First consider Δc where c is a constant function. Then Δc(x,y)=0.

Now consider the identity function i(x)=x. Then Δi(x,y)=1.

Δ is linear in the sense that if f and g are functions, Δ(f+g) = Δf+Δg.

Now consider the product of two functions, f and g.

Δ(fg)(x,y) = (f(x)g(x)-f(y)g(y))/(x-y)

= (f(x)g(x)-f(x)g(y)+f(x)g(y)-f(y)g(y))/(x-y)

= f(x)Δg(x,y)+g(y)Δf(x,y)

So now we have a Leibniz-like rule. We can compute finite differences of polynomials without using subtraction or division! What's more, we can use these formulae to difference algebraic functions defined implicitly by polynomials. For example consider f(x)=1/(1-x). We can rewrite this implicitly, using only addition and multiplication, as

f(x) = 1+x f(x)

Differencing both sides we get

Δf(x,y) = x Δf(x,y)+f(y)

That tells us that Δf(x,y) = f(x)f(y).

Finite Differences of Types

We're now ready to apply our operator to types. Instead of functions on the reals we work with functors on the set of types. A good first example container is the functor F(X)=X^N for an integer N. This is basically just an array of N elements of type X. We could apply the Leibniz rule repeatedly, but we expect to get the same result as if we'd worked over the reals. So setting f(x)=x^N we get

Δf(x,y) = (x^N-y^N)/(x-y) = x^N-1+x^N-2y+x^N-3y²+...+y^N-1

So we know that on types, ΔF(X,Y) = X^N-1+X^N-2Y+...+Y^N-1.

There's a straightforward interpretation we can give this. Differentiating a type makes a hole in it. Finite differencing makes a hole in it, but everything to the left of the hole is of one type and everything on the right is another. For example, for F(X)=X³, ΔF(X,Y)=X²+XY+Y² can be drawn as:

If you've been reading the right papers then at this point it should all become familiar. Finite differencing is none other than dissection, as described by Conor in his Jokers and Clowns paper. I don't know if he was aware that he was talking about finite differences - the paper itself talks about this being a type of derivative. It's sort of implicit when he writes the isomorphism:

right :: p j + (Δp c j , c) → (j , Δp c j ) + p c

With a little rearrangement this becomes the definition of finite difference.

Now that we've recognised dissection as finite difference we can reason informally about dissection using high school algebra. For example, we already know that lists, defined by L(X) = 1+X L(X) can be informally thought of as L(X)=1/(1-X). So using the example I gave above we see that ΔL(X,Y)=1/((1-X)(1-Y)) = L(X)L(Y). So the dissection of a list is a pair of lists, one for the left elements, and one for the right elements. Just what we'd expect.

Another example. Consider the trees defined by T(X)=X+T(X)². Informally we can interpret this as T(X)=(1+√(1-4X))/2. A little algebraic manipulation, using (√x-√y)(√x+√y) = x-y shows that

ΔT(X,Y) = 1/(1-(T(X)+T(Y))

In other words, a dissection of a tree is a list of trees, each of which is a tree of X or a tree of Y. This corresponds to the fact that if you dissect a tree at some element, and then follow the path from the root to the hole left behind, then all of the left branches (in blue) are trees of type X and all of the right branches (in red) are trees of type Y.

If you're geometrically inclined then you can think of types with holes in them as being a kind of tangent to the space of types. Along those lines, dissections become secants. I think this geometric analogy can be taken a lot further and that in fact a non-trivial piece of differential geometry can be made to work with types. But that's for another day.

Oh, I almost forgot. Derivatives are what you get when you compute finite differences for points really close to each other. So I hope you can see that Δf(x,x)=df/dx giving us holes in terms of dissections. Conor mentions this in his paper.

We should also be able to use this approach to compute finite differences in other algebraic structures that don't have subtraction or division.

I can leave you with some exercises:

1. What does finite differencing mean when applied to both ordinary and exponential generating functions?

2. Can you derive the "chain rule" for finite differences? This can be useful when you compute dissections of types defined by sets of mutually recursive definitions.

3. Why is right, defined above, a massaged version of the definition of finite difference? (Hint: define d=((f(x)-f(y))/(x-y). In this equation, eliminate the division by a suitable multiplication and eliminate the subtraction by a suitable addition. And remember that (,) is Haskell notation for the product of types.)

More Parsing With Best First Search

2009-09-13T17:27:00.000-07:00


> {-# LANGUAGE NoMonomorphismRestriction,GeneralizedNewtypeDeriving #-}

I have three goals in this post:

1. Refactoring the technique in my previous post so that building the search tree is entirely separate from searching the tree.
2. Making it work with real-valued weights, not just integers
3. Applying it to an ambiguous parsing problem, making use of a type class to define an abstract grammar.


> import Control.Arrow
> import Control.Monad
> import Control.Monad.Instances
> import Control.Monad.State
> import Data.Either
> import Data.Function
> import Random
> import qualified Data.List as L
> import qualified Data.Map as M

Search Trees

The idea is that I want to search a tree of possibilities where each edge of the tree is marked with a weight. The goal will be to search for leaves that minimise the sum of the weights of the edges down to the leaf.

Here's an example tree:

The minimum weight leaf is at C. If we're working with probabilities then we'll use minus the log of the probability of a branch as the weight. That way multiplication of probabilities becomes additions of weights, and the likeliest leaf has the minimum weight path.

So here's the definition of a search tree. I've given both leaves and edges weights:


> data Search c a = Leaf   { lb::c, leaf::a}
>                 | Choice { lb::c, choices::[Search c a] } deriving Show

(Compare with this.) lb is short for 'lower bound'. It provides a lower bound for the total weight of any option in this subtree (assuming non-negative weights). The tree in the diagram would look like:


> ex1 = Choice 0 [
>   Choice (-log 0.1) [
>       Leaf (-log 0.5) 'A',
>       Leaf (-log 0.5) 'B'],
>   Choice (-log 0.2) [
>       Leaf (-log 0.6) 'C',
>       Leaf (-log 0.4) 'D']]

This tree is a container in a straightforward way and so we can make it an instance of Functor:


> instance Functor (Search c) where
>    fmap f (Leaf   c a ) = Leaf c   $ f a
>    fmap f (Choice c as) = Choice c $ map (fmap f) as

But it's also a monad. >>= maps all of the elements of a tree to trees in their own right, and then grafts those trees into the parent tree:


> instance Num c => Monad (Search c) where
>    return = Leaf 0
>    a >>= f = join $ fmap f a where
>        join (Leaf   c a ) = Choice c [a]
>        join (Choice c as) = Choice c $ map join as

It's easy to make trees into a MonadPlus by simply grafting trees into a new root. MonadPlus is meant to be a monoid, but this operation, as written below, isn't precisely associative. But it's 'morally' associative in that two terms that are meant to be equal describe equivalent search trees. So I'm not going to lose any sleep over it:


> instance Num c => MonadPlus (Search c) where
>    mzero = Choice 0 []
>    a `mplus` b = Choice 0 [a,b]

For our searching we'll need a priority queue. I'll use a skew tree with code I lifted from somewhere I've forgotten:


> data Ord a => Tree a = Null | Fork a (Tree a) (Tree a) deriving Show

> isEmpty :: Ord a => Tree a -> Bool
> isEmpty Null = True
> isEmpty (Fork x a b) = False

> minElem :: Ord a => Tree a -> a
> minElem (Fork x a b) = x

> deleteMin :: Ord a => Tree a -> Tree a
> deleteMin (Fork x a b) = merge a b

> insert :: Ord a => a -> Tree a -> Tree a
> insert x a = merge (Fork x Null Null) a

> merge :: Ord a => Tree a -> Tree a -> Tree a
> merge a Null = a
> merge Null b = b
> merge a b
>  | minElem a <= minElem b = connect a b
>  | otherwise = connect b a

> connect (Fork x a b) c = Fork x b (merge a c)

At each stage in the search we'll pick the unexplored branch with the lowest total weight so far. So when we compare trees we'll compare on their lower bounds. So we need an ordering on the trees as follows:


> instance (Num c) => Eq (Search c a) where
>    (==) = (==) `on` lb

> instance (Num c,Ord c) => Ord (Search c a) where
>    compare = compare `on` lb

The real cost of a choice isn't just the weight immediately visible in the tree but the cost of the journey you took to get there. We use the bumpUp function to put that extra cost into the part of the tree we're currently looking at:


> bumpUp delta (Leaf   c a)  = Leaf   (delta+c) a
> bumpUp delta (Choice c as) = Choice (delta+c) as

The only tricky aspect to this code is that we need to be able to handle infinite trees. We can't have our code simply go off and return when it's found the next match because it might not be possible to do so in a finite time. Instead, the code needs to perform one operation at a time and report what it found at each point, even if that report is just stalling for time. We do this by returning a (possibly infinite) list containing elements that are either (1) the next item found or (2) a new update giving more information about the lower bound of the cost of any item that might be yet to come. This allows the caller to bail out of the search once the cost has passed a certain threshold.

(Returning a useless looking constructor to stall for time is a common design pattern in Haskell. It's an example of how programs that work with codata need to keep being productive and you get something similar with the -|Skip|- in Stream Fusion. First time I write the code I failed to do this and kept wondering why my infinite searches would just hang, despite my great efforts to make it as lazy as possible.)


> runSearch :: (Num c,Ord c) => Tree (Search c a) -> [Either c a]
> runSearch Null = []
> runSearch queue = let
>    m = minElem queue
>    queue' = deleteMin queue
>    in case m of
>        Leaf   c a  -> Left c : Right a : runSearch queue'
>        Choice c as -> Left c           : (runSearch $ foldl (flip insert) queue' $ map (bumpUp c) as)

A quick test of an infinite search: finding Pythagorean triples by brute force. We give each integer as cost one more than the previous one:

I guess this is actually Dijkstra's algorithm, but on a tree rather than a general graph.


> integers m = Choice 1 [Leaf 0 m,integers (m+1)]

> test = do
>    a <- integers 1
>    b <- integers 1
>    c <- integers 1
>    guard $ a*a+b*b==c*c
>    return (a,b,c)

> test1 = runSearch (insert test Null)

If you run test1 you'll notice how the output is noisy because of all those Left w terms. If you'e not worried about non-termination you could just throw out redundant output like so:


> reduce [] = []
> reduce (Left  a : Left b : bs) = reduce (Left b : bs)

Might as well convert weights to probabilities while we're at it:


> reduce (Left  a : bs) = Left (exp (-a)) : reduce bs
> reduce (Right a : bs) = Right a : reduce bs

This version should be a lot less chatty:


> test2 = reduce test1

Grammar

Now that searching works I can turn to an application - a more sophisticated example of what I briefly looked at previously), parsing with ambiguous grammars. So let me first build types to represent parsed sentences in a toy grammar:


> data Noun = Noun String deriving (Show,Eq,Ord)
> data Verb = Verb String deriving (Show,Eq,Ord)
> data Adj  = Adj  String deriving (Show,Eq,Ord)
> data Prep = Prep String deriving (Show,Eq,Ord)

The following two are noun phrase and prepositional phrase:


> data NP = NP [Adj] Noun deriving (Show,Eq,Ord)
> data PP = PP Prep  Noun deriving (Show,Eq,Ord)

And entire sentences:


> data Sentence = Sentence [NP] Verb [NP] [PP] deriving (Show,Eq,Ord)

We want to be able to print parsed sentences so here's a quick 'unparse' type class to recover the underlying string:


> class UnParse a where
>    unParse :: a -> String

> instance UnParse Noun where
>    unParse (Noun a) = a

> instance UnParse Verb where
>    unParse (Verb a) = a

> instance UnParse Adj where
>    unParse (Adj a) = a

> instance UnParse Prep where
>    unParse (Prep a) = a

> instance UnParse NP where
>    unParse (NP a b) = concatMap unParse a ++ unParse b

> instance UnParse PP where
>    unParse (PP a b) = unParse a ++ unParse b

> instance UnParse Sentence where
>    unParse (Sentence a b c d) = concatMap unParse a ++ unParse b ++ concatMap unParse c ++ concatMap unParse d

Now I'm going to approach the problem of parsing ambiguous sentences in two ways. One will be efficient, and one will be inefficient but represent the 'ground truth' against which we'll compare. (This reflects standard practice in graphics publications where authors compare their fancy new algorithm with an ultra-slow but reliable Monte Carlo ray-tracer.)

I'm going to assume that sentences in my language are described by a "context free" probability distribution so that a noun phrase, say, has a fixed probability of being made up of each possible combination of constituents regardless of the context in which it appears.

I need an English word for something that takes a grammar and does something with it but I'm at a loss to think of an example. I'll use 'transducer', even though I don't think that's right.

So a transducer is built from either terminal nodes of one character, or it's one of a choice of transducers, each with a given probability:


> class Transducer t where
>    char :: Char -> t Char
>    choose :: [(Float,t a)] -> t a

And here's our toy grammar. It's nothing like an actual natural language because real grammars take a long time to get right. Note I'm just giving the first couple of type signatures to show that the grammar uses only the Monad and Transducer interfaces:


> string :: (Monad t, Transducer t) => [Char] -> t [Char]
> string ""     = return ""
> string (c:cs) = do {char c; string cs; return (c:cs)}

So, for example, a noun has a 50% chance of being the string ab and a 50% chance of being the string ba:


> noun :: (Monad t, Transducer t) => t Noun
> noun = do
>    a <- choose [(0.5,string "ab"),(0.5,string "ba")]
>    return $ Noun a

> verb :: (Monad t, Transducer t) => t Verb
> verb = do
>    a <- choose [(0.5,string "aa"),(0.5,string "b")]
>    return $ Verb a

> adjective :: (Monad t, Transducer t) => t Adj
> adjective = do
>    a <- choose [(0.5,string "ab"),(0.5,string "aa")]
>    return $ Adj a

> parsePrep = do
>    a <- choose [(0.5,string "a"),(0.5,string "b")]
>    return $ Prep a

Some of our "parts of speech" allow sequences of terms. We need some kind of probabilistic model of how many such terms we can expect. I'm going to assume the probability falls off exponentially with the number of items:


> many :: (Monad t, Transducer t) => Float -> t a -> t [a]
> many p t = choose [
>   (p,return []),
>   (1-p,do
>     a <- t
>     as <- many p t
>     return $ a:as)]

I also have a convenience function for sequences of length at least 1:


> many1 p t = do
>     a <- t
>     as <- many p t
>     return (a:as)

And now the rest of the grammar:


> parseNP = do
>    a <- many 0.5 adjective
>    b <- noun
>    return $ NP a b

> parsePP = do
>    a <- parsePrep
>    b <- noun
>    return $ PP a b

> sentence = do
>    a <- many 0.5 parseNP
>    b <- verb
>    c <- many 0.5 parseNP
>    d <- many 0.5 parsePP
>    return $ Sentence a b c d

We're going to use this grammar with two instances of type Transducer. The first will use the rules of the grammar as production rules to generate random sentences. The second will parse strings using the grammar. So we get two uses from one 'transducer'. This is pretty powerful: we have described the grammar in an abstract way that doesn't asuume any particular use for it.


> newtype Generator a = Generator { unGen :: State StdGen a } deriving Monad
> newtype Parser a    = Parser { runParse :: (String -> Search Float (a,String)) }

Let's implement the generation first:


> instance Transducer Generator where
>    char a = return a
>    choose p = do
>        r <- Generator (State random)
>        case (L.find ((>=r) . fst) $ zip (scanl1 (+) (map fst p)) (map snd p)) of
>            Just opt -> snd opt

We can test it by generating a bunch of random sentences:


> gen = mkStdGen 12343210
> generate n partOfSpeech = (unGen $ sequence (replicate n partOfSpeech)) `evalState` gen

> test3 = mapM_ print $ generate 10 sentence

We can now use generate-and-test to estimate what proportion of randomly generated sentences match a given sentence:


> generateAndTest n partOfSpeech chars = do
>    a <- generate n sentence
>    guard $ unParse a == chars
>    return a

> collectResults n partOfSpeech chars = M.fromListWith (+) $ map (flip (,) 1) $
>   generateAndTest n partOfSpeech chars
> countResults n partOfSpeech chars = mapM_ print $ L.sortBy (flip compare `on` snd) $
>   M.toList $ collectResults n partOfSpeech chars

> test4 = countResults 100000 (noun :: Parser Noun) "abab"

On the other hand we can build a parser, based on Hutton's, just like in my previous post except using this new tree search monad:


> instance Monad Parser where
>   return a = Parser (\cs -> return (a,cs))
>   p >>= f = Parser (\cs -> do
>       (a,cs') <- runParse p cs
>       runParse (f a) cs')

> instance MonadPlus Parser where
>   mzero = Parser (\cs -> mzero)
>   p `mplus` q = Parser (\cs -> runParse p cs `mplus` runParse q cs)

> instance Transducer Parser where
>    char c = Parser $ char' where
>       char' "" = mzero
>       char' (a:as) = if a==c then return (a,as) else mzero
>    choose p = foldl1 mplus $ map (\(p,x) -> prob p >> x) p where
>       prob p = Parser (\cs -> Leaf (-log p) ((),cs))

> goParse (Parser f) x = runSearch $ insert (f x) Null

> end = Parser (\cs -> if cs=="" then return ((),"") else mzero)

> withEnd g = do
>   a <- g
>   end
>   return a

> normalise results = let total = last (lefts results)
>   in map (\x -> case x of
>       Left a -> a / total
>       Right b -> b
>   ) results

> findParse chars = mapM_ print $ reduce $ runSearch $
>   insert (runParse (withEnd sentence) chars) Null

Results

And now we can try running both methods on the same string:


> main = do
>    let string = "ababbbab"
>    findParse string
>    print "-------------------"
>    countResults 1000000 (sentence :: Parser Sentence) string

You should see the parsings from countResults in roughly the same proportion as the relative probabilities given by findParse. Remember that the relative probability of a given parsing is the last Left p term before that parsing. Try playing with string, the number of Monte Carlo runs and the seed. Remember that there is going to be some variation in the randomised algorithm, especially with hard to parse strings, but raising the number of runs will eventually give reasonable numbers. Of course ultimately we don't care about the Monte Carlo method so it's allowed to be slow.

Anyway, none of this is a new algorithm. You can find similar things in papers such as Probabilistic tree transducers and A Generalization of Dijkstra's Algorithm. But what is cool is how easily Haskell allows us to decouple the tree building part from the searching part. (And of course the tree is never fully built, it's built and destroyed lazily as needed.) All of the published algorithms have the parsing and searching hopelessly interleaved so it's hard to see what exactly is going on. Here the search algorithm doesn't need to know anything about grammars, or even that it is searching for parsings.

Semiring Parsing is also easy to implement this way.

BTW If you think my "ab" language is a bit to contrived, check out the last picture here for an example of some natural language that is in a similar spirit :-)

A Monad for Combinatorial Search with Heuristics

2009-07-04T17:35:00.000-07:00

Haskell provides a great way to perform combinatorial searching with backtracking: the list monad. Do-notation provides a nice DSL that makes it easy to express the trying out of different possibilities. But the list monad only performs a simple-minded walk through all of the alternatives giving little opportunity to direct that walk. In particular, it's not easy to provide heuristics to say things like "try this alternative first but if it starts going badly consider this alternative too". This post contains a monad that gives a simple scheme to allow programmers to direct searches in this way.

First the Haskell administrativia...


> import Data.Char
> import Control.Monad
> import Data.Monoid
> import Data.List

When using the list monad, a list is interpreted as a list of candidates in a search. The join function for this monad takes a list of lists of candidates and flattens it into a list of candidates. This is all the list monad really does: you write code that generates new candidates from old, and the >>= function applies this code to all of the candidates it knows about and then flattens this back out to a single list of candidates. Importantly it does this in a lazy way so that you only need look at candidates as they are generated.

This new monad will keep slightly more information: each candidate will have a 'penalty' value attached to it saying how attractive a candidate it is. Candidates with score 0 will be tried first, and those with score n will be tried after those with lower scores. We can represent a collection of candidates and their scores simply as a list of lists. The first list in the list will have those with score 0, the second will have those with score 1 and so on. We'll call these lists penalty lists and the positions within those lists slots.

Here's the definiton of the penalty list type:


> data PList a = P { unO :: [[a]] } deriving (Show,Eq)

It's a functor in a straightforward way:


> instance Functor PList where
>   fmap f (P xs) = P (fmap (fmap f) xs)

The rule we'll adopt is that if you're trying a combination of two candidates then the penalty associated with the combination is the sum of the penalties of the individual objects. To implement this we need an alternative version of the join operation. If we have a penalty list of penalty lists and we have an element in the mth slot in the nth penalty sublist then we want it to end up in the (m+n)th slot in the final penalty list. Within a slot we can just order the elements just like in the original list monad.


> headm :: Monoid m => [m] -> m
> headm (a:as) = a
> headm [] = mempty

> tailm :: Monoid m => [m] -> [m]
> tailm (a:as) = as
> tailm [] = []

> zipm :: Monoid m => [[m]] -> [m]
> zipm ms | all null ms = []
> zipm ms = let
>   heads = map headm ms
>   tails = map tailm ms
>   h = mconcat heads
>   t = zipm (filter (not . null) tails)
>   in h : t

> instance Monad PList where
>   return x = P [[x]]
>   x >>= f = let P xs = (fmap (unO . f) x) in P (join xs) where
>       join []     = []
>       join (m:ms) = let
>           part1 = zipm m
>           part2 = join ms
>           in headm part1 : zipm [tailm part1,part2]

Explaining how join is implemented would take many words so I hope this picture of the computation of an example will do instead. I used the Monoid class simply to avoid directly referring to one level of nesting of brackets. It is intended to be a proper implementation of a Monad satisfying the three monad laws but I haven't proved this and it's possible that it occasionally leaves trailing empty lists around - which have no impact on search results.


> instance MonadPlus PList where
>    mzero = P []
>    mplus (P xs) (P ys) = P (zipm [xs,ys])

We can use this much like the list monad. First it will search for possibilities with zero penalty. When these are exhausted it'll backtrack to the last place where it can start finding possibilities with penalty 1. Then it'll try penalty 2 and so on. Importantly it manages to do this lazily so that we don't explore penalty n+1 until we've finished penalty n.

So now we can start using it. We'll hunt for Pythagorean triples by simply hunting through all of the triples of integers. But we'll try to find solutions where the sum of the integers is as small as possible. So as list of candidate integers we use P [[1],[2],[3]...]. In other words, the integer n has penalty n-1. Here's the code:


> ex1 = do
>   x <- P $ map (\x -> [x]) [1..]
>   y <- P $ map (\y -> [y]) [1..]
>   z <- P $ map (\z -> [z]) [1..]
>   guard $ x*x+y*y==z*z
>   return $ (x,y,z)

Of course we wouldn't really search for Pythagorean triples this way. This is just an illustration of how to use the code. But note, crucially, that the equivalent code using the regular list monad would give us back no solutions. It'd start with x=1 and y=1 and then go off to infinity finding candidates for z. So as a side effect the penalty list allows us to tame some infinite searches.

Anyway, that was a simple numerical example. But this monad can be used with much more complex kinds of search. In fact it almost serves as a drop-in replacement for the list monad. This is a really nice example of the way separation of concerns is easy in Haskell. The task of generating candidates for search can easily be separated from the task of selecting from those candidates, even though the operations are highly interleaved during execution.

So here's a more complex example: writing a parser that can tolerate errors without running into combinatorial explosion. The idea is that we associate a penalty with each error. The penalty will make the parser run on the assumption of no errors until it can no longer parse, and then it'll backtrack on the assumption of one error until that assumption is no longer tenable and so on. We can liberally sprinkle 'erroneous' parsings throughout our code confident that these branches will only be taken in the event that an error-free parsing can't be found.

Firstly, here's a penalty list that we can use to introduce a penalty of just 1.


> penalty :: PList ()
> penalty = P [[],[()]]

If we stick that in the code path then anything following acquires a penalty of 1.

Now we can write a parser. We can implement Hutton's parser in his monad parsers paper with very little modification. We simply replace the usual list with the penalty list and do away with the +++ operator to allow it to be a bit more liberal about backtracking. Here's the parser type:


> newtype Parser a = Parser (String -> PList (a,String))

We could have parameterised that with the underlying monad so that we could have parsers with a choice of search strategy.

The rest is a lot like in Hutton's paper:


> parse (Parser f) x = f x

> instance Monad Parser where
>    return a = Parser (\cs -> P [[(a,cs)]])
>    p >>= f = Parser (\cs -> do
>        (a,cs') <- parse p cs
>        parse (f a) cs')

> instance MonadPlus Parser where
>    mzero = Parser (\cs -> mzero)
>    p `mplus` q = Parser (\cs -> parse p cs `mplus` parse q cs)

> item :: Parser Char
> item = Parser (\cs -> case cs of
>    "" -> mzero
>    (c:cs) -> P [[(c,cs)]])

> sat :: (Char -> Bool) -> Parser Char
> sat p = do
>    c <- item
>    if p c then return c else mzero

> char :: Char -> Parser Char
> char c = sat (c ==)

Now for a simple parsing problem. We'll parse simple arithmetical expressions a lot like in Hutton's paper. But I'm going to tolerate two kinds of error:
1. The shift key doesn't always work so occasionally a shifted or unshifted version of a character may appear and
2. parentheses are occasionally left out by the clumsy user.

Now we can code up a simple grammar for this. First the mapping between shifted and unshifted characters (on a Mac US keyboard):


> lowers = "1234567890-=/"
> uppers = "!@#$%^&*()_+?"
> lower x = lookup x (zip uppers lowers)
> upper x = lookup x (zip lowers uppers)

> upperChar x = case upper x of
>     Nothing -> mzero
>     Just y -> char y >> return x

> lowerChar x = case lower x of
>     Nothing -> mzero
>     Just y -> char y >> return x

A version of penalty wrapped for the parser monad:


> avoid :: Parser ()
> avoid = Parser $ \cs -> do
>    penalty
>    return ((),cs)

Reading keys on the assumption that the shift key may have failed:


> keyChar x = char x `mplus` (avoid >> upperChar x) `mplus` (avoid >> lowerChar x)

> digit = do
>     x <- foldl mplus mzero (map keyChar "0123456789")
>     return (fromIntegral (ord x-ord '0'))

> number1 :: Integer -> Parser Integer
> number1 m = return m `mplus` do
>     n <- digit
>     number1 (10*m+n)

> number :: Parser Integer
> number = do
>     n <- digit
>     number1 n

> chainl :: Parser a -> Parser (a -> a -> a) -> a -> Parser a 
> chainl p op a = (p `chainl1` op) `mplus` return a 
> chainl1 :: Parser a -> Parser (a -> a -> a) -> Parser a 
> p `chainl1` op = do {a <- p; rest a} 
>     where 
>     rest a = (do
>         f <- op
>         b <- p
>         rest (f a b)) `mplus` return a

Optional parentheses:


> shouldHave c = keyChar c `mplus` (avoid >> return c)

And the main part of the expression grammar:


> expr = term `chainl1` addop 
> term = monomial `chainl1` mulop 
> monomial = factor `chainl1` powop
> factor = number `mplus` do {shouldHave '('; n <- expr; shouldHave ')'; return n} 
> powop = keyChar '^' >> return (^)
> addop = do {keyChar '+'; return (+)} `mplus` do {keyChar '-'; return (-)} 
> mulop = do {keyChar '*'; return (*)} `mplus` do {keyChar '/'; return (div)}

Match the end of a string:


> end :: Parser ()
> end = Parser $ \cs ->
>    if null cs then P [[((),"")]] else mzero

We can test it out with:


> completeExpr = do
>   n <- expr
>   end
>   return n

> ex2 = parse completeExpr "2^(1+3"

When we run this we get no error-free parsing but we do get 3 readings with one error. One comes from reading the '(' as 9, one comes from inserting the missing ')' at the end and one comes from inserting ')' after '1'. Note that even for complex expressions we'll quickly find a 1- or 2-error parsing. For the regular list monad we might never get a parsing because there are an infinite number of ways of inserting parentheses.

Anyway, that was just a toy parsing problem. But a more complex application comes to my mind. Some written languages are tricky to parse because their orthography doesn't fully capture the phonetics of the original language, because there are few or no indicators of sentence or even word breaks, and because they have numerous optional orthographic and grammatical rules and use a script whose individual characters are occasionally hard to reliably identify. In such a situation it's good to have a parser driven by heuristics about what is likely to be intended and the penalty list monad might serve well. Here's an example of such a language.

Update: I forgot to add some connections to previous monads I've talked about:

PList is a variation of the convolution monad I described here. It deals with the "wrong category" aspect so it is a true Haskell monad. Penalty lists form some kind of dual to the convolution comonad.
It has much in common with this monad. That monad doesn't do anything smart about ordering searches but it does have the neat ability to 'fuse' different branches of a search so that different ways to arrive at the same place don't add to the combinatorial explosion. It's good for searches where you want to know what the minimum penalty is to get somewhere, but don't care what the best path actually is.

Also, in response to a comment on #haskell I've made the join example more complex so it's easier to generalise from it.

Automata and the A-D-E classification.

2009-06-20T13:01:00.000-07:00

Introduction

The A-D-E classification is a strange ubiquitous pattern that appears in many branches of mathematics. Typically it appears when you try to classify certain types of mathematical construction. If the A-D-E classification applies then you end up with two infinite sequences of cryptically named objects (A₁,A₂,A₃,...) and (D₁,D₂,D₃,...) as well as three leftover objects called E₆, E₇ and E₈. Unfortunately, most of these objects and their classifications are tricky to define using only elementary mathematics. However, there is one type of object that is classified in this way that can be given a relatively straightforward computational description involving a little linear algebra and assuming you know a tiny bit about automata.

But first: why care about the A-D-E classification? Well I tried to say a little bit about how symmetries relate to nature a while back. Certain types of possible symmetry of particle physics can be classified the A-D-E way. The symmetry group corresponding to E8 is the now famous exceptional group E8. I won't be able to get to an explanation of how groups are involved. But at least I'll be able to give a hint about the bigger picture that E8 is part of.

Non-deterministic Finite State Automata

Here's a diagram representing a very simple non-deterministic finite state automaton (NDA):

It can be in one of two states. When in state A it can transition to state B and in state B it can only transition back to state B, but it can do so in two different ways.

Vector Automata

Now I'll introduce a more general kind of automaton: a vector automaton (VA). (I made that term up, it's not meant to correspond with anyone else's terminology.) Every vector automaton is built from an NDA. But each state corresponds to a finite dimensional vector space and each transition corresponds to a linear function mapping from the vector space of the source state to the vector space of the destination state. We could turn the above example into a VA by assigning a 1D vector space V_A to A, a 2D vector space to V_B and defining linear functions:

f : (x) -> (x,0)
g : (x,y) -> (-y,x)
h : (x,y) -> (y,-x)

A VA is just like an NDA in that it transitions from state to state according to the given transitions. But additionally it keeps track of a vector in the vector space corresponding to the current state. Each time it makes a transition the linear function corresponding to that function is applied to the vector. So in the example above, the NDA might start in state A with a scalar value x (ie. a 1D vector). When it makes its first transition its vector becomes the 2D vector (x,0) and after that each transition rotates the vector through 90 degrees clockwise or anticlockwise.

There's a lot of freedom in defining a VA given its underlying NDA. For each node you can pick any vector space you like of any finite dimension, and for each transition you can pick any linear function you like mapping between the source and target vector spaces.

Let's make this a little more formal. An NDA is a finite set of states combined with a finite set of transitions. Each transition has a source and destination, each of which is a state. That's it. You're allowed any finite number of transitions between states and a transition can have the same state as source and destination.

A VA is an NDA combined with a finite dimensional vector space attached to each state and a linear function for each transition such that the function maps from the vector space of the of the source to the vector space of the target.

Now consider this really simple NDA:

When we build a VA from this NDA, for convenience I'll call the vector spaces corresponding to A and B, A and B. And I'll call the linear function from A to B, f. How many VAs can be made from this VA? Clearly an infinite number. But a lot of them are very similar.

Suppose we assume A and B are 2-dimensional and write their elements as pairs (x,y). Suppose f(1,0)=u and f(0,1)=v and that u is not a multiple of v and both are non-zero. Then we can use u and v as a basis for B. If we write f as f' in this new basis we get f'(1,0)=(1,0) and f'(0,1)=(0,1). So by relabelling the basis of B we have actually revealed that an infinite number of choices for f reduce to the same thing apart from a change of basis in B.

Up to change of basis in A and B we find there are only three possibilities:
(1) u and v are distinct, non-zero, and not multiples of each other.
(2) u is non-zero but v is zero
(3) both u and v are zero

We started with an infinite number of possibilities for our particular choice of dimensions and ended up with just 3. We'll say that two VAs are equivalent if we can get one from the other by changing basis like this.

On the other hand, consider this NDA:

Let's choose A, B and C to be 1-dimensional vector spaces and define f, g and h as:

f(x) = x
g(x) = x
h(x) = λx

where λ is any real number. Then h(g(f(x)))=λx. We can choose any λ we like so we have an infinity of possibilities. No amount of basis change is going to change this fact. This is different from the case above because now we're comparing x and h(g(f(x))) which lie in the same vector space. So when our NDAs have loops in them, the space of possible VAs, for most choices of dimension, is infinite.

The sum of two VAs

If we have two VAs corresponding to the same NDA we can combine them together to make a single machine. The state is given by a state in the shared underlying NDA, but we now have a pair of vectors. Each time there is a transition we apply the pair of transforms to transform the two vector simultaneously. But we can encode a pair of vectors as a single vector simply by concatenating together the vector components in some basis. So this machine is just another VA. The dimension of the vector space for each state is simply the sum of the dimensions of the vector spaces in the original VAs. So, given two VAs we can sum them to get a third.

Given a VA for an NDA it may or may not be equivalent to the sum of two simpler VAs. If it isn't, it's said to be irreducible. If it is, then we can ask the same question about the simpler VAs. In this way, every VA is the sum of irreducible VAs.

Main Theorem

Now comes the main result I want to give:

If the graph underlying the NDA is one of the following list then (and only then) there is a finite number of inequivalent irreducible VAs for the NDA. All other VAs for that NDA are simply finite sums of machines from this finite set. Here's the list:

That's it! Weird huh? (Note that diagram X_n has n nodes.)

Strangely, those same diagrams (and a few more) appear (in a quite different way) when classifying the possible symmetries of fundamental particles. The symmetries are given the same names as these diagrams. And in just the same way we get those 'sporadic' symmetries leading up to E₈. Those same diagrams also arise in Catastrophe Theory.

I'd like to sketch the proof, at least in one direction, but I seem to have run out of time. Another day perhaps. But note that none of these diagrams have loops for the reasons I gave above, and I've already shown that A₂ gives only a finite number of choices for a certain choice of dimension.

Meanwhile, I should give the proper mathematician names for the things above. The diagrams listed above are examples of Dynkin diagrams. A non-deterministic automaton as described above is known as a quiver. A vector automaton is normally known as a representation of a quiver. And the theorem above is half of Gabriel's theorem.

Hashing Molecules

2009-06-07T10:03:00.000-07:00

Twenty or so years ago I worked for a pharmaceutical company that had a large database of compounds. That got me thinking about the problem of how to perform lookups based on molecular structures. If you can find a bunch of numbers that encapsulate the molecular structure then you can use them as database keys. But you need to ensure that the same molecule entered in two different ways gets mapped to the same numbers, and you'd like different molecules, such as stereoisomers, or even enantiomers to get mapped to different values.

That got me thinking and around then I came up with an idea for hashing molecules inspired by some of the mathematics I'd been doing not long before. I never got around to coding it up but twenty years later it dawned on me that I could easily do it using a similar technique to what I used for untangling tangles and translating trace diagrams. Anyway, as I've already given examples of how to translate diagrams to monadic expressions I'm going to skimp on the details in this post and just talk about things specific to molecular structures. For this post to make sense you have to understand those earlier posts.

So first comes the usual Haskell preamble:



> {-# LANGUAGE MultiParamTypeClasses,FlexibleInstances,FunctionalDependencies,GeneralizedNewtypeDeriving #-}

> module Main where

> import Data.HashTable
> import Data.List as L
> import Data.Int
> import Data.Array
> import GHC.Arr
> import qualified Data.Map as M
> import Control.Monad
> infixl 5 .+
> infixl 6 .*

I'm going to be using the vector space monad with a 4-dimensional vector space. This type labels the dimensions:



> data I = A | B | C | D deriving (Eq,Ord,Show,Ix,Enum)

n-valent atoms will be represented by functions that consume n-tuples. We'll start with a simple hash:



> c' (a,b,c,d) = hashString ("C" ++ show (a,b,c,d))

My prime motivation for using Haskell for this problem is that the code was super-easy to write and investigate. But it's inefficient. I'll talk about how to remedy this properly at the end. But for now I'm going to memoise many functions as arrays using a Memoisable type class with a memo method. So I'll be using c instead of c':

> c = memo $ \x -> symmetrise a4 c' x .* return ()

Note the use of the symmetrise function. The idea is that the 4 bonds coming out from (singly-bonded) carbon can be thought of as lying at the corners of a tetrahedron.

They have tetrahedral symmetry. So I'd like my hash to also have this symmetry so that, for example, c (i,j,k,l) == c (j,i,l,k). We can enforce this by summing over all 24 permutations of the arguments compatible with this symmetry, also known as A4. So a4 lists all of these permutations:



> a4 (i,j,k,l) = [
>    (i,j,k,l),
>    (j,i,l,k),
>    (k,l,i,j),
>    (l,k,j,i),

>    (i,k,l,j),
>    (i,l,j,k),

>    (k,j,l,i),
>    (l,j,i,k),

>    (j,l,k,i),
>    (l,i,k,j),

>    (j,k,i,l),
>    (k,i,j,l)
>    ]

And symmetrise performs the summation:



> symmetrise group f x = sum (map f (group x))

We can define other molecules too. Hydrogen is easy:



> h = memo $ \a -> hashString ("H" ++ show a) .* return ()

Oxygen has S2 symmetry:



> s2 (i,j) = [ (i,j), (j,i) ]
> o' (a,b) = hashString ("O" ++ show (a,b))
> o = memo $ \x -> symmetrise s2 o' x .* return ()

You can think of a carbon atom and a hydrogen atom, say, as a pair of arrays c_ijkl and h_i, and bonding them together as summation over a shared index. So methane would be the sum over i,j,k,l = A..D of c_ijklh_ih_jh_kh_l. To this end, define a bond as:



> bond :: V Int32 I
> bond = return A .+ return B .+ return C .+ return D

We can make H₂ like so:



> h2 = simp $ do
>    i <- bond
>    h ! i
>    h ! i

i <- bond makes i a bond which we then attach to two hydrogen atoms. Evaluating h2 will give us the hash for hydrogen gas.

Rather than dive straight into CH₄ we can construct some useful building blocks. Hydrogen with a bond already attached:



> h_ :: V Int32 I
> h_ = do
>    m <- bond
>    h ! m
>    return m

I'm using a trailing underscore _ to indicate a free bond.

ch2_ accepts one bond and returns another. Once memoised it is, in effect, just a 4x4 matrix and can be used to rapidly build chains.



> ch2_ = memo $ \i -> simp $ do
>    k <- bond
>    l <- bond
>    m <- bond
>    h ! l
>    h ! m
>    c ! (i,k,l,m)
>    return k

So here's code to make an alkyl chain with a free bind at the end.



> alkyl_ 0 = h_

> alkyl_ n = simp $ do
>    i <- alkyl_ (n-1)
>    ch2_ ! i

We can now make alkanes by attaching a hydrogen atom at the end and make methane as a special case:



> alkane n = simp $ alkyl_ n >>= (h ! )
> methane = alkane 1

Now a hydroxyl group:



> oh = memo $ \i -> simp $ do
>    j <- bond
>    h ! j
>    o ! (i,j)

And you can have a drink on me:



> ethanol = simp $ do
>    i <- alkyl_ 2
>    oh ! i

Carbon double bonds turn out to be straightforward. We can simply use a pair of bonds:



> doubleBond :: V Int32 (I,I)
> doubleBond = simp $ do
>    i <- bond
>    j <- bond
>    return (i,j)

Carbon-carbon double bonds tend not to twist and this is reflected in the hash. We'd have to apply symmetrise if we wanted twistable bonds.

We can make a pre-canned doubly bonded carbon atom pair:



> c_c = memo $ \(i,j,k,l) -> simp $ do
>    (m,n) <- doubleBond
>    c ! (i,j,m,n)
>    c ! (m,n,k,l)

So we can build ethene like this



> ethene = simp $ do
>    i <- h_
>    j <- h_
>    k <- h_
>    l <- h_
>    c_c ! (i,j,k,l)

Here's a methyl group with a free bond



> ch3_ = simp $ do
>    l <- h_
>    ch2_ ! l

So we can build a bunch more compounds



> propene = simp $ do
>    j <- ch3_
>    i <- h_
>    k <- h_
>    l <- h_
>    c_c ! (i,j,k,l)

> cisbut2ene = simp $ do
>    i <- ch3_
>    j <- h_
>    k <- ch3_
>    l <- h_
>    c_c ! (i,j,k,l)

> transbut2ene = simp $ do
>    i <- h_
>    j <- ch3_
>    k <- ch3_
>    l <- h_
>    c_c ! (i,j,k,l)

> cisbut2ene' = simp $ do
>    i <- h_
>    j <- ch3_
>    k <- h_
>    l <- ch3_
>    c_c ! (i,j,k,l)

> _2methylpropene = simp $ do
>    i <- ch3_
>    j <- ch3_
>    k <- h_
>    l <- h_
>    c_c ! (i,j,k,l)

> _2methylpropene' = simp $ do
>    i <- h_
>    j <- h_
>    k <- ch3_
>    l <- ch3_
>    c_c ! (i,j,k,l)

An interesting problem is building a benzene ring. Here's a first attempt with six free bonds:



> ring1 (p,q,r,s,t,u) = simp $ do
>    i <- bond
>    j <- bond
>    k <- bond
>    c_c ! (j,q,i,p)
>    c_c ! (i,u,k,t)
>    c_c ! (k,s,j,r)

The problem with that is that benzene rings are special. The electrons are 'delocalised' so that the ring has rotational symmetry. We need to sum over the two consistent ways to assign single and double bonds around the ring. For the more general case of interlocking benzene rings we must still sum over all consistent assignments.



> ring = memo $ \(p,q,r,s,t,u) -> simp $ ring1 (p,q,r,s,t,u) .+ ring1 (q,r,s,t,u,p)

And some more compunds:



> phenyl = memo $ \p -> simp $ do
>    q <- h_
>    r <- h_
>    s <- h_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> benzene :: V Int32 ()
> benzene = simp $ do
>    i <- bond
>    phenyl ! i
>    h ! i

> phenol :: V Int32 ()
> phenol = simp $ do
>    i <- bond
>    phenyl ! i
>    oh ! i

> toluene :: V Int32 ()
> toluene = simp $ do
>    i <- ch3_
>    phenyl ! i

> toluene' :: V Int32 ()
> toluene' = simp $ do
>    p <- h_
>    q <- h_
>    r <- h_
>    s <- ch3_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> toluene'' :: V Int32 ()
> toluene'' = simp $ do
>    p <- h_
>    q <- h_
>    r <- ch3_
>    s <- h_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> _1_2_dimethylbenzene = simp $ do
>    p <- ch3_
>    q <- ch3_
>    r <- h_
>    s <- h_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> _1_3_dimethylbenzene = simp $ do
>    p <- ch3_
>    q <- h_
>    r <- ch3_
>    s <- h_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> _1_4_dimethylbenzene = simp $ do
>    p <- ch3_
>    q <- h_
>    r <- h_
>    s <- ch3_
>    t <- h_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

> _1_5_dimethylbenzene = simp $ do
>    p <- ch3_
>    q <- h_
>    r <- h_
>    s <- h_
>    t <- ch3_
>    u <- h_
>    ring ! (p,q,r,s,t,u)

As we might hope, the three different ways to define toluene give the same result. We also discover that 1,3- and 1,5-dimethylbenzene are the same compound (or at least probably are, this isn't a perfect hash).



> main = do
>    print $ "toluene = " ++ show toluene
>    print $ "toluene' = " ++ show toluene'
>    print $ "toluene'' = " ++ show toluene''
>    print $ "_1_2_dimethylbenzene = " ++ show _1_2_dimethylbenzene
>    print $ "_1_3_dimethylbenzene = " ++ show _1_3_dimethylbenzene
>    print $ "_1_4_dimethylbenzene = " ++ show _1_4_dimethylbenzene
>    print $ "_1_5_dimethylbenzene = " ++ show _1_5_dimethylbenzene

Now I need to say something about performance. The above code is naive and performs many unnecessary summations. For example, hashing a long chain should only take time linear in its length. But using the above code indiscriminately could give you exponential time. A good implementation might take a divide and conquer approach: slice the molecule in half through a bunch of bonds, compute partial hashes for each half and then sew the halves together in time exponential in the number of bonds you sliced through. For the types of molecules I've seen in real pharmaceutical databases (say) this is actually pretty cheap if you're smart about the slicing. The hashes in the above code could probably be computed many thousands of times faster. As it is, you'll probably need to compile the above with optimisation.

I'm willing to bet that with small changes, and with suitable choice of matrices over the reals, we can get invariants of molecules that predict physical properties. These calulations are reminiscent of algorithms for various types of counting algorithm so at the very least they probably compute things that are meaningful from a statistical mechanics perspective.

Incidentally, this approach to stitching together atoms was inspired by an old paper by R. C. Penner on fatgraphs - nothing to do with chemistry. A few days ago he put a paper online about an application to modelling proteins.

Update: Since writing this I've found a name for what I'm doing. I'm converting a chemical structure diagram into a tensor network. There seems to be lots of literature on how to evaluate these efficiently. In effect, everything I've been doing in this blog in terms of converting diagrams to code is an example of evaluating a tensor network.

Now comes the memoisation class:



> class Memoisable ix where
>    memo :: (ix -> a) -> Array ix a

> instance Memoisable I where
>    memo f = array (A,D) [(i,f i) | i <- [A .. D]]

> instance Memoisable (I,I) where
>    memo f = array ((A,A),(D,D)) [(i,f i) |
>        p <- [A .. D],
>        q <- [A .. D],
>        let i = (p,q) ]

> instance Memoisable (I,I,I,I) where
>    memo f = array ((A,A,A,A),(D,D,D,D)) [(i,f i) |
>        p <- [A .. D],
>        q <- [A .. D],
>        r <- [A .. D],
>        s <- [A .. D],
>        let i = (p,q,r,s) ]

> instance Memoisable (I,I,I,I,I,I) where
>    memo f = array ((A,A,A,A,A,A),(D,D,D,D,D,D)) [(i,f i) |
>        p <- [A .. D],
>        q <- [A .. D],
>        r <- [A .. D],
>        s <- [A .. D],
>        t <- [A .. D],
>        u <- [A .. D],
>        let i = (p,q,r,s,t,u) ]

Missing instances from Data.Array:



> instance (Ix a1, Ix a2, Ix a3, Ix a4, Ix a5, Ix a6) => Ix (a1,a2,a3,a4,a5,a6)  where
>    range ((l1,l2,l3,l4,l5,l6),(u1,u2,u3,u4,u5,u6)) =
>      [(i1,i2,i3,i4,i5,i6) | i1 <- range (l1,u1),
>                             i2 <- range (l2,u2),
>                             i3 <- range (l3,u3),
>                             i4 <- range (l4,u4),
>                             i5 <- range (l5,u5),
>                             i6 <- range (l6,u6)]

>    unsafeIndex ((l1,l2,l3,l4,l5,l6),(u1,u2,u3,u4,u5,u6)) (i1,i2,i3,i4,i5,i6) =
>      unsafeIndex (l6,u6) i6 + unsafeRangeSize (l6,u6) * (
>      unsafeIndex (l5,u5) i5 + unsafeRangeSize (l5,u5) * (
>      unsafeIndex (l4,u4) i4 + unsafeRangeSize (l4,u4) * (
>      unsafeIndex (l3,u3) i3 + unsafeRangeSize (l3,u3) * (
>      unsafeIndex (l2,u2) i2 + unsafeRangeSize (l2,u2) * (
>      unsafeIndex (l1,u1) i1)))))

>    inRange ((l1,l2,l3,l4,l5,l6),(u1,u2,u3,u4,u5,u6)) (i1,i2,i3,i4,i5,i6) =
>      inRange (l1,u1) i1 && inRange (l2,u2) i2 &&
>      inRange (l3,u3) i3 && inRange (l4,u4) i4 &&
>      inRange (l5,u5) i5 && inRange (l6,u6) i6

And the same vector space monad I've used many times before. Strictly speaking, it's more like a semiring module monad as Int32 isn't a field.



> swap (x,y) = (y,x)

> class Num k => VectorSpace k v | v -> k where
>    zero :: v
>    (.+) :: v -> v -> v
>    (.*) :: k -> v -> v
>    (.-) :: v -> v -> v
>    v1 .- v2 = v1 .+ ((-1).*v2)

> data V k a = V { unV :: [(k,a)] } deriving (Show)

> reduce x = filter ((/=0) . fst) $ fmap swap $ M.toList $ M.fromListWith (+) $ fmap swap $ x
> simp (V x) = V (reduce x)

> instance (Ord a,Num k) => Eq (V k a) where
>  V x==V y = reduce x==reduce y

> instance (Ord a,Num k,Ord k) => Ord (V k a) where
>  compare (V x) (V y) = compare (reduce x) (reduce y)

> instance Num k => Functor (V k) where
>    fmap f (V as) = V $ map (\(k,a) -> (k,f a)) as

> instance Num k => Monad (V k) where
>    return a = V [(1,a)]
>    x >>= f = join (fmap f x)
>        where join x = V $ concat $ fmap (uncurry scale) $ unV $ fmap unV x
>              scale k1 as = map (\(k2,a) -> (k1*k2,a)) as

> instance Num r => MonadPlus (V r) where
>    mzero = V []
>    mplus (V x) (V y) = V (x++y)

> instance (Num k,Ord a) => VectorSpace k (V k a) where
>    zero = V []
>    V x .+ V y = V (x ++ y)
>    (.*) k = (>>= (\a -> V [(k,a)]))

> e = return :: Num k => a -> V k a
> coefficient b (V bs) = maybe 0 id (L.lookup b (map swap (reduce bs)))

Trace Diagrams with Monads

2009-05-16T13:34:00.000-07:00

Knot diagrams aren't the only kind of diagram that can be translated nicely into Haskell monad notation. Other types of diagram include Penrose Spin Networks, many kinds of Feynman Diagram, Penrose Tensor Notation, birdtracks and a type of closely related diagram I hadn't met before: Trace Diagrams.

I encourage readers to check out the Wikipedia page and associated papers on trace diagrams as they give a better tutorial than I could write. My aim here is to show how those diagrams can be translated directly into working code just like with knots.

As usual, this is literate Haskell so I need these lines:


> {-# LANGUAGE MultiParamTypeClasses,FunctionalDependencies,FlexibleInstances #-}

> module Main where

> import qualified Data.Map as M
> import Control.Monad
> infixl 5 .+
> infixl 6 .*

I'll reuse my vector space monad code from before and work in a 3D space with the axes labelled X, Y and Z.


> data Space = X | Y | Z deriving (Eq,Show,Ord)

We draw vectors as little boxes with connections emerging from them:

Now recall from my knot posts that we represent a diagram with m legs at the top and n legs at the bottom as an expression that takes an m-tuple as input and returns an n-tuple "in the monad" as output.

Vectors can be represented as elements of V Float Space, for example:


> u,v,w :: V Float Space
> u = return X .- return Y
> v = return X .+ 2.* return Y
> w = return Y .- return Z

I could have emphasised that there are zero inputs at the top by using type signature () -> V Float Space instead.

Given two vectors we can form their dot product. The dot product itself is represented by a little u-shaped curve:

So the dot product of v and w is drawn as:

(The i and j are just so you can see what corresponds to what in the code below. They're not really part of the diagram.)

We can implement the dot product as


> cup :: (Space,Space) -> V Float ()
> cup (i,j) = case (i,j) of
>    (X,X) -> return ()
>    (Y,Y) -> return ()
>    (Z,Z) -> return ()
>    otherwise -> 0 .* return ()

and compute an example using


> vdotw = do
>     i <- v
>     j <- w
>     cup (i,j)

We hook up legs of the diagram using corresponding inputs and outputs in the code.

Now consider this diagram:

If we attach another vector to the free leg then we get the dot product. So this object is a thing that maps vectors to scalars. Ie. it's a dual vector. So dual vectors are represented by diagrams with a free leg at the top. We can redraw this diagram:
In other words, turning a vector v upside down turns it into a dual vector that takes w to the dot product of v and w. Here's some code for making the dual of v.


> dual :: V Float Space -> Space -> V Float ()
> dual v i = do
>     j <- v
>     cup (i,j)

We can also consider cross products. These take two vectors as input and output one. So we're looking at a diagram with two legs at the top and one at the bottom. We'll use a bold dot to represent one of these:

Here's the implementation:


> cross :: (Space,Space) -> V Float Space
> cross (X,Y) = return Z
> cross (Y,Z) = return X
> cross (Z,X) = return Y

> cross (Y,X) = (-1) .* return Z
> cross (Z,Y) = (-1) .* return X
> cross (X,Z) = (-1) .* return Y

> cross _ = mzero

We can form a triple product u.(v×w) like this:

We can then abstract out the triple product bit that looks like this:

Implementing it as:


> trident :: (Space,Space,Space) -> V Float ()
> trident (i,j,k) = do
>    l <- cross (i,j)
>    cup (l,k)

Remember that if u, v and w give the rows of a 3x3 matrix, then u.(v×w) is the determinant of that matrix.

We can also define a dot product for dual vectors that we can draw like this:

The code looks like this:


> cap :: () -> V Float (Space,Space)
> cap () = return (X,X) .+ return (Y,Y) .+ return (Z,Z)

We can now combine the two dot products in a diagram like this:

We can write that as:


> cupcap i = do
>      (j,k) <- cap ()
>      cup (i,j)
>      return k

We'd hope that we could pull this diagram taut and get the identity linear map. If you try applying cupcap to X, Y and Z you'll see it has exactly the same effect as return, which does indeed represent the identity.

(If you allow me to digress, I'll point out that there's something really deep going on with this almost trivial looking identity. It represents the identity map in the sense that it copies the input i to the output k. Imagine we were dealing with the trivial monad, ie. the one that just wraps values. Then no matter how cup and cap were implemented it would be impossible for k to be a copy of i. If you follow the flow of information through that code then i disappears into cup and k is read from cap without it seeing i. If we read from top to bottom we can think of cap as emitting a pair of objects and of cup as absorbing two. There is no way that any information about i can be communicated to k. But in the vector space monad, k can depend on i. As I've mentioned a few times over the years, the universe is described by quantum mechanics which can be described using the vector space monad. Amazingly the above piece of code, or at least something like it, can be physically realised in terms of particles. It describes a process that is fundamentally quantum, and not classical. In fact, Coecke shows that this is a precursor to quantum teleportation in section 3c of this paper. You could also think in terms of information about i being sent back in time through the cap. That's the idea behind this paper on Effective Quantum Time Travel.)

Now we can make a fork by bending down the tines of the cross product:


> fork () = do
>    (i,j) <- cap ()
>    (k,l) <- cap ()
>    m <- cross (j,k)
>    return (i,l,m)

We can write matrices as boxes with a leg for input and a leg for output. Here's an example matrix:


> a :: Space -> V Float Space
> a X = 2 .* return X
> a Y = return Z
> a Z = (-1) .* return Y

It rotates by 90 degrees around the X axis and scales the X axis by a factor of two.

With the help of our two dot products we can turn a matrix upside-down:

The corresponding code is:


> b :: Space -> V Float Space
> b l = do
>     (i,j) <- cap ()
>     k <- a j
>     cup (k,l)
>     return i

Turning a matrix upside down gives its transpose. You'll find that matrix B rotates in the opposite direction to A but still scales by a factor of two.

Surprisingly, 3! times the determinant of a 3x3 matrix A can be represented by this diagram:

So we can write a determinant routine as follows:


> det a = do
>     (i,j,k) <- fork ()
>     i' <- a i
>     j' <- a j
>     k' <- a k
>     (1/6.0) .* trident (i',j',k')

(Again I've labelled the diagram so you can easily see what corresponds where in the code.)

I could keep going, but at this point I'll just defer to Elisha Peterson's paper. I hope that I've given you enough clues to be able to translate his diagrams into Haskell code, in effect giving semantics for his syntax. As an exercise, try writing code to compute the adjugate of a 3x3 matrix.

And a reminder: none of the above is intended as production-worthy code for working with 3-vectors. It is intended purely as a way to give a practical realisation of trace diagrams allow people to experimentally investigate their properties and make testable conjectures.

And now comes the library code needed to make the above code work:


> swap (x,y) = (y,x)

> class Num k => VectorSpace k v | v -> k where
>    zero :: v
>    (.+) :: v -> v -> v
>    (.*) :: k -> v -> v
>    (.-) :: v -> v -> v
>    v1 .- v2 = v1 .+ ((-1).*v2)

> data V k a = V { unV :: [(k,a)] }
> instance (Num k,Ord a,Show a) => Show (V k a) where
>   show (V x) = show (reduce x)

> reduce x = filter ((/=0) . fst) $ fmap swap $ M.toList $ M.fromListWith (+) $ fmap swap $ x

> instance (Ord a,Num k) => Eq (V k a) where
>  V x==V y = reduce x==reduce y

> instance (Ord a,Num k,Ord k) => Ord (V k a) where
>  compare (V x) (V y) = compare (reduce x) (reduce y)

> instance Num k => Functor (V k) where
>    fmap f (V as) = V $ map (\(k,a) -> (k,f a)) as

> instance Num k => Monad (V k) where
>    return a = V [(1,a)]
>    x >>= f = join (fmap f x)
>        where join x = V $ concat $ fmap (uncurry scale) $ unV $ fmap unV x
>              scale k1 as = map (\(k2,a) -> (k1*k2,a)) as

> instance Num r => MonadPlus (V r) where
>    mzero = V []
>    mplus (V x) (V y) = V (x++y)

> instance (Num k,Ord a) => VectorSpace k (V k a) where
>    zero = V []
>    V x .+ V y = V (x ++ y)
>    (.*) k = (>>= (\a -> V [(k,a)]))

> e = return :: Num k => a -> V k a
> coefficient b (V bs) = maybe 0 id (lookup b (map swap (reduce bs)))

The Three Projections of Doctor Futamura

2009-05-03T08:30:00.000-07:00

Introduction

The Three Projections of Futamura are a sequence of applications of a programming technique called 'partial evaluation' or 'specialisation', each one more mind-bending than the previous one. But it shouldn't be programmers who have all the fun. So I'm going to try to explain the three projections in a way that non-programmers can maybe understand too. But whether you're a programmer or not, this kind of self-referential reasoning can hurt your brain. At least it hurts mine. But it's a good pain, right?

So rather than talk about computer programs, I'll talk about machines of the mechanical variety. A bit like computer programs, these machines will have some kind of slot for inputting stuff, and some kind of slot where output will come out. But unlike computer programs, I'll be able to draw pictures of them to show what I'm talking about. I'll also assume these machines have access to an infinite supply of raw materials for manufacturing purposes and I'll also assume that these machines can replicate stuff - because in a computer we can freely make copies of data, until we run out of memory at least.

Minting coins

A really simple example of a machine is one that has a slot for inputting blanks, and outputs newly minted coins:

That's a dedicated $1 manufacturing machine. We could imagine that internally it stamps the appropriate design onto the blank and spits out the result.

It'd be more interesting if we could make a machine with another input slot that allowed us to input the description of the coin. By providing different inputs we could mint a variety of different coins with one machine. I'm going to adopt the convention that when we want to input a description we input a picture of the result we want. I'll draw pictures as rectangles with the subject inside them. Here's a general purpose machine manufacturing pound coins:

The same machine could make dollars, zlotys or yen. You could imagine this machine works by taking the description and then milling the coin CNC style. We call such a machine an interpreter. It interprets the instructions and produces its result.

The interpreter has a great advantage over the dedicated dollar mint. You make make any kind of coin. But it's going to run a lot slower. The dedicated minter can just stamp a coin in one go. The interpreter can't do this because every input might be different. It has to custom mill each coin individually. Is there a way to get the benefits of both types of machine? We could do this: take the coin description and instead of milling the coin directly we mill negative reliefs for both sides of the coin. We then build a new dedicated minting machine that uses these negatives to stamp out the coin. In other words we could make a machine that takes as input a coin description and outputs a dedicated machine to make that type of coin. This kind of machine making machine is called a compiler. It takes a set of instructions, but instead of executing them one by one, it makes a dedicated machine to perform them. Here's one in action:

So here are the two important concepts so far:

Interpreters: these take descriptions or instructions and use them to make the thing described.
Compilers: these take descriptions or instructions and use them to make a machine dedicated to making the thing described. The process of making such a machine from a set of instructions is known as compiling.

The Projections of Doctor Futamura help make clear the relationship between these kinds of things.

Specialisation

We need one more important concept: the specialiser. Suppose we have a machine that has two inputs slots, A and B. But now suppose that when we use the machine we find that we vary the kind of thing we put into slot B, but always end up putting the same thing into slot A. If we know that slot A will always get the same input then we could streamline the machine using our knowledge of the properties of A. This is similar to the minting situation - if we know we're always going to make $1 coins then we can dedicate our machine to that purpose. In fact, if we know that we're always going to input the same thing into slot A we don't even need slot A any more. We could just stick an A inside the machine and whenever the user inputs something to slot B, the machine would then replicate the A and then use it just as if it had been input.

In summary, given a machine with two slots A and B, and given some input suitable for slot A, we could redesign it as a machine with just a B slot that automatically, internally self-feeds the chosen item to A. But we can often do better than this. We don't need to self-feed stuff to slot A. We might be able to redesign the way the machine works based on the assumption that we always get the same stuff going into slot A. For example, in the minting example a dedicate $1 minter was more specialised than just a general purpose minter that interpreted the instructions for making a $1 coin. This process of customising a machine for a particular input to slot A is called specialisation or partial evaluation.

Now imagine we have a machine for automatically specialising designs for machines. It might have two slots: one for inputting a description for a two slot machine with slots A and B, and one for inputting stuff suitable for slot A. It would then print out a description for a customised machine with just a slot B. We could call it a specialisation machine. Here is one at work:

It's converting a description of a two input machine into a description of a one input machine.

The First Projection

The process of specialisation is similar to what I was talking about with dedicated minting machines. Rather than just have a similarity we can completely formalise this. Note that the interpreter above takes two inputs. So the design for an interpreter could be fed into the first input of a specialiser. Now we feed a description the coin we want into slot B. The specialiser whirrs away and eventually outputs a description of a machine that is an interpreter that is dedicated to making that one particular coin. The result will describe a machine with only one input suitable for blanks. In other words, we can use a specialiser as a compiler. This is the first of Doctor Futamura's Projections. Here's a picture of the process at work:

What this shows is that you don't need to make compilers. You can make specialisers instead. This is actually a very practical thing to do in the computing world. For example there are commercial products (I'm not connected with that product in any way) that can specialise code to run on a specific architecture like CUDA. It's entirely practical to convert an interpreter to a compiler with such a tool. By writing a specialiser, the purveyors of such tools allow third parties to develop their own compilers and so this is more useful than just writing a dedicated compiler.

The Second Projection

Time to kick it up a notch. The first input to the specialiser is a description of a two input machine. But the specialiser is itself a two input machine. Are you thinking what I'm thinking yet? We could stuff a description of a specialiser into the specialiser's own first input! In the first projection we provided an interpreter as input to the specialiser. If we know we're always going to want to use the same interpreter then we could streamline the specialiser to work specifically with this input. So we can specialise the specialiser to work with our interpreter like this:

But what is that machine whose description it has output? An interpreter takes as input a description of how to operate on some stuff, like turning blanks into coins. In effect, the output machine has the interpreter built into it. So it takes descriptions and outputs a machine for performing those instructions. In other words it's a compiler. If the specialiser is any good then the compiler will be good too. It won't just hide an interpreter in a box and feed it your description, it will make dedicated parts to ensure your compiler produces a fast dedicated machine. And that is Doctor Futamura's Second Projection.

The Third Projection

But we can go further. The specialiser can accept a description of a specialiser as its first input. That means we can specialise it specifically for this input. And to do that, we use a specialiser. In other words we can feed a descrption of a specialiser into both inputs of the specialiser! Here we go:

But what is the X machine that it outputs? In the second projection we pass in an interpreter as the second argument and get back a compiler. So the third projection gives us a dedicated machine for this task. The X machine accepts the description of an interpreter as input and outputs the description of a compiler. So the X machine is a dedicated interpreter-to-compiler converter. And that is the Third Projection of Doctor Futamura.

If we have a specialiser we never need to make a compiler again. We need only design interpreters that we can automatically convert to compilers. In general it's easier to write interpreters than compilers and so in principle this makes life easier for programmers. It also allows us to compartmentalise the building of compilers. We can separate the interpreter bit from the bit that fashions specific parts for a task. The specialiser does the latter so our would-be compiler writer can concentrate on the former. But who would have guessed that passing a specialiser to itself twice would give us something so useful?

Summary

So here are the projections:

Compiling specific programs to dedicated machines.
Making a dedicated machine for compilation from an interpreter.
Making a machine dedicated to the task of converting interpreters into compilers.

There are lots of variations we can play with. I've just talked about descriptions of things without saying much about what those descriptions look like. In practice there are lots of different 'languages' we can use to express our descriptions. So variations on these projections can generate descriptions in different languages, possibly converting between them. We might also have lots of different specialisers that are themselves optimised for specific types of specialisation. The Futamura projections give interesting ways to combine these. And there are also variations for generating dedicating machines for other tasks related to compiling - like parsing the descriptions we might feed in as input.

If you want to read more on this subject there's a whole book online with example code. They're not easy things to design.

I think that specialisation is a killer feature that I'd like to see more of. Present day compilers (and here I'm talking about computers, not machines in general) are hard-coded black boxes for the task of compilation. They're not very good at allowing you to get in there and tweak the way compilation occurs - for example if you want to generate code according to a strategy you know. Specialisation is a nice alternative to simply bolting an API onto a compiler. It would make it easy for anyone to write optimising and optimised compilers for their own languages and combine such compilers with interpreters for interactive instead of offline compilation.

I learnt about this stuff, as well as lots of other stuff in my blog, from the excellent Vicious Circles. The theory is closely related to the theory of writing quines that I used for my three language quine.

And if you keep your ears to the ground you can hear rumours of a fabled fourth projection...

Faster than a speeding photon

2009-04-04T09:38:00.000-07:00

How to outrun a photon

I thought it would be fun to try to give a readable account of Unruh effect. It's a surprising phenomenon, and there isn't universal agreement over what exactly the theory predicts, let alone whether the effect has ever been observed. It has important implications for physics and philosophy and may even give a way to test some aspects of quantum gravity in the lab.

One way to start the story is consideration of this problem: if a photon is speeding towards you, can you outrun it? Let's simplify things a bit so that we're considering motion in one dimension.

If we're confined to one dimension, we can't dodge the photon, we can only hope to remain ahead of it. As the only things that can travel at the speed of light, c, are massless things like photons, it seems that there is no hope for a massive thing like a person in a spaceship to avoid it. The photon will always be faster than you, and so it'll catch you.

But in theory you can outrun a photon! Do you see the flaw in the above reasoning that made it seem impossible?

The best way to make things clear is to draw a diagram. We'll plot some graphs of position vs. time for some photons and spaceships. We'll have time going up the vertical axis and position along the horizontal axis. Here's an example:

I've chosen units so that one second on the vertical axis is drawn the same size as one light-second on the horizontal axis. The net result is that photons always travel at 45 degree angles to the axes. Massive objects, that travel slower than light, are confined to travel on courses that have angles of smaller than 45 degrees with respect to the vertical axis. The path of the photon is the diagonal black line and the path of a spaceship is in red. It starts to the right of the photon but as we move up the time axis the photon eventually catches up with it.

If the spaceship travels faster then it will follow an angle closer to 45 degrees. Here are a pair of paths corresponding to faster spaceships:

The faster the ship is, the further it gets before the photon catches up. But we're just putting off the inevitable. It seems that whatever we do, the photon will always catch up.

But there's a hidden assumption in the above. By drawing straight lines for the spaceship I was assuming it was travelling at a constant velocity. But there's no reason for that to be true. Here's a different path the spaceship could follow:

At no point does the red path of the spaceship meet the black path of the photon. And yet at no point does the red path reach 45 degrees to the vertical axis. In other words, the spaceship never travels at the speed of light, and yet the photon never catches up with it. Spaceships can outrun photons!

So what kind of path is that? It's actually a hyperbola and it corresponds to a spaceship accelerating at a constant rate. You might wonder how it can be constant acceleration when the speed of the spaceship never exceeds that of light. From an external observer's point of view, after a while it does look like the ship is travelling at a more or less constant velocity close to the speed of light. But from the point of view of someone on the spaceship it feels exactly like constant acceleration. So that is the path that would be taken by a spaceship with its thrusters firing at a constant rate.

An event horizon!

I chose that path so that the spaceship stays just in front of the photon. A photon that starts slightly to the right will eventually catch up with the ship. But photons starting further to the left of the ship will never reach it. This means that absolutely nothing starting to the left of the black photon path can ever catch the ship. That should sound familiar. It's exactly like a black hole. From the point of view of someone on the ship, the diagonal black line is exactly like the event horizon of a black hole. Nothing to the left of it can ever be seen by observers on the ship.

What does it look like if an observer in the ship watches something that crosses the event horizon? Here's another diagram:

Again, the black diagonal line is the path of a photon, which we know is a bit like an event horizon. The blue line is the path of an object at rest. In effect, it's falling over our apparent event horizon. Of course the blue object doesn't see any unusual phenomenon on approaching the event horizon because there's nothing really there - it's only something seen by observers in the ship. The blue object emits a series of photons (shown in green) at equal intervals. As long as these photons are emitted before the event horizon they eventually catch up with the red spaceship. But notice how they arrive at more and more widely spaced intervals. A photon released exactly at the event horizon never reaches the ship. So the viewers on the ship see the spacing between the photons get longer and longer. They'll never see the blue object cross the event horizon they'll just see it getting closer and closer until eventually it appears to freeze. Again, this is just like a black hole event horizon.

The universe looks pretty weird from the point of view of a constantly accelerating observer. Half of it is simply missing behind an event horizon. But that's just the start of the weirdness. When we throw Quantum Mechanics into the mix something much weirder happens.

Matter from a vacuum

It's well known that physicists expect black holes to emit particles as Hawking radiation. But our accelerating observer sees something like a black hole, so we might expect them to see something like Hawking radiation. We also know that an observer at rest sees no event horizon. Which means that we might predict that accelerating observers see particles that observers at rest don't. Can we take such a prediction seriously?

Let's look a bit more closely at this. According to a popular view of quantum mechanics, the vacuum is teeming with vacuum fluctuations - ephemeral particle-antiparticle pairs that briefly come into existence and then annihilate each other. In the diagram below I've drawn some of these events:

As we follow up the time axis, pairs of (complementary coloured) particles come into existence and then annihilate each other. These events are so fleeting that they have no effect on our particle detectors and we see a vacuum. But note that I've drawn one of these events straddling our apparent event horizon. From the point of view of an accelerating observer this looks like a pair of particles coming into existence, but because of the argument I sketched above, they seem to freeze near the event horizon. In other words, to an accelerating observer these fleeting events are no longer fleeting, they look like real particles coming into existence and sticking around forever. Accelerating observers appear to see particles in a vacuum!

What I've described above is absolutely not a rigourous argument. But amazingly, when you use the machinery of quantum field theory you end up making exactly the same prediction: accelerating observers see particles. This is known as the Unruh effect. When you do this properly you can compute a bit more detail. It turns out that the energies of the particles are random with exactly the same distribution as black body radiation. In other words, the vacuum looks like it has a glow corresponding to a particular temperature that is proportional to the acceleration. But it's not a bright glow. You need to accelerate at about 10²⁰ m/s² before the temperature appears to be 1K. Building a thermometer that can survive such accelerations is no mean feat. So it looks like the Unruh effect is a curiosity that might never be observed in the lab.

But it has been suggested that the Unruh effect has already been observed. There aren't many things that can survive that kind of acceleration, but an electron can, and an electron can behave like a thermometer. Electrons in circular particle accelerators routinely undergo the kinds of accelerations we're talking about. They do so because they are driven by a magnetic field. Now an electron has spin, so you can think of it as a bit like a little electric current running round in a loop. That means an electron is like a little dipole electromagnet. Magnets in magnetic fields tend to want to line up along the field - that's how a compass works. So electrons that spend long enough in a particle accelerator, eg. those in a storage ring, should eventually line up with the field. Lining up like this is known as polarisation, and in this particular case it's known as he Sokolov Ternov effect. But when we look at electrons in a storage ring it turns out they're not quite completely lined up, they're slightly depolarised. This is easily explained by Unruh radiation - they're constantly accelerating and so they feel themselves to be in a hot environment. The continual interaction with this hot environment causes the electron spins to be a bit randomised, so they don't all line up nicely.

Unfortunately this isn't definitive evidence for Unruh radiation because when we carry out the full calculation of the Sokolov-Ternov effect it turns out that it predicts partial depolarisation anyway. Now it looks like we don't have evidence for the Unruh effect. But it's not that simple. The Unruh effect isn't a new effect made up by a physicist. It's a prediction based on a new way of looking at fairly conventional physics. The Sokolov-Ternov effect is also predicted from standard physics, just in a different frame of reference. So maybe the partial depolarisation predicted by this effect is in fact the very same thing as the Unruh effect, just looked at from a different point of view.

What does this mean philosophically? We're used to the idea that looking at things from different angles changes how they look. Einstein extended this notion to spacetime so space and time seem different to moving observers. The Unruh effect goes one step further. Whether or not an individual particle exists depends on your point of view. Do particles not have any kind of existence independently of how we look at them? And how does it look to an observer at rest watching an accelerating observer fly by with a thermometer. Do they see a thermometer apparently responding to nothing? Or do they also see the particles once they've interacted with a thermometer? It's all so weird that some physicists take the view that the notion of the particle is outdated and we should only be talking about quantum mechanical wavefunctions.

There's another reason why the Unruh effect is important. The Unruh effect doesn't involve General Relativity, which is all about curved spacetime. But it does use the same mathematical machinery so it gives a way to test out that mathematics. So even if we don't have laboratory black holes to play with, we may still be able to investigate the mathematical framework that predicts phenomena like the Hawking effect. But also note that at least one paper claims it's all a mathematical error and there is no Unruh effect in reality.

In recent years there has been a dramatic increase in the number of papers on the Unruh effect. I expect this trend is going to keep going for a while

References

I learnt about the Unruh effect from Wald's book Quantum field theory in curved spacetime and black hole thermodynamics.
The diagram showing particle creation/annihilation events straddling the event horizon came from Susskind and Lindesay's An Introduction to Black Holes, Information and the String Theory Revolution.
I tried to catch up with recent developments by reading some of the paper The Unruh effect and its applications.

Dinatural Transformations and Coends

2009-03-07T13:13:00.000-08:00

Abstract nonsense warning: my goal here is simply to understand what the definition of a coend means in the context of Haskell.


> {-# LANGUAGE Rank2Types,ExistentialQuantification #-}

Pick a type. Let's call it X. Now consider the type


forall a . (a -> X,a) -> X

What can we say about elements of this type? There's an obvious example of such a function, eval that applies the first element of the pair to the second.


> eval (f,x) = f x

It's also easy to construct others. What property do they all share? Well there's a nice way to answer this question. We can use Janis Voigtländer's Free Theorem Generator. After a bit of simplification it tells us that if f is in


forall a . (a -> X,a) -> X

then for any compatible g


f (x,g y) = f (x . g,y)

We can apply any function we like to the second argument and we get exactly the same result by pre-composing with the first argument.

This is a pretty strong property - it puts a big constraint on what f can do. It holds because any function of type forall a . (a -> X,a) -> X maps to a type X that makes no reference to the quantified a. It can't let any information about the type a escape. And this means that f has to somehow eliminate an element of a and and element of a->X. There's only one non-trivial way to do that: provide the former to the latter as a function argument. So f must factor as h . eval for some function h. The free theorem comes from the fact that eval (x,g y) = eval (x.g,y). Or, f could be the constant function, but that still factors in the same way.

Let's try something slightly more complex. Consider the type


forall a . (a -> X,[a]) -> X

Again we can ask the free theorem generator to give us a property. We're told that if f is of this type, then for any compatible g


f(x,map g y) = f(x . g,y)

The crucial point this time is that in order to eliminate the [a] we have to use fmap to apply the function of type a -> X. So any function of this type must factor through fmap, and that's why the free theorem follows.

Now it's time to put this in a more general framework. We can write both of the above examples as elements of type forall a . S a a -> X where S is some type constructor. In both cases we also have that S a is a functor. But it's also a cofunctor in its first argument. Intuitively this says that S a b contains or produces b's but consumes a's. Here's a class to express this:


> class Difunctor h where
>   lmap :: (b -> a) -> h a c -> h b c
>   rmap :: (c -> d) -> h a c -> h a d

(I made up the word Difunctor.)

For (co)functoriality we insist on the laws


lmap (f . g) = lmap g . lmap f
rmap (f . g) = rmap f . rmap g

We can make both of the examples above instances:


> data Ex1 x a b = Ex1 (a -> x) b

> instance Difunctor (Ex1 x) where
>   lmap f (Ex1 g x) = Ex1 (g . f) x
>   rmap f (Ex1 g x) = Ex1 g (f x)

> data Ex2 x a b = Ex2 (a -> x) [b]

> instance Difunctor (Ex2 x) where
>   lmap f (Ex2 g x) = Ex2 (g . f) x
>   rmap f (Ex2 g x) = Ex2 g (map f x)

Our free theorems were essentially about this type:


> type DiNatural s x = forall a . (s a a -> x)

and in both cases the free theorem could be written as


f . lmap g == f . rmap g

This is a highly non-trivial result. Look at the preconditions for this theorem: they say something about each of the arguments to s individually. And yet we deduce that the two arguments are in fact intimately connected. We can apply g using either lmap or rmap and get the same result. This property is known as dinaturality. More precisely, if for some X, f :: s a a -> X, where s is a difunctor, then f is dinatural if for all compatible g,


f . lmap g == f . rmap g

We have this theorem: if s is an instance of the Difunctor class (obeying its laws) then for any X, elements of DiNatural s X are dinatural.

(Compare with the definition of naturality: f . fmap g = fmap g . f.)

I'm tempted to call this "Dinaturality for Free" but there's already a paper by that name and I don't know what it's about. And note that I only claim it's a theorem. I don't actually know how to prove this uniformly for all difunctors, but I'd stake a beer on it. At least as long as we don't do any weird haskell coding (so our free theorems are always valid) and as long as we restrict ourselves to functions that always terminate.

In the examples I gave above the reason why this holds is related to the fact that any function of the given type can be factored as something following an evaluation type function. More generally, if s is a dinatural then if there is a single Y, and function i :: s a a -> Y, such that every function f :: s a a -> X, for any X, can be factored as f = h . i, then Y is said to be the coend of s.

There are different approaches to computing a coend. Above I've used inspection. The coend of (a -> X,a) is X. But there's also a kind of cheating approach where we can use an existential type to get the type uniformly for all dinaturals:


> data Coend s = forall a . Coend (s a a)

For an explanation of why that works, see Edward Kmett's post on Kan Extensions. For this example, this means we expect the types t and Coend (Ex1 t) to be isomorphic. Here is the isomorphism and its inverse:


> iso :: Coend (Ex1 t) -> t
> iso (Coend (Ex1 f x)) = f x

> iso' :: t -> Coend (Ex1 t)
> iso' x = Coend (Ex1 (const x) ())

I'll leave the proof that iso and iso' are mutual inverses to you.

I hope that gives some idea of what a coend is. Informally it captures the method by which a dinatural transformation of type forall a . s a a -> X is able to eliminate the quantified a. If you look through the Haskell libraries you'll find many dinaturals (or at least things that can be made dinatural through the use of uncurry).

This code and description was inspired by the discussion at nLab.

Beyond Monads

2009-02-14T18:09:00.000-08:00

The state monad gives an elegant way to thread state information through Haskell code. Unfortunately it has an annoying limitation: the state must have the same type throughout the monadic expression. In this post I want to look at how to fix this. Unfortunately, fixing State means it's no longer a monad, but we'll discover a new abstraction that replaces monads. And then we can look at what else this abstraction is good for. The cool bit is that we have to write virtually no new code, and we'll even coax the compiler into doing the hard work of figuring out what the new abstraction should be.

This is all based on an idea that has been invented by a bunch of people independently, although in slightly different forms. I'm being chiefly guided by the paper Parameterized Notions of Computation.

The problem with the state monad is that it is defined by


newtype State s a = State { runState :: s -> (a, s) }

The state going into and out of one of these values is the same, s. We can't vary the type of the state as we pass through our code. But that's really easy to fix, just define:


> import Prelude hiding (return,(>>=),(>>),(.),id,drop)
> import Control.Category

> newtype State s1 s2 a = State { runState :: s1 -> (a, s2) }

I can now just copy and paste the definitions (with name changes to avoid clashes) out of the ghc prelude source code


> return' a = State $ \s -> (a, s)
> m >>>= k  = State $ \s -> let
>   (a, s') = runState m s
>   in runState (k a) s'

> get   = State $ \s -> (s, s)
> put s = State $ \_ -> ((), s)

We don't have to change a thing! The old code exactly matches the new type. We can now write code using the new State:


> test1 = return' 1 >>>= \x ->
>           return' 2 >>>= \y ->
>           get >>>= \z ->
>           put (x+y*z) >>>= \_ ->
>           return' z           
> go1 = runState test1 10

But we're now also able to write code like:


> test2 = return' 1 >>>= \x ->
>           return' 2 >>>= \y ->
>           get >>>= \z ->
>           put (show (x+y*z)) >>>= \_ ->
>           return' z           
> go2 = runState test2 10

The state starts of as an Integer but ends up as a String.

Problem solved! Except that this definition of State doesn't give us a monad and so we lose the benefits of having an interface shared by many monads. Is there a new more appropriate abstraction we can use? Rather than scratch our heads over it, we can just ask ghci to tell us what's going on.


*Main> :t return'
return' :: a -> State s1 s1 a
*Main> :t (>>>=)
(>>>=) :: State s1 s11 t -> (t -> State s11 s2 a) -> State s1 s2 a

This immediately suggests a new abstraction:


> class ParameterisedMonad m where
>   return :: a -> m s s a
>   (>>=) :: m s1 s2 t -> (t -> m s2 s3 a) -> m s1 s3 a

> x >> f = x >>= \_ -> f

It's a lot like the usual Monad class except that we're now parameterising uses of this class with a pair of types. Our new >>= operator also has a compatibility condition on it. We can think of an element of m s1 s2 as having a 'tail' and 'head' living in s1 and s2 respectively. In order to use >>= we require the head of the first argument to match the tail given by the second argument.

Anyway, we have:


> instance ParameterisedMonad State where
>   return = return'
>   (>>=) = (>>>=)

We didn't really design this class, we just used what ghci told us. Will it turn out to be a useful abstraction?

First a category theoretical aside: in this post I talked about how monads were really a kind of abstract monoid. Well ParameterisedMonad is a kind of abstract category. If we were to implement join for this class it would play a role analogous to composition of arrows in a category. In a monoid you can multiply any old elements together to get a new element. In a category, you can't multiply two arrows together unless the tail of the second matches the head of the first.

Now we can generalise the writer monad to a ParameterisedMonad. But there's a twist: every monoid gives rise to a writer. This time we'll find that every category gives rise to a ParameterisedMonad. Here's the definition. Again, it was lifted straight out of the source for the usual Writer monad. The main change is replacing mempty and mappend with id and flip (.).


> data Writer cat s1 s2 a = Writer { runWriter :: (a,cat s1 s2) }
> instance (Category cat) => ParameterisedMonad (Writer cat) where
>   return a = Writer (a,id)
>   m >>= k  = Writer $ let
>      (a, w)  = runWriter m
>      (b, w') = runWriter (k a)
>      in (b, w' . w)
> tell w = Writer ((),w)
> execWriter m = snd (runWriter m)

It's just like the usual Writer monad except that the type of the 'written' data may change. I'll borrow an example (modified a bit) from the paper. Define some type safe stack machine operations that are guaranteed not to blow your stack:


> push n x = (n,x)
> drop (_,x) = x
> dup (n,x) = (n,(n,x))
> add (m,(n,x)) = (m+n,x)
> swap (m,(n,x)) = (n,(m,x))

We can now 'write' the composition of a bunch of these operations as a 'side effect':


> test3 = tell (push 1) >>
>         tell (push 2) >>
>         tell dup >>
>         tell add >>
>         tell swap >>
>         tell drop
> go3 = execWriter test3 ()

I guess there's one last thing I have to find. The mother of all parameterised monads. Again, we lift code from the ghc libraries, this time from Control.Monad.Cont. I just tweak the definition ever so slightly. Normally when you hand a continuation to an element of the Cont type it gives you back an element of the continuation's range. We allow the return of any type. This time the implementations of return and (>>=) remain completely unchanged:


> newtype Cont r1 r2 a = Cont { runCont :: (a -> r2) -> r1 }
> instance ParameterisedMonad Cont where
>   return a = Cont ($ a)
>   m >>= k  = Cont $ \c -> runCont m $ \a -> runCont (k a) c

> i x = Cont (\fred -> x >>= fred)
> run m = runCont m return

> test4 = run $ i (tell (push 1)) >>
>               i (tell (push 2)) >>
>               i (tell dup) >>
>               i (tell add) >>
>               i (tell swap) >>
>               i (tell drop)

> go4 = execWriter test4 ()

So what's going on here? The implementations of these instances require almost trivial changes to the original monads, or in two cases no changes at all apart from the type signature. I have my opinion: Haskell programmers have been using the wrong type class all along. In each case the type signature for return and >>= was too strict and so the functionality was being unnecessarily shackled. By writing the code without a signature, ghci tells us what the correct signature should have been all along. I think it might just possibly be time to consider making ParameterisedMonad as important as Monad to Haskell programming. At the very least, do-notation needs to be adapted to support ParameterisedMonad.

Update: You *can* use do-notation with ParameterisedMonad if you use the NoImplicitPrelude flag.

Update2: Some credits and links:

The Polystate Monad is one of the independent discoveries I mentioned above.
A more general approach to Parameterized Monads in Haskell.
A comment on Parameterized Monads that shows explicitly how to make this work with NoImplicitPrelude.
Oleg's Variable (type)state `monad'.
Wadler discovered this design pattern back in 1993 in Monads and composable continuations.

I didn't contribute anything, this article is just advocacy.

Beyond Regular Expressions: More Incremental String Matching

2009-01-31T14:19:00.000-08:00

In my last post I showed how to incrementally match long strings against regular expressions. I now want to apply similar methods to matching languages that can't be described by regular expressions. (Note that 'language' is just jargon for a set of strings that meet some criterion.) In particular, regular expressions can't be used to test a string for balanced parentheses. This is because we need some kind of mechanism to count how many open parentheses are still pending and a finite state machine can't represent arbitrary integers.

So let's start with a slightly more abstract description of what was going on last time so we can see the bigger picture. We were storing strings in balanced trees with a kind of 'measurement' or 'digest' of the string stored in the nodes of the tree. Each character mapped to an element of a monoid via a function called measure and you can think of the measurement function as acting on entire strings if you mappend together all of the measurements for each of the characters. So what we have is a function f :: String -> M taking strings to some type M (in the last post M was a type of array) with the properties


f (a ++ b) == f a `mappend` f b
f [] == mempty

By noticing that String is itself a monoid we can write this as


f (a `mappend` b) == f a `mappend` f b
f mempty == mempty

Anything satisfying these laws is called a monoid homomorphism, or just homomorphism for short.

So the technique I used worked like this: I found a homomorphism from String to some type with the useful property that for any string s, f s still contains all the information required to figure out if we're dealing with a member of our language. If f turns a string into something more efficient to work with then we can make our string matching more efficient.

Now I want to make the notion of "contains all the information required" more precise by considering an example. Consider strings that consist only of the characters ( and ). Our language will be the set of strings whose parentheses balance. In other words the total number of ( must match the total number of ), and as we scan from left to right we must never see more ) than (. For example, ()()() and ((()))() are in our language, but )()()( isn't. This language is called the Dyck language.

Suppose we're testing whether or not some string is in the Dyck language. If we see () as a substring then if we delete it from the string, it makes no difference to whether or not the string is in the Dyck language. In fact, if we see (()), ((())), (((()))) and so on they can all be deleted. On the other hand, you can't delete )( without knowing about the rest of the string. Deleting it from ()() makes no difference to its membership in the Dyck language, but deleting it from )(() certainly does.

So given a language L, we can say that two strings, x and y, are interchangeable with respect to L if any time we see x as a substring of another string we can replace it with y, and vice versa, without making any difference to whether the string is in the language. Interchangeable strings are a kind of waste of memory. If we're testing for membership of L there's no need to distinguish between them. So we'd like our measurement homomorphism to map all interchangeable strings to the same values. But we don't want to map any more strings to the same value because then we lose information that tells us if a string is an element of L. A homomorphism that strikes this balance perfectly is called the 'canonical homomorphism' and the image of the set of all strings under this homomorphisms is called the syntactic monoid. By 'image', I simply mean all the possible values that could arise from applying the homomorphism to all possible strings.

So lets go back to the Dyck language. Any time we see () we can delete it. But if we delete every occurence of () from a string then all we have left is a bunch of ) followed by a bunch of (. Let's say it's p of the former, and q of the latter. Every string of parentheses can be distilled down to a pair of integers ≥0, (p,q). But does this go far enough? Could we distill any further? Well for any choice of (p,q) it's a good exercise to see that for any other choice of (p',q') there's always a string in the Dyck language where if you have )^p(^q as a substring, replacing it with (p',q') gives you something not in the language. So you can't distill any further. Which means we have our syntactic monoid and canonical homomorphism. In this case the monoid is called the bicyclic monoid and we can implement it as follows:


> {-# LANGUAGE TypeSynonymInstances,FlexibleInstances,MultiParamTypeClasses #-}
> import Data.Foldable
> import Data.Monoid
> import Data.FingerTree hiding (fromList)
> import qualified Data.List as L

> data Bicyclic = B Int Int deriving (Eq,Show)

> hom '(' = B 0 1
> hom ')' = B 1 0

> instance Monoid Bicyclic where
>   mempty = B 0 0
>   B a b `mappend` B c d = B (a-b+max b c) (d-c+max b c)

Where did that code for mappend come from? Consider )^a(^b)^c(^d. We can delete () from the middle many times over.

Now we can more or less reproduce the code of last week and get a Dyck language tester. Once we've distilled a string down to (p,q) we only need to test whether or not p=q=0 to see whether or not our parentheses are balanced:


> matches' s = x==B 0 0 where
>   x = mconcat (map hom s)

> data Elem a = Elem { getElem :: a } deriving Show
> data Size = Size { getSize :: Int } deriving (Eq,Ord,Show)

> instance Monoid Size where
>    mempty = Size 0
>    Size m `mappend` Size n = Size (m+n)

> instance Measured (Size,Bicyclic) (Elem Char) where
>    measure (Elem a) = (Size 1,hom a)

> type FingerString = FingerTree (Size,Bicyclic) (Elem Char)

> insert :: Int -> Char -> FingerString -> FingerString
> insert i c z = l >< (Elem c <| r) where (l,r) = split (\(Size n,_) -> n>i) z

> string = empty :: FingerString

> matchesDyck string = snd (measure string)==B 0 0

> loop string = do
>   print $ map getElem (toList string)
>   print $ "matches? " ++ show (matchesDyck string)
>   print "(Position,Character)"
>   r <- getLine
>   let (i,c) = read r
>   loop $ insert i c string

> main = do
>   loop string

There's a completely different way to test membership of the Dyck language. Replace each ( with 1 and ) with -1. Now scan from left to right keeping track of (1) the sum of all the numbers so far and (2) the minimum value taken by this sum. If the final sum and the final minimal sum are zero, then we have matching parentheses. But we need to do this on substrings without scanning from the beginning in one go. That's an example of a parallel prefix sum problem and it's what I talked about here.

So here's an extended exercise: adapt the parallel prefix sum approach to implement incremental Dyck language testing with fingertrees. You should end up with a canonical homomorphism that's similar to the one above. It'll probably be slightly different but ultimately equivalent.

And here's an even more extended exercise: protein sequences are sequences from a 20 letter alphabet. Each letter can be assigned a hydrophobicity value from certain tables. (Pick whichever table you want.) The hydrophobicity of a string is the sum of the hydrophobicities of its letters. Given a string, we can give it a score corresponding to the largest hydrophobicity of any contiguous substring in it. Use fingertrees and a suitable monoid to track this score as the string is incrementally edited. Note how widely separated substrings can suddenly combine together as stuff between them is adjusted.

If you're interested in Dyck languages with multiple types of parenthesis that need to match you need something much more fiendish.

Fast incremental regular expression matching with monoids

2009-01-24T13:01:00.000-08:00

The Problem

Consider this problem: Fix a regular expression R. Suppose you have a string of length N. There's not much you can do about it, you'll likely have to scan all N characters to test so see if the string matches R. But once you've performed the test, how fast can you test the string again if a small edit is made to it? It seems that in general you'd have to rescan the entire string, or at least rescan from where the edit was made. But it turns out that you can do regular expression matching incrementally so that for many changes you might make to the string, you only require O(log N) time to recompute whether the string matches. This is true even if characters at opposite ends of the string interact to make a succcessful match. What's more, it's remarkably straightforward to implement if we make use of fingertrees and monoids.

I'm going to assume a bit of background for which resources can be found on the web: understanding some basics about monoids, understanding apfelmus's inspiring introduction to fingertrees or the original paper, and you'll need to be completely comfortable with the idea of compiling regular expressions to finite state machines.

As usual, this is literate Haskell. But you need to have the fingertree package installed and we need a bunch of imports.


> {-# LANGUAGE TypeSynonymInstances,FlexibleInstances,MultiParamTypeClasses #-}
> import qualified Data.Array as B
> import Data.Array.Unboxed as U
> import Data.Foldable
> import Data.Monoid
> import Data.FingerTree hiding (fromList)
> import qualified Data.List as L

A Finite State Machine

So let's start with an example regular expression: .*(.*007.*).*. We're looking for "007" enclosed between parentheses, but the parentheses could be millions of characters apart.

A standard technique for finding regular expressions is to compile them to a finite state automaton. It takes quite a bit of code to do that, but it is completely standard. So rather than do that, here's a finite state machine I constructed by hand for this regular expression:

I've used the convention that an unlabelled edge corresponds to any input that isn't matched by another labelled edge. We can express the transitions as a function fsm.


> fsm 0 '(' = 1
> fsm 0 _ = 0

> fsm 1 '0' = 2
> fsm 1 _ = 1

> fsm 2 '0' = 3
> fsm 2 _ = 1

> fsm 3 '7' = 4
> fsm 3 '0' = 3
> fsm 3 _ = 1

> fsm 4 ')' = 5
> fsm 4 _ = 4

> fsm 5 _ = 5

The initial state is 0 and a match corresponds to state 5.

We can test a string in the standard way using


> matches s = Prelude.foldl fsm 0 s==5

Try matches "(00 7)" and matches "He(007xxxxxxxxxxxx)llo".

We can thing of the inputs as being functions acting on the automaton. Each input character is a function that maps the automaton from one state to another. We could use the Haskell composition function, (.), to compose these functions. But (.) doesn't really do anything, f . g is just a closure that says "when the time comes, apply g and then f.

On the other hand we could tabulate our transition functions as follows:


> tabulate f = array (0,5) [(x,f x) | x <- range (0,5)]

We have one such tabulated function for each letter in our alphabet:


> letters = array (' ','z') [(i,tabulate (flip fsm i)) | i <- range (' ','z')]

Given two tabulated functions we can easily form the table of the composition function. In fact, our tabulated functions form a monoid with mappend for composition. I used unboxed arrays for performance:


> type Table = UArray Int Int

> instance Monoid Table where
>    mempty = tabulate id
>    f `mappend` g = tabulate (\state -> (U.!) g ((U.!) f state))

Note that we've cheated a bit. An object of type Table could be any array of Ints indexed by Int. But if we promise to only build arrays indexed by elements of [0..5] and containing elements of the same range then our claim to monoidhood is valid.

Given any string, we can compute whether it matches our regular expression by looking up the corresponding Table in our letters array, composing them, and then checking if the tabulated function maps the initial state 0 to the final state 5:


> matches' s = table!0==5 where
>   table = mconcat (map ((B.!) letters) s)

This is slower and more complex than our original implementation of matches. But what we've now done is 'tease out' a monoid structure from the problem. If we store a string as a sequence of characters represented by a fingertree, we can store in each subtree the element of Table corresponding to the substring it represents. Every time the tree is rebalanced we need to recompute the corresponding Tables. But that's fine, it typically involves only O(log N) operations, and we don't need to write any code, the fingertree will do it for us automatically. Once we've done this, we end up with a representation of strings with the property that we always know what the corresponding Table is. We can freely split and join such trees knowing that the Table will always be up to date.

The only slight complication is that I want to be able to randomly access the nth character of the tree. apfelmus explains that in his post. The change I need to make is that I'm going to use both the Size monoid and the Table monoid, so I need the product monoid.


> data Elem a = Elem { getElem :: a } deriving Show
> data Size = Size { getSize :: Int } deriving (Eq,Ord,Show)

> instance Monoid Size where
>    mempty = Size 0
>    Size m `mappend` Size n = Size (m+n)

We need to implement measure as in the fingertree paper:


> instance Measured (Size,Table) (Elem Char) where
>    measure (Elem a) = (Size 1,(B.!) letters a)

A Fingertree

And now we can define our strings as:


> type FingerString = FingerTree (Size,Table) (Elem Char)

The insertion routine is more or less what's in the paper:


> insert :: Int -> Char -> FingerString -> FingerString
> insert i c z = l >< (Elem c <| r) where (l,r) = split (\(Size n,_) -> n>i) z

Note how I project out the size from the product monoid in order to insert at the correct position.

Here's an example string. Adjust the length to suit your memory and CPU horsepower:


> fromList = L.foldl' (|>) empty
> string = fromList (map Elem $ take 1000000 $ cycle "the quick brown fox jumped over the lazy dog")

(I use a strict form of fromList to ensure the tree actually gets built.)

The actual match function simply projects out the second component of the monoid and again tests to see if it maps the initial state to the final state:


> matches007 string = ((U.!) (snd (measure string)) 0)==5

An Interactive Loop

I recommend compiling with optimisation, something like ghc --make -O5 -o regexp regexp.lhs:


> loop string = do
>   print $ "matches? " ++ show (matches007 string)
>   print "(Position,Character)"
>   r <- getLine
>   let (i,c) = read r
>   loop $ insert i c string

> main = do
>   loop string

Now you can run this interactively. Input values like (100,'f') to insert an 'f' at position 100. It can take a good few seconds to compute the initial tree, but after that the matching process is instantaneous. (Actually, the second match might take a few seconds, that's because despite the foldl' the tree hasn't been fully built.)

A suitable sample input is:


(100,'(')
(900000,')')
(20105,'0')
(20106,'0')
(20107,'7')

Discussion

Note there is quite an overhead for this example. I'm storing an entire Table for each character. But you can easily store chunks of string (like in a rope). This means that some chunks will be rescanned when a string is edited - but rescanning a 1K chunk, say, is a lot less expensive than scanning a gigabyte file in its entirety. Working in blocks will probably speed up the initial scan too, a much smaller tree needs to be built.

When Hinze and Patterson originally wrote the fingertree paper they were motivated by parallel prefix sum methods. Just about any parallel prefix algorithm can be converted to an incremental algorithm using fingertrees. This article is based on the idea of doing this with the parallel lexing scheme described by Hillis and Steele in their classic Connection Machine paper.

So why would you want to match against a fixed regular expression like this? Well this method extends to a full blown incremental lexer. This will lex quickly even if placing a character in a string changes the type of lexemes billions of characters away. See the Hillis and Steele paper for details.

Note there's nothing especially Haskelly about this code except that Haskell made it easy to prototype. You can do this in C++, say, using mutable red-black trees.

Haskell Monoids and their Uses

2009-01-17T13:47:00.000-08:00

Haskell is a great language for constructing code modularly from small but orthogonal building blocks. One of these small blocks is the monoid. Although monoids come from mathematics (algebra in particular) they are found everywhere in computing. You probably use one or two monoids implicitly with every line of code you write, whatever the language, but you might not know it yet. By making them explicit we find interesting new ways of constructing those lines of code. In particular, ways that are often easier to both read and write. So the following is an intro to monoids in Haskell. I'm assuming familiarity with type classes, because Haskell monoids form a type class. I also assume some familiarity with monads, though nothing too complex.

This post is literate Haskell so you can play with the examples directly.

Defining Monoids

In Haskell, a monoid is a type with a rule for how two elements of that type can be combined to make another element of the same type. To be a monoid there also needs to be an element that you can think of as representing 'nothing' in the sense that when it's combined with other elements it leaves the other element unchanged.

A great example is lists. Given two lists, say [1,2] and [3,4], you can join them together using ++ to get [1,2,3,4]. There's also the empty list []. Using ++ to combine [] with any list gives you back the same list, for example []++[1,2,3,4]==[1,2,3,4].

Another example is the type of integers, Integer. Given two elements, say 3 and 4, we can combine them with + to get 7. We also have the element 0 which when added to any other integer leaves it unchanged.

So here is a possible definition for the monoid type class:


class Monoid m where
    mappend :: m -> m -> m
    mempty :: m

The function mappend is the function we use to combine pairs of elements, and mempty is the 'nothing' element. We can make lists an instance like this:


instance Monoid [a] where
    mappend = (++)
    mempty = []

Because we want mempty to do nothing when combined with other elements we also require monoids to obey these two rules


a `mappend` mempty = a

and


mempty `mappend` a = a.

Notice how there are two ways to combine a and b using mappend. We can write a `mappend` b or b `mappend` a. There is no requirement on a monoid that these be equal to each other. (But see below.) But there is another property that monoids are required to have. Suppose we start with the list [3,4]. And now suppose we want to concatenate it with [1,2] on the left and [5,6] on the right. We could do the left concatenation first to get [1,2]++[3,4] and then form ([1,2]++[3,4])++[5,6]. But we could do the right one first and get [1,2]++([3,4]++[5,6]). Because we're concatenating at opposite ends the two operations don't interfere and it doesn't matter which we do first. This gives rise to the third and last requirement we have of monoids:


(a `mappend` b) `mappend` c == a `mappend` (b `mappend` c)

and you can summarise it with the slogan 'combining on the left doesn't interfere with combining on the right'. Notice how the integers, combined with +, also have this property. It's such a useful property it has a name: associativity.

That's a complete specification of what a monoid is. Haskell doesn't enfore the three laws I've given, but anyone reading code using a monoid will expect these laws to hold.

Some Uses of Monoids

But given that we already have individual functions like ++ and +, why would we ever want to use mappend instead?

One reason is that with a monoid we get another function called mconcat for free. mconcat takes a list of values in a monoid and combines them all together. For example mconcat [a,b,c] is equal to a `mappend` (b `mappend` c). Any time you have a monoid you have this quick and easy way to combine a whole list together. But note that there is some ambiguity in the idea behind mconcat. To compute mconcat [a,b,...,c,d] which order should we work in? Should we work from left to right and compute a `mappend` b first? Or should we start with c `mappend` d. That's one place where the associativity law comes in: it makes no difference.

Another place where you might want to use a monoid is in code that is agnostic about how you want to combine elements. Just as mconcat works with any monoid, you might want to write your own code that works with any monoid.

Explicitly using the Monoid type class for a function also tells the reader of your code what your intentions are. If a function has signature [a] -> b you know it takes a list and constructs an object of type b from it. But it has considerable freedom in what it can do with your list. But if you see a function of type (Monoid a) => a -> b, even if it is only used with lists, we know what kind of things the function will do with the list. For example, we know that the function might add things to your list, but it's never going to pull any elements out of your list.

The same type can give rise to a monoid in different ways. For example, I've already mentions that the integers form a monoid. So we could define:


instance Monoid Integer where
    mappend = (+)
    mempty = 0

But there's a good reason not to do that: there's another natural way to make integers into a monoid:



instance Monoid Integer where
    mappend = (*)
    mempty = 1

We can't have both of these definitions at the same time. So the Data.Monoid library doesn't make Integer into a Monoid directly. Instead, it wraps them with Sum and Product. It also does so more generally so that you can make any Num type into a monoid in two different ways. We have both


Num a => Monoid (Sum a)

and


Num a => Monoid (Product a)

To use these we wrap our values in the appropriate wrapper and we can then use the monoid functions. For example mconcat [Sum 2,Sum 3,Sum 4] is Sum 9, but mconcat [Product 2,Product 3,Product 4] is [Product 24].

Using Sum and Product looks like a complicated way to do ordinary addition and multiplication. Why do things that way?

The Writer Monad

You can think of monoids as being accumulators. Given a running total, n, we can add in a new value a to get a new running total n' = n `mappend` a. Accumulating totals is a very common design pattern in real code so it's useful to abstract this idea. This is exactly what the Writer monad allows. We can write monadic code that accumulates values as a "side effect". The function to perform the accumulation is (somewhat confusingly) called tell. Here's an example where we're logging a trace of what we're doing.


> import Data.Monoid
> import Data.Foldable
> import Control.Monad.Writer
> import Control.Monad.State

> fact1 :: Integer -> Writer String Integer
> fact1 0 = return 1
> fact1 n = do
>   let n' = n-1
>   tell $ "We've taken one away from " ++ show n ++ "\n"
>   m <- fact1 n'
>   tell $ "We've called f " ++ show m ++ "\n"
>   let r = n*m
>   tell $ "We've multiplied " ++ show n ++ " and " ++ show m ++ "\n"
>   return r

This is an implementation of the factorial function that tells us what it did. Each time we call tell we combine its argument with the running log of all of the strings that we've 'told' so far. We use runWriter to extract the results back out. If we run


> ex1 = runWriter (fact1 10)

we get back both 10! and a list of what it took to compute this.

But Writer allows us to accumulate more than just strings. We can use it with any monoid. For example, we can use it to count how many multiplications and subtractions were required to compute a given factorial. To do this we simply tell a value of the appropriate type. In this case we want to add values, and the monoid for addition is Sum. So instead we could implement:


> fact2 :: Integer -> Writer (Sum Integer) Integer
> fact2 0 = return 1
> fact2 n = do
>   let n' = n-1
>   tell $ Sum 1
>   m <- fact2 n'
>   let r = n*m
>   tell $ Sum 1
>   return r


> ex2 = runWriter (fact2 10)

There's another way we could have written this, using the state monad:


> fact3 :: Integer -> State Integer Integer
> fact3 0 = return 1
> fact3 n = do
>   let n' = n-1
>   modify (+1)
>   m <- fact3 n'
>   let r = n*m
>   modify (+1)
>   return r

> ex3 = runState (fact3 10) 0

It works just as well, but there is a big advantage to using the Writer version. It has type signature f :: Integer -> Writer (Sum Integer) Integer. We can immediately read from this that our function has a side effect that involves accumulating a number in a purely additive way. It's never going to, for example, multiply the accumulated value. The type information tells us a lot about what is going on inside the function without us having to read a single line of the implementation. The version written with State is free to do whatever it likes with the accumulated value and so it's harder to discern its purpose.

Data.Monoid also provides an Any monoid. This is the Bool type with the disjunction operator, better known as ||. The idea behind the name is that if you combine together any collection of elements of type Any then the result is Any True precisely when at least any one of the original elements is Any True. If we think of these values as accumulators then they provide a kind of one way switch. We start accumulating with mempty, ie. Any False, and we can think of this as being the switch being off. Any time we accumulate Any True into our running 'total' the switch is turned on. This switch can never be switched off again by accumulating any more values. This models a pattern we often see in code: a flag that we want to switch on, as a side effect, if a certain condition is met at any point.


> fact4 :: Integer -> Writer Any Integer
> fact4 0 = return 1
> fact4 n = do
>   let n' = n-1
>   m <- fact4 n'
>   let r = n*m
>   tell (Any (r==120))
>   return r

> ex4 = runWriter (fact4 10)

At the end of our calculation we get n!, but we are also told if at any stage in the calculation two numbers were multiplied to give 120. We can almost read the tell line as if it were English: "tell my caller if any value of r is ever 120". Not only do we get the plumbing for this flag with a minimal amount of code. If we look at the type for this version of f it tells us exactly what's going on. We can read off immediately that this function, as a "side effect", computes a flag that can be turned on, but never turned off. That's a lot of useful information from just a type signature. In many other programming languages we might expect to see a boolean in the type signature, but we'd be forced to read the code to get any idea of how it will be used.

Commutative Monoids, Non-Commutative Monoids and Dual Monoids

Two elements of a monoid, x and y, are said to commute if x `mappend` y == y `mappend` x. The monoid itself is said to be commutative if all of its elements commute with each other. A good example of a commutative monoid is the type of integers. For any pair of integers, a+b==b+a.

If a monoid isn't commutative, it's said to be non-commutative. If it's non-comuutative it means that for some x and y, x `mappend` y isn't the same as y `mappend` x, so mappend and flip mappend are not the same function. For example [1,2] ++ [3,4] is different from [3,4] ++ [1,2]. This has the interesting consequence that we can make another monoid in which the combination function is flip mappend. We can still use the same mempty element, so the first two monoid laws hold. Additionally, it's a nice exercise to prove that the third monoid law still holds. This flipped monoid is called the dual monoid and Data.Monoid provides the Dual type constructor to build the dual of a monoid. We can use this to reverse the order in which the writer monad accumulates values. For example the following code collects the execution trace in reverse order:


> fact5 :: Integer -> Writer (Dual String) Integer
> fact5 0 = return 1
> fact5 n = do
>   let n' = n-1
>   tell $ Dual $ "We've taken one away from " ++ show n ++ "\n"
>   m <- fact5 n'
>   tell $ Dual $ "We've called f " ++ show m ++ "\n"
>   let r = n*m
>   tell $ Dual $ "We've multiplied " ++ show n ++ " and " ++ show m ++ "\n"
>   return r

> ex5 = runWriter (fact5 10)

The Product Monoid

Suppose we want to accumulate two side effects at the same time. For example, maybe we want to both count instructions and leave a readable trace of our computation. We could use monad transformers to combine two writer monads. But there is a slightly easier way - we can combine two monoids into one 'product' monoid. It's defined like this:


instance (Monoid a,Monoid b) => Monoid (a,b) where
    mempty = (mempty,mempty)
    mappend (u,v) (w,x) = (u `mappend` w,v `mappend` x)

Each time we use mappend on the product we actually perform a pair of mappends on each of the elements of the pair. With these small helper functions:


> tellFst a = tell $ (a,mempty)
> tellSnd b = tell $ (mempty,b)

we can now use two monoids simultaneously:


> fact6 :: Integer -> Writer (String,Sum Integer) Integer
> fact6 0 = return 1
> fact6 n = do
>   let n' = n-1
>   tellSnd (Sum 1)
>   tellFst $ "We've taken one away from " ++ show n ++ "\n"
>   m <- fact6 n'
>   let r = n*m
>   tellSnd (Sum 1)
>   tellFst $ "We've multiplied " ++ show n ++ " and " ++ show m ++ "\n"
>   return r

> ex6 = runWriter (fact6 5)

If we had simply implemented our code using one specific monoid, like lists, our code would be very limited in its application. But by using the general Monoid type class we ensure that users of our code can use not just any individual monoid, but even multiple monoids. This can make for more efficient code because it means we can perform multiple accumulations while traversing a data structure once. And yet we still ensure readability because our code is written using the interface to a single monoid making our algorithms simpler to read.

Foldable Data

One last application to mention is the Data.Foldable library. This provides a generic approach to walking through a datastructure, accumulating values as we go. The foldMap function applies a function to each element of our structure and then accumulates the return values of each of these applications. An implementation of foldMap for a tree structure might be:


> data Tree a = Empty | Leaf a | Node (Tree a) a (Tree a)

> instance Foldable Tree where
>   foldMap f Empty = mempty
>   foldMap f (Leaf x) = f x
>   foldMap f (Node l k r) = foldMap f l `mappend` f k `mappend` foldMap f r

We can now use any of the monoids discussed above to compute properties of our trees. For example, we can use the function (== 1) to test whether each element is equal to 1 and then use the Any monoid to find out if any element of the tree is equal to 1. Here are a pair of examples: one to compute whether or not an element is equal to 1, and another to test if every element is greater than 5:


> tree = Node (Leaf 1) 7 (Leaf 2)

> ex7 = foldMap (Any . (== 1)) tree
> ex8 = foldMap (All . (> 5)) tree

Note, of course, that these expressions can be used, unmodified, with any foldable type, not just trees.

I hope you agree that this expresses our intentions in a way that is easy to read.

That suggests another exercise: write something similar to find the minimum or maximum element in a tree. You may need to construct a new monoid along the lines of Any and All. Try finding both in one traversal of the tree using the product monoid.

The foldable example also illustrates another point. The implementor of foldMap for the tree doesn't need to worry about whether the left tree should be combined with the central element before the right tree. Associativity means it can be implemented either way and give the same results.

Recap

Monoids provide a general approach to combining and accumulating values. They allow us to write code that is agnostic about the method we will use to combine values, and that makes our code more reusable. By using named monoids we can write type signatures that express our intentions to people reading our code: for example by using Any instead of Bool we make it clear just how our boolean value is to be used. And we can combine the monoid-based building blocks provided by Haskell libraries to build useful and readable algorithms with a minimum of effort.

Some final notes: mathematicians often refer to mappend as a 'binary operator' and often it's called 'multiplication'. Just like in ordinary algebra, it's often also written with abuttal or using the star operator, ie. ab and a*b might both represent a `mappend` b. You can read more about monoids at Wikipedia. And I wish I had time to talk about monoid morphisms, and why the list monoid is free (and what consequences that might have for how you can your write code), and how compositing gives you monoids and a whole lot more.

Rewriting Monadic Expressions with Template Haskell

2009-01-02T14:24:00.000-08:00

The goal today is to implement an impossible Haskell function. But as this is a literate Haskell post we need to get some boilerplate out of the way:


> module Test where

> import Language.Haskell.TH
> import Control.Monad
> import Control.Monad.Reader
> import Control.Monad.Cont
> import IO

> infixl 1 #

So back to the impossible function: all monads m come equipped with a function of type a -> m a. But it's well known that you can't "extract elements back out of the monad" because there is no function of type m a -> a. So my goal today will be to write such a function. Clearly the first line of code ought to be:


> extract :: Monad m => m a -> a

The rest of the post will fill in the details.

The type declaration tells us what role it plays syntactically, ie. it tells us how we can write code with extract that type checks. But what should the semantics be?

For the IO monad an answer is easy to guess. 1 + extract (readLn :: IO Int) should execute readLn by reading n integer from stdin, strip of the IO part of the the return value and then add 1 to the result. In fact, Haskell already has a function that does exactly that, unsafePerformIO. The goal here is to implement extract in a way that works with any monad.

What might we expect the value of 1 + extract [1,2,3] to be? The value of extract [1,2,3] surely must be 1, 2 or 3. But which one? And what happens if the list is empty? There really isn't any way of answering this while remaining *purely* functional. But if we were running code on a suitable machine we could fork three threads returning one of 1, 2 and 3 in each thread, and then collecing th results together in a list. In other words, we'd expect the final result to be [2,3,4]. This would make extract a lot like McCarthy's ambiguous operator.

So it's clear that we can't interpret extract as a pure function. But we could try implementing it on a new abstract machine. But there's another approach we could take. A couple of times I've talked about a function of type ~~a -> a where ~a is the type a -> Void and Void is the type with no elements. This corresponds to double negation elimination in classical logic. The Curry-Howard isomorphism tells us no such function can be implemented in a pure functional language, but we can translate expressions containing references to such a function into expressions that are completely pure. This is the so-called CPS translation. Anyway, I had this idea that we could do something limilar with monads so that we could translate expressions containing extract into ordinary Haskell code. Turns out there's already a paper on how to do this: Representing Monads by Andrzej Filinski.

To translate all of Haskell this way would be a messy business. But just for fun I thought I'd implement a simpler translation for a small subset of Haskell. It's simply this:

For a choice of monad m denote the translation of both types and values by T. So x::a becomes T(x)::T(a). The translation is simply:

T(a) = m a on types
T(f x) = T(f) `ap` T(x)
T(extract x) = join T(x)

The important thing is that extract of type m a -> a is replaced by join of type m (m a) -> m a.

A great way to translate Haskell code is with Template Haskell. So here's some code:


> (#) x y = liftM2 AppE x y

> rewrite :: Exp -> ExpQ

> rewrite (AppE f x) = do
>     e <- [| extract |]
>     if f==e
>         then [| join |] # rewrite x
>         else [| ap |] # rewrite f # rewrite x

Most of the rest is support for some forms of syntactic sugar. First the infix operators:


> rewrite (InfixE (Just x) f Nothing) =
>     [| fmap |] # return f # rewrite x

> rewrite (InfixE (Just x) f (Just y)) =
>     [| liftM2 |] # return f # rewrite x # rewrite y

> rewrite (InfixE Nothing f (Just y)) =
>     [| fmap |] # ([| flip |] # return f) # rewrite y

And list operations. For example [a,b,c] is sugar for a : b : c : [].


> rewrite (ListE []) = [| return [] |]

> rewrite (ListE (x:xs)) =
>     [| liftM2 |] # [| (:) |] # rewrite x # rewrite (ListE xs)

> rewrite x =
>     [| return |] # return x

> test :: ExpQ -> ExpQ
> test = (>>= rewrite)

extract itself is just a placeholder that is supposd to be translated away:


> extract = undefined

If the above is placed in a file called Test.lhs then you can try compiling the following code.


{-# LANGUAGE TemplateHaskell #-}

import Test
import Control.Monad
import Control.Monad.Reader
import Control.Monad.Cont

ex1 = $(test [|
        extract [1,2] + extract [10,20]
    |])
ex2 = runReader $(test [|
        extract ask + 7
    |]) 10
ex3 = $(test [|
        1+extract (readLn :: IO Int)
    |] )

The big omission is the lack of translation for lambda abstractions. I think that to get this right might requires translating all of the code using -|extract|- from the ground up, not just isolated expressions like those above. And like with CPS, you lose referential transparency and the order in which expressions are evaluated makes a difference.

Anyway, this is a partial answer to the question posed here on "automonadization".

The Mother of all Monads

2008-12-24T13:19:00.001-08:00

Suppose someone stole all the monads but one, which monad would you want it to be? If you're a Haskell programmer you wouldn't be too bothered, you could just roll your own monads using nothing more than functions.

But suppose someone stole do-notation leaving you with a version that only supported one type of monad. Which one would you choose? Rolling your own Haskell syntax is hard so you really want to choose wisely. Is there a universal monad that encompasses the functionality of all other monads?

I often find I learn more computer science by trying to decode random isolated sentences than from reading entire papers. About a year ago I must have skimmed this post because the line "the continuation monad is in some sense the mother of all monads" became stuck in my head. So maybe Cont is the monad we should choose. This post is my investigation of why exactly it's the best choice. Along the way I'll also try to give some insight into how you can make practical use the continuation monad. I'm deliberately going to avoid discussing the underlying mechanism that makes continuations work.

So let's start with this simple piece of code


> import Control.Monad.Cont

> ex1 = do
>   a <- return 1
>   b <- return 10
>   return $ a+b

I haven't specified the monad but in almost every case we'd expect the result to have something to do with the number 11. For the list monad we get [11], for the Maybe monad we get Just 11 and so on. For the Cont monad we get something that takes a function, and applies it to 11. Here's an example of its use:


> test1 = runCont ex1 show

ex1 is just a function that takes as argument show and applies it to 11 to give the string "11". Cont and runCont are just wrapping and unwrapping functions that we can mostly ignore.

We could have done that without continuations. So what exactly does the Cont monad give us here? Well let's make a 'hole' in the code above:

Whatever integer we place in the hole, the value of test1 will be the result of adding one and applying show. So we can think of that picture as being a function whose argument we shove in the hole. Now Haskell is a functional programming language so we expect that we can somehow reify that function and get our hands on it. That's exactly what the continuation monad Cont does. Let's call the function we're talking about by the name fred. How can we get our hands on it? It's with this piece code:


ex1 = do
  a <- return 1
  b <- Cont (\fred -> ...)
  return $ a+b

The ... is a context in which fred represents "the entire surrounding computation". Such a computaton is known as a "continuation". It's a bit hard to get your head around but the Cont monad allows you to write subexpressions that are able to "capture" the entirety of the code around them, as far as the function provided to runCont. To show that this is the case let's apply fred to the number 10:


> ex2 = do
>   a <- return 1
>   b <- Cont (\fred -> fred 10)
>   return $ a+b

> test2 = runCont ex2 show

The entire computation is applied to 10 and we get "11". Now you know what return does in this monad. But that's a convoluted way of doing things. What other advantages do we get? Well the expression for b can do whatever it wants with fred as long as it returns the same type, ie. a string. So we can write this:


> ex3 = do
>   a <- return 1
>   b <- Cont (\fred -> "escape")
>   return $ a+b

> test3 = runCont ex3 show

fred is completely ignored. The entire computation is thrown away and instead of applying show to a number, we simply return "escape". In other words, we have a mechanism for throwing values out of a computation. So continuations provide, among other things, an exception handling mechanism. But that's curious, because that's exactly what the Maybe monad provides. It looks like we might be able to simulate Maybe this way. But rather than do that, let's do something even more radical.


> ex4 = do
>   a <- return 1
>   b <- Cont (\fred -> fred 10 ++ fred 20)
>   return $ a+b

> test4 = runCont ex4 show

We've used fred twice. We've made the code around our "hole" run twice, each time executing with a different starting value. Continuations allow mere subexpressions to take complete control of the expressions within which they lie. That should remind you of something. It's just like the list monad. The above code is a lot like


> test5 = do
>   a <- return 1
>   b <- [10,20]
>   return $ a+b

So can we emulate the list monad? Well instead of converting our integer to a string at the end we want to convert it to a list. So this will work:


> ex6 = do
>   a <- return 1
>   b <- Cont (\fred -> fred 10 ++ fred 20)
>   return $ a+b

> test6 = runCont ex6 (\x -> [x])

We can avoid those ++ operators by using concat:


> ex7 = do
>   a <- return 1
>   b <- Cont (\fred -> concat [fred 10,fred 20])
>   return $ a+b

> test7 = runCont ex7 (\x -> [x])

But now you may notice we can remove almost every depepndence on the list type to get:


> ex8 = do
>   a <- return 1
>   b <- Cont (\fred -> [10,20] >>= fred)
>   return $ a+b

> test8 = runCont ex8 return

Note, we're using monad related functions, but when we do so we're not using do-notation. We can now do one last thing to tidy this up:


> i x = Cont (\fred -> x >>= fred)
> run m = runCont m return

And now we have something close to do-notation for the list monad at our disposal again:


> test9 = run $ do
>   a <- i [1,2]
>   b <- i [10,20]
>   return $ a+b

I hope you can see how this works. i x says that the continuation should be applied to x, not as an ordinary function, but with >>=. But that's just business as usual for monads. So the above should work for any monad.


> test10 = run $ do
>   i $ print "What is your name?"
>   name <- i getLine
>   i $ print $ "Merry Xmas " ++ name

The Grinch has been foiled and we see that the continuation monad really is the mother of all monads.

There are some interesting consequences of this beyond Haskell. Many languages with support for continuations should be extensible to support monads. In particular, if there is an elegant notation for continuations, there should be one for monads too. This is why I didn't want to talk about the underlying mechanism of the Cont monad. Different languages can implement continuations in different ways. An extreme example is (non-portable) C where you can reify continuations by literally flushing out all registers to memory and grabbing the stack. In fact, I've used this to implement something like the list monad for searching in C. (Just for fun, not for real work.) Scheme has call-with-current-continuation which can be used similarly. And even Python's yield does something a little like reifying a continuation and might be usable this way. (Is that's what's going on here? I haven't read that yet.).

This post was also inspired by this paper by Filinski. I haven't followed the details yet (it's tricky) but the gist is similar. I was actually looking at Filinski's paper because of something I'll mention in my next post.

An Approach to Algorithm Parallelisation

2008-11-29T16:44:00.000-08:00

The other day I came across the paper Parallelizing Complex Scans and Reductions lying on a colleague's desk. The first part of the paper discussed how to make a certain family of algorithms run faster on parallel machines and the second half of the paper went on to show how, with some work, the method could be stretched to a wider class of algorithm. What the authors seemed to miss was that the extra work really wasn't necessary and the methods of the first half apply, with no change, to the second half. But don't take this as a criticism! I learnt a whole new way to approach algorithm design, and the trick to making the second half easy uses methods that have become popular in more recent years. Doing a web search I found lots of papers describing something similar to what I did.

This is also a nice example of how the notion of abstraction in computing and the notion of abstraction in mathematics are exactly the same thing. But I'm getting ahead of myself.

So here a direct translation from he paper of some procedural code to find the largest sum that can be made from a subsequence of a sequence. This will be our first implementation of the problem examined in the second half of the paper:


> f1 [] (sofar,max) = (sofar,max)
> f1 (b:bs) (sofar,max) =
>     let sofar' = if sofar+b<0
>             then 0
>             else sofar+b
>         max' = if max<sofar'
>             then sofar'
>             else max
>     in f1 bs (sofar',max')

sofar is a running sum that is reset each time it dips below zero, and max keeps track of the best sum so far. We initialise sofar and sum with 0 and -infinity. Here's an example of how to ue it:


> b :: [Double]
> b = [1..5] ++ [5,4..(-10)] ++ [(-2)..6]

> infinity :: Double
> infinity = 1/0
> test1 b = snd $ f1 b (0,-infinity)

Notice how we prime the algorithm with a starting vector. The 0 corresponds to the fact that at the start we've summed over 0 elements and the -infinity corresponds to the fact that we want the first sum computed to be the highest so far at that point.

Test the code with test1 b. We'll use a similar pattern all the way through this code.

You may see the problem with making this parallelisable: we are maintaining running sums so that the final values of sofar and max all depend on what was computed earlier. It's not obvious that we can break this up into pieces. sofar computes sums of subsequences between resets but chopping the array b into pieces might split such subsequences between processors. How can we handle this cleanly?

The first step is to write version two of the above function using max instead of conditionals:


> f2 [] (sofar,max) = (sofar,max)
> f2 (b:bs) (sofar,max) =
>     let sofar' = Prelude.max (sofar+b) 0
>         max'   = Prelude.max max sofar'
>     in f2 bs (sofar',max')

> test2 b = snd $ f2 b (0,-infinity)

But that doesn't appear to make things any easier, we've just buried the conditionals inside max, it doesn't make the serial dependency go away.

So let's solve another problem instead. In f2 I'll replace max with addition and addition with multiplication. 0 is the identity for addition so we should replace it with the identity for multiplication, 1. Similarly, -infinity is the identity for max so we should replace that with 0. We get:


> f3 [] (sofar,max) = (sofar,max)
> f3 (b:bs) (sofar,max) =
>     let sofar' = sofar*b+1
>         max'   = max+sofar'
>     in f3 bs (sofar',max')

> test3 b = snd $ f3 b (1,0)

That's all very well but (1) it looks no more parallelisable and (2) it's solving the wrong problem. Let's ignore problem (2) for now.

The thing that makes f3 easier to work with is that it's almost a linear function acting on the vector (sofar,max). Linear functions have one very nice property. If f and g are linear then we can compute f(g(x)) by acting with g first, and then applying f. But we can also compose f and g without reference to x giving us another linear function. We only have to know how f and g act on basis elements and we can immediately compute how f . g acts on basis elements. This is usually expressed by writing f and g as matrices. So let's tweak f3 so it's linear in its last argument:


> f4 [] (sofar,max,i) = (sofar,max,i)
> f4 (b:bs) (sofar,max,i) =
>     let sofar' = (sofar * b) + i
>         max'   = max + sofar'
>         i'     = i
>     in f4 bs (sofar',max',i')

> test4 b = let (_,max,_) = f4 b (1,0,1) in max

So now I need to write some code to work with linear functions. I'll do it in a very direct style. Here are some tuples representing basis vectors. (I could have written some fancy vector/matrix code but I don't want to distract from the problem in hand.)


> x,y,z :: Num a => (a,a,a)
> x = (1,0,0)
> y = (0,1,0)
> z = (0,0,1)

And here's some code that computes how a function acts on a basis, in effect finding the matrix for our function with respect to the basis x,y,z:


> matrix f = (f x,f y,f z)

Some simple operations on vectors:


> (a,b,c) .+ (d,e,f) = (a + d,b + e,c + f)
> a .* (b,c,d) = (a * b,a * c,a * d)

And now a little function that, given how f acts on basis elements, can apply f to any vector:


> apply (mx,my,mz) (sofar,max,i) = (sofar .* mx) .+ (max .* my) .+ (i .* mz)

Now we can redo the calculation with f4 by first making the matrix for f4, and then applying that to our starting vector.


> test5 b = let m = matrix (f4 b)
>               (_,max,_) = apply m (1,0,1)
>           in max

Note how by time we've computed m we've done almost all of the work even though the code hasn't yet touched (1,0,1).

But now comes the first bit of magic. We can split our list b into pieces. Compute the corresponding matrix for each piece on a separate processor, and then apply the resulting matrices to our starting vector.

Let's chop our list of reals into pieces of size n:


> chop n [] = []
> chop n l = let (a,b) = splitAt n l in a : chop n b

We'll use pieces of size 10:


> test6 b = max where
>    (_,max,_) = foldr ($) (1,0,1) (reverse b_functions) where
>    b_pieces = chop 10 b

The following maps should be replaced with a parallel version. It's easy to do this.


>    b_matrices = map (matrix . f4) b_pieces
>    b_functions = map apply b_matrices

Great, we've successfully parallelised the code, but it's the wrong algorithm. How can we use this to solve the correct problem? Remember how we replaced max with addition and addition with multiplication. We just have to undo that. That's all! Everything required to prove that the above parallelisation is valid applies over any semiring. At no point did we divide or subtract, and we only used elementary properties of numbers like a*(b+c) = a*b+a*c. That property holds for max and addition. In fact a+max(b,c) = max(a+b,a+c). We don't even have to modify the code. We can just define the max-plus semiring as follows:


> newtype MaxPlus = M Double deriving (Eq,Show)

> instance Num MaxPlus where
>     fromInteger 0 = M (-infinity)
>     fromInteger 1 = M 0
>     fromInteger _ = error "no conversion from general integer"
>     M a + M b = M (max a b)
>     M a * M b = M (a+b)
>     abs _ = error "no abs"
>     signum _ = error "no signum"
>     negate _ = error "no negate"

And now all we need is


> test7 b = test6 (fmap M b)

(I wonder if ghc is smart enough to completely eliminate that fmap M, after all, on a newtype, M should do zero work.)

And that's a completely parallelised version of the original algorithm.

There is a ton of optimisation that can be performed here. In particular, matrix applies a function to a fixed basis. For the particular function we're using here there's a big win from constant folding. The same constant folding applies in any semiring.

And back to the point I made at the beginning. By abstracting from the reals to a general semiring we are able to make the same code perform multiple duties: it can work on functions linear over many semirings, not just the reals. Mathematicians don't work with abstractions just to make life hell for students - they do so because working with abstract entities allows the same words to be reused in a valid way in many different contexts. This benefits both mathematicians and computer scientists.

Here's a link you can use to find out more on this technique.

Note that in real world usage you wouldn't use lists. -|chop|- would take longer than the rest of the algorithm.

PS A curious aside. I spent ages trying to get ghc to compile this code and getting my homebrew HTML markup code to work reliably on it. But eventually I located the problem. I've just returned from Hawai`i where I wrote most of this code. Just for fun I'd put my keyboard into Hawai`ian mode and forgot about it. When I did that, the single quote key started generating the unicode symbol for the Hawai`ian glottal stop. It looks just like a regular single quote in my current terminal font so it was hard to see anything wrong with the code. But, quite reasonably, ghc and many other tools barf if you try to use one where a regular quote is expected.

Some thoughts on reasoning and monads

2008-11-15T12:52:00.000-08:00

Just a short note about some half-formed thoughts on the subject of monads and reasoning.

Haskell types and total functions, at least for a suitable subset of Haskell, form a category, with types and functions as the objects and arrows. Haskell code is usually full of expressions whose types don't correspond to arrows (e.g.. 1::Int). But in category theory we can only really talk of objects and arrows, not elements of objects. Many categories are made up of objects that have no reasonable notion of an element. By writing code in pointfree style we can eliminate reference to individual elements and talk only about arrows.

For example, suppose f :: A -> B and g :: B -> C. Then we can define h :: A -> C in a couple of different ways. Pointfree:


h = g . f

or pointful:


h x = let y = f x
          z = g y
      in z

For examples as simple as this it's not hard to write code that will translate between one style and the other. But the former definition has an advantage, it works in any category. We'll consider a very simple example. Let R be the category of real numbers, where there is an arrow from x to y if x≤y. We can think of f and g as arrows in this category. If we know A≤B and B≤C then the first definition above tells us how to construct an arrow from A to C and hence it proves that A≤C. The second definition makes no sense because it relies on the notion of x as an element of A and uses the functions f and g to construct elements of B and C. These words hold no meaning in the category R where the arrows aren't functions and the objects aren't containers of elements.

Except that's not quite true.

Because there is a scheme for translating the second pointful definition to pointfree form, the second definition does in fact provide a proof that A≤C. We just have to bear in mind that the proof needs to be translated into pointfree form first. In fact, we can happily spend our day using pointful style to generate proofs about R, as long as at the end of the day we translate our proofs to pointfree notation. In fact, Haskell programmers know that it's often much easier to write programs in pointful style so it seems reasonable to guess that there are many proofs that are easier to write in pointful style even though they can't be interpreted literally. Philosophically this is a bit weird. As long as we restrict ourselves to chains of reasoning that can be translated, we can use intuitions about elements of objects to make valid deductions about domains where these notions make no sense.

Part I of Lambek and Scott is about the correct pointful language to use when talking about cartesian closed categories (CCCs). They use a form of typed lambda calculus. Every arrow in a cartesian category can be written as a pointful lambda expression. Even though there isn't a meaningful way to assign meaning to the 'points' individually, every lambda abstraction gives an arrow in a cartesian closed category that can be built using the standard parts that come with a CCC, and vice versa. The pointful definition of h above is an example. Here's a (very) slightly less trivial example. Given f :: A -> C and g :: B -> D we can define h :: (A,C) -> Either B D as


f x = let y = fst x
          z = Left y
      in z

Again, that code can be translated into pointfree form using only the standard parts available in a cartesian closed category. Actually, we only need a category with products and coproducts for this example. So consider a lattice. The above code can be translated into a proof that A∩C≤B∪D. (I'm using cap and cup for join and meet.) Again, it makes no sense to talk of elements of a general lattice as containers of elements, and yet after translation to pointfree notation we have a valid proof. If you can restrict yourself to a suitable subset of set theory then your proofs involving elements of objects can carry over to many other categories.

Part II of Lambek and Scott is about taking this idea to extremes. It's about categories known as toposes. A topos is a category that is even more like the category of sets than a CCC. It's still general enough that there are many kinds of toposes, but you can use a sizable portion of first order logic and ZF to make deductions about them. Again, the literal notion of membership of the objects of a topos might make no sense, but the proofs have a translation to pointfree notation. In fact, it's possible to write entire papers in what looks like conventional set theory language, and have them be valid for other toposes. Anders Kock, for example, writes such papers. Chris Isham has been arguing that topos theory is the correct framework for physics. If you interpret your propositions as being in the category Set then you get classical physics. But those same propositions can be interpreted in other categories, such as one for quantum mechanics, giving a way to use and extend classical language to reason about quantum systems. This set theory-like language is known as the "internal language" of a topos.

Anyway, I'm interested in the notion that Haskell do-notation provides another kind of pointful language that can be used to reason about situations where points don't seem at first to make sense. Consider the lattice of subsets of a topological space, ordered by inclusion. Let Cl(U) be the closure of U. Cl satisfies these properties:

X⊂Y implies Cl(X)⊂Cl(Y)
X⊂Cl(X)
and
Cl(CL(X))⊂Cl(X)

Look familiar? If we make this lattice into a category, with arrows being inclusions, then the first property states that Cl is a functor and the next two say that Cl is a monad. In fact, monads are a kind of generalised closure. So now suppose we're given A⊂Cl(B) and B⊂Cl(C) and wish to prove that A⊂Cl(C). We can rephrase this by saying that if f :: A -> Cl B and g :: B -> Cl C we need an arrow h :: A -> Cl C. We can write one like this:


h x = do
    y <- f x
    z <- g y
    return z

Now it's tempting to interpret the inclusions as functions with y <- f x saying that y is the image of x under the inclusion. (I don’t know about you, but when I write x <- getChar I think to myself “x is the return value from calling getChar”.) But that doesn't really work because the type of f x is Cl B but y is of type B. On the other hand, we can radically reinterpret the above as something like this: when arguing about chains of inclusions of subsets of a topological space, as long as at the end of the chain of inclusions you always end up in the closure of some subset, you're allowed to cheat and nudge a generic point in the closure of a subset back into the original subset. This is exactly parallel to the way do-notation seems to allow us to extract elements out of monadic objects as long as at the end of the do-block we always return an element of monadic type. I'm sure that with a bit of work we could produce a rigorous metatheorem from this. I also expect we can also produce something similar for comonads.

Anyway, the moral is that when working with categories with monads there are may be some interesting and unusual ways to reason. The example of the lattice of subsets is fairly trivial but I'm sure there are other interesting examples. I also expect there's a nice connection with modal logic. I now think of Haskell do-notation as the "internal language" of a category with monads.

Update: I left out the crucial sentences I meant to write. It's easy to see do-notation as a kind of Haskell specific trick for making things like IO heavy code look like traditional procedural code. But comparison with the theory in Lambek and Scott, and Topos Theory in general, makes it clear that do-notation is a member of a family of related languages.

From Monoids to Monads

2008-11-08T10:13:00.000-08:00

Generalising Monoids

The word 'monad' is derived from the word 'monoid'. The explanation usually given is that there is an analogy between monoids and monads. On the surface, this seems a bit unlikely. The join operation in a monad is supposed to correspond to the binary operator in the monoid, but join is a completely different kind of thing, certainly not a binary operator in any usual sense.

I'm going to make this analogy precise so that it's clear that both monoids and monads are examples of the same construction. In fact, I'm going to write some Haskell code to define monoids and monads in almost exactly the same way. I was surprised to find I could do this because instances of Haskell's Monoid and Monad aren't even the same kind of thing (where I'm using 'kind' in its technical sense). But it can be done.

So let's start thinking about monoids. They are traditionally sets equipped with a special element and a binary operator so that the special element acts as an identity for the binary operator, and where the binary operator is associative. We expect type signatures something like one :: m and mult :: m -> m -> m so that, for example, m (m a b) c == m a (m b c). That's fine as it stands, but it doesn't generalise easily. In particular it'd be nice to generalise this definition to other categories. To do that we need to rephrase the definitions so that they are completely point-free and are written purely as the compositions of arrows.

Let's start by thinking about the rule that says multiplying by the identity on the left should leave a value unchanged. Ie. we want mult one x == x for all x. We already have a problem, it refers to a couple of 'points', both one, the identity, and x the unknown. We can deal with the first one easily, we just replace one with an arrow () -> m. But we also need some plumbing to provide two arguments to the mult function. Rather than belabour the point, I'll just give the full code:


> {-# OPTIONS_GHC -fglasgow-exts #-}

> import Control.Monad
> import Test.QuickCheck

> class Monoid m where
>     one :: () -> m
>     mult :: (m,m) -> m

The law for multiplication on the left is then given by law1_left == law1_middle and so on:


>     law1_left,law1_middle,law1_right :: m -> m
>     law1_left   = mult . (one <#> id) . lambda
>     law1_middle = id
>     law1_right  = mult . (id <#> one) . rho

The associativity law is then given by law2_left x = law2_right x.


>     law2_left,law2_right :: ((m,m),m) -> m
>     law2_left = mult . (mult <#> id)
>     law2_right = mult . (id <#> mult) . alpha

The left and right hand sides of the laws are now point-free. But in order to do this I've had to write some auxiliary functions:


> lambda :: a -> ((),a)
> lambda x = ((),x)

> rho :: a -> (a,())
> rho x = (x,())

> alpha :: ((a,b),c) -> (a,(b,c))
> alpha ((x,y),z) = (x,(y,z))

I've also used the fact that (,) is a bifunctor, ie. it's a functor in each of its arguments so (,) doesn't just give a way to generate a new type (a,b) from types a and b. It also combines pairs of arrows to make new arrows. I'll call the part of (,) that acts on arrows by the name <#>:


> (<#>) :: (a -> c) -> (b -> d) -> (a, b) -> (c, d)
> (f <#> g) (x,y) = (f x,g y)

Intuitively, f <#> g maps f on the left and g on the right of (a,b)

Try unpacking those definitions to see that we get the usual monoid laws. For example


law2_left ((x,y),z) == mult $ (mult <#> id) ((x,y),z) == mult (mult (x,y),z)
law2_right ((x,y),z) == mult $ (id <#> mult) $ alpha ((x,y),z) == mult (x,mult (y,z))

So we get


mult (mult (x,y),z) ==  mult (x,mult (y,z))

the usual associativity law.

Now that this definition is point-free it seems we could just carry it over to any category. In fact, we've implicitly done this because we've carried over the definition of a monoid from Set to Hask, the category of Haskell types and total functions. We're so used to treating Hask as a proxy for Set we can forget they are actually different categories. But this definition of a monoid works in both. But what about that lambda, rho and alpha? Well they're easy to define in any Cartesian closed category (CCC). But we don't need all of the features of a CCC to define a monoid, we just need lambda, rho and alpha and some kind of 'product' on the set of objects that also acts like a bifunctor. In fact, there's a bunch of 'obvious' laws that these functions satisfy in a CCC. Any category with these functions satisfying these same laws is called a monoidal category. The above definitions allow us to transfer the definition of a monoid to any such category. For the full definition, see the wikipedia entry.

Anyway, let's check to see if the type Int might be a monoid:


> instance Monoid Int where
>     one _ = 1
>     mult (a,b) = a*b

> check1 = quickCheck $ \n -> law1_left n == law1_middle (n :: Int)
> check2 = quickCheck $ \n -> law1_left n == law1_right (n :: Int)
> check3 = quickCheck $ \n -> law2_left n == (law2_right n :: Int)

Of course that's no proof, but it gives us some confidence. (Eg. what about large numbers close to 2³¹...?)

Another Category in the World of Haskell

Most people around here will be familiar with how Haskell types and functions form a category Hask in the obvious way. But it's less well known that there is another category lurking in Haskell. Consider the set of all Haskell functors. These are endofunctors, ie. functors Hask→Hask. Between any two functors is the set of natural transformations between those functors. (If f and g are functors, then the polymorphic functions f a -> g a are the natural transformations.) In addition, you can compose natural transformations and the composition is associative. In other words, Haskell functors form a category. (See the appendix for more abstract nonsense relating to this.)

We'll call the category of endofunctors on Hask by the name E. Functors can be composed like so:


> type (f :<*> g) x = f (g x)

It'd be cool if this was a product in the usual categorical sense, but it isn't. There isn't a natural way to map to both f a and g a from f (g a) with the universal property of products. Instead it's a weaker type of product which is still a bifunctor. Here's how it acts on arrows (remembering that in E, arrows are natural transformations):


> (<*>) f g = f. fmap g

(We could also have used fmap g . f.) If you think of f <#> g as making f act on the left and g act on the right, then you can think of f <*> g as making f and g act on the outer and inner containers of a container of containers.

Here's the identity element for this product, the identity functor. It plays a role similar to () in Hask:


> data Id x = Id x deriving Show
> instance Functor Id where
>     fmap f (Id x) = Id (f x)

We can now define some familiar looking natural transformations:


> lambda' :: Functor f => f a -> (Id :<*> f) a
> lambda' x = Id x

> rho' :: Functor f => f a -> (f :<*> Id) a
> rho' x = fmap Id x

> alpha' :: f (g (h a)) -> f (g (h a))
> alpha' = id

With these we have turned E into a monoidal category.

(OK, this may be confusing. We're 'multiplying' functors and we have associativity and a left- and right-identity. So functors form a monoid (modulo isomorphism). But that's a distraction, these are not the monoids that you're looking for. See the appendix for more on this.)

So now we can define monoids in this category using almost the same code as above:


> class Functor m => Monoid' m where
>     one' :: Id a -> m a
>     mult' :: (m :<*> m) a -> m a

>     law1_left',law1_middle',law1_right' :: m a -> m a
>     law1_left'   = mult' . (one' <*> id) . lambda'
>     law1_middle' = id
>     law1_right'  = mult' . (id <*> one') . rho'

>     law2_left',law2_right' :: ((m :<*> m) :<*> m) a -> m a
>     law2_left'  = mult' . (mult' <*> id)
>     law2_right' = mult' . (id <*> mult') . alpha'

And here's the punch line. That's precisely a monad, laws 'n' all.

If you want, you can unpack the definitions above and see that they correspond to the usual notion of a monad. We can write code to do the translation:


> data Monad m => TranslateMonad m a = TM { unTM :: m a } deriving (Eq,Show)

> translate :: Monad m => m a -> TranslateMonad m a
> translate x = TM x

> instance (Monad m,Functor m) => Functor (TranslateMonad m) where
>     fmap f (TM x) = TM (fmap f x)

> instance (Functor m,Monad m) => Monoid' (TranslateMonad m) where
>     one' (Id x) = TM $ return x
>     mult' (TM x) = TM $ fmap unTM x >>= id

In other words, any instance of Monad gives an instance of Monoid'. I'll let you write the reverse mapping. We can even check (not prove!) the laws are satisfied by a particular monad by using QuickCheck:


> instance Arbitrary a => Arbitrary (Id a) where
>     arbitrary = liftM Id arbitrary

> instance (Monad m,Eq (m a),Arbitrary (m a)) => Arbitrary (TranslateMonad m a) where
>     arbitrary = liftM TM arbitrary

> check4 = quickCheck $ \n -> law1_left' n == law1_middle' (n :: TranslateMonad [] Int)
> check5 = quickCheck $ \n -> law1_left' n == law1_right' (n :: TranslateMonad [] Int)
> check6 = quickCheck $ \n -> law2_left' n == (law2_right' n :: TranslateMonad [] Int)

I don't know about you, but I find this mind-blowing. We took a definition of a simple concept, the monoid. We then abstractificated its definition so that it applied to any monoidal category, and in doing so we completely changed the meaning of one and mult. And yet the resulting object, including its laws, are precisely what you need to solve a whole host of problems in algebra and Haskell programming. This is an amazing example of the unreasonable effectiveness of mathematics. The concept might be a little tricky to grasp: monads are like ordinary monoids, but with outer and inner replacing left and right. But the payoff from this is that intuitions about monoids carry over to monads. I hope to say more about this in future episodes.

Appendix

Consider the category, H, with one object, Hask, and whose arrows are the endofunctors on Hask. We've shown how the arrows on this category aren't just a set, they themselves form a category. This makes H a 2-category. A category with one object is a monoid. But this is a 2-category, so we have a 2-monoid. In fact, there are some extra details required to show we have a 2-category. For example, we need to show that composition of arrows (which now form a category, not a set) is a functor (not just a function). That follows from the Interchange Law which I've already talked about. But note that H is a monoid in a completely conventional way, its arrows form a set with a binary operator on them. This is not the structure that corresponds to monads, although it plays a part in constructing them. Also, don't confuse this monoid with the one that appears in a MonadPlus.

Confession

I'm having trouble giving <*> the correct type signature. I think it should be something like (forall x.a x -> c x) -> (forall x.b x -> d x) -> (forall x.a (b x) -> c (d x)). GHC doesn't like it. Can anyone come up with the correct thing? The code still works.

Links

Another account of the same subject, modulo some details.

Operads and their Monads

2008-10-25T18:45:00.000-07:00

Hardly a day goes by at the n-Category Cafe without operads being mentioned. So it's time to write some code illustrating them.


> {-# LANGUAGE FlexibleInstances #-}

> import Data.Monoid
> import Control.Monad.Writer
> import Control.Monad.State

Let's define a simple tree type:


> data Tree a = Leaf a | Tree [Tree a] deriving (Eq,Show)

Sometimes we want to apply a function to every element of the tree. That's provided by the fmap member of the Functor type class.


> instance Functor Tree where
>     fmap f (Leaf x) = Leaf (f x)
>     fmap f (Tree ts) = Tree $ map (fmap f) ts

But just as we can't use map to apply monadic functions to a list (we'd write mapM print [1,2,3]), we can't use fmap to apply them to our tree. What we need is a monadic version of fmap. Here's a suitable type class:


> class FunctorM c where
>     fmapM :: Monad m => (a -> m b) -> c a -> m (c b)

(I could have used Data.Traversable but that entails Data.Foldable and that's too much work.)

And now we can implement a monadic version of Tree's Functor instance:


> instance FunctorM Tree where
>     fmapM f (Leaf x) = do
>         y <- f x
>         return (Leaf y)
>     fmapM f (Tree ts) = do
>         ts' <- mapM (fmapM f) ts
>         return (Tree ts')

We can use fmapM to extract the list of elements of a container type:


> serialise a = runWriter $ fmapM putElement a
>               where putElement x = tell [x]

Not only does serialise suck out the elements of a container, it also spits out an empty husk in which all of the elements have been replaced by (). We can think of the latter as the 'shape' of the original structure with the original elements removed. We can formalise this as


> type Shape t = t ()

So we have:


> serialise :: FunctorM t => t a -> (Shape t,[a])

Keeping the shape around allows is to invert serialise to give:


> deserialise :: FunctorM t => Shape t -> [a] -> (t a, [a])
> deserialise t = runState (fmapM getElement t) where
>   getElement () = do
>       x:xs <- get
>       put xs
>       return x

(That's a bit like using the supply monad. This function also returns the leftovers.)

We can also write (apologies for the slightly cryptic use of the writer monad):


> size :: FunctorM t => t a -> Int
> size a = getSum $ execWriter $ fmapM incCounter a
>          where incCounter _ = tell (Sum 1)

Let's try an example. Here's an empty tree:


> tree1 = Tree [Tree [Leaf (),Leaf ()],Leaf ()]

We can pack a bunch of integers into it:


> ex1 = fst $ deserialise tree1 [1,2,3]

And get them back out again:


> ex2 = serialise ex1

serialise separates the shape from the data, something you can read lots more about at Barry Jay's web site.

Remember that trees are also monads:


> instance Monad Tree where
>     return x = Leaf x
>     t >>= f = join (fmap f t) where
>         join (Leaf t) = t
>         join (Tree ts) = Tree (fmap join ts)

The join operation for a tree grafts the elements of a tree of trees back into the tree.

Right, that's enough about trees and shapes for now.

Operads are a bit like the plumbing involved in installing a sprinkler system. Suppose you have a piece, A that splits a single pipe into two:


       ______
      /
     /   ____
____/   /
     A /
____   \
    \   \____
     \
      \______

And you have two more pipes B and C that look like this:


        ______
       /
      /   ____
     /   /
____/    \____
          
    B or C
____      ____
    \    /
     \   \____
      \
       \______

then you can 'compose' them to make a larger system like this:


                       ______
                      /
                     /   ____
                    /   /
          _________/    \____
         /            B  
        /   _______      ____
       /   /       \    /
      /   /         \   \____
     /   /           \
____/   /             \______
     A /
____   \               ______
    \   \             /
     \   \           /   ____
      \   \         /   /
       \   \_______/    \____
        \             C  
         \_________      ____
                   \    /
                    \   \____
                     \
                      \______

(Vim rectangular mode FTW!)

The important thing to note here is that as A had two outputs (or inputs, depending on your point of view) you can attach two more pieces, like B and C, directly to it.

Call the number of outlets the 'degree' of the system. If A has degree n then we can attach n more systems, A₁...A_n to it and the resulting system will have degree degree(A₁)+...+degree(A_n). We can write the result as A(A₁,...,A_n).

We also have the 'identity' pipe that looks like this:


_____________

  identity
_____________

Formally, an operad is a collection of objects, each of which has a 'degree' that's an integer n, n≥0, depending on your application), an identity element of degree 1, and a composition law:


> class Operad a where
>     degree :: a -> Int
>     identity :: a
>     o :: a -> [a] -> a

o is the composition operation. If f has degree n then we can apply it to a list of n more objects. So we only expect to evaluate f `o` fs successfully if degree f == length fs.

There are many identities we'd expect to hold. For example f `o` [identities,...,identity] == f, because adding plain sections of pipe has no effect. We also expect some associativity conditions coming from the fact that it doesn't matter what order we build a pipe assembly in, it'll still function the same way.

We can follow this pipe metaphor quite closely to define what I think of as the prototype Operad. A Fn a is a function that takes n inputs of type a and returns one of type a. As we can't easily introspect and find out how many arguments such a function expects, we also store the degree of the function with the function:


> data Fn a = F { deg::Int, fn::[a] -> a }

> instance Show a => Show (Fn a) where
>     show (F n _) = "<degree " ++ show n ++ " function>"

unconcat is a kind of inverse to concat. You give a list of integers and it chops up a list into pieces with lengths corresponding to your integers. We use this to unpack the arguments to f `o` gs into pieces suitable for the elements of gs to consume.


> unconcat [] [] = []
> unconcat (n:ns) xs = take n xs : unconcat ns (drop n xs)

> instance Operad (Fn a) where
>     degree = deg
>     f `o` gs = let n = sum (map degree gs)
>                in F n (fn f . zipWith fn gs . unconcat (map degree gs))
>     identity = F 1 head

Now compute an example, f(f₁,f₂,f₃) applied to [1,2,3]. It should give 1+1+2*3=8.


> ex3 = fn (f `o`  [f1,f2]) [1,2,3] where
>   f = F 2 (\[x,y] -> x+y)
>   f1 = F 1 (\[x] -> x+1)
>   f2 = F 2 (\[x,y] -> x*y)

(That's a lot like lambda calculus without names. Operads are a bit like n-ary combinators.)

Now I'm going to introduce a different looking operad. Think of V as representing schemes for dicing the real line. Here are some examples:


|0.25|0.25|   0.5   |

|0.1|0.1|    0.8    |

If A divides up the real line into n pieces then you could divide up each of the n pieces using their own schemes. This means that dicing schemes compose. So if we define A, B and C as:


A = |0.5|0.5|

B = |0.75|0.25|
C = |0.1|0.1|0.8|

Then A(B,C) is:


|0.375|0.125|0.05|0.05|0.4|

We could implement V as a list of real numbers, but it's more fun to generalise to any monoid and not worry about divisions summing to 1:


> data V m = V { unV :: [m] } deriving (Eq,Show)

This becomes an operad by allowing the monoid value in a 'parent' scheme multiply the values in a 'child'.


> instance Monoid m => Operad (V m) where
>     degree (V ps) = length ps
>     (V as) `o` bs = V $ op as (map unV bs) where
>         op [] [] = []
>         op (a:as) (b:bs) = map (mappend a) b ++ op as bs
>     identity = V [mempty]

For example, if d₁ cuts the real line in half, and d₂ cuts it into thirds, then d₁(d₁,d₂) will cut it into five pieces of lengths 1/4,1/4,1/6,1/6,1/6:


> ex4 = d1 `o` [d1,d2] where
>   d1 = V [Product (1/2),Product (1/2)]
>   d2 = V [Product (1/3),Product (1/3),Product (1/3)]

If the elements in a V are non-negative and sum to 1 we can think of them as probability distributions. The composition A(A₁,...,A_n) is the distribution of all possible outcomes you can get by selecting a value i in the range {1..n} using distribution A and then selecting a second
value conditionally from distribution A_i. We connect with the recent n-category post on entropy.

In fact we can compute the entropy of a distrbution as follows:


> h (V ps) = - (sum $ map (\(Product x) -> xlogx x) ps) where
>   xlogx 0 = 0
>   xlogx x = x*log x/log 2

We can now look at the 'aside' in that post. From an element of V we can produce a function that computes a corresponding linear combination (at least for Num types):


> linear (V ps) xs = sum $ zipWith (*) (map (\(Product x) -> x) ps) xs

We can now compute the entropy of a distribution in two different ways:


> (ex5,ex6) = (h (d1 `o` [d1,d2]),h d1 + linear d1 (map h [d1,d2])) where
>   d1 = V [Product 0.5,Product 0.5]
>   d2 = V [Product 0.25,Product 0.75]

Now according to this paper on operads we can build a monad from an operad. Here's the construction:


> data MonadWrapper op a = M { shape::op, value::[a] } deriving (Eq,Show)

(The field names aren't from the paper but they do give away what's actually going on...)

The idea is that an element of this construction consists of an element of the operad of degree n, and an n element list. It's a functor in an obvious way:


> instance Functor (MonadWrapper o) where
>     fmap f (M o xs) = M o (map f xs)

It's also a FunctorM:


> instance FunctorM (MonadWrapper o) where
>   fmapM f (M s c) = do
>       c' <- mapM f c
>       return $ M s c'

We can make the construction a monad as follows:


> instance Operad o => Monad (MonadWrapper o) where
>     return x = M identity [x]
>     p >>= f = join (fmap f p) where
>         join (M p xs) = M (p `o` map shape xs) (concatMap value xs)

Now for something to be a monad there are various laws that needs to be satisfied. These follow from the rules (which I haven't explicitly stated) for an operad. When I first looked at that paper I was confused - it seemed that the operad part and the list part didn't interact with each other. And then I suddenly realised what was happening. But hang on for a moment...

Tree shapes make nice operads. The composition rule just grafts child trees into the leaves of the parent:


> instance Operad (Tree ()) where
>     degree t = length (snd (serialise t))
>     identity = Leaf ()
>     t `o` ts = let (r,[]) = deserialise t ts in r >>= id

We can write that more generically so it works with more than trees:

> data OperadWrapper m = O { unO::Shape m }


> instance (FunctorM m,Monad m) => Operad (OperadWrapper m) where
>     degree (O t) = size t
>     identity = O (return ())
>     (O t) `o` ts = let (r,[]) = deserialise t (map unO ts) in O (r >>= id)

So let's use the construction above to make a monad. But what actually is this monad? Each element is a pair with (1) a tree shape of degree n and (2) an n-element list. In other words, it's just a serialised tree. We can define these isomorphisms to make that clearer:


> iso1 :: FunctorM t => t x -> MonadWrapper (t ()) x
> iso1 t = uncurry M (serialise t)

> iso2 :: FunctorM t => MonadWrapper (t ()) x -> t x
> iso2 (M shape contents) = let (tree,[]) = deserialise shape contents in tree

So, for example:


> ex7 = iso2 (iso1 tree) where
>   tree = Tree [Tree [Leaf "Birch",Leaf "Oak"],Leaf "Cypress",Leaf "Binary"]

That construction won't work for all monads, just those monads that come from operads. I'll leave you to characterise those.

And now we have it: a way to think about operads from a computational perspective. They're the shapes of certain monadic serialisable containers. Operadic composition is the just the same grafting operation used in the join operation, using values () as the graft points.

I have a few moments spare so let's actually do something with an operad. First we need the notion of a free operad. This is basically just a set of 'pipe' parts that we can stick together with no equations holding apart from those inherent in the definition of an operad. This is different from the V operad where many different ways of apply the o operator can result in the same result. We can use any set of parts, as long as we can associate an integer with each part:


> class Graded a where
>   grade :: a -> Int

The free operad structure is just a tree:


> data FreeOperad a = I | B a [FreeOperad a] deriving Show

I will be the identity, but it will also serve as a 'terminator' like the way there's always a [] at the end of a list.

An easy way to make a single part an element of an operad:


> b n = let d = grade n in B n (replicate d I)

Here's the instance:


> instance Graded a => Operad (FreeOperad a) where
>   degree I = 1
>   degree (B _ xs) = sum (map degree xs)
>   identity = I
>   I `o` [x] = x
>   B a bs `o` xs = let arities = map degree bs
>                   in B a $ zipWith o bs (unconcat arities xs)

Now I'm going to use this to make an operad and then a monad:


> instance Graded [a] where
>   grade = length

> type DecisionTree = MonadWrapper (FreeOperad [Float])

What we get is a lot like the probability monad except it doesn't give the final probabilities. Instead, it gives the actual tree of possibilities. (I must also point out this hpaste by wli.)


> test = do
>   a <- M (b [Product 0.5,Product 0.5]) [1,2]
>   b <- M (b [Product (1/3.0),Product (1/3.0),Product (1/3.0)]) [1,2,3]
>   return $ a+b

Now we can 'flattten' this tree so that the leaves have the final probabilities:


> flatten :: Monoid m => FreeOperad [m] -> V m
> flatten I = V [mempty]
> flatten (B ms fs) = V $ concat $ zipWith (map . mappend) ms (map (unV . flatten) fs)

This is a morphism of operads. (You can probably guess the definition of such a thing.) It induces a morphism of monads:


> liftOp :: (Operad a,Operad b) => (a -> b) -> MonadWrapper a x -> MonadWrapper b x
> liftOp f (M shape values) = M (f shape) values

liftOp flatten test will give you the final probabilities.

There may just be a possible application of this stuff. The point of separating shape from data is performance. You can store all of your data in flat arrays and do most of your work there. It means you can write fast tight loops and only rebuild the original datatype if needed at the end. If you're lucky you can precompute the shape of the result, allowing you to preallocate a suitable chunk of memory for your final answer to go into. What the operad does is allow you to extend this idea to monadic computations, for suitable monads. If the 'shape' of the computation is independent of the details of the computation, you can use an operad to compute that shape, and then compute the contents of the corresponding array separately. If you look at the instance for MonadWrapper you'll see that the part of the computation that deals with the data is simply a concatMap.

BTW In some papers the definition restricts the degree to ≥1. But that's less convenient for computer science applications. If it really bothers you then you can limit yourself to thinking about containers that contain at least one element.

What's the use of a transfinite ordinal?

2008-10-18T15:14:00.000-07:00

This post is going to talk about one player games (which I'll call 1PGs, or just games). They seem pretty stupid at first, but it turns out that they have a lot of rich and applicable theory. The idea is that at any position in a 1PG there is a set of new positions that you can move to. You lose when you have no legal moves and the goal is to avoid losing for as long as you can. For convenience, let's identify your position in a 1PG with the set of positions you could move to. A position in a 1PG is just another 1PG. So we can give the following recursive formal definition of a 1PG:

A one player game is either the losing position with no legal moves, denoted by ∅, or a set of one player games.

A sequence of legal moves in a 1PG, say x, is then a sequence of the form x ∋x₀ ∋x₁ ∋x₂ ∋....

Did you notice that I sneaked in a really big restriction on 1PGs in that definition? My definition means that 1PGs are the same thing as sets. And sets are well-founded meaning that any sequence of the form x ∋x_0 ∋x₁ ∋x₂ ∋... eventually terminates. So my definition implicitly contains the restriction that you always eventually lose. That's fine for my purposes.

Now suppose that x ∋y ∋z is a legal sequence of plays in some 1PG. You might as well assume that z is in x. The reason is that if z were a legal move from x, you wouldn't take it, because in order to delay the inevitable you're going to prefer to go via y. So we'll assume that any game x is "transitively closed", ie. if x∋y and y∋z then x∋z.

Now I want to restrict things even further to the most boring kind of 1PG of all, games where the successor moves are all totally ordered. Intuitively what I mean is that you can think of every move in the game as simply moving you along a 'line'. Given any two positions in the game, x and y, either x∋y, or y∈x, or x=y. In these games there are no branches, the only thing you can do in these games is delay the inevitable. A bit like life really. We'll also use the notation a<b to mean a∈b.

I'll call these games ordinal 1PGs. You may have noticed that this is precisely the same as the ordinary definition of an ordinal.

Let's call the game ∅ by the name 0. 0 is sudden death. Recursively define the game n={0,1,2,3,n-1}. In the game of n you have up to n steps to play before you die (though you could hurry things along if you wanted to).

But what about this game: ω={0,1,2,...}. There's nothing mysterious about it. In your first move you choose an integer, and then you get to delay the end for up to n moves. This is our first transfinite ordinal, but it just describes a 1PG that anyone can play and which is guaranteed to terminate after a finite number of moves. We can't say in advance how many moves it will take, and we can't even put an upper bound on the game length, but we do know it's finite.

Given two games, a and b, we can define the sum a+b. This is simply the game where you play through b first, and then play through a. Because of the transitivity condition above you could simply abandon b at any stage and then jump into a, but that's not a good thing to do if you're trying to stay alive for as long as possible. The first example of such a 1PG is the game ω+1. You get to make one move, then you choose an n, and then you get to play n.

What about 1+ω? Well you get to choose n, play n moves, and then play an extra move. But that's exactly like playing ω and making the move n+1. So 1+ω and ω are the same game. Clearly addition of 1PGs isn't commutative.

We can multiply two games as well. Given games a and b, ab is defined like this: we have game a and game b in front of us at the same time, a on the left and b on the right. At any turn we have a choice: make a move in the game on the left, or make a move in the game on the right, resetting the left game back to a again. It's like playing through b, but playing through a copy of a at each turn. Note how for integers a and b the game ab is just ab, where the former is game multiplication, and the latter is ordinary arithmetical multiplication.

We can also define exponentiation. We can define aⁿ, for finite integer n, in the obvious way as ω·ω·...·ω n times. Before getting onto more general exponentiation I need to classify ordinal games. There are three kinds. There's

(1) the zero game 0=∅.
(2) the games where there is one 'largest' next move in the game, the one you should take. These are games of the form a+1 and the next position in the game, if you're playing to survive, will be a. For example, the game 7 is 6+1.
(3) the games where there is no single best choice of next move. These are games where there is no largest next move, but instead an ever increasing sequence of possible next moves and you have to pick one. The first example is ω where you're forced to pick an n. These are called limit ordinals.

We can use this to recursively define a^b using the classification of b.

(1) If b=0 then by definition a^b=1 and that's that.
(2) if b=c+1, then a^b=a^c·a.
(3) if b is a limit ordinal then a move in a^b works like this: you pick a move in b, say c, and now you get to play a^c (or, by transitive closure, any smaller game).

For example, you play ω^ω like this: in your first move you pick an n, and then you play ωⁿ. Note how the state of play never becomes infinite, you're always playing with finite rules on a finite 'board'. This is an important point. Even though ω^ω3+22+ω+4 seems like it must be a very big kind of thing, it describes a finite process that always terminates.

So what use are these things?

Suppose you have an algorithm that runs one step at a time, with each step taking a finite time, and where the internal state of the program at step t is s_t. Suppose that you can find a map f from program states to the ordinals such that f(s_t+1)<f(s_t), ie. the program always decreases the ordinal associated to a state. Then the program must eventually terminate.

For example, suppose the program iterates 10 times. We can simply assign f(s_t)=10-t. Suppose instead your algorithm computes some number n at the first step, though you can't figure out what that number is (or maybe it accepts n as an input from a user), then we can use f(0)=ω and so on. We don't know in advance what n is, but that doesn't matter. Whatever n is entered, the program will terminate.

This technique is particularly useful for proving termination of rewrite rules. For example, with the rewrite rules generated by the Buchberger algorithm we can map the first step to ωⁿ, for some n. And if we don't know what n is, we can just start with ω^ω. If you take a peek at the papers of Nachum Dershowitz you'll see how he applies this technique to many other types of rewrite rule.

Transfinite ordinals are not as esoteric as they may first appear to be. Using them to prove termination goes back to Turing and Floyd but I learnt about the method from reading Dershowitz. (I think Floyd suggested using the method to prove termination of the standard rewrite rules for computing derivatives of expressions. It's tricky because differentiating a large product, say, can result in many new terms, so there's a race between the derivative 'knocking down' powers, and products causing growth.)

Probably the best example of using ordinals to prove the termination of a process is in the Hydra game. What's amazing about this game is that you need quite powerful axioms of arithmetic to prove termination.

One more thing. All of this theory generalises to two player games.

Update: I think mjd may have been about to write on the same thing. On the other hand, though he mentions program termination, he says of it "But none of that is what I was planning to discuss". So maybe what I've written can be seen as complementary.