ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add smoothed plot to chart windows

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-08-18 00:53:10 -04:00
parent e4b9da147a
commit 98a7e045bb
4 changed files with 565 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

552
tools/lammps-gui/chartviewer.cpp
View File

@ -18,6 +18,7 @@
#include <QAction>
#include <QApplication>
#include <QChart>
#include <QCheckBox>
#include <QCloseEvent>
#include <QComboBox>
#include <QEvent>
@ -32,6 +33,7 @@
#include <QPushButton>
#include <QSettings>
#include <QSpacerItem>
#include <QSpinBox>
#include <QTextStream>
#include <QTime>
#include <QVBoxLayout>
@ -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<QPointF> calc_sgsmooth(const QList<QPointF> &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<double> float_vect;
//! comfortable array of ints;
typedef std::vector<int> int_vect;
// savitzky golay smoothing.
static float_vect sg_smooth(const float_vect &v, const int w, const int deg);
QList<QPointF> calc_sgsmooth(const QList<QPointF> &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<QPointF> 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<float_vect> {
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<float_vect>(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<float_vect>(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<float_vect>(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:

14
tools/lammps-gui/chartviewer.h
View File

@ -17,14 +17,16 @@
#include <QComboBox>
#include <QList>
#include <QString>
#include <QWidget>
#include <QTime>
#include <QWidget>
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<QtCharts::ChartViewer *> 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

1
tools/lammps-gui/lammpsgui.qrc
View File

@ -13,6 +13,7 @@
<file>icons/application-plot.png</file>
<file>icons/application-yaml.png</file>
<file>icons/axes-img.png</file>
<file>icons/chart-smooth.png</file>
<file>icons/document-new.png</file>
<file>icons/document-open-recent.png</file>
<file>icons/document-open.png</file>