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