leonarski_fandClaude Opus 5
8493095912
The error model says sigma -> b*I for strong reflections, so merged I/sigma flattens off at 1/b - the number reported as ISa. Plotting I/sigma against I with that asymptote drawn on it is what shows whether the reported ISa describes the data or comes from a degenerate fit, which nothing in the window could show before. A third page next to the per-shell plot and table. The merge carries a few thousand strided (I, sigma) pairs to the viewer for it - a shape, not a reflection list; the reflections themselves are in the .mtz/.cif. The hero row gains the Wilson B and the radiation-damage Delta-B. Both were already computed and already in MergeStatistics, so they only needed showing. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
389 lines
18 KiB
C++
389 lines
18 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||
// SPDX-License-Identifier: GPL-3.0-only
|
||
|
||
#include "JFJochMergeStatsWindow.h"
|
||
|
||
#include <cmath>
|
||
#include <limits>
|
||
|
||
#include <QButtonGroup>
|
||
#include <QComboBox>
|
||
#include <QFrame>
|
||
#include <QHBoxLayout>
|
||
#include <QHeaderView>
|
||
#include <QLabel>
|
||
#include <QStackedWidget>
|
||
#include <QStringList>
|
||
#include <QTableWidget>
|
||
#include <QToolButton>
|
||
#include <QVBoxLayout>
|
||
#include <QtCharts/QChart>
|
||
#include <QtCharts/QChartView>
|
||
#include <QtCharts/QLineSeries>
|
||
#include <QtCharts/QLogValueAxis>
|
||
#include <QtCharts/QScatterSeries>
|
||
#include <QtCharts/QValueAxis>
|
||
|
||
#include "gemmi/symmetry.hpp" // find_spacegroup_by_number
|
||
|
||
#include "../charts/JFJochSimpleChartView.h"
|
||
#include "../widgets/ToolbarIcons.h"
|
||
|
||
namespace {
|
||
// Per-shell metric value (NaN -> dropped by the chart). Keep in sync with the combo entries.
|
||
double ShellMetric(const MergeStatisticsShell &s, int metric) {
|
||
switch (metric) {
|
||
case 0: return s.cc_half * 100.0; // CC1/2 (%)
|
||
case 1: return s.r_meas * 100.0; // R-meas (%)
|
||
case 2: return s.mean_i_over_sigma; // <I/sigma>
|
||
case 3: return s.possible_unique_reflections > 0 // Completeness (%)
|
||
? 100.0 * s.unique_reflections / s.possible_unique_reflections : NAN;
|
||
case 4: return s.unique_reflections > 0 // Multiplicity
|
||
? static_cast<double>(s.total_observations) / s.unique_reflections : NAN;
|
||
case 5: return s.cc_ref * 100.0; // CCref (%)
|
||
default: return NAN;
|
||
}
|
||
}
|
||
|
||
// A big-number "card" for the hero row: large navy value over a small grey caption.
|
||
QWidget *MakeCard(const QString &value, const QString &caption, QWidget *parent) {
|
||
auto *card = new QFrame(parent);
|
||
card->setFrameShape(QFrame::StyledPanel);
|
||
auto *l = new QVBoxLayout(card);
|
||
l->setContentsMargins(12, 8, 12, 8);
|
||
l->setSpacing(0);
|
||
auto *v = new QLabel(value, card);
|
||
QFont f = v->font();
|
||
f.setPointSizeF(f.pointSizeF() * 2.2);
|
||
f.setBold(true);
|
||
v->setFont(f);
|
||
v->setStyleSheet("color: #1F3A5F;");
|
||
v->setAlignment(Qt::AlignCenter);
|
||
auto *c = new QLabel(caption, card);
|
||
c->setStyleSheet("color: gray;");
|
||
c->setAlignment(Qt::AlignCenter);
|
||
l->addWidget(v);
|
||
l->addWidget(c);
|
||
return card;
|
||
}
|
||
}
|
||
|
||
JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats,
|
||
double isa, bool has_reference,
|
||
const TwinningAnalysisResult &twinning,
|
||
std::optional<int64_t> space_group_number,
|
||
const std::optional<SearchSpaceGroupResult> &space_group_search,
|
||
const std::vector<std::pair<float, float>> &merged_i_sigma,
|
||
QWidget *parent)
|
||
: QWidget(parent, Qt::Window), stats_(stats), has_reference_(has_reference) {
|
||
setWindowTitle("Merge statistics — " + title);
|
||
setAttribute(Qt::WA_DeleteOnClose);
|
||
resize(740, 560);
|
||
|
||
auto *layout = new QVBoxLayout(this);
|
||
|
||
// --- Hero row: overall numbers ---
|
||
const auto &o = stats_.overall;
|
||
const double completeness = o.possible_unique_reflections > 0
|
||
? 100.0 * o.unique_reflections / o.possible_unique_reflections : NAN;
|
||
const double multiplicity = o.unique_reflections > 0
|
||
? static_cast<double>(o.total_observations) / o.unique_reflections : NAN;
|
||
auto pct = [](double v) { return std::isfinite(v) ? QString::number(v, 'f', 1) + "%" : QStringLiteral("—"); };
|
||
auto num = [](double v, int d) { return std::isfinite(v) ? QString::number(v, 'f', d) : QStringLiteral("—"); };
|
||
|
||
auto *hero = new QHBoxLayout();
|
||
// High-resolution limit of the merge (d_min). When no scaling resolution limit is given it is set
|
||
// by the CC1/2 logistic cutoff, so it varies per run - show it first.
|
||
if (o.d_min < o.d_max)
|
||
hero->addWidget(MakeCard(QString::asprintf("%.2f", o.d_min), "Resolution [Å]", this));
|
||
hero->addWidget(MakeCard(num(isa, 1), "ISa", this));
|
||
hero->addWidget(MakeCard(pct(o.cc_half * 100.0), "CC½", this));
|
||
hero->addWidget(MakeCard(pct(o.r_meas * 100.0), "R-meas", this));
|
||
hero->addWidget(MakeCard(pct(completeness), "Completeness", this));
|
||
hero->addWidget(MakeCard(num(multiplicity, 1), "Multiplicity", this));
|
||
// Wilson B describes how fast the intensities fall off with resolution; the radiation-damage
|
||
// Delta-B is how much MORE they fall off by the end of the run than at the start (rotation only).
|
||
if (std::isfinite(stats_.wilson_b))
|
||
hero->addWidget(MakeCard(num(stats_.wilson_b, 1), "Wilson B [Ų]", this));
|
||
if (std::isfinite(stats_.radiation_damage_delta_b)) {
|
||
const double db = stats_.radiation_damage_delta_b;
|
||
hero->addWidget(MakeCard((db > 0 ? "+" : "") + num(db, 1), "ΔB, damage [Ų]", this));
|
||
}
|
||
if (has_reference_)
|
||
hero->addWidget(MakeCard(pct(o.cc_ref * 100.0), "CCref", this));
|
||
if (space_group_number.has_value()) {
|
||
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(static_cast<int>(*space_group_number));
|
||
const QString sg_name = sg ? QString::fromStdString(sg->short_name())
|
||
: QString::number(*space_group_number);
|
||
hero->addWidget(MakeCard(sg_name, "Space group", this));
|
||
}
|
||
layout->addLayout(hero);
|
||
|
||
// De-novo space-group search: the point group + the ranked candidates it scored (only present when
|
||
// the space group was not fixed by the user). A summary line above a compact table.
|
||
if (space_group_search.has_value() && !space_group_search->candidates.empty()) {
|
||
const auto &sgs = *space_group_search;
|
||
QString summary = "Point group " + QString::fromStdString(sgs.point_group_hm.empty() ? "?" : sgs.point_group_hm);
|
||
if (!sgs.alternatives.empty()) {
|
||
QStringList alts;
|
||
for (const auto &alt : sgs.alternatives)
|
||
alts << QString::fromStdString(alt.short_name());
|
||
summary += " · indistinguishable from these data: " + alts.join(", ");
|
||
}
|
||
auto *sgLabel = new QLabel(summary, this);
|
||
sgLabel->setStyleSheet("color: gray;");
|
||
layout->addWidget(sgLabel);
|
||
|
||
auto *sgTable = new QTableWidget(this);
|
||
sgTable->setEditTriggers(QAbstractItemView::NoEditTriggers);
|
||
sgTable->verticalHeader()->setVisible(false);
|
||
const QStringList sgHeaders{"Space group", "Absent", "Viol.", "⟨I/σ⟩ abs", "⟨I/σ⟩ pres", "Consistent"};
|
||
sgTable->setColumnCount(sgHeaders.size());
|
||
sgTable->setHorizontalHeaderLabels(sgHeaders);
|
||
sgTable->setRowCount(static_cast<int>(sgs.candidates.size()));
|
||
for (int i = 0; i < static_cast<int>(sgs.candidates.size()); ++i) {
|
||
const auto &c = sgs.candidates[i];
|
||
int col = 0;
|
||
auto set = [&](const QString &t, bool leftAlign = false) {
|
||
auto *it = new QTableWidgetItem(t);
|
||
it->setTextAlignment((leftAlign ? Qt::AlignLeft : Qt::AlignRight) | Qt::AlignVCenter);
|
||
if (c.selected) { QFont f = it->font(); f.setBold(true); it->setFont(f); }
|
||
sgTable->setItem(i, col++, it);
|
||
};
|
||
set(QString::fromStdString(c.space_group.short_name()) + (c.selected ? " *" : ""), true);
|
||
set(QString::number(c.absent_observed));
|
||
set(QString::number(c.absent_violations));
|
||
set(QString::number(c.absent_mean_i_over_sigma, 'f', 2));
|
||
set(QString::number(c.present_mean_i_over_sigma, 'f', 2));
|
||
set(c.consistent ? QStringLiteral("yes") : QStringLiteral("no"));
|
||
}
|
||
sgTable->horizontalHeader()->setSectionResizeMode(QHeaderView::Stretch);
|
||
sgTable->setMaximumHeight(160); // a secondary detail — keep the shells plot the focus
|
||
layout->addWidget(sgTable);
|
||
}
|
||
|
||
// --- Twinning test (Padilla-Yeates L-test + second moment) ---
|
||
// A red-outlined banner flags suspected twinning; otherwise a quiet grey "all clear".
|
||
if (twinning.l_test_pairs > 0) {
|
||
QString text = QString("Twinning: ⟨|L|⟩ = %1 ⟨I²⟩/⟨I⟩² = %2")
|
||
.arg(twinning.mean_abs_l, 0, 'f', 3)
|
||
.arg(twinning.second_moment, 0, 'f', 2);
|
||
if (twinning.twinning_suspected)
|
||
text += QString(" — twinning suspected (estimated fraction ≈ %1)")
|
||
.arg(twinning.estimated_twin_fraction, 0, 'f', 2);
|
||
else
|
||
text += " — no twinning indicated";
|
||
|
||
auto *banner = new QLabel(text, this);
|
||
banner->setContentsMargins(10, 6, 10, 6);
|
||
if (twinning.twinning_suspected)
|
||
banner->setStyleSheet("color: rgb(160,0,0); border: 1px solid red; border-radius: 4px;");
|
||
else
|
||
banner->setStyleSheet("color: gray;");
|
||
layout->addWidget(banner);
|
||
}
|
||
|
||
// --- Controls: plot/table view toggle (icon buttons) + per-resolution metric selector ---
|
||
auto *row = new QHBoxLayout();
|
||
auto *plotBtn = new QToolButton(this);
|
||
plotBtn->setIcon(ToolbarIcons::linePlot());
|
||
plotBtn->setCheckable(true);
|
||
plotBtn->setChecked(true);
|
||
plotBtn->setToolTip("Per-resolution plot");
|
||
plotBtn->setStyleSheet(ToolbarIcons::buttonStyle());
|
||
auto *tableBtn = new QToolButton(this);
|
||
tableBtn->setIcon(ToolbarIcons::table());
|
||
tableBtn->setCheckable(true);
|
||
tableBtn->setToolTip("Per-shell table");
|
||
tableBtn->setStyleSheet(ToolbarIcons::buttonStyle());
|
||
auto *isaBtn = new QToolButton(this);
|
||
isaBtn->setIcon(ToolbarIcons::linePlot());
|
||
isaBtn->setCheckable(true);
|
||
isaBtn->setToolTip("ISa diagnostic: merged I/σ against I");
|
||
isaBtn->setStyleSheet(ToolbarIcons::buttonStyle());
|
||
auto *viewGroup = new QButtonGroup(this);
|
||
viewGroup->setExclusive(true);
|
||
viewGroup->addButton(plotBtn, 0);
|
||
viewGroup->addButton(tableBtn, 1);
|
||
viewGroup->addButton(isaBtn, 2);
|
||
row->addWidget(plotBtn);
|
||
row->addWidget(tableBtn);
|
||
row->addWidget(isaBtn);
|
||
row->addSpacing(12);
|
||
metricLabel_ = new QLabel("Metric:", this);
|
||
row->addWidget(metricLabel_);
|
||
metric_ = new QComboBox(this);
|
||
metric_->addItem("CC½ (%)", 0);
|
||
metric_->addItem("R-meas (%)", 1);
|
||
metric_->addItem("⟨I/σ⟩", 2);
|
||
metric_->addItem("Completeness (%)", 3);
|
||
metric_->addItem("Multiplicity", 4);
|
||
if (has_reference_)
|
||
metric_->addItem("CCref (%)", 5);
|
||
row->addWidget(metric_);
|
||
row->addStretch();
|
||
layout->addLayout(row);
|
||
|
||
// --- Stacked plot / table ---
|
||
stack_ = new QStackedWidget(this);
|
||
chart_ = new JFJochSimpleChartView(this);
|
||
table_ = new QTableWidget(this);
|
||
table_->setEditTriggers(QAbstractItemView::NoEditTriggers);
|
||
table_->verticalHeader()->setVisible(false);
|
||
stack_->addWidget(chart_); // page 0
|
||
stack_->addWidget(table_); // page 1
|
||
isa_page_ = buildIsaDiagnostic(merged_i_sigma, isa);
|
||
if (isa_page_)
|
||
stack_->addWidget(isa_page_); // page 2
|
||
else
|
||
isaBtn->setEnabled(false);
|
||
layout->addWidget(stack_, 1);
|
||
|
||
connect(metric_, &QComboBox::currentIndexChanged, this, [this](int) { updatePlot(); });
|
||
connect(viewGroup, &QButtonGroup::idClicked, this, [this](int idx) {
|
||
stack_->setCurrentIndex(idx);
|
||
metric_->setVisible(idx == 0); // the metric selector + its label only apply to the plot
|
||
metricLabel_->setVisible(idx == 0);
|
||
});
|
||
|
||
buildTable();
|
||
updatePlot();
|
||
}
|
||
|
||
// Merged I/sigma against I. The error model says sigma -> b*I for strong reflections, so I/sigma
|
||
// flattens off at 1/b = ISa: the plot shows whether the reported number is where the data actually
|
||
// level out. Log intensity axis, because the reflections span several decades and the asymptote lives
|
||
// at the strong end.
|
||
QWidget *JFJochMergeStatsWindow::buildIsaDiagnostic(const std::vector<std::pair<float, float>> &i_sigma,
|
||
double isa) {
|
||
auto *points = new QScatterSeries();
|
||
points->setName("merged reflections");
|
||
points->setMarkerSize(4.0);
|
||
points->setBorderColor(Qt::transparent);
|
||
points->setColor(QColor(0x1F, 0x3A, 0x5F, 0x60));
|
||
|
||
double i_lo = std::numeric_limits<double>::max(), i_hi = 0.0, r_hi = 0.0;
|
||
for (const auto &[I, sigma] : i_sigma) {
|
||
if (!(I > 0.0f) || !(sigma > 0.0f))
|
||
continue; // the log axis cannot show a non-positive intensity
|
||
const double ratio = I / sigma;
|
||
points->append(I, ratio);
|
||
i_lo = std::min(i_lo, static_cast<double>(I));
|
||
i_hi = std::max(i_hi, static_cast<double>(I));
|
||
r_hi = std::max(r_hi, ratio);
|
||
}
|
||
if (points->count() < 2 || !(i_hi > i_lo)) {
|
||
delete points;
|
||
return nullptr;
|
||
}
|
||
|
||
auto *chart = new QChart();
|
||
chart->legend()->setVisible(false);
|
||
chart->addSeries(points);
|
||
|
||
auto *ax = new QLogValueAxis();
|
||
ax->setBase(10.0);
|
||
ax->setTitleText("merged I");
|
||
ax->setRange(i_lo, i_hi);
|
||
auto *ay = new QValueAxis();
|
||
ay->setTitleText("merged I/σ");
|
||
ay->setRange(0.0, std::max(r_hi, std::isfinite(isa) ? isa * 1.2 : 0.0) * 1.05);
|
||
chart->addAxis(ax, Qt::AlignBottom);
|
||
chart->addAxis(ay, Qt::AlignLeft);
|
||
points->attachAxis(ax);
|
||
points->attachAxis(ay);
|
||
|
||
if (std::isfinite(isa) && isa > 0.0) {
|
||
auto *line = new QLineSeries();
|
||
line->append(i_lo, isa);
|
||
line->append(i_hi, isa);
|
||
QPen pen(QColor(0xC0, 0x39, 0x2B));
|
||
pen.setStyle(Qt::DashLine);
|
||
pen.setWidth(2);
|
||
line->setPen(pen);
|
||
chart->addSeries(line);
|
||
line->attachAxis(ax);
|
||
line->attachAxis(ay);
|
||
}
|
||
|
||
auto *view = new QChartView(chart, this);
|
||
view->setRenderHint(QPainter::Antialiasing);
|
||
view->setToolTip(std::isfinite(isa) && isa > 0.0
|
||
? QString("Dashed line: reported ISa = %1. The strong reflections should level "
|
||
"off there.").arg(isa, 0, 'f', 1)
|
||
: QString("No ISa was determined for this merge."));
|
||
return view;
|
||
}
|
||
|
||
void JFJochMergeStatsWindow::updatePlot() {
|
||
const int metric = metric_->currentData().toInt();
|
||
std::vector<float> x, y;
|
||
double dmax = 0.0;
|
||
x.reserve(stats_.shells.size());
|
||
y.reserve(stats_.shells.size());
|
||
for (const auto &s: stats_.shells) {
|
||
x.push_back(s.mean_one_over_d2);
|
||
const double v = ShellMetric(s, metric);
|
||
y.push_back(static_cast<float>(v));
|
||
if (std::isfinite(v))
|
||
dmax = std::max(dmax, v);
|
||
}
|
||
chart_->UpdateData(x, y, "Resolution [Å]", metric_->currentText(), true);
|
||
|
||
// Absolute y-axis so positions are comparable: CC1/2 + CCref are 0..100, completeness 0..(>=100),
|
||
// multiplicity / R-meas / I-sigma start at 0. Show the labels (UpdateData hides them by default).
|
||
double ymin = 0.0, ymax = 1.0;
|
||
switch (metric) {
|
||
case 0: case 5: ymax = 100.0; break; // CC1/2, CCref
|
||
case 3: ymax = std::max(100.0, dmax * 1.05); break; // Completeness
|
||
default: ymax = dmax > 0.0 ? dmax * 1.05 : 1.0; break; // R-meas, I/sigma, multiplicity
|
||
}
|
||
const auto axes = chart_->chart()->axes(Qt::Vertical);
|
||
if (!axes.isEmpty()) {
|
||
if (auto *vax = qobject_cast<QValueAxis *>(axes.first())) {
|
||
vax->setRange(ymin, ymax);
|
||
vax->setLabelsVisible(true);
|
||
}
|
||
}
|
||
}
|
||
|
||
void JFJochMergeStatsWindow::buildTable() {
|
||
QStringList headers{"Res. [Å]", "N_obs", "Unique", "Compl. %", "Mult.", "⟨I/σ⟩", "R-meas %", "CC½ %"};
|
||
if (has_reference_)
|
||
headers << "CCref %";
|
||
table_->setColumnCount(headers.size());
|
||
table_->setHorizontalHeaderLabels(headers);
|
||
table_->setRowCount(static_cast<int>(stats_.shells.size()) + 1); // shells + overall
|
||
|
||
auto fillRow = [&](int rowIdx, const MergeStatisticsShell &s, const QString &resLabel) {
|
||
int c = 0;
|
||
auto set = [&](const QString &t) {
|
||
auto *it = new QTableWidgetItem(t);
|
||
it->setTextAlignment(Qt::AlignRight | Qt::AlignVCenter); // numbers right-aligned
|
||
table_->setItem(rowIdx, c++, it);
|
||
};
|
||
const double compl_ = s.possible_unique_reflections > 0
|
||
? 100.0 * s.unique_reflections / s.possible_unique_reflections : NAN;
|
||
const double mult = s.unique_reflections > 0
|
||
? static_cast<double>(s.total_observations) / s.unique_reflections : NAN;
|
||
set(resLabel);
|
||
set(QString::number(s.total_observations));
|
||
set(QString::number(s.unique_reflections));
|
||
set(std::isfinite(compl_) ? QString::number(compl_, 'f', 1) : QStringLiteral("—"));
|
||
set(std::isfinite(mult) ? QString::number(mult, 'f', 1) : QStringLiteral("—"));
|
||
set(QString::number(s.mean_i_over_sigma, 'f', 1));
|
||
set(QString::number(s.r_meas * 100.0, 'f', 1));
|
||
set(QString::number(s.cc_half * 100.0, 'f', 1));
|
||
if (has_reference_)
|
||
set(QString::number(s.cc_ref * 100.0, 'f', 1));
|
||
};
|
||
|
||
int rowIdx = 0;
|
||
for (const auto &s: stats_.shells)
|
||
fillRow(rowIdx++, s, QString::number(s.d_min, 'f', 2));
|
||
fillRow(rowIdx, stats_.overall, QStringLiteral("Overall"));
|
||
|
||
// Equal-width columns that fill the table (no super-wide last column).
|
||
table_->horizontalHeader()->setSectionResizeMode(QHeaderView::Stretch);
|
||
}
|