Jungfraujoch/common/IndexingSettings.cpp

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

#include "IndexingSettings.h"
#include "JFJochException.h"
#include "CUDAWrapper.h"

#include <cmath>

#define check_max(param, val, max) if ((val) > (max)) throw JFJochException(JFJochExceptionCategory::InputParameterAboveMax, param)
#define check_min(param, val, min) if ((val) < (min)) throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, param)
#define check_finite(param, val) if (!std::isfinite(val)) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, param)

IndexingSettings::IndexingSettings() {
    if (get_gpu_count() > 0)
        algorithm = IndexingAlgorithmEnum::FFBIDX;
    else
        algorithm = IndexingAlgorithmEnum::None;
}

IndexingSettings &IndexingSettings::ViableCellMinSpots(int64_t input) {
    // Four, not the default six: six is the stills frame gate's number, and a rotation run is no longer
    // accepted or refused on it (the pooled spot test in the first pass decides that), so it has no
    // business being a hard minimum. Four is where the number stops being a convention and becomes a
    // fact - three spots fit any lattice at all (see MIN_SPOT_COUNT), and the cell and orientation
    // refinements this same setting bounds need more observations than parameters.
    check_min("ViableCellMinSpots", input, 4);
    viable_cell_min_spots = input;
    return *this;
}

int64_t IndexingSettings::GetViableCellMinSpots() const {
    return viable_cell_min_spots;
}

IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) {
    switch (input) {
        case IndexingAlgorithmEnum::Auto:
        case IndexingAlgorithmEnum::FFBIDX:
        case IndexingAlgorithmEnum::FFT:
        case IndexingAlgorithmEnum::FFTW:
        case IndexingAlgorithmEnum::None:
            algorithm = input;
            break;
        default:
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  "Invalid value for indexing algorithm enum parameter");
    }
    return *this;
}

// The accepted range, not the default (which stays 500 A - see IndexingSettings.h). The FFT can
// only recover a basis vector up to this length, so the ceiling is exactly the longest cell the
// indexer can ever find; at 500 it excluded 1091 of the PDB's 206950 X-ray entries (0.53%) outright.
// 1200 A leaves 8. Nothing pays for the wider range: the histogram is sized from the VALUE in use,
// and only a caller that asks for more - a given cell that needs it, or the long-axis rescue - gets
// a longer transform.
IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) {
    check_finite("FFT indexing max unit cell (A)", input);
    check_min("FFT indexing max unit cell (A)", input, 50);
    check_max("FFT indexing max unit cell (A)", input, fft_max_unit_cell_limit_A);
    fft_max_unit_cell_A = input;
    return *this;
}

IndexingSettings &IndexingSettings::FFT_MinUnitCell_A(float input) {
    check_finite("FFT indexing min unit cell (A)", input);
    check_min("FFT indexing min unit cell (A)", input, fft_min_unit_cell_limit_A);
    check_max("FFT indexing min unit cell (A)", input, 40);
    fft_min_unit_cell_A = input;
    return *this;
}

IndexingSettings & IndexingSettings::FFT_MaxAngle_deg(float input) {
    check_finite("FFT indexing max angle (deg)", input);
    check_min("FFT indexing max angle (deg)", input, 0);
    check_max("FFT indexing max angle (deg)", input, 180);
    fft_max_angle_deg = input;
    return *this;
}

IndexingSettings & IndexingSettings::FFT_MinAngle_deg(float input) {
    check_finite("FFT indexing min angle (deg)", input);
    check_min("FFT indexing min angle (deg)", input, 0);
    check_max("FFT indexing min angle (deg)", input, 180);
    fft_min_angle_deg = input;
    return *this;
}

float IndexingSettings::GetFFT_MinAngle_deg() const {
    return fft_min_angle_deg;
}

float IndexingSettings::GetFFT_MaxAngle_deg() const {
    return fft_max_angle_deg;
}

IndexingSettings &IndexingSettings::FFT_NumVectors(int64_t input) {
    check_min("FFT indexing number of search vectors", input, 128);
    fft_num_vectors = input;
    return *this;
}

IndexingSettings &IndexingSettings::FFT_HighResolution_A(float input) {
    check_finite("FFT indexing high resolution (A)", input);
    check_min("FFT indexing high resolution (A)", input, 0.5);
    check_max("FFT indexing high resolution (A)", input, 6.0);

    fft_high_resolution_A = input;
    return *this;
}

IndexingAlgorithmEnum IndexingSettings::GetAlgorithm() const {
    return algorithm;
}

float IndexingSettings::GetFFT_MaxUnitCell_A() const {
    return fft_max_unit_cell_A;
}

float IndexingSettings::GetFFT_MinUnitCell_A() const {
    return fft_min_unit_cell_A;
}

int64_t IndexingSettings::GetFFT_NumVectors() const {
    return fft_num_vectors;
}

float IndexingSettings::GetFFT_HighResolution_A() const {
    return fft_high_resolution_A;
}

IndexingSettings &IndexingSettings::Tolerance(float input) {
    check_min("Indexing tolerance", input, 0.0);
    check_max("Indexing tolerance", input, 0.5);
    indexing_tolerance = input;
    return *this;
}

float IndexingSettings::GetTolerance() const {
    return indexing_tolerance;
}

int64_t IndexingSettings::GetIndexingThreads() const {
    return indexing_threads;
}

IndexingSettings &IndexingSettings::IndexingThreads(int64_t input) {
    check_min("Indexing thread count", input, 1);
    check_max("Indexing thread count", input, 64);
    indexing_threads = input;
    return *this;
}

int64_t IndexingSettings::GetRefineThreads() const {
    return refine_threads;
}

IndexingSettings &IndexingSettings::RefineThreads(int64_t input) {
    check_min("Candidate-cell refinement thread count", input, 1);
    check_max("Candidate-cell refinement thread count", input, MAX_REFINE_THREADS);
    refine_threads = input;
    return *this;
}

IndexingSettings &IndexingSettings::UnitCellDistTolerance(float input) {
    check_min("Relative unit cell distance tolerance vs. reference", input, 0.0001);
    check_max("Relative unit cell distance tolerance vs. reference", input, 0.2001);
    unit_cell_dist_tolerance_vs_reference = input;
    return *this;
}

float IndexingSettings::GetUnitCellDistTolerance() const {
    return unit_cell_dist_tolerance_vs_reference;
}

float IndexingSettings::GetUnitCellAngleTolerance_deg() const {
    return unit_cell_angle_tolerance_deg;
}

GeomRefinementAlgorithmEnum IndexingSettings::GetGeomRefinementAlgorithm() const {
    return refinement;
}

IndexingSettings &IndexingSettings::GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum input) {
    refinement = input;
    return *this;
}

IndexingSettings & IndexingSettings::IndexIceRings(bool input) {
    index_ice_rings = input;
    return *this;
}

IndexingSettings & IndexingSettings::RotationIndexing(bool input) {
    enable_rotation_indexing = input;
    return*this;
}

IndexingSettings & IndexingSettings::RotationIndexingMinAngularRange_deg(float input) {
    check_finite("Rotation indexing minimum angular range (deg.)", input);
    check_min("Rotation indexing minimum angular range (deg.)", input, 1.0);
    rotation_indexing_min_angular_range_deg = input;
    return *this;
}

IndexingSettings & IndexingSettings::RotationIndexingAngularStride_deg(float input) {
    check_finite("Rotation indexing angular stride (deg.)", input);
    check_min("Rotation indexing angular stride (deg.)", input, 0.0);
    rotation_indexing_angular_stride_deg = input;
    return *this;
}

bool IndexingSettings::GetRotationIndexing() const {
    return enable_rotation_indexing;
}

float IndexingSettings::GetRotationIndexingMinAngularRange_deg() const {
    return rotation_indexing_min_angular_range_deg;
}

float IndexingSettings::GetRotationIndexingAngularStride_deg() const {
    return rotation_indexing_angular_stride_deg;
}

bool IndexingSettings::GetIndexIceRings() const {
    return index_ice_rings;
}

bool IndexingSettings::GetBlockingBehavior() const {
    return blocking_behavior;
}

IndexingSettings &IndexingSettings::BlockingBehavior(bool input) {
    blocking_behavior = input;
    return *this;
}

int64_t IndexingSettings::GetMaxExtraLattices() const {
    return max_extra_lattices;
}

IndexingSettings &IndexingSettings::MaxExtraLattices(int64_t input) {
    check_min("Max extra lattices", input, 0);
    check_max("Max extra lattices", input, 10);
    max_extra_lattices = input;
    return *this;
}