tensor4all_core/defaults/

idx_tensor.rs

use crate::defaults::DynIndex;
use crate::index_like::IndexLike;
use crate::index_ops::{common_ind_positions, prepare_contraction, prepare_contraction_pairs};
use crate::tensor_like::LinearizationOrder;
use crate::AnyScalar;
use anyhow::{Context, Result};
use num_complex::{Complex32, Complex64};
use num_traits::Zero;
use rand::Rng;
use rand_distr::{Distribution, StandardNormal};
use std::cell::RefCell;
use std::cmp::Reverse;
use std::collections::{HashMap, HashSet};
use std::env;
use std::sync::{Arc, OnceLock};
use std::time::{Duration, Instant};
use tenferro::{
    DType, DotGeneralConfig, Tensor as NativeTensor, TensorRead, TensorValue, TensorView,
};
use tenferro_ad::{extension::adopt_untracked_eager_value, EagerRuntime, EagerTensor};
use tenferro_einsum::{EagerEinsumExt, EinsumSubscripts};
use tenferro_linalg::{EagerTensorLinalgExt, RankRevealingQrOptions};
use tensor4all_tensorbackend::{
    contract_native_tensor, default_eager_ctx, dense_native_tensor_from_col_major,
    dense_native_tensor_from_col_major_owned, diag_native_tensor_from_col_major,
    native_tensor_primal_to_diag, storage_payload_native_read_input, storage_to_native_tensor,
    ExecutionContext, NativeTensorReadInput, TensorElement,
};
use tensor4all_tensorbackend::{
    src_error_estimate as backend_src_error_estimate,
    src_error_estimate_general as backend_src_error_estimate_general, Matrix,
};
use tensor4all_tensorbackend::{IncrementalQr, IncrementalQrScalar};
use tensor4all_tensorbackend::{Storage, StorageKind};

use super::contract::{ContractionOptions, PairwiseContractionOptions};
use super::structured_contraction::{
    normalize_payload_read_for_roots, storage_from_payload_native, storage_payload_native,
    OperandLayout, StructuredContractionPlan, StructuredContractionSpec,
};

fn conjugate_eager(
    inner: &EagerTensor,
) -> std::result::Result<EagerTensor, Arc<dyn std::error::Error + Send + Sync + 'static>> {
    inner.conj().map_err(|source| Arc::new(source) as _)
}

#[derive(Debug, Default, Clone)]
struct PairwiseContractProfileEntry {
    calls: usize,
    total_time: Duration,
    total_bytes: usize,
}

/// Hermitian eigendecomposition of a rank-2 [`IdxTensor`].
/// Eigenvectors are returned as a rank-2 tensor whose first index is the input
/// matrix row index and whose second index labels eigenvector columns. The
/// eigenvalues are detached primal values intended for nonsmooth selection
/// logic such as truncation cutoffs.
/// # Examples
/// ```
/// use tensor4all_core::{DynIndex, IdxTensor};
/// let row = DynIndex::new_dyn(2);
/// let col = DynIndex::new_dyn(2);
/// let matrix = IdxTensor::from_dense(
///     vec![row.clone(), col],
///     vec![1.0_f64, 0.0, 0.0, 2.0],
/// ).unwrap();
/// let decomp = matrix.hermitian_eigendecomposition(1.0e-12).unwrap();
/// assert_eq!(decomp.eigenvalues, vec![1.0, 2.0]);
/// assert_eq!(
///     decomp.eigenvectors.indices(),
///     &[row, decomp.eigenvector_index.clone()]
/// );
/// ```
#[derive(Debug, Clone)]
pub struct TensorHermitianEigendecomposition {
    /// Real eigenvalues in backend Hermitian eigensolver order.
    pub eigenvalues: Vec<f64>,
    /// Eigenvector matrix with one eigenvector in each column.
    pub eigenvectors: IdxTensor,
    /// Index labeling the eigenvector columns.
    pub eigenvector_index: DynIndex,
}

thread_local! {
    static PAIRWISE_CONTRACT_PROFILE_STATE: RefCell<HashMap<&'static str, PairwiseContractProfileEntry>> =
        RefCell::new(HashMap::new());
}

fn pairwise_contract_profile_enabled() -> bool {
    static ENABLED: OnceLock<bool> = OnceLock::new();
    *ENABLED.get_or_init(|| env::var("T4A_PROFILE_PAIRWISE_CONTRACT").is_ok())
}

fn record_pairwise_contract_profile(section: &'static str, elapsed: Duration) {
    if !pairwise_contract_profile_enabled() {
        return;
    }
    PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
        let mut state = state.borrow_mut();
        let entry = state.entry(section).or_default();
        entry.calls += 1;
        entry.total_time += elapsed;
    });
}

fn record_pairwise_contract_profile_bytes(section: &'static str, bytes: usize) {
    if !pairwise_contract_profile_enabled() {
        return;
    }
    PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
        let mut state = state.borrow_mut();
        let entry = state.entry(section).or_default();
        entry.total_bytes += bytes;
    });
}

fn profile_pairwise_contract_section<T>(section: &'static str, f: impl FnOnce() -> T) -> T {
    if !pairwise_contract_profile_enabled() {
        return f();
    }
    let started = Instant::now();
    let result = f();
    record_pairwise_contract_profile(section, started.elapsed());
    result
}

/// Reset the aggregated pairwise `IdxTensor` contraction profile.
pub fn reset_pairwise_contract_profile() {
    PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| state.borrow_mut().clear());
}

/// Print and clear the aggregated pairwise `IdxTensor` contraction profile.
pub fn print_and_reset_pairwise_contract_profile() {
    if !pairwise_contract_profile_enabled() {
        return;
    }
    PAIRWISE_CONTRACT_PROFILE_STATE.with(|state| {
        let mut entries: Vec<_> = state
            .borrow()
            .iter()
            .map(|(section, entry)| (*section, entry.clone()))
            .collect();
        state.borrow_mut().clear();
        entries.sort_by_key(|(_, entry)| Reverse(entry.total_time));

        eprintln!("=== IdxTensor pairwise contract profile ===");
        for (section, entry) in entries {
            let per_call_us = if entry.calls == 0 {
                0.0
            } else {
                entry.total_time.as_secs_f64() * 1.0e6 / entry.calls as f64
            };
            eprintln!(
                "{section}: calls={} total={:.6}ms per_call={:.3}us bytes={}",
                entry.calls,
                entry.total_time.as_secs_f64() * 1.0e3,
                per_call_us,
                entry.total_bytes,
            );
        }
    });
}

fn tensor_profile_bytes(dtype: DType, shape: &[usize]) -> usize {
    let element_size = match dtype {
        DType::F32 => 4,
        DType::F64 => 8,
        DType::C32 => 8,
        DType::C64 => 16,
        DType::I32 => 4,
        DType::I64 => 8,
        DType::Bool => 1,
    };
    shape
        .iter()
        .try_fold(1usize, |bytes, &dim| bytes.checked_mul(dim))
        .and_then(|elements| elements.checked_mul(element_size))
        .unwrap_or(usize::MAX)
}

/// Trait for scalar types that can generate random values from a standard
/// normal distribution.
/// This enables the generic [`IdxTensor::random`] constructor.
pub trait RandomScalar: TensorElement {
    /// Generate a random value from the standard normal distribution.
    fn random_value<R: Rng>(rng: &mut R) -> Self;
}

impl RandomScalar for f64 {
    fn random_value<R: Rng>(rng: &mut R) -> Self {
        StandardNormal.sample(rng)
    }
}

impl RandomScalar for Complex64 {
    fn random_value<R: Rng>(rng: &mut R) -> Self {
        Complex64::new(StandardNormal.sample(rng), StandardNormal.sample(rng))
    }
}

/// Compute the permutation array from original indices to new indices.
/// This function finds the mapping from new indices to original indices by
/// matching index IDs. The result is a permutation array `perm` such that
/// `new_indices[i]` corresponds to `original_indices[perm[i]]`.
/// # Arguments
/// * `original_indices` - The original indices in their current order
/// * `new_indices` - The desired new indices order (must be a permutation of original_indices)
/// # Returns
/// A `Vec<usize>` representing the permutation: `perm[i]` is the position in
/// `original_indices` of the index that should be at position `i` in `new_indices`.
/// # Errors
/// Returns an error when `new_order` contains indices not present in
/// `original` (a missing-index failure) or the two lists differ in length
/// (a length mismatch).
/// # Example
/// ```
/// use tensor4all_core::tensor::compute_permutation_from_indices;
/// use tensor4all_core::DynIndex;
/// let i = DynIndex::new_dyn(2);
/// let j = DynIndex::new_dyn(3);
/// let original = vec![i.clone(), j.clone()];
/// let new_order = vec![j.clone(), i.clone()];
/// let perm = compute_permutation_from_indices(&original, &new_order).unwrap();
/// assert_eq!(perm, vec![1, 0]);  // j is at position 1, i is at position 0
/// ```
pub fn compute_permutation_from_indices(
    original_indices: &[DynIndex],
    new_indices: &[DynIndex],
) -> std::result::Result<Vec<usize>, IdxTensorError> {
    if !(new_indices.len() == original_indices.len()) {
        return Err(
            anyhow::anyhow!("new_indices length must match original_indices length").into(),
        );
    };

    let mut perm = Vec::with_capacity(new_indices.len());
    let mut used = std::collections::HashSet::new();

    for new_idx in new_indices {
        // Find the position of this index in the original indices
        // DynIndex implements Eq, so we can compare directly
        let pos = original_indices
            .iter()
            .position(|old_idx| old_idx == new_idx)
            .ok_or_else(|| {
                anyhow::anyhow!("new_indices must be a permutation of original_indices")
            })?;

        if !(used.insert(pos)) {
            return Err(anyhow::anyhow!("duplicate index in new_indices").into());
        };
        perm.push(pos);
    }

    Ok(perm)
}

/// Compact structured payload kept in the authoritative eager representation.
/// The payload may use any supported eager dtype (`f32`, `f64`, `c32`, or
/// `c64`) and may be either tracked or untracked. Tracking is a property of
/// `payload`, never of the presence of this metadata container.
#[derive(Clone)]
pub(crate) struct StructuredPayload {
    payload: Arc<EagerTensor>,
    payload_dims: Vec<usize>,
    axis_classes: Vec<usize>,
}

/// Error returned when [`IdxTensor::storage`] or
/// [`IdxTensor::to_storage`] cannot produce a compact `f64`/`Complex64`
/// storage snapshot from the authoritative payload.
/// Backend diagnostics remain available through [`std::error::Error::source`]
/// instead of being erased into a display string. The error is cloneable so a
/// deferred failure can be retained by cloned tensors without rebuilding a
/// detached primal value.
/// # Examples
/// ```
/// use std::error::Error;
/// use std::sync::Arc;
/// use tensor4all_core::TensorStorageError;
/// let error = TensorStorageError::Materialization {
///     source: Arc::new(std::io::Error::other("backend unavailable")),
/// };
/// assert!(error.source().is_some());
/// assert!(error.to_string().contains("backend unavailable"));
/// ```
#[derive(Debug, Clone, thiserror::Error)]
pub enum TensorStorageError {
    /// The eager or structured payload could not be converted to compact storage.
    #[error("failed to materialize IdxTensor storage: {source}")]
    Materialization {
        /// Original diagnostic returned by the backend or storage conversion seam.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
    /// An eager payload uses a scalar dtype that compact [`Storage`] cannot hold.
    #[error(
        "compact IdxTensor storage does not support dtype {dtype}; the eager payload remains authoritative"
    )]
    UnsupportedDtype {
        /// Native scalar dtype retained by the eager representation.
        dtype: &'static str,
    },
    /// An eager conjugation operation failed and was deferred by the infallible
    /// [`IdxTensor::conj`] API.
    #[error("failed to conjugate IdxTensor storage: {source}")]
    Conjugation {
        /// Original diagnostic returned by the eager AD backend.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
}

/// Errors returned by the fallible numerical and comparison methods on
/// [`IdxTensor`].
/// The enum is intentionally owned by `tensor4all-core`: callers can match
/// storage, shape, scalar, subtraction, and invalid-value failures without
/// depending on the internal `anyhow` plumbing. Wrapped backend diagnostics
/// retain their complete [`std::error::Error::source`] chain.
///
/// # Remedies
/// - Storage failures: check the operation against the storage kind
///   (structured vs dense) and dtype before calling; the eager payload may
///   remain authoritative for unsupported dtypes.
/// - Shape/index failures: validate indices and dimensions at the call site
///   (`dims`, `indices`, external index sets) before the operation.
/// - NaN/invalid-value failures: inspect the payload for non-finite entries
///   before numerical comparisons.
/// - Backend failures: the wrapped source chain identifies the backend stage;
///   re-run with the backend diagnostic visible (see `source`).
/// # Examples
/// ```
/// use tensor4all_core::IdxTensorError;
/// let error = IdxTensorError::NaNInput {
///     operation: "norm_squared",
/// };
/// assert!(error.to_string().contains("NaN"));
/// ```
#[derive(Debug, Clone, thiserror::Error)]
pub enum IdxTensorError {
    /// Compact storage or deferred storage materialization failed.
    #[error("IdxTensor storage operation failed: {source}")]
    Storage {
        /// Original storage diagnostic, including its backend source chain.
        #[source]
        source: TensorStorageError,
    },
    /// A native eager payload could not be materialized for a numerical operation.
    #[error("IdxTensor materialization failed: {source}")]
    Materialization {
        /// Original backend or eager-runtime diagnostic.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
    /// A rank-zero scalar could not be extracted from a reduction result.
    #[error("IdxTensor scalar extraction failed: {source}")]
    ScalarExtraction {
        /// Original scalar-wrapper or backend diagnostic.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
    /// The reduction result has a scalar dtype that this real-valued operation
    /// cannot interpret.
    #[error("IdxTensor scalar type mismatch: expected {expected}, got {actual}")]
    ScalarTypeMismatch {
        /// Scalar dtype required by the operation.
        expected: &'static str,
        /// Scalar dtype returned by the reduction.
        actual: String,
    },
    /// Tensor shapes, index spaces, or dimension metadata cannot be aligned
    /// for an operation (comparison, index replacement, or other shape-sensitive
    /// transformations).
    #[error("IdxTensor shape mismatch during {operation}: expected {expected}, got {actual}")]
    ShapeMismatch {
        /// Operation that attempted the alignment.
        operation: &'static str,
        /// Expected index/dimension description.
        expected: String,
        /// Actual index/dimension description.
        actual: String,
    },
    /// Tensor subtraction failed while evaluating a comparison.
    #[error("IdxTensor subtraction failed: {source}")]
    Subtraction {
        /// Original arithmetic or backend diagnostic.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
    /// An input contained a NaN and the operation rejected it rather than
    /// silently converting it to zero.
    #[error("IdxTensor {operation} received NaN input")]
    NaNInput {
        /// Numerical operation that observed the NaN.
        operation: &'static str,
    },
    /// A comparison tolerance was NaN, infinite, or negative.
    #[error("IdxTensor tolerance {name} is invalid: {value}")]
    InvalidTolerance {
        /// Name of the invalid tolerance.
        name: &'static str,
        /// Supplied tolerance value.
        value: f64,
    },
    /// Another eager tensor operation failed while preparing a comparison.
    #[error("IdxTensor {operation} failed: {source}")]
    Operation {
        /// Name of the eager operation that failed.
        operation: &'static str,
        /// Original backend or tensor diagnostic.
        #[source]
        source: Arc<dyn std::error::Error + Send + Sync + 'static>,
    },
}

impl From<anyhow::Error> for IdxTensorError {
    fn from(source: anyhow::Error) -> Self {
        Self::operation("IdxTensor", source)
    }
}

impl From<TensorStorageError> for IdxTensorError {
    fn from(source: TensorStorageError) -> Self {
        Self::Storage { source }
    }
}

impl From<tenferro_ad::Error> for IdxTensorError {
    fn from(source: tenferro_ad::Error) -> Self {
        Self::Materialization {
            source: Arc::from(anyhow::Error::new(source).into_boxed_dyn_error()),
        }
    }
}

impl From<tensor4all_tensorbackend::EagerContextError> for IdxTensorError {
    fn from(source: tensor4all_tensorbackend::EagerContextError) -> Self {
        Self::Materialization {
            source: Arc::from(anyhow::Error::new(source).into_boxed_dyn_error()),
        }
    }
}

impl From<tensor4all_tensorbackend::BridgeError> for IdxTensorError {
    fn from(source: tensor4all_tensorbackend::BridgeError) -> Self {
        Self::Materialization {
            source: Arc::new(source),
        }
    }
}

impl IdxTensorError {
    fn boxed(error: anyhow::Error) -> Arc<dyn std::error::Error + Send + Sync + 'static> {
        Arc::from(error.into_boxed_dyn_error())
    }

    fn materialization(error: anyhow::Error) -> Self {
        Self::Materialization {
            source: Self::boxed(error),
        }
    }

    fn scalar_extraction(error: anyhow::Error) -> Self {
        Self::ScalarExtraction {
            source: Self::boxed(error),
        }
    }

    fn operation(operation: &'static str, error: anyhow::Error) -> Self {
        Self::Operation {
            operation,
            source: Self::boxed(error),
        }
    }
}

#[derive(Clone)]
pub(crate) enum IdxTensorStorage {
    Materialized(Arc<Storage>),
    Eager {
        inner: Arc<EagerTensor>,
        axis_classes: Vec<usize>,
    },
    /// One authoritative compact eager payload and its logical layout.
    Compact(Arc<StructuredPayload>),
    /// A storage representation whose eager operation failed before the
    /// infallible tensor API could return an error.
    Deferred {
        source: Box<Self>,
        error: Arc<TensorStorageError>,
    },
}

impl IdxTensorStorage {
    fn from_storage(storage: Arc<Storage>) -> Self {
        Self::Materialized(storage)
    }

    fn from_eager_dense(inner: EagerTensor, rank: usize) -> Self {
        Self::Eager {
            inner: Arc::new(inner),
            axis_classes: IdxTensor::dense_axis_classes(rank),
        }
    }

    fn eager(&self) -> Option<&EagerTensor> {
        match self {
            Self::Materialized(_) => None,
            Self::Eager { inner, .. } => Some(inner.as_ref()),
            Self::Compact(payload) => Some(payload.payload.as_ref()),
            Self::Deferred { source, .. } => source.eager(),
        }
    }

    fn deferred_error(&self) -> Option<&TensorStorageError> {
        match self {
            Self::Deferred { error, .. } => Some(error.as_ref()),
            _ => None,
        }
    }

    fn with_deferred_error(self, error: TensorStorageError) -> Self {
        if self.deferred_error().is_some() {
            self
        } else {
            Self::Deferred {
                source: Box::new(self),
                error: Arc::new(error),
            }
        }
    }

    fn axis_classes(&self) -> &[usize] {
        match self {
            Self::Materialized(storage) => storage.axis_classes(),
            Self::Eager { axis_classes, .. } => axis_classes,
            Self::Compact(payload) => &payload.axis_classes,
            Self::Deferred { source, .. } => source.axis_classes(),
        }
    }

    fn payload_dims(&self) -> &[usize] {
        match self {
            Self::Materialized(storage) => storage.payload_dims(),
            Self::Eager { inner, .. } => inner.shape(),
            Self::Compact(payload) => &payload.payload_dims,
            Self::Deferred { source, .. } => source.payload_dims(),
        }
    }

    fn payload_strides_vec(&self) -> Vec<isize> {
        match self {
            Self::Materialized(storage) => storage.payload_strides().to_vec(),
            Self::Eager { inner, .. } => {
                IdxTensor::col_major_strides(inner.shape()).unwrap_or_default()
            }
            Self::Compact(payload) => {
                IdxTensor::col_major_strides(&payload.payload_dims).unwrap_or_default()
            }
            Self::Deferred { source, .. } => source.payload_strides_vec(),
        }
    }

    fn is_f64(&self) -> bool {
        match self {
            Self::Materialized(storage) => storage.is_f64(),
            Self::Eager { inner, .. } => inner.dtype() == DType::F64,
            Self::Compact(payload) => payload.payload.dtype() == DType::F64,
            Self::Deferred { source, .. } => source.is_f64(),
        }
    }

    fn is_c64(&self) -> bool {
        match self {
            Self::Materialized(storage) => storage.is_c64(),
            Self::Eager { inner, .. } => inner.dtype() == DType::C64,
            Self::Compact(payload) => payload.payload.dtype() == DType::C64,
            Self::Deferred { source, .. } => source.is_c64(),
        }
    }

    fn dtype(&self) -> Option<DType> {
        match self {
            Self::Materialized(storage) => Some(if storage.is_c64() {
                DType::C64
            } else {
                DType::F64
            }),
            Self::Eager { inner, .. } => Some(inner.dtype()),
            Self::Compact(payload) => Some(payload.payload.dtype()),
            Self::Deferred { source, .. } => source.dtype(),
        }
    }

    fn is_complex(&self) -> bool {
        match self {
            Self::Materialized(storage) => storage.is_complex(),
            Self::Eager { inner, .. } => matches!(inner.dtype(), DType::C32 | DType::C64),
            Self::Compact(payload) => {
                matches!(payload.payload.dtype(), DType::C32 | DType::C64)
            }
            Self::Deferred { source, .. } => source.is_complex(),
        }
    }

    fn is_diag(&self) -> bool {
        match self {
            Self::Materialized(storage) => storage.is_diag(),
            Self::Eager { axis_classes, .. } => IdxTensor::is_diag_axis_classes(axis_classes),
            Self::Compact(payload) => IdxTensor::is_diag_axis_classes(&payload.axis_classes),
            Self::Deferred { source, .. } => source.is_diag(),
        }
    }

    fn storage_kind(&self) -> StorageKind {
        match self {
            Self::Materialized(storage) => storage.storage_kind(),
            Self::Eager { axis_classes, .. } => {
                if axis_classes.iter().copied().eq(0..axis_classes.len()) {
                    StorageKind::Dense
                } else if IdxTensor::is_diag_axis_classes(axis_classes) {
                    StorageKind::Diagonal
                } else {
                    StorageKind::Structured
                }
            }
            Self::Compact(payload) => {
                if payload
                    .axis_classes
                    .iter()
                    .copied()
                    .eq(0..payload.axis_classes.len())
                {
                    StorageKind::Dense
                } else if IdxTensor::is_diag_axis_classes(&payload.axis_classes) {
                    StorageKind::Diagonal
                } else {
                    StorageKind::Structured
                }
            }
            Self::Deferred { source, .. } => source.storage_kind(),
        }
    }

    fn materialize_eager_payload(inner: &EagerTensor) -> Result<NativeTensor> {
        let native = inner.duplicate_value()?;
        if native.is_col_major_contiguous()? {
            return Ok(native);
        }
        let read = inner.tensor_read();
        Ok(inner.runtime().with_execution_session(|session| {
            session.to_contiguous_read(TensorRead::from_view(read.tensor_view()))
        })??)
    }

    fn materialize(
        &self,
        logical_rank: usize,
    ) -> std::result::Result<Arc<Storage>, TensorStorageError> {
        match self {
            Self::Materialized(storage) => Ok(Arc::clone(storage)),
            Self::Eager {
                inner,
                axis_classes,
            } => {
                let native = Self::materialize_eager_payload(inner).map_err(|source| {
                    TensorStorageError::Materialization {
                        source: Arc::from(source.into_boxed_dyn_error()),
                    }
                })?;
                let dtype = native.dtype();
                if matches!(dtype, DType::F32 | DType::C32) {
                    return Err(TensorStorageError::UnsupportedDtype {
                        dtype: IdxTensor::dtype_name(dtype),
                    });
                }
                IdxTensor::storage_from_native_with_axis_classes(
                    &native,
                    axis_classes,
                    logical_rank,
                )
                .map(Arc::new)
                .map_err(|source| TensorStorageError::Materialization {
                    source: Arc::from(source.into_boxed_dyn_error()),
                })
            }
            Self::Compact(payload) => {
                let native =
                    Self::materialize_eager_payload(&payload.payload).map_err(|source| {
                        TensorStorageError::Materialization {
                            source: Arc::from(source.into_boxed_dyn_error()),
                        }
                    })?;
                let dtype = native.dtype();
                if matches!(dtype, DType::F32 | DType::C32) {
                    return Err(TensorStorageError::UnsupportedDtype {
                        dtype: IdxTensor::dtype_name(dtype),
                    });
                }
                IdxTensor::storage_from_native_with_axis_classes(
                    &native,
                    &payload.axis_classes,
                    logical_rank,
                )
                .map(Arc::new)
                .map_err(|source| TensorStorageError::Materialization {
                    source: Arc::from(source.into_boxed_dyn_error()),
                })
            }
            Self::Deferred { error, .. } => Err((**error).clone()),
        }
    }

    fn scale_eager_payload(&self, scalar: &AnyScalar) -> Result<Self> {
        let (payload, payload_dims, axis_classes) = match self {
            Self::Materialized(storage) => {
                let native = if storage.is_f64() {
                    let values = storage
                        .payload_f64_col_major_vec()
                        .map_err(anyhow::Error::new)?;
                    dense_native_tensor_from_col_major(&values, storage.payload_dims())?
                } else {
                    let values = storage
                        .payload_c64_col_major_vec()
                        .map_err(anyhow::Error::new)?;
                    dense_native_tensor_from_col_major(&values, storage.payload_dims())?
                };
                (
                    EagerTensor::from_tensor_in(native, default_eager_ctx()?)?,
                    storage.payload_dims().to_vec(),
                    storage.axis_classes().to_vec(),
                )
            }
            Self::Eager {
                inner,
                axis_classes,
            } => (
                (**inner).clone(),
                inner.shape().to_vec(),
                axis_classes.clone(),
            ),
            Self::Compact(compact) => (
                (*compact.payload).clone(),
                compact.payload_dims.clone(),
                compact.axis_classes.clone(),
            ),
            Self::Deferred { error, .. } => return Err(anyhow::Error::new((**error).clone())),
        };
        let scalar_inner = scalar.as_tensor()?.try_materialized_inner()?;
        let target_dtype = IdxTensor::scale_target_dtype(payload.dtype(), scalar_inner.dtype())?;
        let payload = if payload.dtype() == target_dtype {
            payload
        } else {
            payload.cast(target_dtype)?
        };
        let scalar_inner = if scalar_inner.dtype() == target_dtype {
            scalar_inner.clone()
        } else {
            scalar_inner.cast(target_dtype)?
        };
        let scaled = if payload.shape().is_empty() {
            payload.mul(&scalar_inner)?
        } else {
            let subscripts = IdxTensor::scale_subscripts(payload.shape().len())?;
            [&payload, &scalar_inner].einsum_subscripts(&subscripts)?
        };
        match self {
            Self::Eager { .. } => Ok(Self::Eager {
                inner: Arc::new(scaled),
                axis_classes,
            }),
            Self::Compact(_) | Self::Materialized(_) => {
                Ok(Self::Compact(Arc::new(StructuredPayload {
                    payload: Arc::new(scaled),
                    payload_dims,
                    axis_classes,
                })))
            }
            Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
        }
    }

    fn conjugate_with<F>(&self, conjugate: &F) -> std::result::Result<Self, TensorStorageError>
    where
        F: Fn(
            &EagerTensor,
        ) -> std::result::Result<
            EagerTensor,
            Arc<dyn std::error::Error + Send + Sync + 'static>,
        >,
    {
        match self {
            Self::Materialized(storage) => Ok(Self::Materialized(Arc::new(storage.conj()))),
            Self::Eager {
                inner,
                axis_classes,
            } => conjugate(inner)
                .map(|conjugated| Self::Eager {
                    inner: Arc::new(conjugated),
                    axis_classes: axis_classes.clone(),
                })
                .map_err(|source| TensorStorageError::Conjugation { source }),
            Self::Compact(payload) => conjugate(payload.payload.as_ref())
                .map(|conjugated| {
                    Self::Compact(Arc::new(StructuredPayload {
                        payload: Arc::new(conjugated),
                        payload_dims: payload.payload_dims.clone(),
                        axis_classes: payload.axis_classes.clone(),
                    }))
                })
                .map_err(|source| TensorStorageError::Conjugation { source }),
            Self::Deferred { error, .. } => Err((**error).clone()),
        }
    }

    fn sum_scalar(&self) -> Result<AnyScalar> {
        match self {
            Self::Materialized(storage) => {
                if storage.is_f64() {
                    Ok(AnyScalar::new_real(storage.sum::<f64>()))
                } else {
                    let value = storage.sum::<Complex64>();
                    Ok(AnyScalar::new_complex(value.re, value.im))
                }
            }
            Self::Eager { inner, .. } => IdxTensor::native_sum_scalar(inner),
            Self::Compact(payload) => IdxTensor::native_sum_scalar(&payload.payload),
            Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
        }
    }

    fn nonfinite_flags(&self) -> Result<(bool, bool)> {
        match self {
            Self::Materialized(storage) => Ok(storage.payload_nonfinite_flags()),
            Self::Eager { inner, .. } => IdxTensor::native_nonfinite_flags(inner),
            Self::Compact(payload) => IdxTensor::native_nonfinite_flags(&payload.payload),
            Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
        }
    }

    fn payload_value_at(&self, payload_coords: &[usize]) -> Result<Complex64> {
        match self {
            Self::Materialized(storage) => storage
                .scalar_at(payload_coords)
                .map(Complex64::from)
                .map_err(anyhow::Error::new),
            Self::Eager { inner, .. } => {
                IdxTensor::native_complex_payload_value_at(inner, payload_coords)
            }
            Self::Compact(payload) => {
                IdxTensor::native_complex_payload_value_at(&payload.payload, payload_coords)
            }
            Self::Deferred { error, .. } => Err(anyhow::Error::new((**error).clone())),
        }
    }

    fn for_each_payload_value(&self, mut f: impl FnMut(Complex64)) -> Result<()> {
        let payload_dims = self.payload_dims();
        let payload_len = checked_product(payload_dims)?;
        let mut payload_coords = vec![0usize; payload_dims.len()];
        for _ in 0..payload_len {
            f(self.payload_value_at(&payload_coords)?);
            let mut carry = true;
            for (coordinate, &dim) in payload_coords.iter_mut().zip(payload_dims.iter()) {
                if !carry {
                    break;
                }
                *coordinate += 1;
                if *coordinate == dim {
                    *coordinate = 0;
                } else {
                    carry = false;
                }
            }
        }
        Ok(())
    }

    fn compact_payload(&self) -> Option<&StructuredPayload> {
        match self {
            Self::Compact(payload) => Some(payload.as_ref()),
            Self::Deferred { source, .. } => source.compact_payload(),
            _ => None,
        }
    }
}

/// Errors returned when constructing a compact copy-selector tensor.
/// A copy-selector has logical values
/// `scale * delta(left, right) * delta(site, selected_value)` and is used to
/// carry a bond through a fixed physical site without dense bond-squared storage.
/// # Examples
/// ```
/// use tensor4all_core::{DynIndex, StructuredSelectorError, IdxTensor};
/// let left = DynIndex::new_dyn(2);
/// let site = DynIndex::new_dyn(3);
/// let right = DynIndex::new_dyn(4);
/// let error = IdxTensor::from_copy_selector(left, site, right, 1, 1.0_f64)
///     .unwrap_err();
/// assert!(matches!(error, StructuredSelectorError::BondDimensionMismatch { .. }));
/// ```
#[derive(Debug, thiserror::Error)]
pub enum StructuredSelectorError {
    /// The two logical copy axes have different dimensions.
    #[error("copy-selector bond dimensions differ: left={left}, right={right}")]
    BondDimensionMismatch {
        /// Dimension of the left copy axis.
        left: usize,
        /// Dimension of the right copy axis.
        right: usize,
    },
    /// One of the logical axes has dimension zero.
    #[error("copy-selector {axis} dimension must be positive")]
    ZeroDimension {
        /// Name of the zero-dimensional axis.
        axis: &'static str,
    },
    /// The selected physical coordinate is outside the site dimension.
    #[error("selected site value {value} is outside 0..{site_dim}")]
    SelectedValueOutOfBounds {
        /// Requested zero-based physical coordinate.
        value: usize,
        /// Dimension of the physical site.
        site_dim: usize,
    },
    /// The compact payload element count cannot be represented by `usize`.
    #[error("copy-selector payload size overflows usize for dimensions {bond_dim} x {site_dim}")]
    PayloadSizeOverflow {
        /// Dimension shared by the copy axes.
        bond_dim: usize,
        /// Dimension of the physical site.
        site_dim: usize,
    },
    /// A compact payload stride cannot be represented by `isize`.
    #[error("copy-selector bond stride {bond_dim} exceeds isize::MAX")]
    StrideOverflow {
        /// Bond dimension that could not be converted to a stride.
        bond_dim: usize,
    },
    /// Reserving the compact payload failed.
    #[error("could not allocate copy-selector payload with {elements} elements")]
    AllocationFailed {
        /// Number of compact payload elements requested.
        elements: usize,
    },
    /// Backend structured-storage validation failed.
    #[error("invalid copy-selector storage: {message}")]
    InvalidStorage {
        /// Diagnostic returned by structured-storage validation.
        message: String,
    },
}

/// Dynamic-rank tensor with structured payload storage -- the central data type
/// of tensor4all.
/// `IdxTensor` stores a logical multi-dimensional tensor of supported scalar
/// values (`f32`, `f64`, `Complex32`, or `Complex64`) together with a list of
/// [`DynIndex`] labels. `f64`/`Complex64` tensors may use compact [`Storage`]
/// snapshots; `f32`/`Complex32` tensors retain an eager payload as the
/// authoritative representation because compact storage supports only the
/// 64-bit dtypes. The logical layout may be dense, diagonal, or explicitly
/// structured. The indices carry unique
/// identities (UUIDs) so that contraction, addition, and other binary
/// operations can automatically match legs by identity rather than position.
/// # Key Operations
/// | Operation | Method |
/// |-----------|--------|
/// | Create from data | [`from_dense`](Self::from_dense), [`from_diag`](Self::from_diag), [`zeros`](Self::zeros) |
/// | Extract data | [`to_vec`](Self::to_vec), [`into_dense_col_major_parts`](Self::into_dense_col_major_parts), [`sum`](Self::sum), [`only`](Self::only) |
/// | Contraction | [`contract`](Self::contract) |
/// | Arithmetic | [`add`](Self::add), [`scale`](Self::scale), [`axpby`](Self::axpby) |
/// | Factorization | via [`TensorFactorizationLike::factorize`](crate::TensorFactorizationLike::factorize) |
/// | Norms | [`norm`](Self::norm), [`norm_squared`](Self::norm_squared), [`maxabs`](Self::maxabs) |
/// | Index ops | [`replaceind`](Self::replaceind), [`permute_indices`](Self::permute_indices) |
/// # Data Layout
/// Logical dense extraction uses **column-major** order (first index varies
/// fastest), matching Fortran, Julia, and ITensors.jl conventions. Compact
/// structured payloads additionally carry explicit payload dimensions, strides,
/// and logical-axis classes.
/// # Examples
/// ```
/// use tensor4all_core::{IdxTensor, DynIndex};
/// // Create a 2x3 real tensor
/// let i = DynIndex::new_dyn(2);
/// let j = DynIndex::new_dyn(3);
/// let data = vec![1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0];
/// let t = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
/// assert_eq!(t.dims(), vec![2, 3]);
/// assert!(t.is_f64());
/// // Sum all elements: 1+2+3+4+5+6 = 21
/// let s = t.sum().unwrap();
/// assert!((s.real() - 21.0).abs() < 1e-12);
/// // Extract data back out
/// let data_out = t.to_vec::<f64>().unwrap();
/// assert_eq!(data_out, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
/// ```
#[derive(Clone)]
pub struct IdxTensor {
    /// Full index information (includes tags and other metadata).
    pub indices: Vec<DynIndex>,
    /// Authoritative payload representation. Compact storage is used when the
    /// dtype is supported by [`Storage`]; otherwise this retains an eager
    /// payload without promotion.
    pub(crate) storage: IdxTensorStorage,
    /// Lazily materialized logical-dense eager payload for native execution and AD.
    pub(crate) eager_cache: Arc<OnceLock<Arc<EagerTensor>>>,
}

impl IdxTensor {
    fn dense_axis_classes(rank: usize) -> Vec<usize> {
        (0..rank).collect()
    }

    fn dtype_name(dtype: DType) -> &'static str {
        match dtype {
            DType::F32 => "f32",
            DType::F64 => "f64",
            DType::C32 => "c32",
            DType::C64 => "c64",
            DType::I32 => "i32",
            DType::I64 => "i64",
            DType::Bool => "bool",
        }
    }

    fn scalar_dtype(&self) -> Result<DType> {
        if let Some(inner) = self.storage.eager() {
            return Ok(inner.dtype());
        }
        if self.storage.is_f64() {
            Ok(DType::F64)
        } else if self.storage.is_c64() {
            Ok(DType::C64)
        } else {
            Err(anyhow::anyhow!(
                "unable to determine IdxTensor scalar dtype"
            ))
        }
    }

    fn diag_axis_classes(rank: usize) -> Vec<usize> {
        if rank == 0 {
            vec![]
        } else {
            vec![0; rank]
        }
    }

    fn canonicalize_axis_classes(axis_classes: &[usize]) -> Vec<usize> {
        let mut map = std::collections::HashMap::new();
        let mut next = 0usize;
        axis_classes
            .iter()
            .map(|&class_id| {
                *map.entry(class_id).or_insert_with(|| {
                    let canonical = next;
                    next += 1;
                    canonical
                })
            })
            .collect()
    }

    fn permute_axis_classes(&self, perm: &[usize]) -> Vec<usize> {
        let axis_classes = self.storage.axis_classes();
        let permuted: Vec<usize> = perm.iter().map(|&index| axis_classes[index]).collect();
        Self::canonicalize_axis_classes(&permuted)
    }

    fn normalize_insert_axis(op: &str, axis: isize, rank: usize) -> Result<usize> {
        let normalized = if axis < 0 {
            rank as isize + 1 + axis
        } else {
            axis
        };
        if !(normalized >= 0 && normalized <= rank as isize) {
            return Err(anyhow::anyhow!(
                "{op}: axis {axis} is out of bounds for inserting into rank {rank}"
            ));
        };
        Ok(normalized as usize)
    }

    fn is_diag_axis_classes(axis_classes: &[usize]) -> bool {
        axis_classes.len() >= 2 && axis_classes.iter().all(|&class_id| class_id == 0)
    }

    fn validate_axis_classes(axis_classes: &[usize], rank: usize) -> Result<()> {
        if axis_classes.len() != rank {
            return Err(anyhow::anyhow!(
                "axis-class rank {} does not match tensor rank {rank}",
                axis_classes.len()
            ));
        }
        if Self::canonicalize_axis_classes(axis_classes) != axis_classes {
            return Err(anyhow::anyhow!(
                "axis classes must be canonical first-occurrence labels: {axis_classes:?}"
            ));
        }
        Ok(())
    }

    fn einsum_subscripts_from_usize_ids(
        inputs: &[Vec<usize>],
        output: &[usize],
    ) -> Result<EinsumSubscripts> {
        let input_labels = inputs
            .iter()
            .map(|ids| {
                ids.iter()
                    .map(|&id| {
                        u32::try_from(id)
                            .map_err(|_| anyhow::anyhow!("einsum label {id} exceeds u32 range"))
                    })
                    .collect::<Result<Vec<_>>>()
            })
            .collect::<Result<Vec<_>>>()?;
        let output_labels = output
            .iter()
            .map(|&id| {
                u32::try_from(id)
                    .map_err(|_| anyhow::anyhow!("einsum label {id} exceeds u32 range"))
            })
            .collect::<Result<Vec<_>>>()?;
        let input_refs = input_labels.iter().map(Vec::as_slice).collect::<Vec<_>>();
        Ok(EinsumSubscripts::new(&input_refs, &output_labels))
    }

    fn build_binary_einsum_subscripts(
        lhs_rank: usize,
        axes_a: &[usize],
        rhs_rank: usize,
        axes_b: &[usize],
    ) -> Result<EinsumSubscripts> {
        if !(axes_a.len() == axes_b.len()) {
            return Err(anyhow::anyhow!(
                "contract axis length mismatch: lhs {:?}, rhs {:?}",
                axes_a,
                axes_b
            ));
        };

        let mut lhs_ids = vec![usize::MAX; lhs_rank];
        let mut rhs_ids = vec![usize::MAX; rhs_rank];
        let mut next_id = 0usize;

        let mut seen_lhs = vec![false; lhs_rank];
        let mut seen_rhs = vec![false; rhs_rank];

        for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
            if !(lhs_axis < lhs_rank) {
                return Err(anyhow::anyhow!("lhs contract axis {lhs_axis} out of range"));
            };
            if !(rhs_axis < rhs_rank) {
                return Err(anyhow::anyhow!("rhs contract axis {rhs_axis} out of range"));
            };
            if !(!seen_lhs[lhs_axis]) {
                return Err(anyhow::anyhow!("duplicate lhs contract axis {lhs_axis}"));
            };
            if !(!seen_rhs[rhs_axis]) {
                return Err(anyhow::anyhow!("duplicate rhs contract axis {rhs_axis}"));
            };
            seen_lhs[lhs_axis] = true;
            seen_rhs[rhs_axis] = true;
            lhs_ids[lhs_axis] = next_id;
            rhs_ids[rhs_axis] = next_id;
            next_id += 1;
        }

        let mut output_ids = Vec::with_capacity(lhs_rank + rhs_rank - 2 * axes_a.len());
        for id in &mut lhs_ids {
            if *id == usize::MAX {
                *id = next_id;
                output_ids.push(next_id);
                next_id += 1;
            }
        }
        for id in &mut rhs_ids {
            if *id == usize::MAX {
                *id = next_id;
                output_ids.push(next_id);
                next_id += 1;
            }
        }

        Self::einsum_subscripts_from_usize_ids(&[lhs_ids, rhs_ids], &output_ids)
    }

    fn binary_dot_general_config(axes_a: &[usize], axes_b: &[usize]) -> Result<DotGeneralConfig> {
        if !(axes_a.len() == axes_b.len()) {
            return Err(anyhow::anyhow!(
                "contract axis length mismatch: lhs {:?}, rhs {:?}",
                axes_a,
                axes_b
            ));
        };
        Ok(DotGeneralConfig {
            lhs_contracting_dims: axes_a.to_vec(),
            rhs_contracting_dims: axes_b.to_vec(),
            lhs_batch_dims: vec![],
            rhs_batch_dims: vec![],
        })
    }

    fn binary_contraction_axis_classes(
        lhs_axis_classes: &[usize],
        axes_a: &[usize],
        rhs_axis_classes: &[usize],
        axes_b: &[usize],
    ) -> Result<Vec<usize>> {
        debug_assert_eq!(axes_a.len(), axes_b.len());

        fn find(parent: &mut [usize], value: usize) -> usize {
            if parent[value] != value {
                parent[value] = find(parent, parent[value]);
            }
            parent[value]
        }

        fn union(parent: &mut [usize], lhs: usize, rhs: usize) {
            let lhs_root = find(parent, lhs);
            let rhs_root = find(parent, rhs);
            if lhs_root != rhs_root {
                parent[rhs_root] = lhs_root;
            }
        }

        let lhs_payload_rank = match lhs_axis_classes.iter().copied().max() {
            Some(value) => value
                .checked_add(1)
                .ok_or_else(|| anyhow::anyhow!("left payload rank overflows usize"))?,
            None => 0,
        };
        let rhs_payload_rank = match rhs_axis_classes.iter().copied().max() {
            Some(value) => value
                .checked_add(1)
                .ok_or_else(|| anyhow::anyhow!("right payload rank overflows usize"))?,
            None => 0,
        };
        let rhs_offset = lhs_payload_rank;
        let parent_len = lhs_payload_rank
            .checked_add(rhs_payload_rank)
            .ok_or_else(|| anyhow::anyhow!("payload rank sum overflows usize"))?;
        let mut parent: Vec<usize> = (0..parent_len).collect();

        for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
            union(
                &mut parent,
                lhs_axis_classes[lhs_axis],
                rhs_offset
                    .checked_add(rhs_axis_classes[rhs_axis])
                    .ok_or_else(|| anyhow::anyhow!("rhs axis-class offset overflows usize"))?,
            );
        }

        let mut lhs_contracted = vec![false; lhs_axis_classes.len()];
        for &axis in axes_a {
            lhs_contracted[axis] = true;
        }
        let mut rhs_contracted = vec![false; rhs_axis_classes.len()];
        for &axis in axes_b {
            rhs_contracted[axis] = true;
        }

        let mut root_to_class = std::collections::HashMap::new();
        let mut next_class = 0usize;
        let mut axis_classes = Vec::with_capacity(
            lhs_axis_classes
                .len()
                .saturating_add(rhs_axis_classes.len()),
        );

        for (axis, &class_id) in lhs_axis_classes.iter().enumerate() {
            if !lhs_contracted[axis] {
                let root = find(&mut parent, class_id);
                let class = *root_to_class.entry(root).or_insert_with(|| {
                    let value = next_class;
                    next_class += 1;
                    value
                });
                axis_classes.push(class);
            }
        }
        for (axis, &class_id) in rhs_axis_classes.iter().enumerate() {
            if !rhs_contracted[axis] {
                let rhs_class = rhs_offset
                    .checked_add(class_id)
                    .ok_or_else(|| anyhow::anyhow!("rhs axis-class offset overflows usize"))?;
                let root = find(&mut parent, rhs_class);
                let class = *root_to_class.entry(root).or_insert_with(|| {
                    let value = next_class;
                    next_class += 1;
                    value
                });
                axis_classes.push(class);
            }
        }

        Ok(axis_classes)
    }

    fn scale_subscripts(rank: usize) -> Result<EinsumSubscripts> {
        let ids: Vec<usize> = (0..rank).collect();
        Self::einsum_subscripts_from_usize_ids(&[ids.clone(), Vec::new()], &ids)
    }

    fn scale_target_dtype(payload: DType, scalar: DType) -> Result<DType> {
        let target = match payload {
            DType::F32 => match scalar {
                DType::C32 | DType::C64 => DType::C32,
                DType::F32 | DType::F64 => DType::F32,
                dtype => {
                    return Err(anyhow::anyhow!(
                        "unsupported scalar dtype {dtype:?} for f32 scaling"
                    ));
                }
            },
            DType::C32 => match scalar {
                DType::F32 | DType::F64 | DType::C32 | DType::C64 => DType::C32,
                dtype => {
                    return Err(anyhow::anyhow!(
                        "unsupported scalar dtype {dtype:?} for c32 scaling"
                    ));
                }
            },
            DType::F64 => match scalar {
                DType::C32 | DType::C64 => DType::C64,
                DType::F32 | DType::F64 => DType::F64,
                dtype => {
                    return Err(anyhow::anyhow!(
                        "unsupported scalar dtype {dtype:?} for f64 scaling"
                    ));
                }
            },
            DType::C64 => match scalar {
                DType::F32 | DType::F64 | DType::C32 | DType::C64 => DType::C64,
                dtype => {
                    return Err(anyhow::anyhow!(
                        "unsupported scalar dtype {dtype:?} for c64 scaling"
                    ));
                }
            },
            dtype => {
                return Err(anyhow::anyhow!(
                    "unsupported tensor dtype {dtype:?} for scaling"
                ));
            }
        };
        Ok(target)
    }

    fn validate_indices(indices: &[DynIndex]) -> Result<()> {
        let mut seen = HashSet::new();
        for idx in indices {
            if !(seen.insert(idx.clone())) {
                return Err(anyhow::anyhow!("Tensor indices must all be unique"));
            };
        }
        Ok(())
    }

    fn validate_diag_dims(dims: &[usize]) -> Result<()> {
        if !dims.is_empty() {
            let first_dim = dims[0];
            for (i, &dim) in dims.iter().enumerate() {
                if !(dim == first_dim) {
                    return Err(anyhow::anyhow!("DiagTensor requires all indices to have the same dimension, but dims[{i}] = {dim} != dims[0] = {first_dim}"));
                };
            }
        }
        Ok(())
    }

    fn seed_native_payload(storage: &Storage, dims: &[usize]) -> Result<NativeTensor> {
        Ok(storage_to_native_tensor(storage, dims)?)
    }

    fn empty_eager_cache() -> Arc<OnceLock<Arc<EagerTensor>>> {
        Arc::new(OnceLock::new())
    }

    fn eager_cache_with(inner: EagerTensor) -> Arc<OnceLock<Arc<EagerTensor>>> {
        let cache = Arc::new(OnceLock::new());
        let _ = cache.set(Arc::new(inner));
        cache
    }

    fn compact_payload_inner(&self) -> Result<EagerTensor> {
        self.ensure_storage_ready()?;
        if let Some(inner) = self.storage.eager() {
            return Ok(inner.clone());
        }
        Ok(EagerTensor::from_tensor_in(
            storage_payload_native(self.storage.materialize(self.indices.len())?.as_ref())?,
            default_eager_ctx()?,
        )?)
    }

    fn dense_inner_from_payload(
        payload: &EagerTensor,
        axis_classes: &[usize],
        logical_dims: &[usize],
    ) -> Result<EagerTensor> {
        let payload_rank = match axis_classes.iter().copied().max() {
            Some(class_id) => class_id
                .checked_add(1)
                .ok_or_else(|| anyhow::anyhow!("structured payload class rank overflows usize"))?,
            None => 0,
        };
        if !(payload.shape().len() == payload_rank) {
            return Err(anyhow::anyhow!(
                "structured payload rank {} does not match axis classes {:?}",
                payload.shape().len(),
                axis_classes
            ));
        };
        if !(logical_dims.len() == axis_classes.len()) {
            return Err(anyhow::anyhow!(
                "logical rank {} does not match axis class rank {}",
                logical_dims.len(),
                axis_classes.len()
            ));
        };

        if axis_classes == Self::dense_axis_classes(logical_dims.len()) {
            if !(payload.shape() == logical_dims) {
                return Err(anyhow::anyhow!(
                    "dense payload dims {:?} do not match logical dims {:?}",
                    payload.shape(),
                    logical_dims
                ));
            };
            return Ok(payload.clone());
        }

        let mut first_axis_by_class = vec![None; payload_rank];
        let mut dense = payload.clone();
        for (logical_axis, &class_id) in axis_classes.iter().enumerate() {
            let first_axis = match first_axis_by_class[class_id] {
                Some(first_axis) => first_axis,
                None => {
                    first_axis_by_class[class_id] = Some(logical_axis);
                    continue;
                }
            };
            dense = dense.embed_diag(first_axis, logical_axis)?;
        }
        if !(dense.shape() == logical_dims) {
            return Err(anyhow::anyhow!(
                "expanded structured payload dims {:?} do not match logical dims {:?}",
                dense.shape(),
                logical_dims
            ));
        };
        Ok(dense)
    }

    fn tracked_compact_payload_value(&self) -> Option<&StructuredPayload> {
        self.storage
            .deferred_error()
            .is_none()
            .then_some(self.storage.compact_payload())
            .flatten()
            .filter(|value| value.payload.tracks_grad())
    }

    fn ensure_storage_ready(&self) -> Result<()> {
        if let Some(error) = self.storage.deferred_error() {
            return Err(anyhow::Error::new(error.clone()));
        }
        Ok(())
    }

    fn compact_payload_is_logical_dense(&self, payload_dims: &[usize]) -> bool {
        self.storage.axis_classes() == Self::dense_axis_classes(self.indices.len())
            && payload_dims == self.dims()
    }

    fn uses_tracked_compact_storage(&self) -> bool {
        self.tracked_compact_payload_value()
            .is_some_and(|value| !self.compact_payload_is_logical_dense(&value.payload_dims))
    }

    fn ensure_shape_packing_preserves_ad(&self, op_name: &str) -> Result<()> {
        self.ensure_storage_ready()?;
        if !(!self.uses_tracked_compact_storage()) {
            return Err(anyhow::anyhow!("{op_name}: structured AD tensors with compact storage are not supported because materializing compact storage would detach gradients"));
        };
        Ok(())
    }

    fn operand_indices_for_contraction(&self, conjugate: bool) -> Vec<DynIndex> {
        if conjugate {
            self.indices.iter().map(|index| index.conj()).collect()
        } else {
            self.indices.clone()
        }
    }

    fn build_binary_contraction_labels(
        lhs_rank: usize,
        axes_a: &[usize],
        rhs_rank: usize,
        axes_b: &[usize],
    ) -> Result<(Vec<usize>, Vec<usize>, Vec<usize>)> {
        if !(axes_a.len() == axes_b.len()) {
            return Err(anyhow::anyhow!(
                "contract axis length mismatch: lhs {:?}, rhs {:?}",
                axes_a,
                axes_b
            ));
        };

        let mut lhs_ids = vec![usize::MAX; lhs_rank];
        let mut rhs_ids = vec![usize::MAX; rhs_rank];
        let mut next_id = 0usize;

        let mut seen_lhs = vec![false; lhs_rank];
        let mut seen_rhs = vec![false; rhs_rank];

        for (&lhs_axis, &rhs_axis) in axes_a.iter().zip(axes_b.iter()) {
            if !(lhs_axis < lhs_rank) {
                return Err(anyhow::anyhow!("lhs contract axis {lhs_axis} out of range"));
            };
            if !(rhs_axis < rhs_rank) {
                return Err(anyhow::anyhow!("rhs contract axis {rhs_axis} out of range"));
            };
            if !(!seen_lhs[lhs_axis]) {
                return Err(anyhow::anyhow!("duplicate lhs contract axis {lhs_axis}"));
            };
            if !(!seen_rhs[rhs_axis]) {
                return Err(anyhow::anyhow!("duplicate rhs contract axis {rhs_axis}"));
            };
            seen_lhs[lhs_axis] = true;
            seen_rhs[rhs_axis] = true;
            lhs_ids[lhs_axis] = next_id;
            rhs_ids[rhs_axis] = next_id;
            next_id += 1;
        }

        let mut output_ids = Vec::with_capacity(lhs_rank + rhs_rank - 2 * axes_a.len());
        for id in &mut lhs_ids {
            if *id == usize::MAX {
                *id = next_id;
                output_ids.push(next_id);
                next_id += 1;
            }
        }
        for id in &mut rhs_ids {
            if *id == usize::MAX {
                *id = next_id;
                output_ids.push(next_id);
                next_id += 1;
            }
        }

        Ok((lhs_ids, rhs_ids, output_ids))
    }

    fn build_payload_einsum_subscripts(
        input_roots: &[Vec<usize>],
        output_roots: &[usize],
    ) -> Result<EinsumSubscripts> {
        Self::einsum_subscripts_from_usize_ids(input_roots, output_roots)
    }

    fn normalize_eager_payload_for_roots(
        payload: &EagerTensor,
        roots: &[usize],
    ) -> Result<(Option<EagerTensor>, Vec<usize>)> {
        if !(payload.shape().len() == roots.len()) {
            return Err(anyhow::anyhow!(
                "payload rank {} does not match root label count {}",
                payload.shape().len(),
                roots.len()
            ));
        };

        let mut current_payload = None;
        let mut current_roots = roots.to_vec();
        while let Some((axis_a, axis_b)) = Self::first_duplicate_pair(&current_roots) {
            let source = current_payload.as_ref().unwrap_or(payload);
            current_payload = Some(source.extract_diag(axis_a, axis_b)?);
            current_roots.remove(axis_b);
        }

        Ok((current_payload, current_roots))
    }

    fn first_duplicate_pair(values: &[usize]) -> Option<(usize, usize)> {
        let mut first_axis_by_value = std::collections::HashMap::new();
        for (axis, &value) in values.iter().enumerate() {
            if let Some(&first_axis) = first_axis_by_value.get(&value) {
                return Some((first_axis, axis));
            }
            first_axis_by_value.insert(value, axis);
        }
        None
    }

    fn from_structured_payload_inner(
        indices: Vec<DynIndex>,
        payload_inner: EagerTensor,
        payload_dims: Vec<usize>,
        axis_classes: Vec<usize>,
    ) -> Result<Self> {
        Self::validate_indices(&indices)?;
        if payload_inner.shape() != payload_dims {
            return Err(anyhow::anyhow!(
                "structured payload dims {:?} do not match planned payload dims {:?}",
                payload_inner.shape(),
                payload_dims
            ));
        }
        if axis_classes == Self::dense_axis_classes(indices.len()) {
            return Self::from_inner_with_axis_classes(indices, payload_inner, axis_classes);
        }
        let structured_payload = Arc::new(StructuredPayload {
            payload: Arc::new(payload_inner),
            payload_dims,
            axis_classes,
        });
        Ok(Self {
            indices,
            storage: IdxTensorStorage::Compact(structured_payload),
            eager_cache: Self::empty_eager_cache(),
        })
    }

    fn contract_structured_payloads(
        &self,
        other: &Self,
        result_indices: Vec<DynIndex>,
        axes_a: &[usize],
        axes_b: &[usize],
    ) -> Result<Self> {
        let (lhs_labels, rhs_labels, output_labels) = Self::build_binary_contraction_labels(
            self.indices.len(),
            axes_a,
            other.indices.len(),
            axes_b,
        )?;
        Self::contract_structured_payloads_nary(
            &[self, other],
            result_indices,
            vec![lhs_labels, rhs_labels],
            output_labels,
        )
    }

    pub(crate) fn contract_structured_payloads_nary(
        operands: &[&Self],
        result_indices: Vec<DynIndex>,
        input_labels: Vec<Vec<usize>>,
        output_labels: Vec<usize>,
    ) -> Result<Self> {
        if !(!operands.is_empty()) {
            return Err(anyhow::anyhow!("structured contraction needs operands"));
        };
        for operand in operands {
            operand.ensure_storage_ready()?;
        }
        let layouts = operands
            .iter()
            .map(|operand| {
                OperandLayout::new(operand.dims(), operand.storage.axis_classes().to_vec())
            })
            .collect::<Result<Vec<_>>>()?;
        let spec = StructuredContractionSpec {
            input_labels,
            output_labels,
            retained_labels: Default::default(),
        };
        let plan = StructuredContractionPlan::new(&layouts, &spec)?;
        let any_grad = operands.iter().any(|operand| operand.tracks_grad());

        if any_grad {
            let dtypes = operands
                .iter()
                .map(|operand| {
                    operand.storage.dtype().ok_or_else(|| {
                        anyhow::anyhow!("structured contraction operand has no scalar dtype")
                    })
                })
                .collect::<Result<Vec<_>>>()?;
            let target = Self::common_eager_dtype(&dtypes)?;
            let mut payloads = Vec::with_capacity(operands.len());
            for operand in operands {
                let payload = operand.compact_payload_inner()?;
                payloads.push(if payload.dtype() == target {
                    payload
                } else {
                    payload.cast(target)?
                });
            }

            let mut normalized = Vec::with_capacity(payloads.len());
            let mut labels = Vec::with_capacity(payloads.len());
            for (operand_idx, (payload, operand_plan)) in
                payloads.iter().zip(plan.operand_plans.iter()).enumerate()
            {
                let (payload, roots) =
                    Self::normalize_eager_payload_for_roots(payload, &operand_plan.class_roots)?;
                normalized.push(payload.unwrap_or_else(|| payloads[operand_idx].clone()));
                labels.push(roots);
            }
            let refs = normalized.iter().collect::<Vec<_>>();
            let subscripts =
                Self::build_payload_einsum_subscripts(&labels, &plan.output_payload_roots)?;
            let payload = refs.as_slice().einsum_subscripts(&subscripts)?;
            return Self::from_structured_payload_inner(
                result_indices,
                payload,
                plan.output_payload_dims,
                plan.output_axis_classes,
            );
        }

        // The native backend borrows contiguous compact payloads and promotes
        // only operands whose compact dtype differs. No logical dense tensor is
        // constructed on this path.
        let storage_owners = operands
            .iter()
            .map(|operand| match &operand.storage {
                IdxTensorStorage::Materialized(storage) => Some(Arc::clone(storage)),
                IdxTensorStorage::Deferred { source, .. } => match source.as_ref() {
                    IdxTensorStorage::Materialized(storage) => Some(Arc::clone(storage)),
                    _ => None,
                },
                _ => None,
            })
            .collect::<Vec<_>>();
        let mut inputs = Vec::with_capacity(operands.len());
        for (operand_idx, operand) in operands.iter().enumerate() {
            if let Some(storage) = storage_owners[operand_idx].as_ref() {
                inputs.push(storage_payload_native_read_input(storage.as_ref())?);
            } else {
                let inner = operand
                    .storage
                    .eager()
                    .ok_or_else(|| anyhow::anyhow!("structured operand has no compact payload"))?;
                inputs.push(NativeTensorReadInput::Borrowed(inner.tensor_read()));
            }
        }
        let mut normalized = Vec::with_capacity(inputs.len());
        let mut labels = Vec::with_capacity(inputs.len());
        for (input, operand_plan) in inputs.into_iter().zip(plan.operand_plans.iter()) {
            let (input, roots) =
                normalize_payload_read_for_roots(input, &operand_plan.class_roots)?;
            normalized.push(input);
            labels.push(roots);
        }
        let refs = normalized
            .iter()
            .zip(labels.iter())
            .map(|(input, labels)| (input, labels.as_slice()))
            .collect::<Vec<_>>();
        let payload = tensor4all_tensorbackend::einsum_native_tensor_reads(
            &refs,
            &plan.output_payload_roots,
        )?;
        // Wrap the result in the operands' common eager runtime when they share
        // one, so explicit-context tensors (e.g. SVD factors) stay in their
        // context instead of falling back to the process-global default.
        // Mixed runtimes keep the legacy default wrap; rejecting them is #623
        // scope, not this seam's.
        let mut common: Option<Arc<EagerRuntime>> = None;
        let mut mixed = false;
        for operand in operands {
            if let Some(inner) = operand.storage.eager() {
                match &common {
                    None => common = Some(Arc::clone(inner.runtime())),
                    Some(runtime) => {
                        if runtime.id() != inner.ctx_id() {
                            mixed = true;
                        }
                    }
                }
            }
        }
        let payload_inner = match (common, mixed) {
            (Some(runtime), false) => EagerTensor::from_tensor_in(payload, runtime)?,
            _ => EagerTensor::from_tensor_in(payload, default_eager_ctx()?)?,
        };
        Self::from_structured_payload_inner(
            result_indices,
            payload_inner,
            plan.output_payload_dims,
            plan.output_axis_classes,
        )
    }

    fn common_eager_dtype(dtypes: &[DType]) -> Result<DType> {
        let target = if dtypes.contains(&DType::C64)
            || (dtypes.contains(&DType::C32) && dtypes.contains(&DType::F64))
        {
            DType::C64
        } else if dtypes.contains(&DType::F64) || dtypes.contains(&DType::C32) {
            if dtypes.contains(&DType::C32) {
                DType::C32
            } else {
                DType::F64
            }
        } else {
            DType::F32
        };
        if !(dtypes
            .iter()
            .all(|dtype| matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64)))
        {
            return Err(anyhow::anyhow!(
                "structured contraction supports only f32, f64, c32, and c64 operands"
            ));
        };
        Ok(target)
    }

    fn should_use_structured_payload_contract(&self, other: &Self) -> bool {
        self.tracks_grad()
            || other.tracks_grad()
            || self.storage.axis_classes() != Self::dense_axis_classes(self.indices.len())
            || other.storage.axis_classes() != Self::dense_axis_classes(other.indices.len())
    }

    fn storage_from_native_with_axis_classes(
        native: &NativeTensor,
        axis_classes: &[usize],
        logical_rank: usize,
    ) -> Result<Storage> {
        if matches!(native.dtype(), DType::F32 | DType::C32) {
            return Err(anyhow::anyhow!(
                "compact IdxTensor storage does not support dtype {:?}; retain the eager payload",
                native.dtype()
            ));
        }
        if Self::is_diag_axis_classes(axis_classes) {
            match native.dtype() {
                DType::F64 | DType::I32 | DType::I64 | DType::Bool => Storage::from_diag_col_major(
                    native_tensor_primal_to_diag::<f64>(native)?,
                    logical_rank,
                ),
                DType::C64 => Storage::from_diag_col_major(
                    native_tensor_primal_to_diag::<Complex64>(native)?,
                    logical_rank,
                ),
                DType::F32 | DType::C32 => Err(anyhow::anyhow!(
                    "compact IdxTensor storage does not support dtype {:?}",
                    native.dtype()
                )),
            }
        } else {
            storage_from_payload_native(native.duplicate()?, native.shape(), axis_classes.to_vec())
        }
    }

    fn dense_selected_diag_payload<T: TensorElement + Copy + Zero>(
        payload: Vec<T>,
        kept_dims: &[usize],
        selected_positions: &[usize],
    ) -> Result<Vec<T>> {
        let output_len = checked_product(kept_dims)?;
        let mut data = vec![T::zero(); output_len];
        if output_len == 0 {
            return Ok(data);
        }

        let Some((&first_position, rest)) = selected_positions.split_first() else {
            return Ok(data);
        };
        if rest.iter().any(|&position| position != first_position) {
            return Ok(data);
        }

        let value = payload[first_position];
        if kept_dims.is_empty() {
            data[0] = value;
            return Ok(data);
        }

        let mut offset = 0usize;
        let mut stride = 1usize;
        for &dim in kept_dims {
            let term = first_position
                .checked_mul(stride)
                .ok_or_else(|| anyhow::anyhow!("diagonal selection offset overflow"))?;
            offset = offset
                .checked_add(term)
                .ok_or_else(|| anyhow::anyhow!("diagonal selection offset overflow"))?;
            stride = stride
                .checked_mul(dim)
                .ok_or_else(|| anyhow::anyhow!("diagonal selection stride overflow"))?;
        }
        data[offset] = value;
        Ok(data)
    }

    fn select_diag_indices(
        &self,
        kept_indices: Vec<DynIndex>,
        kept_dims: Vec<usize>,
        positions: &[usize],
    ) -> Result<Self> {
        if self.storage.is_f64() {
            let storage = self.storage.materialize(self.indices.len())?;
            let payload = storage
                .payload_f64_col_major_vec()
                .map_err(anyhow::Error::new)?;
            let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
            Self::from_dense(kept_indices, data).map_err(anyhow::Error::from)
        } else if self.storage.is_c64() {
            let storage = self.storage.materialize(self.indices.len())?;
            let payload = storage
                .payload_c64_col_major_vec()
                .map_err(anyhow::Error::new)?;
            let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
            Self::from_dense(kept_indices, data).map_err(anyhow::Error::from)
        } else if self.storage.dtype() == Some(DType::F32) {
            let inner = self
                .storage
                .eager()
                .ok_or_else(|| anyhow::anyhow!("failed to read f32 diagonal payload"))?;
            let payload = inner.value()?.as_slice::<f32>()?.to_vec();
            let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
            Self::from_dense(kept_indices, data).map_err(anyhow::Error::from)
        } else if self.storage.dtype() == Some(DType::C32) {
            let inner = self
                .storage
                .eager()
                .ok_or_else(|| anyhow::anyhow!("failed to read c32 diagonal payload"))?;
            let payload = inner.value()?.as_slice::<Complex32>()?.to_vec();
            let data = Self::dense_selected_diag_payload(payload, &kept_dims, positions)?;
            Self::from_dense(kept_indices, data).map_err(anyhow::Error::from)
        } else {
            Err(anyhow::anyhow!("unsupported diagonal storage scalar type"))
        }
    }

    fn col_major_strides(dims: &[usize]) -> Result<Vec<isize>> {
        let mut strides = Vec::with_capacity(dims.len());
        let mut stride = 1isize;
        for &dim in dims {
            strides.push(stride);
            let dim = isize::try_from(dim)
                .map_err(|_| anyhow::anyhow!("dimension does not fit in isize"))?;
            stride = stride
                .checked_mul(dim)
                .ok_or_else(|| anyhow::anyhow!("column-major stride overflow"))?;
        }
        Ok(strides)
    }

    fn zero_structured_selection<T>(
        kept_indices: Vec<DynIndex>,
        kept_dims: &[usize],
    ) -> Result<Self>
    where
        T: TensorElement + Zero,
    {
        let output_len = checked_product(kept_dims)?;
        Self::from_dense(kept_indices, vec![T::zero(); output_len]).map_err(anyhow::Error::from)
    }

    fn selected_structured_class_positions(
        axis_classes: &[usize],
        payload_rank: usize,
        selected_axes: &[usize],
        positions: &[usize],
    ) -> Option<Vec<Option<usize>>> {
        let mut selected_class_positions = vec![None; payload_rank];
        for (&axis, &position) in selected_axes.iter().zip(positions.iter()) {
            let class_id = axis_classes[axis];
            if let Some(existing) = selected_class_positions[class_id] {
                if existing != position {
                    return None;
                }
            } else {
                selected_class_positions[class_id] = Some(position);
            }
        }
        Some(selected_class_positions)
    }

    fn select_structured_indices_typed<T, F>(
        &self,
        payload: Vec<T>,
        kept_axes: &[usize],
        kept_indices: Vec<DynIndex>,
        kept_dims: Vec<usize>,
        selected: (&[usize], &[usize]),
        make_output: F,
    ) -> Result<Self>
    where
        T: TensorElement + Zero,
        F: FnOnce(Vec<T>, Vec<usize>, Vec<isize>, Vec<usize>) -> Result<Self>,
    {
        let (selected_axes, positions) = selected;
        let payload_dims = self.storage.payload_dims();
        let axis_classes = self.storage.axis_classes();
        let payload_rank = payload_dims.len();
        let Some(selected_class_positions) = Self::selected_structured_class_positions(
            axis_classes,
            payload_rank,
            selected_axes,
            positions,
        ) else {
            return Self::zero_structured_selection::<T>(kept_indices, &kept_dims);
        };

        let selected_class_kept = kept_axes
            .iter()
            .any(|&axis| selected_class_positions[axis_classes[axis]].is_some());
        if selected_class_kept {
            return self.select_structured_indices_dense(
                payload,
                kept_axes,
                kept_indices,
                kept_dims,
                &selected_class_positions,
            );
        }

        let mut old_to_new_class = vec![None; payload_rank];
        let mut output_payload_dims = Vec::new();
        let mut output_axis_classes = Vec::with_capacity(kept_axes.len());
        for &axis in kept_axes {
            let class_id = axis_classes[axis];
            let new_class = match old_to_new_class[class_id] {
                Some(new_class) => new_class,
                None => {
                    let new_class = output_payload_dims.len();
                    old_to_new_class[class_id] = Some(new_class);
                    output_payload_dims.push(payload_dims[class_id]);
                    new_class
                }
            };
            output_axis_classes.push(new_class);
        }

        let output_len = checked_product(&output_payload_dims)?;
        let mut output_payload = Vec::with_capacity(output_len);
        for linear in 0..output_len {
            let output_payload_index = decode_col_major_linear(linear, &output_payload_dims)?;
            let mut input_payload_index = vec![0usize; payload_rank];
            for class_id in 0..payload_rank {
                input_payload_index[class_id] =
                    if let Some(position) = selected_class_positions[class_id] {
                        position
                    } else if let Some(new_class) = old_to_new_class[class_id] {
                        output_payload_index[new_class]
                    } else {
                        return Err(anyhow::anyhow!(
                            "structured payload class {class_id} is neither selected nor kept"
                        ));
                    };
            }
            let input_linear = encode_col_major_linear(&input_payload_index, payload_dims)?;
            output_payload.push(payload[input_linear]);
        }

        let output_strides = Self::col_major_strides(&output_payload_dims)?;
        make_output(
            output_payload,
            output_payload_dims,
            output_strides,
            output_axis_classes,
        )
    }

    fn select_structured_indices_dense<T>(
        &self,
        payload: Vec<T>,
        kept_axes: &[usize],
        kept_indices: Vec<DynIndex>,
        kept_dims: Vec<usize>,
        selected_class_positions: &[Option<usize>],
    ) -> Result<Self>
    where
        T: TensorElement + Zero,
    {
        let payload_dims = self.storage.payload_dims();
        let axis_classes = self.storage.axis_classes();
        let output_len = checked_product(&kept_dims)?;
        let mut output = Vec::with_capacity(output_len);

        for linear in 0..output_len {
            let kept_position = decode_col_major_linear(linear, &kept_dims)?;
            let mut input_payload_index = selected_class_positions.to_vec();
            let mut is_structural_zero = false;

            for (&axis, &position) in kept_axes.iter().zip(kept_position.iter()) {
                let class_id = axis_classes[axis];
                match input_payload_index[class_id] {
                    Some(existing) if existing != position => {
                        is_structural_zero = true;
                        break;
                    }
                    Some(_) => {}
                    None => input_payload_index[class_id] = Some(position),
                }
            }

            if is_structural_zero {
                output.push(T::zero());
                continue;
            }

            let input_payload_index = input_payload_index
                .into_iter()
                .enumerate()
                .map(|(class_id, position)| {
                    position.ok_or_else(|| {
                        anyhow::anyhow!(
                            "structured payload class {class_id} is neither selected nor kept"
                        )
                    })
                })
                .collect::<Result<Vec<_>>>()?;
            let input_linear = encode_col_major_linear(&input_payload_index, payload_dims)?;
            output.push(payload[input_linear]);
        }

        Self::from_dense(kept_indices, output).map_err(anyhow::Error::from)
    }

    fn select_structured_indices(
        &self,
        kept_axes: &[usize],
        kept_indices: Vec<DynIndex>,
        kept_dims: Vec<usize>,
        selected_axes: &[usize],
        positions: &[usize],
    ) -> Result<Self> {
        if self.storage.is_f64() {
            let storage = self.storage.materialize(self.indices.len())?;
            let payload = storage
                .payload_f64_col_major_vec()
                .map_err(anyhow::Error::new)?;
            let output_indices = kept_indices.clone();
            self.select_structured_indices_typed(
                payload,
                kept_axes,
                kept_indices,
                kept_dims,
                (selected_axes, positions),
                move |payload, dims, strides, classes| {
                    let storage = Storage::new_structured(payload, dims, strides, classes)?;
                    Self::from_storage(output_indices, Arc::new(storage))
                        .map_err(anyhow::Error::from)
                },
            )
        } else if self.storage.is_c64() {
            let storage = self.storage.materialize(self.indices.len())?;
            let payload = storage
                .payload_c64_col_major_vec()
                .map_err(anyhow::Error::new)?;
            let output_indices = kept_indices.clone();
            self.select_structured_indices_typed(
                payload,
                kept_axes,
                kept_indices,
                kept_dims,
                (selected_axes, positions),
                move |payload, dims, strides, classes| {
                    let storage = Storage::new_structured(payload, dims, strides, classes)?;
                    Self::from_storage(output_indices, Arc::new(storage))
                        .map_err(anyhow::Error::from)
                },
            )
        } else if self.storage.dtype() == Some(DType::F32) {
            let inner = self
                .storage
                .eager()
                .ok_or_else(|| anyhow::anyhow!("failed to read f32 structured payload"))?;
            let payload = inner.value()?.as_slice::<f32>()?.to_vec();
            let output_indices = kept_indices.clone();
            self.select_structured_indices_typed(
                payload,
                kept_axes,
                kept_indices,
                kept_dims,
                (selected_axes, positions),
                move |payload, dims, _strides, classes| {
                    let native = dense_native_tensor_from_col_major(&payload, &dims)?;
                    let inner = EagerTensor::from_tensor_in(native, default_eager_ctx()?)?;
                    Self::from_structured_payload_inner(output_indices, inner, dims, classes)
                },
            )
        } else if self.storage.dtype() == Some(DType::C32) {
            let inner = self
                .storage
                .eager()
                .ok_or_else(|| anyhow::anyhow!("failed to read c32 structured payload"))?;
            let payload = inner.value()?.as_slice::<Complex32>()?.to_vec();
            let output_indices = kept_indices.clone();
            self.select_structured_indices_typed(
                payload,
                kept_axes,
                kept_indices,
                kept_dims,
                (selected_axes, positions),
                move |payload, dims, _strides, classes| {
                    let native = dense_native_tensor_from_col_major(&payload, &dims)?;
                    let inner = EagerTensor::from_tensor_in(native, default_eager_ctx()?)?;
                    Self::from_structured_payload_inner(output_indices, inner, dims, classes)
                },
            )
        } else {
            Err(anyhow::anyhow!(
                "unsupported structured storage scalar type"
            ))
        }
    }

    fn validate_storage_matches_indices(indices: &[DynIndex], storage: &Storage) -> Result<()> {
        let dims = Self::expected_dims_from_indices(indices);
        let storage_dims = storage.logical_dims();
        if storage_dims != dims {
            return Err(anyhow::anyhow!(
                "storage logical dims {:?} do not match indices dims {:?}",
                storage_dims,
                dims
            ));
        }
        if storage.is_diag() {
            Self::validate_diag_dims(&dims)?;
        }
        Ok(())
    }

    fn try_materialized_inner(&self) -> Result<&EagerTensor> {
        self.ensure_storage_ready()?;
        let logical_dims = self.dims();
        if let Some(value) = self.tracked_compact_payload_value() {
            if self.compact_payload_is_logical_dense(&value.payload_dims) {
                return Ok(value.payload.as_ref());
            }
            if self.eager_cache.get().is_none() {
                let dense = Self::dense_inner_from_payload(
                    value.payload.as_ref(),
                    &value.axis_classes,
                    &logical_dims,
                )?;
                let _ = self.eager_cache.set(Arc::new(dense));
            }
            return self
                .eager_cache
                .get()
                .map(|inner| inner.as_ref())
                .ok_or_else(|| {
                    anyhow::anyhow!("IdxTensor structured AD cache was not initialized")
                });
        }
        if let Some(inner) = self.storage.eager() {
            if self.storage.axis_classes() == Self::dense_axis_classes(self.indices.len()) {
                return Ok(inner);
            }
            if self.eager_cache.get().is_none() {
                let dense = Self::dense_inner_from_payload(
                    inner,
                    self.storage.axis_classes(),
                    &logical_dims,
                )?;
                let _ = self.eager_cache.set(Arc::new(dense));
            }
            return self
                .eager_cache
                .get()
                .map(|inner| inner.as_ref())
                .ok_or_else(|| {
                    anyhow::anyhow!("IdxTensor structured eager cache was not initialized")
                });
        }
        if self.eager_cache.get().is_none() {
            let native = profile_pairwise_contract_section("materialize_storage_to_native", || {
                let storage = self.storage.materialize(self.indices.len())?;
                Self::seed_native_payload(storage.as_ref(), &logical_dims)
            })
            .context("IdxTensor materialization failed")?;
            record_pairwise_contract_profile_bytes(
                "materialize_storage_to_native",
                tensor_profile_bytes(native.dtype(), native.shape()),
            );
            let _ = self.eager_cache.set(Arc::new(EagerTensor::from_tensor_in(
                native,
                default_eager_ctx()?,
            )?));
        }
        self.eager_cache
            .get()
            .map(|inner| inner.as_ref())
            .ok_or_else(|| anyhow::anyhow!("IdxTensor materialization cache was not initialized"))
    }

    pub(crate) fn as_inner(&self) -> Result<&EagerTensor> {
        self.try_materialized_inner()
    }

    /// Compute dims from `indices` order.
    #[inline]
    fn expected_dims_from_indices(indices: &[DynIndex]) -> Vec<usize> {
        indices.iter().map(|idx| idx.dim()).collect()
    }

    /// Get dims in the current `indices` order.
    ///
    /// This is computed on-demand from `indices` (single source of truth).
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(3);
    /// let k = DynIndex::new_dyn(4);
    /// let t = IdxTensor::from_dense(
    ///     vec![i, j, k],
    ///     vec![0.0; 24],
    /// ).unwrap();
    /// assert_eq!(t.dims(), vec![2, 3, 4]);
    /// ```
    pub fn dims(&self) -> Vec<usize> {
        Self::expected_dims_from_indices(&self.indices)
    }

    /// Select fixed coordinates for tensor indices and drop those axes.
    ///
    /// The `selected_indices` slice identifies tensor axes by index identity,
    /// and `positions` gives the zero-based coordinate to take on each
    /// selected axis. Unselected indices are preserved in their original order.
    ///
    /// # Arguments
    ///
    /// * `selected_indices` - Indices to fix and remove from the result. Each
    ///
    ///   index must appear exactly once in this tensor.
    /// * `positions` - Coordinates for `selected_indices`. Each coordinate must
    ///
    ///   be less than the corresponding index dimension.
    ///
    /// # Returns
    ///
    /// A tensor over the unselected indices. Selecting no indices returns a
    /// clone of the original tensor. Selecting all indices returns a rank-0
    /// scalar tensor. Diagonal and structured tensors are sliced from their
    /// compact payload without materializing the original full tensor; the
    /// result keeps structured storage when the remaining logical axes can
    /// still be represented by axis classes.
    ///
    /// # Errors
    /// Returns an error when a selected coordinate is out of range for its index
    /// (an out of bounds failure) or when `selected_indices` and `positions`
    /// differ in length (a length mismatch).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(3);
    /// let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
    /// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
    ///
    /// let selected = tensor.select_indices(&[j], &[1]).unwrap();
    /// assert_eq!(selected.dims(), vec![2]);
    /// assert_eq!(selected.to_vec::<f64>().unwrap(), vec![3.0, 4.0]);
    /// ```
    pub fn select_indices(
        &self,
        selected_indices: &[DynIndex],
        positions: &[usize],
    ) -> std::result::Result<Self, IdxTensorError> {
        if selected_indices.len() != positions.len() {
            return Err(anyhow::anyhow!(
                "selected_indices length {} does not match positions length {}",
                selected_indices.len(),
                positions.len()
            )
            .into());
        }
        if selected_indices.is_empty() {
            return Ok(self.clone());
        }

        let mut selected_axes = Vec::with_capacity(selected_indices.len());
        let mut seen_axes = HashSet::with_capacity(selected_indices.len());
        for (selected, &position) in selected_indices.iter().zip(positions.iter()) {
            let axis = self
                .indices
                .iter()
                .position(|index| index == selected)
                .ok_or_else(|| anyhow::anyhow!("selected index is not present in tensor"))?;
            if !seen_axes.insert(axis) {
                return Err(anyhow::anyhow!("selected index appears more than once").into());
            }
            let dim = self.indices[axis].dim();
            if position >= dim {
                return Err(anyhow::anyhow!(
                    "selected coordinate {position} is out of range for axis {axis} with dim {dim}"
                )
                .into());
            }
            selected_axes.push(axis);
        }

        let kept_axes = self
            .indices
            .iter()
            .enumerate()
            .filter(|(axis, _)| !seen_axes.contains(axis))
            .map(|(axis, _)| axis)
            .collect::<Vec<_>>();
        let kept_indices = kept_axes
            .iter()
            .map(|&axis| self.indices[axis].clone())
            .collect::<Vec<_>>();
        let kept_dims = kept_axes
            .iter()
            .map(|&axis| self.indices[axis].dim())
            .collect::<Vec<_>>();

        if matches!(
            self.storage.storage_kind(),
            StorageKind::Diagonal | StorageKind::Structured
        ) {
            self.ensure_shape_packing_preserves_ad("select_indices")?;
        }
        if self.storage.storage_kind() == StorageKind::Diagonal {
            return self
                .select_diag_indices(kept_indices, kept_dims, positions)
                .map_err(IdxTensorError::from);
        }
        if self.storage.storage_kind() == StorageKind::Structured {
            return self
                .select_structured_indices(
                    &kept_axes,
                    kept_indices,
                    kept_dims,
                    &selected_axes,
                    positions,
                )
                .map_err(IdxTensorError::from);
        }
        if self.storage.storage_kind() != StorageKind::Dense {
            return Err(anyhow::anyhow!(
                "select_indices got unsupported storage kind {:?}",
                self.storage.storage_kind()
            )
            .into());
        }

        // Slice positionally axis by axis in the operand's own runtime.
        // (A starts-tensor + dynamic_slice build would need an explicit
        // host-to-device upload for resident operands, which this
        // context-free helper cannot name.)
        let mut sliced = self.try_materialized_inner()?.clone();
        for (&axis, &position) in selected_axes.iter().zip(positions.iter()) {
            sliced = sliced
                .slice_axis(axis, position..position + 1)
                .map_err(|error| anyhow::anyhow!("select_indices slicing failed: {error}"))?;
        }
        Self::from_inner(kept_indices, sliced.reshape(&kept_dims)?).map_err(IdxTensorError::from)
    }

    /// Stack tensors along a newly inserted index.
    ///
    /// Each input must have exactly the same index order and dimensions. The
    /// `new_index` dimension must match the number of input tensors. The
    /// `axis` argument follows tenferro/PyTorch-style insertion semantics:
    /// `0` inserts before the first existing axis and `-1` appends a trailing
    /// axis. Use `axis = -1` for batched contractions because tenferro uses
    /// trailing batch dimensions as the canonical batched-GEMM layout.
    ///
    /// # Errors
    ///
    /// Returns an error if no tensors are provided, the new index dimension
    /// does not match the number of tensors, an input has a different index
    /// order, `axis` is outside the valid insertion range, or a tracked
    /// structured-AD tensor uses compact storage that would need dense
    /// materialization.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let batch = DynIndex::new_dyn(2);
    /// let a = IdxTensor::from_dense(vec![i.clone()], vec![1.0_f64, 2.0]).unwrap();
    /// let b = IdxTensor::from_dense(vec![i.clone()], vec![3.0_f64, 4.0]).unwrap();
    ///
    /// let stacked = IdxTensor::stack_along_new_index(&[&a, &b], batch.clone(), -1).unwrap();
    ///
    /// assert_eq!(stacked.indices(), &[i, batch]);
    /// assert_eq!(stacked.to_vec::<f64>().unwrap(), vec![1.0, 2.0, 3.0, 4.0]);
    /// ```
    pub fn stack_along_new_index(
        tensors: &[&Self],
        new_index: DynIndex,
        axis: isize,
    ) -> std::result::Result<Self, IdxTensorError> {
        let first = tensors
            .first()
            .copied()
            .ok_or_else(|| anyhow::anyhow!("stack_along_new_index requires at least one tensor"))?;
        if !(new_index.dim() == tensors.len()) {
            return Err(anyhow::anyhow!(
                "stack_along_new_index: new index dim {} does not match tensor count {}",
                new_index.dim(),
                tensors.len()
            )
            .into());
        };

        let base_indices = first.indices.clone();
        for tensor in tensors.iter().copied().skip(1) {
            if !(tensor.indices == base_indices) {
                return Err(anyhow::anyhow!(
                    "stack_along_new_index: input tensors must have identical index order"
                )
                .into());
            };
        }
        for &tensor in tensors {
            tensor.ensure_shape_packing_preserves_ad("stack_along_new_index")?;
        }

        let insert_axis =
            Self::normalize_insert_axis("stack_along_new_index", axis, base_indices.len())?;
        let mut result_indices = base_indices;
        result_indices.insert(insert_axis, new_index);

        let inner_refs = tensors
            .iter()
            .map(|tensor| tensor.try_materialized_inner())
            .collect::<Result<Vec<_>>>()?;
        let stacked = EagerTensor::stack(&inner_refs, axis)?;
        Self::from_inner(result_indices, stacked).map_err(IdxTensorError::from)
    }

    /// Select positions along one index and replace it with a new index.
    ///
    /// This is the retained-axis counterpart to [`Self::select_indices`]:
    /// instead of fixing one coordinate and removing the index, it gathers a
    /// list of positions and keeps the gathered axis under `target_index`.
    /// Repeated positions are allowed; reverse-mode AD accumulates repeated
    /// cotangents through tenferro's scatter-add gather transpose.
    ///
    /// # Errors
    /// Returns an error when a selected position is out of range for the source
    /// index (an out of bounds failure) or when the source and target index
    /// dimensions are incompatible (a shape mismatch).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let source = DynIndex::new_dyn(3);
    /// let target = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(
    ///     vec![source.clone()],
    ///     vec![10.0_f64, 20.0, 30.0],
    /// ).unwrap();
    ///
    /// let selected = tensor.index_select(&source, target.clone(), &[2, 0]).unwrap();
    ///
    /// assert_eq!(selected.indices(), &[target]);
    /// assert_eq!(selected.to_vec::<f64>().unwrap(), vec![30.0, 10.0]);
    /// ```
    pub fn index_select(
        &self,
        source_index: &DynIndex,
        target_index: DynIndex,
        positions: &[usize],
    ) -> std::result::Result<Self, IdxTensorError> {
        if !(target_index.dim() == positions.len()) {
            return Err(anyhow::anyhow!(
                "index_select: target index dim {} does not match position count {}",
                target_index.dim(),
                positions.len()
            )
            .into());
        };
        let axis = self
            .indices
            .iter()
            .position(|index| index == source_index)
            .ok_or_else(|| anyhow::anyhow!("index_select: source index is not present"))?;
        let source_dim = self.indices[axis].dim();
        for &position in positions {
            if !(position < source_dim) {
                return Err(anyhow::anyhow!(
                    "index_select: position {position} is out of range for source dim {source_dim}"
                )
                .into());
            };
        }
        self.ensure_shape_packing_preserves_ad("index_select")?;

        let axis = isize::try_from(axis)
            .map_err(|_| anyhow::anyhow!("index_select: axis does not fit in isize"))?;
        let selected = self
            .try_materialized_inner()?
            .index_select(axis, positions)?;
        let mut result_indices = self.indices.clone();
        result_indices[axis as usize] = target_index;
        Self::from_inner(result_indices, selected).map_err(IdxTensorError::from)
    }

    /// Create a new tensor with dynamic rank.
    ///
    /// # Errors
    /// Returns an error when the storage logical dimension does not match the
    /// index dimension product (a shape mismatch) or when duplicate
    /// indices are provided.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::Storage;
    /// use std::sync::Arc;
    ///
    /// let i = DynIndex::new_dyn(3);
    /// let storage = Arc::new(Storage::new_dense::<f64>(3).unwrap());
    /// let t = IdxTensor::new(vec![i], storage).unwrap();
    /// assert_eq!(t.dims(), vec![3]);
    /// ```
    pub fn new(
        indices: Vec<DynIndex>,
        storage: Arc<Storage>,
    ) -> std::result::Result<Self, IdxTensorError> {
        Self::from_storage(indices, storage)
    }

    /// Create a new tensor with dynamic rank, automatically computing dimensions from indices.
    ///
    /// This is a convenience constructor that extracts dimensions from indices using `IndexLike::dim()`.
    ///
    /// # Errors
    /// Returns an error when the storage logical dimension does not match the
    /// index dimension product (a shape mismatch) or when duplicate
    /// indices are provided.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::Storage;
    /// use std::sync::Arc;
    ///
    /// let i = DynIndex::new_dyn(4);
    /// let storage = Arc::new(Storage::new_dense::<f64>(4).unwrap());
    /// let t = IdxTensor::from_indices(vec![i], storage).unwrap();
    /// assert_eq!(t.dims(), vec![4]);
    /// ```
    pub fn from_indices(
        indices: Vec<DynIndex>,
        storage: Arc<Storage>,
    ) -> std::result::Result<Self, IdxTensorError> {
        Self::new(indices, storage)
    }

    /// Create a tensor from explicit compact storage.
    ///
    /// # Errors
    /// Returns an error when the storage scalar kind is incompatible with the
    /// requested operations (a scalar-kind mismatch) or the storage cannot
    /// represent the given index space (a shape mismatch).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::Storage;
    /// use std::sync::Arc;
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let storage = Arc::new(Storage::new_diag(vec![1.0_f64, 2.0]).unwrap());
    /// let t = IdxTensor::from_storage(vec![i, j], storage).unwrap();
    /// assert_eq!(t.dims(), vec![2, 2]);
    /// ```
    pub fn from_storage(
        indices: Vec<DynIndex>,
        storage: Arc<Storage>,
    ) -> std::result::Result<Self, IdxTensorError> {
        Self::validate_indices(&indices)?;
        Self::validate_storage_matches_indices(&indices, storage.as_ref())?;
        Ok(Self {
            indices,
            storage: IdxTensorStorage::from_storage(storage),
            eager_cache: Self::empty_eager_cache(),
        })
    }

    /// Create a tensor from explicit structured storage.
    ///
    /// This is an alias for [`IdxTensor::from_storage`] with a name that
    /// emphasizes that compact structured metadata is preserved.
    ///
    /// # Errors
    /// Returns an error when the structured storage is invalid (an invalid-storage
    /// failure) or the index space is incompatible (a shape mismatch).
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::{Storage, StorageKind};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let storage = Arc::new(Storage::from_diag_col_major(vec![1.0_f64, 2.0], 2).unwrap());
    /// let tensor = IdxTensor::from_structured_storage(vec![i, j], storage).unwrap();
    /// assert_eq!(tensor.storage().unwrap().storage_kind(), StorageKind::Diagonal);
    /// ```
    pub fn from_structured_storage(
        indices: Vec<DynIndex>,
        storage: Arc<Storage>,
    ) -> std::result::Result<Self, IdxTensorError> {
        Self::from_storage(indices, storage)
    }

    /// Construct a compact copy tensor that selects one physical-site value.
    ///
    /// The returned rank-3 tensor has logical indices `[left, site, right]` and
    /// value `scale` exactly when `left == right` and `site == selected_value`;
    /// every other entry is zero. Its payload has `left.dim * site.dim`
    /// elements rather than `left.dim * site.dim * right.dim` dense elements.
    ///
    /// # Arguments
    ///
    /// - `left`: left copy axis; its dimension must be positive and equal to
    ///
    ///   `right.dim`.
    /// - `site`: physical axis whose selected coordinate remains active.
    /// - `right`: right copy axis paired with `left`.
    /// - `selected_value`: zero-based coordinate in `0..site.dim`.
    /// - `scale`: value stored on the selected copy diagonal.
    ///
    /// # Returns
    ///
    /// A structured tensor with axis classes `[0, 1, 0]`. For `f64` and
    /// `Complex64`, compact storage is retained; `f32` and `Complex32` keep
    /// an eager authoritative payload because compact storage has no 32-bit
    /// scalar representation.
    ///
    /// # Errors
    ///
    /// Returns [`StructuredSelectorError`] when dimensions are zero or
    /// inconsistent, the selected value is out of bounds, checked size or
    /// stride arithmetic overflows, allocation fails, or backend structured
    /// storage validation rejects the metadata.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::StorageKind;
    ///
    /// let left = DynIndex::new_dyn(2);
    /// let site = DynIndex::new_dyn(3);
    /// let right = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_copy_selector(
    ///     left,
    ///     site,
    ///     right,
    ///     1,
    ///     2.5_f64,
    /// ).unwrap();
    ///
    /// assert_eq!(tensor.storage().unwrap().storage_kind(), StorageKind::Structured);
    /// assert_eq!(tensor.storage().unwrap().payload_len(), 6);
    /// assert_eq!(
    ///     tensor.to_vec::<f64>().unwrap(),
    ///     vec![0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.5, 0.0, 0.0],
    /// );
    /// ```
    pub fn from_copy_selector<T>(
        left: DynIndex,
        site: DynIndex,
        right: DynIndex,
        selected_value: usize,
        scale: T,
    ) -> std::result::Result<Self, StructuredSelectorError>
    where
        T: TensorElement + Copy + Zero,
    {
        if left.dim == 0 {
            return Err(StructuredSelectorError::ZeroDimension { axis: "left" });
        }
        if site.dim == 0 {
            return Err(StructuredSelectorError::ZeroDimension { axis: "site" });
        }
        if right.dim == 0 {
            return Err(StructuredSelectorError::ZeroDimension { axis: "right" });
        }
        if left.dim != right.dim {
            return Err(StructuredSelectorError::BondDimensionMismatch {
                left: left.dim,
                right: right.dim,
            });
        }
        if selected_value >= site.dim {
            return Err(StructuredSelectorError::SelectedValueOutOfBounds {
                value: selected_value,
                site_dim: site.dim,
            });
        }

        let payload_len =
            left.dim
                .checked_mul(site.dim)
                .ok_or(StructuredSelectorError::PayloadSizeOverflow {
                    bond_dim: left.dim,
                    site_dim: site.dim,
                })?;
        let _site_stride = isize::try_from(left.dim)
            .map_err(|_| StructuredSelectorError::StrideOverflow { bond_dim: left.dim })?;
        let selected_offset = left.dim.checked_mul(selected_value).ok_or(
            StructuredSelectorError::PayloadSizeOverflow {
                bond_dim: left.dim,
                site_dim: site.dim,
            },
        )?;

        let mut payload = Vec::new();
        payload.try_reserve_exact(payload_len).map_err(|_| {
            StructuredSelectorError::AllocationFailed {
                elements: payload_len,
            }
        })?;
        payload.resize(payload_len, T::zero());
        for bond in 0..left.dim {
            payload[selected_offset + bond] = scale;
        }
        let payload_native = dense_native_tensor_from_col_major(&payload, &[left.dim, site.dim])
            .map_err(|error| StructuredSelectorError::InvalidStorage {
                message: error.to_string(),
            })?;
        let payload_dtype = payload_native.dtype();
        let payload_inner = EagerTensor::from_tensor_in(
            payload_native,
            default_eager_ctx().map_err(|error| StructuredSelectorError::InvalidStorage {
                message: error.to_string(),
            })?,
        )
        .map_err(|error| StructuredSelectorError::InvalidStorage {
            message: error.to_string(),
        })?;
        let payload_dims = vec![left.dim, site.dim];
        let indices = vec![left, site, right];
        if !matches!(
            payload_dtype,
            DType::F32 | DType::F64 | DType::C32 | DType::C64
        ) {
            return Err(StructuredSelectorError::InvalidStorage {
                message: format!("unsupported selector dtype {:?}", payload_dtype),
            });
        }
        Self::from_structured_payload_inner(indices, payload_inner, payload_dims, vec![0, 1, 0])
            .map_err(|error| StructuredSelectorError::InvalidStorage {
                message: error.to_string(),
            })
    }

    /// Create a tensor from a native tenferro payload.
    pub(crate) fn from_native(indices: Vec<DynIndex>, native: NativeTensor) -> Result<Self> {
        let axis_classes = Self::dense_axis_classes(indices.len());
        Self::from_native_with_axis_classes(indices, native, axis_classes)
    }

    pub(crate) fn from_native_with_axis_classes(
        indices: Vec<DynIndex>,
        native: NativeTensor,
        axis_classes: Vec<usize>,
    ) -> Result<Self> {
        Self::from_inner_with_axis_classes(
            indices,
            EagerTensor::from_tensor_in(native, default_eager_ctx()?)?,
            axis_classes,
        )
    }

    /// Adopt a backend-derived value whose result indices and axis classes were
    /// validated by the contraction planner that produced it.
    pub(crate) fn from_untracked_native_with_axis_classes(
        indices: Vec<DynIndex>,
        native: NativeTensor,
        axis_classes: Vec<usize>,
    ) -> Result<Self> {
        Self::from_inner_with_validated_metadata(
            indices,
            adopt_untracked_eager_value(default_eager_ctx()?, TensorValue::from_tensor(native))?,
            axis_classes,
        )
    }

    pub(crate) fn from_inner(indices: Vec<DynIndex>, inner: EagerTensor) -> Result<Self> {
        let axis_classes = Self::dense_axis_classes(indices.len());
        Self::from_inner_with_axis_classes(indices, inner, axis_classes)
    }

    /// Compute the Hermitian eigendecomposition of a rank-2 tensor.
    ///
    /// The tensor must have two square matrix axes. The returned eigenvectors
    /// stay in [`IdxTensor`] form so downstream tensor algebra can preserve
    /// AD metadata where the backend supports it. Eigenvalues are returned as
    /// detached real primal values because truncation and rank selection are
    /// nonsmooth control-flow decisions.
    ///
    /// `hermitian_tol` controls the allowed imaginary part of complex
    /// eigenvalues after the backend solve; use a small non-negative value such
    /// as `1e-12` for numerically Hermitian inputs.
    ///
    /// # Errors
    /// Returns an error when the tensor is not rank-2, when the two indices have
    /// unequal dimensions (a shape or shape mismatch), or when the
    /// eigensolver fails to converge (a non-convergence failure).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{AnyScalar, DynIndex, TensorContractionLike, IdxTensor};
    ///
    /// let row = DynIndex::new_dyn(2);
    /// let col = DynIndex::new_dyn(2);
    /// let matrix = IdxTensor::from_dense(
    ///     vec![row.clone(), col.clone()],
    ///     vec![3.0_f64, 0.0, 0.0, 5.0],
    /// ).unwrap();
    ///
    /// let decomp = matrix.hermitian_eigendecomposition(1.0e-12).unwrap();
    /// let eigenvector = decomp
    ///     .eigenvectors
    ///     .select_indices(&[decomp.eigenvector_index.clone()], &[0])
    ///     .unwrap();
    /// let eigenvector_as_col = eigenvector.replaceind(&row, &col).unwrap();
    /// let applied = IdxTensor::contract(&[&matrix, &eigenvector_as_col]).unwrap();
    /// let expected = eigenvector.scale(AnyScalar::new_real(decomp.eigenvalues[0])).unwrap();
    ///
    /// assert!(applied.isapprox(&expected, 1.0e-12, 0.0).unwrap());
    /// ```
    pub fn hermitian_eigendecomposition(
        &self,
        hermitian_tol: f64,
    ) -> std::result::Result<TensorHermitianEigendecomposition, IdxTensorError> {
        if !(self.indices.len() == 2) {
            return Err(anyhow::anyhow!(
                "IdxTensor::hermitian_eigendecomposition requires a rank-2 tensor, got rank {}",
                self.indices.len()
            )
            .into());
        };
        let dims = self.dims();
        if !(dims[0] == dims[1]) {
            return Err(anyhow::anyhow!(
                "IdxTensor::hermitian_eigendecomposition requires a square matrix, got {}x{}",
                dims[0],
                dims[1]
            )
            .into());
        };
        if !(dims[0] > 0) {
            return Err(anyhow::anyhow!(
                "IdxTensor::hermitian_eigendecomposition requires a non-empty matrix"
            )
            .into());
        };
        if !(hermitian_tol.is_finite() && hermitian_tol >= 0.0) {
            return Err(anyhow::anyhow!(
                "IdxTensor::hermitian_eigendecomposition requires a finite non-negative tolerance"
            )
            .into());
        };

        let input = self.try_materialized_inner()?;
        let (values, vectors) = input
            .eigh()
            .map_err(|source| anyhow::anyhow!("Hermitian eigendecomposition failed: {source}"))?;

        let eigenvalue_index = DynIndex::new_dyn(dims[0]);
        let eigenvector_index = DynIndex::new_dyn(dims[0]);
        let eigenvalue_tensor = Self::from_inner(vec![eigenvalue_index], values)?;
        let eigenvalues = Self::read_real_eigenvalues(&eigenvalue_tensor, hermitian_tol)
            .with_context(|| {
                "IdxTensor::hermitian_eigendecomposition failed to read eigenvalues"
            })?;
        let eigenvectors = Self::from_inner(
            vec![self.indices[0].clone(), eigenvector_index.clone()],
            vectors,
        )?;

        Ok(TensorHermitianEigendecomposition {
            eigenvalues,
            eigenvectors,
            eigenvector_index,
        })
    }

    fn read_real_eigenvalues(values: &Self, hermitian_tol: f64) -> Result<Vec<f64>> {
        if values.is_complex() {
            values
                .to_vec::<Complex64>()?
                .into_iter()
                .enumerate()
                .map(|(index, value)| {
                    let imaginary = value.im.abs();
                    let allowed = hermitian_tol * value.norm().max(1.0);
                    if !matches!(
                        imaginary.partial_cmp(&allowed),
                        Some(std::cmp::Ordering::Less) | Some(std::cmp::Ordering::Equal)
                    ) {
            return Err(anyhow::anyhow!("Hermitian eigenvalue {index} has imaginary part {imaginary}, exceeding tolerance {allowed}"));
        };
                    Ok(value.re)
                })
                .collect()
        } else {
            values.to_vec::<f64>().map_err(anyhow::Error::from)
        }
    }

    pub(crate) fn from_diag_inner(
        indices: Vec<DynIndex>,
        payload_inner: EagerTensor,
    ) -> Result<Self> {
        let dims = Self::expected_dims_from_indices(&indices);
        Self::validate_indices(&indices)?;
        Self::validate_diag_dims(&dims)?;
        let payload_len = checked_product(payload_inner.shape())?;
        Self::validate_diag_payload_len(payload_len, &dims)?;
        let axis_classes = Self::diag_axis_classes(dims.len());
        let diag_inner = payload_inner.embed_diag(0, 1)?;
        Self::from_inner_with_axis_classes(indices, diag_inner, axis_classes)
    }

    fn compact_inner_from_logical(
        inner: &EagerTensor,
        axis_classes: &[usize],
    ) -> Result<EagerTensor> {
        let mut payload = inner.clone();
        let mut classes = axis_classes.to_vec();
        while let Some((axis_a, axis_b)) = Self::first_duplicate_pair(&classes) {
            payload = payload.extract_diag(axis_a, axis_b)?;
            classes.remove(axis_b);
        }
        Ok(payload)
    }

    pub(crate) fn from_inner_with_axis_classes(
        indices: Vec<DynIndex>,
        inner: EagerTensor,
        axis_classes: Vec<usize>,
    ) -> Result<Self> {
        Self::from_inner_with_axis_classes_impl(indices, inner, axis_classes, false)
    }

    fn from_inner_with_validated_metadata(
        indices: Vec<DynIndex>,
        inner: EagerTensor,
        axis_classes: Vec<usize>,
    ) -> Result<Self> {
        Self::from_inner_with_axis_classes_impl(indices, inner, axis_classes, true)
    }

    fn from_inner_with_axis_classes_impl(
        indices: Vec<DynIndex>,
        inner: EagerTensor,
        axis_classes: Vec<usize>,
        validated: bool,
    ) -> Result<Self> {
        let dims = profile_pairwise_contract_section("from_inner_expected_dims", || {
            Self::expected_dims_from_indices(&indices)
        });
        let dense_axis_classes = axis_classes.iter().copied().eq(0..indices.len());
        if validated {
            debug_assert!(Self::validate_indices(&indices).is_ok());
            if !dense_axis_classes {
                Self::validate_axis_classes(&axis_classes, indices.len())?;
            }
        } else {
            profile_pairwise_contract_section("from_inner_validate_indices", || {
                Self::validate_indices(&indices)
            })?;
            Self::validate_axis_classes(&axis_classes, indices.len())?;
        }
        if dims != inner.shape() {
            return Err(anyhow::anyhow!(
                "native payload dims {:?} do not match indices dims {:?}",
                inner.shape(),
                dims
            ));
        }
        if Self::is_diag_axis_classes(&axis_classes) {
            profile_pairwise_contract_section("from_inner_validate_diag_dims", || {
                Self::validate_diag_dims(&dims)
            })?;
        }
        let storage = if dense_axis_classes {
            IdxTensorStorage::from_eager_dense(inner, indices.len())
        } else {
            let payload = Self::compact_inner_from_logical(&inner, &axis_classes)?;
            let payload_dims = payload.shape().to_vec();
            IdxTensorStorage::Compact(Arc::new(StructuredPayload {
                payload: Arc::new(payload),
                payload_dims,
                axis_classes,
            }))
        };
        Ok(Self {
            indices,
            storage,
            eager_cache: Self::empty_eager_cache(),
        })
    }

    /// Borrow the indices.
    pub fn indices(&self) -> &[DynIndex] {
        &self.indices
    }

    pub(crate) fn axis_classes(&self) -> &[usize] {
        self.storage.axis_classes()
    }

    #[cfg(feature = "tenferro-cuda")]
    pub(crate) fn deferred_storage_error(&self) -> Option<&TensorStorageError> {
        self.storage.deferred_error()
    }

    #[cfg(feature = "tenferro-cuda")]
    pub(crate) fn cuda_eager_inner(&self) -> Option<&EagerTensor> {
        self.storage
            .eager()
            .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
    }

    /// Borrow the owning eager runtime without materializing or transferring.
    ///
    /// Returns `None` for tensors with no eager value. Used to keep
    /// contraction results in their operands' context.
    pub(crate) fn eager_runtime(&self) -> Option<Arc<EagerRuntime>> {
        self.storage
            .eager()
            .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
            .map(|inner| Arc::clone(inner.runtime()))
    }

    /// Check whether this tensor's eager value is device-resident.
    ///
    /// Metadata-only: inspects placement without reading data, transferring,
    /// or consulting any execution context. Used to route factorization to the
    /// context-scoped path and to reject context-free calls on resident inputs
    /// before algorithm work begins.
    #[cfg(feature = "tenferro-cuda")]
    pub(crate) fn is_cuda_resident(&self) -> bool {
        self.cuda_eager_inner().is_some_and(|inner| {
            inner.tensor_read().placement().memory_kind == tenferro_tensor::MemoryKind::Device
        })
    }

    /// Check whether this tensor's eager value is device-resident.
    ///
    /// Without the CUDA feature no tensor can be device-resident.
    #[cfg(not(feature = "tenferro-cuda"))]
    pub(crate) fn is_cuda_resident(&self) -> bool {
        false
    }

    #[cfg(feature = "tenferro-cuda")]
    pub(crate) fn cuda_duplicate_native(&self) -> Result<NativeTensor> {
        Ok(self.try_materialized_inner()?.duplicate_value()?)
    }

    /// Enable reverse-mode AD tracking on this tensor by creating a tracked leaf.
    /// # Errors
    /// Returns an error when the tensor is not a scalar (a rank mismatch) or the
    /// AD backend cannot track the tensor's dtype.
    ///
    pub fn enable_grad(self) -> std::result::Result<Self, IdxTensorError> {
        self.ensure_storage_ready()?;
        // Keep the eager payload when available: compact Storage currently
        // stores only f64/C64 and must not promote f32/C32 leaves before AD.
        let eager_payload = self
            .storage
            .eager()
            .or_else(|| self.eager_cache.get().map(AsRef::as_ref))
            .filter(|inner| inner.shape() == self.storage.payload_dims());
        let payload = match eager_payload {
            Some(inner) => inner.duplicate_value()?,
            None => {
                let materialized = self.storage.materialize(self.indices.len())?;
                storage_payload_native(materialized.as_ref())
                    .context("IdxTensor::enable_grad failed")?
            }
        };
        let payload_dims = self.storage.payload_dims().to_vec();
        let axis_classes = self.storage.axis_classes().to_vec();
        let tracked = Arc::new(EagerTensor::requires_grad_in(
            payload,
            default_eager_ctx()?,
        )?);
        let storage = if axis_classes == Self::dense_axis_classes(self.indices.len()) {
            IdxTensorStorage::Eager {
                inner: tracked,
                axis_classes,
            }
        } else {
            IdxTensorStorage::Compact(Arc::new(StructuredPayload {
                payload: tracked,
                payload_dims,
                axis_classes,
            }))
        };
        Ok(Self {
            indices: self.indices,
            storage,
            eager_cache: Self::empty_eager_cache(),
        })
    }

    /// Report whether this tensor participates in gradient tracking.
    pub fn tracks_grad(&self) -> bool {
        self.storage.eager().is_some_and(EagerTensor::tracks_grad)
            || self
                .eager_cache
                .get()
                .is_some_and(|inner| inner.tracks_grad())
    }

    /// Return the accumulated gradient, if one has been stored.
    /// # Errors
    /// Returns an error when the tensor is not a tracked leaf or the gradient is
    /// unavailable for the tensor's dtype (an unavailable-gradient failure).
    ///
    pub fn grad(&self) -> std::result::Result<Option<Self>, IdxTensorError> {
        if let Some(value) = self.tracked_compact_payload_value() {
            let Some(gradient) = value.payload.grad()? else {
                return Ok(None);
            };
            let gradient_shape = gradient.shape().to_vec();
            let gradient_tensor = gradient.to_tensor()?;
            if self.compact_payload_is_logical_dense(&value.payload_dims) {
                return Ok(Some(Self::from_native_with_axis_classes(
                    self.indices.clone(),
                    gradient_tensor,
                    value.axis_classes.clone(),
                )?));
            }
            if gradient_shape != value.payload_dims {
                return Err(anyhow::anyhow!(
                    "gradient payload dims {:?} do not match {:?}",
                    gradient_shape,
                    value.payload_dims
                )
                .into());
            }
            let gradient = EagerTensor::from_tensor_in(gradient_tensor, default_eager_ctx()?)?;
            return Ok(Some(Self::from_structured_payload_inner(
                self.indices.clone(),
                gradient,
                value.payload_dims.clone(),
                value.axis_classes.clone(),
            )?));
        }

        let Some(gradient) = self.try_materialized_inner()?.grad()? else {
            return Ok(None);
        };
        Ok(Some(Self::from_native_with_axis_classes(
            self.indices.clone(),
            gradient.to_tensor()?,
            self.storage.axis_classes().to_vec(),
        )?))
    }

    /// Clear the accumulated gradient stored for this tensor.
    /// # Errors
    /// Returns an error when the tensor is not a tracked leaf (a missing-graph
    /// failure).
    ///
    pub fn clear_grad(&self) -> std::result::Result<(), IdxTensorError> {
        self.ensure_storage_ready()?;
        if let Some(value) = self.tracked_compact_payload_value() {
            value.payload.clear_grad()?;
        }
        if let Some(inner) = self.storage.eager() {
            inner.clear_grad()?;
        }
        if let Some(inner) = self.eager_cache.get() {
            inner.clear_grad()?;
        }
        Ok(())
    }

    /// Run reverse-mode autodiff from this scalar tensor.
    /// # Errors
    /// Returns an error when the tensor is not a scalar (a rank mismatch) or the
    /// reverse pass fails (a graph failure).
    ///
    pub fn backward(&self) -> std::result::Result<(), IdxTensorError> {
        if let Some(value) = self.tracked_compact_payload_value() {
            return value.payload.backward().map(|_| ()).map_err(|e| {
                IdxTensorError::from(anyhow::anyhow!("IdxTensor::backward failed: {e}"))
            });
        }
        self.try_materialized_inner()?
            .backward()
            .map(|_| ())
            .map_err(|e| IdxTensorError::from(anyhow::anyhow!("IdxTensor::backward failed: {e}")))
    }

    /// Detach this tensor from the reverse graph.
    /// # Errors
    /// Returns an error when the tensor is not a tracked leaf (a missing-graph
    /// failure).
    ///
    pub fn detach(&self) -> std::result::Result<Self, IdxTensorError> {
        Self::from_inner_with_axis_classes(
            self.indices.clone(),
            self.try_materialized_inner()?.detach(),
            self.storage.axis_classes().to_vec(),
        )
        .map_err(IdxTensorError::from)
    }

    /// Check if this tensor is already in canonical form.
    pub fn is_simple(&self) -> bool {
        true
    }

    /// Materialize the primal payload as a compact storage snapshot.
    ///
    /// The eager payload remains authoritative for `f32`/`c32` and tracked
    /// structured tensors; this method is a fallible bridge to compact storage.
    ///
    /// # Errors
    ///
    /// Returns [`TensorStorageError`] when an eager backend payload cannot be
    /// converted to compact storage, when its dtype is `f32`/`c32`, or when a
    /// deferred eager operation failed.
    pub fn to_storage(&self) -> std::result::Result<Arc<Storage>, TensorStorageError> {
        self.storage.materialize(self.indices.len())
    }

    /// Materializes and returns a compact storage snapshot.
    ///
    /// # Errors
    /// Returns an error when the compact storage cannot be materialized (a
    /// backend failure).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::StorageKind;
    ///
    /// let tensor = IdxTensor::from_dense(
    ///     vec![DynIndex::new_dyn(2)],
    ///     vec![1.0_f64, 2.0],
    /// )
    /// .unwrap();
    /// assert_eq!(tensor.storage().unwrap().storage_kind(), StorageKind::Dense);
    /// ```
    pub fn storage(&self) -> std::result::Result<Arc<Storage>, TensorStorageError> {
        self.storage.materialize(self.indices.len())
    }

    /// Return the logical storage layout without materializing compact storage.
    ///
    /// For `f32` and `c32`, the eager representation is authoritative because
    /// compact [`Storage`] supports only `f64` and `c64`.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::StorageKind;
    ///
    /// let tensor = IdxTensor::from_diag(
    ///     vec![DynIndex::new_dyn(2), DynIndex::new_dyn(2)],
    ///     vec![1.0_f32, 2.0],
    /// )
    /// .unwrap();
    /// assert_eq!(tensor.storage_kind(), StorageKind::Diagonal);
    /// ```
    pub fn storage_kind(&self) -> StorageKind {
        self.storage.storage_kind()
    }

    /// Return the exact eager scalar dtype without reading tensor values.
    ///
    /// This feature-gated accessor is used by the CUDA TreeTN boundary to
    /// reject mixed dtypes before the first device contraction.
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the dtype cannot be determined from the
    /// tensor metadata.
    ///
    /// # Examples
    ///
    /// ```
    /// # #[cfg(feature = "tenferro-cuda")]
    /// # {
    /// use tensor4all_core::{DynIndex, IdxTensor, IdxTensorError};
    ///
    /// let tensor = IdxTensor::from_dense(vec![DynIndex::new_dyn(1)], vec![1.0_f32]).unwrap();
    /// let dtype = tensor.cuda_dtype().unwrap();
    /// assert_eq!(dtype, tenferro::DType::F32);
    /// let accessor: fn(&IdxTensor) -> Result<tenferro::DType, IdxTensorError> =
    ///     IdxTensor::cuda_dtype;
    /// assert_eq!(
    ///     std::mem::size_of_val(&accessor),
    ///     std::mem::size_of::<fn(&IdxTensor) -> Result<tenferro::DType, IdxTensorError>>(),
    /// );
    /// # }
    /// ```
    #[cfg(feature = "tenferro-cuda")]
    pub fn cuda_dtype(&self) -> std::result::Result<DType, IdxTensorError> {
        self.scalar_dtype().map_err(IdxTensorError::from)
    }

    /// Sum all elements, returning `AnyScalar`.
    ///
    /// # Errors
    /// Returns an error when the reduction fails (a backend or scalar-extraction
    /// failure).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(3);
    /// let t = IdxTensor::from_dense(vec![i], vec![1.0, 2.0, 3.0]).unwrap();
    /// let s = t.sum().unwrap();
    /// assert!((s.real() - 6.0).abs() < 1e-12);
    /// ```
    pub fn sum(&self) -> std::result::Result<AnyScalar, IdxTensorError> {
        self.ensure_storage_ready()?;
        if self.indices.is_empty() {
            return AnyScalar::from_tensor(self.clone()).map_err(IdxTensorError::from);
        }
        if let Some(payload) = self.storage.eager().filter(|payload| payload.tracks_grad()) {
            let axes: Vec<usize> = (0..payload.shape().len()).collect();
            let reduced = payload.reduce_sum(Some(&axes))?;
            return AnyScalar::from_tensor(Self::from_inner(Vec::new(), reduced)?)
                .map_err(IdxTensorError::from);
        }
        self.storage.sum_scalar().map_err(IdxTensorError::from)
    }

    /// Extract the scalar value from a 0-dimensional tensor (or 1-element tensor).
    ///
    /// This is similar to Julia's `only()` function.
    ///
    /// # Errors
    /// Returns an error when the tensor is not rank-0 and does not contain exactly
    /// one element (a rank mismatch).
    /// # Panics
    ///
    /// Panics if the tensor has more than one element.
    ///
    /// # Example
    ///
    /// ```
    /// use tensor4all_core::{IdxTensor, AnyScalar};
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// // Create a scalar tensor (0 dimensions, 1 element)
    /// let indices: Vec<Index<DynId>> = vec![];
    /// let tensor: IdxTensor = IdxTensor::from_dense(indices, vec![42.0]).unwrap();
    ///
    /// assert_eq!(tensor.only().unwrap().real(), 42.0);
    /// ```
    pub fn only(&self) -> std::result::Result<AnyScalar, IdxTensorError> {
        let dims = self.dims();
        let total_size = checked_product(&dims)?;
        if !(total_size == 1 || dims.is_empty()) {
            return Err(anyhow::anyhow!(
                "only() requires a scalar tensor (1 element), got {} elements with dims {:?}",
                if dims.is_empty() { 1 } else { total_size },
                dims
            )
            .into());
        };
        self.sum()
    }

    /// Permute the tensor dimensions using the given new indices order.
    ///
    /// This is the main permutation method that takes the desired new indices
    /// and automatically computes the corresponding permutation of dimensions
    /// and data. The new indices must be a permutation of the original indices
    /// (matched by full index identity).
    ///
    /// # Arguments
    /// * `new_indices` - The desired new indices order. Must be a permutation
    ///
    ///   of `self.indices` (matched by full index identity).
    ///
    /// # Errors
    /// Returns an error when `new_order` does not contain exactly the tensor's
    /// indices (an index-set mismatch or a missing-index failure).
    /// # Panics
    /// Panics if `new_indices.len() != self.indices.len()`, if any full index
    /// identity doesn't match, or if there are duplicate indices.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// // Create a 2×3 tensor
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let indices = vec![i.clone(), j.clone()];
    /// let tensor: IdxTensor = IdxTensor::from_dense(indices, vec![0.0; 6]).unwrap();
    ///
    /// // Permute to 3×2: swap the two dimensions by providing new indices order
    /// let permuted = tensor.permute_indices(&[j, i]).unwrap();
    /// assert_eq!(permuted.dims(), vec![3, 2]);
    /// ```
    pub fn permute_indices(
        &self,
        new_indices: &[DynIndex],
    ) -> std::result::Result<Self, IdxTensorError> {
        // Compute permutation by full index equality
        let perm = compute_permutation_from_indices(&self.indices, new_indices)?;
        if perm.iter().copied().eq(0..perm.len()) {
            return Ok(Self {
                indices: new_indices.to_vec(),
                storage: self.storage.clone(),
                eager_cache: Arc::clone(&self.eager_cache),
            });
        }

        let permuted = self.try_materialized_inner()?.transpose(&perm)?;
        let axis_classes = self.permute_axis_classes(&perm);
        Self::from_inner_with_axis_classes(new_indices.to_vec(), permuted, axis_classes)
            .map_err(IdxTensorError::from)
    }

    /// Permute the tensor dimensions, returning a new tensor.
    ///
    /// This method reorders the indices, dimensions, and data according to the
    /// given permutation. The permutation specifies which old axis each new
    /// axis corresponds to: `new_axis[i] = old_axis[perm[i]]`.
    ///
    /// # Arguments
    /// * `perm` - The permutation: `perm[i]` is the old axis index for new axis `i`
    ///
    /// # Errors
    /// Returns an error when `new_order` does not contain exactly the tensor's
    /// indices (an index-set mismatch or a missing-index failure).
    /// # Panics
    /// Panics if `perm.len() != self.indices.len()` or if the permutation is invalid.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// // Create a 2×3 tensor
    /// let indices = vec![
    ///     Index::new_dyn(2),
    ///     Index::new_dyn(3),
    /// ];
    /// let tensor: IdxTensor = IdxTensor::from_dense(indices, vec![0.0; 6]).unwrap();
    ///
    /// // Permute to 3×2: swap the two dimensions
    /// let permuted = tensor.permute(&[1, 0]).unwrap();
    /// assert_eq!(permuted.dims(), vec![3, 2]);
    /// ```
    pub fn permute(&self, perm: &[usize]) -> std::result::Result<Self, IdxTensorError> {
        if !(perm.len() == self.indices.len()) {
            return Err(anyhow::anyhow!("permutation length must match tensor rank").into());
        };
        let mut seen = HashSet::new();
        for &axis in perm {
            if !(axis < self.indices.len()) {
                return Err(anyhow::anyhow!("permutation axis {axis} out of range").into());
            };
            if !(seen.insert(axis)) {
                return Err(anyhow::anyhow!("duplicate axis {axis} in permutation").into());
            };
        }
        if perm.iter().copied().eq(0..perm.len()) {
            return Ok(self.clone());
        }

        // Permute indices
        let new_indices: Vec<DynIndex> = perm.iter().map(|&i| self.indices[i].clone()).collect();
        let permuted = self.try_materialized_inner()?.transpose(perm)?;
        let axis_classes = self.permute_axis_classes(perm);
        Self::from_inner_with_axis_classes(new_indices, permuted, axis_classes)
            .map_err(IdxTensorError::from)
    }

    pub(crate) fn try_contract_pairwise_default(&self, other: &Self) -> Result<Self> {
        self.try_contract_pairwise_default_with_options(other, PairwiseContractionOptions::new())
    }

    pub(crate) fn try_contract_pairwise_default_with_options(
        &self,
        other: &Self,
        options: PairwiseContractionOptions,
    ) -> Result<Self> {
        // Device-resident structured/diagonal operands cannot enter the native
        // host session used below; route them through the owning-runtime eager
        // plan executor (which also serves N-ary resident contraction).
        // Dense-dense resident pairs keep the pairwise eager path, and host
        // operands keep this function unchanged.
        #[cfg(feature = "tenferro-cuda")]
        if self.is_cuda_resident()
            && other.is_cuda_resident()
            && self.should_use_structured_payload_contract(other)
        {
            return super::contract::contract_pair_via_plan(self, other, options);
        }
        let self_indices = profile_pairwise_contract_section("operand_indices", || {
            self.operand_indices_for_contraction(options.lhs_conj)
        });
        let other_indices = profile_pairwise_contract_section("operand_indices", || {
            other.operand_indices_for_contraction(options.rhs_conj)
        });
        let self_dims = profile_pairwise_contract_section("expected_dims", || {
            Self::expected_dims_from_indices(&self_indices)
        });
        let other_dims = profile_pairwise_contract_section("expected_dims", || {
            Self::expected_dims_from_indices(&other_indices)
        });
        let spec = profile_pairwise_contract_section("prepare_contraction", || {
            prepare_contraction(&self_indices, &self_dims, &other_indices, &other_dims)
        })
        .context("contraction preparation failed")?;
        let result_axis_classes = profile_pairwise_contract_section("result_axis_classes", || {
            Self::binary_contraction_axis_classes(
                self.storage.axis_classes(),
                &spec.axes_a,
                other.storage.axis_classes(),
                &spec.axes_b,
            )
        })?;

        if profile_pairwise_contract_section("structured_check", || {
            self.should_use_structured_payload_contract(other)
        }) {
            if options.has_conj() {
                let lhs = if options.lhs_conj {
                    self.conj()
                } else {
                    self.clone()
                };
                let rhs = if options.rhs_conj {
                    other.conj()
                } else {
                    other.clone()
                };
                return profile_pairwise_contract_section("structured_conj_fallback", || {
                    lhs.try_contract_pairwise_default(&rhs)
                });
            }
            return profile_pairwise_contract_section("structured_payload_contract", || {
                self.contract_structured_payloads(
                    other,
                    spec.result_indices.into_vec(),
                    &spec.axes_a,
                    &spec.axes_b,
                )
            });
        }

        if self.indices.is_empty() && other.indices.is_empty() {
            if options.has_conj() {
                let lhs = if options.lhs_conj {
                    self.conj()
                } else {
                    self.clone()
                };
                let rhs = if options.rhs_conj {
                    other.conj()
                } else {
                    other.clone()
                };
                return lhs.try_contract_pairwise_default(&rhs);
            }
            let result = profile_pairwise_contract_section("scalar_mul", || {
                Ok::<_, anyhow::Error>(
                    self.try_materialized_inner()?
                        .mul(other.try_materialized_inner()?)?,
                )
            })?;
            return profile_pairwise_contract_section("from_inner", || {
                Self::from_inner(spec.result_indices.into_vec(), result)
            });
        }

        let self_dtype = self.try_materialized_inner()?.dtype();
        let other_dtype = other.try_materialized_inner()?.dtype();
        if self_dtype != other_dtype {
            if options.has_conj() {
                let lhs = if options.lhs_conj {
                    self.conj()
                } else {
                    self.clone()
                };
                let rhs = if options.rhs_conj {
                    other.conj()
                } else {
                    other.clone()
                };
                return lhs.try_contract_pairwise_default(&rhs);
            }
            let self_native = self.try_materialized_inner()?.duplicate_value()?;
            let other_native = other.try_materialized_inner()?.duplicate_value()?;
            let result_native = profile_pairwise_contract_section("native_contract", || {
                contract_native_tensor(&self_native, &spec.axes_a, &other_native, &spec.axes_b)
            })?;
            return profile_pairwise_contract_section("from_native", || {
                Self::from_untracked_native_with_axis_classes(
                    spec.result_indices.into_vec(),
                    result_native,
                    result_axis_classes,
                )
            });
        }

        let config = profile_pairwise_contract_section("build_dot_general_config", || {
            Self::binary_dot_general_config(&spec.axes_a, &spec.axes_b)
        })?;
        let result = profile_pairwise_contract_section("dot_general_with_conj", || {
            let lhs = profile_pairwise_contract_section("lhs_try_materialized_inner", || {
                self.try_materialized_inner()
            })?;
            let rhs = profile_pairwise_contract_section("rhs_try_materialized_inner", || {
                other.try_materialized_inner()
            })?;
            profile_pairwise_contract_section("dot_general_execute", || {
                lhs.dot_general_with_conj(rhs, config, options.lhs_conj, options.rhs_conj)
            })
            .map_err(anyhow::Error::from)
        })?;
        record_pairwise_contract_profile_bytes(
            "dot_general_output",
            tensor_profile_bytes(result.dtype(), result.shape()),
        );
        profile_pairwise_contract_section("from_inner_axis_classes", || {
            Self::from_inner_with_axis_classes(
                spec.result_indices.into_vec(),
                result,
                result_axis_classes,
            )
        })
    }

    // Provenance: this retained-axis dense path is a tensor4all-specific
    // implementation of the batch needed by SRC. Its output ordering follows
    // `defaults::contract::build_contraction_plan`; it is not present in
    // RandomMPOMPS and is labelled `[AI-Supplied]` in the SRC audit.
    pub(crate) fn try_contract_pairwise_retaining(
        &self,
        other: &Self,
        retained_indices: &[DynIndex],
    ) -> Result<Self> {
        if retained_indices.is_empty()
            || self.should_use_structured_payload_contract(other)
            || self.try_materialized_inner()?.dtype() != other.try_materialized_inner()?.dtype()
        {
            let options = ContractionOptions::new().with_retain_indices(retained_indices);
            return super::contract::contract_with_options(&[self, other], options)
                .map_err(anyhow::Error::from);
        }

        let common = common_ind_positions(&self.indices, &other.indices);
        if common.is_empty() {
            let options = ContractionOptions::new().with_retain_indices(retained_indices);
            return super::contract::contract_with_options(&[self, other], options)
                .map_err(anyhow::Error::from);
        }

        let mut retained_pairs = Vec::with_capacity(retained_indices.len());
        for retained in retained_indices {
            let Some(pos_a) = self.indices.iter().position(|index| index == retained) else {
                let options = ContractionOptions::new().with_retain_indices(retained_indices);
                return super::contract::contract_with_options(&[self, other], options)
                    .map_err(anyhow::Error::from);
            };
            let Some(pos_b) = other.indices.iter().position(|index| index == retained) else {
                let options = ContractionOptions::new().with_retain_indices(retained_indices);
                return super::contract::contract_with_options(&[self, other], options)
                    .map_err(anyhow::Error::from);
            };
            if !self.indices[pos_a].is_contractable(&other.indices[pos_b]) {
                let options = ContractionOptions::new().with_retain_indices(retained_indices);
                return super::contract::contract_with_options(&[self, other], options)
                    .map_err(anyhow::Error::from);
            }
            retained_pairs.push((pos_a, pos_b));
        }

        let retained_pairs_set: HashSet<(usize, usize)> = retained_pairs.iter().copied().collect();
        let mut contracting_a = Vec::with_capacity(common.len());
        let mut contracting_b = Vec::with_capacity(common.len());
        for &(pos_a, pos_b) in &common {
            if !retained_pairs_set.contains(&(pos_a, pos_b)) {
                contracting_a.push(pos_a);
                contracting_b.push(pos_b);
            }
        }

        let config = DotGeneralConfig {
            lhs_contracting_dims: contracting_a.clone(),
            rhs_contracting_dims: contracting_b.clone(),
            lhs_batch_dims: retained_pairs.iter().map(|&(pos_a, _)| pos_a).collect(),
            rhs_batch_dims: retained_pairs.iter().map(|&(_, pos_b)| pos_b).collect(),
        };
        let result = self
            .try_materialized_inner()?
            .dot_general_with_conj(other.try_materialized_inner()?, config, false, false)
            .map_err(|error| anyhow::anyhow!("retained pairwise contraction failed: {error}"))?;

        // dot_general emits [lhs_free, rhs_free, batch]. Restore the index
        // order used by contract_with_options, where retained axes occur at
        // their first operand position. The permutation is a single dense
        // copy after the batched GEMM and is much cheaper than one einsum per
        // probe column.
        let contracting_a_set: HashSet<usize> = contracting_a.into_iter().collect();
        let contracting_b_set: HashSet<usize> = contracting_b.into_iter().collect();
        let retained_a_set: HashSet<usize> = retained_pairs.iter().map(|&(a, _)| a).collect();
        let retained_b_set: HashSet<usize> = retained_pairs.iter().map(|&(_, b)| b).collect();
        let mut current_indices = self
            .indices
            .iter()
            .enumerate()
            .filter(|(axis, _)| !retained_a_set.contains(axis) && !contracting_a_set.contains(axis))
            .map(|(_, index)| index.clone())
            .chain(
                other
                    .indices
                    .iter()
                    .enumerate()
                    .filter(|(axis, _)| {
                        !retained_b_set.contains(axis) && !contracting_b_set.contains(axis)
                    })
                    .map(|(_, index)| index.clone()),
            )
            .collect::<Vec<_>>();
        current_indices.extend(
            retained_pairs
                .iter()
                .map(|&(pos_a, _)| self.indices[pos_a].clone()),
        );

        let mut desired_indices = self
            .indices
            .iter()
            .enumerate()
            .filter(|(axis, _)| !contracting_a_set.contains(axis))
            .map(|(_, index)| index.clone())
            .chain(
                other
                    .indices
                    .iter()
                    .enumerate()
                    .filter(|(axis, _)| {
                        !contracting_b_set.contains(axis) && !retained_b_set.contains(axis)
                    })
                    .map(|(_, index)| index.clone()),
            )
            .collect::<Vec<_>>();
        if desired_indices.is_empty() {
            desired_indices = current_indices.clone();
        }
        let result = if current_indices == desired_indices {
            result
        } else {
            let permutation = desired_indices
                .iter()
                .map(|desired| {
                    current_indices
                        .iter()
                        .position(|current| current == desired)
                        .ok_or_else(|| {
                            anyhow::anyhow!(
                                "retained pairwise result is missing index {:?}",
                                desired
                            )
                        })
                })
                .collect::<Result<Vec<_>>>()?;
            result.transpose(&permutation)?
        };
        Self::from_inner(desired_indices, result)
    }

    pub(crate) fn try_tensordot_pairwise_explicit(
        &self,
        other: &Self,
        pairs: &[(DynIndex, DynIndex)],
    ) -> Result<Self> {
        use crate::index_ops::ContractionError;

        let self_dims = Self::expected_dims_from_indices(&self.indices);
        let other_dims = Self::expected_dims_from_indices(&other.indices);
        let spec = prepare_contraction_pairs(
            &self.indices,
            &self_dims,
            &other.indices,
            &other_dims,
            pairs,
        )
        .map_err(|e| match e {
            ContractionError::NoCommonIndices => {
                anyhow::anyhow!("tensordot: No pairs specified for contraction")
            }
            ContractionError::BatchContractionNotImplemented => anyhow::anyhow!(
                "tensordot: Common index found but not in contraction pairs. \
                         Batch contraction is not yet implemented."
            ),
            ContractionError::IndexNotFound { tensor } => {
                anyhow::anyhow!("tensordot: Index not found in {} tensor", tensor)
            }
            ContractionError::DimensionMismatch {
                pos_a,
                pos_b,
                dim_a,
                dim_b,
            } => anyhow::anyhow!(
                "tensordot: Dimension mismatch: self[{}]={} != other[{}]={}",
                pos_a,
                dim_a,
                pos_b,
                dim_b
            ),
            ContractionError::DuplicateAxis { tensor, pos } => {
                anyhow::anyhow!("tensordot: Duplicate axis {} in {} tensor", pos, tensor)
            }
        })?;
        let result_axis_classes = Self::binary_contraction_axis_classes(
            self.storage.axis_classes(),
            &spec.axes_a,
            other.storage.axis_classes(),
            &spec.axes_b,
        )?;

        if self.should_use_structured_payload_contract(other) {
            return self.contract_structured_payloads(
                other,
                spec.result_indices.into_vec(),
                &spec.axes_a,
                &spec.axes_b,
            );
        }

        if self.indices.is_empty() && other.indices.is_empty() {
            let result = self
                .try_materialized_inner()?
                .mul(other.try_materialized_inner()?)
                .map_err(|e| anyhow::anyhow!("tensordot scalar multiply failed: {e}"))?;
            return Self::from_inner(spec.result_indices.into_vec(), result);
        }

        let self_dtype = self.try_materialized_inner()?.dtype();
        let other_dtype = other.try_materialized_inner()?.dtype();
        if self_dtype != other_dtype {
            let self_native = self.try_materialized_inner()?.duplicate_value()?;
            let other_native = other.try_materialized_inner()?.duplicate_value()?;
            let result_native =
                contract_native_tensor(&self_native, &spec.axes_a, &other_native, &spec.axes_b)?;
            return Self::from_untracked_native_with_axis_classes(
                spec.result_indices.into_vec(),
                result_native,
                result_axis_classes,
            );
        }

        let subscripts = Self::build_binary_einsum_subscripts(
            self.indices.len(),
            &spec.axes_a,
            other.indices.len(),
            &spec.axes_b,
        )?;
        let result = [
            self.try_materialized_inner()?,
            other.try_materialized_inner()?,
        ]
        .einsum_subscripts(&subscripts)
        .map_err(|e| anyhow::anyhow!("tensordot failed: {e}"))?;
        Self::from_inner_with_axis_classes(
            spec.result_indices.into_vec(),
            result,
            result_axis_classes,
        )
    }

    pub(crate) fn try_outer_product_pairwise(&self, other: &Self) -> Result<Self> {
        use anyhow::Context;

        // Check for common indices - outer product should have none
        let common_positions = common_ind_positions(&self.indices, &other.indices);
        if !common_positions.is_empty() {
            let common_ids: Vec<_> = common_positions
                .iter()
                .map(|(pos_a, _)| self.indices[*pos_a].id())
                .collect();
            return Err(anyhow::anyhow!(
                "outer_product: tensors have common indices {:?}. \
                 Use tensordot to contract common indices, or use sim() to replace \
                 indices with fresh IDs before computing outer product.",
                common_ids
            ))
            .context("outer_product: common indices found");
        }

        // Build result indices and dimensions
        let mut result_indices = self.indices.clone();
        result_indices.extend(other.indices.iter().cloned());
        let result_axis_classes = Self::binary_contraction_axis_classes(
            self.storage.axis_classes(),
            &[],
            other.storage.axis_classes(),
            &[],
        )?;
        if self.should_use_structured_payload_contract(other) {
            return self.contract_structured_payloads(other, result_indices, &[], &[]);
        }
        let self_dtype = self.try_materialized_inner()?.dtype();
        let other_dtype = other.try_materialized_inner()?.dtype();
        if self_dtype != other_dtype {
            let self_native = self.try_materialized_inner()?.duplicate_value()?;
            let other_native = other.try_materialized_inner()?.duplicate_value()?;
            let result_native = contract_native_tensor(&self_native, &[], &other_native, &[])?;
            return Self::from_untracked_native_with_axis_classes(
                result_indices,
                result_native,
                result_axis_classes,
            );
        }

        let subscripts = Self::build_binary_einsum_subscripts(
            self.indices.len(),
            &[],
            other.indices.len(),
            &[],
        )?;
        let result = [
            self.try_materialized_inner()?,
            other.try_materialized_inner()?,
        ]
        .einsum_subscripts(&subscripts)
        .map_err(|e| anyhow::anyhow!("outer_product failed: {e}"))?;
        Self::from_inner_with_axis_classes(result_indices, result, result_axis_classes)
    }
}

// ============================================================================
// Random tensor generation
// ============================================================================

impl IdxTensor {
    /// Create a random tensor with values from standard normal distribution (generic over scalar type).
    ///
    /// For `f64`, each element is drawn from the standard normal distribution.
    /// For `Complex64`, both real and imaginary parts are drawn independently.
    ///
    /// # Type Parameters
    /// * `T` - The scalar element type (must implement [`RandomScalar`])
    /// * `R` - The random number generator type
    ///
    /// # Arguments
    /// * `rng` - Random number generator
    /// * `indices` - The indices for the tensor
    ///
    /// # Errors
    /// Returns an error when the dimension product overflows (an overflow failure)
    /// or the backend cannot generate the requested scalar type.
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    /// use rand::SeedableRng;
    /// use rand_chacha::ChaCha8Rng;
    ///
    /// let mut rng = ChaCha8Rng::seed_from_u64(42);
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let tensor: IdxTensor = IdxTensor::random::<f64, _>(&mut rng, vec![i, j]).unwrap();
    /// assert_eq!(tensor.dims(), vec![2, 3]);
    /// ```
    pub fn random<T: RandomScalar, R: Rng>(
        rng: &mut R,
        indices: Vec<DynIndex>,
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims: Vec<usize> = indices.iter().map(|idx| idx.dim()).collect();
        let size = checked_product(&dims)?;
        let data: Vec<T> = (0..size).map(|_| T::random_value(rng)).collect();
        Self::from_dense(indices, data)
    }
}

impl IdxTensor {
    /// Add two tensors element-wise.
    ///
    /// The tensors must have the same full index set (including tags and prime
    /// levels). If the indices are in a different order, the other tensor will
    /// be permuted to match `self`.
    ///
    /// # Arguments
    /// * `other` - The tensor to add
    ///
    /// # Returns
    /// A new tensor representing `self + other`, or an error if:
    /// - The tensors have different index sets
    /// - The dimensions don't match
    /// - Storage types are incompatible
    ///
    /// # Errors
    /// Returns an error when the two tensors have different index sets (an
    /// index-set mismatch) or the arithmetic reports a failure.
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    ///
    /// let indices_a = vec![i.clone(), j.clone()];
    /// let data_a = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
    /// let tensor_a: IdxTensor = IdxTensor::from_dense(indices_a, data_a).unwrap();
    ///
    /// let indices_b = vec![i.clone(), j.clone()];
    /// let data_b = vec![1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
    /// let tensor_b: IdxTensor = IdxTensor::from_dense(indices_b, data_b).unwrap();
    ///
    /// let sum = tensor_a.add(&tensor_b).unwrap();
    /// // sum = [[2, 3, 4], [5, 6, 7]]
    /// ```
    pub fn add(&self, other: &Self) -> std::result::Result<Self, IdxTensorError> {
        // Validate that both tensors have the same number of indices
        if self.indices.len() != other.indices.len() {
            return Err(anyhow::anyhow!(
                "Index count mismatch: self has {} indices, other has {}",
                self.indices.len(),
                other.indices.len()
            )
            .into());
        }

        // Validate that both tensors have the same set of indices
        let self_set: HashSet<_> = self.indices.iter().collect();
        let other_set: HashSet<_> = other.indices.iter().collect();

        if self_set != other_set {
            return Err(
                anyhow::anyhow!("Index set mismatch: tensors must have the same indices").into(),
            );
        }

        // Permute other to match self's index order (no-op if already aligned)
        let other_aligned = other.permute_indices(&self.indices)?;

        // Validate dimensions match after alignment
        let self_expected_dims = Self::expected_dims_from_indices(&self.indices);
        let other_expected_dims = Self::expected_dims_from_indices(&other_aligned.indices);
        if self_expected_dims != other_expected_dims {
            use crate::TagSetLike;
            let fmt = |indices: &[DynIndex]| -> Vec<String> {
                indices
                    .iter()
                    .map(|idx| {
                        let tags: Vec<String> = idx.tags().iter().collect();
                        format!("{:?}(dim={},tags={:?})", idx.id(), idx.dim(), tags)
                    })
                    .collect()
            };
            return Err(anyhow::anyhow!(
                "Dimension mismatch after alignment.\n\
                 self: dims={:?}, indices(order)={:?}\n\
                 other_aligned: dims={:?}, indices(order)={:?}",
                self_expected_dims,
                fmt(&self.indices),
                other_expected_dims,
                fmt(&other_aligned.indices)
            )
            .into());
        }

        self.axpby(
            AnyScalar::new_real(1.0),
            &other_aligned,
            AnyScalar::new_real(1.0),
        )
    }

    /// Compute a linear combination: `a * self + b * other`.
    ///
    /// Both tensors must have the same full set of indices (including tags and
    /// prime levels). If indices are in a different order, `other` is automatically permuted
    /// to match `self`.
    ///
    /// # Errors
    /// Returns an error when the tensors have different index sets (an index-set
    /// mismatch) or the arithmetic reports a failure.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{AnyScalar, DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let a = IdxTensor::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
    /// let b = IdxTensor::from_dense(vec![i.clone()], vec![3.0, 4.0]).unwrap();
    ///
    /// // 2*a + 3*b = [2+9, 4+12] = [11, 16]
    /// let result = a.axpby(AnyScalar::new_real(2.0), &b, AnyScalar::new_real(3.0)).unwrap();
    /// let data = result.to_vec::<f64>().unwrap();
    /// assert!((data[0] - 11.0).abs() < 1e-12);
    /// assert!((data[1] - 16.0).abs() < 1e-12);
    /// ```
    pub fn axpby(
        &self,
        a: AnyScalar,
        other: &Self,
        b: AnyScalar,
    ) -> std::result::Result<Self, IdxTensorError> {
        // Validate that both tensors have the same number of indices.
        if self.indices.len() != other.indices.len() {
            return Err(anyhow::anyhow!(
                "Index count mismatch: self has {} indices, other has {}",
                self.indices.len(),
                other.indices.len()
            )
            .into());
        }

        // Validate that both tensors have the same set of indices.
        let self_set: HashSet<_> = self.indices.iter().collect();
        let other_set: HashSet<_> = other.indices.iter().collect();
        if self_set != other_set {
            return Err(
                anyhow::anyhow!("Index set mismatch: tensors must have the same indices").into(),
            );
        }

        // Align other tensor axis order to self.
        let other_aligned = other.permute_indices(&self.indices)?;

        // Validate dimensions match after alignment.
        let self_expected_dims = Self::expected_dims_from_indices(&self.indices);
        let other_expected_dims = Self::expected_dims_from_indices(&other_aligned.indices);
        if self_expected_dims != other_expected_dims {
            return Err(anyhow::anyhow!(
                "Dimension mismatch after alignment: self={:?}, other_aligned={:?}",
                self_expected_dims,
                other_expected_dims
            )
            .into());
        }

        let axis_classes = if self.storage.axis_classes() == other_aligned.storage.axis_classes() {
            self.storage.axis_classes().to_vec()
        } else {
            Self::dense_axis_classes(self.indices.len())
        };

        let same_compact_layout = self.storage.payload_dims()
            == other_aligned.storage.payload_dims()
            && self.storage.payload_strides_vec() == other_aligned.storage.payload_strides_vec()
            && self.storage.axis_classes() == other_aligned.storage.axis_classes();
        if same_compact_layout
            && matches!(&self.storage, IdxTensorStorage::Materialized(_))
            && matches!(&other_aligned.storage, IdxTensorStorage::Materialized(_))
            && !self.tracks_grad()
            && !other_aligned.tracks_grad()
            && !a.tracks_grad()
            && !b.tracks_grad()
        {
            let lhs_storage = self.storage.materialize(self.indices.len())?;
            let rhs_storage = other_aligned
                .storage
                .materialize(other_aligned.indices.len())?;
            let combined = lhs_storage
                .axpby(
                    &a.to_backend_scalar(),
                    rhs_storage.as_ref(),
                    &b.to_backend_scalar(),
                )
                .map_err(|e| anyhow::anyhow!("storage axpby failed: {e}"))?;
            return Self::from_storage(self.indices.clone(), Arc::new(combined));
        }

        let lhs = self.scale(a)?;
        let rhs = other_aligned.scale(b)?;
        let combined = lhs
            .try_materialized_inner()?
            .add(rhs.try_materialized_inner()?)
            .map_err(|e| anyhow::anyhow!("tensor addition failed: {e}"))?;
        Self::from_inner_with_axis_classes(self.indices.clone(), combined, axis_classes)
            .map_err(IdxTensorError::from)
    }

    /// Scalar multiplication.
    ///
    /// Multiplies every element by `scalar`.
    ///
    /// # Errors
    /// Returns an error when the scalar coefficient is invalid for the tensor's
    /// scalar type (an invalid scalar dtype) or the backend reports a
    /// failure.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{AnyScalar, DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(3);
    /// let t = IdxTensor::from_dense(vec![i], vec![1.0, 2.0, 3.0]).unwrap();
    /// let scaled = t.scale(AnyScalar::new_real(2.0)).unwrap();
    /// assert_eq!(scaled.to_vec::<f64>().unwrap(), vec![2.0, 4.0, 6.0]);
    /// ```
    pub fn scale(&self, scalar: AnyScalar) -> std::result::Result<Self, IdxTensorError> {
        if matches!(
            &self.storage,
            IdxTensorStorage::Eager { .. }
                | IdxTensorStorage::Compact(_)
                | IdxTensorStorage::Materialized(_)
        ) {
            // Scale via the compact payload only. Materialized structured
            // storage is converted payload-coordinate by payload-coordinate
            // (never the logical domain) and returned as compact storage, so
            // scaling never touches unreferenced strided-gap backing entries.
            let storage = self.storage.scale_eager_payload(&scalar)?;
            return Ok(Self {
                indices: self.indices.clone(),
                storage,
                eager_cache: Self::empty_eager_cache(),
            });
        }

        let self_dtype = self.try_materialized_inner()?.dtype();
        let scalar_dtype = scalar.as_tensor()?.try_materialized_inner()?.dtype();
        if self_dtype != scalar_dtype {
            let target_dtype = Self::scale_target_dtype(self_dtype, scalar_dtype)?;
            let self_inner = self.try_materialized_inner()?;
            let self_inner = if self_inner.dtype() == target_dtype {
                self_inner.clone()
            } else {
                self_inner.cast(target_dtype)?
            };
            let scalar_inner = scalar.as_tensor()?.try_materialized_inner()?;
            let scalar_inner = if scalar_inner.dtype() == target_dtype {
                scalar_inner.clone()
            } else {
                scalar_inner.cast(target_dtype)?
            };
            let scaled = if self.indices.is_empty() {
                self_inner
                    .mul(&scalar_inner)
                    .map_err(|e| anyhow::anyhow!("scalar multiplication failed: {e}"))?
            } else {
                let subscripts = Self::scale_subscripts(self.indices.len())?;
                [&self_inner, &scalar_inner]
                    .einsum_subscripts(&subscripts)
                    .map_err(|e| anyhow::anyhow!("tensor scaling failed: {e}"))?
            };
            return Self::from_inner_with_axis_classes(
                self.indices.clone(),
                scaled,
                self.storage.axis_classes().to_vec(),
            )
            .map_err(IdxTensorError::from);
        }
        let scaled = if self.indices.is_empty() {
            self.try_materialized_inner()?
                .mul(scalar.as_tensor()?.try_materialized_inner()?)
                .map_err(|e| anyhow::anyhow!("scalar multiplication failed: {e}"))?
        } else {
            let subscripts = Self::scale_subscripts(self.indices.len())?;
            [
                self.try_materialized_inner()?,
                scalar.as_tensor()?.try_materialized_inner()?,
            ]
            .einsum_subscripts(&subscripts)
            .map_err(|e| anyhow::anyhow!("tensor scaling failed: {e}"))?
        };
        Self::from_inner_with_axis_classes(
            self.indices.clone(),
            scaled,
            self.storage.axis_classes().to_vec(),
        )
        .map_err(IdxTensorError::from)
    }

    /// Inner product (dot product) of two tensors.
    ///
    /// Computes `⟨self, other⟩ = Σ conj(self)_i * other_i`.
    ///
    /// # Errors
    /// Returns an error when the tensors have different index sets (an index-set
    /// mismatch) or the contraction reports a failure.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(3);
    /// let a = IdxTensor::from_dense(vec![i.clone()], vec![1.0, 2.0, 3.0]).unwrap();
    /// let b = IdxTensor::from_dense(vec![i.clone()], vec![4.0, 5.0, 6.0]).unwrap();
    ///
    /// // <a, b> = 1*4 + 2*5 + 3*6 = 32
    /// let ip = a.inner_product(&b).unwrap();
    /// assert!((ip.real() - 32.0).abs() < 1e-12);
    /// ```
    pub fn inner_product(&self, other: &Self) -> std::result::Result<AnyScalar, IdxTensorError> {
        if self.indices.len() == other.indices.len() {
            let self_set: HashSet<_> = self.indices.iter().collect();
            let other_set: HashSet<_> = other.indices.iter().collect();
            if self_set == other_set {
                let other_aligned = other.permute_indices(&self.indices)?;
                let result = super::contract::contract_pair_with_operand_options(
                    self,
                    &other_aligned,
                    PairwiseContractionOptions::new().with_lhs_conj(true),
                )?;
                return result.sum();
            }
        }

        // Contract self.conj() with other over all indices
        let result = super::contract::contract_pair_with_operand_options(
            self,
            other,
            PairwiseContractionOptions::new().with_lhs_conj(true),
        )?;
        // Result should be a scalar (no indices)
        result.sum()
    }
}

// ============================================================================
// Index Replacement Methods
// ============================================================================

impl IdxTensor {
    /// Replace an index in the tensor with a new index.
    ///
    /// This replaces every index equal to `old_index` (full index equality,
    /// including id, prime level, and tags) with `new_index`.
    /// The storage data is not modified, only the index metadata is changed.
    ///
    /// # Arguments
    /// * `old_index` - The index to replace (matched by full index equality)
    /// * `new_index` - The new index to use
    ///
    /// # Returns
    /// A new tensor with the index replaced. If no index matches `old_index`,
    /// returns a clone of the original tensor.
    ///
    /// # Errors
    /// Returns an error when the new index has an incompatible dimension
    /// (a shape mismatch: the replacement dimension must equal the original).
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let new_i = Index::new_dyn(2);  // Same dimension, different ID
    ///
    /// let indices = vec![i.clone(), j.clone()];
    /// let tensor: IdxTensor = IdxTensor::from_dense(indices, vec![0.0; 6]).unwrap();
    ///
    /// // Replace index i with new_i
    /// let replaced = tensor.replaceind(&i, &new_i).unwrap();
    /// assert_eq!(replaced.indices[0].id, new_i.id);
    /// assert_eq!(replaced.indices[1].id, j.id);
    /// ```
    pub fn replaceind(
        &self,
        old_index: &DynIndex,
        new_index: &DynIndex,
    ) -> std::result::Result<Self, IdxTensorError> {
        // Validate dimension match
        if old_index.dim() != new_index.dim() {
            return Err(IdxTensorError::ShapeMismatch {
                operation: "replaceind",
                expected: format!("dimension {}", old_index.dim()),
                actual: format!("dimension {}", new_index.dim()),
            });
        }

        let new_indices: Vec<_> = self
            .indices
            .iter()
            .map(|idx| {
                if *idx == *old_index {
                    new_index.clone()
                } else {
                    idx.clone()
                }
            })
            .collect();

        Ok(Self {
            indices: new_indices,
            storage: self.storage.clone(),
            eager_cache: Arc::clone(&self.eager_cache),
        })
    }

    /// Replace multiple indices in the tensor.
    ///
    /// This replaces each index in `old_indices` (matched by full index equality,
    /// including id, prime level, and tags) with the corresponding index in
    /// `new_indices`. The storage data is not modified.
    ///
    /// # Arguments
    /// * `old_indices` - The indices to replace (matched by full index equality)
    /// * `new_indices` - The new indices to use
    ///
    /// # Returns
    /// A new tensor with the indices replaced. Indices not found in `old_indices`
    /// are kept unchanged.
    ///
    /// # Errors
    /// Returns an error when `old_indices` and `new_indices` differ in length
    /// (a shape mismatch), or when any replacement index has an incompatible
    /// dimension (a shape mismatch: the replacement dimension must equal the
    /// original).
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let new_i = Index::new_dyn(2);
    /// let new_j = Index::new_dyn(3);
    ///
    /// let indices = vec![i.clone(), j.clone()];
    /// let tensor: IdxTensor = IdxTensor::from_dense(indices, vec![0.0; 6]).unwrap();
    ///
    /// // Replace both indices
    /// let replaced = tensor
    ///     .replace_indices(&[i.clone(), j.clone()], &[new_i.clone(), new_j.clone()])
    ///     .unwrap();
    /// assert_eq!(replaced.indices[0].id, new_i.id);
    /// assert_eq!(replaced.indices[1].id, new_j.id);
    /// ```
    pub fn replace_indices(
        &self,
        old_indices: &[DynIndex],
        new_indices: &[DynIndex],
    ) -> std::result::Result<Self, IdxTensorError> {
        if old_indices.len() != new_indices.len() {
            return Err(IdxTensorError::ShapeMismatch {
                operation: "replace_indices",
                expected: format!("{} indices", old_indices.len()),
                actual: format!("{} indices", new_indices.len()),
            });
        }

        // Validate dimension matches for all replacements
        for (old, new) in old_indices.iter().zip(new_indices.iter()) {
            if old.dim() != new.dim() {
                return Err(IdxTensorError::ShapeMismatch {
                    operation: "replace_indices",
                    expected: format!("dimension {}", old.dim()),
                    actual: format!("dimension {}", new.dim()),
                });
            }
        }

        // Build a map from old indices to new indices
        let replacement_map: std::collections::HashMap<_, _> =
            old_indices.iter().zip(new_indices.iter()).collect();

        let new_indices_vec: Vec<_> = self
            .indices
            .iter()
            .map(|idx| {
                if let Some(new_idx) = replacement_map.get(idx) {
                    (*new_idx).clone()
                } else {
                    idx.clone()
                }
            })
            .collect();

        Ok(Self {
            indices: new_indices_vec,
            storage: self.storage.clone(),
            eager_cache: Arc::clone(&self.eager_cache),
        })
    }
}

// ============================================================================
// Complex Conjugation
// ============================================================================

impl IdxTensor {
    /// Complex conjugate of all tensor elements.
    ///
    /// For real (`f32`/`f64`) tensors, returns a copy (conjugate of real is
    /// identity). For complex (`Complex32`/`Complex64`) tensors, conjugates
    /// each element.
    ///
    /// The indices and dimensions remain unchanged. If an eager backend cannot
    /// perform the conjugation, the failure is retained and reported by the
    /// next fallible materialization or AD-sensitive operation.
    ///
    /// This is inspired by the `conj` operation in ITensorMPS.jl.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    /// use num_complex::Complex64;
    ///
    /// let i = Index::new_dyn(2);
    /// let data = vec![Complex64::new(1.0, 2.0), Complex64::new(3.0, -4.0)];
    /// let tensor: IdxTensor = IdxTensor::from_dense(vec![i], data).unwrap();
    ///
    /// let conj_tensor = tensor.conj();
    /// assert_eq!(
    ///     conj_tensor.to_vec::<Complex64>().unwrap(),
    ///     vec![Complex64::new(1.0, -2.0), Complex64::new(3.0, 4.0)]
    /// );
    /// ```
    pub fn conj(&self) -> Self {
        self.conj_with(&conjugate_eager)
    }

    fn conj_with<F>(&self, conjugate: &F) -> Self
    where
        F: Fn(
            &EagerTensor,
        ) -> std::result::Result<
            EagerTensor,
            Arc<dyn std::error::Error + Send + Sync + 'static>,
        >,
    {
        // Conjugate tensor storage and map indices via IndexLike::conj(). For
        // default undirected indices, conj() is a no-op; this remains future-
        // proof for QSpace-compatible directed indices.
        let new_indices: Vec<DynIndex> = self.indices.iter().map(|idx| idx.conj()).collect();
        let mut storage = match self.storage.conjugate_with(conjugate) {
            Ok(storage) => storage,
            Err(error) => self.storage.clone().with_deferred_error(error),
        };
        let mut eager_cache = if storage.deferred_error().is_some() {
            Arc::clone(&self.eager_cache)
        } else {
            Self::empty_eager_cache()
        };

        if storage.deferred_error().is_none() {
            if let Some(inner) = self.eager_cache.get() {
                match conjugate(inner.as_ref()) {
                    Ok(conjugated) => eager_cache = Self::eager_cache_with(conjugated),
                    Err(source) => {
                        storage =
                            storage.with_deferred_error(TensorStorageError::Conjugation { source });
                        // Keep the original cache alive so a deferred failure
                        // retains its graph until a fallible consumer reports it.
                        eager_cache = Arc::clone(&self.eager_cache);
                    }
                }
            }
        }

        Self {
            indices: new_indices,
            storage,
            eager_cache,
        }
    }
}

#[derive(Debug, Default)]
struct Lassq {
    scale: f64,
    sumsq: f64,
    infinite: bool,
}

impl Lassq {
    fn add_component(&mut self, value: f64) {
        let value = value.abs();
        if value == 0.0 {
            return;
        }
        if value.is_infinite() {
            self.infinite = true;
            return;
        }
        if self.scale < value {
            if self.scale == 0.0 {
                self.sumsq = 1.0;
            } else {
                let ratio = self.scale / value;
                self.sumsq = 1.0 + self.sumsq * ratio * ratio;
            }
            self.scale = value;
        } else {
            let ratio = value / self.scale;
            self.sumsq += ratio * ratio;
        }
    }

    fn add_complex(&mut self, value: Complex64) {
        self.add_component(value.re);
        self.add_component(value.im);
    }

    fn add_scaled(&mut self, scale: f64, coefficient: f64) {
        if scale == 0.0 || coefficient == 0.0 {
            return;
        }
        if self.scale < scale {
            if self.scale == 0.0 {
                self.sumsq = coefficient * coefficient;
            } else {
                let ratio = self.scale / scale;
                self.sumsq = coefficient * coefficient + self.sumsq * ratio * ratio;
            }
            self.scale = scale;
        } else {
            let ratio = scale / self.scale * coefficient;
            self.sumsq += ratio * ratio;
        }
    }

    fn add_component_difference(&mut self, lhs: f64, rhs: f64) {
        if lhs == rhs {
            return;
        }
        let scale = lhs.abs().max(rhs.abs());
        if scale != 0.0 {
            self.add_scaled(scale, (lhs / scale - rhs / scale).abs());
        }
    }

    fn add_complex_difference(&mut self, lhs: Complex64, rhs: Complex64) {
        self.add_component_difference(lhs.re, rhs.re);
        self.add_component_difference(lhs.im, rhs.im);
    }

    fn is_zero(&self) -> bool {
        !self.infinite && self.scale == 0.0
    }

    fn norm(&self) -> f64 {
        if self.infinite {
            f64::INFINITY
        } else if self.scale == 0.0 {
            0.0
        } else {
            self.scale * self.sumsq.sqrt()
        }
    }

    fn norm_squared(&self) -> f64 {
        let norm = self.norm();
        norm * norm
    }

    fn log_norm(&self) -> f64 {
        if self.infinite {
            f64::INFINITY
        } else if self.scale == 0.0 {
            f64::NEG_INFINITY
        } else {
            self.scale.ln() + 0.5 * self.sumsq.ln()
        }
    }
}

// ============================================================================
// Norm Computation
// ============================================================================

impl IdxTensor {
    /// Compute the squared Frobenius norm of the tensor: ||T||² = Σ|T_ijk...|²
    ///
    /// For real tensors: sum of squares of all elements.
    /// For complex tensors: sum of `|z|²` over the compact payload.
    /// The reduction promotes source values to `f64` and uses a stable LASSQ
    /// accumulator, so it does not form a source-dtype `self * conj(self)`.
    ///
    /// # Errors
    /// Returns [`IdxTensorError`] when storage/materialization or scalar
    /// extraction fails, or when the input produces NaN. The result is
    /// accumulated from squared magnitudes, so it is never negative;
    /// positive infinity is preserved.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];  // 1² + 2² + ... + 6² = 91
    /// let tensor: IdxTensor = IdxTensor::from_dense(vec![i, j], data).unwrap();
    ///
    /// assert!((tensor.norm_squared().unwrap() - 91.0).abs() < 1e-10);
    /// ```
    pub fn norm_squared(&self) -> std::result::Result<f64, IdxTensorError> {
        let dtype = self
            .scalar_dtype()
            .map_err(IdxTensorError::scalar_extraction)?;
        if !matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64) {
            return Err(IdxTensorError::ScalarTypeMismatch {
                expected: "f32, f64, c32, or c64",
                actual: Self::dtype_name(dtype).to_string(),
            });
        }
        let (has_nan, _) = self
            .compact_nonfinite_flags()
            .map_err(IdxTensorError::materialization)?;
        if has_nan {
            return Err(IdxTensorError::NaNInput {
                operation: "norm_squared",
            });
        }

        let mut norm = Lassq::default();
        self.storage
            .for_each_payload_value(|value| norm.add_complex(value))
            .map_err(IdxTensorError::materialization)?;
        let value = norm.norm_squared();
        if value.is_nan() {
            return Err(IdxTensorError::NaNInput {
                operation: "norm_squared",
            });
        }
        Ok(value)
    }

    /// Compute the Frobenius norm of the tensor: ||T|| = sqrt(Σ|T_ijk...|²)
    ///
    /// # Errors
    /// Returns [`IdxTensorError`] when norm evaluation fails or when the
    /// input contains NaN. Positive infinity is preserved.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let data = vec![3.0, 4.0];  // sqrt(9 + 16) = 5
    /// let tensor: IdxTensor = IdxTensor::from_dense(vec![i], data).unwrap();
    ///
    /// assert!((tensor.norm().unwrap() - 5.0).abs() < 1e-10);
    /// ```
    pub fn norm(&self) -> std::result::Result<f64, IdxTensorError> {
        Ok(self.norm_squared()?.sqrt())
    }

    /// Maximum absolute value of all elements (L-infinity norm).
    ///
    /// # Errors
    /// Returns [`IdxTensorError`] when authoritative storage or eager
    /// materialization cannot be read, or when the input contains NaN.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(4);
    /// let t = IdxTensor::from_dense(vec![i], vec![-5.0, 1.0, 3.0, -2.0]).unwrap();
    /// assert!((t.maxabs().unwrap() - 5.0).abs() < 1e-12);
    /// ```
    pub fn maxabs(&self) -> std::result::Result<f64, IdxTensorError> {
        if let Some(error) = self.storage.deferred_error() {
            return Err(IdxTensorError::Storage {
                source: error.clone(),
            });
        }
        let dtype = self
            .storage
            .dtype()
            .ok_or_else(|| IdxTensorError::ScalarTypeMismatch {
                expected: "f32, f64, c32, or c64",
                actual: "unknown".to_string(),
            })?;
        if !matches!(dtype, DType::F32 | DType::F64 | DType::C32 | DType::C64) {
            return Err(IdxTensorError::ScalarTypeMismatch {
                expected: "f32, f64, c32, or c64",
                actual: Self::dtype_name(dtype).to_string(),
            });
        }
        let (has_nan, _) = self
            .compact_nonfinite_flags()
            .map_err(IdxTensorError::materialization)?;
        if has_nan {
            return Err(IdxTensorError::NaNInput {
                operation: "maxabs",
            });
        }
        let mut value = 0.0_f64;
        self.storage
            .for_each_payload_value(|scalar| {
                let magnitude = scalar.re.hypot(scalar.im);
                value = value.max(magnitude);
            })
            .map_err(IdxTensorError::materialization)?;
        Ok(value)
    }

    fn native_complex_payload_value_at(
        native: &EagerTensor,
        payload_coords: &[usize],
    ) -> Result<Complex64> {
        if !(payload_coords.len() == native.shape().len()) {
            return Err(anyhow::anyhow!(
                "payload coordinate rank {} does not match payload rank {}",
                payload_coords.len(),
                native.shape().len()
            ));
        };
        for (&coordinate, &dim) in payload_coords.iter().zip(native.shape().iter()) {
            if coordinate >= dim {
                return Err(anyhow::anyhow!(
                    "payload coordinate {coordinate} is out of bounds for dim {dim}"
                ));
            }
        }
        let value = native.value()?;
        match value.as_tensor_view() {
            TensorView::F32(view) => view
                .get(payload_coords)
                .copied()
                .map(|value| Complex64::new(f64::from(value), 0.0))
                .ok_or_else(|| anyhow::anyhow!("failed to read f32 payload value")),
            TensorView::F64(view) => view
                .get(payload_coords)
                .copied()
                .map(|value| Complex64::new(value, 0.0))
                .ok_or_else(|| anyhow::anyhow!("failed to read f64 payload value")),
            TensorView::C32(view) => view
                .get(payload_coords)
                .copied()
                .map(|value| Complex64::new(f64::from(value.re), f64::from(value.im)))
                .ok_or_else(|| anyhow::anyhow!("failed to read c32 payload value")),
            TensorView::C64(view) => view
                .get(payload_coords)
                .copied()
                .ok_or_else(|| anyhow::anyhow!("failed to read c64 payload value")),
            view => Err(anyhow::anyhow!(
                "unsupported payload dtype {:?}",
                view.dtype()
            )),
        }
    }

    fn native_sum_scalar(native: &EagerTensor) -> Result<AnyScalar> {
        let value = native.value()?;
        match native.dtype() {
            DType::F32 => Ok(AnyScalar::from_value(
                value.as_slice::<f32>()?.iter().copied().sum::<f32>(),
            )),
            DType::F64 => Ok(AnyScalar::from_value(
                value.as_slice::<f64>()?.iter().copied().sum::<f64>(),
            )),
            DType::C32 => Ok(AnyScalar::from_value(
                value
                    .as_slice::<Complex32>()?
                    .iter()
                    .copied()
                    .sum::<Complex32>(),
            )),
            DType::C64 => Ok(AnyScalar::from_value(
                value
                    .as_slice::<Complex64>()?
                    .iter()
                    .copied()
                    .sum::<Complex64>(),
            )),
            dtype => Err(anyhow::anyhow!("unsupported dtype {dtype:?}")),
        }
    }

    fn native_nonfinite_flags_typed<T: TensorElement>(
        native: &EagerTensor,
        classify: impl Fn(T) -> (bool, bool),
    ) -> Result<(bool, bool)> {
        let value = native.value()?;
        if let Ok(values) = value.as_slice::<T>() {
            return Ok(values.iter().copied().map(&classify).fold(
                (false, false),
                |(has_nan, has_infinity), (is_nan, is_infinite)| {
                    (has_nan | is_nan, has_infinity | is_infinite)
                },
            ));
        }
        drop(value);
        let value = IdxTensorStorage::materialize_eager_payload(native)?;
        Ok(value.as_slice::<T>()?.iter().copied().map(classify).fold(
            (false, false),
            |(has_nan, has_infinity), (is_nan, is_infinite)| {
                (has_nan | is_nan, has_infinity | is_infinite)
            },
        ))
    }

    fn native_nonfinite_flags(native: &EagerTensor) -> Result<(bool, bool)> {
        match native.dtype() {
            DType::F32 => Self::native_nonfinite_flags_typed(native, |value: f32| {
                (value.is_nan(), value.is_infinite())
            }),
            DType::F64 => Self::native_nonfinite_flags_typed(native, |value: f64| {
                (value.is_nan(), value.is_infinite())
            }),
            DType::C32 => Self::native_nonfinite_flags_typed(native, |value: Complex32| {
                (
                    value.re.is_nan() || value.im.is_nan(),
                    value.re.is_infinite() || value.im.is_infinite(),
                )
            }),
            DType::C64 => Self::native_nonfinite_flags_typed(native, |value: Complex64| {
                (
                    value.re.is_nan() || value.im.is_nan(),
                    value.re.is_infinite() || value.im.is_infinite(),
                )
            }),
            dtype => Err(anyhow::anyhow!("unsupported dtype {dtype:?}")),
        }
    }

    fn compact_nonfinite_flags(&self) -> Result<(bool, bool)> {
        self.storage.nonfinite_flags()
    }

    /// Element-wise subtraction with index alignment.
    ///
    /// This computes `self - other` using the same vector-space semantics as
    /// [`TensorVectorSpace`](crate::TensorVectorSpace).
    ///
    /// # Errors
    /// Returns an error when the tensors have different index sets (an index-set
    /// mismatch) or the arithmetic reports a failure.
    ///
    pub fn sub(&self, other: &Self) -> std::result::Result<Self, IdxTensorError> {
        self.axpby(AnyScalar::new_real(1.0), other, AnyScalar::new_real(-1.0))
    }

    /// Negate all elements.
    ///
    /// # Errors
    /// Returns an error when scalar multiplication fails for the tensor storage
    /// (a dtype mismatch) or the backend reports a failure.
    ///
    pub fn neg(&self) -> std::result::Result<Self, IdxTensorError> {
        self.scale(AnyScalar::new_real(-1.0))
    }

    /// Approximate equality check using Julia `isapprox`-style semantics.
    ///
    /// Values are aligned by index identity and streamed from each tensor's
    /// compact support. Exact zero-tolerance comparisons use exact scalar
    /// equality; nonzero tolerances use scaled sum-of-squares accumulation,
    /// avoiding logical-dense traversal and avoidable underflow/overflow.
    ///
    /// # Errors
    /// Returns [`IdxTensorError`] when tolerances are invalid, the index
    /// spaces cannot be aligned, storage cannot be read, or either input
    /// contains NaN.
    pub fn isapprox(
        &self,
        other: &Self,
        atol: f64,
        rtol: f64,
    ) -> std::result::Result<bool, IdxTensorError> {
        for (name, value) in [("atol", atol), ("rtol", rtol)] {
            if !value.is_finite() || value < 0.0 {
                return Err(IdxTensorError::InvalidTolerance { name, value });
            }
        }
        if self.indices.len() != other.indices.len() {
            return Err(IdxTensorError::ShapeMismatch {
                operation: "isapprox",
                expected: format!("indices {:?}", self.indices),
                actual: format!("indices {:?}", other.indices),
            });
        }

        let other_axis_by_index = other
            .indices
            .iter()
            .cloned()
            .enumerate()
            .map(|(axis, index)| (index, axis))
            .collect::<HashMap<_, _>>();
        let other_positions = self
            .indices
            .iter()
            .map(|index| {
                other_axis_by_index.get(index).copied().ok_or_else(|| {
                    IdxTensorError::ShapeMismatch {
                        operation: "isapprox",
                        expected: format!("indices {:?}", self.indices),
                        actual: format!("indices {:?}", other.indices),
                    }
                })
            })
            .collect::<std::result::Result<Vec<_>, _>>()?;
        let self_dims = self.dims();
        let other_dims = other.dims();
        for (axis, &other_axis) in other_positions.iter().enumerate() {
            if self_dims[axis] != other_dims[other_axis] {
                return Err(IdxTensorError::ShapeMismatch {
                    operation: "isapprox",
                    expected: format!("dims {:?}", self_dims),
                    actual: format!("dims {:?}", other_dims),
                });
            }
        }
        for tensor in [self, other] {
            if tensor
                .compact_nonfinite_flags()
                .map_err(IdxTensorError::materialization)?
                .0
            {
                return Err(IdxTensorError::NaNInput {
                    operation: "isapprox",
                });
            }
        }

        let exact = atol == 0.0 && rtol == 0.0;
        let lhs_payload_dims = self.storage.payload_dims().to_vec();
        let rhs_payload_dims = other.storage.payload_dims().to_vec();
        let lhs_payload_len =
            checked_product(&lhs_payload_dims).map_err(IdxTensorError::materialization)?;
        let rhs_payload_len =
            checked_product(&rhs_payload_dims).map_err(IdxTensorError::materialization)?;
        let lhs_axis_classes = self.storage.axis_classes();
        let rhs_axis_classes = other.storage.axis_classes();
        let self_to_other = other_positions.clone();
        let mut other_to_self = vec![0usize; self_to_other.len()];
        for (self_axis, &other_axis) in self_to_other.iter().enumerate() {
            other_to_self[other_axis] = self_axis;
        }

        let mut diff = Lassq::default();
        let mut lhs_norm = Lassq::default();
        let mut rhs_norm = Lassq::default();
        let mut compare = |lhs: Complex64, rhs: Complex64| -> bool {
            if exact {
                return lhs == rhs;
            }
            let lhs_infinite = lhs.re.is_infinite() || lhs.im.is_infinite();
            let rhs_infinite = rhs.re.is_infinite() || rhs.im.is_infinite();
            if lhs_infinite || rhs_infinite {
                return lhs == rhs;
            }
            lhs_norm.add_complex(lhs);
            rhs_norm.add_complex(rhs);
            diff.add_complex_difference(lhs, rhs);
            true
        };

        // Each compact payload coordinate identifies exactly one logical
        // support point. Map it through the aligned logical axes instead of
        // traversing structural zeros in the logical tensor.
        let mut lhs_coords = vec![0usize; lhs_payload_dims.len()];
        let mut rhs_from_lhs = vec![0usize; rhs_payload_dims.len()];
        let mut rhs_seen = vec![false; rhs_payload_dims.len()];
        for _ in 0..lhs_payload_len {
            let lhs = self
                .storage
                .payload_value_at(&lhs_coords)
                .map_err(IdxTensorError::materialization)?;
            if map_payload_support_coordinate(
                lhs_axis_classes,
                rhs_axis_classes,
                &self_to_other,
                &lhs_coords,
                &mut rhs_from_lhs,
                &mut rhs_seen,
            ) {
                let rhs = other
                    .storage
                    .payload_value_at(&rhs_from_lhs)
                    .map_err(IdxTensorError::materialization)?;
                if !compare(lhs, rhs) {
                    return Ok(false);
                }
            } else if !compare(lhs, Complex64::new(0.0, 0.0)) {
                return Ok(false);
            }

            increment_col_major_coordinate(&mut lhs_coords, &lhs_payload_dims);
        }

        // The reverse pass accounts for support points that exist only in the
        // right tensor. Overlap points were compared in the first pass and are
        // therefore skipped here without a payload-sized visited set.
        let mut rhs_coords = vec![0usize; rhs_payload_dims.len()];
        let mut lhs_from_rhs = vec![0usize; lhs_payload_dims.len()];
        let mut lhs_seen = vec![false; lhs_payload_dims.len()];
        for _ in 0..rhs_payload_len {
            let rhs = other
                .storage
                .payload_value_at(&rhs_coords)
                .map_err(IdxTensorError::materialization)?;
            if !map_payload_support_coordinate(
                rhs_axis_classes,
                lhs_axis_classes,
                &other_to_self,
                &rhs_coords,
                &mut lhs_from_rhs,
                &mut lhs_seen,
            ) && !compare(Complex64::new(0.0, 0.0), rhs)
            {
                return Ok(false);
            }
            increment_col_major_coordinate(&mut rhs_coords, &rhs_payload_dims);
        }

        if exact {
            return Ok(true);
        }
        let absolute_ok = if diff.is_zero() {
            true
        } else if atol == 0.0 || diff.infinite {
            false
        } else {
            diff.log_norm() <= atol.ln()
        };
        let relative_ok = if rtol == 0.0 || diff.is_zero() {
            diff.is_zero()
        } else if diff.infinite {
            lhs_norm.infinite || rhs_norm.infinite
        } else if lhs_norm.infinite || rhs_norm.infinite {
            true
        } else {
            let reference_log = lhs_norm.log_norm().max(rhs_norm.log_norm());
            diff.log_norm() <= rtol.ln() + reference_log
        };
        Ok(absolute_ok || relative_ok)
    }

    /// Create a diagonal Kronecker-delta tensor for one input/output index pair.
    ///
    /// # Errors
    /// Returns an error when the two indices have different dimensions (an
    /// index shape mismatch).
    ///
    pub fn diagonal(
        input_index: &DynIndex,
        output_index: &DynIndex,
    ) -> std::result::Result<Self, IdxTensorError> {
        <Self as TensorConstructionLike>::diagonal(input_index, output_index)
    }

    /// Create a product of Kronecker-delta tensors for paired index lists.
    ///
    /// # Errors
    /// Returns an error if the index lists have different lengths or paired
    /// dimensions do not match.
    pub fn delta(
        input_indices: &[DynIndex],
        output_indices: &[DynIndex],
    ) -> std::result::Result<Self, IdxTensorError> {
        <Self as TensorConstructionLike>::delta(input_indices, output_indices)
    }

    /// Create a scalar tensor equal to one.
    ///
    /// # Errors
    /// Returns an error when dense scalar construction fails for the element type
    /// (an invalid scalar dtype or a construction failure).
    ///
    pub fn scalar_one() -> std::result::Result<Self, IdxTensorError> {
        <Self as TensorConstructionLike>::scalar_one()
    }

    /// Create a tensor filled with ones over the given indices.
    ///
    /// # Errors
    /// Returns an error when the tensor size overflows (an overflow failure) or
    /// dense construction fails.
    ///
    pub fn ones(indices: &[DynIndex]) -> std::result::Result<Self, IdxTensorError> {
        <Self as TensorConstructionLike>::ones(indices)
    }

    /// Create a one-hot tensor with value one at the specified index positions.
    ///
    /// # Errors
    /// Returns an error when any coordinate is outside its index dimension (an
    /// out of bounds failure).
    ///
    pub fn onehot(index_vals: &[(DynIndex, usize)]) -> std::result::Result<Self, IdxTensorError> {
        <Self as TensorConstructionLike>::onehot(index_vals)
    }

    /// Keep one coordinate along an index while retaining that index axis.
    ///
    /// This is the differentiable masking counterpart to [`Self::select_indices`].
    /// It selects the requested slice, forms a one-hot tensor over the removed
    /// axis in the source dtype, and takes an explicit tensor product to restore
    /// the original index order. The implementation stays in the tensor backend,
    /// so structured storage and reverse-mode metadata are preserved whenever
    /// the backend can represent the operation.
    ///
    /// # Arguments
    ///
    /// * `index` - Existing tensor index to mask.
    /// * `position` - Zero-based coordinate to keep; all other coordinates become
    ///
    ///   zero.
    ///
    /// # Errors
    /// Returns an error when the coordinate is outside the index dimension (an
    /// out of bounds failure) or the mask construction fails.
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(vec![i.clone()], vec![3.0_f64, 4.0]).unwrap();
    /// let masked = tensor.mask_index(&i, 1).unwrap();
    ///
    /// assert_eq!(masked.indices(), &[i]);
    /// assert_eq!(masked.to_vec::<f64>().unwrap(), vec![0.0, 4.0]);
    /// assert!(IdxTensor::from_dense(
    ///     vec![DynIndex::new_dyn(2)],
    ///     vec![1.0_f64, 2.0],
    /// )
    /// .unwrap()
    /// .mask_index(&DynIndex::new_dyn(2), 0)
    /// .is_err());
    /// ```
    pub fn mask_index(
        &self,
        index: &DynIndex,
        position: usize,
    ) -> std::result::Result<Self, IdxTensorError> {
        if !(self.indices.iter().any(|candidate| candidate == index)) {
            return Err(anyhow::anyhow!("mask_index: index is not present in tensor").into());
        };
        if !(position < index.dim()) {
            return Err(anyhow::anyhow!(
                "mask_index: position {position} is out of range for dimension {}",
                index.dim()
            )
            .into());
        };

        // Retaining the shared index turns contraction into a backend-level
        // elementwise product instead of materializing a host mask. Construct
        // the constant mask in the input dtype so f32/c32 values and AD graphs
        // are not promoted or detached.
        let mask = match self.scalar_dtype()? {
            DType::F32 => Self::from_dense(
                vec![index.clone()],
                (0..index.dim())
                    .map(|value| if value == position { 1.0_f32 } else { 0.0 })
                    .collect(),
            ),
            DType::F64 => Self::from_dense(
                vec![index.clone()],
                (0..index.dim())
                    .map(|value| if value == position { 1.0_f64 } else { 0.0 })
                    .collect(),
            ),
            DType::C32 => Self::from_dense(
                vec![index.clone()],
                (0..index.dim())
                    .map(|value| {
                        if value == position {
                            num_complex::Complex32::new(1.0, 0.0)
                        } else {
                            num_complex::Complex32::new(0.0, 0.0)
                        }
                    })
                    .collect(),
            ),
            DType::C64 => Self::from_dense(
                vec![index.clone()],
                (0..index.dim())
                    .map(|value| {
                        if value == position {
                            Complex64::new(1.0, 0.0)
                        } else {
                            Complex64::new(0.0, 0.0)
                        }
                    })
                    .collect(),
            ),
            dtype => {
                return Err(anyhow::anyhow!("mask_index does not support dtype {dtype:?}").into())
            }
        }?;
        super::contract::contract_pair_with_options(
            self,
            &mask,
            super::contract::ContractionOptions::new()
                .with_retain_indices(std::slice::from_ref(index)),
        )
    }

    /// Compute the relative distance between two tensors.
    ///
    /// Returns `||A - B|| / ||A||` (Frobenius norm).
    /// If `||A|| = 0`, returns `||B||` instead to avoid division by zero.
    ///
    /// This is the ITensor-style distance function useful for comparing tensors.
    ///
    /// # Arguments
    /// * `other` - The other tensor to compare with
    ///
    /// # Errors
    /// Returns [`IdxTensorError`] when either norm contains NaN, or when
    /// scaling and subtracting the tensors fails.
    ///
    /// # Returns
    /// The relative distance as a f64 value.
    ///
    /// # Note
    /// The indices of both tensors must be permutable to each other.
    /// The result tensor (A - B) uses the index ordering from self.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let data_a = vec![1.0, 0.0];
    /// let data_b = vec![1.0, 0.0];  // Same tensor
    /// let tensor_a: IdxTensor = IdxTensor::from_dense(vec![i.clone()], data_a).unwrap();
    /// let tensor_b: IdxTensor = IdxTensor::from_dense(vec![i.clone()], data_b).unwrap();
    ///
    /// assert!(tensor_a.distance(&tensor_b).unwrap() < 1e-10);  // Zero distance
    /// ```
    pub fn distance(&self, other: &Self) -> std::result::Result<f64, IdxTensorError> {
        let norm_self = self.norm()?;

        // Compute A - B = A + (-1) * B
        let neg_other = other.scale(AnyScalar::new_real(-1.0))?;
        let diff = self.add(&neg_other)?;
        let norm_diff = diff.norm()?;

        if norm_self > 0.0 {
            Ok(norm_diff / norm_self)
        } else {
            Ok(norm_diff)
        }
    }
}

impl std::fmt::Debug for IdxTensor {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("IdxTensor")
            .field("indices", &self.indices)
            .field("dims", &self.dims())
            .field("is_diag", &self.is_diag())
            .finish()
    }
}

/// Create a diagonal tensor with dynamic rank from diagonal data.
/// # Arguments
/// * `indices` - The indices for the tensor (all must have the same dimension)
/// * `diag_data` - The diagonal elements (length must equal the dimension of indices)
///
/// The returned tensor preserves compact diagonal payload metadata; use
/// [`IdxTensor::is_diag`] or [`IdxTensor::storage`] to inspect that
/// representation.
///
/// # Errors
/// Returns an error when the index dimensions are unequal (a dimension
/// mismatch) or the diagonal construction fails.
/// # Panics
/// Panics if indices have different dimensions, or if diag_data length doesn't match.
/// # Examples
/// ```
/// use tensor4all_core::{DynIndex, diag_idx_tensor};
/// let i = DynIndex::new_dyn(3);
/// let j = DynIndex::new_dyn(3);
/// let t = diag_idx_tensor(vec![i, j], vec![1.0, 2.0, 3.0]).unwrap();
/// assert_eq!(t.dims(), vec![3, 3]);
/// assert!(t.is_diag());
/// ```
pub fn diag_idx_tensor(
    indices: Vec<DynIndex>,
    diag_data: Vec<f64>,
) -> std::result::Result<IdxTensor, IdxTensorError> {
    IdxTensor::from_diag(indices, diag_data)
}

#[allow(clippy::type_complexity)]
pub(crate) type UnfoldSplitInnerResult = (
    EagerTensor,
    usize,
    usize,
    usize,
    Vec<DynIndex>,
    Vec<DynIndex>,
);

/// Unfold a tensor into a matrix by splitting indices into left and right groups.
/// This function validates the split, permutes the tensor so that left indices
/// come first, and returns a rank-2 native tenferro tensor along with metadata.
/// # Arguments
/// * `t` - Input tensor
/// * `left_inds` - Indices to place on the left (row) side of the matrix
/// # Returns
/// A tuple `(matrix_tensor, left_len, m, n, left_indices, right_indices)` where:
/// - `matrix_tensor` is a rank-2 `tenferro::Tensor` with shape `[m, n]`
/// - `left_len` is the number of left indices
/// - `m` is the product of left index dimensions
/// - `n` is the product of right index dimensions
/// - `left_indices` is the vector of left indices (cloned)
/// - `right_indices` is the vector of right indices (cloned)
/// # Errors
///
/// Returns an error when the tensor rank is less than 2 (a rank mismatch),
/// when `left_inds` is empty or contains all indices (an invalid split), when
/// `left_inds` contains indices not present in the tensor (a missing-index
/// failure) or duplicates, or when the native reshape fails (a backend
/// failure).
/// # Examples
/// ```
/// use tensor4all_core::{DynIndex, IdxTensor, unfold_split};
/// let i = DynIndex::new_dyn(2);
/// let j = DynIndex::new_dyn(3);
/// // 2x3 dense tensor with data [1..6]
/// let t = IdxTensor::from_dense(
///     vec![i.clone(), j.clone()],
///     vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
/// ).unwrap();
/// let (matrix, left_len, m, n, left_indices, right_indices) =
///     unfold_split(&t, &[i]).unwrap();
/// assert_eq!(left_len, 1);
/// assert_eq!(m, 2);
/// assert_eq!(n, 3);
/// assert_eq!(left_indices.len(), 1);
/// assert_eq!(right_indices.len(), 1);
/// ```
#[allow(clippy::type_complexity)]
pub fn unfold_split(
    t: &IdxTensor,
    left_inds: &[DynIndex],
) -> std::result::Result<
    (
        NativeTensor,
        usize,
        usize,
        usize,
        Vec<DynIndex>,
        Vec<DynIndex>,
    ),
    IdxTensorError,
> {
    let (matrix_inner, left_len, m, n, left_indices, right_indices) =
        unfold_split_inner(t, left_inds)?;

    Ok((
        matrix_inner.duplicate_value()?,
        left_len,
        m,
        n,
        left_indices,
        right_indices,
    ))
}

pub(crate) fn unfold_split_inner(
    t: &IdxTensor,
    left_inds: &[DynIndex],
) -> Result<UnfoldSplitInnerResult> {
    let rank = t.indices.len();

    // Validate rank
    if !(rank >= 2) {
        return Err(anyhow::anyhow!(
            "Tensor must have rank >= 2, got rank {}",
            rank
        ));
    };

    let left_len = left_inds.len();

    // Validate split: must be a proper subset
    if !(left_len > 0 && left_len < rank) {
        return Err(anyhow::anyhow!("Left indices must be a non-empty proper subset of tensor indices (0 < left_len < rank), got left_len={}, rank={}",
        left_len,
        rank));
    };

    // Validate that all left_inds are in the tensor and there are no duplicates
    let tensor_set: HashSet<_> = t.indices.iter().collect();
    let mut left_set = HashSet::new();

    for left_idx in left_inds {
        if !(tensor_set.contains(left_idx)) {
            return Err(anyhow::anyhow!("Index in left_inds not found in tensor"));
        };
        if !(left_set.insert(left_idx)) {
            return Err(anyhow::anyhow!("Duplicate index in left_inds"));
        };
    }

    // Build right_inds: all indices not in left_inds, in original order
    let mut right_inds = Vec::new();
    for idx in &t.indices {
        if !left_set.contains(idx) {
            right_inds.push(idx.clone());
        }
    }

    // Build new_indices: left_inds first, then right_inds
    let mut new_indices = Vec::with_capacity(rank);
    new_indices.extend_from_slice(left_inds);
    new_indices.extend_from_slice(&right_inds);

    // Permute tensor to have left indices first, then right indices
    let unfolded = t.permute_indices(&new_indices)?;

    // Compute matrix dimensions
    let unfolded_dims = unfolded.dims();
    let m = checked_product(&unfolded_dims[..left_len])?;
    let n = checked_product(&unfolded_dims[left_len..])?;

    let matrix_tensor = unfolded.try_materialized_inner()?.reshape(&[m, n])?;

    Ok((
        matrix_tensor,
        left_len,
        m,
        n,
        left_inds.to_vec(),
        right_inds,
    ))
}

// ============================================================================
// TensorIndex implementation for IdxTensor
// ============================================================================

use crate::tensor_index::TensorIndex;

impl TensorIndex for IdxTensor {
    type Index = DynIndex;
    type Error = IdxTensorError;

    fn external_indices(&self) -> Vec<DynIndex> {
        // For IdxTensor, all indices are external.
        self.indices.clone()
    }

    fn num_external_indices(&self) -> usize {
        self.indices.len()
    }

    fn replaceind(
        &self,
        old_index: &DynIndex,
        new_index: &DynIndex,
    ) -> std::result::Result<Self, Self::Error> {
        // Delegate to the inherent method.
        IdxTensor::replaceind(self, old_index, new_index)
    }

    fn replace_indices(
        &self,
        old_indices: &[DynIndex],
        new_indices: &[DynIndex],
    ) -> std::result::Result<Self, Self::Error> {
        // Delegate to the inherent method.
        IdxTensor::replace_indices(self, old_indices, new_indices)
    }
}

// ============================================================================
// TensorLike implementation for IdxTensor
// ============================================================================

use crate::tensor_like::{
    FactorizeError, FactorizeOptions, FactorizeResult, IncrementalQrState, TensorConstructionLike,
    TensorContractionLike, TensorFactorizationLike, TensorVectorSpace,
};

impl TensorVectorSpace for IdxTensor {
    fn norm_squared(&self) -> std::result::Result<f64, Self::Error> {
        IdxTensor::norm_squared(self)
    }

    fn maxabs(&self) -> std::result::Result<f64, Self::Error> {
        IdxTensor::maxabs(self)
    }

    fn isapprox(
        &self,
        other: &Self,
        atol: f64,
        rtol: f64,
    ) -> std::result::Result<bool, Self::Error> {
        IdxTensor::isapprox(self, other, atol, rtol)
    }

    fn axpby(
        &self,
        a: crate::AnyScalar,
        other: &Self,
        b: crate::AnyScalar,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::axpby(self, a, other, b)
    }

    fn scale(&self, scalar: crate::AnyScalar) -> std::result::Result<Self, Self::Error> {
        IdxTensor::scale(self, scalar)
    }

    fn scale_in(
        &self,
        factor: f64,
        context: &ExecutionContext,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::scale_in(self, factor, context)
    }

    fn norm_in(&self, context: &ExecutionContext) -> std::result::Result<f64, Self::Error> {
        IdxTensor::norm_in(self, context)
    }

    fn inner_product(&self, other: &Self) -> std::result::Result<crate::AnyScalar, Self::Error> {
        IdxTensor::inner_product(self, other)
    }
}

impl TensorFactorizationLike for IdxTensor {
    fn factorize(
        &self,
        left_inds: &[DynIndex],
        options: &FactorizeOptions,
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        crate::factorize::factorize(self, left_inds, options)
    }

    fn factorize_in(
        &self,
        left_inds: &[DynIndex],
        options: &FactorizeOptions,
        context: &ExecutionContext,
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        crate::factorize::factorize_in(self, left_inds, options, context)
    }

    fn factorize_auto(
        &self,
        left_inds: &[DynIndex],
        options: &FactorizeOptions,
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        crate::factorize::factorize_auto(self, left_inds, options)
    }

    fn factorize_full_rank(
        &self,
        left_inds: &[DynIndex],
        alg: crate::FactorizeAlg,
        canonical: crate::Canonical,
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        crate::factorize::factorize_full_rank(self, left_inds, alg, canonical)
    }

    fn factorize_full_rank_in(
        &self,
        left_inds: &[DynIndex],
        alg: crate::FactorizeAlg,
        canonical: crate::Canonical,
        context: &ExecutionContext,
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        crate::factorize::factorize_full_rank_in(self, left_inds, alg, canonical, context)
    }

    fn factorize_probe_columns_incremental(
        previous: Option<&FactorizeResult<Self>>,
        all_columns: &[&Self],
        appended_columns: &[&Self],
        left_inds: &[DynIndex],
    ) -> std::result::Result<FactorizeResult<Self>, FactorizeError> {
        factorize_probe_columns_incremental(previous, all_columns, appended_columns, left_inds)
    }

    fn factorize_probe_batch_incremental(
        previous: Option<&FactorizeResult<IdxTensor>>,
        batch_tensor: &IdxTensor,
        batch_index: &DynIndex,
        left_inds: &[DynIndex],
    ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
        factorize_probe_batch_incremental_impl(previous, batch_tensor, batch_index, left_inds)
    }

    fn factorize_probe_batch_incremental_in(
        previous: Option<&FactorizeResult<IdxTensor>>,
        batch_tensor: &IdxTensor,
        batch_index: &DynIndex,
        left_inds: &[DynIndex],
        context: &ExecutionContext,
    ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
        if let Some(previous) = previous {
            previous
                .left
                .validate_context(context)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            previous
                .right
                .validate_context(context)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        }
        batch_tensor
            .validate_context(context)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        // Compatibility boundary: the process-global CPU context keeps the
        // historical host `IncrementalQr` matrix path (bitwise-identical
        // numerics for legacy callers). Every other explicit context uses the
        // scoped recompute path with no host round-trip.
        if context.is_global_default_cpu() {
            return factorize_probe_batch_incremental_impl(
                previous,
                batch_tensor,
                batch_index,
                left_inds,
            );
        }
        Self::resident_probe_batch_qr(previous, batch_tensor, batch_index, left_inds, context)
    }

    fn src_error_estimate(
        &self,
    ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
        if self.is_cuda_resident() {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "CUDA-resident QR factor requires src_error_estimate_in with the owning execution context"
            )));
        }
        let indices = self.indices();
        if indices.len() != 2 {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC QR factor must have rank 2, got {}",
                indices.len()
            )));
        }
        let nrows = indices[0].dim();
        let ncols = indices[1].dim();
        let result = if self.is_f64() {
            let matrix = Matrix::from_col_major_vec(
                nrows,
                ncols,
                self.to_vec::<f64>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );
            backend_src_error_estimate(&matrix)
        } else if self.is_c64() {
            let matrix = Matrix::from_col_major_vec(
                nrows,
                ncols,
                self.to_vec::<Complex64>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );
            backend_src_error_estimate(&matrix)
        } else {
            return Err(FactorizeError::UnsupportedStorage(
                "SRC adaptive estimation currently supports f64 and Complex64 QR factors",
            ));
        };
        result.map_err(|error| FactorizeError::ComputationError(anyhow::anyhow!(error)))
    }

    fn src_error_estimate_in(
        &self,
        context: &ExecutionContext,
    ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
        self.validate_context(context)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        #[cfg(feature = "tenferro-cuda")]
        if matches!(context, ExecutionContext::Cuda(_)) {
            return self.resident_src_error_estimate(context);
        }
        if context.is_global_default_cpu() {
            return TensorFactorizationLike::src_error_estimate(self);
        }
        self.host_src_error_estimate_general()
    }
}

impl IdxTensor {
    /// Evaluate the SRC estimator for an explicit CPU context's restored factor.
    fn host_src_error_estimate_general(
        &self,
    ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
        let indices = self.indices();
        if indices.len() != 2 {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC QR factor must have rank 2, got {}",
                indices.len()
            )));
        }
        let nrows = indices[0].dim();
        let ncols = indices[1].dim();
        let result = if self.is_f64() {
            let matrix = Matrix::from_col_major_vec(
                nrows,
                ncols,
                self.to_vec::<f64>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );
            backend_src_error_estimate_general(&matrix)
        } else if self.is_c64() {
            let matrix = Matrix::from_col_major_vec(
                nrows,
                ncols,
                self.to_vec::<Complex64>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );
            backend_src_error_estimate_general(&matrix)
        } else {
            return Err(FactorizeError::UnsupportedStorage(
                "SRC adaptive estimation currently supports f64 and Complex64 QR factors",
            ));
        };
        result.map_err(|error| FactorizeError::ComputationError(anyhow::anyhow!(error)))
    }

    /// Resident batch-native probe factorization by full RRQR recompute.
    ///
    /// Reconstructs the full sketch block `[Q_prev R_prev, A_new]` with resident
    /// matmul/concatenation and runs column-pivoted rank-revealing QR in the
    /// owning runtime. Only the scalar rank metadata is explicitly read back;
    /// matrix payloads and permutation metadata remain resident. The returned
    /// right factor is projected back into original sketch-column order, so
    /// `left * right` reconstructs the unpermuted sketch. Since that right
    /// factor is not generally triangular, the adaptive estimator uses a
    /// general solve of the small square sketch factor.
    fn resident_probe_batch_qr(
        previous: Option<&FactorizeResult<IdxTensor>>,
        batch_tensor: &IdxTensor,
        batch_index: &DynIndex,
        left_inds: &[DynIndex],
        context: &ExecutionContext,
    ) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
        use crate::defaults::idx_tensor::unfold_split_inner;

        // The host path consumes `batch_index` as the appended-axis label;
        // reject a mismatched label before any arithmetic, as the device path
        // otherwise never observes it.
        if !batch_tensor.indices().contains(batch_index) {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "resident probe batch is missing its batch index"
            )));
        }
        let (appended_inner, _, m, b, _, _) = unfold_split_inner(batch_tensor, left_inds)
            .context("resident probe batch unfold failed")
            .map_err(FactorizeError::ComputationError)?;
        if !matches!(appended_inner.dtype(), DType::F64 | DType::C64) {
            return Err(FactorizeError::UnsupportedStorage(
                "incremental SRC factorization currently supports f64 and Complex64 tensors",
            ));
        }
        if b == 0 {
            if let Some(previous) = previous {
                return Ok(previous.clone());
            }
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "incremental SRC factorization requires at least one column"
            )));
        }
        let full_inner = if let Some(previous) = previous {
            let (q_inner, _, previous_rows, _, _, _) =
                unfold_split_inner(&previous.left, left_inds)
                    .context("resident previous Q unfold failed")
                    .map_err(FactorizeError::ComputationError)?;
            if previous_rows != m {
                return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                    "previous and appended sketch row dimensions differ: {previous_rows} vs {m}"
                )));
            }
            let (r_inner, _, _, _, _, _) =
                unfold_split_inner(&previous.right, std::slice::from_ref(&previous.bond_index))
                    .context("resident previous R unfold failed")
                    .map_err(FactorizeError::ComputationError)?;
            let reconstructed = q_inner.matmul(&r_inner).map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident sketch reconstruction failed"),
                )
            })?;
            EagerTensor::concatenate(&[&reconstructed, &appended_inner], 1).map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident sketch concatenation failed"),
                )
            })?
        } else {
            appended_inner
        };
        let total_width = *full_inner.shape().get(1).ok_or_else(|| {
            FactorizeError::ComputationError(anyhow::anyhow!(
                "resident sketch factor is not rank-2"
            ))
        })?;
        let rank_tolerance = 32.0 * f64::EPSILON * m.max(total_width) as f64;
        let decomposition = full_inner
            .rank_revealing_qr(RankRevealingQrOptions::default().rtol(rank_tolerance))
            .map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident probe batch RRQR failed"),
                )
            })?;
        let rank = Self::read_resident_rank(&decomposition.rank, context)?;
        let maximum_rank = m.min(total_width);
        if rank > maximum_rank {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "resident RRQR returned rank {rank}, exceeding maximum {maximum_rank}"
            )));
        }
        if let Some(previous) = previous {
            if rank < previous.rank {
                return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                    "resident SRC RRQR rank decreased from {} to {rank}",
                    previous.rank
                )));
            }
        }
        let q_full = decomposition
            .q
            .slice_axis(1, 0..rank)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        // RRQR factors the pivoted sketch A[:, permutation]. Restore the
        // original sketch-column order without reading permutation metadata:
        // Q_rank^H A computes the equivalent right factor entirely resident.
        let q_adjoint = q_full
            .transpose(&[1, 0])
            .and_then(|transposed| transposed.conj())
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let r_full = q_adjoint.matmul(&full_inner).map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("resident RRQR right-factor restoration failed"),
            )
        })?;
        let cap = DynIndex::new_bond(rank)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let batch = DynIndex::new_link(total_width)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let mut q_indices = left_inds.to_vec();
        q_indices.push(cap.clone());
        let q_dims: Vec<usize> = q_indices.iter().map(|index| index.dim()).collect();
        let left = Self::from_inner(
            q_indices,
            q_full.reshape(&q_dims).map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident probe batch Q reshape failed"),
                )
            })?,
        )
        .map_err(FactorizeError::ComputationError)?;
        let right = Self::from_inner(
            vec![cap.clone(), batch],
            r_full.reshape(&[rank, total_width]).map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident probe batch R reshape failed"),
                )
            })?,
        )
        .map_err(FactorizeError::ComputationError)?;
        Ok(FactorizeResult::new(left, right, cap, None, rank))
    }

    /// Device-side Appendix-C estimator: solve and reduce on-device, read back
    /// `ncols + 1` decision scalars.
    ///
    /// Singularity and non-finiteness surface as typed solve/validation errors;
    /// unlike the host path there is no separate diagonal pre-check, but every
    /// invalid input still fails instead of producing an estimate.
    #[cfg(feature = "tenferro-cuda")]
    fn resident_src_error_estimate(
        &self,
        context: &ExecutionContext,
    ) -> std::result::Result<tensor4all_tensorbackend::SrcErrorEstimate, FactorizeError> {
        use tensor4all_tensorbackend::SrcErrorEstimate;

        let indices = self.indices();
        if indices.len() != 2 {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC QR factor must have rank 2, got {}",
                indices.len()
            )));
        }
        let nrows = indices[0].dim();
        let ncols = indices[1].dim();
        if nrows != ncols {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC estimator requires a square R, got {nrows}x{ncols}"
            )));
        }
        if nrows == 0 {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC estimator requires a non-empty R"
            )));
        }
        let inner = self.cuda_eager_inner().ok_or_else(|| {
            FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC estimator requires a resident eager QR factor"
            ))
        })?;
        let ExecutionContext::Cuda(cuda) = context else {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "SRC device estimator requires a CUDA execution context"
            )));
        };
        // RRQR restores original sketch-column order, so the factor is square
        // but not generally triangular. Solve its adjoint as a general system.
        let adjoint = inner
            .transpose(&[1, 0])
            .and_then(|transposed| transposed.conj())
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let identity = Self::resident_identity(cuda, nrows, inner.dtype())?;
        // Solve F^H X = I for X = F^{-†}.
        let solved = adjoint.solve(&identity).map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("SRC inverse-adjoint general solve failed"),
            )
        })?;
        let column_sums = solved
            .abs()
            .and_then(|magnitudes| magnitudes.reduce_sum_squares(&[0]))
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let total = inner
            .abs()
            .and_then(|magnitudes| magnitudes.reduce_sum_squares(&[0, 1]))
            .and_then(|scalar| scalar.reshape(&[1]))
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
        let column_norms_sq = Self::read_resident_vector(&column_sums, ncols, context)?;
        let total_norm_sq = Self::read_resident_vector(&total, 1, context)?;
        Self::src_estimate_from_decision_scalars(&column_norms_sq, total_norm_sq[0], ncols)
            .map(|(error, norm)| SrcErrorEstimate { error, norm })
            .map_err(FactorizeError::ComputationError)
    }

    /// Read the scalar integer rank returned by RRQR through its owning context.
    fn read_resident_rank(
        rank: &EagerTensor,
        context: &ExecutionContext,
    ) -> std::result::Result<usize, FactorizeError> {
        if !rank.shape().is_empty() {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "resident RRQR rank metadata has shape {:?}, expected []",
                rank.shape()
            )));
        }
        let resident = rank.to_tensor().map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("resident RRQR rank materialization failed"),
            )
        })?;
        #[cfg(feature = "tenferro-cuda")]
        let host = match context {
            ExecutionContext::Cpu(_) => resident,
            ExecutionContext::Cuda(cuda) => cuda.download(&resident).map_err(|error| {
                FactorizeError::ComputationError(
                    anyhow::Error::new(error).context("resident RRQR rank download failed"),
                )
            })?,
        };
        #[cfg(not(feature = "tenferro-cuda"))]
        let host = match context {
            ExecutionContext::Cpu(_) => resident,
        };
        let values = host.as_slice::<i64>().map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("resident RRQR rank decode failed"),
            )
        })?;
        let value = *values.first().ok_or_else(|| {
            FactorizeError::ComputationError(anyhow::anyhow!(
                "resident RRQR rank metadata is empty"
            ))
        })?;
        usize::try_from(value).map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("resident RRQR rank is negative"),
            )
        })
    }

    /// Build an `n x n` identity in the operand's owning runtime.
    ///
    /// Host-originated construction data takes the same explicit upload path
    /// as [`IdxTensor::from_dense_in`]; the transfer is part of the estimator's
    /// decision payload, not hidden arithmetic.
    #[cfg(feature = "tenferro-cuda")]
    fn resident_identity(
        cuda: &Arc<tensor4all_tensorbackend::CudaExecutionContext>,
        n: usize,
        dtype: DType,
    ) -> std::result::Result<EagerTensor, FactorizeError> {
        let native = match dtype {
            DType::F64 => {
                let mut data = vec![0.0_f64; n * n];
                for diagonal in 0..n {
                    data[diagonal + diagonal * n] = 1.0;
                }
                NativeTensor::from_vec_col_major(vec![n, n], data)
            }
            DType::C64 => {
                let mut data = vec![Complex64::new(0.0, 0.0); n * n];
                for diagonal in 0..n {
                    data[diagonal + diagonal * n] = Complex64::new(1.0, 0.0);
                }
                NativeTensor::from_vec_col_major(vec![n, n], data)
            }
            _other => {
                return Err(FactorizeError::UnsupportedStorage(
                    "SRC adaptive estimation currently supports f64 and Complex64 QR factors",
                ));
            }
        }
        .map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("SRC estimator identity construction failed"),
            )
        })?;
        let uploaded = cuda.upload_cuda(&native).map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("SRC estimator identity upload failed"),
            )
        })?;
        let runtime = cuda.eager_runtime().map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("SRC estimator eager runtime failed"),
            )
        })?;
        // `validate_context` already pinned the factor to this exact CUDA
        // context, so the uploaded identity lands in the factor's runtime.
        EagerTensor::from_tensor_in(uploaded, runtime).map_err(|error| {
            FactorizeError::ComputationError(
                anyhow::Error::new(error).context("SRC estimator identity wrapping failed"),
            )
        })
    }

    /// Wrap a resident decision vector and read it back through `context`.
    #[cfg(feature = "tenferro-cuda")]
    fn read_resident_vector(
        eager: &EagerTensor,
        len: usize,
        context: &ExecutionContext,
    ) -> std::result::Result<Vec<f64>, FactorizeError> {
        if eager.shape() != [len] {
            return Err(FactorizeError::ComputationError(anyhow::anyhow!(
                "device SRC decision vector has shape {:?}, expected [{len}]",
                eager.shape()
            )));
        }
        let decision = Self::from_inner(vec![DynIndex::new_dyn(len)], eager.clone())
            .map_err(FactorizeError::ComputationError)?;
        decision
            .read_decision_data(context)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))
    }

    /// Scalar Appendix-C formula from the `ncols + 1` decision values.
    ///
    /// Mirrors `src_error_estimate_from_inverse_adjoint`: column norms must be
    /// finite and nonzero, and both estimates must be finite.
    #[cfg(feature = "tenferro-cuda")]
    fn src_estimate_from_decision_scalars(
        column_norms_sq: &[f64],
        total_norm_sq: f64,
        ncols: usize,
    ) -> std::result::Result<(f64, f64), anyhow::Error> {
        if !total_norm_sq.is_finite() {
            return Err(anyhow::anyhow!(
                "SRC estimator requires finite entries in R"
            ));
        }
        let mut inverse_column_error_sq = 0.0_f64;
        for (col, &column_norm_sq) in column_norms_sq.iter().enumerate() {
            if !column_norm_sq.is_finite() || column_norm_sq == 0.0 {
                return Err(anyhow::anyhow!(
                    "SRC inverse-adjoint solve returned an invalid column norm at column {col}"
                ));
            }
            inverse_column_error_sq += 1.0 / column_norm_sq;
        }
        let sketch_width = ncols as f64;
        let error_sq = inverse_column_error_sq / sketch_width;
        let norm_estimate_sq = total_norm_sq / sketch_width;
        if !error_sq.is_finite() || !norm_estimate_sq.is_finite() {
            return Err(anyhow::anyhow!(
                "SRC estimator produced a non-finite estimate"
            ));
        }
        Ok((error_sq.sqrt(), norm_estimate_sq.sqrt()))
    }
}

// Provenance: the public seam mirrors the append contract in
// `chriscamano/RandomMPOMPS/code/tensornetwork/incrementalqr.py::IncrementalQR.append`
// and `incrementalqr.cpp::add_cols`; the tensor/index bridge and private state
// are tensor4all-specific and are labelled `[AI-Supplied]` in the audit.
fn factorize_probe_columns_incremental(
    previous: Option<&FactorizeResult<IdxTensor>>,
    all_columns: &[&IdxTensor],
    appended_columns: &[&IdxTensor],
    left_inds: &[DynIndex],
) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
    let first = all_columns.first().ok_or_else(|| {
        FactorizeError::ComputationError(anyhow::anyhow!(
            "incremental SRC factorization requires at least one column"
        ))
    })?;
    if first.is_f64() {
        incremental_probe_factorize_typed::<f64>(previous, all_columns, appended_columns, left_inds)
    } else if first.is_c64() {
        incremental_probe_factorize_typed::<Complex64>(
            previous,
            all_columns,
            appended_columns,
            left_inds,
        )
    } else {
        Err(FactorizeError::UnsupportedStorage(
            "incremental SRC factorization currently supports f64 and Complex64 tensors",
        ))
    }
}

trait IncrementalQrStateScalar: IncrementalQrScalar + TensorElement {
    fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>>;

    fn store(state: IncrementalQr<Self>) -> IncrementalQrState;
}

impl IncrementalQrStateScalar for f64 {
    fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>> {
        match previous.incremental_qr_state()? {
            IncrementalQrState::F64(state) => Some(state.clone()),
            IncrementalQrState::C64(_) => None,
        }
    }

    fn store(state: IncrementalQr<Self>) -> IncrementalQrState {
        IncrementalQrState::F64(state)
    }
}

impl IncrementalQrStateScalar for Complex64 {
    fn resume(previous: &FactorizeResult<IdxTensor>) -> Option<IncrementalQr<Self>> {
        match previous.incremental_qr_state()? {
            IncrementalQrState::F64(_) => None,
            IncrementalQrState::C64(state) => Some(state.clone()),
        }
    }

    fn store(state: IncrementalQr<Self>) -> IncrementalQrState {
        IncrementalQrState::C64(state)
    }
}

fn incremental_probe_factorize_from_matrices<S>(
    previous: Option<&FactorizeResult<IdxTensor>>,
    appended_is_empty: bool,
    left_inds: &[DynIndex],
    initial_matrix: impl FnOnce() -> std::result::Result<Matrix<S>, FactorizeError>,
    appended_matrix: impl FnOnce() -> std::result::Result<Matrix<S>, FactorizeError>,
) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
where
    S: IncrementalQrStateScalar,
{
    let (mut state, previous_left, previous_cap, previous_rank) = if let Some(previous) = previous {
        if appended_is_empty {
            return Ok(previous.clone());
        }
        let previous_rank = previous.rank;
        let resumed = S::resume(previous);
        let reused_previous_q = resumed.is_some();
        let state = if let Some(state) = resumed {
            state
        } else {
            let q = previous.left.clone();
            let cap = previous.bond_index.clone();
            let mut q_indices = left_inds.to_vec();
            q_indices.push(cap.clone());
            let q = q
                .permute_indices(&q_indices)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            let q_dims = q_indices.iter().map(IndexLike::dim).collect::<Vec<_>>();
            let q_rows = q_dims[..q_dims.len() - 1]
                .iter()
                .try_fold(1usize, |size, &dim| size.checked_mul(dim))
                .ok_or_else(|| {
                    FactorizeError::ComputationError(anyhow::anyhow!(
                        "incremental SRC Q row dimension overflows usize"
                    ))
                })?;
            let q_matrix = Matrix::from_col_major_vec(
                q_rows,
                cap.dim(),
                q.to_vec::<S>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );

            let batch = previous.right.indices().last().cloned().ok_or_else(|| {
                FactorizeError::ComputationError(anyhow::anyhow!(
                    "incremental SRC R factor has no batch index"
                ))
            })?;
            let r = previous
                .right
                .permute_indices(&[cap.clone(), batch.clone()])
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            let r_matrix = Matrix::from_col_major_vec(
                cap.dim(),
                batch.dim(),
                r.to_vec::<S>()
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            );
            IncrementalQr::from_factors(q_matrix, r_matrix)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
        };
        (
            state,
            if reused_previous_q {
                Some(previous.left.clone())
            } else {
                None
            },
            if reused_previous_q {
                Some(previous.bond_index.clone())
            } else {
                None
            },
            previous_rank,
        )
    } else {
        let initial_matrix = initial_matrix()?;
        (
            IncrementalQr::new(initial_matrix)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
            None,
            None,
            0,
        )
    };
    if !appended_is_empty && previous.is_some() {
        let appended_matrix = appended_matrix()?;
        state
            .append(&appended_matrix)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    }

    let rank = state.rank();
    if previous_left.is_some() && rank < previous_rank {
        return Err(FactorizeError::ComputationError(anyhow::anyhow!(
            "incremental SRC QR rank decreased from {previous_rank} to {rank}"
        )));
    }
    let r = state.r();
    let cap = if rank == previous_rank {
        match previous_cap.as_ref() {
            Some(cap) => cap.clone(),
            None => DynIndex::new_bond(rank)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
        }
    } else {
        DynIndex::new_bond(rank)
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
    };
    let left = if let Some(previous_left) = previous_left {
        if rank == previous_rank {
            previous_left
        } else {
            let appended_rank = rank - previous_rank;
            let appended_cap = DynIndex::new_bond(appended_rank)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            let appended_q = state
                .q_columns(previous_rank, appended_rank)
                .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            let mut appended_indices = left_inds.to_vec();
            appended_indices.push(appended_cap.clone());
            let appended_left =
                IdxTensor::from_dense(appended_indices, appended_q.as_col_major_slice().to_vec())
                    .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
            IdxTensor::concatenate_along_new_index(
                &[&previous_left, &appended_left],
                &[
                    previous_cap.ok_or_else(|| {
                        FactorizeError::ComputationError(anyhow::anyhow!(
                            "incremental SRC previous QR state has no bond index"
                        ))
                    })?,
                    appended_cap,
                ],
                cap.clone(),
            )
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
        }
    } else {
        let q = state.q();
        let mut q_indices = left_inds.to_vec();
        q_indices.push(cap.clone());
        IdxTensor::from_dense(q_indices, q.as_col_major_slice().to_vec())
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?
    };
    let batch = DynIndex::new_link(r.ncols())
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    let right = IdxTensor::from_dense(vec![cap.clone(), batch], r.as_col_major_slice().to_vec())
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    Ok(FactorizeResult::new(left, right, cap, None, rank)
        .with_incremental_qr_state(S::store(state)))
}

fn incremental_probe_factorize_typed<S>(
    previous: Option<&FactorizeResult<IdxTensor>>,
    all_columns: &[&IdxTensor],
    appended_columns: &[&IdxTensor],
    left_inds: &[DynIndex],
) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
where
    S: IncrementalQrStateScalar,
{
    incremental_probe_factorize_from_matrices::<S>(
        previous,
        appended_columns.is_empty(),
        left_inds,
        || probe_columns_matrix::<S>(all_columns, left_inds),
        || probe_columns_matrix::<S>(appended_columns, left_inds),
    )
}

fn incremental_probe_factorize_batch_typed<S>(
    previous: Option<&FactorizeResult<IdxTensor>>,
    batch_tensor: &IdxTensor,
    batch_index: &DynIndex,
    left_inds: &[DynIndex],
) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError>
where
    S: IncrementalQrStateScalar,
{
    incremental_probe_factorize_from_matrices::<S>(
        previous,
        batch_index.dim() == 0,
        left_inds,
        || probe_batch_matrix::<S>(batch_tensor, batch_index, left_inds),
        || probe_batch_matrix::<S>(batch_tensor, batch_index, left_inds),
    )
}

fn factorize_probe_batch_incremental_impl(
    previous: Option<&FactorizeResult<IdxTensor>>,
    batch_tensor: &IdxTensor,
    batch_index: &DynIndex,
    left_inds: &[DynIndex],
) -> std::result::Result<FactorizeResult<IdxTensor>, FactorizeError> {
    if batch_tensor.is_f64() {
        incremental_probe_factorize_batch_typed::<f64>(
            previous,
            batch_tensor,
            batch_index,
            left_inds,
        )
    } else if batch_tensor.is_c64() {
        incremental_probe_factorize_batch_typed::<Complex64>(
            previous,
            batch_tensor,
            batch_index,
            left_inds,
        )
    } else {
        Err(FactorizeError::UnsupportedStorage(
            "incremental SRC factorization currently supports f64 and Complex64 tensors",
        ))
    }
}

fn probe_columns_matrix<S>(
    columns: &[&IdxTensor],
    left_inds: &[DynIndex],
) -> std::result::Result<Matrix<S>, FactorizeError>
where
    S: TensorElement,
{
    let batch = DynIndex::new_link(columns.len())
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    let stacked = IdxTensor::stack_along_new_index(columns, batch.clone(), -1)
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    let mut ordered_indices = left_inds.to_vec();
    ordered_indices.push(batch);
    let ordered = stacked
        .permute_indices(&ordered_indices)
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    let nrows = left_inds
        .iter()
        .try_fold(1usize, |size, index| size.checked_mul(index.dim()))
        .ok_or_else(|| {
            FactorizeError::ComputationError(anyhow::anyhow!(
                "incremental SRC sketch row dimension overflows usize"
            ))
        })?;
    Ok(Matrix::from_col_major_vec(
        nrows,
        columns.len(),
        ordered
            .to_vec::<S>()
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
    ))
}

/// Like [`probe_columns_matrix`], but for a tensor that already carries a
/// batch axis (`batch_index`) instead of being split into separate column
/// tensors — skips the `stack_along_new_index` call `probe_columns_matrix`
/// needs to re-assemble one.
fn probe_batch_matrix<S>(
    batch_tensor: &IdxTensor,
    batch_index: &DynIndex,
    left_inds: &[DynIndex],
) -> std::result::Result<Matrix<S>, FactorizeError>
where
    S: TensorElement,
{
    let mut ordered_indices = left_inds.to_vec();
    ordered_indices.push(batch_index.clone());
    let ordered = batch_tensor
        .permute_indices(&ordered_indices)
        .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?;
    let nrows = left_inds
        .iter()
        .try_fold(1usize, |size, index| size.checked_mul(index.dim()))
        .ok_or_else(|| {
            FactorizeError::ComputationError(anyhow::anyhow!(
                "incremental SRC sketch row dimension overflows usize"
            ))
        })?;
    Ok(Matrix::from_col_major_vec(
        nrows,
        batch_index.dim(),
        ordered
            .to_vec::<S>()
            .map_err(|error| FactorizeError::ComputationError(anyhow::Error::new(error)))?,
    ))
}

impl TensorContractionLike for IdxTensor {
    fn conj(&self) -> Self {
        // Delegate to the inherent method (complex conjugate for dense tensors)
        IdxTensor::conj(self)
    }

    fn direct_sum(
        &self,
        other: &Self,
        pairs: &[(DynIndex, DynIndex)],
    ) -> std::result::Result<crate::tensor_like::DirectSumResult<Self>, Self::Error> {
        let (tensor, new_indices) = crate::direct_sum::direct_sum(self, other, pairs)?;
        Ok(crate::tensor_like::DirectSumResult {
            tensor,
            new_indices,
        })
    }

    fn outer_product(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
        super::contract::outer_product(self, other)
    }

    fn permuteinds(&self, new_order: &[DynIndex]) -> std::result::Result<Self, Self::Error> {
        // Delegate to the inherent method
        IdxTensor::permute_indices(self, new_order)
    }

    fn fuse_indices(
        &self,
        old_indices: &[DynIndex],
        new_index: DynIndex,
        order: LinearizationOrder,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::fuse_indices(self, old_indices, new_index, order)
    }

    fn contract(tensors: &[&Self]) -> std::result::Result<Self, Self::Error> {
        super::contract::contract(tensors)
    }

    fn contract_retaining_indices(
        tensors: &[&Self],
        retained_indices: &[DynIndex],
    ) -> std::result::Result<Self, Self::Error> {
        if tensors.len() == 2 {
            return tensors[0]
                .try_contract_pairwise_retaining(tensors[1], retained_indices)
                .map_err(IdxTensorError::from);
        }
        let options = ContractionOptions::new().with_retain_indices(retained_indices);
        super::contract::contract_with_options(tensors, options)
    }

    fn contract_pair(&self, other: &Self) -> std::result::Result<Self, Self::Error> {
        super::contract::contract_pair(self, other)
    }
}

impl TensorConstructionLike for IdxTensor {
    fn select_indices(
        &self,
        selected_indices: &[DynIndex],
        positions: &[usize],
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::select_indices(self, selected_indices, positions)
    }

    fn diagonal(
        input_index: &DynIndex,
        output_index: &DynIndex,
    ) -> std::result::Result<Self, Self::Error> {
        let dim = input_index.dim();
        if dim != output_index.dim() {
            return Err(anyhow::anyhow!(
                "Dimension mismatch: input index has dim {}, output has dim {}",
                dim,
                output_index.dim(),
            )
            .into());
        }

        IdxTensor::from_diag(
            vec![input_index.clone(), output_index.clone()],
            vec![1.0_f64; dim],
        )
    }

    fn scalar_one() -> std::result::Result<Self, Self::Error> {
        IdxTensor::from_dense(vec![], vec![1.0_f64])
    }

    fn ones(indices: &[DynIndex]) -> std::result::Result<Self, Self::Error> {
        if indices.is_empty() {
            return <Self as TensorConstructionLike>::scalar_one();
        }
        let dims: Vec<usize> = indices.iter().map(|idx| idx.size()).collect();
        let total_size = checked_total_size(&dims)?;
        IdxTensor::from_dense(indices.to_vec(), vec![1.0_f64; total_size])
    }

    fn ones_in(
        context: &tensor4all_tensorbackend::ExecutionContext,
        indices: &[DynIndex],
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::ones_in(context, indices)
    }

    fn validate_context(
        &self,
        context: &tensor4all_tensorbackend::ExecutionContext,
    ) -> std::result::Result<(), Self::Error> {
        IdxTensor::validate_context(self, context)
    }

    fn from_dense_any(
        indices: Vec<DynIndex>,
        data: Vec<AnyScalar>,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::from_dense_any(indices, data)
    }

    fn from_dense<T: TensorElement + Into<AnyScalar>>(
        indices: Vec<DynIndex>,
        data: Vec<T>,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::from_dense(indices, data)
    }

    fn from_dense_in<T: TensorElement + Into<AnyScalar>>(
        context: &tensor4all_tensorbackend::ExecutionContext,
        indices: Vec<DynIndex>,
        data: Vec<T>,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::from_dense_in(context, indices, data)
    }

    fn stack_along_new_index(
        tensors: &[&Self],
        new_index: DynIndex,
        axis: isize,
    ) -> std::result::Result<Self, Self::Error> {
        IdxTensor::stack_along_new_index(tensors, new_index, axis)
    }

    fn concatenate_along_new_index(
        tensors: &[&Self],
        source_indices: &[DynIndex],
        new_index: DynIndex,
    ) -> std::result::Result<Self, Self::Error> {
        let first = tensors
            .first()
            .ok_or_else(|| anyhow::anyhow!("concatenate requires at least one tensor"))?;
        if tensors.len() != source_indices.len() {
            return Err(anyhow::anyhow!(
                "concatenate tensor count {} does not match source-index count {}",
                tensors.len(),
                source_indices.len()
            )
            .into());
        }
        let first_axis = first
            .indices
            .iter()
            .position(|index| index == &source_indices[0])
            .ok_or_else(|| anyhow::anyhow!("concatenate source index is not present"))?;
        let mut total_dim = 0usize;
        let first_indices = first.indices();
        let mut inners = Vec::with_capacity(tensors.len());
        for (tensor, source) in tensors.iter().zip(source_indices) {
            let axis = tensor
                .indices
                .iter()
                .position(|index| index == source)
                .ok_or_else(|| anyhow::anyhow!("concatenate source index is not present"))?;
            if axis != first_axis {
                return Err(
                    anyhow::anyhow!("concatenate source axes must have the same position").into(),
                );
            }
            if tensor
                .indices()
                .iter()
                .zip(first_indices)
                .enumerate()
                .any(|(position, (actual, expected))| position != first_axis && actual != expected)
            {
                return Err(anyhow::anyhow!(
                    "concatenate tensors must match away from the source axis"
                )
                .into());
            }
            total_dim = total_dim
                .checked_add(source.dim())
                .ok_or_else(|| anyhow::anyhow!("concatenate dimension overflow"))?;
            tensor.ensure_shape_packing_preserves_ad("concatenate_along_new_index")?;
            inners.push(tensor.try_materialized_inner()?);
        }
        if new_index.dim() != total_dim {
            return Err(anyhow::anyhow!(
                "concatenate output dimension {} does not match source dimension sum {}",
                new_index.dim(),
                total_dim
            )
            .into());
        }
        let concatenated = EagerTensor::concatenate(&inners, first_axis)?;
        let mut indices = first_indices.to_vec();
        indices[first_axis] = new_index;
        Self::from_inner(indices, concatenated).map_err(IdxTensorError::from)
    }

    fn onehot(index_vals: &[(DynIndex, usize)]) -> std::result::Result<Self, Self::Error> {
        if index_vals.is_empty() {
            return <Self as TensorConstructionLike>::scalar_one();
        }
        let indices: Vec<DynIndex> = index_vals.iter().map(|(idx, _)| idx.clone()).collect();
        let vals: Vec<usize> = index_vals.iter().map(|(_, v)| *v).collect();
        let dims: Vec<usize> = indices.iter().map(|idx| idx.size()).collect();

        for (k, (&v, &d)) in vals.iter().zip(dims.iter()).enumerate() {
            if v >= d {
                return Err(anyhow::anyhow!(
                    "onehot: value {} at position {} is >= dimension {}",
                    v,
                    k,
                    d
                )
                .into());
            }
        }

        let total_size = checked_total_size(&dims).map_err(Self::Error::from)?;
        let mut data = vec![0.0_f64; total_size];

        let offset = column_major_offset(&dims, &vals).map_err(Self::Error::from)?;
        data[offset] = 1.0;

        Self::from_dense(indices, data)
    }

    // delta() uses the default implementation via diagonal() and outer_product()
}

fn checked_total_size(dims: &[usize]) -> Result<usize> {
    dims.iter().try_fold(1_usize, |acc, &d| {
        if d == 0 {
            return Err(anyhow::anyhow!("invalid dimension 0"));
        }
        acc.checked_mul(d)
            .ok_or_else(|| anyhow::anyhow!("tensor size overflow"))
    })
}

fn column_major_offset(dims: &[usize], vals: &[usize]) -> Result<usize> {
    if dims.len() != vals.len() {
        return Err(anyhow::anyhow!(
            "column_major_offset: dims.len() != vals.len()"
        ));
    }
    checked_total_size(dims)?;

    let mut offset = 0usize;
    let mut stride = 1usize;
    for (k, (&v, &d)) in vals.iter().zip(dims.iter()).enumerate() {
        if d == 0 {
            return Err(anyhow::anyhow!("invalid dimension 0 at position {}", k));
        }
        if v >= d {
            return Err(anyhow::anyhow!(
                "column_major_offset: value {} at position {} is >= dimension {}",
                v,
                k,
                d
            ));
        }
        let term = v
            .checked_mul(stride)
            .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
        offset = offset
            .checked_add(term)
            .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
        stride = stride
            .checked_mul(d)
            .ok_or_else(|| anyhow::anyhow!("column_major_offset: overflow"))?;
    }
    Ok(offset)
}

// ============================================================================
// High-level API for tensor construction (avoids direct Storage access)
// ============================================================================

impl IdxTensor {
    fn any_scalar_payload_to_complex(data: Vec<AnyScalar>) -> Vec<Complex64> {
        data.into_iter()
            .map(|value| {
                value
                    .as_c64()
                    .unwrap_or_else(|| Complex64::new(value.real(), 0.0))
            })
            .collect()
    }

    fn any_scalar_payload_to_real(data: Vec<AnyScalar>) -> Vec<f64> {
        data.into_iter().map(|value| value.real()).collect()
    }

    fn validate_dense_payload_len(data_len: usize, dims: &[usize]) -> Result<()> {
        let expected_len = checked_total_size(dims)?;
        if !(data_len == expected_len) {
            return Err(anyhow::anyhow!(
                "dense payload length {} does not match dims {:?} (expected {})",
                data_len,
                dims,
                expected_len
            ));
        };
        Ok(())
    }

    fn validate_diag_payload_len(data_len: usize, dims: &[usize]) -> Result<()> {
        if !(!dims.is_empty()) {
            return Err(anyhow::anyhow!(
                "diagonal tensor construction requires at least one index"
            ));
        };
        Self::validate_diag_dims(dims)?;
        if !(data_len == dims[0]) {
            return Err(anyhow::anyhow!(
                "diagonal payload length {} does not match diagonal dimension {}",
                data_len,
                dims[0]
            ));
        };
        Ok(())
    }

    /// Create a tensor from dense data with explicit indices.
    ///
    /// This is the recommended high-level API for creating tensors from raw data.
    /// It avoids direct access to `Storage` internals.
    ///
    /// # Type Parameters
    /// * `T` - Scalar type (`f32`, `f64`, `Complex32`, or `Complex64`)
    ///
    /// # Arguments
    /// * `indices` - Vector of indices for the tensor
    /// * `data` - Tensor data in column-major order
    ///
    /// # Errors
    /// Returns an error when the data length does not match the index dimension
    /// product (a shape mismatch).
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
    /// let tensor: IdxTensor = IdxTensor::from_dense(vec![i, j], data).unwrap();
    /// assert_eq!(tensor.dims(), vec![2, 3]);
    /// ```
    pub fn from_dense<T: TensorElement>(
        indices: Vec<DynIndex>,
        data: Vec<T>,
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims = Self::expected_dims_from_indices(&indices);
        Self::validate_indices(&indices)?;
        Self::validate_dense_payload_len(data.len(), &dims)?;
        let native = dense_native_tensor_from_col_major_owned(data, &dims)?;
        Self::from_native(indices, native).map_err(IdxTensorError::from)
    }

    /// Create a tensor from dense payload data provided as [`AnyScalar`] values.
    ///
    /// This is the preferred public API when the caller only knows the scalar
    /// type at runtime.
    ///
    /// # Errors
    /// Returns an error when the payload length does not match the index dimension
    /// product (a shape mismatch) or a scalar conversion fails.
    /// # Examples
    /// ```
    /// use tensor4all_core::{AnyScalar, IdxTensor};
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(2);
    /// let tensor = IdxTensor::from_dense_any(
    ///     vec![i, j],
    ///     vec![
    ///         AnyScalar::new_real(1.0),
    ///         AnyScalar::new_complex(0.0, 1.0),
    ///         AnyScalar::new_real(2.0),
    ///         AnyScalar::new_real(3.0),
    ///     ],
    /// ).unwrap();
    ///
    /// assert!(tensor.is_complex());
    /// assert_eq!(tensor.dims(), vec![2, 2]);
    /// ```
    pub fn from_dense_any(
        indices: Vec<DynIndex>,
        data: Vec<AnyScalar>,
    ) -> std::result::Result<Self, IdxTensorError> {
        if data.iter().any(AnyScalar::is_complex) {
            Self::from_dense(indices, Self::any_scalar_payload_to_complex(data))
        } else {
            Self::from_dense(indices, Self::any_scalar_payload_to_real(data))
        }
    }

    /// Create a diagonal tensor from diagonal payload data with explicit indices.
    ///
    /// All indices must have the same dimension, and `data.len()` must equal
    /// that dimension. The resulting tensor has nonzero entries only on
    /// the multi-index diagonal (`T[i,i,...,i] = data[i]`).
    ///
    /// The returned tensor preserves diagonal metadata; use
    /// [`IdxTensor::is_diag`] or [`IdxTensor::storage_kind`] to inspect
    /// that representation. `f32` and `Complex32` values remain eager and are
    /// never promoted into compact `f64`/`Complex64` storage.
    ///
    /// # Errors
    /// Returns an error when the index dimensions are unequal or the payload
    /// length does not match the diagonal dimension (a shape mismatch).
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(3);
    /// let j = DynIndex::new_dyn(3);
    /// let diag = IdxTensor::from_diag(vec![i, j], vec![1.0, 2.0, 3.0]).unwrap();
    /// assert!(diag.is_diag());
    ///
    /// let data = diag.to_vec::<f64>().unwrap();
    /// // 3x3 identity-like: [1,0,0, 0,2,0, 0,0,3] in column-major
    /// assert!((data[0] - 1.0).abs() < 1e-12);
    /// assert!((data[4] - 2.0).abs() < 1e-12);
    /// assert!((data[8] - 3.0).abs() < 1e-12);
    /// assert!((data[1]).abs() < 1e-12);  // off-diagonal is zero
    /// ```
    pub fn from_diag<T: TensorElement>(
        indices: Vec<DynIndex>,
        data: Vec<T>,
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims = Self::expected_dims_from_indices(&indices);
        Self::validate_indices(&indices)?;
        Self::validate_diag_payload_len(data.len(), &dims)?;
        let native = diag_native_tensor_from_col_major(&data, dims.len())?;
        Self::from_native_with_axis_classes(indices, native, Self::diag_axis_classes(dims.len()))
            .map_err(IdxTensorError::from)
    }

    /// Create a diagonal tensor from diagonal payload data provided as
    /// [`AnyScalar`] values.
    ///
    /// This is the preferred public API when the caller only knows the scalar
    /// type at runtime.
    ///
    /// # Errors
    /// Returns an error when the index dimensions are unequal or the payload
    /// length does not match (a shape mismatch), or a scalar conversion
    /// fails.
    /// # Examples
    /// ```
    /// use tensor4all_core::{AnyScalar, IdxTensor};
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(2);
    /// let tensor = IdxTensor::from_diag_any(
    ///     vec![i, j],
    ///     vec![AnyScalar::new_real(1.0), AnyScalar::new_complex(2.0, -1.0)],
    /// ).unwrap();
    ///
    /// assert!(tensor.is_complex());
    /// assert_eq!(tensor.dims(), vec![2, 2]);
    /// ```
    pub fn from_diag_any(
        indices: Vec<DynIndex>,
        data: Vec<AnyScalar>,
    ) -> std::result::Result<Self, IdxTensorError> {
        if data.iter().any(AnyScalar::is_complex) {
            Self::from_diag(indices, Self::any_scalar_payload_to_complex(data))
        } else {
            Self::from_diag(indices, Self::any_scalar_payload_to_real(data))
        }
    }

    /// Create a copy tensor whose nonzero entries are `value` on the diagonal.
    ///
    /// For indices `[i, j, k]`, the returned tensor satisfies
    /// `T[i, j, k] = value` when `i = j = k`, and zero otherwise.
    ///
    /// # Errors
    /// Returns an error when the index dimensions are unequal (a dimension
    /// mismatch) or the construction fails.
    /// # Examples
    /// ```
    /// use tensor4all_core::{AnyScalar, IdxTensor};
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(2);
    /// let k = Index::new_dyn(2);
    /// let tensor = IdxTensor::copy_tensor(
    ///     vec![i, j, k],
    ///     AnyScalar::new_real(1.0),
    /// ).unwrap();
    ///
    /// assert_eq!(tensor.dims(), vec![2, 2, 2]);
    /// ```
    pub fn copy_tensor(
        indices: Vec<DynIndex>,
        value: AnyScalar,
    ) -> std::result::Result<Self, IdxTensorError> {
        if indices.is_empty() {
            return Self::from_dense_any(vec![], vec![value]);
        }
        let dim = indices[0].dim();
        let data = vec![value; dim];
        Self::from_diag_any(indices, data)
    }

    /// Replace multiple tensor indices with one fused index using an exact local reshape.
    ///
    /// The full indices in `old_indices` identify the axes to fuse and also
    /// define the coordinate order used inside `new_index`. The new fused index
    /// is inserted at the earliest axis position among the fused axes; all
    /// other axes keep their original relative order. Use
    /// [`LinearizationOrder::ColumnMajor`] to match tensor4all's dense vector
    /// layout, or [`LinearizationOrder::RowMajor`] when interoperating with
    /// row-major fused coordinates.
    ///
    /// # Arguments
    /// * `old_indices` - Non-empty list of existing tensor indices to replace.
    ///
    ///   Each index is matched by full identity, must appear exactly once in
    ///   the tensor, must have the same dimension as the matched tensor axis,
    ///   and must not be duplicated in this list.
    /// * `new_index` - Replacement index whose dimension must equal the product
    ///
    ///   of the dimensions in `old_indices`.
    /// * `order` - Linearization convention used to encode the old coordinates
    ///
    ///   into the single coordinate of `new_index`.
    ///
    /// # Returns
    /// A tensor with the same element type and values, but with `old_indices`
    /// replaced by `new_index`.
    ///
    /// # Errors
    /// Returns an error if `old_indices` is empty, contains duplicate IDs,
    /// references an index not present in the tensor, if the fused dimension
    /// does not match the product of the old dimensions, if the replacement
    /// would duplicate a kept index, or if the dense reshape cannot be
    /// represented without overflow.
    ///
    /// # Examples
    /// ```
    /// use tensor4all_core::{DynIndex, LinearizationOrder, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let fused = DynIndex::new_link(4).unwrap();
    /// let tensor = IdxTensor::from_dense(
    ///     vec![i.clone(), j.clone()],
    ///     vec![1.0, 2.0, 3.0, 4.0],
    /// ).unwrap();
    ///
    /// let fused_tensor = tensor
    ///     .fuse_indices(&[i.clone(), j.clone()], fused.clone(), LinearizationOrder::ColumnMajor)
    ///     .unwrap();
    /// assert_eq!(fused_tensor.dims(), vec![4]);
    ///
    /// let roundtrip = fused_tensor
    ///     .unfuse_index(&fused, &[i, j], LinearizationOrder::ColumnMajor)
    ///     .unwrap();
    /// assert!(roundtrip.isapprox(&tensor, 1e-12, 0.0).unwrap());
    /// ```
    pub fn fuse_indices(
        &self,
        old_indices: &[DynIndex],
        new_index: DynIndex,
        order: LinearizationOrder,
    ) -> std::result::Result<Self, IdxTensorError> {
        if !(!old_indices.is_empty()) {
            return Err(anyhow::anyhow!("fuse_indices requires at least one index to fuse").into());
        };

        let old_dims = self.dims();
        let mut seen_indices = HashSet::new();
        let mut old_axes = Vec::with_capacity(old_indices.len());
        for old_index in old_indices {
            if !(seen_indices.insert(old_index)) {
                return Err(anyhow::anyhow!("duplicate index in old_indices").into());
            };
            let axis = self
                .indices
                .iter()
                .position(|idx| idx == old_index)
                .ok_or_else(|| anyhow::anyhow!("index {:?} not found in tensor", old_index))?;
            if !(old_index.dim() == old_dims[axis]) {
                return Err(anyhow::anyhow!(
                    "old index dimension does not match tensor axis dimension"
                )
                .into());
            };
            old_axes.push(axis);
        }

        let fused_dims: Vec<usize> = old_axes.iter().map(|&axis| old_dims[axis]).collect();
        let fused_product = checked_product(&fused_dims)?;
        if !(fused_product == new_index.dim()) {
            return Err(anyhow::anyhow!(
                "product of old index dimensions must match the replacement index dimension"
            )
            .into());
        };

        let insertion_axis =
            old_axes.iter().copied().min().ok_or_else(|| {
                anyhow::anyhow!("fuse_indices requires at least one index to fuse")
            })?;
        let old_axis_set: HashSet<usize> = old_axes.iter().copied().collect();

        let mut result_indices =
            Vec::with_capacity(self.indices.len() - old_indices.len() + 1usize);
        for (axis, index) in self.indices.iter().enumerate() {
            if axis == insertion_axis {
                result_indices.push(new_index.clone());
            }
            if !old_axis_set.contains(&axis) {
                result_indices.push(index.clone());
            }
        }
        let mut result_seen = HashSet::new();
        for index in &result_indices {
            if !(result_seen.insert(index)) {
                return Err(
                    anyhow::anyhow!("fuse_indices result would contain duplicate index").into(),
                );
            };
        }
        Self::validate_indices(&result_indices)?;

        let mut new_dims = Vec::with_capacity(old_dims.len() - old_indices.len() + 1usize);
        for (axis, dim) in old_dims.iter().copied().enumerate() {
            if axis == insertion_axis {
                new_dims.push(new_index.dim());
            }
            if !old_axis_set.contains(&axis) {
                new_dims.push(dim);
            }
        }

        self.ensure_shape_packing_preserves_ad("fuse_indices")?;

        let mut grouped_axes = old_axes.clone();
        if matches!(order, LinearizationOrder::RowMajor) {
            grouped_axes.reverse();
        }
        let mut perm = Vec::with_capacity(self.indices.len());
        perm.extend((0..insertion_axis).filter(|axis| !old_axis_set.contains(axis)));
        perm.extend(grouped_axes);
        perm.extend(
            ((insertion_axis + 1)..self.indices.len()).filter(|axis| !old_axis_set.contains(axis)),
        );
        debug_assert_eq!(perm.len(), self.indices.len());

        let packed = self.permute(&perm)?;
        let reshaped = packed.try_materialized_inner()?.reshape(&new_dims)?;
        Self::from_inner(result_indices, reshaped).map_err(IdxTensorError::from)
    }

    /// Replace one fused index with multiple indices using an exact reshape.
    ///
    /// The caller must specify how the old fused index should be decoded into
    /// the new indices via `order`.
    ///
    /// # Errors
    /// Returns an error when the fused dimension does not equal the product of
    /// the new index dimensions (a shape mismatch).
    /// # Examples
    /// ```
    /// use tensor4all_core::{DynIndex, LinearizationOrder, IdxTensor};
    ///
    /// let fused = DynIndex::new_dyn(4);
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(vec![fused.clone()], vec![1.0, 2.0, 3.0, 4.0]).unwrap();
    ///
    /// let unfused = tensor
    ///     .unfuse_index(&fused, &[i.clone(), j.clone()], LinearizationOrder::ColumnMajor)
    ///     .unwrap();
    ///
    /// let expected = IdxTensor::from_dense(vec![i, j], vec![1.0, 2.0, 3.0, 4.0]).unwrap();
    /// assert!(unfused.isapprox(&expected, 1e-12, 0.0).unwrap());
    /// ```
    pub fn unfuse_index(
        &self,
        old_index: &DynIndex,
        new_indices: &[DynIndex],
        order: LinearizationOrder,
    ) -> std::result::Result<Self, IdxTensorError> {
        if !(!new_indices.is_empty()) {
            return Err(
                anyhow::anyhow!("unfuse_index requires at least one replacement index").into(),
            );
        };

        let axis = self
            .indices
            .iter()
            .position(|idx| idx == old_index)
            .ok_or_else(|| anyhow::anyhow!("index {:?} not found in tensor", old_index))?;

        let replacement_dims: Vec<usize> = new_indices.iter().map(DynIndex::dim).collect();
        let replacement_product = checked_product(&replacement_dims)?;
        if !(replacement_product == old_index.dim()) {
            return Err(anyhow::anyhow!(
                "product of new index dimensions must match the replaced index dimension"
            )
            .into());
        };

        let mut result_indices =
            Vec::with_capacity(self.indices.len() - 1usize + new_indices.len());
        result_indices.extend_from_slice(&self.indices[..axis]);
        result_indices.extend(new_indices.iter().cloned());
        result_indices.extend_from_slice(&self.indices[axis + 1..]);
        Self::validate_indices(&result_indices)?;

        let old_dims = self.dims();
        let mut new_dims = Vec::with_capacity(old_dims.len() - 1usize + replacement_dims.len());
        new_dims.extend_from_slice(&old_dims[..axis]);
        new_dims.extend_from_slice(&replacement_dims);
        new_dims.extend_from_slice(&old_dims[axis + 1..]);

        self.ensure_shape_packing_preserves_ad("unfuse_index")?;

        let mut grouped_indices = new_indices.to_vec();
        let mut grouped_dims = replacement_dims.clone();
        if matches!(order, LinearizationOrder::RowMajor) {
            grouped_indices.reverse();
            grouped_dims.reverse();
        }
        let mut packed_indices =
            Vec::with_capacity(self.indices.len() - 1usize + grouped_indices.len());
        packed_indices.extend_from_slice(&self.indices[..axis]);
        packed_indices.extend(grouped_indices);
        packed_indices.extend_from_slice(&self.indices[axis + 1..]);

        let mut packed_dims = Vec::with_capacity(old_dims.len() - 1usize + grouped_dims.len());
        packed_dims.extend_from_slice(&old_dims[..axis]);
        packed_dims.extend_from_slice(&grouped_dims);
        packed_dims.extend_from_slice(&old_dims[axis + 1..]);

        let reshaped = self.try_materialized_inner()?.reshape(&packed_dims)?;
        let packed = Self::from_inner(packed_indices, reshaped)?;
        if matches!(order, LinearizationOrder::ColumnMajor) {
            Ok(packed)
        } else {
            packed.permute_indices(&result_indices)
        }
    }

    /// Create a scalar (0-dimensional) tensor from a supported element value.
    ///
    /// # Errors
    /// Returns an error when the element type is not supported (an
    /// unsupported-dtype failure).
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    ///
    /// let scalar = IdxTensor::scalar(42.0).unwrap();
    /// assert_eq!(scalar.dims(), Vec::<usize>::new());
    /// assert_eq!(scalar.only().unwrap().real(), 42.0);
    /// ```
    pub fn scalar<T: TensorElement>(value: T) -> std::result::Result<Self, IdxTensorError> {
        Self::from_dense(vec![], vec![value])
    }

    /// Create a tensor filled with zeros of a supported element type.
    ///
    /// # Errors
    /// Returns an error when the dimension product overflows (an overflow failure)
    /// or the element type is unsupported.
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let j = Index::new_dyn(3);
    /// let tensor = IdxTensor::zeros::<f64>(vec![i, j]).unwrap();
    /// assert_eq!(tensor.dims(), vec![2, 3]);
    /// ```
    pub fn zeros<T: TensorElement + Zero + Clone>(
        indices: Vec<DynIndex>,
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims: Vec<usize> = indices.iter().map(|idx| idx.dim()).collect();
        let size = checked_product(&dims)?;
        Self::from_dense(indices, vec![T::zero(); size])
    }
}

// ============================================================================
// High-level API for data extraction (avoids direct .storage() access)
// ============================================================================

impl IdxTensor {
    /// Extract tensor data as a column-major `Vec<T>`.
    ///
    /// # Type Parameters
    /// * `T` - The scalar element type (`f32`, `f64`, `Complex32`, or
    ///
    ///   `Complex64`).
    ///
    /// # Returns
    /// A vector of the tensor data in column-major order.
    ///
    /// # Errors
    /// Returns an error when the tensor dtype does not match the requested element
    /// type (a scalar-kind mismatch) or materialization fails.
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(vec![i], vec![1.0, 2.0]).unwrap();
    /// let data = tensor.to_vec::<f64>().unwrap();
    /// assert_eq!(data, &[1.0, 2.0]);
    /// ```
    pub fn to_vec<T: TensorElement>(&self) -> std::result::Result<Vec<T>, IdxTensorError> {
        self.as_inner()?
            .duplicate_value()?
            .as_slice::<T>()
            .map(|values| values.to_vec())
            .map_err(|source| IdxTensorError::materialization(anyhow::Error::new(source)))
    }

    /// Reads dense column-major tensor values without copying them into a new vector.
    ///
    /// The callback receives the same logical dense values and ordering as
    /// [`Self::to_vec`], with the first tensor index varying fastest. Use this
    /// for read-only kernels that can finish while the callback is active;
    /// use [`Self::to_vec`] when the values must outlive the call.
    ///
    /// # Arguments
    ///
    /// * `read` - A callback that consumes the temporary dense value slice and
    ///   returns the caller's result.
    ///
    /// # Returns
    ///
    /// Returns the callback's result without allocating a result vector.
    ///
    /// # Errors
    ///
    /// Returns an error when the tensor dtype does not match `T` or dense
    /// materialization fails.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(
    ///     vec![i, j],
    ///     vec![1.0_f64, 2.0, 3.0, 4.0],
    /// )?;
    /// let weighted_sum = tensor.with_dense_slice::<f64, _>(|values| {
    ///     values.iter().enumerate().map(|(i, value)| (i + 1) as f64 * value).sum::<f64>()
    /// })?;
    /// assert_eq!(weighted_sum, 30.0);
    /// # Ok::<(), anyhow::Error>(())
    /// ```
    pub fn with_dense_slice<T: TensorElement, R>(
        &self,
        read: impl FnOnce(&[T]) -> R,
    ) -> std::result::Result<R, IdxTensorError> {
        let inner = self.as_inner()?;
        let value = inner.value()?;
        if let Ok(values) = value.as_slice::<T>() {
            return Ok(read(values));
        }

        // Backend-resident and non-contiguous values cannot be borrowed as a
        // host slice. Preserve the previous materializing behavior for those
        // cases while avoiding a full payload copy for the ordinary
        // host-contiguous path above.
        drop(value);
        let values = inner.duplicate_value()?;
        let values = values
            .as_slice::<T>()
            .map_err(|source| IdxTensorError::materialization(anyhow::Error::new(source)))?;
        Ok(read(values))
    }

    /// Consume the tensor and return its indices with dense column-major values.
    ///
    /// Use this when a caller needs to move index metadata and dense payload
    /// values across an API boundary. The returned values are ordered with the
    /// first tensor index varying fastest. Compact diagonal or structured
    /// storage is materialized into dense logical values.
    ///
    /// # Type Parameters
    /// * `T` - The scalar element type to extract: `f32`, `f64`, `Complex32`,
    ///
    ///   or `Complex64`.
    ///
    /// # Returns
    /// The tensor's original indices and dense column-major flat data.
    ///
    /// # Errors
    /// Returns an error when the tensor dtype does not match the requested element
    /// type (a scalar-kind mismatch) or materialization fails.
    /// # Examples
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(
    ///     vec![i.clone(), j.clone()],
    ///     vec![1.0_f64, 2.0, 3.0, 4.0],
    /// ).unwrap();
    ///
    /// let (indices, data) = tensor.into_dense_col_major_parts::<f64>().unwrap();
    ///
    /// assert_eq!(indices, vec![i, j]);
    /// assert_eq!(data, vec![1.0, 2.0, 3.0, 4.0]);
    /// ```
    pub fn into_dense_col_major_parts<T: TensorElement>(
        self,
    ) -> std::result::Result<(Vec<DynIndex>, Vec<T>), IdxTensorError> {
        if !(!self.tracks_grad()) {
            return Err(anyhow::anyhow!("IdxTensor::into_dense_col_major_parts cannot consume tensors with tracked autodiff state").into());
        };
        let data = self.to_vec::<T>()?;
        Ok((self.indices, data))
    }

    /// Check if the tensor has `f64` storage.
    ///
    /// # Example
    /// ```
    /// use tensor4all_core::IdxTensor;
    /// use tensor4all_core::index::{DefaultIndex as Index, DynId};
    ///
    /// let i = Index::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(vec![i], vec![1.0, 2.0]).unwrap();
    /// assert!(tensor.is_f64());
    /// assert!(!tensor.is_complex());
    /// ```
    pub fn is_f64(&self) -> bool {
        self.storage.dtype() == Some(DType::F64)
    }

    /// Check if the tensor has `f32` storage.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let tensor = IdxTensor::from_dense(
    ///     vec![DynIndex::new_dyn(2)],
    ///     vec![1.0_f32, 2.0],
    /// )
    /// .unwrap();
    /// assert!(tensor.is_f32());
    /// ```
    pub fn is_f32(&self) -> bool {
        self.storage.dtype() == Some(DType::F32)
    }

    /// Check if the tensor has complex-32 storage.
    ///
    /// # Examples
    ///
    /// ```
    /// use num_complex::Complex32;
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let tensor = IdxTensor::from_dense(
    ///     vec![DynIndex::new_dyn(2)],
    ///     vec![Complex32::new(1.0, 0.0), Complex32::new(0.0, 1.0)],
    /// )
    /// .unwrap();
    /// assert!(tensor.is_c32());
    /// ```
    pub fn is_c32(&self) -> bool {
        self.storage.dtype() == Some(DType::C32)
    }

    /// Check if the tensor has complex-64 storage.
    ///
    /// # Example
    /// ```
    /// use num_complex::Complex64;
    /// use tensor4all_core::{DynIndex, IdxTensor};
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense(
    ///     vec![i],
    ///     vec![Complex64::new(1.0, 0.0), Complex64::new(0.0, 1.0)],
    /// )
    /// .unwrap();
    /// assert!(tensor.is_c64());
    /// ```
    pub fn is_c64(&self) -> bool {
        self.storage.is_c64()
    }

    /// Check whether the tensor carries diagonal logical axis metadata.
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use tensor4all_tensorbackend::Storage;
    ///
    /// // Tensors from `from_dense` use dense storage
    /// let i = DynIndex::new_dyn(2);
    /// let j = DynIndex::new_dyn(2);
    /// let dense = IdxTensor::from_dense(vec![i, j], vec![1.0, 0.0, 0.0, 1.0]).unwrap();
    /// assert!(!dense.is_diag());
    ///
    /// // Diagonal metadata is preserved when constructing from diagonal storage.
    /// let k = DynIndex::new_dyn(2);
    /// let l = DynIndex::new_dyn(2);
    /// let diag = IdxTensor::from_storage(
    ///     vec![k, l],
    ///     Storage::from_diag_col_major(vec![1.0, 2.0], 2)
    ///         .map(std::sync::Arc::new)
    ///         .unwrap(),
    /// )
    /// .unwrap();
    /// assert!(diag.is_diag());
    /// ```
    pub fn is_diag(&self) -> bool {
        self.storage.is_diag()
    }

    /// Check if the tensor has complex storage (C64).
    ///
    /// # Examples
    ///
    /// ```
    /// use tensor4all_core::{DynIndex, IdxTensor};
    /// use num_complex::Complex64;
    ///
    /// let i = DynIndex::new_dyn(2);
    /// let real_t = IdxTensor::from_dense(vec![i.clone()], vec![1.0, 2.0]).unwrap();
    /// assert!(!real_t.is_complex());
    ///
    /// let complex_t = IdxTensor::from_dense(
    ///     vec![i],
    ///     vec![Complex64::new(1.0, 0.0), Complex64::new(0.0, 1.0)],
    /// ).unwrap();
    /// assert!(complex_t.is_complex());
    /// ```
    pub fn is_complex(&self) -> bool {
        self.storage.is_complex()
    }
    /// Create a dense tensor in a caller-owned execution context.
    ///
    /// Host data is explicitly uploaded when `context` is CUDA-resident; no
    /// global default context is consulted.
    ///
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor};
    /// use tensor4all_tensorbackend::CpuExecutionContext;
    /// use tenferro_cpu::CpuBackend;
    ///
    /// let context = ExecutionContext::Cpu(Arc::new(
    ///     CpuExecutionContext::from_backend(CpuBackend::new()),
    /// ));
    /// let index = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense_in(&context, vec![index], vec![1.0_f64, 2.0])?;
    /// assert_eq!(tensor.to_vec::<f64>()?, vec![1.0, 2.0]);
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the indices, payload, context, or
    /// explicit transfer is invalid.
    pub fn from_dense_in<T: TensorElement>(
        context: &ExecutionContext,
        indices: Vec<DynIndex>,
        data: Vec<T>,
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims = Self::expected_dims_from_indices(&indices);
        Self::validate_indices(&indices)?;
        Self::validate_dense_payload_len(data.len(), &dims)?;
        let native = dense_native_tensor_from_col_major(&data, &dims)
            .map_err(|error| IdxTensorError::operation("context-scoped construction", error))?;
        let inner = match context {
            ExecutionContext::Cpu(context) => {
                let runtime = context.eager_runtime().map_err(|error| {
                    IdxTensorError::operation("CPU context construction", anyhow::Error::new(error))
                })?;
                EagerTensor::from_tensor_in(native, runtime).map_err(|error| {
                    IdxTensorError::operation(
                        "CPU context eager wrapping",
                        anyhow::Error::new(error),
                    )
                })?
            }
            #[cfg(feature = "tenferro-cuda")]
            ExecutionContext::Cuda(context) => {
                let uploaded = context.upload_cuda(&native).map_err(|error| {
                    IdxTensorError::operation(
                        "CUDA context construction",
                        anyhow::Error::new(error),
                    )
                })?;
                EagerTensor::from_tensor_in(
                    uploaded,
                    context.eager_runtime().map_err(|error| {
                        IdxTensorError::operation(
                            "CUDA context eager runtime",
                            anyhow::Error::new(error),
                        )
                    })?,
                )
                .map_err(|error| {
                    IdxTensorError::operation(
                        "CUDA context eager wrapping",
                        anyhow::Error::new(error),
                    )
                })?
            }
        };
        Self::from_inner(indices, inner).map_err(IdxTensorError::from)
    }

    /// Create an all-ones tensor in a caller-owned execution context.
    ///
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor};
    /// use tensor4all_tensorbackend::CpuExecutionContext;
    /// use tenferro_cpu::CpuBackend;
    ///
    /// let context = ExecutionContext::Cpu(Arc::new(
    ///     CpuExecutionContext::from_backend(CpuBackend::new()),
    /// ));
    /// let index = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::ones_in(&context, &[index])?;
    /// assert_eq!(tensor.to_vec::<f64>()?, vec![1.0, 1.0]);
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the index dimensions or context
    /// construction is invalid.
    pub fn ones_in(
        context: &ExecutionContext,
        indices: &[DynIndex],
    ) -> std::result::Result<Self, IdxTensorError> {
        let dims = Self::expected_dims_from_indices(indices);
        let total_size = checked_total_size(&dims)?;
        Self::from_dense_in(context, indices.to_vec(), vec![1.0_f64; total_size])
    }

    /// Validate that this tensor belongs to the supplied execution context.
    ///
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor};
    /// use tensor4all_tensorbackend::CpuExecutionContext;
    /// use tenferro_cpu::CpuBackend;
    ///
    /// let context = ExecutionContext::Cpu(Arc::new(
    ///     CpuExecutionContext::from_backend(CpuBackend::new()),
    /// ));
    /// let index = DynIndex::new_dyn(2);
    /// let tensor = IdxTensor::from_dense_in(&context, vec![index], vec![1.0_f64, 2.0])?;
    /// tensor.validate_context(&context)?;
    /// assert_eq!(tensor.dims(), vec![2]);
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the tensor does not belong to `context`.
    pub fn validate_context(
        &self,
        context: &ExecutionContext,
    ) -> std::result::Result<(), IdxTensorError> {
        let Some(inner) = self.storage.eager() else {
            if matches!(context, ExecutionContext::Cpu(_)) {
                return Ok(());
            }
            return Err(IdxTensorError::operation(
                "context validation",
                anyhow::anyhow!("tensor has no CUDA eager runtime"),
            ));
        };
        match context {
            ExecutionContext::Cpu(context) => {
                let expected = context.eager_runtime().map_err(|error| {
                    IdxTensorError::operation("CPU context validation", anyhow::Error::new(error))
                })?;
                if inner.ctx_id() != expected.id() {
                    return Err(IdxTensorError::operation(
                        "context validation",
                        anyhow::anyhow!("tensor belongs to a different CPU eager context"),
                    ));
                }
            }
            #[cfg(feature = "tenferro-cuda")]
            ExecutionContext::Cuda(context) => {
                self.validate_cuda_residency(context).map_err(|error| {
                    IdxTensorError::operation("CUDA context validation", anyhow::Error::new(error))
                })?;
            }
        }
        Ok(())
    }

    /// Read back a small algorithm-decision tensor through its owning context.
    ///
    /// This is the explicit synchronization/readback boundary for rank
    /// selection and adaptive stopping payloads (singular values, norm
    /// vectors, scalar estimates). With a CUDA context the resident value is
    /// downloaded through that exact context; no arithmetic operand is ever
    /// reconstructed on the host from this payload.
    ///
    /// # Arguments
    ///
    /// * `context` - Caller-owned execution context the tensor must belong to.
    ///
    /// # Returns
    ///
    /// The decision values in column-major order (`f64`, or the real part for
    /// `Complex64`, matching the existing rank-selection semantics).
    ///
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor};
    /// use tensor4all_tensorbackend::CpuExecutionContext;
    /// use tenferro_cpu::CpuBackend;
    ///
    /// let context = ExecutionContext::Cpu(Arc::new(
    ///     CpuExecutionContext::from_backend(CpuBackend::new()),
    /// ));
    /// let tensor = IdxTensor::from_dense_in(
    ///     &context,
    ///     vec![DynIndex::new_dyn(3)],
    ///     vec![0.5_f64, 2.0, 1.0],
    /// )?;
    /// assert_eq!(tensor.read_decision_data(&context)?, vec![0.5, 2.0, 1.0]);
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the tensor does not belong to `context`,
    /// when the explicit transfer fails, or when the storage is neither `f64`
    /// nor `Complex64`.
    pub fn read_decision_data(
        &self,
        context: &ExecutionContext,
    ) -> std::result::Result<Vec<f64>, IdxTensorError> {
        self.validate_context(context)?;
        #[cfg(feature = "tenferro-cuda")]
        if let ExecutionContext::Cuda(cuda) = context {
            let host = self.download(cuda).map_err(|error| {
                IdxTensorError::operation("decision readback download", anyhow::Error::new(error))
            })?;
            return Self::decision_values(&host);
        }
        Self::decision_values(self)
    }

    fn decision_values(host: &Self) -> std::result::Result<Vec<f64>, IdxTensorError> {
        if host.is_f64() {
            host.to_vec::<f64>().map_err(|error| {
                IdxTensorError::operation("decision readback", anyhow::Error::new(error))
            })
        } else if host.is_c64() {
            host.to_vec::<Complex64>()
                .map(|values| values.into_iter().map(|value| value.re).collect())
                .map_err(|error| {
                    IdxTensorError::operation("decision readback", anyhow::Error::new(error))
                })
        } else {
            Err(IdxTensorError::operation(
                "decision readback",
                anyhow::anyhow!("decision payloads support f64 and Complex64 storage only"),
            ))
        }
    }

    /// Scale every element by a real factor in a caller-owned context.
    ///
    /// The scalar is built in the tensor's owning runtime (explicit upload
    /// for CUDA) and multiplied there; host-storage CPU tensors are adopted
    /// into `context` first via explicit construction. With a CPU context
    /// this matches [`IdxTensor::scale`] bit-for-bit.
    ///
    /// # Arguments
    ///
    /// * `factor` - Real scaling factor.
    /// * `context` - Caller-owned execution context the tensor belongs to.
    ///
    /// # Examples
    ///
    /// ```
    /// use std::sync::Arc;
    /// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor};
    /// use tensor4all_tensorbackend::CpuExecutionContext;
    /// use tenferro_cpu::CpuBackend;
    ///
    /// let context = ExecutionContext::Cpu(Arc::new(
    ///     CpuExecutionContext::from_backend(CpuBackend::new()),
    /// ));
    /// let tensor = IdxTensor::from_dense_in(
    ///     &context,
    ///     vec![DynIndex::new_dyn(2)],
    ///     vec![1.0_f64, 2.0],
    /// )?;
    /// let scaled = tensor.scale_in(3.0, &context)?;
    /// assert_eq!(scaled.to_vec::<f64>()?, vec![3.0, 6.0]);
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the tensor does not belong to
    /// `context`, when the storage is neither `f64` nor `Complex64`, or when
    /// the scalar construction, upload, or multiplication fails.
    pub fn scale_in(
        &self,
        factor: f64,
        context: &ExecutionContext,
    ) -> std::result::Result<Self, IdxTensorError> {
        self.validate_context(context)?;
        let adopted_storage;
        let operand: &EagerTensor = match self.eager_runtime() {
            Some(_) => self.as_inner()?,
            // Host storage (CPU only; CUDA validated above): adopt into the
            // context through explicit construction.
            None => {
                adopted_storage = Self::adopt_host_into_context(self, context)?;
                adopted_storage.as_inner()?
            }
        };
        let scalar = Self::context_scalar_in(operand, factor, context)?;
        let scaled = operand.mul(&scalar).map_err(|error| {
            IdxTensorError::operation(
                "context-scoped scaling",
                anyhow::Error::new(error).context("eager multiplication failed"),
            )
        })?;
        Self::from_inner(self.indices.clone(), scaled).map_err(IdxTensorError::from)
    }

    /// Adopt host storage into a context through explicit construction.
    fn adopt_host_into_context(
        tensor: &Self,
        context: &ExecutionContext,
    ) -> std::result::Result<Self, IdxTensorError> {
        if tensor.is_f64() {
            let values = tensor.to_vec::<f64>().map_err(|error| {
                IdxTensorError::operation("context adoption", anyhow::Error::new(error))
            })?;
            Self::from_dense_in(context, tensor.indices.clone(), values)
        } else if tensor.is_c64() {
            let values = tensor.to_vec::<Complex64>().map_err(|error| {
                IdxTensorError::operation("context adoption", anyhow::Error::new(error))
            })?;
            Self::from_dense_in(context, tensor.indices.clone(), values)
        } else {
            Err(IdxTensorError::operation(
                "context adoption",
                anyhow::anyhow!("adoption supports f64 and Complex64 storage only"),
            ))
        }
    }

    /// Build a scalar in the operand's owning runtime (explicit upload for CUDA).
    fn context_scalar_in(
        operand: &EagerTensor,
        factor: f64,
        context: &ExecutionContext,
    ) -> std::result::Result<EagerTensor, IdxTensorError> {
        let native: NativeTensor = match operand.dtype() {
            DType::F64 => NativeTensor::from_vec_col_major(vec![], vec![factor]),
            DType::C64 => {
                NativeTensor::from_vec_col_major(vec![], vec![Complex64::new(factor, 0.0)])
            }
            other => {
                return Err(IdxTensorError::operation(
                    "context-scoped scaling",
                    anyhow::anyhow!("unsupported scalar dtype {other:?}"),
                ));
            }
        }
        .map_err(|error| {
            IdxTensorError::operation("context-scoped scaling", anyhow::anyhow!("{error}"))
        })?;
        let runtime = Self::context_eager_runtime(context)?;
        let resident = match context {
            ExecutionContext::Cpu(_) => native,
            #[cfg(feature = "tenferro-cuda")]
            ExecutionContext::Cuda(cuda) => cuda.upload_cuda(&native).map_err(|error| {
                IdxTensorError::operation(
                    "context-scoped scaling",
                    anyhow::Error::new(error).context("scalar upload failed"),
                )
            })?,
        };
        EagerTensor::from_tensor_in(resident, runtime).map_err(|error| {
            IdxTensorError::operation(
                "context-scoped scaling",
                anyhow::Error::new(error).context("scalar wrapping failed"),
            )
        })
    }

    /// Resolve the eager runtime owned by a context.
    fn context_eager_runtime(
        context: &ExecutionContext,
    ) -> std::result::Result<Arc<EagerRuntime>, IdxTensorError> {
        match context {
            ExecutionContext::Cpu(context) => context.eager_runtime().map_err(|error| {
                IdxTensorError::operation("context-scoped scaling", anyhow::Error::new(error))
            }),
            #[cfg(feature = "tenferro-cuda")]
            ExecutionContext::Cuda(context) => context.eager_runtime().map_err(|error| {
                IdxTensorError::operation("context-scoped scaling", anyhow::Error::new(error))
            }),
        }
    }

    /// Frobenius norm in a caller-owned execution context.
    ///
    /// Reductions run in the tensor's owning runtime; only the scalar result
    /// crosses the explicit readback boundary on CUDA. With a CPU context
    /// this delegates to [`IdxTensor::norm`] unchanged.
    ///
    /// # Errors
    ///
    /// Returns [`IdxTensorError`] when the tensor does not belong to
    /// `context`, when the storage is neither `f64` nor `Complex64`, or when
    /// the reductions or explicit readback fail.
    pub fn norm_in(&self, context: &ExecutionContext) -> std::result::Result<f64, IdxTensorError> {
        self.validate_context(context)?;
        #[cfg(feature = "tenferro-cuda")]
        if matches!(context, ExecutionContext::Cuda(_)) {
            let inner = self.cuda_eager_inner().ok_or_else(|| {
                IdxTensorError::operation(
                    "context-scoped norm",
                    anyhow::anyhow!("tensor has no resident eager value"),
                )
            })?;
            let rank = inner.shape().len();
            let axes: Vec<usize> = (0..rank).collect();
            let squared = inner
                .abs()
                .and_then(|magnitudes| magnitudes.reduce_sum_squares(&axes))
                .map_err(|error| {
                    IdxTensorError::operation("context-scoped norm", anyhow::anyhow!("{error}"))
                })?;
            let scalar = squared.reshape(&[1]).map_err(|error| {
                IdxTensorError::operation("context-scoped norm", anyhow::anyhow!("{error}"))
            })?;
            let values = Self::read_resident_vector(&scalar, 1, context).map_err(|error| {
                IdxTensorError::operation("context-scoped norm", anyhow::Error::new(error))
            })?;
            return Ok(values[0].sqrt());
        }
        #[cfg(not(feature = "tenferro-cuda"))]
        let _ = context;
        self.norm()
    }
}

fn checked_product(dims: &[usize]) -> Result<usize> {
    dims.iter().try_fold(1usize, |acc, &dim| {
        acc.checked_mul(dim)
            .ok_or_else(|| anyhow::anyhow!("dimension product overflow"))
    })
}

fn increment_col_major_coordinate(coords: &mut [usize], dims: &[usize]) {
    let mut carry = true;
    for (coordinate, &dim) in coords.iter_mut().zip(dims.iter()) {
        if !carry {
            break;
        }
        *coordinate += 1;
        if *coordinate == dim {
            *coordinate = 0;
        } else {
            carry = false;
        }
    }
}

fn map_payload_support_coordinate(
    source_axis_classes: &[usize],
    target_axis_classes: &[usize],
    source_to_target_axes: &[usize],
    source_coords: &[usize],
    target_coords: &mut [usize],
    target_seen: &mut [bool],
) -> bool {
    if source_axis_classes.len() != source_to_target_axes.len()
        || target_coords.len() != target_seen.len()
    {
        return false;
    }
    target_seen.fill(false);
    for (source_axis, &target_axis) in source_to_target_axes.iter().enumerate() {
        let Some(&source_class) = source_axis_classes.get(source_axis) else {
            return false;
        };
        let Some(&target_class) = target_axis_classes.get(target_axis) else {
            return false;
        };
        let Some(&source_value) = source_coords.get(source_class) else {
            return false;
        };
        let Some(target_value) = target_coords.get_mut(target_class) else {
            return false;
        };
        if target_seen[target_class] {
            if *target_value != source_value {
                return false;
            }
        } else {
            *target_value = source_value;
            target_seen[target_class] = true;
        }
    }
    target_seen.iter().all(|&seen| seen)
}

fn decode_col_major_linear(linear: usize, dims: &[usize]) -> Result<Vec<usize>> {
    let total = checked_product(dims)?;
    if !(linear < total) {
        return Err(anyhow::anyhow!(
            "linear offset {} out of bounds for dims {:?}",
            linear,
            dims
        ));
    };
    let mut remaining = linear;
    let mut out = Vec::with_capacity(dims.len());
    for &dim in dims {
        out.push(remaining % dim);
        remaining /= dim;
    }
    Ok(out)
}

fn encode_col_major_linear(indices: &[usize], dims: &[usize]) -> Result<usize> {
    if !(indices.len() == dims.len()) {
        return Err(anyhow::anyhow!(
            "index rank {} does not match dims {:?}",
            indices.len(),
            dims
        ));
    };
    let mut linear = 0usize;
    let mut stride = 1usize;
    for (&index, &dim) in indices.iter().zip(dims.iter()) {
        if !(index < dim) {
            return Err(anyhow::anyhow!(
                "index {} out of bounds for dimension {}",
                index,
                dim
            ));
        };
        let term = index
            .checked_mul(stride)
            .ok_or_else(|| anyhow::anyhow!("linear offset overflow"))?;
        linear = linear
            .checked_add(term)
            .ok_or_else(|| anyhow::anyhow!("linear offset overflow"))?;
        stride = stride
            .checked_mul(dim)
            .ok_or_else(|| anyhow::anyhow!("stride overflow"))?;
    }
    Ok(linear)
}

#[cfg(test)]
mod tests {
    use super::*;
    use num_complex::{Complex32, Complex64};
    use std::cell::Cell;
    use tensor4all_tensorbackend::StorageError;

    #[test]
    fn structured_contraction_does_not_install_logical_dense_cache() {
        let n = 8;
        let left = DynIndex::new_dyn(n);
        let site = DynIndex::new_dyn(3);
        let right = DynIndex::new_dyn(n);
        let far = DynIndex::new_dyn(n);
        let end = DynIndex::new_dyn(n);
        let a =
            IdxTensor::from_copy_selector(left, site.clone(), right.clone(), 1, 1.0_f64).unwrap();
        let b =
            IdxTensor::from_copy_selector(right, site.clone(), far.clone(), 1, 2.0_f64).unwrap();
        let c = IdxTensor::from_copy_selector(far, site.clone(), end, 1, 3.0_f64).unwrap();
        let result = crate::defaults::contract::contract_with_options(
            &[&a, &b, &c],
            crate::defaults::contract::ContractionOptions::new()
                .with_retain_indices(std::slice::from_ref(&site)),
        )
        .unwrap();

        assert_eq!(result.storage_kind(), StorageKind::Structured);
        assert!(result.eager_cache.get().is_none());
        assert_eq!(result.storage().unwrap().payload_len(), n * 3);
    }

    #[test]
    fn binary_contraction_axis_classes_preserve_uncontracted_order() {
        let axis_classes =
            IdxTensor::binary_contraction_axis_classes(&[0, 1, 0], &[1], &[0, 1], &[0]).unwrap();

        assert_eq!(axis_classes, vec![0, 0, 1]);
    }

    #[test]
    fn structured_metrics_use_authoritative_compact_payload_for_all_dtypes() {
        fn check(tensor: IdxTensor, expected_sum: f64, expected_norm_squared: f64) {
            assert!(matches!(tensor.storage, IdxTensorStorage::Compact(_)));
            assert!(tensor.eager_cache.get().is_none());
            assert!((tensor.sum().unwrap().real() - expected_sum).abs() < 1.0e-6);
            assert!((tensor.norm_squared().unwrap() - expected_norm_squared).abs() < 1.0e-6);
            assert!((tensor.maxabs().unwrap() - 2.0).abs() < 1.0e-6);
            assert!(tensor.isapprox(&tensor, 0.0, 0.0).unwrap());
            assert!(tensor.eager_cache.get().is_none());
        }

        let indices = || vec![DynIndex::new_dyn(2), DynIndex::new_dyn(2)];
        check(
            IdxTensor::from_diag(indices(), vec![1.0_f32, 2.0]).unwrap(),
            3.0,
            5.0,
        );
        check(
            IdxTensor::from_diag(indices(), vec![1.0_f64, 2.0]).unwrap(),
            3.0,
            5.0,
        );
        check(
            IdxTensor::from_diag(
                indices(),
                vec![Complex32::new(1.0, 0.0), Complex32::new(2.0, 0.0)],
            )
            .unwrap(),
            3.0,
            5.0,
        );
        check(
            IdxTensor::from_diag(
                indices(),
                vec![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)],
            )
            .unwrap(),
            3.0,
            5.0,
        );
    }

    #[test]
    fn payload_storage_error_retains_typed_source() {
        let storage = Storage::from_dense_col_major(vec![1.0_f64, 2.0], &[2]).unwrap();
        let error = storage.scalar_at(&[2]).unwrap_err();
        assert!(matches!(error, StorageError::InvalidStructuredStorage(_)));
    }

    #[test]
    fn materialization_error_retains_backend_source() {
        let native = NativeTensor::from_vec_col_major(vec![2, 3], vec![1.0_f64; 6]).unwrap();
        let inner = EagerTensor::from_tensor_in(native, default_eager_ctx().unwrap()).unwrap();
        let storage = IdxTensorStorage::Eager {
            inner: Arc::new(inner),
            axis_classes: vec![0, 0],
        };

        let error = storage.materialize(2).unwrap_err();
        assert!(matches!(error, TensorStorageError::Materialization { .. }));
        assert!(std::error::Error::source(&error).is_some());
    }

    fn conjugate_with_injected_failure(
        tensor: &IdxTensor,
        target: *const EagerTensor,
        message: &'static str,
    ) -> IdxTensor {
        let calls = Cell::new(0usize);
        let conjugated = tensor.conj_with(&|inner| {
            calls.set(calls.get() + 1);
            if std::ptr::eq(inner, target) {
                Err(Arc::new(std::io::Error::other(message)) as _)
            } else {
                conjugate_eager(inner)
            }
        });
        assert!(calls.get() > 0, "injected closure was not reached");
        conjugated
    }

    fn assert_unwrapped_conjugation_error(
        tensor: IdxTensor,
        target: *const EagerTensor,
        message: &'static str,
    ) -> IdxTensor {
        let conjugated = conjugate_with_injected_failure(&tensor, target, message);
        let error = conjugated.to_storage().unwrap_err();
        assert!(matches!(error, TensorStorageError::Conjugation { .. }));
        let source = std::error::Error::source(&error).unwrap();
        assert_eq!(source.to_string(), message);
        assert!(
            source.source().is_none(),
            "source was wrapped more than once"
        );
        conjugated
    }

    #[test]
    fn authoritative_storage_conjugation_failure_is_deferred_without_detaching() {
        let i = DynIndex::new_dyn(2);
        let native = NativeTensor::from_vec_col_major(vec![2], vec![1.0_f64, 2.0]).unwrap();
        let inner = EagerTensor::requires_grad_in(native, default_eager_ctx().unwrap()).unwrap();
        let tensor = IdxTensor::from_inner(vec![i], inner).unwrap();
        let source = match &tensor.storage {
            IdxTensorStorage::Eager { inner, .. } => Arc::clone(inner),
            IdxTensorStorage::Compact(payload) => Arc::clone(&payload.payload),
            IdxTensorStorage::Materialized(_) | IdxTensorStorage::Deferred { .. } => {
                panic!("tracked eager source expected")
            }
        };

        let conjugated = conjugate_with_injected_failure(
            &tensor,
            Arc::as_ptr(&source),
            "forced authoritative eager conjugation failure",
        );
        assert!(conjugated.tracks_grad());
        assert!(conjugated.is_f64());
        assert!(!conjugated.is_complex());
        assert!(!conjugated.is_diag());
        assert_eq!(conjugated.dims(), vec![2]);

        let error = conjugated.to_storage().unwrap_err();
        let source = std::error::Error::source(&error).unwrap();
        assert_eq!(
            source.to_string(),
            "forced authoritative eager conjugation failure"
        );
        assert!(conjugated.detach().is_err());
    }

    #[test]
    fn structured_payload_conjugation_failure_retains_graph_and_blocks_detached_primal() {
        let i = DynIndex::new_dyn(2);
        let j = DynIndex::new_dyn(2);
        let tensor = IdxTensor::from_diag(
            vec![i, j],
            vec![Complex64::new(1.0, 2.0), Complex64::new(3.0, -4.0)],
        )
        .unwrap()
        .enable_grad()
        .unwrap();

        let target = Arc::as_ptr(
            &tensor
                .storage
                .compact_payload()
                .expect("tracked compact payload")
                .payload,
        );
        let conjugated = assert_unwrapped_conjugation_error(
            tensor,
            target,
            "forced structured AD conjugation failure",
        );
        assert!(conjugated.tracks_grad());
        assert!(conjugated.detach().is_err());
        assert!(conjugated.clone().enable_grad().is_err());
        assert!(conjugated.sum().is_err());
        assert!(conjugated.grad().is_err());
        assert!(conjugated.clear_grad().is_err());
        assert!(conjugated.maxabs().is_err());
        assert!(conjugated.norm_squared().is_err());

        let twice_conjugated = conjugated.conj();
        let error = twice_conjugated.to_storage().unwrap_err();
        assert_eq!(
            std::error::Error::source(&error).unwrap().to_string(),
            "forced structured AD conjugation failure"
        );
    }

    #[test]
    fn eager_cache_conjugation_failure_is_deferred_with_original_diagnostic() {
        let i = DynIndex::new_dyn(2);
        let j = DynIndex::new_dyn(2);
        let tensor = IdxTensor::from_diag(vec![i, j], vec![1.0_f64, 2.0]).unwrap();
        tensor.as_inner().unwrap();
        let target = Arc::as_ptr(tensor.eager_cache.get().unwrap());

        let conjugated = assert_unwrapped_conjugation_error(
            tensor,
            target,
            "forced eager cache conjugation failure",
        );
        assert!(!conjugated.tracks_grad());
        assert!(conjugated.detach().is_err());
    }

    #[test]
    fn encode_col_major_linear_rejects_offset_overflow() {
        let error =
            encode_col_major_linear(&[usize::MAX - 1, usize::MAX - 1], &[usize::MAX, usize::MAX])
                .unwrap_err();
        assert!(error.to_string().contains("linear offset overflow"));
    }

    #[test]
    fn factorize_probe_batch_incremental_matches_the_column_based_path() {
        let row = DynIndex::new_dyn(6);
        let batch = DynIndex::new_dyn(4);
        let data: Vec<f64> = (0..24).map(|i| i as f64 * 0.37 - 1.5).collect();
        let batch_tensor =
            IdxTensor::from_dense(vec![row.clone(), batch.clone()], data.clone()).unwrap();

        let columns: Vec<IdxTensor> = (0..4)
            .map(|position| {
                batch_tensor
                    .select_indices(std::slice::from_ref(&batch), &[position])
                    .unwrap()
            })
            .collect();
        let column_refs: Vec<&IdxTensor> = columns.iter().collect();

        let from_batch = IdxTensor::factorize_probe_batch_incremental(
            None,
            &batch_tensor,
            &batch,
            std::slice::from_ref(&row),
        )
        .unwrap();
        let from_columns = IdxTensor::factorize_probe_columns_incremental(
            None,
            &column_refs,
            &column_refs,
            &[row],
        )
        .unwrap();

        assert_eq!(from_batch.rank, from_columns.rank);
        let batch_data = from_batch.left.to_vec::<f64>().unwrap();
        let columns_data = from_columns.left.to_vec::<f64>().unwrap();
        assert_eq!(batch_data.len(), columns_data.len());
        let max_diff = batch_data
            .iter()
            .zip(columns_data.iter())
            .map(|(a, b)| (a - b).abs())
            .fold(0.0_f64, f64::max);
        assert!(
            max_diff < 1e-12,
            "batch-native and column-based factorizations disagree: max diff = {}",
            max_diff
        );
    }
}