libpappsomspp
Library for mass spectrometry
massspectrumpluscombiner.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 MassSpectrumPlusCombiner &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 = 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  auto bin_iter_for_mz = lower_bound(bin_iter, bin_end_iter, trace_iter->x);
122 
123  if(bin_iter_for_mz != bin_end_iter)
124  {
125  if(bin_iter_for_mz != m_bins.begin())
126  bin_iter = --bin_iter_for_mz;
127  }
128  else
129  throw(ExceptionNotPossible("The bin vector must match the mz value."));
130 
131  while(bin_iter != bin_end_iter)
132  {
133  current_bin = *bin_iter;
134 
135  // qDebug() << "Current bin:" << QString("%1").arg(current_bin, 0, 'f',
136  // 15)
137  //<< "at a distance of:"
138  //<< std::distance(m_bins.begin(), bin_iter);
139 
140  // For the current bin, we start by instantiating a new DataPoint. By
141  // essence, each bin will have at most one corresponding DataPoint.
142 
143  DataPoint new_data_point;
144 
145  // Do not set the y value to 0 so that we can actually test if the
146  // data point is valid later on (try not to push back y=0 data
147  // points).
148  new_data_point.x = current_bin;
149 
150  // Now perform a loop over the data points in the mass spectrum.
151 
152  // qDebug() << "trace_iter:" << trace_iter->toString()
153  //<< "data point distance:"
154  //<< std::distance(trace_iter_begin, trace_iter);
155 
156  while(trace_iter != trace_iter_end)
157  {
158 
159  bool trace_matched = false;
160 
161  // If trace is not to the end and the y value is not 0
162  // apply the shift, perform the rounding and check if the obtained
163  // x value is in the current bin, that is if it is less or equal
164  // to the current bin.
165 
166  // qDebug() << "Thread:" << QThread::currentThreadId();
167  // qDebug() << "trace_iter:" << trace_iter->toString()
168  //<< "data point distance:"
169  //<< std::distance(trace_iter_begin, trace_iter);
170 
171  // if(!Utils::almostEqual(trace_iter->y, 0.0, 10))
172  if(trace_iter->y)
173  {
174  trace_x = trace_iter->x;
175  trace_y = trace_iter->y;
176 
177  // trace_x is the m/z value that we need to combine,
178  // so make sure we check if there is a mz shift to apply.
179 
180  // if(m_mzIntegrationParams.m_applyMzShift)
181  // trace_x += m_mzIntegrationParams.m_mzShift;
182 
183  // Now apply the rounding (if any).
184  if(m_decimalPlaces != -1)
185  trace_x = Utils::roundToDecimals(trace_x, m_decimalPlaces);
186 
187  if(trace_x <= current_bin)
188  {
189 
190  // qDebug() << "Matched, increment trace_iter";
191  new_data_point.y += trace_y;
192 
193  // Let's record that we matched.
194  trace_matched = true;
195 
196  // Because we matched, we can step-up with the
197  // iterator.
198  ++trace_iter;
199  }
200  // else
201  //{
202  // We did have a non-0 y value, but that did not
203  // match. So we do not step-up with the iterator.
204  //}
205  }
206  // End of
207  // if(trace_iter->y)
208  else
209  {
210  // We iterated into a y=0 data point, so just skip it. Let
211  // the below code think that we have matched the point and
212  // iterate one step up.
213 
214  // qDebug() << "The y value is almost equal to 0, increment the "
215  //"trace iter but do nothing else.";
216 
217  trace_matched = true;
218  ++trace_iter;
219  }
220  // At this point, check if none of them matched.
221 
222  if(!trace_matched)
223  {
224  // None of the first and trace mass spectra data
225  // points were found to match the current bin. All we
226  // have to do is go to the next bin. We break and the
227  // bin vector iterator will be incremented.
228 
229  // However, if we had a valid new data point, that
230  // data point needs to be pushed back in the new mass
231  // spectrum.
232 
233  if(new_data_point.isValid())
234  {
235 
236  // We need to check if that bin value is present already in
237  // the map_trace object passed as parameter.
238 
239  std::pair<std::map<pappso_double, pappso_double>::iterator,
240  bool>
241  result =
242  map_trace.insert(std::pair<pappso_double, pappso_double>(
243  new_data_point.x, new_data_point.y));
244 
245  if(!result.second)
246  {
247  // The key already existed! The item was not inserted. We
248  // need to update the value.
249 
250  result.first->second += new_data_point.y;
251 
252  // qDebug() << "Incremented the data point in the map
253  // trace.";
254  }
255  // else
256  //{
257  // qDebug()
258  //<< "Inserted a new data point into the map trace.";
259  //}
260  }
261 
262  // We need to break this loop! That will increment the
263  // bin iterator.
264 
265  break;
266  }
267  }
268  // End of
269  // while(1)
270 
271  // Each time we get here, that means that we have consumed all
272  // the mass spectra data points that matched the current bin.
273  // So go to the next bin.
274 
275  if(trace_iter == trace_iter_end)
276  {
277 
278  // Make sure a last check is done in case one data point was
279  // cooking...
280 
281  if(new_data_point.isValid())
282  {
283 
284  std::pair<std::map<pappso_double, pappso_double>::iterator, bool>
285  result =
286  map_trace.insert(std::pair<pappso_double, pappso_double>(
287  new_data_point.x, new_data_point.y));
288 
289  if(!result.second)
290  {
291  result.first->second += new_data_point.y;
292  }
293  }
294 
295  // Now truly exit the loops...
296  break;
297  }
298 
299  ++bin_iter;
300  }
301  // End of
302  // while(bin_iter != bin_end_iter)
303 
304  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
305  //<< "The combination result mass spectrum being returned has TIC:"
306  //<< new_trace.sumY();
307 
308  return map_trace;
309 }
310 
311 
312 } // namespace pappso
313 
pappso::MassSpectrumPlusCombinerCstSPtr
std::shared_ptr< const MassSpectrumPlusCombiner > MassSpectrumPlusCombinerCstSPtr
Definition: massspectrumpluscombiner.h:16
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition: types.h:48
pappso::MassSpectrumPlusCombiner::combineNoFilteringStep
virtual MapTrace & combineNoFilteringStep(MapTrace &map_trace, const Trace &trace) const
Definition: massspectrumpluscombiner.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::DataPoint
Definition: datapoint.h:21
pappso::MapTrace
Definition: maptrace.h:33
pappso::ExceptionNotPossible
Definition: exceptionnotpossible.h:29
pappso::MassSpectrumPlusCombiner
Definition: massspectrumpluscombiner.h:25
pappso::MassSpectrumPlusCombiner::~MassSpectrumPlusCombiner
virtual ~MassSpectrumPlusCombiner()
Destruct the instance.
Definition: massspectrumpluscombiner.cpp:65
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
massspectrumpluscombiner.h
pappso::Utils::roundToDecimals
static pappso_double roundToDecimals(pappso_double value, int decimal_places)
Definition: utils.cpp:104
pappso::MassSpectrumPlusCombiner::MassSpectrumPlusCombiner
MassSpectrumPlusCombiner()
Construct an uninitialized instance.
Definition: massspectrumpluscombiner.cpp:35
pappso::MassSpectrumCombiner::m_bins
std::vector< pappso_double > m_bins
Definition: massspectrumcombiner.h:58