/* Roger Jang, 20091205 */

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

#define	INPUT	prhs[0]		/* Input Arguments */
#define	OUTPUT	plhs[0]		/* Output Arguments */

/* Mean computation */
double vecNanMean(double *vec, int vecLen){
	int i, validCount=0;
	double total=0;
	for (i=0; i<vecLen; i++)
		if (!mxIsNaN(vec[i])){
			total+=vec[i];
			validCount++;
		}
	return(total/validCount);
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
	double *input, *output;
	int i, m, n;

	/* Check for proper number of arguments */
	if (nrhs<1) {
		char message[100];
		strcpy(message, mexFunctionName());
		strcat(message, " requires 1 input vector.\n");
		mexErrMsgTxt(message);
	}
	
	/* Assign pointers to the various parameters */
	input = mxGetPr(INPUT);
	m = mxGetM(INPUT);
	n = mxGetN(INPUT);

	/* Create a matrix for the return argument */
	if (m==1 | n==1)
		OUTPUT = mxCreateDoubleMatrix(1, 1, mxREAL);
	else
		OUTPUT = mxCreateDoubleMatrix(1, n, mxREAL);
	output = mxGetPr(OUTPUT);
	
	/* Do the actual computation */
	if (m==1 | n==1)	/* input is a vector */
		output[0]=vecNanMean(input, m*n);
	else			/* input is a matrix */
		for (i=0; i<n; i++)
			output[i]=vecNanMean(input+i*m, m);
}