cleaned up; ready for release
[RBC.git] / driver.cu
1 /* This file is part of the Random Ball Cover (RBC) library.
2 * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
3 */
4
5 #include<stdio.h>
6 #include<stdlib.h>
7 #include<cuda.h>
8 #include<sys/time.h>
9 #include<math.h>
10 #include "defs.h"
11 #include "utils.h"
12 #include "utilsGPU.h"
13 #include "rbc.h"
14 #include "brute.h"
15 #include "sKernel.h"
16
17 void parseInput(int,char**);
18 void readData(char*,unint,unint,real*);
19 void readDataText(char*,unint,unint,real*);
20 void orgData(real*,unint,unint,matrix,matrix);
21 void evalNNerror(matrix, matrix, unint*);
22 void evalKNNerror(matrix,matrix,intMatrix);
23
24 char *dataFile, *outFile;
25 unint n=0, m=0, d=0, numReps=0, s=0;
26 unint deviceNum=0;
27 int main(int argc, char**argv){
28 real *data;
29 matrix x, q;
30 unint *NNs;
31 intMatrix NNsK, kNNsRBC;
32 unint i;
33 struct timeval tvB,tvE;
34 cudaError_t cE;
35 rbcStruct rbcS;
36
37 printf("*****************\n");
38 printf("RANDOM BALL COVER\n");
39 printf("*****************\n");
40
41 parseInput(argc,argv);
42
43 cuInit(0);
44 printf("Using GPU #%d\n",deviceNum);
45 if(cudaSetDevice(deviceNum) != cudaSuccess){
46 printf("Unable to select device %d.. exiting. \n",deviceNum);
47 exit(1);
48 }
49
50 unsigned int memFree, memTot;
51 CUcontext pctx;
52 unsigned int flags=0;
53 int device;
54 cudaGetDevice(&device);
55 cuCtxCreate(&pctx,flags,device);
56 cuMemGetInfo(&memFree, &memTot);
57 printf("GPU memory free = %u/%u (MB) \n",memFree/(1024*1024),memTot/(1024*1024));
58
59 data = (real*)calloc( (n+m)*d, sizeof(*data) );
60 x.mat = (real*)calloc( PAD(n)*PAD(d), sizeof(*(x.mat)) );
61
62 //Need to allocate extra space, as each group of q will be padded later.
63 q.mat = (real*)calloc( PAD(m)*PAD(d), sizeof(*(q.mat)) );
64 x.r = n; x.c = d; x.pr = PAD(n); x.pc = PAD(d); x.ld = x.pc;
65 q.r = m; q.c = d; q.pr = PAD(m); q.pc = PAD(d); q.ld = q.pc;
66
67 NNs = (unint*)calloc( m, sizeof(*NNs) );
68 for(i=0; i<m; i++)
69 NNs[i]=DUMMY_IDX;
70
71 readData(dataFile, (n+m), d, data);
72 orgData(data, (n+m), d, x, q);
73 free(data);
74
75 NNsK.r=q.r; NNsK.pr=q.pr; NNsK.pc=NNsK.c=K; NNsK.ld=NNsK.pc;
76 kNNsRBC.r=q.r; kNNsRBC.pr=q.pr; kNNsRBC.pc=kNNsRBC.c=K; kNNsRBC.ld=kNNsRBC.pc;
77 kNNsRBC.mat = (unint*)calloc(kNNsRBC.pr*kNNsRBC.pc, sizeof(*kNNsRBC.mat));
78 NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
79
80 /* printf("running k-brute force..\n"); */
81 /* gettimeofday(&tvB,NULL); */
82 /* bruteK(x,q,NNsK); */
83 /* gettimeofday(&tvE,NULL); */
84 /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
85
86 printf("\nrunning rbc..\n");
87 gettimeofday(&tvB,NULL);
88 buildRBC(x, &rbcS, numReps, s);
89 gettimeofday(&tvE,NULL);
90 printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
91
92 //This finds the 32-NN; if you are only interested in the 1-NN, use queryRBC(..) instead
93 gettimeofday(&tvB,NULL);
94 kqueryRBC(q, rbcS, kNNsRBC);
95 gettimeofday(&tvE,NULL);
96 printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
97
98 destroyRBC(&rbcS);
99 printf("finished \n");
100
101 cE = cudaGetLastError();
102 if( cE != cudaSuccess ){
103 printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
104 }
105
106 evalKNNerror(x,q,kNNsRBC);
107
108 free(NNs);
109 free(NNsK.mat);
110 free(kNNsRBC.mat);
111 free(x.mat);
112 free(q.mat);
113 }
114
115
116 void parseInput(int argc, char **argv){
117 int i=1;
118 if(argc <= 1){
119 printf("\nusage: \n testRBC -f datafile (bin) -n numPts (DB) -m numQueries -d dim -r numReps -s numPtsPerRep [-o outFile] [-g GPU num]\n\n");
120 printf("\tdatafile = binary file containing the data\n");
121 printf("\tnumPts = size of database\n");
122 printf("\tnumQueries = number of queries\n");
123 printf("\tdim = dimensionailty\n");
124 printf("\tnumReps = number of representatives\n");
125 printf("\tnumPtsPerRep = number of points assigned to each representative\n");
126 printf("\toutFile = output file (optional); stored in text format\n");
127 printf("\tGPU num = ID # of the GPU to use (optional) for multi-GPU machines\n");
128 printf("\n\n");
129 exit(0);
130 }
131
132 while(i<argc){
133 if(!strcmp(argv[i], "-f"))
134 dataFile = argv[++i];
135 else if(!strcmp(argv[i], "-n"))
136 n = atoi(argv[++i]);
137 else if(!strcmp(argv[i], "-m"))
138 m = atoi(argv[++i]);
139 else if(!strcmp(argv[i], "-d"))
140 d = atoi(argv[++i]);
141 else if(!strcmp(argv[i], "-r"))
142 numReps = atoi(argv[++i]);
143 else if(!strcmp(argv[i], "-s"))
144 s = atoi(argv[++i]);
145 else if(!strcmp(argv[i], "-o"))
146 outFile = argv[++i];
147 else if(!strcmp(argv[i], "-g"))
148 deviceNum = atoi(argv[++i]);
149 else{
150 fprintf(stderr,"%s : unrecognized option.. exiting\n",argv[i]);
151 exit(1);
152 }
153 i++;
154 }
155
156 if( !n || !m || !d || !numReps || !s || !dataFile){
157 fprintf(stderr,"more arguments needed.. exiting\n");
158 exit(1);
159 }
160
161 if(numReps>n){
162 fprintf(stderr,"can't have more representatives than points.. exiting\n");
163 exit(1);
164 }
165 }
166
167
168 void readData(char *dataFile, unint rows, unint cols, real *data){
169 FILE *fp;
170 unint numRead;
171
172 fp = fopen(dataFile,"r");
173 if(fp==NULL){
174 fprintf(stderr,"error opening file.. exiting\n");
175 exit(1);
176 }
177
178 numRead = fread(data,sizeof(real),rows*cols,fp);
179 if(numRead != rows*cols){
180 fprintf(stderr,"error reading file.. exiting \n");
181 exit(1);
182 }
183 fclose(fp);
184 }
185
186
187 void readDataText(char *dataFile, unint rows, unint cols, real *data){
188 FILE *fp;
189 real t;
190
191 fp = fopen(dataFile,"r");
192 if(fp==NULL){
193 fprintf(stderr,"error opening file.. exiting\n");
194 exit(1);
195 }
196
197 for(int i=0; i<rows; i++){
198 for(int j=0; j<cols; j++){
199 if(fscanf(fp,"%f ", &t)==EOF){
200 fprintf(stderr,"error reading file.. exiting \n");
201 exit(1);
202 }
203 data[IDX( i, j, cols )]=(real)t;
204 }
205 }
206 fclose(fp);
207 }
208
209 //This function splits the data into two matrices, x and q, of
210 //their specified dimensions. The data is split randomly.
211 //It is assumed that the number of rows of data (the parameter n)
212 //is at least as large as x.r+q.r
213 void orgData(real *data, unint n, unint d, matrix x, matrix q){
214
215 unint i,fi,j;
216 unint *p;
217 p = (unint*)calloc(n,sizeof(*p));
218
219 randPerm(n,p);
220
221 for(i=0,fi=0 ; i<x.r ; i++,fi++){
222 for(j=0;j<x.c;j++){
223 x.mat[IDX(i,j,x.ld)] = data[IDX(p[fi],j,d)];
224 }
225 }
226
227 for(i=0 ; i<q.r ; i++,fi++){
228 for(j=0;j<q.c;j++){
229 q.mat[IDX(i,j,q.ld)] = data[IDX(p[fi],j,d)];
230 }
231 }
232
233 free(p);
234 }
235
236
237 //find the error rate of a set of NNs, then print it.
238 void evalNNerror(matrix x, matrix q, unint *NNs){
239 struct timeval tvB, tvE;
240 unint i;
241
242 printf("\nComputing error rates (this might take a while)\n");
243 real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
244 for(i=0;i<q.r;i++){
245 if(NNs[i]>n) printf("error");
246 ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
247 }
248
249 unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
250 gettimeofday(&tvB,NULL);
251 bruteRangeCount(x,q,ranges,cnts);
252 gettimeofday(&tvE,NULL);
253
254 long int nc=0;
255 for(i=0;i<m;i++){
256 nc += cnts[i];
257 }
258 double mean = ((double)nc)/((double)m);
259 double var = 0.0;
260 for(i=0;i<m;i++) {
261 var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
262 }
263 printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
264 printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
265
266 if(outFile){
267 FILE* fp = fopen(outFile, "a");
268 fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) );
269 fclose(fp);
270 }
271
272 free(ranges);
273 free(cnts);
274 }
275
276
277 //evals the error rate of k-nns
278 void evalKNNerror(matrix x, matrix q, intMatrix NNs){
279 struct timeval tvB, tvE;
280 unint i,j,k;
281
282 unint m = q.r;
283 printf("\nComputing error rates (this might take a while)\n");
284
285 unint *ol = (unint*)calloc( q.r, sizeof(*ol) );
286
287 intMatrix NNsB;
288 NNsB.r=q.r; NNsB.pr=q.pr; NNsB.c=NNsB.pc=32; NNsB.ld=NNsB.pc;
289 NNsB.mat = (unint*)calloc( NNsB.pr*NNsB.pc, sizeof(*NNsB.mat) );
290
291 gettimeofday(&tvB,NULL);
292 bruteK(x,q,NNsB);
293 gettimeofday(&tvE,NULL);
294
295 //calc overlap
296 for(i=0; i<m; i++){
297 for(j=0; j<K; j++){
298 for(k=0; k<K; k++){
299 ol[i] += ( NNs.mat[IDX(i, j, NNs.ld)] == NNsB.mat[IDX(i, k, NNsB.ld)] );
300 }
301 }
302 }
303
304 long int nc=0;
305 for(i=0;i<m;i++){
306 nc += ol[i];
307 }
308
309 double mean = ((double)nc)/((double)m);
310 double var = 0.0;
311 for(i=0;i<m;i++) {
312 var += (((double)ol[i])-mean)*(((double)ol[i])-mean)/((double)m);
313 }
314 printf("\tavg overlap = %6.4f/%d; std dev = %6.4f \n", mean, K, sqrt(var));
315
316 FILE* fp;
317 if(outFile){
318 fp = fopen(outFile, "a");
319 fprintf( fp, "%d %d %6.5f %6.5f ", numReps, s, mean, sqrt(var) );
320 }
321
322 real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
323 for(i=0;i<q.r;i++){
324 ranges[i] = distVec(q,x,i,NNs.mat[IDX(i, K-1, NNs.ld)]);
325 }
326
327 unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
328 bruteRangeCount(x,q,ranges,cnts);
329
330 nc=0;
331 for(i=0;i<m;i++){
332 nc += cnts[i];
333 }
334 mean = ((double)nc)/((double)m);
335 var = 0.0;
336 for(i=0;i<m;i++) {
337 var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
338 }
339 printf("\tavg actual rank of 32nd NN returned by the RBC = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
340 printf("(brute k-nn took %6.4f) \n", timeDiff(tvB, tvE));
341
342 if(outFile){
343 fprintf( fp, "%6.5f %6.5f \n", mean, sqrt(var) );
344 fclose(fp);
345 }
346
347 free(cnts);
348 free(ol);
349 free(NNsB.mat);
350 }