// GenMod.cpp

#include <ctime>#include <list>#include <vector>#include <cmath>#include <float.h>#include "LocalPCA.h"#include "HCSampler.h"#include "GenMod.h"

//#define SD_DIS#if defined(WIN32)    #define wxFinite(n) _finite(n)#elif defined(_LINUX)    #define wxFinite(n) finite(n)#else    #define wxFinite(n) ((n)==(n))#endif

namespace az
{
namespace mea
{
namespace gen
{
namespace mod
{

//constructor
ModelBase::ModelBase()
{
	mDataSize	= mDataDim	= 0;
	pData		= 0;
	mExtension	= 0.2;
}

//destructor
ModelBase::~ModelBase()
{
	Clear();
}
//clear data pool
void ModelBase::Clear()
{
	unsigned int i;
	if( pData != 0 )
	{
		for( i=0; i<mDataSize; i++ ) delete []pData[i];
		delete []pData;
		pData	  = 0;
		mDataSize = 0;
	}
}

// convert a real number to a nearest positive integer number 
unsigned int ModelBase::ROUND(double X) 
{ 
	if(X-(unsigned int)(X)<0.5 && (unsigned int)(X)>0 ) 
		return (unsigned int)(X); 
	else 
		return (unsigned int)(X)+1; 
}

// pertubation
double ModelBase::PERTUBATION()	
{ 
	double r;
	do{ r = 1.5 + rnd::gaussian()*0.5; }while(r<0.8 || r>2.2);
	return r;
}

//////////////////////////////////////////////
//Local PCA based EDA generator//constructorModelLocalPCAU::ModelLocalPCAU(){	mLatentDim	= 0;	mMaxCluster	= 0;}//initialize the LPCAvoid ModelLocalPCAU::Set(unsigned int latent, unsigned int cluster, unsigned int trainsteps, double extension){	mLatentDim	= latent;	mMaxCluster	= cluster;	mTrainSteps	= trainsteps;	mExtension	= extension;}//model-based generatorCPopulationMO& ModelLocalPCAU::Generate(CPopulationMO& popnew, CPopulationMO& popref){	unsigned int c,i,j,k;		//Step 1: assign new data	if( popref.Size() != mDataSize ) 	{		Clear();			mDataSize	= popref.Size();		mDataDim	= popref.P().XSize();		pData		= new double*[mDataSize];		for( i=0; i<mDataSize; i++ ) pData[i] = new double[mDataDim];	}	for( i=0; i<mDataSize; i++ ) for( j=0; j<mDataDim; j++ ) pData[i][j] = popref[i][j];		//Step 2: train with Local PCA	alg::LocalPCA lpca;	//alg::Kmeans lpca;	lpca.Set(mMaxCluster, mDataSize, mDataDim, mLatentDim, mTrainSteps);	for( i=0; i<mDataSize; i++ ) for( j=0; j<mDataDim; j++ ) lpca.mvX[i][j] = popref[i][j];	lpca.Train();	//Step 3: calculate the probability of each cluster, i.e the size to create in each cluster	unsigned int nt = 0;	std::vector<unsigned int> nc(mMaxCluster);	for(i=0; i<mMaxCluster; i++) { nc[i] = (lpca.mvNo[i] >= 1) ? 2:0; nt += nc[i];}	double vt = 0.0;	std::vector<double> vc(mMaxCluster);	for(i=0; i<mMaxCluster; i++) 	{ 		vc[i] = 0.0;		if(lpca.mvNo[i] >  1)		{			vc[i] = 1.0;			for(j=0; j<mLatentDim; j++) vc[i] *= lpca.mvProMax[i][j] - lpca.mvProMin[i][j]; 		}		vt += vc[i];	}	double ns = mDataSize - nt + 0.0;	nt	= 0;	for(i=0; i<mMaxCluster; i++) { nc[i] += (unsigned int)(ns*vc[i]/vt); nt += nc[i];}	while(nt<mDataSize)	{		i = rnd::rand((unsigned int)0, (unsigned int)mMaxCluster);		if(lpca.mvNo[i]>1) {nc[i]++; nt++;}	}#ifdef AZ_MODEL_OUT	std::ofstream fhand("model.set");	fhand<<"PCA"<<std::endl;	fhand<<mMaxCluster<<std::endl;#endif	//Step 4: create new trial solutions in each cluster	popnew.Resize(mDataSize);	unsigned int np=0; 	double sd;	for(c=0; c<mMaxCluster; c++)	{		// only one point in the cluster		if(lpca.mvNo[c] == 1)		{			sd  = 0.0;			for(i=0; i<mDataDim; i++) sd += popref.P().XUpp(i)-popref.P().XLow(i);			sd /= mDataDim;			for(i=0; i<nc[c]; i++) for( j=0; j<mDataDim; j++ ) popnew[np+i][j] = lpca.mvMean[c][j] + sd*rnd::gaussian();		}		// more than one point in the cluster		else if(lpca.mvNo[c] > 1)		{			std::vector< std::vector<double> > t(mLatentDim);			for(i=0; i<mLatentDim; i++) t[i].resize(nc[c]);			std::vector<double> low(mLatentDim),upp(mLatentDim);			for(i=0; i<mLatentDim; i++) 			{				low[i] = lpca.mvProMin[c][i]-mExtension*(lpca.mvProMax[c][i]-lpca.mvProMin[c][i]);				upp[i] = lpca.mvProMax[c][i]+mExtension*(lpca.mvProMax[c][i]-lpca.mvProMin[c][i]);			}			alg::LHC(t, low, upp);			sd = 0.0;			for(i=mLatentDim; i<mDataDim; i++) sd += lpca.mvEigenvalue[c][i];			sd = sqrt(sd / double(mDataDim-mLatentDim));						for(i=0; i<nc[c]; i++) for( j=0; j<mDataDim; j++ ) 			{				popnew[np+i][j] = lpca.mvMean[c][j] + sd*rnd::gaussian();				for(k=0; k<mLatentDim; k++) popnew[np+i][j] += t[k][i]*lpca.mvEigenvector[c][k][j];			}		}		np += nc[c];#ifdef AZ_MODEL_OUT		fhand<<lpca.mvProMin[c][0] * lpca.mvEigenvector[c][0][0] + lpca.mvMean[c][0] <<"\t"			 <<lpca.mvProMin[c][0] * lpca.mvEigenvector[c][0][1] + lpca.mvMean[c][1] <<"\t"			 <<lpca.mvProMax[c][0] * lpca.mvEigenvector[c][0][0] + lpca.mvMean[c][0] <<"\t"			 <<lpca.mvProMax[c][0] * lpca.mvEigenvector[c][0][1] + lpca.mvMean[c][1] <<"\t"			 <<sd<<std::endl;#endif	}#ifdef AZ_MODEL_OUT	fhand.close();#endif		// Step 7: check the range	for(i=0; i<mDataSize; i++)		for(j=0; j<mDataDim; j++)		{			if( wxFinite(popnew[i][j]) == 0 ) popnew[i][j] = rnd::rand(popnew.P().XLow(j),popnew.P().XUpp(j));			// border check strategy 1: DE strategy			if(popnew[i][j] < popnew.P().XLow(j))		popnew[i][j] = 0.5*(popnew.P().XLow(j)+popref[i][j]);			else if(popnew[i][j] > popnew.P().XUpp(j))	popnew[i][j] = 0.5*(popnew.P().XUpp(j)+popref[i][j]);	}	return popnew;}

} //namespace mod
} //namespace gen
} //namespace mea
} //namespace az