text: import operator
  from collections import defaultdict, Counter
  from functools import reduce
  import itertools
  from itertools import accumulate
  from typing import Optional, List, Dict as tDict, Tuple as tTuple
  
  import typing
  
  from sympy.core.basic import Basic
  from sympy.core.containers import Tuple
  from sympy.core.expr import Expr
  from sympy.core.function import (Function, Lambda)
  from sympy.core.mul import Mul
  from sympy.core.singleton import S
  from sympy.core.sorting import default_sort_key
  from sympy.core.symbol import (Dummy, Symbol)
  from sympy.matrices.expressions.diagonal import diagonalize_vector
  from sympy.matrices.expressions.matexpr import MatrixExpr
  from sympy.matrices.expressions.special import ZeroMatrix
  from sympy.tensor.array.arrayop import (permutedims, tensorcontraction, tensordiagonal, tensorproduct)
  from sympy.tensor.array.dense_ndim_array import ImmutableDenseNDimArray
  from sympy.tensor.array.ndim_array import NDimArray
  from sympy.tensor.indexed import (Indexed, IndexedBase)
  from sympy.matrices.expressions.matexpr import MatrixElement
  from sympy.tensor.array.expressions.utils import _apply_recursively_over_nested_lists, _sort_contraction_indices, \
      _get_mapping_from_subranks, _build_push_indices_up_func_transformation, _get_contraction_links, \
      _build_push_indices_down_func_transformation
  from sympy.combinatorics import Permutation
  from sympy.combinatorics.permutations import _af_invert
  from sympy.core.sympify import _sympify
  
  
  class _ArrayExpr(Expr):
      pass
  
  
  class ArraySymbol(_ArrayExpr):
      """
      Symbol representing an array expression
      """
  
      def __new__(cls, symbol, shape: typing.Iterable) -> "ArraySymbol":
          if isinstance(symbol, str):
              symbol = Symbol(symbol)
          # symbol = _sympify(symbol)
          shape = Tuple(*map(_sympify, shape))
          obj = Expr.__new__(cls, symbol, shape)
          return obj
  
      @property
      def name(self):
          return self._args[0]
  
      @property
      def shape(self):
          return self._args[1]
  
      def __getitem__(self, item):
          return ArrayElement(self, item)
  
      def as_explicit(self):
          if not all(i.is_Integer for i in self.shape):
              raise ValueError("cannot express explicit array with symbolic shape")
          data = [self[i] for i in itertools.product(*[range(j) for j in self.shape])]
          return ImmutableDenseNDimArray(data).reshape(*self.shape)
  
  
  class ArrayElement(_ArrayExpr):
      """
      An element of an array.
      """
      def __new__(cls, name, indices):
          if isinstance(name, str):
              name = Symbol(name)
          name = _sympify(name)
          indices = _sympify(tuple(indices))
          if hasattr(name, "shape"):
              if any((i >= s) == True for i, s in zip(indices, name.shape)):
                  raise ValueError("shape is out of bounds")
          if any((i < 0) == True for i in indices):
              raise ValueError("shape contains negative values")
          obj = Expr.__new__(cls, name, indices)
          return obj
  
      @property
      def name(self):
          return self._args[0]
  
      @property
      def indices(self):
          return self._args[1]
  
  
  class ZeroArray(_ArrayExpr):
      """
      Symbolic array of zeros. Equivalent to ``ZeroMatrix`` for matrices.
      """
  
      def __new__(cls, *shape):
          if len(shape) == 0:
              return S.Zero
          shape = map(_sympify, shape)
          obj = Expr.__new__(cls, *shape)
          return obj
  
      @property
      def shape(self):
          return self._args
  
      def as_explicit(self):
          if not all(i.is_Integer for i in self.shape):
              raise ValueError("Cannot return explicit form for symbolic shape.")
          return ImmutableDenseNDimArray.zeros(*self.shape)
  
  
  class OneArray(_ArrayExpr):
      """
      Symbolic array of ones.
      """
  
      def __new__(cls, *shape):
          if len(shape) == 0:
              return S.One
          shape = map(_sympify, shape)
          obj = Expr.__new__(cls, *shape)
          return obj
  
      @property
      def shape(self):
          return self._args
  
      def as_explicit(self):
          if not all(i.is_Integer for i in self.shape):
              raise ValueError("Cannot return explicit form for symbolic shape.")
          return ImmutableDenseNDimArray([S.One for i in range(reduce(operator.mul, self.shape))]).reshape(*self.shape)
  
  
  class _CodegenArrayAbstract(Basic):
  
      @property
      def subranks(self):
          """
          Returns the ranks of the objects in the uppermost tensor product inside
          the current object.  In case no tensor products are contained, return
          the atomic ranks.
  
          Examples
          ========
  
          >>> from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct, ArrayContraction
          >>> from sympy import MatrixSymbol
          >>> M = MatrixSymbol("M", 3, 3)
          >>> N = MatrixSymbol("N", 3, 3)
          >>> P = MatrixSymbol("P", 3, 3)
  
          Important: do not confuse the rank of the matrix with the rank of an array.
  
          >>> tp = ArrayTensorProduct(M, N, P)
          >>> tp.subranks
          [2, 2, 2]
  
          >>> co = ArrayContraction(tp, (1, 2), (3, 4))
          >>> co.subranks
          [2, 2, 2]
          """
          return self._subranks[:]
  
      def subrank(self):
          """
          The sum of ``subranks``.
          """
          return sum(self.subranks)
  
      @property
      def shape(self):
          return self._shape
  
  
  class ArrayTensorProduct(_CodegenArrayAbstract):
      r"""
      Class to represent the tensor product of array-like objects.
      """
  
      def __new__(cls, *args, **kwargs):
          args = [_sympify(arg) for arg in args]
  
          normalize = kwargs.pop("normalize", True)
  
          ranks = [get_rank(arg) for arg in args]
  
          obj = Basic.__new__(cls, *args)
          obj._subranks = ranks
          shapes = [get_shape(i) for i in args]
  
          if any(i is None for i in shapes):
              obj._shape = None
          else:
              obj._shape = tuple(j for i in shapes for j in i)
          if normalize:
              return obj._normalize()
          return obj
  
      def _normalize(self):
          args = self.args
          args = self._flatten(args)
  
          ranks = [get_rank(arg) for arg in args]
  
          # Check if there are nested permutation and lift them up:
          permutation_cycles = []
          for i, arg in enumerate(args):
              if not isinstance(arg, PermuteDims):
                  continue
              permutation_cycles.extend([[k + sum(ranks[:i]) for k in j] for j in arg.permutation.cyclic_form])
              args[i] = arg.expr
          if permutation_cycles:
              return PermuteDims(ArrayTensorProduct(*args), Permutation(sum(ranks)-1)*Permutation(permutation_cycles))
  
          if len(args) == 1:
              return args[0]
  
          # If any object is a ZeroArray, return a ZeroArray:
          if any(isinstance(arg, (ZeroArray, ZeroMatrix)) for arg in args):
              shapes = reduce(operator.add, [get_shape(i) for i in args], ())
              return ZeroArray(*shapes)
  
          # If there are contraction objects inside, transform the whole
          # expression into `ArrayContraction`:
          contractions = {i: arg for i, arg in enumerate(args) if isinstance(arg, ArrayContraction)}
          if contractions:
              ranks = [_get_subrank(arg) if isinstance(arg, ArrayContraction) else get_rank(arg) for arg in args]
              cumulative_ranks = list(accumulate([0] + ranks))[:-1]
              tp = self.func(*[arg.expr if isinstance(arg, ArrayContraction) else arg for arg in args])
              contraction_indices = [tuple(cumulative_ranks[i] + k for k in j) for i, arg in contractions.items() for j in arg.contraction_indices]
              return ArrayContraction(tp, *contraction_indices)
  
          diagonals = {i: arg for i, arg in enumerate(args) if isinstance(arg, ArrayDiagonal)}
          if diagonals:
              inverse_permutation = []
              last_perm = []
              ranks = [get_rank(arg) for arg in args]
              cumulative_ranks = list(accumulate([0] + ranks))[:-1]
              for i, arg in enumerate(args):
                  if isinstance(arg, ArrayDiagonal):
                      i1 = get_rank(arg) - len(arg.diagonal_indices)
                      i2 = len(arg.diagonal_indices)
                      inverse_permutation.extend([cumulative_ranks[i] + j for j in range(i1)])
                      last_perm.extend([cumulative_ranks[i] + j for j in range(i1, i1 + i2)])
                  else:
                      inverse_permutation.extend([cumulative_ranks[i] + j for j in range(get_rank(arg))])
              inverse_permutation.extend(last_perm)
              tp = self.func(*[arg.expr if isinstance(arg, ArrayDiagonal) else arg for arg in args])
              ranks2 = [_get_subrank(arg) if isinstance(arg, ArrayDiagonal) else get_rank(arg) for arg in args]
              cumulative_ranks2 = list(accumulate([0] + ranks2))[:-1]
              diagonal_indices = [tuple(cumulative_ranks2[i] + k for k in j) for i, arg in diagonals.items() for j in arg.diagonal_indices]
              return PermuteDims(ArrayDiagonal(tp, *diagonal_indices), _af_invert(inverse_permutation))
  
          return self.func(*args, normalize=False)
  
      @classmethod
      def _flatten(cls, args):
          args = [i for arg in args for i in (arg.args if isinstance(arg, cls) else [arg])]
          return args
  
      def as_explicit(self):
          return tensorproduct(*[arg.as_explicit() if hasattr(arg, "as_explicit") else arg for arg in self.args])
  
  
  class ArrayAdd(_CodegenArrayAbstract):
      r"""
      Class for elementwise array additions.
      """
  
      def __new__(cls, *args, **kwargs):
          args = [_sympify(arg) for arg in args]
          ranks = [get_rank(arg) for arg in args]
          ranks = list(set(ranks))
          if len(ranks) != 1:
              raise ValueError("summing arrays of different ranks")
          shapes = [arg.shape for arg in args]
          if len({i for i in shapes if i is not None}) > 1:
              raise ValueError("mismatching shapes in addition")
  
          normalize = kwargs.pop("normalize", True)
  
          obj = Basic.__new__(cls, *args)
          obj._subranks = ranks
          if any(i is None for i in shapes):
              obj._shape = None
          else:
              obj._shape = shapes[0]
          if normalize:
              return obj._normalize()
          return obj
  
      def _normalize(self):
          args = self.args
  
          # Flatten:
          args = self._flatten_args(args)
  
          shapes = [get_shape(arg) for arg in args]
          args = [arg for arg in args if not isinstance(arg, (ZeroArray, ZeroMatrix))]
          if len(args) == 0:
              if any(i for i in shapes if i is None):
                  raise NotImplementedError("cannot handle addition of ZeroMatrix/ZeroArray and undefined shape object")
              return ZeroArray(*shapes[0])
          elif len(args) == 1:
              return args[0]
          return self.func(*args, normalize=False)
  
      @classmethod
      def _flatten_args(cls, args):
          new_args = []
          for arg in args:
              if isinstance(arg, ArrayAdd):
                  new_args.extend(arg.args)
              else:
                  new_args.append(arg)
          return new_args
  
      def as_explicit(self):
          return reduce(operator.add, [arg.as_explicit() for arg in self.args])
  
  
  class PermuteDims(_CodegenArrayAbstract):
      r"""
      Class to represent permutation of axes of arrays.
  
      Examples
      ========
  
      >>> from sympy.tensor.array.expressions.array_expressions import PermuteDims
      >>> from sympy import MatrixSymbol
      >>> M = MatrixSymbol("M", 3, 3)
      >>> cg = PermuteDims(M, [1, 0])
  
      The object ``cg`` represents the transposition of ``M``, as the permutation
      ``[1, 0]`` will act on its indices by switching them:
  
      `M_{ij} \Rightarrow M_{ji}`
  
      This is evident when transforming back to matrix form:
  
      >>> from sympy.tensor.array.expressions.conv_array_to_matrix import convert_array_to_matrix
      >>> convert_array_to_matrix(cg)
      M.T
  
      >>> N = MatrixSymbol("N", 3, 2)
      >>> cg = PermuteDims(N, [1, 0])
      >>> cg.shape
      (2, 3)
  
      Permutations of tensor products are simplified in order to achieve a
      standard form:
  
      >>> from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
      >>> M = MatrixSymbol("M", 4, 5)
      >>> tp = ArrayTensorProduct(M, N)
      >>> tp.shape
      (4, 5, 3, 2)
      >>> perm1 = PermuteDims(tp, [2, 3, 1, 0])
  
      The args ``(M, N)`` have been sorted and the permutation has been
      simplified, the expression is equivalent:
  
      >>> perm1.expr.args
      (N, M)
      >>> perm1.shape
      (3, 2, 5, 4)
      >>> perm1.permutation
      (2 3)
  
      The permutation in its array form has been simplified from
      ``[2, 3, 1, 0]`` to ``[0, 1, 3, 2]``, as the arguments of the tensor
      product `M` and `N` have been switched:
  
      >>> perm1.permutation.array_form
      [0, 1, 3, 2]
  
      We can nest a second permutation:
  
      >>> perm2 = PermuteDims(perm1, [1, 0, 2, 3])
      >>> perm2.shape
      (2, 3, 5, 4)
      >>> perm2.permutation.array_form
      [1, 0, 3, 2]
      """
  
      def __new__(cls, expr, permutation, **kwargs):
          from sympy.combinatorics import Permutation
          expr = _sympify(expr)
          permutation = Permutation(permutation)
          permutation_size = permutation.size
          expr_rank = get_rank(expr)
          if permutation_size != expr_rank:
              raise ValueError("Permutation size must be the length of the shape of expr")
  
          normalize = kwargs.pop("normalize", True)
  
          obj = Basic.__new__(cls, expr, permutation)
          obj._subranks = [get_rank(expr)]
          shape = get_shape(expr)
          if shape is None:
              obj._shape = None
          else:
              obj._shape = tuple(shape[permutation(i)] for i in range(len(shape)))
          if normalize:
              return obj._normalize()
          return obj
  
      def _normalize(self):
          expr = self.expr
          permutation = self.permutation
          if isinstance(expr, PermuteDims):
              subexpr = expr.expr
              subperm = expr.permutation
              permutation = permutation * subperm
              expr = subexpr
          if isinstance(expr, ArrayContraction):
              expr, permutation = self._PermuteDims_denestarg_ArrayContraction(expr, permutation)
          if isinstance(expr, ArrayTensorProduct):
              expr, permutation = self._PermuteDims_denestarg_ArrayTensorProduct(expr, permutation)
          if isinstance(expr, (ZeroArray, ZeroMatrix)):
              return ZeroArray(*[expr.shape[i] for i in permutation.array_form])
          plist = permutation.array_form
          if plist == sorted(plist):
              return expr
          return self.func(expr, permutation, normalize=False)
  
      @property
      def expr(self):
          return self.args[0]
  
      @property
      def permutation(self):
          return self.args[1]
  
      @classmethod
      def _PermuteDims_denestarg_ArrayTensorProduct(cls, expr, permutation):
          # Get the permutation in its image-form:
          perm_image_form = _af_invert(permutation.array_form)
          args = list(expr.args)
          # Starting index global position for every arg:
          cumul = list(accumulate([0] + expr.subranks))
          # Split `perm_image_form` into a list of list corresponding to the indices
          # of every argument:
          perm_image_form_in_components = [perm_image_form[cumul[i]:cumul[i+1]] for i in range(len(args))]
          # Create an index, target-position-key array:
          ps = [(i, sorted(comp)) for i, comp in enumerate(perm_image_form_in_components)]
          # Sort the array according to the target-position-key:
          # In this way, we define a canonical way to sort the arguments according
          # to the permutation.
          ps.sort(key=lambda x: x[1])
          # Read the inverse-permutation (i.e. image-form) of the args:
          perm_args_image_form = [i[0] for i in ps]
          # Apply the args-permutation to the `args`:
          args_sorted = [args[i] for i in perm_args_image_form]
          # Apply the args-permutation to the array-form of the permutation of the axes (of `expr`):
          perm_image_form_sorted_args = [perm_image_form_in_components[i] for i in perm_args_image_form]
          new_permutation = Permutation(_af_invert([j for i in perm_image_form_sorted_args for j in i]))
          return ArrayTensorProduct(*args_sorted), new_permutation
  
      @classmethod
      def _PermuteDims_denestarg_ArrayContraction(cls, expr, permutation):
          if not isinstance(expr, ArrayContraction):
              return expr, permutation
          if not isinstance(expr.expr, ArrayTensorProduct):
              return expr, permutation
          args = expr.expr.args
          subranks = [get_rank(arg) for arg in expr.expr.args]
  
          contraction_indices = expr.contraction_indices
          contraction_indices_flat = [j for i in contraction_indices for j in i]
          cumul = list(accumulate([0] + subranks))
  
          # Spread the permutation in its array form across the args in the corresponding
          # tensor-product arguments with free indices:
          permutation_array_blocks_up = []
          image_form = _af_invert(permutation.array_form)
          counter = 0
          for i, e in enumerate(subranks):
              current = []
              for j in range(cumul[i], cumul[i+1]):
                  if j in contraction_indices_flat:
                      continue
                  current.append(image_form[counter])
                  counter += 1
              permutation_array_blocks_up.append(current)
  
          # Get the map of axis repositioning for every argument of tensor-product:
          index_blocks = [[j for j in range(cumul[i], cumul[i+1])] for i, e in enumerate(expr.subranks)]
          index_blocks_up = expr._push_indices_up(expr.contraction_indices, index_blocks)
          inverse_permutation = permutation**(-1)
          index_blocks_up_permuted = [[inverse_permutation(j) for j in i if j is not None] for i in index_blocks_up]
  
          # Sorting key is a list of tuple, first element is the index of `args`, second element of
          # the tuple is the sorting key to sort `args` of the tensor product:
          sorting_keys = list(enumerate(index_blocks_up_permuted))
          sorting_keys.sort(key=lambda x: x[1])
  
          # Now we can get the permutation acting on the args in its image-form:
          new_perm_image_form = [i[0] for i in sorting_keys]
          # Apply the args-level permutation to various elements:
          new_index_blocks = [index_blocks[i] for i in new_perm_image_form]
          new_index_perm_array_form = _af_invert([j for i in new_index_blocks for j in i])
          new_args = [args[i] for i in new_perm_image_form]
          new_contraction_indices = [tuple(new_index_perm_array_form[j] for j in i) for i in contraction_indices]
          new_expr = ArrayContraction(ArrayTensorProduct(*new_args), *new_contraction_indices)
          new_permutation = Permutation(_af_invert([j for i in [permutation_array_blocks_up[k] for k in new_perm_image_form] for j in i]))
          return new_expr, new_permutation
  
      @classmethod
      def _check_permutation_mapping(cls, expr, permutation):
          subranks = expr.subranks
          index2arg = [i for i, arg in enumerate(expr.args) for j in range(expr.subranks[i])]
          permuted_indices = [permutation(i) for i in range(expr.subrank())]
          new_args = list(expr.args)
          arg_candidate_index = index2arg[permuted_indices[0]]
          current_indices = []
          new_permutation = []
          inserted_arg_cand_indices = set([])
          for i, idx in enumerate(permuted_indices):
              if index2arg[idx] != arg_candidate_index:
                  new_permutation.extend(current_indices)
                  current_indices = []
                  arg_candidate_index = index2arg[idx]
              current_indices.append(idx)
              arg_candidate_rank = subranks[arg_candidate_index]
              if len(current_indices) == arg_candidate_rank:
                  new_permutation.extend(sorted(current_indices))
                  local_current_indices = [j - min(current_indices) for j in current_indices]
                  i1 = index2arg[i]
                  new_args[i1] = PermuteDims(new_args[i1], Permutation(local_current_indices))
                  inserted_arg_cand_indices.add(arg_candidate_index)
                  current_indices = []
          new_permutation.extend(current_indices)
  
          # TODO: swap args positions in order to simplify the expression:
          # TODO: this should be in a function
          args_positions = list(range(len(new_args)))
          # Get possible shifts:
          maps = {}
          cumulative_subranks = [0] + list(accumulate(subranks))
          for i in range(0, len(subranks)):
              s = set([index2arg[new_permutation[j]] for j in range(cumulative_subranks[i], cumulative_subranks[i+1])])
              if len(s) != 1:
                  continue
              elem = next(iter(s))
              if i != elem:
                  maps[i] = elem
  
          # Find cycles in the map:
          lines = []
          current_line = []
          while maps:
              if len(current_line) == 0:
                  k, v = maps.popitem()
                  current_line.append(k)
              else:
                  k = current_line[-1]
                  if k not in maps:
                      current_line = []
                      continue
                  v = maps.pop(k)
              if v in current_line:
                  lines.append(current_line)
                  current_line = []
                  continue
              current_line.append(v)
          for line in lines:
              for i, e in enumerate(line):
                  args_positions[line[(i + 1) % len(line)]] = e
  
          # TODO: function in order to permute the args:
          permutation_blocks = [[new_permutation[cumulative_subranks[i] + j] for j in range(e)] for i, e in enumerate(subranks)]
          new_args = [new_args[i] for i in args_positions]
          new_permutation_blocks = [permutation_blocks[i] for i in args_positions]
          new_permutation2 = [j for i in new_permutation_blocks for j in i]
          return ArrayTensorProduct(*new_args), Permutation(new_permutation2)  # **(-1)
  
      @classmethod
      def _check_if_there_are_closed_cycles(cls, expr, permutation):
          args = list(expr.args)
          subranks = expr.subranks
          cyclic_form = permutation.cyclic_form
          cumulative_subranks = [0] + list(accumulate(subranks))
          cyclic_min = [min(i) for i in cyclic_form]
          cyclic_max = [max(i) for i in cyclic_form]
          cyclic_keep = []
          for i, cycle in enumerate(cyclic_form):
              flag = True
              for j in range(0, len(cumulative_subranks) - 1):
                  if cyclic_min[i] >= cumulative_subranks[j] and cyclic_max[i] < cumulative_subranks[j+1]:
                      # Found a sinkable cycle.
                      args[j] = PermuteDims(args[j], Permutation([[k - cumulative_subranks[j] for k in cyclic_form[i]]]))
                      flag = False
                      break
              if flag:
                  cyclic_keep.append(cyclic_form[i])
          return ArrayTensorProduct(*args), Permutation(cyclic_keep, size=permutation.size)
  
      def nest_permutation(self):
          r"""
          DEPRECATED.
          """
          ret = self._nest_permutation(self.expr, self.permutation)
          if ret is None:
              return self
          return ret
  
      @classmethod
      def _nest_permutation(cls, expr, permutation):
          if isinstance(expr, ArrayTensorProduct):
              return PermuteDims(*cls._check_if_there_are_closed_cycles(expr, permutation))
          elif isinstance(expr, ArrayContraction):
              # Invert tree hierarchy: put the contraction above.
              cycles = permutation.cyclic_form
              newcycles = ArrayContraction._convert_outer_indices_to_inner_indices(expr, *cycles)
              newpermutation = Permutation(newcycles)
              new_contr_indices = [tuple(newpermutation(j) for j in i) for i in expr.contraction_indices]
              return ArrayContraction(PermuteDims(expr.expr, newpermutation), *new_contr_indices)
          elif isinstance(expr, ArrayAdd):
              return ArrayAdd(*[PermuteDims(arg, permutation) for arg in expr.args])
          return None
  
      def as_explicit(self):
          return permutedims(self.expr.as_explicit(), self.permutation)
  
  
  class ArrayDiagonal(_CodegenArrayAbstract):
      r"""
      Class to represent the diagonal operator.
  
      Explanation
      ===========
  
      In a 2-dimensional array it returns the diagonal, this looks like the
      operation:
  
      `A_{ij} \rightarrow A_{ii}`
  
      The diagonal over axes 1 and 2 (the second and third) of the tensor product
      of two 2-dimensional arrays `A \otimes B` is
  
      `\Big[ A_{ab} B_{cd} \Big]_{abcd} \rightarrow \Big[ A_{ai} B_{id} \Big]_{adi}`
  
      In this last example the array expression has been reduced from
      4-dimensional to 3-dimensional. Notice that no contraction has occurred,
      rather there is a new index `i` for the diagonal, contraction would have
      reduced the array to 2 dimensions.
  
      Notice that the diagonalized out dimensions are added as new dimensions at
      the end of the indices.
      """
  
      def __new__(cls, expr, *diagonal_indices, **kwargs):
          expr = _sympify(expr)
          diagonal_indices = [Tuple(*sorted(i)) for i in diagonal_indices]
          normalize = kwargs.get("normalize", True)
  
          shape = get_shape(expr)
          if shape is not None:
              cls._validate(expr, *diagonal_indices)
              # Get new shape:
              positions, shape = cls._get_positions_shape(shape, diagonal_indices)
          else:
              positions = None
          if len(diagonal_indices) == 0:
              return expr
          obj = Basic.__new__(cls, expr, *diagonal_indices)
          obj._positions = positions
          obj._subranks = _get_subranks(expr)
          obj._shape = shape
          if normalize:
              return obj._normalize()
          return obj
  
      def _normalize(self):
          expr = self.expr
          diagonal_indices = self.diagonal_indices
          if isinstance(expr, ArrayAdd):
              return self._ArrayDiagonal_denest_ArrayAdd(expr, *diagonal_indices)
          if isinstance(expr, ArrayDiagonal):
              return self._ArrayDiagonal_denest_ArrayDiagonal(expr, *diagonal_indices)
          if isinstance(expr, PermuteDims):
              return self._ArrayDiagonal_denest_PermuteDims(expr, *diagonal_indices)
          if isinstance(expr, (ZeroArray, ZeroMatrix)):
              positions, shape = self._get_positions_shape(expr.shape, diagonal_indices)
              return ZeroArray(*shape)
          return self.func(expr, *diagonal_indices, normalize=False)
  
      @staticmethod
      def _validate(expr, *diagonal_indices):
          # Check that no diagonalization happens on indices with mismatched
          # dimensions:
          shape = get_shape(expr)
          for i in diagonal_indices:
              if any(j >= len(shape) for j in i):
                  raise ValueError("index is larger than expression shape")
              if len({shape[j] for j in i}) != 1:
                  raise ValueError("diagonalizing indices of different dimensions")
              if len(i) <= 1:
                  raise ValueError("need at least two axes to diagonalize")
  
      @staticmethod
      def _remove_trivial_dimensions(shape, *diagonal_indices):
          return [tuple(j for j in i) for i in diagonal_indices if shape[i[0]] != 1]
  
      @property
      def expr(self):
          return self.args[0]
  
      @property
      def diagonal_indices(self):
          return self.args[1:]
  
      @staticmethod
      def _flatten(expr, *outer_diagonal_indices):
          inner_diagonal_indices = expr.diagonal_indices
          all_inner = [j for i in inner_diagonal_indices for j in i]
          all_inner.sort()
          # TODO: add API for total rank and cumulative rank:
          total_rank = _get_subrank(expr)
          inner_rank = len(all_inner)
          outer_rank = total_rank - inner_rank
          shifts = [0 for i in range(outer_rank)]
          counter = 0
          pointer = 0
          for i in range(outer_rank):
              while pointer < inner_rank and counter >= all_inner[pointer]:
                  counter += 1
                  pointer += 1
              shifts[i] += pointer
              counter += 1
          outer_diagonal_indices = tuple(tuple(shifts[j] + j for j in i) for i in outer_diagonal_indices)
          diagonal_indices = inner_diagonal_indices + outer_diagonal_indices
          return ArrayDiagonal(expr.expr, *diagonal_indices)
  
      @classmethod
      def _ArrayDiagonal_denest_ArrayAdd(cls, expr, *diagonal_indices):
          return ArrayAdd(*[ArrayDiagonal(arg, *diagonal_indices) for arg in expr.args])
  
      @classmethod
      def _ArrayDiagonal_denest_ArrayDiagonal(cls, expr, *diagonal_indices):
          return cls._flatten(expr, *diagonal_indices)
  
      @classmethod
      def _ArrayDiagonal_denest_PermuteDims(cls, expr: PermuteDims, *diagonal_indices):
          back_diagonal_indices = [[expr.permutation(j) for j in i] for i in diagonal_indices]
          nondiag = [i for i in range(get_rank(expr)) if not any(i in j for j in diagonal_indices)]
          back_nondiag = [expr.permutation(i) for i in nondiag]
          remap = {e: i for i, e in enumerate(sorted(back_nondiag))}
          new_permutation1 = [remap[i] for i in back_nondiag]
          shift = len(new_permutation1)
          diag_block_perm = [i + shift for i in range(len(back_diagonal_indices))]
          new_permutation = new_permutation1 + diag_block_perm
          return PermuteDims(
              ArrayDiagonal(
                  expr.expr,
                  *back_diagonal_indices
              ),
              new_permutation
          )
  
      def _push_indices_down_nonstatic(self, indices):
          transform = lambda x: self._positions[x] if x < len(self._positions) else None
          return _apply_recursively_over_nested_lists(transform, indices)
  
      def _push_indices_up_nonstatic(self, indices):
  
          def transform(x):
              for i, e in enumerate(self._positions):
                  if (isinstance(e, int) and x == e) or (isinstance(e, tuple) and x in e):
                      return i
  
          return _apply_recursively_over_nested_lists(transform, indices)
  
      @classmethod
      def _push_indices_down(cls, diagonal_indices, indices, rank):
          positions, shape = cls._get_positions_shape(range(rank), diagonal_indices)
          transform = lambda x: positions[x] if x < len(positions) else None
          return _apply_recursively_over_nested_lists(transform, indices)
  
      @classmethod
      def _push_indices_up(cls, diagonal_indices, indices, rank):
          positions, shape = cls._get_positions_shape(range(rank), diagonal_indices)
  
          def transform(x):
              for i, e in enumerate(positions):
                  if (isinstance(e, int) and x == e) or (isinstance(e, (tuple, Tuple)) and (x in e)):
                      return i
  
          return _apply_recursively_over_nested_lists(transform, indices)
  
      @classmethod
      def _get_positions_shape(cls, shape, diagonal_indices):
          data1 = tuple((i, shp) for i, shp in enumerate(shape) if not any(i in j for j in diagonal_indices))
          pos1, shp1 = zip(*data1) if data1 else ((), ())
          data2 = tuple((i, shape[i[0]]) for i in diagonal_indices)
          pos2, shp2 = zip(*data2) if data2 else ((), ())
          positions = pos1 + pos2
          shape = shp1 + shp2
          return positions, shape
  
      def as_explicit(self):
          return tensordiagonal(self.expr.as_explicit(), *self.diagonal_indices)
  
  
  class ArrayElementwiseApplyFunc(_CodegenArrayAbstract):
  
      def __new__(cls, function, element):
  
          if not isinstance(function, Lambda):
              d = Dummy('d')
              function = Lambda(d, function(d))
  
          obj = _CodegenArrayAbstract.__new__(cls, function, element)
          obj._subranks = _get_subranks(element)
          return obj
  
      @property
      def function(self):
          return self.args[0]
  
      @property
      def expr(self):
          return self.args[1]
  
      @property
      def shape(self):
          return self.expr.shape
  
      def _get_function_fdiff(self):
          d = Dummy("d")
          function = self.function(d)
          fdiff = function.diff(d)
          if isinstance(fdiff, Function):
              fdiff = type(fdiff)
          else:
              fdiff = Lambda(d, fdiff)
          return fdiff
  
  
  class ArrayContraction(_CodegenArrayAbstract):
      r"""
      This class is meant to represent contractions of arrays in a form easily
      processable by the code printers.
      """
  
      def __new__(cls, expr, *contraction_indices, **kwargs):
          contraction_indices = _sort_contraction_indices(contraction_indices)
          expr = _sympify(expr)
  
          normalize = kwargs.get("normalize", True)
  
          obj = Basic.__new__(cls, expr, *contraction_indices)
          obj._subranks = _get_subranks(expr)
          obj._mapping = _get_mapping_from_subranks(obj._subranks)
  
          free_indices_to_position = {i: i for i in range(sum(obj._subranks)) if all(i not in cind for cind in contraction_indices)}
          obj._free_indices_to_position = free_indices_to_position
  
          shape = get_shape(expr)
          cls._validate(expr, *contraction_indices)
          if shape:
              shape = tuple(shp for i, shp in enumerate(shape) if not any(i in j for j in contraction_indices))
          obj._shape = shape
          if normalize:
              return obj._normalize()
          return obj
  
      def _normalize(self):
          expr = self.expr
          contraction_indices = self.contraction_indices
  
          if len(contraction_indices) == 0:
              return expr
  
          if isinstance(expr, ArrayContraction):
              return self._ArrayContraction_denest_ArrayContraction(expr, *contraction_indices)
  
          if isinstance(expr, (ZeroArray, ZeroMatrix)):
              return self._ArrayContraction_denest_ZeroArray(expr, *contraction_indices)
  
          if isinstance(expr, PermuteDims):
              return self._ArrayContraction_denest_PermuteDims(expr, *contraction_indices)
  
          if isinstance(expr, ArrayTensorProduct):
              expr, contraction_indices = self._sort_fully_contracted_args(expr, contraction_indices)
              expr, contraction_indices = self._lower_contraction_to_addends(expr, contraction_indices)
              if len(contraction_indices) == 0:
                  return expr
  
          if isinstance(expr, ArrayDiagonal):
              return self._ArrayContraction_denest_ArrayDiagonal(expr, *contraction_indices)
  
          if isinstance(expr, ArrayAdd):
              return self._ArrayContraction_denest_ArrayAdd(expr, *contraction_indices)
  
          return self.func(expr, *contraction_indices, normalize=False)
  
      def __mul__(self, other):
          if other == 1:
              return self
          else:
              raise NotImplementedError("Product of N-dim arrays is not uniquely defined. Use another method.")
  
      def __rmul__(self, other):
          if other == 1:
              return self
          else:
              raise NotImplementedError("Product of N-dim arrays is not uniquely defined. Use another method.")
  
      @staticmethod
      def _validate(expr, *contraction_indices):
          shape = get_shape(expr)
          if shape is None:
              return
  
          # Check that no contraction happens when the shape is mismatched:
          for i in contraction_indices:
              if len({shape[j] for j in i if shape[j] != -1}) != 1:
                  raise ValueError("contracting indices of different dimensions")
  
      @classmethod
      def _push_indices_down(cls, contraction_indices, indices):
          flattened_contraction_indices = [j for i in contraction_indices for j in i]
          flattened_contraction_indices.sort()
          transform = _build_push_indices_down_func_transformation(flattened_contraction_indices)
          return _apply_recursively_over_nested_lists(transform, indices)
  
      @classmethod
      def _push_indices_up(cls, contraction_indices, indices):
          flattened_contraction_indices = [j for i in contraction_indices for j in i]
          flattened_contraction_indices.sort()
          transform = _build_push_indices_up_func_transformation(flattened_contraction_indices)
          return _apply_recursively_over_nested_lists(transform, indices)
  
      @classmethod
      def _lower_contraction_to_addends(cls, expr, contraction_indices):
          if isinstance(expr, ArrayAdd):
              raise NotImplementedError()
          if not isinstance(expr, ArrayTensorProduct):
              return expr, contraction_indices
          subranks = expr.subranks
          cumranks = list(accumulate([0] + subranks))
          contraction_indices_remaining = []
          contraction_indices_args = [[] for i in expr.args]
          backshift = set([])
          for i, contraction_group in enumerate(contraction_indices):
              for j in range(len(expr.args)):
                  if not isinstance(expr.args[j], ArrayAdd):
                      continue
                  if all(cumranks[j] <= k < cumranks[j+1] for k in contraction_group):
                      contraction_indices_args[j].append([k - cumranks[j] for k in contraction_group])
                      backshift.update(contraction_group)
                      break
              else:
                  contraction_indices_remaining.append(contraction_group)
          if len(contraction_indices_remaining) == len(contraction_indices):
              return expr, contraction_indices
          total_rank = get_rank(expr)
          shifts = list(accumulate([1 if i in backshift else 0 for i in range(total_rank)]))
          contraction_indices_remaining = [Tuple.fromiter(j - shifts[j] for j in i) for i in contraction_indices_remaining]
          ret = ArrayTensorProduct(*[
              ArrayContraction(arg, *contr) for arg, contr in zip(expr.args, contraction_indices_args)
          ])
          return ret, contraction_indices_remaining
  
      def split_multiple_contractions(self):
          """
          Recognize multiple contractions and attempt at rewriting them as paired-contractions.
  
          This allows some contractions involving more than two indices to be
          rewritten as multiple contractions involving two indices, thus allowing
          the expression to be rewritten as a matrix multiplication line.
  
          Examples:
  
          * `A_ij b_j0 C_jk` ===> `A*DiagMatrix(b)*C`
  
          Care for:
          - matrix being diagonalized (i.e. `A_ii`)
          - vectors being diagonalized (i.e. `a_i0`)
  
          Multiple contractions can be split into matrix multiplications if
          not more than two arguments are non-diagonals or non-vectors.
          Vectors get diagonalized while diagonal matrices remain diagonal.
          The non-diagonal matrices can be at the beginning or at the end
          of the final matrix multiplication line.
          """
  
          editor = _EditArrayContraction(self)
  
          contraction_indices = self.contraction_indices
  
          onearray_insert = []
  
          for indl, links in enumerate(contraction_indices):
              if len(links) <= 2:
                  continue
  
              # Check multiple contractions:
              #
              # Examples:
              #
              # * `A_ij b_j0 C_jk` ===> `A*DiagMatrix(b)*C \otimes OneArray(1)` with permutation (1 2)
              #
              # Care for:
              # - matrix being diagonalized (i.e. `A_ii`)
              # - vectors being diagonalized (i.e. `a_i0`)
  
              # Multiple contractions can be split into matrix multiplications if
              # not more than three arguments are non-diagonals or non-vectors.
              #
              # Vectors get diagonalized while diagonal matrices remain diagonal.
              # The non-diagonal matrices can be at the beginning or at the end
              # of the final matrix multiplication line.
  
              positions = editor.get_mapping_for_index(indl)
  
              # Also consider the case of diagonal matrices being contracted:
              current_dimension = self.expr.shape[links[0]]
  
              not_vectors: tTuple[_ArgE, int] = []
              vectors: tTuple[_ArgE, int] = []
              for arg_ind, rel_ind in positions:
                  arg = editor.args_with_ind[arg_ind]
                  mat = arg.element
                  abs_arg_start, abs_arg_end = editor.get_absolute_range(arg)
                  other_arg_pos = 1-rel_ind
                  other_arg_abs = abs_arg_start + other_arg_pos
                  if ((1 not in mat.shape) or
                      ((current_dimension == 1) is True and mat.shape != (1, 1)) or
                      any(other_arg_abs in l for li, l in enumerate(contraction_indices) if li != indl)
                  ):
                      not_vectors.append((arg, rel_ind))
                  else:
                      vectors.append((arg, rel_ind))
              if len(not_vectors) > 2:
                  # If more than two arguments in the multiple contraction are
                  # non-vectors and non-diagonal matrices, we cannot find a way
                  # to split this contraction into a matrix multiplication line:
                  continue
              # Three cases to handle:
              # - zero non-vectors
              # - one non-vector
              # - two non-vectors
              for v, rel_ind in vectors:
                  v.element = diagonalize_vector(v.element)
              vectors_to_loop = not_vectors[:1] + vectors + not_vectors[1:]
              first_not_vector, rel_ind = vectors_to_loop[0]
              new_index = first_not_vector.indices[rel_ind]
  
              for v, rel_ind in vectors_to_loop[1:-1]:
                  v.indices[rel_ind] = new_index
                  new_index = editor.get_new_contraction_index()
                  assert v.indices.index(None) == 1 - rel_ind
                  v.indices[v.indices.index(None)] = new_index
                  onearray_insert.append(v)
  
              last_vec, rel_ind = vectors_to_loop[-1]
              last_vec.indices[rel_ind] = new_index
  
          for v in onearray_insert:
              editor.insert_after(v, _ArgE(OneArray(1), [None]))
  
          return editor.to_array_contraction()
  
      def flatten_contraction_of_diagonal(self):
          if not isinstance(self.expr, ArrayDiagonal):
              return self
          contraction_down = self.expr._push_indices_down(self.expr.diagonal_indices, self.contraction_indices)
          new_contraction_indices = []
          diagonal_indices = self.expr.diagonal_indices[:]
          for i in contraction_down:
              contraction_group = list(i)
              for j in i:
                  diagonal_with = [k for k in diagonal_indices if j in k]
                  contraction_group.extend([l for k in diagonal_with for l in k])
                  diagonal_indices = [k for k in diagonal_indices if k not in diagonal_with]
              new_contraction_indices.append(sorted(set(contraction_group)))
  
          new_contraction_indices = ArrayDiagonal._push_indices_up(diagonal_indices, new_contraction_indices)
          return ArrayContraction(
              ArrayDiagonal(
                  self.expr.expr,
                  *diagonal_indices
              ),
              *new_contraction_indices
          )
  
      @staticmethod
      def _get_free_indices_to_position_map(free_indices, contraction_indices):
          free_indices_to_position = {}
          flattened_contraction_indices = [j for i in contraction_indices for j in i]
          counter = 0
          for ind in free_indices:
              while counter in flattened_contraction_indices:
                  counter += 1
              free_indices_to_position[ind] = counter
              counter += 1
          return free_indices_to_position
  
      @staticmethod
      def _get_index_shifts(expr):
          """
          Get the mapping of indices at the positions before the contraction
          occurs.
  
          Examples
          ========
  
          >>> from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
          >>> from sympy.tensor.array.expressions.array_expressions import ArrayContraction
          >>> from sympy import MatrixSymbol
          >>> M = MatrixSymbol("M", 3, 3)
          >>> N = MatrixSymbol("N", 3, 3)
          >>> cg = ArrayContraction(ArrayTensorProduct(M, N), [1, 2])
          >>> cg._get_index_shifts(cg)
          [0, 2]
  
          Indeed, ``cg`` after the contraction has two dimensions, 0 and 1. They
          need to be shifted by 0 and 2 to get the corresponding positions before
          the contraction (that is, 0 and 3).
          """
          inner_contraction_indices = expr.contraction_indices
          all_inner = [j for i in inner_contraction_indices for j in i]
          all_inner.sort()
          # TODO: add API for total rank and cumulative rank:
          total_rank = _get_subrank(expr)
          inner_rank = len(all_inner)
          outer_rank = total_rank - inner_rank
          shifts = [0 for i in range(outer_rank)]
          counter = 0
          pointer = 0
          for i in range(outer_rank):
              while pointer < inner_rank and counter >= all_inner[pointer]:
                  counter += 1
                  pointer += 1
              shifts[i] += pointer
              counter += 1
          return shifts
  
      @staticmethod
      def _convert_outer_indices_to_inner_indices(expr, *outer_contraction_indices):
          shifts = ArrayContraction._get_index_shifts(expr)
          outer_contraction_indices = tuple(tuple(shifts[j] + j for j in i) for i in outer_contraction_indices)
          return outer_contraction_indices
  
      @staticmethod
      def _flatten(expr, *outer_contraction_indices):
          inner_contraction_indices = expr.contraction_indices
          outer_contraction_indices = ArrayContraction._convert_outer_indices_to_inner_indices(expr, *outer_contraction_indices)
          contraction_indices = inner_contraction_indices + outer_contraction_indices
          return ArrayContraction(expr.expr, *contraction_indices)
  
      @classmethod
      def _ArrayContraction_denest_ArrayContraction(cls, expr, *contraction_indices):
          return cls._flatten(expr, *contraction_indices)
  
      @classmethod
      def _ArrayContraction_denest_ZeroArray(cls, expr, *contraction_indices):
          contraction_indices_flat = [j for i in contraction_indices for j in i]
          shape = [e for i, e in enumerate(expr.shape) if i not in contraction_indices_flat]
          return ZeroArray(*shape)
  
      @classmethod
      def _ArrayContraction_denest_ArrayAdd(cls, expr, *contraction_indices):
          return ArrayAdd(*[ArrayContraction(i, *contraction_indices) for i in expr.args])
  
      @classmethod
      def _ArrayContraction_denest_PermuteDims(cls, expr, *contraction_indices):
          permutation = expr.permutation
          plist = permutation.array_form
          new_contraction_indices = [tuple(permutation(j) for j in i) for i in contraction_indices]
          new_plist = [i for i in plist if not any(i in j for j in new_contraction_indices)]
          new_plist = cls._push_indices_up(new_contraction_indices, new_plist)
          return PermuteDims(
              ArrayContraction(expr.expr, *new_contraction_indices),
              Permutation(new_plist)
          )
  
      @classmethod
      def _ArrayContraction_denest_ArrayDiagonal(cls, expr: 'ArrayDiagonal', *contraction_indices):
          diagonal_indices = list(expr.diagonal_indices)
          down_contraction_indices = expr._push_indices_down(expr.diagonal_indices, contraction_indices, get_rank(expr.expr))
          # Flatten diagonally contracted indices:
          down_contraction_indices = [[k for j in i for k in (j if isinstance(j, (tuple, Tuple)) else [j])] for i in down_contraction_indices]
          new_contraction_indices = []
          for contr_indgrp in down_contraction_indices:
              ind = contr_indgrp[:]
              for j, diag_indgrp in enumerate(diagonal_indices):
                  if diag_indgrp is None:
                      continue
                  if any(i in diag_indgrp for i in contr_indgrp):
                      ind.extend(diag_indgrp)
                      diagonal_indices[j] = None
              new_contraction_indices.append(sorted(set(ind)))
  
          new_diagonal_indices_down = [i for i in diagonal_indices if i is not None]
          new_diagonal_indices = ArrayContraction._push_indices_up(new_contraction_indices, new_diagonal_indices_down)
          return ArrayDiagonal(
              ArrayContraction(expr.expr, *new_contraction_indices),
              *new_diagonal_indices
          )
  
      @classmethod
      def _sort_fully_contracted_args(cls, expr, contraction_indices):
          if expr.shape is None:
              return expr, contraction_indices
          cumul = list(accumulate([0] + expr.subranks))
          index_blocks = [list(range(cumul[i], cumul[i+1])) for i in range(len(expr.args))]
          contraction_indices_flat = {j for i in contraction_indices for j in i}
          fully_contracted = [all(j in contraction_indices_flat for j in range(cumul[i], cumul[i+1])) for i, arg in enumerate(expr.args)]
          new_pos = sorted(range(len(expr.args)), key=lambda x: (0, default_sort_key(expr.args[x])) if fully_contracted[x] else (1,))
          new_args = [expr.args[i] for i in new_pos]
          new_index_blocks_flat = [j for i in new_pos for j in index_blocks[i]]
          index_permutation_array_form = _af_invert(new_index_blocks_flat)
          new_contraction_indices = [tuple(index_permutation_array_form[j] for j in i) for i in contraction_indices]
          new_contraction_indices = _sort_contraction_indices(new_contraction_indices)
          return ArrayTensorProduct(*new_args), new_contraction_indices
  
      def _get_contraction_tuples(self):
          r"""
          Return tuples containing the argument index and position within the
          argument of the index position.
  
          Examples
          ========
  
          >>> from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
          >>> from sympy import MatrixSymbol
          >>> from sympy.abc import N
          >>> from sympy.tensor.array.expressions.array_expressions import ArrayContraction
          >>> A = MatrixSymbol("A", N, N)
          >>> B = MatrixSymbol("B", N, N)
  
          >>> cg = ArrayContraction(ArrayTensorProduct(A, B), (1, 2))
          >>> cg._get_contraction_tuples()
          [[(0, 1), (1, 0)]]
  
          Notes
          =====
  
          Here the contraction pair `(1, 2)` meaning that the 2nd and 3rd indices
          of the tensor product `A\otimes B` are contracted, has been transformed
          into `(0, 1)` and `(1, 0)`, identifying the same indices in a different
          notation. `(0, 1)` is the second index (1) of the first argument (i.e.
                  0 or `A`). `(1, 0)` is the first index (i.e. 0) of the second
          argument (i.e. 1 or `B`).
          """
          mapping = self._mapping
          return [[mapping[j] for j in i] for i in self.contraction_indices]
  
      @staticmethod
      def _contraction_tuples_to_contraction_indices(expr, contraction_tuples):
          # TODO: check that `expr` has `.subranks`:
          ranks = expr.subranks
          cumulative_ranks = [0] + list(accumulate(ranks))
          return [tuple(cumulative_ranks[j]+k for j, k in i) for i in contraction_tuples]
  
      @property
      def free_indices(self):
          return self._free_indices[:]
  
      @property
      def free_indices_to_position(self):
          return dict(self._free_indices_to_position)
  
      @property
      def expr(self):
          return self.args[0]
  
      @property
      def contraction_indices(self):
          return self.args[1:]
  
      def _contraction_indices_to_components(self):
          expr = self.expr
          if not isinstance(expr, ArrayTensorProduct):
              raise NotImplementedError("only for contractions of tensor products")
          ranks = expr.subranks
          mapping = {}
          counter = 0
          for i, rank in enumerate(ranks):
              for j in range(rank):
                  mapping[counter] = (i, j)
                  counter += 1
          return mapping
  
      def sort_args_by_name(self):
          """
          Sort arguments in the tensor product so that their order is lexicographical.
  
          Examples
          ========
  
          >>> from sympy.tensor.array.expressions.conv_matrix_to_array import convert_matrix_to_array
          >>> from sympy import MatrixSymbol
          >>> from sympy.abc import N
          >>> A = MatrixSymbol("A", N, N)
          >>> B = MatrixSymbol("B", N, N)
          >>> C = MatrixSymbol("C", N, N)
          >>> D = MatrixSymbol("D", N, N)
  
          >>> cg = convert_matrix_to_array(C*D*A*B)
          >>> cg
          ArrayContraction(ArrayTensorProduct(A, D, C, B), (0, 3), (1, 6), (2, 5))
          >>> cg.sort_args_by_name()
          ArrayContraction(ArrayTensorProduct(A, D, B, C), (0, 3), (1, 4), (2, 7))
          """
          expr = self.expr
          if not isinstance(expr, ArrayTensorProduct):
              return self
          args = expr.args
          sorted_data = sorted(enumerate(args), key=lambda x: default_sort_key(x[1]))
          pos_sorted, args_sorted = zip(*sorted_data)
          reordering_map = {i: pos_sorted.index(i) for i, arg in enumerate(args)}
          contraction_tuples = self._get_contraction_tuples()
          contraction_tuples = [[(reordering_map[j], k) for j, k in i] for i in contraction_tuples]
          c_tp = ArrayTensorProduct(*args_sorted)
          new_contr_indices = self._contraction_tuples_to_contraction_indices(
                  c_tp,
                  contraction_tuples
          )
          return ArrayContraction(c_tp, *new_contr_indices)
  
      def _get_contraction_links(self):
          r"""
          Returns a dictionary of links between arguments in the tensor product
          being contracted.
  
          See the example for an explanation of the values.
  
          Examples
          ========
  
          >>> from sympy import MatrixSymbol
          >>> from sympy.abc import N
          >>> from sympy.tensor.array.expressions.conv_matrix_to_array import convert_matrix_to_array
          >>> A = MatrixSymbol("A", N, N)
          >>> B = MatrixSymbol("B", N, N)
          >>> C = MatrixSymbol("C", N, N)
          >>> D = MatrixSymbol("D", N, N)
  
          Matrix multiplications are pairwise contractions between neighboring
          matrices:
  
          `A_{ij} B_{jk} C_{kl} D_{lm}`
  
          >>> cg = convert_matrix_to_array(A*B*C*D)
          >>> cg
          ArrayContraction(ArrayTensorProduct(B, C, A, D), (0, 5), (1, 2), (3, 6))
  
          >>> cg._get_contraction_links()
          {0: {0: (2, 1), 1: (1, 0)}, 1: {0: (0, 1), 1: (3, 0)}, 2: {1: (0, 0)}, 3: {0: (1, 1)}}
  
          This dictionary is interpreted as follows: argument in position 0 (i.e.
          matrix `A`) has its second index (i.e. 1) contracted to `(1, 0)`, that
          is argument in position 1 (matrix `B`) on the first index slot of `B`,
          this is the contraction provided by the index `j` from `A`.
  
          The argument in position 1 (that is, matrix `B`) has two contractions,
          the ones provided by the indices `j` and `k`, respectively the first
          and second indices (0 and 1 in the sub-dict).  The link `(0, 1)` and
          `(2, 0)` respectively. `(0, 1)` is the index slot 1 (the 2nd) of
          argument in position 0 (that is, `A_{\ldot j}`), and so on.
          """
          args, dlinks = _get_contraction_links([self], self.subranks, *self.contraction_indices)
          return dlinks
  
      def as_explicit(self):
          return tensorcontraction(self.expr.as_explicit(), *self.contraction_indices)
  
  
  class _ArgE:
      """
      The ``_ArgE`` object contains references to the array expression
      (``.element``) and a list containing the information about index
      contractions (``.indices``).
  
      Index contractions are numbered and contracted indices show the number of
      the contraction. Uncontracted indices have ``None`` value.
  
      For example:
      ``_ArgE(M, [None, 3])``
      This object means that expression ``M`` is part of an array contraction
      and has two indices, the first is not contracted (value ``None``),
      the second index is contracted to the 4th (i.e. number ``3``) group of the
      array contraction object.
      """
      indices: List[Optional[int]]
  
      def __init__(self, element, indices: Optional[List[Optional[int]]] = None):
          self.element = element
          if indices is None:
              self.indices = [None for i in range(get_rank(element))]
          else:
              self.indices = indices
  
      def __str__(self):
          return "_ArgE(%s, %s)" % (self.element, self.indices)
  
      __repr__ = __str__
  
  
  class _IndPos:
      """
      Index position, requiring two integers in the constructor:
  
      - arg: the position of the argument in the tensor product,
      - rel: the relative position of the index inside the argument.
      """
      def __init__(self, arg: int, rel: int):
          self.arg = arg
          self.rel = rel
  
      def __str__(self):
          return "_IndPos(%i, %i)" % (self.arg, self.rel)
  
      __repr__ = __str__
  
      def __iter__(self):
          yield from [self.arg, self.rel]
  
  
  class _EditArrayContraction:
      """
      Utility class to help manipulate array contraction objects.
  
      This class takes as input an ``ArrayContraction`` object and turns it into
      an editable object.
  
      The field ``args_with_ind`` of this class is a list of ``_ArgE`` objects
      which can be used to easily edit the contraction structure of the
      expression.
  
      Once editing is finished, the ``ArrayContraction`` object may be recreated
      by calling the ``.to_array_contraction()`` method.
      """
  
      def __init__(self, base_array: typing.Union[ArrayContraction, ArrayDiagonal, ArrayTensorProduct]):
  
          expr: Basic
          diagonalized: tTuple[tTuple[int, ...], ...]
          contraction_indices: List[tTuple[int]]
          if isinstance(base_array, ArrayContraction):
              mapping = _get_mapping_from_subranks(base_array.subranks)
              expr = base_array.expr
              contraction_indices = base_array.contraction_indices
              diagonalized = ()
          elif isinstance(base_array, ArrayDiagonal):
  
              if isinstance(base_array.expr, ArrayContraction):
                  mapping = _get_mapping_from_subranks(base_array.expr.subranks)
                  expr = base_array.expr.expr
                  diagonalized = ArrayContraction._push_indices_down(base_array.expr.contraction_indices, base_array.diagonal_indices)
                  contraction_indices = base_array.expr.contraction_indices
              elif isinstance(base_array.expr, ArrayTensorProduct):
                  mapping = {}
                  expr = base_array.expr
                  diagonalized = base_array.diagonal_indices
                  contraction_indices = []
              else:
                  mapping = {}
                  expr = base_array.expr
                  diagonalized = base_array.diagonal_indices
                  contraction_indices = []
  
          elif isinstance(base_array, ArrayTensorProduct):
              expr = base_array
              contraction_indices = []
              diagonalized = ()
          else:
              raise NotImplementedError()
  
          if isinstance(expr, ArrayTensorProduct):
              args = list(expr.args)
          else:
              args = [expr]
  
          args_with_ind: List[_ArgE] = [_ArgE(arg) for arg in args]
          for i, contraction_tuple in enumerate(contraction_indices):
              for j in contraction_tuple:
                  arg_pos, rel_pos = mapping[j]
                  args_with_ind[arg_pos].indices[rel_pos] = i
          self.args_with_ind: List[_ArgE] = args_with_ind
          self.number_of_contraction_indices: int = len(contraction_indices)
          self._track_permutation: Optional[List[List[int]]] = None
  
          mapping = _get_mapping_from_subranks(base_array.subranks)
  
          # Trick: add diagonalized indices as negative indices into the editor object:
          for i, e in enumerate(diagonalized):
              for j in e:
                  arg_pos, rel_pos = mapping[j]
                  self.args_with_ind[arg_pos].indices[rel_pos] = -1 - i
  
      def insert_after(self, arg: _ArgE, new_arg: _ArgE):
          pos = self.args_with_ind.index(arg)
          self.args_with_ind.insert(pos + 1, new_arg)
  
      def get_new_contraction_index(self):
          self.number_of_contraction_indices += 1
          return self.number_of_contraction_indices - 1
  
      def refresh_indices(self):
          updates: tDict[int, int] = {}
          for arg_with_ind in self.args_with_ind:
              updates.update({i: -1 for i in arg_with_ind.indices if i is not None})
          for i, e in enumerate(sorted(updates)):
              updates[e] = i
          self.number_of_contraction_indices: int = len(updates)
          for arg_with_ind in self.args_with_ind:
              arg_with_ind.indices = [updates.get(i, None) for i in arg_with_ind.indices]
  
      def merge_scalars(self):
          scalars = []
          for arg_with_ind in self.args_with_ind:
              if len(arg_with_ind.indices) == 0:
                  scalars.append(arg_with_ind)
          for i in scalars:
              self.args_with_ind.remove(i)
          scalar = Mul.fromiter([i.element for i in scalars])
          if len(self.args_with_ind) == 0:
              self.args_with_ind.append(_ArgE(scalar))
          else:
              from sympy.tensor.array.expressions.conv_array_to_matrix import _a2m_tensor_product
              self.args_with_ind[0].element = _a2m_tensor_product(scalar, self.args_with_ind[0].element)
  
      def to_array_contraction(self):
  
          # Count the ranks of the arguments:
          counter = 0
          # Create a collector for the new diagonal indices:
          diag_indices = defaultdict(list)
  
          count_index_freq = Counter()
          for arg_with_ind in self.args_with_ind:
              count_index_freq.update(Counter(arg_with_ind.indices))
  
          free_index_count = count_index_freq[None]
  
          # Construct the inverse permutation:
          inv_perm1 = []
          inv_perm2 = []
          # Keep track of which diagonal indices have already been processed:
          done = set([])
  
          # Counter for the diagonal indices:
          counter4 = 0
  
          for arg_with_ind in self.args_with_ind:
              # If some diagonalization axes have been removed, they should be
              # permuted in order to keep the permutation.
              # Add permutation here
              counter2 = 0  # counter for the indices
              for i in arg_with_ind.indices:
                  if i is None:
                      inv_perm1.append(counter4)
                      counter2 += 1
                      counter4 += 1
                      continue
                  if i >= 0:
                      continue
                  # Reconstruct the diagonal indices:
                  diag_indices[-1 - i].append(counter + counter2)
                  if count_index_freq[i] == 1 and i not in done:
                      inv_perm1.append(free_index_count - 1 - i)
                      done.add(i)
                  elif i not in done:
                      inv_perm2.append(free_index_count - 1 - i)
                      done.add(i)
                  counter2 += 1
              # Remove negative indices to restore a proper editor object:
              arg_with_ind.indices = [i if i is not None and i >= 0 else None for i in arg_with_ind.indices]
              counter += len([i for i in arg_with_ind.indices if i is None or i < 0])
  
          inverse_permutation = inv_perm1 + inv_perm2
          permutation = _af_invert(inverse_permutation)
  
          # Get the diagonal indices after the detection of HadamardProduct in the expression:
          diag_indices_filtered = [tuple(v) for v in diag_indices.values() if len(v) > 1]
  
          self.merge_scalars()
          self.refresh_indices()
          args = [arg.element for arg in self.args_with_ind]
          contraction_indices = self.get_contraction_indices()
          expr = ArrayContraction(ArrayTensorProduct(*args), *contraction_indices)
          expr2 = ArrayDiagonal(expr, *diag_indices_filtered)
          if self._track_permutation is not None:
              permutation2 = _af_invert([j for i in self._track_permutation for j in i])
              expr2 = PermuteDims(expr2, permutation2)
  
          expr3 = PermuteDims(expr2, permutation)
          return expr3
  
      def get_contraction_indices(self) -> List[List[int]]:
          contraction_indices: List[List[int]] = [[] for i in range(self.number_of_contraction_indices)]
          current_position: int = 0
          for i, arg_with_ind in enumerate(self.args_with_ind):
              for j in arg_with_ind.indices:
                  if j is not None:
                      contraction_indices[j].append(current_position)
                  current_position += 1
          return contraction_indices
  
      def get_mapping_for_index(self, ind) -> List[_IndPos]:
          if ind >= self.number_of_contraction_indices:
              raise ValueError("index value exceeding the index range")
          positions: List[_IndPos] = []
          for i, arg_with_ind in enumerate(self.args_with_ind):
              for j, arg_ind in enumerate(arg_with_ind.indices):
                  if ind == arg_ind:
                      positions.append(_IndPos(i, j))
          return positions
  
      def get_contraction_indices_to_ind_rel_pos(self) -> List[List[_IndPos]]:
          contraction_indices: List[List[_IndPos]] = [[] for i in range(self.number_of_contraction_indices)]
          for i, arg_with_ind in enumerate(self.args_with_ind):
              for j, ind in enumerate(arg_with_ind.indices):
                  if ind is not None:
                      contraction_indices[ind].append(_IndPos(i, j))
          return contraction_indices
  
      def count_args_with_index(self, index: int) -> int:
          """
          Count the number of arguments that have the given index.
          """
          counter: int = 0
          for arg_with_ind in self.args_with_ind:
              if index in arg_with_ind.indices:
                  counter += 1
          return counter
  
      def get_args_with_index(self, index: int) -> List[_ArgE]:
          """
          Get a list of arguments having the given index.
          """
          ret: List[_ArgE] = [i for i in self.args_with_ind if index in i.indices]
          return ret
  
      @property
      def number_of_diagonal_indices(self):
          data = set([])
          for arg in self.args_with_ind:
              data.update({i for i in arg.indices if i is not None and i < 0})
          return len(data)
  
      def track_permutation_start(self):
          permutation = []
          perm_diag = []
          counter: int = 0
          counter2: int = -1
          for arg_with_ind in self.args_with_ind:
              perm = []
              for i in arg_with_ind.indices:
                  if i is not None:
                      if i < 0:
                          perm_diag.append(counter2)
                          counter2 -= 1
                      continue
                  perm.append(counter)
                  counter += 1
              permutation.append(perm)
          max_ind = max([max(i) if i else -1 for i in permutation]) if permutation else -1
          perm_diag = [max_ind - i for i in perm_diag]
          self._track_permutation = permutation + [perm_diag]
  
      def track_permutation_merge(self, destination: _ArgE, from_element: _ArgE):
          index_destination = self.args_with_ind.index(destination)
          index_element = self.args_with_ind.index(from_element)
          self._track_permutation[index_destination].extend(self._track_permutation[index_element]) # type: ignore
          self._track_permutation.pop(index_element) # type: ignore
  
      def get_absolute_free_range(self, arg: _ArgE) -> typing.Tuple[int, int]:
          """
          Return the range of the free indices of the arg as absolute positions
          among all free indices.
          """
          counter = 0
          for arg_with_ind in self.args_with_ind:
              number_free_indices = len([i for i in arg_with_ind.indices if i is None])
              if arg_with_ind == arg:
                  return counter, counter + number_free_indices
              counter += number_free_indices
          raise IndexError("argument not found")
  
      def get_absolute_range(self, arg: _ArgE) -> typing.Tuple[int, int]:
          """
          Return the absolute range of indices for arg, disregarding dummy
          indices.
          """
          counter = 0
          for arg_with_ind in self.args_with_ind:
              number_indices = len(arg_with_ind.indices)
              if arg_with_ind == arg:
                  return counter, counter + number_indices
              counter += number_indices
          raise IndexError("argument not found")
  
  
  def get_rank(expr):
      if isinstance(expr, (MatrixExpr, MatrixElement)):
          return 2
      if isinstance(expr, _CodegenArrayAbstract):
          return len(expr.shape)
      if isinstance(expr, NDimArray):
          return expr.rank()
      if isinstance(expr, Indexed):
          return expr.rank
      if isinstance(expr, IndexedBase):
          shape = expr.shape
          if shape is None:
              return -1
          else:
              return len(shape)
      if hasattr(expr, "shape"):
          return len(expr.shape)
      return 0
  
  
  def _get_subrank(expr):
      if isinstance(expr, _CodegenArrayAbstract):
          return expr.subrank()
      return get_rank(expr)
  
  
  def _get_subranks(expr):
      if isinstance(expr, _CodegenArrayAbstract):
          return expr.subranks
      else:
          return [get_rank(expr)]
  
  
  def get_shape(expr):
      if hasattr(expr, "shape"):
          return expr.shape
      return ()
  
  
  def nest_permutation(expr):
      if isinstance(expr, PermuteDims):
          return expr.nest_permutation()
      else:
          return expr
  