gwseq_io/cram/

reader.rs

//! The CRAM reader.
//!
//! Four read paths over the same machinery, matching the BAM reader's: the
//! index turns a locus into slices, a slice decodes into BAM records, and the
//! filter decides which of them the caller asked for.
//!
//! # Why the unit of work is a slice
//!
//! A CRAM slice cannot be entered partway. `AP` is delta-coded from the
//! slice's own start, the entropy coders run the length of a block, and `NF`
//! reaches from one record to another within the slice — so ten thousand
//! records come back or none do. Unlike a BGZF chunk there is nothing smaller
//! to seek to, which makes the slice the unit of caching as well as of work:
//! two loci landing in one slice decode it once.
//!
//! # Two reads per query, and what is memoised
//!
//! §12: a slice cannot be decoded without its container's compression header,
//! and the index gives a slice's offset *relative to the end of its container
//! header* — so reaching a slice means reading its container header first.
//! Both the container header and the compression header parsed out of it are
//! memoised per container, because every slice in a container wants the same
//! ones and parsing a compression header is two maps of encodings.

use std::sync::Arc;

use bytes::Bytes;
use parking_lot::Mutex;

use crate::bam::header::{HeaderLine, SamHeader};
use crate::bam::record::{decode_block, BamRecord, EntryFilter};
use crate::bam::EntriesRequest;
use crate::error::{Error, Result};
use crate::genomic::{ChrMap, Locs};
use crate::parallel::Executor;
use crate::progress::ProgressTracker;
use crate::source::ByteSource;

use super::compression::CompressionHeader;
use super::container::{
    Block, BlockContentType, ContainerHeader, FileDefinition, EOF_CONTAINER, FILE_DEFINITION_SIZE,
};
use super::crai::{CramIndex, IndexEntry};
use super::record::{decode_slice, ReferenceBases, References};
use super::reference::{resolve, Reference, ReferenceSource};
use super::slice::Slice;

/// Container headers kept parsed. A container is a few hundred kilobytes, so
/// this is a handful of megabytes of reach for a few kilobytes of state.
const CONTAINER_CACHE: usize = 8;

/// Decoded slices kept, in bytes of BAM records. A slice of ten thousand
/// hundred-base reads is roughly 2 MB, so this holds a few dozen.
const SLICE_CACHE_BYTES: usize = 64 << 20;

#[derive(Debug)]
struct Inner {
    source: Arc<dyn ByteSource>,
    executor: Executor,
    /// The `.crai`, when one was found beside the file.
    index: Option<Arc<CramIndex>>,
    /// One built by walking container headers, when there was no `.crai`.
    /// Built on first use rather than at open, so opening a file for its header
    /// costs one read.
    built: Mutex<Option<Arc<CramIndex>>>,
    reference: Option<ReferenceSource>,
    caches: Mutex<Caches>,
    /// Woken whenever a slice decode finishes, so a worker that found one
    /// already in flight can take the result instead of decoding it again.
    slice_ready: parking_lot::Condvar,
}

/// What a reader holds between calls.
#[derive(Debug, Default)]
struct Caches {
    /// `(container offset, header, compression header)`, most recent last.
    containers: Vec<(u64, Arc<ContainerHeader>, Arc<CompressionHeader>)>,
    /// `(container offset, landmark, records)`, most recent last.
    slices: Vec<(u64, u64, Bytes)>,
    slice_bytes: usize,
    /// Slices being decoded right now, by `(container offset, landmark)`.
    ///
    /// Without this, `parallel` workers that miss the same slice at the same
    /// moment all decode it and all insert a copy — each one counted against
    /// the cache's byte budget, so the duplicates then evict what the next
    /// locus needs. Adjacent loci share slices routinely, so this is the
    /// common case rather than a race worth ignoring.
    in_flight: Vec<(u64, u64)>,
}

pub struct CramReader {
    inner: Option<Inner>,
    path: String,
    index_path: String,
    header: SamHeader,
    chr_map: ChrMap,
    /// Reference names by index, shared with every record.
    chr_names: Arc<Vec<String>>,
    /// `@RG` identifiers by their position in the header, for rebuilding the
    /// `RG` tag the encoder was free to drop.
    read_groups: Vec<String>,
    version: (u8, u8),
    index_error: String,
    indexed: bool,
    reference: Reference,
    reference_error: String,
}

impl std::fmt::Debug for CramReader {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("CramReader")
            .field("path", &self.path)
            .field("version", &self.version)
            .field("references", &self.chr_map.len())
            .field("reference", &self.reference.path())
            .field("closed", &self.is_closed())
            .finish()
    }
}

impl CramReader {
    pub fn open(
        path: &str,
        index_path: Option<&str>,
        reference: Option<&str>,
        parallel: i64,
        block_size: Option<u64>,
        max_blocks: Option<usize>,
    ) -> Result<Self> {
        let source = crate::source::open(path, block_size, max_blocks)?;
        Self::from_source(
            source, path, index_path, reference, parallel, block_size, max_blocks,
        )
    }

    pub(crate) fn from_source(
        source: Arc<dyn ByteSource>,
        path: &str,
        index_path: Option<&str>,
        reference: Option<&str>,
        parallel: i64,
        block_size: Option<u64>,
        max_blocks: Option<usize>,
    ) -> Result<Self> {
        let definition = FileDefinition::parse(&source.read_at(0, FILE_DEFINITION_SIZE)?, path)?;
        let (header, chr_map, read_groups) = read_header(source.as_ref(), path)?;

        let mut names = vec![String::new(); chr_map.iter().map(|e| e.index + 1).max().unwrap_or(0)];
        for entry in chr_map.iter() {
            names[entry.index] = entry.id.clone();
        }

        // The index, as the BAM reader treats its own: absent is no error, and
        // unreadable is a different failure worth reporting separately.
        let index_path_given = index_path;
        let index_path = index_path
            .map(str::to_string)
            .unwrap_or_else(|| format!("{path}.crai"));
        // A named index that is not there is the caller's mistake and is said
        // so; the default `<file>.crai` simply not existing is not, and falls
        // through to building one. The BAM reader draws the same line.
        let named = index_path_given.is_some();
        let (index, index_error) =
            if !crate::source::is_url(&index_path) && !std::path::Path::new(&index_path).exists() {
                (
                    None,
                    if named {
                        format!("{index_path} is not there")
                    } else {
                        String::new()
                    },
                )
            } else {
                match crate::source::open(&index_path, block_size, max_blocks)
                    .and_then(|s| s.read_to_end(0))
                    .and_then(|data| CramIndex::parse(&data, &index_path))
                {
                    Ok(index) => (Some(Arc::new(index)), String::new()),
                    Err(e) => (None, e.to_string()),
                }
            };

        // §9's EOF container is the only thing that says a file is whole.
        // Without it the file was cut, and the symptom is a query that finds
        // nothing where it should have found the missing tail — reported here
        // as information, the way `samtools` warns "EOF marker is absent",
        // rather than as a refusal, because everything before the cut still
        // reads.
        let truncated = match source.len() {
            Ok(length) if length >= EOF_CONTAINER.len() as u64 => {
                let at = length - EOF_CONTAINER.len() as u64;
                !matches!(source.read_at(at, EOF_CONTAINER.len()), Ok(tail) if tail[..] == EOF_CONTAINER[..])
            }
            _ => true,
        };

        // The reference. A file that cannot find one still opens and still
        // answers everything but `sequence`, which is what `reference_error`
        // is for — the same shape `index_error` already has.
        // One reference window per worker, at least: `parallel` workers on
        // different loci otherwise evict each other's windows on every locus.
        let slots = if parallel <= 0 {
            std::thread::available_parallelism().map_or(4, |n| n.get())
        } else {
            parallel as usize
        };
        let sequences = sequence_details(&header);
        let mut resolved = resolve(reference, &sequences);
        let (mut reference_source, mut reference_error) = match &resolved {
            // The cache is a directory of bare sequences named by checksum,
            // not a FASTA, so it needs the `@SQ` lines to know which file is
            // which sequence.
            Reference::FromCache(root) => (
                Some(ReferenceSource::open_cache(root, &sequences).with_slots(slots)),
                String::new(),
            ),
            Reference::Given(path) | Reference::FromHeader(path) => {
                match ReferenceSource::open(path) {
                    Ok(source) => (Some(source.with_slots(slots)), String::new()),
                    Err(e) => (None, e.to_string()),
                }
            }
            Reference::None(why) => (None, why.clone()),
        };

        // A reference that opened is not yet a reference that *fits*. Without
        // this check, pointing at the wrong FASTA gives a populated
        // `reference_path`, an empty `reference_error` and every base `N` —
        // which is exactly what a read with no substitutions looks like, so a
        // caller has no way to learn the reference was never used. `samtools`
        // fails the read outright.
        if let Some(source) = &reference_source {
            if let Some(complaint) = reference_complaint(source, &chr_map, &sequences) {
                reference_source = None;
                // `reference_path` has to go with it: a path *and* an error is
                // a contradiction, and a caller reading only the path would
                // conclude a reference was used.
                resolved = Reference::None(complaint.clone());
                reference_error = complaint;
            }
        }

        // A `.crai` beside the file means indexed; so does a local file, whose
        // index can be built from its container headers on demand. A remote
        // file with no `.crai` cannot, because building means a request per
        // container.
        let indexed = index.is_some() || !crate::source::is_url(path);
        // A local file whose `.crai` would not parse still has an index — the
        // one `build` makes from its container headers — so the parse failure
        // is information rather than a refusal. Reported and then set aside;
        // before this it was returned from every read, which made a corrupt
        // index worse than no index at all.
        let index_error = if truncated && index_error.is_empty() {
            format!(
                "{path} has no EOF container, so it was cut short; anything past the cut is \
                 missing and a locus there will find nothing"
            )
        } else if index.is_none() && !index_error.is_empty() && !crate::source::is_url(path) {
            format!("{index_error} — reading it from the file's containers instead")
        } else if indexed || !index_error.is_empty() {
            index_error
        } else {
            format!(
                "{index_path} is not there, and this file is remote, so an index cannot be \
                 built by walking its containers"
            )
        };

        Ok(Self {
            inner: Some(Inner {
                source,
                executor: Executor::new(parallel)?,
                index,
                built: Mutex::new(None),
                reference: reference_source,
                caches: Mutex::new(Caches::default()),
                slice_ready: parking_lot::Condvar::new(),
            }),
            path: path.to_string(),
            index_path,
            header,
            chr_map,
            chr_names: Arc::new(names),
            read_groups,
            version: (definition.major, definition.minor),
            index_error,
            indexed,
            reference: resolved,
            reference_error,
        })
    }

    pub fn header(&self) -> &SamHeader {
        &self.header
    }
    pub fn chr_sizes(&self) -> &ChrMap {
        &self.chr_map
    }
    pub fn index_error(&self) -> &str {
        &self.index_error
    }
    pub fn is_indexed(&self) -> bool {
        self.indexed
    }
    pub fn is_closed(&self) -> bool {
        self.inner.is_none()
    }
    pub fn path(&self) -> &str {
        &self.path
    }
    pub fn parallel(&self) -> usize {
        self.inner.as_ref().map_or(0, |i| i.executor.parallel())
    }
    /// The CRAM version, as `(major, minor)`.
    pub fn version(&self) -> (u8, u8) {
        self.version
    }
    /// The reference in use, or `None` when none resolved.
    pub fn reference_path(&self) -> Option<&str> {
        self.reference.path()
    }
    /// Why there is no reference, when there is none. Empty when one resolved.
    pub fn reference_error(&self) -> &str {
        &self.reference_error
    }
    pub fn close(&mut self) {
        if let Some(inner) = self.inner.take() {
            inner.source.close();
        }
    }

    fn inner(&self) -> Result<&Inner> {
        self.inner.as_ref().ok_or_else(|| Error::Closed {
            path: self.path.clone(),
        })
    }

    /// The index: the `.crai`, or one built from the container headers.
    fn index(&self, inner: &Inner) -> Result<Arc<CramIndex>> {
        if let Some(index) = &inner.index {
            return Ok(index.clone());
        }
        // Only a remote file has nothing to fall back on: building an index
        // means a request per container. A local file with an unreadable
        // `.crai` builds one, and `index_error` says that is what happened.
        if !self.index_error.is_empty() && crate::source::is_url(&self.path) {
            return Err(Error::invalid(format!(
                "cram index {} could not be read: {}",
                self.index_path, self.index_error
            )));
        }
        let mut built = inner.built.lock();
        if let Some(index) = built.as_ref() {
            return Ok(index.clone());
        }
        let index = Arc::new(CramIndex::build(inner.source.as_ref())?);
        *built = Some(index.clone());
        Ok(index)
    }

    /// One list per locus, in the order the loci were given.
    pub fn read_entries(&self, req: &EntriesRequest) -> Result<Vec<Vec<BamRecord>>> {
        let inner = self.inner()?;
        let resolved = self.resolve(&req.locs)?;
        let coverage = resolved.iter().map(|(_, s, e)| (e - s).max(0) as u64).sum();
        let tracker = ProgressTracker::with_callback(coverage, req.progress.clone());
        let index = self.index(inner)?;

        // Loci are spread over the workers in runs rather than round-robin, so
        // each worker walks the request forward and the slices it decoded for
        // one locus serve the ones after it.
        let workers = inner.executor.parallel().min(resolved.len().max(1)).max(1);
        let per_worker = resolved.len().div_ceil(workers).max(1);
        let batches: Vec<(usize, usize)> = (0..workers)
            .map(|w| (w * per_worker, ((w + 1) * per_worker).min(resolved.len())))
            .filter(|(from, to)| from < to)
            .collect();

        let lists = inner.executor.map_batches(&batches, |_, (from, to)| {
            let mut out = Vec::with_capacity(to - from);
            for locus in &resolved[*from..*to] {
                out.push(self.read_locus(inner, &index, *locus, req)?);
                tracker.add((locus.2 - locus.1).max(0) as u64);
            }
            Ok(out)
        })?;
        tracker.done_report();
        Ok(lists.into_iter().flatten().collect())
    }

    /// Every alignment on the named references, in reference then coordinate
    /// order.
    pub fn read_all_entries(&self, req: &EntriesRequest) -> Result<Vec<BamRecord>> {
        let locs = Locs::whole_chromosomes(&self.chr_map, &req.locs.chr_ids)?;
        let whole = EntriesRequest {
            locs,
            ..req.clone()
        };
        Ok(self.read_entries(&whole)?.into_iter().flatten().collect())
    }

    /// Per-locus, lazily.
    pub fn iter_entries(&self, req: &EntriesRequest) -> Result<LocusEntries<'_>> {
        LocusEntries::plan(self, req)
    }

    /// Successive windows over whole references.
    pub fn iter_all_entries(&self, req: &EntriesRequest, window: i64) -> Result<WindowEntries<'_>> {
        WindowEntries::plan(self, req, window)
    }

    fn resolve(&self, locs: &Locs) -> Result<Vec<(usize, i64, i64)>> {
        (0..locs.len())
            .map(|i| {
                let entry = self.chr_map.resolve(&locs.chr_ids[i])?;
                Ok((entry.index, locs.starts[i], locs.ends[i]))
            })
            .collect()
    }

    /// The alignments of one locus.
    fn read_locus(
        &self,
        inner: &Inner,
        index: &CramIndex,
        locus: (usize, i64, i64),
        req: &EntriesRequest,
    ) -> Result<Vec<BamRecord>> {
        let (chr, start, end) = locus;
        let filter = EntryFilter {
            chr_index: Some(chr as i32),
            start,
            end: Some(end),
            standard_flags: req.filter.enabled,
        };
        let mut out = Vec::new();
        for entry in index.slices(chr as i32, start, end) {
            let records = self.slice_records(inner, &entry)?;
            out.extend(decode_block(
                &records,
                req.parse_tags,
                &filter,
                &self.chr_names,
                &self.path,
            )?);
        }
        Ok(out)
    }

    /// One slice's records as BAM bytes, from the cache or from the file.
    fn slice_records(&self, inner: &Inner, entry: &IndexEntry) -> Result<Bytes> {
        let key = (entry.container_offset, entry.landmark);
        {
            let mut caches = inner.caches.lock();
            loop {
                if let Some(position) = caches.slices.iter().position(|(c, l, _)| (*c, *l) == key) {
                    // Move to the back, so the cache evicts what has gone
                    // longest unused rather than what was read longest ago.
                    let hit = caches.slices.remove(position);
                    let records = hit.2.clone();
                    caches.slices.push(hit);
                    return Ok(records);
                }
                if !caches.in_flight.contains(&key) {
                    caches.in_flight.push(key);
                    break;
                }
                // Another worker is decoding this one. Wait for it rather than
                // doing the same work: a slice decode is milliseconds and the
                // whole point of the cache.
                inner.slice_ready.wait(&mut caches);
            }
        }
        // From here the marker is ours and must come off whichever way this
        // goes, or every later reader of this slice waits forever.
        let outcome = self.decode_one_slice(inner, entry);
        let mut caches = inner.caches.lock();
        caches.in_flight.retain(|held| *held != key);
        if let Ok(records) = &outcome {
            // Never twice: a key already present would be counted against the
            // byte budget twice over.
            if !caches.slices.iter().any(|(c, l, _)| (*c, *l) == key) {
                caches.slice_bytes += records.len();
                caches.slices.push((key.0, key.1, records.clone()));
                while caches.slice_bytes > SLICE_CACHE_BYTES && caches.slices.len() > 1 {
                    let (_, _, dropped) = caches.slices.remove(0);
                    caches.slice_bytes -= dropped.len();
                }
            }
        }
        drop(caches);
        inner.slice_ready.notify_all();
        outcome
    }

    /// Read and decode one slice, with no cache involved.
    fn decode_one_slice(&self, inner: &Inner, entry: &IndexEntry) -> Result<Bytes> {
        let (container, compression) = self.container(inner, entry.container_offset)?;
        let offset = container.blocks_offset() + entry.landmark;
        let slice = Slice::read(inner.source.as_ref(), offset, entry.size)?;

        // A slice on one reference gets its window fetched once, up front. A
        // multi-reference slice gets the source itself, and looks a reference
        // up per record — see `References`.
        let window = self.reference_window(inner, &slice)?;
        let references = match (&window, &inner.reference) {
            (Some((bases, start)), _) => References::Fixed(ReferenceBases {
                bases,
                // 1-based, as the record decoder works in.
                start: start + 1,
            }),
            (None, Some(source)) if slice.header.is_multi_ref() => References::ByRefId {
                source,
                names: &self.chr_names,
            },
            _ => References::None,
        };
        decode_slice(
            &slice,
            &compression,
            references,
            &self.read_groups,
            &self.path,
        )
    }

    /// The container header at `offset` and its compression header, memoised.
    fn container(
        &self,
        inner: &Inner,
        offset: u64,
    ) -> Result<(Arc<ContainerHeader>, Arc<CompressionHeader>)> {
        {
            let mut caches = inner.caches.lock();
            if let Some(position) = caches.containers.iter().position(|(at, ..)| *at == offset) {
                let hit = caches.containers.remove(position);
                let out = (hit.1.clone(), hit.2.clone());
                caches.containers.push(hit);
                return Ok(out);
            }
        }

        let header = ContainerHeader::read(inner.source.as_ref(), offset)?;
        // §5: the first block of a container is its compression header.
        let first = header
            .landmarks
            .first()
            .copied()
            .map(|landmark| landmark.max(0) as usize)
            .unwrap_or(header.length.max(0) as usize);
        let data = inner
            .source
            .read_exact_at(header.blocks_offset(), first.max(1))?;
        let block = Block::parse(&data, header.blocks_offset(), &self.path)?;
        if block.content_type != BlockContentType::CompressionHeader {
            return Err(Error::corrupt(
                &self.path,
                header.blocks_offset(),
                format!(
                    "the first block of a container is {:?}, not its compression header",
                    block.content_type
                ),
            ));
        }
        let compression = Arc::new(CompressionHeader::parse(&block.data, &self.path)?);
        let header = Arc::new(header);

        let mut caches = inner.caches.lock();
        caches
            .containers
            .push((offset, header.clone(), compression.clone()));
        while caches.containers.len() > CONTAINER_CACHE {
            caches.containers.remove(0);
        }
        Ok((header, compression))
    }

    /// The one window a single-reference slice needs, fetched before the
    /// decode.
    ///
    /// `None` means there is no single window to fetch — the slice is unmapped,
    /// or it is multi-reference and its records are resolved one at a time.
    fn reference_window(
        &self,
        inner: &Inner,
        slice: &Slice,
    ) -> Result<Option<(Arc<Vec<u8>>, i64)>> {
        // A slice carrying its own reference needs nothing from outside.
        if let Some(embedded) = &slice.embedded_reference {
            let start = slice.header.range().map(|(from, _)| from).unwrap_or(0);
            return Ok(Some((Arc::new(embedded.to_vec()), start)));
        }
        let Some(reference) = &inner.reference else {
            return Ok(None);
        };
        let Some((start, end)) = slice.header.range() else {
            return Ok(None);
        };
        let Some(name) = self.chr_names.get(slice.header.ref_id.max(0) as usize) else {
            return Ok(None);
        };
        if !reference.has(name) {
            return Ok(None);
        }
        Ok(Some(reference.window(name, start, end)?))
    }
}

/// Read the CRAM header container: a SAM text header, and the references and
/// read groups read out of it.
///
/// Unlike BAM, there is no binary reference list — the `@SQ` lines are the
/// only statement of what reference a record's id means.
fn read_header(source: &dyn ByteSource, path: &str) -> Result<(SamHeader, ChrMap, Vec<String>)> {
    let container = ContainerHeader::read(source, FILE_DEFINITION_SIZE as u64)?;
    let data = source.read_exact_at(container.blocks_offset(), container.length.max(0) as usize)?;
    let block = Block::parse(&data, container.blocks_offset(), path)?;
    if block.content_type != BlockContentType::FileHeader {
        return Err(Error::corrupt(
            path,
            container.blocks_offset(),
            format!(
                "the first container holds a {:?} block where its header should be",
                block.content_type
            ),
        ));
    }
    // §8.3: the block opens with a 32-bit length and then the SAM text.
    if block.data.len() < 4 {
        return Err(Error::corrupt(
            path,
            container.blocks_offset(),
            "a header block too short to hold its own length",
        ));
    }
    let length =
        i32::from_le_bytes(block.data[..4].try_into().expect("four bytes")).max(0) as usize;
    let text = &block.data[4..(4 + length).min(block.data.len())];
    let header = SamHeader::parse(&String::from_utf8_lossy(text));

    let mut entries = Vec::new();
    let mut read_groups = Vec::new();
    for line in &header.lines {
        match line.kind.as_str() {
            "SQ" => {
                // Reference ids are positional over the `@SQ` lines, so
                // skipping one renumbers every reference after it and every
                // record then names the wrong chromosome. `LN` is mandatory in
                // SAM, so a line without it is a broken header rather than a
                // line to step over quietly.
                let name = field(line, "SN");
                let size = field(line, "LN").and_then(|v| v.parse::<i64>().ok());
                match (name, size) {
                    (Some(name), Some(size)) => entries.push((name.to_string(), size)),
                    (name, _) => {
                        return Err(Error::format(
                            path,
                            format!(
                                "an @SQ line for {} has no usable LN, and reference ids are \
                                 counted over these lines, so every later one would shift",
                                name.unwrap_or("an unnamed sequence")
                            ),
                        ))
                    }
                }
            }
            "RG" => {
                if let Some(id) = field(line, "ID") {
                    read_groups.push(id.to_string());
                }
            }
            _ => {}
        }
    }
    let chr_map = ChrMap::from_entries(entries);
    Ok((header, chr_map, read_groups))
}

/// `(name, UR, M5)` for each `@SQ` line, which is what a reference is resolved
/// from.
/// Why this reference does not fit this file, or `None` if it does.
///
/// Two questions, and the second is the cheap half of the `M5` check §11 asks
/// for and this reader deliberately skips — hashing whole chromosomes at open
/// is not worth it, but comparing their lengths costs nothing and catches the
/// wrong assembly as well as the wrong names.
///
/// A partial match is not an error: a CRAM aligned against a full assembly may
/// name scaffolds a trimmed FASTA leaves out, and every read on a sequence the
/// FASTA does have is still perfectly readable. What is an error is *none* of
/// them matching, which means the wrong file, and a length disagreement, which
/// means the wrong version of the right file.
fn reference_complaint(
    source: &ReferenceSource,
    chr_map: &ChrMap,
    sequences: &[(String, Option<String>, Option<String>)],
) -> Option<String> {
    if sequences.is_empty() {
        return None;
    }
    let mut found = 0usize;
    for (name, ..) in sequences {
        let Some(length) = source.length(name) else {
            continue;
        };
        found += 1;
        let declared = chr_map.get(name).map(|entry| entry.size);
        if let Some(declared) = declared {
            if declared != length {
                return Some(format!(
                    "this reference has {name} at {length} bases where the file's header says                      {declared}, so it is not the assembly these reads were aligned to"
                ));
            }
        }
    }
    if found == 0 {
        let named = sequences
            .iter()
            .take(3)
            .map(|(name, ..)| name.as_str())
            .collect::<Vec<_>>()
            .join(", ");
        return Some(format!(
            "this reference holds none of the {} sequences this file names ({named}...), so              every sequence would read as N",
            sequences.len()
        ));
    }
    None
}

fn sequence_details(header: &SamHeader) -> Vec<(String, Option<String>, Option<String>)> {
    header
        .lines
        .iter()
        .filter(|line| line.kind == "SQ")
        .filter_map(|line| {
            Some((
                field(line, "SN")?.to_string(),
                field(line, "UR").map(str::to_string),
                field(line, "M5").map(str::to_string),
            ))
        })
        .collect()
}

fn field<'a>(line: &'a HeaderLine, tag: &str) -> Option<&'a str> {
    line.fields
        .iter()
        .find(|f| f.tag == tag)
        .map(|f| f.value.as_str())
}

/// What a walk holds: the loci, and where it is among them.
struct WalkPlan {
    loci: Vec<(usize, i64, i64)>,
    order: Vec<usize>,
    /// For a windowed walk, the start each window reports from — so an
    /// alignment reaching over a boundary belongs to one window only.
    min_starts: Vec<Option<i64>>,
    request: EntriesRequest,
    coverage: u64,
}

/// A walk over loci, one step per locus.
pub struct LocusWalk {
    plan: Arc<WalkPlan>,
    next: usize,
    tracker: Arc<ProgressTracker>,
}

impl std::fmt::Debug for LocusWalk {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("LocusWalk")
            .field("loci", &self.plan.loci.len())
            .field("next", &self.next)
            .finish()
    }
}

impl LocusWalk {
    pub fn plan(reader: &CramReader, req: &EntriesRequest) -> Result<Self> {
        Self::plan_with(reader, req, false)
    }

    fn plan_with(
        reader: &CramReader,
        req: &EntriesRequest,
        from_locus_start: bool,
    ) -> Result<Self> {
        let inner = reader.inner()?;
        // Planned up front, so a request that cannot be served says so before
        // the first step rather than partway through the walk.
        reader.index(inner)?;
        let resolved = reader.resolve(&req.locs)?;
        let mut order: Vec<usize> = (0..resolved.len()).collect();
        if req.sort_locations {
            order.sort_by_key(|i| resolved[*i]);
        }
        let loci: Vec<(usize, i64, i64)> = order.iter().map(|i| resolved[*i]).collect();
        let coverage = loci.iter().map(|(_, s, e)| (e - s).max(0) as u64).sum();
        let min_starts = if from_locus_start {
            loci.iter().map(|(_, start, _)| Some(*start)).collect()
        } else {
            vec![None; loci.len()]
        };
        Ok(Self {
            tracker: Arc::new(ProgressTracker::with_callback(
                coverage,
                req.progress.clone(),
            )),
            plan: Arc::new(WalkPlan {
                min_starts,
                loci,
                order,
                request: req.clone(),
                coverage,
            }),
            next: 0,
        })
    }

    /// The same walk, back at its first locus. Shares the plan.
    pub fn restarted(&self) -> Self {
        Self {
            plan: self.plan.clone(),
            next: 0,
            tracker: Arc::new(ProgressTracker::with_callback(
                self.plan.coverage,
                self.plan.request.progress.clone(),
            )),
        }
    }

    pub fn plan_windows(reader: &CramReader, req: &EntriesRequest, span: i64) -> Result<Self> {
        if span < 1 {
            return Err(Error::invalid(format!(
                "span must be positive (got {span})"
            )));
        }
        let locs = crate::bam::window_locs(&reader.chr_map, &req.locs.chr_ids, span)?;
        let windowed = EntriesRequest {
            locs,
            sort_locations: false,
            ..req.clone()
        };
        Self::plan_with(reader, &windowed, true)
    }

    pub fn len(&self) -> usize {
        self.plan.loci.len()
    }
    pub fn is_empty(&self) -> bool {
        self.plan.loci.is_empty()
    }
    pub fn order(&self) -> &[usize] {
        &self.plan.order
    }

    pub fn next_window(&mut self, reader: &CramReader) -> Option<Result<Vec<BamRecord>>> {
        if self.next >= self.plan.loci.len() {
            self.tracker.done_report();
            return None;
        }
        let index = self.next;
        let locus = self.plan.loci[index];
        let outcome = reader.inner().and_then(|inner| {
            let cram_index = reader.index(inner)?;
            reader.read_locus(inner, &cram_index, locus, &self.plan.request)
        });
        match outcome {
            // A failed step is not a step: the walk stays where it was, so the
            // call after it raises the same error rather than reading past a
            // locus nothing was read for.
            Err(e) => Some(Err(e)),
            Ok(mut records) => {
                self.next += 1;
                if let Some(min_start) = self.plan.min_starts[index] {
                    records.retain(|r| r.start() >= min_start);
                }
                let (_, start, end) = locus;
                self.tracker.add((end - start).max(0) as u64);
                Some(Ok(records))
            }
        }
    }
}

/// [`LocusWalk`] as a plain [`Iterator`], for Rust callers.
#[derive(Debug)]
pub struct LocusEntries<'a> {
    reader: &'a CramReader,
    walk: LocusWalk,
}

impl<'a> LocusEntries<'a> {
    fn plan(reader: &'a CramReader, req: &EntriesRequest) -> Result<Self> {
        Ok(Self {
            reader,
            walk: LocusWalk::plan(reader, req)?,
        })
    }
    pub fn len(&self) -> usize {
        self.walk.len()
    }
    pub fn is_empty(&self) -> bool {
        self.walk.is_empty()
    }
    pub fn order(&self) -> &[usize] {
        self.walk.order()
    }
}

impl Iterator for LocusEntries<'_> {
    type Item = Result<Vec<BamRecord>>;
    fn next(&mut self) -> Option<Self::Item> {
        self.walk.next_window(self.reader)
    }
}

/// [`LocusWalk`] over windows tiling whole references.
#[derive(Debug)]
pub struct WindowEntries<'a> {
    reader: &'a CramReader,
    walk: LocusWalk,
}

impl<'a> WindowEntries<'a> {
    fn plan(reader: &'a CramReader, req: &EntriesRequest, span: i64) -> Result<Self> {
        Ok(Self {
            reader,
            walk: LocusWalk::plan_windows(reader, req, span)?,
        })
    }
    pub fn len(&self) -> usize {
        self.walk.len()
    }
    pub fn is_empty(&self) -> bool {
        self.walk.is_empty()
    }
}

impl Iterator for WindowEntries<'_> {
    type Item = Result<Vec<BamRecord>>;
    fn next(&mut self) -> Option<Self::Item> {
        self.walk.next_window(self.reader)
    }
}