Jungfraujoch/image_analysis/indexing/AnalyzeIndexing.h

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

#pragma once

#include "../../common/CrystalLattice.h"
#include "../../common/DiffractionExperiment.h"
#include "../../common/JFJochMessages.h"

// Minimum fraction of a frame's in-resolution spots that must lie on a candidate lattice for the
// frame to be that crystal's. See the frame gate in AnalyzeIndexing, which is where it is applied.
constexpr float LATTICE_MIN_INDEXED_FRACTION = 0.20f;

// Tally one image's spots by their rank in its intensity-ordered spot list: how many images had a spot
// at that rank at all (`counted`) and on how many of them it lay on the lattice (`indexed`). Ice spots
// are skipped, as they are in the frame gate. Both are added to, and their length bounds the ranks
// considered. Counts rather than weights so that the tally is exact whatever order the images are
// summed in, which is what makes the budget below independent of the thread schedule.
void AddSpotBudgetEvidence(const std::vector<SpotToSave> &spots, bool index_ice_rings,
                           std::vector<int64_t> &indexed, std::vector<int64_t> &counted);

// How far the fall from the peak must exceed the counting noise of the spots for the peak to be one.
//
// Under the null - the spots lie on the lattice at the same rate at every depth - the running sum
// below is a driftless random walk in the counted spots: each is worth 1 - g with probability g and
// -g otherwise, so its step has mean zero and variance g(1-g). The MAXIMUM of such a walk is positive
// whatever the data, so an argmax taken on its own cuts every dataset, including one with nothing to
// cut. What is acted on is the FALL from the peak to the end of the list, which is the maximum of the
// same walk read backwards from the end, and the reflection principle gives that maximum's null law
// exactly: P(fall > z sqrt(g(1-g)T)) = 2(1 - Phi(z)) over T counted spots in all. The search over the
// ranks is therefore already paid for and no further multiple-comparison correction is due. z is set
// for one false cut in a thousand measurements, which over a corpus the size of a rotation test set
// (tens of crystals, a measurement per pass) expects none at all.
constexpr float SPOT_BUDGET_SIGNIFICANCE_Z = 3.29f;   // 2(1 - Phi(z)) = 0.001

// The spot budget those tallies support: the rank at which indexed - LATTICE_MIN_INDEXED_FRACTION *
// counted, summed over the ranks down to it, peaks. Each spot that lies on the lattice is worth
// 1 - LATTICE_MIN_INDEXED_FRACTION and each one that does not costs LATTICE_MIN_INDEXED_FRACTION - the
// same weighing the frame gate applies to a spot list as a whole - so the sum rises exactly while the
// spots at that depth are on the lattice more often than the gate's floor. Deeper than the peak they
// are not: they are no longer this crystal's reflections, and they can only push a frame towards
// rejection while adding nothing the lattice recognises.
//
// Zero - keep the whole list - when the fall from that peak to the end of the list is no larger than
// the counting noise above, which is the case whenever the spots go on lying on the lattice at the
// same rate all the way down, and the case a bare argmax gets wrong.
int64_t SpotBudgetFromEvidence(const std::vector<int64_t> &indexed, const std::vector<int64_t> &counted);

bool AnalyzeIndexing(DataMessage &message,
                     const DiffractionExperiment &experiment,
                     const CrystalLattice &latt,
                     const std::vector<CrystalLattice> &extra_lattices = {});