Mercurial > hg > audiodb
changeset 292:d9a88cfd4ab6
Completed merge of lshlib back to current version of the trunk.
| author | mas01mc |
|---|---|
| date | Tue, 29 Jul 2008 22:01:17 +0000 |
| parents | 63ae0dfc1767 |
| children | 9fd5340faffd |
| files | Makefile MakefileLSH ReporterBase.h UNIT_TEST_LSH.cpp audioDB.cpp audioDB.h common.cpp gengetopt.in index.cpp lshlib.cpp lshlib.h mt19937/mt19937ar.c mt19937/mt19937ar.h mt19937/mt19937ar.out mt19937/mtTest.c mt19937/readme-mt.txt query.cpp reporter.h |
| diffstat | 18 files changed, 4090 insertions(+), 318 deletions(-) [+] |
line wrap: on
line diff
--- a/Makefile Tue Jul 22 20:09:31 2008 +0000 +++ b/Makefile Tue Jul 29 22:01:17 2008 +0000 @@ -36,10 +36,10 @@ soapServer.cpp soapClient.cpp soapC.cpp adb.nsmap: audioDBws.h ${SOAPCPP2} audioDBws.h -%.o: %.cpp audioDB.h adb.nsmap cmdline.h reporter.h +%.o: %.cpp audioDB.h adb.nsmap cmdline.h reporter.h ReporterBase.h lshlib.h g++ -c ${CFLAGS} ${GSOAP_INCLUDE} ${GSL_INCLUDE} -Wall $< -OBJS=insert.o create.o common.o dump.o query.o soap.o sample.o audioDB.o +OBJS=insert.o create.o common.o dump.o query.o soap.o sample.o audioDB.o index.o lshlib.o ${EXECUTABLE}: ${OBJS} soapServer.cpp soapClient.cpp soapC.cpp cmdline.c g++ -o ${EXECUTABLE} ${CFLAGS} ${GSL_INCLUDE} ${LIBGSL} ${GSOAP_INCLUDE} $^ ${GSOAP_CPP}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MakefileLSH Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,34 @@ +CFLAGS= -O3 -g +#CFLAGS= -ggdb --save-temps +MTSOURCES=mt19937/mt19937ar.c +CFLAGS+=-DMT19937 +GENGETOPT=gengetopt +SOAPCPP2=soapcpp2 +GSOAP_CPP=-lgsoap++ +GSOAP_INCLUDE= + +TARGET=TEST_LSH_LIB +TARGET2=UNIT_TEST_LSH +LSHOBS=TEST_LSH_LIB.cpp lshlib.cpp +LSHOBS2=UNIT_TEST_LSH.cpp lshlib.cpp +LSHHDRS=TEST_LSH_LIB.h lshlib.h reporter.h audioDB.h ReporterBase.h +LSHHDRS2=lshlib.h + +ifeq ($(shell uname),Linux) +override CFLAGS+=-D_LARGEFILE64_SOURCE -D_FILE_OFFSET_BITS=64 +endif + +all: ${TARGET} + +unit: ${TARGET2} + +${TARGET}: ${LSHOBS} ${LSHHDRS} ${MTSOURCES} MakefileLSH + ${GENGETOPT} -e <gengetopt.in + ${SOAPCPP2} audioDBws.h + g++ -o ${TARGET} ${CFLAGS} ${LSHOBS} ${MTSOURCES} + +${TARGET2}: ${LSHOBS2} ${LSHHDRS2} ${MTSOURCES} MakefileLSH + ${GENGETOPT} -e <gengetopt.in + ${SOAPCPP2} audioDBws.h + g++ -o ${TARGET2} ${CFLAGS} ${LSHOBS2} ${MTSOURCES} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ReporterBase.h Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,12 @@ + +#ifndef __REPORTERBASE_H +#define __REPORTERBASE_H + +class ReporterBase { +public: + virtual ~ReporterBase(){}; + virtual void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) = 0; + virtual void report(char*,void*) = 0; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UNIT_TEST_LSH.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,140 @@ +// UNIT_TEST_LSH.cpp + +#include <vector> +#include "lshlib.h" +#include "reporter.h" + +#define LSH_IN_CORE + + +#define N_POINT_BITS 14 +#define POINT_BIT_MASK 0x00003FFF + +// Callback method for LSH point retrieval +void add_point(void* reporter, Uns32T pointID, Uns32T qpos, float dist) +{ + ReporterBase* pr = (ReporterBase*)reporter; + pr->add_point(pointID>>N_POINT_BITS, qpos, pointID&POINT_BIT_MASK, dist); +} + +int main(int argc, char* argv[]){ + + int nT = 100; // num tracks + int nP = 1000; // num points-per-track + float w = 4.0;// LSH bucket width + int k = 10; + int m = 2; + int d = 10; + int N = 100000; + int C = 200; + + float radius = 0.001; + char FILENAME[] = "foo.lsh"; + + assert(nP>=nT); + + int fid = open(FILENAME,O_RDONLY); + LSH* lsh; + bool serialized = false; + Uns32T trackBase = 0; + + if(fid< 0){ // Make a new serial LSH file + lsh = new LSH(w,k,m,d,N,C,radius); + assert(lsh); + cout << "NEW LSH:" << endl; + } + else{ + close(fid); // Load LSH structures from disk + lsh = new LSH(FILENAME); + assert(lsh); + cout << "MERGE WITH EXISTING LSH:" << FILENAME << endl; + serialized=true; + trackBase = (lsh->get_maxp()>>N_POINT_BITS)+1; // Our encoding of tracks and points + } + cout << "k:" << lsh->k << " "; + cout << "m:" << lsh->m << "(L:" << lsh->L << ") "; + cout << "d:" << lsh->d << " "; + cout << "N:" << lsh->N << " "; + cout << "C:" << lsh->C << " "; + cout << "R:" << lsh->get_radius() << endl; + cout << "p:" << lsh->p << endl; + cout.flush(); + + cout << endl << "Constructing " << nT << " tracks with " << nP << " vectors of dimension " << d << endl; + cout.flush(); + // Construct sets of database vectors, use one point from each set for testing + vector< vector<float> > vv = vector< vector<float> >(nP); // track vectors + vector< vector<float> > qq = vector< vector<float> >(nP);// query vectors + for(int i=0; i< nP ; i++){ + vv[i]=vector<float>(d); // allocate vector + qq[i]=vector<float>(d); // allocate vector + } + + for(int k = 0 ; k < nT ; k ++){ + cout << "[" << k << "]"; + cout.flush(); + for(int i = 0 ; i< nP ; i++) + for(int j=0; j< d ; j++) + vv[i][j] = genrand_real2() / radius; // MT_19937 random numbers + lsh->insert_point_set(vv, (trackBase+k)<<N_POINT_BITS); + qq[k] = vv[k]; // One identity query per set of database vectors + } + cout << endl; + cout.flush(); + + cout << "Writing serialized LSH tables..." << endl; + // TEST SERIALIZED LSH RETRIEVAL + lsh->serialize(FILENAME); + + // TEST LSH RETRIEVAL IN CORE + printf("\n********** In-core LSH retrieval from %d track%c **********\n", + (lsh->get_maxp()>>N_POINT_BITS)+1,(lsh->get_maxp()>>N_POINT_BITS)>0?'s':' '); + fflush(stdout); + for(int i = 0; i < nT ; i++ ){ + trackSequenceQueryRadNNReporter* pr = new trackSequenceQueryRadNNReporter(nP,nT,(lsh->get_maxp()>>N_POINT_BITS)+1); + lsh->retrieve_point(qq[i], i, &add_point, (void*)pr); // LSH point retrieval from core + printf("query vector %d] t1:%u t2:%0X\n", i, lsh->get_t1(), lsh->get_t2()); + fflush(stdout); + pr->report(0,0); + delete pr; + } + delete lsh; + + cout << "Loading Serialized LSH functions from disk ..." << endl; + cout.flush(); + lsh = new LSH(FILENAME); + assert(lsh); + // lsh->serial_dump_tables(FILENAME); + printf("\n********** Serialized LSH retrieval from %d track%c **********\n", (lsh->get_maxp()>>N_POINT_BITS)+1,(lsh->get_maxp()>>N_POINT_BITS)>1?'s':' '); + fflush(stdout); + for(int i= 0; i < nT ; i++ ){ + trackSequenceQueryRadNNReporter* pr = new trackSequenceQueryRadNNReporter(nP,nT,(lsh->get_maxp()>>N_POINT_BITS)+1); + lsh->serial_retrieve_point(FILENAME, qq[i], i, &add_point, (void*) pr); // LSH serialized point retrieval method + printf("query vector %d] t1:%u t2:%0X\n", i, lsh->get_t1(), lsh->get_t2()); + fflush(stdout); + pr->report(0,0); + delete pr; + } + delete lsh; + +#ifdef LSH_IN_CORE + cout << "Loading Serialized LSH functions and tables from disk ..." << endl; + cout.flush(); + // Unserialize entire lsh tree to core + lsh = new LSH(FILENAME,1); + + // TEST UNSERIALIZED LSH RETRIEVAL IN CORE + printf("\n********** Unserialized LSH in-core retrieval from %d track%c **********\n", (lsh->get_maxp()>>N_POINT_BITS)+1,(lsh->get_maxp()>>N_POINT_BITS)>1?'s':' '); + fflush(stdout); + for(int i = 0; i < nT ; i++ ){ + trackSequenceQueryRadNNReporter* pr = new trackSequenceQueryRadNNReporter(nP,nT,(lsh->get_maxp()>>N_POINT_BITS)+1); + lsh->retrieve_point(qq[i], i, &add_point, (void*) pr); // LSH point retrieval from core + printf("query vector %d] t1:%u t2:%0X\n", i, lsh->get_t1(), lsh->get_t2()); + fflush(stdout); + pr->report(0,0); + delete pr; + } + delete lsh; +#endif + +}
--- a/audioDB.cpp Tue Jul 22 20:09:31 2008 +0000 +++ b/audioDB.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -1,5 +1,23 @@ #include "audioDB.h" +PointPair::PointPair(Uns32T a, Uns32T b, Uns32T c):trackID(a),qpos(b),spos(c){}; + +bool operator<(const PointPair& a, const PointPair& b){ + return ( (a.qpos<b.qpos) || + ((a.qpos==b.qpos) && + ( (a.trackID<b.trackID)) || ((a.trackID==b.trackID)&&(a.spos<b.spos)) ) ); +} + +bool operator>(const PointPair& a, const PointPair& b){ + return ( (a.qpos>b.qpos) || + ((a.qpos==b.qpos) && + ( (a.trackID>b.trackID)) || ((a.trackID==b.trackID)&&(a.spos>b.spos)) ) ); +} + +bool operator==(const PointPair& a, const PointPair& b){ + return ( (a.trackID==b.trackID) && (a.qpos==b.qpos) && (a.spos==b.spos) ); +} + audioDB::audioDB(const unsigned argc, char* const argv[]): O2_AUDIODB_INITIALIZERS { if(processArgs(argc, argv)<0){ @@ -49,6 +67,9 @@ else if(O2_ACTION(COM_DUMP)) dump(dbName); + + else if(O2_ACTION(COM_INDEX)) + index_index_db(dbName); else error("Unrecognized command",command); @@ -96,10 +117,16 @@ munmap(timesTable, timesTableLength); if(l2normTable) munmap(l2normTable, l2normTableLength); - + if(trackOffsetTable) + delete trackOffsetTable; + if(reporter) + delete reporter; + if(exact_evaluation_queue) + delete exact_evaluation_queue; if(rng) gsl_rng_free(rng); - + if(vv) + delete vv; if(dbfid>0) close(dbfid); if(infid>0) @@ -179,6 +206,27 @@ } } + sequenceLength = args_info.sequencelength_arg; + if(sequenceLength < 1 || sequenceLength > 1000) { + error("seqlen out of range: 1 <= seqlen <= 1000"); + } + sequenceHop = args_info.sequencehop_arg; + if(sequenceHop < 1 || sequenceHop > 1000) { + error("seqhop out of range: 1 <= seqhop <= 1000"); + } + + if (args_info.absolute_threshold_given) { + if (args_info.absolute_threshold_arg >= 0) { + error("absolute threshold out of range: should be negative"); + } + use_absolute_threshold = true; + absolute_threshold = args_info.absolute_threshold_arg; + } + if (args_info.relative_threshold_given) { + use_relative_threshold = true; + relative_threshold = args_info.relative_threshold_arg; + } + if(args_info.SERVER_given){ command=COM_SERVER; port=args_info.SERVER_arg; @@ -251,7 +299,10 @@ dbName=args_info.database_arg; inFile=args_info.features_arg; if(args_info.key_given) - key=args_info.key_arg; + if(!args_info.features_given) + error("INSERT: '-k key' argument depends on '-f features'"); + else + key=args_info.key_arg; if(args_info.times_given){ timesFileName=args_info.times_arg; if(strlen(timesFileName)>0){ @@ -277,7 +328,10 @@ dbName=args_info.database_arg; inFile=args_info.featureList_arg; if(args_info.keyList_given) - key=args_info.keyList_arg; // INCONSISTENT NO CHECK + if(!args_info.features_given) + error("INSERT: '-k key' argument depends on '-f features'"); + else + key=args_info.key_arg; // INCONSISTENT NO CHECK /* TO DO: REPLACE WITH if(args_info.keyList_given){ @@ -306,13 +360,64 @@ } return 0; } - + + // Set no_unit_norm flag + no_unit_norming = args_info.no_unit_norming_flag; + lsh_use_u_functions = args_info.lsh_use_u_functions_flag; + + // LSH Index Command + if(args_info.INDEX_given){ + if(radius <= 0 ) + error("INDEXing requires a Radius argument"); + if(!(sequenceLength>0 && sequenceLength <= O2_MAXSEQLEN)) + error("INDEXing requires 1 <= sequenceLength <= 1000"); + command=COM_INDEX; + dbName=args_info.database_arg; + + // Whether to store LSH hash tables for query in core (FORMAT2) + lsh_in_core = args_info.lsh_inCore_flag; + + lsh_param_w = args_info.lsh_w_arg; + if(!(lsh_param_w>0 && lsh_param_w<=O2_SERIAL_MAX_BINWIDTH)) + error("Indexing parameter w out of range (0.0 < w <= 100.0)"); + + lsh_param_k = args_info.lsh_k_arg; + if(!(lsh_param_k>0 && lsh_param_k<=O2_SERIAL_MAX_FUNS)) + error("Indexing parameter k out of range (1 <= k <= 100)"); + + lsh_param_m = args_info.lsh_m_arg; + if(!(lsh_param_m>0 && lsh_param_m<= (1 + (sqrt(1 + O2_SERIAL_MAX_TABLES*8.0)))/2.0)) + error("Indexing parameter m out of range (1 <= m <= 20)"); + + lsh_param_N = args_info.lsh_N_arg; + if(!(lsh_param_N>0 && lsh_param_N<=O2_SERIAL_MAX_ROWS)) + error("Indexing parameter N out of range (1 <= N <= 1000000)"); + + lsh_param_b = args_info.lsh_b_arg; + if(!(lsh_param_b>0 && lsh_param_b<=O2_SERIAL_MAX_TRACKBATCH)) + error("Indexing parameter b out of range (1 <= b <= 10000)"); + + lsh_param_ncols = args_info.lsh_ncols_arg; + if( !(lsh_param_ncols>0 && lsh_param_ncols<=O2_SERIAL_MAX_COLS)) + error("Indexing parameter ncols out of range (1 <= ncols <= 1000"); + + return 0; + } + // Query command and arguments if(args_info.QUERY_given){ command=COM_QUERY; dbName=args_info.database_arg; - inFile=args_info.features_arg; - + // XOR features and key search + if(!args_info.features_given && !args_info.key_given || (args_info.features_given && args_info.key_given)) + error("QUERY requires exactly one of either -f features or -k key"); + if(args_info.features_given) + inFile=args_info.features_arg; // query from file + else{ + query_from_key = true; + key=args_info.key_arg; // query from key + } + if(args_info.keyList_given){ trackFileName=args_info.keyList_arg; if(strlen(trackFileName)>0 && !(trackFile = new std::ifstream(trackFileName,std::ios::in))) @@ -358,7 +463,13 @@ if(queryPoint<0 || queryPoint >10000) error("queryPoint out of range: 0 <= queryPoint <= 10000"); } - + + // Whether to pre-load LSH hash tables for query + lsh_in_core = args_info.lsh_inCore_flag; + + // Whether to perform exact evaluation of points returned by LSH + lsh_exact = args_info.lsh_exact_flag; + pointNN = args_info.pointnn_arg; if(pointNN < 1 || pointNN > O2_MAXNN) { error("pointNN out of range: 1 <= pointNN <= 1000000"); @@ -367,26 +478,6 @@ if(trackNN < 1 || trackNN > O2_MAXNN) { error("resultlength out of range: 1 <= resultlength <= 1000000"); } - sequenceLength = args_info.sequencelength_arg; - if(sequenceLength < 1 || sequenceLength > 1000) { - error("seqlen out of range: 1 <= seqlen <= 1000"); - } - sequenceHop = args_info.sequencehop_arg; - if(sequenceHop < 1 || sequenceHop > 1000) { - error("seqhop out of range: 1 <= seqhop <= 1000"); - } - - if (args_info.absolute_threshold_given) { - if (args_info.absolute_threshold_arg >= 0) { - error("absolute threshold out of range: should be negative"); - } - use_absolute_threshold = true; - absolute_threshold = args_info.absolute_threshold_arg; - } - if (args_info.relative_threshold_given) { - use_relative_threshold = true; - relative_threshold = args_info.relative_threshold_arg; - } return 0; } return -1; // no command found
--- a/audioDB.h Tue Jul 22 20:09:31 2008 +0000 +++ b/audioDB.h Tue Jul 29 22:01:17 2008 +0000 @@ -1,3 +1,6 @@ +#ifndef __AUDIODB_H_ +#define __AUDIODB_H_ + #include <stdio.h> #include <stdlib.h> #include <sys/types.h> @@ -14,6 +17,10 @@ #include <signal.h> #include <gsl/gsl_rng.h> +// includes for LSH indexing +#include "ReporterBase.h" +#include "lshlib.h" + // includes for web services #include "soapH.h" #include "cmdline.h" @@ -30,6 +37,7 @@ #define COM_POWER "--POWER" #define COM_DUMP "--DUMP" #define COM_SERVER "--SERVER" +#define COM_INDEX "--INDEX" #define COM_SAMPLE "--SAMPLE" // parameters @@ -62,15 +70,21 @@ #define O2_DEFAULT_DATASIZE (1355U) // in MB #define O2_DEFAULT_NTRACKS (20000U) #define O2_DEFAULT_DATADIM (9U) - +#define O2_REALTYPE (double) #define O2_MAXFILES (20000U) #define O2_MAXFILESTR (256U) #define O2_FILETABLE_ENTRY_SIZE (O2_MAXFILESTR) #define O2_TRACKTABLE_ENTRY_SIZE (sizeof(unsigned)) #define O2_HEADERSIZE (sizeof(dbTableHeaderT)) #define O2_MEANNUMVECTORS (1000U) -#define O2_MAXDIM (1000U) +#define O2_MAXDIM (2000U) #define O2_MAXNN (1000000U) +#define O2_MAXSEQLEN (8000U) // maximum feature vectors in a sequence +#define O2_MAXTRACKS (10000U) // maximum number of tracks +#define O2_MAXTRACKLEN (1000000U) // maximum shingles in a track +#define O2_MAXDOTPRODUCTMEMORY (sizeof(O2_REALTYPE)*O2_MAXSEQLEN*O2_MAXSEQLEN) // 512MB +#define O2_DISTANCE_TOLERANCE (1e-6) +#define O2_SERIAL_MAX_TRACKBATCH (10000) // Flags #define O2_FLAG_L2NORM (0x1U) @@ -85,6 +99,10 @@ #define O2_N_SEQUENCE_QUERY (0x20U) #define O2_ONE_TO_ONE_N_SEQUENCE_QUERY (0x40U) +// Bit masks for packing (trackID,pointID) into 32-bit unsigned int +#define LSH_N_POINT_BITS 14 +#define LSH_TRACK_MASK 0xFFFFC000 // 2^18 = 262144 tracks +#define LSH_POINT_MASK 0x00003FFF // 2^14 = 16384 points per track // Error Codes #define O2_ERR_KEYNOTFOUND (0xFFFFFF00) @@ -131,7 +149,16 @@ off_t dbSize; } dbTableHeaderT, *dbTableHeaderPtr; -class Reporter; + +class PointPair{ + public: + Uns32T trackID; + Uns32T qpos; + Uns32T spos; + PointPair(Uns32T a, Uns32T b, Uns32T c); +}; + +bool operator<(const PointPair& a, const PointPair& b); class audioDB{ @@ -153,6 +180,7 @@ int powerfd; int dbfid; + int lshfid; bool forWrite; int infid; char* db; @@ -164,6 +192,7 @@ char *fileTable; unsigned* trackTable; + off_t *trackOffsetTable; double* dataBuf; double* inBuf; double* l2normTable; @@ -193,6 +222,7 @@ unsigned sequenceLength; unsigned sequenceHop; bool normalizedDistance; + bool no_unit_norming; unsigned queryPoint; unsigned usingQueryPoint; unsigned usingTimes; @@ -202,13 +232,16 @@ unsigned port; double timesTol; double radius; - + bool query_from_key; + Uns32T query_from_key_index; bool use_absolute_threshold; double absolute_threshold; bool use_relative_threshold; double relative_threshold; - + ReporterBase* reporter; // track/point reporter + priority_queue<PointPair, std::vector<PointPair>, std::less<PointPair> >* exact_evaluation_queue; + // Timers struct timeval tv1; struct timeval tv2; @@ -223,20 +256,22 @@ void delete_arrays(int track, unsigned int numVectors, double **D, double **DD); void read_data(int track, double **data_buffer_p, size_t *data_buffer_size_p); void set_up_query(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned int *nvp); + void set_up_query_from_key(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned *nvp, Uns32T queryIndex); void set_up_db(double **snp, double **vsnp, double **spp, double **vspp, double **mddp, unsigned int *dvp); - void query_loop(const char* dbName, const char* inFile, Reporter *reporter); - + void query_loop(const char* dbName, Uns32T queryIndex); + void query_loop_points(double* query, double* qnPtr, double* qpPtr, double meanQdur, Uns32T numVectors); + double dot_product_points(double* q, double* p, Uns32T L); void initRNG(); void initDBHeader(const char *dbName); void initInputFile(const char *inFile); - void initTables(const char* dbName, const char* inFile); + void initTables(const char* dbName, const char* inFile = 0); + void initTablesFromKey(const char* dbName, const Uns32T queryIndex); void unitNorm(double* X, unsigned d, unsigned n, double* qNorm); void unitNormAndInsertL2(double* X, unsigned dim, unsigned n, unsigned append); void insertTimeStamps(unsigned n, std::ifstream* timesFile, double* timesdata); void insertPowerData(unsigned n, int powerfd, double *powerdata); unsigned getKeyPos(char* key); public: - audioDB(const unsigned argc, char* const argv[]); audioDB(const unsigned argc, char* const argv[], adb__queryResponse *adbQueryResponse); audioDB(const unsigned argc, char* const argv[], adb__statusResponse *adbStatusResponse); @@ -263,64 +298,119 @@ bool powers_acceptable(double p1, double p2); void dump(const char* dbName); - // web services + // LSH indexing parameters and data structures + LSH* lsh; + bool lsh_in_core; // load LSH tables for query into core (true) or keep on disk (false) + bool lsh_use_u_functions; + bool lsh_exact; // flag to indicate use exact evaluation of points returned by LSH + double lsh_param_w; // Width of LSH hash-function bins + Uns32T lsh_param_k; // Number of independent hash functions + Uns32T lsh_param_m; // Combinatorial parameter for m(m-1)/2 hash tables + Uns32T lsh_param_N; // Number of rows per hash table + Uns32T lsh_param_b; // Batch size, in number of tracks, per indexing iteration + Uns32T lsh_param_ncols; // Maximum number of collision in a hash-table row + + // LSH vector<> containers for one in-core copy of a set of feature vectors + vector<float>::iterator vi; // feature vector iterator + vector<vector<float> > *vv; // one-track's worth data + + // LSH indexing and retrieval methods + void index_index_db(const char* dbName); + void index_initialize(double**,double**,double**,double**,unsigned int*); + void index_insert_tracks(Uns32T start_track, Uns32T end_track, double** fvpp, double** sNormpp,double** snPtrp, double** sPowerp, double** spPtrp); + int index_insert_track(Uns32T trackID, double** fvpp, double** snpp, double** sppp); + Uns32T index_insert_shingles(vector<vector<float> >*, Uns32T trackID, double* spp); + void index_make_shingle(vector<vector<float> >*, Uns32T idx, double* fvp, Uns32T dim, Uns32T seqLen); + int index_norm_shingles(vector<vector<float> >*, double* snp, double* spp); + int index_query_loop(const char* dbName, Uns32T queryIndex); + vector<vector<float> >* index_initialize_shingles(Uns32T sz); + int index_init_query(const char* dbName); + int index_exists(const char* dbName, double radius, Uns32T sequenceLength); + char* index_get_name(const char*dbName, double radius, Uns32T sequenceLength); + static void index_add_point_approximate(void* instance, Uns32T pointID, Uns32T qpos, float dist); // static point reporter callback method + static void index_add_point_exact(void* instance, Uns32T pointID, Uns32T qpos, float dist); // static point reporter callback method + static Uns32T index_to_trackID(Uns32T lshID); // Convert lsh point index to audioDB trackID + static Uns32T index_to_trackPos(Uns32T lshID); // Convert lsh point index to audioDB trackPos (spos) + static Uns32T index_from_trackInfo(Uns32T, Uns32T); // Convert audioDB trackID and trackPos to an lsh point index + void initialize_exact_evalutation_queue(); + void index_insert_exact_evaluation_queue(Uns32T trackID, Uns32T qpos, Uns32T spos); + + // Web Services void startServer(); }; -#define O2_AUDIODB_INITIALIZERS \ - dim(0), \ - dbName(0), \ - inFile(0), \ - key(0), \ - trackFileName(0), \ - trackFile(0), \ - command(0), \ - output(0), \ - timesFileName(0), \ - timesFile(0), \ - powerFileName(0), \ - powerFile(0), \ - powerfd(0), \ - dbfid(0), \ - forWrite(false), \ - infid(0), \ - db(0), \ - indata(0), \ - dbH(0), \ - rng(0), \ - fileTable(0), \ - trackTable(0), \ - dataBuf(0), \ - l2normTable(0), \ - timesTable(0), \ - fileTableLength(0), \ - trackTableLength(0), \ - dataBufLength(0), \ - timesTableLength(0), \ - powerTableLength(0), \ - l2normTableLength(0), \ - verbosity(1), \ - nsamples(2000), \ - datasize(O2_DEFAULT_DATASIZE), \ - ntracks(O2_DEFAULT_NTRACKS), \ - datadim(O2_DEFAULT_DATADIM), \ - queryType(O2_POINT_QUERY), \ - pointNN(O2_DEFAULT_POINTNN), \ - trackNN(O2_DEFAULT_TRACKNN), \ - sequenceLength(16), \ - sequenceHop(1), \ - normalizedDistance(true), \ - queryPoint(0), \ - usingQueryPoint(0), \ - usingTimes(0), \ - usingPower(0), \ - isClient(0), \ - isServer(0), \ - port(0), \ - timesTol(0.1), \ - radius(0), \ - use_absolute_threshold(false), \ - absolute_threshold(0.0), \ - use_relative_threshold(false), \ - relative_threshold(0.0) +#define O2_AUDIODB_INITIALIZERS \ + dim(0), \ + dbName(0), \ + inFile(0), \ + key(0), \ + trackFileName(0), \ + trackFile(0), \ + command(0), \ + output(0), \ + timesFileName(0), \ + timesFile(0), \ + powerFileName(0), \ + powerFile(0), \ + powerfd(0), \ + dbfid(0), \ + lshfid(0), \ + forWrite(false), \ + infid(0), \ + db(0), \ + indata(0), \ + dbH(0), \ + rng(0), \ + fileTable(0), \ + trackTable(0), \ + trackOffsetTable(0), \ + dataBuf(0), \ + l2normTable(0), \ + timesTable(0), \ + fileTableLength(0), \ + trackTableLength(0), \ + dataBufLength(0), \ + timesTableLength(0), \ + powerTableLength(0), \ + l2normTableLength(0), \ + verbosity(1), \ + nsamples(2000), \ + datasize(O2_DEFAULT_DATASIZE), \ + ntracks(O2_DEFAULT_NTRACKS), \ + datadim(O2_DEFAULT_DATADIM), \ + queryType(O2_POINT_QUERY), \ + pointNN(O2_DEFAULT_POINTNN), \ + trackNN(O2_DEFAULT_TRACKNN), \ + sequenceLength(16), \ + sequenceHop(1), \ + normalizedDistance(true), \ + no_unit_norming(false), \ + queryPoint(0), \ + usingQueryPoint(0), \ + usingTimes(0), \ + usingPower(0), \ + isClient(0), \ + isServer(0), \ + port(0), \ + timesTol(0.1), \ + radius(0), \ + query_from_key(false), \ + query_from_key_index(O2_ERR_KEYNOTFOUND), \ + use_absolute_threshold(false), \ + absolute_threshold(0.0), \ + use_relative_threshold(false), \ + relative_threshold(0.0), \ + reporter(0), \ + exact_evaluation_queue(0), \ + lsh(0), \ + lsh_in_core(false), \ + lsh_use_u_functions(false), \ + lsh_exact(false), \ + lsh_param_k(0), \ + lsh_param_m(0), \ + lsh_param_N(0), \ + lsh_param_b(0), \ + lsh_param_ncols(0), \ + vv(0) +#endif
--- a/common.cpp Tue Jul 22 20:09:31 2008 +0000 +++ b/common.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -147,6 +147,14 @@ CHECKED_MMAP(double *, powerTable, dbH->powerTableOffset, powerTableLength); CHECKED_MMAP(double *, l2normTable, dbH->l2normTableOffset, l2normTableLength); } + + // build track offset table + trackOffsetTable = new off_t[dbH->numFiles]; + Uns32T cumTrack=0; + for(Uns32T k = 0; k < dbH->numFiles; k++){ + trackOffsetTable[k] = cumTrack; + cumTrack += trackTable[k] * dbH->dim; + } } void audioDB::initInputFile (const char *inFile) { @@ -187,7 +195,7 @@ } } -void audioDB::initTables(const char* dbName, const char* inFile = 0) { +void audioDB::initTables(const char* dbName, const char* inFile) { /* FIXME: initRNG() really logically belongs in the audioDB contructor. However, there are of the order of four constructors at the moment, and more to come from API implementation. Given @@ -196,5 +204,7 @@ database will need an RNG. -- CSR, 2008-07-02 */ initRNG(); initDBHeader(dbName); - initInputFile(inFile); + if(inFile) + initInputFile(inFile); } +
--- a/gengetopt.in Tue Jul 22 20:09:31 2008 +0000 +++ b/gengetopt.in Tue Jul 29 22:01:17 2008 +0000 @@ -23,7 +23,7 @@ option "output" - "output directory" string dependon="DUMP" default="audioDB.dump" optional option "L2NORM" L "unit norm vectors and norm all future inserts." dependon="database" optional option "POWER" P "turn on power flag for database." dependon="database" optional - +option "INDEX" X "build an index for -d database at -R radius" dependon="database" dependon="radius" optional section "Database Information" sectiondesc="Information about databases." option "STATUS" S "output database information to stdout." dependon="database" optional @@ -37,7 +37,7 @@ option "features" f "binary series of vectors file {int sz:ieee double[][sz]:eof}." string typestr="filename" dependon="database" optional option "times" t "list of time points (ascii) for feature vectors." string typestr="filename" dependon="features" optional option "power" w "binary power feature file." string typestr="filename" dependon="database" optional -option "key" k "unique identifier associated with features." string typestr="identifier" dependon="features" optional +option "key" k "unique identifier associated with features." string typestr="identifier" optional text "" option "BATCHINSERT" B "add feature vectors named in a --featureList file (with optional keys in a --keyList file) to the named database." dependon="featureList" optional option "featureList" F "text file containing list of binary feature vector files to process, one per track" string typestr="filename" dependon="database" optional @@ -45,12 +45,9 @@ option "powerList" W "text file containing list of binary power feature file." string typestr="filename" dependon="database" optional option "keyList" K "text file containing list of unique identifiers associated with --features." string typestr="filename" optional -option "sequencelength" l "length of sequences for sequence search." int typestr="length" default="16" optional +section "Database Search" sectiondesc="These commands control the retrieval behaviour.\n" - -section "Database Search" sectiondesc="Thse commands control the retrieval behaviour.\n" - -option "QUERY" Q "content-based search on --database using --features as a query. Optionally restrict the search to those tracks identified in a --keyList." values="point","track","sequence","nsequence","onetoonensequence" typestr="searchtype" dependon="database" dependon="features" optional +option "QUERY" Q "content-based search on --database using --features as a query. Optionally restrict the search to those tracks identified in a --keyList." values="point","track","sequence","nsequence","onetoonensequence" typestr="searchtype" dependon="database" optional option "qpoint" p "ordinal position of query start point in --features file." int typestr="position" default="0" optional option "exhaustive" e "exhaustive search: iterate through all query vectors in search. Overrides --qpoint." flag off optional hidden option "pointnn" n "number of point nearest neighbours to use in retrieval." int typestr="numpoints" default="10" optional @@ -58,15 +55,33 @@ option "expandfactor" x "time compress/expand factor of result length to query length [1.0 .. 100.0]." double default="1.1" optional hidden option "rotate" o "rotate query vectors for rotationally invariant search." flag off optional hidden option "resultlength" r "maximum length of the result list." int typestr="length" default="10" optional +option "sequencelength" l "length of sequences for sequence search." int typestr="length" default="16" optional option "sequencehop" h "hop size of sequence window for sequence search." int typestr="hop" default="1" dependon="QUERY" optional -option "absolute-threshold" - "absolute power threshold for consideration of query or target sequence (in Bels)" double dependon="QUERY" optional +option "absolute-threshold" - "absolute power threshold for consideration of query or target sequence (in Bels)" double optional option "relative-threshold" - "relative power threshold between query and target sequence (in Bels)" double dependon="QUERY" optional +section "Locality-sensitive hashing (LSH) parameters" sectiondesc="These parameters control LSH indexing and retrieval\n" + +option "lsh_w" - "width of LSH hash-function bins. " double default="4.0" dependon="INDEX" optional hidden +option "lsh_k" - "even number of independent hash functions to employ with LSH" int typestr="size" default="8" dependon="INDEX" optional +option "lsh_m" - "number of hash tables is m(m-1)/2" int typestr="size" default="5" dependon="INDEX" optional +option "lsh_N" - "number of rows per hash tables" int typestr="size" default="100000" dependon="INDEX" optional +option "lsh_b" - "number of tracks per indexing iteration" int typestr="size" default="500" dependon="INDEX" optional +option "lsh_ncols" - "number of columns (collisions) to allocate in LSH serialization" int typestr="size" default="250" dependon="INDEX" optional +option "lsh_exact" - "use exact evaluation of points retrieved by LSH." flag off dependon="QUERY" optional +option "lsh_inCore" - "LSH hash tables resident in core (INDEX/QUERY)" flag on optional +option "lsh_use_u_functions" - "use m independent hash functions combinatorically to approximate L independent hash functions." flag off optional + +section "Normalization control parameters" sectiondesc="These parameters control the normalization of feaures at query time\n" + +option "no_unit_norming" - "do not unit norm features when querying an L2Norm databases." flag off optional + section "Web Services" sectiondesc="These commands enable the database process to establish a connection via the internet and operate as separate client and server processes.\n" + option "SERVER" s "run as standalone web service on named port." int typestr="port" default="14475" optional option "client" c "run as a client using named host service." string typestr="hostname:port" optional text " Copyright (c) 2007 Michael Casey, Christophe Rhodes - Goldsmiths, University of London" + Goldsmiths, University of London"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/index.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,504 @@ +// LSH indexing +// +// Construct a persistent LSH table structure +// Store at the same location as dbName +// Naming convention: +// dbName.lsh.${radius}.${sequenceLength} +// +// +// Author: Michael Casey +// Date: 23 June 2008 + +#include "audioDB.h" +#include "ReporterBase.h" + + +/************************* LSH point index to audioDB conversion *****************/ +Uns32T audioDB::index_to_trackID(Uns32T lshID){ + return lshID>>LSH_N_POINT_BITS; +} + +Uns32T audioDB::index_to_trackPos(Uns32T lshID){ + return lshID&LSH_POINT_MASK; +} + +Uns32T audioDB::index_from_trackInfo(Uns32T trackID, Uns32T spos){ + return (trackID << LSH_N_POINT_BITS) | spos; +} + +/************************* LSH indexing and query initialization *****************/ + +char* audioDB::index_get_name(const char*dbName, double radius, Uns32T sequenceLength){ + char* indexName = new char[MAXSTR]; + // Attempt to make new file + if(strlen(dbName) > (MAXSTR - 32)) + error("dbName is too long for LSH index filename appendages"); + strncpy(indexName, dbName, MAXSTR); + sprintf(indexName+strlen(dbName), ".lsh.%019.9f.%d", radius, sequenceLength); + return indexName; +} + +// return true if index exists else return false +int audioDB::index_exists(const char* dbName, double radius, Uns32T sequenceLength){ + // Test to see if file exists + char* indexName = index_get_name(dbName, radius, sequenceLength); + lshfid = open (indexName, O_RDONLY); + delete[] indexName; + close(lshfid); + + if(lshfid<0) + return false; + else + return true; +} + +vector<vector<float> >* audioDB::index_initialize_shingles(Uns32T sz){ + if(vv) + delete vv; + vv = new vector<vector<float> >(sz); + for(Uns32T i=0 ; i < sz ; i++) + (*vv)[i]=vector<float>(dbH->dim*sequenceLength); // allocate shingle storage + return vv; +} + +/******************** LSH indexing audioDB database access forall s \in {S} ***********************/ + +// Prepare the AudioDB database for read access and allocate auxillary memory +void audioDB::index_initialize(double **snp, double **vsnp, double **spp, double **vspp, Uns32T *dvp) { + *dvp = dbH->length / (dbH->dim * sizeof(double)); // number of database vectors + *snp = new double[*dvp]; // songs norm pointer: L2 norm table for each vector + + double *snpp = *snp, *sppp = 0; + memcpy(*snp, l2normTable, *dvp * sizeof(double)); + + if (!(dbH->flags & O2_FLAG_POWER)) { + error("database not power-enabled", dbName); + } + *spp = new double[*dvp]; // song powertable pointer + sppp = *spp; + memcpy(*spp, powerTable, *dvp * sizeof(double)); + + for(Uns32T i = 0; i < dbH->numFiles; i++){ + if(trackTable[i] >= sequenceLength) { + sequence_sum(snpp, trackTable[i], sequenceLength); + sequence_sqrt(snpp, trackTable[i], sequenceLength); + + sequence_sum(sppp, trackTable[i], sequenceLength); + sequence_average(sppp, trackTable[i], sequenceLength); + } + snpp += trackTable[i]; + sppp += trackTable[i]; + } + + *vsnp = *snp; + *vspp = *spp; + + // Move the feature vector read pointer to start of fetures in database + lseek(dbfid, dbH->dataOffset, SEEK_SET); +} + + +/************************ LSH indexing ***********************************/ +void audioDB::index_index_db(const char* dbName){ + + char* newIndexName; + double *fvp = 0, *sNorm = 0, *snPtr = 0, *sPower = 0, *spPtr = 0; + Uns32T dbVectors = 0; + + printf("INDEX: initializing header\n"); + // Check if audioDB exists, initialize header and open database for read + forWrite = false; + initDBHeader(dbName); + + newIndexName = index_get_name(dbName, radius, sequenceLength); + + // Set unit norming flag override + audioDB::normalizedDistance = !audioDB::no_unit_norming; + + printf("INDEX: dim %d\n", dbH->dim); + printf("INDEX: R %f\n", radius); + printf("INDEX: seqlen %d\n", sequenceLength); + printf("INDEX: lsh_w %f\n", lsh_param_w); + printf("INDEX: lsh_k %d\n", lsh_param_k); + printf("INDEX: lsh_m %d\n", lsh_param_m); + printf("INDEX: lsh_N %d\n", lsh_param_N); + printf("INDEX: lsh_b %d\n", lsh_param_b); + printf("INDEX: normalized? %s\n", normalizedDistance?"true":"false"); + fflush(stdout); + + + index_initialize(&sNorm, &snPtr, &sPower, &spPtr, &dbVectors); + + if((lshfid = open(newIndexName,O_RDONLY))<0){ + printf("INDEX: constructing new LSH index\n"); + printf("INDEX: making index file %s\n", newIndexName); + fflush(stdout); + // Construct new LSH index + lsh = new LSH((float)lsh_param_w, lsh_param_k, + lsh_param_m, + (Uns32T)(sequenceLength*dbH->dim), + lsh_param_N, + lsh_param_ncols, + (float)radius); + assert(lsh); + + Uns32T endTrack = lsh_param_b; + if( endTrack > dbH->numFiles) + endTrack = dbH->numFiles; + // Insert up to lsh_param_b tracks + index_insert_tracks(0, endTrack, &fvp, &sNorm, &snPtr, &sPower, &spPtr); + lsh->serialize(newIndexName, lsh_in_core?O2_SERIAL_FILEFORMAT2:O2_SERIAL_FILEFORMAT1); + + // Clean up + delete lsh; + close(lshfid); + } + + // Attempt to open LSH file + if((lshfid = open(newIndexName,O_RDONLY))>0){ + printf("INDEX: merging with existing LSH index\n"); + fflush(stdout); + + // Get the lsh header info and find how many tracks are inserted already + lsh = new LSH(newIndexName, false); // lshInCore=false to avoid loading hashTables here + assert(lsh); + Uns32T maxs = index_to_trackID(lsh->get_maxp())+1; + delete lsh; + + // This allows for updating index after more tracks are inserted into audioDB + for(Uns32T startTrack = maxs; startTrack < dbH->numFiles; startTrack+=lsh_param_b){ + + Uns32T endTrack = startTrack + lsh_param_b; + if( endTrack > dbH->numFiles) + endTrack = dbH->numFiles; + printf("Indexing track range: %d - %d\n", startTrack, endTrack); + fflush(stdout); + lsh = new LSH(newIndexName, lsh_in_core); // Initialize core memory for LSH tables + assert(lsh); + + // Insert up to lsh_param_b database tracks + index_insert_tracks(startTrack, endTrack, &fvp, &sNorm, &snPtr, &sPower, &spPtr); + + // Serialize to file + lsh->serialize(newIndexName, lsh_in_core?O2_SERIAL_FILEFORMAT2:O2_SERIAL_FILEFORMAT1); // Serialize core LSH heap to disk + delete lsh; + } + + close(lshfid); + printf("INDEX: done constructing LSH index.\n"); + fflush(stdout); + + } + else{ + error("Something's wrong with LSH index file"); + exit(1); + } + + + delete[] newIndexName; + if(sNorm) + delete[] sNorm; + if(sPower) + delete[] sPower; + + +} + +void audioDB::index_insert_tracks(Uns32T start_track, Uns32T end_track, + double** fvpp, double** sNormpp,double** snPtrp, + double** sPowerp, double** spPtrp){ + size_t nfv = 0; + double* fvp = 0; // Keep pointer for memory allocation and free() for track data + Uns32T trackID = 0; + + VERB_LOG(1, "indexing tracks..."); + + + for(trackID = start_track ; trackID < end_track ; trackID++ ){ + read_data(trackID, &fvp, &nfv); // over-writes fvp and nfv + *fvpp = fvp; // Protect memory allocation and free() for track data + if(!index_insert_track(trackID, fvpp, snPtrp, spPtrp)) + break; + } + std::cout << "finished inserting." << endl; +} + +int audioDB::index_insert_track(Uns32T trackID, double** fvpp, double** snpp, double** sppp){ + // Loop over the current input track's vectors + Uns32T numVecs = (trackTable[trackID]>O2_MAXTRACKLEN?O2_MAXTRACKLEN:trackTable[trackID]) - sequenceLength + 1 ; + vv = index_initialize_shingles(numVecs); + + for( Uns32T pointID = 0 ; pointID < numVecs; pointID++ ) + index_make_shingle(vv, pointID, *fvpp, dbH->dim, sequenceLength); + + Uns32T numVecsAboveThreshold = index_norm_shingles(vv, *snpp, *sppp); + Uns32T collisionCount = index_insert_shingles(vv, trackID, *sppp); + float meanCollisionCount = numVecsAboveThreshold?(float)collisionCount/numVecsAboveThreshold:0; + + /* index_norm_shingles() only goes as far as the end of the + sequence, which is right, but the space allocated is for the + whole track. */ + + /* But numVecs will be <trackTable[track] if trackTable[track]>O2_MAXTRACKLEN + * So let's be certain the pointers are in the correct place + */ + + *snpp += trackTable[trackID]; + *sppp += trackTable[trackID]; + *fvpp += trackTable[trackID] * dbH->dim; + + std::cout << " n=" << trackTable[trackID] << " n'=" << numVecsAboveThreshold << " E[#c]=" << lsh->get_mean_collision_rate() << " E[#p]=" << meanCollisionCount << endl; + std::cout.flush(); + return true; +} + +Uns32T audioDB::index_insert_shingles(vector<vector<float> >* vv, Uns32T trackID, double* spp){ + Uns32T collisionCount = 0; + cout << "[" << trackID << "]" << fileTable+trackID*O2_FILETABLE_ENTRY_SIZE; + for( Uns32T pointID=0 ; pointID < (*vv).size(); pointID++) + if(!use_absolute_threshold || (use_absolute_threshold && (*spp++ >= absolute_threshold))) + collisionCount += lsh->insert_point((*vv)[pointID], index_from_trackInfo(trackID, pointID)); + return collisionCount; +} + +/********************* LSH shingle construction ***************************/ + +// Construct shingles out of a feature matrix +// inputs: +// idx is vector index in feature matrix +// fvp is base feature matrix pointer double* [numVecs x dbH->dim] +// +// pre-conditions: +// dbH->dim +// sequenceLength +// idx < numVectors - sequenceLength + 1 +// +// post-conditions: +// (*vv)[idx] contains a shingle with dbH->dim*sequenceLength float values + +void audioDB::index_make_shingle(vector<vector<float> >* vv, Uns32T idx, double* fvp, Uns32T dim, Uns32T seqLen){ + assert(idx<(*vv).size()); + vector<float>::iterator ve = (*vv)[idx].end(); + vi=(*vv)[idx].begin(); // shingle iterator + // First feature vector in shingle + if(idx==0){ + while(vi!=ve) + *vi++ = (float)(*fvp++); + } + // Not first feature vector in shingle + else{ + vector<float>::iterator ui=(*vv)[idx-1].begin() + dim; // previous shingle iterator + // Previous seqLen-1 dim-vectors + while(vi!=ve-dim) + *vi++=*ui++; + // Move data pointer to next feature vector + fvp += ( seqLen + idx - 1 ) * dim ; + // New d-vector + while(vi!=ve) + *vi++ = (float)(*fvp++); + } +} + +// norm shingles +// in-place norming, no deletions +// If using power, return number of shingles above power threshold +int audioDB::index_norm_shingles(vector<vector<float> >* vv, double* snp, double* spp){ + int z = 0; // number of above-threshold shingles + float l2norm; + double power; + float oneOverRadius = 1./(float)sqrt(radius); // Passed radius is really radius^2 + float oneOverSqrtl2NormDivRad = oneOverRadius; + if(!spp) + error("LSH indexing and query requires a power feature using -w or -W"); + Uns32T shingleSize = sequenceLength*dbH->dim; + for(Uns32T a=0; a<(*vv).size(); a++){ + l2norm = (float)(*snp++); + if(audioDB::normalizedDistance) + oneOverSqrtl2NormDivRad = (1./l2norm)*oneOverRadius; + + for(Uns32T b=0; b < shingleSize ; b++) + (*vv)[a][b]*=oneOverSqrtl2NormDivRad; + + power = *spp++; + if(use_absolute_threshold){ + if ( power >= absolute_threshold ) + z++; + } + else + z++; + } + return z; +} + + +/*********************** LSH retrieval ****************************/ + + +// return true if indexed query performed else return false +int audioDB::index_init_query(const char* dbName){ + + if(!(index_exists(dbName, radius, sequenceLength))) + return false; + + char* indexName = index_get_name(dbName, radius, sequenceLength); + + // Test to see if file exists + if((lshfid = open (indexName, O_RDONLY)) < 0){ + delete[] indexName; + return false; + } + + printf("INDEX: initializing header\n"); + + lsh = new LSH(indexName, false); // Get the header only here + assert(lsh); + sequenceLength = lsh->get_lshHeader()->dataDim / dbH->dim; // shingleDim / vectorDim + + + if( fabs(radius - lsh->get_radius())>fabs(O2_DISTANCE_TOLERANCE)) + printf("*** Warning: adb_radius (%f) != lsh_radius (%f) ***\n", radius, lsh->get_radius()); + + printf("INDEX: dim %d\n", dbH->dim); + printf("INDEX: R %f\n", lsh->get_radius()); + printf("INDEX: seqlen %d\n", sequenceLength); + printf("INDEX: w %f\n", lsh->get_lshHeader()->get_binWidth()); + printf("INDEX: k %d\n", lsh->get_lshHeader()->get_numFuns()); + printf("INDEX: L (m*(m-1))/2 %d\n", lsh->get_lshHeader()->get_numTables()); + printf("INDEX: N %d\n", lsh->get_lshHeader()->get_numRows()); + printf("INDEX: s %d\n", index_to_trackID(lsh->get_maxp())); + printf("INDEX: Opened LSH index file %s\n", indexName); + fflush(stdout); + + // Check to see if we are loading hash tables into core, and do so if true + if((lsh->get_lshHeader()->flags&O2_SERIAL_FILEFORMAT2) || lsh_in_core){ + printf("INDEX: loading hash tables into core %s\n", (lsh->get_lshHeader()->flags&O2_SERIAL_FILEFORMAT2)?"FORMAT2":"FORMAT1"); + delete lsh; + lsh = new LSH(indexName, true); + } + + delete[] indexName; + return true; +} + +// *Static* approximate NN point reporter callback method for lshlib +void audioDB::index_add_point_approximate(void* instancePtr, Uns32T pointID, Uns32T qpos, float dist){ + assert(instancePtr); // We need an instance for this callback + audioDB* myself = (audioDB*) instancePtr; // Use explicit cast to recover "this" instance + Uns32T trackID = index_to_trackID(pointID); + Uns32T spos = index_to_trackPos(pointID); + // Skip identity in query_from_key + if( !myself->query_from_key || (myself->query_from_key && ( trackID != myself->query_from_key_index )) ) + myself->reporter->add_point(trackID, qpos, spos, dist); +} + +// *Static* exact NN point reporter callback method for lshlib +// Maintain a queue of points to pass to query_points() for exact evaluation +void audioDB::index_add_point_exact(void* instancePtr, Uns32T pointID, Uns32T qpos, float dist){ + assert(instancePtr); // We need an instance for this callback + audioDB* myself = (audioDB*) instancePtr; // Use explicit cast to recover "this" instance + Uns32T trackID = index_to_trackID(pointID); + Uns32T spos = index_to_trackPos(pointID); + // Skip identity in query_from_key + if( !myself->query_from_key || (myself->query_from_key && ( trackID != myself->query_from_key_index )) ) + myself->index_insert_exact_evaluation_queue(trackID, qpos, spos); +} + +void audioDB::initialize_exact_evalutation_queue(){ + if(exact_evaluation_queue) + delete exact_evaluation_queue; + exact_evaluation_queue = new priority_queue<PointPair, std::vector<PointPair>, std::less<PointPair> >; +} + +void audioDB::index_insert_exact_evaluation_queue(Uns32T trackID, Uns32T qpos, Uns32T spos){ + PointPair p(trackID, qpos, spos); + exact_evaluation_queue->push(p); +} + +// return 0: if index does not exist +// return nqv: if index exists +int audioDB::index_query_loop(const char* dbName, Uns32T queryIndex) { + + unsigned int numVectors; + double *query, *query_data; + double *qNorm, *qnPtr, *qPower = 0, *qpPtr = 0; + double meanQdur; + void (*add_point_func)(void*,Uns32T,Uns32T,float); + + // Set the point-reporter callback based on the value of lsh_exact + if(lsh_exact){ + initialize_exact_evalutation_queue(); + add_point_func = &index_add_point_exact; + } + else + add_point_func = &index_add_point_approximate; + + if(!index_init_query(dbName)) // sets-up LSH index structures for querying + return 0; + + char* database = index_get_name(dbName, radius, sequenceLength); + + if(query_from_key) + set_up_query_from_key(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors, queryIndex); + else + set_up_query(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors); // get query vectors + + VERB_LOG(1, "retrieving tracks..."); + + assert(pointNN>0 && pointNN<=O2_MAXNN); + assert(trackNN>0 && trackNN<=O2_MAXNN); + + gettimeofday(&tv1, NULL); + // query vector index + Uns32T Nq = (numVectors>O2_MAXTRACKLEN?O2_MAXTRACKLEN:numVectors) - sequenceLength + 1; + vv = index_initialize_shingles(Nq); // allocate memory to copy query vectors to shingles + cout << "Nq=" << Nq; cout.flush(); + // Construct shingles from query features + for( Uns32T pointID = 0 ; pointID < Nq ; pointID++ ) + index_make_shingle(vv, pointID, query, dbH->dim, sequenceLength); + + // Normalize query vectors + Uns32T numVecsAboveThreshold = index_norm_shingles( vv, qnPtr, qpPtr ); + cout << " Nq'=" << numVecsAboveThreshold << endl; cout.flush(); + + // Nq contains number of inspected points in query file, + // numVecsAboveThreshold is number of points with power >= absolute_threshold + double* qpp = qpPtr; // Keep original qpPtr for possible exact evaluation + if(usingQueryPoint && numVecsAboveThreshold){ + if((lsh->get_lshHeader()->flags&O2_SERIAL_FILEFORMAT2) || lsh_in_core) + lsh->retrieve_point((*vv)[0], queryPoint, add_point_func, (void*)this); + else + lsh->serial_retrieve_point(database, (*vv)[0], queryPoint, add_point_func, (void*)this); + } + else if(numVecsAboveThreshold) + for( Uns32T pointID = 0 ; pointID < Nq; pointID++ ) + if(!use_absolute_threshold || (use_absolute_threshold && (*qpp++ >= absolute_threshold))) + if((lsh->get_lshHeader()->flags&O2_SERIAL_FILEFORMAT2) || lsh_in_core) + lsh->retrieve_point((*vv)[pointID], pointID, add_point_func, (void*)this); + else + lsh->serial_retrieve_point(database, (*vv)[pointID], pointID, add_point_func, (void*)this); + + if(lsh_exact) + // Perform exact distance computation on point pairs in exact_evaluation_queue + query_loop_points(query, qnPtr, qpPtr, meanQdur, numVectors); + + 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)) + + // Close the index file + close(lshfid); + + // Clean up + if(query_data) + delete[] query_data; + if(qNorm) + delete[] qNorm; + if(qPower) + delete[] qPower; + if(database) + delete[] database; + + return Nq; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lshlib.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,1429 @@ +#include "lshlib.h" + +//#define __LSH_DUMP_CORE_TABLES__ +//#define USE_U_FUNCTIONS +//#define LSH_BLOCK_FULL_ROWS + +void err(char*s){cout << s << endl;exit(2);} + +Uns32T get_page_logn(){ + int pagesz = (int)sysconf(_SC_PAGESIZE); + return (Uns32T)log2((double)pagesz); +} + +unsigned align_up(unsigned x, unsigned w){ return ((x) + ((1<<w)-1) & ~((1<<w)-1)); } + +void H::error(const char* a, const char* b, const char *sysFunc) { + cerr << a << ": " << b << endl; + if (sysFunc) { + perror(sysFunc); + } + exit(1); +} + +H::H(Uns32T kk, Uns32T mm, Uns32T dd, Uns32T NN, Uns32T CC): +#ifdef USE_U_FUNCTIONS + use_u_functions(true), +#else + use_u_functions(false), +#endif + bucketCount(0), + pointCount(0), + N(NN), + C(CC), + k(kk), + m(mm), + L(mm*(mm-1)/2), + d(dd) +{ + Uns32T j; + cout << "file size: ~" << (((unsigned long long)L*N*C*sizeof(SerialElementT))/1000000UL) << "MB" << endl; + if(((unsigned long long)L*N*C*sizeof(SerialElementT))>4000000000UL) + error("Maximum size of LSH file exceded: 12*L*N*C > 4000MB"); + else if(((unsigned long long)N*C*sizeof(SerialElementT))>1000000000UL) + cout << "warning: hash tables exceed 1000MB." << endl; + + if(m<2){ + m=2; + L=1; // check value of L + cout << "warning: setting m=2, L=1" << endl; + } + if(use_u_functions && k%2){ + k++; // make sure k is even + cout << "warning: setting k even" << endl; + } + __initialize_data_structures(); + for(j=0; j<L; j++) + for(kk=0; kk<k; kk++) { + r1[j][kk]=__randr(); // random 1..2^29 + r2[j][kk]=__randr(); // random 1..2^29 + } +} + +// Post constructor initialization +void H::__initialize_data_structures(){ + H::P = UH_PRIME_DEFAULT; + + /* FIXME: don't use time(); instead use /dev/random or similar */ + /* FIXME: write out the seed somewhere, so that we can get + repeatability */ +#ifdef MT19937 + init_genrand(time(NULL)); +#else + srand(time(NULL)); // seed random number generator +#endif + Uns32T i,j; + H::h = new bucket**[ H::L ]; + H::r1 = new Uns32T*[ H::L ]; + H::r2 = new Uns32T*[ H::L ]; + assert( H::h && H::r1 && H::r2 ); // failure + for( j = 0 ; j < H::L ; j++ ){ + H::r1[ j ] = new Uns32T[ H::k ]; + H::r2[ j ] = new Uns32T[ H::k ]; + assert( H::r1[j] && H::r2[j] ); // failure + } + + for( j = 0 ; j < H::L ; j++ ){ + H::h[j] = new bucket*[ H::N ]; + assert( H::h[j] ); + for( i = 0 ; i < H::N ; i++) + H::h[j][i] = 0; + } +} + +// Destruct hash tables +H::~H(){ + Uns32T i,j; + for( j=0 ; j < H::L ; j++ ){ + delete[] H::r1[ j ]; + delete[] H::r2[ j ]; + for(i = 0; i< H::N ; i++) + delete H::h[ j ][ i ]; + delete[] H::h[ j ]; + } + delete[] H::r1; + delete[] H::r2; + delete[] H::h; +} + + +// make hash value \in Z +void H::__generate_hash_keys(Uns32T*g,Uns32T* r1, Uns32T* r2){ + H::t1 = __computeProductModDefaultPrime( g, r1, H::k ) % H::N; + H::t2 = __computeProductModDefaultPrime( g, r2, H::k ); + +} + +#define CR_ASSERT(b){if(!(b)){fprintf(stderr, "ASSERT failed on line %d, file %s.\n", __LINE__, __FILE__); exit(1);}} + +// Computes (a.b) mod UH_PRIME_DEFAULT +inline Uns32T H::__computeProductModDefaultPrime(Uns32T *a, Uns32T *b, IntT size){ + LongUns64T h = 0; + + for(IntT i = 0; i < size; i++){ + h = h + (LongUns64T)a[i] * (LongUns64T)b[i]; + h = (h & TWO_TO_32_MINUS_1) + 5 * (h >> 32); + if (h >= UH_PRIME_DEFAULT) { + h = h - UH_PRIME_DEFAULT; + } + CR_ASSERT(h < UH_PRIME_DEFAULT); + } + return h; +} + +Uns32T H::bucket_insert_point(bucket **pp){ + Uns32T collisionCount = 0; + if(!*pp){ + *pp = new bucket(); +#ifdef LSH_BLOCK_FULL_ROWS + (*pp)->t2 = 0; // Use t2 as a collision counter for the row + (*pp)->next = new bucket(); +#endif + } +#ifdef LSH_BLOCK_FULL_ROWS + collisionCount = (*pp)->t2; + if(collisionCount < H::C){ // Block if row is full + (*pp)->t2++; // Increment collision counter + pointCount++; + collisionCount++; + __bucket_insert_point((*pp)->next); // First bucket holds collision count + } +#else + pointCount++; + __bucket_insert_point(*pp); // No collision count storage +#endif + return collisionCount; +} + +void H::__bucket_insert_point(bucket* p){ + if(p->t2 == IFLAG){ // initialization flag, is it in the domain of t2? + p->t2 = H::t2; + bucketCount++; // Record start of new point-locale collision chain + p->snext = new sbucket(); + __sbucket_insert_point(p->snext); + return; + } + + if(p->t2 == H::t2){ + __sbucket_insert_point(p->snext); + return; + } + + if(p->next){ + __bucket_insert_point(p->next); + } + + else{ + p->next = new bucket(); + __bucket_insert_point(p->next); + } + +} + +void H::__sbucket_insert_point(sbucket* p){ + if(p->pointID==IFLAG){ + p->pointID = H::p; + return; + } + + // Search for pointID + if(p->snext){ + __sbucket_insert_point(p->snext); + } + else{ + // Make new point collision bucket at end of list + p->snext = new sbucket(); + __sbucket_insert_point(p->snext); + } +} + +inline bucket** H::__get_bucket(int j){ + return *(h+j); +} + +// hash functions G +G::G(float ww, Uns32T kk,Uns32T mm, Uns32T dd, Uns32T NN, Uns32T CC, float r): + H(kk,mm,dd,NN,CC), + w(ww), + radius(r), + maxp(0), + calling_instance(0), + add_point_callback(0), + lshHeader(0) +{ + Uns32T j; +#ifdef USE_U_FUNCTIONS + G::A = new float**[ H::m ]; // m x k x d random projectors + G::b = new float*[ H::m ]; // m x k random biases +#else + G::A = new float**[ H::L ]; // m x k x d random projectors + G::b = new float*[ H::L ]; // m x k random biases +#endif + G::g = new Uns32T*[ H::L ]; // L x k random projections + assert( G::g && G::A && G::b ); // failure +#ifdef USE_U_FUNCTIONS + // Use m \times u_i functions \in R^{(k/2) \times (d)} + // Combine to make L=m(m-1)/2 hash functions \in R^{k \times d} + for( j = 0; j < H::m ; j++ ){ // m functions u_i(v) + G::A[j] = new float*[ H::k/2 ]; // k/2 x d 2-stable distribution coefficients + G::b[j] = new float[ H::k/2 ]; // bias + assert( G::A[j] && G::b[j] ); // failure + for( kk = 0; kk < H::k/2 ; kk++ ){ + G::A[j][kk] = new float[ H::d ]; + assert( G::A[j][kk] ); // failure + for(Uns32T i = 0 ; i < H::d ; i++ ) + G::A[j][kk][i] = randn(); // Normal + G::b[j][kk] = ranf()*G::w; // Uniform + } + } +#else + // Use m \times u_i functions \in R^{k \times (d)} + // Combine to make L=m(m-1)/2 hash functions \in R^{k \times d} + for( j = 0; j < H::L ; j++ ){ // m functions u_i(v) + G::A[j] = new float*[ H::k ]; // k x d 2-stable distribution coefficients + G::b[j] = new float[ H::k ]; // bias + assert( G::A[j] && G::b[j] ); // failure + for( kk = 0; kk < H::k ; kk++ ){ + G::A[j][kk] = new float[ H::d ]; + assert( G::A[j][kk] ); // failure + for(Uns32T i = 0 ; i < H::d ; i++ ) + G::A[j][kk][i] = randn(); // Normal + G::b[j][kk] = ranf()*G::w; // Uniform + } + } +#endif + + for( j = 0 ; j < H::L ; j++ ){ // L functions g_j(u_a, u_b) a,b \in nchoosek(m,2) + G::g[j] = new Uns32T[ H::k ]; // k x 32-bit hash values, gj(v)=[x0 x1 ... xk-1] xk \in Z + assert( G::g[j] ); + } + + initialize_partial_functions(); // m partially evaluated hash functions +} + +// Serialize from file LSH constructor +// Read parameters from database file +// Load the hash functions, close the database +// Optionally load the LSH tables into head-allocated lists in core +G::G(char* filename, bool lshInCoreFlag): + calling_instance(0), + add_point_callback(0) +{ + int dbfid = unserialize_lsh_header(filename); + unserialize_lsh_functions(dbfid); + initialize_partial_functions(); + + // Format1 only needs unserializing if specifically requested + if(!(lshHeader->flags&O2_SERIAL_FILEFORMAT2) && lshInCoreFlag){ + unserialize_lsh_hashtables_format1(dbfid); + } + + // Format2 always needs unserializing + if(lshHeader->flags&O2_SERIAL_FILEFORMAT2 && lshInCoreFlag){ + unserialize_lsh_hashtables_format2(dbfid); + } + + close(dbfid); +} + +void G::initialize_partial_functions(){ + +#ifdef USE_U_FUNCTIONS + uu = vector<vector<Uns32T> >(H::m); + for( Uns32T aa=0 ; aa < H::m ; aa++ ) + uu[aa] = vector<Uns32T>( H::k/2 ); +#else + uu = vector<vector<Uns32T> >(H::L); + for( Uns32T aa=0 ; aa < H::L ; aa++ ) + uu[aa] = vector<Uns32T>( H::k ); +#endif +} + + +// Generate z ~ N(0,1) +float G::randn(){ +// Box-Muller + float x1, x2; + do{ + x1 = ranf(); + } while (x1 == 0); // cannot take log of 0 + x2 = ranf(); + float z; + z = sqrtf(-2.0 * logf(x1)) * cosf(2.0 * M_PI * x2); + return z; +} + +float G::ranf(){ +#ifdef MT19937 + return (float) genrand_real2(); +#else + return (float)( (double)rand() / ((double)(RAND_MAX)+(double)(1)) ); +#endif +} + +// range is 1..2^29 +/* FIXME: that looks like an ... odd range. Still. */ +Uns32T H::__randr(){ +#ifdef MT19937 + return (Uns32T)((genrand_int32() >> 3) + 1); +#else + return (Uns32T) ((rand() >> 2) + 1); +#endif +} + +G::~G(){ + Uns32T j,kk; +#ifdef USE_U_FUNCTIONS + for( j = 0 ; j < H::m ; j++ ){ + for( kk = 0 ; kk < H::k/2 ; kk++ ) + delete[] A[j][kk]; + delete[] A[j]; + } + delete[] A; + for( j = 0 ; j < H::m ; j++ ) + delete[] b[j]; + delete[] b; +#else + for( j = 0 ; j < H::L ; j++ ){ + for( kk = 0 ; kk < H::k ; kk++ ) + delete[] A[j][kk]; + delete[] A[j]; + } + delete[] A; + for( j = 0 ; j < H::L ; j++ ) + delete[] b[j]; + delete[] b; +#endif + + for( j = 0 ; j < H::L ; j++ ) + delete[] g[j]; + delete[] g; + delete lshHeader; +} + +// Compute all hash functions for vector v +// #ifdef USE_U_FUNCTIONS use Combination of m \times h_i \in R^{(k/2) \times d} +// to make L \times g_j functions \in Z^k +void G::compute_hash_functions(vector<float>& v){ // v \in R^d + float iw = 1. / G::w; // hash bucket width + Uns32T aa, kk; + if( v.size() != H::d ) + error("v.size != H::d","","compute_hash_functions"); // check input vector dimensionality + double tmp = 0; + float *pA, *pb; + Uns32T *pg; + int dd; + vector<float>::iterator vi; + vector<Uns32T>::iterator ui; + +#ifdef USE_U_FUNCTIONS + Uns32T bb; + // Store m dot products to expand + for( aa=0; aa < H::m ; aa++ ){ + ui = uu[aa].begin(); + for( kk = 0 ; kk < H::k/2 ; kk++ ){ + pb = *( G::b + aa ) + kk; + pA = * ( * ( G::A + aa ) + kk ); + dd = H::d; + tmp = 0.; + vi = v.begin(); + while( dd-- ) + tmp += *pA++ * *vi++; // project + tmp += *pb; // translate + tmp *= iw; // scale + *ui++ = (Uns32T) floor(tmp); // floor + } + } + // Binomial combinations of functions u_{a,b} \in Z^{(k/2) \times d} + Uns32T j; + for( aa=0, j=0 ; aa < H::m-1 ; aa++ ) + for( bb = aa + 1 ; bb < H::m ; bb++, j++ ){ + pg= *( G::g + j ); // L \times functions g_j(v) \in Z^k + // u_1 \in Z^{(k/2) \times d} + ui = uu[aa].begin(); + kk=H::k/2; + while( kk-- ) + *pg++ = *ui++; // hash function g_j(v)=[x1 x2 ... x(k/2)]; xk \in Z + // u_2 \in Z^{(k/2) \times d} + ui = uu[bb].begin(); + kk=H::k/2; + while( kk--) + *pg++ = *ui++; // hash function g_j(v)=[x(k/2+1) x(k/2+2) ... xk]; xk \in Z + } +#else + for( aa=0; aa < H::L ; aa++ ){ + ui = uu[aa].begin(); + for( kk = 0 ; kk < H::k ; kk++ ){ + pb = *( G::b + aa ) + kk; + pA = * ( * ( G::A + aa ) + kk ); + dd = H::d; + tmp = 0.; + vi = v.begin(); + while( dd-- ) + tmp += *pA++ * *vi++; // project + tmp += *pb; // translate + tmp *= iw; // scale + *ui++ = (Uns32T) (floor(tmp)); // floor + } + } + // Compute hash functions + for( aa=0 ; aa < H::L ; aa++ ){ + pg= *( G::g + aa ); // L \times functions g_j(v) \in Z^k + // u_1 \in Z^{k \times d} + ui = uu[aa].begin(); + kk=H::k; + while( kk-- ) + *pg++ = *ui++; // hash function g_j(v)=[x1 x2 ... xk]; xk \in Z + } +#endif + +} + + +// single point insertion; inserted values are hash value and pointID +Uns32T G::insert_point(vector<float>& v, Uns32T pp){ + Uns32T collisionCount = 0; + H::p = pp; + if(pp>G::maxp) + G::maxp=pp; // Store highest pointID in database + compute_hash_functions( v ); + for(Uns32T j = 0 ; j < H::L ; j++ ){ // insertion + __generate_hash_keys( *( G::g + j ), *( H::r1 + j ), *( H::r2 + j ) ); + collisionCount += bucket_insert_point( *(h + j) + t1 ); + } + return collisionCount; +} + + +// batch insert for a point set +// inserted values are vector hash value and pointID starting at basePointID +void G::insert_point_set(vector<vector<float> >& vv, Uns32T basePointID){ + for(Uns32T point=0; point<vv.size(); point++) + insert_point(vv[point], basePointID+point); +} + +// point retrieval routine +void G::retrieve_point(vector<float>& v, Uns32T qpos, ReporterCallbackPtr add_point, void* caller){ + calling_instance = caller; + add_point_callback = add_point; + compute_hash_functions( v ); + for(Uns32T j = 0 ; j < H::L ; j++ ){ + __generate_hash_keys( *( G::g + j ), *( H::r1 + j ), *( H::r2 + j ) ); + if( bucket* bPtr = *(__get_bucket(j) + get_t1()) ) +#ifdef LSH_BLOCK_FULL_ROWS + bucket_chain_point( bPtr->next, qpos); +#else + bucket_chain_point( bPtr , qpos); +#endif + } +} + +void G::retrieve_point_set(vector<vector<float> >& vv, ReporterCallbackPtr add_point, void* caller){ + for(Uns32T qpos = 0 ; qpos < vv.size() ; qpos++ ) + retrieve_point(vv[qpos], qpos, add_point, caller); +} + +// export lsh tables to table structure on disk +// +// LSH TABLE STRUCTURE +// ---header 64 bytes --- +// [magic #tables #rows #cols elementSize databaseSize version flags dim #funs 0 0 0 0 0 0] +// +// ---random projections L x k x d float --- +// A[0][0][0] A[0][0][1] ... A[0][0][d-1] +// A[0][1][0] A[0][1][1] ... A[1][1][d-1] +// ... +// A[0][K-1][0] A[0][1][1] ... A[0][k-1][d-1] +// ... +// ... +// A[L-1][0][0] A[M-1][0][1] ... A[L-1][0][d-1] +// A[L-1][1][0] A[M-1][1][1] ... A[L-1][1][d-1] +// ... +// A[L-1][k-1][0] A[M-1][1][1] ... A[L-1][k-1][d-1] +// +// ---bias L x k float --- +// b[0][0] b[0][1] ... b[0][k-1] +// b[1][0] b[1][1] ... b[1][k-1] +// ... +// b[L-1][0] b[L-1][1] ... b[L-1][k-1] +// +// ---random r1 L x k float --- +// r1[0][0] r1[0][1] ... r1[0][k-1] +// r1[1][0] r1[1][1] ... r1[1][k-1] +// ... +// r1[L-1][0] r1[L-1][1] ... r1[L-1][k-1] +// +// ---random r2 L x k float --- +// r2[0][0] r2[0][1] ... r2[0][k-1] +// r2[1][0] r2[1][1] ... r2[1][k-1] +// ... +// r2[L-1][0] r2[L-1][1] ... r2[L-1][k-1] +// +// ---hash table 0: N x C x 8 --- +// [t2 pointID][t2 pointID]...[t2 pointID] +// [t2 pointID][t2 pointID]...[t2 pointID] +// ... +// [t2 pointID][t2 pointID]...[t2 pointID] +// +// ---hash table 1: N x C x 8 --- +// [t2 pointID][t2 pointID]...[t2 pointID] +// [t2 pointID][t2 pointID]...[t2 pointID] +// ... +// [t2 pointID][t2 pointID]...[t2 pointID] +// +// ... +// +// ---hash table L-1: N x C x 8 --- +// [t2 pointID][t2 pointID]...[t2 pointID] +// [t2 pointID][t2 pointID]...[t2 pointID] +// ... +// [t2 pointID][t2 pointID]...[t2 pointID] +// + +// Serial header constructors +SerialHeader::SerialHeader(){;} +SerialHeader::SerialHeader(float W, Uns32T L, Uns32T N, Uns32T C, Uns32T k, Uns32T d, float r, Uns32T p, Uns32T FMT): + lshMagic(O2_SERIAL_MAGIC), + binWidth(W), + numTables(L), + numRows(N), + numCols(C), + elementSize(O2_SERIAL_ELEMENT_SIZE), + version(O2_SERIAL_VERSION), + size(L * align_up(N * C * O2_SERIAL_ELEMENT_SIZE, get_page_logn()) // hash tables + + align_up(O2_SERIAL_HEADER_SIZE + // header + hash functions + L*k*( sizeof(float)*d+2*sizeof(Uns32T)+sizeof(float)),get_page_logn())), + flags(FMT), + dataDim(d), + numFuns(k), + radius(r), + maxp(p){;} // header + +float* G::get_serial_hashfunction_base(char* db){ + if(db&&lshHeader) + return (float*)(db+O2_SERIAL_HEADER_SIZE); + else return NULL; +} + +SerialElementT* G::get_serial_hashtable_base(char* db){ + if(db&&lshHeader) + return (SerialElementT*)(db+get_serial_hashtable_offset()); + else + return NULL; +} + +Uns32T G::get_serial_hashtable_offset(){ + if(lshHeader) + return align_up(O2_SERIAL_HEADER_SIZE + + L*lshHeader->numFuns*( sizeof(float)*lshHeader->dataDim+2*sizeof(Uns32T)+sizeof(float)),get_page_logn()); + else + return 0; +} + +void G::serialize(char* filename, Uns32T serialFormat){ + int dbfid; + char* db; + int dbIsNew=0; + + // Check requested serialFormat + if(!(serialFormat==O2_SERIAL_FILEFORMAT1 || serialFormat==O2_SERIAL_FILEFORMAT2)) + error("Unrecognized serial file format request: ", "serialize()"); + + // Test to see if file exists + if((dbfid = open (filename, O_RDONLY)) < 0) + // If it doesn't, then create the file (CREATE) + if(errno == ENOENT){ + // Create the file + std::cout << "Creating new serialized LSH database:" << filename << "..."; + std::cout.flush(); + serial_create(filename, serialFormat); + dbIsNew=1; + } + else + // The file can't be opened + error("Can't open the file", filename, "open"); + + // Load the on-disk header into core + dbfid = serial_open(filename, 1); // open for write + db = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 1);// get database pointer + serial_get_header(db); // read header + serial_munmap(db, O2_SERIAL_HEADER_SIZE); // drop mmap + + // Check compatibility of core and disk data structures + if( !serial_can_merge(serialFormat) ) + error("Incompatible core and serial LSH, data structure dimensions mismatch."); + + // For new LSH databases write the hashfunctions + if(dbIsNew) + serialize_lsh_hashfunctions(dbfid); + // Write the hashtables in the requested format + if(serialFormat == O2_SERIAL_FILEFORMAT1) + serialize_lsh_hashtables_format1(dbfid, !dbIsNew); + else + serialize_lsh_hashtables_format2(dbfid, !dbIsNew); + + if(!dbIsNew){ + db = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 1);// get database pointer + //serial_get_header(db); // read header + cout << "maxp = " << G::maxp << endl; + lshHeader->maxp=G::maxp; + // Default to FILEFORMAT1 + if(!(lshHeader->flags&O2_SERIAL_FILEFORMAT2)) + lshHeader->flags|=O2_SERIAL_FILEFORMAT2; + memcpy((char*)db, (char*)lshHeader, sizeof(SerialHeaderT)); + serial_munmap(db, O2_SERIAL_HEADER_SIZE); // drop mmap + } + + serial_close(dbfid); +} + +// Test to see if core structure and requested format is +// compatible with currently opened database +int G::serial_can_merge(Uns32T format){ + SerialHeaderT* that = lshHeader; + if( (format==O2_SERIAL_FILEFORMAT2 && !that->flags&O2_SERIAL_FILEFORMAT2) + || (format!=O2_SERIAL_FILEFORMAT2 && that->flags&O2_SERIAL_FILEFORMAT2) + || !( this->w == that->binWidth && + this->L == that->numTables && + this->N == that->numRows && + this->k == that->numFuns && + this->d == that->dataDim && + sizeof(SerialElementT) == that->elementSize && + this->radius == that->radius)){ + serial_print_header(format); + return 0; + } + else + return 1; +} + +// Used as an error message for serial_can_merge() +void G::serial_print_header(Uns32T format){ + std::cout << "Fc:" << format << " Fs:" << lshHeader->flags << endl; + std::cout << "Wc:" << w << " Ls:" << lshHeader->binWidth << endl; + std::cout << "Lc:" << L << " Ls:" << lshHeader->numTables << endl; + std::cout << "Nc:" << N << " Ns:" << lshHeader->numRows << endl; + std::cout << "kc:" << k << " ks:" << lshHeader->numFuns << endl; + std::cout << "dc:" << d << " ds:" << lshHeader->dataDim << endl; + std::cout << "sc:" << sizeof(SerialElementT) << " ss:" << lshHeader->elementSize << endl; + std::cout << "rc:" << this->radius << " rs:" << lshHeader->radius << endl; +} + +int G::serialize_lsh_hashfunctions(int fid){ + float* pf; + Uns32T *pu; + Uns32T x,y,z; + + db = serial_mmap(fid, get_serial_hashtable_offset(), 1);// get database pointer + pf = get_serial_hashfunction_base(db); + + // HASH FUNCTIONS + // Write the random projectors A[][][] +#ifdef USE_U_FUNCTIONS + for( x = 0 ; x < H::m ; x++ ) + for( y = 0 ; y < H::k/2 ; y++ ) +#else + for( x = 0 ; x < H::L ; x++ ) + for( y = 0 ; y < H::k ; y++ ) +#endif + for( z = 0 ; z < d ; z++ ) + *pf++ = A[x][y][z]; + + // Write the random biases b[][] +#ifdef USE_U_FUNCTIONS + for( x = 0 ; x < H::m ; x++ ) + for( y = 0 ; y < H::k/2 ; y++ ) +#else + for( x = 0 ; x < H::L ; x++ ) + for( y = 0 ; y < H::k ; y++ ) +#endif + *pf++=b[x][y]; + + pu = (Uns32T*)pf; + + // Write the Z projectors r1[][] + for( x = 0 ; x < H::L ; x++) + for( y = 0 ; y < H::k ; y++) + *pu++ = r1[x][y]; + + // Write the Z projectors r2[][] + for( x = 0 ; x < H::L ; x++) + for( y = 0; y < H::k ; y++) + *pu++ = r2[x][y]; + + serial_munmap(db, get_serial_hashtable_offset()); + return 1; +} + +int G::serialize_lsh_hashtables_format1(int fid, int merge){ + SerialElementT *pe, *pt; + Uns32T x,y; + + if( merge && !serial_can_merge(O2_SERIAL_FILEFORMAT1) ) + error("Cannot merge core and serial LSH, data structure dimensions mismatch."); + + Uns32T hashTableSize=sizeof(SerialElementT)*lshHeader->numRows*lshHeader->numCols; + Uns32T colCount, meanColCount, colCountN, maxColCount, minColCount; + // Write the hash tables + for( x = 0 ; x < H::L ; x++ ){ + std::cout << (merge ? "merging":"writing") << " hash table " << x << " FORMAT1..."; + std::cout.flush(); + // memory map a single hash table for sequential access + // Align each hash table to page boundary + char* dbtable = serial_mmap(fid, hashTableSize, 1, + align_up(get_serial_hashtable_offset()+x*hashTableSize, get_page_logn())); + if(madvise(dbtable, hashTableSize, MADV_SEQUENTIAL)<0) + error("could not advise hashtable memory","","madvise"); + + maxColCount=0; + minColCount=O2_SERIAL_MAX_COLS; + meanColCount=0; + colCountN=0; + pt=(SerialElementT*)dbtable; + for( y = 0 ; y < H::N ; y++ ){ + // Move disk pointer to beginning of row + pe=pt+y*lshHeader->numCols; + + colCount=0; + if(bucket* bPtr = h[x][y]) + if(merge) +#ifdef LSH_BLOCK_FULL_ROWS + serial_merge_hashtable_row_format1(pe, bPtr->next, colCount); // skip collision counter bucket + else + serial_write_hashtable_row_format1(pe, bPtr->next, colCount); // skip collision counter bucket +#else + serial_merge_hashtable_row_format1(pe, bPtr, colCount); + else + serial_write_hashtable_row_format1(pe, bPtr, colCount); +#endif + if(colCount){ + if(colCount<minColCount) + minColCount=colCount; + if(colCount>maxColCount) + maxColCount=colCount; + meanColCount+=colCount; + colCountN++; + } + } + if(colCountN) + std::cout << "#rows with collisions =" << colCountN << ", mean = " << meanColCount/(float)colCountN + << ", min = " << minColCount << ", max = " << maxColCount + << endl; + serial_munmap(dbtable, hashTableSize); + } + + // We're done writing + return 1; +} + +void G::serial_merge_hashtable_row_format1(SerialElementT* pr, bucket* b, Uns32T& colCount){ + while(b && b->t2!=IFLAG){ + SerialElementT*pe=pr; // reset disk pointer to beginning of row + serial_merge_element_format1(pe, b->snext, b->t2, colCount); + b=b->next; + } +} + +void G::serial_merge_element_format1(SerialElementT* pe, sbucket* sb, Uns32T t2, Uns32T& colCount){ + while(sb){ + if(colCount==lshHeader->numCols){ + std::cout << "!point-chain full " << endl; + return; + } + Uns32T c=0; + // Merge collision chains + while(c<lshHeader->numCols){ + if( (pe+c)->hashValue==IFLAG){ + (pe+c)->hashValue=t2; + (pe+c)->pointID=sb->pointID; + colCount=c+1; + if(c+1<lshHeader->numCols) + (pe+c+1)->hashValue=IFLAG; + break; + } + c++; + } + sb=sb->snext; + } + return; +} + +void G::serial_write_hashtable_row_format1(SerialElementT*& pe, bucket* b, Uns32T& colCount){ + pe->hashValue=IFLAG; + while(b && b->t2!=IFLAG){ + serial_write_element_format1(pe, b->snext, b->t2, colCount); + b=b->next; + } +} + +void G::serial_write_element_format1(SerialElementT*& pe, sbucket* sb, Uns32T t2, Uns32T& colCount){ + while(sb){ + if(colCount==lshHeader->numCols){ + std::cout << "!point-chain full " << endl; + return; + } + pe->hashValue=t2; + pe->pointID=sb->pointID; + pe++; + colCount++; + sb=sb->snext; + } + pe->hashValue=IFLAG; + return; +} + +int G::serialize_lsh_hashtables_format2(int fid, int merge){ + Uns32T x,y; + + if( merge && !serial_can_merge(O2_SERIAL_FILEFORMAT2) ) + error("Cannot merge core and serial LSH, data structure dimensions mismatch."); + + // We must pereform FORMAT1 merges in core + if(merge) + unserialize_lsh_hashtables_format2(fid); + + Uns32T colCount, meanColCount, colCountN, maxColCount, minColCount, t1; + lseek(fid, get_serial_hashtable_offset(), SEEK_SET); + + // Write the hash tables + for( x = 0 ; x < H::L ; x++ ){ + std::cout << (merge ? "merging":"writing") << " hash table " << x << " FORMAT2..."; + std::cout.flush(); + maxColCount=0; + minColCount=O2_SERIAL_MAX_COLS; + meanColCount=0; + colCountN=0; + for( y = 0 ; y < H::N ; y++ ){ + colCount=0; + if(bucket* bPtr = h[x][y]){ + t1 = y | O2_SERIAL_FLAGS_T1_BIT; + if( write(fid, &t1, sizeof(Uns32T)) != sizeof(Uns32T) ){ + close(fid); + error("write error in serial_write_hashtable_format2() [t1]"); + } +#ifdef LSH_BLOCK_FULL_ROWS + serial_write_hashtable_row_format2(fid, bPtr->next, colCount); // skip collision counter bucket +#else + serial_write_hashtable_row_format2(fid, bPtr, colCount); +#endif + } + if(colCount){ + if(colCount<minColCount) + minColCount=colCount; + if(colCount>maxColCount) + maxColCount=colCount; + meanColCount+=colCount; + colCountN++; + } + } + // Write END of table marker + t1 = O2_SERIAL_FLAGS_END_BIT; + if( write(fid, &t1, sizeof(Uns32T)) != sizeof(Uns32T) ){ + close(fid); + error("write error in serial_write_hashtable_format2() [end]"); + } + + if(colCountN) + std::cout << "#rows with collisions =" << colCountN << ", mean = " << meanColCount/(float)colCountN + << ", min = " << minColCount << ", max = " << maxColCount + << endl; + } + + // We're done writing + return 1; +} + +void G::serial_write_hashtable_row_format2(int fid, bucket* b, Uns32T& colCount){ + while(b && b->t2!=IFLAG){ + t2 = O2_SERIAL_FLAGS_T2_BIT; + if( write(fid, &t2, sizeof(Uns32T)) != sizeof(Uns32T) ){ + close(fid); + error("write error in serial_write_hashtable_row_format2()"); + } + t2 = b->t2; + if( write(fid, &t2, sizeof(Uns32T)) != sizeof(Uns32T) ){ + close(fid); + error("write error in serial_write_hashtable_row_format2()"); + } + serial_write_element_format2(fid, b->snext, colCount); + b=b->next; + } +} + +void G::serial_write_element_format2(int fid, sbucket* sb, Uns32T& colCount){ + while(sb){ + if(write(fid, &sb->pointID, sizeof(Uns32T))!=sizeof(Uns32T)){ + close(fid); + error("Write error in serial_write_element_format2()"); + } + colCount++; + sb=sb->snext; + } +} + + +int G::serial_create(char* filename, Uns32T FMT){ + return serial_create(filename, w, L, N, C, k, d, FMT); +} + + +int G::serial_create(char* filename, float binWidth, Uns32T numTables, Uns32T numRows, Uns32T numCols, + Uns32T numFuns, Uns32T dim, Uns32T FMT){ + + if(numTables > O2_SERIAL_MAX_TABLES || numRows > O2_SERIAL_MAX_ROWS + || numCols > O2_SERIAL_MAX_COLS || numFuns > O2_SERIAL_MAX_FUNS + || dim>O2_SERIAL_MAX_DIM){ + error("LSH parameters out of bounds for serialization"); + } + + int dbfid; + if ((dbfid = open (filename, O_RDWR|O_CREAT|O_EXCL, S_IRUSR|S_IWUSR|S_IRGRP|S_IWGRP|S_IROTH|S_IWOTH)) < 0) + error("Can't create serial file", filename, "open"); + get_lock(dbfid, 1); + + // Make header first to get size of serialized database + lshHeader = new SerialHeaderT(binWidth, numTables, numRows, numCols, numFuns, dim, radius, maxp, FMT); + + // go to the location corresponding to the last byte + if (lseek (dbfid, lshHeader->get_size() - 1, SEEK_SET) == -1) + error("lseek error in db file", "", "lseek"); + + // write a dummy byte at the last location + if (write (dbfid, "", 1) != 1) + error("write error", "", "write"); + + db = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 1); + + memcpy (db, lshHeader, O2_SERIAL_HEADER_SIZE); + + serial_munmap(db, O2_SERIAL_HEADER_SIZE); + + close(dbfid); + + std::cout << "done initializing tables." << endl; + + return 1; +} + +char* G::serial_mmap(int dbfid, Uns32T memSize, Uns32T forWrite, off_t offset){ + if(forWrite){ + if ((db = (char*) mmap(0, memSize, PROT_READ | PROT_WRITE, + MAP_SHARED, dbfid, offset)) == (caddr_t) -1) + error("mmap error in request for writable serialized database", "", "mmap"); + } + else if ((db = (char*) mmap(0, memSize, PROT_READ, MAP_SHARED, dbfid, offset)) == (caddr_t) -1) + error("mmap error in read-only serialized database", "", "mmap"); + + return db; +} + +SerialHeaderT* G::serial_get_header(char* db){ + lshHeader = new SerialHeaderT(); + memcpy((char*)lshHeader, db, sizeof(SerialHeaderT)); + + if(lshHeader->lshMagic!=O2_SERIAL_MAGIC) + error("Not an LSH database file"); + + return lshHeader; +} + +void G::serial_munmap(char* db, Uns32T N){ + munmap(db, N); +} + +int G::serial_open(char* filename, int writeFlag){ + int dbfid; + if(writeFlag){ + if ((dbfid = open (filename, O_RDWR)) < 0) + error("Can't open serial file for read/write", filename, "open"); + get_lock(dbfid, writeFlag); + } + else{ + if ((dbfid = open (filename, O_RDONLY)) < 0) + error("Can't open serial file for read", filename, "open"); + get_lock(dbfid, 0); + } + + return dbfid; +} + +void G::serial_close(int dbfid){ + + release_lock(dbfid); + close(dbfid); +} + +int G::unserialize_lsh_header(char* filename){ + + int dbfid; + char* db; + // Test to see if file exists + if((dbfid = open (filename, O_RDONLY)) < 0) + error("Can't open the file", filename, "open"); + close(dbfid); + dbfid = serial_open(filename, 0); // open for read + db = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 0);// get database pointer + serial_get_header(db); // read header + serial_munmap(db, O2_SERIAL_HEADER_SIZE); // drop mmap + + // Unserialize header parameters + H::L = lshHeader->numTables; + H::m = (Uns32T)( (1.0 + sqrt(1 + 8.0*(int)H::L)) / 2.0); + H::N = lshHeader->numRows; + H::C = lshHeader->numCols; + H::k = lshHeader->numFuns; + H::d = lshHeader->dataDim; + G::w = lshHeader->binWidth; + G::radius = lshHeader->radius; + G::maxp = lshHeader->maxp; + + return dbfid; +} + +// unserialize the LSH parameters +// we leave the LSH tree on disk as a flat file +// it is this flat file that we search by memory mapping +void G::unserialize_lsh_functions(int dbfid){ + Uns32T j, kk; + float* pf; + Uns32T* pu; + + // Load the hash functions into core + char* db = serial_mmap(dbfid, get_serial_hashtable_offset(), 0);// get database pointer again + +#ifdef USE_U_FUNCTIONS + G::A = new float**[ H::m ]; // m x k x d random projectors + G::b = new float*[ H::m ]; // m x k random biases +#else + G::A = new float**[ H::L ]; // m x k x d random projectors + G::b = new float*[ H::L ]; // m x k random biases +#endif + G::g = new Uns32T*[ H::L ]; // L x k random projections + assert(g&&A&&b); // failure + + pf = get_serial_hashfunction_base(db); + +#ifdef USE_U_FUNCTIONS + for( j = 0 ; j < H::m ; j++ ){ // L functions gj(v) + G::A[j] = new float*[ H::k/2 ]; // k x d 2-stable distribution coefficients + G::b[j] = new float[ H::k/2 ]; // bias + assert( G::A[j] && G::b[j] ); // failure + for( kk = 0 ; kk < H::k/2 ; kk++ ){ // Normally distributed hash functions +#else + for( j = 0 ; j < H::L ; j++ ){ // L functions gj(v) + G::A[j] = new float*[ H::k ]; // k x d 2-stable distribution coefficients + G::b[j] = new float[ H::k ]; // bias + assert( G::A[j] && G::b[j] ); // failure + for( kk = 0 ; kk < H::k ; kk++ ){ // Normally distributed hash functions +#endif + G::A[j][kk] = new float[ H::d ]; + assert( G::A[j][kk] ); // failure + for(Uns32T i = 0 ; i < H::d ; i++ ) + G::A[j][kk][i] = *pf++; // Normally distributed random vectors + } + } +#ifdef USE_U_FUNCTIONS + for( j = 0 ; j < H::m ; j++ ) // biases b + for( kk = 0 ; kk < H::k/2 ; kk++ ) +#else + for( j = 0 ; j < H::L ; j++ ) // biases b + for( kk = 0 ; kk < H::k ; kk++ ) +#endif + G::b[j][kk] = *pf++; + + for( j = 0 ; j < H::L ; j++){ // 32-bit hash values, gj(v)=[x0 x1 ... xk-1] xk \in Z + G::g[j] = new Uns32T[ H::k ]; + assert( G::g[j] ); + } + + + H::__initialize_data_structures(); + + pu = (Uns32T*)pf; + for( j = 0 ; j < H::L ; j++ ) // Z projectors r1 + for( kk = 0 ; kk < H::k ; kk++ ) + H::r1[j][kk] = *pu++; + + for( j = 0 ; j < H::L ; j++ ) // Z projectors r2 + for( kk = 0 ; kk < H::k ; kk++ ) + H::r2[j][kk] = *pu++; + + serial_munmap(db, get_serial_hashtable_offset()); +} + +void G::unserialize_lsh_hashtables_format1(int fid){ + SerialElementT *pe, *pt; + Uns32T x,y; + Uns32T hashTableSize=sizeof(SerialElementT)*lshHeader->numRows*lshHeader->numCols; + // Read the hash tables into core + for( x = 0 ; x < H::L ; x++ ){ + // memory map a single hash table + // Align each hash table to page boundary + char* dbtable = serial_mmap(fid, hashTableSize, 0, + align_up(get_serial_hashtable_offset()+x*hashTableSize, get_page_logn())); + if(madvise(dbtable, hashTableSize, MADV_SEQUENTIAL)<0) + error("could not advise hashtable memory","","madvise"); + pt=(SerialElementT*)dbtable; + for( y = 0 ; y < H::N ; y++ ){ + // Move disk pointer to beginning of row + pe=pt+y*lshHeader->numCols; + unserialize_hashtable_row_format1(pe, h[x]+y); +#ifdef __LSH_DUMP_CORE_TABLES__ + printf("S[%d,%d]", x, y); + serial_bucket_dump(pe); + printf("C[%d,%d]", x, y); + dump_hashtable_row(h[x][y]); +#endif + } + serial_munmap(dbtable, hashTableSize); + } +} + +void G::unserialize_hashtable_row_format1(SerialElementT* pe, bucket** b){ + Uns32T colCount = 0; + while(colCount!=lshHeader->numCols && pe->hashValue !=IFLAG){ + H::p = pe->pointID; // current point ID + t2 = pe->hashValue; + bucket_insert_point(b); + pe++; + colCount++; + } +} + +void G::unserialize_lsh_hashtables_format2(int fid){ + Uns32T x=0,y=0; + + // Seek to hashtable base offset + if(lseek(fid, get_serial_hashtable_offset(), SEEK_SET)!=get_serial_hashtable_offset()){ + close(fid); + error("Seek error in unserialize_lsh_hashtables_format2"); + } + + // Read the hash tables into core (structure is given in header) + while( x < H::L){ + if(read(fid, &(H::t1), sizeof(Uns32T))!=sizeof(Uns32T)){ + close(fid); + error("Read error","unserialize_lsh_hashtables_format2()"); + } + if(H::t1&O2_SERIAL_FLAGS_END_BIT) + x++; // End of table + else + while(y < H::N){ + // Read a row and move file pointer to beginning of next row or table + if(!(H::t1&O2_SERIAL_FLAGS_T1_BIT)){ + close(fid); + error("State matchine error t1","unserialize_lsh_hashtables_format2()"); + } + y = H::t1 ^ O2_SERIAL_FLAGS_T1_BIT; + if(y>=H::N){ + close(fid); + error("Unserialized hashtable row pointer out of range","unserialize_lsh_hashtables_format2()"); + } + Uns32T token = unserialize_hashtable_row_format2(fid, h[x]+y); + +#ifdef __LSH_DUMP_CORE_TABLES__ + printf("C[%d,%d]", x, y); + dump_hashtable_row(h[x][y]); +#endif + // Check that token is valid + if( !(token&O2_SERIAL_FLAGS_T1_BIT || token&O2_SERIAL_FLAGS_END_BIT) ){ + close(fid); + error("State machine error end of row/table", "unserialize_lsh_hashtables_format2()"); + } + // Check for end of table flag + if(token&O2_SERIAL_FLAGS_END_BIT){ + x++; + break; + } + // Check for new row flag + if(token&O2_SERIAL_FLAGS_T1_BIT) + H::t1 = token; + } + } +} + +Uns32T G::unserialize_hashtable_row_format2(int fid, bucket** b){ + bool pointFound = false; + if(read(fid, &(H::t2), sizeof(Uns32T)) != sizeof(Uns32T)){ + close(fid); + error("Read error T2 token","unserialize_hashtable_row_format2"); + } + if( !(H::t2==O2_SERIAL_FLAGS_END_BIT || H::t2==O2_SERIAL_FLAGS_T2_BIT)){ + close(fid); + error("State machine error: expected E or T2"); + } + while(!(H::t2==O2_SERIAL_FLAGS_END_BIT || H::t2&O2_SERIAL_FLAGS_T1_BIT)){ + pointFound=false; + // Check for T2 token + if(H::t2!=O2_SERIAL_FLAGS_T2_BIT) + error("State machine error T2 token", "unserialize_hashtable_row_format2()"); + // Read t2 value + if(read(fid, &(H::t2), sizeof(Uns32T)) != sizeof(Uns32T)){ + close(fid); + error("Read error t2","unserialize_hashtable_row_format2"); + } + if(read(fid, &(H::p), sizeof(Uns32T)) != sizeof(Uns32T)){ + close(fid); + error("Read error H::p","unserialize_hashtable_row_format2"); + } + while(!(H::p==O2_SERIAL_FLAGS_END_BIT || H::p&O2_SERIAL_FLAGS_T1_BIT || H::p==O2_SERIAL_FLAGS_T2_BIT )){ + pointFound=true; + bucket_insert_point(b); + if(read(fid, &(H::p), sizeof(Uns32T)) != sizeof(Uns32T)){ + close(fid); + error("Read error H::p","unserialize_hashtable_row_format2"); + } + } + H::t2 = H::p; // Copy last found token to t2 + if(!pointFound) + error("State machine error: point", "unserialize_hashtable_row_format2()"); + } + return H::t2; // holds current token +} + +void G::dump_hashtable_row(bucket* p){ + while(p && p->t2!=IFLAG){ + sbucket* sbp = p->snext; + while(sbp){ + printf("(%0X,%u)", p->t2, sbp->pointID); + fflush(stdout); + sbp=sbp->snext; + } + p=p->next; + } + printf("\n"); +} + + +// G::serial_retrieve_point( ... ) +// retrieves (pointID) from a serialized LSH database +// +// inputs: +// filename - file name of serialized LSH database +// vv - query point set +// +// outputs: +// inserts retrieved points into add_point() callback method +void G::serial_retrieve_point_set(char* filename, vector<vector<float> >& vv, ReporterCallbackPtr add_point, void* caller) +{ + int dbfid = serial_open(filename, 0); // open for read + char* dbheader = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 0);// get database pointer + serial_get_header(dbheader); // read header + serial_munmap(dbheader, O2_SERIAL_HEADER_SIZE); // drop header mmap + + if((lshHeader->flags & O2_SERIAL_FILEFORMAT2)){ + close(dbfid); + error("serial_retrieve_point_set is for SERIAL_FILEFORMAT1 only"); + } + + // size of each hash table + Uns32T hashTableSize=sizeof(SerialElementT)*lshHeader->numRows*lshHeader->numCols; + calling_instance = caller; // class instance variable used in ...bucket_chain_point() + add_point_callback = add_point; + + for(Uns32T j=0; j<L; j++){ + // memory map a single hash table for random access + char* db = serial_mmap(dbfid, hashTableSize, 0, + align_up(get_serial_hashtable_offset()+j*hashTableSize,get_page_logn())); + if(madvise(db, hashTableSize, MADV_RANDOM)<0) + error("could not advise local hashtable memory","","madvise"); + SerialElementT* pe = (SerialElementT*)db ; + for(Uns32T qpos=0; qpos<vv.size(); qpos++){ + compute_hash_functions(vv[qpos]); + __generate_hash_keys(*(g+j),*(r1+j),*(r2+j)); + serial_bucket_chain_point(pe+t1*lshHeader->numCols, qpos); // Point to correct row + } + serial_munmap(db, hashTableSize); // drop hashtable mmap + } + serial_close(dbfid); +} + +void G::serial_retrieve_point(char* filename, vector<float>& v, Uns32T qpos, ReporterCallbackPtr add_point, void* caller){ + int dbfid = serial_open(filename, 0); // open for read + char* dbheader = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 0);// get database pointer + serial_get_header(dbheader); // read header + serial_munmap(dbheader, O2_SERIAL_HEADER_SIZE); // drop header mmap + + if((lshHeader->flags & O2_SERIAL_FILEFORMAT2)){ + close(dbfid); + error("serial_retrieve_point is for SERIAL_FILEFORMAT1 only"); + } + + // size of each hash table + Uns32T hashTableSize=sizeof(SerialElementT)*lshHeader->numRows*lshHeader->numCols; + calling_instance = caller; + add_point_callback = add_point; + compute_hash_functions(v); + for(Uns32T j=0; j<L; j++){ + // memory map a single hash table for random access + char* db = serial_mmap(dbfid, hashTableSize, 0, + align_up(get_serial_hashtable_offset()+j*hashTableSize,get_page_logn())); + if(madvise(db, hashTableSize, MADV_RANDOM)<0) + error("could not advise local hashtable memory","","madvise"); + SerialElementT* pe = (SerialElementT*)db ; + __generate_hash_keys(*(g+j),*(r1+j),*(r2+j)); + serial_bucket_chain_point(pe+t1*lshHeader->numCols, qpos); // Point to correct row + serial_munmap(db, hashTableSize); // drop hashtable mmap + } + serial_close(dbfid); +} + +void G::serial_dump_tables(char* filename){ + int dbfid = serial_open(filename, 0); // open for read + char* dbheader = serial_mmap(dbfid, O2_SERIAL_HEADER_SIZE, 0);// get database pointer + serial_get_header(dbheader); // read header + serial_munmap(dbheader, O2_SERIAL_HEADER_SIZE); // drop header mmap + Uns32T hashTableSize=sizeof(SerialElementT)*lshHeader->numRows*lshHeader->numCols; + for(Uns32T j=0; j<L; j++){ + // memory map a single hash table for random access + char* db = serial_mmap(dbfid, hashTableSize, 0, + align_up(get_serial_hashtable_offset()+j*hashTableSize,get_page_logn())); + if(madvise(db, hashTableSize, MADV_SEQUENTIAL)<0) + error("could not advise local hashtable memory","","madvise"); + SerialElementT* pe = (SerialElementT*)db ; + printf("*********** TABLE %d ***************\n", j); + fflush(stdout); + int count=0; + do{ + printf("[%d,%d]", j, count++); + fflush(stdout); + serial_bucket_dump(pe); + pe+=lshHeader->numCols; + }while(pe<(SerialElementT*)db+lshHeader->numRows*lshHeader->numCols); + } + +} + +void G::serial_bucket_dump(SerialElementT* pe){ + SerialElementT* pend = pe+lshHeader->numCols; + while( !(pe->hashValue==IFLAG || pe==pend ) ){ + printf("(%0X,%u)",pe->hashValue,pe->pointID); + pe++; + } + printf("\n"); + fflush(stdout); +} + +void G::serial_bucket_chain_point(SerialElementT* pe, Uns32T qpos){ + SerialElementT* pend = pe+lshHeader->numCols; + while( !(pe->hashValue==IFLAG || pe==pend ) ){ + if(pe->hashValue==t2){ // new match + add_point_callback(calling_instance, pe->pointID, qpos, radius); + } + pe++; + } +} + +void G::bucket_chain_point(bucket* p, Uns32T qpos){ + if(!p || p->t2==IFLAG) + return; + if(p->t2==t2){ // match + sbucket_chain_point(p->snext, qpos); // add to reporter + } + if(p->next){ + bucket_chain_point(p->next, qpos); // recurse + } +} + +void G::sbucket_chain_point(sbucket* p, Uns32T qpos){ + add_point_callback(calling_instance, p->pointID, qpos, radius); + if(p->snext){ + sbucket_chain_point(p->snext, qpos); + } +} + +void G::get_lock(int fd, bool exclusive) { + struct flock lock; + int status; + lock.l_type = exclusive ? F_WRLCK : F_RDLCK; + lock.l_whence = SEEK_SET; + lock.l_start = 0; + lock.l_len = 0; /* "the whole file" */ + retry: + do { + status = fcntl(fd, F_SETLKW, &lock); + } while (status != 0 && errno == EINTR); + if (status) { + if (errno == EAGAIN) { + sleep(1); + goto retry; + } else { + error("fcntl lock error", "", "fcntl"); + } + } +} + +void G::release_lock(int fd) { + struct flock lock; + int status; + + lock.l_type = F_UNLCK; + lock.l_whence = SEEK_SET; + lock.l_start = 0; + lock.l_len = 0; + + status = fcntl(fd, F_SETLKW, &lock); + + if (status) + error("fcntl unlock error", "", "fcntl"); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lshlib.h Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,336 @@ +// lshlib.h - a library for locality sensitive hashtable insertion and retrieval +// +// Author: Michael Casey +// Copyright (c) 2008 Michael Casey, All Rights Reserved + +/* GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + See LICENSE.txt +*/ + +#ifndef __LSHLIB_H +#define __LSHLIB_H + +#include <vector> +#include <queue> +#include <stdio.h> +#include <stdlib.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <sys/mman.h> +#include <fcntl.h> +#include <string.h> +#include <iostream> +#include <fstream> +#include <math.h> +#include <sys/time.h> +#include <assert.h> +#include <float.h> +#include <signal.h> +#include <time.h> +#include <limits.h> +#include <errno.h> +#ifdef MT19937 +#include "mt19937/mt19937ar.h" +#endif + +#define IntT int +#define LongUns64T long long unsigned +#define Uns32T unsigned +#define Int32T int +#define BooleanT int +#define TRUE 1 +#define FALSE 0 + +// A big number (>> max # of points) +#define INDEX_START_EMPTY 1000000000U + +// 4294967291 = 2^32-5 +#define UH_PRIME_DEFAULT 4294967291U + +// 2^29 +#define MAX_HASH_RND 536870912U + +// 2^32-1 +#define TWO_TO_32_MINUS_1 4294967295U + +#define O2_SERIAL_VERSION 1 // Sync with SVN version +#define O2_SERIAL_HEADER_SIZE sizeof(SerialHeaderT) +#define O2_SERIAL_ELEMENT_SIZE sizeof(SerialElementT) +#define O2_SERIAL_MAX_TABLES (200) +#define O2_SERIAL_MAX_ROWS (1000000) +#define O2_SERIAL_MAX_COLS (1000) +#define O2_SERIAL_MAX_DIM (2000) +#define O2_SERIAL_MAX_FUNS (100) +#define O2_SERIAL_MAX_BINWIDTH (200) + +// Flags for Serial Header +#define O2_SERIAL_FILEFORMAT1 (0x1U) // Optimize for on-disk search +#define O2_SERIAL_FILEFORMAT2 (0x2U) // Optimize for in-core search + +// Flags for serialization fileformat2: use high 3 bits of Uns32T +#define O2_SERIAL_FLAGS_T1_BIT (0x80000000U) +#define O2_SERIAL_FLAGS_T2_BIT (0x40000000U) +#define O2_SERIAL_FLAGS_END_BIT (0x20000000U) + +unsigned align_up(unsigned x, unsigned w); + +#define O2_SERIAL_FUNCTIONS_SIZE (align_up(sizeof(float) * O2_SERIAL_MAX_TABLES * O2_SERIAL_MAX_FUNS * O2_SERIAL_MAX_DIM \ ++ sizeof(float) * O2_SERIAL_MAX_TABLES * O2_SERIAL_MAX_FUNS + \ ++ sizeof(Uns32T) * O2_SERIAL_MAX_TABLES * O2_SERIAL_MAX_FUNS * 2 \ ++ O2_SERIAL_HEADER_SIZE,get_page_logn())) + +#define O2_SERIAL_MAX_LSH_SIZE (O2_SERIAL_ELEMENT_SIZE * O2_SERIAL_MAX_TABLES \ + * O2_SERIAL_MAX_ROWS * O2_SERIAL_MAX_COLS + O2_SERIAL_FUNCTIONS_SIZE) + +#define O2_SERIAL_MAGIC ('o'|'2'<<8|'l'<<16|'s'<<24) + +using namespace std; + +Uns32T get_page_logn(); + +// Disk table entry +typedef class SerialElement SerialElementT; +class SerialElement { +public: + Uns32T hashValue; + Uns32T pointID; + + SerialElement(Uns32T h, Uns32T pID): + hashValue(h), + pointID(pID){} +}; + +// Disk header +typedef class SerialHeader SerialHeaderT; +class SerialHeader { +public: + Uns32T lshMagic; // unique identifier for file header + float binWidth; // hash-function bin width + Uns32T numTables; // number of hash tables in file + Uns32T numRows; // size of each hash table + Uns32T numCols; // max collisions in each hash table + Uns32T elementSize; // size of a hash bucket + Uns32T version; // version number of file format + Uns32T size; // total size of database (bytes) + Uns32T flags; // 32 bits of useful information + Uns32T dataDim; // vector dimensionality + Uns32T numFuns; // number of independent hash functions + float radius; // 32-bit floating point radius + Uns32T maxp; // number of unique IDs in the database + Uns32T value_14; // spare value + Uns32T value_15; // spare value + Uns32T value_16; // spare value + + SerialHeader(); + SerialHeader(float W, Uns32T L, Uns32T N, Uns32T C, Uns32T k, Uns32T d, float radius, Uns32T p, Uns32T FMT); + + float get_binWidth(){return binWidth;} + Uns32T get_numTables(){return numTables;} + Uns32T get_numRows(){return numRows;} + Uns32T get_numCols(){return numCols;} + Uns32T get_elementSize(){return elementSize;} + Uns32T get_version(){return version;} + Uns32T get_flags(){return flags;} + Uns32T get_size(){return size;} + Uns32T get_dataDim(){return dataDim;} + Uns32T get_numFuns(){return numFuns;} + Uns32T get_maxp(){return maxp;} +}; + +#define IFLAG 0xFFFFFFFF + +// Point-set collision bucket (sbucket). +// sbuckets form a collision chain that identifies PointIDs falling in the same locale. +// sbuckets are chained from a bucket containing the collision list's t2 identifier +class sbucket { + friend class bucket; + friend class H; + friend class G; + + public: + class sbucket* snext; + unsigned int pointID; + + sbucket(){ + snext=0; + pointID=IFLAG; + } + ~sbucket(){delete snext;} + sbucket* get_snext(){return snext;} +}; + +// bucket structure for a linked list of locales that collide with the same hash value t1 +// different buckets represent different locales, collisions within a locale are chained +// in sbuckets +class bucket { + friend class H; + friend class G; + bucket* next; + sbucket* snext; + public: + unsigned int t2; + bucket(){ + next=0; + snext=0; + t2=IFLAG; + } + ~bucket(){delete next;delete snext;} + bucket* get_next(){return next;} +}; + + +// The hash_functions for locality-sensitive hashing +class H{ + friend class G; + private: + bucket*** h; // hash functions + Uns32T** r1; // random ints for hashing + Uns32T** r2; // random ints for hashing + Uns32T t1; // first hash table key + Uns32T t2; // second hash table key + Uns32T P; // hash table prime number + bool use_u_functions; // flag to optimize computation of hashes + Uns32T bucketCount; // count of number of point buckets allocated + Uns32T pointCount; // count of number of points inserted + + void __initialize_data_structures(); + void __bucket_insert_point(bucket*); + void __sbucket_insert_point(sbucket*); + Uns32T __computeProductModDefaultPrime(Uns32T*,Uns32T*,IntT); + Uns32T __randr(); + bucket** __get_bucket(int j); + void __generate_hash_keys(Uns32T*,Uns32T*,Uns32T*); + void error(const char* a, const char* b = "", const char *sysFunc = 0); + + public: + Uns32T N; // num rows per table + Uns32T C; // num collision per row + Uns32T k; // num projections per hash function + Uns32T m; // ~sqrt num hash tables + Uns32T L; // L = m*(m-1)/2, conversely, m = (1 + sqrt(1 + 8.0*L)) / 2.0 + Uns32T d; // dimensions + Uns32T p; // current point + + H(){;} + H(Uns32T k, Uns32T m, Uns32T d, Uns32T N, Uns32T C); + ~H(); + + + Uns32T bucket_insert_point(bucket**); + + Uns32T get_t1(){return t1;} + Uns32T get_t2(){return t2;} + +}; + + +typedef void (*ReporterCallbackPtr)(void* objPtr, Uns32T pointID, Uns32T queryIndex, float squaredDistance); + +// Interface for indexing and retrieval +class G: public H{ + private: + float *** A; // m x k x d random projectors from R^d->R^k + float ** b; // m x k uniform additive constants + Uns32T ** g; // L x k random hash projections \in Z^k + float w; // width of hash slots (relative to normalized feature space) + float radius;// scaling coefficient for data (1./radius) + vector<vector<Uns32T> > uu; // Storage for m patial hash evaluations + Uns32T maxp; // highest pointID stored in database + void* calling_instance; // store calling object instance for member-function callback + void (*add_point_callback)(void*, Uns32T, Uns32T, float); + + void initialize_partial_functions(); + + void get_lock(int fd, bool exclusive); + void release_lock(int fd); + int serial_create(char* filename, Uns32T FMT); + int serial_create(char* filename, float binWidth, Uns32T nTables, Uns32T nRows, Uns32T nCols, Uns32T k, Uns32T d, Uns32T FMT); + char* serial_mmap(int dbfid, Uns32T sz, Uns32T w, off_t offset = 0); + void serial_munmap(char* db, Uns32T N); + int serial_open(char* filename,int writeFlag); + void serial_close(int dbfid); + + // Function to write hashfunctions to disk + int serialize_lsh_hashfunctions(int fid); + + // Functions to write hashtables to disk in format1 (optimized for on-disk retrieval) + int serialize_lsh_hashtables_format1(int fid, int merge); + void serial_write_hashtable_row_format1(SerialElementT*& pe, bucket* h, Uns32T& colCount); + void serial_write_element_format1(SerialElementT*& pe, sbucket* sb, Uns32T t2, Uns32T& colCount); + void serial_merge_hashtable_row_format1(SerialElementT* pr, bucket* h, Uns32T& colCount); + void serial_merge_element_format1(SerialElementT* pe, sbucket* sb, Uns32T t2, Uns32T& colCount); + int serial_can_merge(Uns32T requestedFormat); // Test to see whether core and on-disk structures are compatible + + // Functions to write hashtables to disk in format2 (optimized for in-core retrieval) + int serialize_lsh_hashtables_format2(int fid, int merge); + void serial_write_hashtable_row_format2(int fid, bucket* h, Uns32T& colCount); + void serial_write_element_format2(int fid, sbucket* sb, Uns32T& colCount); + + // Functions to read serial header and hash functions (format1 and format2) + int unserialize_lsh_header(char* filename); // read lsh header from disk into core + void unserialize_lsh_functions(int fid); // read the lsh hash functions into core + + // Functions to read hashtables in format1 + void unserialize_lsh_hashtables_format1(int fid); // read FORMAT1 hash tables into core (disk format) + void unserialize_hashtable_row_format1(SerialElementT* pe, bucket** b); // read lsh hash table row into core + + // Functions to read hashtables in format2 + void unserialize_lsh_hashtables_format2(int fid); // read FORMAT2 hash tables into core (core format) + Uns32T unserialize_hashtable_row_format2(int fid, bucket** b); // read lsh hash table row into core + + // Helper functions + void serial_print_header(Uns32T requestedFormat); + float* get_serial_hashfunction_base(char* db); + SerialElementT* get_serial_hashtable_base(char* db); + Uns32T get_serial_hashtable_offset(); // Size of SerialHeader + HashFunctions + SerialHeaderT* serial_get_header(char* db); + SerialHeaderT* lshHeader; + + // Core Retrieval/Inspections Functions + void bucket_chain_point(bucket* p, Uns32T qpos); + void sbucket_chain_point(sbucket* p, Uns32T qpos); + void dump_hashtable_row(bucket* p); + + // Serial (Format 1) Retrieval/Inspection Functions + void serial_bucket_chain_point(SerialElementT* pe, Uns32T qpos); + void serial_bucket_dump(SerialElementT* pe); + + // Hash functions + void compute_hash_functions(vector<float>& v); + float randn(); + float ranf(); + + char* db; // pointer to serialized structure + + public: + G(char* lshFile, bool lshInCore = false); // unserialize constructor + G(float w, Uns32T k,Uns32T m, Uns32T d, Uns32T N, Uns32T C, float r); // core constructor + ~G(); + + Uns32T insert_point(vector<float>&, Uns32T pointID); + void insert_point_set(vector<vector<float> >& vv, Uns32T basePointID); + + // point retrieval from core + void retrieve_point(vector<float>& v, Uns32T qpos, ReporterCallbackPtr, void* me=NULL); + // point set retrieval from core + void retrieve_point_set(vector<vector<float> >& vv, ReporterCallbackPtr, void* me=NULL); + // serial point set retrieval + void serial_retrieve_point_set(char* filename, vector<vector<float> >& vv, ReporterCallbackPtr, void* me=NULL); + // serial point retrieval + void serial_retrieve_point(char* filename, vector<float>& vv, Uns32T qpos, ReporterCallbackPtr, void* me=NULL); + + void serialize(char* filename, Uns32T serialFormat = O2_SERIAL_FILEFORMAT1); // write hashfunctions and hashtables to disk + + SerialHeaderT* get_lshHeader(){return lshHeader;} + float get_radius(){return radius;} + Uns32T get_maxp(){return maxp;} + void serial_dump_tables(char* filename); + float get_mean_collision_rate(){ return (float) pointCount / bucketCount ; } +}; + +typedef class G LSH; + + + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mt19937/mt19937ar.c Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,175 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + Copyright (C) 2005, Mutsuo Saito, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +#include <stdio.h> +#include "mt19937ar.h" + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti<N; mti++) { + mt[mti] = + (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length) +{ + int i, j, k; + init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk<N-M;kk++) { + y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); + mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL]; + } + for (;kk<N-1;kk++) { + y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); + mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,0x7fffffff]-interval */ +long genrand_int31(void) +{ + return (long)(genrand_int32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +double genrand_real1(void) +{ + return genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +double genrand_real2(void) +{ + return genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +double genrand_real3(void) +{ + return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void) +{ + unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} +/* These real versions are due to Isaku Wada, 2002/01/09 added */
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mt19937/mt19937ar.h Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,72 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + Copyright (C) 2005, Mutsuo Saito + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s); + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length); + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void); + +/* generates a random number on [0,0x7fffffff]-interval */ +long genrand_int31(void); + +/* These real versions are due to Isaku Wada, 2002/01/09 added */ +/* generates a random number on [0,1]-real-interval */ +double genrand_real1(void); + +/* generates a random number on [0,1)-real-interval */ +double genrand_real2(void); + +/* generates a random number on (0,1)-real-interval */ +double genrand_real3(void); + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mt19937/mt19937ar.out Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,403 @@ +1000 outputs of genrand_int32() +1067595299 955945823 477289528 4107218783 4228976476 +3344332714 3355579695 227628506 810200273 2591290167 +2560260675 3242736208 646746669 1479517882 4245472273 +1143372638 3863670494 3221021970 1773610557 1138697238 +1421897700 1269916527 2859934041 1764463362 3874892047 +3965319921 72549643 2383988930 2600218693 3237492380 +2792901476 725331109 605841842 271258942 715137098 +3297999536 1322965544 4229579109 1395091102 3735697720 +2101727825 3730287744 2950434330 1661921839 2895579582 +2370511479 1004092106 2247096681 2111242379 3237345263 +4082424759 219785033 2454039889 3709582971 835606218 +2411949883 2735205030 756421180 2175209704 1873865952 +2762534237 4161807854 3351099340 181129879 3269891896 + 776029799 2218161979 3001745796 1866825872 2133627728 + 34862734 1191934573 3102311354 2916517763 1012402762 +2184831317 4257399449 2899497138 3818095062 3030756734 +1282161629 420003642 2326421477 2741455717 1278020671 +3744179621 271777016 2626330018 2560563991 3055977700 +4233527566 1228397661 3595579322 1077915006 2395931898 +1851927286 3013683506 1999971931 3006888962 1049781534 +1488758959 3491776230 104418065 2448267297 3075614115 +3872332600 891912190 3936547759 2269180963 2633455084 +1047636807 2604612377 2709305729 1952216715 207593580 +2849898034 670771757 2210471108 467711165 263046873 +3569667915 1042291111 3863517079 1464270005 2758321352 +3790799816 2301278724 3106281430 7974801 2792461636 + 555991332 621766759 1322453093 853629228 686962251 +1455120532 957753161 1802033300 1021534190 3486047311 +1902128914 3701138056 4176424663 1795608698 560858864 +3737752754 3141170998 1553553385 3367807274 711546358 +2475125503 262969859 251416325 2980076994 1806565895 + 969527843 3529327173 2736343040 2987196734 1649016367 +2206175811 3048174801 3662503553 3138851612 2660143804 +1663017612 1816683231 411916003 3887461314 2347044079 +1015311755 1203592432 2170947766 2569420716 813872093 +1105387678 1431142475 220570551 4243632715 4179591855 +2607469131 3090613241 282341803 1734241730 1391822177 +1001254810 827927915 1886687171 3935097347 2631788714 +3905163266 110554195 2447955646 3717202975 3304793075 +3739614479 3059127468 953919171 2590123714 1132511021 +3795593679 2788030429 982155079 3472349556 859942552 +2681007391 2299624053 647443547 233600422 608168955 +3689327453 1849778220 1608438222 3968158357 2692977776 +2851872572 246750393 3582818628 3329652309 4036366910 +1012970930 950780808 3959768744 2538550045 191422718 +2658142375 3276369011 2927737484 1234200027 1920815603 +3536074689 1535612501 2184142071 3276955054 428488088 +2378411984 4059769550 3913744741 2732139246 64369859 +3755670074 842839565 2819894466 2414718973 1010060670 +1839715346 2410311136 152774329 3485009480 4102101512 +2852724304 879944024 1785007662 2748284463 1354768064 +3267784736 2269127717 3001240761 3179796763 895723219 + 865924942 4291570937 89355264 1471026971 4114180745 +3201939751 2867476999 2460866060 3603874571 2238880432 +3308416168 2072246611 2755653839 3773737248 1709066580 +4282731467 2746170170 2832568330 433439009 3175778732 + 26248366 2551382801 183214346 3893339516 1928168445 +1337157619 3429096554 3275170900 1782047316 4264403756 +1876594403 4289659572 3223834894 1728705513 4068244734 +2867840287 1147798696 302879820 1730407747 1923824407 +1180597908 1569786639 198796327 560793173 2107345620 +2705990316 3448772106 3678374155 758635715 884524671 + 486356516 1774865603 3881226226 2635213607 1181121587 +1508809820 3178988241 1594193633 1235154121 326117244 +2304031425 937054774 2687415945 3192389340 2003740439 +1823766188 2759543402 10067710 1533252662 4132494984 + 82378136 420615890 3467563163 541562091 3535949864 +2277319197 3330822853 3215654174 4113831979 4204996991 +2162248333 3255093522 2219088909 2978279037 255818579 +2859348628 3097280311 2569721123 1861951120 2907080079 +2719467166 998319094 2521935127 2404125338 259456032 +2086860995 1839848496 1893547357 2527997525 1489393124 +2860855349 76448234 2264934035 744914583 2586791259 +1385380501 66529922 1819103258 1899300332 2098173828 +1793831094 276463159 360132945 4178212058 595015228 + 177071838 2800080290 1573557746 1548998935 378454223 +1460534296 1116274283 3112385063 3709761796 827999348 +3580042847 1913901014 614021289 4278528023 1905177404 + 45407939 3298183234 1184848810 3644926330 3923635459 +1627046213 3677876759 969772772 1160524753 1522441192 + 452369933 1527502551 832490847 1003299676 1071381111 +2891255476 973747308 4086897108 1847554542 3895651598 +2227820339 1621250941 2881344691 3583565821 3510404498 + 849362119 862871471 797858058 2867774932 2821282612 +3272403146 3997979905 209178708 1805135652 6783381 +2823361423 792580494 4263749770 776439581 3798193823 +2853444094 2729507474 1071873341 1329010206 1289336450 +3327680758 2011491779 80157208 922428856 1158943220 +1667230961 2461022820 2608845159 387516115 3345351910 +1495629111 4098154157 3156649613 3525698599 4134908037 + 446713264 2137537399 3617403512 813966752 1157943946 +3734692965 1680301658 3180398473 3509854711 2228114612 +1008102291 486805123 863791847 3189125290 1050308116 +3777341526 4291726501 844061465 1347461791 2826481581 + 745465012 2055805750 4260209475 2386693097 2980646741 + 447229436 2077782664 1232942813 4023002732 1399011509 +3140569849 2579909222 3794857471 900758066 2887199683 +1720257997 3367494931 2668921229 955539029 3818726432 +1105704962 3889207255 2277369307 2746484505 1761846513 +2413916784 2685127085 4240257943 1166726899 4215215715 +3082092067 3960461946 1663304043 2087473241 4162589986 +2507310778 1579665506 767234210 970676017 492207530 +1441679602 1314785090 3262202570 3417091742 1561989210 +3011406780 1146609202 3262321040 1374872171 1634688712 +1280458888 2230023982 419323804 3262899800 39783310 +1641619040 1700368658 2207946628 2571300939 2424079766 + 780290914 2715195096 3390957695 163151474 2309534542 +1860018424 555755123 280320104 1604831083 2713022383 +1728987441 3639955502 623065489 3828630947 4275479050 +3516347383 2343951195 2430677756 635534992 3868699749 + 808442435 3070644069 4282166003 2093181383 2023555632 +1568662086 3422372620 4134522350 3016979543 3259320234 +2888030729 3185253876 4258779643 1267304371 1022517473 + 815943045 929020012 2995251018 3371283296 3608029049 +2018485115 122123397 2810669150 1411365618 1238391329 +1186786476 3155969091 2242941310 1765554882 279121160 +4279838515 1641578514 3796324015 13351065 103516986 +1609694427 551411743 2493771609 1316337047 3932650856 +4189700203 463397996 2937735066 1855616529 2626847990 + 55091862 3823351211 753448970 4045045500 1274127772 +1124182256 92039808 2126345552 425973257 386287896 +2589870191 1987762798 4084826973 2172456685 3366583455 +3602966653 2378803535 2901764433 3716929006 3710159000 +2653449155 3469742630 3096444476 3932564653 2595257433 + 318974657 3146202484 853571438 144400272 3768408841 + 782634401 2161109003 570039522 1886241521 14249488 +2230804228 1604941699 3928713335 3921942509 2155806892 + 134366254 430507376 1924011722 276713377 196481886 +3614810992 1610021185 1785757066 851346168 3761148643 +2918835642 3364422385 3012284466 3735958851 2643153892 +3778608231 1164289832 205853021 2876112231 3503398282 +3078397001 3472037921 1748894853 2740861475 316056182 +1660426908 168885906 956005527 3984354789 566521563 +1001109523 1216710575 2952284757 3834433081 3842608301 +2467352408 3974441264 3256601745 1409353924 1329904859 +2307560293 3125217879 3622920184 3832785684 3882365951 +2308537115 2659155028 1450441945 3532257603 3186324194 +1225603425 1124246549 175808705 3009142319 2796710159 +3651990107 160762750 1902254979 1698648476 1134980669 + 497144426 3302689335 4057485630 3603530763 4087252587 + 427812652 286876201 823134128 1627554964 3745564327 +2589226092 4202024494 62878473 3275585894 3987124064 +2791777159 1916869511 2585861905 1375038919 1403421920 + 60249114 3811870450 3021498009 2612993202 528933105 +2757361321 3341402964 2621861700 273128190 4015252178 +3094781002 1621621288 2337611177 1796718448 1258965619 +4241913140 2138560392 3022190223 4174180924 450094611 +3274724580 617150026 2704660665 1469700689 1341616587 + 356715071 1188789960 2278869135 1766569160 2795896635 + 57824704 2893496380 1235723989 1630694347 3927960522 + 428891364 1814070806 2287999787 4125941184 3968103889 +3548724050 1025597707 1404281500 2002212197 92429143 +2313943944 2403086080 3006180634 3561981764 1671860914 +1768520622 1803542985 844848113 3006139921 1410888995 +1157749833 2125704913 1789979528 1799263423 741157179 +2405862309 767040434 2655241390 3663420179 2172009096 +2511931187 1680542666 231857466 1154981000 157168255 +1454112128 3505872099 1929775046 2309422350 2143329496 +2960716902 407610648 2938108129 2581749599 538837155 +2342628867 430543915 740188568 1937713272 3315215132 +2085587024 4030765687 766054429 3517641839 689721775 +1294158986 1753287754 4202601348 1974852792 33459103 +3568087535 3144677435 1686130825 4134943013 3005738435 +3599293386 426570142 754104406 3660892564 1964545167 + 829466833 821587464 1746693036 1006492428 1595312919 +1256599985 1024482560 1897312280 2902903201 691790057 +1037515867 3176831208 1968401055 2173506824 1089055278 +1748401123 2941380082 968412354 1818753861 2973200866 +3875951774 1119354008 3988604139 1647155589 2232450826 +3486058011 3655784043 3759258462 847163678 1082052057 + 989516446 2871541755 3196311070 3929963078 658187585 +3664944641 2175149170 2203709147 2756014689 2456473919 +3890267390 1293787864 2830347984 3059280931 4158802520 +1561677400 2586570938 783570352 1355506163 31495586 +3789437343 3340549429 2092501630 896419368 671715824 +3530450081 3603554138 1055991716 3442308219 1499434728 +3130288473 3639507000 17769680 2259741420 487032199 +4227143402 3693771256 1880482820 3924810796 381462353 +4017855991 2452034943 2736680833 2209866385 2128986379 + 437874044 595759426 641721026 1636065708 3899136933 + 629879088 3591174506 351984326 2638783544 2348444281 +2341604660 2123933692 143443325 1525942256 364660499 + 599149312 939093251 1523003209 106601097 376589484 +1346282236 1297387043 764598052 3741218111 933457002 +1886424424 3219631016 525405256 3014235619 323149677 +2038881721 4100129043 2851715101 2984028078 1888574695 +2014194741 3515193880 4180573530 3461824363 2641995497 +3179230245 2902294983 2217320456 4040852155 1784656905 +3311906931 87498458 2752971818 2635474297 2831215366 +3682231106 2920043893 3772929704 2816374944 309949752 +2383758854 154870719 385111597 1191604312 1840700563 + 872191186 2925548701 1310412747 2102066999 1504727249 +3574298750 1191230036 3330575266 3180292097 3539347721 + 681369118 3305125752 3648233597 950049240 4173257693 +1760124957 512151405 681175196 580563018 1169662867 +4015033554 2687781101 699691603 2673494188 1137221356 + 123599888 472658308 1053598179 1012713758 3481064843 +3759461013 3981457956 3830587662 1877191791 3650996736 + 988064871 3515461600 4089077232 2225147448 1249609188 +2643151863 3896204135 2416995901 1397735321 3460025646 + +1000 outputs of genrand_real2() +0.76275443 0.99000644 0.98670464 0.10143112 0.27933125 +0.69867227 0.94218740 0.03427201 0.78842173 0.28180608 +0.92179002 0.20785655 0.54534773 0.69644020 0.38107718 +0.23978165 0.65286910 0.07514568 0.22765211 0.94872929 +0.74557914 0.62664415 0.54708246 0.90959343 0.42043116 +0.86334511 0.19189126 0.14718544 0.70259889 0.63426346 +0.77408121 0.04531601 0.04605807 0.88595519 0.69398270 +0.05377184 0.61711170 0.05565708 0.10133577 0.41500776 +0.91810699 0.22320679 0.23353705 0.92871862 0.98897234 +0.19786706 0.80558809 0.06961067 0.55840445 0.90479405 +0.63288060 0.95009721 0.54948447 0.20645042 0.45000959 +0.87050869 0.70806991 0.19406895 0.79286390 0.49332866 +0.78483914 0.75145146 0.12341941 0.42030252 0.16728160 +0.59906494 0.37575460 0.97815160 0.39815952 0.43595080 +0.04952478 0.33917805 0.76509902 0.61034321 0.90654701 +0.92915732 0.85365931 0.18812377 0.65913428 0.28814566 +0.59476081 0.27835931 0.60722542 0.68310435 0.69387186 +0.03699800 0.65897714 0.17527003 0.02889304 0.86777366 +0.12352068 0.91439461 0.32022990 0.44445731 0.34903686 +0.74639273 0.65918367 0.92492794 0.31872642 0.77749724 +0.85413832 0.76385624 0.32744211 0.91326300 0.27458185 +0.22190155 0.19865383 0.31227402 0.85321225 0.84243342 +0.78544200 0.71854080 0.92503892 0.82703064 0.88306297 +0.47284073 0.70059042 0.48003761 0.38671694 0.60465770 +0.41747204 0.47163243 0.72750808 0.65830223 0.10955369 +0.64215401 0.23456345 0.95944940 0.72822249 0.40888451 +0.69980355 0.26677428 0.57333635 0.39791582 0.85377858 +0.76962816 0.72004885 0.90903087 0.51376506 0.37732665 +0.12691640 0.71249738 0.81217908 0.37037313 0.32772374 +0.14238259 0.05614811 0.74363008 0.39773267 0.94859135 +0.31452454 0.11730313 0.62962618 0.33334237 0.45547255 +0.10089665 0.56550662 0.60539371 0.16027624 0.13245301 +0.60959939 0.04671662 0.99356286 0.57660859 0.40269560 +0.45274629 0.06699735 0.85064246 0.87742744 0.54508392 +0.87242982 0.29321385 0.67660627 0.68230715 0.79052073 +0.48592054 0.25186266 0.93769755 0.28565487 0.47219067 +0.99054882 0.13155240 0.47110470 0.98556600 0.84397623 +0.12875246 0.90953202 0.49129015 0.23792727 0.79481194 +0.44337770 0.96564297 0.67749118 0.55684872 0.27286897 +0.79538393 0.61965356 0.22487929 0.02226018 0.49248200 +0.42247006 0.91797788 0.99250134 0.23449967 0.52531508 +0.10246337 0.78685622 0.34310922 0.89892996 0.40454552 +0.68608407 0.30752487 0.83601319 0.54956031 0.63777550 +0.82199797 0.24890696 0.48801123 0.48661910 0.51223987 +0.32969635 0.31075073 0.21393155 0.73453207 0.15565705 +0.58584522 0.28976728 0.97621478 0.61498701 0.23891470 +0.28518540 0.46809591 0.18371914 0.37597910 0.13492176 +0.66849449 0.82811466 0.56240330 0.37548956 0.27562998 +0.27521910 0.74096121 0.77176757 0.13748143 0.99747138 +0.92504502 0.09175241 0.21389176 0.21766512 0.31183245 +0.23271221 0.21207367 0.57903312 0.77523344 0.13242613 +0.31037988 0.01204835 0.71652949 0.84487594 0.14982178 +0.57423142 0.45677888 0.48420169 0.53465428 0.52667473 +0.46880526 0.49849733 0.05670710 0.79022476 0.03872047 +0.21697212 0.20443086 0.28949326 0.81678186 0.87629474 +0.92297064 0.27373097 0.84625273 0.51505586 0.00582792 +0.33295971 0.91848412 0.92537226 0.91760033 0.07541125 +0.71745848 0.61158698 0.00941650 0.03135554 0.71527471 +0.24821915 0.63636652 0.86159918 0.26450229 0.60160194 +0.35557725 0.24477500 0.07186456 0.51757096 0.62120362 +0.97981062 0.69954667 0.21065616 0.13382753 0.27693186 +0.59644095 0.71500764 0.04110751 0.95730081 0.91600724 +0.47704678 0.26183479 0.34706971 0.07545431 0.29398385 +0.93236070 0.60486023 0.48015011 0.08870451 0.45548581 +0.91872718 0.38142712 0.10668643 0.01397541 0.04520355 +0.93822273 0.18011940 0.57577277 0.91427606 0.30911399 +0.95853475 0.23611214 0.69619891 0.69601980 0.76765372 +0.58515930 0.49479057 0.11288752 0.97187699 0.32095365 +0.57563608 0.40760618 0.78703383 0.43261152 0.90877651 +0.84686346 0.10599030 0.72872803 0.19315490 0.66152912 +0.10210518 0.06257876 0.47950688 0.47062066 0.72701157 +0.48915116 0.66110261 0.60170685 0.24516994 0.12726050 +0.03451185 0.90864994 0.83494878 0.94800035 0.91035206 +0.14480751 0.88458997 0.53498312 0.15963215 0.55378627 +0.35171349 0.28719791 0.09097957 0.00667896 0.32309622 +0.87561479 0.42534520 0.91748977 0.73908457 0.41793223 +0.99279792 0.87908370 0.28458072 0.59132853 0.98672190 +0.28547393 0.09452165 0.89910674 0.53681109 0.37931425 +0.62683489 0.56609740 0.24801549 0.52948179 0.98328855 +0.66403523 0.55523786 0.75886666 0.84784685 0.86829981 +0.71448906 0.84670080 0.43922919 0.20771016 0.64157936 +0.25664246 0.73055695 0.86395782 0.65852932 0.99061803 +0.40280575 0.39146298 0.07291005 0.97200603 0.20555729 +0.59616495 0.08138254 0.45796388 0.33681125 0.33989127 +0.18717090 0.53545811 0.60550838 0.86520709 0.34290701 +0.72743276 0.73023855 0.34195926 0.65019733 0.02765254 +0.72575740 0.32709576 0.03420866 0.26061893 0.56997511 +0.28439072 0.84422744 0.77637570 0.55982168 0.06720327 +0.58449067 0.71657369 0.15819609 0.58042821 0.07947911 +0.40193792 0.11376012 0.88762938 0.67532159 0.71223735 +0.27829114 0.04806073 0.21144026 0.58830274 0.04140071 +0.43215628 0.12952729 0.94668759 0.87391019 0.98382450 +0.27750768 0.90849647 0.90962737 0.59269720 0.96102026 +0.49544979 0.32007095 0.62585546 0.03119821 0.85953001 +0.22017528 0.05834068 0.80731217 0.53799961 0.74166948 +0.77426600 0.43938444 0.54862081 0.58575513 0.15886492 +0.73214332 0.11649057 0.77463977 0.85788827 0.17061997 +0.66838056 0.96076133 0.07949296 0.68521946 0.89986254 +0.05667410 0.12741385 0.83470977 0.63969104 0.46612929 +0.10200126 0.01194925 0.10476340 0.90285217 0.31221221 +0.32980614 0.46041971 0.52024973 0.05425470 0.28330912 +0.60426543 0.00598243 0.97244013 0.21135841 0.78561597 +0.78428734 0.63422849 0.32909934 0.44771136 0.27380750 +0.14966697 0.18156268 0.65686758 0.28726350 0.97074787 +0.63676171 0.96649494 0.24526295 0.08297372 0.54257548 +0.03166785 0.33735355 0.15946671 0.02102971 0.46228045 +0.11892296 0.33408336 0.29875681 0.29847692 0.73767569 +0.02080745 0.62980060 0.08082293 0.22993106 0.25031439 +0.87787525 0.45150053 0.13673441 0.63407612 0.97907688 +0.52241942 0.50580158 0.06273902 0.05270283 0.77031811 +0.05113352 0.24393329 0.75036441 0.37436336 0.22877652 +0.59975358 0.85707591 0.88691457 0.85547165 0.36641027 +0.58720133 0.45462835 0.09243817 0.32981586 0.07820411 +0.25421519 0.36004706 0.60092307 0.46192412 0.36758683 +0.98424170 0.08019934 0.68594024 0.45826386 0.29962317 +0.79365413 0.89231296 0.49478547 0.87645944 0.23590734 +0.28106737 0.75026285 0.08136314 0.79582424 0.76010628 +0.82792971 0.27947652 0.72482861 0.82191216 0.46171689 +0.79189752 0.96043686 0.51609668 0.88995725 0.28998963 +0.55191845 0.03934737 0.83033700 0.49553013 0.98009549 +0.19017594 0.98347750 0.33452066 0.87144372 0.72106301 +0.71272114 0.71465963 0.88361677 0.85571283 0.73782329 +0.20920458 0.34855153 0.46766817 0.02780062 0.74898344 +0.03680650 0.44866557 0.77426312 0.91025891 0.25195236 +0.87319953 0.63265037 0.25552148 0.27422476 0.95217406 +0.39281839 0.66441573 0.09158900 0.94515992 0.07800798 +0.02507888 0.39901462 0.17382573 0.12141278 0.85502334 +0.19902911 0.02160210 0.44460522 0.14688742 0.68020336 +0.71323733 0.60922473 0.95400380 0.99611159 0.90897777 +0.41073520 0.66206647 0.32064685 0.62805003 0.50677209 +0.52690101 0.87473387 0.73918362 0.39826974 0.43683919 +0.80459118 0.32422684 0.01958019 0.95319576 0.98326137 +0.83931735 0.69060863 0.33671416 0.68062550 0.65152380 +0.33392969 0.03451730 0.95227244 0.68200635 0.85074171 +0.64721009 0.51234433 0.73402047 0.00969637 0.93835057 +0.80803854 0.31485260 0.20089527 0.01323282 0.59933780 +0.31584602 0.20209563 0.33754800 0.68604181 0.24443049 +0.19952227 0.78162632 0.10336988 0.11360736 0.23536740 +0.23262256 0.67803776 0.48749791 0.74658435 0.92156640 +0.56706407 0.36683221 0.99157136 0.23421374 0.45183767 +0.91609720 0.85573315 0.37706276 0.77042618 0.30891908 +0.40709595 0.06944866 0.61342849 0.88817388 0.58734506 +0.98711323 0.14744128 0.63242656 0.87704136 0.68347125 +0.84446569 0.43265239 0.25146321 0.04130111 0.34259839 +0.92697368 0.40878778 0.56990338 0.76204273 0.19820348 +0.66314909 0.02482844 0.06669207 0.50205581 0.26084093 +0.65139159 0.41650223 0.09733904 0.56344203 0.62651696 +0.67332139 0.58037374 0.47258086 0.21010758 0.05713135 +0.89390629 0.10781246 0.32037450 0.07628388 0.34227964 +0.42190597 0.58201860 0.77363549 0.49595133 0.86031236 +0.83906769 0.81098161 0.26694195 0.14215941 0.88210306 +0.53634237 0.12090720 0.82480459 0.75930318 0.31847147 +0.92768077 0.01037616 0.56201727 0.88107122 0.35925856 +0.85860762 0.61109408 0.70408301 0.58434977 0.92192494 +0.62667915 0.75988365 0.06858761 0.36156496 0.58057195 +0.13636150 0.57719713 0.59340255 0.63530602 0.22976282 +0.71915530 0.41162531 0.63979565 0.09931342 0.79344045 +0.10893790 0.84450224 0.23122236 0.99485593 0.73637397 +0.17276368 0.13357764 0.74965804 0.64991737 0.61990341 +0.41523170 0.05878239 0.05687301 0.05497131 0.42868366 +0.42571090 0.25810502 0.89642955 0.30439758 0.39310223 +0.11357431 0.04288255 0.23397550 0.11200634 0.85621396 +0.89733974 0.37508865 0.42077265 0.68597384 0.72781399 +0.19296476 0.61699087 0.31667128 0.67756410 0.00177323 +0.05725176 0.79474693 0.18885238 0.06724856 0.68193156 +0.42202167 0.22082041 0.28554673 0.64995708 0.87851940 +0.29124547 0.61009521 0.87374537 0.05743712 0.69902994 +0.81925115 0.45653873 0.37236821 0.31118709 0.52734307 +0.39672836 0.38185294 0.30163915 0.17374510 0.04913278 +0.90404879 0.25742801 0.58266467 0.97663209 0.79823377 +0.36437958 0.15206043 0.26529938 0.22690047 0.05839021 +0.84721160 0.18622435 0.37809403 0.55706977 0.49828704 +0.47659049 0.24289680 0.88477595 0.07807463 0.56245739 +0.73490635 0.21099431 0.13164942 0.75840044 0.66877037 +0.28988183 0.44046090 0.24967434 0.80048356 0.26029740 +0.30416821 0.64151867 0.52067892 0.12880774 0.85465381 +0.02690525 0.19149288 0.49630295 0.79682619 0.43566145 +0.00288078 0.81484193 0.03763639 0.68529083 0.01339574 +0.38405386 0.30537067 0.22994703 0.44000045 0.27217985 +0.53831243 0.02870435 0.86282045 0.61831306 0.09164956 +0.25609707 0.07445781 0.72185784 0.90058883 0.30070608 +0.94476583 0.56822213 0.21933909 0.96772793 0.80063440 +0.26307906 0.31183306 0.16501252 0.55436179 0.68562285 +0.23829083 0.86511559 0.57868991 0.81888344 0.20126869 +0.93172350 0.66028129 0.21786948 0.78515828 0.10262106 +0.35390326 0.79303876 0.63427924 0.90479631 0.31024934 +0.60635447 0.56198079 0.63573813 0.91854197 0.99701497 +0.83085849 0.31692291 0.01925964 0.97446405 0.98751283 +0.60944293 0.13751018 0.69519957 0.68956636 0.56969015 +0.46440193 0.88341765 0.36754434 0.89223647 0.39786427 +0.85055280 0.12749961 0.79452122 0.89449784 0.14567830 +0.45716830 0.74822309 0.28200437 0.42546044 0.17464886 +0.68308746 0.65496587 0.52935411 0.12736159 0.61523955 +0.81590528 0.63107864 0.39786553 0.20102294 0.53292914 +0.75485590 0.59847044 0.32861691 0.12125866 0.58917183 +0.07638293 0.86845380 0.29192617 0.03989733 0.52180460 +0.32503407 0.64071852 0.69516575 0.74254998 0.54587026 +0.48713246 0.32920155 0.08719954 0.63497059 0.54328459 +0.64178757 0.45583809 0.70694291 0.85212760 0.86074305 +0.33163422 0.85739792 0.59908488 0.74566046 0.72157152
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mt19937/mtTest.c Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,65 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + Copyright (C) 2005, Mutsuo Saito, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +#include <stdio.h> +#include "mt19937ar.h" + +int main(void) +{ + int i; + unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4; + init_by_array(init, length); + printf("1000 outputs of genrand_int32()\n"); + for (i=0; i<1000; i++) { + printf("%10lu ", genrand_int32()); + if (i%5==4) printf("\n"); + } + printf("\n1000 outputs of genrand_real2()\n"); + for (i=0; i<1000; i++) { + printf("%10.8f ", genrand_real2()); + if (i%5==4) printf("\n"); + } + return 0; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mt19937/readme-mt.txt Tue Jul 29 22:01:17 2008 +0000 @@ -0,0 +1,79 @@ +This is a Mersenne Twister pseudorandom number generator +with period 2^19937-1 with improved initialization scheme, +modified on 2002/1/26 by Takuji Nishimura and Makoto Matsumoto. +modified on 2005/4/26 by Mutsuo Saito + +Contents of this tar ball: +readme-mt.txt this file +mt19937ar.c the C source (ar: initialize by ARray) +mt19937ar.h the C header file for mt19937ar +mtTest.c the C test main program of mt19937ar.c +mt19937ar.out Test outputs of six types generators. 1000 for each + +1. Initialization + The initialization scheme for the previous versions of MT +(e.g. 1999/10/28 version or earlier) has a tiny problem, that +the most significant bits of the seed is not well reflected +to the state vector of MT. + +This version (2002/1/26) has two initialization schemes: +init_genrand(seed) and init_by_array(init_key, key_length). + +init_genrand(seed) initializes the state vector by using +one unsigned 32-bit integer "seed", which may be zero. + +init_by_array(init_key, key_length) initializes the state vector +by using an array init_key[] of unsigned 32-bit integers +of length key_kength. If key_length is smaller than 624, +then each array of 32-bit integers gives distinct initial +state vector. This is useful if you want a larger seed space +than 32-bit word. + +2. Generation +After initialization, the following type of pseudorandom numbers +are available. + +genrand_int32() generates unsigned 32-bit integers. +genrand_int31() generates unsigned 31-bit integers. +genrand_real1() generates uniform real in [0,1] (32-bit resolution). +genrand_real2() generates uniform real in [0,1) (32-bit resolution). +genrand_real3() generates uniform real in (0,1) (32-bit resolution). +genrand_res53() generates uniform real in [0,1) with 53-bit resolution. + +Note: the last five functions call the first one. +if you need more speed for these five functions, you may +suppress the function call by copying genrand_int32() and +replacing the last return(), following to these five functions. + +3. main() +main() is an example to initialize with an array of length 4, +then 1000 outputs of unsigned 32-bit integers, +then 1000 outputs of real [0,1) numbers. + +4. The outputs +The output of the mt19937ar.c is in the file mt19937ar.out. +If you revise or translate the code, check the output +by using this file. + +5. Cryptography +This generator is not cryptoraphically secure. +You need to use a one-way (or hash) function to obtain +a secure random sequence. + +6. Correspondence +See: +URL http://www.math.keio.ac.jp/matumoto/emt.html +email matumoto@math.keio.ac.jp, nisimura@sci.kj.yamagata-u.ac.jp + +7. Reference +M. Matsumoto and T. Nishimura, +"Mersenne Twister: A 623-Dimensionally Equidistributed Uniform +Pseudo-Random Number Generator", +ACM Transactions on Modeling and Computer Simulation, +Vol. 8, No. 1, January 1998, pp 3--30. + +------- +Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, +All rights reserved. +Copyright (C) 2005, Mutsuo Saito +All rights reserved.
--- a/query.cpp Tue Jul 22 20:09:31 2008 +0000 +++ b/query.cpp Tue Jul 29 22:01:17 2008 +0000 @@ -1,5 +1,4 @@ #include "audioDB.h" - #include "reporter.h" bool audioDB::powers_acceptable(double p1, double p2) { @@ -17,50 +16,88 @@ } void audioDB::query(const char* dbName, const char* inFile, adb__queryResponse *adbQueryResponse) { - initTables(dbName, inFile); - Reporter *r = 0; + // init database tables and dbH first + if(query_from_key) + initTables(dbName); + else + initTables(dbName, inFile); + + // keyKeyPos requires dbH to be initialized + if(query_from_key && (!key || (query_from_key_index = getKeyPos((char*)key))==O2_ERR_KEYNOTFOUND)) + error("Query key not found :",key); + switch (queryType) { case O2_POINT_QUERY: sequenceLength = 1; normalizedDistance = false; - r = new pointQueryReporter<std::greater < NNresult > >(pointNN); + reporter = new pointQueryReporter< std::greater < NNresult > >(pointNN); break; case O2_TRACK_QUERY: sequenceLength = 1; normalizedDistance = false; - r = new trackAveragingReporter<std::greater < NNresult > >(pointNN, trackNN, dbH->numFiles); + reporter = new trackAveragingReporter< std::greater< NNresult > >(pointNN, trackNN, dbH->numFiles); break; - case O2_SEQUENCE_QUERY: + case O2_SEQUENCE_QUERY: + if(no_unit_norming) + normalizedDistance = false; if(radius == 0) { - r = new trackAveragingReporter<std::less < NNresult > >(pointNN, trackNN, dbH->numFiles); + reporter = new trackAveragingReporter< std::less< NNresult > >(pointNN, trackNN, dbH->numFiles); } else { - r = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles); + if(index_exists(dbName, radius, sequenceLength)){ + char* indexName = index_get_name(dbName, radius, sequenceLength); + lsh = new LSH(indexName); + assert(lsh); + reporter = new trackSequenceQueryRadReporter(trackNN, index_to_trackID(lsh->get_maxp())+1); + delete lsh; + delete[] indexName; + } + else + reporter = new trackSequenceQueryRadReporter(trackNN, dbH->numFiles); } break; - case O2_N_SEQUENCE_QUERY : + case O2_N_SEQUENCE_QUERY: + if(no_unit_norming) + normalizedDistance = false; if(radius == 0) { - r = new trackSequenceQueryNNReporter<std::less < NNresult > >(pointNN, trackNN, dbH->numFiles); + reporter = new trackSequenceQueryNNReporter< std::less < NNresult > >(pointNN, trackNN, dbH->numFiles); } else { - r = new trackSequenceQueryRadNNReporter(pointNN,trackNN, dbH->numFiles); + if(index_exists(dbName, radius, sequenceLength)){ + char* indexName = index_get_name(dbName, radius, sequenceLength); + lsh = new LSH(indexName); + assert(lsh); + reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, index_to_trackID(lsh->get_maxp())+1); + delete lsh; + delete[] indexName; + } + else + reporter = new trackSequenceQueryRadNNReporter(pointNN,trackNN, dbH->numFiles); } break; case O2_ONE_TO_ONE_N_SEQUENCE_QUERY : if(radius == 0) { error("query-type not yet supported"); } else { - r = new trackSequenceQueryRadNNReporterOneToOne(pointNN,trackNN, dbH->numFiles); + reporter = new trackSequenceQueryRadNNReporterOneToOne(pointNN,trackNN, dbH->numFiles); } break; default: error("unrecognized queryType in query()"); } - query_loop(dbName, inFile, r); - r->report(fileTable, adbQueryResponse); - delete r; + + // Test for index (again) here + if(radius && index_exists(dbName, radius, sequenceLength)) + index_query_loop(dbName, query_from_key_index); + else + query_loop(dbName, query_from_key_index); + + reporter->report(fileTable, adbQueryResponse); } // return ordinal position of key in keyTable +// this should really be a STL hash map search unsigned audioDB::getKeyPos(char* key){ + if(!dbH) + error("dbH not initialized","getKeyPos"); for(unsigned k=0; k<dbH->numFiles; k++) if(strncmp(fileTable + k*O2_FILETABLE_ENTRY_SIZE, key, strlen(key))==0) return k; @@ -155,8 +192,8 @@ sp = DD[j]; spd = D[j+w] + w; k = trackTable[track] - w; - while(k--) - *sp++ += *spd++; + while(k--) + *sp++ += *spd++; } } } else { // HOP_SIZE != 1 @@ -211,7 +248,7 @@ // actual query point's information. -- CSR, 2007-12-05 void audioDB::set_up_query(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned *nvp) { *nvp = (statbuf.st_size - sizeof(int)) / (dbH->dim * sizeof(double)); - + if(!(dbH->flags & O2_FLAG_L2NORM)) { error("Database must be L2 normed for sequence query","use -L2NORM"); } @@ -285,6 +322,92 @@ } } +// 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(double **qp, double **vqp, double **qnp, double **vqnp, double **qpp, double **vqpp, double *mqdp, unsigned *nvp, Uns32T queryIndex) { + 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; + + *nvp = trackTable[queryIndex]; + if(*nvp < sequenceLength) { + error("Query shorter than requested sequence length", "maybe use -l"); + } + + VERB_LOG(1, "performing norms... "); + + // Read query feature vectors from database + *qp = NULL; + lseek(dbfid, dbH->dataOffset + trackOffsetTable[queryIndex] * sizeof(double), SEEK_SET); + size_t allocatedSize = 0; + read_data(queryIndex, qp, &allocatedSize); + // Consistency check on allocated memory and query feature size + if(*nvp*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(*qnp = new double[*nvp]); + memcpy(*qnp, l2normTable+trackIndexOffset, *nvp*sizeof(double)); + sequence_sum(*qnp, *nvp, sequenceLength); + sequence_sqrt(*qnp, *nvp, sequenceLength); + + if( usingPower ){ + // Copy Power partial-sum coefficients + assert(*qpp = new double[*nvp]); + memcpy(*qpp, powerTable+trackIndexOffset, *nvp*sizeof(double)); + sequence_sum(*qpp, *nvp, sequenceLength); + sequence_average(*qpp, *nvp, sequenceLength); + } + + if (usingTimes) { + unsigned int k; + *mqdp = 0.0; + double *querydurs = new double[*nvp]; + double *timesdata = new double[*nvp*2]; + assert(querydurs && timesdata); + memcpy(timesdata, timesTable+trackIndexOffset, *nvp*sizeof(double)); + for(k = 0; k < *nvp; k++) { + querydurs[k] = timesdata[2*k+1] - timesdata[2*k]; + *mqdp += querydurs[k]; + } + *mqdp /= k; + + VERB_LOG(1, "mean query file duration: %f\n", *mqdp); + + delete [] querydurs; + delete [] timesdata; + } + // Defaults, for exhaustive search (!usingQueryPoint) + *vqp = *qp; + *vqnp = *qnp; + *vqpp = *qpp; + + if(usingQueryPoint) { + if(queryPoint > *nvp || queryPoint > *nvp - sequenceLength + 1) { + error("queryPoint > numVectors-wL+1 in query"); + } else { + VERB_LOG(1, "query point: %u\n", queryPoint); + *vqp = *qp + queryPoint * dbH->dim; + *vqnp = *qnp + queryPoint; + if (usingPower) { + *vqpp = *qpp + queryPoint; + } + *nvp = sequenceLength; + } + } +} + + // 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. @@ -341,14 +464,98 @@ *vspp = *spp; } -void audioDB::query_loop(const char* dbName, const char* inFile, Reporter *reporter) { +// 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(double* query, double* qnPtr, double* qpPtr, double meanQdur, Uns32T numVectors){ + unsigned int dbVectors; + double *sNorm, *snPtr, *sPower = 0, *spPtr = 0; + double *meanDBdur = 0; + + // 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 + set_up_db(&sNorm, &snPtr, &sPower, &spPtr, &meanDBdur, &dbVectors); + + VERB_LOG(1, "matching points..."); + + assert(pointNN>0 && pointNN<=O2_MAXNN); + assert(trackNN>0 && trackNN<=O2_MAXNN); + + // We are guaranteed that the order of points is sorted by: + // qpos, trackID, spos + // 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 over the sequence + double dist; + size_t data_buffer_size = 0; + double *data_buffer = 0; + Uns32T trackOffset; + Uns32T trackIndexOffset; + 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(); + trackOffset=trackOffsetTable[pp.trackID]; // num data elements offset + trackIndexOffset=trackOffset/dbH->dim; // num vectors offset + if((!(usingPower) || powers_acceptable(qpPtr[usingQueryPoint?0:pp.qpos], sPower[trackIndexOffset+pp.spos])) && + ((usingQueryPoint?0:pp.qpos) < numVectors-sequenceLength+1 && pp.spos < trackTable[pp.trackID]-sequenceLength+1)){ + if(currentTrack!=pp.trackID){ + currentTrack=pp.trackID; + lseek(dbfid, dbH->dataOffset + trackOffset * sizeof(double), SEEK_SET); + read_data(currentTrack, &data_buffer, &data_buffer_size); + } + dist = dot_product_points(query+(usingQueryPoint?0:pp.qpos*dbH->dim), data_buffer+pp.spos*dbH->dim, dbH->dim*sequenceLength); + if(normalizedDistance) + dist = 2-(2/(qnPtr[usingQueryPoint?0:pp.qpos]*sNorm[trackIndexOffset+pp.spos]))*dist; + else + if(no_unit_norming) + dist = qnPtr[usingQueryPoint?0:pp.qpos]*qnPtr[usingQueryPoint?0:pp.qpos]+sNorm[trackIndexOffset+pp.spos]*sNorm[trackIndexOffset+pp.spos] - 2*dist; + // else + // dist = dist; + if((!radius) || dist <= (radius+O2_DISTANCE_TOLERANCE)) + reporter->add_point(pp.trackID, pp.qpos, pp.spos, dist); + } + exact_evaluation_queue->pop(); + } +} + +// A completely unprotected dot-product method +// Caller is responsible for ensuring that memory is within bounds +inline double audioDB::dot_product_points(double* q, double* p, Uns32T L){ + double dist = 0.0; + while(L--) + dist += *q++ * *p++; + return dist; +} + +void audioDB::query_loop(const char* dbName, Uns32T queryIndex) { unsigned int numVectors; double *query, *query_data; double *qNorm, *qnPtr, *qPower = 0, *qpPtr = 0; double meanQdur; - set_up_query(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors); + if(query_from_key) + set_up_query_from_key(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors, queryIndex); + else + set_up_query(&query_data, &query, &qNorm, &qnPtr, &qPower, &qpPtr, &meanQdur, &numVectors); unsigned int dbVectors; double *sNorm, *snPtr, *sPower = 0, *spPtr = 0; @@ -365,21 +572,12 @@ double **D = 0; // Differences query and target double **DD = 0; // Matched filter distance - D = new double*[numVectors]; + D = new double*[numVectors]; // pre-allocate DD = new double*[numVectors]; gettimeofday(&tv1, NULL); unsigned processedTracks = 0; - - // build track offset table - off_t *trackOffsetTable = new off_t[dbH->numFiles]; - unsigned cumTrack=0; off_t trackIndexOffset; - for(k = 0; k < dbH->numFiles; k++){ - trackOffsetTable[k] = cumTrack; - cumTrack += trackTable[k] * dbH->dim; - } - char nextKey[MAXSTR]; // Track loop @@ -403,6 +601,18 @@ } } + // 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; // numVectors offset read_data(track, &data_buffer, &data_buffer_size); @@ -425,15 +635,18 @@ for(j = 0; j <= numVectors - wL; j += HOP_SIZE) { for(k = 0; k <= trackTable[track] - wL; k += HOP_SIZE) { double thisDist; - if(normalizedDistance) { + if(normalizedDistance) thisDist = 2-(2/(qnPtr[j]*sNorm[trackIndexOffset+k]))*DD[j][k]; - } else { - thisDist = DD[j][k]; - } + else + if(no_unit_norming) + thisDist = qnPtr[j]*qnPtr[j]+sNorm[trackIndexOffset+k]*sNorm[trackIndexOffset+k] - 2*DD[j][k]; + else + thisDist = DD[j][k]; + // Power test if ((!usingPower) || powers_acceptable(qpPtr[j], sPower[trackIndexOffset + k])) { // radius test - if((!radius) || thisDist < radius) { + if((!radius) || thisDist <= (radius+O2_DISTANCE_TOLERANCE)) { reporter->add_point(track, usingQueryPoint ? queryPoint : j, k, thisDist); } } @@ -452,8 +665,6 @@ (tv1.tv_sec*1000 + tv1.tv_usec/1000)) // Clean up - if(trackOffsetTable) - delete[] trackOffsetTable; if(query_data) delete[] query_data; if(qNorm) @@ -493,3 +704,5 @@ } VERB_LOG(2, "done.\n"); } + +
--- a/reporter.h Tue Jul 22 20:09:31 2008 +0000 +++ b/reporter.h Tue Jul 29 22:01:17 2008 +0000 @@ -1,8 +1,13 @@ +#ifndef __REPORTER_H +#define __REPORTER_H + #include <utility> #include <queue> -#include <deque> #include <set> #include <functional> +#include <iostream> +#include "ReporterBase.h" +#include "audioDB.h" typedef struct nnresult { unsigned int trackID; @@ -32,7 +37,7 @@ return a.count > b.count; } -class Reporter { +class Reporter : public ReporterBase { public: virtual ~Reporter() {}; virtual void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) = 0; @@ -42,7 +47,7 @@ // adbQueryResponse thing everywhere; the fileTable argument is // there solely for convertion trackIDs into names. -- CSR, // 2007-12-10. - virtual void report(char *fileTable, adb__queryResponse *adbQueryResponse) = 0; + virtual void report(char *fileTable, void* adbQueryResponse) = 0; }; template <class T> class pointQueryReporter : public Reporter { @@ -50,7 +55,7 @@ pointQueryReporter(unsigned int pointNN); ~pointQueryReporter(); void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); + void report(char *fileTable, void* adbQueryResponse); private: unsigned int pointNN; std::priority_queue< NNresult, std::vector< NNresult >, T> *queue; @@ -79,7 +84,7 @@ } } -template <class T> void pointQueryReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) { +template <class T> void pointQueryReporter<T>::report(char *fileTable, void *adbQueryResponse) { NNresult r; std::vector<NNresult> v; unsigned int size = queue->size(); @@ -93,26 +98,32 @@ if(adbQueryResponse==0) { for(rit = v.rbegin(); rit < v.rend(); rit++) { r = *rit; - std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + if(fileTable) + std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << r.trackID << " "; std::cout << r.dist << " " << r.qpos << " " << r.spos << std::endl; } } else { - adbQueryResponse->result.__sizeRlist=size; - adbQueryResponse->result.__sizeDist=size; - adbQueryResponse->result.__sizeQpos=size; - adbQueryResponse->result.__sizeSpos=size; - adbQueryResponse->result.Rlist= new char*[size]; - adbQueryResponse->result.Dist = new double[size]; - adbQueryResponse->result.Qpos = new unsigned int[size]; - adbQueryResponse->result.Spos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeRlist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeDist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeQpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeSpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.Rlist= new char*[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist = new double[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Spos = new unsigned int[size]; unsigned int k = 0; for(rit = v.rbegin(); rit < v.rend(); rit++, k++) { r = *rit; - adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR]; - adbQueryResponse->result.Dist[k] = r.dist; - adbQueryResponse->result.Qpos[k] = r.qpos; - adbQueryResponse->result.Spos[k] = r.spos; - snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + ((adb__queryResponse*)adbQueryResponse)->result.Rlist[k] = new char[O2_MAXFILESTR]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist[k] = r.dist; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos[k] = r.qpos; + ((adb__queryResponse*)adbQueryResponse)->result.Spos[k] = r.spos; + if(fileTable) + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + else + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%d", r.trackID); } } } @@ -122,7 +133,7 @@ trackAveragingReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles); ~trackAveragingReporter(); void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); + void report(char *fileTable, void *adbQueryResponse); protected: unsigned int pointNN; unsigned int trackNN; @@ -153,7 +164,7 @@ } } -template <class T> void trackAveragingReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) { +template <class T> void trackAveragingReporter<T>::report(char *fileTable, void *adbQueryResponse) { std::priority_queue < NNresult, std::vector< NNresult>, T> result; for (int i = numFiles-1; i >= 0; i--) { unsigned int size = queues[i].size(); @@ -194,105 +205,36 @@ if(adbQueryResponse==0) { for(rit = v.rbegin(); rit < v.rend(); rit++) { r = *rit; - std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + if(fileTable) + std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << r.trackID << " "; std::cout << r.dist << " " << r.qpos << " " << r.spos << std::endl; } } else { - adbQueryResponse->result.__sizeRlist=size; - adbQueryResponse->result.__sizeDist=size; - adbQueryResponse->result.__sizeQpos=size; - adbQueryResponse->result.__sizeSpos=size; - adbQueryResponse->result.Rlist= new char*[size]; - adbQueryResponse->result.Dist = new double[size]; - adbQueryResponse->result.Qpos = new unsigned int[size]; - adbQueryResponse->result.Spos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeRlist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeDist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeQpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeSpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.Rlist= new char*[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist = new double[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Spos = new unsigned int[size]; unsigned int k = 0; for(rit = v.rbegin(); rit < v.rend(); rit++, k++) { r = *rit; - adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR]; - adbQueryResponse->result.Dist[k] = r.dist; - adbQueryResponse->result.Qpos[k] = r.qpos; - adbQueryResponse->result.Spos[k] = r.spos; - snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + ((adb__queryResponse*)adbQueryResponse)->result.Rlist[k] = new char[O2_MAXFILESTR]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist[k] = r.dist; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos[k] = r.qpos; + ((adb__queryResponse*)adbQueryResponse)->result.Spos[k] = r.spos; + if(fileTable) + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + else + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%d", r.trackID); } } } -// track Sequence Query Radius Reporter -// only return tracks and retrieved point counts -class trackSequenceQueryRadReporter : public Reporter { -public: - trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles); - ~trackSequenceQueryRadReporter(); - void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); - protected: - unsigned int trackNN; - unsigned int numFiles; - std::set<std::pair<unsigned int, unsigned int> > *set; - unsigned int *count; -}; - -trackSequenceQueryRadReporter::trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles): - trackNN(trackNN), numFiles(numFiles) { - set = new std::set<std::pair<unsigned int, unsigned int> >; - count = new unsigned int[numFiles]; - for (unsigned i = 0; i < numFiles; i++) { - count[i] = 0; - } -} - -trackSequenceQueryRadReporter::~trackSequenceQueryRadReporter() { - delete set; - delete [] count; -} - -void trackSequenceQueryRadReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) { - std::set<std::pair<unsigned int, unsigned int> >::iterator it; - std::pair<unsigned int, unsigned int> pair = std::make_pair(trackID, qpos); // only count this once - it = set->find(pair); - if (it == set->end()) { - set->insert(pair); - count[trackID]++; // only count if <tackID,qpos> pair is unique - } -} - -void trackSequenceQueryRadReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) { - std::priority_queue < Radresult, std::vector<Radresult>, std::greater<Radresult> > result; - // KLUDGE: doing this backwards in an attempt to get the same - // tiebreak behaviour as before. - for (int i = numFiles-1; i >= 0; i--) { - Radresult r; - r.trackID = i; - r.count = count[i]; - if(r.count > 0) { - result.push(r); - if (result.size() > trackNN) { - result.pop(); - } - } - } - - Radresult r; - std::vector<Radresult> v; - unsigned int size = result.size(); - for(unsigned int k = 0; k < size; k++) { - r = result.top(); - v.push_back(r); - result.pop(); - } - std::vector<Radresult>::reverse_iterator rit; - - if(adbQueryResponse==0) { - for(rit = v.rbegin(); rit < v.rend(); rit++) { - r = *rit; - std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.count << std::endl; - } - } else { - // FIXME - } -} - // Another type of trackAveragingReporter that reports all pointNN nearest neighbours template <class T> class trackSequenceQueryNNReporter : public trackAveragingReporter<T> { protected: @@ -302,15 +244,16 @@ using trackAveragingReporter<T>::pointNN; public: trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); + void report(char *fileTable, void *adbQueryResponse); }; template <class T> trackSequenceQueryNNReporter<T>::trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles) :trackAveragingReporter<T>(pointNN, trackNN, numFiles){} -template <class T> void trackSequenceQueryNNReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) { +template <class T> void trackSequenceQueryNNReporter<T>::report(char *fileTable, void *adbQueryResponse) { std::priority_queue < NNresult, std::vector< NNresult>, T> result; - std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles]; + std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues + = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles]; for (int i = numFiles-1; i >= 0; i--) { unsigned int size = queues[i].size(); @@ -353,7 +296,11 @@ if(adbQueryResponse==0) { for(rit = v.rbegin(); rit < v.rend(); rit++) { r = *rit; - std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.dist << std::endl; + if(fileTable) + std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << r.trackID << " "; + std::cout << r.dist << std::endl; unsigned int qsize = point_queues[r.trackID].size(); // Reverse the order of the points stored in point_queues for(unsigned int k=0; k < qsize; k++){ @@ -368,90 +315,126 @@ } } } else { - adbQueryResponse->result.__sizeRlist=size; - adbQueryResponse->result.__sizeDist=size; - adbQueryResponse->result.__sizeQpos=size; - adbQueryResponse->result.__sizeSpos=size; - adbQueryResponse->result.Rlist= new char*[size]; - adbQueryResponse->result.Dist = new double[size]; - adbQueryResponse->result.Qpos = new unsigned int[size]; - adbQueryResponse->result.Spos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeRlist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeDist=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeQpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.__sizeSpos=size; + ((adb__queryResponse*)adbQueryResponse)->result.Rlist= new char*[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist = new double[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos = new unsigned int[size]; + ((adb__queryResponse*)adbQueryResponse)->result.Spos = new unsigned int[size]; unsigned int k = 0; for(rit = v.rbegin(); rit < v.rend(); rit++, k++) { r = *rit; - adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR]; - adbQueryResponse->result.Dist[k] = r.dist; - adbQueryResponse->result.Qpos[k] = r.qpos; - adbQueryResponse->result.Spos[k] = r.spos; - snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + ((adb__queryResponse*)adbQueryResponse)->result.Rlist[k] = new char[O2_MAXFILESTR]; + ((adb__queryResponse*)adbQueryResponse)->result.Dist[k] = r.dist; + ((adb__queryResponse*)adbQueryResponse)->result.Qpos[k] = r.qpos; + ((adb__queryResponse*)adbQueryResponse)->result.Spos[k] = r.spos; + if(fileTable) + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE); + else + snprintf(((adb__queryResponse*)adbQueryResponse)->result.Rlist[k], O2_MAXFILESTR, "%d", r.trackID); } } // clean up delete[] point_queues; } +/********************** Radius Reporters **************************/ -// track Sequence Query Radius NN Reporter -// retrieve tracks ordered by query-point matches (one per track per query point) -// -// as well as sorted n-NN points per retrieved track -class trackSequenceQueryRadNNReporter : public Reporter { +class triple{ + public: + unsigned int a; + unsigned int b; + unsigned int c; + + triple(void); + triple(unsigned int, unsigned int, unsigned int); + unsigned int first(); + unsigned int second(); + unsigned int third(); +}; + +triple& make_triple(unsigned int, unsigned int, unsigned int); + +triple::triple(unsigned int a, unsigned int b, unsigned int c):a(a),b(b),c(c){} + +unsigned int triple::first(){return a;} +unsigned int triple::second(){return b;} +unsigned int triple::third(){return c;} + +triple::triple():a(0),b(0),c(0){} + +bool operator< (const triple &t1, const triple &t2) { + return ((t1.a < t2.a) || + ((t1.a == t2.a) && ((t1.b < t2.b) || + ((t1.b == t2.b) && (t1.c < t2.c))))); +} +bool operator== (const triple &t1, const triple &t2) { + return ((t1.a == t2.a) && (t1.b == t2.b) && (t1.c == t2.c)); +} + +triple& make_triple(unsigned int a, unsigned int b, unsigned int c){ + triple* t = new triple(a,b,c); + return *t; +} + + +// track Sequence Query Radius Reporter +// only return tracks and retrieved point counts +class trackSequenceQueryRadReporter : public Reporter { public: - trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles); - ~trackSequenceQueryRadNNReporter(); + trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles); + ~trackSequenceQueryRadReporter(); void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); + void report(char *fileTable, void *adbQueryResponse); protected: - unsigned int pointNN; unsigned int trackNN; unsigned int numFiles; - std::set< NNresult > *set; - std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues; + std::set<std::pair<unsigned int, unsigned int> > *set; + std::set< triple > *set_triple; unsigned int *count; }; -trackSequenceQueryRadNNReporter::trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles): -pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) { - // Where to count Radius track matches (one-to-one) - set = new std::set< NNresult >; - // Where to insert individual point matches (one-to-many) - point_queues = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles]; - +trackSequenceQueryRadReporter::trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles): + trackNN(trackNN), numFiles(numFiles) { + set = new std::set<std::pair<unsigned int, unsigned int> >; + set_triple = new std::set<triple, std::less<triple> >; count = new unsigned int[numFiles]; for (unsigned i = 0; i < numFiles; i++) { count[i] = 0; } } -trackSequenceQueryRadNNReporter::~trackSequenceQueryRadNNReporter() { +trackSequenceQueryRadReporter::~trackSequenceQueryRadReporter() { delete set; + delete set_triple; delete [] count; } -void trackSequenceQueryRadNNReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) { - std::set< NNresult >::iterator it; - NNresult r; - r.trackID = trackID; - r.qpos = qpos; - r.dist = dist; - r.spos = spos; +void trackSequenceQueryRadReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) { + std::set<std::pair<unsigned int, unsigned int> >::iterator it; + std::pair<unsigned int, unsigned int> pair = std::make_pair(trackID, qpos); // only count this once + std::set<triple>::iterator it2; + triple triple; - // Record all matching points (within radius) - if (!isnan(dist)) { - point_queues[trackID].push(r); - if(point_queues[trackID].size() > pointNN) - point_queues[trackID].pop(); - } + triple = make_triple(trackID, qpos, spos); // only count this once + + // Record unique <trackID,qpos,spos> triples (record one collision from all hash tables) + it2 = set_triple->find(triple); - // Record counts of <trackID,qpos> pairs - it = set->find(r); - if (it == set->end()) { - set->insert(r); - count[trackID]++; + if(it2 == set_triple->end()){ + set_triple->insert(triple); + + it = set->find(pair); + if (it == set->end()) { + set->insert(pair); + count[trackID]++; // only count if <tackID,qpos> pair is unique + } } } -void trackSequenceQueryRadNNReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) { +void trackSequenceQueryRadReporter::report(char *fileTable, void *adbQueryResponse) { std::priority_queue < Radresult, std::vector<Radresult>, std::greater<Radresult> > result; // KLUDGE: doing this backwards in an attempt to get the same // tiebreak behaviour as before. @@ -475,6 +458,120 @@ v.push_back(r); result.pop(); } + std::vector<Radresult>::reverse_iterator rit; + + if(adbQueryResponse==0) { + for(rit = v.rbegin(); rit < v.rend(); rit++) { + r = *rit; + if(fileTable) + std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << r.trackID << " "; + std::cout << r.count << std::endl; + } + } else { + // FIXME + } +} + +// track Sequence Query Radius NN Reporter +// retrieve tracks ordered by query-point matches (one per track per query point) +// +// as well as sorted n-NN points per retrieved track +class trackSequenceQueryRadNNReporter : public Reporter { +public: + trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles); + ~trackSequenceQueryRadNNReporter(); + void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); + void report(char *fileTable, void *adbQueryResponse); + protected: + unsigned int pointNN; + unsigned int trackNN; + unsigned int numFiles; + std::set<std::pair<unsigned int, unsigned int> > *set; + std::set< triple > *set_triple; + std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues; + unsigned int *count; +}; + +trackSequenceQueryRadNNReporter::trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles): +pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) { + // Where to count Radius track matches (one-to-one) + set = new std::set<std::pair<unsigned int, unsigned int> >; + set_triple = new std::set<triple, std::less<triple> >; + // Where to insert individual point matches (one-to-many) + point_queues = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles]; + + count = new unsigned int[numFiles]; + for (unsigned i = 0; i < numFiles; i++) { + count[i] = 0; + } +} + +trackSequenceQueryRadNNReporter::~trackSequenceQueryRadNNReporter() { + delete set; + delete set_triple; + delete [] count; +} + +void trackSequenceQueryRadNNReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) { + std::set<std::pair<unsigned int, unsigned int> >::iterator it; + std::set<triple>::iterator it2; + std::pair<unsigned int, unsigned int> pair; + triple triple; + NNresult r; + + pair = std::make_pair(trackID, qpos); // only count this once + triple = make_triple(trackID, qpos, spos); // only count this once + // Record unique <trackID,qpos,spos> triples (record one collision from all hash tables) + it2 = set_triple->find(triple); + if(it2 == set_triple->end()){ + set_triple->insert(triple); + // Record all matching points (within radius) + // Record counts of <trackID,qpos> pairs + it = set->find(pair); + if (it == set->end()) { + set->insert(pair); + count[trackID]++; + if (!isnan(dist)) { + r.trackID = trackID; + r.qpos = qpos; + r.dist = dist; + r.spos = spos; + point_queues[trackID].push(r); + if(point_queues[trackID].size() > pointNN) + point_queues[trackID].pop(); + } + } + } +} + +void trackSequenceQueryRadNNReporter::report(char *fileTable, void *adbQueryResponse) { + std::priority_queue < Radresult, std::vector<Radresult>, std::greater<Radresult> > result; + // KLUDGE: doing this backwards in an attempt to get the same + // tiebreak behaviour as before. + Radresult r; + NNresult rk; + + for (int i = numFiles-1; i >= 0; i--) { + r.trackID = i; + r.count = count[i]; + if(r.count > 0) { + cout.flush(); + result.push(r); + if (result.size() > trackNN) { + result.pop(); + } + } + } + + std::vector<Radresult> v; + unsigned int size = result.size(); + for(unsigned int k = 0; k < size; k++) { + r = result.top(); + v.push_back(r); + result.pop(); + } // Traverse tracks in descending order of count cardinality @@ -484,7 +581,11 @@ if(adbQueryResponse==0) { for(rit = v.rbegin(); rit < v.rend(); rit++) { r = *rit; - std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.count << std::endl; + if(fileTable) + std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << r.trackID << " "; + std::cout << r.count << std::endl; // Reverse the order of the points stored in point_queues unsigned int qsize=point_queues[r.trackID].size(); @@ -492,9 +593,8 @@ point_queue.push(point_queues[r.trackID].top()); point_queues[r.trackID].pop(); } - for(unsigned int k=0; k < qsize; k++){ - NNresult rk = point_queue.top(); + rk = point_queue.top(); std::cout << rk.dist << " " << rk.qpos << " " << rk.spos << std::endl; point_queue.pop(); } @@ -505,7 +605,6 @@ delete[] point_queues; } - /********** ONE-TO-ONE REPORTERS *****************/ // track Sequence Query Radius NN Reporter One-to-One @@ -516,7 +615,7 @@ trackSequenceQueryRadNNReporterOneToOne(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles); ~trackSequenceQueryRadNNReporterOneToOne(); void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist); - void report(char *fileTable, adb__queryResponse *adbQueryResponse); + void report(char *fileTable, void *adbQueryResponse); protected: unsigned int pointNN; unsigned int trackNN; @@ -564,7 +663,7 @@ } -void trackSequenceQueryRadNNReporterOneToOne::report(char *fileTable, adb__queryResponse *adbQueryResponse) { +void trackSequenceQueryRadNNReporterOneToOne::report(char *fileTable, void *adbQueryResponse) { if(adbQueryResponse==0) { std::vector< NNresult >::iterator vit; NNresult rk; @@ -572,11 +671,16 @@ rk = *vit; std::cout << rk.dist << " " << rk.qpos << " " - << rk.spos << " " - << fileTable + rk.trackID*O2_FILETABLE_ENTRY_SIZE - << std::endl; + << rk.spos << " "; + if(fileTable) + std::cout << fileTable + rk.trackID*O2_FILETABLE_ENTRY_SIZE << " "; + else + std::cout << rk.trackID << " "; + std::cout << std::endl; } } else { // FIXME } } + +#endif