Postmodern Haskell and OpenGL: Introducing vinyl-gl

July 2nd, 2013 3 comments

This article is for those familiar with Haskell and, at least passingly, with Jon Sterling’s vinyl library (that packs quite a nice introduction). Further, it is assumed that the reader is familiar with the basics of computer graphics and OpenGL. I have written another article that introduces the use of relatively modern OpenGL practice with Haskell that may serve as a primer for this article.

This post is a guided tour through the Game2D example in the vinyl-gl repository. I think (hope) the examples directory there is a useful starting point to see how this library may be used. If you want to skip the examples and just get coding, the library is available on hackage.

What’s the Problem?

A drawback to the changes undergone by OpenGL over the past few years as it has shifted from a fixed-function to a programmable pipeline is that it now seems much harder to get started. Part of the problem is that, while moving virtually all processing into shader programs provides better performance and vastly more flexibility, programmers must now connect their application code to their GLSL code. This is a bit of a pain as it typically requires thinking about detailed aspects of memory layout, and rather pedantic interactions with OpenGL in order to get past the fact that you must hand it a void* with the promise that this blob of memory is full of useful goodies.

A Worked Example

A first example of connecting Haskell OpenGL code with GLSL is the skeleton of a 2D side-scroller game (think Super Mario Bros.). In these games, levels are often drawn as a set of reusable tiles, and player movement causes a camera to translate across the level. Here’s what it will look like,

I will begin by presenting the GLSL code to be used. I am going to be drawing textured 2D triangles to make up the game level. Our vertex shader will take in a 3×3 matrix that captures the 2D transformations our geometry may undergo. A vertex coordinate and a texture coordinate will comprise the vertex input for the shader, which will then output clip-space coordinates for the vertex data (gl_Position), and pass the texture coordinates through to the fragment shader.

 #version 150
uniform mat3 cam;
in vec2 vertexCoord;
in vec2 texCoord;
out vec2 texCoordFrag;

void main() {
  texCoordFrag = texCoord;
  gl_Position = vec4(cam * (vec3(vertexCoord, 1) * 2 - 1), 1);

The fragment shader simply samples the texture.

 #version 150
uniform sampler2D tex;
in vec2 texCoordFrag;
out vec4 fragColor;

void main() {
  fragColor = texture(tex, texCoordFrag);

In summary, we have uniform inputs cam and tex. We also have per-vertex inputs vertexCoord and texCoord.

Now we can write some Haskell to feed the GLSL monster. We start by pulling in vinyl, GLUtil, and vinyl-gl in addition to some of the usual suspects.

{-# LANGUAGE DataKinds, TypeOperators #-}
import Control.Applicative
import Data.Foldable (foldMap, traverse_)
import Data.Vinyl
import Graphics.GLUtil
import Graphics.GLUtil.Camera2D
import Graphics.Rendering.OpenGL
import Graphics.UI.GLFW (Key(KeyEsc))
import Graphics.VinylGL
import Linear (V2(..), _x, M33)
import System.FilePath ((</>))

We also have some helper modules used in these examples. These modules are tied to the GLFW-b package’s windowing and input facilities. They serve to kick off our main window with an event loop, and map keyboard input to camera movements.

import Keyboard2D (moveCamera)
import Window (initGL, UI(..))

The rendering parts of our application need to refer to some common rendering state. In our case, we just have a 2D transformation matrix that operates on homogenous 2D points, but we use a vinyl record to get used to techniques that let us avoid coupling the rendering state across all rendering functions. Take a look Demo3D.hs example and its associated rendering functions in Geometry.hs to see this in action. Note the association between the field named "cam" and the vertex shader uniform named cam.

type AppInfo = PlainRec '["cam" ::: M33 GLfloat]

Our game is just going to have ground in front of a blue sky background. The ground is made up of columns of dirt tiles, each topped by a grass tile. Given that design, a level in our game is simply a list of column heights. Here we work with game maps made up of columns whose height ranges from 0 to 10.

gameLevel :: [Int]
gameLevel = [3,3,3,4,5,4,3,3,3,4,5,5,6,7,6,6,6,7,6,5,4,3,3]

We will need to convert a tile height to a 2D square in order to draw the tile on the screen. We define this mapping once so we have a compact game level representation and a consistent-by-construction tile size. Since we are only given the height of the tile to produce, we set its left edge to have an X coordinate of 0. We also here convert our game level’s vertical coordinate system of [0,10] to normalized coordinates that range from 0 to 1. The tile function gives us the four vertices of a square for a tile.

tile :: Int -> [V2 GLfloat]
tile h = let h' = fromIntegral h / 10 in V2 <$> [0,0.2] <*> [h', h' - 0.2]

The tile function can be used to build a column of tiles whose left edge has an X coordinate of 0. We now lay these columns out side-by-side extending from a leftmost edge at the Y axis. Each successive column is placed farther along the positive X axis.

spaceColumns :: [[V2 GLfloat]] -> [[V2 GLfloat]]
spaceColumns = zipWith (map . (_x +~)) [0, 0.2 ..]

Building a Bridge to GLSL

The 2D vertices we produce will consist of not only 2D positions, but also texture coordinates that define how our tile images are mapped to 2D tiles. Having this Tex field means that we have the flexibility to switch to a sprite atlas technique, or to repeat an image across a tile. Once again, we are arranging a coincidence of naming here: the vinyl Field pos has a type with name "vertexCoord" that matches up with the vertexCoord our vertex shader expects as input.

type Pos = "vertexCoord" ::: V2 GLfloat
type Tex = "texCoord"    ::: V2 GLfloat

pos :: Pos
pos = Field

tex :: Tex
tex = Field

We’re not going to do anything fancy with texture coordinates for now, so we will simply assign each vertex texture coordinates that map the entire texture to the tile square.

tileTex :: [[V2 GLfloat]] -> [PlainRec [Pos,Tex]]
tileTex = foldMap (flip (zipWith (<+>)) (cycle coords) . map (pos =:))
  where coords = map (tex =:) $ V2 <$> [0,1] <*> [0,1]

We compute all the grassy tiles from our gameLevel data structure together because they share the same texture. By grouping our tiles by their textures, we can simplify eventual rendering.

The grassTiles action computes a list of vertices with position (Pos) and texture coordinate (Tex) fields, then loads all that data into OpenGL. The fields used here will correspond to the inputs expected by our GLSL vertex shader. Data is loaded into OpenGL with the bufferVertices function that takes a list (or a Data.Vector.Storable.Vector) of PlainRecs and feeds the data into OpenGL. Note that the BufferedVertices type is tagged with the fields of the buffered vertex data. This lets us compare our buffered data with the expectations of our GLSL program when possible.

grassTiles :: IO (BufferedVertices [Pos,Tex])
grassTiles = bufferVertices . tileTex . spaceColumns $ map tile gameLevel

Producing the columns of dirt that support our grassy tiles is slightly more complicated. Conceptually, we generate a new tile at every possible tile height down to 1 (remember our game level has tiles going up from 0). This means that rather than producing one tile for each element of the gameLevel structure, we instead produce a list of tiles corresponding to a dirt column in our game world.

dirtTiles :: IO (BufferedVertices [Pos,Tex])
dirtTiles = bufferVertices . tileTex . spaceColumns $ map col gameLevel
  where col :: Int -> [V2 GLfloat]
        col h = foldMap tile [h - 1, h - 2 .. 1]

So much for geometry! We have now defined how to produce a bunch of textured 2D squares from our gameLevel data structure. The next step is to load the images we want to use into OpenGL TextureObjects. Here, we load images from the art directory, and set each to use nearest-neighbor filtering.

loadTextures :: [FilePath] -> IO [TextureObject]
loadTextures = fmap (either error id . sequence) . mapM aux
  where aux f = do img <- readTexture ("art" </> f)
                   traverse_ (const texFilter) img
                   return img
        texFilter = do textureFilter Texture2D $= ((Nearest, Nothing), Nearest)
                       texture2DWrap $= (Repeated, ClampToEdge)

Ready to Render

Finally, we can define how to draw our game world! We will make use of the geometry and texture loading pieces defined above to produce a function that will draw our game level. I obtained the images used here from; with thanks to the artist,

background :: IO (AppInfo -> IO ())
background = 
  do [grass,dirt] <- loadTextures [ "ground.png", "ground_dirt.png" ]

With our textures loaded, we are now going to load our GLSL program by specifying a vertex shader and a fragment shader.

     s <- loadShaderProgram ("etc"</>"game2d.vert") ("etc"</>"game2d.frag")

Our fragment shader has a sampler2D uniform input to support texturing. We are always going to use the first texture unit, denoted 0, so we can set this uniform parameter of our shader program now.

     setUniforms s (texSampler =: 0)

We are now ready to prepare our geometry. Recall that we have code that loads a bunch of 2D vertices into OpenGL: four vertices for each tile. We are going to render each tile as a pair of triangles, but we don’t want to make separate rendering calls for every tile, or even every triangle. To this end, we will also provide OpenGL with vertex indices that describe how to index into the buffered vertices to produce triangles.

We have defined our square tile vertices like this:

0 ------- 2
|         |
|         |
|         |

These four vertices may be drawn as two triangles with counterclockwise winding by considering one triangle consisting of vertices [0,1,2] and another consisting of vertices [2,1,3]. Thus the indices [0,1,2,2,1,3] may be used to index into a buffer of four vertices to render two triangles making up one square tile. We’ve buffered all our tile vertices contiguously, so all we need to do is repeat the [0,1,2,2,1,3] pattern for each tile, shifting the sequence by 4 for each tile (for example, the second tile indexes into the buffered vertex data at [4,5,6,6,5,7]). We build an infinite list of such shifted sequences in the definition for inds below, then take 6 indices for each tile in our game world.

     grassVerts <- grassTiles
     eb <- bufferIndices inds
     grassVAO <- makeVAO $ do enableVertices' s grassVerts
                              bindVertices grassVerts
                              bindBuffer ElementArrayBuffer $= Just eb
     dirtVerts <- dirtTiles

OpenGL is a stateful system, and Vertex Array Objects (VAOs) were added to help encapsulate some of this state. In particular, we need to bind each field of our vertices to its corresponding GLSL shader input. We do this by giving OpenGL offsets and strides into our vertex buffer that define how to pick out each field of each vertex. The vinyl-gl library handles all of this automatically. The convention that drives the entire arrangement is that the symbolic name attached to each field is identical to the corresponding GLSL input’s name.

     dirtVAO <- makeVAO $ do enableVertices' s dirtVerts
                             bindVertices dirtVerts
                             bindBuffer ElementArrayBuffer $= Just eb
     return $ \i -> do currentProgram $= Just (program s)
                       setUniforms s i
                       withVAO grassVAO . withTextures2D [grass] $
                         drawIndexedTris numGrassTris
                       withVAO dirtVAO . withTextures2D [dirt] $
                         drawIndexedTris numDirtTris
  where numGrassTris = fromIntegral $ 2 * length gameLevel
        numDirtTris = fromIntegral . sum $ map (*2) gameLevel
        texSampler = Field :: "tex" ::: GLint
        inds = take (sum $ map (*6) gameLevel) $
               foldMap (flip map [0,1,2,2,1,3] . (+)) [0,4..]

The Window and Main Loop

Our application is driven by a single initialization function. We set the color to which each frame will be cleared, turn on alpha blending, as our textures have alpha channels, and initialize our background renderer.

setup :: IO (AppInfo -> IO ())
setup = do clearColor $= Color4 0.812 0.957 0.969 1
           blend $= Enabled
           blendFunc $= (SrcAlpha, OneMinusSrcAlpha)

The game loop makes use of a frame-tick action that provides an updated UI value and the drawing function returned from setup.

loop :: IO UI -> IO ()
loop tick = setup >>= go camera2D
  where go :: Camera GLfloat -> (AppInfo -> IO ()) -> IO ()
        go c draw = 
          do ui <- tick
             clear [ColorBuffer, DepthBuffer]
             let mCam = camMatrix c
                 info = Field =: mCam
             draw info
             if keysPressed ui ^. contains KeyEsc
             then return () -- terminate
             else go (moveCamera ui c) draw

Open the window and kick off the loop!

main :: IO ()
main = usage >> initGL "2D Platformer" 640 480 >>= loop
usage :: IO ()
usage = putStrLn "Arrow keys to translate, shift+arrow to rotate, esc to exit!"

And that’s it! Since walking through working code isn’t entirely satisfying, what would happen if we screwed up the synchronization between Haskell and GLSL?

Safety Net

Suppose we had written:

type Pos = "vertexCoord" ::: V2 GLint

and treated vertex position as an integer value throughout.

This would result in a runtime error at startup:

Game2D: Type mismatch in vertexCoord

Similarly, had we written:

type Pos = "vertexCoord" ::: V3 GLfloat

We would see the same error at runtime. We are not inspecting the GLSL program at compile time, as it may change between the time we typecheck our Haskell program and run our executable. However, once we have loaded GLSL code at runtime, and tried to enable some BufferedVertices with enableVertices', vinyl-gl can detect mismatches between the Haskell and GLSL sides.

Another potential mishap: had we let our names get out of sync and written,

type Pos = "vCoords" ::: V2 GLfloat

We would get a different error at runtime:

Game2D: GLSL expecting vertexCoord

These hypothetical errors may seem a bit contrived, but easily arise when composing larger applications from small pieces, as we always do in Haskell. The extensible records provided by vinyl combined with the vinyl-gl machinery let us,

  • Append fields to records to grow complex vertex descriptors
  • Project out simpler supertypes from more richly defined descriptors
  • Load complex vertex data into OpenGL, but only map parts of the data to GLSL parameters

This flexibility helps the programmer avoid getting boxed in by specific record types while maintaining good relations between Haskell and GLSL code.

Short and Sweet

To bring things full circle, here is a version of the code from the older article linked at the top compared to a version ported to use vinyl-gl and the newest GLUtil. We have eliminated nearly half the code by reducing OpenGL boilerplate (GLSL ↔ Haskell synchronization, VertexArrayDescriptor specifications, etc.), and better encapsulating rendering actions.

(Old version on the left, new version on the right.)

Categories: Programming Tags: ,

Fast Factory, Video 1

April 13th, 2013 Comments off

I spent the past few months developing the perception system behind this video of a PR2 grasping moving objects.

Perception is driven entirely by the Kinect mounted on the robot’s head, and all processing of the Kinect data is done in Haskell. I will write up a report on that side of things in the near future, but suffice to say that Haskell is more than ready for top flight performance and reliability (we’ve done about a month of testing, where a test day generally involves launching the perception software, then letting it run all day). It’s multithreaded, networked, and GPU accelerated. Free monads were involved. Lenses were involved. Burritos played a crucial role.

GLUtil on Hackage

August 31st, 2012 Comments off

The GLUtil Haskell package I wrote about before is now available on hackage.

Categories: Uncategorized Tags:

Introducing CLUtil

December 22nd, 2011 Comments off

Edited: May 4, 2013

Edited: September 6, 2013

  • Change CLUtil links to a tagged version compatible with the code presented here.


OpenCL is a cross-platform parallel programming standard with support for execution on both CPUs and GPUs. The OpenCL package on hackage provides a direct binding to the API with just enough Haskellosity to make invoking those API functions borderline pleasant. That said, there remains a certain amount of boilerplate that is rather offputting.

For context, we are trying to write program fragments that may be executed in parallel to improve efficiency. Multi-core CPUs are good at this, GPUs can be even better. We can bend GLSL to our needs (see my modern OpenGL with GLSL in Haskell tutorial), but are then left to contend with various graphics and GUI considerations that are completely irrelevant to our goals. OpenCL represents a much more familiar looking programming environment (a kind of C with support for small vectors), with less setup overhead than a GPGPU program using GLSL on the backend.

So let’s get started: How can we use OpenCL from a Haskell program?

First we’ll write a simple OpenCL kernel that adds two vectors, a and b, and stores the element-wise square of the result in c. To make it slightly interesting we’ll operate on four floating point numbers at a time.

__kernel void vecAdd(__global float *a, 
                     __global float *b,
                     __global float *c) {
  size_t i = get_global_id(0);
  float4 r = vload4(i,a) + vload4(i,b);

I’ve saved this code in a file called

To recap: we have an OpenCL kernel that takes two read-only input vectors and a write-only output vector. This sort of configuration fits rather well into a functional setting, so it is where we will focus our efforts as we endeavour to make the common things easy. Now we’re ready for some Haskell:

> import Control.Parallel.CLUtil
> import qualified Data.Vector.Storable as V

I’ve imported my CLUtil library that provides some simple helpers for common Haskell-OpenCL operations. Now let’s write a program that initializes OpenCL, loads and compiles an OpenCL kernel stored in a separate file, and run that kernel on some data.

> test1 = do s <- ezInit CL_DEVICE_TYPE_CPU
>            k <- kernelFromFile s "" "vecAdd"
>            let v1 = V.fromList [1,2,3,4::Float]
>                v2 = V.fromList [5,6,7,8::Float]
>            v3 <- runKernel s k v1 v2 (Out 4) (Work1D 1)
>            print (v3::Vector Float)

And that’s it! The ezInit function intializes an OpenCL device and sets up a context and queue, all of which are returned in a record. The kernelFromFile function reads OpenCL source code from a file, then builds the named kernel using the previously initialized OpenCL state record. Finally, we produce some data in Haskell and call runKernel.

The runKernel action is variadic in an attempt to meet most of your kernel-running needs. Its arguments are an OpenCL state record, a kernel to run, the parameters to be passed to the OpenCL kernel, and the number of global work items (i.e. how many times to invoke the kernel). Recall that our kernel had three parameters, and operated on four elements at once. So, we pass the kernel the two input vectors, a note requesting an output vector of length four, and a specification that our kernel should be run exactly once. We bind the output to another vector, and provide a type annotation on the next line to let runKernel figure out how much memory to allocate for the output vector. This is potentially confusing, hopefully more examples will clarify usage.

Next, let’s consider another kernel that demonstrates a bit more variety in the types.

__kernel void funnyBusiness(__global double *a,
                            __global double *b,
                            __global double *c,
                            __global int *s) {
  size_t i = get_global_id(0);
  double4 r = vload4(i,a) + vload4(i,b);
  double4 rsq = r * r;
  vstore4(rsq, i, c);
  s[i] = (int)(rsq.x + rsq.y + rsq.z + rsq.w);

This kernel is similar to the first, but uses the double scalar type, and also takes an int pointer as an argument that has an int cast of the horizontal sum of the result vector written to it. The Haskell side looks very similar to the last time,

> test2 = do s <- ezInit CL_DEVICE_TYPE_CPU
>            k <- kernelFromFile s "" "funnyBusiness"
>            let v1 = V.enumFromN  0 12 :: Vector Double
>                v2 = V.enumFromN 12 12 :: Vector Double
>            (v3,v4) <- runKernel s k v1 v2
>                                 (Out 12) -- v3 is a 12-element output
>                                 (Out 3)  -- v4 is a 3-element output
>                                 (Work1D $ 12 `div` 4)
>            print $ V.sum (v3::Vector Double)
>            print $ V.sum (v4::Vector CInt)

A few things to note:

  • v3 and v4 have different types. The runKernel action always returns Vectors, be they individual as in test1, or tupled as here. The CLUtil library currently supports up to three output vectors, but this is an arbitrary limit that will be raised as the library matures.

  • Once again we operate on four input elements at a time. Given that our source vectors have twelve elements each, we must run three work items.

  • The type annotation on v4 is CInt rather than Int. This is important as the size of Haskell’s Int may be 32 or 64 bits, while OpenCL’s int is always 32 bits.

One last bit of mystery that remains is the Work1D specification. OpenCL supports up to three dimensions of work items, which can ease the handling of multi-dimensional data. If we were adding a matrix, for example, we might have the elements of the matrices stored in Vectors, but the elements are intuitively addressed by a pair of coordinates. A naive OpenCL kernel for adding matrices might go something like this,

__kernel void floaty(__global float *a, 
                     __global float *b,
                     __global float *c) {
  size_t x = get_global_id(0);
  size_t y = get_global_id(1);
  size_t i = y * get_global_size(1) + x;
  c[i] = a[i] + b[i];

On the OpenCL side, we’re not gaining anything by working with two dimensions of work item indices, but sometimes these coordinates can inform coordinate-dependent logic in the OpenCL program.

The Haskell side is once again very similar, with the exception of the new Work2D specification.

> test3 = do s <- ezInit CL_DEVICE_TYPE_GPU
>            k <- kernelFromFile s "" "floaty"
>            let v1 = V.enumFromN  0 16 :: Vector Float
>                v2 = V.enumFromN 16 16 :: Vector Float
>            v3 <- runKernel s k v1 v2 (Out 16) (Work2D 4 4)
>                  :: IO (Vector Float)
>            mapM_ (\i -> print (V.slice i 4 v3)) [0,4..12]

Note that here we explicitly request a GPU to run our program. The GPU in my laptop supports floats, but not doubles, so I put this kernel in a separate file from the others.

Hopefully usage of CLUtil is becoming clear. So let’s crank things up and take advantage of the Haskell ecosystem.

A blog post describing quasicrystal figures serves as an interesting test of OpenCL performance. As a first pass, we can port the pixel computations from that post to OpenCL without thinking about optimization…

float wave(float,float,float,float);
float wrap(float);

__kernel void quasiCrystal(int pixels, float scale, float phase, int numAngles, __constant float *angles, __global uchar *img) { int x = get_global_id(0); int y = get_global_id(1); float denom = (float)pixels - 1; float u = scale * ((2 * (float)x / denom) - 1); float v = scale * ((2 * (float)y / denom) - 1); float sum = 0.0f; for(int i = 0; i < numAngles; ++i) { sum += wave(phase, angles[i], u, v); } uchar r = (uchar)(255.0f * clamp(wrap(sum), 0.0f, 1.0f)); vstore4((uchar4)(255,r,128,r), y*pixels + x, img); }
float wrap(float n) { int k = (int)n; if(n < 0) k -= 1; float v = n - (float)k; if(k%2 == 1) v = 1 - v; return v; }
float wave(float phase, float theta, float x, float y) { return (cos(cos(theta) * x + sin(theta)*y + phase) + 1) / 2; }

…then write the control logic in Haskell. We also tie our image generation into the gloss1 library for a real-time display.

> import Control.Parallel.CLUtil
> import qualified Data.Vector.Storable as V
> import Data.Word (Word8)
> import Graphics.Gloss hiding (Vector, scale)
> import Graphics.Gloss.Interface.IO.Animate hiding (Vector, scale)
> pixels :: Int
> pixels = 800
> scale :: Float
> scale = 128
> angles :: Int -> Vector Float
> angles n = (* (pi / fromIntegral n)) $ V.enumFromN 0 n
> mkPicture :: Vector Word8 -> Picture
> mkPicture = flip (bitmapOfForeignPtr pixels pixels) False
>           . (\(x,_) -> x)
>           . V.unsafeToForeignPtr0
> main = do s <- ezInit CL_DEVICE_TYPE_ALL
>           k <- kernelFromFile s "" "quasiCrystal"
>           let numPix = pixels * pixels
>               pixels' = fromIntegral pixels :: CInt
>               numAngles = 7
>               allAngles = angles numAngles
>               frame :: Float -> IO Picture
>               frame phase = mkPicture `fmap`
>                             runKernel s k pixels' scale phase
>                                       numAngles allAngles
>                                       (Out (numPix*4)) (Work2D pixels pixels)
>           animateIO (InWindow "Quasicrystal" (pixels,pixels) (0,0)) black frame

On my laptop, running on the CPU consumes 300% of a dual core i5 (with hyperthreading), while running on the GPU consumes 85% CPU. While this example is somewhat more complicated, the ability to switch between CPU and GPU computation driving a visualization, all from an interactive interpreter is highly compelling. The animation runs at a reasonable clip on a dual core CPU, so try it out yourself! It should look something like this,

Quasicrystal animation screenshot

Quasicrystal animation screenshot

These programs may all be found in the examples directory of the CLUtil package.

P.S. When building executables on your system, you may need to supply GHC with an -lOpenCL option. This is not necessary in Mac OS X, where we use the installed OpenCL framework, but is likely needed elsewhere.

  1. Visit the gloss homepage to get the latest version via darcs. You may need to manually relax version constraints on OpenGL and GLUT to install gloss with the latest versions of those packages. Note also that on Mac OS X, installing gloss with the GLFW-b backend (cabal install gloss --flags="-GLUT GLFW") allows you to start and stop the visualization from an interactive GHCi session (see the FAQ on the gloss page for further instructions). ?

Categories: Programming Tags:

Haskell Genetic Algorithm “Hello, world!”

April 15th, 2011 2 comments

A post on reddit linked to several implementations of a cute "Hello, world!" program demonstrating a genetic algorithm that evolves towards a target string. Example programs were written in several languages, and I thought a Haskell version could be worthwhile as it demonstrates the use of random numbers, an issue that frustrates many newcomers to the language.

I use the MonadRandom package to deal with plumbing a random number generator through the code, but otherwise stick to lists and the usual cast of characters from the Prelude.

> {-# LANGUAGE NamedFieldPuns, TupleSections #-}
> import Control.Applicative
> import Control.Arrow (second)
> import Control.Monad (liftM, replicateM)
> import Control.Monad.Random
> import Data.Function (on)
> import Data.List (minimumBy, sortBy, nub, (\\))
> import Data.Ord (comparing)
> import Text.Printf (printf)

The example programs associated with the original post all come with tests. Some of these are standard unit tests that verify hand-picked examples, while others verify properties of random processes. To support both kinds of tests, we use HUnit for unit tests, and QuickCheck for property testing.

> import Test.QuickCheck
> import Test.QuickCheck.Monadic
> import Test.HUnit hiding (assert)

For clarity, we create a type alias for the raw data we will be working with.

> type Gene = String
> target :: Gene
> target = "Hello, world!"

Next are the functions for working with individual Genes.

> mate :: RandomGen g => Gene -> Gene -> Rand g Gene
> mate g1 g2 = (++) <$> flip take g1 <*> flip drop g2 <$> pivot
> where pivot = getRandomR (0, length g1 - 1)
> mutate :: RandomGen g => Gene -> Rand g Gene
> mutate g = (uncurry (++) .) . second . (:) <$> delta <*> parts
> where delta = getRandomR (' ', 'z')
> idx = getRandomR (0, length g - 1)
> parts = second tail . flip splitAt g <$> idx
> fitness :: Gene -> Int
> fitness = sum . map abs . zipWith ((-) `on` fromEnum) target

We also want a random Gene generator.

> randomGene :: RandomGen g => Rand g Gene
> randomGene = replicateM (length target) $ getRandomR (' ', 'z')

Finally we can start thinking about populations. The example Python program I looked at had several parameters for initializing the population and controlling how it evolves. We capture those parameters in a record.

> data PopInfo = PopInfo { size      :: Int
> , crossover :: Float
> , elitism :: Float
> , mutation :: Float }

A Population is a pair of a record describing the population and a collection of Genes.

> type Population = (PopInfo, [Gene])
> defaultPop :: PopInfo
> defaultPop = PopInfo 1024 0.8 0.1 0.03

We will use a helper function to produce a randomized initial population.

> randomPop :: RandomGen g => PopInfo -> Rand g Population
> randomPop = liftA2 (,) <$> pure <*> flip replicateM randomGene . size

A tournament selection method is used to select parent Genes from a Population.

> tournamentSize :: Int
> tournamentSize = 3
> tournamentSelection :: RandomGen g => Population -> Rand g Gene
> tournamentSelection (info, genes) =
> minimumBy (comparing fitness) . map (genes !!) <$>
> replicateM tournamentSize (getRandomR (0, size info - 1))
> twoM :: Monad m => m a -> m (a, a)
> twoM = liftM (\[x,y] -> (x,y)) . replicateM 2
> selectParents :: RandomGen g => Population -> Rand g (Gene, Gene)
> selectParents = twoM . tournamentSelection

The meat of the algorithm is the evolve function that, appropriately enough, evolves a Population.

> evolve :: RandomGen g => Population -> Rand g Population
> evolve p@(info@(PopInfo {size, crossover, elitism, mutation}), genes) =
> (info,) . sortBy (comparing fitness) . (take idx genes ++) <$>
> replicateM (size - idx) (twoM getRandom >>= go)
> where idx = round (fromIntegral size * elitism)
> go (r1,r2) | r1 <= crossover =
> selectParents p >>= uncurry mate >>= addChild r2
> | otherwise = addMutation r2
> addChild r c
> | r <= mutation = mutate c
> | otherwise = return c
> addMutation r
> | r <= mutation = mutate . (genes !!) =<< getRandomR (idx, size - 1)
> | otherwise = (genes !!) <$> getRandomR (idx, size - 1)

Now we’re ready to kick off a multi-generation quest for the string we already wrote down up top.

> iterateUntil :: Monad m => (a -> Bool) -> (a -> m a) -> a -> m a
> iterateUntil stop f = go
> where go x | stop x = return x
> | otherwise = f x >>= go
> maxGenerations :: Int
> maxGenerations = 16384
> main = evalRandIO (randomPop defaultPop >>= iterateUntil done step . (, 0))
> >>= result
> where step (p,gen) = (,) <$> evolve p <*> pure (gen+1)
> done ((_, g:_), generation) =
> generation == maxGenerations || fitness g == 0
> result ((_, g:_), generation)
> | generation == maxGenerations =
> putStrLn "Maximum generations reached without success."
> | fitness g == 0 = printf "Reached target (%d): %s\n" generation g
> | otherwise = putStrLn "Evolution is hard. Let's go shopping."

In keeping with the spirit of the original programs, we include tests to exercise some of the component functionality we’ve put together.

> testGen = run (evalRandIO randomGene) >>= assert . check
> where check g = and $ map ($ g) [ (>= 0) . fitness
> , (== 13) . length
> , all (between 32 122 . fromEnum) ]
> between l r x = l <= x && x <= r
> testMut = run (evalRandIO $ randomGene >>= pairWithMutant) >>= assert . check
> where pairWithMutant = liftA2 (,) <$> pure <*> mutate
> check (g,m) = length g == length m && length (nub g \\ nub m) <= 1
> testMate = run (evalRandIO $ twoM randomGene >>= pairWithChild) >>=
> assert . check
> where pairWithChild (mom,dad) = (mom,dad,) <$> mate mom dad
> check (m,d,c) = length c == 13 &&
> (and . map (\(_,y,z) -> y == z) .
> dropWhile (\(x,y,_) -> x == y) $ zip3 m c d)
> unitTests = test [ "fitness1" ~: 0 ~=? fitness "Hello, world!"
> , "fitness2" ~: 399 ~=? fitness "H5p&J;!l<X\\7l"
> , "fitness3" ~: 297 ~=? fitness "Vc;fx#QRP8V\\$"
> , "fitness4" ~: 415 ~=? fitness "t\\O`E_Jx$n=NF" ]
> runTests = do mapM_ (quickCheck . monadicIO) [testGen, testMut, testMate]
> runTestTT unitTests

You can compile this post with GHC, or just run it in GHCi to poke around and ensure that the tests all pass.

Categories: Programming Tags: