From 210664f45453a02d30c4b6098bec2a6637459c2c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 17 Oct 2026 10:48:36 +0200 Subject: [PATCH] Model scale fit: no log ratio at |Fcalc| = 0, and no LM step once the model has underflowed Two ways the scale fit's normal matrix went singular and gemmi's jordan_solve threw "Trying to reverse singular matrix. Column 0 is zeroed.", ending the run: 1. The isotropic start (fit_isotropic_b_approximately) takes log(Fobs/|Fcalc|) over the strong reflections. At a phase-restricted reflection whose Fcalc and Fmask share a phase, some solvent pair of the k_sol/b_sol grid can cancel the model exactly in float, |Fcalc + k_sol exp(-b_sol s^2) Fmask| = 0. log(Fobs/0) = inf made that grid point's k_overall and B NaN, every sum of the Levenberg-Marquardt NaN, and the solve failed. This is what stopped three open-arm sets in the output pass after MR on the GPU (the scale fit after model validation's rigid body): one strong reflection, one grid point. The CPU path has the same exposure; which grid point hits an exact zero depends on the last bits of Fcalc, so the --no-gpu run of the same set went through. Such a reflection has no log ratio and is now left out of the start - in the GPU kernel, in FixedSolventFit::Start and, as a marked local change, in gemmi's own fit_isotropic_b_approximately, so the three stay equal and the held-solvent fits are covered too. 2. In FitModelScale's CPU Levenberg-Marquardt (ModelValidation_CCModelUnrelatedModelConfirmsNothing in a CUDA=OFF build): on 66 reflections of an unrelated model the fit walked to k_overall -424 and a B where k exp(-B s^2) |F| underflows at every reflection; the squares in J^T J were all 0 while J^T r kept a residue near 1e-210, which jordan_solve rejects. The fit now ends there, at its best so far, as it would at the lambda limit; that grid point then loses on R. The GPU fit takes the derivatives in float and their squares in double, so it cannot reach this state. Both changes only touch fits that threw before. Byte identity checked against an rc176 build on the standard and inhouse timing benchmarks and on 4nwv, 8tyy, 6oel (p.hkl, p_model.pdb, MTZ data, MR_* report keys). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- gemmi_gph/gemmi/scaling.hpp | 5 +++ .../structure_refinement/ModelScaleGPU.cu | 3 +- .../structure_refinement/ModelScaling.cpp | 7 ++++ tests/ModelScaleGPUTest.cpp | 20 +++++++++++ tests/ModelScalingTest.cpp | 33 +++++++++++++++++++ 5 files changed, 67 insertions(+), 1 deletion(-) diff --git a/gemmi_gph/gemmi/scaling.hpp b/gemmi_gph/gemmi/scaling.hpp index 630774c0b..519cec6e7 100644 --- a/gemmi_gph/gemmi/scaling.hpp +++ b/gemmi_gph/gemmi/scaling.hpp @@ -232,6 +232,11 @@ struct Scaling { if (p.fobs < 1 || p.fobs < p.sigma) // skip weak reflections continue; double fcalc = std::abs(get_fcalc(p)); + // Local change (upstream fits it): a reflection where the solvent exactly cancels the model, + // |Fcalc| = 0 in float, has no log ratio - its log(Fobs/0) = inf made k_overall and B NaN, + // and the solver that starts from them failed on a singular matrix. + if (fcalc == 0) + continue; double x = p.stol2; double y = std::log(static_cast(p.fobs / fcalc)); sx += x; diff --git a/image_analysis/structure_refinement/ModelScaleGPU.cu b/image_analysis/structure_refinement/ModelScaleGPU.cu index d77e4a6ed..f26a72744 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.cu +++ b/image_analysis/structure_refinement/ModelScaleGPU.cu @@ -71,7 +71,8 @@ __global__ void scale_sums(const int *__restrict__ hkl, const double *__restrict const double fo = fobs[i]; if (f.mode == MODE_ISOTROPIC) { - if (!strong[i]) + // fit_isotropic_b_approximately() (local change in gemmi_gph): no log ratio at |F| = 0 + if (!strong[i] || f_abs == 0) continue; const double x = stol2[i]; const double y = logf(static_cast(fo / f_abs)); diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp index d7b941c5a..6a099189b 100644 --- a/image_analysis/structure_refinement/ModelScaling.cpp +++ b/image_analysis/structure_refinement/ModelScaling.cpp @@ -65,6 +65,8 @@ struct FixedSolventFit { if (p.fobs < 1 || p.fobs < p.sigma) // skip weak reflections continue; double fcalc_abs = points[i].fcalc_abs; + if (fcalc_abs == 0) // local change in gemmi_gph: no log ratio + continue; double x = p.stol2; double y = std::log(static_cast(p.fobs / fcalc_abs)); sx += x; @@ -228,6 +230,11 @@ void LevMarFit(FixedSolventFit &target) { for (;;) { if (eval_limit > 0 && eval_count >= eval_limit) break; + // Not in gemmi: a fit that has walked to where k_overall exp(-B s^2) |F| underflows at every + // reflection has no step left to take - the squares in alpha are all 0 while beta keeps a + // residue near 1e-210, and jordan_solve() throws on that. It ends here, at its best so far. + if (alpha[0] == 0) + break; // prepare next parameters -> temp_beta std::copy(alpha, alpha + NA * NA, temp_alpha); diff --git a/tests/ModelScaleGPUTest.cpp b/tests/ModelScaleGPUTest.cpp index 86c167f83..67a52067f 100644 --- a/tests/ModelScaleGPUTest.cpp +++ b/tests/ModelScaleGPUTest.cpp @@ -382,4 +382,24 @@ TEST_CASE("ModelScaleGPUPool_MatchesFitModelScale", "[ModelValidation][gpu]") { } } +// ModelScaling fits past a reflection the solvent cancels exactly (ModelScalingTest.cpp), on the card. +TEST_CASE("ModelScaleGPUPool_FitsPastAReflectionTheSolventCancelsExactly", "[ModelValidation][gpu]") { + if (get_gpu_count() == 0) + SKIP("No GPU"); + gemmi::UnitCell cell(72, 44, 51, 90, 112, 90); + gemmi::Scaling scaling = MakeScaling(cell, *gemmi::find_spacegroup_by_name("C 1 2 1"), 2.5); + // exp(-b_sol s^2) = 1 at s = 0, so at the grid's k_sol = 0.10 the sum is -0.1f + 0.1f = 0 for every b_sol + auto &p = *std::find_if(scaling.points.begin(), scaling.points.end(), + [](const auto &q) { return q.fobs >= 1 && q.fobs >= q.sigma; }); + p.stol2 = 0; + p.fcmol = {-0.1f, 0.f}; + p.fmask = {1.f, 0.f}; + + ModelScaleGPUPool pool(0, scaling.points.size(), 1); + ModelScaleReport report; + REQUIRE_NOTHROW(report = pool.Fit(scaling, 80.0, 2)); + CHECK(std::isfinite(report.r_work_fit)); + CHECK(std::isfinite(scaling.k_overall)); +} + #endif diff --git a/tests/ModelScalingTest.cpp b/tests/ModelScalingTest.cpp index b5b5247f7..3898b033d 100644 --- a/tests/ModelScalingTest.cpp +++ b/tests/ModelScalingTest.cpp @@ -271,3 +271,36 @@ TEST_CASE("ModelScaling fits a grid point exactly as gemmi's Scaling does") { CHECK(ours.calculate_r_factor() == gemmis.calculate_r_factor()); } } + +// A reflection where the solvent term cancels the model exactly - |Fcalc + k_sol exp(-b_sol s^2) Fmask| = 0 +// in float - at some grid points. Its log(Fobs/|Fcalc|) was inf, which made the starting scale there NaN and +// the solver fail on a singular matrix. It has no log ratio, so the starting fit leaves it out. +TEST_CASE("ModelScaling fits past a reflection the solvent cancels exactly") { + const gemmi::UnitCell cell = Cell(); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21"); + REQUIRE(sg != nullptr); + gemmi::Scaling scaling(cell, sg); + scaling.use_solvent = true; + MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/20.0, + /*k_sol=*/0.35, /*b_sol=*/45.0, /*d_min=*/2.0, /*d_max=*/50.0); + // exp(-b_sol s^2) = 1 at s = 0, so at the grid's k_sol = 0.10 the sum is -0.1f + 0.1f = 0 for every b_sol + auto &p = scaling.points[0]; + p.stol2 = 0; + p.fcmol = {-0.1f, 0.f}; + p.fmask = {1.f, 0.f}; + REQUIRE(p.fobs >= 1); + + gemmi::Scaling gemmis = scaling; + ModelScaleReport report; + REQUIRE_NOTHROW(report = FitModelScale(scaling)); + CHECK(std::isfinite(report.r_work_fit)); + + gemmis.fix_k_sol = true; + gemmis.fix_b_sol = true; + gemmis.k_sol = 0.10; + gemmis.b_sol = 10.0; + REQUIRE(std::abs(gemmis.get_fcalc(gemmis.points[0])) == 0); + gemmis.fit_isotropic_b_approximately(); + CHECK(std::isfinite(gemmis.k_overall)); + REQUIRE_NOTHROW(gemmis.fit_parameters()); +}