/*========================================================================
  PISA  (https://sop.tik.ee.ethz.ch/pisa/)
  ========================================================================
  Computer Engineering (TIK)
  ETH Zurich
  ========================================================================
  SHV - Sampling-based HyperVolume-oriented Algorithm

  authors: Johannes Bader, johannes.bader@tik.ee.ethz.ch

  last change: 20.04.2009
  ========================================================================
 */


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <float.h>

#include "shv.h"

/* common parameters */
int alpha;  /* number of individuals in initial population */
int mu;     /* number of individuals selected as parents */
int lambda; /* number of offspring individuals */
int dim;    /* number of objectives */


/* local parameters from paramfile*/
int seed;   /* seed for random number generator */
double bound;
int nrOfJunks;
int nrOfSamplesPerJunk;

/* other variables */
char cfgfile[FILE_NAME_LENGTH];  /* 'cfg' file (common parameters) */
char inifile[FILE_NAME_LENGTH];  /* 'ini' file (initial population) */
char selfile[FILE_NAME_LENGTH];  /* 'sel' file (parents) */
char arcfile[FILE_NAME_LENGTH];  /* 'arc' file (archive) */
char varfile[FILE_NAME_LENGTH];  /* 'var' file (offspring) */

/* population containers */
pop *pp_all = NULL;
pop *pp_new = NULL;
pop *pp_sel = NULL;
/* front partitioning */
fpart front_part;

/*-----------------------| initialization |------------------------------*/
void initialize(char *paramfile, char *filenamebase)
/* Performs the necessary initialization to start in state 0. Reads all
 * parameters
 */
{
	FILE *fp;
	int result;
	char str[CFG_ENTRY_LENGTH];

	/* reading parameter file with parameters for selection */
	fp = fopen(paramfile, "r");
	assert(fp != NULL);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "seed") == 0);
	result = fscanf(fp, "%d", &seed);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "bound") == 0);
	result = fscanf(fp, "%lf", &bound);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "junks") == 0);
	result = fscanf(fp, "%d", &nrOfJunks);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "junksize") == 0);
	result = fscanf(fp, "%d", &nrOfSamplesPerJunk );

	/* fscanf() returns EOF if reading fails. */
	assert(result != EOF); /* no EOF, parameters correctly read */
	fclose(fp);

	srand(seed); /* seeding random number generator */

	sprintf(varfile, "%svar", filenamebase);
	sprintf(selfile, "%ssel", filenamebase);
	sprintf(cfgfile, "%scfg", filenamebase);
	sprintf(inifile, "%sini", filenamebase);
	sprintf(arcfile, "%sarc", filenamebase);

	/* reading cfg file with common configurations for both parts */
	fp = fopen(cfgfile, "r");
	assert(fp != NULL);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "alpha") == 0);
	fscanf(fp, "%d", &alpha);
	assert(alpha > 0);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "mu") == 0);
	fscanf(fp, "%d", &mu);
	assert(mu > 0);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "lambda") == 0);
	fscanf(fp, "%d", &lambda);
	assert(lambda > 0);

	fscanf(fp, "%s", str);
	assert(strcmp(str, "dim") == 0);
	result = fscanf(fp, "%d", &dim);
	assert(result != EOF); /* no EOF, 'dim' correctly read */
	assert(dim > 0);

	fclose(fp);

	/* create individual and archive pop */
	pp_all = create_pop(alpha + lambda, dim);
	pp_sel = create_pop(mu, dim);
	/* create temporary populations (indices only) */
	create_front_part(&front_part, alpha + lambda);
}


/*-------------------| memory allocation functions |---------------------*/

void* chk_malloc(size_t size)
/* Wrapper function for malloc(). Checks for failed allocations. */
{
	void *return_value = malloc(size);
	if(return_value == NULL)
		PISA_ERROR("Selector: Out of memory.");
	return (return_value);
}

void create_front_part(fpart*  partp, int  max_pop_size)
{
	int i;

	partp->max_fronts = (max_pop_size > 0 ? max_pop_size : 0);
	partp->fronts = 0;
	partp->size = 0;
	partp->front_array = (front*) chk_malloc(max_pop_size * sizeof(front));
	for (i = 0; i < max_pop_size; i++)
	{
		partp->front_array[i].size = 0;
		partp->front_array[i].members = chk_malloc(max_pop_size * sizeof(int));
	}
}

pop* create_pop(int maxsize, int dim)
/* Allocates memory for a population. */
{
	int i;
	pop *pp;

	assert(dim >= 0);
	assert(maxsize >= 0);

	pp = (pop*) chk_malloc(sizeof(pop));
	pp->size = 0;
	pp->maxsize = maxsize;
	pp->ind_array = (ind**) chk_malloc(maxsize * sizeof(ind*));

	for (i = 0; i < maxsize; i++)
		pp->ind_array[i] = NULL;

	return (pp);
}


ind* create_ind(int dim)
/* Allocates memory for one individual. */
{
	ind *p_ind;

	assert(dim >= 0);

	p_ind = (ind*) chk_malloc(sizeof(ind));

	p_ind->index = -1;
	p_ind->fitness = -1;
	p_ind->f = (double*) chk_malloc(dim * sizeof(double));
	return (p_ind);
}

void free_front_part(fpart* partp)
{
	int i;

	if (partp->front_array != NULL)
	{
		for (i = 0; i < partp->max_fronts; i++)
			if (partp->front_array[i].members != NULL)
				chk_free(partp->front_array[i].members);
		chk_free(partp->front_array);
		partp->front_array = NULL;
	}
	partp->fronts = 0;
	partp->max_fronts = 0;
	partp->size = 0;
}

void free_memory()
/* Frees all memory. */
{
	free_pop(pp_sel);
	complete_free_pop(pp_all);
	free_pop(pp_new);
	pp_sel = NULL;
	pp_all = NULL;
	pp_new = NULL;
	free_front_part(&front_part);
}

void free_pop(pop *pp)
/* Frees memory for given population. */
{
	if (pp != NULL)
	{
		free(pp->ind_array);
		free(pp);
	}
}

void complete_free_pop(pop *pp)
/* Frees memory for given population and for all individuals in the
   population. */
{
	int i = 0;
	if (pp != NULL)
	{
		if(pp->ind_array != NULL)
		{
			for (i = 0; i < pp->size; i++)
			{
				if (pp->ind_array[i] != NULL)
				{
					free_ind(pp->ind_array[i]);
					pp->ind_array[i] = NULL;
				}
			}

			free(pp->ind_array);
		}

		free(pp);
	}
}

void free_ind(ind *p_ind)
/* Frees memory for given individual. */
{
	assert(p_ind != NULL);

	free(p_ind->f);
	free(p_ind);
}

/*-----------------------| selection functions|--------------------------*/

void selection()
/**
 *  Do environmental and mating selection
 */
{
	mergeOffspring();
	generateFrontPartition( &front_part );
	environmentalSelection();
	cleanUpArchive(&(front_part));
	matingSelection();

	return;
}

void cleanUpArchive(fpart* partp)
/** Removes all individuals from pp_all that are not referenced in partp
 *
 * @param[in] partp partition of the population
 */
{
	pop* new_pop;
	int i, j;

	new_pop = create_pop(alpha + lambda, dim);
	new_pop->size = 0;
	for (i = partp->fronts - 1; i >= 0; i--)
	{
		for (j = partp->front_array[i].size - 1; j >= 0; j--)
		{
			int index = partp->front_array[i].members[j];
			new_pop->ind_array[new_pop->size] = pp_all->ind_array[index];
			pp_all->ind_array[index] = NULL;
			partp->front_array[i].members[j] = new_pop->size;
			new_pop->size++;
		}
	}
	complete_free_pop(pp_all);
	pp_all = new_pop;
}

void mergeOffspring()
/**
 * Merge the offspring individuals (pp_new) with the archive population pp_new
 *
 * @pre pp_all->maxsize >= pp_all->size + pp_new->size
 * @remark the population pp_new is added to pp_all, the population pp_new
 * 	is emptied.
 */
{
	int i;

	assert(pp_all->size + pp_new->size <= pp_all->maxsize);

	for (i = 0; i < pp_new->size; i++)
	{
		pp_all->ind_array[pp_all->size + i] = pp_new->ind_array[i];
	}

	pp_all->size += pp_new->size;

	free_pop(pp_new);
	pp_new = NULL;
}

int dominates(int a, int b)
/* Determines if one individual dominates another.
   Minimizing fitness values. */
{
	int i;
	int a_is_worse = 0;
	int equal = 1;

	for (i = 0; i < dim && !a_is_worse; i++)
	{
		a_is_worse = pp_all->ind_array[a]->f[i] > pp_all->ind_array[b]->f[i];
		equal = (pp_all->ind_array[a]->f[i] == pp_all->ind_array[b]->f[i]) && equal;
	}

	return (!equal && !a_is_worse);
}

void generateFrontPartition( fpart* front_part )
/**
 * Partition the population pp_all into fronts
 *
 * @param[out] front_part partition of the front
 * @pre the population pp_all needs to be valid
 */
{
	int i,j;
	int actFront = 0;
	int notDominated;
	int checked[ pp_all->size ];
	int actchecked[ pp_all->size ];
	int added;

	cleanPart( front_part );

	for( i = 0; i < pp_all->size; i++ )
		checked[i] = 0;

	added = 0;
	actFront = 0;
	while( added < pp_all->size )
	{
		for( i = 0; i < pp_all->size; i++ )
		{
			if( checked[i] == 1 )
				continue;
			actchecked[i] = 0;
			notDominated = 1;
			j = 0;
			while( notDominated == 1 && j < pp_all->size )
			{
				if( i != j && checked[j] == 0 && dominates(j, i) )
					notDominated = 0;
				j++;
			}
			if( notDominated )
			{
				actchecked[i] = 1;
				insertInPart(front_part, actFront, i );
				added++;
			}
		}
		for( j = 0; j < pp_all->size; j++ )
			if( actchecked[j] == 1 )
				checked[j] = 1;
		actFront++;
	}
}

void cleanPart(fpart* partp)
{
	int i;

	for (i = 0; i < partp->max_fronts; i++)
		partp->front_array[i].size = 0;
	partp->size = 0;
	partp->fronts = 0;
}

void insertInPart( fpart* partp, int nr, int index )
{
	if( partp->fronts < nr + 1 )
		partp->fronts = nr + 1;
	partp->size++;
	partp->front_array[ nr ].members[ partp->front_array[ nr ].size ] = index;
	partp->front_array[ nr ].size++;
}

void removeIndividual( int sel, fpart* partp, front *fp )
{
	assert( sel >= 0 );
	assert( sel < fp->size );
	fp->size--;
	fp->members[sel] = fp->members[fp->size];
	partp->size--;
}

void  swap(double  *front, int  i, int  j)
{
	int  k;
	double  temp;

	for (k = 0; k < dim; k++) {
		temp = front[i * dim + k];
		front[i * dim + k] = front[j * dim + k];
		front[j * dim + k] = temp;
	}
}

int  dominatesPoints(double  *point1, double  *point2, int  no_objectives)
/* returns true if 'point1' dominates 'points2' with respect to the
       to the first 'no_objectives' objectives
 */
{
	int  i;
	int  better_in_any_objective, worse_in_any_objective;

	better_in_any_objective = 0;
	worse_in_any_objective = 0;
	for (i = 0; i < no_objectives && !worse_in_any_objective; i++)
		if (point1[i] > point2[i])
			better_in_any_objective = 1;
		else if (point1[i] < point2[i])
			worse_in_any_objective = 1;
	return (!worse_in_any_objective && better_in_any_objective);
}

int  filter_nondominated_set(double  *front, int  no_points,
		int  no_objectives)
/* all nondominated points regarding the first 'no_objectives' dimensions
       are collected; the points 0..no_points-1 in 'front' are
       considered; the points in 'front' are resorted, such that points
       [0..n-1] represent the nondominated points; n is returned
 */
{
	int  i, j;
	int  n;

	n = no_points;
	i = 0;
	while (i < n) {
		j = i + 1;
		while (j < n) {
			if (dominatesPoints(&(front[i * dim]), &(front[j * dim]),
					no_objectives)) {
				/* remove point 'j' */
				n--;
				swap(front, j, n);
			}
			else if (dominatesPoints(&(front[j * dim]), &(front[i * dim]),
					no_objectives)) {
				/* remove point 'i'; ensure that the point copied to index 'i'
		   is considered in the next outer loop (thus, decrement i) */
				n--;
				swap(front, i, n);
				i--;
				break;
			}
			else
				j++;
		}
		i++;
	}
	return n;
}

double  surface_unchanged_to(double  *front, int  no_points, int  objective)
/* calculate next value regarding dimension 'objective'; consider
	points 0..no_points-1 in 'front'
 */
{
	int     i;
	double  min, value;

	min = front[objective];
	for (i = 1; i < no_points; i++) {
		value = front[i * dim + objective];
		if (value < min)  min = value;
	}

	return min;
}

int  reduce_nondominated_set(double  *front, int  no_points, int  objective,
		double  threshold)
/* remove all points which have a value <= 'threshold' regarding the
       dimension 'objective'; the points [0..no_points-1] in 'front' are
       considered; 'front' is resorted, such that points [0..n-1] represent
       the remaining points; 'n' is returned
 */
{
	int  n;
	int  i;

	n = no_points;
	for (i = 0; i < n; i++)
		if (front[i * dim + objective] <= threshold) {
			n--;
			swap(front, i, n);
		}

	return n;
}

double  calc_hypervolume(double  *front, int  no_points, int  no_objectives)
{
	int     n;
	double  volume, distance;

	volume = 0;
	distance = 0;
	n = no_points;
	while (n > 0) {
		int     no_nondominated_points;
		double  temp_vol, temp_dist;

		no_nondominated_points = filter_nondominated_set(front, n,
				no_objectives - 1);
		temp_vol = 0;
		if (no_objectives < 3) {
			temp_vol = front[0];
		}
		else
			temp_vol = calc_hypervolume(front, no_nondominated_points,
					no_objectives - 1);
		temp_dist = surface_unchanged_to(front, n, no_objectives - 1);
		volume += temp_vol * (temp_dist - distance);
		distance = temp_dist;
		n = reduce_nondominated_set(front, n, no_objectives - 1, distance);
	}

	return volume;
}

double unaryHypervolumeIndicator(int *A, int sizea, double bound  )
{
	double pointArray[ sizea*dim ];
	int i, k;
	for (i = 0; i < sizea; i++)
		for (k = 0; k < dim; k++)
			pointArray[i * dim + k] = bound - pp_all->ind_array[ A[i] ]->f[k];

	return calc_hypervolume(pointArray, sizea, dim );
}

void assignFitness(front* fp )
/* fitness is to be maximized */
{
	int i;
	front* fpTemp1;
	fpart dum1_part;  /* temporary storage */
	create_front_part(&dum1_part, alpha + lambda);

	double indicatorValue;

	fpTemp1 = &(dum1_part.front_array[0]);
	for (i = 0; i < fp->size; i++)
		fpTemp1->members[i] = fp->members[i];
	dum1_part.size = fp->size;
	dum1_part.fronts = 1;
	dum1_part.front_array[0].size = fp->size;


	fpTemp1->size--;
	for (i = fpTemp1->size; i >= 0; i--)
	{
		int temp = fpTemp1->members[i];
		fpTemp1->members[i] = fpTemp1->members[fpTemp1->size];
		indicatorValue = unaryHypervolumeIndicator( fpTemp1->members, fpTemp1->size,
				bound );
		pp_all->ind_array[temp]->fitness = -indicatorValue;
		fpTemp1->members[fpTemp1->size] = temp;
	}


	free_front_part(&dum1_part);
}

void removeTheWorstIndividual(front* fp, fpart* partp )
{
	int j;
	int sel = 0;
	double min = -1;

	for (j = fp->size - 1; j >= 0; j--)
	{
		if ( j == fp->size - 1 || min >= pp_all->ind_array[fp->members[j]]->fitness)
		{
			min = pp_all->ind_array[fp->members[j]]->fitness;
			sel = j;
		}
	}
	removeIndividual( sel, partp, fp );
}
void getObjectiveVector( int* A, int sizea, int i, double* objVect  )
{
	assert( i >= 0);
	assert( i < sizea );
	int k;
	for( k = 0; k < dim; k++ ) {
		objVect[ k ] = pp_all->ind_array[ A[i] ]->f[k];
	}
}

void getObjectiveArray( int* A, int sizea, double* pointArray )
/* Returns an array of all objectives referenced by A */
{
	int i,k;
	for( i = 0; i < sizea; i++ )
	{
		for( k = 0; k < dim; k++ )
		{
			pointArray[ i*dim + k ] = pp_all->ind_array[ A[i] ]->f[k];
		}
	}
}

int getNextSamplingIndex( int* h, int* N, double* V, int sizea, int* save,
		int decreaseDependent )
/* Get the next sampling index. All save points ignored, all points with 0
 * sample sizes prioritized */
{
	double sel = -1;
	double maxDecrease = -1;

	int i;
	double pt, actDecrease;
	for( i = 0; i < sizea; i++ )
	{
		if( save[i] > 0 )
			continue;
		if( N[i] == 0 )
		{
			sel = i;
			break;
		}
		if( decreaseDependent == 0 )
			actDecrease = rand();
		else {
			pt = ( h[i] + 1 ) / ( (double)N[i] + 2.0 );
			actDecrease = V[i]*V[i]*( pt*(1-pt) / ( (double)N[i] + 3.0 ) -
					pt*(1-pt) / ( (double)N[i] + 4.0 ) );
		}
		if( sel == -1 || actDecrease > maxDecrease )
		{
			sel = i;
			maxDecrease = actDecrease;
		}
	}
	return sel;
}

void quick_sort_indices(double values[], int indices[], int left, int right )
{
	double pivot;
	int pivot_ind;
	int l_hold, r_hold;

	l_hold = left;
	r_hold = right;
	pivot = values[left];
	pivot_ind = indices[left];
	while (left < right)
	{
		while ((values[right] >= pivot) && (left < right))
			right--;
		if (left != right)
		{
			values[left] = values[right];
			indices[left] = indices[right];
			left++;
		}
		while ((values[left] <= pivot) && (left < right))
			left++;
		if (left != right)
		{
			values[right] = values[left];
			indices[right] = indices[left];
			right--;
		}
	}
	values[left] = pivot;
	indices[left] = pivot_ind;
	pivot = left;
	left = l_hold;
	right = r_hold;
	if (left < pivot)
		quick_sort_indices(values, indices, left, pivot-1 );
	if (right > pivot)
		quick_sort_indices(values, indices, pivot+1, right );
}

void calcSamplingHypercubes( int *A, int sizea, double* SHC_lower,	double* SHC_upper )
/* Calculate the Sampling hypercubes for all individuals in A
 * Assumptions: Minimizing */
{
	getObjectiveArray(A, sizea, SHC_lower);

	int i,k,j;
	int Ind[ sizea ];
	double objVect[ sizea ];

	/* All dimensions */
	for( k = 0; k < dim; k++ )
	{
		for( i = 0; i < sizea; i++ )
			Ind[i] = i;
		getObjectiveVector( A, sizea, k, objVect  );
		quick_sort_indices( objVect, Ind, 0, sizea-1 );
		/* All individuals */
		for( j = 0; j < sizea; j++ )
		{
			int found = 0;
			int l = j;
			/* If there are multiple individuals with the same kth coordinate,
			 * check all */
			while( l > 0 && objVect[l-1] == objVect[l] )
				l--;

			do {
				if( l == j )
					continue;
				int o = 0;
				found = 1;
				while( found == 1 && o < dim)
				{
					if( o != k )
						found = (pp_all->ind_array[A[Ind[j]]]->f[o] -
								pp_all->ind_array[A[Ind[l]]]->f[o] ) >= 0;
					o++;
				}
			} while( found == 0 && ++l < sizea );
			assert( (Ind[j]*dim + k) >= 0 );
			assert( (Ind[j]*dim + k) < dim*sizea );
			if( found )
				SHC_upper[ Ind[j]*dim + k ] = pp_all->ind_array[ A[ Ind[l] ] ]->f[k];
			else
				SHC_upper[ Ind[j]*dim + k ] = bound;

			/* Debug Checks */
			if( SHC_upper[ Ind[j]*dim + k ] < SHC_lower[ Ind[j]*dim + k ] )
			{
				fprintf(stderr,"Error: Upper bound (%lf) smaller than lower "
						"bound (%lf)\n", SHC_upper[ Ind[j]*dim + k ],
						SHC_lower[ Ind[j]*dim + k ] ); fflush(stderr);
			}
		}
	}
}

void resetSamplingCounts( int* h, int* N, int size )
{
	while( --size >= 0 )
	{
		h[ size ] = 0;
		N[ size ] = 0;
	}
}


void calcSamplingHypercubesVolume( double* SHCVol, double* SHC_lower, double* SHC_upper,
		int size )
{
	int i,k;

	for( i = 0; i < size; i++ )
	{
		SHCVol[i] = 1.0;
		for( k = 0; k < dim; k++ )
			SHCVol[i] *= SHC_upper[ i*dim + k ] - SHC_lower[ i*dim + k ];
	}
}

int weaklyDominatesMin( double *point1, double *point2, int no_objectives )
{
	int better;
	better = 1;
	int i = 0;

	while( i < no_objectives && better )
	{
		better = point1[i] <= point2[i];
		i++;
	}
	return better;
}

void monteCarloSampling( int* h, int* N, int* A, int sizea, double* SHC_lower,
		double* SHC_upper, int nrOfJunks, int nrOfSamplesPerJunk, double* SHCVol )
{
	double sample[ dim ];
	int save[ sizea ];
	int i,c,s,k,j;
	for( i = 0; i < sizea; i++ )
		save[i] = ( SHCVol[i] == 0);

	for( c = 0; c < nrOfJunks; c++ )
	{
		i = getNextSamplingIndex( h, N, SHCVol, sizea, save, 1 );

		assert( i >= 0 );
		assert( i < sizea );

		/* In this case, sampling makes no sense */
		if( SHCVol[i] == 0 )
			continue;

		N[i] += nrOfSamplesPerJunk;

		for( s = 0; s < nrOfSamplesPerJunk; s++ )
		{
			for( k = 0; k < dim; k++ )
				sample[ k ] = drand( SHC_lower[ i*dim + k ],
						SHC_upper[ i*dim + k ] );

			int single = 1;
			j = 0;
			while( single == 1 && j < sizea )
			{
				if( j != i )
					single = !weaklyDominatesMin( SHC_lower + (j*dim), sample, dim);
				j++;
			}
			if( single )
				h[i]++;
		}
	}
}

void monteCarloSampleNewAfterRemoval( fpart* partp, front* fp, int nrOfJunks,
		int nrOfSamplesPerJunk )
{
	int totalNrOfSamples = 0;
	double* SHCVol = NULL;
	double* SHC_upper = NULL;
	double* SHC_lower = NULL;
	int h[ fp->size];
	int N[ fp->size];
	double minVol = DBL_MAX;
	double checkVol;
	int sel = -1;
	int i;

	SHC_lower = (double*)chk_malloc( sizeof(double)*(fp->size)*dim );
	SHC_upper  = (double*)chk_malloc( sizeof(double)*(fp->size)*dim );
	SHCVol = (double*)chk_malloc( sizeof(double)*fp->size );

	while( partp->size > alpha )
	{

		calcSamplingHypercubes( fp->members, fp->size, SHC_lower, SHC_upper);
		calcSamplingHypercubesVolume( SHCVol, SHC_lower, SHC_upper, fp->size );
		resetSamplingCounts( h, N, fp->size);
		monteCarloSampling( h, N, fp->members, fp->size, SHC_lower, SHC_upper,
				nrOfJunks, nrOfSamplesPerJunk, SHCVol );
		totalNrOfSamples += nrOfJunks*nrOfSamplesPerJunk;

		sel = -1;
		for( i = 0; i < fp->size; i++ )
		{
			if( N[i] > 0 )
				checkVol = (h[i]/(double)N[i] )*SHCVol[i];
			else
				checkVol = 0;
			if( sel == -1  || checkVol < minVol ) {
				minVol = checkVol;
				sel = i;
			}
		}
		assert( sel >= 0 );
		assert( sel < fp->size );
		removeIndividual( sel, partp, fp );
	}

	chk_free(SHCVol); SHCVol = NULL;
	chk_free(SHC_lower); SHC_lower = NULL;
	chk_free(SHC_upper); SHC_upper = NULL;
}

void environmentalSelection()
/**
 * Selects alpha individuals of the population pp_all
 *
 * @pre pp_all and front_part set
 * @post front_part.size == alpha
 * @remark operates only on the front partition front_part. No individuals
 * 		are actually removed.
 */
{
	int i;
	/** Start with front wise reduction */
	for( i = front_part.fronts - 1; i >= 0; i-- )
	{
		if( front_part.size - front_part.front_array[i].size >= alpha )
		{
			front_part.size -= front_part.front_array[i].size;
			front_part.front_array[i].size = 0;
			front_part.fronts--;
		}
		else
			break;
	}
	/** Then remove from worst front */
	front* fp = &( front_part.front_array[ front_part.fronts-1 ] );
	while( front_part.size > alpha )
	{
		monteCarloSampleNewAfterRemoval(&(front_part), fp,	nrOfJunks, nrOfSamplesPerJunk );
	}

	assert( front_part.size == alpha );
}

void matingSelection()
/**
 * Select parents individuals uniformly from pp_all and add them to pp_sel
 */
{
	int i;

	pp_sel->size = mu;
	for (i = 0; i < mu; i++)
	{
		addToSelection( i, irand( pp_all->size ) );
	}
}

void addToSelection(int i, int c )
/**
 * adds the cth individual of pp_all to the ith place of pp_sel
 *
 * @param[in] c index of pp_all to be added
 * @param[in] i index of pp_sel to be replaced
 * @pre 0 <= i < pp_sel->size
 * @pre 0 <= c <= pp_all->size
 */
{
	assert( 0 <= i );
	assert( i < pp_sel->size );
	assert( 0 <= c );
	assert( c < pp_all->size );
	pp_sel->ind_array[i] = pp_all->ind_array[c];
}


void select_initial()
/* Performs initial selection. */
{
	selection();
}


void select_normal()
/* Performs normal selection.*/
{
	selection();
}

int irand(int range)
/* Generate a random integer. */
{
	int j;
	j=(int) ((double)range * (double) rand() / (RAND_MAX+1.0));
	return (j);
}

double drand( double from, double to )
{
	double j;
	j = from + (double)( (to-from)*rand() / ( RAND_MAX + 1.0) );
	return (j);
}

/*--------------------| data exchange functions |------------------------*/

int read_ini()
{
	int i;
	pp_new = create_pop(alpha, dim);

	for (i = 0; i < alpha; i++)
		pp_new->ind_array[i] = create_ind(dim);
	pp_new->size = alpha;

	return (read_pop(inifile, pp_new, alpha, dim));
}


int read_var()
{
	int i;
	pp_new = create_pop(lambda, dim);

	for (i = 0; i < lambda; i++)
		pp_new->ind_array[i] = create_ind(dim);

	pp_new->size = lambda;
	return (read_pop(varfile, pp_new, lambda, dim));
}


void write_sel()
{
	assert( pp_sel->size == mu );
	write_pop(selfile, pp_sel, mu);
}


void write_arc()
{
	assert( pp_all->size == alpha );
	write_pop(arcfile, pp_all, pp_all->size);
}


int check_sel()
{
	return (check_file(selfile));
}


int check_arc()
{
	return (check_file(arcfile));
}

void chk_free( void* handle )
{
	if( handle != NULL ) {
		free( handle );
	}
}