mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 88 Wiki Activity

performance: the French-Wilson quadrature stops recomputing what its grid already fixes

Browse Source 
The posterior is integrated on j = (i + 1/2) * dj, so sqrt(j) and log(j) each
split into a per-reflection factor and a term that depends only on i:
sqrt(j) = sqrt(dj) * sqrt(i + 1/2) and log(j) = log(dj) + log(i + 1/2). The
i-dependent halves are the same for every reflection and every intensity, so
they are tabulated once for the run rather than recomputed at each of the 400
grid points. The exponential in the second pass is left alone, but a point more
than 37 below the peak carries a weight under 1e-16 and cannot move a
normalised sum of doubles, so it is skipped instead of exponentiated - and the
posterior is sharply peaked, so most of the grid is skipped. Finally the Wilson
prior and the centric flag belong to the reflection rather than to one of its
three intensities, so the two symmetry lookups are made once and shared instead
of three times over.

None of this changes the arithmetic: on five datasets spanning the corpus's
speed range the space group, the unique-reflection count, R_meas and CC1/2 are
identical to the digit.

Measured over those five, each run twice and the two binaries interleaved so
machine load cannot favour one arm: 771.2 s to 730.7 s, 5.3 % of total wall,
between 3.2 % and 7.8 % per dataset and in the same direction on all five. The
work removed is largest where the run merges most - the quadrature runs on
every weak intensity of every merge, and a large cell at high resolution reaches
a few million of them.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_fandClaude Opus 5 committed 2026-09-09 14:13:31 +02:00
1 parent 11c15837c6
commit e74517c2fb
1 file changed
+42 -13
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/FrenchWilson.cpp
+42 -13
View File
@@ -24,8 +24,23 @@ struct Posterior {
// so the exponentials never overflow/underflow. acentric: p(J) ~ exp(-J/S); centric:
// p(J) ~ exp(-J/2S)/sqrt(J).
// `logw` is caller-owned scratch of npts doubles (one per worker), so the integration allocates nothing.
// The quadrature grid is j = (i + 1/2) * dj, so sqrt(j) and log(j) separate into a per-reflection
// factor and a term that depends only on i: sqrt(j) = sqrt(dj) * sqrt(i + 1/2) and
// log(j) = log(dj) + log(i + 1/2). The i-dependent halves are the same for every reflection, so
// they are tabulated once instead of being recomputed npts times per intensity.
struct QuadratureGrid {
std::vector<double> sqrt_i, log_i;
explicit QuadratureGrid(int npts) : sqrt_i(npts), log_i(npts) {
for (int i = 0; i < npts; ++i) {
const double x = i + 0.5;
sqrt_i[i] = std::sqrt(x);
log_i[i] = std::log(x);
}
}
};
Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool centric, int npts,
std::vector<double> &logw) {
std::vector<double> &logw, const QuadratureGrid &grid) {
const double inv_2s2 = 1.0 / (2.0 * sigma * sigma);
// The posterior is the Gaussian likelihood tilted by the exponential prior, so it peaks at
// I - sigma^2/S and decays over whichever of sigma and S is TIGHTER. Ranging to I + 8 sigma
@@ -39,25 +54,34 @@ Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool
const double j_max = peak + 10.0 * width;
const double dj = j_max / npts;
const double log_dj = centric ? std::log(dj) : 0.0;
double max_logw = -std::numeric_limits<double>::infinity();
for (int i = 0; i < npts; ++i) {
const double j = (i + 0.5) * dj;
const double diff = I - j;
const double log_prior = centric ? (-j / (2.0 * sigma_wilson) - 0.5 * std::log(j))
: (-j / sigma_wilson);
const double log_prior = centric
? (-j / (2.0 * sigma_wilson) - 0.5 * (log_dj + grid.log_i[i]))
: (-j / sigma_wilson);
logw[i] = log_prior - diff * diff * inv_2s2;
max_logw = std::max(max_logw, logw[i]);
}
// exp(-37) is 8e-17: a point that far below the peak cannot change a normalised sum of
// doubles, and the posterior is sharply peaked, so most of the grid is skipped outright.
constexpr double LOG_NEGLIGIBLE = -37.0;
const double sqrt_dj = std::sqrt(dj);
double sum_w = 0, sum_wI = 0, sum_wF = 0;
for (int i = 0; i < npts; ++i) {
const double shifted = logw[i] - max_logw;
if (shifted < LOG_NEGLIGIBLE)
continue;
const double j = (i + 0.5) * dj;
const double w = std::exp(logw[i] - max_logw);
const double w = std::exp(shifted);
if (!std::isfinite(w))
continue;
sum_w += w;
sum_wI += w * j;
sum_wF += w * std::sqrt(j);
sum_wF += w * (sqrt_dj * grid.sqrt_i[i]);
}
if (sum_w <= 0.0) {
const double j = std::max(I, 0.0);
@@ -128,16 +152,17 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, const gemmi::Space
// French-Wilson |F| for one intensity of reflection r (its mean, or one Bijvoet hand); the shell
// Wilson prior, epsilon and centric flag are the reflection's, shared by all three.
// The Wilson prior and the centric flag belong to the reflection, not to the intensity, so they
// are looked up once and shared by its mean and both Bijvoet hands - three symmetry lookups
// became one.
auto fw_one = [&](const MergedReflection &r, float I, float sigma, float &F, float &sigF,
std::vector<double> &logw) {
std::vector<double> &logw, const QuadratureGrid &grid, double sigma_wilson,
bool centric) {
if (!std::isfinite(I) || !std::isfinite(sigma) || sigma <= 0.0f) { naive_one(I, sigma, F, sigF); return; }
// Strong reflections: the FW correction is negligible, <|F|> = sqrt(I).
if (I > opts.strong_cutoff * sigma) { naive_one(I, sigma, F, sigF); return; }
const auto s = shells.GetShell(r.d);
const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean);
const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}});
const Posterior post = integrate_posterior(I, sigma, sigma_wilson, centric,
opts.integration_points, logw);
opts.integration_points, logw, grid);
F = static_cast<float>(post.mean_F);
sigF = static_cast<float>(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F)));
};
@@ -146,13 +171,17 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, const gemmi::Space
const int n = static_cast<int>(merged.size());
const int nt = std::clamp(opts.num_threads, 1, n);
const int chunk = (n + nt - 1) / nt;
const QuadratureGrid grid(opts.integration_points);
auto do_chunk = [&](int lo, int hi) {
std::vector<double> logw(opts.integration_points);
for (int i = lo; i < hi; ++i) {
MergedReflection &r = merged[i];
fw_one(r, r.I, r.sigma, r.F, r.sigmaF, logw);
fw_one(r, r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus, logw);
fw_one(r, r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus, logw);
const auto s = shells.GetShell(r.d);
const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean);
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);
fw_one(r, r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus, logw, grid, sigma_wilson, centric);
}
};
if (nt == 1) {