mas01cr@243: #include <utility>
mas01cr@243: #include <queue>
mas01cr@277: #include <deque>
mas01cr@243: #include <set>
mas01cr@243: #include <functional>
mas01cr@243: 
mas01cr@243: typedef struct nnresult {
mas01cr@243:   unsigned int trackID;
mas01cr@243:   double dist;
mas01cr@243:   unsigned int qpos;
mas01cr@243:   unsigned int spos;
mas01cr@243: } NNresult;
mas01cr@243: 
mas01cr@243: typedef struct radresult {
mas01cr@243:   unsigned int trackID;
mas01cr@243:   unsigned int count;
mas01cr@243: } Radresult;
mas01cr@243: 
mas01cr@243: bool operator< (const NNresult &a, const NNresult &b) {
mas01cr@243:   return a.dist < b.dist;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: bool operator> (const NNresult &a, const NNresult &b) {
mas01cr@243:   return a.dist > b.dist;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: bool operator< (const Radresult &a, const Radresult &b) {
mas01cr@243:   return a.count < b.count;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: class Reporter {
mas01cr@243: public:
mas01cr@243:   virtual ~Reporter() {};
mas01cr@243:   virtual void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) = 0;
mas01cr@243:   // FIXME: this interface is a bit wacky: a relic of previous, more
mas01cr@243:   // confused times.  Really it might make sense to have separate
mas01cr@243:   // reporter classes for WS and for stdout, rather than passing this
mas01cr@243:   // adbQueryResponse thing everywhere; the fileTable argument is
mas01cr@243:   // there solely for convertion trackIDs into names.  -- CSR,
mas01cr@243:   // 2007-12-10.
mas01cr@243:   virtual void report(char *fileTable, adb__queryResponse *adbQueryResponse) = 0;
mas01cr@243: };
mas01cr@243: 
mas01cr@243: template <class T> class pointQueryReporter : public Reporter {
mas01cr@243: public:
mas01cr@243:   pointQueryReporter(unsigned int pointNN);
mas01cr@243:   ~pointQueryReporter();
mas01cr@243:   void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
mas01cr@243:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@243: private:
mas01cr@243:   unsigned int pointNN;
mas01cr@243:   std::priority_queue< NNresult, std::vector< NNresult >, T> *queue;
mas01cr@243: };
mas01cr@243: 
mas01cr@243: template <class T> pointQueryReporter<T>::pointQueryReporter(unsigned int pointNN)
mas01cr@243:   : pointNN(pointNN) {
mas01cr@243:   queue = new std::priority_queue< NNresult, std::vector< NNresult >, T>;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> pointQueryReporter<T>::~pointQueryReporter() {
mas01cr@243:   delete queue;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> void pointQueryReporter<T>::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
mas01cr@246:   if (!isnan(dist)) {
mas01cr@246:     NNresult r;
mas01cr@246:     r.trackID = trackID;
mas01cr@246:     r.qpos = qpos;
mas01cr@246:     r.spos = spos;
mas01cr@246:     r.dist = dist;
mas01cr@246:     queue->push(r);
mas01cr@246:     if(queue->size() > pointNN) {
mas01cr@246:       queue->pop();
mas01cr@246:     }
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> void pointQueryReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@243:   NNresult r;
mas01cr@243:   std::vector<NNresult> v;
mas01cr@243:   unsigned int size = queue->size();
mas01cr@243:   for(unsigned int k = 0; k < size; k++) {
mas01cr@243:     r = queue->top();
mas01cr@243:     v.push_back(r);
mas01cr@243:     queue->pop();
mas01cr@243:   }
mas01cr@243:   std::vector<NNresult>::reverse_iterator rit;
mas01cr@243:       
mas01cr@243:   if(adbQueryResponse==0) {
mas01cr@243:     for(rit = v.rbegin(); rit < v.rend(); rit++) {
mas01cr@243:       r = *rit;
mas01cr@277:       std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " ";
mas01cr@243:       std::cout << r.dist << " " << r.qpos << " " << r.spos << std::endl;
mas01cr@243:     }
mas01cr@243:   } else {
mas01cr@243:     adbQueryResponse->result.__sizeRlist=size;
mas01cr@243:     adbQueryResponse->result.__sizeDist=size;
mas01cr@243:     adbQueryResponse->result.__sizeQpos=size;
mas01cr@243:     adbQueryResponse->result.__sizeSpos=size;
mas01cr@243:     adbQueryResponse->result.Rlist= new char*[size];
mas01cr@243:     adbQueryResponse->result.Dist = new double[size];
mas01cr@243:     adbQueryResponse->result.Qpos = new unsigned int[size];
mas01cr@243:     adbQueryResponse->result.Spos = new unsigned int[size];
mas01cr@243:     unsigned int k = 0;
mas01cr@243:     for(rit = v.rbegin(); rit < v.rend(); rit++, k++) {
mas01cr@243:       r = *rit;
mas01cr@243:       adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR];
mas01cr@243:       adbQueryResponse->result.Dist[k] = r.dist;
mas01cr@243:       adbQueryResponse->result.Qpos[k] = r.qpos;
mas01cr@243:       adbQueryResponse->result.Spos[k] = r.spos;
mas01cr@277:       snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE);
mas01cr@243:     }
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> class trackAveragingReporter : public Reporter {
mas01cr@243:  public:
mas01cr@243:   trackAveragingReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles);
mas01cr@243:   ~trackAveragingReporter();
mas01cr@243:   void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
mas01cr@243:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@277:  protected:
mas01cr@243:   unsigned int pointNN;
mas01cr@243:   unsigned int trackNN;
mas01cr@243:   unsigned int numFiles;
mas01cr@243:   std::priority_queue< NNresult, std::vector< NNresult>, T > *queues;  
mas01cr@243: };
mas01cr@243: 
mas01cr@243: template <class T> trackAveragingReporter<T>::trackAveragingReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles) 
mas01cr@243:   : pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) {
mas01cr@243:   queues = new std::priority_queue< NNresult, std::vector< NNresult>, T >[numFiles];
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> trackAveragingReporter<T>::~trackAveragingReporter() {
mas01cr@243:   delete [] queues;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> void trackAveragingReporter<T>::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
mas01cr@246:   if (!isnan(dist)) {
mas01cr@246:     NNresult r;
mas01cr@246:     r.trackID = trackID;
mas01cr@246:     r.qpos = qpos;
mas01cr@246:     r.spos = spos;
mas01cr@246:     r.dist = dist;
mas01cr@246:     queues[trackID].push(r);
mas01cr@246:     if(queues[trackID].size() > pointNN) {
mas01cr@246:       queues[trackID].pop();
mas01cr@246:     }
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@243: template <class T> void trackAveragingReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@243:   std::priority_queue < NNresult, std::vector< NNresult>, T> result;
mas01cr@243:   for (int i = numFiles-1; i >= 0; i--) {
mas01cr@243:     unsigned int size = queues[i].size();
mas01cr@243:     if (size > 0) {
mas01cr@243:       NNresult r;
mas01cr@243:       double dist = 0;
mas01cr@243:       NNresult oldr = queues[i].top();
mas01cr@243:       for (unsigned int j = 0; j < size; j++) {
mas01cr@243:         r = queues[i].top();
mas01cr@243:         dist += r.dist;
mas01cr@243:         queues[i].pop();
mas01cr@243:         if (r.dist == oldr.dist) {
mas01cr@243:           r.qpos = oldr.qpos;
mas01cr@243:           r.spos = oldr.spos;
mas01cr@243:         } else {
mas01cr@243:           oldr = r;
mas01cr@243:         }
mas01cr@243:       }
mas01cr@243:       dist /= size;
mas01cr@243:       r.dist = dist; // trackID, qpos and spos are magically right already.
mas01cr@243:       result.push(r);
mas01cr@243:       if (result.size() > trackNN) {
mas01cr@243:         result.pop();
mas01cr@243:       }
mas01cr@243:     }
mas01cr@243:   }
mas01cr@243: 
mas01cr@243:   NNresult r;
mas01cr@243:   std::vector<NNresult> v;
mas01cr@243:   unsigned int size = result.size();
mas01cr@243:   for(unsigned int k = 0; k < size; k++) {
mas01cr@243:     r = result.top();
mas01cr@243:     v.push_back(r);
mas01cr@243:     result.pop();
mas01cr@243:   }
mas01cr@243:   std::vector<NNresult>::reverse_iterator rit;
mas01cr@243:       
mas01cr@243:   if(adbQueryResponse==0) {
mas01cr@243:     for(rit = v.rbegin(); rit < v.rend(); rit++) {
mas01cr@243:       r = *rit;
mas01cr@277:       std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " ";
mas01cr@243:       std::cout << r.dist << " " << r.qpos << " " << r.spos << std::endl;
mas01cr@243:     }
mas01cr@243:   } else {
mas01cr@243:     adbQueryResponse->result.__sizeRlist=size;
mas01cr@243:     adbQueryResponse->result.__sizeDist=size;
mas01cr@243:     adbQueryResponse->result.__sizeQpos=size;
mas01cr@243:     adbQueryResponse->result.__sizeSpos=size;
mas01cr@243:     adbQueryResponse->result.Rlist= new char*[size];
mas01cr@243:     adbQueryResponse->result.Dist = new double[size];
mas01cr@243:     adbQueryResponse->result.Qpos = new unsigned int[size];
mas01cr@243:     adbQueryResponse->result.Spos = new unsigned int[size];
mas01cr@243:     unsigned int k = 0;
mas01cr@243:     for(rit = v.rbegin(); rit < v.rend(); rit++, k++) {
mas01cr@243:       r = *rit;
mas01cr@243:       adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR];
mas01cr@243:       adbQueryResponse->result.Dist[k] = r.dist;
mas01cr@243:       adbQueryResponse->result.Qpos[k] = r.qpos;
mas01cr@243:       adbQueryResponse->result.Spos[k] = r.spos;
mas01cr@277:       snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE);
mas01cr@243:     }
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@277: // track Sequence Query Radius Reporter
mas01cr@277: // only return tracks and retrieved point counts
mas01cr@243: class trackSequenceQueryRadReporter : public Reporter { 
mas01cr@243: public:
mas01cr@243:   trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles);
mas01cr@243:   ~trackSequenceQueryRadReporter();
mas01cr@243:   void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
mas01cr@243:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@277:  protected:
mas01cr@243:   unsigned int trackNN;
mas01cr@243:   unsigned int numFiles;
mas01cr@243:   std::set<std::pair<unsigned int, unsigned int> > *set;
mas01cr@243:   unsigned int *count;
mas01cr@243: };
mas01cr@243: 
mas01cr@243: trackSequenceQueryRadReporter::trackSequenceQueryRadReporter(unsigned int trackNN, unsigned int numFiles):
mas01cr@243:   trackNN(trackNN), numFiles(numFiles) {
mas01cr@243:   set = new std::set<std::pair<unsigned int, unsigned int> >;
mas01cr@243:   count = new unsigned int[numFiles];
mas01cr@243:   for (unsigned i = 0; i < numFiles; i++) {
mas01cr@243:     count[i] = 0;
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@243: trackSequenceQueryRadReporter::~trackSequenceQueryRadReporter() {
mas01cr@243:   delete set;
mas01cr@243:   delete [] count;
mas01cr@243: }
mas01cr@243: 
mas01cr@243: void trackSequenceQueryRadReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
mas01cr@243:   std::set<std::pair<unsigned int, unsigned int> >::iterator it;
mas01cr@277:   std::pair<unsigned int, unsigned int> pair = std::make_pair(trackID, qpos); // only count this once
mas01cr@243:   it = set->find(pair);
mas01cr@243:   if (it == set->end()) {
mas01cr@243:     set->insert(pair);
mas01cr@277:     count[trackID]++; // only count if <tackID,qpos> pair is unique
mas01cr@243:   }
mas01cr@243: }
mas01cr@243: 
mas01cr@243: void trackSequenceQueryRadReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@243:   std::priority_queue < Radresult > result;
mas01cr@243:   // KLUDGE: doing this backwards in an attempt to get the same
mas01cr@243:   // tiebreak behaviour as before.
mas01cr@243:   for (int i = numFiles-1; i >= 0; i--) {
mas01cr@243:     Radresult r;
mas01cr@243:     r.trackID = i;
mas01cr@243:     r.count = count[i];
mas01cr@243:     if(r.count > 0) {
mas01cr@243:       result.push(r);
mas01cr@243:       if (result.size() > trackNN) {
mas01cr@243:         result.pop();
mas01cr@243:       }
mas01cr@243:     }
mas01cr@243:   }
mas01cr@243: 
mas01cr@243:   Radresult r;
mas01cr@243:   std::vector<Radresult> v;
mas01cr@243:   unsigned int size = result.size();
mas01cr@243:   for(unsigned int k = 0; k < size; k++) {
mas01cr@243:     r = result.top();
mas01cr@243:     v.push_back(r);
mas01cr@243:     result.pop();
mas01cr@243:   }
mas01cr@243:   std::vector<Radresult>::reverse_iterator rit;
mas01cr@243:       
mas01cr@243:   if(adbQueryResponse==0) {
mas01cr@243:     for(rit = v.rbegin(); rit < v.rend(); rit++) {
mas01cr@243:       r = *rit;
mas01cr@277:       std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.count << std::endl;
mas01cr@243:     }
mas01cr@243:   } else {
mas01cr@243:     // FIXME
mas01cr@243:   }
mas01cr@243: }
mas01cr@277: 
mas01cr@277: // Another type of trackAveragingReporter that reports all pointNN nearest neighbours
mas01cr@277: template <class T> class trackSequenceQueryNNReporter : public trackAveragingReporter<T> {
mas01cr@277:  protected:
mas01cr@277:   using trackAveragingReporter<T>::numFiles;
mas01cr@277:   using trackAveragingReporter<T>::queues;
mas01cr@277:   using trackAveragingReporter<T>::trackNN;
mas01cr@277:   using trackAveragingReporter<T>::pointNN;
mas01cr@277:  public:
mas01cr@277:   trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles);
mas01cr@277:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@277: };
mas01cr@277: 
mas01cr@277: template <class T> trackSequenceQueryNNReporter<T>::trackSequenceQueryNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles)
mas01cr@277: :trackAveragingReporter<T>(pointNN, trackNN, numFiles){}
mas01cr@277: 
mas01cr@277: template <class T> void trackSequenceQueryNNReporter<T>::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@277:   std::priority_queue < NNresult, std::vector< NNresult>, T> result;
mas01cr@277:   std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles];
mas01cr@277:   
mas01cr@277:   for (int i = numFiles-1; i >= 0; i--) {
mas01cr@277:     unsigned int size = queues[i].size();
mas01cr@277:     if (size > 0) {
mas01cr@277:       NNresult r;
mas01cr@277:       double dist = 0;
mas01cr@277:       NNresult oldr = queues[i].top();
mas01cr@277:       for (unsigned int j = 0; j < size; j++) {
mas01cr@277:         r = queues[i].top();
mas01cr@277:         dist += r.dist;
mas01cr@277: 	point_queues[i].push(r);
mas01cr@277: 	queues[i].pop();
mas01cr@277:         if (r.dist == oldr.dist) {
mas01cr@277:           r.qpos = oldr.qpos;
mas01cr@277:           r.spos = oldr.spos;
mas01cr@277:         } else {
mas01cr@277:           oldr = r;
mas01cr@277:         }
mas01cr@277:       }
mas01cr@277:       dist /= size;
mas01cr@277:       r.dist = dist; // trackID, qpos and spos are magically right already.
mas01cr@277:       result.push(r);
mas01cr@277:       if (result.size() > trackNN) {
mas01cr@277:         result.pop();
mas01cr@277:       }
mas01cr@277:     }
mas01cr@277:   }
mas01cr@277: 
mas01cr@277:   NNresult r;
mas01cr@277:   std::vector<NNresult> v;
mas01cr@277:   unsigned int size = result.size();
mas01cr@277:   for(unsigned int k = 0; k < size; k++) {
mas01cr@277:     r = result.top();
mas01cr@277:     v.push_back(r);
mas01cr@277:     result.pop();
mas01cr@277:   }
mas01cr@277:   std::vector<NNresult>::reverse_iterator rit;
mas01cr@277:   std::priority_queue< NNresult, std::vector< NNresult>, std::greater<NNresult> > point_queue;      
mas01cr@277: 
mas01cr@277:   if(adbQueryResponse==0) {
mas01cr@277:     for(rit = v.rbegin(); rit < v.rend(); rit++) {
mas01cr@277:       r = *rit;
mas01cr@277:       std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.dist << std::endl;
mas01cr@277:       unsigned int qsize = point_queues[r.trackID].size();
mas01cr@277:       // Reverse the order of the points stored in point_queues
mas01cr@277:       for(unsigned int k=0; k < qsize; k++){
mas01cr@277: 	point_queue.push( point_queues[r.trackID].top() );
mas01cr@277: 	point_queues[r.trackID].pop();
mas01cr@277:       }
mas01cr@277:       
mas01cr@277:       for(unsigned int k = 0; k < qsize; k++) {
mas01cr@277: 	NNresult rk = point_queue.top();
mas01cr@277: 	std::cout << rk.dist << " " << rk.qpos << " " << rk.spos << std::endl;
mas01cr@277: 	point_queue.pop();
mas01cr@277:       }
mas01cr@277:     }
mas01cr@277:   } else {
mas01cr@277:     adbQueryResponse->result.__sizeRlist=size;
mas01cr@277:     adbQueryResponse->result.__sizeDist=size;
mas01cr@277:     adbQueryResponse->result.__sizeQpos=size;
mas01cr@277:     adbQueryResponse->result.__sizeSpos=size;
mas01cr@277:     adbQueryResponse->result.Rlist= new char*[size];
mas01cr@277:     adbQueryResponse->result.Dist = new double[size];
mas01cr@277:     adbQueryResponse->result.Qpos = new unsigned int[size];
mas01cr@277:     adbQueryResponse->result.Spos = new unsigned int[size];
mas01cr@277:     unsigned int k = 0;
mas01cr@277:     for(rit = v.rbegin(); rit < v.rend(); rit++, k++) {
mas01cr@277:       r = *rit;
mas01cr@277:       adbQueryResponse->result.Rlist[k] = new char[O2_MAXFILESTR];
mas01cr@277:       adbQueryResponse->result.Dist[k] = r.dist;
mas01cr@277:       adbQueryResponse->result.Qpos[k] = r.qpos;
mas01cr@277:       adbQueryResponse->result.Spos[k] = r.spos;
mas01cr@277:       snprintf(adbQueryResponse->result.Rlist[k], O2_MAXFILESTR, "%s", fileTable+r.trackID*O2_FILETABLE_ENTRY_SIZE);
mas01cr@277:     }
mas01cr@277:   }
mas01cr@277:   // clean up
mas01cr@277:   delete[] point_queues;
mas01cr@277: }
mas01cr@277: 
mas01cr@277: 
mas01cr@277: // track Sequence Query Radius NN Reporter
mas01cr@277: // retrieve tracks ordered by query-point matches (one per track per query point)
mas01cr@277: //
mas01cr@277: // as well as sorted n-NN points per retrieved track
mas01cr@277: class trackSequenceQueryRadNNReporter : public Reporter { 
mas01cr@277: public:
mas01cr@277:   trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles);
mas01cr@277:   ~trackSequenceQueryRadNNReporter();
mas01cr@277:   void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
mas01cr@277:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@277:  protected:
mas01cr@277:   unsigned int pointNN;
mas01cr@277:   unsigned int trackNN;
mas01cr@277:   unsigned int numFiles;
mas01cr@277:   std::set< NNresult > *set;
mas01cr@277:   std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> > *point_queues;
mas01cr@277:   unsigned int *count;
mas01cr@277: };
mas01cr@277: 
mas01cr@277: trackSequenceQueryRadNNReporter::trackSequenceQueryRadNNReporter(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles):
mas01cr@277: pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) {
mas01cr@277:   // Where to count Radius track matches (one-to-one)
mas01cr@277:   set = new std::set< NNresult >; 
mas01cr@277:   // Where to insert individual point matches (one-to-many)
mas01cr@277:   point_queues = new std::priority_queue< NNresult, std::vector< NNresult>, std::less<NNresult> >[numFiles];
mas01cr@277:   
mas01cr@277:   count = new unsigned int[numFiles];
mas01cr@277:   for (unsigned i = 0; i < numFiles; i++) {
mas01cr@277:     count[i] = 0;
mas01cr@277:   }
mas01cr@277: }
mas01cr@277: 
mas01cr@277: trackSequenceQueryRadNNReporter::~trackSequenceQueryRadNNReporter() {
mas01cr@277:   delete set;
mas01cr@277:   delete [] count;
mas01cr@277: }
mas01cr@277: 
mas01cr@277: void trackSequenceQueryRadNNReporter::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
mas01cr@277:   std::set< NNresult >::iterator it;
mas01cr@277:   NNresult r;
mas01cr@277:   r.trackID = trackID;
mas01cr@277:   r.qpos = qpos;
mas01cr@277:   r.dist = dist;
mas01cr@277:   r.spos = spos;
mas01cr@277: 
mas01cr@277:   // Record all matching points (within radius)
mas01cr@277:   if (!isnan(dist)) {
mas01cr@277:     point_queues[trackID].push(r);
mas01cr@277:     if(point_queues[trackID].size() > pointNN)
mas01cr@277:       point_queues[trackID].pop();
mas01cr@277:   }
mas01cr@277: 
mas01cr@277:   // Record counts of <trackID,qpos> pairs
mas01cr@277:   it = set->find(r);
mas01cr@277:   if (it == set->end()) {
mas01cr@277:     set->insert(r);
mas01cr@277:     count[trackID]++;
mas01cr@277:   }
mas01cr@277: }
mas01cr@277: 
mas01cr@277: void trackSequenceQueryRadNNReporter::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@277:   std::priority_queue < Radresult > result;
mas01cr@277:   // KLUDGE: doing this backwards in an attempt to get the same
mas01cr@277:   // tiebreak behaviour as before.
mas01cr@277:   for (int i = numFiles-1; i >= 0; i--) {
mas01cr@277:     Radresult r;
mas01cr@277:     r.trackID = i;
mas01cr@277:     r.count = count[i];
mas01cr@277:     if(r.count > 0) {
mas01cr@277:       result.push(r);
mas01cr@277:       if (result.size() > trackNN) {
mas01cr@277:         result.pop();
mas01cr@277:       }
mas01cr@277:     }
mas01cr@277:   }
mas01cr@277: 
mas01cr@277:   Radresult r;
mas01cr@277:   std::vector<Radresult> v;
mas01cr@277:   unsigned int size = result.size();
mas01cr@277:   for(unsigned int k = 0; k < size; k++) {
mas01cr@277:     r = result.top();
mas01cr@277:     v.push_back(r);
mas01cr@277:     result.pop();
mas01cr@277:   }
mas01cr@277: 
mas01cr@277: 
mas01cr@277:   // Traverse tracks in descending order of count cardinality
mas01cr@277:   std::vector<Radresult>::reverse_iterator rit;
mas01cr@277:   std::priority_queue< NNresult, std::vector< NNresult>, std::greater<NNresult> > point_queue;
mas01cr@277: 
mas01cr@277:   if(adbQueryResponse==0) {
mas01cr@277:     for(rit = v.rbegin(); rit < v.rend(); rit++) {
mas01cr@277:       r = *rit;
mas01cr@277:       std::cout << fileTable + r.trackID*O2_FILETABLE_ENTRY_SIZE << " " << r.count << std::endl;
mas01cr@277: 
mas01cr@277:       // Reverse the order of the points stored in point_queues
mas01cr@277:       unsigned int qsize=point_queues[r.trackID].size();
mas01cr@277:       for(unsigned int k=0; k < qsize; k++){
mas01cr@277: 	point_queue.push(point_queues[r.trackID].top());
mas01cr@277: 	point_queues[r.trackID].pop();
mas01cr@277:       }
mas01cr@277: 
mas01cr@277:       for(unsigned int k=0; k < qsize; k++){
mas01cr@277: 	NNresult rk = point_queue.top();
mas01cr@277: 	std::cout << rk.dist << " " << rk.qpos << " " << rk.spos << std::endl;
mas01cr@277: 	point_queue.pop();
mas01cr@277:       }
mas01cr@277:     }
mas01cr@277:   } else {
mas01cr@277:     // FIXME
mas01cr@277:   }
mas01cr@277:   delete[] point_queues;
mas01cr@277: }
mas01cr@277: 
mas01cr@277: 
mas01cr@277: /********** ONE-TO-ONE REPORTERS *****************/
mas01cr@277: 
mas01cr@277: // track Sequence Query Radius NN Reporter One-to-One
mas01cr@277: // for each query point find the single best matching target point in all database
mas01cr@277: // report qpos, spos and trackID
mas01cr@277: class trackSequenceQueryRadNNReporterOneToOne : public Reporter { 
mas01cr@277: public:
mas01cr@277:   trackSequenceQueryRadNNReporterOneToOne(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles);
mas01cr@277:   ~trackSequenceQueryRadNNReporterOneToOne();
mas01cr@277:   void add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist);
mas01cr@277:   void report(char *fileTable, adb__queryResponse *adbQueryResponse);
mas01cr@277:  protected:
mas01cr@277:   unsigned int pointNN;
mas01cr@277:   unsigned int trackNN;
mas01cr@277:   unsigned int numFiles;
mas01cr@277:   std::set< NNresult > *set;
mas01cr@277:   std::vector< NNresult> *point_queue;
mas01cr@277:   unsigned int *count;
mas01cr@277: 
mas01cr@277: };
mas01cr@277: 
mas01cr@277: trackSequenceQueryRadNNReporterOneToOne::trackSequenceQueryRadNNReporterOneToOne(unsigned int pointNN, unsigned int trackNN, unsigned int numFiles):
mas01cr@277: pointNN(pointNN), trackNN(trackNN), numFiles(numFiles) {
mas01cr@277:   // Where to count Radius track matches (one-to-one)
mas01cr@277:   set = new std::set< NNresult >; 
mas01cr@277:   // Where to insert individual point matches (one-to-many)
mas01cr@277:   point_queue = new std::vector< NNresult >;
mas01cr@277:   
mas01cr@277:   count = new unsigned int[numFiles];
mas01cr@277:   for (unsigned i = 0; i < numFiles; i++) {
mas01cr@277:     count[i] = 0;
mas01cr@277:   }
mas01cr@277: }
mas01cr@277: 
mas01cr@277: trackSequenceQueryRadNNReporterOneToOne::~trackSequenceQueryRadNNReporterOneToOne() {
mas01cr@277:   delete set;
mas01cr@277:   delete [] count;
mas01cr@277: }
mas01cr@277: 
mas01cr@277: void trackSequenceQueryRadNNReporterOneToOne::add_point(unsigned int trackID, unsigned int qpos, unsigned int spos, double dist) {
mas01cr@277:   std::set< NNresult >::iterator it;
mas01cr@277:   NNresult r;
mas01cr@277: 
mas01cr@277:   r.qpos = qpos;
mas01cr@277:   r.trackID = trackID;
mas01cr@277:   r.spos = spos;
mas01cr@277:   r.dist = dist;
mas01cr@277: 
mas01cr@277:   if(point_queue->size() < r.qpos + 1){
mas01cr@277:     point_queue->resize( r.qpos + 1 );
mas01cr@277:     (*point_queue)[r.qpos].dist = 1e6;
mas01cr@277:   }
mas01cr@277: 
mas01cr@277:   if (r.dist < (*point_queue)[r.qpos].dist)
mas01cr@277:     (*point_queue)[r.qpos] = r;
mas01cr@277: 
mas01cr@277: }
mas01cr@277: 
mas01cr@277: void trackSequenceQueryRadNNReporterOneToOne::report(char *fileTable, adb__queryResponse *adbQueryResponse) {
mas01cr@277:   if(adbQueryResponse==0) {
mas01cr@277:     std::vector< NNresult >::iterator vit;
mas01cr@277:     NNresult rk;
mas01cr@277:     for( vit = point_queue->begin() ; vit < point_queue->end() ; vit++ ){
mas01cr@277:       rk = *vit;
mas01cr@277:       std::cout << rk.dist << " " 
mas01cr@277: 		<< rk.qpos << " " 
mas01cr@277: 		<< rk.spos << " " 
mas01cr@277: 		<< fileTable + rk.trackID*O2_FILETABLE_ENTRY_SIZE 
mas01cr@277: 		<< std::endl;
mas01cr@277:       }
mas01cr@277:   } else {
mas01cr@277:     // FIXME
mas01cr@277:   }
mas01cr@277: }