dawn@0: /* 
dawn@0: * rdiff_inplace.c
dawn@0: * running difference function
dawn@0: * 
dawn@0: * Alain de Cheveigné, CNRS/Ircam	Dec 2001
dawn@0: * (c) 2001 CNRS
dawn@0: */
dawn@0: 
dawn@0: #include <math.h>
dawn@0: #include "mex.h"
dawn@0: 
dawn@0: /* #define MWRAP */
dawn@0: #include "mwrap_check.h"
dawn@0: 
dawn@0: /* Input Arguments */
dawn@0: #define X_IN	prhs[0]
dawn@0: #define Y_IN	prhs[1]
dawn@0: #define R_IN	prhs[2]
dawn@0: #define L_IN 	prhs[3]
dawn@0: #define N_IN	prhs[4]
dawn@0: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: static void rdiff_inplace(
dawn@0: 				double *xp,		/* input vector 1 */
dawn@0: 				double *yp,		/* input vector 2 */
dawn@0: 				double *rp,		/* df matrix (time X lag) */
dawn@0: 				double *lp,		/* lag matrix */
dawn@0: 				int N,			/* size of summation window */
dawn@0: 				int m, 			/* rows in corr matrix */
dawn@0: 				int n	 		/* columns in corr matrix */
dawn@0: 			) 
dawn@0: {
dawn@0: 	int j,k,i, xlag, ylag;
dawn@0: 	double z, d, *r, *x, *y, *bump;
dawn@0: 		
dawn@0: 	
dawn@0: 	for (j=0;j<n;j++) {									/* for each column (lag pair) */
dawn@0: 		xlag = GET(lp[2*j]);
dawn@0: 		ylag = GET(lp[2*j+1]);
dawn@0: 		if (N==1) {
dawn@0: 			r=&rp[j*m];
dawn@0: 			x=&xp[xlag];
dawn@0: 			y=&yp[ylag];
dawn@0: 			bump=r+m;
dawn@0: 			while(r<bump) {								/* for each row */
dawn@0: 				d = GET(*x)-GET(*y);
dawn@0: 				SET(*r) = d*d;	
dawn@0: 				r++; x++; y++;
dawn@0: 			}
dawn@0: 		} else {
dawn@0: 			for (k=0; k<m; k++) {						/* for each row */
dawn@0: 				z=0;
dawn@0: 				x=&xp[xlag+k*N];
dawn@0: 				y=&yp[ylag+k*N];
dawn@0: 				bump=x+N;
dawn@0: 				while(x<bump) {
dawn@0: 					d = GET(*x)-GET(*y);
dawn@0: 					z+=d*d;
dawn@0: 					x++; y++;
dawn@0: 				}
dawn@0: 				SET(rp[j*m+k]) = z;
dawn@0: 			}
dawn@0: 		}
dawn@0: 	}
dawn@0: 	
dawn@0: 	return;						
dawn@0: }
dawn@0: 
dawn@0: void mexFunction(
dawn@0: 				int nlhs, mxArray *plhs[],
dawn@0: 				int nrhs, const mxArray *prhs[]
dawn@0: 				)
dawn@0: {
dawn@0: 	double	*xp, *yp, *rp, *lp, *up, *np; 
dawn@0: 	int nx, mx, ny, my, nr, mr, nl, ml, xmax, ymax, j, N, xlag, ylag;
dawn@0: 	
dawn@0: 	/* Check for proper number of arguments */
dawn@0: 	if (nrhs < 4 || nrhs > 5) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE takes 4 or 5 input arguments");
dawn@0: 	}
dawn@0: 	
dawn@0: 	/* Check type of input */
dawn@0: 	if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || 
dawn@0: 		mxIsSparse(X_IN)	|| !mxIsDouble(X_IN) ) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: X should be doubles");
dawn@0: 	}
dawn@0: 	if (!mxIsNumeric(Y_IN) || mxIsComplex(Y_IN) || 
dawn@0: 		mxIsSparse(Y_IN)	|| !mxIsDouble(Y_IN) ) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: Y should be doubles");
dawn@0: 	}
dawn@0: 	if (!mxIsNumeric(R_IN) || mxIsComplex(R_IN) || 
dawn@0: 		mxIsSparse(R_IN)	|| !mxIsDouble(R_IN) ) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: R should be doubles");
dawn@0: 	}
dawn@0: 	if (nrhs==5) {
dawn@0: 		if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || 
dawn@0: 			mxIsSparse(N_IN)	 || !mxIsDouble(N_IN) ) {
dawn@0: 			mexErrMsgTxt("RDIFF_INPLACE: N should double");
dawn@0: 		}
dawn@0: 		np = mxGetPr(N_IN);
dawn@0: 		N=*np;
dawn@0: 	} else {
dawn@0: 		N=1;
dawn@0: 	}
dawn@0: 	mx=mxGetM(X_IN);		/* rows */
dawn@0: 	nx=mxGetN(X_IN);		/* columns */
dawn@0: 	my=mxGetM(Y_IN);		
dawn@0: 	ny=mxGetN(Y_IN);		
dawn@0: 	mr=mxGetM(R_IN);		
dawn@0: 	nr=mxGetN(R_IN);		
dawn@0: 	ml=mxGetM(L_IN);		
dawn@0: 	nl=mxGetN(L_IN);		
dawn@0: 
dawn@0: 	if (nx>1 || ny>1) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: X and Y should be column vectors");
dawn@0: 	}
dawn@0: 	if (nr !=nl) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: result and lag matrices should have same number of columns");
dawn@0: 	}
dawn@0: 	if ( ml != 2) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: lags should be a two row matrix");
dawn@0: 	}
dawn@0: 	if (!nx || !mx || !ny || !my || !nr || !mr || !nl || !ml) {
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: one of the input matrices is empty");
dawn@0: 	}	
dawn@0: 		
dawn@0: 	xp = mxGetPr(X_IN);
dawn@0: 	yp = mxGetPr(Y_IN);
dawn@0: 	rp = mxGetPr(R_IN);
dawn@0: 	lp = mxGetPr(L_IN);
dawn@0: 	
dawn@0: 	checkin_matrix((mxArray *) L_IN);
dawn@0: 			
dawn@0: 	/* find maximum lag for x and y and check that no lag is negative */
dawn@0: 	xmax=GET(lp[0]);
dawn@0: 	ymax=GET(lp[1]);
dawn@0: 	for (j=1;j<nl;j++) {
dawn@0: 		xlag = GET(lp[2*j]);
dawn@0: 		ylag = GET(lp[2*j+1]);
dawn@0: 		if ( xlag<0 || ylag<0 ) {
dawn@0: 			mexErrMsgTxt("RDIFF_INPLACE: LAGS should be non-negative");
dawn@0: 		}
dawn@0: 		if (xlag>xmax ) xmax=xlag;
dawn@0: 		if (ylag>ymax ) ymax=ylag;
dawn@0: 	}
dawn@0: 
dawn@0: 	if ( (N*mr+ymax)> my) {
dawn@0: 		mexPrintf("mr*N: %d, ymax: %d, my: %d\n", mr*N, ymax, my);
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: Y data too short");
dawn@0: 	}
dawn@0: 	if ( (N*mr+xmax)> mx) {
dawn@0: 		mexPrintf("mr*N: %d, xmax: %d, mx: %d\n", mr*N, xmax, mx);
dawn@0: 		mexErrMsgTxt("RDIFF_INPLACE: X data too short");
dawn@0: 	}
dawn@0: 	
dawn@0: 	checkin_matrix((mxArray *) X_IN);
dawn@0: 	if (yp!=xp) checkin_matrix((mxArray *) Y_IN);
dawn@0: 	checkin_matrix((mxArray *) R_IN);
dawn@0: 	
dawn@0: 	/* Do the actual computations in a subroutine */
dawn@0: 	rdiff_inplace(xp,yp,rp,lp,N,mr,nr);
dawn@0: 
dawn@0: 	checkout_matrix((mxArray *) X_IN);
dawn@0: 	if (yp!=xp) checkout_matrix((mxArray *) Y_IN); 
dawn@0: 	checkout_matrix((mxArray *) R_IN);
dawn@0: 	checkout_matrix((mxArray *) L_IN);
dawn@0: 	return;
dawn@0: }
dawn@0: