fix(build): make graph construction invariant under input FASTA order

.gitignore

@@ -23,4 +23,7 @@
 Thumbs.db
 __pycache__/
 ehthumbs.db
-/debug
+/debug
+
+# compiled typst notes (source .typ is tracked)
+/notes/*.pdf

notes/assets/n00/dual_benchmark.csv

@@ -0,0 +1,21 @@
+dataset,variant,seconds,n_paths,n_blocks,pangenome_bp,core_bp,core_blocks
+example.fa,nodual,0.0,2,1,1000,1000,1
+example.fa,dual,0.0,2,1,1000,1000,1
+russian_doll_plasmids.fa.gz,nodual,0.2,4,11,50790,33610,5
+russian_doll_plasmids.fa.gz,dual,0.2,4,11,50790,33610,5
+flu-h3.fa,nodual,0.2,51,1,1737,1737,1
+flu-h3.fa,dual,0.2,51,1,1737,1737,1
+flu-h1.fa,nodual,0.8,171,10,2332,762,1
+flu-h1.fa,dual,0.9,171,10,2332,762,1
+ges-1.fa,nodual,7.8,33,218,596005,1016,1
+ges-1.fa,dual,8.2,33,218,596005,1016,1
+mpox.fa,nodual,19.3,13,63,207072,165129,19
+mpox.fa,dual,21.9,13,63,207072,165129,19
+sc2.fa,nodual,14.4,169,404,145717,1094,3
+sc2.fa,dual,17.0,169,404,145717,1094,3
+campylobacter-3.fa,nodual,17.3,3,512,3445626,102659,49
+campylobacter-3.fa,dual,21.8,3,514,3445993,102867,50
+klebs.fa.gz,nodual,142,9,1381,7644985,4458249,268
+klebs.fa.gz,dual,157,9,1372,7644453,4457866,268
+ecoli.fa.gz,nodual,157,10,2914,7830290,3782006,498
+ecoli.fa.gz,dual,187,10,2908,7827250,3782120,498

packages/pangraph/src/align/block_names.rs

@@ -0,0 +1,203 @@
+use crate::make_internal_report;
+use crate::pangraph::pangraph_block::{BlockId, PangraphBlock};
+use eyre::Report;
+use std::collections::{BTreeMap, HashMap};
+use std::hash::{Hash, Hasher};
+use twox_hash::XxHash64;
+
+/// Number of hex digits used to render the consensus hash in a canonical name.
+const HASH_WIDTH: usize = 16;
+
+/// Number of decimal digits used to render the block id in a canonical name.
+/// `usize::MAX` is 20 digits, so this never truncates.
+const ID_WIDTH: usize = 20;
+
+/// Content-derived ordering key for a block: hash of its consensus sequence, tie-broken by block id.
+///
+/// The tie-break only engages for blocks whose consensus is byte-identical, in which case the
+/// alignment between them is symmetric and the choice cannot bias the result.
+pub type BlockKey = (u64, BlockId);
+
+/// Content-derived naming and ordering for the blocks handed to an alignment backend.
+///
+/// `BlockId`s are assigned from the position of a record in the input (see `Pangraph::singleton`),
+/// so using them to name sequences leaks the input order into the aligner. `minimap2` run with
+/// `-X` keeps only one direction of each pair, chosen by `strcmp` of the sequence names, which
+/// makes the query/reference roles - and therefore the merged consensus - depend on the order the
+/// FASTA files were listed on the command line.
+///
+/// This type replaces those names with a key derived from the block consensus, so that both the
+/// order blocks are fed to the aligner and the names they carry are a pure function of the block
+/// content. Names are zero-padded to a fixed width, so byte-wise (`strcmp`) comparison agrees with
+/// the numeric `(consensus_hash, block_id)` order that [`BlockNames::canonical_order`] uses.
+pub struct BlockNames {
+  names: BTreeMap<BlockId, String>,
+  ids: HashMap<String, BlockId>,
+  keys: BTreeMap<BlockId, BlockKey>,
+  order: Vec<BlockId>,
+}
+
+/// Hashes a block's consensus sequence.
+///
+/// Content-derived, and therefore independent of how `BlockId`s were assigned.
+pub fn consensus_hash(block: &PangraphBlock) -> u64 {
+  let mut hasher = XxHash64::with_seed(0);
+  block.consensus().hash(&mut hasher);
+  hasher.finish()
+}
+
+/// Computes the content-derived ordering key of a block.
+fn consensus_key(id: BlockId, block: &PangraphBlock) -> BlockKey {
+  (consensus_hash(block), id)
+}
+
+/// Renders a canonical, fixed-width name for a block.
+fn canonical_name((hash, id): BlockKey) -> String {
+  format!("{hash:0HASH_WIDTH$x}_{:0ID_WIDTH$}", id.0)
+}
+
+impl BlockNames {
+  /// Builds the canonical naming and ordering for a set of blocks.
+  pub fn from_blocks(blocks: &BTreeMap<BlockId, PangraphBlock>) -> Self {
+    let keys: BTreeMap<BlockId, BlockKey> = blocks
+      .iter()
+      .map(|(&id, block)| (id, consensus_key(id, block)))
+      .collect();
+
+    let mut order: Vec<BlockId> = keys.keys().copied().collect();
+    order.sort_unstable_by_key(|id| keys[id]);
+
+    let names: BTreeMap<BlockId, String> = keys.iter().map(|(&id, &key)| (id, canonical_name(key))).collect();
+    let ids: HashMap<String, BlockId> = names.iter().map(|(&id, name)| (name.clone(), id)).collect();
+
+    Self {
+      names,
+      ids,
+      keys,
+      order,
+    }
+  }
+
+  /// Blocks in canonical (content-derived) order.
+  pub fn canonical_order(&self) -> impl Iterator<Item = BlockId> + '_ {
+    self.order.iter().copied()
+  }
+
+  /// Canonical name of a block, as handed to the alignment backend.
+  pub fn name(&self, id: BlockId) -> &str {
+    &self.names[&id]
+  }
+
+  /// Recovers the block a canonical name refers to.
+  ///
+  /// Fails if the aligner returned a name we never fed it, which would otherwise surface as a
+  /// silently wrong block id.
+  pub fn id_of(&self, name: &str) -> Result<BlockId, Report> {
+    self
+      .ids
+      .get(name)
+      .copied()
+      .ok_or_else(|| make_internal_report!("Aligner returned unknown sequence name '{name}'"))
+  }
+
+  /// Content-derived sort key of a block, for order-independent tie-breaking.
+  pub fn sort_key(&self, id: BlockId) -> BlockKey {
+    self.keys[&id]
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+  use crate::pangraph::pangraph_node::NodeId;
+  use crate::representation::seq::Seq;
+  use itertools::Itertools;
+  use pretty_assertions::assert_eq;
+  use rstest::rstest;
+
+  /// Builds a block map from `(id, consensus)` pairs.
+  fn blocks_of(spec: &[(usize, &str)]) -> BTreeMap<BlockId, PangraphBlock> {
+    spec
+      .iter()
+      .map(|&(id, seq)| {
+        let bid = BlockId(id);
+        (bid, PangraphBlock::from_consensus(Seq::from_str(seq), bid, NodeId(id)))
+      })
+      .collect()
+  }
+
+  #[rstest]
+  fn canonical_names_are_fixed_width_and_unique() {
+    let blocks = blocks_of(&[(0, "ACGT"), (1, "TTTT"), (12345, "GGGG")]);
+    let names = BlockNames::from_blocks(&blocks);
+
+    let rendered = blocks.keys().map(|&id| names.name(id).to_owned()).collect_vec();
+    assert!(rendered.iter().all(|n| n.len() == HASH_WIDTH + 1 + ID_WIDTH));
+    assert_eq!(rendered.iter().unique().count(), rendered.len());
+  }
+
+  #[rstest]
+  fn canonical_order_matches_lexicographic_name_order() {
+    let blocks = blocks_of(&[(0, "ACGT"), (1, "TTTT"), (2, "GGGG"), (3, "CCCC")]);
+    let names = BlockNames::from_blocks(&blocks);
+
+    let by_order = names.canonical_order().map(|id| names.name(id)).collect_vec();
+    let sorted = by_order.iter().copied().sorted().collect_vec();
+    assert_eq!(by_order, sorted);
+  }
+
+  /// The property the fix exists for: permuting block ids over the same consensus set must not
+  /// change the order or the relative naming the aligner sees.
+  #[rstest]
+  fn canonical_order_is_invariant_under_id_permutation() {
+    let seqs = ["ACGTACGT", "TTTTGGGG", "GGGGCCCC", "CACACACA"];
+
+    let forward = blocks_of(&seqs.iter().enumerate().map(|(i, s)| (i, *s)).collect_vec());
+    let reversed = blocks_of(
+      &seqs
+        .iter()
+        .enumerate()
+        .map(|(i, s)| (seqs.len() - 1 - i, *s))
+        .collect_vec(),
+    );
+
+    let n_fwd = BlockNames::from_blocks(&forward);
+    let n_rev = BlockNames::from_blocks(&reversed);
+
+    // The sequences appear in the same order regardless of how ids were assigned.
+    let seq_of = |blocks: &BTreeMap<BlockId, PangraphBlock>, id: BlockId| blocks[&id].consensus().as_str().to_owned();
+    let order_fwd = n_fwd.canonical_order().map(|id| seq_of(&forward, id)).collect_vec();
+    let order_rev = n_rev.canonical_order().map(|id| seq_of(&reversed, id)).collect_vec();
+    assert_eq!(order_fwd, order_rev);
+  }
+
+  #[rstest]
+  fn names_round_trip_to_block_ids() {
+    let blocks = blocks_of(&[(0, "ACGT"), (7, "TTTT")]);
+    let names = BlockNames::from_blocks(&blocks);
+
+    for &id in blocks.keys() {
+      assert_eq!(names.id_of(names.name(id)).unwrap(), id);
+    }
+  }
+
+  #[rstest]
+  fn unknown_name_is_rejected() {
+    let blocks = blocks_of(&[(0, "ACGT")]);
+    let names = BlockNames::from_blocks(&blocks);
+    let err = names.id_of("not-a-name").unwrap_err();
+    assert!(err.to_string().contains("unknown sequence name"));
+  }
+
+  /// Blocks sharing a consensus must still receive distinct names: identical names would make
+  /// `minimap2` treat the pair as a self-comparison and discard the diagonal anchors that carry
+  /// their (perfect) alignment, so the two blocks could never merge.
+  #[rstest]
+  fn identical_consensus_blocks_get_distinct_names() {
+    let blocks = blocks_of(&[(0, "ACGTACGT"), (1, "ACGTACGT")]);
+    let names = BlockNames::from_blocks(&blocks);
+
+    assert_ne!(names.name(BlockId(0)), names.name(BlockId(1)));
+    assert_eq!(names.sort_key(BlockId(0)).0, names.sort_key(BlockId(1)).0);
+  }
+}

packages/pangraph/src/align/minimap2_lib/align_with_minimap2_lib.rs

@@ -1,34 +1,45 @@
 use crate::align::alignment::{Alignment, Hit};
 use crate::align::alignment_args::AlignmentArgs;
+use crate::align::block_names::BlockNames;
 use crate::pangraph::pangraph_block::{BlockId, PangraphBlock};
 use crate::pangraph::strand::Strand;
 use crate::{make_error, make_internal_error};
 use eyre::{Report, WrapErr};
-use itertools::{Itertools, izip};
+use itertools::Itertools;
 use minimap2::{Minimap2Args, Minimap2Index, Minimap2Mapper, Minimap2Preset, Minimap2Result};
 use noodles::sam::record::Cigar;
 use num_traits::clamp_min;
 use rayon::prelude::*;
 use std::collections::BTreeMap;
 use std::str::FromStr;
 
+/// Aligns the consensus sequences of `blocks` against each other, all-vs-all.
+///
+/// Blocks are fed in canonical (content-derived) order under canonical names, so that neither the
+/// set of alignments nor the query/reference role of each pair depends on how `BlockId`s were
+/// assigned - and therefore on the order the input FASTA files were listed. See [`BlockNames`].
 pub fn align_with_minimap2_lib(
   blocks: &BTreeMap<BlockId, PangraphBlock>,
+  names: &BlockNames,
   params: &AlignmentArgs,
 ) -> Result<Vec<Alignment>, Report> {
-  let (names, seqs): (Vec<String>, Vec<&str>) = blocks
-    .iter()
-    .map(|(id, block)| (id.to_string(), block.consensus().as_str()))
+  let (seq_names, seqs): (Vec<&str>, Vec<&str>) = names
+    .canonical_order()
+    .map(|id| (names.name(id), blocks[&id].consensus().as_str()))
     .unzip();
 
-  let alns: Vec<Alignment> = align_with_minimap2_lib_impl(&seqs, &names, params)?;
+  let alns: Vec<Alignment> = align_with_minimap2_lib_impl(&seqs, &seq_names, &|n| names.id_of(n), params)?;
 
   Ok(alns)
 }
 
+/// Aligner mechanics, decoupled from how sequence names map back to blocks.
+///
+/// `resolve` recovers the [`BlockId`] a name refers to; the caller decides the naming scheme.
 fn align_with_minimap2_lib_impl(
   seqs: &[impl AsRef<str>],
   names: &[impl AsRef<str>],
+  resolve: &(dyn Fn(&str) -> Result<BlockId, Report> + Sync),
   params: &AlignmentArgs,
 ) -> Result<Vec<Alignment>, Report> {
   if names.len() != seqs.len() {
@@ -61,11 +72,15 @@ fn align_with_minimap2_lib_impl(
 
   let idx = Minimap2Index::new(&seqs, &names, &args)?;
 
-  let results: Vec<Minimap2Result> = izip!(&seqs, &names)
-    .par_bridge()
+  // `par_iter().zip()` over slices is an indexed parallel iterator, for which `collect` preserves
+  // input order by construction. `par_bridge()` does not guarantee order, which would leave the
+  // energy-sort tie-breaking in `filter_matches` at the mercy of thread scheduling.
+  let results: Vec<Minimap2Result> = seqs
+    .par_iter()
+    .zip(names.par_iter())
     .map_init(
       || Minimap2Mapper::new(&idx).unwrap(),
-      move |mapper, (seq, name)| {
+      |mapper, (seq, name)| {
         mapper
           .run_map(seq, name)
           .wrap_err_with(|| format!("When aligning sequence '{name}'"))
@@ -75,7 +90,7 @@ fn align_with_minimap2_lib_impl(
 
   let alns = results
     .into_iter()
-    .map(Alignment::from_minimap_paf_obj)
+    .map(|res| Alignment::from_minimap_paf_obj(res, resolve))
     .collect::<Result<Vec<Vec<_>>, Report>>()?
     .into_iter()
     .flatten()
@@ -86,20 +101,23 @@ fn align_with_minimap2_lib_impl(
 
 #[allow(clippy::multiple_inherent_impl)]
 impl Alignment {
-  pub fn from_minimap_paf_obj(res: Minimap2Result) -> Result<Vec<Self>, Report> {
+  pub fn from_minimap_paf_obj(
+    res: Minimap2Result,
+    resolve: &(dyn Fn(&str) -> Result<BlockId, Report> + Sync),
+  ) -> Result<Vec<Self>, Report> {
     let Minimap2Result { pafs, .. } = res;
     pafs
       .into_iter()
       .map(|paf| {
         if let Some(cg) = &paf.cg {
           Ok(Alignment {
             qry: Hit::new(
-              BlockId::from_str(&paf.q.name)?,
+              resolve(&paf.q.name)?,
               paf.q.len,
               (paf.q.start as usize, paf.q.end as usize),
             ),
             reff: Hit::new(
-              BlockId::from_str(&paf.t.name)?,
+              resolve(&paf.t.name)?,
               paf.t.len,
               (paf.t.start as usize, paf.t.end as usize),
             ),
@@ -183,7 +201,8 @@ mod tests {
       ..AlignmentArgs::default()
     };
 
-    let actual = align_with_minimap2_lib_impl(&seqs, &names, &params)?;
+    // Names here are the FASTA record ids, so resolve them as plain integers.
+    let actual = align_with_minimap2_lib_impl(&seqs, &names, &|n| BlockId::from_str(&n.to_owned()), &params)?;
 
     let expected = vec![Alignment {
       qry: Hit::new(BlockId(0), 998, (0, 996)),

packages/pangraph/src/align/mmseqs/align_with_mmseqs.rs

@@ -1,5 +1,6 @@
 use crate::align::alignment::Alignment;
 use crate::align::alignment_args::AlignmentArgs;
+use crate::align::block_names::BlockNames;
 use crate::align::mmseqs::paf::PafTsvRecord;
 use crate::io::fasta::FastaWriter;
 use crate::io::file::open_file_or_stdin;
@@ -15,8 +16,14 @@ use std::io::Read;
 use std::process::Command;
 use tempfile::Builder as TempDirBuilder;
 
+/// Aligns the consensus sequences of `blocks` against each other using mmseqs.
+///
+/// Like the minimap2 backend, blocks are written in canonical (content-derived) order under
+/// canonical names, so the result does not depend on how `BlockId`s were assigned. See
+/// [`BlockNames`].
 pub fn align_with_mmseqs(
   blocks: &BTreeMap<BlockId, PangraphBlock>,
+  names: &BlockNames,
   params: &AlignmentArgs,
 ) -> Result<Vec<Alignment>, Report> {
   // TODO: This uses a global resource - filesystem.
@@ -29,9 +36,9 @@ pub fn align_with_mmseqs(
 
   {
     let mut writer = FastaWriter::from_path(&input_path)?;
-    blocks
-      .iter()
-      .try_for_each(|(id, block)| writer.write(id.to_string(), &None, block.consensus()))?;
+    names
+      .canonical_order()
+      .try_for_each(|id| writer.write(names.name(id), &None, blocks[&id].consensus()))?;
   }
 
   let output_column_names = PafTsvRecord::fields_names().join(",");
@@ -66,7 +73,7 @@ pub fn align_with_mmseqs(
   let mut paf_str = String::new();
   open_file_or_stdin(&Some(output_path))?.read_to_string(&mut paf_str)?;
 
-  Alignment::from_paf_str(&paf_str)
+  Alignment::from_paf_str(&paf_str, &|n| names.id_of(n))
 }
 
 // FIXME: This test is failing after commit a62b19b018b4b2f9602bc75335d4ab5ddbc7abf5