view query.cpp @ 462:f689510baaf4 api-inversion

Simplify audioDB::query_loop_points.

Using the new functions audiodb_track_id_datum() and
audiodb_datum_qpointers(), much of the body of the method disappears.
Of course, we've probably introduced some inefficiencies and extra
memory copies, but I'm fairly sure that this method is going to be
dominated by disk i/o time anyway, so it doesn't matter.

#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 (audiodb_index_exists(adb->path, qspec.refine.radius, qspec.qid.sequence_length)) {
	char *indexName = audiodb_index_get_name(adb->path, qspec.refine.radius, qspec.qid.sequence_length);
	lsh = index_allocate(indexName, false);
	reporter = new trackSequenceQueryRadReporter(trackNN, audiodb_index_to_track_id(lsh->get_maxp(), audiodb_lsh_n_point_bits(adb))+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 (audiodb_index_exists(adb->path, qspec.refine.radius, qspec.qid.sequence_length)){
	char *indexName = audiodb_index_get_name(adb->path, qspec.refine.radius, qspec.qid.sequence_length);
	lsh = index_allocate(indexName, false);
	reporter = new trackSequenceQueryRadNNReporter(pointNN, trackNN, audiodb_index_to_track_id(lsh->get_maxp(), audiodb_lsh_n_point_bits(adb))+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");
  }

  /* Somewhere around here is where the implementation of
   * audiodb_query_spec() starts. */

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

  switch(qspec.params.distance) {
  case ADB_DISTANCE_DOT_PRODUCT:
    switch(qspec.params.accumulation) {
    case ADB_ACCUMULATION_DB:
      qstate.accumulator = new DBAccumulator<adb_result_dist_gt>(qspec.params.npoints);
      break;
    case ADB_ACCUMULATION_PER_TRACK:
      qstate.accumulator = new PerTrackAccumulator<adb_result_dist_gt>(qspec.params.npoints, qspec.params.ntracks);
      break;
    case ADB_ACCUMULATION_ONE_TO_ONE:
      qstate.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:
      qstate.accumulator = new DBAccumulator<adb_result_dist_lt>(qspec.params.npoints);
      break;
    case ADB_ACCUMULATION_PER_TRACK:
      qstate.accumulator = new PerTrackAccumulator<adb_result_dist_lt>(qspec.params.npoints, qspec.params.ntracks);
      break;
    case ADB_ACCUMULATION_ONE_TO_ONE:
      qstate.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) && audiodb_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, &qstate);
  }
  else{
    VERB_LOG(1, "Calling brute-force query on database %s\n", dbName);
    if(query_loop(adb, &qspec, &qstate)) {
      error("query_loop failed");
    }
  }

  adb_query_results_t *rs = qstate.accumulator->get_points();

  delete qstate.accumulator;
  delete qstate.allowed_keys;

  /* End of audiodb_query_spec() function */

  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);
  }
  audiodb_query_free_results(adb, &qspec, rs);

  reporter->report(fileTable, adbQueryResponse);
}

int audiodb_query_free_results(adb_t *adb, adb_query_spec_t *spec, adb_query_results_t *rs) {
  free(rs->results);
  free(rs);
  return 0;
}

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_track_id_datum(adb_t *adb, uint32_t track_id, adb_datum_t *d) {
  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 = (*adb->keys)[track_id].c_str();
    /* 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));
    }
  }
  return 0;
}

int audiodb_datum_qpointers(adb_datum_t *d, uint32_t sequence_length, double **vector_data, double **vector, adb_qpointers_internal_t *qpointers) {
  uint32_t nvectors = d->nvectors;

  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;
  }

  *vector = *vector_data;
  qpointers->l2norm = qpointers->l2norm_data;
  qpointers->power = qpointers->power_data;
  return 0;
}

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 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;
    }
    audiodb_track_id_datum(adb, track_id, &d);
  } else {
    return 1;
  }

  /* FIXME: check the overflow logic here */
  if(sequence_start + sequence_length > d.nvectors) {
    if(datum != &d) {
      audiodb_free_datum(&d);
    }
    return 1;
  }

  audiodb_datum_qpointers(&d, sequence_length, vector_data, vector, qpointers);

  /* Finally, if applicable, set up the moving qpointers. */
  if(spec->qid.flags & ADB_QUERY_ID_FLAG_EXHAUSTIVE) {
    /* the qpointers are already at the start, and so correct. */
  } else {
    /* adjust the qpointers to point to the correct place in the sequence */
    *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;
    }
    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, adb_qstate_internal_t *qstate, 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(qstate->exact_evaluation_queue->size() == 0) {
    return;
  }

  /* 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, as in query_loop()). */
  double dist;
  double *dbdata = 0, *dbdata_pointer;
  Uns32T currentTrack = 0x80000000; // KLUDGE: Initialize with a value outside of track index range
  Uns32T npairs = qstate->exact_evaluation_queue->size();
  while(npairs--) {
    PointPair pp = qstate->exact_evaluation_queue->top();
    if(currentTrack != pp.trackID) {
      SAFE_DELETE_ARRAY(dbdata);
      SAFE_DELETE_ARRAY(dbpointers.l2norm_data);
      SAFE_DELETE_ARRAY(dbpointers.power_data);
      SAFE_DELETE_ARRAY(dbpointers.mean_duration);
      currentTrack = pp.trackID;
      adb_datum_t d = {0};
      if(audiodb_track_id_datum(adb, pp.trackID, &d)) {
        error("failed to get datum");
      }
      if(audiodb_datum_qpointers(&d, sequence_length, &dbdata, &dbdata_pointer, &dbpointers)) {
        audiodb_free_datum(&d);
        error("failed to get dbpointers");
      }
      audiodb_free_datum(&d);
    }
    Uns32T qPos = usingQueryPoint?0:pp.qpos;// index for query point
    Uns32T sPos = 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 && sPos<(*adb->track_lengths)[pp.trackID]-sequence_length+1 ) ){
      // Compute distance
      dist = audiodb_dot_product(query + qPos*adb->header->dim, dbdata + 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;
        qstate->accumulator->add_point(&r);
      }
    }
    qstate->exact_evaluation_queue->pop();
  }
  // Cleanup
  SAFE_DELETE_ARRAY(dbdata);
  SAFE_DELETE_ARRAY(dbpointers.l2norm_data);
  SAFE_DELETE_ARRAY(dbpointers.power_data);
  SAFE_DELETE_ARRAY(dbpointers.mean_duration);
  delete qstate->exact_evaluation_queue;
}

int audioDB::query_loop(adb_t *adb, adb_query_spec_t *spec, adb_qstate_internal_t *qstate) {
  
  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);

  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);

  std::set<std::string>::iterator keys_end = qstate->allowed_keys->end();
  for(track = 0; track < adb->header->numFiles; track++) {
    unsigned t = track;
    
    while (qstate->allowed_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;
                qstate->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;
}