1. Genetic Algorithm problem

Hello guys i have coded a GA which optimizes airfoil shapes. The problem i have is that the code runs fine untill generation 50 and then on generation 51 it crashes..! Actually the code is still running but is produces zeros and some infinite numbers. Can you please help me out to find out what am i doing wrong..??

The code works like this...

Aifoils are produced from x-y coordinates. I have the coordinates on the x-axis fixed and each new shape is produced by altering the values on the y-axis. The code uses Xfoil as an evaluation tool.

These are the files needed for the code to run

Gadata stores the limits of each control point. in a form lower upper
ControlPoints2.txt Contains the control points coordinate
airfoilcoordinates.dat Contains the airfoil's coordinates
airfoil.txt Contains the commands for the Xfoil
xfoil.def Contains the parameters for the Xfoil
Ga_log.txt Contains the results of the algorithm
run.bat Isrfoil.txt and execute the commands
pplot.exe and pxplot.exe are needed in order for Xfoil to plot the airfoils they come together with Xfoil when you download it.. you can download it from here http://web.mit.edu/drela/Public/web/xfoil/!

The coordinates for the spline are on a form of
x1
y1
x2
y2 etc

In order for the code to run you have to change the extension of xfoil.txt to xfoil.def and run.txt to run.bat i have renamed them in order to be able to upload them.

Thank you again and i would be happy to answer any questions you may have about the code and genetic algorithms and evolution strategies...!!!

2. Just so we're all clear on this, we can comment out the calls to run external processes and just watch the code you posted go through it's iterations.

3. OK, but Xfoil is used in order to calculate the value of fitness function...!

4. Meh - doesn't look like it.

In non_uniform_mutation() and crossover() (possibly others)
> double x = ((rand() % 1000)/1000);
> double r = (rand()%1); /* r = (randval(0,1) >= 0.5); */
These are always 0. The calculations are done in integer arithmetic, and then only at the last moment, converted to doubles.
You have a randval() function, perhaps you should use it.

Are you attempting to just grab the last line of the file?
A while loop of
Code:
```while ( fgets(buff,sizeof(buff),fd) != NULL ) {
}
if ( feof(fd) ) {  // note: if ferror() returns true, buff is garbage.
if ( sscanf(buff, "%lf %lf %lf",&n,&cl,&cd) == 3 ) {
printf("cl: %1.5f  cd: %1.5f\n",cl,cd);
}
}```
would suffice.
Also, calling something which is a FILE* an 'fd' is a misnomer.
'fd' is typically used to denote file descriptors, which are integers, not FILE* pointers.

> float rbasis(int c,float u,int npts,int x[60],float h[28],float r[28])
28 has a #define constant (so use it), but what is 60?
The best I can tell is that it's a bit more than 28*2

> for (i = 1; i <= npts; i++)
There are several places where you seem to assume that arrays start at 1.

> for (ispline = 1; ispline <= 2*npts; ispline=ispline+2)
This might explain why some of your arrays are 60, as being a bit more than 28*2, to take into account all of this 1-based array stuff.
Except SOME of your 1-based arrays are not padded at the end, and this will cause grief.

5. Salem Thanks a lot for reviewing my code...! I took your comments into consideration and i did this...

Actually my randval values doesn't produce so random values so i did some changes into the functions of crossover and non_uniform_mutation() now my functions are like this...

Code:
```void non_uniform_mutation(void)
{
int i,j;
double lbound,hbound;
double r,x;

for(i = 0; i < POPSIZE; i++)
{
for(j = 0; j < NVARS; j++)
{
x = (float)(rand() % 1000)/1000;
if(x < PMUTATION)
{
lbound = population[i].lower[j];
hbound = population[i].upper[j];

population[i].gene[j] = randval(lbound,hbound);

r = rand() % 2;        /* r = (randval(0,1) >= 0.5); */

if(r == 0)
{
population[i].gene[j] = population[i].gene[j] + delta2(generation,hbound-population[i].gene[j]);

}

if(r == 1)
{
population[i].gene[j] = population[i].gene[j] - delta2(generation,population[i].gene[j]-lbound);

}

}
}
}
}```
Code:
```void crossover(void)
{
int mem,one=0;
int first = 0;
double x;

for(mem =0; mem < POPSIZE; mem++)
{
x = (float)(rand() % 1000)/1000;

if(x < PXOVER)
{
first++;

if(first%2 == 0)
{
random_Xover(one,mem);
}

else
{
one = mem;
}
}
}
}```
Now as far as the float rbasis(int c,float u,int npts,int x[60],float h[28],float r[28])

r is an array which contains the values of the basis function of a b-spline. Yes it has a constant value but it is by change if you change the degree of the spline then you should change the size of r[] as more basis functions will be needed.

Yes my arrays start from one because they represent curve points in a form of (x1,y1),(x2,y2) etc....

Actually you are right but i usually assign values at the beginning to the arrays for example

Code:
```for (ispline=1; ispline <= npts; ispline++)
{
h[ispline] = 1.0;
}```

Code:
```for (ispline = 1; ispline <= 2*p1; ispline++)
{
p[ispline] = 0.;
}

nplusc = npts + k;

/*  zero and redimension the knot vector and the basis array */

for(i = 0; i <= npts; i++){
nbasis[i] = 0.;
}

for(i = 0; i <= nplusc; i++){
x[i] = 0.;
}```
Unfortunately i still don't know why it produces infinite numbers after generation 50...

6. > for(i = 0; i <= npts; i++)
Well you see, stuff like this is a red rag to a bull when it comes to C programming, because it screams array overrun.

If you've got an array[N], then having <= N in your loop test is always a bad thing.

7. Hello Salem..!!

I changed my for loops as you proposed but now my spline function is calculating nan values..!!

This is my updated code can you please tell me why is that happenning cause i cannot figure it out...!!

Code:
```int knot(int n,int c,int x[60])

{
int nplusc,nplus2,i;

nplusc = n + c;
nplus2 = n + 2;

x[0] = 0;
for (i = 0; i < nplusc; i++){
if ( (i > c) && (i < nplus2) )
x[i+1] = x[i] + 1;
else
x[i] = x[i+1];
}

return 0;
}

float rbasis(int c,float u,int npts,int x[60],float h[28],float r[30])

{
int nplusc;
int i,k;
float d,e;
float sum;
float temp[36];

nplusc = npts + c;

/*        printf("knot vector is \n");
for (i = 1; i <= nplusc; i++){
printf(" %d %d \n", i,x[i]);
}
printf("u is %lf \n", u);
*/

/* calculate the first order nonrational basis functions n[i]    */

for (i = 0; i < nplusc; i++){
if (( u >= x[i]) && (u < x[i+1]))
temp[i] = 1;
else
temp[i] = 0;
}

/* calculate the higher order nonrational basis functions */

for (k = 0; k < c+1; k++){
for (i = 0; i < nplusc-k+1; i++){
if (temp[i] != 0) {   /* if the lower order basis function is zero skip the calculation */
d = ((u-x[i])*temp[i])/(x[i+k]-x[i]);}
else{
d = 0;}

if (temp[i+1] != 0) {    /* if the lower order basis function is zero skip the calculation */
e = ((x[i+k]-u)*temp[i+1])/(x[i+k]-x[i]);}
else {
e = 0;}   /*line for check */

temp[i] = d + e;
}
}

if (u == (float)x[nplusc]){        /*    pick up last point    */
temp[npts] = 1;
}
/*
printf("Nonrational basis functions are \n");
for (i=1; i<= npts; i++){
printf("%lf ", temp[i]);
}
printf("\n");
*/
/* calculate sum for denominator of rational basis functions */

sum = 0.;
for (i = 0; i < npts; i++){
sum = sum + temp[i]*h[i];
}

/* form rational basis functions and put in r vector */

for (i = 0; i < npts; i++){
if (sum != 0){
r[i] = (temp[i]*h[i])/(sum);}
else
r[i] = 0;
}

return 0;
}

float rbspline()

{
int i,j,icount,jcount;
int i1;
int x[60];        /* data points store vector */
int nplusc;
int ispline;
int npts,k;
int p1;

float step;
float u;
float nbasis[30];
float temp;

float b[60];  /* control vertices matrix*/
float p[203]; /* allows for up to 200 points on curve */
float h[28];

FILE *splineinput;
FILE *airfoilcoordinates;

splineinput = fopen("controlpoints2.txt","r");
airfoilcoordinates = fopen("airfoilcoordinates.dat","w");

npts = 28;
k = 6;     /* third ordech=getc();r, change for other orders */
p1 = 100;   /* eleven points on curve */

for (ispline = 0; ispline < 2*npts; ispline++)
{
fscanf(splineinput,"%f\n ",&b[ispline]);
}

printf("\nPolygon points\n");

for (ispline = 0; ispline < 2*npts; ispline=ispline+2)
{
printf(" %1.5f %1.5f  \n",b[ispline],b[ispline+1]);
}

for (ispline=0; ispline < npts; ispline++)
{
h[ispline] = 1.0;
}

for (ispline = 0; ispline < 2*p1; ispline++)
{
p[ispline] = 0.;
}

nplusc = npts + k;

/*  zero and redimension the knot vector and the basis array */

for(i = 0; i < npts; i++){
nbasis[i] = 0.;
}

for(i = 0; i < nplusc; i++){
x[i] = 0.;
}

/* generate the uniform open knot vector */

knot(npts,k,x);

/*
printf("The knot vector is ");
for (i = 1; i <= nplusc; i++){
printf(" %d ", x[i]);
}
printf("\n");
*/

icount = 0;

/*    calculate the points on the rational B-spline curve */

u = 0;
step = ((float)x[nplusc])/((float)(p1-1));

for (i1 = 0; i1 < p1; i1++){

if ((float)x[nplusc] - u < 5e-6){
u = (float)x[nplusc];
}

rbasis(k,u,npts,x,h,nbasis);      /* generate the basis function for this value of t */
/*
printf("u = %f \n",u);
printf("nbasis = ");
for (i = 1; i <= npts; i++){
printf("%f  ",nbasis[i]);
}
printf("\n");
*/
for (j = 0; j < 2; j++){      /* generate a point on the curve */
jcount = j;
p[icount+j] = 0.;

for (i = 0; i < npts; i++){ /* Do local matrix multiplication */
temp = nbasis[i]*b[jcount];
p[icount + j] = p[icount + j] + temp;
/*
printf("jcount,nbasis,b,nbasis*b,p = %d %f %f %f %f\n",jcount,nbasis[i],b[jcount],temp,p[icount+j]);
*/
jcount = jcount + 2;
}
}
/*
printf("icount, p %d %f %f %l \n",icount,p[icount+1],p[icount+2],p[icount+3]);
*/
icount = icount + 2;
u = u + step;
}

printf("\nHomogeneous weighting vector is \n");
for (ispline = 0; ispline < npts; ispline++){
printf(" %1.4f ", h[ispline]);
}

printf("\n");

/*printf("\nCurve points\n"); */

for (ispline = 0; ispline < 2*p1; ispline=ispline+2)
{
fprintf(airfoilcoordinates,"%1.5f %1.5f \n",p[ispline],p[ispline+1]);
}

fclose(splineinput);
fclose(airfoilcoordinates);

return 0;
}```

8. In knot(), you changed
Code:
```	x[1] = 0;
for (i = 2; i <= nplusc; i++){
if ( (i > c) && (i < nplus2) )
x[i] = x[i-1] + 1;

to
x[0] = 0;
for (i = 0; i < nplusc; i++) {
if ((i > c) && (i < nplus2))
x[i + 1] = x[i] + 1;```
Why have you changed things over from x[i-1] to x[i+1] ?

And
Code:
```/* calculate the higher order nonrational basis functions */
for (k = 2; k <= c; k++){
for (i = 1; i <= nplusc-k; i++){
if (temp[i] != 0) {   /* if the lower order basis function is zero skip the calculation */
became
/* calculate the higher order nonrational basis functions */
for (k = 0; k < c + 1; k++) {
for (i = 0; i < nplusc - k + 1; i++) {```
Sliding the array down one index would have made
for (k = 2; k <= c; k++)
into
for (k = 1; k < c; k++)

And
for (i = 1; i <= nplusc-k; i++)
into
for (i = 0; i < nplusc-k; i++)

It's basically
- subtract 1 from the start condition
- turn the comparison from <= into <

Rewriting the code and then complaining it doesn't work doesn't help.

I would tell you to learn how to use a debugger, then you can put breakpoints in the code, single step the execution, examine your variables until you home in on where the calculation is going wrong.
Which will typically be down to either uninitialised variables or array overruns.

If you want to stick with your 1-based arrays, which did at least seem to work until the last iteration, then perhaps the easier alternative is to just make the arrays a bit longer than necessary.
It wastes a bit of memory, but it might be an easier code transformation than the one you've basically managed to mess up.

9. Salem thank you for your support...!! Since i believe that the problem might lie within the arrays of the rbspline function maybe... I wanted to ask you is there a way to dynamically change the size of these arrays...?