comparison audioDB.cpp @ 20:0519fc406b29

New major version, mostly tested: both sequence queries (-Q seq  --pointnn N and -Q seq --radius R) now work, all reported distances are Euclidean.

   case O2_FLAG_POINT_QUERY:
     pointQuery(dbName, inFile, adbQueryResult);
     break;
   case O2_FLAG_SEQUENCE_QUERY:
     if(radius==0)
-      trackSequenceQuery(dbName, inFile, adbQueryResult);
+      trackSequenceQueryNN(dbName, inFile, adbQueryResult);
     else
-      trackSequenceQueryEuc(dbName, inFile, adbQueryResult);
+      trackSequenceQueryRad(dbName, inFile, adbQueryResult);
     break;
   case O2_FLAG_TRACK_QUERY:
     trackPointQuery(dbName, inFile, adbQueryResult);
     break;
   default:

     delete qNorm;
   if(timesdata)
     delete timesdata;
   if(dbdurs)
     delete dbdurs;
-}
-
-void audioDB::sequenceQuery(const char* dbName, const char* inFile, adb__queryResult *adbQueryResult){  
-
 }
 
 // trackPointQuery  
 // return the trackNN closest tracks to the query track
 // uses average of pointNN points per track 

     delete timesdata;
   if(meanDBdur)
     delete meanDBdur;
 
 }
-  
-void audioDB::deleteDB(const char* dbName, const char* inFile){
-
-}
-
-// NBest matched filter distance between query and target tracks
-// efficient implementation
-// outputs average of N minimum matched filter distances
-void audioDB::trackSequenceQuery(const char* dbName, const char* inFile, adb__queryResult *adbQueryResult){
+
+
+// k nearest-neighbor (k-NN) search between query and target tracks
+// efficient implementation based on matched filter
+// assumes normed shingles
+// outputs distances of retrieved shingles, max retreived = pointNN shingles per per track
+void audioDB::trackSequenceQueryNN(const char* dbName, const char* inFile, adb__queryResult *adbQueryResult){
   
   initTables(dbName, inFile);
   
   // For each input vector, find the closest pointNN matching output vectors and report
   // we use stdout in this stub version

     error("Database must be L2 normed for sequence query","use -l2norm");
   
   if(verbosity>1)
     cerr << "performing norms ... "; cerr.flush();
   unsigned dbVectors = dbH->length/(sizeof(double)*dbH->dim);
+
   // 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);    
   unitNorm(queryCopy, dbH->dim, numVectors, qNorm);
   query = queryCopy;
+
   // Make norm measurements relative to sequenceLength
   unsigned w = sequenceLength-1;
   unsigned i,j;
   double* ps;
   double tmp1,tmp2;
+
   // Copy the L2 norm values to core to avoid disk random access later on
   memcpy(sNorm, l2normTable, dbVectors*sizeof(double));
   double* snPtr = sNorm;
   for(i=0; i<dbH->numFiles; i++){
-    if(trackTable[i]>sequenceLength){
+    if(trackTable[i]>=sequenceLength){
       tmp1=*snPtr;
       j=1;
       w=sequenceLength-1;
       while(w--)
 	*snPtr+=snPtr[j++];
       ps = snPtr+1;
       w=trackTable[i]-sequenceLength; // +1 - 1
       while(w--){
 	tmp2=*ps;
-	*ps=*(ps-1)-tmp1+*(ps+sequenceLength);
+	*ps=*(ps-1)-tmp1+*(ps+sequenceLength-1);
 	tmp1=tmp2;
+	ps++;
+      }
+      ps = snPtr;
+      w=trackTable[i]-sequenceLength+1;
+      while(w--){
+	*ps=sqrt(*ps);
 	ps++;
       }
     }
     snPtr+=trackTable[i];
   }

     *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;
+      w = trackTable[i]-sequenceLength;
       pn = sMeanL2+i;
+      *pn=0;
       while(w--)
-	*pn+=*ps++;
-      *pn/=trackTable[i]-sequenceLength+1;
+	if(*ps>0)
+	  *pn+=*ps++;
+      *pn/=trackTable[i]-sequenceLength;
       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/2; // 50% of the mean shingle power
-  SILENCE_THRESH/=10; // 10% of the mean shingle power is SILENCE
-  
+  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;
   w=sequenceLength-1;
   i=1;
   tmp1=*qNorm;
   while(w--)
     *qNorm+=qNorm[i++];
   ps = qNorm+1;
-  qMeanL2 = *qNorm;
-  w=numVectors-sequenceLength;
+  w=numVectors-sequenceLength; // +1 -1
   while(w--){
     tmp2=*ps;
-    *ps=*(ps-1)-tmp1+*(ps+sequenceLength);
+    *ps=*(ps-1)-tmp1+*(ps+sequenceLength-1);
     tmp1=tmp2;
-    qMeanL2+=*ps;
-    *ps++;
+    ps++;
+  }
+  ps = qNorm;
+  qMeanL2 = 0;
+  w=numVectors-sequenceLength+1;
+  while(w--){
+    *ps=sqrt(*ps);
+    qMeanL2+=*ps++;
   }
   qMeanL2 /= numVectors-sequenceLength+1;
+
   if(verbosity>1)
     cerr << "done." << endl;    
   
   
   if(verbosity>1)

   unsigned k,l,m,n,track,trackOffset=0, HOP_SIZE=sequenceHop, wL=sequenceLength;
   double thisDist;
   double oneOverWL=1.0/wL;
   
   for(k=0; k<pointNN; k++){
-    distances[k]=0.0;
+    distances[k]=1.0e6;
     qIndexes[k]=~0;
     sIndexes[k]=~0;    
   }
   
   for(k=0; k<trackNN; k++){
-    trackDistances[k]=0.0;
+    trackDistances[k]=1.0e6;
     trackQIndexes[k]=~0;
     trackSIndexes[k]=~0;
     trackIDs[k]=~0;
   }
 

       query=query+queryPoint*dbH->dim;
       qNorm=qNorm+queryPoint;
       numVectors=wL;
     }
   
-  double ** D = 0;    // Cross-correlation between query and target 
+  double ** D = 0;    // Differences query and target 
   double ** DD = 0;   // Matched filter distance
 
   D = new double*[numVectors];
   assert(D);
   DD = new double*[numVectors];

     trackOffsetTable[k]=cumTrack;
     cumTrack+=trackTable[k]*dbH->dim;
   }
 
   char nextKey [MAXSTR];
+
+  // chi^2 statistics
+  double sampleCount = 0;
+  double sampleSum = 0;
+  double logSampleSum = 0;
+  double minSample = 1e9;
+  double maxSample = 0;
+
+  // Track loop 
   for(processedTracks=0, track=0 ; processedTracks < dbH->numFiles ; track++, processedTracks++){
 
     // get trackID from file if using a control file
     if(trackFile){
       if(!trackFile->eof()){

     if(sequenceLength<trackTable[track]){  // test for short sequences
       
       if(verbosity>7)
 	cerr << track << "." << trackIndexOffset << "." << trackTable[track] << " | ";cerr.flush();
 		
-      // Cross-correlation matrix
+      // Sum products matrix
       for(j=0; j<numVectors;j++){
 	D[j]=new double[trackTable[track]]; 
 	assert(D[j]);
 
       }

       for(j=0; j<numVectors;j++){
 	DD[j]=new double[trackTable[track]];
 	assert(DD[j]);
       }
 
-      // Cross Correlation
+      double tmp;
+      // Dot product
       for(j=0; j<numVectors; j++)
 	for(k=0; k<trackTable[track]; k++){
 	  qp=query+j*dbH->dim;
 	  sp=dataBuf+trackOffset+k*dbH->dim;
 	  DD[j][k]=0.0; // Initialize matched filter array

 	    k=trackTable[track]-w;
 	    while(k--)
 	      *sp+++=*spd++;
 	  }
       }
+
       else{ // HOP_SIZE != 1
 	for(w=0; w<wL; w++)
 	  for(j=0; j<numVectors-w; j+=HOP_SIZE){
 	    sp=DD[j];
 	    spd=D[j+w]+w;

 	  cerr << "within duration tolerance." << endl;
 	  cerr.flush();
 	}
 
 	// Search for minimum distance by shingles (concatenated vectors)
-	for(j=0;j<numVectors-wL+1;j+=HOP_SIZE)
-	  for(k=0;k<trackTable[track]-wL+1;k+=HOP_SIZE){
-	    
-	    diffL2 = fabs(qNorm[j] - sNorm[k]);
+	for(j=0;j<numVectors-wL;j+=HOP_SIZE)
+	  for(k=0;k<trackTable[track]-wL;k+=HOP_SIZE){
+	    thisDist=2-(2/(qNorm[j]*sNorm[trackIndexOffset+k]))*DD[j][k];
+	    if(verbosity>10)
+	      cerr << thisDist << " " << qNorm[j] << " " << sNorm[trackIndexOffset+k] << endl;
+	    // Gather chi^2 statistics
+	    if(thisDist<minSample)
+	      minSample=thisDist;
+	    else if(thisDist>maxSample)
+	      maxSample=thisDist;
+	    if(thisDist>1e-9){
+	      sampleCount++;
+	      sampleSum+=thisDist;
+	      logSampleSum+=log(thisDist);
+	    }
+
+	    // diffL2 = fabs(qNorm[j] - sNorm[trackIndexOffset+k]);
 	    // Power test
 	    if(!USE_THRESH || 
 	       // Threshold on mean L2 of Q and S sequences
-	       (USE_THRESH && qNorm[j]>SILENCE_THRESH && sNorm[k]>SILENCE_THRESH && 
+	       (USE_THRESH && qNorm[j]>SILENCE_THRESH && sNorm[trackIndexOffset+k]>SILENCE_THRESH && 
 		// Are both query and target windows above mean energy?
-		(qNorm[j]>qMeanL2 && sNorm[k]>sMeanL2[track] &&  diffL2 < DIFF_THRESH )))
-	      thisDist=DD[j][k]*oneOverWL;
+		(qNorm[j]>qMeanL2*.25 && sNorm[trackIndexOffset+k]>sMeanL2[track]*.25))) // &&  diffL2 < DIFF_THRESH )))
+	      thisDist=thisDist; // Computed above
 	    else
-	      thisDist=0.0;
-	    
-	    // NBest match algorithm
-	    for(m=0; m<pointNN; m++){
-	      if(thisDist>=distances[m]){	  
+	      thisDist=1000000.0;
+
+	    // k-NN match algorithm
+	    m=pointNN;
+	    while(m--){
+	      if(thisDist<=distances[m])
+		if(m==0 || thisDist>=distances[m-1]){
 		// Shuffle distances up the list
 		for(l=pointNN-1; l>m; l--){
 		  distances[l]=distances[l-1];
 		  qIndexes[l]=qIndexes[l-1];
 		  sIndexes[l]=sIndexes[l-1];

 		  qIndexes[m]=queryPoint;
 		else
 		  qIndexes[m]=j;
 		sIndexes[m]=k;
 		break;
-	      }
+		}
 	    }
 	  }
 	// Calculate the mean of the N-Best matches
 	thisDist=0.0;
 	for(m=0; m<pointNN; m++)
-	  if(distances[m]<0.000001){ // Stop rubbish songs getting good scores
-	    thisDist=0.0;
-	    break;
-	  }
-	  else
-	    thisDist+=distances[m];
+	  thisDist+=distances[m];
 	thisDist/=pointNN;
 	
 	// Let's see the distances then...
 	if(verbosity>3)
 	  cerr << fileTable+track*O2_FILETABLESIZE << " " << thisDist << endl;
 
+
 	// All the track stuff goes here
 	n=trackNN;
 	while(n--){
-	  if(thisDist>=trackDistances[n]){
-	    if((n==0 || thisDist<=trackDistances[n-1])){
+	  if(thisDist<=trackDistances[n]){
+	    if((n==0 || thisDist>=trackDistances[n-1])){
 	      // Copy all values above up the queue
 	      for( l=trackNN-1 ; l > n ; l--){
 		trackDistances[l]=trackDistances[l-1];
 		trackQIndexes[l]=trackQIndexes[l-1];
 		trackSIndexes[l]=trackSIndexes[l-1];

 	  }
 	  else
 	    break;
 	}
       } // Duration match
-      
-      // per-track reset array values
-      for(unsigned k=0; k<pointNN; k++){
-	distances[k]=0.0;
-	qIndexes[k]=~0;
-	sIndexes[k]=~0;    
-      }
-      
+            
       // Clean up current track
       if(D!=NULL){
 	for(j=0; j<numVectors; j++)
 	  delete[] D[j];
       }

       if(DD!=NULL){
 	for(j=0; j<numVectors; j++)
 	  delete[] DD[j];
       }
     }
+    // per-track reset array values
+    for(unsigned k=0; k<pointNN; k++){
+      distances[k]=1.0e6;
+      qIndexes[k]=~0;
+      sIndexes[k]=~0;    
+    }
   }
 
   gettimeofday(&tv2,NULL);
-  if(verbosity>1)
+  if(verbosity>1){
     cerr << endl << "processed tracks :" << processedTracks << " matched tracks: " << successfulTracks << " elapsed time:" 
 	 << ( tv2.tv_sec*1000 + tv2.tv_usec/1000 ) - ( tv1.tv_sec*1000+tv1.tv_usec/1000 ) << " msec" << endl;
-  
+    cerr << "sampleCount: " << sampleCount << " sampleSum: " << sampleSum << " logSampleSum: " << logSampleSum 
+	 << " minSample: " << minSample << " maxSample: " << maxSample << endl;
+  }  
   if(adbQueryResult==0){
     if(verbosity>1)
       cerr<<endl;
     // Output answer
     // Loop over nearest neighbours
     for(k=0; k < min(trackNN,successfulTracks); k++)
-      cout << fileTable+trackIDs[k]*O2_FILETABLESIZE << " " << trackDistances[k] << " " << trackQIndexes[k] << " " << trackSIndexes[k] << endl;
+      cout << fileTable+trackIDs[k]*O2_FILETABLESIZE << " " << trackDistances[k] << " " 
+	   << trackQIndexes[k] << " " << trackSIndexes[k] << endl;
   }
   else{ // Process Web Services Query
     int listLen = min(trackNN, processedTracks);
     adbQueryResult->__sizeRlist=listLen;
     adbQueryResult->__sizeDist=listLen;

   }
 
 
   // Clean up
   if(trackOffsetTable)
-    delete trackOffsetTable;
+    delete[] trackOffsetTable;
   if(queryCopy)
-    delete queryCopy;
+    delete[] queryCopy;
   //if(qNorm)
   //delete qNorm;
   if(D)
     delete[] D;
   if(DD)
     delete[] DD;
   if(timesdata)
-    delete timesdata;
+    delete[] timesdata;
   if(meanDBdur)
-    delete meanDBdur;
-
-
-}
-
-// NBest matched filter distance between query and target tracks
-// efficient implementation
-// outputs average of N minimum matched filter distances
-void audioDB::trackSequenceQueryEuc(const char* dbName, const char* inFile, adb__queryResult *adbQueryResult){
+    delete[] meanDBdur;
+
+
+}
+
+// Radius search between query and target tracks
+// efficient implementation based on matched filter
+// assumes normed shingles
+// outputs count of retrieved shingles, max retreived = one shingle per query shingle per track
+void audioDB::trackSequenceQueryRad(const char* dbName, const char* inFile, adb__queryResult *adbQueryResult){
   
   initTables(dbName, inFile);
   
   // For each input vector, find the closest pointNN matching output vectors and report
   // we use stdout in this stub version