1 // ==========================================================================
2 // SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2018, Knut Reinert, FU Berlin
5 // All rights reserved.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are met:
9 //
10 // * Redistributions of source code must retain the above copyright
11 // notice, this list of conditions and the following disclaimer.
12 // * Redistributions in binary form must reproduce the above copyright
13 // notice, this list of conditions and the following disclaimer in the
14 // documentation and/or other materials provided with the distribution.
15 // * Neither the name of Knut Reinert or the FU Berlin nor the names of
16 // its contributors may be used to endorse or promote products derived
17 // from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
29 // DAMAGE.
30 //
31 // ==========================================================================
32 // Author: Rene Rahn <rene.rahn@fu-berlin.de>
33 // ==========================================================================
34 // Defines the recursion formula for the dp-alignment algorithms.
35 // ==========================================================================
36
37 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_FORMULA_H_
38 #define SEQAN_INCLUDE_SEQAN_ALIGN_DP_FORMULA_H_
39
40 namespace seqan {
41
42 // ============================================================================
43 // Forwards
44 // ============================================================================
45
46 // ============================================================================
47 // Tags, Classes, Enums
48 // ============================================================================
49
50 // ----------------------------------------------------------------------------
51 // Tag RecursionDirectionDiagonal
52 // ----------------------------------------------------------------------------
53
54 struct RecursionDirectionDiagonal_;
55 typedef Tag<RecursionDirectionDiagonal_> RecursionDirectionDiagonal;
56
57 // ----------------------------------------------------------------------------
58 // Tag RecursionDirectionHorizontal
59 // ----------------------------------------------------------------------------
60
61 struct RecursionDirectionHorizontal_;
62 typedef Tag<RecursionDirectionHorizontal_> RecursionDirectionHorizontal;
63
64 // ----------------------------------------------------------------------------
65 // Tag RecursionDirectionVertical
66 // ----------------------------------------------------------------------------
67
68 struct RecursionDirectionVertical_;
69 typedef Tag<RecursionDirectionVertical_> RecursionDirectionVertical;
70
71 // ----------------------------------------------------------------------------
72 // Tag RecursionDirectionAll
73 // ----------------------------------------------------------------------------
74
75 struct RecursionDirectionAll_;
76 typedef Tag<RecursionDirectionAll_> RecursionDirectionAll;
77
78 // ----------------------------------------------------------------------------
79 // Tag RecursionDirectionUpperDiagonal
80 // ----------------------------------------------------------------------------
81
82 struct RecursionDirectionUpperDiagonal_;
83 typedef Tag<RecursionDirectionUpperDiagonal_> RecursionDirectionUpperDiagonal;
84
85 // ----------------------------------------------------------------------------
86 // Tag RecursionDirectionLowerDiagonal
87 // ----------------------------------------------------------------------------
88
89 struct RecursionDirectionLowerDiagonal_;
90 typedef Tag<RecursionDirectionLowerDiagonal_> RecursionDirectionLowerDiagonal;
91
92 // ----------------------------------------------------------------------------
93 // Tag RecursionDirectionZero
94 // ----------------------------------------------------------------------------
95
96 struct RecursionDirectionZero_;
97 typedef Tag<RecursionDirectionZero_> RecursionDirectionZero;
98
99 // ============================================================================
100 // Metafunctions
101 // ============================================================================
102
103 // Helper typedef to get the correct score value type from the score-matrix navigator.
104 template <typename TCellTuple>
105 using ExtractedScoreValueType_ = std::decay_t<decltype(_scoreOfCell(std::get<0>(std::declval<TCellTuple>())))>;
106
107 // ============================================================================
108 // Functions
109 // ============================================================================
110
111 template <typename TTarget,
112 typename TSourceLeft,
113 typename TSourceRight,
114 typename TTraceLeft,
115 typename TTraceRight>
SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>,typename TraceBitMap_<TTarget>::Type)116 inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>, typename TraceBitMap_<TTarget>::Type)
117 _maxScore(TTarget & target,
118 TSourceLeft const & srcLeft,
119 TSourceRight const & srcRight,
120 TTraceRight const /**/,
121 TTraceLeft const /**/,
122 TracebackOff const & /*tag*/)
123 {
124 using std::max;
125 target = max(static_cast<TTarget>(srcLeft), static_cast<TTarget>(srcRight));
126 return TraceBitMap_<TTarget>::NONE;
127 }
128
129 template <typename TTarget,
130 typename TSourceLeft,
131 typename TSourceRight,
132 typename TTraceLeft,
133 typename TTraceRight>
SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>,typename TraceBitMap_<TTarget>::Type)134 inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>, typename TraceBitMap_<TTarget>::Type)
135 _maxScore(TTarget & target,
136 TSourceLeft const & srcLeft,
137 TSourceRight const & srcRight,
138 TTraceRight const /**/,
139 TTraceLeft const /**/,
140 TracebackOff const & /*tag*/)
141 {
142 target = max(srcLeft, srcRight);
143 return TraceBitMap_<TTarget>::NONE;
144 }
145
146 template <typename TTarget,
147 typename TSourceLeft,
148 typename TSourceRight,
149 typename TTraceLeft,
150 typename TTraceRight,
151 typename TGapsPlacement>
SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>,typename TraceBitMap_<TTarget>::Type)152 inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>, typename TraceBitMap_<TTarget>::Type)
153 _maxScore(TTarget & target,
154 TSourceLeft const & srcLeft,
155 TSourceRight const & srcRight,
156 TTraceLeft const traceLeft,
157 TTraceRight const traceRight,
158 TracebackOn<TracebackConfig_<SingleTrace, TGapsPlacement> > const & /*tag*/)
159 {
160 return (srcLeft < srcRight)
161 ? (target = srcRight, traceRight)
162 : (target = srcLeft, traceLeft);
163 }
164
165 template <typename TTarget,
166 typename TSourceLeft,
167 typename TSourceRight,
168 typename TTraceLeft,
169 typename TTraceRight,
170 typename TGapsPlacement>
SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>,typename TraceBitMap_<TTarget>::Type)171 inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>, typename TraceBitMap_<TTarget>::Type)
172 _maxScore(TTarget & target,
173 TSourceLeft const & srcLeft,
174 TSourceRight const & srcRight,
175 TTraceLeft const traceLeft,
176 TTraceRight const traceRight,
177 TracebackOn<TracebackConfig_<SingleTrace, TGapsPlacement> > const & /*tag*/)
178 {
179 auto cmp = cmpGt(srcRight, srcLeft);
180 target = blend(srcLeft, srcRight, cmp);
181 return blend(traceLeft, traceRight, cmp);
182 }
183
184 template <typename TTarget,
185 typename TSourceLeft,
186 typename TSourceRight,
187 typename TTraceLeft,
188 typename TTraceRight,
189 typename TGapsPlacement>
SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>,typename TraceBitMap_<TTarget>::Type)190 inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TTarget>>>, typename TraceBitMap_<TTarget>::Type)
191 _maxScore(TTarget & target,
192 TSourceLeft const & srcLeft,
193 TSourceRight const & srcRight,
194 TTraceLeft const traceLeft,
195 TTraceRight const traceRight,
196 TracebackOn<TracebackConfig_<CompleteTrace, TGapsPlacement> > const & /*tag*/)
197 {
198 return (srcLeft == srcRight)
199 ? (target = srcLeft, traceLeft | traceRight)
200 : (srcLeft < srcRight)
201 ? (target = srcRight, traceRight)
202 : (target = srcLeft, traceLeft);
203 }
204
205 template <typename TTarget,
206 typename TSourceLeft,
207 typename TSourceRight,
208 typename TTraceLeft,
209 typename TTraceRight,
210 typename TGapsPlacement>
SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>,typename TraceBitMap_<TTarget>::Type)211 inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TTarget>>, typename TraceBitMap_<TTarget>::Type)
212 _maxScore(TTarget & target,
213 TSourceLeft const & srcLeft,
214 TSourceRight const & srcRight,
215 TTraceLeft const traceLeft,
216 TTraceRight const traceRight,
217 TracebackOn<TracebackConfig_<CompleteTrace, TGapsPlacement> > const & /*tag*/)
218 {
219 auto cmpG = cmpGt(srcRight, srcLeft);
220 auto cmpE = cmpEq(srcRight, srcLeft);
221 target = blend(srcLeft, srcRight, cmpG);
222 auto result = blend(traceLeft, traceRight, cmpG);
223 return blend(result, traceLeft | traceRight, cmpE);
224 }
225
226 // ----------------------------------------------------------------------------
227 // Function _computeScore [RecursionDirectionDiagonal]
228 // ----------------------------------------------------------------------------
229 // Independent of gap cost model.
230 template <typename TScoreValue, typename TGapCosts,
231 typename TSequenceHValue,
232 typename TSequenceVValue,
233 typename TScoringScheme,
234 typename TAlgorithm, typename TTracebackConfig, typename TExecPolicy>
235 inline auto
_computeScore(DPCell_<TScoreValue,TGapCosts> & current,DPCell_<TScoreValue,TGapCosts> & diagonal,DPCell_<TScoreValue,TGapCosts> const &,DPCell_<TScoreValue,TGapCosts> const &,TSequenceHValue const & seqHVal,TSequenceVValue const & seqVVal,TScoringScheme const & scoringScheme,RecursionDirectionDiagonal const &,DPProfile_<TAlgorithm,TGapCosts,TTracebackConfig,TExecPolicy> const &)236 _computeScore(DPCell_<TScoreValue, TGapCosts> & current,
237 DPCell_<TScoreValue, TGapCosts> & diagonal,
238 DPCell_<TScoreValue, TGapCosts> const & /*horizontal*/,
239 DPCell_<TScoreValue, TGapCosts> const & /*vertical*/,
240 TSequenceHValue const & seqHVal,
241 TSequenceVValue const & seqVVal,
242 TScoringScheme const & scoringScheme,
243 RecursionDirectionDiagonal const &,
244 DPProfile_<TAlgorithm, TGapCosts, TTracebackConfig, TExecPolicy> const &)
245 {
246 _scoreOfCell(current) = _scoreOfCell(diagonal) + score(scoringScheme, seqHVal, seqVVal);
247 setGapExtension(current, False(), False(), createVector<TScoreValue>(-1));
248
249 if (IsLocalAlignment_<TAlgorithm>::VALUE)
250 {
251 return _maxScore(_scoreOfCell(current),
252 TraceBitMap_<TScoreValue>::NONE,
253 _scoreOfCell(current),
254 TraceBitMap_<TScoreValue>::NONE,
255 TraceBitMap_<TScoreValue>::DIAGONAL,
256 TTracebackConfig{});
257 }
258 else
259 {
260 return TraceBitMap_<TScoreValue>::DIAGONAL;
261 }
262 }
263
264 // ----------------------------------------------------------------------------
265 // Function _computeScore [RecursionDirectionZero]
266 // ----------------------------------------------------------------------------
267
268 // Independent of gap cost model.
269 template <typename TScoreValue, typename TGapCosts,
270 typename TSequenceHValue,
271 typename TSequenceVValue,
272 typename TScoringScheme,
273 typename TAlgoTag, typename TTraceFlag, typename TExecPolicy>
274 inline auto
_computeScore(DPCell_<TScoreValue,TGapCosts> & current,DPCell_<TScoreValue,TGapCosts> & diagonal,DPCell_<TScoreValue,TGapCosts> const &,DPCell_<TScoreValue,TGapCosts> & vertical,TSequenceHValue const &,TSequenceVValue const &,TScoringScheme const &,RecursionDirectionZero const &,DPProfile_<TAlgoTag,TGapCosts,TTraceFlag,TExecPolicy> const &)275 _computeScore(DPCell_<TScoreValue, TGapCosts> & current,
276 DPCell_<TScoreValue, TGapCosts> & diagonal,
277 DPCell_<TScoreValue, TGapCosts> const & /*horizontal*/,
278 DPCell_<TScoreValue, TGapCosts> & vertical,
279 TSequenceHValue const & /*seqHVal*/,
280 TSequenceVValue const & /*seqVVal*/,
281 TScoringScheme const & /*scoringScheme*/,
282 RecursionDirectionZero const &,
283 DPProfile_<TAlgoTag, TGapCosts, TTraceFlag, TExecPolicy> const &)
284 {
285 _scoreOfCell(current) = TraceBitMap_<TScoreValue>::NONE;
286 _scoreOfCell(diagonal) = _scoreOfCell(current);
287 _scoreOfCell(vertical) = _scoreOfCell(current);
288 return TraceBitMap_<TScoreValue>::NONE;
289 }
290
291 } // namespace seqan
292
293 #endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_FORMULA_H_
294