Imported Debian patch 0.12.7-3
[bowtie.git] / diff_sample.h
1 #ifndef DIFF_SAMPLE_H_
2 #define DIFF_SAMPLE_H_
3
4 #include <stdint.h>
5 #include <seqan/sequence.h>
6 #include <seqan/index.h> // for LarssonSadakane
7 #include "assert_helpers.h"
8 #include "multikey_qsort.h"
9 #include "timer.h"
10 #include "auto_array.h"
11
12 using namespace std;
13 using namespace seqan;
14
15 #ifndef VMSG_NL
16 #define VMSG_NL(args...) \
17 if(this->verbose()) { \
18         stringstream tmp; \
19         tmp << args << endl; \
20         this->verbose(tmp.str()); \
21 }
22 #endif
23
24 #ifndef VMSG
25 #define VMSG(args...) \
26 if(this->verbose()) { \
27         stringstream tmp; \
28         tmp << args; \
29         this->verbose(tmp.str()); \
30 }
31 #endif
32
33 /**
34  * Routines for calculating, sanity-checking, and dispensing difference
35  * cover samples to clients.
36  */
37
38 /**
39  *
40  */
41 struct sampleEntry {
42         uint32_t maxV;
43         uint32_t numSamples;
44         uint32_t samples[128];
45 };
46
47 /// Array of Colbourn and Ling calculated difference covers up to
48 /// r = 16 (maxV = 5953)
49 static struct sampleEntry clDCs[16];
50 static bool clDCs_calced = false; /// have clDCs been calculated?
51
52 /**
53  * Check that the given difference cover 'ds' actually covers all
54  * differences for a periodicity of v.
55  */
56 template<typename T>
57 static bool dcRepOk(T v, String<T>& ds) {
58         // diffs[] records all the differences observed
59         bool *covered = new bool[v];
60         for(T i = 1; i < v; i++) {
61                 covered[i] = false;
62         }
63         for(T di = T(); di < length(ds); di++) {
64                 for(T dj = di+1; dj < length(ds); dj++) {
65                         assert_lt(ds[di], ds[dj]);
66                         T d1 = (ds[dj] - ds[di]);
67                         T d2 = (ds[di] + v - ds[dj]);
68                         assert_lt(d1, v);
69                         assert_lt(d2, v);
70                         covered[d1] = true;
71                         covered[d2] = true;
72                 }
73         }
74         bool ok = true;
75         for(T i = 1; i < v; i++) {
76                 if(covered[i] == false) {
77                         ok = false;
78                         break;
79                 }
80         }
81         delete[] covered;
82         return ok;
83 }
84
85 /**
86  * Return true iff each element of ts (with length 'limit') is greater
87  * than the last.
88  */
89 template<typename T>
90 static bool increasing(T* ts, size_t limit) {
91         for(size_t i = 0; i < limit-1; i++) {
92                 if(ts[i+1] <= ts[i]) return false;
93         }
94         return true;
95 }
96
97 /**
98  * Return true iff the given difference cover covers difference 'diff'
99  * mod 'v'.
100  */
101 template<typename T>
102 static inline bool hasDifference(T *ds, T d, T v, T diff) {
103         // diffs[] records all the differences observed
104         for(T di = T(); di < d; di++) {
105                 for(T dj = di+1; dj < d; dj++) {
106                         assert_lt(ds[di], ds[dj]);
107                         T d1 = (ds[dj] - ds[di]);
108                         T d2 = (ds[di] + v - ds[dj]);
109                         assert_lt(d1, v);
110                         assert_lt(d2, v);
111                         if(d1 == diff || d2 == diff) return true;
112                 }
113         }
114         return false;
115 }
116
117 /**
118  * Exhaustively calculate optimal difference cover samples for v = 4,
119  * 8, 16, 32, 64, 128, 256 and store results in p2DCs[]
120  */
121 template<typename T>
122 void calcExhaustiveDC(T i, bool verbose = false, bool sanityCheck = false) {
123         T v = i;
124         bool *diffs = new bool[v];
125         // v is the target period
126         T ld = (T)ceil(sqrt(v));
127         // ud is the upper bound on |D|
128         T ud = v / 2;
129         // for all possible |D|s
130         bool ok = true;
131         T *ds = NULL;
132         T d;
133         for(d = ld; d <= ud+1; d++) {
134                 // for all possible |D| samples
135                 ds = new T[d];
136                 for(T j = 0; j < d; j++) {
137                         ds[j] = j;
138                 }
139                 assert(increasing(ds, d));
140                 while(true) {
141                         // reset diffs[]
142                         for(T t = 1; t < v; t++) {
143                                 diffs[t] = false;
144                         }
145                         T diffCnt = 0;
146                         // diffs[] records all the differences observed
147                         for(T di = 0; di < d; di++) {
148                                 for(T dj = di+1; dj < d; dj++) {
149                                         assert_lt(ds[di], ds[dj]);
150                                         T d1 = (ds[dj] - ds[di]);
151                                         T d2 = (ds[di] + v - ds[dj]);
152                                         assert_lt(d1, v);
153                                         assert_lt(d2, v);
154                                         assert_gt(d1, 0);
155                                         assert_gt(d2, 0);
156                                         if(!diffs[d1]) diffCnt++; diffs[d1] = true;
157                                         if(!diffs[d2]) diffCnt++; diffs[d2] = true;
158                                 }
159                         }
160                         // Do we observe all possible differences (except 0)
161                         ok = diffCnt == v-1;
162                         if(ok) {
163                                 // Yes, all differences are covered
164                                 break;
165                         } else {
166                                 // Advance ds
167                                 // (Following is commented out because it turns out
168                                 // it's slow)
169                                 // Find a missing difference
170                                 //uint32_t missing = 0xffffffff;
171                                 //for(uint32_t t = 1; t < v; t++) {
172                                 //      if(diffs[t] == false) {
173                                 //              missing = diffs[t];
174                                 //              break;
175                                 //      }
176                                 //}
177                                 //assert_neq(missing, 0xffffffff);
178                                 assert(increasing(ds, d));
179                                 bool advanced = false;
180                                 bool keepGoing = false;
181                                 do {
182                                         keepGoing = false;
183                                         for(T bd = d-1; bd > 1; bd--) {
184                                                 T dif = (d-1)-bd;
185                                                 if(ds[bd] < v-1-dif) {
186                                                         ds[bd]++;
187                                                         assert_neq(0, ds[bd]);
188                                                         // Reset subsequent ones
189                                                         for(T bdi = bd+1; bdi < d; bdi++) {
190                                                                 assert_eq(0, ds[bdi]);
191                                                                 ds[bdi] = ds[bdi-1]+1;
192                                                                 assert_gt(ds[bdi], ds[bdi-1]);
193                                                         }
194                                                         assert(increasing(ds, d));
195                                                         // (Following is commented out because
196                                                         // it turns out it's slow)
197                                                         // See if the new DC has the missing value
198                                                         //if(!hasDifference(ds, d, v, missing)) {
199                                                         //      keepGoing = true;
200                                                         //      break;
201                                                         //}
202                                                         advanced = true;
203                                                         break;
204                                                 } else {
205                                                         ds[bd] = 0;
206                                                         // keep going
207                                                 }
208                                         }
209                                 } while(keepGoing);
210                                 // No solution for this |D|
211                                 if(!advanced) break;
212                                 assert(increasing(ds, d));
213                         }
214                 } // next sample assignment
215                 if(ok) {
216                         break;
217                 }
218                 delete[] ds;
219         } // next |D|
220         assert(ok);
221         delete[] diffs;
222         cout << "Did exhaustive v=" << v << " |D|=" << d << endl;
223         cout << "  ";
224         for(T i = 0; i < d; i++) {
225                 cout << ds[i];
226                 if(i < d-1) cout << ",";
227         }
228         cout << endl;
229         delete[] ds;
230 }
231
232 /**
233  * Routune for calculating the elements of clDCs up to r = 16 using the
234  * technique of Colbourn and Ling.
235  *
236  * See http://citeseer.ist.psu.edu/211575.html
237  */
238 template <typename T>
239 void calcColbournAndLingDCs(bool verbose = false, bool sanityCheck = false) {
240         for(T r = 0; r < 16; r++) {
241                 T maxv = 24*r*r + 36*r + 13; // Corollary 2.3
242                 T numsamp = 6*r + 4;
243                 clDCs[r].maxV = maxv;
244                 clDCs[r].numSamples = numsamp;
245                 memset(clDCs[r].samples, 0, 4 * 128);
246                 T i;
247                 // clDCs[r].samples[0] = 0;
248                 // Fill in the 1^r part of the B series
249                 for(i = 1; i < r+1; i++) {
250                         clDCs[r].samples[i] = clDCs[r].samples[i-1] + 1;
251                 }
252                 // Fill in the (r + 1)^1 part
253                 clDCs[r].samples[r+1] = clDCs[r].samples[r] + r + 1;
254                 // Fill in the (2r + 1)^r part
255                 for(i = r+2; i < r+2+r; i++) {
256                         clDCs[r].samples[i] = clDCs[r].samples[i-1] + 2*r + 1;
257                 }
258                 // Fill in the (4r + 3)^(2r + 1) part
259                 for(i = r+2+r; i < r+2+r+2*r+1; i++) {
260                         clDCs[r].samples[i] = clDCs[r].samples[i-1] + 4*r + 3;
261                 }
262                 // Fill in the (2r + 2)^(r + 1) part
263                 for(i = r+2+r+2*r+1; i < r+2+r+2*r+1+r+1; i++) {
264                         clDCs[r].samples[i] = clDCs[r].samples[i-1] + 2*r + 2;
265                 }
266                 // Fill in the last 1^r part
267                 for(i = r+2+r+2*r+1+r+1; i < r+2+r+2*r+1+r+1+r; i++) {
268                         clDCs[r].samples[i] = clDCs[r].samples[i-1] + 1;
269                 }
270                 assert_eq(i, numsamp);
271                 assert_lt(i, 128);
272                 if(sanityCheck) {
273                         // diffs[] records all the differences observed
274                         bool *diffs = new bool[maxv];
275                         for(T i = 0; i < numsamp; i++) {
276                                 for(T j = i+1; j < numsamp; j++) {
277                                         T d1 = (clDCs[r].samples[j] - clDCs[r].samples[i]);
278                                         T d2 = (clDCs[r].samples[i] + maxv - clDCs[r].samples[j]);
279                                         assert_lt(d1, maxv);
280                                         assert_lt(d2, maxv);
281                                         diffs[d1] = true;
282                                         diffs[d2] = true;
283                                 }
284                         }
285                         // Should have observed all possible differences (except 0)
286                         for(T i = 1; i < maxv; i++) {
287                                 if(diffs[i] == false) cout << r << ", " << i << endl;
288                                 assert(diffs[i] == true);
289                         }
290                         delete[] diffs;
291                 }
292         }
293         clDCs_calced = true;
294 }
295
296 /**
297  * Entries 4-57 are transcribed from page 6 of Luk and Wong's paper
298  * "Two New Quorum Based Algorithms for Distributed Mutual Exclusion",
299  * which is also used and cited in the Burkhardt and Karkkainen's
300  * papers on difference covers for sorting.  These samples are optimal
301  * according to Luk and Wong.
302  *
303  * All other entries are generated via the exhaustive algorithm in
304  * calcExhaustiveDC().
305  *
306  * The 0 is stored at the end of the sample as an end-of-list marker,
307  * but 0 is also an element of each.
308  *
309  * Note that every difference cover has a 0 and a 1.  Intuitively,
310  * any optimal difference cover sample can be oriented (i.e. rotated)
311  * such that it includes 0 and 1 as elements.
312  *
313  * All samples in this list have been verified to be complete covers.
314  *
315  * A value of 0xffffffff in the first column indicates that there is no
316  * sample for that value of v.  We do not keep samples for values of v
317  * less than 3, since they are trivial (and the caller probably didn't
318  * mean to ask for it).
319  */
320 static uint32_t dc0to64[65][10] = {
321         {0xffffffff},                     // 0
322         {0xffffffff},                     // 1
323         {0xffffffff},                     // 2
324         {1, 0},                           // 3
325         {1, 2, 0},                        // 4
326         {1, 2, 0},                        // 5
327         {1, 3, 0},                        // 6
328         {1, 3, 0},                        // 7
329         {1, 2, 4, 0},                     // 8
330         {1, 2, 4, 0},                     // 9
331         {1, 2, 5, 0},                     // 10
332         {1, 2, 5, 0},                     // 11
333         {1, 3, 7, 0},                     // 12
334         {1, 3, 9, 0},                     // 13
335         {1, 2, 3, 7, 0},                  // 14
336         {1, 2, 3, 7, 0},                  // 15
337         {1, 2, 5, 8, 0},                  // 16
338         {1, 2, 4, 12, 0},                 // 17
339         {1, 2, 5, 11, 0},                 // 18
340         {1, 2, 6, 9, 0},                  // 19
341         {1, 2, 3, 6, 10, 0},              // 20
342         {1, 4, 14, 16, 0},                // 21
343         {1, 2, 3, 7, 11, 0},              // 22
344         {1, 2, 3, 7, 11, 0},              // 23
345         {1, 2, 3, 7, 15, 0},              // 24
346         {1, 2, 3, 8, 12, 0},              // 25
347         {1, 2, 5, 9, 15, 0},              // 26
348         {1, 2, 5, 13, 22, 0},             // 27
349         {1, 4, 15, 20, 22, 0},            // 28
350         {1, 2, 3, 4, 9, 14, 0},           // 29
351         {1, 2, 3, 4, 9, 19, 0},           // 30
352         {1, 3, 8, 12, 18, 0},             // 31
353         {1, 2, 3, 7, 11, 19, 0},          // 32
354         {1, 2, 3, 6, 16, 27, 0},          // 33
355         {1, 2, 3, 7, 12, 20, 0},          // 34
356         {1, 2, 3, 8, 12, 21, 0},          // 35
357         {1, 2, 5, 12, 14, 20, 0},         // 36
358         {1, 2, 4, 10, 15, 22, 0},         // 37
359         {1, 2, 3, 4, 8, 14, 23, 0},       // 38
360         {1, 2, 4, 13, 18, 33, 0},         // 39
361         {1, 2, 3, 4, 9, 14, 24, 0},       // 40
362         {1, 2, 3, 4, 9, 15, 25, 0},       // 41
363         {1, 2, 3, 4, 9, 15, 25, 0},       // 42
364         {1, 2, 3, 4, 10, 15, 26, 0},      // 43
365         {1, 2, 3, 6, 16, 27, 38, 0},      // 44
366         {1, 2, 3, 5, 12, 18, 26, 0},      // 45
367         {1, 2, 3, 6, 18, 25, 38, 0},      // 46
368         {1, 2, 3, 5, 16, 22, 40, 0},      // 47
369         {1, 2, 5, 9, 20, 26, 36, 0},      // 48
370         {1, 2, 5, 24, 33, 36, 44, 0},     // 49
371         {1, 3, 8, 17, 28, 32, 38, 0},     // 50
372         {1, 2, 5, 11, 18, 30, 38, 0},     // 51
373         {1, 2, 3, 4, 6, 14, 21, 30, 0},   // 52
374         {1, 2, 3, 4, 7, 21, 29, 44, 0},   // 53
375         {1, 2, 3, 4, 9, 15, 21, 31, 0},   // 54
376         {1, 2, 3, 4, 6, 19, 26, 47, 0},   // 55
377         {1, 2, 3, 4, 11, 16, 33, 39, 0},  // 56
378         {1, 3, 13, 32, 36, 43, 52, 0},    // 57
379
380         // Generated by calcExhaustiveDC()
381         {1, 2, 3, 7, 21, 33, 37, 50, 0},  // 58
382         {1, 2, 3, 6, 13, 21, 35, 44, 0},  // 59
383         {1, 2, 4, 9, 15, 25, 30, 42, 0},  // 60
384         {1, 2, 3, 7, 15, 25, 36, 45, 0},  // 61
385         {1, 2, 4, 10, 32, 39, 46, 51, 0}, // 62
386         {1, 2, 6, 8, 20, 38, 41, 54, 0},  // 63
387         {1, 2, 5, 14, 16, 34, 42, 59, 0}  // 64
388 };
389
390 /**
391  * Get a difference cover for the requested periodicity v.
392  */
393 template <typename T>
394 static String<T> getDiffCover(T v,
395                               bool verbose = false,
396                               bool sanityCheck = false)
397 {
398         assert_gt(v, 2);
399         String<T> ret;
400
401         // Can we look it up in our hardcoded array?
402         if(v <= 64 && dc0to64[v][0] == 0xffffffff) {
403                 if(verbose) cout << "v in hardcoded area, but hardcoded entry was all-fs" << endl;
404                 return ret;
405         } else if(v <= 64) {
406                 append(ret, 0);
407                 for(size_t i = 0; i < 10; i++) {
408                         if(dc0to64[v][i] == 0) break;
409                         append(ret, dc0to64[v][i]);
410                 }
411                 if(sanityCheck) assert(dcRepOk(v, ret));
412                 return ret;
413         }
414
415         // Can we look it up in our calcColbournAndLingDCs array?
416         if(!clDCs_calced) {
417                 calcColbournAndLingDCs<uint32_t>(verbose, sanityCheck);
418                 assert(clDCs_calced);
419         }
420         for(size_t i = 0; i < 16; i++) {
421                 if(v <= clDCs[i].maxV) {
422                         for(size_t j = 0; j < clDCs[i].numSamples; j++) {
423                                 T s = clDCs[i].samples[j];
424                                 if(s >= v) {
425                                         s %= v;
426                                         for(size_t k = 0; k < length(ret); k++) {
427                                                 if(s == ret[k]) break;
428                                                 if(s < ret[k]) {
429                                                         insertValue(ret, k, s);
430                                                         break;
431                                                 }
432                                         }
433                                 } else {
434                                         append(ret, s % v);
435                                 }
436                         }
437                         if(sanityCheck) assert(dcRepOk(v, ret));
438                         return ret;
439                 }
440         }
441         cerr << "Error: Could not find a difference cover sample for v=" << v << endl;
442         throw 1;
443 }
444
445 /**
446  * Calculate and return a delta map based on the given difference cover
447  * and periodicity v.
448  */
449 template <typename T>
450 static String<T> getDeltaMap(T v, const String<T>& dc) {
451         // Declare anchor-map-related items
452         String<T> amap;
453         size_t amapEnts = 1;
454         fill(amap, v, 0xffffffff, Exact());
455         amap[0] = 0;
456         // Print out difference cover (and optionally calculate
457         // anchor map)
458         for(size_t i = 0; i < length(dc); i++) {
459                 for(size_t j = i+1; j < length(dc); j++) {
460                         assert_gt(dc[j], dc[i]);
461                         T diffLeft  = dc[j] - dc[i];
462                         T diffRight = dc[i] + v - dc[j];
463                         assert_lt(diffLeft, v);
464                         assert_lt(diffRight, v);
465                         if(amap[diffLeft] == 0xffffffff) {
466                                 amap[diffLeft] = dc[i];
467                                 amapEnts++;
468                         }
469                         if(amap[diffRight] == 0xffffffff) {
470                                 amap[diffRight] = dc[j];
471                                 amapEnts++;
472                         }
473                 }
474         }
475         return amap;
476 }
477
478 /**
479  * Return population count (count of all bits set to 1) of i.
480  */
481 template<typename T>
482 static unsigned int popCount(T i) {
483         unsigned int cnt = 0;
484         for(size_t j = 0; j < sizeof(T)*8; j++) {
485                 if(i & 1) cnt++;
486                 i >>= 1;
487         }
488         return cnt;
489 }
490
491 /**
492  * Calculate log-base-2 of i
493  */
494 template<typename T>
495 static unsigned int myLog2(T i) {
496         assert_eq(1, popCount(i)); // must be power of 2
497         for(size_t j = 0; j < sizeof(T)*8; j++) {
498                 if(i & 1) return j;
499                 i >>= 1;
500         }
501         assert(false);
502         return 0xffffffff;
503 }
504
505 /**
506  *
507  */
508 template<typename TStr>
509 class DifferenceCoverSample {
510 public:
511
512         DifferenceCoverSample(const TStr& __text,
513                               uint32_t __v,
514                               bool __verbose = false,
515                               bool __sanity = false,
516                               ostream& __logger = cout) :
517                 _text(__text),
518                 _v(__v),
519                 _verbose(__verbose),
520                 _sanity(__sanity),
521                 _ds(getDiffCover(_v, _verbose, _sanity)),
522                 _dmap(getDeltaMap(_v, _ds)),
523                 _d(length(_ds)),
524                 _doffs(),
525                 _isaPrime(),
526                 _dInv(),
527                 _log2v(myLog2(_v)),
528                 _vmask(0xffffffff << _log2v),
529                 _logger(__logger)
530         {
531                 assert_gt(_d, 0);
532                 assert_eq(1, popCount(_v)); // must be power of 2
533                 // Build map from d's to idx's
534                 fill(_dInv, _v, 0xffffffff, Exact());
535                 for(size_t i = 0; i < length(_ds); i++) _dInv[_ds[i]] = i;
536         }
537         
538         /**
539          * Allocate an amount of memory that simulates the peak memory
540          * usage of the DifferenceCoverSample with the given text and v.
541          * Throws bad_alloc if it's not going to fit in memory.  Returns
542          * the approximate number of bytes the Cover takes at all times.
543          */
544         static size_t simulateAllocs(const TStr& text, uint32_t v) {
545                 String<uint32_t> ds = getDiffCover(v, false /*verbose*/, false /*sanity*/);
546                 size_t len = length(text);
547                 size_t sPrimeSz = (len / v) * length(ds);
548                 // sPrime, sPrimeOrder, _isaPrime all exist in memory at
549                 // once and that's the peak
550                 AutoArray<uint32_t> aa(sPrimeSz * 3 + (1024 * 1024 /*out of caution*/));
551                 return sPrimeSz * 4; // sPrime array
552         }
553
554         uint32_t v() const                   { return _v; }
555         uint32_t log2v() const               { return _log2v; }
556         uint32_t vmask() const               { return _vmask; }
557         uint32_t modv(uint32_t i) const      { return i & ~_vmask; }
558         uint32_t divv(uint32_t i) const      { return i >> _log2v; }
559         uint32_t d() const                   { return _d; }
560         bool verbose() const                 { return _verbose; }
561         bool sanityCheck() const             { return _sanity; }
562         const TStr& text() const             { return _text; }
563         const String<uint32_t>& ds() const   { return _ds; }
564         const String<uint32_t>& dmap() const { return _dmap; }
565         ostream& log() const                 { return _logger; }
566
567         void     build();
568         uint32_t tieBreakOff(uint32_t i, uint32_t j) const;
569         int64_t  breakTie(uint32_t i, uint32_t j) const;
570         bool     isCovered(uint32_t i) const;
571         uint32_t rank(uint32_t i) const;
572
573         /**
574          * Print out the suffix array such that every sample offset has its
575          * rank filled in and every non-sample offset is shown as '-'.
576          */
577         void print(ostream& out) {
578                 for(size_t i = 0; i < length(_text); i++) {
579                         if(isCovered(i)) {
580                                 out << rank(i);
581                         } else {
582                                 out << "-";
583                         }
584                         if(i < length(_text)-1) {
585                                 out << ",";
586                         }
587                 }
588                 out << endl;
589         }
590
591 private:
592
593         void doBuiltSanityCheck() const;
594         void buildSPrime(String<uint32_t>& sPrime);
595
596         bool built() const {
597                 return length(_isaPrime) > 0;
598         }
599
600         void verbose(const string& s) const {
601                 if(this->verbose()) {
602                         this->log() << s;
603                         this->log().flush();
604                 }
605         }
606
607         const TStr&      _text;     // text to sample
608         uint32_t         _v;        // periodicity of sample
609         bool             _verbose;  //
610         bool             _sanity;   //
611         String<uint32_t> _ds;       // samples: idx -> d
612         String<uint32_t> _dmap;     // delta map
613         uint32_t         _d;        // |D| - size of sample
614         String<uint32_t> _doffs;    // offsets into sPrime/isaPrime for each d idx
615         String<uint32_t> _isaPrime; // ISA' array
616         String<uint32_t> _dInv;     // Map from d -> idx
617         uint32_t         _log2v;
618         uint32_t         _vmask;
619         ostream&         _logger;
620 };
621
622 /**
623  * Return true iff suffixes with offsets suf1 and suf2 out of host
624  * string 'host' are identical up to depth 'v'.
625  */
626 template <typename TStr>
627 static inline bool suffixLt(const TStr& host, uint32_t suf1, uint32_t suf2) {
628         uint32_t hlen = length(host);
629         assert_neq(suf1, suf2);
630         uint32_t i = 0;
631         while(suf1 + i < hlen && suf2 + i < hlen) {
632                 if(host[suf1+i] < host[suf2+i]) return true;
633                 if(host[suf1+i] > host[suf2+i]) return false;
634                 i++;
635         }
636         if(suf1 + i == hlen) {
637                 assert_lt(suf2 + i, hlen);
638                 return false;
639         }
640         assert_eq(suf2 + i, hlen);
641         return true;
642 }
643
644 /**
645  * Sanity-check the difference cover by first inverting _isaPrime then
646  * checking that each successive suffix really is less than the next.
647  */
648 template <typename TStr>
649 void DifferenceCoverSample<TStr>::doBuiltSanityCheck() const {
650         uint32_t v = this->v();
651         assert(built());
652         VMSG_NL("  Doing sanity check");
653         uint32_t added = 0;
654         String<uint32_t> sorted;
655         fill(sorted, length(_isaPrime), 0xffffffff, Exact());
656         for(size_t di = 0; di < this->d(); di++) {
657                 uint32_t d = _ds[di];
658                 size_t i = 0;
659                 for(size_t doi = _doffs[di]; doi < _doffs[di+1]; doi++, i++) {
660                         assert_eq(0xffffffff, sorted[_isaPrime[doi]]);
661                         // Maps the offset of the suffix to its rank
662                         sorted[_isaPrime[doi]] = v*i + d;
663                         added++;
664                 }
665         }
666         assert_eq(added, length(_isaPrime));
667         for(size_t i = 0; i < length(sorted)-1; i++) {
668                 assert(suffixLt(this->text(), sorted[i], sorted[i+1]));
669         }
670 }
671
672 /**
673  * Build the s' array by sampling suffixes (suffix offsets, actually)
674  * from t according to the difference-cover sample and pack them into
675  * an array of machine words in the order dictated by the "mu" mapping
676  * described in Burkhardt.
677  *
678  * Also builds _doffs map.
679  */
680 template <typename TStr>
681 void DifferenceCoverSample<TStr>::buildSPrime(String<uint32_t>& sPrime) {
682         const TStr& t = this->text();
683         const String<uint32_t>& ds = this->ds();
684         uint32_t tlen = length(t);
685         uint32_t v = this->v();
686         uint32_t d = this->d();
687         assert_gt(v, 2);
688         assert_lt(d, v);
689         // Record where each d section should begin in sPrime
690         uint32_t tlenDivV = this->divv(tlen);
691         uint32_t tlenModV = this->modv(tlen);
692         uint32_t sPrimeSz = 0;
693         assert(empty(_doffs));
694         reserve(_doffs, d+1, Exact());
695         assert_eq(capacity(_doffs), d+1);
696         for(uint32_t di = 0; di < d; di++) {
697                 // mu mapping
698                 uint32_t sz = tlenDivV + ((ds[di] <= tlenModV) ? 1 : 0);
699                 assert_geq(sz, 0);
700                 appendValue(_doffs, sPrimeSz);
701                 sPrimeSz += sz;
702         }
703         appendValue(_doffs, sPrimeSz);
704         #ifndef NDEBUG
705         if(tlenDivV > 0) {
706                 for(size_t i = 0; i < d; i++) {
707                         assert_gt(_doffs[i+1], _doffs[i]);
708                         uint32_t diff = _doffs[i+1] - _doffs[i];
709                         assert(diff == tlenDivV || diff == tlenDivV+1);
710                 }
711         }
712         #endif
713         assert_eq(length(_doffs), d+1);
714         // Size sPrime appropriately
715         reserve(sPrime, sPrimeSz+1, Exact()); // reserve extra slot for LS
716         fill(sPrime, sPrimeSz, 0xffffffff, Exact());
717         // Slot suffixes from text into sPrime according to the mu
718         // mapping; where the mapping would leave a blank, insert a 0
719         uint32_t added = 0;
720         uint32_t i = 0;
721         for(uint32_t ti = 0; ti <= tlen; ti += v) {
722                 for(uint32_t di = 0; di < d; di++) {
723                         uint32_t tti = ti + ds[di];
724                         if(tti > tlen) break;
725                         uint32_t spi = _doffs[di] + i;
726                         assert_lt(spi, _doffs[di+1]);
727                         assert_leq(tti, tlen);
728                         assert_lt(spi, sPrimeSz);
729                         assert_eq(0xffffffff, sPrime[spi]);
730                         sPrime[spi] = tti; added++;
731                 }
732                 i++;
733         }
734         assert_eq(added, sPrimeSz);
735 }
736
737 /**
738  * Return true iff suffixes with offsets suf1 and suf2 out of host
739  * string 'host' are identical up to depth 'v'.
740  */
741 template <typename TStr>
742 static inline bool suffixSameUpTo(const TStr& host,
743                                   uint32_t suf1,
744                                   uint32_t suf2,
745                                   uint32_t v)
746 {
747         for(uint32_t i = 0; i < v; i++) {
748                 bool endSuf1 = suf1+i >= length(host);
749                 bool endSuf2 = suf2+i >= length(host);
750                 if((endSuf1 && !endSuf2) || (!endSuf1 && endSuf2)) return false;
751                 if(endSuf1 && endSuf2) return true;
752                 if(host[suf1+i] != host[suf2+i]) return false;
753         }
754         return true;
755 }
756
757 /**
758  * Calculates a ranking of all suffixes in the sample and stores them,
759  * packed according to the mu mapping, in _isaPrime.
760  */
761 template <typename TStr>
762 void DifferenceCoverSample<TStr>::build() {
763         // Local names for relevant types
764         typedef typename Value<TStr>::Type TAlphabet;
765         VMSG_NL("Building DifferenceCoverSample");
766         // Local names for relevant data
767         const TStr& t = this->text();
768         uint32_t v = this->v();
769         assert_gt(v, 2);
770         // Build s'
771         String<uint32_t> sPrime;
772         VMSG_NL("  Building sPrime");
773         buildSPrime(sPrime);
774         assert_gt(length(sPrime), 0);
775         assert_leq(length(sPrime), length(t)+1); // +1 is because of the end-cap
776         uint32_t nextRank = 0;
777         {
778                 VMSG_NL("  Building sPrimeOrder");
779                 String<uint32_t> sPrimeOrder;
780                 reserve(sPrimeOrder, length(sPrime)+1, Exact()); // reserve extra slot for LS
781                 resize(sPrimeOrder, length(sPrime), Exact());
782                 for(size_t i = 0; i < length(sPrimeOrder); i++) {
783                         sPrimeOrder[i] = i;
784                 }
785                 // sPrime now holds suffix-offsets for DC samples.
786                 {
787                         Timer timer(cout, "  V-Sorting samples time: ", this->verbose());
788                         VMSG_NL("  V-Sorting samples");
789                         // Extract backing-store array from sPrime and sPrimeOrder;
790                         // the mkeyQSortSuf2 routine works on the array for maximum
791                         // efficiency
792                         uint32_t *sPrimeArr = (uint32_t*)begin(sPrime);
793                         size_t slen = length(sPrime);
794                         assert_eq(sPrimeArr[0], sPrime[0]);
795                         assert_eq(sPrimeArr[slen-1], sPrime[slen-1]);
796                         uint32_t *sPrimeOrderArr = (uint32_t*)begin(sPrimeOrder);
797                         assert_eq(sPrimeOrderArr[0], sPrimeOrder[0]);
798                         assert_eq(sPrimeOrderArr[slen-1], sPrimeOrder[slen-1]);
799                         // Sort sample suffixes up to the vth character using a
800                         // multikey quicksort.  Sort time is proportional to the
801                         // number of samples times v.  It isn't quadratic.
802                         // sPrimeOrder is passed in as a swapping partner for
803                         // sPrimeArr, i.e., every time the multikey qsort swaps
804                         // elements in sPrime, it swaps the same elements in
805                         // sPrimeOrder too.  This allows us to easily reconstruct
806                         // what the sort did.
807                         mkeyQSortSuf2(t, sPrimeArr, slen, sPrimeOrderArr,
808                                       ValueSize<TAlphabet>::VALUE,
809                                       this->verbose(), this->sanityCheck(), v);
810                         // Make sure sPrime and sPrimeOrder are consistent with
811                         // their respective backing-store arrays
812                         assert_eq(sPrimeArr[0], sPrime[0]);
813                         assert_eq(sPrimeArr[slen-1], sPrime[slen-1]);
814                         assert_eq(sPrimeOrderArr[0], sPrimeOrder[0]);
815                         assert_eq(sPrimeOrderArr[slen-1], sPrimeOrder[slen-1]);
816                 }
817                 // Now assign the ranking implied by the sorted sPrime/sPrimeOrder
818                 // arrays back into sPrime.
819                 VMSG_NL("  Allocating rank array");
820                 reserve(_isaPrime, length(sPrime)+1, Exact());
821                 fill(_isaPrime, length(sPrime), 0xffffffff, Exact());
822                 assert_gt(length(_isaPrime), 0);
823                 {
824                         Timer timer(cout, "  Ranking v-sort output time: ", this->verbose());
825                         VMSG_NL("  Ranking v-sort output");
826                         for(size_t i = 0; i < length(sPrime)-1; i++) {
827                                 // Place the appropriate ranking
828                                 _isaPrime[sPrimeOrder[i]] = nextRank;
829                                 // If sPrime[i] and sPrime[i+1] are identical up to v, then we
830                                 // should give the next suffix the same rank
831                                 if(!suffixSameUpTo(t, sPrime[i], sPrime[i+1], v)) nextRank++;
832                         }
833                         _isaPrime[sPrimeOrder[length(sPrime)-1]] = nextRank; // finish off
834                 }
835                 // sPrimeOrder is destroyed
836                 // All the information we need is now in _isaPrime
837         }
838         #ifndef NDEBUG
839         // Check that all ranks are sane
840         for(size_t i = 0; i < length(_isaPrime); i++) {
841                 assert_neq(_isaPrime[i], 0xffffffff);
842                 assert_lt(_isaPrime[i], length(_isaPrime));
843         }
844         #endif
845         // Now pass the sPrimeRanks[] array to LarssonSadakane (in SeqAn).
846         append(_isaPrime, length(_isaPrime));
847         append(sPrime, length(sPrime));
848         {
849                 Timer timer(cout, "  Invoking Larsson-Sadakane on ranks time: ", this->verbose());
850                 VMSG_NL("  Invoking Larsson-Sadakane on ranks");
851                 createSuffixArray(sPrime, _isaPrime, LarssonSadakane(), length(_isaPrime));
852         }
853         // sPrime now contains the suffix array (which we ignore)
854         assert_eq(length(_isaPrime), length(sPrime));
855         assert_gt(length(_isaPrime), 0);
856         // chop off final character of _isaPrime
857         resize(_isaPrime, length(_isaPrime)-1);
858         // Subtract 1 from each isaPrime (to adjust for LarssonSadakane
859         // always ranking the final suffix as lexicographically first)
860         for(size_t i = 0; i < length(_isaPrime); i++) {
861                 _isaPrime[i]--;
862         }
863         VMSG_NL("  Sanity-checking and returning");
864
865         // done!
866         //if(this->verbose()) print(cout);
867         if(this->sanityCheck()) doBuiltSanityCheck();
868 }
869
870 /**
871  * Return true iff index i within the text is covered by the difference
872  * cover sample.  Allow i to be off the end of the text; simplifies
873  * logic elsewhere.
874  */
875 template <typename TStr>
876 bool DifferenceCoverSample<TStr>::isCovered(uint32_t i) const {
877         assert(built());
878         uint32_t modi = this->modv(i);
879         assert_lt(modi, length(_dInv));
880         return _dInv[modi] != 0xffffffff;
881 }
882
883 /**
884  * Given a text offset that's covered, return its lexicographical rank
885  * among the sample suffixes.
886  */
887 template <typename TStr>
888 uint32_t DifferenceCoverSample<TStr>::rank(uint32_t i) const {
889         assert(built());
890         assert_lt(i, length(this->text()));
891         uint32_t imodv = this->modv(i);
892         assert_neq(0xffffffff, _dInv[imodv]); // must be in the sample
893         uint32_t ioff = this->divv(i);
894         assert_lt(ioff, _doffs[_dInv[imodv]+1] - _doffs[_dInv[imodv]]);
895         uint32_t isaIIdx = _doffs[_dInv[imodv]] + ioff;
896         assert_lt(isaIIdx, length(_isaPrime));
897         uint32_t isaPrimeI = _isaPrime[isaIIdx];
898         assert_leq(isaPrimeI, length(_isaPrime));
899         return isaPrimeI;
900 }
901
902 /**
903  * Return: < 0 if suffix i is lexicographically less than suffix j; > 0
904  * if suffix j is lexicographically greater.
905  */
906 template <typename TStr>
907 int64_t DifferenceCoverSample<TStr>::breakTie(uint32_t i, uint32_t j) const {
908         assert(built());
909         assert_neq(i, j);
910         assert_lt(i, length(this->text()));
911         assert_lt(j, length(this->text()));
912         uint32_t imodv = this->modv(i);
913         uint32_t jmodv = this->modv(j);
914         assert_neq(0xffffffff, _dInv[imodv]); // must be in the sample
915         assert_neq(0xffffffff, _dInv[jmodv]); // must be in the sample
916         uint32_t dimodv = _dInv[imodv];
917         uint32_t djmodv = _dInv[jmodv];
918         uint32_t ioff = this->divv(i);
919         uint32_t joff = this->divv(j);
920         assert_lt(dimodv+1, length(_doffs));
921         assert_lt(djmodv+1, length(_doffs));
922         // assert_lt: expected (32024) < (0)
923         assert_lt(ioff, _doffs[dimodv+1] - _doffs[dimodv]);
924         assert_lt(joff, _doffs[djmodv+1] - _doffs[djmodv]);
925         uint32_t isaIIdx = _doffs[dimodv] + ioff;
926         uint32_t isaJIdx = _doffs[djmodv] + joff;
927         assert_lt(isaIIdx, length(_isaPrime));
928         assert_lt(isaJIdx, length(_isaPrime));
929         assert_neq(isaIIdx, isaJIdx); // ranks must be unique
930         uint32_t isaPrimeI = _isaPrime[isaIIdx];
931         uint32_t isaPrimeJ = _isaPrime[isaJIdx];
932         assert_neq(isaPrimeI, isaPrimeJ); // ranks must be unique
933         assert_leq(isaPrimeI, length(_isaPrime));
934         assert_leq(isaPrimeJ, length(_isaPrime));
935         return (int64_t)isaPrimeI - (int64_t)isaPrimeJ;
936 }
937
938 /**
939  * Given i, j, return the number of additional characters that need to
940  * be compared before the difference cover can break the tie.
941  */
942 template <typename TStr>
943 uint32_t DifferenceCoverSample<TStr>::tieBreakOff(uint32_t i, uint32_t j) const {
944         const TStr& t = this->text();
945         const String<uint32_t>& dmap = this->dmap();
946         assert(built());
947         // It's actually convenient to allow this, but we're permitted to
948         // return nonsense in that case
949         if(t[i] != t[j]) return 0xffffffff;
950         //assert_eq(t[i], t[j]); // if they're unequal, there's no tie to break
951         uint32_t v = this->v();
952         assert_neq(i, j);
953         assert_lt(i, length(t));
954         assert_lt(j, length(t));
955         uint32_t imod = this->modv(i);
956         uint32_t jmod = this->modv(j);
957         uint32_t diffLeft = (jmod >= imod)? (jmod - imod) : (jmod + v - imod);
958         uint32_t diffRight = (imod >= jmod)? (imod - jmod) : (imod + v - jmod);
959         assert_lt(diffLeft, length(dmap));
960         assert_lt(diffRight, length(dmap));
961         uint32_t destLeft = dmap[diffLeft];   // offset where i needs to be
962         uint32_t destRight = dmap[diffRight]; // offset where i needs to be
963         assert(isCovered(destLeft));
964         assert(isCovered(destLeft+diffLeft));
965         assert(isCovered(destRight));
966         assert(isCovered(destRight+diffRight));
967         assert_lt(destLeft, v);
968         assert_lt(destRight, v);
969         uint32_t deltaLeft = (destLeft >= imod)? (destLeft - imod) : (destLeft + v - imod);
970         if(deltaLeft == v) deltaLeft = 0;
971         uint32_t deltaRight = (destRight >= jmod)? (destRight - jmod) : (destRight + v - jmod);
972         if(deltaRight == v) deltaRight = 0;
973         assert_lt(deltaLeft, v);
974         assert_lt(deltaRight, v);
975         assert(isCovered(i+deltaLeft));
976         assert(isCovered(j+deltaLeft));
977         assert(isCovered(i+deltaRight));
978         assert(isCovered(j+deltaRight));
979         return min(deltaLeft, deltaRight);
980 }
981
982 #endif /*DIFF_SAMPLE_H_*/