k-RBC implemented and debugged; error routines added
[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, NNsKCPU, 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 for(i=0;i<m;i++)
76 NNs[i]=DUMMY_IDX;
77
78 NNsK.r=q.r; NNsK.pr=q.pr; NNsK.pc=NNsK.c=K; NNsK.ld=NNsK.pc;
79 NNsKCPU.r=q.r; NNsKCPU.pr=q.pr; NNsKCPU.pc=NNsKCPU.c=K; NNsKCPU.ld=NNsKCPU.pc;
80 kNNsRBC.r=q.r; kNNsRBC.pr=q.pr; kNNsRBC.pc=kNNsRBC.c=K; kNNsRBC.ld=kNNsRBC.pc;
81 kNNsRBC.mat = (unint*)calloc(kNNsRBC.pr*kNNsRBC.pc, sizeof(*kNNsRBC.mat));
82 NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
83 NNsKCPU.mat = (unint*)calloc(NNsKCPU.pr*NNsKCPU.pc, sizeof(*NNsKCPU.mat));
84
85 /* printf("running k-brute force..\n"); */
86 /* gettimeofday(&tvB,NULL); */
87 /* bruteK(x,q,NNsK); */
88 /* gettimeofday(&tvE,NULL); */
89 /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
90
91 /* printf("running regular brute force..\n"); */
92 /* gettimeofday(&tvB,NULL); */
93 /* bruteSearch(x,q,NNs); */
94 /* gettimeofday(&tvE,NULL); */
95 /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
96
97
98 /* printf("running cpu version..\n"); */
99 /* gettimeofday(&tvB,NULL); */
100 /* bruteKCPU(x,q,NNsKCPU); */
101 /* gettimeofday(&tvE,NULL); */
102 /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
103
104
105 printf("\nrunning rbc..\n");
106 gettimeofday(&tvB,NULL);
107 buildRBC(x, &rbcS, numReps, s);
108 gettimeofday(&tvE,NULL);
109 printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
110
111 /* gettimeofday(&tvB,NULL); */
112 /* queryRBC(q, rbcS, NNs); */
113 /* gettimeofday(&tvE,NULL); */
114 /* printf("\t.. query time for rbc = %6.4f \n",timeDiff(tvB,tvE)); */
115
116 gettimeofday(&tvB,NULL);
117 kqueryRBC(q, rbcS, kNNsRBC);
118 gettimeofday(&tvE,NULL);
119 printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
120
121 destroyRBC(&rbcS);
122 printf("finished \n");
123
124 cE = cudaGetLastError();
125 if( cE != cudaSuccess ){
126 printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
127 }
128
129 evalKNNerror(x,q,kNNsRBC);
130
131 /* for(i=0;i<q.r;i++){ */
132 /* int j; */
133 /* for(j=0;j<K;j++){ */
134 /* if(NNsK.mat[IDX(i,j,NNsK.ld)] != kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)]) */
135 /* printf("%d %d: (%d: %6.2f %d: %6.2f) ",i,j,NNsK.mat[IDX(i,j,NNsK.ld)],distVec(q,x,i,NNsK.mat[IDX(i,j,NNsK.ld)]),kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)] ,distVec(q,x,i,kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)])); */
136 /* } */
137 /* /\* printf("\n"); *\/ */
138 /* } */
139
140 /* if(outFile){ */
141 /* FILE* fp = fopen(outFile, "a"); */
142 /* fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) ); */
143 /* fclose(fp); */
144 /* } */
145
146 /* free(ranges); */
147 /* free(cnts); */
148 free(NNs);
149 free(NNsK.mat);
150 free(NNsKCPU.mat);
151 free(x.mat);
152 free(q.mat);
153 }
154
155
156 void parseInput(int argc, char **argv){
157 int i=1;
158 if(argc <= 1){
159 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");
160 printf("\tdatafile = binary file containing the data\n");
161 printf("\tnumPts = size of database\n");
162 printf("\tnumQueries = number of queries\n");
163 printf("\tdim = dimensionailty\n");
164 printf("\tnumReps = number of representatives\n");
165 printf("\tnumPtsPerRep = number of points assigned to each representative\n");
166 printf("\toutFile = output file (optional); stored in text format\n");
167 printf("\tGPU num = ID # of the GPU to use (optional) for multi-GPU machines\n");
168 printf("\n\n");
169 exit(0);
170 }
171
172 while(i<argc){
173 if(!strcmp(argv[i], "-f"))
174 dataFile = argv[++i];
175 else if(!strcmp(argv[i], "-n"))
176 n = atoi(argv[++i]);
177 else if(!strcmp(argv[i], "-m"))
178 m = atoi(argv[++i]);
179 else if(!strcmp(argv[i], "-d"))
180 d = atoi(argv[++i]);
181 else if(!strcmp(argv[i], "-r"))
182 numReps = atoi(argv[++i]);
183 else if(!strcmp(argv[i], "-s"))
184 s = atoi(argv[++i]);
185 else if(!strcmp(argv[i], "-o"))
186 outFile = argv[++i];
187 else if(!strcmp(argv[i], "-g"))
188 deviceNum = atoi(argv[++i]);
189 else{
190 fprintf(stderr,"%s : unrecognized option.. exiting\n",argv[i]);
191 exit(1);
192 }
193 i++;
194 }
195
196 if( !n || !m || !d || !numReps || !s || !dataFile){
197 fprintf(stderr,"more arguments needed.. exiting\n");
198 exit(1);
199 }
200
201 if(numReps>n){
202 fprintf(stderr,"can't have more representatives than points.. exiting\n");
203 exit(1);
204 }
205 }
206
207
208 void readData(char *dataFile, unint rows, unint cols, real *data){
209 FILE *fp;
210 unint numRead;
211
212 fp = fopen(dataFile,"r");
213 if(fp==NULL){
214 fprintf(stderr,"error opening file.. exiting\n");
215 exit(1);
216 }
217
218 numRead = fread(data,sizeof(real),rows*cols,fp);
219 if(numRead != rows*cols){
220 fprintf(stderr,"error reading file.. exiting \n");
221 exit(1);
222 }
223 fclose(fp);
224 }
225
226
227 void readDataText(char *dataFile, unint rows, unint cols, real *data){
228 FILE *fp;
229 real t;
230
231 fp = fopen(dataFile,"r");
232 if(fp==NULL){
233 fprintf(stderr,"error opening file.. exiting\n");
234 exit(1);
235 }
236
237 for(int i=0; i<rows; i++){
238 for(int j=0; j<cols; j++){
239 if(fscanf(fp,"%f ", &t)==EOF){
240 fprintf(stderr,"error reading file.. exiting \n");
241 exit(1);
242 }
243 data[IDX( i, j, cols )]=(real)t;
244 }
245 }
246 fclose(fp);
247 }
248
249 //This function splits the data into two matrices, x and q, of
250 //their specified dimensions. The data is split randomly.
251 //It is assumed that the number of rows of data (the parameter n)
252 //is at least as large as x.r+q.r
253 void orgData(real *data, unint n, unint d, matrix x, matrix q){
254
255 unint i,fi,j;
256 unint *p;
257 p = (unint*)calloc(n,sizeof(*p));
258
259 randPerm(n,p);
260
261 for(i=0,fi=0 ; i<x.r ; i++,fi++){
262 for(j=0;j<x.c;j++){
263 x.mat[IDX(i,j,x.ld)] = data[IDX(p[fi],j,d)];
264 }
265 }
266
267 for(i=0 ; i<q.r ; i++,fi++){
268 for(j=0;j<q.c;j++){
269 q.mat[IDX(i,j,q.ld)] = data[IDX(p[fi],j,d)];
270 }
271 }
272
273 free(p);
274 }
275
276
277 //find the error rate of a set of NNs, then print it.
278 void evalNNerror(matrix x, matrix q, unint *NNs){
279 struct timeval tvB, tvE;
280 unint i;
281
282 printf("\nComputing error rates (this might take a while)\n");
283 real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
284 for(i=0;i<q.r;i++){
285 if(NNs[i]>n) printf("error");
286 ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
287 }
288
289 unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
290 gettimeofday(&tvB,NULL);
291 bruteRangeCount(x,q,ranges,cnts);
292 gettimeofday(&tvE,NULL);
293
294 long int nc=0;
295 for(i=0;i<m;i++){
296 nc += cnts[i];
297 }
298 double mean = ((double)nc)/((double)m);
299 double var = 0.0;
300 for(i=0;i<m;i++) {
301 var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
302 }
303 printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
304 printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
305 }
306
307
308 //evals the error rate of k-nns
309 void evalKNNerror(matrix x, matrix q, intMatrix NNs){
310 struct timeval tvB, tvE;
311 unint i,j,k;
312
313 unint m = q.r;
314 printf("\nComputing error rates (this might take a while)\n");
315
316 unint *ol = (unint*)calloc( q.r, sizeof(*ol) );
317
318
319 intMatrix NNsB;
320 NNsB.r=q.r; NNsB.pr=q.pr; NNsB.c=NNsB.pc=32; NNsB.ld=NNsB.pc;
321 NNsB.mat = (unint*)calloc( NNsB.pr*NNsB.pc, sizeof(*NNsB.mat) );
322
323 gettimeofday(&tvB,NULL);
324 bruteK(x,q,NNsB);
325 gettimeofday(&tvE,NULL);
326
327 //calc overlap
328 for(i=0; i<m; i++){
329 for(j=0; j<K; j++){
330 for(k=0; k<K; k++){
331 ol[i] += ( NNs.mat[IDX(i, j, NNs.ld)] == NNsB.mat[IDX(i, k, NNsB.ld)] );
332 }
333 }
334 }
335
336 long int nc=0;
337 for(i=0;i<m;i++){
338 nc += ol[i];
339 }
340
341 double mean = ((double)nc)/((double)m);
342 double var = 0.0;
343 for(i=0;i<m;i++) {
344 var += (((double)ol[i])-mean)*(((double)ol[i])-mean)/((double)m);
345 }
346 printf("\tavg overlap = %6.4f/%d; std dev = %6.4f \n", mean, K, sqrt(var));
347
348
349 real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
350 for(i=0;i<q.r;i++){
351 ranges[i] = distVec(q,x,i,NNs.mat[IDX(i, K-1, NNs.ld)]);
352 }
353
354 unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
355 bruteRangeCount(x,q,ranges,cnts);
356
357 nc=0;
358 for(i=0;i<m;i++){
359 nc += cnts[i];
360 }
361 mean = ((double)nc)/((double)m);
362 var = 0.0;
363 for(i=0;i<m;i++) {
364 var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
365 }
366 printf("\tavg actual rank of 32nd NN returned by the RBC = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
367
368
369 printf("(brute k-nn took %6.4f) \n", timeDiff(tvB, tvE));
370
371 free(cnts);
372 free(ol);
373 free(NNsB.mat);
374 }