diff --git a/.gitmodules b/.gitmodules index a80f848..ab9bf2b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "lib/simde"] path = lib/simde url = https://github.com/nemequ/simde.git +[submodule "ext/TAL"] + path = ext/TAL + url = https://github.com/IntelLabs/Trans-Omics-Acceleration-Library.git diff --git a/Makefile b/Makefile index 18622f5..a5b1806 100644 --- a/Makefile +++ b/Makefile @@ -1,14 +1,44 @@ -CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra -CPPFLAGS= -DHAVE_KALLOC -INCLUDES= +CFLAGS= -Wall -O2 -Wc++-compat #-Wextra +CPPFLAGS= -DHAVE_KALLOC -march=native + +OPT_FLAGS= -DVECTORIZED_CHAINING -DALIGN_AVX + +ifeq ($(lhash), 1) + OPT_FLAGS+= -DLISA_HASH -DUINT64 -DVECTORIZE +endif + + +ifeq ($(manual_profile), 1) + CPPFLAGS+= -DMANUAL_PROFILING +endif + +ifeq ($(arch), avx512) + CPPFLAGS+= -xCORE-AVX512 +endif + +ifeq ($(disable_output), 1) + CPPFLAGS+= -DDISABLE_OUTPUT +endif + +ifeq ($(no_opt),) + CPPFLAGS+= $(OPT_FLAGS) +endif + +INCLUDES= -I./ext/TAL/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 chain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o ksw2_ll_sse.o PROG= minimap2 PROG_EXTRA= sdust minimap2-lite -LIBS= -lm -lz -lpthread +LIBS= -lm -lz -lpthread + +CC=$(CXX) +ifeq ($(CC), g++) + CC=g++ -std=c++11 +endif 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 diff --git a/README.md b/README.md index 9b4dfe9..30e47ad 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,32 @@ +## mm2-fast +### Introduction +mm2-fast is an accelerated implementation of minimap2 on modern CPUs. mm2-fast accelerates all three major modules of minimap2: Seeding, Chaining, and Alignment, achieving upto 3.5x speedup over minimap2. +mm2-fast is a drop-in replacement of minimap2, providing the same functionality with the exact same output. +In the current version, all the modules are optimized using **AVX-512** vectorization. For non-AVX-512 machines, mm2-fast runs as minimap2 code. + +### Usage +```sh +git clone --recursive https://github.com/lh3/minimap2.git -b fast-contrib mm2-fast +cd mm2-fast + +# Compile and run mm2-fast (without seeding module optimizations). +make clean && make +./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam + +# Compile and run mm2-fast (with all three optimized modules) +1. Build learned hash table index for optimized seeding module + Pre-requisite: Install "Rust" and add path to .bashrc file. For Rust installation, visit https://rustup.rs/ + ./build_rmi.sh ref.fa map-ont ##takes two arguments: 1. path-to-reference-seq-file 2. preset + +2. Compile and run +make clean && make lhash=1 +./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam + +# Compile with all optimizations disabled (runs as minimap2) +make clean && make no_opt=1 +``` +The original README content of minimap2 follows. + [![GitHub Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=social&logo=github&label=Download)](https://github.com/lh3/minimap2/releases) [![BioConda Install](https://img.shields.io/conda/dn/bioconda/minimap2.svg?style=flag&label=BioConda%20install)](https://anaconda.org/bioconda/minimap2) [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy) diff --git a/align.c b/align.c index 6028abb..c95512c 100644 --- a/align.c +++ b/align.c @@ -1,3 +1,33 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ + #include #include #include @@ -5,7 +35,9 @@ #include "minimap.h" #include "mmpriv.h" #include "ksw2.h" - +#include "ksw2_extd2_avx.h" +#include +extern uint64_t alignment_time; static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi) { int i, j; @@ -312,6 +344,10 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // 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); @@ -327,8 +363,17 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint 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 + else{ +#ifdef ALIGN_AVX +#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); +#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 +#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); @@ -336,6 +381,9 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint 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 } static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x) diff --git a/build_rmi.sh b/build_rmi.sh new file mode 100755 index 0000000..d162799 --- /dev/null +++ b/build_rmi.sh @@ -0,0 +1,17 @@ +ref_data=$1 +preset=$2 + +make clean && make no_opt=1 +touch temp_read.fastq +./minimap2 -ax $2 $1 temp_read.fastq -Z 1 >/dev/null + + +kv_file=$1"_"$2"_minimizers_key_value_sorted" + +full_path=`readlink -f $kv_file` + +cd ./ext/TAL +make lisa_hash +./build-lisa-hash-index $kv_file + +rm ../../temp_read.fastq diff --git a/chain.c b/chain.c index a2f7ac5..94e9be0 100644 --- a/chain.c +++ b/chain.c @@ -1,10 +1,39 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ #include #include #include #include "minimap.h" #include "mmpriv.h" #include "kalloc.h" - +#include "parallel_chaining_32_bit.h" static const char LogTable256[256] = { #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, @@ -19,9 +48,11 @@ static inline int ilog2_32(uint32_t v) return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v]; } + mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale, int is_cdna, int n_segs, 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 k, *f, *p, *t, *v, n_u, n_v; + int32_t k, *p, *t, *v, n_u, n_v; + uint32_t *f; int64_t i, j, st = 0; uint64_t *u, *u2, sum_qspan = 0; float avg_qspan; @@ -32,11 +63,33 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m kfree(km, a); return 0; } - f = (int32_t*)kmalloc(km, n * 4); + f = (uint32_t*)kmalloc(km, n * 4); p = (int32_t*)kmalloc(km, n * 4); t = (int32_t*)kmalloc(km, n * 4); v = (int32_t*)kmalloc(km, n * 4); memset(t, 0, n * 4); +#if VECTORIZED_CHAINING + 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, gap_scale, is_cdna, n_segs); + + obj.mm_dp_vectorized(n, &anchors[0], anchor_r, anchor_q, anchor_l, f, p, v, 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); +#else for (i = 0; i < n; ++i) sum_qspan += a[i].y>>32&0xff; avg_qspan = (float)sum_qspan / n; @@ -85,7 +138,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m 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 } - +#endif // find the ending positions of chains memset(t, 0, n * 4); for (i = 0; i < n; ++i) diff --git a/ext/TAL b/ext/TAL new file mode 160000 index 0000000..75ce461 --- /dev/null +++ b/ext/TAL @@ -0,0 +1 @@ +Subproject commit 75ce461846bbadeef0c22fb727394349e4b1d89a diff --git a/index.c b/index.c index b24a505..cd78061 100644 --- a/index.c +++ b/index.c @@ -1,4 +1,37 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ #include +#include +#include +#include +using namespace std; #include #if defined(WIN32) || defined(_WIN32) #include // for open(2) @@ -53,6 +86,37 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag) return mi; } + +void mm_idx_destroy_mm_hash(mm_idx_t *mi) +{ + 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<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) +{ + uint32_t i; + 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; @@ -97,6 +161,81 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n) } } +//Output minimap2's hash table entries +void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi) +{ + std::map> m; + + ofstream f(f_name); + fprintf(stderr, "Building sorted key-val map"); + + uint32_t i,j; + uint64_t num_values = 0; + for (i = 0; i < 1U<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<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, "Storing hash to %s \n", f_name); + vector key_list; + key_list.push_back(m.size()); + for(auto k : m){ + key_list.push_back(k.first); + f<; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ +#include "ksw2_extd2_avx.h" + +#ifdef __AVX512BW__ +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) +{ + __m512i ind512_srli = _mm512_set1_epi8(15); + __mmask8 msk_srli = 0x0003; + __m512i bt32_ = _mm512_setr_epi32(0,0,0,0,4,4,4,4,8,8,8,8,12,12,12,12); + + int8_t index[64] __attribute((aligned(64))); + + for (int i=0; i<64; i++) + index[i] = i%16 - 1; + + index[0] = 15; + index[16] = 31; + index[32] = 47; + index[48] = 63; + + __m512i shf512a, shf512b, slli512; + __m512i ind512_slli = _mm512_load_si512((__m512i*) index); + __mmask8 mska = 0x90; + __mmask64 mskb = 0x0001000100010000; + __mmask64 mskc = 0x1; + __mmask64 mskc_ar[4] = {0x1, 0x10000, 0x100000000, 0x1000000000000}; + + + + #define __dp_code_block1_pcl \ + /*__mmask64 mskc_ = mskc; */ \ + /*if (t == st_) */ \ + __mmask64 mskc_ = (t == st_) ? mskc_ar[(st0 - t*64)/16]:mskc; \ + z = _mm512_load_si512(&s[t]); \ + xt1 = _mm512_load_si512(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(xt1, 15); */ /* tmp <- x[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)x)[t*64 + 63]); \ + /* xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_);*/ /* xt1 <- x[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(xt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + xt1 = _mm512_mask_blend_epi8(mskc_, slli512, x1_); \ + x1_ = tmp; \ + vt1 = _mm512_load_si512(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ + /* tmp = _mm_srli_si128(vt1, 15); */ /* tmp <- v[r-1][t+15] */ \ + tmp = _mm512_set1_epi8(((int8_t*)v)[t*64 + 63]); \ + /* vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); *//* vt1 <- v[r-1][t-1..t+14] */ \ + shf512a = _mm512_shuffle_epi8(vt1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + vt1 = _mm512_mask_blend_epi8(mskc_, slli512, v1_); \ + v1_ = tmp; \ + a = _mm512_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ + ut = _mm512_load_si512(&u[t]); /* ut <- u[t..t+15] */ \ + b = _mm512_add_epi8(_mm512_load_si512(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ \ + x2t1= _mm512_load_si512(&x2[t]); \ + /* tmp = _mm_srli_si128(x2t1, 15);*/ \ + tmp = _mm512_set1_epi8(((int8_t*)x2)[t*64 + 63]); \ + /* x2t1= _mm_or_si128(_mm_slli_si128(x2t1, 1), x21_); */ \ + shf512a = _mm512_shuffle_epi8(x2t1, ind512_slli); \ + shf512b = _mm512_shuffle_i32x4(shf512a, shf512a, mska); \ + slli512 = _mm512_mask_blend_epi8(mskb, shf512a, shf512b); \ + x2t1 = _mm512_mask_blend_epi8(mskc_, slli512, x21_); \ + x21_= tmp; \ + a2= _mm512_add_epi8(x2t1, vt1); \ + b2= _mm512_add_epi8(_mm512_load_si512(&y2[t]), ut); + + + #define __dp_code_block2_pcl \ + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ + tmp = _mm512_sub_epi8(z, q_); \ + a = _mm512_sub_epi8(a, tmp); \ + b = _mm512_sub_epi8(b, tmp); \ + tmp = _mm512_sub_epi8(z, q2_); \ + a2= _mm512_sub_epi8(a2, tmp); \ + b2= _mm512_sub_epi8(b2, tmp); + + __mmask64 msk_ar[5] = {0xFFFF, 0xFFFFFFFF, 0xFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF}; + __mmask64 msk_ar2[5] = {0xFFFF, 0xFFFF, 0xFFFFFFFF, 0xFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF}; + + int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff; + int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); + int32_t *H = 0, H0 = 0, last_H0_t = 0; + uint8_t *qr, *sf, *mem, *mem2 = 0; + + __m512i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; + __m512i *u, *v, *x, *y, *x2, *y2, *s, *p = 0; + __m512i one_, two_, three_, four_, s1_, s2_, s3_, s4_; + + ksw_reset_extz(ez); + if (m <= 1 || qlen <= 0 || tlen <= 0) return; + + if (q2 + e2 < q + e) t = q, q = q2, q2 = t, t = e, e = e2, e2 = t; // make sure q+e no larger than q2+e2 + s1_ = _mm512_set1_epi8(0x08); + s2_ = _mm512_set1_epi8(0x10); + s3_ = _mm512_set1_epi8(0x20); + s4_ = _mm512_set1_epi8(0x40); + + one_ = _mm512_set1_epi8(1); + two_ = _mm512_set1_epi8(2); + three_ = _mm512_set1_epi8(3); + four_ = _mm512_set1_epi8(4); + + zero_ = _mm512_set1_epi8(0); + q_ = _mm512_set1_epi8(q); + q2_ = _mm512_set1_epi8(q2); + qe_ = _mm512_set1_epi8(q + e); + qe2_ = _mm512_set1_epi8(q2 + e2); + sc_mch_ = _mm512_set1_epi8(mat[0]); + sc_mis_ = _mm512_set1_epi8(mat[1]); + sc_N_ = mat[m*m-1] == 0? _mm512_set1_epi8(-e2) : _mm512_set1_epi8(mat[m*m-1]); + m1_ = _mm512_set1_epi8(m - 1); // wildcard + + if (w < 0) w = tlen > qlen? tlen : qlen; + wl = wr = w; + tlen_ = (tlen + 63) / 64; + n_col_ = qlen < tlen? qlen : tlen; + n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 63) / 64 + 1; + qlen_ = (qlen + 63) / 64; + for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { + max_sc = max_sc > mat[t]? max_sc : mat[t]; + min_sc = min_sc < mat[t]? min_sc : mat[t]; + } + if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches + + long_thres = e != e2? (q2 - q) / (e - e2) - 1 : 0; + if (q2 + e2 + long_thres * e2 > q + e + long_thres * e) + ++long_thres; + long_diff = long_thres * (e - e2) - (q2 - q) - e2; + + mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1 + 63, 64); + u = (__m512i*)(((size_t)mem + 63) >> 6 << 6); // 16-byte aligned + v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_; + s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 64; + memset(u, -q - e, tlen_ * 64); + memset(v, -q - e, tlen_ * 64); + memset(x, -q - e, tlen_ * 64); + memset(y, -q - e, tlen_ * 64); + memset(x2, -q2 - e2, tlen_ * 64); + memset(y2, -q2 - e2, tlen_ * 64); + if (!approx_max) { + H = (int32_t*)kmalloc(km, tlen_ * 64 * 4); + for (t = 0; t < tlen_ * 64; ++t) H[t] = KSW_NEG_INF; + } + if (with_cigar) { + mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 64); + p = (__m512i*)(((size_t)mem2 + 63) >> 6 << 6); + off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); + off_end = off + qlen + tlen - 1; + } + + for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t]; + memcpy(sf, target, tlen); + + + for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { + int st = 0, en = tlen - 1, st0, en0, st_, en_; + int8_t x1, x21, v1; + uint8_t *qrr = qr + (qlen - 1 - r); + int8_t *u8 = (int8_t*)u, *v8 = (int8_t*)v, *x8 = (int8_t*)x, *x28 = (int8_t*)x2; + __m512i x1_, x21_, v1_; + // find the boundaries + if (st < r - qlen + 1) st = r - qlen + 1; + if (en > r) en = r; + if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil + if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor + if (st > en) { + ez->zdropped = 1; + break; + } + st0 = st, en0 = en; + int st_new = st / 16 * 16, en_new = (en + 16) / 16 * 16 - 1; + // int st_new = st / 64 * 64, en_new = (en + 64) / 64 * 64 - 1; + int stb = st, enb = en; + st = st / 64 * 64, en = (en + 64) / 64 * 64 - 1; + int stn = stb / 16 * 16, enn = (enb + 16) / 16 * 16 - 1; + // set boundary conditions + if (st_new > 0) { + if (st_new - 1 >= last_st && st_new - 1 <= last_en) { + x1 = x8[st_new - 1], x21 = x28[st_new - 1], v1 = v8[st_new - 1]; // (r-1,s-1) calculated in the last round + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = -q - e; + } + } else { + x1 = -q - e, x21 = -q2 - e2; + v1 = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + if (en_new >= r) { + ((int8_t*)y)[r] = -q - e, ((int8_t*)y2)[r] = -q2 - e2; + u8[r] = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; + } + + // loop fission: set scores first + if (!(flag & KSW_EZ_GENERIC_SC)) { + for (t = st0; t <= en0; t += 64) { + __m512i sq, st, tmp_512, mask_512; + __mmask64 tmp, mask; + sq = _mm512_loadu_si512((__m512i*)&sf[t]); + st = _mm512_loadu_si512((__m512i*)&qrr[t]); +// mask = _mm512_or_si512(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_)); + mask = (_mm512_cmpeq_epi8_mask(sq, m1_) | _mm512_cmpeq_epi8_mask(st, m1_)); + tmp = _mm512_cmpeq_epi8_mask(sq, st); + + tmp_512 = _mm512_mask_blend_epi8(tmp, sc_mis_, sc_mch_); + tmp_512 = _mm512_mask_blend_epi8(mask, tmp_512, sc_N_); + if (t + 64 > en0) + { + __mmask64 msk; + int ind = (en0 - t + 16)/16; + //assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + _mm512_mask_storeu_epi8((__m512i*)((int8_t*)s + t), msk, tmp_512); + } + else + _mm512_storeu_si512((__m512i*)((int8_t*)s + t), tmp_512); + + } + } else { + for (t = st0; t <= en0; ++t) + ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]]; + } + + // core loop + // fprintf(stderr, "- r: %d, x1: %d, x21: %d, v1: %d, en_new: %d, e: %d, q: %d\n", + //r, x1, x21, v1, en_new, e, q); + x1_ = _mm512_set1_epi8((uint8_t)x1); + x21_ = _mm512_set1_epi8((uint8_t)x21); + v1_ = _mm512_set1_epi8((uint8_t)v1); + + //st_ = st / 16, en_ = en / 16; + st_ = st / 64, en_ = en / 64; + //assert(en_ - st_ + 1 <= n_col_); + if (!with_cigar) { // score only + for (t = st_; t <= en_; ++t) { + __m512i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __dp_code_block1_pcl; + + z = _mm512_max_epi8(z, a); + z = _mm512_max_epi8(z, b); + z = _mm512_max_epi8(z, a2); + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; // save u[] and v[]; update a, b, a2 and b2 + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + + } + else { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + if (t == en_) { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + //assert(ind >= 0); + msk = msk_ar2[ind]; + // fprintf(stderr, "en0: %d, t: %d, ind: %d, msk: %d\n", en0, t, ind, msk); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + else + { + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_max_epi8(a, zero_), qe_)); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_max_epi8(b, zero_), qe_)); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_max_epi8(a2, zero_), qe2_)); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_max_epi8(b2, zero_), qe2_)); + } + // for (int l=0; l<64; l++) + // fprintf(stderr, "%d ", ((int8_t*)x)[l]); + + } + + } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + + + for (t = st_; t < en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 4 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + } + { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a, z),zero_, one_); // d = a > z? 1 : 0 + z = _mm512_max_epi8(z, a); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b, z), d, two_); // d = b > z? 2 : d + z = _mm512_max_epi8(z, b); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(a2, z), d, three_); // d = a2 > z? 3 : d + z = _mm512_max_epi8(z, a2); + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(b2, z), d, four_); // d = b2 > z? 3 : d + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + // __dp_code_block2_pcl; + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + // //assert(ind >= 0 && ind < 5); + msk = msk_ar2[ind]; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + tmp_mask = _mm512_cmpgt_epi8_mask(a, zero_); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s1_)); // d = a > 0? 1<<3 : 0 + + tmp_mask = _mm512_cmpgt_epi8_mask(b, zero_); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s2_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(a2, zero_); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, a2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s3_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(b2, zero_); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, zero_, b2), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, zero_, s4_)); // d = b > 0? 1<<6 : 0 + //_mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + } + } + + + } else { // gap right-alignment + __m512i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + // off[r] = stn, off_end[r] = enn; + // fprintf(stderr, "t: %d, st0: %d\n", st_, st0); + for (t = st_; t < en_; ++t) { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a), one_, zero_); + z = _mm512_max_epi8(z, a); + // d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), _mm512_set1_epi8(2), d); + // d = z > b? d : 2 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), two_, d); // d = z > b? d : 2 + z = _mm512_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), three_, d); // d = z > a2? d : 3 + z = _mm512_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), four_, d); // d = z > b2? d : 4 + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + __mmask64 msk; + { + _mm512_store_si512(&u[t], _mm512_sub_epi8(z, vt1)); + _mm512_store_si512(&v[t], _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_store_si512(&x[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + // d = a > 0? 1<<3 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s1_, zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_store_si512(&y[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + // d = b > 0? 1<<4 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s2_, zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_store_si512(&x2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + // d = a > 0? 1<<5 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s3_, zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_store_si512(&y2[t], _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + // d = b > 0? 1<<6 : 0 + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s4_, zero_)); // d = b > 0? 1<<6 : 0 + _mm512_store_si512(&pr[t], d); + + } + } + //for (t = st_; t <= en_; ++t)// Last iteration unrolled + { + __m512i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __mmask64 tmp_mask; + + __dp_code_block1_pcl; + + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a), one_, zero_); + z = _mm512_max_epi8(z, a); + // d = z > b? d : 2 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b), two_, d); // d = z > b? d : 2 + z = _mm512_max_epi8(z, b); + // d = z > a2? d : 3 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, a2), three_, d); // d = z > a2? d : 3 + z = _mm512_max_epi8(z, a2); + // d = z > b2? d : 4 + d = _mm512_mask_blend_epi8(_mm512_cmpgt_epi8_mask(z, b2), four_, d); // d = z > b2? d : 4 + z = _mm512_max_epi8(z, b2); + z = _mm512_min_epi8(z, sc_mch_); + + // __dp_code_block2_pcl; + + __mmask64 msk; + { + // __mmask64 msk; + int ind = (en0 - t*64 + 16)/16; + msk = msk_ar2[ind]; + off_end[r] -= (4-ind)*16; + + _mm512_mask_storeu_epi8(&u[t], msk, _mm512_sub_epi8(z, vt1)); + _mm512_mask_storeu_epi8(&v[t], msk, _mm512_sub_epi8(z, ut)); + tmp = _mm512_sub_epi8(z, q_); + a = _mm512_sub_epi8(a, tmp); + b = _mm512_sub_epi8(b, tmp); + tmp = _mm512_sub_epi8(z, q2_); + a2= _mm512_sub_epi8(a2, tmp); + b2= _mm512_sub_epi8(b2, tmp); + } + + { + + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a); + _mm512_mask_storeu_epi8(&x[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s1_, zero_)); // d = a > 0? 1<<3 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b); + _mm512_mask_storeu_epi8(&y[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b, zero_), qe_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s2_, zero_)); // d = b > 0? 1<<4 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, a2); + _mm512_mask_storeu_epi8(&x2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, a2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s3_, zero_)); // d = a > 0? 1<<5 : 0 + tmp_mask = _mm512_cmpgt_epi8_mask(zero_, b2); + _mm512_mask_storeu_epi8(&y2[t], msk, _mm512_sub_epi8(_mm512_mask_blend_epi8(tmp_mask, b2, zero_), qe2_)); + d = _mm512_or_si512(d, _mm512_mask_blend_epi8(tmp_mask, s4_, zero_)); // d = b > 0? 1<<6 : 0 + // _mm512_store_si512(&pr[t], d); + _mm512_mask_storeu_epi8(&pr[t], msk, d); + + } + } + } + + if (!approx_max) { // find the exact max with a 32-bit score array + int32_t max_H, max_t; + // compute H[], max_H and max_t + if (r > 0) { + int32_t HH[16], tt[16], en1 = st0 + (en0 - st0) / 16 * 16, i; + __m512i max_H_, max_t_; + max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] : H[en0] + v8[en0]; // special casing the last element + + max_t = en0; + max_H_ = _mm512_set1_epi32(max_H); + max_t_ = _mm512_set1_epi32(max_t); + for (t = st0; t < en1; t += /*4*/16) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + __m512i H1, tmp, t_; + __mmask16 tmp_mask; + H1 = _mm512_loadu_si512((__m512i*)&H[t]); + __m128i t__ = _mm_load_si128((__m128i*) &v8[t]); + t_ = _mm512_cvtepi8_epi32(t__); + H1 = _mm512_add_epi32(H1, t_); + _mm512_storeu_si512((__m512i*)&H[t], H1); + // making it 4 lanes to match accuracy + + __m512i shfH, shft, max1, max2; + t_ = _mm512_set1_epi32(t); + t_ = _mm512_add_epi32(t_, bt32_); + shfH = _mm512_shuffle_i32x4(H1, H1, 0x31); + shft = _mm512_shuffle_i32x4(t_, t_, 0x31); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, H1); + max1 = _mm512_mask_blend_epi32(tmp_mask, H1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, t_, shft); + shfH = _mm512_shuffle_i32x4(max1, max1, 0x2); + shft = _mm512_shuffle_i32x4(max2, max2, 0x2); + tmp_mask = _mm512_cmpgt_epi32_mask(shfH, max1); + max1 = _mm512_mask_blend_epi32(tmp_mask, max1, shfH); + max2 = _mm512_mask_blend_epi32(tmp_mask, max2, shft); + tmp_mask = _mm512_cmpgt_epi32_mask(max1, max_H_); + max_H_ = _mm512_mask_blend_epi32(tmp_mask, max_H_, max1); + max_t_ = _mm512_mask_blend_epi32(tmp_mask, max_t_, max2); + + + } + _mm512_storeu_si512((__m512i*)HH, max_H_); + _mm512_storeu_si512((__m512i*)tt, max_t_); + + int rem = (en0 - t) / 4; + for (int l=0; l HH[j]) { + HH[j] = H[t]; + tt[j] = bt; + } + t++; + } + } + + for (i = 0; i < 4; ++i) + if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i; + + for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + H[t] += (int32_t)v8[t]; + if (H[t] > max_H) { + max_H = H[t], max_t = t; + } + } + } else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0 + + // update ez + if (en0 == tlen - 1 && H[en0] > ez->mte) { + ez->mte = H[en0], ez->mte_q = r - en_new; + } + if (r - st0 == qlen - 1 && H[st0] > ez->mqe) { + ez->mqe = H[st0], ez->mqe_t = st0; + } + + if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H[tlen - 1]; + } + } else { // find approximate max; Z-drop might be inaccurate, too. + if (r > 0) { + if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { + int32_t d0 = v8[last_H0_t]; + int32_t d1 = u8[last_H0_t + 1]; + if (d0 > d1) H0 += d0; + else H0 += d1, ++last_H0_t; + } else if (last_H0_t >= st0 && last_H0_t <= en0) { + H0 += v8[last_H0_t]; + } else { + ++last_H0_t, H0 += u8[last_H0_t]; + } + } else H0 = v8[0] - qe, last_H0_t = 0; + if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e2)) { + break; + } + if (r == qlen + tlen - 2 && en0 == tlen - 1) { + ez->score = H0; + } + } + // last_st = st, last_en = en; + last_st = st_new, last_en = en_new; + + } + + kfree(km, mem); + if (!approx_max) kfree(km, H); + if (with_cigar) { // backtrack + int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); + if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) { + ez->reach_end = 1; + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (ez->max_t >= 0 && ez->max_q >= 0) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*64, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } + kfree(km, mem2); kfree(km, off); + } +} + +#endif diff --git a/ksw2_extd2_avx.h b/ksw2_extd2_avx.h new file mode 100644 index 0000000..abcea0a --- /dev/null +++ b/ksw2_extd2_avx.h @@ -0,0 +1,40 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ +#include +#include +#include +#include "ksw2.h" +#include +#include +#include +#include + +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); diff --git a/main.c b/main.c index 395c0d9..26d615f 100644 --- a/main.c +++ b/main.c @@ -1,13 +1,54 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ #include #include #include +#include #include #include "bseq.h" #include "minimap.h" #include "mmpriv.h" #include "ketopt.h" +#include #define MM_VERSION "2.18-r1015" +using namespace std; + +#ifdef MANUAL_PROFILING +uint64_t num_reads = 0, minimizer_hit_time = 0, dp_chaining_time = 0, alignment_time = 0; +#endif + +#ifdef LISA_HASH +#include "lisa_hash.h" +lisa_hash *lh; +#endif #ifdef __linux__ #include @@ -108,11 +149,12 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha int main(int argc, char *argv[]) { - 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:"; + 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:Z:"; ketopt_t o = KETOPT_INIT; mm_mapopt_t opt; mm_idxopt_t ipt; int i, c, n_threads = 3, n_parts, old_best_n = -1; + uint64_t total_time = 0; char *fnw = 0, *rg = 0, *junc_bed = 0, *s, *alt_list = 0; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; @@ -122,9 +164,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; @@ -141,6 +185,7 @@ int main(int argc, char *argv[]) while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { if (c == 'w') ipt.w = atoi(o.arg); + else if (c == 'Z') opt.L_hash = atoi(o.arg); else if (c == 'k') ipt.k = atoi(o.arg); else if (c == 'H') ipt.flag |= MM_I_HPC; else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I @@ -345,7 +390,12 @@ 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); + total_time = __rdtsc(); if (idx_rdr == 0) { fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno)); return 1; @@ -387,9 +437,23 @@ 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); + if(opt.L_hash == 1) { + fprintf(stderr, "Generating lisa-hash..\n"); + mm_idx_dump_hash(preset_arg.c_str(), mi); + fprintf(stderr, "Lisa-hash saving done.. \n"); + exit(0); + } if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1); if (alt_list) mm_idx_alt_read(mi, alt_list); ret = 0; +#ifdef LISA_HASH + fprintf(stderr, "Using LISA_HASH..\n"); + mm_idx_destroy_mm_hash(mi); + char* prefix; + lh = new lisa_hash(preset_arg, prefix); + fprintf(stderr, "Loading done.\n"); + total_time = __rdtsc(); +#endif 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); @@ -398,11 +462,15 @@ 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); 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); @@ -422,5 +490,13 @@ 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); } - return 0; + +#ifdef MANUAL_PROFILING + fprintf(stderr, "\n Number of reads = %lld Minimizer hit time = %lld dp_chaining time = %lld alignment time = %lld total time = %lld \n", num_reads, minimizer_hit_time, dp_chaining_time, alignment_time, __rdtsc() - total_time); +#endif + fprintf(stderr, "Total ticks: %lld \n",__rdtsc() - total_time); +#ifdef LISA_HASH + delete lh; +#endif +return 0; } diff --git a/map.c b/map.c index 0812783..9f3a100 100644 --- a/map.c +++ b/map.c @@ -1,3 +1,32 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ #include #include #include @@ -9,6 +38,17 @@ #include "mmpriv.h" #include "bseq.h" #include "khash.h" +#include + +#ifdef LISA_HASH +#include "lisa_hash.h" +extern lisa_hash *lh; +#endif + +#ifdef MANUAL_PROFILING +extern uint64_t num_reads, minimizer_hit_time, dp_chaining_time, alignment_time; +#endif + struct mm_tbuf_s { void *km; @@ -87,6 +127,77 @@ typedef struct { const uint64_t *cr; } mm_match_t; + +#ifdef LISA_HASH +static mm_match_t *collect_matches_lisa_hash(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos) +{ + 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)); + + int rep_st = 0, rep_en = 0, n_m; + size_t i; + mm_match_t *m; + *n_mini_pos = 0; + *mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t)); + m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t)); + + for (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); + for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) { + const uint64_t *cr; + mm128_t *p = &mv->a[i]; + uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff; + int t; + cr = cr_batch[i]; t = t_batch[i]; +/*Correctness check for lisa_hash*/ +#ifdef LISA_HASH_ASSERT + int t_minimap2_original; + const uint64_t *cr_minimap2_hash = mm_idx_get(mi, p->x>>8, &t); + cr_minimap2_hash = mm_idx_get(mi, p->x>>8, &t_minimap2_original); + assert(t == t_minimap2_original); +#endif + if (t >= max_occ) { + + int en = (q_pos >> 1) + 1, st = en - q_span; + if (st > rep_en) { + *rep_len += rep_en - rep_st; + rep_st = st, rep_en = en; + } else rep_en = en; + } else { +#ifdef LISA_HASH_ASSERT + //Correctness assertion + for(int itr = 0; itr < t; itr++){ + assert((cr[itr] == cr_minimap2_hash[itr])); + } +#endif + mm_match_t *q = &m[n_m++]; + q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32; + q->is_tandem = 0; + 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; + *n_a += q->n; + (*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q_pos>>1; + } + } + + free(cr_batch); + free(t_batch); + free(minimizers); + free(lisa_pos); + + *rep_len += rep_en - rep_st; + *_n_m = n_m; + return m; +} +#endif + + static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos) { int rep_st = 0, rep_en = 0, n_m; @@ -146,6 +257,51 @@ static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const cha return 0; } + +static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len, + int *n_mini_pos, uint64_t **mini_pos) +{ + int i, n_m; + mm_match_t *m; + mm128_t *a; +#ifndef LISA_HASH + m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos); +#else + m = collect_matches_lisa_hash(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos); +#endif + uint64_t mapping_count = 0; + mapping_count = 0; + + a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t)); + for (i = 0, *n_a = 0; i < n_m; ++i) { + + mm_match_t *q = &m[i]; + const uint64_t *r = q->cr; + uint32_t k; + for (k = 0; k < q->n; ++k) { + uint64_t r_k = r[k]; + int32_t is_self, rpos = (uint32_t)r_k >> 1; + mm128_t *p; + if (skip_seed(opt->flag, r_k, q, qname, qlen, mi, &is_self)) continue; + p = &a[(*n_a)++]; + if ((r_k&1) == (q->q_pos&1)) { // forward strand + p->x = (r_k & 0xffffffff00000000ULL) | rpos; + p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1; + } else { // reverse strand + p->x = 1ULL<<63 | (r_k & 0xffffffff00000000ULL) | rpos; + p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1); + } + p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT; + if (q->is_tandem) p->y |= MM_SEED_TANDEM; + if (is_self) p->y |= MM_SEED_SELF; + } + } + + kfree(km, m); + radix_sort_128x(a, a + (*n_a)); + return a; +} + static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos) { @@ -212,40 +368,6 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max return a; } -static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len, - int *n_mini_pos, uint64_t **mini_pos) -{ - int i, n_m; - mm_match_t *m; - mm128_t *a; - m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos); - a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t)); - for (i = 0, *n_a = 0; i < n_m; ++i) { - mm_match_t *q = &m[i]; - const uint64_t *r = q->cr; - uint32_t k; - for (k = 0; k < q->n; ++k) { - int32_t is_self, rpos = (uint32_t)r[k] >> 1; - mm128_t *p; - if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue; - p = &a[(*n_a)++]; - if ((r[k]&1) == (q->q_pos&1)) { // forward strand - p->x = (r[k]&0xffffffff00000000ULL) | rpos; - p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1; - } else { // reverse strand - p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos; - p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1); - } - p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT; - if (q->is_tandem) p->y |= MM_SEED_TANDEM; - if (is_self) p->y |= MM_SEED_SELF; - } - } - kfree(km, m); - radix_sort_128x(a, a + (*n_a)); - return a; -} - static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a) { if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s) @@ -271,6 +393,11 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) { + +#ifdef MANUAL_PROFILING + num_reads++; +#endif + int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos; int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR); uint32_t hash; @@ -293,8 +420,18 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv); if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); - else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); + else { +#ifdef MANUAL_PROFILING + uint64_t mm_hit_start = __rdtsc(); +#endif + + a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); + +#ifdef MANUAL_PROFILING + minimizer_hit_time += (__rdtsc() - mm_hit_start); +#endif + } if (mm_dbg_flag & MM_DBG_PRINT_SEED) { fprintf(stderr, "RS\t%d\n", rep_len); for (i = 0; i < n_a; ++i) @@ -312,9 +449,16 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** max_chain_gap_ref = opt->max_frag_len - qlen_sum; if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap; } else max_chain_gap_ref = opt->max_gap; - +#ifdef MANUAL_PROFILING + uint64_t dp_start = __rdtsc(); +#endif a = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); + +#ifdef MANUAL_PROFILING + dp_chaining_time += (__rdtsc() - dp_start); +#endif + if (opt->max_occ > opt->mid_occ && rep_len > 0) { int rechain = 0; if (n_regs0 > 0) { // test if the best chain has all the segments @@ -335,6 +479,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** kfree(b->km, mini_pos); if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); + a = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); } } @@ -430,6 +575,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori; const char *qseqs[MM_MAX_SEG]; mm_tbuf_t *b = s->buf[tid]; + assert(s->n_seg[i] <= MM_MAX_SEG); if (mm_dbg_flag & MM_DBG_PRINT_QNAME) fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq); @@ -561,6 +707,7 @@ static void *worker_pipeline(void *shared, int step, void *in) else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag); return in; } else if (step == 2) { // step 2: output + void *km = 0; step_t *s = (step_t*)in; const mm_idx_t *mi = p->mi; @@ -569,6 +716,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 @@ -603,6 +751,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]); diff --git a/minimap.h b/minimap.h index a00869e..2b8c0e8 100644 --- a/minimap.h +++ b/minimap.h @@ -1,3 +1,32 @@ +/* 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 ; + Vasimuddin Md ; Sanchit Misra ; + Chirag Jain ; Heng Li +*/ #ifndef MINIMAP2_H #define MINIMAP2_H @@ -153,6 +182,8 @@ typedef struct { int64_t max_sw_mat; const char *split_prefix; + // Store minimizer hash to a file as key and list of values + int L_hash; } mm_mapopt_t; // index reader @@ -267,6 +298,14 @@ 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 * @@ -296,6 +335,21 @@ void mm_idx_stat(const mm_idx_t *idx); */ 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 *