Build Yourself a Bluetooth Controlled Six-Legged Robot

Introduction

I don't know why you can't easily find toy robots that can be controlled by Bluetooth. But that lack can be remedied by a few hours work. So you can see exactly what I mean, here is the robot I'm going to describe being controlled by me via a program running on a Mac:

As it's controlled from a much more powerful computer you can use sophisticated algorithms to direct it. And if you don't want to do that, it can easily be programmed to run autonomously. Either way, you can make use of a pair of infrared proximity sensors as input to your algorithms.

The first thing I must do is give credit where it's due. This is a modification of someone else's design. In fact, it's a sample project at Pololu Electronics and Robotics modified by addition of an off-the-shelf Bluetooth board.

You'll need:

1. Everything required to build the Pololu sample project. That includes Windows. I was developing on a Mac so I made use of VMWare. I can't vouch for any other Windows emulation because you'll need a USB port and I've experienced troubles with emulation and USB in the past. You'll only need Window for initial configuration of the robot. After that you can use the development platform of your choice. The list of parts is here. Order the partial kit for the motor controller. You'll need the header pins it comes with.
2. One BlueSMiRF Gold Bluetooth modem from Sparkfun.

The first step is to build the Pololu project exactly as described. At this point, you'll have a complete robot that can run autonomously. You can program it from Windows using Pololu's scripting language described in the manual.

Note that the battery used is this. Its connector is different to that described in the project. You simply need to solder a pair of pins rather than a special type of connector:

By the way, there is a minor flaw in Pololu's design. The sensor and servos are glued to the battery using hot glue. When the battery is charged it heats. You can guess the rest. So once everything works, you may want to switch to a different glue. I used epoxy.

The strands of wire from the servos are very fine. At one point I think some came astray causing a short-circuit. I was surprised the controller board still worked after I saw smoke rise up from it. So after that I used some dabs of hot glue to act as insulation around my soldering.

I found that I couldn't use the motor controller via an external USB hub, I always had to use a cable directly from my Mac.

In the process of making the robot you cut the connectors off the servos. Keep two of these with 2 or 3 inches of wire attached. We'll reuse them later.

Setting up Bluetooth

In its default mode, the robot is scriptable and controllable via the USB port. However, two pins on the board are attached to a serial port (running at 5V, not the usual RS232 voltage). These can be connected directly to a pair of pins on the BlueSMiRF board.

First we'll get the modem working independently of the robot. So go ahead and solder 6 header pins to the board like so:

Now we can connect it directly to the robot's battery using its connector to make the following circuit:

Conveniently, VCC and GND are next to each other, just like in the battery connector. The little red LED marked Stat should start flashing. This means it's ready to go. Take care not to reverse the polarity (this isn't Star Trek) or make an off by one pin error.

Now go ahead and pair it with your main computer. On the Mac it's the usual process. It'll probably appear as a device called Firefly-something. There have been various revisions of documentation for BlueSMiRFs. I guessed, correctly, that the correct passkey was 1234.

Once it's connected you should get a new serial port. Mine appeared as /dev/tty.FireFly-DB29-SPP. SPP is Serial Port Profile. Now you need a terminal application that can talk to a serial port. On Unix machines you can use screen followed by the path to the dev file. You'll find help on the web for other OSes. On Windows I expect you'll get a new COM port. It didn't seem to matter what baud rate I chose at this stage. I guess that SPP presents an interface that looks like a serial port but that the connection rate doesn't really mean anything. But I may be wrong and you may need to experiment and/or read documentation. It may take several attempt to successfully connect. (It likes to prove it has the power to say no if it wants.) When you do, you'll see the Conn LED light up green.

If you connect within one minute of powering up the BlueSMiRF you can configure it by entering its command mode. One revision of the documentation says you should type +++. My version used $$$. You won't see your keys echoed on the screen but it should reply CMD. Now type d and return to get a summary of its status. Now type SU,19200 (and return) to set the baud rate on the serial port side. It'll go much faster than this but it's worth being cautious to start with. It acknowledges this but it doesn't actually change speed until you restart the device (eg. by power cycling it). Quit command mode by typing ---. At this point, anything you type in the terminal will come out at 19200 baud on the TX pin. Quit screen by hitting ^A^X.

I'm not sure how you can easily test this short of completing the robot. I tested it using a program on a MSP430 Launchpad. I wish I had a fancy diagnostic device to read serial transmissions. I ought to build one.

By the way, talking to the BlueSMiRF as if it's a serial port is a form of backward compatibility. It should be possible to talk to it directly via the RFCOMM API. Among other things, this would allow better diagnostics.

Configuring the Robot

Now you need to configure the robot to use 19200 baud on its TX/RX pins. Connect via USB to the Windows Control Center application (like you've done before if you built the original design) and set it to use the UART at 19200 baud. Make sure CRC is off. Like so:

Once that is done you don't need Windows again unless you change the baud rate.

Connecting Modem to Robot

Now you need to connect the robot to the modem. You have 6 pins on the modem board so you need a 6 pin connector. I glued together the two three pin connectors I mentioned above to make a connector like this:  

The two outer wires can be soldered together because we're not using flow control, and so the BlueSMiRF is clear to send (CTS) whenever it is ready to send (RTS), The ends of the other four wires can be soldered directly to the motor controller board to make this circuit:

Now comes the iffy bit. I hot glued the BlueSMiRF board directly to the side of the servos like this:

That's iffy because it places electric motors directly next to a device using RF. Motors emit lots of RF noise. But they can't be that bad, after all servos have been used for years on radio controlled models, and Bluetooth probably has some error correction. It seemed to work for me. (But see note below.)

Talking to the Robot

Now we need to talk to the robot. Just about every development platform has a way of talking to serial ports. I used Haskell with the serialport package. As I had fixed the baud rate, I used the Pololu compact protocol. As the focus in this article is on hardware, I'll just give the code at the end. Once it's running you should be able to control the robot using

• f - take a step forward
• b - take a step back
• l - turn left
• r - turn right
• e - read the state of the two proximity sensors

Now you're free to code up whatever you want. I experimented with localization using a method I learnt from Eric Kidd.

Notes

I occasionally had dropped bytes. I haven't yet worked out a way to 100% reliably recover when this happens and I'm not sure if it's due to the RF noise I mentioned above. You can probably mess with the timeout parameter in the Pololu Control Center to ensure the robot returns to a known state if communications cease for a bit.

The Apple SPP driver is pretty crappy. It generally works (apart from the annoying need to make multiple attempts to connect), but if you interrupt your code at the wrong time with ^C (say) you can find it locks up so hard that you end up with an unkillable zombie process. The OS didn't even seem able to kill it on shutdown (maybe it would have eventually timed out) and I had to power down the Mac if I wanted to use the robot again. It only happened a couple of times in two weeks of intensive use, but it's annoying.

And if the robot stops working - remember these batteries are pretty small and are doing a lot, so you may just need a recharge.

Code

The code uses two threads, a server that talks to the serial port and a client for users. This means that independent threads can control the legs and eyes, say, without the serial port transactions becoming interleaved. This isn't very polished but should be good enough to start experimentation. I've only tested it under Mac OS X but it looks platform independent to me.

> module Robot where

> import Prelude hiding (Left, Right)
> import Control.Concurrent
> import Control.Monad.Trans
> import Control.Monad.Trans.Maybe
> import Control.Monad.Trans.Reader
> import Control.Exception
> import System.Hardware.Serialport


Server side using Pololu compact protocol.


> sendByte :: Int -> SerialPort -> IO ()
> sendByte b s = sendChar s (toEnum b)


> getByte :: SerialPort -> IO (Maybe Int)
> getByte s = do
>     b <- recvChar s
>     return $ fmap fromEnum b


> setTarget :: Int -> Int -> SerialPort -> IO ()
> setTarget channel target port = do
>     sendByte 0x84 port
>     sendByte channel port
>     let (b1, b0) = target `divMod` 0x80
>     sendByte b0 port
>     sendByte b1 port


> getPosition channel port = runMaybeT $ do
>     liftIO $ sendByte 0x90 port
>     liftIO $ sendByte channel port
>     lo <- MaybeT $ liftIO $ getByte port
>     hi <- MaybeT $ liftIO $ getByte port
>     return $ lo + 0x100*hi


> data SerialCommand = SetTarget Int Int
>                    | GetPosition Int (Maybe Int -> IO ())
>                    | End


> serialExec command port = 
>     case command of
>         SetTarget channel target -> setTarget channel target port
>         GetPosition channel continuation ->
>              getPosition channel port >>= continuation
>         otherwise -> return ()


> serialThread commandChannel port = do
>     command <- liftIO $ readChan commandChannel
>     serialExec command port
>     case command of
>         End -> do
>             closeSerial port
>             return ()
>         otherwise -> serialThread commandChannel port

Client side

> data Servo = Servo { channel :: Int, loLimit :: Int, hiLimit :: Int }


> data Direction = Left | Right | Forward | Back


> start tty = try (openSerial tty defaultSerialSettings
>        { baudRate = B19200 }) >>=
>     either
>         (\ex -> do
>             print (ex :: IOException)
>             threadDelay 250000
>             start tty)
>         (\port -> do
>             commandChannel <- newChan
>             forkIO $ serialThread commandChannel port
>             return commandChannel)


> type RobotM a = ReaderT (Chan SerialCommand) IO a


> writeChan' = flip writeChan


> end = ReaderT $ writeChan' End


> setServo limit servo = ReaderT $ writeChan' $
>    SetTarget (channel servo) (limit servo)


> setLo = setServo loLimit
> setHi = setServo hiLimit


> readEye channel cont = ReaderT $ writeChan' $ GetPosition channel cont


> delay = 100*1000


> raise Left    = setLo midLegs
> raise Right   = setHi midLegs
> forward Left  = setLo leftLegs
> forward Right = setHi rightLegs
> back Left     = setHi leftLegs
> back Right    = setLo rightLegs


> opposite Left    = Right
> opposite Right   = Left
> opposite Forward = Back
> opposite Back    = Forward


> move Forward = forward
> move Back    = back


> halfCycle side direction0 direction1 = do
>     raise side
>     liftIO $ threadDelay delay


>     move direction0 Left
>     move direction1 Right
>     liftIO $ threadDelay delay


> motion direction0 direction1 = do
>     halfCycle Left  direction0 direction1
>     halfCycle Right (opposite direction0) (opposite direction1)


> walkCycle    = motion Forward Back
> reverseCycle = motion Back Forward
> turn Right   = motion Forward Forward
> turn Left    = motion Back Back


> withRobot tty cmds = do
>     commandChannel <- start tty
>     flip runReaderT commandChannel (cmds >> end)


> eyes = do
>     readEye 4 $ print . (("Left eye  = " ++) . show . fmap (<512))
>     readEye 3 $ print . (("Right eye = " ++) . show . fmap (<512))


> commandLoop = do
>     key <- liftIO $ getChar
>     case key of
>         'f' -> walkCycle
>         'b' -> reverseCycle
>         'l' -> turn Left
>         'r' -> turn Right
>         'e' -> eyes
>         otherwise -> return ()
>     if key=='q'
>         then return ()
>         else commandLoop 


> testRobot = do
>     withRobot tty $ do
>         liftIO $ print "Ready"
>         commandLoop


Here's some stuff for you to configure. In particular, set the limits of the servos so that the legs don't whack into the body of the robot, stripping the gears. These numbers will be a little different depending on exactly how you built the robot. You may want to start with the upper and lower limits closer to 6000, the middle of the range of the servo motion.


> tty = "/dev/tty.FireFly-DB29-SPP"


> rightLegs = Servo 0 5000 6500
> midLegs   = Servo 1 5000 6500
> leftLegs  = Servo 2 5000 6500


One Last Thing
Double check everything I say above. I could have easily made a mistake. In particular, make sure that I haven't inadvertently introduced short circuits by comparing my diagrams against the online documentation for the parts.

---

Generalising Gödel's Theorem with Multiple Worlds. Part III.

Fixed Points
A number of branches of mathematics have some sort of implicit function theorem. These guarantee that we can define functions implicitly, rather than explicitly. For example, we can define the function f from the positive reals to the positive reals by the relation f(x)^2 = x. In this case we can write an explicit formula, f(x) = x^1/2, but the implicit function theorem gives quite general conditions when the implicit equation defines a function uniquely, even when it is very hard to write an explicit formula. We have something similar in the theory of datatypes. We define the list type implicitly in terms of lists: L(X) = 1+X L(X). There is a unique smallest type (and unique largest type) that satisfies this equation and we use such equations in Haskell programs to uniquely specify our types. (Haskell uses the largest type.)

Provability logic also has a kind of implicit function theorem. Consider the following function:

> f0 p = Neg (Box p)

Suppose the equation p <-> f0 p had a solution. Then we would have a proposition that asserts its own unprovability. In fact, we know such a proposition: the Godel sentence from Godel's first incompleteness theorem. However, the Godel sentence uses a clever Quining technique to make a proposition that makes assertions about its own Godel number. We can't express such a thing in the language of provability logic. However, consider Godel's second incompleteness theorem. This tells us that if PA can prove its consistency then it is in fact inconsistent. Or to turn it around, if PA is consistent, then PA can't prove its consistency. Asserting that PA is consistent is the same as ◊⊤ (because asserting that ⊤ is consistent with the rest of PA is the same as saying the rest of PA is consistent). So Godel's second incompleteness theorem can be written in provability logic as ¬◻(◊⊤)→◊⊤. Using the methods of part one we can show that ◊⊤→¬◻(◊⊤) is also valid. So Dia T is the solution to p <-> f0 p:

> test0 = valid $ let p = Dia T in p <-> f0 p

Considered as a data structure, notice how f0 places its argument inside a Box. Consider just the functions Prop -> Prop that (1) simply copy p into various places in its return value (ie. that don't analyse p in any way) and that (2) put all copies of its arguments somewhere inside a Box (possibly in a deeply nested way). We'll call these 'modalising' functions. Then there is an amazing result due to Solovay:

The Fixed Point Theorem
If f is modalizing then it has a unique fixed point in provability logic. (Unique in the sense that if p and q are fixed points for f, then p <-> q is valid.

Our goal in the next post in this series will be to implement an algorithm to find these fixed points.

Think about what this means. We can make up any old function that has this property and find a corresponding Godel-like theorem. For example, pick

> f1 p = Box p --> Box (Neg p)

This has a fixed point p = Box (Box F) -> Box F which we can test with

> test1 = valid $ let p = Box (Box F) --> Box F in p <-> f1 p

Or in English (quoting verbatim from Boolos) "a sentence is equivalent to the assertion that it is disprovable-if-provable if and only if it is equivalent to the assertion that arithmetic is inconsistent if the inconsistency of arithmetic is provable". (That's a good argument for using provability logic instead of English!) These are the generalisations of Godel's theorem in my title.

But there are many other important consequences. The solution to this fixed point equation doesn't involve any self-reference (because we can't do self-reference in provability logic as we have no way to talk about Godel numbers). Self-reference is less essential than you might think - we can solve these equations without using it.

At first the fixed point theorem seems like a positive thing: we can crank out as many of these theorems as we like. But there's also a flip side: it says that once we've proved Lob's theorem, there aren't any more techniques we need in order to construct fixed points. So there aren't any funky new variations on what Godel did, waiting to be discovered, that will give us a bunch of new fixed points. We have them all already.

In order to find these fixed points I'm going to use the runTableau function we wrote last time. But that's for the next article in the series. In this post I want to put it to a simpler use:

Automatically Drawing Tableau Diagrams
First we'll need a very simple library of functions for drawing ASCII art diagrams. We'll represent a drawing simply as a list of ASCII strings, one for each row. Think of diagrams as forming a box. The width of a box is the length of its longest row:

> width box = foldr max 0 (map length box)

Sometimes we'll want to pad the length of the rows so that they all have the same length:

> pad len b a = a ++ replicate (len - length a) b

The aside function allows us to 'typeset' one box next to another:

> aside a b = let w = width a
>                 h = max (length a) (length b)
>                 a' = pad h [] a
>                 b' = pad h [] b
>                 a'' = map (pad w ' ') a'
>             in zipWith (++) a'' b'

We can draw nice frames around our boxes:

> frame a = let w = width a
>               h = length a
>               strut = replicate h "|"
>               rule = "+" ++ replicate w '-' ++ "+"
>           in [rule] ++ strut `aside` a `aside` strut ++ [rule]

And the rule for disjunctions requires us to draw side by side boxes with a vertical line between them:

> alt a b = let h = max (length a) (length b)
>               strut = replicate h " | "
>           in a `aside` strut `aside` b

Generating diagrams is now a matter of generating a diagram for each of the 'hooks' in our algorithm:

> diagramRules = TableauRules {

We'll modify the rules so that instead of simply returning a Bool to indicate closure we return a pair. The first element is a list of strings representing the rows of the ASCII art, but the second element plays the same role as before:

>     closes             = snd,

These are the functions that handle discovery of a contradiction. They draw a "X" in a circle to indicate this:

>     foundF             = \_ -> (["(X)"], True),
>     foundContradiction = \_ -> (["(X)"], True),

>     open               = \_ -> ([], False),

We deal with conjunctions simply by listing the two subpropositions that went into them:

>     conjRule           = \a b -> first ([show a, show b] ++),

For disjunctions we draw the diagrams:

>     disjRule           = \p a b (tl, lb) (tr, rb) ->
>                             (([show a] ++ tl) `alt` ([show b] ++ tr), lb && rb),

>     doubleNegation     = \q -> first ([show q] ++),
>     combineWorlds      = \(t0, b0) (t1, b1) -> (t0 ++ t1, b0 || b1),
>     processWorld       = \_ t -> t,

When we create a box for a new world we 'import' some propositions from the parent world. We mark these with an asterisk:

>     tableau            = \br (t, b) -> (frame (map (("* " ++) . show) br ++ map ("  "++) t), b)
> }

> diagram p = do
>     let (t, b) = runTableau diagramRules [Neg p]
>     mapM_ putStrLn t
>     print b

Now we can prove that we really do have a fixed point of f1:

> diagram1 = diagram $ let p = Box (Box F) --> Box F in p <-> f1 p

You may want to pick a really small font and stretch your terminal wide before you run that:

See how every box ends in an (X), just as we want.

note
To run the code above, just append this article to the previous one to make a single literate Haskell program.

Quine Central

Quines for Everyone
This is an interruption to the sequence on provability logic (though it's not entirely unrelated).

The code below spits out a Haskell program that prints out a Perl program that prints out a Python program that prints out a Ruby program that prints out a C program that prints out a Java program that prints out the original program. Nothing new, an obvious generalisation of this. (Well, truth be told, it's a block of HTML that generates some Haskell code...)

But there's one big difference. To allow everyone else to join in the fun you can configure it yourself! Any non-empty list of languages will do, including length one. You may start hitting language line length limits if you make your list too long. Note that you can repeat languages so you can give someone some Haskell which decays into Perl after n iterations.

The code is geared towards generating tightly packed code but it's easily adapted to generate something slightly more readable. Apologies for the C warnings. Trivial to fix.

It's easily extended to support many more languages. C++, C#, obj-C, ocaml, Prolog, Lisp, Scheme and Go, say, should all be trivial apart from maybe a tiny bit of work with delimiting strings. (Eg. for Go you may need to tweak the import statement slightly so it doesnt't use double quotation marks.)

The code leaves many opportunities for refactoring but it's not like anyone is actually going to use this code for real production so I'm leaving it as is now.

I've only tested it under MacOS X. I don't know if there are carriage return/linefeed issues with other OSes. The shell script at the end is a regression test.

There was a little bit of theory involved which I learnt from Vicious Circles.

Here's a challenge for you: write a quine that takes as input the name of a language and outputs the same thing implemented in the input language. Much harder than what I just wrote.

> import Data.List

Here's the bit you can easily play with:

> langs = [Haskell, Perl, Python, Ruby, C, Java]

Implementation

> data Languages = Haskell | Ruby | Perl | C | Python | Java

> sequenceFromString Haskell s = "map toEnum[" ++ (intercalate "," $
>     map (\c -> show (fromEnum c)) s) ++ "]"
> sequenceFromString Perl s    = (intercalate "," $
>     map (\c -> "chr(" ++ show (fromEnum c) ++ ")") s)
> sequenceFromString Python s  = (intercalate "+" $
>     map (\c -> "chr(" ++ show (fromEnum c) ++ ")") s)
> sequenceFromString Ruby s    = (intercalate "+" $
>     map (\c -> show (fromEnum c) ++ ".chr") s)
> sequenceFromString C s       = concatMap
>     (\c -> "putchar(" ++ show (fromEnum c) ++ ");") s
> sequenceFromString Java s    = concatMap
>     (\c -> "o.write(" ++ show (fromEnum c) ++ ");") s

> paramList' Haskell = intercalate " " . map (\n -> "a" ++ show n)
> paramList' C       = intercalate "," . map (\n -> "char *a" ++ show n)
> paramList' Python  = intercalate "," . map (\n -> "a" ++ show n)
> paramList' Ruby    = intercalate "," . map (\n -> "a" ++ show n)
> paramList' Java    = intercalate "," . map (\n -> "String a" ++ show n)

> paramList Perl    _ = ""
> paramList lang n = paramList' lang [0..n-1]

> driver l args = defn l ++ intercalate (divider l) args ++ endDefn l

> divider C       = "\",\""
> divider Perl    = "','"
> divider Ruby    = "\",\""
> divider Python  = "\",\""
> divider Haskell = "\" \""
> divider Java    = "\",\""

> defn C       = "main(){q(\""
> defn Perl    = "&q('"
> defn Python  = "q(\""
> defn Ruby    = "q(\""
> defn Haskell = "main=q \""
> defn Java    = "public static void main(String[]args){q(\""

> endDefn C       = "\");}"
> endDefn Perl    = "')"
> endDefn Python  = "\")"
> endDefn Ruby    = "\")"
> endDefn Haskell = "\""
> endDefn Java    = "\");}}"

> arg Haskell n = "a" ++ show n
> arg Perl n    = "$_[" ++ show n ++ "]"
> arg C n       = "printf(a" ++ show n ++ ");"
> arg Python n  = "a" ++ show n
> arg Ruby n    = "a" ++ show n
> arg Java n    = "o.print(a" ++ show n ++ ");"

> argDivide Haskell l = "++" ++ sequenceFromString Haskell (divider l) ++ "++"
> argDivide Perl l    = ","    ++ sequenceFromString Perl (divider l) ++ ","
> argDivide C l       = sequenceFromString C (divider l)
> argDivide Python l  = "+" ++ sequenceFromString Python (divider l) ++ "+"
> argDivide Ruby l    = "+" ++ sequenceFromString Ruby (divider l) ++ "+"
> argDivide Java l    = sequenceFromString Java (divider l)

> argList lang1 lang2 n = intercalate (argDivide lang1 lang2) $
>     map (arg lang1) ([1..n-1] ++ [0])

> fromTo Haskell l n = "q " ++ paramList Haskell n ++ "=putStrLn$a0++" ++
>     sequenceFromString Haskell ("\n" ++ defn l) ++ "++" ++
>     argList Haskell l n ++ "++" ++ sequenceFromString Haskell (endDefn l)
> fromTo Perl    l n = "sub q {" ++ "print $_[0]," ++
>     sequenceFromString Perl ("\n" ++ defn l) ++ "," ++ argList Perl l n ++ "," ++
>     sequenceFromString Perl (endDefn l ++ "\n") ++ "}"
> fromTo Python  l n = "def q(" ++ paramList Python n ++
>     "): print a0+" ++ sequenceFromString Python ("\n" ++ defn l) ++
>     "+" ++ argList Python l n ++ "+" ++ sequenceFromString Python (endDefn l)
> fromTo Ruby    l n = "def q(" ++ paramList Ruby n ++
>     ") print a0+" ++ sequenceFromString Ruby ("\n" ++ defn l) ++
>     "+" ++ argList Ruby l n ++ "+" ++ sequenceFromString Ruby (endDefn l ++ "\n") ++ " end"
> fromTo C       l n = "q(" ++ paramList C n ++ "){" ++ "printf(a0);" ++
>     sequenceFromString C ("\n" ++ defn l) ++ argList C l n ++
>     sequenceFromString C (endDefn l ++ "\n") ++ "}"
> fromTo Java    l n = "public class quine{public static void q(" ++
>     paramList Java n ++ "){java.io.PrintStream o=System.out;o.print(a0);" ++
>     sequenceFromString Java ("\n" ++ defn l) ++ argList Java l n ++
>     sequenceFromString Java (endDefn l ++ "\n") ++ "}"

> main = do
>     let n = length langs
>     let langs' = cycle langs
>     putStrLn $ fromTo (head langs') (head (tail langs')) n
>     putStrLn $ driver (head langs') $ zipWith (\lang1 lang2 -> fromTo lang1 lang2 n)
>         (take n (tail langs')) (tail (tail langs'))

Regression Test
Assuming this article is stored in quineCentral.lhs.

runghc quineCentral.lhs>1.hs
cat 1.hs
echo "---------------------------------"
runghc 1.hs>2.pl
cat 2.pl
echo "---------------------------------"
perl 2.pl>3.py
cat 3.py
echo "---------------------------------"
python 3.py>4.ruby
cat 4.ruby
echo "---------------------------------"
ruby 4.ruby>5.c
cat 5.c
echo "---------------------------------"
gcc -o 5 5.c
./5>quine.java
cat quine.java
echo "---------------------------------"
javac quine.java
java quine>7.hs
cat 7.hs
diff 1.hs 7.hs

Generalising Gödel's Theorem with Multiple Worlds. Part II.

Introduction
Last time we looked at a method for testing whether propositions of provability logic were valid by looking at the consequences of propositions within nested collections of worlds. This lends itself naturally to an algorithm that we can implement in Haskell. The diagrams I used last time are a variant on what are known as tableaux. But tableaux can be used in a number of ways and so we need code that is suitably generalised. In the code that follows I've ensured that the core algorithm has an interface that is suitable for carrying out the four tasks I'll demand of it. This means that the code leans more towards practicality than mathematical elegance.

And I apologise in advance: this post is mostly a bunch of implementation details. We'll get back to some mathematics next time. Nonetheless, by the end of this post you'll have working code to test the validty of propositions of provability logic.

Central to tableau algorithms is pattern matching. I'd like the Haskell pattern matcher do much of the work, and to increase its flexibility I'll need this extension:

> {-# LANGUAGE ViewPatterns #-}

We'll need these libraries too:

> import Control.Applicative
> import Control.Arrow
> import Control.Monad
> import Data.Function
> import List
> import Maybe
> import Text.Show

Logical Propositions
Now we'll need a bunch of logical operators. The first three are constructors for a proposition type and the rest are sugar to make expressions look nicer:

> infixr 1 :->
> infixr 2 :\/
> infixr 3 :/\

> infixr 1 -->
> infixr 1 <--
> infixr 1 <->
> infixr 2 \/
> infixr 3 /\

> (\/)    = (:\/)
> (/\)    = (:/\)
> (-->)   = (:->)
> (<--)   = flip (:->)
> p <-> q = (p :-> q) :/\ (q :-> p)

Here's our basic proposition type:

> data Prop = Letter String
>           | Prop :\/ Prop
>           | Prop :/\ Prop
>           | Prop :-> Prop
>           | Box Prop
>           | Dia Prop
>           | F
>           | T
>           | Neg Prop deriving (Eq, Ord)

I want show to know about operator precedence for propositions:

> instance Show Prop where
>     showsPrec p (a :/\ b) = showParen (p>3) $ showsPrec 3 a . showString " /\\ " . showsPrec 3 b
>     showsPrec p (a :\/ b) = showParen (p>2) $ showsPrec 2 a . showString " \\/ " . showsPrec 2 b
>     showsPrec p (a :-> b) = showParen (p>1) $ showsPrec 1 a . showString " --> " . showsPrec 1 b
>     showsPrec p (Neg r)   = showParen (p>4) $ showString "Neg " . showsPrec 5 r
>     showsPrec p (Box r)   = showParen (p>4) $ showString "Box " . showsPrec 5 r
>     showsPrec p (Dia r)   = showParen (p>4) $ showString "Dia " . showsPrec 5 r
>     showsPrec p (Letter n)= showParen (p>5) $ showsPrec 6 n
>     showsPrec p T         = showString "T"
>     showsPrec p F         = showString "F"

Some simple rules for simplification of some logical expressions:

> simplify p = let simplify' (a :\/ F) = a
>                  simplify' (F :\/ b) = b
>                  simplify' (a :/\ T) = a
>                  simplify' (T :/\ b) = b
>                  simplify' (a :\/ T) = T
>                  simplify' (T :\/ b) = T
>                  simplify' (a :/\ F) = F
>                  simplify' (F :/\ b) = F
>                  simplify' (F :-> b) = T
>                  simplify' (T :-> b) = b
>                  simplify' (a :-> F) = Neg a
>                  simplify' (a :-> T) = T
>                  simplify' (Neg T) = F
>                  simplify' (Neg F) = T
>                  simplify' (Box T) = T
>                  simplify' (Dia F) = F
>                  simplify' z = z
>    in case p of
>        a :/\ b -> let a' = simplify a
>                       b' = simplify b
>                   in simplify' (a' :/\ b')
>        a :\/ b -> let a' = simplify a
>                       b' = simplify b
>                   in simplify' (a' :\/ b')
>        a :-> b -> simplify' (simplify a :-> simplify b)
>        Box a   -> simplify' (Box (simplify a))
>        Dia a   -> simplify' (Dia (simplify a))
>        Neg (Neg a) -> simplify a
>        a           -> a

Kinds of Proposition
I'm actually going to use more than one proposition type. So when we do case analysis I need to make my patterns more abstract so they can work with multiple types. I'm going to use the PropType type to represent the ways we're going to classify logical propositions. The types are:

1. Atomic: A single letter or its negation
2. Constant: Simply T or F or a negation thereof.
3. DoubleNegation.
4. Disjunction: used to represent things like a∧b or ¬(a∨b).
5. Conjunction: used to represent things like a∨b or ¬(a∧b).
6. Provability: These are statements about provability like those starting with ◻ or ¬◊.
7. Consistency: These are statements about consistency like those starting with ◊ or ¬◻.

As this type is a simple container we can make it a functor too.

> data PropType a = Atomic a
>                 | Constant a
>                 | DoubleNegation a
>                 | Disjunction a a
>                 | Conjunction a a
>                 | Provability a
>                 | Consistency a

> instance Functor PropType where
>     fmap f (Atomic a)         = Atomic (f a)
>     fmap f (Constant a)       = Constant (f a)
>     fmap f (DoubleNegation a) = DoubleNegation (f a)
>     fmap f (Provability a)    = Provability (f a)
>     fmap f (Consistency a)    = Consistency (f a)
>     fmap f (Conjunction a b)  = Conjunction (f a) (f b)
>     fmap f (Disjunction a b)  = Disjunction (f a) (f b)

I'll introduce a typeclass that will allow us to use the PropType view to query propositions:

> class PropTypeable a where
>     propType :: a -> PropType a
>     neg      :: a -> a
>     isF      :: a -> Bool
>     negative :: a -> Bool
>     positiveComponent :: a -> Prop

And now here we have the cases that I summarised in English above:

> instance PropTypeable Prop where
>     propType (a :\/ b)        = Disjunction a b
>     propType (a :/\ b)        = Conjunction a b
>     propType (Neg (a :\/ b))  = Conjunction (Neg a) (Neg b)
>     propType (Neg (a :/\ b))  = Disjunction (Neg a) (Neg b)
>     propType (a :-> b)        = Disjunction (Neg a) b
>     propType (Neg (a :-> b))  = Conjunction a (Neg b)
>     propType (Neg (Neg a))    = DoubleNegation a
>     propType (Box a)          = Provability a
>     propType (Neg (Box a))    = Consistency (Neg a)
>     propType (Dia a)          = Consistency a
>     propType (Neg (Dia a))    = Provability (Neg a)
>     propType (Letter a)       = Atomic (Letter a)
>     propType (Neg (Letter a)) = Atomic (Neg (Letter a))
>     propType T                = Constant T
>     propType F                = Constant F
>     propType (Neg F)          = Constant T
>     propType (Neg T)          = Constant F
>     neg                       = Neg
>     isF F                     = True
>     isF (Neg T)               = True
>     isF _                     = False
>     positiveComponent (Neg a) = a
>     positiveComponent a       = a
>     negative (Neg _)          = True
>     negative _                = False

It'll be a while before we need the full generality so it's going to seem like overkill for the moment!

And some pre-packaged letters for convenience:

> [a, b, c, d, p, q, r, s, t] = map (Letter . return) "abcdpqrst"

We're going to need some operations that act on lists.

placesWhere finds all of the elements of a list for which some predicate holds. Instead of just listing the elements that match, it lists the elements paired with the rest of the list after the matching element is removed. We can think of these pairs as elements and their surrounding context:

> placesWhere p []     = []
> placesWhere p (x:xs) = let r = map (second (x:)) $ placesWhere p xs
>                        in if p x then ((x, xs) : r) else r

This finds something in the intersection of two sets using a given 'equality' predicate for matching. As we may be using a predicate different from == we need to see both of the (possibly different) elements that satisfy the predicate.

> findIntersection eq a b = listToMaybe [(x, y) | x <- a, y <- b, x `eq` y]

Sometimes we'll meet propositions that we can start to reason about using ordinary propositional calculus. These will match the propositional predicate:

> propositional (propType -> DoubleNegation _) = True
> propositional (propType -> Conjunction _ _)  = True
> propositional (propType -> Disjunction _ _)  = True
> propositional _                              = False

On the other hand the provability predicate is used to identify propositions that need rules pertaining to provability and consistency:

> provability (propType -> Provability _) = True
> provability (propType -> Consistency _) = True
> provability _                           = False

The Algorithm
And now we're almost ready to implement the tableau rules. Because we'll be using tableaux in a number of different ways I need lots of hooks into the algorithm that can perform different operations. I've collected all of these hooks into a single type. The algorithm will take a proposition of type prop (which will be Prop for the first three cases) and produce something of type result.

> data TableauRules prop result = TableauRules {

If our result corresponds to a world that is self-contradictory it is said to close. Here's how we indicate a closed (and hence not really existing) world. In the simplest case we won't actually store any information about a world, just whether or not it closes. So closes will be the identity function:

>     closes :: result -> Bool,

Occasionally we'll find a world with something obviously equivalent to F in it. It closes. Here's what we want to return in that case. The argument is the offending proposition:

>     foundF :: prop -> result,

Sometimes we'll find a pair that obviosuly contradict, like a and ¬a:

>     foundContradiction :: (prop, prop) -> result,

And sometimes we'll find an open world (ie. a real non-closed one). This function gets handed the list of propositions that have been found to hold in it:

>     open :: [prop] -> result,

Here's what we do when we find a conjunction. I hope you remember that when we meet a conjunction we can delete it and replace it with the two subpropositions. conjRule is handed the subpropositions as well as the result from proceeding with the tableau rule for conjunctions:

>     conjRule :: prop -> prop -> result -> result,

Disjunctions work a little differently. When handling a∨b, say, we need to handle two subtableaux, one with a and one with b. The first argument to disjRule rules is the disjunction itself, the next two are the left and right subpropositions, and the last two arguments are the results of continuing the two subtableaux.

>     disjRule :: prop -> prop -> prop -> result -> result -> result,

With doubleNegation we get to see propositions that have undergone double negation elimination.

>     doubleNegation :: prop -> result -> result,

When we use Dia to open new worlds we need to ensure that each of these subworlds is valid. Each subworld is processed with processWorld . For example, when we're drawing tableau diagrams we can use this book to draw a frame around the subtableaux. We then fold together these subworlds using combineWorlds:

>     processWorld  :: prop -> result -> result,
>     combineWorlds :: result -> result -> result,

Lastly we have our driver function that kicks off the whole tableau algorithm on a list of propositions:

>     tableau :: [prop] -> result -> result
> }

We can now implement a tableau computation algorithm that supports these hooks. We'll start with the check for whether or not we have an immediate contradiction:

> simpleClosure rules ps = case find isF ps of
>    Just a  -> foundF rules a

Split the propositions into those that are negated and those that are not. We're looking for propositions in one part that directly contradict propositions in the other:

>    Nothing ->
>      let (neg, pos) = partition negative ps
>          maybePair = findIntersection ((==) `on` positiveComponent) neg pos
>      in case maybePair of
>          Just pair -> foundContradiction rules pair
>          Nothing   -> open rules ps

Double negation elimination is straightforward to apply. Delete the original proposition and replace it with the version without double negation:

> applyDNeg rules p a props = doubleNegation rules a $
>   applyPropositional rules (a : delete p props)

The rule for handling conjunctions. We delete the conjunction from the current list of propositions and replace it with the two subpropositions:

> applyConj rules p a b props = conjRule rules a b $
>   applyPropositional rules (a : b : delete p props)

Disjunctions require running two separate subtableaux:

> applyDisj rules p a b props =
>    let props' = delete p props
>        left   = applyPropositional rules (a : props')
>        right  = applyPropositional rules (b : props')
>    in disjRule rules p a b left right

Here we tie together the rules for propositional calculus. We use a bit of case analysis to decide which rule to apply, and if no rule applies we try the provability rules instead:

> applyPropositional rules props =
>     let t = simpleClosure rules props in if closes rules t
>         then t
>         else case find propositional props of
>             Nothing -> applyProvability t rules props
>             Just p  -> case p of
>                 (propType -> DoubleNegation q) -> applyDNeg rules p q props
>                 (propType -> Conjunction a b)  -> applyConj rules p a b props
>                 (propType -> Disjunction a b)  -> applyDisj rules p a b props

When we've exhausted all possible rules from propositional calculus we scan for propositions like Dia p or Neg (Box p). These may imply the existence of subworlds. We then try to instantiate these subworlds, seeded according to the rules I gave in the previous article:

> applyProvability t rules props =
>     let impliedWorlds = placesWhere consistency props

>         consistency (propType -> Consistency _) = True
>         consistency _ = False

>         testWorld (p@(propType -> Consistency q), props) =

In the following line, neg p corresponds to the application of Löb's Theorem. provabilities is the list of propositions inherited by a subworld from statements about provability in the parent:

>              let tableau = runTableau rules (q : neg p : provabilities)
>                  provabilities = do 
>                      p@(propType -> Provability q) <- props
>                      [p, q]
>              in processWorld rules p tableau

>     in foldr (combineWorlds rules) t (map testWorld impliedWorlds)

And finally, here's where we kick our algorithm off:

> runTableau rules props = tableau rules props $ applyPropositional rules props

Testing Validity
We can now use a set of simple rules to test the validity of propositions. As mentioned above, the return data is a Bool used to indicate whether a subworld was invalid. By and large, these rules do the trivial thing. Note how the rule for disjunction requires both of the alternatives to be invalid in order to completely invalidate, so we use (&&). But when considering subworlds, just one bad subworlds is enough to invalidate a world:

> validRules = TableauRules {
>     closes = id,
>     open   = \_ -> False,

>     foundF             = \_ -> True,
>     foundContradiction = \_ -> True,

>     conjRule       = \_ _ t -> t,
>     disjRule       = \_ _ _ -> (&&),
>     doubleNegation = \_ t -> t,

>     combineWorlds = (||),
>     processWorld  = \_ t -> t,

>     tableau = \_ t -> t
> }

We can now write a simple validty test. We negate the proposition we're interested in and test whether the implied world closes:

> valid p = runTableau validRules [neg p]

Here's a small regression test to ensure everything works. It's just a bunch of examples that I worked out by hand or lifted from Boolos's book:

> valids = [
>         T,
>         a :-> a,
>         Box a :-> Box a,
>         Box a :-> Box (Box a),
>         Box (Box a :-> a) :-> Box a,
>         Box F <-> Box (Dia T),
>         let x = p :/\ q :-> r :-> a in Box (Box x :-> x) :-> Box x,
>         F :-> Dia p,
>         Box (Dia p) :-> Box (Box F :-> F),
>         (Box F \/ q /\ Dia (Box F /\ Neg q)) <->
>           (Dia (Box F \/ q /\ Dia (Box F /\ Neg q))
>           --> q /\ Neg (Box (Box F \/ q /\ Dia (Box F /\ Neg q)
>           --> q)))
>     ]

> invalids = [
>         F,
>         a :-> Box a,
>         Box a :-> a,
>         Box (Box a :-> a) :-> a,
>         Dia T,
>         Box (Dia T),
>         Neg (Box F),
>         (Box F \/ p /\ Dia (Box F /\ Neg q)) <->
>           (Dia (Box F \/ q /\ Dia (Box F /\ Neg q))
>           --> q /\ Neg (Box (Box F \/ q /\ Dia (Box F /\ Neg q)
>           --> q)))
>     ]

If everything is working, regress1 should give the result True.

> regress1 = do
>     print $ (and $ map valid valids) &&
>             (and $ map (not . valid) invalids)

That's enough implementation. In the next installment we'll start putting this code to work.

References
This code is an implementation of the algorithm described in Chapter 10 of The Logic of Provability. See that book if you want (1) a proof that the above algorithm always terminates and (2) that it really does correctly decide the validity of propositions of provability logic.

Generalising Gödel's Theorem with Multiple Worlds. Part I.

Introduction
In his first incompleteness theorem, Gödel showed us that we can construct a sentence that denies its own provability. In his second incompleteness theorem he showed that an example of such a sentence is the one that asserts the inconsistency of arithmetic. If arithmetic is consistent then it can't prove its own consistency. On the other hand, if arithmetic is inconsistent then we can prove anything, and hence we can prove its consistency.

Can we generalise what Gödel did? For example, can we construct sentences that we can prove assert their own provability? What about sentences that deny that their provability is provable? Or what about sentences that assert that if they're provable then it's not provable that it's inconsistent that they imply that they're inconsistent with the rest of arithmetic?

Not only can we do these things, we can also write a computer program that generates such theorems for us. We can do so by working with the idea that a consistent set of axioms describes a world, and that a set of axioms able to talk about sets of axioms describes a world within a world, and a set of axioms that...well you can guess how it goes. A bit like Inception really.

Gödel's Theorem
Briefly, Gödel's first incompleteness theorem goes like this: we work with PA, the logical system built from Peano's Axioms. Within PA we can state and prove theorems about arithmetic like the fact that 1+2=3 or that for any prime there is always another greater prime. But even though PA is about numbers, we can talk about other things if we can encode them as numbers. In particular, Gödel came up with a scheme to encode propositions of PA as numbers. We use [P] to represent the number for the proposition P in Gödel's scheme. A proof is essentially just a list of propositions where each one is connected to some earlier ones by simple mechanical rules. These rules can be turned into arithmetical statements about the Gödel numbers for the propositions. So in the language of PA it is possible to assert that a particular number is the Gödel number of a provable proposition. Let Prov(n) denote the proposition of PA that says that n is the Gödel number of a provable proposition. Prov(n) is a proposition of PA and so is ¬Prov(n). Imagine we could find a proposition G with the property that G↔¬Prov([G]). It would assert its own unprovability. But it also appears to involve stuffing a representation of the Gödel number of G within G as well as all the rules for determining provability. Amazingly Gödel figured out how to do this (using tricks not unlike those used for quining). If G were false, G would be provable. So if we trust PA as a means of reasoning about proofs in PA, then G is true, though it can't be proved using the axioms of PA.

Provability Logic
We're going to be interested specifically in provability, so we don't need all of the power of PA. So we're going to work with a simplified domain specific logic called GL (for Gödel-Löb), otherwise known as Provability Logic. The idea is that sentences of GL will be shorthand for classes of statement in PA. GL will contain propositional calculus. So here's an example statement in GL: p∧q. The unknowns p and q represent propositions of PA and statements of GL are considered valid if they're provable whatever statements of PA they represent. For example, we could assign p="there are at least 10 primes" and q="7>2", in which case p∧q holds. But we could assign p="there are just 10 primes" in which case it's false. So p∧q isn't a valid proposition of GL as it doesn't hold for all p and q. On the other hand p→p∨q is valid because no matter what crazy propositions we assign to p and q, p→p∨q is true. (Of course, when I say "there are at least 10 primes", I mean a long and complicated sentence of PA that amounts to the same thing.)

But there's more to GL than propositional calculus. It also has the one-argument predicate ◻ which asserts that its argument is provable. More precisely, ◻p says that whatever proposition of PA p represents, let's say it's P, we have Prov([P]). It's just shorthand. Here's another example: ◻(q→◻p) says that whatever P we assign to p, and whatever Q we assign to q, Prov([Q →Prov([P])]). Or in English it says "for any propositions P and Q, it is provable that Q implies that P is provable".

We'll use ⊤ and ⊥ from propositional calculus as the always true and always false propositions in the usual way. ◻⊥ is the assertion that we can prove ⊥. In other words it's the assertion that PA is inconsistent. So now we can state Gödel's second incompleteness theorem: ¬◻⊥→¬◻¬◻⊥. If PA is consistent, we can't prove it.

We can also introduce the symbol ◊. This is just shorthand for ¬◻¬. ◊p says that it's not provable that p is false. In other words, it says that p is consistent with the rest of PA. So ◻ is provability and ◊ is consistency.

A set of assignments of propositions of PA to a bunch of letters in GL can be thought of as a world. For example, we can imagine a world in which p="2>1" and q="2<1". In this world, we have p and ¬q. The GL proposition ◻p says that we have a proof of p so it must be true in all worlds. Conversely, ◊p says that it's not true there is no world where p holds. In other words, it asserts the existence of a world where p holds. So worlds can talk about other worlds. If we have ◊◊p in some world, then it's asserting that there's another world in which ◊p holds. In other words, ◊◊p asserts there is a world in which it is asserted that there is another world where p holds. We can draw pictures to represent this. If a world has propositions that talk about another world then we draw the talked about world as a kind of subworld. Here's how we can picture ◊◊p:

Worlds are assignments of propositions of PA to letters of GL. But most of the time we won't particularly care about specific propositions themselves like "2>1". We'll be more interested in what truth assignments we can make to the propositions represented by the letters. So we can think of a world as a place where the letters of GL have been assigned truth values, true or false. And we can think of each world as containing subworlds consisting of propositions that can be proved or disproved by their parent worlds.

I'm going to spend most of my time looking at how we can explore and unfold all of the implications contained in a world.

The Rules of the Game
We can give some rules. If a world contains p∧q like this:

then it must contain both p and q as well:

I crossed out the p∧q as it became redundant. It's important to note that when I wrote p∧q we could have any propositions of GL standing in for p and q. So a world containing (p∨q)∧r must also contain p∨q and r. Other rules may apply too. If we had p∧q and r∧s we can unfold both to get p, q, r and s. This rule also kicks in if it applies to negated propositions that become conjunctions if we "push down" the negation using de Morgan's laws. For example if a world contains ¬(p→q) then it also contains p and ¬q.

If a world contains p∨q then we don't know exactly what it looks like.

There are two possibilities. We can draw both together like this:

The line in the middle means that the world either looks like what's on the left or it looks like what's on the right. Once we've indicated there are two alternatives then the original proposition became redundant again. When we have a vertical line like this then everything above the vertical line applies in both the left and right possibilities. This saves us having to split our world into two and copy all of our formulae to both sides. Like in the case of conjunctions this rule also kicks in for other kinds of disjunctions. So a world containing p→q splits into two separate worlds headed by ¬p and q.

If a world contains ⊥, or a contradiction, in means that it wasn't really a valid world after all. If we meet a world like this, we know that our starting point must have been been self-contradictory. If a world implies the existence of a subworld that isn't valid, then it must itself be invalid:

On the other hand, if we've split a world into two possibilities because we found p∨q then even if one branch is invalid the world might still be valid if the other branch is valid.

If we have a world containing ◻p:

then we know that p must be in every subworld. Actually, we know it must also hold in any subsubworld too, all the way down. This is because if we can prove something we can then use that proof directly to form a constructive proof that we can prove it. So that means ◻p holds in every subworld too:

We treat negation as above. So if we see ¬◻p we treat it like ◊¬p.

As we have also said, if we have a world with a ◊p in it then it must contain a subworld with p in it. Note that ◊ implies the existence of subworlds, but ◻ doesn't. It just tells us what must be in them if they exist.

But there's one last rule we'll need. It's Lob's theorem. This is the big theorem on which everything else I say depends. It states that ◻(◻p →p) →◻p. If we can prove that proving something implies its truth then we can prove it. I could sketch a proof here, but I highly recommend the cartoon proof here. In a way, Löb's theorem is a bit sad. The raison d'etre of mathematics is that we can use proofs to be sure of things. In other words, we take for granted that ◻p→p. But if we could prove this then we could prove ◻p, even if p were false! ◻p→p is not valid!

(Philosophical digression: Mathematicians assume ◻p→p. They don't assume ◻p→p because there's a proof. They assume it because experience has shown it to work. So I claim that ◻p→p is an empirical fact that we learn by scientific induction. That's controversial because I'm basically saying that much of mathematics is empirical - at least it is if you talk about truth rather than proof. If you disagree, don't worry about it. The rest of what I say here is independent of this digression.)

Anyway, we can flip Löb's theorem around to get ◊p →◊(p ∧¬◊p). So if we have a world in which ◊p holds:

then we have a subworld in which both p and ¬◊p hold:

And that's all we'll need apart from the obvious rule that we can remove double negation. Suppose we have a proposition of GL. We can use the rules above to extract as many implications as we can. Eventually there will come a point where there is nothing more we can do. If we can do this without hitting a contradiction then we've found a bunch of possible worlds for the proposition. If we can't find a valid world, however, the the original proposition must have been false. You might ask whether or not there any other rules we need in addition to the ones above. Maybe there are other theorems like Löb's theorem that we need to use. Amazingly it can be shown that no other rules are needed. This is a sure-fire terminating algorithm for determining whether or not a proposition of GL is valid! This is a powerful tool. We can now start constructing wild propositions like those I started with and find out whether they are valid. (Note that I've not in any way proved this procedure always works. You'll need to look at Boolos's book to see why.)

Some Proofs
Let's work through some examples. First I'll prove something from propositional calculus: (p∧q)∨(p∧r)→p. We start by drawing a world with the negation:

Now ¬(a→b) is the same as a∧¬b. So we get:

Now the only way to proceed is to consider the disjunction and consider two alternatives:

On both sides of the vertical line we find p. But that contradicts the ¬p we discovered earlier. So there is no way the negation of our original proposition can hold in any world. And therefore the original proposition must be valid.

How about showing Gödel's second incompleteness theorem this way. ◊⊤ says that ⊤ is consistent with the rest of arithmetic. Ie. it expresses the consistency of arithmetic. The theorem is then ◊⊤→¬◻◊⊤, ie. if PA is consistent, then PA can't prove it is consistent. We'll start with this:

We can deal with the negated implication like before. We can also deal straightforwardly with the double negation:

There's now no way to proceed except to use the ◊⊤ to open up a new world. Remember that when we use this we must use Löb's theorem as well as inheriting p and ◻p from the parent world's ◻p:

And we get a contradiction because we have both ◊⊤ and ¬◊⊤. So Gödel's second incompleteness theorem is indeed valid. Note that this isn't a demonstration from scratch. We've shown its validity from Löb's theorem. So this isn't really a useful way to show it's valid. But it *is* a useful way to show the validity of generalisations of Gödel's theorem. Unfortunately, I have to stop for now.

In the coming posts I'll implement all of the above as a computer program. If there's any ambiguity in what I've said I hope the source code to the program will resolve those ambiguities. What's more, our program won't just test the validity of a proposition but it will draw out a nice picture of our world with its subworlds. So you'll be able to trace through every step to make sure you understand!

But that's not all. We'll also see how to write a program to illustrate a bunch more theorems like Craig's Interpolation Lemma and Beth's Definability Theorem and then we'll finish with a program that is designed to construct self-referential propositions. In particular, given any self-referential description like "p is a proposition that is equivalent to the proposition that denies that it's provable that p's provability is consistent with p itself" it will solve to find p, even though GL doesn't allow us to directly construct self-referential propositions.

So let me recap: a world is an assignment of consistent truth values to the letters (and consequently propositions) of GL. Some of these propositions imply the existence of other subworlds with different truth values for these propositions. We draw these worlds as subworlds of the original world. For a world to be valid it mustn't contain any contradictions (unless the world contains a bunch of alternatives in which case just one alternative needs to be valid.) A proposition of GL is valid if unfolding its negation doesn't result in an invalid world.

Exercises
Which are the following are valid? If you report back any difficulties you have I can incorporate any needed revisions into the description above.
1. p→◻p
2. ◻p→◻◻p
3. ◊p→◊◻◻p
4. ◊(p→◻q) →◊(◻(◊(p∧q)))
5. ◻(p∧q∧r∧s) →◻p∧◻q∧◻r∧◻s

References
Most of this stuff is based on Chapter 10 of The Logic of Provability by Boolos. The idea of using multiple worlds to prove theorems is due to Kripke. I believe the procedure of unfolding the implications of a proposition in the tree-like way I describe above is due to Smullyan.

Update
Solutions to problems:
1. Not valid. Diagram:

2. Valid
3. Valid.
4. Valid.
5. Valid. Diagram:

Beating the odds with entangled qubits

Quantum mechanics allows the possibility of "spooky action at a distance", correlations between widely separated but simultaneous random events that can't be explained by probability theory. These events look like they secretly communicate with each other, but we also know that quantum mechanics prevents us sending messages faster than the speed of light. Nonetheless, even though we can't exploit non-locality to send messages faster than the speed of light, two cooperating parties can exploit non-locality to perform tasks better than would be possible without non-locality. The CHSH game is one such example.

My goal here is to write code to emulate the CHSH game. It will require reusing the probability and quantum mechanics monads I've used here many times before. So I won't be explaining how these work.

> import Data.Map (toList, fromListWith)
> import Complex
> infixl 7 .*

The CHSH game is cooperative in the sense that the two players, A and B, are attempting to work together to win. The two players are widely separated. Between the two players is the game show host. The host randomly generates a pair of bits, s and t. s is sent to A and t is sent to B. Neither A nor B gets to see the message sent to their partner. A and B must now simultaneously make their moves, stating a choice of bit.

Call A's move a and B's move b. A and B win if they can arrange that a XOR b equals s AND t. So, for example, if A receives a true bit, and thinks B has also received a true bit, then A wants to make a move that differs from B's, otherwise A wants to make the same move as B. A and B are allowed to plan as much as they like beforehand but it should be pretty clear that they can't possibly guarantee a win.

We can formally write the victory condition as:

> victory a b s t = (a `xor` b) == (s && t)

Now there's a pretty good strategy A and B can adopt: three quarters of the time, s AND t will be false. In that case, A and B want their answers to match. So they could simply choose false, regardless of what message the game show host sends them. We can give their strategies as a function of the host's message:

> astrategy s = False
> bstrategy t = False

Now we can simulate our game:

> game = do

The host picks a random bit to send to each player:

>    s <- 0.5 .* return False + 0.5 .* return True
>    t <- 0.5 .* return False + 0.5 .* return True

The players now respond to each of their messages:

>    let a = astrategy s
>    let b = bstrategy t

Now we can collect the replies and score the result:

>    let score = victory a b s t
>    return score

> play1 = collect game

Running play1 gives the expected result that A and B have a 3/4 chance of winning. It's not hard to prove classically that they can do no better than this.

But in a quantum universe it is possible to do better! We now allow A and B to adopt strategies that involve making measurements of a quantum system. To describe their strategies we need to use the quantum monad. Here is the previous strategy rewritten for this monad. In this case, the argument b is the state of the quantum system they observe. The first element of the pair is the move in the game, the second element is the state the physical system is left in by the player:

> aqstrategy b s = qreturn (False, b)
> bqstrategy b t = qreturn (False, b)

Now we can rewrite game to support quantum processes:

> game' aqstrategy bqstrategy bits = do
>    s <- 0.5 .* preturn False + 0.5 .* preturn True
>    t <- 0.5 .* preturn False + 0.5 .* preturn True
>    (score, _, _) <- collect $ observe $ do
>        (abit, bbit) <- bits
>        (a, abit') <- aqstrategy abit s
>        (b, bbit') <- bqstrategy bbit t
>        let score = victory a b s t

Note that we have to return abit' and bbit' because quantum processes are reversible and can't erase information about a state. Also note that abit' and bbit' can be widely separated in space.

>        qreturn (score, abit', bbit')

Probabilistic processes can erase whatever they like:

>    preturn score

The quantum version of a coin toss to generate a random bit, ie. an equal superposition of False and True:

> coin' = (1/sqrt 2) .* qreturn False + (1/sqrt 2) .* qreturn True

> bits = do
>   a <- coin'
>   b <- coin'
>   return (a, b)

So now we play exactly as in play1 except that we give each player an independent qubit, which they ignore:

> play2 = collect $ game' aqstrategy bqstrategy bits

Unsurprisingly, the probability of winning is just the same as before.

But now we can try something impossible in the classical case. We give each player half of a perfectly correlated pair of qubits. The players are now entangled and can exploit non-locality. Of course if their strategies ignore the qubits we get the same result as before:

> bell = (1/sqrt 2) .* qreturn (False, False) + (1/sqrt 2) .* qreturn (True, True)
> play3 = collect $ game' aqstrategy bqstrategy bell

Now comes the surprising bit. The players can each look at the (classical) bit given to them by the game host. Depending what it is they rotate the qubit's state through some angle in state space. (In the case of the qubit being an electron state, this is an actual physical rotation of the electron.) In these strategies, the choice of move is the same as the state the qubit is left in after observation, hence the qreturn (b', b') bit.

> aqstrategy' b False = rotate (0) b >>= \b' -> qreturn (b', b')
> aqstrategy' b True  = rotate (pi/2) b >>= \b' -> qreturn (b', b')

> bqstrategy' b False = rotate (pi/4) b >>= \b' -> qreturn (b', b')
> bqstrategy' b True  = rotate (-pi/4) b >>= \b' -> qreturn (b', b')

And now when we play, the probability of winning is greater than 3/4.

> play4 = collect $ game' aqstrategy' bqstrategy' bell

All of the 'communication' took place before the game started. A and B didn't communicate s and t to each other. And yet they can beat the classical odds.

So in conclusion: Quantum mechanics gives opportunities for collusion that are impossible classically. Sadly we don't yet know how to maintain the state of separated entangled qubits for extended periods of time. But I remember seeing recently that people are managing to maintain qubit states for nanoseconds.

By the way, there are games where it is possible to achieve a 100% success rate with the help of quantum states. These give examples of what is known as quantum pseudo-telepathy. I presume the "pseudo" is because despite the 100% success rate, it still doesn't give a way to send messages instantly.

A last thought from me: one reason why humans send messages is to allow them to coordinate strategies. But quantum game theory shows that we can coordinate strategies without sending messages. In other words, even though non-locality doesn't give us faster-than-light communication, it does allow us to do things that were previously thought to require FTL. I think this may have some profound consequences.

And an example from a different domain: in biochemistry one could imagine remote parts of ligands coordinating the way they bind to receptors, something that would be completely missed by the kind of quasi-classical simulation I've seen biochemists use.

---

My standard quantum mechanics code:

> data W b a = W { runW :: [(a, b)] } deriving (Eq, Show, Ord)

> mapW f (W l) = W $ map (\(a, b) -> (a, f b)) l

> instance Functor (W b) where
>  fmap f (W a) = W $ map (\(a, p) -> (f a, p)) a

> instance Num b => Monad (W b) where
>  return x = W [(x, 1)]
>  l >>= f = W $ concatMap (\(W d, p) -> map (\(x, q)->(x, p*q)) d) (runW $ fmap f l)

> a .* b = mapW (a*) b

> instance (Eq a, Show a, Num b) => Num (W b a) where
>   W a + W b = W $ (a ++ b)
>   a - b = a + (-1) .* b
>   _ * _ = error "Num is annoying"
>   abs _ = error "Num is annoying"
>   signum _ = error "Num is annoying"
>   fromInteger a = if a==0 then W [] else error "fromInteger can only take zero argument"

> collect :: (Ord a, Num b) => W b a -> W b a
> collect = W . filter ((/= 0) . snd) . toList . fromListWith (+) . runW

> type P a = W Double a
> type Q a = W (Complex Double) a

> a `xor` b = a/=b

> rotate :: Double -> Bool -> Q Bool
> rotate theta True = let theta' = theta :+ 0
>   in cos (theta'/2) .* return True - sin (theta'/2) .* return False
> rotate theta False = let theta' = theta :+ 0
>   in cos (theta'/2) .* return False + sin (theta'/2) .* return True

> observe :: Ord a => Q a -> P a
> observe = W . map (\(a, w) -> (a, magnitude (w*w))) . runW . collect

Some help for the compiler (and maybe humans too):

> preturn = return :: a -> P a
> qreturn = return :: a -> Q a

Statistical Fingertrees

I have no time to post a proper article. But I have time to post a mini-article with details in the links.

> {-# LANGUAGE MultiParamTypeClasses #-}
> import Data.FingerTree
> import Data.Monoid

I'm a bad person. I often use this method to compute the mean and variance of some samples. It's not robust. It performs badly if the variance is small compared to the size of the samples. I shouldn't use it. As penance, this is a quick article on dynamically computing variances robustly and efficiently for datasets that are frequently manipulated.

For convenience I'll talk about the unscaled variance: the sum, not the average of the square deviation from the mean. You can easily compute the variance from this if you know the dataset size.

If we know the size, mean and unscaled variance of two sets (more properly, multisets) we can find the mean and unscaled variance of the their union using the formula here. This method is much more robust that the naive algorithm.

The rule for combining two datasets gives us a monoid:

> data Stats = Stats { n :: Float, mean :: Float, unscaledVariance :: Float }
>              deriving Show

> instance Monoid Stats where
>     mempty = Stats 0 0 undefined
>     Stats n m v `mappend` Stats 0 _ _ = Stats n m v
>     Stats 0 _ _ `mappend` Stats n m v = Stats n m v
>     Stats n m v `mappend` Stats n' m' v' = 
>       let delta = m' - m
>       in Stats (n + n') ((n*m+n'*m')/(n+n'))
>                (v + v' + delta*delta*n*n'/(n+n'))

Given a single sample, we can compute its stats:

> instance Measured Stats Float where
>     measure x = Stats 1 x 0

Now we need just one more line of code:

> type StatsTree = FingerTree Stats Float

We now have a data structure that allows us to freely split, join, delete elements from and add elements to sequences of samples, all the while robustly keeping track of their mean and unscaled variance (and hence their variance).

For example:

> example = fromList [1..10] :: StatsTree
> test = let (_, b) = split ((>=4) . n) example
>            (c, _) = split ((>3)  . n) b
>        in measure c

computes the stats for the 3 elements starting at the 4th element of example.

An example application might be maintaining rolling averages and variances for a sliding window.

This was inspired an article by John D Cook somewhere around here.

On Removing Singularities from Rational Functions

Introduction
Suppose we have the function

> f x = 1/(x+x^2) - 1/(x+2*x^2)

Some basic algebraic manipulation shows that in the limit as x→0, f(x)→1. But we can't simply compute f 0 because this computation involves division by zero at intermediate stages. How can we automate the process of computing the limit without implementing symbolic algebra?

I've already described one way to remove singularities from a function. But that approach is very limited in its applicability.

This article is about a variation on the approach to formal power series that nicely showcases some advantages of lazy lists. It will allow us to form Laurent series of functions so we can keep track of the singularities.

The usual Haskell approach to power series allows you to examine the coefficients of any term in the power series of the functions you can form. These series can't be used, however, to evaluate the function. Doing so requires summing an infinite series, but we can't do so reliably because no matter how many terms in a power series we add, we can never be sure that there aren't more large terms further downstream that we haven't reached yet. And if we want to perform computations completely over the rationals, say, we don't want to be dealing with infinite sums.

I'd like to look at a way of working with power series that allows us to perform exact computations making it possible to answer questions like "what is the sum of all the terms in this power series starting with the x^n term?" By extending to Laurent series, and implementing the ability to selectively sum over just the terms with non-negative powers, we can compute functions like f above at 0 and simply skip over the troublesome http://en.wikipedia.org/wiki/Pole_%28complex_analysis%29.

Power Series
When I previously discussed power series I used code that worked with the coefficients of the power series. This time we want to work with values of the function so it makes sense to store, not the coefficients a_i but the terms themselves, a_ixⁱ. So instead of a list of coefficients, Num a => [a] we need a representation that looks a little like:

> data Power a = Power (a -> [a])

where we pass x in as an argument to the function contained in a Power. But we also want to allow Laurent series so we need to also store an offset to say which (possibly negative) term our series starts with:

> data Laurent a = Laurent (a -> (Int, [a]))

But this fails us for at least two reasons:

1. We have the individual terms, but to evaluate the function requires summing all of the terms in an infinite list.
2. If we have a Laurent series, then we need to store values of a_ixⁱ for x=0 and i<0. We'll end up with division by zero errors.

Partial Sum Series
So here's what we'll do instead. Suppose our power series is Σ_i=n^∞a_ixⁱ. We'll store the terms s_j=Σ_i=j^∞a_ix^i-j. Our type will look like:

> data Partial a = Partial (a -> (Int, [a]))

> instance Eq (Partial a)
> instance Show (Partial a)

It's straightforward to add two functions in this form. We just add them term by term after first aligning them so that the xⁱ term in one is lined up with the xⁱ term in the other:

> instance Num a => Num (Partial a) where

>  Partial f + Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         pad 0 _ ys = ys
>         pad n x ys = let z:zs = pad (n-1) x ys
>                      in x*z : z : zs
>         l = min m n
>     in (l, zipWith (+) (pad (m-l) x xs) (pad (n-l) x ys))

Notice the slight subtlety in the alignment routine pad. By the definition above, the jth term has a factor of x^j built into it. So we need to multiply by x each time we pad our list on the left.

Now we need to multiply series. We know from ordinary power series that we need some sort of convolution. But it looks like for this case we have an extra complication. We appear to need to difference our representation to get back the original terms, convolve, and then resum. Amazingly, we don't need to do this at all. We can convolve 'in place' so to speak.

Here's what an ordinary convolution looks like when we want to multiply the sequence of terms (a_i) by (b_i):

In this example, the blue diagonal corresponds to the terms that are summed to get the 4th term in the result.

However, we wish to work with partial sums s_j=Σ_i=j^∞a_ix^i-j and t_j=Σ_i=j^∞t_ix^i-j, constructing the partial sums of the convolution of a and b from s and t. The partial sums of the convolution can be derived from the partial sums by tweaking the convolution so it looks like this:

The blue terms work just like before and need to be summed. But we also need to subtract off the red terms, weighted by a factor of x. That's it! (I'll leave that as an exercise to prove. The inclusion-exclusion principle helps.)

The neat thing is that the red terms for each sum are a subset of the blue terms needed for the next element. We don't need to perform two separate sums. We can share much of the computation between the red and blue terms. All we need to do is write an ordinary convolution routine that additionally returns not just the blue terms, but a pair containing the blue sum and the red sum.

>  Partial f * Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         (outer, inner) = convolve xs ys
>         f' a b = a-x*b -- (the subtraction I mentioned above)
>     in (m+n, zipWith f' outer inner)
>  fromInteger n = let n' = fromInteger n in Partial $ \_ -> (0, n' : repeat 0)

>  negate (Partial f) = Partial $ \x -> let (m, xs) = f x
>                                 in (m, map negate xs)
>  signum = error "signum not implemented"
>  abs    = error "signum not implemented"

This is an ordinary convolution routine tweaked to return the partial sum inner.

> convolve (a0:ar@(a1:as)) ~(b0:br@(b1:bs)) =
>  let (inner, _) = convolve ar br
>      ab = map (a0 *) bs
>      ba = map (* b0) as
>  in (a0*b0 : a0*b1+a1*b0
>            : zipWith3 (\a b c -> a+b+c) inner ab ba, 0 : inner)

The code is very similar to the usual power series multiplication routine. We can also use the same method described by McIlroy to divide our series.

As our series are a munged up version of the usual power series it's pretty surprising that it's possible to divide with so little code:

> instance Fractional a => Fractional (Partial a) where
>   fromRational n = let n' = fromRational n in Partial $ \_ -> (0, n' : repeat 0)
>   recip (Partial f) = Partial $ \x ->
>      let nibble (n, x:xs) | x==0      = nibble (n+1, xs)
>                           | otherwise = (n, (x:xs))
>          (n, xs) = nibble (f x)
>      in (-n, rconvolve x xs)

In effect, rconvolve solves the equation convolve a b==1:

> rconvolve x (a0:ar@(a1:as)) =
>   let (outer, inner) = convolve ar result
>       f a b = x*b-a
>       r = -1/f a0 a1
>       result = recip a0 : (map (r *) $ zipWith f outer inner)
>   in result

Note one ugly quirk of this code. I need to 'nibble' off leading zeroes from the series. This requires our underlying type a to have computable equality. (In principle we can work around this using parallel or).

That's it. We can now write a function to compute the positive part of a rational function. (By positive part, I mean all of the terms using non-negative powers of x.)

> pos f z = let Partial g = f $ Partial $ \x -> (1, 1 : repeat 0)
>               (n, xs) = g z
>           in if n>0
>               then z^n*head xs
>               else xs!!(-n)

Here are some examples:

> test1 = let f x = (1+2*x)/(3-4*x*x)
>         in pos (\x -> 1/(f x-f 0)/x) (0::Rational)

> test2 = pos (\x -> 1/(1+4*x+3*x^2+x^3) - 1/(1+x)) (1::Rational)

The original example I started with:

> test3 = pos (\x -> 1/(x+x^2) - 1/(x+2*x^2)) (0::Rational)

No division by zero anywhere!

Conclusions
The code works. But it does have limitations. As written it only supports rational functions. It's not hard to extend to square roots. (Try writing the code - it makes a nice exercise.) Unfortunately, any implementation of square root will (I think) require a division by x. This means that you'll be able to compute the positive part away from zero, but not at zero.

This method can't be extended fully to transcendental functions. But it is possible to add partial support for them. In fact, So with a little work we can still compute the positive part of functions like 1/sqrt(cos x-1) away from x==0. But applying cos to an arbitrary rational function may need more complex methods. I encourage you to experiment.

Note that this code makes good use of laziness. If your function has no singularities then you might find it performs no computations beyond what is required to compute the ordinary numerical value.

Automatic even/odd splitting

Statement of the Problem
Suppose you have a real valued function on the reals, say f. We can split it into the sum of an even and odd part:

f(x) = f_odd(x)+f_even(x)

where

f_odd(x) = (f(x)-f(-x))/2, f_even(x) = (f(x)+f(-x))/2

If f_odd has a power series around zero, then all of its terms must have odd powers in x. So f_odd(x)/x must have all even powers and it becomes natural to 'compress' down the terms to form f'(x) = f_odd(sqrt(x))/sqrt(x). (Take that as a definition of f' for this article.) We can implement this operation as a higher order function:

> import Prelude hiding (odd)

> odd f x = let s = sqrt x in 0.5*(f s - f (-s))/s

Here's a simple example:

> f x = x*(1+x*(2+x*(3+7*x)))

> test0 = odd f 3

But there's something not quite right about this. If f is rational, then so is odd f. But the implementation of odd involves square roots. Among other things, square roots introduce inaccuracy. As square roots don't appear in the final result, can we eliminate them from the intermediate steps of the computation too?

A Better Solution
Let's use the results of last week's article to compute odd another way. We want a linear function that maps as follows:

1	->	0
x	->	1
x2	->	0
x3	->	x
x2n	->	0
x2n+1	->	xn

Here's an automaton:

If we start at 0, end at 1, and take exactly n steps, then the product of the factors we collect up along the way is given by the second column of that table. In n is even there is no such path so we collect up 0. As we're working with polynomials over the reals, rather than types, we have x1=1x and so on. We can construct a transition matrix:

0	1
x	0

We now do something like we did last time. Any time we have a function of some variable x, we replace x with the transition matrix. Our functions now take matrix values like

a	b
c	d

Any polynomial of our transition matrix always gives us equal elements along the diagonal. This is true even if we form the inverse of the transition matrix. So we don't need to store d. So now we implement a simple matrix type needing only to store three elements instead of four:

> data O a = O a a a deriving (Show, Eq)

> instance Num a => Num (O a) where
>    O a b c + O a' b' c' = O (a+a') (b+b') (c+c')
>    O a b c * O a' b' c' = O (a*a'+b*c') (a*b'+b*a') (c*a'+a*c')
>    fromInteger n = let i = fromInteger n in O i 0 0
>    negate (O a b c) = O (negate a) (negate b) (negate c)

Notice how similar this is to automatic differentiation. We can extend to reciprocals too:

> instance Fractional a => Fractional (O a) where
>    fromRational n = let i = fromRational n in O i 0 0
>    recip (O a b c) = let idet = recip (a*a-b*c) in O (idet*a) (idet*negate b) (idet*negate c)

And now we can implement a replacement for odd:

> transition x = O 0 1 x
> odd' f x = let O _ fx _ = f (transition x) in fx

> test1 = odd' f 3

Another example:

> test3 = odd' (\x -> (x+3*x*x-1/x)/(x*x)) 2
> test4 = odd  (\x -> (x+3*x*x-1/x)/(x*x)) 2

This new version has many advantages: it uses only rational functions, it's more accurate, and it's well defined at zero.

Conclusion
Automatic differentiation is just one of a family of methods that can be used to compute a wide variety of functions of real-valued functions. Essentially we're just working over real-valued matrices instead of real numbers. By using automata we can simplify the process of working out which matrices to use. (Though for the simple example above, you may have been able to guess the matrix without any other help).

(BTW I think there are hidden automata lurking in a few places in mathematics. For example, in Umbral calculus.)