1 /**** 2 DIAMOND protein aligner 3 Copyright (C) 2016-2021 Max Planck Society for the Advancement of Science e.V. 4 Benjamin Buchfink 5 6 Code developed by Benjamin Buchfink <benjamin.buchfink@tue.mpg.de> 7 8 This program is free software: you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation, either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see <http://www.gnu.org/licenses/>. 20 ****/ 21 22 namespace DP { namespace Swipe { 23 namespace DISPATCH_ARCH { 24 25 template<typename _sv> 26 struct Matrix 27 { 28 struct ColumnIterator 29 { ColumnIteratorMatrix::ColumnIterator30 ColumnIterator(_sv* hgap_front, _sv* score_front) : 31 hgap_ptr_(hgap_front), 32 score_ptr_(score_front) 33 { } 34 inline void operator++() 35 { 36 ++hgap_ptr_; ++score_ptr_; 37 } hgapMatrix::ColumnIterator38 inline _sv hgap() const 39 { 40 return *hgap_ptr_; 41 } diagMatrix::ColumnIterator42 inline _sv diag() const 43 { 44 return *score_ptr_; 45 } set_hgapMatrix::ColumnIterator46 inline void set_hgap(const _sv& x) 47 { 48 *hgap_ptr_ = x; 49 } set_scoreMatrix::ColumnIterator50 inline void set_score(const _sv& x) 51 { 52 *score_ptr_ = x; 53 } trace_maskMatrix::ColumnIterator54 std::nullptr_t trace_mask() { 55 return nullptr; 56 } 57 _sv *hgap_ptr_, *score_ptr_; 58 }; MatrixMatrix59 Matrix(int rows, int) 60 { 61 hgap_.resize(rows); 62 score_.resize(rows + 1); 63 const auto z = _sv(); 64 std::fill(hgap_.begin(), hgap_.end(), z); 65 std::fill(score_.begin(), score_.end(), z); 66 } beginMatrix67 inline ColumnIterator begin(int) 68 { 69 return ColumnIterator(hgap_.begin(), score_.begin()); 70 } set_zeroMatrix71 void set_zero(int c) 72 { 73 const int l = (int)hgap_.size(); 74 const auto z = extract_channel(_sv(), 0); 75 for (int i = 0; i < l; ++i) { 76 hgap_[i] = set_channel(hgap_[i], c, z); 77 score_[i] = set_channel(score_[i], c, z); 78 } 79 score_[l] = set_channel(score_[l], c, z); 80 } colsMatrix81 constexpr int cols() const { 82 return 1; 83 } 84 _sv operator[](int i) const { 85 return score_[i + 1]; 86 } 87 private: 88 #ifdef __APPLE__ 89 MemBuffer<_sv> hgap_, score_; 90 #else 91 static thread_local MemBuffer<_sv> hgap_, score_; 92 #endif 93 }; 94 95 #ifndef __APPLE__ 96 template<typename _sv> thread_local MemBuffer<_sv> Matrix<_sv>::hgap_; 97 template<typename _sv> thread_local MemBuffer<_sv> Matrix<_sv>::score_; 98 #endif 99 100 template<typename _sv> 101 struct TracebackVectorMatrix 102 { 103 typedef typename ::DISPATCH_ARCH::ScoreTraits<_sv>::TraceMask TraceMask; 104 typedef void* Stat; 105 struct ColumnIterator 106 { ColumnIteratorTracebackVectorMatrix::ColumnIterator107 ColumnIterator(_sv* hgap_front, _sv* score_front, TraceMask* trace_mask_front) : 108 hgap_ptr_(hgap_front), 109 score_ptr_(score_front), 110 trace_mask_ptr_(trace_mask_front) 111 { } 112 inline void operator++() 113 { 114 ++hgap_ptr_; ++score_ptr_; ++trace_mask_ptr_; 115 } hgapTracebackVectorMatrix::ColumnIterator116 inline _sv hgap() const 117 { 118 return *hgap_ptr_; 119 } diagTracebackVectorMatrix::ColumnIterator120 inline _sv diag() const 121 { 122 return *score_ptr_; 123 } trace_maskTracebackVectorMatrix::ColumnIterator124 inline TraceMask* trace_mask() { 125 return trace_mask_ptr_; 126 } set_hgapTracebackVectorMatrix::ColumnIterator127 inline void set_hgap(const _sv& x) 128 { 129 *hgap_ptr_ = x; 130 } set_scoreTracebackVectorMatrix::ColumnIterator131 inline void set_score(const _sv& x) 132 { 133 *score_ptr_ = x; 134 } statTracebackVectorMatrix::ColumnIterator135 std::nullptr_t stat() { 136 return nullptr; 137 } hstatTracebackVectorMatrix::ColumnIterator138 std::nullptr_t hstat() { 139 return nullptr; 140 } set_hstatTracebackVectorMatrix::ColumnIterator141 void set_hstat(std::nullptr_t) {} set_zeroTracebackVectorMatrix::ColumnIterator142 inline void set_zero() {} 143 _sv *hgap_ptr_, *score_ptr_; 144 TraceMask* trace_mask_ptr_; 145 }; 146 147 struct TracebackIterator 148 { TracebackIteratorTracebackVectorMatrix::TracebackIterator149 TracebackIterator(const TraceMask *mask, const TraceMask* mask_begin, const TraceMask* mask_end, int rows, int i, int j, size_t channel) : 150 rows_(rows), 151 mask_(mask), 152 mask_begin_(mask_begin), 153 mask_end_(mask_end), 154 channel_mask_vgap(TraceMask::vmask(channel)), 155 channel_mask_hgap(TraceMask::hmask(channel)), 156 i(i), 157 j(j) 158 { 159 assert(i >= 0 && j >= 0); 160 } wrap_maskTracebackVectorMatrix::TracebackIterator161 void wrap_mask() { 162 if (mask_ < mask_begin_) 163 mask_ = mask_end_ - (mask_begin_ - mask_); 164 } maskTracebackVectorMatrix::TracebackIterator165 TraceMask mask() const { 166 return *mask_; 167 } walk_diagonalTracebackVectorMatrix::TracebackIterator168 void walk_diagonal() 169 { 170 mask_ -= rows_ + 1; 171 wrap_mask(); 172 --i; 173 --j; 174 assert(i >= -1 && j >= -1); 175 } walk_gapTracebackVectorMatrix::TracebackIterator176 pair<Edit_operation, int> walk_gap() 177 { 178 if (mask_->gap & channel_mask_vgap) { 179 int l = 0; 180 do { 181 ++l; 182 --i; 183 --mask_; 184 } while (((mask_->open & channel_mask_vgap) == 0) && (i > 0)); 185 return std::make_pair(op_insertion, l); 186 } 187 else { 188 int l = 0; 189 do { 190 ++l; 191 --j; 192 mask_ -= rows_; 193 wrap_mask(); 194 } while (((mask_->open & channel_mask_hgap) == 0) && (j > 0)); 195 return std::make_pair(op_deletion, l); 196 } 197 } 198 const int rows_; 199 const TraceMask* mask_, *mask_begin_, *mask_end_; 200 const decltype(TraceMask::gap) channel_mask_vgap, channel_mask_hgap; 201 int i, j; 202 }; 203 tracebackTracebackVectorMatrix204 TracebackIterator traceback(int col, int i, int j, size_t channel) const 205 { 206 return TracebackIterator(&trace_mask_[col*rows_ + i], trace_mask_.begin(), trace_mask_.end(), rows_, i, j, channel); 207 } 208 TracebackVectorMatrixTracebackVectorMatrix209 TracebackVectorMatrix(int rows, int cols) : 210 rows_(rows), 211 cols_(cols) 212 { 213 hgap_.resize(rows); 214 score_.resize(rows + 1); 215 trace_mask_.resize(cols * rows); 216 std::fill(hgap_.begin(), hgap_.end(), _sv()); 217 std::fill(score_.begin(), score_.end(), _sv()); 218 } 219 beginTracebackVectorMatrix220 inline ColumnIterator begin(int col) 221 { 222 return ColumnIterator(hgap_.begin(), score_.begin(), &trace_mask_[col*rows_]); 223 } 224 set_zeroTracebackVectorMatrix225 void set_zero(int c) 226 { 227 const int l = (int)hgap_.size(); 228 for (int i = 0; i < l; ++i) { 229 hgap_[i] = set_channel(hgap_[i], c, ::DISPATCH_ARCH::ScoreTraits<_sv>::zero_score()); 230 score_[i] = set_channel(score_[i], c, ::DISPATCH_ARCH::ScoreTraits<_sv>::zero_score()); 231 } 232 score_[l] = set_channel(score_[l], c, ::DISPATCH_ARCH::ScoreTraits<_sv>::zero_score()); 233 } 234 colsTracebackVectorMatrix235 int cols() const { 236 return cols_; 237 } 238 239 _sv operator[](int i) const { 240 return _sv(); 241 } 242 243 #ifdef __APPLE__ 244 MemBuffer<_sv> hgap_, score_; 245 #else 246 static thread_local MemBuffer<_sv> hgap_, score_; 247 #endif 248 MemBuffer<TraceMask> trace_mask_; 249 private: 250 int rows_, cols_; 251 }; 252 253 #ifndef __APPLE__ 254 template<typename _sv> thread_local MemBuffer<_sv> TracebackVectorMatrix<_sv>::hgap_; 255 template<typename _sv> thread_local MemBuffer<_sv> TracebackVectorMatrix<_sv>::score_; 256 #endif 257 258 template<typename Sv, bool Traceback> 259 struct SelectMatrix { 260 }; 261 262 template<typename Sv> 263 struct SelectMatrix<Sv, true> { 264 using Type = TracebackVectorMatrix<Sv>; 265 }; 266 267 template<typename Sv> 268 struct SelectMatrix<Sv, false> { 269 using Type = Matrix<Sv>; 270 }; 271 272 }}}