massspectrumpluscombiner.cpp

Go to the documentation of this file.

/////////////////////// StdLib includes
#include <vector>
#include <map>
#include <cmath>
 
 
/////////////////////// Qt includes
#include <QDebug>
#include <QFile>
#include <QThread>
#if 0
// For debugging purposes.
#include <QFile>
#endif
 
 
/////////////////////// Local includes
#include "massspectrumpluscombiner.h"
#include "pappsomspp/core/utils.h"
#include "pappsomspp/core/exception/exceptionnotpossible.h"
 
namespace pappso
{
 
 
//! Construct an uninitialized instance.
MassSpectrumPlusCombiner::MassSpectrumPlusCombiner(QObject *parent_p)
  : MassSpectrumCombiner(parent_p)
{
}
 
MassSpectrumPlusCombiner::MassSpectrumPlusCombiner(int decimal_places,
                                                   QObject *parent_p)
  : MassSpectrumCombiner(decimal_places, parent_p)
{
}
 
MassSpectrumPlusCombiner::MassSpectrumPlusCombiner(
  const MassSpectrumPlusCombiner &other, QObject *parent_p)
  : MassSpectrumCombiner(other, parent_p)
 
{
}
 
MassSpectrumPlusCombiner::MassSpectrumPlusCombiner(
  MassSpectrumPlusCombinerCstSPtr other, QObject *parent_p)
  : MassSpectrumCombiner(other, parent_p)
 
{
}
 
//! Destruct the instance.
MassSpectrumPlusCombiner::~MassSpectrumPlusCombiner()
{
}
 
MassSpectrumPlusCombiner &
MassSpectrumPlusCombiner::operator=(const MassSpectrumPlusCombiner &other)
{
  if(this == &other)
    return *this;
 
  MassSpectrumCombiner::m_decimalPlaces = other.m_decimalPlaces;
 
  m_bins.assign(other.m_bins.begin(), other.m_bins.end());
 
  return *this;
}
 
// We get a trace as parameter that this combiner needs to combine into the
// map_trace we also get as parameter. This combiner is for a mass spectrum, and
// thus it needs to have been configured to hold a vector of m/z values that are
// bins into which DataPoint_s from trace are combined into map_trace. It is the
// bins that drive the traversal of trace. For each bin, check if there is
// (or are, if bins are big) data points in the trace that fit in it.
 
// One important consideration is that the bins are formed before the call to
// this combine() function and thus the distance between one m/z value and the
// next is already determined. This means that we do not need to bother with the
// conversion of m/z values to new decimal_places-defined values.
 
// Another idea is that, upon the very first combination of the first spectrum
// of the the thousands of spectra that are going to potentially be combined
// into map_trace object, that map_trace object is *empty*. This means that we
// can set the bins in that map_trace object and then iterated directly in the
// destination map_trace and do the intensity incrementation calculation right
// into the .second member of the pair that is being iterated into in the
// map_trace object.
MapTrace &
MassSpectrumPlusCombiner::combine(MapTrace &map_trace, const Trace &trace) const
{
  // qDebug();
 
  // If there is no single point in the trace we need to combine into the
  // map_trace, then there is nothing to do.
  if(!trace.size())
    {
      // qDebug() << "Thread:" << QThread::currentThreadId()
      //<< "Returning right away because trace is empty.";
      return map_trace;
    }
 
  // The destination map trace will be filled-in with the result of the
  // combination.
 
  // Sanity check:
  if(!m_bins.size())
    throw(ExceptionNotPossible("The bin vector cannot be empty."));
 
  // qDebug() << qSetRealNumberPrecision(6) << "There are " << m_bins.size() <<
  // "bins running from"
  //          << m_bins.front() << "to" << m_bins.back();
 
  // Since the map_trace argument is used all along the combination of
  // potentially thousands of mass spectra, we need to be careful. But
  // clearly, if that map_trace object is empty, then, one way to accelerate the
  // whole process is to copy the bins in there and later assume that
  // the bin is already in it. So that we need not test each time if we have
  // inserted the new data point or if we simply need to increment its
  // intensity.
 
  if(!map_trace.size())
    {
      std::transform(m_bins.begin(),
                     m_bins.end(),
                     std::inserter(map_trace, map_trace.end()),
                     [](double key) {
                       return std::make_pair(key, 0.0);
                     });
    }
 
  // We will need to only use these iterator variables if we do not want to
  // loose consistency.
 
  using TraceIter            = std::vector<DataPoint>::const_iterator;
  TraceIter trace_iter_begin = trace.begin();
  TraceIter trace_iter       = trace_iter_begin;
  TraceIter trace_iter_end   = trace.end();
 
  // qDebug() << qSetRealNumberPrecision(6) << "The trace to combine has " <<
  // trace.size()
  //          << "data points, with first a m/z" << trace.front().x << "and last
  //          at m/z"
  //          << trace.back().x;
 
  // We will iterate in the destination map_trace that has been filled with
  // the bins initially computed for this combiner.
  using BinIter                  = std::map<double, double>::iterator;
  BinIter map_trace_bin_iter     = map_trace.begin();
  BinIter map_trace_bin_end_iter = map_trace.end();
 
  // We do not want to iterate in too many map_trace items before arriving
  // to the map_trace item that matches the m/z value of the first data
  // point in trace.
 
  // Get the first map_trace bin that would match the very first data point in
  // the trace that we need to combine into the map_trace. We want the first
  // map_trace element having key >= trace_iter->x.
  BinIter map_trace_bin_iter_for_mz = map_trace.upper_bound(trace_iter->x);
 
  // If it's not the beginning, the previous element is what we want
  if(map_trace_bin_iter_for_mz != map_trace.begin())
    {
      map_trace_bin_iter =
        --map_trace_bin_iter_for_mz; // This is the greatest map_trace key <
                                     // trace_iter->x
 
      // qDebug() << "Found map_trace key=" << map_trace_bin_iter_for_mz->first
      //          << ", value=" << map_trace_bin_iter_for_mz->second;
    }
  else
    {
      // qDebug() << "The first bin found to match the first data point in trace
      // is:"
      // << map_trace.begin()->first;
    }
 
  // qDebug() << "The first bin found to match the first m/z value in the trace
  // to combine ("
  //          << trace_iter->x
  //          << ") is at index:" << std::distance(map_trace.begin(),
  //          map_trace_bin_iter);
 
  while(map_trace_bin_iter != map_trace_bin_end_iter)
    {
      // qDebug() << qSetRealNumberPrecision(6)
      //          << "Now iterating in new bin with mz:" <<
      //          map_trace_bin_iter->first
      //          << "at a distance of the first bin:"
      //          << std::distance(map_trace.begin(), map_trace_bin_iter);
 
      // Now perform a loop over the data points in the mass spectrum.
 
      // qDebug() << "Start looping in the trace to check if its DataPoint.x "
      //"value matches that of the current bin."
      //<< "trace_iter:" << trace_iter->toString()
      //<< "at data point distance from first point in trace:"
 
      while(trace_iter != trace_iter_end)
        {
          // qDebug() << qSetRealNumberPrecision(6)
          //          << "trace data point being iterated into: " <<
          //          trace_iter->x << " - "
          //          << trace_iter->y;
 
          // If trace is not to the end and the y value is not 0, perform the
          // rounding and check if the obtained x value is in the current bin,
          // that is if it is less or equal to the current bin.
 
          // qDebug() << "Thread:" << QThread::currentThreadId();
          // qDebug() << "trace_iter:" << trace_iter->toString()
          //<< "data point distance:"
          //<< std::distance(trace_iter_begin, trace_iter);
 
 
          if(trace_iter->x <= map_trace_bin_iter->first)
            {
              // qDebug() << qSetRealNumberPrecision(6) << "Matched because
              // trace x" << trace_iter->x
              //          << "is <= current bin m/z" <<
              //          map_trace_bin_iter->first;
 
              // qDebug() << "Incrementing the y value by " << trace_iter->y;
 
              map_trace_bin_iter->second += trace_iter->y;
 
              // qDebug() << qSetRealNumberPrecision(6)
              //          << "After incrementation of the intensity: " <<
              //          map_trace_bin_iter->second;
 
              // Because we matched, we can step-up with the
              // trace iterator.  The idea is that, if the bins are large, then,
              // multiple data points from a given mass spec acquisition' trace
              // may fall inside that bin.
 
              // qDebug() << "Now increment the trace iterator.";
 
              ++trace_iter;
            }
          else
            {
              // The trace data point's x did not match the current map_trace
              // bin m/z value because the former was greater than the latter.
 
              // Because the DataPoint instances in the trace are sorted in
              // increasing x order, we know we won't have to add any data point
              // to this bin. So we break this loop, which will increment the
              // map_trace bin iterator.
 
              break;
            }
        }
      // End of
      // while(trace_iter != trace_iter_end)
 
      // Each time we get here, that means that we have consumed all
      // the mass spectra data points that matched the current bin.
      // So go to the next bin.
 
      if(trace_iter == trace_iter_end)
        {
          // Well, in fact we truly need to exit the loops...
          break;
        }
 
      // qDebug() << "Now increment the bin iterator.";
      ++map_trace_bin_iter;
    }
  // End of
  // while(bin_iter != bin_end_iter)
 
  return map_trace;
}
 
MapTrace &
MassSpectrumPlusCombiner::combine(MapTrace &map_trace_out,
                                  const MapTrace &map_trace_in) const
{
  // qDebug();
 
  if(!map_trace_in.size())
    {
      // qDebug() << "Thread:" << QThread::currentThreadId()
      //<< "Returning right away because map_trace_in is empty.";
      return map_trace_out;
    }
 
  Trace trace(map_trace_in);
 
  return combine(map_trace_out, trace);
}
 
 
} // namespace pappso