k-RBC implemented and debugged; error routines added
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Thu, 18 Nov 2010 16:38:04 +0000 (17:38 +0100)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Thu, 18 Nov 2010 16:38:04 +0000 (17:38 +0100)
driver.cu
kernelWrap.cu
kernelWrap.h
kernels.cu
rbc.cu
rbc.h

index 117e3a6..72b1495 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 evalNNerror(matrix, matrix, unint*);
+void evalKNNerror(matrix,matrix,intMatrix);
 
 char *dataFile, *outFile;
 unint n=0, m=0, d=0, numReps=0, s=0;
@@ -27,7 +28,7 @@ int main(int argc, char**argv){
   real *data;
   matrix x, q;
   unint *NNs;
-  intMatrix NNsK, NNsKCPU;
+  intMatrix NNsK, NNsKCPU, kNNsRBC;
   unint i;
   struct timeval tvB,tvE;
   cudaError_t cE;
@@ -76,92 +77,66 @@ int main(int argc, char**argv){
   
   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;
+  kNNsRBC.r=q.r; kNNsRBC.pr=q.pr; kNNsRBC.pc=kNNsRBC.c=K; kNNsRBC.ld=kNNsRBC.pc;
+  kNNsRBC.mat = (unint*)calloc(kNNsRBC.pr*kNNsRBC.pc, sizeof(*kNNsRBC.mat));
   NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
   NNsKCPU.mat = (unint*)calloc(NNsKCPU.pr*NNsKCPU.pc, sizeof(*NNsKCPU.mat));
-  
-  int tester[5][5];
-  int k=0;
-  int j;
-  for(i=0;i<5;i++){
-    for(j=0;j<5;j++){
-      tester[i][j]=k++;
-    }
-  }
 
+  /* 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 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"); */
+  /* gettimeofday(&tvB,NULL); */
+  /* bruteSearch(x,q,NNs); */
+  /* gettimeofday(&tvE,NULL); */
+  /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
 
-  printf("running regular brute force..\n");
-  gettimeofday(&tvB,NULL);
-  bruteSearch(x,q,NNs);
-  gettimeofday(&tvE,NULL);
-  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
+
+  /* printf("running cpu version..\n"); */
+  /* gettimeofday(&tvB,NULL); */
+  /* bruteKCPU(x,q,NNsKCPU); */
+  /* gettimeofday(&tvE,NULL); */
+  /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
 
 
-  printf("running cpu version..\n");
+  printf("\nrunning rbc..\n");
   gettimeofday(&tvB,NULL);
-  bruteKCPU(x,q,NNsKCPU);
+  buildRBC(x, &rbcS, numReps, s);
   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("\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") */;
+  gettimeofday(&tvB,NULL);
+  kqueryRBC(q, rbcS, kNNsRBC);
+  gettimeofday(&tvE,NULL);
+  printf("\t.. query time for krbc = %6.4f \n",timeDiff(tvB,tvE));
+  
+  destroyRBC(&rbcS);
+  printf("finished \n");
   
   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)); */
-
-
+  evalKNNerror(x,q,kNNsRBC);
+
+/* for(i=0;i<q.r;i++){ */
+    /* int j; */
+ /*    for(j=0;j<K;j++){ */
+ /*      if(NNsK.mat[IDX(i,j,NNsK.ld)] != kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)]) */
+ /*    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)])); */
+ /*    } */
+ /* /\*    printf("\n"); *\/ */
+ /*  } */
   /* if(outFile){ */
   /*   FILE* fp = fopen(outFile, "a"); */
   /*   fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) ); */
@@ -298,3 +273,102 @@ void orgData(real *data, unint n, unint d, matrix x, matrix q){
   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));
+}
+
+
+//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) );
+  
+  gettimeofday(&tvB,NULL);
+  bruteK(x,q,NNsB);
+  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));
+  
+  
+  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));
+
+  free(cnts);
+  free(ol);
+  free(NNsB.mat);
+}
index 1cee847..d46f2cb 100644 (file)
@@ -128,6 +128,24 @@ void planNNWrap(const matrix dq, const unint *dqMap, const matrix dx, const intM
 }
 
 
+void planKNNWrap(const matrix dq, const unint *dqMap, const matrix dx, const intMatrix dxMap, matrix dMins, intMatrix dMinIDs, compPlan dcP, unint compLength){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid;
+  unint todo;
+
+  grid.x = 1;
+  unint numDone = 0;
+  while( numDone<compLength ){
+    todo = MIN( (compLength-numDone) , MAX_BS*BLOCK_SIZE );
+    grid.y = todo/BLOCK_SIZE;
+    planKNNKernel<<<grid,block>>>(dq,dqMap,dx,dxMap,dMins,dMinIDs,dcP,numDone);
+    numDone += todo;
+  }
+  cudaThreadSynchronize();
+}
+
+
+
 void rangeCountWrap(const matrix dq, const matrix dx, real *dranges, unint *dcounts){
   dim3 block(BLOCK_SIZE,BLOCK_SIZE);
   dim3 grid;
index b30bd76..96f31c3 100644 (file)
@@ -14,5 +14,6 @@ void nnWrap(const matrix,const matrix,real*,unint*);
 void knnWrap(const matrix,const matrix,matrix,intMatrix);
 void rangeCountWrap(const matrix,const matrix,real*,unint*);
 void planNNWrap(const matrix,const unint*,const matrix,const intMatrix,real*,unint*,compPlan,unint);
+void planKNNWrap(const matrix,const unint*,const matrix,const intMatrix,matrix,intMatrix,compPlan,unint);
 
 #endif
index 4924638..3524b06 100644 (file)
@@ -138,7 +138,7 @@ __global__ void planKNNKernel(const matrix Q, const unint *qMap, const matrix X,
       }
      
       dNN[offQ][offX+32] = (xB+offX<groupCount)? ans:MAX_REAL;
-      idNN[offQ][offX+32] = xB + offX;
+      idNN[offQ][offX+32] = (xB+offX<groupCount)? xMap.mat[IDX( xG, xB+offX, xMap.ld )]: DUMMY_IDX; 
       __syncthreads();
 
       sort16off( dNN, idNN );
diff --git a/rbc.cu b/rbc.cu
index 51607f9..4cfef0f 100644 (file)
--- a/rbc.cu
+++ b/rbc.cu
@@ -67,6 +67,57 @@ void queryRBC(const matrix q, const rbcStruct rbcS, unint *NNs){
 }
 
 
+void kqueryRBC(const matrix q, const rbcStruct rbcS, intMatrix NNs){
+  unint m = q.r;
+  unint numReps = rbcS.dr.r;
+  unint compLength;
+  compPlan dcP;
+  unint *qMap, *dqMap;
+  qMap = (unint*)calloc(PAD(m+(BLOCK_SIZE-1)*PAD(numReps)),sizeof(*qMap));
+  matrix dq;
+  copyAndMove(&dq, &q);
+  
+  charMatrix cM;
+  cM.r=cM.c=numReps; cM.pr=cM.pc=cM.ld=PAD(numReps);
+  cM.mat = (char*)calloc( cM.pr*cM.pc, sizeof(*cM.mat) );
+  
+  unint *repIDsQ;
+  repIDsQ = (unint*)calloc( m, sizeof(*repIDsQ) );
+  real *distToRepsQ;
+  distToRepsQ = (real*)calloc( m, sizeof(*distToRepsQ) );
+  unint *groupCountQ;
+  groupCountQ = (unint*)calloc( PAD(numReps), sizeof(*groupCountQ) );
+  
+  computeReps(dq, rbcS.dr, repIDsQ, distToRepsQ);
+
+  //How many points are assigned to each group?
+  computeCounts(repIDsQ, m, groupCountQ);
+  
+  //Set up the mapping from groups to queries (qMap).
+  buildQMap(q, qMap, repIDsQ, numReps, &compLength);
+  
+  // Setup the computation matrix.  Currently, the computation matrix is 
+  // just the identity matrix: each query assigned to a particular 
+  // representative is compared only to that representative's points.  
+  idIntersection(cM);
+
+  initCompPlan(&dcP, cM, groupCountQ, rbcS.groupCount, numReps);
+
+  checkErr( cudaMalloc( (void**)&dqMap, compLength*sizeof(*dqMap) ) );
+  cudaMemcpy( dqMap, qMap, compLength*sizeof(*dqMap), cudaMemcpyHostToDevice );
+  
+  computeKNNs(rbcS.dx, rbcS.dxMap, dq, dqMap, dcP, NNs, compLength);
+
+  free(qMap);
+  freeCompPlan(&dcP);
+  cudaFree(dq.mat);
+  free(cM.mat);
+  free(repIDsQ);
+  free(distToRepsQ);
+  free(groupCountQ);
+}
+
+
 void buildRBC(const matrix x, rbcStruct *rbcS, unint numReps, unint s){
   //const matrix dr, intMatrix xmap, unint *counts, unint s){
   unint n = x.pr;
@@ -247,7 +298,6 @@ void idIntersection(charMatrix cM){
 }
 
 
-
 void fullIntersection(charMatrix cM){
   unint i,j;
   for(i=0;i<cM.r;i++){
@@ -268,12 +318,29 @@ void computeNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan dc
   planNNWrap(dq, dqMap, dx, dxMap, dMins, dMinIDs, dcP, compLength);
   cudaMemcpy( NNs, dMinIDs, dq.r*sizeof(*NNs), cudaMemcpyDeviceToHost);
   
-      
   cudaFree(dMins);
   cudaFree(dMinIDs);
 }
 
 
+void computeKNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan dcP, intMatrix NNs, unint compLength){
+  matrix dMins;
+  intMatrix dMinIDs;
+  dMins.r=compLength; dMins.pr=compLength; dMins.c=K; dMins.pc=K; dMins.ld=dMins.pc;
+  dMinIDs.r=compLength; dMinIDs.pr=compLength; dMinIDs.c=K; dMinIDs.pc=K; dMinIDs.ld=dMinIDs.pc;
+
+  checkErr( cudaMalloc((void**)&dMins.mat,dMins.pr*dMins.pc*sizeof(*dMins.mat)) );
+  checkErr( cudaMalloc((void**)&dMinIDs.mat,dMinIDs.pr*dMinIDs.pc*sizeof(*dMinIDs.mat)) );
+
+  planKNNWrap(dq, dqMap, dx, dxMap, dMins, dMinIDs, dcP, compLength);
+  cudaMemcpy( NNs.mat, dMinIDs.mat, dq.r*K*sizeof(*NNs.mat), cudaMemcpyDeviceToHost);
+
+  
+  cudaFree(dMins.mat);
+  cudaFree(dMinIDs.mat);
+}
+
+
 //This calls the dist1Kernel wrapper, but has it compute only 
 //a submatrix of the all-pairs distance matrix.  In particular,
 //only distances from dr[start,:].. dr[start+length-1] to all of x
diff --git a/rbc.h b/rbc.h
index bee5afb..eaf19a1 100644 (file)
--- a/rbc.h
+++ b/rbc.h
@@ -10,8 +10,8 @@
 
 void buildRBC(const matrix,rbcStruct*,unint, unint);
 void queryRBC(const matrix,const rbcStruct,unint*);
+void kqueryRBC(const matrix,const rbcStruct,intMatrix);
 void destroyRBC(rbcStruct*);
-
 void distSubMat(matrix,matrix,matrix,unint,unint);
 void computeReps(matrix,matrix,unint*,real*);
 void computeRadii(unint*,real*,real*,unint,unint);
@@ -22,6 +22,7 @@ void fullIntersection(charMatrix);
 void initCompPlan(compPlan*,charMatrix,unint*,unint*,unint);
 void freeCompPlan(compPlan*);
 void computeNNs(matrix,intMatrix,matrix,unint*,compPlan,unint*,unint);
+void computeKNNs(matrix,intMatrix,matrix,unint*,compPlan,intMatrix,unint);
 void setupReps(matrix,rbcStruct*,int);
 
 #endif