// Genetic Algorithm for nonlinear programming
// Written by Microsoft Visual C++
// Copyright by UTLab @ Tsinghua University
// http://orsc.edu.cn/UTLab

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "UTLab.h"

static void  initialization(void);
static void  evaluation(int gen);
static void  selection(void);
static void  crossover(void);
static void  mutation(void);
static void  objective_function(void);
static int   constraint_check(double x[]);

#define N 3  // number of variables
#define M 1  // number of objectives
#define TYPE 1 // 1=max;-1=min
#define GEN 400 // maximum generation number
#define POP_SIZE 30
#define P_MUTATION 0.2
#define P_CROSSOVER 0.3

double  CHROMOSOME[POP_SIZE+1][N+1];
double  OBJECTIVE[POP_SIZE+1][M+1];
double  q[POP_SIZE+1];

static void objective_function(void)
{
	double x1,x2,x3;
	int i;
	for(i = 1; i <= POP_SIZE; i++) {
		x1 = CHROMOSOME[i][1];
		x2 = CHROMOSOME[i][2];
		x3 = CHROMOSOME[i][3];
		OBJECTIVE[i][1] = sqrt(x1)+sqrt(x2)+sqrt(x3);
	}
	for(i=1;i<=POP_SIZE;i++)
	  OBJECTIVE[i][0]= OBJECTIVE[i][1];
}

static int constraint_check(double x[])
{
	double a;
	int n; 
	for(n=1;n<=N;n++) if(x[n]<0) return 0;
	a = x[1]*x[1]+2*x[2]*x[2]+3*x[3]*x[3];
	if(a>1) return 0;
	return 1;
}

static void initialization(void)
{
  double x[N+1]; // N is the number of variables
  int i,j;
  for(i=1; i<=POP_SIZE; i++){
	  mark:
	  for(j=1; j<=N; j++) x[j]=myu(0,1);
	  if(constraint_check(x)==0) goto mark;
	  for(j=1; j<=N; j++) CHROMOSOME[i][j]=x[j];
  }
}

main()
{
  int i, j;
  double a;

  q[0]=0.05; a=0.05;
  for(i=1; i<=POP_SIZE; i++) {a=a*0.95; q[i]=q[i-1]+a;}
  initialization();
  evaluation(0);
  for(i=1; i<=GEN; i++) {
	  selection();
	  crossover();
	  mutation();
	  evaluation(i);
	  printf("\nGeneration NO.%d\n", i);
	  printf("x=(");
	  for(j=1; j<=N; j++) {
		  if(j<N) printf("%3.4f,",CHROMOSOME[0][j]);
		  else printf("%3.4f",CHROMOSOME[0][j]);
	  }
	  if(M==1) printf(")\nf=%3.4f\n", OBJECTIVE[0][1]);
	  else {
	      printf(")\nf=(");
	      for(j=1; j<=M; j++) {
		     if(j<M) printf("%3.4f,", OBJECTIVE[0][j]);
		     else printf("%3.4f", OBJECTIVE[0][j]);
		  }
          printf(")  Aggregating Value=%3.4f\n",OBJECTIVE[0][0]);
	  }
  }
  printf("\n");
  return 1;
}

static void evaluation(int gen)
{
  double a;
  int   i, j, k, label;
  objective_function();
  if(gen==0){
	 for(k=0; k<=M; k++) OBJECTIVE[0][k]=OBJECTIVE[1][k];
	 for(j = 1; j <= N; j++) CHROMOSOME[0][j]=CHROMOSOME[1][j];
  }
  for(i=0; i<POP_SIZE; i++){
	  label=0;  a=OBJECTIVE[i][0];
	  for(j=i+1; j<=POP_SIZE; j++)
		 if((TYPE*a)<(TYPE*OBJECTIVE[j][0])) {
			 a=OBJECTIVE[j][0];
			 label=j;
		 }
	  if(label!=0) {
		 for(k=0; k<=M; k++) {
			 a=OBJECTIVE[i][k];
			 OBJECTIVE[i][k]=OBJECTIVE[label][k];
			 OBJECTIVE[label][k]=a;
		 }
		 for(j=1; j<=N; j++) {
			 a=CHROMOSOME[i][j];
			 CHROMOSOME[i][j]=CHROMOSOME[label][j];
			 CHROMOSOME[label][j]=a;
		 }
	  }
  }
}

static void selection()
{
  double r, temp[POP_SIZE+1][N+1];
  int   i, j, k;
  for(i=1; i<=POP_SIZE; i++) {
	  r=myu(0, q[POP_SIZE]);
	  for(j=0; j<=POP_SIZE; j++) {
		  if(r<=q[j]) {
			  for(k=1; k<=N; k++) temp[i][k]=CHROMOSOME[j][k];
			  break;
		  }
	  }
  }
  for(i=1; i<=POP_SIZE; i++)
	 for(k=1; k<=N; k++)
		 CHROMOSOME[i][k]=temp[i][k];
}

static void crossover()
{
  int   i, j, jj, k, pop;
  double r, x[N+1], y[N+1];
  pop=POP_SIZE/2;
  for(i=1; i<=pop; i++) {
	 if(myu(0,1)>P_CROSSOVER) continue;
	 j=(int)myu(1,POP_SIZE);
	 jj=(int)myu(1,POP_SIZE);
	 r=myu(0,1);
	 for(k=1; k<=N; k++) {
		 x[k]=r*CHROMOSOME[j][k]+(1-r)*CHROMOSOME[jj][k];
		 y[k]=r*CHROMOSOME[jj][k]+(1-r)*CHROMOSOME[j][k];
	 }
	 if(constraint_check(x)==1)
		 for(k=1; k<=N; k++) CHROMOSOME[j][k]=x[k];
	 if(constraint_check(y)==1)
		 for(k=1; k<=N; k++) CHROMOSOME[jj][k]=y[k];
  }
}

static void mutation(void)
{
  int i, j, k;
  double x[N+1], y[N+1], infty, direction[N+1];
  double INFTY=10, precision=0.0001;
  for(i=1; i<=POP_SIZE; i++) {
	  if(myu(0,1)>P_MUTATION) continue;
	  for(k=1; k<=N; k++) x[k] = CHROMOSOME[i][k];
	  for(k=1; k<=N; k++)
		  if(myu(0,1)<0.5) direction[k]=myu(-1,1);
		  else direction[k]=0;
	  infty=myu(0,INFTY);
	  while(infty>precision) {
		  for(j=1; j<=N; j++) y[j]=x[j]+infty*direction[j];
		  if(constraint_check(y)==1) {
			 for(k=1; k<=N; k++) CHROMOSOME[i][k]=y[k];
			 break;
		  }
		  infty=myu(0,infty);
	  }
  }
}