/*========================================================================
  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 <stdlib.h>
#include <float.h>
#include <math.h>

/*Debugging*/
#include <stdio.h>

#include "prefrel.h"
#include "population.h"
#include "hypervolume.h"
#include "epsilon.h"
#include "dominance.h"
#include "utils.h"
#include "assert.h"

double equalDistributionMeasure(int *A, int sizea)
{
	int i, j, k;
	double  diversity;
	diversity = 0.0;

	double minDist1, minDist2;
	double dist;
	/* Calculate the mean minimal distance */
	double mindist = -1;

	/* Calculate Deviation */
	mindist = -1;
	if (sizea > 1)
	{
		for (i = 0; i < sizea; i++)
		{
			minDist1 = DBL_MAX;
			minDist2 = DBL_MAX;
			for (j = 0; j < sizea; j++)
			{
				if (i != j)
				{
					dist = 0;
					for (k = 0; k < getNrOfObjectives(); k++)
						dist += getObjectiveValueDifference( A[i],A[j],k)*
						getObjectiveValueDifference( A[i],A[j],k);
					if (minDist1 > dist)
					{
						minDist2 = minDist1;
						minDist1 = dist;
					}
					else if (minDist2 > dist)
						minDist2 = dist;
				}
			}
			minDist1 = sqrt(minDist1);
			minDist2 = sqrt(minDist2);

			diversity -= 1.0/(minDist1+0.0001) + 1.0/(minDist2+0.0001);
		}
	}
	return diversity;
}

double unaryTschebyscheff(int *A, int sizea )
/* The smaller v the better */
{
	int k;
	int i = 0;
	double w1;
	double wr;
	double v = 0.0;
	double pr;
	int nrOfObj = getNrOfObjectives();

	assert( nrOfObj > 1 );

	/* Iterate over all weight vectors */
	for( pr = 0; pr < nrOfObj; pr++)
	{
		for( w1 = 0.7; w1 <= 1; w1 += 0.01 )
		{
			wr = (1-w1)/(double)(nrOfObj - 1.0);
			double u = 0.0;
			/* Iterate over all individuals */
			for( i = 0; i < sizea; i++ )
			{
				double max = 0.0;
				double tmp;
				/* Calculate the utility function */
				for( k = 0; k < nrOfObj; k++ )
				{
					if( k == pr)
						tmp = w1*getObjectiveValue( A[i], k);
					else
						tmp = wr*getObjectiveValue( A[i], k);
					if( k == 0 || tmp > max )
						max = tmp;
				}
				if( i == 0 || max < u ) {
					u = max;
				}
			}
			v += u;
		}
	}
	return v;
}

double unaryIndicator(int type, int* A, int sizea, double bound, int scaling )
/* indicator values are to be maximized */
{
	double v;

	switch (type)
	{
	case 3:
		v = -1*unaryTschebyscheff(A,sizea);
		break;
	case 4:
		v = unaryHypervolumeIndicator(A, sizea, bound, scaling);
		break;
	case 11:
		v = -1*unaryEpsilonFrontRef(A,sizea);
		break;
	case 14:
		v = equalDistributionMeasure(A,sizea);
		break;
	case 30:
		v = -1*unaryWeightedEpsilonReferencePoint(A,sizea,0);
		break;
	case 31:
		v = -1*unaryWeightedEpsilonReferencePoint(A,sizea,1);
		break;
	case 0:
	default:
		v = 0.0;
		break;
	}
	return v;
}

int prefDomPlusUnaryIndicator(int indicatorType, int sorting, int ranking,
		fpart* partpA, fpart* partpB, double bound,
		int scaling)
/* SIDE EFFECTS: partpA and partpB may be temporarily rearranged */
{
	int i;
	int aBetter, bBetter;
	double indicatorValueA, indicatorValueB;

	if (sorting == 0)
	{
		performDominanceSorting(ranking, partpA);
		performDominanceSorting(ranking, partpB);
	}
	aBetter = 0;
	bBetter = 0;
	int indiff = 0;
	i = 0;
	while (i < partpA->fronts && i < partpB->fronts && !aBetter && !bBetter)
	{
		switch (dominanceCheckPop(partpA->front_array[i].members,
				partpA->front_array[i].size,
				partpB->front_array[i].members,
				partpB->front_array[i].size))
		{
		case a_better_b:
			aBetter = 1;
			break;
		case b_better_a:
			if( b_better_a )
				bBetter = 1;
			break;
		case incomparable:
			indicatorValueA = unaryIndicator(indicatorType, partpA->front_array[i].members,
					partpA->front_array[i].size, bound, scaling);
			indicatorValueB = unaryIndicator(indicatorType, partpB->front_array[i].members,
					partpB->front_array[i].size,  bound, scaling);
			if (indicatorValueA > indicatorValueB)
				aBetter = 1;
			else if (indicatorValueB > indicatorValueA)
				bBetter = 1;
			break;
		default:
			indiff = 1;
			break;
		}
		if (sorting == 0)
			i = partpA->fronts;
		else
			i++;
	}
	if (sorting == 0)
	{
		undoDominanceSorting(partpA, ranking);
		undoDominanceSorting(partpB, ranking);
	}
	if( aBetter == 1 && bBetter == 0 )
		return 1;
	if( indiff == 1 )
		return 1;
	return 0;
}

comp prefDomPlusIndicatorFront(int indicatorType,
		front* frontA, front* frontB,
		double bound,
		int scaling )
{
	double indicatorValueA;
	double indicatorValueB;


	comp test = incomparable;

	test = dominanceCheckPop(frontA->members, frontA->size,
			frontB->members, frontB->size);


	switch ( test )
	{
	case a_better_b:
		return a_better_b;
	case b_better_a:
		return b_better_a;
	case incomparable:
		if( indicatorType < 50 ) {
			indicatorValueA = unaryIndicator(indicatorType, frontA->members,
					frontA->size, bound, scaling);
			indicatorValueB = unaryIndicator(indicatorType, frontB->members,
					frontB->size, bound, scaling);
		}
		if (indicatorValueA > indicatorValueB)
			return a_better_b;
		else if (indicatorValueB > indicatorValueA)
			return b_better_a;
		break;
	default:
		break;
	}
	return indifferent;
}

int prefDomPlus( int prefrel_type, fpart* partA, fpart* partB, int bound, int scaling )
{
	int remainingType; /*???xx*/
	int paretoFrontType; /*?xx??*/
	int procedure; /*x????*/

	getReductionTypes(prefrel_type, &procedure, &paretoFrontType, &remainingType);


	comp status = incomparable;

	front* currFrontA;
	front* currFrontB;

	int currFrontIndexA =  0;
	int currFrontIndexB =  0;

	int step = 1;

	while( status != a_better_b && status != b_better_a &&
			currFrontIndexA < partA->fronts && currFrontIndexB < partB->fronts )
	{
		currFrontA = &(partA->front_array[ currFrontIndexA ]);
		currFrontB = &(partB->front_array[ currFrontIndexB ]);
		if( step == 1 ) {
			status = prefDomPlusIndicatorFront(paretoFrontType,
					currFrontA, currFrontB, bound, scaling );
			if( procedure == 2 && status != a_better_b && status != b_better_a )
			{
				status = prefDomPlusIndicatorFront(remainingType,
						currFrontA, currFrontB, bound, scaling );
			}
		}
		else {
			status = prefDomPlusIndicatorFront(remainingType,
					currFrontA, currFrontB, bound, scaling );
		}
		currFrontIndexA++;
		currFrontIndexB++;
		step++;
	}

	if( status == a_better_b || status == indifferent )
		return 1;
	else
		return 0;
}

int prefrel(fpart* partpA, fpart* partpB, int type, int sorting, int ranking,
		double bound, int scaling)
/** Returns 1 if A is better than B and 0 otherwise */
{
	int result;
	/* using preferences on populations */
	if( type >= 0 )
	{
		result = prefDomPlus( type, partpA, partpB, bound, scaling );
		return result;
	}
	type = -type;
	/* using preferences on fronts */
	result = prefDomPlusUnaryIndicator(type, sorting, ranking,
				partpA, partpB, bound, scaling);

	return result;
}