mas01cr@239: #include "audioDB.h"
mas01cr@239: #include "reporter.h"
mas01cr@239: 
mas01cr@239: bool audioDB::powers_acceptable(double p1, double p2) {
mas01cr@239:   if (use_absolute_threshold) {
mas01cr@239:     if ((p1 < absolute_threshold) || (p2 < absolute_threshold)) {
mas01cr@239:       return false;
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239:   if (use_relative_threshold) {
mas01cr@239:     if (fabs(p1-p2) > fabs(relative_threshold)) {
mas01cr@239:       return false;
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239:   return true;
mas01cr@239: }
mas01cr@239: 
mas01cr@239: void audioDB::query(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse) {
mas01mc@292:   // init database tables and dbH first
mas01mc@292:   if(query_from_key)
mas01mc@292:     initTables(dbName);
mas01mc@292:   else
mas01mc@292:     initTables(dbName, inFile);
mas01mc@292: 
mas01mc@292:   // keyKeyPos requires dbH to be initialized
mas01mc@292:   if(query_from_key && (!key || (query_from_key_index = getKeyPos((char*)key))==O2_ERR_KEYNOTFOUND))
mas01mc@292:     error("Query key not found :",key);  
mas01mc@292:   
mas01cr@239:   switch (queryType) {
mas01cr@239:   case O2_POINT_QUERY:
mas01cr@239:     sequenceLength = 1;
mas01cr@239:     normalizedDistance = false;
mas01mc@292:     reporter = new pointQueryReporter< std::greater < NNresult > >(pointNN);
mas01cr@239:     break;
mas01cr@239:   case O2_TRACK_QUERY:
mas01cr@239:     sequenceLength = 1;
mas01cr@239:     normalizedDistance = false;
mas01mc@292:     reporter = new trackAveragingReporter< std::greater< NNresult > >(pointNN, trackNN, dbH->numFiles);
mas01cr@239:     break;
mas01mc@292:   case O2_SEQUENCE_QUERY:    
mas01mc@292:     if(no_unit_norming)
mas01mc@292:       normalizedDistance = false;
mas01cr@239:     if(radius == 0) {
mas01mc@292:       reporter = new trackAveragingReporter< std::less< NNresult > >(pointNN, trackNN, dbH->numFiles);
mas01cr@239:     } else {
mas01mc@292:       if(index_exists(dbName, radius, sequenceLength)){
mas01mc@292: 	char* indexName = index_get_name(dbName, radius, sequenceLength);
mas01mc@308: 	lsh = index_allocate(indexName, false);
mas01mc@319: 	reporter = new trackSequenceQueryRadReporter(trackNN, index_to_trackID(lsh->get_maxp(), lsh_n_point_bits)+1);
mas01mc@292: 	delete[] indexName;
mas01mc@292:       }
mas01mc@292:       else
mas01mc@292: 	reporter = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles);
mas01cr@239:     }
mas01cr@239:     break;
mas01mc@292:   case O2_N_SEQUENCE_QUERY:
mas01mc@292:     if(no_unit_norming)
mas01mc@292:       normalizedDistance = false;
mas01mc@248:     if(radius == 0) {
mas01mc@292:       reporter = new trackSequenceQueryNNReporter< std::less < NNresult > >(pointNN, trackNN, dbH->numFiles);
mas01mc@248:     } else {
mas01mc@292:       if(index_exists(dbName, radius, sequenceLength)){
mas01mc@292: 	char* indexName = index_get_name(dbName, radius, sequenceLength);
mas01mc@308: 	lsh = index_allocate(indexName, false);
mas01mc@319: 	reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, index_to_trackID(lsh->get_maxp(), lsh_n_point_bits)+1);
mas01mc@292: 	delete[] indexName;
mas01mc@292:       }
mas01mc@292:       else
mas01mc@292: 	reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, dbH->numFiles);
mas01mc@248:     }
mas01mc@248:     break;
mas01mc@263:   case O2_ONE_TO_ONE_N_SEQUENCE_QUERY :
mas01mc@263:     if(radius == 0) {
mas01mc@263:       error("query-type not yet supported");
mas01mc@263:     } else {
mas01mc@292:       reporter = new trackSequenceQueryRadNNReporterOneToOne(pointNN,trackNN, dbH->numFiles);
mas01mc@263:     }
mas01mc@263:     break;
mas01cr@239:   default:
mas01cr@239:     error("unrecognized queryType in query()");
mas01cr@239:   }  
mas01mc@292: 
mas01mc@292:   // Test for index (again) here
mas01mc@292:   if(radius && index_exists(dbName, radius, sequenceLength)) 
mas01mc@292:     index_query_loop(dbName, query_from_key_index);
mas01mc@292:   else
mas01mc@292:     query_loop(dbName, query_from_key_index);
mas01mc@292: 
mas01mc@292:   reporter->report(fileTable, adbQueryResponse);
mas01cr@239: }
mas01cr@239: 
mas01cr@239: // return ordinal position of key in keyTable
mas01mc@292: // this should really be a STL hash map search
mas01cr@239: unsigned audioDB::getKeyPos(char* key){  
mas01mc@292:   if(!dbH)
mas01mc@292:     error("dbH not initialized","getKeyPos");
mas01cr@239:   for(unsigned k=0; k<dbH->numFiles; k++)
mas01cr@256:     if(strncmp(fileTable + k*O2_FILETABLE_ENTRY_SIZE, key, strlen(key))==0)
mas01cr@239:       return k;
mas01cr@239:   error("Key not found",key);
mas01cr@239:   return O2_ERR_KEYNOTFOUND;
mas01cr@239: }
mas01cr@239: 
mas01cr@239: // This is a common pattern in sequence queries: what we are doing is
mas01cr@239: // taking a window of length seqlen over a buffer of length length,
mas01cr@239: // and placing the sum of the elements in that window in the first
mas01cr@239: // element of the window: thus replacing all but the last seqlen
mas01cr@239: // elements in the buffer with the corresponding windowed sum.
mas01cr@239: void audioDB::sequence_sum(double *buffer, int length, int seqlen) {
mas01cr@239:   double tmp1, tmp2, *ps;
mas01cr@239:   int j, w;
mas01cr@239: 
mas01cr@239:   tmp1 = *buffer;
mas01cr@239:   j = 1;
mas01cr@239:   w = seqlen - 1;
mas01cr@239:   while(w--) {
mas01cr@239:     *buffer += buffer[j++];
mas01cr@239:   }
mas01cr@239:   ps = buffer + 1;
mas01cr@239:   w = length - seqlen; // +1 - 1
mas01cr@239:   while(w--) {
mas01cr@239:     tmp2 = *ps;
mas01cr@239:     if(isfinite(tmp1)) {
mas01cr@239:       *ps = *(ps - 1) - tmp1 + *(ps + seqlen - 1);
mas01cr@239:     } else {
mas01cr@239:       for(int i = 1; i < seqlen; i++) {
mas01cr@239:         *ps += *(ps + i);
mas01cr@239:       }
mas01cr@239:     }
mas01cr@239:     tmp1 = tmp2;
mas01cr@239:     ps++;
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01cr@239: // In contrast to sequence_sum() above, sequence_sqrt() and
mas01cr@239: // sequence_average() below are simple mappers across the sequence.
mas01cr@239: void audioDB::sequence_sqrt(double *buffer, int length, int seqlen) {
mas01cr@239:   int w = length - seqlen + 1;
mas01cr@239:   while(w--) {
mas01cr@239:     *buffer = sqrt(*buffer);
mas01cr@239:     buffer++;
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01cr@239: void audioDB::sequence_average(double *buffer, int length, int seqlen) {
mas01cr@239:   int w = length - seqlen + 1;
mas01cr@239:   while(w--) {
mas01cr@239:     *buffer /= seqlen;
mas01cr@239:     buffer++;
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01cr@239: void audioDB::initialize_arrays(int track, unsigned int numVectors, double *query, double *data_buffer, double **D, double **DD) {
mas01cr@239:   unsigned int j, k, l, w;
mas01cr@239:   double *dp, *qp, *sp;
mas01cr@239: 
mas01cr@239:   const unsigned HOP_SIZE = sequenceHop;
mas01cr@239:   const unsigned wL = sequenceLength;
mas01cr@239: 
mas01cr@239:   for(j = 0; j < numVectors; j++) {
mas01cr@239:     // Sum products matrix
mas01cr@239:     D[j] = new double[trackTable[track]]; 
mas01cr@239:     assert(D[j]);
mas01cr@239:     // Matched filter matrix
mas01cr@239:     DD[j]=new double[trackTable[track]];
mas01cr@239:     assert(DD[j]);
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   // Dot product
mas01cr@239:   for(j = 0; j < numVectors; j++)
mas01cr@239:     for(k = 0; k < trackTable[track]; k++){
mas01cr@239:       qp = query + j * dbH->dim;
mas01cr@239:       sp = data_buffer + k * dbH->dim;
mas01cr@239:       DD[j][k] = 0.0; // Initialize matched filter array
mas01cr@239:       dp = &D[j][k];  // point to correlation cell j,k
mas01cr@239:       *dp = 0.0;      // initialize correlation cell
mas01cr@239:       l = dbH->dim;         // size of vectors
mas01cr@239:       while(l--)
mas01cr@239:         *dp += *qp++ * *sp++;
mas01cr@239:     }
mas01cr@239:   
mas01cr@239:   // Matched Filter
mas01cr@239:   // HOP SIZE == 1
mas01cr@239:   double* spd;
mas01cr@239:   if(HOP_SIZE == 1) { // HOP_SIZE = shingleHop
mas01cr@239:     for(w = 0; w < wL; w++) {
mas01cr@239:       for(j = 0; j < numVectors - w; j++) { 
mas01cr@239:         sp = DD[j];
mas01cr@239:         spd = D[j+w] + w;
mas01cr@239:         k = trackTable[track] - w;
mas01mc@292: 	while(k--)
mas01mc@292: 	  *sp++ += *spd++;
mas01cr@239:       }
mas01cr@239:     }
mas01cr@239:   } else { // HOP_SIZE != 1
mas01cr@239:     for(w = 0; w < wL; w++) {
mas01cr@239:       for(j = 0; j < numVectors - w; j += HOP_SIZE) {
mas01cr@239:         sp = DD[j];
mas01cr@239:         spd = D[j+w]+w;
mas01cr@239:         for(k = 0; k < trackTable[track] - w; k += HOP_SIZE) {
mas01cr@239:           *sp += *spd;
mas01cr@239:           sp += HOP_SIZE;
mas01cr@239:           spd += HOP_SIZE;
mas01cr@239:         }
mas01cr@239:       }
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01cr@239: void audioDB::delete_arrays(int track, unsigned int numVectors, double **D, double **DD) {
mas01cr@239:   if(D != NULL) {
mas01cr@239:     for(unsigned int j = 0; j < numVectors; j++) {
mas01cr@239:       delete[] D[j];
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239:   if(DD != NULL) {
mas01cr@239:     for(unsigned int j = 0; j < numVectors; j++) {
mas01cr@239:       delete[] DD[j];
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01mc@319: void audioDB::read_data(int trkfid, int track, double **data_buffer_p, size_t *data_buffer_size_p) {
mas01cr@239:   if (trackTable[track] * sizeof(double) * dbH->dim > *data_buffer_size_p) {
mas01cr@239:     if(*data_buffer_p) {
mas01cr@239:       free(*data_buffer_p);
mas01cr@239:     }
mas01cr@239:     { 
mas01cr@239:       *data_buffer_size_p = trackTable[track] * sizeof(double) * dbH->dim;
mas01cr@239:       void *tmp = malloc(*data_buffer_size_p);
mas01cr@239:       if (tmp == NULL) {
mas01cr@239:         error("error allocating data buffer");
mas01cr@239:       }
mas01cr@239:       *data_buffer_p = (double *) tmp;
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: 
mas01mc@319:   read(trkfid, *data_buffer_p, trackTable[track] * sizeof(double) * dbH->dim);
mas01cr@239: }
mas01cr@239: 
mas01cr@239: // These names deserve some unpicking.  The names starting with a "q"
mas01cr@239: // are pointers to the query, norm and power vectors; the names
mas01cr@239: // starting with "v" are things that will end up pointing to the
mas01cr@239: // actual query point's information.  -- CSR, 2007-12-05
mas01cr@239: void audioDB::set_up_query(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned *nvp) {
mas01cr@239:   *nvp = (statbuf.st_size - sizeof(int)) / (dbH->dim * sizeof(double));
mas01mc@292: 
mas01cr@239:   if(!(dbH->flags & O2_FLAG_L2NORM)) {
mas01cr@239:     error("Database must be L2 normed for sequence query","use -L2NORM");
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   if(*nvp < sequenceLength) {
mas01cr@239:     error("Query shorter than requested sequence length", "maybe use -l");
mas01cr@239:   }
mas01cr@239:   
mas01cr@239:   VERB_LOG(1, "performing norms... ");
mas01cr@239: 
mas01cr@239:   *qp = new double[*nvp * dbH->dim];
mas01cr@239:   memcpy(*qp, indata+sizeof(int), *nvp * dbH->dim * sizeof(double));
mas01cr@239:   *qnp = new double[*nvp];
mas01cr@239:   unitNorm(*qp, dbH->dim, *nvp, *qnp);
mas01cr@239: 
mas01cr@239:   sequence_sum(*qnp, *nvp, sequenceLength);
mas01cr@239:   sequence_sqrt(*qnp, *nvp, sequenceLength);
mas01cr@239: 
mas01cr@239:   if (usingPower) {
mas01cr@239:     *qpp = new double[*nvp];
mas01cr@239:     if (lseek(powerfd, sizeof(int), SEEK_SET) == (off_t) -1) {
mas01cr@239:       error("error seeking to data", powerFileName, "lseek");
mas01cr@239:     }
mas01cr@239:     int count = read(powerfd, *qpp, *nvp * sizeof(double));
mas01cr@239:     if (count == -1) {
mas01cr@239:       error("error reading data", powerFileName, "read");
mas01cr@239:     }
mas01cr@239:     if ((unsigned) count != *nvp * sizeof(double)) {
mas01cr@239:       error("short read", powerFileName);
mas01cr@239:     }
mas01cr@239: 
mas01cr@239:     sequence_sum(*qpp, *nvp, sequenceLength);
mas01cr@239:     sequence_average(*qpp, *nvp, sequenceLength);
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   if (usingTimes) {
mas01cr@239:     unsigned int k;
mas01cr@239:     *mqdp = 0.0;
mas01cr@239:     double *querydurs = new double[*nvp];
mas01cr@239:     double *timesdata = new double[*nvp*2];
mas01cr@239:     insertTimeStamps(*nvp, timesFile, timesdata);
mas01cr@239:     for(k = 0; k < *nvp; k++) {
mas01cr@239:       querydurs[k] = timesdata[2*k+1] - timesdata[2*k];
mas01cr@239:       *mqdp += querydurs[k];
mas01cr@239:     }
mas01cr@239:     *mqdp /= k;
mas01cr@239: 
mas01cr@239:     VERB_LOG(1, "mean query file duration: %f\n", *mqdp);
mas01cr@239: 
mas01cr@239:     delete [] querydurs;
mas01cr@239:     delete [] timesdata;
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   // Defaults, for exhaustive search (!usingQueryPoint)
mas01cr@239:   *vqp = *qp;
mas01cr@239:   *vqnp = *qnp;
mas01cr@239:   *vqpp = *qpp;
mas01cr@239: 
mas01cr@239:   if(usingQueryPoint) {
mas01cr@239:     if(queryPoint > *nvp || queryPoint > *nvp - sequenceLength + 1) {
mas01cr@239:       error("queryPoint > numVectors-wL+1 in query");
mas01cr@239:     } else {
mas01cr@239:       VERB_LOG(1, "query point: %u\n", queryPoint);
mas01cr@239:       *vqp = *qp + queryPoint * dbH->dim;
mas01cr@239:       *vqnp = *qnp + queryPoint;
mas01cr@239:       if (usingPower) {
mas01cr@239:         *vqpp = *qpp + queryPoint;
mas01cr@239:       }
mas01cr@239:       *nvp = sequenceLength;
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: }
mas01cr@239: 
mas01mc@292: // Does the same as set_up_query(...) but from database features instead of from a file
mas01mc@292: // Constructs the same outputs as set_up_query
mas01mc@292: void audioDB::set_up_query_from_key(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned *nvp, Uns32T queryIndex) {
mas01mc@292:   if(!trackTable)
mas01mc@292:     error("trackTable not initialized","set_up_query_from_key");
mas01mc@292: 
mas01mc@292:   if(!(dbH->flags & O2_FLAG_L2NORM)) {
mas01mc@292:     error("Database must be L2 normed for sequence query","use -L2NORM");
mas01mc@292:   }
mas01mc@292:   
mas01mc@292:   if(dbH->flags & O2_FLAG_POWER)
mas01mc@292:     usingPower = true;
mas01mc@292:   
mas01mc@292:   if(dbH->flags & O2_FLAG_TIMES)
mas01mc@292:     usingTimes = true;
mas01mc@292: 
mas01mc@292:   *nvp = trackTable[queryIndex];  
mas01mc@292:   if(*nvp < sequenceLength) {
mas01mc@292:     error("Query shorter than requested sequence length", "maybe use -l");
mas01mc@292:   }
mas01mc@292:   
mas01mc@292:   VERB_LOG(1, "performing norms... ");
mas01mc@292: 
mas01mc@320:   // For LARGE_ADB load query features from file
mas01mc@320:   if( dbH->flags & O2_FLAG_LARGE_ADB ){
mas01mc@320:     if(infid>0)
mas01mc@320:       close(infid);
mas01mc@321:     char* prefixedString = new char[O2_MAXFILESTR];
mas01mc@321:     char* tmpStr = prefixedString;
mas01mc@321:     strncpy(prefixedString, featureFileNameTable+queryIndex*O2_FILETABLE_ENTRY_SIZE, O2_MAXFILESTR);
mas01mc@321:     prefix_name(&prefixedString, adb_feature_root);
mas01mc@321:     if(tmpStr!=prefixedString)
mas01mc@321:       delete[] tmpStr;
mas01mc@321:     initInputFile(prefixedString, false); // nommap, file pointer at correct position
mas01mc@320:     size_t allocatedSize = 0;
mas01mc@320:     read_data(infid, queryIndex, qp, &allocatedSize); // over-writes qp and allocatedSize
mas01mc@320:     // Consistency check on allocated memory and query feature size
mas01mc@320:     if(*nvp*sizeof(double)*dbH->dim != allocatedSize)
mas01mc@320:       error("Query memory allocation failed consitency check","set_up_query_from_key");
mas01mc@320:     // Allocated and calculate auxillary sequences: l2norm and power
mas01mc@320:     init_track_aux_data(queryIndex, *qp, qnp, vqnp, qpp, vqpp);
mas01mc@320:   }
mas01mc@320:   else{ // Load from self-contained ADB database
mas01mc@320:     // Read query feature vectors from database
mas01mc@320:     *qp = NULL;
mas01mc@320:     lseek(dbfid, dbH->dataOffset + trackOffsetTable[queryIndex] * sizeof(double), SEEK_SET);
mas01mc@320:     size_t allocatedSize = 0;
mas01mc@320:     read_data(dbfid, queryIndex, qp, &allocatedSize);
mas01mc@320:     // Consistency check on allocated memory and query feature size
mas01mc@320:     if(*nvp*sizeof(double)*dbH->dim != allocatedSize)
mas01mc@320:       error("Query memory allocation failed consitency check","set_up_query_from_key");
mas01mc@320:     
mas01mc@320:     Uns32T trackIndexOffset = trackOffsetTable[queryIndex]/dbH->dim; // Convert num data elements to num vectors
mas01mc@320:     // Copy L2 norm partial-sum coefficients
mas01mc@320:     assert(*qnp = new double[*nvp]);
mas01mc@320:     memcpy(*qnp, l2normTable+trackIndexOffset, *nvp*sizeof(double));
mas01mc@320:     sequence_sum(*qnp, *nvp, sequenceLength);
mas01mc@320:     sequence_sqrt(*qnp, *nvp, sequenceLength);
mas01mc@320:     
mas01mc@320:     if( usingPower ){
mas01mc@320:       // Copy Power partial-sum coefficients
mas01mc@320:       assert(*qpp = new double[*nvp]);
mas01mc@320:       memcpy(*qpp, powerTable+trackIndexOffset, *nvp*sizeof(double));
mas01mc@320:       sequence_sum(*qpp, *nvp, sequenceLength);
mas01mc@320:       sequence_average(*qpp, *nvp, sequenceLength);
mas01mc@320:     }
mas01mc@320:     
mas01mc@320:     if (usingTimes) {
mas01mc@320:       unsigned int k;
mas01mc@320:       *mqdp = 0.0;
mas01mc@320:       double *querydurs = new double[*nvp];
mas01mc@320:       double *timesdata = new double[*nvp*2];
mas01mc@320:       assert(querydurs && timesdata);
mas01mc@320:       memcpy(timesdata, timesTable+trackIndexOffset, *nvp*sizeof(double));    
mas01mc@320:       for(k = 0; k < *nvp; k++) {
mas01mc@320: 	querydurs[k] = timesdata[2*k+1] - timesdata[2*k];
mas01mc@320: 	*mqdp += querydurs[k];
mas01mc@320:       }
mas01mc@320:       *mqdp /= k;
mas01mc@320:       
mas01mc@320:       VERB_LOG(1, "mean query file duration: %f\n", *mqdp);
mas01mc@320:       
mas01mc@320:       delete [] querydurs;
mas01mc@320:       delete [] timesdata;
mas01mc@320:     }
mas01mc@292:   }
mas01mc@292: 
mas01mc@292:   // Defaults, for exhaustive search (!usingQueryPoint)
mas01mc@292:   *vqp = *qp;
mas01mc@292:   *vqnp = *qnp;
mas01mc@292:   *vqpp = *qpp;
mas01mc@292: 
mas01mc@292:   if(usingQueryPoint) {
mas01mc@292:     if(queryPoint > *nvp || queryPoint > *nvp - sequenceLength + 1) {
mas01mc@292:       error("queryPoint > numVectors-wL+1 in query");
mas01mc@292:     } else {
mas01mc@292:       VERB_LOG(1, "query point: %u\n", queryPoint);
mas01mc@292:       *vqp = *qp + queryPoint * dbH->dim;
mas01mc@292:       *vqnp = *qnp + queryPoint;
mas01mc@292:       if (usingPower) {
mas01mc@292:         *vqpp = *qpp + queryPoint;
mas01mc@292:       }
mas01mc@292:       *nvp = sequenceLength;
mas01mc@292:     }
mas01mc@292:   }
mas01mc@292: }
mas01mc@292: 
mas01mc@292: 
mas01cr@239: // FIXME: this is not the right name; we're not actually setting up
mas01cr@239: // the database, but copying various bits of it out of mmap()ed tables
mas01cr@239: // in order to reduce seeks.
mas01cr@239: void audioDB::set_up_db(double **snp, double **vsnp, double **spp, double **vspp, double **mddp, unsigned int *dvp) {
mas01cr@239:   *dvp = dbH->length / (dbH->dim * sizeof(double));
mas01cr@239:   *snp = new double[*dvp];
mas01cr@239: 
mas01cr@239:   double *snpp = *snp, *sppp = 0;
mas01cr@239:   memcpy(*snp, l2normTable, *dvp * sizeof(double));
mas01cr@239: 
mas01cr@239:   if (usingPower) {
mas01cr@239:     if (!(dbH->flags & O2_FLAG_POWER)) {
mas01cr@239:       error("database not power-enabled", dbName);
mas01cr@239:     }
mas01cr@239:     *spp = new double[*dvp];
mas01cr@239:     sppp = *spp;
mas01cr@239:     memcpy(*spp, powerTable, *dvp * sizeof(double));
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   for(unsigned int i = 0; i < dbH->numFiles; i++){
mas01cr@239:     if(trackTable[i] >= sequenceLength) {
mas01cr@239:       sequence_sum(snpp, trackTable[i], sequenceLength);
mas01cr@239:       sequence_sqrt(snpp, trackTable[i], sequenceLength);
mas01cr@239: 
mas01cr@239:       if (usingPower) {
mas01cr@239: 	sequence_sum(sppp, trackTable[i], sequenceLength);
mas01cr@239:         sequence_average(sppp, trackTable[i], sequenceLength);
mas01cr@239:       }
mas01cr@239:     }
mas01cr@239:     snpp += trackTable[i];
mas01cr@239:     if (usingPower) {
mas01cr@239:       sppp += trackTable[i];
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   if (usingTimes) {
mas01cr@239:     if(!(dbH->flags & O2_FLAG_TIMES)) {
mas01cr@239:       error("query timestamps provided for non-timed database", dbName);
mas01cr@239:     }
mas01cr@239: 
mas01cr@239:     *mddp = new double[dbH->numFiles];
mas01cr@239: 
mas01cr@239:     for(unsigned int k = 0; k < dbH->numFiles; k++) {
mas01cr@239:       unsigned int j;
mas01cr@239:       (*mddp)[k] = 0.0;
mas01cr@239:       for(j = 0; j < trackTable[k]; j++) {
mas01cr@239: 	(*mddp)[k] += timesTable[2*j+1] - timesTable[2*j];
mas01cr@239:       }
mas01cr@239:       (*mddp)[k] /= j;
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   *vsnp = *snp;
mas01cr@239:   *vspp = *spp;
mas01cr@239: }
mas01cr@239: 
mas01mc@292: // query_points()
mas01mc@292: //
mas01mc@292: // using PointPairs held in the exact_evaluation_queue compute squared distance for each PointPair
mas01mc@292: // and insert result into the current reporter.
mas01mc@292: //
mas01mc@292: // Preconditions:
mas01mc@292: // A query inFile has been opened with setup_query(...) and query pointers initialized
mas01mc@292: // The database contains some points
mas01mc@292: // An exact_evaluation_queue has been allocated and populated
mas01mc@292: // A reporter has been allocated
mas01mc@292: //
mas01mc@292: // Postconditions:
mas01mc@292: // reporter contains the points and distances that meet the reporter constraints 
mas01mc@292: 
mas01mc@292: void audioDB::query_loop_points(double* query, double* qnPtr, double* qpPtr, double meanQdur, Uns32T numVectors){ 
mas01mc@292:   unsigned int dbVectors;
mas01mc@315:   double *sNorm = 0, *snPtr, *sPower = 0, *spPtr = 0;
mas01mc@292:   double *meanDBdur = 0;
mas01mc@292: 
mas01mc@292:   // check pre-conditions
mas01mc@292:   assert(exact_evaluation_queue&&reporter);
mas01mc@292:   if(!exact_evaluation_queue->size()) // Exit if no points to evaluate
mas01mc@292:     return;
mas01mc@292: 
mas01mc@292:   // Compute database info
mas01mc@292:   // FIXME: we more than likely don't need very much of the database
mas01mc@292:   // so make a new method to build these values per-track or, even better, per-point
mas01mc@320:   if( !( dbH->flags & O2_FLAG_LARGE_ADB) )
mas01mc@320:     set_up_db(&sNorm, &snPtr, &sPower, &spPtr, &meanDBdur, &dbVectors);
mas01mc@292: 
mas01mc@292:   VERB_LOG(1, "matching points...");
mas01mc@292: 
mas01mc@292:   assert(pointNN>0 && pointNN<=O2_MAXNN);
mas01mc@292:   assert(trackNN>0 && trackNN<=O2_MAXNN);
mas01mc@292: 
mas01mc@292:   // We are guaranteed that the order of points is sorted by:
mas01mc@320:   // trackID, spos, qpos
mas01mc@292:   // so we can be relatively efficient in initialization of track data.
mas01mc@292:   // Here we assume that points don't overlap, so we will use exhaustive dot
mas01mc@320:   // product evaluation instead of memoization of partial sums which is used
mas01mc@320:   // for exhaustive brute-force evaluation from smaller databases: e.g. query_loop()
mas01mc@292:   double dist;
mas01mc@292:   size_t data_buffer_size = 0;
mas01mc@292:   double *data_buffer = 0;
mas01mc@320:   Uns32T trackOffset = 0;
mas01mc@320:   Uns32T trackIndexOffset = 0;
mas01mc@292:   Uns32T currentTrack = 0x80000000; // Initialize with a value outside of track index range
mas01mc@292:   Uns32T npairs = exact_evaluation_queue->size();
mas01mc@292:   while(npairs--){
mas01mc@292:     PointPair pp = exact_evaluation_queue->top();
mas01mc@320:     // Large ADB track data must be loaded here for sPower
mas01mc@320:     if(dbH->flags & O2_FLAG_LARGE_ADB){
mas01mc@320:       trackOffset=0;
mas01mc@320:       trackIndexOffset=0;
mas01mc@292:       if(currentTrack!=pp.trackID){
mas01mc@321: 	char* prefixedString = new char[O2_MAXFILESTR];
mas01mc@321: 	char* tmpStr = prefixedString;
mas01mc@320: 	// On currentTrack change, allocate and load track data
mas01mc@292: 	currentTrack=pp.trackID;
mas01mc@320: 	SAFE_DELETE_ARRAY(sNorm);
mas01mc@320: 	SAFE_DELETE_ARRAY(sPower);
mas01mc@320: 	if(infid>0)
mas01mc@320: 	  close(infid);
mas01mc@320: 	// Open and check dimensions of feature file
mas01mc@321: 	strncpy(prefixedString, featureFileNameTable+pp.trackID*O2_FILETABLE_ENTRY_SIZE, O2_MAXFILESTR);
mas01mc@321: 	prefix_name((char ** const) &prefixedString, adb_feature_root);
mas01mc@321: 	if (prefixedString!=tmpStr)
mas01mc@321: 	  delete[] tmpStr;
mas01mc@321: 	initInputFile(prefixedString, false); // nommap, file pointer at correct position
mas01mc@320: 	// Load the feature vector data for current track into data_buffer
mas01mc@320: 	read_data(infid, pp.trackID, &data_buffer, &data_buffer_size);	
mas01mc@320: 	// Load power and calculate power and l2norm sequence sums
mas01mc@320: 	init_track_aux_data(pp.trackID, data_buffer, &sNorm, &snPtr, &sPower, &spPtr);
mas01mc@320:       }
mas01mc@320:     }
mas01mc@320:     else{
mas01mc@320:       // These offsets are w.r.t. the entire database of feature vectors and auxillary variables
mas01mc@320:       trackOffset=trackOffsetTable[pp.trackID]; // num data elements offset
mas01mc@320:       trackIndexOffset=trackOffset/dbH->dim;    // num vectors offset
mas01mc@320:     }    
mas01mc@320:     Uns32T qPos = usingQueryPoint?0:pp.qpos;// index for query point
mas01mc@320:     Uns32T sPos = trackIndexOffset+pp.spos; // index into l2norm table
mas01mc@320:     // Test power thresholds before computing distance
mas01mc@320:     if( ( !usingPower || powers_acceptable(qpPtr[qPos], sPower[sPos])) &&
mas01mc@320: 	( qPos<numVectors-sequenceLength+1 && pp.spos<trackTable[pp.trackID]-sequenceLength+1 ) ){
mas01mc@320:       // Non-large ADB track data is loaded inside power test for efficiency
mas01mc@320:       if( !(dbH->flags & O2_FLAG_LARGE_ADB) && (currentTrack!=pp.trackID) ){
mas01mc@320: 	// On currentTrack change, allocate and load track data
mas01mc@320: 	currentTrack=pp.trackID;
mas01mc@320: 	lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
mas01mc@319: 	read_data(dbfid, currentTrack, &data_buffer, &data_buffer_size);
mas01mc@292:       }
mas01mc@320:       // Compute distance
mas01mc@320:       dist = dot_product_points(query+qPos*dbH->dim, data_buffer+pp.spos*dbH->dim, dbH->dim*sequenceLength);
mas01mc@320:       double qn = qnPtr[qPos];
mas01mc@320:       double sn = sNorm[sPos];
mas01mc@292:       if(normalizedDistance) 
mas01mc@320: 	dist = 2 - (2/(qn*sn))*dist;
mas01mc@292:       else 
mas01mc@292: 	if(no_unit_norming)
mas01mc@320: 	  dist = qn*qn + sn*sn - 2*dist;
mas01mc@292:       // else
mas01mc@292:       // dist = dist;      
mas01mc@314:       if((!radius) || dist <= (O2_LSH_EXACT_MULT*radius+O2_DISTANCE_TOLERANCE)) 
mas01mc@320: 	reporter->add_point(pp.trackID, pp.qpos, pp.spos, dist);    
mas01mc@292:     }
mas01mc@292:     exact_evaluation_queue->pop();
mas01mc@292:   }
mas01mc@315:   // Cleanup
mas01mc@320:   SAFE_DELETE_ARRAY(sNorm);
mas01mc@320:   SAFE_DELETE_ARRAY(sPower);
mas01mc@320:   SAFE_DELETE_ARRAY(meanDBdur);
mas01mc@292: }
mas01mc@292: 
mas01mc@292: // A completely unprotected dot-product method
mas01mc@292: // Caller is responsible for ensuring that memory is within bounds
mas01mc@292: inline double audioDB::dot_product_points(double* q, double* p, Uns32T  L){
mas01mc@292:   double dist = 0.0;
mas01mc@292:   while(L--)
mas01mc@292:     dist += *q++ * *p++;
mas01mc@292:   return dist;
mas01mc@292: }
mas01mc@292: 
mas01mc@292: void audioDB::query_loop(const char* dbName, Uns32T queryIndex) {
mas01cr@239:   
mas01cr@239:   unsigned int numVectors;
mas01cr@239:   double *query, *query_data;
mas01cr@239:   double *qNorm, *qnPtr, *qPower = 0, *qpPtr = 0;
mas01cr@239:   double meanQdur;
mas01cr@239: 
mas01mc@319:   if( dbH->flags & O2_FLAG_LARGE_ADB )
mas01mc@319:     error("error: LARGE_ADB requires indexed query");
mas01mc@319: 
mas01mc@292:   if(query_from_key)
mas01mc@292:     set_up_query_from_key(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors, queryIndex);
mas01mc@292:   else
mas01mc@292:     set_up_query(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors);
mas01cr@239: 
mas01cr@239:   unsigned int dbVectors;
mas01cr@239:   double *sNorm, *snPtr, *sPower = 0, *spPtr = 0;
mas01cr@239:   double *meanDBdur = 0;
mas01cr@239: 
mas01cr@239:   set_up_db(&sNorm, &snPtr, &sPower, &spPtr, &meanDBdur, &dbVectors);
mas01cr@239: 
mas01cr@239:   VERB_LOG(1, "matching tracks...");
mas01cr@239:   
mas01cr@239:   assert(pointNN>0 && pointNN<=O2_MAXNN);
mas01cr@239:   assert(trackNN>0 && trackNN<=O2_MAXNN);
mas01cr@239: 
mas01cr@239:   unsigned j,k,track,trackOffset=0, HOP_SIZE=sequenceHop, wL=sequenceLength;
mas01cr@239:   double **D = 0;    // Differences query and target 
mas01cr@239:   double **DD = 0;   // Matched filter distance
mas01cr@239: 
mas01mc@292:   D = new double*[numVectors]; // pre-allocate 
mas01cr@239:   DD = new double*[numVectors];
mas01cr@239: 
mas01cr@239:   gettimeofday(&tv1, NULL); 
mas01cr@239:   unsigned processedTracks = 0;
mas01cr@239:   off_t trackIndexOffset;
mas01cr@239:   char nextKey[MAXSTR];
mas01cr@239: 
mas01cr@239:   // Track loop 
mas01cr@239:   size_t data_buffer_size = 0;
mas01cr@239:   double *data_buffer = 0;
mas01cr@239:   lseek(dbfid, dbH->dataOffset, SEEK_SET);
mas01cr@239: 
mas01cr@239:   for(processedTracks=0, track=0 ; processedTracks < dbH->numFiles ; track++, processedTracks++) {
mas01cr@239: 
mas01cr@239:     trackOffset = trackOffsetTable[track];     // numDoubles offset
mas01cr@239: 
mas01cr@239:     // get trackID from file if using a control file
mas01cr@239:     if(trackFile) {
mas01cr@239:       trackFile->getline(nextKey,MAXSTR);
mas01cr@239:       if(!trackFile->eof()) {
mas01cr@239: 	track = getKeyPos(nextKey);
mas01cr@239:         trackOffset = trackOffsetTable[track];
mas01cr@239:         lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
mas01cr@239:       } else {
mas01cr@239: 	break;
mas01cr@239:       }
mas01cr@239:     }
mas01cr@239: 
mas01mc@292:     // skip identity on query_from_key
mas01mc@292:     if( query_from_key && (track == queryIndex) ) {
mas01mc@292:       if(queryIndex!=dbH->numFiles-1){
mas01mc@292: 	track++;
mas01mc@292: 	trackOffset = trackOffsetTable[track];
mas01mc@292: 	lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET);
mas01mc@292:       }
mas01mc@292:       else{
mas01mc@292: 	break;
mas01mc@292:       }
mas01mc@292:     }
mas01mc@292: 
mas01cr@239:     trackIndexOffset=trackOffset/dbH->dim; // numVectors offset
mas01cr@239: 
mas01mc@319:     read_data(dbfid, track, &data_buffer, &data_buffer_size);
mas01cr@239:     if(sequenceLength <= trackTable[track]) {  // test for short sequences
mas01cr@239:       
mas01cr@239:       VERB_LOG(7,"%u.%jd.%u | ", track, (intmax_t) trackIndexOffset, trackTable[track]);
mas01cr@239:       
mas01cr@239:       initialize_arrays(track, numVectors, query, data_buffer, D, DD);
mas01cr@239: 
mas01cr@239:       if(usingTimes) {
mas01cr@239:         VERB_LOG(3,"meanQdur=%f meanDBdur=%f\n", meanQdur, meanDBdur[track]);
mas01cr@239:       }
mas01cr@239: 
mas01cr@239:       if((!usingTimes) || fabs(meanDBdur[track]-meanQdur) < meanQdur*timesTol) {
mas01cr@239:         if(usingTimes) {
mas01cr@239:           VERB_LOG(3,"within duration tolerance.\n");
mas01cr@239:         }
mas01cr@239: 
mas01cr@239: 	// Search for minimum distance by shingles (concatenated vectors)
mas01cr@239: 	for(j = 0; j <= numVectors - wL; j += HOP_SIZE) {
mas01cr@239: 	  for(k = 0; k <= trackTable[track] - wL; k += HOP_SIZE) {
mas01cr@239:             double thisDist;
mas01mc@292:             if(normalizedDistance) 
mas01cr@239:               thisDist = 2-(2/(qnPtr[j]*sNorm[trackIndexOffset+k]))*DD[j][k];
mas01mc@292: 	    else 
mas01mc@292: 	      if(no_unit_norming)
mas01mc@292: 		thisDist = qnPtr[j]*qnPtr[j]+sNorm[trackIndexOffset+k]*sNorm[trackIndexOffset+k] - 2*DD[j][k];
mas01mc@292: 	      else
mas01mc@292: 		thisDist = DD[j][k];
mas01mc@292: 
mas01cr@239: 	    // Power test
mas01cr@239: 	    if ((!usingPower) || powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k])) {
mas01cr@239:               // radius test
mas01mc@292:               if((!radius) || thisDist <= (radius+O2_DISTANCE_TOLERANCE)) {
mas01cr@239:                 reporter->add_point(track, usingQueryPoint ? queryPoint : j, k, thisDist);
mas01cr@239:               }
mas01cr@239:             }
mas01cr@239:           }
mas01cr@239:         }
mas01cr@239:       } // Duration match            
mas01cr@239:       delete_arrays(track, numVectors, D, DD);
mas01cr@239:     }
mas01cr@239:   }
mas01cr@239: 
mas01cr@239:   free(data_buffer);
mas01cr@239: 
mas01cr@239:   gettimeofday(&tv2,NULL);
mas01cr@239:   VERB_LOG(1,"elapsed time: %ld msec\n",
mas01cr@239:            (tv2.tv_sec*1000 + tv2.tv_usec/1000) - 
mas01cr@239:            (tv1.tv_sec*1000 + tv1.tv_usec/1000))
mas01cr@239: 
mas01cr@239:   // Clean up
mas01cr@239:   if(query_data)
mas01cr@239:     delete[] query_data;
mas01cr@239:   if(qNorm)
mas01cr@239:     delete[] qNorm;
mas01cr@239:   if(sNorm)
mas01cr@239:     delete[] sNorm;
mas01cr@239:   if(qPower)
mas01cr@239:     delete[] qPower;
mas01cr@239:   if(sPower)
mas01cr@239:     delete[] sPower;
mas01cr@239:   if(D)
mas01cr@239:     delete[] D;
mas01cr@239:   if(DD)
mas01cr@239:     delete[] DD;
mas01cr@239:   if(meanDBdur)
mas01cr@239:     delete[] meanDBdur;
mas01cr@239: }
mas01cr@239: 
mas01cr@239: // Unit norm block of features
mas01cr@239: void audioDB::unitNorm(double* X, unsigned dim, unsigned n, double* qNorm){
mas01cr@239:   unsigned d;
mas01cr@239:   double L2, *p;
mas01cr@239: 
mas01cr@239:   VERB_LOG(2, "norming %u vectors...", n);
mas01cr@239:   while(n--) {
mas01cr@239:     p = X;
mas01cr@239:     L2 = 0.0;
mas01cr@239:     d = dim;
mas01cr@239:     while(d--) {
mas01cr@239:       L2 += *p * *p;
mas01cr@239:       p++;
mas01cr@239:     }
mas01cr@239:     if(qNorm) {
mas01cr@239:       *qNorm++=L2;
mas01cr@239:     }
mas01cr@239:     X += dim;
mas01cr@239:   }
mas01cr@239:   VERB_LOG(2, "done.\n");
mas01cr@239: }
mas01mc@292: 
mas01mc@292: