pairwisesqrdistance.cpp

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

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

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

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

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

void mexFunction(
	int nlhs,	mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	double	*out, *mat1, *mat2, diff, squaredSum;
	int dim1, dataNum1, dim2, dataNum2, i, j, k;
	
	/* Check for proper number of arguments */
	if ((nrhs<1)||(nrhs>2)||(nlhs>1)) {
		char message[100];
		sprintf(message, "Usage:\n\tdistMat=%s(mat1, mat2)\n\tdistMat=%s(mat)", mexFunctionName(), mexFunctionName());
		mexErrMsgTxt(message);
	}

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

	/* 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 = nrhs==1? mat1:mxGetPr(MAT2);

	/* Do the actual computations */
	if (nrhs==1)	// Symmetric matrix
		for (i=0; i<dataNum1; i++)
			for (j=i+1; j<dataNum2; j++){
				squaredSum = 0;
				for (k=0; k<dim1; k++){
					diff = mat1[i*dim1+k]-mat2[j*dim2+k];
					squaredSum += diff*diff;
				}
				out[j*dataNum1+i] = out[i*dataNum1+j] = squaredSum;
			}
	else
		for (i=0; i<dataNum1; i++)
			for (j=0; j<dataNum2; j++) {
				squaredSum = 0;
				for (k=0; k<dim1; k++){
					diff = mat1[i*dim1+k]-mat2[j*dim2+k];
					squaredSum += diff*diff;
				}
				out[j*dataNum1+i] = squaredSum;
			}
}