comparison audioDB.cpp @ 155:b251c49379d1 powertable

Implement power thresholding for radius search.

(Also clean up a bit more after oneself, deleting power vectors if
appropriate)

   // Copy the L2 norm values to core to avoid disk random access later on
   memcpy(sNorm, l2normTable, dbVectors*sizeof(double));
   double* qnPtr = qNorm;
   double* snPtr = sNorm;
 
-  double *sPower = 0, *qPower;
+  double *sPower = 0, *qPower = 0;
   double *spPtr = 0, *qpPtr = 0;
 
   if (usingPower) {
     if (!(dbH->flags & O2_FLAG_POWER)) {
       error("database not power-enabled", dbName);

     delete[] queryCopy;
   if(qNorm)
     delete[] qNorm;
   if(sNorm)
     delete[] sNorm;
+  if(qPower)
+    delete[] qPower;
+  if(sPower) 
+    delete[] sPower;
   if(D)
     delete[] D;
   if(DD)
     delete[] DD;
   if(timesdata)

   // we use stdout in this stub version
   unsigned numVectors = (statbuf.st_size-sizeof(int))/(sizeof(double)*dbH->dim);
   double* query = (double*)(indata+sizeof(int));
   double* queryCopy = 0;
 
-  double qMeanL2;
-  double* sMeanL2;
-
-  unsigned USE_THRESH=0;
-  double SILENCE_THRESH=0;
-  double DIFF_THRESH=0;
-
   if(!(dbH->flags & O2_FLAG_L2NORM) )
     error("Database must be L2 normed for sequence query","use -l2norm");
   
   if(verbosity>1) {
     cerr << "performing norms ... "; cerr.flush();

   // Make a copy of the query
   queryCopy = new double[numVectors*dbH->dim];
   memcpy(queryCopy, query, numVectors*dbH->dim*sizeof(double));
   qNorm = new double[numVectors];
   sNorm = new double[dbVectors];
-  sMeanL2=new double[dbH->numFiles];
-  assert(qNorm&&sNorm&&queryCopy&&sMeanL2&&sequenceLength);    
+  assert(qNorm&&sNorm&&queryCopy&&sequenceLength);    
   unitNorm(queryCopy, dbH->dim, numVectors, qNorm);
   query = queryCopy;
 
   // Make norm measurements relative to sequenceLength
   unsigned w = sequenceLength-1;

 
   // Copy the L2 norm values to core to avoid disk random access later on
   memcpy(sNorm, l2normTable, dbVectors*sizeof(double));
   double* snPtr = sNorm;
   double* qnPtr = qNorm;
+
+  double *sPower = 0, *qPower = 0;
+  double *spPtr = 0, *qpPtr = 0;
+
+  if (usingPower) {
+    if(!(dbH->flags & O2_FLAG_POWER)) {
+      error("database not power-enabled", dbName);
+    }
+    sPower = new double[dbVectors];
+    spPtr = sPower;
+    memcpy(sPower, powerTable, dbVectors * sizeof(double));
+  }
+
   for(i=0; i<dbH->numFiles; i++){
-    if(trackTable[i]>=sequenceLength){
+    if(trackTable[i]>=sequenceLength) {
       sequence_sum(snPtr, trackTable[i], sequenceLength);
+      if (usingPower) {
+        sequence_sum(spPtr, trackTable[i], sequenceLength);
+      }
 
       ps = snPtr;
       w=trackTable[i]-sequenceLength+1;
       while(w--){
 	*ps=sqrt(*ps);
 	ps++;
       }
-    }
-    snPtr+=trackTable[i];
-  }
-  
-  double* pn = sMeanL2;
-  w=dbH->numFiles;
-  while(w--)
-    *pn++=0.0;
-  ps=sNorm;
-  unsigned processedTracks=0;
-  for(i=0; i<dbH->numFiles; i++){
-    if(trackTable[i]>sequenceLength-1){
-      w = trackTable[i]-sequenceLength+1;
-      pn = sMeanL2+i;
-      *pn=0;
-      while(w--)
-	if(*ps>0)
-	  *pn+=*ps++;
-      *pn/=trackTable[i]-sequenceLength+1;
-      SILENCE_THRESH+=*pn;
-      processedTracks++;
-    }
-    ps = sNorm + trackTable[i];
-  }
-  if(verbosity>1) {
-    cerr << "processedTracks: " << processedTracks << endl;
-  }
-    
-  SILENCE_THRESH/=processedTracks;
-  USE_THRESH=1; // Turn thresholding on
-  DIFF_THRESH=SILENCE_THRESH; //  mean shingle power
-  SILENCE_THRESH/=5; // 20% of the mean shingle power is SILENCE
-  if(verbosity>4) {
-    cerr << "silence thresh: " << SILENCE_THRESH;
-  }
-
+      if (usingPower) {
+        ps = spPtr;
+        w = trackTable[i] - sequenceLength + 1;
+        while(w--) {
+          *ps = *ps / sequenceLength;
+          ps++;
+        }
+      }
+    }
+    snPtr += trackTable[i];
+    if (usingPower) {
+      spPtr += trackTable[i];
+    }
+  }
+  
   sequence_sum(qnPtr, numVectors, sequenceLength);
+  if (usingPower) {
+    qPower = new double[numVectors];
+    qpPtr = qPower;
+    if (lseek(powerfd, sizeof(int), SEEK_SET) == (off_t) -1) {
+      error("error seeking to data", powerFileName, "lseek");
+    }
+    int count = read(powerfd, qPower, numVectors * sizeof(double));
+    if (count == -1) {
+      error("error reading data", powerFileName, "read");
+    }
+    if ((unsigned) count != numVectors * sizeof(double)) {
+      error("short read", powerFileName);
+    }
+
+    sequence_sum(qpPtr, numVectors, sequenceLength);
+    ps = qpPtr;
+    w = numVectors - sequenceLength + 1;
+    while(w--) {
+      *ps = *ps / sequenceLength;
+      ps++;
+    }
+  }
 
   ps = qnPtr;
-  qMeanL2 = 0;
-  w=numVectors-sequenceLength+1;
+  w = numVectors - sequenceLength + 1;
   while(w--){
     *ps=sqrt(*ps);
-    qMeanL2+=*ps++;
-  }
-  qMeanL2 /= numVectors-sequenceLength+1;
+  }
 
   if(verbosity>1) {
     cerr << "done." << endl;    
   }
   

   assert(D);
   DD = new double*[numVectors];
   assert(DD);
 
   gettimeofday(&tv1, NULL); 
-  processedTracks=0;
+  unsigned processedTracks = 0;
   unsigned successfulTracks=0;
 
   double* qp;
   double* sp;
   double* dp;

 
 	// Search for minimum distance by shingles (concatenated vectors)
 	for(j=0;j<=numVectors-wL;j+=HOP_SIZE)
 	  for(k=0;k<=trackTable[track]-wL;k+=HOP_SIZE){
 	    thisDist=2-(2/(qnPtr[j]*sNorm[trackIndexOffset+k]))*DD[j][k];
-	    if(verbosity>10) {
+	    if(verbosity>9) {
 	      cerr << thisDist << " " << qnPtr[j] << " " << sNorm[trackIndexOffset+k] << endl;
             }
 	    // Gather chi^2 statistics
 	    if(thisDist<minSample)
 	      minSample=thisDist;

 	      logSampleSum+=log(thisDist);
 	    }
 
 	    // diffL2 = fabs(qnPtr[j] - sNorm[trackIndexOffset+k]);
 	    // Power test
-	    if(!USE_THRESH || 
-	       // Threshold on mean L2 of Q and S sequences
-	       (USE_THRESH && qnPtr[j]>SILENCE_THRESH && sNorm[trackIndexOffset+k]>SILENCE_THRESH && 
-		// Are both query and target windows above mean energy?
-		(qnPtr[j]>qMeanL2*.25 && sNorm[trackIndexOffset+k]>sMeanL2[track]*.25))) // &&  diffL2 < DIFF_THRESH )))
-	      thisDist=thisDist; // Computed above
-	    else
-	      thisDist=1000000.0;
+            if (usingPower) {
+              if (!(powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k]))) {
+                thisDist = 1000000.0;
+              }
+            }
+
 	    if(thisDist>=0 && thisDist<=radius){
 	      distances[0]++; // increment count
 	      break; // only need one track point per query point
 	    }
 	  }

     delete[] queryCopy;
   if(qNorm)
     delete[] qNorm;
   if(sNorm)
     delete[] sNorm;
-  if(sMeanL2)
-    delete[] sMeanL2;
+  if(qPower)
+    delete[] qPower;
+  if(sPower) 
+    delete[] sPower;
   if(D)
     delete[] D;
   if(DD)
     delete[] DD;
   if(timesdata)
     delete[] timesdata;
   if(meanDBdur)
     delete[] meanDBdur;
-
-
 }
 
 // Unit norm block of features
 void audioDB::unitNorm(double* X, unsigned dim, unsigned n, double* qNorm){
   unsigned d;