Jungfraujoch/image_analysis/scale_merge/ReferenceFreeReindex.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <cmath>
#include <cstdint>
#include <vector>

#include "gemmi/symmetry.hpp"

#include "../../common/DiffractionExperiment.h"
#include "../IntegrationOutcome.h"

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;

// 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<double> metric_preference;   // per operator: share of the images whose refined cell is
                                             // closer to `cell` as indexed than after the operator
    std::vector<uint8_t> allowed;            // allowed[o * (ops.size() + 1) + g]: indexing g open to image o
};
MetricAnchors AnchorOnCells(const std::vector<IntegrationOutcome> &outcomes, const UnitCell &cell,
                            const std::vector<gemmi::Op> &ops);

// 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<double> 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
};
ReferenceAnchorResult AnchorReferenceReindex(std::vector<IntegrationOutcome> &outcomes, const UnitCell &cell,
                                             const std::vector<gemmi::Op> &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 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<double> metric_preference;   // per operator: share of crystals whose refined cell is
                                             // closer to `cell` as indexed than after the operator
    int sweeps = 0;                // passes over the crystals until (almost) none changed (the start kept)
    std::vector<int> 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<IntegrationOutcome> &outcomes,
                                                       const DiffractionExperiment &x,
                                                       const UnitCell &cell,
                                                       const std::vector<gemmi::Op> &ops,
                                                       size_t nthreads, bool apply);

// 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<gemmi::Op> &ops, const std::vector<double> &metric_preference);

// 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.
size_t RelabelRunIndexing(std::vector<IntegrationOutcome> &outcomes, std::vector<int> &choice,
                          const DiffractionExperiment &x, const std::vector<gemmi::Op> &ops,
                          const std::vector<double> &metric_preference);