typedef T Type;

SAtomic(T t = 0) { m_atomic.store(t); }

SAtomic(const atomic<T> &a) { m_atomic.store(a.load()); }

SAtomic(const SAtomic &other) { m_atomic.store(other.m_atomic.load()); }

SAtomic& operator=(const SAtomic &other) {

m_atomic.store(other.m_atomic.load());

return *this;

}

SAtomic& operator=(T t) {

m_atomic.store(t);

return *this;

}

bool Set(T value, T expected = 0) { return m_atomic.compare_exchange_strong(expected, value); }

operator T() const { return m_atomic.load(); }

T Load() const { return m_atomic.load(); }

atomic<T> m_atomic;

public:

void Restart() { start(); }

string Elapsed() const { return format(); }

void Stop() { stop (); }

void Resume() { resume(); }

typedef list<future<void>> ThreadsStatus;

ThreadsStatus active_threads_status;

// int total_jobs = jobs.size();

// cerr << "Remaining " << total_jobs << " jobs from " << total_jobs << endl;

//create ncores threads

for(int i = 0; i < ncores && !jobs.empty(); ++i) {

active_threads_status.push_front(async(launch::async, jobs.front()));

jobs.pop_front();

}

//for each finished thread create a new one until done

chrono::milliseconds span (1);

while(!active_threads_status.empty()) {

for(auto iloop = active_threads_status.begin(); iloop != active_threads_status.end(); ) {

auto done = iloop++;

if(done->wait_for(span) == future_status::timeout) // not ready

continue;

done->get();

active_threads_status.erase(done);

if(!jobs.empty()) {

active_threads_status.push_front(async(launch::async, jobs.front()));

jobs.pop_front();

}

// cerr << "Remaining jobs " << jobs.size()+active_threads_status.size() << " from " << total_jobs << endl;

}

}

// cerr << endl;

public:

CReadHolder(bool contains_paired) : m_total_seq(0), m_contains_paired(contains_paired) {};

// inserts read at the end

template <typename Container>

void PushBack(const Container& read) {

int shift = (m_total_seq*2)%64;

int read_len = 0;

for(auto it = read.rbegin(); it != read.rend(); ++it) { // put backward for kmer compatibility

if(shift == 0)

m_storage.push_back(0);

m_storage.back() += ((find(bin2NT.begin(), bin2NT.end(), *it) - bin2NT.begin()) << shift);

shift = (shift+2)%64;

++read_len;

}

m_read_length.push_back(read_len);

m_total_seq += read_len;

}

template <typename RandomIterator>

void PushBack(RandomIterator begin, uint32_t len) {

int shift = (m_total_seq*2)%64;

for(RandomIterator it = begin+len-1; ; --it) {

if(shift == 0)

m_storage.push_back(0);

m_storage.back() += ((find(bin2NT.begin(), bin2NT.end(), *it) - bin2NT.begin()) << shift);

shift = (shift+2)%64;

if(it == begin)

break;

}

m_read_length.push_back(len);

m_total_seq += len;

}

// insert sequence from other container

class string_iterator;

void PushBack(const string_iterator& is) {

size_t read_len = is.ReadLen();

m_read_length.push_back(read_len);

size_t destination_first_bit = 2*m_total_seq;

m_total_seq += read_len;

m_storage.resize((2*m_total_seq+63)/64);

const CReadHolder& other_holder = *is.m_readholderp;

size_t bit_from = is.m_position;

size_t bit_to = bit_from+2*read_len;

other_holder.CopyBits(bit_from, bit_to, m_storage, destination_first_bit, m_storage.size());

}

// swaps contents with other

void Swap(CReadHolder& other) {

swap(m_storage, other.m_storage);

swap(m_read_length, other.m_read_length);

swap(m_total_seq, other.m_total_seq);

}

// deletes all sequences and releases memory

void Clear() { CReadHolder(m_contains_paired).Swap(*this); }

// Total nucleotide count of the sequnce

size_t TotalSeq() const { return m_total_seq; }

// Maximal length of included sequences

size_t MaxLength() const {

if(m_read_length.empty())

return 0;

else

return *max_element(m_read_length.begin(), m_read_length.end());

}

// the number of kmers of give length that could be generated

size_t KmerNum(unsigned kmer_len) const {

size_t num = 0;

if(m_read_length.empty())

return num;

for(auto l : m_read_length) {

if(l >= kmer_len)

num += l-kmer_len+1;

}

return num;

}

// total number of sequences

size_t ReadNum() const { return m_read_length.size(); }

size_t MemoryFootprint() const { return 8*m_storage.capacity()+4*m_read_length.capacity(); } // memory in bytes

void Reserve(size_t seq, size_t num = 0) {

m_storage.reserve(seq/32+1);

if(num > 0)

m_read_length.reserve(num);

}

// shortest sequence length at xx% of total length

size_t NXX(double xx) const {

vector<uint32_t> read_length(m_read_length.begin(), m_read_length.end());

sort(read_length.begin(), read_length.end());

size_t nxx = 0;

size_t len = 0;

for(int j = (int)read_length.size()-1; j >= 0 && len < xx*m_total_seq; --j) {

nxx = read_length[j];

len += read_length[j];

}

return nxx;

}

// shortest sequence length at 50% of total length

size_t N50() const { return NXX(0.5); }

// iterator-type clas to access kmers

class kmer_iterator;

kmer_iterator kend() const { return kmer_iterator(0, *this, 2*m_total_seq); }

kmer_iterator kbegin(int kmer_len) const { return kmer_iterator(kmer_len, *this); }

class kmer_iterator {

public:

// dereference operator; returns value!

TKmer operator*() const {

TKmer kmer(m_kmer_len, 0);

uint64_t* guts = kmer.getPointer();

size_t bit_from = m_position;

size_t bit_to = bit_from+2*m_kmer_len;

m_readholderp->CopyBits(bit_from, bit_to, guts, 0, (2*m_kmer_len+63)/64);

return kmer;

}

// iterator advance

kmer_iterator& operator++() {

if(m_position == 2*(m_readholderp->m_total_seq-m_kmer_len)) {

m_position = 2*m_readholderp->m_total_seq;

return *this;

}

m_position += 2;

if(++m_position_in_read == m_readholderp->m_read_length[m_read]-m_kmer_len+1) {

m_position += 2*(m_kmer_len-1);

++m_read;

m_position_in_read = 0;

SkipShortReads();

}

return *this;

}

// doesn't check read boundaries - should be used only if landing in the SAME read

kmer_iterator& operator+=(int l) {

m_position += 2*l;

m_position_in_read += l;

return *this;

}

friend bool operator==(kmer_iterator const& li, kmer_iterator const& ri) { return li.m_position == ri.m_position && li.m_readholderp == ri.m_readholderp; }

friend bool operator!=(kmer_iterator const& li, kmer_iterator const& ri) { return li.m_position != ri.m_position || li.m_readholderp != ri.m_readholderp; }

friend class CReadHolder;

private:

kmer_iterator(int kmer_len, const CReadHolder& rholder, size_t position = 0, size_t position_in_read = 0, size_t read = 0) : m_readholderp(&rholder), m_read(read), m_position(position), m_kmer_len(kmer_len), m_position_in_read(position_in_read) {

SkipShortReads();

}

void SkipShortReads() {

while(m_position < 2*m_readholderp->m_total_seq && m_read < m_readholderp->m_read_length.size() && m_readholderp->m_read_length[m_read] < m_kmer_len)

m_position += 2*m_readholderp->m_read_length[m_read++];

}

const CReadHolder* m_readholderp;

size_t m_read; // read number

size_t m_position; // BIT num in concatenated sequence

uint32_t m_kmer_len;

uint32_t m_position_in_read; // SYMBOL in read

};

// iterator-type clas to access reads

string_iterator send() const { return string_iterator(*this, 2*m_total_seq, m_read_length.size()); }

string_iterator sbegin() const { return string_iterator(*this); }

enum {eSingle = 0, eFirstMate = 1, eSecondMate = 2};

class string_iterator {

public:

string_iterator() : m_readholderp(nullptr), m_position(0), m_read(0) {}

string operator*() const {

int read_length = m_readholderp->m_read_length[m_read];

string read;

read.reserve(read_length);

size_t position = m_position+2*(read_length-1);

for(int i = 0; i < read_length; ++i) {

read.push_back(bin2NT[(m_readholderp->m_storage[position/64] >> position%64) & 3]);

position -= 2;

}

return read;

}

// returns inversed binary sequence (not complemented)

// assumes that destination is extended properly and filled with 0s

void BSeq(int shift, uint64_t* destination) const {

size_t position = m_position+2*shift;

size_t len = 2*(ReadLen()-shift);

m_readholderp->CopyBits(position, position+len, destination, 0, (len+63)/64);

}

// returns clipped binary sequence in correct order

// assumes that destination is extended properly and filled with 0s

// left/right refer to the original sequence

void TrueBSeq(size_t left_clip, size_t right_clip, bool reverse_complement, uint64_t* destination) const {

auto Reverse = [](uint64_t& word) {

word = ((word & 0x3333333333333333) << 2) | ((word >> 2) & 0x3333333333333333); // swap adjacent pairs

word = ((word & 0x0F0F0F0F0F0F0F0F) << 4) | ((word >> 4) & 0x0F0F0F0F0F0F0F0F); // swap nibbles

word = ((word & 0x00FF00FF00FF00FF) << 8) | ((word >> 8) & 0x00FF00FF00FF00FF); // swap bytes

word = ((word & 0x0000FFFF0000FFFF) << 16) | ((word >> 16) & 0x0000FFFF0000FFFF); // swap 16 bit chunks

word = ((word & 0x00000000FFFFFFFF) << 32) | ((word >> 32) & 0x00000000FFFFFFFF); // swap 32 bit chunks

};

size_t position = m_position+2*right_clip; // sequence stored reversed

size_t len = 2*(ReadLen()-right_clip-left_clip);

size_t destination_size = (len+63)/64;

if(reverse_complement) {

m_readholderp->CopyBits(position, position+len, destination, 0, destination_size); // already reversed; not complemented

for(size_t p = 0; p < destination_size; ++p) // complement (will also convert trailing As into Ts)

destination[p] ^= 0xAAAAAAAAAAAAAAAA;

int partial_bits = len%64;

if(partial_bits > 0) // remove trailing Ts

destination[destination_size-1] &= (1ULL << partial_bits) - 1;

} else {

int shift_to_right_end = 64*destination_size-len;

m_readholderp->CopyBits(position, position+len, destination, shift_to_right_end, destination_size); // reversed and shifted to the end of the destination

for(size_t p = 0; p < destination_size/2; ++p) {

swap(destination[p], destination[destination_size-1-p]);

Reverse(destination[p]);

Reverse(destination[destination_size-1-p]);

}

if(destination_size%2)

Reverse(destination[destination_size/2]);

}

}

// returns number of equal nucleotides (2bit) from the beginning

// could be longer than actual sequence length if sequence is not multiple of 32

static size_t CommomSeqLen(const uint64_t* seq1p, const uint64_t* seq2p, size_t word_len) {

auto last = seq1p+word_len;

auto mism = mismatch(seq1p, last, seq2p);

size_t extend = 32*(mism.first-seq1p);

if(mism.first != last)

extend += (ffsll(*mism.first ^ *mism.second)-1)/2; // after ^ all matches are 0s; ffs returns 1-based position of the first bit set to 1

return extend;

}

string_iterator& operator++() {

if(m_read == m_readholderp->m_read_length.size())

return *this;

m_position += 2*m_readholderp->m_read_length[m_read++];

return *this;

}

size_t ReadLen() const { return m_readholderp->m_read_length[m_read]; }

kmer_iterator KmersForRead(int kmer_len) const {

if(kmer_len <= (int)m_readholderp->m_read_length[m_read])

return kmer_iterator(kmer_len, *m_readholderp, m_position, 0, m_read);

else

return m_readholderp->kend();

}

size_t Hash() const {

hash<const CReadHolder*> h1;

hash<size_t> h2;

return h1(m_readholderp)^h2(m_position);

}

struct SHash { size_t operator()(const string_iterator& is) const { return is.Hash(); } };

bool HasMate() const { return m_readholderp->m_contains_paired; }

int PairType() const {

if(!m_readholderp->m_contains_paired)

return eSingle;

else if(m_read%2) // odd

return eSecondMate;

else // even

return eFirstMate;

}

string_iterator GetMate() const { // undefined behavior if not paired container

if(m_read%2) // odd

return string_iterator(*m_readholderp, m_position-2*m_readholderp->m_read_length[m_read-1], m_read-1);

else // even

return string_iterator(*m_readholderp, m_position+2*m_readholderp->m_read_length[m_read], m_read+1);

}

friend bool operator==(const string_iterator& li, const string_iterator& ri) { return li.m_read == ri.m_read && li.m_readholderp == ri.m_readholderp; }

friend bool operator!=(const string_iterator& li, const string_iterator& ri) { return li.m_read != ri.m_read || li.m_readholderp != ri.m_readholderp; }

friend class CReadHolder;

private:

string_iterator(const CReadHolder& rholder, size_t position = 0, size_t read = 0) : m_readholderp(&rholder), m_position(position), m_read(read) {}

const CReadHolder* m_readholderp;

size_t m_position;

size_t m_read;

};

private:

// efficiently copies sequence to destination without converting it to string

// assumes that destination is extended properly and filled with 0; destination_size - number of 'used' 8-byte words in destination after copy

template <typename Dest>

void CopyBits(size_t bit_from, size_t bit_to, Dest& destination, size_t destination_bit_from, size_t destination_size) const {

if(bit_to <= bit_from)

return;

size_t word = bit_from/64;

size_t last_word = (bit_to-1)/64;

unsigned shift = bit_from%64;

size_t destination_word = destination_bit_from/64;

unsigned destination_shift = destination_bit_from%64;

if(shift > 0) { // first word partial

uint64_t chunk = (m_storage[word++] >> shift);

if(destination_shift > 0) { // destination word partial

destination[destination_word] += (chunk << destination_shift);

if(shift <= destination_shift) // we used all remaining destination word

++destination_word;

if(shift < destination_shift && destination_word < destination_size) // first word spills out

destination[destination_word] += (chunk >> (64-destination_shift));

} else { // desination word is not partial - it is bigger than chunk

destination[destination_word] = chunk;

}

destination_shift = (destination_shift+64-shift)%64;

}

for( ; word <= last_word; ++word, ++destination_word) {

if(destination_shift > 0) {

destination[destination_word] += (m_storage[word] << destination_shift);

if(destination_word+1 < destination_size)

destination[destination_word+1] += (m_storage[word] >> (64-destination_shift));

} else {

destination[destination_word] = m_storage[word];

}

}

int partial_bits = (destination_bit_from+bit_to-bit_from)%64;

if(partial_bits > 0) {

uint64_t mask = (1ULL << partial_bits) - 1;

destination[destination_size-1] &= mask;

}

}

vector<uint64_t> m_storage;

vector<uint32_t> m_read_length;

size_t m_total_seq;

bool m_contains_paired;

int SLOPE_LEN = 5;

int peak = min(rlimit,(int)bins.size()-SLOPE_LEN-1);

while(peak >= SLOPE_LEN) {

bool maxim = true;

for(int i = 1; i <= SLOPE_LEN && maxim; ++i)

maxim = bins[peak+i].second < bins[peak].second;

for(int i = 1; i <= SLOPE_LEN && maxim; ++i)

maxim = bins[peak-i].second < bins[peak].second;

if(maxim)

break;

--peak;

}

if(peak < SLOPE_LEN)

return -1;

int valley = 0;

for(int i = 1; i <= peak; ++i) {

if(bins[i].second < bins[valley].second)

valley = i;

}

if(valley == peak)

return -1;

for(int i = valley; i < (int)bins.size(); ++i) {

if(bins[i].second > bins[peak].second)

peak = i;

}

if(bins[valley].second < 0.7*bins[peak].second)

return valley;

else

return -1;

unsigned MIN_NUM = 100;

size_t gsize = 0;

for(auto& bin : bins) {

if(bin.second >= MIN_NUM)

gsize += bin.first*bin.second;

}

// step back over repeats and plasmids that are not likely to be more than 20 percent of the genome

int rl = 0;

size_t gs = 0;

for(auto& bin : bins) {

gs += bin.first*bin.second;

if(rl < (int)bins.size()-1)

++rl;

if(gs > 0.8*gsize)

break;

}

// find histogram portion with biggest volume and estimate genome size as number of kmers in the portion

int valley = -1;

int rlimit = rl;

size_t genome = 0;

size_t genome_vol = 0;

while(true) {

int v = FindValleyAndPeak(bins, rl);

size_t g = 0;

size_t g_vol = 0;

for(int i = max(0, v); i <= rl; ++i)

{

g_vol += (bins[i].first*bins[i].second);

g += bins[i].second;

}

if((v >= 0 && g > genome) || g_vol > genome_vol) {

valley = v;

rlimit = rl;

genome = g;

genome_vol = g_vol;

// cerr << valley << " " << rlimit << " " << genome << endl;

}

if(v < 0)

break;

rl = v;

}

return make_pair(valley, rlimit);

pair<int,int> grange = HistogramRange(bins);

if(grange.first < 0)

grange.first = 0;

size_t genome = 0;

size_t kmers = 0;

for(int i = grange.first; i <= grange.second; ++i) {

genome += bins[i].second;

kmers += bins[i].first*bins[i].second;

}

if(genome > 0)

return double(kmers)/genome;

else

return 0.;

pair<int,int> grange = HistogramRange(bins);

if(grange.first < 0)

grange.first = 0;

size_t genome = 0;

for(int i = grange.first; i <= grange.second; ++i)

genome += bins[i].second;

return genome;

// A hash with kmer as a key

// Implemented using a boost::variant of unordered_map<<LargeInt<N>,V> with maximal N = 16 which allows kmer size up to 512

public:

typedef V MappedType;

typedef TKmerMapN<V> Type;

CKmerMap(int kmer_len = 0) : m_kmer_len(kmer_len) {

if(m_kmer_len > 0)

m_container = CreateVariant<TKmerMapN<V>, TLargeIntMap, V>((m_kmer_len+31)/32);

}

size_t Size() const { return apply_visitor(container_size(), m_container); } // number of elements in the container

void Reserve(size_t rsrv) { apply_visitor(reserve(rsrv), m_container); } // reserves hash table for rsrv elements

void Clear() { apply_visitor(clear(), m_container); } // clear hash table

V& operator[] (const TKmer& kmer) {

if(m_kmer_len == 0)

throw runtime_error("Can't insert in uninitialized container");

return apply_visitor(mapper(kmer), m_container);

}

V* Find(const TKmer& kmer) { return apply_visitor(find(kmer), m_container); } // returns nullptr if not found

int KmerLen() const { return m_kmer_len; }

template <typename Prob>

void GetInfo(Prob& prob) { apply_visitor(get_info<Prob>(prob), m_container); } // scans the containier and calls prob(k, v) for each mapped element

private:

template <typename Prob>

struct get_info : public boost::static_visitor<> {

get_info(Prob& p) : prob(p) {}

template <typename T> void operator()(T& v) const {

for(auto& val : v)

prob(TKmer(val.first), val.second);

}

Prob& prob;

};

struct container_size : public boost::static_visitor<size_t> { template <typename T> size_t operator()(const T& v) const { return v.size();} };

struct clear : public boost::static_visitor<> { template <typename T> void operator()(const T& v) const { v.clear();} };

struct reserve : public boost::static_visitor<> {

reserve(size_t r) : rsrv(r) {}

template <typename T> void operator() (T& v) const { v.reserve(rsrv); }

size_t rsrv;

};

struct mapper : public boost::static_visitor<V&> {

mapper(const TKmer& k) : kmer(k) {}

template <typename T> V& operator()(T& v) const {

typedef typename T::key_type large_t;

return v[kmer.get<large_t>()];

}

const TKmer& kmer;

};

struct find : public boost::static_visitor<V*> {

find(const TKmer& k) : kmer(k) {}

template <typename T> V* operator()(T& v) const {

typedef typename T::key_type large_t;

typename T::iterator it = v.find(kmer.get<large_t>());

if(it != v.end())

return &(it->second);

else

return 0;

}

const TKmer& kmer;

};

Type m_container;

int m_kmer_len;