Monads for Markov chains

This is my first literate Haskell post, and actually my first post on Haskell, and my first attempt at literate programming. I spent (wasted?) a lot of time in trying to make the code segments render nicely on WordPress, with a process that was only partially automated and very tedious. I think I will cook up some script to automate it completely, if I ever get to write another such post.

Suppose you need to model a finite Markov chain in code. There are essentially two ways of doing that: one is to simply run a simulation of the Markov chain using a random number generator to obtain dice rolls and random cards from the decks, the other is to create a stochastic matrix containing the transition probabilities for each pair of states.

In this post I will show how a single monadic description of the Markov chain dynamics can be used to obtain both a simulator and the transition matrix.

> {-# LANGUAGE MultiParamTypeClasses, 
>     FlexibleInstances, 
>     GeneralizedNewtypeDeriving #-}
>
> import Control.Arrow
> import Control.Monad
> import Control.Monad.State.Strict
> import Data.Array
> import Random

Let’s start with an example of Markov chain and how we would like to be able to implement in Haskell. Consider a simplified version of the familiar Monopoly game: there are 40 squares (numbered 0 to 39), you throw two 6-sided dice each turn, some special squares have particular effects (see below), if you get a double roll three times in a row, you go to jail. The special squares are:

30: go to jail
2, 17, 33: Community Chest
7, 22, 36: Chance

Community Chest (CC) and Chance (CH) make you take a card from a deck and move to some other place depending on what’s written on the card. You will find the details on the code, so I won’t explain them here.

This is of course a Markov chain, where the states can be represented by:

> type Square = Int
> data GameState = GS {
>       position :: Square,
>       doubles :: Int } deriving (Eq, Ord, Show)

and a description of the game can be given in a monadic style like this:

> sGO :: Square
> sGO = 0
> 
> sJAIL :: Square
> sJAIL = 10
>
> finalize :: Square -> Game Square
> finalize n
>     | n == 2 || n == 17 || n == 33 = cc n
>     | n == 7 || n == 22 || n == 36 = ch n
>     | n == 30 = return sJAIL
>     | otherwise = return n
> 
> cc :: Square -> Game Square
> cc n = do i <- choose (1 :: Int, 16)
>           return $ case i of
>                      1 -> sGO
>                      2 -> sJAIL
>                      _ -> n
> 
> ch :: Square -> Game Square
> ch n = do i <- choose (1 :: Int, 16)
>           return $ case i of
>                      1 -> sGO
>                      2 -> sJAIL
>                      3 -> 11
>                      4 -> 24
>                      5 -> 39
>                      6 -> 5
>                      7 -> nextR n
>                      8 -> nextR n
>                      9 -> nextU n
>                      10 -> n - 3
>                      _ -> n
>     where
>       nextR n = let n' = n + 5
>                 in n' - (n' `mod` 5)
>       nextU n
>           | n >= 12 && n < 28 = 28
>           | otherwise = 12
> 
> roll :: Game (Int, Int)
> roll = let r1 = choose (1, 6)
>        in liftM2 (,) r1 r1
> 
> markDouble :: Bool -> Game ()
> markDouble True = modify $ \s -> s {
>                     doubles = doubles s + 1 }
> markDouble False = modify $ \s -> s {
>                      doubles = 0
>                    }
>
> goTo :: Square -> Game ()
> goTo n = let n' = n `mod` 40
>          in modify $ \s -> s { position = n' }
> 
> game :: Game ()
> game = do n <- liftM position get
>           (a, b) <- roll
>           markDouble (a == b)
>           d <- liftM doubles get
>           if d == 3
>            then do markDouble False
>                    goTo sJAIL
>            else do let n' = n + a + b
>                    n'' <- finalize n'
>                    goTo n''
> 

As you can see, Game is a state monad, with an additional function choose that gives us a random element of a range:

> class MonadState s m => MonadMC s m where
>     choose :: (Enum a) => (a, a) -> m a

This can be implemented very easily using the (strict) state monad and a random generator:

> newtype MCSim s a = MCSim (State ([s], StdGen) a)
>     deriving Monad
> 
> instance MonadState s (MCSim s) where
>     get = MCSim $ liftM (head . fst) get
>     put x = MCSim . modify $ \(xs, g) -> (x : xs, g)
> 
> instance MonadMC s (MCSim s) where
>     choose (a, b) = MCSim $
>                     do (xs, g) <- get
>                        let bnds = (fromEnum a, fromEnum b)
>                        let (y, g') = randomR bnds g
>                        put (xs, g')
>                        return . toEnum $ y
> 
> -- type Game a = MCSim GameState a
>
> runSim :: StdGen -> Int -> s -> MCSim s () -> [s]
> runSim g n start m = fst $ execState m' ([start], g)
>     where
>       (MCSim m') = foldr (>>) (return ()) $ replicate n m

The runSim function runs the simulation and returns the list of visited states. This is already quite nice, but the best thing is that the same code can be used to create the transition matrix, just swapping in a new implementation of the Game type alias:

> newtype MC s a = MC (s -> [(s, Double, a)])
> 
> instance Monad (MC s) where
>     return x = MC $ \s -> return (s, 1.0, x)
>     (MC m) >>= f = MC $ \s ->
>                    do (s, p, x) <- m s
>                       let (MC m') = f x
>                       (s', q, y) <- m' s
>                       return (s', p * q, y)
> 
> instance MonadState s (MC s) where
>     get = MC $ \s -> return (s, 1.0, s)
>     put x = MC $ \s -> return (x, 1.0, ())
> 
> instance MonadMC s (MC s) where
>     choose (a, b) = let r = [a..b]
>                         p = recip . fromIntegral . length $ r
>                     in MC $ \s -> map (\x -> (s, p, x)) r
> 
> type Game a = MC GameState a

The idea is that we keep track of all possible destination states for a given state, with associated conditional probabilities. For those familiar with Eric Kidd’s series on probability monads, this is basically:

type MC s a = StateT s (PerhapsT [] a)

Now, how to get a transition matrix from such a monad? Of course, we have to require that the states are indexable:

> markov :: Ix s =>
>           MC s () -> (s, s) -> Array (s, s) Double
> markov (MC m) r = accumArray (+) 0.0 (double r) $
>  	     	    range r >>= transitions
>     where
>       mkAssoc s (s', p, _) = ((s, s'), p)
>       transitions s = map (mkAssoc s) $ m s
> 	double (a, b) = ((a, a), (b, b))

So we iterate over all states and use the probability values contained in the monad to fill in the array cells corresponding to the selected state pair.

To actually apply this to our Monopoly example, we need to make GameState indexable:

> nextState :: GameState -> GameState
> nextState (GS p d) = if d == 2
>                      then GS (p + 1) 0
>                      else GS p (d + 1)
> 
> instance Ix GameState where
>     range (s1, s2) = takeWhile (<= s2) .
>                      iterate nextState $ s1
>     index (s1, s2) s =
>         let poss = (position s1, position s2)
>         in index poss (position s) * 3 +
>            doubles s - doubles s1
>     inRange (s1, s2) s = s1 <= s && s <= s2
>     rangeSize (s1, s2) = index (s1, s2) s2 + 1

then finally we can try:

> monopoly :: (GameState, GameState)
> monopoly = (GS 0 0, GS 39 2)
> 
> initialState :: Array GameState Double
> initialState = let n = rangeSize monopoly
>                    p = recip $ fromIntegral n
>                in listArray monopoly $ replicate n p
> 
> statDistr :: Int -> [(GameState, Double)]
> statDistr n = let mat = markov game monopoly
>                   distributions = iterate (.* mat)
>                         initialState
>                   st = distributions !! n
>               in assocs st

where .* is a simple vector-matrix multiplication function:

> infixl 5 .*
> (.*) :: (Ix i, Num a) =>
>            Array i a -> Array (i, i) a -> Array i a
> (.*) x y = array resultBounds
>               [(i, sum [x!k * y!(k,i) | k <- range (l,u)])
>                | i <- range (l'',u'') ]
>         where (l, u) = bounds x
>               ((l', l''), (u', u'')) = bounds y
>               resultBounds
>                 | (l,u)==(l',u') = (l'', u'')
>                 | otherwise = error ".*: incompatible bounds"

Calling statDistr 100 will return an association list of states with corresponding probability in an approximation of the stationary distribution, computed by applying the power method to the transition matrix. The number 100 is a pure guess, I don’t know how to estimate the number of iterations necessary for convergence, but that is out of the scope of this post, anyway.

Rewriting from scratch

I have recently come across this article (and those linked by it), while I was considering a rewrite for my own pet project.

While I tend to think that rewriting from scratch is usually a bad idea, if we are talking about open source projects, there’s another element you should take into account: how pleasurable it is to work on it. If you end up with a big and hairy code base you cannot touch without suddenly feeling all sad and pessimistic about the fate of the world, you should realize that something has gone deeply wrong: you are supposed to work on your project for fun (yeah, to scratch an itch, maybe, too…), and when the fun stops, most of the reasons to continue disappear. You have no (paying) customers, you have no deadlines, you will be able (or unable, depending on what you do in the rest of your time) to pay your bills anyway.

At this point, either you throw everything into the garbage bin – which may be sad, disappointing and frustrating, both for you and your users – or you take whatever lessons you may have learned in the process, and use your extra wisdom to start from scratch a hopefully better and more mature design and implementation that won’t have those shortcomings that brought you to this miserable situation.

It will probably have other defects and issues, though, more subtle and harder to spot in an upfront design. You have to accept it. And you have to accept the fact that you may need to start from scratch once again in the future.

Here is a simple form you can use to announce your decision to start from scratch. Just copy-paste it into your blog, filling in the required fields. It may save you precious time that you will be able to spend on your n-th rewrite

I'll rewrite my pet project from scratch because:

[ ] It's not modular enough
[ ] The original design didn't take __________________ into consideration
[ ] It's too slow
[ ] There is too much legacy code which I can't seem to get rid of
[ ] I want to write it in ____________________

----

[ ] And this time I'll use TDD/BDD

[ ] And this time it will be done right!

Signed (the developer)

_____________________________

Fighting zombies with a double fork

No, this post is not about a special weapon against Romero’s undeads (though that would be cool , but about a trick I’ve just discovered after a lot of struggling with child processes that aren’t properly reaped by their parent.
Let me try to be a bit more specific: suppose you have a parent that needs to spawn two kinds of child processes. Processes of the first kind perform a computation and return to the parent immediately (so they are handled via a fork+exec+waitpid combination), while those of the second kind should be started, their pid recorded so that they may be killed later, and then forgotten.
The question is: how to handle the second kind of subprocesses? If you just fork+exec, when you later kill them, you get zombie processes. If you ignore SIGCHLD (so auto-reaping is enabled) you can’t perform a blocking wait on the children of the first kind anymore. Handling SIGCHLD is perhaps possible, but I think you need to keep track of the PIDs of the spawned processes and do some magic in the handler. I didn’t manage to make this working.

What I did manage to get working is the following: just spawn the process in a temporary child process. This way the grandchild will be orphaned and won’t require a wait from its grandfather.

There are however some problems with this approach. Namely, the child needs to send the PID of its own child to the parent. I’ve accomplished that through a pipe. Here is some example code to better illustrate this technique:

from subprocess import Popen, PIPE, STDOUT
import sys

def spawn(cmd):
    cmdline = ('python', __file__) + cmd
    return int(Popen(cmdline, stdout = PIPE).stdout.read())

if __name__ == '__main__':
    print Popen(sys.argv[1:], stdout = open('/dev/null', 'w'), stderr = STDOUT).pid
    sys.exit(0)

skema 0.1

I have just released skema, a very simple tool I wrote to automate some repetitive coding tasks. Nothing fancy, just a ruby script that expands ERB templates in a directory, reading variable values from the command line, configuration files, or interactively.
Its purpose is similar to kapptemplate (it already includes a minimal KDE application template), but it may be faster to use for little things, and it is probably easier to create templates for it.
Here is a tarball for those of you who want to try it out. You can find a short tutorial in the README file that should be enough to get you started.
Here is a screenshot of skema in interactive mode while expanding the kapp template:

Some thoughts on git

As I mentioned in a previous post, we started using git for managing Tagua source code. There’s currently a lot of controversy in the topic of distributed versus centralized VCS’s, and I feel like expressing my own ideas, too.
I’m no git guru (yet), so please don’t get offended if I missed something in my criticism. I hope to generate a fruitful discussion, because I really think that the VCS is an important element in the development of an open source project.

Why git rules

• It’s decentralized. I fully agree with Linus on this point: the decentralized model is superior, for various reasons I won’t talk about in detail here.
• It’s simple. I don’t understand why people keep complaining about git being overly complicated. The object model is fairly straightforward, as soon as you stop thinking about revisions the cvs/svn way, and branching/merging is trivial.
• The object model is very well designed, flexible and useful.

Why git sucks

Well, git has a number of minor (?) defects, like portability issues, impossibility of doing partial checkouts (actually I think this is a problem with all the decentralized VCS’s) and stuff like that, but let me concentrate on those issues which are fundamental to git as an open source project, and are unlikely to be addressed in the future:

• It is written in (unreadable) C. I don’t want to start a flame war (really!), but browsing through git’s code is a painful experience. Huge .c files with tons of big functions with no apparent structure and almost no comment. Ubiquitous micro-management stuff for handling manually allocated buffers. In other words… a mess.
• It (ab)uses the UNIX philosophy. Do one thing, and do it well. Yeah, of course it’s a great idea, but maybe they shouldn’t have taken it too literally. I mean, if we talk about modularity and reusable components, I’m all for it, but it doesn’t mean you have to make a different executable for each task (even low level ones that should be invisible to the end user) and chain them together via text pipes and bash scripts!
• It’s not reusable by third party applications, at least not in any efficient/convenient way. If you want to build a higher-level abstraction over git (a so-called porcelain), you need to fork out every now and then to get input from git, and then parse its textual output. Yes, I know about libgit.a, but that’s not really usable as a library (half of git functions call exit on failure) and has no API documentation at all.

That said, I have to admit that most of its defects aren’t that important if you are just going to use it, and not develop it. However, using something you have problems inspecting and modifying feels a bit like using a closed source product…

KBoard is dead. Long live Tagua!

We finally made a decision, and named the project Tagua. Thanks to all the people who suggested possible names, and expecially to Riccardo Iaconelli who came up with the definitive one.
I’m now doing a global s/kboard/tagua/g on webpages, code, tracker, cron jobs, git repository…

Ah, yes, I forgot to mention that we moved the project to git when we resumed the development at Akademy. There has been a long thread on kde-core-devel discussing the adoption of git by subprojects, and there emerged that moving away from the KDE svn (though apparently discouraged) does not mean parting from the KDE project. I’ll probably talk about our experience with git in a dedicated post.

So, if you want to try the development version of Tagua, either grab a nightly snapshot (temporarely located here), or clone our main repository:

git clone http://repo.or.cz/r/tagua.git

KBoard: an update

KBoard is a project by Maurizio and myself for a generic board game application. It started back in October 2005, and progressed slowly but steady until the end of 2006. After then it was left starving on the playground, basically untouched until a few days ago.
After a prolific discussion at Akademy, we are now back developing KBoard at full speed. Here are a couple of screenshots showing our progresses so far:

Compare them with the old screenshots and you’ll see how much has been done in just a few weeks: the main layout and the clock are now themable via Lua scripts, just like the chessboard and pieces.
I shall thank Maurizio for this astonishing work, and generally for his invaluable Lua theme loader, which is one of the finest pieces of code in KBoard.

By the way, most of the recent changes are actually invisible to the user, and concern the API used by the graphical interface to communicate with game plugins (called variants). I recently worked (not that much) on the animation stuff, and ported three games to the new API. Maurizio made the whole refactoring of the code using the old API at akademy, and revised the concepts around the pool (i.e. the place where captured pieces end up in games like Shogi and Crazyhouse).

We plan to release KBoard 1.0 by the end of the year. It might be an optimistic estimate, but if we continue to work with this speed, it should be actually possible. There is a lot of things to fix, but the planned features are almost all there. Needless to say, any help (developers, artists) would be really appreciated.

At the moment, our biggest dilemma is the name. We decided that the name should be changed (KBoard is ugly and doesn’t even make clear what the application is about), but after a whole day of awfully poor proposals, we still have to find a decent alternative. So you if have a nice name in mind (not necessarily starting with k, think about Phonon, Plasma, Solid… those are good names!) please tell us!

Canvas adapter and plasma

I’ve finally finished polishing the code of KGameCanvasAdapter, a hundred lines of straightforward code that stand as a bridge between KGameCanvas code and any QPainter based drawing system, including of course QGV.

The only test case for the adapter is a port of my kollision game to plasma. Porting was very easy: it was just a matter of having the MainArea class inherit from Plasma::Applet and KGameCanvasAdapter instead of QWidget. The resulting applet is not really a game, at the moment: it is just a box with bouncing red balls.
However, I think it’s a good test for both my adapter and plasma, since it can show how well plasma performs with CPU intensive plasmoids. Here is a short screencast of the kollision plasmoid in action.

The code is available in a git repository:

git clone http://kollision.ath.cx/plasma.git plasma

will download the whole patched plasma directory. To compile it, you should move it inside a kdebase/workspace working copy, so that it uses the existing cmake build system of kdebase.

KBattleship: almost ready

This week I’ve been hacking quite a lot on KBattleship, adding almost all required features: a decent AI, network play and sounds. While Riccardo is working on KWelcomeScreen, a library class inspired by this blog post by Johann, I’m going to make some minor adjustments, add small missing features like scoring, and if everything goes well I guess we can move the application to kdereview within next week.
This is a screenshot showing a game between the old KBattleship and my rewrite:

Dynamic types and virtual inner types

The concept based polymorphism I tried to explain in my previous post can be more conveniently expressed with another syntax which avoids concepts at all.

The key observation is that runtime polymorphism and F-bounded polymorphism are very similar in nature, and, to some extent, the latter can be implemented using the former. If I have understood correctly, that is exactly how the Java virtual machine implements generics.

For example, suppose you have an abstract class Shape with a virtual function draw, and several concrete classes like Square, Circle, Triangle, which override draw.
Client code using runtime polymorphism would look like:

void render(Shape* shape)
{
  shape->draw();
}

while using generic programming one could write:

template <typename S>
where Inherits<S, Shape>
void render(S* s)
{
  s->draw();
}

assuming that the concept Inherits<T, U> is verified when T inherits from U (possibly indirectly). An approximation for Inherits is

auto concept Inherits<typename T, typename U>
{
  // T can be implicitly converted to U
  T::operator U();
}

Using the generic code in such situations should probably be considered a mistake, because if one wants to use template based polymorphism, defining the abstract class Shape is pointless. It would be better to directly define a ShapeConcept having a member function draw.

However, the example suggests that a syntax resembling templates could be used to express runtime polymorphism:

template <Shape! S>
void render(S* s)
{
  s->draw();
}

the idea being that S is the runtime type of s. The compiler could instantiate the template immediately and implement the body just like a Shape pointer were used instead of S. Here S is not really a type, since for example expressions like

new S

should be rejected; let’s call things like S dynamic types. They behave much like incomplete types (cannot be instantiated directy, sizeof cannot be called on them), but are a different beast: for the purpose of typechecking they are just subtypes of their parent concept (Shape, in this example), but for code generation they are considered equivalent to it.

Why could all of this syntax be possibly useful? Well, when using ordinary abstract classes, it does not really add anything to the language, but it shows its power when another (more than syntactical) extension is introduced: virtual inner types.

Just like the name suggests, virtual inner types are inner types which behave polymorphically. The typical (and motivating) example here is still that of the abstract factory. So suppose you have the following abstract factory:

class AbstractAnimal;
class AbstractFood;

class AbstractFactory
{
public:
  virtual AbstractAnimal* createAnimal() = 0;
  virtual AbstractFood* createFood() = 0;
};

and the following abstract products:

class AbstractAnimal
{
public:
  virtual void eat(Food*) = 0;
};

class AbstractFood
{
public:
  virtual int quality() const = 0;
};

As explained in my previous posts [1] [2], there’s no way to explain to the compiler that the Food that an animal is able to eat is only that created by the same factory which created that very animal.
Virtual inner types could be used to solve the problem:

class AbstractAnimal;
class AbstractFood;

class AbstractFactory
{
public:
  virtual typedef AbstractAnimal Animal;
  virtual typedef AbstractFood Food;
  virtual Animal* createAnimal() = 0;
  virtual Food* createFood() = 0;
};

class AbstractAnimal
{
public:
  virtual typedef AbstractFactory Factory;
  virtual void eat(Factory::Food* food) = 0;
};

class AbstractFood
{
public:
  virtual typedef AbstractFactory Factory; // unused, only here for completeness
  virtual int quality() const = 0;
};

The idea is that virtual typedef Base X forces inherited class to define an inner type X inheriting from Base. Furthermore, this type behaves polymorphically (the type itself, not only instances!). That means that if A is a dynamic type with parent AbstractAnimal, then A::Factory is the correct dynamic type with parent AbstractFactory.

To explain this fact in more detail, let’s continue the example:

class Cow;
class Grass;

class CowFactory : public AbstractFactory
{
public:
  typedef Cow Animal;
  typedef Grass Food;
  Cow* createAnimal();
  Grass* createFood();
};

class Cow : public Animal
{
public:
  typedef CowFactory Factory;
  void eat(Grass* grass);
};

class Grass : public Food
{
public:
  typedef CowFactory Factory;
  int quality() const;
};

notice that we have defined inner types corresponding to the virtual typedefs of parent classes. Besides, the argument to eat is Grass and not AbstractFood.
Client code, in fact, cannot simply pass an AbstractFood object to the eat member function. The following function should not compile:

void feed(AbstractAnimal* animal, AbstractFood* food)
{
  animal->eat(food);
}

but can be fixed with the use of dynamic types

template <AbstractFactory! Factory>
void feed(Factory::Animal* animal, Factory::Food* food)
{
  animal->eat(food);
}

expressing the fact that both animal and food come from the same factory of dynamic type Factory. The food parameter (actually stored in an AbstractFood variable), is statically downcasted to Grass before being passed to the eat member function. The static typechecking via dynamic types ensures that such a downcast is always safe, so there is no need to perform a dynamic cast.

These ideas for a language extension are essentially equivalent to the ones expressed in the previous post, but seem to require less drastic changes in compiler implementations. Furthermore, they don’t use concepts, so it could be even implemented on top of the existing gcc.