twitcher 0.6.9

//! Startup checks for setting combinations that are contradictory, unreachable, or silently
//! ignored.
//!
//! None of these ever fail a run: an invocation may take hours, so it must not die on a
//! heuristic. Every finding is emitted as a `tracing::warn!`.
//!
//! The checks are split into pure functions returning the messages and thin wrappers that log
//! them, so that the logic is unit-testable without capturing logs.

use generic_a_star::cost::AStarCost;
use lib_tsalign::config::TemplateSwitchConfig;
use lib_tsalign::costs::U64Cost;
use lib_tsalign::costs::cost_function::CostFunction;
use lib_tsalign::costs::gap_affine::GapAffineAlignmentCostTable;

use compact_genome::implementation::alphabets::dna_alphabet_or_n::DnaAlphabetOrN;

use crate::common::aligner::cli::{CliAlignmentArgs, MLSSelector};
use crate::common::cluster_settings::{ClusterStrategy, ClusteringSettings};

type Costs = TemplateSwitchConfig<DnaAlphabetOrN, U64Cost>;
type EditCosts = GapAffineAlignmentCostTable<DnaAlphabetOrN, U64Cost>;

fn infinite() -> U64Cost {
    U64Cost::from_primitive(u64::MAX)
}

/// Emits a warning for every questionable combination of alignment settings and costs.
pub fn warn_about_alignment_settings(args: &CliAlignmentArgs, costs: &Costs) {
    for message in alignment_setting_warnings(args, costs) {
        tracing::warn!("{message}");
    }
}

/// Emits a warning for every questionable combination of clustering settings.
pub fn warn_about_clustering_settings(settings: &ClusteringSettings) {
    for message in clustering_setting_warnings(settings) {
        tracing::warn!("{message}");
    }
}

/// Renders a finite bound of a cost table, or `unbounded` if there is none.
fn render_bound(bound: Option<isize>) -> String {
    bound.map_or_else(|| "unbounded".to_owned(), |bound| bound.to_string())
}

/// The finite bounds of an offset cost table, where `None` means unbounded in that direction.
///
/// Returns `None` if the table is infinite everywhere, i.e. if it permits no jump at all.
fn offset_bounds(costs: &CostFunction<isize, U64Cost>) -> Option<(Option<isize>, Option<isize>)> {
    let lower = costs.minimum_finite_input()?;
    Some((
        (lower != isize::MIN).then_some(lower),
        costs.maximum_finite_input(),
    ))
}

/// Returns true if the cost table is zero everywhere, i.e. if it imposes no preference and no
/// limit on the values it prices.
fn is_free<Input: Clone>(costs: &CostFunction<Input, U64Cost>) -> bool {
    let points: Vec<(Input, U64Cost)> = costs.clone().into();
    points
        .iter()
        .all(|(_, cost)| *cost == U64Cost::from_primitive(0))
}

/// Returns true if any base or gap in the table costs something to open.
fn charges_gap_open(costs: &EditCosts) -> bool {
    costs.max_gap_open_cost() != U64Cost::from_primitive(0)
}

/// Which template switch geometries the cost model leaves open, and whether any of them can be
/// realised at all.
#[expect(
    clippy::struct_excessive_bools,
    reason = "independent predicates about the cost model, not a state machine"
)]
struct Geometries {
    /// Some geometry with a different ancestor and descendant has a finite base cost.
    inter_strand: bool,
    /// Some geometry with the same ancestor and descendant has a finite base cost.
    intra_strand: bool,
    /// Some forward geometry has a finite base cost.
    forward: bool,
    /// The shortest inner segment the `Length` table prices, if any.
    min_length: Option<usize>,
    /// A template switch can be entered: some geometry has both a finite base cost and an offset
    /// table that is finite somewhere, and inner segments have a finite length cost.
    possible: bool,
}

impl Geometries {
    fn of(costs: &Costs) -> Self {
        let infinite = infinite();
        let base = &costs.base_cost;

        let inter_strand = [base.rqf, base.rqr, base.qrf, base.qrr]
            .iter()
            .any(|cost| *cost < infinite);
        let intra_strand = [base.rrf, base.rrr, base.qqf, base.qqr]
            .iter()
            .any(|cost| *cost < infinite);
        let min_length = costs.length_costs.minimum_finite_input();

        let inter_usable =
            inter_strand && costs.rq_qr_offset_costs.minimum_finite_input().is_some();
        let intra_usable =
            intra_strand && costs.rr_qq_offset_costs.minimum_finite_input().is_some();

        Self {
            inter_strand,
            intra_strand,
            forward: [base.rrf, base.rqf, base.qrf, base.qqf]
                .iter()
                .any(|cost| *cost < infinite),
            min_length,
            possible: min_length.is_some() && (inter_usable || intra_usable),
        }
    }
}

/// Warnings about alignment settings that contradict each other or the loaded cost model.
///
/// `padding` bounds where a template switch may jump to, because it is the only ancestor
/// sequence the aligner ever sees. Cost tables that price jumps, inner lengths or geometries
/// beyond what `padding` can reach are therefore partly dead weight.
pub fn alignment_setting_warnings(args: &CliAlignmentArgs, costs: &Costs) -> Vec<String> {
    let mut warnings = Vec::new();
    let geometries = Geometries::of(costs);

    // A cost model under which no template switch can ever be found. `--no-ts` does the same
    // thing but skips the machinery, so it is strictly faster.
    if !args.no_ts && !geometries.possible {
        let reason = if geometries.min_length.is_none() {
            "the `Length` table is infinite everywhere"
        } else if !geometries.inter_strand && !geometries.intra_strand {
            "all eight base costs are infinite"
        } else {
            "every offset table for a geometry with a finite base cost is infinite everywhere"
        };
        warnings.push(format!(
            "No template switch can ever be found with these costs, because {reason}. Pass --no-ts to skip the template switch alignment entirely."
        ));
    }

    warnings.extend(reachability_warnings(args, costs, &geometries));

    // `preprocess-filter` is implemented for reverse template switches only, so with forward
    // geometries enabled it prunes nothing at all.
    if geometries.possible
        && geometries.forward
        && !args.no_ts
        && !args.use_fpa
        && matches!(
            args.min_length_strategy.unwrap_or_default(),
            MLSSelector::PreprocessFilter
        )
    {
        warnings.push(
            "Forward template switches are enabled, but the minimum length strategy `preprocess-filter` only prunes reverse template switches. Consider --min-length-strategy lookahead."
                .to_owned(),
        );
    }

    if args.use_fpa {
        warnings.extend(fpa_warnings(costs, geometries.intra_strand));
    }

    // Flanks are documented as unsupported, and they additionally disable the equal-cost
    // extension of the alignment beyond its range.
    if costs.left_flank_length != 0 || costs.right_flank_length != 0 {
        warnings.push(format!(
            "The cost model sets left_flank_length = {} and right_flank_length = {}. Flanks are not supported and additionally disable extending the alignment beyond its range; set both to 0.",
            costs.left_flank_length, costs.right_flank_length,
        ));
    }

    if args.no_ts {
        let flags: Vec<_> = [
            ("--costs", args.costs.is_some()),
            ("--min-length-strategy", args.min_length_strategy.is_some()),
            ("--allow-mixed-descendants", args.allow_mixed_descendants),
        ]
        .into_iter()
        .filter_map(|(flag, given)| given.then_some(flag))
        .collect();

        if !flags.is_empty() {
            warnings.push(format!(
                "--no-ts disables template switch alignment, so {} has no effect.",
                flags.join(", "),
            ));
        }
    }

    warnings
}

/// Warnings about cost tables that price template switches the padded window can never contain.
fn reachability_warnings(
    args: &CliAlignmentArgs,
    costs: &Costs,
    geometries: &Geometries,
) -> Vec<String> {
    let mut warnings = Vec::new();
    let padding = isize::try_from(args.padding).unwrap_or(isize::MAX);

    for (name, table, reachable) in [
        (
            "RQQROffset",
            &costs.rq_qr_offset_costs,
            geometries.inter_strand && !args.no_ts,
        ),
        (
            "RRQQOffset",
            &costs.rr_qq_offset_costs,
            geometries.intra_strand && !args.no_ts && !args.use_fpa,
        ),
    ] {
        if !reachable {
            continue;
        }
        let Some((lower, upper)) = offset_bounds(table) else {
            continue;
        };
        if lower.is_none_or(|lower| lower < -padding) || upper.is_none_or(|upper| upper > padding) {
            warnings.push(format!(
                "The `{name}` table prices template switch jumps in [{}, {}], but --padding {} limits reachable jumps to ±{}. Jumps beyond that can never be found; raise --padding or narrow the table.",
                render_bound(lower),
                render_bound(upper),
                args.padding,
                args.padding,
            ));
        }
    }

    // The inner segment is copied from the ancestor, which spans the cluster plus the padding on
    // both sides. Pricing lengths far beyond that only widens the search.
    if !args.no_ts
        && let Some(max_length) = costs.length_costs.maximum_finite_input()
        && max_length > args.padding.saturating_mul(2)
    {
        warnings.push(format!(
            "The `Length` table permits template switches of up to {max_length} bases, but with --padding {} only some more than {} bases of template are available. Probably you wnat to lower the last `Length` breakpoint or raise --padding.",
            args.padding,
            args.padding.saturating_mul(2),
        ));
    }

    warnings
}

/// A single warning naming the configured parts of the cost model that the four-point aligner
/// silently disregards.
fn fpa_warnings(costs: &Costs, intra_strand: bool) -> Option<String> {
    let mut ignored = Vec::new();

    if charges_gap_open(&costs.primary_edit_costs)
        || charges_gap_open(&costs.secondary_forward_edit_costs)
        || charges_gap_open(&costs.secondary_reverse_edit_costs)
    {
        ignored.push("the gap open costs (the four-point aligner is not gap-affine)");
    }
    if !is_free(&costs.length_costs) {
        ignored.push("the `Length` table");
    }
    if !is_free(&costs.length_difference_costs) {
        ignored.push("the `LengthDifference` table");
    }
    if !is_free(&costs.forward_anti_descendant_gap_costs) {
        ignored.push("the `ForwardAntiDescendantGap` table");
    }
    if !is_free(&costs.reverse_anti_descendant_gap_costs) {
        ignored.push("the `ReverseAntiDescendantGap` table");
    }
    if intra_strand {
        ignored.push(
            "the intra-strand base costs `rrf`/`rrr`/`qqf`/`qqr` and the `RRQQOffset` table (the four-point aligner only finds inter-strand template switches)",
        );
    }

    (!ignored.is_empty()).then(|| {
        format!(
            "--fpa ignores parts of the cost model that are configured here: {}.",
            ignored.join(", "),
        )
    })
}

/// Warnings about clustering settings whose gates are either inert or so tight that almost
/// nothing can pass them.
pub fn clustering_setting_warnings(settings: &ClusteringSettings) -> Vec<String> {
    let mut warnings = Vec::new();

    // A chain of events spaced `max_gap` apart approaches a density of `1 / (max_gap + 1)` from
    // above, so no legal cluster can ever fall below that.
    #[expect(clippy::cast_precision_loss)]
    let minimum_reachable_density = 1.0 / (settings.max_gap as f64 + 1.0);
    if settings.min_density != 0.0 {
        warnings.push("--cluster-min-density is not recommended to be used. See https://version.helsinki.fi/kraujasp/twitcher/-/work_items/76 for more info.".to_string());
        if settings.min_density <= minimum_reachable_density {
            warnings.push(format!(
            "--cluster-min-density {} can never reject a cluster, because --cluster-max-gap {} already forces a density above {minimum_reachable_density:.4}.",
            settings.min_density, settings.max_gap,
        ));
        } else if settings.strategy == ClusterStrategy::Legacy && settings.min_density > 1.0 {
            warnings.push(format!(
            "--cluster-min-density {} exceeds the density of a run of single nucleotide variants under --cluster-strategy legacy, so almost nothing will be clustered.",
            settings.min_density,
        ));
        }
    }

    if settings.strategy == ClusterStrategy::Legacy && settings.min_records == Some(1) {
        warnings.push(
            "--cluster-min-records 1 under --cluster-strategy legacy turns every single variant into a cluster to be realigned, which is slow and mostly produces noise."
                .to_owned(),
        );
    }

    warnings
}

#[cfg(test)]
mod tests {
    use super::*;

    const DEFAULT_TSA: &str = include_str!("../../default_costs.tsa");

    fn costs_from(tsa: &str) -> Costs {
        TemplateSwitchConfig::read_plain(tsa.as_bytes()).unwrap()
    }

    fn default_costs() -> Costs {
        costs_from(DEFAULT_TSA)
    }

    /// The default cost file with `from` replaced by `to`, to vary a single table per test.
    fn costs_with(from: &str, to: &str) -> Costs {
        assert!(
            DEFAULT_TSA.contains(from),
            "cost file does not contain {from:?}"
        );
        costs_from(&DEFAULT_TSA.replace(from, to))
    }

    fn assert_mentions(warnings: &[String], needle: &str) {
        assert!(
            warnings.iter().any(|w| w.contains(needle)),
            "expected a warning mentioning {needle:?}, got {warnings:#?}"
        );
    }

    /// The bundled costs must be self-consistent with the default padding, otherwise every plain
    /// invocation of twitcher emits warnings.
    #[test]
    fn default_settings_and_default_costs_are_silent() {
        let warnings = alignment_setting_warnings(&CliAlignmentArgs::default(), &default_costs());
        assert!(warnings.is_empty(), "{warnings:#?}");
        assert!(clustering_setting_warnings(&ClusteringSettings::default()).is_empty());
    }

    #[test]
    fn offset_window_wider_than_padding_warns() {
        let args = CliAlignmentArgs {
            padding: 30,
            ..Default::default()
        };
        assert_mentions(
            &alignment_setting_warnings(&args, &default_costs()),
            "`RQQROffset` table prices template switch jumps in [-200, 200]",
        );
    }

    /// A table that is finite for every jump is unbounded, not merely too wide.
    #[test]
    fn unbounded_offset_window_warns() {
        let costs = costs_with(
            "RQQROffset\n -inf -200 201\n  inf    0 inf",
            "RQQROffset\n -inf\n    0",
        );
        assert_mentions(
            &alignment_setting_warnings(&CliAlignmentArgs::default(), &costs),
            "jumps in [unbounded, unbounded]",
        );
    }

    /// The intra-strand table is only reachable if some intra-strand geometry has a finite base
    /// cost, which the defaults do not have.
    #[test]
    fn unreachable_offset_table_is_not_warned_about() {
        let args = CliAlignmentArgs {
            padding: 30,
            ..Default::default()
        };
        let warnings = alignment_setting_warnings(&args, &default_costs());
        assert!(
            !warnings.iter().any(|w| w.contains("RRQQOffset")),
            "{warnings:#?}"
        );
    }

    #[test]
    fn inner_length_beyond_the_available_template_warns() {
        let args = CliAlignmentArgs {
            padding: 30,
            ..Default::default()
        };
        assert_mentions(
            &alignment_setting_warnings(&args, &default_costs()),
            "permits template switches of up to 400 bases",
        );
    }

    /// An unbounded `Length` table is taken as deliberate and not warned about.
    #[test]
    fn unbounded_inner_length_is_not_warned_about() {
        let costs = costs_with(
            "Length\n   0  5  6 7 8 401\n inf 15 10 5 0 inf",
            "Length\n   0 5\n inf 0",
        );
        let warnings = alignment_setting_warnings(&CliAlignmentArgs::default(), &costs);
        assert!(
            !warnings
                .iter()
                .any(|w| w.contains("permits template switches")),
            "{warnings:#?}"
        );
    }

    #[test]
    fn cost_model_without_any_usable_geometry_warns() {
        let costs = costs_with(
            "rqr_cost = 2\nqrr_cost = 2",
            "rqr_cost = inf\nqrr_cost = inf",
        );
        assert_mentions(
            &alignment_setting_warnings(&CliAlignmentArgs::default(), &costs),
            "all eight base costs are infinite",
        );
    }

    #[test]
    fn cost_model_without_any_finite_length_warns() {
        let costs = costs_with(
            "Length\n   0  5  6 7 8 401\n inf 15 10 5 0 inf",
            "Length\n   0\n inf",
        );
        assert_mentions(
            &alignment_setting_warnings(&CliAlignmentArgs::default(), &costs),
            "the `Length` table is infinite everywhere",
        );
    }

    #[test]
    fn no_template_switch_possible_is_not_warned_about_with_no_ts() {
        let costs = costs_with(
            "rqr_cost = 2\nqrr_cost = 2",
            "rqr_cost = inf\nqrr_cost = inf",
        );
        let args = CliAlignmentArgs {
            no_ts: true,
            ..Default::default()
        };
        let warnings = alignment_setting_warnings(&args, &costs);
        assert!(
            !warnings.iter().any(|w| w.contains("No template switch")),
            "{warnings:#?}"
        );
    }

    #[test]
    fn forward_geometry_with_preprocess_filter_warns() {
        let costs = costs_with("rqf_cost = inf", "rqf_cost = 8");
        assert_mentions(
            &alignment_setting_warnings(&CliAlignmentArgs::default(), &costs),
            "only prunes reverse template switches",
        );
    }

    #[test]
    fn forward_geometry_with_lookahead_is_silent() {
        let costs = costs_with("rqf_cost = inf", "rqf_cost = 8");
        let args = CliAlignmentArgs {
            min_length_strategy: Some(MLSSelector::Lookahead),
            ..Default::default()
        };
        let warnings = alignment_setting_warnings(&args, &costs);
        assert!(
            !warnings.iter().any(|w| w.contains("prunes reverse")),
            "{warnings:#?}"
        );
    }

    #[test]
    fn fpa_reports_the_cost_tables_it_ignores() {
        let args = CliAlignmentArgs {
            use_fpa: true,
            ..Default::default()
        };
        let warnings = alignment_setting_warnings(&args, &default_costs());
        assert_mentions(&warnings, "--fpa ignores parts of the cost model");
        assert_mentions(&warnings, "gap open costs");
        assert_mentions(&warnings, "the `Length` table");
        // Free everywhere in the defaults, so it must not be listed.
        assert!(
            !warnings
                .iter()
                .any(|w| w.contains("`ReverseAntiDescendantGap`")),
            "{warnings:#?}"
        );
    }

    #[test]
    fn fpa_reports_unsupported_intra_strand_geometries() {
        let costs = costs_with("rrr_cost = inf", "rrr_cost = 3");
        let args = CliAlignmentArgs {
            use_fpa: true,
            ..Default::default()
        };
        assert_mentions(
            &alignment_setting_warnings(&args, &costs),
            "only finds inter-strand template switches",
        );
    }

    #[test]
    fn nonzero_flanks_warn() {
        let costs = costs_with("left_flank_length = 0", "left_flank_length = 5");
        assert_mentions(
            &alignment_setting_warnings(&CliAlignmentArgs::default(), &costs),
            "Flanks are not supported",
        );
    }

    #[test]
    fn no_ts_with_template_switch_only_flags_warns() {
        let args = CliAlignmentArgs {
            no_ts: true,
            costs: Some("some_costs.tsa".to_owned()),
            allow_mixed_descendants: true,
            ..Default::default()
        };
        assert_mentions(
            &alignment_setting_warnings(&args, &default_costs()),
            "--costs, --allow-mixed-descendants has no effect",
        );
    }

    #[test]
    fn density_gate_that_can_never_reject_warns() {
        let settings = ClusteringSettings {
            min_density: 0.01,
            max_gap: 20,
            ..Default::default()
        };
        assert_mentions(
            &clustering_setting_warnings(&settings),
            "can never reject a cluster",
        );
    }

    #[test]
    fn density_gate_above_one_warns_under_legacy() {
        let settings = ClusteringSettings {
            min_density: 1.5,
            ..Default::default()
        };
        assert_mentions(
            &clustering_setting_warnings(&settings),
            "almost nothing will be clustered",
        );
    }

    #[test]
    fn single_record_clusters_warn_under_legacy() {
        let settings = ClusteringSettings {
            min_records: Some(1),
            ..Default::default()
        };
        assert_mentions(
            &clustering_setting_warnings(&settings),
            "--cluster-min-records 1 under --cluster-strategy legacy",
        );
    }

    #[test]
    fn single_record_clusters_are_expected_under_edit_mass() {
        let settings = ClusteringSettings {
            strategy: ClusterStrategy::EditMass,
            min_records: Some(1),
            ..Default::default()
        };
        assert!(clustering_setting_warnings(&settings).is_empty());
    }
}