From 2bae65c72e8d62044b0ff2c2a5a969e71a0180dc Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 31 Oct 2026 12:54:02 +0100 Subject: [PATCH] Stills partiality: each crystal refines its radial scale and rocking width, held to the run's spread The post-refinement now fits, per crystal and beside the orientation tilt, a radial scale eta of the reciprocal lattice (q -> (1 + eta) q: the photon energy of an XFEL shot, or a cell the spot positions did not pin) and a scale of the run's rocking width (the crystal's own mosaicity and size). Each changes the partialities differently across resolution, so neither is degenerate with the scale G. Both are held by Gaussian priors about the run whose widths come from the data: a first cycle refines against loose priors, and the spread of those estimates over the crystals, less the uncertainty each carries from its own reflections, is the prior of the second. A run whose crystals agree (a synthetic test, or a stable beam) is held to the run's curve; one whose crystals differ is not. On the quick panel (one run per set, merged MTZ): XFEL serial 18k frames, S anomalous peaks mean 6.80 -> 7.35 (weakest site 2.64 -> 1.65), CCref last shell 0.657 -> 0.653; XFEL serial 1.5k frames CCref last shell 0.222 -> 0.252; monochromatic jet CCref last 0.405 -> 0.458 (overall 0.447 -> 0.396); pink jet 0.299 -> 0.316; rotation sweeps as stills: lysozyme 5 keV peaks 3.36 -> 3.13 (CCref last 0.674 -> 0.689), thaumatin unchanged, a Zn site sweep: peaks 18.8 -> 20.6. On the full XFEL run CC_anom between independent halves 0.334 -> 0.39 and R-free 0.243 -> 0.235. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/ACKNOWLEDGEMENT.md | 10 + .../scale_merge/StillsPartialityRefine.cpp | 273 +++++++++++------- .../scale_merge/StillsPartialityRefine.h | 35 ++- 3 files changed, 209 insertions(+), 109 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 80e94cd33..9330bf040 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -78,6 +78,16 @@ high-throughput reduction and analysis of serial femtosecond X-ray diffraction d J. Appl. Cryst. 47, 1118-1131 [doi:10.1107/S1600576714007626](https://doi.org/10.1107/S1600576714007626). +**Per-crystal post-refinement of serial stills** - refining, for each still, the partiality +parameters the spot positions cannot pin (a small reorientation, the crystal's own rocking width and the +radial scale of its reciprocal lattice, which a shot's photon energy or a crystal's cell moves) against +the merged intensities follows the post-refinement of PRIME and the mosaic-block model of Sauter et al.: +M. Uervirojnangkoorn, O. B. Zeldin, A. Y. Lyubimov, J. Hattne, A. S. Brewster, N. K. Sauter, A. T. Brunger +and W. I. Weis, "Enabling X-ray free electron laser crystallography for challenging biological systems +from a limited number of crystals" (2015), eLife 4, e05421; N. K. Sauter, J. Hattne, A. S. Brewster, +N. Echols, P. H. Zwart and P. D. Adams, "Improved crystal orientation and physical properties from +single-shot XFEL stills" (2014), Acta Cryst. D70, 3299-3309 (DOIs not verified here, so not given). + Spot extraction groups strong pixels into spots with the sparse connected-component labelling of the ACTS traccc project: P. Gessinger, H. M. Gray, A. Krasznahorkay, C. Leggett, J. Niermann, A. Salzburger, S. N. Swatman and B. Yeo, "traccc: GPU track reconstruction library for HEP diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 6096da7f4..89d5413db 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -25,6 +25,10 @@ namespace { // Model error of the per-crystal scale fit, as a fraction of the predicted intensity - ScaleOnTheFly's. constexpr double SCALE_MODEL_ERROR = 0.3; + // The smallest priors on a crystal's radial scale and log width scale (RunPrior). + constexpr double MIN_PRIOR_ETA = 1e-5; + constexpr double MIN_PRIOR_LN_WIDTH = 0.01; + double SafeInv(double x, double fallback) { if (!std::isfinite(x) || x == 0.0) return fallback; @@ -79,30 +83,6 @@ namespace { } }; - // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y). - // Mirrors BraggPrediction: dist_ewald = |S| - 1/lambda with S = q_rot + S0, and - // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = gamma0^2 + (gamma_e*d*)^2 + (bw*|q_z|)^2 - three - // independent broadenings added in quadrature: the reciprocal-lattice point's own radius (gamma0, - // ~1/domain size, resolution-INdependent), the mosaic/divergence spread (gamma_e*d*, proportional to - // d*) and the bandwidth smear along the beam. In practice the fit below returns gamma_e ~ 0 and the - // width is essentially gamma0 - see there. - double ComputeP(double qx, double qy, double qz, - double psi_x, double psi_y, - double s0x, double s0y, double s0z, double inv_lambda, - const Width &width) { - double aa[3] = {psi_x, psi_y, 0.0}; - double q[3] = {qx, qy, qz}; - double qr[3]; - ceres::AngleAxisRotatePoint(aa, q, qr); - const double Sx = qr[0] + s0x, Sy = qr[1] + s0y, Sz = qr[2] + s0z; - const double de = std::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - inv_lambda; - const double dstar = std::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); - const double sig2 = width.Sigma2(dstar, qr[2]); - if (!(sig2 > 0.0)) - return 1.0; - return std::exp(-0.5 * de * de / sig2); - } - // The prediction window cuts the rocking curve. Reflections fill a window of half-width W evenly, so // their plain moment of de is u = W^2 / 3, and the intensity-weighted moment m of a Gaussian of width // s cut at +-W is u * g(a), a = W / s, with g(a) = 3/a^2 (1 - 2 a phi(a) / erf(a / sqrt 2)): 1 for a @@ -170,79 +150,118 @@ namespace { -s * ky, s * kx, c}; } - // ComputeP with the tilt given as its rotation matrix. - struct TiltModel { + // A crystal's model of its partialities: the orientation tilt (psi_x, psi_y) about the two axes + // perpendicular to the beam, a radial scale eta of its reciprocal lattice (q -> (1 + eta) q: a + // wavelength or a cell a little off for this crystal - on an XFEL the photon energy of its shot) and a + // scale exp(ln_width) of the run's rocking width (a crystal's own mosaicity or domain size). + // Per-crystal mosaicity and lattice scale in post-refinement, following Uervirojnangkoorn et al. (2015) + // eLife 4, e05421 and Sauter et al. (2014) Acta Cryst. D70, 3299-3309. + struct CrystalModel { + double psi_x = 0.0, psi_y = 0.0, eta = 0.0, ln_width = 0.0; + }; + + // The partiality of a reflection whose base reciprocal vector is q under a crystal model. Mirrors + // BraggPrediction: dist_ewald = |S| - 1/lambda with S = q' + S0, q' = (1 + eta) R(psi) q, and + // p = exp(-dist_ewald^2 / 2 sigma^2), sigma = exp(ln_width) sigma_run(d*). + + struct PartialityModel { double s0x, s0y, s0z, inv_lambda; Width width; - double P(const std::array &R, const FitObs &o) const { - const double qx = R[0] * o.qx + R[1] * o.qy + R[2] * o.qz; - const double qy = R[3] * o.qx + R[4] * o.qy + R[5] * o.qz; - const double qz = R[6] * o.qx + R[7] * o.qy + R[8] * o.qz; + double At(const std::array &R, const CrystalModel &m, double qx0, double qy0, double qz0) const { + const double f = 1.0 + m.eta; + const double qx = f * (R[0] * qx0 + R[1] * qy0 + R[2] * qz0); + const double qy = f * (R[3] * qx0 + R[4] * qy0 + R[5] * qz0); + const double qz = f * (R[6] * qx0 + R[7] * qy0 + R[8] * qz0); const double de = VecLen(qx + s0x, qy + s0y, qz + s0z) - inv_lambda; - const double sig2 = width.Sigma2(VecLen(qx, qy, qz), qz); + const double sig2 = width.Sigma2(VecLen(qx, qy, qz), qz) * std::exp(2.0 * m.ln_width); if (!(sig2 > 0.0)) return 1.0; return std::exp(-0.5 * de * de / sig2); } + + double P(const std::array &R, const CrystalModel &m, const FitObs &o) const { + return At(R, m, o.qx, o.qy, o.qz); + } }; - // Refine the orientation tilt (psi_x, psi_y) holding the scale G fixed: minimise - // sum_i rho(r_i^2) / 2 + |psi / sigma_prior|^2 / 2, r_i = weight_i (G p_i(psi) Iref_i / corr_i - Iobs_i), - // with the Cauchy rho(s) = k^2 log(1 + s / k^2) - the intensity-space objective of ScaleOnTheFly, and - // the Gaussian prior a proper MAP term since r_i is already in sigma units. Levenberg-Marquardt on the - // normal equations reweighted by rho', the Jacobian by central differences, every trial tilt held inside - // the bound. Returns false when the cost cannot be evaluated, leaving psi as it came in. - bool RefineTilt(const std::vector &obs, const TiltModel &model, double G, double robust_k, - double inv_sigma_prior, double max_tilt, double psi[2]) { + // Refine the crystal model holding the scale G fixed: minimise + // sum_i rho(r_i^2) / 2 + sum_k (x_k / sigma_k)^2 / 2, r_i = weight_i (G p_i(x) Iref_i / corr_i - Iobs_i), + // over x = (psi_x, psi_y, eta, ln_width), with the Cauchy rho(s) = k^2 log(1 + s / k^2) - the intensity- + // space objective of ScaleOnTheFly - and Gaussian priors about the run (x = 0), proper MAP terms since + // r_i is already in sigma units. Levenberg-Marquardt on the normal equations reweighted by rho', the + // Jacobian by central differences, every trial held inside the bounds. Returns false when the cost + // cannot be evaluated, leaving the model as it came in. + bool RefineCrystal(const std::vector &obs, const PartialityModel &model, double G, double robust_k, + const std::array &inv_sigma_prior, const std::array &bound, + CrystalModel &crystal, std::array &variance) { constexpr int MAX_ITERATIONS = 25; - constexpr double STEP = 1e-7; // rad, for the central differences + constexpr std::array STEP = {1e-7, 1e-7, 1e-7, 1e-5}; // for the central differences const double k2 = robust_k * robust_k; - const double prior2 = inv_sigma_prior * inv_sigma_prior; + auto as_model = [](const Eigen::Vector4d &x) { return CrystalModel{x(0), x(1), x(2), x(3)}; }; - auto cost = [&](const double p[2]) { - const auto R = TiltMatrix(p[0], p[1]); - double c = 0.5 * prior2 * (p[0] * p[0] + p[1] * p[1]); + auto cost = [&](const Eigen::Vector4d &x) { + const auto R = TiltMatrix(x(0), x(1)); + const CrystalModel m = as_model(x); + double c = 0.0; + for (int k = 0; k < 4; ++k) + c += 0.5 * inv_sigma_prior[k] * inv_sigma_prior[k] * x(k) * x(k); for (const auto &o: obs) { - const double r = o.weight * (G * model.P(R, o) * o.inv_corr * o.Iref - o.Iobs); + const double r = o.weight * (G * model.P(R, m, o) * o.inv_corr * o.Iref - o.Iobs); c += 0.5 * k2 * std::log1p(r * r / k2); } return c; }; - double current = cost(psi); + Eigen::Vector4d x(crystal.psi_x, crystal.psi_y, crystal.eta, crystal.ln_width); + double current = cost(x); if (!std::isfinite(current)) return false; double lambda = 1e-4; + Eigen::Matrix4d data_A = Eigen::Matrix4d::Identity(); // the normal matrix of the data alone for (int iter = 0; iter < MAX_ITERATIONS; ++iter) { - const auto R = TiltMatrix(psi[0], psi[1]); - const auto Rxp = TiltMatrix(psi[0] + STEP, psi[1]), Rxm = TiltMatrix(psi[0] - STEP, psi[1]); - const auto Ryp = TiltMatrix(psi[0], psi[1] + STEP), Rym = TiltMatrix(psi[0], psi[1] - STEP); - double A00 = prior2, A01 = 0.0, A11 = prior2; - double b0 = prior2 * psi[0], b1 = prior2 * psi[1]; + Eigen::Matrix4d A = Eigen::Matrix4d::Zero(); + Eigen::Vector4d b = Eigen::Vector4d::Zero(); + for (int k = 0; k < 4; ++k) { + A(k, k) = inv_sigma_prior[k] * inv_sigma_prior[k]; + b(k) = A(k, k) * x(k); + } + // The model at x and at x +- one step in each parameter. + std::array, 9> R; + std::array m; + for (int j = 0; j < 9; ++j) { + Eigen::Vector4d xj = x; + if (j > 0) + xj((j - 1) / 2) += (j % 2 ? 1.0 : -1.0) * STEP[(j - 1) / 2]; + R[j] = TiltMatrix(xj(0), xj(1)); + m[j] = as_model(xj); + } for (const auto &o: obs) { const double a = o.weight * G * o.inv_corr * o.Iref; - const double r = a * model.P(R, o) - o.weight * o.Iobs; - const double jx = a * (model.P(Rxp, o) - model.P(Rxm, o)) / (2.0 * STEP); - const double jy = a * (model.P(Ryp, o) - model.P(Rym, o)) / (2.0 * STEP); + const double r = a * model.P(R[0], m[0], o) - o.weight * o.Iobs; + Eigen::Vector4d J; + for (int k = 0; k < 4; ++k) + J(k) = a * (model.P(R[2 * k + 1], m[2 * k + 1], o) - model.P(R[2 * k + 2], m[2 * k + 2], o)) + / (2.0 * STEP[k]); const double w = 1.0 / (1.0 + r * r / k2); // rho' - A00 += w * jx * jx; - A01 += w * jx * jy; - A11 += w * jy * jy; - b0 += w * jx * r; - b1 += w * jy * r; + A += w * J * J.transpose(); + b += w * r * J; } + data_A = A; + for (int k = 0; k < 4; ++k) + data_A(k, k) -= inv_sigma_prior[k] * inv_sigma_prior[k]; bool accepted = false; - double trial[2] = {psi[0], psi[1]}; + Eigen::Vector4d trial = x; double trial_cost = current; while (lambda < 1e10) { - const double M00 = A00 * (1.0 + lambda), M11 = A11 * (1.0 + lambda); - const double det = M00 * M11 - A01 * A01; - if (!(det > 0.0)) + Eigen::Matrix4d M = A; + M.diagonal() *= 1.0 + lambda; + const Eigen::Vector4d dx = M.ldlt().solve(b); + if (!dx.allFinite()) break; - trial[0] = std::clamp(psi[0] - (M11 * b0 - A01 * b1) / det, -max_tilt, max_tilt); - trial[1] = std::clamp(psi[1] - (M00 * b1 - A01 * b0) / det, -max_tilt, max_tilt); + for (int k = 0; k < 4; ++k) + trial(k) = std::clamp(x(k) - dx(k), -bound[k], bound[k]); trial_cost = cost(trial); if (std::isfinite(trial_cost) && trial_cost < current) { accepted = true; @@ -253,14 +272,17 @@ namespace { } if (!accepted) break; - const double step = std::hypot(trial[0] - psi[0], trial[1] - psi[1]); const double decrease = current - trial_cost; - psi[0] = trial[0]; - psi[1] = trial[1]; + const double step = (trial - x).norm(); + x = trial; current = trial_cost; if (decrease <= 1e-6 * current || step <= 1e-10) break; } + crystal = as_model(x); + const Eigen::Matrix4d cov = data_A.inverse(); + for (int k = 0; k < 4; ++k) + variance[k] = cov(k, k); return true; } } @@ -273,12 +295,12 @@ StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) min_partiality_(x.GetScalingSettings().GetMinPartiality()), bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} -double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, +StillsPartialityRefine::CrystalFit StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, const std::map &reference, const std::optional &run_width, - double radial_scale) const { + double radial_scale, const CrystalPrior &prior) const { if (outcome.reflections.empty()) - return 0.0; + return {}; const Coord Astar = outcome.latt.Astar(); const Coord Bstar = outcome.latt.Bstar(); @@ -337,7 +359,7 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, } if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0) - return 0.0; + return {}; // Without the run's measured width (too few reflections in the run to measure it), the width comes // from this crystal's own reflections. Solve the 2x2 normal equations for de^2 = A + B d*^2. A degenerate spread in d* (all reflections in @@ -372,44 +394,47 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, .table_dstar = run_width ? &run_width->table_dstar : nullptr, .table_s = run_width ? &run_width->table_s : nullptr}; - double psi[2] = {0.0, 0.0}; + CrystalModel crystal; double G = 1.0; - const TiltModel model{.s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda, - .width = width}; - const bool refine_tilt = obs.size() >= settings_.min_reflections; - const int inner = refine_tilt ? settings_.inner_iterations : 1; + const PartialityModel model{.s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda, .width = width}; + const bool refine_crystal = obs.size() >= settings_.min_reflections; + const int inner = refine_crystal ? settings_.inner_iterations : 1; + const std::array inv_sigma_prior = {kRadToDeg / settings_.prior_sigma_deg, + kRadToDeg / settings_.prior_sigma_deg, + 1.0 / prior.sigma_eta, 1.0 / prior.sigma_ln_width}; + std::array variance{NAN, NAN, NAN, NAN}; + const std::array bound = {settings_.max_tilt_rad, settings_.max_tilt_rad, + settings_.max_eta, settings_.max_ln_width}; for (int it = 0; it < inner; ++it) { - // (1) Solve G given the current partialities, and weight the tilt fit by the variances at that model. + // (1) Solve G given the current partialities, and weight the model fit by the variances at that model. + const auto R = TiltMatrix(crystal.psi_x, crystal.psi_y); std::vector coeff(obs.size()); - for (size_t j = 0; j < obs.size(); ++j) { - const double p = ComputeP(obs[j].qx, obs[j].qy, obs[j].qz, psi[0], psi[1], - S0.x, S0.y, S0.z, inv_lambda, width); - coeff[j] = p * obs[j].inv_corr * obs[j].Iref; - } + for (size_t j = 0; j < obs.size(); ++j) + coeff[j] = model.P(R, crystal, obs[j]) * obs[j].inv_corr * obs[j].Iref; G = SolveScaleIRLS(obs, coeff, settings_.robust_k); if (!(G > 0.0) || !std::isfinite(G)) - return 0.0; + return {}; for (size_t j = 0; j < obs.size(); ++j) obs[j].weight = 1.0 / std::sqrt(ModelVariance(obs[j], G * coeff[j])); - if (!refine_tilt) + if (!refine_crystal) break; - // (2) Refine the tilt holding G fixed. A cost that cannot be evaluated keeps the tilt this crystal - // came in with and stops refining it - G alone is still a usable model. - const double inv_sigma_prior = settings_.prior_sigma_deg > 0.0 ? kRadToDeg / settings_.prior_sigma_deg : 0.0; - if (!RefineTilt(obs, model, G, settings_.robust_k, inv_sigma_prior, settings_.max_tilt_rad, psi)) + // (2) Refine the crystal model holding G fixed. A cost that cannot be evaluated keeps the model this + // crystal came in with and stops refining it - G alone is still a usable model. + if (!RefineCrystal(obs, model, G, settings_.robust_k, inv_sigma_prior, bound, crystal, variance)) break; } // Write the refined partiality + scale correction onto every reflection of the crystal (not only the // fit subset), so the merge sees a consistent model. image_scale_corr = prescaling_corr * qe_corr * flight_corr / // (partiality * G), the same composition ScaleOnTheFly writes. + const auto R = TiltMatrix(crystal.psi_x, crystal.psi_y); for (auto &r: outcome.reflections) { const Coord q = base_q(r); - const double p = ComputeP(q.x, q.y, q.z, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda, width); + const double p = model.At(R, crystal, q.x, q.y, q.z); r.partiality = static_cast(p); const float corr = r.prescaling_corr * r.qe_corr * r.flight_corr; const double denom = p * G; @@ -428,8 +453,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, outcome.image_scale_cc = cc; outcome.image_scale_cc_n = cc_n; - const double tilt_deg = std::sqrt(psi[0] * psi[0] + psi[1] * psi[1]) * kRadToDeg; - return tilt_deg; + return {.refined = refine_crystal, .tilt_deg = std::hypot(crystal.psi_x, crystal.psi_y) * kRadToDeg, + .eta = crystal.eta, .ln_width = crystal.ln_width, .var_eta = variance[2], .var_ln_width = variance[3]}; } double StillsPartialityRefine::MeasureRadialScale(const std::vector &outcomes, @@ -718,6 +743,40 @@ std::optional StillsPartialityRefine::MeasureR return width; } +StillsPartialityRefine::CrystalPrior StillsPartialityRefine::RunPrior(const std::vector &fits) const { + // The spread of a parameter over the crystals: the variance of the estimates less the median variance + // they carry from their own data, floored at a small value (a crystal is never held exactly). + auto spread = [&](auto value, auto variance, double floor) { + std::vector x, v; + for (const auto &f: fits) + if (f.refined && std::isfinite(variance(f)) && variance(f) > 0.0) { + x.push_back(value(f)); + v.push_back(variance(f)); + } + if (x.size() < 10) + return floor; + double mean = 0.0; + for (double e: x) + mean += e; + mean /= static_cast(x.size()); + double var = 0.0; + for (double e: x) + var += (e - mean) * (e - mean); + var /= static_cast(x.size() - 1); + std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); + return std::max(floor, std::sqrt(std::max(0.0, var - v[v.size() / 2]))); + }; + const CrystalPrior prior{ + .sigma_eta = spread([](const CrystalFit &f) { return f.eta; }, + [](const CrystalFit &f) { return f.var_eta; }, MIN_PRIOR_ETA), + .sigma_ln_width = spread([](const CrystalFit &f) { return f.ln_width; }, + [](const CrystalFit &f) { return f.var_ln_width; }, MIN_PRIOR_LN_WIDTH)}; + Logger("StillsPartialityRefine").Info("Stills crystals differ in radial scale by {:.1e} and in rocking width " + "by {:.0f}% (sd over the crystals, net of their own uncertainty)", + prior.sigma_eta, 100.0 * prior.sigma_ln_width); + return prior; +} + double StillsPartialityRefine::Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell, const std::optional &given_width, @@ -732,7 +791,6 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si nthreads = AvailableCpus(); nthreads = std::max(1, nthreads); - double last_mean_tilt = 0.0; // Where the frames were predicted through a curve measured before the image pass, the partialities use // that same curve. Measured again here, through windows that follow it, the width reads narrower at low // resolution (its constant term falls to about zero), and the merge it gave moved R-free against the @@ -750,6 +808,13 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si if (used_width) *used_width = run_width; + // The priors on a crystal's radial scale and width are the run's own spread of them: the first cycle + // refines every crystal against loose priors, and the spread of those estimates less their own + // uncertainty (method of moments) is how much crystals really differ. On an XFEL each shot has its + // photon energy and the spread is large; on a monochromatic or a stable pink beam it comes out small, + // and the crystals are held to the run's curve. + CrystalPrior prior{.sigma_eta = settings_.max_eta, .sigma_ln_width = settings_.max_ln_width}; + std::vector fits(outcomes.size()); for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. const std::vector merged = MergeAll(experiment_, outcomes, nthreads, reference_cell); @@ -757,19 +822,12 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si for (const auto &m: merged) reference[hkl_key_generator_(m)] = m.I; - std::atomic tilt_sum{0.0}; - std::atomic tilt_n{0}; std::atomic next{0}; - auto worker = [&]() { size_t i = next.fetch_add(1); while (i < outcomes.size()) { - const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width, radial_scale) : 0.0; - if (t > 0.0) { - double prev = tilt_sum.load(); - while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {} - tilt_n.fetch_add(1); - } + fits[i] = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width, radial_scale, prior) + : CrystalFit{}; i = next.fetch_add(1); } }; @@ -786,8 +844,17 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si f.get(); } - last_mean_tilt = tilt_n > 0 ? tilt_sum.load() / static_cast(tilt_n.load()) : 0.0; + if (outer == 0) + prior = RunPrior(fits); } + double tilt_sum = 0.0; + size_t tilt_n = 0; + for (const auto &f: fits) + if (f.refined) { + tilt_sum += f.tilt_deg; + ++tilt_n; + } + const double last_mean_tilt = tilt_n > 0 ? tilt_sum / static_cast(tilt_n) : 0.0; return last_mean_tilt; } diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 5b99d6d36..6f644fb8c 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -23,10 +23,18 @@ // reshapes the SPATIAL pattern of partialities in a way a single G cannot mimic. That breaks the // degeneracy that killed the scalar-sigma fit. // +// Two more parameters per crystal, each held by a prior about the run: a radial scale eta of its +// reciprocal lattice (q -> (1 + eta) q, which moves the excitation error by about eta lambda d*^2 / 2 - +// on an XFEL the photon energy of the crystal's own shot, which spot positions cannot pin) and a scale of +// the run's rocking width (the crystal's own mosaicity and size). Both change p differently across +// resolution, so neither is degenerate with G. Their priors are not constants: the first cycle refines +// against loose ones, and the run's spread of the estimates, net of their own uncertainty, becomes the +// prior of the next (RunPrior) - large on SASE XFEL data, small on a monochromatic or stable pink beam. +// // Because the integrated intensity I_obs is fixed, refining the tilt only recomputes the partiality // analytically from the stored per-image lattice/geometry (q = A*.h + B*.k + C*.l, then // dist_ewald = |q + S0| - 1/lambda, signed) - NO pixel re-integration. The loop is: merge -> per-crystal -// refine tilt (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and +// refine tilt, eta and width (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and // image_scale_corr -> re-merge, iterated a few times. Mutates each reflection's `partiality` and // `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections. The // width of the Gaussian partiality is the run's, measured once from the intensity-weighted excitation @@ -38,8 +46,10 @@ public: int inner_iterations = 3; // (solve G) <-> (refine tilt) alternations per crystal size_t min_reflections = 40; // skip tilt refinement below this (anti-overfit on sparse crystals) double max_tilt_rad = 0.0175; // hard bound on |dpsi| (~1 deg); indexing already refined orientation - double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0 (0 = off); tames weak-data overfit, + double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0; tames weak-data overfit, // inert on strong data (well-supported tilts overcome it) + double max_eta = 5e-3; // bound on a crystal's radial scale |eta|, and its prior in the first cycle + double max_ln_width = 1.1; // the same for the log of its width scale (a factor of 3) double robust_k = 3.0; // Cauchy loss scale (sigma units) double gamma_e = 0.0; // angular width sigma(d*) = gamma_e*d* (0 = measure from data) }; @@ -84,8 +94,21 @@ private: const double min_partiality_; const float bandwidth_sigma_; - // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped). - double RefineOne(IntegrationOutcome &outcome, - const std::map &reference, - const std::optional &run_width, double radial_scale) const; + // The priors of a crystal's radial scale eta and log width scale (Gaussian, about 0). + struct CrystalPrior { + double sigma_eta, sigma_ln_width; + }; + // One crystal's refined model: |dpsi| in degrees, eta, ln width and the variances of the last two from + // its own data (no prior). refined = false where only G (or nothing) was fitted. + struct CrystalFit { + bool refined = false; + double tilt_deg = 0.0, eta = 0.0, ln_width = 0.0, var_eta = NAN, var_ln_width = NAN; + }; + + // Refine one crystal against the reference map. + CrystalFit RefineOne(IntegrationOutcome &outcome, const std::map &reference, + const std::optional &run_width, double radial_scale, + const CrystalPrior &prior) const; + // The run's spread of eta and ln width over the crystals, net of their own uncertainty. + CrystalPrior RunPrior(const std::vector &fits) const; };