//! Virtual Dataset (VDS) assembly, following libhdf5's `H5Dvirtual.c`.
//!
//! A virtual dataset stores no data of its own: a list of mappings (kept in
//! the global heap) pairs a selection of the virtual dataspace with a
//! selection of a *source* dataset, in the same file (`"."`) or another one.
//! Reading it means reading each source and scattering the selected source
//! elements into the virtual buffer, pairing the two selections element by
//! element in row-major order. Elements no mapping supplies — unmapped
//! regions, and mappings whose source file or dataset does not exist — read
//! as the virtual dataset's **fill value**, as in libhdf5.
//!
//! Source files other than the virtual file itself are obtained through a
//! caller-supplied [`VdsFileResolver`], since this crate has no filesystem.

#[cfg(not(feature = "std"))]
use alloc::{format, string::String, vec, vec::Vec};

use crate::data_layout::{DataLayout, VdsMapping, parse_vds_mappings};
use crate::dataspace::Dataspace;
use crate::datatype::Datatype;
use crate::error::FormatError;
use crate::selection::{SerializedSelection, UNLIMITED};

/// Resolves the name of an external VDS source file, as stored in the
/// mapping, to that file's bytes.
///
/// `Ok(None)` means the file does not exist; its mappings then read as the
/// fill value, as libhdf5 does for a missing source. `Err` refuses the name
/// (e.g. a path the caller will not follow) and fails the read, so that a
/// refused source is never passed off as fill.
pub type VdsFileResolver<'a> = dyn Fn(&str) -> Result<Option<Vec<u8>>, FormatError> + 'a;

/// A fully assembled virtual dataset.
#[derive(Debug, Clone, PartialEq)]
pub struct VirtualData {
    /// The virtual dataset's extent.
    pub dims: Vec<u64>,
    /// Raw element bytes, row-major, in the virtual dataset's datatype.
    pub data: Vec<u8>,
    /// Number of elements no mapping supplied; they hold the fill value.
    pub unmapped: u64,
}

fn vds_err(msg: impl Into<String>) -> FormatError {
    FormatError::ChunkedReadError(msg.into())
}

/// Upper bound on the printf-style source datasets probed for one mapping.
const MAX_PRINTF_BLOCKS: u64 = 1 << 20;

/// A source file or dataset name, parsed for printf-style `%b` block-number
/// substitutions (`H5D_virtual_parse_source_name`): `%b` is the block
/// number, `%%` a literal `%`, and any other `%` sequence is invalid.
#[derive(Debug, Clone, PartialEq)]
struct SourceName {
    /// Literal text around the substitutions: `segments.len() == subs + 1`.
    segments: Vec<String>,
}

impl SourceName {
    fn parse(name: &str) -> Result<SourceName, FormatError> {
        let mut segments = vec![String::new()];
        let mut chars = name.chars();
        while let Some(c) = chars.next() {
            if c != '%' {
                segments.last_mut().expect("never empty").push(c);
                continue;
            }
            match chars.next() {
                Some('b') => segments.push(String::new()),
                Some('%') => segments.last_mut().expect("never empty").push('%'),
                _ => {
                    return Err(vds_err(format!(
                        "invalid format specifier in VDS source name {name:?}"
                    )));
                }
            }
        }
        Ok(SourceName { segments })
    }

    /// Number of `%b` substitutions.
    fn subs(&self) -> usize {
        self.segments.len() - 1
    }

    /// The name with every `%b` replaced by `block`.
    fn build(&self, block: u64) -> String {
        let mut out = String::new();
        for (i, seg) in self.segments.iter().enumerate() {
            if i > 0 {
                out.push_str(&format!("{block}"));
            }
            out.push_str(seg);
        }
        out
    }
}

/// How a mapping's selections relate, as libhdf5 classifies them.
#[derive(Debug, Clone, Copy, PartialEq)]
enum Kind {
    /// Both selections have a fixed size.
    Fixed,
    /// Both are unlimited (in `vdim` / `sdim`): the mapping grows with the
    /// source dataset's extent.
    Unlimited { vdim: usize, sdim: usize },
    /// The virtual selection repeats a block without limit in `vdim`; block
    /// `j` comes from the source named by substituting `j` for `%b`.
    Printf { vdim: usize },
}

/// One mapping with its selections decoded.
struct Mapping {
    file: SourceName,
    dataset: SourceName,
    vsel: SerializedSelection,
    ssel: SerializedSelection,
    kind: Kind,
}

impl Mapping {
    fn new(m: VdsMapping) -> Result<Mapping, FormatError> {
        let (vsel, _) = SerializedSelection::decode(&m.virtual_selection)?;
        let (ssel, _) = SerializedSelection::decode(&m.source_selection)?;
        let file = SourceName::parse(&m.source_file)?;
        let dataset = SourceName::parse(&m.source_dataset)?;
        let subs = file.subs() + dataset.subs();
        // The checks of H5D_virtual_check_mapping_pre/_post.
        let kind = match (vsel.unlimited_dim(), ssel.unlimited_dim()) {
            (Some(vdim), None) => {
                if subs == 0 {
                    return Err(vds_err(
                        "unlimited virtual selection with a limited source selection \
                         and no %b in the source names",
                    ));
                }
                match &vsel {
                    SerializedSelection::Regular { count, block, .. }
                        if count[vdim] == UNLIMITED && block[vdim] != UNLIMITED => {}
                    _ => {
                        return Err(vds_err(
                            "printf VDS mapping needs a virtual selection with an unlimited count",
                        ));
                    }
                }
                Kind::Printf { vdim }
            }
            (Some(vdim), Some(sdim)) => {
                if non_unlimited_elements(&vsel, vdim) != non_unlimited_elements(&ssel, sdim) {
                    return Err(vds_err(
                        "unlimited VDS mapping: virtual and source selections differ \
                         outside the unlimited dimension",
                    ));
                }
                Kind::Unlimited { vdim, sdim }
            }
            (None, Some(_)) => {
                return Err(vds_err(
                    "VDS mapping with an unlimited source selection and a limited \
                     virtual selection is not supported",
                ));
            }
            (None, None) => Kind::Fixed,
        };
        if subs > 0 && !matches!(kind, Kind::Printf { .. }) {
            return Err(vds_err(
                "%b in a VDS source name without an unlimited virtual selection",
            ));
        }
        Ok(Mapping {
            file,
            dataset,
            vsel,
            ssel,
            kind,
        })
    }
}

/// Elements a regular selection selects outside dimension `skip`.
fn non_unlimited_elements(sel: &SerializedSelection, skip: usize) -> Option<u64> {
    let SerializedSelection::Regular { count, block, .. } = sel else {
        return None;
    };
    (0..count.len())
        .filter(|&d| d != skip)
        .try_fold(1u64, |acc, d| {
            acc.checked_mul(count[d].checked_mul(block[d])?)
        })
}

/// Load and decode the mapping list of a virtual layout.
fn load_mappings(
    file_data: &[u8],
    layout: &DataLayout,
    length_size: u8,
) -> Result<Vec<Mapping>, FormatError> {
    let DataLayout::Virtual {
        global_heap_address,
        global_heap_index,
        ..
    } = layout
    else {
        return Err(vds_err("not a virtual dataset layout"));
    };
    let Some(addr) = *global_heap_address else {
        return Ok(Vec::new());
    };
    let coll =
        crate::global_heap::GlobalHeapCollection::parse(file_data, addr as usize, length_size)?;
    let index = u16::try_from(*global_heap_index)
        .map_err(|_| vds_err("VDS mapping heap index out of range"))?;
    let obj = coll
        .get_object(index)
        .ok_or(FormatError::GlobalHeapObjectNotFound {
            collection_address: addr,
            index,
        })?;
    parse_vds_mappings(&obj.data, length_size)?
        .into_iter()
        .map(Mapping::new)
        .collect()
}

/// `H5S__hyper_get_clip_diminfo`: the count and block a regular selection has
/// in its unlimited dimension once clipped to `clip`.
fn clip_diminfo(start: u64, stride: u64, count: u64, block: u64, clip: u64) -> (u64, u64) {
    if start >= clip {
        if block == UNLIMITED {
            (count, 0)
        } else {
            (0, block)
        }
    } else if block == UNLIMITED || block == stride {
        (1, clip - start)
    } else {
        ((clip - start).div_ceil(stride.max(1)), block)
    }
}

/// The unlimited-dimension parameters (start, stride, count, block) of a
/// regular selection.
fn unlim_diminfo(sel: &SerializedSelection, dim: usize) -> Result<[u64; 4], FormatError> {
    match sel {
        SerializedSelection::Regular {
            start,
            stride,
            count,
            block,
        } => Ok([start[dim], stride[dim], count[dim], block[dim]]),
        _ => Err(vds_err(
            "unlimited VDS selection is not a regular hyperslab",
        )),
    }
}

/// `H5S_hyper_get_clip_extent_match` with `incl_trail = false` (the
/// "last available" view): the extent to clip `clip_sel` (unlimited in
/// `clip_dim`) to so that it holds as many slices as `match_sel` (unlimited
/// in `match_dim`) holds when clipped to `match_clip`.
fn clip_extent_match(
    clip_sel: &SerializedSelection,
    clip_dim: usize,
    match_sel: &SerializedSelection,
    match_dim: usize,
    match_clip: u64,
) -> Result<u64, FormatError> {
    let overflow = || FormatError::Overflow("VDS clip extent overflow".into());
    let [mstart, mstride, mcount, mblock] = unlim_diminfo(match_sel, match_dim)?;
    let (count, block) = clip_diminfo(mstart, mstride, mcount, mblock, match_clip);
    let slices = if block == 0 || count == 0 {
        0
    } else if count == 1 {
        block
    } else {
        let mut n = block.checked_mul(count).ok_or_else(overflow)?;
        let span = mstride
            .checked_mul(count - 1)
            .and_then(|s| s.checked_add(block))
            .ok_or_else(overflow)?;
        let room = match_clip - mstart;
        if span > room {
            n -= span - room;
        }
        n
    };

    // H5S__hyper_get_clip_extent_real
    let [start, stride, _, block] = unlim_diminfo(clip_sel, clip_dim)?;
    if slices == 0 {
        return Ok(0);
    }
    let extent = if block == UNLIMITED || block == stride {
        start.checked_add(slices)
    } else {
        let full = slices / block;
        let rem = slices - full * block;
        if rem > 0 {
            full.checked_mul(stride)
                .and_then(|o| start.checked_add(o))
                .and_then(|e| e.checked_add(rem))
        } else {
            (full - 1)
                .checked_mul(stride)
                .and_then(|o| start.checked_add(o))
                .and_then(|e| e.checked_add(block))
        }
    };
    extent.ok_or_else(overflow)
}

/// How each mapping is read, and the resulting extent.
struct Plan {
    dims: Vec<u64>,
    steps: Vec<Step>,
}

#[derive(Clone, Copy)]
enum Step {
    Fixed,
    /// Read with the virtual selection clipped to `vclip` and the source
    /// selection to the source's extent; `None` if the source is missing.
    Unlimited(Option<u64>),
    /// Read blocks `0..blocks`.
    Printf(u64),
}

/// Work out the extent libhdf5 gives the virtual dataset
/// (`H5D__virtual_set_extent_unlim`, default view `H5D_VDS_LAST_AVAILABLE`
/// with a printf gap of 0) and how much of each unlimited mapping is read.
fn plan(mappings: &[Mapping], stored: &[u64], sources: &mut Sources) -> Result<Plan, FormatError> {
    let overflow = || FormatError::Overflow("VDS extent overflow".into());
    let rank = stored.len();
    let mut new_dims: Vec<Option<u64>> = vec![None; rank];
    // Minimum extent needed by the limited parts of every virtual selection
    // (H5D_virtual_update_min_dims).
    let mut min_dims = vec![0u64; rank];
    let mut steps = Vec::with_capacity(mappings.len());

    for m in mappings {
        let skip = match m.kind {
            Kind::Unlimited { vdim, .. } | Kind::Printf { vdim } => Some(vdim),
            Kind::Fixed => None,
        };
        if let Some(ends) = selection_bounds_end(&m.vsel)? {
            if ends.len() != rank {
                return Err(vds_err("VDS selection rank does not match dataspace rank"));
            }
            for d in (0..rank).filter(|&d| Some(d) != skip) {
                min_dims[d] = min_dims[d].max(ends[d].checked_add(1).ok_or_else(overflow)?);
            }
        }

        let (vdim, clip, step) = match m.kind {
            Kind::Fixed => {
                steps.push(Step::Fixed);
                continue;
            }
            Kind::Unlimited { vdim, sdim } => {
                match sources.dims(&m.file.build(0), &m.dataset.build(0))? {
                    Some(src_dims) => {
                        let extent = *src_dims.get(sdim).ok_or_else(|| {
                            vds_err("VDS source rank does not match its selection")
                        })?;
                        let clip = clip_extent_match(&m.vsel, vdim, &m.ssel, sdim, extent)?;
                        (vdim, clip, Step::Unlimited(Some(clip)))
                    }
                    None => (vdim, 0, Step::Unlimited(None)),
                }
            }
            Kind::Printf { vdim } => {
                // With a gap of 0 the search stops at the first missing
                // source dataset.
                let mut found = 0u64;
                while sources
                    .dims(&m.file.build(found), &m.dataset.build(found))?
                    .is_some()
                {
                    found += 1;
                    if found > MAX_PRINTF_BLOCKS {
                        return Err(vds_err("too many printf-style VDS source datasets"));
                    }
                }
                let clip = if found == 0 {
                    0
                } else {
                    // End of block `found - 1` in the unlimited dimension.
                    let [start, stride, _, block] = unlim_diminfo(&m.vsel, vdim)?;
                    (found - 1)
                        .checked_mul(stride)
                        .and_then(|o| start.checked_add(o))
                        .and_then(|e| e.checked_add(block))
                        .ok_or_else(overflow)?
                };
                (vdim, clip, Step::Printf(found))
            }
        };
        if vdim >= rank {
            return Err(vds_err("VDS selection rank does not match dataspace rank"));
        }
        new_dims[vdim] = Some(new_dims[vdim].map_or(clip, |n| n.max(clip)));
        steps.push(step);
    }

    let dims = (0..rank)
        .map(|d| match new_dims[d] {
            None => stored[d],
            Some(n) => n.max(min_dims[d]),
        })
        .collect();
    Ok(Plan { dims, steps })
}

/// The last selected coordinate in each dimension (`H5S_SELECT_BOUNDS`),
/// ignoring any unlimited dimension; `None` for ALL/NONE and empty selections.
fn selection_bounds_end(sel: &SerializedSelection) -> Result<Option<Vec<u64>>, FormatError> {
    let overflow = || FormatError::Overflow("VDS selection bounds overflow".into());
    match sel {
        SerializedSelection::All | SerializedSelection::None => Ok(None),
        SerializedSelection::Regular {
            start,
            stride,
            count,
            block,
        } => {
            if count.contains(&0) || block.contains(&0) {
                return Ok(None);
            }
            let mut ends = Vec::with_capacity(start.len());
            for d in 0..start.len() {
                if count[d] == UNLIMITED || block[d] == UNLIMITED {
                    ends.push(0);
                    continue;
                }
                let end = (count[d] - 1)
                    .checked_mul(stride[d])
                    .and_then(|o| start[d].checked_add(o))
                    .and_then(|e| e.checked_add(block[d] - 1))
                    .ok_or_else(overflow)?;
                ends.push(end);
            }
            Ok(Some(ends))
        }
        SerializedSelection::Blocks { rank, ends, .. } => {
            if ends.is_empty() {
                return Ok(None);
            }
            let mut max = vec![0u64; *rank];
            for e in ends.chunks_exact(*rank) {
                for d in 0..*rank {
                    max[d] = max[d].max(e[d]);
                }
            }
            Ok(Some(max))
        }
    }
}

/// The virtual dataset's extent as libhdf5 reports it (`H5Dget_space`).
///
/// For a virtual dataset whose mappings are all of fixed size this is the
/// stored dataspace. With unlimited or printf-style mappings libhdf5
/// recomputes the unlimited dimension from the sources present (the default
/// "last available" view: the largest extent any mapping can fill), which
/// needs the source files, read through `resolver`.
pub fn virtual_dataset_extent(
    file_data: &[u8],
    layout: &DataLayout,
    dataspace: &Dataspace,
    _offset_size: u8,
    length_size: u8,
    resolver: Option<&VdsFileResolver>,
) -> Result<Vec<u64>, FormatError> {
    let mappings = load_mappings(file_data, layout, length_size)?;
    if mappings.iter().all(|m| m.kind == Kind::Fixed) {
        return Ok(dataspace.dimensions.clone());
    }
    let mut sources = Sources::new(file_data, resolver);
    Ok(plan(&mappings, &dataspace.dimensions, &mut sources)?.dims)
}

/// Read a whole virtual dataset, at the extent
/// [`virtual_dataset_extent`] reports.
///
/// `fill` is the virtual dataset's fill value (from its fill value message;
/// `None` for the default of zeros); every element no mapping supplies holds
/// it. External source files are read through `resolver`; without one, a
/// mapping to another file is an error.
#[allow(clippy::too_many_arguments)]
pub fn read_virtual_dataset(
    file_data: &[u8],
    layout: &DataLayout,
    dataspace: &Dataspace,
    datatype: &Datatype,
    fill: Option<&[u8]>,
    _offset_size: u8,
    length_size: u8,
    resolver: Option<&VdsFileResolver>,
) -> Result<VirtualData, FormatError> {
    let mappings = load_mappings(file_data, layout, length_size)?;
    let mut sources = Sources::new(file_data, resolver);
    let Plan { dims, steps } = plan(&mappings, &dataspace.dimensions, &mut sources)?;

    let elem_size = datatype.type_size() as usize;
    let total = dims
        .iter()
        .try_fold(1u64, |acc, &d| acc.checked_mul(d))
        .ok_or_else(|| FormatError::Overflow("virtual dataset extent".into()))?;
    let mut data = crate::chunked_read::alloc_output(crate::chunked_read::checked_byte_len(
        total, elem_size,
    )?)?;
    if let Some(fill) = fill.filter(|f| f.len() == elem_size && f.iter().any(|&b| b != 0)) {
        for element in data.chunks_exact_mut(elem_size) {
            element.copy_from_slice(fill);
        }
    }
    let mut mapped = vec![false; usize::try_from(total).map_err(|_| vds_err("VDS too large"))?];

    for (m, step) in mappings.iter().zip(&steps) {
        match (*step, m.kind) {
            (Step::Fixed, _) => {
                let (file, dset) = (m.file.build(0), m.dataset.build(0));
                let Some(src) = sources.dataset(&file, &dset, datatype)? else {
                    continue; // missing source file or dataset: fill
                };
                let vidx = selection_indices(&m.vsel, &dims, None)?;
                let sidx = selection_indices(&m.ssel, &src.dims, None)?;
                scatter(&mut data, &mut mapped, &src.raw, &vidx, &sidx, elem_size)?;
            }
            (Step::Unlimited(Some(vclip)), Kind::Unlimited { vdim, sdim }) => {
                let (file, dset) = (m.file.build(0), m.dataset.build(0));
                let Some(src) = sources.dataset(&file, &dset, datatype)? else {
                    continue;
                };
                let vidx = selection_indices(&m.vsel, &dims, Some((vdim, vclip)))?;
                let extent = *src
                    .dims
                    .get(sdim)
                    .ok_or_else(|| vds_err("VDS source rank does not match its selection"))?;
                let sidx = selection_indices(&m.ssel, &src.dims, Some((sdim, extent)))?;
                scatter(&mut data, &mut mapped, &src.raw, &vidx, &sidx, elem_size)?;
            }
            (Step::Unlimited(None), _) => {}
            (Step::Printf(blocks), Kind::Printf { vdim }) => {
                for j in 0..blocks {
                    let Some(src) =
                        sources.dataset(&m.file.build(j), &m.dataset.build(j), datatype)?
                    else {
                        continue;
                    };
                    let vblock = unlim_block(&m.vsel, vdim, j)?;
                    let vidx = selection_indices(&vblock, &dims, None)?;
                    let sidx = selection_indices(&m.ssel, &src.dims, None)?;
                    scatter(&mut data, &mut mapped, &src.raw, &vidx, &sidx, elem_size)?;
                }
            }
            _ => return Err(vds_err("internal error: VDS plan does not match mapping")),
        }
    }

    let unmapped = mapped.iter().filter(|&&m| !m).count() as u64;
    Ok(VirtualData {
        dims,
        data,
        unmapped,
    })
}

/// `H5S_hyper_get_unlim_block`: block `j` of a selection whose count is
/// unlimited in `dim`.
fn unlim_block(
    sel: &SerializedSelection,
    dim: usize,
    j: u64,
) -> Result<SerializedSelection, FormatError> {
    let SerializedSelection::Regular {
        start,
        stride,
        count,
        block,
    } = sel
    else {
        return Err(vds_err("printf VDS selection is not a regular hyperslab"));
    };
    let mut start = start.clone();
    let mut count = count.clone();
    start[dim] = j
        .checked_mul(stride[dim])
        .and_then(|o| start[dim].checked_add(o))
        .ok_or_else(|| FormatError::Overflow("VDS block start overflow".into()))?;
    count[dim] = 1;
    Ok(SerializedSelection::Regular {
        start,
        stride: stride.clone(),
        count,
        block: block.clone(),
    })
}

/// Copy source element `sidx[i]` to virtual element `vidx[i]` for every `i`.
fn scatter(
    out: &mut [u8],
    mapped: &mut [bool],
    src: &[u8],
    vidx: &[u64],
    sidx: &[u64],
    elem_size: usize,
) -> Result<(), FormatError> {
    if vidx.len() != sidx.len() {
        return Err(vds_err("virtual/source selection element counts differ"));
    }
    for (&v, &s) in vidx.iter().zip(sidx) {
        let (vo, so) = (v as usize * elem_size, s as usize * elem_size);
        if vo + elem_size > out.len() || so + elem_size > src.len() {
            return Err(vds_err("virtual dataset selection out of bounds"));
        }
        out[vo..vo + elem_size].copy_from_slice(&src[so..so + elem_size]);
        mapped[v as usize] = true;
    }
    Ok(())
}

/// Row-major linear indices of the elements `sel` selects in a dataspace of
/// shape `dims`, in the order libhdf5 iterates them (row-major).
///
/// `clip` = `(dim, limit)` drops every coordinate `>= limit` in `dim` — the
/// clipping libhdf5 applies to an unlimited selection
/// (`H5S_hyper_clip_unlim`); an unlimited selection must be clipped.
fn selection_indices(
    sel: &SerializedSelection,
    dims: &[u64],
    clip: Option<(usize, u64)>,
) -> Result<Vec<u64>, FormatError> {
    let overflow = || FormatError::Overflow("VDS selection index overflow".into());
    let rank = dims.len();
    let total = dims
        .iter()
        .try_fold(1u64, |acc, &d| acc.checked_mul(d))
        .ok_or_else(overflow)?;
    let mut row_stride = vec![1u64; rank];
    for d in (0..rank.saturating_sub(1)).rev() {
        row_stride[d] = row_stride[d + 1]
            .checked_mul(dims[d + 1])
            .ok_or_else(overflow)?;
    }
    if sel.rank().is_some_and(|r| r != rank) {
        return Err(vds_err("VDS selection rank does not match dataspace rank"));
    }
    let limit = |d: usize| match clip {
        Some((cd, l)) if cd == d => l,
        _ => u64::MAX,
    };

    match sel {
        SerializedSelection::All => Ok((0..total).collect()),
        SerializedSelection::None => Ok(Vec::new()),
        SerializedSelection::Regular {
            start,
            stride,
            count,
            block,
        } => {
            // Selected coordinates along each dimension, in order.
            let mut per_dim: Vec<Vec<u64>> = Vec::with_capacity(rank);
            for d in 0..rank {
                let lim = limit(d);
                if (count[d] == UNLIMITED || block[d] == UNLIMITED) && lim == u64::MAX {
                    return Err(vds_err("unlimited VDS selection was not clipped"));
                }
                let mut coords = Vec::new();
                let mut ci = 0u64;
                'blocks: while ci < count[d] {
                    let base = ci
                        .checked_mul(stride[d])
                        .and_then(|o| start[d].checked_add(o))
                        .ok_or_else(overflow)?;
                    if base >= lim {
                        break;
                    }
                    let mut bi = 0u64;
                    while bi < block[d] {
                        let coord = base.checked_add(bi).ok_or_else(overflow)?;
                        if coord >= lim {
                            break 'blocks;
                        }
                        // Past the extent is malformed; bail before the list
                        // can grow without bound.
                        if coord >= dims[d] {
                            return Err(vds_err("VDS selection exceeds the dataspace extent"));
                        }
                        coords.push(coord);
                        bi += 1;
                    }
                    ci += 1;
                }
                per_dim.push(coords);
            }
            if per_dim.iter().any(|c| c.is_empty()) {
                return Ok(Vec::new());
            }
            let n = per_dim
                .iter()
                .try_fold(1usize, |acc, c| acc.checked_mul(c.len()))
                .ok_or_else(overflow)?;
            let mut out = Vec::with_capacity(n);
            let mut idx = vec![0usize; rank];
            loop {
                let lin: u64 = (0..rank).map(|d| per_dim[d][idx[d]] * row_stride[d]).sum();
                out.push(lin);
                // Mixed-radix increment, last dimension fastest.
                let mut d = rank;
                loop {
                    if d == 0 {
                        return Ok(out);
                    }
                    d -= 1;
                    idx[d] += 1;
                    if idx[d] < per_dim[d].len() {
                        break;
                    }
                    idx[d] = 0;
                }
            }
        }
        SerializedSelection::Blocks {
            rank: _,
            starts,
            ends,
        } => {
            // libhdf5 serializes the union as disjoint blocks, so their volumes
            // never add up to more than the dataspace.
            let mut volume = 0u64;
            for (s, e) in starts.chunks_exact(rank).zip(ends.chunks_exact(rank)) {
                for d in 0..rank {
                    if e[d] >= dims[d] {
                        return Err(vds_err("VDS selection exceeds the dataspace extent"));
                    }
                }
                let v = s
                    .iter()
                    .zip(e)
                    .try_fold(1u64, |acc, (&s, &e)| acc.checked_mul(e - s + 1))
                    .ok_or_else(overflow)?;
                volume = volume.checked_add(v).ok_or_else(overflow)?;
                if volume > total {
                    return Err(vds_err("VDS selection blocks overlap"));
                }
            }
            let mut out = Vec::with_capacity(volume as usize);
            for (s, e) in starts.chunks_exact(rank).zip(ends.chunks_exact(rank)) {
                let mut cur = s.to_vec();
                'block: loop {
                    if (0..rank).all(|d| cur[d] < limit(d)) {
                        out.push((0..rank).map(|d| cur[d] * row_stride[d]).sum());
                    }
                    for d in (0..rank).rev() {
                        if cur[d] < e[d] {
                            cur[d] += 1;
                            continue 'block;
                        }
                        cur[d] = s[d];
                    }
                    break;
                }
            }
            out.sort_unstable();
            out.dedup();
            Ok(out)
        }
    }
}

/// A source dataset's decoded contents.
struct SourceData {
    dims: Vec<u64>,
    raw: Vec<u8>,
}

/// Source files and datasets, fetched on demand. The most recently used
/// external file is kept, since consecutive mappings usually share one.
struct Sources<'a, 'r> {
    file_data: &'a [u8],
    resolver: Option<&'r VdsFileResolver<'r>>,
    cached_file: Option<(String, Option<Vec<u8>>)>,
}

impl<'a, 'r> Sources<'a, 'r> {
    fn new(file_data: &'a [u8], resolver: Option<&'r VdsFileResolver<'r>>) -> Self {
        Sources {
            file_data,
            resolver,
            cached_file: None,
        }
    }

    /// The bytes of source file `name`, or `None` if it does not exist.
    fn file(&mut self, name: &str) -> Result<Option<&[u8]>, FormatError> {
        if name == "." {
            return Ok(Some(self.file_data));
        }
        if self.cached_file.as_ref().is_none_or(|(n, _)| n != name) {
            let resolver = self.resolver.ok_or_else(|| {
                vds_err("external-file virtual dataset sources require a file resolver")
            })?;
            self.cached_file = Some((String::from(name), resolver(name)?));
        }
        // An external file is handed over whole; its addresses are relative
        // to its superblock, so skip any user block.
        match self.cached_file.as_ref().and_then(|(_, b)| b.as_deref()) {
            Some(bytes) => Ok(Some(crate::signature::split_user_block(bytes)?.1)),
            None => Ok(None),
        }
    }

    /// The extent of source dataset `path` in file `file`, or `None` when
    /// either does not exist.
    fn dims(&mut self, file: &str, path: &str) -> Result<Option<Vec<u64>>, FormatError> {
        let Some(bytes) = self.file(file)? else {
            return Ok(None);
        };
        Ok(open_source(bytes, path)?.map(|s| s.dataspace.dimensions))
    }

    /// Read source dataset `path` from file `file`, or `None` when either
    /// does not exist. Its datatype must be the virtual dataset's: libhdf5
    /// converts between types here, which is not supported.
    fn dataset(
        &mut self,
        file: &str,
        path: &str,
        datatype: &Datatype,
    ) -> Result<Option<SourceData>, FormatError> {
        // Variable-length and reference elements are addresses into the file
        // that holds them (global-heap IDs, object addresses). Copied out of
        // another file they would be decoded against the virtual dataset's
        // file and name some other object, so refuse rather than return them.
        if file != "." && holds_file_addresses(datatype) {
            return Err(vds_err(format!(
                "VDS source {path} in {file}: variable-length and reference data \
                 from another file is not supported"
            )));
        }
        let Some(bytes) = self.file(file)? else {
            return Ok(None);
        };
        let Some(src) = open_source(bytes, path)? else {
            return Ok(None);
        };
        read_source(bytes, src, path, datatype).map(Some)
    }
}

/// Whether elements of `dt` contain addresses into their own file:
/// variable-length data (global-heap IDs) or references.
fn holds_file_addresses(dt: &Datatype) -> bool {
    match dt {
        Datatype::VariableLength { .. } | Datatype::Reference { .. } => true,
        Datatype::Compound { members, .. } => {
            members.iter().any(|m| holds_file_addresses(&m.datatype))
        }
        Datatype::Array { base_type, .. } | Datatype::Enumeration { base_type, .. } => {
            holds_file_addresses(base_type)
        }
        _ => false,
    }
}

/// An opened source dataset's object header.
struct OpenSource {
    offset_size: u8,
    length_size: u8,
    header: crate::object_header::ObjectHeader,
    dataspace: Dataspace,
}

fn source_message<'h>(
    src: &'h OpenSource,
    path: &str,
    t: crate::message_type::MessageType,
) -> Result<&'h crate::object_header::HeaderMessage, FormatError> {
    src.header
        .messages
        .iter()
        .find(|m| m.msg_type == t)
        .ok_or_else(|| vds_err(format!("VDS source {path} has no {t:?} message")))
}

/// Open source dataset `path` of the file in `file_data`, or `None` if there
/// is no such object (libhdf5 reads a missing source as fill).
fn open_source(file_data: &[u8], path: &str) -> Result<Option<OpenSource>, FormatError> {
    use crate::message_type::MessageType;
    use crate::shared_message::message_data_with_sohm;

    // `file_data` starts at the superblock (see `Sources::file`).
    let sb = crate::superblock::Superblock::parse(file_data, 0)?;
    let (os, ls) = (sb.offset_size, sb.length_size);
    let addr = match crate::group_v2::resolve_path_any(file_data, &sb, path) {
        Ok(a) => a,
        Err(FormatError::PathNotFound(_)) => return Ok(None),
        Err(e) => return Err(e),
    };
    let header = crate::object_header::ObjectHeader::parse(file_data, addr as usize, os, ls)?;
    let mut src = OpenSource {
        offset_size: os,
        length_size: ls,
        header,
        dataspace: Dataspace {
            space_type: crate::dataspace::DataspaceType::Null,
            rank: 0,
            dimensions: Vec::new(),
            max_dimensions: None,
        },
    };
    let ds_msg = source_message(&src, path, MessageType::Dataspace)?;
    src.dataspace = Dataspace::parse(&message_data_with_sohm(file_data, ds_msg, os, ls)?, ls)?;
    Ok(Some(src))
}

/// Read an opened source dataset in full (its own fill value applied to
/// unallocated chunks).
fn read_source(
    file_data: &[u8],
    src: OpenSource,
    path: &str,
    datatype: &Datatype,
) -> Result<SourceData, FormatError> {
    use crate::filter_pipeline::FilterPipeline;
    use crate::message_type::MessageType;
    use crate::shared_message::message_data_with_sohm;

    let (os, ls) = (src.offset_size, src.length_size);
    let dt_msg = source_message(&src, path, MessageType::Datatype)?;
    let (src_type, _) = Datatype::parse(&message_data_with_sohm(file_data, dt_msg, os, ls)?)?;
    // libhdf5 converts each source to the virtual dataset's type. Only the
    // conversion that is a pure byte swap is done here.
    let swap = if &src_type == datatype {
        false
    } else if differs_only_in_byte_order(&src_type, datatype) {
        true
    } else {
        return Err(vds_err(format!(
            "VDS source {path} has a different datatype from the virtual dataset \
             (only a byte-order conversion is supported)"
        )));
    };
    let layout = DataLayout::parse(
        &source_message(&src, path, MessageType::DataLayout)?.data,
        os,
        ls,
    )?;
    // A source that is itself virtual could form a cycle (A -> B -> A) and
    // recurse without bound. Nested virtual sources are not supported.
    if matches!(layout, DataLayout::Virtual { .. }) {
        return Err(vds_err(
            "virtual dataset source is itself virtual (unsupported)",
        ));
    }
    let pipeline = src
        .header
        .messages
        .iter()
        .find(|m| m.msg_type == MessageType::FilterPipeline)
        .map(|m| {
            message_data_with_sohm(file_data, m, os, ls).and_then(|d| FilterPipeline::parse(&d))
        })
        .transpose()?;
    let raw = crate::fill_value::read_full_with_fill(
        &src.header.messages,
        file_data,
        &layout,
        &src.dataspace,
        src_type.type_size() as usize,
        os,
        ls,
        || {
            crate::data_read::read_raw_data_full(
                file_data,
                &layout,
                &src.dataspace,
                &src_type,
                pipeline.as_ref(),
                os,
                ls,
            )
        },
    )?;
    let mut raw = raw;
    if swap {
        let size = datatype.type_size() as usize;
        for element in raw.chunks_exact_mut(size) {
            element.reverse();
        }
    }
    Ok(SourceData {
        dims: src.dataspace.dimensions,
        raw,
    })
}

/// Whether two numeric types are identical apart from their byte order (both
/// little- or big-endian), so converting one to the other is a byte swap.
fn differs_only_in_byte_order(a: &Datatype, b: &Datatype) -> bool {
    use crate::datatype::DatatypeByteOrder::{BigEndian, LittleEndian};
    let mut a = a.clone();
    match &mut a {
        Datatype::FixedPoint { byte_order, .. }
        | Datatype::FloatingPoint { byte_order, .. }
        | Datatype::BitField { byte_order, .. } => {
            *byte_order = match byte_order {
                LittleEndian => BigEndian,
                BigEndian => LittleEndian,
                _ => return false,
            }
        }
        _ => return false,
    }
    &a == b
}

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

    fn regular(start: u64, stride: u64, count: u64, block: u64) -> SerializedSelection {
        SerializedSelection::Regular {
            start: vec![start],
            stride: vec![stride],
            count: vec![count],
            block: vec![block],
        }
    }

    #[test]
    fn printf_names_follow_libhdf5() {
        let n = SourceName::parse("f-%b.h5").unwrap();
        assert_eq!(n.subs(), 1);
        assert_eq!(n.build(12), "f-12.h5");
        let n = SourceName::parse("100%%_%b_%b").unwrap();
        assert_eq!(n.build(3), "100%_3_3");
        let n = SourceName::parse("plain%%name").unwrap();
        assert_eq!((n.subs(), n.build(7)), (0, "plain%name".to_string()));
        // Anything else after '%' (or a trailing '%') is invalid.
        assert!(SourceName::parse("a%d").is_err());
        assert!(SourceName::parse("a%").is_err());
    }

    #[test]
    fn clip_extent_matches_libhdf5_arithmetic() {
        // 7 source slices (contiguous unlimited source) into blocks of 3
        // every 4: two full blocks and one slice of a third -> extent 9.
        let src = regular(0, 1, UNLIMITED, 1);
        let v = regular(0, 4, UNLIMITED, 3);
        assert_eq!(clip_extent_match(&v, 0, &src, 0, 7).unwrap(), 9);
        // Exactly two blocks: the extent ends at the end of the last block.
        assert_eq!(clip_extent_match(&v, 0, &src, 0, 6).unwrap(), 7);
        // An empty source gives an empty mapping.
        assert_eq!(clip_extent_match(&v, 0, &src, 0, 0).unwrap(), 0);
        // Unlimited block: the extent is start + slices.
        let vb = regular(2, 1, 1, UNLIMITED);
        assert_eq!(clip_extent_match(&vb, 0, &src, 0, 5).unwrap(), 7);
        // A strided source clipped mid-block counts only the selected slices.
        let src2 = regular(1, 4, UNLIMITED, 2); // 1,2, 5,6, 9,10 ...
        let dense = regular(0, 1, UNLIMITED, 1);
        assert_eq!(clip_extent_match(&dense, 0, &src2, 0, 6).unwrap(), 3);
    }

    #[test]
    fn clipped_selection_drops_the_partial_tail() {
        let v = regular(0, 4, UNLIMITED, 3);
        assert_eq!(
            selection_indices(&v, &[20], Some((0, 9))).unwrap(),
            vec![0, 1, 2, 4, 5, 6, 8]
        );
        // Unclipped unlimited selections cannot be enumerated.
        assert!(selection_indices(&v, &[20], None).is_err());
    }

    #[test]
    fn unlimited_mapping_rules() {
        let sel = |s: &SerializedSelection| -> Vec<u8> {
            // Serialize as version 2 (8-byte regular).
            let SerializedSelection::Regular {
                start,
                stride,
                count,
                block,
            } = s
            else {
                unreachable!()
            };
            let mut b = Vec::new();
            b.extend_from_slice(&2u32.to_le_bytes());
            b.extend_from_slice(&2u32.to_le_bytes());
            b.push(1);
            b.extend_from_slice(&0u32.to_le_bytes());
            b.extend_from_slice(&(start.len() as u32).to_le_bytes());
            for d in 0..start.len() {
                for v in [start[d], stride[d], count[d], block[d]] {
                    b.extend_from_slice(&v.to_le_bytes());
                }
            }
            b
        };
        let mapping = |file: &str, v: &SerializedSelection, s: &SerializedSelection| VdsMapping {
            source_file: file.into(),
            source_dataset: "d".into(),
            source_selection: sel(s),
            virtual_selection: sel(v),
        };
        let unlim = regular(0, 10, UNLIMITED, 10);
        let fixed = regular(0, 1, 1, 10);
        // Unlimited virtual + limited source needs %b in a name ...
        assert!(Mapping::new(mapping("f.h5", &unlim, &fixed)).is_err());
        let m = Mapping::new(mapping("f%b.h5", &unlim, &fixed)).unwrap();
        assert_eq!(m.kind, Kind::Printf { vdim: 0 });
        // ... and %b is refused anywhere else.
        assert!(Mapping::new(mapping("f%b.h5", &fixed, &fixed)).is_err());
        let src_unlim = regular(0, 1, UNLIMITED, 1);
        let m = Mapping::new(mapping("f.h5", &unlim, &src_unlim)).unwrap();
        assert_eq!(m.kind, Kind::Unlimited { vdim: 0, sdim: 0 });
        // An unlimited source into a limited virtual selection is refused.
        assert!(Mapping::new(mapping("f.h5", &fixed, &src_unlim)).is_err());
    }
}