MikkelSchubert / adapterremoval / #36

Committed 22 Jul 2024 09:33AM UTC coverage: 87.26% (-12.7%) from 100.0%

Build # #36

Build Type

push

travis-ci

Committed by

MikkelSchubert

Commit Message

remove duplicate tests

Run Details

2185 of 2504 relevant lines covered (87.26%)

16293.15 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/simd_avx512bw.cpp

/*************************************************************************\
 * AdapterRemoval - cleaning next-generation sequencing reads            *
 *                                                                       *
 * Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com           *
 *                                                                       *
 * This program is free software: you can redistribute it and/or modify  *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * This program is distributed in the hope that it will be useful,       *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have received a copy of the GNU General Public License     *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>. *
\*************************************************************************/
#include "simd.hpp"    // declarations
#include <bitset>      // for bitset
#include <cstddef>     // for size_t
#include <immintrin.h> // for _mm512_set1_epi8, __m512i, _mm512_lo...

namespace adapterremoval {

namespace simd {

namespace {

/** Counts the number of masked bytes **/
inline auto
count_masked_avx512(__mmask64 value)
{
  // Generate 64 bit mask from most significant bits of each byte and count bits
  return std::bitset<64>(value).count();
}

} // namespace

bool
compare_subsequences_avx512(size_t& n_mismatches,
                            size_t& n_ambiguous,
                            const char* seq_1,
                            const char* seq_2,
                            size_t length,
                            size_t max_penalty)
{
  //! Mask of all Ns
  const auto n_mask = _mm512_set1_epi8('N');

  while (length) {
    __m512i s1;
    __m512i s2;

    if (length < 64) {
      const auto mask = 0xFFFFFFFFFFFFFFFFLLU >> (64 - length);
      s1 = _mm512_maskz_loadu_epi8(mask, seq_1);
      s2 = _mm512_maskz_loadu_epi8(mask, seq_2);

      length = 0;
    } else {
      s1 = _mm512_loadu_epi8(seq_1);
      s2 = _mm512_loadu_epi8(seq_2);

      seq_1 += 64;
      seq_2 += 64;
      length -= 64;
    }

    // Sets 0xFF for every byte where one or both nts is N
    const auto ns_mask =
      _mm512_cmpeq_epu8_mask(s1, n_mask) | _mm512_cmpeq_epu8_mask(s2, n_mask);

    // Sets 0xFF for every byte where bytes are equal or N
    const auto eq_mask = _mm512_cmpeq_epu8_mask(s1, s2) | ns_mask;

    n_mismatches += 64 - count_masked_avx512(eq_mask);
    if (2 * n_mismatches + n_ambiguous > max_penalty) {
      return false;
    }

    // Early termination is almost always due to mismatches, so updating the
    // number of Ns after the above check saves time in the common case.
    n_ambiguous += count_masked_avx512(ns_mask);
  }

  return true;
}

} // namespace simd

} // namespace adapterremoval

1	/*************************************************************************\
2	* AdapterRemoval - cleaning next-generation sequencing reads *
3	* *
4	* Copyright (C) 2024 by Mikkel Schubert - mikkelsch@gmail.com *
5	* *
6	* This program is free software: you can redistribute it and/or modify *
7	* it under the terms of the GNU General Public License as published by *
8	* the Free Software Foundation, either version 3 of the License, or *
9	* (at your option) any later version. *
10	* *
11	* This program is distributed in the hope that it will be useful, *
12	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14	* GNU General Public License for more details. *
15	* *
16	* You should have received a copy of the GNU General Public License *
17	* along with this program. If not, see <http://www.gnu.org/licenses/>. *
18	\*************************************************************************/
19	#include "simd.hpp" // declarations
20	#include <bitset> // for bitset
21	#include <cstddef> // for size_t
22	#include <immintrin.h> // for _mm512_set1_epi8, __m512i, _mm512_lo...
23
24	namespace adapterremoval {
25
26	namespace simd {
27
28	namespace {
29
30	/ Counts the number of masked bytes /
31	inline auto
32	count_masked_avx512(__mmask64 value)	×
33	{
34	// Generate 64 bit mask from most significant bits of each byte and count bits
35	return std::bitset<64>(value).count();	×
36	}
37
38	} // namespace
39
40	bool
41	compare_subsequences_avx512(size_t& n_mismatches,	×
42	size_t& n_ambiguous,
43	const char* seq_1,
44	const char* seq_2,
45	size_t length,
46	size_t max_penalty)
47	{
48	//! Mask of all Ns
49	const auto n_mask = _mm512_set1_epi8('N');	×
50
51	while (length) {	×
52	__m512i s1;	×
53	__m512i s2;	×
54
55	if (length < 64) {	×
56	const auto mask = 0xFFFFFFFFFFFFFFFFLLU >> (64 - length);	×
57	s1 = _mm512_maskz_loadu_epi8(mask, seq_1);	×
58	s2 = _mm512_maskz_loadu_epi8(mask, seq_2);	×
59
60	length = 0;	×
61	} else {
62	s1 = _mm512_loadu_epi8(seq_1);	×
63	s2 = _mm512_loadu_epi8(seq_2);	×
64
65	seq_1 += 64;	×
66	seq_2 += 64;	×
67	length -= 64;	×
68	}
69
70	// Sets 0xFF for every byte where one or both nts is N
71	const auto ns_mask =	×
72	_mm512_cmpeq_epu8_mask(s1, n_mask) \| _mm512_cmpeq_epu8_mask(s2, n_mask);	×
73
74	// Sets 0xFF for every byte where bytes are equal or N
75	const auto eq_mask = _mm512_cmpeq_epu8_mask(s1, s2) \| ns_mask;	×
76
77	n_mismatches += 64 - count_masked_avx512(eq_mask);	×
78	if (2 * n_mismatches + n_ambiguous > max_penalty) {	×
79	return false;
80	}
81
82	// Early termination is almost always due to mismatches, so updating the
83	// number of Ns after the above check saves time in the common case.
84	n_ambiguous += count_masked_avx512(ns_mask);	×
85	}
86
87	return true;
88	}
89
90	} // namespace simd
91
92	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #36

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous