BitMagic-C++
xsample03.cpp
Go to the documentation of this file.
1/*
2Copyright(c) 2018 Anatoliy Kuznetsov(anatoliy_kuznetsov at yahoo.com)
3
4Licensed under the Apache License, Version 2.0 (the "License");
5you may not use this file except in compliance with the License.
6You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10Unless required by applicable law or agreed to in writing, software
11distributed under the License is distributed on an "AS IS" BASIS,
12WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13See the License for the specific language governing permissions and
14limitations under the License.
15
16For more information please visit: http://bitmagic.io
17*/
18
19/** \example xsample03.cpp
20 Seach for human mutation (SNP) in within chr1.
21 Benchmark comaprison of different methods
22
23 \sa bm::sparse_vector
24 \sa bm::rsc_sparse_vector
25 \sa bm::sparse_vector_scanner
26*/
27
28/*! \file xsample03.cpp
29 \brief Example: SNP search in human genome
30
31 Brief description of used method:
32 1. Parse SNP chromosome report and extract information about SNP number and
33 location in the chromosome
34 2. Use this information to build bit-transposed sparse_vector<>
35 where vector position matches chromosome position and SNP ids (aka rsid)
36 is kept as a bit-transposed matrix
37 3. Build rank-select compressed sparse vector, dropping all NULL columns
38 (this data format is pretty sparse, since number of SNPs is significantly
39 less than number of chromosome bases (1:5 or less)
40 Use memory report to understand memory footprint for each form of storage
41 4. Run benchmarks searching for 500 randomly selected SNPs using
42 - bm::sparse_vector<>
43 - bm::rsc_sparse_vector<>
44 - std::vector<pair<unsigned, unsigned> >
45
46 This example should be useful for construction of compressed columnar
47 tables with parallel search capabilities.
48
49*/
50
51#include <iostream>
52#include <sstream>
53#include <chrono>
54#include <regex>
55#include <time.h>
56#include <stdio.h>
57#include <vector>
58#include <chrono>
59#include <map>
60#include <utility>
61
62using namespace std;
63
64#include "bm.h"
65#include "bmalgo.h"
66#include "bmserial.h"
67#include "bmrandom.h"
68#include "bmsparsevec.h"
69#include "bmsparsevec_compr.h"
70#include "bmsparsevec_algo.h"
71#include "bmsparsevec_serial.h"
72#include "bmalgo_similarity.h"
73#include "bmsparsevec_util.h"
74
75
76#include "bmdbg.h"
77#include "bmtimer.h"
78#include "bmundef.h" /* clear the pre-proc defines from BM */
79
80static
82{
83 std::cerr
84 << "BitMagic SNP Search Sample Utility (c) 2018" << std::endl
85 << "-isnp file-name -- input set file (SNP FASTA) to parse" << std::endl
86 << "-svout spase vector output -- sparse vector name to save" << std::endl
87 << "-rscout rs-compressed spase vector output -- name to save" << std::endl
88 << "-svin sparse vector input -- sparse vector file name to load " << std::endl
89 << "-rscin rs-compressed sparse vector input -- file name to load " << std::endl
90 << "-diag -- run diagnostics" << std::endl
91 << "-timing -- collect timings" << std::endl
92 ;
93}
94
95
96
97
98// Arguments
99//
100std::string sv_out_name;
101std::string rsc_out_name;
102std::string sv_in_name;
103std::string rsc_in_name;
104std::string isnp_name;
105bool is_diag = false;
106bool is_timing = false;
107bool is_bench = false;
108
109static
110int parse_args(int argc, char *argv[])
111{
112 for (int i = 1; i < argc; ++i)
113 {
114 std::string arg = argv[i];
115 if ((arg == "-h") || (arg == "--help"))
116 {
117 show_help();
118 return 0;
119 }
120
121 if (arg == "-svout" || arg == "--svout")
122 {
123 if (i + 1 < argc)
124 {
125 sv_out_name = argv[++i];
126 }
127 else
128 {
129 std::cerr << "Error: -svout requires file name" << std::endl;
130 return 1;
131 }
132 continue;
133 }
134 if (arg == "-rscout" || arg == "--rscout")
135 {
136 if (i + 1 < argc)
137 {
138 rsc_out_name = argv[++i];
139 }
140 else
141 {
142 std::cerr << "Error: -rscout requires file name" << std::endl;
143 return 1;
144 }
145 continue;
146 }
147
148 if (arg == "-svin" || arg == "--svin")
149 {
150 if (i + 1 < argc)
151 {
152 sv_in_name = argv[++i];
153 }
154 else
155 {
156 std::cerr << "Error: -svin requires file name" << std::endl;
157 return 1;
158 }
159 continue;
160 }
161
162 if (arg == "-rscin" || arg == "--rscin")
163 {
164 if (i + 1 < argc)
165 {
166 rsc_in_name = argv[++i];
167 }
168 else
169 {
170 std::cerr << "Error: -rscin requires file name" << std::endl;
171 return 1;
172 }
173 continue;
174 }
175
176 if (arg == "-isnp" || arg == "--isnp" || arg == "-snp" || arg == "--snp")
177 {
178 if (i + 1 < argc)
179 {
180 isnp_name = argv[++i];
181 }
182 else
183 {
184 std::cerr << "Error: -isnp requires file name" << std::endl;
185 return 1;
186 }
187 continue;
188 }
189
190 if (arg == "-diag" || arg == "--diag" || arg == "-d" || arg == "--d")
191 is_diag = true;
192 if (arg == "-timing" || arg == "--timing" || arg == "-t" || arg == "--t")
193 is_timing = true;
194 if (arg == "-bench" || arg == "--bench" || arg == "-b" || arg == "--b")
195 is_bench = true;
196
197 } // for i
198 return 0;
199}
200
201
202// Global types
203//
206typedef std::vector<std::pair<unsigned, unsigned> > vector_pairs;
207
208// ----------------------------------------------------------------------------
209
211
212// SNP report format parser (extraction and transformation)
213// Parser extracts SNP id (rsid) and positio on genome to build
214// sparse vector where index (position in vector) metches position on the
215// chromosome 1.
216//
217static
218int load_snp_report(const std::string& fname, sparse_vector_u32& sv)
219{
220 bm::chrono_taker tt1(cout, "1. parse input SNP chr report", 1, &timing_map);
221
222 std::ifstream fin(fname.c_str(), std::ios::in);
223 if (!fin.good())
224 return -1;
225
226 unsigned rs_id, rs_pos;
227 size_t idx;
228
229 std::string line;
230 std::string delim = " \t";
231
232 std::regex reg("\\s+");
233 std::sregex_token_iterator it_end;
234
235 bm::bvector<> bv_rs;
236 bv_rs.init();
237
238 unsigned rs_cnt = 0;
239 for (unsigned i = 0; std::getline(fin, line); ++i)
240 {
241 if (line.empty() ||
242 !isdigit(line.front())
243 )
244 continue;
245
246 // regex based tokenizer
247 std::sregex_token_iterator it(line.begin(), line.end(), reg, -1);
248 std::vector<std::string> line_vec(it, it_end);
249 if (line_vec.empty())
250 continue;
251
252 // parse columns of interest
253 try
254 {
255 rs_id = unsigned(std::stoul(line_vec.at(0), &idx));
256
257 if (bv_rs.test(rs_id))
258 {
259 continue;
260 }
261 rs_pos = unsigned(std::stoul(line_vec.at(11), &idx));
262
263 bv_rs.set_bit_no_check(rs_id);
264 sv.set(rs_pos, rs_id);
265
266 ++rs_cnt;
267 }
268 catch (std::exception& /*ex*/)
269 {
270 continue; // detailed disgnostics commented out
271 // error detected, because some columns are sometimes missing
272 // just ignore it
273 //
274 /*
275 std::cerr << ex.what() << "; ";
276 std::cerr << "rs=" << line_vec.at(0) << " pos=" << line_vec.at(11) << std::endl;
277 continue;
278 */
279 }
280 if (rs_cnt % (4 * 1024) == 0)
281 std::cout << "\r" << rs_cnt << " / " << i; // PROGRESS report
282 } // for i
283
284 std::cout << std::endl;
285 std::cout << "SNP count=" << rs_cnt << std::endl;
286 return 0;
287}
288
289// Generate random subset of random values from a sparse vector
290//
291static
292void generate_random_subset(const sparse_vector_u32& sv, std::vector<unsigned>& vect, unsigned count)
293{
295
296 bm::random_subset<bm::bvector<> > rand_sampler;
297 bm::bvector<> bv_sample;
298 rand_sampler.sample(bv_sample, *bv_null, count);
299
300 bm::bvector<>::enumerator en = bv_sample.first();
301 for (; en.valid(); ++en)
302 {
303 unsigned idx = *en;
304 unsigned v = sv[idx];
305 vect.push_back(v);
306 }
307}
308
309// build std::vector of pairs (position to rs)
310//
311static
313{
316
317 for (; it != it_end; ++it)
318 {
319 if (!it.is_null())
320 {
321 std::pair<unsigned, unsigned> pos2rs = std::make_pair(it.pos(), it.value());
322 vp.push_back(pos2rs);
323 }
324 }
325}
326
327// O(N) -- O(N/2) linear search in vector of pairs (position - rsid)
328//
329static
330bool search_vector_pairs(const vector_pairs& vp, unsigned rs_id, unsigned& pos)
331{
332 for (unsigned i = 0; i < vp.size(); ++i)
333 {
334 if (vp[i].second == rs_id)
335 {
336 pos = vp[i].first;
337 return true;
338 }
339 }
340 return false;
341}
342
343// SNP search benchmark
344// Search for SNPs using different data structures (Bitmagic and STL)
345//
346// An extra step is verification, to make sure all methods are consistent
347//
348static
350{
351 const unsigned rs_sample_count = 2000;
352
353 std::vector<unsigned> rs_vect;
354 generate_random_subset(sv, rs_vect, rs_sample_count);
355 if (rs_vect.empty())
356 {
357 std::cerr << "Benchmark subset empty!" << std::endl;
358 return;
359 }
360
361 // build traditional sparse vector
362 vector_pairs vp;
363 build_vector_pairs(sv, vp);
364
365 // search result bit-vectors
366 //
367 bm::bvector<> bv_found1;
368 bm::bvector<> bv_found2;
369 bm::bvector<> bv_found3;
370
371 bv_found1.init(); bv_found2.init(); bv_found3.init();// pre-initialize vectors
372
373 if (!sv.empty())
374 {
375 bm::chrono_taker tt1(cout, "09. rs search (sv)", unsigned(rs_vect.size()), &timing_map);
376
377 // scanner class
378 // it's important to keep scanner class outside the loop to avoid
379 // unnecessary re-allocs and construction costs
380 //
381
383
384 for (unsigned i = 0; i < rs_vect.size(); ++i)
385 {
386 unsigned rs_id = rs_vect[i];
387 unsigned rs_pos;
388 bool found = scanner.find_eq(sv, rs_id, rs_pos);
389
390 if (found)
391 {
392 bv_found1.set_bit_no_check(rs_pos);
393 }
394 else
395 {
396 std::cout << "Error: rs_id = " << rs_id << " not found!" << std::endl;
397 }
398 } // for
399 }
400
401 if (!csv.empty())
402 {
403 bm::chrono_taker tt1(cout, "10. rs search (rsc_sv)", unsigned(rs_vect.size()), &timing_map);
404
406
407 for (unsigned i = 0; i < rs_vect.size(); ++i)
408 {
409 unsigned rs_id = rs_vect[i];
410 unsigned rs_pos;
411 bool found = scanner.find_eq(csv, rs_id, rs_pos);
412
413 if (found)
414 {
415 bv_found2.set_bit_no_check(rs_pos);
416 }
417 else
418 {
419 std::cout << "rs_id = " << rs_id << " not found!" << std::endl;
420 }
421 } // for
422 }
423
424 if (vp.size())
425 {
426 bm::chrono_taker tt1(cout, "11. rs search (std::vector<>)", unsigned(rs_vect.size()), &timing_map);
427
428 for (unsigned i = 0; i < rs_vect.size(); ++i)
429 {
430 unsigned rs_id = rs_vect[i];
431 unsigned rs_pos;
432 bool found = search_vector_pairs(vp, rs_id, rs_pos);
433
434 if (found)
435 {
436 bv_found3.set_bit_no_check(rs_pos);
437 }
438 else
439 {
440 std::cout << "rs_id = " << rs_id << " not found!" << std::endl;
441 }
442 } // for
443 }
444
445 // compare results from various methods (check integrity)
446 int res = bv_found1.compare(bv_found2);
447 if (res != 0)
448 {
449 std::cerr << "Error: search discrepancy (sparse search) detected!" << std::endl;
450 }
451 res = bv_found1.compare(bv_found3);
452 if (res != 0)
453 {
454 std::cerr << "Error: search discrepancy (std::vector<>) detected!" << std::endl;
455 }
456
457}
458
459
460int main(int argc, char *argv[])
461{
462 if (argc < 3)
463 {
464 show_help();
465 return 1;
466 }
467
470
471 try
472 {
473 auto ret = parse_args(argc, argv);
474 if (ret != 0)
475 return ret;
476
477 if (!isnp_name.empty())
478 {
479 auto res = load_snp_report(isnp_name, sv);
480 if (res != 0)
481 {
482 return res;
483 }
484 }
485 if (!sv_in_name.empty())
486 {
487 bm::chrono_taker tt1(cout, "02. Load sparse vector", 1, &timing_map);
488 file_load_svector(sv, sv_in_name);
489 }
490
491 // load rs-compressed sparse vector
492 if (!rsc_in_name.empty())
493 {
494 bm::chrono_taker tt1(cout, "03. Load rsc sparse vector", 1, &timing_map);
495 file_load_svector(csv, rsc_in_name);
496 }
497
498 // if rs-compressed vector is not available - build it on the fly
499 if (csv.empty() && !sv.empty())
500 {
502 {
503 bm::chrono_taker tt1(cout, "04. rs compress sparse vector", 1, &timing_map);
504 csv.load_from(sv);
505 }
506 {
507 bm::chrono_taker tt1(cout, "05. rs de-compress sparse vector", 1, &timing_map);
508 csv.load_to(sv2);
509 }
510
511 if (!sv.equal(sv2)) // diagnostics check (just in case)
512 {
513 std::cerr << "Error: rs-compressed vector check failed!" << std::endl;
514 return 1;
515 }
516 }
517
518 // save SV vector
519 if (!sv_out_name.empty() && !sv.empty())
520 {
521 bm::chrono_taker tt1(cout, "07. Save sparse vector", 1, &timing_map);
522 sv.optimize();
523 file_save_svector(sv, sv_out_name);
524 }
525
526 // save RS sparse vector
527 if (!rsc_out_name.empty() && !csv.empty())
528 {
529 bm::chrono_taker tt1(cout, "08. Save RS sparse vector", 1, &timing_map);
530 csv.optimize();
531 file_save_svector(csv, rsc_out_name);
532 }
533
534 if (is_diag) // print memory diagnostics
535 {
536 if (!sv.empty())
537 {
538 std::cout << std::endl
539 << "sparse vector statistics:"
540 << std::endl;
541 bm::print_svector_stat(std::cout, sv, true);
542 }
543 if (!csv.empty())
544 {
545 std::cout << std::endl
546 << "RS compressed sparse vector statistics:"
547 << std::endl;
548 bm::print_svector_stat(std::cout, csv, true);
549 }
550 }
551
552 if (is_bench) // run set of benchmarks
553 {
554 run_benchmark(sv, csv);
555 }
556
557 if (is_timing) // print all collected timings
558 {
559 std::cout << std::endl << "Performance:" << std::endl;
561 }
562 }
563 catch (std::exception& ex)
564 {
565 std::cerr << "Error:" << ex.what() << std::endl;
566 return 1;
567 }
568
569 return 0;
570}
571
Compressed bit-vector bvector<> container, set algebraic methods, traversal iterators.
Algorithms for bvector<> (main include)
Generation of random subset.
Serialization / compression of bvector<>. Set theoretical operations on compressed BLOBs.
Sparse constainer sparse_vector<> for integer types using bit-transposition transform.
Algorithms for bm::sparse_vector.
Compressed sparse container rsc_sparse_vector<> for integer types.
Serialization for sparse_vector<>
Timing utilities for benchmarking (internal)
pre-processor un-defines to avoid global space pollution (internal)
const bvector_type * get_null_bvector() const BMNOEXCEPT
Get bit-vector of assigned values or NULL (if not constructed that way)
Definition: bmbmatrix.h:430
Constant iterator designed to enumerate "ON" bits.
Definition: bm.h:603
bool valid() const BMNOEXCEPT
Checks if iterator is still valid.
Definition: bm.h:283
Bitvector Bit-vector container with runtime compression of bits.
Definition: bm.h:115
bool test(size_type n) const BMNOEXCEPT
returns true if bit n is set and false is bit n is 0.
Definition: bm.h:1480
enumerator first() const
Returns enumerator pointing on the first non-zero bit.
Definition: bm.h:1849
void init()
Explicit post-construction initialization. Must be caled to make sure safe use of *_no_check() method...
Definition: bm.h:2258
int compare(const bvector< Alloc > &bvect) const BMNOEXCEPT
Lexicographical comparison with a bitvector.
Definition: bm.h:3709
void set_bit_no_check(size_type n)
Set bit without checking preconditions (size, etc)
Definition: bm.h:4405
Utility class to collect performance measurements and statistics.
Definition: bmtimer.h:41
std::map< std::string, statistics > duration_map_type
test name to duration map
Definition: bmtimer.h:66
static void print_duration_map(TOut &tout, const duration_map_type &dmap, format fmt=ct_time)
Definition: bmtimer.h:150
void sample(BV &bv_out, const BV &bv_in, size_type sample_count)
Get random subset of input vector.
Definition: bmrandom.h:140
void load_from(const sparse_vector_type &sv_src)
Load compressed vector from a sparse vector (with NULLs)
void load_to(sparse_vector_type &sv) const
Exort compressed vector to a sparse vector (with NULLs)
void optimize(bm::word_t *temp_block=0, typename bvector_type::optmode opt_mode=bvector_type::opt_compress, statistics *stat=0)
run memory optimization for all vector slices
bool empty() const BMNOEXCEPT
return true if vector is empty
algorithms for sparse_vector scan/search
void find_eq(const SV &sv, value_type value, bvector_type &bv_out)
find all sparse vector elements EQ to search value
succinct sparse vector with runtime compression using bit-slicing / transposition method
Definition: bmsparsevec.h:87
bool equal(const sparse_vector< Val, BV > &sv, bm::null_support null_able=bm::use_null) const BMNOEXCEPT
check if another sparse vector has the same content and size
Definition: bmsparsevec.h:2246
void push_back(value_type v)
push value back into vector
Definition: bmsparsevec.h:1838
void set(size_type idx, value_type v)
set specified element with bounds checking and automatic resize
Definition: bmsparsevec.h:1804
const_iterator end() const BMNOEXCEPT
Provide const iterator access to the end
Definition: bmsparsevec.h:559
bool empty() const BMNOEXCEPT
return true if vector is empty
Definition: bmsparsevec.h:716
const_iterator begin() const BMNOEXCEPT
Provide const iterator access to container content
Definition: bmsparsevec.h:2256
void optimize(bm::word_t *temp_block=0, typename bvector_type::optmode opt_mode=bvector_type::opt_compress, typename sparse_vector< Val, BV >::statistics *stat=0)
run memory optimization for all vector planes
Definition: bmsparsevec.h:2094
@ use_null
support "non-assigned" or "NULL" logic
Definition: bmconst.h:229
int main(int argc, char *argv[])
Definition: xsample03.cpp:460
static int parse_args(int argc, char *argv[])
Definition: xsample03.cpp:110
std::vector< std::pair< unsigned, unsigned > > vector_pairs
Definition: xsample03.cpp:206
bool is_bench
Definition: xsample03.cpp:107
std::string sv_in_name
Definition: xsample03.cpp:102
std::string rsc_in_name
Definition: xsample03.cpp:103
static void build_vector_pairs(const sparse_vector_u32 &sv, vector_pairs &vp)
Definition: xsample03.cpp:312
bm::sparse_vector< unsigned, bm::bvector<> > sparse_vector_u32
Definition: xsample03.cpp:204
static bool search_vector_pairs(const vector_pairs &vp, unsigned rs_id, unsigned &pos)
Definition: xsample03.cpp:330
static void show_help()
Definition: xsample03.cpp:81
std::string isnp_name
Definition: xsample03.cpp:104
bm::chrono_taker ::duration_map_type timing_map
Definition: xsample03.cpp:210
bool is_diag
Definition: xsample03.cpp:105
bm::rsc_sparse_vector< unsigned, sparse_vector_u32 > rsc_sparse_vector_u32
Definition: xsample03.cpp:205
static void generate_random_subset(const sparse_vector_u32 &sv, std::vector< unsigned > &vect, unsigned count)
Definition: xsample03.cpp:292
std::string sv_out_name
Definition: xsample03.cpp:100
std::string rsc_out_name
Definition: xsample03.cpp:101
static int load_snp_report(const std::string &fname, sparse_vector_u32 &sv)
Definition: xsample03.cpp:218
static void run_benchmark(const sparse_vector_u32 &sv, const rsc_sparse_vector_u32 &csv)
Definition: xsample03.cpp:349
bool is_timing
Definition: xsample03.cpp:106