#ifndef EBWT_SEARCH_BACKTRACK_H_
#define EBWT_SEARCH_BACKTRACK_H_

#include <stdexcept>
#include <seqan/sequence.h>
#include "pat.h"
#include "qual.h"
#include "ebwt_search_util.h"
#include "range.h"
#include "range_source.h"
#include "aligner_metrics.h"
#include "search_globals.h"

/**
 * Class that coordinates quality- and quantity-aware backtracking over
 * some range of a read sequence.
 *
 * The creator can configure the BacktrackManager to treat different
 * stretches of the read differently.
 */
class GreedyDFSRangeSource {

	typedef std::pair<int, int> TIntPair;
	typedef seqan::String<seqan::Dna> DnaString;

public:
	GreedyDFSRangeSource(
			const Ebwt<DnaString>* ebwt,
			const EbwtSearchParams<DnaString>& params,
			const BitPairReference* refs,
			uint32_t qualThresh,  /// max acceptable q-distance
			const int maxBts, /// maximum # backtracks allowed
			uint32_t reportPartials = 0,
			bool reportExacts = true,
			bool reportRanges = false,
			PartialAlignmentManager* partials = NULL,
			String<QueryMutation>* muts = NULL,
			bool verbose = true,
			vector<String<Dna5> >* os = NULL,
			bool considerQuals = true,  // whether to consider quality values when making backtracking decisions
			bool halfAndHalf = false, // hacky way of supporting separate revisitable regions
			bool maqPenalty = true) :
		_refs(refs),
		_qry(NULL),
		_qlen(0),
		_qual(NULL),
		_name(NULL),
		_ebwt(ebwt),
		_params(params),
		_unrevOff(0),
		_1revOff(0),
		_2revOff(0),
		_3revOff(0),
		_maqPenalty(maqPenalty),
		_qualThresh(qualThresh),
		_pairs(NULL),
		_elims(NULL),
		_mms(),
		_refcs(),
		_chars(NULL),
		_reportPartials(reportPartials),
		_reportExacts(reportExacts),
		_reportRanges(reportRanges),
		_partials(partials),
		_muts(muts),
		_os(os),
		_sanity(_os != NULL && _os->size() > 0),
		_considerQuals(considerQuals),
		_halfAndHalf(halfAndHalf),
		_5depth(0),
		_3depth(0),
		_numBts(0),
		_totNumBts(0),
		_maxBts(maxBts),
		_precalcedSideLocus(false),
		_preLtop(),
		_preLbot(),
		_verbose(verbose),
		_ihits(0llu)
	{ }

	~GreedyDFSRangeSource() {
		if(_pairs != NULL) delete[] _pairs;
		if(_elims != NULL) delete[] _elims;
		if(_chars != NULL) delete[] _chars;
	}

	/**
	 * Set a new query read.
	 */
	void setQuery(Read& r) {
		const bool fw = _params.fw();
		const bool ebwtFw = _ebwt->fw();
		if(ebwtFw) {
			_qry  = fw ? &r.patFw : &r.patRc;
			_qual = fw ? &r.qual  : &r.qualRev;
		} else {
			_qry  = fw ? &r.patFwRev : &r.patRcRev;
			_qual = fw ? &r.qualRev  : &r.qual;
		}
		_name = &r.name;
		// Reset _qlen
		if(length(*_qry) > _qlen) {
			try {
				_qlen = length(*_qry);
				// Resize _pairs
				if(_pairs != NULL) { delete[] _pairs; }
				_pairs = new TIndexOffU[_qlen*_qlen*8];
				// Resize _elims
				if(_elims != NULL) { delete[] _elims; }
				_elims = new uint8_t[_qlen*_qlen];
				memset(_elims, 0, _qlen*_qlen);
				// Resize _chars
				if(_chars != NULL) { delete[] _chars; }
				_chars = new char[_qlen];
				assert(_pairs != NULL && _elims != NULL && _chars != NULL);
			} catch(std::bad_alloc& e) {
				ThreadSafe _ts(&gLock);
				cerr << "Unable to allocate memory for depth-first "
				     << "backtracking search; new length = " << length(*_qry)
				     << endl;
				throw 1;
			}
		} else {
			// New length is less than old length, so there's no need
			// to resize any data structures.
			assert(_pairs != NULL && _elims != NULL && _chars != NULL);
			_qlen = length(*_qry);
		}
		_mms.clear();
		_refcs.clear();
		assert_geq(length(*_qual), _qlen);
		if(_verbose) {
			cout << "setQuery(_qry=" << (*_qry) << ", _qual=" << (*_qual) << ")" << endl;
		}
		// Initialize the random source using new read as part of the
		// seed.
		_color = r.color;
		_seed = r.seed;
		_patid = r.patid;
		_primer = r.primer;
		_trimc = r.trimc;
		_rand.init(r.seed);
	}

	/**
	 * Apply a batch of mutations to this read, possibly displacing a
	 * previous batch of mutations.
	 */
	void setMuts(String<QueryMutation>* muts) {
		if(_muts != NULL) {
			// Undo previous mutations
			assert_gt(length(*_muts), 0);
			undoPartialMutations();
		}
		_muts = muts;
		if(_muts != NULL) {
			assert_gt(length(*_muts), 0);
			applyPartialMutations();
		}
	}

	/**
	 * Set backtracking constraints.
	 */
	void setOffs(uint32_t depth5,   // depth of far edge of hi-half
	             uint32_t depth3,   // depth of far edge of lo-half
	             uint32_t unrevOff, // depth above which we cannot backtrack
	             uint32_t revOff1,  // depth above which we may backtrack just once
	             uint32_t revOff2,  // depth above which we may backtrack just twice
	             uint32_t revOff3)  // depth above which we may backtrack just three times
	{
		_5depth   = depth5;
		_3depth   = depth3;
		assert_geq(depth3, depth5);
		_unrevOff = unrevOff;
		_1revOff  = revOff1;
		_2revOff  = revOff2;
		_3revOff  = revOff3;
	}

	/**
	 * Reset number of backtracks to 0.
	 */
	void resetNumBacktracks() {
		_totNumBts = 0;
	}

	/**
	 * Return number of backtracks since the last time the count was
	 * reset.
	 */
	uint32_t numBacktracks() {
		return _totNumBts;
	}

	/**
	 * Set whether to report exact hits.
	 */
	void setReportExacts(int stratum) {
		_reportExacts = stratum;
	}

	/**
	 * Set the Bowtie index to search against.
	 */
	void setEbwt(const Ebwt<String<Dna> >* ebwt) {
		_ebwt = ebwt;
	}

	/**
	 * Return the current range
	 */
	Range& range() {
		return _curRange;
	}

	/**
	 * Set _qlen.  Don't let it exceed length of query.
	 */
	void setQlen(uint32_t qlen) {
		assert(_qry != NULL);
		_qlen = min<uint32_t>((uint32_t)length(*_qry), qlen);
	}

	/// Return the maximum number of allowed backtracks in a given call
	/// to backtrack()
	uint32_t maxBacktracks() {
		return _maxBts;
	}

	/**
	 * Initiate the recursive backtracking routine starting at the
	 * extreme right-hand side of the pattern.  Use the ftab to match
	 * the first several characters in one chomp, as long as doing so
	 * does not "jump over" any legal backtracking targets.
	 *
	 * Return true iff the HitSink has indicated that we're done with
	 * this read.
	 */
	bool backtrack(uint32_t ham = 0) {
		assert_gt(length(*_qry), 0);
		assert_leq(_qlen, length(*_qry));
		assert_geq(length(*_qual), length(*_qry));
		const Ebwt<String<Dna> >& ebwt = *_ebwt;
		int ftabChars = ebwt._eh._ftabChars;
		int nsInSeed = 0; int nsInFtab = 0;
		if(!tallyNs(nsInSeed, nsInFtab)) {
			// No alignments are possible because of the distribution
			// of Ns in the read in combination with the backtracking
			// constraints.
			return false;
		}
		bool ret;
		// m = depth beyond which ftab must not extend or else we might
		// miss some legitimate paths
		uint32_t m = min<uint32_t>(_unrevOff, (uint32_t)_qlen);
		if(nsInFtab == 0 && m >= (uint32_t)ftabChars) {
			uint32_t ftabOff = calcFtabOff();
			TIndexOffU top = ebwt.ftabHi(ftabOff);
			TIndexOffU bot = ebwt.ftabLo(ftabOff+1);
			if(_qlen == (TIndexOffU)ftabChars && bot > top) {
				// We have a match!
				if(_reportPartials > 0) {
					// Oops - we're trying to find seedlings, so we've
					// gone too far; start again
					ret = backtrack(0,   // depth
					                0,   // top
					                0,   // bot
					                ham,
					                nsInFtab > 0);
				} else {
					// We have a match!
					ret = reportAlignment(0, top, bot, ham);
				}
			} else if (bot > top) {
				// We have an arrow pair from which we can backtrack
				ret = backtrack(ftabChars, // depth
				                top,       // top
				                bot,       // bot
				                ham,
				                nsInFtab > 0);
			} else {
				// The arrows are already closed; give up
				ret = false;
			}
		} else {
			// The ftab *does* extend past the unrevisitable portion;
			// we can't use it in this case, because we might jump past
			// a legitimate mismatch
			ret = backtrack(0,   // depth
			                0,   // top
			                0,   // bot
			                ham,
			                // disable ftab jumping if there is more
			                // than 1 N in it
			                nsInFtab > 0);
		}
		if(finalize()) ret = true;
		return ret;
	}

	/**
	 * If there are any buffered results that have yet to be committed,
	 * commit them.  This happens when looking for partial alignments.
	 */
	bool finalize() {
		bool ret = false;
		if(_reportPartials > 0) {
			// We're in partial alignment mode; take elements of the
			// _partialBuf and install them in the _partials database
			assert(_partials != NULL);
			if(_partialsBuf.size() > 0) {
#ifndef NDEBUG
				for(size_t i = 0; i < _partialsBuf.size(); i++) {
					assert(_partialsBuf[i].repOk(_qualThresh, (uint32_t)_qlen, (*_qual), _maqPenalty));
				}
#endif
				_partials->addPartials(_params.patId(), _partialsBuf);
				_partialsBuf.clear();
				ret = true;
			} else {
				assert(!ret);
			}
		}
		assert_eq(0, _partialsBuf.size());
		return ret;
	}

	/**
	 * Starting at the given "depth" relative to the 5' end, and the
	 * given top and bot indexes (where top=0 and bot=0 means it's up
	 * to us to calculate the initial range), and initial weighted
	 * hamming distance iham, find a hit using randomized, quality-
	 * aware backtracking.
	 */
	bool backtrack(uint32_t depth,
	               TIndexOffU top,
	               TIndexOffU bot,
	               uint32_t iham = 0,
	               bool disableFtab = false)
	{
		HitSinkPerThread& sink = _params.sink();
		_ihits = sink.retainedHits().size();

		// Initiate the recursive, randomized quality-aware backtracker
		// with a stack depth of 0 (no backtracks so far)
		_bailedOnBacktracks = false;
		bool done = backtrack(0, depth, _unrevOff, _1revOff, _2revOff, _3revOff,
		                      top, bot, iham, iham, _pairs, _elims, disableFtab);

		_totNumBts += _numBts;
		_numBts = 0;
		_precalcedSideLocus = false;
		_bailedOnBacktracks = false;
		return done;
	}

	/**
	 * Recursive routine for progressing to the next backtracking
	 * decision given some initial conditions.  If a hit is found, it
	 * is recorded and true is returned.  Otherwise, if there are more
	 * backtracking opportunities, the function will call itself
	 * recursively and return the result.  As soon as there is a
	 * mismatch and no backtracking opportunities, false is returned.
	 */
	bool backtrack(uint32_t  stackDepth, // depth of the recursion stack; = # mismatches so far
	               uint32_t  depth,    // next depth where a post-pair needs to be calculated
	               uint32_t  unrevOff, // depths < unrevOff are unrevisitable
	               uint32_t  oneRevOff,// depths < oneRevOff are 1-revisitable
	               uint32_t  twoRevOff,// depths < twoRevOff are 2-revisitable
	               uint32_t  threeRevOff,// depths < threeRevOff are 3-revisitable
	               TIndexOffU  top,      // top arrow in pair prior to 'depth'
	               TIndexOffU  bot,      // bottom arrow in pair prior to 'depth'
	               uint32_t  ham,      // weighted hamming distance so far
	               uint32_t  iham,     // initial weighted hamming distance
	               TIndexOffU* pairs,    // portion of pairs array to be used for this backtrack frame
	               uint8_t*  elims,    // portion of elims array to be used for this backtrack frame
	               bool disableFtab = false)
	{
		// Can't have already exceeded weighted hamming distance threshold
		assert_leq(stackDepth, depth);
		assert_gt(length(*_qry), 0);
		assert_leq(_qlen, length(*_qry));
		assert_geq(length(*_qual), length(*_qry));
		assert(_qry != NULL);
		assert(_qual != NULL);
		assert(_name != NULL);
		assert(_qlen != 0);
		assert_leq(ham, _qualThresh);
		assert_lt(depth, _qlen); // can't have run off the end of qry
		assert_geq(bot, top);    // could be that both are 0
		assert(pairs != NULL);
		assert(elims != NULL);
		assert_leq(stackDepth, _qlen);
		const Ebwt<String<Dna> >& ebwt = *_ebwt;
		HitSinkPerThread& sink = _params.sink();
		uint64_t prehits = sink.numValidHits();
		if(_halfAndHalf) {
			assert_eq(0, _reportPartials);
			assert_gt(_3depth, _5depth);
		}
		if(_reportPartials) {
			assert(!_halfAndHalf);
		}
		if(_verbose) {
			cout << "  backtrack(stackDepth=" << stackDepth << ", "
			     << "depth=" << depth << ", "
			     << "top=" << top << ", "
			     << "bot=" << bot << ", "
			     << "ham=" << ham << ", "
			     << "iham=" << iham << ", "
			     << "pairs=" << pairs << ", "
			     << "elims=" << (void*)elims << "): \"";
			for(int i = (int)depth - 1; i >= 0; i--) {
				cout << _chars[i];
			}
			cout << "\"" << endl;
		}
		// Do this early on so that we can clear _precalcedSideLocus
		// before we have too many opportunities to bail and leave it
		// 'true'
		SideLocus ltop, lbot;
		if(_precalcedSideLocus) {
			ltop = _preLtop;
			lbot = _preLbot;
			_precalcedSideLocus = false;
		} else if(top != 0 || bot != 0) {
			SideLocus::initFromTopBot(top, bot, ebwt._eh, ebwt._ebwt, ltop, lbot);
		}
		// Check whether we've exceeded any backtracking limit
		if(_halfAndHalf) {
			if(_maxBts > 0 && _numBts == _maxBts) {
				_bailedOnBacktracks = true;
				return false;
			}
			_numBts++;
		}
		// # positions with at least one legal outgoing path
		uint32_t altNum = 0;
		// # positions tied for "best" outgoing qual
		uint32_t eligibleNum = 0;
		// total range-size for all eligibles
		TIndexOffU eligibleSz = 0;
		// If there is just one eligible slot at the moment (a common
		// case), these are its parameters
		uint32_t eli = 0;
		bool     elignore = true; // ignore the el values because they didn't come from a recent override
		TIndexOffU eltop = 0;
		TIndexOffU elbot = 0;
		uint32_t elham = ham;
		char     elchar = 0;
		int      elcint = 0;
		// The lowest quality value associated with any alternative
		// ranges; all alternative ranges with this quality are
		// eligible
		uint8_t lowAltQual = 0xff;
		uint32_t d = depth;
		uint32_t cur = (uint32_t)_qlen - d - 1; // current offset into _qry
		while(cur < _qlen) {
			// Try to advance further given that
			if(_verbose) {
				cout << "    cur=" << cur << " \"";
				for(int i = (int)d - 1; i >= 0; i--) {
					cout << _chars[i];
				}
				cout << "\"";
			}

			// If we're searching for a half-and-half solution, then
			// enforce the boundary-crossing constraints here.
			if(_halfAndHalf && !hhCheckTop(stackDepth, d, iham, _mms, prehits)) {
				return false;
			}

			bool curIsEligible = false;
			// Reset eligibleNum and eligibleSz if there are any
			// eligible pairs discovered at this spot
			bool curOverridesEligible = false;
			// Determine whether ranges at this location are
			// candidates for backtracking
			int c = (int)(*_qry)[cur];
			assert_leq(c, 4);
			uint8_t q = qualAt(cur);
			// The current query position is a legit alternative if it a) is
			// not in the unrevisitable region, and b) the quality ceiling (if
			// one exists) is not exceeded
			bool curIsAlternative =
				(d >= unrevOff) &&
			    (!_considerQuals ||
			     (ham + mmPenalty(_maqPenalty, q) <= _qualThresh));
			if(curIsAlternative) {
				if(_considerQuals) {
					// Is it the best alternative?
					if(q < lowAltQual) {
						// Ranges at this depth in this backtracking frame are
						// eligible, unless we learn otherwise.  Ranges previously
						// thought to be eligible are not any longer.
						curIsEligible = true;
						curOverridesEligible = true;
					} else if(q == lowAltQual) {
						// Ranges at this depth in this backtracking frame
						// are eligible, unless we learn otherwise
						curIsEligible = true;
					}
				} else {
					// When quality values are not considered, all positions
					// are eligible
					curIsEligible = true;
				}
			}
			if(curIsEligible) assert(curIsAlternative);
			if(curOverridesEligible) assert(curIsEligible);
			if(curIsAlternative && !curIsEligible) {
				assert_gt(eligibleSz, 0);
				assert_gt(eligibleNum, 0);
			}
			if(_verbose) {
				cout << " alternative: " << curIsAlternative;
				cout << ", eligible: " << curIsEligible;
				if(curOverridesEligible) cout << "(overrides)";
				cout << endl;
			}

			// If c is 'N', then it's guaranteed to be a mismatch
			if(c == 4 && d > 0) {
				// Force the 'else if(curIsAlternative)' branch below
				top = bot = 1;
			} else if(c == 4) {
				// We'll take the 'if(top == 0 && bot == 0)' branch below
				assert_eq(0, top);
				assert_eq(0, bot);
			}
			// Calculate the ranges for this position
			if(top == 0 && bot == 0) {
				// Calculate first quartet of ranges using the _fchr[]
				// array
				               pairs[0 + 0] = ebwt._fchr[0];
				pairs[0 + 4] = pairs[1 + 0] = ebwt._fchr[1];
				pairs[1 + 4] = pairs[2 + 0] = ebwt._fchr[2];
				pairs[2 + 4] = pairs[3 + 0] = ebwt._fchr[3];
				pairs[3 + 4]                = ebwt._fchr[4];
				// Update top and bot
				if(c < 4) {
					top = pairTop(pairs, d, c); bot = pairBot(pairs, d, c);
					assert_geq(bot, top);
				}
			} else if(curIsAlternative) {
				// Clear pairs
				memset(&pairs[d*8], 0, 8 * OFF_SIZE);
				// Calculate next quartet of ranges
				ebwt.mapLFEx(ltop, lbot, &pairs[d*8], &pairs[(d*8)+4]);
				// Update top and bot
				if(c < 4) {
					top = pairTop(pairs, d, c); bot = pairBot(pairs, d, c);
					assert_geq(bot, top);
				}
			} else {
				// This query character is not even a legitimate
				// alternative (because backtracking here would blow
				// our mismatch quality budget), so no need to do the
				// bookkeeping for the entire quartet, just do c
				if(c < 4) {
					if(top+1 == bot) {
						bot = top = ebwt.mapLF1(top, ltop, c);
						if(bot != OFF_MASK) bot++;
					} else {
						top = ebwt.mapLF(ltop, c); bot = ebwt.mapLF(lbot, c);
						assert_geq(bot, top);
					}
				}
			}
			if(top != bot) {
				// Calculate loci from row indices; do it now so that
				// those prefetches are fired off as soon as possible.
				// This eventually calls SideLocus.initfromRow().
				SideLocus::initFromTopBot(top, bot, ebwt._eh, ebwt._ebwt, ltop, lbot);
			}
			// Update the elim array
			eliminate(elims, d, c);

			if(curIsAlternative) {
				// Given the just-calculated range quartet, update
				// elims, altNum, eligibleNum, eligibleSz
				for(int i = 0; i < 4; i++) {
					if(i == c) continue;
					assert_leq(pairTop(pairs, d, i), pairBot(pairs, d, i));
					TIndexOffU spread = pairSpread(pairs, d, i);
					if(spread == 0) {
						// Indicate this char at this position is
						// eliminated as far as this backtracking frame is
						// concerned, since its range is empty
						elims[d] |= (1 << i);
						assert_lt(elims[d], 16);
					}
					if(spread > 0 && ((elims[d] & (1 << i)) == 0)) {
						// This char at this position is an alternative
						if(curIsEligible) {
							if(curOverridesEligible) {
								// Only now that we know there is at least
								// one potential backtrack target at this
								// most-eligible position should we reset
								// these eligibility parameters
								lowAltQual = q;
								eligibleNum = 0;
								eligibleSz = 0;
								curOverridesEligible = false;
								// Remember these parameters in case
								// this turns out to be the only
								// eligible target
								eli = d;
								eltop = pairTop(pairs, d, i);
								elbot = pairBot(pairs, d, i);
								assert_eq(elbot-eltop, spread);
								elham = mmPenalty(_maqPenalty, q);
								elchar = "acgt"[i];
								elcint = i;
								elignore = false;
							}
							eligibleSz += spread;
							eligibleNum++;
						}
						assert_gt(eligibleSz, 0);
						assert_gt(eligibleNum, 0);
						altNum++;
					}
				}
			}
			if(altNum > 0) {
				assert_gt(eligibleSz, 0);
				assert_gt(eligibleNum, 0);
			}
			assert_leq(eligibleNum, eligibleSz);
			assert_leq(eligibleNum, altNum);
			assert_lt(elims[d], 16);
			assert(sanityCheckEligibility(depth, d, unrevOff, lowAltQual, eligibleSz, eligibleNum, pairs, elims));

			// Achieved a match, but need to keep going
			bool backtrackDespiteMatch = false;
			bool reportedPartial = false;
			if(cur == 0 &&  // we've consumed the entire pattern
			   top < bot && // there's a hit to report
			   stackDepth < _reportPartials && // not yet used up our mismatches
			   _reportPartials > 0)  // there are still legel backtracking targets
			{
				assert(!_halfAndHalf);
				if(altNum > 0) backtrackDespiteMatch = true;
				if(stackDepth > 0) {
					// This is a legit seedling; report it
					reportPartial(stackDepth);
					reportedPartial = true;
				}
				// Now continue on to find legitimate seedlings with
				// more mismatches than this one
			}
			// Check whether we've obtained an exact alignment when
			// we've been instructed not to report exact alignments
			bool invalidExact = false;
			if(cur == 0 && stackDepth == 0 && bot > top && !_reportExacts) {
				invalidExact = true;
				backtrackDespiteMatch = true;
			}
			// Set this to true if the only way to make legal progress
			// is via one or more additional backtracks.  This is
			// helpful in half-and-half mode.
			bool mustBacktrack = false;
			bool invalidHalfAndHalf = false;
			if(_halfAndHalf) {
				ASSERT_ONLY(uint32_t lim = (_3revOff == _2revOff)? 2 : 3);
				if((d == (_5depth-1)) && top < bot) {
					// We're crossing the boundary separating the hi-half
					// from the non-seed portion of the read.
					// We should induce a mismatch if we haven't mismatched
					// yet, so that we don't waste time pursuing a match
					// that was covered by a previous phase
					assert_eq(0, _reportPartials);
					assert_leq(stackDepth, lim-1);
					invalidHalfAndHalf = (stackDepth == 0);
					if(stackDepth == 0 && altNum > 0) {
						backtrackDespiteMatch = true;
						mustBacktrack = true;
					} else if(stackDepth == 0) {
						// We're returning from the bottommost frame
						// without having found any hits; let's
						// sanity-check that there really aren't any
						return false;
					}
				}
				else if((d == (_3depth-1)) && top < bot) {
					// We're crossing the boundary separating the lo-half
					// from the non-seed portion of the read
					assert_eq(0, _reportPartials);
					assert_leq(stackDepth, lim);
					assert_gt(stackDepth, 0);
					// Count the mismatches in the lo and hi halves
					uint32_t loHalfMms = 0, hiHalfMms = 0;
					assert_geq(_mms.size(), stackDepth);
					for(size_t i = 0; i < stackDepth; i++) {
						uint32_t d = (uint32_t)_qlen - _mms[i] - 1;
						if     (d < _5depth) hiHalfMms++;
						else if(d < _3depth) loHalfMms++;
						else assert(false);
					}
					assert_leq(loHalfMms + hiHalfMms, lim);
					invalidHalfAndHalf = (loHalfMms == 0 || hiHalfMms == 0);
					if((stackDepth < 2 || invalidHalfAndHalf) && altNum > 0) {
						// We backtracked fewer times than necessary;
						// force a backtrack
						mustBacktrack = true;
						backtrackDespiteMatch = true;
					} else if(stackDepth < 2) {
						return false;
					}
				}
				if(d < _5depth-1) {
					assert_leq(stackDepth, lim-1);
				}
				else if(d >= _5depth && d < _3depth-1) {
					assert_gt(stackDepth, 0);
					assert_leq(stackDepth, lim);
				}
			}
			// This is necessary for the rare case where we're about
			// to declare success because bot > top and we've consumed
			// the final character, but all hits between top and bot
			// are spurious.  This check ensures that we keep looking
			// for non-spurious hits in that case.
			if(cur == 0 &&            // we made it to the left-hand-side of the read
			   bot > top &&           // there are alignments to report
			   !invalidHalfAndHalf && // alignment isn't disqualified by half-and-half requirement
			   !invalidExact &&       // alignment isn't disqualified by no-exact-hits setting
			   !reportedPartial)      // for when it's a partial alignment we've already reported
			{
				bool ret = reportAlignment(stackDepth, top, bot, ham);
				if(!ret) {
					// reportAlignment returned false, so enter the
					// backtrack loop and keep going
					top = bot;
				} else {
					// reportAlignment returned true, so stop
					return true;
				}
			}
			//
			// Mismatch with alternatives
			//
			while((top == bot || backtrackDespiteMatch) && altNum > 0) {
				if(_verbose) cout << "    top (" << top << "), bot ("
				                 << bot << ") with " << altNum
				                 << " alternatives, eligible: "
				                 << eligibleNum << ", " << eligibleSz
				                 << endl;
				assert_gt(eligibleSz, 0);
				assert_gt(eligibleNum, 0);
				// Mismatch!  Must now choose where we are going to
				// take our quality penalty.  We can only look as far
				// back as our last decision point.
				assert(sanityCheckEligibility(depth, d, unrevOff, lowAltQual, eligibleSz, eligibleNum, pairs, elims));
				// Pick out the arrow pair we selected and target it
				// for backtracking
				ASSERT_ONLY(uint32_t eligiblesVisited = 0);
				size_t i = d, j = 0;
				assert_geq(i, depth);
				TIndexOffU bttop = 0;
				TIndexOffU btbot = 0;
				uint32_t btham = ham;
				char     btchar = 0;
				int      btcint = 0;
				uint32_t icur = 0;
				// The common case is that eligibleSz == 1
				if(eligibleNum > 1 || elignore) {
					ASSERT_ONLY(bool foundTarget = false);
					// Walk from left to right
					for(; i >= depth; i--) {
						assert_geq(i, unrevOff);
						icur = (uint32_t)(_qlen - i - 1); // current offset into _qry
						uint8_t qi = qualAt(icur);
						assert_lt(elims[i], 16);
						if((qi == lowAltQual || !_considerQuals) && elims[i] != 15) {
							// This is the leftmost eligible position with at
							// least one remaining backtrack target
							TIndexOffU posSz = 0;
							// Add up the spreads for A, C, G, T
							for(j = 0; j < 4; j++) {
								if((elims[i] & (1 << j)) == 0) {
									assert_gt(pairSpread(pairs, i, j), 0);
									posSz += pairSpread(pairs, i, j);
								}
							}
							// Generate a random number
							assert_gt(posSz, 0);
							uint32_t r = _rand.nextU32() % posSz;
							for(j = 0; j < 4; j++) {
								if((elims[i] & (1 << j)) == 0) {
									// This range has not been eliminated
									ASSERT_ONLY(eligiblesVisited++);
									uint32_t spread = pairSpread(pairs, i, j);
									if(r < spread) {
										// This is our randomly-selected
										// backtrack target
										ASSERT_ONLY(foundTarget = true);
										bttop = pairTop(pairs, i, j);
										btbot = pairBot(pairs, i, j);
										btham += mmPenalty(_maqPenalty, qi);
										btcint = (uint32_t)j;
										btchar = "acgt"[j];
										assert_leq(btham, _qualThresh);
										break; // found our target; we can stop
									}
									r -= spread;
								}
							}
							assert(foundTarget);
							break; // escape left-to-right walk
						}
					}
					assert_leq(i, d);
					assert_lt(j, 4);
					assert_leq(eligiblesVisited, eligibleNum);
					assert(foundTarget);
					assert_neq(0, btchar);
					assert_gt(btbot, bttop);
					assert_leq(btbot-bttop, eligibleSz);
				} else {
					// There was only one eligible target; we can just
					// copy its parameters
					assert_eq(1, eligibleNum);
					assert(!elignore);
					i = eli;
					bttop = eltop;
					btbot = elbot;
					btham += elham;
					j = btcint = elcint;
					btchar = elchar;
					assert_neq(0, btchar);
					assert_gt(btbot, bttop);
					assert_leq(btbot-bttop, eligibleSz);
				}
				// This is the earliest that we know what the next top/
				// bot combo is going to be
				SideLocus::initFromTopBot(bttop, btbot,
				                          ebwt._eh, ebwt._ebwt,
				                          _preLtop, _preLbot);
				icur = (uint32_t)(_qlen - i - 1); // current offset into _qry
				// Slide over to the next backtacking frame within
				// pairs and elims; won't interfere with our frame or
				// any of our parents' frames
				TIndexOffU *newPairs = pairs + (_qlen*8);
				uint8_t  *newElims = elims + (_qlen);
				// If we've selected a backtracking target that's in
				// the 1-revisitable region, then we ask the recursive
				// callee to consider the 1-revisitable region as also
				// being unrevisitable (since we just "used up" all of
				// our visits)
				uint32_t btUnrevOff  = unrevOff;
				uint32_t btOneRevOff = oneRevOff;
				uint32_t btTwoRevOff = twoRevOff;
				uint32_t btThreeRevOff = threeRevOff;
				assert_geq(i, unrevOff);
				assert_geq(oneRevOff, unrevOff);
				assert_geq(twoRevOff, oneRevOff);
				assert_geq(threeRevOff, twoRevOff);
				if(i < oneRevOff) {
					// Extend unrevisitable region to include former 1-
					// revisitable region
					btUnrevOff = oneRevOff;
					// Extend 1-revisitable region to include former 2-
					// revisitable region
					btOneRevOff = twoRevOff;
					// Extend 2-revisitable region to include former 3-
					// revisitable region
					btTwoRevOff = threeRevOff;
				}
				else if(i < twoRevOff) {
					// Extend 1-revisitable region to include former 2-
					// revisitable region
					btOneRevOff = twoRevOff;
					// Extend 2-revisitable region to include former 3-
					// revisitable region
					btTwoRevOff = threeRevOff;
				}
				else if(i < threeRevOff) {
					// Extend 2-revisitable region to include former 3-
					// revisitable region
					btTwoRevOff = threeRevOff;
				}
				// Note the character that we're backtracking on in the
				// mm array:
				if(_mms.size() <= stackDepth) {
					assert_eq(_mms.size(), stackDepth);
					_mms.push_back(icur);
				} else {
					_mms[stackDepth] = icur;
				}
				assert_eq(1, dna4Cat[(int)btchar]);
				if(_refcs.size() <= stackDepth) {
					assert_eq(_refcs.size(), stackDepth);
					_refcs.push_back(btchar);
				} else {
					_refcs[stackDepth] = btchar;
				}
#ifndef NDEBUG
				for(uint32_t j = 0; j < stackDepth; j++) {
					assert_neq(_mms[j], icur);
				}
#endif
				_chars[i] = btchar;
				assert_leq(i+1, _qlen);
				bool ret;
				if(i+1 == _qlen) {
					ret = reportAlignment(stackDepth+1, bttop, btbot, btham);
				} else if(_halfAndHalf &&
				          !disableFtab &&
				          _2revOff == _3revOff &&
				          i+1 < (uint32_t)ebwt._eh._ftabChars &&
				          (uint32_t)ebwt._eh._ftabChars <= _5depth)
				{
					// The ftab doesn't extend past the unrevisitable portion,
					// so we can go ahead and use it
					// Rightmost char gets least significant bit-pairs
					int ftabChars = ebwt._eh._ftabChars;
					TIndexOffU ftabOff = (TIndexOffU)(int)(*_qry)[_qlen - ftabChars];
					assert_lt(ftabOff, 4);
					assert_lt(ftabOff, ebwt._eh._ftabLen-1);
					for(int j = ftabChars - 1; j > 0; j--) {
						ftabOff <<= 2;
						if(_qlen-j == icur) {
							ftabOff |= btcint;
						} else {
							assert_lt((int)(*_qry)[_qlen-j], 4);
							ftabOff |= (int)(*_qry)[_qlen-j];
						}
						assert_lt(ftabOff, ebwt._eh._ftabLen-1);
					}
					assert_lt(ftabOff, ebwt._eh._ftabLen-1);
					TIndexOffU ftabTop = ebwt.ftabHi(ftabOff);
					TIndexOffU ftabBot = ebwt.ftabLo(ftabOff+1);
					assert_geq(ftabBot, ftabTop);
					if(ftabTop == ftabBot) {
						ret = false;
					} else {
						assert(!_precalcedSideLocus);
						assert_leq(iham, _qualThresh);
						ret = backtrack(stackDepth+1,
						                ebwt._eh._ftabChars,
						                btUnrevOff,  // new unrevisitable boundary
						                btOneRevOff, // new 1-revisitable boundary
						                btTwoRevOff, // new 2-revisitable boundary
						                btThreeRevOff, // new 3-revisitable boundary
						                ftabTop,  // top arrow in range prior to 'depth'
						                ftabBot,  // bottom arrow in range prior to 'depth'
						                btham,  // weighted hamming distance so far
						                iham,   // initial weighted hamming distance
						                newPairs,
						                newElims);
					}
				} else {
					// We already called initFromTopBot for the range
					// we're going to continue from
					_precalcedSideLocus = true;
					assert_leq(iham, _qualThresh);
					// Continue from selected alternative range
					ret = backtrack(stackDepth+1,// added 1 mismatch to alignment
					                (uint32_t)i+1, // start from next position after
					                btUnrevOff,  // new unrevisitable boundary
					                btOneRevOff, // new 1-revisitable boundary
					                btTwoRevOff, // new 2-revisitable boundary
					                btThreeRevOff, // new 3-revisitable boundary
					                bttop,  // top arrow in range prior to 'depth'
					                btbot,  // bottom arrow in range prior to 'depth'
					                btham,  // weighted hamming distance so far
					                iham,   // initial weighted hamming distance
					                newPairs,
					                newElims);
				}
				if(ret) {
					assert_gt(sink.numValidHits(), prehits);
					return true; // return, signaling that we're done
				}
				if(_bailedOnBacktracks ||
				  (_halfAndHalf && (_maxBts > 0) && (_numBts >= _maxBts)))
				{
					_bailedOnBacktracks = true;
					return false;
				}
				// No hit was reported; update elims[], eligibleSz,
				// eligibleNum, altNum
				_chars[i] = (*_qry)[icur];
				assert_neq(15, elims[i]);
				ASSERT_ONLY(uint8_t oldElim = elims[i]);
				elims[i] |= (1 << j);
				assert_lt(elims[i], 16);
				assert_gt(elims[i], oldElim);
				eligibleSz -= (btbot-bttop);
				eligibleNum--;
				elignore = true;
				assert_geq(eligibleNum, 0);
				altNum--;
				assert_geq(altNum, 0);
				if(altNum == 0) {
					// No alternative backtracking points; all legal
					// backtracking targets have been exhausted
					assert_eq(0, altNum);
					assert_eq(0, eligibleSz);
					assert_eq(0, eligibleNum);
					return false;
				}
				else if(eligibleNum == 0 && _considerQuals) {
					// Find the next set of eligible backtrack points
					// by re-scanning this backtracking frame (from
					// 'depth' up to 'd')
					lowAltQual = 0xff;
					for(size_t k = d; k >= depth && k <= _qlen; k--) {
						size_t kcur = _qlen - k - 1; // current offset into _qry
						uint8_t kq = qualAt(kcur);
						if(k < unrevOff) break; // already visited all revisitable positions
						bool kCurIsAlternative = (ham + mmPenalty(_maqPenalty, kq) <= _qualThresh);
						bool kCurOverridesEligible = false;
						if(kCurIsAlternative) {
							if(kq < lowAltQual) {
								// This target is more eligible than
								// any targets that came before, so we
								// set it to supplant/override them
								kCurOverridesEligible = true;
							}
							if(kq <= lowAltQual) {
								// Position is eligible
								for(int l = 0; l < 4; l++) {
									if((elims[k] & (1 << l)) == 0) {
										// Not yet eliminated
										TIndexOffU spread = pairSpread(pairs, k, l);
										if(kCurOverridesEligible) {
											// Clear previous eligible results;
											// this one's better
											lowAltQual = kq;
											kCurOverridesEligible = false;
											// Keep these parameters in
											// case this target turns
											// out to be the only
											// eligible target and we
											// can avoid having to
											// recalculate them
											eligibleNum = 0;
											eligibleSz = 0;
											eli = (uint32_t)k;
											eltop = pairTop(pairs, k, l);
											elbot = pairBot(pairs, k, l);
											assert_eq(elbot-eltop, spread);
											elham = mmPenalty(_maqPenalty, kq);
											elchar = "acgt"[l];
											elcint = l;
											elignore = false;
										}
										eligibleNum++;
										assert_gt(spread, 0);
										eligibleSz += spread;
									}
								}
							}
						}
					}
				}
				assert_gt(eligibleNum, 0);
				assert_leq(eligibleNum, altNum);
				assert_gt(eligibleSz, 0);
				assert_geq(eligibleSz, eligibleNum);
				assert(sanityCheckEligibility(depth, d, unrevOff, lowAltQual, eligibleSz, eligibleNum, pairs, elims));
				// Try again
			} // while(top == bot && altNum > 0)
			if(mustBacktrack || invalidHalfAndHalf || invalidExact) {
				return false;
			}
			// Mismatch with no alternatives
			if(top == bot && altNum == 0) {
				assert_eq(0, altNum);
				assert_eq(0, eligibleSz);
				assert_eq(0, eligibleNum);
				return false;
			}
			// Match!
			_chars[d] = (*_qry)[cur];
			d++; cur--;
		} // while(cur < _qlen)
		assert_eq(0xffffffff, cur);
		assert_gt(bot, top);
		if(_reportPartials > 0) {
			// Stack depth should not exceed given hamming distance
			assert_leq(stackDepth, _reportPartials);
		}
		bool ret = false;
		if(stackDepth >= _reportPartials) {
			ret = reportAlignment(stackDepth, top, bot, ham);
		}
		return ret;
	}

	/**
	 * Pretty print a hit along with the backtracking constraints.
	 */
	void printHit(const Hit& h) {
		::printHit(*_os, h, *_qry, _qlen, _unrevOff, _1revOff, _2revOff, _3revOff, _ebwt->fw());
	}

	/**
	 * Return true iff we're enforcing a half-and-half constraint
	 * (forced edits in both seed halves).
	 */
	bool halfAndHalf() const {
		return _halfAndHalf;
	}

protected:

	/**
	 * Return true iff we're OK to continue after considering which
	 * half-seed boundary we're passing through, together with the
	 * number of mismatches accumulated so far.  Return false if we
	 * should stop because a half-and-half constraint is violated.  If
	 * we're not currently passing a half-seed boundary, just return
	 * true.
	 */
	bool hhCheck(uint32_t stackDepth, uint32_t depth,
	             const std::vector<uint32_t>& mms, bool empty)
	{
		ASSERT_ONLY(uint32_t lim = (_3revOff == _2revOff)? 2 : 3);
		if((depth == (_5depth-1)) && !empty) {
			// We're crossing the boundary separating the hi-half
			// from the non-seed portion of the read.
			// We should induce a mismatch if we haven't mismatched
			// yet, so that we don't waste time pursuing a match
			// that was covered by a previous phase
			assert_eq(0, _reportPartials);
			assert_leq(stackDepth, lim-1);
			return stackDepth > 0;
		} else if((depth == (_3depth-1)) && !empty) {
			// We're crossing the boundary separating the lo-half
			// from the non-seed portion of the read
			assert_eq(0, _reportPartials);
			assert_leq(stackDepth, lim);
			assert_gt(stackDepth, 0);
			// Count the mismatches in the lo and hi halves
			uint32_t loHalfMms = 0, hiHalfMms = 0;
			for(size_t i = 0; i < stackDepth; i++) {
				uint32_t depth = (uint32_t)(_qlen - mms[i] - 1);
				if     (depth < _5depth) hiHalfMms++;
				else if(depth < _3depth) loHalfMms++;
				else assert(false);
			}
			assert_leq(loHalfMms + hiHalfMms, lim);
			bool invalidHalfAndHalf = (loHalfMms == 0 || hiHalfMms == 0);
			return (stackDepth >= 2 && !invalidHalfAndHalf);
		}
		if(depth < _5depth-1) {
			assert_leq(stackDepth, lim-1);
		}
		else if(depth >= _5depth && depth < _3depth-1) {
			assert_gt(stackDepth, 0);
			assert_leq(stackDepth, lim);
		}
		return true;
	}

	/**
	 * Calculate the stratum of the partial (or full) alignment
	 * currently under consideration.  Stratum is equal to the number
	 * of mismatches in the seed portion of the alignment.
	 */
	int calcStratum(const std::vector<TIndexOffU>& mms, uint32_t stackDepth) {
		int stratum = 0;
		for(size_t i = 0; i < stackDepth; i++) {
			if(mms[i] >= (_qlen - _3revOff)) {
				// This mismatch falls within the seed; count it
				// toward the stratum to report
				stratum++;
				// Don't currently support more than 3
				// mismatches in the seed
				assert_leq(stratum, 3);
			}
		}
		return stratum;
	}

	/**
	 * Mark character c at depth d as being eliminated with respect to
	 * future backtracks.
	 */
	void eliminate(uint8_t *elims, uint32_t d, int c) {
		if(c < 4) {
			elims[d] = (1 << c);
			assert_gt(elims[d], 0);
			assert_lt(elims[d], 16);
		} else {
			elims[d] = 0;
		}
		assert_lt(elims[d], 16);
	}

	/**
	 * Return true iff the state of the backtracker as encoded by
	 * stackDepth, d and iham is compatible with the current half-and-
	 * half alignment mode.  prehits is for sanity checking when
	 * bailing.
	 */
	bool hhCheckTop(uint32_t stackDepth,
	                uint32_t d,
	                uint32_t iham,
	                const std::vector<TIndexOffU>& mms,
	                uint64_t prehits = 0xffffffffffffffffllu)
	{
		assert_eq(0, _reportPartials);
		// Crossing from the hi-half into the lo-half
		if(d == _5depth) {
			if(_3revOff == _2revOff) {
				// Total of 2 mismatches allowed: 1 hi, 1 lo
				// The backtracking logic should have prevented us from
				// backtracking more than once into this region
				assert_leq(stackDepth, 1);
				// Reject if we haven't encountered mismatch by this point
				if(stackDepth == 0) {
					return false;
				}
			} else { // if(_3revOff != _2revOff)
				// Total of 3 mismatches allowed: 1 hi, 1 or 2 lo
				// The backtracking logic should have prevented us from
				// backtracking more than twice into this region
				assert_leq(stackDepth, 2);
				// Reject if we haven't encountered mismatch by this point
				if(stackDepth < 1) {
					return false;
				}
			}
		} else if(d == _3depth) {
			// Crossing from lo-half to outside of the seed
			if(_3revOff == _2revOff) {
				// Total of 2 mismatches allowed: 1 hi, 1 lo
				// The backtracking logic should have prevented us from
				// backtracking more than twice within this region
				assert_leq(stackDepth, 2);
				// Must have encountered two mismatches by this point
				if(stackDepth < 2) {
					// We're returning from the bottommost frame
					// without having found any hits; let's
					// sanity-check that there really aren't any
					return false;
				}
			} else { // if(_3revOff != _2revOff)
				// Total of 3 mismatches allowed: 1 hi, 1 or 2 lo
				// Count the mismatches in the lo and hi halves
				int loHalfMms = 0, hiHalfMms = 0;
				assert_geq(mms.size(), stackDepth);
				for(size_t i = 0; i < stackDepth; i++) {
					TIndexOffU d = (TIndexOffU)(_qlen - mms[i] - 1);
					if     (d < _5depth) hiHalfMms++;
					else if(d < _3depth) loHalfMms++;
					else assert(false);
				}
				assert_leq(loHalfMms + hiHalfMms, 3);
				assert_gt(hiHalfMms, 0);
				if(loHalfMms == 0) {
					// We're returning from the bottommost frame
					// without having found any hits; let's
					// sanity-check that there really aren't any
					return false;
				}
				assert_geq(stackDepth, 2);
				// The backtracking logic should have prevented us from
				// backtracking more than twice within this region
				assert_leq(stackDepth, 3);
			}
		} else {
			// We didn't just cross a boundary, so do an in-between check
			if(d >= _5depth) {
				assert_geq(stackDepth, 1);
			} else if(d >= _3depth) {
				assert_geq(stackDepth, 2);
			}
		}
		return true;
	}

	/**
	 * Return the Phred quality value for the most likely base at
	 * offset 'off' in the read.
	 */
	inline uint8_t qualAt(size_t off) {
		return phredCharToPhredQual((*_qual)[off]);
	}

	/// Get the top offset for character c at depth d
	inline TIndexOffU pairTop(TIndexOffU* pairs, size_t d, size_t c) {
		return pairs[d*8 + c + 0];
	}

	/// Get the bot offset for character c at depth d
	inline TIndexOffU pairBot(TIndexOffU* pairs, size_t d, size_t c) {
		return pairs[d*8 + c + 4];
	}

	/// Get the spread between the bot and top offsets for character c
	/// at depth d
	inline TIndexOffU pairSpread(TIndexOffU* pairs, size_t d, size_t c) {
		assert_geq(pairBot(pairs, d, c), pairTop(pairs, d, c));
		return pairBot(pairs, d, c) - pairTop(pairs, d, c);
	}

	/**
	 * Tally how many Ns occur in the seed region and in the ftab-
	 * jumpable region of the read.  Check whether the mismatches
	 * induced by the Ns already violates the current policy.  Return
	 * false if the policy is already violated, true otherwise.
	 */
	bool tallyNs(int& nsInSeed, int& nsInFtab) {
		const Ebwt<String<Dna> >& ebwt = *_ebwt;
		int ftabChars = ebwt._eh._ftabChars;
		// Count Ns in the seed region of the read and short-circuit if
		// the configuration of Ns guarantees that there will be no
		// valid alignments given the backtracking constraints.
		for(size_t i = 0; i < _3revOff; i++) {
			if((int)(*_qry)[_qlen-i-1] == 4) {
				nsInSeed++;
				if(nsInSeed == 1) {
					if(i < _unrevOff) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else if(nsInSeed == 2) {
					if(i < _1revOff) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else if(nsInSeed == 3) {
					if(i < _2revOff) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else {
					assert_gt(nsInSeed, 3);
					return false;     // Exceeded mm budget on Ns alone
				}
			}
		}
		// Calculate the number of Ns there are in the region that
		// would get jumped over if the ftab were used.
		for(size_t i = 0; i < (size_t)ftabChars && i < _qlen; i++) {
			if((int)(*_qry)[_qlen-i-1] == 4) nsInFtab++;
		}
		return true;
	}

	/**
	 * Calculate the offset into the ftab for the rightmost 'ftabChars'
	 * characters of the current query. Rightmost char gets least
	 * significant bit-pair.
	 */
	uint32_t calcFtabOff() {
		const Ebwt<String<Dna> >& ebwt = *_ebwt;
		int ftabChars = ebwt._eh._ftabChars;
		uint32_t ftabOff = (*_qry)[_qlen - ftabChars];
		assert_lt(ftabOff, 4);
		assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		for(int i = ftabChars - 1; i > 0; i--) {
			ftabOff <<= 2;
			assert_lt((uint32_t)(*_qry)[_qlen-i], 4);
			ftabOff |= (uint32_t)(*_qry)[_qlen-i];
			assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		}
		assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		return ftabOff;
	}

	/**
	 * Mutate the _qry string according to the contents of the _muts
	 * array, which represents a partial alignment.
	 */
	void applyPartialMutations() {
		if(_muts == NULL) {
			// No mutations to apply
			return;
		}
		for(size_t i = 0; i < length(*_muts); i++) {
			const QueryMutation& m = (*_muts)[i];
			assert_lt(m.pos, _qlen);
			assert_leq(m.oldBase, 4);
			assert_lt(m.newBase, 4);
			assert_neq(m.oldBase, m.newBase);
			assert_eq((uint32_t)((*_qry)[m.pos]), (uint32_t)m.oldBase);
			(*_qry)[m.pos] = (Dna5)(int)m.newBase; // apply it
		}
	}

	/**
	 * Take partial-alignment mutations present in the _muts list and
	 * place them on the _mm list so that they become part of the
	 * reported alignment.
	 */
	void promotePartialMutations(int stackDepth) {
		if(_muts == NULL) {
			// No mutations to undo
			return;
		}
		size_t numMuts = length(*_muts);
		assert_leq(numMuts, _qlen);
		for(size_t i = 0; i < numMuts; i++) {
			// Entries in _mms[] are in terms of offset into
			// _qry - not in terms of offset from 3' or 5' end
			assert_lt(stackDepth + i, _qlen);
			// All partial-alignment mutations should fall
			// within bounds
			assert_lt((*_muts)[i].pos, _qlen);
			// All partial-alignment mutations should fall
			// within unrevisitable region
			assert_lt(_qlen - (*_muts)[i].pos - 1, _unrevOff);
#ifndef NDEBUG
			// Shouldn't be any overlap between mismatched positions
			// and positions that mismatched in the partial alignment.
			for(size_t j = 0; j < stackDepth + i; j++) {
				assert_neq(_mms[j], (uint32_t)(*_muts)[i].pos);
			}
#endif
			if(_mms.size() <= stackDepth + i) {
				assert_eq(_mms.size(), stackDepth + i);
				_mms.push_back((*_muts)[i].pos);
			} else {
				_mms[stackDepth + i] = (*_muts)[i].pos;
			}
			if(_refcs.size() <= stackDepth + i) {
				assert_eq(_refcs.size(), stackDepth + i);
				_refcs.push_back("ACGT"[(*_muts)[i].newBase]);
			} else {
				_refcs[stackDepth + i] = "ACGT"[(*_muts)[i].newBase];
			}
		}
	}

	/**
	 * Undo mutations to the _qry string, returning it to the original
	 * read.
	 */
	void undoPartialMutations() {
		if(_muts == NULL) {
			// No mutations to undo
			return;
		}
		for(size_t i = 0; i < length(*_muts); i++) {
			const QueryMutation& m = (*_muts)[i];
			assert_lt(m.pos, _qlen);
			assert_leq(m.oldBase, 4);
			assert_lt(m.newBase, 4);
			assert_neq(m.oldBase, m.newBase);
			assert_eq((uint32_t)((*_qry)[m.pos]), (uint32_t)m.newBase);
			(*_qry)[m.pos] = (Dna5)(int)m.oldBase; // undo it
		}
	}

	/**
	 * Report a range of alignments with # mismatches = stackDepth and
	 * with the mutations (also mismatches) contained in _muts.  The
	 * range is delimited by top and bot.  Returns true iff one or more
	 * full alignments were successfully reported and the caller can
	 * stop searching.
	 */
	bool reportAlignment(uint32_t stackDepth, TIndexOffU top,
			TIndexOffU bot, uint16_t cost)
	{
#ifndef NDEBUG
		// No two elements of _mms[] should be the same
		assert_geq(_mms.size(), stackDepth);
		for(size_t i = 0; i < stackDepth; i++) {
			for(size_t j = i+1; j < stackDepth; j++) {
				assert_neq(_mms[j], _mms[i]);
			}
			// All elements of _mms[] should fall within bounds
			assert_lt(_mms[i], _qlen);
		}
#endif
		if(_reportPartials) {
			assert_leq(stackDepth, _reportPartials);
			if(stackDepth > 0) {
				// Report this partial alignment.  A partial alignment
				// is defined purely by its mismatches; top and bot are
				// ignored.
				reportPartial(stackDepth);
			}
			return false; // keep going - we want to find all partial alignments
		}
		int stratum = 0;
		if(stackDepth > 0) {
			stratum = calcStratum(_mms, stackDepth);
		}
		assert_lt(stratum, 4);
		assert_geq(stratum, 0);
		bool hit;
		// If _muts != NULL then this alignment extends a partial
		// alignment, so we have to account for the differences present
		// in the partial.
		if(_muts != NULL) {
			// Undo partial-alignment mutations to get original _qry
			ASSERT_ONLY(String<Dna5> tmp = (*_qry));
			undoPartialMutations();
			assert_neq(tmp, (*_qry));
			// Add the partial-alignment mutations to the _mms[] array
			promotePartialMutations(stackDepth);
			// All muts are in the seed, so they count toward the stratum
			size_t numMuts = length(*_muts);
			stratum += numMuts;
			cost |= (stratum << 14);
			assert_geq(cost, (uint32_t)(stratum << 14));
			// Report the range of full alignments
			hit = reportFullAlignment((uint32_t)(stackDepth + numMuts), top, bot, stratum, cost);
			// Re-apply partial-alignment mutations
			applyPartialMutations();
			assert_eq(tmp, (*_qry));
		} else {
			// Report the range of full alignments
			cost |= (stratum << 14);
			assert_geq(cost, (uint32_t)(stratum << 14));
			hit = reportFullAlignment(stackDepth, top, bot, stratum, cost);
		}
		return hit;
	}

	/**
	 * Report a range of full alignments with # mismatches = stackDepth.
	 * The range is delimited by top and bot.  Returns true if one or
	 * more alignments were successfully reported.  Returns true iff
	 * one or more full alignments were successfully reported and the
	 * caller can stop searching.
	 */
	bool reportFullAlignment(uint32_t stackDepth,
			TIndexOffU top,
			TIndexOffU bot,
	                         int stratum,
	                         uint16_t cost)
	{
		assert_gt(bot, top);
		if(stackDepth == 0 && !_reportExacts) {
			// We are not reporting exact hits (usually because we've
			// already reported them as part of a previous invocation
			// of the backtracker)
			return false;
		}
		assert(!_reportRanges);
		TIndexOffU spread = bot - top;
		// Pick a random spot in the range to begin report
		TIndexOffU r = top + (_rand.nextU<TIndexOffU>() % spread);
		for(TIndexOffU i = 0; i < spread; i++) {
			TIndexOffU ri = r + i;
			if(ri >= bot) ri -= spread;
			// reportChaseOne takes the _mms[] list in terms of
			// their indices into the query string; not in terms
			// of their offset from the 3' or 5' end.
			assert_geq(cost, (uint32_t)(stratum << 14));
			if(_ebwt->reportChaseOne((*_qry), _qual, _name,
			                         _color, _primer, _trimc, colorExEnds,
			                         snpPhred, _refs, _mms, _refcs,
			                         stackDepth, ri, top, bot,
			                         (uint32_t)_qlen, stratum, cost, _patid,
			                         _seed, _params))
			{
				// Return value of true means that we can stop
				return true;
			}
			// Return value of false means that we should continue
			// searching.  This could happen if we the call to
			// reportChaseOne() reported a hit, but the user asked for
			// multiple hits and we haven't reached the ceiling yet.
			// This might also happen if the call to reportChaseOne()
			// didn't report a hit because the alignment was spurious
			// (i.e. overlapped some padding).
		}
		// All range elements were examined and we should keep going
		return false;
	}

	/**
	 * Report the partial alignment represented by the current stack
	 * state (_mms[] and stackDepth).
	 */
	bool reportPartial(uint32_t stackDepth) {
		// Sanity-check stack depth
		if(_3revOff != _2revOff) {
			assert_leq(stackDepth, 3);
		} else if(_2revOff != _1revOff) {
			assert_leq(stackDepth, 2);
		} else {
			assert_leq(stackDepth, 1);
		}

		// Possibly report
		assert_gt(_reportPartials, 0);
		assert(_partials != NULL);
		ASSERT_ONLY(uint32_t qualTot = 0);
		PartialAlignment al;
		al.u64.u64 = 0xffffffffffffffffllu;
		assert_leq(stackDepth, 3);
		assert_gt(stackDepth, 0);

		// First mismatch
		assert_gt(_mms.size(), 0);
		assert_lt(_mms[0], _qlen);
		// First, append the mismatch position in the read
		al.entry.pos0 = (uint16_t)_mms[0]; // pos
		ASSERT_ONLY(uint8_t qual0 = mmPenalty(_maqPenalty, phredCharToPhredQual((*_qual)[_mms[0]])));
		ASSERT_ONLY(qualTot += qual0);
		uint32_t ci = (uint32_t)(_qlen - _mms[0] - 1);
		// _chars[] is index in terms of RHS-relative depth
		int c = (int)(Dna5)_chars[ci];
		assert_lt(c, 4);
		assert_neq(c, (int)(*_qry)[_mms[0]]);
		// Second, append the substituted character for the position
		al.entry.char0 = c;

		if(stackDepth > 1) {
			assert_gt(_mms.size(), 1);
			// Second mismatch
			assert_lt(_mms[1], _qlen);
			// First, append the mismatch position in the read
			al.entry.pos1 = (uint16_t)_mms[1]; // pos
			ASSERT_ONLY(uint8_t qual1 = mmPenalty(_maqPenalty, phredCharToPhredQual((*_qual)[_mms[1]])));
			ASSERT_ONLY(qualTot += qual1);
			ci = (uint32_t)(_qlen - _mms[1] - 1);
			// _chars[] is index in terms of RHS-relative depth
			c = (int)(Dna5)_chars[ci];
			assert_lt(c, 4);
			assert_neq(c, (int)(*_qry)[_mms[1]]);
			// Second, append the substituted character for the position
			al.entry.char1 = c;
			if(stackDepth > 2) {
				assert_gt(_mms.size(), 2);
				// Second mismatch
				assert_lt(_mms[2], _qlen);
				// First, append the mismatch position in the read
				al.entry.pos2 = (uint16_t)_mms[2]; // pos
				ASSERT_ONLY(uint8_t qual2 = mmPenalty(_maqPenalty, phredCharToPhredQual((*_qual)[_mms[2]])));
				ASSERT_ONLY(qualTot += qual2);
				ci = (uint32_t)(_qlen - _mms[2] - 1);
				// _chars[] is index in terms of RHS-relative depth
				c = (int)(Dna5)_chars[ci];
				assert_lt(c, 4);
				assert_neq(c, (int)(*_qry)[_mms[2]]);
				// Second, append the substituted character for the position
				al.entry.char2 = c;
			} else {
				// Signal that the '2' slot is empty
				al.entry.pos2 = 0xffff;
			}
		} else {
			// Signal that the '1' slot is empty
			al.entry.pos1 = 0xffff;
		}

		assert_leq(qualTot, _qualThresh);
		assert(validPartialAlignment(al));
#ifndef NDEBUG
		assert(al.repOk(_qualThresh, (uint32_t)_qlen, (*_qual), _maqPenalty));
		for(size_t i = 0; i < _partialsBuf.size(); i++) {
			assert(validPartialAlignment(_partialsBuf[i]));
			assert(!samePartialAlignment(_partialsBuf[i], al));
		}
#endif
		_partialsBuf.push_back(al);
		return true;
	}

	/**
	 * Check that the given eligibility parameters (lowAltQual,
	 * eligibleSz, eligibleNum) are correct, given the appropriate
	 * inputs (pairs, elims, depth, d, unrevOff)
	 */
	bool sanityCheckEligibility(uint32_t  depth,
	                            uint32_t  d,
	                            uint32_t  unrevOff,
	                            uint32_t  lowAltQual,
	                            uint32_t  eligibleSz,
	                            uint32_t  eligibleNum,
	                            TIndexOffU* pairs,
	                            uint8_t*  elims)
	{
		// Sanity check that the lay of the land is as we
		// expect given eligibleNum and eligibleSz
		size_t i = max(depth, unrevOff), j = 0;
		uint32_t cumSz = 0;
		uint32_t eligiblesVisited = 0;
		for(; i <= d; i++) {
			uint32_t icur = (uint32_t)(_qlen - i - 1); // current offset into _qry
			uint8_t qi = qualAt(icur);
			assert_lt(elims[i], 16);
			if((qi == lowAltQual || !_considerQuals) && elims[i] != 15) {
				// This is an eligible position with at least
				// one remaining backtrack target
				for(j = 0; j < 4; j++) {
					if((elims[i] & (1 << j)) == 0) {
						// This pair has not been eliminated
						assert_gt(pairBot(pairs, i, j), pairTop(pairs, i, j));
						cumSz += pairSpread(pairs, i, j);
						eligiblesVisited++;
					}
				}
			}
		}
		assert_eq(cumSz, eligibleSz);
		assert_eq(eligiblesVisited, eligibleNum);
		return true;
	}

	const BitPairReference* _refs; // reference sequences (or NULL if not colorspace)
	String<Dna5>*       _qry;    // query (read) sequence
	size_t              _qlen;   // length of _qry
	String<char>*       _qual;   // quality values for _qry
	String<char>*       _name;   // name of _qry
	bool                _color;  // whether read is colorspace
	const Ebwt<String<Dna> >* _ebwt;   // Ebwt to search in
	const EbwtSearchParams<String<Dna> >& _params;   // Ebwt to search in
	uint32_t            _unrevOff; // unrevisitable chunk
	uint32_t            _1revOff;  // 1-revisitable chunk
	uint32_t            _2revOff;  // 2-revisitable chunk
	uint32_t            _3revOff;  // 3-revisitable chunk
	/// Whether to round qualities off Maq-style when calculating penalties
	bool                _maqPenalty;
	uint32_t            _qualThresh; // only accept hits with weighted
	                             // hamming distance <= _qualThresh
	TIndexOffU           *_pairs;  // ranges, leveled in parallel
	                             // with decision stack
	uint8_t            *_elims;  // which ranges have been
	                             // eliminated, leveled in parallel
	                             // with decision stack
	std::vector<TIndexOffU> _mms;  // array for holding mismatches
	std::vector<uint8_t> _refcs;  // array for holding mismatches
	// Entries in _mms[] are in terms of offset into
	// _qry - not in terms of offset from 3' or 5' end
	char               *_chars;  // characters selected so far
	// If > 0, report partial alignments up to this many mismatches
	uint32_t            _reportPartials;
	/// Do not report alignments with stratum < this limit
	bool                _reportExacts;
	/// When reporting a full alignment, report top/bot; don't chase
	/// any of the results
	bool                _reportRanges;
	/// Append partial alignments here
	PartialAlignmentManager *_partials;
	/// Set of mutations that apply for a partial alignment
	String<QueryMutation> *_muts;
	/// Reference texts (NULL if they are unavailable
	vector<String<Dna5> >* _os;
	/// Whether to use the _os array together with a naive matching
	/// algorithm to double-check reported alignments (or the lack
	/// thereof)
	bool                _sanity;
	/// Whether to consider quality values when deciding where to
	/// backtrack
	bool                _considerQuals;
	bool                _halfAndHalf;
	/// Depth of 5'-seed-half border
	uint32_t            _5depth;
	/// Depth of 3'-seed-half border
	uint32_t            _3depth;
	/// Default quals
	String<char>        _qualDefault;
	/// Number of backtracks in last call to backtrack()
	uint32_t            _numBts;
	/// Number of backtracks since last reset
	uint32_t            _totNumBts;
	/// Max # of backtracks to allow before giving up
	uint32_t            _maxBts;
	/// Whether we precalcualted the Ebwt locus information for the
	/// next top/bot pair
	bool    _precalcedSideLocus;
	/// Precalculated top locus
	SideLocus           _preLtop;
	/// Precalculated bot locus
	SideLocus           _preLbot;
	/// Flag to record whether a 'false' return from backtracker is due
	/// to having exceeded one or more backrtacking limits
	bool                _bailedOnBacktracks;
	/// Source of pseudo-random numbers
	RandomSource        _rand;
	/// Be talkative
	bool                _verbose;
	uint64_t            _ihits;
	// Holding area for partial alignments
	vector<PartialAlignment> _partialsBuf;
	// Current range to expose to consumers
	Range               _curRange;
	uint32_t            _patid;
	char                _primer;
	char                _trimc;
	uint32_t            _seed;
#ifndef NDEBUG
	std::set<TIndexOff> allTops_;
#endif
};

/**
 * Class that coordinates quality- and quantity-aware backtracking over
 * some range of a read sequence.
 *
 * The creator can configure the BacktrackManager to treat different
 * stretches of the read differently.
 */
class EbwtRangeSource : public RangeSource {
	typedef Ebwt<String<Dna> > TEbwt;
	typedef std::pair<int, int> TIntPair;
public:
	EbwtRangeSource(
			const TEbwt* ebwt,
			bool         fw,
			TIndexOffU   qualLim,
			bool         reportExacts,
			bool         verbose,
			bool         quiet,
			int          halfAndHalf,
			bool         partial,
			bool         maqPenalty,
			bool         qualOrder,
			AlignerMetrics *metrics = NULL) :
		RangeSource(),
		qry_(NULL),
		qlen_(0),
		qual_(NULL),
		name_(NULL),
		ebwt_(ebwt),
		fw_(fw),
		offRev0_(0),
		offRev1_(0),
		offRev2_(0),
		offRev3_(0),
		maqPenalty_(maqPenalty),
		qualOrder_(qualOrder),
		qualLim_(qualLim),
		reportExacts_(reportExacts),
		halfAndHalf_(halfAndHalf),
		partial_(partial),
		depth5_(0),
		depth3_(0),
		verbose_(verbose),
		quiet_(quiet),
		skippingThisRead_(false),
		metrics_(metrics)
	{ curEbwt_ = ebwt_; }

	/**
	 * Set a new query read.
	 */
	virtual void setQuery(Read& r, Range *seedRange) {
		const bool ebwtFw = ebwt_->fw();
		if(ebwtFw) {
			qry_  = fw_ ? &r.patFw : &r.patRc;
			qual_ = fw_ ? &r.qual  : &r.qualRev;
		} else {
			qry_  = fw_ ? &r.patFwRev : &r.patRcRev;
			qual_ = fw_ ? &r.qualRev  : &r.qual;
		}
		name_ = &r.name;
		if(seedRange != NULL) seedRange_ = *seedRange;
		else                  seedRange_.invalidate();
		qlen_ = length(*qry_);
		skippingThisRead_ = false;
		// Apply edits from the partial alignment to the query pattern
		if(seedRange_.valid()) {
			qryBuf_ = *qry_;
			const size_t srSz = seedRange_.mms.size();
			assert_gt(srSz, 0);
			assert_eq(srSz, seedRange_.refcs.size());
			for(size_t i = 0; i < srSz; i++) {
				assert_lt(seedRange_.mms[i], qlen_);
				char rc = (char)seedRange_.refcs[i];
				assert(rc == 'A' || rc == 'C' || rc == 'G' || rc == 'T');
				ASSERT_ONLY(char oc = (char)qryBuf_[qlen_ - seedRange_.mms[i] - 1]);
				assert_neq(rc, oc);
				qryBuf_[qlen_ - seedRange_.mms[i] - 1] = (Dna5)rc;
				assert_neq((Dna5)rc, (*qry_)[qlen_ - seedRange_.mms[i] - 1]);
			}
			qry_ = &qryBuf_;
		}
		// Make sure every qual is a valid qual ASCII character (>= 33)
		for(size_t i = 0; i < length(*qual_); i++) {
			assert_geq((*qual_)[i], 33);
		}
		assert_geq(length(*qual_), qlen_);
		this->done = false;
		this->foundRange = false;
		color_ = r.color;
		rand_.init(r.seed);
	}

	/**
	 * Set backtracking constraints.
	 */
	void setOffs(uint32_t depth5,   // depth of far edge of hi-half
	             uint32_t depth3,   // depth of far edge of lo-half
	             uint32_t unrevOff, // depth above which we cannot backtrack
	             uint32_t revOff1,  // depth above which we may backtrack just once
	             uint32_t revOff2,  // depth above which we may backtrack just twice
	             uint32_t revOff3)  // depth above which we may backtrack just three times
	{
		depth5_   = depth5;
		depth3_   = depth3;
		assert_geq(depth3_, depth5_);
		offRev0_ = unrevOff;
		offRev1_  = revOff1;
		offRev2_  = revOff2;
		offRev3_  = revOff3;
	}

	/**
	 * Return true iff this RangeSource is allowed to report exact
	 * alignments (exact = no edits).
	 */
	bool reportExacts() const {
		return reportExacts_;
	}

	/// Return the current range
	virtual Range& range() {
		return curRange_;
	}

	/**
	 * Set qlen_ according to parameter, except don't let it fall below
	 * the length of the query.
	 */
	void setQlen(uint32_t qlen) {
		assert(qry_ != NULL);
		qlen_ = min<uint32_t>((uint32_t)length(*qry_), qlen);
	}

	/**
	 * Initiate continuations so that the next call to advance() begins
	 * a new search.  Note that contMan is empty upon return if there
	 * are no valid continuations to begin with.  Also note that
	 * calling initConts() may result in finding a range (i.e., if we
	 * immediately jump to a valid range using the ftab).
	 */
	virtual void
	initBranch(PathManager& pm) {
		assert(curEbwt_ != NULL);
		assert_gt(length(*qry_), 0);
		assert_leq(qlen_, length(*qry_));
		assert_geq(length(*qual_), length(*qry_));
		const Ebwt<String<Dna> >& ebwt = *ebwt_;
		int ftabChars = ebwt._eh._ftabChars;
		this->foundRange = false;
		int nsInSeed = 0; int nsInFtab = 0;
		ASSERT_ONLY(allTops_.clear());
		if(skippingThisRead_) {
			this->done = true;
			return;
		}
		if(qlen_ < 4) {
			uint32_t maxmms = 0;
			if(offRev0_ != offRev1_) maxmms = 1;
			if(offRev1_ != offRev2_) maxmms = 2;
			if(offRev2_ != offRev3_) maxmms = 3;
			if(qlen_ <= maxmms) {
				if(!quiet_) {
					ThreadSafe _ts(&gLock);
					cerr << "Warning: Read (" << (*name_) << ") is less than " << (maxmms+1) << " characters long; skipping..." << endl;
				}
				this->done = true;
				skippingThisRead_ = true;
				return;
			}
		}
		if(!tallyNs(nsInSeed, nsInFtab)) {
			// No alignments are possible because of the distribution
			// of Ns in the read in combination with the backtracking
			// constraints.
			return;
		}
		// icost = total cost penalty (major bits = stratum, minor bits =
		// quality penalty) incurred so far by partial alignment
		uint16_t icost = (seedRange_.valid()) ? seedRange_.cost : 0;
		// iham = total quality penalty incurred so far by partial alignment
		uint16_t iham = (seedRange_.valid() && qualOrder_) ? (seedRange_.cost & ~0xc000): 0;
		assert_leq(iham, qualLim_);
		// m = depth beyond which ftab must not extend or else we might
		// miss some legitimate paths
		uint32_t m = min<uint32_t>(offRev0_, (uint32_t)qlen_);
		// Let skipInvalidExact = true if using the ftab would be a
		// waste because it would jump directly to an alignment we
		// couldn't use.
		bool ftabSkipsToEnd = (qlen_ == (uint32_t)ftabChars);
		bool skipInvalidExact = (!reportExacts_ && ftabSkipsToEnd);

		// If it's OK to use the ftab...
		if(nsInFtab == 0 && m >= (uint32_t)ftabChars && !skipInvalidExact) {
			// Use the ftab to jump 'ftabChars' chars into the read
			// from the right
			uint32_t ftabOff = calcFtabOff();
			TIndexOffU top = ebwt.ftabHi(ftabOff);
			TIndexOffU bot = ebwt.ftabLo(ftabOff+1);
			if(qlen_ == (uint32_t)ftabChars && bot > top) {
				// We found a range with 0 mismatches immediately.  Set
				// fields to indicate we found a range.
				assert(reportExacts_);
				curRange_.top     = top;
				curRange_.bot     = bot;
				curRange_.stratum = (icost >> 14);
				curRange_.cost    = icost;
				curRange_.numMms  = 0;
				curRange_.ebwt    = ebwt_;
				curRange_.fw      = fw_;
				curRange_.mms.clear(); // no mismatches
				curRange_.refcs.clear(); // no mismatches
				// Lump in the edits from the partial alignment
				addPartialEdits();
				assert(curRange_.repOk());
				// no need to do anything with curRange_.refcs
				this->foundRange  = true;
				//this->done = true;
				return;
			} else if (bot > top) {
				// We have a range to extend
				assert_leq(top, ebwt._eh._len);
				assert_leq(bot, ebwt._eh._len);
				Branch *b = pm.bpool.alloc();
				if(b == NULL) {
					assert(pm.empty());
					return;
				}
				if(!b->init(
				        pm.rpool, pm.epool, pm.bpool.lastId(), (uint32_t)qlen_,
				        offRev0_, offRev1_, offRev2_, offRev3_,
				        0, ftabChars, icost, iham, top, bot,
				        ebwt._eh, ebwt._ebwt))
				{
					// Negative result from b->init() indicates we ran
					// out of best-first chunk memory
					assert(pm.empty());
					return;
				}
				assert(!b->curtailed_);
				assert(!b->exhausted_);
				assert_gt(b->depth3_, 0);
				pm.push(b); // insert into priority queue
				assert(!pm.empty());
			} else {
				// The arrows are already closed within the
				// unrevisitable region; give up
			}
		} else {
			// We can't use the ftab, so we start from the rightmost
			// position and use _fchr
			Branch *b = pm.bpool.alloc();
			if(b == NULL) {
				assert(pm.empty());
				return;
			}
			if(!b->init(pm.rpool, pm.epool, pm.bpool.lastId(), (uint32_t)qlen_,
			        offRev0_, offRev1_, offRev2_, offRev3_,
			        0, 0, icost, iham, 0, 0, ebwt._eh, ebwt._ebwt))
			{
				// Negative result from b->init() indicates we ran
				// out of best-first chunk memory
				assert(pm.empty());
				return;
			}
			assert(!b->curtailed_);
			assert(!b->exhausted_);
			assert_gt(b->depth3_, 0);
			pm.push(b); // insert into priority queue
			assert(!pm.empty());
		}
		return;
	}

	/**
	 * Advance along the lowest-cost branch managed by the given
	 * PathManager.  Keep advancing until condition 'until' is
	 * satisfied.  Typically, the stopping condition 'until' is
	 * set to stop whenever pm's minCost changes.
	 */
	virtual void
	advanceBranch(int until, uint16_t minCost, PathManager& pm) {
		assert(curEbwt_ != NULL);

		// Let this->foundRange = false; we'll set it to true iff this call
		// to advance yielded a new valid-alignment range.
		this->foundRange = false;

		// Can't have already exceeded weighted hamming distance threshold
		assert_gt(length(*qry_), 0);
		assert_leq(qlen_, length(*qry_));
		assert_geq(length(*qual_), length(*qry_));
		assert(!pm.empty());

		do {
			assert(pm.repOk());
			// Get the highest-priority branch according to the priority
			// queue in 'pm'
			Branch* br = pm.front();
			// Shouldn't be curtailed or exhausted
			assert(!br->exhausted_);
			assert(!br->curtailed_);
			assert_gt(br->depth3_, 0);
			assert_leq(br->ham_, qualLim_);
			if(verbose_) {
				br->print((*qry_), (*qual_), minCost, cout, (halfAndHalf_>0), partial_, fw_, ebwt_->fw());
				if(!br->edits_.empty()) {
					cout << "Edit: ";
					for(size_t i = 0; i < br->edits_.size(); i++) {
						Edit e = br->edits_.get(i);
						cout << (curEbwt_->fw() ? (qlen_ - e.pos - 1) : e.pos)
							 << (char)e.chr;
						if(i < br->edits_.size()-1) cout << " ";
					}
					cout << endl;
				}

			}
			assert(br->repOk((uint32_t)qlen_));

			ASSERT_ONLY(int stratum = br->cost_ >> 14); // shift the stratum over
			assert_lt(stratum, 4);
			// Not necessarily true with rounding
			uint32_t depth = br->tipDepth();

			const Ebwt<String<Dna> >& ebwt = *ebwt_;

			if(halfAndHalf_ > 0) assert_gt(depth3_, depth5_);

			bool reportedPartial = false;
			bool invalidExact = false;
			bool empty = false;
			bool hit = false;
			uint16_t cost = br->cost_;
			uint32_t cur = 0;
			uint32_t nedits = 0;

			if(halfAndHalf_ && !hhCheckTop(br, depth, 0)) {
				// Stop extending this branch because it violates a half-
				// and-half constraint
				if(metrics_ != NULL) metrics_->curBacktracks_++;
				pm.curtail(br, (uint32_t)qlen_, depth3_, qualOrder_);
				goto bail;
			}

			cur = (uint32_t)(qlen_ - depth - 1); // current offset into qry_
			if(depth < qlen_) {
				// Determine whether ranges at this location are candidates
				// for backtracking
				int c = (int)(*qry_)[cur]; // get char at this position
				int nextc = -1;
				if(cur < qlen_-1) nextc = (int)(*qry_)[cur+1];
				assert_leq(c, 4);
				// If any uncalled base's penalty is still under
				// the ceiling, then this position is an alternative
				uint8_t q[4] = {'!', '!', '!', '!'};
				uint8_t bestq;
				// get unrounded penalties at this position
				bestq = q[0] = q[1] = q[2] = q[3] =
					mmPenalty(maqPenalty_, qualAt(cur));

				// The current query position is a legit alternative if it a) is
				// not in the unrevisitable region, and b) its selection would
				// not necessarily cause the quality ceiling (if one exists) to
				// be exceeded
				bool curIsAlternative = (depth >= br->depth0_) &&
				                        (br->ham_ + bestq <= qualLim_);
				ASSERT_ONLY(TIndexOffU obot = br->bot_);
				TIndexOffU otop = br->top_;

				// If c is 'N', then it's a mismatch
				if(c == 4 && depth > 0) {
					// Force the 'else if(curIsAlternative)' or 'else'
					// branches below
					br->top_ = br->bot_ = 1;
				} else if(c == 4) {
					// We'll take the 'if(br->top == 0 && br->bot == 0)'
					// branch below
					assert_eq(0, br->top_);
					assert_eq(0, br->bot_);
				}

				// Get the range state for the current position
				RangeState *rs = br->rangeState();
				assert(rs != NULL);
				// Calculate the ranges for this position
				if(br->top_ == 0 && br->bot_ == 0) {
					// Calculate first quartet of ranges using the _fchr[]
					// array
								  rs->tops[0] = ebwt._fchr[0];
					rs->bots[0] = rs->tops[1] = ebwt._fchr[1];
					rs->bots[1] = rs->tops[2] = ebwt._fchr[2];
					rs->bots[2] = rs->tops[3] = ebwt._fchr[3];
					rs->bots[3]               = ebwt._fchr[4];
					ASSERT_ONLY(int r =)
					br->installRanges(c, nextc, qualLim_ - br->ham_, q);
					assert(r < 4 || c == 4);
					// Update top and bot
					if(c < 4) {
						br->top_ = rs->tops[c];
						br->bot_ = rs->bots[c];
					}
				} else if(curIsAlternative && (br->bot_ > br->top_ || c == 4)) {
					// Calculate next quartet of ranges.  We hope that the
					// appropriate cache lines are prefetched.
					assert(br->ltop_.valid());
								  rs->tops[0] =
					rs->bots[0] = rs->tops[1] =
					rs->bots[1] = rs->tops[2] =
					rs->bots[2] = rs->tops[3] =
					rs->bots[3]               = 0;
					if(br->lbot_.valid()) {
						if(metrics_ != NULL) metrics_->curBwtOps_++;
						ebwt.mapLFEx(br->ltop_, br->lbot_, (TIndexOffU*)rs->tops, (TIndexOffU*)rs->bots);
					} else {
#ifndef NDEBUG
						TIndexOffU tmptops[] = {0, 0, 0, 0};
						TIndexOffU tmpbots[] = {0, 0, 0, 0};
						SideLocus ltop, lbot;
						ltop.initFromRow(otop, ebwt_->_eh, ebwt_->_ebwt);
						lbot.initFromRow(obot, ebwt_->_eh, ebwt_->_ebwt);
						ebwt.mapLFEx(ltop, lbot, tmptops, tmpbots);
#endif
						if(metrics_ != NULL) metrics_->curBwtOps_++;
						int cc = ebwt.mapLF1((TIndexOffU&)otop, br->ltop_);
						br->top_ = otop;
						assert(cc == -1 || (cc >= 0 && cc < 4));
						if(cc >= 0) {
							assert_lt(cc, 4);
							rs->tops[cc] = br->top_;
							rs->bots[cc] = (br->top_ + 1);
						}
#ifndef NDEBUG
						for(int i = 0; i < 4; i++) {
							assert_eq(tmpbots[i] - tmptops[i],
							          rs->bots[i] - rs->tops[i]);
						}
#endif
					}
					ASSERT_ONLY(int r =)
					br->installRanges(c, nextc, qualLim_ - br->ham_, q);
					assert(r < 4 || c == 4);
					// Update top and bot
					if(c < 4) {
						br->top_ = rs->tops[c];
						br->bot_ = rs->bots[c];
					} else {
						br->top_ = br->bot_ = 1;
					}
				} else if(br->bot_ > br->top_) {
					// This read position is not a legitimate backtracking
					// alternative.  No need to do the bookkeeping for the
					// entire quartet, just do c.  We hope that the
					// appropriate cache lines are prefetched before now;
					// otherwise, we're about to take an expensive cache
					// miss.
					assert(br->ltop_.valid());
					rs->eliminated_ = true; // eliminate all alternatives leaving this node
					assert(br->eliminated(br->len_));
					if(c < 4) {
						if(br->top_ + 1 == br->bot_) {
							if(metrics_ != NULL) metrics_->curBwtOps_++;
							br->bot_ = br->top_ = ebwt.mapLF1(br->top_, br->ltop_, c);
							if(br->bot_ != OFF_MASK) br->bot_++;
						} else {
							if(metrics_ != NULL) metrics_->curBwtOps_++;
							br->top_ = ebwt.mapLF(br->ltop_, c);
							assert(br->lbot_.valid());
							if(metrics_ != NULL) metrics_->curBwtOps_++;
							br->bot_ = ebwt.mapLF(br->lbot_, c);
						}
					}
				} else {
					rs->eliminated_ = true;
				}
				assert(rs->repOk());
				// br->top_ and br->bot_ now contain the next top and bot
			} else {
				// The continuation had already processed the whole read
				assert_eq(qlen_, depth);
				cur = 0;
			}
			empty = (br->top_ == br->bot_);
			hit = (cur == 0 && !empty);

			// Check whether we've obtained an exact alignment when
			// we've been instructed not to report exact alignments
			nedits = (uint32_t)br->edits_.size();
			invalidExact = (hit && nedits == 0 && !reportExacts_);
			assert_leq(br->ham_, qualLim_);

			// Set this to true if the only way to make legal progress
			// is via one or more additional backtracks.
			if(halfAndHalf_ && !hhCheck(br, depth, empty)) {
				// This alignment doesn't satisfy the half-and-half
				// requirements; reject it
				if(metrics_ != NULL) metrics_->curBacktracks_++;
				pm.curtail(br, (uint32_t)qlen_, depth3_, qualOrder_);
				goto bail;
			}

			if(hit &&            // there is a range to report
			   !invalidExact &&  // not disqualified by no-exact-hits setting
			   !reportedPartial) // not an already-reported partial alignment
			{
				if(verbose_) {
					if(partial_) {
						cout << " Partial alignment:" << endl;
					} else {
						cout << " Final alignment:" << endl;
					}
					br->len_++;
					br->print((*qry_), (*qual_), minCost, cout, halfAndHalf_ > 0, partial_, fw_, ebwt_->fw());
					br->len_--;
					cout << endl;
				}
				assert_gt(br->bot_, br->top_);
				assert_leq(br->ham_, qualLim_);
				assert_leq((uint32_t)(br->cost_ & ~0xc000), qualLim_);
				if(metrics_ != NULL) metrics_->setReadHasRange();
				curRange_.top     = br->top_;
				curRange_.bot     = br->bot_;
				curRange_.cost    = br->cost_;
				curRange_.stratum = (br->cost_ >> 14);
				curRange_.numMms  = nedits;
				curRange_.fw      = fw_;
				curRange_.mms.clear();
				curRange_.refcs.clear();
				for(size_t i = 0; i < nedits; i++) {
					curRange_.mms.push_back((uint32_t)(qlen_ - br->edits_.get(i).pos - 1));
					curRange_.refcs.push_back((char)br->edits_.get(i).chr);
				}
				addPartialEdits();
				curRange_.ebwt    = ebwt_;
				this->foundRange  = true;
	#ifndef NDEBUG
				TIndexOff top2 = (TIndexOff)br->top_;
				top2++; // ensure it's not 0
				if(ebwt_->fw()) top2 = -top2;
				assert(allTops_.find(top2) == allTops_.end());
				allTops_.insert(top2);
	#endif
				assert(curRange_.repOk());
				// Must curtail because we've consumed the whole pattern
				if(metrics_ != NULL) metrics_->curBacktracks_++;
				pm.curtail(br, (uint32_t)qlen_, depth3_, qualOrder_);
			} else if(empty || cur == 0) {
				// The branch couldn't be extended further
				if(metrics_ != NULL) metrics_->curBacktracks_++;
				pm.curtail(br, (uint32_t)qlen_, depth3_, qualOrder_);
			} else {
				// Extend the branch by one position; no change to its cost
				// so there's no need to reconsider where it lies in the
				// priority queue
				assert_neq(0, cur);
				br->extend();
			}
		bail:
			// Make sure the front element of the priority queue is
			// extendable (i.e. not curtailed) and then prep it.
			if(!pm.splitAndPrep(rand_, (uint32_t)qlen_, qualLim_, depth3_,
			                    qualOrder_,
			                    ebwt_->_eh, ebwt_->_ebwt, ebwt_->_fw))
			{
				pm.reset(0);
				assert(pm.empty());
			}
			if(pm.empty()) {
				// No more branches
				break;
			}
			assert(!pm.front()->curtailed_);
			assert(!pm.front()->exhausted_);

			if(until == ADV_COST_CHANGES && pm.front()->cost_ != cost) break;
			else if(until == ADV_STEP) break;

		} while(!this->foundRange);
		if(!pm.empty()) {
			assert(!pm.front()->curtailed_);
			assert(!pm.front()->exhausted_);
		}
	}

	/**
	 * Return true iff we're enforcing a half-and-half constraint
	 * (forced edits in both seed halves).
	 */
	int halfAndHalf() const {
		return halfAndHalf_;
	}

protected:

	/**
	 * Lump all the seed-alignment edits from the seedRange_ range
	 * found previously to the curRange_ range just found.
	 */
	void addPartialEdits() {
		// Lump in the edits from the partial alignment
		if(seedRange_.valid()) {
			const size_t srSz = seedRange_.mms.size();
			for(size_t i = 0; i < srSz; i++) {
				curRange_.mms.push_back((uint32_t)(qlen_ - seedRange_.mms[i] - 1));
				curRange_.refcs.push_back(seedRange_.refcs[i]);
			}
			curRange_.numMms += srSz;
		}
	}

	/**
	 * Return true iff we're OK to continue after considering which
	 * half-seed boundary we're passing through, together with the
	 * number of mismatches accumulated so far.  Return false if we
	 * should stop because a half-and-half constraint is violated.  If
	 * we're not currently passing a half-seed boundary, just return
	 * true.
	 */
	bool hhCheck(Branch *b, uint32_t depth, bool empty) {
		const uint32_t nedits = (uint32_t)b->edits_.size();
		ASSERT_ONLY(uint32_t lim3 = (offRev3_ == offRev2_)? 2 : 3);
		ASSERT_ONLY(uint32_t lim5 = (offRev1_ == offRev0_)? 2 : 1);
		if((depth == (depth5_-1)) && !empty) {
			// We're crossing the boundary separating the hi-half
			// from the non-seed portion of the read.
			// We should induce a mismatch if we haven't mismatched
			// yet, so that we don't waste time pursuing a match
			// that was covered by a previous phase
			assert_leq(nedits, lim5);
			return nedits > 0;
		} else if((depth == (depth3_-1)) && !empty) {
			// We're crossing the boundary separating the lo-half
			// from the non-seed portion of the read
			assert_leq(nedits, lim3);
			assert_gt(nedits, 0);
			// Count the mismatches in the lo and hi halves
			uint32_t loHalfMms = 0, hiHalfMms = 0;
			for(size_t i = 0; i < nedits; i++) {
				uint32_t depth = b->edits_.get(i).pos;
				if     (depth < depth5_) hiHalfMms++;
				else if(depth < depth3_) loHalfMms++;
				else assert(false);
			}
			assert_leq(loHalfMms + hiHalfMms, lim3);
			bool invalidHalfAndHalf = (loHalfMms == 0 || hiHalfMms == 0);
			return (nedits >= (uint32_t)halfAndHalf_ && !invalidHalfAndHalf);
		}
#ifndef NDEBUG
		if(depth < depth5_-1) {
			assert_leq(nedits, lim5);
		}
		else if(depth >= depth5_ && depth < depth3_-1) {
			assert_gt(nedits, 0);
			assert_leq(nedits, lim3);
		}
#endif
		return true;
	}

	/**
	 * Return true iff the state of the backtracker as encoded by
	 * stackDepth, d and iham is compatible with the current half-and-
	 * half alignment mode.  prehits is for sanity checking when
	 * bailing.
	 */
	bool hhCheckTop(Branch* b,
	                uint32_t d,
	                uint32_t iham,
	                uint64_t prehits = 0xffffffffffffffffllu)
	{
		// Crossing from the hi-half into the lo-half
		ASSERT_ONLY(uint32_t lim3 = (offRev3_ == offRev2_)? 2 : 3);
		ASSERT_ONLY(uint32_t lim5 = (offRev1_ == offRev0_)? 2 : 1);
		const uint32_t nedits = (uint32_t)b->edits_.size();
		if(d == depth5_) {
			assert_leq(nedits, lim5);
			if(nedits == 0) {
				return false;
			}
		} else if(d == depth3_) {
			assert_leq(nedits, lim3);
			if(nedits < (uint32_t)halfAndHalf_) {
				return false;
			}
		}
#ifndef NDEBUG
		else {
			// We didn't just cross a boundary, so do an in-between check
			if(d >= depth5_) {
				assert_geq(nedits, 1);
			} else if(d >= depth3_) {
				assert_geq(nedits, lim3);
			}
		}
#endif
		return true;
	}

	/**
	 * Return the Phred-scale quality value at position 'off'
	 */
	inline uint8_t qualAt(size_t off) {
		return phredCharToPhredQual((*qual_)[off]);
	}

	/**
	 * Tally how many Ns occur in the seed region and in the ftab-
	 * jumpable region of the read.  Check whether the mismatches
	 * induced by the Ns already violates the current policy.  Return
	 * false if the policy is already violated, true otherwise.
	 */
	bool tallyNs(int& nsInSeed, int& nsInFtab) {
		const Ebwt<String<Dna> >& ebwt = *ebwt_;
		int ftabChars = ebwt._eh._ftabChars;
		// Count Ns in the seed region of the read and short-circuit if
		// the configuration of Ns guarantees that there will be no
		// valid alignments given the backtracking constraints.
		for(size_t i = 0; i < offRev3_; i++) {
			if((int)(*qry_)[qlen_-i-1] == 4) {
				nsInSeed++;
				if(nsInSeed == 1) {
					if(i < offRev0_) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else if(nsInSeed == 2) {
					if(i < offRev1_) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else if(nsInSeed == 3) {
					if(i < offRev2_) {
						return false; // Exceeded mm budget on Ns alone
					}
				} else {
					assert_gt(nsInSeed, 3);
					return false;     // Exceeded mm budget on Ns alone
				}
			}
		}
		// Calculate the number of Ns there are in the region that
		// would get jumped over if the ftab were used.
		for(size_t i = 0; i < (size_t)ftabChars && i < qlen_; i++) {
			if((int)(*qry_)[qlen_-i-1] == 4) nsInFtab++;
		}
		return true;
	}

	/**
	 * Calculate the offset into the ftab for the rightmost 'ftabChars'
	 * characters of the current query. Rightmost char gets least
	 * significant bit-pair.
	 */
	uint32_t calcFtabOff() {
		const Ebwt<String<Dna> >& ebwt = *ebwt_;
		int ftabChars = ebwt._eh._ftabChars;
		uint32_t ftabOff = (*qry_)[qlen_ - ftabChars];
		assert_lt(ftabOff, 4);
		assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		for(int i = ftabChars - 1; i > 0; i--) {
			ftabOff <<= 2;
			assert_lt((uint32_t)(*qry_)[qlen_-i], 4);
			ftabOff |= (uint32_t)(*qry_)[qlen_-i];
			assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		}
		assert_lt(ftabOff, ebwt._eh._ftabLen-1);
		return ftabOff;
	}

	String<Dna5>*       qry_;    // query (read) sequence
	String<Dna5>        qryBuf_; // for composing modified qry_ strings
	size_t              qlen_;   // length of _qry
	String<char>*       qual_;   // quality values for _qry
	String<char>*       name_;   // name of _qry
	bool                color_;  // true -> read is colorspace
	const Ebwt<String<Dna> >*   ebwt_;   // Ebwt to search in
	bool                fw_;
	uint32_t            offRev0_; // unrevisitable chunk
	uint32_t            offRev1_;  // 1-revisitable chunk
	uint32_t            offRev2_;  // 2-revisitable chunk
	uint32_t            offRev3_;  // 3-revisitable chunk
	/// Whether to round qualities off Maq-style when calculating penalties
	bool                maqPenalty_;
	/// Whether to order paths on our search in a way that takes
	/// qualities into account.  If this is false, the effect is that
	/// the first path reported is guaranteed to be in the best
	/// stratum, but it's not guaranteed to have the best quals.
	bool                qualOrder_;
	/// Reject alignments where sum of qualities at mismatched
	/// positions is greater than qualLim_
	uint32_t            qualLim_;
	/// Report exact alignments iff this is true
	bool                reportExacts_;
	/// Whether to use the _os array together with a naive matching
	/// algorithm to double-check reported alignments (or the lack
	/// thereof)
	int                 halfAndHalf_;
	/// Whether we're generating partial alignments for a longer
	/// alignment in the opposite index.
	bool                partial_;
	/// Depth of 5'-seed-half border
	uint32_t            depth5_;
	/// Depth of 3'-seed-half border
	uint32_t            depth3_;
	/// Source of pseudo-random numbers
	RandomSource        rand_;
	/// Be talkative
	bool                verbose_;
	/// Suppress unnecessary output
	bool                quiet_;
	// Current range to expose to consumers
	Range               curRange_;
	// Range for the partial alignment we're extending (NULL if we
	// aren't extending a partial)
	Range               seedRange_;
	// Starts as false; set to true as soon as we know we want to skip
	// all further processing of this read
	bool                skippingThisRead_;
	// Object encapsulating metrics
	AlignerMetrics*     metrics_;
#ifndef NDEBUG
	std::set<TIndexOff> allTops_;
#endif
};

/**
 * Concrete factory for EbwtRangeSource objects.
 */
class EbwtRangeSourceFactory {
	typedef Ebwt<String<Dna> > TEbwt;
public:
	EbwtRangeSourceFactory(
			const TEbwt* ebwt,
			bool         fw,
			uint32_t     qualThresh,
			bool         reportExacts,
			bool         verbose,
			bool         quiet,
			bool         halfAndHalf,
			bool         seeded,
			bool         maqPenalty,
			bool         qualOrder,
			AlignerMetrics *metrics = NULL) :
			ebwt_(ebwt),
			fw_(fw),
			qualThresh_(qualThresh),
			reportExacts_(reportExacts),
			verbose_(verbose),
			quiet_(quiet),
			halfAndHalf_(halfAndHalf),
			seeded_(seeded),
			maqPenalty_(maqPenalty),
			qualOrder_(qualOrder),
			metrics_(metrics) { }

	/**
	 * Return new EbwtRangeSource with predefined params.s
	 */
	EbwtRangeSource *create() {
		return new EbwtRangeSource(ebwt_, fw_, qualThresh_,
		                           reportExacts_, verbose_, quiet_,
		                           halfAndHalf_, seeded_, maqPenalty_,
		                           qualOrder_, metrics_);
	}

protected:
	const TEbwt* ebwt_;
	bool         fw_;
	uint32_t     qualThresh_;
	bool         reportExacts_;
	bool         verbose_;
	bool         quiet_;
	bool         halfAndHalf_;
	bool         seeded_;
	bool         maqPenalty_;
	bool         qualOrder_;
	AlignerMetrics *metrics_;
};

/**
 * What boundary within the alignment to "pin" a particular
 * backtracking constraint to.
 */
enum SearchConstraintExtent {
	PIN_TO_BEGINNING = 1, // depth 0; i.e., constraint is inactive
	PIN_TO_LEN,           // constraint applies to while alignment
	PIN_TO_HI_HALF_EDGE,  // constraint applies to hi-half of seed region
	PIN_TO_SEED_EDGE      // constraint applies to entire seed region
};

/**
 * Concrete RangeSourceDriver that deals properly with
 * GreedyDFSRangeSource by calling setOffs() with the appropriate
 * parameters when initializing it;
 */
class EbwtRangeSourceDriver :
	public SingleRangeSourceDriver<EbwtRangeSource>
{
public:
	EbwtRangeSourceDriver(
			EbwtSearchParams<String<Dna> >& params,
			EbwtRangeSource* rs,
			bool fw,
			bool seed,
			bool maqPenalty,
			bool qualOrder,
			HitSink& sink,
			HitSinkPerThread* sinkPt,
			uint32_t seedLen,
			bool nudgeLeft,
			SearchConstraintExtent rev0Off,
			SearchConstraintExtent rev1Off,
			SearchConstraintExtent rev2Off,
			SearchConstraintExtent rev3Off,
			vector<String<Dna5> >& os,
			bool verbose,
			bool quiet,
			bool mate1,
			ChunkPool* pool,
			int *btCnt) :
			SingleRangeSourceDriver<EbwtRangeSource>(
					params, rs, fw, sink, sinkPt, os, verbose,
					quiet, mate1, 0, pool, btCnt),
			seed_(seed),
			maqPenalty_(maqPenalty),
			qualOrder_(qualOrder),
			rs_(rs), seedLen_(seedLen),
			nudgeLeft_(nudgeLeft),
			rev0Off_(rev0Off), rev1Off_(rev1Off),
			rev2Off_(rev2Off), rev3Off_(rev3Off),
			verbose_(verbose), quiet_(quiet)
	{
		if(seed_) assert_gt(seedLen, 0);
	}

	virtual ~EbwtRangeSourceDriver() { }

	bool seed() const { return seed_; }

	bool ebwtFw() const { return rs_->curEbwt()->fw(); }

	/**
	 * Called every time setQuery() is called in the parent class,
	 * after setQuery() has been called on the RangeSource but before
	 * initConts() has been called.
	 */
	virtual void initRangeSource(const String<char>& qual) {
		// If seedLen_ is huge, then it will always cover the whole
		// alignment
		assert_eq(len_, seqan::length(qual));
		uint32_t s = (seedLen_ > 0 ? min(seedLen_, len_) : len_);
		uint32_t sLeft  = s >> 1;
		uint32_t sRight = s >> 1;
		// If seed has odd length, then nudge appropriate half up by 1
		if((s & 1) != 0) { if(nudgeLeft_) sLeft++; else sRight++; }
		uint32_t rev0Off = cextToDepth(rev0Off_, sRight, s, len_);
		uint32_t rev1Off = cextToDepth(rev1Off_, sRight, s, len_);
		uint32_t rev2Off = cextToDepth(rev2Off_, sRight, s, len_);
		uint32_t rev3Off = cextToDepth(rev3Off_, sRight, s, len_);
		// Truncate the pattern if necessary
		uint32_t qlen = (uint32_t)seqan::length(qual);
		if(seed_) {
			if(len_ > s) {
				rs_->setQlen(s);
				qlen = s;
			}
			assert(!rs_->reportExacts());
		}
		// If there are any Ns in the unrevisitable region, then this
		// driver is guaranteed to yield no fruit.
		uint16_t minCost = 0;
		if(rs_->reportExacts()) {
			// Keep minCost at 0
		} else if (!rs_->halfAndHalf() && rev0Off < s) {
			// Exacts not allowed, so there must be at least 1 mismatch
			// outside of the unrevisitable area
			minCost = 1 << 14;
			if(qualOrder_) {
				uint8_t lowQual = 0xff;
				for(uint32_t d = rev0Off; d < s; d++) {
					uint8_t lowAtPos;
					lowAtPos = qual[qlen - d - 1];
					if(lowAtPos < lowQual) lowQual = lowAtPos;
				}
				assert_lt(lowQual, 0xff);
				minCost += mmPenalty(maqPenalty_, phredCharToPhredQual(lowQual));
			}
		} else if(rs_->halfAndHalf() && sRight > 0 && sRight < (s-1)) {
			// Half-and-half constraints are active, so there must be
			// at least 1 mismatch in both halves of the seed
			assert(rs_->halfAndHalf());
			minCost = (seed_ ? 3 : 2) << 14;
			if(qualOrder_) {
				assert(rs_->halfAndHalf() == 2 || rs_->halfAndHalf() == 3);
				uint8_t lowQual1 = 0xff;
				for(uint32_t d = 0; d < sRight; d++) {
					uint8_t lowAtPos;
					lowAtPos = qual[qlen - d - 1];
					if(lowAtPos < lowQual1) lowQual1 = lowAtPos;
				}
				assert_lt(lowQual1, 0xff);
				minCost += mmPenalty(maqPenalty_, phredCharToPhredQual(lowQual1));
				uint8_t lowQual2_1 = 0xff;
				uint8_t lowQual2_2 = 0xff;
				for(uint32_t d = sRight; d < s; d++) {
					uint8_t lowAtPos;
					lowAtPos = qual[qlen - d - 1];
					if(lowAtPos < lowQual2_1) {
						if(lowQual2_1 != 0xff) {
							lowQual2_2 = lowQual2_1;
						}
						lowQual2_1 = lowAtPos;
					} else if(lowAtPos < lowQual2_2) {
						lowQual2_2 = lowAtPos;
					}
				}
				assert_lt(lowQual2_1, 0xff);
				minCost += mmPenalty(maqPenalty_, phredCharToPhredQual(lowQual2_1));
				if(rs_->halfAndHalf() > 2 && lowQual2_2 != 0xff) {
					minCost += mmPenalty(maqPenalty_, phredCharToPhredQual(lowQual2_2));
				}
			}
		}
		if(verbose_) cout << "initRangeSource minCost: " << minCost << endl;
		this->minCostAdjustment_ = minCost;
		rs_->setOffs(sRight,   // depth of far edge of hi-half (only matters where half-and-half is possible)
		             s,        // depth of far edge of lo-half (only matters where half-and-half is possible)
		             rev0Off,  // depth above which we cannot backtrack
		             rev1Off,  // depth above which we may backtrack just once
		             rev2Off,  // depth above which we may backtrack just twice
		             rev3Off); // depth above which we may backtrack just three times
	}

protected:

	/**
	 * Convert a search constraint extent to an actual depth into the
	 * read.
	 */
	inline uint32_t cextToDepth(SearchConstraintExtent cext,
	                            uint32_t sRight,
	                            uint32_t s,
	                            uint32_t len)
	{
		if(cext == PIN_TO_SEED_EDGE)    return s;
		if(cext == PIN_TO_HI_HALF_EDGE) return sRight;
		if(cext == PIN_TO_BEGINNING)    return 0;
		if(cext == PIN_TO_LEN)          return len;
		cerr << "Bad SearchConstraintExtent: " << cext;
		throw 1;
	}

	bool seed_;
	bool maqPenalty_;
	bool qualOrder_;
	EbwtRangeSource* rs_;
	uint32_t seedLen_;
	bool nudgeLeft_;
	SearchConstraintExtent rev0Off_;
	SearchConstraintExtent rev1Off_;
	SearchConstraintExtent rev2Off_;
	SearchConstraintExtent rev3Off_;
	bool verbose_;
	bool quiet_;
};

/**
 * Concrete RangeSourceDriver that deals properly with
 * GreedyDFSRangeSource by calling setOffs() with the appropriate
 * parameters when initializing it;
 */
class EbwtRangeSourceDriverFactory {
public:
	EbwtRangeSourceDriverFactory(
			EbwtSearchParams<String<Dna> >& params,
			EbwtRangeSourceFactory* rs,
			bool fw,
			bool seed,
			bool maqPenalty,
			bool qualOrder,
			HitSink& sink,
			HitSinkPerThread* sinkPt,
			uint32_t seedLen,
			bool nudgeLeft,
			SearchConstraintExtent rev0Off,
			SearchConstraintExtent rev1Off,
			SearchConstraintExtent rev2Off,
			SearchConstraintExtent rev3Off,
			vector<String<Dna5> >& os,
			bool verbose,
			bool quiet,
			bool mate1,
			ChunkPool* pool,
			int *btCnt = NULL) :
			params_(params),
			rs_(rs),
			fw_(fw),
			seed_(seed),
			maqPenalty_(maqPenalty),
			qualOrder_(qualOrder),
			sink_(sink),
			sinkPt_(sinkPt),
			seedLen_(seedLen),
			nudgeLeft_(nudgeLeft),
			rev0Off_(rev0Off),
			rev1Off_(rev1Off),
			rev2Off_(rev2Off),
			rev3Off_(rev3Off),
			os_(os),
			verbose_(verbose),
			quiet_(quiet),
			mate1_(mate1),
			pool_(pool),
			btCnt_(btCnt)
	{ }

	~EbwtRangeSourceDriverFactory() {
		delete rs_; rs_ = NULL;
	}

	/**
	 * Return a newly-allocated EbwtRangeSourceDriver with the given
	 * parameters.
	 */
	EbwtRangeSourceDriver *create() const {
		return new EbwtRangeSourceDriver(
				params_, rs_->create(), fw_, seed_, maqPenalty_,
				qualOrder_, sink_, sinkPt_, seedLen_, nudgeLeft_,
				rev0Off_, rev1Off_, rev2Off_, rev3Off_, os_, verbose_,
				quiet_, mate1_, pool_, btCnt_);
	}

protected:
	EbwtSearchParams<String<Dna> >& params_;
	EbwtRangeSourceFactory* rs_;
	bool fw_;
	bool seed_;
	bool maqPenalty_;
	bool qualOrder_;
	HitSink& sink_;
	HitSinkPerThread* sinkPt_;
	uint32_t seedLen_;
	bool nudgeLeft_;
	SearchConstraintExtent rev0Off_;
	SearchConstraintExtent rev1Off_;
	SearchConstraintExtent rev2Off_;
	SearchConstraintExtent rev3Off_;
	vector<String<Dna5> >& os_;
	bool verbose_;
	bool quiet_;
	bool mate1_;
	ChunkPool* pool_;
	int *btCnt_;
};

/**
 * A RangeSourceDriver that manages two child EbwtRangeSourceDrivers,
 * one for searching for seed strings with mismatches in the hi-half,
 * and one for extending those seed strings toward the 3' end.
 */
class EbwtSeededRangeSourceDriver : public RangeSourceDriver<EbwtRangeSource> {
	typedef RangeSourceDriver<EbwtRangeSourceDriver>* TRangeSrcDrPtr;
	typedef CostAwareRangeSourceDriver<EbwtRangeSource> TCostAwareRangeSrcDr;
public:
	EbwtSeededRangeSourceDriver(
			EbwtRangeSourceDriverFactory* rsFact,
			EbwtRangeSourceDriver* rsSeed,
			bool fw,
			uint32_t seedLen,
			bool verbose,
			bool quiet,
			bool mate1) :
			RangeSourceDriver<EbwtRangeSource>(true, 0),
			rsFact_(rsFact), rsFull_(false, NULL, verbose, quiet, true),
			rsSeed_(rsSeed), patsrc_(NULL), seedLen_(seedLen), fw_(fw),
			mate1_(mate1), seedRange_(0)
	{
		assert(rsSeed_->seed());
	}

	virtual ~EbwtSeededRangeSourceDriver() {
		delete rsFact_; rsFact_ = NULL;
		delete rsSeed_; rsSeed_ = NULL;
	}

	/**
	 * Prepare this aligner for the next read.
	 */
	virtual void setQueryImpl(PatternSourcePerThread* patsrc, Range *partial) {
		this->done = false;
		rsSeed_->setQuery(patsrc, partial);
		this->minCostAdjustment_ = max(rsSeed_->minCostAdjustment_, rsSeed_->minCost);
		this->minCost = this->minCostAdjustment_;
		rsFull_.clearSources();
		rsFull_.setQuery(patsrc, partial);
		rsFull_.minCost = this->minCost;
		assert_gt(rsFull_.minCost, 0);
		patsrc_ = patsrc;
		// The minCostAdjustment comes from the seed range source
		// driver, based on Ns and quals in the hi-half
		this->foundRange = false;
		ASSERT_ONLY(allTops_.clear());
		assert_eq(this->minCost, min<uint16_t>(rsSeed_->minCost, rsFull_.minCost));
	}

	/**
	 * Advance the aligner by one memory op.  Return true iff we're
	 * done with this read.
	 */
	virtual void advance(int until) {
		assert(!this->foundRange);
		until = max<int>(until, ADV_COST_CHANGES);
		ASSERT_ONLY(uint16_t preCost = this->minCost);
		advanceImpl(until);
		if(this->foundRange) {
			assert_eq(range().cost, preCost);
		}
#ifndef NDEBUG
		if(this->foundRange) {
			// Assert that we have not yet dished out a range with this
			// top offset
			assert_gt(range().bot, range().top);
			assert(range().ebwt != NULL);
			TIndexOff top = (TIndexOff)range().top;
			top++; // ensure it's not 0
			if(!range().ebwt->fw()) top = -top;
			assert(allTops_.find(top) == allTops_.end());
			allTops_.insert(top);
		}
#endif
	}

	/**
	 * Advance the aligner by one memory op.  Return true iff we're
	 * done with this read.
	 */
	virtual void advanceImpl(int until) {
		assert(!this->done);
		assert(!this->foundRange);
		assert_gt(rsFull_.minCost, 0);
		// Advance the seed range source
		if(rsSeed_->done && rsFull_.done &&
		   !rsSeed_->foundRange && !rsFull_.foundRange)
		{
			this->done = true;
			return;
		}
		if(rsSeed_->done && !rsSeed_->foundRange) {
			rsSeed_->minCost = 0xffff;
			if(rsFull_.minCost > this->minCost) {
				this->minCost = rsFull_.minCost;
				// Cost changed, so return
				return;
			}
		}
		if(rsFull_.done && !rsFull_.foundRange) {
			rsFull_.minCost = 0xffff;
			if(rsSeed_->minCost > this->minCost) {
				this->minCost = rsSeed_->minCost;
				// Cost changed, so return
				return;
			}
		}
		assert(rsSeed_->minCost != 0xffff || rsFull_.minCost != 0xffff);
		// Extend a partial alignment
		ASSERT_ONLY(uint16_t oldMinCost = this->minCost);
		assert_eq(this->minCost, min<uint16_t>(rsSeed_->minCost, rsFull_.minCost));
		bool doFull = rsFull_.minCost <= rsSeed_->minCost;
		if(!doFull) {
			// Advance the partial-alignment generator
			assert_eq(rsSeed_->minCost, this->minCost);
			if(!rsSeed_->foundRange) {
				rsSeed_->advance(until);
			}
			if(rsSeed_->foundRange) {
				assert_eq(this->minCost, rsSeed_->range().cost);
				assert_eq(oldMinCost, rsSeed_->range().cost);
				seedRange_ = &rsSeed_->range();
				rsSeed_->foundRange = false;
				assert_geq(seedRange_->cost, this->minCostAdjustment_);
				this->minCostAdjustment_ = seedRange_->cost;
				assert_gt(seedRange_->numMms, 0);
				// Keep the range for the hi-half partial alignment so
				// that the driver can (a) modify the pattern string
				// and (b) modify results from the RangeSource to
				// include these edits.
				EbwtRangeSourceDriver *partial = rsFact_->create();
				partial->minCost = seedRange_->cost;
				rsFull_.minCost = seedRange_->cost;
				rsFull_.addSource(partial, seedRange_);
				if(rsFull_.foundRange) {
					this->foundRange = true;
					rsFull_.foundRange = false;
					assert(rsFull_.range().repOk());
					assert_eq(range().cost, oldMinCost);
				}
			}
			if(rsSeed_->minCost > this->minCost) {
				this->minCost = rsSeed_->minCost;
				if(!rsFull_.done) {
					this->minCost = min(this->minCost, rsFull_.minCost);
					assert_eq(this->minCost, min<uint16_t>(rsSeed_->minCost, rsFull_.minCost));
				}
			}
		} else {
			// Extend a full alignment
			assert(!rsFull_.done);
			assert(!rsFull_.foundRange);
			uint16_t oldFullCost = rsFull_.minCost;
			if(!rsFull_.foundRange) {
				rsFull_.advance(until);
			}
			// Found a minimum-cost range
			if(rsFull_.foundRange) {
				this->foundRange = true;
				rsFull_.foundRange = false;
				assert(rsFull_.range().repOk());
				assert_eq(range().cost, oldMinCost);
			}
			assert_geq(rsFull_.minCost, oldFullCost);
			// Did the min cost change?
			if(rsFull_.minCost > oldFullCost) {
				// If a range was found, hold on to it and save it for
				// later.  Update the minCost.
				assert(!rsSeed_->done || rsSeed_->minCost == 0xffff);
				this->minCost = min(rsFull_.minCost, rsSeed_->minCost);
			}
		}
	}

	/**
	 * Return the range found.
	 */
	virtual Range& range() {
		Range& r = rsFull_.range();
		r.fw = fw_;
		r.mate1 = mate1_;
		return r;
	}

	/**
	 * Return whether we're generating ranges for the first or the
	 * second mate.
	 */
	virtual bool mate1() const {
		return mate1_;
	}

	/**
	 * Return true iff current pattern is forward-oriented.
	 */
	virtual bool fw() const {
		return fw_;
	}

protected:

	EbwtRangeSourceDriverFactory* rsFact_;
	TCostAwareRangeSrcDr rsFull_;
	EbwtRangeSourceDriver* rsSeed_;
	PatternSourcePerThread* patsrc_;
	uint32_t seedLen_;
	bool fw_;
	bool mate1_;
	bool generating_;
	Range *seedRange_;
};

#endif /*EBWT_SEARCH_BACKTRACK_H_*/