localpca.cpp

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

// LocalPCA.cpp

#include <algorithm>
#include <cmath>
#include "Random.h"
#include "LocalPCA.h"

namespace az
{

namespace alg
{

void LocalPCA::Train()
{
	unsigned int	c,	// cluster index
					m,	// row index
					n;	// col index
	unsigned int i;

	// initialize the PI matrix
	mvPI.resize(mNClu);
	for(c=0; c<mNClu; c++)
	{
		mvPI[c].resize(mDX);
		for(m=0; m<mDX; m++) mvPI[c][m].resize(mDX);
		for(m=0; m<mDX; m++) for(n=m; n<mDX; n++) mvPI[c][m][n] = mvPI[c][m][n] = (m==n) ? 1.0 : 0.0;
	}

	// initialize the center points
	std::vector<unsigned int> index(mNX);
	for(i=0; i<mNX; i++) index[i] = i;
	for(c=0; c<mNClu; c++)
	{
		std::swap(index[c], index[rnd::rand(c, mNX)]);
		mvMean[c] = mvX[index[c]];
	}

	// trainning
	unsigned int failupdate = 0;
	double	dis, mindis, pro, err;
	std::vector<double> meanold(mDX);
	while( mIter++ < mMaxIter && failupdate<mNClu )
	{
		// partition trainning data
		for(c=0; c<mNClu; c++) mvNo[c] = 0;	// number point in each cluster
		for(m=0; m<mNX; m++)
		{
			mindis		= 1.0E100;
			mvIndex[m]	= 0;				// to which cluster it belongs
			for(c=0; c<mNClu; c++)
				if( (dis = Distance(m,c))< mindis)
				{
					mindis		= dis;
					mvIndex[m]	= c;
				}
			mvNo[mvIndex[m]]++;
		}

		// update parameters
		failupdate	= 0;
		for(c=0; c<mNClu; c++) 
		{
			// save old mean
			meanold = mvMean[c];

			// no data assigned to cluster c
			if(mvNo[c] < 1)
			{
				mvMean[c] = mvX[rnd::rand((unsigned int)0, (unsigned int)mNX)];
				for(m=0; m<mDX; m++) 
				{
					mvPI[c][m][m] = 1.0;
					for(n=m+1; n<mDX; n++) mvPI[c][m][n] = mvPI[c][n][m] = 0.0;
				}
			}
			// only one data assigned to cluster c
			else if(mvNo[c] == 1)
			{
				for(m=0; m<mNX; m++) if(mvIndex[m] == c) break;
				mvMean[c] = mvX[m];
				for(m=0; m<mDX; m++) 
				{
					mvPI[c][m][m] = 1.0;
					for(n=m+1; n<mDX; n++) mvPI[c][m][n] = mvPI[c][n][m] = 0.0;
				}
			}
			// more than one data assigned to cluster c
			else if(mvNo[c] > 1)
			{
				// find data index
				index.resize(mvNo[c]); n=0;
				for(m=0; m<mNX; m++) if(mvIndex[m] == c) index[n++] = m;
				
				// calculate the eigenvalues and eigenvectors
				Eigen(mvMean[c], mvEigenvalue[c], mvEigenvector[c], index);

				// calculate pi
				for(m=0; m<mDX; m++) for(n=0; n<mDX; n++)
				{
					mvPI[c][m][n] = 0.0;
					//for(i=0; i<mDLat; i++) mvPI[c][m][n] += mvEigenvector[c][i][m]*mvEigenvector[c][i][n];
					for(i=mDLat; i<mDX; i++) mvPI[c][m][n] += mvEigenvector[c][i][m]*mvEigenvector[c][i][n];
				}
			}// else if(mvNo[c] > 1)

			// calculate the error
			err = 0.0;
			for(m=0; m<mDX; m++) err += (meanold[m] - mvMean[c][m])*(meanold[m] - mvMean[c][m]);
			failupdate += (sqrt(err)<mErrTol) ? 1:0;
		}// for(c=0; c<mNClu; c++)
	}// while( mIter++ < mMaxIter && failupdate<mNClu )

	// calculate the projection
	for(c=0; c<mNClu; c++) if(mvNo[c] > 1)
	{
		// copy data
		index.resize(mvNo[c]); n=0;
		for(m=0; m<mNX;   m++) if(mvIndex[m] == c) index[n++] = m;
		for(m=0; m<mDLat; m++) {mvProMin[c][m] = 1.0E100; mvProMax[c][m] = -1.0E100;}
		for(m=0; m<index.size(); m++) for(n=0; n<mDLat; n++)
		{
			pro = 0.0;
			for(i=0; i<mDX; i++) pro += (mvX[index[m]][i]-mvMean[c][i])*mvEigenvector[c][n][i];
			if(pro>mvProMax[c][n]) mvProMax[c][n] = pro;
			if(pro<mvProMin[c][n]) mvProMin[c][n] = pro;
		}
	}
}

double LocalPCA::Distance(unsigned int m, unsigned int c)
{
	unsigned int i,j;
	double dis = 0.0;
	std::vector<double> tmpdis(mDX);
	for(i=0; i<mDX; i++)
	{
		tmpdis[i] = 0.0;
		for(j=0; j<mDX; j++) tmpdis[i] += (mvX[m][j]-mvMean[c][j])*mvPI[c][j][i];
	}
	for(i=0; i<mDX; i++) dis += tmpdis[i]*(mvX[m][i]-mvMean[c][i]);
	return dis;

}

} //namespace alg

} //namespace az