On representing some real numbers exactly

There are uncountably many real numbers, and only countably many finite strings of symbols, so we know for sure that there is no scheme to represent all real numbers as finite strings of symbols in such a way that different reals get different representations. However, there are useful subsets of the real numbers that *are* equipped with finite encoding schemes. I hope to give some examples of such sets, including one, the set of periods, that isn't currently well known outside the mathematical world.

It's traditional to think of the real numbers as divided up into three types of number: the rational numbers, the algebraic numbers, and the leftovers known collectively as the transcendental numbers. If you're not familiar with those terms I'll be explaining them below. We can represent this as follows:

ℚ⊂A⊂ℝ

I'm using unicode for blackboard bold R and Q there. I hope your browser supports them. ℝ is the set of real numbers, ℚ is the set of rational numbers, and A is the algebraic numbers.

But there are many different kinds of transcendental number in existence, so the subject of this article will be three more levels of classification we could insert between the set of algebraic numbers and the set of real numbers in the sequence above.

But first, let's revisit ℚ and A.

The Rational Numbers, ℚ
The rational numbers are fairly straightforward to deal with. They are nothing more than the ratios of integers. Given two such ratios we can easily compare them to see whether they are equal. For example, we're taught at a young age how to tell if 2/3 is the same as 4/6.

However, it didn't take long for the Ancient Greeks to realise that they had some equations they could solve approximately using rationals, but couldn't solve exactly. The best known example is solving x²=2.

Assume x is rational, so that x=p/q, with p and q integers with highest factors cancelled out. We have p²=2q², and so p must be divisible by 2. But then p² is divisible by 4, and so q is divisible by 2. This means that we can cancel 2 from the top and bottom of p/q. That contradicts our assumption about cancellation. So the equation has no rational solution.

If we expect x²=2 to have a solution we must work with numbers that can't be represented as fractions.

Notice how the rational numbers come equipped with a notation to describe them. We can just write one integer over another, with a horizontal line between them. Given such a description we can immediately tell if it's valid (it's only invalid if the denominator is 0) and whether it equals some other rational number.

One interesting property of the rational numbers is that they are *dense* in the real numbers. This is just another way of saying that any real number can be approximated as well as we like by a rational number. For example, we can approximate √2 to within 1/1000 by the rational number 141421/10000. This is an important property that will come up again later.

We also know that there are only countably many rational numbers because there are only countably many pairs of integers.

The Real Algebraic Numbers, A
The real algebraic numbers are all the real numbers that we can obtain by finding roots of polynomial equations whose coefficients are rational numbers. For example, the real solutions to x⁵-2x²+½=0 are all algebraic numbers. Among other things, this includes the solution to x²=2 that wasn't contained in the rationals. The algebraic real numbers also include the solutions to "insoluble" algebraic equations. (By 'algebraic' equations, I just mean polynomial equations.) Even though Galois proved we can't write expressions for the solutions to all algebraic equations using the four arithmetic operations and nth roots, these solutions still exist in the real numbers, and they are all, by definition, algebraic numbers.

We can describe an algebraic number by writing down the equation it solves, and additionally providing some description to say which root of the equation we're interested in. Unfortunately, given any algebraic number there is always an infinite number of algebraic equations that it satisfies. So like with the rationals there is some redundancy. The good news, however, is that given two such descriptions of algebraic numbers there is an algorithm to tell whether or not they describe equal numbers.

Even though there are vastly more algebraic numbers than rational numbers, they still only form a countable set. As mentioned in the previous paragraph, every algebraic number can be described by a finite string of symbols, and there are only countably many such strings.

The Computable Real Numbers, C
Now I want to insert a set between A and ℝ. These are the computable real numbers. But first think back to what I said about the rational numbers.

Pick a real number x. Given any ε we can always find a rational number that comes within ε of x. We can define a function f on the rationals with the property that f(ε) is rational and within ε of x for any rational ε. So f gives rational approximations to x to any desired accuracy. Such a function uniquely specifies x. In fact x is the limit of f(ε) as ε tends to zero. The computable real numbers are the real numbers specified in this way by *computable* functions. Notice how there are no infinities involved here. A computable real number is represented by nothing more than a finite string of symbols forming a computer program that processes finite integers. We can find out what the real number is to any desired accuracy in a finite time.

Because we can write computer programs to find roots of algebraic equations to any desired accuracy we know that the algebraic numbers are contained in the computables. So now we have:

ℚ⊂A ⊂C ⊂ℝ

Where I'm using C to represent the computables.

Almost any number we ever want to work with is computable. Whether it's √2, or π, or the value of some humongous integral that's appeared in your engineering problem, chances are there's an algorithm somewhere that can approximate it as accurately as you like, given enough CPU time and RAM.

Note that every computable is described by a finite string of symbols, and so the computables form a countable set.

There are some problems with computable numbers. Given two representations of computable numbers we have no way of telling whether or not they are equal. It's not just that this is hard to do. There simply is no algorithm to prove that the numbers are equal. The problem is that the two computer programs may give exactly the same approximations to every degree of accuracy down to a certain ε. But for higher accuracy they may give different results. We have no way of knowing in advance at what value of ε they'll start differing. And anyway, there simply is no algorithm for telling if two computer programs will always generate the same results. So paradoxically, testing the equality of two computable numbers is itself uncomputable.

Actually, the situation is much worse. We can't even tell if we have a valid representation of a computable real. In order to do that, we need to know that our computer program to generate approximations terminates. But solving that would solve the halting problem.

By the way, despite these issues this representation of computable real numbers as functions leads to practical ways to do arbitrary precision arithmetic on a computer in a way that gives us guaranteed bounds on the accuracy of our results.

From here we could go two ways. We could try to rein in our class of number to something more manageable. Or we could try to push further and try to find ways to represent real numbers that aren't even computable. Let's do the former first.

The Periods, P
The periods are an interesting class of real number that doesn't seem to be well known. The definition may be old, but it was only recently that they were being promoted in mathematical circles as an interesting thing to study. So I thought I'd contribute to their promotion.

Consider the number π. It arises straightforwardly in geometric problems. An example is computing the area of a circle. Back in the 18th century Lambert proved that it was irrational and in the 19th century Lindemann proved that it was transcendental. But in a sense the transcendental numbers are simply a rubbish heap into which the leftover numbers have been discarded. Can we get some kind of handle on at least some of these numbers in a way that puts π back on the map?

The real number π is the area of a unit circle. So it is the area of the region of the plane given by the equation x²+y²<1.



Here's a similar construction of another real number:



The area of that region is log(2). (Natural logarithm of course!) We can construct the logarithm of any positive rational number in a similar way. These are all transcendental numbers. But in a way they are nice transcendental numbers. They arise from considering areas of the plane described by straightforward algebraic inequalities. (As in the previous section, I'm using 'algebraic inequalities' to mean 'polynomial inequalities'.)

That suggests a class of real number: those numbers that can be represented as the volume of a region defined by a bunch of algebraic inequalities with rational coefficients. These are known as the periods. I'm talking about generalised volumes in n-dimensions, not just areas in the plane. Clearly π and log(2) are periods.

But does this representation solve the problem we had with computable numbers where we were unable to guarantee we could check the equality of two numbers? Now it gets interesting. The answer is: we don't know. It is conjectured that if we have the same number represented in two different ways as a period that we can transform one representation into another simply by using a small set of elementary operations. It is also conjectured that there is a terminating algorithm for finding such a sequence of operations, or if one doesn't exist, demonstrating this.

One obvious question now is this: are there any computable reals that aren't periods? A recent paper, Periods and elementary real numbers claims to exhibit one by means of a kind of diagonalisation argument. But in general it's hard to prove that a number isn't a period. I don't believe there is a proof yet that e = 2.718... isn't a period, but mathematicians expect that it isn't.

I think the name 'period' comes from the fact that the period of a pendulum of rational length, in a gravitational field of rational strength, is a period. Computing this requires elliptic functions, and these often result in periods when applied to rational numbers.

Although mathematics is about numbers, few mathematical publications talk much about specific real numbers. Curiously, when they do, they often talk about numbers that are periods. For example there have been many recent papers on the values of the Riemann zeta function for interesting arguments. These are often periods.

It has been suggested that the study of periods is actually the study of algebraic geometry in disguise. There is certainly no end to the interesting mathematics we can do using only periods.

We now have:

Q ⊂A ⊂P ⊂C ⊂R

The Definables
In the section on computable numbers I mentioned that we could go the other way and try to find a bigger class of real numbers that could be represented by finite strings of symbols. We could simply go "all the way" and consider those real numbers that can be defined, by any means possible, using the symbols of mathematics. We can try to pin this down a bit better. We'll work with the language of set theory. In this language we can write strings of symbols like S = "x>0 and x²=2". Such a string uniquely defines a real number if when we glue the string "there exists a unique x such that" onto the beginning of it we get a true proposition. We can now represent the number √2 as the string S.

It looks like we now have:

Q ⊂A ⊂P ⊂C ⊂D ⊂R

with D the set of definables.

Except there's a problem. The set D is not definable! The problem is this: we can talk about strings of symbols representing real numbers. To talk about these in set theory we'd encode these strings of symbols as sets so that we can apply the language of set theory. In order to use our attempted definition of definable we need some way to say when a string of symbols represents a true proposition. But to define a set of definables we need to talk about true propositions in the language of set theory. Godel showed us how to talk about the provability of a proposition within set theory. He did this by showing that provability is about mechanical operations we can perform on strings. But there's nothing analogous for talking about the truth of propositions. In fact, Tarski showed us this is impossible. So while we can talk about all kinds of individual numbers as being definable, we can't construct the set of definable numbers.

You may be interested to see an example of a definable number that isn't computable. Probably the most publicised example is Chaitin's constant. It represents the probability that a randomly generated string of symbols is a computer program that terminates in a finite time. We can't actually compute this number because it requires us to solve the halting problem. Nonetheless, it's perfectly well defined.

You can find the original paper on periods by Zagier and Kontsevich here. I first found out about periods from the book Mathematics Unlimited.

Update: Jared asked if I intended the real algebraic numbers. That's what I said in the text, but I confusingly wrote ℚ for this set. But that usually means the *complex* algebraic numbers. So I've changed notation and now use A for the real algebraic numbers.

Similarly, the periods are usually defined to be complex numbers whose real and imaginary parts are given by algebraically specified volumes. I was trying to avoid mention of the complex numbers to keep prerequisites to a minimum and the essential ideas here work without mention of the complex numbers.

A Partial Ordering of some Category Theory applied to Haskell

I've had a few requests from people wanting to teach themselves applications of Category Theory to Haskell based on posts I've made. I've made things difficult by posting stuff at random levels of difficulty and without any kind of organising thread through them. So here's an attempt to list a bunch of posts related to aspects of category theory. I've grouped them by themes and within each theme I've tried to list the articles in order of difficulty. Unfortunately there can be big gaps between one article and the next as none of this material was intended to be linked together continuously. Nonetheless, I hope this is of some help.

First a warning: Category Theory Screws You Up!.

Monads
The first theme has to be monads of course. But don't forget: monads don't do anything. They're simply an interface to something that you must already have implemented some other way. So don't believe all that hype about how monads are what allow Haskell to use side effects and I/O.

1. If you just want to use IO and don't care about monads: The IO Monad for People who Simply Don't Care


2. This gets more hits than anything else on my blog: You Could Have Invented Monads! (And Maybe You Already Have.)


3. The Trivial Monad


4. Homeland Security Threat Level Monad


5. Monads, Kleisli Arrows, Comonads and other Rambling Thoughts


6. The idea that monads give a way to describe substitutions in a tree forms the basis for this and the next two posts: Variable substitution gives a...


7. Monads are Trees with Grafting


8. Where do monads come from?


9. From Monoids to Monads


10. Some thoughts on reasoning and monads


11. The Mother of all Monads


12. Beyond Monads
    

Fold and Unfold
I think this is one of the great applications of Category Theory to Computer Science. Structural recursion can be characterised really nicely in terms of F-algebras. That's cool. But even cooler is that when you dualise the definitions you get a great way to look at non-terminating computations on things like streams.

1. Data and Codata


2. Coalgebras and Automata


3. Purely functional recursive types in Haskell and Python

Commutative Monads and Vector Spaces
Trying to order these is tricky. I'm not sure I define the term commutative monad until the talk.

1. I realised a while back that the operation to build a vector space from a basis is a monad. In fact, like many well known algebraic structures, we get a commutative monad. Monads for vector spaces, probability and quantum mechanics pt. I


2. Monads, Vector Spaces and Quantum Mechanics pt. II


3. Trace Diagrams with Monads


4. A talk this time: Commutative Monads, Diagrams and Knots


5. And what happens when you try to use a non-commutative monad when a commutative monad is expected: Why isn't ListT a monad?

Comonads
I'm not sure the killer application for comonads has been found yet. But I do think they're good for things like dataflow and cellular automata fit the comonad model very well:

1. Evaluating cellular automata is comonadic


2. Comonadic Arrays


3. Comonads and reading from the future


4. The Monads Hidden Behind Every Zipper


5. Categories of polynomials and comonadic plumbing

Category Theory
And these are generally categorical articles

1. You Could Have Defined Natural Transformations


2. Reverse Engineering Machines with the Yoneda Lemma


3. A Yonedic Addendum


4. Applying the Yoneda lemma to memoize functions that at first seem unmemoizable: Memoizing Polymorphic Functions with High School Algebra and Quantifiers


5. Products, Limits and Parametric Polymorphism


6. Dinatural Transformations and Coends


7. The goal here was not to understand what Category Theory applies to Haskell, but how Haskell code can be interpreted in other categories: "What Category do Haskell Types and Functions Live In?"


8. A first step towards 2-category theory: The Interchange Law

---

Products, Limits and Parametric Polymorphism

> {-# LANGUAGE RankNTypes #-}

When I wrote about memoizing polymorphic types I mentioned that you can think of forall a. F(a) as the product over all types a of F(a), where F is some type level function. For example F might be a type constructor like []. That's not completely accurate, as I hope to now explain. Along the way we should get some insight into the meaning of the limit of a functor in a category.

Suppose we have two types, A and B. We can form their product (A, B). We have the two projections fst and snd and if we have an element x in (A, B) we know that there is no necessary relationship between fst x and snd x. We can freely choose x so that each of fst x and snd x can take on any values we like in A and B.

But now consider an element of forall a. F(a). For each concrete type X we have a projection π_X::forall a. F(a) -> F(X). So it looks like a product of all types. However, we can't freely choose elements of forall a. F(a) so as to get any element of F(X) we like for each choice of X. To demonstrate this, consider an element x of the type forall a. [a]. For any choice of X we get a projection. For example, picking X to be Int or String gives:

> p1 :: (forall a. [a]) -> [Int]
> p2 :: (forall a. [a]) -> [String]
> p1 x = x
> p2 x = x

We can draw a simple category diagram representing this:




However, forall a. [a] comes with a free theorem. For this particular type we have a free theorem that says that for any f :: X -> Y, fmap f (p1 x) == p2 x. For example consider the well known function show :: Int -> String. The free theorem tells us that this diagram commutes:




So if p1 x == [3] then p2 x == ["3"]. We have lost free choice. But we have lost a lot more freedom than this. We have a commuting triangle like this for absolutely any function f :: X -> Y. It should be clear that there is no way we can pick elements of our list to satisfy all of these constraints. So x must be the empty list.

Cones
This scenario of having one projection for each type has a name. It's an example of a cone. Let's borrow the definition from Wikipedia:

Let F:J→C be a functor. Let N be an object of C. A cone from N to F is a family of morphisms with one morphism for each X,

π_X:N→F(X)

so that for every morphism f:X→Y, the following diagram commutes:



For any Haskell functor F, the free theorem tells us that we have exactly these diagrams with forall a. F(a) playing the role of N. So forall a. F(a), with its projections to F(X), forms a cone.

Limits
If F is an instance of the Haskell Functor type class, ie. an endofunctor on Hask, then the type forall a. F(a) gives us a cone. But not just any old cone. I don't know how to prove this but I'm pretty sure it's true that the free theorems for a functor are the only non-trivial relations between F(X) and F(Y) that we're forced to obey. If that's true, then, in a sense, forall a. F(a) is the "biggest" type satisfying the free theorems. We can make this more precise by saying that for any cone N with associated projections π_X, we can can map it uniquely to forall a. F(a) so that the following diagram commutes:



This special kind of cone has a name. It's called a limit. In other words, ∀a. F(a) = lim F. In fact, this is exactly how Limit is defined in category-extras.

In a sense you can think of a limit of a functor in any category as being like a product for which a version of the free theorems for a functor holds.

Colimits
A dual story can be told for existential types and colimits. But to do this we need free theorems for existential types. I'll leave that until I've figured out a nice way to derive these free theorems...

Final words
I should have written this article before I wrote one on coends. Think of it as a prequel.

The fact that we don't have complete freedom of choice when defining polymorphic elements in Haskell is what we mean by 'parametric' polymorphism. Instead of specifying one individual value for each type we define the elements in a uniform way. In a language like C++ we can use template specialisation to freely construct a rule for getting a value from a type using 'ad hoc' polymorphism. That freedom comes at a price - it becomes harder to reason about polymorphism.

It's amazing that many definitions in category theory emerge naturally (pun fully intended) from the free theorems. I keep hoping that one day I'll find a paper on exactly what is going on here that I understand. The original free theorems paper is very uncategorical in its language.

---

The Categorification of the Naturals

A heavyweight looking title, but this post is really about nothing more than doing arithmetic.

Peano Arithmetic
I've seen many articles on type level arithmetic. They all seem to share the idea that the Haskell type system can be made to perform computations by treating types as symbols that can be manipulated according to rules. But every article I have seen seems to miss the important idea that the naturals don't have to simply be empty symbols - that they are perfectly good types with elements and that the basic operations of arithmetic have nice a interpretation as functions between types. Implementing these missing pieces will also give an example of categorification.

As usual, some Haskell administration first because this post is runnable Haskell code:

> {-# LANGUAGE ScopedTypeVariables, UndecidableInstances #-}
> {-# OPTIONS -fglasgow-exts #-}

Here are what are commonly called (some of) the Peano axioms defining addition and multiplication:

1. 0+b = b
2. Sa+b = S(a+b)
3. 0.b = 0
4. Sa.b = b+a.b

The idea is that S represents the "successor" function maping n to n+1. Using just these definitions, and induction, we can define addition and multiplication for all natural numbers. For example, 3 is represented by SSS0 and 2 by SS0 and we can compute 3+2 using

2+3
= SS0+SSS0 by definition
= S(S0+SSS0) by 2
= S(S(0+SSS0)) by 2
= SSSSS0 by 1
= 5 by definition

But where do addition and multiplication come from? One point of view is that the natural numbers are what we get when we take finite sets but consider sets of the same size to be equal. We can do the same with finite types. The type Bool and Maybe () both have two elements (ignoring bottoms) and are isomorpic. We can just consider these to be the same type, called 2. Given two types A and B we can form Either A B. The number of elements in this new type is the sum of the number of elements in A and B. If we blur the distinction between isomorphic types we can think Either as being the addition operator. Similarly, (,) can be thought of as multiplication. The Peano axioms now describe the properties of addition and multiplication defined in this way.

When we consider different types to be equal we lose some information. In particular, we lose that fact that given two types of the same size, we can construct an explicit isomorphism between them. But there's no need to do this. We can go back to the Peano axioms and reinterpret them as a recipe for constructing the isomorphism. If we do this, then any theorem we prove (constructively) using the Peano axioms can be interpreted as explicitly constructing an isomorphism between types. We normally just forget about the isomorphism. This 'forgetting' is so common that it has a name: decategorification. Putting the structure back is called categorification.

Type Level Naturals
We will represent the natural number n as a type with precisely n elements. We'll start with the type representing zero. Obviously it must have no elements. It's traditionally called Void.

> data Void
> instance Show Void where
>     show _ = undefined

That undefined will cause no problems as we can never pass an argument into show.

If Void is playing the role of 0 we need something to play the role of S. That's Maybe. Given a type A, Maybe A is the type with one more element. So we can mimic the definitions of the natural numbers:

> type One   = Maybe Void
> type Two   = Maybe One
> type Three = Maybe Two
> type Four  = Maybe Three
> type Five  = Maybe Four
> type Six   = Maybe Five

and so on. I'll call these the natural number types. We can also label the elements of these types. Here are some elements:

> zero  = Nothing
> one   = Just zero
> two   = Just one
> three = Just two
> four  = Just three
> five  = Just four

Addition
Now we can define addition. We want to be able to take a pair of natural number types A and B and construct an explicit isomorphism between Either A B and a natural number type which I'll label Plus A B. I'll call the isomorphisms one way plus and the other way plus'. Here's a suitable type class:

> class Plussable a b where
>    type Plus a b
>    plus  :: Either a b -> Plus a b
>    plus' :: Plus a b -> Either a b

From axiom 1 we want 0+b=b. This immediately gives:

> instance Plussable Void b where
>    type Plus Void b = b

>    plus (Right b) = b
>    plus'       b  = Right b

We can view axiom 2, Sa+b = S(a+b), as:



The implementation of plus implements the mapping of the shaded square directly. If we ignore the shaded square and consider only the unshaded ones, then we are left with another simpler addition. We can implement the isomorphism for that by using plus recursively. Here's the code:

> instance Plussable a b => Plussable (Maybe a) b where
>    type Plus (Maybe a) b = Maybe (Plus a b)
>    plus (Left Nothing)  = Nothing
>    plus (Left (Just a)) = Just ((plus :: Either a b -> Plus a b) (Left a))
>    plus (Right b)       = Just ((plus :: Either a b -> Plus a b) (Right b))

>    plus' Nothing = Left Nothing
>    plus' (Just x) =
>        let i' = plus' :: Plus a b -> Either a b
>        in  case i' x of
>            Left  a -> Left (Just a)
>            Right b -> Right b

Multiplication
Now we can implement multiplication similarly. First the type class:

> class Timesable a b where
>    type Times a b
>    times  :: (a, b) -> Times a b
>    times' :: Times a b -> (a, b)

Multiplication by zero gives zero. This is straightforward to implement for the simply reason that we don't actually have to implement isomorphisms for the empty type:

> instance Timesable Void b where
>    type Times Void b = Void
>    times  _ = undefined
>    times' _ = undefined

(That's not quite true, Haskell, for some reason, forces us to write a line of code that can never be used. I think this ought to be fixed.)

> instance (Timesable a b, Plussable b (Times a b)) => Timesable (Maybe a) b where
>     type Times (Maybe a) b = Plus b (Times a b)

>     times (Nothing, b) =
>         let i = plus :: Either b (Times a b) -> Plus b (Times a b)
>         in  i (Left b)

>     times (Just a, b) =
>         let i = plus :: Either b (Times a b) -> Plus b (Times a b)
>         in  i (Right (times ((a, b))))

>     times' b = 
>         let i' = plus' :: Plus b (Times a b) -> Either b (Times a b)
>         in case i' b of
>             Left  b  -> (Nothing, b)
>             Right ab -> let (a, b) = times' ab in (Just a, b)

That's it. We've decategorified type level arithmetic. Given an equality like 2*3=6 we automatically get an isomorphism like

times :: (Two, Three) -> Six

Isomorphisms from Equations
But what about more general equation like 2*2+5 = 3*3? Can we automatically construct the isomorphism?

One approach is simply to reduce each side of the equation to its canonical form, in this case 9, and then use this to construct a pair of isomorphisms, one from the left hand side to the 9 element natural number type, and one from 9 element natural number type to the right hand side. We'll use a type class to indicate that a type can be reduced to canonical form. The map doing the reduction will be called canonical:

> class Canonicable a where
>     type Canonical a

>     canonical  :: a -> Canonical a
>     canonical' :: Canonical a -> a

Void is already in canonical form so there's nothing to do in this case:

> instance Canonicable Void where
>     type Canonical Void = Void

>     canonical  = id
>     canonical' = id

If something is of type Maybe A, and A is reducible to canonical form, then we can simply reduce Maybe A in two steps:

> instance Canonicable a => Canonicable (Maybe a) where
>     type Canonical (Maybe a) = Maybe (Canonical a)

>     canonical  Nothing  = Nothing
>     canonical  (Just n) = Just (canonical n)
>     canonical' Nothing  = Nothing
>     canonical' (Just n) = Just (canonical' n)

Now I give the rule for reducing Either A B to canonical form. We just have to reduce A and B to canonical form and then apply plus:

> instance (Canonicable m, Canonicable n, Plussable (Canonical m) (Canonical n)) => Canonicable (Either m n) where
>     type Canonical (Either m n) = Plus (Canonical m) (Canonical n)
>     canonical (Left m) =
>         let i = plus :: Either (Canonical m) (Canonical n) -> Plus (Canonical m) (Canonical n)
>         in  i (Left (canonical m))
>     canonical (Right n) =
>         let i = plus :: Either (Canonical m) (Canonical n) -> Plus (Canonical m) (Canonical n)
>         in  i (Right (canonical n))
>     canonical' x =
>         let i' = plus' :: Plus (Canonical m) (Canonical n) -> Either (Canonical m) (Canonical n)
>         in case i' x of
>             Left  m -> Left  (canonical' m)
>             Right n -> Right (canonical' n)

Now we need to do the same for multiplication. I'm beginning to feel sorry for the stress we're putting the compiler through:

> instance (Canonicable m, Canonicable n, Timesable (Canonical m) (Canonical n)) => Canonicable (m, n) where
>     type Canonical (m, n) = Times (Canonical m) (Canonical n)
>     canonical (m, n) =
>         let i = times :: (Canonical m, Canonical n) -> Times (Canonical m) (Canonical n)
>         in  i (canonical m, canonical n)
>     canonical' x =
>         let i' = times' :: Times (Canonical m) (Canonical n) -> (Canonical m, Canonical n)
>             (m, n) = times' x
>         in (canonical' m, canonical' n)

Now using the canonical forms we can build the isomorphism for any equation:

> iso :: (Canonical m ~ Canonical n, Canonicable m, Canonicable n) => m -> n
> iso m = canonical' (canonical m)

So let's return to 2*2+5=3*3. The isomorphism should be:

> test = iso :: Either (Two, Two) Five -> (Three, Three)

If we've done our job correctly, the compiler won't complain that it can't build the isomorphism.

If you really want you can try running this code for a few values:

> go1 = Left (zero, one)
> go2 = Left (one, zero)
> go3 = Right four

Try writing code to implement the inverse, checking that it does give the inverse for these three cases.

Conclusions
So there you have it, categorified arithmetic. Of course categorifying the naturals isn't so hard. But what does it mean to categorify the number π? You'll have to read some John Baez to find out more.

There sort of is an application of the operations defined above. The type Three say is the type of indices into a three element type. More generally, these natural number types give indices into fixed length containers and the addition and multiplication operations give type safe ways to map between containers that have the same size. This could be used to pack n-dimensional fixed size arrays into 1-dimensional arrays and vice-versa with compile-time checking of array indices. In practice, however, the compiler would need to be smart enough to realise it could use integers internally rather than the more complex structures it's probably using. But it's curious to see similar operations appear in some OpenCL array manipulation code I've been playing with.

The code above isn't all that pretty. As I've said before: Haskell is two languages. There's the value level language and the type level one. The former is much prettier than the latter, especially if you can use type inference to eliminate the latter.

By the way, you can view iso as a command to trigger the Haskell compiler to prove there is an isomorphism between two types of a certain class. This is very similar to what a tactic in Coq does. In fact, the code I've written above is very similar to what a proof in Coq might look like. The main difference is that Coq gives you a helping hand and can fill in details whereas Haskell forces us to do all of the work ourselves.

An Irrelevant Aside
When I was still at high school a friend returned to visit after a few months at university. He'd been playing with Prolog and showed me how to define Peano arithmetic in that language. Since then, I've sort of been obsessed with squeezing Peano arithmetic out of every computational system that can do it. Hence my C++ code here. I looked him up on the web and it turns out he also wrote the original BSD automounter. Small world.

Tagging Monad Transformer Layers

A quick post extracted from some code I was writing at the weekend.

> {-# OPTIONS_GHC -fglasgow-exts #-}
> {-# LANGUAGE ScopedTypeVariables, OverlappingInstances #-}

> import Control.Monad.Trans
> import Control.Monad.State
> import Control.Monad.Writer
> import Control.Monad.Identity

Monad transformers can get a little ugly. Here's a toy example that looks pretty bad:

> test1 :: StateT Int (StateT Int (StateT Int (WriterT String Identity))) Int
> test1 = do
>    put 1
>    lift $ put 2
>    lift $ lift $ put 3
>    a <- get
>    b <- lift $ get
>    c <- lift $ lift $ get
>    lift $ lift $ lift $ tell $ show $ a+b+c
>    return $ a*b*c

> go1 = runIdentity (runWriterT (runStateT (runStateT (runStateT test1 0) 0) 0))

There are obvious ways to make it prettier, like the suggestions in RWH. But despite what it says there, the monad "layout" is still "hardwired" and the code is fragile if you decide to insert more layers into your transformer stack. It's no way to program.

So here's an alternative I came up with. First we make a bunch of tags:

> data A = A
> data B = B
> data C = C
> data D = D

We can now label each of the monad transformers with a tag:

> test2 :: TStateT A Int (TStateT B Int (TStateT C Int (TWriterT D String Identity))) Int

And now we can have everything lifted to the appropriate layer automatically:

> test2 = do
>    tput A 1
>    tput B 2
>    tput C 3
>    a <- tget A
>    b <- tget B
>    c <- tget C
>    ttell D $ show $ a+b+c
>    return $ a*b*c

> go2 = runIdentity (runTWriterT (runTStateT (runTStateT (runTStateT test2 0) 0) 0))

Much more readable and much more robust. Change the order of the layers, or insert new ones, and the code still works.

I've tried to make this minimally invasive. It just introduces one new monad transformer that can be used to tag any other. The definitions like TStateT and tput are just trivial wrapped versions of their originals.

Anyway, this is just the first thing that came to mind and I threw it together quickly. Surely nobody else likes all those lifts. So what other solutions already exist? I'd rather use someone else's well tested library than my hastily erected solution:

> data T tag m a = T { runTag :: m a } deriving Show

> instance Monad m => Monad (T tag m) where
>    return a = T (return a)
>    T x >>= f = T $ x >>= (runTag . f)

> instance MonadTrans (T tag) where
>    lift m = T m

> class TWith tag (m :: * -> *) (n :: * -> *) where
>    taggedLift :: tag -> m a -> n a

> instance (Monad m, m ~ n) => TWith tag m (T tag n) where
>    taggedLift _ x = lift x

> instance (Monad m, Monad n, TWith tag m n, MonadTrans t) => TWith tag m (t n) where
>    taggedLift tag x = lift (taggedLift tag x)

> type TStateT tag s m = T tag (StateT s m)
> runTStateT = runStateT . runTag

> tput tag x = taggedLift tag (put x)
> tget tag = taggedLift tag get

> type TWriterT tag w m = T tag (WriterT w m)
> runTWriterT = runWriterT . runTag

> ttell tag x = taggedLift tag (tell x)

---

Target Enumeration with the Euler Characteristic. Parts 1 & 2

Part 1

> import Data.Char

A Statement of the Problem
The problem I ultimately want to solve, and its solution, is described in the paper Target Enumeration via Euler Characteristic Integrals. My goal here is to show how to implement that solution on a computer and make it accessible to a wider audience.

Suppose we have a set of targets we want to count. This could be anything from enemy tanks rolling over the plains to electronically tagged wildlife roaming the countryside. Each target has a region of influence which might simply be circular in shape, or might be more complex and depend on the target. Now suppose that we have a high density of sensors scattered over our domain and that each sensor can tell us how many regions of influence it lies in. Roughly speaking, each sensor counts how many targets are nearby. How do we compute how many targets we have in total?

Here's an illustration:





There are four targets. The region of influence for each one is coloured making it easy to see which region is which. I've labelled each region with an integer showing how many targets can be detected in that region. The idea is that we'd have a very dense scattering of sensors in our domain, each sensor reporting an integer. In effect we'd be getting an image like a rasterised version of that picture. But we wouldn't be getting the convenient colours, just an integer per pixel.

At first it seems like a trivial problem. The sensors can all count, and if every target is in range of a sensor, every target will be counted. But we can't simply add the numbers from all of the sensors as many sensors will be in the domain of influence of the same target. If we sum all of the numbers we'll be counting each target many times over. We need to be able to subtract off the targets that are counted twice. But some targets will be counted three times and so on. And how do we tell when a target has been counted twice when all we have are counts?

We'll make one simplifying assumption in solving this problem: that the regions of influence are simply connected. In other words, they are basically some kind of shape that doesn't have holes in it. That could mean anything from a square or disk to a shape like the letter 'W'. But it excludes shapes like annuli or the letter 'B'. If we make this assumption then we can solve this problem with a very simple algorithm that will work in almost all cases. In fact, the only time it fails will be situations where no algorithm could possibly work. But there's a little ground to cover before getting to the solution.

We'll make another simplifying assumption for now. That the sensors are arranged in a rectangular grid. So the data we get back from the sensors will be a grid filled with integers. That essentially turns our problem into one of image processing and we can think of sensor values as pixels. Here's a picture where I've drawn one domain of influence and I've indicated the values returned for three of the sensors.




Simple Grids
So lets assume the sensors have coordinates given by pairs of integers and that they return integer counts. The state of all the sensors can be represented by a function of this type:

> type Field = Int -> Int -> Int

We'll assume that we get zero if we try to read from beyond our domain. We can represent a grid of sensors, including the grid's width and height, using:

> data Grid = Grid Int Int Field

For efficiency something of type Field ought to read data from an array, but I'll not be assuming arrays here.

We can define display and addition of two grids:

> instance Eq Grid

> instance Show Grid where
>   show (Grid w h f) = concat
>       [[digit (f x y) | x <- [0..w-1]] ++ "\n" | y <- [0..h-1]]
> digit 0 = '.'
> digit n = chr (48+n)

> instance Num Grid where
>   Grid w0 h0 f0 + Grid w1 h1 f1 = Grid
>       (w0 `max` w1) (h0 `max` h1)
>       (\x y -> f0 x y + f1 x y)

Ourultimate goal is to define some kind of count function with signature:

count :: Grid -> Int

Now suppose the function f gives the counts corresponding to one set of targets and g is the count corresponding to another. If the region of influence of these two sets of targets is separated by at least one 'pixel' then it should be clear that

count f + count g == count (f+g)

So at least approximately, count is additive. We also need it to be translation invariant. There's only one function that has these properties, summing up the values at all pixels:

> gsum (Grid w h f) = sum [f x y | x <- [0..w-1], y <- [0..h-1]]

We can implement functions to make some example grids:

> point x y = Grid  (x+1) (y+1)
>                   (\x0 y0 -> if (x0, y0) == (x,y) then 1 else 0)
> circle x y r = Grid  (x+r+1) (y+r+1)
>                      (\x0 y0 -> if (x-x0)^2+(y-y0)^2<r^2 then 1 else 0)

And now we can build and display some examples:

> test1 = circle 10 10 5+circle 7 13 4+point 5 5+point 9 12
> test2 = gsum test1

Here's a typical output:

*Main> test1
................
................
................
................
................
.....1..........
........11111...
.......1111111..
......111111111.
......111111111.
.....1222211111.
....11222221111.
....11222321111.
....1112222111..
....111122211...
....1111111.....
.....11111......
................
*Main> test2
116

It should be pretty clear that this doesn't count the number of targets. So how can we implement something additive and yet count targets?

Another operation we can perform on grids is scale them. Here's an implementation of scaling a grid:

> scale n (Grid w h f) =  Grid (w*n) (h*n)
>                         (\x y -> f (x `div` n) (y `div` n))

Scaling up an `image' shouldn't change the number of targets detected. It should only correspond to the same number of targets with double-sized regions of influence. So we'd also like the following property:

count (n `scale` f) = count f

It's easy to see that that gsum actually has the following property for n>0:

gsum (n `scale` f) = n^2 * gsum f

(^ is the power function. For some reason lhs2TeX displays it as an up arrow.) These requirements are pretty tough to meet with an additive operation. But there's an amazing transformation we can perform on the data first. Instead of working on a grid with one value for each pixel we'll also store values for the 'edges' between pixels and for the 'vertices' at the corners of pixels.

Euler Grids
So lets define a new kind of grid to be a tuple of Field functions, one for faces (ie. the pixels), one for horizontal edges, one for vertical edges, and one for vertices.

> data EGrid = EGrid {
>       eWidth::Int, eHeight::Int,
>       faces::Field, hedges::Field, vedges::Field, vertices::Field
>   }

The lower left vertex is (0,0) but we need to add an extra row and column of vertices on the right. Similarly we'll need an extra row and an extra column of edges. We can now `resample' our original grid onto one of these new style grids:

> g2e (Grid w h f) = EGrid w h
>   f
>   (\x y  -> f (x-1)  y      `max` f x y)
>   (\x y  -> f x      (y-1)  `max` f x y)
>   (\x y  -> f (x-1)  (y-1)  `max` f (x-1) y `max` f x (y-1) `max` f x y)

I'm using the rule that the value along an edge will be the maximum of the values on the two impinging faces. Similarly, the vertices acquire the maximum of the four faces they meet.

I'll try to illustrate that here:



I hope you can see, from the placement of the labels, how I've attached values to edges and vertices as well as faces.

We now have a bit more freedom. We have three different types of sum we can carry out:

> fsum (EGrid w h f _ _ _) = gsum (Grid w h f)
> esum (EGrid w h _ e f _) = gsum (Grid (w+1) h e)+gsum (Grid w (h+1) f)
> vsum (EGrid w h _ _ _ v) = gsum (Grid (w+1) (h+1) v)

(We could sum over horizontal and vertical edges separately too, but if we did that then a 90 degree rotation would give a different target count.)

Now we can define a measurement function that takes three `weights' and gives us back a weighted sum:

> measure a b c g = let e = g2e g in a*vsum e+b*esum e+c*fsum e

We can reproduce the gsum function as

> area = measure 0 0 1

Try some examples to test that

((n^2) *) . area = area . scale n

Exercises
Now I can leave you with some challenges:

1. Find some suitable arguments a, b and c to measure so that we get:

mystery_property1 = measure a b c
((n^1) *) . mystery_property1 = mystery_property1 . scale n

I'll let you assume that there is some choice of values that works.

(Hint: you just need to try applying mystery_property1 to a few scalings of some chosen shape. You'll quickly find some simultaneous equations in a, b and c to solve. Solve them.)

2. Can you find a simple geometric interpretation for mystery_property1? Assume that the original input grid simply consists of zeros and ones, so that it's a binary image. It shouldn't be hard to find a good interpretation. It's a little harder if it isn't a binary image so don't worry too much about that case.

3. Now find some suitable arguments a, b and c to measure so that we get:

mystery_property2 = measure a b c
((n^0) *) . mystery_property2 = mystery_property2 . scale n

4. Can you find a simple interpretation for binary images? You might think you have it immediately so work hard to find counterexamples. Have a solid interpretation yet? And can you extend it to images that consist of more general integers?

5. Optimise the code for mystery_property2 assuming the image is binary and that the input is on a 2D array. Ultimately you should get some code that walks a 2D array doing something very simple at each pixel. Can you understand how it's managing to compute something that fits the interpretation you gave?

6. Define a version of scale called gscale that works on EGrids. Among other things we should expect:

g2e . (scale n) == gscale n . g2e

and that the invariance properties of the mystery properties should hold with respect to gscale.

I'll answer most of these questions in my next post. If you find mystery_property2 you've rediscovered one of the deepest notions in mathematics. Something that appears in fields as diverse as combinatorics, algebraic geometry, algebraic topology, graph theory, group theory, geometric probability, fluid dynamics, and, of course, image processing.

---

Part 2
The Semi-perimeter
Let's start with exercise 1 from my previous post. I allowed you to assume there was a solution. Knowing this we only need to try scaling up one shape. Here are three scalings of a single pixel:





For the 1x1 square we have: 1 face, 4 edges, 4 vertices.
For the 2x2 square we have: 4 faces, 12 edges, 9 vertices.
For the 3x3 square we have: 9 faces, 24 edges, 16 vertices.

(If you don't feel like counting these yourself you can use code like:

measure 1 0 0 $ point 0 0
measure 1 0 0 $ 2 `scale` point 0 0
measure 1 0 0 $ 3 `scale` point 0 0

)
and so on.

So now we have some equations:

 4a+ 4b+ c =  x
 9a+12b+4c = 2x
16a+24b+9c = 3x

We find a=0, c=-2b. I'll pick b = 1, c = -2. There's a straightfoward interpretation of measure 0 1 (-2) for binary images. It computes half of the perimeter, the semi-perimeter. There's a nice way to see this. We can try to count the number of edges in a shape starting from the number of faces in it. You can think of each face as being surrounded by 4 "half-thickness" edges. Where two faces meet we get a full thickness edge, so using half the number of faces counts internal edges correctly. But around the border of a shape we are left with contributions from only a face on one side. So we're only counting the perimeter edges by half. We get a shortfall of the semi-perimeter. Working backwards tells us how to compute the semi-perimeter from the total number of edges and faces.

For more general images, not just binary ones, we can roughly think of measure 0 1 (-2) as computing the sum of the semi-perimeters of all of the isocontours of our image.

The Euler Characteristic
The interesting case is now exercise 3. This time our equation is:

 4a+ 4b+ c = x
 9a+12b+4c = x
16a+24b+9c = x

Now we get a solution a = 1, b = -1, c = 1. If you try computing this for a few shapes it looks like it's counting the number of connected components of a binary image. However, once you realise the possibility of some holes in your image you find that it always turns out to be the total number of connected components minus the total number of holes.

Here's an example. Treat this as a single complete image:




It has two components and one hole so we expect measure 1 (-1) 1 to give us 1. We can count:

14 faces
42 edges
29 vertices

Giving 29-42+14 = 1.

I'm not sure which is the easiest proof that vertices-edges+faces counts the number of components minus the number of holes. One approach is this: we only need to consider one connected component at a time. Remove its holes. Now build up the shape one pixel at a time starting with one pixel and ensuring that you have a hole-free shape at each stage. It's not hard to enumerate all of the possible ways you can add one pixel to an existing shape and show that each such addition leaves measure 1 (-1) 1 unchanged. If you now make a single pixel hole in your shape you'll see that it lowers measure 1 (-1) 1 by 1. If you now continue to add pixels to the hole, in a way that doesn't change the number of holes, you'll see that measure 1 (-1) 1 remains unchanged again.

measure 1 (-1) 1 computes what is known as the Euler characteristic of a shape. I talked a little about this in one context earlier and showed how to compute it in another context here. The Euler characteristic is a topological invariant of a shape in the sense that a rubber sheet deformation of a shape leaves the number of holes and the number of components unchanged.

The above description shows that the Euler characteristic is particularly easy to compute. It simply requires a map-reduce operation over the entire grid. But what about the separate terms: the number of components and the number of holes? These seem like simpler notions and you might expect them to be just as easy to compute. Actually they are harder to compute. Compare also with flood fill algorithms which solve a related problem. Minsky and Papert show in their book Perceptrons that any topological invariant that can be learnt by a one layer neural network (with certain resonable restrictions) must be a function of the Euler characteristic. I find it quite amazing that this notion from topology is connected (no pun intended) to learnability.

We can define

> euler = measure 1 (-1) 1

I have sketched an argument that euler counts #components-#holes. If we assume that each of our connected components has no holes then it counts the number of components. But here's a neat thing: if we *add* two images that contain strictly overlapping shapes (ie. not just touching each other along their boundaries) then because of additivity, euler will still count the number of shapes. In other words, if you did the exercises then you solved the target enumeration problem. It's pretty miraculous. You could splat down thousands of geometric shapes into an image. They can overlap as much as you like. But as long as they don't touch along a boundary you can still compute the total number of shapes. If two shapes do touch along a common boundary then no algorithm can work, after all they'll be indistinguishable from a single connected shape. For a quick example, notice how

> test3 = euler test1

recovers that test1 is the sum of 4 shapes that don't touch.

Generalised Integrals
Consider our original measurement area. This sums the values at each pixel. It is a numerical approximation to the integral of a function sampled at each pixel. Likewise each of the measure functions is numerical approximation to a generalised type of integral. The original paper uses these integrals to solve its problem. I have simply used a discrete version.

I apologise for only sketching proofs. It takes considerably more work to provide rigorous proofs. But I encourage you to experiment with the code and attempt to find counterexamples. The history of the Euler characteristic is itself characteristed by a kind of back and forth between attempted proofs and counterexamples than in a strange way mirrors the innocent looking definition: vertices-edges+faces.

By the way, some people have propsed that the Euler characteristic is a kind of generalisation of the idea of counting. It shares many properties with the usual notion of cardinality.

One last thing: I have implicitly shown that target counting is learnable by a certain type of one-layer neural network.

And thanks to @alpheccar for pointing out the original target enumeration paper.

---

Historical Addendum
I'm repeating this as a possibly apocryphal story I have heard from other parties: Minsky and Papert demonstrated that the only learnable topological invariants for single layer network are functions of the Euler characteristic. In particular, they demonstrated the unlearnability of connectedness. This was a precisely stated no-go theorem that discouraged and slowed investment in neural network research for many years and helped contribute to the AI Winter. Is there a good historical book on this period of AI research?

---

Counting Targets using the Euler Characteristic, Part 1

Using LaTeX was so much easier than HTML that I'm doing it again. The PDF is here and the literate Haskell source is here.

Abstract

The problem I ultimately want to solve, and its solution, is described in the paper Target Enumeration via Euler Characteristic Integrals by Baryshnikov and Ghrist. My goal here is to show how to implement that solution on a computer, and by doing so make it accessible to a wider audience.

Thanks to @alpheccar for tweeting about the original paper.

---

Monads are Trees with Grafting

This is an attempt to collect together, in tutorial form, a few of the things I've said about monads going back as far as my field guide. It's probably not a good first tutorial, but it contains things I'd wish I'd known immediately after reading my first tutorial.

Unfortunately, I used a few more LaTeX features than I ought to have while drafting it making it hard to get back into HTML form cleanly. So here is a PDF: Monads are Trees with Grafting.

Probably the main inspiration was this post by Derek Elkins.

---

Where do monads come from?

Introduction
Monads play an important part in implementing computational "effects" in Haskell, and yet monads come out of pure mathematics where they are used to describe operations in algebraic structures. What do these two applications have to do with each other? In this post I hope to answer that question and maybe give another view on monads.

> {-# LANGUAGE GADTs #-}

Algebra for Flowcharts
You might not have seen one of these in a long time. A flowchart describing a subroutine:



For now I'm considering flowcharts only with nodes to set or reset just one flag, test the flag, and return values. I'm not allowing any paths to reconverge. So these flowcharts are trees describing the possible paths of control flow.

We can describe such flowcharts with a more compact notation. Use x_n to mean:



Think of x_n being an variable of unknown value.
Use the + operator to mean tests of the flag. So x+y represents:



Note that this '+' isn't commutative because the left branch will always mean 'no' and the right branch will mean 'yes'.

And we'll use the symbols '0' and '1', applied on the left, to mean these two nodes respectively:



The flowchart above can now be written compactly as 1x₃+0(x₄+x₅). (Note that I'm overloading the notation for addition and multiplication, this isn't intended to be ordinary addition and multiplication of numbers.)

It's clear that the 'return 5' never gets executed. So the above flowchart is equivalent to this one:



The rule describing this simplification can be expressed with our compact notation as the rule:

0(a+b) = 0a
1(a+b) = 1b

So we already have a strong clue that there is a connection between algebra and computation. Here are some more rules that you can check for yourself:

abc = bc (here a and b are each either 0 or 1)

a+(b+c) = a+c = (a+b)+c

a+a = a = 0a+1a

Substitution
Now suppose we wish to extend this subroutine so it does a bit more work. Instead of simply returning n we might want it to set a flag according to whether the previous result is even or odd, and subtract 1 if it is odd. It's a little awkward because I've insisted that these trees don't rejoin so we need to patch the tree in three places using the rule:



Note that the descriptions of the two cases and the n here aren't part of the flowchart, they're part of the language I'm using to talk about flowcharts.

We can do something similar for our algebraic expressions. We can define a mathematical function taking integers to algebraic expressions such that

f(n)
= 0x_n if n is even
= 1x_n-1 if n is odd

and now define an operation * by

f*u = u with x_n replaced by f(n) throughout

Our updated flowchart can now be represented by f*(1x₃+0(x₄+x₅)). Now you might see why I used the weird x_n notation. The * operation is ordinary substitution of an expression for a variable that we learn in high school, so I used x_n notation to make it clear that x_n is a variable. You can think of f as a dictionary saying how the values of the variables are to be expanded.

Monads
Let TN be the set of expressions in our algebra where the subscripts of the x's come from the set N. Then we have the functions:

η:N→TN defined by η(n)=x_n
and
μ:T(T(N))→T(N) defined by μ(a)=i*a, where i is the identity.

We have a monad. (Exercise: check the monad laws.) Note that this definition isn't specific to our particular algebra. This definition applies to a very broad class of algebraic expression. Monads in mathematics are used to describe substitution of subexpressions into expressions. But I hope you can now see why they have a role to play in computing: the act of substituting an expression for a variable is essentially the same thing as grafting more operations on to the end of our tree of computations.

But if you look at the definitions of typical monads in the Haskell libraries, they don't look much like trees, and the implementation of >>= doesn't look much like grafting trees. The rest of this post will be about working through a specific example to show that it really is about tree grafting after all.

A Type for Expression Trees
We'll need to define a type to represent our algebraic expressions in Haskell. The n is going to be used to label our variables, and s will be explained shortly:

> data Expr s n where

Now we need constructors to represent 0 and 1. We can simplify this by combining them into one constructor and using the type Bool to signify whether it is 0 or 1 that we want. But there's no need to restrict ourselves to Bool so I'll make the code more general and use s instead:

>    Put :: s -> Expr s n -> Expr s n

So 0 is False and 1 is True and `Put` is our multiplication operator. (Maybe you've already guessed why I called it that.)

The next constructor will be for addition. I'm going to give it the weird name Get. The obvious thing to write would be:

    Get :: Expr s n -> Expr s n -> Expr s n

We could uncurry that instead and use

    Get :: (Expr s n, Expr s n) -> Expr s n

But we know that a pair (a, a) is the same as a function Bool -> a. So I'm going to actually use:

>    Get :: (s -> Expr s n) -> Expr s n

So instead of giving us a binary node in a tree, Get gives us a possible infinite s-ary node. But we'll be working with s = Bool for now.

Now we need variables. We'll just use:

>    Var :: n -> Expr s n

This is a recursive type so we can write a fold for it. If we'd been using generic library we'd get the fold for free but I'm giving it explicitly:

> fold :: (s -> w -> w) -> ((s -> w) -> w) -> (n -> w) -> Expr s n -> w
> fold put get var (Put s x) = put s (fold put get var x)
> fold put get var (Get c) = get (\s -> fold put get var (c s))
> fold put get var (Var n) = var n

We can implement the monad for substitution:

> instance Monad (Expr s) where
>   return n = Var n
>   x >>= f = fold Put Get f x where

Note how the fold recursively replaces Get with Get, Put with Put and Var n with f n.

By now the game's up. If you read my post from last week then you know that I'm building something like the State monad. But although this monad has a similar interface, internally it's nothing like State. For one thing, Expr is a tree type and internally it could end up with an arbitrarily large number of nodes. How can we extract the familar State monad from this?

Canonical Form
Expressions in various algebraic languages can often be reduced to a canonical form. For example, we have disjunctive normal form for (a fragment of) logic. Expressions in our language above can be reduced to a canonical form. In the language of flowcharts, every one of our flowcharts can be reduced to something of the form:



Algebraically the claim is that any expression reduces to something of the form au+bv. In the general case the items a and b are of type s and u and v can be represented by elements of n. So this is something of type ((s, n), (s, n)). Again we replace the pair by a function to get the type s -> (n, s). So we want a function to convert to canonical form:

> canonical :: Expr s n -> (s -> (n, s))
> canonical = fold put get var where

We can get a return box into canonical form like this:



Algebraically this is x_n→0x_n+1x_n and we can implement it as

>   var n = \s -> (n, s)

Now we need to simplify the application of 0 or 1 to something already in canonical form. 0(au+bv) = 0au = au and 1(au+bv) = 1bv = bv. Essentially it just picks the first or second term in the sum according to whether we're multiplying on the left by 0 or 1. But that's exactly what the function type in s -> (n, s) does. So we can just write:

>   put s c = \t -> c s

Now we need to consider sums of two terms in canonical form. We can derive this just by looking at the sum of two sums: (p+q)+(r+s). This is just p+s. We get this by picking the left branch twice or the right branch twice, so we're just applying our function twice to the same thing:

>   get c s = c s s

Comparing with my last post you should see that canonical is what I called an interpreter. So the act of interpreting these tree expressions is one and the same thing as reducing them to a canonical form. This shouldn't be so surprising - executing a functional program is itself a process of reduction to a normal form. Here's an example of using the above reduction to get our flowchart executed:

> test1 =
>   Get $ \c -> if c
>       then Put False $
>       Get $ \d -> if d
>           then Var 5
>           else Var 4
>       else
>           Put True $
>           Var 3
> go1 = canonical test1 False

canonical returns a different tree to the original tree. But we can easily map back to trees again:

> inject :: (s -> (n,s)) -> Expr s n
> inject f = Get $ (\(n, s) -> Put s (Var n)) . f

So the composition inject . canonical takes a tree and reduces it to canonical form keeping it as a tree.

We can view things like this: the State s n type is just the Expr s n type modulo an equivalence relation. The usual State type is a tree but this fact is obscured because we only work with a representative of each equivalence class, not entire trees. But the Monad instance is still simply tree grafting. We can prove this as follows:

> data State s n = State { runState :: s -> (n, s) }

> instance Monad (State s) where
>   return n = State $ canonical $ Var n
>   State f >>= g = State $ canonical $ inject f >>= (inject . runState . g)

Note how the implementations consist solely of isomorphisms or wrapping/unwrapping. The type State is simply inheriting the Monad instance from the tree type. We can now put Put and Get into a form usable in the monad. Again, we're just wrapping and "canonicalising" near-trivial operations.

> get = State $ canonical (Get (\s -> Var s))
> put n = State $ canonical $ Put n (Var ())

> test2 = do
>   c <- get
>   if c
>       then do
>           put False
>           d <- get
>           if d
>               then return 5
>               else return 4
>       else do
>           put True
>           return 3

> go2 = runState test2 False

Summary
Like in the previous post We're starting with signatures for functions. We can form an expression tree type. With the substitution operation this naturally becomes a monad.

We also specified equations for our functions. This leads to the idea of a reduced form for our expression trees. The state monad (and, in fact, most of the standard monads) are examples of reduced forms. If we implement canonical for an algebra, the monad follows for free. More precisely: given a signature and equations, the tree type and monad can be derived automatically. It seems tricky to derive the rules for reducing expressions to canonical form, but once you've done that the monad for the reduced form can be derived automatically.

The act of evaluating a monadic expression is the same as reducing an expression tree to its canonical form.

Various Thoughts
Given two algebraic theories with their own signatures and equations they can be combined together. One obvious way is to simply "freely" throw all of the signatures and equations together. This gives the "sum" of the two theories. If we additionally impose some new equations, a certain type of compatibility condition between the theories, then we get the product. So starting with two theories we get the monads for the two theories and the monad for the sum or product. This might provide an alternative approach to monad transformers for combining effects. (Actually, I'm doubtful of this plan for automatically combining effects because surprising things can happen, but I still think there's something useful to be extracted from the theory.)

Could this form the basis for more monad tutorials? Programmers from imperative, functional and other backgrounds are familiar with the ideas of trees and substitutions.

By the way, it's sometimes useful to implement a monad by building a tree first and then interpreting it separately. That's what I did a few posts ago and I found it much easier than deriving a complete monad in one go.

These trees are related to the diagrams I drew here.

You can repeat much of what I said above with all kinds of algebraic structures. My ICFP talk was on what you can do with the monad that arises from vector spaces and a spent a slide or two on the rules for reduction to canonical form.

An earlier post of mine on substitution.

You can interpret the canonical form of the expressions above as saying that when you have a bunch of switches (in certain types of code) you can replace them with just one. This is reminiscent of the old folk theorem about loops. In fact, the algebra I describe above is a lot like Kleene Algebra with Tests which can be used to prove that result.

This was all motivated by this post on Lambda the Ultimate. Some comments by Derek Elkins (that he's probably long forgotten) and pieces of this by Conor McBride helped me write the code.

I originally planned to say more about forming sums and products of theories but the code is going to need more work...

Update: Following one of the links I gave above I just found Unimo which is closely related.

---