trace.cpp

Go to the documentation of this file.

#include <numeric>
#include <limits>
#include <vector>
#include <map>
#include <cmath>
#include <algorithm>
 
#include <QDebug>
 
#include "trace.h"
#include "maptrace.h"
#include "../processing/combiners/tracepluscombiner.h"
#include "../processing/combiners/traceminuscombiner.h"
#include "../types.h"
#include "../pappsoexception.h"
#include "../exception/exceptionoutofrange.h"
#include "../processing/filters/filterresample.h"
#include "../processing/filters/filterpass.h"
 
 
int traceMetaTypeId    = qRegisterMetaType<pappso::Trace>("pappso::Trace");
int tracePtrMetaTypeId = qRegisterMetaType<pappso::Trace *>("pappso::Trace *");
 
 
namespace pappso
{
 
std::vector<DataPoint>::iterator
findFirstEqualOrGreaterX(std::vector<DataPoint>::iterator begin,
                         std::vector<DataPoint>::iterator end,
                         const double &value)
{
  return std::find_if(begin, end, [value](const DataPoint &to_compare) {
    if(to_compare.x < value)
      {
        return false;
      }
    return true;
  });
}
 
std::vector<DataPoint>::const_iterator
findFirstEqualOrGreaterX(std::vector<DataPoint>::const_iterator begin,
                         std::vector<DataPoint>::const_iterator end,
                         const double &value)
{
  return std::find_if(begin, end, [value](const DataPoint &to_compare) {
    if(to_compare.x < value)
      {
        return false;
      }
    return true;
  });
}
 
std::vector<DataPoint>::iterator
findFirstGreaterX(std::vector<DataPoint>::iterator begin,
                  std::vector<DataPoint>::iterator end,
                  const double &value)
{
  return std::find_if(begin, end, [value](const DataPoint &to_compare) {
    if(to_compare.x > value)
      {
        return true;
      }
    return false;
  });
}
 
std::vector<DataPoint>::const_iterator
findFirstGreaterX(std::vector<DataPoint>::const_iterator begin,
                  std::vector<DataPoint>::const_iterator end,
                  const double &value)
{
  return std::find_if(begin, end, [value](const DataPoint &to_compare) {
    if(to_compare.x > value)
      {
        return true;
      }
    return false;
  });
}
 
 
std::vector<DataPoint>::iterator
findDifferentYvalue(std::vector<DataPoint>::iterator begin,
                    std::vector<DataPoint>::iterator end,
                    const double &y_value)
{
  return std::find_if(begin, end, [y_value](const DataPoint &to_compare) {
    if(to_compare.y != y_value)
      {
        return true;
      }
    return false;
  });
}
 
std::vector<DataPoint>::const_iterator
findDifferentYvalue(std::vector<DataPoint>::const_iterator begin,
                    std::vector<DataPoint>::const_iterator end,
                    const double &y_value)
{
  return std::find_if(begin, end, [y_value](const DataPoint &to_compare) {
    if(to_compare.y != y_value)
      {
        return true;
      }
    return false;
  });
}
 
 
std::vector<DataPoint>::const_iterator
minYDataPoint(std::vector<DataPoint>::const_iterator begin,
              std::vector<DataPoint>::const_iterator end)
{
  return std::min_element(
    begin, end, [](const DataPoint &a, const DataPoint &b) {
      return a.y < b.y;
    });
}
 
 
std::vector<DataPoint>::const_iterator
maxYDataPoint(std::vector<DataPoint>::const_iterator begin,
              std::vector<DataPoint>::const_iterator end)
{
  return std::max_element(
    begin, end, [](const DataPoint &a, const DataPoint &b) {
      return a.y < b.y;
    });
}
 
 
std::vector<DataPoint>::iterator
maxYDataPoint(std::vector<DataPoint>::iterator begin,
              std::vector<DataPoint>::iterator end)
{
  return std::max_element(
    begin, end, [](const DataPoint &a, const DataPoint &b) {
      return a.y < b.y;
    });
}
 
 
std::vector<DataPoint>::const_iterator
moveLowerYRigthDataPoint(const Trace &trace,
                         std::vector<DataPoint>::const_iterator begin)
{
  if(begin == trace.end())
    return begin;
  auto it     = begin + 1;
  auto result = begin;
  while((it != trace.end()) && (it->y <= result->y))
    {
      it++;
      result++;
    }
  return result;
}
 
std::vector<DataPoint>::const_iterator
moveLowerYLeftDataPoint(const Trace &trace,
                        std::vector<DataPoint>::const_iterator begin)
{
  if(begin == trace.begin())
    return begin;
  auto it     = begin - 1;
  auto result = begin;
  while((it != trace.begin()) && (it->y <= result->y))
    {
      it--;
      result--;
    }
  return result;
}
 
 
double
sumYTrace(std::vector<DataPoint>::const_iterator begin,
          std::vector<DataPoint>::const_iterator end,
          double init)
{
  return std::accumulate(
    begin, end, init, [](double a, const DataPoint &b) { return a + b.y; });
}
 
double
meanYTrace(std::vector<DataPoint>::const_iterator begin,
           std::vector<DataPoint>::const_iterator end)
{
  pappso_double nb_element = distance(begin, end);
  if(nb_element == 0)
    throw ExceptionOutOfRange(
      QObject::tr("unable to compute mean on a trace of size 0"));
  return (sumYTrace(begin, end, 0) / nb_element);
}
 
double
medianYTrace(std::vector<DataPoint>::const_iterator begin,
             std::vector<DataPoint>::const_iterator end)
{
  pappso_double nb_element = distance(begin, end);
  if(nb_element == 0)
    throw ExceptionOutOfRange(
      QObject::tr("unable to compute median on a trace of size 0"));
 
  std::vector<DataPoint> data(begin, end);
  std::nth_element(
    data.begin(),
    data.begin() + data.size() / 2,
    data.end(),
    [](const DataPoint &a, const DataPoint &b) { return a.y < b.y; });
  return data[data.size() / 2].y;
}
 
double
areaTrace(std::vector<DataPoint>::const_iterator begin,
          std::vector<DataPoint>::const_iterator end)
{
 
  if(begin == end)
    return 0;
  auto previous = begin;
  auto next     = begin + 1;
  double area   = 0;
  while(next != end)
    {
      area += ((next->x - previous->x) * (previous->y + next->y)) / (double)2;
      previous++;
      next++;
    }
  return area;
}
 
 
Trace
flooredLocalMaxima(std::vector<DataPoint>::const_iterator begin,
                   std::vector<DataPoint>::const_iterator end,
                   double y_floor)
{
  Trace local_maxima_trace;
 
  Trace single_peak_trace;
 
  DataPoint previous_data_point;
 
  for(auto iter = begin; iter != end; ++iter)
    {
      DataPoint iterated_data_point(iter->x, iter->y);
 
      // qDebug().noquote() << "Current data point:"
      //<< iterated_data_point.toString();
 
      if(iterated_data_point.y < y_floor)
        {
          // qDebug() << "under the floor";
 
          if(single_peak_trace.size())
            {
              // qDebug() << "There was a single peak trace cooking";
 
              local_maxima_trace.push_back(single_peak_trace.maxYDataPoint());
 
              // qDebug().noquote() << "pushed back local maximum point:"
              //<< local_maxima_trace.back().toString();
 
              // Clean and set the context.
              single_peak_trace.clear();
 
              previous_data_point = iterated_data_point;
 
              continue;
            }
          else
            {
              // qDebug() << "no single peak trace cooking";
 
              previous_data_point = iterated_data_point;
 
              continue;
            }
        }
      else
        {
          // qDebug() << "over the floor";
 
          // The iterated value is greater than the y_floor value, so we need to
          // handle it.
 
          if(iterated_data_point.y == previous_data_point.y)
            {
              // We are in a flat region, no need to change anything to the
              // context, just skip the point.
              continue;
            }
          else if(iterated_data_point.y > previous_data_point.y)
            {
              // qDebug().noquote() << "ascending in a peak";
 
              // The previously iterated y value was smaller than the presently
              // iterated one, so we are ascending in a peak.
 
              // All we need to do is set the context.
 
              single_peak_trace.push_back(iterated_data_point);
 
              // qDebug().noquote() << "pushed back normal point:"
              //<< single_peak_trace.back().toString();
 
              previous_data_point = iterated_data_point;
 
              continue;
            }
          else
            {
              // qDebug().noquote() << "started descending in a peak";
 
              // No, the currently iterated y value is less than the previously
              // iterated value.
 
              single_peak_trace.push_back(iterated_data_point);
 
              // qDebug().noquote() << "pushed back normal point:"
              //<< single_peak_trace.back().toString();
 
              previous_data_point = iterated_data_point;
 
              continue;
            }
        }
    }
  // End of
  // for(auto iter = begin; iter != end; ++iter)
 
  // Attention, we might arrive here with a peak being created, we need to get
  // its maximum if that peak is non-empty;
 
  if(single_peak_trace.size())
    {
 
      local_maxima_trace.push_back(single_peak_trace.maxYDataPoint());
 
      // qDebug().noquote()
      //<< "was cooking a peak: pushed back local maximum point:"
      //<< local_maxima_trace.back().toString();
    }
 
  return local_maxima_trace;
}
 
 
Trace::Trace()
{
}
 
 
Trace::Trace(
  const std::vector<std::pair<pappso_double, pappso_double>> &dataPoints)
{
  reserve(dataPoints.size());
 
  for(auto &dataPoint : dataPoints)
    {
      push_back(DataPoint(dataPoint));
    }
 
  std::sort(begin(), end(), [](const DataPoint &a, const DataPoint &b) {
    return (a.y < b.y);
  });
}
 
 
Trace::Trace(const std::vector<DataPoint> &dataPoints)
  : std::vector<DataPoint>(dataPoints)
{
}
 
 
Trace::Trace(const std::vector<DataPoint> &&dataPoints)
  : std::vector<DataPoint>(std::move(dataPoints))
{
  // This constructor used by the MassSpectrum && constructor.
}
 
 
Trace::Trace(const MapTrace &map_trace)
{
  for(auto &&item : map_trace)
    push_back(DataPoint(item.first, item.second));
}
 
Trace::Trace(const Trace &other) : std::vector<DataPoint>(other)
{
}
 
 
Trace::Trace(const Trace &&other) : std::vector<DataPoint>(std::move(other))
{
  // This constructor used by the MassSpectrum && constructor.
}
 
 
Trace::~Trace()
{
  // Calls the destructor for each DataPoint object in the vector.
  clear();
}
 
 
size_t
Trace::initialize(const std::vector<pappso_double> &xVector,
                  const std::vector<pappso_double> &yVector)
{
  // Sanity check
  if(xVector.size() != yVector.size())
    qFatal(
      "Fatal error at %s@%d -- %s(). "
      "xVector and yVector must have the same size."
      "Program aborted.",
      __FILE__,
      __LINE__,
      __FUNCTION__);
 
  // Do not force the release of the all the vector space, because we prefer
  // resizing.  clear();
 
  resize(xVector.size());
 
  for(std::size_t iter = 0; iter < xVector.size(); ++iter)
    {
      push_back(DataPoint(xVector.at(iter), yVector.at(iter)));
    }
 
  return size();
}
 
 
size_t
Trace::initialize(const std::map<pappso_double, pappso_double> &map)
{
 
  // Do not force the release of the all the vector space, because we prefer
  // resizing.  clear(false);
 
  resize(map.size());
 
  for(auto &&item : map)
    {
      push_back(DataPoint(item.first, item.second));
    }
 
  return size();
}
 
 
size_t
Trace::initialize(const Trace &other)
{
  *this = other;
 
  return size();
}
 
 
Trace &
Trace::operator=(const Trace &other)
{
  assign(other.begin(), other.end());
 
  return *this;
}
 
 
Trace &
Trace::operator=(Trace &&other)
{
  vector<DataPoint>::operator=(std::move(other));
  return *this;
}
 
 
TraceSPtr
Trace::makeTraceSPtr() const
{
  return std::make_shared<Trace>(*this);
}
 
 
TraceCstSPtr
Trace::makeTraceCstSPtr() const
{
  return std::make_shared<const Trace>(*this);
}
 
 
std::vector<pappso_double>
Trace::xToVector() const
{
  std::vector<pappso_double> vector;
 
  for(auto &&dataPoint : *this)
    vector.push_back(dataPoint.x);
 
  return vector;
}
 
 
std::vector<pappso_double>
Trace::yToVector() const
{
  std::vector<pappso_double> vector;
 
  for(auto &&dataPoint : *this)
    vector.push_back(dataPoint.y);
 
  return vector;
}
 
 
std::map<pappso_double, pappso_double>
Trace::toMap() const
{
  std::map<pappso_double, pappso_double> map;
 
  std::pair<std::map<pappso_double, pappso_double>::iterator, bool> ret;
 
  for(auto &&dataPoint : *this)
    {
      ret = map.insert(
        std::pair<pappso_double, pappso_double>(dataPoint.x, dataPoint.y));
 
      if(ret.second == false)
        {
          qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << "()"
                   << "It is odd that the Trace contains multiple same keys.";
 
          // No insertion, then increment the y value.
          ret.first->second += dataPoint.y;
        }
    }
 
  return map;
}
 
 
// const DataPoint &
// Trace::dataPointWithX(pappso_double value) const
//{
// auto iterator =
// std::find_if(begin(), end(), [value](const DataPoint &dataPoint) {
// return (dataPoint.x == value);
//});
 
// if(iterator != end())
//{
//// The returned data point is valid.
// return *iterator;
//}
// else
//{
//// The returned data point is invalid because it is not initialized.
// return DataPoint();
//}
//}
 
 
std::vector<DataPoint>::iterator
Trace::dataPointIteratorxWithX(pappso_double value)
{
  auto iterator =
    std::find_if(begin(), end(), [value](const DataPoint &dataPoint) {
      return (dataPoint.x == value);
    });
 
  return iterator;
}
 
 
std::vector<DataPoint>::const_iterator
Trace::dataPointCstIteratorxWithX(pappso_double value) const
{
  auto iterator =
    std::find_if(begin(), end(), [value](const DataPoint &dataPoint) {
      return (dataPoint.x == value);
    });
 
  return iterator;
}
 
 
std::size_t
Trace::dataPointIndexWithX(pappso_double value) const
{
  std::vector<DataPoint>::const_iterator iterator =
    dataPointCstIteratorxWithX(value);
 
  if(iterator != end())
    return std::distance(begin(), iterator);
 
  return std::numeric_limits<std::size_t>::max();
}
 
 
DataPoint
Trace::containsX(pappso_double value, PrecisionPtr precision_p) const
{
  auto iterator = std::find_if(
    begin(), end(), [value, precision_p](const DataPoint &data_point) {
      if(precision_p)
        {
          pappso_double delta = precision_p->delta(value);
 
          if(data_point.x >= (value - delta) && data_point.x <= (value + delta))
            return true;
          else
            return false;
        }
      else
        {
          return (data_point.x == value);
        }
    });
 
  if(iterator != end())
    {
      // The returned data point is valid.
      return *iterator;
    }
  else
    {
      // The returned data point is invalid because it is not initialized.
      return DataPoint();
    }
}
 
 
const DataPoint &
Trace::minYDataPoint() const
{
  auto dataPoint = std::min_element(
    begin(), end(), [](const DataPoint &a, const DataPoint &b) {
      return (b.y < a.y);
    });
 
  if(dataPoint == end())
    {
      throw ExceptionOutOfRange(
        QObject::tr("unable to get min peak intensity on spectrum size %1")
          .arg(size()));
    }
 
  return (*dataPoint);
}
 
 
const DataPoint &
Trace::maxYDataPoint() const
{
  auto dataPoint = std::max_element(
    begin(), end(), [](const DataPoint &a, const DataPoint &b) {
      return (a.y < b.y);
    });
 
  if(dataPoint == end())
    {
      throw ExceptionOutOfRange(
        QObject::tr("unable to get max peak intensity on spectrum size %1")
          .arg(size()));
    }
 
  return (*dataPoint);
}
 
 
pappso_double
Trace::minY() const
{
  return minYDataPoint().y;
}
 
 
pappso_double
Trace::maxY() const
{
  return maxYDataPoint().y;
}
 
 
pappso_double
Trace::sumY() const
{
  // double sum = 0;
 
  // for(auto &&dp : m_dataPoints)
  // sum += dp.y;
 
  // qDebug() << __FILE__ << "@" << __LINE__ << __FUNCTION__ << " ()"
  //<< "Returning sum/tic:" << sum;
 
  // return sum;
 
  return std::accumulate(begin(),
                         end(),
                         (double)0,
                         [](pappso_double sum, const DataPoint &dataPoint) {
                           return (sum + dataPoint.y);
                         });
}
 
 
pappso_double
Trace::sumY(double mzStart, double mzEnd) const
{
  auto begin_it = findFirstEqualOrGreaterX(this->begin(), this->end(), mzStart);
  return sumYTrace(
    begin_it, findFirstGreaterX(begin_it, this->end(), mzEnd), 0);
}
 
 
pappso_double
Trace::maxY(double mzStart, double mzEnd) const
{
  std::vector<DataPoint>::const_iterator begin_it =
    findFirstEqualOrGreaterX(this->begin(), this->end(), mzStart);
 
  double max_y = 0;
 
  while(begin_it != this->end())
    {
      if(begin_it->y > max_y)
        max_y = begin_it->y;
      begin_it++;
    }
  return max_y;
}
 
 
void
Trace::sortX()
{
  std::sort(begin(), end(), [](const DataPoint &a, const DataPoint &b) {
    return (a.x < b.x);
  });
}
 
 
void
Trace::unique()
{
  auto last =
    std::unique(begin(), end(), [](const DataPoint &a, const DataPoint &b) {
      return (a.x == b.x);
    });
 
  erase(last, end());
}
 
 
std::vector<pappso_double>
Trace::xValues()
{
  std::vector<pappso_double> values;
 
  for(auto &&dataPoint : *this)
    {
      values.push_back(dataPoint.x);
    }
 
  return values;
}
 
 
std::vector<pappso_double>
Trace::yValues()
{
  std::vector<pappso_double> values;
 
  for(auto &&dataPoint : *this)
    {
      values.push_back(dataPoint.y);
    }
 
  return values;
}
 
 
QString
Trace::toString() const
{
  // Even if the spectrum is empty, we should return an empty string.
  QString text;
 
  for(auto &&dataPoint : *this)
    {
      text.append(QString("%1 %2\n")
                    .arg(dataPoint.x, 0, 'f', 10)
                    .arg(dataPoint.y, 0, 'f', 10));
    }
 
  return text;
}
 
 
Trace &
Trace::filter(const FilterInterface &filter)
{
  return filter.filter(*this);
}
 
} // namespace pappso