annotate query.cpp @ 443:cb44e57a96fa api-inversion

New (internal) function audiodb_query_spec_qpointers()

Intended to replace audioDB::set_up_query and
audioDB::set_up_query_from_key; it isn't yet called, but I want to save
my work now before messing with stuff, because the function itself is
long and complicated: it handles three separate ways of getting data
together (datum provided, key provided + LARGE_ADB, key provided +
!LARGE_ADB).

There is insufficient error checking as yet, but hey.

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

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

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

  qspec.refine.flags = 0;
  /* FIXME: trackFile / ADB_REFINE_KEYLIST */
  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;
  }

  /* FIXME qspec.qid.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(&qspec, dbName, query_from_key_index);
  }
  else{
    VERB_LOG(1, "Calling brute-force query on database %s\n", dbName);
    query_loop(&qspec, query_from_key_index);
  }

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

static 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) {
    std::map<std::string,uint32_t>::iterator it;
    it = adb->keys->find(datum->key);
    if(it == adb->keys->end()) {
      return 1;
    }
    uint32_t track_id = (*it).second;
    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;
  if(sequence_start > nvectors - sequence_length) {
    /* is there something to free?  goto error */
    return 1;
  }

  qpointers->nvectors = nvectors;

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

  qpointers->l2norm_data = (double *) malloc(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 = (double *) malloc(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 = (double *) calloc(1, sizeof(double));
    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;
    qpointers->power = qpointers->power_data + spec->qid.sequence_start;
  }


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

// These names deserve some unpicking.  The names starting with a "q"
// are pointers to the query, norm and power vectors; the names
// starting with "v" are things that will end up pointing to the
// actual query point's information.  -- CSR, 2007-12-05
void audioDB::set_up_query(adb_query_spec_t *spec, double **qp, double **vqp, adb_qpointers_internal_t *qpointers) {
  uint32_t nvectors = (statbuf.st_size - sizeof(int)) / (dbH->dim * sizeof(double));
  qpointers->nvectors = nvectors;

  uint32_t sequence_length = spec->qid.sequence_length;
  if(!(dbH->flags & O2_FLAG_L2NORM)) {
    error("Database must be L2 normed for sequence query","use -L2NORM");
  }

  if(nvectors < sequence_length) {
    error("Query shorter than requested sequence length", "maybe use -l");
  }
  
  VERB_LOG(1, "performing norms... ");

  *qp = new double[nvectors * dbH->dim];
  memcpy(*qp, indata+sizeof(int), nvectors * dbH->dim * sizeof(double));
  qpointers->l2norm_data = new double[nvectors];
  audiodb_l2norm_buffer(*qp, dbH->dim, nvectors, qpointers->l2norm_data);

  audiodb_sequence_sum(qpointers->l2norm_data, nvectors, sequence_length);
  audiodb_sequence_sqrt(qpointers->l2norm_data, nvectors, sequence_length);

  if (usingPower) {
    qpointers->power_data = new double[nvectors];
    if (lseek(powerfd, sizeof(int), SEEK_SET) == (off_t) -1) {
      error("error seeking to data", powerFileName, "lseek");
    }
    int count = read(powerfd, qpointers->power_data, nvectors * sizeof(double));
    if (count == -1) {
      error("error reading data", powerFileName, "read");
    }
    if ((unsigned) count != nvectors * sizeof(double)) {
      error("short read", powerFileName);
    }

    audiodb_sequence_sum(qpointers->power_data, nvectors, sequence_length);
    audiodb_sequence_average(qpointers->power_data, nvectors, sequence_length);
  }

  if (usingTimes) {
    unsigned int k;
    qpointers->mean_duration = new double[1];
    *qpointers->mean_duration = 0.0;
    double *querydurs = new double[nvectors];
    double *timesdata = new double[2*nvectors];
    insertTimeStamps(nvectors, timesFile, timesdata);
    for(k = 0; k < nvectors; k++) {
      querydurs[k] = timesdata[2*k+1] - timesdata[2*k];
      *qpointers->mean_duration += querydurs[k];
    }
    *qpointers->mean_duration /= k;

    VERB_LOG(1, "mean query file duration: %f\n", *qpointers->mean_duration);

    delete [] querydurs;
    delete [] timesdata;
  }

  // Defaults, for exhaustive search (!usingQueryPoint)
  *vqp = *qp;
  qpointers->l2norm = qpointers->l2norm_data;
  qpointers->power = qpointers->power_data;

  if(usingQueryPoint) {
    if( !(queryPoint < nvectors && queryPoint < nvectors - sequence_length + 1) ) {
      error("queryPoint >= nvectors-sequence_length+1 in query");
    } else {
      VERB_LOG(1, "query point: %u\n", queryPoint);
      *vqp = *qp + queryPoint * dbH->dim;
      qpointers->l2norm = qpointers->l2norm_data + queryPoint;
      if (usingPower) {
        qpointers->power = qpointers->power_data + queryPoint;
      }
      qpointers->nvectors = sequence_length;
    }
  }
}

// Does the same as set_up_query(...) but from database features instead of from a file
// Constructs the same outputs as set_up_query
void audioDB::set_up_query_from_key(adb_query_spec_t *spec, double **qp, double **vqp, adb_qpointers_internal_t *qpointers, Uns32T queryIndex) {
  uint32_t sequence_length = spec->qid.sequence_length;
  if(!trackTable)
    error("trackTable not initialized","set_up_query_from_key");

  if(!(dbH->flags & O2_FLAG_L2NORM)) {
    error("Database must be L2 normed for sequence query","use -L2NORM");
  }
  
  if(dbH->flags & O2_FLAG_POWER)
    usingPower = true;
  
  if(dbH->flags & O2_FLAG_TIMES)
    usingTimes = true;

  uint32_t nvectors = trackTable[queryIndex];
  qpointers->nvectors = nvectors;
  if(nvectors < sequence_length) {
    error("Query shorter than requested sequence length", "maybe use -l");
  }
  
  VERB_LOG(1, "performing norms... ");

  // For LARGE_ADB load query features from file
  if( dbH->flags & O2_FLAG_LARGE_ADB ){
    if(infid>0)
      close(infid);
    char* prefixedString = new char[O2_MAXFILESTR];
    char* tmpStr = prefixedString;
    strncpy(prefixedString, featureFileNameTable+queryIndex*O2_FILETABLE_ENTRY_SIZE, O2_MAXFILESTR);
    prefix_name(&prefixedString, adb_feature_root);
    if(tmpStr!=prefixedString)
      delete[] tmpStr;
    initInputFile(prefixedString, false); // nommap, file pointer at correct position
    size_t allocatedSize = 0;
    if(audiodb_read_data(adb, infid, queryIndex, qp, &allocatedSize))
      error("failed to read data"); // over-writes qp and allocatedSize
    // Consistency check on allocated memory and query feature size
    if(nvectors*sizeof(double)*dbH->dim != allocatedSize)
      error("Query memory allocation failed consitency check","set_up_query_from_key");
    // Allocated and calculate auxillary sequences: l2norm and power
    init_track_aux_data(queryIndex, *qp, &qpointers->l2norm_data, &qpointers->l2norm, &qpointers->power_data, &qpointers->power);
  }
  else{ // Load from self-contained ADB database
    // Read query feature vectors from database
    *qp = NULL;
    lseek(dbfid, dbH->dataOffset + trackOffsetTable[queryIndex] * sizeof(double), SEEK_SET);
    size_t allocatedSize = 0;
    if(audiodb_read_data(adb, dbfid, queryIndex, qp, &allocatedSize))
      error("failed to read data");
    // Consistency check on allocated memory and query feature size
    if(nvectors*sizeof(double)*dbH->dim != allocatedSize)
      error("Query memory allocation failed consitency check","set_up_query_from_key");
    
    Uns32T trackIndexOffset = trackOffsetTable[queryIndex]/dbH->dim; // Convert num data elements to num vectors
    // Copy L2 norm partial-sum coefficients
    assert(qpointers->l2norm_data = new double[nvectors]);
    memcpy(qpointers->l2norm_data, l2normTable+trackIndexOffset, nvectors*sizeof(double));
    audiodb_sequence_sum(qpointers->l2norm_data, nvectors, sequence_length);
    audiodb_sequence_sqrt(qpointers->l2norm_data, nvectors, sequence_length);
    
    if( usingPower ){
      // Copy Power partial-sum coefficients
      assert(qpointers->power_data = new double[nvectors]);
      memcpy(qpointers->power_data, powerTable+trackIndexOffset, nvectors*sizeof(double));
      audiodb_sequence_sum(qpointers->power_data, nvectors, sequence_length);
      audiodb_sequence_average(qpointers->power_data, nvectors, sequence_length);
    }
    
    if (usingTimes) {
      unsigned int k;
      qpointers->mean_duration = new double[1];
      *qpointers->mean_duration = 0.0;
      double *querydurs = new double[nvectors];
      double *timesdata = new double[nvectors*2];
      assert(querydurs && timesdata);
      memcpy(timesdata, timesTable+trackIndexOffset, nvectors*sizeof(double));    
      for(k = 0; k < nvectors; k++) {
	querydurs[k] = timesdata[2*k+1] - timesdata[2*k];
	*qpointers->mean_duration += querydurs[k];
      }
      *qpointers->mean_duration /= k;
      
      VERB_LOG(1, "mean query file duration: %f\n", *qpointers->mean_duration);
      
      delete [] querydurs;
      delete [] timesdata;
    }
  }

  // Defaults, for exhaustive search (!usingQueryPoint)
  *vqp = *qp;
  qpointers->l2norm = qpointers->l2norm_data;
  qpointers->power = qpointers->power_data;

  if(usingQueryPoint) {
    if( !(queryPoint < nvectors && queryPoint < nvectors - sequence_length + 1) ) {
      error("queryPoint >= nvectors-sequence_length+1 in query");
    } else {
      VERB_LOG(1, "query point: %u\n", queryPoint);
      *vqp = *qp + queryPoint * dbH->dim;
      qpointers->l2norm = qpointers->l2norm_data + queryPoint;
      if (usingPower) {
        qpointers->power = qpointers->power_data + queryPoint;
      }
      qpointers->nvectors = sequence_length;
    }
  }
}


// FIXME: this is not the right name; we're not actually setting up
// the database, but copying various bits of it out of mmap()ed tables
// in order to reduce seeks.
static int audiodb_set_up_db(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_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);

  // check pre-conditions
  assert(exact_evaluation_queue&&reporter);
  if(!exact_evaluation_queue->size()) // Exit if no points to evaluate
    return;

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

  VERB_LOG(1, "matching points...");

  // 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.
  // Here we assume that points 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: e.g. 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(dbH->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, false); // nommap, file pointer at correct position
	// 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=trackOffsetTable[pp.trackID]; // num data elements offset
      trackIndexOffset=trackOffset/dbH->dim;    // 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<trackTable[pp.trackID]-sequence_length+1 ) ){
      // Non-large ADB track data is loaded inside power test for efficiency
      if( !(dbH->flags & O2_FLAG_LARGE_ADB) && (currentTrack!=pp.trackID) ){
	// On currentTrack change, allocate and load track data
	currentTrack=pp.trackID;
	lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), 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*dbH->dim, data_buffer+pp.spos*dbH->dim, dbH->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 = fileTable + pp.trackID * O2_FILETABLE_ENTRY_SIZE;
        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);
}

void audioDB::query_loop(adb_query_spec_t *spec, Uns32T queryIndex) {
  
  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( dbH->flags & O2_FLAG_LARGE_ADB )
    error("error: LARGE_ADB requires indexed query");

  if(query_from_key)
    set_up_query_from_key(spec, &query_data, &query, &qpointers, queryIndex);
  else
    set_up_query(spec, &query_data, &query, &qpointers);

  if(audiodb_set_up_db(adb, spec, &dbpointers)) {
    error("failed to set up db");
  }

  VERB_LOG(1, "matching tracks...");
  
  unsigned j,k,track,trackOffset=0, HOP_SIZE=sequenceHop;
  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];

  gettimeofday(&tv1, NULL); 
  unsigned processedTracks = 0;
  off_t trackIndexOffset;
  char nextKey[MAXSTR];

  // 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 = audiodb_key_index(adb, nextKey);
        if(track == (uint32_t) -1) {
          error("key not found", nextKey);
        }
        trackOffset = trackOffsetTable[track];
        lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
      } else {
	break;
      }
    }

    // skip identity on query_from_key
    if( query_from_key && (track == queryIndex) ) {
      if(queryIndex!=dbH->numFiles-1){
	track++;
	trackOffset = trackOffsetTable[track];
	lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
      }
      else{
	break;
      }
    }

    trackIndexOffset=trackOffset/dbH->dim; // qpointers.nvectors offset

    if(audiodb_read_data(adb, dbfid, track, &data_buffer, &data_buffer_size))
      error("failed to read data");
    if(wL <= trackTable[track]) {  // test for short sequences
      
      VERB_LOG(7,"%u.%jd.%u | ", track, (intmax_t) trackIndexOffset, trackTable[track]);
      
      audiodb_initialize_arrays(adb, spec, track, qpointers.nvectors, query, data_buffer, D, DD);

      if(spec->refine.flags & ADB_REFINE_DURATION_RATIO) {
        VERB_LOG(3,"meanQdur=%f meanDBdur=%f\n", qpointers.mean_duration[0], dbpointers.mean_duration[track]);
      }

      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) {
        if(spec->refine.flags & ADB_REFINE_DURATION_RATIO) {
          VERB_LOG(3,"within duration tolerance.\n");
        }

	// Search for minimum distance by shingles (concatenated vectors)
	for(j = 0; j <= qpointers.nvectors - wL; j += HOP_SIZE) {
	  for(k = 0; k <= trackTable[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 = fileTable + track * O2_FILETABLE_ENTRY_SIZE;
                r.dist = thisDist;
                r.qpos = usingQueryPoint ? queryPoint : j;
                r.ipos = k;
                accumulator->add_point(&r);
              }
            }
          }
        }
      } // Duration match            
      audiodb_delete_arrays(track, qpointers.nvectors, D, DD);
    }
  }

  free(data_buffer);

  gettimeofday(&tv2,NULL);
  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(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;
}