Merge branch 'perf-modelval' into rc175

image_analysis/structure_refinement/ModelMaskGPU.cu

@@ -189,6 +189,49 @@ ModelMaskGPU::ModelMaskGPU(cudaStream_t stream, size_t max_points)
                               "Solvent mask grid too large for int32 labels");
 }
 
+namespace {
+
+// gemmi Grid::calculate_spacing(): 1 / (n |a*|), where a* is row k of orth's inverse, which is the
+// cross product of the other two columns over the determinant (= volume).
+void GridSpacing(const ModelMaskGrid &grid, double spacing[3]) {
+    const int n[3] = {grid.nu, grid.nv, grid.nw};
+    const double *o = grid.orth;
+    const double col[3][3] = {{o[0], o[3], o[6]}, {o[1], o[4], o[7]}, {o[2], o[5], o[8]}};
+    for (int k = 0; k < 3; k++) {
+        const double *p = col[(k + 1) % 3], *q = col[(k + 2) % 3];
+        const double cx = p[1] * q[2] - p[2] * q[1], cy = p[2] * q[0] - p[0] * q[2], cz = p[0] * q[1] - p[1] * q[0];
+        spacing[k] = grid.volume / (n[k] * std::sqrt(cx * cx + cy * cy + cz * cz));
+    }
+}
+
+}  // namespace
+
+bool ModelMaskGPU::ShrinkIsNoOp(const ModelMaskGrid &grid) {
+    const int n[3] = {grid.nu, grid.nv, grid.nw};
+    double spacing[3];
+    GridSpacing(grid, spacing);
+    const double *o = grid.orth;
+    // gemmi's shrink, set_margin_around(r_shrink), only ever looks at the grid offsets other than 0
+    // within r_shrink, in the box of floor(r_shrink / spacing). No such offset, no change.
+    const int d[3] = {static_cast<int>(std::floor(R_SHRINK / spacing[0])),
+                      static_cast<int>(std::floor(R_SHRINK / spacing[1])),
+                      static_cast<int>(std::floor(R_SHRINK / spacing[2]))};
+    for (int w = -d[2]; w <= d[2]; w++)
+        for (int v = -d[1]; v <= d[1]; v++)
+            for (int u = -d[0]; u <= d[0]; u++) {
+                const double f[3] = {static_cast<double>(u) / n[0], static_cast<double>(v) / n[1],
+                                     static_cast<double>(w) / n[2]};
+                double r2 = 0;
+                for (int i = 0; i < 3; i++) {
+                    const double x = o[3 * i] * f[0] + o[3 * i + 1] * f[1] + o[3 * i + 2] * f[2];
+                    r2 += x * x;
+                }
+                if (r2 <= R_SHRINK * R_SHRINK && r2 != 0.0)
+                    return false;
+            }
+    return true;
+}
+
 void ModelMaskGPU::SetGrid(const ModelMaskGrid &grid, const std::vector<ModelMaskOp> &ops) {
     npoints = static_cast<size_t>(grid.nu) * grid.nv * grid.nw;
     if (npoints > max_points)
@@ -206,35 +249,10 @@ void ModelMaskGPU::SetGrid(const ModelMaskGrid &grid, const std::vector<ModelMas
         for (int j = 0; j < 3; j++)
             geom.orth_n[3 * i + j] = static_cast<float>(grid.orth[3 * i + j] / n[j]);
 
-    // gemmi Grid::calculate_spacing(): 1 / (n |a*|), where a* is row k of orth's inverse, which is the
-    // cross product of the other two columns over the determinant (= volume).
-    const double *o = grid.orth;
-    const double col[3][3] = {{o[0], o[3], o[6]}, {o[1], o[4], o[7]}, {o[2], o[5], o[8]}};
-    for (int k = 0; k < 3; k++) {
-        const double *p = col[(k + 1) % 3], *q = col[(k + 2) % 3];
-        const double cx = p[1] * q[2] - p[2] * q[1], cy = p[2] * q[0] - p[0] * q[2], cz = p[0] * q[1] - p[1] * q[0];
-        geom.spacing[k] = grid.volume / (n[k] * std::sqrt(cx * cx + cy * cy + cz * cz));
-    }
-
-    // gemmi's shrink, set_margin_around(r_shrink), only ever looks at the grid offsets other than 0
-    // within r_shrink, in the box of floor(r_shrink / spacing). No such offset, no change.
-    const int d[3] = {static_cast<int>(std::floor(R_SHRINK / geom.spacing[0])),
-                      static_cast<int>(std::floor(R_SHRINK / geom.spacing[1])),
-                      static_cast<int>(std::floor(R_SHRINK / geom.spacing[2]))};
-    for (int w = -d[2]; w <= d[2]; w++)
-        for (int v = -d[1]; v <= d[1]; v++)
-            for (int u = -d[0]; u <= d[0]; u++) {
-                const double f[3] = {static_cast<double>(u) / n[0], static_cast<double>(v) / n[1],
-                                     static_cast<double>(w) / n[2]};
-                double r2 = 0;
-                for (int i = 0; i < 3; i++) {
-                    const double x = o[3 * i] * f[0] + o[3 * i + 1] * f[1] + o[3 * i + 2] * f[2];
-                    r2 += x * x;
-                }
-                if (r2 <= R_SHRINK * R_SHRINK && r2 != 0.0)
-                    throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
-                                          "Solvent mask grid too fine: the GPU mask has no shrink step");
-            }
+    GridSpacing(grid, geom.spacing);
+    if (!ShrinkIsNoOp(grid))
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "Solvent mask grid too fine: the GPU mask has no shrink step");
 
     // gemmi SolventMasker::remove_islands(), rounded the same way.
     island_limit = static_cast<int>(static_cast<size_t>(ISLAND_MIN_VOLUME * npoints / grid.volume));

image_analysis/structure_refinement/ModelMaskGPU.h

@@ -50,8 +50,12 @@ public:
     static size_t DeviceBytes(size_t max_points);
     ModelMaskGPU(cudaStream_t stream, size_t max_points);
 
+    // Whether gemmi's shrink (r_shrink = 0.8 A) leaves this grid as it is - the only grids this mask is
+    // right on. Host code; asked before any engine is reserved, so a grid it fails goes to the CPU first.
+    static bool ShrinkIsNoOp(const ModelMaskGrid &grid);
+
     // Per zone. Throws if the grid has more than max_points points, if there are more than MAX_OPS
-    // operators, or if gemmi's shrink (r_shrink = 0.8 A) would change anything on this grid.
+    // operators, or if ShrinkIsNoOp() fails on it.
     void SetGrid(const ModelMaskGrid &grid, const std::vector<ModelMaskOp> &ops);
 
     // d_atoms: hydrogens and unoccupied atoms already left out. d_mask: the whole grid, 1 = solvent, 0 = macromolecule. Queued on the stream.

image_analysis/structure_refinement/ModelScaling.cpp

@@ -11,12 +11,93 @@
 
 namespace {
 
+// gemmi::Scaling<float> as a grid point fits it: the same parameters, model values, derivatives
+// and R. The solvent pair is fixed there, so |Fcalc + k_sol exp(-b_sol s^2) Fmask| of a reflection
+// is the same at every evaluation of the solver; it is taken once per grid point here rather than
+// at every evaluation, where it was most of the cost of the fit. Every expression is the one
+// scaling.hpp evaluates, in the same types, so every number is the same to the bit.
+struct FixedSolventFit {
+    struct Point {
+        gemmi::Miller hkl;
+        float fobs;
+        float fcalc_abs;   // std::abs(Scaling::get_fcalc(p))
+        double get_y() const { return fobs; }
+        double get_weight() const { return 1.0; }
+    };
+    const std::vector<gemmi::Vec6> &constraint_matrix;
+    double k_overall = 1.0;
+    gemmi::SMat33<double> b_star{0, 0, 0, 0, 0, 0};
+    std::vector<Point> points;
+
+    explicit FixedSolventFit(const gemmi::Scaling<float> &s) : constraint_matrix(s.constraint_matrix) {}
+
+    // The reflections at the solvent pair of `s`, and its scale as the start of the fit.
+    void Start(const gemmi::Scaling<float> &s) {
+        k_overall = s.k_overall;
+        b_star = s.b_star;
+        points.resize(s.points.size());
+        for (size_t i = 0; i < points.size(); ++i)
+            points[i] = {s.points[i].hkl, s.points[i].fobs, std::abs(s.get_fcalc(s.points[i]))};
+    }
+
+    // Scaling::get_parameters(), set_parameters(), get_overall_scale_factor(), compute_value() and
+    // compute_value_and_derivatives(), with k_sol and b_sol fixed.
+    std::vector<double> get_parameters() const {
+        std::vector<double> ret;
+        ret.push_back(k_overall);
+        for (const gemmi::Vec6 &v : constraint_matrix)
+            ret.push_back(gemmi::vec6_dot(v, b_star));
+        return ret;
+    }
+    void set_parameters(const double *p) {
+        k_overall = p[0];
+        int n = 0;
+        b_star = {0, 0, 0, 0, 0, 0};
+        for (const gemmi::Vec6 &row : constraint_matrix) {
+            double d = p[++n];
+            b_star.u11 += row[0] * d;
+            b_star.u22 += row[1] * d;
+            b_star.u33 += row[2] * d;
+            b_star.u12 += row[3] * d;
+            b_star.u13 += row[4] * d;
+            b_star.u23 += row[5] * d;
+        }
+    }
+    void set_parameters(const std::vector<double> &p) { set_parameters(p.data()); }
+    double get_overall_scale_factor(const gemmi::Miller &hkl) const {
+        return k_overall * std::exp(-0.25 * b_star.r_u_r(hkl));
+    }
+    double compute_value(const Point &p) const {
+        return p.fcalc_abs * (float) get_overall_scale_factor(p.hkl);
+    }
+    double compute_value_and_derivatives(const Point &p, std::vector<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;
+        double y = k_overall * fe;
+        dy_da[0] = fe;
+        gemmi::SMat33<double> du = {
+            -0.25 * y * (h.x * h.x),
+            -0.25 * y * (h.y * h.y),
+            -0.25 * y * (h.z * h.z),
+            -0.5 * y * (h.x * h.y),
+            -0.5 * y * (h.x * h.z),
+            -0.5 * y * (h.y * h.z),
+        };
+        for (size_t j = 0; j < constraint_matrix.size(); ++j)
+            dy_da[n + j] = gemmi::vec6_dot(constraint_matrix[j], du);
+        return y;
+    }
+};
+
 // R-factor of the current parameters over the fitted reflections. This is what the grid is
 // selected on, and it is the quantity the scale exists to make small.
-double RFactor(const gemmi::Scaling<float> &scaling) {
+double RFactor(const FixedSolventFit &fit) {
     double num = 0, den = 0;
-    for (const auto &p : scaling.points) {
-        num += std::fabs(p.fobs - scaling.compute_value(p));
+    for (const auto &p : fit.points) {
+        num += std::fabs(p.fobs - fit.compute_value(p));
         den += p.fobs;
     }
     return den > 0 ? num / den : 1.0;
@@ -66,25 +147,32 @@ ModelScaleReport FitModelScale(gemmi::Scaling<float> &scaling, ModelScaleBox box
     gemmi::SMat33<double> best_b_star{0, 0, 0, 0, 0, 0};
 
     struct PointFit { double k_sol, b_sol, r = NAN, k_overall = 1.0; gemmi::SMat33<double> b_star{0, 0, 0, 0, 0, 0}; };
-    auto fit_point = [](gemmi::Scaling<float> &s, PointFit &pf) {
+    auto fit_point = [](gemmi::Scaling<float> &s, FixedSolventFit &fit, PointFit &pf) {
         s.k_sol = pf.k_sol;
         s.b_sol = pf.b_sol;
         s.fit_isotropic_b_approximately();   // a fresh starting point for this solvent pair
-        s.fit_parameters();                  // k_overall + anisotropic B only
-        pf.r = RFactor(s);
-        pf.k_overall = s.k_overall;
-        pf.b_star = s.b_star;
+        fit.Start(s);
+        gemmi::LevMar levmar;
+        levmar.fit(fit);                     // k_overall + anisotropic B only
+        s.k_overall = fit.k_overall;
+        s.b_star = fit.b_star;
+        pf.r = RFactor(fit);
+        pf.k_overall = fit.k_overall;
+        pf.b_star = fit.b_star;
     };
     auto try_points = [&](std::vector<PointFit> &pts) {
         if (independent)
             ParallelChunks(static_cast<int>(pts.size()), nthreads, [&](int lo, int hi) {
                 gemmi::Scaling<float> local = scaling;
+                FixedSolventFit fit(local);
                 for (int i = lo; i < hi; ++i)
-                    fit_point(local, pts[i]);
+                    fit_point(local, fit, pts[i]);
             });
-        else
+        else {
+            FixedSolventFit fit(scaling);
             for (auto &pf : pts)
-                fit_point(scaling, pf);
+                fit_point(scaling, fit, pf);
+        }
         for (const auto &pf : pts) {
             ++report.n_grid;
             // A diverged fit gives r = NaN; latched as best_r it wins every later r < best_r.

image_analysis/structure_refinement/ModelValidation.cpp

@@ -537,20 +537,29 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         dc.rate = 1.5;
         dc.set_grid_cell_and_spacegroup(ms.st);
         dc.set_refmac_compatible_blur(ms.st.models[0]);
+
+        // The mask is made on a thread of its own, beside the density: both only read the model, each
+        // fills a grid of its own, and between them they are the whole cost of this step.
+        // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
+        std::future<gemmi::AsuData<std::complex<float>>> fmask = std::async(std::launch::async,
+                [&model = ms.st.models[0], unit_cell = dc.grid.unit_cell, spacegroup = dc.grid.spacegroup,
+                 spacing = dc.requested_grid_spacing(), to_d] {
+            gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
+            gemmi::Grid<float> mask_grid;
+            mask_grid.unit_cell = unit_cell;
+            mask_grid.spacegroup = spacegroup;
+            mask_grid.set_size_from_spacing(spacing, gemmi::GridSizeRounding::Up);
+            masker.put_mask_on_grid(mask_grid, model);
+            return MapToFPhi(mask_grid).prepare_asu_data(to_d, 0);
+        });
         dc.put_model_density_on_grid(ms.st.models[0]);
         ms.fcalc = MapToFPhi(dc.grid)
                        .prepare_asu_data(dc.d_min, dc.blur, false, false, false);
-
-        // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
-        gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
-        gemmi::Grid<float> mask_grid;
-        mask_grid.unit_cell = dc.grid.unit_cell;
-        mask_grid.spacegroup = dc.grid.spacegroup;
-        mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
-        masker.put_mask_on_grid(mask_grid, ms.st.models[0]);
-        ms.fmask = MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0);
+        ms.fmask = fmask.get();
     };
-    compute_model_factors(mdl, d_min);
+    // The model's structure factors are made beside the reservation of the GPU engines below, which
+    // does not need them.
+    std::future<void> model_factors = std::async(std::launch::async, [&] { compute_model_factors(mdl, d_min); });
 
     // The rigid body's device, decided once for the whole validation - the real fit and every replicate
     // of the null on the same one, so that a verdict never depends on how full the card was when a
@@ -564,13 +573,16 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         for (const MergedReflection &r : obs)
             if (!std::isnan(r.F) && r.d >= finest_zone)
                 ++zone_observations;
+        // One engine for the real model and one for each replicate of the null placed beside it.
         rigid_body_engines = RigidBodyGPUPool::Create(st.models[0], ucell, *sg, d_min, zone_observations,
-                                                      std::min<size_t>(std::max<size_t>(nthreads, 1), 4), logger);
+                                                      std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES) + 1,
+                                                      logger);
         rigid_body_pool = rigid_body_engines.get();
     }
 #else
     (void) rigid_body_gpu;
 #endif
+    model_factors.get();
 
     gemmi::GroupOps gops = sg->operations();
     gemmi::ReciprocalAsu asu(sg);
@@ -726,11 +738,18 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     };
     auto probe_indexing = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in) {
         IndexingProbe out;
+        // Each relabelling is fitted on a thread of its own: the fits only read the structure factors.
+        std::vector<std::future<double>> relabelled;
+        for (const auto &op : reindex_ops)
+            relabelled.push_back(std::async(std::launch::async, [&, op] {
+                return fit_model(ms, ReindexReflections(obs_in, op)).r_free;
+            }));
         out.identity_fit = fit_model(ms, obs_in);
         std::vector<double> r_free{out.identity_fit->r_free};   // identity first, then the twin laws
         double best_r_free = r_free.front();
-        for (const auto &op : reindex_ops) {
-            const double cand = fit_model(ms, ReindexReflections(obs_in, op)).r_free;
+        for (size_t i = 0; i < reindex_ops.size(); i++) {
+            const gemmi::Op &op = reindex_ops[i];
+            const double cand = relabelled[i].get();
             r_free.push_back(cand);
             if (cand < best_r_free) {
                 best_r_free = cand;
@@ -762,7 +781,6 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     if (obs_model == &obs)
         identity_fit = std::move(indexing.identity_fit);
     indexing.identity_fit.reset();
-    Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit), d_min);
 
     // --- is there anything for the model to decide? ---
     // There are three: the space-group label, where the model asserts the other enantiomorph; the
@@ -777,12 +795,11 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     const bool decision_pending = result.model_enantiomorph_candidate
                                   || !(indexing.op == gemmi::Op::identity())
                                   || !(result.change_of_basis_op == gemmi::Op::identity());
-    // A null that cannot be built is a question that could not be put, and this design already has a
-    // state for that. Without this the run dies here - model validation runs BEFORE the reflection
-    // files are written, so one failed replicate would take the .mtz, .cif, .hkl and _unmerged.mtz
-    // with it, after the merge has already been paid for. The replicates are rotated models fed to a
-    // scaling path that throws on data it cannot pair up, so this is the adversarial input for it.
-    try {
+    std::vector<gemmi::Mat33> equivalent;      // the orientations the crystal cannot tell from the model's
+    double reach_deg = 0;
+    std::vector<gemmi::Mat33> null_rotation;
+    double null_d_min = d_min;
+    std::vector<MergedReflection> obs_null;
     if (decision_pending) {
         // --- the null: this model, this data, in random orientations ---
         // An R-factor on its own says nothing about whether a model belongs to a crystal. What a model
@@ -821,12 +838,10 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         std::vector<gemmi::Op> lattice_ops{gemmi::Op::identity()};
         for (const gemmi::Op &law : ReindexAmbiguityOperators(cell, ambiguity_sg))
             lattice_ops.push_back(law);
-        std::vector<gemmi::Mat33> equivalent;
         for (const gemmi::Op &law : lattice_ops)
             for (const gemmi::Op &s : gops.sym_ops)
                 equivalent.push_back(cartesian_rotation(s.combine(law), ucell));
-        const double reach_deg = RigidBodyReachDeg(st.models[0], d_min);
-        std::vector<gemmi::Mat33> null_rotation;
+        reach_deg = RigidBodyReachDeg(st.models[0], d_min);
         {
             std::mt19937 rng(NULL_SEED);
             int redraws = 0;
@@ -857,15 +872,79 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
             if (r.rfree_flag && !std::isnan(r.F))
                 free_d.push_back(r.d);
         std::sort(free_d.begin(), free_d.end(), std::greater<>());
-        double null_d_min = std::max(d_min, NULL_D_MIN);
+        null_d_min = std::max(d_min, NULL_D_MIN);
         if (free_d.size() < NULL_MIN_FREE)
             null_d_min = d_min;
         else
             null_d_min = std::min<double>(null_d_min, free_d[NULL_MIN_FREE - 1]);
-        std::vector<MergedReflection> obs_null;
         for (const MergedReflection &r : obs)
             if (r.d >= null_d_min)
                 obs_null.push_back(r);
+    }
+    const gemmi::Structure null_model = st;   // the model as read, which the replicates are made from
+    // Each replicate works on its own copy of the model as read and of its structure factors, and the
+    // replicates share nothing else, so they run concurrently, and beside the real model's placement
+    // below, which moves the real model - the one the maps and the atom-density readout are still
+    // taken from - and nothing a replicate reads.
+    //
+    // The check before the fit is only a filter: the rigid body can walk further than its nominal
+    // reach, and a replicate it carries onto an orientation equivalent to the model's has been
+    // refined into the model's own solution (measured: 24 deg, R-free 0.56 -> 0.39, which took the
+    // null's spread from 0.016 to 0.066 and a fitting model below the gate). So where the
+    // replicate ENDS is checked as well, and one that ended within reach is drawn again and
+    // placed again - from a generator of its own, so the draws a replicate takes do not depend on
+    // the order the concurrent replicates finish in.
+    std::vector<double> null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES);
+    std::atomic<int> placed_redraws{0};
+    auto run_replicate = [&](int i) {
+        std::mt19937 rng(NULL_SEED + 1 + i);
+        gemmi::Mat33 rot = null_rotation[i];
+        for (int redraw = 0;; redraw++) {
+            ModelState rep;
+            rep.st = null_model;
+            reorient_about_centroid(rep.st.models[0], rot);
+            compute_model_factors(rep, null_d_min);
+            std::optional<Fit> identity_fit;
+            if (!reindex_ops.empty()) {
+                IndexingProbe probe = probe_indexing(rep, obs_null);
+                null_margin[i] = probe.margin;
+                identity_fit = std::move(probe.identity_fit);   // replicates are placed against obs
+            }
+            const Placement placed = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min);
+            null_r_work[i] = placed.fit.r_work;
+            const gemmi::Mat33 ended = placed.rb_rotation.multiply(rot);
+            if (redraw >= NULL_MAX_REDRAWS || angle_to_nearest_deg(ended, equivalent) >= reach_deg)
+                break;
+            ++placed_redraws;
+            rot = random_rotation(rng);
+            for (int k = 0; k < NULL_MAX_REDRAWS && angle_to_nearest_deg(rot, equivalent) < reach_deg; k++)
+                rot = random_rotation(rng);
+        }
+    };
+    // On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool
+    // worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of
+    // every replicate but the one on the calling thread ran serially, and the null lasted as long
+    // as its slowest serial replicate. From threads outside the pool those passes spread over the
+    // pool as the real fit's do. Each pass splits its work the same way wherever it runs, so the
+    // numbers are the same.
+    std::atomic<int> next_replicate{0};
+    std::vector<std::future<void>> replicate_threads;
+    if (decision_pending)
+        for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
+            replicate_threads.push_back(std::async(std::launch::async, [&] {
+                for (int i = next_replicate++; i < NULL_REPLICATES; i = next_replicate++)
+                    run_replicate(i);
+            }));
+    // The real model is placed while they run; what is compared with them is read off it below.
+    Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit), d_min);
+
+    // A null that cannot be built is a question that could not be put, and this design already has a
+    // state for that. Without this the run dies here - model validation runs BEFORE the reflection
+    // files are written, so one failed replicate would take the .mtz, .cif, .hkl and _unmerged.mtz
+    // with it, after the merge has already been paid for. The replicates are rotated models fed to a
+    // scaling path that throws on data it cannot pair up, so this is the adversarial input for it.
+    try {
+    if (decision_pending) {
         // The real model's side of it: the indexing margin from the model as read, as a replicate's
         // is, and the R-work of the model where its placement left it, against the reflections that
         // placement was fitted to.
@@ -889,58 +968,6 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         }
         double real_r_work = coarse_r_work(*obs_model);
 
-        // Each replicate works on its own copy of the model and of its structure factors, and the
-        // replicates share nothing else, so they run concurrently and the real model - which the
-        // maps and the atom-density readout below are still taken from - is never moved.
-        //
-        // The check before the fit is only a filter: the rigid body can walk further than its nominal
-        // reach, and a replicate it carries onto an orientation equivalent to the model's has been
-        // refined into the model's own solution (measured: 24 deg, R-free 0.56 -> 0.39, which took the
-        // null's spread from 0.016 to 0.066 and a fitting model below the gate). So where the
-        // replicate ENDS is checked as well, and one that ended within reach is drawn again and
-        // placed again - from a generator of its own, so the draws a replicate takes do not depend on
-        // the order the concurrent replicates finish in.
-        std::vector<double> null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES);
-        std::atomic<int> placed_redraws{0};
-        auto run_replicate = [&](int i) {
-            std::mt19937 rng(NULL_SEED + 1 + i);
-            gemmi::Mat33 rot = null_rotation[i];
-            for (int redraw = 0;; redraw++) {
-                ModelState rep;
-                rep.st = mdl.st;
-                SetModelPositions(rep.st.models[0], as_read);
-                reorient_about_centroid(rep.st.models[0], rot);
-                compute_model_factors(rep, null_d_min);
-                std::optional<Fit> identity_fit;
-                if (!reindex_ops.empty()) {
-                    IndexingProbe probe = probe_indexing(rep, obs_null);
-                    null_margin[i] = probe.margin;
-                    identity_fit = std::move(probe.identity_fit);   // replicates are placed against obs
-                }
-                const Placement placed = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min);
-                null_r_work[i] = placed.fit.r_work;
-                const gemmi::Mat33 ended = placed.rb_rotation.multiply(rot);
-                if (redraw >= NULL_MAX_REDRAWS || angle_to_nearest_deg(ended, equivalent) >= reach_deg)
-                    break;
-                ++placed_redraws;
-                rot = random_rotation(rng);
-                for (int k = 0; k < NULL_MAX_REDRAWS && angle_to_nearest_deg(rot, equivalent) < reach_deg; k++)
-                    rot = random_rotation(rng);
-            }
-        };
-        // On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool
-        // worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of
-        // every replicate but the one on the calling thread ran serially, and the null lasted as long
-        // as its slowest serial replicate. From threads outside the pool those passes spread over the
-        // pool as the real fit's do. Each pass splits its work the same way wherever it runs, so the
-        // numbers are the same.
-        std::atomic<int> next_replicate{0};
-        std::vector<std::future<void>> replicate_threads;
-        for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
-            replicate_threads.push_back(std::async(std::launch::async, [&] {
-                for (int i = next_replicate++; i < NULL_REPLICATES; i = next_replicate++)
-                    run_replicate(i);
-            }));
         for (std::future<void> &f : replicate_threads)
             f.get();
         if (placed_redraws > 0)
@@ -1235,24 +1262,18 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     }
 
     // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) ---
-    const gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc);
-    const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4");
-    write_ccp4(map_from_coefficients(mapfofc), output_prefix + "_fofc.ccp4");
-    {
-        double s = 0; int n = 0;
-        for (gemmi::Model &m : st.models)
-            for (gemmi::Chain &ch : m.chains)
-                for (gemmi::Residue &r : ch.residues)
-                    for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; }
-        result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
-    }
+    // The three maps - this one, the difference map and the anomalous map - share nothing but what
+    // they read, so the other two are made and written on threads of their own beside this one.
+    std::future<void> fofc_map = std::async(std::launch::async, [&] {
+        write_ccp4(map_from_coefficients(mapfofc), output_prefix + "_fofc.ccp4");
+    });
 
     // --- anomalous difference map, where the merge kept the Bijvoet split ---
     // Coefficients F(+) - F(-) carried on the model phase turned back by 90 degrees. Its peaks sit on
     // the anomalous scatterers, so reading the map at each of the model's own atoms names them,
     // rather than leaving a list of coordinates for someone to look up.
     // Following ANODE, Thorn & Sheldrick (2011) J. Appl. Cryst. 44, 1285-1287
-    {
+    std::future<void> anomalous_map = std::async(std::launch::async, [&] {
         // Read from the merged reflections as they came in, and carry each one into the model's frame
         // here: the hand each Bijvoet difference belongs to is a property of the frame the merge was
         // made in, and both operators can change it. F(+) and F(-) are attached to the + index of the
@@ -1332,7 +1353,20 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
                 sites.resize(MAX_ANOMALOUS_SITES);
             result.anomalous_sites = std::move(sites);
         }
+    });
+
+    const gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc);
+    const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4");
+    {
+        double s = 0; int n = 0;
+        for (gemmi::Model &m : st.models)
+            for (gemmi::Chain &ch : m.chains)
+                for (gemmi::Residue &r : ch.residues)
+                    for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; }
+        result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
     }
+    fofc_map.get();
+    anomalous_map.get();
 
     // --- the enantiomorph, now that there is something to decide it on ---
     // Two things have to hold before the model's hand is written on these data. The model has to

image_analysis/structure_refinement/RigidBodyGPU.cpp

@@ -186,6 +186,11 @@ std::unique_ptr<RigidBodyGPUPool> RigidBodyGPUPool::Create(const gemmi::Model &m
                             "GPU's gridding at {:.1f} A", zone_d);
                 return nullptr;
             }
+            if (!RigidBodyGPUEngine::MaskSupports(zone)) {
+                logger.Info("Model validation: the rigid body runs on the CPU - the solvent mask grid at {:.1f} A "
+                            "is finer than the GPU's mask, which has no shrink step, can take", zone_d);
+                return nullptr;
+            }
             cap.atoms = std::max(cap.atoms, zone.atoms.size());
             cap.grid_points = std::max(cap.grid_points, static_cast<size_t>(zone.nu) * zone.nv * zone.nw);
             cap.complex_points = std::max(cap.complex_points, static_cast<size_t>(zone.nu / 2 + 1) * zone.nv * zone.nw);
@@ -203,7 +208,7 @@ std::unique_ptr<RigidBodyGPUPool> RigidBodyGPUPool::Create(const gemmi::Model &m
         constexpr size_t HEADROOM = 1ull << 30;
         const size_t budget = std::min(total / 4, free > HEADROOM ? free - HEADROOM : 0);
         const size_t fit = bytes > 0 ? budget / bytes : 0;
-        const size_t want = std::min<size_t>({fit, std::max<size_t>(max_engines, 1), 4});
+        const size_t want = std::min<size_t>(fit, std::max<size_t>(max_engines, 1));
         const int device = RigidBodyGPUEngine::CurrentDevice();
         for (size_t i = 0; i < want; i++)
             pool->engines_.push_back(std::make_unique<RigidBodyGPUEngine>(cap, device));

image_analysis/structure_refinement/RigidBodyGPU.cu

@@ -731,6 +731,16 @@ bool RigidBodyGPUEngine::Supports(const RigidBodyGPUZone &zone) {
     return true;
 }
 
+bool RigidBodyGPUEngine::MaskSupports(const RigidBodyGPUZone &zone) {
+    ModelMaskGrid mask_grid;
+    mask_grid.nu = zone.nu;
+    mask_grid.nv = zone.nv;
+    mask_grid.nw = zone.nw;
+    std::copy(zone.orth, zone.orth + 9, mask_grid.orth);
+    mask_grid.volume = zone.volume;
+    return ModelMaskGPU::ShrinkIsNoOp(mask_grid);
+}
+
 void RigidBodyGPUEngine::MemoryInfo(size_t &free, size_t &total) {
     cuda_err(cudaMemGetInfo(&free, &total));
 }

image_analysis/structure_refinement/RigidBodyGPU.h

@@ -42,7 +42,7 @@ public:
     // and never the last gigabyte of what is free, since the merge may be running beside it - or where
     // the cell is too small for the gather (an atom's box wider than the cell). Logged either way.
     // `max_observations`: the most working reflections any fit will be given in its finest zone.
-    // `max_engines`: at most this many, however much memory there is (1 to 4).
+    // `max_engines`: at most this many, however much memory there is.
     static std::unique_ptr<RigidBodyGPUPool> Create(const gemmi::Model &model, const gemmi::UnitCell &cell,
                                                     const gemmi::SpaceGroup &sg, double d_min,
                                                     size_t max_observations, size_t max_engines, Logger &logger);

image_analysis/structure_refinement/RigidBodyGPUEngine.h

@@ -101,6 +101,8 @@ public:
     // Whether the gather reproduces gemmi's box walk on this zone: every atom's box narrower than the cell,
     // so that no point is reached by two images of one atom.
     static bool Supports(const RigidBodyGPUZone &zone);
+    // Whether the GPU's solvent mask, which has no shrink step, is gemmi's on this zone's grid.
+    static bool MaskSupports(const RigidBodyGPUZone &zone);
     static void MemoryInfo(size_t &free, size_t &total);
     static int CurrentDevice();
 

rugnux/RugnuxScaleMerge.cpp

@@ -3234,6 +3234,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
     // with no reference MTZ - the alternative indexing. Both are relabelings of the same
     // measurements, so the statistics computed above are unaffected; what changes is that the file,
     // the R-factors and the maps then all describe one indexing instead of two.
+    std::future<void> placed_model_file;
     if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass
         && !superseded) {
         phase("Validating against model");
@@ -3338,9 +3339,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
         // relabelled to the model's enantiomorph, so "did not move" is not "identical to the input
         // file". A rejected model is scored, placed and mapped like any other, and that is the case
         // where someone most wants to see the model in the density, so it gets its file too.
+        // Written beside the reflection files below: nothing after this reads it.
         if (validation.placed_model)
-            WritePlacedModel(*validation.placed_model, *result.consensus_cell,
-                             experiment_.GetSpaceGroupOrP1(), config_.output_prefix, logger);
+            placed_model_file = std::async(std::launch::async,
+                    [placed = validation.placed_model, cell = *result.consensus_cell,
+                     sg = experiment_.GetSpaceGroupOrP1(), prefix = config_.output_prefix, &logger] {
+                WritePlacedModel(*placed, cell, sg, prefix, logger);
+            });
     }
 
     // The unmerged MTZ (written with the other output files, see there) is built from here on,
@@ -3454,5 +3459,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                         p1.size(), FinalistLedgerToText(led));
         }
     }
+    if (placed_model_file.valid())
+        placed_model_file.get();
     return false;
 }

tests/ModelMaskGPUTest.cpp

@@ -214,6 +214,9 @@ TEST_CASE("ModelMaskGPU_RefusesAGridTheShrinkWouldChange", "[ModelValidation][gp
     ModelMaskGPU mask(stream, fine.point_count());
     CHECK_NOTHROW(mask.SetGrid(MaskGridDescription(coarse), ops));
     CHECK_THROWS(mask.SetGrid(MaskGridDescription(fine), ops));
+    // The same answer ahead of time, which is what sends such a grid to the CPU before any engine is made.
+    CHECK(ModelMaskGPU::ShrinkIsNoOp(MaskGridDescription(coarse)));
+    CHECK_FALSE(ModelMaskGPU::ShrinkIsNoOp(MaskGridDescription(fine)));
 }
 
 #endif

tests/ModelScalingTest.cpp

@@ -6,8 +6,11 @@
 #include <cmath>
 #include <algorithm>
 #include <complex>
+#include <random>
+#include <utility>
 
 #include "../image_analysis/structure_refinement/ModelScaling.h"
+#include "gemmi/asumask.hpp"
 #include "gemmi/scaling.hpp"
 #include "gemmi/symmetry.hpp"
 #include "gemmi/unitcell.hpp"
@@ -164,3 +167,69 @@ TEST_CASE("ModelScaling is the same on any number of threads") {
     CHECK(serial.b_star.u22 == parallel.b_star.u22);
     CHECK(serial.b_star.u33 == parallel.b_star.u33);
 }
+
+// FitModelScale fits each solvent pair of its grid on a copy of gemmi::Scaling's target that takes the
+// solvent term once per pair (ModelScaling.cpp). On a box of a single pair it has to give exactly what
+// gemmi's own Scaling gives for that pair - to the bit, since the validation's numbers are made of it.
+TEST_CASE("ModelScaling fits a grid point exactly as gemmi's Scaling does") {
+    const gemmi::UnitCell cell(61.0, 72.0, 83.0, 90.0, 101.0, 90.0);
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 1 21 1");
+    REQUIRE(sg != nullptr);
+    const gemmi::GroupOps gops = sg->operations();
+    const gemmi::ReciprocalAsu asu(sg);
+    gemmi::Scaling<float> scaling(cell, sg);
+    scaling.use_solvent = true;
+    std::mt19937 rng(20261008);
+    std::uniform_real_distribution<double> uniform(0.0, 1.0);
+    const gemmi::SMat33<double> b_star = gemmi::SMat33<double>{22, 31, 27, 0, 3, 0}.transformed_by(cell.frac.mat);
+    const gemmi::Miller lim = cell.get_hkl_limits(2.5);
+    for (int h = -lim[0]; h <= lim[0]; h++)
+        for (int k = -lim[1]; k <= lim[1]; k++)
+            for (int l = -lim[2]; l <= lim[2]; l++) {
+                const gemmi::Miller m{{h, k, l}};
+                if ((h == 0 && k == 0 && l == 0) || cell.calculate_d(m) < 2.5 || !asu.is_in(m)
+                    || gops.is_systematically_absent(m))
+                    continue;
+                const double stol2 = cell.calculate_stol_sq(m);
+                const std::complex<float> fc(std::polar(100.0 * std::exp(-8.0 * stol2) * (0.2 + uniform(rng)),
+                                                        6.283 * uniform(rng)));
+                const std::complex<float> fm(std::polar(900.0 * std::exp(-20.0 * stol2), 6.283 * uniform(rng)));
+                const std::complex<double> total = std::complex<double>(fc)
+                                                   + 0.4 * std::exp(-50.0 * stol2) * std::complex<double>(fm);
+                const double fobs = 2.5 * std::exp(-0.25 * b_star.r_u_r(m)) * std::abs(total)
+                                    * (0.9 + 0.2 * uniform(rng));
+                gemmi::Scaling<float>::Point p{};
+                p.hkl = m;
+                p.stol2 = stol2;
+                p.fcmol = fc;
+                p.fmask = fm;
+                p.fobs = static_cast<float>(fobs);
+                p.sigma = static_cast<float>(0.05 * fobs + 0.5);
+                scaling.points.push_back(p);
+            }
+    REQUIRE(scaling.points.size() > 1000);
+
+    for (const auto &[k_sol, b_sol] : {std::pair{0.35, 40.0}, std::pair{0.20, 70.0}}) {
+        gemmi::Scaling<float> ours = scaling;
+        FitModelScale(ours, ModelScaleBox{k_sol, k_sol, b_sol, b_sol}, 4);
+
+        gemmi::Scaling<float> gemmis = scaling;
+        gemmis.fix_k_sol = true;
+        gemmis.fix_b_sol = true;
+        gemmis.k_sol = k_sol;
+        gemmis.b_sol = b_sol;
+        gemmis.fit_isotropic_b_approximately();
+        gemmis.fit_parameters();
+
+        CHECK(ours.k_sol == k_sol);
+        CHECK(ours.b_sol == b_sol);
+        CHECK(ours.k_overall == gemmis.k_overall);
+        CHECK(ours.b_star.u11 == gemmis.b_star.u11);
+        CHECK(ours.b_star.u22 == gemmis.b_star.u22);
+        CHECK(ours.b_star.u33 == gemmis.b_star.u33);
+        CHECK(ours.b_star.u12 == gemmis.b_star.u12);
+        CHECK(ours.b_star.u13 == gemmis.b_star.u13);
+        CHECK(ours.b_star.u23 == gemmis.b_star.u23);
+        CHECK(ours.calculate_r_factor() == gemmis.calculate_r_factor());
+    }
+}