changeset 227:083956accbce refactoring

And now, the big one: two into one will go.

* Reporter abstract base class;
* Two instantiations of same, replicating the point accumulating and
  report functionality of trackSequenceQueryNN() and
  trackSequenceQueryRad();
* adjust trackSequenceQueryNN() to
  - include radius threshold test if radius is non-zero;
  - use new Reporter argument for accumulating matches;
  - move report() call into query();
* new VERB_LOG macro to take away zillions of std::cerr << things;
* we can now delete trackSequenceQueryRad() entirely.

WOOHOO!

--- a/audioDB.h	Wed Dec 05 13:51:34 2007 +0000
+++ b/audioDB.h	Thu Dec 06 14:41:07 2007 +0000
@@ -114,6 +114,7 @@
   off_t dbSize;
 } dbTableHeaderT, *dbTableHeaderPtr;
 
+class Reporter;
 
 class audioDB{
   
@@ -196,8 +197,7 @@
   void read_data(int track, double **data_buffer_p, size_t *data_buffer_size_p);
   void set_up_query(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned int *nvp);
   void set_up_db(double **snp, double **vsnp, double **spp, double **vspp, double **mddp, unsigned int *dvp);
-  void trackSequenceQueryNN(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse=0);
-  void trackSequenceQueryRad(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse=0);
+  void trackSequenceQueryNN(const char* dbName, const char* inFile, Reporter *reporter);
 
   void initDBHeader(const char *dbName);
   void initInputFile(const char *inFile);

--- a/query.cpp	Wed Dec 05 13:51:34 2007 +0000
+++ b/query.cpp	Thu Dec 06 14:41:07 2007 +0000
@@ -1,5 +1,218 @@
 #include "audioDB.h"
 
+#include <utility>
+#include <queue>
+#include <set>
+#include <functional>
+
+#define VERB_LOG(vv, ...) \
+  if(verbosity > vv) { \
+    fprintf(stderr, __VA_ARGS__); \
+    fflush(stderr); \
+  }
+
+typedef struct nnresult {
+  unsigned int trackID;
+  double dist;
+  unsigned int qpos;
+  unsigned int spos;
+} NNresult;
+
+typedef struct radresult {
+  unsigned int trackID;
+  unsigned int count;
+} Radresult;
+
+bool operator< (const NNresult &a, const NNresult &b) {
+  return a.dist < b.dist;
+}
+
+bool operator<= (const NNresult &a, const NNresult &b) {
+  return a.dist <= b.dist;
+}
+
+bool operator< (const Radresult &a, const Radresult &b) {
+  return a.count < b.count;
+}
+
+class Reporter {
+public:
+  virtual ~Reporter() {};
+  virtual void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) = 0;
+  virtual void report(char *fileTable, adb__queryResponse *adbQueryResponse) = 0;
+};
+
+class trackSequenceQueryNNReporter : public Reporter {
+public:
+  trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles);
+  ~trackSequenceQueryNNReporter();
+  void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
+  void report(char *fileTable, adb__queryResponse *adbQueryResponse);
+private:
+  unsigned int pointNN;
+  unsigned int trackNN;
+  unsigned int numFiles;
+  std::priority_queue< NNresult > *queues;
+};
+
+trackSequenceQueryNNReporter::trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles)
+  : pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) {
+  queues = new std::priority_queue< NNresult >[numFiles];
+}
+
+trackSequenceQueryNNReporter::~trackSequenceQueryNNReporter() {
+  delete [] queues;
+}
+
+void trackSequenceQueryNNReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
+  NNresult r;
+  r.trackID = trackID;
+  r.qpos = qpos;
+  r.spos = spos;
+  r.dist = dist;
+  queues[trackID].push(r);
+  if(queues[trackID].size() > pointNN) {
+    queues[trackID].pop();
+  }
+}
+
+void trackSequenceQueryNNReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
+  std::priority_queue < NNresult > result;
+  for (int i = numFiles-1; i >= 0; i--) {
+    unsigned int size = queues[i].size();
+    if (size > 0) {
+      NNresult r;
+      double dist = 0;
+      NNresult oldr = queues[i].top();
+      for (unsigned int j = 0; j < size; j++) {
+        r = queues[i].top();
+        dist += r.dist;
+        queues[i].pop();
+        if (r.dist == oldr.dist) {
+          r.qpos = oldr.qpos;
+          r.spos = oldr.spos;
+        } else {
+          oldr = r;
+        }
+      }
+      dist /= size;
+      r.dist = dist; // trackID, qpos and spos are magically right already.
+      result.push(r);
+      if (result.size() > trackNN) {
+        result.pop();
+      }
+    }
+  }
+
+  NNresult r;
+  std::vector<NNresult> v;
+  unsigned int size = result.size();
+  for(unsigned int k = 0; k < size; k++) {
+    r = result.top();
+    v.push_back(r);
+    result.pop();
+  }
+  std::vector<NNresult>::reverse_iterator rit;
+      
+  if(adbQueryResponse==0) {
+    for(rit = v.rbegin(); rit < v.rend(); rit++) {
+      r = *rit;
+      std::cout << fileTable + r.trackID*O2_FILETABLESIZE << " ";
+      std::cout << r.dist << " " << r.qpos << " " << r.spos << std::endl;
+    }
+  } else {
+    adbQueryResponse->result.__sizeRlist=size;
+    adbQueryResponse->result.__sizeDist=size;
+    adbQueryResponse->result.__sizeQpos=size;
+    adbQueryResponse->result.__sizeSpos=size;
+    adbQueryResponse->result.Rlist= new char*[size];
+    adbQueryResponse->result.Dist = new double[size];
+    adbQueryResponse->result.Qpos = new unsigned int[size];
+    adbQueryResponse->result.Spos = new unsigned int[size];
+    unsigned int k = 0;
+    for(rit = v.rbegin(); rit < v.rend(); rit++, k++) {
+      r = *rit;
+      adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR];
+      adbQueryResponse->result.Dist[k] = r.dist;
+      adbQueryResponse->result.Qpos[k] = r.qpos;
+      adbQueryResponse->result.Spos[k] = r.spos;
+      snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLESIZE);
+    }
+  }
+}
+
+class trackSequenceQueryRadReporter : public Reporter { 
+public:
+  trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles);
+  ~trackSequenceQueryRadReporter();
+  void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
+  void report(char *fileTable, adb__queryResponse *adbQueryResponse);
+private:
+  unsigned int trackNN;
+  unsigned int numFiles;
+  std::set<std::pair<unsigned int, unsigned int> > *set;
+  unsigned int *count;
+};
+
+trackSequenceQueryRadReporter::trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles):
+  trackNN(trackNN), numFiles(numFiles) {
+  set = new std::set<std::pair<unsigned int, unsigned int> >;
+  count = new unsigned int[numFiles];
+  for (unsigned i = 0; i < numFiles; i++) {
+    count[i] = 0;
+  }
+}
+
+trackSequenceQueryRadReporter::~trackSequenceQueryRadReporter() {
+  delete set;
+  delete [] count;
+}
+
+void trackSequenceQueryRadReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
+  std::set<std::pair<unsigned int, unsigned int> >::iterator it;
+  std::pair<unsigned int, unsigned int> pair = std::make_pair(trackID, qpos);
+  it = set->find(pair);
+  if (it == set->end()) {
+    set->insert(pair);
+    count[trackID]++;
+  }
+}
+
+void trackSequenceQueryRadReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
+  std::priority_queue < Radresult > result;
+  // KLUDGE: doing this backwards in an attempt to get the same
+  // tiebreak behaviour as before.
+  for (int i = numFiles-1; i >= 0; i--) {
+    Radresult r;
+    r.trackID = i;
+    r.count = count[i];
+    if(r.count > 0) {
+      result.push(r);
+      if (result.size() > trackNN) {
+        result.pop();
+      }
+    }
+  }
+
+  Radresult r;
+  std::vector<Radresult> v;
+  unsigned int size = result.size();
+  for(unsigned int k = 0; k < size; k++) {
+    r = result.top();
+    v.push_back(r);
+    result.pop();
+  }
+  std::vector<Radresult>::reverse_iterator rit;
+      
+  if(adbQueryResponse==0) {
+    for(rit = v.rbegin(); rit < v.rend(); rit++) {
+      r = *rit;
+      std::cout << fileTable + r.trackID*O2_FILETABLESIZE << " " << r.count << std::endl;
+    }
+  } else {
+  }
+}
+
 bool audioDB::powers_acceptable(double p1, double p2) {
   if (use_absolute_threshold) {
     if ((p1 < absolute_threshold) || (p2 < absolute_threshold)) {
@@ -15,12 +228,19 @@
 }
 
 void audioDB::query(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse) {
-  switch(queryType) {
+  initTables(dbName, inFile);
+
+  switch (queryType) {
   case O2_SEQUENCE_QUERY:
-    if(radius==0)
-      trackSequenceQueryNN(dbName, inFile, adbQueryResponse);
-    else
-      trackSequenceQueryRad(dbName, inFile, adbQueryResponse);
+    Reporter *r;
+    if(radius == 0) {
+      r = new trackSequenceQueryNNReporter(pointNN, trackNN, dbH->numFiles);
+    } else {
+      r = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles);
+    }
+    trackSequenceQueryNN(dbName, inFile, r);
+    r->report(fileTable, adbQueryResponse);
+    delete r;
     break;
   default:
     error("unrecognized queryType in query()");
@@ -182,9 +402,7 @@
     error("Query shorter than requested sequence length", "maybe use -l");
   }
   
-  if(verbosity>1) {
-    std::cerr << "performing norms ... "; std::cerr.flush();
-  }
+  VERB_LOG(1, "performing norms... ");
 
   *qp = new double[*nvp * dbH->dim];
   memcpy(*qp, indata+sizeof(int), *nvp * dbH->dim * sizeof(double));
@@ -222,9 +440,9 @@
       *mqdp += querydurs[k];
     }
     *mqdp /= k;
-    if(verbosity>1) {
-      std::cerr << "mean query file duration: " << *mqdp << std::endl;
-    }
+
+    VERB_LOG(1, "mean query file duration: %f\n", *mqdp);
+
     delete [] querydurs;
     delete [] timesdata;
   }
@@ -238,9 +456,7 @@
     if(queryPoint > *nvp || queryPoint > *nvp - sequenceLength + 1) {
       error("queryPoint > numVectors-wL+1 in query");
     } else {
-      if(verbosity>1) {
-	std::cerr << "query point: " << queryPoint << std::endl; std::cerr.flush();
-      }
+      VERB_LOG(1, "query point: %ud\n", queryPoint);
       *vqp = *qp + queryPoint * dbH->dim;
       *vqnp = *qnp + queryPoint;
       if (usingPower) {
@@ -307,9 +523,7 @@
   *vspp = *spp;
 }
 
-void audioDB::trackSequenceQueryNN(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse){
-  
-  initTables(dbName, inFile);
+void audioDB::trackSequenceQueryNN(const char* dbName, const char* inFile, Reporter *reporter) {
   
   unsigned int numVectors;
   double *query, *query_data;
@@ -323,51 +537,21 @@
   double *meanDBdur = 0;
 
   set_up_db(&sNorm, &snPtr, &sPower, &spPtr, &meanDBdur, &dbVectors);
-  
-  if(verbosity>1) {
-    std::cerr << "matching tracks..." << std::endl;
-  }
+
+  VERB_LOG(1, "matching tracks...");
   
   assert(pointNN>0 && pointNN<=O2_MAXNN);
   assert(trackNN>0 && trackNN<=O2_MAXNN);
-  
-  // Make temporary dynamic memory for results
-  double trackDistances[trackNN];
-  unsigned trackIDs[trackNN];
-  unsigned trackQIndexes[trackNN];
-  unsigned trackSIndexes[trackNN];
-  
-  double distances[pointNN];
-  unsigned qIndexes[pointNN];
-  unsigned sIndexes[pointNN];
-  
-  unsigned j,k,l,m,n,track,trackOffset=0, HOP_SIZE=sequenceHop, wL=sequenceLength;
-  double thisDist;
-  
-  for(k=0; k<pointNN; k++){
-    distances[k]=1.0e6;
-    qIndexes[k]=~0;
-    sIndexes[k]=~0;    
-  }
-  
-  for(k=0; k<trackNN; k++){
-    trackDistances[k]=1.0e6;
-    trackQIndexes[k]=~0;
-    trackSIndexes[k]=~0;
-    trackIDs[k]=~0;
-  }
 
-  double ** D = 0;    // Differences query and target 
-  double ** DD = 0;   // Matched filter distance
+  unsigned j,k,track,trackOffset=0, HOP_SIZE=sequenceHop, wL=sequenceLength;
+  double **D = 0;    // Differences query and target 
+  double **DD = 0;   // Matched filter distance
 
   D = new double*[numVectors];
-  assert(D);
   DD = new double*[numVectors];
-  assert(DD);
 
   gettimeofday(&tv1, NULL); 
   unsigned processedTracks = 0;
-  unsigned successfulTracks=0;
 
   // build track offset table
   off_t *trackOffsetTable = new off_t[dbH->numFiles];
@@ -378,14 +562,7 @@
     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;
+  char nextKey[MAXSTR];
 
   // Track loop 
   size_t data_buffer_size = 0;
@@ -411,160 +588,45 @@
     trackIndexOffset=trackOffset/dbH->dim; // numVectors offset
 
     read_data(track, &data_buffer, &data_buffer_size);
-    if(sequenceLength<=trackTable[track]){  // test for short sequences
+    if(sequenceLength <= trackTable[track]) {  // test for short sequences
       
-      if(verbosity>7) {
-	std::cerr << track << "." << trackIndexOffset << "." << trackTable[track] << " | ";std::cerr.flush();
-      }
+      VERB_LOG(7,"%u.%jd.%u | ", track, (intmax_t) trackIndexOffset, trackTable[track]);
       
       initialize_arrays(track, numVectors, query, data_buffer, D, DD);
 
-      if(verbosity>3 && usingTimes) {
-	std::cerr << "meanQdur=" << meanQdur << " meanDBdur=" << meanDBdur[track] << std::endl;
-	std::cerr.flush();
+      if(usingTimes) {
+        VERB_LOG(3,"meanQdur=%f meanDBdur=%f\n", meanQdur, meanDBdur[track]);
       }
 
-      if(!usingTimes || 
-	 (usingTimes 
-	  && fabs(meanDBdur[track]-meanQdur)<meanQdur*timesTol)){
-
-	if(verbosity>3 && usingTimes) {
-	  std::cerr << "within duration tolerance." << std::endl;
-	  std::cerr.flush();
-	}
+      if((!usingTimes) || fabs(meanDBdur[track]-meanQdur) < meanQdur*timesTol) {
+        if(usingTimes) {
+          VERB_LOG(3,"within duration tolerance.\n");
+        }
 
 	// Search for minimum distance by shingles (concatenated vectors)
-	for(j=0;j<=numVectors-wL;j+=HOP_SIZE)
+	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>9) {
-	      std::cerr << thisDist << " " << qnPtr[j] << " " << sNorm[trackIndexOffset+k] << std::endl;
+	    double thisDist=2-(2/(qnPtr[j]*sNorm[trackIndexOffset+k]))*DD[j][k];
+	    // Power test
+	    if ((!usingPower) || powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k])) {
+              // radius test
+              if((!radius) || thisDist < radius) {
+                reporter->add_point(track, usingQueryPoint ? queryPoint : j, k, thisDist);
+              }
             }
-	    // 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(qnPtr[j] - sNorm[trackIndexOffset+k]);
-	    // Power test
-	    if (usingPower) {
-	      if (!(powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k]))) {
-		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];
-		}
-		distances[m]=thisDist;
-		if(usingQueryPoint)
-		  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] == 1000000.0) break;
-	  thisDist+=distances[m];
+          }
         }
-	thisDist/=m;
-	
-	// Let's see the distances then...
-	if(verbosity>3) {
-	  std::cerr << fileTable+track*O2_FILETABLESIZE << " " << thisDist << std::endl;
-        }
-
-
-	// All the track stuff goes here
-	n=trackNN;
-	while(n--){
-	  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];
-		trackIDs[l]=trackIDs[l-1];
-	      }
-	      trackDistances[n]=thisDist;
-	      trackQIndexes[n]=qIndexes[0];
-	      trackSIndexes[n]=sIndexes[0];
-	      successfulTracks++;
-	      trackIDs[n]=track;
-	      break;
-	    }
-	  }
-	  else
-	    break;
-	}
       } // Duration match            
       delete_arrays(track, numVectors, D, DD);
     }
-    // per-track reset array values
-    for(unsigned k=0; k<pointNN; k++){
-      distances[k]=1.0e6;
-      qIndexes[k]=~0;
-      sIndexes[k]=~0;    
-    }
   }
 
   free(data_buffer);
 
   gettimeofday(&tv2,NULL);
-  if(verbosity>1) {
-    std::cerr << std::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" << std::endl;
-    std::cerr << "sampleCount: " << sampleCount << " sampleSum: " << sampleSum << " logSampleSum: " << logSampleSum 
-	 << " minSample: " << minSample << " maxSample: " << maxSample << std::endl;
-  }  
-  if(adbQueryResponse==0){
-    if(verbosity>1) {
-      std::cerr<<std::endl;
-    }
-    // Output answer
-    // Loop over nearest neighbours
-    for(k=0; k < std::min(trackNN,successfulTracks); k++)
-      std::cout << fileTable+trackIDs[k]*O2_FILETABLESIZE << " " << trackDistances[k] << " " 
-	   << trackQIndexes[k] << " " << trackSIndexes[k] << std::endl;
-  }
-  else{ // Process Web Services Query
-    int listLen = std::min(trackNN, processedTracks);
-    adbQueryResponse->result.__sizeRlist=listLen;
-    adbQueryResponse->result.__sizeDist=listLen;
-    adbQueryResponse->result.__sizeQpos=listLen;
-    adbQueryResponse->result.__sizeSpos=listLen;
-    adbQueryResponse->result.Rlist= new char*[listLen];
-    adbQueryResponse->result.Dist = new double[listLen];
-    adbQueryResponse->result.Qpos = new unsigned int[listLen];
-    adbQueryResponse->result.Spos = new unsigned int[listLen];
-    for(k=0; k<(unsigned)adbQueryResponse->result.__sizeRlist; k++){
-      adbQueryResponse->result.Rlist[k]=new char[O2_MAXFILESTR];
-      adbQueryResponse->result.Dist[k]=trackDistances[k];
-      adbQueryResponse->result.Qpos[k]=trackQIndexes[k];
-      adbQueryResponse->result.Spos[k]=trackSIndexes[k];
-      sprintf(adbQueryResponse->result.Rlist[k], "%s", fileTable+trackIDs[k]*O2_FILETABLESIZE);
-    }
-  }
+  VERB_LOG(1,"elapsed time: %ld msec\n",
+           (tv2.tv_sec*1000 + tv2.tv_usec/1000) - 
+           (tv1.tv_sec*1000 + tv1.tv_usec/1000))
 
   // Clean up
   if(trackOffsetTable)
@@ -577,267 +639,7 @@
     delete[] sNorm;
   if(qPower)
     delete[] qPower;
-  if(sPower) 
-    delete[] sPower;
-  if(D)
-    delete[] D;
-  if(DD)
-    delete[] DD;
-  if(meanDBdur)
-    delete[] meanDBdur;
-}
-
-void audioDB::trackSequenceQueryRad(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse){
-  
-  initTables(dbName, inFile);
-  
-  unsigned int numVectors;
-  double *query, *query_data;
-  double *qNorm, *qnPtr, *qPower = 0, *qpPtr = 0;
-  double meanQdur;
-
-  set_up_query(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors);
-
-  unsigned int dbVectors;
-  double *sNorm, *snPtr, *sPower = 0, *spPtr = 0;
-  double *meanDBdur = 0;
-
-  set_up_db(&sNorm, &snPtr, &sPower, &spPtr, &meanDBdur, &dbVectors);
-  
-  if(verbosity>1) {
-    std::cerr << "matching tracks..." << std::endl;
-  }
-  
-  assert(pointNN>0 && pointNN<=O2_MAXNN);
-  assert(trackNN>0 && trackNN<=O2_MAXNN);
-  
-  // Make temporary dynamic memory for results
-  double trackDistances[trackNN];
-  unsigned trackIDs[trackNN];
-  unsigned trackQIndexes[trackNN];
-  unsigned trackSIndexes[trackNN];
-  
-  double distances[pointNN];
-  unsigned qIndexes[pointNN];
-  unsigned sIndexes[pointNN];
-  
-  unsigned j,k,l,n,track,trackOffset=0;
-  unsigned const HOP_SIZE=sequenceHop;
-  unsigned const wL=sequenceLength;
-  double thisDist;
-  
-  for(k=0; k<pointNN; k++){
-    distances[k]=0.0;
-    qIndexes[k]=~0;
-    sIndexes[k]=~0;    
-  }
-  
-  for(k=0; k<trackNN; k++){
-    trackDistances[k]=0.0;
-    trackQIndexes[k]=~0;
-    trackSIndexes[k]=~0;
-    trackIDs[k]=~0;
-  }
-
-  double ** D = 0;    // Differences query and target 
-  double ** DD = 0;   // Matched filter distance
-
-  D = new double*[numVectors];
-  assert(D);
-  DD = new double*[numVectors];
-  assert(DD);
-
-  gettimeofday(&tv1, NULL); 
-  unsigned processedTracks = 0;
-  unsigned successfulTracks=0;
-
-  // build track offset table
-  off_t *trackOffsetTable = new off_t[dbH->numFiles];
-  unsigned cumTrack=0;
-  off_t trackIndexOffset;
-  for(k=0; k<dbH->numFiles;k++){
-    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 
-  size_t data_buffer_size = 0;
-  double *data_buffer = 0;
-  lseek(dbfid, dbH->dataOffset, SEEK_SET);
-
-  for(processedTracks=0, track=0 ; processedTracks < dbH->numFiles ; track++, processedTracks++){
-
-    trackOffset = trackOffsetTable[track];     // numDoubles offset
-
-    // get trackID from file if using a control file
-    if(trackFile) {
-      trackFile->getline(nextKey,MAXSTR);
-      if(!trackFile->eof()) {
-	track = getKeyPos(nextKey);
-        trackOffset = trackOffsetTable[track];
-        lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
-      } else {
-	break;
-      }
-    }
-
-    trackIndexOffset=trackOffset/dbH->dim; // numVectors offset
-
-    read_data(track, &data_buffer, &data_buffer_size);
-    if(sequenceLength<=trackTable[track]){  // test for short sequences
-      
-      if(verbosity>7) {
-	std::cerr << track << "." << trackIndexOffset << "." << trackTable[track] << " | ";std::cerr.flush();
-      }
-
-      initialize_arrays(track, numVectors, query, data_buffer, D, DD);
-
-      if(verbosity>3 && usingTimes) {
-	std::cerr << "meanQdur=" << meanQdur << " meanDBdur=" << meanDBdur[track] << std::endl;
-	std::cerr.flush();
-      }
-
-      if(!usingTimes || 
-	 (usingTimes 
-	  && fabs(meanDBdur[track]-meanQdur)<meanQdur*timesTol)){
-
-	if(verbosity>3 && usingTimes) {
-	  std::cerr << "within duration tolerance." << std::endl;
-	  std::cerr.flush();
-	}
-
-	// 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>9) {
-	      std::cerr << thisDist << " " << qnPtr[j] << " " << sNorm[trackIndexOffset+k] << std::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(qnPtr[j] - sNorm[trackIndexOffset+k]);
-	    // Power test
-            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
-	    }
-	  }
-	// How many points were below threshold ?
-	thisDist=distances[0];
-	
-	// Let's see the distances then...
-	if(verbosity>3) {
-	  std::cerr << fileTable+track*O2_FILETABLESIZE << " " << thisDist << std::endl;
-        }
-
-	// All the track stuff goes here
-	n=trackNN;
-	while(n--){
-	  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];
-		trackIDs[l]=trackIDs[l-1];
-	      }
-	      trackDistances[n]=thisDist;
-	      trackQIndexes[n]=qIndexes[0];
-	      trackSIndexes[n]=sIndexes[0];
-	      successfulTracks++;
-	      trackIDs[n]=track;
-	      break;
-	    }
-	  }
-	  else
-	    break;
-	}
-      } // Duration match
-      delete_arrays(track, numVectors, D, DD);
-    }
-    // per-track reset array values
-    for(unsigned k=0; k<pointNN; k++){
-      distances[k]=0.0;
-      qIndexes[k]=~0;
-      sIndexes[k]=~0;    
-    }
-  }
-
-  free(data_buffer);
-
-  gettimeofday(&tv2,NULL);
-  if(verbosity>1) {
-    std::cerr << std::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" << std::endl;
-    std::cerr << "sampleCount: " << sampleCount << " sampleSum: " << sampleSum << " logSampleSum: " << logSampleSum 
-	 << " minSample: " << minSample << " maxSample: " << maxSample << std::endl;
-  }
-  
-  if(adbQueryResponse==0){
-    if(verbosity>1) {
-      std::cerr<<std::endl;
-    }
-    // Output answer
-    // Loop over nearest neighbours
-    for(k=0; k < std::min(trackNN,successfulTracks); k++)
-      std::cout << fileTable+trackIDs[k]*O2_FILETABLESIZE << " " << trackDistances[k] << std::endl;
-  }
-  else{ // Process Web Services Query
-    int listLen = std::min(trackNN, processedTracks);
-    adbQueryResponse->result.__sizeRlist=listLen;
-    adbQueryResponse->result.__sizeDist=listLen;
-    adbQueryResponse->result.__sizeQpos=listLen;
-    adbQueryResponse->result.__sizeSpos=listLen;
-    adbQueryResponse->result.Rlist= new char*[listLen];
-    adbQueryResponse->result.Dist = new double[listLen];
-    adbQueryResponse->result.Qpos = new unsigned int[listLen];
-    adbQueryResponse->result.Spos = new unsigned int[listLen];
-    for(k=0; k<(unsigned)adbQueryResponse->result.__sizeRlist; k++){
-      adbQueryResponse->result.Rlist[k]=new char[O2_MAXFILESTR];
-      adbQueryResponse->result.Dist[k]=trackDistances[k];
-      adbQueryResponse->result.Qpos[k]=trackQIndexes[k];
-      adbQueryResponse->result.Spos[k]=trackSIndexes[k];
-      sprintf(adbQueryResponse->result.Rlist[k], "%s", fileTable+trackIDs[k]*O2_FILETABLESIZE);
-    }
-  }
-
-  // Clean up
-  if(trackOffsetTable)
-    delete[] trackOffsetTable;
-  if(query_data)
-    delete[] query_data;
-  if(qNorm)
-    delete[] qNorm;
-  if(sNorm)
-    delete[] sNorm;
-  if(qPower)
-    delete[] qPower;
-  if(sPower) 
+  if(sPower)
     delete[] sPower;
   if(D)
     delete[] D;
@@ -851,9 +653,8 @@
 void audioDB::unitNorm(double* X, unsigned dim, unsigned n, double* qNorm){
   unsigned d;
   double L2, *p;
-  if(verbosity>2) {
-    std::cerr << "norming " << n << " vectors...";std::cerr.flush();
-  }
+
+  VERB_LOG(2, "norming %u vectors...", n);
   while(n--) {
     p = X;
     L2 = 0.0;
@@ -867,7 +668,5 @@
     }
     X += dim;
   }
-  if(verbosity>2) {
-    std::cerr << "done..." << std::endl;
-  }
+  VERB_LOG(2, "done.\n");
 }