From a8d3a8dbc23dc66362e24f611a063954c7e80982 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Mon, 26 Oct 2026 07:33:16 +0100
Subject: [PATCH] Refinement: the mask grid sized without allocating it; the
 L-BFGS pairs keep their s.y

BulkSolventMask::MaskGrid, documented as "sizes and spacing set, values not allocated", allocated and zeroed the
whole fine mask grid (up to 2^26 floats) - in BulkSolventMask's constructor and in every GPU engine's, i.e. once
more for each weight-trial slot, only for the grid's size. It now takes good_grid_size() as set_size_from_spacing()
does, and the CPU engine's Fmask() allocates its copy. rugnux_refine_test, RTX 5080, -N 8: 9yl4 8.6 -> 7.8 s,
7d98 4.85 -> 4.0 s (engine 0.66 -> 0.16 s, slots 0.60 -> 0.23 s on 9yl4); the high-resolution sets, whose mask
grid is small, are unchanged.

RefinementMinimiser keeps s.y of every pair it stores and Direction() takes rho and gamma from it instead of two
more dot products per pair per call (the same Dot() of the same vectors): Direction's CPU time -22% on 9yl4; no
wall change with the trials side by side, where other slots fill the card meanwhile.

Bitwise: refined CIF identical on 5reo, 5f6m, 9yl4, 7d98, a 5f6m MD start (GPU) and 6o2h (CPU engine).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 6f5caf5e1b1183a752acc5cca78bfb7ef002b12a)
---
 .../structure_refinement/RefinementMinimiser.h         |  8 ++++++--
 image_analysis/structure_refinement/RefinementXray.cpp | 10 +++++++++-
 2 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/image_analysis/structure_refinement/RefinementMinimiser.h b/image_analysis/structure_refinement/RefinementMinimiser.h
index 42bb3afe1..82ac9e118 100644
--- a/image_analysis/structure_refinement/RefinementMinimiser.h
+++ b/image_analysis/structure_refinement/RefinementMinimiser.h
@@ -27,6 +27,7 @@ public:
     void Reset() {
         s_.clear();
         y_.clear();
+        sy_.clear();
     }
     // The estimated diagonal of the Hessian; every element > 0.
     void SetDiagonal(std::vector<double> h) { h_ = std::move(h); }
@@ -94,9 +95,11 @@ public:
             if (sy > 1e-10 * std::sqrt(Dot(s, s) * Dot(y, y))) {
                 s_.push_back(std::move(s));
                 y_.push_back(std::move(y));
+                sy_.push_back(sy);
                 if (s_.size() > memory_) {
                     s_.pop_front();
                     y_.pop_front();
+                    sy_.pop_front();
                 }
             }
             const double decrease = f - f_new;
@@ -124,7 +127,7 @@ private:
         const size_t n = g.size(), m = s_.size();
         std::vector<double> q = g, a(m), rho(m);
         for (size_t i = m; i-- > 0;) {
-            rho[i] = 1.0 / Dot(s_[i], y_[i]);
+            rho[i] = 1.0 / sy_[i];
             a[i] = rho[i] * Dot(s_[i], q);
             for (size_t k = 0; k < n; k++)
                 q[k] -= a[i] * y_[i][k];
@@ -134,7 +137,7 @@ private:
             double yhy = 0;
             for (size_t k = 0; k < n; k++)
                 yhy += y_[m - 1][k] * y_[m - 1][k] / h_[k];
-            gamma = Dot(s_[m - 1], y_[m - 1]) / yhy;
+            gamma = sy_[m - 1] / yhy;
         }
         for (size_t k = 0; k < n; k++)
             q[k] *= gamma / h_[k];
@@ -149,5 +152,6 @@ private:
 
     size_t memory_;
     std::deque<std::vector<double>> s_, y_;
+    std::deque<double> sy_;   // s.y of each pair, as taken when it was kept
     std::vector<double> h_;
 };
diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp
index 1ce87e07d..3cf0506bb 100644
--- a/image_analysis/structure_refinement/RefinementXray.cpp
+++ b/image_analysis/structure_refinement/RefinementXray.cpp
@@ -282,7 +282,14 @@ gemmi::Grid<float> BulkSolventMask::MaskGrid(const gemmi::UnitCell &cell, const
     grid.unit_cell = cell;
     grid.spacegroup = &sg;
     const double spacing = std::max(std::min(d_min / 3.0, MASK_SPACING_MAX), std::cbrt(cell.volume / MASK_POINTS_MAX));
-    grid.set_size_from_spacing(spacing, gemmi::GridSizeRounding::Up);
+    // set_size_from_spacing()'s size, without allocating the values
+    const std::array<int, 3> n = gemmi::good_grid_size({{cell.a / spacing, cell.b / spacing, cell.c / spacing}},
+                                                       gemmi::GridSizeRounding::Up, &sg);
+    grid.nu = n[0];
+    grid.nv = n[1];
+    grid.nw = n[2];
+    grid.calculate_spacing();
+    grid.axis_order = gemmi::AxisOrder::XYZ;
     return grid;
 }
 
@@ -291,6 +298,7 @@ void BulkSolventMask::Fmask(const gemmi::Model &model, const std::vector<gemmi::
     if (orbit_leaders_.empty())   // at the first mask: a GPU engine's refinement never draws one here
         orbit_leaders_ = OrbitLeaders(grid_, nthreads_);
     gemmi::Grid<float> mask = grid_;
+    mask.data.resize(mask.point_count());
     PutMaskOnGrid(mask, model, orbit_leaders_, nthreads_);
     const gemmi::FPhiGrid<float> fm = MapToFPhi(mask);
     f.clear();