19 #define ALIGNMENT_DEBUG
20 #undef ALIGNMENT_DEBUG
61 template <
typename SpectrumType1,
typename SpectrumType2>
62 void getSpectrumAlignment(std::vector<std::pair<Size, Size> >& alignment,
const SpectrumType1& s1,
const SpectrumType2& s2)
const
64 if (!s1.isSorted() || !s2.isSorted())
71 double tolerance = (double)param_.getValue(
"tolerance");
73 if (!param_.getValue(
"is_relative_tolerance").toBool() )
75 std::map<Size, std::map<Size, std::pair<Size, Size> > > traceback;
76 std::map<Size, std::map<Size, double> > matrix;
80 for (
Size i = 1; i <= s1.size(); ++i)
82 matrix[i][0] = i * tolerance;
83 traceback[i][0] = std::make_pair(i - 1, 0);
85 for (
Size j = 1; j <= s2.size(); ++j)
87 matrix[0][j] = j * tolerance;
88 traceback[0][j] = std::make_pair(0, j - 1);
93 Size last_i(0), last_j(0);
96 for (
Size i = 1; i <= s1.size(); ++i)
98 double pos1(s1[i - 1].getMZ());
100 for (
Size j = left_ptr; j <= s2.size(); ++j)
102 bool off_band(
false);
104 double pos2(s2[j - 1].getMZ());
105 double diff_align = fabs(pos1 - pos2);
108 if (pos2 > pos1 && diff_align > tolerance)
110 if (i < s1.size() && j < s2.size() && s1[i].getMZ() < pos2)
117 if (pos1 > pos2 && diff_align > tolerance && j > left_ptr + 1)
122 double score_align = diff_align;
124 if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j - 1) != matrix[i - 1].end())
126 score_align += matrix[i - 1][j - 1];
130 score_align += (i - 1 + j - 1) * tolerance;
133 double score_up = tolerance;
134 if (matrix.find(i) != matrix.end() && matrix[i].find(j - 1) != matrix[i].end())
136 score_up += matrix[i][j - 1];
140 score_up += (i + j - 1) * tolerance;
143 double score_left = tolerance;
144 if (matrix.find(i - 1) != matrix.end() && matrix[i - 1].find(j) != matrix[i - 1].end())
146 score_left += matrix[i - 1][j];
150 score_left += (i - 1 + j) * tolerance;
153 #ifdef ALIGNMENT_DEBUG
154 cerr << i <<
" " << j <<
" " << left_ptr <<
" " << pos1 <<
" " << pos2 <<
" " << score_align <<
" " << score_left <<
" " << score_up << endl;
157 if (score_align <= score_up && score_align <= score_left && diff_align <= tolerance)
159 matrix[i][j] = score_align;
160 traceback[i][j] = std::make_pair(i - 1, j - 1);
166 if (score_up <= score_left)
168 matrix[i][j] = score_up;
169 traceback[i][j] = std::make_pair(i, j - 1);
173 matrix[i][j] = score_left;
174 traceback[i][j] = std::make_pair(i - 1, j);
190 #ifdef ALIGNMENT_DEBUG
192 cerr <<
"TheMatrix: " << endl <<
" \t \t";
193 for (
Size j = 0; j != s2.size(); ++j)
195 cerr << s2[j].getPosition()[0] <<
" \t";
198 for (
Size i = 0; i <= s1.size(); ++i)
202 cerr << s1[i - 1].getPosition()[0] <<
" \t";
208 for (
Size j = 0; j <= s2.size(); ++j)
210 if (matrix.has(i) && matrix[i].has(j))
212 if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
218 if (traceback[i][j].first == i - 1 && traceback[i][j].second == j)
228 cerr << matrix[i][j] <<
" \t";
244 while (i >= 1 && j >= 1)
246 if (traceback[i][j].first == i - 1 && traceback[i][j].second == j - 1)
248 alignment.push_back(std::make_pair(i - 1, j - 1));
250 Size new_i = traceback[i][j].first;
251 Size new_j = traceback[i][j].second;
259 #ifdef ALIGNMENT_DEBUG
262 cerr <<
"Alignment (size=" << alignment.size() <<
"): " << endl;
264 Size i_s1(0), i_s2(0);
265 for (vector<pair<Size, Size> >::const_reverse_iterator it = alignment.rbegin(); it != alignment.rend(); ++it, ++i_s1, ++i_s2)
267 while (i_s1 < it->first - 1)
269 cerr << i_s1 <<
" " << s1[i_s1].getPosition()[0] <<
" " << s1[i_s1].getIntensity() << endl;
272 while (i_s2 < it->second - 1)
274 cerr <<
" \t " << i_s2 <<
" " << s2[i_s2].getPosition()[0] <<
" " << s2[i_s2].getIntensity() << endl;
277 cerr <<
"(" << s1[it->first - 1].getPosition()[0] <<
" <-> " << s2[it->second - 1].getPosition()[0] <<
") ("
278 << it->first <<
"|" << it->second <<
") (" << s1[it->first - 1].getIntensity() <<
"|" << s2[it->second - 1].getIntensity() <<
")" << endl;
287 for (; it != it.
end(); ++it) alignment.emplace_back(it.
refIdx(), it.
tgtIdx());
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:66
A method or algorithm argument contains illegal values.
Definition: Exception.h:629
For each element in the reference container the closest peak in the target will be searched....
Definition: MatchedIterator.h:46
size_t refIdx() const
index into reference container
Definition: MatchedIterator.h:162
size_t tgtIdx() const
index into target container
Definition: MatchedIterator.h:168
static MatchedIterator end()
the end iterator
Definition: MatchedIterator.h:196
Aligns the peaks of two sorted spectra Method 1: Using a banded (width via 'tolerance' parameter) ali...
Definition: SpectrumAlignment.h:43
void getSpectrumAlignment(std::vector< std::pair< Size, Size > > &alignment, const SpectrumType1 &s1, const SpectrumType2 &s2) const
Definition: SpectrumAlignment.h:62
SpectrumAlignment & operator=(const SpectrumAlignment &source)
assignment operator
SpectrumAlignment(const SpectrumAlignment &source)
copy constructor
SpectrumAlignment()
default constructor
~SpectrumAlignment() override
destructor
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:97
static String & reverse(String &this_s)
Definition: StringUtilsSimple.h:330
Main OpenMS namespace.
Definition: openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:19