oxmera_tensor/

tensor.rs

//! The tensor value: a layout over shared storage, with an optional
//! autograd tape node.

use std::sync::{Arc, Mutex};

use oxmera_core::layout::contiguous_strides;
use oxmera_core::{DType, Device, Error, Layout, Result, Shape, Strides};
use rand::SeedableRng;
use rand_distr::Distribution;

use crate::autograd::{AutogradMeta, GradFn};
use crate::backend::backend_for;
use crate::storage::{CpuStorage, Storage};

/// A tensor: shared storage viewed through a layout.
///
/// Cloning a tensor is cheap — it clones the layout and bumps the storage
/// refcount, never the data. View operations (`reshape`, `permute`,
/// `narrow`, …) produce new tensors over the same storage whenever the
/// layout arithmetic allows it.
#[derive(Clone)]
pub struct Tensor {
    storage: Arc<Storage>,
    layout: Layout,
    autograd: Option<Arc<AutogradMeta>>,
}

impl std::fmt::Debug for Tensor {
    /// Prints only the tensor's metadata — shape, dtype, device and
    /// whether it tracks gradients. Never the storage contents: a tensor
    /// can hold gigabytes, and a derived `Debug` dumped all of it into
    /// every log line and panic message.
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("Tensor")
            .field("shape", self.shape())
            .field("dtype", &self.dtype())
            .field("device", &self.device())
            .field("requires_grad", &self.requires_grad())
            .finish()
    }
}

impl Tensor {
    // ---- construction ---------------------------------------------------

    /// A tensor over existing storage with an explicit layout.
    ///
    /// Errors when the layout addresses elements outside the storage.
    pub fn from_storage(storage: Arc<Storage>, layout: Layout) -> Result<Self> {
        let available = storage_len(&storage) as isize;
        let Some((lo, hi)) = addressed_bounds(&layout) else {
            return Err(Error::InvalidArgument {
                op: "Tensor::from_storage",
                detail: "layout extents overflow the address space, so it cannot be proven \
                         in bounds"
                    .into(),
            });
        };
        if lo < 0 || hi > available {
            return Err(Error::InvalidArgument {
                op: "Tensor::from_storage",
                detail: format!(
                    "layout addresses [{lo}, {hi}) but storage holds {available} elements"
                ),
            });
        }
        Ok(Self {
            storage,
            layout,
            autograd: None,
        })
    }

    /// A contiguous CPU tensor holding `data` with shape `shape`.
    ///
    /// Errors when `data.len()` does not equal `shape.numel()`.
    pub fn from_vec_f32(data: Vec<f32>, shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::from_vec_f32")?;
        if data.len() != numel {
            return Err(Error::ShapeMismatch {
                expected: Shape::from([data.len()]),
                got: shape,
                op: "Tensor::from_vec_f32",
            });
        }
        Ok(Self {
            storage: Arc::new(Storage::from_f32_vec(data)),
            layout: Layout::contiguous(shape),
            autograd: None,
        })
    }

    /// A contiguous CPU tensor copying `data` with shape `shape`.
    pub fn from_slice(data: &[f32], shape: impl Into<Shape>) -> Result<Self> {
        Self::from_vec_f32(data.to_vec(), shape)
    }

    /// A contiguous CPU `F64` tensor holding `data`. `f64` tensors live on
    /// the CPU (the GPU backends carry `f32`); every op the CPU backend
    /// implements accepts them, and [`Tensor::to_dtype`] converts.
    pub fn from_vec_f64(data: Vec<f64>, shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::from_vec_f64")?;
        if data.len() != numel {
            return Err(Error::ShapeMismatch {
                expected: Shape::from([data.len()]),
                got: shape,
                op: "Tensor::from_vec_f64",
            });
        }
        Ok(Self {
            storage: Arc::new(Storage::from_f64_vec(data)),
            layout: Layout::contiguous(shape),
            autograd: None,
        })
    }

    /// A contiguous CPU `I64` tensor holding `data` (indices, targets).
    pub fn from_vec_i64(data: Vec<i64>, shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::from_vec_i64")?;
        if data.len() != numel {
            return Err(Error::ShapeMismatch {
                expected: Shape::from([data.len()]),
                got: shape,
                op: "Tensor::from_vec_i64",
            });
        }
        Ok(Self {
            storage: Arc::new(Storage::from_i64_vec(data)),
            layout: Layout::contiguous(shape),
            autograd: None,
        })
    }

    /// A CPU tensor of zeros.
    ///
    /// # Panics
    /// Panics if the shape's element count overflows `usize`. Build the
    /// shape from untrusted input through [`Tensor::try_zeros`] for a
    /// typed error instead.
    pub fn zeros(shape: impl Into<Shape>) -> Self {
        Self::try_zeros(shape).expect("shape element count overflows usize")
    }

    /// A CPU tensor of zeros, or [`Error::InvalidArgument`] when the
    /// shape's element count overflows `usize`.
    pub fn try_zeros(shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::try_zeros")?;
        Self::from_vec_f32(try_filled(numel, 0.0, "Tensor::try_zeros")?, shape)
    }

    /// A CPU tensor of ones.
    ///
    /// # Panics
    /// Panics if the shape's element count overflows `usize`; see
    /// [`Tensor::try_ones`].
    pub fn ones(shape: impl Into<Shape>) -> Self {
        Self::try_ones(shape).expect("shape element count overflows usize")
    }

    /// A CPU tensor of ones, or [`Error::InvalidArgument`] when the
    /// shape's element count overflows `usize`.
    pub fn try_ones(shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::try_ones")?;
        Self::from_vec_f32(try_filled(numel, 1.0, "Tensor::try_ones")?, shape)
    }

    /// A CPU tensor filled with `value`.
    ///
    /// # Panics
    /// Panics if the shape's element count overflows `usize`; see
    /// [`Tensor::try_full`].
    pub fn full(shape: impl Into<Shape>, value: f32) -> Self {
        Self::try_full(shape, value).expect("shape element count overflows usize")
    }

    /// A CPU tensor filled with `value`, or [`Error::InvalidArgument`]
    /// when the shape's element count overflows `usize`.
    pub fn try_full(shape: impl Into<Shape>, value: f32) -> Result<Self> {
        let shape = shape.into();
        let numel = checked_numel(&shape, "Tensor::try_full")?;
        Self::from_vec_f32(try_filled(numel, value, "Tensor::try_full")?, shape)
    }

    /// A rank-0 scalar tensor.
    pub fn scalar(value: f32) -> Self {
        Self::from_vec_f32(vec![value], Shape::from([])).expect("scalar always fits")
    }

    /// Standard-normal random CPU tensor, seeded from the OS.
    pub fn randn(shape: impl Into<Shape>) -> Self {
        Self::randn_with_seed(shape, rand::random())
    }

    /// Standard-normal random CPU tensor with a fixed seed, for
    /// reproducible tests and examples.
    pub fn randn_with_seed(shape: impl Into<Shape>, seed: u64) -> Self {
        let shape = shape.into();
        let numel = shape.numel();
        let mut rng = rand::rngs::StdRng::seed_from_u64(seed);
        let normal = rand_distr::StandardNormal;
        let data: Vec<f32> = (0..numel).map(|_| normal.sample(&mut rng)).collect();
        Self::from_vec_f32(data, shape).expect("lengths match by construction")
    }

    // ---- accessors -------------------------------------------------------

    /// The shape of this view.
    pub fn shape(&self) -> &Shape {
        &self.layout.shape
    }

    /// The dimension extents, outermost first.
    pub fn dims(&self) -> &[usize] {
        self.layout.shape.dims()
    }

    /// The rank (number of dimensions).
    pub fn ndim(&self) -> usize {
        self.layout.shape.ndim()
    }

    /// The total number of elements.
    pub fn numel(&self) -> usize {
        self.layout.shape.numel()
    }

    /// The full layout of this view.
    pub fn layout(&self) -> &Layout {
        &self.layout
    }

    /// The element type.
    pub fn dtype(&self) -> DType {
        self.storage.dtype()
    }

    /// The device the storage lives on.
    pub fn device(&self) -> Device {
        self.storage.device()
    }

    /// The shared storage behind this view.
    pub fn storage(&self) -> &Arc<Storage> {
        &self.storage
    }

    // ---- element access (CPU) --------------------------------------------

    /// The element at a logical index, as `f32`.
    ///
    /// Errors on rank mismatch, out-of-bounds, non-float dtype, or non-CPU
    /// storage.
    pub fn get_f32(&self, index: &[usize]) -> Result<f32> {
        let offset = self.layout.offset_of(index)?;
        Ok(self.storage.cpu()?.f32s()?[offset])
    }

    /// The element at a logical index, as `f64` (from an `F64` tensor).
    pub fn get_f64(&self, index: &[usize]) -> Result<f64> {
        let offset = self.layout.offset_of(index)?;
        Ok(self.storage.cpu()?.f64s()?[offset])
    }

    /// The element at a logical index, as `i64`.
    pub fn get_i64(&self, index: &[usize]) -> Result<i64> {
        let offset = self.layout.offset_of(index)?;
        Ok(self.storage.cpu()?.i64s()?[offset])
    }

    /// Every element in logical (row-major) order, as `f32`, from CPU
    /// storage.
    pub fn to_vec_f32(&self) -> Result<Vec<f32>> {
        let src = self.storage.cpu()?.f32s()?;
        Ok(gather_logical(src, &self.layout))
    }

    /// Every element in logical (row-major) order, as `f64`, from an `F64`
    /// CPU tensor (use [`Tensor::to_dtype`] first for an `f32` one).
    pub fn to_vec_f64(&self) -> Result<Vec<f64>> {
        let src = self.storage.cpu()?.f64s()?;
        Ok(gather_logical(src, &self.layout))
    }

    /// Every element in logical (row-major) order, as `i64`.
    pub fn to_vec_i64(&self) -> Result<Vec<i64>> {
        let src = self.storage.cpu()?.i64s()?;
        Ok(gather_logical(src, &self.layout))
    }

    // ---- views -----------------------------------------------------------

    fn view(&self, layout: Layout) -> Self {
        Self {
            storage: Arc::clone(&self.storage),
            layout,
            autograd: None,
        }
    }

    /// A view (or copy, when this view is not contiguous) with the same
    /// elements in a new shape.
    pub fn reshape(&self, shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        if checked_numel(&shape, "reshape")? != self.numel() {
            return Err(Error::ShapeMismatch {
                expected: self.shape().clone(),
                got: shape,
                op: "reshape",
            });
        }
        let base = if self.layout.is_contiguous() {
            self.clone()
        } else {
            self.contiguous_data()?
        };
        let layout = Layout {
            strides: contiguous_strides(&shape),
            shape,
            offset: base.layout.offset,
        };
        let out = base.view(layout);
        Ok(crate::ops::record_view(self, out, ViewKind::Reshape))
    }

    /// A view with dimensions reordered by `perm` (a permutation of
    /// `0..ndim`).
    pub fn permute(&self, perm: &[usize]) -> Result<Self> {
        let n = self.ndim();
        if perm.len() != n || {
            let mut seen = vec![false; n];
            perm.iter()
                .any(|&p| p >= n || std::mem::replace(&mut seen[p], true))
        } {
            return Err(Error::InvalidArgument {
                op: "permute",
                detail: format!("{perm:?} is not a permutation of 0..{n}"),
            });
        }
        let dims = self.dims();
        let strides = self.layout.strides.values();
        let new_dims: Vec<usize> = perm.iter().map(|&p| dims[p]).collect();
        let new_strides: Vec<isize> = perm.iter().map(|&p| strides[p]).collect();
        let layout = Layout {
            shape: Shape::new(new_dims),
            strides: Strides::new(new_strides),
            offset: self.layout.offset,
        };
        let out = self.view(layout);
        Ok(crate::ops::record_view(
            self,
            out,
            ViewKind::Permute(perm.to_vec()),
        ))
    }

    /// A view with dimensions `d0` and `d1` swapped.
    pub fn transpose(&self, d0: usize, d1: usize) -> Result<Self> {
        let mut perm: Vec<usize> = (0..self.ndim()).collect();
        if d0 >= perm.len() || d1 >= perm.len() {
            return Err(Error::InvalidArgument {
                op: "transpose",
                detail: format!("dims ({d0}, {d1}) out of range for rank {}", perm.len()),
            });
        }
        perm.swap(d0, d1);
        self.permute(&perm)
    }

    /// The matrix transpose: the last two dimensions swapped.
    pub fn t(&self) -> Result<Self> {
        let n = self.ndim();
        if n < 2 {
            return Err(Error::InvalidArgument {
                op: "t",
                detail: format!("needs rank >= 2, got {n}"),
            });
        }
        self.transpose(n - 2, n - 1)
    }

    /// A view of `len` elements of dimension `dim` starting at `start`.
    pub fn narrow(&self, dim: usize, start: usize, len: usize) -> Result<Self> {
        let dims = self.dims();
        let end = start.checked_add(len);
        if dim >= dims.len() || end.is_none_or(|e| e > dims[dim]) {
            return Err(Error::InvalidArgument {
                op: "narrow",
                detail: format!(
                    "dim {dim}, start {start} len {len} against shape {:?}",
                    self.shape()
                ),
            });
        }
        let mut new_dims = dims.to_vec();
        new_dims[dim] = len;
        let strides = self.layout.strides.values().to_vec();
        let offset = (self.layout.offset as isize + start as isize * strides[dim]) as usize;
        let layout = Layout {
            shape: Shape::new(new_dims),
            strides: Strides::new(strides),
            offset,
        };
        let out = self.view(layout);
        Ok(crate::ops::record_view(
            self,
            out,
            ViewKind::Narrow { dim, start, len },
        ))
    }

    /// A view of `range` along `dim` — sugar over [`Tensor::narrow`].
    pub fn slice(&self, dim: usize, range: std::ops::Range<usize>) -> Result<Self> {
        if range.end < range.start {
            // Saturating to an empty view turned a caller mistake into a
            // silently empty tensor that fails much later.
            return Err(Error::InvalidArgument {
                op: "slice",
                detail: format!("range {}..{} is reversed", range.start, range.end),
            });
        }
        self.narrow(dim, range.start, range.end - range.start)
    }

    /// A zero-copy broadcast view to `shape` (stride 0 on expanded axes).
    pub fn broadcast_to(&self, shape: impl Into<Shape>) -> Result<Self> {
        let shape = shape.into();
        checked_numel(&shape, "broadcast_to")?;
        let layout = broadcast_layout(&self.layout, &shape)?;
        let out = self.view(layout);
        Ok(crate::ops::record_view(self, out, ViewKind::Broadcast))
    }

    /// A broadcast view that records nothing on the tape — backend
    /// plumbing; prefer [`Tensor::broadcast_to`] in user code.
    pub fn broadcast_view(&self, shape: &Shape) -> Result<Self> {
        checked_numel(shape, "broadcast_view")?;
        let layout = broadcast_layout(&self.layout, shape)?;
        Ok(self.view(layout))
    }

    /// A view with a new size-1 dimension inserted at `dim`.
    pub fn unsqueeze(&self, dim: usize) -> Result<Self> {
        let mut dims = self.dims().to_vec();
        if dim > dims.len() {
            return Err(Error::InvalidArgument {
                op: "unsqueeze",
                detail: format!("dim {dim} out of range for rank {}", dims.len()),
            });
        }
        dims.insert(dim, 1);
        let mut strides = self.layout.strides.values().to_vec();
        strides.insert(dim, 0);
        let layout = Layout {
            shape: Shape::new(dims),
            strides: Strides::new(strides),
            offset: self.layout.offset,
        };
        let out = self.view(layout);
        Ok(crate::ops::record_view(self, out, ViewKind::Reshape))
    }

    /// This tensor's elements, in logical order, in fresh contiguous
    /// storage on the same device. A no-op clone when already contiguous.
    pub fn contiguous(&self) -> Result<Self> {
        if self.layout.is_contiguous() && self.layout.offset == 0 {
            return Ok(self.clone());
        }
        let out = self.contiguous_data()?;
        Ok(crate::ops::record_view(self, out, ViewKind::Contiguous))
    }

    /// The contiguous copy without autograd recording — backend plumbing;
    /// prefer [`Tensor::contiguous`] in user code.
    pub fn contiguous_untracked(&self) -> Result<Self> {
        self.contiguous_data()
    }

    /// The contiguous copy without autograd recording (plumbing).
    pub(crate) fn contiguous_data_crate(&self) -> Result<Self> {
        self.contiguous_data()
    }

    /// The contiguous copy without autograd recording (plumbing).
    pub(crate) fn contiguous_data(&self) -> Result<Self> {
        match self.device() {
            Device::Cpu => {
                let shape = self.shape().clone();
                match self.storage.cpu()? {
                    CpuStorage::F32(_) => Tensor::from_vec_f32(self.to_vec_f32()?, shape),
                    CpuStorage::F64(_) => Tensor::from_vec_f64(self.to_vec_f64()?, shape),
                    CpuStorage::I64(_) => Tensor::from_vec_i64(self.to_vec_i64()?, shape),
                    CpuStorage::U8(_) => Err(Error::UnsupportedDType {
                        dtype: DType::U8,
                        op: "contiguous",
                    }),
                }
            }
            device => backend_for(device)?.contiguous(self),
        }
    }

    // ---- autograd surface -------------------------------------------------

    /// Mark (or unmark) this tensor as a gradient-accumulating leaf, in
    /// place, returning it for chaining.
    pub fn requires_grad_(mut self, requires: bool) -> Self {
        match (&self.autograd, requires) {
            (Some(meta), _) if meta.grad_fn.is_none() => {
                // Leaf: rebuild the meta with the new flag.
                self.autograd = requires.then(|| {
                    Arc::new(AutogradMeta {
                        requires_grad: true,
                        grad: Mutex::new(None),
                        grad_fn: None,
                    })
                });
            }
            (Some(_), true) => { /* non-leaf already tracked; nothing to do */ }
            (Some(_), false) => self.autograd = None,
            (None, true) => {
                self.autograd = Some(Arc::new(AutogradMeta {
                    requires_grad: true,
                    grad: Mutex::new(None),
                    grad_fn: None,
                }));
            }
            (None, false) => {}
        }
        self
    }

    /// Whether gradients **accumulate on this tensor** during `backward`
    /// — true for leaves marked with [`requires_grad_`](Self::requires_grad_).
    ///
    /// This answers a narrower question than PyTorch's `requires_grad`:
    /// a tensor *computed from* such a leaf is on the tape but does not
    /// accumulate a gradient of its own, so it reports `false` here. Ask
    /// [`is_tracked`](Self::is_tracked) for "is this on the graph at all".
    pub fn requires_grad(&self) -> bool {
        self.autograd.as_ref().is_some_and(|m| m.requires_grad)
    }

    /// Whether this tensor participates in the autograd tape at all —
    /// true for a leaf that requires grad and for anything computed from
    /// one while recording was enabled; false for constants and for
    /// everything produced under [`no_grad`](crate::autograd::no_grad).
    ///
    /// This is the predicate that observes `no_grad`:
    ///
    /// ```
    /// use oxmera_tensor::tensor::Tensor;
    /// use oxmera_tensor::autograd::no_grad;
    ///
    /// let a = Tensor::from_slice(&[1.0, 2.0], [2]).unwrap().requires_grad_(true);
    /// assert!(a.mul_scalar(3.0).unwrap().is_tracked());
    /// assert!(!no_grad(|| a.mul_scalar(3.0).unwrap()).is_tracked());
    /// ```
    pub fn is_tracked(&self) -> bool {
        self.autograd.is_some()
    }

    /// The accumulated gradient, if a backward pass has produced one.
    pub fn grad(&self) -> Option<Tensor> {
        self.autograd
            .as_ref()
            .and_then(|m| m.grad.lock().expect("grad mutex poisoned").clone())
    }

    /// Clear this tensor's accumulated gradient.
    pub fn zero_grad(&self) {
        if let Some(meta) = &self.autograd {
            *meta.grad.lock().expect("grad mutex poisoned") = None;
        }
    }

    /// The same view without any tape connection.
    pub fn detach(&self) -> Self {
        Self {
            storage: Arc::clone(&self.storage),
            layout: self.layout.clone(),
            autograd: None,
        }
    }

    /// Propagate gradients from this scalar through the recorded tape.
    ///
    /// Errors when the tensor is not a scalar; use
    /// [`Tensor::backward_with`] to seed a non-scalar output.
    pub fn backward(&self) -> Result<()> {
        if self.numel() != 1 {
            return Err(Error::InvalidArgument {
                op: "backward",
                detail: format!(
                    "output has {} elements; seed a non-scalar with backward_with",
                    self.numel()
                ),
            });
        }
        // The seed lives where the loss lives: a CPU seed against a
        // device-resident graph failed at the first VJP with a
        // DeviceMismatch (found by the CUDA backend's end-to-end test, and
        // latent on Metal).
        let seed = Tensor::ones(self.shape().clone())
            .to_dtype(self.dtype())?
            .to_device(self.device())?;
        self.backward_with(seed)
    }

    /// Propagate gradients seeding this tensor's gradient with `seed`.
    pub fn backward_with(&self, seed: Tensor) -> Result<()> {
        let seed = if seed.device() == self.device() {
            seed
        } else {
            seed.to_device(self.device())?
        };
        crate::autograd::run_backward(self, seed)
    }

    pub(crate) fn autograd_meta(&self) -> Option<Arc<AutogradMeta>> {
        self.autograd.clone()
    }

    /// Attach a tape node to this tensor (used by the op layer).
    /// Detach this tensor's tape node and hand it back, leaving the tensor a
    /// leaf. `AutogradMeta::drop` uses this to dismantle a deep tape
    /// iteratively rather than recursing once per node.
    pub(crate) fn take_autograd(&mut self) -> Option<Arc<AutogradMeta>> {
        self.autograd.take()
    }

    pub(crate) fn with_grad_fn(mut self, grad_fn: GradFn) -> Self {
        self.autograd = Some(Arc::new(AutogradMeta {
            requires_grad: false,
            grad: Mutex::new(None),
            grad_fn: Some(grad_fn),
        }));
        self
    }
}

/// The internal view taxonomy the op layer uses to build view VJPs.
#[derive(Debug, Clone)]
pub(crate) enum ViewKind {
    /// Reshape/unsqueeze: gradient reshapes back to the input shape.
    Reshape,
    /// Permute by this permutation: gradient permutes by the inverse.
    Permute(Vec<usize>),
    /// Narrow: gradient scatters back into zeros of the input shape.
    Narrow {
        /// Narrowed dimension.
        dim: usize,
        /// Range start.
        start: usize,
        /// Range length.
        len: usize,
    },
    /// Broadcast view: gradient sum-reduces back to the input shape.
    Broadcast,
    /// Contiguous copy: gradient passes through (reshaped if needed).
    Contiguous,
}

fn storage_len(storage: &Storage) -> usize {
    match storage.data() {
        crate::storage::StorageData::Cpu(c) => c.len(),
        #[cfg(target_os = "macos")]
        crate::storage::StorageData::Metal(b) => {
            b.buffer().length() as usize / storage.dtype().size_in_bytes()
        }
        crate::storage::StorageData::Opaque(b) => b.len(),
    }
}

/// The half-open range of storage indices a layout can address, as
/// `(min, max_exclusive)`. Negative strides lower the minimum below the
/// offset, so a valid layout needs `min >= 0` as well as
/// `max_exclusive <= storage length`; the pre-0.5 check accounted for
/// positive strides only and let an underflowing negative-stride layout
/// through to a panic on the first read.
fn addressed_bounds(layout: &Layout) -> Option<(isize, isize)> {
    if layout.shape.checked_numel()? == 0 {
        return Some((0, 0));
    }
    let base = isize::try_from(layout.offset).ok()?;
    let mut lo = base;
    let mut hi = base;
    for (&d, &s) in layout.shape.dims().iter().zip(layout.strides.values()) {
        if d > 1 {
            let span = isize::try_from(d - 1).ok()?.checked_mul(s)?;
            if s >= 0 {
                hi = hi.checked_add(span)?;
            } else {
                lo = lo.checked_add(span)?;
            }
        }
    }
    Some((lo, hi.checked_add(1)?))
}

/// Broadcast `layout` to `target`, stride 0 on expanded axes.
pub(crate) fn broadcast_layout(layout: &Layout, target: &Shape) -> Result<Layout> {
    let src = layout.shape.dims();
    let dst = target.dims();
    if dst.len() < src.len() {
        return Err(Error::BroadcastIncompatible {
            lhs: layout.shape.clone(),
            rhs: target.clone(),
        });
    }
    let lead = dst.len() - src.len();
    let mut strides = vec![0isize; dst.len()];
    for i in 0..src.len() {
        let (s, d) = (src[i], dst[lead + i]);
        if s == d {
            strides[lead + i] = layout.strides.values()[i];
        } else if s == 1 {
            strides[lead + i] = 0;
        } else {
            return Err(Error::BroadcastIncompatible {
                lhs: layout.shape.clone(),
                rhs: target.clone(),
            });
        }
    }
    Ok(Layout {
        shape: target.clone(),
        strides: Strides::new(strides),
        offset: layout.offset,
    })
}

/// Gather a strided view's elements into logical row-major order.
pub(crate) fn gather_logical<T: Copy>(src: &[T], layout: &Layout) -> Vec<T> {
    let numel = layout.shape.numel();
    let mut out = Vec::with_capacity(numel);
    if numel == 0 {
        return out;
    }
    let dims = layout.shape.dims();
    if layout.is_contiguous() {
        let start = layout.offset;
        out.extend_from_slice(&src[start..start + numel]);
        return out;
    }
    let strides = layout.strides.values();
    // Row fast path: when the innermost stride is 1 the row is a slice
    // copy; when it is 0 (a broadcast dimension being materialized) the
    // row is one value repeated. Only a genuinely strided innermost
    // dimension — a transposed view — falls through to the odometer.
    let ndim = dims.len();
    let inner = dims[ndim - 1];
    let inner_stride = strides[ndim - 1];
    if inner > 0 && (inner_stride == 0 || inner_stride == 1) {
        let outer = Layout {
            shape: Shape::new(dims[..ndim - 1].to_vec()),
            strides: Strides::new(strides[..ndim - 1].to_vec()),
            offset: layout.offset,
        };
        let mut walker = crate::cpu_iter::OffsetWalker::at(&outer, 0);
        for _ in 0..numel / inner {
            let base = walker.next_offset();
            if inner_stride == 1 {
                out.extend_from_slice(&src[base..base + inner]);
            } else {
                out.resize(out.len() + inner, src[base]);
            }
        }
        return out;
    }
    let mut index = vec![0usize; dims.len()];
    let mut offset = layout.offset as isize;
    loop {
        out.push(src[offset as usize]);
        // Odometer increment, last dimension fastest, offset updated
        // incrementally.
        let mut d = dims.len();
        loop {
            if d == 0 {
                return out;
            }
            d -= 1;
            index[d] += 1;
            offset += strides[d];
            if index[d] < dims[d] {
                break;
            }
            offset -= dims[d] as isize * strides[d];
            index[d] = 0;
        }
        if out.len() == numel {
            return out;
        }
    }
}

/// The element count of a caller-supplied shape, or a typed error when it
/// does not fit in `usize` (see [`Shape::checked_numel`]).
/// `numel` copies of `value`, reporting an allocation failure as a typed
/// error. `vec![value; numel]` *aborts the process* when the allocator
/// refuses, which a `try_` constructor must never do — that is the whole
/// reason the caller reached for the fallible form.
fn try_filled(numel: usize, value: f32, op: &'static str) -> Result<Vec<f32>> {
    let mut data: Vec<f32> = Vec::new();
    data.try_reserve_exact(numel)
        .map_err(|_| Error::InvalidArgument {
            op,
            detail: format!("cannot allocate {numel} f32 elements"),
        })?;
    data.resize(numel, value);
    Ok(data)
}

fn checked_numel(shape: &Shape, op: &'static str) -> Result<usize> {
    shape.checked_numel().ok_or_else(|| Error::InvalidArgument {
        op,
        detail: format!(
            "shape {:?} has more elements than fit in usize",
            shape.dims()
        ),
    })
}