[Review only] Resolve indexing ambiguities by correlation

[tjlane]

! This PR was vibe-coded.

Important

This is a review-only aggregate view of the existing stacked PRs. It adds no new commits and should not be merged directly. The smaller PRs remain the merge path.

Stack context

This draft collects only the reindexing layers of #31 in one diff:

• check when alternative indexing is possible: Check when alternative indexing is possible #351
• match and normalize reflections for comparison: Match and normalize reflections for reindexing comparisons #352
• choose and apply the best reindexing by correlation: Choose and apply the best reindexing by correlation #353
• reject unsuitable data and unclear results: Reject unsuitable data and unclear reindexing results #354

phase_alignment.py, has_origin_shift_ambiguity(), and their tests are deliberately excluded from this diff. They begin in #364 and are reviewed in #362.

What this implements

• detects whether a dataset has a metric-compatible indexing ambiguity
• enumerates proper-hand alternative indexing operators using Gemmi
• scores identity and alternative operators by correlation on one common reflection set
• converts amplitudes to intensities and applies equal-count resolution-bin normalization before scoring
• returns the transformed DataSet, selected operator, candidate scores, runner-up score, and correlation gap
• rejects malformed, unmerged, non-isomorphous, anomalous, duplicated, or otherwise unscorable inputs with explicit errors
• warns for weak correlations and rejects low or non-unique solutions

Reviewer focus

• operator direction and the deliberate exclusion of hand inversion
• whether every operator is scored on the same finite common-HKL intersection
• amplitude-to-intensity conversion and resolution normalization
• candidate/result semantics, especially runner-up and gap handling
• conservative validation for duplicate ASU reflections, constant data, and insufficient overlap

Please leave line-level review on the constituent PR where possible; this draft is intended to make the complete reindexing story easy to inspect in one place.

Related: #31, #174

[codecov]

Codecov Report

❌ Patch coverage is 98.36735% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 87.77%. Comparing base (9b8b06a) to head (0ec98c7).

Files with missing lines	Patch %	Lines
reciprocalspaceship/algorithms/reindexing.py	98.30%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #361      +/-   ##
==========================================
+ Coverage   86.86%   87.77%   +0.90%     
==========================================
  Files          40       42       +2     
  Lines        2856     3100     +244     
==========================================
+ Hits         2481     2721     +240     
- Misses        375      379       +4     

Flag	Coverage Δ
unittests	87.77% <98.36%> (+0.90%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[kmdalton]

please refactor phase_alignment.py and corresponding tests out of this review. it belongs in #362

[tjlane]

Thanks—addressed. phase_alignment.py, has_origin_shift_ambiguity(), and their tests now enter the stack in #364, so they are absent from this reindexing-only diff and included in #362. I also updated both aggregate summaries to make that boundary explicit. Please take another look when convenient.

[minhuanli]

Thanks for putting the aggregate view together, it made the full reindexing story much easier to follow. Below is a high-level review of the whole stack. Each item is tagged with the constituent PR where it was introduced, so fixes can land in the right place.

Summary of what the stack does

• has_reindexing_ambiguity(ds) reports whether the dataset's space group and unit cell allow an alternative (proper) indexing, using Gemmi's twin-law search.
• reindex_by_correlation(dataset, reference, ...):
  1. lists the identity plus the Gemmi operators;
  2. maps the moving Miller indices through each operator into the ASU;
  3. scores every candidate on one shared set of finite reflections, after squaring F to I and normalizing by resolution bin;
  4. picks the highest Pearson CC, and rejects a result with a low CC or a small gap to the runner-up.

What works

For merged, non-anomalous data with exact (merohedral) ambiguities, it works well. I took the tests/data/fmodel structures, scrambled each with its known operator and asked the code to recover it. With the gap check off, it picked the correct operator in all 24 structure/operator combinations, including a non-involutive y,z,x case, so the operator direction is right. 6OVT (P6₁) and 5W79 (P3₂21) pass with the default thresholds, and 6OVT still passes with moderate noise. The test that dataset.py's concat is not shadowed at import is a good catch.

Bugs that need fixing before merge

A. The pseudo-merohedral option (max_obliquity > 0) fails (#351, #353)

ds = rs.read_mtz("tests/data/fmodel/6GL4.mtz")[["FMODEL"]].dropna()   # P 1 2 1, a≈b
rs.algorithms.has_reindexing_ambiguity(ds, max_obliquity=5)            # -> True
rs.algorithms.reindex_by_correlation(ds, ds, data_key="FMODEL",
                                     reference_key="FMODEL", max_obliquity=5)
# PhaseAlignmentInputError: merged data must contain unique Miller indices after mapping to the ASU

It fails even when a dataset is compared with itself. Gemmi returns twin laws, and a twin law is not always an indexing ambiguity. Two of the three operators here (-y,-x,-z, y,-x,z) do not preserve P121's rotations, so applying them sends distinct reflections to the same ASU index. For the same reason, has_reindexing_ambiguity answers "yes" when there is no usable ambiguity.

Suggested fix: keep only operators op where op·G·op⁻¹ = G for the space group's rotation part, or drop pseudo-merohedral support from the docstrings for now. Please also add an end-to-end test with max_obliquity > 0; none exists today.

B. Coarse bins make wrong operators score too high, so clean data gets rejected (#352)

_resolution_normalize puts at least 100 reflections in each bin and caps the count at 20 bins. Inside each wide bin the intensity still falls with resolution, and every candidate shares that trend. It adds correlation to right and wrong operators alike and shrinks the gap between them. Runner-up CC with the default bins compared with about 20 reflections per bin:

data	default bins	~20 refl/bin
1CTJ	0.976	0.649
6S34	0.993	0.828
5W79	0.437	0.039

Even random Wilson-distributed intensities with the same HKLs give about 0.2 on the wrong operator, where it should be about 0. As a result, noise-free copies of 1CTJ, 6ITG, 6DWF and 6H64, reindexed by their twin operator, raise NoClearSolutionError: the true operator scores 1.000, but the gap is only 0.01–0.03. These test files go only to 8 Å, which exaggerates the effect. Some of the gap is also genuine pseudosymmetry (6ITG, 6DWF, 6H64). Still, the binning weakness is real. Suggested fix: use finer bins or a smooth normalization. rs.algorithms.scale_merged_intensities.mean_intensity_by_resolution already provides a kernel-smoothed ⟨I⟩(d).

C. Negative intensities crash the call (#352)

Merged, unscaled intensities (e.g. Aimless IMEAN) often have an outer shell whose mean is ≤ 0. That currently raises intensities must have a positive finite mean in every resolution bin. To reproduce: take 6OVT and replace the outer 10% of reflections with near-zero-mean noise, or add heavy noise to the moving intensities. Please skip or down-weight such bins, or clearly require amplitudes or French–Wilson-treated intensities.

Caveats to decide on or document

• Thresholds (Choose and apply the best reindexing by correlation #353): the 0.05 gap and 0.2 minimum are said to be calibrated on 6OVT only, and that calibration isn't in the repo. A fixed gap works against exactly the pseudosymmetric crystals where ambiguity matters most. Please add the calibration as a test or script, document its limits, or consider a relative measure.
• Anomalous columns (Reject unsuitable data and unclear reindexing results #354): any F(+)/I(+) column makes the call fail, even when data_key is IMEAN, so a typical merged MTZ like tests/data/algorithms/HEWL_SSAD_24IDC.mtz can't be used as-is. Either document this or ignore or transform those columns. Rejecting HL columns makes sense.
• Hand inversion is excluded on purpose. That is correct for chiral macromolecules, but please say so explicitly as a design decision.
• Phase origin: the returned dataset may still be origin-shifted relative to the reference, and that is resolved only in Check when alternative origin shifts are possible #364. Anyone who uses Check when alternative indexing is possible #351–Reject unsuitable data and unclear reindexing results #354 alone gets unaligned phases.
• The result quietly takes the reference's cell and space group, and the moving dataset's own cell is discarded.
• The isomorphism check runs on the unreindexed cells with the default 5% tolerance.
• With many candidates, the shared reflection set can shrink a lot. For example, cubic-metric P1 has 23 alternative operators.
• The correlation has no sigma weighting.

API, structure and docs, for consistency with the rest of the repo

1. Error classes (Check when alternative indexing is possible #351, Choose and apply the best reindexing by correlation #353). The repo doesn't define custom exceptions elsewhere; merge, compute_intensity_from_structurefactor and scale_merged_intensities all raise ValueError. If we keep new classes, they shouldn't be named PhaseAlignment* in a reindexing API, and _errors.py shouldn't describe itself as being for phase alignment. Also, has_reindexing_ambiguity raises a plain ValueError for a non-DataSet input while everything else raises PhaseAlignmentInputError. The error classes are in __all__ but not in docs/api.
2. has_reindexing_ambiguity overlaps existing API (Check when alternative indexing is possible #351). It is essentially bool(ds.reindexing_ops) / ds.find_twin_laws(...) plus validation. Consider dropping it or making it a DataSet method, since the repo's convention is that user-facing helpers live on DataSet.
3. Reuse existing helpers (Match and normalize reflections for reindexing comparisons #352, Choose and apply the best reindexing by correlation #353): DataSet.assign_resolution_bins / bin_by_percentile for binning, mean_intensity_by_resolution for smooth normalization, and DataSet.compute_dHKL instead of calling cell.calculate_d_array directly.
4. Signature style (Choose and apply the best reindexing by correlation #353). Existing algorithms take positional keys with defaults (intensity_key="I", sf_key) and return a DataSet or a tuple. Required keyword-only data_key/reference_key and a dataclass result would be new patterns. They're fine, but let's agree on them explicitly since they set the precedent for [Review only] Align datasets to a common origin #362 and Check when alternative origin shifts are possible #364. ReindexingCandidate might work as an internal type rather than public API.
5. Typing and validation weight (Check when alternative indexing is possible #351–Reject unsuitable data and unclear reindexing results #354). Final, TypeAlias, the TYPE_CHECKING import (the one line tests never run) and the checks for complex or int32-overflowing Miller indices go well beyond the rest of algorithms/, which is lightly typed and lightly validated. Trimming them would make the code easier to maintain.
6. Docstrings (Choose and apply the best reindexing by correlation #353):
   • The link to align_phases points at code that doesn't exist in this part of the stack.
   • Raises should list the actual exception types.
   • Please cite a reference for "following the strategy used by Pointless" (Evans 2006/2011).
   • max_obliquity shouldn't promise pseudo-merohedral support until A is fixed.
   • Minor: the repo's docstrings write rs.DataSet and "(default: X)".

Status

Requesting changes on the stack, but not on this aggregate, which isn't meant to be merged. A, B and C are blocking. The caveats should be decided or documented, and the API points in the last section can be settled together.

Note: a line-by-line review on the constituent PRs (#351–#354) will follow once the bugs above (A, B and C) are fixed.