mafft_core/

engine.rs

/// High-level MAFFT alignment engine.

use rayon::prelude::*;

use mafft_io::read_fasta;
use mafft_scoring::{build_context, build_context_with_kimura};
use mafft_tree::{DistanceMatrix, musclesupg, ktuple_distance, scoring_matrix_distance};
use mafft_tree::parttree_split::{build_parttree_topology};
use mafft_tree::parttree_pivot::PtSeqKind;
use mafft_align::{build_local_homology_table, GapModel};
use mafft_types::{ScoringModel, SequenceSet, LocalHomologyTable};

use crate::progressive::MultipleAlignment;
use crate::refinement::{iterative_refine, RefinementParams};
use crate::add::{add_sequences, add_sequences_keeplength};

/// Alignment mode (strategy).
#[derive(Debug, Clone)]
pub enum AlignmentMode {
    /// FFT-NS-2: fast progressive (default).
    FftNs2,
    /// FFT-NS-i: progressive + limited iterative refinement.
    FftNsi { iterations: usize },
    /// G-INS-i: global alignment + iterative refinement.
    GInsi { iterations: usize },
    /// L-INS-i: local alignment + iterative refinement.
    LInsi { iterations: usize },
    /// E-INS-i: generalized affine + iterative refinement.
    EInsi { iterations: usize },
    /// Q-INS-i: RNA alignment with McCaskill base-pair probabilities.
    QInsi { iterations: usize },
    /// X-INS-i: RNA alignment with CONTRAfold structure predictions.
    XInsi { iterations: usize },
}

impl Default for AlignmentMode {
    fn default() -> Self {
        Self::FftNs2
    }
}

/// Matrix shift applied when rebuilding the refinement-tree distances the
/// way C's `dndpre` does, for the modes that have no `pairlocalalign` step
/// (FFT-NS-i and friends).
///
/// `scripts/mafft` does NOT pass `-h` to the `dndpre` invocation that writes
/// `hat2` for `dvtditr`, so C's `constants()` falls back to its per-alphabet
/// DEFAULT `poffset` — and the two alphabets do not share one:
///
/// | alphabet   | default `poffset`          | `offset = (int)(600/1000 * poffset + 0.5)` | shift |
/// |------------|----------------------------|--------------------------------------------|-------|
/// | nucleotide | `DEFAULTOFS_N = -369` (`DNA.h:3`)    | -220                             | 220   |
/// | protein    | `DEFAULTOFS_B = -123` (`blosum.c:3`) | -73                              | 73    |
///
/// The DP matrix is shifted by `-offset`. Using the protein 73 for DNA made
/// the refinement distances differ from C's `hat2` outright — on
/// `mtb_cds_120x1400` the leaf pair (0,58) came out 0.307 against C's 0.253 —
/// which reordered UPGMA merges (steps 10 and 11 swapped), and a swapped
/// merge changes the group on 131 of 237 refinement branches. DNA FFT-NS-i
/// byte-parity with C over BAliBASE `bali2dna` went 67/141 → 132/141 when
/// this was corrected.
pub fn dndpre_offset_shift(is_nucleotide: bool) -> i32 {
    if is_nucleotide { 220 } else { 73 }
}

/// Scale factors turning the pair-phase `ppenalty`-style integers
/// (`lgop * 1000`, …) into DP units: `(gap_scale, offset_scale)`.
///
/// Mirrors C `constants.c`. Nucleotide (`:316-322`):
/// `penalty = (int)( 3 * 600.0/1000.0 * ppenalty + 0.5 )` and likewise
/// `penalty_ex` / `penalty_OP` / `penalty_dist`, but
/// `offset = (int)( 1 * 600.0/1000.0 * poffset + 0.5 )`. Protein
/// (`:672-677`): `600.0/1000.0` for all of them. The `3 *` on nucleotide
/// gap penalties is load-bearing — without it DNA pairwise gaps cost a
/// third of C's and L-INS-i / G-INS-i / E-INS-i diverge — so it lives in
/// one named place with a test pinning it.
pub fn pair_penalty_scales(is_nucleotide: bool) -> (f64, f64) {
    if is_nucleotide {
        (3.0 * 600.0 / 1000.0, 1.0 * 600.0 / 1000.0)
    } else {
        (600.0 / 1000.0, 600.0 / 1000.0)
    }
}

/// The main MAFFT alignment engine.
#[derive(Debug, Clone)]
pub struct MafftEngine {
    pub mode: AlignmentMode,
    pub scoring_model: ScoringModel,
    /// Number of guide tree rebuilds. C's FFT-NS-2 default is 2.
    pub retree: usize,
    /// Gap opening penalty override (positive float, e.g. 1.53 → internal -1530).
    /// None = use default for the scoring model.
    pub gap_open: Option<f64>,
    /// Offset/extension penalty override (positive float, e.g. 0.123 → internal -123).
    /// None = use default.
    pub gap_offset: Option<f64>,
    /// Gap extension penalty override (`--exp`). User passes a positive
    /// float (e.g. `--exp 0.1`); C negates internally (`gexp = -1.0 * arg`)
    /// and then `constants()` scales by `(int)(scale * gexp + 0.5)`
    /// (`scale = 600/1000` for protein, `3*600/1000` for DNA). When
    /// `Some`, overrides `scoring.gap.extend` in the same pattern as
    /// `gap_offset` overrides `scoring.gap.offset`. Default `None` keeps
    /// the model default (0 for protein/DNA).
    pub gap_extend: Option<f64>,
    /// Per-class pairwise gap params for L-INS-i (`--lop` / `--lep` /
    /// `--lexp`). Override the hardcoded `lgop=-2.00 / laof=0.100 /
    /// lexp=-0.100` C defaults applied in the local pairwise alignment
    /// stage. Each `None` = use the C default.
    pub pair_lop: Option<f64>,
    pub pair_lep: Option<f64>,
    pub pair_lexp: Option<f64>,
    /// Per-class pairwise gap params for E-INS-i generalized affine
    /// (`--gop` / `--gep` / `--gexp`). Override the C defaults
    /// `pggop=-1.53 / pgaof=0.10 / pgexp=-0.00`. Each `None` = use the
    /// C default.
    pub pair_gop: Option<f64>,
    pub pair_gep: Option<f64>,
    pub pair_gexp: Option<f64>,
    /// `--shiftpenalty` factor for `--allowshift`. C's `spfactor`
    /// multiplies the gap-open penalty to derive the per-cell shift
    /// cost: `penalty_shift = (int)(spfactor * penalty)`. Default 2.0
    /// (the C `--allowshift` baseline). `None` keeps the default.
    pub shift_penalty_factor: Option<f64>,
    /// `--minimumweight` floor applied to per-sequence weights in the
    /// intergroup-score accumulation. Threaded into
    /// `RefinementParams.minimum_weight`. `None` keeps the C default
    /// (`0.00001`).
    pub minimum_weight: Option<f64>,
    /// `--nwildcard`: fill the DNA scoring matrix's `'n'` row
    /// with `round(0.25 * self_score)` per residue, matching C's
    /// `constants.c::nscore`. When false (default and `--nzero`),
    /// the N row stays at the build-time defaults (effectively
    /// zero for the unscored entries). DNA-only — protein inputs
    /// silently bypass.
    pub nwildcard: bool,
    /// `--skipiterate F` threshold — mirrors C's `dvtditr -E
    /// $fixthreshold` → `autosubalignment = F`. When F exceeds the
    /// max distance-from-tip in the guide tree, refinement is
    /// skipped entirely (matches C's `generatesubalignmentstable`
    /// returning 1, which prints the "WARNING: Iterative refinement
    /// was not done" diagnostic and exits). For small F values, the
    /// port (R-3, closed 2026-06-03) generates sub-alignment
    /// clusters via `mafft_tree::generate_subalignments_table` and
    /// sets per-(step, side) skip flags in
    /// `RefinementParams::skip_branches`, mirroring C's
    /// `dvtditr.c:997-1006` `includemember && !samemember` gate.
    pub skipiterate: Option<f64>,
    /// `--bestfirst` refinement strategy. Default false (BAATARI2,
    /// matches C MAFFT's default). When true, refinement evaluates
    /// every branch against the same baseline alignment per
    /// iteration, picks the one with the largest gain, applies,
    /// repeats — mirroring C's `parallelizationstrategy = BESTFIRST`.
    pub bestfirst: bool,
    /// `--thread N`. C selects a different refinement implementation on
    /// `nthread > 0` (`tditeration.c:1433`): `athread` walks the tree in
    /// a fixed order and converges/stops by different rules than the
    /// single-threaded loop — see `RefinementParams::per_cycle_convergence`.
    /// `0` (the default, and what C's script passes for both no `--thread`
    /// and `--thread 0`) selects the single-threaded rules.
    pub nthread: usize,
    /// `--oneiteration` "one-vs-others" refinement (C's
    /// `disttbfast -r` → `dooneiteration` in
    /// `mafft-upstream/core/disttbfast.c:2217`). Runs once after the
    /// progressive merge and before regular refinement. ONLY
    /// triggered in the disttbfast-path modes (FFT-NS-2, FFT-NS-i);
    /// L/G/E-INS-i bypass it because `scripts/mafft:2673` only
    /// passes `-r` to `disttbfast`, never to `tbfast` or `dvtditr`.
    pub oneiteration: bool,
    /// `--pileup`: build a comb-tree guide via
    /// [`mafft_tree::Topology::pileup_chain`] instead of UPGMA,
    /// then run progressive merge once with no refinement (C
    /// strategy name "Pileup-NS-1", `scripts/mafft:2169`). Skips
    /// distance computation entirely. Forces `retree = 1` and
    /// disables `--maxiterate` refinement.
    pub pileup: bool,
    /// Tree-linkage method for UPGMA cluster joining. Mirrors C's
    /// `tbfast -X $sueff` (`scripts/mafft:264,419,422,427`). Default
    /// = `Mix { sueff: 0.1 }` (C default). `--averagelinkage` →
    /// `Mix { sueff: 1.0 }` ≡ `Average`; `--minimumlinkage` → `Mix
    /// { sueff: 0.0 }` ≡ `Minimum`; `--mixedlinkage F` → `Mix { sueff:
    /// F }`. C's `--youngestlinkage` is a separate algorithm
    /// (memory-saving k-mer tree builder with on-demand cluster
    /// distance recompute); rust wires it via `self.memsavetree`,
    /// which uses the same algorithm family; the dedicated
    /// `--youngestlinkage` port (`mafft_tree::youngestlinkage_tree`)
    /// is byte-identical to C (see `TODO.md` R-8).
    pub cluster_method: mafft_tree::ClusterMethod,
    /// Disable FFT: force pure DP for all alignment steps.
    pub nofft: bool,
    /// Enable long-range gap shift penalty (--allowshift). In MAFFT 7.526 the
    /// warp DP itself is dead code (`defs.c:54 trywarp = 0` and never set);
    /// the actual `--allowshift` effect is to set `unalign_level = 0.8` which
    /// triggers per-step `makedynamicmtx` (`disttbfast.c:2304`). Kept as a
    /// boolean for CLI symmetry; only `unalign_level > 0` has runtime effect.
    pub allowshift: bool,
    /// Per-step substitution-score offset = (distfromtip - unalign_level) * 600
    /// (clamped at 0). Mirrors C `specificityconsideration` + `dist2offset`
    /// + `makedynamicmtx`. 0 = disabled, 0.8 = `--allowshift` default.
    pub unalign_level: f64,
    /// Kimura R parameter for DNA distance model (--kimura).
    pub kimura_r: Option<i32>,
    /// Use PartTree for guide tree construction (--parttree).
    pub parttree: bool,
    /// Use DP-based PartTree (--dpparttree).
    pub dpparttree: bool,
    /// Group size for PartTree partitioning (--groupsize).
    pub groupsize: Option<usize>,
    /// Reorder output sequences in guide-tree DFS order (--reorder). Default
    /// is input order (--inputorder), matching C MAFFT 7.526.
    pub reorder_output: bool,
    /// Use a user-supplied guide tree (`--treein FILE`). Format matches C
    /// MAFFT's `_guidetree`: nseq-1 lines of `im jm len0 len1` (1-indexed,
    /// im < jm), as produced by `newick2mafft.rb`. When `Some`, distance
    /// computation and tree building are skipped — the loaded tree is
    /// used for every progressive pass (mirrors C `tbfast.c:2072-2078`).
    pub treein_path: Option<std::path::PathBuf>,
    /// Use the memory-saving guide-tree algorithm (`--memsavetree`).
    /// Mirrors C MAFFT `compacttree_memsaveselectable` with `howcompact=2`
    /// (`mltaln9.c:5491`) — k-mer-based distances computed on the fly with
    /// no full distance matrix. Enabled by `--auto` for the 100k+ bracket.
    pub memsavetree: bool,
    /// `--youngestlinkage` — same family as memsavetree but with per-step
    /// recomputation of cluster distances after each join. C MAFFT
    /// `mltaln9.c::compacttree_memsaveselectable(howcompact=2, memsave=1)`.
    pub youngestlinkage: bool,
    /// `--leavegappyregion` / `--legacygappenalty` — disable the
    /// gap-aware DP reweighting (`legacygapcost = 1`,
    /// `Salignmm.c:1604-1610`). Restores pre-7.110 behaviour where
    /// gappy columns are scored as if fully nongap.
    pub legacy_gap_cost: bool,
    /// Seed local-homology table (`--seed FILE` constraints). Mirrors
    /// C MAFFT's `hat3.seed` produced by `multi2hat3s` — pairwise
    /// `korh = 'k'` regions between seed sequences with `opt`
    /// pre-multiplied by `tsuyosa = user_nseq² * 100`. The table is
    /// sized to the full (seeds + user input) `nseq`. When `Some`,
    /// the engine folds these entries into its pairwise homology
    /// table (or uses them directly for non-INS-i modes) and forces
    /// `iterate ≥ 2` so the refinement step picks them up
    /// (`scripts/mafft:1911-1923`).
    pub seed_homology: Option<LocalHomologyTable>,
    /// `--memsave` Hirschberg DP routing. When true and the non-FFT
    /// progressive merge would call `profile_align`, route through
    /// `mafft_align::msalignmm` instead (linear-space DP — mirrors C
    /// MAFFT's `MSalignmm` in `tbfast.c:1159-1161` under `alg='M'`).
    /// For inputs that fit in memory the alignment is the same as
    /// `profile_align`; only memory usage differs.
    pub memsave_dp: bool,
    /// `--c-compat` opt-in: replicate C MAFFT's `Salignmm.c::A__align`
    /// static-TLS memoization (`reuseprofiles` / `cpmx_calc_add`) so
    /// tied-DP-cell choices match C bit-for-bit at the cost of
    /// carrying per-thread cross-call state. Default false (stateless,
    /// pure progressive engine). Enable to reproduce C's output on
    /// inputs where the §B.2 / BALIBASE-corpus residual divergences
    /// matter for downstream byte-equality requirements. See the
    /// BB20027 deep-dive in `balibase_parity_run.md` for the diagnosis
    /// (historical: every BALIBASE fixture is now byte-identical without it).
    pub c_compat: bool,
}

impl Default for MafftEngine {
    fn default() -> Self {
        Self {
            mode: AlignmentMode::FftNs2,
            scoring_model: ScoringModel::Blosum(62),
            retree: 2,
            gap_open: None,
            gap_offset: None,
            gap_extend: None,
            pair_lop: None,
            pair_lep: None,
            pair_lexp: None,
            pair_gop: None,
            pair_gep: None,
            pair_gexp: None,
            shift_penalty_factor: None,
            minimum_weight: None,
            skipiterate: None,
            bestfirst: false,
            nthread: 0,
            oneiteration: false,
            nwildcard: false,
            pileup: false,
            cluster_method: mafft_tree::ClusterMethod::default(),
            nofft: false,
            allowshift: false,
            unalign_level: 0.0,
            kimura_r: None,
            parttree: false,
            dpparttree: false,
            groupsize: None,
            reorder_output: false,
            treein_path: None,
            memsavetree: false,
            youngestlinkage: false,
            legacy_gap_cost: false,
            seed_homology: None,
            memsave_dp: false,
            c_compat: false,
        }
    }
}

impl MafftEngine {
    pub fn new(mode: AlignmentMode) -> Self {
        Self { mode, scoring_model: ScoringModel::Blosum(62), retree: 2,
            gap_open: None, gap_offset: None, gap_extend: None,
            pair_lop: None, pair_lep: None, pair_lexp: None,
            pair_gop: None, pair_gep: None, pair_gexp: None,
            shift_penalty_factor: None, minimum_weight: None,
            skipiterate: None, bestfirst: false, nthread: 0, oneiteration: false, nwildcard: false,
            pileup: false,
            cluster_method: mafft_tree::ClusterMethod::default(),
            nofft: false, allowshift: false, unalign_level: 0.0,
            kimura_r: None, parttree: false, dpparttree: false,
            groupsize: None, reorder_output: false, treein_path: None,
            memsavetree: false, youngestlinkage: false, legacy_gap_cost: false,
            seed_homology: None, memsave_dp: false, c_compat: false }
    }

    /// Enable `--c-compat`: replicate C MAFFT's static-TLS cpmx
    /// memoization so tied-DP-cell choices match C bit-for-bit.
    pub fn with_c_compat(mut self, c_compat: bool) -> Self {
        self.c_compat = c_compat;
        self
    }

    /// Set the number of guide tree rebuilds.
    pub fn with_retree(mut self, retree: usize) -> Self {
        self.retree = retree;
        self
    }

    /// Set gap opening penalty (positive float, e.g. 1.53).
    pub fn with_gap_open(mut self, op: f64) -> Self {
        self.gap_open = Some(op);
        self
    }

    /// Set offset/extension penalty (positive float, e.g. 0.123).
    pub fn with_gap_offset(mut self, ep: f64) -> Self {
        self.gap_offset = Some(ep);
        self
    }

    /// Use PartTree for guide tree (--parttree).
    pub fn with_parttree(mut self, parttree: bool) -> Self {
        self.parttree = parttree;
        self
    }

    /// Use DP-based PartTree (--dpparttree).
    pub fn with_dpparttree(mut self, dpparttree: bool) -> Self {
        self.dpparttree = dpparttree;
        self
    }

    /// Set group size for PartTree (--groupsize).
    pub fn with_groupsize(mut self, groupsize: usize) -> Self {
        self.groupsize = Some(groupsize);
        self
    }

    /// Set Kimura R parameter for DNA distance model (default 2).
    pub fn with_kimura(mut self, kimura_r: i32) -> Self {
        self.kimura_r = Some(kimura_r);
        self
    }

    /// Enable long-range gap shift penalty (CLI symmetry only — see field
    /// docstring; only `unalign_level > 0` has runtime effect).
    pub fn with_allowshift(mut self, allowshift: bool) -> Self {
        self.allowshift = allowshift;
        self
    }

    /// Set per-step dynamic-matrix offset (`specificityconsideration`).
    /// Disabled at 0.0; `--allowshift` defaults to 0.8.
    pub fn with_unalign_level(mut self, level: f64) -> Self {
        self.unalign_level = level;
        self
    }

    /// Disable FFT: force pure DP for all alignment steps.
    pub fn with_nofft(mut self, nofft: bool) -> Self {
        self.nofft = nofft;
        self
    }

    /// Emit output sequences in guide-tree DFS order (`--reorder`). When
    /// `false` (default), output stays in input order (`--inputorder`).
    pub fn with_reorder(mut self, reorder: bool) -> Self {
        self.reorder_output = reorder;
        self
    }

    /// Set scoring model (e.g. BLOSUM with specific number).
    pub fn with_scoring_model(mut self, model: ScoringModel) -> Self {
        self.scoring_model = model;
        self
    }

    /// Align a set of sequences.
    pub fn align(&self, input: &SequenceSet) -> MultipleAlignment {
        let seq_type = input.seq_type;
        let scoring_model = if seq_type.is_nucleotide() {
            ScoringModel::Dna
        } else {
            self.scoring_model
        };

        let mut scoring = if let Some(kr) = self.kimura_r {
            build_context_with_kimura(scoring_model, seq_type, kr)
        } else {
            build_context(scoring_model, seq_type)
        };

        // `--nwildcard` (and the implicit case from `unalignlevel != 0`
        // — see `scripts/mafft:1437` setting `nmodel=" -: "`). Fills
        // the DNA scoring matrix's `'n'` row with 25%-self-score
        // values. DNA-only — protein inputs no-op. Apply BEFORE any
        // gap penalty overrides so they still take effect on the
        // residue submatrix.
        let want_nwildcard = self.nwildcard || self.unalign_level > 0.0;
        if want_nwildcard {
            mafft_scoring::apply_nwildcard(&mut scoring);
        }

        // Apply gap penalty overrides if set.
        // C convention: --op 1.53 means ppenalty = -1530 (multiply by -1000).
        // After scaling: penalty = (int)(600/1000 * ppenalty + 0.5).
        if let Some(op) = self.gap_open {
            let ppenalty = -(op * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 3.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            scoring.gap.open = (scale * ppenalty as f64 + 0.5) as i32;
        }
        if let Some(ep) = self.gap_offset {
            let poffset = -(ep * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 1.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            let new_offset = (scale * poffset as f64 + 0.5) as i32;
            // C's constants() bakes the offset into the scoring matrix during
            // construction: `n_distmp[i][j] -= offset`. Our build_context()
            // builds the matrix with offset=0 (matching C's default aof=0 from
            // the shell script), NOT with gap_params.offset. So the "old"
            // offset baked into the matrix is 0, regardless of what
            // scoring.gap.offset says.
            let matrix_offset = 0i32;
            let delta = new_offset - matrix_offset;
            if delta != 0 {
                let nscored = scoring.nscoredalphabets;
                for i in 0..nscored {
                    for j in 0..nscored {
                        scoring.substitution_matrix[i][j] -= delta;
                        scoring.consweight_matrix[i][j] = scoring.substitution_matrix[i][j] as f64;
                        scoring.fft_matrix[i][j] = scoring.substitution_matrix[i][j] + new_offset;
                    }
                }
            }
            scoring.gap.offset = new_offset;
        }
        // `--exp` (gap extension penalty). C: `gexp = -1.0 * arg` then
        // `penalty_ex = (int)(scale * gexp + 0.5)` where `scale =
        // 600/1000` (protein) or `3*600/1000` (DNA). Same shape as the
        // `gap_open` override above.
        if let Some(exp) = self.gap_extend {
            let pgexp = -(exp * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 3.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            scoring.gap.extend = (scale * pgexp as f64 + 0.5) as i32;
        }

        let nseq = input.nseq();
        let quiet_mode = false;
        // C `disttbfast.c:4450` (`tbfast.c:3320`, `pairlocalalign.c:3372`)
        // calls `gappick0(bseq[i], seq[i])` for every sequence before
        // anything is aligned — the engine works on RESIDUE-ONLY
        // sequences regardless of whether the input FASTA had gaps.
        // Without this, feeding a previously-aligned FASTA (gaps in
        // input) produces a different alignment than C because the rust
        // progressive sees the gapped form (closes R-6: 16+1 adversarial
        // fixture diverged by 369 lines on combined_17.fa, byte-identical
        // on the residue-only c17_ungapped.fa).
        //
        // Only `-` is a gap here. `gappick0` (`mltaln9.c:10537`) leaves
        // `.` alone: for protein it is a legal residue (`amino_n['.']` =
        // 23, `blosum.c:12`), for nucleotide it is illegal and C exits
        // before reaching this point (`seqcheck`, `mltaln9.c:60`).
        let sequences: Vec<Vec<u8>> = input.sequences.iter().map(|s| {
            s.data.iter().copied().filter(|&c| c != b'-').collect()
        }).collect();
        let names: Vec<String> = input.sequences.iter().map(|s| s.name.clone()).collect();

        let use_fft = !self.nofft && matches!(
            self.mode,
            AlignmentMode::FftNs2 | AlignmentMode::FftNsi { .. }
        );

        // Step 1: Initial guide tree
        // For PartTree mode, route through the C-equivalent splittbfast
        // pipeline (`crates/mafft-tree/src/parttree_split.rs`). Both
        // `--parttree` and `--dpparttree` route through it — the C
        // distinction (`partdist="ktuples"` vs `partdist="localalign"`,
        // `scripts/mafft:392/395`) is a distance-metric variation
        // inside C's `splittbfast`; our Rust pipeline currently uses
        // k-tuple distance for both. The 36-seq sample is below
        // PartTree's recursion threshold so both modes byte-identical
        // C MAFFT 7.526. A true DP-based distance for `--dpparttree`
        // larger-input runs is not yet ported (would slot into
        // `parttree_dist.rs`).
        let use_parttree = self.parttree || self.dpparttree;
        let parttree_topo = if use_parttree {
            let kind = if scoring.seq_type.is_nucleotide() {
                PtSeqKind::Dna
            } else {
                PtSeqKind::Protein
            };
            let picksize = 50;
            Some(build_parttree_topology(&sequences, kind, picksize))
        } else {
            None
        };

        // For L-INS-i / E-INS-i, MAFFT's `pairlocalalign` (then `tbfast`) replaces
        // the 6-mer initial distance with a distance derived from all-vs-all
        // pairwise local alignments. We piggyback on `build_local_homology_table`,
        // which already runs the same pairwise alignments to populate the
        // homology constraint table — we keep both outputs (distance + table).
        let pair_kind = match self.mode {
            AlignmentMode::LInsi { .. } => Some(mafft_align::PairAligner::Local),
            AlignmentMode::GInsi { .. } => Some(mafft_align::PairAligner::Global),
            AlignmentMode::EInsi { .. } => {
                Some(mafft_align::PairAligner::GeneralizedAffine)
            }
            _ => None,
        };
        let mut pairwise_for_constraints = if let Some(aligner) = pair_kind {
            let seq_refs: Vec<&[u8]> = sequences.iter()
                .map(|s| s.as_slice()).collect();
            // C's `pairlocalalign` uses pairwise-specific gap penalties,
            // NOT the progressive ones (`scripts/mafft:91-92,201-203`).
            // For L-INS-i (`-L`): lgop=-2.00, lexp=-0.100, laof=0.100.
            // For G-INS-i (`-A`): pgop=$pggop, pgexp=$pggexp, pgaof=$pgaof.
            //   Defaults match L-INS-i values for protein
            //   (`scripts/mafft:91-92`). Same numbers below.
            // C's argument parser (`pairlocalalign.c:1671-1675`) does
            //   ppenalty = (int)( atof(arg) * 1000 - 0.5 )
            // which is C-style truncation toward zero (`as i32` in Rust).
            // For `-f -2.00` this gives -2000 (not -2001). Then
            // `constants.c:1014-1016` does
            //   penalty = (int)( 600/1000 * ppenalty + 0.5 )
            // which is round-half-up for positive and round-half-up-toward-
            // zero for negative — also `as i32` truncation in Rust because
            // for negative numbers like -1199.5, `(int)` gives -1199 not
            // -1200.
            //
            // For `-2.00 / -0.100 / 0.100`: C gets penalty = -1199,
            // penalty_ex = -59, offset = 59. Round-naively in Rust we'd
            // get -1200 / -60 / 60 — off by 1, which propagates through
            // `iscore` and yields a 1-per-residue gap in `opt` (~0.01
            // off vs C across all pairs).
            let cc_int = |x: f64, mul: f64| -> i32 {
                ((x * mul) - 0.5) as i32
            };
            let cc_scale = |ppen: i32, scale: f64| -> i32 {
                ((scale * ppen as f64) + 0.5) as i32
            };
            // L-INS-i / G-INS-i defaults (script:91-92,201-203).
            // E-INS-i overrides (`scripts/mafft:1940-1948`): when
            // distance="localgenaf" (and `oldgenafparam != 1`), the script
            // resets `lexp="0.0"` and `laof="0.0"` so the regular gap-extend
            // and matrix-offset are zeroed out, leaving only the gen-affine
            // skip-gap (LGOP) as the long-range penalty.
            let is_einsi = matches!(self.mode, AlignmentMode::EInsi { .. });
            // `scripts/mafft:1469-1473`: when `unalignlevel > 0` zero
            // `lexp=laof=pgexp=pgaof=0` for the pair phase.
            let unalign_active = self.unalign_level > 0.0;
            // Pair-phase gap defaults. C scripts/mafft assigns these per
            // mode. Both L-INS-i and E-INS-i pairwise alignment use
            // `lgop`/`lexp`/`laof` (the "L-INS" gap params); `--lop` /
            // `--lep` / `--lexp` override. The `pggop`/`pgaof`/`pgexp`
            // family (set by `--gop`/`--gep`/`--gexp` in C) is reserved
            // for the X-INS-i / Q-INS-i RNA pipelines which we don't
            // exercise — surfacing them at the CLI but they are
            // currently inert for protein/DNA workflows. For E-INS-i,
            // `lexp` and `laof` are forced to 0 (the generalized-affine
            // skip-gap cost — `lgop_op = LGOP = -6.00` below — covers
            // long-range gaps instead).
            let lgop: f64 = self.pair_lop.unwrap_or(-2.00);
            let lexp: f64 = if is_einsi || unalign_active {
                self.pair_lexp.unwrap_or(0.0)
            } else {
                self.pair_lexp.unwrap_or(-0.100)
            };
            let laof: f64 = if is_einsi || unalign_active {
                self.pair_lep.unwrap_or(0.0)
            } else {
                self.pair_lep.unwrap_or(0.100)
            };
            // E-INS-i extras (`scripts/mafft:198-199`):
            //   LGOP=-6.00 → ppenalty_OP (skip-gap open).
            //   LEXP= 0.0 → ppenalty_EX (skip-gap extend, unused; C
            //               comments out the extension increments).
            let lgop_op: f64 = -6.00;
            // C scales the pair-phase penalties differently for nucleotide
            // and protein (`constants.c:316-322` vs `:672-677`):
            //   nucleotide: penalty/penalty_ex/penalty_OP = 3 * 600/1000 * pp
            //               offset                        = 1 * 600/1000 * po
            //   protein:    everything                    =     600/1000 * pp
            // The gap penalties carry the `3 *`; the offset does not (C
            // writes the `1 *` out explicitly beside the `3 *`s). Applying
            // the protein factor to DNA made pairwise gaps a third of C's
            // cost, so L-INS-i / G-INS-i / E-INS-i opened gaps C refused.
            // Same idiom as the progressive-phase penalties above.
            let (gap_scale, offset_scale) = pair_penalty_scales(seq_type.is_nucleotide());
            let p_open = cc_int(lgop, 1000.0);
            let p_ext = cc_int(lexp, 1000.0);
            let p_offset = cc_int(laof, 1000.0);
            let p_op = cc_int(lgop_op, 1000.0);
            let mut pair_gap = GapModel::new(
                cc_scale(p_open, gap_scale) as f64,
                cc_scale(p_ext, gap_scale) as f64,
            );
            // C `constants.c:277-278`: `if (penalty_shift_factor < 10) trywarp = 1`.
            // With `--allowshift`, `spfactor = 2.0` (< 10) → warp DP fires.
            // `penalty_shift = (int)(penalty_shift_factor * penalty)`
            // (`constants.c:318`). For pair phase: penalty = -1199, sp = 2.0,
            // so penalty_shift = -2398.
            if self.unalign_level > 0.0 {
                let spfactor = self.shift_penalty_factor.unwrap_or(2.0);
                let penalty_shift = (spfactor * pair_gap.open) as i32 as f64;
                pair_gap.shift = Some(penalty_shift);
            }
            let pair_op = cc_scale(p_op, gap_scale) as f64;
            let pair_offset_int: i32 = cc_scale(p_offset, offset_scale);
            let nscored = scoring.nscoredalphabets;
            // The DP layer takes f64 matrices (post §9c migration). Build
            // the shifted matrix from `consweight_matrix` (= f64 view of
            // `substitution_matrix`) and subtract the pair offset there.
            let pair_offset_f64 = pair_offset_int as f64;
            let mut shifted: Vec<Vec<f64>> = scoring.consweight_matrix.clone();
            for i in 0..nscored {
                for j in 0..nscored {
                    shifted[i][j] -= pair_offset_f64;
                }
            }
            // C's `L__align11` sets `localthr = -offset + scoreoffset * 600`
            // (Lalign11.c:248-249). With `scoreoffset = 0` and the
            // `offset = (int)(0.6 * poffset + 0.5)` computed above
            // (= `pair_offset_int`), C uses `localthr = -pair_offset_int`.
            // Our `local_align` computes `localthr = -score_offset * 600`,
            // so to reach `localthr = -pair_offset_int` we pass
            // `score_offset = pair_offset_int / 600`.
            let score_offset_for_local = pair_offset_int as f64 / 600.0;
            let (table, dist) = mafft_align::build_homology_table_with_unalign(
                &seq_refs,
                &shifted,
                &scoring.amino_map,
                &pair_gap,
                score_offset_for_local,
                aligner,
                pair_op,
                self.unalign_level,
            );
            // For L-INS-i, tbfast computes pairwise alignments in-memory
            // via `callpairlocalalign=1`. The `iscore` distance matrix
            // is passed straight to `fixed_musclesupg_double_realloc_…`
            // without a hat2 file round-trip (the script does not invoke
            // a separate pairlocalalign + hat2 read), so we keep full
            // double precision here. The 3-decimal hat2 rounding only
            // applies on paths that genuinely write/read `hat2` (e.g.
            // FFT-NS-i + dndpre).
            Some((table, DistanceMatrix::from_full(&dist)))
        } else {
            None
        };

        let mut dm = if use_parttree {
            // Skip full distance matrix — PartTree builds tree directly
            DistanceMatrix::new(nseq)
        } else if let Some((_, ref pre_dm)) = pairwise_for_constraints {
            pre_dm.clone()
        } else {
            compute_distance_matrix_from_seqs(&sequences)
        };

        // Load user-supplied guide tree early (`--treein`), so the
        // importance-recomputation below uses the same topology C does
        // (`tbfast.c:2967 counteff_simple_double_nostatic_memsave( njob, topol, len, dep, eff )`
        // where `topol`/`len` come from `loadtree` when `treein=1`).
        let user_topo: Option<mafft_tree::Topology> = self.treein_path.as_ref().map(|path| {
            mafft_tree::parse_mafft_tree(path, nseq).unwrap_or_else(|e| {
                eprintln!("--treein: {e}");
                std::process::exit(1);
            })
        });

        // `--seed`: fold seed-derived `hat3.seed` entries into the
        // pairwise homology table BEFORE `recompute_importance`, so the
        // position-vote pass weighs seed regions together with pairwise
        // ones. Mirrors C MAFFT's `cat hat3.seed hat3 > hat3`
        // (`scripts/mafft:2523-2540`) — tbfast then reads the merged
        // file before calling `calcimportance_half`.
        //
        // When the mode doesn't run pairwise homology (FFT-NS-i with
        // `--seed`), promote the seed table to be the constraint table
        // outright; the pairwise distance matrix is left alone (FFT-NS-i
        // uses ktuple distances by default — same as without `--seed`).
        if let Some(ref seed_lh) = self.seed_homology {
            if let Some((ref mut table, _)) = pairwise_for_constraints {
                mafft_align::merge_homology_tables(table, seed_lh);
            }
        }

        // C's `tbfast` calls `calcimportance_half` (mltaln9.c:11756) AFTER
        // the initial tree to replace each region's provisional importance
        // with `mean(position-vote support over region) * region.opt`,
        // then symmetrize across (i,j)/(j,i). `region.opt` is stored in
        // C's post-`tbfast.c:2202` scale (`isumscore / sumoverlap`), so
        // the impmtx contributions match C's numerically.
        if pairwise_for_constraints.is_some() && !use_parttree {
            // C computes the constraint `importance` TWICE, with DIFFERENT
            // trees, and the two phases must be mirrored separately:
            //
            //  1. tbfast (progressive) calls `calcimportance` using weights
            //     from the FULL-PRECISION in-memory `iscore` UPGMA tree
            //     (`tbfast.c:2926` builds it from the un-rounded distance
            //     matrix). This importance drives the progressive merge.
            //  2. dvtditr (refinement) RE-reads the original hat3 and calls
            //     `calcimportance` again, this time with weights from the
            //     3-decimal `hat2` tree (`readhat2_pointer`, dvtditr.c:753).
            //
            // This block is phase (1): use the full-precision `dm` tree.
            // Phase (2) is mirrored just before `iterative_refine` below,
            // where `local_hom`'s importance is recomputed from the rounded
            // refinement tree. Using the rounded tree here instead would
            // regress E-INS BB40004 (progressive merge diverges, 588 lines);
            // using the full tree for refinement would regress BB50001
            // (244 lines). Verified both directions empirically.
            //
            // With `--treein`, C uses the loaded user tree for both phases.
            let initial_topo = user_topo.clone()
                .unwrap_or_else(|| musclesupg(&dm, self.cluster_method));
            let weights = mafft_tree::sequence_weights(&initial_topo);
            let seq_refs: Vec<&[u8]> = sequences.iter()
                .map(|s| s.as_slice()).collect();
            if let Some((ref mut table, _)) = pairwise_for_constraints {
                mafft_align::recompute_importance(table, &seq_refs, &weights);
            }
        }

        // Step 2: Build guide tree and progressive align, repeating `retree` times.
        // Each iteration after the first computes distances from the ALIGNMENT
        // (not the raw sequences), producing a better tree.
        //
        // C's `scripts/mafft:142-156` sets `defaultcycle=1` for L/G/E/Q/X-INS-i
        // (vs `defaultcycle=2` for FFT-NS-2 and FFT-NS-i). The script then
        // applies two post-`--retree` rewrites that override the user value:
        //
        //   1. `scripts/mafft:1840-1842` clamps `cycle = min(cycle, 3)` for
        //      ALL distance modes — `--retree 5` becomes 3 in C.
        //   2. `scripts/mafft:1934-1936` forces `cycle = 1` in the
        //      `distance ∈ {local, global, localgenaf, globalgenaf, scarna}`
        //      branch (= L/G/E/Q/X-INS-i) regardless of `--retree`. So
        //      `--retree 3 --localpair` runs with cycle=1 in C, not 3.
        //
        // Mirror both: cap at 3 for FFT-NS-*, force 1 for INS-i regardless
        // of the user value. Without this, `--retree N --localpair` (N > 1)
        // silently runs extra passes vs C and diverges by hundreds of lines
        // (898 lines on the 36-seq sample with `--retree 3 --localpair`).
        let retree = if matches!(
            self.mode,
            AlignmentMode::LInsi { .. }
                | AlignmentMode::GInsi { .. }
                | AlignmentMode::EInsi { .. }
                | AlignmentMode::QInsi { .. }
                | AlignmentMode::XInsi { .. }
        ) {
            1
        } else if self.pileup {
            // `--pileup` is "Pileup-NS-1" — single progressive pass,
            // no second tree-rebuild (`scripts/mafft:2169` strategy
            // name; C disttbfast forces `cycle = 1` for the pileup
            // guidetree variant at line 1991).
            1
        } else {
            self.retree.clamp(1, 3)
        };
        let mut msa = MultipleAlignment {
            sequences: sequences.clone(),
            names: names.clone(),
            score: 0.0,
            step_trace: Vec::new(), guide_tree: None, first_pass_sequences: None, distance_matrix: None,
        };
        let mut accumulated_trace = Vec::new();
        let penalty_dist = scoring.gap.open;
        // Final progressive guide tree — used for `--reorder` output ordering
        // (mirrors C `tbfast.c:2928` writing the order file from the
        // post-UPGMA topology, BEFORE any iterative refinement).
        let mut final_progressive_topo: Option<mafft_tree::Topology> = None;
        // Intermediate alignment after pass 0 of the retree loop. C's
        // `--parttree` script runs `splittbfast` TWICE: CALL 1 produces
        // `pre_1` (this is what we want to capture here) and CALL 2 reads
        // `pre_1` as `orialn` for its `naivepairscore11`-based distance
        // computation. Without this we'd feed CALL 2 the FINAL alignment
        // (`pre_2`) and the scores would diverge.
        let mut first_pass_msa: Option<Vec<Vec<u8>>> = None;

        // `user_topo` was loaded above (before recompute_importance) so the
        // LH table weights and the progressive merge tree are derived from
        // the SAME topology — matching C's `tbfast.c:2072` (loadtree) +
        // `tbfast.c:2967` (counteff_simple from loaded topol) sequencing.

        // `--memsavetree`: build the guide tree with the C MAFFT
        // `compacttree_memsaveselectable` algorithm. C uses k-mer-based
        // `distcompact` in disttbfast (pass 0, raw sequences) and switches
        // to MSA-based `distcompact_msa` in tbfast (pass 1+, aligned
        // sequences). Mirrors `disttbfast.c:4018` then
        // `tbfast.c:2538`.
        //
        // Pass 0 uses k-mer; pass 1+ uses MSA. The MSA tree is rebuilt
        // INSIDE the retree loop from `msa.sequences` after each pass.
        let memsavetree_kmer_topo: Option<mafft_tree::Topology> = if self.memsavetree || self.youngestlinkage {
            let seq_refs: Vec<&[u8]> = sequences.iter()
                .map(|s| s.as_slice()).collect();
            let is_dna = scoring.seq_type.is_nucleotide();
            if self.youngestlinkage {
                Some(mafft_tree::memsavetree::youngestlinkage_tree(&seq_refs, is_dna))
            } else {
                Some(mafft_tree::memsavetree::memsavetree(&seq_refs, is_dna))
            }
        } else {
            None
        };

        for pass in 0..retree {
            let topo = if let Some(ref t) = user_topo {
                t.clone()
            } else if self.pileup {
                // `--pileup`: skip the distance matrix entirely and
                // build a comb-tree directly from the input order.
                // C `mltaln9.c::createchain` with `shuffle=0`.
                mafft_tree::Topology::pileup_chain(nseq)
            } else if self.memsavetree || self.youngestlinkage {
                if pass == 0 {
                    memsavetree_kmer_topo.clone().expect("memsavetree topology must be cached")
                } else {
                    // MSA-based rebuild from the prior pass's alignment.
                    // For --youngestlinkage, use the compacttree=4 MSA
                    // variant (`youngestlinkage_tree_msa`); for
                    // --memsavetree, use the compacttree=3 MSA variant
                    // (`memsavetree_msa`).
                    let aligned_refs: Vec<&[u8]> = msa.sequences.iter()
                        .map(|s| s.as_slice()).collect();
                    if self.youngestlinkage {
                        mafft_tree::memsavetree::youngestlinkage_tree_msa(
                            &aligned_refs,
                            &scoring.consweight_matrix,
                            &scoring.amino_map,
                            scoring.gap.open as f64,
                        )
                    } else {
                        mafft_tree::memsavetree::memsavetree_msa(
                            &aligned_refs,
                            &scoring.consweight_matrix,
                            &scoring.amino_map,
                            scoring.gap.open as f64,
                        )
                    }
                }
            } else if pass == 0 && use_parttree {
                parttree_topo.clone().unwrap()
            } else {
                musclesupg(&dm, self.cluster_method)
            };

            // C always progresses from raw input on each retree pass.
            let input_seqs = sequences.clone();

            // Shift penalty: penalty_shift = penalty_shift_factor * penalty
            // (`constants.c:318`). Default factor = 100 (trywarp = 0). With
            // `--allowshift`, `scripts/mafft:1428` sets `spfactor=2.00`, which
            // triggers `trywarp = 1` (constants.c:277-278: `if (factor < 10)`)
            // and gives `penalty_shift = 2.0 * penalty`. The previous "0.8"
            // here was confused with `unalignlevel = 0.8` — different knob.
            let shift = if self.allowshift {
                let spfactor = self.shift_penalty_factor.unwrap_or(2.0);
                Some(spfactor * scoring.gap.open as f64)
            } else {
                None
            };

            // For L-INS-i / E-INS-i, thread the local-homology table through
            // the progressive merges so they pick up the same per-cell
            // importance bonuses the refinement DP already uses. C's tbfast
            // does this via Falign_localhom (FFT) or partA__align (per
            // segment). We currently only handle the non-FFT branch in
            // `progressive_align_with_constraints` — that's the path
            // L-INS-i takes since the engine sets `use_fft = false` for
            // any non-FftNs2/FftNsi mode.
            let progress_constraints = pairwise_for_constraints
                .as_ref().map(|(t, _)| t);
            // C's `tbfast` is invoked with different `outgap` settings per
            // mode (`scripts/mafft:2584,2593,2601`): G-INS-i omits the
            // `$termgapopt = -O` flag so `outgap = 1` (head/tail gap
            // penalized). L-INS-i and E-INS-i pass `-O` so `outgap = 0`.
            // The progressive A__align/profile_align_imp call propagates
            // this as `headgp = tailgp = outgap`.
            // C `outgap=1` (terminal gaps penalized) is the global default
            // (`splittbfast.c:560`, `disttbfast.c:185`) and is overridden to
            // 0 by the `-O` flag (`scripts/mafft:291 termgapopt=" -O "`).
            // The mafft script passes `-O` to disttbfast/tbfast for most
            // modes but withholds it for:
            //   - `--globalpair` (G-INS-i / G-INS-1) — `scripts/mafft:2584`
            //     vs L-INS-i/E-INS-i which include termgapopt.
            //   - `--parttree` / `--dpparttree` — `scripts/mafft:2655` does
            //     not include `$termgapopt` in the splittbfast call.
            let penalize_term_gaps = matches!(self.mode, AlignmentMode::GInsi { .. })
                || use_parttree;
            // C `splittbfast.c:6` `#define WEIGHT 0` makes `--parttree` use
            // `fastconjuction_noweight` (uniform per-cluster weights) for
            // its internal `pairalign`. We mirror that by passing a
            // uniform-1.0 weight vector when `use_parttree`. `--pileup`
            // also disables weighting (`disttbfast.c:3962` sets
            // `weight = 0; tbrweight = 0` → `eff[i] = 1.0` for all i at
            // `disttbfast.c:4131`). All other modes derive weights from
            // the guide tree's branch lengths.
            let weights_override: Option<Vec<f64>> = if use_parttree || self.pileup {
                Some(vec![1.0; sequences.len()])
            } else {
                None
            };
            // Disable the cached profile blend (`blend_profiles_exact`).
            // C's `createcpmxresult` (Salignmm.c:608) omits the eff*1.0 gap
            // contribution at gap-insertion positions ("tsukawanai" comment
            // at line 624). C also only uses the cache in pass 0 (`treebase`)
            // and forces fresh `cpmx_calc_new` in pass 1+ (`dooneiteration`
            // — disttbfast.c:2288-2289 sets `cpmxchild0/1 = NULL`). Rust's
            // blend mirrors C's createcpmxresult math-faithfully, but
            // matching when the cache fires across all BB tests is delicate
            // — disabling the blend entirely (always `cpmx_calc_new`) gives
            // byte-identical output across the BB20018 / BB40046 set without
            // perf impact at typical alignment sizes.
            msa = crate::progressive::progressive_align_full_c_compat_ex(
                &input_seqs, &names, &topo, &scoring, use_fft, shift,
                progress_constraints, penalize_term_gaps,
                weights_override.as_deref(), self.unalign_level,
                self.legacy_gap_cost, self.memsave_dp, self.c_compat,
                false,
            );
            accumulated_trace.extend(msa.step_trace.iter().copied());
            final_progressive_topo = Some(topo.clone());
            if pass == 0 && first_pass_msa.is_none() {
                first_pass_msa = Some(msa.sequences.clone());
            }

            // For the next retree pass, recompute distances from the now-aligned
            // sequences (matching C's disttbfast behavior in the second iteration
            // of `iguidetree`). The refinement-tree distance matrix is built
            // separately further down with a different (offset-shifted) matrix
            // mirroring dndpre's invocation.
            if pass + 1 < retree {
                dm = compute_distance_matrix_scoring(
                    &msa.sequences, &scoring.substitution_matrix,
                    &scoring.amino_map, penalty_dist,
                );
            }
        }
        msa.step_trace = accumulated_trace;

        // Step 3: Build local homology table (for constrained modes)
        let uses_constraints = matches!(
            self.mode,
            AlignmentMode::LInsi { .. }
                | AlignmentMode::EInsi { .. }
                | AlignmentMode::GInsi { .. }
        );
        let uses_rna_constraints = matches!(
            self.mode,
            AlignmentMode::QInsi { .. } | AlignmentMode::XInsi { .. }
        );
        let mut local_hom = if uses_rna_constraints {
            // RNA modes: compute base-pair probabilities using external tools,
            // then use them as constraints for iterative refinement.
            let bpp_result = match &self.mode {
                AlignmentMode::QInsi { .. } => {
                    crate::external::compute_bpp_mccaskill(
                        &sequences,
                    )
                }
                AlignmentMode::XInsi { .. } => {
                    crate::external::compute_bpp_contrafold(
                        &sequences,
                    )
                }
                _ => unreachable!(),
            };
            match bpp_result {
                Ok(bpp_tables) => {
                    if !quiet_mode {
                        eprintln!("RNA structure: computed BPP for {} sequences", bpp_tables.len());
                    }
                    // For now, use standard local homology as fallback.
                    // Full BPP→constraint integration would convert base-pair
                    // probabilities into pairwise constraints here.
                    let seq_refs: Vec<&[u8]> = sequences.iter().map(|s| s.as_slice()).collect();
                    let gap = GapModel::new(scoring.gap.open as f64, scoring.gap.extend as f64);
                    let (table, _dist) = build_local_homology_table(
                        &seq_refs,
                        &scoring.consweight_matrix,
                        &scoring.amino_map,
                        &gap,
                        0.0,
                    );
                    Some(table)
                }
                Err(e) => {
                    eprintln!("{e}");
                    std::process::exit(1);
                }
            }
        } else if uses_constraints {
            // Reuse the table built up-front for the initial distance matrix.
            // Keep the initial pairwise distance matrix alongside for the
            // refinement tree (mirrors C's `dvtditr` reading the `hat2` file
            // written by initial pairlocalalign instead of recomputing from
            // the progressive alignment).
            pairwise_for_constraints.as_ref().map(|(t, _)| t.clone())
        } else if let Some(ref seed_lh) = self.seed_homology {
            // `--seed` with a non-INS-i mode (e.g. FFT-NS-i + `--seed`):
            // no pairwise homology was built, but seed entries still
            // need to drive refinement. Run `recompute_importance` on
            // the seed-only table here (using a sequence-weight vector
            // derived from the final progressive guide tree, mirroring
            // `tbfast.c:2967` calling `calcimportance` after the post-
            // progressive tree is in hand).
            let mut table = seed_lh.clone();
            let seq_refs: Vec<&[u8]> = msa.sequences.iter()
                .map(|s| s.as_slice()).collect();
            let weights = final_progressive_topo.as_ref()
                .map(mafft_tree::sequence_weights)
                .unwrap_or_else(|| vec![1.0; nseq]);
            mafft_align::recompute_importance(&mut table, &seq_refs, &weights);
            Some(table)
        } else {
            None
        };
        // Stash the initial pairwise distance matrix for the refinement tree
        // (modes that ran pairlocalalign, where C's dvtditr reads `hat2`).
        let initial_pairwise_dm: Option<DistanceMatrix> = pairwise_for_constraints
            .as_ref()
            .map(|(_, dm)| dm.clone());

        // `--oneiteration`: C's `disttbfast -r` → `dooneiteration`
        // (`disttbfast.c:2217-2538`). Runs AFTER progressive merge,
        // BEFORE regular refinement. Gated on the disttbfast-path
        // modes ONLY (FFT-NS-2, FFT-NS-i) — `scripts/mafft:2673`
        // passes `-r` only to `disttbfast`, not to `tbfast`/`dvtditr`.
        // L/G/E-INS-i use pairlocalalign+tbfast, so their pipeline
        // never sees `-r` even when `--oneiteration` is given (C
        // confirmed: `--localpair --oneiteration` ≡ `--localpair`
        // alone, byte-identical).
        if self.oneiteration
            && matches!(self.mode,
                AlignmentMode::FftNs2 | AlignmentMode::FftNsi { .. })
        {
            // `final_progressive_topo` is the guide tree from the
            // last `progressive_align_full_c_compat_ex` pass
            // (set at line 802 above). Always Some at this point
            // for the disttbfast-path modes (FFT-NS-2/i).
            let oneiter_topo = final_progressive_topo
                .as_ref()
                .expect("FFT-NS-2/i path always sets final_progressive_topo")
                .clone();
            let refine_shift = if self.allowshift {
                let spfactor = self.shift_penalty_factor.unwrap_or(2.0);
                Some((spfactor * scoring.gap.open as f64) as i32 as f64)
            } else {
                None
            };
            let one_iter_params = RefinementParams {
                max_iterations: 0, // unused by one_vs_others_refine
                use_fft: true,
                legacy_gap_cost: self.legacy_gap_cost,
                shift: refine_shift,
                unalign_level: self.unalign_level,
                minimum_weight: self.minimum_weight.unwrap_or(0.00001),
                ..Default::default()
            };
            crate::refinement::one_vs_others_refine(
                &mut msa, &oneiter_topo, &scoring, &one_iter_params,
            );
        }

        // Step 4: Iterative refinement (if mode requires it)
        match &self.mode {
            AlignmentMode::FftNs2 => {}
            AlignmentMode::FftNsi { iterations }
            | AlignmentMode::GInsi { iterations }
            | AlignmentMode::LInsi { iterations }
            | AlignmentMode::EInsi { iterations }
            | AlignmentMode::QInsi { iterations }
            | AlignmentMode::XInsi { iterations } => {
                // C's mafft script does NOT pass `-h` to dndpre (the second
                // invocation that writes hat2 for dvtditr). dndpre therefore
                // uses the BLOSUM62 default `poffset = -123` → offset = -73,
                // shifting every cell of the scoring matrix by +73 relative
                // to the offset=0 matrix disttbfast and dvtditr use for DP.
                // The refinement tree dvtditr builds reads this hat2, so the
                // distance matrix we feed `musclesupg` here must use the same
                // shifted matrix and operate on the FINAL progressive
                // alignment (`msa.sequences`).
                // For modes that ran an initial pairlocalalign step (L-INS-i,
                // G-INS-i, E-INS-i, ...), C's `dvtditr` reads `hat2` written
                // by tbfast's pairlocalalign — initial pairwise distances at
                // 3-decimal precision. We mirror that exactly.
                //
                // For modes without pairlocalalign (FFT-NS-i, distance="ktuples"),
                // C's script invokes `dndpre` between tbfast and dvtditr to
                // recompute distances from the progressive alignment with
                // `dndpre`'s DEFAULT poffset shift. We mirror that path here.
                let dm = if let Some(ref initial_dm) = initial_pairwise_dm {
                    // Mimic hat2 file's `%.3f` rounding so musclesupg sees the
                    // same distances dvtditr sees.
                    let n = initial_dm.nseq;
                    let mut rounded = DistanceMatrix::new(n);
                    for i in 0..n {
                        for j in (i + 1)..n {
                            let d = (initial_dm.get(i, j) * 1000.0).round() / 1000.0;
                            rounded.set(i, j, d);
                        }
                    }
                    rounded
                } else {
                    let dndpre_offset_shift = dndpre_offset_shift(seq_type.is_nucleotide());
                    let mut shifted_matrix: Vec<Vec<i32>> = scoring.substitution_matrix
                        .iter()
                        .map(|row| row.iter().map(|&v| v + dndpre_offset_shift).collect())
                        .collect();
                    let nscored = scoring.nscoredalphabets;
                    for i in 0..shifted_matrix.len() {
                        for j in 0..shifted_matrix[i].len() {
                            if i >= nscored || j >= nscored {
                                shifted_matrix[i][j] = 0;
                            }
                        }
                    }
                    let penalty_dist = scoring.gap.open;
                    let raw = compute_distance_matrix_scoring(
                        &msa.sequences, &shifted_matrix,
                        &scoring.amino_map, penalty_dist,
                    );
                    // C's dndpre writes hat2 with `%.3f` precision
                    // (`io.c::write_hat2`). dvtditr then reads back these
                    // 3-decimal values. Round here so musclesupg sees the
                    // same distances dvtditr sees — needed for `--treeout`
                    // branch-length parity with C.
                    let n = raw.nseq;
                    let mut rounded = DistanceMatrix::new(n);
                    for i in 0..n {
                        for j in (i + 1)..n {
                            let d = (raw.get(i, j) * 1000.0).round() / 1000.0;
                            rounded.set(i, j, d);
                        }
                    }
                    rounded
                };
                // With `--treein`, C's dvtditr also loads the user tree
                // (`dvtditr.c:766-768` if(intree) ...
                // veryfastsupg_double_loadtree). Override the distance-based
                // rebuild so refinement runs against the same user-supplied
                // topology as the progressive pass.
                let topo = if let Some(ref t) = user_topo {
                    t.clone()
                } else {
                    musclesupg(&dm, self.cluster_method)
                };
                // C's `--treeout` writes the refinement tree built inside
                // `dvtditr` (not the progressive tbfast tree). Override
                // `final_progressive_topo` so the Newick we emit for
                // FFT-NS-i / *-INS-i modes matches what C writes.
                final_progressive_topo = Some(topo.clone());
                // C's mafft script (scripts/mafft:1512) sets `iteratelimit=254`
                // for BESTFIRST/BAATARI0, else `iteratelimit=16`, then caps
                // `iterate` to that. Mirror that strategy-dependent cap so
                // `--bestfirst --maxiterate 100` actually runs 100 best-move
                // iterations (C width 713) instead of stopping at 16 (rust
                // width 725 before this fix).
                let iterate_limit = if self.bestfirst { 254 } else { 16 };
                // `dvtditr.c:704-708`: `if( njob == 2 ) { weight = 0; niter = 1; }`
                // — a pair is refined exactly once, unweighted (uniform
                // weights are what `BranchWeights` already returns at 2).
                let iterate_limit = if nseq == 2 { 1 } else { iterate_limit };
                let capped_iterations = (*iterations).min(iterate_limit);
                // C's mafft script always passes -F (use_fft=1) to dvtditr
                // for refinement (scripts/mafft line 1531: rnaoptit=" -F "),
                // regardless of whether progressive alignment used FFT.
                // `--allowshift`: C's dvtditr gets `-Q 2.0` →
                // `penalty_shift_factor = 2.0` → `trywarp = 1`, with
                // `penalty_shift = (int)(penalty_shift_factor * penalty)`
                // (constants.c:318). The refinement `penalty` is the scaled
                // gap-open (`scoring.gap.open`). Pass it so the refinement
                // profile DP enables the warp/shift state (already ported in
                // `profile_align_imp_with_boundary` via `gap.shift`).
                let refine_shift = if self.allowshift {
                    let spfactor = self.shift_penalty_factor.unwrap_or(2.0);
                    Some((spfactor * scoring.gap.open as f64) as i32 as f64)
                } else {
                    None
                };
                let params = RefinementParams {
                    max_iterations: capped_iterations,
                    // C picks the refinement implementation on `nthread > 0`
                    // (`tditeration.c:1433`) and the two converge differently.
                    per_cycle_convergence: self.nthread > 0,
                    use_fft: true,
                    legacy_gap_cost: self.legacy_gap_cost,
                    shift: refine_shift,
                    // Multi-distance-class DP gates on `unalign_level > 0`, not
                    // on `--allowshift`: C `scripts/mafft:1436` enables the
                    // `-s #` (specificityconsideration) path whenever
                    // `unalignlevel != 0.0`, whether set via `--allowshift`
                    // (→ 0.8) or `--unalignlevel #` directly. Only the warp/
                    // shift DP (`refine_shift` above) is allowshift-specific
                    // (it needs `spfactor < 10`, which `--unalignlevel` alone
                    // does not set).
                    unalign_level: self.unalign_level,
                    minimum_weight: self.minimum_weight.unwrap_or(0.00001),
                    bestfirst: self.bestfirst,
                    ..Default::default()
                };
                // C `dvtditr.c:882` switches to segmented refinement (split
                // the alignment at high-conservation anchors, refine each
                // segment independently) whenever `constraint == 0` and
                // `bunkatsu != 0` — i.e. FFT-NS-i but not the *-INS-i modes,
                // which set `constraint=2` and stay single-segment.
                // `--seed` adds local homology constraints (constraint != 0),
                // so seeded FFT-NS-i must also stay on the single-segment
                // refinement path — gate on `local_hom.is_none()`.
                // Phase (2) of the C importance computation (see the long
                // comment at the progressive-phase `recompute_importance`
                // above): dvtditr RE-computes `importance` from the 3-decimal
                // `hat2` refinement tree before iterating. `topo` here is that
                // rounded tree (built from the hat2-rounded `dm` / dndpre
                // distances just above, or the user tree). Recompute
                // `local_hom`'s importance from it so refinement sees the
                // same constraint weights C's dvtditr does. The progressive
                // merge already consumed the full-precision-tree importance.
                // Skipped under `--treein` (user_topo) since both phases use
                // the same loaded tree, and for RNA modes (importance there
                // is not distance-tree-derived).
                if local_hom.is_some() && pairwise_for_constraints.is_some()
                    && user_topo.is_none() && !uses_rna_constraints
                {
                    let weights = mafft_tree::sequence_weights(&topo);
                    let seq_refs: Vec<&[u8]> = msa.sequences.iter()
                        .map(|s| s.as_slice()).collect();
                    if let Some(ref mut lh) = local_hom {
                        mafft_align::recompute_importance(lh, &seq_refs, &weights);
                    }
                }
                // `--skipiterate F`: C's `dvtditr -E $fixthreshold` →
                // `autosubalignment = F`. The function
                // `generatesubalignmentstable` (mltaln9.c:15330-15407)
                // walks the tree and identifies sub-alignment clusters
                // whose internal merges are all ≤ F. Two outcomes:
                //   - Whole tree below threshold → skip refinement
                //     entirely (`distfromtip[0] <= threshold`, returns 1).
                //   - Otherwise → sub-alignments are recorded, and
                //     `dvtditr.c:997-1006` marks topology branches that
                //     are STRICT SUBSETS of any sub-alignment as
                //     `skipthisbranch[step][side]=1` so refinement
                //     skips them.
                let (skip_refinement, skip_branches_vec) = if let Some(f) = self.skipiterate {
                    let (sub_alignments, all_below) =
                        mafft_tree::generate_subalignments_table(&topo, f);
                    if all_below {
                        eprintln!(
                            "\n#################################################################\n\
                             # WARNING: Iterative refinment was not done because you gave a\n\
                             # large --skipiterate value ({f:.3}).\n\
                             #################################################################\n"
                        );
                        (true, Vec::new())
                    } else {
                        // Build the per-(step, side) skip mask from
                        // sub-alignments. A branch is skipped iff its
                        // subtree (step.left or step.right) is a STRICT
                        // SUBSET of any sub-alignment cluster (C's
                        // `includemember && !samemember` at
                        // `dvtditr.c:997-1006`).
                        let nsteps = topo.steps.len();
                        let sub_sets: Vec<std::collections::BTreeSet<usize>> =
                            sub_alignments.iter().map(|s| s.iter().copied().collect()).collect();
                        let skip_branches: Vec<(bool, bool)> = topo.steps.iter().map(|step| {
                            let l: std::collections::BTreeSet<usize> = step.left.iter().copied().collect();
                            let r: std::collections::BTreeSet<usize> = step.right.iter().copied().collect();
                            let mut skip_l = false;
                            let mut skip_r = false;
                            for s in &sub_sets {
                                if !skip_l && l.is_subset(s) && l != *s { skip_l = true; }
                                if !skip_r && r.is_subset(s) && r != *s { skip_r = true; }
                                if skip_l && skip_r { break; }
                            }
                            (skip_l, skip_r)
                        }).collect();
                        let _ = (nsteps, sub_alignments);
                        (false, skip_branches)
                    }
                } else {
                    (false, Vec::new())
                };
                let params = RefinementParams {
                    skip_branches: skip_branches_vec,
                    ..params
                };
                // `--skipiterate F` small-F: skip-branches honored
                // only by the standard `iterative_refine` BAATARI2
                // loop. C's dvtditr does this regardless of the
                // segmented mode; for parity we route the small-F
                // case through the un-segmented refinement.
                let has_skip_branches = !params.skip_branches.is_empty();
                let use_segmented = matches!(self.mode, AlignmentMode::FftNsi { .. })
                    && local_hom.is_none()
                    && !has_skip_branches;
                if !skip_refinement {
                    if params.bestfirst {
                        // `--bestfirst` (C `parallelizationstrategy=BESTFIRST`):
                        // pick best-gain branch per iteration. Bypasses both
                        // the BAATARI2 walk and the FFT-segmented variant.
                        crate::refinement::bestfirst_refine(
                            &mut msa, &topo, &scoring, &params, local_hom.as_ref(),
                        );
                    } else if use_segmented {
                        crate::refinement::segmented_iterative_refine(
                            &mut msa, &topo, &scoring, &params, local_hom.as_ref(),
                        );
                    } else {
                        iterative_refine(
                            &mut msa, &topo, &scoring, &params, local_hom.as_ref(),
                        );
                    }
                }
            }
        }

        // `--reorder`: permute output to the C-equivalent reorder ordering.
        //
        // - Non-PartTree: tree-DFS over the final progressive guide tree
        //   (`tbfast.c:2928` calls `topolorderz` on the post-UPGMA topology).
        // - PartTree (`--parttree` / `--dpparttree`): C runs `splittbfast`
        //   TWICE (`scripts/mafft:2655` and `:2681`). CALL 1 uses raw 6-mer
        //   distances; CALL 2 passes `-Z` (`fromaln=1`) and recomputes
        //   distances via `naivepairscore11` on the aligned sequences. The
        //   final output order is the COMPOSITION:
        //     `final_order[k] = call1_order[call2_order[k]]`
        //   We mirror both passes to reach byte-identity.
        if self.reorder_output {
            let order: Option<Vec<usize>> = if use_parttree {
                let kind = if scoring.seq_type.is_nucleotide() {
                    PtSeqKind::Dna
                } else {
                    PtSeqKind::Protein
                };
                // CALL 1: parttree pivot pipeline on raw sequences.
                let call1_order = mafft_tree::parttree_split::compute_parttree_order(
                    &sequences, kind, 50,
                );
                // Reorder the FIRST-PASS aligned MSA into CALL 1's order so
                // CALL 2 sees `pre_1` (C's intermediate alignment), not the
                // final `pre_2`. Without using `first_pass_msa` here, our
                // CALL 2 distances would diverge from C's because the two
                // passes produce subtly different alignments.
                let source_msa: &Vec<Vec<u8>> = first_pass_msa
                    .as_ref().unwrap_or(&msa.sequences);
                let aligned_reordered: Vec<Vec<u8>> = call1_order
                    .iter().map(|&i| source_msa[i].clone()).collect();
                // CALL 2: parttree pivot pipeline with `fromaln=1` scoring
                // on the reordered aligned MSA. Uses the progressive-phase
                // substitution matrix and gap penalty (matches C's `penalty`
                // global set by `constants()`).
                let call2_order = mafft_tree::parttree_split::compute_parttree_order_fromaln(
                    &aligned_reordered,
                    &scoring.consweight_matrix,
                    &scoring.amino_map,
                    scoring.gap.open as f64,
                );
                // Compose: final_order[k] = call1_order[call2_order[k]].
                Some(call2_order.iter().map(|&k| call1_order[k]).collect())
            } else {
                final_progressive_topo.as_ref().map(|t| t.dfs_order())
            };
            if let Some(order) = order {
                if order.len() == msa.sequences.len() {
                    msa.sequences = order.iter().map(|&i| msa.sequences[i].clone()).collect();
                    msa.names = order.iter().map(|&i| msa.names[i].clone()).collect();
                }
            }
        }

        // Expose the final progressive guide tree to callers (used for
        // `--treeout` Newick serialization by the CLI binary).
        msa.guide_tree = final_progressive_topo;
        // Expose the first-pass alignment too — `--parttree --treeout`
        // and `--parttree --reorder` both need C MAFFT's `pre_1` to
        // reproduce CALL 2's tree / order generation.
        msa.first_pass_sequences = first_pass_msa;
        // Expose the (post-musclesupg) distance matrix used by the
        // progressive merge. `--distout` writes this to `<input>.hat2`,
        // and `--scoreout` derives the unweighted SP score from it.
        // Empty / None for paths that didn't compute a full dm
        // (PartTree, `--treein` with a user tree, etc.).
        if dm.nseq == nseq && dm.nseq > 0 {
            msa.distance_matrix = Some(dm.clone());
        }

        msa
    }

    /// Add new sequences to an existing alignment.
    ///
    /// `existing_input` is the already-aligned MSA (FASTA with gaps).
    /// `new_input` contains the new unaligned sequences to add.
    /// `keeplength` if true, preserves the existing alignment's column structure.
    pub fn add_to_alignment(
        &self,
        existing_input: &SequenceSet,
        new_input: &SequenceSet,
        keeplength: bool,
    ) -> MultipleAlignment {
        let seq_type = existing_input.seq_type;
        let scoring_model = if seq_type.is_nucleotide() {
            ScoringModel::Dna
        } else {
            self.scoring_model
        };

        let mut scoring = build_context(scoring_model, seq_type);

        if let Some(op) = self.gap_open {
            let ppenalty = -(op * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 3.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            scoring.gap.open = (scale * ppenalty as f64 + 0.5) as i32;
        }
        if let Some(ep) = self.gap_offset {
            let poffset = -(ep * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 1.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            let new_offset = (scale * poffset as f64 + 0.5) as i32;
            let matrix_offset = 0i32;
            let delta = new_offset - matrix_offset;
            if delta != 0 {
                let nscored = scoring.nscoredalphabets;
                for i in 0..nscored {
                    for j in 0..nscored {
                        scoring.substitution_matrix[i][j] -= delta;
                        scoring.consweight_matrix[i][j] = scoring.substitution_matrix[i][j] as f64;
                        scoring.fft_matrix[i][j] = scoring.substitution_matrix[i][j] + new_offset;
                    }
                }
            }
            scoring.gap.offset = new_offset;
        }

        let use_fft = !self.nofft && matches!(
            self.mode,
            AlignmentMode::FftNs2 | AlignmentMode::FftNsi { .. }
        );

        let existing = MultipleAlignment {
            sequences: existing_input.sequences.iter().map(|s| s.data.clone()).collect(),
            names: existing_input.sequences.iter().map(|s| s.name.clone()).collect(),
            score: 0.0,
            step_trace: Vec::new(), guide_tree: None, first_pass_sequences: None, distance_matrix: None,
        };

        let new_sequences: Vec<Vec<u8>> = new_input.sequences.iter().map(|s| s.data.clone()).collect();
        let new_names: Vec<String> = new_input.sequences.iter().map(|s| s.name.clone()).collect();

        if keeplength {
            add_sequences_keeplength(&existing, &new_sequences, &new_names, &scoring, use_fft)
        } else {
            add_sequences(&existing, &new_sequences, &new_names, &scoring, use_fft)
        }
    }

    /// Same as [`add_to_alignment`] with `keeplength = true`, but also
    /// returns the per-added-sequence list of dropped insertion runs.
    /// Used by `--mapout` / `--compactmapout` to emit the `.map`
    /// file. Each entry is `(start_pos_in_addbk_0based, run_length)`.
    pub fn add_to_alignment_with_map(
        &self,
        existing_input: &SequenceSet,
        new_input: &SequenceSet,
    ) -> (MultipleAlignment, Vec<Vec<(usize, usize)>>) {
        let seq_type = existing_input.seq_type;
        let scoring_model = if seq_type.is_nucleotide() {
            ScoringModel::Dna
        } else {
            self.scoring_model
        };
        let mut scoring = build_context(scoring_model, seq_type);
        if let Some(op) = self.gap_open {
            let ppenalty = -(op * 1000.0) as i32;
            let scale = if seq_type.is_nucleotide() { 3.0 * 600.0 / 1000.0 } else { 600.0 / 1000.0 };
            scoring.gap.open = (scale * ppenalty as f64 + 0.5) as i32;
        }
        let use_fft = !self.nofft && matches!(
            self.mode,
            AlignmentMode::FftNs2 | AlignmentMode::FftNsi { .. }
        );
        let existing = MultipleAlignment {
            sequences: existing_input.sequences.iter().map(|s| s.data.clone()).collect(),
            names: existing_input.sequences.iter().map(|s| s.name.clone()).collect(),
            score: 0.0,
            step_trace: Vec::new(), guide_tree: None, first_pass_sequences: None, distance_matrix: None,
        };
        let new_sequences: Vec<Vec<u8>> = new_input.sequences.iter().map(|s| s.data.clone()).collect();
        let new_names: Vec<String> = new_input.sequences.iter().map(|s| s.name.clone()).collect();
        crate::add::add_sequences_keeplength_with_map(
            &existing, &new_sequences, &new_names, &scoring, use_fft,
        )
    }

    /// Convenience: read FASTA file and align.
    pub fn align_file(&self, path: &std::path::Path) -> Result<MultipleAlignment, mafft_io::IoError> {
        let input = read_fasta(path)?;
        Ok(self.align(&input))
    }
}

/// Compute pairwise 6-tuple distances from raw (unaligned) sequences.
/// Matches C's default distance computation using `commonsextet_p`.
fn compute_distance_matrix_from_seqs(sequences: &[Vec<u8>]) -> DistanceMatrix {
    let nseq = sequences.len();
    let pairs: Vec<(usize, usize, f64)> = (0..nseq)
        .into_par_iter()
        .flat_map(|i| {
            let seqs = sequences;
            ((i + 1)..nseq).into_par_iter().map(move |j| {
                let d = ktuple_distance(&seqs[i], &seqs[j], 6);
                (i, j, d)
            })
        })
        .collect();

    let mut dm = DistanceMatrix::new(nseq);
    for (i, j, d) in pairs {
        dm.set(i, j, d);
    }
    dm
}

/// Compute pairwise distances from aligned sequences using scoring matrix.
///
/// Ports C's `msadistmtxthread` which uses `naivepairscorefast` for the
/// retree distance computation. This produces different distances from
/// simple identity distance and thus a different guide tree.
fn compute_distance_matrix_scoring(
    sequences: &[Vec<u8>],
    matrix: &[Vec<i32>],
    amino_map: &[u8; 256],
    penalty_dist: i32,
) -> DistanceMatrix {
    let nseq = sequences.len();
    let pairs: Vec<(usize, usize, f64)> = (0..nseq)
        .into_par_iter()
        .flat_map(|i| {
            let seqs = sequences;
            ((i + 1)..nseq).into_par_iter().map(move |j| {
                let d = scoring_matrix_distance(&seqs[i], &seqs[j], matrix, amino_map, penalty_dist);
                (i, j, d)
            })
        })
        .collect();

    let mut dm = DistanceMatrix::new(nseq);
    for (i, j, d) in pairs {
        dm.set(i, j, d);
    }
    dm
}


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

    /// C `constants.c:316-322` (nucleotide) vs `:672-677` (protein). The
    /// nucleotide `3 *` on the gap penalties is what keeps DNA L-INS-i /
    /// G-INS-i / E-INS-i byte-identical to C; the offset stays at `1 *`.
    /// C's `dndpre` default `poffset` differs by alphabet: `DEFAULTOFS_N`
    /// (`DNA.h:3`) vs `DEFAULTOFS_B` (`blosum.c:3`). Using the protein value
    /// for DNA silently reorders the refinement guide tree.
    /// C selects its refinement implementation on `nthread > 0`
    /// (`tditeration.c:1433`), and its script maps both *no* `--thread` and
    /// `--thread 0` to `dvtditr -C 0`. So 0 must mean the single-threaded
    /// convergence rule, and anything >= 1 the `athread` per-cycle rule.
    #[test]
    fn nthread_selects_the_convergence_rule_like_c() {
        let per_cycle = |nthread: usize| nthread > 0;
        assert!(!per_cycle(0), "no --thread / --thread 0 => C -C 0 => single-threaded rule");
        assert!(per_cycle(1), "--thread 1 => C -C 1 => athread rule");
        assert!(per_cycle(4));
        // Default engine must not opt into the athread rule.
        assert_eq!(MafftEngine::new(AlignmentMode::FftNs2).nthread, 0);
    }

    #[test]
    fn dndpre_offset_shift_mirrors_constants_c() {
        // offset = (int)( 600/1000 * poffset + 0.5 ); shift = -offset.
        let c_offset = |poffset: i32| (0.6 * poffset as f64 + 0.5) as i32;
        assert_eq!(c_offset(-369), -220, "nucleotide DEFAULTOFS_N");
        assert_eq!(c_offset(-123), -73, "protein DEFAULTOFS_B");
        assert_eq!(dndpre_offset_shift(true), 220, "DNA must not use the protein shift");
        assert_eq!(dndpre_offset_shift(false), 73);
    }

    #[test]
    fn pair_penalty_scales_mirror_constants_c() {
        let (gap, off) = pair_penalty_scales(true);
        assert_eq!(gap, 3.0 * 600.0 / 1000.0, "nucleotide gap scale must carry C's `3 *`");
        assert_eq!(off, 600.0 / 1000.0, "nucleotide offset scale is `1 *`, not `3 *`");
        let (gap, off) = pair_penalty_scales(false);
        assert_eq!(gap, 600.0 / 1000.0);
        assert_eq!(off, 600.0 / 1000.0);
        // Worked values for the L-INS-i defaults (lgop=-2.00 → ppenalty=-2000,
        // laof=0.100 → poffset=100), rounded the way C's `(int)(x + 0.5)` does:
        let cc = |pp: i32, sc: f64| ((sc * pp as f64) + 0.5) as i32;
        assert_eq!(cc(-2000, pair_penalty_scales(true).0), -3599);  // C: -3600+0.5 → -3599
        assert_eq!(cc(-2000, pair_penalty_scales(false).0), -1199);
        assert_eq!(cc(100, pair_penalty_scales(true).1), 60);
    }
    use mafft_types::{Sequence, SeqType};

    fn make_test_input() -> SequenceSet {
        SequenceSet {
            sequences: vec![
                Sequence { name: "s1".into(), data: b"ACDEFGHIKLMNP".to_vec() },
                Sequence { name: "s2".into(), data: b"ACDEFHIKLMNP".to_vec() },
                Sequence { name: "s3".into(), data: b"ACDEHIKLMNP".to_vec() },
            ],
            seq_type: SeqType::Protein,
        }
    }

    #[test]
    fn engine_fftns2() {
        let engine = MafftEngine::new(AlignmentMode::FftNs2);
        let msa = engine.align(&make_test_input());
        assert_eq!(msa.nseq(), 3);
        let w = msa.width();
        assert!(w >= 13);
        for seq in &msa.sequences { assert_eq!(seq.len(), w); }
    }

    #[test]
    fn engine_with_refinement() {
        let engine = MafftEngine::new(AlignmentMode::FftNsi { iterations: 5 });
        let msa = engine.align(&make_test_input());
        assert_eq!(msa.nseq(), 3);
        let w = msa.width();
        for seq in &msa.sequences { assert_eq!(seq.len(), w); }
    }

    #[test]
    fn engine_two_sequences() {
        let input = SequenceSet {
            sequences: vec![
                Sequence { name: "a".into(), data: b"ACDEFGHIK".to_vec() },
                Sequence { name: "b".into(), data: b"ACDEFGHIK".to_vec() },
            ],
            seq_type: SeqType::Protein,
        };
        let engine = MafftEngine::default();
        let msa = engine.align(&input);
        assert_eq!(msa.nseq(), 2);
        assert_eq!(msa.sequences[0], msa.sequences[1]);
    }

    #[test]
    fn engine_linsi_mode() {
        let engine = MafftEngine::new(AlignmentMode::LInsi { iterations: 2 });
        let msa = engine.align(&make_test_input());
        assert_eq!(msa.nseq(), 3);
        let w = msa.width();
        for seq in &msa.sequences { assert_eq!(seq.len(), w); }
    }

    #[test]
    fn engine_retree_1_vs_2() {
        let input = make_test_input();
        let msa1 = MafftEngine::new(AlignmentMode::FftNs2).with_retree(1).align(&input);
        let msa2 = MafftEngine::new(AlignmentMode::FftNs2).with_retree(2).align(&input);
        // Both should produce valid alignments
        assert_eq!(msa1.nseq(), 3);
        assert_eq!(msa2.nseq(), 3);
        let w1 = msa1.width();
        let w2 = msa2.width();
        for seq in &msa1.sequences { assert_eq!(seq.len(), w1); }
        for seq in &msa2.sequences { assert_eq!(seq.len(), w2); }
    }

    /// Guard: refinement tree uses scoring-matrix distance, not identity distance.
    ///
    /// C's dvtditr.c reads scoring-matrix-based distances from the hat2 file
    /// (written during the retree pass) to build the UPGMA tree for refinement.
    /// This test verifies that the engine calls `compute_distance_matrix_scoring`
    /// (not `compute_distance_matrix_from_alignment`) by checking that the
    /// refinement path in the source uses `scoring.substitution_matrix`.
    ///
    /// Functional check: FFT-NS-i with refinement produces a valid alignment
    /// on a dataset where scoring-matrix vs identity distance would yield
    /// different guide trees (sequences with varying conservation levels).
    #[test]
    fn engine_refinement_uses_scoring_matrix_distance() {
        let input = SequenceSet {
            sequences: vec![
                Sequence { name: "s1".into(), data: b"ACDEFGHIKLMNPQRSTVWY".to_vec() },
                Sequence { name: "s2".into(), data: b"ACDEFGHIKLMNPQRSTVWY".to_vec() },
                Sequence { name: "s3".into(), data: b"WWWWWWWWWWWWWWWWWWWW".to_vec() },
                Sequence { name: "s4".into(), data: b"ACDHIKLMNP".to_vec() },
                Sequence { name: "s5".into(), data: b"ACDEHIKLMNPQR".to_vec() },
            ],
            seq_type: SeqType::Protein,
        };
        let engine = MafftEngine::new(AlignmentMode::FftNsi { iterations: 5 });
        let msa = engine.align(&input);
        assert_eq!(msa.nseq(), 5);
        let w = msa.width();
        for (i, seq) in msa.sequences.iter().enumerate() {
            assert_eq!(seq.len(), w, "sequence {i} has wrong width after refinement");
            let residues = seq.iter().filter(|&&c| c != b'-').count();
            assert_eq!(residues, input.sequences[i].data.len(),
                "sequence {i} lost residues during refinement");
        }
    }

    /// Guard: engine passes cut=0.0 (default) to refinement params.
    ///
    /// Verifies the engine uses `..Default::default()` which has cut=0.0,
    /// not a hardcoded nonzero value.
    #[test]
    fn engine_refinement_params_use_default_cut() {
        // The RefinementParams default must have cut=0.0.
        // The engine constructs params with `..Default::default()`,
        // so this transitively guards the engine's behavior.
        let params = crate::refinement::RefinementParams::default();
        assert_eq!(params.cut, 0.0);
    }

}