diff --git a/tools/lammps-gui/chartviewer.cpp b/tools/lammps-gui/chartviewer.cpp index 413a169468..c8fe04ffef 100644 --- a/tools/lammps-gui/chartviewer.cpp +++ b/tools/lammps-gui/chartviewer.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -32,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -54,6 +56,20 @@ ChartWindow::ChartWindow(const QString &_filename, QWidget *parent) : // workaround for incorrect highlight bug on macOS auto *dummy = new QPushButton(QIcon(), ""); dummy->hide(); + + smooth = new QCheckBox("Show smooth graph"); + smooth->setCheckState(Qt::Unchecked); + window = new QSpinBox; + window->setRange(5, 100); + window->setValue(10); + window->setEnabled(false); + window->setToolTip("Smoothing Window Size"); + order = new QSpinBox; + order->setRange(1, 10); + order->setValue(4); + order->setEnabled(false); + order->setToolTip("Smoothing Order"); + auto *normal = new QPushButton(QIcon(":/icons/gtk-zoom-fit.png"), ""); normal->setToolTip("Reset zoom to normal"); @@ -61,7 +77,12 @@ ChartWindow::ChartWindow(const QString &_filename, QWidget *parent) : top->addWidget(menu); top->addSpacerItem(new QSpacerItem(1, 1, QSizePolicy::Expanding, QSizePolicy::Minimum)); top->addWidget(dummy); + top->addWidget(smooth); + top->addWidget(window); + top->addWidget(order); + top->addWidget(new QLabel(" ")); top->addWidget(normal); + top->addWidget(new QLabel(" ")); top->addWidget(new QLabel("Select data:")); top->addWidget(columns); saveAsAct = file->addAction("&Save Graph As...", this, &ChartWindow::saveAs); @@ -86,6 +107,11 @@ ChartWindow::ChartWindow(const QString &_filename, QWidget *parent) : layout->addLayout(top); setLayout(layout); + connect(smooth, &QPushButton::released, this, &ChartWindow::update_smooth); + connect(window, &QAbstractSpinBox::editingFinished, this, &ChartWindow::update_smooth); + connect(order, &QAbstractSpinBox::editingFinished, this, &ChartWindow::update_smooth); + connect(window, &QSpinBox::valueChanged, this, &ChartWindow::update_smooth); + connect(order, &QSpinBox::valueChanged, this, &ChartWindow::update_smooth); connect(normal, &QPushButton::released, this, &ChartWindow::reset_zoom); connect(columns, SIGNAL(currentIndexChanged(int)), this, SLOT(change_chart(int))); installEventFilter(this); @@ -160,6 +186,19 @@ void ChartWindow::stop_run() if (main) main->stop_run(); } +void ChartWindow::update_smooth() +{ + bool do_smooth = smooth->isChecked(); + window->setEnabled(do_smooth); + order->setEnabled(do_smooth); + + int wval = window->value(); + int oval = order->value(); + + for (auto &c : charts) + c->smooth_param(do_smooth, wval, oval); +} + void ChartWindow::saveAs() { if (charts.empty()) return; @@ -317,8 +356,9 @@ bool ChartWindow::eventFilter(QObject *watched, QEvent *event) /* -------------------------------------------------------------------- */ ChartViewer::ChartViewer(const QString &title, int _index, QWidget *parent) : - QChartView(parent), last_step(-1), index(_index), chart(new QChart), series(new QLineSeries), - xaxis(new QValueAxis), yaxis(new QValueAxis) + QChartView(parent), last_step(-1), index(_index), window(10), order(4), chart(new QChart), + series(new QLineSeries), smooth(nullptr), xaxis(new QValueAxis), yaxis(new QValueAxis), + do_smooth(false) { chart->legend()->hide(); chart->addAxis(xaxis, Qt::AlignBottom); @@ -343,6 +383,17 @@ ChartViewer::ChartViewer(const QString &title, int _index, QWidget *parent) : /* -------------------------------------------------------------------- */ +ChartViewer::~ChartViewer() +{ + delete xaxis; + delete yaxis; + delete smooth; + delete series; + delete chart; +} + +/* -------------------------------------------------------------------- */ + void ChartViewer::add_data(int step, double data) { if (last_step < step) { @@ -375,6 +426,7 @@ void ChartViewer::add_data(int step, double data) if (last_update.msecsTo(QTime::currentTime()) > settings.value("updchart", "500").toInt()) { last_update = QTime::currentTime(); reset_zoom(); + update_smooth(); } } } @@ -385,6 +437,7 @@ void ChartViewer::reset_zoom() { auto points = series->points(); + // get min/max for plot qreal xmin = 1.0e100; qreal xmax = -1.0e100; qreal ymin = 1.0e100; @@ -396,6 +449,17 @@ void ChartViewer::reset_zoom() ymax = qMax(ymax, p.y()); } + // if plotting the smoothed plot, check for its min/max values, too + if (smooth) { + auto spoints = smooth->points(); + for (auto &p : spoints) { + xmin = qMin(xmin, p.x()); + xmax = qMax(xmax, p.x()); + ymin = qMin(ymin, p.y()); + ymax = qMax(ymax, p.y()); + } + } + // avoid (nearly) empty ranges double deltax = xmax - xmin; if ((deltax / ((xmax == 0.0) ? 1.0 : xmax)) < 1.0e-10) { @@ -423,6 +487,490 @@ void ChartViewer::reset_zoom() yaxis->setRange(ymin, ymax); } +/* -------------------------------------------------------------------- */ + +void ChartViewer::smooth_param(bool _do_smooth, int _window, int _order) +{ + // turn off smooth plot + if (!_do_smooth) { + if (smooth) { + chart->removeSeries(smooth); + delete smooth; + smooth = nullptr; + } + } + do_smooth = _do_smooth; + window = _window; + order = _order; + update_smooth(); +} + +/* -------------------------------------------------------------------- */ + +// update smooth plot data + +static QList calc_sgsmooth(const QList &input, const int window, const int order); + +void ChartViewer::update_smooth() +{ + if (do_smooth) { + if (series->count() > (2 * window)) { + if (!smooth) { + smooth = new QLineSeries; + chart->addSeries(smooth); + smooth->attachAxis(xaxis); + smooth->attachAxis(yaxis); + } + smooth->clear(); + smooth->append(calc_sgsmooth(series->points(), window, order)); + } + } +} + +//! default convergence +static constexpr double TINY_FLOAT = 1.0e-300; + +//! comfortable array of doubles +typedef std::vector float_vect; +//! comfortable array of ints; +typedef std::vector int_vect; + +// savitzky golay smoothing. +static float_vect sg_smooth(const float_vect &v, const int w, const int deg); + +QList calc_sgsmooth(const QList &input, int window, int order) +{ + const int ndat = input.count(); + if (ndat < 2 * window + 2) window = ndat / 2 - 1; + + if (window > 1) { + float_vect in(ndat); + QList rv; + + for (int i = 0; i < ndat; ++i) { + in[i] = input[i].y(); + } + float_vect out = sg_smooth(in, window, order); + + for (int i = 0; i < ndat; ++i) { + rv.append(QPointF(input[i].x(), out[i])); + } + return rv; + } else { + return input; + } +} + +/*! matrix class. + * + * This is a matrix class derived from a vector of float_vects. Note that + * the matrix elements indexed [row][column] with indices starting at 0 (c + * style). Also note that because of its design looping through rows should + * be faster than looping through columns. + * + * \brief two dimensional floating point array + */ +class float_mat : public std::vector { +private: + //! disable the default constructor + explicit float_mat() {}; + //! disable assignment operator until it is implemented. + float_mat &operator=(const float_mat &) { return *this; }; + +public: + //! constructor with sizes + float_mat(const std::size_t rows, const std::size_t cols, const double def = 0.0); + //! copy constructor for matrix + float_mat(const float_mat &m); + //! copy constructor for vector + float_mat(const float_vect &v); + + //! use default destructor + // ~float_mat() {}; + + //! get size + int nr_rows(void) const { return size(); }; + //! get size + int nr_cols(void) const { return front().size(); }; +}; + +// constructor with sizes +float_mat::float_mat(const std::size_t rows, const std::size_t cols, const double defval) : + std::vector(rows) +{ + for (std::size_t i = 0; i < rows; ++i) { + (*this)[i].resize(cols, defval); + } +#if 0 + if ((rows < 1) || (cols < 1)) { + char buffer[1024]; + + sprintf(buffer, "cannot build matrix with %d rows and %d columns\n", + rows, cols); + sgs_error(buffer); + } +#endif +} + +// copy constructor for matrix +float_mat::float_mat(const float_mat &m) : std::vector(m.size()) +{ + + float_mat::iterator inew = begin(); + float_mat::const_iterator iold = m.begin(); + for (/* empty */; iold < m.end(); ++inew, ++iold) { + const auto oldsz = iold->size(); + inew->resize(oldsz); + const float_vect oldvec(*iold); + *inew = oldvec; + } +} + +// copy constructor for vector +float_mat::float_mat(const float_vect &v) : std::vector(1) +{ + + const auto oldsz = v.size(); + front().resize(oldsz); + front() = v; +} + +////////////////////// +// Helper functions // +////////////////////// + +//! permute() orders the rows of A to match the integers in the index array. +void permute(float_mat &A, int_vect &idx) +{ + int_vect i(idx.size()); + + for (int j = 0; j < A.nr_rows(); ++j) { + i[j] = j; + } + + // loop over permuted indices + for (int j = 0; j < A.nr_rows(); ++j) { + if (i[j] != idx[j]) { + + // search only the remaining indices + for (int k = j + 1; k < A.nr_rows(); ++k) { + if (i[k] == idx[j]) { + std::swap(A[j], A[k]); // swap the rows and + i[k] = i[j]; // the elements of + i[j] = idx[j]; // the ordered index. + break; // next j + } + } + } + } +} + +/*! \brief Implicit partial pivoting. + * + * The function looks for pivot element only in rows below the current + * element, A[idx[row]][column], then swaps that row with the current one in + * the index map. The algorithm is for implicit pivoting (i.e., the pivot is + * chosen as if the max coefficient in each row is set to 1) based on the + * scaling information in the vector scale. The map of swapped indices is + * recorded in swp. The return value is +1 or -1 depending on whether the + * number of row swaps was even or odd respectively. */ +static int partial_pivot(float_mat &A, const std::size_t row, const std::size_t col, + float_vect &scale, int_vect &idx, double tol) +{ + if (tol <= 0.0) tol = TINY_FLOAT; + + int swapNum = 1; + + // default pivot is the current position, [row,col] + std::size_t pivot = row; + double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]]; + + // loop over possible pivots below current + for (int j = row + 1; j < A.nr_rows(); ++j) { + + const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]]; + + // if this elem is larger, then it becomes the pivot + if (tmp > piv_elem) { + pivot = j; + piv_elem = tmp; + } + } + +#if 0 + if(piv_elem < tol) { + sgs_error("partial_pivot(): Zero pivot encountered.\n") +#endif + + if (pivot > row) { // bring the pivot to the diagonal + int j = idx[row]; // reorder swap array + idx[row] = idx[pivot]; + idx[pivot] = j; + swapNum = -swapNum; // keeping track of odd or even swap + } + return swapNum; +} + +/*! \brief Perform backward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is upper + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the lower triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_backsubst(float_mat &A, float_mat &a, bool diag = false) +{ + for (int r = (A.nr_rows() - 1); r >= 0; --r) { + for (int c = (A.nr_cols() - 1); c > r; --c) { + for (int k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if (!diag) { + for (int k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Perform forward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is lower + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the upper triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_forwsubst(float_mat &A, float_mat &a, bool diag = true) +{ + for (int r = 0; r < A.nr_rows(); ++r) { + for (int c = 0; c < r; ++c) { + for (int k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if (!diag) { + for (int k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Performs LU factorization in place. + * + * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of + * swapped indeces is recorded in idx. The return value is +1 or -1 + * depending on whether the number of row swaps was even or odd + * respectively. idx must be preinitialized to a valid set of indices + * (e.g., {1,2, ... ,A.nr_rows()}). */ +static int lu_factorize(float_mat &A, int_vect &idx, double tol = TINY_FLOAT) +{ + if (tol <= 0.0) tol = TINY_FLOAT; +#if 0 + if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) { + sgs_error("lu_factorize(): cannot handle empty " + "or nonsquare matrices.\n"); + + return 0; + } +#endif + float_vect scale(A.nr_rows()); // implicit pivot scaling + for (int i = 0; i < A.nr_rows(); ++i) { + double maxval = 0.0; + for (int j = 0; j < A.nr_cols(); ++j) { + if (fabs(A[i][j]) > maxval) maxval = fabs(A[i][j]); + } + if (maxval == 0.0) { +#if 0 + sgs_error("lu_factorize(): zero pivot found.\n"); +#endif + return 0; + } + scale[i] = 1.0 / maxval; + } + + int swapNum = 1; + for (int c = 0; c < A.nr_cols(); ++c) { // loop over columns + swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal + for (int r = 0; r < A.nr_rows(); ++r) { // loop over rows + int lim = (r < c) ? r : c; + for (int j = 0; j < lim; ++j) { + A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c]; + } + if (r > c) A[idx[r]][c] /= A[idx[c]][c]; + } + } + permute(A, idx); + return swapNum; +} + +/*! \brief Solve a system of linear equations. + * Solves the inhomogeneous matrix problem with lu-decomposition. Note that + * inversion may be accomplished by setting a to the identity_matrix. */ +static float_mat lin_solve(const float_mat &A, const float_mat &a, double tol = TINY_FLOAT) +{ + float_mat B(A); + float_mat b(a); + int_vect idx(B.nr_rows()); + + for (int j = 0; j < B.nr_rows(); ++j) { + idx[j] = j; // init row swap label array + } + lu_factorize(B, idx, tol); // get the lu-decomp. + permute(b, idx); // sort the inhomogeneity to match the lu-decomp + lu_forwsubst(B, b); // solve the forward problem + lu_backsubst(B, b); // solve the backward problem + return b; +} + +/////////////////////// +// related functions // +/////////////////////// + +//! Returns the inverse of a matrix using LU-decomposition. +static float_mat invert(const float_mat &A) +{ + const int n = A.size(); + float_mat E(n, n, 0.0); + float_mat B(A); + + for (int i = 0; i < n; ++i) { + E[i][i] = 1.0; + } + + return lin_solve(B, E); +} + +//! returns the transposed matrix. +static float_mat transpose(const float_mat &a) +{ + float_mat res(a.nr_cols(), a.nr_rows()); + + for (int i = 0; i < a.nr_rows(); ++i) { + for (int j = 0; j < a.nr_cols(); ++j) { + res[j][i] = a[i][j]; + } + } + return res; +} + +//! matrix multiplication. +float_mat operator*(const float_mat &a, const float_mat &b) +{ + float_mat res(a.nr_rows(), b.nr_cols()); +#if 0 + if (a.nr_cols() != b.nr_rows()) { + sgs_error("incompatible matrices in multiplication\n"); + return res; + } +#endif + for (int i = 0; i < a.nr_rows(); ++i) { + for (int j = 0; j < b.nr_cols(); ++j) { + double sum(0.0); + for (int k = 0; k < a.nr_cols(); ++k) { + sum += a[i][k] * b[k][j]; + } + res[i][j] = sum; + } + } + return res; +} + +//! calculate savitzky golay coefficients. +static float_vect sg_coeff(const float_vect &b, const std::size_t deg) +{ + const std::size_t rows(b.size()); + const std::size_t cols(deg + 1); + float_mat A(rows, cols); + float_vect res(rows); + + // generate input matrix for least squares fit + for (std::size_t i = 0; i < rows; ++i) { + for (std::size_t j = 0; j < cols; ++j) { + A[i][j] = pow(double(i), double(j)); + } + } + + float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); + + for (std::size_t i = 0; i < b.size(); ++i) { + res[i] = c[0][0]; + for (std::size_t j = 1; j <= deg; ++j) { + res[i] += c[j][0] * pow(double(i), double(j)); + } + } + return res; +} + +/*! \brief savitzky golay smoothing. + * + * This method means fitting a polynome of degree 'deg' to a sliding window + * of width 2w+1 throughout the data. The needed coefficients are + * generated dynamically by doing a least squares fit on a "symmetric" unit + * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome + * yields the sg-coefficients. at the border non symmectric vectors b are + * used. */ +float_vect sg_smooth(const float_vect &v, const int width, const int deg) +{ + float_vect res(v.size(), 0.0); +#if 0 + if ((width < 1) || (deg < 0) || (v.size() < (2 * width + 2))) { + sgs_error("sgsmooth: parameter error.\n"); + return res; + } +#endif + const int window = 2 * width + 1; + const int endidx = v.size() - 1; + + // do a regular sliding window average + if (deg == 0) { + // handle border cases first because we need different coefficients + for (int i = 0; i < width; ++i) { + const double scale = 1.0 / double(i + 1); + const float_vect c1(width, scale); + for (int j = 0; j <= i; ++j) { + res[i] += c1[j] * v[j]; + res[endidx - i] += c1[j] * v[endidx - j]; + } + } + + // now loop over rest of data. reusing the "symmetric" coefficients. + const double scale = 1.0 / double(window); + const float_vect c2(window, scale); + for (std::size_t i = 0; i <= (v.size() - window); ++i) { + for (int j = 0; j < window; ++j) { + res[i + width] += c2[j] * v[i + j]; + } + } + return res; + } + + // handle border cases first because we need different coefficients + for (int i = 0; i < width; ++i) { + float_vect b1(window, 0.0); + b1[i] = 1.0; + + const float_vect c1(sg_coeff(b1, deg)); + for (int j = 0; j < window; ++j) { + res[i] += c1[j] * v[j]; + res[endidx - i] += c1[j] * v[endidx - j]; + } + } + + // now loop over rest of data. reusing the "symmetric" coefficients. + float_vect b2(window, 0.0); + b2[width] = 1.0; + const float_vect c2(sg_coeff(b2, deg)); + + for (std::size_t i = 0; i <= (v.size() - window); ++i) { + for (int j = 0; j < window; ++j) { + res[i + width] += c2[j] * v[i + j]; + } + } + return res; +} + // Local Variables: // c-basic-offset: 4 // End: diff --git a/tools/lammps-gui/chartviewer.h b/tools/lammps-gui/chartviewer.h index ca1c33979f..58f65ba215 100644 --- a/tools/lammps-gui/chartviewer.h +++ b/tools/lammps-gui/chartviewer.h @@ -17,14 +17,16 @@ #include #include #include -#include #include +#include class QAction; class QCloseEvent; class QEvent; class QMenuBar; class QMenu; +class QCheckBox; +class QSpinBox; namespace QtCharts { class ChartViewer; } @@ -49,6 +51,7 @@ private slots: void quit(); void reset_zoom(); void stop_run(); + void update_smooth(); void saveAs(); void exportDat(); @@ -67,6 +70,8 @@ private: QComboBox *columns; QAction *saveAsAct, *exportCsvAct, *exportDatAct, *exportYamlAct; QAction *closeAct, *stopAct, *quitAct; + QCheckBox *smooth; + QSpinBox *window, *order; QString filename; QList charts; @@ -85,9 +90,12 @@ class ChartViewer : public QChartView { public: explicit ChartViewer(const QString &title, int index, QWidget *parent = nullptr); + ~ChartViewer(); void add_data(int step, double data); void reset_zoom(); + void smooth_param(bool _do_smooth, int _window, int _order); + void update_smooth(); int get_index() const { return index; }; int get_count() const { return series->count(); } @@ -97,11 +105,13 @@ public: private: int last_step, index; + int window, order; QChart *chart; - QLineSeries *series; + QLineSeries *series, *smooth; QValueAxis *xaxis; QValueAxis *yaxis; QTime last_update; + bool do_smooth; }; } // namespace QtCharts #endif diff --git a/tools/lammps-gui/lammpsgui.qrc b/tools/lammps-gui/lammpsgui.qrc index 5de1a5efba..cf9dd20dda 100644 --- a/tools/lammps-gui/lammpsgui.qrc +++ b/tools/lammps-gui/lammpsgui.qrc @@ -13,6 +13,7 @@ icons/application-plot.png icons/application-yaml.png icons/axes-img.png + icons/chart-smooth.png icons/document-new.png icons/document-open-recent.png icons/document-open.png diff --git a/tools/lammps-gui/preferences.cpp b/tools/lammps-gui/preferences.cpp index bce4ed46b0..2955d37eb1 100644 --- a/tools/lammps-gui/preferences.cpp +++ b/tools/lammps-gui/preferences.cpp @@ -260,8 +260,8 @@ GeneralTab::GeneralTab(QSettings *_settings, LammpsWrapper *_lammps, QWidget *pa connect(getallfont, &QPushButton::released, this, &GeneralTab::newallfont); connect(gettextfont, &QPushButton::released, this, &GeneralTab::newtextfont); - auto *freqlabel = new QLabel("Data update interval (ms)"); - auto *freqval = new QSpinBox; + auto *freqlabel = new QLabel("Data update interval (ms)"); + auto *freqval = new QSpinBox; freqval->setRange(1, 1000); freqval->setStepType(QAbstractSpinBox::AdaptiveDecimalStepType); freqval->setValue(settings->value("updfreq", "10").toInt());