mas01cr@239: #include "audioDB.h" mas01cr@239: 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) { mas01cr@239: initTables(dbName, inFile); mas01cr@239: Reporter *r = 0; mas01cr@239: switch (queryType) { mas01cr@239: case O2_POINT_QUERY: mas01cr@239: sequenceLength = 1; mas01cr@239: normalizedDistance = false; mas01cr@239: r = new pointQueryReporter >(pointNN); mas01cr@239: break; mas01cr@239: case O2_TRACK_QUERY: mas01cr@239: sequenceLength = 1; mas01cr@239: normalizedDistance = false; mas01cr@239: r = new trackAveragingReporter >(pointNN, trackNN, dbH->numFiles); mas01cr@239: break; mas01cr@239: case O2_SEQUENCE_QUERY: mas01cr@239: if(radius == 0) { mas01cr@239: r = new trackAveragingReporter >(pointNN, trackNN, dbH->numFiles); mas01cr@239: } else { mas01cr@239: r = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles); mas01cr@239: } mas01cr@239: break; mas01mc@248: case O2_N_SEQUENCE_QUERY : mas01mc@248: if(radius == 0) { mas01mc@248: r = new trackSequenceQueryNNReporter >(pointNN, trackNN, dbH->numFiles); mas01mc@248: } else { mas01mc@250: r = 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@263: r = new trackSequenceQueryRadNNReporterOneToOne(pointNN,trackNN, dbH->numFiles); mas01mc@263: } mas01mc@263: break; mas01cr@239: default: mas01cr@239: error("unrecognized queryType in query()"); mas01cr@239: } mas01cr@240: query_loop(dbName, inFile, r); mas01cr@239: r->report(fileTable, adbQueryResponse); mas01cr@239: delete r; mas01cr@239: } mas01cr@239: mas01cr@239: // return ordinal position of key in keyTable mas01cr@239: unsigned audioDB::getKeyPos(char* key){ mas01cr@239: for(unsigned k=0; knumFiles; 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; mas01cr@239: while(k--) mas01cr@239: *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: mas01cr@239: void audioDB::read_data(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: mas01cr@239: read(dbfid, *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)); mas01cr@239: 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: 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: mas01cr@240: void audioDB::query_loop(const char* dbName, const char* inFile, Reporter *reporter) { 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: mas01cr@239: 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: mas01cr@239: D = new double*[numVectors]; mas01cr@239: DD = new double*[numVectors]; mas01cr@239: mas01cr@239: gettimeofday(&tv1, NULL); mas01cr@239: unsigned processedTracks = 0; mas01cr@239: mas01cr@239: // build track offset table mas01cr@239: off_t *trackOffsetTable = new off_t[dbH->numFiles]; mas01cr@239: unsigned cumTrack=0; mas01cr@239: off_t trackIndexOffset; mas01cr@239: for(k = 0; k < dbH->numFiles; k++){ mas01cr@239: trackOffsetTable[k] = cumTrack; mas01cr@239: cumTrack += trackTable[k] * dbH->dim; mas01cr@239: } mas01cr@239: 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: mas01cr@239: trackIndexOffset=trackOffset/dbH->dim; // numVectors offset mas01cr@239: mas01cr@239: read_data(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; mas01cr@239: if(normalizedDistance) { mas01cr@239: thisDist = 2-(2/(qnPtr[j]*sNorm[trackIndexOffset+k]))*DD[j][k]; mas01cr@239: } else { mas01cr@239: thisDist = DD[j][k]; mas01cr@239: } mas01cr@239: // Power test mas01cr@239: if ((!usingPower) || powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k])) { mas01cr@239: // radius test mas01cr@239: if((!radius) || thisDist < radius) { 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(trackOffsetTable) mas01cr@239: delete[] trackOffsetTable; 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: }