xsample07a.cpp
Go to the documentation of this file.
/*
Copyright(c) 2002-2020 Anatoliy Kuznetsov(anatoliy_kuznetsov at yahoo.com)
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
For more information please visit: http://bitmagic.io
*/
/** \example xsample07a.cpp
Use of bvector<> for k-mer fingerprint K should be short,
no minimizers here
*/
/*! \file xsample07a.cpp
\brief Example: Use of bvector<> for k-mer fingerprint
K should be short, no minimizers here (k < 24)
Example loads FASTA file (large multi-molecule file is expected,
builds a collection of k-mers for each molecule and runs
clusterization algorithm on the input collection using
set intersect (logical AND) as a similarity measure.
Example uses std::async for running parallel jobs.
*/
#include <assert.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <vector>
#include <list>
#include <map>
#include <algorithm>
#include <utility>
#include <memory>
#include <future>
#include <thread>
#include <mutex>
#include <atomic>
#include "bm64.h" // use 48-bit vectors
#include "bmalgo.h"
#include "bmserial.h"
#include "bmaggregator.h"
#include "bmsparsevec_compr.h"
#include "bmsparsevec_algo.h"
#include "bmrandom.h"
// BitMagic utilities for debug and timings
#include "bmdbg.h"
#include "bmtimer.h"
#include "bmundef.h" /* clear the pre-proc defines from BM */
#include "dna_finger.h"
using namespace std;
// Arguments
//
std::string ifa_name;
std::string ikd_name;
std::string ikd_counts_name;
std::string kh_name;
std::string ikd_rep_name;
std::string ikd_freq_name;
bool is_diag = false;
bool is_timing = false;
bool is_bench = false;
unsigned ik_size = 8;
unsigned parallel_jobs = 4;
unsigned f_percent = 5; // percent of k-mers we try to clear as over-represented
#include "cmd_args.h"
// Global types
//
typedef bm::bvector<> bvector_type;
typedef std::vector<char> vector_char_type;
typedef bm::dynamic_heap_matrix<unsigned, bm::bvector<>::allocator_type> distance_matrix_type;
typedef std::vector<std::unique_ptr<bvector_type> > bvector_ptr_vector_type;
typedef bvector_type::size_type bv_size_type;
typedef std::vector<std::pair<bv_size_type, bv_size_type> > bv_ranges_vector;
// Global vars
//
bm::chrono_taker<>::duration_map_type timing_map;
/// wait for any opening in a list of futures
/// used to schedule parallel tasks with CPU overbooking control
///
template<typename FV>
void wait_for_slot(FV& futures, unsigned* parallel_cnt, unsigned concurrency)
{
do
{
for (auto e = futures.begin(); e != futures.end(); ++e)
{
std::future_status status = e->wait_for(std::chrono::milliseconds(100));
if (status == std::future_status::ready)
{
(*parallel_cnt) -= 1;
futures.erase(e);
break;
}
} // for e
} while (*parallel_cnt >= concurrency);
}
/// Collection of sequences and k-mer fingerprint vectors
///
class CSequenceColl
{
public:
typedef std::vector<unsigned char> buffer_type;
public:
CSequenceColl()
{}
CSequenceColl(const CSequenceColl&) = delete;
void add_sequence(const string& acc, vector_char_type* seq_ptr)
{
m_acc.push_back(acc);
m_seqs.emplace_back(seq_ptr);
}
void set_buffer(size_t i, const buffer_type& buf)
{
unique_ptr<buffer_type> buf_ptr(new buffer_type(buf));
{
static std::mutex mtx_counts_lock;
std::lock_guard<std::mutex> guard(mtx_counts_lock);
if (m_kmer_bufs.size() <= i)
m_kmer_bufs.resize(i+1);
m_kmer_bufs[i].reset(buf_ptr.release());
}
}
void sync_buffers_size()
{
m_kmer_bufs.resize(this->size());
}
size_t size() const
{ assert(m_seqs.size() == m_acc.size()); return m_seqs.size(); }
const string& get_acc(size_t i) const { return m_acc[i]; }
const vector_char_type& get_sequence(size_t i) const { return *(m_seqs[i]); }
size_t seq_size(size_t i) const { return m_seqs[i]->size(); }
size_t total_seq_size() const
{
size_t sum = 0;
for (size_t i = 0; i < m_seqs.size(); ++i)
sum += seq_size(i);
return sum;
}
///
size_t buf_size() const { return m_kmer_bufs.size(); }
/// Get k-mer vector BLOB size
size_t get_buf_size(size_t i) const { return m_kmer_bufs[i]->size(); }
/// Get k-mer BLOB pointer
const unsigned char* get_buf(size_t i) const
{
const buffer_type* p = m_kmer_bufs[i].get();
if (!p)
return 0;
return p->data();
}
/// Deserialize group of k-mer fingerprint vectors
///
void deserialize_k_mers(bvector_ptr_vector_type& k_mers_vect,
const bm::bvector<>& bv_req,
bm::bvector<>::size_type bv_req_count) const;
private:
vector<unique_ptr<vector_char_type> > m_seqs;
vector<string> m_acc;
vector<unique_ptr<buffer_type> > m_kmer_bufs;
};
void CSequenceColl::deserialize_k_mers(bvector_ptr_vector_type& k_mers_vect,
const bm::bvector<>& bv_req,
bm::bvector<>::size_type bv_req_count) const
{
std::list<std::future<void> > futures;
k_mers_vect.resize(0);
k_mers_vect.reserve(bv_req_count);
bm::bvector<>::enumerator en(bv_req.first());
for (; en.valid(); ++en)
{
auto i_idx = *en;
bm::bvector<>* bv = new bm::bvector<>(); // k-mer fingerprint
k_mers_vect.emplace_back(bv);
const unsigned char* buf = this->get_buf(i_idx);
if (!buf)
continue;
futures.emplace_back( // async decompress
std::async(std::launch::async,
[bv, buf]()
{ BM_DECLARE_TEMP_BLOCK(tb); bm::deserialize(*bv, buf, tb); }
));
} // for en
for (auto& e : futures)
e.wait();
}
// -----------------------------------------------------------------------
//
/// Load multi-sequence FASTA
///
static
int load_FASTA(const std::string& fname, CSequenceColl& seq_coll)
{
unique_ptr<vector_char_type> seq_vect(new vector_char_type());
std::string line, acc;
std::ifstream fin(fname.c_str(), std::ios::in);
if (!fin.good())
return -1;
for (size_t i = 0; std::getline(fin, line); ++i)
{
if (line.empty())
continue;
if (line.front() == '>') // defline
{
if (!acc.empty())
{
seq_vect->shrink_to_fit();
seq_coll.add_sequence(acc, seq_vect.release());
acc.resize(0);
seq_vect.reset(new vector_char_type());
}
std::size_t pos = line.find_first_of(":");
if (pos == std::string::npos) // not found
{
acc = line;
}
else
{
acc = line.substr(1, pos-1);
}
continue;
}
for (std::string::iterator it = line.begin(); it != line.end(); ++it)
seq_vect->push_back(*it);
} // for
if (!acc.empty())
{
seq_vect->shrink_to_fit();
seq_coll.add_sequence(acc, seq_vect.release());
}
cout << "\r \r" << endl;
return 0;
}
/// save k-mer vectors to a file
static
void save_kmer_buffers(const std::string& fname, const CSequenceColl& seq_coll)
{
char magic_ch = '\t';
std::ofstream bfile (fname, std::ios::out | std::ios::binary);
if (!bfile.good())
{
std::cerr << "Cannot open file for write: " << fname << std::endl;
exit(1);
}
// save collection size
size_t sz = seq_coll.size();
bfile.write((char*)&sz, std::streamsize(sizeof(sz)));
// save the collection elements
//
for (size_t i = 0; i < sz; ++i)
{
size_t buf_size = 0;
const unsigned char* buf = seq_coll.get_buf(i);
if (!buf)
{
bfile.write((char*)&buf_size, std::streamsize(sizeof(buf_size)));
continue;
}
buf_size = seq_coll.get_buf_size(i);
bfile.write((char*)&buf_size, std::streamsize(sizeof(buf_size)));
if (buf_size)
{
bfile.write((char*)buf, std::streamsize(buf_size));
bfile.write((char*)&magic_ch, 1);
}
} // for i
}
/// Load k-mer vectors
///
static
void load_kmer_buffers(const std::string& fname, CSequenceColl& seq_coll)
{
char magic_ch = '\t';
std::ifstream bfile (fname, std::ios::in | std::ios::binary);
if (!bfile.good())
{
std::cerr << "Cannot open file for read: " << fname << std::endl;
exit(1);
}
// save collection size
size_t sz;
bfile.read((char*)&sz, std::streamsize(sizeof(sz)));
CSequenceColl::buffer_type buf;
// load the collection elements
//
for (size_t i = 0; i < sz; ++i)
{
size_t buf_size = 0;
bfile.read((char*)&buf_size, std::streamsize(sizeof(buf_size)));
if (buf_size)
{
buf.resize(buf_size);
bfile.read((char*) buf.data(), std::streamsize(buf_size));
char control_ch = 0;
bfile.read((char*)&control_ch, 1);
if (control_ch != magic_ch)
{
cerr << "Error: read failure!" << endl;
exit(1);
}
seq_coll.set_buffer(i, buf);
}
} // for i
}
inline
bool get_DNA_code(char bp, bm::id64_t& dna_code)
{
switch (bp)
{
case 'A':
dna_code = 0; // 00
break;
case 'T':
dna_code = 1; // 01
break;
case 'G':
dna_code = 2; // 10
break;
case 'C':
dna_code = 3; // 11
break;
default: // ambiguity codes are ignored (for simplicity)
return false;
}
return true;
}
/// Calculate k-mer as an unsigned long integer
///
///
/// @return true - if k-mer is "true" (not 'NNNNNN')
///
inline
bool get_kmer_code(const char* dna,
size_t pos, unsigned k_size,
bm::id64_t& k_mer)
{
// generate k-mer
//
bm::id64_t k_acc = 0;
unsigned shift = 0;
dna += pos;
for (size_t i = 0; i < k_size; ++i)
{
char bp = dna[i];
bm::id64_t dna_code;
bool valid = get_DNA_code(bp, dna_code);
if (!valid)
return false;
k_acc |= (dna_code << shift); // accumulate new code within 64-bit accum
shift += 2; // each DNA base pair needs 2-bits to store
} // for i
k_mer = k_acc;
return true;
}
/// Translate integer code to DNA letter
///
inline
char int2DNA(unsigned code)
{
static char lut[] = { 'A', 'T', 'G', 'C', 'N', 'M', '$' };
if (code < 5)
return lut[code];
assert(0);
return 'N';
}
/// Translate k-mer code into ATGC DNA string
///
/// @param dna - target string
/// @param k_mer - k-mer code
/// @param k_size -
inline
void translate_kmer(std::string& dna, bm::id64_t kmer_code, unsigned k_size)
{
dna.resize(k_size);
for (size_t i = 0; i < k_size; ++i)
{
unsigned dna_code = unsigned(kmer_code & 3);
char bp = int2DNA(dna_code);
dna[i] = bp;
kmer_code >>= 2;
} // for i
assert(!kmer_code);
}
/**
This function turns each k-mer into an integer number and encodes it
in a bit-vector (presense vector)
The natural limitation here is that integer has to be less tha 48-bits
(limitations of bm::bvector<>)
This method build a presense k-mer fingerprint vector which can be
used for Jaccard distance comparison.
@param bv - [out] - target bit-vector
@param seq_vect - [out] DNA sequence vector
@param k-size - dimention for k-mer generation
@param k_buf - sort buffer for generated k-mers
@param chunk_size - sort buffer size (number of k-mers per sort)
*/
template<typename BV>
void generate_k_mer_bvector(BV& bv,
const vector_char_type& seq_vect,
unsigned k_size,
std::vector<bm::id64_t>& k_buf,
const bm::id64_t chunk_size = 400000000
)
{
bv.clear();
bv.init(); // need to explicitly init to use bvector<>::set_bit_no_check()
if (seq_vect.empty())
return;
const char* dna_str = &seq_vect[0];
k_buf.reserve(size_t(chunk_size));
k_buf.resize(0);
{
bm::id64_t k_mer_code;
vector_char_type::size_type dna_sz = seq_vect.size()-(k_size-1);
vector_char_type::size_type pos = 0;
bool valid = false;
for (; pos < dna_sz; ++pos)
{
valid = get_kmer_code(dna_str, pos, k_size, k_mer_code);
if (valid)
{
k_buf.push_back(k_mer_code);
break;
}
} // for pos
const unsigned k_shift = (k_size-1) * 2;
if (valid)
{
for (++pos; pos < dna_sz; ++pos)
{
bm::id64_t bp_code;
valid = get_DNA_code(dna_str[pos + (k_size - 1)], bp_code);
if (!valid)
{
pos += k_size; // wind fwrd to the next BP char
for (; pos < dna_sz; ++pos) // search for the next valid k-mer
{
valid = get_kmer_code(dna_str, pos, k_size, k_mer_code);
if (valid)
{
k_buf.push_back(k_mer_code);
break;
}
}
continue;
}
// shift out the previous base pair code, OR the new arrival
k_mer_code = ((k_mer_code >> 2) | (bp_code << k_shift));
// generated k-mer codes are accumulated in buffer for sorting
k_buf.push_back(k_mer_code);
if (k_buf.size() == chunk_size) // soring check.point
{
std::sort(k_buf.begin(), k_buf.end());
if (k_buf.size())
{
bv.set(&k_buf[0], k_buf.size(), bm::BM_SORTED); // fast bulk set
k_buf.resize(0);
bv.optimize(); // periodically re-optimize to save memory
}
float pcnt = float(pos) / float(dna_sz);
pcnt *= 100;
cout << "\r" << unsigned(pcnt) << "% of " << dna_sz
<< " (" << (pos+1) <<") "
<< flush;
}
} // for pos
}
if (k_buf.size()) // add last incomplete chunk here
{
std::sort(k_buf.begin(), k_buf.end());
bv.set(&k_buf[0], k_buf.size(), bm::BM_SORTED); // fast bulk set
}
}
}
std::atomic_ullong k_mer_progress_count(0);
static
void generate_k_mers(CSequenceColl& seq_coll, unsigned k_size,
size_t from, size_t to)
{
assert(from <= to);
if (!seq_coll.size() || (from >= seq_coll.size()))
return;
std::vector<bm::id64_t> k_buf; // sort buffer
BM_DECLARE_TEMP_BLOCK(tb)
CSequenceColl::buffer_type buf;
typedef bm::bvector<>::allocator_type allocator_type;
typedef allocator_type::allocator_pool_type allocator_pool_type;
allocator_pool_type pool; // local pool for blocks
bm::bvector<> bv;
bm::bvector<>::mem_pool_guard mp_guard_bv; // memory pool reduces allocation calls to heap
mp_guard_bv.assign_if_not_set(pool, bv);
if (!to || to >= seq_coll.size())
to = seq_coll.size()-1;
bm::serializer<bm::bvector<> > bvs; // serializer object
bvs.set_bookmarks(false);
unsigned cnt = 0;
for (size_t i = from; i <= to; ++i)
{
const vector_char_type& seq_vect = seq_coll.get_sequence(i);
generate_k_mer_bvector(bv, seq_vect, k_size, k_buf);
// serialize the vector
//
typename bm::bvector<>::statistics st;
bv.optimize(tb, bm::bvector<>::opt_compress, &st);
buf.resize(st.max_serialize_mem);
size_t blob_size = bvs.serialize(bv, &buf[0], buf.size());
buf.resize(blob_size);
seq_coll.set_buffer(i, buf);
// local progress report counter is just to avoid atomic too often
++cnt;
if (cnt >= 100)
{
k_mer_progress_count.fetch_add(cnt);
cnt = 0;
}
} // for i
}
static
void generate_k_mers_parallel(CSequenceColl& seq_coll, unsigned k_size,
unsigned concurrency)
{
size_t total_seq_size = seq_coll.total_seq_size();
if (!concurrency)
concurrency = 1;
size_t batch_size = total_seq_size / concurrency;
if (!batch_size)
batch_size = total_seq_size;
std::list<std::future<void> > futures;
for (size_t from = 0; from <= seq_coll.size(); )
{
size_t to = from;
for (size_t to_pick = 0; to < seq_coll.size(); ++to)
{
to_pick += seq_coll.seq_size(to);
if (to_pick >= batch_size)
break;
} // for
futures.emplace_back(
std::async(std::launch::async,
[&seq_coll, k_size, from, to]() { generate_k_mers(seq_coll, k_size, from, to); }
));
from = to+1;
} // for from
// wait for all jobs to finish, print progress report
//
unsigned long long cnt = seq_coll.size();
for (auto& e : futures)
{
unsigned long long c_prev = 0;
while(1)
{
std::future_status status = e.wait_for(std::chrono::seconds(60));
if (status == std::future_status::ready)
break;
// progress report (entertainment)
//
unsigned long long c = k_mer_progress_count;
auto delta = c - c_prev;
c_prev = c;
auto remain_cnt = cnt - c;
auto remain_min = remain_cnt / delta;
cout << "\r" << c << ": progress per minute=" << delta;
if (remain_min < 120)
{
cout << " wait for " << remain_min << "m " << flush;
}
else
{
auto remain_h = remain_min / 60;
cout << " wait for " << remain_h << "h " << flush;
}
} // while
} // for
cout << endl;
}
// -----------------------------------------------------------------------
/// Group (clustrer) of sequences
///
class CSeqGroup
{
public:
CSeqGroup(bm::id64_t lead_id = ~0ull)
: m_lead_id(lead_id), m_bv_members(bm::BM_GAP)
{
m_bv_members.set(lead_id);
}
/// set id for the group representative
void set_lead(bm::id64_t lead_id)
{ add_member(m_lead_id = lead_id); }
/// Get lead id
bm::id64_t get_lead() const { return m_lead_id; }
/// check is cluster is non-empty
bool is_assigned() { return m_lead_id != ~0ull; }
/// add a member to the group
void add_member(bm::id64_t id) { m_bv_members.set_bit_no_check(id); }
void add_member(bm::id64_t id, const bm::bvector<>& bv_kmer)
{
m_bv_members.set_bit_no_check(id);
m_bv_kmer_union |= bv_kmer;
}
void add_member_sync(bm::id64_t id, const bm::bvector<>& bv_kmer)
{
std::lock_guard<std::mutex> guard(mtx_add_member_lock);
m_bv_members.set_bit_no_check(id); // a bit faster than set()
m_bv_kmer_union |= bv_kmer;
}
void merge_member_sync(bm::bvector<>& bv_seq, bm::bvector<>& bv_kmer)
{
std::lock_guard<std::mutex> guard(mtx_add_member_lock);
m_bv_members.merge(bv_seq);
m_bv_kmer_union.merge(bv_kmer);
}
bm::id64_t count_and_union_sync(const bm::bvector<>& bv)
{
std::lock_guard<std::mutex> guard(mtx_add_member_lock);
return bm::count_and(bv, m_bv_kmer_union);
}
void clear_member(bm::id64_t id) { m_bv_members.set(id, false); }
bm::bvector<>& get_rep() { return m_bv_rep; }
const bm::bvector<>& get_rep() const { return m_bv_rep; }
const bm::bvector<>& get_members() const { return m_bv_members; }
bm::bvector<>& get_members() { return m_bv_members; }
bm::bvector<>& get_kmer_union() { return m_bv_kmer_union; }
const bm::bvector<>& get_kmer_union() const { return m_bv_kmer_union; }
friend class CSeqClusters;
private:
bm::id64_t m_lead_id; ///< groups lead vector ID
bm::bvector<> m_bv_members; ///< vector of all group member ids
std::mutex mtx_add_member_lock; ///< concurrency protection for add_member_sync()
bm::bvector<> m_bv_rep; ///< group's representative vector
bm::bvector<> m_bv_kmer_union; ///< union of all k-mers of members
};
// -----------------------------------------------------------------------
//
//
class CSeqClusters
{
public:
typedef std::vector<std::unique_ptr<CSeqGroup> > groups_vector_type;
public:
CSeqClusters()
{}
CSeqClusters(const CSeqClusters&)=delete;
void add_group(CSeqGroup* sg) { m_seq_groups.emplace_back(sg); }
/// memebers moved into their own group
void take_group(bm::bvector<>& bv_members);
/// Acquire all groups from another cluster collection
///
void merge_from(CSeqClusters& sc);
/// Remove groups which turned empty after clusterization
void clear_empty_groups();
/// Compute union of all cluster group members
const bm::bvector<>& union_all_groups();
/// Resolve duplicate membership between groups
///
void resolve_duplicates(const CSequenceColl& seq_coll);
/// Find the best representatives in all cluster groups
/// the criteria is maximum absolute similarity to all members
///
void elect_leaders(const CSequenceColl& seq_coll,
unsigned concurrency);
/// calculate avg cluster population count
bm::id64_t compute_avg_count() const;
size_t groups_size() const { return m_seq_groups.size(); }
CSeqGroup* get_group(size_t idx) { return m_seq_groups[idx].get(); }
/// print clusterization report
void print_summary(const char* title) const;
private:
bm::bvector<> m_all_members; ///< Union of all group members
groups_vector_type m_seq_groups; ///< vector of all formed clusters
bm::aggregator<bvector_type> agg; ///< fast aggregator for set UNION ops
};
void CSeqClusters::clear_empty_groups()
{
for (groups_vector_type::iterator it = m_seq_groups.begin();
it != m_seq_groups.end(); )
{
CSeqGroup* sg = it->get();
bm::bvector<>& bv_mem = sg->get_members();
auto cnt = bv_mem.count();
if (cnt < 2) // practically empty group
it = m_seq_groups.erase(it);
else
++it;
} // for
}
void CSeqClusters::take_group(bm::bvector<>& bv_members)
{
bm::id64_t lead_id = bv_members.get_first();
CSeqGroup* sg = new CSeqGroup(lead_id);
sg->get_members().swap(bv_members); // move members to the cluster
add_group(sg);
}
const bm::bvector<>& CSeqClusters::union_all_groups()
{
agg.set_optimization();
m_all_members.clear();
for (groups_vector_type::const_iterator it = m_seq_groups.begin();
it != m_seq_groups.end(); ++it)
{
const CSeqGroup* sg = it->get();
agg.add(&sg->get_members());
} // for
agg.combine_or(m_all_members); // run UNION of all member vectors
agg.reset();
return m_all_members;
}
bm::id64_t CSeqClusters::compute_avg_count() const
{
bm::id64_t cnt = 0;
for (groups_vector_type::const_iterator it = m_seq_groups.begin();
it != m_seq_groups.end(); ++it)
{
const CSeqGroup* sg = it->get();
cnt += sg->get_members().count();
}
cnt = cnt / m_seq_groups.size();
return cnt;
}
void CSeqClusters::merge_from(CSeqClusters& sc)
{
for (auto it = sc.m_seq_groups.begin(); it != sc.m_seq_groups.end(); ++it)
m_seq_groups.emplace_back(it->release());
sc.m_seq_groups.clear();
}
void resolve_duplicates(CSeqGroup& seq_group1,
CSeqGroup& seq_group2,
const CSequenceColl& seq_coll);
void CSeqClusters::resolve_duplicates(const CSequenceColl& seq_coll)
{
for (size_t i = 0; i < m_seq_groups.size(); ++i)
{
CSeqGroup* sg1 = m_seq_groups[i].get();
for (size_t j = 0; j < i; ++j)
{
CSeqGroup* sg2 = m_seq_groups[j].get();
// resolve pairwise conflicts between cluster groups
::resolve_duplicates(*sg1, *sg2, seq_coll);
} // for j
} // for i
}
/// Compute similarity distances for one row/vector (1:N) of distance matrix
///
static
void compute_and_sim_row(
unsigned* row,
const bm::bvector<>* bv_i,
size_t i,
const std::vector<std::unique_ptr<bvector_type> >& k_mers_vect)
{
size_t j;
for (j = 0; j < i; ++j)
{
const bm::bvector<>* bv_j = k_mers_vect[j].get();
bm::id64_t and_cnt = bm::count_and(*bv_i, *bv_j);
row[j] = unsigned(and_cnt);
}
auto cnt = bv_i->count();
row[j] = unsigned(cnt);
}
/// Compute similarity distances matrix (COUNT(AND(a, b))
///
static
void compute_and_sim(distance_matrix_type& dm,
const CSequenceColl& seq_coll,
const bm::bvector<>& bv_mem,
bm::bvector<>::size_type bv_mem_cnt,
unsigned concurrency)
{
if (concurrency < 1)
concurrency = 1;
auto N = bv_mem_cnt;
const unsigned k_max_electors = 500;
bm::bvector<> bv_sub; // subset vector
if (N > k_max_electors) // TODO: parameterize the
{
bm::random_subset<bm::bvector<> > rsub;
rsub.sample(bv_sub, bv_mem, k_max_electors); // pick random sunset
}
bvector_ptr_vector_type k_mers_vect;
// materialize list of vectors used in distance calculation
//
seq_coll.deserialize_k_mers(k_mers_vect, bv_mem, N);
std::list<std::future<void> > futures;
unsigned parallel_cnt = 0; // number of jobs in flight at a time
size_t i;
for (i = 0; i < N; ++i)
{
const bm::bvector<>* bv_i = k_mers_vect[i].get();
unsigned* row = dm.row(i);
do
{
if (parallel_cnt < concurrency)
{
futures.emplace_back(
std::async(std::launch::async,
[row, bv_i, i, &k_mers_vect]()
{ compute_and_sim_row(row, bv_i, i, k_mers_vect); }
));
++parallel_cnt;
break;
}
else
{
// wait for an async() slot to open (overbooking control)
::wait_for_slot(futures, &parallel_cnt, concurrency);
}
} while(1);
} // for i
// sync point
for (auto& e : futures)
e.wait();
dm.replicate_triange(); // copy to full simetrical distances from triangle
}
/// Compute union (Universe) of all k-mers in the cluster group
/// Implemented as a OR of all k-mer fingerprints
static
void compute_seq_group_union(CSeqGroup& seq_group,
const CSequenceColl& seq_coll)
{
BM_DECLARE_TEMP_BLOCK(tb)
bm::bvector<>& bv_kmer_union = seq_group.get_kmer_union();
bv_kmer_union.clear();
bm::bvector<>& bv_all_members = seq_group.get_members();
for(bm::bvector<>::enumerator en(bv_all_members) ;en.valid(); ++en)
{
auto idx = *en;
const unsigned char* buf = seq_coll.get_buf(idx);
if (!buf)
continue;
bm::deserialize(bv_kmer_union, buf, tb); // deserialize is an OR operation
} // for en
bv_kmer_union.optimize(tb);
}
void CSeqClusters::elect_leaders(const CSequenceColl& seq_coll,
unsigned concurrency)
{
bm::operation_deserializer<bm::bvector<> > od;
bm::random_subset<bm::bvector<> > rsub;
std::list<std::future<void> > futures;
for (size_t k = 0; k < m_seq_groups.size(); ++k)
{
CSeqGroup* sg = m_seq_groups[k].get();
bm::bvector<>& bv_all_members = sg->get_members();
bv_all_members.set(sg->get_lead());
auto N = bv_all_members.count();
auto all_members_count = N; (void) all_members_count;
// determine the size of the "electoral colledge" on available concurrency
unsigned k_max_electors = 200 * unsigned(log2(concurrency));
if (k_max_electors < 500)
k_max_electors = 500;
const bm::bvector<>* bv_mem = &bv_all_members; // vector of electors
bm::bvector<> bv_sub; // subset vector
if (N > k_max_electors) // TODO: parameterize the
{
// pick a random sub-set as the "electoral colledge"
rsub.sample(bv_sub, bv_all_members, k_max_electors);
bv_sub.set(sg->get_lead()); // current leader always takes part
bv_mem = &bv_sub;
N = bv_sub.count();
}
// NxN distance matrix between members
//
distance_matrix_type dm(N, N);
dm.init();
dm.set_zero();
// compute triangular distance matrix
//
compute_and_sim(dm, seq_coll, *bv_mem, N, concurrency);
// leader election based on maximum similarity to other cluster
// elements
//
bm::id64_t best_score = 0;
bm::id64_t leader_idx = sg->get_lead();
assert(bv_mem->test(leader_idx));
auto rank = bv_mem->count_range(0, leader_idx);
assert(rank);
dm.sum(best_score, rank-1); // use sum or all similarities as a score for the leader
bm::id64_t old_leader_idx = leader_idx;
bm::bvector<>::enumerator en(bv_mem->first());
for (size_t i = 0; en.valid(); ++en, ++i)
{
bm::id64_t cand_score;
dm.sum(cand_score, i); // use sum or all similarities as a score
if (cand_score > best_score) // better candidate for a leader
{
best_score = cand_score;
leader_idx = *en;
}
} // for en
if (leader_idx != old_leader_idx) // found a new leader
{
sg->set_lead(leader_idx);
const unsigned char* buf = seq_coll.get_buf(leader_idx);
assert(buf);
bm::bvector<>* bv = &sg->get_rep();
// async replace the leader representative k-mer vector
//
futures.emplace_back(
std::async(std::launch::async,
[bv, buf]() { bv->clear(); bm::deserialize(*bv, buf); }
));
// re-compute k-mer UNION of all vectors
futures.emplace_back(
std::async(std::launch::async,
[sg, &seq_coll]() { compute_seq_group_union(*sg, seq_coll); }
));
}
} // for i
for (auto& e : futures) // wait for all forked tasks
e.wait();
}
void CSeqClusters::print_summary(const char* title) const
{
cout << title << endl;
for (size_t i = 0; i < m_seq_groups.size(); ++i)
{
const CSeqGroup* sg = m_seq_groups[i].get();
const bm::bvector<>& bv_mem = sg->get_members();
cout << sg->get_lead() << ": "
<< bv_mem.count() << endl;
}
cout << "-----------\nTotal: " << m_seq_groups.size() << endl << endl;
}
///
static
void compute_group(CSeqGroup& seq_group,
const CSequenceColl& seq_coll,
const bm::bvector<>& bv_exceptions,
float similarity_cut_off)
{
assert(similarity_cut_off < 1);
assert(seq_group.is_assigned());
auto sz = seq_coll.buf_size();
size_t lead_id = seq_group.get_lead();
if (lead_id >= sz)
return;
const unsigned char* buf = seq_coll.get_buf(lead_id);
if (!buf)
return;
bm::bvector<>& bv = seq_group.get_rep();
bm::deserialize(bv, buf);
auto i_cnt = bv.count();
// approximate number of k-mers we consider similar
float similarity_target = float(i_cnt * float(similarity_cut_off));
bm::operation_deserializer<bm::bvector<> > od;
bool found = false;
for (size_t i = 0; i < sz; ++i)
{
bool is_except = bv_exceptions.test(i);
if (is_except)
continue;
buf = seq_coll.get_buf(i);
if (!buf)
continue;
// constant deserializer AND just to count the product
// without actual deserialization (from the compressed BLOB)
//
bm::id64_t and_cnt = od.deserialize(bv, buf, 0, bm::set_COUNT_AND);
if (and_cnt && (float(and_cnt) > similarity_target)) // similar enough to be a candidate
{
seq_group.add_member(i);
found = true;
}
} // for i
if (!found)
seq_group.clear_member(lead_id);
}
/// Resolve duplicate members between two groups
void resolve_duplicates(CSeqGroup& seq_group1,
CSeqGroup& seq_group2,
const CSequenceColl& seq_coll)
{
if (&seq_group1 == & seq_group2) // self check
return;
bm::bvector<>& bv1 = seq_group1.get_members();
bm::bvector<>& bv2 = seq_group2.get_members();
// build intersect between group members
bm::bvector<> bv_and;
bv_and.bit_and(bv1, bv2, bm::bvector<>::opt_none);
if (bv_and.any()) // double membership detected
{
bm::bvector<>& bv_rep1 = seq_group1.get_rep();
bm::bvector<>& bv_rep2 = seq_group2.get_rep();
bm::operation_deserializer<bm::bvector<> > od;
// evaluate each double member for best membership placement
//
for (bm::bvector<>::enumerator en(bv_and); en.valid(); ++en)
{
auto idx = *en;
auto lead_idx1 = seq_group1.get_lead();
auto lead_idx2 = seq_group2.get_lead();
assert(lead_idx1 != lead_idx2);
if (idx == lead_idx1)
{
seq_group2.clear_member(idx);
continue;
}
if (idx == lead_idx2)
{
seq_group1.clear_member(idx);
continue;
}
const unsigned char* buf = seq_coll.get_buf(idx);
assert(buf);
// resolve conflict based on max.absolute similarity
//
bm::id64_t and_cnt1 = od.deserialize(bv_rep1, buf, 0, bm::set_COUNT_AND);
bm::id64_t and_cnt2 = od.deserialize(bv_rep2, buf, 0, bm::set_COUNT_AND);
if (and_cnt1 >= and_cnt2)
seq_group2.clear_member(idx);
else
seq_group1.clear_member(idx);
} // for
}
}
/// Utility class to accumulate cahnges to cluster
/// before commiting it (mutex syncronous operation)
///
struct CKMerAcc
{
CKMerAcc(size_t sz)
: bv_members(sz), bv_kmers(sz)
{}
void add(size_t cluster_id,
bm::bvector<>::size_type m_id, bm::bvector<>& bv_kmer)
{
bm::bvector<>* bv_m = bv_members[cluster_id].get();
bm::bvector<>* bv_k = bv_kmers[cluster_id].get();
if (!bv_m)
{
assert(!bv_kmers[cluster_id].get());
bv_members[cluster_id].reset(new bm::bvector<>(bm::BM_GAP));
bv_kmers[cluster_id].reset(new bm::bvector<>(bm::BM_GAP));
bv_m = bv_members[cluster_id].get();
bv_k = bv_kmers[cluster_id].get();
}
bv_m->set(m_id);
bv_k->merge(bv_kmer);
}
bvector_ptr_vector_type bv_members;
bvector_ptr_vector_type bv_kmers;
};
/// Compute AND similarity to all available clusters assign to the most similar
/// using cluster representative
///
static
void assign_to_best_cluster(CSeqClusters& cluster_groups,
const CSequenceColl& seq_coll,
const bm::bvector<>& bv_seq_ids,
bm::bvector<>::size_type seq_id_from,
bm::bvector<>::size_type seq_id_to)
{
BM_DECLARE_TEMP_BLOCK(tb)
bm::bvector<>::allocator_pool_type pool;
CKMerAcc acc(cluster_groups.groups_size());
bm::bvector<>::enumerator en(bv_seq_ids);
en.go_to(seq_id_from);
for ( ;en.valid(); ++en)
{
auto seq_id = *en;
if (seq_id > seq_id_to)
break;
const unsigned char* buf = seq_coll.get_buf(seq_id);
if (!buf)
continue;
bm::bvector<> bv_k_mer;
bv_k_mer.set_allocator_pool(&pool); // for faster memory recycling
bm::deserialize(bv_k_mer, buf, tb);
bm::bvector<>::size_type best_score(0);
bm::bvector<>::size_type cluster_idx(~0ull);
// analyse candidate's similarity to all clusters via representative
//
for (size_t i = 0; i < cluster_groups.groups_size(); ++i)
{
CSeqGroup* sg = cluster_groups.get_group(i);
// - COUNT(AND) similarity to the representative of each cluster
//
bm::bvector<>& bv_rep_k_mer = sg->get_rep();
auto rep_and_cnt = bm::count_and(bv_k_mer, bv_rep_k_mer);
if (rep_and_cnt > best_score)
{
cluster_idx = i; best_score = rep_and_cnt;
}
} // for i
if (cluster_idx != ~0ull) // cluster association via representative
{
acc.add(cluster_idx, seq_id, bv_k_mer);
}
} // for all seq-ids in the range
// merge all accumulated cluster assignmnets all at once
//
for (size_t i = 0; i < cluster_groups.groups_size(); ++i)
{
bm::bvector<>* bv_m = acc.bv_members[i].get();
if (!bv_m)
continue;
bm::bvector<>* bv_k = acc.bv_kmers[i].get();
CSeqGroup* sg = cluster_groups.get_group(i);
sg->merge_member_sync(*bv_m, *bv_k);
} // for
}
/// Compute AND similarity to all available clusters assign to the most similar
/// using UNION of k-mers in the cluster
/// This is a more relaxed assignmnet, used when representative does not work
static
void assign_to_best_cluster_union(CSeqClusters& cluster_groups,
const CSequenceColl& seq_coll,
const bm::bvector<>& bv_seq_ids,
bm::bvector<>::size_type seq_id_from,
bm::bvector<>::size_type seq_id_to)
{
BM_DECLARE_TEMP_BLOCK(tb)
bm::bvector<>::allocator_pool_type pool;
bm::bvector<>::enumerator en(bv_seq_ids);
en.go_to(seq_id_from);
for ( ;en.valid(); ++en)
{
auto seq_id = *en;
if (seq_id > seq_id_to)
break;
const unsigned char* buf = seq_coll.get_buf(seq_id);
if (!buf)
continue;
bm::bvector<> bv_k_mer;
bv_k_mer.set_allocator_pool(&pool); // for faster memory recycling
bm::deserialize(bv_k_mer, buf, tb);
bm::bvector<>::size_type best_score(0);
bm::bvector<>::size_type cluster_idx(~0ull);
{
// try to extend search using UNION of all cluster k-mers
//
best_score = 0;
for (size_t i = 0; i < cluster_groups.groups_size(); ++i)
{
CSeqGroup* sg = cluster_groups.get_group(i);
// - COUNT(AND) similarity to the representative of each cluster
//
bm::id64_t uni_and_cnt = sg->count_and_union_sync(bv_k_mer);
if (uni_and_cnt > best_score)
{
cluster_idx = i; best_score = uni_and_cnt;
}
} // for i
if (cluster_idx != ~0ull) // cluster association via representative
{
CSeqGroup* sg = cluster_groups.get_group(cluster_idx);
sg->add_member_sync(seq_id, bv_k_mer);
}
}
} // for all seq-ids in the range
}
/// Pick random sequences as cluster seed elements, try attach
/// initial sequences based on weighted similarity measure
///
static
void compute_random_clusters(CSeqClusters& cluster_groups,
const CSequenceColl& seq_coll,
const bm::bvector<>& bv_total,
bm::random_subset<bvector_type>& rsub,
unsigned num_clust,
float similarity_cut_off,
unsigned concurrency)
{
bm::bvector<> bv_rsub; // random subset of sequences
rsub.sample(bv_rsub, bv_total, num_clust); // pick random sequences as seeds
std::list<std::future<void> > futures;
unsigned parallel_cnt = 0;
for (bm::bvector<>::enumerator en=bv_rsub.first(); en.valid(); ++en)
{
auto idx = *en;
CSeqGroup* sg = new CSeqGroup(idx);
cluster_groups.add_group(sg);
do
{
if (parallel_cnt < concurrency)
{
futures.emplace_back(
std::async(std::launch::async,
[&seq_coll, sg, &bv_rsub, similarity_cut_off]()
{ compute_group(*sg, seq_coll, bv_rsub, similarity_cut_off); }
));
++parallel_cnt;
break;
}
else
{
// wait for an async() slot to open
::wait_for_slot(futures, &parallel_cnt, concurrency);
}
} while(1);
} // for en
// wait for completion of initial cluster group formation
for (auto& e : futures)
e.wait();
}
static
void compute_jaccard_clusters(CSeqClusters& seq_clusters,
const CSequenceColl& seq_coll,
unsigned num_clust,
float similarity_cut_off,
unsigned concurrency)
{
assert(similarity_cut_off < 1);
if (!seq_coll.buf_size())
return; // nothing to do
bm::random_subset<bm::bvector<> > rsub; // sub-set getter
bm::bvector<> bv_total;
bv_total.set_range(0, seq_coll.buf_size());
bm::id64_t rcount = 0;
const unsigned max_pass = 3;
for (unsigned pass = 0; pass < max_pass; ++pass)
{
CSeqClusters cluster_groups;
compute_random_clusters(cluster_groups, seq_coll, bv_total, rsub,
num_clust, similarity_cut_off, concurrency);
// remove possible empty clusters (inital seeds were picked at random)
//
cluster_groups.clear_empty_groups();
cluster_groups.resolve_duplicates(seq_coll);
cluster_groups.clear_empty_groups();
// print summary after the initial formation of cluster groups
//
cluster_groups.print_summary("Inital cluster formations:");
// re-elect representatives
cluster_groups.elect_leaders(seq_coll, concurrency);
cluster_groups.print_summary("After lead re-election:");
// sub-set of sequence ids already distributed into clusters
bm::bvector<>::size_type total_count = bv_total.count();
{
cout << " total = " << total_count << endl;
const bm::bvector<>& bv_clust = cluster_groups.union_all_groups();
cout << " clustered = " << bv_clust.count() << endl;
bv_total -= bv_clust; // exlude all already clustered
total_count = bv_total.count();
cout << " remain = " << total_count << endl;
}
if (!total_count)
break;
std::list<std::future<void> > futures;
// run split algorithm to determine approximately equal ranges
// for parallel processing
bv_ranges_vector pair_vect;
bm::bvector<>::size_type split_rank = total_count / concurrency; // target population count per job
bm::rank_range_split(bv_total, split_rank, pair_vect);
assert(pair_vect.size());
for (size_t k = 0; k < pair_vect.size(); ++k)
{
auto seq_id_from = pair_vect[k].first;
auto seq_id_to = pair_vect[k].second;
futures.emplace_back(
std::async(std::launch::async,
[&cluster_groups, &seq_coll, &bv_total, seq_id_from, seq_id_to]()
{ assign_to_best_cluster(cluster_groups, seq_coll, bv_total, seq_id_from, seq_id_to); }
));
}
for (auto& e : futures) // sync point
e.wait();
cluster_groups.print_summary("Clusters after phase 2 recruitment");
// check if there are sequences not yet belonging to any cluster
{
const bm::bvector<>& bv_clust = cluster_groups.union_all_groups();
bv_total -= bv_clust; // exlude all already clustered
rcount = bv_total.count();
if (rcount)
{
cout << "Undistributed sequences = " << rcount << endl;
}
else
{
seq_clusters.merge_from(cluster_groups);
break;
}
}
seq_clusters.merge_from(cluster_groups);
bm::id64_t avg_group_count = seq_clusters.compute_avg_count();
if (rcount < avg_group_count) // not worth another pass ?
{
seq_clusters.take_group(bv_total);
break;
}
cout << "PASS=" << (pass+1) << endl << endl;
} // for pass
// try to assign to the global pool of clusters using UNION
// which relaxes assignmnet
if (rcount)
{
{
const bm::bvector<>& bv_clust = seq_clusters.union_all_groups();
cout << endl << " clustered = " << bv_clust.count() << endl;
bv_total -= bv_clust; // exlude all already clustered
rcount = bv_total.count();
cout << " remain = " << rcount << endl;
}
if (rcount)
{
bv_ranges_vector pair_vect;
bm::bvector<>::size_type split_rank = rcount / concurrency;
bm::rank_range_split(bv_total, split_rank, pair_vect);
std::list<std::future<void> > futures;
for (size_t k = 0; k < pair_vect.size(); ++k)
{
auto seq_id_from = pair_vect[k].first;
auto seq_id_to = pair_vect[k].second;
futures.emplace_back(
std::async(std::launch::async,
[&seq_clusters, &seq_coll, &bv_total, seq_id_from, seq_id_to]()
{ assign_to_best_cluster_union(seq_clusters, seq_coll, bv_total, seq_id_from, seq_id_to); }
));
}
for (auto& e : futures) // sync point
e.wait();
{
const bm::bvector<>& bv_clust = seq_clusters.union_all_groups();
cout << endl << " clustered = " << bv_clust.count() << endl;
bv_total -= bv_clust; // exlude all already clustered
rcount = bv_total.count();
cout << " remain = " << rcount << endl;
}
}
}
if (rcount)
{
seq_clusters.take_group(bv_total);
}
seq_clusters.clear_empty_groups();
seq_clusters.resolve_duplicates(seq_coll);
seq_clusters.clear_empty_groups();
seq_clusters.print_summary("Final clusters summary:");
}
int main(int argc, char *argv[])
{
CSequenceColl seq_coll;
try
{
auto ret = parse_args(argc, argv);
if (ret != 0)
{
cerr << "cmd-line parse error. " << endl;
return ret;
}
if (!ifa_name.empty()) // FASTA file load
{
bm::chrono_taker tt1(cout, "1. Load FASTA", 1, &timing_map);
// limitation: loads a single molecule only
//
auto res = load_FASTA(ifa_name, seq_coll);
if (res != 0)
return res;
}
cout << "Sequences size = " << seq_coll.size() << endl;
if (ik_size && !ifa_name.empty())
{
{
bm::chrono_taker tt1(cout, "2. Generate k-mers", 1, &timing_map);
seq_coll.sync_buffers_size();
generate_k_mers_parallel(seq_coll, ik_size, parallel_jobs);
}
if (!ikd_name.empty())
{
bm::chrono_taker tt1(cout, "3. Save k-mers", 1, &timing_map);
save_kmer_buffers(ikd_name, seq_coll);
}
}
if (ik_size && ifa_name.empty() && !ikd_name.empty())
{
{
bm::chrono_taker tt1(cout, "4. Load k-mers", 1, &timing_map);
load_kmer_buffers(ikd_name, seq_coll);
}
if (seq_coll.buf_size())
{
CSeqClusters seq_clusters;
bm::chrono_taker tt1(cout, "5. k-mer similarity clustering", 1, &timing_map);
compute_jaccard_clusters(seq_clusters, seq_coll, 10, 0.05f, parallel_jobs);
}
}
if (is_timing)
{
std::cout << std::endl << "Performance:" << std::endl;
bm::chrono_taker<>::print_duration_map(cout, timing_map, bm::chrono_taker<>::ct_time);
}
}
catch(std::exception& ex)
{
std::cerr << ex.what() << std::endl;
return 1;
}
return 0;
}
BM_DECLARE_TEMP_BLOCK
#define BM_DECLARE_TEMP_BLOCK(x)
Definition bm.h:47
bmaggregator.h
Algorithms for fast aggregation of N bvectors.
bmalgo.h
Algorithms for bvector<> (main include).
bmrandom.h
Generation of random subset.
bmserial.h
Serialization / compression of bvector<>. Set theoretical operations on compressed BLOBs.
bmsparsevec_algo.h
Algorithms for bm::sparse_vector.
bmsparsevec_compr.h
Compressed sparse container rsc_sparse_vector<> for integer types.
bmtimer.h
Timing utilities for benchmarking (internal).
bmundef.h
pre-processor un-defines to avoid global space pollution (internal)
CSeqClusters
Definition xsample07a.cpp:758
CSeqClusters::merge_from
void merge_from(CSeqClusters &sc)
Acquire all groups from another cluster collection.
Definition xsample07a.cpp:861
CSeqClusters::groups_size
size_t groups_size() const
Definition xsample07a.cpp:795
CSeqClusters::elect_leaders
void elect_leaders(const CSequenceColl &seq_coll, unsigned concurrency)
Find the best representatives in all cluster groups the criteria is maximum absolute similarity to al...
Definition xsample07a.cpp:994
CSeqClusters::resolve_duplicates
void resolve_duplicates(const CSequenceColl &seq_coll)
Resolve duplicate membership between groups.
Definition xsample07a.cpp:873
CSeqClusters::CSeqClusters
CSeqClusters(const CSeqClusters &)=delete
CSeqClusters::clear_empty_groups
void clear_empty_groups()
Remove groups which turned empty after clusterization.
Definition xsample07a.cpp:807
CSeqClusters::print_summary
void print_summary(const char *title) const
print clusterization report
Definition xsample07a.cpp:1087
CSeqClusters::add_group
void add_group(CSeqGroup *sg)
Definition xsample07a.cpp:766
CSeqClusters::groups_vector_type
std::vector< std::unique_ptr< CSeqGroup > > groups_vector_type
Definition xsample07a.cpp:760
CSeqClusters::get_group
CSeqGroup * get_group(size_t idx)
Definition xsample07a.cpp:796
CSeqClusters::union_all_groups
const bm::bvector & union_all_groups()
Compute union of all cluster group members.
Definition xsample07a.cpp:830
CSeqClusters::compute_avg_count
bm::id64_t compute_avg_count() const
calculate avg cluster population count
Definition xsample07a.cpp:847
CSeqClusters::take_group
void take_group(bm::bvector<> &bv_members)
memebers moved into their own group
Definition xsample07a.cpp:822
CSeqClusters::CSeqClusters
CSeqClusters()
Definition xsample07a.cpp:762
CSeqGroup
Group (clustrer) of sequences.
Definition xsample07a.cpp:687
CSeqGroup::add_member
void add_member(bm::id64_t id, const bm::bvector<> &bv_kmer)
Definition xsample07a.cpp:706
CSeqGroup::get_rep
bm::bvector & get_rep()
Definition xsample07a.cpp:733
CSeqGroup::get_kmer_union
const bm::bvector & get_kmer_union() const
Definition xsample07a.cpp:740
CSeqGroup::count_and_union_sync
bm::id64_t count_and_union_sync(const bm::bvector<> &bv)
Definition xsample07a.cpp:725
CSeqGroup::CSeqGroup
CSeqGroup(bm::id64_t lead_id=~0ull)
Definition xsample07a.cpp:689
CSeqGroup::add_member_sync
void add_member_sync(bm::id64_t id, const bm::bvector<> &bv_kmer)
Definition xsample07a.cpp:711
CSeqGroup::get_members
const bm::bvector & get_members() const
Definition xsample07a.cpp:736
CSeqGroup::set_lead
void set_lead(bm::id64_t lead_id)
set id for the group representative
Definition xsample07a.cpp:696
CSeqGroup::get_rep
const bm::bvector & get_rep() const
Definition xsample07a.cpp:734
CSeqGroup::add_member
void add_member(bm::id64_t id)
add a member to the group
Definition xsample07a.cpp:705
CSeqGroup::CSeqClusters
friend class CSeqClusters
Definition xsample07a.cpp:742
CSeqGroup::get_kmer_union
bm::bvector & get_kmer_union()
Definition xsample07a.cpp:739
CSeqGroup::get_members
bm::bvector & get_members()
Definition xsample07a.cpp:737
CSeqGroup::merge_member_sync
void merge_member_sync(bm::bvector<> &bv_seq, bm::bvector<> &bv_kmer)
Definition xsample07a.cpp:718
CSeqGroup::is_assigned
bool is_assigned()
check is cluster is non-empty
Definition xsample07a.cpp:702
CSeqGroup::get_lead
bm::id64_t get_lead() const
Get lead id.
Definition xsample07a.cpp:699
CSeqGroup::clear_member
void clear_member(bm::id64_t id)
Definition xsample07a.cpp:731
CSequenceColl
Collection of sequences and k-mer fingerprint vectors.
Definition xsample07a.cpp:131
CSequenceColl::add_sequence
void add_sequence(const string &acc, vector_char_type *seq_ptr)
Definition xsample07a.cpp:140
CSequenceColl::CSequenceColl
CSequenceColl(const CSequenceColl &)=delete
CSequenceColl::deserialize_k_mers
void deserialize_k_mers(bvector_ptr_vector_type &k_mers_vect, const bm::bvector<> &bv_req, bm::bvector<>::size_type bv_req_count) const
Deserialize group of k-mer fingerprint vectors.
Definition xsample07a.cpp:208
CSequenceColl::total_seq_size
size_t total_seq_size() const
Definition xsample07a.cpp:171
CSequenceColl::get_sequence
const vector_char_type & get_sequence(size_t i) const
Definition xsample07a.cpp:167
CSequenceColl::buffer_type
std::vector< unsigned char > buffer_type
Definition xsample07a.cpp:133
CSequenceColl::get_buf
const unsigned char * get_buf(size_t i) const
Get k-mer BLOB pointer.
Definition xsample07a.cpp:186
CSequenceColl::get_buf_size
size_t get_buf_size(size_t i) const
Get k-mer vector BLOB size.
Definition xsample07a.cpp:183
CSequenceColl::size
size_t size() const
Definition xsample07a.cpp:163
CSequenceColl::buf_size
size_t buf_size() const
Definition xsample07a.cpp:180
CSequenceColl::seq_size
size_t seq_size(size_t i) const
Definition xsample07a.cpp:169
CSequenceColl::CSequenceColl
CSequenceColl()
Definition xsample07a.cpp:136
CSequenceColl::get_acc
const string & get_acc(size_t i) const
Definition xsample07a.cpp:166
CSequenceColl::sync_buffers_size
void sync_buffers_size()
Definition xsample07a.cpp:158
CSequenceColl::set_buffer
void set_buffer(size_t i, const buffer_type &buf)
Definition xsample07a.cpp:146
bm::aggregator
Algorithms for fast aggregation of a group of bit-vectors.
Definition bmaggregator.h:122
bm::bvector::enumerator
Constant iterator designed to enumerate "ON" bits.
Definition bm.h:603
bm::bvector::enumerator::go_to
bool go_to(size_type pos) BMNOEXCEPT
go to a specific position in the bit-vector (or next)
Definition bm.h:8209
bm::bvector::iterator_base::valid
bool valid() const BMNOEXCEPT
Checks if iterator is still valid.
Definition bm.h:283
bm::bvector
Bitvector Bit-vector container with runtime compression of bits.
Definition bm.h:115
bm::bvector::opt_compress
@ opt_compress
compress blocks when possible (GAP/prefix sum)
Definition bm.h:137
bm::bvector::opt_none
@ opt_none
no optimization
Definition bm.h:134
bm::bvector::test
bool test(size_type n) const BMNOEXCEPT
returns true if bit n is set and false is bit n is 0.
Definition bm.h:1502
bm::bvector::merge
void merge(bm::bvector< Alloc > &bvect)
Merge/move content from another vector.
Definition bm.h:5816
bm::bvector::allocator_pool_type
allocator_type::allocator_pool_type allocator_pool_type
Definition bm.h:118
bm::bvector::size
size_type size() const BMNOEXCEPT
Returns bvector's capacity (number of bits it can store).
Definition bm.h:1300
bm::bvector::count
size_type count() const BMNOEXCEPT
population count (count of ON bits)
Definition bm.h:2401
bm::bvector::bit_and
bm::bvector< Alloc > & bit_and(const bm::bvector< Alloc > &bv1, const bm::bvector< Alloc > &bv2, typename bm::bvector< Alloc >::optmode opt_mode=opt_none)
3-operand AND : this := bv1 AND bv2
Definition bm.h:6118
bm::bvector::set
bvector< Alloc > & set(size_type n, bool val=true)
Sets bit n if val is true, clears bit n if val is false.
Definition bm.h:4188
bm::bvector::set_allocator_pool
void set_allocator_pool(allocator_pool_type *pool_ptr) BMNOEXCEPT
Set allocator pool for local (non-th readed) memory cyclic(lots of alloc-free ops) opertations.
Definition bm.h:1040
bm::bvector::optimize
void optimize(bm::word_t *temp_block=0, optmode opt_mode=opt_compress, statistics *stat=0)
Optimize memory bitvector's memory allocation.
Definition bm.h:3635
bm::bvector::any
bool any() const BMNOEXCEPT
Returns true if any bits in this bitset are set, and otherwise returns false.
Definition bm.h:2451
bm::bvector::size_type
bvector_size_type size_type
Definition bm.h:121
bm::bvector::first
enumerator first() const
Returns enumerator pointing on the first non-zero bit.
Definition bm.h:1871
bm::bvector::swap
void swap(bvector< Alloc > &bvect) BMNOEXCEPT
Exchanges content of bv and this bvector.
Definition bm.h:3966
bm::bvector::count_range
size_type count_range(size_type left, size_type right, const rs_index_type &rs_idx) const BMNOEXCEPT
Returns count of 1 bits in the given range [left..right] Uses rank-select index to accelerate the sea...
Definition bm.h:3516
bm::bvector::set_range
bvector< Alloc > & set_range(size_type left, size_type right, bool value=true)
Sets all bits in the specified closed interval [left,right] Interval must be inside the bvector's siz...
Definition bm.h:2368
bm::bvector::get_first
size_type get_first() const BMNOEXCEPT
find first 1 bit in vector. Function may return 0 and this requires an extra check if bit 0 is actual...
Definition bm.h:1600
bm::bvector::clear
void clear(const size_type *ids, size_type ids_size, bm::sort_order so=bm::BM_UNKNOWN)
clear list of bits in this bitset
Definition bm.h:4149
bm::bvector::allocator_type
Alloc allocator_type
Definition bm.h:117
bm::chrono_taker
Utility class to collect performance measurements and statistics.
Definition bmtimer.h:41
bm::chrono_taker::duration_map_type
std::map< std::string, statistics > duration_map_type
test name to duration map
Definition bmtimer.h:66
bm::chrono_taker::ct_time
@ ct_time
Definition bmtimer.h:59
bm::chrono_taker::print_duration_map
static void print_duration_map(TOut &tout, const duration_map_type &dmap, format fmt=ct_time)
Definition bmtimer.h:150
bm::operation_deserializer
Deserializer, performs logical operations between bit-vector and serialized bit-vector.
Definition bmserial.h:927
bm::operation_deserializer::deserialize
size_type deserialize(bvector_type &bv, const unsigned char *buf, set_operation op, bool exit_on_one=false)
Deserialize bvector using buffer as set operation argument.
Definition bmserial.h:6581
bm::random_subset
Definition bmrandom.h:57
bm::random_subset::sample
void sample(BV &bv_out, const BV &bv_in, size_type sample_count)
Get random subset of input vector.
Definition bmrandom.h:140
bm::serializer
Bit-vector serialization class.
Definition bmserial.h:76
bm::serializer::set_bookmarks
void set_bookmarks(bool enable, unsigned bm_interval=256) BMNOEXCEPT
Add skip-markers to serialization BLOB for faster range decode at the expense of some BLOB size incre...
Definition bmserial.h:1284
bm::serializer::serialize
size_type serialize(const BV &bv, unsigned char *buf, size_t buf_size)
Bitvector serialization into memory block.
Definition bmserial.h:2703
bm::bvector::mem_pool_guard
bm::alloc_pool_guard< allocator_pool_type, bvector< Alloc > > mem_pool_guard
Definition bm.h:790
bm::BM_SORTED
@ BM_SORTED
input set is sorted (ascending order)
Definition bmconst.h:206
bm::set_COUNT_AND
@ set_COUNT_AND
Definition bmconst.h:175
bm::BM_GAP
@ BM_GAP
GAP compression is ON.
Definition bmconst.h:148
bm::deserialize
size_t deserialize(BV &bv, const unsigned char *buf, bm::word_t *temp_block=0, const bm::bv_ref_vector< BV > *ref_vect=0)
Bitvector deserialization from a memory BLOB.
Definition bmserial.h:3137
bm::rank_range_split
void rank_range_split(const BV &bv, typename BV::size_type rank, PairVect &target_v)
Algorithm to identify bit-vector ranges (splits) for the rank.
Definition bmalgo.h:394
bm::count_and
BV::size_type count_and(const BV &bv1, const BV &bv2) BMNOEXCEPT
Computes bitcount of AND operation of two bitsets.
Definition bmalgo.h:49
bm
Definition bm.h:78
bm::id64_t
unsigned long long int id64_t
Definition bmconst.h:35
timing_map
bm::chrono_taker ::duration_map_type timing_map
Definition sample11.cpp:46
main
int main(void)
Definition sample1.cpp:43
bv_size_type
bm::bvector ::size_type bv_size_type
Definition sample24.cpp:54
bv_ranges_vector
std::vector< std::pair< bv_size_type, bv_size_type > > bv_ranges_vector
Definition sample24.cpp:55
CKMerAcc
Utility class to accumulate cahnges to cluster before commiting it (mutex syncronous operation).
Definition xsample07a.cpp:1215
CKMerAcc::add
void add(size_t cluster_id, bm::bvector<>::size_type m_id, bm::bvector<> &bv_kmer)
Definition xsample07a.cpp:1220
CKMerAcc::bv_kmers
bvector_ptr_vector_type bv_kmers
Definition xsample07a.cpp:1238
CKMerAcc::bv_members
bvector_ptr_vector_type bv_members
Definition xsample07a.cpp:1237
CKMerAcc::CKMerAcc
CKMerAcc(size_t sz)
Definition xsample07a.cpp:1216
bm::bv_statistics::max_serialize_mem
size_t max_serialize_mem
estimated maximum memory for serialization
Definition bmfunc.h:61
bm::bvector::statistics
Statistical information about bitset's memory allocation details.
Definition bm.h:125
parse_args
static int parse_args(int argc, char *argv[])
Definition xsample03.cpp:110
is_bench
bool is_bench
Definition xsample03.cpp:107
is_diag
bool is_diag
Definition xsample03.cpp:105
is_timing
bool is_timing
Definition xsample03.cpp:106
ifa_name
std::string ifa_name
Definition xsample04.cpp:72
parallel_jobs
unsigned parallel_jobs
Definition xsample04a.cpp:68
vector_char_type
std::vector< char > vector_char_type
Definition xsample06.cpp:133
ikd_rep_name
std::string ikd_rep_name
Definition xsample07.cpp:84
ik_size
unsigned ik_size
Definition xsample07.cpp:89
kh_name
std::string kh_name
Definition xsample07.cpp:83
ikd_freq_name
std::string ikd_freq_name
Definition xsample07.cpp:85
ikd_counts_name
std::string ikd_counts_name
Definition xsample07.cpp:82
ikd_name
std::string ikd_name
Definition xsample07.cpp:81
f_percent
unsigned f_percent
Definition xsample07.cpp:91
k_mer_progress_count
std::atomic_ullong k_mer_progress_count(0)
compute_jaccard_clusters
static void compute_jaccard_clusters(CSeqClusters &seq_clusters, const CSequenceColl &seq_coll, unsigned num_clust, float similarity_cut_off, unsigned concurrency)
Definition xsample07a.cpp:1420
assign_to_best_cluster
static void assign_to_best_cluster(CSeqClusters &cluster_groups, const CSequenceColl &seq_coll, const bm::bvector<> &bv_seq_ids, bm::bvector<>::size_type seq_id_from, bm::bvector<>::size_type seq_id_to)
Compute AND similarity to all available clusters assign to the most similar using cluster representat...
Definition xsample07a.cpp:1245
generate_k_mers
static void generate_k_mers(CSequenceColl &seq_coll, unsigned k_size, size_t from, size_t to)
Definition xsample07a.cpp:561
compute_random_clusters
static void compute_random_clusters(CSeqClusters &cluster_groups, const CSequenceColl &seq_coll, const bm::bvector<> &bv_total, bm::random_subset< bvector_type > &rsub, unsigned num_clust, float similarity_cut_off, unsigned concurrency)
Pick random sequences as cluster seed elements, try attach initial sequences based on weighted simila...
Definition xsample07a.cpp:1374
wait_for_slot
void wait_for_slot(FV &futures, unsigned *parallel_cnt, unsigned concurrency)
wait for any opening in a list of futures used to schedule parallel tasks with CPU overbooking contro...
Definition xsample07a.cpp:111
translate_kmer
void translate_kmer(std::string &dna, bm::id64_t kmer_code, unsigned k_size)
Translate k-mer code into ATGC DNA string.
Definition xsample07a.cpp:444
save_kmer_buffers
static void save_kmer_buffers(const std::string &fname, const CSequenceColl &seq_coll)
save k-mer vectors to a file
Definition xsample07a.cpp:294
bvector_type
bm::bvector bvector_type
Definition xsample07a.cpp:94
load_FASTA
static int load_FASTA(const std::string &fname, CSequenceColl &seq_coll)
Load multi-sequence FASTA.
Definition xsample07a.cpp:244
compute_and_sim_row
static void compute_and_sim_row(unsigned *row, const bm::bvector<> *bv_i, size_t i, const std::vector< std::unique_ptr< bvector_type > > &k_mers_vect)
Compute similarity distances for one row/vector (1:N) of distance matrix.
Definition xsample07a.cpp:890
load_kmer_buffers
static void load_kmer_buffers(const std::string &fname, CSequenceColl &seq_coll)
Load k-mer vectors.
Definition xsample07a.cpp:332
distance_matrix_type
bm::dynamic_heap_matrix< unsigned, bm::bvector<>::allocator_type > distance_matrix_type
Definition xsample07a.cpp:96
compute_group
static void compute_group(CSeqGroup &seq_group, const CSequenceColl &seq_coll, const bm::bvector<> &bv_exceptions, float similarity_cut_off)
Definition xsample07a.cpp:1102
get_kmer_code
bool get_kmer_code(const char *dna, size_t pos, unsigned k_size, bm::id64_t &k_mer)
Calculate k-mer as an unsigned long integer.
Definition xsample07a.cpp:402
vector_char_type
std::vector< char > vector_char_type
Definition xsample07a.cpp:95
int2DNA
char int2DNA(unsigned code)
Translate integer code to DNA letter.
Definition xsample07a.cpp:429
generate_k_mer_bvector
void generate_k_mer_bvector(BV &bv, const vector_char_type &seq_vect, unsigned k_size, std::vector< bm::id64_t > &k_buf, const bm::id64_t chunk_size=400000000)
This function turns each k-mer into an integer number and encodes it in a bit-vector (presense vector...
Definition xsample07a.cpp:474
assign_to_best_cluster_union
static void assign_to_best_cluster_union(CSeqClusters &cluster_groups, const CSequenceColl &seq_coll, const bm::bvector<> &bv_seq_ids, bm::bvector<>::size_type seq_id_from, bm::bvector<>::size_type seq_id_to)
Compute AND similarity to all available clusters assign to the most similar using UNION of k-mers in ...
Definition xsample07a.cpp:1315
generate_k_mers_parallel
static void generate_k_mers_parallel(CSequenceColl &seq_coll, unsigned k_size, unsigned concurrency)
Definition xsample07a.cpp:616
get_DNA_code
bool get_DNA_code(char bp, bm::id64_t &dna_code)
Definition xsample07a.cpp:374
compute_seq_group_union
static void compute_seq_group_union(CSeqGroup &seq_group, const CSequenceColl &seq_coll)
Compute union (Universe) of all k-mers in the cluster group Implemented as a OR of all k-mer fingerpr...
Definition xsample07a.cpp:973
resolve_duplicates
void resolve_duplicates(CSeqGroup &seq_group1, CSeqGroup &seq_group2, const CSequenceColl &seq_coll)
Resolve duplicate members between two groups.
Definition xsample07a.cpp:1155
k_mer_progress_count
std::atomic_ullong k_mer_progress_count(0)
compute_and_sim
static void compute_and_sim(distance_matrix_type &dm, const CSequenceColl &seq_coll, const bm::bvector<> &bv_mem, bm::bvector<>::size_type bv_mem_cnt, unsigned concurrency)
Compute similarity distances matrix (COUNT(AND(a, b)).
Definition xsample07a.cpp:910
bvector_ptr_vector_type
std::vector< std::unique_ptr< bvector_type > > bvector_ptr_vector_type
Definition xsample07a.cpp:97