cubicsplinemodel.cpp

Go to the documentation of this file.

// Copyright 2026 Filippo Rusconi
// Inspired by work by Lars Nilse in OpenMS
 
 
/////////////////////// stdlib includes
#include <limits>
#include <algorithm>
 
 
/////////////////////// Qt includes
#include <QList>
#include <QDebug>
 
 
/////////////////////// pappsomspp includes
 
 
/////////////////////// Local includes
#include "pappsomspp/core/processing/detection/cubicsplinemodel.h"
 
namespace pappso
{
 
CubicSplineModel::CubicSplineModel()
{
  // qDebug() << "Allocating:" << this;
}
 
CubicSplineModel::CubicSplineModel(const QList<double> &x_values,
                                   const QList<double> &y_values)
{
  // qDebug() << "Allocating:" << this;
 
  Q_ASSERT(x_values.size() == y_values.size());
  Q_ASSERT(x_values.size() >= 2);
 
  double eps = std::numeric_limits<double>::epsilon();
 
  // Ensure the x_values are sorted increasingly.
 
  Q_ASSERT(std::adjacent_find(
             x_values.cbegin(), x_values.cend(), [=](double a, double b) {
               return !(b > a + eps);
             }) == x_values.cend());
 
  setup(x_values, y_values);
 
  // qDebug() << "Number of knots:" << m_knots.size();
}
 
// Use delegating constructor.
CubicSplineModel::CubicSplineModel(const QMap<double, double> &x_y_values_map)
  : CubicSplineModel(x_y_values_map.keys(), x_y_values_map.values())
{
}
 
CubicSplineModel::CubicSplineModel(const CubicSplineModel &other)
{
  // qDebug() << "Allocating:" << this;
  m_constCoeffs     = other.m_constCoeffs;
  m_linearCoeffs    = other.m_linearCoeffs;
  m_quadraticCoeffs = other.m_quadraticCoeffs;
  m_cubicCoeffs     = other.m_cubicCoeffs;
  m_knots           = other.m_knots;
}
 
CubicSplineModel *
CubicSplineModel::clone(const CubicSplineModel &other)
{
  // qDebug() << "Cloning.";
  CubicSplineModel *copy_p = new CubicSplineModel();
 
  copy_p->m_name            = other.m_name;
  copy_p->m_constCoeffs     = other.m_constCoeffs;
  copy_p->m_linearCoeffs    = other.m_linearCoeffs;
  copy_p->m_quadraticCoeffs = other.m_quadraticCoeffs;
  copy_p->m_cubicCoeffs     = other.m_cubicCoeffs;
  copy_p->m_knots           = other.m_knots;
 
  return copy_p;
}
 
CubicSplineModel::~CubicSplineModel()
{
  // qDebug() << "Destructing the cubic spline model";
}
 
void
CubicSplineModel::setup(const QList<double> &x_values,
                        const QList<double> &y_values)
{
  // qDebug() << "Setting up with x_values size:" << x_values.size();
 
  Q_ASSERT(x_values.size() == y_values.size());
 
  const size_t n = x_values.size() - 1;
 
  std::vector<double> h;
  h.reserve(n);
  m_constCoeffs.reserve(n);
  m_knots.reserve(n + 1);
  // do the 0'th element manually, since the loop below only starts at 1
  h.push_back(x_values[1] - x_values[0]);
  m_knots.push_back(x_values[0]);
  m_constCoeffs.push_back(y_values[0]);
 
  std::vector<double> mu(n, 0.0);
  std::vector<double> z(n, 0.0);
  for(unsigned i = 1; i < n; ++i)
    {
      h.push_back(x_values[i + 1] - x_values[i]);
      const double l =
        2 * (x_values[i + 1] - x_values[i - 1]) - h[i - 1] * mu[i - 1];
      mu[i] = h[i] / l;
      z[i]  = (3 *
                (y_values[i + 1] * h[i - 1] -
                 y_values[i] * (x_values[i + 1] - x_values[i - 1]) +
                 y_values[i - 1] * h[i]) /
                (h[i - 1] * h[i]) -
              h[i - 1] * z[i - 1]) /
             l;
      // store x,y -- required for evaluation later on
      m_knots.push_back(x_values[i]);
      m_constCoeffs.push_back(y_values[i]);
    }
  // 'm_knots' needs to be full length (all other member vectors (except
  // m_quadraticCoeffs) are one element shorter)
  m_knots.push_back(x_values[n]);
 
  m_linearCoeffs.resize(n);
  m_cubicCoeffs.resize(n);
  m_quadraticCoeffs.resize(n + 1);
  m_quadraticCoeffs.back() = 0;
  for(int j = static_cast<int>(n) - 1; j >= 0; --j)
    {
      m_quadraticCoeffs[j] = z[j] - mu[j] * m_quadraticCoeffs[j + 1];
      m_linearCoeffs[j] =
        (y_values[j + 1] - y_values[j]) / h[j] -
        h[j] * (m_quadraticCoeffs[j + 1] + 2 * m_quadraticCoeffs[j]) / 3;
      m_cubicCoeffs[j] =
        (m_quadraticCoeffs[j + 1] - m_quadraticCoeffs[j]) / (3 * h[j]);
    }
  // qDebug() << "m_knots size:" << m_knots.size();
}
 
CubicSplineModel &
CubicSplineModel::operator=(const CubicSplineModel &other)
{
  if(&other == this)
    return *this;
 
  m_name            = other.m_name;
  m_constCoeffs     = other.m_constCoeffs;
  m_linearCoeffs    = other.m_linearCoeffs;
  m_quadraticCoeffs = other.m_quadraticCoeffs;
  m_cubicCoeffs     = other.m_cubicCoeffs;
  m_knots           = other.m_knots;
 
  return *this;
}
 
const QList<double> &
CubicSplineModel::getKnots() const
{
  return m_knots;
}
 
double
CubicSplineModel::evalSplineAt(double x_value) const
{
  Q_ASSERT(x_value >= m_knots.front() || x_value <= m_knots.back());
 
  // What is the index of the closest knot left of x_value (or the same).
  unsigned i = static_cast<unsigned>(
    std::lower_bound(m_knots.begin(), m_knots.end(), x_value) -
    m_knots.begin());
 
  if(m_knots[i] > x_value || m_knots.back() == x_value)
    {
      --i;
    }
 
  const double delta_x = x_value - m_knots[i];
 
  double eval_value =
    ((m_cubicCoeffs[i] * delta_x + m_quadraticCoeffs[i]) * delta_x +
     m_linearCoeffs[i]) *
      delta_x +
    m_constCoeffs[i];
 
  // qDebug() << "Now returning evaluation value:" << eval_value;
 
  return eval_value;
}
 
double
CubicSplineModel::derivative(const double x_value) const
{
  return derivatives(x_value, 1);
}
 
double
CubicSplineModel::derivatives(const double x_value, unsigned order) const
{
  // qDebug() << "Size of m_knots:" << m_knots.size();
 
  Q_ASSERT(x_value >= m_knots.front() && x_value <= m_knots.back());
 
  Q_ASSERT(order >= 1 && order <= 3);
 
  // determine index of closest node left of (or exactly at) x_value
  unsigned i = static_cast<unsigned>(
    std::lower_bound(m_knots.begin(), m_knots.end(), x_value) -
    m_knots.begin());
 
  if(m_knots[i] > x_value ||
     m_knots.back() ==
       x_value) // also, i must not point to last index in 'm_knots',
                // since all other vectors are one element shorter
    {
      --i;
    }
 
  const double delta_x = x_value - m_knots[i];
  if(order == 1)
    {
      return m_linearCoeffs[i] + 2 * m_quadraticCoeffs[i] * delta_x +
             3 * m_cubicCoeffs[i] * delta_x * delta_x;
    }
  else if(order == 2)
    {
      return 2 * m_quadraticCoeffs[i] + 6 * m_cubicCoeffs[i] * delta_x;
    }
  else
    {
      return 6 * m_cubicCoeffs[i];
    }
}
 
void
spline_bisection(const CubicSplineModel &spline_model,
                 double const mz_at_left,
                 double const mz_at_right,
                 double &center_peak_mz,
                 double &center_peak_intensity,
                 double const threshold)
{
  // qDebug() << "spline model has" << spline_model.getKnots().size();
 
  // calculate maximum by evaluating the spline's 1st derivative
  // (bisection method)
  double lefthand  = mz_at_left;
  double righthand = mz_at_right;
 
  bool lefthand_sign = true;
  double eps         = std::numeric_limits<double>::epsilon();
 
  // bisection
  do
    {
      double mid                = (lefthand + righthand) / 2.0;
      double midpoint_deriv_val = spline_model.derivative(mid);
 
      // if deriv nearly zero then maximum already found
      if(!(std::fabs(midpoint_deriv_val) > eps))
        {
          break;
        }
 
      bool midpoint_sign = (midpoint_deriv_val < 0.0) ? false : true;
 
      if(lefthand_sign ^ midpoint_sign)
        {
          righthand = mid;
        }
      else
        {
          lefthand = mid;
        }
    }
  while(righthand - lefthand > threshold);
 
  center_peak_mz = (lefthand + righthand) / 2;
 
  center_peak_intensity = spline_model.evalSplineAt(center_peak_mz);
 
  // qDebug() << "center_peak_intensity:" << center_peak_intensity;
}
 
} // namespace pappso