pairwisedistance.cpp

来自「一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有」· C++ 代码 · 共 87 行

/* Pairwise squared Euclidean distance */
/* Roger Jang, 20030608 */

// How to compile:
// MATLAB 7.1: mex pairwiseDistance.cpp d:/users/jang/c/lib/utility.cpp -output pairwiseDistance.dll
// Others: mex pairwiseDistance.cpp d:/users/jang/c/lib/utility.cpp

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

/* Input Arguments */
#define	MAT1	prhs[0]
#define	MAT2	prhs[1]
#define	METHOD	prhs[2]

/* Output Arguments */
#define	OUT	plhs[0]

void mexFunction(
	int nlhs,	mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	double	*out, *mat1, *mat2;
	int dim, dataNum1, dim2, dataNum2, i, j, k, method, isOneMat;
	
	/* Check for proper number of arguments */
	if ((nrhs<1)||(nrhs>3)||(nlhs>1)){
		char message[100];
		sprintf(message, "Usage:\n\tdistMat=%s(mat1, mat2, method)\n\tdistMat=%s(mat1, mat2)\n\tdistMat=%s(mat, method)\n\tdistMat=%s(mat)",
			mexFunctionName(), mexFunctionName(), mexFunctionName(), mexFunctionName());
		mexErrMsgTxt(message);
	}
	
	method=2;	// Default to L2 norm
	if (nrhs==3){
		method=(int)mxGetScalar(METHOD);
		isOneMat=0;
	} else if (nrhs==2){
		if (mxGetM(MAT2)*mxGetN(MAT2)==1){
			method=(int)mxGetScalar(MAT2);
			isOneMat=1;
		} else {
			isOneMat=0;
		}
	} else {
		isOneMat=1;
	}
	
	if(!mxIsDouble(MAT1) || mxIsComplex(MAT1))
		mexErrMsgTxt("Input MAT1 must be a double real matrix.");
	
//	printf("method=%d\n", method);
//	printf("isOneMat=%d\n", isOneMat);

	/* Check dimensions */
	dim = mxGetM(MAT1);
	dataNum1 = mxGetN(MAT1);
	dim2 = isOneMat==1? dim:mxGetM(MAT2);
	dataNum2 = isOneMat==1? dataNum1:mxGetN(MAT2);
	if (dim != dim2)
		mexErrMsgTxt("Matrix sizes mismatch!");

	/* Check method */
	if ((method<1)||(method>2))
		mexErrMsgTxt("This function only supports L1 and L2 norm!");
	
	/* Create a matrix for the return argument */
	OUT = mxCreateDoubleMatrix(dataNum1, dataNum2, mxREAL);

	/* Assign pointers to the various parameters */
	out = mxGetPr(OUT);
	mat1 = mxGetPr(MAT1);
	mat2 = isOneMat==1? mat1:mxGetPr(MAT2);

	/* Do the actual computations */
	if (isOneMat==1)	// Return a symmetric matrix
		for (i=0; i<dataNum1; i++)
			for (j=i+1; j<dataNum2; j++){
				out[j*dataNum1+i] = out[i*dataNum1+j] = vecDistance(mat1+i*dim, mat2+j*dim, dim, method);
			}
	else
		for (i=0; i<dataNum1; i++)
			for (j=0; j<dataNum2; j++) {
				out[j*dataNum1+i] = vecDistance(mat1+i*dim, mat2+j*dim, dim, method);
			}
}