updated NN functions so that they return the distances as well as the indices
[RBC.git] / driver.cu
index 7bd2711..aee8ee5 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -18,7 +18,8 @@ void parseInput(int,char**);
 void readData(char*,unint,unint,real*);
 void readDataText(char*,unint,unint,real*);
 void orgData(real*,unint,unint,matrix,matrix);
 void readData(char*,unint,unint,real*);
 void readDataText(char*,unint,unint,real*);
 void orgData(real*,unint,unint,matrix,matrix);
-
+void evalNNerror(matrix, matrix, unint*);
+void evalKNNerror(matrix,matrix,intMatrix);
 
 char *dataFile, *outFile;
 unint n=0, m=0, d=0, numReps=0, s=0;
 
 char *dataFile, *outFile;
 unint n=0, m=0, d=0, numReps=0, s=0;
@@ -26,8 +27,8 @@ unint deviceNum=0;
 int main(int argc, char**argv){
   real *data;
   matrix x, q;
 int main(int argc, char**argv){
   real *data;
   matrix x, q;
-  unint *NNs, *NNsBrute;
-  unint i;
+  intMatrix nnsBrute, nnsRBC;
+  matrix distsBrute, distsRBC;
   struct timeval tvB,tvE;
   cudaError_t cE;
   rbcStruct rbcS;
   struct timeval tvB,tvE;
   cudaError_t cE;
   rbcStruct rbcS;
@@ -37,22 +38,16 @@ int main(int argc, char**argv){
   printf("*****************\n");
   
   parseInput(argc,argv);
   printf("*****************\n");
   
   parseInput(argc,argv);
-
-  cuInit(0);
+  
   printf("Using GPU #%d\n",deviceNum);
   if(cudaSetDevice(deviceNum) != cudaSuccess){
     printf("Unable to select device %d.. exiting. \n",deviceNum);
     exit(1);
   }
   
   printf("Using GPU #%d\n",deviceNum);
   if(cudaSetDevice(deviceNum) != cudaSuccess){
     printf("Unable to select device %d.. exiting. \n",deviceNum);
     exit(1);
   }
   
-  unsigned int memFree, memTot;
-  CUcontext pctx;
-  unsigned int flags=0;
-  int device;
-  cudaGetDevice(&device);
-  cuCtxCreate(&pctx,flags,device);
-  cuMemGetInfo(&memFree, &memTot);
-  printf("GPU memory free = %u/%u (MB) \n",memFree/(1024*1024),memTot/(1024*1024));
+  size_t memFree, memTot;
+  cudaMemGetInfo(&memFree, &memTot);
+  printf("GPU memory free = %lu/%lu (MB) \n",(unsigned long)memFree/(1024*1024),(unsigned long)memTot/(1024*1024));
 
   data  = (real*)calloc( (n+m)*d, sizeof(*data) );
   x.mat = (real*)calloc( PAD(n)*PAD(d), sizeof(*(x.mat)) );
 
   data  = (real*)calloc( (n+m)*d, sizeof(*data) );
   x.mat = (real*)calloc( PAD(n)*PAD(d), sizeof(*(x.mat)) );
@@ -62,35 +57,40 @@ int main(int argc, char**argv){
   x.r = n; x.c = d; x.pr = PAD(n); x.pc = PAD(d); x.ld = x.pc;
   q.r = m; q.c = d; q.pr = PAD(m); q.pc = PAD(d); q.ld = q.pc;
 
   x.r = n; x.c = d; x.pr = PAD(n); x.pc = PAD(d); x.ld = x.pc;
   q.r = m; q.c = d; q.pr = PAD(m); q.pc = PAD(d); q.ld = q.pc;
 
-  NNs = (unint*)calloc( m, sizeof(*NNs) );
-  for(i=0; i<m; i++)
-    NNs[i]=DUMMY_IDX;
-  NNsBrute = (unint*)calloc( m, sizeof(*NNsBrute) );
-
+  //Load data 
   readData(dataFile, (n+m), d, data);
   orgData(data, (n+m), d, x, q);
   free(data);
 
   readData(dataFile, (n+m), d, data);
   orgData(data, (n+m), d, x, q);
   free(data);
 
-  for(i=0;i<m;i++)
-    NNs[i]=NNsBrute[i]=DUMMY_IDX;
-  
-  /* printf("running brute force..\n"); */
-  /* gettimeofday(&tvB,NULL); */
-  /* bruteSearch(x,q,NNsBrute); */
-  /* gettimeofday(&tvE,NULL); */
-  /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
+  //Allocate space for NNs and dists
+  nnsBrute.r=q.r; nnsBrute.pr=q.pr; nnsBrute.pc=nnsBrute.c=K; nnsBrute.ld=nnsBrute.pc;
+  nnsBrute.mat = (unint*)calloc(nnsBrute.pr*nnsBrute.pc, sizeof(*nnsBrute.mat));
+  nnsRBC.r=q.r; nnsRBC.pr=q.pr; nnsRBC.pc=nnsRBC.c=K; nnsRBC.ld=nnsRBC.pc;
+  nnsRBC.mat = (unint*)calloc(nnsRBC.pr*nnsRBC.pc, sizeof(*nnsRBC.mat));
   
   
+  distsBrute.r=q.r; distsBrute.pr=q.pr; distsBrute.pc=distsBrute.c=K; distsBrute.ld=distsBrute.pc;
+  distsBrute.mat = (real*)calloc(distsBrute.pr*distsBrute.pc, sizeof(*distsBrute.mat));
+  distsRBC.r=q.r; distsRBC.pr=q.pr; distsRBC.pc=distsRBC.c=K; distsRBC.ld=distsRBC.pc;
+  distsRBC.mat = (real*)calloc(distsRBC.pr*distsRBC.pc, sizeof(*distsRBC.mat));
+
+  printf("running k-brute force..\n");
+  gettimeofday(&tvB,NULL);
+  bruteK(x,q,nnsBrute,distsBrute);
+  gettimeofday(&tvE,NULL);
+  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
+
   printf("\nrunning rbc..\n");
   gettimeofday(&tvB,NULL);
   buildRBC(x, &rbcS, numReps, s);
   gettimeofday(&tvE,NULL);
   printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
 
   printf("\nrunning rbc..\n");
   gettimeofday(&tvB,NULL);
   buildRBC(x, &rbcS, numReps, s);
   gettimeofday(&tvE,NULL);
   printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
 
+  //This finds the 32-NN; if you are only interested in the 1-NN, use queryRBC(..) instead
   gettimeofday(&tvB,NULL);
   gettimeofday(&tvB,NULL);
-  queryRBC(q, rbcS, NNs);
+  kqueryRBC(q, rbcS, nnsRBC, distsRBC);
   gettimeofday(&tvE,NULL);
   gettimeofday(&tvE,NULL);
-  printf("\t.. query time for rbc = %6.4f \n",timeDiff(tvB,tvE));
-
+  printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
+  
   destroyRBC(&rbcS);
   printf("finished \n");
   
   destroyRBC(&rbcS);
   printf("finished \n");
   
@@ -98,42 +98,15 @@ int main(int argc, char**argv){
   if( cE != cudaSuccess ){
     printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
   }
   if( cE != cudaSuccess ){
     printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
   }
-
-  printf("\nComputing error rates (this might take a while)\n");
-  real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
-  for(i=0;i<q.r;i++){
-    if(NNs[i]>n) printf("error");
-    ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
-  }
-
-  unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
-  gettimeofday(&tvB,NULL);
-  bruteRangeCount(x,q,ranges,cnts);
-  gettimeofday(&tvE,NULL);
   
   
-  long int nc=0;
-  for(i=0;i<m;i++){
-    nc += cnts[i];
-  }
-  double mean = ((double)nc)/((double)m);
-  double var = 0.0;
-  for(i=0;i<m;i++) {
-    var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
-  }
-  printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
-  printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
-
-
-  if(outFile){
-    FILE* fp = fopen(outFile, "a");
-    fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) );
-    fclose(fp);
-  }
-
-  free(ranges);
-  free(cnts);
-  free(NNs);
-  free(NNsBrute);
+  evalKNNerror(x,q,nnsRBC);
+  
+  cudaThreadExit();
+  
+  free(nnsBrute.mat);
+  free(nnsRBC.mat);
+  free(distsBrute.mat);
+  free(distsRBC.mat);
   free(x.mat);
   free(q.mat);
 }
   free(x.mat);
   free(q.mat);
 }
@@ -259,3 +232,122 @@ void orgData(real *data, unint n, unint d, matrix x, matrix q){
   free(p);
 }
 
   free(p);
 }
 
+
+//find the error rate of a set of NNs, then print it.
+void evalNNerror(matrix x, matrix q, unint *NNs){
+  struct timeval tvB, tvE;
+  unint i;
+
+  printf("\nComputing error rates (this might take a while)\n");
+  real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
+  for(i=0;i<q.r;i++){
+    if(NNs[i]>n) printf("error");
+    ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
+  }
+
+  unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
+  gettimeofday(&tvB,NULL);
+  bruteRangeCount(x,q,ranges,cnts);
+  gettimeofday(&tvE,NULL);
+  
+  long int nc=0;
+  for(i=0;i<m;i++){
+    nc += cnts[i];
+  }
+  double mean = ((double)nc)/((double)m);
+  double var = 0.0;
+  for(i=0;i<m;i++) {
+    var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
+  }
+  printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
+  printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
+  
+  if(outFile){
+    FILE* fp = fopen(outFile, "a");
+    fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) );
+    fclose(fp);
+  }
+
+  free(ranges);
+  free(cnts);
+}
+
+
+//evals the error rate of k-nns
+void evalKNNerror(matrix x, matrix q, intMatrix NNs){
+  struct timeval tvB, tvE;
+  unint i,j,k;
+
+  unint m = q.r;
+  printf("\nComputing error rates (this might take a while)\n");
+  
+  unint *ol = (unint*)calloc( q.r, sizeof(*ol) );
+  
+  intMatrix NNsB;
+  NNsB.r=q.r; NNsB.pr=q.pr; NNsB.c=NNsB.pc=32; NNsB.ld=NNsB.pc;
+  NNsB.mat = (unint*)calloc( NNsB.pr*NNsB.pc, sizeof(*NNsB.mat) );
+  matrix distsBrute;
+  distsBrute.r=q.r; distsBrute.pr=q.pr; distsBrute.c=distsBrute.pc=K; distsBrute.ld=distsBrute.pc;
+  distsBrute.mat = (real*)calloc( distsBrute.pr*distsBrute.pc, sizeof(*distsBrute.mat) );
+
+  gettimeofday(&tvB,NULL);
+  bruteK(x,q,NNsB,distsBrute);
+  gettimeofday(&tvE,NULL);
+
+   //calc overlap
+  for(i=0; i<m; i++){
+    for(j=0; j<K; j++){
+      for(k=0; k<K; k++){
+       ol[i] += ( NNs.mat[IDX(i, j, NNs.ld)] == NNsB.mat[IDX(i, k, NNsB.ld)] );
+      }
+    }
+  }
+
+  long int nc=0;
+  for(i=0;i<m;i++){
+    nc += ol[i];
+  }
+
+  double mean = ((double)nc)/((double)m);
+  double var = 0.0;
+  for(i=0;i<m;i++) {
+    var += (((double)ol[i])-mean)*(((double)ol[i])-mean)/((double)m);
+  }
+  printf("\tavg overlap = %6.4f/%d; std dev = %6.4f \n", mean, K, sqrt(var));
+
+  FILE* fp;
+  if(outFile){
+    fp = fopen(outFile, "a");
+    fprintf( fp, "%d %d %6.5f %6.5f ", numReps, s, mean, sqrt(var) );
+  }
+
+  real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
+  for(i=0;i<q.r;i++){
+    ranges[i] = distVec(q,x,i,NNs.mat[IDX(i, K-1, NNs.ld)]);
+  }
+  
+  unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
+  bruteRangeCount(x,q,ranges,cnts);
+  
+  nc=0;
+  for(i=0;i<m;i++){
+    nc += cnts[i];
+  }
+  mean = ((double)nc)/((double)m);
+  var = 0.0;
+  for(i=0;i<m;i++) {
+    var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
+  }
+  printf("\tavg actual rank of 32nd NN returned by the RBC = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
+  printf("(brute k-nn took %6.4f) \n", timeDiff(tvB, tvE));
+
+  if(outFile){
+    fprintf( fp, "%6.5f %6.5f \n", mean, sqrt(var) );
+    fclose(fp);
+  }
+
+  free(cnts);
+  free(ol);
+  free(NNsB.mat);
+  free(distsBrute.mat);
+}