diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index d89692ff7..60846728d 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1048,8 +1048,12 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s PlaceRiders(ctx, topo); if (cycle == 0 && fit != nullptr) ctx.Restore(*fit); - else + else { + // a trial's likelihood fit starts from nothing the slot's previous trial left + if (cycle == 0) + ctx.rows.bin_d.clear(); ctx.Refit(); + } if (adp_w > 0) ctx.MoveOverallB(); { diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp index 3cf0506bb..b3a014843 100644 --- a/image_analysis/structure_refinement/RefinementXray.cpp +++ b/image_analysis/structure_refinement/RefinementXray.cpp @@ -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 &fc, - const std::vector &idx, double d, double beta, bool with_sigma) { - double sum = 0, dfc = 0; - for (size_t i : idx) - sum += refine_math::AmplitudeNegLogLikelihood(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 fo, fc, var, log_var; - std::vector centric; - FixedVarianceRows(const RefinementReflections &refl, const std::vector &fc_all, - const std::vector &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 -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 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(bins), nthreads, [&](int b) { const std::vector 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 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++) { diff --git a/image_analysis/structure_refinement/RefinementXray.h b/image_analysis/structure_refinement/RefinementXray.h index 6cd615d70..8f35f9ce2 100644 --- a/image_analysis/structure_refinement/RefinementXray.h +++ b/image_analysis/structure_refinement/RefinementXray.h @@ -51,6 +51,8 @@ struct XrayRowModel { std::vector k_total; std::vector> fmask; std::vector d, var; + // D per likelihood bin, where EstimateLikelihoodParameters() found it - the next fit's start (empty: none) + std::vector bin_d; }; // The flat bulk-solvent model's mask (SolventMasker with the Refmac radii, PutMaskOnGrid) and its structure factors