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

Refinement: the likelihood's D and beta searches by Newton's method; D on the working set from the last fit's

Browse Source 
Deterministic, rounding changes. The per-bin ML fit took 4 rounds x (40 + 40) golden-section steps on the
free set and a 50-step golden section for D on the working set - each step a Bessel I0/I1 per row, ~4.6M per
Refit at 92k reflections, 27-37% of a trial's wall on the high-resolution sets. Each one-dimensional search is
now a safeguarded Newton iteration on the analytic first and second derivatives (bracketed on the sign of the
derivative, bisection where a step leaves the bracket or the curvature is not positive), within the same
brackets; the free-set alternation keeps its 4 rounds from the moments, so beta is the one it was (a converged
joint fit was tried: the 4 rounds are not converged in the low-resolution bins, and converging them moves
6iu9's FOM by -0.003, so it was not kept). D on the working set starts from the bin's D of the context's
previous Refit (XrayRowModel::bin_d), restored with a trial's fit and cleared at a trial's first Refit, so a
trial depends on nothing its slot ran before (side by side = one at a time, md5).

14 parity sets, score only: R-work/R-free identical to 4 digits, mean |FOM - REFMAC| 0.01186 -> 0.01186,
D within 1e-5 and var within 1e-4 (relative) per reflection. refit.ml (TIMING sum, -N 8): 5reo 1.47 -> 0.11 s,
5f6m 3.88 -> 0.26, 9yl4 1.37 -> 0.11, 5f6m MD start 4.97 -> 0.46. Wall: 5reo 2.03 -> 1.75, 5f6m 3.45 -> 2.60,
9yl4 7.99 -> 7.36, 7d98 4.04 -> 3.78, 5f6m MD start 6.10 -> 4.82 s. 5f6m MD starts 01-04 median R-free
0.1693 -> 0.1692.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit a91511945e0ccd7145a93877c0c231e87f86fb44)
This commit is contained in:
leonarski_f committed 2026-10-26 09:16:33 +01:00
1 parent a8d3a8dbc2
commit 2c1df0ebfa
3 files changed
+83 -74
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/structure_refinement/RefinementXray.cpp
+76 -73
View File
@@ -49,76 +49,56 @@ double Variance(const RefinementReflections &refl, size_t i, double beta, bool w
return refl.epsilon[i] * beta + (refl.centric[i] ? s2 : 2.0 * s2);
}
// -ln L summed over rows `idx` at D and beta.
double BinNegLogLikelihood(const RefinementReflections &refl, const std::vector<double> &fc,
const std::vector<size_t> &idx, double d, double beta, bool with_sigma) {
double sum = 0, dfc = 0;
for (size_t i : idx)
sum += refine_math::AmplitudeNegLogLikelihood<double>(refl.fo[i], fc[i], d, Variance(refl, i, beta, with_sigma),
refl.centric[i], dfc);
return sum;
}
// Rows `idx` at a fixed beta, gathered for searches over D alone: the variance and its logarithm taken once rather
// than at every evaluation. FixedVarianceNegLogLikelihood() at D is BinNegLogLikelihood() at D and that beta, to
// the bit: the same expressions in the same order, only log(var) hoisted.
struct FixedVarianceRows {
std::vector<double> fo, fc, var, log_var;
std::vector<char> centric;
FixedVarianceRows(const RefinementReflections &refl, const std::vector<double> &fc_all,
const std::vector<size_t> &idx, double beta, bool with_sigma) {
for (size_t i : idx) {
fo.push_back(refl.fo[i]);
fc.push_back(fc_all[i]);
var.push_back(Variance(refl, i, beta, with_sigma));
log_var.push_back(std::log(var.back()));
centric.push_back(refl.centric[i]);
}
}
// The first and second derivatives of one row's -ln L (refine_math::AmplitudeNegLogLikelihood) in D and in the
// variance, at D and var.
struct RowDerivatives {
double d, dd, v, vv;
};
// refine_math::AmplitudeNegLogLikelihood() summed over the rows at D, with log(var) given.
double FixedVarianceNegLogLikelihood(const FixedVarianceRows &rows, double d) {
double sum = 0;
for (size_t k = 0; k < rows.fo.size(); k++) {
const double fo = rows.fo[k], fc = rows.fc[k], var = rows.var[k];
if (!rows.centric[k]) {
const double x = 2.0 * d * fo * fc / var;
double log_i0, ratio;
refine_math::BesselI0I1(x, log_i0, ratio);
sum += rows.log_var[k] + (fo * fo + d * d * fc * fc) / var - log_i0;
} else {
const double x = d * fo * fc / var;
const double log_cosh = x + std::log(1.0 + std::exp(-2.0 * x)) - 0.6931471805599453;
sum += 0.5 * rows.log_var[k] + (fo * fo + d * d * fc * fc) / (2.0 * var) - log_cosh;
}
RowDerivatives LikelihoodDerivatives(double fo, double fc, double d, double var, bool centric) {
const double q = fo * fo + d * d * fc * fc;
if (!centric) {
// x = 2 D Fo Fc / var, m = I1(x) / I0(x) = d ln I0 / dx, m' = 1 - m / x - m^2
const double a = 2.0 * fo * fc / var, x = a * d;
double log_i0, m;
refine_math::BesselI0I1(x, log_i0, m);
const double m1 = x > 0 ? 1.0 - m / x - m * m : 0.5;
return {2.0 * d * fc * fc / var - a * m,
2.0 * fc * fc / var - a * a * m1,
1.0 / var - q / (var * var) + m * x / var,
(-1.0 + 2.0 * q / var - x * (m1 * x + 2.0 * m)) / (var * var)};
}
return sum;
// x = D Fo Fc / var, t = tanh x = d ln cosh / dx
const double a = fo * fc / var, x = a * d;
const double t = std::tanh(x), t1 = 1.0 - t * t;
return {d * fc * fc / var - a * t,
fc * fc / var - a * a * t1,
0.5 / var - q / (2.0 * var * var) + t * x / var,
(-0.5 + q / var - x * (t1 * x + 2.0 * t)) / (var * var)};
}
// The minimum of f over [lo, hi] by golden-section search.
// The minimum of a function over [lo, hi] from x, given its first and second derivative (derivatives(x, g, h)):
// Newton's method on the first derivative, kept inside the interval where that changes sign and bisecting it
// wherever a step would leave it or the curvature is not positive. A minimum on a bound is the bound.
template <typename F>
double GoldenMinimum(F f, double lo, double hi, int iterations) {
const double g = 0.5 * (std::sqrt(5.0) - 1.0);
double a = lo, b = hi;
double c = b - g * (b - a), d = a + g * (b - a);
double fc = f(c), fd = f(d);
for (int i = 0; i < iterations; i++) {
if (fc < fd) {
b = d;
d = c;
fd = fc;
c = b - g * (b - a);
fc = f(c);
} else {
a = c;
c = d;
fc = fd;
d = a + g * (b - a);
fd = f(d);
}
double NewtonMinimum(F derivatives, double x, double lo, double hi) {
x = std::clamp(x, lo, hi);
for (int iteration = 0; iteration < 200; iteration++) {
double g = 0, h = 0;
derivatives(x, g, h);
if (g < 0)
lo = x;
else
hi = x;
double next = h > 0 ? x - g / h : lo - 1.0;
if (!(next > lo && next < hi))
next = 0.5 * (lo + hi);
const bool converged = std::fabs(next - x) < 1e-10 || hi - lo < 1e-10;
x = next;
if (converged)
break;
}
return 0.5 * (a + b);
return x;
}
class RefinementEngineCPU : public RefinementEngine {
@@ -333,6 +313,7 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
}
std::vector<double> centre(bins), bin_d(bins), bin_beta(bins);
const bool warm = rows.bin_d.size() == bins;
// each bin is fitted on its own
ParallelFor(static_cast<int>(bins), nthreads, [&](int b) {
const std::vector<size_t> idx(free.begin() + b * free.size() / bins, free.begin() + (b + 1) * free.size() / bins);
@@ -349,25 +330,47 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
// beta: the spread of the free amplitudes about D Fc, measurement error included, by maximum likelihood
// (alternating one-dimensional minimisations over ln D and ln beta)
for (int round = 0; round < 4; round++) {
const FixedVarianceRows at_beta(refl, fc, idx, beta, false);
const double ld = GoldenMinimum([&](double x) { return FixedVarianceNegLogLikelihood(at_beta, std::exp(x)); },
std::log(d) - 2.0, std::log(std::min(d * 7.4, 3.0)), 40);
const double ld = NewtonMinimum([&](double x, double &g, double &h) {
const double dx = std::exp(x);
for (size_t i : idx) {
const RowDerivatives r = LikelihoodDerivatives(refl.fo[i], fc[i], dx, refl.epsilon[i] * beta,
refl.centric[i]);
g += dx * r.d;
h += dx * dx * r.dd + dx * r.d;
}
}, std::log(d), std::log(d) - 2.0, std::log(std::min(d * 7.4, 3.0)));
d = std::exp(ld);
const double lb = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, idx, d, std::exp(x), false); },
std::log(beta) - 3.0, std::log(beta) + 3.0, 40);
const double lb = NewtonMinimum([&](double x, double &g, double &h) {
for (size_t i : idx) {
const double var = refl.epsilon[i] * std::exp(x);
const RowDerivatives r = LikelihoodDerivatives(refl.fo[i], fc[i], d, var, refl.centric[i]);
g += var * r.v;
h += var * var * r.vv + var * r.v;
}
}, std::log(beta), std::log(beta) - 3.0, std::log(beta) + 3.0);
beta = std::exp(lb);
}
// D: by maximum likelihood on the working reflections, with the measurement error added to that beta.
// Where the error is a large part of the spread, the variance counts it twice and D falls towards 0 -
// the high-resolution D and figure of merit REFMAC 5.8 reports.
// D: by maximum likelihood on the working reflections, with the measurement error added to that beta, from
// where the last fit found it. Where the error is a large part of the spread, the variance counts it twice
// and D falls towards 0 - the high-resolution D and figure of merit REFMAC 5.8 reports.
// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 (D from the working set, the variance from the free)
if (!work_in_bin[b].empty()) {
const FixedVarianceRows working(refl, fc, work_in_bin[b], beta, true);
d = GoldenMinimum([&](double x) { return FixedVarianceNegLogLikelihood(working, x); }, 0.0, 3.0, 50);
std::vector<double> var;
for (size_t i : work_in_bin[b])
var.push_back(Variance(refl, i, beta, true));
d = NewtonMinimum([&](double x, double &g, double &h) {
for (size_t k = 0; k < var.size(); k++) {
const size_t i = work_in_bin[b][k];
const RowDerivatives r = LikelihoodDerivatives(refl.fo[i], fc[i], x, var[k], refl.centric[i]);
g += r.d;
h += r.dd;
}
}, warm ? rows.bin_d[b] : d, 0.0, 3.0);
}
bin_d[b] = d;
bin_beta[b] = beta;
});
rows.bin_d = bin_d;
rows.d.resize(n);
rows.var.resize(n);
for (size_t i = 0; i < n; i++) {