fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-15 21:17:54 +08:00
Code Issues Packages Projects Releases Wiki Activity

mm2-fast-v2.22 init

Browse Source
This commit is contained in:
Saurabh
2021-11-18 21:14:38 -08:00
committed by Heng Li
parent c7cdb758db
commit ba2f19ba37
12 changed files with 3535 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
Makefile
View File

@@ -1,6 +1,42 @@
CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra
CPPFLAGS= -DHAVE_KALLOC
INCLUDES=
CPPFLAGS= -DHAVE_KALLOC #-march=native #-DALIGN_AVX -DPARALLEL_CHAINING #-DMANUAL_PROFILING
COMP_FLAG = -march=native
ifeq ($(avx2_compile), 1)
COMP_FLAG = -mavx2
endif
#CPPFLAGS= -DHAVE_KALLOC -mavx2 -DALIGN_AVX -DAPPLY_AVX2 -DPARALLEL_CHAINING #-DLISA_HASH -DUINT64 -DVECTORIZE #-DMANUAL_PROFILING
#CPPFLAGS= -DHAVE_KALLOC -mavx2 -DPARALLEL_CHAINING #-DMANUAL_PROFILING
OPT_FLAGS= -DPARALLEL_CHAINING -DALIGN_AVX -DAPPLY_AVX2
OPT_FLAGS+=$(COMP_FLAG)
ifeq ($(lhash_index), 1)
CPPFLAGS+= -DLISA_INDEX
endif
ifeq ($(lhash), 1)
OPT_FLAGS+= -DLISA_HASH -DUINT64 -DVECTORIZE
endif
ifeq ($(manual_profile), 1)
CPPFLAGS+= -DMANUAL_PROFILING
endif
#ifeq ($(use_avx2), 1)
# OPT_FLAGS+= -DAPPLY_AVX2
#endif
ifeq ($(disable_output), 1)
CPPFLAGS+= -DDISABLE_OUTPUT
endif
ifeq ($(no_opt),)
CPPFLAGS+= $(OPT_FLAGS)
endif
#INCLUDES=
#INCLUDES= -I./ext/TAL_offline/src/LISA-hash #-I./ext/TAL/src/dynamic-programming
OBJS= kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \
lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \
ksw2_ll_sse.o
@@ -10,7 +46,7 @@ LIBS= -lm -lz -lpthread
ifeq ($(arm_neon),) # if arm_neon is not defined
ifeq ($(sse2only),) # if sse2only is not defined
OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o ksw2_extd2_avx.o
else # if sse2only is defined
OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
endif

66
align.c
View File

@@ -5,6 +5,13 @@
#include "minimap.h"
#include "mmpriv.h"
#include "ksw2.h"
#include "ksw2_extd2_avx.h"
#include <x86intrin.h>
extern uint64_t avg;
extern uint64_t alignment_time;
extern void *km1;
extern uint64_t km_size;// = 500000000; // 500 MB
extern int km_top;
static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi)
{
@@ -313,6 +320,7 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
}
}
#if 0
static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
{
if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
@@ -340,7 +348,65 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
fprintf(stderr, "\n");
}
}
#endif
#if 1
static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
{
#ifdef MANUAL_PROFILING
uint64_t align_start = __rdtsc();
#endif
if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
int i;
fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr);
fputc('\n', stderr);
for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
fputc('\n', stderr);
}
if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
ksw_reset_extz(ez);
ez->zdropped = 1;
} else if (opt->flag & MM_F_SPLICE)
ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag, junc, ez);
else if (opt->q == opt->q2 && opt->e == opt->e2)
ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
else{
#if defined (ALIGN_AVX) && (defined(__AVX512BW__) || (defined(__AVX2__) && defined(APPLY_AVX2)))
#ifdef __AVX512BW__
ksw_extd2_avx512(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
#elif __AVX2__
avg = 0;
// uint64_t *ptr_km = (uint64_t *) km1;
// for(uint64_t itr = 0; itr < km_size/512; itr++){
// avg+=ptr_km[itr];
// }
//#ifdef MANUAL_PROFILING
// uint64_t align_start = __rdtsc();
//#endif
ksw_extd2_avx2(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
//#ifdef MANUAL_PROFILING
// alignment_time += (__rdtsc() - align_start);
//#endif
#endif
#else
ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
#endif
}
if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
int i;
fprintf(stderr, "score=%d, cigar=", ez->score);
for (i = 0; i < ez->n_cigar; ++i)
fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
fprintf(stderr, "\n");
}
#ifdef MANUAL_PROFILING
alignment_time += (__rdtsc() - align_start);
#endif
}
#endif
static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
{
int64_t i, off0 = mi->seq[rid].offset, off = off0 + x;

16
build_rmi.sh Executable file
View File

@@ -0,0 +1,16 @@
ref_data=$1
preset=$2
make clean && make lhash_index=1
touch temp_read.fastq
./minimap2 -ax $2 $1 temp_read.fastq >/dev/null
kv_file=$1"_"$2"_minimizers_key_value_sorted"
full_path=`readlink -f $kv_file`
cd ./ext/TAL_offline
make lisa_hash
./build-lisa-hash-index $full_path
rm ../../temp_read.fastq

360
index.c
View File

@@ -14,6 +14,21 @@
#include "mmpriv.h"
#include "kvec.h"
#include "khash.h"
#include <map>
#include <fstream>
#include <vector>
#include <algorithm>
#include <x86intrin.h>
using namespace std;
#ifdef LISA_HASH
#include "lisa_hash.h"
extern lisa_hash<uint64_t, uint64_t> *lh;
#endif
#define idx_hash(a) ((a)>>1)
#define idx_eq(a, b) ((a)>>1 == (b)>>1)
@@ -52,9 +67,43 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
if (!(mm_dbg_flag & 1)) mi->km = km_init();
return mi;
}
void mm_idx_destroy_mm_hash(mm_idx_t *mi)
{
//fprintf(stderr, "mm_destroy_hash\n");
uint32_t i;
if (mi == 0) return;
if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
if (mi->B) {
for (i = 0; i < 1U<<mi->b; ++i) {
free(mi->B[i].p);
free(mi->B[i].a.a);
kh_destroy(idx, (idxhash_t*)mi->B[i].h);
}
}
}
void mm_idx_destroy_seq(mm_idx_t *mi)
{
//fprintf(stderr, "mm_destroy_seq\n");
uint32_t i;
if (mi == 0) return;
if (mi->I) {
for (i = 0; i < mi->n_seq; ++i)
free(mi->I[i].a);
free(mi->I);
}
if (!mi->km) {
for (i = 0; i < mi->n_seq; ++i)
free(mi->seq[i].name);
free(mi->seq);
} else km_destroy(mi->km);
free(mi->B); free(mi->S); free(mi);
}
void mm_idx_destroy(mm_idx_t *mi)
{
uint32_t i;
if (mi == 0) return;
if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
@@ -96,6 +145,317 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
return &b->p[kh_val(h, k)>>32];
}
}
//Output minimap2's hash table entries
class hash_entry {
public:
uint64_t key;
uint64_t n;
uint64_t *p;
hash_entry(uint64_t k, uint64_t n_, uint64_t *p_){
key = k;
n = n_;
p = p_;
}
};
bool key_sort( hash_entry i1, hash_entry i2)
{
return (i1.key < i2.key);
}
#if 0
void mm_idx_load_key_value_lisa(const char* f_name, const mm_idx_t *mi)
{
uint64_t tic = __rdtsc();
std::vector<hash_entry> v_hash;
//ofstream f(f_name);
fprintf(stderr, "Building sorted key-val map\n");
uint32_t i,j;
uint64_t num_values = 0;
for (i = 0; i < 1U<<mi->b; ++i) {
//fprintf(stderr, "BucketID %lu \n", i);
idxhash_t *h = (idxhash_t*)mi->B[i].h;
khint_t k;
if (h == 0) continue;
for (k = 0; k < kh_end(h); ++k){
if (kh_exist(h, k)) {
uint64_t key = kh_key(h, k), bucket_id = i;
key = key>>1;
key = key<<mi->b | bucket_id;
if(kh_key(h, k)&1)
{
//print key value
//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
v_hash.push_back(hash_entry(key, kh_val(h, k), NULL));
}
else
{ // print key
uint32_t n = (uint32_t)kh_val(h, k);
//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
// for 0 to lsb 32 val
// print b->p[msb 32 of val]
v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0]));
}
}
}
}
sort(v_hash.begin(), v_hash.end(), key_sort);
fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
fprintf(stderr, "Storing hash to %s \n", f_name);
tic = __rdtsc();
int64_t itr_p = 0;
for( int i = 0; i < v_hash.size(); i++){
if(v_hash[i].p == NULL){
//f<<v_hash[i].key << " "<<1<<"\n"<<v_hash[i].n<<" \n";
lh->p[itr_p++] = v_hash[i].n;
continue;
}
//f<<v_hash[i].key << " "<<v_hash[i].n<<endl;
for(int j = 0; j < v_hash[i].n; j++){
// f<<v_hash[i].p[j]<<" ";
lh->p[itr_p++] = v_hash[i].p[j];
num_values++;
}
//f<<endl;
}
//f.close();
string size_file_name = (string) f_name + "_size";
ofstream size_f(size_file_name);
size_f<<v_hash.size()<<" "<<num_values;
size_f.close();
string prefix = (string)f_name + "_keys";
string keys_bin_file_name = prefix + ".uint64";
ofstream wf(keys_bin_file_name, ios::out | ios::binary);
wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));
wf.close();
key_list.clear();
m.clear();
v_hash.clear();
fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
}
#endif
void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi)
{
uint64_t tic = __rdtsc();
//std::map<uint64_t, vector<uint64_t>> m;
std::vector<hash_entry> v_hash;
//ofstream f(f_name);
fprintf(stderr, "Building sorted key-val map\n");
uint32_t i,j;
uint64_t num_values = 0;
for (i = 0; i < 1U<<mi->b; ++i) {
//fprintf(stderr, "BucketID %lu \n", i);
idxhash_t *h = (idxhash_t*)mi->B[i].h;
khint_t k;
if (h == 0) continue;
for (k = 0; k < kh_end(h); ++k){
if (kh_exist(h, k)) {
uint64_t key = kh_key(h, k), bucket_id = i;
key = key>>1;
key = key<<mi->b | bucket_id;
if(kh_key(h, k)&1)
{
//print key value
//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
//m[key].push_back(kh_val(h, k));
v_hash.push_back(hash_entry(key, kh_val(h, k), NULL));
}
else
{ // print key
uint32_t n = (uint32_t)kh_val(h, k);
//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
// for 0 to lsb 32 val
// print b->p[msb 32 of val]
v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0]));
}
}
}
}
sort(v_hash.begin(), v_hash.end(), key_sort);
fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
fprintf(stderr, "Storing hash to %s \n", f_name);
tic = __rdtsc();
vector<uint64_t> key_list;
vector<uint64_t> val_list;
vector<uint64_t> p_list;
/*
key_list.push_back(m.size());
for(auto k : m){
key_list.push_back(k.first);
f<<k.first << " "<<k.second.size()<<endl;
for(int j = 0; j < k.second.size(); j++){
f<<k.second[j]<<" ";
num_values++;
}
f<<endl;
}
*/
key_list.push_back(v_hash.size());
int64_t itr_p = 0;
uint64_t sum_pos = 0;
string f1_name = (string)f_name + "_pos_bin";
string f2_name = (string)f_name + "_val_bin";
ofstream f1(f1_name, ios::out | ios::binary);
ofstream f2(f2_name, ios::out | ios::binary);
for( int i = 0; i < v_hash.size(); i++){
key_list.push_back(v_hash[i].key);
if(v_hash[i].p == NULL){
//f<<v_hash[i].key << " "<<1<<"\n"<<v_hash[i].n<<" \n";
val_list.push_back(sum_pos<<32|(uint64_t)1);
sum_pos+=1;
p_list.push_back(v_hash[i].n);
num_values++;
continue;
}
//f<<v_hash[i].key << " "<<v_hash[i].n<<endl;
val_list.push_back(sum_pos<<32|(uint64_t)v_hash[i].n);
sum_pos+=v_hash[i].n;
num_values+=v_hash[i].n;
for(int j = 0; j < v_hash[i].n; j++){
//f<<v_hash[i].p[j]<<" ";
p_list.push_back(v_hash[i].p[j]);
}
// f<<endl;
}
f1.write((char*)&val_list[0], (val_list.size())*sizeof(uint64_t));
f2.write((char*)&p_list[0], (p_list.size())*sizeof(uint64_t));
f1.close();
f2.close();
fprintf(stderr, "Index sorted SoA time %lld \n", __rdtsc() - tic);
//f.close();
string size_file_name = (string) f_name + "_size";
ofstream size_f(size_file_name);
size_f<<v_hash.size()<<" "<<num_values;
size_f.close();
string prefix = (string)f_name + "_keys";
string keys_bin_file_name = prefix + ".uint64";
ofstream wf(keys_bin_file_name, ios::out | ios::binary);
wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));
wf.close();
key_list.clear();
//m.clear();
v_hash.clear();
fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
}
void mm_idx_dump_hash_1(const char* f_name, const mm_idx_t *mi)
{
uint64_t tic = __rdtsc();
std::map<uint64_t, vector<uint64_t>> m;
ofstream f(f_name);
fprintf(stderr, "Building sorted key-val map\n");
uint32_t i,j;
uint64_t num_values = 0;
for (i = 0; i < 1U<<mi->b; ++i) {
//fprintf(stderr, "BucketID %lu \n", i);
idxhash_t *h = (idxhash_t*)mi->B[i].h;
khint_t k;
if (h == 0) continue;
for (k = 0; k < kh_end(h); ++k){
if (kh_exist(h, k)) {
uint64_t key = kh_key(h, k), bucket_id = i;
key = key>>1;
key = key<<mi->b | bucket_id;
if(kh_key(h, k)&1)
{
//print key value
//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
m[key].push_back(kh_val(h, k));
}
else
{ // print key
uint32_t n = (uint32_t)kh_val(h, k);
//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
// for 0 to lsb 32 val
// print b->p[msb 32 of val]
for(j = 0; j < n; j++)
{
//fprintf(stderr, "%llu ", mi->B[i].p[(kh_val(h, k)>>32) + j]);
m[key].push_back(mi->B[i].p[(kh_val(h, k)>>32) + j]);
}
}
}
}
}
fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
fprintf(stderr, "Storing hash to %s \n", f_name);
tic = __rdtsc();
vector<uint64_t> key_list;
key_list.push_back(m.size());
for(auto k : m){
key_list.push_back(k.first);
f<<k.first << " "<<k.second.size()<<endl;
for(int j = 0; j < k.second.size(); j++){
f<<k.second[j]<<" ";
num_values++;
}
f<<endl;
}
f.close();
string size_file_name = (string) f_name + "_size";
ofstream size_f(size_file_name);
size_f<<m.size()<<" "<<num_values;
size_f.close();
string prefix = (string)f_name + "_keys";
string keys_bin_file_name = prefix + ".uint64";
ofstream wf(keys_bin_file_name, ios::out | ios::binary);
wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));
wf.close();
key_list.clear();
m.clear();
fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
}
void mm_idx_stat(const mm_idx_t *mi)
{

2319
ksw2_extd2_avx.c Normal file
View File

File diff suppressed because it is too large Load Diff

42
ksw2_extd2_avx.h Normal file
View File

@@ -0,0 +1,42 @@
/* The MIT License
Copyright (c) 2018- Dana-Farber Cancer Institute
2017-2018 Broad Institute, Inc.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Modified Copyright (C) 2021 Intel Corporation
Contacts: Saurabh Kalikar <saurabh.kalikar@intel.com>;
Vasimuddin Md <vasimuddin.md@intel.com>; Sanchit Misra <sanchit.misra@intel.com>;
Chirag Jain <chirag@iisc.ac.in>; Heng Li <hli@jimmy.harvard.edu>
*/
#include <string.h>
#include <stdio.h>
#include <assert.h>
#include "ksw2.h"
#include <immintrin.h>
#include <x86intrin.h>
#include <smmintrin.h>
#include <emmintrin.h>
void ksw_extd2_avx512(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
void ksw_extd2_avx2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);

218
lchain.c
View File

@@ -5,6 +5,11 @@
#include "mmpriv.h"
#include "kalloc.h"
#include "krmq.h"
//#include "simd_chain.h"
#include "parallel_chaining_32_bit.h"
extern uint64_t dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_t *p, int32_t *v, int32_t *t, int32_t min_cnt, int32_t min_sc, int32_t *n_u_, int32_t *n_v_)
{
@@ -88,15 +93,65 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t max_dist_x, int32_t max_dist_y, int32_t bw, float chn_pen_gap, float chn_pen_skip, int is_cdna, int n_seg)
{
int32_t dq = (int32_t)ai->y - (int32_t)aj->y, dr, dd, dg, q_span, sc;
int32_t sidi = (ai->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
int32_t sidj = (aj->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
if (dq <= 0 || dq > max_dist_x) return INT32_MIN;
dr = (int32_t)(ai->x - aj->x);
if (sidi == sidj && (dr == 0 || dq > max_dist_y)) return INT32_MIN;
uint64_t ai_x, ai_y, aj_x, aj_y;
ai_x = ai->x; ai_y = ai->y; aj_x = aj->x; aj_y = aj->y;
#ifdef CHAIN_DEBUG
int32_t sc_vect = obj.comput_sc_vectorized_avx2_caller(ai_x, ai_y, aj_x, aj_y, aj->y>>32&0xff);
#endif
//if (sc_vect == 0) return INT32_MIN;
//else
//return sc_vect;
//fprintf(stderr, "%lld %lld %lld %lld \n", ai_x, ai_y, aj_x, aj_y);
//fprintf(stderr, "%lld %lld %lld %f %f %d %d\n", max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
int32_t dq = (int32_t)ai_y - (int32_t)aj_y, dr, dd, dg, q_span, sc;
int32_t sidi = (ai_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
int32_t sidj = (aj_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
if (dq <= 0 || dq > max_dist_x) {
#ifdef CHAIN_DEBUG
if(INT32_MIN != sc_vect){
//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
}
#endif
return INT32_MIN;
}
dr = (int32_t)(ai_x - aj_x);
if (sidi == sidj && (dr == 0 || dq > max_dist_y)) {
#ifdef CHAIN_DEBUG
if(INT32_MIN != sc_vect){
//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
}
#endif
return INT32_MIN;
}
dd = dr > dq? dr - dq : dq - dr;
if (sidi == sidj && dd > bw) return INT32_MIN;
if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) return INT32_MIN;
if (sidi == sidj && dd > bw) {
#ifdef CHAIN_DEBUG
if(INT32_MIN != sc_vect){
//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
}
#endif
return INT32_MIN;
}
if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) {
#ifdef CHAIN_DEBUG
if(INT32_MIN != sc_vect){
//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
}
#endif
return INT32_MIN;
}
dg = dr < dq? dr : dq;
q_span = aj->y>>32&0xff;
sc = q_span < dg? q_span : dg;
@@ -110,6 +165,13 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
else sc -= (int)(lin_pen + .5f * log_pen);
} else sc -= (int)(lin_pen + .5f * log_pen);
}
#ifdef CHAIN_DEBUG
if(sc != sc_vect ){
//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
fprintf(stderr, "outer: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
}
#endif
return sc;
}
@@ -124,10 +186,18 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
int is_cdna, int n_seg, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
{ // TODO: make sure this works when n has more than 32 bits
int32_t *f, *t, *v, n_u, n_v, mmax_f = 0;
///fprintf(stderr, "chaining called\n");
#ifdef MANUAL_PROFILING
uint64_t align_start = __rdtsc();
#endif
int32_t *f, *t, *v, *v_1, *p_1, n_u, n_v, mmax_f = 0;
int64_t *p, i, j, max_ii, st = 0, n_iter = 0;
uint64_t *u;
uint32_t* f_1;
if (_u) *_u = 0, *n_u_ = 0;
if (n == 0 || a == 0) {
kfree(km, a);
@@ -136,20 +206,67 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
if (max_dist_x < bw) max_dist_x = bw;
if (max_dist_y < bw && !is_cdna) max_dist_y = bw;
KMALLOC(km, p, n);
KMALLOC(km, p_1, n);
KMALLOC(km, f, n);
KMALLOC(km, f_1, n);
KMALLOC(km, v, n);
KMALLOC(km, v_1, n);
KCALLOC(km, t, n);
#ifdef PARALLEL_CHAINING
// Parallel chaining data-structures
anchor_t* anchors = (anchor_t*)malloc(n* sizeof(anchor_t));
for (i = 0; i < n; ++i) {
uint64_t ri = a[i].x;
int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!!
anchors[i].r = ri;
anchors[i].q = qi;
anchors[i].l = q_span;
}
num_bits_t *anchor_r, *anchor_q, *anchor_l;
create_SoA_Anchors_32_bit(anchors, n, anchor_r, anchor_q, anchor_l);
//dp_chain obj(max_dist_x, max_dist_y, bw, max_skip, max_iter, 0, is_cdna, n_seg);
dp_chain obj(max_dist_x, max_dist_y, bw, max_skip, max_iter, min_cnt, min_sc, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
obj.mm_dp_vectorized(n, &anchors[0], anchor_r, anchor_q, anchor_l, f_1, p_1, v_1, max_dist_x, max_dist_y, NULL, NULL);
// -16 is due to extra padding at the start of arrays
anchor_r -= 16; anchor_q -= 16; anchor_l -= 16;
free(anchor_r);
free(anchor_q);
free(anchor_l);
free(anchors);
for(int i = 0; i < n; i++){
// if(f[i] != f_1[i] || p[i] != p_1[i] || v[i] !=v_1[i])
// {
// fprintf(stderr, "i:%d %d %d %d %d %d %d\n",i, f[i], f_1[i], p[i], p_1[i], v[i], v_1[i] );
// }
#if 1
f[i] = f_1[i];
p[i] = p_1[i];
v[i] = v_1[i];
#endif
}
//
#else
// fill the score and backtrack arrays
for (i = 0, max_ii = -1; i < n; ++i) {
int64_t max_j = -1, end_j;
int32_t max_f = a[i].y>>32&0xff, n_skip = 0;
while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist_x)) ++st;
if (i - st > max_iter) st = i - max_iter;
int my_cnt = 0;
for (j = i - 1; j >= st; --j) {
int32_t sc;
sc = comput_sc(&a[i], &a[j], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
++n_iter;
// if(i == 177){
//fprintf(stderr, "args: %d %d %d %d %d\n", a[i].x, a[i].y, a[j].x, a[j].y, a[j].y>>32&0xff);
//fprintf(stderr, "j_th %d score: %d\n", ++my_cnt, sc);
// }
if (sc == INT32_MIN) continue;
sc += f[j];
if (sc > max_f) {
@@ -162,32 +279,72 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
if (p[j] >= 0) t[p[j]] = i;
}
end_j = j;
int debug_iter = 2057329;
//if (i == debug_iter) fprintf(stderr, "mm2 -- endj: %d max_ii: %d max_f: %d \n", end_j, max_ii, max_f);
#if 1
if (max_ii < 0 || a[i].x - a[max_ii].x > (int64_t)max_dist_x) {
int32_t max = INT32_MIN;
max_ii = -1;
for (j = i - 1; j >= st; --j)
if (max < f[j]) max = f[j], max_ii = j;
for (j = i - 1; j >= st; --j) {
if (max < (int32_t)f[j]) max = f[j], max_ii = j;
}
}
#endif
#if 1
if (max_ii >= 0 && max_ii < end_j) {
int32_t tmp;
tmp = comput_sc(&a[i], &a[max_ii], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
// if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d \n", end_j, max_ii, max_f, tmp);
if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) {
// if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d \n", end_j, max_ii, max_f, tmp);
max_f = tmp + f[max_ii], max_j = max_ii;
// if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d sum : %d \n", end_j, max_ii, max_f, tmp, tmp + f[max_ii]);
}
}
#endif
f[i] = max_f, p[i] = max_j;
v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak
#if 1
if (max_ii < 0 || (a[i].x - a[max_ii].x <= (int64_t)max_dist_x && f[max_ii] < f[i]))
max_ii = i;
if (mmax_f < max_f) mmax_f = max_f;
#endif
}
#endif
#ifdef CHAIN_DEBUG
for(int i = 0; i < n; i++){
if(f[i] != f_1[i] || p[i] != p_1[i] || v[i] !=v_1[i])
{
fprintf(stderr, "i:%d %d %d %d %d %d %d\n",i, f[i], f_1[i], p[i], p_1[i], v[i], v_1[i] );
}
#if 0
f[i] = f_1[i];
p[i] = p_1[i];
v[i] = v_1[i];
#endif
}
#endif
u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, &n_u, &n_v);
*n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here
kfree(km, p); kfree(km, f); kfree(km, t);
kfree(km, p); kfree(km, p_1); kfree(km, f); kfree(km, f_1); kfree(km, t); kfree(km, v_1);
if (n_u == 0) {
kfree(km, a); kfree(km, v);
return 0;
}
#ifdef MANUAL_PROFILING
dp_time += __rdtsc() - align_start;
#endif
return compact_a(km, n_u, u, n_v, v, a);
}
@@ -225,6 +382,11 @@ static inline int32_t comput_sc_simple(const mm128_t *ai, const mm128_t *aj, flo
mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
{
#ifdef MANUAL_PROFILING
uint64_t start = __rdtsc();
#endif
uint64_t tim;
//fprintf(stderr, "rmq call \n");
int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0;
int64_t *p, i, i0, st = 0, st_inner = 0, n_iter = 0;
uint64_t *u;
@@ -252,6 +414,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
int32_t q_span = a[i].y>>32&0xff, max_f = q_span;
lc_elem_t s, *q, *r, lo, hi;
// add in-range anchors
#ifdef MANUAL_PROFILING_RMQ
tim = __rdtsc();
#endif
if (i0 < i && a[i0].x != a[i].x) {
int64_t j;
for (j = i0; j < i; ++j) {
@@ -266,7 +431,13 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
}
i0 = i;
}
#ifdef MANUAL_PROFILING_RMQ
rmq_t1 += __rdtsc() - tim;
#endif
// get rid of active chains out of range
#ifdef MANUAL_PROFILING_RMQ
tim = __rdtsc();
#endif
while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist || krmq_size(head, root) > cap_rmq_size)) {
s.y = (int32_t)a[st].y, s.i = st;
if ((q = krmq_find(lc_elem, root, &s, 0)) != 0) {
@@ -275,6 +446,12 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
}
++st;
}
#ifdef MANUAL_PROFILING_RMQ
rmq_t2 += __rdtsc() - tim;
#endif
#ifdef MANUAL_PROFILING_RMQ
tim = __rdtsc();
#endif
if (max_dist_inner > 0) { // similar to the block above, but applied to the inner tree
while (st_inner < i && (a[i].x>>32 != a[st_inner].x>>32 || a[i].x > a[st_inner].x + max_dist_inner || krmq_size(head, root_inner) > cap_rmq_size)) {
s.y = (int32_t)a[st_inner].y, s.i = st_inner;
@@ -285,6 +462,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
++st_inner;
}
}
#ifdef MANUAL_PROFILING_RMQ
rmq_t3 += __rdtsc() - tim;
#endif
// RMQ
lo.i = INT32_MAX, lo.y = (int32_t)a[i].y - max_dist;
hi.i = 0, hi.y = (int32_t)a[i].y;
@@ -304,6 +484,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
krmq_itr_t(lc_elem) itr;
krmq_itr_find(lc_elem, root_inner, lo, &itr);
while ((q = krmq_at(&itr)) != 0) {
#ifdef MANUAL_PROFILING_RMQ
tim = __rdtsc();
#endif
if (q->y < (int32_t)a[i].y - max_dist_inner) break;
++n_rmq_iter;
j = q->i;
@@ -319,11 +502,15 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
if (p[j] >= 0) t[p[j]] = i;
}
if (!krmq_itr_prev(lc_elem, &itr)) break;
#ifdef MANUAL_PROFILING_RMQ
rmq_t4 += __rdtsc() - tim;
#endif
}
n_iter += n_rmq_iter;
}
}
}
// set max
assert(max_j < 0 || (a[max_j].x < a[i].x && (int32_t)a[max_j].y < (int32_t)a[i].y));
f[i] = max_f, p[i] = max_j;
@@ -340,5 +527,8 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
kfree(km, a); kfree(km, v);
return 0;
}
#ifdef MANUAL_PROFILING
rmq_time += __rdtsc() - start;
#endif
return compact_a(km, n_u, u, n_v, v, a);
}

307
lisa_hash.h Normal file
View File

@@ -0,0 +1,307 @@
/*************************************************************************************
MIT License
Copyright (c) 2020 Intel Labs
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Authors: Saurabh Kalikar <saurabh.kalikar@intel.com>; Sanchit Misra <sanchit.misra@intel.com>;
*****************************************************************************************/
#include<map>
#include<vector>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include "RMI.h"
using namespace std;
template<typename rmi_key_t, typename rmi_val_t>
class lisa_hash{
public:
rmi_val_t* p;
rmi_val_t* p_bin;
private:
rmi_key_t* keys;
uint64_t* values_enc;
uint64_t* values_enc_bin;
uint64_t keys_size;
uint64_t p_size;
RMI<rmi_key_t> *rmi;
void mem_alloc(uint64_t hash_size, uint64_t p_size){
// Removed for memory optimization
#if 0
keys = (rmi_key_t*)malloc((1+hash_size)*sizeof(rmi_key_t));
keys = keys + 1; // keys[-1] stores total number of keys
#endif
values_enc = (uint64_t*)malloc(hash_size*sizeof(uint64_t));
p = (rmi_val_t*)malloc(p_size*sizeof(rmi_val_t));
//values_enc_bin = (uint64_t*)malloc(hash_size*sizeof(uint64_t));
//p_bin = (rmi_val_t*)malloc(p_size*sizeof(rmi_val_t));
keys_size = hash_size;
// keys[-1] = keys_size;
this->p_size = p_size;
fprintf(stderr, "Memory allocated %lld \n", p_size);
}
//This function is no longer used.
void store_sorted_keys(string binFileName){
ofstream wf(binFileName, ios::out | ios::binary);
wf.write((char*)&keys[-1], (1+keys_size)*sizeof(rmi_key_t));
wf.close();
rmi_key_t *temp_keys = (rmi_key_t*)malloc((1 + keys_size)*sizeof(rmi_key_t));
ifstream rf(binFileName, ios::out | ios::binary);
rf.read((char*)&temp_keys[0], (1 + keys_size)*sizeof(rmi_key_t));
if(temp_keys[0] == keys[-1]){
for(int i = 0; i < keys_size; i++){
if(temp_keys[i+1] != keys[i]){
cout<<"Error: "<<temp_keys[i+1]<< " " << keys[i]<<endl;
break;
}
}
}
else{
cout<<"File writing error!!\n"<<temp_keys[0] << " "<<keys[-1] << keys_size;
}
free(temp_keys);
}
void set_key_value(uint64_t i, rmi_key_t key, uint64_t p_pos, uint64_t p_size){
// keys[i] = key;//<<1;
values_enc[i] = p_pos << 32 | p_size;
}
void get_val(uint64_t i, uint64_t &pos, uint64_t &p_size){
pos = values_enc[i] >> 32;
p_size = (uint32_t) values_enc[i];
}
void load_bin(string inputFile){
fprintf(stderr, "Loading from bin\n");
string f1_name = (string) inputFile + "_pos_bin";
string f2_name = (string) inputFile + "_val_bin";
ifstream instream_f1(f1_name, ifstream::binary);
ifstream instream_f2(f2_name, ifstream::binary);
instream_f1.seekg(0);
instream_f1.read((char*)&values_enc[0], keys_size*sizeof(uint64_t));
instream_f1.close();
instream_f2.seekg(0);
instream_f2.read((char*)&p[0], p_size*sizeof(uint64_t));
instream_f2.close();
/*
for(uint64_t i = 0; i < p_size; i++){
if(p[i] != p_bin[i]){
fprintf(stderr, "Error!! %lld %lld\n",p[i], p_bin[i] );
}
}
for(uint64_t i = 0; i < keys_size; i++){
if(values_enc[i] != values_enc_bin[i]){
fprintf(stderr, "Error!! %lld %lld\n", values_enc[i], values_enc_bin[i]);
}
}
*/
}
void load(string inputFile){
ifstream f(inputFile);
rmi_key_t key;
rmi_val_t val;
uint64_t n;
int64_t offset = 0;
int64_t i = 0;
while(f>>key){
f>>n;
//values_enc[++key_size] = val_size << 32 | n;
for(int j = 0; j < n; j++){
f>>val;
p[offset + j] = val;
}
set_key_value(i, key, offset, n);
offset +=n;
i++;
}
}
public:
lisa_hash(string inputFile, char* rmi_prefix, long leaf = 0){
uint64_t start_time, load_time, rmi_building_time, rmi_object;
start_time = __rdtsc();
uint64_t val_count = 0;
ifstream f_size(inputFile+"_size");
f_size>>keys_size;
f_size>>p_size;
mem_alloc(keys_size, p_size);
fprintf(stderr, "Num_keys: %lld, num_values = %lld", keys_size, p_size);
//load(inputFile);
load_bin(inputFile);
string prefix = inputFile + "_keys";
#ifdef UINT64
string keys_bin_file_name = prefix + ".uint64";
#else
string keys_bin_file_name = prefix + ".f64";
#endif
ifstream rf(keys_bin_file_name, ios::out | ios::binary);
if(!rf.good()){
cout<<"Error: Binary file with keys not found!!\n";
//call store keys
//store_sorted_keys(keys_bin_file_name);
exit(0);
}
string keys_rmi_file = prefix + ".rmi_PARAMETERS";
ifstream rmi_f(keys_rmi_file, ios::out | ios::binary);
if(!rmi_f.good() || leaf != 0){
cout<<"rmi file not found: "<< keys_rmi_file <<endl;
if(leaf == 0)
{
cout<<"Number of rmi leaf nodes are not provided for "<<keys_size<<" keys\n";
cout<<"Using default number of leaf nodes: "<<keys_size/32<<"\n";
leaf = keys_size/32;
}
load_time = __rdtsc() - start_time;
start_time = __rdtsc();
fprintf(stderr,"TIMER LOG: kay-val files loading time- %lld\n", load_time);
//string script = "./scripts/build_rmi.sh";
string script = "./scripts/build-rmi.linear_spline.linear.sh";
string keys_path = " " + keys_bin_file_name;
string rmi_path = " " + keys_rmi_file;
string num_leaf = " " + std::to_string(leaf);
#ifdef UINT64
string key_type = " UINT64";
#else
string key_type = " F64";
#endif
string cmd = script + keys_path + rmi_path + num_leaf + key_type;
system(cmd.c_str());
rmi_building_time = __rdtsc() - start_time;
fprintf(stderr,"TIMER LOG: rmi building time- %lld\n", rmi_building_time);
}
rmi = new RMI<rmi_key_t>(&prefix[0]);
}
rmi_val_t* get_hash_value(rmi_key_t key, int*n){
uint64_t index = rmi->lookup(key);
if(index == -1){
cout<<"Key not found\n";
return NULL;
}
uint64_t pos, p_size;
get_val(index, pos, p_size);
*n = p_size;
return &p[pos];
}
rmi_val_t* get_hash_values_batched(rmi_key_t *keys, uint64_t num_keys, int* &num_values){
int64_t *pos = (int64_t*) malloc(num_keys*sizeof(int64_t));
rmi->lookup_batched(keys, num_keys, &pos[0]);
num_values = (int*) malloc(num_keys*sizeof(int));
rmi_val_t **p_ptrs = (rmi_val_t**) malloc(num_keys*sizeof(rmi_val_t*));
uint64_t total_num_values = 0;
for(int i = 0; i < num_keys; i++){
if(pos[i] == -1){
num_values[i] = 0;
continue;
}
num_values[i] = (uint32_t) values_enc[pos[i]];
total_num_values+= num_values[i];
p_ptrs[i] = &p[values_enc[pos[i]] >> 32];
}
free(pos);
rmi_val_t* ret_values = (rmi_val_t*) malloc(total_num_values*sizeof(rmi_val_t));
uint64_t cnt= 0;
for(int i = 0; i < num_keys; i++){
rmi_val_t *value_start_ptr = p_ptrs[i];
uint64_t numhit = num_values[i];
for(int j = 0; j < numhit; j++){
ret_values[cnt++] = value_start_ptr[j];
}
}
free(p_ptrs);
return ret_values;
}
void mm_idx_get_batched(uint64_t* &minimizers, uint64_t num_minimizers, int64_t* &pos, uint64_t** &p_ptrs, int* &num_hits){
rmi->lookup_batched(minimizers, num_minimizers, &pos[0]);
for(int i = 0; i < num_minimizers; i++){
int64_t p_i = pos[i];
if(p_i < 0 || p_i > keys_size )
num_hits[i] = 0;
else
num_hits[i] = (uint32_t) values_enc[p_i];
p_ptrs[i] = p + (values_enc[p_i] >> 32);
}
}
~lisa_hash(){
delete rmi;
free(values_enc);
free(p);
}
};

100
main.c
View File

@@ -1,5 +1,3 @@
/* Fast contrib: v.22 */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
@@ -8,6 +6,73 @@
#include "minimap.h"
#include "mmpriv.h"
#include "ketopt.h"
#include <x86intrin.h>
#include <immintrin.h>
#include <sys/time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <string>
#include <map>
#include <errno.h>
#include "bseq.h"
#include "minimap.h"
#include "mmpriv.h"
#include "ketopt.h"
//#include "profile.h"
#include <stdint.h>
#include <unistd.h>
#include <x86intrin.h>
using namespace std;
uint64_t avg;
uint64_t minimizer_lookup_time, alignment_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
#ifdef LISA_HASH
#include "lisa_hash.h"
lisa_hash<uint64_t, uint64_t> *lh;
#endif
// New memory allocation approach for alignment optimizations
//
void *km1;
uint64_t km_size = 500000000; // 500 MB
int km_top;
/*
void *kcalloc_(void* km, int count, int size)
{
assert(count*size < km_size);
km_top += count*size + 1024;
memset(km, 0, count * size);
// printf("km_top: %d\n", km_top);
return km;
}
void *kmalloc_(void* km, int count) {
if(km_top + count >= km_size)
printf("count: %d\n", count);
assert(km_top + count < km_size);
void *mem = (void*) ((int8_t*) km + km_top);
km_top += count + 1024;
// printf("km_top: %d\n", km_top);
return mem;
}
void kfree_all() { km_top = 0;}
*/
// Memory for alignment end
#ifndef __rdtsc
#ifdef _rdtsc
#define __rdtsc _rdtsc
#else
#define __rdtsc __builtin_ia32_rdtsc
#endif
#endif
#define MM_VERSION "2.22-r1101"
@@ -119,6 +184,9 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const
int main(int argc, char *argv[])
{
// Memory allocation for alignment optimizations
//km1 = calloc(km_size, 1); // 10 MB init contg. alloc
const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:";
ketopt_t o = KETOPT_INIT;
mm_mapopt_t opt;
@@ -133,9 +201,11 @@ int main(int argc, char *argv[])
liftrlimit();
mm_realtime0 = realtime();
mm_set_opt(0, &ipt, &opt);
string preset_arg = "";
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
if (c == 'x') {
preset_arg += (string) o.arg;
if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
return 1;
@@ -368,6 +438,7 @@ int main(int argc, char *argv[])
fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
return 1;
}
preset_arg = (string)argv[o.ind] + "_" + preset_arg + "_minimizers_key_value_sorted";
idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
if (idx_rdr == 0) {
fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno));
@@ -410,6 +481,9 @@ int main(int argc, char *argv[])
__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
if (mm_verbose >= 3) mm_idx_stat(mi);
#ifdef LISA_INDEX
mm_idx_dump_hash(preset_arg.c_str(), mi);
#endif
if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
if (alt_list) mm_idx_alt_read(mi, alt_list);
if (argc - (o.ind + 1) == 0) {
@@ -417,6 +491,16 @@ int main(int argc, char *argv[])
continue; // no query files
}
ret = 0;
#ifdef LISA_HASH
fprintf(stderr, "Using LISA_HASH..\n");
mm_idx_destroy_mm_hash(mi);
char* prefix;
lh = new lisa_hash<uint64_t, uint64_t>(preset_arg, prefix);
fprintf(stderr, "Loading done.\n");
// total_time = __rdtsc();
// fprintf(stderr, "\nIndexing Real time: %.3f sec;\n", realtime() - mapping_time);
#endif
mm_realtime0 = realtime();
if (!(opt.flag & MM_F_FRAG_MODE)) {
for (i = o.ind + 1; i < argc; ++i) {
ret = mm_map_file(mi, argv[i], &opt, n_threads);
@@ -425,12 +509,17 @@ int main(int argc, char *argv[])
} else {
ret = mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
}
mm_idx_destroy(mi);
//mm_idx_destroy(mi);
if (ret < 0) {
fprintf(stderr, "ERROR: failed to map the query file\n");
exit(EXIT_FAILURE);
}
}
#ifdef LISA_HASH
mm_idx_destroy_seq(mi);
#else
mm_idx_destroy(mi);
#endif
n_parts = idx_rdr->n_parts;
mm_idx_reader_close(idx_rdr);
@@ -449,5 +538,10 @@ int main(int argc, char *argv[])
fprintf(stderr, " %s", argv[i]);
fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
}
fprintf(stderr, "minimizer-lookup: %lld dp: %lld rmq: %lld rmq_t1: %lld rmq_t2: %lld rmq_t3: %lld rmq_t4: %lld alignment: %lld %lld\n", minimizer_lookup_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4, alignment_time, avg);
#ifdef LISA_HASH
delete lh;
#endif
return 0;
}

20
map.c
View File

@@ -10,6 +10,8 @@
#include "bseq.h"
#include "khash.h"
extern uint64_t rmq_time;
struct mm_tbuf_s {
void *km;
int rep_len, frag_gap;
@@ -276,7 +278,12 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
if (opt->flag & MM_F_RMQ) {
a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
// a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
// opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
} else {
//fprintf(stderr, "dp call - n_a = %lld\n", n_a);
a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
}
@@ -284,12 +291,23 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y;
if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) {
#ifdef MANUAL_PROFILING
// uint64_t tim = __rdtsc();
#endif
// fprintf(stderr, "pre: rmq rechain call - n_a = %lld n_regs = %lld\n",n_a, n_regs0);
int32_t i;
int64_t prev_n_a = n_a;
for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i];
kfree(b->km, u);
radix_sort_128x(a, a + n_a);
// fprintf(stderr, "post: rmq rechain call - prev_n_a = %lld n_a = %lld n_regs = %lld\n",prev_n_a, n_a, n_regs0);
// a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
// opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
#ifdef MANUAL_PROFILING
// rmq_time += __rdtsc() - tim;
#endif
}
} else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads
int rechain = 0;
@@ -563,6 +581,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init();
for (k = 0; k < s->n_frag; ++k) {
int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
#ifndef DISABLE_OUTPUT
for (i = seg_st; i < seg_en; ++i) {
mm_bseq1_t *t = &s->seq[i];
if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
@@ -597,6 +616,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
mm_err_puts(p->str.s);
}
}
#endif
for (i = seg_st; i < seg_en; ++i) {
for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p);
free(s->reg[i]);

20
minimap.h
View File

@@ -285,6 +285,13 @@ mm_idx_t *mm_idx_load(FILE *fp);
*/
void mm_idx_dump(FILE *fp, const mm_idx_t *mi);
/**
* Store hash table from minimap2 index into a file
* @param f_name File name for output file
* @param mi minimap2 index
*/
void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi);
/**
* Create an index from strings in memory
*
@@ -313,6 +320,19 @@ void mm_idx_stat(const mm_idx_t *idx);
* @param r minimap2 index
*/
void mm_idx_destroy(mm_idx_t *mi);
/**
* Destroy/deallocate an hash table index
*
* @param r minimap2 index
*/
void mm_idx_destroy_mm_hash(mm_idx_t *mi);
/**
* Destroy/deallocate target sequences
*
* @param r minimap2 index
*/
void mm_idx_destroy_seq(mm_idx_t *mi);
/**
* Initialize a thread-local buffer for mapping

45
seed.c
View File

@@ -1,9 +1,40 @@
#include "mmpriv.h"
#include "kalloc.h"
#include "ksort.h"
#include <stdlib.h>
#include<algorithm>
#include <x86intrin.h>
#ifdef LISA_HASH
#include "lisa_hash.h"
extern lisa_hash<uint64_t, uint64_t> *lh;
#endif
extern uint64_t minimizer_lookup_time;
mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_)
{
#ifdef MANUAL_PROFILING
uint64_t lookup_start = __rdtsc();
#endif
#ifdef LISA_HASH
//-----------------------------------
uint64_t** cr_batch = (uint64_t**) malloc((mv->n)*sizeof(uint64_t*));
int* t_batch = (int*)malloc((mv->n)*sizeof(int));
uint64_t* minimizers = (uint64_t*) malloc((mv->n)*sizeof(uint64_t));
int64_t* lisa_pos = (int64_t*) malloc((max(32, (int)mv->n))* sizeof(int64_t));
for (size_t i = 0; i < mv->n; i++) {
mm128_t *p = &mv->a[i];
minimizers[i] = p->x>>8;
}
lh->mm_idx_get_batched(minimizers, mv->n, lisa_pos, cr_batch, t_batch);
//-----------------------------------
#endif
mm_seed_t *m;
size_t i;
int32_t k;
@@ -14,7 +45,12 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv,
mm128_t *p = &mv->a[i];
uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
int t;
#ifdef LISA_HASH
t = t_batch[i];
cr = cr_batch[i];
#else
cr = mm_idx_get(mi, p->x>>8, &t);
#endif
if (t == 0) continue;
q = &m[k++];
q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32;
@@ -22,7 +58,16 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv,
if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1;
if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1;
}
#ifdef LISA_HASH
free(cr_batch);
free(t_batch);
free(minimizers);
free(lisa_pos);
#endif
*n_m_ = k;
#ifdef MANUAL_PROFILING
minimizer_lookup_time += __rdtsc() - lookup_start;
#endif
return m;
}