/* 
 *  FINDINDICES.C:
 *  MEX-File version of findhandles.m
 *
 *  Gary Levenson 
 *  Copyright (c) 1984-98 by The MathWorks, Inc.
 *  $Revision: 1.5 $ $Date: 1998/06/03 12:25:38 $
 */

#include "mex.h"

/*
 * syntax for M-File:
 *   Match = findhandles(NewHandleVect, OldHandleVect, FirstSecond)
 * 
 * purpose:
 *   Find all the instances of handles of New in Old and return indeces.
 *   FirstSecond determines whether to search over rows or columns
 *
 * function Match=findhandles(NewHandleVect,OldHandleVect,FirstSecond)
 * 
 * if isempty(OldHandleVect)|isempty(NewHandleVect),
 *  Match=[];  
 * else,  
 *  OldHandleVect=reshape(OldHandleVect,prod(size(OldHandleVect)),1);
 *  NewHandleVect=reshape(NewHandleVect,1,prod(size(NewHandleVect)));
 * 
 *  OldHandleMat=OldHandleVect(:,ones(size(NewHandleVect,2),1));
 *  NewHandleMat=NewHandleVect(ones(size(OldHandleVect,1),1),:);
 *  MatchMat= OldHandleMat==NewHandleMat;
 *  [Row,Column]=find(MatchMat);
 *  Match=(1:length(NewHandleVect));
 *  if FirstSecond==1,  
 *    Match=Column;    
 *  else,
 *    Match=Row;    
 *  end    
 * 
 * end % if isempty  
 * 
 */

#define NEW_HANDLE_VECT prhs[0]
#define OLD_HANDLE_VECT prhs[1]
#define ROW_OR_COL      prhs[2]
#define MATCH           plhs[0]

#define isMatrix(A)     (mxGetM(A)>1 && mxGetN(A)>1)

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

void mexFunction(
		 int           nlhs,
		 mxArray       *plhs[],
		 int           nrhs,
		 const mxArray *prhs[]
		 )
{
    int i, j, m1, m2, n1, n2, rowOrCol;
    int maxDim1, maxDim2, length = 0;
    double *newHandles, *oldHandles, *match, *output;
    bool empty;

#if 0
    /* Check for proper number of input and output args */

    if (nrhs != 3) {
	mexErrMsgTxt("FINDINDICES requires three input arguments");
    } else if (nlhs > 1) {
	mexErrMsgTxt("FINDINDICES requires one output arguments");
    }

    if (!mxIsNumeric(NEW_HANDLE_VECT) || mxIsComplex(NEW_HANDLE_VECT) ||
	mxIsSparse(NEW_HANDLE_VECT)  || !mxIsDouble(NEW_HANDLE_VECT) ||
	isMatrix(NEW_HANDLE_VECT)) {
	mexErrMsgTxt("NewHandleVect must be a real vector");
    }
    if (!mxIsNumeric(OLD_HANDLE_VECT) || mxIsComplex(OLD_HANDLE_VECT) ||
	mxIsSparse(OLD_HANDLE_VECT)  || !mxIsDouble(OLD_HANDLE_VECT) || 
	isMatrix(OLD_HANDLE_VECT)) {
	mexErrMsgTxt("OldHandleVect must be a real vector");
    }
    if (!mxIsNumeric(ROW_OR_COL) || mxIsComplex(ROW_OR_COL) ||
	mxIsSparse(ROW_OR_COL)   || !mxIsDouble(ROW_OR_COL) || 
	mxGetM(ROW_OR_COL) > 1   || mxGetN(ROW_OR_COL) > 1) {
	mexErrMsgTxt("rowOrCol must be a scalar");
    }
#endif
    
    m1 = mxGetM(NEW_HANDLE_VECT);
    n1 = mxGetN(NEW_HANDLE_VECT);
    maxDim1 = max(m1,n1);
    
    m2 = mxGetM(OLD_HANDLE_VECT);
    n2 = mxGetN(OLD_HANDLE_VECT);
    maxDim2 = max(m2,n2);
    
    empty = (m1 == 0 || n1 == 0 || m2 == 0 || n2 == 0);

    if (!empty) {
	
	newHandles  = mxGetPr(NEW_HANDLE_VECT);
	oldHandles  = mxGetPr(OLD_HANDLE_VECT);
	rowOrCol    = (int)mxGetScalar(ROW_OR_COL);
	output      = mxCalloc(max(maxDim1,maxDim2), sizeof(double));
	
	/*
	 * If rowOrCol is one, find the instances of new in old
	 */
	if (rowOrCol==1) {
	    for (j = 0; j<maxDim2; j++) {
		for (i = 0; i<maxDim1; i++) {
		    if (newHandles[i] == oldHandles[j]) {
			*output++ = i+1;
			length++;
		    }
		}
	    }
	} else {
	    for (j = 0; j<maxDim1; j++) {
		for (i = 0; i<maxDim2; i++) {
		    if (oldHandles[i] == newHandles[j]) {
			*output++ = i+1;
			length++;
		    }
		}
	    }
	}
	output -= length;   
	
	MATCH  = mxCreateDoubleMatrix(length,1,0);
	mxFree(mxGetPr(MATCH));
	mxSetPr(MATCH, output);
	return;
    } else {
	/* at least one input is [] */
	MATCH  = mxCreateDoubleMatrix(length,1,0);
	return;
    }	
}