libpappsomspp
Library for mass spectrometry
massspectrumminuscombiner.cpp
Go to the documentation of this file.
1 #include <numeric>
2 #include <limits>
3 #include <vector>
4 #include <map>
5 #include <cmath>
6 
7 #include <QDebug>
8 #include <QFile>
9 
10 #if 0
11 // For debugging purposes.
12 #include <QFile>
13 #endif
14 
16 #include "../../types.h"
17 #include "../../utils.h"
18 #include "../../pappsoexception.h"
19 #include "../../exception/exceptionoutofrange.h"
20 #include "../../exception/exceptionnotpossible.h"
21 #include "traceminuscombiner.h"
22 
23 
24 namespace pappso
25 {
26 
27 
28 //! Construct an uninitialized instance.
30 {
31 }
32 
33 
35  : MassSpectrumCombiner(decimal_places)
36 {
37 }
38 
39 
42  : MassSpectrumCombiner(other)
43 
44 {
45  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()";
46 }
47 
48 
50  const MassSpectrumMinusCombiner &other)
51  : MassSpectrumCombiner(other)
52 
53 {
54  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()";
55 }
56 
57 
58 //! Destruct the instance.
60 {
61 }
62 
63 
64 MapTrace &
66  const Trace &trace) const
67 {
68  if(!trace.size())
69  {
70  // qDebug() << "Thread:" << QThread::currentThreadId()
71  //<< "Returning right away because trace is empty.";
72  return map_trace;
73  }
74 
75  // We will need to only use these iterator variables if we do not want to
76  // loose consistency.
77 
78  using TraceIter = std::vector<DataPoint>::const_iterator;
79  TraceIter trace_iter_begin = trace.begin();
80  TraceIter trace_iter = trace_iter_begin;
81  TraceIter trace_iter_end = trace.end();
82 
83  // The destination map trace will be filled-in with the result of the
84  // combination.
85 
86  // Sanity check:
87  if(!m_bins.size())
88  throw(ExceptionNotPossible("The bin vector cannot be empty."));
89 
90  using BinIter = std::vector<pappso_double>::const_iterator;
91 
92  BinIter bin_iter = m_bins.begin();
93  BinIter bin_end_iter = m_bins.end();
94 
95  // qDebug() << "initial bins iter at a distance of:"
96  //<< std::distance(m_bins.begin(), bin_iter)
97  //<< "bins distance:" << std::distance(m_bins.begin(), m_bins.end())
98  //<< "bins size:" << m_bins.size() << "first bin:" << m_bins.front()
99  //<< "last bin:" << m_bins.back();
100 
101  // Iterate in the vector of bins and for each bin check if there are matching
102  // data points in the trace.
103 
104  pappso_double current_bin = 0;
105 
106  pappso_double trace_x = 0;
107  pappso_double trace_y = 0;
108 
109  // Lower bound returns an iterator pointing to the first element in the
110  // range [first, last) that is not less than (i.e. greater or equal to)
111  // value, or last if no such element is found.
112 
113  auto bin_iter_for_mz = lower_bound(bin_iter, bin_end_iter, trace_iter->x);
114 
115  if(bin_iter_for_mz != bin_end_iter)
116  {
117  if(bin_iter_for_mz != m_bins.begin())
118  bin_iter = --bin_iter_for_mz;
119  }
120  else
121  throw(ExceptionNotPossible("The bin vector must match the mz value."));
122 
123  while(bin_iter != bin_end_iter)
124  {
125  current_bin = *bin_iter;
126 
127  // qDebug() << "Current bin:" << QString("%1").arg(current_bin, 0, 'f',
128  // 15)
129  //<< "at a distance of:"
130  //<< std::distance(m_bins.begin(), bin_iter);
131 
132  // For the current bin, we start by instantiating a new DataPoint. By
133  // essence, each bin will have at most one corresponding DataPoint.
134 
135  DataPoint new_data_point;
136 
137  // Do not set the y value to 0 so that we can actually test if the
138  // data point is valid later on (try not to push back y=0 data
139  // points).
140  new_data_point.x = current_bin;
141 
142  // Now perform a loop over the data points in the mass spectrum.
143 
144  // qDebug() << "trace_iter:" << trace_iter->toString()
145  //<< "data point distance:"
146  //<< std::distance(trace_iter_begin, trace_iter);
147 
148  while(trace_iter != trace_iter_end)
149  {
150 
151  bool trace_matched = false;
152 
153  // If trace is not to the end and the y value is not 0
154  // apply the shift, perform the rounding and check if the obtained
155  // x value is in the current bin, that is if it is less or equal
156  // to the current bin.
157 
158  // qDebug() << "Thread:" << QThread::currentThreadId();
159  // qDebug() << "trace_iter:" << trace_iter->toString()
160  //<< "data point distance:"
161  //<< std::distance(trace_iter_begin, trace_iter);
162 
163  // if(!Utils::almostEqual(trace_iter->y, 0.0, 10))
164  if(trace_iter->y)
165  {
166  trace_x = trace_iter->x;
167  trace_y = trace_iter->y;
168 
169  // trace_x is the m/z value that we need to combine,
170  // so make sure we check if there is a mz shift to apply.
171 
172  // if(m_mzIntegrationParams.m_applyMzShift)
173  // trace_x += m_mzIntegrationParams.m_mzShift;
174 
175  // Now apply the rounding (if any).
176  if(m_decimalPlaces != -1)
177  trace_x = Utils::roundToDecimals(trace_x, m_decimalPlaces);
178 
179  if(trace_x <= current_bin)
180  {
181 
182  // qDebug() << "Matched, increment trace_iter";
183  new_data_point.y -= trace_y;
184 
185  // Let's record that we matched.
186  trace_matched = true;
187 
188  // Because we matched, we can step-up with the
189  // iterator.
190  ++trace_iter;
191  }
192  // else
193  //{
194  // We did have a non-0 y value, but that did not
195  // match. So we do not step-up with the iterator.
196  //}
197  }
198  // End of
199  // if(trace_iter->y)
200  else
201  {
202  // We iterated into a y=0 data point, so just skip it. Let
203  // the below code think that we have matched the point and
204  // iterate one step up.
205 
206  // qDebug() << "The y value is almost equal to 0, increment the "
207  //"trace iter but do nothing else.";
208 
209  trace_matched = true;
210  ++trace_iter;
211  }
212  // At this point, check if none of them matched.
213 
214  if(!trace_matched)
215  {
216  // None of the first and trace mass spectra data
217  // points were found to match the current bin. All we
218  // have to do is go to the next bin. We break and the
219  // bin vector iterator will be incremented.
220 
221  // However, if we had a valid new data point, that
222  // data point needs to be pushed back in the new mass
223  // spectrum.
224 
225  if(new_data_point.isValid())
226  {
227 
228  // We need to check if that bin value is present already in
229  // the map_trace object passed as parameter.
230 
231  std::pair<std::map<pappso_double, pappso_double>::iterator,
232  bool>
233  result =
234  map_trace.insert(std::pair<pappso_double, pappso_double>(
235  new_data_point.x, -new_data_point.y));
236 
237  if(!result.second)
238  {
239  // The key already existed! The item was not inserted. We
240  // need to update the value.
241 
242  result.first->second -= new_data_point.y;
243 
244  // qDebug() << "Incremented the data point in the map
245  // trace.";
246  }
247  // else
248  //{
249  // qDebug()
250  //<< "Inserted a new data point into the map trace.";
251  //}
252  }
253 
254  // We need to break this loop! That will increment the
255  // bin iterator.
256 
257  break;
258  }
259  }
260  // End of
261  // while(1)
262 
263  // Each time we get here, that means that we have consumed all
264  // the mass spectra data points that matched the current bin.
265  // So go to the next bin.
266 
267  if(trace_iter == trace_iter_end)
268  {
269 
270  // Make sure a last check is done in case one data point was
271  // cooking...
272 
273  if(new_data_point.isValid())
274  {
275 
276  std::pair<std::map<pappso_double, pappso_double>::iterator, bool>
277  result =
278  map_trace.insert(std::pair<pappso_double, pappso_double>(
279  new_data_point.x, -new_data_point.y));
280 
281  if(!result.second)
282  {
283  result.first->second -= new_data_point.y;
284  }
285  }
286 
287  // Now truly exit the loops...
288  break;
289  }
290 
291  ++bin_iter;
292  }
293  // End of
294  // while(bin_iter != bin_end_iter)
295 
296  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
297  //<< "The combination result mass spectrum being returned has TIC:"
298  //<< new_trace.sumY();
299 
300  return map_trace;
301 }
302 
303 
304 } // namespace pappso
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition: types.h:69
pappso::MassSpectrumMinusCombiner::combineNoFilteringStep
virtual MapTrace & combineNoFilteringStep(MapTrace &map_trace, const Trace &trace) const
Definition: massspectrumminuscombiner.cpp:65
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
Definition: aa.cpp:39
traceminuscombiner.h
pappso::MassSpectrumMinusCombiner::~MassSpectrumMinusCombiner
virtual ~MassSpectrumMinusCombiner()
Destruct the instance.
Definition: massspectrumminuscombiner.cpp:59
pappso::DataPoint
Definition: datapoint.h:21
pappso::MapTrace
Definition: maptrace.h:33
pappso::ExceptionNotPossible
Definition: exceptionnotpossible.h:53
pappso::Trace
A simple container of DataPoint instances.
Definition: trace.h:132
pappso::DataPoint::isValid
bool isValid() const
Definition: datapoint.cpp:117
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:29
pappso::Utils::roundToDecimals
static pappso_double roundToDecimals(pappso_double value, int decimal_places)
Definition: utils.cpp:125
pappso::MassSpectrumMinusCombiner
Definition: massspectrumminuscombiner.h:26
pappso::MassSpectrumCombiner::m_bins
std::vector< pappso_double > m_bins
Definition: massspectrumcombiner.h:58
massspectrumminuscombiner.h