/*
* Copyright 2011, Ben Langmead <langmea@cs.jhu.edu>
*
* This file is part of Bowtie 2.
*
* Bowtie 2 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Bowtie 2 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Bowtie 2. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef DIFF_SAMPLE_H_
#define DIFF_SAMPLE_H_
#include <stdint.h>
#include <string.h>
#include "assert_helpers.h"
#include "multikey_qsort.h"
#include "timer.h"
#include "ds.h"
#include "mem_ids.h"
#include "ls.h"
#include "btypes.h"
using namespace std;
#ifndef VMSG_NL
#define VMSG_NL(...) \
if(this->verbose()) { \
stringstream tmp; \
tmp << __VA_ARGS__ << endl; \
this->verbose(tmp.str()); \
}
#endif
#ifndef VMSG
#define VMSG(...) \
if(this->verbose()) { \
stringstream tmp; \
tmp << __VA_ARGS__; \
this->verbose(tmp.str()); \
}
#endif
/**
* Routines for calculating, sanity-checking, and dispensing difference
* cover samples to clients.
*/
/**
*
*/
struct sampleEntry {
uint32_t maxV;
uint32_t numSamples;
uint32_t samples[128];
};
/// Array of Colbourn and Ling calculated difference covers up to
/// r = 16 (maxV = 5953)
extern struct sampleEntry clDCs[16];
extern bool clDCs_calced; /// have clDCs been calculated?
/**
* Check that the given difference cover 'ds' actually covers all
* differences for a periodicity of v.
*/
template<typename T>
static bool dcRepOk(T v, EList<T>& ds) {
// diffs[] records all the differences observed
AutoArray<bool> covered(v, EBWT_CAT);
for(T i = 1; i < v; i++) {
covered[i] = false;
}
for(T di = T(); di < ds.size(); di++) {
for(T dj = di+1; dj < ds.size(); dj++) {
assert_lt(ds[di], ds[dj]);
T d1 = (ds[dj] - ds[di]);
T d2 = (ds[di] + v - ds[dj]);
assert_lt(d1, v);
assert_lt(d2, v);
covered[d1] = true;
covered[d2] = true;
}
}
bool ok = true;
for(T i = 1; i < v; i++) {
if(covered[i] == false) {
ok = false;
break;
}
}
return ok;
}
/**
* Return true iff each element of ts (with length 'limit') is greater
* than the last.
*/
template<typename T>
static bool increasing(T* ts, size_t limit) {
for(size_t i = 0; i < limit-1; i++) {
if(ts[i+1] <= ts[i]) return false;
}
return true;
}
/**
* Return true iff the given difference cover covers difference 'diff'
* mod 'v'.
*/
template<typename T>
static inline bool hasDifference(T *ds, T d, T v, T diff) {
// diffs[] records all the differences observed
for(T di = T(); di < d; di++) {
for(T dj = di+1; dj < d; dj++) {
assert_lt(ds[di], ds[dj]);
T d1 = (ds[dj] - ds[di]);
T d2 = (ds[di] + v - ds[dj]);
assert_lt(d1, v);
assert_lt(d2, v);
if(d1 == diff || d2 == diff) return true;
}
}
return false;
}
/**
* Exhaustively calculate optimal difference cover samples for v = 4,
* 8, 16, 32, 64, 128, 256 and store results in p2DCs[]
*/
template<typename T>
void calcExhaustiveDC(T i, bool verbose = false, bool sanityCheck = false) {
T v = i;
AutoArray<bool> diffs(v, EBWT_CAT);
// v is the target period
T ld = (T)ceil(sqrt(v));
// ud is the upper bound on |D|
T ud = v / 2;
// for all possible |D|s
bool ok = true;
T *ds = NULL;
T d;
for(d = ld; d <= ud+1; d++) {
// for all possible |D| samples
AutoArray<T> ds(d, EBWT_CAT);
for(T j = 0; j < d; j++) {
ds[j] = j;
}
assert(increasing(ds, d));
while(true) {
// reset diffs[]
for(T t = 1; t < v; t++) {
diffs[t] = false;
}
T diffCnt = 0;
// diffs[] records all the differences observed
for(T di = 0; di < d; di++) {
for(T dj = di+1; dj < d; dj++) {
assert_lt(ds[di], ds[dj]);
T d1 = (ds[dj] - ds[di]);
T d2 = (ds[di] + v - ds[dj]);
assert_lt(d1, v);
assert_lt(d2, v);
assert_gt(d1, 0);
assert_gt(d2, 0);
if(!diffs[d1]) diffCnt++; diffs[d1] = true;
if(!diffs[d2]) diffCnt++; diffs[d2] = true;
}
}
// Do we observe all possible differences (except 0)
ok = diffCnt == v-1;
if(ok) {
// Yes, all differences are covered
break;
} else {
// Advance ds
// (Following is commented out because it turns out
// it's slow)
// Find a missing difference
//uint32_t missing = 0xffffffff;
//for(uint32_t t = 1; t < v; t++) {
// if(diffs[t] == false) {
// missing = diffs[t];
// break;
// }
//}
//assert_neq(missing, 0xffffffff);
assert(increasing(ds, d));
bool advanced = false;
bool keepGoing = false;
do {
keepGoing = false;
for(T bd = d-1; bd > 1; bd--) {
T dif = (d-1)-bd;
if(ds[bd] < v-1-dif) {
ds[bd]++;
assert_neq(0, ds[bd]);
// Reset subsequent ones
for(T bdi = bd+1; bdi < d; bdi++) {
assert_eq(0, ds[bdi]);
ds[bdi] = ds[bdi-1]+1;
assert_gt(ds[bdi], ds[bdi-1]);
}
assert(increasing(ds, d));
// (Following is commented out because
// it turns out it's slow)
// See if the new DC has the missing value
//if(!hasDifference(ds, d, v, missing)) {
// keepGoing = true;
// break;
//}
advanced = true;
break;
} else {
ds[bd] = 0;
// keep going
}
}
} while(keepGoing);
// No solution for this |D|
if(!advanced) break;
assert(increasing(ds, d));
}
} // next sample assignment
if(ok) {
break;
}
} // next |D|
assert(ok);
cout << "Did exhaustive v=" << v << " |D|=" << d << endl;
cout << " ";
for(T i = 0; i < d; i++) {
cout << ds[i];
if(i < d-1) cout << ",";
}
cout << endl;
}
/**
* Routune for calculating the elements of clDCs up to r = 16 using the
* technique of Colbourn and Ling.
*
* See http://citeseer.ist.psu.edu/211575.html
*/
template <typename T>
void calcColbournAndLingDCs(bool verbose = false, bool sanityCheck = false) {
for(T r = 0; r < 16; r++) {
T maxv = 24*r*r + 36*r + 13; // Corollary 2.3
T numsamp = 6*r + 4;
clDCs[r].maxV = maxv;
clDCs[r].numSamples = numsamp;
memset(clDCs[r].samples, 0, 4 * 128);
T i;
// clDCs[r].samples[0] = 0;
// Fill in the 1^r part of the B series
for(i = 1; i < r+1; i++) {
clDCs[r].samples[i] = clDCs[r].samples[i-1] + 1;
}
// Fill in the (r + 1)^1 part
clDCs[r].samples[r+1] = clDCs[r].samples[r] + r + 1;
// Fill in the (2r + 1)^r part
for(i = r+2; i < r+2+r; i++) {
clDCs[r].samples[i] = clDCs[r].samples[i-1] + 2*r + 1;
}
// Fill in the (4r + 3)^(2r + 1) part
for(i = r+2+r; i < r+2+r+2*r+1; i++) {
clDCs[r].samples[i] = clDCs[r].samples[i-1] + 4*r + 3;
}
// Fill in the (2r + 2)^(r + 1) part
for(i = r+2+r+2*r+1; i < r+2+r+2*r+1+r+1; i++) {
clDCs[r].samples[i] = clDCs[r].samples[i-1] + 2*r + 2;
}
// Fill in the last 1^r part
for(i = r+2+r+2*r+1+r+1; i < r+2+r+2*r+1+r+1+r; i++) {
clDCs[r].samples[i] = clDCs[r].samples[i-1] + 1;
}
assert_eq(i, numsamp);
assert_lt(i, 128);
if(sanityCheck) {
// diffs[] records all the differences observed
AutoArray<bool> diffs(maxv, EBWT_CAT);
for(T i = 0; i < numsamp; i++) {
for(T j = i+1; j < numsamp; j++) {
T d1 = (clDCs[r].samples[j] - clDCs[r].samples[i]);
T d2 = (clDCs[r].samples[i] + maxv - clDCs[r].samples[j]);
assert_lt(d1, maxv);
assert_lt(d2, maxv);
diffs[d1] = true;
diffs[d2] = true;
}
}
// Should have observed all possible differences (except 0)
for(T i = 1; i < maxv; i++) {
if(diffs[i] == false) cout << r << ", " << i << endl;
assert(diffs[i] == true);
}
}
}
clDCs_calced = true;
}
/**
* A precalculated list of difference covers.
*/
extern uint32_t dc0to64[65][10];
/**
* Get a difference cover for the requested periodicity v.
*/
template <typename T>
static EList<T> getDiffCover(
T v,
bool verbose = false,
bool sanityCheck = false)
{
assert_gt(v, 2);
EList<T> ret;
ret.clear();
// Can we look it up in our hardcoded array?
if(v <= 64 && dc0to64[v][0] == 0xffffffff) {
if(verbose) cout << "v in hardcoded area, but hardcoded entry was all-fs" << endl;
return ret;
} else if(v <= 64) {
ret.push_back(0);
for(size_t i = 0; i < 10; i++) {
if(dc0to64[v][i] == 0) break;
ret.push_back(dc0to64[v][i]);
}
if(sanityCheck) assert(dcRepOk(v, ret));
return ret;
}
// Can we look it up in our calcColbournAndLingDCs array?
if(!clDCs_calced) {
calcColbournAndLingDCs<uint32_t>(verbose, sanityCheck);
assert(clDCs_calced);
}
for(size_t i = 0; i < 16; i++) {
if(v <= clDCs[i].maxV) {
for(size_t j = 0; j < clDCs[i].numSamples; j++) {
T s = clDCs[i].samples[j];
if(s >= v) {
s %= v;
for(size_t k = 0; k < ret.size(); k++) {
if(s == ret[k]) break;
if(s < ret[k]) {
ret.insert(s, k);
break;
}
}
} else {
ret.push_back(s % v);
}
}
if(sanityCheck) assert(dcRepOk(v, ret));
return ret;
}
}
cerr << "Error: Could not find a difference cover sample for v=" << v << endl;
throw 1;
}
/**
* Calculate and return a delta map based on the given difference cover
* and periodicity v.
*/
template <typename T>
static EList<T> getDeltaMap(T v, const EList<T>& dc) {
// Declare anchor-map-related items
EList<T> amap;
size_t amapEnts = 1;
amap.resizeExact((size_t)v);
amap.fill(0xffffffff);
amap[0] = 0;
// Print out difference cover (and optionally calculate
// anchor map)
for(size_t i = 0; i < dc.size(); i++) {
for(size_t j = i+1; j < dc.size(); j++) {
assert_gt(dc[j], dc[i]);
T diffLeft = dc[j] - dc[i];
T diffRight = dc[i] + v - dc[j];
assert_lt(diffLeft, v);
assert_lt(diffRight, v);
if(amap[diffLeft] == 0xffffffff) {
amap[diffLeft] = dc[i];
amapEnts++;
}
if(amap[diffRight] == 0xffffffff) {
amap[diffRight] = dc[j];
amapEnts++;
}
}
}
return amap;
}
/**
* Return population count (count of all bits set to 1) of i.
*/
template<typename T>
static unsigned int popCount(T i) {
unsigned int cnt = 0;
for(size_t j = 0; j < sizeof(T)*8; j++) {
if(i & 1) cnt++;
i >>= 1;
}
return cnt;
}
/**
* Calculate log-base-2 of i
*/
template<typename T>
static unsigned int myLog2(T i) {
assert_eq(1, popCount(i)); // must be power of 2
for(size_t j = 0; j < sizeof(T)*8; j++) {
if(i & 1) return (int)j;
i >>= 1;
}
assert(false);
return 0xffffffff;
}
/**
*
*/
template<typename TStr>
class DifferenceCoverSample {
public:
DifferenceCoverSample(const TStr& __text,
uint32_t __v,
bool __verbose = false,
bool __sanity = false,
ostream& __logger = cout) :
_text(__text),
_v(__v),
_verbose(__verbose),
_sanity(__sanity),
_ds(getDiffCover(_v, _verbose, _sanity)),
_dmap(getDeltaMap(_v, _ds)),
_d((uint32_t)_ds.size()),
_doffs(),
_isaPrime(),
_dInv(),
_log2v(myLog2(_v)),
_vmask(OFF_MASK << _log2v),
_logger(__logger)
{
assert_gt(_d, 0);
assert_eq(1, popCount(_v)); // must be power of 2
// Build map from d's to idx's
_dInv.resizeExact((size_t)v());
_dInv.fill(0xffffffff);
uint32_t lim = (uint32_t)_ds.size();
for(uint32_t i = 0; i < lim; i++) {
_dInv[_ds[i]] = i;
}
}
/**
* Allocate an amount of memory that simulates the peak memory
* usage of the DifferenceCoverSample with the given text and v.
* Throws bad_alloc if it's not going to fit in memory. Returns
* the approximate number of bytes the Cover takes at all times.
*/
static size_t simulateAllocs(const TStr& text, uint32_t v) {
EList<uint32_t> ds(getDiffCover(v, false /*verbose*/, false /*sanity*/));
size_t len = text.length();
size_t sPrimeSz = (len / v) * ds.size();
// sPrime, sPrimeOrder, _isaPrime all exist in memory at
// once and that's the peak
AutoArray<TIndexOffU> aa(sPrimeSz * 3 + (1024 * 1024 /*out of caution*/), EBWT_CAT);
return sPrimeSz * 4; // sPrime array
}
uint32_t v() const { return _v; }
uint32_t log2v() const { return _log2v; }
uint32_t vmask() const { return _vmask; }
uint32_t modv(TIndexOffU i) const { return (uint32_t)(i & ~_vmask); }
TIndexOffU divv(TIndexOffU i) const { return i >> _log2v; }
uint32_t d() const { return _d; }
bool verbose() const { return _verbose; }
bool sanityCheck() const { return _sanity; }
const TStr& text() const { return _text; }
const EList<uint32_t>& ds() const { return _ds; }
const EList<uint32_t>& dmap() const { return _dmap; }
ostream& log() const { return _logger; }
void build(int nthreads);
uint32_t tieBreakOff(TIndexOffU i, TIndexOffU j) const;
int64_t breakTie(TIndexOffU i, TIndexOffU j) const;
bool isCovered(TIndexOffU i) const;
TIndexOffU rank(TIndexOffU i) const;
/**
* Print out the suffix array such that every sample offset has its
* rank filled in and every non-sample offset is shown as '-'.
*/
void print(ostream& out) {
for(size_t i = 0; i < _text.length(); i++) {
if(isCovered(i)) {
out << rank(i);
} else {
out << "-";
}
if(i < _text.length()-1) {
out << ",";
}
}
out << endl;
}
private:
void doBuiltSanityCheck() const;
void buildSPrime(EList<TIndexOffU>& sPrime, size_t padding);
bool built() const {
return _isaPrime.size() > 0;
}
void verbose(const string& s) const {
if(this->verbose()) {
this->log() << s.c_str();
this->log().flush();
}
}
const TStr& _text; // text to sample
uint32_t _v; // periodicity of sample
bool _verbose; //
bool _sanity; //
EList<uint32_t> _ds; // samples: idx -> d
EList<uint32_t> _dmap; // delta map
uint32_t _d; // |D| - size of sample
EList<TIndexOffU> _doffs; // offsets into sPrime/isaPrime for each d idx
EList<TIndexOffU> _isaPrime; // ISA' array
EList<uint32_t> _dInv; // Map from d -> idx
uint32_t _log2v;
TIndexOffU _vmask;
ostream& _logger;
};
/**
* Sanity-check the difference cover by first inverting _isaPrime then
* checking that each successive suffix really is less than the next.
*/
template <typename TStr>
void DifferenceCoverSample<TStr>::doBuiltSanityCheck() const {
uint32_t v = this->v();
assert(built());
VMSG_NL(" Doing sanity check");
TIndexOffU added = 0;
EList<TIndexOffU> sorted;
sorted.resizeExact(_isaPrime.size());
sorted.fill(OFF_MASK);
for(size_t di = 0; di < this->d(); di++) {
uint32_t d = _ds[di];
size_t i = 0;
for(size_t doi = _doffs[di]; doi < _doffs[di+1]; doi++, i++) {
assert_eq(OFF_MASK, sorted[_isaPrime[doi]]);
// Maps the offset of the suffix to its rank
sorted[_isaPrime[doi]] = (TIndexOffU)(v*i + d);
added++;
}
}
assert_eq(added, _isaPrime.size());
#ifndef NDEBUG
for(size_t i = 0; i < sorted.size()-1; i++) {
assert(sstr_suf_lt(this->text(), sorted[i], this->text(), sorted[i+1], false));
}
#endif
}
/**
* Build the s' array by sampling suffixes (suffix offsets, actually)
* from t according to the difference-cover sample and pack them into
* an array of machine words in the order dictated by the "mu" mapping
* described in Burkhardt.
*
* Also builds _doffs map.
*/
template <typename TStr>
void DifferenceCoverSample<TStr>::buildSPrime(
EList<TIndexOffU>& sPrime,
size_t padding)
{
const TStr& t = this->text();
const EList<uint32_t>& ds = this->ds();
TIndexOffU tlen = (TIndexOffU)t.length();
uint32_t v = this->v();
uint32_t d = this->d();
assert_gt(v, 2);
assert_lt(d, v);
// Record where each d section should begin in sPrime
TIndexOffU tlenDivV = this->divv(tlen);
uint32_t tlenModV = this->modv(tlen);
TIndexOffU sPrimeSz = 0;
assert(_doffs.empty());
_doffs.resizeExact((size_t)d+1);
for(uint32_t di = 0; di < d; di++) {
// mu mapping
TIndexOffU sz = tlenDivV + ((ds[di] <= tlenModV) ? 1 : 0);
assert_geq(sz, 0);
_doffs[di] = sPrimeSz;
sPrimeSz += sz;
}
_doffs[d] = sPrimeSz;
#ifndef NDEBUG
if(tlenDivV > 0) {
for(size_t i = 0; i < d; i++) {
assert_gt(_doffs[i+1], _doffs[i]);
TIndexOffU diff = _doffs[i+1] - _doffs[i];
assert(diff == tlenDivV || diff == tlenDivV+1);
}
}
#endif
assert_eq(_doffs.size(), d+1);
// Size sPrime appropriately
sPrime.resizeExact((size_t)sPrimeSz + padding);
sPrime.fill(OFF_MASK);
// Slot suffixes from text into sPrime according to the mu
// mapping; where the mapping would leave a blank, insert a 0
TIndexOffU added = 0;
TIndexOffU i = 0;
for(uint64_t ti = 0; ti <= tlen; ti += v) {
for(uint32_t di = 0; di < d; di++) {
TIndexOffU tti = ti + ds[di];
if(tti > tlen) break;
TIndexOffU spi = _doffs[di] + i;
assert_lt(spi, _doffs[di+1]);
assert_leq(tti, tlen);
assert_lt(spi, sPrimeSz);
assert_eq(OFF_MASK, sPrime[spi]);
sPrime[spi] = tti; added++;
}
i++;
}
assert_eq(added, sPrimeSz);
}
/**
* Return true iff suffixes with offsets suf1 and suf2 out of host
* string 'host' are identical up to depth 'v'.
*/
template <typename TStr>
static inline bool suffixSameUpTo(
const TStr& host,
TIndexOffU suf1,
TIndexOffU suf2,
TIndexOffU v)
{
for(TIndexOffU i = 0; i < v; i++) {
bool endSuf1 = suf1+i >= host.length();
bool endSuf2 = suf2+i >= host.length();
if((endSuf1 && !endSuf2) || (!endSuf1 && endSuf2)) return false;
if(endSuf1 && endSuf2) return true;
if(host[suf1+i] != host[suf2+i]) return false;
}
return true;
}
template<typename TStr>
struct VSortingParam {
DifferenceCoverSample<TStr>* dcs;
TIndexOffU* sPrimeArr;
size_t sPrimeSz;
TIndexOffU* sPrimeOrderArr;
size_t depth;
const EList<size_t>* boundaries;
size_t* cur;
MUTEX_T* mutex;
};
template<typename TStr>
static void VSorting_worker(void *vp)
{
VSortingParam<TStr>* param = (VSortingParam<TStr>*)vp;
DifferenceCoverSample<TStr>* dcs = param->dcs;
const TStr& host = dcs->text();
const size_t hlen = host.length();
uint32_t v = dcs->v();
while(true) {
size_t cur = 0;
{
ThreadSafe ts(param->mutex, true);
cur = *(param->cur);
(*param->cur)++;
}
if(cur >= param->boundaries->size()) return;
size_t begin = (cur == 0 ? 0 : (*param->boundaries)[cur-1]);
size_t end = (*param->boundaries)[cur];
assert_leq(begin, end);
if(end - begin <= 1) continue;
mkeyQSortSuf2(
host,
hlen,
param->sPrimeArr,
param->sPrimeSz,
param->sPrimeOrderArr,
4,
begin,
end,
param->depth,
v);
}
}
/**
* Calculates a ranking of all suffixes in the sample and stores them,
* packed according to the mu mapping, in _isaPrime.
*/
template <typename TStr>
void DifferenceCoverSample<TStr>::build(int nthreads) {
// Local names for relevant types
VMSG_NL("Building DifferenceCoverSample");
// Local names for relevant data
const TStr& t = this->text();
uint32_t v = this->v();
assert_gt(v, 2);
// Build s'
EList<TIndexOffU> sPrime;
// Need to allocate 2 extra elements at the end of the sPrime and _isaPrime
// arrays. One element that's less than all others, and another that acts
// as needed padding for the Larsson-Sadakane sorting code.
size_t padding = 1;
VMSG_NL(" Building sPrime");
buildSPrime(sPrime, padding);
size_t sPrimeSz = sPrime.size() - padding;
assert_gt(sPrime.size(), padding);
assert_leq(sPrime.size(), t.length() + padding + 1);
TIndexOffU nextRank = 0;
{
VMSG_NL(" Building sPrimeOrder");
EList<TIndexOffU> sPrimeOrder;
sPrimeOrder.resizeExact(sPrimeSz);
for(TIndexOffU i = 0; i < sPrimeSz; i++) {
sPrimeOrder[i] = i;
}
// sPrime now holds suffix-offsets for DC samples.
{
Timer timer(cout, " V-Sorting samples time: ", this->verbose());
VMSG_NL(" V-Sorting samples");
// Extract backing-store array from sPrime and sPrimeOrder;
// the mkeyQSortSuf2 routine works on the array for maximum
// efficiency
TIndexOffU *sPrimeArr = (TIndexOffU*)sPrime.ptr();
assert_eq(sPrimeArr[0], sPrime[0]);
assert_eq(sPrimeArr[sPrimeSz-1], sPrime[sPrimeSz-1]);
TIndexOffU *sPrimeOrderArr = (TIndexOffU*)sPrimeOrder.ptr();
assert_eq(sPrimeOrderArr[0], sPrimeOrder[0]);
assert_eq(sPrimeOrderArr[sPrimeSz-1], sPrimeOrder[sPrimeSz-1]);
// Sort sample suffixes up to the vth character using a
// multikey quicksort. Sort time is proportional to the
// number of samples times v. It isn't quadratic.
// sPrimeOrder is passed in as a swapping partner for
// sPrimeArr, i.e., every time the multikey qsort swaps
// elements in sPrime, it swaps the same elements in
// sPrimeOrder too. This allows us to easily reconstruct
// what the sort did.
if(nthreads == 1) {
mkeyQSortSuf2(t, sPrimeArr, sPrimeSz, sPrimeOrderArr, 4,
this->verbose(), this->sanityCheck(), v);
} else {
int query_depth = 0;
int tmp_nthreads = nthreads;
while(tmp_nthreads > 0) {
query_depth++;
tmp_nthreads >>= 1;
}
EList<size_t> boundaries; // bucket boundaries for parallelization
TIndexOffU *sOrig = NULL;
if(this->sanityCheck()) {
sOrig = new TIndexOffU[sPrimeSz];
memcpy(sOrig, sPrimeArr, OFF_SIZE * sPrimeSz);
}
mkeyQSortSuf2(t, sPrimeArr, sPrimeSz, sPrimeOrderArr, 4,
this->verbose(), false, query_depth, &boundaries);
if(boundaries.size() > 0) {
AutoArray<tthread::thread*> threads(nthreads);
EList<VSortingParam<TStr> > tparams;
size_t cur = 0;
MUTEX_T mutex;
for(int tid = 0; tid < nthreads; tid++) {
// Calculate bucket sizes by doing a binary search for each
// suffix and noting where it lands
tparams.expand();
tparams.back().dcs = this;
tparams.back().sPrimeArr = sPrimeArr;
tparams.back().sPrimeSz = sPrimeSz;
tparams.back().sPrimeOrderArr = sPrimeOrderArr;
tparams.back().depth = query_depth;
tparams.back().boundaries = &boundaries;
tparams.back().cur = &cur;
tparams.back().mutex = &mutex;
threads[tid] = new tthread::thread(VSorting_worker<TStr>, (void*)&tparams.back());
}
for (int tid = 0; tid < nthreads; tid++) {
threads[tid]->join();
}
}
if(this->sanityCheck()) {
sanityCheckOrderedSufs(t, t.length(), sPrimeArr, sPrimeSz, v);
for(size_t i = 0; i < sPrimeSz; i++) {
assert_eq(sPrimeArr[i], sOrig[sPrimeOrderArr[i]]);
}
delete[] sOrig;
}
}
// Make sure sPrime and sPrimeOrder are consistent with
// their respective backing-store arrays
assert_eq(sPrimeArr[0], sPrime[0]);
assert_eq(sPrimeArr[sPrimeSz-1], sPrime[sPrimeSz-1]);
assert_eq(sPrimeOrderArr[0], sPrimeOrder[0]);
assert_eq(sPrimeOrderArr[sPrimeSz-1], sPrimeOrder[sPrimeSz-1]);
}
// Now assign the ranking implied by the sorted sPrime/sPrimeOrder
// arrays back into sPrime.
VMSG_NL(" Allocating rank array");
_isaPrime.resizeExact(sPrime.size());
ASSERT_ONLY(_isaPrime.fill(OFF_MASK));
assert_gt(_isaPrime.size(), 0);
{
Timer timer(cout, " Ranking v-sort output time: ", this->verbose());
VMSG_NL(" Ranking v-sort output");
for(size_t i = 0; i < sPrimeSz-1; i++) {
// Place the appropriate ranking
_isaPrime[sPrimeOrder[i]] = nextRank;
// If sPrime[i] and sPrime[i+1] are identical up to v, then we
// should give the next suffix the same rank
if(!suffixSameUpTo(t, sPrime[i], sPrime[i+1], v)) nextRank++;
}
_isaPrime[sPrimeOrder[sPrimeSz-1]] = nextRank; // finish off
#ifndef NDEBUG
for(size_t i = 0; i < sPrimeSz; i++) {
assert_neq(OFF_MASK, _isaPrime[i]);
assert_lt(_isaPrime[i], sPrimeSz);
}
#endif
}
// sPrimeOrder is destroyed
// All the information we need is now in _isaPrime
}
_isaPrime[_isaPrime.size()-1] = (TIndexOffU)sPrimeSz;
sPrime[sPrime.size()-1] = (TIndexOffU)sPrimeSz;
// _isaPrime[_isaPrime.size()-1] and sPrime[sPrime.size()-1] are just
// spacer for the Larsson-Sadakane routine to use
{
Timer timer(cout, " Invoking Larsson-Sadakane on ranks time: ", this->verbose());
VMSG_NL(" Invoking Larsson-Sadakane on ranks");
if(sPrime.size() >= LS_SIZE) {
cerr << "Error; sPrime array has so many elements that it can't be converted to a signed array without overflow." << endl;
throw 1;
}
LarssonSadakane<TIndexOff> ls;
ls.suffixsort(
(TIndexOff*)_isaPrime.ptr(),
(TIndexOff*)sPrime.ptr(),
(TIndexOff)sPrimeSz,
(TIndexOff)sPrime.size(),
0);
}
// chop off final character of _isaPrime
_isaPrime.resizeExact(sPrimeSz);
for(size_t i = 0; i < _isaPrime.size(); i++) {
_isaPrime[i]--;
}
#ifndef NDEBUG
for(size_t i = 0; i < sPrimeSz-1; i++) {
assert_lt(_isaPrime[i], sPrimeSz);
assert(i == 0 || _isaPrime[i] != _isaPrime[i-1]);
}
#endif
VMSG_NL(" Sanity-checking and returning");
if(this->sanityCheck()) doBuiltSanityCheck();
}
/**
* Return true iff index i within the text is covered by the difference
* cover sample. Allow i to be off the end of the text; simplifies
* logic elsewhere.
*/
template <typename TStr>
bool DifferenceCoverSample<TStr>::isCovered(TIndexOffU i) const {
assert(built());
uint32_t modi = this->modv(i);
assert_lt(modi, _dInv.size());
return _dInv[modi] != 0xffffffff;
}
/**
* Given a text offset that's covered, return its lexicographical rank
* among the sample suffixes.
*/
template <typename TStr>
TIndexOffU DifferenceCoverSample<TStr>::rank(TIndexOffU i) const {
assert(built());
assert_lt(i, this->text().length());
uint32_t imodv = this->modv(i);
assert_neq(0xffffffff, _dInv[imodv]); // must be in the sample
TIndexOffU ioff = this->divv(i);
assert_lt(ioff, _doffs[_dInv[imodv]+1] - _doffs[_dInv[imodv]]);
TIndexOffU isaIIdx = _doffs[_dInv[imodv]] + ioff;
assert_lt(isaIIdx, _isaPrime.size());
TIndexOffU isaPrimeI = _isaPrime[isaIIdx];
assert_leq(isaPrimeI, _isaPrime.size());
return isaPrimeI;
}
/**
* Return: < 0 if suffix i is lexicographically less than suffix j; > 0
* if suffix j is lexicographically greater.
*/
template <typename TStr>
int64_t DifferenceCoverSample<TStr>::breakTie(TIndexOffU i, TIndexOffU j) const {
assert(built());
assert_neq(i, j);
assert_lt(i, this->text().length());
assert_lt(j, this->text().length());
uint32_t imodv = this->modv(i);
uint32_t jmodv = this->modv(j);
assert_neq(0xffffffff, _dInv[imodv]); // must be in the sample
assert_neq(0xffffffff, _dInv[jmodv]); // must be in the sample
uint32_t dimodv = _dInv[imodv];
uint32_t djmodv = _dInv[jmodv];
TIndexOffU ioff = this->divv(i);
TIndexOffU joff = this->divv(j);
assert_lt(dimodv+1, _doffs.size());
assert_lt(djmodv+1, _doffs.size());
// assert_lt: expected (32024) < (0)
assert_lt(ioff, _doffs[dimodv+1] - _doffs[dimodv]);
assert_lt(joff, _doffs[djmodv+1] - _doffs[djmodv]);
TIndexOffU isaIIdx = _doffs[dimodv] + ioff;
TIndexOffU isaJIdx = _doffs[djmodv] + joff;
assert_lt(isaIIdx, _isaPrime.size());
assert_lt(isaJIdx, _isaPrime.size());
assert_neq(isaIIdx, isaJIdx); // ranks must be unique
TIndexOffU isaPrimeI = _isaPrime[isaIIdx];
TIndexOffU isaPrimeJ = _isaPrime[isaJIdx];
assert_neq(isaPrimeI, isaPrimeJ); // ranks must be unique
assert_leq(isaPrimeI, _isaPrime.size());
assert_leq(isaPrimeJ, _isaPrime.size());
return (int64_t)isaPrimeI - (int64_t)isaPrimeJ;
}
/**
* Given i, j, return the number of additional characters that need to
* be compared before the difference cover can break the tie.
*/
template <typename TStr>
uint32_t DifferenceCoverSample<TStr>::tieBreakOff(TIndexOffU i, TIndexOffU j) const {
const TStr& t = this->text();
const EList<uint32_t>& dmap = this->dmap();
assert(built());
// It's actually convenient to allow this, but we're permitted to
// return nonsense in that case
if(t[i] != t[j]) return 0xffffffff;
//assert_eq(t[i], t[j]); // if they're unequal, there's no tie to break
uint32_t v = this->v();
assert_neq(i, j);
assert_lt(i, t.length());
assert_lt(j, t.length());
uint32_t imod = this->modv(i);
uint32_t jmod = this->modv(j);
uint32_t diffLeft = (jmod >= imod)? (jmod - imod) : (jmod + v - imod);
uint32_t diffRight = (imod >= jmod)? (imod - jmod) : (imod + v - jmod);
assert_lt(diffLeft, dmap.size());
assert_lt(diffRight, dmap.size());
uint32_t destLeft = dmap[diffLeft]; // offset where i needs to be
uint32_t destRight = dmap[diffRight]; // offset where i needs to be
assert(isCovered(destLeft));
assert(isCovered(destLeft+diffLeft));
assert(isCovered(destRight));
assert(isCovered(destRight+diffRight));
assert_lt(destLeft, v);
assert_lt(destRight, v);
uint32_t deltaLeft = (destLeft >= imod)? (destLeft - imod) : (destLeft + v - imod);
if(deltaLeft == v) deltaLeft = 0;
uint32_t deltaRight = (destRight >= jmod)? (destRight - jmod) : (destRight + v - jmod);
if(deltaRight == v) deltaRight = 0;
assert_lt(deltaLeft, v);
assert_lt(deltaRight, v);
assert(isCovered(i+deltaLeft));
assert(isCovered(j+deltaLeft));
assert(isCovered(i+deltaRight));
assert(isCovered(j+deltaRight));
return min(deltaLeft, deltaRight);
}
#endif /*DIFF_SAMPLE_H_*/