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

ModelScaling: the LM matrices after an accepted step reuse the step's exp(-B/4 h.h)

Browse Source 
LevMarFit takes compute_wssr() in a copy (ComputeWSSR) that keeps each
reflection's exp(-0.25 * b_star.r_u_r(h)); when the step is accepted,
ComputeLMMatrices is evaluated at exactly those parameters and takes the same
values from there instead of calling exp again. Same expressions, same
parameters, same numbers.

Bench against the base FitModelScale (single thread, 48k reflections, P2_1):
1.65 -> 1.10 s for this and the previous commit together; P1, P2_12_12_1,
P6_122, I23 and the dependent n<=5 path identical to the bit and 30-40% faster.
Validation outputs (p.hkl md5, p.mtz / p_maps.mtz data, maps, p_model.pdb, log)
identical to the base: CPU build --mode scale on 5reo, myob_x10sa, 8sa8, 8xtg,
5ebi at -N 1, default and -N 8; full CPU runs of 9ea5, 6iu5, 7dkp; the GPU build
with no visible device on the five --mode scale sets.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_f
2026-10-10 11:22:26 +02:00
co-authored by Claude Opus 5.5
parent 76f9192fb4
commit 7a6b194e09
1 changed files with 30 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/structure_refinement/ModelScaling.cpp
+30 -11
View File
@@ -105,11 +105,11 @@ struct FixedSolventFit {
double compute_value(const Point &p) const {
return p.fcalc_abs * (float) get_overall_scale_factor(p.hkl);
}
// For NA parameters: the B components are NA - 1 constraint rows.
// For NA parameters: the B components are NA - 1 constraint rows. `kaniso` is
// exp(-0.25 * b_star.r_u_r(h)) at the current parameters, as compute_value() takes it too.
template <int NA>
double compute_value_and_derivatives(const Point &p, double *dy_da) const {
double compute_value_and_derivatives(const Point &p, double kaniso, double *dy_da) const {
gemmi::Vec3 h(p.hkl);
double kaniso = std::exp(-0.25 * b_star.r_u_r(h));
double fcalc_abs = p.fcalc_abs;
int n = 1;
double fe = fcalc_abs * kaniso;
@@ -149,13 +149,18 @@ double RFactor(const FixedSolventFit &fit) {
// adding +0 would change, since it starts at +0 and a sum is only -0 when both its terms are. And
// where gemmi fills the lower half of alpha, the whole square is summed here, a rectangle the compiler
// can vectorise; the lower half is gemmi's, and the upper is replaced by it at the end, as in gemmi.
//
// `kaniso`: each reflection's exp(-0.25 * b_star.r_u_r(h)) at the current parameters where ComputeWSSR()
// has just taken them there (nullptr: computed here) - the same value, which need not be taken twice.
template <int NA>
double ComputeLMMatrices(const FixedSolventFit &fit, double *alpha_out, double *beta_out) {
double ComputeLMMatrices(const FixedSolventFit &fit, const double *kaniso, double *alpha_out, double *beta_out) {
long double wssr = 0; // long double here notably increases the accuracy
double alpha[NA * NA] = {}, beta[NA] = {};
double dy_da[NA];
for (const auto &p : fit.points) {
double y = fit.compute_value_and_derivatives<NA>(p, dy_da);
for (size_t i = 0; i < fit.points.size(); ++i) {
const auto &p = fit.points[i];
const double k = kaniso ? kaniso[i] : std::exp(-0.25 * fit.b_star.r_u_r(gemmi::Vec3(p.hkl)));
double y = fit.compute_value_and_derivatives<NA>(p, k, dy_da);
double dy_sig = p.get_y() - y;
for (int j = 0; j != NA; ++j) {
const bool used = dy_da[j] != 0;
@@ -177,8 +182,21 @@ double ComputeLMMatrices(const FixedSolventFit &fit, double *alpha_out, double *
return (double) wssr;
}
// gemmi's LevMar::fit() (levmar.hpp) with its default settings, on ComputeLMMatrices<NA>() - the same
// iterations to the same answer, bit for bit.
// gemmi's compute_wssr() (levmar.hpp) on a FixedSolventFit, keeping each reflection's
// exp(-0.25 * b_star.r_u_r(h)) in `kaniso` for ComputeLMMatrices() at the same parameters.
double ComputeWSSR(const FixedSolventFit &fit, std::vector<double> &kaniso) {
long double wssr = 0; // long double here notably increases the accuracy
for (size_t i = 0; i < fit.points.size(); ++i) {
const auto &p = fit.points[i];
kaniso[i] = std::exp(-0.25 * fit.b_star.r_u_r(p.hkl));
// FixedSolventFit::compute_value(), weight 1
wssr += gemmi::sq(p.get_y() - p.fcalc_abs * (float) (fit.k_overall * kaniso[i]));
}
return (double) wssr;
}
// gemmi's LevMar::fit() (levmar.hpp) with its default settings, on ComputeLMMatrices<NA>() and
// ComputeWSSR() - the same iterations to the same answer, bit for bit.
template <int NA>
void LevMarFit(FixedSolventFit &target) {
const int eval_limit = 100;
@@ -194,7 +212,8 @@ void LevMarFit(FixedSolventFit &target) {
double alpha[NA * NA], beta[NA], temp_alpha[NA * NA];
std::vector<double> temp_beta(NA);
const double initial_wssr = ComputeLMMatrices<NA>(target, alpha, beta);
std::vector<double> kaniso(target.points.size());
const double initial_wssr = ComputeLMMatrices<NA>(target, nullptr, alpha, beta);
double wssr = initial_wssr;
int small_change_counter = 0;
@@ -219,7 +238,7 @@ void LevMarFit(FixedSolventFit &target) {
temp_beta[i] += best_a[i];
target.set_parameters(temp_beta);
double new_wssr = gemmi::compute_wssr(target);
double new_wssr = ComputeWSSR(target, kaniso);
++eval_count;
if (new_wssr < wssr) {
double rel_change = (wssr - new_wssr) / wssr;
@@ -235,7 +254,7 @@ void LevMarFit(FixedSolventFit &target) {
} else {
small_change_counter = 0;
}
ComputeLMMatrices<NA>(target, alpha, beta);
ComputeLMMatrices<NA>(target, kaniso.data(), alpha, beta); // at the parameters just tried
++eval_count;
lambda *= lambda_down_factor;
} else { // worse fitting