pub fn shape_of<T: MatrixValue>(matrix: &AbstractLiteral<T>) -> Option<MatrixShape<T::Dom>> {

    let mut res = shape_of_inner(matrix)?;

    res.strides.reverse();

    res.dims.reverse();

    res.idx_doms.reverse();

    Some(res)

}

fn shape_of_inner<T: MatrixValue>(matrix: &AbstractLiteral<T>) -> Option<MatrixShape<T::Dom>> {

    let AbstractLiteral::Matrix(elems, dom) = matrix else {

    let sz = elems.len();

    if sz == 0 {

    };

    let fst = elems[0].as_nested_matrix().and_then(shape_of_inner);

    for elem in elems.iter().skip(1) {

        debug_assert_eq!(

            fst,

            elem.as_nested_matrix().and_then(shape_of_inner),

    Some(match fst {

        None => MatrixShape {

            size: sz,

            strides: vec![1],

            dims: vec![sz],

            idx_doms: vec![dom.clone()],

        },

            res.strides.push(res.size);

            res.dims.push(sz);

            res.idx_doms.push(dom.clone());

            res.size = if res.size == 0 { sz } else { sz * res.size };

            res

}

pub fn shape_of_matrix_expr(expr: &Expr) -> Option<MatrixShape<DomainPtr>> {

    match expr {

        _ => None,

}

pub fn shape_of_dom(

    matrix_dom_gd: &GroundDomain,

) -> Result<MatrixShape<Moo<GroundDomain>>, DomainOpError> {

    let GroundDomain::Matrix(_, idx_doms, _) = matrix_dom_gd else {

    let len = idx_doms.len();

    let mut strides = VecDeque::with_capacity(len);

    let mut dimensions = VecDeque::with_capacity(len);

    let mut size: usize = 1;

    for gd in idx_doms.iter().rev() {

        let gd_sz = gd.len_usize()?;

        strides.push_front(size);

        dimensions.push_front(gd_sz);

        size = size.checked_mul(gd_sz).ok_or(DomainOpError::TooLarge)?;

    Ok(MatrixShape {

        size,

        dims: dimensions.into(),

        strides: strides.into(),

        idx_doms: idx_doms.clone(),

    })

}

pub fn partial_flatten<T: MatrixValue>(n: usize, matrix: AbstractLiteral<T>) -> AbstractLiteral<T> {

    if n == 0 {

        return matrix;

    }

        n > 0 && n < shape.dims.len(),

    let new_strides = Vec::from(&shape.strides[n..]);

    let flattened = flatten_owned(matrix).collect_vec();

    let res = unflatten_list(&flattened, &new_strides);

    res.into_nested_matrix()

}

pub fn flatten<T: MatrixValue>(matrix: &AbstractLiteral<T>) -> impl Iterator<Item = &T> {

    let AbstractLiteral::Matrix(elems, _) = matrix else {

    flatten_inner(elems)

}

fn flatten_inner<'a, T: MatrixValue>(elems: &'a [T]) -> impl Iterator<Item = &'a T> {

    elems.iter().flat_map(|elem| {

        if let Some(m) = elem.as_nested_matrix() {

            Box::new(flatten(m)) as Box<dyn Iterator<Item = &'a T>>

            Box::new(std::iter::once(elem)) as Box<dyn Iterator<Item = &'a T>>

    })

}

pub fn flatten_owned<T: MatrixValue>(matrix: AbstractLiteral<T>) -> impl Iterator<Item = T> {

    let AbstractLiteral::Matrix(elems, _) = matrix else {

    flatten_owned_inner(elems)

}

fn flatten_owned_inner<T: MatrixValue>(elems: Vec<T>) -> impl Iterator<Item = T> {

    elems

        .into_iter()

        .flat_map(|elem| match elem.into_nested_matrix() {

            Ok(m) => Box::new(flatten_owned(m)) as Box<dyn Iterator<Item = T>>,

            Err(leaf) => Box::new(std::iter::once(leaf)) as Box<dyn Iterator<Item = T>>,

        })

}

pub fn unflatten_matrix<T: MatrixValue>(

    elems: &[T],

    index_domains: &[T::Dom],

    strides: &[usize],

) -> T {

    let dom = index_domains.first().expect("no index domains").clone();

    let stride = *strides.first().expect("no strides");

    if index_domains.len() == 1 {

        return T::from(AbstractLiteral::Matrix(Vec::from(elems), dom));

    }

    let mut inners = Vec::<T>::with_capacity(stride);

    let mut i_start: usize = 0;

    while i_start < elems.len() {

        let next = i_start + stride;

        let elem = unflatten_matrix(&elems[i_start..next], &index_domains[1..], &strides[1..]);

        inners.push(elem);

        i_start = next;

    T::from(AbstractLiteral::Matrix(inners, dom))

}

pub fn unflatten_list<T: MatrixValue>(elems: &[T], strides: &[usize]) -> T {

    let stride = *strides.first().expect("no strides");

    if strides.len() == 1 {

        return AbstractLiteral::matrix_implied_indices(Vec::from(elems)).into();

    }

    let mut inners = Vec::<T>::with_capacity(stride);

    let mut i_start: usize = 0;

    while i_start < elems.len() {

        let next = i_start + stride;

        let elem = unflatten_list(&elems[i_start..next], &strides[1..]);

        inners.push(elem);

        i_start = next;

    AbstractLiteral::matrix_implied_indices(inners).into()

}

pub fn index_domains<T: MatrixValue>(matrix: &AbstractLiteral<T>) -> Vec<T::Dom> {

    shape_of(matrix)

}

pub fn try_enumerate_indices(

    index_domains: Vec<Moo<GroundDomain>>,

) -> Result<impl Iterator<Item = Vec<Literal>>, DomainOpError> {

    let domains = index_domains

        .into_iter()

        .map(|x| x.values().map(|values| values.collect_vec()))

    Ok(domains.into_iter().multi_cartesian_product())

}

pub fn enumerate_indices(

    index_domains: Vec<Moo<GroundDomain>>,

) -> impl Iterator<Item = Vec<Literal>> {

    try_enumerate_indices(index_domains).expect("index domain should be enumerable with .values()")

}

pub fn num_elements(index_domains: &[Moo<GroundDomain>]) -> Result<u64, DomainOpError> {

    let idx_dom_lengths = index_domains

        .iter()

        .map(|d| d.length())

    Ok(idx_dom_lengths.iter().product())

}

pub fn flatten_enumerate(

    matrix: AbstractLiteral<Literal>,

) -> impl Iterator<Item = (Vec<Literal>, Literal)> {

    let index_domains: Vec<Moo<GroundDomain>> = shape

        .idx_doms

        .into_iter()

        .zip(shape.dims)

        .map(|(domain, len)| bound_index_domain_from_length(domain, len))

    izip!(enumerate_indices(index_domains), flatten_owned(matrix))

}

pub fn enumerate_index_union_indices(

    a_domains: &[Moo<GroundDomain>],

    b_domains: &[Moo<GroundDomain>],

) -> Result<impl Iterator<Item = Vec<Literal>>, DomainOpError> {

    if a_domains.len() != b_domains.len() {

    }

    let idx_domains: Result<Vec<_>, _> = a_domains

        .iter()

        .zip(b_domains.iter())

        .map(|(a, b)| a.union(b))

    let idx_domains = idx_domains?.into_iter().map(Moo::new).collect();

    try_enumerate_indices(idx_domains)

}

pub fn flat_index_to_full_index(index_domains: &[Moo<GroundDomain>], index: u64) -> Vec<Literal> {

    let mut remaining = index;

    let mut multipliers = vec![1; index_domains.len()];

    for i in (1..index_domains.len()).rev() {

        multipliers[i - 1] = multipliers[i] * index_domains[i].as_ref().length().unwrap();

    let mut coords = Vec::new();

    for m in multipliers.iter() {

        coords.push(((remaining / m + 1) as i32).into());

        remaining %= *m;

    coords

}

pub fn bound_index_domains_of_expr(expr: &Expr) -> Option<Vec<Moo<GroundDomain>>> {

    let dom = expr.domain_of().and_then(|dom| dom.resolve().ok())?;

    let GroundDomain::Matrix(_, index_domains, _) = dom.as_ref() else {

        return None;

    let Some(dimension_lengths) = expr_matrix_dimension_lengths(expr) else {

        return Some(index_domains.clone());

}

pub fn safe_index_optimised(m: Expr, idx: Literal) -> Option<Expr> {

    match m {

            let mut idxs = idxs;

            let (slice_idx, _) = idxs.iter().find_position(|opt| opt.is_none())?;

            let _ = idxs[slice_idx].replace(idx.into());

            let Some(idxs) = idxs.into_iter().collect::<Option<Vec<_>>>() else {

            Some(Expr::SafeIndex(Metadata::new(), mat, idxs))

            let Literal::Int(index) = idx else {

            let index_domains = bound_index_domains_of_expr(inner.as_ref())?;

            if index_domains.iter().any(|domain| domain.length().is_err()) {

            }

            let flat_index = flat_index_to_full_index(&index_domains, (index - 1) as u64);

            let flat_index: Vec<Expr> = flat_index.into_iter().map(Into::into).collect();

            Some(Expr::SafeIndex(Metadata::new(), inner, flat_index))

        _ => Some(Expr::SafeIndex(

            Metadata::new(),

            Moo::new(m),

            vec![idx.into()],

        )),

}

fn expr_matrix_dimension_lengths(expr: &Expr) -> Option<Vec<usize>> {

    Some(shape_of_matrix_expr(expr)?.dims)

}

fn bound_index_domain_from_length(mut domain: Moo<GroundDomain>, len: usize) -> Moo<GroundDomain> {

    match Moo::make_mut(&mut domain) {

        GroundDomain::Int(ranges, _) if ranges.len() == 1 && len > 0 => {

            ranges[0] = match &ranges[0] {

                Range::UnboundedR(start) => Range::Bounded(*start, start + (len as i32 - 1)),

                Range::Bounded(low, high) => Range::Bounded(*low, *high),

                Range::Single(value) => Range::Single(*value),

            domain

        _ => domain,

}

    fn as_nested_matrix(&self) -> Option<&AbstractLiteral<Literal>> {

        match self {

            Literal::AbstractLiteral(m @ AbstractLiteral::Matrix(..)) => Some(m),

            _ => None,

    }

    fn into_nested_matrix(self) -> Result<AbstractLiteral<Literal>, Self> {

        match self {

            Literal::AbstractLiteral(m @ AbstractLiteral::Matrix(..)) => Ok(m),

            other => Err(other),

    }