libpappsomspp
Library for mass spectrometry
massspectrumminuscombiner.cpp
Go to the documentation of this file.
1 /////////////////////// StdLib includes
2 #include <numeric>
3 #include <limits>
4 #include <vector>
5 #include <map>
6 #include <cmath>
7 #include <iostream>
8 #include <iomanip>
9 
10 
11 /////////////////////// Qt includes
12 #include <QDebug>
13 #include <QFile>
14 #include <QThread>
15 #if 0
16 // For debugging purposes.
17 #include <QFile>
18 #endif
19 
20 
21 /////////////////////// Local includes
23 #include "../../types.h"
24 #include "../../utils.h"
25 #include "../../pappsoexception.h"
26 #include "../../exception/exceptionoutofrange.h"
27 #include "../../exception/exceptionnotpossible.h"
28 
29 
30 namespace pappso
31 {
32 
33 
34 //! Construct an uninitialized instance.
36 {
37 }
38 
39 
41  : MassSpectrumCombiner(decimal_places)
42 {
43 }
44 
45 
47  const MassSpectrumMinusCombiner &other)
48  : MassSpectrumCombiner(other)
49 
50 {
51 }
52 
53 
56  : MassSpectrumCombiner(other)
57 
58 {
59 }
60 
61 
62 //! Destruct the instance.
64 {
65 }
66 
67 
70 {
71  if(this == &other)
72  return *this;
73 
75 
76  m_bins.assign(other.m_bins.begin(), other.m_bins.end());
77 
78  return *this;
79 }
80 
81 
82 MapTrace &
84  const Trace &trace) const
85 {
86  // qDebug();
87 
88  if(!trace.size())
89  {
90  // qDebug() << "Thread:" << QThread::currentThreadId()
91  //<< "Returning right away because trace is empty.";
92  return map_trace;
93  }
94 
95  // We will need to only use these iterator variables if we do not want to
96  // loose consistency.
97 
98  using TraceIter = std::vector<DataPoint>::const_iterator;
99  TraceIter trace_iter_begin = trace.begin();
100  TraceIter trace_iter = trace_iter_begin;
101  TraceIter trace_iter_end = trace.end();
102 
103  // The destination map trace will be filled-in with the result of the
104  // combination.
105 
106  // Sanity check:
107  if(!m_bins.size())
108  throw(ExceptionNotPossible("The bin vector cannot be empty."));
109 
110  using BinIter = std::vector<pappso_double>::const_iterator;
111 
112  BinIter bin_iter = m_bins.begin();
113  BinIter bin_end_iter = m_bins.end();
114 
115  // qDebug() << "initial bins iter at a distance of:"
116  //<< std::distance(m_bins.begin(), bin_iter)
117  //<< "bins distance:" << std::distance(m_bins.begin(), m_bins.end())
118  //<< "bins size:" << m_bins.size() << "first bin:" << m_bins.front()
119  //<< "last bin:" << m_bins.back();
120 
121  // Iterate in the vector of bins and for each bin check if there are matching
122  // data points in the trace.
123 
124  pappso_double current_bin_mz_value = 0;
125 
126  pappso_double trace_x = 0;
127  pappso_double trace_y = 0;
128 
129  // Lower bound returns an iterator pointing to the first element in the
130  // range [first, last) that is not less than (i.e. greater or equal to)
131  // value, or last if no such element is found.
132 
133  // Get the first bin that would match the very first data point in the trace
134  // that we need to combine into the map_trace.
135 
136  auto bin_iter_for_mz = lower_bound(bin_iter, bin_end_iter, trace_iter->x);
137 
138  if(bin_iter_for_mz != bin_end_iter)
139  {
140  // If we found a bin and that the bin is not the first bin of the bin
141  // vector, then go back one bin to be sure we do not miss data points.
142 
143  if(bin_iter_for_mz != m_bins.begin())
144  bin_iter = --bin_iter_for_mz;
145  }
146  else
147  throw(ExceptionNotPossible("The bin vector must match the mz value."));
148 
149  // Now iterate in the bin vector starting from the std::prev(found bin).
150 
151  while(bin_iter != bin_end_iter)
152  {
153  current_bin_mz_value = *bin_iter;
154 
155  // qDebug() << "Iterating in new bin:"
156  //<< QString("%1").arg(current_bin_mz_value, 0, 'f', 15)
157  //<< "at a distance from the first bin of:"
158  //<< std::distance(m_bins.begin(), bin_iter);
159 
160  // For the current bin, we start by instantiating a new DataPoint. By
161  // essence, each bin will have at most one corresponding DataPoint.
162 
163  DataPoint bin_data_point;
164 
165  // Do not set the y value to 0 so that we can actually test if the
166  // data point is valid later on (try not to push back y=0 data
167  // points).
168 
169  bin_data_point.x = current_bin_mz_value;
170 
171  // qDebug() << "Seed bin_data_point.x with current_bin_mz_value value:"
172  //<< QString("%1").arg(bin_data_point.x, 0, 'f', 15);
173 
174  // Now perform a loop over the data points in the mass spectrum.
175 
176  while(trace_iter != trace_iter_end)
177  {
178  bool trace_matched = false;
179 
180  // If we are not at the end of trace and if the y value of the
181  // currently iterated trace datapoint is not 0, perform the rounding
182  // and check if the obtained x value is in the current bin, that is if
183  // it is less or equal to the current bin.
184 
185  // qDebug() << "Thread:" << QThread::currentThreadId();
186 
187  // qDebug() << "Iterating in new trace datapoint:"
188  //<< trace_iter->toString(15) << "at distance:"
189  //<< std::distance(trace_iter_begin, trace_iter);
190 
191  if(trace_iter->y)
192  {
193  // qDebug() << "The y value of trace iterated value is not 0.";
194 
195  // trace_x is the m/z value that we need to combine,
196 
197  trace_x = trace_iter->x;
198  trace_y = trace_iter->y;
199 
200  // Now apply the rounding (if any).
201  if(m_decimalPlaces != -1)
202  trace_x = Utils::roundToDecimals(trace_x, m_decimalPlaces);
203 
204  if(trace_x <= current_bin_mz_value)
205  {
206  // qDebug()
207  //<< "trace_x <= current_bin_mz_value: the match happened.";
208 
209  // qDebug() << "Going to compute" << bin_data_point.y << "-"
210  //<< trace_y << "with current_bin_mz_value:"
211  //<< QString("%1").arg(
212  // current_bin_mz_value, 0, 'f', 15);
213 
214  // This is where the bin_data_point actually gets decremented
215  // by the y value of the currently iterated trace data point.
216  // So bin_data_point has already the MINUS'ed value in it.
217  // Relate this code line with the result.first->second +=
218  // bin_data_point.y; code line elsewhere in this function.
219  bin_data_point.y -= trace_y;
220 
221  // qDebug() << "New bin_data_point.y value:" <<
222  // bin_data_point.y;
223 
224  // Let's record that we matched.
225  trace_matched = true;
226 
227  // Because we matched, we can step-up with the
228  // iterator.
229 
230  // qDebug() << "The values matched, increment the trace "
231  //"itereator but not the bin iterator, as maybe "
232  //"the next trace datapoint also matches the bin.";
233 
234  ++trace_iter;
235  }
236  // else
237  //{
238  // We did have a non-0 y value, but that did not
239  // match. So we do not step-up with the iterator because that
240  // trace data point will fit into another bin not yet iterated
241  // into.
242  //}
243  }
244  // End of
245  // if(trace_iter->y)
246  else
247  {
248  // We iterated into a y=0 data point, so just skip it. Let the
249  // below code think that we have matched the point and iterate one
250  // step up.
251 
252  // qDebug() << "The y value of the currently iterated trace data "
253  //"point is almost equal to 0, increment the "
254  //"trace iter but do nothing else.";
255 
256  trace_matched = true;
257  ++trace_iter;
258  }
259 
260  // At this point, check if the trace data point matched the currently
261  // iterated bin.
262 
263  if(!trace_matched)
264  {
265 
266  // qDebug() << "The currently iterated trace datapoint did not "
267  //"match the current bin.";
268 
269  // The trace data point did not match the currently iterated bin.
270  // All we have to do is go to the next bin. We break and the bin
271  // vector iterator will be incremented.
272 
273  // However, if we had a valid new data point, that
274  // data point needs to be pushed back in the new mass
275  // spectrum.
276 
277  if(bin_data_point.isValid())
278  {
279 
280  // qDebug() << "But we had a valid bin data point cooking, so
281  // " "we have to take that into account:"
282  //<< bin_data_point.toString(15);
283 
284  // We need to check if that bin value is present already in
285  // the map_trace object passed as parameter.
286 
287  std::pair<std::map<pappso_double, pappso_double>::iterator,
288  bool>
289  result =
290  map_trace.insert(std::pair<pappso_double, pappso_double>(
291  bin_data_point.x, -bin_data_point.y));
292 
293  if(!result.second)
294  {
295  // qDebug()
296  //<< "The x value (mz) was already present in the "
297  //"map trace, simply update its y value (intensity).";
298 
299  // The key already existed! The item was not inserted. We
300  // need to update the y value.
301 
302  // qDebug() << "Going to compute" << result.first->second
303  //<< "+" << bin_data_point.y;
304 
305  // This might be counterintuitive: one would expect the
306  // use of the '-=' operator because we are
307  // MINUS-combining. But no! bin_data_point only get y
308  // values already minus'ed. So we need to add that
309  // minus'ed value, otherwise we add the value (- -value is
310  // the same as +value).
311  result.first->second += bin_data_point.y;
312 
313  // qDebug() << "New y value for the item in the map trace
314  // "
315  //"(result.first->second value):"
316  //<< result.first->second;
317  }
318  // else
319  //{
320  // qDebug()
321  //<< "The data point has a mz value that was not present "
322  //"in the map trace."
323  //<< " Inserted a new data point into the map trace:"
324  //<< bin_data_point.x << "," << bin_data_point.y;
325  //}
326  }
327 
328  // We need to break this loop! That will increment the
329  // bin iterator.
330 
331  break;
332  }
333  }
334  // End of
335  // while(trace_iter != trace_iter_end)
336 
337  // Each time we get here, that means that we have consumed all
338  // the mass spectra data points that matched the current bin.
339  // So go to the next bin.
340 
341  if(trace_iter == trace_iter_end)
342  {
343 
344  // Make sure a last check is done in case one data point was
345  // cooking...
346 
347  if(bin_data_point.isValid())
348  {
349 
350  std::pair<std::map<pappso_double, pappso_double>::iterator, bool>
351  result =
352  map_trace.insert(std::pair<pappso_double, pappso_double>(
353  bin_data_point.x, -bin_data_point.y));
354 
355  if(!result.second)
356  {
357  // qDebug()
358  //<< "The x value (mz) was already present in the "
359  //"map trace, simply update its y value (intensity).";
360 
361  // The key already existed! The item was not inserted. We
362  // need to update the y value.
363 
364  // qDebug() << "Going to compute" << result.first->second <<
365  // "+"
366  //<< bin_data_point.y;
367 
368 
369  // This might be counterintuitive: one would expect the
370  // use of the '-=' operator because we are
371  // MINUS-combining. But no! bin_data_point only get y
372  // values already minus'ed. So we need to add that
373  // minus'ed value, otherwise we add the value (- -value is
374  // the same as +value).
375  result.first->second += bin_data_point.y;
376 
377  // qDebug() << "New y value for the item in the map trace "
378  //"(result.first->second value):"
379  //<< result.first->second;
380  }
381  // else
382  //{
383  // qDebug()
384  //<< "The data point has a mz value that was not present "
385  //"in the map trace."
386  //<< " Inserted a new data point into the map trace:"
387  //<< bin_data_point.x << "," << bin_data_point.y;
388  //}
389  }
390 
391  // Now truly exit the loops...
392  break;
393  }
394 
395  ++bin_iter;
396  }
397  // End of
398  // while(bin_iter != bin_end_iter)
399 
400  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
401  //<< "The combination result mass spectrum being returned has TIC:"
402  //<< new_trace.sumY();
403 
404  return map_trace;
405 }
406 
407 
408 MapTrace &
410  const MapTrace &map_trace_in) const
411 {
412  // qDebug();
413 
414  if(!map_trace_in.size())
415  {
416  // qDebug() << "Thread:" << QThread::currentThreadId()
417  //<< "Returning right away because map_trace_in is empty.";
418  return map_trace_out;
419  }
420 
421  Trace trace(map_trace_in);
422 
423  return combine(map_trace_out, trace);
424 }
425 
426 } // namespace pappso
int m_decimalPlaces
Number of decimals to use for the keys (x values)
std::vector< pappso_double > m_bins
virtual MapTrace & combine(MapTrace &map_trace, const Trace &trace) const override
MassSpectrumMinusCombiner()
Construct an uninitialized instance.
virtual ~MassSpectrumMinusCombiner()
Destruct the instance.
MassSpectrumMinusCombiner & operator=(const MassSpectrumMinusCombiner &other)
A simple container of DataPoint instances.
Definition: trace.h:148
static pappso_double roundToDecimals(pappso_double value, int decimal_places)
Definition: utils.cpp:120
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
std::shared_ptr< const MassSpectrumMinusCombiner > MassSpectrumMinusCombinerCstSPtr
double pappso_double
A type definition for doubles.
Definition: types.h:49
pappso_double x
Definition: datapoint.h:23
bool isValid() const
Definition: datapoint.cpp:132
pappso_double y
Definition: datapoint.h:24