// ==========================================================================
// SeqAn - The Library for Sequence Analysis
// ==========================================================================
// Copyright (c) 2006-2018, Knut Reinert, FU Berlin
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the name of Knut Reinert or the FU Berlin nor the names of
// its contributors may be used to endorse or promote products derived
// from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
// DAMAGE.
//
// ==========================================================================
// Author: Rene Rahn <rene.rahn@fu-berlin.de>
// ==========================================================================
// Implements the core of the dp algorithms.
// This is the crucial part of the refactoring of the alignment algorithms.
// It implements - at the moment only a column wise approach - the core
// loop structure for all alignment profiles. We generally differ between an
// unbanded alignment which is very easy, a banded alignment and a special
// case of the banded alignment the Hamming distance, where upper diagonal
// equals lower diagonal.
//
// The unbanded alignment:
// The computation of the unbanded alignment is divided into three parts.
// In the following we refer to a track as the part where the inner loop
// is iterating (in case of column wise navigation a track is equivalent
// to a column).
// First we compute the initial track. Afterwards we continue with all
// inner tracks of the dp matrix and in the end we compute the last track
// separately. This is because all three types have a special property that
// is different from the other track types.
// Each track itself is further divided into three parts, namely the first cell
// the inner cell and the last cell, corresponding to the initial row,
// all inner rows and the last row of a typical dp matrix. This partition of
// the dp matrix allows us to easily change the behavior of different cells
// according to the chosen dp profile at compile time.
// See alignment_dp_meta_info.h to learn about the different meta objects
// that manage the characteristics of each cell of a particular track type.
//
// The banded alignment:
// In the banded alignment we generally divide the dp matrix into the same
// partition as for the unbanded alignment. The only difference is that we,
// additionally add a specific information of how the current track is
// located within the dp matrix. Since we only consider a band we do not
// necessarily span over the full matrix size for a particular column.
// We distinguish between the locations: PartialColumnTop,
// PartialColumnMiddle, PartialColumnBottom and FullColumn (which is the
// default location for unbanded alignments). Each location of the column
// implies a different composition of the cells contained within a
// particular track. Thus, we are able to set different recursion
// directions and tracking informations for each cell independent from the
// runtime. The only difference is that the outer for-loop (iterating over
// the tracks) is split up into three loops. The first loop then only
// iterates over these tracks that are located at the top of the matrix.
// The second for-loop iterates over the tracks that either are of type
// PartialColumnMiddle or FullColumn (wide bands, where the upper diagonal
// begins behind the track where the lower diagonal crosses the last row of
// the dp matrix). And the last for-loop iterates over the tail of the band
// which is located at the PartialColumnBottom.
//
// The Hamming distance:
// In the special case where upper diagonal equals lower diagonal we only
// have to parse one diagonal of the matrix so we have a special
// implementation for that, though it works for all dp profiles.
//
// Restricitons:
// At the moment we have implemented a restriction such that not all bands
// are accepted. If the dp profile consists of the standard global alignment
// algorithm (NeedlemanWunsch or Gotoh), the band is required to go through
// the sink and the source of the dp matrix. If this is not given the
// alignment algorithm is aborted and the score std::numeric_limits<TScoreValue>::min()
// is returned.
// There are no further restrictions.
//
// GapCosts:
// Another detail of the new module is the selection of the gap functions,
// which is also now part of the compile time configuration. Whenever an
// algorithm is implemented it would automatically work for both gap
// functions (linear gaps and affine gaps).
//
// Navigation:
// It is possible to a certain degree to change the behavior of how to parse
// through the dp matrix. Using the new navigators one can implement
// different behaviors for different matrices. At the moment we only support
// column wise navigation for full and sparse score matrices and for full
// traceback matrices. Another detail of this navigators comes into account,
// when we want to compute only the score. We actually create a navigator
// for the dp matrix but implemented it this way that it gets never actually
// called when the traceback is disabled. Thus we do not store the traceback
// matrix if it is not necessary.
//
// Traceback:
// The traceback is now implemented as a single function that is used by all
// alignment profiles. Here we prefer the diagonal direction before the
// vertical before the horizontal direction.
// All tracebacks are first stored within the String<TraceSegment> object
// and afterwards, when the traceback is finished adapted to its given
// target object such as Align, Graph, Fragments, etc.
//
// Tracing:
// We use now an object called DPScout to keep track of the maximal score.
// This object scouts for the best value and can be overloaded to implement
// different strategies of how the score should be traced. Togehter with
// the meta_info file it only traces thus cells that are allowed to be
// traced given the current dp profile. Since this is also a compile time
// property we do not need to track every cell for the global alignment,
// while we do in the local alignment.
//
// Structure:
// The sequences within the matrix are marked as horizontal and vertical
// sequence to determine there orientation within the matrix.
// ==========================================================================
#ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_ALGORITHM_IMPL_H_
#define SEQAN_INCLUDE_SEQAN_ALIGN_DP_ALGORITHM_IMPL_H_
namespace seqan {
// ============================================================================
// Forwards
// ============================================================================
// ============================================================================
// Tags, Classes, Enums
// ============================================================================
// ============================================================================
// Metafunctions
// ============================================================================
// ============================================================================
// Functions
// ============================================================================
// ----------------------------------------------------------------------------
// Function prepareAlign()
// ----------------------------------------------------------------------------
template<typename TSequence, typename TAlignSpec>
inline void
prepareAlign(StringSet<Align<TSequence, TAlignSpec> > & align,
TSequence const & strH,
StringSet<TSequence> const & setV)
{
size_t numAlignments = length(setV);
SEQAN_ASSERT_EQ(length(align), 0u);
SEQAN_ASSERT_GT(numAlignments, 0u);
resize(align, numAlignments);
for(size_t i = 0; i < numAlignments; ++i)
{
resize(rows(align[i]), 2);
assignSource(row(align[i], 0), strH);
assignSource(row(align[i], 1), setV[i]);
}
}
// ----------------------------------------------------------------------------
// Function _checkBandProperties()
// ----------------------------------------------------------------------------
// Checks whether the chosen band fits the dp profile.
template <typename TSequenceH, typename TSequenceV, typename TAlignmentProfile>
inline bool _checkBandProperties(TSequenceH const & /*seqH*/,
TSequenceV const & /*seqV*/,
DPBandConfig<BandOff> const & /*band*/,
TAlignmentProfile const & /*alignProfile*/)
{
return true;
}
template <typename TSequenceH, typename TSequenceV, typename TAlignmentProfile>
inline bool _checkBandProperties(TSequenceH const & seqH,
TSequenceV const & seqV,
DPBandConfig<BandOn> const & band,
TAlignmentProfile const & /*alignProfile*/)
{
typedef typename MakeSigned<typename Size<TSequenceH>::Type>::Type TSignedSize;
// Check if the intersection between band and DP matrix is empty.
if (upperDiagonal(band) < (0 - static_cast<TSignedSize>(length(seqV))) ||
lowerDiagonal(band) > static_cast<TSignedSize>(length(seqH)))
{
return false;
}
// If the band begins before the beginning of the horizontal sequence
// then check if free end-gaps are enabled at the beginning of the vertical sequence.
if (upperDiagonal(band) < 0 && !IsFreeEndGap_<TAlignmentProfile, DPFirstColumn>::VALUE)
return false;
// If the band begins before the beginning of the vertical sequence
// then check if free end-gaps are enabled at the beginning of the horizontal sequence.
if (lowerDiagonal(band) > 0 && !IsFreeEndGap_<TAlignmentProfile, DPFirstRow>::VALUE)
return false;
// If the band ends behind the end of the vertical sequence
// then check if free end-gaps are enabled at the end of the horizontal sequence.
if (upperDiagonal(band) + static_cast<TSignedSize>(length(seqV)) < static_cast<TSignedSize>(length(seqH)) &&
!IsFreeEndGap_<TAlignmentProfile, DPLastRow>::VALUE)
{
return false;
}
// If the band ends behind the end of the horizontal sequence
// then check if free end-gaps are enabled at the end of the vertical sequence.
if (lowerDiagonal(band) + static_cast<TSignedSize>(length(seqV)) > static_cast<TSignedSize>(length(seqH)) &&
!IsFreeEndGap_<TAlignmentProfile, DPLastColumn>::VALUE)
{
return false;
}
return true;
}
// ----------------------------------------------------------------------------
// Function _invalidDPSettings()
// ----------------------------------------------------------------------------
// Checks if the settings for the dp algorithm are valid.
// Returns true if they are valid, false otherwise.
template <typename TSequenceH, typename TSequenceV, typename TBand, typename TAlignmentProfile>
inline bool _isValidDPSettings(TSequenceH const & seqH,
TSequenceV const & seqV,
TBand const & band,
TAlignmentProfile const & alignProfile)
{
// Check if the sequences are empty.
if (empty(seqH) || empty(seqV))
{
return false;
}
return _checkBandProperties(seqH, seqV, band, alignProfile);
}
// ----------------------------------------------------------------------------
// Function _isBandEnabled()
// ----------------------------------------------------------------------------
// Returns true if a band is selected, otherwise false.
template <typename TBandSpec>
inline bool
_isBandEnabled(DPBandConfig<TBandSpec> const & /*band*/)
{
return IsSameType<TBandSpec, BandOn>::VALUE;
}
// ----------------------------------------------------------------------------
// Function _computeCell()
// ----------------------------------------------------------------------------
// Computes the score and tracks it if enabled.
template <typename TDPScout,
typename TTraceMatrixNavigator,
typename TDPCell,
typename TSequenceHValue, typename TSequenceVValue, typename TScoringScheme, typename TColumnDescriptor,
typename TCellDescriptor, typename TDPProfile>
inline void
_computeCell(TDPScout & scout,
TTraceMatrixNavigator & traceMatrixNavigator,
TDPCell & current,
TDPCell & diagonal,
TDPCell const & horizontal,
TDPCell & vertical,
TSequenceHValue const & seqHVal,
TSequenceVValue const & seqVVal,
TScoringScheme const & scoringScheme,
TColumnDescriptor const &,
TCellDescriptor const &, // One of FirstCell, InnerCell or LastCell.
TDPProfile const &)
{
typedef DPMetaColumn_<TDPProfile, TColumnDescriptor> TMetaColumn;
assignValue(traceMatrixNavigator,
_computeScore(current, diagonal, horizontal, vertical, seqHVal, seqVVal, scoringScheme,
typename RecursionDirection_<TMetaColumn, TCellDescriptor>::Type(),
TDPProfile()));
if (TrackingEnabled_<TMetaColumn, TCellDescriptor>::VALUE)
{
typedef typename LastColumnEnabled_<TDPProfile, TColumnDescriptor>::Type TIsLastColumn;
typedef typename LastRowEnabled_<TDPProfile, TCellDescriptor, TColumnDescriptor>::Type TIsLastRow;
// TODO(rrahn): Refactor to set vertical score only when max is updated.
if (IsTracebackEnabled_<TDPProfile>::VALUE)
{
_setVerticalScoreOfCell(current, _verticalScoreOfCell(vertical));
}
_scoutBestScore(scout, current, traceMatrixNavigator,
TIsLastColumn(), TIsLastRow());
}
}
// ----------------------------------------------------------------------------
// Function _precomputeScoreMatrixOffset()
// ----------------------------------------------------------------------------
// Default fallback if scoring scheme is not a matrix.
template <typename TSeqValue,
typename TScoringScheme>
inline TSeqValue const &
_precomputeScoreMatrixOffset(TSeqValue const & seqVal,
TScoringScheme const & /*score*/)
{
return seqVal;
}
// ----------------------------------------------------------------------------
// Function _computeTrack()
// ----------------------------------------------------------------------------
template <typename TDPScout,
typename TDPScoreMatrixNavigator,
typename TDPTraceMatrixNavigator,
typename TSeqHValue,
typename TSeqVValue,
typename TSeqVIterator,
typename TScoringScheme,
typename TDPCell,
typename TColumnDescriptor,
typename TDPProfile>
inline void
_computeTrack(TDPScout & scout,
TDPScoreMatrixNavigator & dpScoreMatrixNavigator,
TDPTraceMatrixNavigator & dpTraceMatrixNavigator,
TSeqHValue const & seqHValue,
TSeqVValue const & seqVValue,
TSeqVIterator const & seqBegin,
TSeqVIterator const & seqEnd,
TScoringScheme const & scoringScheme,
TDPCell & cacheDiag,
TDPCell & cacheVert,
TColumnDescriptor const &,
TDPProfile const &)
{
_goNextCell(dpScoreMatrixNavigator, TColumnDescriptor(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, TColumnDescriptor(), FirstCell());
_preInitCacheDiagonal(cacheDiag, dpScoreMatrixNavigator, TColumnDescriptor());
// Precompute the row of the scoring matrix for future look-ups.
TSeqHValue tmpSeqH = _precomputeScoreMatrixOffset(seqHValue, scoringScheme);
// Initilaize SIMD version with multiple end points.
_preInitScoutVertical(scout);
// Compute the first cell.
_computeCell(scout,
dpTraceMatrixNavigator,
value(dpScoreMatrixNavigator),
cacheDiag,
previousCellHorizontal(dpScoreMatrixNavigator),
cacheVert,
tmpSeqH,
seqVValue,
scoringScheme,
TColumnDescriptor(), FirstCell(), TDPProfile());
TSeqVIterator iter = seqBegin;
for (; iter != seqEnd - 1; ++iter)
{
_goNextCell(dpScoreMatrixNavigator, TColumnDescriptor(), InnerCell());
_goNextCell(dpTraceMatrixNavigator, TColumnDescriptor(), InnerCell());
_incVerticalPos(scout);
// Compute the inner cell.
if (SEQAN_UNLIKELY(_reachedVerticalEndPoint(scout, iter)))
{
_computeCell(scout,
dpTraceMatrixNavigator,
value(dpScoreMatrixNavigator),
cacheDiag,
previousCellHorizontal(dpScoreMatrixNavigator),
cacheVert,
tmpSeqH, sequenceEntryForScore(scoringScheme, container(iter), position(iter)),
scoringScheme, TColumnDescriptor(), LastCell(), TDPProfile());
_nextVerticalEndPos(scout);
}
else
{
_computeCell(scout,
dpTraceMatrixNavigator,
value(dpScoreMatrixNavigator),
cacheDiag,
previousCellHorizontal(dpScoreMatrixNavigator),
cacheVert,
tmpSeqH, sequenceEntryForScore(scoringScheme, container(iter), position(iter)),
scoringScheme, TColumnDescriptor(), InnerCell(), TDPProfile());
}
}
_goNextCell(dpScoreMatrixNavigator, TColumnDescriptor(), LastCell());
_goNextCell(dpTraceMatrixNavigator, TColumnDescriptor(), LastCell());
_incVerticalPos(scout);
_computeCell(scout,
dpTraceMatrixNavigator,
value(dpScoreMatrixNavigator),
cacheDiag,
previousCellHorizontal(dpScoreMatrixNavigator),
cacheVert,
tmpSeqH,
sequenceEntryForScore(scoringScheme, container(iter), position(iter)),
scoringScheme,
TColumnDescriptor(), LastCell(), TDPProfile());
}
template <typename TDPScout,
typename TDPScoreMatrixNavigator,
typename TDPTraceMatrixNavigator,
typename TSeqHValue,
typename TSeqVValue,
typename TSeqVIterator,
typename TScoringScheme,
typename TColumnDescriptor,
typename TDPProfile>
inline void
_computeTrack(TDPScout & scout,
TDPScoreMatrixNavigator & dpScoreMatrixNavigator,
TDPTraceMatrixNavigator & dpTraceMatrixNavigator,
TSeqHValue const & seqHValue,
TSeqVValue const & seqVValue,
TSeqVIterator const & seqBegin,
TSeqVIterator const & seqEnd,
TScoringScheme const & scoringScheme,
TColumnDescriptor const &,
TDPProfile const &)
{
using TDPCell = std::decay_t<decltype(value(dpScoreMatrixNavigator))>;
TDPCell cacheDiag;
TDPCell cacheVert;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator, seqHValue, seqVValue, seqBegin, seqEnd,
scoringScheme, cacheDiag, cacheVert, TColumnDescriptor{}, TDPProfile{});
}
// ----------------------------------------------------------------------------
// Function _computeUnbandedAlignmentHelperTerminate()
// ----------------------------------------------------------------------------
template <typename TDPCell, typename TSpec>
inline bool //TODO(C++11) constexpr
_computeAlignmentHelperCheckTerminate(DPScout_<TDPCell, TSpec > const & /**/)
{
return false;
}
template <typename TDPCell, typename TSpec>
inline bool
_computeAlignmentHelperCheckTerminate(DPScout_<TDPCell,Terminator_<TSpec> > const & s)
{
return _terminationCriteriumIsMet(s);
}
// ----------------------------------------------------------------------------
// Function _computeUnbandedAlignment()
// ----------------------------------------------------------------------------
// Computes the standard DP-algorithm.
template <typename TDPScout,
typename TDPScoreMatrixNavigator,
typename TDPTraceMatrixNavigator,
typename TSequenceH,
typename TSequenceV,
typename TScoringScheme,
typename TBand,
typename TAlignmentAlgo, typename TGapCosts, typename TTraceFlag, typename TExecPolicy>
inline void
_computeAlignmentImpl(TDPScout & scout,
TDPScoreMatrixNavigator & dpScoreMatrixNavigator,
TDPTraceMatrixNavigator & dpTraceMatrixNavigator,
TSequenceH const & seqH,
TSequenceV const & seqV,
TScoringScheme const & scoringScheme,
TBand const & /*band*/,
DPProfile_<TAlignmentAlgo, TGapCosts, TTraceFlag, TExecPolicy> const & dpProfile,
NavigateColumnWise const & /*tag*/)
{
typedef typename Iterator<TSequenceH const, Rooted>::Type TConstSeqHIterator;
typedef typename Iterator<TSequenceV const, Rooted>::Type TConstSeqVIterator;
// Initilaize SIMD version with multiple end points.
_preInitScoutHorizontal(scout);
// ============================================================================
// PREPROCESSING
// ============================================================================
TConstSeqVIterator seqVBegin = begin(seqV, Rooted());
TConstSeqVIterator seqVEnd = end(seqV, Rooted());
SEQAN_ASSERT_GT(length(seqH), 0u);
SEQAN_ASSERT_GT(length(seqV), 0u);
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInitialColumn, FullColumn>(), dpProfile);
// ============================================================================
// MAIN DP
// ============================================================================
TConstSeqHIterator seqHIter = begin(seqH, Rooted());
TConstSeqHIterator seqHIterEnd = end(seqH, Rooted()) - 1;
for (; seqHIter != seqHIterEnd; ++seqHIter)
{
_incHorizontalPos(scout);
// We might only select it if SIMD version is available.
if (SEQAN_UNLIKELY(_reachedHorizontalEndPoint(scout, seqHIter)))
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, FullColumn>(), dpProfile);
_nextHorizontalEndPos(scout);
}
else
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInnerColumn, FullColumn>(), dpProfile);
}
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
}
// ============================================================================
// POSTPROCESSING
// ============================================================================
_incHorizontalPos(scout);
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, FullColumn>(), dpProfile);
// If we compute only the single option. we need to check if there are other possibilities at the end.
// Traceback only from Diagonal, but could also come from vertical or horizontal.
// for (unsigned i = 0; i < length(recMatrix); ++i)
// {
// std::cout << recMatrix[i]._score << "\t";
// }
// std::cout << std::endl;
}
// ----------------------------------------------------------------------------
// Function _computeAlignment() banded
// ----------------------------------------------------------------------------
// Computes the banded DP-algorithm.
template <typename TDPScout,
typename TDPScoreMatrixNavigator,
typename TDPTraceMatrixNavigator,
typename TSequenceH,
typename TSequenceV,
typename TScoringScheme,
typename TBand,
typename TAlignmentAlgo, typename TGapCosts, typename TTraceFlag, typename TExecPolicy>
inline void
_computeAlignmentImpl(TDPScout & scout,
TDPScoreMatrixNavigator & dpScoreMatrixNavigator,
TDPTraceMatrixNavigator & dpTraceMatrixNavigator,
TSequenceH const & seqH,
TSequenceV const & seqV,
TScoringScheme const & scoringScheme,
TBand const & band,
DPProfile_<TAlignmentAlgo, TGapCosts, TTraceFlag, TExecPolicy> const & dpProfile,
NavigateColumnWiseBanded const & /*tag*/)
{
typedef DPProfile_<TAlignmentAlgo, TGapCosts, TTraceFlag> TDPProfile;
typedef typename MakeSigned<typename Size<TSequenceH>::Type>::Type TSignedSizeSeqH;
typedef typename MakeSigned<typename Size<TSequenceV>::Type>::Type TSignedSizeSeqV;
typedef typename Iterator<TSequenceH const, Rooted>::Type TConstSeqHIterator;
typedef typename Iterator<TSequenceV const, Rooted>::Type TConstSeqVIterator;
using TDPScoreValue = std::decay_t<decltype(value(dpScoreMatrixNavigator))>;
// Caching these cells improves performance significantly.
TDPScoreValue cacheDiag;
TDPScoreValue cacheVert;
if (upperDiagonal(band) == lowerDiagonal(band))
{
_computeHammingDistance(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator, seqH, seqV, scoringScheme, band,
dpProfile);
return;
}
// Now we have the problem of not knowing when we are in the last cell.
// ============================================================================
// PREPROCESSING
// ============================================================================
TSignedSizeSeqH seqHlength = static_cast<TSignedSizeSeqH>(length(seqH));
TSignedSizeSeqH seqVlength = static_cast<TSignedSizeSeqV>(length(seqV));
TConstSeqVIterator seqVBegin = begin(seqV, Rooted()) - _min(0, 1 + upperDiagonal(band));
TConstSeqVIterator seqVEnd = begin(seqV, Rooted()) - _min(0, _max(-seqVlength, lowerDiagonal(band)));
// We have to distinguish two band sizes. Some which spans the whole matrix in between and thus who not.
// This can be distinguished, if UpperDiagonal > length(seqV) + LowerDiagonal
// We start at least at the first position of the horizontal sequence or wherever the lower diagonal begins first.
TConstSeqHIterator seqHIterBegin = begin(seqH, Rooted()) + _max(0, _min(seqHlength - 1, lowerDiagonal(band)));
// The horizontal initial phase ends after the upper diagonal but at most after the horizontal sequence, or there is no horizontal initialization phase.
TConstSeqHIterator seqHIterEndColumnTop = begin(seqH, Rooted()) + _min(seqHlength - 1, _max(0, upperDiagonal(band)));
// The middle band phase ends after the lower diagonal crosses the bottom of the alignment matrix or after the horizontal sequence if it is smaller.
TConstSeqHIterator seqHIterEndColumnMiddle = begin(seqH, Rooted()) + _min(seqHlength - 1, _max(0, seqVlength + lowerDiagonal(band)));
// Swap the two iterators if we are in a band that spans over the full column.
if (upperDiagonal(band) > seqVlength + lowerDiagonal(band))
std::swap(seqHIterEndColumnTop, seqHIterEndColumnMiddle);
// The bottom band phase ends after the upper diagonal of the band crosses the bottom of the matrix or after the horizontal sequence if it is smaller.
TConstSeqHIterator seqHIterEndColumnBottom = begin(seqH, Rooted()) + _max(0, _min(seqHlength,
upperDiagonal(band) + seqVlength) - 1);
// The Initial column can be PartialColumnTop which is given if the upper diagonal is >= 0,
// otherwise it only can be PartialColumnMiddle or PartialColumnBottom depending where the lower diagonal is.
// Check for single initialization cells in InitialColumn and FinalColumn.
if (seqHIterBegin == end(seqH, Rooted()) - 1)
{
// Set the iterator to the begin of the track.
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell());
// Only one cell
_computeCell(scout, dpTraceMatrixNavigator, value(dpScoreMatrixNavigator),
cacheDiag, previousCellHorizontal(dpScoreMatrixNavigator), cacheVert,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIterBegin)),
sequenceEntryForScore(scoringScheme, seqV, 0), scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell(), TDPProfile());
// we might need to additionally track this point.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, PartialColumnTop> >, FirstCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, True(), False());
return;
}
if (seqHIterEndColumnBottom == begin(seqH, Rooted()))
{
// Set the iterator to the begin of the track.
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnBottom>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnBottom>(), FirstCell());
// Only one cell
_computeCell(scout, dpTraceMatrixNavigator, value(dpScoreMatrixNavigator),
cacheDiag, previousCellHorizontal(dpScoreMatrixNavigator), cacheVert,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin)), scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnBottom>(), FirstCell(), TDPProfile());
// We might need to additionally track this point.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInitialColumn, PartialColumnBottom> >, LastCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, False(), True());
return;
}
if (upperDiagonal(band) < 0)
{
++seqVBegin;
if (lowerDiagonal(band) > -seqVlength)
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnMiddle>(), dpProfile);
else
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnBottom>(), dpProfile);
}
else if (lowerDiagonal(band) >= 0)
{
// Set the iterator to the begin of the track.
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell());
//TODO(rrahn): We possibly need to set the cache values here?
// Should we not just compute the cell?
_computeCell(scout, dpTraceMatrixNavigator, value(dpScoreMatrixNavigator),
cacheDiag, previousCellHorizontal(dpScoreMatrixNavigator), cacheVert,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIterBegin)),
sequenceEntryForScore(scoringScheme, seqV, 0),
scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), FirstCell(), TDPProfile());
// we might need to additionally track this point.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, PartialColumnTop> >, FirstCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, False(), False());
}
else // Upper diagonal >= 0 and lower Diagonal < 0
{
if (lowerDiagonal(band) <= -seqVlength) // The band is bounded by the top and bottom of the matrix.
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInitialColumn, FullColumn>(), dpProfile);
}
else // The band is bounded by the top but not the bottom of the matrix.
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, 0),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInitialColumn, PartialColumnTop>(), dpProfile);
}
}
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
// ============================================================================
// MAIN DP
// ============================================================================
TConstSeqHIterator seqHIter = seqHIterBegin;
// Compute the first part of the band, where the band is bounded by the top but not by the bottom of the matrix.
for (; seqHIter != seqHIterEndColumnTop; ++seqHIter)
{
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInnerColumn, PartialColumnTop>(), dpProfile);
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
}
// TODO(rmaerker): Check if putting the if-statement before the actual algorithm can speedup the code.
// Check whether the band spans over the full column or not at some point.
if (upperDiagonal(band) > seqVlength + lowerDiagonal(band))
{
// Compute the second part of the band, where the band is bounded by the top and the bottom of the matrix.
// We might want to track the current cell here, since this is the first cell that crosses the bottom but is
// not part of the FullColumn tracks.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, FullColumn> >, LastCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, False(), True());
for (; seqHIter != seqHIterEndColumnMiddle; ++seqHIter)
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInnerColumn, FullColumn>(), dpProfile);
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
}
}
else // Compute the second part of the band, where the band is not bounded by the top and bottom of the matrix
{
for (; seqHIter != seqHIterEndColumnMiddle; ++seqHIter)
{
++seqVBegin;
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInnerColumn, PartialColumnMiddle>(), dpProfile);
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
} // We might want to track the current cell here, since this is the first cell that crosses the bottom.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom> >, LastCell>::VALUE)
{
if (lowerDiagonal(band) + seqVlength < seqHlength)
{
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, False(), True());
}
}
}
// Compute the third part of the band, where the band, is bounded by the bottom but not by the top of the matrix.
for (; seqHIter != seqHIterEndColumnBottom; ++seqHIter)
{
++seqVBegin;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom>(), dpProfile);
if (_computeAlignmentHelperCheckTerminate(scout))
{
return;
}
}
// ============================================================================
// POSTPROCESSING
// ============================================================================
// Check where the last track of the column is located.
if (seqHIter - begin(seqH, Rooted()) < seqHlength - 1) // Case 1: The band ends before the final column is reached.
{
// Set the iterator to the begin of the track.
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom>(), FirstCell());
_preInitCacheDiagonal(cacheDiag, dpScoreMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom>());
_computeCell(scout, dpTraceMatrixNavigator, value(dpScoreMatrixNavigator),
cacheDiag, previousCellHorizontal(dpScoreMatrixNavigator), cacheVert,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin)),
scoringScheme,
MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom>(), FirstCell(), TDPProfile());
// We might need to additionally track this point.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, PartialColumnBottom> >, LastCell>::VALUE)
{
_setVerticalScoreOfCell(value(dpScoreMatrixNavigator), _verticalScoreOfCell(cacheVert));
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, False(), True());
}
}
else if (seqHIter == end(seqH, Rooted()) - 1) // Case 2: The band ends somewhere in the final column of the matrix.
{
// Case2a: The band ends in the last cell of the final column.
if (upperDiagonal(band) == seqHlength - seqVlength)
{
// Set the iterator to the begin of the track.
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom>(), FirstCell());
_preInitCacheDiagonal(cacheDiag, dpScoreMatrixNavigator, MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom>());
_computeCell(scout, dpTraceMatrixNavigator, value(dpScoreMatrixNavigator),
cacheDiag, previousCellHorizontal(dpScoreMatrixNavigator), cacheVert,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin)),
scoringScheme,
MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom>(), FirstCell(), TDPProfile());
// we might need to additionally track this point.
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom> >, LastCell>::VALUE)
{
_setVerticalScoreOfCell(value(dpScoreMatrixNavigator), _verticalScoreOfCell(cacheVert));
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, True(), True());
}
}
else // Case2b: At least two cells intersect between the band and the matrix in the final column of the matrix.
{
if (upperDiagonal(band) >= seqHlength) // The band is bounded by the top of the matrix only or by the top and the bottom.
{
if (lowerDiagonal(band) + seqVlength > seqHlength) // The band is bounded by the top of the matrix
{
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, PartialColumnTop>(), dpProfile);
}
else // The band is bounded by the top and the bottom of the matrix.
{
if (lowerDiagonal(band) + seqVlength + 1 > seqHlength) // We have to go into the last cell.
{
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme, cacheDiag, cacheVert,
MetaColumnDescriptor<DPFinalColumn, PartialColumnTop>(), dpProfile);
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, FullColumn> >, LastCell>::VALUE)
{
_setVerticalScoreOfCell(value(dpScoreMatrixNavigator), _verticalScoreOfCell(cacheVert));
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, True(), True());
}
}
else
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, 0),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, FullColumn>(), dpProfile);
}
}
else // The band is bounded by bottom of matrix or completely unbounded.
{
++seqVBegin;
if (lowerDiagonal(band) + seqVlength <= seqHlength) // The band is bounded by the bottom of the matrix.
{
if (lowerDiagonal(band) + seqVlength == seqHlength) // We have to go into the last cell.
{
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme, cacheDiag, cacheVert,
MetaColumnDescriptor<DPFinalColumn, PartialColumnMiddle>(), dpProfile);
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom> >, LastCell>::VALUE)
{
_setVerticalScoreOfCell(value(dpScoreMatrixNavigator), _verticalScoreOfCell(cacheVert));
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator, True(), True());
}
}
else
{
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, PartialColumnBottom>(), dpProfile);
}
}
else // The band is unbounded by the matrix.
{
++seqVEnd;
_computeTrack(scout, dpScoreMatrixNavigator, dpTraceMatrixNavigator,
sequenceEntryForScore(scoringScheme, seqH, position(seqHIter)),
sequenceEntryForScore(scoringScheme, seqV, position(seqVBegin) - 1),
seqVBegin, seqVEnd, scoringScheme,
MetaColumnDescriptor<DPFinalColumn, PartialColumnMiddle>(), dpProfile);
}
}
}
}
}
// ----------------------------------------------------------------------------
// Function _computeHammingDistance()
// ----------------------------------------------------------------------------
// Computes the Hamming-Distance if the band-size is 1.
template <typename TDPScout,
typename TDPScoreMatrixNavigator,
typename TDPTraceMatrixNavigator,
typename TSequenceH,
typename TSequenceV,
typename TScoringScheme,
typename TBand,
typename TAlignmentAlgo, typename TGapCosts, typename TTraceFlag, typename TExecPolicy>
inline void
_computeHammingDistance(TDPScout & scout,
TDPScoreMatrixNavigator & dpScoreMatrixNavigator,
TDPTraceMatrixNavigator & dpTraceMatrixNavigator,
TSequenceH const & seqH,
TSequenceV const & seqV,
TScoringScheme const & scoringScheme,
TBand const & band,
DPProfile_<TAlignmentAlgo, TGapCosts, TTraceFlag, TExecPolicy> const &)
{
typedef typename MakeSigned<typename Size<TSequenceH const>::Type>::Type TSignedSizeSeqH;
typedef typename MakeSigned<typename Size<TSequenceV const>::Type>::Type TSignedSizeSeqV;
typedef typename Iterator<TSequenceH const, Rooted>::Type TConstSeqHIterator;
typedef typename Iterator<TSequenceV const, Rooted>::Type TConstSeqVIterator;
typedef typename Value<TDPScoreMatrixNavigator>::Type TDPCell;
typedef DPProfile_<TAlignmentAlgo, TGapCosts, TTraceFlag> TDPProfile;
// ============================================================================
// PREPROCESSING
// ============================================================================
TSignedSizeSeqH seqHlength = static_cast<TSignedSizeSeqH>(length(seqH));
TSignedSizeSeqH seqVlength = static_cast<TSignedSizeSeqV>(length(seqV));
TConstSeqHIterator itH = begin(seqH, Rooted()) + _max(0, _min(seqHlength - 1, upperDiagonal(band)));
TConstSeqHIterator itHEnd = begin(seqH, Rooted()) + _min(seqHlength - 1, upperDiagonal(band) + seqVlength);
TConstSeqVIterator itV = begin(seqV, Rooted()) + _max(0, _min(seqVlength - 1, -lowerDiagonal(band)));
TConstSeqVIterator itVEnd = begin(seqV, Rooted()) + _min(seqVlength - 1, lowerDiagonal(band) + seqHlength);
TDPCell dummy;
assignValue(dpTraceMatrixNavigator,
_computeScore(value(dpScoreMatrixNavigator), dummy, dummy, dummy,
sequenceEntryForScore(scoringScheme, seqH, position(itH)),
sequenceEntryForScore(scoringScheme, seqV, position(itV)),
scoringScheme, RecursionDirectionZero(), TDPProfile()));
if (upperDiagonal(band) < 0)
{
if (upperDiagonal(band) == -seqVlength)
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInitialColumn, FullColumn> >, LastCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
return;
}
else
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInitialColumn, FullColumn> >, InnerCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
}
else if (lowerDiagonal(band) > 0)
{
if (lowerDiagonal(band) == seqHlength)
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, FullColumn> >, FirstCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
return;
}
else
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, FullColumn> >, FirstCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
}
else
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInitialColumn, FullColumn> >, FirstCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
TDPCell prevDiagonal = value(dpScoreMatrixNavigator);
// ============================================================================
// MAIN DP
// ============================================================================
while (itH != itHEnd && itV != itVEnd)
{
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, FullColumn>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, FullColumn>(), FirstCell());
assignValue(dpTraceMatrixNavigator,
_computeScore(value(dpScoreMatrixNavigator), prevDiagonal, dummy, dummy,
sequenceEntryForScore(scoringScheme, seqH, position(itH)),
sequenceEntryForScore(scoringScheme, seqV, position(itV)),
scoringScheme, RecursionDirectionDiagonal(), TDPProfile()));
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, FullColumn> >, InnerCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
prevDiagonal = value(dpScoreMatrixNavigator);
++itH;
++itV;
}
// ============================================================================
// POSTPROCESSING
// ============================================================================
_goNextCell(dpScoreMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, FullColumn>(), FirstCell());
_goNextCell(dpTraceMatrixNavigator, MetaColumnDescriptor<DPInnerColumn, FullColumn>(), FirstCell());
assignValue(dpTraceMatrixNavigator,
_computeScore(value(dpScoreMatrixNavigator), prevDiagonal, dummy, dummy,
sequenceEntryForScore(scoringScheme, seqH, position(itH)),
sequenceEntryForScore(scoringScheme, seqV, position(itV)),
scoringScheme, RecursionDirectionDiagonal(), TDPProfile()));
if (itH == itHEnd)
{
if (itV == itVEnd) // Is in the last cell of final column
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, FullColumn> >, LastCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPInnerColumn, FullColumn> >, LastCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
else
{
if (TrackingEnabled_<DPMetaColumn_<TDPProfile, MetaColumnDescriptor<DPFinalColumn, FullColumn> >, InnerCell>::VALUE)
_scoutBestScore(scout, value(dpScoreMatrixNavigator), dpTraceMatrixNavigator);
}
}
// ----------------------------------------------------------------------------
// Function _printScoreMatrix()
// ----------------------------------------------------------------------------
template <typename TTraceMatrix>
void _printScoreMatrix(TTraceMatrix & scoreMatrix)
{
typedef typename Size<TTraceMatrix>::Type TSize;
TSize dimH = length(scoreMatrix, +DPMatrixDimension_::HORIZONTAL);
TSize dimV = length(scoreMatrix, +DPMatrixDimension_::VERTICAL);
for (unsigned row = 0; row < dimV; ++row)
{
for (unsigned column = 0; column < dimH; ++column)
if (_scoreOfCell(value(scoreMatrix, row + column * dimV)) <= DPCellDefaultInfinity<DPCell_<int, LinearGaps>>::VALUE)
std::cout << "-∞\t";
else
std::cout << _scoreOfCell(value(scoreMatrix, row + column * dimV)) << "\t";
std::cout << std::endl;
}
std::cout << std::endl;
}
// ----------------------------------------------------------------------------
// Function _printTracebackMatrix()
// ----------------------------------------------------------------------------
template <typename TTraceMatrix>
void _printTracebackMatrix(TTraceMatrix & dpTraceMatrix)
{
typedef typename Size<TTraceMatrix>::Type TSize;
TSize dimH = length(dpTraceMatrix, +DPMatrixDimension_::HORIZONTAL);
TSize dimV = length(dpTraceMatrix, +DPMatrixDimension_::VERTICAL);
for (unsigned row = 0; row < dimV; ++row)
{
for (unsigned column = 0; column < dimH; ++column)
std::cout << _translateTraceValue(value(dpTraceMatrix, row + column * dimV)) << "\t";
std::cout << std::endl;
}
std::cout << std::endl;
}
template <typename TTraceMatrix, typename TPosition>
void _printTracebackMatrix(TTraceMatrix & dpTraceMatrix, TPosition const simdLane)
{
typedef typename Size<TTraceMatrix>::Type TSize;
TSize dimH = length(dpTraceMatrix, +DPMatrixDimension_::HORIZONTAL);
TSize dimV = length(dpTraceMatrix, +DPMatrixDimension_::VERTICAL);
for (unsigned row = 0; row < dimV; ++row)
{
for (unsigned column = 0; column < dimH; ++column)
std::cout << _translateTraceValue(value(dpTraceMatrix, row + column * dimV)[simdLane]) << "\t";
std::cout << std::endl;
}
std::cout << std::endl;
}
// ----------------------------------------------------------------------------
// Function _correctTraceValue()
// ----------------------------------------------------------------------------
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TScoreValue> > >, void)
_correctTraceValue(TTraceNavigator &,
DPScout_<DPCell_<TScoreValue, LinearGaps>, TDPScoutSpec> const &)
{
// Nothing to do.
}
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TScoreValue> >, void)
_correctTraceValue(TTraceNavigator &,
DPScout_<DPCell_<TScoreValue, LinearGaps>, TDPScoutSpec> const &)
{
// Nothing to do.
}
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TScoreValue> > >, void)
_correctTraceValue(TTraceNavigator & traceNavigator,
DPScout_<DPCell_<TScoreValue, AffineGaps>, TDPScoutSpec> const & dpScout)
{
_setToPosition(traceNavigator, maxHostPosition(dpScout));
if (_verticalScoreOfCell(dpScout._maxScore) == _scoreOfCell(dpScout._maxScore))
{
value(traceNavigator) &= ~TraceBitMap_<TScoreValue>::DIAGONAL;
value(traceNavigator) |= TraceBitMap_<TScoreValue>::MAX_FROM_VERTICAL_MATRIX;
}
else if (_horizontalScoreOfCell(dpScout._maxScore) == _scoreOfCell(dpScout._maxScore))
{
value(traceNavigator) &= ~TraceBitMap_<TScoreValue>::DIAGONAL;
value(traceNavigator) |= TraceBitMap_<TScoreValue>::MAX_FROM_HORIZONTAL_MATRIX;
}
}
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TScoreValue> >, void)
_correctTraceValue(TTraceNavigator & traceNavigator,
DPScout_<DPCell_<TScoreValue, AffineGaps>, TDPScoutSpec> const & dpScout)
{
using TMaskType = typename SimdMaskVector<TScoreValue>::Type;
_setToPosition(traceNavigator, toGlobalPosition(traceNavigator,
maxHostCoordinate(dpScout, +DPMatrixDimension_::HORIZONTAL),
maxHostCoordinate(dpScout, +DPMatrixDimension_::VERTICAL)));
TMaskType flag = createVector<TMaskType>(0);
assignValue(flag, dpScout._simdLane, -1);
auto cmpV = cmpEq(_verticalScoreOfCell(dpScout._maxScore), _scoreOfCell(dpScout._maxScore)) & flag;
auto cmpH = cmpEq(_horizontalScoreOfCell(dpScout._maxScore), _scoreOfCell(dpScout._maxScore)) & flag;
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) & ~TraceBitMap_<TScoreValue>::DIAGONAL,
cmpV | cmpH);
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) | TraceBitMap_<TScoreValue>::MAX_FROM_VERTICAL_MATRIX,
cmpV);
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) | TraceBitMap_<TScoreValue>::MAX_FROM_HORIZONTAL_MATRIX,
cmpH);
}
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Not<Is<SimdVectorConcept<TScoreValue> > >, void)
_correctTraceValue(TTraceNavigator & traceNavigator,
DPScout_<DPCell_<TScoreValue, DynamicGaps>, TDPScoutSpec> const & dpScout)
{
_setToPosition(traceNavigator, maxHostPosition(dpScout));
if (isGapExtension(dpScout._maxScore, DynamicGapExtensionVertical()))
{
value(traceNavigator) &= ~TraceBitMap_<TScoreValue>::DIAGONAL;
value(traceNavigator) |= TraceBitMap_<TScoreValue>::MAX_FROM_VERTICAL_MATRIX;
}
else if (isGapExtension(dpScout._maxScore, DynamicGapExtensionHorizontal()))
{
value(traceNavigator) &= ~TraceBitMap_<TScoreValue>::DIAGONAL;
value(traceNavigator) |= TraceBitMap_<TScoreValue>::MAX_FROM_HORIZONTAL_MATRIX;
}
}
template <typename TTraceNavigator, typename TScoreValue, typename TDPScoutSpec>
inline SEQAN_FUNC_ENABLE_IF(Is<SimdVectorConcept<TScoreValue> >, void)
_correctTraceValue(TTraceNavigator & traceNavigator,
DPScout_<DPCell_<TScoreValue, DynamicGaps>, TDPScoutSpec> const & dpScout)
{
using TMaskType = typename SimdMaskVector<TScoreValue>::Type;
_setToPosition(traceNavigator, maxHostPosition(dpScout));
TMaskType flag = createVector<TMaskType>(0);
assignValue(flag, dpScout._simdLane, -1);
auto cmpV = isGapExtension(dpScout._maxScore, DynamicGapExtensionVertical()) & flag;
auto cmpH = isGapExtension(dpScout._maxScore, DynamicGapExtensionHorizontal()) & flag;
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) & ~TraceBitMap_<TScoreValue>::DIAGONAL,
cmpV | cmpH);
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) | TraceBitMap_<TScoreValue>::MAX_FROM_VERTICAL_MATRIX,
cmpV);
value(traceNavigator) = blend(value(traceNavigator),
value(traceNavigator) | TraceBitMap_<TScoreValue>::MAX_FROM_HORIZONTAL_MATRIX,
cmpH);
}
// ----------------------------------------------------------------------------
// Function _finishAlignment()
// ----------------------------------------------------------------------------
template <typename TTraceTarget,
typename TTraceMatNavigator,
typename TScoreValue, typename TGapsModel, typename TDPScoutSpec,
typename TSeqH,
typename TSeqV,
typename TBandSwitch,
typename TAlignmentAlgorithm, typename TGapScheme, typename TTraceFlag, typename TExecPolicy>
inline SEQAN_FUNC_ENABLE_IF(Not<IsTracebackEnabled_<TTraceFlag> >, TScoreValue)
_finishAlignment(TTraceTarget & /*traceSegments*/,
TTraceMatNavigator & /*dpTraceMatrixNavigator*/,
DPScout_<DPCell_<TScoreValue, TGapsModel>, TDPScoutSpec> & dpScout,
TSeqH const & /*seqH*/,
TSeqV const & /*seqV*/,
DPBandConfig<TBandSwitch> const & /*band*/,
DPProfile_<TAlignmentAlgorithm, TGapScheme, TTraceFlag, TExecPolicy> const & /*dpProfile*/)
{
return maxScore(dpScout);
}
template <typename TTraceTarget,
typename TTraceMatNavigator,
typename TScoreValue, typename TGapsModel, typename TDPScoutSpec,
typename TSeqH,
typename TSeqV,
typename TBandSwitch,
typename TAlignmentAlgorithm, typename TGapScheme, typename TTraceFlag, typename TExecPolicy>
inline SEQAN_FUNC_ENABLE_IF(And<Is<SimdVectorConcept<TScoreValue> >, IsTracebackEnabled_<TTraceFlag> >, TScoreValue)
_finishAlignment(TTraceTarget & traceSegments,
TTraceMatNavigator & dpTraceMatrixNavigator,
DPScout_<DPCell_<TScoreValue, TGapsModel>, TDPScoutSpec> & scout,
TSeqH const & seqH,
TSeqV const & seqV,
DPBandConfig<TBandSwitch> const & band,
DPProfile_<TAlignmentAlgorithm, TGapScheme, TTraceFlag, TExecPolicy> const & dpProfile)
{
typedef typename Size<TTraceTarget>::Type TSize;
for(TSize i = 0; i < length(traceSegments); ++i)
{
_setSimdLane(dpTraceMatrixNavigator, i);
_setSimdLane(scout, i);
if (IsSingleTrace_<TTraceFlag>::VALUE)
{
_correctTraceValue(dpTraceMatrixNavigator, scout);
}
_computeTraceback(traceSegments[i], dpTraceMatrixNavigator,
toGlobalPosition(dpTraceMatrixNavigator,
maxHostCoordinate(scout, +DPMatrixDimension_::HORIZONTAL),
maxHostCoordinate(scout, +DPMatrixDimension_::VERTICAL)),
_hostLengthH(scout, seqH),
_hostLengthV(scout, seqV), band, dpProfile);
}
return maxScore(scout);
}
template <typename TTraceTarget,
typename TTraceMatNavigator,
typename TScoreValue, typename TGapsModel, typename TDPScoutSpec,
typename TSeqH,
typename TSeqV,
typename TBandSwitch,
typename TAlignmentAlgorithm, typename TGapScheme, typename TTraceFlag, typename TExecPolicy>
inline SEQAN_FUNC_ENABLE_IF(And<Not<Is<SimdVectorConcept<TScoreValue> > >, IsTracebackEnabled_<TTraceFlag> >, TScoreValue)
_finishAlignment(TTraceTarget & traceSegments,
TTraceMatNavigator & dpTraceMatrixNavigator,
DPScout_<DPCell_<TScoreValue, TGapsModel>, TDPScoutSpec> & dpScout,
TSeqH const & seqH,
TSeqV const & seqV,
DPBandConfig<TBandSwitch> const & band,
DPProfile_<TAlignmentAlgorithm, TGapScheme, TTraceFlag, TExecPolicy> const & dpProfile)
{
if (IsSingleTrace_<TTraceFlag>::VALUE)
_correctTraceValue(dpTraceMatrixNavigator, dpScout);
_computeTraceback(traceSegments, dpTraceMatrixNavigator, dpScout, seqH, seqV, band, dpProfile);
return maxScore(dpScout);
}
// ----------------------------------------------------------------------------
// Function _computeAligmnment()
// ----------------------------------------------------------------------------
template <typename TDPScoreValue, typename TTraceValue, typename TScoreMatHost, typename TTraceMatHost,
typename TTraceTarget,
typename TScoutState,
typename TSequenceH,
typename TSequenceV,
typename TScoreScheme,
typename TBandSwitch,
typename TAlignmentAlgorithm, typename TGapScheme, typename TTraceFlag, typename TExecPolicy>
inline typename Value<TScoreScheme>::Type
_computeAlignment(DPContext<TDPScoreValue, TTraceValue, TScoreMatHost, TTraceMatHost> & dpContext,
TTraceTarget & traceSegments,
TScoutState & scoutState,
TSequenceH const & seqH,
TSequenceV const & seqV,
TScoreScheme const & scoreScheme,
DPBandConfig<TBandSwitch> const & band,
DPProfile_<TAlignmentAlgorithm, TGapScheme, TTraceFlag, TExecPolicy> const & dpProfile)
{
typedef typename DefaultScoreMatrixSpec_<TAlignmentAlgorithm>::Type TScoreMatrixSpec;
typedef DPMatrix_<TDPScoreValue, TScoreMatrixSpec, TScoreMatHost> TDPScoreMatrix;
typedef DPMatrix_<TTraceValue, FullDPMatrix, TTraceMatHost> TDPTraceMatrix;
using TNavigationSpec = std::conditional_t<std::is_same<TBandSwitch, BandOff>::value,
NavigateColumnWise,
NavigateColumnWiseBanded>;
typedef DPMatrixNavigator_<TDPScoreMatrix, DPScoreMatrix, TNavigationSpec> TDPScoreMatrixNavigator;
typedef DPMatrixNavigator_<TDPTraceMatrix, DPTraceMatrix<TTraceFlag>, TNavigationSpec> TDPTraceMatrixNavigator;
typedef typename ScoutSpecForAlignmentAlgorithm_<TAlignmentAlgorithm, TScoutState>::Type TDPScoutSpec;
typedef DPScout_<TDPScoreValue, TDPScoutSpec> TDPScout;
typedef typename Value<TScoreScheme>::Type TScoreValue;
// Check if current dp settings are valid. If not return infinity value for dp score value.
if (!_isValidDPSettings(seqH, seqV, band, dpProfile))
return createVector<TScoreValue>(std::numeric_limits<typename Value<TScoreValue>::Type>::min()); // NOTE(rrahn): In case of non-simd version, createVector returns just a scalar.
TDPScoreMatrix dpScoreMatrix;
TDPTraceMatrix dpTraceMatrix;
// TODO(rmaerker): Check whether the matrix allocation can be reduced if upperDiagonal < 0?
setLength(dpScoreMatrix, +DPMatrixDimension_::HORIZONTAL, length(seqH) + 1 - std::max(0, lowerDiagonal(band)));
setLength(dpTraceMatrix, +DPMatrixDimension_::HORIZONTAL, length(seqH) + 1 - std::max(0, lowerDiagonal(band)));
SEQAN_IF_CONSTEXPR (IsSameType<TBandSwitch, BandOff>::VALUE)
{
setLength(dpScoreMatrix, +DPMatrixDimension_::VERTICAL, length(seqV) + 1);
setLength(dpTraceMatrix, +DPMatrixDimension_::VERTICAL, length(seqV) + 1);
}
else
{
int bandSize = _min(static_cast<int>(length(seqH)), upperDiagonal(band)) - _max(lowerDiagonal(band), -static_cast<int>(length(seqV))) + 1;
setLength(dpScoreMatrix, +DPMatrixDimension_::VERTICAL, _min(static_cast<int>(length(seqV)) + 1, bandSize));
setLength(dpTraceMatrix, +DPMatrixDimension_::VERTICAL, _min(static_cast<int>(length(seqV)) + 1, bandSize));
}
// We set the host to the score matrix and the dp matrix.
setHost(dpScoreMatrix, getDpScoreMatrix(dpContext));
setHost(dpTraceMatrix, getDpTraceMatrix(dpContext));
resize(dpScoreMatrix);
// We do not need to allocate the memory for the trace matrix if the traceback is disabled.
if (IsTracebackEnabled_<TTraceFlag>::VALUE)
resize(dpTraceMatrix);
TDPScoreMatrixNavigator dpScoreMatrixNavigator{dpScoreMatrix, band};
TDPTraceMatrixNavigator dpTraceMatrixNavigator{dpTraceMatrix, band};
TDPScout dpScout(scoutState);
#if SEQAN_ALIGN_SIMD_PROFILE
profile.preprTimer += sysTime() - timer;
timer = sysTime();
#endif
// Execute the alignment.
_computeAlignmentImpl(dpScout, dpScoreMatrixNavigator, dpTraceMatrixNavigator, seqH, seqV, scoreScheme, band,
dpProfile, TNavigationSpec{});
#if SEQAN_ALIGN_SIMD_PROFILE
profile.alignTimer += sysTime() - timer;
timer = sysTime();
#endif
return _finishAlignment(traceSegments, dpTraceMatrixNavigator, dpScout, seqH, seqV, band, dpProfile);
}
} // namespace seqan
#endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_DP_ALGORITHM_IMPL_H_
