genmod.cpp

该算法是张青复与周爱民的基于RM-MEDA算法的程序

// 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
//constructor
ModelLocalPCAU::ModelLocalPCAU()
{
	mLatentDim	= 0;
	mMaxCluster	= 0;
}

//initialize the LPCA
void ModelLocalPCAU::Set(unsigned int latent, unsigned int cluster, unsigned int trainsteps, double extension)
{
	mLatentDim	= latent;
	mMaxCluster	= cluster;
	mTrainSteps	= trainsteps;
	mExtension	= extension;
}

//model-based generator
CPopulationMO& 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