wolffd@0: /* CREATED:2010-01-17 08:12:57 by Brian McFee <bmcfee@cs.ucsd.edu> */
wolffd@0: /* binarysearch.c
wolffd@0:  *
wolffd@0:  *  binary search
wolffd@0:  *
wolffd@0:  * Compile: 
wolffd@0:  *  mex binarysearch.c
wolffd@0:  */
wolffd@0: 
wolffd@0: #include "mex.h"
wolffd@0: 
wolffd@0: 
wolffd@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
wolffd@0: {
wolffd@0:     /* Declare variables */
wolffd@0:     int nQuery, nData;          /* number of elements */
wolffd@0:     double *pQuery, *pData;     /* input arrays */
wolffd@0:     double *positions;          /* output array(s) */
wolffd@0: 
wolffd@0:     int i;
wolffd@0: 
wolffd@0:     if (nrhs != 2) {
wolffd@0:         mexErrMsgTxt("Two input arguments are required.");
wolffd@0:     }
wolffd@0:     if (nlhs > 1) {
wolffd@0:         mexErrMsgTxt("Too many output arguments.");
wolffd@0:     }
wolffd@0: 
wolffd@0:     if (!(mxIsDouble(prhs[0])) || !(mxIsDouble(prhs[1]))) {
wolffd@0:         mexErrMsgTxt("Input arrays must be of type double.");
wolffd@0:     }
wolffd@0: 
wolffd@0:     nQuery      = mxGetNumberOfElements(prhs[0]);
wolffd@0:     nData       = mxGetNumberOfElements(prhs[1]);
wolffd@0:     pQuery      = (double *)mxGetPr(prhs[0]);
wolffd@0:     pData       = (double *)mxGetPr(prhs[1]);
wolffd@0:     
wolffd@0:     plhs[0]     = mxCreateDoubleMatrix(1, nQuery, mxREAL);
wolffd@0:     positions   = mxGetPr(plhs[0]);
wolffd@0: 
wolffd@0:     /* Now for the meat */
wolffd@0: 
wolffd@0:     for (i = 0; i < nQuery; i++) {
wolffd@0:         positions[i] = 0;
wolffd@0:         positions[i] = binarysearch(pQuery[i], pData, nData);
wolffd@0:     }
wolffd@0: 
wolffd@0:     mxSetN(plhs[0], nQuery);
wolffd@0: }
wolffd@0: 
wolffd@0: int binarysearch(double q, double *data, int n) {
wolffd@0:     int l = 0;
wolffd@0:     int u = n-1;
wolffd@0:     int pivot;
wolffd@0: 
wolffd@0:     while (l < u) {
wolffd@0:         pivot = l + (u - l) / 2;
wolffd@0: 
wolffd@0:         if (q > data[pivot]) {
wolffd@0:             u = pivot - 1;
wolffd@0:         } else {
wolffd@0:             l = pivot + 1;
wolffd@0:         }
wolffd@0:     }
wolffd@0: 
wolffd@0:     /* Break ties to the right */
wolffd@0:     if (l == n || q > data[l]) {
wolffd@0:         return l;
wolffd@0:     } else {
wolffd@0:         return l + 1;
wolffd@0:     }
wolffd@0: 
wolffd@0: }