view query.cpp @ 457:823bca1e10f5 api-inversion

Sketch out a "query state" structure.

As yet it's completely unused, but the intention is that accumulated
state will be collected into one of these structures for each query, and
then passed around, to help reduce the need for silly arglists.  It's
possible that this structure will also grow a pointer to the adb itself,
and be the thing passed to the LSH callback; we'll see how that
develops.

#include "audioDB.h"
#include "reporter.h"

#include "audioDB-internals.h"
#include "accumulators.h"

bool audiodb_powers_acceptable(adb_query_refine_t *r, double p1, double p2) {
  if (r->flags & ADB_REFINE_ABSOLUTE_THRESHOLD) {
    if ((p1 < r->absolute_threshold) || (p2 < r->absolute_threshold)) {
      return false;
    }
  }
  if (r->flags & ADB_REFINE_RELATIVE_THRESHOLD) {
    if (fabs(p1-p2) > fabs(r->relative_threshold)) {
      return false;
    }
  }
  return true;
}

void audioDB::query(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse) {

  // init database tables and dbH first
  if(query_from_key)
    initTables(dbName);
  else
    initTables(dbName, inFile);

  adb_query_spec_t qspec;
  adb_datum_t datum = {0};

  qspec.refine.flags = 0;
  if(trackFile) {
    qspec.refine.flags |= ADB_REFINE_INCLUDE_KEYLIST;
    std::vector<const char *> v;
    char *k = new char[MAXSTR];
    trackFile->getline(k, MAXSTR);    
    while(!trackFile->eof()) {
      v.push_back(k);
      k = new char[MAXSTR];
      trackFile->getline(k, MAXSTR);    
    }
    delete [] k;
    qspec.refine.include.nkeys = v.size();
    qspec.refine.include.keys = new const char *[qspec.refine.include.nkeys];
    for(unsigned int k = 0; k < qspec.refine.include.nkeys; k++) {
      qspec.refine.include.keys[k] = v[k];
    }
  }
  if(query_from_key) {
    qspec.refine.flags |= ADB_REFINE_EXCLUDE_KEYLIST;
    qspec.refine.exclude.nkeys = 1;
    qspec.refine.exclude.keys = &key;
  }
  if(radius) {
    qspec.refine.flags |= ADB_REFINE_RADIUS;
    qspec.refine.radius = radius;
  }
  if(use_absolute_threshold) {
    qspec.refine.flags |= ADB_REFINE_ABSOLUTE_THRESHOLD;
    qspec.refine.absolute_threshold = absolute_threshold;
  }
  if(use_relative_threshold) {
    qspec.refine.flags |= ADB_REFINE_RELATIVE_THRESHOLD;
    qspec.refine.relative_threshold = relative_threshold;
  }
  if(usingTimes) {
    qspec.refine.flags |= ADB_REFINE_DURATION_RATIO;
    qspec.refine.duration_ratio = timesTol;
  }
  /* FIXME: not sure about this any more; maybe it belongs in
     query_id?  Or maybe we just don't need a flag for it? */
  qspec.refine.hopsize = sequenceHop;
  if(sequenceHop != 1) {
    qspec.refine.flags |= ADB_REFINE_HOP_SIZE;
  }

  if(query_from_key) {
    datum.key = key;
  } else {
    int fd;
    struct stat st;

    /* FIXME: around here there are all sorts of hideous leaks. */
    fd = open(inFile, O_RDONLY);
    if(fd < 0) {
      error("failed to open feature file", inFile);
    }
    fstat(fd, &st);
    read(fd, &datum.dim, sizeof(uint32_t));
    datum.nvectors = (st.st_size - sizeof(uint32_t)) / (datum.dim * sizeof(double));
    datum.data = (double *) malloc(st.st_size - sizeof(uint32_t));
    read(fd, datum.data, st.st_size - sizeof(uint32_t));
    close(fd);
    if(usingPower) {
      uint32_t one;
      fd = open(powerFileName, O_RDONLY);
      if(fd < 0) {
        error("failed to open power file", powerFileName);
      }
      read(fd, &one, sizeof(uint32_t));
      if(one != 1) {
        error("malformed power file dimensionality", powerFileName);
      }
      datum.power = (double *) malloc(datum.nvectors * sizeof(double));
      if(read(fd, datum.power, datum.nvectors * sizeof(double)) != (ssize_t) (datum.nvectors * sizeof(double))) {
        error("malformed power file", powerFileName);
      }
      close(fd);
    }
    if(usingTimes) {
      datum.times = (double *) malloc(2 * datum.nvectors * sizeof(double));
      insertTimeStamps(datum.nvectors, timesFile, datum.times);
    }
  }

  qspec.qid.datum = &datum;
  qspec.qid.sequence_length = sequenceLength;
  qspec.qid.flags = usingQueryPoint ? 0 : ADB_QUERY_ID_FLAG_EXHAUSTIVE;
  qspec.qid.sequence_start = queryPoint;

  switch(queryType) {
  case O2_POINT_QUERY:
    qspec.qid.sequence_length = 1;
    qspec.params.accumulation = ADB_ACCUMULATION_DB;
    qspec.params.distance = ADB_DISTANCE_DOT_PRODUCT;
    qspec.params.npoints = pointNN;
    qspec.params.ntracks = 0;
    reporter = new pointQueryReporter< std::greater < NNresult > >(pointNN);
    break;
  case O2_TRACK_QUERY:
    qspec.qid.sequence_length = 1;
    qspec.params.accumulation = ADB_ACCUMULATION_PER_TRACK;
    qspec.params.distance = ADB_DISTANCE_DOT_PRODUCT;
    qspec.params.npoints = pointNN;
    qspec.params.ntracks = trackNN;
    reporter = new trackAveragingReporter< std::greater< NNresult > >(pointNN, trackNN, dbH->numFiles);
    break;
  case O2_SEQUENCE_QUERY:
  case O2_N_SEQUENCE_QUERY:
    qspec.params.accumulation = ADB_ACCUMULATION_PER_TRACK;
    qspec.params.distance = no_unit_norming ? ADB_DISTANCE_EUCLIDEAN : ADB_DISTANCE_EUCLIDEAN_NORMED;
    qspec.params.npoints = pointNN;
    qspec.params.ntracks = trackNN;
    switch(queryType) {
    case O2_SEQUENCE_QUERY:
      if(!(qspec.refine.flags & ADB_REFINE_RADIUS)) {
        reporter = new trackAveragingReporter< std::less< NNresult > >(pointNN, trackNN, dbH->numFiles);
      } else if (index_exists(adb->path, qspec.refine.radius, qspec.qid.sequence_length)) {
	char* indexName = index_get_name(adb->path, qspec.refine.radius, qspec.qid.sequence_length);
	lsh = index_allocate(indexName, false);
	reporter = new trackSequenceQueryRadReporter(trackNN, index_to_trackID(lsh->get_maxp(), lsh_n_point_bits)+1);
	delete[] indexName;
      } else {
	reporter = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles);
      }
      break;
    case O2_N_SEQUENCE_QUERY:
      if(!(qspec.refine.flags & ADB_REFINE_RADIUS)) {
        reporter = new trackSequenceQueryNNReporter< std::less < NNresult > >(pointNN, trackNN, dbH->numFiles);
      } else if (index_exists(adb->path, qspec.refine.radius, qspec.qid.sequence_length)){
	char* indexName = index_get_name(adb->path, qspec.refine.radius, qspec.qid.sequence_length);
	lsh = index_allocate(indexName, false);
	reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, index_to_trackID(lsh->get_maxp(), lsh_n_point_bits)+1);
	delete[] indexName;
      } else {
	reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, dbH->numFiles);
      }
      break;
    }
    break;
  case O2_ONE_TO_ONE_N_SEQUENCE_QUERY:
    qspec.params.accumulation = ADB_ACCUMULATION_ONE_TO_ONE;
    qspec.params.distance = ADB_DISTANCE_EUCLIDEAN_NORMED;
    qspec.params.npoints = 0;
    qspec.params.ntracks = 0;
    break;
  default:
    error("unrecognized queryType");
  }

  // keyKeyPos requires dbH to be initialized
  if(query_from_key && (!key || (query_from_key_index = audiodb_key_index(adb, key)) == (uint32_t) -1)) 
    error("Query key not found", key);  

  switch(qspec.params.distance) {
  case ADB_DISTANCE_DOT_PRODUCT:
    switch(qspec.params.accumulation) {
    case ADB_ACCUMULATION_DB:
      accumulator = new DBAccumulator<adb_result_dist_gt>(qspec.params.npoints);
      break;
    case ADB_ACCUMULATION_PER_TRACK:
      accumulator = new PerTrackAccumulator<adb_result_dist_gt>(qspec.params.npoints, qspec.params.ntracks);
      break;
    case ADB_ACCUMULATION_ONE_TO_ONE:
      accumulator = new NearestAccumulator<adb_result_dist_gt>();
      break;
    default:
      error("unknown accumulation");
    }
    break;
  case ADB_DISTANCE_EUCLIDEAN_NORMED:
  case ADB_DISTANCE_EUCLIDEAN:
    switch(qspec.params.accumulation) {
    case ADB_ACCUMULATION_DB:
      accumulator = new DBAccumulator<adb_result_dist_lt>(qspec.params.npoints);
      break;
    case ADB_ACCUMULATION_PER_TRACK:
      accumulator = new PerTrackAccumulator<adb_result_dist_lt>(qspec.params.npoints, qspec.params.ntracks);
      break;
    case ADB_ACCUMULATION_ONE_TO_ONE:
      accumulator = new NearestAccumulator<adb_result_dist_lt>();
      break;
    default:
      error("unknown accumulation");
    }
    break;
  default:
    error("unknown distance function");
  }
  
  // Test for index (again) here
  if((qspec.refine.flags & ADB_REFINE_RADIUS) && index_exists(adb->path, qspec.refine.radius, qspec.qid.sequence_length)){ 
    VERB_LOG(1, "Calling indexed query on database %s, radius=%f, sequence_length=%d\n", adb->path, qspec.refine.radius, qspec.qid.sequence_length);
    index_query_loop(adb, &qspec);
  }
  else{
    VERB_LOG(1, "Calling brute-force query on database %s\n", dbName);
    if(query_loop(adb, &qspec)) {
      error("query_loop failed");
    }
  }

  adb_query_results_t *rs = accumulator->get_points();
  for(unsigned int k = 0; k < rs->nresults; k++) {
    adb_result_t r = rs->results[k];
    reporter->add_point(audiodb_key_index(adb, r.key), r.qpos, r.ipos, r.dist);
  }

  reporter->report(fileTable, adbQueryResponse);
}

static void audiodb_initialize_arrays(adb_t *adb, adb_query_spec_t *spec, int track, unsigned int numVectors, double *query, double *data_buffer, double **D, double **DD) {
  unsigned int j, k, l, w;
  double *dp, *qp, *sp;

  const unsigned HOP_SIZE = spec->refine.hopsize;
  const unsigned wL = spec->qid.sequence_length;

  for(j = 0; j < numVectors; j++) {
    // Sum products matrix
    D[j] = new double[(*adb->track_lengths)[track]]; 
    assert(D[j]);
    // Matched filter matrix
    DD[j]=new double[(*adb->track_lengths)[track]];
    assert(DD[j]);
  }

  // Dot product
  for(j = 0; j < numVectors; j++)
    for(k = 0; k < (*adb->track_lengths)[track]; k++){
      qp = query + j * adb->header->dim;
      sp = data_buffer + k * adb->header->dim;
      DD[j][k] = 0.0; // Initialize matched filter array
      dp = &D[j][k];  // point to correlation cell j,k
      *dp = 0.0;      // initialize correlation cell
      l = adb->header->dim;         // size of vectors
      while(l--)
        *dp += *qp++ * *sp++;
    }
  
  // Matched Filter
  // HOP SIZE == 1
  double* spd;
  if(HOP_SIZE == 1) { // HOP_SIZE = shingleHop
    for(w = 0; w < wL; w++) {
      for(j = 0; j < numVectors - w; j++) { 
        sp = DD[j];
        spd = D[j+w] + w;
        k = (*adb->track_lengths)[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;
        for(k = 0; k < (*adb->track_lengths)[track] - w; k += HOP_SIZE) {
          *sp += *spd;
          sp += HOP_SIZE;
          spd += HOP_SIZE;
        }
      }
    }
  }
}

static void audiodb_delete_arrays(int track, unsigned int numVectors, double **D, double **DD) {
  if(D != NULL) {
    for(unsigned int j = 0; j < numVectors; j++) {
      delete[] D[j];
    }
  }
  if(DD != NULL) {
    for(unsigned int j = 0; j < numVectors; j++) {
      delete[] DD[j];
    }
  }
}

int audiodb_read_data(adb_t *adb, int trkfid, int track, double **data_buffer_p, size_t *data_buffer_size_p) {
  uint32_t track_length = (*adb->track_lengths)[track];
  size_t track_size = track_length * sizeof(double) * adb->header->dim;
  if (track_size > *data_buffer_size_p) {
    if(*data_buffer_p) {
      free(*data_buffer_p);
    }
    { 
      *data_buffer_size_p = track_size;
      void *tmp = malloc(track_size);
      if (tmp == NULL) {
        goto error;
      }
      *data_buffer_p = (double *) tmp;
    }
  }

  read_or_goto_error(trkfid, *data_buffer_p, track_size);
  return 0;

 error:
  return 1;
}

void audioDB::insertTimeStamps(unsigned numVectors, std::ifstream *timesFile, double *timesdata) {
  assert(usingTimes);

  unsigned numtimes = 0;

  if(!timesFile->is_open()) {
    error("problem opening times file on timestamped database", timesFileName);
  }

  double timepoint, next;
  *timesFile >> timepoint;
  if (timesFile->eof()) {
    error("no entries in times file", timesFileName);
  }
  numtimes++;
  do {
    *timesFile >> next;
    if (timesFile->eof()) {
      break;
    }
    numtimes++;
    timesdata[0] = timepoint;
    timepoint = (timesdata[1] = next);
    timesdata += 2;
  } while (numtimes < numVectors + 1);

  if (numtimes < numVectors + 1) {
    error("too few timepoints in times file", timesFileName);
  }

  *timesFile >> next;
  if (!timesFile->eof()) {
    error("too many timepoints in times file", timesFileName);
  }
}

int audiodb_query_spec_qpointers(adb_t *adb, adb_query_spec_t *spec, double **vector_data, double **vector, adb_qpointers_internal_t *qpointers) {
  adb_datum_t *datum;
  adb_datum_t d = {0};
  uint32_t nvectors;
  uint32_t sequence_length;
  uint32_t sequence_start;

  datum = spec->qid.datum;
  sequence_length = spec->qid.sequence_length;
  sequence_start = spec->qid.sequence_start;
    
  if(datum->data) {
    if(datum->dim != adb->header->dim) {
      return 1;
    }
    /* initialize d, and mark that nothing needs freeing later. */
    d = *datum;
    datum = &d;
  } else if (datum->key) {
    uint32_t track_id;
    if((track_id = audiodb_key_index(adb, datum->key)) == (uint32_t) -1) {
      return 1;
    }
    off_t track_offset = (*adb->track_offsets)[track_id];
    
    if(adb->header->flags & O2_FLAG_LARGE_ADB) {
      /* create a reference/insert, then use adb_insert_create_datum() */
      adb_reference_t reference = {0};
      char features[MAXSTR], power[MAXSTR], times[MAXSTR];
      lseek(adb->fd, adb->header->dataOffset + track_id * O2_FILETABLE_ENTRY_SIZE, SEEK_SET);
      /* FIXME: learn not to worry and love the bomb^Wbuffer overflow */
      read(adb->fd, features, MAXSTR);
      reference.features = features;
      if(adb->header->flags & O2_FLAG_POWER) {
        lseek(adb->fd, adb->header->powerTableOffset + track_id * O2_FILETABLE_ENTRY_SIZE, SEEK_SET);
        read(adb->fd, power, MAXSTR);
        reference.power = power;
      }
      if(adb->header->flags & O2_FLAG_TIMES) {
        lseek(adb->fd, adb->header->timesTableOffset + track_id * O2_FILETABLE_ENTRY_SIZE, SEEK_SET);
        read(adb->fd, times, MAXSTR);
        reference.times = times;
      }
      audiodb_insert_create_datum(&reference, &d);
    } else {
      /* initialize from sources of data that we already have */
      d.nvectors = (*adb->track_lengths)[track_id];
      d.dim = adb->header->dim;
      d.key = datum->key;
      /* read out stuff from the database tables */
      d.data = (double *) malloc(d.nvectors * d.dim * sizeof(double));
      lseek(adb->fd, adb->header->dataOffset + track_offset, SEEK_SET);
      read(adb->fd, d.data, d.nvectors * d.dim * sizeof(double));
      if(adb->header->flags & O2_FLAG_POWER) {
        d.power = (double *) malloc(d.nvectors * sizeof(double));
        lseek(adb->fd, adb->header->powerTableOffset + track_offset / d.dim, SEEK_SET);
        read(adb->fd, d.power, d.nvectors * sizeof(double));
      }
      if(adb->header->flags & O2_FLAG_TIMES) {
        d.times = (double *) malloc(2 * d.nvectors * sizeof(double));
        lseek(adb->fd, adb->header->timesTableOffset + track_offset / d.dim, SEEK_SET);
        read(adb->fd, d.times, 2 * d.nvectors * sizeof(double));
      }
    }
  } else {
    return 1;
  }

  /* Now we have a full(ish) datum, compute all the qpointery stuff
     that we care about (l2norm/power/mean duration).  (This bit could
     conceivably become a new function) */
  nvectors = d.nvectors;
  /* FIXME: check the overflow logic here */
  if(sequence_start + sequence_length > nvectors) {
    /* is there something to free?  goto error */
    return 1;
  }

  qpointers->nvectors = nvectors;

  size_t vector_size = nvectors * sizeof(double) * d.dim;
  *vector_data = new double[vector_size];
  memcpy(*vector_data, d.data, vector_size);

  qpointers->l2norm_data = new double[vector_size / d.dim];
  audiodb_l2norm_buffer(*vector_data, d.dim, nvectors, qpointers->l2norm_data);
  audiodb_sequence_sum(qpointers->l2norm_data, nvectors, sequence_length);
  audiodb_sequence_sqrt(qpointers->l2norm_data, nvectors, sequence_length);

  if(d.power) {
    qpointers->power_data = new double[vector_size / d.dim];
    memcpy(qpointers->power_data, d.power, vector_size / d.dim);
    audiodb_sequence_sum(qpointers->power_data, nvectors, sequence_length);
    audiodb_sequence_average(qpointers->power_data, nvectors, sequence_length);
  }

  if(d.times) {
    qpointers->mean_duration = new double[1];
    *qpointers->mean_duration = 0;
    for(unsigned int k = 0; k < nvectors; k++) {
      *qpointers->mean_duration += d.times[2*k+1] - d.times[2*k];
    }
    *qpointers->mean_duration /= nvectors;
  }


  /* Finally, set up the moving qpointers. */
  if(spec->qid.flags & ADB_QUERY_ID_FLAG_EXHAUSTIVE) {
    *vector = *vector_data;
    qpointers->l2norm = qpointers->l2norm_data;
    qpointers->power = qpointers->power_data;
  } else {
    *vector = *vector_data + spec->qid.sequence_start * d.dim;
    qpointers->l2norm = qpointers->l2norm_data + spec->qid.sequence_start;
    if(d.power) {
      qpointers->power = qpointers->power_data + spec->qid.sequence_start;
    }
    /* FIXME: this is a little bit ugly.  No, a lot ugly.  But at the
     * moment this is how query_loop() knows when to stop, so for
     * now... */
    qpointers->nvectors = sequence_length;
  }

  /* Clean up: free any bits of datum that we have ourselves
   * allocated. */
  if(datum != &d) {
    audiodb_free_datum(&d);
  }

  return 0;
}

static int audiodb_set_up_dbpointers(adb_t *adb, adb_query_spec_t *spec, adb_qpointers_internal_t *dbpointers) {
  uint32_t nvectors = adb->header->length / (adb->header->dim * sizeof(double));
  uint32_t sequence_length = spec->qid.sequence_length;

  bool using_power = spec->refine.flags & (ADB_REFINE_ABSOLUTE_THRESHOLD|ADB_REFINE_RELATIVE_THRESHOLD);
  bool using_times = spec->refine.flags & ADB_REFINE_DURATION_RATIO;
  double *times_table = NULL;


  dbpointers->nvectors = nvectors;
  dbpointers->l2norm_data = new double[nvectors];

  double *snpp = dbpointers->l2norm_data, *sppp = 0;
  lseek(adb->fd, adb->header->l2normTableOffset, SEEK_SET);
  read_or_goto_error(adb->fd, dbpointers->l2norm_data, nvectors * sizeof(double));

  if (using_power) {
    if (!(adb->header->flags & O2_FLAG_POWER)) {
      goto error;
    }
    dbpointers->power_data = new double[nvectors];
    sppp = dbpointers->power_data;
    lseek(adb->fd, adb->header->powerTableOffset, SEEK_SET);
    read_or_goto_error(adb->fd, dbpointers->power_data, nvectors * sizeof(double));
  }

  for(unsigned int i = 0; i < adb->header->numFiles; i++){
    size_t track_length = (*adb->track_lengths)[i];
    if(track_length >= sequence_length) {
      audiodb_sequence_sum(snpp, track_length, sequence_length);
      audiodb_sequence_sqrt(snpp, track_length, sequence_length);
      if (using_power) {
	audiodb_sequence_sum(sppp, track_length, sequence_length);
        audiodb_sequence_average(sppp, track_length, sequence_length);
      }
    }
    snpp += track_length;
    if (using_power) {
      sppp += track_length;
    }
  }

  if (using_times) {
    if(!(adb->header->flags & O2_FLAG_TIMES)) {
      goto error;
    }

    dbpointers->mean_duration = new double[adb->header->numFiles];

    times_table = (double *) malloc(2 * nvectors * sizeof(double));
    if(!times_table) {
      goto error;
    }
    lseek(adb->fd, adb->header->timesTableOffset, SEEK_SET);
    read_or_goto_error(adb->fd, times_table, 2 * nvectors * sizeof(double));
    for(unsigned int k = 0; k < adb->header->numFiles; k++) {
      size_t track_length = (*adb->track_lengths)[k];
      unsigned int j;
      dbpointers->mean_duration[k] = 0.0;
      for(j = 0; j < track_length; j++) {
	dbpointers->mean_duration[k] += times_table[2*j+1] - times_table[2*j];
      }
      dbpointers->mean_duration[k] /= j;
    }

    free(times_table);
    times_table = NULL;
  }

  dbpointers->l2norm = dbpointers->l2norm_data;
  dbpointers->power = dbpointers->power_data;
  return 0;

 error:
  if(dbpointers->l2norm_data) {
    delete [] dbpointers->l2norm_data;
  }
  if(dbpointers->power_data) {
    delete [] dbpointers->power_data;
  }
  if(dbpointers->mean_duration) {
    delete [] dbpointers->mean_duration;
  }
  if(times_table) {
    free(times_table);
  }
  return 1;
  
}

// query_points()
//
// using PointPairs held in the exact_evaluation_queue compute squared distance for each PointPair
// and insert result into the current reporter.
//
// Preconditions:
// A query inFile has been opened with setup_query(...) and query pointers initialized
// The database contains some points
// An exact_evaluation_queue has been allocated and populated
// A reporter has been allocated
//
// Postconditions:
// reporter contains the points and distances that meet the reporter constraints 

void audioDB::query_loop_points(adb_t *adb, adb_query_spec_t *spec, double* query, adb_qpointers_internal_t *qpointers) {
  adb_qpointers_internal_t dbpointers = {0};

  uint32_t sequence_length = spec->qid.sequence_length;
  bool power_refine = spec->refine.flags & (ADB_REFINE_ABSOLUTE_THRESHOLD|ADB_REFINE_RELATIVE_THRESHOLD);

  if(exact_evaluation_queue->size() == 0) {
    return;
  }

  // Compute database info.  FIXME: we more than likely don't need
  // very much of the database so write a new function to build these
  // values per-track or, even better, per-point
  if(!(adb->header->flags & O2_FLAG_LARGE_ADB)) {
    if(audiodb_set_up_dbpointers(adb, spec, &dbpointers)) {
      error("failed to set up dbpointers");
    }
  }

  // We are guaranteed that the order of points is sorted by:
  // {trackID, spos, qpos} so we can be relatively efficient in
  // initialization of track data.  We assume that points usually
  // don't overlap, so we will use exhaustive dot product evaluation
  // instead of memoization of partial sums which is used for
  // exhaustive brute-force evaluation from smaller databases in
  // query_loop()
  double dist;
  size_t data_buffer_size = 0;
  double *data_buffer = 0;
  Uns32T trackOffset = 0;
  Uns32T trackIndexOffset = 0;
  Uns32T currentTrack = 0x80000000; // Initialize with a value outside of track index range
  Uns32T npairs = exact_evaluation_queue->size();
  while(npairs--){
    PointPair pp = exact_evaluation_queue->top();
    // Large ADB track data must be loaded here for sPower
    if(adb->header->flags & O2_FLAG_LARGE_ADB) {
      trackOffset=0;
      trackIndexOffset=0;
      if(currentTrack!=pp.trackID){
	char* prefixedString = new char[O2_MAXFILESTR];
	char* tmpStr = prefixedString;
	// On currentTrack change, allocate and load track data
	currentTrack=pp.trackID;
	SAFE_DELETE_ARRAY(dbpointers.l2norm_data);
	SAFE_DELETE_ARRAY(dbpointers.power_data);
	if(infid>0)
	  close(infid);
	// Open and check dimensions of feature file
	strncpy(prefixedString, featureFileNameTable+pp.trackID*O2_FILETABLE_ENTRY_SIZE, O2_MAXFILESTR);
	prefix_name((char ** const) &prefixedString, adb_feature_root);
	if (prefixedString!=tmpStr)
	  delete[] tmpStr;
	initInputFile(prefixedString);
	// Load the feature vector data for current track into data_buffer
	if(audiodb_read_data(adb, infid, pp.trackID, &data_buffer, &data_buffer_size))
          error("failed to read data");
	// Load power and calculate power and l2norm sequence sums
	init_track_aux_data(pp.trackID, data_buffer, &dbpointers.l2norm_data, &dbpointers.l2norm, &dbpointers.power_data, &dbpointers.power);
      }
    }
    else{
      // These offsets are w.r.t. the entire database of feature vectors and auxillary variables
      trackOffset = (*adb->track_offsets)[pp.trackID];
      trackIndexOffset = trackOffset/(adb->header->dim * sizeof(double)); // num vectors offset
    }    
    Uns32T qPos = usingQueryPoint?0:pp.qpos;// index for query point
    Uns32T sPos = trackIndexOffset+pp.spos; // index into l2norm table
    // Test power thresholds before computing distance
    if( ( (!power_refine) || audiodb_powers_acceptable(&spec->refine, qpointers->power[qPos], dbpointers.power[sPos])) &&
	( qPos<qpointers->nvectors-sequence_length+1 && pp.spos<(*adb->track_lengths)[pp.trackID]-sequence_length+1 ) ){
      // Non-large ADB track data is loaded inside power test for efficiency
      if(!(adb->header->flags & O2_FLAG_LARGE_ADB) && (currentTrack!=pp.trackID) ){
	// On currentTrack change, allocate and load track data
	currentTrack=pp.trackID;
	lseek(dbfid, adb->header->dataOffset + trackOffset, SEEK_SET);
	if(audiodb_read_data(adb, dbfid, currentTrack, &data_buffer, &data_buffer_size))
          error("failed to read data");
      }
      // Compute distance
      dist = audiodb_dot_product(query + qPos*adb->header->dim, data_buffer + pp.spos*adb->header->dim, adb->header->dim*sequence_length);
      double qn = qpointers->l2norm[qPos];
      double sn = dbpointers.l2norm[sPos];
      switch(spec->params.distance) {
      case ADB_DISTANCE_EUCLIDEAN_NORMED:
	dist = 2 - (2/(qn*sn))*dist;
        break;
      case ADB_DISTANCE_EUCLIDEAN:
        dist = qn*qn + sn*sn - 2*dist;
        break;
      }
      if((!radius) || dist <= (O2_LSH_EXACT_MULT*radius+O2_DISTANCE_TOLERANCE)) {
        adb_result_t r;
        r.key = (*adb->keys)[pp.trackID].c_str();
        r.dist = dist;
        r.qpos = pp.qpos;
        r.ipos = pp.spos;
        accumulator->add_point(&r);
      }
    }
    exact_evaluation_queue->pop();
  }
  // Cleanup
  SAFE_DELETE_ARRAY(dbpointers.l2norm_data);
  SAFE_DELETE_ARRAY(dbpointers.power_data);
  SAFE_DELETE_ARRAY(dbpointers.mean_duration);
}

int audioDB::query_loop(adb_t *adb, adb_query_spec_t *spec) {
  
  double *query, *query_data;
  adb_qpointers_internal_t qpointers = {0}, dbpointers = {0};

  bool power_refine = spec->refine.flags & (ADB_REFINE_ABSOLUTE_THRESHOLD|ADB_REFINE_RELATIVE_THRESHOLD);

  std::set<std::string> keys;
  if(spec->refine.flags & ADB_REFINE_INCLUDE_KEYLIST) {
    for(unsigned int k = 0; k < spec->refine.include.nkeys; k++) {
      keys.insert(spec->refine.include.keys[k]);
    }
  } else {
    for(unsigned int k = 0; k < adb->header->numFiles; k++) {
      keys.insert((*adb->keys)[k]);
    }
  }
  if(spec->refine.flags & ADB_REFINE_EXCLUDE_KEYLIST) {
    for(unsigned int k = 0; k < spec->refine.exclude.nkeys; k++) {
      keys.erase(spec->refine.exclude.keys[k]);
    }
  }

  if(adb->header->flags & O2_FLAG_LARGE_ADB) {
    /* FIXME: actually it would be nice to support this mode of
     * operation, but for now... */
    return 1;
  }

  if(audiodb_query_spec_qpointers(adb, spec, &query_data, &query, &qpointers)) {
    return 1;
  }

  if(audiodb_set_up_dbpointers(adb, spec, &dbpointers)) {
    return 1;
  }

  unsigned j,k,track,trackOffset=0, HOP_SIZE = spec->refine.hopsize;
  unsigned wL = spec->qid.sequence_length;
  double **D = 0;    // Differences query and target 
  double **DD = 0;   // Matched filter distance

  D = new double*[qpointers.nvectors]; // pre-allocate 
  DD = new double*[qpointers.nvectors];

  off_t trackIndexOffset;

  // Track loop 
  size_t data_buffer_size = 0;
  double *data_buffer = 0;
  lseek(adb->fd, adb->header->dataOffset, SEEK_SET);

  for(track = 0; track < adb->header->numFiles; track++) {
    unsigned t = track;

    while (keys.find((*adb->keys)[track]) == keys.end()) {
      track++;
      if(track == adb->header->numFiles) {
        goto loop_finish;
      }
    }
    trackOffset = (*adb->track_offsets)[track];
    if(track != t) {
      lseek(adb->fd, adb->header->dataOffset + trackOffset, SEEK_SET);
    }
    trackIndexOffset = trackOffset / (adb->header->dim * sizeof(double)); // dbpointers.nvectors offset

    if(audiodb_read_data(adb, adb->fd, track, &data_buffer, &data_buffer_size)) {
      return 1;
    }
    if(wL <= (*adb->track_lengths)[track]) {  // test for short sequences
      
      audiodb_initialize_arrays(adb, spec, track, qpointers.nvectors, query, data_buffer, D, DD);

      if((!(spec->refine.flags & ADB_REFINE_DURATION_RATIO)) || 
         fabs(dbpointers.mean_duration[track]-qpointers.mean_duration[0]) < qpointers.mean_duration[0]*spec->refine.duration_ratio) {

	// Search for minimum distance by shingles (concatenated vectors)
	for(j = 0; j <= qpointers.nvectors - wL; j += HOP_SIZE) {
	  for(k = 0; k <= (*adb->track_lengths)[track] - wL; k += HOP_SIZE) {
            double thisDist = 0;
            double qn = qpointers.l2norm[j];
            double sn = dbpointers.l2norm[trackIndexOffset + k];
            switch(spec->params.distance) {
            case ADB_DISTANCE_EUCLIDEAN_NORMED:
              thisDist = 2-(2/(qn*sn))*DD[j][k];
              break;
            case ADB_DISTANCE_EUCLIDEAN:
              thisDist = qn*qn + sn*sn - 2*DD[j][k];
              break;
            case ADB_DISTANCE_DOT_PRODUCT:
              thisDist = DD[j][k];
              break;
            }
	    // Power test
	    if ((!power_refine) || audiodb_powers_acceptable(&spec->refine, qpointers.power[j], dbpointers.power[trackIndexOffset + k])) {
              // radius test
              if((!(spec->refine.flags & ADB_REFINE_RADIUS)) || 
                 thisDist <= (spec->refine.radius+O2_DISTANCE_TOLERANCE)) {
                adb_result_t r;
                r.key = (*adb->keys)[track].c_str();
                r.dist = thisDist;
                if(spec->qid.flags & ADB_QUERY_ID_FLAG_EXHAUSTIVE) {
                  r.qpos = j;
                } else {
                  r.qpos = spec->qid.sequence_start;
                }
                r.ipos = k;
                accumulator->add_point(&r);
              }
            }
          }
        }
      } // Duration match            
      audiodb_delete_arrays(track, qpointers.nvectors, D, DD);
    }
  }

 loop_finish:

  free(data_buffer);

  // Clean up
  if(query_data)
    delete[] query_data;
  if(qpointers.l2norm_data)
    delete[] qpointers.l2norm_data;
  if(qpointers.power_data)
    delete[] qpointers.power_data;
  if(qpointers.mean_duration)
    delete[] qpointers.mean_duration;
  if(dbpointers.power_data)
    delete[] dbpointers.power_data;
  if(dbpointers.l2norm_data)
    delete[] dbpointers.l2norm_data;
  if(D)
    delete[] D;
  if(DD)
    delete[] DD;
  if(dbpointers.mean_duration)
    delete[] dbpointers.mean_duration;

  return 0;
}