lcsmex.c

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

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* Input Arguments */
#define	IN1	prhs[0]
#define	IN2	prhs[1]

/* Output Arguments */
#define	OUT1	plhs[0]
#define	OUT2	plhs[1]

#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif

static int lcs(int m, double *a, int n, double *b, double *result) {
	int **table;
	int count, i, j, temp;

	/* Memory allocation */
	table = (int **) malloc(sizeof(int)*(m+1));
	table[0] = (int *) malloc(sizeof(int)*(m+1)*(n+1));
	for (i=1; i<=m; i++)
		table[i] = table[i-1]+n+1;

	/* Set initial conditions for dynamic programming table */
	for (i=0; i<=m; table[i][0]=0, i++);
	for (j=0; j<=n; table[0][j]=0, j++);

	/* Dynamic programming */
	for (i=1; i<=m; i++)
		for (j=1; j<=n; j++)
			if (a[i-1]==b[j-1])
				table[i][j] = table[i-1][j-1]+1;
			else if (table[i][j-1] > table[i-1][j])
				table[i][j] = table[i][j-1];
			else
				table[i][j] = table[i-1][j];

	/* Print the DP table */
	/*
	for (i=0; i<=m; i++) {
		for (j=0; j<=n; j++)
			printf("%d ", table[i][j]);
		printf("\n");
	}
	*/

	/* Put the result */
	count = table[m][n];
	temp = count;
	for (i=m, j=n; (i!=0)&&(j!=0); )
		if (table[i][j] == table[i-1][j])
			i--;
		else if (table[i][j] == table[i][j-1])
			j--;
		else {
			result[--temp] = a[i-1];
			i--; j--;
		}

	/* Free memory */
	free(table[0]); free(table);

	return(count);
}

void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray*prhs[] ) { 
	double *in1, *in2, *out1, *out2, *result;
	unsigned int m, n, i;
    
	/* Check for proper number of arguments */
	if (nrhs != 2)
		mexErrMsgTxt("Two input arguments required."); 

	m = mxGetM(IN1)*mxGetN(IN1);
	n = mxGetM(IN2)*mxGetN(IN2);
	if (mxIsComplex(IN1) || mxIsComplex(IN2))
		mexErrMsgTxt("Inputs must be real."); 
    
	/* Create a matrix for the return argument */ 
	OUT1 = mxCreateDoubleMatrix(1, 1, mxREAL); 
    
	/* Assign pointers to the various parameters */ 
	out1 = mxGetPr(OUT1);
	in1 = mxGetPr(IN1); 
	in2 = mxGetPr(IN2);

	result = (double *) malloc(sizeof(double)*MAX(m, n));
	
	/* Return the first output, length of LCS. */
	*out1 = (double) lcs(m, in1, n, in2, result);
	/* Return the second output, LCS vector. */
	if (nlhs > 1) {
		OUT2 = mxCreateDoubleMatrix(1, (int) *out1, mxREAL);
		out2 = mxGetPr(OUT2);
		for (i=0; i<*out1; i++)
			out2[i] = result[i];
	}
	free(result);

	return;
}