mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 80 Wiki Activity
Files
5c18de45a35b1c220c8bcc4d3815096f0db85a1d
Jungfraujoch/tests/ModelValidationTest.cpp
T
leonarski_fandClaude Opus 5.5 5c18de45a3 Rigid body GPU: faster gather, zone tables kept per pool, high-priority streams 
The gather stages each brick's atoms cooperatively with the Cartesian position
of the image nearest the brick, so a point no longer wraps and transforms
every atom it visits (a narrow cell keeps the per-point image). What a zone
needs of the atoms is worked out once per pool and zone rather than once per
fit. Engine streams run at the device's highest priority, because the first
validation runs beside the P1 cross-check merge.

8t7r-sized fit (63.6k atoms, 30 evaluations, 14 Jacobians): 0.80 -> 0.45 s;
6oel-sized: 0.50 -> 0.27 s.

Tests: a non-CUDA build compiles the GPU test helpers away; the GPU-vs-CPU
residual bound is 2e-4 of <Fobs>, the resolution of gemmi's own scale fit.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
2026-09-28 19:03:30 +02:00

1386 lines
68 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <map>
#include <random>
#include <sstream>
#include <gemmi/mmread_gz.hpp>
#include <gemmi/fourier.hpp>
#include <gemmi/solmask.hpp>
#include "../common/Logger.h"
#include "../rugnux/ModelFFT.h"
#include "../rugnux/ModelGrid.h"
#include "../rugnux/ModelValidation.h"
#include "../rugnux/RigidBodyRefine.h"
#ifdef JFJOCH_USE_CUDA
#include "../rugnux/RigidBodyGPU.h"
#include "../common/CUDAWrapper.h"
#endif
#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
// atoms, a cell and a space group, which is all these tests are about. Neutral by construction:
// no real specimen's parameters are involved.
const char *kMmcif = R"(data_test
_cell.length_a 40.000
_cell.length_b 50.000
_cell.length_c 60.000
_cell.angle_alpha 90.000
_cell.angle_beta 90.000
_cell.angle_gamma 90.000
_symmetry.space_group_name_H-M 'P 1'
loop_
_atom_site.group_PDB
_atom_site.id
_atom_site.type_symbol
_atom_site.label_alt_id
_atom_site.label_atom_id
_atom_site.label_comp_id
_atom_site.label_asym_id
_atom_site.label_seq_id
_atom_site.Cartn_x
_atom_site.Cartn_y
_atom_site.Cartn_z
_atom_site.occupancy
_atom_site.B_iso_or_equiv
ATOM 1 C . CA GLY A 1 10.000 12.000 14.000 1.00 20.00
ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00
)";
const char *kPdb =
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n"
"ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n"
"ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n"
"END\n";
// A synthetic "protein": carbons filling one asymmetric unit of a small cell. The space group
// matters - in P1 the origin is free in all three directions, so |F| does not change when the
// whole content is translated and there is no translation to recover. No specimen is involved;
// the positions come out of a fixed seed.
const char *kCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 21 21 21 4\n";
// The same cell in a polar group: b is the unique axis, so the origin is free along y alone.
const char *kPolarCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 1 2 1 2\n";
std::string ClusterPdb(const char *cryst = kCryst) {
std::string pdb = cryst;
std::mt19937 rng(20260902);
std::uniform_real_distribution<double> x(2, 14), y(2, 16), z(2, 18);
char line[96];
for (int i = 1; i <= 150; i++) {
std::snprintf(line, sizeof line,
"ATOM %5d C UNK A 1 %8.3f%8.3f%8.3f 1.00 20.00 C\n",
i, x(rng), y(rng), z(rng));
pdb += line;
}
return pdb + "END\n";
}
// The same two atoms as kPdb, plus an anisotropic U, a partial occupancy and a water: a written
// model has to carry what the input carried, not only its coordinates.
const char *kPdbRich =
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n"
"ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n"
"ANISOU 1 CA GLY A 1 1000 1200 1400 100 200 300 C\n"
"ATOM 2 CB GLY A 1 12.000 14.000 16.000 0.60 25.00 C\n"
"HETATM 3 O HOH A 101 20.000 22.000 24.000 1.00 30.00 O\n"
"END\n";
std::string WriteTemp(const std::string &name, const char *content) {
std::ofstream f(name);
f << content;
f.close();
return name;
}
}
// --model used to call gemmi::read_pdb unconditionally, so a deposited model handed over as mmCIF -
// which is how the PDB serves coordinates by default - was refused outright. Both formats now go
// through the same reader, which decides on the file's content.
TEST_CASE("ModelValidation_ReadsPdbAndMmcif", "[ModelValidation]") {
Logger logger("ModelValidation_ReadsPdbAndMmcif");
const auto pdb = WriteTemp("model_validation_test.pdb", kPdb);
const auto cif = WriteTemp("model_validation_test.cif", kMmcif);
const auto from_pdb = ModelReferenceIntensities(pdb, {}, {}, 4.0, logger);
const auto from_cif = ModelReferenceIntensities(cif, {}, {}, 4.0, logger);
REQUIRE_FALSE(from_pdb.empty());
REQUIRE_FALSE(from_cif.empty());
// The same structure either way, so the same reflections come out of it.
CHECK(from_cif.size() == from_pdb.size());
// The extension is not what decides: the same mmCIF under a .pdb name still reads.
const auto misnamed = WriteTemp("model_validation_test_misnamed.pdb", kMmcif);
CHECK_FALSE(ModelReferenceIntensities(misnamed, {}, {}, 4.0, logger).empty());
std::filesystem::remove(pdb);
std::filesystem::remove(cif);
std::filesystem::remove(misnamed);
}
// A model that cannot be used must say why. Returning an empty result and logging was enough to lose
// the fact entirely: the run finished successfully with no R-free and no maps, which is exactly what
// a run without --model looks like.
TEST_CASE("ModelValidation_UnusableModelGivesAReason", "[ModelValidation]") {
Logger logger("ModelValidation_UnusableModelGivesAReason");
const auto empty_file = WriteTemp("model_validation_test_bogus.pdb", "not a coordinate file\n");
const auto result = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60,
.alpha = 90, .beta = 90, .gamma = 90},
empty_file, "", logger);
CHECK_FALSE(result.ok);
CHECK_FALSE(result.failure_reason.empty());
CHECK_THAT(result.failure_reason, Catch::Matchers::ContainsSubstring(empty_file));
const auto missing = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60,
.alpha = 90, .beta = 90, .gamma = 90},
"model_validation_test_does_not_exist.pdb", "", logger);
CHECK_FALSE(missing.ok);
CHECK_FALSE(missing.failure_reason.empty());
std::filesystem::remove(empty_file);
}
// A rigid-body step is only worth taking if it can find a shift it was not told about. The check is
// closed: the "observed" amplitudes are the model's own, so the answer is known to be zero shift.
TEST_CASE("ModelValidation_RigidBodyRecoversASmallShift", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyRecoversASmallShift");
const auto path = WriteTemp("rigid_body_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();
// "Observed" amplitudes: the model's own structure factors, so the target's minimum is exactly
// where the model started.
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();
const std::vector<gemmi::Position> original = ModelPositions(st.models[0]);
std::vector<gemmi::Position> displaced;
for (const gemmi::Position &p : original)
displaced.emplace_back(p.x + 0.40, p.y - 0.30, p.z + 0.20); // 0.54 A off
SetModelPositions(st.models[0], displaced);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
CHECK(result.converged);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
double before = 0, after = 0;
for (size_t i = 0; i < original.size(); i++) {
before += original[i].dist_sq(displaced[i]);
after += original[i].dist_sq(refined[i]);
}
before = std::sqrt(before / original.size());
after = std::sqrt(after / original.size());
logger.Info("Rigid-body test: rmsd from the truth {:.3f} A -> {:.3f} A", before, after);
CHECK(after < 0.2 * before);
std::filesystem::remove(path);
}
// A polar space group leaves the origin free along one direction: moving the whole cell content
// along it multiplies every structure factor by a phase and changes no amplitude, so the data
// cannot say where the body sits along it and the refinement must not pretend otherwise. The check
// is closed - the "observed" amplitudes are the model's own - and the model is displaced in all
// three directions at once, so the same run says both what is recovered and what is left alone.
TEST_CASE("ModelValidation_RigidBodyLeavesThePolarDirectionAlone", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyLeavesThePolarDirectionAlone");
const auto path = WriteTemp("rigid_body_polar_test.pdb", ClusterPdb(kPolarCryst).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();
// b is the unique axis of P 1 2 1, so y is the free direction and x and z are determined.
const std::vector<gemmi::Position> original = ModelPositions(st.models[0]);
std::vector<gemmi::Position> displaced;
for (const gemmi::Position &p : original)
displaced.emplace_back(p.x + 0.35, p.y + 0.50, p.z - 0.30);
SetModelPositions(st.models[0], displaced);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
CHECK(result.converged);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
gemmi::Vec3 left;
for (size_t i = 0; i < original.size(); i++)
left += gemmi::Vec3(refined[i]) - gemmi::Vec3(original[i]);
left *= 1.0 / static_cast<double>(original.size());
logger.Info("Rigid-body polar test: left over ({:.3f}, {:.3f}, {:.3f}) A", left.x, left.y, left.z);
CHECK(std::fabs(left.x) < 0.10);
CHECK(std::fabs(left.z) < 0.10);
// Along b nothing was refined away, because there is nothing there to refine.
CHECK(left.y == Catch::Approx(0.50).margin(0.02));
std::filesystem::remove(path);
}
// The bulk solvent the placement is scored through has to stay inside the range a flat solvent model
// means anything in. gemmi's own scaler is an unbounded Levenberg-Marquardt and reaches b_sol of
// hundreds or thousands of A^2, which does not corrupt a reported number here but distorts the target
// that decides where the model goes, at every one of the hundreds of evaluations.
TEST_CASE("ModelValidation_RigidBodySolventStaysPhysical", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodySolventStaysPhysical");
const auto path = WriteTemp("rigid_body_solvent_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();
// Turned right around, so the model explains nothing: that is where an unbounded solvent fit
// has nothing to hold it and runs away.
std::vector<gemmi::Position> turned = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : turned)
centre += p;
centre *= 1.0 / static_cast<double>(turned.size());
for (gemmi::Position &p : turned)
p = gemmi::Position(centre - (gemmi::Vec3(p) - centre));
SetModelPositions(st.models[0], turned);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
logger.Info("Rigid-body solvent test: k_sol {:.3f}, b_sol {:.1f} A^2", result.k_sol, result.b_sol);
CHECK(result.k_sol >= 0.10);
CHECK(result.k_sol <= 0.60);
CHECK(result.b_sol >= 10.0);
CHECK(result.b_sol <= 80.0);
std::filesystem::remove(path);
}
// sigma_A is what says how much of the model to believe, so the two ends of its range are what the
// weighting has to get right: a model that explains the data completely, and one that explains none
// of it.
TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValidation]") {
gemmi::UnitCell cell(40, 50, 60, 90, 90, 90);
std::mt19937 rng(12345);
std::normal_distribution<double> normal(0.0, 1.0);
auto weights = [&](bool agreeing) {
std::vector<SigmaAReflection> refl;
for (int i = 0; i < 2000; i++) {
SigmaAReflection r;
r.f_calc = std::fabs(normal(rng)) * 100;
r.f_obs = agreeing ? r.f_calc : std::fabs(normal(rng)) * 100;
r.inv_d2 = 0.01 + 0.2 * (i / 2000.0);
r.free = (i % 20) == 0; // the usual 5 %
refl.push_back(r);
}
return EstimateSigmaA(refl, cell);
};
const SigmaAResult perfect = weights(true);
const SigmaAResult useless = weights(false);
CHECK(perfect.mean_fom > 0.85);
CHECK(useless.mean_fom < 0.2);
CHECK(perfect.shells == 2); // 100 free reflections, 50 to a shell
// No free reflections to estimate on: the coefficients are left alone rather than weighted by a
// number that was never measured.
std::vector<SigmaAReflection> no_free;
for (int i = 0; i < 100; i++)
no_free.push_back({100.0, 100.0, 0.05, 1, false, false});
const SigmaAResult unweighted = EstimateSigmaA(no_free, cell);
CHECK(unweighted.weight.size() == no_free.size());
CHECK(unweighted.weight[0].m == 1.0);
CHECK(unweighted.weight[0].d == 1.0);
}
// The model rugnux scored has to reach disk, or a user overlaying their input model on rugnux's maps
// is wrong by the whole rigid-body shift. Two things have to hold: the file carries what the input
// carried, and it carries the cell and space group the reflection files beside it are written in -
// which, once --model has adopted the model's enantiomorph, is neither the data's original label nor
// necessarily the input model's.
// A depositor's model and rugnux's own reduction routinely describe the SAME lattice on different
// axes - I-centred against C-centred, unique axis b against c, a cyclic permutation of an
// orthorhombic cell - and the space-group NUMBER is identical in every one of those, so nothing that
// compares numbers can notice. Re-fractionalizing straight across such a pair scrambles the model,
// and the rigid body cannot undo it: six parameters about a centroid are not a change of basis. The
// check is closed - the data are the model's own structure factors - so a model that reaches the
// data's description of the lattice comes out at R near zero, and one that does not comes out near
// the 0.6 an unrelated structure gives.
TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidation]") {
Logger logger("ModelValidation_FindsTheDatasDescriptionOfTheLattice");
// The data: the cluster's own |F|^2 in its own frame, 30 x 34 x 38, P 21 21 21.
const auto data_model = WriteTemp("cob_data_test.pdb", ClusterPdb().c_str());
auto obs = ModelReferenceIntensities(data_model, {}, {}, 2.5, logger);
REQUIRE(obs.size() > 500);
for (size_t i = 0; i < obs.size(); i++) {
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
obs[i].rfree_flag = (i % 20) == 0;
}
const UnitCell data_cell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *data_sg = gemmi::find_spacegroup_by_name("P 21 21 21");
const std::string prefix = (std::filesystem::temp_directory_path() / "cob_test").string();
// The same structure as a deposition on permuted axes, a' = b, b' = c, c' = a. The cell is
// orthogonal, so the Cartesian coordinates permute with the axes.
std::string permuted = "CRYST1 34.000 38.000 30.000 90.00 90.00 90.00 P 21 21 21 4\n";
{
std::istringstream in(ClusterPdb());
std::string line;
char buf[96];
while (std::getline(in, line)) {
if (line.compare(0, 4, "ATOM") != 0)
continue;
std::snprintf(buf, sizeof buf, "%s%8.3f%8.3f%8.3f%s",
line.substr(0, 30).c_str(), std::stod(line.substr(38, 8)),
std::stod(line.substr(46, 8)), std::stod(line.substr(30, 8)),
line.substr(54).c_str());
permuted += buf;
permuted += "\n";
}
permuted += "END\n";
}
const auto other_setting = WriteTemp("cob_model_test.pdb", permuted.c_str());
const auto moved = ValidateAgainstModel(obs, data_cell, other_setting, prefix, logger, data_sg,
/*probe_indexing_ambiguity=*/false, 1, 1.0);
REQUIRE(moved.ok);
logger.Info("Change of basis test: took {} out of {}, R-work {:.4f}",
moved.change_of_basis_op.triplet(), moved.setting_as_read, moved.r_work);
CHECK_FALSE(moved.change_of_basis_op == gemmi::Op::identity());
CHECK(moved.r_work < 0.15);
// The isomorphous case, which is what nearly every run is: the model already describes the
// lattice the way the data do, the only operators that map its cell are its own symmetry, and
// nothing is scored or moved.
const auto same = ValidateAgainstModel(obs, data_cell, data_model, prefix, logger, data_sg,
false, 1, 1.0);
REQUIRE(same.ok);
CHECK(same.change_of_basis_op == gemmi::Op::identity());
CHECK(same.setting_as_read.empty());
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(data_model);
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");
const auto input = WriteTemp("write_model_test_input.pdb", kPdbRich);
gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
// The model in the data's own group, which is the ordinary case: the frame it is written in is
// then wholly the caller's. The cell it arrives in (40/50/60) is not the data's and is replaced.
st.spacegroup_hm = "P 43 21 2";
// The tetragonal lysozyme cell and one enantiomorph of its group, standing in for what
// AdoptModelFrame settled.
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
REQUIRE(sg != nullptr);
WritePlacedModel(st, data_cell, *sg, "write_model_test", logger);
const std::string written = "write_model_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
CHECK(back.cell.a == Catch::Approx(79.0));
CHECK(back.cell.c == Catch::Approx(38.0));
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 96);
REQUIRE(back.models.size() == 1);
REQUIRE(back.models[0].chains.size() == 1);
const gemmi::Chain &chain = back.models[0].chains[0];
REQUIRE(chain.residues.size() == 2); // GLY and the water
REQUIRE(chain.residues[0].atoms.size() == 2);
REQUIRE(chain.residues[1].atoms.size() == 1);
// The coordinates are written as they stand - WritePlacedModel changes the label, not the atoms.
const gemmi::Atom &ca = chain.residues[0].atoms[0];
const gemmi::Atom &cb = chain.residues[0].atoms[1];
CHECK(ca.pos.x == Catch::Approx(10.0));
CHECK(ca.pos.z == Catch::Approx(14.0));
CHECK(cb.pos.y == Catch::Approx(14.0));
// B-factors, occupancies and the anisotropic U survive the round trip.
CHECK(ca.b_iso == Catch::Approx(20.0));
CHECK(cb.b_iso == Catch::Approx(25.0));
CHECK(cb.occ == Catch::Approx(0.60));
CHECK(ca.aniso.nonzero());
CHECK(ca.aniso.u11 == Catch::Approx(0.1000).margin(1e-4));
CHECK(ca.aniso.u12 == Catch::Approx(0.0100).margin(1e-4));
CHECK(chain.residues[1].name == "HOH");
std::filesystem::remove(input);
std::filesystem::remove(written);
}
// A model cannot be labelled with a group its own coordinates do not obey. Where the data were merged
// in a supergroup of the model's - an over-merge across a pseudo-symmetry operation, or the hand the
// model's fit did not earn - that operation would generate atoms the model does not contain, and a
// refinement program acts on it silently. The model keeps its own group and the caller is warned.
TEST_CASE("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey", "[ModelValidation]") {
Logger logger("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey");
const auto input = WriteTemp("write_model_sg_test_input.pdb", kPdbRich); // P 1
const gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
REQUIRE(sg != nullptr);
WritePlacedModel(st, data_cell, *sg, "write_model_sg_test", logger);
const std::string written = "write_model_sg_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
// The cell is still the data's - the coordinates do sit in it - but the group is the model's own.
CHECK(back.cell.a == Catch::Approx(79.0));
CHECK(back.cell.c == Catch::Approx(38.0));
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 1);
std::filesystem::remove(input);
std::filesystem::remove(written);
std::filesystem::remove("write_model_sg_test_model.pdb");
}
// The other direction is not a mismatch at all: a model whose coordinates obey more symmetry than the
// data were merged in obeys the data's group too - the file simply holds more than one asymmetric unit
// of it, which is what a model in P1 always does. The reflection files' group is kept, as before.
TEST_CASE("WriteModel_KeepsTheDataGroupWhenItIsASubgroup", "[ModelValidation]") {
Logger logger("WriteModel_KeepsTheDataGroupWhenItIsASubgroup");
const auto input = WriteTemp("write_model_sub_test_input.pdb", kPdbRich);
gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
st.spacegroup_hm = "P 43 21 2";
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *p1 = gemmi::find_spacegroup_by_name("P 1");
REQUIRE(p1 != nullptr);
WritePlacedModel(st, data_cell, *p1, "write_model_sub_test", logger);
const std::string written = "write_model_sub_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 1);
std::filesystem::remove(input);
std::filesystem::remove(written);
std::filesystem::remove("write_model_sub_test_model.pdb");
}
// CC(model, data) has to follow where the signal actually is, or it cannot support the one-sided
// claim it exists for. The check is closed: the "observed" intensities are the model's own with
// Gaussian noise added, and how much noise is chosen per shell - almost none in the first, some in
// the second, enough to bury the signal in the third - so the answer is known before the run.
TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") {
Logger logger("ModelValidation_CCModelFollowsTheSignalByShell");
const auto path = WriteTemp("cc_model_test.pdb", ClusterPdb().c_str());
auto obs = ModelReferenceIntensities(path, {}, {}, 2.5, logger);
REQUIRE(obs.size() > 1000);
// The shells the correlation is reported on, coarse to fine, and the noise each one gets as a
// multiple of the r.m.s. intensity of that shell. Nothing coarser than the first shell is kept:
// the reference intensities carry the bulk solvent at fixed constants while the validation fits
// its own, and below about 6 A that difference is a large part of |F| and would decorrelate a
// shell this test needs to be clean.
const std::vector<float> shells{5.0f, 3.2f, 2.5f};
const double noise[3] = {0.02, 1.0, 30.0};
auto shell_of = [&](float d) {
for (size_t s = 0; s < shells.size(); s++)
if (d > shells[s]) return s;
return shells.size();
};
std::erase_if(obs, [&](const MergedReflection &r) { return r.d > 6.0f || shell_of(r.d) >= shells.size(); });
REQUIRE(obs.size() > 500);
std::vector<double> sum_i2(shells.size(), 0.0);
std::vector<int> count(shells.size(), 0);
for (const auto &r : obs) {
sum_i2[shell_of(r.d)] += static_cast<double>(r.I) * r.I;
++count[shell_of(r.d)];
}
std::mt19937 rng(20260907);
std::normal_distribution<double> normal(0.0, 1.0);
for (size_t i = 0; i < obs.size(); i++) {
const size_t bin = shell_of(obs[i].d);
const double sd = noise[bin] * std::sqrt(sum_i2[bin] / count[bin]);
obs[i].I = static_cast<float>(obs[i].I + sd * normal(rng));
obs[i].sigma = static_cast<float>(std::max(1.0, sd));
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
obs[i].rfree_flag = (i % 20) == 0;
}
const std::string prefix = (std::filesystem::temp_directory_path() / "cc_model_test").string();
const auto result =
ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38,
.alpha = 90, .beta = 90, .gamma = 90},
path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"),
/*probe_indexing_ambiguity=*/false, 1, 1.0, shells);
REQUIRE(result.ok);
REQUIRE(result.cc_model_shells.size() == shells.size());
int n_total = 0;
for (size_t s = 0; s < shells.size(); s++) {
const auto &sh = result.cc_model_shells[s];
logger.Info("CC(model,data) {:.2f} A: {:.3f} on {} refl, {:+.1f} sigma",
sh.d_min, sh.cc, sh.n, sh.sigma);
CHECK(sh.d_min == shells[s]);
CHECK(sh.n > 20);
n_total += sh.n;
}
CHECK(n_total == result.cc_model_n);
// Essentially noiseless: the model is the data, so the correlation is high and hugely significant.
CHECK(result.cc_model_shells[0].cc > 0.9);
CHECK(result.cc_model_shells[0].sigma > 10.0);
// Noise at the shell's own r.m.s. still leaves plenty to see.
CHECK(result.cc_model_shells[1].cc > 0.25);
CHECK(result.cc_model_shells[1].sigma > 5.0);
// Buried: the shell must NOT come out significant, or the one-sided test would fire on noise.
CHECK(std::fabs(result.cc_model_shells[2].cc) < 0.15);
CHECK(std::fabs(result.cc_model_shells[2].sigma) < 4.0);
// No shells asked for, none reported: a run that did not measure it writes no key.
const auto no_shells =
ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38,
.alpha = 90, .beta = 90, .gamma = 90},
path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"),
false, 1, 1.0);
CHECK(no_shells.ok);
CHECK(no_shells.cc_model_shells.empty());
for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
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);
}
// 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.
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",
"CRYST1 40.000 50.000 60.000 90.00 100.00 90.00 C 1 2 1 4\n",
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 21 21 21 4\n",
"CRYST1 60.000 60.000 60.000 90.00 90.00 90.00 I 2 3 24\n",
"CRYST1 80.000 80.000 80.000 90.00 90.00 90.00 F 41 3 2 96\n",
};
for (const char *cryst : crysts) {
const auto path = WriteTemp("parallel_gridding_test.pdb", ClusterPdb(cryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
std::filesystem::remove(path);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
int i = 0;
for (gemmi::Chain &ch : st.models[0].chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
if (i++ % 3 == 0)
a.aniso = {0.30f, 0.25f, 0.20f, 0.02f, -0.01f, 0.03f};
const auto same = [](const gemmi::Grid<float> &a, const gemmi::Grid<float> &b) {
return a.data.size() == b.data.size() &&
std::memcmp(a.data.data(), b.data.data(), a.data.size() * sizeof(float)) == 0;
};
for (double d_min : {6.0, 3.5}) {
gemmi::DensityCalculator<gemmi::IT92<float>, float> gemmi_dc;
gemmi_dc.d_min = d_min;
gemmi_dc.rate = 1.5;
gemmi_dc.grid.unit_cell = st.cell;
gemmi_dc.grid.spacegroup = sg;
gemmi_dc.set_refmac_compatible_blur(st.models[0]);
gemmi_dc.put_model_density_on_grid(st.models[0]);
gemmi::Grid<float> gemmi_mask;
gemmi_mask.unit_cell = st.cell;
gemmi_mask.spacegroup = sg;
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]);
for (size_t nthreads : {1, 4}) {
const std::vector<size_t> leaders = OrbitLeaders(gemmi_mask, nthreads);
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.grid.unit_cell = st.cell;
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));
gemmi::Grid<float> mask;
mask.unit_cell = st.cell;
mask.spacegroup = sg;
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
CHECK(same(mask, gemmi_mask));
}
}
}
}
namespace {
// The five groups of ModelValidation_ParallelGriddingMatchesGemmi: none, a centring, screws, a
// cubic body centring and a cubic face centring with a diamond glide.
const char *kRigidBodyCrysts[] = {
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n",
"CRYST1 40.000 50.000 60.000 90.00 100.00 90.00 C 1 2 1 4\n",
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 21 21 21 4\n",
"CRYST1 60.000 60.000 60.000 90.00 90.00 90.00 I 2 3 24\n",
"CRYST1 80.000 80.000 80.000 90.00 90.00 90.00 F 41 3 2 96\n",
};
// ClusterPdb() in `cryst`, every third atom anisotropic.
gemmi::Structure AnisoCluster(const char *cryst) {
const auto path = WriteTemp("rigid_body_composition_test.pdb", ClusterPdb(cryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
std::filesystem::remove(path);
st.setup_cell_images();
int i = 0;
for (gemmi::Chain &ch : st.models[0].chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
if (i++ % 3 == 0)
a.aniso = {0.30f, 0.25f, 0.20f, 0.02f, -0.01f, 0.03f};
return st;
}
// The zone's density calculator, as the rigid body sets it up for `model`.
gemmi::DensityCalculator<gemmi::IT92<float>, float> ZoneDensity(const gemmi::Structure &st,
const gemmi::Model &model, double d_min) {
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.grid.unit_cell = st.cell;
dc.grid.spacegroup = st.find_spacegroup();
dc.set_refmac_compatible_blur(model);
return dc;
}
}
// The rigid body's Fcalc is composed from the transform of one copy of the model instead of being
// taken from the symmetrized grid. The two are the same sum rearranged, so they must agree to
// rounding, and must give a value to exactly the reflections prepare_asu_data() does - the
// systematic absences of the group and a reflection sitting exactly on d_min (the (10 0 0) of the
// 60 A cubic cell at 6 A) are left out by both.
TEST_CASE("ModelValidation_RigidBodyP1FcalcMatchesSymmetrized", "[ModelValidation]") {
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::Model &model = st.models[0];
for (double d_min : {6.0, 3.5}) {
auto sym = ZoneDensity(st, model, d_min);
sym.put_model_density_on_grid(model);
gemmi::FPhiGrid<float> sym_f = MapToFPhi(sym.grid);
const auto ref = sym_f.prepare_asu_data(d_min, sym.blur, false, false, false);
// Every index of the ASU in the grid's box, whatever its resolution, absences and 000 too.
std::vector<gemmi::Miller> candidates;
for (const auto &hv : sym_f.prepare_asu_data(0, 0, true, true, false).v)
candidates.push_back(hv.hkl);
auto copy = ZoneDensity(st, model, d_min);
PutModelDensityOnGrid(copy, model, {}, 4);
const SymmetryComposition composition(copy.grid, d_min, candidates);
std::vector<std::complex<double>> f;
composition.Compose(MapToFPhi(copy.grid), copy.blur, f, nullptr, 4);
REQUIRE(composition.Hkl().size() == ref.v.size());
double mean = 0, worst = 0;
for (size_t i = 0; i < ref.v.size(); i++) {
REQUIRE(composition.Hkl()[i] == ref.v[i].hkl);
mean += std::abs(ref.v[i].value) / static_cast<double>(ref.v.size());
worst = std::max(worst, std::abs(f[i] - std::complex<double>(ref.v[i].value)));
}
INFO(cryst << " at " << d_min << " A: worst " << worst / mean << " of the mean |F|");
CHECK(worst <= 1e-4 * mean);
}
}
}
// The derivative of the composed Fcalc with respect to a translation of the body is a phase factor
// per term, exact - against a central difference of the regridded, recomposed Fcalc. Along a
// direction the origin is free in (all three in P1, b in C2) the amplitude does not change at all.
TEST_CASE("ModelValidation_RigidBodyTranslationDerivativeIsExact", "[ModelValidation]") {
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const double d_min = 3.5;
const gemmi::DensityCalculator<gemmi::IT92<float>, float> zone = ZoneDensity(st, st.models[0], d_min);
gemmi::Grid<float> grid;
grid.unit_cell = zone.grid.unit_cell;
grid.spacegroup = zone.grid.spacegroup;
grid.set_size_from_spacing(zone.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
std::vector<gemmi::Miller> candidates;
const gemmi::ReciprocalAsu asu(grid.spacegroup);
for (int h = -25; h <= 25; h++)
for (int k = -25; k <= 25; k++)
for (int l = -25; l <= 25; l++)
if (asu.is_in({{h, k, l}}))
candidates.push_back({{h, k, l}});
const SymmetryComposition composition(grid, d_min, candidates);
REQUIRE(composition.Hkl().size() > 200);
auto fcalc = [&](const gemmi::Vec3 &shift, std::vector<std::array<std::complex<double>, 3>> *df_dt) {
gemmi::Model model = st.models[0];
for (gemmi::Chain &ch : model.chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
a.pos += gemmi::Position(shift);
auto dc = ZoneDensity(st, model, d_min);
PutModelDensityOnGrid(dc, model, {}, 4);
std::vector<std::complex<double>> f;
composition.Compose(MapToFPhi(dc.grid), dc.blur, f, df_dt, 4);
return f;
};
std::vector<std::array<std::complex<double>, 3>> df_dt;
const std::vector<std::complex<double>> f = fcalc({}, &df_dt);
// A five-point difference, whose truncation error goes as eps^4: the grid is float, and at the
// 1e-3 A of a plain central difference its rounding alone reads as 3e-3 of the derivative. What
// is left, about 5e-4, is the sampling - an atom gridded after a move is the moved atom only to
// the accuracy of the grid - which the exact derivative does not have and the difference does.
const double eps = 0.05;
const bool p1 = st.find_spacegroup()->number == 1, c2 = st.find_spacegroup()->number == 5;
for (int k = 0; k < 3; k++) {
gemmi::Vec3 e;
e.at(k) = eps;
const auto p_1 = fcalc(e, nullptr), m_1 = fcalc(-e, nullptr);
const auto p_2 = fcalc(2 * e, nullptr), m_2 = fcalc(-2 * e, nullptr);
double diff = 0, norm = 0, amplitude = 0;
for (size_t m = 0; m < f.size(); m++) {
const std::complex<double> numeric = (8.0 * (p_1[m] - m_1[m]) - (p_2[m] - m_2[m])) / (12 * eps);
diff += std::norm(numeric - df_dt[m][k]);
norm += std::norm(df_dt[m][k]);
amplitude += gemmi::sq(std::real(std::conj(f[m]) * df_dt[m][k]) / std::abs(f[m]));
}
INFO(cryst << " axis " << k << ": relative error " << std::sqrt(diff / norm)
<< ", amplitude part " << std::sqrt(amplitude / norm));
CHECK(std::sqrt(diff / norm) <= 1e-3);
if (p1 || (c2 && k == 1))
CHECK(std::sqrt(amplitude / norm) <= 1e-5);
}
}
}
namespace {
// The rigid body's Jacobian at q against a central difference of its own residuals - which
// re-fit the scale at every evaluation - with the bulk-solvent mask held, as the Jacobian holds
// it. Per column (those in `columns`): the cosine between the two and the ratio of their norms.
// Also logged, not checked: how far the difference moves once the mask is let move with the body,
// which is what holding it costs.
// The model's own amplitudes to `d_min`: "observed" data whose minimum is where the model is.
gemmi::AsuData<gemmi::ValueSigma<float>> OwnAmplitudes(const char *cryst, const gemmi::Structure &st,
double d_min, Logger &logger) {
const auto path = WriteTemp("rigid_body_own_amplitudes.pdb", ClusterPdb(cryst).c_str());
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
std::filesystem::remove(path);
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = st.find_spacegroup();
for (const auto &r : ref)
if (st.cell.calculate_d({{r.h, r.k, r.l}}) >= d_min)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
return fobs;
}
// The rigid body's target on the CPU, or on an engine of `gpu` where it is given.
std::unique_ptr<RigidBodyTargetBase> MakeTarget(gemmi::Model &model, const gemmi::Structure &st,
RigidBodyGPUPool *gpu) {
#ifdef JFJOCH_USE_CUDA
if (gpu != nullptr)
return std::make_unique<RigidBodyTargetGPU>(*gpu, model, st.cell, *st.find_spacegroup(), 4);
#endif
return std::make_unique<RigidBodyTarget>(model, st.cell, *st.find_spacegroup(), 4);
}
void CheckRigidBodyJacobian(const char *cryst, const double q0[6], const std::vector<int> &columns,
double min_cosine, double max_norm_error, bool gpu = false) {
Logger logger("CheckRigidBodyJacobian");
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
// "Observed" amplitudes: the model's own, as placed at q = 0.
const auto path = WriteTemp("rigid_body_jacobian_test.pdb", ClusterPdb(cryst).c_str());
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
std::filesystem::remove(path);
REQUIRE_FALSE(ref.empty());
for (double zone : {6.0, 3.5}) {
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
if (st.cell.calculate_d({{r.h, r.k, r.l}}) >= zone)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
gemmi::Model model = st.models[0];
RigidBodyGPUPool *pool = nullptr;
#ifdef JFJOCH_USE_CUDA
std::unique_ptr<RigidBodyGPUPool> engines;
if (gpu) {
engines = RigidBodyGPUPool::Create(model, st.cell, *sg, zone, fobs.v.size(), 1, logger);
REQUIRE(engines);
pool = engines.get();
}
#else
(void) gpu;
#endif
const std::unique_ptr<RigidBodyTargetBase> target_backend = MakeTarget(model, st, pool);
RigidBodyTargetBase &target = *target_backend;
target.SetZone(fobs, zone);
const size_t n = target.NumObservations();
std::vector<double> r(n), jacobian(n * 6);
REQUIRE(target.Residuals(q0, r.data()));
REQUIRE(target.Jacobian(q0, jacobian.data()));
auto difference = [&](int j, double h) {
double qp[6], qm[6];
std::copy(q0, q0 + 6, qp);
std::copy(q0, q0 + 6, qm);
qp[j] += h;
qm[j] -= h;
std::vector<double> rp(n), rm(n), d(n);
REQUIRE(target.Residuals(qp, rp.data()));
REQUIRE(target.Residuals(qm, rm.data()));
for (size_t i = 0; i < n; i++)
d[i] = (rp[i] - rm[i]) / (2 * h);
return d;
};
std::vector<std::vector<double>> held(6), moving(6);
target.hold_mask = true;
REQUIRE(target.Residuals(q0, r.data())); // the mask the difference holds is q0's
for (int j : columns)
held[j] = difference(j, 0.005);
// The mask is binary on the grid, so with it moving the difference is taken at the step the
// Jacobian had when every column moved it (0.01 of the zone's resolution): at 0.005 A it
// mostly measures which grid points a few atoms happened to cross.
target.hold_mask = false;
for (int j : columns)
moving[j] = difference(j, 0.01 * zone);
double largest = 0;
for (int j : columns) {
double s = 0;
for (size_t i = 0; i < n; i++)
s += gemmi::sq(held[j][i]);
largest = std::max(largest, std::sqrt(s));
}
for (int j : columns) {
double dot = 0, ours = 0, theirs = 0, dot_moving = 0, norm_moving = 0;
for (size_t i = 0; i < n; i++) {
dot += jacobian[i * 6 + j] * held[j][i];
ours += gemmi::sq(jacobian[i * 6 + j]);
theirs += gemmi::sq(held[j][i]);
dot_moving += moving[j][i] * held[j][i];
norm_moving += gemmi::sq(moving[j][i]);
}
ours = std::sqrt(ours);
theirs = std::sqrt(theirs);
norm_moving = std::sqrt(norm_moving);
const double cosine = dot / (ours * theirs);
logger.Info("{} at {:.1f} A, column {}: cosine {:.5f}, norm {:.4f} of the difference's; "
"with the mask moving, cosine {:.4f} and norm {:.4f}",
std::string(cryst).substr(55, 11), zone, j, cosine, ours / theirs,
dot_moving / (norm_moving * theirs), norm_moving / theirs);
INFO(cryst << " at " << zone << " A, column " << j << ": cosine " << cosine
<< ", norm ratio " << ours / theirs);
if (theirs < 1e-3 * largest) {
// A direction the origin is free in: |F| does not change, so neither may the column.
CHECK(ours < 1e-2 * largest);
continue;
}
CHECK(cosine >= min_cosine);
CHECK(std::fabs(ours / theirs - 1) <= max_norm_error);
}
}
}
}
// The rotation columns at a placement already rotated (about 1.4 deg, so the step is taken from a
// rotated body and not from the model as read) with anisotropic atoms, whose U the placement does not
// turn with the body - the columns must differentiate exactly the function the residuals evaluate.
TEST_CASE("ModelValidation_RigidBodyRotationDerivativeMatchesDifferences", "[ModelValidation]") {
const double q0[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts)
CheckRigidBodyJacobian(cryst, q0, {0, 1, 2}, 0.99, 0.03);
}
// The whole Jacobian - analytic translation, forward-difference rotation, the scale re-fit folded in by
// projection - against the central difference of the residuals it describes. The projection is
// Kaufman's, which leaves out a term that grows with the residuals (the scale's derivatives depend on
// the placement too, weighted by how badly the model fits): with the body 0.07 A off it agrees with the
// difference to about 2e-3 in cosine, with it 0.35 A off to about 2e-2, and the bounds follow.
TEST_CASE("ModelValidation_RigidBodyJacobianMatchesNumericScaleRefit", "[ModelValidation]") {
const double close[6] = {0.0, 0.0, 0.0, 0.05, -0.04, 0.03};
const double far[6] = {0.0, 0.0, 0.0, 0.25, -0.20, 0.15};
for (const char *cryst : kRigidBodyCrysts) {
CheckRigidBodyJacobian(cryst, close, {0, 1, 2, 3, 4, 5}, 0.99, 0.02);
CheckRigidBodyJacobian(cryst, far, {0, 1, 2, 3, 4, 5}, 0.975, 0.03);
}
}
// The null checks where a replicate ENDED against the orientations equivalent to the model's, from the
// rotation the rigid body reports - which must be the rotation it applied: every atom's offset from the
// centroid after the refinement is that rotation of its offset before.
TEST_CASE("ModelValidation_RigidBodyReportsTheRotationItApplied", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyReportsTheRotationItApplied");
const auto path = WriteTemp("rigid_body_rotation_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);
std::filesystem::remove(path);
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();
// Turned 3 deg about z through the centroid, so there is a rotation to take back.
std::vector<gemmi::Position> turned = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : turned)
centre += p;
centre *= 1.0 / static_cast<double>(turned.size());
const double a = 3.0 * PI / 180.0;
const gemmi::Mat33 rz(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1);
for (gemmi::Position &p : turned)
p = gemmi::Position(rz.multiply(gemmi::Vec3(p) - centre) + centre);
SetModelPositions(st.models[0], turned);
const RigidBodyRefineResult result = RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
REQUIRE(result.angle_deg > 1.0);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
gemmi::Vec3 moved_centre;
for (const gemmi::Position &p : refined)
moved_centre += p;
moved_centre *= 1.0 / static_cast<double>(refined.size());
double worst = 0;
for (size_t i = 0; i < refined.size(); i++) {
const gemmi::Vec3 expected = result.rotation.multiply(gemmi::Vec3(turned[i]) - centre);
worst = std::max(worst, (gemmi::Vec3(refined[i]) - moved_centre - expected).length());
}
CHECK(worst < 1e-6);
const double trace = result.rotation[0][0] + result.rotation[1][1] + result.rotation[2][2];
CHECK(std::acos(std::clamp((trace - 1) / 2, -1.0, 1.0)) * 180 / PI ==
Catch::Approx(result.angle_deg).margin(1e-6));
}
#ifdef JFJOCH_USE_CUDA
// The GPU target is the CPU target's function, computed on the device: at the same placement the two
// give the same residuals and the same Jacobian to rounding (float distances on the device, cuFFT for
// FFTW), on every group of the composition tests, both zones, anisotropic atoms included.
TEST_CASE("RigidBodyGPU_MatchesCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0) {
SUCCEED("no GPU");
return;
}
Logger logger("RigidBodyGPU_MatchesCPU");
const double q0[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
for (double zone : {6.0, 3.5}) {
const auto fobs = OwnAmplitudes(cryst, st, zone, logger);
gemmi::Model cpu_model = st.models[0], gpu_model = st.models[0];
auto pool = RigidBodyGPUPool::Create(gpu_model, st.cell, *sg, zone, fobs.v.size(), 1, logger);
REQUIRE(pool);
RigidBodyTarget cpu(cpu_model, st.cell, *sg, 4);
RigidBodyTargetGPU gpu(*pool, gpu_model, st.cell, *sg, 4);
cpu.SetZone(fobs, zone);
gpu.SetZone(fobs, zone);
const size_t n = cpu.NumObservations();
REQUIRE(gpu.NumObservations() == n);
std::vector<double> rc(n), rg(n), jc(6 * n), jg(6 * n);
REQUIRE(cpu.Residuals(q0, rc.data()));
REQUIRE(gpu.Residuals(q0, rg.data()));
REQUIRE(cpu.Jacobian(q0, jc.data()));
REQUIRE(gpu.Jacobian(q0, jg.data()));
CHECK(gpu.unmatched == cpu.unmatched);
CHECK(gpu.k_sol == cpu.k_sol);
CHECK(gpu.b_sol == cpu.b_sol);
double worst = 0, rms = 0;
for (size_t i = 0; i < n; i++) {
worst = std::max(worst, std::fabs(rc[i] - rg[i]));
rms += rc[i] * rc[i] / static_cast<double>(n);
}
INFO(cryst << " at " << zone << " A: residuals differ by at most " << worst << ", rms residual "
<< std::sqrt(rms));
// Most cases agree to a few 1e-6 of <Fobs>. The bound is set by gemmi's scale fit instead: its
// Levenberg-Marquardt stops at a relative change of 1e-5, so a near-tie in its accept or stop
// decision - which a 1e-5 change of Fcalc can flip, on the CPU alone - moves the scale by up
// to about 1e-4 of |F|.
CHECK(worst <= 2e-4);
for (int j = 0; j < 6; j++) {
double diff = 0, norm = 0;
for (size_t i = 0; i < n; i++) {
diff += gemmi::sq(jc[6 * i + j] - jg[6 * i + j]);
norm += gemmi::sq(jc[6 * i + j]);
}
INFO(cryst << " at " << zone << " A, column " << j << ": relative difference "
<< std::sqrt(diff / norm));
CHECK(std::sqrt(diff) <= 2e-3 * std::sqrt(norm) + 1e-9);
}
}
}
}
// The whole Jacobian on the GPU against the central difference of the GPU's own residuals, with the
// bounds of ModelValidation_RigidBodyJacobianMatchesNumericScaleRefit.
TEST_CASE("RigidBodyGPU_JacobianMatchesNumericScaleRefit", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0) {
SUCCEED("no GPU");
return;
}
const double close[6] = {0.0, 0.0, 0.0, 0.05, -0.04, 0.03};
const double rotated[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts) {
CheckRigidBodyJacobian(cryst, close, {0, 1, 2, 3, 4, 5}, 0.99, 0.02, true);
CheckRigidBodyJacobian(cryst, rotated, {0, 1, 2}, 0.99, 0.03, true);
}
}
namespace {
// A whole placement on the GPU from the model displaced by 0.54 A and turned 2 deg, against the
// model's own amplitudes; the model is left where the fit put it.
RigidBodyRefineResult DisplacedFit(const char *cryst, gemmi::Structure &st, RigidBodyGPUPool *pool,
Logger &logger) {
const gemmi::SpaceGroup *sg = st.find_spacegroup();
const auto fobs = OwnAmplitudes(cryst, st, 3.0, logger);
std::vector<gemmi::Position> moved = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : moved)
centre += p;
centre *= 1.0 / static_cast<double>(moved.size());
const double a = 2.0 * PI / 180.0;
const gemmi::Mat33 rz(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1);
for (gemmi::Position &p : moved)
p = gemmi::Position(rz.multiply(gemmi::Vec3(p) - centre) + centre + gemmi::Vec3(0.40, -0.30, 0.20));
SetModelPositions(st.models[0], moved);
return RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger, 4, pool);
}
}
// The GPU fit walks where the CPU fit walks: from the same displaced start, to the same placement
// within a few thousandths of an angstrom.
TEST_CASE("RigidBodyGPU_FitAgreesWithCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0) {
SUCCEED("no GPU");
return;
}
Logger logger("RigidBodyGPU_FitAgreesWithCPU");
for (const char *cryst : {kCryst, kPolarCryst, kRigidBodyCrysts[3]}) {
gemmi::Structure cpu_st = AnisoCluster(cryst), gpu_st = AnisoCluster(cryst);
auto pool = RigidBodyGPUPool::Create(gpu_st.models[0], gpu_st.cell, *gpu_st.find_spacegroup(), 3.0, 100000, 1,
logger);
REQUIRE(pool);
const RigidBodyRefineResult cpu = DisplacedFit(cryst, cpu_st, nullptr, logger);
const RigidBodyRefineResult gpu = DisplacedFit(cryst, gpu_st, pool.get(), logger);
CHECK(gpu.converged == cpu.converged);
const std::vector<gemmi::Position> pc = ModelPositions(cpu_st.models[0]), pg = ModelPositions(gpu_st.models[0]);
double rmsd = 0;
for (size_t i = 0; i < pc.size(); i++)
rmsd += pc[i].dist_sq(pg[i]) / static_cast<double>(pc.size());
INFO(cryst << ": CPU " << cpu.angle_deg << " deg " << cpu.shift_A << " A, GPU " << gpu.angle_deg << " deg "
<< gpu.shift_A << " A, " << std::sqrt(rmsd) << " A apart");
CHECK(std::sqrt(rmsd) < 2e-3);
}
}
// Deterministic: the same fit twice, and on a pool of one engine and of four, gives the same placement
// bit for bit.
TEST_CASE("RigidBodyGPU_Deterministic", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0) {
SUCCEED("no GPU");
return;
}
Logger logger("RigidBodyGPU_Deterministic");
const char *cryst = kRigidBodyCrysts[4];
std::vector<std::vector<gemmi::Position>> placed;
for (size_t engines : {1, 1, 4}) {
gemmi::Structure st = AnisoCluster(cryst);
auto pool = RigidBodyGPUPool::Create(st.models[0], st.cell, *st.find_spacegroup(), 3.0, 100000, engines, logger);
REQUIRE(pool);
const RigidBodyRefineResult r = DisplacedFit(cryst, st, pool.get(), logger);
CHECK(r.evaluations > 0);
placed.push_back(ModelPositions(st.models[0]));
}
for (size_t k = 1; k < placed.size(); k++)
for (size_t i = 0; i < placed[0].size(); i++) {
CHECK(placed[k][i].x == placed[0][i].x);
CHECK(placed[k][i].y == placed[0][i].y);
CHECK(placed[k][i].z == placed[0][i].z);
}
}
#endif
Reference in New Issue View Git Blame Copy Permalink