// Stochastic Dependent-Chance Programming
// Project Investment Optimization Problem 
// Written by Microsoft Visual C++
// Copyright by UTLab @ Tsinghua University
// http://orsc.edu.cn/UTLab

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

#define N 7    // number of input neurons (decision variables)
#define M 2     // number of objectives
#define Cycles 10000     // number of simulation cycles
#define TYPE 1 // 1=max; -1=min
#define GENE 3000
#define POP_SIZE 60
#define P_MUTATION 0.1
#define P_CROSSOVER 0.2
#define rate 0.006   // interest rate per month


static int CHROMOSOME[POP_SIZE+1][N+1];
static double OBJECTIVE[POP_SIZE+1][M+1], q[POP_SIZE+1];
static void   Initialization(void);
static void   evaluation(int gen);
static void   selection(void);
static void   crossover(void);
static void   mutation(void);
static void   Objectives(void);
static void   rearrange(int x[]);
static void   f(int x[N+1], double y[M+1]);
static int    constraint_check(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(int x[N+1], double y[M+1])
{
   int i;
   double mu[Cycles+1],sum[Cycles+1],P[6],xi[12],t,temp1,temp2,temp3,temp4;
  
 static double c[12]={0,1500,1800,430,1600,800,500,2000,2100,550,530,630};
     
     temp1=temp2=temp3=temp4=0;
   for(i=1; i<=Cycles; i++) {
                sum[i]=0;        
                P[1]=x[1];
		P[2]=x[1];
		P[3]=x[1];
		P[4]=x[1];
		P[5]=x[1];

		xi[1]=myu(9-0.75,9+1);
			mu[i]=triangle(xi[1],7,9,12);
            
			xi[2]=myu(5-0.75,5+0.75);
			t=triangle(xi[2],3,5,7);
			if(mu[i]>t) mu[i]=t;
						
			xi[3]=myu(10-1,10+0.75);
            t=triangle(xi[3],7,10,12);
			if(mu[i]>t) mu[i]=t;
            
			xi[4]=myu(6-0.75,6+1);
			t=triangle(xi[4],4,6,9);
			if(mu[i]>t) mu[i]=t;

			xi[5]=myu(8-0.75,8+1);
			t=triangle(xi[5],6,8,11);
			if(mu[i]>t) mu[i]=t;
			
			
			xi[6]=myu(8-0.75,8+0.75);
			t=triangle(xi[6],6,8,10);
			if(mu[i]>t) mu[i]=t;
			

			xi[7]=myu(9-0.75,9+1);
			t=triangle(xi[7],6,9,12);
			if(mu[i]>t) mu[i]=t;
			
			xi[8]=myu(6-0.75,6+0.75);
			t=triangle(xi[8],5,6,8);
			if(mu[i]>t) mu[i]=t;
            
			xi[9]=myu(10-0.75,10+0.75);
			t=triangle(xi[9],8,10,12);
			if(mu[i]>t) mu[i]=t;
            
			xi[10]=myu(16-1,16+0.75);
			t=triangle(xi[10],13,16,18);
			if(mu[i]>t) mu[i]=t;
            
			xi[11]=myu(11-1,11+1);
			t=triangle(xi[11],9,11,13);
			if(mu[i]>t) mu[i]=t;

               
          
            
			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[i]+=c[1]*pow(1+rate,ceil(P[0]-x[1]));
			sum[i]+=c[2]*pow(1+rate,ceil(P[0]-x[1]));
			sum[i]+=c[3]*pow(1+rate,ceil(P[0]-x[1]));
			sum[i]+=c[4]*pow(1+rate,ceil(P[0]-x[2]));
			sum[i]+=c[5]*pow(1+rate,ceil(P[0]-x[3]));
			sum[i]+=c[6]*pow(1+rate,ceil(P[0]-x[3]));
			sum[i]+=c[7]*pow(1+rate,ceil(P[0]-x[3]));
			sum[i]+=c[8]*pow(1+rate,ceil(P[0]-x[4]));
			sum[i]+=c[9]*pow(1+rate,ceil(P[0]-x[5]));
			sum[i]+=c[10]*pow(1+rate,ceil(P[0]-x[6]));
			sum[i]+=c[11]*pow(1+rate,ceil(P[0]-x[7]));
                 if(P[0]-x[1] > 32) {
                if(mu[i]>temp2)
                 temp2=mu[i];
			}
            else {
	           if(mu[i]>temp1)
                 temp1=mu[i];
			} 
                 if(sum[i] > 14370 )  {
                if(mu[i]>temp4)
                 temp4=mu[i];
			}
            else {
	           if(mu[i]>temp3)
                 temp3=mu[i];
			}       
   }
          //printf("%f,%f,%f,%f\n",temp1,temp2,temp3,temp4);
        y[2] = (temp1+1-temp2)/2.;
		y[1] = (temp3+1-temp4)/2.;
      // printf("%f\n", y[2]);
}


static void Objectives(void)
{
	double  y[M+1];
	int i,j,x[N+1];

	for(i=1; i<=POP_SIZE; i++) {
		for(j=1;j<=N;j++) x[j] = CHROMOSOME[i][j];
		f(x,y);
		//for(j=1;j<=M;j++) printf("%5.8f,    ",y[j]);
        //printf("\n");
        for(j=1;j<=M;j++) OBJECTIVE[i][j] = y[j];
		}
	for(i=1; i<=POP_SIZE; i++)
	  OBJECTIVE[i][0]= OBJECTIVE[i][1];
	                    
}


static int constraint_check(int x[])
{
	double y[M+1]; 
	for(int i=1;i<=N;i++)
	{
	if( x[i]<0) return 0;
	}    
     f(x,y);
	if(y[2]<0.9) 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,12));
	  x[3]=int(myu(0,12));
	  x[2]=int(myu(0,12));	  
	  x[1]=int(myu(0,10));
	  
          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\n",i);
  }
}



int main( void )
{
   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<=GENE; 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("\nGeneration 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(")\n f=%8.15f\n", OBJECTIVE[0][1]);
	  else {
	      printf(")\nf=(");
	      for(j=1; j<=M; j++) {
		     if(j<M) printf("%8.15f,", OBJECTIVE[0][j]);
		     else printf("%8.15f", OBJECTIVE[0][j]);
		  }
          printf(")  Aggregating Value=%8.15f\n",OBJECTIVE[0][0]);
	  }
   }
   printf("\nThe final result has been written in RESULT.dat.\n\n");
   fclose(fp);
   return 1;
}








static void evaluation(int gen)
{
  double a;
  int   i, j, k, label, b;

  Objectives();

  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++) {
			 b = CHROMOSOME[i][j];
			 CHROMOSOME[i][j] = CHROMOSOME[label][j];
			 CHROMOSOME[label][j] = b;
		 }
	  }
  }
}

static void selection(void)
{
  double rr;
  int   i, j, k, temp[POP_SIZE+1][N+1];

  for(i=1; i<=POP_SIZE; i++) {
	  rr = myu(0, q[POP_SIZE]);
	  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(void)
{
  int   i, j, jj, k, x[N+1], y[N+1];
  double rr;

  for(i = 1; i <= POP_SIZE/2; i++) {
	 if(myu(0,1)>P_CROSSOVER) continue;
	 j = (int)myu(1,POP_SIZE);
	 jj = (int)myu(1,POP_SIZE);
	 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); // rearrange the schedule x
	 if(constraint_check(x) == 1)
		 for(k=1; k<=N; k++) CHROMOSOME[j][k] = x[k];

        rearrange(y); // rearrange the schedule y
        if(constraint_check(y) == 1)
		 for(k=1; k<=N; k++) CHROMOSOME[jj][k] = y[k];
  }
}



static void mutation(void)
{
  int i, j, y[N+1];
  double  infty, direction[N+1];
  double precision = 0.0001;
  double INFTY = 20;

  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++) y[j] = int(CHROMOSOME[i][j]+infty*direction[j]);

                  rearrange(y); // rearrange the schedule y
		  if(constraint_check(y) == 1) {
			 for(j=1; j<=N; j++) CHROMOSOME[i][j] = y[j];
			 break;
		  }
		  infty = myu(0,infty);
	  }
  }
}