Daniel@0
|
1 /* CREATED:2010-01-17 08:12:57 by Brian McFee <bmcfee@cs.ucsd.edu> */
|
Daniel@0
|
2 /* binarysearch.c
|
Daniel@0
|
3 *
|
Daniel@0
|
4 * binary search
|
Daniel@0
|
5 *
|
Daniel@0
|
6 * Compile:
|
Daniel@0
|
7 * mex -DNAN_EQUALS_ZERO binarysearch.c
|
Daniel@0
|
8 */
|
Daniel@0
|
9
|
Daniel@0
|
10 #include "mex.h"
|
Daniel@0
|
11
|
Daniel@0
|
12
|
Daniel@0
|
13 #if NAN_EQUALS_ZERO
|
Daniel@0
|
14 #define IsNonZero(d) ((d) != 0.0 || mxIsNan(d))
|
Daniel@0
|
15 #else
|
Daniel@0
|
16 #define IsNonZero(d) ((d) != 0.0)
|
Daniel@0
|
17 #endif
|
Daniel@0
|
18
|
Daniel@0
|
19 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
Daniel@0
|
20 {
|
Daniel@0
|
21 /* Declare variables */
|
Daniel@0
|
22 int nQuery, nData; /* number of elements */
|
Daniel@0
|
23 double *pQuery, *pData; /* input arrays */
|
Daniel@0
|
24 double *positions; /* output array(s) */
|
Daniel@0
|
25
|
Daniel@0
|
26 int i;
|
Daniel@0
|
27
|
Daniel@0
|
28 if (nrhs != 2) {
|
Daniel@0
|
29 mexErrMsgTxt("Two input arguments are required.");
|
Daniel@0
|
30 }
|
Daniel@0
|
31 if (nlhs > 1) {
|
Daniel@0
|
32 mexErrMsgTxt("Too many output arguments.");
|
Daniel@0
|
33 }
|
Daniel@0
|
34
|
Daniel@0
|
35 if (!(mxIsDouble(prhs[0])) || !(mxIsDouble(prhs[1]))) {
|
Daniel@0
|
36 mexErrMsgTxt("Input arrays must be of type double.");
|
Daniel@0
|
37 }
|
Daniel@0
|
38
|
Daniel@0
|
39 nQuery = mxGetNumberOfElements(prhs[0]);
|
Daniel@0
|
40 nData = mxGetNumberOfElements(prhs[1]);
|
Daniel@0
|
41 pQuery = (double *)mxGetPr(prhs[0]);
|
Daniel@0
|
42 pData = (double *)mxGetPr(prhs[1]);
|
Daniel@0
|
43
|
Daniel@0
|
44 plhs[0] = mxCreateDoubleMatrix(1, nQuery, mxREAL);
|
Daniel@0
|
45 positions = mxGetPr(plhs[0]);
|
Daniel@0
|
46
|
Daniel@0
|
47 /* Now for the meat */
|
Daniel@0
|
48
|
Daniel@0
|
49 for (i = 0; i < nQuery; i++) {
|
Daniel@0
|
50 positions[i] = 0;
|
Daniel@0
|
51 positions[i] = binarysearch(pQuery[i], pData, nData);
|
Daniel@0
|
52 }
|
Daniel@0
|
53
|
Daniel@0
|
54 mxSetN(plhs[0], nQuery);
|
Daniel@0
|
55 }
|
Daniel@0
|
56
|
Daniel@0
|
57 int binarysearch(double q, double *data, int n) {
|
Daniel@0
|
58 int l = 0;
|
Daniel@0
|
59 int u = n-1;
|
Daniel@0
|
60 int pivot;
|
Daniel@0
|
61
|
Daniel@0
|
62 while (l < u) {
|
Daniel@0
|
63 pivot = l + (u - l) / 2;
|
Daniel@0
|
64
|
Daniel@0
|
65 if (q > data[pivot]) {
|
Daniel@0
|
66 u = pivot - 1;
|
Daniel@0
|
67 } else {
|
Daniel@0
|
68 l = pivot + 1;
|
Daniel@0
|
69 }
|
Daniel@0
|
70 }
|
Daniel@0
|
71
|
Daniel@0
|
72 /* Break ties to the right */
|
Daniel@0
|
73 if (l == n || q > data[l]) {
|
Daniel@0
|
74 return l;
|
Daniel@0
|
75 } else {
|
Daniel@0
|
76 return l + 1;
|
Daniel@0
|
77 }
|
Daniel@0
|
78
|
Daniel@0
|
79 }
|