v1.0.0-rc.172 (#82)

* Fixed `jfjoch_broker` cancelling every data collection with a CUDA "out of memory" error after long operation: GPU memory no longer leaks with each collection.
* Rugnux scales a rotation sweep until the per-frame scales settle instead of for a fixed three rounds, and says so when they did not - merged intensities, and the space group, resolution cut and frame rejection read off them, change accordingly; `--scaling-iterations` is now the cap on that loop (default 100).
* Rugnux places every frame of a marCCD, SMV or miniCBF series at the spindle angle its own header states, so a series with missing frames, or with angles written modulo 360, is no longer read at the wrong geometry or refused.
* Every rotation run writes two diagnostic files beside its reflections: `<prefix>_detector.jpg`, the detector projection with the pixel mask and the detected beam-stop shadow drawn on it, and `<prefix>_plot.txt`, one row per image.

Reviewed-on: #82
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

tests/ModelValidationTest.cpp

@@ -7,16 +7,20 @@
 #include <cstdio>
 #include <filesystem>
 #include <fstream>
+#include <map>
 #include <random>
 #include <sstream>
 
 #include <gemmi/mmread_gz.hpp>
+#include <gemmi/fourier.hpp>
 
 #include "../common/Logger.h"
+#include "../rugnux/ModelFFT.h"
 #include "../rugnux/ModelValidation.h"
 #include "../rugnux/RigidBodyRefine.h"
 #include "../rugnux/SigmaA.h"
 #include "../rugnux/WriteModel.h"
+#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
 
 namespace {
     // A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with
@@ -393,6 +397,82 @@ TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidat
     std::filesystem::remove(other_setting);
 }
 
+// An alternative indexing is settled by relabelling the DATA into the model's indexing, not by moving
+// the model into the data's: that is what puts every dataset of one crystal form in one convention.
+// Only where a reference has already fixed the data's indexing is the model moved instead.
+TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValidation]") {
+    Logger logger("ModelValidation_ReindexesTheDataIntoTheModelsIndexing");
+
+    // Point group 4 on a tetragonal lattice (4/mmm): one alternative indexing.
+    const auto model = WriteTemp("reidx_model_test.pdb",
+        ClusterPdb("CRYST1   34.000   34.000   38.000  90.00  90.00  90.00 P 4           4\n").c_str());
+    auto own = ModelReferenceIntensities(model, {}, {}, 2.5, logger);
+    REQUIRE(own.size() > 500);
+    for (size_t i = 0; i < own.size(); i++) {
+        own[i].F = std::sqrt(std::max(0.0f, own[i].I));
+        own[i].rfree_flag = (i % 20) == 0;
+    }
+    const UnitCell cell{.a = 34, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
+    const gemmi::SpaceGroup *p4 = gemmi::find_spacegroup_by_name("P 4");
+    const auto laws = ReindexAmbiguityOperators(cell, *p4);
+    REQUIRE(laws.size() == 1);
+    const std::string prefix = (std::filesystem::temp_directory_path() / "reidx_test").string();
+
+    // The data as a run that picked the other indexing would have merged them.
+    const auto obs = ReindexReflections(own, laws.front());
+
+    const auto to_model = ValidateAgainstModel(obs, cell, model, prefix, logger, p4,
+                                               /*probe_indexing_ambiguity=*/true, 1, 1.0);
+    REQUIRE(to_model.ok);
+    CHECK(to_model.change_of_basis_op == gemmi::Op::identity());
+    CHECK(to_model.indexing_decided);
+    CHECK(to_model.indexing_op == laws.front());
+    CHECK(to_model.r_work < 0.15);
+
+    // A reference already fixed the data's indexing: the data stay, and the model is moved.
+    const auto fixed = ValidateAgainstModel(obs, cell, model, prefix, logger, p4, false, 1, 1.0);
+    REQUIRE(fixed.ok);
+    CHECK_FALSE(fixed.change_of_basis_op == gemmi::Op::identity());
+    CHECK(fixed.indexing_op == gemmi::Op::identity());
+    CHECK(fixed.r_work < 0.15);
+
+    // A near-perfect twin of the mis-indexed data, 45 % of it in the model's indexing: the model
+    // prefers that indexing by less than a model in a random orientation prefers one, so it has decided
+    // nothing and the data must keep the indexing they were merged in.
+    std::vector<MergedReflection> twinned = obs;
+    {
+        const gemmi::GroupOps gops = p4->operations();
+        const gemmi::ReciprocalAsu asu(p4);
+        auto key = [&](const gemmi::Miller &h) { return asu.to_asu(h, gops).first; };
+        std::map<gemmi::Miller, float> by_hkl;
+        for (const auto &r : obs)
+            by_hkl[key({{r.h, r.k, r.l}})] = r.I;
+        for (auto &r : twinned) {
+            const auto mate = by_hkl.find(key(laws.front().apply_to_hkl({{r.h, r.k, r.l}})));
+            REQUIRE(mate != by_hkl.end());
+            r.I = 0.55f * r.I + 0.45f * mate->second;
+            r.F = std::sqrt(std::max(0.0f, r.I));
+        }
+    }
+    const auto twin = ValidateAgainstModel(twinned, cell, model, prefix, logger, p4, true, 1, 1.0);
+    REQUIRE(twin.ok);
+    CHECK_FALSE(twin.indexing_decided);
+    CHECK(twin.indexing_op == gemmi::Op::identity());
+    CHECK(twin.change_of_basis_op == gemmi::Op::identity());
+
+    // Data already in the model's indexing: nothing moves, and there is nothing to arbitrate.
+    const auto same = ValidateAgainstModel(own, cell, model, prefix, logger, p4, true, 1, 1.0);
+    REQUIRE(same.ok);
+    CHECK(same.change_of_basis_op == gemmi::Op::identity());
+    CHECK(same.indexing_op == gemmi::Op::identity());
+    CHECK_FALSE(same.fit_tested);
+    CHECK(same.r_work < 0.15);
+
+    for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
+        std::filesystem::remove(prefix + suffix);
+    std::filesystem::remove(model);
+}
+
 TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]") {
     Logger logger("WriteModel_KeepsTheContentAndTakesTheGivenFrame");
 
@@ -587,3 +667,147 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]")
         std::filesystem::remove(prefix + suffix);
     std::filesystem::remove(path);
 }
+
+// The model path's structure factors come from FFTW; they must be gemmi's own transform to float
+// precision, in the same layout, so prepare_asu_data() reads the same reflections from either.
+TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") {
+    // An even grid, and an odd one on every axis (FFTW and pocketfft split odd lengths differently).
+    const auto size = GENERATE(std::array<int, 3>{20, 24, 30}, std::array<int, 3>{15, 21, 27});
+    gemmi::Grid<float> map;
+    map.unit_cell.set(40.0, 50.0, 60.0, 90.0, 95.0, 90.0);
+    map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
+    map.set_size(size[0], size[1], size[2]);
+    std::mt19937 rng(7);
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    for (auto &x : map.data)
+        x = dist(rng);
+
+    gemmi::FPhiGrid<float> ref = gemmi::transform_map_to_f_phi(map, true);
+    gemmi::FPhiGrid<float> ours = MapToFPhi(map);
+    REQUIRE(ours.nu == ref.nu);
+    REQUIRE(ours.nv == ref.nv);
+    REQUIRE(ours.nw == ref.nw);
+    REQUIRE(ours.half_l == ref.half_l);
+    REQUIRE(ours.data.size() == ref.data.size());
+    double largest = 0, worst = 0;
+    for (size_t i = 0; i < ref.data.size(); i++) {
+        largest = std::max(largest, static_cast<double>(std::abs(ref.data[i])));
+        worst = std::max(worst, static_cast<double>(std::abs(ours.data[i] - ref.data[i])));
+    }
+    CHECK(worst <= 1e-5 * largest);
+
+    const auto a = ref.prepare_asu_data(4.0, 0, false, false, false);
+    const auto b = ours.prepare_asu_data(4.0, 0, false, false, false);
+    REQUIRE(a.v.size() == b.v.size());
+    for (size_t i = 0; i < a.v.size(); i++)
+        CHECK(a.v[i].hkl == b.v[i].hkl);
+}
+
+// The output maps come from FFTW too; each must be gemmi's own map to float precision, point for point
+// in the same layout. Coefficients with arbitrary phases, negative indices, expanded by symmetry and
+// Friedel, on grids that are odd along u and v (w is even by construction of a half-l grid).
+TEST_CASE("ModelValidation_MapFromFPhiMatchesGemmi", "[ModelValidation]") {
+    const char *sg_name = GENERATE("P 1", "P 1 21 1", "P 21 21 21");
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(sg_name);
+    REQUIRE(sg != nullptr);
+    gemmi::AsuData<std::complex<float>> coef;
+    coef.unit_cell_.set(31.0, 43.0, 57.0, 90.0, sg->number == 4 ? 104.0 : 90.0, 90.0);
+    coef.spacegroup_ = sg;
+    const gemmi::ReciprocalAsu asu(sg);
+    const gemmi::GroupOps gops = sg->operations();
+    std::mt19937 rng(11);
+    std::uniform_real_distribution<float> amp(0.1f, 10.0f), phase(-3.14159f, 3.14159f);
+    for (int h = -7; h <= 7; h++)
+        for (int k = -9; k <= 9; k++)
+            for (int l = -11; l <= 11; l++) {
+                const gemmi::Op::Miller hkl{{h, k, l}};
+                if ((h == 0 && k == 0 && l == 0) || !asu.is_in(hkl) || gops.is_systematically_absent(hkl))
+                    continue;
+                coef.v.push_back({hkl, std::polar(amp(rng), phase(rng))});
+            }
+    REQUIRE(coef.v.size() > 500);
+    // P 1 takes an odd grid on u and v; the screw axes need even factors.
+    const std::array<int, 3> size = sg->number == 1 ? std::array<int, 3>{17, 21, 26}
+                                                     : std::array<int, 3>{18, 24, 26};
+    gemmi::FPhiGrid<float> grid = gemmi::get_f_phi_on_grid<float>(coef, size, true);
+    grid.data[grid.index_n(2, -3, 4)] = std::complex<float>(1.0f, NAN);   // a missing coefficient
+
+    const gemmi::Grid<float> ours = MapFromFPhi(grid);
+    const gemmi::Grid<float> ref = gemmi::transform_f_phi_grid_to_map(gemmi::FPhiGrid<float>(grid));
+    REQUIRE(ours.nu == ref.nu);
+    REQUIRE(ours.nv == ref.nv);
+    REQUIRE(ours.nw == ref.nw);
+    REQUIRE(ours.axis_order == ref.axis_order);
+    REQUIRE(ours.spacegroup == ref.spacegroup);
+    REQUIRE(ours.data.size() == ref.data.size());
+    double largest = 0, worst = 0;
+    for (size_t i = 0; i < ref.data.size(); i++) {
+        REQUIRE(std::isfinite(ours.data[i]));
+        largest = std::max(largest, static_cast<double>(std::abs(ref.data[i])));
+        worst = std::max(worst, static_cast<double>(std::abs(ours.data[i] - ref.data[i])));
+    }
+    CHECK(largest > 0);
+    CHECK(worst <= 1e-5 * largest);
+}
+
+// Map -> coefficients -> map is the identity (the V/N and 1/V scales cancel the unnormalised
+// transforms), on a grid odd along u and v.
+TEST_CASE("ModelValidation_ModelFFTRoundTrip", "[ModelValidation]") {
+    gemmi::Grid<float> map;
+    map.unit_cell.set(35.0, 45.0, 55.0, 80.0, 95.0, 105.0);
+    map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
+    map.set_size(15, 21, 28);
+    std::mt19937 rng(3);
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    for (auto &x : map.data)
+        x = dist(rng);
+
+    const gemmi::Grid<float> back = MapFromFPhi(MapToFPhi(map));
+    REQUIRE(back.nu == map.nu);
+    REQUIRE(back.nv == map.nv);
+    REQUIRE(back.nw == map.nw);
+    double worst = 0;
+    for (size_t i = 0; i < map.data.size(); i++)
+        worst = std::max(worst, static_cast<double>(std::abs(back.data[i] - map.data[i])));
+    CHECK(worst <= 1e-5);
+}
+
+// The Jacobian's six columns are evaluated in parallel; the placement must be the serial one, bit for
+// bit.
+TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation]") {
+    Logger logger("ModelValidation_RigidBodySameOnAnyNumberOfThreads");
+
+    const auto path = WriteTemp("rigid_body_threads_test.pdb", ClusterPdb().c_str());
+    gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
+    const gemmi::SpaceGroup *sg = st.find_spacegroup();
+    REQUIRE(sg != nullptr);
+    st.setup_cell_images();
+
+    const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
+    REQUIRE_FALSE(ref.empty());
+    gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
+    fobs.unit_cell_ = st.cell;
+    fobs.spacegroup_ = sg;
+    for (const auto &r : ref)
+        fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
+    fobs.ensure_sorted();
+
+    std::vector<gemmi::Position> displaced;
+    for (const gemmi::Position &p : ModelPositions(st.models[0]))
+        displaced.emplace_back(p.x + 0.40, p.y - 0.30, p.z + 0.20);
+
+    gemmi::Model serial = st.models[0], parallel = st.models[0];
+    SetModelPositions(serial, displaced);
+    SetModelPositions(parallel, displaced);
+    const RigidBodyRefineResult r1 = RefineRigidBody(serial, st.cell, *sg, fobs, 3.0, logger, 1);
+    const RigidBodyRefineResult r6 = RefineRigidBody(parallel, st.cell, *sg, fobs, 3.0, logger, 6);
+    CHECK(r1.evaluations == r6.evaluations);
+    CHECK(r1.angle_deg == r6.angle_deg);
+    CHECK(r1.shift_A == r6.shift_A);
+    const auto p1 = ModelPositions(serial), p6 = ModelPositions(parallel);
+    REQUIRE(p1.size() == p6.size());
+    for (size_t i = 0; i < p1.size(); i++)
+        CHECK((p1[i].x == p6[i].x && p1[i].y == p6[i].y && p1[i].z == p6[i].z));
+
+    std::filesystem::remove(path);
+}