// Genetic Algorithm for Principle of Least Squares
// 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 2  // 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];

double DATA[6][2]={
0.6, 0.1,
1.0,0.3,
1.5, 0.4,
2.0, 0.6,
2.8, 0.8,
3.6,0.9,
};

static void objective_function(void)
{
	double x1,x2, t;
	int i, j;
	for(i = 1; i <= POP_SIZE; i++) {
		x1 = CHROMOSOME[i][1];
		x2 = CHROMOSOME[i][2];
        OBJECTIVE[i][1] = 0;
		for(j=0; j<=5; j++){
            t=1+exp(3.14159*(x1-log(DATA[j][0]))/(sqrt(3)*x2));
            t=1/t;
			t=(t-DATA[j][1])*(t-DATA[j][1]);
	    	OBJECTIVE[i][1] = OBJECTIVE[i][1]+t;
		}
	}
	for(i=1;i<=POP_SIZE;i++)
	  OBJECTIVE[i][0]= OBJECTIVE[i][1];
}

static int constraint_check(double x[])
{
	if(x[1]>5||x[1]<0||x[2]>5||x[2]<0) 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,5);
	  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]);
	  }
	  printf(")\nLeast Squares = %3.4f\n", OBJECTIVE[0][1]);
 }
  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);
	  }
  }
}