leonarski_fandClaude Opus 5
f706b2d8d2
FitShape set out.shells before running the regressions and returned early if either weighted least-squares call failed, leaving NaN coefficients with flat and convex both false. AnalyzeAnisotropy tested exactly those two, so it read the failure as Linear and marked it stable: the report and _reflns.jfjoch_aniso_shape claimed LINEAR over a shell count that suggested it meant something, with nan printed beside it. ShapeFit now carries an explicit ok, set only once every regression has succeeded, and result.shape_shells reports the shells the signature was actually fitted over - none, where it was not - so the renderer and the rebinning caution skip it and the verdict stays Undetermined. The directional-limit line tested the first and third principal axes for a finite d_min but not the second. ConeLimits fills each independently, so a crystal whose outer two directions fall in cones too sparse to measure while the middle one is well covered suppressed the whole line, the d_min spread with it. All three are tested now. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016L1qig74oYQzfUJJZbbxFh
1476 lines
74 KiB
C++
1476 lines
74 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "AnisotropyAnalysis.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <bit>
|
|
#include <cmath>
|
|
#include <cstdint>
|
|
#include <iomanip>
|
|
#include <numeric>
|
|
#include <tuple>
|
|
#include <sstream>
|
|
|
|
#include <Eigen/Dense>
|
|
|
|
#include "HKLKey.h"
|
|
|
|
#include "gemmi/eig3.hpp"
|
|
#include "gemmi/scaling.hpp"
|
|
|
|
namespace {
|
|
|
|
using Tensor = gemmi::SMat33<double>;
|
|
|
|
// ---------------------------------------------------------------- tuned constants
|
|
// The binning is the one the shape test was validated on: 12 equal-count resolution shells x 60
|
|
// directions on a hemisphere, a cell entering the fit once it holds 12 independent reflections.
|
|
// Shell count, direction count and the cell-population floor were each varied by a factor of two
|
|
// either way without moving a verdict, so none of the three is sharp.
|
|
constexpr int N_SHELLS = 12;
|
|
constexpr int N_DIRECTIONS = 60;
|
|
constexpr int MIN_CELL_REFLECTIONS = 12;
|
|
constexpr int MIN_CELLS = 40; // fewer than this and no tensor is fitted at all
|
|
constexpr int MIN_SHAPE_SHELLS = 5; // fewer than this and the s^2 signature is not fitted
|
|
constexpr int MIN_CELLS_PER_SHAPE_SHELL = 8;
|
|
|
|
// The whole diagnostic is refused below this: at <I/sigma> ~ 1 neither the tensor nor its shape is
|
|
// interpretable, and one battery dataset at 0.65 produced numbers that mean nothing.
|
|
constexpr double MIN_MEAN_I_OVER_SIGMA = 1.0;
|
|
|
|
// Below ~90 deg of OBSERVED rotation a lab-fixed systematic reaches a second and a third tensor
|
|
// direction (canonical correlation 0.74-0.85 against 0.05 at a full sweep), the anisotropy such a
|
|
// systematic manufactures grows 5-9x, and the power to detect a real 10 A^2 falls to <= 0.25.
|
|
constexpr double MIN_ROTATION_DEG = 90.0;
|
|
|
|
// Symmetry-expanded reflections are capped by striding the unique list: the estimator is a fit of a
|
|
// few tens of cell means and does not improve past this, while the memory does grow.
|
|
constexpr size_t MAX_ENTRIES = 4000000;
|
|
|
|
// The unmerged arm strides whole unique reflections down to this many. The systematic-error floor is
|
|
// set by the systematic, not by counting: over a 12x range in the number of unique reflections and a
|
|
// 10x range in <I/sigma> the smallest anisotropy the gate can establish moves by less than 2x, so
|
|
// more than this buys nothing and costs time and memory.
|
|
constexpr int MAX_CLUSTERS = 20000;
|
|
|
|
// The most partials ScaledObservations assembles on rotation data (see there). In partials because
|
|
// partials are what its sort costs, and generous because how many unique reflections that is
|
|
// depends on how many frames a rocking curve spans and how many mates the Laue class gives it,
|
|
// which together vary by an order of magnitude: over the rotation battery it leaves between 25
|
|
// thousand and 183 thousand of them, against the MAX_CLUSTERS the floor above actually reads.
|
|
constexpr size_t MAX_SAMPLE_PARTS = 8u << 20;
|
|
|
|
// splitmix64's finalizer. Sampling on this rather than on the key itself matters: a packed hkl key
|
|
// has h in its low bits, so any direct test on it selects a set of lattice planes, not a sample.
|
|
uint64_t MixKey(uint64_t x) {
|
|
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
|
|
x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
|
|
return x ^ (x >> 31);
|
|
}
|
|
|
|
// The cone half-angle of the directional diffraction limits, and the <I/sigma(I)> the limit is read
|
|
// at. Both are AIMLESS's shipped defaults (ANALYSIS CONE 20, ISIGMINIMUM 1.5 -> here 2.0, the value
|
|
// this project already uses for a resolution statement).
|
|
constexpr double CONE_HALF_ANGLE_DEG = 20.0;
|
|
constexpr double CONE_I_OVER_SIGMA = 2.0;
|
|
|
|
// The s^2 signature. FLAT is declared at c0 > +0.70 with c0/sigma(c0) > 6: over a 38-crystal battery
|
|
// that fires on six of the seven crystals whose deficit does not follow a Debye-Waller B, with no
|
|
// false positive among the other 24, and the false-positive count stays 0 anywhere between z > 5 and
|
|
// z > 8. A pure-B simulation on the same reflection lists returns |c0| <= 0.18, so the machinery
|
|
// cannot invent an intercept of this size. The curvature arm uses the same z for symmetry; it has no
|
|
// false-positive calibration of its own, and it is reported rather than acted on.
|
|
constexpr double FLAT_INTERCEPT = 0.70;
|
|
constexpr double SHAPE_Z = 6.0;
|
|
// CONVEX asks for the s^4 term to carry more than half the deficit at the resolution limit, not
|
|
// merely to be significant: a mild positive curvature is common and is not a different mechanism.
|
|
constexpr double CONVEX_SHARE = 0.5;
|
|
|
|
// Gate bands. The denominator is this dataset's own systematic error scale, not a counting-statistics
|
|
// error bar: against known ground truth a counting-statistics gate calls "real" in 15-58% of clean
|
|
// isotropic datasets, because the symmetry-forbidden directions are determined from contrasts BETWEEN
|
|
// symmetry mates of one reflection (which share their true |F|^2, so Wilson scatter cancels) while
|
|
// the symmetry-allowed ones can only be determined BETWEEN different reflections (which do not).
|
|
// 3.5 is a 2-5% false-positive rate pooled over six artefact classes and 6-14% in the worst of them;
|
|
// 5.0 is 0-0.5% and 0-3%.
|
|
constexpr double GATE_MARGINAL = 2.0;
|
|
constexpr double GATE_ESTABLISHED = 3.5;
|
|
constexpr double GATE_STRONG = 5.0;
|
|
|
|
// The forbidden subspace of a real dataset is dominated by systematics, not by counting noise, so
|
|
// this ratio is above 1 on essentially every dataset: measured over the battery its median is 6.5
|
|
// and only two of 22 exceed 20, so "it is above its counting null" would fire on almost all of
|
|
// them. This is the upper tail. Together with p == 1 and a gate that established nothing it picks
|
|
// out two datasets in the battery - one of them the crystal whose space-group determination is
|
|
// itself known to be unstable - where the unconditional p == 1 caution fired on 17 of 38.
|
|
constexpr double FORBIDDEN_Z_SUSPECT = 20.0;
|
|
|
|
// E[dB] / sigma_per_component for a zero-mean Gaussian tensor confined to a subspace: the eigenvalues
|
|
// of a noisily estimated tensor repel, so dB has a strictly positive expectation even when the tensor
|
|
// is exactly zero. Per-subspace Monte-Carlo medians, indexed by the dimension of the subspace;
|
|
// K_ALLOWED is for the symmetry-allowed subspaces, K_FORBIDDEN for the symmetry-forbidden ones (the
|
|
// two differ at the same dimension because they are different subspaces of the deviatoric space).
|
|
constexpr double K_ALLOWED[6] = {0.0, 0.824, 1.583, 2.075, 0.0, 2.872};
|
|
constexpr double K_FORBIDDEN[6] = {0.0, 0.0, 1.66, 2.13, 2.53, 2.87};
|
|
|
|
// ---------------------------------------------------------------- tensor algebra
|
|
double FrobDot(const Tensor &a, const Tensor &b) {
|
|
return a.u11 * b.u11 + a.u22 * b.u22 + a.u33 * b.u33
|
|
+ 2.0 * (a.u12 * b.u12 + a.u13 * b.u13 + a.u23 * b.u23);
|
|
}
|
|
|
|
Tensor Scale(const Tensor &a, double s) {
|
|
return {a.u11 * s, a.u22 * s, a.u33 * s, a.u12 * s, a.u13 * s, a.u23 * s};
|
|
}
|
|
|
|
Tensor Deviatoric(const Tensor &a) { return a.added_kI(-a.trace() / 3.0); }
|
|
|
|
// s^T E s from the cell-averaged outer product q = <(sx^2, sy^2, sz^2, sx sy, sx sz, sy sz)>.
|
|
double QuadForm(const Tensor &e, const std::array<double, 6> &q) {
|
|
return e.u11 * q[0] + e.u22 * q[1] + e.u33 * q[2]
|
|
+ 2.0 * (e.u12 * q[3] + e.u13 * q[4] + e.u23 * q[5]);
|
|
}
|
|
|
|
// Gram-Schmidt in the Frobenius inner product; directions that the earlier ones already span drop out.
|
|
std::vector<Tensor> Orthonormalise(const std::vector<Tensor> &raw) {
|
|
std::vector<Tensor> out;
|
|
for (Tensor t : raw) {
|
|
for (const auto &o : out)
|
|
t = t - Scale(o, FrobDot(t, o));
|
|
const double n = std::sqrt(FrobDot(t, t));
|
|
if (n > 1e-9)
|
|
out.push_back(Scale(t, 1.0 / n));
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// The cell with its lengths and angles snapped to what the crystal system requires.
|
|
gemmi::UnitCell MetricIdealCell(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) {
|
|
if (!sg)
|
|
return cell;
|
|
double a = cell.a, b = cell.b, c = cell.c;
|
|
double alpha = cell.alpha, beta = cell.beta, gamma = cell.gamma;
|
|
switch (sg->crystal_system()) {
|
|
case gemmi::CrystalSystem::Triclinic:
|
|
break;
|
|
case gemmi::CrystalSystem::Monoclinic:
|
|
if (sg->monoclinic_unique_axis() == 'a') { beta = 90.0; gamma = 90.0; }
|
|
else if (sg->monoclinic_unique_axis() == 'c') { alpha = 90.0; beta = 90.0; }
|
|
else { alpha = 90.0; gamma = 90.0; }
|
|
break;
|
|
case gemmi::CrystalSystem::Orthorhombic:
|
|
alpha = beta = gamma = 90.0;
|
|
break;
|
|
case gemmi::CrystalSystem::Tetragonal:
|
|
a = b = 0.5 * (a + b);
|
|
alpha = beta = gamma = 90.0;
|
|
break;
|
|
case gemmi::CrystalSystem::Trigonal:
|
|
case gemmi::CrystalSystem::Hexagonal:
|
|
if (sg->ext == 'R') {
|
|
a = b = c = (a + b + c) / 3.0;
|
|
alpha = beta = gamma = (alpha + beta + gamma) / 3.0;
|
|
} else {
|
|
a = b = 0.5 * (a + b);
|
|
alpha = beta = 90.0;
|
|
gamma = 120.0;
|
|
}
|
|
break;
|
|
case gemmi::CrystalSystem::Cubic:
|
|
a = b = c = (a + b + c) / 3.0;
|
|
alpha = beta = gamma = 90.0;
|
|
break;
|
|
}
|
|
gemmi::UnitCell out;
|
|
out.set(a, b, c, alpha, beta, gamma);
|
|
return out;
|
|
}
|
|
|
|
// The symmetry-allowed anisotropy directions, deviatoric, in the Cartesian frame.
|
|
//
|
|
// gemmi::adp_symmetry_constraints gives the standard ADP symmetry constraints in the crystal-axis
|
|
// (B*) parameterisation - 6 / 4 / 3 / 2 / 2 / 1 vectors for triclinic / monoclinic / orthorhombic /
|
|
// tetragonal / trigonal-hexagonal / cubic. Transforming each to Cartesian is a bijection on symmetric
|
|
// matrices, and the isotropic direction is always allowed, so removing the trace leaves exactly
|
|
// 5 / 3 / 2 / 1 / 1 / 0 free deviatoric parameters.
|
|
//
|
|
// Following Sheriff & Hendrickson (1987) Acta Cryst. A43, 118-121
|
|
std::vector<Tensor> AllowedDeviatoricBasis(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) {
|
|
std::vector<Tensor> raw;
|
|
// The constraints are a statement about the IDEAL metric, and rugnux writes the unconstrained
|
|
// refined cell: a beta of 89.97 deg in a C2 crystal leaves the isotropic tensor outside the span
|
|
// of the monoclinic constraints, the trace is then not fully removed, and the count comes out
|
|
// one too high. A cell whose angle is a hundredth of a degree off does not create a fourth
|
|
// anisotropy direction, so the basis is built on the idealised metric.
|
|
const gemmi::UnitCell ideal = MetricIdealCell(cell, sg);
|
|
for (const gemmi::Vec6 &v : gemmi::adp_symmetry_constraints(sg)) {
|
|
const Tensor b_star{v[0], v[1], v[2], v[3], v[4], v[5]};
|
|
raw.push_back(Deviatoric(b_star.transformed_by<double>(ideal.orth.mat)));
|
|
}
|
|
return Orthonormalise(raw);
|
|
}
|
|
|
|
std::vector<Tensor> FullDeviatoricBasis() {
|
|
const std::vector<Tensor> raw{{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0},
|
|
{0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}};
|
|
std::vector<Tensor> dev;
|
|
dev.reserve(raw.size());
|
|
for (const auto &t : raw)
|
|
dev.push_back(Deviatoric(t));
|
|
return Orthonormalise(dev);
|
|
}
|
|
|
|
// What is left of the deviatoric space once the symmetry-allowed directions are taken out. In those
|
|
// directions the true tensor is exactly zero by symmetry, whatever the crystal is - so whatever is
|
|
// measured there is this dataset's own systematic error, measured on this dataset.
|
|
std::vector<Tensor> ForbiddenBasis(const std::vector<Tensor> &full, const std::vector<Tensor> &allowed) {
|
|
std::vector<Tensor> raw;
|
|
for (Tensor t : full) {
|
|
for (const auto &a : allowed)
|
|
t = t - Scale(a, FrobDot(t, a));
|
|
raw.push_back(t);
|
|
}
|
|
return Orthonormalise(raw);
|
|
}
|
|
|
|
// ---------------------------------------------------------------- binning
|
|
// A Fibonacci spiral on the hemisphere: near-uniform, and no pole or seam where a lattice direction
|
|
// could pile up.
|
|
std::vector<gemmi::Vec3> HemisphereDirections(int n) {
|
|
std::vector<gemmi::Vec3> d;
|
|
d.reserve(n);
|
|
const double golden = gemmi::pi() * (1.0 + std::sqrt(5.0));
|
|
for (int i = 0; i < n; ++i) {
|
|
const double z = (i + 0.5) / n;
|
|
const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
|
|
const double phi = golden * (i + 0.5);
|
|
d.emplace_back(r * std::cos(phi), r * std::sin(phi), z);
|
|
}
|
|
return d;
|
|
}
|
|
|
|
struct Cells {
|
|
std::vector<double> mu, se, s2;
|
|
std::vector<std::array<double, 6>> q;
|
|
std::vector<int> shell, n;
|
|
// The cluster structure, kept so the covariance can be made robust to it: one entry per
|
|
// measurement that went into a cell, grouped by the unique reflection it belongs to.
|
|
std::vector<int32_t> entry_cell, entry_cluster;
|
|
std::vector<float> entry_value;
|
|
int n_shells = 0;
|
|
};
|
|
|
|
struct Usable {
|
|
std::array<int, 3> hkl;
|
|
double value; // epsilon-corrected merged intensity, negatives kept
|
|
double s2;
|
|
int shell = 0;
|
|
};
|
|
|
|
std::vector<Usable> SelectReflections(const std::vector<MergedReflection> &merged,
|
|
const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg,
|
|
int n_shells) {
|
|
const gemmi::GroupOps gops = sg ? sg->operations() : gemmi::GroupOps{};
|
|
std::vector<Usable> u;
|
|
u.reserve(merged.size());
|
|
for (const auto &r : merged) {
|
|
if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f))
|
|
continue;
|
|
const std::array<int, 3> hkl{r.h, r.k, r.l};
|
|
const double eps = sg ? std::max(1, gops.epsilon_factor(hkl)) : 1;
|
|
u.push_back({hkl, r.I / eps, 1.0 / (static_cast<double>(r.d) * r.d), 0});
|
|
}
|
|
if (u.empty())
|
|
return u;
|
|
// Equal-count shells in s^2. Every symmetry copy of a reflection has the same |s|, so equal-count
|
|
// over the expanded set and over the unique set are the same edges.
|
|
std::vector<int> order(u.size());
|
|
std::iota(order.begin(), order.end(), 0);
|
|
std::sort(order.begin(), order.end(), [&](int a, int b) { return u[a].s2 < u[b].s2; });
|
|
for (size_t i = 0; i < order.size(); ++i)
|
|
u[order[i]].shell = std::min(n_shells - 1, static_cast<int>(i * n_shells / order.size()));
|
|
return u;
|
|
}
|
|
|
|
// Keep the cells that hold enough measurements, renumber them, and carry the entry list over.
|
|
void FinishCells(Cells &c, const std::vector<double> &sum, const std::vector<double> &sumsq,
|
|
const std::vector<double> &s2sum, const std::vector<std::array<double, 6>> &qsum,
|
|
const std::vector<int> &count, const std::vector<int32_t> &raw_cell,
|
|
const std::vector<int32_t> &raw_cluster, const std::vector<float> &raw_value) {
|
|
std::vector<int32_t> remap(count.size(), -1);
|
|
for (size_t cid = 0; cid < count.size(); ++cid) {
|
|
if (count[cid] < MIN_CELL_REFLECTIONS)
|
|
continue;
|
|
const double n = count[cid];
|
|
const double mean = sum[cid] / n;
|
|
const double var = std::max(0.0, (sumsq[cid] - n * mean * mean) / (n - 1.0));
|
|
remap[cid] = static_cast<int32_t>(c.mu.size());
|
|
c.mu.push_back(mean);
|
|
c.se.push_back(std::sqrt(var / n));
|
|
c.s2.push_back(s2sum[cid] / n);
|
|
c.q.push_back({qsum[cid][0] / n, qsum[cid][1] / n, qsum[cid][2] / n,
|
|
qsum[cid][3] / n, qsum[cid][4] / n, qsum[cid][5] / n});
|
|
c.shell.push_back(static_cast<int>(cid) / N_DIRECTIONS);
|
|
c.n.push_back(count[cid]);
|
|
}
|
|
// A cell whose measurements happen to agree exactly has no error bar of its own; give it the
|
|
// typical one rather than an infinite weight.
|
|
std::vector<double> positive;
|
|
for (double s : c.se)
|
|
if (s > 0.0)
|
|
positive.push_back(s);
|
|
if (!positive.empty()) {
|
|
std::nth_element(positive.begin(), positive.begin() + positive.size() / 2, positive.end());
|
|
const double median = positive[positive.size() / 2];
|
|
for (double &s : c.se)
|
|
if (!(s > 0.0))
|
|
s = median;
|
|
}
|
|
for (size_t e = 0; e < raw_cell.size(); ++e)
|
|
if (remap[raw_cell[e]] >= 0) {
|
|
c.entry_cell.push_back(remap[raw_cell[e]]);
|
|
c.entry_cluster.push_back(raw_cluster[e]);
|
|
c.entry_value.push_back(raw_value[e]);
|
|
}
|
|
}
|
|
|
|
Cells BuildCells(const std::vector<Usable> &refl, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup *sg, int n_shells) {
|
|
Cells c;
|
|
c.n_shells = n_shells;
|
|
const std::vector<gemmi::Op> ops = sg ? sg->operations().sym_ops : std::vector<gemmi::Op>{gemmi::Op::identity()};
|
|
const auto dirs = HemisphereDirections(N_DIRECTIONS);
|
|
// Friedel needs no separate operation: -s folds onto the same hemisphere direction as s.
|
|
const size_t stride = std::max<size_t>(1, (refl.size() * ops.size() + MAX_ENTRIES - 1) / MAX_ENTRIES);
|
|
|
|
const int n_cells_total = n_shells * N_DIRECTIONS;
|
|
std::vector<double> sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0);
|
|
std::vector<std::array<double, 6>> qsum(n_cells_total, {0, 0, 0, 0, 0, 0});
|
|
std::vector<int> count(n_cells_total, 0);
|
|
std::vector<int32_t> raw_cell, raw_cluster;
|
|
std::vector<float> raw_value;
|
|
|
|
std::vector<std::pair<int, gemmi::Vec3>> hits; // (cell, s) of this reflection's copies
|
|
for (size_t i = 0; i < refl.size(); i += stride) {
|
|
const auto &r = refl[i];
|
|
hits.clear();
|
|
for (const auto &op : ops) {
|
|
const gemmi::Op::Miller m = op.apply_to_hkl(r.hkl);
|
|
gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2]));
|
|
if (s.z < 0)
|
|
s = gemmi::Vec3(-s.x, -s.y, -s.z);
|
|
const double len = s.length();
|
|
if (!(len > 0))
|
|
continue;
|
|
const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len);
|
|
int best = 0;
|
|
double best_dot = -1.0;
|
|
for (int j = 0; j < N_DIRECTIONS; ++j) {
|
|
const double dot = std::fabs(unit.dot(dirs[j]));
|
|
if (dot > best_dot) { best_dot = dot; best = j; }
|
|
}
|
|
hits.emplace_back(r.shell * N_DIRECTIONS + best, s);
|
|
}
|
|
// A reflection contributes to a cell ONCE however many of its symmetry copies land there, so
|
|
// a cell's standard error counts independent measurements.
|
|
std::sort(hits.begin(), hits.end(), [](const auto &a, const auto &b) { return a.first < b.first; });
|
|
const int32_t cluster = static_cast<int32_t>(i);
|
|
int previous = -1;
|
|
for (const auto &[cid, s] : hits) {
|
|
if (cid == previous)
|
|
continue;
|
|
previous = cid;
|
|
sum[cid] += r.value;
|
|
sumsq[cid] += r.value * r.value;
|
|
s2sum[cid] += r.s2;
|
|
qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z;
|
|
qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z;
|
|
count[cid] += 1;
|
|
raw_cell.push_back(static_cast<int32_t>(cid));
|
|
raw_cluster.push_back(cluster);
|
|
raw_value.push_back(static_cast<float>(r.value));
|
|
}
|
|
}
|
|
|
|
FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value);
|
|
return c;
|
|
}
|
|
|
|
// ---------------------------------------------------------------- the tensor fit
|
|
struct TensorFit {
|
|
bool ok = false;
|
|
Tensor b{0, 0, 0, 0, 0, 0};
|
|
Eigen::MatrixXd cov; // nb x nb, cluster-robust
|
|
double chi2red = NAN;
|
|
};
|
|
|
|
// mu_cell = K_shell * exp(-1/2 s^T B s), fitted on the INTENSITY scale by weighted Gauss-Newton with
|
|
// weights 1/SE^2, one free K per resolution shell profiled alongside the tensor. Because K is free per
|
|
// shell, every isotropic feature - the Wilson curve, an ice ring, a noise floor, a scaling error - is
|
|
// absorbed exactly, and only the l=2 angular part drives the tensor. Nothing is censored: a cell whose
|
|
// mean is zero or negative is kept, which is what keeps the fit honest in a direction that has died.
|
|
//
|
|
// The covariance is a cluster-robust sandwich, clustering by unique reflection. That is not a detail:
|
|
// symmetry mates of one reflection share their true |F|^2, so a symmetry-forbidden tensor direction -
|
|
// determined from contrasts WITHIN a cluster - carries none of the Wilson scatter that a
|
|
// symmetry-allowed direction carries, and a covariance that assumed independent cells would get the
|
|
// ratio of the two backwards by a factor of several.
|
|
//
|
|
// Following Popov & Bourenkov (2003) Acta Cryst. D59, 1145-1153
|
|
TensorFit FitTensor(const Cells &c, const std::vector<Tensor> &basis, bool want_covariance) {
|
|
TensorFit fit;
|
|
const int nb = static_cast<int>(basis.size());
|
|
const int nc = static_cast<int>(c.mu.size());
|
|
if (nb == 0 || nc < MIN_CELLS)
|
|
return fit;
|
|
|
|
// Shells that carry a positive mean intensity: ln K is undefined for one that does not, and a
|
|
// shell with no signal at all cannot say anything about direction either.
|
|
std::vector<double> shell_num(c.n_shells, 0.0), shell_den(c.n_shells, 0.0);
|
|
for (int i = 0; i < nc; ++i) {
|
|
const double w = 1.0 / (c.se[i] * c.se[i]);
|
|
shell_num[c.shell[i]] += w * c.mu[i];
|
|
shell_den[c.shell[i]] += w;
|
|
}
|
|
std::vector<int> slot(c.n_shells, -1);
|
|
int ns = 0;
|
|
std::vector<double> ln_k;
|
|
for (int s = 0; s < c.n_shells; ++s)
|
|
if (shell_den[s] > 0.0 && shell_num[s] / shell_den[s] > 0.0) {
|
|
slot[s] = ns++;
|
|
ln_k.push_back(std::log(shell_num[s] / shell_den[s]));
|
|
}
|
|
std::vector<int> keep;
|
|
for (int i = 0; i < nc; ++i)
|
|
if (slot[c.shell[i]] >= 0)
|
|
keep.push_back(i);
|
|
if (static_cast<int>(keep.size()) < MIN_CELLS || ns == 0)
|
|
return fit;
|
|
|
|
const int n = static_cast<int>(keep.size());
|
|
const int np = ns + nb;
|
|
Eigen::MatrixXd G(n, nb);
|
|
Eigen::VectorXd mu(n), se(n);
|
|
std::vector<int> row_slot(n);
|
|
for (int i = 0; i < n; ++i) {
|
|
const int ci = keep[i];
|
|
for (int k = 0; k < nb; ++k)
|
|
G(i, k) = -0.5 * QuadForm(basis[k], c.q[ci]);
|
|
mu(i) = c.mu[ci];
|
|
se(i) = c.se[ci];
|
|
row_slot[i] = slot[c.shell[ci]];
|
|
}
|
|
|
|
Eigen::VectorXd beta = Eigen::VectorXd::Zero(np);
|
|
for (int s = 0; s < ns; ++s)
|
|
beta(s) = ln_k[s];
|
|
|
|
auto predict = [&](const Eigen::VectorXd &b, Eigen::VectorXd &pred) {
|
|
for (int i = 0; i < n; ++i) {
|
|
double e = b(row_slot[i]);
|
|
for (int k = 0; k < nb; ++k)
|
|
e += G(i, k) * b(ns + k);
|
|
pred(i) = std::exp(std::clamp(e, -60.0, 60.0));
|
|
}
|
|
};
|
|
auto chi2_of = [&](const Eigen::VectorXd &pred) {
|
|
double s = 0.0;
|
|
for (int i = 0; i < n; ++i) {
|
|
const double r = (mu(i) - pred(i)) / se(i);
|
|
s += r * r;
|
|
}
|
|
return s;
|
|
};
|
|
|
|
Eigen::VectorXd pred(n), resid(n);
|
|
predict(beta, pred);
|
|
double chi2 = chi2_of(pred);
|
|
Eigen::MatrixXd J(n, np);
|
|
for (int iter = 0; iter < 60; ++iter) {
|
|
J.setZero();
|
|
for (int i = 0; i < n; ++i) {
|
|
const double f = -pred(i) / se(i);
|
|
J(i, row_slot[i]) = f;
|
|
for (int k = 0; k < nb; ++k)
|
|
J(i, ns + k) = f * G(i, k);
|
|
resid(i) = (mu(i) - pred(i)) / se(i);
|
|
}
|
|
const Eigen::MatrixXd A = J.transpose() * J;
|
|
const Eigen::VectorXd grad = J.transpose() * resid;
|
|
const Eigen::VectorXd step =
|
|
A.ldlt().solve(-grad); // Gauss-Newton on chi2/2
|
|
if (!step.allFinite())
|
|
break;
|
|
// Halve the step until chi2 falls; a fit that cannot improve at all is converged.
|
|
double t = 1.0;
|
|
bool improved = false;
|
|
Eigen::VectorXd trial(np), trial_pred(n);
|
|
for (int back = 0; back < 20; ++back) {
|
|
trial = beta + t * step;
|
|
predict(trial, trial_pred);
|
|
const double trial_chi2 = chi2_of(trial_pred);
|
|
if (trial_chi2 < chi2) {
|
|
const bool converged = chi2 - trial_chi2 < 1e-10 * std::max(1.0, chi2);
|
|
beta = trial;
|
|
pred = trial_pred;
|
|
chi2 = trial_chi2;
|
|
improved = !converged;
|
|
break;
|
|
}
|
|
t *= 0.5;
|
|
}
|
|
if (!improved)
|
|
break;
|
|
}
|
|
|
|
Tensor b{0, 0, 0, 0, 0, 0};
|
|
for (int k = 0; k < nb; ++k)
|
|
b = b + Scale(basis[k], beta(ns + k));
|
|
fit.b = b;
|
|
fit.chi2red = chi2 / std::max(1, n - np);
|
|
fit.ok = true;
|
|
if (!want_covariance)
|
|
return fit;
|
|
|
|
// Sandwich: A^-1 B A^-1 with A the Gauss-Newton normal matrix and B the sum of outer products of
|
|
// per-cluster scores. A cell mean is the mean of its reflections, so a reflection's share of the
|
|
// cell's score is its own share of the cell's residual.
|
|
for (int i = 0; i < n; ++i) {
|
|
const double f = -pred(i) / se(i);
|
|
J(i, row_slot[i]) = f;
|
|
for (int k = 0; k < nb; ++k)
|
|
J(i, ns + k) = f * G(i, k);
|
|
}
|
|
const Eigen::MatrixXd A = J.transpose() * J;
|
|
std::vector<int32_t> cell_row(c.mu.size(), -1);
|
|
for (int i = 0; i < n; ++i)
|
|
cell_row[keep[i]] = i;
|
|
Eigen::MatrixXd B = Eigen::MatrixXd::Zero(np, np);
|
|
Eigen::VectorXd score = Eigen::VectorXd::Zero(np);
|
|
int32_t current = -1;
|
|
auto flush = [&]() {
|
|
if (current >= 0)
|
|
B.noalias() += score * score.transpose();
|
|
score.setZero();
|
|
};
|
|
for (size_t e = 0; e < c.entry_cell.size(); ++e) {
|
|
if (c.entry_cluster[e] != current) {
|
|
flush();
|
|
current = c.entry_cluster[e];
|
|
}
|
|
const int i = cell_row[c.entry_cell[e]];
|
|
if (i < 0)
|
|
continue;
|
|
const int ci = keep[i];
|
|
const double w = (c.entry_value[e] - pred(i)) * pred(i)
|
|
/ (c.n[ci] * c.se[ci] * c.se[ci]);
|
|
score(row_slot[i]) += w;
|
|
for (int k = 0; k < nb; ++k)
|
|
score(ns + k) += w * G(i, k);
|
|
}
|
|
flush();
|
|
const Eigen::MatrixXd ainv = A.ldlt().solve(Eigen::MatrixXd::Identity(np, np));
|
|
const Eigen::MatrixXd full = ainv * B * ainv;
|
|
fit.cov = full.bottomRightCorner(nb, nb);
|
|
return fit;
|
|
}
|
|
|
|
struct Eigen3 {
|
|
double value[3];
|
|
double vec[3][3];
|
|
};
|
|
|
|
Eigen3 Decompose(const Tensor &b) {
|
|
Eigen3 out{};
|
|
double d[3];
|
|
const gemmi::Mat33 v = gemmi::eigen_decomposition(b, d);
|
|
int idx[3] = {0, 1, 2};
|
|
std::sort(idx, idx + 3, [&](int a, int c) { return d[a] > d[c]; });
|
|
for (int n = 0; n < 3; ++n) {
|
|
out.value[n] = d[idx[n]];
|
|
for (int j = 0; j < 3; ++j)
|
|
out.vec[n][j] = v[j][idx[n]]; // eigenvectors are the columns
|
|
}
|
|
return out;
|
|
}
|
|
|
|
double DeltaB(const Tensor &b) {
|
|
double d[3];
|
|
gemmi::eigen_decomposition(b, d);
|
|
return *std::max_element(d, d + 3) - *std::min_element(d, d + 3);
|
|
}
|
|
|
|
// The dB that counting statistics alone would produce in a subspace: E[dB] = k * sigma, with sigma the
|
|
// per-component standard error and k the constant for that subspace's dimension.
|
|
double CountingNull(const Eigen::MatrixXd &cov, double k) {
|
|
if (cov.rows() == 0)
|
|
return 0.0;
|
|
return k * std::sqrt(std::max(0.0, cov.trace() / cov.rows()));
|
|
}
|
|
|
|
// ---------------------------------------------------------------- the s^2 signature
|
|
struct ShapeFit {
|
|
// The regressions below all ran. Distinct from shells: there can be plenty of shells and
|
|
// still no fit, and a ShapeFit that returned early carries NaN coefficients with flat and
|
|
// convex both false - which reads as a clean "Linear" to anyone testing those two alone.
|
|
bool ok = false;
|
|
int shells = 0;
|
|
double c0 = NAN, c0_err = NAN, c1 = NAN, slope_through_origin = NAN;
|
|
double curvature = NAN, curvature_err = NAN, curvature_share = NAN;
|
|
double residual = NAN;
|
|
bool flat = false, convex = false;
|
|
};
|
|
|
|
// One shell's l=2 amplitude: mu_cell = K exp(A u_cell), u = -1/2 (s^T bhat s)/|s|^2, K profiled out,
|
|
// fitted on the intensity scale so no cell is dropped. A genuine tensor gives A = |B_dev| s^2, i.e. a
|
|
// straight line through the origin; the fit makes no assumption at all about how A depends on s.
|
|
bool FitShellAmplitude(const std::vector<double> &u, const std::vector<double> &mu,
|
|
const std::vector<double> &se, double &a_out, double &a_err) {
|
|
auto chi2 = [&](double a) {
|
|
double num = 0.0, den = 0.0;
|
|
std::vector<double> g(u.size());
|
|
for (size_t i = 0; i < u.size(); ++i) {
|
|
g[i] = std::exp(std::clamp(a * u[i], -60.0, 60.0));
|
|
const double w = 1.0 / (se[i] * se[i]);
|
|
num += mu[i] * g[i] * w;
|
|
den += g[i] * g[i] * w;
|
|
}
|
|
const double k = den > 0.0 ? num / den : 0.0;
|
|
double s = 0.0;
|
|
for (size_t i = 0; i < u.size(); ++i) {
|
|
const double r = (mu[i] - k * g[i]) / se[i];
|
|
s += r * r;
|
|
}
|
|
return s;
|
|
};
|
|
// A coarse scan brackets the minimum, golden section refines it. A is |B_dev| s^2, so a few
|
|
// hundred covers any crystal; the search is on a smooth one-parameter curve.
|
|
double lo = -300.0, hi = 300.0, best = 0.0, best_chi2 = chi2(0.0);
|
|
for (int i = 0; i <= 120; ++i) {
|
|
const double a = lo + (hi - lo) * i / 120.0;
|
|
const double v = chi2(a);
|
|
if (v < best_chi2) { best_chi2 = v; best = a; }
|
|
}
|
|
lo = best - 5.0;
|
|
hi = best + 5.0;
|
|
constexpr double INV_PHI = 0.6180339887498949;
|
|
double x1 = hi - INV_PHI * (hi - lo), x2 = lo + INV_PHI * (hi - lo);
|
|
double f1 = chi2(x1), f2 = chi2(x2);
|
|
for (int i = 0; i < 60; ++i) {
|
|
if (f1 < f2) { hi = x2; x2 = x1; f2 = f1; x1 = hi - INV_PHI * (hi - lo); f1 = chi2(x1); }
|
|
else { lo = x1; x1 = x2; f1 = f2; x2 = lo + INV_PHI * (hi - lo); f2 = chi2(x2); }
|
|
}
|
|
a_out = 0.5 * (lo + hi);
|
|
// One-parameter Gauss-Newton error bar, scaled by the residual so a shell the model does not
|
|
// describe reports a large one.
|
|
double num = 0.0, den = 0.0;
|
|
std::vector<double> g(u.size());
|
|
for (size_t i = 0; i < u.size(); ++i) {
|
|
g[i] = std::exp(std::clamp(a_out * u[i], -60.0, 60.0));
|
|
const double w = 1.0 / (se[i] * se[i]);
|
|
num += mu[i] * g[i] * w;
|
|
den += g[i] * g[i] * w;
|
|
}
|
|
const double k = den > 0.0 ? num / den : 0.0;
|
|
double jtj = 0.0, ss = 0.0;
|
|
for (size_t i = 0; i < u.size(); ++i) {
|
|
const double dr = k * g[i] * u[i] / se[i];
|
|
jtj += dr * dr;
|
|
const double r = (mu[i] - k * g[i]) / se[i];
|
|
ss += r * r;
|
|
}
|
|
if (!(jtj > 0.0) || u.size() < 3)
|
|
return false;
|
|
a_err = std::sqrt(ss / (u.size() - 2) / jtj);
|
|
return std::isfinite(a_out) && std::isfinite(a_err) && a_err > 0.0;
|
|
}
|
|
|
|
// Weighted least squares of `y` on the columns of `x`; returns the parameters and the covariance.
|
|
bool Wls(const Eigen::MatrixXd &x, const Eigen::VectorXd &y, const Eigen::VectorXd &w,
|
|
Eigen::VectorXd &p, Eigen::MatrixXd &cov, double &chi2) {
|
|
const Eigen::MatrixXd xtw = x.transpose() * w.asDiagonal();
|
|
const Eigen::MatrixXd normal = xtw * x;
|
|
const Eigen::LDLT<Eigen::MatrixXd> ldlt(normal);
|
|
if (ldlt.info() != Eigen::Success)
|
|
return false;
|
|
p = ldlt.solve(xtw * y);
|
|
const Eigen::VectorXd r = y - x * p;
|
|
chi2 = r.transpose() * w.asDiagonal() * r;
|
|
cov = ldlt.solve(Eigen::MatrixXd::Identity(x.cols(), x.cols()));
|
|
return p.allFinite();
|
|
}
|
|
|
|
ShapeFit FitShape(const Cells &c, const Tensor &tensor) {
|
|
ShapeFit out;
|
|
const double norm = std::sqrt(FrobDot(tensor, tensor));
|
|
if (!(norm > 0.0))
|
|
return out;
|
|
const Tensor bhat = Scale(tensor, 1.0 / norm);
|
|
|
|
std::vector<double> xs, as, errs;
|
|
for (int s = 0; s < c.n_shells; ++s) {
|
|
std::vector<double> u, mu, se;
|
|
double s2sum = 0.0;
|
|
for (size_t i = 0; i < c.mu.size(); ++i) {
|
|
if (c.shell[i] != s)
|
|
continue;
|
|
u.push_back(-0.5 * QuadForm(bhat, c.q[i]) / std::max(c.s2[i], 1e-12));
|
|
mu.push_back(c.mu[i]);
|
|
se.push_back(c.se[i]);
|
|
s2sum += c.s2[i];
|
|
}
|
|
if (static_cast<int>(u.size()) < MIN_CELLS_PER_SHAPE_SHELL)
|
|
continue;
|
|
const auto [lo, hi] = std::minmax_element(u.begin(), u.end());
|
|
if (*hi - *lo < 1e-6)
|
|
continue;
|
|
double a = 0.0, err = 0.0;
|
|
if (!FitShellAmplitude(u, mu, se, a, err))
|
|
continue;
|
|
xs.push_back(s2sum / u.size());
|
|
as.push_back(a);
|
|
errs.push_back(err);
|
|
}
|
|
out.shells = static_cast<int>(xs.size());
|
|
if (out.shells < MIN_SHAPE_SHELLS)
|
|
return out;
|
|
|
|
const int n = out.shells;
|
|
Eigen::VectorXd y(n), w(n);
|
|
for (int i = 0; i < n; ++i) {
|
|
y(i) = as[i];
|
|
w(i) = 1.0 / (errs[i] * errs[i]);
|
|
}
|
|
Eigen::VectorXd p;
|
|
Eigen::MatrixXd cov;
|
|
double chi2 = 0.0;
|
|
Eigen::MatrixXd x1(n, 1), x2(n, 2), xq(n, 2);
|
|
for (int i = 0; i < n; ++i) {
|
|
x1(i, 0) = xs[i];
|
|
x2(i, 0) = 1.0;
|
|
x2(i, 1) = xs[i];
|
|
xq(i, 0) = xs[i];
|
|
xq(i, 1) = xs[i] * xs[i];
|
|
}
|
|
if (!Wls(x1, y, w, p, cov, chi2))
|
|
return out;
|
|
out.slope_through_origin = p(0);
|
|
out.residual = chi2 / std::max(1, n - 1);
|
|
if (!Wls(x2, y, w, p, cov, chi2))
|
|
return out;
|
|
out.c0 = p(0);
|
|
out.c0_err = std::sqrt(std::max(0.0, cov(0, 0)));
|
|
out.c1 = p(1);
|
|
if (Wls(xq, y, w, p, cov, chi2)) {
|
|
out.curvature = p(1);
|
|
out.curvature_err = std::sqrt(std::max(0.0, cov(1, 1)));
|
|
// How much of the deficit at the resolution limit the s^4 term carries.
|
|
const double s2max = *std::max_element(xs.begin(), xs.end());
|
|
const double total = p(0) * s2max + p(1) * s2max * s2max;
|
|
out.curvature_share = total != 0.0 ? p(1) * s2max * s2max / total : NAN;
|
|
}
|
|
out.flat = out.c0 > FLAT_INTERCEPT && out.c0_err > 0.0 && out.c0 / out.c0_err > SHAPE_Z;
|
|
out.convex = !out.flat && out.curvature > 0.0 && out.curvature_err > 0.0
|
|
&& out.curvature / out.curvature_err > SHAPE_Z
|
|
&& out.curvature_share > CONVEX_SHARE;
|
|
out.ok = true;
|
|
return out;
|
|
}
|
|
|
|
// ---------------------------------------------------------------- directional limits
|
|
// <I/sigma(I)> in a cone about each principal direction, read where it falls through 2.0. This uses no
|
|
// model of the fall-off at all, which is why it and the tensor are reported side by side.
|
|
//
|
|
// Following Evans & Murshudov (2013) Acta Cryst. D69, 1204-1214
|
|
void ConeLimits(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup *sg, const Eigen3 &axes, int n_shells,
|
|
double (&d_min)[3], bool (&censored)[3]) {
|
|
const std::vector<gemmi::Op> ops = sg ? sg->operations().sym_ops
|
|
: std::vector<gemmi::Op>{gemmi::Op::identity()};
|
|
std::vector<double> s2;
|
|
s2.reserve(merged.size());
|
|
for (const auto &r : merged)
|
|
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f)
|
|
s2.push_back(1.0 / (static_cast<double>(r.d) * r.d));
|
|
if (s2.size() < static_cast<size_t>(n_shells * MIN_CELL_REFLECTIONS))
|
|
return;
|
|
std::vector<double> sorted = s2;
|
|
std::sort(sorted.begin(), sorted.end());
|
|
std::vector<double> edge(n_shells + 1);
|
|
for (int i = 0; i <= n_shells; ++i)
|
|
edge[i] = sorted[std::min(sorted.size() - 1, sorted.size() * i / n_shells)];
|
|
|
|
const double cos_cone = std::cos(CONE_HALF_ANGLE_DEG * gemmi::pi() / 180.0);
|
|
std::vector<std::vector<double>> isig_sum(3, std::vector<double>(n_shells, 0.0));
|
|
std::vector<std::vector<int>> isig_n(3, std::vector<int>(n_shells, 0));
|
|
size_t iu = 0;
|
|
for (const auto &r : merged) {
|
|
if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f))
|
|
continue;
|
|
const double this_s2 = s2[iu++];
|
|
int shell = 0;
|
|
while (shell + 1 < n_shells && this_s2 > edge[shell + 1])
|
|
++shell;
|
|
bool in_cone[3] = {false, false, false};
|
|
for (const auto &op : ops) {
|
|
const gemmi::Op::Miller m = op.apply_to_hkl({{r.h, r.k, r.l}});
|
|
const gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2]));
|
|
const double len = s.length();
|
|
if (!(len > 0))
|
|
continue;
|
|
for (int n = 0; n < 3; ++n) {
|
|
const double dot = (s.x * axes.vec[n][0] + s.y * axes.vec[n][1] + s.z * axes.vec[n][2]) / len;
|
|
if (std::fabs(dot) >= cos_cone)
|
|
in_cone[n] = true;
|
|
}
|
|
}
|
|
for (int n = 0; n < 3; ++n)
|
|
if (in_cone[n]) {
|
|
isig_sum[n][shell] += r.I / r.sigma;
|
|
isig_n[n][shell] += 1;
|
|
}
|
|
}
|
|
for (int n = 0; n < 3; ++n) {
|
|
double last_s2 = NAN, last_isig = NAN, limit_s2 = NAN;
|
|
bool crossed = false;
|
|
for (int s = 0; s < n_shells; ++s) {
|
|
if (isig_n[n][s] < MIN_CELL_REFLECTIONS)
|
|
continue;
|
|
const double mid = 0.5 * (edge[s] + edge[s + 1]);
|
|
const double isig = isig_sum[n][s] / isig_n[n][s];
|
|
if (isig >= CONE_I_OVER_SIGMA) {
|
|
last_s2 = mid;
|
|
last_isig = isig;
|
|
limit_s2 = edge[s + 1];
|
|
} else if (std::isfinite(last_s2)) {
|
|
// Linear interpolation of <I/sigma> in s^2 between the last shell above the threshold
|
|
// and the first below it.
|
|
const double f = (last_isig - CONE_I_OVER_SIGMA) / (last_isig - isig);
|
|
limit_s2 = last_s2 + f * (mid - last_s2);
|
|
crossed = true;
|
|
break;
|
|
} else {
|
|
break;
|
|
}
|
|
}
|
|
if (std::isfinite(limit_s2) && limit_s2 > 0.0) {
|
|
d_min[n] = 1.0 / std::sqrt(limit_s2);
|
|
// <I/sigma> never fell through the threshold: the outermost shell that carries this cone
|
|
// is still above it, so the value is the edge of the measured data and not this crystal's
|
|
// own limit along that direction.
|
|
censored[n] = !crossed;
|
|
}
|
|
}
|
|
}
|
|
|
|
const char *Band(double g) {
|
|
if (g < GATE_MARGINAL) return "not established";
|
|
if (g < GATE_ESTABLISHED) return "marginal";
|
|
if (g < GATE_STRONG) return "established";
|
|
return "strong";
|
|
}
|
|
|
|
// ---------------------------------------------------------------- the systematic-error floor
|
|
// The same cells, built from UNMERGED observations at the Miller index each was measured at, with
|
|
// the unique reflection they reduce to as the cluster. Nothing is symmetry-expanded and nothing is
|
|
// de-duplicated: two observations of the same reflection recorded on different frames are two
|
|
// measurements, and the difference between them is exactly the signal this arm is after.
|
|
Cells BuildObservationCells(const std::vector<AnisotropyObservation> &obs, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup *sg, int n_shells) {
|
|
Cells c;
|
|
c.n_shells = n_shells;
|
|
const HKLKeyGenerator key(/*merge_friedel=*/true, sg ? *sg : *gemmi::find_spacegroup_by_number(1));
|
|
|
|
// Order the observations by the unique reflection they belong to, so the cluster-robust sum
|
|
// below can walk them one cluster at a time; then stride whole clusters if there are more than
|
|
// the fit needs. The estimate is flat in the number of unique reflections over a factor of 12,
|
|
// so striding costs nothing and bounds the work.
|
|
std::vector<std::pair<uint64_t, int32_t>> order;
|
|
order.reserve(obs.size());
|
|
for (size_t i = 0; i < obs.size(); ++i) {
|
|
const auto &o = obs[i];
|
|
if (!std::isfinite(o.I) || !std::isfinite(o.sigma) || o.sigma <= 0.0f || !(o.d > 0.0f))
|
|
continue;
|
|
order.emplace_back(key(o.h, o.k, o.l).pack(), static_cast<int32_t>(i));
|
|
}
|
|
if (order.empty())
|
|
return c;
|
|
std::sort(order.begin(), order.end());
|
|
int32_t n_clusters = 1;
|
|
for (size_t i = 1; i < order.size(); ++i)
|
|
if (order[i].first != order[i - 1].first)
|
|
++n_clusters;
|
|
const int cluster_stride = std::max<int>(1, n_clusters / MAX_CLUSTERS);
|
|
|
|
const auto dirs = HemisphereDirections(N_DIRECTIONS);
|
|
const int n_cells_total = n_shells * N_DIRECTIONS;
|
|
std::vector<double> sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0);
|
|
std::vector<std::array<double, 6>> qsum(n_cells_total, {0, 0, 0, 0, 0, 0});
|
|
std::vector<int> count(n_cells_total, 0);
|
|
std::vector<int32_t> raw_cell, raw_cluster;
|
|
std::vector<float> raw_value;
|
|
|
|
std::vector<int32_t> kept;
|
|
std::vector<int32_t> cluster_of;
|
|
uint64_t previous_key = 0;
|
|
int32_t seen = -1;
|
|
for (const auto &[k, idx] : order) {
|
|
if (seen < 0 || k != previous_key) {
|
|
previous_key = k;
|
|
++seen;
|
|
}
|
|
if (seen % cluster_stride != 0)
|
|
continue;
|
|
kept.push_back(idx);
|
|
cluster_of.push_back(seen / cluster_stride);
|
|
}
|
|
if (kept.size() < static_cast<size_t>(MIN_CELLS * MIN_CELL_REFLECTIONS))
|
|
return c;
|
|
|
|
// Equal-count shells over the observations that are kept.
|
|
std::vector<double> s2(kept.size());
|
|
for (size_t i = 0; i < kept.size(); ++i) {
|
|
const auto &o = obs[kept[i]];
|
|
s2[i] = 1.0 / (static_cast<double>(o.d) * o.d);
|
|
}
|
|
std::vector<int32_t> by_s2(kept.size());
|
|
std::iota(by_s2.begin(), by_s2.end(), 0);
|
|
std::sort(by_s2.begin(), by_s2.end(), [&](int32_t a, int32_t b) { return s2[a] < s2[b]; });
|
|
std::vector<int> shell(kept.size());
|
|
for (size_t i = 0; i < by_s2.size(); ++i)
|
|
shell[by_s2[i]] = std::min(n_shells - 1, static_cast<int>(i * n_shells / by_s2.size()));
|
|
|
|
raw_cell.reserve(kept.size());
|
|
raw_value.reserve(kept.size());
|
|
raw_cluster.reserve(kept.size());
|
|
for (size_t i = 0; i < kept.size(); ++i) {
|
|
const auto &o = obs[kept[i]];
|
|
gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(o.h, o.k, o.l));
|
|
if (s.z < 0)
|
|
s = gemmi::Vec3(-s.x, -s.y, -s.z);
|
|
const double len = s.length();
|
|
if (!(len > 0))
|
|
continue;
|
|
const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len);
|
|
int best = 0;
|
|
double best_dot = -1.0;
|
|
for (int j = 0; j < N_DIRECTIONS; ++j) {
|
|
const double dot = std::fabs(unit.dot(dirs[j]));
|
|
if (dot > best_dot) { best_dot = dot; best = j; }
|
|
}
|
|
const int cid = shell[i] * N_DIRECTIONS + best;
|
|
sum[cid] += o.I;
|
|
sumsq[cid] += static_cast<double>(o.I) * o.I;
|
|
s2sum[cid] += s2[i];
|
|
qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z;
|
|
qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z;
|
|
count[cid] += 1;
|
|
raw_cell.push_back(cid);
|
|
raw_cluster.push_back(cluster_of[i]);
|
|
raw_value.push_back(o.I);
|
|
}
|
|
FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value);
|
|
return c;
|
|
}
|
|
|
|
// The systematic error scale this dataset carries, and the part of it that is merely counting
|
|
// noise, both measured in the tensor directions the Laue class forbids.
|
|
struct SystematicFloor {
|
|
bool ok = false;
|
|
double sigma_excess = NAN; // per tensor component, A^2
|
|
double forbidden_z = NAN; // that scatter over its own counting null
|
|
int n_observations = 0;
|
|
};
|
|
|
|
SystematicFloor MeasureSystematicFloor(const std::vector<AnisotropyObservation> &obs,
|
|
const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg,
|
|
const std::vector<Tensor> &full,
|
|
const std::vector<Tensor> &forbidden) {
|
|
SystematicFloor out;
|
|
const int q = static_cast<int>(forbidden.size());
|
|
if (q == 0 || obs.empty())
|
|
return out;
|
|
const Cells cells = BuildObservationCells(obs, cell, sg, N_SHELLS);
|
|
if (cells.mu.empty())
|
|
return out;
|
|
const TensorFit fit = FitTensor(cells, full, true);
|
|
if (!fit.ok || fit.cov.rows() != static_cast<int>(full.size()))
|
|
return out;
|
|
|
|
Eigen::MatrixXd basis_map(static_cast<int>(full.size()), q);
|
|
for (int j = 0; j < q; ++j)
|
|
for (size_t k = 0; k < full.size(); ++k)
|
|
basis_map(static_cast<int>(k), j) = FrobDot(full[k], forbidden[j]);
|
|
Tensor forb{0, 0, 0, 0, 0, 0};
|
|
Eigen::VectorXd theta(static_cast<int>(full.size()));
|
|
for (size_t k = 0; k < full.size(); ++k)
|
|
theta(static_cast<int>(k)) = FrobDot(fit.b, full[k]);
|
|
const Eigen::VectorXd theta_forb = basis_map.transpose() * theta;
|
|
for (int j = 0; j < q; ++j)
|
|
forb = forb + Scale(forbidden[j], theta_forb(j));
|
|
const Eigen::MatrixXd cov_forb = basis_map.transpose() * fit.cov * basis_map;
|
|
|
|
const double k_q = K_FORBIDDEN[q];
|
|
if (!(k_q > 0.0))
|
|
return out;
|
|
const double sigma_sys = DeltaB(forb) / k_q;
|
|
const double sigma_stat = CountingNull(cov_forb, k_q) / k_q;
|
|
// Take the counting part out: what is left is the systematic alone.
|
|
out.sigma_excess = std::sqrt(std::max(0.0, sigma_sys * sigma_sys - sigma_stat * sigma_stat));
|
|
if (sigma_stat > 0.0)
|
|
out.forbidden_z = sigma_sys / sigma_stat;
|
|
out.n_observations = static_cast<int>(cells.entry_cell.size());
|
|
out.ok = true;
|
|
return out;
|
|
}
|
|
}
|
|
|
|
std::vector<AnisotropyObservation> ScaledObservations(const std::vector<IntegrationOutcome> &outcomes,
|
|
bool rotation, const gemmi::SpaceGroup *space_group,
|
|
double min_partiality) {
|
|
// Per-image scale, indexed the way the outcomes are.
|
|
std::vector<double> g(outcomes.size(), 0.0);
|
|
for (size_t i = 0; i < outcomes.size(); ++i)
|
|
if (outcomes[i].image_scale_g.has_value() && std::isfinite(*outcomes[i].image_scale_g)
|
|
&& *outcomes[i].image_scale_g > 0.0f)
|
|
g[i] = *outcomes[i].image_scale_g;
|
|
|
|
std::vector<AnisotropyObservation> out;
|
|
// The same acceptance the merge itself uses (Merge.cpp): everything in an IntegrationOutcome came
|
|
// out of the integrator, and Reflection::observed is not set by the _process.h5 reader, so testing
|
|
// it would silently empty this on the --mode scale path.
|
|
auto usable = [](const Reflection &r) {
|
|
return std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f
|
|
&& std::isfinite(r.rlp) && r.rlp > 0.0f && std::isfinite(r.partiality);
|
|
};
|
|
|
|
if (!rotation) {
|
|
// A still is a whole measurement of its reflection, and consecutive stills are different
|
|
// crystals, so there is nothing to assemble.
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
if (!(g[i] > 0.0))
|
|
continue;
|
|
for (const auto &r : outcomes[i].reflections) {
|
|
if (!usable(r) || !(r.partiality >= min_partiality))
|
|
continue;
|
|
const double corr = r.rlp / (r.partiality * g[i]);
|
|
out.push_back({r.h, r.k, r.l, static_cast<float>(r.I * corr),
|
|
static_cast<float>(r.sigma * corr), r.d});
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// A rotation reflection is integrated image by image, so it arrives as a run of partials over
|
|
// consecutive frames. Assemble each run into one full the way the 3D combine and the unmerged
|
|
// export do - same raw hkl, frames no further apart than the combine's own gap, parts added with
|
|
// their variances in quadrature - because a single partial divided by its own partiality carries
|
|
// the rocking-curve model's error as well as its intensity, and that error is a function of the
|
|
// reflection's direction relative to the spindle, which is precisely the direction this
|
|
// diagnostic is measuring.
|
|
constexpr float MAX_FRAME_GAP = 2.0f;
|
|
// The sort key travels with the part instead of being read back through the pointer, the way the
|
|
// merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
|
|
// and an indirect compare is a cache miss on every one of them. The keys are the same values in
|
|
// the same order, so introsort makes the same comparisons and the same swaps and leaves the same
|
|
// order - which matters, because two parts can genuinely share (h,k,l) and image_number and the
|
|
// event sums below are floating point.
|
|
struct Part {
|
|
int32_t h, k, l;
|
|
float image_number;
|
|
const Reflection *r;
|
|
size_t image;
|
|
};
|
|
// Assemble a SAMPLE of the unique reflections, not all of them. Everything below exists to feed
|
|
// MeasureSystematicFloor, which strides whole unique reflections down to MAX_CLUSTERS and on a
|
|
// large sweep keeps about one in seventy of what it is handed; taking the sample here instead of
|
|
// there removes the sort rather than shrinking it. Two things decide what the sample may be:
|
|
// * whole unique reflections. The floor is measured from the contrast BETWEEN a reflection's
|
|
// symmetry mates, so a reflection is taken with all of its mates or not at all - which is why
|
|
// this needs the space group the floor will be measured in, and samples on the ASU key.
|
|
// * coverage in direction and in resolution, which is what the fit's 12 shells x 60 directions
|
|
// read. A hash of the ASU key is a uniform sample of unique reflections and carries both; a
|
|
// stride over the parts, or over hkl order, carries neither.
|
|
// How many unique reflections the partials cover is not known here, and does not need to be: one
|
|
// part in 2^n is one unique reflection in 2^n whatever their number, so the rate comes from the
|
|
// partial count, which is.
|
|
size_t n_parts = 0;
|
|
for (size_t i = 0; i < outcomes.size(); ++i)
|
|
if (g[i] > 0.0)
|
|
n_parts += outcomes[i].reflections.size();
|
|
const HKLKeyGenerator asu_key(/*merge_friedel=*/true,
|
|
space_group ? *space_group : *gemmi::find_spacegroup_by_number(1));
|
|
const uint64_t sample_mask =
|
|
std::bit_ceil((n_parts + MAX_SAMPLE_PARTS - 1) / MAX_SAMPLE_PARTS) - 1;
|
|
std::vector<Part> parts;
|
|
// Reserve what the sample holds rather than letting the vector grow: the growth chain copies a few
|
|
// hundred megabytes on a large sweep.
|
|
parts.reserve(std::min(n_parts, MAX_SAMPLE_PARTS));
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
if (!(g[i] > 0.0))
|
|
continue;
|
|
for (const auto &r : outcomes[i].reflections) {
|
|
if (!usable(r))
|
|
continue;
|
|
if ((MixKey(asu_key(r.h, r.k, r.l).pack()) & sample_mask) != 0)
|
|
continue;
|
|
parts.push_back({r.h, r.k, r.l, r.image_number, &r, i});
|
|
}
|
|
}
|
|
std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) {
|
|
return std::tie(a.h, a.k, a.l, a.image_number) < std::tie(b.h, b.k, b.l, b.image_number);
|
|
});
|
|
for (size_t i = 0; i < parts.size();) {
|
|
size_t j = i + 1;
|
|
while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k
|
|
&& parts[j].l == parts[i].l
|
|
&& parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP)
|
|
++j;
|
|
double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, p_g = 0.0;
|
|
for (size_t m = i; m < j; ++m) {
|
|
const Reflection &r = *parts[m].r;
|
|
sum_p += r.partiality;
|
|
sum_I += static_cast<double>(r.I) * r.rlp;
|
|
sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
|
|
p_g += r.partiality * g[parts[m].image];
|
|
}
|
|
const Reflection &first = *parts[i].r;
|
|
i = j;
|
|
// The scale is the rocking curve's own weighted mean over the frames the event spans; the
|
|
// partiality divisor only corrects an event the sweep cut short, since a complete one sums
|
|
// to 1 by construction.
|
|
if (!(sum_p >= min_partiality) || !(p_g > 0.0))
|
|
continue;
|
|
const double scale = p_g / sum_p; // partiality-weighted mean G
|
|
const double corr = 1.0 / (sum_p * scale);
|
|
out.push_back({first.h, first.k, first.l, static_cast<float>(sum_I * corr),
|
|
static_cast<float>(std::sqrt(sum_var) * corr), first.d});
|
|
}
|
|
return out;
|
|
}
|
|
|
|
const char *AnisotropyShapeCode(AnisotropyShape shape) {
|
|
switch (shape) {
|
|
case AnisotropyShape::Linear: return "LINEAR";
|
|
case AnisotropyShape::Flat: return "FLAT";
|
|
case AnisotropyShape::Convex: return "CONVEX";
|
|
default: return "UNDETERMINED";
|
|
}
|
|
}
|
|
|
|
const char *AnisotropyVerdictCode(AnisotropyVerdict verdict) {
|
|
switch (verdict) {
|
|
case AnisotropyVerdict::NotDetected: return "NOT_DETECTED";
|
|
case AnisotropyVerdict::Detected: return "DETECTED";
|
|
default: return "CANNOT_DETERMINE";
|
|
}
|
|
}
|
|
|
|
AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
|
|
const std::vector<AnisotropyObservation> &unmerged,
|
|
const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup *space_group,
|
|
const AnisotropyRunInfo &run) {
|
|
AnisotropyResult result;
|
|
if (merged.empty() || !cell.is_crystal())
|
|
return result;
|
|
|
|
const auto allowed = AllowedDeviatoricBasis(cell, space_group);
|
|
const auto full = FullDeviatoricBasis();
|
|
const auto forbidden = ForbiddenBasis(full, allowed);
|
|
const int p = static_cast<int>(allowed.size());
|
|
const int q = static_cast<int>(forbidden.size());
|
|
result.n_free_parameters = p;
|
|
|
|
if (p == 0) {
|
|
// Cubic. The anisotropy tensor is forced isotropic by symmetry, so its deviatoric part is exactly
|
|
// zero - there is nothing to fit and nothing to gate.
|
|
result.n_reflections = static_cast<int>(merged.size());
|
|
result.delta_b = 0.0;
|
|
result.delta_b_linear = 0.0;
|
|
result.delta_b_flat = 0.0;
|
|
result.eigenvalue[0] = result.eigenvalue[1] = result.eigenvalue[2] = 0.0;
|
|
result.verdict = AnisotropyVerdict::NotDetected;
|
|
result.band = "no anisotropy is possible in this Laue class";
|
|
return result;
|
|
}
|
|
|
|
double isig_sum = 0.0;
|
|
int isig_n = 0;
|
|
for (const auto &r : merged)
|
|
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
|
|
isig_sum += r.I / r.sigma;
|
|
++isig_n;
|
|
}
|
|
const double mean_isig = isig_n > 0 ? isig_sum / isig_n : 0.0;
|
|
|
|
const auto refl = SelectReflections(merged, cell, space_group, N_SHELLS);
|
|
if (refl.size() < static_cast<size_t>(MIN_CELLS * MIN_CELL_REFLECTIONS))
|
|
return result;
|
|
|
|
const Cells cells = BuildCells(refl, cell, space_group, N_SHELLS);
|
|
const TensorFit constrained = FitTensor(cells, allowed, true);
|
|
if (!constrained.ok)
|
|
return result;
|
|
result.n_reflections = static_cast<int>(refl.size());
|
|
result.n_cells = static_cast<int>(cells.mu.size());
|
|
|
|
const Eigen3 axes = Decompose(constrained.b);
|
|
result.delta_b = axes.value[0] - axes.value[2];
|
|
for (int n = 0; n < 3; ++n) {
|
|
result.eigenvalue[n] = axes.value[n];
|
|
for (int j = 0; j < 3; ++j)
|
|
result.eigenvector[n][j] = axes.vec[n][j];
|
|
}
|
|
double d_min = 0.0;
|
|
for (const auto &r : merged)
|
|
if (r.d > 0.0f && (d_min == 0.0 || r.d < d_min))
|
|
d_min = r.d;
|
|
if (d_min > 0.0)
|
|
result.fold_weakening = std::exp(result.delta_b / (2.0 * d_min * d_min));
|
|
ConeLimits(merged, cell, space_group, axes, N_SHELLS, result.d_min_axis, result.d_min_censored);
|
|
{
|
|
double lo = INFINITY, hi = -INFINITY;
|
|
for (double d : result.d_min_axis)
|
|
if (std::isfinite(d)) { lo = std::min(lo, d); hi = std::max(hi, d); }
|
|
if (std::isfinite(lo) && std::isfinite(hi))
|
|
result.d_min_spread = hi - lo;
|
|
}
|
|
|
|
// --- the resolution signature, and whether the binning decides it ---
|
|
const ShapeFit shape = FitShape(cells, constrained.b);
|
|
// A fit that did not run was fitted over no shells. Reported that way rather than as the shells
|
|
// it had available, so the renderer and the cautions below do not print NaN coefficients beside a
|
|
// shell count that suggests they mean something.
|
|
result.shape_shells = shape.ok ? shape.shells : 0;
|
|
result.shape_intercept = shape.c0;
|
|
result.shape_intercept_z = shape.c0_err > 0.0 ? shape.c0 / shape.c0_err : NAN;
|
|
result.shape_slope = shape.c1;
|
|
result.shape_curvature_z = shape.curvature_err > 0.0 ? shape.curvature / shape.curvature_err : NAN;
|
|
result.shape_curvature_share = shape.curvature_share;
|
|
result.shape_residual = shape.residual;
|
|
if (shape.ok && shape.shells >= MIN_SHAPE_SHELLS) {
|
|
result.shape = shape.flat ? AnisotropyShape::Flat
|
|
: shape.convex ? AnisotropyShape::Convex : AnisotropyShape::Linear;
|
|
if (std::isfinite(shape.slope_through_origin) && shape.slope_through_origin != 0.0) {
|
|
result.delta_b_linear = result.delta_b * shape.c1 / shape.slope_through_origin;
|
|
result.delta_b_flat = result.delta_b - result.delta_b_linear;
|
|
}
|
|
// Rebin at 8 and at 16 shells and refit the tensor each time, so the direction the amplitude is
|
|
// measured along is re-derived too. A verdict that moves is not a measurement.
|
|
result.shape_stable = true;
|
|
for (int alt : {8, 16}) {
|
|
const auto alt_refl = SelectReflections(merged, cell, space_group, alt);
|
|
const Cells alt_cells = BuildCells(alt_refl, cell, space_group, alt);
|
|
const TensorFit alt_fit = FitTensor(alt_cells, allowed, false);
|
|
if (!alt_fit.ok) {
|
|
result.shape_stable = false;
|
|
break;
|
|
}
|
|
const ShapeFit alt_shape = FitShape(alt_cells, alt_fit.b);
|
|
if (!alt_shape.ok || alt_shape.shells < MIN_SHAPE_SHELLS || alt_shape.flat != shape.flat) {
|
|
result.shape_stable = false;
|
|
break;
|
|
}
|
|
}
|
|
if (!result.shape_stable)
|
|
result.shape = AnisotropyShape::Undetermined;
|
|
}
|
|
// The gate is applied to the part of delta_b that behaves as a Debye-Waller B, when that could be
|
|
// separated; to delta_b itself when it could not.
|
|
// A FLAT crystal can fit a negative linear slope; a negative anisotropy is not a smaller one, it is
|
|
// no established Debye-Waller component at all.
|
|
const double magnitude = std::max(0.0, std::isfinite(result.delta_b_linear) ? result.delta_b_linear
|
|
: result.delta_b);
|
|
|
|
// --- the gate ---
|
|
auto refuse = [&](std::string why) {
|
|
result.verdict = AnisotropyVerdict::CannotDetermine;
|
|
result.refusal = std::move(why);
|
|
};
|
|
if (mean_isig < MIN_MEAN_I_OVER_SIGMA) {
|
|
refuse("the merged data are at the noise floor (<I/sigma(I)> below 1), where neither the tensor "
|
|
"nor its resolution signature means anything");
|
|
} else if (q == 0) {
|
|
// Triclinic. Every other Laue class leaves directions in which the true tensor is zero by
|
|
// symmetry, and measuring those is how this dataset's systematic error is established. Triclinic
|
|
// leaves none, and three candidate substitutes (a degree-4 spherical-harmonic block, a
|
|
// low/high-resolution split, a floor borrowed from similar datasets) were built and measured
|
|
// against known truth: the first two are blind to real anisotropy, as a denominator must be, but
|
|
// neither delivers a usable scale, and the borrowed floor spans three orders of magnitude.
|
|
refuse("the Laue class is triclinic, which leaves no symmetry-forbidden tensor direction, so this "
|
|
"dataset carries no measurement of its own systematic error");
|
|
} else if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < MIN_ROTATION_DEG) {
|
|
refuse("the observed rotation range is too short for the tensor to be separated from a lab-fixed "
|
|
"systematic");
|
|
} else if (std::isfinite(run.observed_rotation_deg) && !run.dose_term_in_scale_model) {
|
|
// A dose ramp is not a significance failure - it is a confident false detection. A simulated
|
|
// isotropic crystal given a 25 A^2 relative-B ramp over 180 deg returns dB ~ 9 A^2 with a
|
|
// significance no test can see through, because the ramp really is a smooth quadratic function of
|
|
// direction once the sweep geometry is folded in.
|
|
refuse("the scale model carried no dose term, and an uncorrected dose ramp manufactures anisotropy "
|
|
"that no significance test can see through");
|
|
} else if (unmerged.empty()) {
|
|
// Merged intensities have exact Laue symmetry by construction, so the symmetry-forbidden tensor
|
|
// directions - the only place a dataset can measure its own systematic error - are identically
|
|
// zero there whatever the crystal carries. Gating on what is left would be gating on counting
|
|
// statistics, which against known ground truth calls "real" in 15-58% of clean isotropic datasets.
|
|
refuse("the systematic-error scale can only be measured on unmerged observations, and none were "
|
|
"available for this merge");
|
|
} else {
|
|
const SystematicFloor sys = MeasureSystematicFloor(unmerged, cell, space_group, full, forbidden);
|
|
result.n_observations = sys.n_observations;
|
|
result.forbidden_z = sys.forbidden_z;
|
|
if (!sys.ok) {
|
|
refuse("the symmetry-forbidden tensor directions of the unmerged observations could not be "
|
|
"measured, so this dataset's systematic-error scale is unknown");
|
|
} else {
|
|
const double k_p = K_ALLOWED[p];
|
|
result.sigma_systematic = sys.sigma_excess;
|
|
const double allowed_null = CountingNull(constrained.cov, k_p);
|
|
result.floor = std::sqrt(k_p * result.sigma_systematic * k_p * result.sigma_systematic
|
|
+ allowed_null * allowed_null);
|
|
if (result.floor > 0.0) {
|
|
result.significance = magnitude / result.floor;
|
|
result.detection_limit = GATE_ESTABLISHED * result.floor;
|
|
result.band = Band(result.significance);
|
|
result.verdict = result.significance > GATE_ESTABLISHED ? AnisotropyVerdict::Detected
|
|
: AnisotropyVerdict::NotDetected;
|
|
} else {
|
|
refuse("the systematic-error floor could not be measured");
|
|
}
|
|
}
|
|
}
|
|
|
|
// --- cautions: measured things that bias the verdict, in either direction ---
|
|
if (result.shape == AnisotropyShape::Flat)
|
|
result.cautions.emplace_back(
|
|
"the directional deficit does not follow exp(-1/2 s^T B s), so the fitted deltaB is a fit "
|
|
"of the wrong functional form and may be an under-estimate");
|
|
if (result.shape == AnisotropyShape::Undetermined && result.shape_shells >= MIN_SHAPE_SHELLS)
|
|
result.cautions.emplace_back(
|
|
"the resolution signature changed when the shells were rebinned, so it is reported as "
|
|
"undetermined rather than as a measurement");
|
|
if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < 180.0)
|
|
result.cautions.emplace_back(
|
|
"the observed rotation range is under 180 deg, which leaves a second tensor direction "
|
|
"reachable by a lab-fixed systematic and reduces the power of the test");
|
|
if (std::isfinite(run.radiation_damage_relative_b) && std::fabs(run.radiation_damage_relative_b) > 10.0)
|
|
result.cautions.emplace_back(
|
|
"the run carries substantial radiation damage; the scale model removes its average "
|
|
"monotone part, but a non-monotone dose ramp needs --relative-b to be taken out too");
|
|
// R2's Laue-class caution, fired on its measured signature rather than on the Laue class alone.
|
|
// Too high a symmetry assignment pushes real anisotropy into the directions the systematic error is
|
|
// measured in, so sigma_sys is inflated by the signal itself and the gate collapses (measured: a
|
|
// 42x inflation of sigma_sys, a gate of 106 falling to 3.1). The signature is a forbidden subspace
|
|
// far above its own counting noise TOGETHER with a gate that has established nothing. p == 1 on its
|
|
// own is not the signature - it is true of every tetragonal, trigonal and hexagonal dataset, 17 of
|
|
// 38 in the battery, most of them correctly assigned.
|
|
if (p == 1 && result.verdict == AnisotropyVerdict::NotDetected
|
|
&& std::isfinite(result.forbidden_z) && result.forbidden_z > FORBIDDEN_Z_SUSPECT)
|
|
result.cautions.emplace_back(
|
|
"the space group is worth checking before this verdict is believed: the tensor directions "
|
|
"this Laue class forbids carry far more scatter than their own counting noise, yet nothing "
|
|
"was established in the single direction it allows - which is what a too-high symmetry "
|
|
"assignment looks like, real anisotropy being measured as systematic error");
|
|
return result;
|
|
}
|
|
|
|
std::string AnisotropyToText(const AnisotropyResult &result) {
|
|
std::ostringstream os;
|
|
if (result.n_reflections == 0)
|
|
return os.str();
|
|
os << std::fixed;
|
|
os << "Diffraction anisotropy\n";
|
|
if (result.n_free_parameters == 0) {
|
|
os << " The Laue class forces the anisotropy tensor to be isotropic: deltaB = 0 exactly, with no\n"
|
|
<< " free parameter to fit. This is symmetry, not a measurement.\n";
|
|
return os.str();
|
|
}
|
|
os << std::setprecision(2);
|
|
os << " Anisotropic deltaB (range of principal components) = " << result.delta_b << " A^2"
|
|
<< " [" << result.n_free_parameters << " free direction"
|
|
<< (result.n_free_parameters == 1 ? "" : "s") << " in this Laue class]\n";
|
|
os << " Principal components (A^2, relative):";
|
|
for (double v : result.eigenvalue)
|
|
os << " " << (v - result.eigenvalue[2]);
|
|
os << "\n";
|
|
if (std::isfinite(result.fold_weakening))
|
|
os << " Strongest / weakest direction at the resolution limit: " << result.fold_weakening << "x\n";
|
|
// Every axis, not just the outer two: ConeLimits fills each independently, so a crystal whose
|
|
// first and third directions fall in cones too sparse to measure - but whose second is well
|
|
// covered - would otherwise suppress the whole line, d_min_spread included.
|
|
if (std::isfinite(result.d_min_axis[0]) || std::isfinite(result.d_min_axis[1])
|
|
|| std::isfinite(result.d_min_axis[2])) {
|
|
bool any_censored = false;
|
|
os << " d_min along the principal directions (<I/sigma(I)> = " << CONE_I_OVER_SIGMA
|
|
<< " in a " << CONE_HALF_ANGLE_DEG << " deg cone):";
|
|
for (int n = 0; n < 3; ++n) {
|
|
if (!std::isfinite(result.d_min_axis[n])) {
|
|
os << " -";
|
|
} else if (result.d_min_censored[n]) {
|
|
os << " <" << result.d_min_axis[n];
|
|
any_censored = true;
|
|
} else {
|
|
os << " " << result.d_min_axis[n];
|
|
}
|
|
}
|
|
os << " A\n";
|
|
if (any_censored)
|
|
os << " \"<\" marks a direction in which <I/sigma(I)> had not fallen to "
|
|
<< CONE_I_OVER_SIGMA << " where the measured data\n"
|
|
<< " end: that figure is the edge of the data, not this crystal's limit along that\n"
|
|
<< " direction, and the spread between the directions is therefore a lower bound.\n";
|
|
}
|
|
if (result.shape_shells >= MIN_SHAPE_SHELLS) {
|
|
os << " Resolution signature of the deficit over " << result.shape_shells << " shells: "
|
|
<< AnisotropyShapeCode(result.shape) << " (intercept " << std::showpos << result.shape_intercept
|
|
<< std::noshowpos << ", z = " << result.shape_intercept_z << "; slope " << result.shape_slope
|
|
<< " A^2; through-origin chi2/dof " << result.shape_residual << ")\n";
|
|
if (result.shape == AnisotropyShape::Flat)
|
|
os << " The deficit does not follow exp(-1/2 s^T B s): " << result.delta_b_flat
|
|
<< " A^2 of the deltaB above does not behave as a Debye-Waller B. The fall-off is being\n"
|
|
<< " described by the wrong functional form, so the number may be an UNDER-estimate,\n"
|
|
<< " not an over-estimate.\n";
|
|
else if (result.shape == AnisotropyShape::Convex)
|
|
os << " The deficit grows faster than s^2, which a Debye-Waller B cannot do.\n";
|
|
}
|
|
if (result.verdict == AnisotropyVerdict::CannotDetermine) {
|
|
os << " => CANNOT DETERMINE: " << result.refusal << ".\n";
|
|
} else {
|
|
// The magnitude judged is the Debye-Waller part of the deltaB, which is what the floor is a
|
|
// floor on; on a crystal whose deficit is not a Debye-Waller fall-off the two differ. Which of
|
|
// the two numbers is which is said here rather than left to the reader: deltaB_linear exceeds
|
|
// the deltaB on roughly a third of a real battery, and a verdict line that quoted the larger
|
|
// one without saying so read as a contradiction of the headline above.
|
|
const bool linear = std::isfinite(result.delta_b_linear);
|
|
const double judged = std::max(0.0, linear ? result.delta_b_linear : result.delta_b);
|
|
os << std::setprecision(2)
|
|
<< (result.verdict == AnisotropyVerdict::Detected ? " => DETECTED (" : " => NOT DETECTED (")
|
|
<< result.band << "): ratio " << result.significance << " = "
|
|
<< (linear ? "deltaB_linear " : "deltaB ") << judged
|
|
<< " A^2 / this dataset's own systematic-error floor of " << result.floor << " A^2.\n";
|
|
if (linear) {
|
|
os << " The gate is applied to deltaB_linear, the part of the deltaB that follows\n"
|
|
<< " exp(-1/2 s^T B s)";
|
|
if (result.delta_b_linear > result.delta_b)
|
|
os << " - it exceeds the deltaB here because the fitted line passes below\n"
|
|
<< " the origin";
|
|
else if (result.delta_b_linear < 0.0)
|
|
os << " - the fitted linear slope is negative here, which is no established\n"
|
|
<< " Debye-Waller component at all rather than a small one, so it enters as zero";
|
|
os << ". The anisotropy to quote is the deltaB of " << result.delta_b << " A^2 above.\n";
|
|
}
|
|
}
|
|
if (std::isfinite(result.detection_limit))
|
|
os << " Below about " << result.detection_limit << " A^2 nothing could be established on this\n"
|
|
<< " dataset. That limit is set by the systematic error, not by counting, so it does NOT\n"
|
|
<< " improve with more reflections or a longer exposure.\n";
|
|
for (const auto &c : result.cautions)
|
|
os << " Note: " << c << "\n";
|
|
os << " Reported and not corrected: no intensity is changed and no reflection is removed.\n";
|
|
return os.str();
|
|
}
|