RESULT-CHANGING (owner decision): stills tilt post-refinement solved directly instead of by Ceres
Each crystal's (psi_x, psi_y) tilt was refined by a Ceres problem rebuilt for every inner iteration: one auto-differentiated residual block and one Cauchy loss per reflection. With the scaling reference cleaned up it was 70% of the stills scale/merge. The same objective - Cauchy on the sigma-weighted intensity residuals plus the Gaussian prior on the tilt, the tilt bounded - is now minimised by a small Levenberg-Marquardt on the rho'-reweighted 2x2 normal equations, the tilt applied to all reflections through one rotation matrix and the Jacobian taken by central differences. It reaches the same answer to optimizer precision: on the new synthetic test (StillsPartialityRefine_RecoversTheTilt) partiality error 0.0559 -> 0.03019 with Ceres and 0.03016 with this, mean tilt 0.02065 vs 0.02067 deg; on two 25k-frame XFEL runs (--mode scale) CC1/2, I/sigma, R_pim and the CC against a deposited lysozyme model agree to 0.0005, mean |dpsi| 0.023 deg on both. Scale+merge 6.9 -> 3.0 s (--mode scale, loaded box). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -10,6 +10,8 @@
|
||||
#include "../image_analysis/scale_merge/HKLKey.h"
|
||||
#include "../image_analysis/scale_merge/Merge.h"
|
||||
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
|
||||
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
|
||||
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
|
||||
#include "gemmi/reciproc.hpp"
|
||||
|
||||
TEST_CASE("HKLKey_NoSG_noMergeFriedel") {
|
||||
@@ -650,3 +652,97 @@ TEST_CASE("MergeOnTheFly_AddImagesIsTheSerialMerge") {
|
||||
CHECK(serial_stats.overall.cc_half == parallel_stats.overall.cc_half);
|
||||
CHECK(serial_stats.overall.r_meas == parallel_stats.overall.r_meas);
|
||||
}
|
||||
|
||||
// Stills partiality post-refinement on noiseless synthetic crystals: each crystal's true orientation is
|
||||
// tilted from the stored one by up to 0.05 deg, its intensities are the true partials of a fixed set of
|
||||
// full intensities, and its stored partialities those of the stored orientation. The refinement has to
|
||||
// move the partialities towards the true ones. (The Ceres solve this replaced gave the same numbers: error
|
||||
// 0.0559 -> 0.0302, mean tilt 0.0206 deg.)
|
||||
TEST_CASE("StillsPartialityRefine_RecoversTheTilt") {
|
||||
constexpr float CELL = 50.0f, WAVELENGTH = 1.0f, D_MIN = 2.0f;
|
||||
constexpr double GAMMA0 = 0.001, ACCEPT = 1.7 * GAMMA0; // RLP width; |dist_ewald| accepted at prediction
|
||||
DiffractionExperiment x;
|
||||
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
||||
DiffractionGeometry geom;
|
||||
geom.Wavelength_A(WAVELENGTH);
|
||||
const Coord S0 = geom.GetScatteringVector();
|
||||
|
||||
std::mt19937 rng(11);
|
||||
std::normal_distribution<double> gauss(0.0, 1.0);
|
||||
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
||||
auto full_I = [](int h, int k, int l) { // a fixed, Friedel-symmetric pseudo-random full intensity
|
||||
const uint64_t key = HKLKeyGenerator(true, gemmi::get_spacegroup_by_name("P 1"))(h, k, l).pack();
|
||||
return 10.0 + 990.0 * static_cast<double>((key * 0x9e3779b97f4a7c15ULL) >> 40) / static_cast<double>(1ULL << 24);
|
||||
};
|
||||
auto rotate = [](const double R[9], const Coord &v) {
|
||||
return Coord(static_cast<float>(R[0] * v.x + R[1] * v.y + R[2] * v.z),
|
||||
static_cast<float>(R[3] * v.x + R[4] * v.y + R[5] * v.z),
|
||||
static_cast<float>(R[6] * v.x + R[7] * v.y + R[8] * v.z));
|
||||
};
|
||||
auto rotation = [](double ax, double ay, double az, double angle, double R[9]) { // Rodrigues
|
||||
const double n = std::sqrt(ax * ax + ay * ay + az * az);
|
||||
ax /= n; ay /= n; az /= n;
|
||||
const double c = std::cos(angle), s = std::sin(angle), t = 1.0 - c;
|
||||
const double M[9] = {c + t * ax * ax, t * ax * ay - s * az, t * ax * az + s * ay,
|
||||
t * ax * ay + s * az, c + t * ay * ay, t * ay * az - s * ax,
|
||||
t * ax * az - s * ay, t * ay * az + s * ax, c + t * az * az};
|
||||
std::copy(M, M + 9, R);
|
||||
};
|
||||
auto partiality = [&](const Coord &q) {
|
||||
const double de = (q + S0).Length() - 1.0 / WAVELENGTH;
|
||||
return std::exp(-0.5 * de * de / (GAMMA0 * GAMMA0));
|
||||
};
|
||||
|
||||
std::vector<IntegrationOutcome> outcomes(150);
|
||||
std::vector<std::vector<double>> true_p(outcomes.size());
|
||||
const int hmax = static_cast<int>(CELL / D_MIN);
|
||||
for (size_t c = 0; c < outcomes.size(); ++c) {
|
||||
double R[9], T[9];
|
||||
rotation(gauss(rng), gauss(rng), gauss(rng), 2.0 * M_PI * unit(rng), R);
|
||||
const double tilt = (0.02 + 0.03 * unit(rng)) * M_PI / 180.0, phi = 2.0 * M_PI * unit(rng);
|
||||
rotation(std::cos(phi), std::sin(phi), 0.0, tilt, T); // about an axis normal to the beam
|
||||
auto &o = outcomes[c];
|
||||
o.geom = geom;
|
||||
o.latt = CrystalLattice(rotate(R, Coord(CELL, 0, 0)), rotate(R, Coord(0, CELL, 0)),
|
||||
rotate(R, Coord(0, 0, CELL)));
|
||||
const double G = 0.5 + unit(rng);
|
||||
for (int h = -hmax; h <= hmax; ++h)
|
||||
for (int k = -hmax; k <= hmax; ++k)
|
||||
for (int l = -hmax; l <= hmax; ++l) {
|
||||
const Coord q = o.latt.Astar() * h + o.latt.Bstar() * k + o.latt.Cstar() * l;
|
||||
const float dstar = q.Length();
|
||||
if (dstar == 0.0f || dstar > 1.0f / D_MIN)
|
||||
continue;
|
||||
if (std::fabs((q + S0).Length() - 1.0 / WAVELENGTH) > ACCEPT)
|
||||
continue;
|
||||
const double p_true = partiality(rotate(T, q));
|
||||
Reflection r{};
|
||||
r.h = h; r.k = k; r.l = l;
|
||||
r.d = 1.0f / dstar;
|
||||
r.I = static_cast<float>(G * p_true * full_I(h, k, l));
|
||||
r.sigma = static_cast<float>(1.0 + 0.01 * r.I);
|
||||
r.prescaling_corr = 1.0f;
|
||||
r.partiality = static_cast<float>(partiality(q));
|
||||
r.image_scale_corr = 1.0f / r.partiality;
|
||||
o.reflections.push_back(r);
|
||||
true_p[c].push_back(p_true);
|
||||
}
|
||||
REQUIRE(o.reflections.size() >= 40);
|
||||
}
|
||||
|
||||
auto mean_abs_error = [&]() {
|
||||
double sum = 0.0;
|
||||
size_t n = 0;
|
||||
for (size_t c = 0; c < outcomes.size(); ++c)
|
||||
for (size_t i = 0; i < outcomes[c].reflections.size(); ++i, ++n)
|
||||
sum += std::fabs(outcomes[c].reflections[i].partiality - true_p[c][i]);
|
||||
return sum / static_cast<double>(n);
|
||||
};
|
||||
const double before = mean_abs_error();
|
||||
ScaleOnTheFly(x, MergeAll(x, outcomes, 4, std::nullopt)).Scale(outcomes, 4);
|
||||
const double mean_tilt = StillsPartialityRefine(x).Run(outcomes, 4, std::nullopt);
|
||||
const double after = mean_abs_error();
|
||||
CHECK(mean_tilt > 0.015);
|
||||
CHECK(mean_tilt < 0.06);
|
||||
CHECK(after < 0.6 * before); // measured 0.54: the prior holds the larger tilts back
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user