/************************************************************************
 *
 * Author:	 Luis Gerardo de la Fraga	 fraga@cs.cinvestav.mx
 *
 * Computer Science Deparment, Cinvestav.
 *
 * Copyright (c) 2007, Cinvestav.
 *
 * Permission is granted to copy and distribute this file, for noncommercial
 * use, provided (a) this copyright notice is preserved, and (b) no attempt
 * is made to restrict redistribution of this file, 

 *  All comments concerning this program may be sent to the
 *  e-mail address 'fraga@cs.cinvestav.mx'
 *
 **************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

/**
#define DEBUG 0
   Activating this DEBUG flag produce:
   The best guy is printed every 20 generations
   And this flag can be activated at compiled time as:

   gcc -O2 -DDEBUG -o edellipse edellipse.c -lm
**/
#define NPARAM  5

#define A  0
#define B  1
#define XC 2
#define YC 3
#define ALPHA 4
#define FIT   5

typedef struct param {
	double p[6];
} PARAM;

typedef struct point {
	double x, y;
} POINT;

typedef struct data {
	int n;
	POINT *pts;
	POINT limits[NPARAM]; /** For a, b, xc, yc, and alpha **/
	int ngen, popsize;
	PARAM *guys;
	double Pc, Pf; 
	double s;  /** stop condition **/
} DATA;

float ran0( int * idum, int seed )
{
	#define SIZE 100
	static float maxran, v[SIZE], y;
	static int iff=0;
	int j;

	if ( *idum <0 || iff == 0) {
		srandom( seed );
		iff=1;
		*idum=1;
		maxran = RAND_MAX + 1.0;
		for (j=0;j<100;j++) v[j]=random();
		y=random();
	}
	j= (int)( y*SIZE/maxran );
	y=v[j];
	v[j]=random();
	return (y/maxran);
}

double find_point_in_elipse( double a, double b, double x0, double y0 )
{
	double l0, l1, v, v1, v2;
	double f, df, s, c;
	double p;
	int i;

	/** Correction added Jan 12, 2010 **/
	l0 = sqrt( x0*x0 + y0*y0 );
	if ( l0 < 1e-6 ) return M_PI_2;

	v = a*a - b*b;

	if( fabs(y0) < b && fabs(x0) <= v/a ) {
		if ( fabs(y0) < 1e-6 ) {
			l0 = M_PI_2;
		}
		else {
			p = fabs(x0)/(1 + fabs(y0)/l0);
			if( x0 < 0 ) { p = -p; }
			l0 = atan2( a*y0, b*(x0-p) );
		}
	}
	else {
		l0 = atan2( a*y0, b*x0 );
	}

	v1 = y0*b;
	v2 = x0*a;
	for ( i=0; i<=3; i++ ) {
		s = sin(l0);	
		c = cos(l0);	
		f = v1*c - v2*s + v*c*s;
		df = -v1*s - v2*c + v*(c*c-s*s);

		if ( fabs(df) < 1e-8 ) {
			fprintf( stderr,  "# NO $x0 $y0\n" );
            return ( l0 );
        }

		l1 = l0 - f/df;
		if ( fabs(l1-l0)< 0.00001 ){
			return( l1 );
		}

		l0 = l1;
	}
	return l0;
}


void evaluateGuy( PARAM *ptrg, POINT *ptrp, int n ) 
{
	int i;
	double c, s, d, m;
	POINT ptran, prota, dif;

#ifdef CONSTR
	ptrg->p[B] = 20.0/ptrg->p[A];
#endif

	if ( ptrg->p[A] < ptrg->p[B] ) {
		s = ptrg->p[A];
		ptrg->p[A] = ptrg->p[B];
		ptrg->p[B] = s;
	}	

	c = cos( ptrg->p[ALPHA] );
	s = sin( ptrg->p[ALPHA] );

	d = 0;
	for ( i=0; i<n; i++ ) {
		ptran.x = ptrp[i].x - ptrg->p[XC];  
		ptran.y = ptrp[i].y - ptrg->p[YC];  

		prota.x =  ptran.x*c + ptran.y*s;
		prota.y = -ptran.x*s + ptran.y*c;

		m = find_point_in_elipse( ptrg->p[A], ptrg->p[B], prota.x, prota.y ); 

		dif.x = ptrg->p[A]*cos(m) - prota.x;
		dif.y = ptrg->p[B]*sin(m) - prota.y;

		/** d += fabs(dif.x) + fabs(dif.y); **/

		d += sqrt( dif.x*dif.x + dif.y*dif.y );

		/** d += dif.x*dif.x + dif.y*dif.y; **/
	}
	/** fprintf( stderr, "d = %lf\n", d ); **/
	

	ptrg->p[FIT] = d;
}

/** The main algorithm **/	
void eDifferential( DATA *pd )
{
	int i, i0, i1, i2;
	int n, seed=6;
	int j, k, jrand;
	double v, m, m2;
	PARAM newguy;

	/** Evaluate to all guys **/
	for(i=0; i<pd->popsize; i++)
		evaluateGuy( &pd->guys[i], pd->pts, pd->n ); 

	n = pd->popsize;
	for(i=0; i<pd->ngen; i++){
		for(k=0; k<pd->popsize; k++){
			/** We select three different parents **/
			i0 = (int)(ran0( &seed, k )*n);
			do { 
				i1 = (int)(ran0( &seed, k )*n);
			} while( i1 == i0 );
			do { 
				i2 = (int)(ran0( &seed, k )*n);
			} while( i2==i1 || i2 == i0 );

			/** We apply the crossoverBin **/
#ifndef CONSTR
			jrand = (int)(5.0*ran0( &seed, k ));
#else
			jrand = (int)(3.0*ran0( &seed, k ));
			if ( jrand >= B ) jrand++;
#endif

#ifndef CONSTR
			for ( j=0; j<NPARAM; j++ ) {
#else
			for ( j=0; j<=3; j++ ) {
				/** An example of how to use constraints b=a/20, alpha=pi/2 **/
				if ( j == B ) j++;
#endif
				if ( (ran0( &seed, k ) < pd->Pc) || j==jrand ) {
					newguy.p[j] = pd->guys[i2].p[j] +
						pd->Pf * (pd->guys[i0].p[j] - pd->guys[i1].p[j]); 

				/** We avoid this value be outside its bounds **/
				if ( newguy.p[j] < pd->limits[j].x ||
					 newguy.p[j] > pd->limits[j].y ) {
					newguy.p[j] = (double)(((double)(
					   (long)(fabs(newguy.p[j]*1e8))%
							 (long)((pd->limits[j].y-pd->limits[j].x)*1e8))/1e8)
							+ pd->limits[j].x);
					}
				}
				else{
					newguy.p[j] = pd->guys[k].p[j];
				}
			}
			/** Evaluation of the new guy 
			fprintf( stderr, "%lf %lf %lf %lf %lf\n",
				 pd->guys1[k].p[A], pd->guys1[k].p[B], pd->guys1[k].p[XC],
				 pd->guys1[k].p[YC], pd->guys1[k].p[ALPHA] );
			**/
#ifdef CONSTR
			pd->guys1[k].p[ALPHA] = pd->guys[k].p[ALPHA];
#endif
			evaluateGuy( &newguy, pd->pts, pd->n ); 

			/** Copy the new guy if it is better **/
			if( newguy.p[FIT] < pd->guys[k].p[FIT] ){
				for ( j=0; j<=NPARAM; j++ ) 
					pd->guys[k].p[j] = newguy.p[j];
			}
		} /** End of one generation **/
		/** We find the best and the worst guys **/
		i0 = i1 = 0;
		for ( j=1; j<n; j++ ) {
			if ( pd->guys[j].p[FIT] < pd->guys[i0].p[FIT] )
				i0 = j;
		}
		for ( j=1; j<n; j++ ) {
			if ( pd->guys[j].p[FIT] > pd->guys[i1].p[FIT] )
				i1 = j;
		}
		if ( pd->guys[i1].p[FIT] - pd->guys[i0].p[FIT] < pd->s ) 
			break;
	}
	printf( "%d %lf %lf %lf %lf %lf %lf\n", i,
		pd->guys[i0].p[A], pd->guys[i0].p[B], 
		pd->guys[i0].p[XC], pd->guys[i0].p[YC], 
		pd->guys[i0].p[ALPHA], pd->guys[i0].p[FIT]  ); 
}



int initialization( DATA *pd )
{
	pd->pts = (POINT *)malloc( pd->n * sizeof(POINT) );
	if ( pd->pts == NULL ) return 1;

	pd->guys = (PARAM *)malloc( pd->popsize * sizeof(PARAM) );
	if ( pd->guys == NULL ){
		free( pd->pts );
		return 1;
	}
	
	return 0;
}

int read_input_data( char *namefile, DATA *pd )
{
	FILE *fp;
	double x, y;
	int i, n;
	char line[64];

	if ( (fp=fopen( namefile, "r" ))==NULL ){
		fprintf( stderr, "Error: can't open %s file\n", namefile );
		return 1;
	}	
	i=0;
	while ( fgets( line, 64, fp ) != NULL ) {
		n = sscanf( line, "%lf %lf", &x, &y );

		if ( n == 2 ){
			pd->pts[i].x = x;
			pd->pts[i].y = y;
			i++;
		} 

	}
	fclose( fp );
	/**
	for (i=0; i<pd->n; i++ ) 
		printf ("%d %lf %lf\n", i, pd->pts[i].x, pd->pts[i].y );
	pd->guys[0].p[A] = 6.51872; 
	pd->guys[0].p[B] = 3.03189;
	pd->guys[0].p[XC] = 2.69962;
	pd->guys[0].p[YC] = 3.81596;
	pd->guys[0].p[ALPHA] = 0.35962;

	evaluateGuy( &pd->guys[0], pd->pts, pd->n ); 
	printf( "Error = %lf\n", pd->guys[0].p[FIT] );
**/
	 
	if ( i != pd->n ) {
		fprintf( stderr, "Error: I read %d points\n", i );
		return 1;
	}
	return 0;
}

void free_memory ( DATA *pd )
{
	free( pd->guys );
	free( pd->pts );
}


double halton( int i, int d )
{
	/**
	 hasta 5 dimensiones (en d)
	 para más dimensiones no es buena esta aproximación,
	 hay que usar la version scrambled
	**/
	static int prime[5] = { 2, 3, 5, 7, 11 };
	int x, p1, p2;
	double c;

	p1 = prime[d];
	p2 = p1;
	c = 0.0;
	do {
		x = i % p1;
		c += (double)x/p2;
		i /= p1;
		p2 *= p1;
	} while( i > 0 );

	return( c );
}

void dhalton( int i, double * r, double * u )
{
	double v;
	int d;

	for(  d=0; d<5; d++ ) {
		v = halton( i, d );
		v += r[d];
		if ( v >= 1.0 ) v -= 1.0;
		u[d] = v;
	}
}


int main( int argc, char * argv[] ) 
{
	DATA d;
	int i, j; 
	double range[5], u[5], r[5];
	int initran=-3;   /** A negative value restart the random number generator **/
	int seed;
	char file[64];

	if( argc != 19 ) { 
		fprintf( stderr, "Args: a_min a_max b_min b_max xc_min xc_max\n\tyc_min yc_max alpha_min alpha_max n data_file ngen popsize Pc Pf s seed\n" );
		return 1;
	}

	d.limits[0].x = strtod( argv[1], (char **)NULL);
	d.limits[0].y = strtod( argv[2], (char **)NULL);

	d.limits[1].x = strtod( argv[3], (char **)NULL);
	d.limits[1].y = strtod( argv[4], (char **)NULL);

	d.limits[2].x = strtod( argv[5], (char **)NULL);
	d.limits[2].y = strtod( argv[6], (char **)NULL);

	d.limits[3].x = strtod( argv[7], (char **)NULL);
	d.limits[3].y = strtod( argv[8], (char **)NULL);

	d.limits[4].x = strtod( argv[9], (char **)NULL);
	d.limits[4].y = strtod( argv[10], (char **)NULL);

	d.n = (int)strtol( argv[11], (char **)NULL, 10); 
	strncpy( file, argv[12], 64 );

	d.ngen = (int)strtol( argv[13], (char **)NULL, 10); 
	d.popsize = (int)strtol( argv[14], (char **)NULL, 10);
	if ( d.popsize < 3 ){
		fprintf( stderr, "Population size must be greater than 3\n" );
		return 1;
	}
	d.Pc = strtod( argv[15], (char **)NULL);
	d.Pf = strtod( argv[16], (char **)NULL);

	d.s = strtod( argv[17], (char **)NULL);
	seed = (int)strtol( argv[18], (char **)NULL, 10); 

	if( initialization( &d ) ) {
		fprintf( stderr, "Error allocating memory\n" );
		return 1;
	}
	if( read_input_data( file, &d ) ) {
		free_memory ( &d );
		fprintf( stderr, "Error reading data\n" );
		return 1;
	}
	/** fprintf( stderr, "Read %d data pairs\n", d.n ); **/

	for( j=0; j<NPARAM; j++ ) {
		range[j] = d.limits[j].y - d.limits[j].x;
		r[j] = ran0( &initran, seed );
	}
	/** printf( "I: %lf %lf %lf %lf %lf\n", r[0], r[1], r[2], r[3], r[4] );
		
	     we initialize the guys with random values 
	     Change 2/03/2011 We use now a random Halton sequence **/
	for( i=0; i<d.popsize; i++ ) {
		dhalton( i, r, u );
		/** printf( "%lf %lf %lf %lf %lf\n", u[0], u[1], u[2], u[3], u[4] ); **/
		for( j=0; j<NPARAM; j++ ) 
			d.guys[i].p[j] = d.limits[j].x + range[j]*u[j];
	}

	eDifferential( &d );

	free_memory ( &d );
	return 0;
}