leonarski_fandClaude Opus 5
0e6ece3c38
BuildObservationCells strides whole ASU clusters down to MAX_CLUSTERS = 20000 and keeps 1.4% of what it is handed - so the assembly before it sorted and summed every integrated partial in the sweep to throw 98.6% of the result away on the next line. The comment there says the striding "costs nothing and bounds the work"; it bounded it after the expensive part. Sample first instead: keep a partial iff a splitmix64 hash of its ASU key falls in the low one part in 2^n. Whole ASU clusters are mandatory rather than convenient - the cluster-robust sandwich sums one score per cluster, and the forbidden-direction signal IS the contrast between mates, so a raw-hkl sample would collapse that onto the independent covariance and inflate the floor. The hash rather than the key matters too: HKLKey::pack puts h in the low bits, so masking the key directly would select lattice planes instead of a sample. The rate comes from the partial count, not from merged.size(): the merged list is resolution-cut, and taking the rate from it made the first version never fire at all on the crystal it was written for. One part in 2^n of the partials is one unique reflection in 2^n, so the count already computed for the reserve is the right one. DELTA_B is unchanged by construction, not by luck. AnalyzeAnisotropy touches the unmerged observations in exactly two places - an empty test and MeasureSystematicFloor - so the tensor, the shape verdict and the directional limits come from the merged list alone. Over the battery: 0 verdict flips of 38; DELTA_B, DELTA_B_LINEAR, SHAPE, FREE_DIRECTIONS and the D_MIN keys differ on none of the 38, byte-identical rather than within noise; every reported key identical on 26; the mmCIF diff is exactly _reflns.jfjoch_aniso_floor and _significance; no caution appears or disappears anywhere. The spread those two keys move within was measured rather than assumed, and re-running cannot measure it because the code is deterministic. The shipped code already samples - by striding cluster order - so moving the phase of its own stride gives the sampling spread it has always carried. Against that control the change moves the floor less: ratio range 0.803-1.481 against 0.497-1.343, geomean 1.031 against 0.941. The gate therefore did not get easier to trip - the floor's geomean moves +3.1% here against -5.9% for the phase control, the crystal nearest the gate is not sampled at all and is bit-identical, and the two sampled crystals nearest the gate both move away from it. N_OBSERVATIONS goes up, not down, because a smaller pool lands the downstream integer stride slightly above the cluster cap. The ASU reduction is 14-45 ns a call in isolation, which nearly sank the design; in place it is free, because it overlaps the memory stalls of the Reflection stream it walks. The 24 crystals below the sampling cap run the shipped algorithm plus the reduction and cost +1.79 s in total, median +0.04 s. The block that introduced this measures 8.07 s before and 1.97 s after, which is 5.74 s off the largest crystal (-12.3%) and -11.9 s of anisotropy time over the battery. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
1464 lines
73 KiB
C++
1464 lines
73 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 {
|
|
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;
|
|
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);
|
|
result.shape_shells = shape.shells;
|
|
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.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.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 && 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";
|
|
if (std::isfinite(result.d_min_axis[0]) || 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();
|
|
}
|