From 1f8b5b1986aa3cd478742f2cfe7c3a9ee0065ea5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 26 Oct 2026 08:08:32 +0100 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit ff285fb525ecf4b480698d4aafe212a34301a8fe) --- .../structure_refinement/ModelScaleGPU.cu | 19 +++++++++----- .../structure_refinement/ModelScaleGPU.h | 3 ++- .../structure_refinement/ModelScaling.cpp | 23 +++++++++------- .../structure_refinement/ModelScaling.h | 7 ++++- .../structure_refinement/Refinement.cpp | 26 ++++++++++++++----- 5 files changed, 53 insertions(+), 25 deletions(-) diff --git a/image_analysis/structure_refinement/ModelScaleGPU.cu b/image_analysis/structure_refinement/ModelScaleGPU.cu index 39267021d..d16649fcf 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.cu +++ b/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 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 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 fine; diff --git a/image_analysis/structure_refinement/ModelScaleGPU.h b/image_analysis/structure_refinement/ModelScaleGPU.h index 4d3fec578..c68cf2b3d 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.h +++ b/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); diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp index e62931b62..67f04a930 100644 --- a/image_analysis/structure_refinement/ModelScaling.cpp +++ b/image_analysis/structure_refinement/ModelScaling.cpp @@ -427,12 +427,16 @@ ModelScaleReport FitModelScale(gemmi::Scaling &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 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 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 fine; @@ -511,7 +515,8 @@ ModelScaleGPUPool::~ModelScaleGPUPool() { set_gpu(current); } -ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling &scaling, double b_hi, size_t nthreads, bool robust) { +ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling &scaling, double b_hi, size_t nthreads, bool robust, + double k_start, double b_start) { const size_t n = scaling.points.size(); std::vector> hkl(n); std::vector stol2(n); @@ -554,11 +559,11 @@ ModelScaleReport ModelScaleGPUPool::Fit(gemmi::Scaling &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(); diff --git a/image_analysis/structure_refinement/ModelScaling.h b/image_analysis/structure_refinement/ModelScaling.h index 9ec565a24..8630d62f9 100644 --- a/image_analysis/structure_refinement/ModelScaling.h +++ b/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 &scaling, double b_hi, size_t nthreads, bool robust = false); + // k_start, b_start: ModelScaleBox's. + ModelScaleReport Fit(gemmi::Scaling &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()); } diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 760b19f60..5d532a825 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -126,6 +126,9 @@ struct Context { XrayRowModel rows; std::vector> 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 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 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 &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)