/*========================================================================
  PISA  (www.tik.ee.ethz.ch/pisa/)
  ========================================================================
  Computer Engineering (TIK)
  ETH Zurich
  ========================================================================
  SPAM - Set Preference Algorithm for Multiobjective Optimization

  authors: Johannes Bader, johannes.bader@tik.ee.ethz.ch
           Eckart Zitzler, eckart.zitzler@tik.ee.ethz.ch
           Marco Laumanns, laumanns@tik.ee.ethz.ch

  revision by: Stefan Bleuler, stefan.bleuler@tik.ee.ethz.ch
			   Dimo Brockhoff, dimo.brockhoff@tik.ee.ethz.ch

  last change: 01.12.2008
  ========================================================================
 */

#include <assert.h>

#include "hypervolume.h"
#include "utils.h"
#include "population.h"
#include "dominance.h"
#include "math.h"

double hypervolume(double  *refs, int dim)
/* Usual hypervolume */
{
	int  i;
	double  volume;

	/* calculate hypervolume */
	volume = 1.0;
	for (i = 0; i < dim; i ++)
	{
		volume *= refs[dim + i] - refs[i];
	}

	return volume;
}

double compute_integral(double *refs, double bound, int dim )
{
	double volume=0;

	volume = hypervolume(refs, dim);
	return volume;
}

double calc_hypercube_integral(double  *refs, int dim, int scaling, double bound)
{
	double  volume;
	/* factor, the WDF near front is larger then maximal hypervolume on entire
		search space when scaling == 0 */
	double factor;
	double *scaledrefs;
	int i,j;
	int scalingnecessary;

	scaledrefs = chk_malloc(sizeof(double) * dim * 2);
	/* initialize factor: */
	factor = 1.0;
	for (i=0; i<dim; i++)
	{
		factor *= bound;
	}

	/* get deep copy of hypercube's reference points, the integral
		of which is sought */
	for (i=0; i<2*dim; i++)
	{
		scaledrefs[i] = refs[i];
	}

	/* add hypervolume on entire search space */
	volume = hypervolume(scaledrefs, dim);

	/* if scaling == 0:
			- compute usual hypervolume of hypercube ( done above )
			- add integral of additional weight distribution function
				near Pareto-optimal front, here in [0,1] x ... x [0,1]
	*/
	if (scaling == 0)
	{

		/* scale points for additional kernel near Pareto-optimal front
		   ------------------------------------------------------------ */

		/* scaling is necessary if second point lies in hypercube
			[(bound-1),bound] x [(bound-1),bound] x ... x [(bound-1)-bound] */
		scalingnecessary = 1;
		for (i=0; i<dim; i++)
		{
			if (scaledrefs[dim+i] < (bound-1.0))
			{
				scalingnecessary = 0;
				break;
			}
		}
		if (scalingnecessary)
		{
			/* if scaling is necessary, translate reference points
				to [0,1] x [0,1] x ... x [0,1] and cut - if necessary -
				on bounding box of this hypercube */
			for (i=0; i<2*dim; i++)
			{
				scaledrefs[i] = scaledrefs[i] - (bound-1.0);
			}
			for (i=0; i<dim; i++)
			{
				if (scaledrefs[i] < 0)
				{
					scaledrefs[i] = 0;
				}
			}
		}
		else
		{
			for (i=0; i<2*dim; i++)
			{
				scaledrefs[i] = 0;
			}
		}

	/* compute integral of additional weight distribution function
		near Pareto-optimal front */
		volume += factor * compute_integral(scaledrefs, bound, dim);

	}

	/* if scaling == 3:
			- compute usual indicator for hypercube, already scaled
				to [0,1] x [0,1] ( done above )
			- add line integrals of constant weight distribution functions
				on the axes
	*/
	if (scaling == 3)
	{
		/* add weight distibution function on axis */
		for (i = 0; i < dim; i ++)
		{
			if (scaledrefs[i] == 0.0)
			{
				double  tempVol = 1.0;
				for (j = 0; j < dim; j++)
					if (j != i)
						tempVol *= scaledrefs[dim + j] - scaledrefs[j];
				volume += tempVol;
			}
		}
	}

	chk_free( scaledrefs ); scaledrefs = NULL;

	return volume;
}
void swap(double  *front, int dim, 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 filter_nondominated_set(double  *front, int dim, 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 (dominates(&(front[i * dim]), &(front[j * dim]),
			              no_objectives))
			{
				/* remove point 'j' */
				n--;
				swap(front, dim, j, n);
			}
			else if (dominates(&(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, dim, i, n);
				i--;
				break;
			}
			else
				j++;
		}
		i++;
	}
	return n;
}

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

	assert(no_points > 0);
	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 dim, 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, dim, i, n);
		}

	return n;
}

double  calc_hypervolume_cubewise(double  *front, double  *refs,
                                  int  no_points, int  no_objectives, int dim,
                                  double bound, int scaling)
{
	int     n;
	double  volume, distance;

	volume = 0;
	distance = 0;
	n = no_points;
	refs[no_objectives - 1] = 0.0;
	while (n > 0)
	{
		int     no_nondominated_points;
		double  temp_vol, temp_dist;

		if (no_objectives > 1)
			no_nondominated_points = filter_nondominated_set(front, dim, n,
			                         no_objectives - 1);
		else
			no_nondominated_points = filter_nondominated_set(front, dim, n,
			                         no_objectives);

		assert(no_nondominated_points <= n);
		assert(no_nondominated_points > 0);

		temp_dist = surface_unchanged_to(front, dim, n, no_objectives - 1);
		refs[dim + no_objectives - 1] = temp_dist;

		temp_vol = 0;
		if (no_objectives < 3)
		{
			assert(no_nondominated_points > 0);
			refs[0] = 0.0;
			refs[dim] = front[0];
			temp_vol = calc_hypercube_integral(refs, dim, scaling,  bound);
		}
		else
			temp_vol = calc_hypervolume_cubewise(front, refs, no_nondominated_points,
			                                     no_objectives - 1, dim, bound, scaling);
		volume += temp_vol;

		refs[no_objectives - 1] = temp_dist;
		distance = temp_dist;
		n = reduce_nondominated_set(front, dim, n, no_objectives - 1, distance);
	}

	return volume;
}

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

	assert(no_objectives > 0);
	volume = 0;
	distance = 0;
	n = no_points;
	if (no_objectives == 1) {
		filter_nondominated_set(front, dim, n, 1);
		volume = front[0];
	}
	else {
		 while (n > 0)
		 {
			 int     no_nondominated_points;
			 double  temp_vol, temp_dist;

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

double unaryHypervolumeIndicator(int  *A, int sizea, double bound, int scaling )
{
	int dim = getNrOfObjectives();
	double pointArray[ sizea*dim ];
	double refPoints[ dim*2 ];
	int i, k;
	for (i = 0; i < sizea; i++)
		for (k = 0; k < dim; k++)
			pointArray[i * dim + k] =
			    getReversedObjectiveValue( A[i], k, scaling, bound );

	return calc_hypervolume_cubewise(pointArray, refPoints, sizea, dim, dim, bound, scaling);
}