diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 7a5c5a5e7..4180ce21b 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -197,9 +197,15 @@ P. Evans, "Scaling and assessment of data quality" (2006), Acta Cryst. D62, 72-8 **Amplitudes from intensities** — the posterior-mean amplitude of each merged intensity under the acentric and centric Wilson priors is French and Wilson's; giving no amplitude to an intensity more than 3.7 sigma below zero, and leaving such intensities out of the prior, follows CCP4's -[ctruncate](https://www.ccp4.ac.uk/). S. French and K. Wilson, "On the treatment of negative intensity +[ctruncate](https://www.ccp4.ac.uk/) (C. Ballard and N. Stein). So does scaling each reflection's +Wilson prior by the anisotropy tensor along its direction, which ctruncate has done by default since +its version 1.7 ("use anisotropy in prior for truncate procedure"); rugnux uses its own tensor for it +(see Diffraction anisotropy below). S. French and K. Wilson, "On the treatment of negative intensity observations" (1978), Acta Cryst. A34, 517-525 -[doi:10.1107/S0567739478001114](https://doi.org/10.1107/S0567739478001114). +[doi:10.1107/S0567739478001114](https://doi.org/10.1107/S0567739478001114); ctruncate is cited through +the CCP4 suite: M. D. Winn, C. C. Ballard, K. D. Cowtan et al., +"Overview of the CCP4 suite and current developments" (2011), Acta Cryst. D67, 235-242 +[doi:10.1107/S0907444910045749](https://doi.org/10.1107/S0907444910045749). **Absorption as spherical harmonics** — describing an empirical absorption correction as a series of real spherical harmonics of the beam directions in the crystal frame is Blessing's; its use as a diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a64928c70..2d2102e54 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -6,6 +6,7 @@ * Rugnux: `.hkl` now holds unmerged, scaled full reflections (SHELX HKLF 4) on rotation data, so SHELXL computes Rint itself. * Rugnux reads Rigaku d*TREK SMV images (Saturn CCD), including detector 2theta, image orientation and encoded pixel overflows. * Rugnux integrates spots that grow wider than the integration disk away from the beam (typical of small molecules at high X-ray energy) over their measured footprint. +* Rugnux: French-Wilson amplitudes (`F`/`SIGF`) use an anisotropic Wilson prior from the fitted anisotropy tensor, as ctruncate does; intensities are unchanged. ### 1.0.0-rc.173 diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 147c15f8d..0242995d9 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -50,7 +50,7 @@ The methods draw on, and in places reimplement, solutions from: - R. J. Read, P. D. Adams & A. J. McCoy, "Intensity statistics in the presence of translational noncrystallographic symmetry", *Acta Cryst.* **D69** (2013), 176-183 (the native-Patterson detection of translational pseudo-symmetry, and the intensity modulation it produces, which the axial-zone screw-absence test scores against). - A. Barty, R. A. Kirian, F. R. N. C. Maia et al., "Cheetah: software for high-throughput reduction and analysis of serial femtosecond X-ray diffraction data", *J. Appl. Cryst.* **47** (2014), 1118-1131 (peakfinder8: the per-resolution-ring background statistics of §3.2). - A. Hennequin, B. Couturier, V. V. Gligorov & L. Lacassagne, "SparseCCL: Connected Components Labeling and Analysis for sparse images", DASIP 2019, 65-70 (the connected-component labelling of §3.4, used via ACTS/traccc). -- S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities). +- S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities), and CCP4's ctruncate (C. Ballard & N. Stein), whose anisotropic Wilson prior §10.8 follows, cited through M. D. Winn et al., "Overview of the CCP4 suite and current developments", *Acta Cryst.* **D67** (2011), 235-242. - A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free cross-validation). - R. A. Fisher, "Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population", *Biometrika* **10** (1915), 507-521 (the z-transformation on which the correction surfaces' held-out half-set CC1/2 is compared). - M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index c3586f883..dc01bfad4 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -76,7 +76,7 @@ By default the reported/written high-resolution limit is trimmed where $\mathrm{ ### 13.5 Diffraction anisotropy -How fast the intensity falls off with resolution can depend on direction. Rugnux measures that, reports it, and does nothing else with it: no intensity is corrected, no reflection is removed on a directional criterion, and the merged data and the written files do not depend on direction at all. +How fast the intensity falls off with resolution can depend on direction. Rugnux measures that and reports it. No intensity is corrected and no reflection is removed on a directional criterion; the one use of the tensor is the Wilson prior of the French–Wilson amplitudes (§10.8), so `F`/`SIGF` follow the fall-off along each reflection's direction while the intensities do not depend on direction at all. **The tensor.** A deviatoric anisotropic displacement tensor is fitted to the merged intensities as diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 4a55123ab..de38ee633 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -106,7 +106,9 @@ Only the ring moves. The signal disk $r_1$ stays circular, deliberately: it sets What a circular $r_1$ loses is flux, and that loss is **not** a function of resolution alone: measured per reflection, it carries a directional component worth several Ų with a definite principal axis, on top of the isotropic part. Nor is there anything in the merge to absorb it. There is **no per-shell scale**, and there cannot usefully be one: every scale in §10 is fitted against a reference built from a reflection's own symmetry equivalents, and equivalents share $s^2$ exactly, so any function of $s^2$ lies in the exact null space of the whole scaling model — a per-shell parameter would have zero residual to fit against. (XDS and DIALS have the same null space, for the same reason.) The isotropic part of the loss is instead degenerate with the overall Wilson $B$ and is silently reported as part of it, so **the reported `WILSON_B` / `_reflns.B_iso_Wilson_estimate` carries an $r_1$-dependent contribution**: measured across a constant-ring-area radius sweep it falls monotonically as the disk grows, by 0.5 Ų on sharp strong data and by up to ~10 Ų on weak wide-spot data. What this costs the *data* is much less than what it costs the flux, because most of the loss is matched by a proportional $\sigma$: it moves no CC$_{1/2}$ and no $R_\text{meas}$, and — to within a few hundredths of an ångström — no resolution cut. -**Measured spot footprint (automatic).** The radii above are chosen from spots at 5 Å, which at high X-ray energy sit close to the beam. Away from it a spot can grow several times wider — radially from the sensor's parallax and the obliquity of the incidence, tangentially from the crystal's azimuthal spread, which rotates the diffracted beam about the incident one and smears the spot along its ring. On small-molecule data at 20–25 keV the standard deviation grows from ~1 px near the beam to ~5 px at the detector edge: the $r_1 = 4$ disk holds a quarter of the flux there, the $6\ldots13$ px ring a third of it, and the profile widths learned inside $r_1$ (§9.3) saturate near $r_1^2/4$. So the pre-scan measures every spot it finds with a window that follows the spot — three of its own standard deviations, iterated and re-centred — separately along and across the radius, and tabulates the median widths $\sigma_ ho,\sigma_ au$ against the distance from the beam. Wherever $3\max(\sigma_ ho,\sigma_ au)>r_1$ the integrator then (i) starts the background ring at $3\sigma$ along each axis, (ii) sums the reflection over the $r_1$ disk **and** the $3\sigma$ footprint ellipse, so the summation — the profile fit's seed and its fallback — holds the spot rather than its core, and (iii) builds the per-reflection Gaussian at the measured widths on a grid grown to hold them. Where every spot fits the disk nothing is installed and the integration is unchanged bit for bit, which is the case for compact protein spots; like the measured radius, the footprint applies to the canonical pass and not to the geometry pre-pass, and a canonical pass whose wider rings the neighbours starve falls back to the settings without it. Judged by refining the published structures with SHELXL, it removes the intensity loss that grew with resolution on the small-molecule sets (rugnux/model intensity in the outermost shell 0.81–0.91 → 0.98–1.02). +**Measured spot footprint (automatic).** The radii above are chosen from spots at 5 Å, which at high X-ray energy sit close to the beam. Away from it a spot can grow several times wider — radially from the sensor's parallax and the obliquity of the incidence, tangentially from the crystal's azimuthal spread, which rotates the diffracted beam about the incident one and smears the spot along its ring. On small-molecule data at 20–25 keV the standard deviation grows from ~1 px near the beam to ~5 px at the detector edge: the $r_1 = 4$ disk holds a quarter of the flux there, the $6\ldots13$ px ring a third of it, and the profile widths learned inside $r_1$ (§9.3) saturate near $r_1^2/4$. So the pre-scan measures every spot it finds with a window that follows the spot — three of its own standard deviations, iterated and re-centred — separately along and across the radius, and tabulates the median widths $\sigma_ +ho,\sigma_ au$ against the distance from the beam. Wherever $3\max(\sigma_ +ho,\sigma_ au)>r_1$ the integrator then (i) starts the background ring at $3\sigma$ along each axis, (ii) sums the reflection over the $r_1$ disk **and** the $3\sigma$ footprint ellipse, so the summation — the profile fit's seed and its fallback — holds the spot rather than its core, and (iii) builds the per-reflection Gaussian at the measured widths on a grid grown to hold them. Where every spot fits the disk nothing is installed and the integration is unchanged bit for bit, which is the case for compact protein spots; like the measured radius, the footprint applies to the canonical pass and not to the geometry pre-pass, and a canonical pass whose wider rings the neighbours starve falls back to the settings without it. Judged by refining the published structures with SHELXL, it removes the intensity loss that grew with resolution on the small-molecule sets (rugnux/model intensity in the outermost shell 0.81–0.91 → 0.98–1.02). ### 9.2 Box summation (seed and fallback) @@ -430,7 +432,7 @@ $ \sigma_F = \sqrt{\langle J\rangle - \langle|F|\rangle^2}. $ -The prior mean is $\Sigma = \varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell}$, where $\varepsilon$ is the reflection's epsilon (symmetry-enhancement) multiplicity and $\langle I/\varepsilon\rangle$ is the Wilson mean in its resolution shell (so reflections on symmetry elements, and each shell, are treated correctly); a shell whose mean is not positive takes the mean of the nearest lower-resolution shell. As in `ctruncate`, an intensity more than 3.7σ below zero gets no amplitude (the intensity is kept) and does not enter the shell mean. Strong reflections ($I>20\sigma$) short-circuit to $|F|=\sqrt{I}$, where the French–Wilson bias is below 0.3%; a reflection with an unusable $I/\sigma$ falls back to $\sqrt{\max(I,0)}$. The integral is evaluated numerically with a log-shift for stability. +The prior mean is $\Sigma = \varepsilon\,K_\mathrm{shell}\,a(\mathbf{h})$, where $\varepsilon$ is the reflection's epsilon (symmetry-enhancement) multiplicity, $a(\mathbf{h}) = \exp(-\tfrac12\mathbf{s}^\mathsf{T}B\,\mathbf{s})$ carries the deviatoric anisotropy tensor $B$ of §13.5 along the reflection's own direction, and $K_\mathrm{shell} = \sum I/\varepsilon \,/ \sum a$ over its resolution shell, so the priors of a shell still average to its measured Wilson mean. The amplitudes are first made with $a = 1$ at the end of the merge and made again once the tensor has been fitted; only $F$/$\sigma_F$ change, never an intensity. With an isotropic prior the weak direction of an anisotropic crystal gets a prior set mostly by the strong direction, which turns its noise into amplitude; for isotropic data the tensor is near zero and the prior reduces to the shell mean, so it is used whenever a tensor was fitted. A shell whose $K$ is not positive takes that of the nearest lower-resolution shell. As in `ctruncate`, an intensity more than 3.7σ below zero gets no amplitude (the intensity is kept) and does not enter the shell mean. Strong reflections ($I>20\sigma$) short-circuit to $|F|=\sqrt{I}$, where the French–Wilson bias is below 0.3%; a reflection with an unusable $I/\sigma$ falls back to $\sqrt{\max(I,0)}$. The integral is evaluated numerically with a log-shift for stability. Amplitudes are written as MTZ `F`/`SIGF`, mmCIF `_refln.F_meas_au`/`F_meas_sigma_au`, and appended to the text HKL, alongside the intensity columns. The **same** $|F|$ feed the model-validation step (§14), so the reflection file and the maps use one consistent set of amplitudes. diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index 15813d1dc..51abe2bdb 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -1409,6 +1409,18 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector &merged, return result; } +gemmi::Mat33 AnisotropyTensorHKL(const AnisotropyResult &result, const gemmi::UnitCell &cell) { + if (!std::isfinite(result.eigenvalue[0])) + return gemmi::Mat33(0.0); + gemmi::Mat33 b(0.0); + for (int n = 0; n < 3; ++n) + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + b[i][j] += result.eigenvalue[n] * result.eigenvector[n][i] * result.eigenvector[n][j]; + // s = h^T frac, as everywhere above, so s^T B s = h^T (frac B frac^T) h. + return cell.frac.mat.multiply(b).multiply(cell.frac.mat.transpose()); +} + std::string AnisotropyToText(const AnisotropyResult &result) { std::ostringstream os; if (result.n_reflections == 0) diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.h b/image_analysis/scale_merge/AnisotropyAnalysis.h index 61f5573ca..47021eb69 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.h +++ b/image_analysis/scale_merge/AnisotropyAnalysis.h @@ -11,8 +11,9 @@ #include "gemmi/symmetry.hpp" #include "gemmi/unitcell.hpp" -// Diffraction-anisotropy diagnostic. It DESCRIBES and it REPORTS: no intensity is corrected, no -// reflection is removed, and the written reflection file does not depend on direction in any way. +// Diffraction-anisotropy diagnostic. It DESCRIBES and it REPORTS: no intensity is corrected and no +// reflection is removed. The one use of the tensor is the Wilson prior of the French-Wilson +// amplitudes (AnisotropyTensorHKL), so F follows direction while the intensities do not. // Two quantities are reported, and they are not the same thing: // // - the anisotropic deltaB, the range of the principal components of the anisotropy tensor (the @@ -172,3 +173,7 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector &merged, const AnisotropyRunInfo &run = {}); std::string AnisotropyToText(const AnisotropyResult &result); + +// The fitted tensor in the Miller-index basis of `cell`: Q with s^T B s = h^T Q h, the form +// FrenchWilsonOptions::anisotropy_hkl takes. Zero (the isotropic prior) when no tensor was fitted. +gemmi::Mat33 AnisotropyTensorHKL(const AnisotropyResult &result, const gemmi::UnitCell &cell); diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp index 632890e44..c50c987d2 100644 --- a/image_analysis/scale_merge/FrenchWilson.cpp +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -10,6 +10,7 @@ #include #include "../../common/ResolutionShells.h" +#include "gemmi/math.hpp" #include "gemmi/symmetry.hpp" namespace { @@ -125,15 +126,28 @@ void ApplyFrenchWilson(std::vector &merged, const gemmi::Space return; } - // Wilson mean intensity per resolution shell. + // Wilson mean intensity per resolution shell, and with an anisotropy tensor the + // direction it is expected along. The prior of reflection h is then epsilon * K * a(h), with + // a(h) = exp(-1/2 h^T Q h) - the form the tensor was fitted in - and the shell constant + // K = sum(I/epsilon) / sum(a) over the shell, so the priors of a shell still average to its + // measured mean. A shell's weak direction gets a prior matched to its own intensities instead of + // one set mostly by the strong direction, which turned the weak direction's noise into amplitude. + // Isotropic data give a near-zero tensor, a(h) ~ 1 and the isotropic prior back, so the tensor is + // used whenever one was fitted, with no detection gate in front of it: the prior needs the best + // estimate of along h, whether the anisotropy comes from the crystal or not. + // Following the anisotropic prior of CCP4 ctruncate (Ballard & Stein; Winn et al. (2011) Acta Cryst. D67, 235-242) ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, opts.num_shells); - std::vector shell_sum(opts.num_shells, 0.0); + std::vector shell_sum(opts.num_shells, 0.0), shell_norm(opts.num_shells, 0.0); std::vector shell_count(opts.num_shells, 0); - double global_sum = 0.0; + double global_sum = 0.0, global_norm = 0.0; int global_count = 0; auto epsilon = [&](const MergedReflection &r) { return std::max(1, gops.epsilon_factor_without_centering({{r.h, r.k, r.l}})); }; + auto anisotropy = [&](const MergedReflection &r) { + const gemmi::Vec3 h(r.h, r.k, r.l); + return std::exp(-0.5 * h.dot(opts.anisotropy_hkl.multiply(h))); + }; // An intensity that gets no amplitude (below reject_below, see fw_one) stays out of the prior too, // as ctruncate leaves its outliers out of the norm: kept in, the systematically negative // intensities of a background over-subtracted on a powder ring pull the shell mean down and with @@ -143,18 +157,21 @@ void ApplyFrenchWilson(std::vector &merged, const gemmi::Space || r.I < opts.reject_below * r.sigma) continue; const double i_over_eps = r.I / epsilon(r); + const double a = anisotropy(r); global_sum += i_over_eps; + global_norm += a; ++global_count; if (const auto s = shells.GetShell(r.d)) { shell_sum[*s] += i_over_eps; + shell_norm[*s] += a; ++shell_count[*s]; } } - const double global_mean = global_count > 0 ? std::max(global_sum / global_count, 1e-10) : 1.0; + const double global_mean = global_count > 0 ? std::max(global_sum / global_norm, 1e-10) : 1.0; std::vector shell_mean(opts.num_shells, global_mean); for (int s = 0; s < opts.num_shells; ++s) if (shell_count[s] >= opts.min_reflections_per_shell) - shell_mean[s] = shell_sum[s] / shell_count[s]; + shell_mean[s] = shell_sum[s] / shell_norm[s]; // A shell whose mean intensity is not positive has no measurable signal, but a prior of ~0 would // still take every amplitude in it to ~0 - more confidently than any measurement says. It takes // the nearest lower-resolution shell's mean instead (shell 0 is the lowest resolution), which is @@ -197,7 +214,7 @@ void ApplyFrenchWilson(std::vector &merged, const gemmi::Space for (int i = lo; i < hi; ++i) { MergedReflection &r = merged[i]; const auto s = shells.GetShell(r.d); - const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean); + const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean) * anisotropy(r); const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}}); fw_one(r, r.I, r.sigma, r.F, r.sigmaF, logw, grid, sigma_wilson, centric); fw_one(r, r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus, logw, grid, sigma_wilson, centric); diff --git a/image_analysis/scale_merge/FrenchWilson.h b/image_analysis/scale_merge/FrenchWilson.h index 5a25034ac..e767bc82a 100644 --- a/image_analysis/scale_merge/FrenchWilson.h +++ b/image_analysis/scale_merge/FrenchWilson.h @@ -7,6 +7,7 @@ #include #include "../../common/Reflection.h" +#include "gemmi/math.hpp" #include "gemmi/symmetry.hpp" struct FrenchWilsonOptions { @@ -16,12 +17,17 @@ struct FrenchWilsonOptions { double strong_cutoff = 20.0; // I/sigma above which <|F|> = sqrt(I) (FW bias negligible) double reject_below = -3.7; // I/sigma below which no amplitude is given (as ctruncate) int num_threads = 1; // workers for the per-reflection integration (independent per reflection) + // Anisotropic Wilson prior: the anisotropy tensor B (ln = c - 1/2 s^T B s, A^2) in the + // Miller-index basis, Q with s^T B s = h^T Q h (AnisotropyTensorHKL). The default, zero, is the + // isotropic prior. + gemmi::Mat33 anisotropy_hkl = gemmi::Mat33(0.0); }; // Fill F and sigmaF on every merged reflection with the French-Wilson estimate of the amplitude: // the posterior mean |F| given the measured intensity I and its sigma under the Wilson prior. The -// prior uses the resolution-shell mean intensity, the correct centric/acentric form, and the -// reflection's epsilon (symmetry-enhancement) multiplicity. Strong reflections reduce to sqrt(I); +// prior uses the resolution-shell mean intensity - scaled along the reflection's direction by the +// anisotropy tensor when one is given - the correct centric/acentric form, and the reflection's +// epsilon (symmetry-enhancement) multiplicity. Strong reflections reduce to sqrt(I); // reflections with an unusable I/sigma fall back to sqrt(max(I,0)) with propagated sigma. An // intensity below reject_below sigmas gets no amplitude (F = NaN) and stays out of the prior; the // intensity itself is kept. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index e0fc66f16..14e9848f4 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1932,19 +1932,20 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (an.d_min_censored[worst_i]) p += fmt::format("\n Along {} the crystal reaches AT LEAST {:.2f} A - the measured data" " end before\n the signal does.", along(worst), *worst); - p += "\n Nothing was corrected or removed: the merged data and the written files do not\n" - " depend on direction at all."; + p += "\n No intensity was corrected and no reflection removed; only the French-Wilson\n" + " amplitudes (F) take the direction into account, through their Wilson prior."; Add(s, Prose(p)); } else if (an.verdict == AnisotropyVerdict::Detected) { Add(s, Prose(fmt::format( - " Anisotropy: {} (deltaB {:.1f} A^2). Nothing was corrected or removed: the merged\n" - " data do not depend on direction.", severity, an.delta_b))); + " Anisotropy: {} (deltaB {:.1f} A^2). No intensity was corrected and no reflection\n" + " removed; only the French-Wilson amplitudes (F) take the direction into account.", + severity, an.delta_b))); } else if (an.verdict == AnisotropyVerdict::CannotDetermine) { Add(s, Prose(" Anisotropy: cannot be determined in this Laue class - there is no\n" " symmetry-forbidden direction to calibrate the fit's own noise against.")); } else { Add(s, Prose(" Anisotropy: below this data set's own noise floor. No directional statement can\n" - " be made, and the merged data do not depend on direction.")); + " be made, and no intensity was corrected for direction.")); } Add(s, Prose("\n How much the fall-off depends on direction, and whether that is established above this\n" " data set's own systematic error. ANISOTROPY_DELTA_B is the range of the principal\n" diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 74aa0c628..0b8c64974 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -72,6 +72,7 @@ #include "../image_analysis/lattice_search/LatticeSearch.h" #include "../image_analysis/lattice_search/LePageLattice.h" #include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/TranslationalNCS.h" #include "../image_analysis/scale_merge/HKLKey.h" @@ -8806,8 +8807,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", wilson.b, wilson.correlation, wilson.n_shells); - // Diffraction anisotropy. Report-only, like the twinning and Wilson analyses above: it - // describes the data and corrects nothing. The tensor and its resolution signature come from + // Diffraction anisotropy. It describes the data and corrects no intensity; its tensor is + // used only as the Wilson prior of the amplitudes, below. The tensor and its resolution signature come from // the merged intensities; the error bar the verdict is gated on has to come from the // unmerged observations, because a merge has exact Laue symmetry by construction and the // tensor directions the symmetry forbids - the only place a dataset measures its own @@ -8815,6 +8816,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (anisotropy.valid()) { sm.statistics.anisotropy = anisotropy.get(); stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n"; + // The amplitudes are the one place the tensor is used: French-Wilson again, with the + // Wilson prior of each reflection taken along its own direction. The intensities and + // everything measured on them above are left as they are. + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = config_.nthreads; + fw_opts.anisotropy_hkl = AnisotropyTensorHKL(sm.statistics.anisotropy, + gemmi::UnitCell(*result.consensus_cell)); + ApplyFrenchWilson(sm.merged, experiment_.GetSpaceGroupOrP1(), fw_opts); } } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3a8364894..c8732e70f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -49,6 +49,7 @@ #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/TranslationalNCS.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" @@ -1939,8 +1940,8 @@ static int RunRugnux(int argc, char **argv) { twin_sg ? twin_sg->centring_type() : 'P'); std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; - // Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - report-only, - // nothing is corrected. The stored reflections have just been re-scaled above, so they carry + // Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - no intensity + // is corrected, the tensor is only the amplitudes' Wilson prior. The stored reflections have just been re-scaled above, so they carry // this merge's own per-image scale. if (experiment.GetUnitCell()) { AnisotropyRunInfo aniso_run; @@ -1955,6 +1956,12 @@ static int RunRugnux(int argc, char **argv) { : 0.0f), *experiment.GetUnitCell(), twin_sg, aniso_run); std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl; + // The amplitudes again with the anisotropic Wilson prior, as the full pipeline does. + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = static_cast(nthreads); + fw_opts.anisotropy_hkl = AnisotropyTensorHKL(merged_statistics.anisotropy, + gemmi::UnitCell(*experiment.GetUnitCell())); + ApplyFrenchWilson(merged_reflections, experiment.GetSpaceGroupOrP1(), fw_opts); } // Before the reflection files, as in the full pipeline: the model settles the enantiomorph and, diff --git a/tests/FrenchWilsonTest.cpp b/tests/FrenchWilsonTest.cpp index da15e7296..ed498865b 100644 --- a/tests/FrenchWilsonTest.cpp +++ b/tests/FrenchWilsonTest.cpp @@ -6,6 +6,7 @@ #include #include +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" #include "../image_analysis/scale_merge/FrenchWilson.h" namespace { @@ -167,3 +168,51 @@ TEST_CASE("French-Wilson: a 5 sigma intensity is still pulled toward the prior", CHECK(v.back().F < 0.95f * std::sqrt(100.0f)); CHECK(v.back().F > std::sqrt(60.0f)); } + +TEST_CASE("French-Wilson: the anisotropic prior follows the tensor along each direction", "[french_wilson]") { + // A cubic 10 A cell with a deviatoric tensor of +20 A^2 along x and -10 A^2 along y and z: every + // ordinary reflection measures exactly 1000 * exp(-1/2 s^T B s), one shell. A weak probe along x + // must get the amplitude an isotropic prior at x's own gives, and less than the same + // measurement along y. + const gemmi::UnitCell cell(10, 10, 10, 90, 90, 90); + AnisotropyResult tensor; + const double eigenvalue[3] = {20.0, -10.0, -10.0}; + for (int n = 0; n < 3; ++n) { + tensor.eigenvalue[n] = eigenvalue[n]; + for (int j = 0; j < 3; ++j) + tensor.eigenvector[n][j] = n == j ? 1.0 : 0.0; + } + const gemmi::Mat33 q = AnisotropyTensorHKL(tensor, cell); + auto expected = [&](int h, int k, int l) { + const gemmi::Vec3 v(h, k, l); + return 1000.0 * std::exp(-0.5 * v.dot(q.multiply(v))); + }; + CHECK(expected(5, 0, 0) == Catch::Approx(1000.0 * std::exp(-0.5 * 20.0 * 0.25))); + + auto with_probes = [](std::vector v) { + v.push_back(Refl(5, 0, 0, 2.0f, 30.0f, 40.0f)); + v.push_back(Refl(0, 5, 0, 2.0f, 30.0f, 40.0f)); + return v; + }; + std::vector aniso_set; + for (int h = -8; h <= 8; ++h) + for (int k = -8; k <= 8; ++k) + for (int l = 1; l <= 8; ++l) + aniso_set.push_back(Refl(h, k, l, static_cast(cell.calculate_d({{h, k, l}})), + static_cast(expected(h, k, l)), 5.0f)); + std::vector iso_set; // the same reflections, all at the x probe's + for (const auto &r : aniso_set) + iso_set.push_back(Refl(r.h, r.k, r.l, r.d, static_cast(expected(5, 0, 0)), 5.0f)); + aniso_set = with_probes(aniso_set); + iso_set = with_probes(iso_set); + + FrenchWilsonOptions opts; + opts.num_shells = 1; + ApplyFrenchWilson(iso_set, SG(1), opts); + opts.anisotropy_hkl = q; + ApplyFrenchWilson(aniso_set, SG(1), opts); + const auto &along_x = aniso_set[aniso_set.size() - 2]; + const auto &along_y = aniso_set.back(); + CHECK(along_x.F == Catch::Approx(iso_set[iso_set.size() - 2].F).epsilon(0.01)); + CHECK(along_x.F < 0.9f * along_y.F); +}