BitMagic-C++
xsample04.cpp
/*
Copyright(c) 2018 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 xsample04.cpp
*/
/*! \file xsample04.cpp
\brief Example: DNA substring search
*/
#include <iostream>
#include <sstream>
#include <chrono>
#include <regex>
#include <time.h>
#include <stdio.h>
#include <stdexcept>
#include <memory>
#include <vector>
#include <chrono>
#include <map>
#include <utility>
#include <algorithm>
#include <unordered_map>
#include "bm.h"
#include "bmalgo.h"
#include "bmserial.h"
#include "bmaggregator.h"
#include "bmdbg.h"
#include "bmtimer.h"
using namespace std;
static
void show_help()
{
std::cerr
<< "BitMagic DNA Search Sample (c) 2018" << std::endl
<< "-fa file-name -- input FASTA file" << std::endl
<< "-s hi|lo -- run substring search benchmark" << std::endl
<< "-diag -- run diagnostics" << std::endl
<< "-timing -- collect timings" << std::endl
;
}
// Arguments
//
std::string ifa_name;
bool is_diag = false;
bool is_timing = false;
bool is_bench = false;
bool is_search = false;
bool h_word_set = true;
static
int parse_args(int argc, char *argv[])
{
for (int i = 1; i < argc; ++i)
{
std::string arg = argv[i];
if ((arg == "-h") || (arg == "--help"))
{
return 0;
}
if (arg == "-fa" || arg == "--fa")
{
if (i + 1 < argc)
{
ifa_name = argv[++i];
}
else
{
std::cerr << "Error: -fa requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-diag" || arg == "--diag" || arg == "-d" || arg == "--d")
is_diag = true;
if (arg == "-timing" || arg == "--timing" || arg == "-t" || arg == "--t")
is_timing = true;
if (arg == "-bench" || arg == "--bench" || arg == "-b" || arg == "--b")
is_bench = true;
if (arg == "-search" || arg == "--search" || arg == "-s" || arg == "--s")
{
is_search = true;
if (i + 1 < argc)
{
std::string a = argv[i+1];
if (a != "-")
{
if (a == "l" || a == "lo")
{
h_word_set = false;
++i;
}
else
if (a == "h" || a == "hi")
{
h_word_set = true;
++i;
}
}
}
}
} // for i
return 0;
}
// Global types
//
typedef std::map<std::string, unsigned> freq_map;
typedef std::vector<std::pair<unsigned, std::string> > dict_vect;
// ----------------------------------------------------------------------------
// FASTA format parser
static
int load_FASTA(const std::string& fname, std::vector<char>& seq_vect)
{
bm::chrono_taker tt1("1. Parse FASTA", 1, &timing_map);
seq_vect.resize(0);
std::ifstream fin(fname.c_str(), std::ios::in);
if (!fin.good())
return -1;
std::string line;
for (unsigned i = 0; std::getline(fin, line); ++i)
{
if (line.empty() ||
line.front() == '>')
continue;
for (std::string::iterator it = line.begin(); it != line.end(); ++it)
seq_vect.push_back(*it);
} // for
return 0;
}
/**
Utility for keeping all DNA finger print vectors and search
using various techniques
*/
{
public:
enum { eA = 0, eC, eG, eT, eN, eEnd };
/// Build fingerprint bit-vectors from the original sequence
///
void Build(const vector<char>& sequence)
{
// use bulk insert iterator (faster way to construct a bit-index)
//
for (size_t i = 0; i < sequence.size(); ++i)
{
unsigned pos = unsigned(i);
switch (sequence[i])
{
case 'A':
iA = pos;
break;
case 'C':
iC = pos;
break;
case 'G':
iG = pos;
break;
case 'T':
iT = pos;
break;
case 'N':
iN = pos;
break;
default:
break;
}
}
}
/// Return fingerprint bit-vector
const bm::bvector<>& GetVector(char letter) const
{
switch (letter)
{
case 'A':
return m_FPrintBV[eA];
case 'C':
return m_FPrintBV[eC];
case 'G':
return m_FPrintBV[eG];
case 'T':
return m_FPrintBV[eT];
case 'N':
return m_FPrintBV[eN];
default:
break;
}
throw runtime_error("Error. Invalid letter!");
}
/// Find word strings
/// using shift + and on fingerprint vectors
/// (horizontal, non-fused basic method)
///
void Find(const string& word, vector<unsigned>& res)
{
if (word.empty())
return;
bm::bvector<> bv(GetVector(word[0])); // step 1: copy first vector
// run series of shifts + logical ANDs
for (size_t i = 1; i < word.size(); ++i)
{
bv.shift_right(); // SHIFT the accumulator bit-vector
// get and AND the next fingerprint
const bm::bvector<>& bv_mask = GetVector(word[i]);
bv &= bv_mask;
auto any = bv.any();
if (!any)
break;
}
// translate results from bvector of word ends to result
unsigned ws = unsigned(word.size()) - 1;
TranslateResults(bv, ws, res);
};
/// This method uses cache blocked aggregator with fused SHIFT+AND kernel
///
void FindAggFused(const string& word, vector<unsigned>& res)
{
if (word.empty())
return;
// first we setup aggregator, add a group of vectors to be processed
m_Agg.reset();
for (size_t i = 0; i < word.size(); ++i)
{
const bm::bvector<>& bv_mask = GetVector(word[i]);
m_Agg.add(&bv_mask);
}
// now run the whole algorithm to get benefits of cache blocking
//
m_Agg.combine_shift_right_and(bv);
// translate results from bvector of word ends to result
unsigned ws = unsigned(word.size()) - 1;
TranslateResults(bv, ws, res);
};
/// Find a set of words in one pass using pipeline
/// of aggregators (this is very experimental)
///
void FindCollection(const vector<tuple<string,int> >& words,
vector<vector<unsigned>>& hits)
{
vector<unique_ptr<aggregator_type> > agg_pipeline;
unsigned ws = 0;
for (const auto& w : words)
{
unique_ptr<aggregator_type> agg_ptr(new aggregator_type());
agg_ptr->set_operation(aggregator_type::BM_SHIFT_R_AND);
const string& word = get<0>(w);
for (size_t i = 0; i < word.size(); ++i)
{
const bm::bvector<>& bv_mask = GetVector(word[i]);
agg_ptr->add(&bv_mask);
}
agg_pipeline.emplace_back(agg_ptr.release());
ws = unsigned(word.size()) - 1;
}
// run the pipeline
vector<unique_ptr<aggregator_type> >::iterator>(agg_pipeline.begin(), agg_pipeline.end());
// convert the results
for (size_t i = 0; i < agg_pipeline.size(); ++i)
{
const aggregator_type* agg_ptr = agg_pipeline[i].get();
auto bv = agg_ptr->get_target();
vector<unsigned> res;
res.reserve(12000);
TranslateResults(*bv, ws, res);
hits.emplace_back(res);
}
}
protected:
/// Translate search results vector using (word size) left shift
///
void TranslateResults(const bm::bvector<>& bv,
unsigned left_shift,
vector<unsigned>& res)
{
for (; en.valid(); ++en)
{
auto pos = *en;
res.push_back(pos - left_shift);
}
}
private:
bm::bvector<> m_FPrintBV[eEnd];
};
static const size_t WORD_SIZE = 28;
using THitList = vector<unsigned>;
/// generate the most frequent words of specified length from the input sequence
///
static
void generate_kmers(vector<tuple<string,int>>& top_words,
vector<tuple<string,int>>& lo_words,
const vector<char>& data,
size_t N,
unsigned word_size)
{
cout << "k-mer generation... " << endl;
top_words.clear();
lo_words.clear();
if (data.size() < word_size)
return;
size_t end_pos = data.size() - word_size;
size_t i = 0;
map<string, int> words;
while (i < end_pos)
{
string s(&data[i], word_size);
if (s.find('N') == string::npos)
words[s] += 1;
i += word_size;
if (i % 10000 == 0)
{
cout << "\r" << i << "/" << end_pos << flush;
}
}
cout << endl << "Picking k-mer samples..." << flush;
multimap<int,string, greater<int>> dst;
for_each(words.begin(), words.end(), [&](const std::pair<string,int>& p)
{
dst.emplace(p.second, p.first);
});
{
auto it = dst.begin();
for(size_t count = 0; count < N && it !=dst.end(); ++it,++count)
top_words.emplace_back(it->second, it->first);
}
{
auto it = dst.rbegin();
for(size_t count = 0; count < N && it !=dst.rend(); ++it, ++count)
lo_words.emplace_back(it->second, it->first);
}
cout << "OK" << endl;
}
/// 2-way string matching
///
static
void find_word_2way(vector<char>& data,
const char* word, unsigned word_size,
{
if (data.size() < word_size)
return;
size_t i = 0;
size_t end_pos = data.size() - word_size;
while (i < end_pos)
{
bool found = true;
for (size_t j = i, k = 0, l = word_size - 1; l > k; ++j, ++k, --l)
{
if (data[j] != word[k] || data[i + l] != word[l])
{
found = false;
break;
}
}
if (found)
r.push_back(unsigned(i));
++i;
}
}
/// Find all words in one pass (cache coherent algorithm)
/// (variation of 2-way string matching for collection search)
///
static
void find_words(const vector<char>& data,
vector<const char*> words,
unsigned word_size,
vector<vector<unsigned>>& hits)
{
if (data.size() < word_size)
return;
size_t i = 0;
size_t end_pos = data.size() - word_size;
size_t words_size = words.size();
while (i < end_pos)
{
for (size_t idx = 0; idx < words_size; ++idx)
{
auto& word = words[idx];
bool found = true;
for (size_t j = i, k = 0, l = word_size - 1; l > k; ++j, ++k, --l)
{
if (data[j] != word[k] || data[i + l] != word[l])
{
found = false;
break;
}
} // for
if (found)
{
hits[idx].push_back(unsigned(i));
break;
}
} // for
++i;
} // while
}
/// Check search result match
///
static
bool hitlist_compare(const THitList& h1, const THitList& h2)
{
if (h1.size() != h2.size())
{
cerr << "size1 = " << h1.size() << " size2 = " << h2.size() << endl;
return false;
}
for (size_t i = 0; i < h1.size(); ++i)
{
if (h1[i] != h2[i])
return false;
}
return true;
}
int main(int argc, char *argv[])
{
if (argc < 3)
{
return 1;
}
std::vector<char> seq_vect;
try
{
auto ret = parse_args(argc, argv);
if (ret != 0)
return ret;
if (!ifa_name.empty())
{
auto res = load_FASTA(ifa_name, seq_vect);
if (res != 0)
return res;
std::cout << "FASTA sequence size=" << seq_vect.size() << std::endl;
{
bm::chrono_taker tt1("2. Build DNA index", 1, &timing_map);
idx.Build(seq_vect);
}
}
if (is_search)
{
vector<tuple<string,int> > h_words;
vector<tuple<string,int> > l_words;
vector<tuple<string,int>>& words = h_word_set ? h_words : l_words;
// generate search sets for benchmarking
//
generate_kmers(h_words, l_words, seq_vect, 25, WORD_SIZE);
vector<THitList> word_hits;
vector<THitList> word_hits_agg;
// search all words in one pass and
// store results in list of hits according to the order of words
// (this method uses memory proximity
// of searched words to maximize CPU cache hit rate)
{
vector<const char*> word_list;
for (const auto& w : words)
{
word_list.push_back(get<0>(w).c_str());
}
word_hits.resize(words.size());
for_each(word_hits.begin(), word_hits.end(), [](THitList& ht) {
ht.reserve(12000);
});
bm::chrono_taker tt1("6. String search 2-way single pass",
unsigned(words.size()), &timing_map);
find_words(seq_vect, word_list, unsigned(WORD_SIZE), word_hits);
}
// collection search, runs all hits at once
//
{
bm::chrono_taker tt1("7. Aggregated search single pass",
unsigned(words.size()), &timing_map);
idx.FindCollection(words, word_hits_agg);
}
// a few variants of word-by-word searches
//
for (size_t word_idx = 0; word_idx < words.size(); ++ word_idx)
{
auto& word = get<0>(words[word_idx]);
THitList hits1;
{
bm::chrono_taker tt1("3. String search 2-way", 1, &timing_map);
find_word_2way(seq_vect,
word.c_str(), unsigned(word.size()),
hits1);
}
THitList hits2;
{
bm::chrono_taker tt1("4. Search with bvector SHIFT+AND", 1, &timing_map);
idx.Find(word, hits2);
}
THitList hits4;
{
bm::chrono_taker tt1("5. Search with aggregator fused SHIFT+AND", 1, &timing_map);
idx.FindAggFused(word, hits4);
}
// check correctness
if (!hitlist_compare(hits1, hits2)
|| !hitlist_compare(hits2, hits4))
{
cout << "Mismatch ERROR for: " << word << endl;
}
else
if (!hitlist_compare(word_hits[word_idx], hits1)
|| !hitlist_compare(word_hits_agg[word_idx], hits1))
{
cout << "Sigle pass mismatch ERROR for: " << word << endl;
}
else
{
cout << word_idx << ": " << word << ": " << hits1.size() << " hits " << endl;
}
}
}
if (is_timing) // print all collected timings
{
std::cout << std::endl << "Performance:" << std::endl;
}
}
catch (std::exception& ex)
{
std::cerr << "Error:" << ex.what() << std::endl;
return 1;
}
return 0;
}
bmaggregator.h
Algorithms for fast aggregation of N bvectors.
bm::chrono_taker
Utility class to collect performance measurements and statistics.
Definition: bmtimer.h:39
generate_kmers
static void generate_kmers(vector< tuple< string, int >> &top_words, vector< tuple< string, int >> &lo_words, const vector< char > &data, size_t N, unsigned word_size)
generate the most frequent words of specified length from the input sequence
Definition: xsample04.cpp:367
DNA_FingerprintScanner::Find
void Find(const string &word, vector< unsigned > &res)
Find word strings using shift + and on fingerprint vectors (horizontal, non-fused basic method)
Definition: xsample04.cpp:250
bm::bvector::any
bool any() const BMNOEXCEPT
Returns true if any bits in this bitset are set, and otherwise returns false.
Definition: bm.h:2244
bm::BM_SORTED
input set is sorted (ascending order)
Definition: bmconst.h:191
DNA_FingerprintScanner::Build
void Build(const vector< char > &sequence)
Build fingerprint bit-vectors from the original sequence.
Definition: xsample04.cpp:189
bm::aggregator
Algorithms for fast aggregation of a group of bit-vectors.
Definition: bmaggregator.h:58
bm::bvector::enumerator
Constant iterator designed to enumerate "ON" bits.
Definition: bm.h:599
DNA_FingerprintScanner
Utility for keeping all DNA finger print vectors and search using various techniques.
Definition: xsample04.cpp:180
bm::aggregator::get_target
const bvector_type * get_target() const
Definition: bmaggregator.h:326
DNA_FingerprintScanner::FindCollection
void FindCollection(const vector< tuple< string, int > > &words, vector< vector< unsigned >> &hits)
Find a set of words in one pass using pipeline of aggregators (this is very experimental)
Definition: xsample04.cpp:302
is_timing
bool is_timing
Definition: xsample04.cpp:74
bmalgo.h
Algorithms for bvector<> (main include)
show_help
static void show_help()
Definition: xsample04.cpp:56
aggregator_type
bm::aggregator< bm::bvector<> > aggregator_type
Definition: xsample04.cpp:144
is_diag
bool is_diag
Definition: xsample04.cpp:73
find_words
static void find_words(const vector< char > &data, vector< const char * > words, unsigned word_size, vector< vector< unsigned >> &hits)
Find all words in one pass (cache coherent algorithm) (variation of 2-way string matching for collect...
Definition: xsample04.cpp:451
bm::bvector<>
freq_map
std::map< std::string, unsigned > freq_map
Definition: xsample04.cpp:141
bm::bvector::iterator_base::valid
bool valid() const BMNOEXCEPT
Checks if iterator is still valid. Analog of != 0 comparison for pointers.
Definition: bm.h:280
timing_map
bm::chrono_taker::duration_map_type timing_map
Definition: xsample04.cpp:148
bmtimer.h
Timing utilities for benchmarking (internal)
is_search
bool is_search
Definition: xsample04.cpp:76
bm::chrono_taker::ct_all
Definition: bmtimer.h:60
ifa_name
std::string ifa_name
Definition: xsample04.cpp:72
is_bench
bool is_bench
Definition: xsample04.cpp:75
dict_vect
std::vector< std::pair< unsigned, std::string > > dict_vect
Definition: xsample04.cpp:142
parse_args
static int parse_args(int argc, char *argv[])
Definition: xsample04.cpp:80
h_word_set
bool h_word_set
Definition: xsample04.cpp:77
bm::aggregator_pipeline_execute
void aggregator_pipeline_execute(It first, It last)
Experimental method ro run multiple aggregators in sync.
Definition: bmaggregator.h:481
DNA_FingerprintScanner::FindAggFused
void FindAggFused(const string &word, vector< unsigned > &res)
This method uses cache blocked aggregator with fused SHIFT+AND kernel.
Definition: xsample04.cpp:277
bm::chrono_taker::duration_map_type
std::map< std::string, statistics > duration_map_type
test name to duration map
Definition: bmtimer.h:65
bmserial.h
Serialization / compression of bvector<>. Set theoretical operations on compressed BLOBs.
bm::chrono_taker::print_duration_map
static void print_duration_map(const duration_map_type &dmap, format fmt=ct_time)
Definition: bmtimer.h:127
hitlist_compare
static bool hitlist_compare(const THitList &h1, const THitList &h2)
Check search result match.
Definition: xsample04.cpp:490
load_FASTA
static int load_FASTA(const std::string &fname, std::vector< char > &seq_vect)
Definition: xsample04.cpp:152
THitList
vector< unsigned > THitList
Definition: xsample04.cpp:362
main
int main(int argc, char *argv[])
Definition: xsample04.cpp:508
bm::bvector::first
enumerator first() const
Returns enumerator pointing on the first non-zero bit.
Definition: bm.h:1770
bm::bvector::bulk_insert_iterator
Output iterator iterator designed to set "ON" bits based on input sequence of integers.
Definition: bm.h:461
find_word_2way
static void find_word_2way(vector< char > &data, const char *word, unsigned word_size, THitList &r)
2-way string matching
Definition: xsample04.cpp:421
bm::aggregator< bm::bvector<> >::BM_SHIFT_R_AND
Definition: bmaggregator.h:79
bm.h
Compressed bit-vector bvector<> container, set algebraic methods, traversal iterators.
WORD_SIZE
static const size_t WORD_SIZE
Definition: xsample04.cpp:361