Loading...
Searching...
No Matches
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
pappso::CubicSplineModel
Definition cubicsplinemodel.h:24
pappso::CubicSplineModel::m_quadraticCoeffs
QList< double > m_quadraticCoeffs
Definition cubicsplinemodel.h:51
pappso::CubicSplineModel::setup
void setup(const QList< double > &x_values, const QList< double > &y_values)
Definition cubicsplinemodel.cpp:89
pappso::CubicSplineModel::CubicSplineModel
CubicSplineModel()
Definition cubicsplinemodel.cpp:24
pappso::CubicSplineModel::~CubicSplineModel
virtual ~CubicSplineModel()
Definition cubicsplinemodel.cpp:83
pappso::CubicSplineModel::m_name
QString m_name
Definition cubicsplinemodel.h:48
pappso::CubicSplineModel::operator=
CubicSplineModel & operator=(const CubicSplineModel &other)
Definition cubicsplinemodel.cpp:147
pappso::CubicSplineModel::m_linearCoeffs
QList< double > m_linearCoeffs
Definition cubicsplinemodel.h:50
pappso::CubicSplineModel::derivatives
double derivatives(const double x_value, unsigned order) const
Definition cubicsplinemodel.cpp:203
pappso::CubicSplineModel::derivative
double derivative(const double x_value) const
Definition cubicsplinemodel.cpp:197
pappso::CubicSplineModel::getKnots
const QList< double > & getKnots() const
Definition cubicsplinemodel.cpp:163
pappso::CubicSplineModel::m_cubicCoeffs
QList< double > m_cubicCoeffs
Definition cubicsplinemodel.h:52
pappso::CubicSplineModel::m_constCoeffs
QList< double > m_constCoeffs
Definition cubicsplinemodel.h:49
pappso::CubicSplineModel::m_knots
QList< double > m_knots
Definition cubicsplinemodel.h:53
pappso::CubicSplineModel::clone
CubicSplineModel * clone(const CubicSplineModel &other)
Definition cubicsplinemodel.cpp:68
pappso::CubicSplineModel::evalSplineAt
double evalSplineAt(double x_value) const
Definition cubicsplinemodel.cpp:169
cubicsplinemodel.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappso::PeptideIonNter::a
@ a
Definition peptide.h:77
pappso::PeptideIonNter::b
@ b
Definition peptide.h:74
pappso::spline_bisection
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)
Definition cubicsplinemodel.cpp:241
pappso::PeptideIonCter::z
@ z
Definition peptide.h:90