libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
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
30namespace pappso
31{
32
33
34//! Construct an uninitialized instance.
38
39
41 : MassSpectrumCombiner(decimal_places)
42{
43}
44
45
52
53
60
61
62//! Destruct the instance.
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// We get a trace as parameter that this combiner needs to combine to the
83// map_trace we also get as parameter. This combiner is for a mass spectrum, and
84// thus it needs to have been configured to hold a vector of m/z values that are
85// bins into which DataPoints from trace are combined into map_trace. It is the
86// bins that drive the the traversal of trace. For each bin, check if there is
87// (or are, if bins are big) data points in the trace that fit in it.
88
90MassSpectrumPlusCombiner::combine(MapTrace &map_trace, const Trace &trace) const
91{
92 // qDebug();
93
94 // If there is no single point in the trace we need to combine into the
95 // map_trace, then there is nothing to do.
96 if(!trace.size())
97 {
98 // qDebug() << "Thread:" << QThread::currentThreadId()
99 //<< "Returning right away because trace is empty.";
100 return map_trace;
101 }
102
103 // We will need to only use these iterator variables if we do not want to
104 // loose consistency.
105
106 using TraceIter = std::vector<DataPoint>::const_iterator;
107 TraceIter trace_iter_begin = trace.begin();
108 TraceIter trace_iter = trace_iter_begin;
109 TraceIter trace_iter_end = trace.end();
110
111 // The destination map trace will be filled-in with the result of the
112 // combination.
113
114 // Sanity check:
115 if(!m_bins.size())
116 throw(ExceptionNotPossible("The bin vector cannot be empty."));
117
118 using BinIter = std::vector<pappso_double>::const_iterator;
119
120 BinIter bin_iter = m_bins.begin();
121 BinIter bin_end_iter = m_bins.end();
122
123 // qDebug() << "initial bins iter at a distance of:"
124 //<< std::distance(m_bins.begin(), bin_iter)
125 //<< "bins distance:" << std::distance(m_bins.begin(), m_bins.end())
126 //<< "bins size:" << m_bins.size() << "first bin:" << m_bins.front()
127 //<< "last bin:" << m_bins.back();
128
129 // Iterate in the vector of bins and for each bin check if there are matching
130 // data points in the trace.
131
132 pappso_double current_bin_mz = 0;
133
134 pappso_double trace_x = 0;
135 pappso_double trace_y = 0;
136
137 // Lower bound returns an iterator pointing to the first element in the
138 // range [first, last) that is not less than (i.e. greater or equal to)
139 // value, or last if no such element is found.
140
141 auto bin_iter_for_mz = lower_bound(bin_iter, bin_end_iter, trace_iter->x);
142
143 if(bin_iter_for_mz != bin_end_iter)
144 {
145 if(bin_iter_for_mz != m_bins.begin())
146 bin_iter = --bin_iter_for_mz;
147 }
148 else
149 throw(ExceptionNotPossible("The bin vector must match the mz value."));
150
151 while(bin_iter != bin_end_iter)
152 {
153 current_bin_mz = *bin_iter;
154
155 // std::cout << std::setprecision(10) << "current bin: " << current_bin_mz
156 //<< std::endl;
157
158 // qDebug() << "Current bin:"
159 //<< QString("%1").arg(current_bin_mz, 0, 'f', 15)
160 //<< "at a distance of:"
161 //<< std::distance(m_bins.begin(), bin_iter);
162
163 // For the current bin, we start by instantiating a new DataPoint. By
164 // essence, each bin will have at most one corresponding DataPoint.
165
166 DataPoint new_data_point;
167
168 // Do not set the y value to 0 so that we can actually test if the
169 // data point is valid later on (try not to push back y=0 data
170 // points).
171 new_data_point.x = current_bin_mz;
172
173 // Now perform a loop over the data points in the mass spectrum.
174
175 // qDebug() << "Start looping in the trace to check if its DataPoint.x "
176 //"value matches that of the current bin."
177 //<< "trace_iter:" << trace_iter->toString()
178 //<< "data point distance:"
179 //<< std::distance(trace_iter_begin, trace_iter);
180
181 while(trace_iter != trace_iter_end)
182 {
183 // std::cout << std::setprecision(10)
184 //<< "trace data point being iterated into: " << trace_iter->x
185 //<< " - " << trace_iter->y << std::endl;
186
187 bool trace_matched = false;
188
189 // If trace is not to the end and the y value is not 0
190 // apply the shift, perform the rounding and check if the obtained
191 // x value is in the current bin, that is if it is less or equal
192 // to the current bin.
193
194 // qDebug() << "Thread:" << QThread::currentThreadId();
195 // qDebug() << "trace_iter:" << trace_iter->toString()
196 //<< "data point distance:"
197 //<< std::distance(trace_iter_begin, trace_iter);
198
199 // if(!Utils::almostEqual(trace_iter->y, 0.0, 10))
200
201 trace_x = trace_iter->x;
202 trace_y = trace_iter->y;
203
204 if(trace_y)
205 {
206
207 // FIXME: unbelievable behaviour: when building in release mode
208 // this code, under i386 (but not x86_64), this code fails if the
209 // following%S statement is missing. There are a variety of
210 // variants that work, the common denominator being that trace_x
211 // has to be output along with a string.
212
213 // std::cout << std::setprecision(10)
214 //<< "iterated trace point: " << trace_x << " - "
215 //<< trace_y << std::endl;
216
217 // For this reason I had to deactivate the related tests
218 // for i386 in tests/test_massspectrumcombiner.cpp
219
220 // trace_x is the m/z value that we need to combine,
221 // so make sure we check if there is a mz shift to apply.
222
223 // if(m_mzIntegrationParams.m_applyMzShift)
224 // trace_x += m_mzIntegrationParams.m_mzShift;
225
226 // Now apply the rounding (if any).
227 if(m_decimalPlaces != -1)
228 trace_x = Utils::roundToDecimals(trace_x, m_decimalPlaces);
229
230 if(trace_x <= current_bin_mz)
231 {
232
233 // std::cout
234 //<< std::setprecision(10)
235 //<< "Matched, because trace point.x is <= the current "
236 //"bin; increment trace_iter"
237 //<< std::endl;
238
239 new_data_point.y += trace_y;
240
241 // std::cout << std::setprecision(10)
242 //<< "After incrementation of the intensity: "
243 //<< new_data_point.y << std::endl;
244
245 // Let's record that we matched.
246 trace_matched = true;
247
248 // Because we matched, we can step-up with the
249 // iterator.
250 ++trace_iter;
251 }
252 // else
253 //{
254 // We did have a non-0 y value, but that did not
255 // match. So we do not step-up with the iterator.
256 //}
257 }
258 // End of
259 // if(trace_iter->y)
260 else
261 {
262 // std::cout << std::setprecision(10)
263 //<< "itered trace point intensity is 0: " << trace_x
264 //<< " - " << trace_y << std::endl;
265
266
267 // We iterated into a y=0 data point, so just skip it. Let
268 // the below code think that we have matched the point and
269 // iterate one step up.
270
271 // qDebug() << "The y value is almost equal to 0, increment the "
272 //"trace iter but do nothing else.";
273
274 trace_matched = true;
275 ++trace_iter;
276 }
277 // At this point, check if none of them matched.
278
279 if(!trace_matched)
280 {
281 // None of the first and trace mass spectra data
282 // points were found to match the current bin. All we
283 // have to do is go to the next bin. We break and the
284 // bin vector iterator will be incremented.
285
286 // However, if we had a valid new data point, that
287 // data point needs to be pushed back in the new mass
288 // spectrum.
289
290 if(new_data_point.isValid())
291 {
292
293 // We need to check if that bin value is present already in
294 // the map_trace object passed as parameter.
295
296 std::pair<std::map<pappso_double, pappso_double>::iterator,
297 bool>
298 result =
299 map_trace.insert(std::pair<pappso_double, pappso_double>(
300 new_data_point.x, new_data_point.y));
301
302 if(!result.second)
303 {
304 // The key already existed! The item was not inserted. We
305 // need to update the value.
306
307 result.first->second += new_data_point.y;
308
309 // qDebug()
310 //<< "Incremented the data point in the map trace.";
311 }
312 // else
313 //{
314 // qDebug()
315 //<< "Inserted a new data point into the map trace.";
316 //}
317 }
318
319 // We need to break this loop! That will increment the
320 // bin iterator.
321
322 break;
323 }
324 }
325 // End of
326 // while(trace_iter != trace_iter_end)
327
328 // Each time we get here, that means that we have consumed all
329 // the mass spectra data points that matched the current bin.
330 // So go to the next bin.
331
332 if(trace_iter == trace_iter_end)
333 {
334
335 // Make sure a last check is done in case one data point was
336 // cooking...
337
338 if(new_data_point.isValid())
339 {
340
341 std::pair<std::map<pappso_double, pappso_double>::iterator, bool>
342 result =
343 map_trace.insert(std::pair<pappso_double, pappso_double>(
344 new_data_point.x, new_data_point.y));
345
346 if(!result.second)
347 {
348 result.first->second += new_data_point.y;
349 }
350 }
351
352 // Now truly exit the loops...
353 break;
354 }
355
356 ++bin_iter;
357 }
358 // End of
359 // while(bin_iter != bin_end_iter)
360
361 // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
362 //<< "The combination result mass spectrum being returned has TIC:"
363 //<< new_trace.sumY();
364
365 return map_trace;
366}
367
368
369MapTrace &
371 const MapTrace &map_trace_in) const
372{
373 // qDebug();
374
375 if(!map_trace_in.size())
376 {
377 // qDebug() << "Thread:" << QThread::currentThreadId()
378 //<< "Returning right away because map_trace_in is empty.";
379 return map_trace_out;
380 }
381
382 Trace trace(map_trace_in);
383
384 return combine(map_trace_out, trace);
385}
386
387
388} // namespace pappso
int m_decimalPlaces
Number of decimals to use for the keys (x values)
std::vector< pappso_double > m_bins
MassSpectrumPlusCombiner()
Construct an uninitialized instance.
virtual ~MassSpectrumPlusCombiner()
Destruct the instance.
virtual MapTrace & combine(MapTrace &map_trace, const Trace &trace) const override
MassSpectrumPlusCombiner & operator=(const MassSpectrumPlusCombiner &other)
A simple container of DataPoint instances.
Definition trace.h:148
static pappso_double roundToDecimals(pappso_double value, int decimal_places)
Definition utils.cpp:142
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
double pappso_double
A type definition for doubles.
Definition types.h:50
std::shared_ptr< const MassSpectrumPlusCombiner > MassSpectrumPlusCombinerCstSPtr
pappso_double x
Definition datapoint.h:23
bool isValid() const
pappso_double y
Definition datapoint.h:24