Related posts:

1. Update: Reference for Functional Equations
2. A Reference for Functional Equations

(Original Image by Valerie Everett)

It is sometime necessary to move an object in a physics simulation to a specific point. On the one hand, it can be difficult to analyse the exact force you have to apply; on the other hand it might not look so good if you animate the object’s position directly.

A compromise that works well in many situations is to use a spring-damper system to move the object.

The trick is simple: we apply two forces—the one is proportional to the displacement; the other is proportional to the velocity. Tweaked correctly, they combine to give realistic movement to the desired point.

The spring force is proportional to the difference between the current position and the position where we want the object:

Here, is a positive value called the spring constant.

As you can see, the force gets smaller as our object approaches the desired position, and becomes zero when it reaches that position. Unfortunately, in the absence of friction or drag, the velocity is not zero at this point, so the object overshoots the desired position, and moves past it. The force becomes bigger, but in the opposite direction. The object keeps on moving, slowing down, and finally starts moving in the opposite direction towards the desired position. This goes on indefinitely.

When there is friction or drag, we might be lucky enough for the system to slow down the object sufficiently so that its velocity becomes zero when the object reaches the desired position. This can be tricky to accomplish, though, and might impact the simulation environment in undesirable ways.

It is better to add a counteracting force explicitly. We add a damper force that is propostional to the velocity of the object, again opposite in direction. Here is the viscous damping coefficient, also a positive number.

We then apply the sum of the forces to our object:

The trick is to choose to get the behaviour we want.

Fortunately, this is easy. The following table summarizes how the damping constant affects behaviour. Here, is the mass of the object.

	The object oscillates.
	The object oscillates, but the oscillations die down.
	The object moves to the desired position without oscillating in minimum time.
	The object moves to the desired position without oscillating, and takes longer as c increases.

If you want to see an explanation of how this works, see the Wikipedia article on damping.

By choosing , we are left with only one parameter to tweak (the spring constant), with which we can adjust the time it will take for the object to reach the desired spot.

A simple implementation of this idea is given by the following function. The function should be called for every simulation frame, until we are satisfied that the object reached its spot:

private void MoveTo(Rigidbody rigidbody, Vector3 newPosition,
      float springConstant)
{
  Vector3 desiredDisplacement = rigidbody.position - newPosition;
  Vector3 springForce = -springConstant * desiredDisplacement;

  float viscousDampingCoefficient
    = 2 * sqrt(rigidbody.mass * springConstant);

  Vector3 dampingForce
    = -viscousDampingCoefficient * rigidbody.velocity;

  Vector3 totalForce = springForce + dampingForce;
  rigidbody.AddForce(totalForce);
}

This will work even when there is drag or friction, except that the object will move slower (in this case we can decrease the artificial damping, although it is a bit risky). When there is an external force applied to the object, the object will come to rest at some point away from the desired position. By increasing the spring constant, the distance between this point and the desired point can be made smaller. Thus, we can also use this scheme to maintain objects at a certain height, for example, which can give a rather realistic simulation of a hovercraft or even a drifting object.

Update: This is an example of a PID controller (proportional–integral–derivative controller), for which there is a C++ implementation in my Special Numbers Library.

Related posts:

1. Cellular Automata for Simulation in Games
2. Random Steering – 7 Components for a Toolkit
3. Force Field Editor v1.0

(Original image by GoAwayStupidAI).

Below are four C++ implementations of the region quadtree (the kind used for image compression, for example). The different implementations were made in an attempt to optimise construction of quadtrees. (For a tutorial on implementing region quadtrees, see Issue 26 [6.39 MB zip] of Dev.Mag).

• NaiveQuadtree is the straightforward implementation.
• AreaSumTableQuadtree uses a summed area table to perform fast calculations of the mean and variance of regions in the data grid.
• AugmentedAreaSumTableQuadtree is the same, except that the area sum table has an extra row and column of zeros to prevents if-then logic that slows it down and makes it tricky to understand.
• SimpleQuadtree is the same as AugmentedAreaSumTableQuadtree , except that no distinction is made (at a class level) between different node types.

The interfaces of all quadtrees are the same, but I did not want to extend from a base class. (Instead, a compile time check is performed on the classes, using boost concepts).

The results of the performance (on my machine!) of the trees are as follows:

As it turns out, the different implementations do not differ significantly. Constructing the nodes takes long, and not so much the calculations necessary to determine whether a node should split or not, and what data should be in the node. Had I profiled properly before I started I would not have gone through this exercise…

Between these four implementations, the NaiveQuadtree is the one I recommend; I left in the other implementations for anyone interested.

The one good thing that came from this experiment is that I found out that using a zero-augmented summed area table can increase performance quite a bit. This is useful for max-filters and other algorithms that use these tables.

You can download the source code:

Quadtree.zip (44 KB, Visual Studio).

(It requires the boost library for concept checking, but nothing else. Everything will still work if you remove all references to the boost library. If you already have boost, just hook up to the include path in Visual Studio).

…or read the online documentation.

Related posts:

1. Quadtrees
2. Quadtrees
3. A simple texture algorithm – faster code and more results

When implementing image algorithms, I am prone to make these mistakes:

• swapping x and y;
• working on the wrong channel;
• making off-by-one errors, especially in window algorithms;
• making division-by-zero errors;
• handling borders incorrectly; and
• handling non-power-of-two sized images incorrectly.

Since these types of errors are prevalent in many image-processing algorithms, it would be useful to develop, once and for all, general tests that will catch these errors quickly for any algorithm.

This post is about such tests.

The General Idea

The general idea is to exploit the fact that many algorithms satisfy invariants such as these:

some_transform(algorithm(image))
  == algorithm(some_transform(image))

For example: a 3-by-3 box blur is invariant under a vertical flip:

box_blur3x3(flip_vertical(image))
  == flip_vertical(box_blur3x3 (image))

It does not matter whether we apply the flip before or after applying the box filter—the result should be the same.

What kind of errors can the above test expose? Since it checks whether the algorithm (not the image) is vertically symmetric, the test can catch certain off-by-one errors along the vertical axis. Here is an example of a faulty implementation of the blur algorithm. Here the function sum adds all the pixel values in the rectangle x0..x1 – 1 and y0..y1 – 1.

//C-like pseudo-code
Image & box_filter3x3(Image & image)
{
  forXY(image, x, y)
  {
    x0, y0 = max(0, x – 1), max(0, y – 1)
    x1, y1 = min(x + 1, image.width), min(y + 1, image.height)
    result = sum(image, x0, x1, y0, y1) / ((x1 – x0) * (y1 – y0));
  }

  image = result;
  return image;
}

Can you spot the error? The actual window used is only 2 by 2. But how does the test catch this error? To see why the test will fail, notice the difference in windows for the top-left and bottom-left pixels:

Top Left:

x0, y0 == 0, 0
x1, y1 == 1, 1

Bottom Left

x0, y0 == width – 2, height – 2
x1, y1 == width , height

As you can see, the two windows have different sizes! The bottom-left pixel window is 1 by 1, but the top left pixel is 2 by 2. Thus flipping the image before and after applying the box-blur gives different results.

The correct algorithm is

//C-like pseudo-code
Image & box_filter3x3(Image & image)
{
  forXY(image, x, y)
  {
    x0, y0 = max(0, x – 1), max(0, y – 1)
    x1, y1 = min(x + 2, image.width), min(y + 2, image.height)
    result = sum(image, x0, x1, y0, y1) /
      ((x1 – x0) * (y1 – y0));
  }

  image = result;

  return Image;
}

Now the windows for the top and bottom left pixels are:

x0, y0 == 0, 0
x1, y1 == 2, 2

Bottom Left

x0, y0 == width – 2, height – 2
x1, y1 == width , height

We can see the window sizes are now the same.

Transforms

Here is a list of suggested transforms and the types of errors they can catch:

Many algorithms should (in theory) also be invariant under scaling. However, because of the complexity involved in interpolation and sampling, I do not recommend using scaling for testing. It is quite difficult to determine under some interpolation or sampling scheme whether an algorithm should in fact be (exactly) invariant or not under scaling. This of course also applies to other transforms that rely on sampling or interpolation.

You will probably also select transforms that are more specific to the algorithms you are implementing. Keep these as simple as possible—not only to avoid implementation errors, but also to avoid subtle misconceptions: it must be easy to see (or well established) that a certain algorithm is invariant under a transform. Test your own transform functions aggressively.

Generic Solutions

Writing test code for each transform is easy, but tedious (and hence error prone). Here I show how a generic test can be implemented using C++ macros or functional programming. The generic test can be used to test whether any algorithm is invariant under a given transform.

Macros Implementation

Probably the easiest way to implement a general test-function in C++ is to implement it as a macro. Here is how this looks:

//C++
#define TEST_TRANSFORM_INVARIANCE(image, transform, command) \
  do { \
    Image image; \
    make_test_image(image); \
    \
    Image original_image(image); \
    command; \
    \
    Image image1(image); \
    transform(image1); \
    \
    image = original_image; \
    transform(image); \
    command; \
    \
    if (image != image1) \
      report_useful_failure_message(...); \
  } while(0)

We can now use this macro like this:

TEST_TRANSFORM_INVARIANCE (image, flipVertical, my_algorithm(image));

The nice thing about the macro approach is how easy functions that take any number of parameters can be tested:

TEST_TRANSFORM_INVARIANCE (image, flipVertical, my_algorithm(image, 10, 20));

Notice that the first parameter is a variable name. The macro declares this variable, and the user can use it in the command parameter to pass the image to the algorithm.

Functional Programming Implementation

If macros are not available, or using them is not desirable, you can still use a generic approach. Unfortunately, we need to define functions with a consistent argument list for each test we want to perform, which means we might need to write more code in some languages.

The general test function can be defined as follows in C++:

//C++
test_transform_inavriance(
  Image & (* function)(Image &),
  Image & (* transform)(Image &))
{
  Image image;
  make_test_image(image);

  Image original_image(image);
  function(image);

  Image image1(image);
  transform(image1);

  image = original_image;
  transform(image);
  function(image);

  if (image != image1)
    report_useful_failure_message(...);
}

To use it with a single argument function, we simply call it like this:

test_transform_inavriance(flip_vertical, my_algorithm);

To use it with a function that takes more than one parameter, we need to define a function to call the extra parameters:

// C++
Image & test_my_algorithm_wrapper(Image & image)
{
  my_algorithm(image, 10, 20);
}
...
test_transform_inavriance(flip_vertical, my_algorithm_wrapper);

In languages that support closures, you can wrap functions more cleanly so that you do not need to define a function for each test command. For example, in Python:

# Python
def wrap(fn, args):
  def wrapper(image):
    fn(image, *args)
  return fn
...
test_transform_inavriance(flip_vertical, wrap(my_algorithm, 10, 20));

Test Bundles

Both the macro and functional programming approaches allow you to combine tests of different transforms in convenient functions:

// C++
#define TEST_ALL(image, command)\
  do { \
    TEST_TRANSFORM_INVARIANCE (image, flip_vertical, command); \
    TEST_TRANSFORM_INVARIANCE (image, flip_horizontal, command); \
    ...
  } while(0)

// C++
void test_all(Image & (* fn)(Image &))
{
  test_transform_inavriance(flip_vertical, function);
  test_transform_inavriance(flip_ horizontal, function);
  ...
}

Test Images

It is important that your images do not destroy the very asymmetry that you are trying to expose. For instance, if we ran the test at the very beginning of this post with an all-zero image, the test will have passed. The image should be asymmetric under the transform we are using, that is, the following must hold:

transform(image) != image

For different transforms, this means different things. For example, for flip_diagonal, the image must not be square, for inverse, the image must not be the constant grayscale image 0.5.

It is useful to build an extra test into our test macro or function, to make sure that we do not inadvertently break this requirement. Here is the modified macro:

// C++
#define TEST_TRANSFORM_INVARIANCE(image, transform, command)\
  do { \
    Image image; \
    make_test_image(image); \
    \
    Image original_image(image); \
    transform(image); \
    \
    if(image == original_image) \
      report_unsuitable_test_image(...); \
    \
    image = original_image; \
    command; \
    \
    Image image1(image); \
    transform(image1); \
    \
    image = original_image; \
    transform(image); \
    command; \
    \
    if (image != image1) \
      report_useful_failure_message(...); \
  } while(0)

For many algorithms, using an image of noise (independent for each channel) works well. Make sure the image have unequal dimensions that are prime numbers (for example, 11 and 13). Using prime numbers ensures that the two dimensions will not have any common divisors, which makes certain kind of tests for recursive algorithms (such as quad-tree compression) more robust.

Reporting Failures

In the example implementation, we have the line

if (image != image1) report_useful_failure_message(...);

This is a somewhat simplistic scheme. What is a useful failure message? Of course, the message must report the test that failed. But it should also give information about how the test failed:

• whether it was a mismatch of the number of channels, image dimensions, or pixel values;
• the number of channels and image dimensions;
• the number of pixels for which the test failed; and
• the first pixel location where the test failed, and the expected and actual pixel values.

These bits of information can help us hunt down the cause. For example, if the test data is an 11-by-13 image, and the test fails for 13 pixels, the error is most likely a border problem caused by an off-by-one error. If the first pixels that fails is 0, 0, and the test fails for (almost) all pixels, and the pixel values differ only by sign, we can guess that there is some sign error made that affects the entire image.

The actual test will thus be a bit more complicated, delegating to a function such as this:

test_image_equality(
  Image & image1,
  Image image 2,
  FailureInfo & info)
{
  info.has_failed = false;
  info.failed_pixels_count = 0;
  info.dimensions_image1 = image1.dimensions;
  info.dimensions_image2 = image2.dimensions;
  info.channels_image1 = image1.channels
  info.channels_image2 = image2.channels

  if(image1.channels != image2.channels)
  {
     info.has_failed = true;
     info.failure_type = CHANNEL_MISMATCH;
     return;
  }

  if(image1.dimensions != image2.dimensions)
  {
    info.has_failed = true;
    info.failure_type = DIMENSION_MISMATCH;
    return;
  }

  forXY(image, x, y)
  {
    if(abs(image1(x, y) - image2(x, y)) < THRESHOLD)
    {
      if (!info.has_failed)
      {
        info.has_failed = true;
        info.failure_type = VALUE_MISMATCH;
        info.first_failed_x = x;
        info.first_failed_y = y;
        info_first_failed_image1 = image1(x,y);
        info_first_failed_image2 = image2(x,y);
      }

      info.failed_pixels_count++;
    }
  }
}

In the test macro or function, we call the function like this:

...
FailureInfo info; \
test_image_equality(image, image1, info); \
if (info.has_failed) report_failure(info); \
...

Using Third-Party Image Libraries

It is important that you understand the transforms you use for testing very well. It is therefore a good idea to use you own transforms, and not those of a library. This way, unknown design choices in the library cannot bite you. Many image libraries do not for example specify how they handle borders, division by zero, rounding, etc. These details are often not important when using the algorithms in a task—but they can make a test fail (or succeed) when it should not.

The same danger exists when you use third party algorithms as building blocks for your own. In this case testing your code will be harder. As a starting point, consider testing whether third party algorithms satisfy the invariants you expect to hold. If they do not, you may need to modify your test code to accommodate for the discrepancy. If they do, you can be a little more confident that your tests will be reliable.

(Yes, it is not normally recommended to test other people’s code with unit tests. However, this is a once-off test to make sure that your own tests are reliable, and since it is so easy to implement (since you already have the test macro / function and transforms), I do not see the harm. Just keep the third-party library test separate from your own.)

Defining Correctness—When are two images equal?

Correctness of an image algorithm is a subtle issue, mostly because the discrete, quantized, finite model of image computation is inherently an approximation of the “real” thing. The question is too broad to tackle here (and frankly, I do not have enough knowledge to do it), so I will focus on aspects that are specifically relevant to this kind of testing.

Let us start with an example of an essentially correct algorithm that fails a test that we expect that it should not:

// C++
Image & threshold(Image & image)
{
  forXY(image, x, y)
    image(x, y) = (image(x, y) <= 0.5f) ? 0 : 1;

  return image;
}

We expect this simple algorithm to be invariant under the inverse transform. However, it is not. Consider a pixel value of 0.5: the inverse of this is 0.5. So both the pixel and its inverse will map to zero, and hence the following will not be true:

inverse(threshold(image)) == threshold(inverse(image)) //not true
inverse(threshold(0.5f)) == inverse(0) == 1 //left side
threshold(inverse(0.5f)) == threshold(0.5f) == 0 //right side

It looks like the problem is easily solved by changing the algorithm:

// C++
Image & threshold(Image ℑ)
{
  forXY(image, x, y)
    image(x, y) = (image(x, y) < 0.5f) ? 0 : 1;

  return image;
}

But, for floating point numbers on a machine, there is another number 0.5 – epsilon such that the following algorithm is equivalent to the one above:

// C++
Image & threshold(Image & image)
{
  forXY(image, x, y)
    image(x, y) = (image(x, y) <= 0.5f - epsilon) ? 0 : 1;

  return image;
}

But for this algorithm our invariant does not hold for a pixel with the value 0.5 – epsilon! Since the last two algorithms are equivalent, it means the modified algorithm is also not correct in this strict sense. In fact, it is not possible to implement an algorithm that is correct in this sense (using floating-point numbers).

But we feel that any of the above implementation is correct, and hence our test is wrong. So how do we change it, and should we?

We can change one (or more) of three things:

• the test data
• the test transform
• how we test for equality in images

Because of the simplicity and general usefulness of both the transform and the test data, I feel it is best to leave these as is. That means we have to change the measure of equality. Since we use random data, we might want to use the “is probably equal” measure—that is, two images are equal if a certain percentage of pixels are equal (within some threshold). Since the likelihood of generating 0.5 is low, our test should pass. I do not like this idea for two reasons:

First, it is not generally useful. For example, another test might fail when the borders are not equal (a much larger number of pixels), even though we feel that it must not for some specific algorithm. We then need another equality measure to handle that case.

Second, it is bound to hide a class of errors where a small number of pixels are incorrectly processed, and we do not want that.

Instead, we would like some measure that is generally useful, but can be tweaked for the specific algorithm to overcome those cases where we expect the algorithm to fail.

One way to do this is to use equality masks: a simple bit array that says whether we are interested in that pixel or not. In this case, we need the mask to be constructed from the test data; here is how we do it:

//C++
BitImage & make_ignore_value_mask(
  BitImage & mask,
  Image & image,
  float value)
{
  forXY(image, x, y)
    mask = image(x, y) == 0.5f;
}

We then change our macro to this:

#define TEST_TRANSFORM_INVARIANCE (image, mask, transform, command, mask_command) \
  do { \
    Image image; \
    make_test_image(image); \
    \
    Image original_image(image); \
    transform(image); \
    \
    if(image == original_image)
      report_unsuitable_test_image(...); \
    \
    image = original_image; \
    command; \
    \
    Image image1(image); \
    transform(image1); \
    \
    image = original_image; \
    transform(image); \
    command; \
    \
    BitImage mask(image.width, image.height); \
    mask_command; \
    \
    if (mask * image != mask * image1)
      report_useful_failure_message(...); \
  } while(0)

We call the macro like this:

TEST_TRANSFORM_INVARIANCE (image, mask, inverse, threshold(image), make_ignore_value_mask(mask, image, value));

The functional programming solutions are similarly modified:

// C++
test_transform_inavriance(
  Image & (* function)(Image &),
  BitImage & (*mask_creation_function) (BitImage & mask, Image & image),
  Image & (* transform)(Image &))
{
  Image image;
  make_test_image(image);

  Image original_image(image);

  if(image == original_image)
      report_unsuitable_test_image(...); 

  image = original_image;
  function(image);

  Image image1(image);
  transform(image1);

  image = original_image;
  transform(image);
  function(image);

  BitImage mask;
  mask_creation_function(mask, image);

  if (mask * image != mask * image1)
    report_useful_failure_message(...);
}

We need a wrapping function to call our extra argument:

BitImage & void wrapped_ make_ignore_value_mask(
  BitImage & mask,
  Image & image)
{
  return make_ignore_value_mask(mask, image, 0,5f);
}

Now we can call our test function like this:

test_transform_inavriance(image, mask, inverse, threshold,
  wrapped_make_ignore mask(mask, image));

In languages that support closures such as Python, we can again use the wrapping technique as before:

test_transform_inavriance(image, mask, inverse, threshold,
  wrap(make_ignore mask(mask, image), 0.5))

Here are a few suggestions for useful mask types:

Be careful for inadvertent ignore_all masks. For example, the ignore_border mask might be all zero when the border is thicker than half the smallest image dimension. It is a good idea to build in a test to check that the mask is not all-zero in the test macro or function.

Also watch out when transforms or algorithms change image dimensions—make sure the mask is constructed correctly.

We have now answered how we can change the test to pass for the thresholding example (and more generally); it remains to answer whether we should. My feeling is: only when you have to. Generally, I test without regard for such subtle issues. When a test fails, I carefully try to understand whether it is because the algorithm is fundamentally wrong, or whether it is just an obscure instance that makes the test fail. In the latter case, I will amend the test with the necessary masks.

Update: Here is another example of how an invariance test can fail even when the algorithm is fundamentally correct. Consider a region quadtree compression scheme. Suppose we want to compress a 5&times 5 image, and suppose the algorithm divides the image in four rectangles. Because 5 is odd, the rectangles will differ in size. The biggest rectangle should always be in the same spot, regardless of the reflection; thus the algorithm should fail when testing pixels along the center. We cannot really address this issue with masks. To me it is unclear how to handle this at all (and clearly, we do not want to limit our tests to powers-of-two images only).

Warning: A False sense of Security

When you have a large group of transforms to construct invariance checks from, you can easily test for a whole bunch of errors with just a few lines of code. If the tests pass, it is easy to think that the algorithm is correct. It may not be:

• Perhaps the test image is unsuitable. This may be by coincidence if the test image is random.
• Any particular test does not catch all errors of a certain type.
• All transforms do not cover the entire set of possible errors.

A set of invariance tests does in general not prove the correctness of a particular algorithm. It merely exposes certain kinds of errors. Therefore, additional tests for correctness are necessary.

Invariance Test Checklist / Summary

The basic structure of the test macro / function is as follows:

1. Create test data.
2. Calculate transform(algorithm(test_image)).
3. Calculate algorithm(transform(test_image)).
4. Compare these results and report any failures.

Remember:

• Check that the test data is changed by the transform.
• Use masks to cater for peripheral special cases where certain pixels, values, or regions in an image should be ignored for comparisons.
• Check that your masks are not all zero (sanity check).
• Report useful information when a test fails:
  • what kind of failure (dimension mismatch, channel mismatch, value mismatch);
  • in the case of value mismatch, the number of pixels that are mismatched and the values of the first pixels from each image that fail to match.

More on Unit Testing for Image Processing

• MTEST – A unit test harness for MATLAB code
• Testing Graphics Code

How do You Test Your Image Algorithms?

Let me know in the comments

Related posts:

1. Quadtrees
2. Simple, Fast* Approximate Minimum / Maximum Filters
3. A simple texture algorithm – faster code and more results

*Fast = not toooo slow…

For the image restoration tool I had to implement min and max filters (also erosion and dilation—in this case with a square structuring element). Implementing these efficiently is not so easy. The naive approach is to simply check all the pixels in the window, and select the maximum or minimum value. This algorithm’s run time is quadratic in the window width, which can be a bit slow for the bigger windows that I am interested in. There are some very efficient algorithms available, but they are quite complicated to implement properly (some require esoteric data structures, for example monotonic wedges (PDF)), and many are not suitable for floating point images.

So I came up with this approximation scheme. It uses some expensive floating point operations, but its run time is constant in the window width.

The crux of the algorithm is the following approximation, which works well for , .

I am not going into the mathematical details of why this a fair approximation.

To turn this into a max filter is straightforward:

1. Calculate the power image

Raise each pixel to the power p.

2. Calculate a summed area table

The summed area table (or integral) of an image is a table containing in each cell (x, y), the sum of all the pixels top and left of (x, y); for the power image it is given by:

Update: This can be efficiently be calculated using the following recurrence:

and  if or .

3. Calculate the window mean

Use the cumulative image sum to calculate the mean in the window. Here, is the window radius, A stands for average:

4. Calculate the approximate maximum

That’s it!

The minimum filter is implemented by finding the maximum of the inverse image and then inverting the result. This approach uses two more passes through the image; so far I could not find a direct approximation for the minimum value. (Update: It seems using gives an approximation for finding the minimum of a set of values.)

Notes

• For very large images, you might run into floating point issues when calculating the cumulative sum image. You can reduce this effect somewhat by subtracting the (total image) mean of the power image from each pixel before calculating the image sum.
• For the image sizes I worked on, the value  was suitable. The higher , the more accurate results, as long as you have no underflow.
• Raising powers is a slow operation. By choosing a suitable value of , you can use a suitable, faster approximation for calculating the power.
• Update: It is possible to implement this as a two-pass algorithm, first working on columns, then on rows. This requires constructing two tables, which makes the algorithm slower. However, using this approach allows bigger images to be processed before overflow in the summed are tables becomes a problem.

Initially I thought that once the maximum can be found, we can also find the second maximum, and so on, so that any statistical order filter can be constructed.

The basic idea was that if is the maximum, we could find the second largest value like this:

The idea is that we remove the maximum value from the sum, so that the result must yield the second largest element. Of course, we do not use the actual maximum, but the approximation—and this is where the formula fails. When we use the approximate maximum, it is easy to show that the formula above will yield the exact same approximate maximum, and not the second largest element.

I should have know it seemed too good to be true!

Update: while looking for links for this post, I found this article: Robust Local Max-Min Filters by Normalized Power-Weighted Filtering (PDF) which initially looked like it was essentially the same thing. Although it is very similar, it is in fact different. The approximate maximum is given by the value:

for .

Using gives an approximation to the minimum.

The nice thing about this approximation is that it seems to work well for all real , not just for .

Related posts:

1. Region Quadtrees in C++
2. A simple texture algorithm – faster code and more results
3. A Simple Trick for Moving Objects in a Physics Simulation

Many textures used for 3D art start from photographs. Ideally, such textures should be uniformly lit so that the texture does not interfere with the lighting applied by the 3D software. Often, lighting artefacts must be removed by hand. This can be tedious and time consuming.

The tool provided here aims to automate this process. It is still in an experimental phase, so it is very crude. Below you can see some of the before and after pictures.

I will say more about how these work in a later post.

If you want to give it a try, go here to download (caveats: command-line only; Windows only; requires ImageMagick.)

The tool can also improve artefacts in tileable textures:

BEFORE (distracting horizontal striping)

AFTER (reduced striping)

Related posts:

1. How to Turn XSI Mod Tool into a Level Editor for your XNA Games: Example Updated for XNA 3.0.
2. Guerrilla Tool Development
3. How to Turn XSI Mod Tool into a Level Editor for your XNA Games: Updated for XNA 3.0.

Download

Related posts:

1. A simple texture algorithm – faster code and more results
2. A Simple Procedural Texture Algorithm – More Results and Code

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 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://code.google.com/p/specialnumbers/).

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). Notice that the points do not lie in a straight line, the agent will first move to point (1, 1, 0), and then turn before moving further. An onlooker will not understand this strange behaviour.

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 of 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 a 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: one signal’s 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 a 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!

Related posts:

1. Random Steering – 7 Components for a Toolkit
2. 5 Tips for Prototyping Slow Algorithms
3. Cellular Automata for Simulation in Games

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

Related posts:

1. How to Turn XSI Mod Tool into a Level Editor for your XNA Games: Updated for XNA 3.0.
2. How to Turn XSI Mod Tool into a Level Editor for your XNA Games: Example Updated for XNA 3.0.

(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.

Related posts:

1. Generating Random Numbers with Arbitrary Distributions
2. Cellular Automata for Simulation in Games
3. 5 Tips for Prototyping Slow Algorithms