diff --git a/image_analysis/structure_refinement/ModelMaskGPU.cu b/image_analysis/structure_refinement/ModelMaskGPU.cu index 5a2db4b21..dd06d96dd 100644 --- a/image_analysis/structure_refinement/ModelMaskGPU.cu +++ b/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(std::floor(R_SHRINK / spacing[0])), + static_cast(std::floor(R_SHRINK / spacing[1])), + static_cast(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(u) / n[0], static_cast(v) / n[1], + static_cast(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 &ops) { npoints = static_cast(grid.nu) * grid.nv * grid.nw; if (npoints > max_points) @@ -206,35 +249,10 @@ void ModelMaskGPU::SetGrid(const ModelMaskGrid &grid, const std::vector(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(std::floor(R_SHRINK / geom.spacing[0])), - static_cast(std::floor(R_SHRINK / geom.spacing[1])), - static_cast(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(u) / n[0], static_cast(v) / n[1], - static_cast(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(static_cast(ISLAND_MIN_VOLUME * npoints / grid.volume)); diff --git a/image_analysis/structure_refinement/ModelMaskGPU.h b/image_analysis/structure_refinement/ModelMaskGPU.h index 9f821bfb5..9d1968c63 100644 --- a/image_analysis/structure_refinement/ModelMaskGPU.h +++ b/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 &ops); // d_atoms: hydrogens and unoccupied atoms already left out. d_mask: the whole grid, 1 = solvent, 0 = macromolecule. Queued on the stream. diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp index 4fa5aee59..0f3c7b995 100644 --- a/image_analysis/structure_refinement/ModelScaling.cpp +++ b/image_analysis/structure_refinement/ModelScaling.cpp @@ -11,12 +11,93 @@ namespace { +// gemmi::Scaling 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 &constraint_matrix; + double k_overall = 1.0; + gemmi::SMat33 b_star{0, 0, 0, 0, 0, 0}; + std::vector points; + + explicit FixedSolventFit(const gemmi::Scaling &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 &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 get_parameters() const { + std::vector 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 &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 &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 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 &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 &scaling, ModelScaleBox box gemmi::SMat33 best_b_star{0, 0, 0, 0, 0, 0}; struct PointFit { double k_sol, b_sol, r = NAN, k_overall = 1.0; gemmi::SMat33 b_star{0, 0, 0, 0, 0, 0}; }; - auto fit_point = [](gemmi::Scaling &s, PointFit &pf) { + auto fit_point = [](gemmi::Scaling &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 &pts) { if (independent) ParallelChunks(static_cast(pts.size()), nthreads, [&](int lo, int hi) { gemmi::Scaling 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. diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index f7cf29601..3673483dc 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -537,20 +537,29 @@ ModelValidationResult Validate(const std::vector &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>> 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 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 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 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 &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(std::max(nthreads, 1), 4), logger); + std::min(std::max(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 &merged, }; auto probe_indexing = [&](ModelState &ms, const std::vector &obs_in) { IndexingProbe out; + // Each relabelling is fitted on a thread of its own: the fits only read the structure factors. + std::vector> 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 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 &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 &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 equivalent; // the orientations the crystal cannot tell from the model's + double reach_deg = 0; + std::vector null_rotation; + double null_d_min = d_min; + std::vector 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 &merged, std::vector lattice_ops{gemmi::Op::identity()}; for (const gemmi::Op &law : ReindexAmbiguityOperators(cell, ambiguity_sg)) lattice_ops.push_back(law); - std::vector 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 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 &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(null_d_min, free_d[NULL_MIN_FREE - 1]); - std::vector 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 null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES); + std::atomic 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 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 next_replicate{0}; + std::vector> replicate_threads; + if (decision_pending) + for (size_t t = 0; t < std::min(std::max(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 &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 null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES); - std::atomic 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 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 next_replicate{0}; - std::vector> replicate_threads; - for (size_t t = 0; t < std::min(std::max(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 &f : replicate_threads) f.get(); if (placed_redraws > 0) @@ -1235,24 +1262,18 @@ ModelValidationResult Validate(const std::vector &merged, } // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- - const gemmi::Grid 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 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 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 &merged, sites.resize(MAX_ANOMALOUS_SITES); result.anomalous_sites = std::move(sites); } + }); + + const gemmi::Grid 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 diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cpp b/image_analysis/structure_refinement/RigidBodyGPU.cpp index e6d4cddcc..d85f696fc 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cpp +++ b/image_analysis/structure_refinement/RigidBodyGPU.cpp @@ -186,6 +186,11 @@ std::unique_ptr 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(zone.nu) * zone.nv * zone.nw); cap.complex_points = std::max(cap.complex_points, static_cast(zone.nu / 2 + 1) * zone.nv * zone.nw); @@ -203,7 +208,7 @@ std::unique_ptr 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({fit, std::max(max_engines, 1), 4}); + const size_t want = std::min(fit, std::max(max_engines, 1)); const int device = RigidBodyGPUEngine::CurrentDevice(); for (size_t i = 0; i < want; i++) pool->engines_.push_back(std::make_unique(cap, device)); diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cu b/image_analysis/structure_refinement/RigidBodyGPU.cu index 1db77013a..f7b283f3e 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cu +++ b/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)); } diff --git a/image_analysis/structure_refinement/RigidBodyGPU.h b/image_analysis/structure_refinement/RigidBodyGPU.h index c4a64308a..7f2651fb5 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.h +++ b/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 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); diff --git a/image_analysis/structure_refinement/RigidBodyGPUEngine.h b/image_analysis/structure_refinement/RigidBodyGPUEngine.h index 9b93f1554..648bc00f0 100644 --- a/image_analysis/structure_refinement/RigidBodyGPUEngine.h +++ b/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(); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 969b46947..3197f3765 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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 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; } diff --git a/tests/ModelMaskGPUTest.cpp b/tests/ModelMaskGPUTest.cpp index f464ce36c..d6ca1e728 100644 --- a/tests/ModelMaskGPUTest.cpp +++ b/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 diff --git a/tests/ModelScalingTest.cpp b/tests/ModelScalingTest.cpp index 741f19d39..300f97b0f 100644 --- a/tests/ModelScalingTest.cpp +++ b/tests/ModelScalingTest.cpp @@ -6,8 +6,11 @@ #include #include #include +#include +#include #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 scaling(cell, sg); + scaling.use_solvent = true; + std::mt19937 rng(20261008); + std::uniform_real_distribution uniform(0.0, 1.0); + const gemmi::SMat33 b_star = gemmi::SMat33{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 fc(std::polar(100.0 * std::exp(-8.0 * stol2) * (0.2 + uniform(rng)), + 6.283 * uniform(rng))); + const std::complex fm(std::polar(900.0 * std::exp(-20.0 * stol2), 6.283 * uniform(rng))); + const std::complex total = std::complex(fc) + + 0.4 * std::exp(-50.0 * stol2) * std::complex(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::Point p{}; + p.hkl = m; + p.stol2 = stol2; + p.fcmol = fc; + p.fmask = fm; + p.fobs = static_cast(fobs); + p.sigma = static_cast(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 ours = scaling; + FitModelScale(ours, ModelScaleBox{k_sol, k_sol, b_sol, b_sol}, 4); + + gemmi::Scaling 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()); + } +}