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 }}}