cleaned up driver
[RBC.git] / driver.cu
index aa7c06c..38bc304 100644 (file)
--- a/driver.cu
+++ b/driver.cu
 #include "sKernel.h"
 
 void parseInput(int,char**);
-void readData(char*,unint,unint,real*);
-void readDataText(char*,unint,unint,real*);
-void orgData(real*,unint,unint,matrix,matrix);
-
-
-char *dataFile, *outFile;
-unint n=0, m=0, d=0, numReps=0, s=0;
+void readData(char*,matrix);
+void readDataText(char*,matrix);
+void evalNNerror(matrix, matrix, unint*);
+void evalKNNerror(matrix,matrix,intMatrix);
+
+char *dataFileX, *dataFileQ, *outFile;
+char runBrute=0, runEval=0;
+unint n=0, m=0, d=0, numReps=0;
 unint deviceNum=0;
 int main(int argc, char**argv){
-  real *data;
   matrix x, q;
-  unint *NNs;
-  intMatrix NNsK, NNsKCPU;
-  unint i;
+  intMatrix nnsRBC;
+  matrix distsRBC;
   struct timeval tvB,tvE;
   cudaError_t cE;
   rbcStruct rbcS;
@@ -38,131 +37,75 @@ int main(int argc, char**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);
   }
   
-  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));
-
-  data  = (real*)calloc( (n+m)*d, sizeof(*data) );
-  x.mat = (real*)calloc( PAD(n)*PAD(d), sizeof(*(x.mat)) );
-
-  //Need to allocate extra space, as each group of q will be padded later.
-  q.mat = (real*)calloc( PAD(m)*PAD(d), sizeof(*(q.mat)) );
-  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;
-  
-  readData(dataFile, (n+m), d, data);
-  orgData(data, (n+m), d, x, q);
-  free(data);
-
-  for(i=0;i<m;i++)
-    NNs[i]=DUMMY_IDX;
-  
-  NNsK.r=q.r; NNsK.pr=q.pr; NNsK.pc=NNsK.c=K; NNsK.ld=NNsK.pc;
-  NNsKCPU.r=q.r; NNsKCPU.pr=q.pr; NNsKCPU.pc=NNsKCPU.c=K; NNsKCPU.ld=NNsKCPU.pc;
-  NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
-  NNsKCPU.mat = (unint*)calloc(NNsKCPU.pr*NNsKCPU.pc, sizeof(*NNsKCPU.mat));
-
-  printf("running k-brute force..\n");
-  gettimeofday(&tvB,NULL);
-  bruteK(x,q,NNsK);
-  gettimeofday(&tvE,NULL);
-  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
-
-  printf("running regular brute force..\n");
+  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));
+
+  initMat( &x, n, d );
+  initMat( &q, m, d );
+  x.mat = (real*)calloc( sizeOfMat(x), sizeof(*(x.mat)) );
+  q.mat = (real*)calloc( sizeOfMat(q), sizeof(*(q.mat)) );
+    
+  //Load data 
+  readData( dataFileX, x );
+  readData( dataFileQ, q );
+
+  //Allocate space for NNs and dists
+  initIntMat( &nnsRBC, m, K );
+  initMat( &distsRBC, m, K );
+  nnsRBC.mat = (unint*)calloc( sizeOfIntMat(nnsRBC), sizeof(*nnsRBC.mat) );
+  distsRBC.mat = (real*)calloc( sizeOfMat(distsRBC), sizeof(*distsRBC.mat) );
+
+  printf("\nrunning rbc..\n");
+  //Build the RBC
   gettimeofday(&tvB,NULL);
-  bruteSearch(x,q,NNs);
+  buildRBC(x, &rbcS, numReps, numReps);
   gettimeofday(&tvE,NULL);
-  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
-
-
-  printf("running cpu version..\n");
+  printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
+  
+  //This finds the 32-NNs; if you are only interested in the 1-NN, use queryRBC(..) instead
   gettimeofday(&tvB,NULL);
-  bruteKCPU(x,q,NNsKCPU);
+  kqueryRBC(q, rbcS, nnsRBC, distsRBC);
   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)); */
-
-  /* gettimeofday(&tvB,NULL); */
-  /* queryRBC(q, rbcS, NNs); */
-  /* gettimeofday(&tvE,NULL); */
-  /* printf("\t.. query time for rbc = %6.4f \n",timeDiff(tvB,tvE)); */
-
-  /* destroyRBC(&rbcS); */
-  /* printf("finished \n") */;
+  printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
   
+  if( runBrute ){
+    intMatrix nnsBrute;
+    matrix distsBrute;
+    initIntMat( &nnsBrute, m, K );
+    nnsBrute.mat = (unint*)calloc( sizeOfIntMat(nnsBrute), sizeof(*nnsBrute.mat) );
+    initMat( &distsBrute, m, K );
+    distsBrute.mat = (real*)calloc( sizeOfMat(distsBrute), sizeof(*distsBrute.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));
+    
+    free(nnsBrute.mat);
+    free(distsBrute.mat);
+  }
+
   cE = cudaGetLastError();
   if( cE != cudaSuccess ){
     printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
   }
-
-  for(i=0;i<q.r;i++){
-    int j;
-    for(j=0;j<K;j++){
-      if(NNsK.mat[IDX(i,j,NNsK.ld)] != NNsKCPU.mat[IDX(i,j,NNsKCPU.ld)])
-       printf("%d %d: (%6.2f %6.2f) ",i,j,distVec(q,x,i,NNsK.mat[IDX(i,j,NNsK.ld)]),distVec(q,x,i,NNsKCPU.mat[IDX(i,j,NNsKCPU.ld)]));
-    }
- /*    printf("\n"); */
-  }
-  /* 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(NNsK.mat);
-  free(NNsKCPU.mat);
+  if( runEval )
+    evalKNNerror(x,q,nnsRBC);
+  
+  destroyRBC(&rbcS);
+  cudaThreadExit();
+  free(nnsRBC.mat);
+  free(distsRBC.mat);
   free(x.mat);
   free(q.mat);
 }
@@ -171,22 +114,26 @@ int main(int argc, char**argv){
 void parseInput(int argc, char **argv){
   int i=1;
   if(argc <= 1){
-    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");
-    printf("\tdatafile     = binary file containing the data\n");
+    printf("\nusage: \n  testRBC -x datafileX -q datafileQ  -n numPts (DB) -m numQueries -d dim -r numReps [-o outFile] [-g GPU num] [-b] [-e]\n\n");
+    printf("\tdatafileX    = binary file containing the database\n");
+    printf("\tdatafileQ    = binary file containing the queries\n");
     printf("\tnumPts       = size of database\n");
     printf("\tnumQueries   = number of queries\n");
     printf("\tdim          = dimensionailty\n");
     printf("\tnumReps      = number of representatives\n");
-    printf("\tnumPtsPerRep = number of points assigned to each representative\n");
     printf("\toutFile      = output file (optional); stored in text format\n");
     printf("\tGPU num      = ID # of the GPU to use (optional) for multi-GPU machines\n");
+    printf("\n\tuse -b to run brute force in addition the RBC\n");
+    printf("\tuse -e option to run evaluation routine\n");
     printf("\n\n");
     exit(0);
   }
   
   while(i<argc){
-    if(!strcmp(argv[i], "-f"))
-      dataFile = argv[++i];
+    if(!strcmp(argv[i], "-x"))
+      dataFileX = argv[++i];
+    else if(!strcmp(argv[i], "-q"))
+      dataFileQ = argv[++i];
     else if(!strcmp(argv[i], "-n"))
       n = atoi(argv[++i]);
     else if(!strcmp(argv[i], "-m"))
@@ -195,8 +142,6 @@ void parseInput(int argc, char **argv){
       d = atoi(argv[++i]);
     else if(!strcmp(argv[i], "-r"))
       numReps = atoi(argv[++i]);
-    else if(!strcmp(argv[i], "-s"))
-      s = atoi(argv[++i]);
     else if(!strcmp(argv[i], "-o"))
       outFile = argv[++i];
     else if(!strcmp(argv[i], "-g"))
@@ -208,7 +153,7 @@ void parseInput(int argc, char **argv){
     i++;
   }
 
-  if( !n || !m || !d || !numReps || !s || !dataFile){
+  if( !n || !m || !d || !numReps || !dataFileX || !dataFileQ ){
     fprintf(stderr,"more arguments needed.. exiting\n");
     exit(1);
   }
@@ -220,7 +165,8 @@ void parseInput(int argc, char **argv){
 }
 
 
-void readData(char *dataFile, unint rows, unint cols, real *data){
+void readData(char *dataFile, matrix x){
+  unint i;
   FILE *fp;
   unint numRead;
 
@@ -230,18 +176,22 @@ void readData(char *dataFile, unint rows, unint cols, real *data){
     exit(1);
   }
     
-  numRead = fread(data,sizeof(real),rows*cols,fp);
-  if(numRead != rows*cols){
-    fprintf(stderr,"error reading file.. exiting \n");
-    exit(1);
+  for( i=0; i<x.r; i++ ){ //can't load everything in one fread
+                           //because matrix is padded.
+    numRead = fread( &x.mat[IDX( i, 0, x.ld )], sizeof(real), x.c, fp );
+    if(numRead != x.c){
+      fprintf(stderr,"error reading file.. exiting \n");
+      exit(1);
+    }
   }
   fclose(fp);
 }
 
 
-void readDataText(char *dataFile, unint rows, unint cols, real *data){
+void readDataText(char *dataFile, matrix x){
   FILE *fp;
   real t;
+  int i,j;
 
   fp = fopen(dataFile,"r");
   if(fp==NULL){
@@ -249,42 +199,134 @@ void readDataText(char *dataFile, unint rows, unint cols, real *data){
     exit(1);
   }
     
-  for(int i=0; i<rows; i++){
-    for(int j=0; j<cols; j++){
+  for(i=0; i<x.r; i++){
+    for(j=0; j<x.c; j++){
       if(fscanf(fp,"%f ", &t)==EOF){
        fprintf(stderr,"error reading file.. exiting \n");
        exit(1);
       }
-      data[IDX( i, j, cols )]=(real)t;
+      x.mat[IDX( i, j, x.ld )]=(real)t;
     }
   }
   fclose(fp);
 }
 
-//This function splits the data into two matrices, x and q, of 
-//their specified dimensions.  The data is split randomly.
-//It is assumed that the number of rows of data (the parameter n)
-//is at least as large as x.r+q.r
-void orgData(real *data, unint n, unint d, matrix x, matrix q){
-   
-  unint i,fi,j;
-  unint *p;
-  p = (unint*)calloc(n,sizeof(*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 %6.5f %6.5f \n", numReps, 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");
   
-  randPerm(n,p);
+  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) );
 
-  for(i=0,fi=0 ; i<x.r ; i++,fi++){
-    for(j=0;j<x.c;j++){
-      x.mat[IDX(i,j,x.ld)] = data[IDX(p[fi],j,d)];
+  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)] );
+      }
     }
   }
 
-  for(i=0 ; i<q.r ; i++,fi++){
-    for(j=0;j<q.c;j++){
-      q.mat[IDX(i,j,q.ld)] = data[IDX(p[fi],j,d)];
-    } 
+  long int nc=0;
+  for(i=0;i<m;i++){
+    nc += ol[i];
   }
 
-  free(p);
-}
+  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 %6.5f %6.5f ", numReps, 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);
+}