null program

30 Nov 2007

Traveling Salesman Problem by Genetic Algorithm

Download Source: genetic.tar.gz (6.19KB)

Here is another project for my artificial intelligence class. I wrote a generic genetic algorithm class in C++ and then applied that class to the traveling salesman problem. A genetic algorithm more tuned to the traveling salesman problem would work better.

This particular implementation can use up to 16 points defined in the weight matrix stored in travel.dat. This weight matrix can either be defined by hand or generated using gendat.m from a list of points stored in points.txt. A chromosome is 64 bits wide, which is 16 points with 4 bits each. To make sure that every possible chromosome is a valid solution, the points are selected out of a circular queue. Every 4 bits describes how far along the queue to walk before pulling out a point. With the circular queue, the chromosome could be as short as 50 bits, but I was trying different things and 64 bits is the simplest way to represent a solution. Here are some sample chromosomes:

1101100010100000011001010101010111000110101100111111000010111010, 9315
0000110100000001000010010101010111011000101100010110000010111010, 9339
1010001100010001010010011001010010100000101100011111000010111110, 9410
0101001000000011010010011001010011010100101011100111000111011110, 9355
0101001001000000010110100001010011010110101101000001000101000110, 9349
0101011100010011000010011011011011010101101100011111010010111110, 9311
0000111101000001000010011001010010110100101100011111100010000010, 9350
1000001111100000010110101011011011100000111101010111000010111010, 9428
0000111011000001000010010111011111111110101100010111000001000111, 9448

The second number is the fitness value of the chromosome (10000 - path length). Below is a path found after 20000 iterations:

gnuplot plot

It takes many iterations to find a reasonable solution and it never finds a *really* good solution. This is because each node in the chromosome depends on every single node before it. This is terrible for a genetic algorithm, but it's really the best I could think of when using a generic genetic algorithm class like this.

A much better method would have the genetic algorithm actually know about the problem at hand, working with nodes rather than bits. Breeding would make cuts on nodes. Mutations would swap single nodes. Perhaps this can be written another time.

[] permanent link

20 Nov 2007

Noise Fractals and Clouds

I was reading about fractal terrains and wanted to give it a shot. As usual, I like to try things out in GNU Octave first by writing up a quick prototype, and, if it is still interesting enough, I will move to another language for better performance (Octave can be frustratingly slow sometimes). Additionally, when it comes to matrix manipulation, the Matlab language tends to be the most concise and powerful. Most functions and operators work on entire matrices or vectors at a time, avoiding many iteration loops, and, more notably, matrix notation built right into the language. This is powerful in the same way regular expressions are built into Perl. The problem is when your data shouldn't be in matrix form, and the language, lacking any other kinds of data structures (especially hash tables! (ignoring the struct hack)), forces you to fit your problem into matrices. Also, the implementations of the Matlab language tend to be very slow.

</rant> Anyway, back to noise.

So, the first and easiest noise algorithm I found was the diamond-square algorithm. Basically, it is noisy interpolation applied recursively. All of the noise adds up to provide something that may be good for height maps, possibly providing procedurally generated terrain for a game. Here is an example/pretty picture. Imagine this as being terrain,

heightmap

You can see we have a sort of mountain going on in the middle. Here is the "plasma fractal" view of our noise.

bird's eye view

These were generated using this Octave code. I believe that this is both the fastest and most concise way to do this in Octave. You can see we are throwing away a lot of random values that we pulled, but at the same time avoiding loops. You will need at least version 2.9 of Octave because, as far as I know, Octave 2.1 doesn't have interp2. And, unfortunately, interp2 is much slower than it probably could be.

function m = diamond_square (m, i, c)
  if isempty (m)
    m = zeros (2);
  end
  
  for k = 1:i
    m = interp2 (m, 1);
    
    ## Define points we want to randomize
    gridmap = ones (size (m));
    gridmap(1:2:end, 1:2:end) = 0;
    gridmap = find(gridmap);       # Makes Octave happy
    
    ## Define random values to be added
    randmap = c * randn (size (m));
    
    m(gridmap) += randmap(gridmap);
    
    c = c / 2;
  end
end

Something of note here compared to other implementations of diamond-square, my code above adds Gaussian noise (with Octave's randn) rather than uniform noise. The scaling variable c is actually adjusting the standard deviation of our noise and not the limits of the noise. I imagine this makes the terrain more natural, as Gaussian noise tends to approximate real noise well.

The first argument provides a base terrain to work from. With this you can define a mountain, islands, hillside, etc. When an empty matrix is provided, the terrain will be grown from flat ground. The second argument decides how many iterations we are going to run. The above example performs 6 iterations on flat terrain. The last argument decides how dramatic the terrain is, which doesn't have much meaning when starting from flat terrain. The example above was produced with two calls like this (one for each image),

octave> mesh(m=diamond_square([], 6, 0.1))
octave> imagesc(m); colormap(gray);

We can add some water, or perhaps lava or something, by choosing a water height and chopping off everything below it. I will choose -0.1 as our water level.

octave> w = m;
octave> w(w < -0.01) = -0.01;
octave> mesh(w)

Water

Here is an example of providing some base terrain. Let's put a mountain in the corner of the map,

octave> a = [0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 1 0; 0 0 0 0 0]
octave> mesh (b = diamond_square (a, 4, 0.15))

Mountain base

Becomes,

Mountain noise

Octave-Forge has a surf function that should make this look nicer, but it doesn't. If we want to render this to make it look nice, we will need another program to do it. Let's do that another time.

Another way to make noise is Perlin noise. Generally, this noise will be better and more useful than diamond-square noise (examples later). I personally like this person's approach: Perlin Noise. It is easy to implement and understand, though it doesn't have the all advantages you get from generating Perlin noise the normal way.

If you don't feel like clicking through and see the nice introduction to Perlin noise there, here is the idea of this approach: we create different frequencies of noise and mix them all together, giving more weight to lower frequency noise than high frequency noise. Below, you can see that I generated (Gaussian) noise in different frequencies, smoothing by spline interpolation. Then we add this all together to get our final Perlin noise.

Noise + Noise + Noise + Noise + Noise + Noise = Noise

Here is the Octave code I am using to generate Perlin noise,

function s = perlin (m)
  s = zeros(m);    # output image
  
  w = m;           # width of current layer
  i = 0;           # iterations
  while w > 3
    i = i + 1;
    d = interp2(randn(w), i-1, "spline");

    s = s + i * d(1:m, 1:m);
    w -= ceil(w/2 - 1);
  end
  
end

The first and only argument gives the side length of the image you want to generate - it only generates square images. This is an extremely simple approach, as there are no parameters to adjust to define the nature of the noise you want to generate. Fortunately, what I am going to try next works fine without these parameters (or with the default, hard-coded parameters if you wish).

Continuing with using some ideas from Hugo Elias (I think that's his name), I am going to use this noise to attempt to create some realistic looking clouds.

First, we will generate some noise. Let's generate some Perlin noise,

octave> n = perlin (200);

Perlin noise

Now, this doesn't look too much like clouds. To get clouds, we can apply an exponential function to the data and set anything below 0 to 0 (contrast stretching). If we scale our noise between 0 and 1, the function will look like this,

Exponential function

To adjust the clouds, we can move this function left and right across the x-axis. We can adjust the "time constant" of the function to change the sharpness of the clouds. Here is the function I wrote to convert the diamond-square or Perlin noise into cloud cover. The output is scaled between 0 and 255 to aid in image output.

function a = get_clouds (a)
  ## Scale a between 0 and 1
  a = a - min(a(:));
  a = a / max(a(:));
  
  ## Parameters
  density = 0.5;
  sharpness = 0.1;
  
  a = 1 - e.^(-(a - density) * sharpness);
  a(a < 0) = 0;

  ## Scale between 0 to 255 and quantize
  a = a / max(a(:));
  a = round(a * 255);
end

Now run this on our noise,

octave> c = get_clouds(n);

Perlin Clouds

Well, that looks a bit more like cloud cover. We just need to apply a colormap to this. I wrote this colormap function for this purpose,

function c = cloud_cmap ()
  c = [0.25 0.25 1];
  for i = 2:256
    c(i, :) = (i-2)/256 * (1 - c(2, :)) + c(2, :);
  end
end

Apply the colormap,

octave> imwrite("clouds.png", c+1, cloud_cmap)

Clouds

Wow, that looks pretty good now. We could improve this by fading to transparent rather than blue. Then do a projective transformation on the clouds and lay them over top a blue gradient. You could do this with image editing software such as GIMP, or you can continue to use Octave like I did at the end of this post.

So, in about 30 lines of Octave (including code on the interactive command line) code we could generate some kind-of realistic looking clouds. Here, I will use the cloud demo to show one particular advantage of Perlin noise over diamond-square (or at least my implementation). Here are some diamond-square clouds,

Clouds Clouds Clouds

And here are some Perlin clouds, which I think look a bit better,

Clouds Clouds Clouds

Notice the straight lines in the diamond-square clouds? You can see it right in the middle of the first image. This comes from the way that the 2-dimensional interpolation stretches the noise in the vertical and horizontal directions, drawing these lines out. This may only be apparent in my implementation, as I am probably missing the "diamond" part of the algorithm. Oh well.

Anyway, to take the clouds a bit further, you can use Octave's imperspectivewarp to apply a perspective transformation to the cloud images. I put some code together that does this transformation as well as adds a gradient,

function i = pers_clouds (n)
  w = size(n, 1);
  c = get_clouds (n);
  t = -pi/32;
  P = [cos(t) sin(t) 0; -sin(t) cos(t) 0; 0.001 0.002 1];
  
  ## Perspective transformation
  pc = imperspectivewarp (c/255, P, "cubic");
  pw = size(pc, 1);
  ph = size(pc, 2);
  
  ## Create and combine background gradient
  [dump back] = meshgrid(1:ph, 1:pw);
  i = pc * 4 * pw + back;
  
  ## Fit between 0 to 255 for image
  i = i - min (i(:));
  i = round (i / max (i(:)) * 255);
  i(isnan(i)) = 0;
end

Provide either Perlin noise or diamond-square noise and it will return an image that you can write out to a file,

octave> n = perlin (1000);
octave> imwrite ("clouds.png", pers_clouds(n) + 1, cloud_cmap)

After cropping the image with something like kolourpaint (as I did below), you get,

Clouds by Perspective

Then we could try sticking it in an image,

Eiffel Tower Before => Eiffel Tower After

Because I snatched this image from Wikipedia (they won't miss it), this image is licensed under Creative Commons Attribution 1.0.

[] permanent link

13 Nov 2007

Neural Network Blackjack Game

I have a really nice post coming up soon (it is mostly written up too!), but I want to spend more time on it to make good. Since I need something to post this week, and I still have material to move here, I have dug this up from my old website.

Get the source code (C++): neural.tar.gz (16.27KB) (GPL)

This is a neural network I wrote for an artificial intelligence class I took about a year ago. It includes a stand alone neural network class you can easily use in your own C++ program. Built around this neural network is a simple version of Blackjack (hit or stand only). You can play the neural network at Blackjack after it has finished training, which can take up to a minute or so.

My implementation of the neural network is a really simple one, using only back propagation, but it still has some neat surprises in it. When I was working with it, I would use a simple GNU Octave script to watch what was going on in real time.

Neural
Network Plot

The x-axis is the number of iterations in tens of thousands and the y-axis describes how often the neural network plays exactly the same way a cheater would at the same seat. A cheater is defined as someone who knows what card is next and can play perfectly. Note: the neural network itself is not cheating. At the end, it agrees with the cheater about 83% of the time. This is the script that reads the neural network output. For those who don't recognize the she-bang, save this to the file plotdat and set the execution permission (chmod +x plotdat).

#!/usr/bin/octave -qf
# Usage: 
#  plotdat filename [loop] [last index]

dat = dlmread(argv{1});
x = length(dat);
loop = 0;

if (length(argv) > 1)
  if (argv{2} == 'loop')
    loop = 1;
  else
    x = argv{2};
  end
end

plot(dat(1:x));

while(loop)
  dat = dlmread(argv{1});
  plot(dat);  
  sleep(1);
end

pause;

When you run the program, a blank, dumb neural network is created. It doesn't know how to play blackjack. In fact, just as it seems with my sister at times, it doesn't know how to do anything at all. It's like a newborn baby, but without any instincts. To make this neural network useful, it is taught to play blackjack by running it very quickly through about one million games, all tutored by a teacher who gets to cheat by looking ahead in the deck. This allows the teacher to always know the proper move.

The training works by giving the network an input and the desired output (determined by the cheating). The network adjusts its internal weights depending on the error between what its current output for the input and the desired output. In doing this, the neural network picks up on the statistical nature of its inputs and will pick up on patterns.

When setting up the neural network in your own program, the tricky part is determining how to encode the inputs so that patterns will be found. In this case, I provide 3 integers to the network: its lower bound score, its best score, and the visible opponent score. By lower bound and best score, I am talking about Aces. All aces are treated as 1's in the lower bound and ideal (highest without bust) values in the best score. Scores never exceed 31, so we can encode this in 5 bits each for a total of 15 bits. We only need 1 output bit: hit (0) or stand (1). The network looks something like this, but with 30 hidden layer neurons and a lot more connections.

Neural
Network

Here is an example run of the blackjack game,

Creating network ... created!
Training ...
Win, lose, total, c%, win%: 1536, 7640, 10000, 42.1168% , 15.36%
Win, lose, total, c%, win%: 3728, 14438, 20000, 55.8875% , 21.92%

... (lots of output for a few seconds) ...

Win, lose, total, c%, win%: 203744, 700338, 980000, 82.0455% , 20.3%
Win, lose, total, c%, win%: 205849, 707461, 990000, 82.5488% , 21.05%
Win, lose, total, c%, win%: 207981, 714515, 1000000, 82.3059% , 21.32%

Begin game:
Computer is dealt: (hole) 5 
Computer is dealt: (hole) 5 4 
Computer stands.
Your hand: 2 8 
Hit? (y/n): y
Your hand: 2 8 3 
Hit? (y/n): y
Your hand: 2 8 3 10 
You bust with 23
Computer wins with 19 against your 23
Play again? (y/n) : y

Begin game:
Computer is dealt: (hole) 4 
Computer is dealt: (hole) 4 10 
Computer stands.
Your hand: 10 A 
Hit? (y/n): n
Your hand: 10 A 
You win with 21 against computer's 20
Play again? (y/n) : n

Computer wins, losses, push: 1 (50%), 1 (50%), 0 (0%)

If you decide to try using the neural network in your own program, be sure to play with different sized networks to see what works best. In my implementation, a small network would be 1 layer with 3 nodes and a large network would be 3 or 4 layers with hundreds of nodes. Bigger networks may not work well at all, and there is no way to know what size is best other than trial and error.

[] permanent link

06 Nov 2007

Iterated Prisoner's Dilemma

Animation

I was reading about the prisoner's dilemma game the other day and was inspired to simulate it myself. It would also be a good project to start learning common lisp. All of the common lisp code is available in its original source file here: prison.lisp. I have only tried this code in my favorite common lisp implementation, CLISP, as well CMUCL.

In prisoner's dilemma, two players acting as prisoners are given the option of cooperating with or betraying (defecting) the other player. Each player's decision along with his opponents decision determines the length of his prison sentence. It is bad news for the cooperating player when the other player is defecting.

Prisoner's dilemma becomes more interesting in the iterated version of the game, where the same two players play repeatedly. This allows players to "punish" each other for uncooperative play. Scoring generally works as so (higher is better),

		Player A
		coop	defect
Player B	coop	(3,3)	(0,5)
Player B	defect	(5,0)	(1,1)

The most famous, and strongest individual strategy, is tit-for-tat. This player begins by playing cooperatively, then does whatever the its opponent did last. Here is the common lisp code to run a tit-for-tat strategy,

(defun tit-for-tat ()
  (lambda (x)
    (if (null x) :coop x)))

If you are unfamiliar with common lisp, the lambda part is returning an anonymous function that actually plays the tit-for-tat strategy. The tit-for-tat function generates a tit-for-tat player along with its own closure. The argument to the anonymous function supplies the opponent's last move, which is one of the symbols :coop or :defect. In the case of the first move, nil is passed. These are some really simple strategies that ignore their arguments,

(defun rand-play ()
  (lambda (x)
    (declare (ignore x))
    (if (> (random 2) 0) :coop :defect)))

(defun switcher-coop ()
  (let ((last :coop))
    (lambda (x)
      (declare (ignore x))
      (if (eq last :coop) 
	  (setf last :defect) 
	  (setf last :coop)))))

(defun switcher-defect ()
  (let ((last :defect))
    (lambda (x)
      (declare (ignore x))
      (if (eq last :coop) 
	  (setf last :defect) 
	  (setf last :coop)))))

(defun always-coop ()
  (lambda (x) 
    (declare (ignore x))
    :coop))

(defun always-defect ()
  (lambda (x) 
    (declare (ignore x))
    :defect))

Patrick Grim did an interesting study about ten years ago on iterated prisoner's dilemma involving competing strategies in a 2-dimensional area: Undecidability in the Spatialized Prisoner's Dilemma: Some Philosophical Implications. It is very interesting, but I really wanted to play around with some different configurations myself. So what I did was extend my iterated prisoner's dilemma engine above to run over a 2-dimensional grid.

Grim's idea was this: place different strategies in a 2-dimensional grid. Each strategy competes against its immediate neighbors. (The paper doesn't specify which kind of neighbor, 4-connected or 8-connected, so I went with 4-connected.) The sum of these competitions are added up to make that cell's final score. After scoring, each cell takes on the strategy of its highest neighbor, if any of its neighbors have a higher score than itself. Repeat.

The paper showed some interesting results, where the tit-for-tat strategy would sometimes dominate, and, in other cases, be quickly wiped out, depending on starting conditions. Here was my first real test of my simulation. Three strategies were placed randomly in a 50x50 grid: tit-for-tat, always-cooperate, and always-defect. This is the first twenty iterations. It stabilizes after 16 iterations.

(run-random-matrix 50 100 20 '(tit-for-tat always-coop always-defect))

Animated GIF

White is always-cooperate, black is always-defect, and cyan is tit-for-tat. Notice how the always-defect quickly exploits the always-cooperate and dominates the first few iterations. However, as the always-cooperate resource becomes exhausted, the tit-for-tat cooperative strategy works together with itself, as well as the remaining always-cooperate, to eliminate the always-defect invaders, who have no one left to exploit. In the end, a few always-defect cells are left in equilibrium, feeding off of always-cooperate neighbors, who themselves have enough cooperating neighbors to hold their ground.

The effect can be seen more easily here. Around the outside is tit-for-tat, in the middle is always-cooperate, and a single always-defect cell is placed in the middle.

(run-matrix (create-three-box) 100 30)

Animated GIF

The asymmetric pattern is due to the way that ties are broken.

The lisp code only spits out text, which isn't very easy to follow whats going on. To generate these gifs, I first used this Octave script to convert the text into images. Just dump the lisp output to a text file and remove the hash table dump at the end. Then run this script on that file: pd_plot.m. The text file input should look like this: example.txt. You will need Octave-Forge

The script will make PNGs. You can either change the script to make GIFs (didn't try this myself), or use something like ImageMagick to convert the images afterward. Then, you compile frames into the animated GIF using Gifsicle.

See if you can come up with some different strategies and make some special patterns for them. You may be able to observe some interesting interactions. The image at the beginning of the article uses all of the listed strategies in a random matrix.

I will continue to try out some more to see if I can find something particularly interesting.

[] permanent link

Don't stop here! This isn't everything. Check out the archives (on the left) for more posts. Or just have a look at the index.