Refinement: a Refit's solvent grid starts around where the context's last one ended

Logic change (deterministic), kept on evidence. Every Refit searched the whole k_sol/b_sol box - an 88-point
coarse pass of plain least squares, then the robust fine pass around its winner. ModelScaleBox (and
ModelScaleGPUPool::Fit / ModelScaleGPU::FitSolvent) take an optional start; with one, the coarse pass is
skipped and the fine pass is laid around it. Refinement passes the grid winner of the context's previous
Refit (before MoveOverallB moves the overall B into b_sol); it travels with a trial's fit, and a trial's first
Refit without one (Context::ForgetStarts, which also clears the likelihood's D start) searches the whole box,
so a trial still depends on nothing its slot ran before (side by side = one at a time, md5).

Evidence (GPU, -N 8): full refinement of the 14 parity sets from their placed models - 12 refined CIFs
identical (md5), 9ig7 and 6iu9 differ, final R-work/R-free equal to 4 digits but 6iu9 R-work 0.2911 -> 0.2910;
of the 5 timing cases 4 identical, 7d98 (cif start) R-free 0.3096 -> 0.3095; 5f6m MD starts 01-04 median R-free
0.1692 -> 0.1692; 6o2h on the CPU engine identical (md5). The score-only parity table is not touched (one Refit,
no start). refit.scale (TIMING sum): 5f6m 1.52 -> 1.10 s. Wall: 5reo 1.68 -> 1.53, 5f6m 2.43 -> 2.24, 9yl4 7.19 ->
7.02, 7d98 3.70 -> 3.65, 5f6m MD start 4.57 -> 4.22 s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit ff285fb525ecf4b480698d4aafe212a34301a8fe)

image_analysis/structure_refinement/ModelScaleGPU.cu

@@ -597,7 +597,8 @@ ModelScaleParams ModelScaleGPU::Fit(const float2 *d_fcmol, const float2 *d_fmask
 // FitModelScale's grid (ModelScaling.cpp, default ModelScaleBox with b_sol up to b_hi): a coarse pass
 // over the box and a finer one around its winner, every grid point its own fit, the winner the lowest
 // finite R in grid order. The fits of a pass run together, one launch per Levenberg-Marquardt step.
-ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d_fmask, double b_hi, bool robust) {
+ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d_fmask, double b_hi, bool robust,
+                                          double k_start, double b_start) {
     ModelSolventFit out;
     if (n_ < 20)
         return out;
@@ -631,12 +632,16 @@ ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d
         return ModelScaleFitState{1.0, {0, 0, 0, 0, 0, 0}, ks, bs, 0, 0};
     };
 
-    std::vector<ModelScaleFitState> coarse;
-    for (double ks = k_lo; ks <= k_hi + 1e-9; ks += 0.05)
-        for (double bs = b_lo; bs <= b_hi + 1e-9; bs += b_step)
-            coarse.push_back(point(ks, bs));
-    try_points(coarse, false);
-    const double k0 = out.k_sol, b0 = out.b_sol;
+    double k0 = k_start, b0 = b_start;
+    if (k0 < 0) {
+        std::vector<ModelScaleFitState> coarse;
+        for (double ks = k_lo; ks <= k_hi + 1e-9; ks += 0.05)
+            for (double bs = b_lo; bs <= b_hi + 1e-9; bs += b_step)
+                coarse.push_back(point(ks, bs));
+        try_points(coarse, false);
+        k0 = out.k_sol;
+        b0 = out.b_sol;
+    }
     const double k_hi2 = std::min(k_hi, k0 + 0.05);
     const double b_hi2 = std::min(b_hi, b0 + b_step);
     std::vector<ModelScaleFitState> fine;

image_analysis/structure_refinement/ModelScaleGPU.h

@@ -63,8 +63,9 @@ public:
     // (five or fewer reflections for the isotropic fit).
     // robust: the fine pass by robust least squares, as FitModelScale with ModelScaleBox::robust - the
     // reweighting on the device, the median of each fit's first residuals included.
+    // k_start, b_start: ModelScaleBox's - the fine pass around them, no coarse pass (negative: as FitModelScale).
     ModelSolventFit FitSolvent(const float2 *d_fcmol, const float2 *d_fmask, double b_hi = 80.0,
-                               bool robust = false);
+                               bool robust = false, double k_start = -1, double b_start = -1);
     // What a robust FitSolvent takes on the card beyond DeviceBytes(), reserved at its first call.
     static size_t RobustBytes(size_t max_points);
 

image_analysis/structure_refinement/ModelScaling.cpp

@@ -427,12 +427,16 @@ ModelScaleReport FitModelScale(gemmi::Scaling<float> &scaling, ModelScaleBox box
     // 10 A^2, or in eight values where the box is wider than the default, so a raised ceiling costs
     // no more grid points.
     const double b_step = std::max(10.0, (box.b_hi - box.b_lo) / 7.0);
-    std::vector<PointFit> coarse;
-    for (double ks = box.k_lo; ks <= box.k_hi + 1e-9; ks += 0.05)
-        for (double bs = box.b_lo; bs <= box.b_hi + 1e-9; bs += b_step)
-            coarse.push_back(PointFit{ks, bs});
-    try_points(coarse);
-    const double k0 = best_k_sol, b0 = best_b_sol;
+    double k0 = box.k_start, b0 = box.b_start;
+    if (k0 < 0) {
+        std::vector<PointFit> coarse;
+        for (double ks = box.k_lo; ks <= box.k_hi + 1e-9; ks += 0.05)
+            for (double bs = box.b_lo; bs <= box.b_hi + 1e-9; bs += b_step)
+                coarse.push_back(PointFit{ks, bs});
+        try_points(coarse);
+        k0 = best_k_sol;
+        b0 = best_b_sol;
+    }
     const double k_hi2 = std::min(box.k_hi, k0 + 0.05);
     const double b_hi2 = std::min(box.b_hi, b0 + b_step);
     std::vector<PointFit> fine;
@@ -511,7 +515,8 @@ ModelScaleGPUPool::~ModelScaleGPUPool() {
     set_gpu(current);
 }
 
-ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling<float> &scaling, double b_hi, size_t nthreads, bool robust) {
+ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling<float> &scaling, double b_hi, size_t nthreads, bool robust,
+                                        double k_start, double b_start) {
     const size_t n = scaling.points.size();
     std::vector<std::array<int, 3>> hkl(n);
     std::vector<double> stol2(n);
@@ -554,11 +559,11 @@ ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling<float> &scaling, double b
         if (cudaMemcpyAsync(e->fcalc, fc.data(), n * sizeof(float2), cudaMemcpyHostToDevice, e->stream) != cudaSuccess
             || cudaMemcpyAsync(e->fmask, fm.data(), n * sizeof(float2), cudaMemcpyHostToDevice, e->stream) != cudaSuccess)
             throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "ModelScaleGPUPool: the upload failed");
-        fit = e->scale.FitSolvent(e->fcalc, e->fmask, b_hi, robust);
+        fit = e->scale.FitSolvent(e->fcalc, e->fmask, b_hi, robust, k_start, b_start);
     } catch (const ModelScaleGPUTooFewReflections &) {
         set_gpu(current);
         release();
-        return FitModelScale(scaling, {.b_hi = b_hi, .robust = robust}, nthreads);
+        return FitModelScale(scaling, {.b_hi = b_hi, .robust = robust, .k_start = k_start, .b_start = b_start}, nthreads);
     } catch (...) {
         set_gpu(current);
         release();

image_analysis/structure_refinement/ModelScaling.h

@@ -34,6 +34,9 @@ struct ModelScaleBox {
     // The overall scale and anisotropic B of the points of the fine pass by robust least squares
     // (ModelScaling.cpp, RobustReweight) rather than plain least squares on |F|: what refinement scores a model with.
     bool robust = false;
+    // Where a fit of the model as it was a moment ago found the solvent: the fine pass around (k_start, b_start)
+    // without the coarse pass. Negative: the coarse pass finds where to look.
+    double k_start = -1, b_start = -1;
 };
 
 // The b_sol ceiling for `model`: 80 A^2, or twice the model's mean atomic B where that is higher.
@@ -73,7 +76,9 @@ public:
     // Safe to call from several threads. Throws JFJochException on a CUDA failure. The robust fit's reweighting
     // runs on the device; it is FitModelScale's robust fit to rounding (float residuals and weights, the
     // anisotropic factor of the steps in float), deterministic run to run.
-    ModelScaleReport Fit(gemmi::Scaling<float> &scaling, double b_hi, size_t nthreads, bool robust = false);
+    // k_start, b_start: ModelScaleBox's.
+    ModelScaleReport Fit(gemmi::Scaling<float> &scaling, double b_hi, size_t nthreads, bool robust = false,
+                         double k_start = -1, double b_start = -1);
     // What the pool holds on the card.
     size_t Bytes() const { return DeviceBytes(max_points_, engines_.size()); }
 

image_analysis/structure_refinement/Refinement.cpp

@@ -126,6 +126,9 @@ struct Context {
     XrayRowModel rows;
     std::vector<std::complex<float>> fcalc, fmask;
     double k_sol = 0, b_sol = 0;
+    // where the last Refit()'s solvent grid ended (before MoveOverallB()): the next one's fine pass is around it
+    // (negative: the next one searches the whole box)
+    double grid_k_sol = -1, grid_b_sol = -1;
     double k_overall = 0;
     gemmi::SMat33<double> b_star{0, 0, 0, 0, 0, 0};
     double b_floor = B_FLOOR;   // the cycle's lowest allowed B
@@ -189,8 +192,9 @@ struct Context {
         double k_sol = 0, b_sol = 0, k_overall = 0;
         gemmi::SMat33<double> b_star{0, 0, 0, 0, 0, 0};
         XrayRowModel rows;
+        double grid_k_sol = -1, grid_b_sol = -1;
     };
-    Fit Fitted() const { return {b_floor, k_sol, b_sol, k_overall, b_star, rows}; }
+    Fit Fitted() const { return {b_floor, k_sol, b_sol, k_overall, b_star, rows, grid_k_sol, grid_b_sol}; }
     // Refit() as it was where `fit` was taken, for the atoms as they were there.
     void Restore(const Fit &fit) {
         b_floor = fit.b_floor;
@@ -200,8 +204,15 @@ struct Context {
         k_overall = fit.k_overall;
         b_star = fit.b_star;
         rows = fit.rows;
+        grid_k_sol = fit.grid_k_sol;
+        grid_b_sol = fit.grid_b_sol;
         engine->SetRowModel(rows);
     }
+    // The next Refit() from nothing an earlier one found: the whole solvent box, the likelihood from the moments.
+    void ForgetStarts() {
+        grid_k_sol = grid_b_sol = -1;
+        rows.bin_d.clear();
+    }
 
     // A new macro-cycle: the blur for the B the cycle may reach, the bulk-solvent mask of the model as it is,
     // the scale and the solvent fitted to the working set, and the likelihood's parameters from the free set.
@@ -244,13 +255,14 @@ struct Context {
         // the scale of the score: by robust least squares (ModelScaleBox::robust)
 #ifdef JFJOCH_USE_CUDA
         if (scale_gpu)
-            scale_gpu->Fit(scaling, SolventBMax(Model()), nthreads, true);
+            scale_gpu->Fit(scaling, SolventBMax(Model()), nthreads, true, grid_k_sol, grid_b_sol);
         else
 #endif
-            FitModelScale(scaling, {.b_hi = SolventBMax(Model()), .robust = true}, nthreads);
+            FitModelScale(scaling, {.b_hi = SolventBMax(Model()), .robust = true, .k_start = grid_k_sol,
+                                    .b_start = grid_b_sol}, nthreads);
         s3.reset();
-        k_sol = scaling.k_sol;
-        b_sol = scaling.b_sol;
+        k_sol = grid_k_sol = scaling.k_sol;
+        b_sol = grid_b_sol = scaling.b_sol;
         k_overall = scaling.k_overall;
         b_star = scaling.b_star;
         rows.k_total.resize(refl.Size());
@@ -1053,9 +1065,9 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector<int> &refinable, const s
         if (cycle == 0 && fit != nullptr)
             ctx.Restore(*fit);
         else {
-            // a trial's likelihood fit starts from nothing the slot's previous trial left
+            // a trial starts from nothing the slot's previous trial left
             if (cycle == 0)
-                ctx.rows.bin_d.clear();
+                ctx.ForgetStarts();
             ctx.Refit();
         }
         if (adp_w > 0)