leonarski_fandClaude Opus 5.5
dbc550d29f
On a large-cell merohedral serial run the reference-free resolution left about a fifth of the crystals in the wrong indexing: the merge read twin-pair CC 0.56 against 0.15-0.31 for the deposited data, and a twin refinement of it found a twin fraction of 0.21. Started from the reference-resolved assignment, the resolver itself moved a fifth of the crystals away, and the different starts ended in different mixed states - it was the comparison, not the start. Two defects in how the crystals were put side by side: - Intensities were divided by the run's mean in their resolution shell, and when any shell's mean was zero or below - which every shell past the resolution the run reaches can be, since integration predicts to the detector edge - all shells fell back to the raw intensities. The comparison was then dominated by the fall-off with resolution, which both indexings share. - Intensities follow the Wilson distribution and a stills partial carries a multiplicative error as large as itself; the long tail let a few very strong reflections dominate the product of two crystals. Groups of strong crystals then held one another in the wrong indexing against the weaker majority: single-crystal moves cannot undo such a group, and which groups formed depended on the order of the passes. The crystals are now compared on the square root of the intensity (sign kept), put on each shell's mean and spread rather than divided by its mean. Every start then ends in the same assignment. On the large-cell run 98.8% of the crystals a reference decides clearly agree with the reference (85.7% before); the merge's twin-pair CC is 0.22, its twin refinement finds a twin fraction of 0.025 and calls the data untwinned, R_split 7.6% -> 6.4%, and the battery R-free 0.258 -> 0.247 (deposited data 0.244, published 0.249). A trigonal serial set with the same kind of ambiguity: twin-pair CC 0.19 -> 0.09, R-free 0.196 -> 0.191. Sets whose cells decide the ambiguity are unchanged. Deterministic and thread-count independent as before (same sequential passes; only the per-reflection values changed). New unit test: 2000 weak stills with heavy-tailed intensities, a resolution fall-off and correlated twin mates; the old comparison put 9.9% of them in the wrong indexing, now under 3%. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
626 lines
30 KiB
C++
626 lines
30 KiB
C++
// 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 <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <random>
|
|
#include <unordered_map>
|
|
#include <unordered_set>
|
|
#include <vector>
|
|
|
|
#include "../image_analysis/scale_merge/CrystalSetting.h"
|
|
#include "../image_analysis/scale_merge/HKLKey.h"
|
|
#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
|
|
#include "../image_analysis/scale_merge/ReferenceFreeReindex.h"
|
|
#include "../image_analysis/scale_merge/SublatticeRelabel.h"
|
|
|
|
namespace {
|
|
const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); }
|
|
|
|
UnitCell Tetragonal() { return UnitCell{78.0f, 78.0f, 37.0f, 90.0f, 90.0f, 90.0f}; }
|
|
|
|
// A reference set: one distinct intensity per Laue-ASU reflection, so that reflections related by
|
|
// a twin law (which are NOT Laue-equivalent in P4) carry different intensities.
|
|
std::vector<MergedReflection> DistinctReference(int space_group_number) {
|
|
const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(space_group_number));
|
|
std::vector<MergedReflection> v;
|
|
std::unordered_set<uint64_t> seen;
|
|
for (int h = -8; h <= 8; ++h)
|
|
for (int k = -8; k <= 8; ++k)
|
|
for (int l = 0; l <= 8; ++l) {
|
|
if (h == 0 && k == 0 && l == 0) continue;
|
|
const uint64_t kk = key(h, k, l).pack();
|
|
if (!seen.insert(kk).second) continue; // one row per ASU reflection
|
|
MergedReflection r;
|
|
r.h = h; r.k = k; r.l = l;
|
|
r.I = static_cast<float>(100 + kk % 9973); // distinct, spread
|
|
r.sigma = 1.0f;
|
|
r.d = 60.0f / (1 + h * h + k * k + l * l);
|
|
v.push_back(r);
|
|
}
|
|
return v;
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Reindex: holohedral crystal has no indexing ambiguity", "[reindex]") {
|
|
// P4(3)2(1)2 (422, holohedral for the tetragonal lattice) -> no twin laws.
|
|
CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(96)).empty());
|
|
// P422 likewise.
|
|
CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(89)).empty());
|
|
}
|
|
|
|
TEST_CASE("Reindex: merohedral crystal exposes the ambiguity operators", "[reindex]") {
|
|
// P4 (point group 4) in a tetragonal lattice (422) -> a non-trivial reindexing coset.
|
|
CHECK_FALSE(ReindexAmbiguityOperators(Tetragonal(), SG(75)).empty());
|
|
}
|
|
|
|
TEST_CASE("Reindex: reference agreement recovers a misindexed dataset", "[reindex]") {
|
|
const int sg = 75; // P4
|
|
const auto reference = DistinctReference(sg);
|
|
const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg));
|
|
REQUIRE_FALSE(laws.empty());
|
|
|
|
// Deliberately mis-index the data by one twin law.
|
|
const auto data = ReindexReflections(reference, laws.front());
|
|
|
|
const auto choice = ChooseReindex(
|
|
data, Tetragonal(), SG(sg),
|
|
[&](const std::vector<MergedReflection> &m) { return ReferenceIntensityCC(m, reference, SG(sg)); });
|
|
|
|
CHECK_FALSE(choice.is_identity); // a reindex was needed
|
|
CHECK(choice.score > 0.99); // the winner realigns with the reference
|
|
CHECK(choice.identity_score < 0.9); // leaving it mis-indexed correlates poorly
|
|
CHECK(choice.n_candidates >= 2);
|
|
}
|
|
|
|
TEST_CASE("Reindex: a correctly indexed dataset keeps identity", "[reindex]") {
|
|
const int sg = 75;
|
|
const auto reference = DistinctReference(sg);
|
|
const auto choice = ChooseReindex(
|
|
reference, Tetragonal(), SG(sg),
|
|
[&](const std::vector<MergedReflection> &m) { return ReferenceIntensityCC(m, reference, SG(sg)); });
|
|
CHECK(choice.is_identity);
|
|
CHECK(choice.score > 0.99);
|
|
}
|
|
|
|
TEST_CASE("Reindex into the ASU: a twin law permutes the reflections without losing any", "[reindex]") {
|
|
const int sg = 75; // P4
|
|
const auto reference = DistinctReference(sg);
|
|
const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg));
|
|
REQUIRE_FALSE(laws.empty());
|
|
|
|
// Mis-index by a twin law, then reindex back into the ASU: the labels must land where an export
|
|
// needs them, and the reflection the label carries must be the one the reference has there.
|
|
const auto misindexed = ReindexReflections(reference, laws.front());
|
|
const auto fixed = ReindexMergedIntoAsu(misindexed, laws.front(), SG(sg), /*merge_friedel=*/true);
|
|
REQUIRE(fixed.size() == reference.size());
|
|
|
|
const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(sg));
|
|
std::unordered_map<uint64_t, float> ref_by_key;
|
|
for (const auto &r : reference)
|
|
ref_by_key[key(r).pack()] = r.I;
|
|
|
|
std::unordered_set<uint64_t> seen;
|
|
for (const auto &r : fixed) {
|
|
const HKLKey k = key(r);
|
|
CHECK(k.h == r.h); // already at its own ASU index
|
|
CHECK(k.k == r.k);
|
|
CHECK(k.l == r.l);
|
|
CHECK(seen.insert(k.pack()).second); // no two reflections land on one label
|
|
const auto it = ref_by_key.find(k.pack());
|
|
REQUIRE(it != ref_by_key.end());
|
|
CHECK(it->second == r.I);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Reindex into the ASU: the change of hand swaps the Bijvoet halves", "[reindex]") {
|
|
// The change of hand is the inversion, so it leaves the Laue-ASU label alone and moves the
|
|
// anomalous signal instead - which is the whole of what adopting a model's hand does to intensities.
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(96); // P4(3)2(1)2
|
|
REQUIRE(sg != nullptr);
|
|
const HKLKeyGenerator key(true, *sg);
|
|
const HKLKey asu = key(3, 1, 2);
|
|
|
|
MergedReflection r;
|
|
r.h = asu.h; r.k = asu.k; r.l = asu.l;
|
|
r.I = 100.0f; r.sigma = 1.0f; r.d = 10.0f;
|
|
r.I_plus = 110.0f; r.sigma_plus = 2.0f; r.I_minus = 90.0f; r.sigma_minus = 3.0f;
|
|
r.F_plus = 10.5f; r.F_minus = 9.5f;
|
|
|
|
const auto merged = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/true);
|
|
REQUIRE(merged.size() == 1);
|
|
CHECK(merged[0].h == asu.h);
|
|
CHECK(merged[0].k == asu.k);
|
|
CHECK(merged[0].l == asu.l);
|
|
CHECK(merged[0].I == 100.0f);
|
|
CHECK(merged[0].I_plus == 90.0f);
|
|
CHECK(merged[0].I_minus == 110.0f);
|
|
CHECK(merged[0].sigma_plus == 3.0f);
|
|
CHECK(merged[0].sigma_minus == 2.0f);
|
|
CHECK(merged[0].F_plus == 9.5f);
|
|
CHECK(merged[0].F_minus == 10.5f);
|
|
|
|
// With the mates kept apart the merge stores the minus hand at -hkl, so the row moves there.
|
|
const auto anom = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/false);
|
|
REQUIRE(anom.size() == 1);
|
|
CHECK(anom[0].h == -asu.h);
|
|
CHECK(anom[0].k == -asu.k);
|
|
CHECK(anom[0].l == -asu.l);
|
|
CHECK(anom[0].I == 100.0f);
|
|
}
|
|
|
|
TEST_CASE("Reindex into the ASU: both mates of an anomalous pair follow the same hand", "[reindex]") {
|
|
// With the mates kept apart, the merge stores the plus mate at +hkl_asu and the minus mate at
|
|
// -hkl_asu, while I_plus/I_minus are attached in the plus convention on BOTH rows alike. An
|
|
// alternative-indexing operator is rotation-type, so op(-x) == -op(x): whatever it does to one
|
|
// mate's hand it does to the other's, and the two rows must still agree afterwards about which
|
|
// intensity is which. Deciding the swap from the hand the new index lands on instead of from the
|
|
// CHANGE of hand moves exactly one of the two, and the pair ends up self-contradictory.
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(75); // P4, Laue 4/m
|
|
REQUIRE(sg != nullptr);
|
|
const HKLKeyGenerator key_gen(/*merge_friedel=*/false, *sg);
|
|
const HKLKey asu = key_gen(3, 1, 2);
|
|
REQUIRE(asu.plus);
|
|
|
|
// The alternative indexing of a P4 lattice: 4/mmm holds it, 4/m does not.
|
|
const gemmi::Op op = gemmi::parse_triplet("k,h,-l");
|
|
|
|
MergedReflection plus;
|
|
plus.h = asu.h; plus.k = asu.k; plus.l = asu.l;
|
|
plus.I = 100.0f; plus.sigma = 1.0f; plus.d = 10.0f;
|
|
plus.I_plus = 110.0f; plus.sigma_plus = 2.0f; plus.I_minus = 90.0f; plus.sigma_minus = 3.0f;
|
|
plus.F_plus = 10.5f; plus.F_minus = 9.5f;
|
|
|
|
MergedReflection minus = plus; // same anomalous split, stored at -hkl_asu
|
|
minus.h = -asu.h; minus.k = -asu.k; minus.l = -asu.l;
|
|
|
|
const auto out = ReindexMergedIntoAsu({plus, minus}, op, SG(75), /*merge_friedel=*/false);
|
|
REQUIRE(out.size() == 2);
|
|
|
|
// The two rows are still one pair: same Laue label up to the Friedel sign, opposite hands.
|
|
CHECK(out[0].h == -out[1].h);
|
|
CHECK(out[0].k == -out[1].k);
|
|
CHECK(out[0].l == -out[1].l);
|
|
|
|
// ... and they still tell the same story about which hand carries which intensity.
|
|
CHECK(out[0].I_plus == out[1].I_plus);
|
|
CHECK(out[0].I_minus == out[1].I_minus);
|
|
CHECK(out[0].sigma_plus == out[1].sigma_plus);
|
|
CHECK(out[0].sigma_minus == out[1].sigma_minus);
|
|
CHECK(out[0].F_plus == out[1].F_plus);
|
|
CHECK(out[0].F_minus == out[1].F_minus);
|
|
|
|
// And the swap did have to happen: "k,h,-l" takes (3,1,2) to (1,3,-2), whose Laue class has its
|
|
// ASU representative at l > 0, so the row that was the plus mate of its pair is the minus mate of
|
|
// the new one. The intensity now standing at the plus index is the one that was at -hkl before.
|
|
CHECK(out[0].I_plus == 90.0f);
|
|
CHECK(out[0].I_minus == 110.0f);
|
|
}
|
|
|
|
TEST_CASE("Reindex: every ambiguity operator is a proper rotation", "[reindex]") {
|
|
// A relabelling by a lattice rotation keeps the hand; one with determinant -1 would swap every
|
|
// Bijvoet pair. The candidates come from the lattice symmetry without the inversion, so none may
|
|
// be improper, in any crystal system that has an ambiguity.
|
|
const std::vector<std::pair<UnitCell, int>> cases{
|
|
{Tetragonal(), 75}, // P4
|
|
{UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 143}, // P3
|
|
{UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 149}, // P312
|
|
{UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 168}, // P6
|
|
{UnitCell{77.0f, 77.0f, 77.0f, 90.0f, 90.0f, 90.0f}, 197}, // I23
|
|
{UnitCell{77.0f, 77.0f, 77.0f, 90.0f, 90.0f, 90.0f}, 199}, // I2(1)3
|
|
};
|
|
for (const auto &[cell, sg] : cases) {
|
|
const auto ops = ReindexAmbiguityOperators(cell, SG(sg));
|
|
CHECK_FALSE(ops.empty());
|
|
for (const auto &op : ops)
|
|
CHECK(op.det_rot() == gemmi::Op::DEN * gemmi::Op::DEN * gemmi::Op::DEN);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Reindex: a near-symmetry law moves the cell with the indices", "[reindex]") {
|
|
// Pseudo-merohedry: b within 2% of c in a primitive orthorhombic cell, and beta a degree off 90 in
|
|
// a monoclinic one. Each admits an indexing law that is not a symmetry of the cell, so relabelling
|
|
// the reflections by it permutes the axes (or flips beta), and the cell they are written with has
|
|
// to follow - every reflection then keeps its own d-spacing.
|
|
const std::vector<std::pair<UnitCell, int>> cases{
|
|
{UnitCell{50.0f, 60.0f, 61.2f, 90.0f, 90.0f, 90.0f}, 16}, // P222, b ~ c
|
|
{UnitCell{70.0f, 79.0f, 83.0f, 90.0f, 91.1f, 90.0f}, 3}, // P2, beta ~ 90
|
|
};
|
|
for (const auto &[cell, sg] : cases) {
|
|
const auto laws = ReindexAmbiguityOperators(cell, SG(sg));
|
|
REQUIRE_FALSE(laws.empty());
|
|
bool moved = false;
|
|
for (const gemmi::Op &law : laws) {
|
|
const UnitCell written = CellInBasis(cell, HklOperator(law));
|
|
moved = moved || std::abs(written.b - cell.b) > 0.5f || std::abs(written.beta - cell.beta) > 0.5f;
|
|
const gemmi::UnitCell before = cell, after = written;
|
|
for (const gemmi::Op::Miller hkl : {gemmi::Op::Miller{1, 2, 3}, gemmi::Op::Miller{-4, 0, 7},
|
|
gemmi::Op::Miller{5, -3, 1}})
|
|
CHECK(after.calculate_d(law.apply_to_hkl(hkl)) == Catch::Approx(before.calculate_d(hkl)).epsilon(1e-5));
|
|
}
|
|
CHECK(moved);
|
|
}
|
|
}
|
|
|
|
namespace {
|
|
// A P2(1) cell with a near-symmetry: a+c is within 1% of a, so the lattice is metrically close to
|
|
// a C-centred orthorhombic one (obliquity ~0.7 deg) without being one.
|
|
UnitCell PseudoMonoclinic() { return UnitCell{64.0f, 29.0f, 36.0f, 90.0f, 107.0f, 90.0f}; }
|
|
|
|
Coord Recip(const CrystalLattice &latt, int h, int k, int l) {
|
|
return latt.Astar() * static_cast<float>(h) + latt.Bstar() * static_cast<float>(k)
|
|
+ latt.Cstar() * static_cast<float>(l);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Reindex: an ambiguity operator of a pseudo-symmetric monoclinic cell keeps every reflection's resolution",
|
|
"[reindex]") {
|
|
const UnitCell cell = PseudoMonoclinic();
|
|
const auto laws = ReindexAmbiguityOperators(cell, SG(4));
|
|
REQUIRE_FALSE(laws.empty());
|
|
const CrystalLattice latt(cell);
|
|
for (const auto &op : laws)
|
|
for (int h = -6; h <= 6; ++h)
|
|
for (int l = -6; l <= 6; ++l) {
|
|
if (h == 0 && l == 0) continue;
|
|
// As the resolver applies it: an operator that is a near-symmetry of the metric moves an
|
|
// index to one of nearly the same length, never to a finer one.
|
|
const auto m = op.apply_to_hkl({{h, 3, l}});
|
|
const float d0 = 1.0f / Recip(latt, h, 3, l).Length();
|
|
const float d1 = 1.0f / Recip(latt, m[0], m[1], m[2]).Length();
|
|
CHECK(std::fabs(d1 - d0) / d0 < 0.03f);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Reindex: the reindexed lattice predicts every spot where it was", "[reindex]") {
|
|
const UnitCell cell = PseudoMonoclinic();
|
|
const auto laws = ReindexAmbiguityOperators(cell, SG(4));
|
|
REQUIRE_FALSE(laws.empty());
|
|
const CrystalLattice latt(cell);
|
|
for (const auto &op : laws) {
|
|
const CrystalLattice moved = ReindexLattice(latt, op);
|
|
for (int h = -4; h <= 4; ++h)
|
|
for (int k = -4; k <= 4; ++k)
|
|
for (int l = -4; l <= 4; ++l) {
|
|
const auto m = op.apply_to_hkl({{h, k, l}});
|
|
const Coord before = Recip(latt, h, k, l);
|
|
const Coord after = Recip(moved, m[0], m[1], m[2]);
|
|
CHECK((after - before).Length() < 1e-5f);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Sublattice: a frame on the given cell in another setting is described on the given cell's axes",
|
|
"[reindex]") {
|
|
// The given C222(1)-style cell, and the same lattice with b and c swapped (and a flipped to keep the
|
|
// basis right-handed): not close axis by axis, but the same lattice.
|
|
const UnitCell given{40.0f, 180.0f, 143.0f, 90.0f, 90.0f, 90.0f};
|
|
const UnitCell frame{40.2f, 143.4f, 180.5f, 90.0f, 90.0f, 90.0f};
|
|
IndexingSettings settings;
|
|
REQUIRE_FALSE(frame.is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
|
|
const auto on_given = CellOnReference(frame, given, settings);
|
|
REQUIRE(on_given.has_value());
|
|
CHECK(on_given->is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
|
|
CHECK(on_given->b == Catch::Approx(180.5f).epsilon(0.001));
|
|
CHECK(on_given->c == Catch::Approx(143.4f).epsilon(0.001));
|
|
|
|
// An unrelated cell is not put on it.
|
|
CHECK_FALSE(CellOnReference(UnitCell{55.0f, 66.0f, 77.0f, 90.0f, 90.0f, 90.0f}, given, settings).has_value());
|
|
}
|
|
|
|
TEST_CASE("Reindex: a stills set indexed in random hands is made consistent from the data alone", "[reindex]") {
|
|
// P4 on a tetragonal lattice: two indexings, related by a twin law. Every crystal records a random
|
|
// part of the same synthetic structure with its own scale and noise, and half of them, at random,
|
|
// carry the other indexing.
|
|
const UnitCell cell = Tetragonal();
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(cell).SpaceGroupNumber(75);
|
|
const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1());
|
|
REQUIRE(laws.size() == 1);
|
|
|
|
const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1());
|
|
std::mt19937 rng(7);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> noise(0.0, 0.3);
|
|
std::unordered_map<uint64_t, double> structure;
|
|
std::vector<std::array<int, 3>> hkls;
|
|
for (int h = -12; h <= 12; ++h)
|
|
for (int k = -12; k <= 12; ++k)
|
|
for (int l = -6; l <= 6; ++l) {
|
|
if (h == 0 && k == 0 && l == 0)
|
|
continue;
|
|
hkls.push_back({h, k, l});
|
|
structure.try_emplace(key(h, k, l).pack(), wilson(rng));
|
|
}
|
|
|
|
const int n_crystals = 300;
|
|
std::vector<IntegrationOutcome> outcomes(n_crystals);
|
|
std::vector<bool> flipped(n_crystals);
|
|
for (int c = 0; c < n_crystals; ++c) {
|
|
flipped[c] = rng() % 2;
|
|
const double scale = 0.5 + (rng() % 1000) / 1000.0;
|
|
outcomes[c].latt = CrystalLattice(cell);
|
|
for (int i = 0; i < 150; ++i) {
|
|
const auto &t = hkls[rng() % hkls.size()];
|
|
const auto m = flipped[c] ? laws[0].apply_to_hkl({{t[0], t[1], t[2]}}) : gemmi::Op::Miller{{t[0], t[1], t[2]}};
|
|
Reflection r{};
|
|
r.h = m[0]; r.k = m[1]; r.l = m[2];
|
|
r.d = 1.0f / std::sqrt((t[0] * t[0] + t[1] * t[1]) / (cell.a * cell.a) + t[2] * t[2] / (cell.c * cell.c));
|
|
r.I = static_cast<float>(scale * structure.at(key(t[0], t[1], t[2]).pack()) * (1.0 + noise(rng)));
|
|
r.sigma = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
outcomes[c].reflections.push_back(r);
|
|
}
|
|
}
|
|
|
|
const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/true);
|
|
CHECK(result.crystals == n_crystals);
|
|
// The overall hand is free: either exactly the flipped crystals are reindexed, or exactly the others.
|
|
int agree = 0;
|
|
for (int c = 0; c < n_crystals; ++c)
|
|
agree += (result.choice[c] != 0) == flipped[c];
|
|
CHECK((agree == n_crystals || agree == 0));
|
|
const size_t n_flipped = std::count(flipped.begin(), flipped.end(), true);
|
|
CHECK(result.reindexed == std::min(n_flipped, n_crystals - n_flipped));
|
|
CHECK(result.reindexed <= n_crystals / 2); // the overall hand that moves the fewest crystals
|
|
CHECK(result.decided_fraction > 0.9);
|
|
|
|
// Done again on the consistent set, nothing moves.
|
|
const auto again = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/false);
|
|
CHECK(again.reindexed == 0);
|
|
}
|
|
|
|
TEST_CASE("Reindex: many weak stills with heavy-tailed intensities are made consistent from the data alone",
|
|
"[reindex]") {
|
|
// As above, but in the regime of a serial run: many crystals with few reflections each, a
|
|
// multiplicative error on every intensity as large as the intensity (partiality), a resolution
|
|
// fall-off, and twin-related reflections that often carry similar intensities. Compared on the
|
|
// intensities themselves, a few very strong reflections dominate and a tenth of the crystals end up
|
|
// in the wrong indexing.
|
|
const UnitCell cell = Tetragonal();
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(cell).SpaceGroupNumber(75);
|
|
const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1());
|
|
REQUIRE(laws.size() == 1);
|
|
|
|
const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1());
|
|
std::mt19937 rng(5);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::uniform_real_distribution<double> uniform(0.0, 1.0);
|
|
std::normal_distribution<double> normal(0.0, 1.0);
|
|
std::unordered_map<uint64_t, double> structure;
|
|
std::vector<std::array<int, 3>> hkls;
|
|
for (int h = -14; h <= 14; ++h)
|
|
for (int k = -14; k <= 14; ++k)
|
|
for (int l = -8; l <= 8; ++l)
|
|
if (h != 0 || k != 0 || l != 0) {
|
|
hkls.push_back({h, k, l});
|
|
structure.try_emplace(key(h, k, l).pack(), wilson(rng));
|
|
}
|
|
// Two twin pairs in five carry nearly the same intensity.
|
|
std::unordered_set<uint64_t> paired;
|
|
for (const auto &t : hkls) {
|
|
const auto m = laws[0].apply_to_hkl({{t[0], t[1], t[2]}});
|
|
const uint64_t a = key(t[0], t[1], t[2]).pack(), b = key(m[0], m[1], m[2]).pack();
|
|
if (a == b || !paired.insert(a).second || !paired.insert(b).second)
|
|
continue;
|
|
if (uniform(rng) < 0.4)
|
|
structure[b] = structure[a] * (0.7 + 0.6 * uniform(rng));
|
|
}
|
|
|
|
const int n_crystals = 2000;
|
|
std::vector<IntegrationOutcome> outcomes(n_crystals);
|
|
std::vector<bool> flipped(n_crystals);
|
|
for (int c = 0; c < n_crystals; ++c) {
|
|
flipped[c] = rng() % 2;
|
|
const double scale = std::exp(0.5 * normal(rng));
|
|
outcomes[c].latt = CrystalLattice(cell);
|
|
for (int i = 0; i < 60; ++i) {
|
|
const auto &t = hkls[rng() % hkls.size()];
|
|
const double s2 = (t[0] * t[0] + t[1] * t[1]) / (cell.a * cell.a) + t[2] * t[2] / (cell.c * cell.c);
|
|
const auto m = flipped[c] ? laws[0].apply_to_hkl({{t[0], t[1], t[2]}}) : gemmi::Op::Miller{{t[0], t[1], t[2]}};
|
|
Reflection r{};
|
|
r.h = m[0]; r.k = m[1]; r.l = m[2];
|
|
r.d = static_cast<float>(1.0 / std::sqrt(s2));
|
|
r.I = static_cast<float>(scale * structure.at(key(t[0], t[1], t[2]).pack()) * std::exp(-15.0 * s2)
|
|
* std::exp(normal(rng)) + 0.02 * normal(rng));
|
|
r.sigma = 0.02f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
outcomes[c].reflections.push_back(r);
|
|
}
|
|
}
|
|
|
|
const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/false);
|
|
int agree = 0;
|
|
for (int c = 0; c < n_crystals; ++c)
|
|
agree += (result.choice[c] != 0) == flipped[c];
|
|
CHECK(std::max(agree, n_crystals - agree) > 0.97 * n_crystals);
|
|
}
|
|
|
|
TEST_CASE("Reindex: a holohedral crystal gives the stills resolution nothing to do", "[reindex]") {
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(Tetragonal()).SpaceGroupNumber(89); // P422
|
|
std::vector<IntegrationOutcome> outcomes(3);
|
|
const auto laws = ReindexAmbiguityOperators(Tetragonal(), x.GetSpaceGroupOrP1());
|
|
const auto result = ResolveReindexReferenceFree(outcomes, x, Tetragonal(), laws, 1, true);
|
|
CHECK(laws.empty());
|
|
CHECK(result.reindexed == 0);
|
|
CHECK(result.crystals == 0);
|
|
}
|
|
|
|
TEST_CASE("Reindex: an operator the crystals' own cells decide is left to them", "[reindex]") {
|
|
// Orthorhombic with a and b 3% apart: a twin law of the metric to within the obliquity tolerance,
|
|
// but every crystal's refined cell says which axis is which - nothing for the data to decide.
|
|
const UnitCell cell{78.0f, 80.4f, 37.0f, 90.0f, 90.0f, 90.0f};
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(cell).SpaceGroupNumber(16); // P222
|
|
const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1());
|
|
REQUIRE_FALSE(laws.empty());
|
|
|
|
std::vector<IntegrationOutcome> outcomes(50);
|
|
std::mt19937 rng(3);
|
|
for (auto &o : outcomes) {
|
|
o.latt = CrystalLattice(cell);
|
|
for (int i = 0; i < 100; ++i) {
|
|
Reflection r{};
|
|
r.h = static_cast<int>(rng() % 15) - 7; r.k = static_cast<int>(rng() % 15) - 7; r.l = static_cast<int>(rng() % 7) + 1;
|
|
r.d = 3.0f; r.I = static_cast<float>(rng() % 1000); r.sigma = 1.0f; r.partiality = 1.0f; r.prescaling_corr = 1.0f;
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 1, true);
|
|
REQUIRE(result.metric_preference.size() == laws.size());
|
|
for (double p : result.metric_preference)
|
|
CHECK(p == 1.0);
|
|
CHECK(result.crystals == 0);
|
|
CHECK(result.reindexed == 0);
|
|
}
|
|
|
|
TEST_CASE("Reindex: where the cells decide, the crystals whose cell disagrees are decided by the intensities",
|
|
"[reindex]") {
|
|
// As above, a and b 3% apart, but a few crystals were indexed on swapped axes: their own cells say so,
|
|
// and the intensities, read against the crystals the cells anchor, put them back.
|
|
const UnitCell cell{78.0f, 80.4f, 37.0f, 90.0f, 90.0f, 90.0f};
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(cell).SpaceGroupNumber(16); // P222
|
|
const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1());
|
|
REQUIRE(laws.size() == 1);
|
|
|
|
const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1());
|
|
std::mt19937 rng(11);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> noise(0.0, 0.3);
|
|
std::unordered_map<uint64_t, double> structure;
|
|
std::vector<std::array<int, 3>> hkls;
|
|
for (int h = -12; h <= 12; ++h)
|
|
for (int k = -12; k <= 12; ++k)
|
|
for (int l = -6; l <= 6; ++l)
|
|
if (h != 0 || k != 0 || l != 0) {
|
|
hkls.push_back({h, k, l});
|
|
structure.try_emplace(key(h, k, l).pack(), wilson(rng));
|
|
}
|
|
|
|
const int n_crystals = 300;
|
|
std::vector<IntegrationOutcome> outcomes(n_crystals);
|
|
std::vector<bool> swapped(n_crystals, false);
|
|
for (int c = 0; c < n_crystals; ++c) {
|
|
swapped[c] = c % 12 == 5; // 25 of 300: the cells still decide the law
|
|
outcomes[c].latt = swapped[c] ? ReindexLattice(CrystalLattice(cell), laws[0]) : CrystalLattice(cell);
|
|
for (int i = 0; i < 150; ++i) {
|
|
const auto &t = hkls[rng() % hkls.size()];
|
|
const auto m = swapped[c] ? laws[0].apply_to_hkl({{t[0], t[1], t[2]}}) : gemmi::Op::Miller{{t[0], t[1], t[2]}};
|
|
Reflection r{};
|
|
r.h = m[0]; r.k = m[1]; r.l = m[2];
|
|
r.d = 1.0f / std::sqrt(t[0] * t[0] / (cell.a * cell.a) + t[1] * t[1] / (cell.b * cell.b)
|
|
+ t[2] * t[2] / (cell.c * cell.c));
|
|
r.I = static_cast<float>(structure.at(key(t[0], t[1], t[2]).pack()) * (1.0 + noise(rng)));
|
|
r.sigma = 0.1f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
outcomes[c].reflections.push_back(r);
|
|
}
|
|
}
|
|
|
|
const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/true);
|
|
REQUIRE(result.metric_preference.size() == 1);
|
|
CHECK(result.metric_preference[0] >= METRIC_DECIDED);
|
|
CHECK(result.crystals == 25);
|
|
for (int c = 0; c < n_crystals; ++c)
|
|
CHECK((result.choice[c] != 0) == swapped[c]);
|
|
}
|
|
|
|
TEST_CASE("Reindex: with a reference, a crystal its own cell anchors keeps that indexing", "[reindex]") {
|
|
// a and b 3% apart: the crystals' cells decide the twin law. The reference is written in the other
|
|
// indexing, so image by image it moves every crystal; afterwards, the crystals whose cell says which
|
|
// axis is which go back, and only those whose cell cannot tell (a = b) stay where the reference put them.
|
|
const UnitCell cell{78.0f, 80.4f, 37.0f, 90.0f, 90.0f, 90.0f};
|
|
const UnitCell undecided{79.2f, 79.2f, 37.0f, 90.0f, 90.0f, 90.0f};
|
|
DiffractionExperiment x;
|
|
x.SetUnitCell(cell).SpaceGroupNumber(16); // P222
|
|
const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1());
|
|
REQUIRE(laws.size() == 1);
|
|
|
|
const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1());
|
|
std::mt19937 rng(5);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> noise(0.0, 0.2);
|
|
std::unordered_map<uint64_t, double> structure;
|
|
std::vector<std::array<int, 3>> hkls;
|
|
for (int h = -12; h <= 12; ++h)
|
|
for (int k = -12; k <= 12; ++k)
|
|
for (int l = -6; l <= 6; ++l)
|
|
if (h != 0 || k != 0 || l != 0) {
|
|
hkls.push_back({h, k, l});
|
|
structure.try_emplace(key(h, k, l).pack(), wilson(rng));
|
|
}
|
|
const auto d_of = [&](const std::array<int, 3> &t) {
|
|
return 1.0f / std::sqrt(t[0] * t[0] / (cell.a * cell.a) + t[1] * t[1] / (cell.b * cell.b)
|
|
+ t[2] * t[2] / (cell.c * cell.c));
|
|
};
|
|
std::vector<MergedReflection> reference;
|
|
std::unordered_set<uint64_t> seen;
|
|
for (const auto &t : hkls) {
|
|
if (!seen.insert(key(t[0], t[1], t[2]).pack()).second)
|
|
continue;
|
|
const auto m = laws[0].apply_to_hkl({{t[0], t[1], t[2]}});
|
|
MergedReflection r;
|
|
r.h = t[0]; r.k = t[1]; r.l = t[2];
|
|
r.I = static_cast<float>(structure.at(key(m[0], m[1], m[2]).pack()));
|
|
r.sigma = 1.0f;
|
|
r.d = d_of(t);
|
|
reference.push_back(r);
|
|
}
|
|
|
|
const int n_crystals = 200;
|
|
std::vector<IntegrationOutcome> outcomes(n_crystals);
|
|
std::vector<bool> open(n_crystals);
|
|
for (int c = 0; c < n_crystals; ++c) {
|
|
open[c] = c % 10 == 3;
|
|
outcomes[c].latt = CrystalLattice(open[c] ? undecided : cell);
|
|
for (int i = 0; i < 150; ++i) {
|
|
const auto &t = hkls[rng() % hkls.size()];
|
|
Reflection r{};
|
|
r.h = t[0]; r.k = t[1]; r.l = t[2];
|
|
r.d = d_of(t);
|
|
r.I = static_cast<float>(structure.at(key(t[0], t[1], t[2]).pack()) * (1.0 + noise(rng)));
|
|
r.sigma = 0.1f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
outcomes[c].reflections.push_back(r);
|
|
}
|
|
}
|
|
const auto as_indexed = outcomes;
|
|
const auto as_indexed_now = [&](int c) {
|
|
for (size_t i = 0; i < outcomes[c].reflections.size(); ++i) {
|
|
const auto &r = outcomes[c].reflections[i], &r0 = as_indexed[c].reflections[i];
|
|
if (r.h != r0.h || r.k != r0.k || r.l != r0.l)
|
|
return false;
|
|
}
|
|
return true;
|
|
};
|
|
|
|
const ReindexAmbiguityResolver resolver(x, reference);
|
|
for (auto &o : outcomes) {
|
|
o.reference_cc = resolver.Resolve(o.reflections, o.latt);
|
|
REQUIRE(o.reference_cc.size() == 2);
|
|
}
|
|
for (int c = 0; c < n_crystals; ++c)
|
|
REQUIRE_FALSE(as_indexed_now(c)); // the reference moved them all
|
|
|
|
const auto result = AnchorReferenceReindex(outcomes, cell, laws);
|
|
REQUIRE(result.metric_preference.size() == 1);
|
|
CHECK(result.metric_preference[0] >= METRIC_DECIDED);
|
|
CHECK(result.overruled == static_cast<size_t>(std::count(open.begin(), open.end(), false)));
|
|
for (int c = 0; c < n_crystals; ++c) {
|
|
CHECK(as_indexed_now(c) == !open[c]);
|
|
CHECK(outcomes[c].reference_cc.empty());
|
|
if (!open[c])
|
|
CHECK(outcomes[c].latt.GetUnitCell().b == Catch::Approx(cell.b).epsilon(1e-4));
|
|
}
|
|
}
|