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