This post was a originally published on Luma Labs, now dead.

As old as stimulus-response techniques are, they still form an important part of many AI systems, even if it is a thin layer underneath a sophisticated decision, planning, or learning system. In this tutorial I give some advice to their design and implementation, mostly out of experience gained from implementing the AI for some racing games.

A stimulus response agent (or a reactive agent) is an agent that takes inputs from its world through sensors, and then takes action based on those inputs through actuators. Between the stimulus and response, there is a processing unit that can be arbitrarily complex. An example of such an agent is one that controls a vehicle in a racing game: the agent “looks” at the road and nearby vehicles, and then decides how much to turn and break.

This definition is so broad that it is hard to think of agents that are not stimulus response agents. Some agents base their actions not on inputs, but on their internal state alone. See for example how random steering algorithms can lead to interesting behaviour.

Some agents don’t act. They just digest. The only way they can be useful if we look at what they are “thinking”. In fact, they are often designed with the sole purpose of “telling” us what they think, and they are not usually called agents, but rather decision systems. A face recognition system is an example of this. Of course, if we consider “telling” or “deciding” an action, these decision system become stimulus response agents too – but since we do not think of them as agents, the distinction is useful.

Design

1. Do not design your system around an overly simplistic view

The typical presentation of a stimulus response agent looks like this:

This presentation, while technically correct, is suitable only for flowers blowing in the wind. A real-world stimulus-response agent might look more like this:

Note that there are five factors that influence decisions:

• Game state is the mode of the game. For example, when an agent can be in a PLAY state (actively competing against the player) or a CHILL state (after taking over from the player at the end of the game.
• Global world state is the current conditions of the world, such as the weather or time of day.
• Local world state is what is normally thought of as the “stimulus”. It is the properties of the world around the agent that changes as the agents moves or otherwise acts in the world. For example, the shape of the path near a path following agent is part of the local world state.
• Agent parameters constitute the fixed limits of the agent’s sensory processing, and reactionary systems. For example, how far can the AI see, how quickly will it react, and how high can it jump. These parameters often describe the bounds for a state variable.
• Agent state is the changing parameters of the agent’s behaviour.

Note that all these can be mapped to the simplistic view of our agent. However, it is important to recognise the differences in these groups of variables, because the agent is made aware of them (connected to them) in different ways.

The world state is not sensed, but communicated through passed parameters, or a query to the world state object.

The game state is enforced by the game system, and typically the agent does not query or sense the game state – it is communicated from the outside.

The local world state is sensed by the agent through its sensors. Typically this means the agent must have hooks into the world. This data can be noisy or incomplete, so the agent might use filtering, and hence will store its own processed copy of the data.

The agent parameters of are partially hard-coded, and partially generated on construction, or read from a file. These do not change over the agent’s life time.

The agent state is calculated on updates. It has two parts:

• The discrete state is usually implemented as a state machine. In a vehicle driving agent, the discrete states can be PULL_AWAY, NORMAL, WRONG_WAY, and so on.
• The continuous state of an agent is the state determined by a set of parameters. For a vehicle driving agent, these can include “deviation_from_path”, “distance_from_closest_player”, “aggressiveness” and so on. Note that the continuous state for every discrete state can be described by a different set of parameters, possibly overlapping.In general, a different set of response equations will be used for every discrete state. In fact, the sole purpose of the state machine is to select the correct set of response equations for every situation.

The table below summarises the factors that influence a stimulus response agent’s decisions.

2. Clearly distinguish between objects that have state, and those that do not (or should not)

It is easy to build objects that have states when they should not. It is, for example, reasonable to design a path-following object so that it keeps track of an agent’s progress along the path (such as the last node reached), to simplify communication between the agent and the path. In this scenario, the communication looks like this:

There are at least two problems with this set-up:

• every agent must have its own copy of the path; and
• the agent can’t ask any path information that violates the path object’s assumption about where the agent is.

To illustrate the last point, consider the situation where the path keeps track of the last node reached by the agent, and updates that node as the agent reports its position to the path object. It uses this node and the next node to interpolate a suitable position for the agent to move to. All is fine, until the agent needs information about the future, for example, it’s projection one second from now on the path. One second from now, the agent may have reached another node, something this simple system cannot handle properly.

Another way to do this, is for the agent to store the last node reached and pass it along its position in queries:

This does not solve the problem really – must the agent now store a separate node for queries into the future?

The best solution is to have a path manager. It has two responsibilities:

• to store any agent state information relating to the path – the last node reached, in our example; and
• to handle calculations for queries (based only on data from the path).

Now, when the agent asks the path manager about its projection 1 second from now, the path manager can check against the path whether a node will be reached, and perform the correct calculation. [In practise, the agent will transform time into distance travelled in time at its current velocity, so that the path manager only needs to know about distance, not time].

Every agent must have its own path manager, but the path managers can share the same path.

Consider another, slightly different, example. A certain agent has a “angriness” factor, which is randomised every 10 updates between two fixed extremes. Every type of agent has a different set of extremes (a monster will have higher extremes than a cow, for example), which helps determine the agent’s “personality”. It is easy to dump the extremes and the current value in one object (typically, the agent object). This is a poor design for two reasons:

• There is unnecessary duplicated data, with the associated disadvantages (such as making network packets bulkier, making save files of the game state bigger, and so on).
• The distinction between of agent state and agent parameters gets muddy. This makes it easier to make mistakes, overcomplicate the design, and introduce unnecessary inefficiencies.

A better design is to put the extremes into an “angriness profile” object, which can be shared by agents of the same type; in the example there will be two instances, one for cows and one for monsters.

Note that this issue goes hand-in-hand with preferring composition over inheritance hierarchies.

3. Draw pictures to show the ideal response to a given stimulus situation

It is very important to define clearly what you want the agent to do given some stimulus. The following graphs shows a few examples for a vehicle steering agent.

Steering as a function of angular deviation from path.	Nitrous usage as a function of angular deviation from path (reckless driver).	Nitrous usage as a function of angular deviation from path (conservative driver).

Most variables will depend on more than one other variable. Do not worry too much about that when drawing these diagrams. Focus on typical cases, and jot down any requirements or assumptions.

Also, do not try to write down any equations for your graphs. Any relationship can always be represented by a response curve, something that is much quicker to set up and easy to implement. A response curve is essentially a look-up table with linear interpolation. It can be used to approximate arbitrary functions by specifying the input range, and a number of output samples. The idea is discussed more thoroughly in AI Programming Wisdom 1 (The beauty of Response Curves by Bob Alexander).

4. Draw System Diagrams to show how variables relate

In a complicated system, it is easy to lose track of how variables relate. A forgotten relationship can often lead to behaviour that is hard to understand. Having an easy-to-reference diagram will be a great help in debugging faulty behaviour.

The diagram above distinguishes between relationships under AI control, and relationships that are under physics control – i.e. Relationships that cannot be changed by the agent. The diagram also indicates whether a relationship is positive (i.e. more of the one will result in more of the other), or negative (i.e. more of one will result in less of the other).

5. Use a spreadsheet to graph response curves to typical situations

Once you have your equations, you have to make sure that they actually describe the curves you intend. Also, in complicated formulas, it is important to check how these curves behave on typical data sets – especially if you employ filtering (see below).

6. Eliminate unimportant feedback loops.

Feedback is important for proper control in most systems. However, every loop makes it more difficult to analyse and understand the system – something that becomes apparent when things go wrong. If your system is complex, try to remove some of the loops (by removing some of the dependencies). You can always complicate your system once you have it working.

7. Prototype your agent on a computer algebra system

A computer algebra system (CAS) is the perfect environment to prototype a stimulus response agent:

• The high-level syntax, dynamic variables, and fat libraries allows you develop designs quickly, so that you can try out more designs.
• The built-in graph plotting capabilities are useful for debugging (see below) and understanding complex designs.
• You can simulate scenarios that are hard to obtain in the actual game, leading to a design that is much more robust.

Filtering

8. Filtering input, output, and intermediate values can drastically improve behaviour

There are several reasons to use filtering in your agent:

• To reduce noise that might confuse your agent about the state the world (or it) is in.
• To smooth out (and hence fill in) incomplete data. For example, a path made of line segments can be “filtered” so that it appears smooth to the agent.
• To make the actions of the agent seem more determined. An agent that changes it’s mind 10 times a second is not very believable.
• To prevent abrupt reactions to changes in the world that will make the agent less believable.

There are several approaches to filtering. One approach that works well in many cases is to change a variable only by some maximum amount, like this:

updateFilteredX(x)
{
   filteredX += clamp(x – filteredX, -maxIncX, maxIncX);
}

Here maxIncX is much smaller than typical values of x (using a tenth is a good starting point).

When you have many filtered variables, it becomes unmanageable to filter them all in this way – for every variable you want to filter, you need three extra variables, and you need to remember to use the right updating technique whenever the variable changes. It is better to encapsulate the filtering in a special class (as has been done in the Special Numbers Library – see http://specialnumbers.google.com).

Other approaches to filtering include:

• statistical filtering techniques, such as least squares and Kalman filtering;
• convolutional filters; and
• PID controllers.

9. Make your filtering technique frame-rate independent, and do it from the start

The rate of updates of filters can have a dramatic effect on resulting behaviour. This rate can either fluctuate naturally (as it goes with frame rate), or be the difference between a Debug and Release versions of your application.

The resulting changes in behaviour makes it hard to:

• debug faulty behaviour; and
• spot problems to begin with (because it always works fine on your machine!).

How you make your filtering frame-rate independent depends on the method of filtering used. For the system explained above, it merely means updating with an amount proportional to the elapsed time (between updates).

updateFilteredX(x, elapsedTime)
{
   filteredX += clamp(x – filteredX,
                -maxIncX*elapsedTime, maxIncX*elapsedTime);
}

The updateFilteredX function above will use a smaller value for maxIncX than the one that does not compensate for time. It is for this reason that you should implement frame-rate independence from the start. Otherwise you will need to re-tweak all your values increment values again.

10. Do not filter the channels of n-dimensional data separately

If you filter a position, for example, by the method described above, you might be tempted to do something like the following:

updateFilteredXYZ(point)
{
   filteredPoint.x += clamp(point.x – filteredPoint.x,
                               -maxInc, maxInc);
   filteredPoint.y += clamp(point.y – filteredPoint.y,
                               -maxInc, maxInc);
   filteredPoint.z += clamp(point.z – filteredPoint.z,
                               -maxInc, maxInc);
}

The above code will lead to strange behaviour. Consider for example,

• filteredPoint = (0, 0, 0),
• point = (2, 1, 0) (for the next 4 updates)
• maxInc = 0.5

The variable filteredPoint will now go through these values on the next four updates:

(0.5, 0.5, 0)

(1.0, 1.0, 0)

(1.5, 1.0, 0)

(2.0, 1.0, 0)

If this filteredPoint represents an agent trying to reach (2, 1, 0)’s heading, the agent will first move to point (1, 1, 0), and then turn towards (2, 1, 0), and move to it.

The correct way to implement filtering for multi-dimensional data is to add a clamped amount of difference between the value and the filteredValue to the filtered value.

updateFilteredXYZ(point)
{
   updateAmount = clamp((filteredPoint – point).len(),
                        0, maxInc)
   filteredPoint += (point – filteredPoint) * updateAmount;
}

Now the next four values of filteredPoint lie in a straight line:

(0.316, 0.632, 0)

(0.632, 1.265, 0)

(0.948, 1.897, 0)

(1, 2, 0)

11. Beware filtering that reduces the range of a variable

Some kinds of filtering can reduce the range of a variable – for example, instead of reaching zero, a variable may merely come close to zero. This is important in checks that assumes variables will come to have certain value eventually – these checks must be modified to compensate for the range discrepancies.

12. Beware of accumulated latency when filtering variables through various layers of logic

Filtering usually results in some latency. Where filtering is used in a layered system, these latencies can add up, and produce a system that is slow to respond. In many situations these latencies can be reduced or eliminated if you can make accurate predictions about the future.

You can reduce the filtering by using it only where it is really important – not everywhere. You can also exchange smoothness for a quicker response. How you accomplish this depends on the method of filtering you use; in the filtering described above you can increase the value of the increment.

If you feel adventurous, you can try out a system that filters only after it becomes aware that it is necessary. For example, the amount of filtering you apply can depend on the amount of noise measured over the last number of frames. In many situations you can approximate the amount of noise by the amount energy of the measured signal. Systems like these can be hard to make robust without resorting to difficult mathematics or special software.

You can also switch to a more sophisticated filtering approach, such as a Kalman filter or PID controller. PID controllers uses prediction to reduce latency; Kalman filtering uses a system model that makes it more effective in complicated scenarios.

13. Be Careful When filtering combined signals

You might want to add (or combine otherwise) several signals, filter this combination, and work on this filtered combined signal. This is a mistake when the signals have typical data transitions at different frequencies: ones signals noise is another’s data.

To illustrate the problem, look at the two signals below. Signal 1 is a low frequency signal, with a spike of noise. Signal 2 is a pure high frequency signal. The noise is clearly visible in the combined signal. Also shown are several filtered versions of the combined signal. Notice that the more we filter, the more we reduce the high frequency signal. And if we do not filter enough, the noise is still very present.

The solution is to filter the signals separately – the low frequency signal can be filtered much more than the high frequency signal – and only then combine the filtered sequences. The result is shown in the last figure.

Signal 1	Signal 2	Combined signal
Combined signal filtered (0.1)	Combined signal filtered (0.5)	Combined signal filtered (1.0)
Signals filtered before combined.

14. Avoid if-then logic in response calculations.

If-then logic signifies one of two things:

• a state difference that should be implemented as proper discrete state of the state machine; or
• the need for a step function.

In the first case, your state should be broken up into two states, with proper transitions. A stimulus response agent is complex enough with you having to track trough various conditions to understand behaviour. Every state should have one response that depends continuously on the stimulus. Because state machines scale so badly, introducing more states is not very desirable either. There are of course several ways to organise states properly to alleviate the scaling problem, but this is the topic of another tutorial.

If your if-then logic is against the value of a variable (for example ‘if x > 0′), you can replace it with an appropriately chosen step function. Using step functions instead of if-then logic simplifies the design, and can be replaced with drop-in “smooth” versions when it becomes evident that it is necessary.

Testing and Debugging

15. Rendering the current state-machine state on screen.

It is a small point, but often strange behaviour results not from incorrect response to stimulus, but from incorrect state transitions. When implementing your agent inside the game, having state transitions visible on screen will make that kind of bug that much easier to spot.

16. Graph variables and their internal states on screen.

Graphing variables (filtered and unfiltered) over time on screen will save you hours of debugging time. Whenever your agent does something funny, you can see at a glance if there are variables with unexpected values that could explain the agent’s actions, perhaps using the CAS prototype as a reference.

Unfortunately, this tip becomes less useful as the system gets more complex. Having twenty plots on your screen will do little to spot problems.

17. Use unit tests to ensure correct implementation of formulas.

Stimulus response agents are perfect candidates for unit testing, especially when using the CAS prototype as a reference implementation. Unit testing your assumptions about the system will ensure correctness, and will help to make you aware of changes that break the system early. It will also make it easier to change and optimise formulas, or change them with smoother or sharper versions. Having a constant reference against which you can tweak parameters also allows you to be daring and make aggressive changes.

18. Use generators to mock the agent’s environment.

A generator is a object that generates a sequence of values. Typically it has a getValue() method, which returns a different value every time it is called. Generators are very similar to iterators, except that they usually generate values on the fly, rather than pointing to members of a container (although reading values from a container is one valid way to generate values). Generators are often used to implement iterators.

To implement unit tests and CAS prototype, you need to mock various aspects of the game, the game world, and other agents. Using generators for this purpose makes mocking unbelievably easy and quick. Use generators to generate streams of

• elapsed times;
• path angle deviations;

• positions and velocities of other agents; and
• state transitions.

You can even mock parts of your agent to test sections in isolation. It helps if you build a library of generators and functions (response curves) that can be stacked together. You will do good to implement at least the following functions:

• step, line, and ramp;
• clamped line and sigmoid (or atan);
• cos, square, and saw tooth;
• exponential and exponential decay;

You should implement the following generators:

• a general generator that generates values from a function;
• a constant generator;
• a sum and product generator (generates the sum or product of two generators); and
• generators that generate random values of various distributions.

Conclusion

A complex stimulus response agent can be a challenge to design and implement properly; every approach has its set of headaches, and every solution brings other considerations to the table. Solving problems is just half the fun – not solving them is the other half!

In the Dev.Mag article Guerrilla Tool Development, I give some ideas for getting some useful tools on a tight budget. Check it out!

(Original image by Hljod.Huskona / CC BY-SA 2.0).

I used to hate neural nets. Mostly, I realise now, because I struggled to implement them correctly. Texts explaining the working of neural nets focus heavily on the mathematical mechanics, and this is good for theoretical understanding and correct usage. However, this approach is terrible for the poor implementer, neglecting many of the details that concern him or her.

This tutorial is an implementation guide. It is not an explanation of how or why neural nets work, or when they should or should not be used. This tutorial will tell you step by step how to implement a very basic neural network. It comes with a simple example problem, and I include several results that you can compare with those that you find.

I tried to make the design as straightforward as possible. The training algorithm is simple backpropagation. There are no hidden layers (I will treat that in an upcoming tutorial), no momentum, no adaptive learning rates, and no sophisticated stopping conditions. Those are, in a sense, easy to add once you have a working neural net against which you can benchmark more elaborate designs.

To keep the implementation simple, I did not bother with optimisation. This too can easily be addressed once a working neural net is in place against which you can verify correctness and measure performance improvements.

The brief introduction below is a very superficial explanation of a neural net; it is included mostly to establish terminology and help you map it to the concepts that are explained in more detail in other texts.

Preliminary remarks and overview

What we are doing

The problem we are trying to solve is this: we have some measurements (features of an object), and we have a good idea that these features might tell us in which class the object belongs. For example, if we are dealing with fruit, knowing the size, colour, and “roughness” of the skin, we might deduce which type the fruit is. Of course, we want a program for this.

Now in heaven, there is a function for exactly that task – we give the function the features, and it spits out the class.

Down here on earth, we are not so lucky (well, not for most problems anyway). So instead, we have to do with some approximation. A neural net is one possibility – there are also others. A neural net (one without any hidden layers) is parameterised by a weight matrix. Different problems in general have different weight matrices. To solve our problem, we need to find a suitable matrix.

If both our set of known samples and the problem itself are reasonable, we might expect to find such a matrix. But not directly – we have to use some kind of iterative scheme – we have to “train” our neural net. We use some of the known samples for this.

Behind the scenes, the weight matrix and the feature vector are combined using some matrix operations to give the output vector, which is converted to a class. The mathematical details of this can be found elsewhere (more).

After we have trained the neural net, we can use it in an application to classify objects that we may or may not have encountered before. Generally, the training and application programs are separate. The training program is the difficult part; most of the tutorial deals with training.

It is common to divide the known samples in three sets, called the training set, validation set, and test set, typically divided in the ratio 50%, 25%, and 25%. The training set is used to update the weights iteratively; the validation set is used to stop the training algorithm, and the test set is used to estimate how well our trained neural net will do in the wild.

Each iteration of the training algorithm is called, quite poetically, an epoch. Many different things affect the number of epochs we will need to train the neural net: the problem complexity, how we scale our updates, and how the weight matrix has been initialised.

We can control the speed with a parameter called the learning rate. In general, higher learning rate means faster learning (although, when it is set too high, the network might become unstable and not learn at all).

Representation

We usually represent the features (our measurements) as a vector of numbers, like this:

[0.01, 1.07, 3.60]

In cases of non-numeric values, we convert whatever we have to numbers with some scheme. For instance, when dealing with colours, we might want to use RGB values, or we might want to make a list of colours (red, yellow, green, etc.), assign a number to each of these labels, and use the corresponding number.

Classes, usually not being numbers, are similarly treated: we assign a number (often arbitrarily) to each class, and use that in our computations. For the actual training, however, we use an output vector. The output vector has a 1 in the position of the class number, and 0 everywhere else.

Thus, if we are dealing with three classes of fruit (apple, pear, banana), and number the classes 1, 2, and 3 respectively, we will have the following output vectors:

It is convenient to put all the inputs of a set together in a single matrix, where each row is a sample. Similarly, outputs and classes are also put into matrices, with input sample in a row (say row number n) corresponds to output sample in row n, and also the class in row n.

The implementation below makes use of high-level matrix operations. This avoids many of the errors that can creep in loop-dense code. The following table will give dimensions of all the matrixes involved; this will be helpful during implementation (especially for assert statements). The numbers are given in terms of the number of inputs (features) of the problem, the number of outputs.

Overview

The basic idea of the training algorithm is the following:

First we load in the data: the training samples, the validation samples, and the test samples.

Then we start to train: we run the backpropagation algorithm on random samples. After each iteration, we see how our network is doing so far (on the validation set), and then we decide whether to keep training or not.

After we stopped, we do a final evaluation of our network on the test set – this gives us an indication of whether the neural net will generalise well to samples not originally in the training set.

After we have found a weight matrix that we can live with, we cfan incorporate this in the application were we need the functionality.

The training program has x functions:

• train
• load data
• feed-forward
• evaluate network
• backpropagation
• output to class
• class to output
• activation
• activation derivative

The application program uses these functions:

• feed-forward
• output vector to class

Implementation

1. Gather the necessary libraries (or write them)

You will need the following libraries:

• A library that supports matrix algebra; and
• A library that plots graphs (x versus y).

If you can’t find a matrix library for your implementation language, then you can write a simple library yourself. Since neural nets do not require matrix inverses or long chains of matrix products, you need not worry (much) about numerical stability, thus the implementation is straightforward.

The implementation needs to support the following matrix operations:

• matrix transposition;
• matrix addition;
• matrix multiplication with a scalar;
• ordinary matrix multiplication;
• Hadamard multiplication (component-wise multiplication);
• Kronecker multiplication (only necessary for between row and column vectors); and
• horizontal matrix concatenation.

The first few operations are standard matrix operations, but you might be less familiar with the last three. (Also check out the Wikipedia article on matrix multiplication – it covers all the types of multiplication mentioned here.)

Hadamard multiplication of matrices is defined for two matrices of equal dimensions. Each component of the new matrix is the product of corresponding components in the two multiplicands, that is:

Z[i][j] = X[i][j] * Y[i][j]

The Kronecker product of a row vector and column vector is defined as a matrix whose components are given by:

Z[i][j] = X[0][i] * Y[j][0]

It is possible to define the product for arbitrary matrices, but we don’t need it.

The horizontal concatenation combines two matrices with the same number of rows. For example, the matrices A and B below will be concatenated to form the new matrix C:

A simple implementation simply constructs a new matrix whose components are given by

if j < X_width
  Z[i][j] = X[i][j]
else
   Z[i][j] = Y[i, j – X_width]

If no graph libraries are available, simply write a function that will output a tab-separated list of the input and output sequences to plot. You can then load or paste this into your favourite spreadsheet program to make the necessary plots.

2. Implement Output and Class conversion functions

This is very simple: implement a function that converts an output matrix to a class number vector, and another that converts a class number to an output vector.

For example, the output_to_class function will take the following matrix

1 0 0

0 1 0

0 0 1

1 0 0

0 0 1

and convert it to:

1

2

3

1

3

(The second function will convert the second matrix back to the first matrix).

3. Implement a function to read in data files

For this tutorial you can use the following three files:

• iris_training.dat
• iris_validation.dat
• iris_test.dat

These three files contain samples from the ICU iris dataset, a simple and quite famous dataset. In each file, samples or contained in rows. Each row has seven entries, separated by tabs. The first four entries are features of irises (sepal length, sepal width, petal length, and petal width); the last three is the outputs denoting the species of iris (setosa, versicolor, and virginica). I have preprocessed the values a bit to get them in the appropriate ranges.

You must read in the data so that you can treat the inputs of each set as a single matrix; similarly for the outputs.

I find it useful to store all the data in a structure, like this:

• data_set
  • input_count
  • output_count
  • training_set
    • inputs
    • outputs
    • classes
    • count
    • bias
  • validation_set
    • inputs
    • outputs
    • classes
    • count
    • bias
  • test_set
    • inputs
    • outputs
    • classes
    • count
    • bias

This makes it more useful to send the data as parameters.

4. Implement an activation function and its derivative

The activation function must take in a matrix X, and return a matrix Y. Y is computed by applying a function component-wise to X. For now, use the hyperbolic tangent function:

5. Implement the feed-forward function

The function must take as arguments an input matrix, weight matrix, and a bias node matrix.

The function should return an output matrix, and a net matrix

These are computed as follows:

net = mul(weights, horcat(inputs, bias))
output = activate(net)

The bias matrix is a constant column vector of 1s with as many rows as the input matrix. This vector corresponds to the bias nodes. The implementation here is a bit clumsy, but for now, the approach used here minimises the potential for error.

6. Implement a weight initialisation function

This function must take in a maximum weight, a width and height, and return a matrix of the given width and height, randomly initialised in the range [-max_weight max_weight].

7. Implement a function that evaluates the network error.

The function must take in:

• an input matrix,
• a weight matrix,
• a target output matrix,
• a target class matrix,
• a bias matrix.

The function must return the error e, and the classification error c.

To compute these, first compute the output matrix Z using the feed-forward function (you can ignore the net matrix).

[output net] = feedforward(inputs, weights, bias)

Now subtract the target output matrix from the output matrix, square the components, add together, and normalise:

error = sum_all_components((target_outputs – outputs)^2) ...

   / (sample_count * output_count)

From the output matrix, calculate the classes:

classes = classes_from_output_vectors(outputs)

Count the number of classes that corresponds with the target classes, and divide by the number of samples to normalise:

c = sum_all_components(classes != target_classes)/sample_count

(Here, our inequality returns a matrix of 0s and 1s, with 1s in positions where the corresponding components in classes and target_classes are not equal.)

8. Implement a dummy backpropagation function

The function should take in:

• An input matrix
• A weight matrix
• a learning rate (eta, as in the Greek letter)
• a bias vector

The function must return an updated weight matrix. For now, return W as is.

9. Implement the train function

The training function should take in three sets, the training_set, validation_set, and test_set. Implement a way to limit the maximum number of samples that will actually be used for training (you can also do his in the main program described in the next section). This is very helpful for debugging purposes (especially if you plan to later replace the backpropagation algorithm with something a little faster – and more complicated).

The function should return a weight matrix, and error values as floats.

Initialise a value plot_graphs to true. This is a debug flag, so it is appropriate to implement this as a macro if it is supported by the implementation language.

The function should initialise a weight matrix using initialise weights. For now, use a max_weight of 1/2.

The function should also construct three bias vectors bias_training, bias_validate, and bias_test. Each must contain only 1s, with as many rows as there are inputs in the training, validation and test sets respectively.

Implement a while loop that stops after 500 iterations. (We will change the while condition later to something else, so do not use a for loop).

Inside the loop, call the backpropagation algorithm. Use the training set inputs, the weights, (for now) a fixed learning rate of 0.1, and bias vector bias_train. Assign the result to weights.

Still inside the loop, call the network error function three times: one time for each of the training, validation, and test sets. Use the weight matrix, and the appropriate bias vector. Wrap these calls in an if-statement that tests for a value plot_graphs. (If your language supports it, you can use conditional compilation on the value of plot_graphs).

Store the errors in six arrays (error_train, classification_error_train, etc.), with the current epoch number as index.

After the loop, plot the six error arrays as a function of epoch number. Wrap this in an if-statement (or conditional compilation statement) that tests for the value plot_graphs.

Call the network error function again, on all three sets as before.

Return the weights, and the six errors.

10. Implement the main training program

The program should load in the sets (using the load_sets function), and pass these to the training algorithm.

11. Run the program

The important thing is that everything should run. You should see your error plots; at this stage they should be straight, horizontal lines. Because of the random weight initialisation, we cannot predict where these lines will lie (so do not be alarmed if they do not look exactly the same as below – as long as they are straight and horizontal).

12. Implement the backpropagation function

You have already created the dummy function; now you can put in the actual calculations.

First, select a random sample.

Now, calculate the net matrix and output matrix using the feed-forward function.

[output, net] = feedforward(random_sample, weights, bias)

Calculate the error vector

error_vector = target_outputs - outputs

Calculate the sensitivity.

delta = hammard(error_vector, activation_diff(net))

The corresponding mathematical expression in the textbook might look like this:

Calculate the weights delta:

weights_delta = scalar_mul(eta, kronecker(transpose(outputs), delta))

The corresponding mathematical expression in the textbook might look like this:

Update the weights:

weights = add(weights, weights_delta)

and return the matrix.

13 Run the program (AGAIN)

First, set the debug option to train on only one sample. The curve should like like this:

Notice that the error curves are smooth, the training error rate goes to 0, and the other error rates go to 0.4, plus or minus – this depends on the initial weights. It will also be different if you use a sample other than the first.

If your curve looks right, change the debug option to train on all samples. Your error curves should now look something like this:

The important feature of the error curves is that they should steadily descent on average. If your curve does not resemble the one shown here; there is a mistake in your implementation. It is often helpful to limit the training set to a single sample (thereby eliminating the randomisation of sampling during training) to see what is going on.

Although the errors often follow the relationship training < validation < test, this is not always the case, so do not be alarmed if this is not true on a run (You can see it is not always true in the run above).

14 Implement a proper stopping condition

Change the while loop to stop when the validation error drops below a threshold. Note that this threshold usually depends on the problem. There are better stopping conditions that are less sensitive to the problem at hand, but this one will do for now.

15 Implement a statistical analysis

This part is important for you to get an idea of the robustness of the neural net. In practice, a very simple analysis will suffice.

This part need to train the algorithm 30 times, and then report the mean, standard deviation and maximum of the

• training time,
• regression error, and
• classification error (on the test sets).

In general, you would like all these values to be “low”. Here are some experiments for the iris data set with different learning rates. For each, 30 runs were made; other parameters are as described earlier (max_weight = 1/2, validation_stop_threshold = 0.1).

Using the neural net in a program

To use the neural net in a program after you have trained it, you need to save the weights found by the training program to a file. Your application must then read in this weights, and use it with the feedforward function to calculate the class.

Here is pseudo-code for the program. Note that the training program and classification program needn’t be implemented in the same language. This allows you to take advantage of speed and interface components that might not be available in your target platform. It is a very good idea to implement your training algorithm on a computer algebra system (such as Matlab or Octave) where you can take advantage of both matrix and graphing capabilities (the code provided below works in both).

Using the training algorithm on other problems

When using your neural net for other algorithms, you might need to change the learning rate, stopping threshold, and weight_max for weight the initialisation. The error plots are indispensible for this purpose.

As long as the learning rate is not too high, it should not affect the quality of the solution, only the number of iterations necessary to obtain it. The stopping threshold, however, has does affect the quality of the solution: if it is tool low, the problem might not be solved, or the neural net might train very long; if it is too high, you will get poor performance. There are better stopping conditions available; once you have everything working, you should investigate these.

Remember that you should not base your decisions on a single run, as runs can differ quite drastically from one another. Perform a few runs, and base decisions on these.

The value weight_max can usually be chosen as 1/sqrt(input_count) with good results, but here too it depends on the problem.

For further benchmarking, check out the Proben1 datasets, described here. (Also check out the erratum note on this site – search for proben on the page.)

Download

The following code works in Matlab and Octave. (Included is a randint function; if you are using Matlab you can remove it, because it is already implemented in Matlab).

iris_data_files.zip (3 KB)

basic_neural_net_0_1.zip (10 KB)

The program provided cheats a little – instead of reading in the raw .dat files, it reads in the .mat file that already has the data in the right structure.

Check it out!

The Python Image Code has also been updated with some of the algorithms explained in the article.

Download

XSIModToolLevelEditor3.zip (5.2 MB).

I finally laid my hands on Donald Knuth’s The Art of Computer Programming (what a wonderful set of books!), and found a neat algorithm for generating random integers 0, 1, 2, … , n – 1, with probabilities p_0, p_1, … , p_(n-1).

I have written about generating random numbers (floats) with arbitrary distributions for one dimension and higher dimensions, and indeed that method can be adapted for generating integers with specific probabilities. However, the method described below is much more concise, and efficient (I would guess) for this special case. Moreover, it is also easy to adapt it to generate floats for continuous distributions.

Description of the Algorithm

The basic idea of the algorithm is simple. We have two tables of length n that contains integers (K, L), and a third table that contains probabilities (P). The first table merely contains the integers 0 to n-1 (thus, we need not actually store it explicitly). The other two tables are computed before generation (more about that below).

To generate a random number, we generate a random integer (uniformly distributed between 0 and n-1 inclusive), and a random float (uniformly distributed between 0 and 1). The integer tells us which cells in the table to use. First we lookup the probability in P at that index. If the float is smaller than that value, we return the integer in K, otherwise we return the integer in L. (Note, the integer in K is exactly the random integer itself. Therefore, we do not actually have a table K – we simply use the random integer.)

Now this will only give the desired result if the tables have been constructed for this to work. Before looking at how these tables are generated, let us look at a very simple example. Suppose we want to generate 0, 1, 2, with probabilities 3/18, 7/18, 8/18. Can you see why the following tables will work?

There is only one way to generate 0: if the random integer is 0, and the random float is below 1/2. This will happen with probability 1/3 x 1/2 = 1/6 = 3/18.

There are two ways of generating 1:

• if the random integer is 1 (probability 1/3 = 6/18); or
• if the random integer is 2, and the float is above 5/6 (probability 1/3 x 1/6 = 1/18).

Adding these probabilities, we get 6/18 + 1/18 = 7/18.

There are also two ways to generate 2:

• if the random integer is 0 and the random float is above 1/2 (1/3 x 1/2 = 1/6 = 3/18); or
• if the random integer is 2 and the random float is below 5/6 (1/3 x 5/6 = 5/18).

Adding these probabilities, we get 3/18 + 5/18 = 8/18.

Generating the tables

Consider this problem:

• We have n squares that we want to paint.
• We have five colours of paint, possibly different amounts of paint for each colour.
• Each square has a border in one of the colours; no two borders are the same colour.
• In total, there is just enough paint to cover the n squares exactly.
• We want to paint each square with at most two colours, with one colour matching the border.

To do this, we sort the paint buckets in ascending order. We paint the square that matches the first bucket with the first bucket, and whatever remains with the last bucket. The first bucket is empty (why?), and there might be some paint remaining in the last bucket. We now put this bucket back so that the buckets are sorted according to the new quantities of paint. The first square is completely covered (why?). Note that the painted square’s colour corresponds with the depleted colour.

The situation is now: we have n-1 unpainted squares, and n-1 colours of paint. This is the same problem as the initial problem, with one less square and one less colour. Therefore, we proceed as before, and repeat this process until all the squares have been painted and all the paint has been used.

To answer the two why’s above:

The first bucket is always empty, because the smallest bucket cannot cover more than one square. Since it is the smallest bucket, all other buckets must contain at least as much paint. Thus, we have n colours, and enough paint of each colour to paint more than one square. Thus, in total, we must have more paint than is required to paint n squares. But we said that we have exactly the right amount of paint needed, not more. Therefore, the smallest amount of paint cannot cover more than one square.

The first square is always covered completely, since the last bucket always contains enough paint for at least one square. If it did not, since it is the largest bucket, all other buckets will paint less than one square. In total, we would have n colours, each that can cover less than one square. Thus, all our paint will cover less than n squares: we do not have enough paint. But we said that we do, so this cannot be. Therefore, we must have enough paint in the last bucket to cover at least one square.

Below is an illustration of this process for three colours. Here we assume that 1 litre of paint covers 1 square. We have 1/3 = 3/6 litre of red, 7/6 litre of green, and 8/6 litre of cyan.

Sort Paint

Paint from first bucket

Paint from last bucket

Sort Paint

Paint from first bucket

Paint from last bucket

Sort paint

Paint from first bucket

As you can see, the solution above corresponds with the tables given in the example above. Indeed, the algorithm for calculating the tables is exactly the same as the paint algorithm:

• The n colours correspond to the n integers (from 0 to n-1) that we want to generate.
• The initial amounts of paints corresponds with the (relative) probabilities that we want to generate each integer.
• The amount of paint used to paint a square of the same border is the entry in table P – the probability of using the number associated with that cell (i.e., “the border”).
• The other colour used to paint a square (if any) corresponds to the entry in table K.

A Small Optimisation in the Implementation

Generating two uniform random numbers for the price of one

There is a trick to generate a random integer (0 <= n < k) and a random float (0 <= x < 1) from a single random float (0 <= u < 1) that is used in the implementation of this algorithm (see download below).

The trick is:

• n = floor (uk)
• x = uk – n

This assumes things about the random generator and the accuracy required. (I do not want to get into the details here).

Download

A python implementation of the above algorithm.

non_uniform_random_int.py

I will cover the case for 2D sets here, but the ideas are easily extended to any dimension.

Visual Inspection From a Scatter Plot

The easy way to estimate the distribution is to simply look at a scatter plot of the samples.

	A scatter plot of a 2D sample set. (A pixel was simply drawn at each point of the sample. The images was slightly blurred to give the pixels more substance).Here we can already estimate the distribution, as we can intuitively “see” the shape.

If the distribution is simple, we can often “guess” suitable parameters for a mathematical function. But this method is not suitable when more accuracy is required. And obviously this method does not work for higher dimensions.

I always use a visual inspection as a first step when using the other methods described here. First, it allows you to decide which approach to use. Second, it serves as a rough benchmark to test the results of a better method against.

Averaging Over a Grid

For this method, we divide the domain in cells, and count the number of points in each cell. To normalise, we divide the counts by the total number of points to get a distribution.

Typically, we choose a cell size large enough so that all cells contain a few points.

It is possible to use different cell sizes, so that regions of higher density has more cells, for instance, by using a quad tree. However, this approach makes it difficult to interpret the values (they should be scaled down by a factor of the area they represent), and to generate random numbers from this distribution. I can’t really imagine a situation where this approach will be preferred above using a regular grid or convolution.

Convolution

There are two ways to implement convolution. The first should be used when the sample set is small relative to the size of the domain, otherwise the second method will be more efficient.

Sparse Point Convolution

Generate a very granular empty grid over the domain.

For each point p in the sample set,

• map p to the grid (calculate the coordinates x,y of p in the grid),
• add a number to all the cells in the neighbourhood of p in the grid.

Now normalise the grid.

The neighbourhood is typically a circle or a square. The number you add can be the same for all neighbours (in which case any positive number will do), or it can be scaled depending on the distance from p.

The size of the neighbourhood depends on the density of your sample. In general, the denser it is, the smaller the neighbourhood can be. Smaller neighbourhoods lead to faster execution. To prevent “holes” in the approximation, use a radius that corresponds to the maximum of the distances between a point and its closest neighbour.

	Estimation with a constant square neighbourhood.
	Estimation with a constant circular neighbourhood.
	Estimation with a circular neighbourhood with a falloff.

Traditional Convolution

Generate a very granular empty grid over the domain.

For each point p in the sample set, map p to the grid (calculate the coordinates x,y of p in the grid), and add one to that location in the grid.

Now choose a square, symmetrical convolution matrix, and perform a discrete convolution on the grid:

new_grid[i, j] = sum_{i,j} grid[i][j] * c[i][j]

Here the i, j go over the indices of the convolution matrix. (Normally, the convolution is defined as new_grid[i, j] = sum_{i,j} grid[n - i][n - j] * c[i][j]. However, since we are suing a symmetrical matrix, these definitions are equivalent, and we need not perform the extra calculation).

Now normalise the grid.

The convolution can be a square or circle of 1s, or be filled with numbers that grow smaller outwards. These correspond to the three neighbourhoods described for the sparse convolution.

Note that the new_grid is larger than the original by one less than the size of the convolution matrix in each dimension. The centre of the new grid corresponds with the original grid.

For example, if the original grid was 100×100, and the convolution matrix was 5×5, the new grid will be 104×104. The point (2, 2) in the new grid corresponds to point (0, 0) in the original grid.

About Normalisation

It is customary to normalise the distribution so that all the probabilities add to 1. But for many purposes we need only relative probabilities, and this step can be skipped. For example, the method of generating random numbers described in a previous post uses only relative probabilities.

A Few Tips

• Always test your distribution by generating a random set from it, and comparing it with the original sample. They should match qualitatively.
• The smaller your sample set, the cruder the approximation should be. That is, cells, neighbourhoods or convolution matrices should be big. There is a limit to the accuracy you can obtain from any sample – if you try to exceed it, your results will be poor.
• The most common implementation errors are made at the borders (of any grid or matrix) – watch out for them!

Download

There is an example of implementation in 2D in with the Python Image Code. See the file random_distributions_demo.py.

I have already covered how to generate random numbers from arbitrary distributions in the one-dimensional case. Here we look at a generalisation of that method that works for higher dimensions.

The basic trick, while easy to understand, is hard to put in words (without reverting to mathematical equations). For two dimensions, we divide the plane into slices. Each slice is a 1D distribution. We also calculate a distribution from summing the frequencies in each slice. The latter distribution gives us one coordinate, and the appropriate slice to use. The distribution of that slice then gives the second coordinate. All distributions are put into inverse accumulative response curves as was done to generate one-dimensional random numbers. (You should review that before implementing the 2D case).

In more dimensions, we also slice the space up into 1D distributions. Sums of these give us more distributions, which we can sum again, and again, until we reach a single distribution. This is used for the first coordinate, and to determine which distribution to use for the next coordinate. This goes on, until a 1D slice gives us the final coordinate. Again, all distributions are converted to inverse accumulative response curves.

If the above is unclear, I hope the detailed description below clears things up.

Two Dimensions

Describe the Distribution

Divide the domain into a regular grid. Let us use the domain 10..40 x 10..40. This means, all the points we will generate will fall in the square between the lines x = 10, x = 40, y = 10 and y = 40.

For later use, we will need the minimum x-coordinate of the domain (x0 = 10), the width of the domain (w = x1 – x0 = 40 – 10 = 30).

Assign to each cell the number of random numbers that should be generated for that cell for any suitable total of points. You need not normalise these values.

For our example we will use the 4×4 grid shown below:

For later use, we need the width of the grid (gw = 4).

The method described here won’t work if any of these values are 0.

Calculate Cumulative Grids

Calculate the cumulative sums for the columns of the array.

Calculate a cumulative sum for the last row

Construct Inverse Response Curves

For each column of cumulative sums, append a zero at the beginning, and create an inverse response curve as described for the one dimensional case. Thus the first columns inverse response curve will be created from 0, 1, 3, 7, 15. Call these curves cy[0] … cy[3].

For the row of cumulative sums, append a zero at the beginning, and create an inverse response curve. Call this curve cx.

Determine Random Point

You are now ready to generate random points from the specified distribution.

First, generate two uniform random numbers, urx and ury.

Use urx to do a lookup into curve c[x]. The result rx gives the x-coordinate of your non-uniform random number. Use this number to decide which column curve to use: subtract the minimum domain x-coordinate, divide it by the domains width, multiply it by the number of columns, and floor it to an integer.

ix = (rx – x0)/w * gw

Now use ury to do a lookup in cy[ix]. The result ry is the y-coordinate of the number.

This image shows a plot of samples generated from the distribution specified above (normalised to fit in the image boundaries). Brighter areas indicate that more points occur in that region.

Visually, it looks like the sample mimics the original distribution well.

N Dimensions

Describe the Distribution

Divide the domain in a regular N-dimensional grid, and assign frequencies to cells in the grid. For future calculations, we will need the lowest coordinate of the domain along each dimension (x0_k), the width of the domain for each dimension (w_k), and the number of cells in the grid for each dimension (gw_k).

Calculate Cumulative Grids

Accumulate sums along one dimension in a N-dimensional grid G_N. The last subgrid contains the totals of all “columns” along that dimension, and is a (N-1) dimensional structure.

Accumulate sums of this subgrid into a (N-1) dimensional grid G_(N-1) along another dimension. Again, the last subgrid contains the totals of the columns of that dimension.

This must be repeated until a single row, G_1 is produced. In general, the last subgrid in G_k contains the totals of columns in G_(k-1) along dimension k. This must be accumulated in a (k-1)-dimensional grid G_(k-1).

Construct Inverse Response Curves

Now for each grid G_k, we need to construct inverse response curves from the columns along dimension k. Remember to append 0 at the beginning of each column. The curves must be put into a k-1-dimensional structure, C_k.

Determine Random Point

Generate N uniformly distributed random numbers ur_1…ur_N.

Use ur_1 to do a lookup in C_1. The result r_1 is the first coordinate of your point. Use it to determine to calculate an appropriate index i_1:

i_1 = (ur_1 – x0_1)/w_1 * gw_1

i_1 is the index of the curve to use from C_2: lookup ur_2 in C_2[i_1] to obtain r_2. This is the second coordinate of your point. Determine i_2:

i_2 = (ur_2 – x0_2)/w_2 * gw_2

i_1 and i_2 determine the curve to use from C3, i.e. the curve C3[i_1, i_2].

Repeat this process until all coordinates are determined. In general, use ur_k to do a lookup in C_k[i_1, i_2, …, i_(k-1)] to obtain r_k. Use this to calculate i_k:

i_k = (ur_k – x0_k)/w_k * gw_k

That’s it!

Download

There is an example of implementation in 2D in with the Python Image Code. See the file random_distributions_demo.py. (I know it is annoying that it is coupled with all the other image code… I am working on a better solution!)

A cellular automata system is one of the best demonstrations of emergence. If you do not know what cellular automata (CA) is, then you should go download Conway’s Game of Life immediately:

Conway’s Game of Life

Essentially, CA is a collection of state machines, updated in discrete time intervals. The next state of one of these depends on the current state as well as the states of neighbours. Usually, the state machines correspond to cells in a grid, and the neighbours of a cell are the cells connected to that cell. For a more detailed explanation, see the Wikipedia article.

Even simple update rules can lead to interesting behaviour: patterns that cannot be predicted from the rules except by running them. With suitable rules, CA can simulate many systems:

• Natural phenomena: weather, fire, plant growth, migration patterns, spread of disease.
• Socio-economic phenomena: urbanisation, segregation, construction and property development, traffic, spread of news.

In games, these can be used as game elements (such as fire and disease), or they can also add some decorative life to a game (such as city development in a simulation game).

The great thing about CA systems is that they are not hard to implement, and they have intrinsic fun value (I mean here, for the programmer!).

Before we get into the details, let me add a small disclaimer here:

• This tutorial is not comprehensive. I only touch on a few ideas in a wide field. See some of the resources at the end for more information.
• This tutorial is not representative. I merely chose topics I found interesting or thought would be useful for game developers.
• I am not suggesting that using cellular automata is the best way to simulate any of the examples given here. For any problem, there might be solutions that are faster, easier to implement, or yield better results.

Implementation

To implement a CA system, you need to do the following:

• Use a suitable system to store states, for example, a 2D array.
• Decide on a suitable neighbourhood shape, and work out some preliminary rules. Decide how you will handle edges.
• Implement a suitable way to describe your rules.
• Implement an update function that hooks into the main game loop.
• Implement a visualisation scheme (for example, drawing tiles depending on states in the array).

Some details of these are described below.

Grid

A typical implementation uses two grids for storing states. One is the current world state; the other is the new world state. Once all cells have been updated for the current iteration, pointers to the arrays can be swapped. This way updating the cells one by one will not have an effect on other cells until the next iteration, and you do not have to create a new array in each iteration.

Using grids with square cells are the easiest to implement, but the concept is also applicable to other grids: hexagonal and triangular grids are also common, and even polar grids have been used to simulate spiral star systems.

It is also possible to use irregular systems, although the implementation then becomes a challenge. One example of an irregular system is a presentation based on a Poisson distribution of points. The trick here is to find an efficient way to locate neighbours. Fortunately, this problem has already been solved in creating the points (see the Poisson Sampling tutorial in Dev.Mag and download code from Luma Labs). It is in fact done by using an invisible grid! Of course, visualising a state is also a bit more difficult.

Neighbourhoods

For 2D cellular automata, the most common neighbourhoods are the Von Neumann neighbourhood and the Moore neighbourhood, shown below (the blue cell’s neighbourhood is shown in red). Other neighbourhoods are also possible.

Von Neumann neighbourhood	Moore neighbourhood	Extended Von Neumann neighbourhood

The cell can be included in the neighbourhood or not. It is advisable to keep the neighbourhoods as small as possible: not only does the neighbourhood size impact performance, it also affects the number of rules to be specified and thus the tweaking complexity.

Specifying Update Rules

Specifying the rules concisely is the most challenging part of implementing CA. With very simple models, you might want to use a few if/else statements. This, however, makes it hard to modify (tweak or extend), and is therefore not suitable for more complicated systems.

It is best to implement rules as a table that can be indexed using the neighbourhood states.

Table as an array. In this method, each row represents the combined state of the surrounding cells. The final entry in the row is the new state of the current cell. Each of the remaining columns represents a state. The entries in these columns represent the number of cells surrounding the current cell in that state.

A cell surrounded with three cells in the Dry Bush state and one cell in the Burning Bush state will go into the Burning Bush state.

This table can be quite large, so that you must take care to implement lookups efficiently. Keep the table sorted, and perform binary search to lookup entries. Below are some improvements on this scheme.

Table as a tree. In the table above, it is easy to see that if there are three dry bush nodes, there are two possibilities for the other states (either 1 0, or 0 1). We can use this to implement a tree. Each level in the tree represents a state. A value in a node represents how many cells in the neighbourhood are in that state. A path from root to leaf fully describes the states of the entire neighbourhood. We do not add nodes for the last state, as it is uniquely determined by the other states.

Below is a tree corresponding to the piece of table above.

BB = Burning Bush

DB = Dry Bush

This is a more efficient implementation; unfortunately, it is also more complicated. It is inconvenient to specify the tree in a tree-friendly format (XML or s-expressions, for example). You still need to specify it in table format (space separated text file, for example) and build the tree from that.

Implementation as a hash table. Here you use tuples of state counts as keys. For example, the third entry in the table above can be accessed with the tuple (3, 0, 1). If your implementation language does not support tuples, you might need to implement them yourself. Be sure to make them immutable!

Implementation as an N-dimensional array. This is a very simple method, although it wastes a lot of memory. N is the number of states, and we use the state counts to index into the array. In our example, we will use a 3D array, so that the third element can be accessed using transition_array[3, 0, 1].

The method above can be used even when the implementation language does not provide N dimensional arrays. Simply use a 1D array, and interpret the state counts as the digits of a base-4 number. Thus, the third row will be located at 3*16 + 0*4 + 1*1 = 49.

Note that since the numbers are restricted by the requirement that they add to four, there will be unused array positions in both these last two methods. You can save memory by not using the last state to index (since it is uniquely determined by the other states). Thus, we can use a 2D array, and access the 3,0,1 rule with 3,0 alone: transition_array[3, 0]. Alternatively, if we need to use a 1D array, we can store the rule in the position 3*4 + 0*1 = 12. Although the same scheme is possible for the hash table implementation, I would advise you to use the full indexing scheme instead, since it will not waste more memory if you do, and it is somewhat easier to understand if you use full indices.

Some rules do not need the detailed tables above. For instance, in a system of only two states, where transitions depend on the count of the states of the neighbourhood, we can use a simple array, indexed by one of the states’ counts. Many other schemes are possible: in general, aim for quick lookups and easy specification.

Edges

The edges of your simulation might present some problems: edges often serve as unwanted sources, sinks or stabilisers.

A source is the origin of outward growth (states change progressively to a specific state in an outward direction).

A sink is the exact opposite:

Sources and sinks are emergent properties of the CA system, and are not physical attributes of the cells.

A stabiliser is a section that keeps neighbouring states fixed. For example, in the fire simulation (described below), there is an unnatural wall of fire around the edges at the end of the simulation.

There are several ways you can handle edges:

Ignore them. When the edge anomalies are acceptable, it is best not to complicate the model. This is only possible if the rules are specified in a way that does not require state counts to sum to the size of the neighbourhood. For example, if there is a rule “if num_states_fire > 1, then …”.

Wrap around. This solution is very elegant, easy to implement, and its effects are easy to understand. However, it is usually easy to see that the world has been wrapped in this way, and hence it is not always a visually a satisfactory solution. It is possible, however, to simulate a grid bigger than displayed area (say about double). This will hide the fact that it is wrapped, but still give you the benefits of having a world without any discontinuities.

In the CA system shown left, we want to calculate the next state of the grey cell. Two neighbours fall outside the edges. If we wrap around, we use the four red squares as shown.

Reflection. Use a neighbouring state to perform calculations.

The CA system shown left is the same situation as above. If we use neighbouring states for the cells beyond the edges, we use the state of the grey cell (twice!), in addition to the two red cells shown.

Fix their states. The edges can be given special, constant states that make the CA behave as intended. Finding the right state(s) for the edges might be difficult.

An implementation trick that applies to the last three edge-handling techniques is to include a border in your simulation. For the Von Neumann and Moore neighbourhoods, this border is only one cell thick.

This border is used for updates like all other cells, but is itself updated differently. The update rules for wrapping and reflection will simply copy states from the right places; the border for fixed states will never be updated.

In the figure, hatched cells are updated according to the special rules of the edge-handling scheme – in this case, wrapping is used, so they are simply copied from the appropriate cells. All other cells are updated normally. Note that here we do not need special logic for updating the grey cell.

A Few Examples

(Also see the downloads at the bottom of the post).

Heat

Here, continuous states (floats between 0 and 1) are used. Every cell is updated to a weighted average of its Moore neighbourhood cells: corners have a smaller weight than other cells. For this simulation, heat sinks and sources have been added: they keep the cell’s temperature at 0 and 1 respectively.

City Roads

This simulation use separate neighbourhoods for creating and destroying a road. The creation neighbourhood is an extended Von Neumann neighbourhood; the destruction neighbourhood is a Moor neighbourhood. Creation and destruction is according to probabilities based on the count of neighbouring roads.

Influence

This simulation models human emotion as influenced by others. A character has an internal state, and an externally visible state. Characters can only observer the external state of their neighbours. The internal state is then a modified: a character surrounded by angry people will itself become angry, and so on. The external state is probabilistically determined by the internal state.

Disease Spread

For our game Epidemic, Chris and I used a simple diffusion model. We used a floating-point number for the state, and represented the number of infected people in a cell. We used different colours for different ratios to visualise infection. At each time step, a cell would become infected by an amount proportional to that of the neighbourhood, with some random scaling.

Visualising States

Versions of tiles

A very effective way of making CA simulations more impressive is to have more than one tile for each state, and display them randomly. This is especially helpful for hiding regularity in off-line simulations. Tile sets with different versions for corners, edges, and central regions can drastically improve the look of the game, although it complicates choosing the correct tile version. In addition, it means you need tiles for various combinations of states. Assuming we can rotate and reflect tiles, in the fire example with three states, we need the following tiles:

3 centre tiles
3 straight-edge tiles
3 diagonal line tiles
3 cross tiles
6 corner tiles
3 half-straight tiles

This is 21 tiles for only three states. It is advisable to automate tile creation.

Tile

States can be animated. Instead of displaying a single tile, display an animation sequence. You will need to think about how to start and stop these animations so that they are seamless. For instance, you can make all animation sequences the same length, and only update states between animation sequences.

You can also increase the seamlessness of the simulation by animating transitions between states. This can be as easy as blending between two tiles over a period of time (assuming that updates are made far enough apart).

No Tiles

Instead of using tiles to represent a state of a cell, you can put objects in a cell. For instance, in a plant growth simulation, you can put a number of trees in the cell depending on the state. The advantages of this are:

• You do not have to create dozens of tiles for variety.
• The player can interact with the objects. For instance, the plant placement of trees can hinder a player’s movement.

This shows a simple forest simulation. Trees grow and die depending on the surrounding trees. This system uses a grid, but not tiles: trees can grow in random positions in tiles. For this simulation, three random coordinates have been chosen for each cell at the beginning of the simulation. Depending on the states, one, two, or three trees will be drawn at these points. Choosing a fixed set of points makes it easier to draw trees at the same coordinates every tick. A somewhat more complicated simulation can choose new points whenever a new tree comes into being.

CA Enhancements

The images below show three frames from a simple fire simulation.

The green is grass (Dry Bush), the orange is fire (Burning Bush) and the brown is burnt grass (Charcoal).

As you can see, the simulation captures some features of a fire:

• fire grows outward,
• grass only burns once,
• grass needs fire to catch fire (there is no spontaneous combustion).

Even so, the symmetry and regularity makes this fire very unconvincing. Below we look at some ideas to enhance a simple CA simulation.

Probabilistic State Transitions

Instead changing to some specific state based on the surrounding states, you can change to any state with some probability. These probabilities depend on the surrounding states, and can possibly be zero.

Adding probability to the mix can have several advantages:

It introduces variability in patterns. This means a specific (perhaps the start-up pattern) will not always produce the same results.

It can enhance the organic look of the system. Cross and squares patterns are common in deterministic CA systems. For many purposes, this is unsatisfactory (have you ever seen a fire spreading as a growing square?) Adding probabilistic state changes will make the patterns more organic and varied.

It can make discreet behaviour more continuous. For example, two common “problems” with CA systems are exponential growth (where all cells go to a certain state outwards very quickly), or sudden death (where all cells go to a certain state very quickly inwards). The problem is often that a single change in the update rules results in a exponential growth system to become a sudden death system, with seemingly no way to get a “sweet spot”. If states can only change with a certain probability, tweaking the probabilities allows you to get intermediate patterns (stable growth, equilibrium, or stable death).

It allows you to control the speed of the simulation. It is especially helpful to slow down a simulation without making it more granular (in the time dimension). It is useful to multiply all your probabilities with a global parameter that you can tweak to control the speed.

Note that in the example above, certain cells did not burn at all – per chance, they did not caught fire when their neighbours were burning. Timers (see below) can reduce or eliminate this effect.

Adding randomness can increase the emergence, and hence also unwanted emergence, so care must be taken that the simulation is always valid.

Tile Interpolation

Normally, you want your simulation to be as simple as possible. However, when such simple simulations are displayed (visualised), the results can be unsatisfactory. One way to get better visual quality is to use tile (cell) interpolation to make the simulation seem more complicated than it really is. When combined with a bit of randomisation, the effects can be very convincing.

The basic idea is to simulate only every nth cell with the CA approach. Other cells are interpolated. The number of interpolated cells can be adjusted, depending on the detail level required (balanced with the combinatorial number of tiles required!)

In the fire simulation, using only three states makes the simulation seem very primitive. Suppose that we simulate only every second row, and in these, only every second column, and fill the remaining cells depending on their neighbours. We have to create a function that will return a tile for each combination of surrounding tiles:

get_center_tile(north, east, south, west)

Interpolation rules are similar to update rules of the simulation, and hence can easily be as complicated. If this is the case, you would do better to implement more states and not do any interpolation. Interpolation should only be used if it significantly reduces the implementation complexity.

Timers

Some interesting effects can be obtained if each cell has a timer, that is, each cell only updates once the timer expires. Some rule determines when the timer is reset. For instance, the timer can be set:

• randomly,
• as a function of the state,
• as a function of the surrounding cells,
• as a function of another property of the tile (*1),
• as a function of some external process (*2).

(*1) For instance, if each cell has an associated “burning capacity” parameter, (simulating material types), a fire state timer can be set in proportion to this parameter.

(*2) For example, we can set timers according to a wind factor.

Timers can have the following effects:

• add some randomness to a simulation,
• slow the simulation down,
• smear the simulation (for example, in the fire simulation, the burning border is much thicker, and no cells are skipped).

Although you can often implement your grid cells as classes (or structs, or whatever aggregating construct your language provides) with states and timers embedded, it is simpler (on more efficient) to have separate integer grids – one for states and one for timers.

Note that if timers do not have the same expiration time, you the CA system is asynchronous and the behaviour can differ dramatically from a synchronous CA system. Thus, you should include timers early in the development process, so that their effect does not force you to re-tweak all the parameters.

If your world is big, you cannot run the CA for the whole world at every time step – it will kill performance. Below are some ideas for improving performance.

Instead of simulating the entire world, you can only simulate the part that is relevant to the game, usually the visible region of the game. A simple technique is to use a wrapped grid, similar to the grids used for the visual port in scrollable 2D games. This grid can be slightly larger than the visible grid, so that it allows cells to update naturally for a few iterations before they become visible, and so increase the illusion of a persistent world.

To hide that fact that you are only simulating a very small piece of the world, you must update “new” cells appropriately. There are several options:

Update cells to a natural starting state. In many cases, this solution is adequate, especially if you simulate a slightly larger grid than is visible.

Copy states from the existing simulation. Here you can copy neighbouring cells, or randomly from a nearby region. Copying from nearby cells promotes spatial coherence.

Initialise states statistically. Based on a neighbourhood of cells, each cell is initialised probabilistically. For example, each cell can be initialised randomly to one of three neighbouring cell’s states.

Use hierarchical modelling. Use another CA system that represents the entire world in very low detail, and use this to initialise the cells of the high detail version. Note that both the states and the update rules will differ for the high and low detail simulations. Use probabilities to prevent homogeneous systems. You can even use several levels of detail, and incorporate them as needed. Beware the extra implementation complexity, though!

This technique does not work when the world is should be changed permanently.

Instead of only simulating a part of the world, you can only update a part of the world. This requires enough memory for the whole world, but only processing power for the cells you wish to update. This method has the advantage of real persistence, thus permanent changes can be reliably simulated. Which cells you update will determine how your simulation will play out:

• Update all visible cells (or a slightly larger region). This is the most natural update scheme, and its effects are similar to partial simulation described above.
• Divide the world into discontiguous regions, and update a different region each time step. For example, divide the world into two sets: one contains every other cell, the other the rest, and update these two sets alternately. To save a significant amount of processing, you need many sets. To hide patterns in the updates that might break the illusion of a continuous world, the sets should be irregular. This method has the advantage that real things happen off screen. On the down side, it slows down the simulation significantly.
• Combine the two methods above.

Partial updating, like the inclusion of timers, lead to an asynchronous system, and hence the behaviour might be very different from that of a fully updated system.

In certain cases, simulations can be baked offline and be used in games. This means you perform the simulation before the time, and record each frame. The game then reads these frames, and plays them back.

Various schemes are possible to get more with less. For instance, you can also interpolate the simulation in time. With some thought, you can simulate events in a way that they can be looped, stacked together, and started and stopped unobtrusively. This means you can blend them with the interaction in a way that will not make it obvious that the simulation has been performed off-line.

Unwanted Emergence

Not all emergence is good, and you might find it hard to prevent certain behaviours. If the unwanted behaviour is only visual, and it occurs only rarely, it might be acceptable. However, if behaviour impacts play negatively, it must be prevented.

As with all tweaking, you should not spend too much time fumbling around. To proceed more “scientifically”, ask yourself the following questions:

• Is there a way to detect (measure objectively) the unwanted behaviour?
• Instead of tweaking your CA, can you detect-and-recover system be used?
• Is the behaviour continuous? In other words, is there a degree of badness, or is it either good or bad?
• If it is continuous, is it feasible to kick in a preventative system? (When the badness goes above a threshold, updates are processed differently).
• Is it feasible to find a set of parameters through an optimisation algorithm (for instance, genetic programming)? It might even be possible to search exhaustively for a solution.
• Is my model too complex? Try to eliminate states and simplify transition rules. You might even consider breaking your simulation into two or more independent CA systems.
• If you find yourself struggling to tweak the parameters for too long, you should perhaps look for a better model that is less sensitive to the exact parameters.

Links

General

http://cell-auto.com/ A nice breakdown of various topics related to Cellular Automata. The article on optimisation is especially interesting.

http://www.collidoscope.com/modernca/ Contains many non-traditional variations of cellular automata.

http://www.fourmilab.ch/cellab/ Application for exploring cellular automata.

Applications and Implementation

http://realtimecollisiondetection.net/blog/?p=57 Path finding using cellular automata.

http://www.ibm.com/developerworks/java/library/j-camusic/ Cellular Automata and music (interesting way of creating procedural music).

http://blog.sigfpe.com/2006/12/evaluating-cellular-automata-is.html About using lazy, infinite containers for cellular automata.

http://madeira.cc.hokudai.ac.jp/RD/takai/automa.html Graphics applications for cellular automata – lots of videos!

Game Design

http://www.gamestudies.org/0102/pearce/ An article about game design.

http://www.cp.eng.chula.ac.th/~vishnu/gameResearch/Sweetser_Thesis.pdf A thesis that looks at game design from a emergence perspective.

http://fora.tv/2006/06/26/Will_Wright_and_Brian_Eno#chapter_01 Will Wright and Brian Eno on generative systems.

Downloads

Game Maker Demos

• Requires Game Maker 7.
• Contains 9 demos.

ca_demos.zip (135 KB)