parzen_classif.c

Parzen窗估计法的C语言代码和Matlab应用的动态链接库。

/*


 parzen_classif : Returns the estimated Labels ytest [1 x Ntest] given Xtrain data [d x Ntrain] and train label ytrain [1 x Ntrain]
 

 Usage
 -------


 [ytest , densite]    = parzen_classif(Xtrain , ytrain , Xtest , [sigma] );

 Inputs
 -------

 Xtrain        Train data   (d x Ntrain)

 ytrain        Train labels (1 x Ntrain)

 Xtest         Test  data   (d x Ntest)

 sigma         Noise Standard deviation of the rbf (default sigma = 1.0)




 Ouputs
 -------

 ytest         Estimated labels (1 x Ntest)

 densite       Estimated density (m x Ntest) where m denotes the number of class

	
 To compile
 ----------

 mex -output parzen_classif.dll parzen_classif.c


 mex  -f mexopts_intelamd.bat -output parzen_classif.dll parzen_classif.c


 Example 1
 ---------

 d      = 5;
 Ntrain = 100;
 Ntest  = 20;
 Xtrain = randn(d , Ntrain);
 ytrain = double(rand(1 , Ntrain)>0.5);
 Xtest  = randn(d , Ntest);
 sigma  = 1;

 [ytest , densite]  = parzen_classif(Xtrain , ytrain , Xtest , sigma );


Example 2
---------

load wine

[d , N]    = size(X);
n1         = round(0.7*N);
n2         = N - n1;
nbite      = 1000;
sigma      = 1;
Perf       = zeros(1 , nbite);

for i=1:nbite

 ind        = randperm(length(y));

 ind1       = ind(1:n1);
 ind2       = ind(n1+1:N);

 Xtrain     = X(: , ind1);
 ytrain     = y(ind1);
 Xtest      = X(: , ind2);

 ytest      = y(ind2);
 ytest_est  = parzen_classif(Xtrain , ytrain , Xtest , sigma );

 Perf(i)    = sum(ytest == ytest_est)/n2;

end

 


  
Author : S閎astien PARIS : [e-mail]
-------


Reference  : Vincent, P. and Bengio, Y. (2003). Manifold parzen windows. In Becker, S., Thrun, S., and Obermayer, K., editors,
-----------  Advances in Neural Information Processing Systems 15, Cambridge, MA. MIT Press.


*/


#include <math.h>
#include "mex.h"


/*--------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------*/


void qs( double * , int , int  ); 


void parzen_classif(double * , double * , double * , double  , double * , double *, int  , int , int  , int , double * , double* );

/*--------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------*/




void mexFunction( int nlhs, mxArray *plhs[] , int nrhs, const mxArray *prhs[] )

{
	
	
	double *Xtrain , *ytrain , *Xtest;
	
	double *ytest, *densite;
		
	int  d , Ntrain , Ntest , m=0 , i , currentlabel;

	double sigma = 1.0;

	double *ytrainsorted , *labels;

	
	
	/*--------------------------------------------------------------------------------*/
	/*--------------------------------------------------------------------------------*/
	/* -------------------------- Parse INPUT  -------------------------------------- */
	/*--------------------------------------------------------------------------------*/	
	/*--------------------------------------------------------------------------------*/
	
		
	/* ----- Input 1 ----- */
	
	
	Xtrain           = mxGetPr(prhs[0]);
		
	
	d                = mxGetM(prhs[0]);
	
	Ntrain           = mxGetN(prhs[0]);


	/* ----- Input 2 ----- */
	
	
	ytrain           = mxGetPr(prhs[1]);
		

	if(mxGetN(prhs[1]) != Ntrain)
	{

		 mexErrMsgTxt("ytrain must be (1 x Ntrain)");	

	}
	

	/* ----- Input 3 ----- */

	
	Xtest            = mxGetPr(prhs[2]);


	if(mxGetM(prhs[2]) != d)
	{

		 mexErrMsgTxt("Xtest must be (d x Ntest)");	

	}

    	
	Ntest            = mxGetN(prhs[2]);

	
	/* ----- Input 4 ----- */


	if (nrhs > 3)

	{
		
		sigma       = (double)mxGetScalar(prhs[3]);
		
	}
	


	
	
	
	/* Determine unique Labels */
	
	ytrainsorted  = mxMalloc(Ntrain*sizeof(double));
	
	
	for ( i = 0 ; i < Ntrain; i++ ) 
	{
		
		ytrainsorted[i] = ytrain[i];
		
	}
	
	
	qs( ytrainsorted , 0 , Ntrain - 1 );
	
	
	labels       = mxMalloc(sizeof(double)); 
	
	labels[m]    = ytrainsorted[0];
	
	currentlabel = labels[0];
	
	for (i = 0 ; i < Ntrain ; i++) 
	{ 
		if (currentlabel != ytrainsorted[i]) 
		{ 
			labels       = (double *)mxRealloc(labels , (m+2)*sizeof(double)); 
			
			labels[++m]  = ytrainsorted[i]; 
			
			currentlabel = ytrainsorted[i];
			
		} 
	} 
	
	m++; 


	/*--------------------------------------------------------------------------------*/
	/*---------------------------------------,----------------------------------------*/
	/* -------------------------- Parse OUTPUT  ------------------------------------- */
	/*--------------------------------------------------------------------------------*/
	/*--------------------------------------------------------------------------------*/
	
	/* ----- output 1 ----- */
	
	
	plhs[0]        = mxCreateDoubleMatrix(1 , Ntest , mxREAL);

	ytest          = mxGetPr(plhs[0]);

	plhs[1]        = mxCreateDoubleMatrix(m , Ntest , mxREAL);

	densite        = mxGetPr(plhs[1]);
	
	
	/*---------------------------------------------------------------------------------*/
	/*---------------------------------------------------------------------------------*/
	/* ----------------------- MAIN CALL  -------------------------------------------- */
	/*---------------------------------------------------------------------------------*/
	/*---------------------------------------------------------------------------------*/	
	/*---------------------------------------------------------------------------------*/
	
	
	parzen_classif(Xtrain , ytrain , Xtest , sigma , 
		           ytest , densite ,
				   d , m , Ntrain , Ntest , labels , ytrainsorted);


	mxFree(labels);
	
	mxFree(ytrainsorted);

	
	
	/*-----------------------------------------------*/
	/*-----------------------------------------------*/
	/* ------------ END of Mex File ---------------- */
	/*-----------------------------------------------*/
	/*-----------------------------------------------*/
		
}

/*----------------------------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------------------------*/



void parzen_classif(double *Xtrain , double *ytrain , double *Xtest , double sigma , 
					double *ytest , double *densite,
					int d , int m, int Ntrain , int Ntest , double *labels , double *ytrainsorted)
						
{
	int i , j , l , t , id , jd , im , ind;
	
	double temp , res  , maxi;
	
	double cte = -0.5/(sigma*sigma), cteprior , epsilon = 10e-50 ;
		
	
	double  *prior;
	
	
	
	/*             Prior            */
	
	prior    =  mxMalloc(m*sizeof(double)); 

	for (t = 0 ; t < m ; t++)
	{
		
		prior[t]    = 0.0;
		
	}
	
	for (i = 0 ; i < Ntrain ; i++) 
	{ 
		for (t = 0 ; t < m ; t++)
		{
			
			if (labels[t] == ytrainsorted[i])
				
			{
				
				prior[t]++;
				
			}
				
		}
	}

		
	for (t = 0 ; t < m ; t++)
	{
		
		prior[t]    /= Ntrain;
		
	}


	
	/* Classify */

		
	for (i = 0; i < Ntest ; i++) 
		
	{
		
		id    = i*d;

		im    = i*m;

		for (j = 0 ; j < Ntrain ; j++) 
			
		{
			jd  = j*d;
			
			
            for(t = 0 ; t < m ; t++)
			{
				if(ytrain[j] == labels[t])
					
				{
					
					ind = t;
				} 
				
			}
			
			res = 0.0;
			
			for (l = 0 ; l < d ; l++)
				
			{
				
				temp  = (Xtest[l + id] - Xtrain[l + jd]);
				
				res  += (temp*temp);
			}
			
			densite[ind + im] += exp(cte*res) + epsilon;
			
		}

		        
		cteprior = 0.0;

		for(t = 0 ; t < m ; t++)
			
		{
			densite[t + im] *=prior[t];

			cteprior        += densite[t + im];

		}

		cteprior = 1.0/cteprior;

		ind      = 0;
		
		maxi     = 0.0;
		
		for(t = 0 ; t < m ; t++)
			
		{
			temp            = densite[t + im]*cteprior;

			densite[t + im] = temp;
			
			if(temp > maxi)
			{
				
				maxi = temp;
				
				ind  = t;
			}
			
			
		}
		
		
		ytest[i] = labels[ind];
		
	}
	

	mxFree(prior);

	
}


/*-------------------------------------------------------------------------------------------------*/


void qs( double *array, int left, int right ) 
{
	
	double pivot;	// pivot element.
	
	int holex;	// hole index.
	
	int i;
	
	holex          = left + ( right - left )/2;

	pivot          = array[ holex ];		     // get pivot from middle of array.
	
	array[holex]   = array[ left ];              // move "hole" to beginning of
	
	holex          = left;			             // range we are sorting.
	
	for ( i = left + 1 ; i <= right ; i++ ) 
	{
		if ( array[ i ] <= pivot ) 
		{
			array[ holex ] = array[ i ];

			array[ i ]     = array[ ++holex ];
		}
	}
	
	if ( holex - left > 1 ) 
	{
	
		qs( array, left, holex - 1 );
	
	}
	if ( right - holex > 1 ) 
	{
	
		qs( array, holex + 1, right );
	
	}
	
	array[ holex ] = pivot;
	
}