From 8bf3d16b9049dd71e7c4577df094e3b6df085f52 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 9 Nov 2026 16:14:55 +0100 Subject: [PATCH] Serial stills: the indexing ambiguity weighs each crystal's intensities and its own cell, never a lock The cell anchor (AnchorOnCells) let the crystals' own cells decide an operator once 90% of them fitted the run's cell better as indexed, and then locked every crystal whose cell preferred its indexing. On a P21 crystal with two pseudo-equal axes 0.2% apart that preference was close to 100% by construction: the per-frame cells are refined under the monoclinic constraint in the indexing they were found in (alpha and gamma held at 90) and put on an obtuse beta afterwards, so CellDistance(reindexed) - CellDistance(as indexed) = 2 (beta_f - 90)(beta_run - 90) >= 0 for every frame. Measured: 95-100% "metric preference" on all eleven pseudo-cubic operators, ~95% of the crystals locked in their random indexing, 135 of ~1600 left to the intensities. The cells scatter by ~1.2% per axis, six times the 0.2% between the axes. A -z reference could not overrule the anchor, and --mode scale -z did no reindexing against the reference. Now each crystal's evidence is a log-likelihood of each indexing: - CellMetric: its cell against the run's after each indexing, by shape (log lengths less their mean) and by each angle's distance from 90 deg; an angle the refinement held at 90/120 deg is not measured, a cell several crystals share exactly is the unrefined prior and says nothing. The per-parameter scatter is fitted on the run by EM over the indexings, so an operator is told apart by the cells exactly as far as their scatter allows - no METRIC_TIE / METRIC_DECIDED constants. - Its intensities against all the other crystals: (n - 3)(zeta_right - zeta_wrong) atanh(CC), the Fisher-z log-odds; zeta_wrong is the mean over the other indexings (the best of many noisy CCs reads as high as the right one on weak crystals and weighed the intensities down to nothing). The Brehm-Diederichs sweep makes the crystals consistent first; the cells - and a -z reference, with its CC as a log-likelihood of the same form - then choose the run's indexing as a whole; then crystal by crystal, sequentially as in the sweep (a simultaneous update drifted off the consistent answer on weak data), each takes the indexing of highest summed log-likelihood. A crystal is decided by the data and its cell, not by the reference: weighed in per crystal, a model reference cost accuracy (97% against 99.5% with it fixing the frame), since it reads a weak crystal less clearly than the merge of all the others and shares the crystal's noise with it. The per-image Resolve at integration time stays as the start; the separate anchor step (AnchorReferenceReindex) is gone, and with it IntegrationOutcome::reference_cc. The overall relabelling is restricted to operators that keep the merge's symmetry-equivalent reflections together. Where the crystal's Laue group is not normal in the lattice's (monoclinic on a pseudo-tetragonal lattice), the other coset representatives send crystals differing by the two-fold into different indexings; the former "fewest moved" rule could take one and scramble a consistent answer (unit test: 243 of 300 crystals left consistent). MergeStills (--mode scale, live checkpoints) runs the same resolution, with the reference where given. ResolveReindexReferenceFree -> ResolveIndexingAmbiguity (it takes the reference). The live per-crystal step (ReindexOnReference) and RelabelRunIndexing use the run's fitted cell model in place of the anchors. Against model truth (per-frame CC with model Fc, frames the model decides at > 3 sigma), 4000-frame subsets of pseudo-tetragonal P21 serial runs: given a right-shape cell 57% -> 99.5%; with the refined geometry 99.7% / 98.4% (rc184 48-66%); --mode scale 64-66% -> 99.8-99.9%, and both halves of one run now merge in the same indexing; with -z 80% -> 100%, 94.5% -> 99.1%, --mode scale -z (no reindexing before) -> 99.9-100%; the run whose axes tie metrically stays 100%. Stills battery subset (7 sets, 5 with an ambiguity): no regression; the pseudo-tetragonal P21 set R_split 0.170 -> 0.151, R-free 0.286 -> 0.284. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK (cherry picked from commit 790696641baab3954ffebb71b9a043f54e5d2dc6) --- docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 6 +- docs/RUGNUX_ADVANCED.md | 13 +- docs/RUGNUX_MODEL_AND_REFERENCE.md | 15 +- docs/RUGNUX_OVERVIEW.md | 4 +- docs/RUGNUX_REPORT.md | 17 +- image_analysis/IndexAndRefine.cpp | 2 +- image_analysis/IntegrationOutcome.h | 3 - .../scale_merge/ReferenceFreeReindex.cpp | 587 ++++++++++++------ .../scale_merge/ReferenceFreeReindex.h | 133 ++-- .../scale_merge/ReindexAmbiguity.cpp | 35 +- image_analysis/scale_merge/ReindexAmbiguity.h | 20 +- rugnux/ResultReport.cpp | 6 +- rugnux/Rugnux.h | 14 +- rugnux/RugnuxScaleMerge.cpp | 79 +-- rugnux/StillsLiveMerge.cpp | 24 +- rugnux/StillsLiveMerge.h | 7 +- tests/ReindexAmbiguityTest.cpp | 316 ++++++---- 17 files changed, 788 insertions(+), 493 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 87a5a8c18..1a5a5064d 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -473,11 +473,11 @@ A reference dataset (`--reference`, MTZ or SF-mmCIF) supplies known intensities **Fix the space group and cell.** Unless overridden on the command line (`-S` for the space group, `-C` for the cell), the reference's space group is adopted and its cell is used as the soft reference cell — indexing may still drift the cell within tolerance, so a small mismatch between reference and data is absorbed rather than rejected. This applies to both stills and rotation data. -**Resolve the indexing (merohedral) ambiguity.** When the lattice symmetry is higher than the crystal's Laue symmetry (e.g. $P3$, $P4$, $P6$, $C2$), more than one indexing of the same lattice is geometrically valid, and the two solutions produce *different* merged intensities that a self-consistent scale cannot tell apart — which of them is the reference's own, only an external reference can say. The candidate reindexings are the identity together with the twin-law cosets of the metric symmetry (from the unit-cell metric and the Laue group); each is scored by the intensity correlation $\mathrm{CC}_\mathrm{ref}$ of the reindexed merge against the reference, and the data are re-merged in the best-correlating indexing. The reindex is **metric-preserving** — only the $hkl$ labels change, the cell is unchanged — and it is a no-op for a holohedral crystal, which has no twin laws (the lattice and Laue symmetry coincide). For rotation data this is done once, after the space group is determined, and the whole merge is then repeated in the chosen indexing. For stills, where each crystal is indexed independently and a merge of mixed indexings would average reflections that are not symmetry mates, it is done **per image**, at integration time: a run with a reference resolves the ambiguity image by image (the image's partiality/Lorentz-corrected intensities are correlated with the reference under each candidate operator, which is scale-invariant, and the best-correlating one is adopted; a crystal whose own cell decides the indexing keeps it, as described below). In neither workflow is the reference a **scale** target: both scale against their own data (§10.2), so $\mathrm{ISa}$ and the merging statistics come from the data alone and no cross-dataset systematic is imported. Because the stills choice is made at integration time, a later re-merge of stored reflections (`--mode scale`) cannot repair a dataset integrated without a reference; the data-only resolution below runs in the full pipeline only. Where there is no reference dataset but there is a **model** (`--model`), the reference intensities are computed from it instead - $|F_\mathrm{model}|^2$ from the atomic structure factors with a flat bulk-solvent contribution at the standard constants ($k_\mathrm{sol}=0.35$, $B_\mathrm{sol}=46$ Å$^2$), which are not fitted because there are no observations yet. Nothing is scaled against them; they serve only to rank the candidate indexings, and the correlation that does the ranking is scale-invariant. This needs the cell and space group up front (`-C` / `-S`, as serial indexing wants anyway); on rotation data the same job is done after the merge, in §14.5, where a merge exists to fit the model to properly. +**Resolve the indexing (merohedral) ambiguity.** When the lattice symmetry is higher than the crystal's Laue symmetry (e.g. $P3$, $P4$, $P6$, $C2$), more than one indexing of the same lattice is geometrically valid, and the two solutions produce *different* merged intensities that a self-consistent scale cannot tell apart — which of them is the reference's own, only an external reference can say. The candidate reindexings are the identity together with the twin-law cosets of the metric symmetry (from the unit-cell metric and the Laue group); each is scored by the intensity correlation $\mathrm{CC}_\mathrm{ref}$ of the reindexed merge against the reference, and the data are re-merged in the best-correlating indexing. The reindex is **metric-preserving** — only the $hkl$ labels change, the cell is unchanged — and it is a no-op for a holohedral crystal, which has no twin laws (the lattice and Laue symmetry coincide). For rotation data this is done once, after the space group is determined, and the whole merge is then repeated in the chosen indexing. For stills, where each crystal is indexed independently and a merge of mixed indexings would average reflections that are not symmetry mates, it is decided **per crystal**: at integration time each image is put into the indexing whose partiality/Lorentz-corrected intensities correlate best with the reference (scale-invariant), and before the merge that correlation is weighed, crystal by crystal, beside the crystal's agreement with all the others and its own cell, as described below. In neither workflow is the reference a **scale** target: both scale against their own data (§10.2), so $\mathrm{ISa}$ and the merging statistics come from the data alone and no cross-dataset systematic is imported. A re-merge of stored reflections (`--mode scale`) makes the same choice, with or without a reference. Where there is no reference dataset but there is a **model** (`--model`), the reference intensities are computed from it instead - $|F_\mathrm{model}|^2$ from the atomic structure factors with a flat bulk-solvent contribution at the standard constants ($k_\mathrm{sol}=0.35$, $B_\mathrm{sol}=46$ Å$^2$), which are not fitted because there are no observations yet. Nothing is scaled against them; they serve only to rank the candidate indexings, and the correlation that does the ranking is scale-invariant. This needs the cell and space group up front (`-C` / `-S`, as serial indexing wants anyway); on rotation data the same job is done after the merge, in §14.5, where a merge exists to fit the model to properly. -**Serial stills without a reference or model.** The data cannot say which indexing is the *true* one, but they can say which crystals are in the *same* one, and that is all a merge needs. Where the space group is fixed (`-S`) and its lattice admits alternative indexings, the crystals are made consistent with each other before anything is scaled, following Brehm & Diederichs (2014). Each crystal is compared on the square roots of its partiality- and Lorentz-corrected intensities, the sign kept: intensities follow the Wilson distribution, whose long tail lets a few very strong reflections dominate every product of two crystals, so that the strongest crystals settle among themselves and a group of them can hold one another in the wrong indexing against the weaker majority (on a large-cell serial run this left a fifth of the crystals in the wrong indexing); amplitudes have no such tail. They are put on the run's mean and spread in their resolution shell (20 shells of equal reciprocal volume; a shell past the resolution the run reaches holds only noise and may have a mean of zero or below, so it is not divided by its mean), standardised within the crystal so that its unknown scale and offset drop out, and weighted by $1/(1+(\sigma/s)^2)$, $s$ being the crystal's spread, so that a weak or barely recorded reflection counts for little. The agreement of two crystals is the sum of the products of these values over the reflections they share, and the objective is that agreement summed over every pair of crystals, each in the indexing assigned to it. That sum is $(|M|^2 - \mathrm{const})/2$ with $M$ the per-reflection sum over all crystals, so a crystal's agreement with all the others is one dot product with $M$ less its own contribution: each crystal in turn takes the indexing in which that dot product is largest, $M$ is updated, and the passes repeat until no more than one crystal in a thousand moves. A pass costs one read of the reflections, not one per pair of crystals, which is what keeps tens of thousands of crystals to seconds. The objective only grows, so it converges; it is started four times (once from the indexing as found, three times from a seeded random assignment) and the start that ends highest is kept, and of the equally consistent answers related by an overall relabelling the one that moves the fewest crystals is taken. The order of every floating-point sum is fixed (crystals in image order, sequential passes), so the result does not depend on the thread count. The run reports how clearly the crystals separated: for each crystal the correlation with the mean of all the others in its chosen indexing and in the best other one, compared by Fisher's $z$ in units of its standard error ($\sqrt{2/(n-3)}$ for $n$ shared reflections) — the median of that margin and the share of crystals above $3\sigma$. +**Serial stills without a reference or model.** The data cannot say which indexing is the *true* one, but they can say which crystals are in the *same* one, and that is all a merge needs. Where the space group is fixed (`-S`) and its lattice admits alternative indexings, the crystals are made consistent with each other before anything is scaled, following Brehm & Diederichs (2014). Each crystal is compared on the square roots of its partiality- and Lorentz-corrected intensities, the sign kept: intensities follow the Wilson distribution, whose long tail lets a few very strong reflections dominate every product of two crystals, so that the strongest crystals settle among themselves and a group of them can hold one another in the wrong indexing against the weaker majority (on a large-cell serial run this left a fifth of the crystals in the wrong indexing); amplitudes have no such tail. They are put on the run's mean and spread in their resolution shell (20 shells of equal reciprocal volume; a shell past the resolution the run reaches holds only noise and may have a mean of zero or below, so it is not divided by its mean), standardised within the crystal so that its unknown scale and offset drop out, and weighted by $1/(1+(\sigma/s)^2)$, $s$ being the crystal's spread, so that a weak or barely recorded reflection counts for little. The agreement of two crystals is the sum of the products of these values over the reflections they share, and the objective is that agreement summed over every pair of crystals, each in the indexing assigned to it. That sum is $(|M|^2 - \mathrm{const})/2$ with $M$ the per-reflection sum over all crystals, so a crystal's agreement with all the others is one dot product with $M$ less its own contribution: each crystal in turn takes the indexing in which that dot product is largest, $M$ is updated, and the passes repeat until no more than one crystal in a thousand moves. A pass costs one read of the reflections, not one per pair of crystals, which is what keeps tens of thousands of crystals to seconds. The objective only grows, so it converges; it is started four times (once from the indexing as found, three times from a seeded random assignment) and the start that ends highest is kept, and of the equally consistent answers related by an overall relabelling the one the crystals' cells favour together is taken (below) - where they say nothing, the one that moves the fewest crystals. The order of every floating-point sum is fixed (crystals in image order, sequential passes), so the result does not depend on the thread count. The run reports how clearly the crystals separated: for each crystal the correlation with the mean of all the others in its chosen indexing and in the best other one, compared by Fisher's $z$ in units of its standard error ($\sqrt{2/(n-3)}$ for $n$ shared reflections) — the median of that margin and the share of crystals above $3\sigma$. -A twin law is only an ambiguity for crystals whose own lattice cannot tell the indexings apart. On a pseudo-merohedral cell whose axes differ by a few per cent, each crystal's refined cell already says which axis is which, and the indexer has used that; letting the intensities overrule it there moves crystals that were right. So for each law the share of crystals whose refined cell fits the run's cell better as indexed than after the law is measured (a tie counting half — on an exactly merohedral lattice the share is about one half). Where it is at least 0.9 the law is decided by the cells: a crystal whose cell prefers its indexing stays there, as an anchor, and only the crystals whose cell does not are decided by the intensities, against the anchors. With a reference or a model the per-image choice above stands except where the cells decide: once all images are integrated, a crystal its own cell anchors goes back to that indexing, and the reference chooses only among the indexings left open to the others; the same data-only resolution is then run only to measure it: the report gives the number of crystals on which the data alone would have chosen otherwise. +A twin law is only partly an ambiguity where the lattice is pseudo-merohedral: a crystal's refined cell then also says something about which axis is which - but only as far as the run's cell differs under the law by more than one crystal's cell scatters, and a cell refined in the indexing it was found in carries a trace of that indexing whatever the truth (a cell refined under the monoclinic constraint keeps its right angles where the indexing put them, and the convention of an obtuse $\beta$ is taken after indexing). Measured on a P2$_1$ run with two axes 0.2% apart, a former rule that let the cells decide a law once 90% of them fitted better as indexed locked 95% of the crystals in the indexing they had been found in at random: every cell fitted better as indexed, by construction. So the cells are read as evidence of the same kind as the intensities - a log-likelihood of each indexing, added to the intensities' for each crystal, never a lock. A crystal truly in indexing $g$ has, as indexed, the run's cell after the inverse of $g$. Its cell is compared with that prediction by shape: the logarithm of each length less their mean (the operators keep the volume; the size carries the distance and wavelength errors, which move all axes together) and each angle's distance from 90° (the obtuse convention cannot then carry the indexing); an angle the crystal's refinement held at 90° or 120° was not measured and is left out, and a cell several crystals share to the last digit is the prior the refinement did not move, which says nothing. The scatter of each of these six quantities about the prediction is fitted on the run itself, by expectation-maximisation over the indexings, each equally likely - so a law is told apart by the cells exactly as far as their scatter allows, with no tolerance to set. A crystal's intensities enter as $(n-3)(\zeta_\mathrm{right}-\zeta_\mathrm{wrong})\,\mathrm{atanh}(\mathrm{CC}_g)$ for each indexing $g$, $\mathrm{CC}_g$ being its correlation on $n$ reflections with the mean of all the other crystals (or with the reference): Fisher's $z$ of a correlation has a variance of $1/(n-3)$, so with $z$ near $\zeta_\mathrm{right}$ in the right indexing and $\zeta_\mathrm{wrong}$ in a wrong one the log-odds of two indexings is $(n-3)(\zeta_\mathrm{right}-\zeta_\mathrm{wrong})(z_g-z_h)$; $\zeta_\mathrm{right}$ is the run's median $z$ in the chosen indexing and $\zeta_\mathrm{wrong}$ the median of each crystal's mean $z$ over the other indexings (not the best of them, which on weak crystals reads as high as the right one and would weigh the intensities down to nothing). Without a reference, the intensities alone first make the crystals consistent as above, the cells choose the overall relabelling, and then each crystal in turn - taken out of the running sums, as in the sweep, so that every decision is made against the others as they stand - takes the indexing of highest summed log-likelihood against the mean of all the others, until no more than one in a thousand moves. Only an operator that keeps the merge's symmetry-equivalent reflections together is an overall relabelling: where the crystal's Laue group is not normal in the lattice's (a monoclinic crystal on a pseudo-tetragonal lattice), an operator that puts the two-fold on another axis would send crystals differing by the two-fold into different indexings and scramble a consistent answer, so it is never taken as one. With a reference, each crystal's correlation with it is a log-likelihood of the same form (its weight from the run's own separation of the best indexing from the others against the reference), and the reference and the cells together choose the overall relabelling, so that the merge is written in the reference's indexing; the report gives the number of crystals put into another indexing than their correlation with the reference alone would. A crystal itself is decided by its agreement with the other crystals and its cell, not by the reference: a reference - a model's intensities most of all - reads a weak crystal less clearly than the merge of all the others, and its correlation shares the crystal's own noise with that agreement, so the two are not independent evidence. Measured against model truth on a pseudo-tetragonal P2$_1$ run: 93% of the crystals the model can decide in the right indexing when each was decided on a model reference, 97% with the reference weighed in crystal by crystal, 99.5% with the reference fixing the indexing as a whole. The log gives the fitted scatter, per law the median log-likelihood by which a crystal's cell favours the indexing it came in, and the number of crystals the cells moved against their intensities. ### 10.10 Ice rings at the scale and merge stages diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 7bc48e612..61df3afc4 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -76,9 +76,9 @@ either way, so which comes out is arbitrary. What that costs depends on the work crystals, and the report says how clearly the crystals separated (`AMBIGUITY_SEPARATION_SIGMA=`, `AMBIGUITY_DECIDED_FRACTION=`) and how many were moved (`AMBIGUITY_REINDEXED=`). The indexing the whole merge ends up in is still arbitrary against other data sets, exactly as on rotation data. - On a pseudo-merohedral cell whose axes differ by a few per cent, each crystal's own refined cell - already says which axis is which; there only the crystals whose cell does not are decided by the - intensities, and the others keep their indexing. + On a pseudo-merohedral cell each crystal's own refined cell also says something about which axis is + which, as far as the run's axes differ by more than one crystal's cell scatters; it is weighed + beside the intensities, crystal by crystal, and decides where they cannot. Every run that merges tests for the ambiguity, and where it exists and nothing resolves it (serial stills whose space group was determined by the run rather than given), says so in the log and in the @@ -100,7 +100,7 @@ How to resolve it: | Situation | What to do | | --- | --- | | **Rotation**, a reference dataset exists | `-z reference.mtz`. Once the space group is settled, each candidate reindexing of the merged intensities is correlated with the reference and the best-correlating one is re-merged. Only the *hkl* labels change; the cell does not | -| **Serial stills**, a reference dataset exists | `-z reference.mtz`. Resolved **per image**, at integration time, by correlating each crystal's intensities with the reference — so the merge never mixes hands in the first place. Where the crystals' own cells tell the indexings apart (pseudo-merohedry), a crystal its cell decides keeps that indexing and the reference decides the rest. The data-only resolution below runs as a check: `AMBIGUITY_DATA_DISAGREES=` | +| **Serial stills**, a reference dataset exists | `-z reference.mtz`. The crystals are made consistent with each other as in the rows below, and the reference fixes the indexing the merge is written in. A crystal is decided by its agreement with all the others and its own cell rather than by the reference, which reads a weak crystal less clearly than the merge of all the others does: `AMBIGUITY_DATA_DISAGREES=` counts the crystals put elsewhere than the reference alone would | | **Rotation**, only a model | `--model model.pdb`. The merged data are fitted to the model in each candidate indexing and the lowest R-free wins — provided the model first beats its random-orientation null (`MODEL_FIT= ACCEPTED`) and its lead over the runner-up beats the lead a random placement of the same model takes; the written reflections are then reindexed into it, so the file, the R-factors and the maps agree. Where either bar is missed the data keep the indexing they were merged in, and the report says by how much (`MODEL_INDEXING_MARGIN_SIGMA`). It is the **data** that are relabelled, never the model: every dataset of one crystal form processed against the same model comes out in that model's indexing, so the reflection files and maps of a screening campaign compare directly. The relabelling is a proper rotation of the lattice, so the merging statistics are unchanged and I(+)/I(-) keep their hand | | **Serial stills**, only a model, space group given (`-S`) | The crystals are made consistent with each other from the data alone, as in the row below - the model plays no part in the processing. Validation then fits the merge to the model in each candidate indexing, as for rotation data, and relabels the written reflections into the model's indexing where it wins | | **Serial stills**, neither, space group given (`-S`) | Nothing to do: the crystals are made consistent with each other from the data alone (`INDEXING_AMBIGUITY_RESOLVED= data`), in an indexing that is arbitrary against other data sets. A reference still settles the weakest crystals more reliably than the other crystals can, and puts the merge in its own indexing | @@ -109,8 +109,9 @@ How to resolve it: Two things to know about where the choice lands: - **`--mode scale` resolves it as the full run does.** On stills with no reference MTZ it makes the - crystals consistent from the data alone before merging; the per-image test against a reference - happens at integration time and is kept in the stored reflections. On rotation data, where + crystals consistent from the data alone before merging; with `-z` the reference fixes the indexing of + the whole, as in a full run with the reference - so one reference puts several runs, or both arms of a pump-probe run, in + the same indexing. On rotation data, where the ambiguity is one choice for the whole dataset, `--mode scale --model model.pdb` does resolve it. - **The choice reaches the written reflections**, not only the R-factors and the maps: the merged `.mtz` / `.cif` / `.hkl` (and `_unmerged.mtz`, where it is asked for) are written in the indexing diff --git a/docs/RUGNUX_MODEL_AND_REFERENCE.md b/docs/RUGNUX_MODEL_AND_REFERENCE.md index 7b7dbb471..baa27592f 100644 --- a/docs/RUGNUX_MODEL_AND_REFERENCE.md +++ b/docs/RUGNUX_MODEL_AND_REFERENCE.md @@ -203,19 +203,18 @@ This is the only point at which a model affects which observations are merged to not gated by the null — there is no merge yet to test the model on. It chooses only between indexings that are equally valid for each crystal; a model that does not belong to the crystal carries no information to choose with, so the crystals stay mixed, as they would without a model. -Where the lattice is only pseudo-merohedral and the crystals' own cells already say which indexing is -which, the cells win: a crystal its cell anchors keeps that indexing, and the model decides -only the crystals whose cell cannot tell. +Where the lattice is only pseudo-merohedral, each crystal's own cell is weighed beside its correlation +with the model, as far as the run's cells tell the indexings apart. After the merge, the model is validated as in the steps above. Without `-C` and `-S`, or on rotation data, the model is not used before the merge at all. Without a model or reference, a stills run given `-S` settles the same question from the data: it cannot say which indexing is the true one, but it can say which crystals share one, and puts each -crystal into the indexing in which it agrees best with all the others, a crystal its own cell decides -again keeping that indexing (`INDEXING_AMBIGUITY_RESOLVED= data`). Without `-S` - the space group then -determined by the run, after the merge - nothing resolves it, and the run warns. With a model, that resolution -still runs, as a check only: `AMBIGUITY_DATA_DISAGREES=` in the report counts the crystals on which -the data alone would have chosen differently from the model. +crystal into the indexing in which it agrees best with all the others, its own cell again weighed +beside that (`INDEXING_AMBIGUITY_RESOLVED= data`). Without `-S` - the space group then +determined by the run, after the merge - nothing resolves it, and the run warns. With a reference MTZ, +that same resolution runs and the reference fixes the indexing the merge is written in: `AMBIGUITY_DATA_DISAGREES=` in the report counts the +crystals the data and their cells put into another indexing than the reference alone would. ## Where a model can touch the space group diff --git a/docs/RUGNUX_OVERVIEW.md b/docs/RUGNUX_OVERVIEW.md index 1bbf10108..afd42f5e0 100644 --- a/docs/RUGNUX_OVERVIEW.md +++ b/docs/RUGNUX_OVERVIEW.md @@ -96,8 +96,8 @@ orientation-tilt post-refinement rather than a rocking curve, and where the latt merohedral indexing ambiguity the crystals are put into one indexing before the first merge: from the data alone, each crystal into the indexing in which it agrees best with all the others, where the space group is given (`-S`); or per image against a reference or a model, which then decides the indexing -the merge is written in. Either way a crystal whose own refined cell already tells the indexings apart -keeps its indexing ([Advanced ▸ the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)). The input may be +the merge is written in. Either way a crystal's own refined cell is weighed beside its intensities, as +far as the run's cells tell the indexings apart ([Advanced ▸ the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)). The input may be several files given together - Jungfraujoch masters, or the files an XFEL pipeline writes, read without a geometry file ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)). Before the images are processed the beam centre is decided on a sample of frames, the file's against the measured ones, and the diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index df7b5f3d7..ef6d418e9 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -302,18 +302,19 @@ and a run that produced a textbook data set read identically for their first thr admits alternative indexing: one sweep is indexed consistently, so it matters only against other data; the operators are **`INDEXING_AMBIGUITY_OPERATORS=`** in section 2). On serial stills the same ambiguity mixes the hands in the merge. Where the space group was given (`-S`), the crystals are made - consistent from the data alone - by their own cells where those decide the indexing, and by the - intensities otherwise - and section 2 says so: **`INDEXING_AMBIGUITY_RESOLVED= data`**, + consistent from their intensities against each other and their own cells, weighed together - and + section 2 says so: **`INDEXING_AMBIGUITY_RESOLVED= data`**, **`AMBIGUITY_CRYSTALS=`** (the crystals that took part), **`AMBIGUITY_REINDEXED=`** (of them, the ones moved), **`AMBIGUITY_SEPARATION_SIGMA=`** (per crystal, how far its chosen indexing beats the best other one in correlation with all the other crystals, in sigma; the median) and - **`AMBIGUITY_DECIDED_FRACTION=`** (the share of crystals decided at over 3 sigma). With a reference or - a model the choice is theirs, except for the crystals whose own cell decides - (`INDEXING_AMBIGUITY_RESOLVED= reference`); the data-only resolution then runs as a check, its + **`AMBIGUITY_DECIDED_FRACTION=`** (the share of crystals decided at over 3 sigma). With a reference + MTZ (`-z`) the reference fixes the indexing the merge is written in + (`INDEXING_AMBIGUITY_RESOLVED= reference`); `AMBIGUITY_CRYSTALS=`, `AMBIGUITY_SEPARATION_SIGMA=` and `AMBIGUITY_DECIDED_FRACTION=` are given as - above, and **`AMBIGUITY_DATA_DISAGREES=`** takes the place of `AMBIGUITY_REINDEXED=`: the crystals on - which the data alone would have chosen otherwise. `AMBIGUITY_TIME_S=` (with `--developer`) is the - time the data-only resolution took. None of these keys is written by `--mode scale`, which does not run it. Where the space group was determined by the run instead, the mixed merge is an + above, and **`AMBIGUITY_DATA_DISAGREES=`** takes the place of `AMBIGUITY_REINDEXED=`: the crystals the + data and their cells put into another indexing than the reference alone would. `AMBIGUITY_TIME_S=` + (with `--developer`) is the time the resolution took. `--mode scale` resolves it the same way and logs + it, but writes none of these keys. Where the space group was determined by the run instead, the mixed merge is an `INDEXING_AMBIGUITY` warning. `INDEXING_RATE=` - the share of frames the per-image indexer took - is a stills measure: a rotation run indexes the sweep as a whole, so there it is written only with `--developer`, and `NO_LATTICE` rests on whether the sweep was indexed. diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 2c4cf73b0..7d8eca32c 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -1107,7 +1107,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { return; auto scaling_start_time = std::chrono::steady_clock::now(); - outcome.reference_cc = reindex_resolver->Resolve(outcome.reflections, outcome.latt); + reindex_resolver->Resolve(outcome.reflections, outcome.latt); auto scaling_end_time = std::chrono::steady_clock::now(); msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); } diff --git a/image_analysis/IntegrationOutcome.h b/image_analysis/IntegrationOutcome.h index 2d1b37cc0..074a1f331 100644 --- a/image_analysis/IntegrationOutcome.h +++ b/image_analysis/IntegrationOutcome.h @@ -24,9 +24,6 @@ struct IntegrationOutcome { std::optional image_scale_cc_n; std::optional image_scale_g; std::optional image_scale_wedge_deg; - // Stills with a reference: the image's CC with it in each indexing it was offered, as indexed first - // (ReindexAmbiguityResolver::Resolve); empty otherwise. - std::vector reference_cc; }; // How far apart, in frames, two partials of one raw hkl may sit and still belong to the same rocking diff --git a/image_analysis/scale_merge/ReferenceFreeReindex.cpp b/image_analysis/scale_merge/ReferenceFreeReindex.cpp index 52dc2753f..a3a052eae 100644 --- a/image_analysis/scale_merge/ReferenceFreeReindex.cpp +++ b/image_analysis/scale_merge/ReferenceFreeReindex.cpp @@ -6,6 +6,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -20,16 +23,39 @@ namespace { constexpr int RESOLUTION_BINS = 20; constexpr int RESTARTS = 4; constexpr int MAX_SWEEPS = 100; + constexpr int MAX_ROUNDS = 20; + constexpr int EM_ITERATIONS = 100; constexpr double DECIDED_SIGMA = 3.0; - // Two fits closer than this (a 0.1% difference of one axis, squared) are the same fit. - constexpr double METRIC_TIE = 1e-6; - double CellDistance(const UnitCell &u, const UnitCell &ref) { - const auto rel = [](double v, double r) { return (v - r) / r; }; - const auto ang = [](double v, double r) { return (v - r) * std::numbers::pi / 180.0; }; - return rel(u.a, ref.a) * rel(u.a, ref.a) + rel(u.b, ref.b) * rel(u.b, ref.b) - + rel(u.c, ref.c) * rel(u.c, ref.c) + ang(u.alpha, ref.alpha) * ang(u.alpha, ref.alpha) - + ang(u.beta, ref.beta) * ang(u.beta, ref.beta) + ang(u.gamma, ref.gamma) * ang(u.gamma, ref.gamma); + using CellFeatures = std::array; + + // A cell's shape, as CellMetric compares it: the log of each length less their mean, and each angle's + // distance from 90 degrees (rad). + CellFeatures Features(const UnitCell &u) { + const double la = std::log(u.a), lb = std::log(u.b), lc = std::log(u.c); + const double mean = (la + lb + lc) / 3.0; + const auto fold = [](double deg) { return std::fabs(deg - 90.0) * std::numbers::pi / 180.0; }; + return {la - mean, lb - mean, lc - mean, fold(u.alpha), fold(u.beta), fold(u.gamma)}; + } + + // Which of a crystal's cell parameters its refinement measured: an angle held at 90 or 120 degrees was not. + std::array Measured(const UnitCell &u) { + const auto free = [](double deg) { return std::fabs(deg - 90.0) > 1e-3 && std::fabs(deg - 120.0) > 1e-3; }; + return {true, true, true, free(u.alpha), free(u.beta), free(u.gamma)}; + } + + // A cell to the last digit printed, to find the cells several crystals share. + std::array CellKey(const UnitCell &u) { + const auto r = [](double v) { return static_cast(std::llround(v * 1e4)); }; + return {r(u.a), r(u.b), r(u.c), r(u.alpha), r(u.beta), r(u.gamma)}; + } + + // The cell a crystal truly in indexing g has as indexed: the run's after the inverse of g. + std::vector Predictions(const UnitCell &cell, const std::vector &indexings) { + std::vector out; + for (const auto &op : indexings) + out.push_back(Features(ReindexLattice(CrystalLattice(cell), op.inverse()).GetUnitCell())); + return out; } double Median(std::vector v) { @@ -39,129 +65,180 @@ namespace { return v[v.size() / 2]; } - // Any consistent answer relabelled as a whole by one of the indexings is as consistent: take the - // one that leaves the most crystals where they are - unless the cells decide an indexing, in which - // case the crystals they anchor fix the whole. Two indexings are the same for the merge when they - // send a general reflection to the same reflection of the asymmetric unit. Returns, for each - // indexing g, the indexing a crystal in g is in after that relabel. choice: each crystal's indexing. + size_t ArgMax(const std::vector &v) { + return std::max_element(v.begin(), v.end()) - v.begin(); + } + + // Per operator, the median over the crystals whose cell says anything of how much better it fits as + // indexed than after the operator. + std::vector MetricLogOdds(const std::vector> &ll, size_t K) { + std::vector out; + for (size_t g = 1; g < K; ++g) { + std::vector v; + for (const auto &l : ll) + if (l.size() == K && std::any_of(l.begin(), l.end(), [](double x) { return x != 0.0; })) + v.push_back(l[0] - l[g]); + out.push_back(v.empty() ? 0.0 : Median(v)); + } + return out; + } + + // A consistent answer relabelled as a whole by one of the indexings is as consistent: take the one + // the crystals' cells favour together (metric: per entry of choice, the log-likelihood of each indexing, + // or empty) - and where they say nothing, the one that leaves the most crystals where they are. Two + // indexings are the same for the merge when they send a general reflection to the same reflection of + // the asymmetric unit. Returns, for each indexing g, the indexing a crystal in g is in after that + // relabel. choice: each crystal's indexing. std::vector GlobalRelabel(const std::vector &indexings, const std::vector &choice, - const HKLKeyGenerator &key_gen, bool anchored) { + const gemmi::SpaceGroup &sg, const std::vector> &metric) { + const HKLKeyGenerator key_gen(/*merge_friedel=*/true, sg); const int K = static_cast(indexings.size()); const std::array general{2, 5, 11}; const auto same = [&](const gemmi::Op::Miller &a, const gemmi::Op::Miller &b) { return key_gen(a[0], a[1], a[2]) == key_gen(b[0], b[1], b[2]); }; - const auto composed = [&](int global, int g) { // indexing g, then the overall relabel - return indexings[global].apply_to_hkl(indexings[g].apply_to_hkl(general)); - }; + // becomes[h][c]: a crystal in c after the overall relabel h + std::vector> becomes(K, std::vector(K)); + for (int h = 0; h < K; ++h) + for (int c = 0; c < K; ++c) { + becomes[h][c] = c; + const auto composed = indexings[h].apply_to_hkl(indexings[c].apply_to_hkl(general)); + for (int g = 0; g < K; ++g) + if (same(composed, indexings[g].apply_to_hkl(general))) { + becomes[h][c] = g; + break; + } + } int global = 0; + double global_ll = -INFINITY; size_t global_kept = 0; - for (int h = 0; h < (anchored ? 1 : K); ++h) { + for (int h = 0; h < K; ++h) { + // Only an indexing that keeps the merge's symmetry-equivalent reflections together relabels a merge + // as a whole. Where the crystal's Laue group is not normal in the lattice's (a monoclinic crystal on + // a pseudo-tetragonal lattice, say), some do not - one that puts the two-fold on another axis - and + // taken as one they would split crystals that differ by a symmetry of the crystal into different + // indexings, scrambling a consistent answer. + const auto relabelled = indexings[h].apply_to_hkl(general); + bool keeps_symmetry = true; + for (const auto &op : sg.operations().sym_ops) + keeps_symmetry = keeps_symmetry + && same(indexings[h].apply_to_hkl(op.apply_to_hkl(general)), relabelled); + if (!keeps_symmetry) + continue; + double ll = 0; size_t kept = 0; - for (const int c : choice) - kept += same(composed(h, c), general); - if (kept > global_kept) { + for (size_t i = 0; i < choice.size(); ++i) { + const int g = becomes[h][choice[i]]; + kept += g == 0; + if (!metric.empty() && metric[i].size() == static_cast(K)) + ll += metric[i][g]; + } + if (ll > global_ll || (ll == global_ll && kept > global_kept)) { + global_ll = ll; global_kept = kept; global = h; } } - std::vector becomes(K); - for (int c = 0; c < K; ++c) { - becomes[c] = c; - for (int g = 0; g < K; ++g) - if (same(composed(global, c), indexings[g].apply_to_hkl(general))) { - becomes[c] = g; - break; + return becomes[global]; + } +} + +std::vector CellMetric::LogLikelihood(const CrystalLattice &latt, const std::vector &ops) const { + std::vector indexings{gemmi::Op::identity()}; + indexings.insert(indexings.end(), ops.begin(), ops.end()); + std::vector ll(indexings.size(), 0.0); + const UnitCell u = latt.GetUnitCell(); + if (priors.contains(CellKey(u))) + return ll; + const CellFeatures x = Features(u); + const auto measured = Measured(u); + const auto prediction = Predictions(cell, indexings); + for (size_t g = 0; g < indexings.size(); ++g) + for (int j = 0; j < 6; ++j) + if (measured[j] && variance[j] > 0) { + const double r = x[j] - prediction[g][j]; + ll[g] -= r * r / (2.0 * variance[j]); + } + return ll; +} + +CellMetric FitCellMetric(const std::vector &outcomes, const UnitCell &cell, + const std::vector &ops) { + std::vector indexings{gemmi::Op::identity()}; + indexings.insert(indexings.end(), ops.begin(), ops.end()); + const size_t K = indexings.size(); + CellMetric m; + m.cell = cell; + + std::map, int> seen; + for (const auto &o : outcomes) + if (!o.reflections.empty()) + ++seen[CellKey(o.latt.GetUnitCell())]; + for (const auto &[key, n] : seen) + if (n > 1) + m.priors.insert(key); + + const auto prediction = Predictions(cell, indexings); + std::vector x; + std::vector> measured; + for (const auto &o : outcomes) { + if (o.reflections.empty()) + continue; + const UnitCell u = o.latt.GetUnitCell(); + if (m.priors.contains(CellKey(u))) + continue; + x.push_back(Features(u)); + measured.push_back(Measured(u)); + } + m.crystals = x.size(); + if (x.empty()) + return m; + + // Expectation-maximisation: each crystal in each indexing with the weight its cell gives it, the scatter + // the weighted residuals have; from all indexings equally weighted. + std::vector w(x.size() * K, 1.0 / static_cast(K)); + for (int iter = 0; iter < EM_ITERATIONS; ++iter) { + std::array sum{}, count{}; + for (size_t f = 0; f < x.size(); ++f) + for (int j = 0; j < 6; ++j) + if (measured[f][j]) { + count[j] += 1.0; + for (size_t g = 0; g < K; ++g) { + const double r = x[f][j] - prediction[g][j]; + sum[j] += w[f * K + g] * r * r; + } } + for (int j = 0; j < 6; ++j) + m.variance[j] = count[j] > 0 ? sum[j] / count[j] : 0.0; + for (size_t f = 0; f < x.size(); ++f) { + std::vector ll(K, 0.0); + for (size_t g = 0; g < K; ++g) + for (int j = 0; j < 6; ++j) + if (measured[f][j] && m.variance[j] > 0) { + const double r = x[f][j] - prediction[g][j]; + ll[g] -= r * r / (2.0 * m.variance[j]); + } + const double top = *std::max_element(ll.begin(), ll.end()); + double total = 0; + for (size_t g = 0; g < K; ++g) + total += w[f * K + g] = std::exp(ll[g] - top); + for (size_t g = 0; g < K; ++g) + w[f * K + g] /= total; } - return becomes; } + return m; } -MetricAnchors AnchorOnCells(const std::vector &outcomes, const UnitCell &cell, - const std::vector &ops) { - std::vector indexings{gemmi::Op::identity()}; - indexings.insert(indexings.end(), ops.begin(), ops.end()); - const size_t K = indexings.size(); - - // How much better an image's own refined cell fits the run's cell as indexed than after indexing g. - const auto cell_preference = [&](const IntegrationOutcome &o, size_t g) { - return CellDistance(ReindexLattice(o.latt, indexings[g]).GetUnitCell(), cell) - - CellDistance(o.latt.GetUnitCell(), cell); - }; - // Which indexings the cells decide: those the indexed images' cells prefer the identity over (a tie - // counting half) for at least METRIC_DECIDED of them. - MetricAnchors out; - std::vector decided(K, false); - for (size_t g = 1; g < K; ++g) { - double prefer = 0, n = 0; - for (const auto &o : outcomes) { - if (o.reflections.empty()) - continue; - const double p = cell_preference(o, g); - prefer += std::fabs(p) < METRIC_TIE ? 0.5 : p > 0 ? 1.0 : 0.0; - n += 1; - } - out.metric_preference.push_back(n > 0 ? prefer / n : 0.5); - decided[g] = out.metric_preference.back() >= METRIC_DECIDED; - } - // Where the cells decide an indexing, an image whose own cell prefers the identity stays as indexed. - out.allowed.assign(outcomes.size() * K, 1); - for (size_t o = 0; o < outcomes.size(); ++o) - for (size_t g = 1; g < K; ++g) - if (decided[g] && !outcomes[o].reflections.empty() && cell_preference(outcomes[o], g) >= METRIC_TIE) - out.allowed[o * K + g] = 0; - return out; +double IntensityLogLikelihood(double cc, double n, double weight) { + if (!std::isfinite(cc) || n <= 3) + return 0.0; + return (n - 3.0) * weight * std::atanh(std::clamp(cc, -0.999, 0.999)); } -ReferenceAnchorResult AnchorReferenceReindex(std::vector &outcomes, const UnitCell &cell, - const std::vector &ops) { - std::vector indexings{gemmi::Op::identity()}; - indexings.insert(indexings.end(), ops.begin(), ops.end()); - const size_t K = indexings.size(); - // The indexing best correlated with the reference among those open: the one as indexed unless another - // beats it; a score that is not finite never wins (the rule of the per-image resolver). - const auto best = [&](const std::vector &cc, const uint8_t *open) { - size_t b = 0; - for (size_t g = 1; g < K; ++g) - if (open[g] && std::isfinite(cc[g]) && (!std::isfinite(cc[b]) || cc[g] > cc[b])) - b = g; - return b; - }; - const std::vector all_open(K, 1); - - // Each image back in the indexing it was indexed in. - std::vector by_reference(outcomes.size(), 0); - for (size_t o = 0; o < outcomes.size(); ++o) { - auto &out = outcomes[o]; - if (out.reference_cc.size() != K) - continue; - by_reference[o] = best(out.reference_cc, all_open.data()); - if (by_reference[o] != 0) - ReindexImage(out.reflections, out.latt, indexings[by_reference[o]].inverse()); - } - - const MetricAnchors anchors = AnchorOnCells(outcomes, cell, ops); - ReferenceAnchorResult result; - result.metric_preference = anchors.metric_preference; - for (size_t o = 0; o < outcomes.size(); ++o) { - auto &out = outcomes[o]; - if (out.reference_cc.size() != K) - continue; - const size_t g = best(out.reference_cc, &anchors.allowed[o * K]); - result.overruled += g != by_reference[o]; - if (g != 0) - ReindexImage(out.reflections, out.latt, indexings[g]); - out.reference_cc.clear(); // the scores were for the indexing as indexed, which is no longer known - } - return result; -} - -ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector &outcomes, - const DiffractionExperiment &x, - const UnitCell &cell, - const std::vector &all_ops, - size_t nthreads, bool apply) { +ReferenceFreeReindexResult ResolveIndexingAmbiguity(std::vector &outcomes, + const DiffractionExperiment &x, const UnitCell &cell, + const std::vector &all_ops, + const std::vector &reference, size_t nthreads) { ReferenceFreeReindexResult out; out.choice.assign(outcomes.size(), -1); @@ -173,8 +250,7 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector(indexings.size()); - const MetricAnchors anchors = AnchorOnCells(outcomes, cell, all_ops); - out.metric_preference = anchors.metric_preference; + out.model.metric = FitCellMetric(outcomes, cell, all_ops); const auto accept = [&](const Reflection &r) { return !r.on_ice_ring && AcceptReflection(r, s.GetHighResolutionLimit_A(), s.GetLowResolutionLimit_A()) @@ -236,18 +312,60 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector allowed(static_cast(N) * K); - std::vector free(N, false); - for (int c = 0; c < N; ++c) - for (int g = 0; g < K; ++g) { - allowed[c * K + g] = anchors.allowed[crystal_outcome[c] * K + g]; - free[c] = free[c] || (g > 0 && allowed[c * K + g]); + // What each crystal's own cell says of each indexing (CellMetric). + std::vector> metric(N); + ParallelFor(N, nthreads, [&](int c) { + metric[c] = out.model.metric.LogLikelihood(outcomes[crystal_outcome[c]].latt, all_ops); + }); + out.metric_log_odds = MetricLogOdds(metric, K); + out.crystals = N; + + // ... and its CC with the reference in each indexing, as a log-likelihood weighted by how far the right + // indexing stands out in the run (IntensityLogLikelihood): the best indexing against the mean of the + // others. + std::vector> by_reference(N); + if (!reference.empty()) { + const ReindexAmbiguityResolver resolver(x, reference); + std::vector scored(N); + ParallelFor(N, nthreads, [&](int c) { scored[c] = resolver.Scores(outcomes[crystal_outcome[c]].reflections, all_ops); }); + std::vector z_best, z_rest; + for (const auto &sc : scored) { + std::vector zs; + for (const double r : sc.cc) + if (std::isfinite(r)) + zs.push_back(std::atanh(std::clamp(r, -0.999, 0.999))); + if (zs.size() < 2) + continue; + std::sort(zs.begin(), zs.end(), std::greater<>()); + z_best.push_back(zs[0]); + z_rest.push_back(std::accumulate(zs.begin() + 1, zs.end(), 0.0) / static_cast(zs.size() - 1)); } - out.crystals = std::count(free.begin(), free.end(), true); - if (out.crystals == 0) - return out; + const double weight = std::max(0.0, Median(z_best) - Median(z_rest)); + // An indexing the reference cannot score counts as the crystal's worst scored one; a crystal it cannot + // score at all gets nothing from it. + for (int c = 0; c < N; ++c) { + std::vector ll(K, INFINITY); + for (int g = 0; g < K; ++g) + if (std::isfinite(scored[c].cc[g])) + ll[g] = IntensityLogLikelihood(scored[c].cc[g], static_cast(scored[c].n), weight); + const double worst = *std::min_element(ll.begin(), ll.end()); + if (!std::isfinite(worst)) + continue; + for (double &v : ll) + v = std::isfinite(v) ? v : worst; + by_reference[c] = ll; + } + } + // What chooses the run's indexing as a whole: the cells and the reference. A crystal itself is decided by + // its agreement with all the other crystals and its cell, not by the reference: a reference - a model's + // intensities most of all - reads a weak crystal less clearly than the merge of all the others does, + // and its CC shares the crystal's own noise with that agreement, so the two are not independent. + // Measured on a pseudo-tetragonal P21 serial run against model truth: 97% of the crystals right with + // the reference weighed in crystal by crystal, 99.5% with the reference fixing the indexing as a whole. + std::vector> evidence = metric; + for (int c = 0; c < N; ++c) + for (int g = 0; g < K && !by_reference[c].empty(); ++g) + evidence[c][g] += by_reference[c][g]; const size_t nnz = begin.back(); // The mean and spread of the amplitudes in shells of equal reciprocal volume over the whole run, so @@ -367,8 +485,7 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector 0) for (int c = 0; c < N; ++c) { - const int g = static_cast(rng() % K); - choice[c] = allowed[c * K + g] ? g : 0; + choice[c] = static_cast(rng() % K); } std::fill(M.begin(), M.end(), 0.0); for (int c = 0; c < N; ++c) @@ -383,7 +500,7 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector best_dot) { @@ -406,56 +523,132 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector C(n_keys, 0.0); - for (int c = 0; c < N; ++c) - for (size_t i = begin[c]; i < begin[c + 1]; ++i) { - M[key[best_choice[c]][i]] += z[i]; - C[key[best_choice[c]][i]] += 1.0; - } - std::vector sigma(N, NAN), cc_chosen(N, NAN), cc_other(N, NAN); - ParallelFor(N, nthreads, [&](int c) { - std::vector cc(K, NAN); - std::vector n_used(K, 0.0); - for (int g = 0; g < K; ++g) { - double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0, n = 0; + // Each crystal's CC with the mean of all the other crystals in every indexing, and on how many reflections. + std::vector cc(static_cast(N) * K), n_used(static_cast(N) * K); + std::vector C(n_keys); + const auto correlate = [&] { + std::fill(M.begin(), M.end(), 0.0); + std::fill(C.begin(), C.end(), 0.0); + for (int c = 0; c < N; ++c) for (size_t i = begin[c]; i < begin[c + 1]; ++i) { - const int32_t k = key[g][i]; - const double own = g == best_choice[c] ? 1.0 : 0.0; - const double cnt = C[k] - own; - if (cnt <= 0) - continue; - const double xv = z[i], yv = (M[k] - own * z[i]) / cnt; - sx += xv; sy += yv; sxx += xv * xv; syy += yv * yv; sxy += xv * yv; n += 1; + M[key[best_choice[c]][i]] += z[i]; + C[key[best_choice[c]][i]] += 1.0; } - const double vx = sxx - sx * sx / n, vy = syy - sy * sy / n; - if (n > 3 && vx > 0 && vy > 0) { - cc[g] = std::clamp((sxy - sx * sy / n) / std::sqrt(vx * vy), -0.999, 0.999); - n_used[g] = n; + ParallelFor(N, nthreads, [&](int c) { + for (int g = 0; g < K; ++g) { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0, n = 0; + for (size_t i = begin[c]; i < begin[c + 1]; ++i) { + const int32_t k = key[g][i]; + const double own = g == best_choice[c] ? 1.0 : 0.0; + const double cnt = C[k] - own; + if (cnt <= 0) + continue; + const double xv = z[i], yv = (M[k] - own * z[i]) / cnt; + sx += xv; sy += yv; sxx += xv * xv; syy += yv * yv; sxy += xv * yv; n += 1; + } + const double vx = sxx - sx * sx / n, vy = syy - sy * sy / n; + const bool ok = n > 3 && vx > 0 && vy > 0; + cc[c * K + g] = ok ? std::clamp((sxy - sx * sy / n) / std::sqrt(vx * vy), -0.999, 0.999) : NAN; + n_used[c * K + g] = ok ? n : 0.0; } - } - const int g0 = best_choice[c]; + }); + }; + // The best other indexing of each crystal, -1 where there is none to compare with. + const auto other_of = [&](int c) { int other = -1; for (int g = 0; g < K; ++g) - if (g != g0 && allowed[c * K + g] && std::isfinite(cc[g]) && (other < 0 || cc[g] > cc[other])) + if (g != best_choice[c] && std::isfinite(cc[c * K + g]) && (other < 0 || cc[c * K + g] > cc[c * K + other])) other = g; - if (!std::isfinite(cc[g0]) || other < 0) - return; - cc_chosen[c] = cc[g0]; - cc_other[c] = cc[other]; - sigma[c] = (std::atanh(cc[g0]) - std::atanh(cc[other])) - * std::sqrt((std::min(n_used[g0], n_used[other]) - 3.0) / 2.0); - }); - std::vector sig, ccc, cco; - for (int c = 0; c < N; ++c) - if (free[c] && std::isfinite(sigma[c])) { - sig.push_back(sigma[c]); - ccc.push_back(cc_chosen[c]); - cco.push_back(cc_other[c]); + return other; + }; + + // The intensities alone leave the run's indexing as a whole open: the reference and the cells choose it. + { + const std::vector becomes = GlobalRelabel(indexings, best_choice, sg, evidence); + for (int &g : best_choice) + g = becomes[g]; + } + + // Each crystal into the indexing its intensities and its cell together favour, against the mean of all + // the others as they stand, until (almost) none moves. The weight of the intensities is the run's own + // separation of the chosen indexing from the others (IntensityLogLikelihood). + for (int round = 0; round < MAX_ROUNDS; ++round) { + correlate(); + // A wrong indexing's z is the mean over the other indexings, not the best of them: the best of many + // noisy CCs reads high, and on weak crystals would weigh the intensities down to nothing. + std::vector z_chosen, z_other; + for (int c = 0; c < N; ++c) { + if (!std::isfinite(cc[c * K + best_choice[c]])) + continue; + double sum = 0, n = 0; + for (int g = 0; g < K; ++g) + if (g != best_choice[c] && std::isfinite(cc[c * K + g])) { + sum += std::atanh(cc[c * K + g]); + n += 1; + } + if (n > 0) { + z_chosen.push_back(std::atanh(cc[c * K + best_choice[c]])); + z_other.push_back(sum / n); + } } + out.model.intensity_weight = std::max(0.0, Median(z_chosen) - Median(z_other)); + // Crystal by crystal, as the sweep above, so that every decision is made against the others as they + // now stand: taken out of the sums, its CC with the mean of the others in each indexing, and back in. + int changed = 0; + out.cell_moved = 0; + for (int c = 0; c < N; ++c) { + for (size_t i = begin[c]; i < begin[c + 1]; ++i) { + M[key[best_choice[c]][i]] -= z[i]; + C[key[best_choice[c]][i]] -= 1.0; + } + // An indexing with no reflection in common with the others cannot be told and never wins; where + // there is none, the crystal stays. + std::vector intensity(K, -INFINITY), both(K, -INFINITY); + for (int g = 0; g < K; ++g) { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0, n = 0; + for (size_t i = begin[c]; i < begin[c + 1]; ++i) { + const int32_t k = key[g][i]; + if (C[k] <= 0) + continue; + const double xv = z[i], yv = M[k] / C[k]; + sx += xv; sy += yv; sxx += xv * xv; syy += yv * yv; sxy += xv * yv; n += 1; + } + const double vx = sxx - sx * sx / n, vy = syy - sy * sy / n; + if (n > 3 && vx > 0 && vy > 0) { + const double r = (sxy - sx * sy / n) / std::sqrt(vx * vy); + intensity[g] = IntensityLogLikelihood(r, n, out.model.intensity_weight); + both[g] = intensity[g] + metric[c][g]; + } + } + const size_t b = ArgMax(both); + if (std::isfinite(both[b])) { + changed += static_cast(b) != best_choice[c]; + out.cell_moved += b != ArgMax(intensity); + best_choice[c] = static_cast(b); + } + for (size_t i = begin[c]; i < begin[c + 1]; ++i) { + M[key[best_choice[c]][i]] += z[i]; + C[key[best_choice[c]][i]] += 1.0; + } + } + if (changed <= N / 1000) + break; + } + correlate(); + + // How clearly each crystal is decided by its intensities: its CC with the mean of all the other crystals, + // in the chosen indexing and in the best other one, and the difference of the two in sigma (Fisher's z; + // each CC on n reflections has an error of 1/sqrt(n - 3), and the two are nearly independent). + std::vector sig, ccc, cco; + for (int c = 0; c < N; ++c) { + const int g0 = best_choice[c], other = other_of(c); + if (other < 0 || !std::isfinite(cc[c * K + g0])) + continue; + ccc.push_back(cc[c * K + g0]); + cco.push_back(cc[c * K + other]); + sig.push_back((std::atanh(cc[c * K + g0]) - std::atanh(cc[c * K + other])) + * std::sqrt((std::min(n_used[c * K + g0], n_used[c * K + other]) - 3.0) / 2.0)); + } out.median_sigma = Median(sig); out.cc_chosen = Median(ccc); out.cc_other = Median(cco); @@ -463,40 +656,39 @@ ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector(std::count_if(sig.begin(), sig.end(), [](double v) { return v > DECIDED_SIGMA; })) / static_cast(sig.size()); - const bool anchored = std::any_of(out.metric_preference.begin(), out.metric_preference.end(), - [](double p) { return p >= METRIC_DECIDED; }); - const std::vector becomes = GlobalRelabel(indexings, best_choice, key_gen, anchored); for (int c = 0; c < N; ++c) { - const int final_g = becomes[best_choice[c]]; + const int final_g = best_choice[c]; out.choice[crystal_outcome[c]] = final_g; + if (!by_reference[c].empty()) + out.reference_overruled += final_g != static_cast(ArgMax(by_reference[c])); if (final_g == 0) continue; ++out.reindexed; - if (apply) { - auto &o = outcomes[crystal_outcome[c]]; - ReindexImage(o.reflections, o.latt, indexings[final_g]); - } + auto &o = outcomes[crystal_outcome[c]]; + ReindexImage(o.reflections, o.latt, indexings[final_g]); } return out; } -int ReindexOnReference(IntegrationOutcome &outcome, const ReindexAmbiguityResolver &resolver, const UnitCell &cell, - const std::vector &ops, const std::vector &metric_preference) { +int ReindexOnReference(IntegrationOutcome &outcome, const ReindexAmbiguityResolver &resolver, + const std::vector &ops, const AmbiguityModel &model) { if (ops.empty() || outcome.reflections.empty()) return 0; - const std::vector cc = resolver.Scores(outcome.reflections); - const auto open = [&](size_t g) { - if (g == 0 || metric_preference.size() < g || metric_preference[g - 1] < METRIC_DECIDED) - return true; - const double p = CellDistance(ReindexLattice(outcome.latt, ops[g - 1]).GetUnitCell(), cell) - - CellDistance(outcome.latt.GetUnitCell(), cell); - return p < METRIC_TIE; - }; - // The rule of AnchorReferenceReindex: as indexed unless an open indexing beats it. + const auto scored = resolver.Scores(outcome.reflections); + const std::vector metric = model.metric.LogLikelihood(outcome.latt, ops); + // As indexed unless another indexing is favoured. size_t best = 0; - for (size_t g = 1; g <= ops.size(); ++g) - if (open(g) && std::isfinite(cc[g]) && (!std::isfinite(cc[best]) || cc[g] > cc[best])) + double best_ll = -INFINITY; + for (size_t g = 0; g <= ops.size(); ++g) { + if (!std::isfinite(scored.cc[g])) // cannot be told, so never wins + continue; + const double ll = IntensityLogLikelihood(scored.cc[g], static_cast(scored.n), model.intensity_weight) + + metric[g]; + if (ll > best_ll) { + best_ll = ll; best = g; + } + } if (best != 0) ReindexImage(outcome.reflections, outcome.latt, ops[best - 1]); return static_cast(best); @@ -504,19 +696,20 @@ int ReindexOnReference(IntegrationOutcome &outcome, const ReindexAmbiguityResolv size_t RelabelRunIndexing(std::vector &outcomes, std::vector &choice, const DiffractionExperiment &x, const std::vector &ops, - const std::vector &metric_preference) { + const CellMetric &metric) { if (ops.empty()) return 0; std::vector indexings{gemmi::Op::identity()}; indexings.insert(indexings.end(), ops.begin(), ops.end()); - const HKLKeyGenerator key_gen(/*merge_friedel=*/true, x.GetSpaceGroupOrP1()); - const bool anchored = std::any_of(metric_preference.begin(), metric_preference.end(), - [](double p) { return p >= METRIC_DECIDED; }); + // Each crystal's cell, read as indexed: the lattice is now in choice[o], so it goes back first. std::vector taking_part; + std::vector> ll; for (size_t o = 0; o < outcomes.size(); ++o) - if (choice[o] >= 0 && !outcomes[o].reflections.empty()) + if (choice[o] >= 0 && !outcomes[o].reflections.empty()) { taking_part.push_back(choice[o]); - const std::vector becomes = GlobalRelabel(indexings, taking_part, key_gen, anchored); + ll.push_back(metric.LogLikelihood(ReindexLattice(outcomes[o].latt, indexings[choice[o]].inverse()), ops)); + } + const std::vector becomes = GlobalRelabel(indexings, taking_part, x.GetSpaceGroupOrP1(), ll); size_t moved = 0; for (size_t o = 0; o < outcomes.size(); ++o) { if (choice[o] < 0 || becomes[choice[o]] == choice[o]) diff --git a/image_analysis/scale_merge/ReferenceFreeReindex.h b/image_analysis/scale_merge/ReferenceFreeReindex.h index 2f783979b..e2fba08a7 100644 --- a/image_analysis/scale_merge/ReferenceFreeReindex.h +++ b/image_analysis/scale_merge/ReferenceFreeReindex.h @@ -4,7 +4,9 @@ #pragma once #include +#include #include +#include #include #include "gemmi/symmetry.hpp" @@ -14,84 +16,93 @@ class ReindexAmbiguityResolver; -// Serial stills with no reference: every crystal is indexed in one of the indexings its lattice allows -// (the identity and the twin-law operators of ReindexAmbiguityOperators), at random. Which one each -// crystal is in is decided from the data alone, by putting each crystal into the indexing in which it -// agrees best with all the others. The overall choice among the equivalent consistent answers is the -// one that moves the fewest crystals. -// -// An operator is only an ambiguity where the crystal's own lattice cannot tell the indexings apart. Where -// the cell is merely close to a more symmetric one (two axes a few per cent apart, say), the crystals' -// refined cells already say which axis is which, and the indexer has used that. For such an operator a -// crystal whose cell prefers the indexing it came in stays there, as an anchor, and only the others are -// decided by the intensities. -// An operator whose reindexed cell is the worse fit to the run's cell for at least this share of the -// crystals is one the cells decide. For a true (merohedral) ambiguity the two fits are equal and the -// share is a half. -constexpr double METRIC_DECIDED = 0.9; +// Serial stills: every crystal is indexed in one of the indexings its lattice allows (the identity and +// the twin-law operators of ReindexAmbiguityOperators), at random. Which one each crystal is in is decided +// from its intensities against all the other crystals and its own refined cell (CellMetric), added as +// log-likelihoods; neither locks a crystal: where one cannot tell two indexings apart it adds nothing, and +// where it can, it is weighed against the other. The run's indexing as a whole - which the crystals' +// agreement with each other leaves open - is chosen by the cells and, where one is given, a reference. -// Which indexings each image's own refined cell leaves open. Indexing g is the image as indexed for g = 0 -// and after ops[g - 1] otherwise; `cell` is the run's cell. An operator the cells decide is closed to -// every image whose cell fits `cell` better as indexed; everything else is open. -struct MetricAnchors { - std::vector metric_preference; // per operator: share of the images whose refined cell is - // closer to `cell` as indexed than after the operator - std::vector allowed; // allowed[o * (ops.size() + 1) + g]: indexing g open to image o -}; -MetricAnchors AnchorOnCells(const std::vector &outcomes, const UnitCell &cell, - const std::vector &ops); +// What a crystal's own refined cell says about the indexing it is in. A crystal truly in indexing g has, +// as indexed, the run's cell after the inverse of g, up to the scatter of the refined cells about it. That +// scatter is fitted on the run itself (expectation-maximisation over the indexings, each equally likely), +// so an operator is told apart by the cells only as far as the run's cells differ under it by more than +// one crystal's cell scatters. A cell is compared by its shape - each length after the overall size is +// taken out, since the operators keep the volume and the size carries the distance and wavelength errors +// - and by its angles' distance from 90 degrees, since the monoclinic convention of an obtuse angle is +// taken after indexing, whatever the true indexing. An angle a crystal's refinement held at 90 or 120 +// degrees was not measured, and a cell shared exactly by several crystals is the prior that the +// refinement did not move: neither says anything. +struct CellMetric { + UnitCell cell; // the run's + std::array variance{}; // of a crystal's cell about the prediction (log lengths, angles + // in rad), 0 where no crystal measured it + std::set> priors; // cells more than one crystal has to the last digit + size_t crystals = 0; // crystals the fit used -// With a reference - a reference MTZ or intensities from the model - each image has been put, image by -// image, into the indexing that correlates best with it (ReindexAmbiguityResolver, which leaves the scores in -// IntegrationOutcome::reference_cc, over the same `ops`). The reference knows nothing of the cells, so -// here the cells have their say first, as above: an image its own cell anchors keeps that indexing, -// and the reference decides among the indexings left open to the others. -struct ReferenceAnchorResult { - std::vector metric_preference; // as in MetricAnchors - size_t overruled = 0; // images the reference had put into an indexing their cell - // excludes, now in the best open one instead + // Per indexing (as indexed, then ops[g - 1]), the log-likelihood of the crystal's cell; all zero where + // it says nothing. + [[nodiscard]] std::vector LogLikelihood(const CrystalLattice &latt, + const std::vector &ops) const; +}; +CellMetric FitCellMetric(const std::vector &outcomes, const UnitCell &cell, + const std::vector &ops); + +// A crystal's intensities as a log-likelihood of each indexing, from its CC with the others (or with a +// reference) in it, on n reflections: (n - 3) * weight * atanh(cc). Fisher's z of a CC on n reflections has +// a variance of 1 / (n - 3); with z near zeta_right in the right indexing and near zeta_wrong in a wrong +// one, the log-odds of two indexings is (n - 3) (zeta_right - zeta_wrong) (z_g - z_h), so the weight is the +// run's zeta_right - zeta_wrong: the median z of the chosen indexing less the median of the mean z of the others. +double IntensityLogLikelihood(double cc, double n, double weight); + +// What a run's resolution leaves for the crystals that come after it (StillsLiveMerge). +struct AmbiguityModel { + CellMetric metric; + double intensity_weight = 0.0; // as in IntensityLogLikelihood }; -ReferenceAnchorResult AnchorReferenceReindex(std::vector &outcomes, const UnitCell &cell, - const std::vector &ops); struct ReferenceFreeReindexResult { - size_t crystals = 0; // crystals that took part (enough reflections to say anything, and an - // indexing their own cell leaves open) + size_t crystals = 0; // crystals that took part (enough reflections to say anything) + size_t reference_overruled = 0; // with a reference: crystals put into another indexing than their CC + // with the reference alone would size_t reindexed = 0; // of them, put into an indexing other than the one they came in double median_sigma = NAN; // per crystal, how far its chosen indexing beats the best other one // in agreement with the rest of the merge, in sigma; the median double decided_fraction = NAN; // share of the crystals for which that margin is over 3 sigma double cc_chosen = NAN; // median CC of a crystal with the merge of all the others, as chosen double cc_other = NAN; // ... and in the best other indexing - std::vector metric_preference; // per operator: share of crystals whose refined cell is - // closer to `cell` as indexed than after the operator + size_t cell_moved = 0; // crystals their cell put into another indexing than their intensities alone + std::vector metric_log_odds; // per operator: median over the crystals of how much better their + // cell fits as indexed than after it (log-likelihood) + AmbiguityModel model; int sweeps = 0; // passes over the crystals until (almost) none changed (the start kept) std::vector choice; // per outcome: -1 did not take part, 0 kept, k reindexed by ops[k - 1] }; -// `cell` is the run's cell, `ops` its twin-law operators (identity excluded). With apply = false nothing is changed - the -// result then says how many crystals the data alone would put into another indexing than they are in, -// which, after a reference has resolved them, is the disagreement with the reference. -ReferenceFreeReindexResult ResolveReindexReferenceFree(std::vector &outcomes, - const DiffractionExperiment &x, - const UnitCell &cell, - const std::vector &ops, - size_t nthreads, bool apply); +// `cell` is the run's cell, `ops` its twin-law operators (identity excluded). The intensities first make the +// crystals consistent with each other (Brehm & Diederichs); the run's indexing as a whole is then the one the +// reference and the cells favour together (where neither says anything, the one that moves the fewest +// crystals); then each crystal takes the indexing its intensities and cell favour, crystal by crystal, until +// (almost) none moves. `reference`: merged intensities in the run's cell and group, or empty. The crystals are +// reindexed in place. +ReferenceFreeReindexResult ResolveIndexingAmbiguity(std::vector &outcomes, + const DiffractionExperiment &x, const UnitCell &cell, + const std::vector &ops, + const std::vector &reference, size_t nthreads); -// A crystal that comes after a run has been resolved (ResolveReindexReferenceFree, which gave -// metric_preference): put into the indexing that correlates best with the merge of that run (`resolver`), -// among those its own cell leaves open - an operator the run's cells decided is closed to a crystal whose -// cell fits `cell` better as indexed, as in AnchorOnCells. Returns the indexing it is put in (0 as indexed, -// k after ops[k - 1]). -int ReindexOnReference(IntegrationOutcome &outcome, const ReindexAmbiguityResolver &resolver, const UnitCell &cell, - const std::vector &ops, const std::vector &metric_preference); +// A crystal that comes after a run has been resolved (ResolveIndexingAmbiguity, which gave `model`): put +// into the indexing its CC with the merge of that run (`resolver`) and its own cell together favour. Returns +// the indexing it is put in (0 as indexed, k after ops[k - 1]). +int ReindexOnReference(IntegrationOutcome &outcome, const ReindexAmbiguityResolver &resolver, + const std::vector &ops, const AmbiguityModel &model); // The run's indexing as a whole, once its crystals were put into theirs in steps (a checkpoint's -// ResolveReindexReferenceFree, then crystal by crystal with ReindexOnReference): relabelled by the indexing -// that leaves the most crystals as they were indexed, as ResolveReindexReferenceFree takes it on the whole -// run - unless the run's cells decide an indexing. choice: per outcome, the indexing it is in (0 as indexed, -// k after ops[k - 1], -1 not taking part), updated. Returns the crystals moved. Relabelling every crystal -// by one indexing changes no intensity nor partiality, only the frame the merge is reported in. +// ResolveIndexingAmbiguity, then crystal by crystal with ReindexOnReference): relabelled by the indexing +// the crystals' cells favour together, or, where they say nothing, by the one that leaves the most crystals +// as they were indexed, as ResolveIndexingAmbiguity takes it on the whole run. choice: per outcome, the +// indexing it is in (0 as indexed, k after ops[k - 1], -1 not taking part), updated. Returns the crystals +// moved. Relabelling every crystal by one indexing changes no intensity nor partiality, only the frame the +// merge is reported in. size_t RelabelRunIndexing(std::vector &outcomes, std::vector &choice, const DiffractionExperiment &x, const std::vector &ops, - const std::vector &metric_preference); + const CellMetric &metric); diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index f13869e7a..8089346dc 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -236,9 +236,9 @@ bool ReindexAmbiguityResolver::Accept(const Reflection &r) const { } double ReindexAmbiguityResolver::ReferenceCC(std::span reflections, - const gemmi::Op &op) const { + const gemmi::Op &op, size_t &n) const { double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; - size_t n = 0; + n = 0; for (const auto &r : reflections) { if (!Accept(r) || r.partiality < s.GetMinPartiality()) continue; @@ -266,27 +266,30 @@ double ReindexAmbiguityResolver::ReferenceCC(std::span reflect // Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this // image alone, the reindexing whose intensities correlate best with the external reference and apply it. -std::vector ReindexAmbiguityResolver::Resolve(std::span reflections, CrystalLattice &latt) const { +void ReindexAmbiguityResolver::Resolve(std::span reflections, CrystalLattice &latt) const { if (ops.empty()) - return {}; - - std::vector cc{ReferenceCC(reflections, gemmi::Op::identity())}; - const BestReindex best = PickBestReindex(cc[0], ops, [&](const gemmi::Op &op) { - cc.push_back(ReferenceCC(reflections, op)); - return cc.back(); - }); + return; + size_t n = 0; + const BestReindex best = PickBestReindex(ReferenceCC(reflections, gemmi::Op::identity(), n), ops, + [&](const gemmi::Op &op) { return ReferenceCC(reflections, op, n); }); if (!best.is_identity) ReindexImage(reflections, latt, best.op); - return cc; } -std::vector ReindexAmbiguityResolver::Scores(std::span reflections) const { +ReindexAmbiguityResolver::Scored ReindexAmbiguityResolver::Scores(std::span reflections, + const std::vector &with_ops) const { + Scored out; + out.cc.push_back(ReferenceCC(reflections, gemmi::Op::identity(), out.n)); + size_t n_op = 0; + for (const auto &op : with_ops) + out.cc.push_back(ReferenceCC(reflections, op, n_op)); + return out; +} + +ReindexAmbiguityResolver::Scored ReindexAmbiguityResolver::Scores(std::span reflections) const { if (ops.empty()) return {}; - std::vector cc{ReferenceCC(reflections, gemmi::Op::identity())}; - for (const auto &op : ops) - cc.push_back(ReferenceCC(reflections, op)); - return cc; + return Scores(reflections, ops); } void ReindexImage(std::span reflections, CrystalLattice &latt, const gemmi::Op &op) { diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 5fe6112e5..78ace25b8 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -131,8 +131,8 @@ class ReindexAmbiguityResolver { bool Accept(const Reflection &r) const; // Pearson correlation of this image's partiality/Lorentz-corrected intensities against the reference, // matched by ASU index under `op`. Scale-invariant (the unknown per-image scale cancels). NaN when - // too few reflections match. - double ReferenceCC(std::span reflections, const gemmi::Op &op) const; + // too few reflections match. n: the reflections matched. + double ReferenceCC(std::span reflections, const gemmi::Op &op, size_t &n) const; public: // The reference must be an external, correctly-handed dataset (not the data's own running merge). @@ -141,10 +141,14 @@ public: // Reindex this image's reflections in place into the hand best-correlated with the reference. No-op // when there is no ambiguity. // Relabels the image's reflections and moves its lattice onto the same axes, so that a spot's new - // index still predicts it where it was measured. Returns the CC in each indexing, as indexed first and - // then the operators of ReindexAmbiguityOperators(cell, space group) - empty when there are none - - // for AnchorReferenceReindex to revisit the choice once the whole run's cells are known. - std::vector Resolve(std::span reflections, CrystalLattice &latt) const; - // The CC in each indexing, as Resolve returns them, with nothing changed. - std::vector Scores(std::span reflections) const; + // index still predicts it where it was measured. + void Resolve(std::span reflections, CrystalLattice &latt) const; + // The CC in each indexing - as indexed, then after each of `ops` - with nothing changed. + struct Scored { + std::vector cc; + size_t n = 0; // reflections the CC as indexed was computed on + }; + Scored Scores(std::span reflections, const std::vector &ops) const; + // ... after the operators of ReindexAmbiguityOperators(cell, space group); empty where there are none. + Scored Scores(std::span reflections) const; }; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index c35ad1a0b..cac4f624c 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1136,9 +1136,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyBool("LATTICE_FOUND", result.consensus_cell.has_value())); if (!result.indexing_ambiguity_operators.empty()) Add(s, KeyText("INDEXING_AMBIGUITY_OPERATORS", result.indexing_ambiguity_operators)); - // Serial stills: the indexing ambiguity from the data alone - deciding the merge where no reference - // or model was given, otherwise measured against the reference's answer (REINDEXED is then the - // number of crystals on which the two disagree). + // Serial stills: the indexing ambiguity, resolved from the data - and against the reference where one + // was given (DATA_DISAGREES is then the number of crystals the data and their cells put into another + // indexing than the reference alone would). if (result.ambiguity_reference_free) { Add(s, KeyText("INDEXING_AMBIGUITY_RESOLVED", result.ambiguity_reference_free_applied ? "data" : "reference")); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index b18f18c0e..ade3887ac 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -29,6 +29,7 @@ #include "../image_analysis/spot_finding/SpotFindingSettings.h" #include "../image_analysis/spot_finding/SpotUtils.h" // PowderRings #include "../image_analysis/scale_merge/Merge.h" // MergeStatistics +#include "../image_analysis/scale_merge/ReferenceFreeReindex.h" // CellMetric #include "../image_analysis/scale_merge/TwinningAnalysis.h" // TwinningAnalysisResult #include "../image_analysis/scale_merge/CentreOfSymmetry.h" // CentreOfSymmetryResult #include "../image_analysis/scale_merge/StillsPartialityRefine.h" // StillsPartialityRefine::RunWidth @@ -383,6 +384,11 @@ OtherCellForm LeaveOutOtherCellForms(std::vector &outcomes, OtherCellForm LeaveOutMinorCellForm(std::vector &outcomes, const UnitCell &cell, const IndexingSettings &settings, Logger &logger); +// What the crystals' own cells say of the indexing ambiguity (CellMetric): their scatter, and per operator the +// median log-likelihood of as indexed over the operator (metric_log_odds). +void LogCellMetric(Logger &logger, const CellMetric &metric, const std::vector &ops, + const std::vector &metric_log_odds); + // The stills scale and merge, the one both --mode mx and --mode scale make: each image scaled against the // merge of all, the partialities post-refined against it (unless for_search or --simple-stills), the error // model fitted, the observations merged, cut at the resolution fall-off and the statistics taken. @@ -843,10 +849,10 @@ struct ProcessResult { // it matters only against other data. Empty where there is none, or where it was resolved. std::string indexing_ambiguity_operators; - // Serial stills: the indexing ambiguity resolved from the data alone (ReferenceFreeReindex.h), where - // the cell and space group admit one. Without a reference or model it decides the merge; with one, it - // is only measured against the reference's answer (ambiguity_reference_free_applied = false, and - // ambiguity_reindexed is then the number of crystals on which the two disagree). + // Serial stills: the indexing ambiguity resolved (ReferenceFreeReindex.h), where the cell and space + // group admit one, from the data - and, with a reference MTZ, against it as well + // (ambiguity_reference_free_applied = false; ambiguity_reindexed is then the number of crystals the data + // and their cells put into another indexing than the reference alone would). bool ambiguity_reference_free = false; bool ambiguity_reference_free_applied = false; size_t ambiguity_crystals = 0; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 3e33beedd..1575ccea2 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -660,6 +661,19 @@ OtherCellForm LeaveOutOtherCellForms(std::vector &outcomes, return {other_cells.size(), rel.size(), MeanUnitCell(other_cells), below + above}; } +void LogCellMetric(Logger &logger, const CellMetric &metric, const std::vector &ops, + const std::vector &metric_log_odds) { + const auto sd = [&](int j, double unit) { return metric.variance[j] > 0 ? std::sqrt(metric.variance[j]) * unit : NAN; }; + logger.Info("Indexing ambiguity: the crystals' own cells scatter by {:.2f}/{:.2f}/{:.2f}% in shape and " + "{:.2f}/{:.2f}/{:.2f} deg in angle ({} crystals; {} cell(s) shared exactly by several crystals " + "are the prior and say nothing)", sd(0, 100.0), sd(1, 100.0), sd(2, 100.0), + sd(3, 180.0 / std::numbers::pi), sd(4, 180.0 / std::numbers::pi), sd(5, 180.0 / std::numbers::pi), + metric.crystals, metric.priors.size()); + for (size_t i = 0; i < ops.size() && i < metric_log_odds.size(); ++i) + logger.Info("Indexing ambiguity {}: a crystal's own cell favours as indexed by a median log-likelihood " + "of {:.2f}", IndexTriplet(ops[i]), metric_log_odds[i]); +} + OtherCellForm LeaveOutMinorCellForm(std::vector &outcomes, const UnitCell &cell, const IndexingSettings &settings, Logger &logger) { std::vector frame; @@ -1972,25 +1986,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const bool guard_beside = merge_first_part && guard_ahead.valid(); if (guard_beside) rsm->SetWriteBackPerFrameScale(false); - // Stills with a reference or a model: each image was put into the indexing that agrees best with it, - // image by image. Where the crystals' own cells decide an indexing, those cells come first, and the - // reference decides only the crystals they leave open. Before the relabel below, which changes axes. - if (!rsm && ambiguity_reference && result.consensus_cell && experiment_.GetUnitCell() - && experiment_.GetGemmiSpaceGroup()) { - const auto ops = ReindexAmbiguityOperators(*experiment_.GetUnitCell(), experiment_.GetSpaceGroupOrP1()); - if (!ops.empty()) { - const auto anchored = AnchorReferenceReindex(indexer->GetIntegrationOutcome(), *result.consensus_cell, - ops); - for (size_t i = 0; i < ops.size(); ++i) - logger.Info("Indexing ambiguity {}: {:.1f}% of the crystals' own cells fit better as indexed{}", - IndexTriplet(ops[i]), 100.0 * anchored.metric_preference[i], - anchored.metric_preference[i] >= METRIC_DECIDED - ? " - the cells decide it; the reference decides only crystals whose cell does not" - : ""); - logger.Info("Indexing ambiguity: {} crystal(s) the reference had put against their own cell kept " - "in an indexing the cell allows", anchored.overruled); - } - } // Stills: a crystal indexed in a lattice that contains the run's - an axis multiple, or the run's in // another setting - is put on the run's cell before anything is scaled, rather than left out. if (!rsm && result.consensus_cell) { @@ -2002,43 +1997,35 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } const auto initial_sg = experiment_.GetGemmiSpaceGroup(); // Stills in a group whose lattice admits alternative indexings: each crystal was indexed in one of - // them at random. With a reference or a model each image has already been put into the reference's - // (IndexAndRefine::ScaleImage); without one the indexings are made consistent here, from the data - // alone, before anything is scaled. With one, the same is done only to measure it against the - // reference's answer. + // them at random. They are made consistent before anything is scaled, from their intensities against + // each other, against a reference where one is given (which then also fixes the indexing the merge is + // written in; IndexAndRefine::ScaleImage put each image into its best indexing against it already) and + // from their own cells. if (!rsm && initial_sg && result.consensus_cell) { const auto ops = ReindexAmbiguityOperators(*result.consensus_cell, *initial_sg); if (!ops.empty()) { const auto t0 = std::chrono::steady_clock::now(); - const auto rf = ResolveReindexReferenceFree(indexer->GetIntegrationOutcome(), experiment_, - *result.consensus_cell, ops, config_.nthreads, - /*apply=*/!ambiguity_reference); + const auto rf = ResolveIndexingAmbiguity(indexer->GetIntegrationOutcome(), experiment_, + *result.consensus_cell, ops, config_.reference_data, + config_.nthreads); result.ambiguity_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - for (size_t i = 0; i < ops.size(); ++i) - logger.Info("Indexing ambiguity {}: {:.1f}% of the crystals' own cells fit better as indexed{}", - IndexTriplet(ops[i]), 100.0 * rf.metric_preference[i], - rf.metric_preference[i] >= METRIC_DECIDED - ? " - the cells decide it; only crystals whose cell does not are decided by " - "the intensities" : ""); + LogCellMetric(logger, rf.model.metric, ops, rf.metric_log_odds); result.ambiguity_reference_free = true; result.ambiguity_reference_free_applied = !ambiguity_reference; result.ambiguity_crystals = rf.crystals; - result.ambiguity_reindexed = rf.reindexed; + result.ambiguity_reindexed = ambiguity_reference ? rf.reference_overruled : rf.reindexed; result.ambiguity_separation_sigma = rf.median_sigma; result.ambiguity_decided_fraction = rf.decided_fraction; - if (!ambiguity_reference) - logger.Info("Indexing ambiguity resolved from the data: {} of {} crystals reindexed; median " - "separation {:.1f} sigma, {:.0f}% of crystals decided at over 3 sigma (CC with the " - "rest {:.3f} as chosen, {:.3f} otherwise); {} sweeps, {:.2f} s", - rf.reindexed, rf.crystals, rf.median_sigma, 100.0 * rf.decided_fraction, - rf.cc_chosen, rf.cc_other, rf.sweeps, result.ambiguity_time_s); - else - logger.Info("Indexing ambiguity from the data alone, against the reference's: it would put {} " - "of {} crystals into another indexing; median separation {:.1f} sigma, {:.0f}% of " - "crystals decided at over 3 sigma; {:.2f} s", rf.reindexed, rf.crystals, - rf.median_sigma, 100.0 * rf.decided_fraction, result.ambiguity_time_s); - Note(fmt::format("indexing ambiguity from the data: {} of {} crystals {}", rf.reindexed, rf.crystals, - ambiguity_reference ? "disagree with the reference" : "reindexed")); + logger.Info("Indexing ambiguity resolved {}: {} of {} crystals reindexed; median separation {:.1f} " + "sigma, {:.0f}% of crystals decided at over 3 sigma (CC with the rest {:.3f} as chosen, " + "{:.3f} otherwise); {} crystal(s) where the cell outweighed the intensities; {} sweeps, " + "{:.2f} s", ambiguity_reference ? "against the reference and from the data" : "from the data", + rf.reindexed, rf.crystals, rf.median_sigma, 100.0 * rf.decided_fraction, rf.cc_chosen, + rf.cc_other, rf.cell_moved, rf.sweeps, result.ambiguity_time_s); + if (ambiguity_reference) + logger.Info("Indexing ambiguity: {} of {} crystals put by the data and their cells into another " + "indexing than the reference alone would", rf.reference_overruled, rf.crystals); + Note(fmt::format("indexing ambiguity: {} of {} crystals reindexed", rf.reindexed, rf.crystals)); ambiguity_reference = true; } } diff --git a/rugnux/StillsLiveMerge.cpp b/rugnux/StillsLiveMerge.cpp index 46d317dc0..b4a6c68cb 100644 --- a/rugnux/StillsLiveMerge.cpp +++ b/rugnux/StillsLiveMerge.cpp @@ -64,11 +64,19 @@ StillsMergeOutput MergeStills(DiffractionExperiment experiment, std::vector(out.images)); if (flag_ice) FlagIceRings(outcomes, ice_width); - if (const auto ops = AmbiguityOperators(experiment); !ops.empty() && reference_data.empty()) { - const auto rf = ResolveReindexReferenceFree(outcomes, experiment, cell, ops, nthreads, /*apply=*/true); - logger.Info("Indexing ambiguity resolved from the data: {} of {} crystals reindexed", rf.reindexed, - rf.crystals); - out.metric_preference = rf.metric_preference; + if (const auto ops = AmbiguityOperators(experiment); !ops.empty()) { + const auto rf = ResolveIndexingAmbiguity(outcomes, experiment, cell, ops, reference_data, nthreads); + LogCellMetric(logger, rf.model.metric, ops, rf.metric_log_odds); + logger.Info("Indexing ambiguity resolved {}: {} of {} crystals reindexed; median separation {:.1f} sigma, " + "{:.0f}% of crystals decided at over 3 sigma (CC with the rest {:.3f} as chosen, {:.3f} " + "otherwise); {} crystal(s) where the cell outweighed the intensities", + reference_data.empty() ? "from the data" : "against the reference and from the data", + rf.reindexed, rf.crystals, rf.median_sigma, 100.0 * rf.decided_fraction, rf.cc_chosen, + rf.cc_other, rf.cell_moved); + if (!reference_data.empty()) + logger.Info("Indexing ambiguity: {} of {} crystals put by the data and their cells into another " + "indexing than the reference alone would", rf.reference_overruled, rf.crystals); + out.ambiguity = rf.model; out.choice = rf.choice; } out.merge = StillsScaleMerge(experiment, outcomes, cell, run_width, reference_data, /*for_search=*/false, @@ -107,7 +115,7 @@ StillsMergeOutput StillsLiveMerge::Checkpoint(int nthreads) { // post-refined - and kept so. ops = AmbiguityOperators(experiment); out = MergeStills(experiment, outcomes, given_cell, run_width, {}, flag_ice, ice_width, nthreads, logger); - metric_preference = out.metric_preference; + ambiguity = out.ambiguity; choice = out.choice; choice.resize(outcomes.size(), -1); run_width = out.merge.used_width; @@ -142,7 +150,7 @@ void StillsLiveMerge::Prepare(int nthreads) { if (!ops.empty()) { const ReindexAmbiguityResolver resolver(experiment, reference); ParallelFor(static_cast(batch.size()), nthreads, [&](int i) { - choice[pending[i]] = ReindexOnReference(batch[i], resolver, cell, ops, metric_preference); + choice[pending[i]] = ReindexOnReference(batch[i], resolver, ops, ambiguity); }); } const ScaleOnTheFly scaling(experiment, reference); @@ -177,7 +185,7 @@ void StillsLiveMerge::LeaveOutOtherForms() { } void StillsLiveMerge::Finish(StillsMergeOutput &out, int nthreads) { - if (const size_t moved = RelabelRunIndexing(outcomes, choice, experiment, ops, metric_preference); moved > 0) + if (const size_t moved = RelabelRunIndexing(outcomes, choice, experiment, ops, ambiguity.metric); moved > 0) logger.Info("Indexing ambiguity: the run relabelled as a whole, {} crystals moved", moved); // Every crystal is scaled and post-refined: what is left of StillsScaleMerge is its merge. ScaleOnTheFly::RejectCollapsedScales(outcomes); diff --git a/rugnux/StillsLiveMerge.h b/rugnux/StillsLiveMerge.h index cf20ba9e9..43f2f3a46 100644 --- a/rugnux/StillsLiveMerge.h +++ b/rugnux/StillsLiveMerge.h @@ -14,14 +14,15 @@ struct StillsMergeOutput { StillsMergeResult merge; size_t images = 0; // frames with reflections that went in - std::vector metric_preference; // of the ambiguity operators (ResolveReindexReferenceFree) + AmbiguityModel ambiguity; // what the ambiguity resolution left (ResolveIndexingAmbiguity) std::vector choice; // ... and each frame's indexing (ReferenceFreeReindexResult::choice) double seconds = 0.0; }; // The frames of the run's cell form (the given cell's population where one was given), put on the run's // cell where indexed in a lattice containing it, their ice-ring reflections flagged (flag_ice), the indexing -// ambiguity resolved from the data (where no reference_data is given), then StillsScaleMerge. `experiment` carries the cell and space group +// ambiguity resolved (ResolveIndexingAmbiguity, against reference_data too where it is given), then +// StillsScaleMerge. `experiment` carries the cell and space group // merged in; outcomes are changed in place. StillsMergeOutput MergeStills(DiffractionExperiment experiment, std::vector &outcomes, const std::optional &given_cell, @@ -80,7 +81,7 @@ private: std::optional prior; std::optional> error_model; // of the last checkpoint's merge std::vector ops; - std::vector metric_preference; + AmbiguityModel ambiguity; std::vector choice; // each frame's indexing, as ReferenceFreeReindexResult::choice void LeaveOutOtherForms(); diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp index d29fe66f9..b6b0769c9 100644 --- a/tests/ReindexAmbiguityTest.cpp +++ b/tests/ReindexAmbiguityTest.cpp @@ -357,7 +357,7 @@ TEST_CASE("Reindex: a stills set indexed in random hands is made consistent from } } - const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/true); + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); CHECK(result.crystals == n_crystals); // The overall hand is free: either exactly the flipped crystals are reindexed, or exactly the others. int agree = 0; @@ -370,7 +370,7 @@ TEST_CASE("Reindex: a stills set indexed in random hands is made consistent from 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); + const auto again = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); CHECK(again.reindexed == 0); } @@ -435,7 +435,7 @@ TEST_CASE("Reindex: many weak stills with heavy-tailed intensities are made cons } } - const auto result = ResolveReindexReferenceFree(outcomes, x, cell, laws, 4, /*apply=*/false); + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); int agree = 0; for (int c = 0; c < n_crystals; ++c) agree += (result.choice[c] != 0) == flipped[c]; @@ -447,25 +447,77 @@ TEST_CASE("Reindex: a holohedral crystal gives the stills resolution nothing to x.SetUnitCell(Tetragonal()).SpaceGroupNumber(89); // P422 std::vector outcomes(3); const auto laws = ReindexAmbiguityOperators(Tetragonal(), x.GetSpaceGroupOrP1()); - const auto result = ResolveReindexReferenceFree(outcomes, x, Tetragonal(), laws, 1, true); + const auto result = ResolveIndexingAmbiguity(outcomes, x, Tetragonal(), laws, {}, 1); 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. +namespace { + // A crystal's refined cell: `cell` with each length off by a relative error of `scatter`. + UnitCell Scattered(const UnitCell &cell, double scatter, std::mt19937 &rng) { + std::normal_distribution e(0.0, scatter); + return UnitCell{static_cast(cell.a * (1.0 + e(rng))), static_cast(cell.b * (1.0 + e(rng))), + static_cast(cell.c * (1.0 + e(rng))), cell.alpha, cell.beta, cell.gamma}; + } + + // A synthetic P222 structure on hkl up to (12, 12, 6), one Wilson-distributed intensity per reflection + // of the asymmetric unit. + struct Structure { + std::unordered_map intensity; + std::vector> hkls; + }; + Structure MakeStructure(const HKLKeyGenerator &key, std::mt19937 &rng) { + std::exponential_distribution wilson(1.0); + Structure s; + 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) { + s.hkls.push_back({h, k, l}); + s.intensity.try_emplace(key(h, k, l).pack(), wilson(rng)); + } + return s; + } + + // A crystal of `structure` truly in the other indexing (by `law`) where `flipped`, recording 150 random + // reflections with a relative noise of 0.3 - in its own reported indices - with the lattice `latt`. + IntegrationOutcome Crystal(const Structure &structure, const HKLKeyGenerator &key, const UnitCell &cell, + const gemmi::Op &law, bool flipped, const CrystalLattice &latt, std::mt19937 &rng) { + std::normal_distribution noise(0.0, 0.3); + IntegrationOutcome o; + o.latt = latt; + for (int i = 0; i < 150; ++i) { + const auto &t = structure.hkls[rng() % structure.hkls.size()]; + const auto m = flipped ? law.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(structure.intensity.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; + o.reflections.push_back(r); + } + return o; + } +} + +TEST_CASE("Reindex: a crystal's own cell says which axis is which where the run's axes differ beyond its scatter", + "[reindex]") { + // Orthorhombic with a and b 3% apart: a twin law of the metric to within the obliquity tolerance, but each + // crystal's refined cell, scattering by 0.3%, says which axis is which. The intensities say nothing. 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()); + REQUIRE(laws.size() == 1); std::vector outcomes(50); std::mt19937 rng(3); for (auto &o : outcomes) { - o.latt = CrystalLattice(cell); + o.latt = CrystalLattice(Scattered(cell, 0.003, rng)); for (int i = 0; i < 100; ++i) { Reflection r{}; r.h = static_cast(rng() % 15) - 7; r.k = static_cast(rng() % 15) - 7; r.l = static_cast(rng() % 7) + 1; @@ -473,18 +525,18 @@ TEST_CASE("Reindex: an operator the crystals' own cells decide is left to them", 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); + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 1); + REQUIRE(result.metric_log_odds.size() == 1); + CHECK(result.metric_log_odds[0] > 10.0); + CHECK(result.model.metric.variance[0] == Catch::Approx(0.003 * 0.003).epsilon(0.5)); + CHECK(result.crystals == 50); CHECK(result.reindexed == 0); } -TEST_CASE("Reindex: where the cells decide, the crystals whose cell disagrees are decided by the intensities", +TEST_CASE("Reindex: where the cells tell the axes apart, they and the intensities put back the crystals indexed on swapped axes", "[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. + // and so do their intensities. const UnitCell cell{78.0f, 80.4f, 37.0f, 90.0f, 90.0f, 90.0f}; DiffractionExperiment x; x.SetUnitCell(cell).SpaceGroupNumber(16); // P222 @@ -493,6 +545,123 @@ TEST_CASE("Reindex: where the cells decide, the crystals whose cell disagrees ar const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1()); std::mt19937 rng(11); + const Structure structure = MakeStructure(key, rng); + const int n_crystals = 300; + std::vector outcomes; + std::vector swapped(n_crystals, false); + for (int c = 0; c < n_crystals; ++c) { + swapped[c] = c % 12 == 5; + const CrystalLattice latt(Scattered(cell, 0.003, rng)); + outcomes.push_back(Crystal(structure, key, cell, laws[0], swapped[c], + swapped[c] ? ReindexLattice(latt, laws[0]) : latt, rng)); + } + + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); + CHECK(result.crystals == n_crystals); + for (int c = 0; c < n_crystals; ++c) + CHECK((result.choice[c] != 0) == swapped[c]); +} + +TEST_CASE("Reindex: cells that only echo the hand the indexer chose do not lock the crystals in it", "[reindex]") { + // a and b 0.2% apart, well inside the 1% the crystals' refined cells scatter by. Every crystal is truly + // in one of the two indexings at random, but its cell is reported in whichever setting fits the run's + // cell better - as the indexer or a setting convention picks it - so every cell prefers the indexing it + // came in, whatever the truth. The intensities have to decide, and do. + const UnitCell cell{78.0f, 78.16f, 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(17); + const Structure structure = MakeStructure(key, rng); + const int n_crystals = 300; + std::vector outcomes; + std::vector flipped(n_crystals); + for (int c = 0; c < n_crystals; ++c) { + flipped[c] = rng() % 2; + UnitCell u = Scattered(cell, 0.01, rng); + if (std::fabs(u.a - cell.b) + std::fabs(u.b - cell.a) < std::fabs(u.a - cell.a) + std::fabs(u.b - cell.b)) + std::swap(u.a, u.b); + outcomes.push_back(Crystal(structure, key, cell, laws[0], flipped[c], CrystalLattice(u), rng)); + } + + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); + CHECK(result.crystals == n_crystals); + 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) == n_crystals); + CHECK(result.decided_fraction > 0.9); +} + +TEST_CASE("Reindex: with a reference, cells that cannot tell the axes apart leave the indexing to the intensities", "[reindex]") { + // The crystals of the test above, with a reference written in their true indexing: every crystal ends in + // it - the reference fixes the run's indexing as a whole - and the cells, which only echo the indexer's + // choice, put none of them elsewhere. + const UnitCell cell{78.0f, 78.16f, 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(19); + const Structure structure = MakeStructure(key, rng); + std::vector reference; + std::unordered_set seen; + for (const auto &t : structure.hkls) { + if (!seen.insert(key(t[0], t[1], t[2]).pack()).second) + continue; + MergedReflection r; + r.h = t[0]; r.k = t[1]; r.l = t[2]; + r.I = static_cast(structure.intensity.at(key(t[0], t[1], t[2]).pack())); + r.sigma = 1.0f; + 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)); + reference.push_back(r); + } + + const int n_crystals = 200; + std::vector outcomes; + std::vector flipped(n_crystals); + for (int c = 0; c < n_crystals; ++c) { + flipped[c] = rng() % 10 < 7; // most of them: the indexing that moves the fewest is the wrong one + UnitCell u = Scattered(cell, 0.01, rng); + if (std::fabs(u.a - cell.b) + std::fabs(u.b - cell.a) < std::fabs(u.a - cell.a) + std::fabs(u.b - cell.b)) + std::swap(u.a, u.b); + outcomes.push_back(Crystal(structure, key, cell, laws[0], flipped[c], CrystalLattice(u), rng)); + } + const auto as_indexed = outcomes; + const auto moved = [&](int c) { + const auto &r = outcomes[c].reflections[0], &r0 = as_indexed[c].reflections[0]; + return r.h != r0.h || r.k != r0.k || r.l != r0.l; + }; + + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, reference, 4); + CHECK(result.crystals == n_crystals); + CHECK(result.reference_overruled == 0); + for (int c = 0; c < n_crystals; ++c) + CHECK(moved(c) == flipped[c]); +} + +TEST_CASE("Reindex: a monoclinic crystal on a pseudo-tetragonal lattice is relabelled as a whole only by a symmetry of its merge", + "[reindex]") { + // P21 with a and b equal: four indexings. Relabelling a whole merge by a diagonal operator would put the + // two-fold axis on a, so it is no relabelling of a P21 merge - taken as one, it splits crystals that differ + // by the two-fold into different indexings. Most crystals are truly in a diagonal indexing here, so the + // indexing that moves the fewest of them is such an operator. + const UnitCell cell{67.0f, 67.0f, 71.0f, 90.0f, 90.0f, 90.0f}; + DiffractionExperiment x; + x.SetUnitCell(cell).SpaceGroupNumber(4); // P21 + const auto laws = ReindexAmbiguityOperators(cell, x.GetSpaceGroupOrP1()); + REQUIRE(laws.size() == 3); + std::vector indexings{gemmi::Op::identity()}; + indexings.insert(indexings.end(), laws.begin(), laws.end()); + + const HKLKeyGenerator key(true, x.GetSpaceGroupOrP1()); + std::mt19937 rng(23); std::exponential_distribution wilson(1.0); std::normal_distribution noise(0.0, 0.3); std::unordered_map structure; @@ -507,17 +676,16 @@ TEST_CASE("Reindex: where the cells decide, the crystals whose cell disagrees ar const int n_crystals = 300; std::vector outcomes(n_crystals); - std::vector swapped(n_crystals, false); + std::vector truth(n_crystals); 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); + truth[c] = c % 10 < 7 ? 2 : static_cast(rng() % 4); + outcomes[c].latt = 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]}}; + const auto m = indexings[truth[c]].apply_to_hkl({{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.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(structure.at(key(t[0], t[1], t[2]).pack()) * (1.0 + noise(rng))); r.sigma = 0.1f; r.partiality = 1.0f; @@ -526,100 +694,16 @@ TEST_CASE("Reindex: where the cells decide, the crystals whose cell disagrees ar } } - 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 wilson(1.0); - std::normal_distribution noise(0.0, 0.2); - std::unordered_map structure; - std::vector> 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 &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)); + const auto result = ResolveIndexingAmbiguity(outcomes, x, cell, laws, {}, 4); + // Consistent: a general reflection, put into each crystal's indexing as indexed and then as chosen, lands + // on the same reflection of the merge for every crystal. + const std::array general{2, 5, 11}; + const auto merged = [&](int c) { + const auto m = indexings[result.choice[c]].apply_to_hkl(indexings[truth[c]].apply_to_hkl(general)); + return key(m[0], m[1], m[2]).pack(); }; - std::vector reference; - std::unordered_set 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(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 outcomes(n_crystals); - std::vector 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(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); - } + int consistent = 0; 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(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)); - } + consistent += merged(c) == merged(0); + CHECK(consistent == n_crystals); }