From 8d199cceba5fea85319667812be4b05bb40ffbcd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 18 Oct 2026 17:07:13 +0200 Subject: [PATCH] ModelValidationTest: the CPU density is gemmi's to rounding, the same on any number of threads The density's walk now places its points in single precision (the Metal density series), so it is no longer gemmi's bit for bit: the two gridding tests compare it to gemmi within 1e-5 relative plus four cutoff-sized terms (a point at a radius' edge can fall on the other side), and check that 1 and 4 threads give the same grid bit for bit. The solvent mask is still compared to gemmi's bit for bit. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- tests/ModelValidationTest.cpp | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 7b8e40e4a..1697896d9 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -1076,8 +1076,10 @@ TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation } // The rigid body puts each probe placement on the grid with its own parallel copy of gemmi's gridding -// (density, solvent mask and their symmetrization). It must give gemmi's grids bit for bit, on any -// number of threads, for isotropic and anisotropic atoms and for groups with and without centring. +// (density, solvent mask and their symmetrization). The mask must be gemmi's bit for bit; the density +// gemmi's but for the last bits of its single-precision walk (a point at a radius' edge can also fall on +// the other side, adding or losing a term at the cutoff). Both must be the same on any number of threads, +// for isotropic and anisotropic atoms and for groups with and without centring. TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") { const char *crysts[] = { "CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n", @@ -1103,6 +1105,15 @@ TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") { return a.data.size() == b.data.size() && std::memcmp(a.data.data(), b.data.data(), a.data.size() * sizeof(float)) == 0; }; + // gemmi's density to rounding: relative 1e-5, and the few cutoff-sized terms at a radius' edge. + const auto close = [](const gemmi::Grid &a, const gemmi::Grid &b, float cutoff) { + if (a.data.size() != b.data.size()) + return false; + for (size_t i = 0; i < a.data.size(); i++) + if (std::fabs(a.data[i] - b.data[i]) > 1e-5f * std::fabs(b.data[i]) + 4 * cutoff) + return false; + return true; + }; for (double d_min : {6.0, 3.5}) { gemmi::DensityCalculator, float> gemmi_dc; gemmi_dc.d_min = d_min; @@ -1118,6 +1129,7 @@ TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") { gemmi_mask.set_size_from_spacing(gemmi_dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); gemmi::SolventMasker(gemmi::AtomicRadiiSet::Refmac).put_mask_on_grid(gemmi_mask, st.models[0]); + std::vector one_thread; for (size_t nthreads : {1, 4}) { const std::vector leaders = OrbitLeaders(gemmi_mask, nthreads); @@ -1128,7 +1140,11 @@ TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") { dc.grid.spacegroup = sg; dc.set_refmac_compatible_blur(st.models[0]); PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads); - CHECK(same(dc.grid, gemmi_dc.grid)); + CHECK(close(dc.grid, gemmi_dc.grid, gemmi_dc.cutoff)); + if (one_thread.empty()) + one_thread = dc.grid.data; + else + CHECK(dc.grid.data == one_thread); gemmi::Grid mask; mask.unit_cell = st.cell; @@ -1173,7 +1189,9 @@ TEST_CASE("ModelValidation_GriddingAnAtomWhoseRadiusStartsAtTheCutoff", "[ModelV gemmi_dc.put_model_density_on_grid(st.models[0]); PutModelDensityOnGrid(dc, st.models[0], {}, 4); REQUIRE(dc.grid.data.size() == gemmi_dc.grid.data.size()); - CHECK(std::memcmp(dc.grid.data.data(), gemmi_dc.grid.data.data(), dc.grid.data.size() * sizeof(float)) == 0); + // gemmi's density but for the last bits of the single-precision walk (see above). + for (size_t i = 0; i < dc.grid.data.size(); i++) + CHECK(std::fabs(dc.grid.data[i] - gemmi_dc.grid.data[i]) <= 1e-5f * std::fabs(gemmi_dc.grid.data[i]) + 4 * dc.cutoff); CHECK(std::all_of(dc.grid.data.begin(), dc.grid.data.end(), [](float v) { return std::isfinite(v); })); }