massspectrumminuscombiner.cpp

Go to the documentation of this file.

/////////////////////// StdLib includes
#include <numeric>
#include <limits>
#include <vector>
#include <map>
#include <cmath>
#include <iostream>
#include <iomanip>
 
 
/////////////////////// Qt includes
#include <QDebug>
#include <QFile>
#include <QThread>
#if 0
// For debugging purposes.
#include <QFile>
#endif
 
 
/////////////////////// Local includes
#include "massspectrumminuscombiner.h"
#include "../../types.h"
#include "pappsomspp/core/utils.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "pappsomspp/core/exception/exceptionnotpossible.h"
 
 
namespace pappso
{
 
 
//! Construct an uninitialized instance.
MassSpectrumMinusCombiner::MassSpectrumMinusCombiner()
{
}
 
 
MassSpectrumMinusCombiner::MassSpectrumMinusCombiner(int decimal_places)
  : MassSpectrumCombiner(decimal_places)
{
}
 
 
MassSpectrumMinusCombiner::MassSpectrumMinusCombiner(
  const MassSpectrumMinusCombiner &other)
  : MassSpectrumCombiner(other)
 
{
}
 
 
MassSpectrumMinusCombiner::MassSpectrumMinusCombiner(
  MassSpectrumMinusCombinerCstSPtr other)
  : MassSpectrumCombiner(other)
 
{
}
 
 
//! Destruct the instance.
MassSpectrumMinusCombiner::~MassSpectrumMinusCombiner()
{
}
 
 
MassSpectrumMinusCombiner &
MassSpectrumMinusCombiner::operator=(const MassSpectrumMinusCombiner &other)
{
  if(this == &other)
    return *this;
 
  m_decimalPlaces = other.m_decimalPlaces;
 
  m_bins.assign(other.m_bins.begin(), other.m_bins.end());
 
  return *this;
}
 
 
MapTrace &
MassSpectrumMinusCombiner::combine(MapTrace &map_trace,
                                   const Trace &trace) const
{
  // qDebug();
 
  if(!trace.size())
    {
      // qDebug() << "Thread:" << QThread::currentThreadId()
      //<< "Returning right away because trace is empty.";
      return map_trace;
    }
 
  // 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();
 
  // 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."));
 
  using BinIter = std::vector<pappso_double>::const_iterator;
 
  BinIter bin_iter     = m_bins.begin();
  BinIter bin_end_iter = m_bins.end();
 
  // qDebug() << "initial bins iter at a distance of:"
  //<< std::distance(m_bins.begin(), bin_iter)
  //<< "bins distance:" << std::distance(m_bins.begin(), m_bins.end())
  //<< "bins size:" << m_bins.size() << "first bin:" << m_bins.front()
  //<< "last bin:" << m_bins.back();
 
  // Iterate in the vector of bins and for each bin check if there are matching
  // data points in the trace.
 
  pappso_double current_bin_mz_value = 0;
 
  pappso_double trace_x = 0;
  pappso_double trace_y = 0;
 
  // Lower bound returns an iterator pointing to the first element in the
  // range [first, last) that is not less than (i.e. greater or equal to)
  // value, or last if no such element is found.
 
  // Get the first bin that would match the very first data point in the trace
  // that we need to combine into the map_trace.
 
  auto bin_iter_for_mz = lower_bound(bin_iter, bin_end_iter, trace_iter->x);
 
  if(bin_iter_for_mz != bin_end_iter)
    {
      // If we found a bin and that the bin is not the first bin of the bin
      // vector, then go back one bin to be sure we do not miss data points.
 
      if(bin_iter_for_mz != m_bins.begin())
        bin_iter = --bin_iter_for_mz;
    }
  else
    throw(ExceptionNotPossible("The bin vector must match the mz value."));
 
  // Now iterate in the bin vector starting from the std::prev(found bin).
 
  while(bin_iter != bin_end_iter)
    {
      current_bin_mz_value = *bin_iter;
 
      // qDebug() << "Iterating in new bin:"
      //<< QString("%1").arg(current_bin_mz_value, 0, 'f', 15)
      //<< "at a distance from the first bin of:"
      //<< std::distance(m_bins.begin(), bin_iter);
 
      // For the current bin, we start by instantiating a new DataPoint. By
      // essence, each bin will have at most one corresponding DataPoint.
 
      DataPoint bin_data_point;
 
      // Do not set the y value to 0 so that we can actually test if the
      // data point is valid later on (try not to push back y=0 data
      // points).
 
      bin_data_point.x = current_bin_mz_value;
 
      // qDebug() << "Seed bin_data_point.x with current_bin_mz_value value:"
      //<< QString("%1").arg(bin_data_point.x, 0, 'f', 15);
 
      // Now perform a loop over the data points in the mass spectrum.
 
      while(trace_iter != trace_iter_end)
        {
          bool trace_matched = false;
 
          // If we are not at the end of trace and if the y value of the
          // currently iterated trace datapoint 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() << "Iterating in new trace datapoint:"
          //<< trace_iter->toString(15) << "at distance:"
          //<< std::distance(trace_iter_begin, trace_iter);
 
          if(trace_iter->y)
            {
              // qDebug() << "The y value of trace iterated value is not 0.";
 
              // trace_x is the m/z value that we need to combine,
 
              trace_x = trace_iter->x;
              trace_y = trace_iter->y;
 
              // Now apply the rounding (if any).
              if(m_decimalPlaces != -1)
                trace_x = Utils::roundToDecimals(trace_x, m_decimalPlaces);
 
              if(trace_x <= current_bin_mz_value)
                {
                  // qDebug()
                  //<< "trace_x <= current_bin_mz_value: the match happened.";
 
                  // qDebug() << "Going to compute" << bin_data_point.y << "-"
                  //<< trace_y << "with current_bin_mz_value:"
                  //<< QString("%1").arg(
                  // current_bin_mz_value, 0, 'f', 15);
 
                  // This is where the bin_data_point actually gets decremented
                  // by the y value of the currently iterated trace data point.
                  // So bin_data_point has already the MINUS'ed value in it.
                  // Relate this code line with the result.first->second +=
                  // bin_data_point.y; code line elsewhere in this function.
                  bin_data_point.y -= trace_y;
 
                  // qDebug() << "New bin_data_point.y value:" <<
                  // bin_data_point.y;
 
                  // Let's record that we matched.
                  trace_matched = true;
 
                  // Because we matched, we can step-up with the
                  // iterator.
 
                  // qDebug() << "The values matched, increment the trace "
                  //"itereator but not the bin iterator, as maybe "
                  //"the next trace datapoint also matches the bin.";
 
                  ++trace_iter;
                }
              // else
              //{
              // We did have a non-0 y value, but that did not
              // match. So we do not step-up with the iterator because that
              // trace data point will fit into another bin not yet iterated
              // into.
              //}
            }
          // End of
          // if(trace_iter->y)
          else
            {
              // We iterated into a y=0 data point, so just skip it. Let the
              // below code think that we have matched the point and iterate one
              // step up.
 
              // qDebug() << "The y value of the currently iterated trace data "
              //"point is almost equal to 0, increment the "
              //"trace iter but do nothing else.";
 
              trace_matched = true;
              ++trace_iter;
            }
 
          // At this point, check if the trace data point matched the currently
          // iterated bin.
 
          if(!trace_matched)
            {
 
              // qDebug() << "The currently iterated trace datapoint did not "
              //"match the current bin.";
 
              // The trace data point did not match the currently iterated bin.
              // All we have to do is go to the next bin. We break and the bin
              // vector iterator will be incremented.
 
              // However, if we had a valid new data point, that
              // data point needs to be pushed back in the new mass
              // spectrum.
 
              if(bin_data_point.isValid())
                {
 
                  // qDebug() << "But we had a valid bin data point cooking, so
                  // " "we have to take that into account:"
                  //<< bin_data_point.toString(15);
 
                  // We need to check if that bin value is present already in
                  // the map_trace object passed as parameter.
 
                  std::pair<std::map<pappso_double, pappso_double>::iterator,
                            bool>
                    result =
                      map_trace.insert(std::pair<pappso_double, pappso_double>(
                        bin_data_point.x, -bin_data_point.y));
 
                  if(!result.second)
                    {
                      // qDebug()
                      //<< "The x value (mz) was already present in the "
                      //"map trace, simply update its y value (intensity).";
 
                      // The key already existed! The item was not inserted. We
                      // need to update the y value.
 
                      // qDebug() << "Going to compute" << result.first->second
                      //<< "+" << bin_data_point.y;
 
                      // This might be counterintuitive: one would expect the
                      // use of the '-=' operator because we are
                      // MINUS-combining. But no! bin_data_point only get y
                      // values already minus'ed. So we need to add that
                      // minus'ed value, otherwise we add the value (- -value is
                      // the same as +value).
                      result.first->second += bin_data_point.y;
 
                      // qDebug() << "New y value for the item in the map trace
                      // "
                      //"(result.first->second value):"
                      //<< result.first->second;
                    }
                  // else
                  //{
                  // qDebug()
                  //<< "The data point has a mz value that was not present "
                  //"in the map trace."
                  //<< " Inserted a new data point into the map trace:"
                  //<< bin_data_point.x << "," << bin_data_point.y;
                  //}
                }
 
              // We need to break this loop! That will increment the
              // 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)
        {
 
          // Make sure a last check is done in case one data point was
          // cooking...
 
          if(bin_data_point.isValid())
            {
 
              std::pair<std::map<pappso_double, pappso_double>::iterator, bool>
                result =
                  map_trace.insert(std::pair<pappso_double, pappso_double>(
                    bin_data_point.x, -bin_data_point.y));
 
              if(!result.second)
                {
                  // qDebug()
                  //<< "The x value (mz) was already present in the "
                  //"map trace, simply update its y value (intensity).";
 
                  // The key already existed! The item was not inserted. We
                  // need to update the y value.
 
                  // qDebug() << "Going to compute" << result.first->second <<
                  // "+"
                  //<< bin_data_point.y;
 
 
                  // This might be counterintuitive: one would expect the
                  // use of the '-=' operator because we are
                  // MINUS-combining. But no! bin_data_point only get y
                  // values already minus'ed. So we need to add that
                  // minus'ed value, otherwise we add the value (- -value is
                  // the same as +value).
                  result.first->second += bin_data_point.y;
 
                  // qDebug() << "New y value for the item in the map trace "
                  //"(result.first->second value):"
                  //<< result.first->second;
                }
              // else
              //{
              // qDebug()
              //<< "The data point has a mz value that was not present "
              //"in the map trace."
              //<< " Inserted a new data point into the map trace:"
              //<< bin_data_point.x << "," << bin_data_point.y;
              //}
            }
 
          // Now truly exit the loops...
          break;
        }
 
      ++bin_iter;
    }
  // End of
  // while(bin_iter != bin_end_iter)
 
  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
  //<< "The combination result mass spectrum being returned has TIC:"
  //<< new_trace.sumY();
 
  return map_trace;
}
 
 
MapTrace &
MassSpectrumMinusCombiner::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