A Neighborhood of Infinity

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.

Untangling with Continued Fractions: part 5

2008-10-12T09:15:00.000-07:00


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

> module Main where

> import qualified Data.Map as M
> import Control.Monad
> import Prelude
> import Ratio hiding (numerator,denominator)

> infixl 5 .+
> infixl 6 .*

So far I've described a way to use monad notation to represent knots, links and tangles, and I've provided a way to implement a monad so that the notation gives rise to the bracket polynomial when applied to a knot. What happens when we do the same with a tangle?

A tangle corresponds to something with this type:


> type Tangle = (Bool,Bool) -> KnotMonad (Bool,Bool)

Here are two examples of tangles, ie. expressions with this type:


> infinity (i,j) = return (i,j)
> zero (i,j) = do
>   cup (i,j)
>   cap ()

These correspond to these tangles:

And the names of those expressions are the corresponding extended rational values.

Now here's a surprising but nontrivial fact that holds when T is a rational tangle: if we set A²=i (where i²=-1) then the value of the monadic expression corresponding to the tangle is of the form a<[∞]>+b<[0]> where the [r] are the rational tangles above, and <T> is the bracket polynomial (evaluated at A²=i). The rational associated to the original tangle is then given by -ia/b. I believe this discovery is, like many other things I've talked about in this blog, due to Conway.

We now have the makings of a complete algorithm. Using the polynomial monad defined in the previous installment we evaluate our tangle. We then compute the coefficients a and b, which are polynomials in A. Now compute a and b at A²=-i, and then compute -ia/b and use the operations z→z+1, z→z-1 and z→-1/z to reduce this rational to zero, interpreting these operations as twist, antitwist and rotation respectively (as I described here).

There's a shortcut I'm going to take but I don't have a proof it's guaranteed to work so it's risky. Instead of computing polynomials in A, doing polynomial division, and then evaluating at A²=-i, I'm going to do everything from the first step with the assumption that A²=-1. (The catch is that if a and b are polynomials in A, then a and b may have a common divisor that cancels when computing a/b. If we substitute A²=-1 from the beginning we get division of zero by zero. I don't know for sure if this can't happen. It's not a big deal if it does though, I just need to write a polynomial division routine.)

Now if A²=-1 then we could choose A=(1/2)^1/2+(1/2)^1/2i. That suggests using real arithmetic but we need to extract results that are rational numbers. So instead introduce the field K=Q((1/2)^1/2), ie. the field of numbers of the form u+v(1/2)^1/2. We can implement this as:


> data K = K { re2::Rational, im2::Rational } deriving (Eq,Show)

> instance Num K where
>   K a b + K a' b' = K (a+a') (b+b')
>   K a b * K a' b' = K (a*a'+2*b*b') (a*b'+a'*b)
>   negate (K a b) = K (negate a) (negate b)
>   abs _ = error ""
>   signum _ = error ""
>   fromInteger i = K (fromInteger i) 0

> instance Fractional K where
>   recip (K a b) = let r = recip (a*a-2*b) in K (r*a) (-r*b)
>   fromRational x = K x 0

We actually want to work with the field L=K(i), ie. complex numbers built from elements of K. We can implement this as:


> type L = Complex K

(Notice how it's not at all hard to build up algebraic structures like this in Haskell without a symbolic algebra package. For more powerful code in this vein, see David Amos's code).

Now we have our field, our monad is that of vector spaces over L:


> type KnotMonad = V L

Let's define some useful constants. The square root of 1/2, A, B=1/A and i:


> hr2 = K 0 (1%2)
> a = hr2 :+ hr2
> b = hr2 :+ (-hr2)
> i = 0 :+ 1

Here's the knot monad code I've presented before. (This version uses -|()|- in a way that may appear a bit useless. I'm just making explicit when we're working with a 1-dimensional vector space generated by -|()|-, which is much the same thing as the base field.)


> cup :: (Bool,Bool) -> KnotMonad ()
> cup (u,v) = case (u,v) of
>   (False,True) -> (-i * b) .* return ()
>   (True,False) -> (i * a) .* return ()
>   otherwise -> vzero

> cap :: () -> KnotMonad (Bool,Bool)
> cap () = (-i * b) .* return (False,True) .+ (i * a) .* return (True,False)

> over :: Tangle
> over (u,v) = a .* do
>   () <- cup (u,v)
>   cap ()
>   .+
>   b .* return (u,v)

> under :: Tangle
> under (u,v) = b .* do
>   () <- cup (u,v)
>   cap ()
>   .+
>   a .* return (u,v)

Here's where we work out the coefficients a and b. We treat a tangle as the function that it is and look at what happens when we apply it to two different values in (Bool,Bool). Applying to value' to infinity and zero shows that they do in fact compute a and b correctly:


> value' p =
>   let alpha = coefficient (False,False) $ p (False,False)
>       beta  = coefficient (True,False)  $ p (False,True)
>   in (alpha,beta)

From a and b we extract -ia/b as a rational:


> value :: Tangle -> Rational
> value p =
>   let (a,b) = value' p
>   in re2 $ realPart $ -i*a/b

And now we reduce to 0. This isn't the best algorithm for reducing rationals to zero. But it's what I wrote in my first attempt and I'm sticking with it for consistency with part one:


> steps :: Rational -> [String]
> steps 0 = []
> steps q | q<= -1 = "twist":steps (q+1)
>         | -1<q && q<1 = "rotate":steps (-1/q)
>         | q>=1 = "antitwist":steps (q-1)

> untangle t = steps (value t)

Here's the example from the first installment:


> example :: Tangle
> example (a,b) = do
>   (c,d) <- over (a,b)
>   (e,f) <- cap ()
>   (g,h) <- over (c,e)
>   (i,j) <- over (f,d)
>   (m,n) <- cap ()
>   (k,l) <- cap ()
>   (q,r) <- over (h,k)
>   (s,y) <- over (l,i)
>   (o,p) <- over (n,g)
>   (t,u) <- under (p,q)
>   (v,w) <- under (r,s)
>   (x,z) <- over (y,j)
>   cup (o,t)
>   cup (u,v)
>   cup (w,x)
>   return (m,z)

I'll leave you to type untangle example in ghci and refer back to my original video.

To find out more, and to read about the details I've left out, try the many papers by Kauffman, including this one.

That was hard! I didn't anticipate how much stuff I'd have to introduce to just sketch out how this works. It's one of those cases where the comments need to be about 30 or 40 times larger than the code and done properly, there's enough material to fill a good few chapters of a book. So I apologise if I skimped on some details along the way, and left out some needed explanations. But I hope that I've at least motivated you to take a peek at knot theory and appreciate how nicely tangles and monads play together.

And remember this is only code to illustrate some mathematics - it's inefficient and not how I'd really untangle big rational tangles! (For example, just using the monad associativity laws you can often rearrange the monad expressions for knots into equivalent forms that execute much quicker.)

The rest is similar to code I've used in previous installments:


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

> class Num k => VectorSpace k v | v -> k where
>     vzero :: 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

> 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
>     vzero = V []
>     V x .+ V y = V (x ++ y)
>     (.*) k = (>>= (\a -> V [(k,a)]))

> coefficient b (V bs) = maybe 0 id (lookup b (map swap (reduce bs)))

> data Complex a = (:+) { realPart :: a, imagPart :: a } deriving (Eq,Show)

> instance Num a => Num (Complex a) where
>   (a :+ b) + (a' :+ b') = (a+a') :+ (b+b')
>   (a :+ b) * (a' :+ b') = (a*a'-b*b') :+ (a*b'+a'*b)
>   negate (a :+ b) = (-a) :+ (-b)
>   fromInteger n = fromInteger n :+ 0
>   abs (a :+ b) = undefined
>   signum (a :+ b) = undefined

> instance Fractional a => Fractional (Complex a) where
>   recip (a :+ b) = let r = recip (a*a+b*b) in ((a*r) :+ (-b*r))
>   fromRational q = fromRational q :+ 0

On writing Python one-liners.

2008-09-26T08:36:00.000-07:00


# I'm off to the Oregon Shakespeare Festival for the weekend so no
# time to finish my Untangling series. But I do have time to post this...

# The other day someone asked me if it was possible to implement a
# certain task in one line of Python. My response was "yes, but you
# wouldn't want to." You can think of what follows as the demonstration
# that you really don't want to.

# The problem is, there are lots of Python constructions you can't
# do on one line: assigning variables, defining functions, using
# conditionals, recursion and so on. So it appears there's a limit
# to what you can do. Amazingly, all of these problems can be solved
# with the use of lambdas.

# One place you often want one line lambdas is in the use of the map
# function. So let's stick with that as an example and warm up with
# an easy example. Lets say we want to add 1 to all of the numbers
# from 0 to 19. We can write:

def inc(n):
    return n+1

print "1.",map(inc,range(0,20))

# Lambdas enable you to write that as a one-liner:

print "1.",map(lambda x:x+1,range(0,20))

# Now what happens if we want to apply a function that uses a local
# (constant) variable. For example:

def f1(n):
    m = 2*n+1
    return m+m*m

print "2.",map(f1,range(0,20))

# We can eliminate the assigment by means of a function. m changes from
# being a local to a function argument:

def f2(n):
    def g(m):
        return m+m*m
    return g(2*n+1)

print "2.",map(f2,range(0,20))

# But that seems to make things worse because now we have a function
# definition instead of a variable assigment. But we can eliminate
# that with a lambda:

def f3(n):
    return (lambda m:m+m*m)(2*n+1)

print "2.",map(f3,range(0,20))

# And now we can write the whole thing as one big lambda:

print "2.",map(lambda n:(lambda m:m+m*m)(2*n+1),range(0,20))

# But what happens if we want to define a bunch of local functions
# inside our lambda? We can use the same technique. We (1) use
# lambdas to allow us to fake local variable definitions and
# (2) use lambdas to define the functions. So we can replace this:

def f(n):
    def g(n):
        return 2*n+1

    def h(n):
        return 3*n-2

    return g(h(g(h(n))))

print "3.",map(f,range(0,20))

# with this:

print "3.",map(lambda n:(lambda g,h,n:g(h(g(h(n)))))(lambda n:2*n+1,lambda n:3*n-2,n),range(0,20))

# We're passing in lambdas as arguments to a lambda so that they
# get bound to g and h.

# But what happens if we want a recursive definition like:

def fact(n):
    if n<=1:
        return 1
    else:
        return n*fact(n-1)

print "4.",map(fact,range(0,20))

# Because the function is bound to the name fact, we can refer to
# fact by name inside the defintion of fact. How can we do this in
# a one-line lambda? We can start by eliminating the self-reference
# in a fairly standard way. First rewrite fact so it doesn't call
# itself, but calls a function passed in as an argument:

def urfact(f,n):
    if n<=1:
        return 1
    else:
        return n*f(f,n-1)

# We can then make a factorial function by calling urfact with itself,
# allowing it to carry out the recursion:

def fact1(n):
    return urfact(urfact,n)

print "4.",map(fact1,range(0,20))

# We can use urfact directly without fact1:

print "4.",map(lambda n:urfact(urfact,n),range(0,20))

# Now all we need is a one-liner definition of urfact. It's tempting to try:

def urfact1(f,n):
    return {True:1,False:n*f(f,n-1)}[n<=1]

# But that results in non-termination because it evaluates *both* branches
# of the conditional. In a lazy language like Haskell, the branch not taken
# wouldn't be evaluated. But Python can be made lazy by using...you
# guessed it...lambda. The expression lambda:x is like a lazy version
# of the expression x. It doesn't get evaluated until you apply it as a
# function of no arguments, at which time it takes on the value x. So we
# can rewrite urfact as:

def urfact2(f,n):
    return {True:lambda:1,False:lambda:n*f(f,n-1)}[n<=1]()

# And we get

print "4.",map(lambda n:urfact2(urfact2,n),range(0,20))

# Note how we use urfact2 twice, so we need the trick above for local
# definitions to get it down to one reference:

print "4.",map(lambda n:(lambda f:f(f,n))(urfact2),range(0,20))

# And now we can replace urfact2 with a lambda

print "4.",map(lambda n:(lambda f:f(f,n))(lambda f,n:{True:lambda:1,False:lambda:n*f(f,n-1)}[n<=1]()),range(0,20))

# And that's enough for now. I expect that if we keep going we'll find 
# way to rewrite any Python program as a one liner. You could probably
# even write a program to automatically 'compile' a usable subset of
# Python into one liners.

# But you probably don't want to do that.

Untangling with Continued Fractions: Part 4

2008-09-20T16:07:00.000-07:00


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

> module Main where

> import Test.QuickCheck
> import qualified Data.Map as M
> import Control.Monad
> import Ratio

> infixl 5 .+
> infixl 6 .*

As I've mentioned previously, we need to choose values for our pieces of knot or tangle so that when they're combined together we get something that tells us about the knot or tangle, not the diagram. For example, suppose our diagram contains this piece:

That piece is isotopic to:

We want any kind of evaluation of our diagram to be independent of which one of these we've chosen. In other words, among other things we need this function:


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

to equal this function:


> straight i = do
>   return i

I've already said how in the vector space monad, these correspond to summations, so really we want these two expressions to be equal (using the summation convention):

M^ijN_jk=δⁱ_k

(δⁱ_k is 1 if i=k and zero otherwise, ie. the identity matrix.)
The left hand side is just the rule for multiplying matrices. So we require M to be the inverse of N.

Now we can start on some real knot theory. There is a collection of 'moves' we can perform on a diagram called Reidemeister moves. Here are diagrams illustrating them:

There are three types of move that we call type I, type II and type III.

These diagrams are intended to show changes we can make to a part of a bigger knot or tangle. It's not hard to see that these are isotopies, they're just simple operations we can perform on a knot or tangle by moving the string a little bit. If you perform one of these moves on a knot or tangle then you should end up with merely a different diagram representing the same knot or tangle. Each of these corresponds to an equality of monadic expression, or an equality of summations. But the really important part is that these are all you need to consider. Every isotopy can be constructed out of these moves. So if we can satisfy the corresponding equalities of expressions then we automatically have a knot invariant that tells us about the underlying knot or tangle, not just the diagram.

Unfortunately, there's a catch. Satisfying the equality for a type I move is too onerous. So we'll just concentrate on type II and III. But the good news is that it doesn't matter, there are workarounds.

Let a be an unknown complex number. We're going to build expressions from a but they're all going to be polynomials in a and its reciprocal, like a+a^2-3/a. These are known as Laurent polynomials and I pointed out recently that I can use the Group type class to represent them.


> i = 0 :+ 1
> newtype Z = Z Int deriving (Eq,Ord,Show,Num)
> type Poly = V (Complex Rational) Z

In other words, elements of Poly are weighted linear combinations of integers. We can think of each integer n as representing aⁿ and the weights being coefficients of these powers. Amazingly the group multiplication multiplies these things out correctly.

a and ia represent a and a^-1.


> a = return (Z 1) :: Poly
> ia = return (Z (-1)) :: Poly

We can now satisfy the type II and III Reidemeister moves with these definitions:


> cup :: (Bool,Bool) -> V Poly ()
> cup (u,v) = case (u,v) of
>   (False,True) -> (-i .* ia) .* return ()
>   (True,False) -> (i .* a) .* return ()
>   otherwise -> zero

> cap :: V Poly (Bool,Bool)
> cap = (-i .* ia) .* return (False,True) .+ (i .* a) .* return (True,False)

> over :: (Bool,Bool) -> V Poly (Bool,Bool)
> over (u,v) = a .* do
>   () <- cup (u,v)
>   cap
>   .+
>   ia .* return (u,v)

> under :: (Bool,Bool) -> V Poly (Bool,Bool)
> under (u,v) = ia .* do
>   () <- cup (u,v)
>   cap
>   .+
>   a .* return (u,v)

(I'm making the V monad do double duty here. I'm using it to represent the Laurent polynomials, but I'm also using it to represent vectors over the Laurent polynomials.)

For example, here's what a type III move looks like if we try to arrange things nicely in terms of our standard 'components':

That gives rise to the two functions


> left (i,j,k) = do
>   (l,m) <- under (i,j)
>   (n,o) <- over (m,k)
>   (p,q) <- over (l,n)
>   return (p,q,o)

and


> right (i,j,k) = do
>   (l,m) <- over (j,k)
>   (n,o) <- over (i,l)
>   (p,q) <- under (o,m)
>   return (n,p,q)

whose equality can be tested with


> test1 = quickCheck $ \(x,y,z) -> left (x,y,z)==right (x,y,z)

I'll leave checking the type II rule as an exercise.

With these defintions we can take any knot diagram and compute a value for it. Note that because a not has no loose ends we have no inputs and no outputs so the result won't be a function but a value, a polynomial in a and its reciprocal. This is known as the Bracket Polynomial.

As mentioned above, the bracket polynomial fails to be a knot invariant. Two diagrams related by type I moves will have different bracket polynomials. With a little bit of sleight of hand we can apply a fudge factor to l with those extra loops and we end up with the Jones polynomial. As that's not my main goal here, I'll let you read how to use the writhe get from the bracket polynomial to the Jones polynomial at Wikipedia. But do I want to make another digression.

It's tricky to define knot invariants. The problems is that knots are topologically invariant in the sense that you can deform a knot in lots of ways and it's still the same knot. This rules out most attempts at using the geometry of a knot to define an invariant. Just about any geometrical feature you find in a particular configuration will be ruined by deforming the knot. But here's an idea: suppose you could compute some geometrical property from a particular curve representing a knot and then find the average of that value for all possible curves representing the same knot. Then we'd get a knot invariant as it wouldn't depend on a particular choice of representative curve. But there's a catch. We need to average over an infinite space, the space of all curves representing a knot. On the other hand, there is a perfectly good tool for averaging over infinite sets: integration. But there's another catch: the space of all curves representing a knot is very large, infinite dimensional in fact, and we don't have very good theory for how to do this. At least mathematicians don't. But physicists integrate over spaces of paths all the time, ever since Feynman came up with the path integral. Amazingly, Ed Witten showed that tere was already a physical model that fit the bill: Chern-Simons theory. In fact, Witten showed that the Jones polynomial, which as I've already mentioned came originally out of statistical mechanics, could be derived from this theory.

Anyway, enough of that digression. In the next installment I'll show how the above can be modified for use with tangles and give my final algorithm.

In the last installment (I think) I'll show how we can derive the rational number for a tangle using the bracket polynomial.

The following is the 'library' behind the code above. Most of it is simply built up from what I've said in recent weeks.


> 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

> 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)))

> diag a = (a,a)
> both f g (a,b) = (f a,g b)

> class Group m a where
>     unit :: () -> m a
>     counit :: a -> m ()
>     mult :: (a,a) -> m a
>     comult :: a -> m (a,a)
>     anti :: a -> m a

> instance Monad m => Group m Z where
>     unit _ = return (Z 0)
>     counit _ = return ()
>     mult = return . uncurry (+)
>     comult = return . diag
>     anti = return . negate

> newtype Identity a = I a deriving (Eq,Ord,Show)

> instance Monad Identity where
>     return x = I x
>     I x >>= f = f x

> data K = K Rational Rational deriving (Eq,Show)

> instance Num K where
>   K a b + K a' b' = K (a+a') (b+b')
>   K a b * K a' b' = K (a*a'+2*b*b') (a*b'+a'*b)
>   negate (K a b) = K (negate a) (negate b)
>   abs _ = error ""
>   signum _ = error ""
>   fromInteger i = K (fromInteger i) 0

> instance Num k => Num (V k Z) where
>   a + b = a .+ b
>   a * b = do
>       u <- a
>       v <- b
>       mult (u,v)
>   negate a = (-1) .* a
>   fromInteger n = fromInteger n .* return 0
>   abs a = error ""
>   signum a = error ""

> data Complex a = (:+) { realPart :: a, imagPart :: a } deriving (Eq,Show)

> instance Num a => Num (Complex a) where
>   (a :+ b) + (a' :+ b') = (a+a') :+ (b+b')
>   (a :+ b) * (a' :+ b') = (a*a'-b*b') :+ (a*b'+a'*b)
>   negate (a :+ b) = (-a) :+ (-b)
>   fromInteger n = fromInteger n :+ 0
>   abs (a :+ b) = undefined
>   signum (a :+ b) = undefined

> instance Fractional a => Fractional (Complex a) where
>   recip (a :+ b) = let r = recip (a*a+b*b) in ((a*r) :+ (-b*r))
>   fromRational q = fromRational q :+ 0

> instance Ord (Complex Float) where
>   compare (a :+ b) (c :+ d) = compare (a,b) (c,d)

Two Papers and a Presentation

2008-09-16T09:39:00.000-07:00

I've been slow to write the last couple of parts on my untangling series as I've been pretty busy lately. Among other things I've finally got around to submitting a couple of graphics related papers (that I actually wrote ages ago) for publication. Not sure they'll get accepted but there's no harm trying. I'll probably write blog posts about these in the near future now that I have permission from my employer to talk about them.

Note: Neither of these are anything I actually use at work. Well, I've used the adjoint one a tiny bit. The latter is in the spirit of computational geometry papers: it demonstrates an algorithm that runs (asymptotically) in the fastest possible time, but in practice you'll probably want to use an approximation that looks just as good. I was amazed I was able to find a practical application for a result that came out of toric varieties, but it's going to be tough on the reviewers as there seem to be so few people in the graphics world who know about multidimensional DSP techniques.

Two Tricks for the
Price of One: Linear Filters and their Adjoints

Fast and Exact
Convolution with Polygonal Filters

Of course, straight after submitting I noticed a bunch of pieces missing from some expressions. In places designed to help disambiguate things for the reviewers as well. Ho hum.

I also recently presented a very short presentation at the "Stupid Rat Tricks" part of the Renderman User Group: quine.sl: a self-shading shader. (For more on RSL see Wikipedia.) The source is available here. I don't think the company lawyers fully appreciated that this shader now allows their handiwork to be used like a texture map.

My presentation was completely blown away by the guy who implemented a large part of the Renderman interface completely in bourne shell. I can't think of a stupider thing than that! :-)

Back to untangling in the next post.

Update: I added a Mathematica 5 notebook illustrating the transfer functions in the polygonal filter paper.

Untangling with Continued Fractions: Part 3

2008-08-31T15:40:00.001-07:00

It's time to return to the vector space monad. But first, after some Haskell preamble, I need to talk about the Einstein summation convention.


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

> module Main where

> import qualified Data.Map as M
> import Control.Monad
> import Test.QuickCheck

> infixl 4 .+
> infixl 4 .-
> infixl 5 .*

We're going to be dealing with multidimensional arrays of values that I'll call tensors. We can refer to the individual elements of these arrays by providing tuples of indices drawn from some set of indices. For example we might refer to the elements of a 10 by 10 array by using pairs of indices, each drawn from the set {0,1,2,...,9}. We'll write some indices as subscripts and some as superscripts. Here are some examples: vⁱ, w_j,T^ij,W^ij_kl. vⁱ means the ith value in array v, T^ij means the (i,j)-th element of T and so on.

We can form sums over such indices like this

Σ_jkT^ijW_jkv^k

where the sum is over the corresponding index sets for j and k.
The Einstein summation convention is simply this: if an expression has an index repeated once as a subscript and once as a superscript then we implicitly assume a sum over the index without writing the summation symbol. So the above expression becomes simply

T^ijW_jkv^k

The vector space monad allows us to express something similar using do-notation. Let's use the type Bool for our index set, so we're selecting elements of our arrays using the values {False,True} rather than integers. Let e_i be the (2-dimensional) vector with the property that its j-component is 1 for j=i and 0 for j≠i. In do-notation this is represented by return i. We use superscripts to indicate the components of a vector so the components of e_i are e^j_i. We can write the vector v=3e_False-5e_Truea, with components (3,-5) as


> v :: V Float Bool
> v = 3 .*return False .- 5 .*return True

Dual to vectors we have covectors. You can think of these as linear functions that consume vectors and return reals. We can define a basis of covectors fⁱ with the property that fⁱ(e_j) is 1 for i=j and 0 otherwise. In our Haskell code it'll be more convenient to think of covectors as mapping vectors to the one-dimensional vector space whose basis element is labelled by (). When we make a vector using an expression like return True, that True is a constructor for the type Bool. We expect to use the dual operation for covectors, and that dual is case analysis. So, for example, we can write w = f^False+2*f^True, wih components (1,2), as


> w :: Bool -> V Float ()
> w i = case i of
>   False -> return ()
>   True -> 2 .* return ()

Applying w to v gives the expression vⁱw_i using the summation convention. Using do-notation it looks like


> vw :: V Float ()
> vw = reduce $ do
>   i <- v
>   w i

The point to note is that just like with the summation convention, we don't need to explicitly sum over i. i serves as a kind of iterator object and a vector v automatically supplies i with values in the appropriate range. w then consumes i. So we can translate the summation convention to do-notation by thinking of upper indices corresponding to places that iterators can emerge from, and thinking of lower indices as places where iterators can be consumed.

Another example is matrix representing a transformation of a vector to a vector. It must consume one iterator and emit another one. In the Einstein summation convention such a transformation might be written rⁱ=Mⁱ_jv^j. We could write a simple 2D 2x scale and 90 degree rotation matrix as


> m :: Bool -> V Float Bool
> m j = case j of
>   False -> 2 .*return True
>   True -> -2 .*return False

We can then write the transformation as


> r = do
>   j <- v
>   i <- m j
>   return i

Another simple example is an inner product. We can represent the ordinary dot product as a matrix d_ij and write the square length of v as v . v = d_ijvⁱv^j. We can implement this as


> v2 = do
>   i <- v
>   j <- v
>   d (i,j)

where we define


> d (i,j) = case (i,j) of
>   (False,False) -> return ()
>   (True,True) -> return ()
>   otherwise -> mzero

Similarly, tuples of superscripts correspond to iterators that are tuples of indices.

With that in place, we're now in a position to consider what monadic notation for knots means in the vector-space-over Bool monad. Essentially what we're doing is representing 'components' of a tangle or knot, the cups, caps and crossings, as tensors. When we glue these components with labels we're forming sums. In fact, each labelled edge connecting two components corresponds to an index and so we can think of it as possibly taking two different values, False or True. So the monadic expression corresponds to a sum over all possible ways to assign values from Bool to the indices with a corresponding weight for each term. This should be familiar to some people, it's just like what takes place in statistical thermodynamics. We consider a bunch of particles and form a sum over all of the possible states that those particles could have, with a corresponding weight for each configuration of states. A good example is the Ising model in which electrons sit at the vertices of a regular lattice. Two neighbouring electrons in a ferromagnet tend to want to align the same way, so confgurations in which neighbouring electrons have the same spin will get assigned a higher weight.

In our tangle problem, the trick is to find a magic set of weights with the property that two different summations representing isotopic knots or tangles give the same result. This turns out to be too hard, but we relax it a bit we will be led to one of the most fruitful new approaches in knot theory in recent years - the Jones Polynomial. Vaughan Jones was investigating certain types of model from statistical mechanics and came to realise that some of the equations that were arising were the same ones needed to get suitable weights for knot theory. From this realisation emerged his eponymous polynomial. But, that's all I'll be saying about this polynomial in this post.

Anyway, to make things clearer here's an example:

We interpret cap, cup, over and under as the tensors M^ij, N_ij, S^ij_kl and T^ij_kl. In monadic notation we get


> example (i,j) = do
>   (k,l) <- over (i,j)
>   (m,n) <- cap
>   cup (l,m)
>   return (k,n)

With the summation convention in mind this evaluates to:

S^ij_klN^lmM_m,n

and without it we see that all we're doing is computing

Σ_l,mS^ij_klN^lmM_m,n

We can try an elementary choice for our tensors above. We can define M and N to be given by identity matrices:


> cap :: V Float (Bool,Bool)
> cap = return (False,False) .+ return (True,True)

> cup :: (Bool,Bool) -> V Float ()
> cup (i,j) = case (i,j) of
>   (False,False) -> return ()
>   (True,True) -> return ()
>   otherwise -> mzero

We can then choose S and T to simply swap their arguments:


> over,under :: (Bool,Bool) -> V Float (Bool,Bool)
> over (i,j) = return (j,i)
> under (i,j) = return (j,i)

With these choices, the weighting is such that a non-zero weight is assigned to a knot or tangle if and only if all of the labels along a single strand are assigned the same value, in other words we're now summing over ways of assigning labels to entire strands. For rational tangles we can only end up with three possible values corresponding to the three ways the two outputs can be connected to the two inputs. (Straight through, the tangle I called ∞; straight through with a swap, the tangle called -1; and connecting together both the inputs, and both the outputs, the tangle called 0.)


> tangleInfinity,tangleMinus1,tangleZero :: (Bool,Bool) -> V Float (Bool,Bool)
> tangleInfinity (i,j) = return (i,j)
> tangleMinus1 (i,j) = return (j,i)
> tangleZero (i,j) = do
>   cup (i,j)
>   cap

> test1 = quickCheck $ \x -> (example x == tangleInfinity x)
> test2 = quickCheck $ \x -> (example x == tangleMinus1 x)
> test3 = quickCheck $ \x -> (example x == tangleZero x)

Only test2 succeeds so we know example is a tangle that forms its outputs by swapping the inputs.

So even with this simple assignment we're able to extract some basic connectivity information from our monadic representation. With smarter choices we'll get more information, but that's enough for now. We now have the prerequisites and in the next installment we can start doing some proper knot theory.


> data V r b = V [(r,b)] deriving (Eq,Show)
> unV (V bs) = bs

> instance Num r => Functor (V r) where
>    fmap f (V bs) = V $ map (\(r,b) -> (r,f b)) bs

> scale x (bs) = map (\(r,b) -> (x*r,b)) bs

> instance Num r => Monad (V r) where
>    return b = V [(1,b)]
>    x >>= f = join (fmap f x)
>        where join x = g $ fmap unV x
>              g (V bs) = V $ concat $ map (uncurry scale) bs

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

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

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

> instance (Num r,Ord b) => VectorSpace r (V r b) where
>    zero = V []
>    V v1 .+ V v2 = reduce (V (v1 ++ v2))
>    r .* V v = reduce (V $ scale r v)

Untangling with Continued Fractions: Part 2

2008-08-23T15:00:00.000-07:00

Recall that I defined a rational tangle to be what you get by starting with a pair of untangled strings whose ends are at infinity and then sliding the strings around as much as you like keeping the ends at infinity. If you don't like the infinity bit, just think of strings inside a sphere with the ends constrained to lie on (but free to slide around) the surface of the sphere. We also insist that the ends end up at the four diagonal directions.

Two rational tangles are considered to be isotopic if you can slide the strings about to turn one into the other, this time keeping the ends of the strings fixed. (Such a motion is called an isotopy.)

As I've previously mentioned, you only need twists, antitwists and rotations to untangle a rational tangle. Here are the twist operations again:

I'm using a box to represent some unknown tangle. The antitwist really is an inverse for the twist because one followed by the other gives you what you started with once you pull the strings straight:

We can also draw a rotation like this:

Let's call the twist T (so the antitwist is T^-1) and the rotation S and write products of operations from right to left in the usual way. As both S and T are invertible, they generate a group. It's pretty clear that S⁴=1, where 1 is the identity operation. So we don't need a clockwise rotation, we can just use S³.

But here's a surprising fact: S²=1. If we perform a 180 degree rotation on a rational tangle the result might look different at first, but it is isotopic to what we started with. I'll leave you to convince yourself of this, but you can prove it by induction on the number of twists or antitwists you have in your diagram. Try it for simple rational tangles first, maybe real ones made with real string.

Here's an even more surprising fact, (TS)³=1. I can demonstrate this pictorially. Start with he result of applying T:

The result of ST:

The result of TST:

The result of STST:

After TSTST:

But by sliding one string behind the tangle, allowed in an isotopy, we find that that is isotopic to this:

It's not hard to see that another S brings us back to where we started. (To do that we need to rotate the block labelled A, but that is allowed in an isotopy as we're not moving the ends of the strings.) This means we can write an antitwist in terms of twists and rotations. But it's convenient to keep using antitwists.

So we have S²=1 and (ST)³=1. It can also be proved that any other relationships between S and T can be derived from these using the operations in a group. But this is a well studied group. We can approach it like this. Define two functions s and t by

s(x) = -1/x

t(x) = x+1.

It's easy to show that s²=identity=(st)³. In fact, we have an isomorphism between the group of operations we can perform on rational tangles and the functions we can build with s and t. It goes further.

Take the rational numbers and throw in ∞ to give us the extended rationals. Let s and t act in the usual way on most of the rationals but throw in the extra rules s(∞)=0, s(0)=∞ and t(∞)=∞ to handle infinity. We have defined an action of s and t on the extended rationals. We can exactly parallel this with the rational tangles. Call this tangle 0:

We can reach every rational tangle by applying the operations S and T to 0. But if we replace S with s and T with t then the sequence of s's and t's will give us an extended rational number that labels the tangle. Because S and T obey the same equations as s and t, isotopic knots will end up getting the same extended rational. So we can use extended rationals to label the rational tangles. Here are some examples:

Firstly just an antiwist, T^-1. That corresponds to t^-1(0)=-1.

Two antitwists gives -2:

Composing a 90 degree rotation gives

And applying another antitwist subtracts one again giving:

The name rational tangle is justified.

So, given any any rational tangle we simply need to find the corresponding extended rational and then figure out how to write it as a composition of s's and t's applied to 0. This is more or less the usual technique for finding the continued fraction representation of a rational. And once we know the sequence, we can then apply it in reverse to untangle the tangle.

But there's one catch, how can we go from the monadic representation of a tangle that I described previously to the corresponding extended rational? We'll start on that problem in my next post, but amazingly the vector space monad can be used to do almost all of the work. It's surprising: it's only a short bit of code that works out the untangling operations, but the explication is getting rather long now...

An amazing application of this work is in the study of the untangling of DNA. DNA is a double helix and the two parts of the helix must be separated in order to transcribe it. As you might imagine, the DNA can become incredibly tangled by this process. To deal with this, the topoisomerases snip, manipulate and repair the DNA so as to eliminate these tangles. It was realised that it was possible to work out exactly what operations the topoisomerases were performing by analysing the before and after DNA strings in terms of rational tangles.

The perosn who first figured out this relationship between tangles and rational numbers was the ubiquitous JH Conway.

References

For the proofs missing above, these references should fill in the details:
Rational Tangles, Kauffman and Goldman
Knots and Physics, Kauffman
And if I've lost you above, the following is a nice easy introduction:
Conway's Rational Tangles, Davis
Modeling protein-DNA complexes with tangles

One last thing. I've noticed some, but not all, papers define the twist oppositely from me.

Untangling with Continued Fractions: Part 1

2008-08-16T20:38:00.000-07:00

Continuing from before It's time now do discuss the problem of forming a representation of a knot or tangle in machine readable form.

I've defined rational tangles but I should set that in a wider context. A knot is essentially a single closed loop embedded in 3D space and a link is a non-intersecting union of a bunch of these. Any tangle, knot or link can be projected down to 2D space so that you get a finite collection of 'over' or 'under' crossings connected by arcs. You might have to jiggle things around a bit to ensure that you don't get any degeneracies like two separate parts of a knot being projected to the same segment of curve in 2D, but this is always possible. Here's an example of such a 2D projection, a knot diagram, for the photograph I posted last week:

Roughly speaking, two knots or links or tangles are equivalent if you can slide the strings about to get from one to the other without passing one string through another. In the case of rational tangles we have the extra constraint that the free ends mustn't move and must always remain 'at infinity' so you can't pass loops over the end. Given two diagrams for rational tangles, the task is to tell when they represent something equivalent. I'll start on this in a future posting, but for now I just want to consider the problem of getting diagrams like that above into machine readable from. (Incidentally, the rigorous mathematical definition of knot equivalence, ambient isotopy, relies on the notion of finding a continuous bijection between the space around the two knots, not the knot itself.)

We don't need anything particularly clever to form our representation as there is a fairly commonplace way to represent connectons between components using monads. If we have some kind of block with inputs a and b and outputs c and d we can represent this as a line of do-notation looking like



(c,d) <- block (a,b)

Examples of such notation are in Matthew Naylor's article on Lava. We could also use arrows but in this case there is no need.

So now we need to break up a diagram of a rational tangle into components and hook them up, working down the diagram from top to bottom. There are four essential components that can be discerned and I'm calling them cups, caps, overs and unders. Here are the corresponding diagrams for each one:

Cups have two inputs but no output and caps have two outputs and no inputs.

Now all we have to do is label the inputs and outputs of each arc, write the corresponding lines as above, and collect them together in a do-block. Here I've redrawn the above diagram marking the cups and caps with red circles and the overs and unders with green circles.

I've also labelled the inputs and outputs. We can now write the following block of code.



> example (a,b) = do
>   (c,d) <- over (a,b)
>   (e,f) <- cap
>   (g,h) <- over (c,e)
>   (i,j) <- over (f,d)
>   (m,n) <- cap
>   (k,l) <- cap
>   (q,r) <- over (h,k)
>   (s,y) <- over (l,i)
>   (o,p) <- over (n,g)
>   (t,u) <- under (p,q)
>   (v,w) <- under (r,s)
>   (x,z) <- over (y,j)
>   cup (o,t)
>   cup (u,v)
>   cup (w,x)
>   return (m,z)

Note how this block has a pair of inputs (a,b) and a pair of outputs (m,z) corresponding to the strings at the top and bottom. Clearly knots and links should have no inputs or outputs but rational tangles should have a pair of inputs and a pair of outputs. The precise order shouldn't matter as long as you don't try to use a name that isn't yet in scope.

So what monad should we use? It's easy to imagine some kind of state monad that allows us to generate fresh labels for each of the connections and collects up a graph-like representation of our tangle. But the surprise is this: it turns out we don't need to do anything complicated like this. With suitable definitions of cup, cap, over and under, not only does the vector space monad give us the representation we want, it also does most of our computation for us. But first I need to explain the underlying mathematics in an upcoming installment.