cleaned up driver
[RBC.git] / driver.cu
index 6ec580f..38bc304 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #include<stdio.h>
 #include<stdlib.h>
 #include<cuda.h>
 #include "sKernel.h"
 
 void parseInput(int,char**);
-void readData(char*,int,int,real*);
-void orgData(real*,int,int,matrix,matrix);
-
+void readData(char*,matrix);
+void readDataText(char*,matrix);
+void evalNNerror(matrix, matrix, unint*);
+void evalKNNerror(matrix,matrix,intMatrix);
 
-char *dataFile, *outFile;
-int n=0, m=0, d=0, numReps=0, s=0;
-int deviceNum=0;
+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;
-  int *NNs, *NNsBrute;
-  int i;
+  intMatrix nnsRBC;
+  matrix distsRBC;
   struct timeval tvB,tvE;
   cudaError_t cE;
+  rbcStruct rbcS;
 
   printf("*****************\n");
   printf("RANDOM BALL COVER\n");
   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 = (int*)calloc( m, sizeof(*NNs) );
-  for(i=0; i<m; i++)
-    NNs[i]=-1;
-  NNsBrute = (int*)calloc( m, sizeof(*NNsBrute) );
+  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));
 
-  readData(dataFile, (n+m), d, data);
-  orgData(data, (n+m), d, x, q);
-  free(data);
+  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("db:\n"); */
-  /* printMat(x); */
-  /* printf("\nqueries: \n"); */
-  /* printMat(q); */
-  /* printf("\n\n"); */
-  
-  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)); */
-  
-  
   printf("\nrunning rbc..\n");
+  //Build the RBC
   gettimeofday(&tvB,NULL);
-  rbc(x,q,numReps,s,NNs); 
+  buildRBC(x, &rbcS, numReps, numReps);
   gettimeofday(&tvE,NULL);
-  printf("\t.. total time elapsed for rbc = %6.4f \n",timeDiff(tvB,tvE));
-  printf("finished \n");
-  
-  cE = cudaGetLastError();
-  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++)
-    ranges[i] = distL1(q,x,i,NNs[i]) - 10e-6;
+  printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
   
-
-  int *cnts = (int*)calloc(q.pr,sizeof(*cnts));
+  //This finds the 32-NNs; if you are only interested in the 1-NN, use queryRBC(..) instead
   gettimeofday(&tvB,NULL);
-  bruteRangeCount(x,q,ranges,cnts);
+  kqueryRBC(q, rbcS, nnsRBC, distsRBC);
   gettimeofday(&tvE,NULL);
+  printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
   
-  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);
+  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);
   }
-  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);
+  cE = cudaGetLastError();
+  if( cE != cudaSuccess ){
+    printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
   }
-
-  free(ranges);
-  free(cnts);
-  free(NNs);
-  free(NNsBrute);
+  
+  if( runEval )
+    evalKNNerror(x,q,nnsRBC);
+  
+  destroyRBC(&rbcS);
+  cudaThreadExit();
+  free(nnsRBC.mat);
+  free(distsRBC.mat);
   free(x.mat);
   free(q.mat);
 }
@@ -137,13 +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("\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("\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"))
@@ -152,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"))
@@ -165,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);
   }
@@ -177,9 +165,10 @@ void parseInput(int argc, char **argv){
 }
 
 
-void readData(char *dataFile, int rows, int cols, real *data){
+void readData(char *dataFile, matrix x){
+  unint i;
   FILE *fp;
-  int numRead;
+  unint numRead;
 
   fp = fopen(dataFile,"r");
   if(fp==NULL){
@@ -187,39 +176,157 @@ void readData(char *dataFile, int rows, int cols, real *data){
     exit(1);
   }
     
-  numRead = fread(data,sizeof(real),rows*cols,fp);
-  if(numRead != rows*cols){
-    fprintf(stderr,"error reading file.. exiting \n");
+  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, matrix x){
+  FILE *fp;
+  real t;
+  int i,j;
+
+  fp = fopen(dataFile,"r");
+  if(fp==NULL){
+    fprintf(stderr,"error opening file.. exiting\n");
     exit(1);
   }
+    
+  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);
+      }
+      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, int n, int d, matrix x, matrix q){
-   
-  int i,fi,j;
-  int *p;
-  p = (int*)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");
+  
+  unint *ol = (unint*)calloc( q.r, sizeof(*ol) );
   
-  randPerm(n,p);
+  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);
+}