diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp index 8337c8c68..1ce87e07d 100644 --- a/image_analysis/structure_refinement/RefinementXray.cpp +++ b/image_analysis/structure_refinement/RefinementXray.cpp @@ -59,6 +59,43 @@ double BinNegLogLikelihood(const RefinementReflections &refl, const std::vector< 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]); + } + } +}; + +// 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; + } + } + return sum; +} + // The minimum of f over [lo, hi] by golden-section search. template double GoldenMinimum(F f, double lo, double hi, int iterations) { @@ -304,7 +341,8 @@ 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 double ld = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, idx, std::exp(x), beta, false); }, + 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); d = std::exp(ld); const double lb = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, idx, d, std::exp(x), false); }, @@ -315,9 +353,10 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std:: // 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()) - d = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, work_in_bin[b], x, beta, true); }, - 0.0, 3.0, 50); + 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); + } bin_d[b] = d; bin_beta[b] = beta; });