plots/cartesian/XYSmoothCurve.cpp

0001 /*
0002     File                 : XYSmoothCurve.cpp
0003     Project              : LabPlot
0004     Description          : A xy-curve defined by a smooth
0005     --------------------------------------------------------------------
0006     SPDX-FileCopyrightText: 2016 Stefan Gerlach <stefan.gerlach@uni.kn>
0007     SPDX-FileCopyrightText: 2017 Alexander Semke <alexander.semke@web.de>
0008     SPDX-License-Identifier: GPL-2.0-or-later
0009 */
0010
0011 /*!
0012   \class XYSmoothCurve
0013   \brief A xy-curve defined by a smooth
0014
0015   \ingroup worksheet
0016 */
0017
0018 #include "XYSmoothCurve.h"
0019 #include "XYSmoothCurvePrivate.h"
0020 #include "backend/core/column/Column.h"
0021 #include "backend/lib/XmlStreamReader.h"
0022 #include "backend/lib/commandtemplates.h"
0023 #include "backend/lib/macros.h"
0024
0025 #include <KLocalizedString>
0026
0027 #include <QElapsedTimer>
0028 #include <QIcon>
0029 #include <QThreadPool>
0030
0031 extern "C" {
0032 #include "backend/nsl/nsl_sf_kernel.h"
0033 }
0034 #include "backend/nsl/nsl_stats.h"
0035 #include <gsl/gsl_math.h> // gsl_pow_*
0036
0037 XYSmoothCurve::XYSmoothCurve(const QString& name)
0038     : XYAnalysisCurve(name, new XYSmoothCurvePrivate(this), AspectType::XYSmoothCurve) {
0039 }
0040
0041 XYSmoothCurve::XYSmoothCurve(const QString& name, XYSmoothCurvePrivate* dd)
0042     : XYAnalysisCurve(name, dd, AspectType::XYSmoothCurve) {
0043 }
0044
0045 // no need to delete the d-pointer here - it inherits from QGraphicsItem
0046 // and is deleted during the cleanup in QGraphicsScene
0047 XYSmoothCurve::~XYSmoothCurve() = default;
0048
0049 void XYSmoothCurve::recalculate() {
0050     Q_D(XYSmoothCurve);
0051     d->recalculate();
0052 }
0053
0054 /*!
0055     Returns an icon to be used in the project explorer.
0056 */
0057 QIcon XYSmoothCurve::icon() const {
0058     return QIcon::fromTheme(QStringLiteral("labplot-xy-smoothing-curve"));
0059 }
0060
0061 // ##############################################################################
0062 // ##########################  getter methods  ##################################
0063 // ##############################################################################
0064 const AbstractColumn* XYSmoothCurve::roughsColumn() const {
0065     Q_D(const XYSmoothCurve);
0066     return d->roughColumn;
0067 }
0068
0069 BASIC_SHARED_D_READER_IMPL(XYSmoothCurve, XYSmoothCurve::SmoothData, smoothData, smoothData)
0070
0071 const XYAnalysisCurve::Result& XYSmoothCurve::result() const {
0072     Q_D(const XYSmoothCurve);
0073     return d->smoothResult;
0074 }
0075
0076 // ##############################################################################
0077 // #################  setter methods and undo commands ##########################
0078 // ##############################################################################
0079 STD_SETTER_CMD_IMPL_F_S(XYSmoothCurve, SetSmoothData, XYSmoothCurve::SmoothData, smoothData, recalculate)
0080 void XYSmoothCurve::setSmoothData(const XYSmoothCurve::SmoothData& smoothData) {
0081     Q_D(XYSmoothCurve);
0082     exec(new XYSmoothCurveSetSmoothDataCmd(d, smoothData, ki18n("%1: set options and perform the smooth")));
0083 }
0084
0085 // ##############################################################################
0086 // ######################### Private implementation #############################
0087 // ##############################################################################
0088 XYSmoothCurvePrivate::XYSmoothCurvePrivate(XYSmoothCurve* owner)
0089     : XYAnalysisCurvePrivate(owner)
0090     , q(owner) {
0091 }
0092
0093 // no need to delete xColumn and yColumn, they are deleted
0094 // when the parent aspect is removed
0095 XYSmoothCurvePrivate::~XYSmoothCurvePrivate() = default;
0096
0097 void XYSmoothCurvePrivate::resetResults() {
0098     smoothResult = XYAnalysisCurve::Result();
0099 }
0100
0101 bool XYSmoothCurvePrivate::recalculateSpecific(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn) {
0102     DEBUG(Q_FUNC_INFO)
0103     QElapsedTimer timer;
0104     timer.start();
0105
0106     if (roughVector)
0107         roughVector->clear();
0108
0109     if (!roughColumn) {
0110         roughColumn = new Column(QStringLiteral("rough"), AbstractColumn::ColumnMode::Double);
0111         roughColumn->setFixed(true); // visible in the project explorer but cannot be modified (renamed, deleted, etc.)
0112         roughVector = static_cast<QVector<double>*>(roughColumn->data());
0113         q->addChild(roughColumn);
0114     }
0115
0116     // check column sizes
0117     if (tmpXDataColumn->rowCount() != tmpYDataColumn->rowCount()) {
0118         smoothResult.available = true;
0119         smoothResult.valid = false;
0120         smoothResult.status = i18n("Number of x and y data points must be equal.");
0121         return true;
0122     }
0123
0124     // copy all valid data point for the smooth to temporary vectors
0125     QVector<double> xdataVector;
0126     QVector<double> ydataVector;
0127
0128     double xmin;
0129     double xmax;
0130     if (smoothData.autoRange) {
0131         xmin = tmpXDataColumn->minimum();
0132         xmax = tmpXDataColumn->maximum();
0133     } else {
0134         xmin = smoothData.xRange.first();
0135         xmax = smoothData.xRange.last();
0136     }
0137
0138     XYAnalysisCurve::copyData(xdataVector, ydataVector, tmpXDataColumn, tmpYDataColumn, xmin, xmax);
0139
0140     // number of data points to smooth
0141     const size_t n = (size_t)xdataVector.size();
0142     if (n < 2) {
0143         smoothResult.available = true;
0144         smoothResult.valid = false;
0145         smoothResult.status = i18n("Not enough data points available.");
0146         return true;
0147     }
0148
0149     double* xdata = xdataVector.data();
0150     double* ydata = ydataVector.data();
0151
0152     double* ydataOriginal = new double[n];
0153     memcpy(ydataOriginal, ydata, n * sizeof(double));
0154
0155     // smooth settings
0156     const nsl_smooth_type type = smoothData.type;
0157     const size_t points = smoothData.points;
0158     const nsl_smooth_weight_type weight = smoothData.weight;
0159     const double percentile = smoothData.percentile;
0160     const int order = smoothData.order;
0161     const nsl_smooth_pad_mode padMode = smoothData.mode;
0162     const double lvalue = smoothData.lvalue;
0163     const double rvalue = smoothData.rvalue;
0164
0165     DEBUG(" smooth type:" << nsl_smooth_type_name[type]);
0166     DEBUG(" points = " << points);
0167     DEBUG(" weight: " << nsl_smooth_weight_type_name[weight]);
0168     DEBUG(" percentile = " << percentile);
0169     DEBUG(" order = " << order);
0170     DEBUG(" pad mode =  " << nsl_smooth_pad_mode_name[padMode]);
0171     DEBUG(" const. values = " << lvalue << ' ' << rvalue);
0172
0173     ///////////////////////////////////////////////////////////
0174     int status = 0;
0175     gsl_set_error_handler_off();
0176
0177     switch (type) {
0178     case nsl_smooth_type_moving_average:
0179         status = nsl_smooth_moving_average(ydata, n, points, weight, padMode);
0180         break;
0181     case nsl_smooth_type_moving_average_lagged:
0182         status = nsl_smooth_moving_average_lagged(ydata, n, points, weight, padMode);
0183         break;
0184     case nsl_smooth_type_percentile:
0185         status = nsl_smooth_percentile(ydata, n, points, percentile, padMode);
0186         break;
0187     case nsl_smooth_type_savitzky_golay:
0188         if (padMode == nsl_smooth_pad_constant)
0189             nsl_smooth_pad_constant_set(lvalue, rvalue);
0190         status = nsl_smooth_savgol(ydata, n, points, order, padMode);
0191         break;
0192     }
0193
0194     xVector->resize((int)n);
0195     yVector->resize((int)n);
0196     memcpy(xVector->data(), xdata, n * sizeof(double));
0197     memcpy(yVector->data(), ydata, n * sizeof(double));
0198     ///////////////////////////////////////////////////////////
0199
0200     // write the result
0201     smoothResult.available = true;
0202     smoothResult.valid = (status == 0);
0203     smoothResult.status = QString::number(status);
0204     smoothResult.elapsedTime = timer.elapsed();
0205
0206     // fill rough vector
0207     if (roughVector) {
0208         roughVector->resize((int)n);
0209         for (int i = 0; i < (int)n; ++i)
0210             roughVector->data()[i] = ydataOriginal[i] - ydata[i];
0211         roughColumn->setChanged();
0212     }
0213
0214     delete[] ydataOriginal;
0215
0216     return true;
0217 }
0218
0219 // ##############################################################################
0220 // ##################  Serialization/Deserialization  ###########################
0221 // ##############################################################################
0222 //! Save as XML
0223 void XYSmoothCurve::save(QXmlStreamWriter* writer) const {
0224     Q_D(const XYSmoothCurve);
0225
0226     writer->writeStartElement(QStringLiteral("xySmoothCurve"));
0227
0228     // write the base class
0229     XYAnalysisCurve::save(writer);
0230
0231     // write xy-smooth-curve specific information
0232     //  smooth data
0233     writer->writeStartElement(QStringLiteral("smoothData"));
0234     writer->writeAttribute(QStringLiteral("autoRange"), QString::number(d->smoothData.autoRange));
0235     writer->writeAttribute(QStringLiteral("xRangeMin"), QString::number(d->smoothData.xRange.first()));
0236     writer->writeAttribute(QStringLiteral("xRangeMax"), QString::number(d->smoothData.xRange.last()));
0237     writer->writeAttribute(QStringLiteral("type"), QString::number(d->smoothData.type));
0238     writer->writeAttribute(QStringLiteral("points"), QString::number(d->smoothData.points));
0239     writer->writeAttribute(QStringLiteral("weight"), QString::number(d->smoothData.weight));
0240     writer->writeAttribute(QStringLiteral("percentile"), QString::number(d->smoothData.percentile));
0241     writer->writeAttribute(QStringLiteral("order"), QString::number(d->smoothData.order));
0242     writer->writeAttribute(QStringLiteral("mode"), QString::number(d->smoothData.mode));
0243     writer->writeAttribute(QStringLiteral("lvalue"), QString::number(d->smoothData.lvalue));
0244     writer->writeAttribute(QStringLiteral("rvalue"), QString::number(d->smoothData.rvalue));
0245     writer->writeEndElement(); // smoothData
0246
0247     // smooth results (generated columns)
0248     writer->writeStartElement(QStringLiteral("smoothResult"));
0249     writer->writeAttribute(QStringLiteral("available"), QString::number(d->smoothResult.available));
0250     writer->writeAttribute(QStringLiteral("valid"), QString::number(d->smoothResult.valid));
0251     writer->writeAttribute(QStringLiteral("status"), d->smoothResult.status);
0252     writer->writeAttribute(QStringLiteral("time"), QString::number(d->smoothResult.elapsedTime));
0253
0254     // save calculated columns if available
0255     if (saveCalculations() && d->xColumn) {
0256         d->xColumn->save(writer);
0257         d->yColumn->save(writer);
0258     }
0259
0260     if (d->roughColumn)
0261         d->roughColumn->save(writer);
0262
0263     writer->writeEndElement(); //"smoothResult"
0264
0265     writer->writeEndElement(); //"xySmoothCurve"
0266 }
0267
0268 //! Load from XML
0269 bool XYSmoothCurve::load(XmlStreamReader* reader, bool preview) {
0270     Q_D(XYSmoothCurve);
0271
0272     QXmlStreamAttributes attribs;
0273     QString str;
0274
0275     while (!reader->atEnd()) {
0276         reader->readNext();
0277         if (reader->isEndElement() && reader->name() == QLatin1String("xySmoothCurve"))
0278             break;
0279
0280         if (!reader->isStartElement())
0281             continue;
0282
0283         if (reader->name() == QLatin1String("xyAnalysisCurve")) {
0284             if (!XYAnalysisCurve::load(reader, preview))
0285                 return false;
0286         } else if (!preview && reader->name() == QLatin1String("smoothData")) {
0287             attribs = reader->attributes();
0288             READ_INT_VALUE("autoRange", smoothData.autoRange, bool);
0289             READ_DOUBLE_VALUE("xRangeMin", smoothData.xRange.first());
0290             READ_DOUBLE_VALUE("xRangeMax", smoothData.xRange.last());
0291             READ_INT_VALUE("type", smoothData.type, nsl_smooth_type);
0292             READ_INT_VALUE("points", smoothData.points, size_t);
0293             READ_INT_VALUE("weight", smoothData.weight, nsl_smooth_weight_type);
0294             READ_DOUBLE_VALUE("percentile", smoothData.percentile);
0295             READ_INT_VALUE("order", smoothData.order, int);
0296             READ_INT_VALUE("mode", smoothData.mode, nsl_smooth_pad_mode);
0297             READ_DOUBLE_VALUE("lvalue", smoothData.lvalue);
0298             READ_DOUBLE_VALUE("rvalue", smoothData.rvalue);
0299         } else if (!preview && reader->name() == QLatin1String("smoothResult")) {
0300             attribs = reader->attributes();
0301             READ_INT_VALUE("available", smoothResult.available, int);
0302             READ_INT_VALUE("valid", smoothResult.valid, int);
0303             READ_STRING_VALUE("status", smoothResult.status);
0304             READ_INT_VALUE("time", smoothResult.elapsedTime, int);
0305         } else if (!preview && reader->name() == QLatin1String("column")) {
0306             Column* column = new Column(QString(), AbstractColumn::ColumnMode::Double);
0307             if (!column->load(reader, preview)) {
0308                 delete column;
0309                 return false;
0310             }
0311             if (column->name() == QLatin1String("x"))
0312                 d->xColumn = column;
0313             else if (column->name() == QLatin1String("y"))
0314                 d->yColumn = column;
0315             else
0316                 d->roughColumn = column;
0317         } else { // unknown element
0318             reader->raiseUnknownElementWarning();
0319             if (!reader->skipToEndElement())
0320                 return false;
0321         }
0322     }
0323
0324     if (preview)
0325         return true;
0326
0327     // wait for data to be read before using the pointers
0328     QThreadPool::globalInstance()->waitForDone();
0329
0330     if (d->xColumn && d->yColumn) {
0331         d->xColumn->setHidden(true);
0332         addChild(d->xColumn);
0333
0334         d->yColumn->setHidden(true);
0335         addChild(d->yColumn);
0336
0337         d->xVector = static_cast<QVector<double>*>(d->xColumn->data());
0338         d->yVector = static_cast<QVector<double>*>(d->yColumn->data());
0339
0340         static_cast<XYCurvePrivate*>(d_ptr)->xColumn = d->xColumn;
0341         static_cast<XYCurvePrivate*>(d_ptr)->yColumn = d->yColumn;
0342
0343         recalcLogicalPoints();
0344     }
0345
0346     if (d->roughColumn) {
0347         addChild(d->roughColumn);
0348         d->roughVector = static_cast<QVector<double>*>(d->roughColumn->data());
0349     }
0350
0351     return true;
0352 }