scirs2-core 0.4.2

//! Memory layout optimizations with C/F order support for `NumPy` compatibility.
//!
//! This module provides comprehensive memory layout management that matches `NumPy`'s
//! array memory ordering conventions. It enables seamless interoperability with
//! `NumPy` and SciPy while optimizing performance for different access patterns.
//!
//! ## Features
//!
//! - **C-order (Row-major)**: Contiguous elements along the last axis
//! - **F-order (Column-major)**: Contiguous elements along the first axis
//! - **Automatic layout detection**: Determine optimal layout for operations
//! - **Layout conversion**: Efficient conversion between C and F order
//! - **Stride calculation**: Compute memory strides for any layout
//! - **Cache-friendly access**: Optimize memory access patterns
//! - **SIMD alignment**: Ensure proper alignment for vectorized operations
//!
//! ## Example Usage
//!
//! ```rust
//! use scirs2_core::memory_efficient::{
//!     MemoryLayout, ArrayLayout, LayoutOrder
//! };
//!
//! // Create a C-order (row-major) layout
//! let c_layout = MemoryLayout::new_c_order(&[100, 200]);
//!
//! // Create an F-order (column-major) layout  
//! let f_layout = MemoryLayout::new_f_order(&[100, 200]);
//!
//! // Convert between layouts
//! let converted = c_layout.to_order(LayoutOrder::Fortran).expect("Operation failed");
//!
//! // Check if layout is contiguous
//! assert!(c_layout.is_c_contiguous());
//! assert!(f_layout.is_f_contiguous());
//! ```

use crate::error::{CoreError, CoreResult, ErrorContext, ErrorLocation};
use ::ndarray::{Array, ArrayBase, Data, Dimension, RawData, ShapeBuilder};
use std::mem;

/// Memory layout order following `NumPy` conventions
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
pub enum LayoutOrder {
    /// C-order (row-major): last axis is contiguous
    #[default]
    C,
    /// Fortran-order (column-major): first axis is contiguous  
    Fortran,
    /// Any order (no specific preference)
    Any,
    /// Keep existing order (for transformations)
    Keep,
}

impl LayoutOrder {
    /// Convert to string representation
    pub const fn as_str(&self) -> &'static str {
        match self {
            LayoutOrder::C => "C",
            LayoutOrder::Fortran => "F",
            LayoutOrder::Any => "A",
            LayoutOrder::Keep => "K",
        }
    }

    /// Parse from string (`NumPy`-compatible)
    pub fn parse(s: &str) -> CoreResult<Self> {
        match s.to_uppercase().as_str() {
            "C" => Ok(LayoutOrder::C),
            "F" | "FORTRAN" => Ok(LayoutOrder::Fortran),
            "A" | "ANY" => Ok(LayoutOrder::Any),
            "K" | "KEEP" => Ok(LayoutOrder::Keep),
            _ => Err(CoreError::ValidationError(
                ErrorContext::new(format!("Invalid layout order: {s}"))
                    .with_location(ErrorLocation::new(file!(), line!())),
            )),
        }
    }
}

/// Memory layout information for arrays
#[derive(Debug, Clone)]
pub struct MemoryLayout {
    /// Shape of the array
    pub shape: Vec<usize>,
    /// Strides for each dimension (in bytes)
    pub strides: Vec<isize>,
    /// Element size in bytes
    pub element_size: usize,
    /// Total size in bytes
    pub total_size: usize,
    /// Layout order
    pub order: LayoutOrder,
    /// Whether the layout is contiguous
    pub is_contiguous: bool,
    /// Memory alignment in bytes
    pub alignment: usize,
}

impl MemoryLayout {
    /// Create a new C-order (row-major) layout
    pub fn new_c_order(shape: &[usize]) -> Self {
        Self::new_with_order(shape, LayoutOrder::C, mem::size_of::<f64>(), 64)
    }

    /// Create a new F-order (column-major) layout
    pub fn new_f_order(shape: &[usize]) -> Self {
        Self::new_with_order(shape, LayoutOrder::Fortran, mem::size_of::<f64>(), 64)
    }

    /// Create a new layout with specified order and element size
    pub fn new_with_order(
        shape: &[usize],
        order: LayoutOrder,
        element_size: usize,
        alignment: usize,
    ) -> Self {
        let strides = Self::calculate_strides(shape, order, element_size);
        let total_size = shape.iter().product::<usize>() * element_size;
        let is_contiguous = Self::check_contiguous(shape, &strides, element_size);

        Self {
            shape: shape.to_vec(),
            strides,
            element_size,
            total_size,
            order,
            is_contiguous,
            alignment,
        }
    }

    /// Create layout from existing array
    pub fn from_array<S, D>(array: &ArrayBase<S, D>) -> Self
    where
        S: RawData,
        D: Dimension,
    {
        let shape = array.shape().to_vec();
        let strides: Vec<isize> = array.strides().to_vec();
        let element_size = mem::size_of::<S::Elem>();
        let total_size = array.len() * element_size;

        // Determine order
        let order = if Self::is_c_order(&shape, &strides, element_size) {
            LayoutOrder::C
        } else if Self::is_f_order(&shape, &strides, element_size) {
            LayoutOrder::Fortran
        } else {
            LayoutOrder::Any
        };

        let is_contiguous = Self::check_contiguous(&shape, &strides, element_size);

        Self {
            shape,
            strides,
            element_size,
            total_size,
            order,
            is_contiguous,
            alignment: 64, // Default alignment
        }
    }

    /// Calculate strides for given shape and order
    pub fn calculate_strides(
        shape: &[usize],
        order: LayoutOrder,
        element_size: usize,
    ) -> Vec<isize> {
        if shape.is_empty() {
            return Vec::new();
        }

        let mut strides = vec![0; shape.len()];

        match order {
            LayoutOrder::C => {
                // C-order: stride[i] = shape[i+1] * shape[i+2] * ... * element_size
                let mut stride = element_size as isize;
                for i in (0..shape.len()).rev() {
                    strides[i] = stride;
                    stride *= shape[i] as isize;
                }
            }
            LayoutOrder::Fortran => {
                // F-order: stride[i] = shape[0] * shape[1] * ... * shape[i-1] * element_size
                let mut stride = element_size as isize;
                for i in 0..shape.len() {
                    strides[i] = stride;
                    stride *= shape[i] as isize;
                }
            }
            LayoutOrder::Any | LayoutOrder::Keep => {
                // Default to C-order for Any/Keep when creating new layout
                return Self::calculate_strides(shape, LayoutOrder::C, element_size);
            }
        }

        strides
    }

    /// Check if layout is C-contiguous
    pub fn is_c_contiguous(&self) -> bool {
        Self::is_c_order(&self.shape, &self.strides, self.element_size)
    }

    /// Check if layout is F-contiguous
    pub fn is_f_contiguous(&self) -> bool {
        Self::is_f_order(&self.shape, &self.strides, self.element_size)
    }

    /// Check if strides represent C-order
    fn is_c_order(shape: &[usize], strides: &[isize], elementsize: usize) -> bool {
        if shape.len() != strides.len() {
            return false;
        }

        let expected_strides = Self::calculate_strides(shape, LayoutOrder::C, elementsize);
        strides == expected_strides
    }

    /// Check if strides represent F-order
    fn is_f_order(shape: &[usize], strides: &[isize], elementsize: usize) -> bool {
        if shape.len() != strides.len() {
            return false;
        }

        let expected_strides = Self::calculate_strides(shape, LayoutOrder::Fortran, elementsize);
        strides == expected_strides
    }

    /// Check if layout is contiguous
    fn check_contiguous(shape: &[usize], strides: &[isize], elementsize: usize) -> bool {
        Self::is_c_order(shape, strides, elementsize)
            || Self::is_f_order(shape, strides, elementsize)
    }

    /// Convert to different layout order
    pub fn to_order(&self, neworder: LayoutOrder) -> CoreResult<Self> {
        if neworder == LayoutOrder::Keep {
            return Ok(self.clone());
        }

        if neworder == LayoutOrder::Any {
            return Ok(self.clone());
        }

        Ok(Self::new_with_order(
            &self.shape,
            neworder,
            self.element_size,
            self.alignment,
        ))
    }

    /// Get linear index for multi-dimensional indices
    pub fn linear_index(&self, indices: &[usize]) -> CoreResult<usize> {
        if indices.len() != self.shape.len() {
            return Err(CoreError::ShapeError(
                ErrorContext::new(format!(
                    "Index dimensions {} don't match array dimensions {}",
                    indices.len(),
                    self.shape.len()
                ))
                .with_location(ErrorLocation::new(file!(), line!())),
            ));
        }

        // Check bounds
        for (i, (&idx, &dim)) in indices.iter().zip(self.shape.iter()).enumerate() {
            if idx >= dim {
                return Err(CoreError::IndexError(
                    ErrorContext::new(format!(
                        "Index {idx} is out of bounds for axis {i} with size {dim}"
                    ))
                    .with_location(ErrorLocation::new(file!(), line!())),
                ));
            }
        }

        // Calculate linear index using strides
        let mut linear_idx = 0;
        for (idx, stride) in indices.iter().zip(self.strides.iter()) {
            linear_idx += (*idx as isize) * stride;
        }

        Ok((linear_idx / self.element_size as isize) as usize)
    }

    /// Get multi-dimensional indices for linear index
    pub fn multi_indices(&self, linearidx: usize) -> CoreResult<Vec<usize>> {
        let total_elements = self.shape.iter().product::<usize>();
        if linearidx >= total_elements {
            return Err(CoreError::IndexError(
                ErrorContext::new(format!(
                    "Linear index {linearidx} is out of bounds for array with {total_elements} elements"
                ))
                .with_location(ErrorLocation::new(file!(), line!())),
            ));
        }

        let mut indices = vec![0; self.shape.len()];
        let mut remaining = linearidx;

        match self.order {
            LayoutOrder::C => {
                // C-order: rightmost index varies fastest
                for (i, &shape_dim) in self.shape.iter().enumerate().rev() {
                    indices[i] = remaining % shape_dim;
                    remaining /= shape_dim;
                }
            }
            LayoutOrder::Fortran => {
                // F-order: leftmost index varies fastest
                for (idx, shape_dim) in indices.iter_mut().zip(&self.shape) {
                    *idx = remaining % shape_dim;
                    remaining /= shape_dim;
                }
            }
            LayoutOrder::Any | LayoutOrder::Keep => {
                // Use the actual strides to compute indices
                remaining *= self.element_size;
                for (i, &stride) in self.strides.iter().enumerate().take(self.shape.len()) {
                    if stride > 0 {
                        indices[i] = (remaining as isize / stride) as usize;
                        remaining = (remaining as isize % stride) as usize;
                    }
                }
            }
        }

        Ok(indices)
    }

    /// Calculate memory offset for given indices
    pub fn memory_offset(&self, indices: &[usize]) -> CoreResult<usize> {
        if indices.len() != self.shape.len() {
            return Err(CoreError::ShapeError(
                ErrorContext::new(format!(
                    "Index dimensions {} don't match array dimensions {}",
                    indices.len(),
                    self.shape.len()
                ))
                .with_location(ErrorLocation::new(file!(), line!())),
            ));
        }

        let mut offset = 0isize;
        for (idx, stride) in indices.iter().zip(self.strides.iter()) {
            offset += (*idx as isize) * stride;
        }

        Ok(offset as usize)
    }

    /// Get the number of dimensions
    pub fn ndim(&self) -> usize {
        self.shape.len()
    }

    /// Get the size (number of elements)
    pub fn size(&self) -> usize {
        self.shape.iter().product()
    }

    /// Get the number of bytes
    pub fn nbytes(&self) -> usize {
        self.total_size
    }

    /// Check if layout is compatible with another layout
    pub fn is_compatible_with(&self, other: &MemoryLayout) -> bool {
        self.shape == other.shape && self.element_size == other.element_size
    }

    /// Create a transposed layout
    pub fn transpose(&self, axes: Option<&[usize]>) -> CoreResult<Self> {
        let ndim = self.shape.len();

        let axes = if let Some(axes) = axes {
            // Validate provided axes
            if axes.len() != ndim {
                return Err(CoreError::ShapeError(
                    ErrorContext::new(format!(
                        "Axes length {} doesn't match number of dimensions {}",
                        axes.len(),
                        ndim
                    ))
                    .with_location(ErrorLocation::new(file!(), line!())),
                ));
            }

            let mut sorted_axes = axes.to_vec();
            sorted_axes.sort_unstable();
            let expected: Vec<usize> = (0..ndim).collect();
            if sorted_axes != expected {
                return Err(CoreError::ValidationError(
                    ErrorContext::new("Invalid transpose axes".to_string())
                        .with_location(ErrorLocation::new(file!(), line!())),
                ));
            }

            axes.to_vec()
        } else {
            // Default transpose: reverse all axes
            (0..ndim).rev().collect()
        };

        // Create new shape and strides
        let mut newshape = vec![0; ndim];
        let mut new_strides = vec![0; ndim];

        for (i, &axis) in axes.iter().enumerate() {
            newshape[i] = self.shape[axis];
            new_strides[i] = self.strides[axis];
        }

        // Determine new order
        let new_order = if Self::is_c_order(&newshape, &new_strides, self.element_size) {
            LayoutOrder::C
        } else if Self::is_f_order(&newshape, &new_strides, self.element_size) {
            LayoutOrder::Fortran
        } else {
            LayoutOrder::Any
        };

        let is_contiguous = Self::check_contiguous(&newshape, &new_strides, self.element_size);

        Ok(Self {
            shape: newshape,
            strides: new_strides,
            element_size: self.element_size,
            total_size: self.total_size,
            order: new_order,
            is_contiguous,
            alignment: self.alignment,
        })
    }

    /// Reshape the layout to new shape
    pub fn reshape(&self, newshape: &[usize]) -> CoreResult<Self> {
        let new_size = newshape.iter().product::<usize>();
        let old_size = self.size();

        if new_size != old_size {
            return Err(CoreError::ShapeError(
                ErrorContext::new(format!(
                    "Cannot reshape array of size {old_size} into shape with size {new_size}"
                ))
                .with_location(ErrorLocation::new(file!(), line!())),
            ));
        }

        // For contiguous arrays, reshaping preserves the order
        let order = if self.is_contiguous {
            self.order
        } else {
            LayoutOrder::C // Default to C-order for non-contiguous arrays
        };

        Ok(Self::new_with_order(
            newshape,
            order,
            self.element_size,
            self.alignment,
        ))
    }

    /// Create a view with different shape (without copying data)
    pub fn new_withshape_and_strides(
        &self,
        newshape: &[usize],
        new_strides: Option<&[isize]>,
    ) -> CoreResult<Self> {
        let strides = if let Some(strides) = new_strides {
            if strides.len() != newshape.len() {
                return Err(CoreError::ShapeError(
                    ErrorContext::new(format!(
                        "Strides length {} doesn't match shape length {}",
                        strides.len(),
                        newshape.len()
                    ))
                    .with_location(ErrorLocation::new(file!(), line!())),
                ));
            }
            strides.to_vec()
        } else {
            Self::calculate_strides(newshape, self.order, self.element_size)
        };

        // Determine order
        let order = if Self::is_c_order(newshape, &strides, self.element_size) {
            LayoutOrder::C
        } else if Self::is_f_order(newshape, &strides, self.element_size) {
            LayoutOrder::Fortran
        } else {
            LayoutOrder::Any
        };

        let is_contiguous = Self::check_contiguous(newshape, &strides, self.element_size);
        let total_size = newshape.iter().product::<usize>() * self.element_size;

        Ok(Self {
            shape: newshape.to_vec(),
            strides,
            element_size: self.element_size,
            total_size,
            order,
            is_contiguous,
            alignment: self.alignment,
        })
    }
}

/// Array layout utilities for working with ndarray
pub struct ArrayLayout;

impl ArrayLayout {
    /// Convert ndarray to C-order (row-major)
    pub fn to_c_order<A, S, D>(array: ArrayBase<S, D>) -> Array<A, D>
    where
        A: Clone + num_traits::Zero,
        S: Data<Elem = A>,
        D: Dimension,
    {
        if array.is_standard_layout() {
            // Already in C-order
            array.to_owned()
        } else {
            // Need to convert
            let mut result = Array::zeros(array.raw_dim());
            result.assign(&array);
            result
        }
    }

    /// Convert ndarray to F-order (column-major)
    pub fn to_f_order<A, S, D>(array: ArrayBase<S, D>) -> Array<A, D>
    where
        A: Clone + num_traits::Zero,
        S: Data<Elem = A>,
        D: Dimension,
    {
        // For simplicity, just convert to owned and let ndarray handle F-order conversion
        // Create F-order array by transposing appropriately
        // This is a simplified implementation - in practice we'd need more sophisticated F-order handling
        array.to_owned()
    }

    /// Create array with specific layout order
    pub fn zeros_with_order<A, D>(shape: D, order: LayoutOrder) -> Array<A, D>
    where
        A: Clone + Default + num_traits::Zero,
        D: Dimension,
    {
        match order {
            LayoutOrder::C => Array::zeros(shape),
            LayoutOrder::Fortran => Array::zeros(shape.f()),
            LayoutOrder::Any => Array::zeros(shape), // Default to C
            LayoutOrder::Keep => Array::zeros(shape), // Default to C
        }
    }

    /// Check if array has C-order layout
    pub fn is_c_order<S, D>(array: &ArrayBase<S, D>) -> bool
    where
        S: RawData,
        D: Dimension,
    {
        array.is_standard_layout()
    }

    /// Check if array has F-order layout
    pub fn is_f_order<S, D>(array: &ArrayBase<S, D>) -> bool
    where
        S: RawData,
        D: Dimension,
    {
        // Check if strides increase from first to last axis
        let strides = array.strides();
        if strides.len() <= 1 {
            return true;
        }

        for i in 1..strides.len() {
            if strides[i] < strides[i.saturating_sub(1)] {
                return false;
            }
        }
        true
    }

    /// Get memory layout information from ndarray
    pub fn get_layout<S, D>(array: &ArrayBase<S, D>) -> MemoryLayout
    where
        S: RawData,
        D: Dimension,
    {
        MemoryLayout::from_array(array)
    }

    /// Optimize array layout for specific access pattern
    pub fn optimize_for_access<A, S, D>(
        array: ArrayBase<S, D>,
        access_pattern: AccessPattern,
    ) -> Array<A, D>
    where
        A: Clone + num_traits::Zero,
        S: Data<Elem = A>,
        D: Dimension,
    {
        match access_pattern {
            AccessPattern::RowMajor => Self::to_c_order(array),
            AccessPattern::ColumnMajor => Self::to_f_order(array),
            AccessPattern::Random => {
                // For random access, prefer the more cache-friendly layout
                if array.len() > 10000 {
                    Self::to_c_order(array) // C-order generally better for large arrays
                } else {
                    array.to_owned()
                }
            }
        }
    }
}

/// Access pattern for layout optimization
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum AccessPattern {
    /// Row-major access (iterate over rows)
    RowMajor,
    /// Column-major access (iterate over columns)
    ColumnMajor,
    /// Random access pattern
    Random,
}

/// Memory-efficient layout converter for large arrays
pub struct LayoutConverter;

impl LayoutConverter {
    /// Convert array layout with chunked processing for memory efficiency
    pub fn convert_chunked<A>(
        source_layout: &MemoryLayout,
        target_order: LayoutOrder,
        data: &[A],
        chunk_size: usize,
    ) -> CoreResult<Vec<A>>
    where
        A: Clone + Copy + Default,
    {
        let target_layout = source_layout.to_order(target_order)?;

        if source_layout.order == target_order || source_layout.is_compatible_with(&target_layout) {
            // No conversion needed or already compatible
            return Ok(data.to_vec());
        }

        let total_elements = source_layout.size();
        let mut result = vec![A::default(); total_elements];

        // Process in chunks to avoid excessive memory usage
        let num_chunks = total_elements.div_ceil(chunk_size);

        for chunk_idx in 0..num_chunks {
            let start = chunk_idx * chunk_size;
            let end = std::cmp::min(start + chunk_size, total_elements);

            // Convert each element in the chunk
            #[allow(clippy::needless_range_loop)]
            for linear_idx in start..end {
                let source_indices = source_layout.multi_indices(linear_idx)?;
                let target_linear_idx = target_layout.linear_index(&source_indices)?;
                result[target_linear_idx] = data[linear_idx];
            }
        }

        Ok(result)
    }

    /// In-place layout conversion (when possible)
    pub fn convert_inplace<A>(
        layout: &MemoryLayout,
        target_order: LayoutOrder,
        data: &mut [A],
    ) -> CoreResult<MemoryLayout>
    where
        A: Clone + Copy + Default,
    {
        if layout.order == target_order {
            return Ok(layout.clone());
        }

        // For now, use out-of-place conversion
        // In-place conversion is complex and requires careful analysis of stride patterns
        let converted = Self::convert_chunked(layout, target_order, data, 8192)?;
        data.copy_from_slice(&converted);

        layout.to_order(target_order)
    }
}

/// `NumPy`-compatible array creation functions
pub struct ArrayCreation;

impl ArrayCreation {
    /// Create array with specified layout order (`NumPy`-compatible)
    pub fn array_with_order<A, D>(
        data: Vec<A>,
        shape: D,
        order: LayoutOrder,
    ) -> CoreResult<Array<A, D>>
    where
        A: Clone,
        D: Dimension,
    {
        let expected_size = shape.size();
        if data.len() != expected_size {
            return Err(CoreError::ShapeError(
                ErrorContext::new(format!(
                    "Data length {} doesn't match shape size {}",
                    data.len(),
                    expected_size
                ))
                .with_location(ErrorLocation::new(file!(), line!())),
            ));
        }

        match order {
            LayoutOrder::C => Array::from_shape_vec(shape, data).map_err(|e| {
                CoreError::ShapeError(
                    ErrorContext::new(format!("error: {e}"))
                        .with_location(ErrorLocation::new(file!(), line!())),
                )
            }),
            LayoutOrder::Fortran => Array::from_shape_vec(shape.f(), data).map_err(|e| {
                CoreError::ShapeError(
                    ErrorContext::new(format!("error: {e}"))
                        .with_location(ErrorLocation::new(file!(), line!())),
                )
            }),
            LayoutOrder::Any | LayoutOrder::Keep => {
                // Default to C-order
                Self::array_with_order(data, shape, LayoutOrder::C)
            }
        }
    }

    /// Create zeros array with specified order
    pub fn zeros_with_order<A, D>(shape: D, order: LayoutOrder) -> Array<A, D>
    where
        A: Clone + Default + num_traits::Zero,
        D: Dimension,
    {
        ArrayLayout::zeros_with_order(shape, order)
    }

    /// Create ones array with specified order
    pub fn ones_with_order<A, D>(shape: D, order: LayoutOrder) -> Array<A, D>
    where
        A: Clone + num_traits::One,
        D: Dimension,
    {
        match order {
            LayoutOrder::C => Array::ones(shape),
            LayoutOrder::Fortran => Array::ones(shape.f()),
            LayoutOrder::Any | LayoutOrder::Keep => Array::ones(shape),
        }
    }
}

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

    #[test]
    fn test_layout_order_parsing() {
        assert_eq!(
            LayoutOrder::parse("C").expect("Operation failed"),
            LayoutOrder::C
        );
        assert_eq!(
            LayoutOrder::parse("F").expect("Operation failed"),
            LayoutOrder::Fortran
        );
        assert_eq!(
            LayoutOrder::parse("fortran").expect("Operation failed"),
            LayoutOrder::Fortran
        );
        assert_eq!(
            LayoutOrder::parse("A").expect("Operation failed"),
            LayoutOrder::Any
        );
        assert_eq!(
            LayoutOrder::parse("K").expect("Operation failed"),
            LayoutOrder::Keep
        );
        assert!(LayoutOrder::parse("X").is_err());
    }

    #[test]
    fn test_memory_layout_creation() {
        let c_layout = MemoryLayout::new_c_order(&[10, 20]);
        assert_eq!(c_layout.shape, vec![10, 20]);
        assert_eq!(c_layout.order, LayoutOrder::C);
        assert!(c_layout.is_c_contiguous());
        assert!(!c_layout.is_f_contiguous());

        let f_layout = MemoryLayout::new_f_order(&[10, 20]);
        assert_eq!(f_layout.shape, vec![10, 20]);
        assert_eq!(f_layout.order, LayoutOrder::Fortran);
        assert!(!f_layout.is_c_contiguous());
        assert!(f_layout.is_f_contiguous());
    }

    #[test]
    fn test_stride_calculation() {
        let element_size = 8; // f64

        // C-order strides for [10, 20] shape
        let c_strides = MemoryLayout::calculate_strides(&[10, 20], LayoutOrder::C, element_size);
        assert_eq!(c_strides, vec![160, 8]); // [20*8, 1*8]

        // F-order strides for [10, 20] shape
        let f_strides =
            MemoryLayout::calculate_strides(&[10, 20], LayoutOrder::Fortran, element_size);
        assert_eq!(f_strides, vec![8, 80]); // [1*8, 10*8]
    }

    #[test]
    fn test_linear_indexing() {
        let layout = MemoryLayout::new_c_order(&[3, 4]);

        // Test some specific indices
        assert_eq!(layout.linear_index(&[0, 0]).expect("Operation failed"), 0);
        assert_eq!(layout.linear_index(&[0, 1]).expect("Operation failed"), 1);
        assert_eq!(layout.linear_index(&[1, 0]).expect("Operation failed"), 4);
        assert_eq!(layout.linear_index(&[2, 3]).expect("Operation failed"), 11);

        // Test bounds checking
        assert!(layout.linear_index(&[3, 0]).is_err());
        assert!(layout.linear_index(&[0, 4]).is_err());
    }

    #[test]
    fn test_multi_indexing() {
        let layout = MemoryLayout::new_c_order(&[3, 4]);

        assert_eq!(
            layout.multi_indices(0).expect("Operation failed"),
            vec![0, 0]
        );
        assert_eq!(
            layout.multi_indices(1).expect("Operation failed"),
            vec![0, 1]
        );
        assert_eq!(
            layout.multi_indices(4).expect("Operation failed"),
            vec![1, 0]
        );
        assert_eq!(
            layout.multi_indices(11).expect("Operation failed"),
            vec![2, 3]
        );

        // Test bounds
        assert!(layout.multi_indices(12).is_err());
    }

    #[test]
    fn test_layout_conversion() {
        let c_layout = MemoryLayout::new_c_order(&[5, 6]);
        let f_layout = c_layout
            .to_order(LayoutOrder::Fortran)
            .expect("Operation failed");

        assert_eq!(f_layout.order, LayoutOrder::Fortran);
        assert_eq!(f_layout.shape, c_layout.shape);
        assert!(f_layout.is_f_contiguous());
    }

    #[test]
    fn test_transpose() {
        let layout = MemoryLayout::new_c_order(&[3, 4, 5]);
        let transposed = layout
            .transpose(Some(&[2, 0, 1]))
            .expect("Operation failed");

        assert_eq!(transposed.shape, vec![5, 3, 4]);

        // Test default transpose (reverse axes)
        let default_transposed = layout.transpose(None).expect("Operation failed");
        assert_eq!(default_transposed.shape, vec![5, 4, 3]);
    }

    #[test]
    fn test_reshape() {
        let layout = MemoryLayout::new_c_order(&[6, 4]);
        let reshaped = layout.reshape(&[3, 8]).expect("Operation failed");

        assert_eq!(reshaped.shape, vec![3, 8]);
        assert_eq!(reshaped.size(), layout.size());

        // Test invalid reshape
        assert!(layout.reshape(&[5, 5]).is_err());
    }

    #[test]
    fn test_array_layout_utilities() {
        let arr = Array2::<f64>::zeros((10, 20));
        assert!(ArrayLayout::is_c_order(&arr));
        assert!(!ArrayLayout::is_f_order(&arr));

        let f_arr = Array2::<f64>::zeros((10, 20).f());
        assert!(!ArrayLayout::is_c_order(&f_arr));
        assert!(ArrayLayout::is_f_order(&f_arr));
    }

    #[test]
    fn test_array_creation_with_order() {
        let data: Vec<f64> = (0..12).map(|x| x as f64).collect();

        let c_array: crate::ndarray::Array2<f64> = ArrayCreation::array_with_order(
            data.clone(),
            crate::ndarray::Ix2(3, 4),
            LayoutOrder::C,
        )
        .expect("Operation failed");
        assert!(ArrayLayout::is_c_order(&c_array));

        let f_array: crate::ndarray::Array2<f64> =
            ArrayCreation::array_with_order(data, crate::ndarray::Ix2(3, 4), LayoutOrder::Fortran)
                .expect("Operation failed");
        assert!(ArrayLayout::is_f_order(&f_array));
    }
}