IndelQuery.h

   1 /**
   2  * k-errors queries
   3  *
   4  * Aligment using suffix filters.
   5  */
   6
   7 #ifndef _IndelQuery_H_
   8 #define _IndelQuery_H_
   9
  10 #include "ResultSet.h"
  11 #include "EditDistance.h"
  12 #include "Query.h"
  13
  14 #include <string>
  15 #include <cstdlib>
  16 #include <cstring>
  17
  18 namespace SXSI
  19 {
  20
  21 class IndelQuery : public Query
  22 {
  23 public:
  24     IndelQuery(TextCollection const *tc_)
  25         : Query(tc_), checked(0), splitPos(0), ed(0)
  26     {
  27         checked = new ResultSet(tc->getLength());
  28     }
  29
  30     virtual ~IndelQuery()
  31     {
  32         delete checked; checked = 0;
  33         delete [] splitPos; splitPos = 0;
  34         delete ed; ed = 0;
  35     }
  36
  37 protected:
  38     void firstStep(void)
  39     {
  40         checked->clear();
  41         constructPartition();
  42
  43         uchar *reversepat = new uchar[patlen+1];
  44         for (unsigned i = klimit+1; i > 1; --i)
  45         {
  46             if (debug) std::cout << "new try: " << splitPos[i] << ", split = " << i << std::endl;
  47
  48             for (unsigned j = 0; j < splitPos[i-1]; j++)
  49                 reversepat[j] = pat[splitPos[i-1]-j-1];
  50             reversepat[splitPos[i-1]] = 0;
  51             ed = new MyersEditDistanceIncremental(reversepat, splitPos[i-1], i-1);
  52
  53             exactMatch(i);
  54
  55             delete ed;
  56             ed = 0;
  57         }
  58         delete [] reversepat;
  59         reversepat = 0;
  60
  61         if (debug) std::cout << "new try: " << splitPos[1] << ", split = 1" << std::endl;
  62         ed = 0; // Not used, exact search only for the first block
  63         exactMatch(1);
  64     }
  65
  66 private:
  67
  68     void exactMatch(const unsigned split)
  69     {
  70         unsigned i = splitPos[split];
  71         const unsigned border = splitPos[split-1];
  72         if (debug) std::cerr << "Exact match: i = " << i << ", border = " << border << std::endl;
  73
  74         while (i != border  && pushChar(pat[i-1]))
  75             --i;
  76
  77         if (debug) std::cerr << "matched upto: i = " << i << std::endl;
  78         if (i == 0)
  79             checkMatch();
  80         else
  81             if (i == border)
  82                 nextStep(split);
  83
  84         // Rewind stack
  85         match.clear();
  86         while (!smin.empty()) smin.pop();
  87         while (!smax.empty()) smax.pop();
  88         smin.push(0);
  89         smax.push(textlen-1);
  90     }
  91
  92     void nextStep(const unsigned split)
  93     {
  94         // Maximum edit distance for the current filter is split-1
  95         if (ed->dist() <= split-1) checkMatch();
  96
  97         if (ed->end()) return;
  98
  99         // ALPHABET has been defined in Query
 100         for (const char *cp = ALPHABET; cp < ALPHABET + ALPHABET_SIZE; ++cp) {
 101             char c = *cp;
 102             if (!pushChar(c)) continue;
 103             const unsigned colmin = ed->pushChar(c);
 104
 105             // Edit distance is for reversed pattern, compute prefix lenght:
 106             const unsigned prefixl = ed->length - ed->getSuffixLength();
 107             // How many errors are we allowing on a prefix of given length?
 108             unsigned i = 1;
 109             while (splitPos[i] < prefixl)
 110                 ++i;
 111
 112             assert(i < klimit+1);
 113             assert(split-1 > i-1);
 114             const unsigned maxk = (split-1) - (i-1);
 115
 116             if (debug)
 117                 std::cerr << "c = " << c << ", split = " << split << ", maxk = " << maxk << ", colmin = " << colmin
 118                           << ", prefixl = " << prefixl << ", dist = " << ed->dist() << std::endl;
 119
 120             if (colmin <= maxk) nextStep(split);
 121
 122             //if (debug) std::cerr << "pop c = " << c << std::endl;
 123             ed->popChar();
 124             popChar();
 125         }
 126     }
 127
 128     inline void checkMatch(void)
 129     {
 130 //        if (debug) std::cerr << "min,max: " << smin.top() << ", " << smax.top() << std::endl;
 131         ulong min = smin.top();
 132         ulong max = smax.top();
 133         MyersEditDistance med(patlen, klimit);
 134         for (ulong i = min; i <= max; i++) {
 135             if (!checked->get(i))
 136             {
 137                 checked->set(i);
 138                 med.init((uchar const *)pat);
 139                 validateCandidate(med, i);
 140             }
 141         }
 142     }
 143
 144     void validateCandidate(MyersEditDistance &med, ulong sapos)
 145     {
 146         uchar *candidate = tc->getSuffix(sapos, patlen+klimit);
 147         unsigned l = strlen((char*)candidate);
 148         if (l < patlen)
 149         {
 150             delete [] candidate;
 151             return;
 152         }
 153
 154         std::pair<unsigned, unsigned> dist(0, patlen);
 155         if (klimit == 0)
 156         {  // MyersEditDistance class does not support klimit == 0
 157             unsigned i = 0;
 158             while (i < patlen && pat[i] == candidate[i])
 159                 ++i;
 160             if (i < patlen)
 161             {
 162                 delete [] candidate;
 163                 return; // There is at least one mismatch (and klimit == 0)
 164             }
 165         }
 166         else // Using prefixDist() when klimit > 0
 167             dist = med.prefixDist(candidate, l);
 168
 169         if (dist.first <= klimit)
 170             newMatch(sapos);
 171         delete [] candidate;
 172     }
 173
 174
 175     inline void newMatch(ulong sapos)
 176     {
 177         TextCollection::DocId di = tc->getDoc(sapos);
 178         result.insert(di);
 179     }
 180
 181
 182     void constructPartition(void)
 183     {
 184         delete [] splitPos;
 185         splitPos = new unsigned[klimit+2];
 186
 187         if (klimit == 0)
 188         {
 189             splitPos[0] = 0;
 190             splitPos[1]=patlen;
 191             return;
 192         }
 193
 194         // First piece size
 195         unsigned psize = patlen/(klimit+2);
 196         if (psize == 0)
 197         {
 198             std::cerr << "error: First psize == 0, for patlen = " << patlen << ", klimit = " << klimit << std::endl;
 199             std::exit(1);
 200         }
 201         splitPos[0] = 0;
 202         splitPos[1] = 2*psize;
 203
 204         // Last klimit are of size:
 205         psize = (patlen-splitPos[1])/klimit;
 206         if (psize == 0)
 207         {
 208             std::cerr << "error: psize == 0, for patlen = " << patlen << ", klimit = " << klimit << std::endl;
 209             std::exit(1);
 210         }
 211
 212         for (unsigned i = 2; i <= klimit; i++)
 213             splitPos[i]=splitPos[i-1]+psize;
 214         splitPos[klimit+1]=patlen;
 215
 216         if (splitPos[klimit] >= patlen)
 217         {
 218             std::cerr << "error: splitPos[" << klimit << "] == " << splitPos[klimit] << " >= patlen, for patlen = " << patlen << ", klimit = " << klimit << std::endl;
 219             std::exit(1);
 220         }
 221
 222         if (debug)
 223         {
 224             std::cerr << "splitpos: ";
 225             for (unsigned i= 0;i<=klimit+1; i++)
 226                 std::cerr << splitPos[i] << ", ";
 227             std::cerr << std::endl;
 228         }
 229     }
 230
 231
 232     ResultSet *checked;
 233     unsigned *splitPos;
 234     MyersEditDistanceIncremental *ed;
 235
 236     IndelQuery();
 237     // No copy constructor or assignment
 238     IndelQuery(IndelQuery const&);
 239     IndelQuery& operator = (IndelQuery const&);
 240 };
 241
 242 } // namespace
 243
 244 #endif // _IndelQuery_H_