plots/cartesian/XYConvolutionCurve.cpp

0001 /*
0002     File                 : XYConvolutionCurve.cpp
0003     Project              : LabPlot
0004     Description          : A xy-curve defined by a convolution
0005     --------------------------------------------------------------------
0006     SPDX-FileCopyrightText: 2018 Stefan Gerlach <stefan.gerlach@uni.kn>
0007     SPDX-License-Identifier: GPL-2.0-or-later
0008 */
0009
0010 /*!
0011   \class XYConvolutionCurve
0012   \brief A xy-curve defined by a convolution
0013
0014   \ingroup worksheet
0015 */
0016
0017 #include "XYConvolutionCurve.h"
0018 #include "XYConvolutionCurvePrivate.h"
0019 #include "backend/core/column/Column.h"
0020 #include "backend/lib/XmlStreamReader.h"
0021 #include "backend/lib/commandtemplates.h"
0022 #include "backend/lib/macros.h"
0023
0024 #include <KLocalizedString>
0025 #include <QElapsedTimer>
0026 #include <QIcon>
0027 #include <QThreadPool>
0028
0029 #include <gsl/gsl_math.h>
0030
0031 XYConvolutionCurve::XYConvolutionCurve(const QString& name)
0032     : XYAnalysisCurve(name, new XYConvolutionCurvePrivate(this), AspectType::XYConvolutionCurve) {
0033 }
0034
0035 XYConvolutionCurve::XYConvolutionCurve(const QString& name, XYConvolutionCurvePrivate* dd)
0036     : XYAnalysisCurve(name, dd, AspectType::XYConvolutionCurve) {
0037 }
0038
0039 // no need to delete the d-pointer here - it inherits from QGraphicsItem
0040 // and is deleted during the cleanup in QGraphicsScene
0041 XYConvolutionCurve::~XYConvolutionCurve() = default;
0042
0043 void XYConvolutionCurve::recalculate() {
0044     Q_D(XYConvolutionCurve);
0045     d->recalculate();
0046 }
0047
0048 const XYAnalysisCurve::Result& XYConvolutionCurve::result() const {
0049     Q_D(const XYConvolutionCurve);
0050     return d->convolutionResult;
0051 }
0052
0053 /*!
0054     Returns an icon to be used in the project explorer.
0055 */
0056 QIcon XYConvolutionCurve::icon() const {
0057     //  return QIcon::fromTheme("labplot-xy-convolution-curve");//not available yet
0058     return QIcon::fromTheme(QStringLiteral("labplot-xy-curve"));
0059 }
0060
0061 // ##############################################################################
0062 // ##########################  getter methods  ##################################
0063 // ##############################################################################
0064 BASIC_SHARED_D_READER_IMPL(XYConvolutionCurve, XYConvolutionCurve::ConvolutionData, convolutionData, convolutionData)
0065
0066 const XYConvolutionCurve::ConvolutionResult& XYConvolutionCurve::convolutionResult() const {
0067     Q_D(const XYConvolutionCurve);
0068     return d->convolutionResult;
0069 }
0070
0071 // ##############################################################################
0072 // #################  setter methods and undo commands ##########################
0073 // ##############################################################################
0074 STD_SETTER_CMD_IMPL_F_S(XYConvolutionCurve, SetConvolutionData, XYConvolutionCurve::ConvolutionData, convolutionData, recalculate)
0075 void XYConvolutionCurve::setConvolutionData(const XYConvolutionCurve::ConvolutionData& convolutionData) {
0076     Q_D(XYConvolutionCurve);
0077     exec(new XYConvolutionCurveSetConvolutionDataCmd(d, convolutionData, ki18n("%1: set options and perform the convolution")));
0078 }
0079
0080 // ##############################################################################
0081 // ######################### Private implementation #############################
0082 // ##############################################################################
0083 XYConvolutionCurvePrivate::XYConvolutionCurvePrivate(XYConvolutionCurve* owner)
0084     : XYAnalysisCurvePrivate(owner)
0085     , q(owner) {
0086 }
0087
0088 // no need to delete xColumn and yColumn, they are deleted
0089 // when the parent aspect is removed
0090 XYConvolutionCurvePrivate::~XYConvolutionCurvePrivate() = default;
0091
0092 void XYConvolutionCurvePrivate::resetResults() {
0093     convolutionResult = XYConvolutionCurve::ConvolutionResult();
0094 }
0095
0096 bool XYConvolutionCurvePrivate::preparationValid(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn) {
0097     Q_UNUSED(tmpXDataColumn);
0098     return tmpYDataColumn != nullptr;
0099 }
0100
0101 bool XYConvolutionCurvePrivate::recalculateSpecific(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn) {
0102     QElapsedTimer timer;
0103     timer.start();
0104
0105     // determine the data source columns
0106     const AbstractColumn* tmpY2DataColumn = nullptr;
0107     if (dataSourceType == XYAnalysisCurve::DataSourceType::Spreadsheet) {
0108         // spreadsheet columns as data source
0109         tmpY2DataColumn = y2DataColumn;
0110     } else {
0111         // curve columns as data source
0112         // no y2 column: use standard kernel
0113     }
0114
0115     // copy all valid data point for the convolution to temporary vectors
0116     QVector<double> xdataVector;
0117     QVector<double> ydataVector;
0118     QVector<double> y2dataVector;
0119
0120     double xmin, xmax;
0121     if (tmpXDataColumn && convolutionData.autoRange) {
0122         xmin = tmpXDataColumn->minimum();
0123         xmax = tmpXDataColumn->maximum();
0124     } else {
0125         xmin = convolutionData.xRange.first();
0126         xmax = convolutionData.xRange.last();
0127     }
0128
0129     // only copy those data where values are valid and in range
0130     if (tmpXDataColumn) { // x-axis present (with possible range)
0131         for (int row = 0; row < tmpXDataColumn->rowCount(); ++row) {
0132             if (tmpXDataColumn->isValid(row) && !tmpXDataColumn->isMasked(row) && tmpYDataColumn->isValid(row) && !tmpYDataColumn->isMasked(row)) {
0133                 if (tmpXDataColumn->valueAt(row) >= xmin && tmpXDataColumn->valueAt(row) <= xmax) {
0134                     xdataVector.append(tmpXDataColumn->valueAt(row));
0135                     ydataVector.append(tmpYDataColumn->valueAt(row));
0136                 }
0137             }
0138         }
0139     } else { // no x-axis: take all valid values
0140         for (int row = 0; row < tmpYDataColumn->rowCount(); ++row)
0141             if (tmpYDataColumn->isValid(row) && !tmpYDataColumn->isMasked(row))
0142                 ydataVector.append(tmpYDataColumn->valueAt(row));
0143     }
0144
0145     const nsl_conv_kernel_type kernel = convolutionData.kernel;
0146     const size_t kernelSize = convolutionData.kernelSize;
0147     if (tmpY2DataColumn != nullptr) {
0148         for (int row = 0; row < tmpY2DataColumn->rowCount(); ++row)
0149             if (tmpY2DataColumn->isValid(row) && !tmpY2DataColumn->isMasked(row))
0150                 y2dataVector.append(tmpY2DataColumn->valueAt(row));
0151         DEBUG("kernel = given response");
0152     } else {
0153         DEBUG("kernel = " << nsl_conv_kernel_name[kernel] << ", size = " << kernelSize);
0154         double* k = new double[kernelSize];
0155         nsl_conv_standard_kernel(k, kernelSize, kernel);
0156         for (size_t i = 0; i < kernelSize; i++)
0157             y2dataVector.append(k[i]);
0158         delete[] k;
0159     }
0160
0161     const size_t n = (size_t)ydataVector.size(); // number of points for signal
0162     const size_t m = (size_t)y2dataVector.size(); // number of points for response
0163     if (n < 1 || m < 1) {
0164         convolutionResult.available = true;
0165         convolutionResult.valid = false;
0166         convolutionResult.status = i18n("Not enough data points available.");
0167         return true;
0168     }
0169
0170     double* xdata = xdataVector.data();
0171     double* ydata = ydataVector.data();
0172     double* y2data = y2dataVector.data();
0173
0174     // convolution settings
0175     const double samplingInterval = convolutionData.samplingInterval;
0176     const nsl_conv_direction_type direction = convolutionData.direction;
0177     const nsl_conv_type_type type = convolutionData.type;
0178     const nsl_conv_method_type method = convolutionData.method;
0179     const nsl_conv_norm_type norm = convolutionData.normalize;
0180     const nsl_conv_wrap_type wrap = convolutionData.wrap;
0181
0182     DEBUG("signal n = " << n << ", response m = " << m);
0183     DEBUG("sampling interval = " << samplingInterval);
0184     DEBUG("direction = " << nsl_conv_direction_name[direction]);
0185     DEBUG("type = " << nsl_conv_type_name[type]);
0186     DEBUG("method = " << nsl_conv_method_name[method]);
0187     DEBUG("norm = " << nsl_conv_norm_name[norm]);
0188     DEBUG("wrap = " << nsl_conv_wrap_name[wrap]);
0189
0190     ///////////////////////////////////////////////////////////
0191     size_t np;
0192     if (type == nsl_conv_type_linear)
0193         np = n + m - 1;
0194     else
0195         np = GSL_MAX(n, m);
0196
0197     double* out = (double*)malloc(np * sizeof(double));
0198     int status = nsl_conv_convolution_direction(ydata, n, y2data, m, direction, type, method, norm, wrap, out);
0199
0200     if (direction == nsl_conv_direction_backward)
0201         if (type == nsl_conv_type_linear)
0202             np = abs((int)(n - m)) + 1;
0203
0204     xVector->resize((int)np);
0205     yVector->resize((int)np);
0206     // take given x-axis values or use index
0207     if (tmpXDataColumn != nullptr) {
0208         int size = GSL_MIN(xdataVector.size(), (int)np);
0209         memcpy(xVector->data(), xdata, size * sizeof(double));
0210         double sampleInterval = (xVector->data()[size - 1] - xVector->data()[0]) / (xdataVector.size() - 1);
0211         DEBUG("xdata size = " << xdataVector.size() << ", np = " << np << ", sample interval = " << sampleInterval);
0212         for (int i = size; i < (int)np; i++) // fill missing values
0213             xVector->data()[i] = xVector->data()[size - 1] + (i - size + 1) * sampleInterval;
0214     } else { // fill with index (starting with 0)
0215         for (size_t i = 0; i < np; i++)
0216             xVector->data()[i] = i * samplingInterval;
0217     }
0218
0219     memcpy(yVector->data(), out, np * sizeof(double));
0220     free(out);
0221     ///////////////////////////////////////////////////////////
0222
0223     // write the result
0224     convolutionResult.available = true;
0225     convolutionResult.valid = (status == 0);
0226     convolutionResult.status = QString::number(status);
0227     convolutionResult.elapsedTime = timer.elapsed();
0228
0229     return true;
0230 }
0231
0232 // ##############################################################################
0233 // ##################  Serialization/Deserialization  ###########################
0234 // ##############################################################################
0235 //! Save as XML
0236 void XYConvolutionCurve::save(QXmlStreamWriter* writer) const {
0237     Q_D(const XYConvolutionCurve);
0238
0239     writer->writeStartElement(QStringLiteral("xyConvolutionCurve"));
0240
0241     // write the base class
0242     XYAnalysisCurve::save(writer);
0243
0244     // write xy-convolution-curve specific information
0245     //  convolution data
0246     writer->writeStartElement(QStringLiteral("convolutionData"));
0247     writer->writeAttribute(QStringLiteral("samplingInterval"), QString::number(d->convolutionData.samplingInterval));
0248     writer->writeAttribute(QStringLiteral("kernel"), QString::number(d->convolutionData.kernel));
0249     writer->writeAttribute(QStringLiteral("kernelSize"), QString::number(d->convolutionData.kernelSize));
0250     writer->writeAttribute(QStringLiteral("autoRange"), QString::number(d->convolutionData.autoRange));
0251     writer->writeAttribute(QStringLiteral("xRangeMin"), QString::number(d->convolutionData.xRange.first()));
0252     writer->writeAttribute(QStringLiteral("xRangeMax"), QString::number(d->convolutionData.xRange.last()));
0253     writer->writeAttribute(QStringLiteral("direction"), QString::number(d->convolutionData.direction));
0254     writer->writeAttribute(QStringLiteral("type"), QString::number(d->convolutionData.type));
0255     writer->writeAttribute(QStringLiteral("method"), QString::number(d->convolutionData.method));
0256     writer->writeAttribute(QStringLiteral("normalize"), QString::number(d->convolutionData.normalize));
0257     writer->writeAttribute(QStringLiteral("wrap"), QString::number(d->convolutionData.wrap));
0258     writer->writeEndElement(); // convolutionData
0259
0260     // convolution results (generated columns)
0261     writer->writeStartElement(QStringLiteral("convolutionResult"));
0262     writer->writeAttribute(QStringLiteral("available"), QString::number(d->convolutionResult.available));
0263     writer->writeAttribute(QStringLiteral("valid"), QString::number(d->convolutionResult.valid));
0264     writer->writeAttribute(QStringLiteral("status"), d->convolutionResult.status);
0265     writer->writeAttribute(QStringLiteral("time"), QString::number(d->convolutionResult.elapsedTime));
0266
0267     // save calculated columns if available
0268     if (saveCalculations() && d->xColumn) {
0269         d->xColumn->save(writer);
0270         d->yColumn->save(writer);
0271     }
0272     writer->writeEndElement(); //"convolutionResult"
0273
0274     writer->writeEndElement(); //"xyConvolutionCurve"
0275 }
0276
0277 //! Load from XML
0278 bool XYConvolutionCurve::load(XmlStreamReader* reader, bool preview) {
0279     Q_D(XYConvolutionCurve);
0280
0281     QXmlStreamAttributes attribs;
0282     QString str;
0283
0284     while (!reader->atEnd()) {
0285         reader->readNext();
0286         if (reader->isEndElement() && reader->name() == QLatin1String("xyConvolutionCurve"))
0287             break;
0288
0289         if (!reader->isStartElement())
0290             continue;
0291
0292         if (reader->name() == QLatin1String("xyAnalysisCurve")) {
0293             if (!XYAnalysisCurve::load(reader, preview))
0294                 return false;
0295         } else if (!preview && reader->name() == QLatin1String("convolutionData")) {
0296             attribs = reader->attributes();
0297             READ_DOUBLE_VALUE("samplingInterval", convolutionData.samplingInterval);
0298             READ_INT_VALUE("kernel", convolutionData.kernel, nsl_conv_kernel_type);
0299             READ_INT_VALUE("kernelSize", convolutionData.kernelSize, size_t);
0300             READ_INT_VALUE("autoRange", convolutionData.autoRange, bool);
0301             READ_DOUBLE_VALUE("xRangeMin", convolutionData.xRange.first());
0302             READ_DOUBLE_VALUE("xRangeMax", convolutionData.xRange.last());
0303             READ_INT_VALUE("direction", convolutionData.direction, nsl_conv_direction_type);
0304             READ_INT_VALUE("type", convolutionData.type, nsl_conv_type_type);
0305             READ_INT_VALUE("method", convolutionData.method, nsl_conv_method_type);
0306             READ_INT_VALUE("normalize", convolutionData.normalize, nsl_conv_norm_type);
0307             READ_INT_VALUE("wrap", convolutionData.wrap, nsl_conv_wrap_type);
0308         } else if (!preview && reader->name() == QLatin1String("convolutionResult")) {
0309             attribs = reader->attributes();
0310             READ_INT_VALUE("available", convolutionResult.available, int);
0311             READ_INT_VALUE("valid", convolutionResult.valid, int);
0312             READ_STRING_VALUE("status", convolutionResult.status);
0313             READ_INT_VALUE("time", convolutionResult.elapsedTime, int);
0314         } else if (!preview && reader->name() == QLatin1String("column")) {
0315             Column* column = new Column(QString(), AbstractColumn::ColumnMode::Double);
0316             if (!column->load(reader, preview)) {
0317                 delete column;
0318                 return false;
0319             }
0320             if (column->name() == QLatin1String("x"))
0321                 d->xColumn = column;
0322             else if (column->name() == QLatin1String("y"))
0323                 d->yColumn = column;
0324         } else { // unknown element
0325             reader->raiseUnknownElementWarning();
0326             if (!reader->skipToEndElement())
0327                 return false;
0328         }
0329     }
0330
0331     if (preview)
0332         return true;
0333
0334     // wait for data to be read before using the pointers
0335     QThreadPool::globalInstance()->waitForDone();
0336
0337     if (d->xColumn && d->yColumn) {
0338         d->xColumn->setHidden(true);
0339         addChild(d->xColumn);
0340
0341         d->yColumn->setHidden(true);
0342         addChild(d->yColumn);
0343
0344         d->xVector = static_cast<QVector<double>*>(d->xColumn->data());
0345         d->yVector = static_cast<QVector<double>*>(d->yColumn->data());
0346
0347         static_cast<XYCurvePrivate*>(d_ptr)->xColumn = d->xColumn;
0348         static_cast<XYCurvePrivate*>(d_ptr)->yColumn = d->yColumn;
0349
0350         recalcLogicalPoints();
0351     }
0352
0353     return true;
0354 }