diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 9eb12b27..1c172a2f 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -17,6 +17,8 @@ ADD_LIBRARY(JFJochScaleMerge RfreeFlags.h FrenchWilson.cpp FrenchWilson.h + ReindexAmbiguity.cpp + ReindexAmbiguity.h ScalingResult.h ScalingResult.cpp) TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon) diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp new file mode 100644 index 00000000..929f641e --- /dev/null +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ReindexAmbiguity.h" + +#include +#include + +#include "gemmi/twin.hpp" + +#include "HKLKey.h" + +std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, + double max_obliquity_deg) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); + if (sg == nullptr) + return {}; + return gemmi::find_twin_laws(static_cast(cell), sg, max_obliquity_deg, + /*all_ops=*/false); +} + +std::vector ReindexReflections(const std::vector &merged, + const gemmi::Op &op) { + std::vector out = merged; + for (auto &r : out) { + const gemmi::Op::Miller h = op.apply_to_hkl({{static_cast(r.h), + static_cast(r.k), + static_cast(r.l)}}); + r.h = h[0]; + r.k = h[1]; + r.l = h[2]; + } + return out; +} + +ReindexChoice ChooseReindex(const std::vector &merged, + const UnitCell &cell, int space_group_number, + const std::function &)> &score, + double max_obliquity_deg) { + ReindexChoice best; + best.identity_score = score(merged); + best.score = best.identity_score; + + const auto ops = ReindexAmbiguityOperators(cell, space_group_number, max_obliquity_deg); + best.n_candidates = 1 + static_cast(ops.size()); + for (const auto &op : ops) { + const double s = score(ReindexReflections(merged, op)); + if (s > best.score) { + best.score = s; + best.op = op; + best.is_identity = false; + } + } + return best; +} + +double ReferenceIntensityCC(const std::vector &merged, + const std::vector &reference, + int space_group_number) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); + if (sg == nullptr) + return 0.0; + const HKLKeyGenerator key(/*merge_friedel=*/true, *sg); + + std::unordered_map ref; + ref.reserve(reference.size()); + for (const auto &r : reference) + if (std::isfinite(r.I)) + ref[key(r).pack()] = r.I; + + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + int n = 0; + for (const auto &m : merged) { + if (!std::isfinite(m.I)) + continue; + const auto it = ref.find(key(m).pack()); + if (it == ref.end()) + continue; + const double x = m.I, y = it->second; + sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; + ++n; + } + if (n < 10) + return 0.0; + const double cov = n * sxy - sx * sy; + const double vx = n * sxx - sx * sx; + const double vy = n * syy - sy * sy; + return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : 0.0; +} diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h new file mode 100644 index 00000000..300276d3 --- /dev/null +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -0,0 +1,52 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "gemmi/symmetry.hpp" + +#include "../../common/Reflection.h" +#include "../../common/UnitCell.h" + +// Resolving the indexing ("alternative indexing" / reindexing) ambiguity: when the lattice symmetry +// is higher than the crystal Laue symmetry (merohedral / pseudo-merohedral), a crystal can be validly +// indexed in several ways related by reindexing operators. The right one is picked by agreement with a +// reference (reference-MTZ intensities or a model), so no reference-free ambigator (cosym) is needed. +// +// This is distinct from the screw-axis / enantiomorph ambiguity (P4_1 vs P4_3): that leaves the merged +// intensities unchanged, so it is never resolved here — it is adopted from the reference/model. + +// The candidate reindexing operators = twin laws = cosets of the crystal point group in the lattice +// symmetry. Empty when the lattice symmetry equals the crystal Laue symmetry (holohedral crystal), +// i.e. there is no indexing ambiguity. Identity is NOT included in the returned list. +std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, + double max_obliquity_deg = 2.0); + +// A copy of `merged` with each Miller index transformed by `op` (the intensities are unchanged). +std::vector ReindexReflections(const std::vector &merged, + const gemmi::Op &op); + +struct ReindexChoice { + gemmi::Op op = gemmi::Op::identity(); // the winning operator (identity = keep the current indexing) + bool is_identity = true; + double score = 0.0; // score of the winning operator + double identity_score = 0.0; // score with no reindex, for comparison / logging + int n_candidates = 1; // identity + twin laws that were tried +}; + +// Pick the reindexing operator (identity or one of the twin laws) that MAXIMIZES `score(reindexed)`. +// `score` is higher-is-better (e.g. correlation against a reference, or minus the R-free against a +// model). Returns identity when there is no ambiguity or identity scores best. +ReindexChoice ChooseReindex(const std::vector &merged, + const UnitCell &cell, int space_group_number, + const std::function &)> &score, + double max_obliquity_deg = 2.0); + +// Ready-made scorer: Pearson correlation of merged intensities against a reference set, matched by the +// Friedel-merged (Laue) ASU index in `space_group_number`. Returns 0 when too few reflections match. +double ReferenceIntensityCC(const std::vector &merged, + const std::vector &reference, + int space_group_number); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e3e41957..5b1d2565 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -76,6 +76,7 @@ ADD_EXECUTABLE(jfjoch_test MergeScaleTest.cpp RfreeFlagsTest.cpp FrenchWilsonTest.cpp + ReindexAmbiguityTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp new file mode 100644 index 00000000..bfc7953f --- /dev/null +++ b/tests/ReindexAmbiguityTest.cpp @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" + +namespace { + 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 DistinctReference(int space_group_number) { + const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(space_group_number)); + std::vector v; + std::unordered_set 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(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(), 96).empty()); + // P422 likewise. + CHECK(ReindexAmbiguityOperators(Tetragonal(), 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(), 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); + 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, + [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, 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, + [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, sg); }); + CHECK(choice.is_identity); + CHECK(choice.score > 0.99); +}