// For Project Investment Optimization Problem in Stochastic Environment (EVM)
// Written by Microsoft Visual C++
// Copyright by UTLab at Tsinghua University

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <time.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  f(void);
static int   constraint_check(int x[]);

#define N 7  // number of variables
#define M 2  // number of objectives
#define   TYPE -1 // 1=max;-1=min
#define   GEN 3000 // maximum generation number
#define   POP_SIZE 60
#define r 0.006   // interest rate per month
double P_MUTATION=0.2;
double P_CROSSOVER=0.3;
double OBJECTIVE[POP_SIZE+1][M+1], q[POP_SIZE+1];
int CHROMOSOME[POP_SIZE+1][N+1];
static void   rearrange(int x[]);

static void rearrange(int x[]) // to rearrange the schedule x!!
{
double temp[4];

//Path 1
temp[1]=x[1];
temp[2]=x[2];
temp[3]=x[5];
sortminmax(temp,1,3);
x[1]=temp[1];
x[2]=temp[2];
x[5]=temp[3];

//Path 2
temp[1]=x[1];
temp[2]=x[3];
temp[3]=x[5];
sortminmax(temp,1,3);
x[1]=temp[1];
x[3]=temp[2];
x[5]=temp[3];

//Path 3
temp[1]=x[1];
temp[2]=x[3];
temp[3]=x[6];
sortminmax(temp,1,3);
x[1]=temp[1];
x[3]=temp[2];
x[6]=temp[3];

//Path 4
temp[1]=x[1];
temp[2]=x[3];
temp[3]=x[7];
sortminmax(temp,1,3);
x[1]=temp[1];
x[3]=temp[2];
x[7]=temp[3];

//Path 5
temp[1]=x[1];
temp[2]=x[4];
temp[3]=x[7];
sortminmax(temp,1,3);
x[1]=temp[1];
x[4]=temp[2];
x[7]=temp[3];
}


static void f(void)
{
    int i,j;
	double sum,xi[12],P[6],timesum;
	static double c[12]= {0,1500,1800,430,1600,800,500,2000,2100,550,530,630};

	for(i=1;i<=POP_SIZE;i++)
	{
		sum=0;
		timesum=0;
		for(j=1;j<=5000;j++) 
		{
			P[1]=CHROMOSOME[i][1];
			P[2]=CHROMOSOME[i][1];
			P[3]=CHROMOSOME[i][1];
			P[4]=CHROMOSOME[i][1];
			P[5]=CHROMOSOME[i][1];
			
			xi[1]=myn(9,3);
			
			xi[2]=myn(5,2);
						
			xi[3]=myn(10,3);
            
			xi[4]=myn(6,2);
			
			xi[5]=myn(8,2);
						
			xi[6]=myn(8,2);
			
			xi[7]=myn(9,3);
			
			xi[8]=myn(6,1);
			
			xi[9]=myn(10,2);
			
			xi[10]=myn(15,3);
			
			xi[11]=myn(11,2);
			//Compute the length of the eight paths
			P[1]+=xi[1]; 
			if(P[1]<CHROMOSOME[i][2]) P[1]=CHROMOSOME[i][2];
            P[1]+=xi[4]; 
			if(P[1]<CHROMOSOME[i][5]) P[1]=CHROMOSOME[i][5];
            
            P[2]+=xi[2]; 
			if(P[2]<CHROMOSOME[i][3]) P[2]=CHROMOSOME[i][3];
                        P[3]=P[2];
                        P[4]=P[2];
            P[2]+=xi[5]; 
			if(P[1]<P[2]) P[1]=P[2];

             P[1]+=xi[9]; 
			
			// P 1 represents P 1,2
			
	    			
            P[3]+=xi[6]; 
			if(P[3]<CHROMOSOME[i][6]) P[3]=CHROMOSOME[i][6];
            P[3]+=xi[10]; 

			// P 3 represents P 3

            P[4]+=xi[7]; 
			if(P[4]<CHROMOSOME[i][7]) P[4]=CHROMOSOME[i][7];
                        
                        
            P[5]+=xi[3];
			if(P[5]<CHROMOSOME[i][4]) P[5]=CHROMOSOME[i][4];
			P[5]+=xi[8];
			
            if(P[4]<P[5]) P[4]=P[5];
            P[4]+=xi[11];
			//P 4 represents P 4,5			
			
                     

                        if(P[1]<P[3]) P[1]=P[3];
						if(P[1]<P[4]) P[1]=P[4];
                        

                        if(P[1]-int(P[1])==0)            
			P[0]=P[1];
						else P[0]=int(P[1])+1;

           timesum+=P[0];
			
			sum+=c[1]*pow(1+r,ceil(P[0]-CHROMOSOME[i][1]));
			sum+=c[2]*pow(1+r,ceil(P[0]-CHROMOSOME[i][1]));
			sum+=c[3]*pow(1+r,ceil(P[0]-CHROMOSOME[i][1]));
			sum+=c[4]*pow(1+r,ceil(P[0]-CHROMOSOME[i][2]));
			sum+=c[5]*pow(1+r,ceil(P[0]-CHROMOSOME[i][3]));
			sum+=c[6]*pow(1+r,ceil(P[0]-CHROMOSOME[i][3]));
			sum+=c[7]*pow(1+r,ceil(P[0]-CHROMOSOME[i][3]));
			sum+=c[8]*pow(1+r,ceil(P[0]-CHROMOSOME[i][4]));
			sum+=c[9]*pow(1+r,ceil(P[0]-CHROMOSOME[i][5]));
			sum+=c[10]*pow(1+r,ceil(P[0]-CHROMOSOME[i][6]));
			sum+=c[11]*pow(1+r,ceil(P[0]-CHROMOSOME[i][7]));
			
		}
		OBJECTIVE[i][1]=sum/5000;
        OBJECTIVE[i][2]=timesum/5000.;
	}
    for(i=1;i<=POP_SIZE;i++) OBJECTIVE[i][0]= OBJECTIVE[i][1];   
}



static int constraint_check(int x[])
{
	double P[6],sum,xi[12];
    int i,j;
	sum=0;
	for(i=1;i<=N;i++) {if(x[i]<0) return 0;}
	for(j=1;j<=5000;j++)
	{
		P[1]=x[1];
		P[2]=x[1];
		P[3]=x[1];
		P[4]=x[1];
		P[5]=x[1];

	        xi[1]=myn(9,3);
			
			xi[2]=myn(5,2);
						
			xi[3]=myn(10,3);
            
			xi[4]=myn(6,2);
			
			xi[5]=myn(8,2);
						
			xi[6]=myn(8,2);
			
			xi[7]=myn(9,3);
			
			xi[8]=myn(6,1);
			
			xi[9]=myn(10,2);
			
			xi[10]=myn(15,3);
			
			xi[11]=myn(11,2);
			
                
            P[1]+=xi[1]; 
			if(P[1]<x[2]) P[1]=x[2];
            P[1]+=xi[4]; 
			if(P[1]<x[5]) P[1]=x[5];
            
            P[2]+=xi[2]; 
			if(P[2]<x[3]) P[2]=x[3];
                        P[3]=P[2];
                        P[4]=P[2];
            P[2]+=xi[5]; 
			if(P[1]<P[2]) P[1]=P[2];

             P[1]+=xi[9]; 
			
			// P 1 represents P 1,2
			
	    			
            P[3]+=xi[6]; 
			if(P[3]<x[6]) P[3]=x[6];
            P[3]+=xi[10]; 

			// P 3 represents P 3

            P[4]+=xi[7]; 
			if(P[4]<x[7]) P[4]=x[7];
                        
                        
            P[5]+=xi[3];
			if(P[5]<x[4]) P[5]=x[4];
			P[5]+=xi[8];
			
            if(P[4]<P[5]) P[4]=P[5];
            P[4]+=xi[11];
			//P 4 represents P 4,5			
			
                     

                        if(P[1]<P[3]) P[1]=P[3];
						if(P[1]<P[4]) P[1]=P[4];
                        

                        if(P[1]-int(P[1])==0)            
			P[0]=P[1];
						else P[0]=int(P[1])+1;
		sum+=P[0];
	}
	sum=sum/5000.;
	
	sum=sum-x[1];
	if(sum>32) return 0; 
    return 1;
}

static void initialization(void)
{
  int x[N+1];
  int i,j;
  for(i=1; i<=POP_SIZE; i++)
  {
      mark:  //****************************************************
	  x[7]=int(myu(0,20));
	  x[6]=int(myu(0,20));
	  x[5]=int(myu(0,20));
	  x[4]=int(myu(0,20));
	  x[3]=int(myu(0,20));
	  x[2]=int(myu(0,20));	  
	  x[1]=int(myu(0,20));

          rearrange(x); // rearrange the schedule x
      if(constraint_check(x)==0) goto mark;
	  for(j=1; j<=N; j++) CHROMOSOME[i][j]=x[j];
	  printf("N=%d",i);
          
  }
}

main()
{
  int i, j;
  double a;
  FILE *fp;

  srand((unsigned)time(NULL));
  q[0]=0.05; a=0.05;
  for(i=1; i<=POP_SIZE; i++) {a=a*0.95; q[i]=q[i-1]+a;}
  fp=fopen("RESULT.dat","w");
  initialization();
  evaluation(0);
  for(i=1; i<=GEN; i++) {
      selection();
      crossover();
      mutation();
      evaluation(i);
           fprintf(fp,"No. %d,",i);
	  for(j=1; j<=N; j++) fprintf(fp,"%d,",CHROMOSOME[0][j]);
	  for(j=1; j<=M; j++) fprintf(fp,"%8.15f, ",OBJECTIVE[0][j]);
      fprintf(fp,"%8.15f\n",OBJECTIVE[0][0]);
      printf("\n Generation NO.%d\n", i);
      printf("x=(");
      for(j=1; j<=N; j++) {
          if(j<N) printf("%d,",CHROMOSOME[0][j]);
          else printf("%d",CHROMOSOME[0][j]);
           
      }
      if(M==1) printf(")\nf=%lf\n", OBJECTIVE[0][1]);
      else {
          printf(")\nf=(");
          for(j=1; j<=M; j++) {
             printf("%lf,", OBJECTIVE[0][j]);
          }
          printf(")  Aggregating Value=%lf\n",OBJECTIVE[0][0]);
      }
  }
  printf("\n");
  printf("\nThe final result has been written in RESULT.dat.\n\n");
   fclose(fp);
  return 1;
}

static void evaluation(int gen)
{
  int i, j, k, l,label,b;double a;
  f();
  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(l=1; l<=N; l++) {
             b=CHROMOSOME[i][l];
             CHROMOSOME[i][l]=CHROMOSOME[label][l];
             CHROMOSOME[label][l]=b;
         }
      }
  }
}

static void selection()
{
  double rr;
  int   i, j, k,temp[POP_SIZE+1][N+1];;
  for(i=1; i<=POP_SIZE; i++) {
      rr=myu(0, 0.785);
      for(j=0; j<=POP_SIZE; j++) {
          if(rr<=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;
  int x[N+1], y[N+1]; double rr;
  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);
     for(k=1; k<N; k++) x[k]=CHROMOSOME[j][k];
     for(k=1; k<N; k++) y[k]=CHROMOSOME[jj][k];
	 
		 rr=myu(0,1);
      
	 for(k = 1; k <= N; k++) {
		 x[k] = int(rr*CHROMOSOME[j][k]+(1-rr)*CHROMOSOME[jj][k]);
		 y[k] = int(rr*CHROMOSOME[jj][k]+(1-rr)*CHROMOSOME[j][k]);
	 }
     rearrange(x);
     if(constraint_check(x)==1)
         for(k=1; k<=N; k++) CHROMOSOME[j][k]=x[k];
		 rearrange(y);
     if(constraint_check(y)==1)
         for(k=1; k<=N; k++) CHROMOSOME[jj][k]=y[k];
  }
}

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