# HG changeset patch
# User mas01cr
# Date 1193927828 0
# Node ID b251c49379d1cc2637dbfbfae000c4a6bc63ad91
# Parent  97f704a32bf28766ce26cb5b65d08c7aa289f1db
Implement power thresholding for radius search.

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

diff -r 97f704a32bf2 -r b251c49379d1 audioDB.cpp
--- a/audioDB.cpp	Thu Nov 01 11:50:34 2007 +0000
+++ b/audioDB.cpp	Thu Nov 01 14:37:08 2007 +0000
@@ -1694,7 +1694,7 @@
   double* qnPtr = qNorm;
   double* snPtr = sNorm;
 
-  double *sPower = 0, *qPower;
+  double *sPower = 0, *qPower = 0;
   double *spPtr = 0, *qpPtr = 0;
 
   if (usingPower) {
@@ -2125,6 +2125,10 @@
     delete[] qNorm;
   if(sNorm)
     delete[] sNorm;
+  if(qPower)
+    delete[] qPower;
+  if(sPower) 
+    delete[] sPower;
   if(D)
     delete[] D;
   if(DD)
@@ -2149,13 +2153,6 @@
   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");
   
@@ -2169,8 +2166,7 @@
   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;
 
@@ -2183,9 +2179,25 @@
   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;
@@ -2193,52 +2205,50 @@
 	*ps=sqrt(*ps);
 	ps++;
       }
+      if (usingPower) {
+        ps = spPtr;
+        w = trackTable[i] - sequenceLength + 1;
+        while(w--) {
+          *ps = *ps / sequenceLength;
+          ps++;
+        }
+      }
     }
-    snPtr+=trackTable[i];
+    snPtr += trackTable[i];
+    if (usingPower) {
+      spPtr += 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++;
+  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");
     }
-    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;
+    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++;
+    }
   }
 
-  sequence_sum(qnPtr, numVectors, sequenceLength);
-
   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;    
@@ -2335,7 +2345,7 @@
   assert(DD);
 
   gettimeofday(&tv1, NULL); 
-  processedTracks=0;
+  unsigned processedTracks = 0;
   unsigned successfulTracks=0;
 
   double* qp;
@@ -2453,7 +2463,7 @@
 	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
@@ -2469,14 +2479,12 @@
 
 	    // 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
@@ -2579,8 +2587,10 @@
     delete[] qNorm;
   if(sNorm)
     delete[] sNorm;
-  if(sMeanL2)
-    delete[] sMeanL2;
+  if(qPower)
+    delete[] qPower;
+  if(sPower) 
+    delete[] sPower;
   if(D)
     delete[] D;
   if(DD)
@@ -2589,8 +2599,6 @@
     delete[] timesdata;
   if(meanDBdur)
     delete[] meanDBdur;
-
-
 }
 
 // Unit norm block of features