ebwt.h

   1 #ifndef EBWT_H_
   2 #define EBWT_H_
   3
   4 #include <stdint.h>
   5 #include <string.h>
   6 #include <algorithm>
   7 #include <iostream>
   8 #include <fstream>
   9 #include <vector>
  10 #include <set>
  11 #include <sstream>
  12 #include <fcntl.h>
  13 #include <errno.h>
  14 #include <stdexcept>
  15 #include <seqan/sequence.h>
  16 #include <seqan/index.h>
  17 #include <sys/stat.h>
  18 #ifdef BOWTIE_MM
  19 #include <sys/mman.h>
  20 #include <sys/shm.h>
  21 #endif
  22 #include "auto_array.h"
  23 #include "shmem.h"
  24 #include "alphabet.h"
  25 #include "assert_helpers.h"
  26 #include "bitpack.h"
  27 #include "blockwise_sa.h"
  28 #include "endian_swap.h"
  29 #include "word_io.h"
  30 #include "random_source.h"
  31 #include "hit.h"
  32 #include "ref_read.h"
  33 #include "threading.h"
  34 #include "bitset.h"
  35 #include "str_util.h"
  36 #include "mm.h"
  37 #include "timer.h"
  38 #include "refmap.h"
  39 #include "color_dec.h"
  40 #include "reference.h"
  41
  42 using namespace std;
  43 using namespace seqan;
  44
  45 #ifndef PREFETCH_LOCALITY
  46 // No locality by default
  47 #define PREFETCH_LOCALITY 2
  48 #endif
  49
  50 // From ccnt_lut.cpp, automatically generated by gen_lookup_tables.pl
  51 extern uint8_t cCntLUT_4[4][4][256];
  52
  53 #ifndef VMSG_NL
  54 #define VMSG_NL(args...) \
  55 if(this->verbose()) { \
  56         stringstream tmp; \
  57         tmp << args << endl; \
  58         this->verbose(tmp.str()); \
  59 }
  60 #endif
  61
  62 #ifndef VMSG
  63 #define VMSG(args...) \
  64 if(this->verbose()) { \
  65         stringstream tmp; \
  66         tmp << args; \
  67         this->verbose(tmp.str()); \
  68 }
  69 #endif
  70
  71 /**
  72  * Flags describing type of Ebwt.
  73  */
  74 enum EBWT_FLAGS {
  75         EBWT_COLOR = 2,     // true -> Ebwt is colorspace
  76         EBWT_ENTIRE_REV = 4 // true -> reverse Ebwt is the whole
  77                             // concatenated string reversed, rather than
  78                             // each stretch reversed
  79 };
  80
  81 /**
  82  * Extended Burrows-Wheeler transform header.  This together with the
  83  * actual data arrays and other text-specific parameters defined in
  84  * class Ebwt constitute the entire Ebwt.
  85  */
  86 class EbwtParams {
  87
  88 public:
  89         EbwtParams() { }
  90
  91         EbwtParams(uint32_t len,
  92                    int32_t lineRate,
  93                    int32_t linesPerSide,
  94                    int32_t offRate,
  95                    int32_t isaRate,
  96                    int32_t ftabChars,
  97                    bool color,
  98                    bool entireReverse)
  99         {
 100                 init(len, lineRate, linesPerSide, offRate, isaRate, ftabChars, color, entireReverse);
 101         }
 102
 103         EbwtParams(const EbwtParams& eh) {
 104                 init(eh._len, eh._lineRate, eh._linesPerSide, eh._offRate,
 105                      eh._isaRate, eh._ftabChars, eh._color, eh._entireReverse);
 106         }
 107
 108         void init(uint32_t len, int32_t lineRate, int32_t linesPerSide,
 109                   int32_t offRate, int32_t isaRate, int32_t ftabChars,
 110                   bool color, bool entireReverse)
 111         {
 112                 _color = color;
 113                 _entireReverse = entireReverse;
 114                 _len = len;
 115                 _bwtLen = _len + 1;
 116                 _sz = (len+3)/4;
 117                 _bwtSz = (len/4 + 1);
 118                 _lineRate = lineRate;
 119                 _linesPerSide = linesPerSide;
 120                 _origOffRate = offRate;
 121                 _offRate = offRate;
 122                 _offMask = 0xffffffff << _offRate;
 123                 _isaRate = isaRate;
 124                 _isaMask = 0xffffffff << ((_isaRate >= 0) ? _isaRate : 0);
 125                 _ftabChars = ftabChars;
 126                 _eftabLen = _ftabChars*2;
 127                 _eftabSz = _eftabLen*4;
 128                 _ftabLen = (1 << (_ftabChars*2))+1;
 129                 _ftabSz = _ftabLen*4;
 130                 _offsLen = (_bwtLen + (1 << _offRate) - 1) >> _offRate;
 131                 _offsSz = _offsLen*4;
 132                 _isaLen = (_isaRate == -1)? 0 : ((_bwtLen + (1 << _isaRate) - 1) >> _isaRate);
 133                 _isaSz = _isaLen*4;
 134                 _lineSz = 1 << _lineRate;
 135                 _sideSz = _lineSz * _linesPerSide;
 136                 _sideBwtSz = _sideSz - 8;
 137                 _sideBwtLen = _sideBwtSz*4;
 138                 _numSidePairs = (_bwtSz+(2*_sideBwtSz)-1)/(2*_sideBwtSz);
 139                 _numSides = _numSidePairs*2;
 140                 _numLines = _numSides * _linesPerSide;
 141                 _ebwtTotLen = _numSidePairs * (2*_sideSz);
 142                 _ebwtTotSz = _ebwtTotLen;
 143                 assert(repOk());
 144         }
 145
 146         uint32_t len() const           { return _len; }
 147         uint32_t bwtLen() const        { return _bwtLen; }
 148         uint32_t sz() const            { return _sz; }
 149         uint32_t bwtSz() const         { return _bwtSz; }
 150         int32_t  lineRate() const      { return _lineRate; }
 151         int32_t  linesPerSide() const  { return _linesPerSide; }
 152         int32_t  origOffRate() const   { return _origOffRate; }
 153         int32_t  offRate() const       { return _offRate; }
 154         uint32_t offMask() const       { return _offMask; }
 155         int32_t  isaRate() const       { return _isaRate; }
 156         uint32_t isaMask() const       { return _isaMask; }
 157         int32_t  ftabChars() const     { return _ftabChars; }
 158         uint32_t eftabLen() const      { return _eftabLen; }
 159         uint32_t eftabSz() const       { return _eftabSz; }
 160         uint32_t ftabLen() const       { return _ftabLen; }
 161         uint32_t ftabSz() const        { return _ftabSz; }
 162         uint32_t offsLen() const       { return _offsLen; }
 163         uint32_t offsSz() const        { return _offsSz; }
 164         uint32_t isaLen() const        { return _isaLen; }
 165         uint32_t isaSz() const         { return _isaSz; }
 166         uint32_t lineSz() const        { return _lineSz; }
 167         uint32_t sideSz() const        { return _sideSz; }
 168         uint32_t sideBwtSz() const     { return _sideBwtSz; }
 169         uint32_t sideBwtLen() const    { return _sideBwtLen; }
 170         uint32_t numSidePairs() const  { return _numSidePairs; }
 171         uint32_t numSides() const      { return _numSides; }
 172         uint32_t numLines() const      { return _numLines; }
 173         uint32_t ebwtTotLen() const    { return _ebwtTotLen; }
 174         uint32_t ebwtTotSz() const     { return _ebwtTotSz; }
 175         bool color() const             { return _color; }
 176         bool entireReverse() const     { return _entireReverse; }
 177
 178         /**
 179          * Set a new suffix-array sampling rate, which involves updating
 180          * rate, mask, sample length, and sample size.
 181          */
 182         void setOffRate(int __offRate) {
 183                 _offRate = __offRate;
 184                 _offMask = 0xffffffff << _offRate;
 185                 _offsLen = (_bwtLen + (1 << _offRate) - 1) >> _offRate;
 186                 _offsSz = _offsLen*4;
 187         }
 188
 189         /**
 190          * Set a new inverse suffix-array sampling rate, which involves
 191          * updating rate, mask, sample length, and sample size.
 192          */
 193         void setIsaRate(int __isaRate) {
 194                 _isaRate = __isaRate;
 195                 _isaMask = 0xffffffff << _isaRate;
 196                 _isaLen = (_bwtLen + (1 << _isaRate) - 1) >> _isaRate;
 197                 _isaSz = _isaLen*4;
 198         }
 199
 200         /// Check that this EbwtParams is internally consistent
 201         bool repOk() const {
 202                 assert_gt(_len, 0);
 203                 assert_gt(_lineRate, 3);
 204                 assert_geq(_offRate, 0);
 205                 assert_leq(_ftabChars, 16);
 206                 assert_geq(_ftabChars, 1);
 207                 assert_lt(_lineRate, 32);
 208                 assert_lt(_linesPerSide, 32);
 209                 assert_lt(_ftabChars, 32);
 210                 assert_eq(0, _ebwtTotSz % (2*_lineSz));
 211                 return true;
 212         }
 213
 214         /**
 215          * Pretty-print the header contents to the given output stream.
 216          */
 217         void print(ostream& out) const {
 218                 out << "Headers:" << endl
 219                     << "    len: "          << _len << endl
 220                     << "    bwtLen: "       << _bwtLen << endl
 221                     << "    sz: "           << _sz << endl
 222                     << "    bwtSz: "        << _bwtSz << endl
 223                     << "    lineRate: "     << _lineRate << endl
 224                     << "    linesPerSide: " << _linesPerSide << endl
 225                     << "    offRate: "      << _offRate << endl
 226                     << "    offMask: 0x"    << hex << _offMask << dec << endl
 227                     << "    isaRate: "      << _isaRate << endl
 228                     << "    isaMask: 0x"    << hex << _isaMask << dec << endl
 229                     << "    ftabChars: "    << _ftabChars << endl
 230                     << "    eftabLen: "     << _eftabLen << endl
 231                     << "    eftabSz: "      << _eftabSz << endl
 232                     << "    ftabLen: "      << _ftabLen << endl
 233                     << "    ftabSz: "       << _ftabSz << endl
 234                     << "    offsLen: "      << _offsLen << endl
 235                     << "    offsSz: "       << _offsSz << endl
 236                     << "    isaLen: "       << _isaLen << endl
 237                     << "    isaSz: "        << _isaSz << endl
 238                     << "    lineSz: "       << _lineSz << endl
 239                     << "    sideSz: "       << _sideSz << endl
 240                     << "    sideBwtSz: "    << _sideBwtSz << endl
 241                     << "    sideBwtLen: "   << _sideBwtLen << endl
 242                     << "    numSidePairs: " << _numSidePairs << endl
 243                     << "    numSides: "     << _numSides << endl
 244                     << "    numLines: "     << _numLines << endl
 245                     << "    ebwtTotLen: "   << _ebwtTotLen << endl
 246                     << "    ebwtTotSz: "    << _ebwtTotSz << endl
 247                     << "    reverse: "      << _entireReverse << endl;
 248         }
 249
 250         uint32_t _len;
 251         uint32_t _bwtLen;
 252         uint32_t _sz;
 253         uint32_t _bwtSz;
 254         int32_t  _lineRate;
 255         int32_t  _linesPerSide;
 256         int32_t  _origOffRate;
 257         int32_t  _offRate;
 258         uint32_t _offMask;
 259         int32_t  _isaRate;
 260         uint32_t _isaMask;
 261         int32_t  _ftabChars;
 262         uint32_t _eftabLen;
 263         uint32_t _eftabSz;
 264         uint32_t _ftabLen;
 265         uint32_t _ftabSz;
 266         uint32_t _offsLen;
 267         uint32_t _offsSz;
 268         uint32_t _isaLen;
 269         uint32_t _isaSz;
 270         uint32_t _lineSz;
 271         uint32_t _sideSz;
 272         uint32_t _sideBwtSz;
 273         uint32_t _sideBwtLen;
 274         uint32_t _numSidePairs;
 275         uint32_t _numSides;
 276         uint32_t _numLines;
 277         uint32_t _ebwtTotLen;
 278         uint32_t _ebwtTotSz;
 279         bool     _color;
 280         bool     _entireReverse;
 281 };
 282
 283 /**
 284  * Exception to throw when a file-realted error occurs.
 285  */
 286 class EbwtFileOpenException : public std::runtime_error {
 287 public:
 288         EbwtFileOpenException(const std::string& msg = "") :
 289                 std::runtime_error(msg) { }
 290 };
 291
 292 /**
 293  * Calculate size of file with given name.
 294  */
 295 static inline int64_t fileSize(const char* name) {
 296         std::ifstream f;
 297         f.open(name, std::ios_base::binary | std::ios_base::in);
 298         if (!f.good() || f.eof() || !f.is_open()) { return 0; }
 299         f.seekg(0, std::ios_base::beg);
 300         std::ifstream::pos_type begin_pos = f.tellg();
 301         f.seekg(0, std::ios_base::end);
 302         return static_cast<int64_t>(f.tellg() - begin_pos);
 303 }
 304
 305 // Forward declarations for Ebwt class
 306 class SideLocus;
 307 template<typename TStr> class EbwtSearchParams;
 308
 309 /**
 310  * Extended Burrows-Wheeler transform data.
 311  *
 312  * An Ebwt may be transferred to and from RAM with calls to
 313  * evictFromMemory() and loadIntoMemory().  By default, a newly-created
 314  * Ebwt is not loaded into memory; if the user would like to use a
 315  * newly-created Ebwt to answer queries, they must first call
 316  * loadIntoMemory().
 317  */
 318 template <typename TStr>
 319 class Ebwt {
 320 public:
 321         typedef typename Value<TStr>::Type TAlphabet;
 322
 323         #define Ebwt_INITS \
 324             _toBigEndian(currentlyBigEndian()), \
 325             _overrideOffRate(__overrideOffRate), \
 326             _overrideIsaRate(__overrideIsaRate), \
 327             _verbose(verbose), \
 328             _passMemExc(passMemExc), \
 329             _sanity(sanityCheck), \
 330             _fw(__fw), \
 331             _in1(MM_FILE_INIT), \
 332             _in2(MM_FILE_INIT), \
 333             _zOff(0xffffffff), \
 334             _zEbwtByteOff(0xffffffff), \
 335             _zEbwtBpOff(-1), \
 336             _nPat(0), \
 337             _nFrag(0), \
 338             _plen(NULL), \
 339             _rstarts(NULL), \
 340             _fchr(NULL), \
 341             _ftab(NULL), \
 342             _eftab(NULL), \
 343             _offs(NULL), \
 344             _isa(NULL), \
 345             _ebwt(NULL), \
 346             _useMm(false), \
 347             useShmem_(false), \
 348             _refnames(), \
 349             rmap_(NULL), \
 350             mmFile1_(NULL), \
 351             mmFile2_(NULL)
 352
 353 #ifdef EBWT_STATS
 354 #define Ebwt_STAT_INITS \
 355         ,mapLFExs_(0llu), \
 356         mapLFs_(0llu), \
 357         mapLFcs_(0llu), \
 358         mapLF1cs_(0llu), \
 359         mapLF1s_(0llu)
 360 #else
 361 #define Ebwt_STAT_INITS
 362 #endif
 363
 364         /// Construct an Ebwt from the given input file
 365         Ebwt(const string& in,
 366              int color,
 367              int needEntireReverse,
 368              bool __fw,
 369              int32_t __overrideOffRate = -1,
 370              int32_t __overrideIsaRate = -1,
 371              bool useMm = false,
 372              bool useShmem = false,
 373              bool mmSweep = false,
 374              bool loadNames = false,
 375              const ReferenceMap* rmap = NULL,
 376              bool verbose = false,
 377              bool startVerbose = false,
 378              bool passMemExc = false,
 379              bool sanityCheck = false) :
 380              Ebwt_INITS
 381              Ebwt_STAT_INITS
 382         {
 383                 assert(!useMm || !useShmem);
 384                 rmap_ = rmap;
 385                 _useMm = useMm;
 386                 useShmem_ = useShmem;
 387                 _in1Str = in + ".1.ebwt";
 388                 _in2Str = in + ".2.ebwt";
 389                 readIntoMemory(
 390                         color,         // expect colorspace reference?
 391                         __fw ? -1 : needEntireReverse, // need REF_READ_REVERSE
 392                         true,          // stop after loading the header portion?
 393                         &_eh,          // params structure to fill in
 394                         mmSweep,       // mmSweep
 395                         loadNames,     // loadNames
 396                         startVerbose); // startVerbose
 397                 // If the offRate has been overridden, reflect that in the
 398                 // _eh._offRate field
 399                 if(_overrideOffRate > _eh._offRate) {
 400                         _eh.setOffRate(_overrideOffRate);
 401                         assert_eq(_overrideOffRate, _eh._offRate);
 402                 }
 403                 // Same with isaRate
 404                 if(_overrideIsaRate > _eh._isaRate) {
 405                         _eh.setIsaRate(_overrideIsaRate);
 406                         assert_eq(_overrideIsaRate, _eh._isaRate);
 407                 }
 408                 assert(repOk());
 409         }
 410
 411         /// Construct an Ebwt from the given header parameters and string
 412         /// vector, optionally using a blockwise suffix sorter with the
 413         /// given 'bmax' and 'dcv' parameters.  The string vector is
 414         /// ultimately joined and the joined string is passed to buildToDisk().
 415         Ebwt(int color,
 416              int32_t lineRate,
 417              int32_t linesPerSide,
 418              int32_t offRate,
 419              int32_t isaRate,
 420              int32_t ftabChars,
 421              const string& file,   // base filename for EBWT files
 422              bool __fw,
 423              bool useBlockwise,
 424              uint32_t bmax,
 425              uint32_t bmaxSqrtMult,
 426              uint32_t bmaxDivN,
 427              int dcv,
 428              vector<FileBuf*>& is,
 429              vector<RefRecord>& szs,
 430              vector<uint32_t>& plens,
 431              uint32_t sztot,
 432              const RefReadInParams& refparams,
 433              uint32_t seed,
 434              int32_t __overrideOffRate = -1,
 435              int32_t __overrideIsaRate = -1,
 436              bool verbose = false,
 437              bool passMemExc = false,
 438              bool sanityCheck = false) :
 439              Ebwt_INITS
 440              Ebwt_STAT_INITS,
 441              _eh(joinedLen(szs),
 442                  lineRate,
 443                  linesPerSide,
 444                  offRate,
 445                  isaRate,
 446                  ftabChars,
 447                  color,
 448                  refparams.reverse == REF_READ_REVERSE)
 449         {
 450                 _in1Str = file + ".1.ebwt";
 451                 _in2Str = file + ".2.ebwt";
 452                 // Open output files
 453                 ofstream fout1(_in1Str.c_str(), ios::binary);
 454                 if(!fout1.good()) {
 455                         cerr << "Could not open index file for writing: \"" << _in1Str << "\"" << endl
 456                              << "Please make sure the directory exists and that permissions allow writing by" << endl
 457                              << "Bowtie." << endl;
 458                         throw 1;
 459                 }
 460                 ofstream fout2(_in2Str.c_str(), ios::binary);
 461                 if(!fout2.good()) {
 462                         cerr << "Could not open index file for writing: \"" << _in2Str << "\"" << endl
 463                              << "Please make sure the directory exists and that permissions allow writing by" << endl
 464                              << "Bowtie." << endl;
 465                         throw 1;
 466                 }
 467                 // Build
 468                 initFromVector(
 469                         is,
 470                         szs,
 471                         plens,
 472                         sztot,
 473                         refparams,
 474                         fout1,
 475                         fout2,
 476                         useBlockwise,
 477                         bmax,
 478                         bmaxSqrtMult,
 479                         bmaxDivN,
 480                         dcv,
 481                         seed);
 482                 // Close output files
 483                 fout1.flush();
 484                 int64_t tellpSz1 = (int64_t)fout1.tellp();
 485                 VMSG_NL("Wrote " << fout1.tellp() << " bytes to primary EBWT file: " << _in1Str);
 486                 fout1.close();
 487                 bool err = false;
 488                 if(tellpSz1 > fileSize(_in1Str.c_str())) {
 489                         err = true;
 490                         cerr << "Index is corrupt: File size for " << _in1Str << " should have been " << tellpSz1
 491                              << " but is actually " << fileSize(_in1Str.c_str()) << "." << endl;
 492                 }
 493                 fout2.flush();
 494                 int64_t tellpSz2 = (int64_t)fout2.tellp();
 495                 VMSG_NL("Wrote " << fout2.tellp() << " bytes to secondary EBWT file: " << _in2Str);
 496                 fout2.close();
 497                 if(tellpSz2 > fileSize(_in2Str.c_str())) {
 498                         err = true;
 499                         cerr << "Index is corrupt: File size for " << _in2Str << " should have been " << tellpSz2
 500                              << " but is actually " << fileSize(_in2Str.c_str()) << "." << endl;
 501                 }
 502                 if(err) {
 503                         cerr << "Please check if there is a problem with the disk or if disk is full." << endl;
 504                         throw 1;
 505                 }
 506                 // Reopen as input streams
 507                 VMSG_NL("Re-opening _in1 and _in2 as input streams");
 508                 if(_sanity) {
 509                         VMSG_NL("Sanity-checking Ebwt");
 510                         assert(!isInMemory());
 511                         readIntoMemory(
 512                                 color,
 513                                 __fw ? -1 : refparams.reverse == REF_READ_REVERSE,
 514                                 false,
 515                                 NULL,
 516                                 false,
 517                                 true,
 518                                 false);
 519                         sanityCheckAll(refparams.reverse);
 520                         evictFromMemory();
 521                         assert(!isInMemory());
 522                 }
 523                 VMSG_NL("Returning from Ebwt constructor");
 524         }
 525
 526         bool isPacked();
 527
 528         /**
 529          * Write the rstarts array given the szs array for the reference.
 530          */
 531         void szsToDisk(const vector<RefRecord>& szs, ostream& os, int reverse) {
 532                 size_t seq = 0;
 533                 uint32_t off = 0;
 534                 uint32_t totlen = 0;
 535                 for(unsigned int i = 0; i < szs.size(); i++) {
 536                         if(szs[i].len == 0) continue;
 537                         if(szs[i].first) off = 0;
 538                         off += szs[i].off;
 539 #ifdef ACCOUNT_FOR_ALL_GAP_REFS
 540                         if(szs[i].first && szs[i].len > 0) seq++;
 541 #else
 542                         if(szs[i].first) seq++;
 543 #endif
 544                         size_t seqm1 = seq-1;
 545                         assert_lt(seqm1, _nPat);
 546                         size_t fwoff = off;
 547                         if(reverse == REF_READ_REVERSE) {
 548                                 // Invert pattern idxs
 549                                 seqm1 = _nPat - seqm1 - 1;
 550                                 // Invert pattern idxs
 551                                 assert_leq(off + szs[i].len, _plen[seqm1]);
 552                                 fwoff = _plen[seqm1] - (off + szs[i].len);
 553                         }
 554                         writeU32(os, totlen, this->toBe()); // offset from beginning of joined string
 555                         writeU32(os, seqm1,  this->toBe()); // sequence id
 556                         writeU32(os, fwoff,  this->toBe()); // offset into sequence
 557                         totlen += szs[i].len;
 558                         off += szs[i].len;
 559                 }
 560         }
 561
 562         /**
 563          * Helper for the constructors above.  Takes a vector of text
 564          * strings and joins them into a single string with a call to
 565          * joinToDisk, which does a join (with padding) and writes some of
 566          * the resulting data directly to disk rather than keep it in
 567          * memory.  It then constructs a suffix-array producer (what kind
 568          * depends on 'useBlockwise') for the resulting sequence.  The
 569          * suffix-array producer can then be used to obtain chunks of the
 570          * joined string's suffix array.
 571          */
 572         void initFromVector(
 573                 vector<FileBuf*>& is,
 574                 vector<RefRecord>& szs,
 575                 vector<uint32_t>& plens,
 576                 uint32_t sztot,
 577                 const RefReadInParams& refparams,
 578                 ofstream& out1,
 579                 ofstream& out2,
 580                 bool useBlockwise,
 581                 uint32_t bmax,
 582                 uint32_t bmaxSqrtMult,
 583                 uint32_t bmaxDivN,
 584                 int dcv,
 585                 uint32_t seed)
 586         {
 587                 // Compose text strings into single string
 588                 VMSG_NL("Calculating joined length");
 589                 TStr s; // holds the entire joined reference after call to joinToDisk
 590                 uint32_t jlen;
 591                 jlen = joinedLen(szs);
 592                 assert_geq(jlen, sztot);
 593                 VMSG_NL("Writing header");
 594                 writeFromMemory(true, out1, out2);
 595                 try {
 596                         VMSG_NL("Reserving space for joined string");
 597                         seqan::reserve(s, jlen, Exact());
 598                         VMSG_NL("Joining reference sequences");
 599                         if(refparams.reverse == REF_READ_REVERSE) {
 600                                 {
 601                                         Timer timer(cout, "  Time to join reference sequences: ", _verbose);
 602                                         joinToDisk(is, szs, plens, sztot, refparams, s, out1, out2, seed);
 603                                 } {
 604                                         Timer timer(cout, "  Time to reverse reference sequence: ", _verbose);
 605                                         vector<RefRecord> tmp;
 606                                         reverseInPlace(s);
 607                                         reverseRefRecords(szs, tmp, false, false);
 608                                         szsToDisk(tmp, out1, refparams.reverse);
 609                                 }
 610                         } else {
 611                                 Timer timer(cout, "  Time to join reference sequences: ", _verbose);
 612                                 joinToDisk(is, szs, plens, sztot, refparams, s, out1, out2, seed);
 613                                 szsToDisk(szs, out1, refparams.reverse);
 614                         }
 615                         // Joined reference sequence now in 's'
 616                 } catch(bad_alloc& e) {
 617                         // If we throw an allocation exception in the try block,
 618                         // that means that the joined version of the reference
 619                         // string itself is too larger to fit in memory.  The only
 620                         // alternatives are to tell the user to give us more memory
 621                         // or to try again with a packed representation of the
 622                         // reference (if we haven't tried that already).
 623                         cerr << "Could not allocate space for a joined string of " << jlen << " elements." << endl;
 624                         if(!isPacked() && _passMemExc) {
 625                                 // Pass the exception up so that we can retry using a
 626                                 // packed string representation
 627                                 throw e;
 628                         }
 629                         // There's no point passing this exception on.  The fact
 630                         // that we couldn't allocate the joined string means that
 631                         // --bmax is irrelevant - the user should re-run with
 632                         // ebwt-build-packed
 633                         if(isPacked()) {
 634                                 cerr << "Please try running bowtie-build on a computer with more memory." << endl;
 635                         } else {
 636                                 cerr << "Please try running bowtie-build in packed mode (-p/--packed) or in automatic" << endl
 637                                      << "mode (-a/--auto), or try again on a computer with more memory." << endl;
 638                         }
 639                         if(sizeof(void*) == 4) {
 640                                 cerr << "If this computer has more than 4 GB of memory, try using a 64-bit executable;" << endl
 641                                      << "this executable is 32-bit." << endl;
 642                         }
 643                         throw 1;
 644                 }
 645                 // Succesfully obtained joined reference string
 646                 assert_geq(length(s), jlen);
 647                 if(bmax != 0xffffffff) {
 648                         VMSG_NL("bmax according to bmax setting: " << bmax);
 649                 }
 650                 else if(bmaxSqrtMult != 0xffffffff) {
 651                         bmax *= bmaxSqrtMult;
 652                         VMSG_NL("bmax according to bmaxSqrtMult setting: " << bmax);
 653                 }
 654                 else if(bmaxDivN != 0xffffffff) {
 655                         bmax = max<uint32_t>(jlen / bmaxDivN, 1);
 656                         VMSG_NL("bmax according to bmaxDivN setting: " << bmax);
 657                 }
 658                 else {
 659                         bmax = (uint32_t)sqrt(length(s));
 660                         VMSG_NL("bmax defaulted to: " << bmax);
 661                 }
 662                 int iter = 0;
 663                 bool first = true;
 664                 // Look for bmax/dcv parameters that work.
 665                 while(true) {
 666                         if(!first && bmax < 40 && _passMemExc) {
 667                                 cerr << "Could not find approrpiate bmax/dcv settings for building this index." << endl;
 668                                 if(!isPacked()) {
 669                                         // Throw an exception exception so that we can
 670                                         // retry using a packed string representation
 671                                         throw bad_alloc();
 672                                 } else {
 673                                         cerr << "Already tried a packed string representation." << endl;
 674                                 }
 675                                 cerr << "Please try indexing this reference on a computer with more memory." << endl;
 676                                 if(sizeof(void*) == 4) {
 677                                         cerr << "If this computer has more than 4 GB of memory, try using a 64-bit executable;" << endl
 678                                                  << "this executable is 32-bit." << endl;
 679                                 }
 680                                 throw 1;
 681                         }
 682                         if(dcv > 4096) dcv = 4096;
 683                         if((iter % 6) == 5 && dcv < 4096 && dcv != 0) {
 684                                 dcv <<= 1; // double difference-cover period
 685                         } else {
 686                                 bmax -= (bmax >> 2); // reduce by 25%
 687                         }
 688                         VMSG("Using parameters --bmax " << bmax);
 689                         if(dcv == 0) {
 690                                 VMSG_NL(" and *no difference cover*");
 691                         } else {
 692                                 VMSG_NL(" --dcv " << dcv);
 693                         }
 694                         iter++;
 695                         try {
 696                                 {
 697                                         VMSG_NL("  Doing ahead-of-time memory usage test");
 698                                         // Make a quick-and-dirty attempt to force a bad_alloc iff
 699                                         // we would have thrown one eventually as part of
 700                                         // constructing the DifferenceCoverSample
 701                                         dcv <<= 1;
 702                                         size_t sz = DifferenceCoverSample<TStr>::simulateAllocs(s, dcv >> 1);
 703                                         AutoArray<uint8_t> tmp(sz);
 704                                         dcv >>= 1;
 705                                         // Likewise with the KarkkainenBlockwiseSA
 706                                         sz = KarkkainenBlockwiseSA<TStr>::simulateAllocs(s, bmax);
 707                                         AutoArray<uint8_t> tmp2(sz);
 708                                         // Now throw in the 'ftab' and 'isaSample' structures
 709                                         // that we'll eventually allocate in buildToDisk
 710                                         AutoArray<uint32_t> ftab(_eh._ftabLen * 2);
 711                                         AutoArray<uint8_t> side(_eh._sideSz);
 712                                         // Grab another 20 MB out of caution
 713                                         AutoArray<uint32_t> extra(20*1024*1024);
 714                                         // If we made it here without throwing bad_alloc, then we
 715                                         // passed the memory-usage stress test
 716                                         VMSG("  Passed!  Constructing with these parameters: --bmax " << bmax << " --dcv " << dcv);
 717                                         if(isPacked()) {
 718                                                 VMSG(" --packed");
 719                                         }
 720                                         VMSG_NL("");
 721                                 }
 722                                 VMSG_NL("Constructing suffix-array element generator");
 723                                 KarkkainenBlockwiseSA<TStr> bsa(s, bmax, dcv, seed, _sanity, _passMemExc, _verbose);
 724                                 assert(bsa.suffixItrIsReset());
 725                                 assert_eq(bsa.size(), length(s)+1);
 726                                 VMSG_NL("Converting suffix-array elements to index image");
 727                                 buildToDisk(bsa, s, out1, out2);
 728                                 out1.flush(); out2.flush();
 729                                 if(out1.fail() || out2.fail()) {
 730                                         cerr << "An error occurred writing the index to disk.  Please check if the disk is full." << endl;
 731                                         throw 1;
 732                                 }
 733                                 break;
 734                         } catch(bad_alloc& e) {
 735                                 if(_passMemExc) {
 736                                         VMSG_NL("  Ran out of memory; automatically trying more memory-economical parameters.");
 737                                 } else {
 738                                         cerr << "Out of memory while constructing suffix array.  Please try using a smaller" << endl
 739                                                  << "number of blocks by specifying a smaller --bmax or a larger --bmaxdivn" << endl;
 740                                         throw 1;
 741                                 }
 742                         }
 743                         first = false;
 744                 }
 745                 assert(repOk());
 746                 // Now write reference sequence names on the end
 747 #ifdef ACCOUNT_FOR_ALL_GAP_REFS
 748                 assert_geq(this->_refnames.size(), this->_nPat);
 749 #else
 750                 assert_eq(this->_refnames.size(), this->_nPat);
 751 #endif
 752                 for(size_t i = 0; i < this->_refnames.size(); i++) {
 753                         out1 << this->_refnames[i] << endl;
 754                 }
 755                 out1 << '\0';
 756                 out1.flush(); out2.flush();
 757                 if(out1.fail() || out2.fail()) {
 758                         cerr << "An error occurred writing the index to disk.  Please check if the disk is full." << endl;
 759                         throw 1;
 760                 }
 761                 VMSG_NL("Returning from initFromVector");
 762         }
 763
 764         /**
 765          * Return the length that the joined string of the given string
 766          * list will have.  Note that this is indifferent to how the text
 767          * fragments correspond to input sequences - it just cares about
 768          * the lengths of the fragments.
 769          */
 770         uint32_t joinedLen(vector<RefRecord>& szs) {
 771                 uint32_t ret = 0;
 772                 for(unsigned int i = 0; i < szs.size(); i++) {
 773                         ret += szs[i].len;
 774                 }
 775                 return ret;
 776         }
 777
 778         /// Destruct an Ebwt
 779         ~Ebwt() {
 780                 // Only free buffers if we're *not* using memory-mapped files
 781                 if(!_useMm) {
 782                         // Delete everything that was allocated in read(false, ...)
 783                         if(_fchr    != NULL) delete[] _fchr;    _fchr    = NULL;
 784                         if(_ftab    != NULL) delete[] _ftab;    _ftab    = NULL;
 785                         if(_eftab   != NULL) delete[] _eftab;   _eftab   = NULL;
 786                         if(_offs != NULL && !useShmem_) {
 787                                 delete[] _offs; _offs = NULL;
 788                         } else if(_offs != NULL && useShmem_) {
 789                                 FREE_SHARED(_offs);
 790                         }
 791                         if(_isa     != NULL) delete[] _isa;     _isa     = NULL;
 792                         if(_plen    != NULL) delete[] _plen;    _plen    = NULL;
 793                         if(_rstarts != NULL) delete[] _rstarts; _rstarts = NULL;
 794                         if(_ebwt != NULL && !useShmem_) {
 795                                 delete[] _ebwt; _ebwt = NULL;
 796                         } else if(_ebwt != NULL && useShmem_) {
 797                                 FREE_SHARED(_ebwt);
 798                         }
 799                 }
 800                 MM_FILE_CLOSE(_in1);
 801                 MM_FILE_CLOSE(_in2);
 802 #ifdef EBWT_STATS
 803                 cout << (_fw ? "Forward index:" : "Mirror index:") << endl;
 804                 cout << "  mapLFEx:   " << mapLFExs_ << endl;
 805                 cout << "  mapLF:     " << mapLFs_   << endl;
 806                 cout << "  mapLF(c):  " << mapLFcs_  << endl;
 807                 cout << "  mapLF1(c): " << mapLF1cs_ << endl;
 808                 cout << "  mapLF(c):  " << mapLF1s_  << endl;
 809 #endif
 810         }
 811
 812         /// Accessors
 813         const EbwtParams& eh() const     { return _eh; }
 814         uint32_t    zOff() const         { return _zOff; }
 815         uint32_t    zEbwtByteOff() const { return _zEbwtByteOff; }
 816         int         zEbwtBpOff() const   { return _zEbwtBpOff; }
 817         uint32_t    nPat() const         { return _nPat; }
 818         uint32_t    nFrag() const        { return _nFrag; }
 819         uint32_t*   fchr() const         { return _fchr; }
 820         uint32_t*   ftab() const         { return _ftab; }
 821         uint32_t*   eftab() const        { return _eftab; }
 822         uint32_t*   offs() const         { return _offs; }
 823         uint32_t*   isa() const          { return _isa; }
 824         uint32_t*   plen() const         { return _plen; }
 825         uint32_t*   rstarts() const      { return _rstarts; }
 826         uint8_t*    ebwt() const         { return _ebwt; }
 827         const ReferenceMap* rmap() const { return rmap_; }
 828         bool        toBe() const         { return _toBigEndian; }
 829         bool        verbose() const      { return _verbose; }
 830         bool        sanityCheck() const  { return _sanity; }
 831         vector<string>& refnames()       { return _refnames; }
 832         bool        fw() const           { return _fw; }
 833
 834         /// Return true iff the Ebwt is currently in memory
 835         bool isInMemory() const {
 836                 if(_ebwt != NULL) {
 837                         assert(_eh.repOk());
 838                         assert(_ftab != NULL);
 839                         assert(_eftab != NULL);
 840                         assert(_fchr != NULL);
 841                         assert(_offs != NULL);
 842                         assert(_isa != NULL);
 843                         assert(_rstarts != NULL);
 844                         assert_neq(_zEbwtByteOff, 0xffffffff);
 845                         assert_neq(_zEbwtBpOff, -1);
 846                         return true;
 847                 } else {
 848                         assert(_ftab == NULL);
 849                         assert(_eftab == NULL);
 850                         assert(_fchr == NULL);
 851                         assert(_offs == NULL);
 852                         assert(_rstarts == NULL);
 853                         assert_eq(_zEbwtByteOff, 0xffffffff);
 854                         assert_eq(_zEbwtBpOff, -1);
 855                         return false;
 856                 }
 857         }
 858
 859         /// Return true iff the Ebwt is currently stored on disk
 860         bool isEvicted() const {
 861                 return !isInMemory();
 862         }
 863
 864         /**
 865          * Load this Ebwt into memory by reading it in from the _in1 and
 866          * _in2 streams.
 867          */
 868         void loadIntoMemory(
 869                 int color,
 870                 int needEntireReverse,
 871                 bool loadNames,
 872                 bool verbose)
 873         {
 874                 readIntoMemory(
 875                         color,      // expect index to be colorspace?
 876                         needEntireReverse, // require reverse index to be concatenated reference reversed
 877                         false,      // stop after loading the header portion?
 878                         NULL,       // params
 879                         false,      // mmSweep
 880                         loadNames,  // loadNames
 881                         verbose);   // startVerbose
 882         }
 883
 884         /**
 885          * Frees memory associated with the Ebwt.
 886          */
 887         void evictFromMemory() {
 888                 assert(isInMemory());
 889                 if(!_useMm) {
 890                         delete[] _fchr;
 891                         delete[] _ftab;
 892                         delete[] _eftab;
 893                         if(!useShmem_) delete[] _offs;
 894                         delete[] _isa;
 895                         // Keep plen; it's small and the client may want to query it
 896                         // even when the others are evicted.
 897                         //delete[] _plen;
 898                         delete[] _rstarts;
 899                         if(!useShmem_) delete[] _ebwt;
 900                 }
 901                 _fchr  = NULL;
 902                 _ftab  = NULL;
 903                 _eftab = NULL;
 904                 _offs  = NULL;
 905                 _isa   = NULL;
 906                 // Keep plen; it's small and the client may want to query it
 907                 // even when the others are evicted.
 908                 //_plen  = NULL;
 909                 _rstarts = NULL;
 910                 _ebwt    = NULL;
 911                 _zEbwtByteOff = 0xffffffff;
 912                 _zEbwtBpOff = -1;
 913         }
 914
 915         /**
 916          * Non-static facade for static function ftabHi.
 917          */
 918         uint32_t ftabHi(uint32_t i) const {
 919                 return Ebwt::ftabHi(_ftab, _eftab, _eh._len, _eh._ftabLen,
 920                                     _eh._eftabLen, i);
 921         }
 922
 923         /**
 924          * Get "high interpretation" of ftab entry at index i.  The high
 925          * interpretation of a regular ftab entry is just the entry
 926          * itself.  The high interpretation of an extended entry is the
 927          * second correpsonding ui32 in the eftab.
 928          *
 929          * It's a static member because it's convenient to ask this
 930          * question before the Ebwt is fully initialized.
 931          */
 932         static uint32_t ftabHi(uint32_t *ftab,
 933                                uint32_t *eftab,
 934                                uint32_t len,
 935                                uint32_t ftabLen,
 936                                uint32_t eftabLen,
 937                                uint32_t i)
 938         {
 939                 assert_lt(i, ftabLen);
 940                 if(ftab[i] <= len) {
 941                         return ftab[i];
 942                 } else {
 943                         uint32_t efIdx = ftab[i] ^ 0xffffffff;
 944                         assert_lt(efIdx*2+1, eftabLen);
 945                         return eftab[efIdx*2+1];
 946                 }
 947         }
 948
 949         /**
 950          * Non-static facade for static function ftabLo.
 951          */
 952         uint32_t ftabLo(uint32_t i) const {
 953                 return Ebwt::ftabLo(_ftab, _eftab, _eh._len, _eh._ftabLen,
 954                                     _eh._eftabLen, i);
 955         }
 956
 957         /**
 958          * Get "low interpretation" of ftab entry at index i.  The low
 959          * interpretation of a regular ftab entry is just the entry
 960          * itself.  The low interpretation of an extended entry is the
 961          * first correpsonding ui32 in the eftab.
 962          *
 963          * It's a static member because it's convenient to ask this
 964          * question before the Ebwt is fully initialized.
 965          */
 966         static uint32_t ftabLo(uint32_t *ftab,
 967                                uint32_t *eftab,
 968                                uint32_t len,
 969                                uint32_t ftabLen,
 970                                uint32_t eftabLen,
 971                                uint32_t i)
 972         {
 973                 assert_lt(i, ftabLen);
 974                 if(ftab[i] <= len) {
 975                         return ftab[i];
 976                 } else {
 977                         uint32_t efIdx = ftab[i] ^ 0xffffffff;
 978                         assert_lt(efIdx*2+1, eftabLen);
 979                         return eftab[efIdx*2];
 980                 }
 981         }
 982
 983         /**
 984          * When using read() to create an Ebwt, we have to set a couple of
 985          * additional fields in the Ebwt object that aren't part of the
 986          * parameter list and are not stored explicitly in the file.  Right
 987          * now, this just involves initializing _zEbwtByteOff and
 988          * _zEbwtBpOff from _zOff.
 989          */
 990         void postReadInit(EbwtParams& eh) {
 991                 uint32_t sideNum     = _zOff / eh._sideBwtLen;
 992                 uint32_t sideCharOff = _zOff % eh._sideBwtLen;
 993                 uint32_t sideByteOff = sideNum * eh._sideSz;
 994                 _zEbwtByteOff = sideCharOff >> 2;
 995                 assert_lt(_zEbwtByteOff, eh._sideBwtSz);
 996                 _zEbwtBpOff = sideCharOff & 3;
 997                 assert_lt(_zEbwtBpOff, 4);
 998                 if((sideNum & 1) == 0) {
 999                         // This is an even (backward) side
1000                         _zEbwtByteOff = eh._sideBwtSz - _zEbwtByteOff - 1;
1001                         _zEbwtBpOff = 3 - _zEbwtBpOff;
1002                         assert_lt(_zEbwtBpOff, 4);
1003                 }
1004                 _zEbwtByteOff += sideByteOff;
1005                 assert(repOk(eh)); // Ebwt should be fully initialized now
1006         }
1007
1008         /**
1009          * Pretty-print the Ebwt to the given output stream.
1010          */
1011         void print(ostream& out) const {
1012                 print(out, _eh);
1013         }
1014
1015         /**
1016          * Pretty-print the Ebwt and given EbwtParams to the given output
1017          * stream.
1018          */
1019         void print(ostream& out, const EbwtParams& eh) const {
1020                 eh.print(out); // print params
1021                 out << "Ebwt (" << (isInMemory()? "memory" : "disk") << "):" << endl
1022                     << "    zOff: "         << _zOff << endl
1023                     << "    zEbwtByteOff: " << _zEbwtByteOff << endl
1024                     << "    zEbwtBpOff: "   << _zEbwtBpOff << endl
1025                     << "    nPat: "  << _nPat << endl
1026                     << "    plen: ";
1027                 if(_plen == NULL) {
1028                         out << "NULL" << endl;
1029                 } else {
1030                         out << "non-NULL, [0] = " << _plen[0] << endl;
1031                 }
1032                 out << "    rstarts: ";
1033                 if(_rstarts == NULL) {
1034                         out << "NULL" << endl;
1035                 } else {
1036                         out << "non-NULL, [0] = " << _rstarts[0] << endl;
1037                 }
1038                 out << "    ebwt: ";
1039                 if(_ebwt == NULL) {
1040                         out << "NULL" << endl;
1041                 } else {
1042                         out << "non-NULL, [0] = " << _ebwt[0] << endl;
1043                 }
1044                 out << "    fchr: ";
1045                 if(_fchr == NULL) {
1046                         out << "NULL" << endl;
1047                 } else {
1048                         out << "non-NULL, [0] = " << _fchr[0] << endl;
1049                 }
1050                 out << "    ftab: ";
1051                 if(_ftab == NULL) {
1052                         out << "NULL" << endl;
1053                 } else {
1054                         out << "non-NULL, [0] = " << _ftab[0] << endl;
1055                 }
1056                 out << "    eftab: ";
1057                 if(_eftab == NULL) {
1058                         out << "NULL" << endl;
1059                 } else {
1060                         out << "non-NULL, [0] = " << _eftab[0] << endl;
1061                 }
1062                 out << "    offs: ";
1063                 if(_offs == NULL) {
1064                         out << "NULL" << endl;
1065                 } else {
1066                         out << "non-NULL, [0] = " << _offs[0] << endl;
1067                 }
1068         }
1069
1070         // Building
1071         static TStr join(vector<TStr>& l, uint32_t seed);
1072         static TStr join(vector<FileBuf*>& l, vector<RefRecord>& szs, uint32_t sztot, const RefReadInParams& refparams, uint32_t seed);
1073         void joinToDisk(vector<FileBuf*>& l, vector<RefRecord>& szs, vector<uint32_t>& plens, uint32_t sztot, const RefReadInParams& refparams, TStr& ret, ostream& out1, ostream& out2, uint32_t seed);
1074         void buildToDisk(InorderBlockwiseSA<TStr>& sa, const TStr& s, ostream& out1, ostream& out2);
1075
1076         // I/O
1077         void readIntoMemory(int color, int needEntireReverse, bool justHeader, EbwtParams *params, bool mmSweep, bool loadNames, bool startVerbose);
1078         void writeFromMemory(bool justHeader, ostream& out1, ostream& out2) const;
1079         void writeFromMemory(bool justHeader, const string& out1, const string& out2) const;
1080
1081         // Sanity checking
1082         void printRangeFw(uint32_t begin, uint32_t end) const;
1083         void printRangeBw(uint32_t begin, uint32_t end) const;
1084         void sanityCheckUpToSide(int upToSide) const;
1085         void sanityCheckAll(int reverse) const;
1086         void restore(TStr& s) const;
1087         void checkOrigs(const vector<String<Dna5> >& os, bool color, bool mirror) const;
1088
1089         // Searching and reporting
1090         void joinedToTextOff(uint32_t qlen, uint32_t off, uint32_t& tidx, uint32_t& textoff, uint32_t& tlen) const;
1091         inline bool report(const String<Dna5>& query, String<char>* quals, String<char>* name, bool color, bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector<uint32_t>& mmui32, const std::vector<uint8_t>& refcs, size_t numMms, uint32_t off, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, uint32_t patid, uint32_t seed, const EbwtSearchParams<TStr>& params) const;
1092         inline bool reportChaseOne(const String<Dna5>& query, String<char>* quals, String<char>* name, bool color, bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector<uint32_t>& mmui32, const std::vector<uint8_t>& refcs, size_t numMms, uint32_t i, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, uint32_t patid, uint32_t seed, const EbwtSearchParams<TStr>& params, SideLocus *l = NULL) const;
1093         inline bool reportReconstruct(const String<Dna5>& query, String<char>* quals, String<char>* name, String<Dna5>& lbuf, String<Dna5>& rbuf, const uint32_t *mmui32, const char* refcs, size_t numMms, uint32_t i, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, const EbwtSearchParams<TStr>& params, SideLocus *l = NULL) const;
1094         inline int rowL(const SideLocus& l) const;
1095         inline uint32_t countUpTo(const SideLocus& l, int c) const;
1096         inline void countUpToEx(const SideLocus& l, uint32_t* pairs) const;
1097         inline uint32_t countFwSide(const SideLocus& l, int c) const;
1098         inline void countFwSideEx(const SideLocus& l, uint32_t *pairs) const;
1099         inline uint32_t countBwSide(const SideLocus& l, int c) const;
1100         inline void countBwSideEx(const SideLocus& l, uint32_t *pairs) const;
1101         inline uint32_t mapLF(const SideLocus& l ASSERT_ONLY(, bool overrideSanity = false)) const;
1102         inline void mapLFEx(const SideLocus& l, uint32_t *pairs ASSERT_ONLY(, bool overrideSanity = false)) const;
1103         inline void mapLFEx(const SideLocus& ltop, const SideLocus& lbot, uint32_t *tops, uint32_t *bots ASSERT_ONLY(, bool overrideSanity = false)) const;
1104         inline uint32_t mapLF(const SideLocus& l, int c ASSERT_ONLY(, bool overrideSanity = false)) const;
1105         inline uint32_t mapLF1(uint32_t row, const SideLocus& l, int c ASSERT_ONLY(, bool overrideSanity = false)) const;
1106         inline int mapLF1(uint32_t& row, const SideLocus& l ASSERT_ONLY(, bool overrideSanity = false)) const;
1107         /// Check that in-memory Ebwt is internally consistent with respect
1108         /// to given EbwtParams; assert if not
1109         bool inMemoryRepOk(const EbwtParams& eh) const {
1110                 assert_leq(ValueSize<TAlphabet>::VALUE, 4);
1111                 assert_geq(_zEbwtBpOff, 0);
1112                 assert_lt(_zEbwtBpOff, 4);
1113                 assert_lt(_zEbwtByteOff, eh._ebwtTotSz);
1114                 assert_lt(_zOff, eh._bwtLen);
1115                 assert(_rstarts != NULL);
1116                 assert_geq(_nFrag, _nPat);
1117                 return true;
1118         }
1119
1120         /// Check that in-memory Ebwt is internally consistent; assert if
1121         /// not
1122         bool inMemoryRepOk() const {
1123                 return repOk(_eh);
1124         }
1125
1126         /// Check that Ebwt is internally consistent with respect to given
1127         /// EbwtParams; assert if not
1128         bool repOk(const EbwtParams& eh) const {
1129                 assert(_eh.repOk());
1130                 if(isInMemory()) {
1131                         return inMemoryRepOk(eh);
1132                 }
1133                 return true;
1134         }
1135
1136         /// Check that Ebwt is internally consistent; assert if not
1137         bool repOk() const {
1138                 return repOk(_eh);
1139         }
1140
1141         bool       _toBigEndian;
1142         int32_t    _overrideOffRate;
1143         int32_t    _overrideIsaRate;
1144         bool       _verbose;
1145         bool       _passMemExc;
1146         bool       _sanity;
1147         bool       _fw;     // true iff this is a forward index
1148         MM_FILE  _in1;    // input fd for primary index file
1149         MM_FILE  _in2;    // input fd for secondary index file
1150         string     _in1Str; // filename for primary index file
1151         string     _in2Str; // filename for secondary index file
1152         uint32_t   _zOff;
1153         uint32_t   _zEbwtByteOff;
1154         int        _zEbwtBpOff;
1155         uint32_t   _nPat;  /// number of reference texts
1156         uint32_t   _nFrag; /// number of fragments
1157         uint32_t*  _plen;
1158         uint32_t*  _rstarts; // starting offset of fragments / text indexes
1159         // _fchr, _ftab and _eftab are expected to be relatively small
1160         // (usually < 1MB, perhaps a few MB if _fchr is particularly large
1161         // - like, say, 11).  For this reason, we don't bother with writing
1162         // them to disk through separate output streams; we
1163         uint32_t*  _fchr;
1164         uint32_t*  _ftab;
1165         uint32_t*  _eftab; // "extended" entries for _ftab
1166         // _offs may be extremely large.  E.g. for DNA w/ offRate=4 (one
1167         // offset every 16 rows), the total size of _offs is the same as
1168         // the total size of the input sequence
1169         uint32_t*  _offs;
1170         uint32_t*  _isa;
1171         // _ebwt is the Extended Burrows-Wheeler Transform itself, and thus
1172         // is at least as large as the input sequence.
1173         uint8_t*   _ebwt;
1174         bool       _useMm;        /// use memory-mapped files to hold the index
1175         bool       useShmem_;     /// use shared memory to hold large parts of the index
1176         vector<string> _refnames; /// names of the reference sequences
1177         const ReferenceMap* rmap_; /// mapping into another reference coordinate space
1178         char *mmFile1_;
1179         char *mmFile2_;
1180         EbwtParams _eh;
1181
1182 #ifdef EBWT_STATS
1183         uint64_t   mapLFExs_;
1184         uint64_t   mapLFs_;
1185         uint64_t   mapLFcs_;
1186 #endif
1187
1188 private:
1189
1190         ostream& log() const {
1191                 return cout; // TODO: turn this into a parameter
1192         }
1193
1194         /// Print a verbose message and flush (flushing is helpful for
1195         /// debugging)
1196         void verbose(const string& s) const {
1197                 if(this->verbose()) {
1198                         this->log() << s;
1199                         this->log().flush();
1200                 }
1201         }
1202 };
1203
1204 /// Specialization for packed Ebwts - return true
1205 template<>
1206 bool Ebwt<String<Dna, Packed<> > >::isPacked() {
1207         return true;
1208 }
1209
1210 /// By default, Ebwts are not packed
1211 template<typename TStr>
1212 bool Ebwt<TStr>::isPacked() {
1213         return false;
1214 }
1215
1216 /**
1217  * Structure encapsulating search parameters, such as whether and how
1218  * to backtrack and how to deal with multiple equally-good hits.
1219  */
1220 template<typename TStr>
1221 class EbwtSearchParams {
1222 public:
1223         EbwtSearchParams(HitSinkPerThread& sink,
1224                          const vector<String<Dna5> >& texts,
1225                          bool fw = true,
1226                          bool ebwtFw = true) :
1227                 _sink(sink),
1228                 _texts(texts),
1229                 _patid(0xffffffff),
1230                 _fw(fw) { }
1231
1232         HitSinkPerThread& sink() const { return _sink; }
1233         void setPatId(uint32_t patid)  { _patid = patid; }
1234         uint32_t patId() const         { return _patid; }
1235         void setFw(bool fw)            { _fw = fw; }
1236         bool fw() const                { return _fw; }
1237         /**
1238          * Report a hit.  Returns true iff caller can call off the search.
1239          */
1240         bool reportHit(const String<Dna5>& query, // read sequence
1241                        String<char>* quals, // read quality values
1242                        String<char>* name,  // read name
1243                        bool color,          // true -> read is colorspace
1244                        bool colExEnds,      // true -> exclude nucleotides at extreme ends after decoding
1245                        int snpPhred,        // penalty for a SNP
1246                        const BitPairReference* ref, // reference (= NULL if not necessary)
1247                        const ReferenceMap* rmap, // map to another reference coordinate system
1248                        bool ebwtFw,         // whether index is forward (true) or mirror (false)
1249                        const std::vector<uint32_t>& mmui32, // mismatch list
1250                        const std::vector<uint8_t>& refcs,  // reference characters
1251                        size_t numMms,      // # mismatches
1252                        U32Pair h,          // ref coords
1253                        U32Pair mh,         // mate's ref coords
1254                        bool mfw,           // mate's orientation
1255                        uint16_t mlen,      // mate length
1256                        U32Pair a,          // arrow pair
1257                        uint32_t tlen,      // length of text
1258                        uint32_t qlen,      // length of query
1259                        int stratum,        // alignment stratum
1260                        uint16_t cost,      // cost of alignment
1261                        uint32_t oms,       // approx. # other valid alignments
1262                        uint32_t patid,
1263                        uint32_t seed,
1264                        uint8_t mate) const
1265         {
1266 #ifndef NDEBUG
1267                 // Check that no two elements of the mms array are the same
1268                 for(size_t i = 0; i < numMms; i++) {
1269                         for(size_t j = i+1; j < numMms; j++) {
1270                                 assert_neq(mmui32[i], mmui32[j]);
1271                         }
1272                 }
1273 #endif
1274                 // If ebwtFw is true, then 'query' and 'quals' are reversed
1275                 // If _fw is false, then 'query' and 'quals' are reverse complemented
1276                 assert(!color || ref != NULL);
1277                 assert(quals != NULL);
1278                 assert(name != NULL);
1279                 assert_eq(mmui32.size(), refcs.size());
1280                 assert_leq(numMms, mmui32.size());
1281                 assert_gt(qlen, 0);
1282                 Hit hit;
1283                 hit.stratum = stratum;
1284                 hit.cost = cost;
1285                 hit.patSeq = query;
1286                 hit.quals = *quals;
1287                 if(!ebwtFw) {
1288                         // Re-reverse the pattern and the quals back to how they
1289                         // appeared in the read file
1290                         ::reverseInPlace(hit.patSeq);
1291                         ::reverseInPlace(hit.quals);
1292                 }
1293                 if(color) {
1294                         hit.colSeq = hit.patSeq;
1295                         hit.colQuals = hit.quals;
1296                         hit.crefcs.resize(qlen, 0);
1297                         // Turn the mmui32 and refcs arrays into the mm FixedBitset and
1298                         // the refc vector
1299                         for(size_t i = 0; i < numMms; i++) {
1300                                 if (ebwtFw != _fw) {
1301                                         // The 3' end is on the left but the mm vector encodes
1302                                         // mismatches w/r/t the 5' end, so we flip
1303                                         uint32_t off = qlen - mmui32[i] - 1;
1304                                         hit.cmms.set(off);
1305                                         hit.crefcs[off] = refcs[i];
1306                                 } else {
1307                                         hit.cmms.set(mmui32[i]);
1308                                         hit.crefcs[i] = refcs[i];
1309                                 }
1310                         }
1311                         assert(ref != NULL);
1312                         char read[1024];
1313                         uint32_t rfbuf[(1024+16)/4];
1314                         ASSERT_ONLY(char rfbuf2[1024]);
1315                         char qual[1024];
1316                         char ns[1024];
1317                         char cmm[1024];
1318                         char nmm[1024];
1319                         int cmms = 0;
1320                         int nmms = 0;
1321                         // TODO: account for indels when calculating these bounds
1322                         size_t readi = 0;
1323                         size_t readf = seqan::length(hit.patSeq);
1324                         size_t refi = 0;
1325                         size_t reff = readf + 1;
1326                         bool maqRound = false;
1327                         for(size_t i = 0; i < qlen + 1; i++) {
1328                                 if(i < qlen) {
1329                                         read[i] = (int)hit.patSeq[i];
1330                                         qual[i] = mmPenalty(maqRound, phredCharToPhredQual(hit.quals[i]));
1331                                 }
1332                                 ASSERT_ONLY(rfbuf2[i] = ref->getBase(h.first, h.second + i));
1333                         }
1334                         int offset = ref->getStretch(rfbuf, h.first, h.second, qlen + 1);
1335                         char *rf = (char*)rfbuf + offset;
1336                         for(size_t i = 0; i < qlen + 1; i++) {
1337                                 assert_eq(rf[i], rfbuf2[i]);
1338                                 rf[i] = (1 << rf[i]);
1339                         }
1340                         decodeHit(
1341                                 read,  // ASCII colors, '0', '1', '2', '3', '.'
1342                                 qual,  // ASCII quals, Phred+33 encoded
1343                                 readi, // offset of first character within 'read' to consider
1344                                 readf, // offset of last char (exclusive) in 'read' to consider
1345                                 rf,    // reference sequence, as masks
1346                                 refi, // offset of first character within 'ref' to consider
1347                                 reff, // offset of last char (exclusive) in 'ref' to consider
1348                                 snpPhred, // penalty incurred by a SNP
1349                                 ns,  // decoded nucleotides are appended here
1350                                 cmm, // where the color mismatches are in the string
1351                                 nmm, // where nucleotide mismatches are in the string
1352                                 cmms, // number of color mismatches
1353                                 nmms);// number of nucleotide mismatches
1354                         size_t nqlen = qlen + (colExEnds ? -1 : 1);
1355                         seqan::resize(hit.patSeq, nqlen);
1356                         seqan::resize(hit.quals, nqlen);
1357                         hit.refcs.resize(nqlen);
1358                         size_t lo = colExEnds ? 1 : 0;
1359                         size_t hi = colExEnds ? qlen : qlen+1;
1360                         size_t destpos = 0;
1361                         for(size_t i = lo; i < hi; i++, destpos++) {
1362                                 // Set sequence character
1363                                 assert_leq(ns[i], 4);
1364                                 assert_geq(ns[i], 0);
1365                                 hit.patSeq[destpos] = (Dna5)(int)ns[i];
1366                                 // Set initial quality
1367                                 hit.quals[destpos] = '!';
1368                                 // Color mismatches penalize quality
1369                                 if(i > 0) {
1370                                         if(cmm[i-1] == 'M') {
1371                                                 if((int)hit.quals[destpos] + (int)qual[i-1] > 126) {
1372                                                         hit.quals[destpos] = 126;
1373                                                 } else {
1374                                                         hit.quals[destpos] += qual[i-1];
1375                                                 }
1376                                         } else if((int)hit.colSeq[i-1] != 4) {
1377                                                 hit.quals[destpos] -= qual[i-1];
1378                                         }
1379                                 }
1380                                 if(i < qlen) {
1381                                         if(cmm[i] == 'M') {
1382                                                 if((int)hit.quals[destpos] + (int)qual[i] > 126) {
1383                                                         hit.quals[destpos] = 126;
1384                                                 } else {
1385                                                         hit.quals[destpos] += qual[i];
1386                                                 }
1387                                         } else if((int)hit.patSeq[i] != 4) {
1388                                                 hit.quals[destpos] -= qual[i];
1389                                         }
1390                                 }
1391                                 if(hit.quals[destpos] < '!') {
1392                                         hit.quals[destpos] = '!';
1393                                 }
1394                                 if(nmm[i] != 'M') {
1395                                         uint32_t off = i - (colExEnds? 1:0);
1396                                         if(!_fw) off = nqlen - off - 1;
1397                                         assert_lt(off, nqlen);
1398                                         hit.mms.set(off);
1399                                         hit.refcs[off] = "ACGT"[ref->getBase(h.first, h.second+i)];
1400                                 }
1401                         }
1402                         if(colExEnds) {
1403                                 // Extreme bases have been removed; that makes the
1404                                 // nucleotide alignment one character shorter than the
1405                                 // color alignment
1406                                 qlen--; mlen--;
1407                                 // It also shifts the alignment's offset up by 1
1408                                 h.second++;
1409                         } else {
1410                                 // Extreme bases are included; that makes the
1411                                 // nucleotide alignment one character longer than the
1412                                 // color alignment
1413                                 qlen++; mlen++;
1414                         }
1415                 } else {
1416                         // Turn the mmui32 and refcs arrays into the mm FixedBitset and
1417                         // the refc vector
1418                         hit.refcs.resize(qlen, 0);
1419                         for(size_t i = 0; i < numMms; i++) {
1420                                 if (ebwtFw != _fw) {
1421                                         // The 3' end is on the left but the mm vector encodes
1422                                         // mismatches w/r/t the 5' end, so we flip
1423                                         uint32_t off = qlen - mmui32[i] - 1;
1424                                         hit.mms.set(off);
1425                                         hit.refcs[off] = refcs[i];
1426                                 } else {
1427                                         hit.mms.set(mmui32[i]);
1428                                         hit.refcs[mmui32[i]] = refcs[i];
1429                                 }
1430                         }
1431                 }
1432                 // Check the hit against the original text, if it's available
1433                 if(_texts.size() > 0) {
1434                         assert_lt(h.first, _texts.size());
1435                         FixedBitset<1024> diffs;
1436                         // This type of check assumes that only mismatches are
1437                         // possible.  If indels are possible, then we either need
1438                         // the caller to provide information about indel locations,
1439                         // or we need to extend this to a more complicated check.
1440                         assert_leq(h.second + qlen, length(_texts[h.first]));
1441                         for(size_t i = 0; i < qlen; i++) {
1442                                 assert_neq(4, (int)_texts[h.first][h.second + i]);
1443                                 // Forward pattern appears at h
1444                                 if((int)hit.patSeq[i] != (int)_texts[h.first][h.second + i]) {
1445                                         uint32_t qoff = i;
1446                                         // if ebwtFw != _fw the 3' end is on on the
1447                                         // left end of the pattern, but the diff vector
1448                                         // should encode mismatches w/r/t the 5' end,
1449                                         // so we flip
1450                                         if (_fw) diffs.set(qoff);
1451                                         else     diffs.set(qlen - qoff - 1);
1452                                 }
1453                         }
1454                         if(diffs != hit.mms) {
1455                                 // Oops, mismatches were not where we expected them;
1456                                 // print a diagnostic message before asserting
1457                                 cerr << "Expected " << hit.mms.str() << " mismatches, got " << diffs.str() << endl;
1458                                 cerr << "  Pat:  " << hit.patSeq << endl;
1459                                 cerr << "  Tseg: ";
1460                                 for(size_t i = 0; i < qlen; i++) {
1461                                         cerr << _texts[h.first][h.second + i];
1462                                 }
1463                                 cerr << endl;
1464                                 cerr << "  mmui32: ";
1465                                 for(size_t i = 0; i < numMms; i++) {
1466                                         cerr << mmui32[i] << " ";
1467                                 }
1468                                 cerr << endl;
1469                                 cerr << "  FW: " << _fw << endl;
1470                                 cerr << "  Ebwt FW: " << ebwtFw << endl;
1471                         }
1472                         if(diffs != hit.mms) assert(false);
1473                 }
1474                 hit.h = h;
1475                 if(rmap != NULL) rmap->map(hit.h);
1476                 hit.patId = ((patid == 0xffffffff) ? _patid : patid);
1477                 hit.patName = *name;
1478                 hit.mh = mh;
1479                 hit.fw = _fw;
1480                 hit.mfw = mfw;
1481                 hit.mlen = mlen;
1482                 hit.oms = oms;
1483                 hit.mate = mate;
1484                 hit.color = color;
1485                 hit.seed = seed;
1486                 assert(hit.repOk());
1487                 return sink().reportHit(hit, stratum);
1488         }
1489 private:
1490         HitSinkPerThread& _sink;
1491         const vector<String<Dna5> >& _texts; // original texts, if available (if not
1492                                     // available, _texts.size() == 0)
1493         uint32_t _patid;      // id of current read
1494         bool _fw;             // current read is forward-oriented
1495 };
1496
1497 /**
1498  * Encapsulates a location in the bwt text in terms of the side it
1499  * occurs in and its offset within the side.
1500  */
1501 struct SideLocus {
1502         SideLocus() :
1503         _sideByteOff(0),
1504         _sideNum(0),
1505         _charOff(0),
1506         _fw(true),
1507         _by(-1),
1508         _bp(-1) { }
1509
1510         /**
1511          * Construct from row and other relevant information about the Ebwt.
1512          */
1513         SideLocus(uint32_t row, const EbwtParams& ep, const uint8_t* ebwt) {
1514                 initFromRow(row, ep, ebwt);
1515         }
1516
1517         /**
1518          * Init two SideLocus objects from a top/bot pair, using the result
1519          * from one call to initFromRow to possibly avoid a second call.
1520          */
1521         static void initFromTopBot(uint32_t top,
1522                                    uint32_t bot,
1523                                    const EbwtParams& ep,
1524                                    const uint8_t* ebwt,
1525                                    SideLocus& ltop,
1526                                    SideLocus& lbot)
1527         {
1528                 const uint32_t sideBwtLen = ep._sideBwtLen;
1529                 const uint32_t sideBwtSz  = ep._sideBwtSz;
1530                 assert_gt(bot, top);
1531                 ltop.initFromRow(top, ep, ebwt);
1532                 uint32_t spread = bot - top;
1533                 if(ltop._charOff + spread < sideBwtLen) {
1534                         lbot._charOff = ltop._charOff + spread;
1535                         lbot._sideNum = ltop._sideNum;
1536                         lbot._sideByteOff = ltop._sideByteOff;
1537                         lbot._fw = ltop._fw;
1538                         lbot._by = lbot._charOff >> 2;
1539                         assert_lt(lbot._by, (int)sideBwtSz);
1540                         if(!lbot._fw) lbot._by = sideBwtSz - lbot._by - 1;
1541                         lbot._bp = lbot._charOff & 3;
1542                         if(!lbot._fw) lbot._bp ^= 3;
1543                 } else {
1544                         lbot.initFromRow(bot, ep, ebwt);
1545                 }
1546         }
1547
1548         /**
1549          * Calculate SideLocus based on a row and other relevant
1550          * information about the shape of the Ebwt.
1551          */
1552         void initFromRow(uint32_t row, const EbwtParams& ep, const uint8_t* ebwt) {
1553                 const uint32_t sideSz     = ep._sideSz;
1554                 // Side length is hard-coded for now; this allows the compiler
1555                 // to do clever things to accelerate / and %.
1556                 _sideNum                  = row / 224;
1557                 _charOff                  = row % 224;
1558                 _sideByteOff              = _sideNum * sideSz;
1559                 assert_leq(row, ep._len);
1560                 assert_leq(_sideByteOff + sideSz, ep._ebwtTotSz);
1561 #ifndef NO_PREFETCH
1562                 __builtin_prefetch((const void *)(ebwt + _sideByteOff),
1563                                    0 /* prepare for read */,
1564                                    PREFETCH_LOCALITY);
1565 #endif
1566                 // prefetch tjside too
1567                 _fw = (_sideNum & 1) != 0; // odd-numbered sides are forward
1568                 _by = _charOff >> 2; // byte within side
1569                 assert_lt(_by, (int)ep._sideBwtSz);
1570                 _bp = _charOff & 3;  // bit-pair within byte
1571                 if(!_fw) {
1572                         _by = ep._sideBwtSz - _by - 1;
1573                         _bp ^= 3;
1574                 }
1575         }
1576
1577         /// Return true iff this is an initialized SideLocus
1578         bool valid() {
1579                 return _bp != -1;
1580         }
1581
1582         /// Make this look like an invalid SideLocus
1583         void invalidate() {
1584                 _bp = -1;
1585         }
1586
1587         const uint8_t *side(const uint8_t* ebwt) const {
1588                 return ebwt + _sideByteOff;
1589         }
1590
1591         const uint8_t *oside(const uint8_t* ebwt) const {
1592                 return ebwt + _sideByteOff + (_fw? (-128) : (128));
1593         }
1594
1595         uint32_t _sideByteOff; // offset of top side within ebwt[]
1596         uint32_t _sideNum;     // index of side
1597         uint16_t _charOff;     // character offset within side
1598         bool _fw;              // side is forward or backward?
1599         int16_t _by;           // byte within side (not adjusted for bw sides)
1600         int8_t _bp;            // bitpair within byte (not adjusted for bw sides)
1601 };
1602
1603 #include "ebwt_search_backtrack.h"
1604
1605 ///////////////////////////////////////////////////////////////////////
1606 //
1607 // Functions for printing and sanity-checking Ebwts
1608 //
1609 ///////////////////////////////////////////////////////////////////////
1610
1611 /**
1612  * Given a range of positions in the EBWT array within the BWT portion
1613  * of a forward side, print the characters at those positions along
1614  * with a summary occ[] array.
1615  */
1616 template<typename TStr>
1617 void Ebwt<TStr>::printRangeFw(uint32_t begin, uint32_t end) const {
1618         assert(isInMemory());
1619         uint32_t occ[] = {0, 0, 0, 0};
1620         assert_gt(end, begin);
1621         for(uint32_t i = begin; i < end; i++) {
1622                 uint8_t by = this->_ebwt[i];
1623                 for(int j = 0; j < 4; j++) {
1624                         // Unpack from lowest to highest bit pair
1625                         int twoBit = unpack_2b_from_8b(by, j);
1626                         occ[twoBit]++;
1627                         cout << "ACGT"[twoBit];
1628                 }
1629                 assert_eq(0, (occ[0] + occ[1] + occ[2] + occ[3]) & 3);
1630         }
1631         cout << ":{" << occ[0] << "," << occ[1] << "," << occ[2] << "," << occ[3] << "}" << endl;
1632 }
1633
1634 /**
1635  * Given a range of positions in the EBWT array within the BWT portion
1636  * of a backward side, print the characters at those positions along
1637  * with a summary occ[] array.
1638  */
1639 template<typename TStr>
1640 void Ebwt<TStr>::printRangeBw(uint32_t begin, uint32_t end) const {
1641         assert(isInMemory());
1642         uint32_t occ[] = {0, 0, 0, 0};
1643         assert_gt(end, begin);
1644         for(uint32_t i = end-1; i >= begin; i--) {
1645                 uint8_t by = this->_ebwt[i];
1646                 for(int j = 3; j >= 0; j--) {
1647                         // Unpack from lowest to highest bit pair
1648                         int twoBit = unpack_2b_from_8b(by, j);
1649                         occ[twoBit]++;
1650                         cout << "ACGT"[twoBit];
1651                 }
1652                 assert_eq(0, (occ[0] + occ[1] + occ[2] + occ[3]) & 3);
1653                 if(i == 0) break;
1654         }
1655         cout << ":{" << occ[0] << "," << occ[1] << "," << occ[2] << "," << occ[3] << "}" << endl;
1656 }
1657
1658 /**
1659  * Check that the ebwt array is internally consistent up to (and not
1660  * including) the given side index by re-counting the chars and
1661  * comparing against the embedded occ[] arrays.
1662  */
1663 template<typename TStr>
1664 void Ebwt<TStr>::sanityCheckUpToSide(int upToSide) const {
1665         assert(isInMemory());
1666         uint32_t occ[] = {0, 0, 0, 0};
1667         uint32_t occ_save[] = {0, 0};
1668         uint32_t cur = 0; // byte pointer
1669         const EbwtParams& eh = this->_eh;
1670         bool fw = false;
1671         while(cur < (upToSide * eh._sideSz)) {
1672                 assert_leq(cur + eh._sideSz, eh._ebwtTotLen);
1673                 for(uint32_t i = 0; i < eh._sideBwtSz; i++) {
1674                         uint8_t by = this->_ebwt[cur + (fw ? i : eh._sideBwtSz-i-1)];
1675                         for(int j = 0; j < 4; j++) {
1676                                 // Unpack from lowest to highest bit pair
1677                                 int twoBit = unpack_2b_from_8b(by, fw ? j : 3-j);
1678                                 occ[twoBit]++;
1679                                 //if(_verbose) cout << "ACGT"[twoBit];
1680                         }
1681                         assert_eq(0, (occ[0] + occ[1] + occ[2] + occ[3]) % 4);
1682                 }
1683                 assert_eq(0, (occ[0] + occ[1] + occ[2] + occ[3]) % eh._sideBwtLen);
1684                 if(fw) {
1685                         // Finished forward bucket; check saved [G] and [T]
1686                         // against the two uint32_ts encoded here
1687                         ASSERT_ONLY(uint32_t *u32ebwt = reinterpret_cast<uint32_t*>(&this->_ebwt[cur + eh._sideBwtSz]));
1688                         ASSERT_ONLY(uint32_t gs = u32ebwt[0]);
1689                         ASSERT_ONLY(uint32_t ts = u32ebwt[1]);
1690                         assert_eq(gs, occ_save[0]);
1691                         assert_eq(ts, occ_save[1]);
1692                         fw = false;
1693                 } else {
1694                         // Finished backward bucket; check current [A] and [C]
1695                         // against the two uint32_ts encoded here
1696                         ASSERT_ONLY(uint32_t *u32ebwt = reinterpret_cast<uint32_t*>(&this->_ebwt[cur + eh._sideBwtSz]));
1697                         ASSERT_ONLY(uint32_t as = u32ebwt[0]);
1698                         ASSERT_ONLY(uint32_t cs = u32ebwt[1]);
1699                         assert(as == occ[0] || as == occ[0]-1); // one 'a' is a skipped '$' and doesn't count toward occ[]
1700                         assert_eq(cs, occ[1]);
1701                         occ_save[0] = occ[2]; // save gs
1702                         occ_save[1] = occ[3]; // save ts
1703                         fw = true;
1704                 }
1705                 cur += eh._sideSz;
1706         }
1707 }
1708
1709 /**
1710  * Sanity-check various pieces of the Ebwt
1711  */
1712 template<typename TStr>
1713 void Ebwt<TStr>::sanityCheckAll(int reverse) const {
1714         const EbwtParams& eh = this->_eh;
1715         assert(isInMemory());
1716         // Check ftab
1717         for(uint32_t i = 1; i < eh._ftabLen; i++) {
1718                 assert_geq(this->ftabHi(i), this->ftabLo(i-1));
1719                 assert_geq(this->ftabLo(i), this->ftabHi(i-1));
1720                 assert_leq(this->ftabHi(i), eh._bwtLen+1);
1721         }
1722         assert_eq(this->ftabHi(eh._ftabLen-1), eh._bwtLen);
1723
1724         // Check offs
1725         int seenLen = (eh._bwtLen + 31) >> 5;
1726         uint32_t *seen;
1727         try {
1728                 seen = new uint32_t[seenLen]; // bitvector marking seen offsets
1729         } catch(bad_alloc& e) {
1730                 cerr << "Out of memory allocating seen[] at " << __FILE__ << ":" << __LINE__ << endl;
1731                 throw e;
1732         }
1733         memset(seen, 0, 4 * seenLen);
1734         uint32_t offsLen = eh._offsLen;
1735         for(uint32_t i = 0; i < offsLen; i++) {
1736                 assert_lt(this->_offs[i], eh._bwtLen);
1737                 int w = this->_offs[i] >> 5;
1738                 int r = this->_offs[i] & 31;
1739                 assert_eq(0, (seen[w] >> r) & 1); // shouldn't have been seen before
1740                 seen[w] |= (1 << r);
1741         }
1742         delete[] seen;
1743
1744         // Check nPat
1745         assert_gt(this->_nPat, 0);
1746
1747         // Check plen, flen
1748         for(uint32_t i = 0; i < this->_nPat; i++) {
1749                 assert_geq(this->_plen[i], 0);
1750         }
1751
1752         // Check rstarts
1753         for(uint32_t i = 0; i < this->_nFrag-1; i++) {
1754                 assert_gt(this->_rstarts[(i+1)*3], this->_rstarts[i*3]);
1755                 if(reverse == REF_READ_REVERSE) {
1756                         assert(this->_rstarts[(i*3)+1] >= this->_rstarts[((i+1)*3)+1]);
1757                 } else {
1758                         assert(this->_rstarts[(i*3)+1] <= this->_rstarts[((i+1)*3)+1]);
1759                 }
1760         }
1761
1762         // Check ebwt
1763         sanityCheckUpToSide(eh._numSides);
1764         VMSG_NL("Ebwt::sanityCheck passed");
1765 }
1766
1767 ///////////////////////////////////////////////////////////////////////
1768 //
1769 // Functions for searching Ebwts
1770 //
1771 ///////////////////////////////////////////////////////////////////////
1772
1773 /**
1774  * Return the final character in row i (i.e. the i'th character in the
1775  * BWT transform).  Note that the 'L' in the name of the function
1776  * stands for 'last', as in the literature.
1777  */
1778 template<typename TStr>
1779 inline int Ebwt<TStr>::rowL(const SideLocus& l) const {
1780         // Extract and return appropriate bit-pair
1781 #ifdef SIXTY4_FORMAT
1782         return (((uint64_t*)l.side(this->_ebwt))[l._by >> 3] >> ((((l._by & 7) << 2) + l._bp) << 1)) & 3;
1783 #else
1784         return unpack_2b_from_8b(l.side(this->_ebwt)[l._by], l._bp);
1785 #endif
1786 }
1787
1788 /**
1789  * Inline-function version of the above.  This does not always seem to
1790  * be inlined
1791  */
1792 #if 0
1793 // Use gcc's intrinsic popcountll.  I don't recommend it because it
1794 // seems to be somewhat slower than the bit-bashing pop64 routine both
1795 // on an AMD server and on an Intel workstation.  On the other hand,
1796 // perhaps when the builtin is used GCC is smart enough to insert a
1797 // pop-count instruction on architectures that have one (e.g. Itanium).
1798 // For now, it's disabled.
1799 #define pop64(x) __builtin_popcountll(x)
1800 #elif 0
1801 __declspec naked int __stdcall pop64
1802 (uint64_t v)
1803 {
1804 static const uint64_t C55 = 0x5555555555555555ll;
1805 static const uint64_t C33 = 0x3333333333333333ll;
1806 static const uint64_t C0F = 0x0F0F0F0F0F0F0F0Fll;
1807 __asm {
1808    MOVD      MM0, [ESP+4] ;v_low
1809    PUNPCKLDQ MM0, [ESP+8] ;v
1810    MOVQ      MM1, MM0     ;v
1811    PSRLD     MM0, 1       ;v >> 1
1812    PAND      MM0, [C55]   ;(v >> 1) & 0x55555555
1813    PSUBD     MM1, MM0     ;w = v - ((v >> 1) & 0x55555555)
1814    MOVQ      MM0, MM1     ;w
1815    PSRLD     MM1, 2       ;w >> 2
1816    PAND      MM0, [C33]   ;w & 0x33333333
1817    PAND      MM1, [C33]   ;(w >> 2)  & 0x33333333
1818    PADDD     MM0, MM1     ;x = (w & 0x33333333) +
1819                           ; ((w >> 2) & 0x33333333)
1820    MOVQ      MM1, MM0     ;x
1821    PSRLD     MM0, 4       ;x >> 4
1822    PADDD     MM0, MM1     ;x + (x >> 4)
1823    PAND      MM0, [C0F]   ;y = (x + (x >> 4) & 0x0F0F0F0F)
1824    PXOR      MM1, MM1     ;0
1825    PSADBW    (MM0, MM1)   ;sum across all 8 bytes
1826    MOVD      EAX, MM0     ;result in EAX per calling
1827                           ; convention
1828    EMMS                   ;clear MMX state
1829    RET  8                 ;pop 8-byte argument off stack
1830                           ; and return
1831    }
1832 }
1833 #elif 0
1834 // Use a bytewise LUT version of popcount.  This is slower than the
1835 // bit-bashing pop64 routine both on an AMD server and on an Intel
1836 // workstation.  It seems to be about the same speed as the GCC builtin
1837 // on Intel, and a bit faster than it on AMD.  For now, it's disabled.
1838 const int popcntU8Table[256] = {
1839     0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
1840     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
1841     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
1842     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
1843     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
1844     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
1845     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
1846     3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
1847 };
1848
1849 // Use this bytewise population count table
1850 inline static int pop64(uint64_t x) {
1851         const unsigned char * p = (const unsigned char *) &x;
1852         return popcntU8Table[p[0]] +
1853                popcntU8Table[p[1]] +
1854                popcntU8Table[p[2]] +
1855                popcntU8Table[p[3]] +
1856                popcntU8Table[p[4]] +
1857                popcntU8Table[p[5]] +
1858                popcntU8Table[p[6]] +
1859                popcntU8Table[p[7]];
1860 }
1861 #else
1862 // Use this standard bit-bashing population count
1863 inline static int pop64(uint64_t x) {
1864    x = x - ((x >> 1) & 0x5555555555555555llu);
1865    x = (x & 0x3333333333333333llu) + ((x >> 2) & 0x3333333333333333llu);
1866    x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0Fllu;
1867    x = x + (x >> 8);
1868    x = x + (x >> 16);
1869    x = x + (x >> 32);
1870    return x & 0x3F;
1871 }
1872 #endif
1873
1874 /**
1875  * Tricky-bit-bashing bitpair counting for given two-bit value (0-3)
1876  * within a 64-bit argument.
1877  */
1878 inline static int countInU64(int c, uint64_t dw) {
1879         uint64_t dwA  = dw &  0xAAAAAAAAAAAAAAAAllu;
1880         uint64_t dwNA = dw & ~0xAAAAAAAAAAAAAAAAllu;
1881         uint64_t tmp;
1882         switch(c) {
1883         case 0:
1884                 tmp = (dwA >> 1) | dwNA;
1885                 break;
1886         case 1:
1887                 tmp = ~(dwA >> 1) & dwNA;
1888                 break;
1889         case 2:
1890                 tmp = (dwA >> 1) & ~dwNA;
1891                 break;
1892         case 3:
1893                 tmp = (dwA >> 1) & dwNA;
1894                 break;
1895         default:
1896                 throw;
1897         }
1898         tmp = pop64(tmp); // Gets 7.62% in profile
1899         if(c == 0) {
1900                 tmp = 32 - tmp;
1901         }
1902         assert_leq(tmp, 32);
1903         assert_geq(tmp, 0);
1904         return (int)tmp;
1905 }
1906
1907 /**
1908  * Tricky-bit-bashing bitpair counting for given two-bit value (0-3)
1909  * within a 64-bit argument.
1910  *
1911  * Function gets 2.32% in profile
1912  */
1913 inline static void countInU64Ex(uint64_t dw, uint32_t* arrs) {
1914         uint64_t dwA  = dw &  0xAAAAAAAAAAAAAAAAllu;
1915         uint64_t dwNA = dw & ~0xAAAAAAAAAAAAAAAAllu;
1916         arrs[0] += (32 - pop64((dwA >> 1) | dwNA));
1917         arrs[1] += pop64(~(dwA >> 1) & dwNA);
1918         arrs[2] += pop64((dwA >> 1) & ~dwNA);
1919         arrs[3] += pop64((dwA >> 1) & dwNA);
1920 }
1921
1922 /**
1923  * Counts the number of occurrences of character 'c' in the given Ebwt
1924  * side up to (but not including) the given byte/bitpair (by/bp).
1925  *
1926  * This is a performance-critical function.  This is the top search-
1927  * related hit in the time profile.
1928  *
1929  * Function gets 11.09% in profile
1930  */
1931 template<typename TStr>
1932 inline uint32_t Ebwt<TStr>::countUpTo(const SideLocus& l, int c) const {
1933         // Count occurrences of c in each 64-bit (using bit trickery);
1934         // Someday countInU64() and pop() functions should be
1935         // vectorized/SSE-ized in case that helps.
1936         uint32_t cCnt = 0;
1937         const uint8_t *side = l.side(this->_ebwt);
1938         int i = 0;
1939 #if 1
1940         for(; i + 7 < l._by; i += 8) {
1941                 cCnt += countInU64(c, *(uint64_t*)&side[i]);
1942         }
1943 #else
1944         for(; i + 2 < l._by; i += 2) {
1945                 cCnt += cCntLUT_16b_4[c][*(uint16_t*)&side[i]];
1946         }
1947 #endif
1948 #ifdef SIXTY4_FORMAT
1949         // Calculate number of bit pairs to shift off the end
1950         const int bpShiftoff = 32 - (((l._by & 7) << 2) + l._bp);
1951         if(bpShiftoff < 32) {
1952                 assert_lt(bpShiftoff, 32);
1953                 const uint64_t sw = (*(uint64_t*)&side[i]) << (bpShiftoff << 1);
1954                 cCnt += countInU64(c, sw);
1955                 if(c == 0) cCnt -= bpShiftoff; // we turned these into As
1956         }
1957 #else
1958         // Count occurences of c in the rest of the side (using LUT)
1959         for(; i < l._by; i++) {
1960                 cCnt += cCntLUT_4[0][c][side[i]];
1961         }
1962         // Count occurences of c in the rest of the byte
1963         if(l._bp > 0) {
1964                 cCnt += cCntLUT_4[(int)l._bp][c][side[i]];
1965         }
1966 #endif
1967         return cCnt;
1968 }
1969
1970 /**
1971  * Counts the number of occurrences of character 'c' in the given Ebwt
1972  * side up to (but not including) the given byte/bitpair (by/bp).
1973  */
1974 template<typename TStr>
1975 inline void Ebwt<TStr>::countUpToEx(const SideLocus& l, uint32_t* arrs) const {
1976         int i = 0;
1977         // Count occurrences of c in each 64-bit (using bit trickery);
1978         // note: this seems does not seem to lend a significant boost to
1979         // performance.  If you comment out this whole loop (which won't
1980         // affect correctness - it will just cause the following loop to
1981         // take up the slack) then runtime does not change noticeably.
1982         // Someday the countInU64() and pop() functions should be
1983         // vectorized/SSE-ized in case that helps.
1984         const uint8_t *side = l.side(this->_ebwt);
1985         for(; i+7 < l._by; i += 8) {
1986                 countInU64Ex(*(uint64_t*)&side[i], arrs);
1987         }
1988 #ifdef SIXTY4_FORMAT
1989         // Calculate number of bit pairs to shift off the end
1990         const int bpShiftoff = 32 - (((l._by & 7) << 2) + l._bp);
1991         assert_leq(bpShiftoff, 32);
1992         if(bpShiftoff < 32) {
1993                 const uint64_t sw = (*(uint64_t*)&l.side(this->_ebwt)[i]) << (bpShiftoff << 1);
1994                 countInU64Ex(sw, arrs);
1995                 arrs[0] -= bpShiftoff;
1996         }
1997 #else
1998         // Count occurences of c in the rest of the side (using LUT)
1999         for(; i < l._by; i++) {
2000                 arrs[0] += cCntLUT_4[0][0][side[i]];
2001                 arrs[1] += cCntLUT_4[0][1][side[i]];
2002                 arrs[2] += cCntLUT_4[0][2][side[i]];
2003                 arrs[3] += cCntLUT_4[0][3][side[i]];
2004         }
2005         // Count occurences of c in the rest of the byte
2006         if(l._bp > 0) {
2007                 arrs[0] += cCntLUT_4[(int)l._bp][0][side[i]];
2008                 arrs[1] += cCntLUT_4[(int)l._bp][1][side[i]];
2009                 arrs[2] += cCntLUT_4[(int)l._bp][2][side[i]];
2010                 arrs[3] += cCntLUT_4[(int)l._bp][3][side[i]];
2011         }
2012 #endif
2013 }
2014
2015 /**
2016  * Count all occurrences of character c from the beginning of the
2017  * forward side to <by,bp> and add in the occ[] count up to the side
2018  * break just prior to the side.
2019  */
2020 template<typename TStr>
2021 inline uint32_t Ebwt<TStr>::countFwSide(const SideLocus& l, int c) const {
2022         assert_lt(c, 4);
2023         assert_geq(c, 0);
2024         assert_lt(l._by, (int)this->_eh._sideBwtSz);
2025         assert_geq(l._by, 0);
2026         assert_lt(l._bp, 4);
2027         assert_geq(l._bp, 0);
2028         const uint8_t *side = l.side(this->_ebwt);
2029         uint32_t cCnt = countUpTo(l, c);
2030         assert_leq(cCnt, this->_eh._sideBwtLen);
2031         if(c == 0 && l._sideByteOff <= _zEbwtByteOff && l._sideByteOff + l._by >= _zEbwtByteOff) {
2032                 // Adjust for the fact that we represented $ with an 'A', but
2033                 // shouldn't count it as an 'A' here
2034                 if((l._sideByteOff + l._by > _zEbwtByteOff) ||
2035                    (l._sideByteOff + l._by == _zEbwtByteOff && l._bp > _zEbwtBpOff))
2036                 {
2037                         cCnt--; // Adjust for '$' looking like an 'A'
2038                 }
2039         }
2040         uint32_t ret;
2041         // Now factor in the occ[] count at the side break
2042         if(c < 2) {
2043                 const uint32_t *ac = reinterpret_cast<const uint32_t*>(side - 8);
2044                 assert_leq(ac[0], this->_eh._numSides * this->_eh._sideBwtLen); // b/c it's used as padding
2045                 assert_leq(ac[1], this->_eh._len);
2046                 ret = ac[c] + cCnt + this->_fchr[c];
2047         } else {
2048                 const uint32_t *gt = reinterpret_cast<const uint32_t*>(side + this->_eh._sideSz - 8); // next
2049                 assert_leq(gt[0], this->_eh._len); assert_leq(gt[1], this->_eh._len);
2050                 ret = gt[c-2] + cCnt + this->_fchr[c];
2051         }
2052 #ifndef NDEBUG
2053         assert_leq(ret, this->_fchr[c+1]); // can't have jumpded into next char's section
2054         if(c == 0) {
2055                 assert_leq(cCnt, this->_eh._sideBwtLen);
2056         } else {
2057                 assert_leq(ret, this->_eh._bwtLen);
2058         }
2059 #endif
2060         return ret;
2061 }
2062
2063 /**
2064  * Count all occurrences of character c from the beginning of the
2065  * forward side to <by,bp> and add in the occ[] count up to the side
2066  * break just prior to the side.
2067  */
2068 template<typename TStr>
2069 inline void Ebwt<TStr>::countFwSideEx(const SideLocus& l, uint32_t* arrs) const
2070 {
2071         assert_lt(l._by, (int)this->_eh._sideBwtSz);
2072         assert_geq(l._by, 0);
2073         assert_lt(l._bp, 4);
2074         assert_geq(l._bp, 0);
2075         countUpToEx(l, arrs);
2076 #ifndef NDEBUG
2077         assert_leq(arrs[0], this->_fchr[1]); // can't have jumped into next char's section
2078         assert_leq(arrs[1], this->_fchr[2]); // can't have jumped into next char's section
2079         assert_leq(arrs[2], this->_fchr[3]); // can't have jumped into next char's section
2080         assert_leq(arrs[3], this->_fchr[4]); // can't have jumped into next char's section
2081 #endif
2082         assert_leq(arrs[0], this->_eh._sideBwtLen);
2083         assert_leq(arrs[1], this->_eh._sideBwtLen);
2084         assert_leq(arrs[2], this->_eh._sideBwtLen);
2085         assert_leq(arrs[3], this->_eh._sideBwtLen);
2086         const uint8_t *side = l.side(this->_ebwt);
2087         if(l._sideByteOff <= _zEbwtByteOff && l._sideByteOff + l._by >= _zEbwtByteOff) {
2088                 // Adjust for the fact that we represented $ with an 'A', but
2089                 // shouldn't count it as an 'A' here
2090                 if((l._sideByteOff + l._by > _zEbwtByteOff) ||
2091                    (l._sideByteOff + l._by == _zEbwtByteOff && l._bp > _zEbwtBpOff))
2092                 {
2093                         arrs[0]--; // Adjust for '$' looking like an 'A'
2094                 }
2095         }
2096         // Now factor in the occ[] count at the side break
2097         const uint32_t *ac = reinterpret_cast<const uint32_t*>(side - 8);
2098         const uint32_t *gt = reinterpret_cast<const uint32_t*>(side + this->_eh._sideSz - 8);
2099 #ifndef NDEBUG
2100         assert_leq(ac[0], this->_fchr[1] + this->_eh.sideBwtLen());
2101         assert_leq(ac[1], this->_fchr[2]-this->_fchr[1]);
2102         assert_leq(gt[0], this->_fchr[3]-this->_fchr[2]);
2103         assert_leq(gt[1], this->_fchr[4]-this->_fchr[3]);
2104 #endif
2105         assert_leq(ac[0], this->_eh._len + this->_eh.sideBwtLen()); assert_leq(ac[1], this->_eh._len);
2106         assert_leq(gt[0], this->_eh._len); assert_leq(gt[1], this->_eh._len);
2107         arrs[0] += (ac[0] + this->_fchr[0]);
2108         arrs[1] += (ac[1] + this->_fchr[1]);
2109         arrs[2] += (gt[0] + this->_fchr[2]);
2110         arrs[3] += (gt[1] + this->_fchr[3]);
2111 #ifndef NDEBUG
2112         assert_leq(arrs[0], this->_fchr[1]); // can't have jumpded into next char's section
2113         assert_leq(arrs[1], this->_fchr[2]); // can't have jumpded into next char's section
2114         assert_leq(arrs[2], this->_fchr[3]); // can't have jumpded into next char's section
2115         assert_leq(arrs[3], this->_fchr[4]); // can't have jumpded into next char's section
2116 #endif
2117 }
2118
2119 /**
2120  * Count all instances of character c from <by,bp> to the logical end
2121  * (actual beginning) of the backward side, and subtract that from the
2122  * occ[] count up to the side break.
2123  */
2124 template<typename TStr>
2125 inline uint32_t Ebwt<TStr>::countBwSide(const SideLocus& l, int c) const {
2126         assert_lt(c, 4);
2127         assert_geq(c, 0);
2128         assert_lt(l._by, (int)this->_eh._sideBwtSz);
2129         assert_geq(l._by, 0);
2130         assert_lt(l._bp, 4);
2131         assert_geq(l._bp, 0);
2132         const uint8_t *side = l.side(this->_ebwt);
2133         uint32_t cCnt = countUpTo(l, c);
2134         if(rowL(l) == c) cCnt++;
2135         assert_leq(cCnt, this->_eh._sideBwtLen);
2136         if(c == 0 && l._sideByteOff <= _zEbwtByteOff && l._sideByteOff + l._by >= _zEbwtByteOff) {
2137                 // Adjust for the fact that we represented $ with an 'A', but
2138                 // shouldn't count it as an 'A' here
2139                 if((l._sideByteOff + l._by > _zEbwtByteOff) ||
2140                    (l._sideByteOff + l._by == _zEbwtByteOff && l._bp >= _zEbwtBpOff))
2141                 {
2142                         cCnt--;
2143                 }
2144         }
2145         uint32_t ret;
2146         // Now factor in the occ[] count at the side break
2147         if(c < 2) {
2148                 const uint32_t *ac = reinterpret_cast<const uint32_t*>(side + this->_eh._sideSz - 8);
2149                 assert_leq(ac[0], this->_eh._numSides * this->_eh._sideBwtLen); // b/c it's used as padding
2150                 assert_leq(ac[1], this->_eh._len);
2151                 ret = ac[c] - cCnt + this->_fchr[c];
2152         } else {
2153                 const uint32_t *gt = reinterpret_cast<const uint32_t*>(side + (2*this->_eh._sideSz) - 8); // next
2154                 assert_leq(gt[0], this->_eh._len); assert_leq(gt[1], this->_eh._len);
2155                 ret = gt[c-2] - cCnt + this->_fchr[c];
2156         }
2157 #ifndef NDEBUG
2158         assert_leq(ret, this->_fchr[c+1]); // can't have jumped into next char's section
2159         if(c == 0) {
2160                 assert_leq(cCnt, this->_eh._sideBwtLen);
2161         } else {
2162                 assert_lt(ret, this->_eh._bwtLen);
2163         }
2164 #endif
2165         return ret;
2166 }
2167
2168 /**
2169  * Count all instances of character c from <by,bp> to the logical end
2170  * (actual beginning) of the backward side, and subtract that from the
2171  * occ[] count up to the side break.
2172  */
2173 template<typename TStr>
2174 inline void Ebwt<TStr>::countBwSideEx(const SideLocus& l, uint32_t* arrs) const {
2175         assert_lt(l._by, (int)this->_eh._sideBwtSz);
2176         assert_geq(l._by, 0);
2177         assert_lt(l._bp, 4);
2178         assert_geq(l._bp, 0);
2179         const uint8_t *side = l.side(this->_ebwt);
2180         countUpToEx(l, arrs);
2181         arrs[rowL(l)]++;
2182         assert_leq(arrs[0], this->_eh._sideBwtLen);
2183         assert_leq(arrs[1], this->_eh._sideBwtLen);
2184         assert_leq(arrs[2], this->_eh._sideBwtLen);
2185         assert_leq(arrs[3], this->_eh._sideBwtLen);
2186         if(l._sideByteOff <= _zEbwtByteOff && l._sideByteOff + l._by >= _zEbwtByteOff) {
2187                 // Adjust for the fact that we represented $ with an 'A', but
2188                 // shouldn't count it as an 'A' here
2189                 if((l._sideByteOff + l._by > _zEbwtByteOff) ||
2190                    (l._sideByteOff + l._by == _zEbwtByteOff && l._bp >= _zEbwtBpOff))
2191                 {
2192                         arrs[0]--; // Adjust for '$' looking like an 'A'
2193                 }
2194         }
2195         // Now factor in the occ[] count at the side break
2196         const uint32_t *ac = reinterpret_cast<const uint32_t*>(side + this->_eh._sideSz - 8);
2197         const uint32_t *gt = reinterpret_cast<const uint32_t*>(side + (2*this->_eh._sideSz) - 8);
2198 #ifndef NDEBUG
2199         assert_leq(ac[0], this->_fchr[1] + this->_eh.sideBwtLen());
2200         assert_leq(ac[1], this->_fchr[2]-this->_fchr[1]);
2201         assert_leq(gt[0], this->_fchr[3]-this->_fchr[2]);
2202         assert_leq(gt[1], this->_fchr[4]-this->_fchr[3]);
2203 #endif
2204         assert_leq(ac[0], this->_eh._len + this->_eh.sideBwtLen()); assert_leq(ac[1], this->_eh._len);
2205         assert_leq(gt[0], this->_eh._len); assert_leq(gt[1], this->_eh._len);
2206         arrs[0] = (ac[0] - arrs[0] + this->_fchr[0]);
2207         arrs[1] = (ac[1] - arrs[1] + this->_fchr[1]);
2208         arrs[2] = (gt[0] - arrs[2] + this->_fchr[2]);
2209         arrs[3] = (gt[1] - arrs[3] + this->_fchr[3]);
2210 #ifndef NDEBUG
2211         assert_leq(arrs[0], this->_fchr[1]); // can't have jumped into next char's section
2212         assert_leq(arrs[1], this->_fchr[2]); // can't have jumped into next char's section
2213         assert_leq(arrs[2], this->_fchr[3]); // can't have jumped into next char's section
2214         assert_leq(arrs[3], this->_fchr[4]); // can't have jumped into next char's section
2215 #endif
2216 }
2217
2218 /**
2219  * Given top and bot loci, calculate counts of all four DNA chars up to
2220  * those loci.  Used for more advanced backtracking-search.
2221  */
2222 template<typename TStr>
2223 inline void Ebwt<TStr>::mapLFEx(const SideLocus& ltop,
2224                                 const SideLocus& lbot,
2225                                 uint32_t *tops,
2226                                 uint32_t *bots
2227                                 ASSERT_ONLY(, bool overrideSanity)
2228                                 ) const
2229 {
2230         // TODO: Where there's overlap, reuse the count for the overlapping
2231         // portion
2232 #ifdef EBWT_STATS
2233         const_cast<Ebwt<TStr>*>(this)->mapLFExs_++;
2234 #endif
2235         assert_eq(0, tops[0]); assert_eq(0, bots[0]);
2236         assert_eq(0, tops[1]); assert_eq(0, bots[1]);
2237         assert_eq(0, tops[2]); assert_eq(0, bots[2]);
2238         assert_eq(0, tops[3]); assert_eq(0, bots[3]);
2239         if(ltop._fw) countFwSideEx(ltop, tops); // Forward side
2240         else         countBwSideEx(ltop, tops); // Backward side
2241         if(lbot._fw) countFwSideEx(lbot, bots); // Forward side
2242         else         countBwSideEx(lbot, bots); // Backward side
2243 #ifndef NDEBUG
2244         if(_sanity && !overrideSanity) {
2245                 // Make sure results match up with individual calls to mapLF;
2246                 // be sure to override sanity-checking in the callee, or we'll
2247                 // have infinite recursion
2248                 assert_eq(mapLF(ltop, 0, true), tops[0]);
2249                 assert_eq(mapLF(ltop, 1, true), tops[1]);
2250                 assert_eq(mapLF(ltop, 2, true), tops[2]);
2251                 assert_eq(mapLF(ltop, 3, true), tops[3]);
2252                 assert_eq(mapLF(lbot, 0, true), bots[0]);
2253                 assert_eq(mapLF(lbot, 1, true), bots[1]);
2254                 assert_eq(mapLF(lbot, 2, true), bots[2]);
2255                 assert_eq(mapLF(lbot, 3, true), bots[3]);
2256         }
2257 #endif
2258 }
2259
2260 #ifndef NDEBUG
2261 /**
2262  * Given top and bot loci, calculate counts of all four DNA chars up to
2263  * those loci.  Used for more advanced backtracking-search.
2264  */
2265 template<typename TStr>
2266 inline void Ebwt<TStr>::mapLFEx(const SideLocus& l,
2267                                 uint32_t *arrs
2268                                 ASSERT_ONLY(, bool overrideSanity)
2269                                 ) const
2270 {
2271         assert_eq(0, arrs[0]);
2272         assert_eq(0, arrs[1]);
2273         assert_eq(0, arrs[2]);
2274         assert_eq(0, arrs[3]);
2275         if(l._fw) countFwSideEx(l, arrs); // Forward side
2276         else      countBwSideEx(l, arrs); // Backward side
2277 #ifndef NDEBUG
2278         if(_sanity && !overrideSanity) {
2279                 // Make sure results match up with individual calls to mapLF;
2280                 // be sure to override sanity-checking in the callee, or we'll
2281                 // have infinite recursion
2282                 assert_eq(mapLF(l, 0, true), arrs[0]);
2283                 assert_eq(mapLF(l, 1, true), arrs[1]);
2284                 assert_eq(mapLF(l, 2, true), arrs[2]);
2285                 assert_eq(mapLF(l, 3, true), arrs[3]);
2286         }
2287 #endif
2288 }
2289 #endif
2290
2291 /**
2292  * Given row i, return the row that the LF mapping maps i to.
2293  */
2294 template<typename TStr>
2295 inline uint32_t Ebwt<TStr>::mapLF(const SideLocus& l
2296                                   ASSERT_ONLY(, bool overrideSanity)
2297                                   ) const
2298 {
2299 #ifdef EBWT_STATS
2300         const_cast<Ebwt<TStr>*>(this)->mapLFs_++;
2301 #endif
2302         uint32_t ret;
2303         assert(l.side(this->_ebwt) != NULL);
2304         int c = rowL(l);
2305         assert_lt(c, 4);
2306         assert_geq(c, 0);
2307         if(l._fw) ret = countFwSide(l, c); // Forward side
2308         else      ret = countBwSide(l, c); // Backward side
2309         assert_lt(ret, this->_eh._bwtLen);
2310 #ifndef NDEBUG
2311         if(_sanity && !overrideSanity) {
2312                 // Make sure results match up with results from mapLFEx;
2313                 // be sure to override sanity-checking in the callee, or we'll
2314                 // have infinite recursion
2315                 uint32_t arrs[] = { 0, 0, 0, 0 };
2316                 mapLFEx(l, arrs, true);
2317                 assert_eq(arrs[c], ret);
2318         }
2319 #endif
2320         return ret;
2321 }
2322
2323 /**
2324  * Given row i and character c, return the row that the LF mapping maps
2325  * i to on character c.
2326  */
2327 template<typename TStr>
2328 inline uint32_t Ebwt<TStr>::mapLF(const SideLocus& l, int c
2329                                   ASSERT_ONLY(, bool overrideSanity)
2330                                   ) const
2331 {
2332 #ifdef EBWT_STATS
2333         const_cast<Ebwt<TStr>*>(this)->mapLFcs_++;
2334 #endif
2335         uint32_t ret;
2336         assert_lt(c, 4);
2337         assert_geq(c, 0);
2338         if(l._fw) ret = countFwSide(l, c); // Forward side
2339         else      ret = countBwSide(l, c); // Backward side
2340         assert_lt(ret, this->_eh._bwtLen);
2341 #ifndef NDEBUG
2342         if(_sanity && !overrideSanity) {
2343                 // Make sure results match up with results from mapLFEx;
2344                 // be sure to override sanity-checking in the callee, or we'll
2345                 // have infinite recursion
2346                 uint32_t arrs[] = { 0, 0, 0, 0 };
2347                 mapLFEx(l, arrs, true);
2348                 assert_eq(arrs[c], ret);
2349         }
2350 #endif
2351         return ret;
2352 }
2353
2354 /**
2355  * Given row i and character c, return the row that the LF mapping maps
2356  * i to on character c.
2357  */
2358 template<typename TStr>
2359 inline uint32_t Ebwt<TStr>::mapLF1(uint32_t row, const SideLocus& l, int c
2360                                    ASSERT_ONLY(, bool overrideSanity)
2361                                    ) const
2362 {
2363 #ifdef EBWT_STATS
2364         const_cast<Ebwt<TStr>*>(this)->mapLF1cs_++;
2365 #endif
2366         if(rowL(l) != c || row == _zOff) return 0xffffffff;
2367         uint32_t ret;
2368         assert_lt(c, 4);
2369         assert_geq(c, 0);
2370         if(l._fw) ret = countFwSide(l, c); // Forward side
2371         else      ret = countBwSide(l, c); // Backward side
2372         assert_lt(ret, this->_eh._bwtLen);
2373 #ifndef NDEBUG
2374         if(_sanity && !overrideSanity) {
2375                 // Make sure results match up with results from mapLFEx;
2376                 // be sure to override sanity-checking in the callee, or we'll
2377                 // have infinite recursion
2378                 uint32_t arrs[] = { 0, 0, 0, 0 };
2379                 mapLFEx(l, arrs, true);
2380                 assert_eq(arrs[c], ret);
2381         }
2382 #endif
2383         return ret;
2384 }
2385
2386 /**
2387  * Given row i and character c, return the row that the LF mapping maps
2388  * i to on character c.
2389  */
2390 template<typename TStr>
2391 inline int Ebwt<TStr>::mapLF1(uint32_t& row, const SideLocus& l
2392                               ASSERT_ONLY(, bool overrideSanity)
2393                               ) const
2394 {
2395 #ifdef EBWT_STATS
2396         const_cast<Ebwt<TStr>*>(this)->mapLF1s_++;
2397 #endif
2398         if(row == _zOff) return -1;
2399         int c = rowL(l);
2400         assert_lt(c, 4);
2401         assert_geq(c, 0);
2402         if(l._fw) row = countFwSide(l, c); // Forward side
2403         else      row = countBwSide(l, c); // Backward side
2404         assert_lt(row, this->_eh._bwtLen);
2405 #ifndef NDEBUG
2406         if(_sanity && !overrideSanity) {
2407                 // Make sure results match up with results from mapLFEx;
2408                 // be sure to override sanity-checking in the callee, or we'll
2409                 // have infinite recursion
2410                 uint32_t arrs[] = { 0, 0, 0, 0 };
2411                 mapLFEx(l, arrs, true);
2412                 assert_eq(arrs[c], row);
2413         }
2414 #endif
2415         return c;
2416 }
2417
2418 /**
2419  * Take an offset into the joined text and translate it into the
2420  * reference of the index it falls on, the offset into the reference,
2421  * and the length of the reference.  Use a binary search through the
2422  * sorted list of reference fragment ranges t
2423  */
2424 template<typename TStr>
2425 void Ebwt<TStr>::joinedToTextOff(uint32_t qlen, uint32_t off,
2426                                  uint32_t& tidx,
2427                                  uint32_t& textoff,
2428                                  uint32_t& tlen) const
2429 {
2430         uint32_t top = 0;
2431         uint32_t bot = _nFrag; // 1 greater than largest addressable element
2432         uint32_t elt = 0xffffffff;
2433         // Begin binary search
2434         while(true) {
2435                 ASSERT_ONLY(uint32_t oldelt = elt);
2436                 elt = top + ((bot - top) >> 1);
2437                 assert_neq(oldelt, elt); // must have made progress
2438                 uint32_t lower = _rstarts[elt*3];
2439                 uint32_t upper;
2440                 if(elt == _nFrag-1) {
2441                         upper = _eh._len;
2442                 } else {
2443                         upper = _rstarts[((elt+1)*3)];
2444                 }
2445                 assert_gt(upper, lower);
2446                 uint32_t fraglen = upper - lower;
2447                 if(lower <= off) {
2448                         if(upper > off) { // not last element, but it's within
2449                                 // off is in this range; check if it falls off
2450                                 if(off + qlen > upper) {
2451                                         // it falls off; signal no-go and return
2452                                         tidx = 0xffffffff;
2453                                         assert_lt(elt, _nFrag-1);
2454                                         return;
2455                                 }
2456                                 tidx = _rstarts[(elt*3)+1];
2457                                 assert_lt(tidx, this->_nPat);
2458                                 assert_leq(fraglen, this->_plen[tidx]);
2459                                 // it doesn't fall off; now calculate textoff.
2460                                 // Initially it's the number of characters that precede
2461                                 // the alignment in the fragment
2462                                 uint32_t fragoff = off - _rstarts[(elt*3)];
2463                                 if(!this->_fw) {
2464                                         fragoff = fraglen - fragoff - 1;
2465                                         fragoff -= (qlen-1);
2466                                 }
2467                                 // Add the alignment's offset into the fragment
2468                                 // ('fragoff') to the fragment's offset within the text
2469                                 textoff = fragoff + _rstarts[(elt*3)+2];
2470                                 assert_lt(textoff, this->_plen[tidx]);
2471                                 break; // done with binary search
2472                         } else {
2473                                 // 'off' belongs somewhere in the region between elt
2474                                 // and bot
2475                                 top = elt;
2476                         }
2477                 } else {
2478                         // 'off' belongs somewhere in the region between top and
2479                         // elt
2480                         bot = elt;
2481                 }
2482                 // continue with binary search
2483         }
2484         tlen = this->_plen[tidx];
2485 }
2486
2487 /**
2488  * Report a potential match at offset 'off' with pattern length
2489  * 'qlen'.  Filter out spurious matches that span texts.
2490  */
2491 template<typename TStr>
2492 inline bool Ebwt<TStr>::report(const String<Dna5>& query,
2493                                String<char>* quals,
2494                                String<char>* name,
2495                                bool color,
2496                                bool colExEnds,
2497                                int snpPhred,
2498                                const BitPairReference* ref,
2499                                const std::vector<uint32_t>& mmui32,
2500                                const std::vector<uint8_t>& refcs,
2501                                size_t numMms,
2502                                uint32_t off,
2503                                uint32_t top,
2504                                uint32_t bot,
2505                                uint32_t qlen,
2506                                int stratum,
2507                                uint16_t cost,
2508                                uint32_t patid,
2509                                uint32_t seed,
2510                                const EbwtSearchParams<TStr>& params) const
2511 {
2512         VMSG_NL("In report");
2513         assert_geq(cost, (uint32_t)(stratum << 14));
2514         assert_lt(off, this->_eh._len);
2515         uint32_t tidx;
2516         uint32_t textoff;
2517         uint32_t tlen;
2518         joinedToTextOff(qlen, off, tidx, textoff, tlen);
2519         if(tidx == 0xffffffff) {
2520                 return false;
2521         }
2522         return params.reportHit(
2523                         query,                    // read sequence
2524                         quals,                    // read quality values
2525                         name,                     // read name
2526                         color,                    // true -> read is colorspace
2527                         colExEnds,                // true -> exclude nucleotides on ends
2528                         snpPhred,                 // phred probability of SNP
2529                         ref,                      // reference sequence
2530                         rmap_,                    // map to another reference coordinate system
2531                         _fw,                      // true = index is forward; false = mirror
2532                         mmui32,                   // mismatch positions
2533                         refcs,                    // reference characters for mms
2534                         numMms,                   // # mismatches
2535                         make_pair(tidx, textoff), // position
2536                         make_pair(0, 0),          // (bogus) mate position
2537                         true,                     // (bogus) mate orientation
2538                         0,                        // (bogus) mate length
2539                         make_pair(top, bot),      // arrows
2540                         tlen,                     // textlen
2541                         qlen,                     // qlen
2542                         stratum,                  // alignment stratum
2543                         cost,                     // cost, including stratum & quality penalty
2544                         bot-top-1,                // # other hits
2545                         patid,                    // pattern id
2546                         seed,                     // pseudo-random seed
2547                         0);                       // mate (0 = unpaired)
2548 }
2549
2550 #include "row_chaser.h"
2551
2552 /**
2553  * Report a result.  Involves walking backwards along the original
2554  * string by way of the LF-mapping until we reach a marked SA row or
2555  * the row corresponding to the 0th suffix.  A marked row's offset
2556  * into the original string can be read directly from the this->_offs[]
2557  * array.
2558  */
2559 template<typename TStr>
2560 inline bool Ebwt<TStr>::reportChaseOne(const String<Dna5>& query,
2561                                        String<char>* quals,
2562                                        String<char>* name,
2563                                        bool color,
2564                                        bool colExEnds,
2565                                        int snpPhred,
2566                                        const BitPairReference* ref,
2567                                        const std::vector<uint32_t>& mmui32,
2568                                        const std::vector<uint8_t>& refcs,
2569                                        size_t numMms,
2570                                        uint32_t i,
2571                                        uint32_t top,
2572                                        uint32_t bot,
2573                                        uint32_t qlen,
2574                                        int stratum,
2575                                        uint16_t cost,
2576                                        uint32_t patid,
2577                                        uint32_t seed,
2578                                        const EbwtSearchParams<TStr>& params,
2579                                        SideLocus *l) const
2580 {
2581         VMSG_NL("In reportChaseOne");
2582         uint32_t off;
2583         uint32_t jumps = 0;
2584         ASSERT_ONLY(uint32_t origi = i);
2585         SideLocus myl;
2586         const uint32_t offMask = this->_eh._offMask;
2587         const uint32_t offRate = this->_eh._offRate;
2588         const uint32_t* offs = this->_offs;
2589         // If the caller didn't give us a pre-calculated (and prefetched)
2590         // locus, then we have to do that now
2591         if(l == NULL) {
2592                 l = &myl;
2593                 l->initFromRow(i, this->_eh, this->_ebwt);
2594         }
2595         assert(l != NULL);
2596         assert(l->valid());
2597         // Walk along until we reach the next marked row to the left
2598         while(((i & offMask) != i) && i != _zOff) {
2599                 // Not a marked row; walk left one more char
2600                 uint32_t newi = mapLF(*l); // calc next row
2601                 assert_neq(newi, i);
2602                 i = newi;                                  // update row
2603                 l->initFromRow(i, this->_eh, this->_ebwt); // update locus
2604                 jumps++;
2605         }
2606         // This is a marked row
2607         if(i == _zOff) {
2608                 // Special case: it's the row corresponding to the
2609                 // lexicographically smallest suffix, which is implicitly
2610                 // marked 0
2611                 off = jumps;
2612                 VMSG_NL("reportChaseOne found zoff off=" << off << " (jumps=" << jumps << ")");
2613         } else {
2614                 // Normal marked row, calculate offset of row i
2615                 off = offs[i >> offRate] + jumps;
2616                 VMSG_NL("reportChaseOne found off=" << off << " (jumps=" << jumps << ")");
2617         }
2618 #ifndef NDEBUG
2619         {
2620                 uint32_t rcoff = RowChaser<TStr>::toFlatRefOff(this, qlen, origi);
2621                 assert_eq(rcoff, off);
2622         }
2623 #endif
2624         return report(query, quals, name, color, colExEnds, snpPhred, ref,
2625                       mmui32, refcs, numMms, off, top, bot, qlen, stratum,
2626                       cost, patid, seed, params);
2627 }
2628
2629 /**
2630  * Report a result.  Involves walking backwards along the original
2631  * string by way of the LF-mapping until we reach a marked SA row or
2632  * the row corresponding to the 0th suffix.  A marked row's offset
2633  * into the original string can be read directly from the this->_offs[]
2634  * array.
2635  */
2636 template<typename TStr>
2637 inline bool Ebwt<TStr>::reportReconstruct(const String<Dna5>& query,
2638                                           String<char>* quals,
2639                                           String<char>* name,
2640                                           String<Dna5>& lbuf,
2641                                           String<Dna5>& rbuf,
2642                                           const uint32_t *mmui32,
2643                                           const char* refcs,
2644                                           size_t numMms,
2645                                           uint32_t i,
2646                                           uint32_t top,
2647                                           uint32_t bot,
2648                                           uint32_t qlen,
2649                                           int stratum,
2650                                           const EbwtSearchParams<TStr>& params,
2651                                           SideLocus *l) const
2652 {
2653         VMSG_NL("In reportReconstruct");
2654         assert_gt(_eh._isaLen, 0); // Must have inverse suffix array to reconstruct
2655         uint32_t off;
2656         uint32_t jumps = 0;
2657         SideLocus myl;
2658         const uint32_t offMask = this->_eh._offMask;
2659         const uint32_t offRate = this->_eh._offRate;
2660         const uint32_t* offs = this->_offs;
2661         const uint32_t* isa = this->_isa;
2662         assert(isa != NULL);
2663         if(l == NULL) {
2664                 l = &myl;
2665                 myl.initFromRow(i, this->_eh, this->_ebwt);
2666         }
2667         assert(l != NULL);
2668         clear(lbuf); clear(rbuf);
2669         // Walk along until we reach the next marked row to the left
2670         while(((i & offMask) != i) && i != _zOff) {
2671                 // Not a marked row; walk left one more char
2672                 int c = rowL(*l);
2673                 appendValue(lbuf, (Dna5)c);
2674                 uint32_t newi;
2675                 assert_lt(c, 4);
2676                 assert_geq(c, 0);
2677                 if(l->_fw) newi = countFwSide(*l, c); // Forward side
2678                 else       newi = countBwSide(*l, c); // Backward side
2679                 assert_lt(newi, this->_eh._bwtLen);
2680                 assert_neq(newi, i);
2681                 i = newi;                                  // update row
2682                 l->initFromRow(i, this->_eh, this->_ebwt); // update locus
2683                 jumps++;
2684         }
2685         // This is a marked row
2686         if(i == _zOff) {
2687                 // Special case: it's the row corresponding to the
2688                 // lexicographically smallest suffix, which is implicitly
2689                 // marked 0
2690                 off = jumps;
2691                 VMSG_NL("reportChaseOne found zoff off=" << off << " (jumps=" << jumps << ")");
2692         } else {
2693                 // Normal marked row, calculate offset of row i
2694                 off = offs[i >> offRate] + jumps;
2695                 VMSG_NL("reportChaseOne found off=" << off << " (jumps=" << jumps << ")");
2696         }
2697         // 'off' now holds the text offset of the first (leftmost) position
2698         // involved in the alignment.  Next we call joinedToTextOff to
2699         // check whether the seed is valid (i.e., does not straddle a
2700         // boundary between two reference seuqences) and to obtain its
2701         // extents
2702         uint32_t tidx;    // the index (id) of the reference we hit in
2703         uint32_t textoff; // the offset of the alignment within the reference
2704         uint32_t tlen;    // length of reference seed hit in
2705         joinedToTextOff(qlen, off, tidx, textoff, tlen);
2706         if(tidx == 0xffffffff) {
2707                 // The seed straddled a reference boundary, and so is spurious.
2708                 // Return false, indicating that we shouldn't stop.
2709                 return false;
2710         }
2711         if(jumps > textoff) {
2712                 // In our progress toward a marked row, we passed the boundary
2713                 // between the reference sequence containing the seed and the
2714                 // reference sequence to the left of it.  That's OK, we just
2715                 // need to knock off the extra characters we added to 'lbuf'.
2716                 assert_eq(jumps, length(lbuf));
2717                 _setLength(lbuf, textoff);
2718                 jumps = textoff;
2719                 assert_eq(textoff, length(lbuf));
2720         } else if(jumps < textoff) {
2721                 // Keep walking until we reach the end of the reference
2722                 assert_neq(i, _zOff);
2723                 uint32_t diff = textoff-jumps;
2724                 for(size_t j = 0; j < diff; j++) {
2725                         // Not a marked row; walk left one more char
2726                         int c = rowL(*l);
2727                         appendValue(lbuf, (Dna5)c);
2728                         uint32_t newi;
2729                         assert_lt(c, 4);
2730                         assert_geq(c, 0);
2731                         if(l->_fw) newi = countFwSide(*l, c); // Forward side
2732                         else       newi = countBwSide(*l, c); // Backward side
2733                         assert_lt(newi, this->_eh._bwtLen);
2734                         assert_neq(newi, i);
2735                         i = newi;                                  // update row
2736                         assert_neq(i, _zOff);
2737                         l->initFromRow(i, this->_eh, this->_ebwt); // update locus
2738                         jumps++;
2739                 }
2740                 assert_eq(textoff, jumps);
2741                 assert_eq(textoff, length(lbuf));
2742         }
2743         assert_eq(textoff, jumps);
2744         assert_eq(textoff, length(lbuf));
2745         // Calculate the right-hand extent of the reference
2746         uint32_t ref_right = off - textoff + tlen;
2747         // Round the right-hand extent to the nearest ISA element that maps
2748         // to it or a character to its right
2749         uint32_t ref_right_rounded = ref_right;
2750         if((ref_right_rounded & _eh._isaMask) != ref_right_rounded) {
2751                 ref_right_rounded = ((ref_right_rounded >> _eh._isaRate)+1) << _eh._isaRate;
2752         }
2753         // TODO: handle case where ref_right_rounded is off the end of _isa
2754         // Let the current suffix-array elt be determined by the ISA
2755         if((ref_right_rounded >> _eh._isaRate) >= _eh._isaLen) {
2756                 i = _eh._len;
2757                 ref_right_rounded = _eh._len;
2758         } else {
2759                 i = isa[ref_right_rounded >> _eh._isaRate];
2760         }
2761         uint32_t right_steps_rounded = ref_right_rounded - (off + qlen);
2762         uint32_t right_steps = ref_right - (off + qlen);
2763         l->initFromRow(i, this->_eh, this->_ebwt); // update locus
2764         for(size_t j = 0; j < right_steps_rounded; j++) {
2765                 // Not a marked row; walk left one more char
2766                 int c = rowL(*l);
2767                 appendValue(rbuf, (Dna5)c);
2768                 uint32_t newi;
2769                 assert_lt(c, 4); assert_geq(c, 0);
2770                 if(l->_fw) newi = countFwSide(*l, c); // Forward side
2771                 else       newi = countBwSide(*l, c); // Backward side
2772                 assert_lt(newi, this->_eh._bwtLen);
2773                 assert_neq(newi, i);
2774                 i = newi;                                  // update row
2775                 assert_neq(i, _zOff);
2776                 l->initFromRow(i, this->_eh, this->_ebwt); // update locus
2777                 jumps++;
2778         }
2779         if(right_steps_rounded > right_steps) {
2780                 jumps -= (right_steps_rounded - right_steps);
2781                 _setLength(rbuf, right_steps);
2782         }
2783         assert_eq(right_steps, length(rbuf));
2784         assert_eq(tlen, jumps + qlen);
2785         ::reverseInPlace(lbuf);
2786         ::reverseInPlace(rbuf);
2787         {
2788                 cout << "reportReconstruct:" << endl
2789                          << "  " << lbuf << query << rbuf << endl;
2790                 cout << "  ";
2791                 for(size_t i = 0; i < length(lbuf); i++) cout << " ";
2792                 cout << query << endl;
2793         }
2794         // Now we've reconstructed the
2795         return false;
2796 }
2797
2798 /**
2799  * Transform this Ebwt into the original string in linear time by using
2800  * the LF mapping to walk backwards starting at the row correpsonding
2801  * to the end of the string.  The result is written to s.  The Ebwt
2802  * must be in memory.
2803  */
2804 template<typename TStr>
2805 void Ebwt<TStr>::restore(TStr& s) const {
2806         assert(isInMemory());
2807         resize(s, this->_eh._len, Exact());
2808         uint32_t jumps = 0;
2809         uint32_t i = this->_eh._len; // should point to final SA elt (starting with '$')
2810         SideLocus l(i, this->_eh, this->_ebwt);
2811         while(i != _zOff) {
2812                 assert_lt(jumps, this->_eh._len);
2813                 //if(_verbose) cout << "restore: i: " << i << endl;
2814                 // Not a marked row; go back a char in the original string
2815                 uint32_t newi = mapLF(l);
2816                 assert_neq(newi, i);
2817                 s[this->_eh._len - jumps - 1] = rowL(l);
2818                 i = newi;
2819                 l.initFromRow(i, this->_eh, this->_ebwt);
2820                 jumps++;
2821         }
2822         assert_eq(jumps, this->_eh._len);
2823 }
2824
2825 /**
2826  * Check that this Ebwt, when restored via restore(), matches up with
2827  * the given array of reference sequences.  For sanity checking.
2828  */
2829 template <typename TStr>
2830 void Ebwt<TStr>::checkOrigs(const vector<String<Dna5> >& os,
2831                             bool color, bool mirror) const
2832 {
2833         TStr rest;
2834         restore(rest);
2835         uint32_t restOff = 0;
2836         size_t i = 0, j = 0;
2837         if(mirror) {
2838                 // TODO: FIXME
2839                 return;
2840         }
2841         while(i < os.size()) {
2842                 size_t olen = length(os[i]);
2843                 int lastorig = -1;
2844                 for(; j < olen; j++) {
2845                         size_t joff = j;
2846                         if(mirror) joff = olen - j - 1;
2847                         if((int)os[i][joff] == 4) {
2848                                 // Skip over Ns
2849                                 lastorig = -1;
2850                                 if(!mirror) {
2851                                         while(j < olen && (int)os[i][j] == 4) j++;
2852                                 } else {
2853                                         while(j < olen && (int)os[i][olen-j-1] == 4) j++;
2854                                 }
2855                                 j--;
2856                                 continue;
2857                         }
2858                         if(lastorig == -1 && color) {
2859                                 lastorig = os[i][joff];
2860                                 continue;
2861                         }
2862                         if(color) {
2863                                 assert_neq(-1, lastorig);
2864                                 assert_eq(dinuc2color[(int)os[i][joff]][lastorig], rest[restOff]);
2865                         } else {
2866                                 assert_eq(os[i][joff], rest[restOff]);
2867                         }
2868                         lastorig = (int)os[i][joff];
2869                         restOff++;
2870                 }
2871                 if(j == length(os[i])) {
2872                         // Moved to next sequence
2873                         i++;
2874                         j = 0;
2875                 } else {
2876                         // Just jumped over a gap
2877                 }
2878         }
2879 }
2880
2881 ///////////////////////////////////////////////////////////////////////
2882 //
2883 // Functions for reading and writing Ebwts
2884 //
2885 ///////////////////////////////////////////////////////////////////////
2886
2887 /**
2888  * Read an Ebwt from file with given filename.
2889  */
2890 template<typename TStr>
2891 void Ebwt<TStr>::readIntoMemory(
2892         int color,
2893         int needEntireRev,
2894         bool justHeader,
2895         EbwtParams *params,
2896         bool mmSweep,
2897         bool loadNames,
2898         bool startVerbose)
2899 {
2900         bool switchEndian; // dummy; caller doesn't care
2901 #ifdef BOWTIE_MM
2902         char *mmFile[] = { NULL, NULL };
2903 #endif
2904         if(_in1Str.length() > 0) {
2905                 if(_verbose || startVerbose) {
2906                         cerr << "  About to open input files: ";
2907                         logTime(cerr);
2908                 }
2909 #ifdef BOWTIE_MM
2910                 // Initialize our primary and secondary input-stream fields
2911                 if(_in1 != -1) close(_in1);
2912                 if(_verbose || startVerbose) {
2913                         cerr << "Opening \"" << _in1Str << "\"" << endl;
2914                 }
2915                 if((_in1 = open(_in1Str.c_str(), O_RDONLY)) < 0) {
2916                         cerr << "Could not open index file " << _in1Str << endl;
2917                 }
2918                 if(_in2 != -1) close(_in2);
2919                 if(_verbose || startVerbose) {
2920                         cerr << "Opening \"" << _in2Str << "\"" << endl;
2921                 }
2922                 if((_in2 = open(_in2Str.c_str(), O_RDONLY)) < 0) {
2923                         cerr << "Could not open index file " << _in2Str << endl;
2924                 }
2925 #else
2926                 // Initialize our primary and secondary input-stream fields
2927                 if(_in1 != NULL) fclose(_in1);
2928                 if(_verbose || startVerbose) cerr << "Opening \"" << _in1Str << "\"" << endl;
2929                 if((_in1 = fopen(_in1Str.c_str(), "rb")) == NULL) {
2930                         cerr << "Could not open index file " << _in1Str << endl;
2931                 }
2932                 if(_in2 != NULL) fclose(_in2);
2933                 if(_verbose || startVerbose) cerr << "Opening \"" << _in2Str << "\"" << endl;
2934                 if((_in2 = fopen(_in2Str.c_str(), "rb")) == NULL) {
2935                         cerr << "Could not open index file " << _in2Str << endl;
2936                 }
2937 #endif
2938                 if(_verbose || startVerbose) {
2939                         cerr << "  Finished opening input files: ";
2940                         logTime(cerr);
2941                 }
2942
2943 #ifdef BOWTIE_MM
2944                 if(_useMm && !justHeader) {
2945                         const char *names[] = {_in1Str.c_str(), _in2Str.c_str()};
2946                         int fds[] = { _in1, _in2 };
2947                         for(int i = 0; i < 2; i++) {
2948                                 if(_verbose || startVerbose) {
2949                                         cerr << "  Memory-mapping input file " << (i+1) << ": ";
2950                                         logTime(cerr);
2951                                 }
2952                                 struct stat sbuf;
2953                                 if (stat(names[i], &sbuf) == -1) {
2954                                         perror("stat");
2955                                         cerr << "Error: Could not stat index file " << names[i] << " prior to memory-mapping" << endl;
2956                                         throw 1;
2957                                 }
2958                                 mmFile[i] = (char*)mmap((void *)0, sbuf.st_size,
2959                                                                                 PROT_READ, MAP_SHARED, fds[i], 0);
2960                                 if(mmFile == (void *)(-1)) {
2961                                         perror("mmap");
2962                                         cerr << "Error: Could not memory-map the index file " << names[i] << endl;
2963                                         throw 1;
2964                                 }
2965                                 if(mmSweep) {
2966                                         int sum = 0;
2967                                         for(off_t j = 0; j < sbuf.st_size; j += 1024) {
2968                                                 sum += (int) mmFile[i][j];
2969                                         }
2970                                         if(startVerbose) {
2971                                                 cerr << "  Swept the memory-mapped ebwt index file 1; checksum: " << sum << ": ";
2972                                                 logTime(cerr);
2973                                         }
2974                                 }
2975                         }
2976                         mmFile1_ = mmFile[0];
2977                         mmFile2_ = mmFile[1];
2978                 }
2979 #endif
2980         }
2981 #ifdef BOWTIE_MM
2982         else if(_useMm && !justHeader) {
2983                 mmFile[0] = mmFile1_;
2984                 mmFile[1] = mmFile2_;
2985         }
2986         if(_useMm && !justHeader) {
2987                 assert(mmFile[0] == mmFile1_);
2988                 assert(mmFile[1] == mmFile2_);
2989         }
2990 #endif
2991
2992         if(_verbose || startVerbose) {
2993                 cerr << "  Reading header: ";
2994                 logTime(cerr);
2995         }
2996
2997         // Read endianness hints from both streams
2998         size_t bytesRead = 0;
2999         switchEndian = false;
3000         uint32_t one = readU32(_in1, switchEndian); // 1st word of primary stream
3001         bytesRead += 4;
3002         #ifndef NDEBUG
3003         assert_eq(one, readU32(_in2, switchEndian)); // should match!
3004         #else
3005         readU32(_in2, switchEndian);
3006         #endif
3007         if(one != 1) {
3008                 assert_eq((1u<<24), one);
3009                 assert_eq(1, endianSwapU32(one));
3010                 switchEndian = true;
3011         }
3012
3013         // Can't switch endianness and use memory-mapped files; in order to
3014         // support this, someone has to modify the file to switch
3015         // endiannesses appropriately, and we can't do this inside Bowtie
3016         // or we might be setting up a race condition with other processes.
3017         if(switchEndian && _useMm) {
3018                 cerr << "Error: Can't use memory-mapped files when the index is the opposite endianness" << endl;
3019                 throw 1;
3020         }
3021
3022         // Reads header entries one by one from primary stream
3023         uint32_t len          = readU32(_in1, switchEndian);
3024         bytesRead += 4;
3025         int32_t  lineRate     = readI32(_in1, switchEndian);
3026         bytesRead += 4;
3027         int32_t  linesPerSide = readI32(_in1, switchEndian);
3028         bytesRead += 4;
3029         int32_t  offRate      = readI32(_in1, switchEndian);
3030         bytesRead += 4;
3031         // TODO: add isaRate to the actual file format (right now, the
3032         // user has to tell us whether there's an ISA sample and what the
3033         // sampling rate is.
3034         int32_t  isaRate      = _overrideIsaRate;
3035         int32_t  ftabChars    = readI32(_in1, switchEndian);
3036         bytesRead += 4;
3037         // chunkRate was deprecated in an earlier version of Bowtie; now
3038         // we use it to hold flags.
3039         int32_t flags = readI32(_in1, switchEndian);
3040         bool entireRev = false;
3041         if(flags < 0 && (((-flags) & EBWT_COLOR) != 0)) {
3042                 if(color != -1 && !color) {
3043                         cerr << "Error: -C was not specified when running bowtie, but index is in colorspace.  If" << endl
3044                              << "your reads are in colorspace, please use the -C option.  If your reads are not" << endl
3045                              << "in colorspace, please use a normal index (one built without specifying -C to" << endl
3046                              << "bowtie-build)." << endl;
3047                         throw 1;
3048                 }
3049                 color = 1;
3050         } else if(flags < 0) {
3051                 if(color != -1 && color) {
3052                         cerr << "Error: -C was specified when running bowtie, but index is not in colorspace.  If" << endl
3053                              << "your reads are in colorspace, please use a colorspace index (one built using" << endl
3054                              << "bowtie-build -C).  If your reads are not in colorspace, don't specify -C when" << endl
3055                              << "running bowtie." << endl;
3056                         throw 1;
3057                 }
3058                 color = 0;
3059         }
3060         if(flags < 0 && (((-flags) & EBWT_ENTIRE_REV) == 0)) {
3061                 if(needEntireRev != -1 && needEntireRev != 0) {
3062                         cerr << "Error: This index is not compatible with this version of bowtie.  Please use a" << endl
3063                              << "current version of bowtie-build." << endl;
3064                         throw 1;
3065                 }
3066         } else entireRev = true;
3067         bytesRead += 4;
3068
3069         // Create a new EbwtParams from the entries read from primary stream
3070         EbwtParams *eh;
3071         bool deleteEh = false;
3072         if(params != NULL) {
3073                 params->init(len, lineRate, linesPerSide, offRate, isaRate, ftabChars, color, entireRev);
3074                 if(_verbose || startVerbose) params->print(cerr);
3075                 eh = params;
3076         } else {
3077                 eh = new EbwtParams(len, lineRate, linesPerSide, offRate, isaRate, ftabChars, color, entireRev);
3078                 deleteEh = true;
3079         }
3080
3081         // Set up overridden suffix-array-sample parameters
3082         uint32_t offsLen = eh->_offsLen;
3083         uint32_t offRateDiff = 0;
3084         uint32_t offsLenSampled = offsLen;
3085         if(_overrideOffRate > offRate) {
3086                 offRateDiff = _overrideOffRate - offRate;
3087         }
3088         if(offRateDiff > 0) {
3089                 offsLenSampled >>= offRateDiff;
3090                 if((offsLen & ~(0xffffffff << offRateDiff)) != 0) {
3091                         offsLenSampled++;
3092                 }
3093         }
3094
3095         // Set up overridden inverted-suffix-array-sample parameters
3096         uint32_t isaLen = eh->_isaLen;
3097         uint32_t isaRateDiff = 0;
3098         uint32_t isaLenSampled = isaLen;
3099         if(_overrideIsaRate > isaRate) {
3100                 isaRateDiff = _overrideIsaRate - isaRate;
3101         }
3102         if(isaRateDiff > 0) {
3103                 isaLenSampled >>= isaRateDiff;
3104                 if((isaLen & ~(0xffffffff << isaRateDiff)) != 0) {
3105                         isaLenSampled++;
3106                 }
3107         }
3108
3109         // Can't override the offrate or isarate and use memory-mapped
3110         // files; ultimately, all processes need to copy the sparser sample
3111         // into their own memory spaces.
3112         if(_useMm && (offRateDiff || isaRateDiff)) {
3113                 cerr << "Error: Can't use memory-mapped files when the offrate or isarate is overridden" << endl;
3114                 throw 1;
3115         }
3116
3117         // Read nPat from primary stream
3118         this->_nPat = readI32(_in1, switchEndian);
3119         bytesRead += 4;
3120         if(this->_plen != NULL && !_useMm) {
3121                 // Delete it so that we can re-read it
3122                 delete[] this->_plen;
3123                 this->_plen = NULL;
3124         }
3125
3126         // Read plen from primary stream
3127         if(_useMm) {
3128 #ifdef BOWTIE_MM
3129                 this->_plen = (uint32_t*)(mmFile[0] + bytesRead);
3130                 bytesRead += this->_nPat*4;
3131                 lseek(_in1, this->_nPat*4, SEEK_CUR);
3132 #endif
3133         } else {
3134                 try {
3135                         if(_verbose || startVerbose) {
3136                                 cerr << "Reading plen (" << this->_nPat << "): ";
3137                                 logTime(cerr);
3138                         }
3139                         this->_plen = new uint32_t[this->_nPat];
3140                         if(switchEndian) {
3141                                 for(uint32_t i = 0; i < this->_nPat; i++) {
3142                                         this->_plen[i] = readU32(_in1, switchEndian);
3143                                 }
3144                         } else {
3145                                 MM_READ_RET r = MM_READ(_in1, (void*)this->_plen, this->_nPat*4);
3146                                 if(r != (MM_READ_RET)(this->_nPat*4)) {
3147                                         cerr << "Error reading _plen[] array: " << r << ", " << (this->_nPat*4) << endl;
3148                                         throw 1;
3149                                 }
3150                         }
3151                 } catch(bad_alloc& e) {
3152                         cerr << "Out of memory allocating plen[] in Ebwt::read()"
3153                                  << " at " << __FILE__ << ":" << __LINE__ << endl;
3154                         throw e;
3155                 }
3156         }
3157
3158         bool shmemLeader;
3159
3160         // TODO: I'm not consistent on what "header" means.  Here I'm using
3161         // "header" to mean everything that would exist in memory if we
3162         // started to build the Ebwt but stopped short of the build*() step
3163         // (i.e. everything up to and including join()).
3164         if(justHeader) goto done;
3165
3166         this->_nFrag = readU32(_in1, switchEndian);
3167         bytesRead += 4;
3168         if(_verbose || startVerbose) {
3169                 cerr << "Reading rstarts (" << this->_nFrag*3 << "): ";
3170                 logTime(cerr);
3171         }
3172         assert_geq(this->_nFrag, this->_nPat);
3173         if(_useMm) {
3174 #ifdef BOWTIE_MM
3175                 this->_rstarts = (uint32_t*)(mmFile[0] + bytesRead);
3176                 bytesRead += this->_nFrag*4*3;
3177                 lseek(_in1, this->_nFrag*4*3, SEEK_CUR);
3178 #endif
3179         } else {
3180                 this->_rstarts = new uint32_t[this->_nFrag*3];
3181                 if(switchEndian) {
3182                         for(uint32_t i = 0; i < this->_nFrag*3; i += 3) {
3183                                 // fragment starting position in joined reference
3184                                 // string, text id, and fragment offset within text
3185                                 this->_rstarts[i]   = readU32(_in1, switchEndian);
3186                                 this->_rstarts[i+1] = readU32(_in1, switchEndian);
3187                                 this->_rstarts[i+2] = readU32(_in1, switchEndian);
3188                         }
3189                 } else {
3190                         MM_READ_RET r = MM_READ(_in1, (void *)this->_rstarts, this->_nFrag*4*3);
3191                         if(r != (MM_READ_RET)(this->_nFrag*4*3)) {
3192                                 cerr << "Error reading _rstarts[] array: " << r << ", " << (this->_nFrag*4*3) << endl;
3193                                 throw 1;
3194                         }
3195                 }
3196         }
3197
3198         if(_useMm) {
3199 #ifdef BOWTIE_MM
3200                 this->_ebwt = (uint8_t*)(mmFile[0] + bytesRead);
3201                 bytesRead += eh->_ebwtTotLen;
3202                 lseek(_in1, eh->_ebwtTotLen, SEEK_CUR);
3203 #endif
3204         } else {
3205                 // Allocate ebwt (big allocation)
3206                 if(_verbose || startVerbose) {
3207                         cerr << "Reading ebwt (" << eh->_ebwtTotLen << "): ";
3208                         logTime(cerr);
3209                 }
3210                 bool shmemLeader = true;
3211                 if(useShmem_) {
3212                         shmemLeader = ALLOC_SHARED_U8(
3213                                 (_in1Str + "[ebwt]"), eh->_ebwtTotLen, &this->_ebwt,
3214                                 "ebwt[]", (_verbose || startVerbose));
3215                         if(_verbose || startVerbose) {
3216                                 cerr << "  shared-mem " << (shmemLeader ? "leader" : "follower") << endl;
3217                         }
3218                 } else {
3219                         try {
3220                                 this->_ebwt = new uint8_t[eh->_ebwtTotLen];
3221                         } catch(bad_alloc& e) {
3222                                 cerr << "Out of memory allocating the ebwt[] array for the Bowtie index.  Please try" << endl
3223                                      << "again on a computer with more memory." << endl;
3224                                 throw 1;
3225                         }
3226                 }
3227                 if(shmemLeader) {
3228                         // Read ebwt from primary stream
3229                         MM_READ_RET r = MM_READ(_in1, (void *)this->_ebwt, eh->_ebwtTotLen);
3230                         if(r != (MM_READ_RET)eh->_ebwtTotLen) {
3231                                 cerr << "Error reading ebwt array: returned " << r << ", length was " << (eh->_ebwtTotLen) << endl
3232                                      << "Your index files may be corrupt; please try re-building or re-downloading." << endl
3233                                      << "A complete index consists of 6 files: XYZ.1.ebwt, XYZ.2.ebwt, XYZ.3.ebwt," << endl
3234                                      << "XYZ.4.ebwt, XYZ.rev.1.ebwt, and XYZ.rev.2.ebwt.  The XYZ.1.ebwt and " << endl
3235                                      << "XYZ.rev.1.ebwt files should have the same size, as should the XYZ.2.ebwt and" << endl
3236                                      << "XYZ.rev.2.ebwt files." << endl;
3237                                 throw 1;
3238                         }
3239                         if(switchEndian) {
3240                                 uint8_t *side = this->_ebwt;
3241                                 for(size_t i = 0; i < eh->_numSides; i++) {
3242                                         uint32_t *cums = reinterpret_cast<uint32_t*>(side + eh->_sideSz - 8);
3243                                         cums[0] = endianSwapU32(cums[0]);
3244                                         cums[1] = endianSwapU32(cums[1]);
3245                                         side += this->_eh._sideSz;
3246                                 }
3247                         }
3248                         if(useShmem_) NOTIFY_SHARED(this->_ebwt, eh->_ebwtTotLen);
3249                 } else {
3250                         // Seek past the data and wait until master is finished
3251                         MM_SEEK(_in1, eh->_ebwtTotLen, SEEK_CUR);
3252                         if(useShmem_) WAIT_SHARED(this->_ebwt, eh->_ebwtTotLen);
3253                 }
3254         }
3255
3256         // Read zOff from primary stream
3257         _zOff = readU32(_in1, switchEndian);
3258         bytesRead += 4;
3259         assert_lt(_zOff, len);
3260
3261         try {
3262                 // Read fchr from primary stream
3263                 if(_verbose || startVerbose) cerr << "Reading fchr (5)" << endl;
3264                 if(_useMm) {
3265 #ifdef BOWTIE_MM
3266                         this->_fchr = (uint32_t*)(mmFile[0] + bytesRead);
3267                         bytesRead += 5*4;
3268                         lseek(_in1, 5*4, SEEK_CUR);
3269 #endif
3270                 } else {
3271                         this->_fchr = new uint32_t[5];
3272                         for(int i = 0; i < 5; i++) {
3273                                 this->_fchr[i] = readU32(_in1, switchEndian);
3274                                 assert_leq(this->_fchr[i], len);
3275                                 if(i > 0) assert_geq(this->_fchr[i], this->_fchr[i-1]);
3276                         }
3277                 }
3278                 assert_gt(this->_fchr[4], this->_fchr[0]);
3279                 // Read ftab from primary stream
3280                 if(_verbose || startVerbose) {
3281                         cerr << "Reading ftab (" << eh->_ftabLen << "): ";
3282                         logTime(cerr);
3283                 }
3284                 if(_useMm) {
3285 #ifdef BOWTIE_MM
3286                         this->_ftab = (uint32_t*)(mmFile[0] + bytesRead);
3287                         bytesRead += eh->_ftabLen*4;
3288                         lseek(_in1, eh->_ftabLen*4, SEEK_CUR);
3289 #endif
3290                 } else {
3291                         this->_ftab = new uint32_t[eh->_ftabLen];
3292                         if(switchEndian) {
3293                                 for(uint32_t i = 0; i < eh->_ftabLen; i++)
3294                                         this->_ftab[i] = readU32(_in1, switchEndian);
3295                         } else {
3296                                 MM_READ_RET r = MM_READ(_in1, (void *)this->_ftab, eh->_ftabLen*4);
3297                                 if(r != (MM_READ_RET)(eh->_ftabLen*4)) {
3298                                         cerr << "Error reading _ftab[] array: " << r << ", " << (eh->_ftabLen*4) << endl;
3299                                         throw 1;
3300                                 }
3301                         }
3302                 }
3303                 // Read etab from primary stream
3304                 if(_verbose || startVerbose) {
3305                         cerr << "Reading eftab (" << eh->_eftabLen << "): ";
3306                         logTime(cerr);
3307                 }
3308                 if(_useMm) {
3309 #ifdef BOWTIE_MM
3310                         this->_eftab = (uint32_t*)(mmFile[0] + bytesRead);
3311                         bytesRead += eh->_eftabLen*4;
3312                         lseek(_in1, eh->_eftabLen*4, SEEK_CUR);
3313 #endif
3314                 } else {
3315                         this->_eftab = new uint32_t[eh->_eftabLen];
3316                         if(switchEndian) {
3317                                 for(uint32_t i = 0; i < eh->_eftabLen; i++)
3318                                         this->_eftab[i] = readU32(_in1, switchEndian);
3319                         } else {
3320                                 MM_READ_RET r = MM_READ(_in1, (void *)this->_eftab, eh->_eftabLen*4);
3321                                 if(r != (MM_READ_RET)(eh->_eftabLen*4)) {
3322                                         cerr << "Error reading _eftab[] array: " << r << ", " << (eh->_eftabLen*4) << endl;
3323                                         throw 1;
3324                                 }
3325                         }
3326                 }
3327                 for(uint32_t i = 0; i < eh->_eftabLen; i++) {
3328                         if(i > 0 && this->_eftab[i] > 0) {
3329                                 assert_geq(this->_eftab[i], this->_eftab[i-1]);
3330                         } else if(i > 0 && this->_eftab[i-1] == 0) {
3331                                 assert_eq(0, this->_eftab[i]);
3332                         }
3333                 }
3334         } catch(bad_alloc& e) {
3335                 cerr << "Out of memory allocating fchr[], ftab[] or eftab[] arrays for the Bowtie index." << endl
3336                      << "Please try again on a computer with more memory." << endl;
3337                 throw 1;
3338         }
3339
3340         // Read reference sequence names from primary index file (or not,
3341         // if --refidx is specified)
3342         if(loadNames) {
3343                 while(true) {
3344                         char c = '\0';
3345                         if(MM_READ(_in1, (void *)(&c), (size_t)1) != (MM_READ_RET)1) break;
3346                         bytesRead++;
3347                         if(c == '\0') break;
3348                         else if(c == '\n') {
3349                                 this->_refnames.push_back("");
3350                         } else {
3351                                 if(this->_refnames.size() == 0) {
3352                                         this->_refnames.push_back("");
3353                                 }
3354                                 this->_refnames.back().push_back(c);
3355                         }
3356                 }
3357         }
3358
3359         bytesRead = 4; // reset for secondary index file (already read 1-sentinel)
3360
3361         shmemLeader = true;
3362         if(_verbose || startVerbose) {
3363                 cerr << "Reading offs (" << offsLenSampled << " 32-bit words): ";
3364                 logTime(cerr);
3365         }
3366         if(!_useMm) {
3367                 if(!useShmem_) {
3368                         // Allocate offs_
3369                         try {
3370                                 this->_offs = new uint32_t[offsLenSampled];
3371                         } catch(bad_alloc& e) {
3372                                 cerr << "Out of memory allocating the offs[] array  for the Bowtie index." << endl
3373                                          << "Please try again on a computer with more memory." << endl;
3374                                 throw 1;
3375                         }
3376                 } else {
3377                         shmemLeader = ALLOC_SHARED_U32(
3378                                 (_in2Str + "[offs]"), offsLenSampled*4, &this->_offs,
3379                                 "offs", (_verbose || startVerbose));
3380                 }
3381         }
3382
3383         if(_overrideOffRate < 32) {
3384                 if(shmemLeader) {
3385                         // Allocate offs (big allocation)
3386                         if(switchEndian || offRateDiff > 0) {
3387                                 assert(!_useMm);
3388                                 const uint32_t blockMaxSz = (2 * 1024 * 1024); // 2 MB block size
3389                                 const uint32_t blockMaxSzU32 = (blockMaxSz >> 2); // # U32s per block
3390                                 char *buf = new char[blockMaxSz];
3391                                 for(uint32_t i = 0; i < offsLen; i += blockMaxSzU32) {
3392                                         uint32_t block = min<uint32_t>(blockMaxSzU32, offsLen - i);
3393                                         MM_READ_RET r = MM_READ(_in2, (void *)buf, block << 2);
3394                                         if(r != (MM_READ_RET)(block << 2)) {
3395                                                 cerr << "Error reading block of offs array: " << r << ", " << (block << 2) << endl
3396                                                      << "Your index files may be corrupt; please try re-building or re-downloading." << endl
3397                                                      << "A complete index consists of 6 files: XYZ.1.ebwt, XYZ.2.ebwt, XYZ.3.ebwt," << endl
3398                                                      << "XYZ.4.ebwt, XYZ.rev.1.ebwt, and XYZ.rev.2.ebwt.  The XYZ.1.ebwt and " << endl
3399                                                      << "XYZ.rev.1.ebwt files should have the same size, as should the XYZ.2.ebwt and" << endl
3400                                                      << "XYZ.rev.2.ebwt files." << endl;
3401                                                 throw 1;
3402                                         }
3403                                         uint32_t idx = i >> offRateDiff;
3404                                         for(uint32_t j = 0; j < block; j += (1 << offRateDiff)) {
3405                                                 assert_lt(idx, offsLenSampled);
3406                                                 this->_offs[idx] = ((uint32_t*)buf)[j];
3407                                                 if(switchEndian) {
3408                                                         this->_offs[idx] = endianSwapU32(this->_offs[idx]);
3409                                                 }
3410                                                 idx++;
3411                                         }
3412                                 }
3413                                 delete[] buf;
3414                         } else {
3415                                 if(_useMm) {
3416 #ifdef BOWTIE_MM
3417                                         this->_offs = (uint32_t*)(mmFile[1] + bytesRead);
3418                                         bytesRead += (offsLen << 2);
3419                                         lseek(_in2, (offsLen << 2), SEEK_CUR);
3420 #endif
3421                                 } else {
3422                                         // If any of the high two bits are set
3423                                         if((offsLen & 0xf0000000) != 0) {
3424                                                 if(sizeof(char *) <= 4) {
3425                                                         cerr << "Sanity error: sizeof(char *) <= 4 but offsLen is " << hex << offsLen << endl;
3426                                                         throw 1;
3427                                                 }
3428                                                 // offsLen << 4 overflows sometimes, so do it in four reads
3429                                                 char *offs = (char *)this->_offs;
3430                                                 for(int i = 0; i < 16; i++) {
3431                                                         MM_READ_RET r = MM_READ(_in2, (void*)offs, offsLen >> 2);
3432                                                         if(r != (MM_READ_RET)(offsLen >> 2)) {
3433                                                                 cerr << "Error reading block of _offs[] array: " << r << ", " << (offsLen >> 2) << endl;
3434                                                                 throw 1;
3435                                                         }
3436                                                         offs += (offsLen >> 2);
3437                                                 }
3438                                         } else {
3439                                                 // Do it all in one read
3440                                                 MM_READ_RET r = MM_READ(_in2, (void*)this->_offs, offsLen << 2);
3441                                                 if(r != (MM_READ_RET)(offsLen << 2)) {
3442                                                         cerr << "Error reading _offs[] array: " << r << ", " << (offsLen << 2) << endl;
3443                                                         throw 1;
3444                                                 }
3445                                         }
3446                                 }
3447                         }
3448
3449                         {
3450                                 ASSERT_ONLY(Bitset offsSeen(len+1));
3451                                 for(uint32_t i = 0; i < offsLenSampled; i++) {
3452                                         assert(!offsSeen.test(this->_offs[i]));
3453                                         ASSERT_ONLY(offsSeen.set(this->_offs[i]));
3454                                         assert_leq(this->_offs[i], len);
3455                                 }
3456                         }
3457
3458                         if(useShmem_) NOTIFY_SHARED(this->_offs, offsLenSampled*4);
3459                 } else {
3460                         // Not the shmem leader
3461                         MM_SEEK(_in2, offsLenSampled*4, SEEK_CUR);
3462                         if(useShmem_) WAIT_SHARED(this->_offs, offsLenSampled*4);
3463                 }
3464         }
3465
3466         // Allocate _isa[] (big allocation)
3467         if(_verbose || startVerbose) {
3468                 cerr << "Reading isa (" << isaLenSampled << "): ";
3469                 logTime(cerr);
3470         }
3471         if(!_useMm) {
3472                 try {
3473                         this->_isa = new uint32_t[isaLenSampled];
3474                 } catch(bad_alloc& e) {
3475                         cerr << "Out of memory allocating the isa[] array  for the Bowtie index." << endl
3476                                  << "Please try again on a computer with more memory." << endl;
3477                         throw 1;
3478                 }
3479         }
3480         // Read _isa[]
3481         if(switchEndian || isaRateDiff > 0) {
3482                 assert(!_useMm);
3483                 for(uint32_t i = 0; i < isaLen; i++) {
3484                         if((i & ~(0xffffffff << isaRateDiff)) != 0) {
3485                                 char tmp[4];
3486                                 MM_READ_RET r = MM_READ(_in2, (void *)tmp, 4);
3487                                 if(r != (MM_READ_RET)4) {
3488                                         cerr << "Error reading a word of the _isa[] array: " << r << ", 4" << endl;
3489                                         throw 1;
3490                                 }
3491                         } else {
3492                                 uint32_t idx = i >> isaRateDiff;
3493                                 assert_lt(idx, isaLenSampled);
3494                                 this->_isa[idx] = readU32(_in2, switchEndian);
3495                         }
3496                 }
3497         } else {
3498                 if(_useMm) {
3499 #ifdef BOWTIE_MM
3500                         this->_isa = (uint32_t*)(mmFile[1] + bytesRead);
3501                         bytesRead += (isaLen << 2);
3502                         lseek(_in2, (isaLen << 2), SEEK_CUR);
3503 #endif
3504                 } else {
3505                         MM_READ_RET r = MM_READ(_in2, (void *)this->_isa, isaLen*4);
3506                         if(r != (MM_READ_RET)(isaLen*4)) {
3507                                 cerr << "Error reading _isa[] array: " << r << ", " << (isaLen*4) << endl;
3508                                 throw 1;
3509                         }
3510                 }
3511         }
3512
3513         {
3514                 ASSERT_ONLY(Bitset isasSeen(len+1));
3515                 for(uint32_t i = 0; i < isaLenSampled; i++) {
3516                         assert(!isasSeen.test(this->_isa[i]));
3517                         ASSERT_ONLY(isasSeen.set(this->_isa[i]));
3518                         assert_leq(this->_isa[i], len);
3519                 }
3520         }
3521
3522         this->postReadInit(*eh); // Initialize fields of Ebwt not read from file
3523         if(_verbose || startVerbose) print(cerr, *eh);
3524
3525         // The fact that _ebwt and friends actually point to something
3526         // (other than NULL) now signals to other member functions that the
3527         // Ebwt is loaded into memory.
3528
3529   done: // Exit hatch for both justHeader and !justHeader
3530
3531         // Be kind
3532         if(deleteEh) delete eh;
3533 #ifdef BOWTIE_MM
3534         lseek(_in1, 0, SEEK_SET);
3535         lseek(_in2, 0, SEEK_SET);
3536 #else
3537         rewind(_in1); rewind(_in2);
3538 #endif
3539 }
3540
3541 /**
3542  * Read reference names from an input stream 'in' for an Ebwt primary
3543  * file and store them in 'refnames'.
3544  */
3545 static inline void
3546 readEbwtRefnames(istream& in, vector<string>& refnames) {
3547         // _in1 must already be open with the get cursor at the
3548         // beginning and no error flags set.
3549         assert(in.good());
3550         assert_eq((streamoff)in.tellg(), ios::beg);
3551
3552         // Read endianness hints from both streams
3553         bool switchEndian = false;
3554         uint32_t one = readU32(in, switchEndian); // 1st word of primary stream
3555         if(one != 1) {
3556                 assert_eq((1u<<24), one);
3557                 switchEndian = true;
3558         }
3559
3560         // Reads header entries one by one from primary stream
3561         uint32_t len          = readU32(in, switchEndian);
3562         int32_t  lineRate     = readI32(in, switchEndian);
3563         int32_t  linesPerSide = readI32(in, switchEndian);
3564         int32_t  offRate      = readI32(in, switchEndian);
3565         int32_t  ftabChars    = readI32(in, switchEndian);
3566         // BTL: chunkRate is now deprecated
3567         int32_t flags = readI32(in, switchEndian);
3568         bool color = false;
3569         bool entireReverse = false;
3570         if(flags < 0) {
3571                 color = (((-flags) & EBWT_COLOR) != 0);
3572                 entireReverse = (((-flags) & EBWT_ENTIRE_REV) != 0);
3573         }
3574
3575         // Create a new EbwtParams from the entries read from primary stream
3576         EbwtParams eh(len, lineRate, linesPerSide, offRate, -1, ftabChars, color, entireReverse);
3577
3578         uint32_t nPat = readI32(in, switchEndian); // nPat
3579         in.seekg(nPat*4, ios_base::cur); // skip plen
3580
3581         // Skip rstarts
3582         uint32_t nFrag = readU32(in, switchEndian);
3583         in.seekg(nFrag*4*3, ios_base::cur);
3584
3585         // Skip ebwt
3586         in.seekg(eh._ebwtTotLen, ios_base::cur);
3587
3588         // Skip zOff from primary stream
3589         readU32(in, switchEndian);
3590
3591         // Skip fchr
3592         in.seekg(5 * 4, ios_base::cur);
3593
3594         // Skip ftab
3595         in.seekg(eh._ftabLen*4, ios_base::cur);
3596
3597         // Skip eftab
3598         in.seekg(eh._eftabLen*4, ios_base::cur);
3599
3600         // Read reference sequence names from primary index file
3601         while(true) {
3602                 char c = '\0';
3603                 in.read(&c, 1);
3604                 if(in.eof()) break;
3605                 if(c == '\0') break;
3606                 else if(c == '\n') {
3607                         refnames.push_back("");
3608                 } else {
3609                         if(refnames.size() == 0) {
3610                                 refnames.push_back("");
3611                         }
3612                         refnames.back().push_back(c);
3613                 }
3614         }
3615         if(refnames.back().empty()) {
3616                 refnames.pop_back();
3617         }
3618
3619         // Be kind
3620         in.clear(); in.seekg(0, ios::beg);
3621         assert(in.good());
3622 }
3623
3624 /**
3625  * Read reference names from the index with basename 'in' and store
3626  * them in 'refnames'.
3627  */
3628 static inline void
3629 readEbwtRefnames(const string& instr, vector<string>& refnames) {
3630         ifstream in;
3631         // Initialize our primary and secondary input-stream fields
3632         in.open((instr + ".1.ebwt").c_str(), ios_base::in | ios::binary);
3633         if(!in.is_open()) {
3634                 throw EbwtFileOpenException("Cannot open file " + instr);
3635         }
3636         assert(in.is_open());
3637         assert(in.good());
3638         assert_eq((streamoff)in.tellg(), ios::beg);
3639         readEbwtRefnames(in, refnames);
3640 }
3641
3642 /**
3643  * Read just enough of the Ebwt's header to get its flags
3644  */
3645 static inline int32_t readFlags(const string& instr) {
3646         ifstream in;
3647         // Initialize our primary and secondary input-stream fields
3648         in.open((instr + ".1.ebwt").c_str(), ios_base::in | ios::binary);
3649         if(!in.is_open()) {
3650                 throw EbwtFileOpenException("Cannot open file " + instr);
3651         }
3652         assert(in.is_open());
3653         assert(in.good());
3654         bool switchEndian = false;
3655         uint32_t one = readU32(in, switchEndian); // 1st word of primary stream
3656         if(one != 1) {
3657                 assert_eq((1u<<24), one);
3658                 assert_eq(1, endianSwapU32(one));
3659                 switchEndian = true;
3660         }
3661         readU32(in, switchEndian);
3662         readI32(in, switchEndian);
3663         readI32(in, switchEndian);
3664         readI32(in, switchEndian);
3665         readI32(in, switchEndian);
3666         int32_t flags = readI32(in, switchEndian);
3667         return flags;
3668 }
3669
3670 /**
3671  * Read just enough of the Ebwt's header to determine whether it's
3672  * colorspace.
3673  */
3674 static inline bool
3675 readEbwtColor(const string& instr) {
3676         int32_t flags = readFlags(instr);
3677         if(flags < 0 && (((-flags) & EBWT_COLOR) != 0)) {
3678                 return true;
3679         } else {
3680                 return false;
3681         }
3682 }
3683
3684 /**
3685  * Read just enough of the Ebwt's header to determine whether it's
3686  * entirely reversed.
3687  */
3688 static inline bool
3689 readEntireReverse(const string& instr) {
3690         int32_t flags = readFlags(instr);
3691         if(flags < 0 && (((-flags) & EBWT_ENTIRE_REV) != 0)) {
3692                 return true;
3693         } else {
3694                 return false;
3695         }
3696 }
3697
3698 /**
3699  * Write an extended Burrows-Wheeler transform to a pair of output
3700  * streams.
3701  *
3702  * @param out1 output stream to primary file
3703  * @param out2 output stream to secondary file
3704  * @param be   write in big endian?
3705  */
3706 template<typename TStr>
3707 void Ebwt<TStr>::writeFromMemory(bool justHeader,
3708                                  ostream& out1,
3709                                  ostream& out2) const
3710 {
3711         const EbwtParams& eh = this->_eh;
3712         assert(eh.repOk());
3713         uint32_t be = this->toBe();
3714         assert(out1.good());
3715         assert(out2.good());
3716
3717         // When building an Ebwt, these header parameters are known
3718         // "up-front", i.e., they can be written to disk immediately,
3719         // before we join() or buildToDisk()
3720         writeI32(out1, 1, be); // endian hint for priamry stream
3721         writeI32(out2, 1, be); // endian hint for secondary stream
3722         writeU32(out1, eh._len,          be); // length of string (and bwt and suffix array)
3723         writeI32(out1, eh._lineRate,     be); // 2^lineRate = size in bytes of 1 line
3724         writeI32(out1, eh._linesPerSide, be); // not used
3725         writeI32(out1, eh._offRate,      be); // every 2^offRate chars is "marked"
3726         writeI32(out1, eh._ftabChars,    be); // number of 2-bit chars used to address ftab
3727         int32_t flags = 1;
3728         if(eh._color) flags |= EBWT_COLOR;
3729         if(eh._entireReverse) flags |= EBWT_ENTIRE_REV;
3730         writeI32(out1, -flags, be); // BTL: chunkRate is now deprecated
3731
3732         if(!justHeader) {
3733                 assert(isInMemory());
3734                 // These Ebwt parameters are known after the inputs strings have
3735                 // been joined() but before they have been built().  These can
3736                 // written to the disk next and then discarded from memory.
3737                 writeU32(out1, this->_nPat,      be);
3738                 for(uint32_t i = 0; i < this->_nPat; i++)
3739                 writeU32(out1, this->_plen[i], be);
3740                 assert_geq(this->_nFrag, this->_nPat);
3741                 writeU32(out1, this->_nFrag, be);
3742                 for(uint32_t i = 0; i < this->_nFrag*3; i++)
3743                         writeU32(out1, this->_rstarts[i], be);
3744
3745                 // These Ebwt parameters are discovered only as the Ebwt is being
3746                 // built (in buildToDisk()).  Of these, only 'offs' and 'ebwt' are
3747                 // terribly large.  'ebwt' is written to the primary file and then
3748                 // discarded from memory as it is built; 'offs' is similarly
3749                 // written to the secondary file and discarded.
3750                 out1.write((const char *)this->ebwt(), eh._ebwtTotLen);
3751                 writeU32(out1, this->zOff(), be);
3752                 uint32_t offsLen = eh._offsLen;
3753                 for(uint32_t i = 0; i < offsLen; i++)
3754                         writeU32(out2, this->_offs[i], be);
3755                 uint32_t isaLen = eh._isaLen;
3756                 for(uint32_t i = 0; i < isaLen; i++)
3757                         writeU32(out2, this->_isa[i], be);
3758
3759                 // 'fchr', 'ftab' and 'eftab' are not fully determined until the
3760                 // loop is finished, so they are written to the primary file after
3761                 // all of 'ebwt' has already been written and only then discarded
3762                 // from memory.
3763                 for(int i = 0; i < 5; i++)
3764                         writeU32(out1, this->_fchr[i], be);
3765                 for(uint32_t i = 0; i < eh._ftabLen; i++)
3766                         writeU32(out1, this->ftab()[i], be);
3767                 for(uint32_t i = 0; i < eh._eftabLen; i++)
3768                         writeU32(out1, this->eftab()[i], be);
3769         }
3770 }
3771
3772 /**
3773  * Given a pair of strings representing output filenames, and assuming
3774  * this Ebwt object is currently in memory, write out this Ebwt to the
3775  * specified files.
3776  *
3777  * If sanity-checking is enabled, then once the streams have been
3778  * fully written and closed, we reopen them and read them into a
3779  * (hopefully) exact copy of this Ebwt.  We then assert that the
3780  * current Ebwt and the copy match in all of their fields.
3781  */
3782 template<typename TStr>
3783 void Ebwt<TStr>::writeFromMemory(bool justHeader,
3784                                  const string& out1,
3785                                  const string& out2) const
3786 {
3787         const EbwtParams& eh = this->_eh;
3788         assert(isInMemory());
3789         assert(eh.repOk());
3790
3791         ofstream fout1(out1.c_str(), ios::binary);
3792         ofstream fout2(out2.c_str(), ios::binary);
3793         writeFromMemory(justHeader, fout1, fout2);
3794         fout1.close();
3795         fout2.close();
3796
3797         // Read the file back in and assert that all components match
3798         if(_sanity) {
3799                 if(_verbose)
3800                         cout << "Re-reading \"" << out1 << "\"/\"" << out2 << "\" for sanity check" << endl;
3801                 Ebwt copy(out1, out2, _verbose, _sanity);
3802                 assert(!isInMemory());
3803                 copy.loadIntoMemory(eh._color ? 1 : 0, false, false);
3804                 assert(isInMemory());
3805             assert_eq(eh._lineRate,     copy.eh()._lineRate);
3806             assert_eq(eh._linesPerSide, copy.eh()._linesPerSide);
3807             assert_eq(eh._offRate,      copy.eh()._offRate);
3808             assert_eq(eh._isaRate,      copy.eh()._isaRate);
3809             assert_eq(eh._ftabChars,    copy.eh()._ftabChars);
3810             assert_eq(eh._len,          copy.eh()._len);
3811             assert_eq(_zOff,             copy.zOff());
3812             assert_eq(_zEbwtBpOff,       copy.zEbwtBpOff());
3813             assert_eq(_zEbwtByteOff,     copy.zEbwtByteOff());
3814                 assert_eq(_nPat,             copy.nPat());
3815                 for(uint32_t i = 0; i < _nPat; i++)
3816                         assert_eq(this->_plen[i], copy.plen()[i]);
3817                 assert_eq(this->_nFrag, copy.nFrag());
3818                 for(uint32_t i = 0; i < this->nFrag*3; i++) {
3819                         assert_eq(this->_rstarts[i], copy.rstarts()[i]);
3820                 }
3821                 for(uint32_t i = 0; i < 5; i++)
3822                         assert_eq(this->_fchr[i], copy.fchr()[i]);
3823                 for(uint32_t i = 0; i < eh._ftabLen; i++)
3824                         assert_eq(this->ftab()[i], copy.ftab()[i]);
3825                 for(uint32_t i = 0; i < eh._eftabLen; i++)
3826                         assert_eq(this->eftab()[i], copy.eftab()[i]);
3827                 for(uint32_t i = 0; i < eh._offsLen; i++)
3828                         assert_eq(this->_offs[i], copy.offs()[i]);
3829                 for(uint32_t i = 0; i < eh._isaLen; i++)
3830                         assert_eq(this->_isa[i], copy.isa()[i]);
3831                 for(uint32_t i = 0; i < eh._ebwtTotLen; i++)
3832                         assert_eq(this->ebwt()[i], copy.ebwt()[i]);
3833                 copy.sanityCheckAll();
3834                 if(_verbose)
3835                         cout << "Read-in check passed for \"" << out1 << "\"/\"" << out2 << "\"" << endl;
3836         }
3837 }
3838
3839 ///////////////////////////////////////////////////////////////////////
3840 //
3841 // Functions for building Ebwts
3842 //
3843 ///////////////////////////////////////////////////////////////////////
3844
3845 /**
3846  * Join several text strings together in a way that's compatible with
3847  * the text-chunking scheme dictated by chunkRate parameter.
3848  *
3849  * The non-static member Ebwt::join additionally builds auxilliary
3850  * arrays that maintain a mapping between chunks in the joined string
3851  * and the original text strings.
3852  */
3853 template<typename TStr>
3854 TStr Ebwt<TStr>::join(vector<TStr>& l, uint32_t seed) {
3855         RandomSource rand; // reproducible given same seed
3856         rand.init(seed);
3857         TStr ret;
3858         size_t guessLen = 0;
3859         for(size_t i = 0; i < l.size(); i++) {
3860                 guessLen += length(l[i]);
3861         }
3862         reserve(ret, guessLen, Exact());
3863         for(size_t i = 0; i < l.size(); i++) {
3864                 TStr& s = l[i];
3865                 assert_gt(length(s), 0);
3866                 append(ret, s);
3867         }
3868         return ret;
3869 }
3870
3871 /**
3872  * Join several text strings together in a way that's compatible with
3873  * the text-chunking scheme dictated by chunkRate parameter.
3874  *
3875  * The non-static member Ebwt::join additionally builds auxilliary
3876  * arrays that maintain a mapping between chunks in the joined string
3877  * and the original text strings.
3878  */
3879 template<typename TStr>
3880 TStr Ebwt<TStr>::join(vector<FileBuf*>& l,
3881                       vector<RefRecord>& szs,
3882                       uint32_t sztot,
3883                       const RefReadInParams& refparams,
3884                       uint32_t seed)
3885 {
3886         RandomSource rand; // reproducible given same seed
3887         rand.init(seed);
3888         RefReadInParams rpcp = refparams;
3889         TStr ret;
3890         size_t guessLen = sztot;
3891         reserve(ret, guessLen, Exact());
3892         ASSERT_ONLY(size_t szsi = 0);
3893         for(size_t i = 0; i < l.size(); i++) {
3894                 // For each sequence we can pull out of istream l[i]...
3895                 assert(!l[i]->eof());
3896                 bool first = true;
3897                 while(!l[i]->eof()) {
3898                         RefRecord rec = fastaRefReadAppend(*l[i], first, ret, rpcp);
3899 #ifndef ACCOUNT_FOR_ALL_GAP_REFS
3900                         if(rec.first && rec.len == 0) rec.first = false;
3901 #endif
3902                         first = false;
3903                         size_t bases = rec.len;
3904                         assert_eq(rec.off, szs[szsi].off);
3905                         assert_eq(rec.len, szs[szsi].len);
3906                         assert_eq(rec.first, szs[szsi].first);
3907                         ASSERT_ONLY(szsi++);
3908                         if(bases == 0) continue;
3909                 }
3910         }
3911         return ret;
3912 }
3913
3914 /**
3915  * Join several text strings together according to the text-chunking
3916  * scheme specified in the EbwtParams.  Ebwt fields calculated in this
3917  * function are written directly to disk.
3918  *
3919  * It is assumed, but not required, that the header values have already
3920  * been written to 'out1' before this function is called.
3921  *
3922  * The static member Ebwt::join just returns a joined version of a
3923  * list of strings without building any of the auxiliary arrays.
3924  * Because the pseudo-random number generator is the same, we expect
3925  * this function and the static function to give the same result given
3926  * the same seed.
3927  */
3928 template<typename TStr>
3929 void Ebwt<TStr>::joinToDisk(
3930         vector<FileBuf*>& l,
3931         vector<RefRecord>& szs,
3932         vector<uint32_t>& plens,
3933         uint32_t sztot,
3934         const RefReadInParams& refparams,
3935         TStr& ret,
3936         ostream& out1,
3937         ostream& out2,
3938         uint32_t seed = 0)
3939 {
3940         RandomSource rand; // reproducible given same seed
3941         rand.init(seed);
3942         RefReadInParams rpcp = refparams;
3943         assert_gt(szs.size(), 0);
3944         assert_gt(l.size(), 0);
3945         assert_gt(sztot, 0);
3946         // Not every fragment represents a distinct sequence - many
3947         // fragments may correspond to a single sequence.  Count the
3948         // number of sequences here by counting the number of "first"
3949         // fragments.
3950         this->_nPat = 0;
3951         this->_nFrag = 0;
3952 #ifdef ACCOUNT_FOR_ALL_GAP_REFS
3953         int nGapFrag = 0;
3954 #endif
3955         for(size_t i = 0; i < szs.size(); i++) {
3956                 if(szs[i].len > 0) this->_nFrag++;
3957 #ifdef ACCOUNT_FOR_ALL_GAP_REFS
3958                 if(szs[i].len == 0 && szs[i].off > 0) nGapFrag++;
3959                 if(szs[i].first && szs[i].len > 0) this->_nPat++;
3960 #else
3961                 // For all records where len=0 and first=1, set first=0
3962                 assert(szs[i].len > 0 || !szs[i].first);
3963                 if(szs[i].first) this->_nPat++;
3964 #endif
3965         }
3966         assert_gt(this->_nPat, 0);
3967         assert_geq(this->_nFrag, this->_nPat);
3968         this->_rstarts = NULL;
3969         writeU32(out1, this->_nPat, this->toBe());
3970         assert_eq(plens.size(), this->_nPat);
3971         // Allocate plen[]
3972         try {
3973                 this->_plen = new uint32_t[this->_nPat];
3974         } catch(bad_alloc& e) {
3975                 cerr << "Out of memory allocating plen[] in Ebwt::join()"
3976                      << " at " << __FILE__ << ":" << __LINE__ << endl;
3977                 throw e;
3978         }
3979         // For each pattern, set plen
3980         for(size_t i = 0; i < plens.size(); i++) {
3981                 this->_plen[i] = plens[i];
3982                 writeU32(out1, this->_plen[i], this->toBe());
3983         }
3984         // Write the number of fragments
3985         writeU32(out1, this->_nFrag, this->toBe());
3986         size_t seqsRead = 0;
3987         ASSERT_ONLY(uint32_t szsi = 0);
3988         ASSERT_ONLY(uint32_t entsWritten = 0);
3989         // For each filebuf
3990         for(unsigned int i = 0; i < l.size(); i++) {
3991                 assert(!l[i]->eof());
3992                 bool first = true;
3993                 uint32_t patoff = 0;
3994                 // For each *fragment* (not necessary an entire sequence) we
3995                 // can pull out of istream l[i]...
3996                 while(!l[i]->eof()) {
3997                         string name;
3998                         // Push a new name onto our vector
3999                         _refnames.push_back("");
4000                         //uint32_t oldRetLen = length(ret);
4001                         RefRecord rec = fastaRefReadAppend(*l[i], first, ret, rpcp, &_refnames.back());
4002 #ifndef ACCOUNT_FOR_ALL_GAP_REFS
4003                         if(rec.first && rec.len == 0) rec.first = false;
4004 #endif
4005                         first = false;
4006                         if(rec.first) {
4007                                 if(_refnames.back().length() == 0) {
4008                                         // If name was empty, replace with an index
4009                                         ostringstream stm;
4010                                         stm << (_refnames.size()-1);
4011                                         _refnames.back() = stm.str();
4012                                 }
4013                         } else {
4014                                 // This record didn't actually start a new sequence so
4015                                 // no need to add a name
4016                                 //assert_eq(0, _refnames.back().length());
4017                                 _refnames.pop_back();
4018                         }
4019                         assert_lt(szsi, szs.size());
4020                         assert(szs[szsi].first == 0 || szs[szsi].first == 1);
4021                         assert_eq(rec.off, szs[szsi].off);
4022                         assert_eq(rec.len, szs[szsi].len);
4023                         // szs[szsi].first == 2 sometimes?!?!  g++ is unable to do
4024                         // the following correctly, regardless of how I write it
4025                         //assert((rec.first == 0) == (szs[szsi].first == 0));
4026                         assert(rec.first || rec.off > 0);
4027                         ASSERT_ONLY(szsi++);
4028 #ifdef ACCOUNT_FOR_ALL_GAP_REFS
4029                         if(rec.len == 0) continue;
4030                         if(rec.first && rec.len > 0) seqsRead++;
4031                         assert_leq(rec.len, this->_plen[seqsRead-1]);
4032 #else
4033                         if(rec.first) seqsRead++;
4034                         if(rec.len == 0) continue;
4035                         assert_leq(rec.len, this->_plen[seqsRead-1]);
4036 #endif
4037                         // Reset the patoff if this is the first fragment
4038                         if(rec.first) patoff = 0;
4039                         patoff += rec.off; // add fragment's offset from end of last frag.
4040                         // Adjust rpcps
4041                         //uint32_t seq = seqsRead-1;
4042                         ASSERT_ONLY(entsWritten++);
4043                         // This is where rstarts elements are written to the output stream
4044                         //writeU32(out1, oldRetLen, this->toBe()); // offset from beginning of joined string
4045                         //writeU32(out1, seq,       this->toBe()); // sequence id
4046                         //writeU32(out1, patoff,    this->toBe()); // offset into sequence
4047                         patoff += rec.len;
4048                 }
4049                 assert_gt(szsi, 0);
4050                 l[i]->reset();
4051                 assert(!l[i]->eof());
4052                 #ifndef NDEBUG
4053                 int c = l[i]->get();
4054                 assert_eq('>', c);
4055                 assert(!l[i]->eof());
4056                 l[i]->reset();
4057                 assert(!l[i]->eof());
4058                 #endif
4059         }
4060         assert_eq(entsWritten, this->_nFrag);
4061 }
4062
4063
4064 /**
4065  * Build an Ebwt from a string 's' and its suffix array 'sa' (which
4066  * might actually be a suffix array *builder* that builds blocks of the
4067  * array on demand).  The bulk of the Ebwt, i.e. the ebwt and offs
4068  * arrays, is written directly to disk.  This is by design: keeping
4069  * those arrays in memory needlessly increases the footprint of the
4070  * building process.  Instead, we prefer to build the Ebwt directly
4071  * "to disk" and then read it back into memory later as necessary.
4072  *
4073  * It is assumed that the header values and join-related values (nPat,
4074  * plen) have already been written to 'out1' before this function
4075  * is called.  When this function is finished, it will have
4076  * additionally written ebwt, zOff, fchr, ftab and eftab to the primary
4077  * file and offs to the secondary file.
4078  *
4079  * Assume DNA/RNA/any alphabet with 4 or fewer elements.
4080  * Assume occ array entries are 32 bits each.
4081  *
4082  * @param sa            the suffix array to convert to a Ebwt
4083  * @param buildISA      whether to output an ISA sample into out2 after
4084  *                      the SA sample
4085  * @param s             the original string
4086  * @param out
4087  */
4088 template<typename TStr>
4089 void Ebwt<TStr>::buildToDisk(InorderBlockwiseSA<TStr>& sa,
4090                              const TStr& s,
4091                              ostream& out1,
4092                              ostream& out2)
4093 {
4094         const EbwtParams& eh = this->_eh;
4095
4096         assert(eh.repOk());
4097         assert_eq(length(s)+1, sa.size());
4098         assert_eq(length(s), eh._len);
4099         assert_gt(eh._lineRate, 3);
4100         assert(sa.suffixItrIsReset());
4101         assert_leq((int)ValueSize<Dna>::VALUE, 4);
4102
4103         uint32_t  len = eh._len;
4104         uint32_t  ftabLen = eh._ftabLen;
4105         uint32_t  sideSz = eh._sideSz;
4106         uint32_t  ebwtTotSz = eh._ebwtTotSz;
4107         uint32_t  fchr[] = {0, 0, 0, 0, 0};
4108         uint32_t* ftab = NULL;
4109         uint32_t  zOff = 0xffffffff;
4110
4111         // Save # of occurrences of each character as we walk along the bwt
4112         uint32_t occ[4] = {0, 0, 0, 0};
4113         // Save 'G' and 'T' occurrences between backward and forward buckets
4114         uint32_t occSave[2] = {0, 0};
4115
4116         // Record rows that should "absorb" adjacent rows in the ftab.
4117         // The absorbed rows represent suffixes shorter than the ftabChars
4118         // cutoff.
4119         uint8_t absorbCnt = 0;
4120         uint8_t *absorbFtab;
4121         try {
4122                 VMSG_NL("Allocating ftab, absorbFtab");
4123                 ftab = new uint32_t[ftabLen];
4124                 memset(ftab, 0, 4 * ftabLen);
4125                 absorbFtab = new uint8_t[ftabLen];
4126                 memset(absorbFtab, 0, ftabLen);
4127         } catch(bad_alloc &e) {
4128                 cerr << "Out of memory allocating ftab[] or absorbFtab[] "
4129                      << "in Ebwt::buildToDisk() at " << __FILE__ << ":"
4130                      << __LINE__ << endl;
4131                 throw e;
4132         }
4133         assert(ftab != NULL);
4134         assert(absorbFtab != NULL);
4135
4136         // Allocate the side buffer; holds a single side as its being
4137         // constructed and then written to disk.  Reused across all sides.
4138 #ifdef SIXTY4_FORMAT
4139         uint64_t *ebwtSide = NULL;
4140 #else
4141         uint8_t *ebwtSide = NULL;
4142 #endif
4143         try {
4144 #ifdef SIXTY4_FORMAT
4145                 ebwtSide = new uint64_t[sideSz >> 3];
4146 #else
4147                 ebwtSide = new uint8_t[sideSz];
4148 #endif
4149         } catch(bad_alloc &e) {
4150                 cerr << "Out of memory allocating ebwtSide[] in "
4151                      << "Ebwt::buildToDisk() at " << __FILE__ << ":"
4152                      << __LINE__ << endl;
4153                 throw e;
4154         }
4155         assert(ebwtSide != NULL);
4156
4157         // Allocate a buffer to hold the ISA sample, which we accumulate in
4158         // the loop and then output at the end.  We can't write output the
4159         // ISA right away because the order in which we calculate its
4160         // elements is based on the suffix array, which we only see bit by
4161         // bit
4162         uint32_t *isaSample = NULL;
4163         if(eh._isaRate >= 0) {
4164                 try {
4165                         isaSample = new uint32_t[eh._isaLen];
4166                 } catch(bad_alloc &e) {
4167                         cerr << "Out of memory allocating isaSample[] in "
4168                              << "Ebwt::buildToDisk() at " << __FILE__ << ":"
4169                              << __LINE__ << endl;
4170                         throw e;
4171                 }
4172                 assert(isaSample != NULL);
4173         }
4174
4175         // Points to the base offset within ebwt for the side currently
4176         // being written
4177         uint32_t side = 0;
4178         // Points to a byte offset from 'side' within ebwt[] where next
4179         // char should be written
4180 #ifdef SIXTY4_FORMAT
4181         int sideCur = (eh._sideBwtSz >> 3) - 1;
4182 #else
4183         int sideCur = eh._sideBwtSz - 1;
4184 #endif
4185
4186         // Whether we're assembling a forward or a reverse bucket
4187         bool fw = false;
4188
4189         // Did we just finish writing a forward bucket?  (Must be true when
4190         // we exit the loop.)
4191         bool wroteFwBucket = false;
4192
4193         // Have we skipped the '$' in the last column yet?
4194         ASSERT_ONLY(bool dollarSkipped = false);
4195
4196         uint32_t si = 0;   // string offset (chars)
4197         ASSERT_ONLY(uint32_t lastSufInt = 0);
4198         ASSERT_ONLY(bool inSA = true); // true iff saI still points inside suffix
4199                                        // array (as opposed to the padding at the
4200                                        // end)
4201         // Iterate over packed bwt bytes
4202         VMSG_NL("Entering Ebwt loop");
4203         ASSERT_ONLY(uint32_t beforeEbwtOff = (uint32_t)out1.tellp());
4204         while(side < ebwtTotSz) {
4205                 wroteFwBucket = false;
4206                 // Sanity-check our cursor into the side buffer
4207                 assert_geq(sideCur, 0);
4208                 assert_lt(sideCur, (int)eh._sideBwtSz);
4209                 assert_eq(0, side % sideSz); // 'side' must be on side boundary
4210                 ebwtSide[sideCur] = 0; // clear
4211                 assert_lt(side + sideCur, ebwtTotSz);
4212                 // Iterate over bit-pairs in the si'th character of the BWT
4213 #ifdef SIXTY4_FORMAT
4214                 for(int bpi = 0; bpi < 32; bpi++, si++)
4215 #else
4216                 for(int bpi = 0; bpi < 4; bpi++, si++)
4217 #endif
4218                 {
4219                         int bwtChar;
4220                         bool count = true;
4221                         if(si <= len) {
4222                                 // Still in the SA; extract the bwtChar
4223                                 uint32_t saElt = sa.nextSuffix();
4224                                 // (that might have triggered sa to calc next suf block)
4225                                 if(isaSample != NULL && (saElt & eh._isaMask) == saElt) {
4226                                         // This element belongs in the ISA sample.  Add
4227                                         // an entry mapping the text offset to the offset
4228                                         // into the suffix array that holds the suffix
4229                                         // beginning with the character at that text offset
4230                                         assert_lt((saElt >> eh._isaRate), eh._isaLen);
4231                                         isaSample[saElt >> eh._isaRate] = si;
4232                                 }
4233                                 if(saElt == 0) {
4234                                         // Don't add the '$' in the last column to the BWT
4235                                         // transform; we can't encode a $ (only A C T or G)
4236                                         // and counting it as, say, an A, will mess up the
4237                                         // LR mapping
4238                                         bwtChar = 0; count = false;
4239                                         ASSERT_ONLY(dollarSkipped = true);
4240                                         zOff = si; // remember the SA row that
4241                                                    // corresponds to the 0th suffix
4242                                 } else {
4243                                         bwtChar = (int)(Dna)(s[saElt-1]);
4244                                         assert_lt(bwtChar, 4);
4245                                         // Update the fchr
4246                                         fchr[bwtChar]++;
4247                                 }
4248                                 // Update ftab
4249                                 if((len-saElt) >= (uint32_t)eh._ftabChars) {
4250                                         // Turn the first ftabChars characters of the
4251                                         // suffix into an integer index into ftab
4252                                         uint32_t sufInt = 0;
4253                                         for(int i = 0; i < eh._ftabChars; i++) {
4254                                                 sufInt <<= 2;
4255                                                 assert_lt(i, (int)(len-saElt));
4256                                                 sufInt |= (unsigned char)(Dna)(s[saElt+i]);
4257                                         }
4258                                         // Assert that this prefix-of-suffix is greater
4259                                         // than or equal to the last one (true b/c the
4260                                         // suffix array is sorted)
4261                                         #ifndef NDEBUG
4262                                         if(lastSufInt > 0) assert_geq(sufInt, lastSufInt);
4263                                         lastSufInt = sufInt;
4264                                         #endif
4265                                         // Update ftab
4266                                         assert_lt(sufInt+1, ftabLen);
4267                                         ftab[sufInt+1]++;
4268                                         if(absorbCnt > 0) {
4269                                                 // Absorb all short suffixes since the last
4270                                                 // transition into this transition
4271                                                 absorbFtab[sufInt] = absorbCnt;
4272                                                 absorbCnt = 0;
4273                                         }
4274                                 } else {
4275                                         // Otherwise if suffix is fewer than ftabChars
4276                                         // characters long, then add it to the 'absorbCnt';
4277                                         // it will be absorbed into the next transition
4278                                         assert_lt(absorbCnt, 255);
4279                                         absorbCnt++;
4280                                 }
4281                                 // Suffix array offset boundary? - update offset array
4282                                 if((si & eh._offMask) == si) {
4283                                         assert_lt((si >> eh._offRate), eh._offsLen);
4284                                         // Write offsets directly to the secondary output
4285                                         // stream, thereby avoiding keeping them in memory
4286                                         writeU32(out2, saElt, this->toBe());
4287                                 }
4288                         } else {
4289                                 // Strayed off the end of the SA, now we're just
4290                                 // padding out a bucket
4291                                 #ifndef NDEBUG
4292                                 if(inSA) {
4293                                         // Assert that we wrote all the characters in the
4294                                         // string before now
4295                                         assert_eq(si, len+1);
4296                                         inSA = false;
4297                                 }
4298                                 #endif
4299                                 // 'A' used for padding; important that padding be
4300                                 // counted in the occ[] array
4301                                 bwtChar = 0;
4302                         }
4303                         if(count) occ[bwtChar]++;
4304                         // Append BWT char to bwt section of current side
4305                         if(fw) {
4306                                 // Forward bucket: fill from least to most
4307 #ifdef SIXTY4_FORMAT
4308                                 ebwtSide[sideCur] |= ((uint64_t)bwtChar << (bpi << 1));
4309                                 if(bwtChar > 0) assert_gt(ebwtSide[sideCur], 0);
4310 #else
4311                                 pack_2b_in_8b(bwtChar, ebwtSide[sideCur], bpi);
4312                                 assert_eq((ebwtSide[sideCur] >> (bpi*2)) & 3, bwtChar);
4313 #endif
4314                         } else {
4315                                 // Backward bucket: fill from most to least
4316 #ifdef SIXTY4_FORMAT
4317                                 ebwtSide[sideCur] |= ((uint64_t)bwtChar << ((31 - bpi) << 1));
4318                                 if(bwtChar > 0) assert_gt(ebwtSide[sideCur], 0);
4319 #else
4320                                 pack_2b_in_8b(bwtChar, ebwtSide[sideCur], 3-bpi);
4321                                 assert_eq((ebwtSide[sideCur] >> ((3-bpi)*2)) & 3, bwtChar);
4322 #endif
4323                         }
4324                 } // end loop over bit-pairs
4325                 assert_eq(dollarSkipped ? 3 : 0, (occ[0] + occ[1] + occ[2] + occ[3]) & 3);
4326 #ifdef SIXTY4_FORMAT
4327                 assert_eq(0, si & 31);
4328 #else
4329                 assert_eq(0, si & 3);
4330 #endif
4331                 if(fw) sideCur++;
4332                 else   sideCur--;
4333 #ifdef SIXTY4_FORMAT
4334                 if(sideCur == (int)eh._sideBwtSz >> 3)
4335 #else
4336                 if(sideCur == (int)eh._sideBwtSz)
4337 #endif
4338                 {
4339                         // Forward side boundary
4340                         assert_eq(0, si % eh._sideBwtLen);
4341 #ifdef SIXTY4_FORMAT
4342                         sideCur = (eh._sideBwtSz >> 3) - 1;
4343 #else
4344                         sideCur = eh._sideBwtSz - 1;
4345 #endif
4346                         assert(fw); fw = false; wroteFwBucket = true;
4347                         // Write 'G' and 'T'
4348                         assert_leq(occSave[0], occ[2]);
4349                         assert_leq(occSave[1], occ[3]);
4350                         uint32_t *u32side = reinterpret_cast<uint32_t*>(ebwtSide);
4351                         side += sideSz;
4352                         assert_leq(side, eh._ebwtTotSz);
4353                         u32side[(sideSz >> 2)-2] = endianizeU32(occSave[0], this->toBe());
4354                         u32side[(sideSz >> 2)-1] = endianizeU32(occSave[1], this->toBe());
4355                         // Write forward side to primary file
4356                         out1.write((const char *)ebwtSide, sideSz);
4357                 } else if (sideCur == -1) {
4358                         // Backward side boundary
4359                         assert_eq(0, si % eh._sideBwtLen);
4360                         sideCur = 0;
4361                         assert(!fw); fw = true;
4362                         // Write 'A' and 'C'
4363                         uint32_t *u32side = reinterpret_cast<uint32_t*>(ebwtSide);
4364                         side += sideSz;
4365                         assert_leq(side, eh._ebwtTotSz);
4366                         u32side[(sideSz >> 2)-2] = endianizeU32(occ[0], this->toBe());
4367                         u32side[(sideSz >> 2)-1] = endianizeU32(occ[1], this->toBe());
4368                         occSave[0] = occ[2]; // save 'G' count
4369                         occSave[1] = occ[3]; // save 'T' count
4370                         // Write backward side to primary file
4371                         out1.write((const char *)ebwtSide, sideSz);
4372                 }
4373         }
4374         VMSG_NL("Exited Ebwt loop");
4375         assert(ftab != NULL);
4376         assert_neq(zOff, 0xffffffff);
4377         if(absorbCnt > 0) {
4378                 // Absorb any trailing, as-yet-unabsorbed short suffixes into
4379                 // the last element of ftab
4380                 absorbFtab[ftabLen-1] = absorbCnt;
4381         }
4382         // Assert that our loop counter got incremented right to the end
4383         assert_eq(side, eh._ebwtTotSz);
4384         // Assert that we wrote the expected amount to out1
4385         assert_eq(((uint32_t)out1.tellp() - beforeEbwtOff), eh._ebwtTotSz);
4386         // assert that the last thing we did was write a forward bucket
4387         assert(wroteFwBucket);
4388
4389         //
4390         // Write zOff to primary stream
4391         //
4392         writeU32(out1, zOff, this->toBe());
4393
4394         //
4395         // Finish building fchr
4396         //
4397         // Exclusive prefix sum on fchr
4398         for(int i = 1; i < 4; i++) {
4399                 fchr[i] += fchr[i-1];
4400         }
4401         assert_eq(fchr[3], len);
4402         // Shift everybody up by one
4403         for(int i = 4; i >= 1; i--) {
4404                 fchr[i] = fchr[i-1];
4405         }
4406         fchr[0] = 0;
4407         if(_verbose) {
4408                 for(int i = 0; i < 5; i++)
4409                         cout << "fchr[" << "ACGT$"[i] << "]: " << fchr[i] << endl;
4410         }
4411         // Write fchr to primary file
4412         for(int i = 0; i < 5; i++) {
4413                 writeU32(out1, fchr[i], this->toBe());
4414         }
4415
4416         //
4417         // Finish building ftab and build eftab
4418         //
4419         // Prefix sum on ftable
4420         uint32_t eftabLen = 0;
4421         assert_eq(0, absorbFtab[0]);
4422         for(uint32_t i = 1; i < ftabLen; i++) {
4423                 if(absorbFtab[i] > 0) eftabLen += 2;
4424         }
4425         assert_leq(eftabLen, (uint32_t)eh._ftabChars*2);
4426         eftabLen = eh._ftabChars*2;
4427         uint32_t *eftab = NULL;
4428         try {
4429                 eftab = new uint32_t[eftabLen];
4430                 memset(eftab, 0, 4 * eftabLen);
4431         } catch(bad_alloc &e) {
4432                 cerr << "Out of memory allocating eftab[] "
4433                      << "in Ebwt::buildToDisk() at " << __FILE__ << ":"
4434                      << __LINE__ << endl;
4435                 throw e;
4436         }
4437         assert(eftab != NULL);
4438         uint32_t eftabCur = 0;
4439         for(uint32_t i = 1; i < ftabLen; i++) {
4440                 uint32_t lo = ftab[i] + Ebwt::ftabHi(ftab, eftab, len, ftabLen, eftabLen, i-1);
4441                 if(absorbFtab[i] > 0) {
4442                         // Skip a number of short pattern indicated by absorbFtab[i]
4443                         uint32_t hi = lo + absorbFtab[i];
4444                         assert_lt(eftabCur*2+1, eftabLen);
4445                         eftab[eftabCur*2] = lo;
4446                         eftab[eftabCur*2+1] = hi;
4447                         ftab[i] = (eftabCur++) ^ 0xffffffff; // insert pointer into eftab
4448                         assert_eq(lo, Ebwt::ftabLo(ftab, eftab, len, ftabLen, eftabLen, i));
4449                         assert_eq(hi, Ebwt::ftabHi(ftab, eftab, len, ftabLen, eftabLen, i));
4450                 } else {
4451                         ftab[i] = lo;
4452                 }
4453         }
4454         assert_eq(Ebwt::ftabHi(ftab, eftab, len, ftabLen, eftabLen, ftabLen-1), len+1);
4455         // Write ftab to primary file
4456         for(uint32_t i = 0; i < ftabLen; i++) {
4457                 writeU32(out1, ftab[i], this->toBe());
4458         }
4459         // Write eftab to primary file
4460         for(uint32_t i = 0; i < eftabLen; i++) {
4461                 writeU32(out1, eftab[i], this->toBe());
4462         }
4463         // Write isa to primary file
4464         if(isaSample != NULL) {
4465                 ASSERT_ONLY(Bitset sawISA(eh._len+1));
4466                 for(uint32_t i = 0; i < eh._isaLen; i++) {
4467                         uint32_t s = isaSample[i];
4468                         assert_leq(s, eh._len);
4469                         assert(!sawISA.test(s));
4470                         ASSERT_ONLY(sawISA.set(s));
4471                         writeU32(out2, s, this->toBe());
4472                 }
4473                 delete[] isaSample;
4474         }
4475         delete[] ftab;
4476         delete[] eftab;
4477         delete[] absorbFtab;
4478
4479         // Note: if you'd like to sanity-check the Ebwt, you'll have to
4480         // read it back into memory first!
4481         assert(!isInMemory());
4482         VMSG_NL("Exiting Ebwt::buildToDisk()");
4483 }
4484
4485 /**
4486  * Try to find the Bowtie index specified by the user.  First try the
4487  * exact path given by the user.  Then try the user-provided string
4488  * appended onto the path of the "indexes" subdirectory below this
4489  * executable, then try the provided string appended onto
4490  * "$BOWTIE_INDEXES/".
4491  */
4492 string adjustEbwtBase(const string& cmdline,
4493                                           const string& ebwtFileBase,
4494                                           bool verbose = false);
4495
4496 #endif /*EBWT_H_*/