Month: September 2012

Simulated Annealing

Recently, I have been fascinated by algorithms generally used for AI. Simulated Annealing is one of the basic algorithms for finding solutions a constrained problem. I have been greatly helped along this path by a book : Artificial Intelligence Application Programming(1st edition) by M. Tim Jones.

Annealing is a process in which a metal is first heated to a high temperature and then slowly cooled down resulting in formation of large, ordered crystal structures. The output of this process is dependent on the time taken for the cooling. Faster cooling times result in formation of brittle, fragile structures while slow cooling results in the formation of large stable crystals.

Simulated annealing as the name suggests, is creating an algorithm which mimics this cooling process. Basically, several problems of varied nature can be solved by creating an analogy with the annealing process. Many problems can thus be equated to the annealing process. For example, given a minimization problem, we can equate the ‘Energy’ of the solution to the value of the objective function and then randomly tweak the values within the given constraints to find the minimum solution.

The algorithm for simulated annealing is as follows:

1. Create an initial solution.

2. Assess it’s value

3. Randomly tweak the solution.

4. Assess the new solution

5.Check acceptance criteria

6. Reduce Temperature

7. Go to step 3

When the temperature reaches a certain predetermined value, the current solution which satisfies all criteria is outputted as the best solution.

To model this, let us assume that the lower the ‘energy’ of a solution, the better it is. Now for the algorithm to converge on a solution, we must propagate the solutions having lower values of energy. Thus we use the acceptance criteria for that. So, if the new solution is lower in energy than the current solution, we move on to reducing the temperature. But, if the working solution is worse, then we have an acceptance criteria which we can follow. In simulated annealing, we use an acceptance criteria based on the laws of thermodynamics.

The following example is a basic example of simulated annealing.

N-Queens problem:

A very famous problem, it was first solved by Carl Friedrich Gauss in 1850 by trial and error. Since then, several methods have been used including divide-and-conquer, depth first, genetic algorithms etc. The problem in it’s conception is very simple. We have an N x N chess board and we have N queens. Now, we have to place the queens such that there is no conflict between any two queens as per the familiar rules of chess. As a queen can move diagonally, vertically and horizontally, we have to ensure that for an acceptable solution, there must be only one queen in each column and row and it’s diagonals must not have any other queen.

The formulation of this problem can be done as follows:

Lets take an 1 x N matrix. The index of the matrix is the row index of the position of the queen while the value at that index is the column index of the queen. Thus we have defined the position of each of the N – queens. Now, we randomly allocate each of the queens on the board such that no 2 queens occupy the same row or column. This can be easily ensured by checking the condition that the value at each index is not equal to any before it.

We define the energy of each solution as the no. of conflicts i.e. the no. of queens threatening each other. We also define several macros which which help us in the coding process


 #define MAX_LENGTH 30
 #define INIT_TEMP 30.0
 #define FIN_TEMP 0.5
 #define ALPHA 0.98
 #define STEPS 100

MAX_LENGTH is the no. of queens i.e. N.

The temperature parameter basically defines the region in which the solution is being searched for. As the temperature decreases, the range of the solution space also decreases thus resulting in an optimum solution. The higher the temperature, the more time it will take for convergence but also there is a greater probability of finding a solution. Now, the final temperature cannot be zero as the value of the exponent then cannot be calculated. Thus, we try and take the value as close to zero as possible but at the cost of time. I tried several values and found out that at 0.5, the algorithm converges to a solution for even large N’s in a fairly decent amount of time.

ALPHA is the value of the constant in the linear function

T_{t+1} = \alpha T_{t}

This function can be any decreasing function of T including non linear functions. I used the above one because it is the easiest to analyse. Now, at each temperature, we allow for 100 iterations for a best solution to emerge. This solution is then propagated to the next temperature. This can be imagined as moving up or down a hill until we find a ditch which represents over optimum value. Then we can start following that ditch until we reach the actual global minimum. To ensure, we stay within the constraints of the ditch, we decrease the temperature and keep decreasing it with every 100 steps until at last we find a solution.

Thus we can see that STEPS must be sufficiently large while ALPHA must be as close to 1 as possible so as to generate the optimal solution. But we must also keep in mind the timin constraints for the program to find a solution.

Now let us encode the problem:

We create a struct for each solution in which we create an array representing the solution and the energy associated with it.

typedef struct {
solutionType solution;
float energy;
} memberType;

Now, we create several functions to help us with the program

void tweakSolution( memberType *member)
{
int temp, x, y;

x=rand() % MAX_LENGTH;

do
{
y=rand() % MAX_LENGTH;
}while(x == y );

temp=member->solution[x];
member->solution[x]=member->solution[y];
member->solution[y]=temp;

}

This is a simple function which will help us create a randomly perturbed solution.


void initSolution( memberType *member)
{
int i;

for(i=0;i<MAX_LENGTH;i++)
{
member->solution[i] = i;
}

for(i=0;i<MAX_LENGTH;i++)
{
tweakSolution(member);
}

}

This function initialises the first solution.

void computeEnergy(memberType *member)
{
int i,j,x,y,tempx,tempy;

char board[MAX_LENGTH][MAX_LENGTH];

int conflicts;

const int dx[4] = { -1,1,-1,1 };

const int dy[4] = { -1, 1, 1 , -1};

bzero((void*)board, MAX_LENGTH*MAX_LENGTH);

for(i=0;i<MAX_LENGTH; i++)
{
board[i][member->solution[i]] = 'Q';
}

conflicts = 0;
for(i=0;i<MAX_LENGTH;i++)
{
x=i;
y=member->solution[i];

for(j=0;j<4;j++)
{
tempx=x;
tempy=y;

while(1)
{
tempx+=dx[j];
tempy+=dy[j];

if(tempx<0 || tempx>=MAX_LENGTH || tempy<0 || tempy >= MAX_LENGTH)
{
break;
}

if(board[tempx][tempy] == 'Q' )
conflicts++;
}
}
}

member->energy = (float) conflicts;
}

This function is used to compute the energy of each solution. We define energy as the no. of conflicts. Since, we have ensured that while tweaking, we avoid placing 2 queens on the same row or column, we can easily reduce the complexity of finding the no. of conflicts to checking the diagonals of each queen.

There are a few more functions for just simplifying the coding process and emitting the putput for better readbility and understanding

void copySolution( memberType *dest, memberType *src)
{
int i;

for(i=0;i<MAX_LENGTH;i++)
{
dest->solution[i] = src->solution[i];
}

dest->energy = src->energy;
}

void emitSolution( memberType *member)
{
char board[MAX_LENGTH][MAX_LENGTH];

int x,y;
bzero((void*)board, MAX_LENGTH*MAX_LENGTH);

for(x=0;x<MAX_LENGTH;x++)
{
board[x][member->solution[x]] = 'Q';
}

printf("board: \n");

for(y=0;y<MAX_LENGTH;y++)
{
for(x=0;x<MAX_LENGTH;x++)
{
if(board[x][y] == 'Q')
printf("Q ");
else
printf(". ");
}
printf("\n");

}
printf("\n\n");

}

Now, we implement the actual simulated annealing algorithm in the main()

int main()
{
int timer = 0,step,solution=0, useNew, accepted;

float temp = INIT_TEMP;
memberType current,working,best;

FILE *fp;

fp=fopen("stats.txt", "w");

srand(time(NULL));

initSolution(&current);
computeEnergy(&current);

best.energy=100.0;

copySolution(&working, &current);

while( temp > FIN_TEMP )
{
printf("\n Temperature = %f", temp);

accepted = 0;

/* Monte Carlo step*/

for( step = 0; step < STEPS; step++);
{
useNew=0;

tweakSolution(&working);
computeEnergy(&working);

if(working.energy <= current.energy)
{
useNew = 1;
}
else
{
float test = rand() % 1;
float delta = working.energy - current.energy;
float calc = exp(-delta/temp);

if(calc > test)
{
accepted++;
useNew = 1;
}
}
}

if(useNew)
{
useNew = 0;
copySolution(&current, &working);

if(current.energy < best.energy)
{
copySolution(&best, &current);
solution = 1;
}

else
{
copySolution(&working, &current);
}

}

fprintf(fp, "%d %f %f %d \n", timer++, temp, best.energy, accepted);
printf("Best Energy = %f\n", best.energy);

temp *= ALPHA;
}

fclose(fp);

if(solution)
{
emitSolution(&best);
}

return 0;
}

In the main program, we have created two loops. The outer loop deals with the temperature and reduces the temperature after every STEPS steps. The inner loop on the other hand is an implementation of a Metropolis algorithm. I will surely discuss the Metropolis algorithm in detail in some other post. But, for now, what it basically does is that it picks up a tweaked solution and computes it’s energy. If the energy of the solution is lower than the energy of the current solution, then it replaces the current solution by the tweaked solution. But if the energy is higher,  it calculates the difference between the current and the tweaked energy, calculates it’s probability as according to Formula 1.1 and then compares it to a randomly generated threshold value.

The reason for doing this is that if the algorithm encounters a local optimum, it may get stuck in the same place thus causing the algorithm to fail. Due to the Metropolis Algorithm, we avoid this by allowing the solution to take values which may not be the best for that iteration.

The solution that we have thus far obtained is then compared to the best solution we have so far. If it’s energy is lower, we copy the new solution into the best one at the end of each temperature change.

The program will terminate when we reach our chosen final temperature thus giving us the optimal solution.

Now that we know ‘Simulated Annealing’, we can apply it to solve problems in various domains,

  • Path generation
  • Transportation problems
  • Image Reconstruction
  • Circuit layout

Here’s a screen-shot of the output.

Much of the source code is provided in the above mentioned book. The following post will proceed on with Ant Algorithms. Thanks for reading and please contact me for any questions.

Color Predicates

It’s been a long time since I posted something here. But here it is.

I started my current project with a simple objective- I wanted to build a completely natural human- computer interaction system. Basically, my final aim is to just have a few cameras and microphones work as the input devices instead of a mouse and keyboard. So, with this in mind, I started out with my first objective- to segment out and track the hand in a video stream.
The first problem that I encountered was that skin as such is not a single color but a heady mixture of reds, yellows and white.The problem is simplified a bit in the HSV colorspace, but not enough. A simple thresholding program will not work at all given such few constraints. So, while researching on how to solve this problem, I got a paper mentioning color predicates.

Now, what are color predicates? Color predicates, very simply put are a map in the colorspace. A simple region within which every color included is taken as being part of the colors to be segmented and outside of which are the rejected colors. So, I decided that this was the way to go. Took me a couple of days but I was able to build a simple color predicate training program.

My way of implementing color predicates:

I used a matrix of 255×255 to represent the Hue and Saturation values possible. I then provided several images of my hand cropped onto a white background ( Used GIMP – awesome piece of software ). Then I scanned each image pixel by pixel and rejecting white, I simply incremented the value at each H,S index of the matrix if I found a color other than white. Now using this simple algorithm, I trained a color predicate for around 20 images taken in various illumination schemes.  With a few manipulations for rejecting low values, I then used the color predicate for simple thresholding of images in my video stream.

            

The best part of this exercise was when I looked at the HSV map thus created in a text editor….amazingly, there were two different regions which could be clearly seen, corresponding to the red spectrum and the yellow spectrum.

The code for the above exercise:


//Color Predicates ver 1.2

#include<cv.h>
#include<highgui.h>
#include<stdio.h>

int colorp[256][256];
IplImage *img, *imghsv;
int h,s,v;
int main()
{

int i,j;
for(i=0;i<256;i++)
{
for(j=0;j<256;j++)
{

colorp[i][j]=0;
}
}
FILE *fp = fopen("cp2.bin", "r");
for(i=0;i<256;i++)
fread(colorp[i], sizeof(int), 256, fp);

CvScalar Data;
img=cvLoadImage("train8.jpg");
imghsv=cvCreateImage(cvGetSize(img),8,3);
cvCvtColor(img,imghsv,CV_BGR2HSV);
cvNamedWindow("Image");
cvNamedWindow("HSV");

for(i=0;i<img->height;i++)
{
for(j=0;j<img->width;j++)
{
Data=cvGet2D(imghsv,i,j);
h=Data.val[0];
s=Data.val[1];
v=Data.val[2];
if(s!=0)
{
colorp[h][s]+=1;
}

}
}
fclose(fp);
FILE* fp2=fopen("cp2.bin","w");
for(i=0;i<256;i++)
{
fwrite(colorp[i], sizeof(int), 256*256, fp2);
}

fclose(fp2);

return 0;

}

In the next post, I will be describing my attempt to use the created color predicate.