Repository: sympy/sympy

Issue:
NumPyPrinter incorrectly prints ArrayTensorProduct
`NumPyPrinter` assumes that all arguments of `ArrayTensorProduct` are 2 dimensional arrays (matrices). `TensorflowPrinter` correctly implements this for arbitrary shapes (including scalars):
```python
>>> from sympy.printing.numpy import NumPyPrinter
>>> from sympy.printing.tensorflow import TensorflowPrinter
>>> from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct, ArraySymbol
>>> print(NumPyPrinter().doprint(ArrayTensorProduct(-1, ArraySymbol('A',(3,3,3)))))
numpy.einsum(-1, [0, 1], A, [2, 3])
>>> print(TensorflowPrinter().doprint(ArrayTensorProduct(-1, ArraySymbol('A',(3,3,3)))))
tensorflow.linalg.einsum(",abc", -1, A)
```

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 416399b4c2..ef894343dd 100644
--- a/.mailmap
+++ b/.mailmap
@@ -915,9 +915,7 @@ Raymond Wong <rayman_407@yahoo.com> <rayman@CoolGuy.(none)>
 Raymond Wong <rayman_407@yahoo.com> <rayman@CoolGuyJr.(none)>
 Raymond Wong <rayman_407@yahoo.com> <rayman@CoolGuyJr>
 Rehas Sachdeva <aquannie@gmail.com>
-Remco de Boer <29308176+redeboer@users.noreply.github.com>
 Remco de Boer <redeboer@gmx.com>
-Remco de Boer <redeboer@outlook.com>
 Renato Coutinho <renato.coutinho@gmail.com>
 Renato Orsino <renato.orsino@gmail.com>
 Rhea Parekh <rheaparekh12@gmail.com>
diff --git a/doc/src/modules/stats.rst b/doc/src/modules/stats.rst
index dd2c2d1d9d..dca5c3d038 100644
--- a/doc/src/modules/stats.rst
+++ b/doc/src/modules/stats.rst
@@ -146,7 +146,6 @@ Interface
 .. autofunction:: coskewness
 .. autofunction:: median
 .. autofunction:: std
-.. autofunction:: quantile
 .. autofunction:: sample
 .. autofunction:: sample_iter
 .. autofunction:: factorial_moment
diff --git a/setup.py b/setup.py
index d6cbea4205..1bc82b26a2 100755
--- a/setup.py
+++ b/setup.py
@@ -446,7 +446,6 @@ def run(self):
               'sympy.parsing.latex': ['*.txt', '*.g4'],
               'sympy.integrals.rubi.parsetools': ['header.py.txt'],
               'sympy.plotting.tests': ['test_region_*.png'],
-              'sympy': ['py.typed']
               },
           data_files=[('share/man/man1', ['doc/man/isympy.1'])],
           cmdclass={'test': test_sympy,
diff --git a/sympy/printing/numpy.py b/sympy/printing/numpy.py
index 0a49b031ac..0234f85f3b 100644
--- a/sympy/printing/numpy.py
+++ b/sympy/printing/numpy.py
@@ -1,5 +1,5 @@
 from sympy.core import S
-from .pycode import PythonCodePrinter, _known_functions_math, _print_known_const, _print_known_func, _unpack_integral_limits
+from .pycode import PythonCodePrinter, _known_functions_math, _print_known_const, _print_known_func, _unpack_integral_limits, ArrayPrinter
 from .codeprinter import CodePrinter
 
 
@@ -30,7 +30,7 @@
 _numpy_known_functions = {k: 'numpy.' + v for k, v in _known_functions_numpy.items()}
 _numpy_known_constants = {k: 'numpy.' + v for k, v in _known_constants_numpy.items()}
 
-class NumPyPrinter(PythonCodePrinter):
+class NumPyPrinter(ArrayPrinter, PythonCodePrinter):
     """
     Numpy printer which handles vectorized piecewise functions,
     logical operators, etc.
@@ -250,64 +250,6 @@ def _print_BlockMatrix(self, expr):
         return '{}({})'.format(self._module_format(self._module + '.block'),
                                  self._print(expr.args[0].tolist()))
 
-    def _print_ArrayTensorProduct(self, expr):
-        array_list = [j for i, arg in enumerate(expr.args) for j in
-                (self._print(arg), "[%i, %i]" % (2*i, 2*i+1))]
-        return "%s(%s)" % (self._module_format(self._module + '.einsum'), ", ".join(array_list))
-
-    def _print_ArrayContraction(self, expr):
-        from ..tensor.array.expressions.array_expressions import ArrayTensorProduct
-        base = expr.expr
-        contraction_indices = expr.contraction_indices
-        if not contraction_indices:
-            return self._print(base)
-        if isinstance(base, ArrayTensorProduct):
-            counter = 0
-            d = {j: min(i) for i in contraction_indices for j in i}
-            indices = []
-            for rank_arg in base.subranks:
-                lindices = []
-                for i in range(rank_arg):
-                    if counter in d:
-                        lindices.append(d[counter])
-                    else:
-                        lindices.append(counter)
-                    counter += 1
-                indices.append(lindices)
-            elems = ["%s, %s" % (self._print(arg), ind) for arg, ind in zip(base.args, indices)]
-            return "%s(%s)" % (
-                self._module_format(self._module + '.einsum'),
-                ", ".join(elems)
-            )
-        raise NotImplementedError()
-
-    def _print_ArrayDiagonal(self, expr):
-        diagonal_indices = list(expr.diagonal_indices)
-        if len(diagonal_indices) > 1:
-            # TODO: this should be handled in sympy.codegen.array_utils,
-            # possibly by creating the possibility of unfolding the
-            # ArrayDiagonal object into nested ones. Same reasoning for
-            # the array contraction.
-            raise NotImplementedError
-        if len(diagonal_indices[0]) != 2:
-            raise NotImplementedError
-        return "%s(%s, 0, axis1=%s, axis2=%s)" % (
-            self._module_format("numpy.diagonal"),
-            self._print(expr.expr),
-            diagonal_indices[0][0],
-            diagonal_indices[0][1],
-        )
-
-    def _print_PermuteDims(self, expr):
-        return "%s(%s, %s)" % (
-            self._module_format("numpy.transpose"),
-            self._print(expr.expr),
-            self._print(expr.permutation.array_form),
-        )
-
-    def _print_ArrayAdd(self, expr):
-        return self._expand_fold_binary_op(self._module + '.add', expr.args)
-
     def _print_NDimArray(self, expr):
         if len(expr.shape) == 1:
             return self._module + '.array(' + self._print(expr.args[0]) + ')'
@@ -316,6 +258,12 @@ def _print_NDimArray(self, expr):
         # Should be possible to extend to more dimensions
         return CodePrinter._print_not_supported(self, expr)
 
+    _add = "add"
+    _einsum = "einsum"
+    _transpose = "transpose"
+    _ones = "ones"
+    _zeros = "zeros"
+
     _print_lowergamma = CodePrinter._print_not_supported
     _print_uppergamma = CodePrinter._print_not_supported
     _print_fresnelc = CodePrinter._print_not_supported
diff --git a/sympy/printing/pycode.py b/sympy/printing/pycode.py
index 3a72ab606c..70e714503a 100644
--- a/sympy/printing/pycode.py
+++ b/sympy/printing/pycode.py
@@ -183,41 +183,6 @@ def _expand_reduce_binary_op(self, op, args):
                 self._expand_reduce_binary_op(args[Nhalf:]),
             )
 
-    def _get_einsum_string(self, subranks, contraction_indices):
-        letters = self._get_letter_generator_for_einsum()
-        contraction_string = ""
-        counter = 0
-        d = {j: min(i) for i in contraction_indices for j in i}
-        indices = []
-        for rank_arg in subranks:
-            lindices = []
-            for i in range(rank_arg):
-                if counter in d:
-                    lindices.append(d[counter])
-                else:
-                    lindices.append(counter)
-                counter += 1
-            indices.append(lindices)
-        mapping = {}
-        letters_free = []
-        letters_dum = []
-        for i in indices:
-            for j in i:
-                if j not in mapping:
-                    l = next(letters)
-                    mapping[j] = l
-                else:
-                    l = mapping[j]
-                contraction_string += l
-                if j in d:
-                    if l not in letters_dum:
-                        letters_dum.append(l)
-                else:
-                    letters_free.append(l)
-            contraction_string += ","
-        contraction_string = contraction_string[:-1]
-        return contraction_string, letters_free, letters_dum
-
     def _print_NaN(self, expr):
         return "float('nan')"
 
@@ -422,6 +387,142 @@ def _hprint_Pow(self, expr, rational=False, sqrt='math.sqrt'):
         return "{}**{}".format(base_str, exp_str)
 
 
+class ArrayPrinter:
+
+    def _arrayify(self, indexed):
+        from sympy.tensor.array.expressions.conv_indexed_to_array import convert_indexed_to_array
+        try:
+            return convert_indexed_to_array(indexed)
+        except Exception:
+            return indexed
+
+    def _get_einsum_string(self, subranks, contraction_indices):
+        letters = self._get_letter_generator_for_einsum()
+        contraction_string = ""
+        counter = 0
+        d = {j: min(i) for i in contraction_indices for j in i}
+        indices = []
+        for rank_arg in subranks:
+            lindices = []
+            for i in range(rank_arg):
+                if counter in d:
+                    lindices.append(d[counter])
+                else:
+                    lindices.append(counter)
+                counter += 1
+            indices.append(lindices)
+        mapping = {}
+        letters_free = []
+        letters_dum = []
+        for i in indices:
+            for j in i:
+                if j not in mapping:
+                    l = next(letters)
+                    mapping[j] = l
+                else:
+                    l = mapping[j]
+                contraction_string += l
+                if j in d:
+                    if l not in letters_dum:
+                        letters_dum.append(l)
+                else:
+                    letters_free.append(l)
+            contraction_string += ","
+        contraction_string = contraction_string[:-1]
+        return contraction_string, letters_free, letters_dum
+
+    def _get_letter_generator_for_einsum(self):
+        for i in range(97, 123):
+            yield chr(i)
+        for i in range(65, 91):
+            yield chr(i)
+        raise ValueError("out of letters")
+
+    def _print_ArrayTensorProduct(self, expr):
+        letters = self._get_letter_generator_for_einsum()
+        contraction_string = ",".join(["".join([next(letters) for j in range(i)]) for i in expr.subranks])
+        return '%s("%s", %s)' % (
+                self._module_format(self._module + "." + self._einsum),
+                contraction_string,
+                ", ".join([self._print(arg) for arg in expr.args])
+        )
+
+    def _print_ArrayContraction(self, expr):
+        from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
+        base = expr.expr
+        contraction_indices = expr.contraction_indices
+
+        if isinstance(base, ArrayTensorProduct):
+            elems = ",".join(["%s" % (self._print(arg)) for arg in base.args])
+            ranks = base.subranks
+        else:
+            elems = self._print(base)
+            ranks = [len(base.shape)]
+
+        contraction_string, letters_free, letters_dum = self._get_einsum_string(ranks, contraction_indices)
+
+        if not contraction_indices:
+            return self._print(base)
+        if isinstance(base, ArrayTensorProduct):
+            elems = ",".join(["%s" % (self._print(arg)) for arg in base.args])
+        else:
+            elems = self._print(base)
+        return "%s(\"%s\", %s)" % (
+            self._module_format(self._module + "." + self._einsum),
+            "{}->{}".format(contraction_string, "".join(sorted(letters_free))),
+            elems,
+        )
+
+    def _print_ArrayDiagonal(self, expr):
+        from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
+        diagonal_indices = list(expr.diagonal_indices)
+        if isinstance(expr.expr, ArrayTensorProduct):
+            subranks = expr.expr.subranks
+            elems = expr.expr.args
+        else:
+            subranks = expr.subranks
+            elems = [expr.expr]
+        diagonal_string, letters_free, letters_dum = self._get_einsum_string(subranks, diagonal_indices)
+        elems = [self._print(i) for i in elems]
+        return '%s("%s", %s)' % (
+            self._module_format(self._module + "." + self._einsum),
+            "{}->{}".format(diagonal_string, "".join(letters_free+letters_dum)),
+            ", ".join(elems)
+        )
+
+    def _print_PermuteDims(self, expr):
+        return "%s(%s, %s)" % (
+            self._module_format(self._module + "." + self._transpose),
+            self._print(expr.expr),
+            self._print(expr.permutation.array_form),
+        )
+
+    def _print_ArrayAdd(self, expr):
+        return self._expand_fold_binary_op(self._module + "." + self._add, expr.args)
+
+    def _print_OneArray(self, expr):
+        return "%s((%s,))" % (
+            self._module_format(self._module+ "." + self._ones),
+            ','.join(map(self._print,expr.args))
+        )
+
+    def _print_ZeroArray(self, expr):
+        return "%s((%s,))" % (
+            self._module_format(self._module+ "." + self._zeros),
+            ','.join(map(self._print,expr.args))
+        )
+
+    def _print_Assignment(self, expr):
+        #XXX: maybe this needs to happen at a higher level e.g. at _print or
+        #doprint?
+        lhs = self._print(self._arrayify(expr.lhs))
+        rhs = self._print(self._arrayify(expr.rhs))
+        return "%s = %s" % ( lhs, rhs )
+
+    def _print_IndexedBase(self, expr):
+        return self._print_ArraySymbol(expr)
+
+
 class PythonCodePrinter(AbstractPythonCodePrinter):
 
     def _print_sign(self, e):
diff --git a/sympy/printing/tensorflow.py b/sympy/printing/tensorflow.py
index 0643e7728c..5c3a342bef 100644
--- a/sympy/printing/tensorflow.py
+++ b/sympy/printing/tensorflow.py
@@ -6,12 +6,12 @@
 from sympy.codegen.cfunctions import Sqrt
 from sympy.external import import_module
 from sympy.printing.precedence import PRECEDENCE
-from sympy.printing.pycode import AbstractPythonCodePrinter
+from sympy.printing.pycode import AbstractPythonCodePrinter, ArrayPrinter
 import sympy
 
 tensorflow = import_module('tensorflow')
 
-class TensorflowPrinter(AbstractPythonCodePrinter):
+class TensorflowPrinter(ArrayPrinter, AbstractPythonCodePrinter):
     """
     Tensorflow printer which handles vectorized piecewise functions,
     logical operators, max/min, and relational operators.
@@ -196,13 +196,6 @@ def _print_MatPow(self, expr):
         return self._expand_fold_binary_op(
             "tensorflow.linalg.matmul", [expr.base]*expr.exp)
 
-    def _print_Assignment(self, expr):
-        # TODO: is this necessary?
-        return "%s = %s" % (
-            self._print(expr.lhs),
-            self._print(expr.rhs),
-        )
-
     def _print_CodeBlock(self, expr):
         # TODO: is this necessary?
         ret = []
@@ -210,73 +203,12 @@ def _print_CodeBlock(self, expr):
             ret.append(self._print(subexpr))
         return "\n".join(ret)
 
-    def _get_letter_generator_for_einsum(self):
-        for i in range(97, 123):
-            yield chr(i)
-        for i in range(65, 91):
-            yield chr(i)
-        raise ValueError("out of letters")
-
-    def _print_ArrayTensorProduct(self, expr):
-        letters = self._get_letter_generator_for_einsum()
-        contraction_string = ",".join(["".join([next(letters) for j in range(i)]) for i in expr.subranks])
-        return '%s("%s", %s)' % (
-                self._module_format('tensorflow.linalg.einsum'),
-                contraction_string,
-                ", ".join([self._print(arg) for arg in expr.args])
-        )
-
-    def _print_ArrayContraction(self, expr):
-        from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
-        base = expr.expr
-        contraction_indices = expr.contraction_indices
-        contraction_string, letters_free, letters_dum = self._get_einsum_string(base.subranks, contraction_indices)
-
-        if not contraction_indices:
-            return self._print(base)
-        if isinstance(base, ArrayTensorProduct):
-            elems = ["%s" % (self._print(arg)) for arg in base.args]
-            return "%s(\"%s\", %s)" % (
-                self._module_format("tensorflow.linalg.einsum"),
-                contraction_string,
-                ", ".join(elems)
-            )
-        raise NotImplementedError()
-
-    def _print_ArrayDiagonal(self, expr):
-        from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct
-        diagonal_indices = list(expr.diagonal_indices)
-        if len(diagonal_indices) > 1:
-            # TODO: this should be handled in sympy.codegen.array_utils,
-            # possibly by creating the possibility of unfolding the
-            # ArrayDiagonal object into nested ones. Same reasoning for
-            # the array contraction.
-            raise NotImplementedError
-        if len(diagonal_indices[0]) != 2:
-            raise NotImplementedError
-        if isinstance(expr.expr, ArrayTensorProduct):
-            subranks = expr.expr.subranks
-            elems = expr.expr.args
-        else:
-            subranks = expr.subranks
-            elems = [expr.expr]
-        diagonal_string, letters_free, letters_dum = self._get_einsum_string(subranks, diagonal_indices)
-        elems = [self._print(i) for i in elems]
-        return '%s("%s", %s)' % (
-            self._module_format("tensorflow.linalg.einsum"),
-            "{}->{}{}".format(diagonal_string, "".join(letters_free), "".join(letters_dum)),
-            ", ".join(elems)
-        )
-
-    def _print_PermuteDims(self, expr):
-        return "%s(%s, %s)" % (
-            self._module_format("tensorflow.transpose"),
-            self._print(expr.expr),
-            self._print(expr.permutation.array_form),
-        )
-
-    def _print_ArrayAdd(self, expr):
-        return self._expand_fold_binary_op('tensorflow.math.add', expr.args)
+    _module = "tensorflow"
+    _einsum = "linalg.einsum"
+    _add = "math.add"
+    _transpose = "transpose"
+    _ones = "ones"
+    _zeros = "zeros"
 
 
 def tensorflow_code(expr, **settings):
diff --git a/sympy/py.typed b/sympy/py.typed
deleted file mode 100644
index e69de29bb2..0000000000
diff --git a/sympy/tensor/array/expressions/conv_indexed_to_array.py b/sympy/tensor/array/expressions/conv_indexed_to_array.py
index e88ca2bbd6..94ea8d2c7e 100644
--- a/sympy/tensor/array/expressions/conv_indexed_to_array.py
+++ b/sympy/tensor/array/expressions/conv_indexed_to_array.py
@@ -1,6 +1,5 @@
 from collections import defaultdict
 
-from sympy.combinatorics.permutations import _af_invert
 from sympy.concrete.summations import Sum
 from sympy.core.add import Add
 from sympy.core.mul import Mul
@@ -13,7 +12,7 @@
 from sympy.matrices.expressions.matexpr import MatrixElement
 from sympy.tensor.array.expressions.array_expressions import ArrayDiagonal, \
     get_shape, ArrayElement, _array_tensor_product, _array_diagonal, _array_contraction, _array_add, \
-    _permute_dims, OneArray, ArrayAdd
+    _permute_dims
 from sympy.tensor.array.expressions.utils import _get_argindex, _get_diagonal_indices
 
 
@@ -94,11 +93,8 @@ def _get_pos(elem, indices):
                 return i
         raise ValueError("not found")
 
-    permutation = _af_invert([_get_pos(i, first_indices) for i in indices])
-    if isinstance(result, ArrayAdd):
-        return _array_add(*[_permute_dims(arg, permutation) for arg in result.args])
-    else:
-        return _permute_dims(result, permutation)
+    permutation = [_get_pos(i, first_indices) for i in indices]
+    return _permute_dims(result, permutation)
 
 
 def _convert_indexed_to_array(expr):
@@ -225,24 +221,15 @@ def _convert_indexed_to_array(expr):
         args, indices = zip(*[_convert_indexed_to_array(arg) for arg in expr.args])
         args = list(args)
         # Check if all indices are compatible. Otherwise expand the dimensions:
-        index0 = []
-        shape0 = []
-        for arg, arg_indices in zip(args, indices):
-            arg_indices_set = set(arg_indices)
-            arg_indices_missing = arg_indices_set.difference(index0)
-            index0.extend([i for i in arg_indices if i in arg_indices_missing])
-            arg_shape = get_shape(arg)
-            shape0.extend([arg_shape[i] for i, e in enumerate(arg_indices) if e in arg_indices_missing])
-        for i, (arg, arg_indices) in enumerate(zip(args, indices)):
-            if len(arg_indices) < len(index0):
-                missing_indices_pos = [i for i, e in enumerate(index0) if e not in arg_indices]
-                missing_shape = [shape0[i] for i in missing_indices_pos]
-                arg_indices = tuple(index0[j] for j in missing_indices_pos) + arg_indices
-                args[i] = _array_tensor_product(OneArray(*missing_shape), args[i])
-            permutation = Permutation([arg_indices.index(j) for j in index0])
+        index0set = set(indices[0])
+        index0 = indices[0]
+        for i in range(1, len(args)):
+            if set(indices[i]) != index0set:
+                raise NotImplementedError("indices must be the same")
+            permutation = Permutation([index0.index(j) for j in indices[i]])
             # Perform index permutations:
             args[i] = _permute_dims(args[i], permutation)
-        return _array_add(*args), tuple(index0)
+        return _array_add(*args), index0
     if isinstance(expr, Pow):
         subexpr, subindices = _convert_indexed_to_array(expr.base)
         if isinstance(expr.exp, (int, Integer)):
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index cfb156b7c9..0df60a2cd5 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -9,12 +9,14 @@
 from operator import gt
 
 
-
 # this is the logical location of these functions
+# from sympy.core.compatibility import ordered
+# from sympy.core.compatibility import default_sort_key  # noqa: F401
+from sympy.utilities.misc import as_int
+
 from sympy.utilities.enumerative import (
     multiset_partitions_taocp, list_visitor, MultisetPartitionTraverser)
 
-from sympy.utilities.misc import as_int
 from sympy.utilities.decorator import deprecated
 
 
@@ -1943,10 +1945,7 @@ def has_dups(seq):
     if isinstance(seq, (dict, set, Dict, Set)):
         return False
     unique = set()
-    try:
-        return any(True for s in seq if s in unique or unique.add(s))
-    except TypeError:
-        return len(seq) != len(list(uniq(seq)))
+    return any(True for s in seq if s in unique or unique.add(s))
 
 
 def has_variety(seq):
@@ -2182,105 +2181,64 @@ def multiset_derangements(s):
     of derangements is desired or else all values will be the same:
 
     >>> list(uniq([i for i in multiset_derangements('1233')]))
-    [[None, None, None, None]]
+    [['3', '3', '2', '1']]
     >>> [i.copy() for i in multiset_derangements('1233')]
     [['3', '3', '1', '2'], ['3', '3', '2', '1']]
-    >>> [''.join(i) for i in multiset_derangements('1233')]
-    ['3312', '3321']
     """
-    from sympy.core.sorting import ordered
-    # create multiset dictionary of hashable elements or else
-    # remap elements to integers
-    try:
-        ms = multiset(s)
-    except TypeError:
-        # give each element a canonical integer value
-        key = dict(enumerate(ordered(uniq(s))))
-        h = []
-        for si in s:
-            for k in key:
-                if key[k] == si:
-                    h.append(k)
-                    break
-        for i in multiset_derangements(h):
-            yield [key[j] for j in i]
-        return
-
-    mx = max(ms.values())  # max repetition of any element
-    n = len(s)  # the number of elements
+    ms = multiset(s)
+    mx = max(ms.values())
+    n = len(s)
+    # special cases
 
-    ## special cases
-
-    # 1) one element has more than half the total cardinality of s: no
-    # derangements are possible.
+    # 0) impossible case
     if mx*2 > n:
         return
 
-    # 2) all elements appear once: singletons
+    # 1) singletons
     if len(ms) == n:
-        yield from _set_derangements(s)
+        for p in generate_derangements(s):
+            yield p
         return
 
-    # find the first element that is repeated the most to place
-    # in the following two special cases where the selection
-    # is unambiguous: either there are two elements with multiplicity
-    # of mx or else there is only one with multiplicity mx
     for M in ms:
         if ms[M] == mx:
             break
-
-    inonM = [i for i in range(n) if s[i] != M]  # location of non-M
-    iM = [i for i in range(n) if s[i] == M]  # locations of M
+    inonM = [i for i in range(n) if s[i] != M]
+    iM = [i for i in range(n) if s[i] == M]
     rv = [None]*n
 
-    # 3) half are the same
+    # 2) half are the same
     if 2*mx == n:
-        # M goes into non-M locations
         for i in inonM:
             rv[i] = M
-        # permutations of non-M go to M locations
         for p in multiset_permutations([s[i] for i in inonM]):
             for i, pi in zip(iM, p):
                 rv[i] = pi
             yield rv
-        # clean-up (and encourages proper use of routine)
-        rv[:] = [None]*n
         return
 
-    # 4) single repeat covers all but 1 of the non-repeats:
-    # if there is one repeat then the multiset of the values
-    # of ms would be {mx: 1, 1: n - mx}, i.e. there would
-    # be n - mx + 1 values with the condition that n - 2*mx = 1
-    if n - 2*mx == 1 and len(ms.values()) == n - mx + 1:
+    # 3) single repeat covers all but 1 of the non-repeats
+    if n - 2*mx == 1 and len(ms.values()) - 1 == n - mx:
         for i in range(len(inonM)):
             i1 = inonM[i]
             ifill = inonM[:i] + inonM[i+1:]
             for j in ifill:
                 rv[j] = M
+            rv[i1] = s[i1]
             for p in permutations([s[j] for j in ifill]):
-                rv[i1] = s[i1]
                 for j, pi in zip(iM, p):
                     rv[j] = pi
-                k = i1
                 for j in iM:
-                    rv[j], rv[k] = rv[k], rv[j]
+                    rv[j], rv[i1] = rv[i1], rv[j]
                     yield rv
-                    k = j
-        # clean-up (and encourages proper use of routine)
-        rv[:] = [None]*n
+                    i1 = j
         return
 
-    ## general case is handled with 3 helpers:
-    #    1) `finish_derangements` will place the last two elements
-    #       which have arbitrary multiplicities, e.g. for multiset
-    #       {c: 3, a: 2, b: 2}, the last two elements are a and b
-    #    2) `iopen` will tell where a given element can be placed
-    #    3) `do` will recursively place elements into subsets of
-    #        valid locations
-
     def finish_derangements():
         """Place the last two elements into the partially completed
         derangement, and yield the results.
+        In non-recursive version, this will be inlined, but a little
+        easier to understand as a function for now.
         """
 
         a = take[1][0]  # penultimate element
@@ -2305,12 +2263,11 @@ def finish_derangements():
         if len(forced_a) > a_ct or len(forced_b) > b_ct:
             # No derangement possible
             return
-
         for i in forced_a:
             rv[i] = a
         for i in forced_b:
             rv[i] = b
-        for a_place in combinations(open_free, a_ct - len(forced_a)):
+        for a_place in subsets(open_free, a_ct - len(forced_a)):
             for a_pos in a_place:
                 rv[a_pos] = a
             for i in open_free:
@@ -2320,7 +2277,6 @@ def finish_derangements():
             # Clean up/undo the final placements
             for i in open_free:
                 rv[i] = None
-
         # additional cleanup - clear forced_a, forced_b
         for i in forced_a:
             rv[i] = None
@@ -2328,35 +2284,21 @@ def finish_derangements():
             rv[i] = None
 
     def iopen(v):
-        # return indices at which element v can be placed in rv:
-        # locations which are not already occupied if that location
-        # does not already contain v in the same location of s
         return [i for i in range(n) if rv[i] is None and s[i] != v]
 
     def do(j):
-        if j == 1:
-            # handle the last two elements (regardless of multiplicity)
-            # with a special method
-            yield from finish_derangements()
+        if j == -1:
+            yield rv
         else:
-            # place the mx elements of M into a subset of places
-            # into which it can be replaced
             M, mx = take[j]
-            for i in combinations(iopen(M), mx):
-                # place M
+            for i in subsets(iopen(M), mx):
                 for ii in i:
                     rv[ii] = M
-                # recursively place the next element
                 yield from do(j - 1)
-                # mark positions where M was placed as once again
-                # open for placement of other elements
                 for ii in i:
                     rv[ii] = None
-
-    # process elements in order of canonically decreasing multiplicity
     take = sorted(ms.items(), key=lambda x:(x[1], x[0]))
     yield from do(len(take) - 1)
-    rv[:] = [None]*n
 
 
 def random_derangement(t, choice=None, strict=True):
@@ -2456,28 +2398,9 @@ def pick(rv, n):
     return rv
 
 
-def _set_derangements(s):
-    """
-    yield derangements of items in ``s`` which are assumed to contain
-    no repeated elements
+def generate_derangements(perm):
     """
-    if len(s) < 2:
-        return
-    if len(s) == 2:
-        yield [s[1], s[0]]
-        return
-    if len(s) == 3:
-        yield [s[1], s[2], s[0]]
-        yield [s[2], s[0], s[1]]
-        return
-    for p in permutations(s):
-        if not any(i == j for i, j in zip(p, s)):
-            yield list(p)
-
-
-def generate_derangements(s):
-    """
-    Return unique derangements of the elements of iterable ``s``.
+    Routine to generate unique derangements or sets or multisets.
 
     Examples
     ========
@@ -2485,8 +2408,10 @@ def generate_derangements(s):
     >>> from sympy.utilities.iterables import generate_derangements
     >>> list(generate_derangements([0, 1, 2]))
     [[1, 2, 0], [2, 0, 1]]
-    >>> list(generate_derangements([0, 1, 2, 2]))
-    [[2, 2, 0, 1], [2, 2, 1, 0]]
+    >>> list(generate_derangements([0, 1, 2, 3]))
+    [[1, 0, 3, 2], [1, 2, 3, 0], [1, 3, 0, 2], [2, 0, 3, 1], \
+    [2, 3, 0, 1], [2, 3, 1, 0], [3, 0, 1, 2], [3, 2, 0, 1], \
+    [3, 2, 1, 0]]
     >>> list(generate_derangements([0, 1, 1]))
     []
 
@@ -2496,10 +2421,20 @@ def generate_derangements(s):
     sympy.functions.combinatorial.factorials.subfactorial
 
     """
-    if not has_dups(s):
-        yield from _set_derangements(s)
+    if not has_dups(perm):
+        s = perm
+        if len(perm) == 2:
+            yield [s[1],s[0]]
+            return
+        if len(perm) == 3:
+            yield [s[1],s[2],s[0]]
+            yield [s[2],s[0],s[1]]
+            return
+        for p in permutations(s):
+            if not any(i == j for i, j in zip(p, s)):
+                yield list(p)
     else:
-        for p in multiset_derangements(s):
+        for p in multiset_derangements(perm):
             yield list(p)
 
 


Test patch:
diff --git a/sympy/printing/tests/test_pycode.py b/sympy/printing/tests/test_pycode.py
index 9c1285d7dd..4f02ea38de 100644
--- a/sympy/printing/tests/test_pycode.py
+++ b/sympy/printing/tests/test_pycode.py
@@ -12,9 +12,11 @@
 from sympy.printing.pycode import (
     MpmathPrinter, PythonCodePrinter, pycode, SymPyPrinter
 )
+from sympy.printing.tensorflow import TensorflowPrinter
 from sympy.printing.numpy import NumPyPrinter, SciPyPrinter
 from sympy.testing.pytest import raises, skip
-from sympy.tensor import IndexedBase
+from sympy.tensor import IndexedBase, Idx
+from sympy.tensor.array.expressions.array_expressions import ArraySymbol, ArrayDiagonal, ArrayContraction, ZeroArray, OneArray
 from sympy.external import import_module
 from sympy.functions.special.gamma_functions import loggamma
 from sympy.parsing.latex import parse_latex
@@ -217,7 +219,6 @@ def test_frac():
     from sympy.functions.elementary.integers import frac
 
     expr = frac(x)
-
     prntr = NumPyPrinter()
     assert prntr.doprint(expr) == 'numpy.mod(x, 1)'
 
@@ -392,3 +393,28 @@ def test_numerical_accuracy_functions():
     assert prntr.doprint(expm1(x)) == 'numpy.expm1(x)'
     assert prntr.doprint(log1p(x)) == 'numpy.log1p(x)'
     assert prntr.doprint(cosm1(x)) == 'scipy.special.cosm1(x)'
+
+def test_array_printer():
+    A = ArraySymbol('A', (4,4,6,6,6))
+    I = IndexedBase('I')
+    i,j,k = Idx('i', (0,1)), Idx('j', (2,3)), Idx('k', (4,5))
+
+    prntr = NumPyPrinter()
+    assert prntr.doprint(ZeroArray(5)) == 'numpy.zeros((5,))'
+    assert prntr.doprint(OneArray(5)) == 'numpy.ones((5,))'
+    assert prntr.doprint(ArrayContraction(A, [2,3])) == 'numpy.einsum("abccd->abd", A)'
+    assert prntr.doprint(I) == 'I'
+    assert prntr.doprint(ArrayDiagonal(A, [2,3,4])) == 'numpy.einsum("abccc->abc", A)'
+    assert prntr.doprint(ArrayDiagonal(A, [0,1], [2,3])) == 'numpy.einsum("aabbc->cab", A)'
+    assert prntr.doprint(ArrayContraction(A, [2], [3])) == 'numpy.einsum("abcde->abe", A)'
+    assert prntr.doprint(Assignment(I[i,j,k], I[i,j,k])) == 'I = I'
+
+    prntr = TensorflowPrinter()
+    assert prntr.doprint(ZeroArray(5)) == 'tensorflow.zeros((5,))'
+    assert prntr.doprint(OneArray(5)) == 'tensorflow.ones((5,))'
+    assert prntr.doprint(ArrayContraction(A, [2,3])) == 'tensorflow.linalg.einsum("abccd->abd", A)'
+    assert prntr.doprint(I) == 'I'
+    assert prntr.doprint(ArrayDiagonal(A, [2,3,4])) == 'tensorflow.linalg.einsum("abccc->abc", A)'
+    assert prntr.doprint(ArrayDiagonal(A, [0,1], [2,3])) == 'tensorflow.linalg.einsum("aabbc->cab", A)'
+    assert prntr.doprint(ArrayContraction(A, [2], [3])) == 'tensorflow.linalg.einsum("abcde->abe", A)'
+    assert prntr.doprint(Assignment(I[i,j,k], I[i,j,k])) == 'I = I'
diff --git a/sympy/tensor/array/expressions/tests/test_convert_index_to_array.py b/sympy/tensor/array/expressions/tests/test_convert_index_to_array.py
index a9cc2f382e..419c335ecd 100644
--- a/sympy/tensor/array/expressions/tests/test_convert_index_to_array.py
+++ b/sympy/tensor/array/expressions/tests/test_convert_index_to_array.py
@@ -7,7 +7,7 @@
 from sympy.combinatorics import Permutation
 from sympy.tensor.array.expressions.array_expressions import ArrayContraction, ArrayTensorProduct, \
     ArrayDiagonal, ArrayAdd, PermuteDims, ArrayElement, _array_tensor_product, _array_contraction, _array_diagonal, \
-    _array_add, _permute_dims, ArraySymbol, OneArray
+    _array_add, _permute_dims, ArraySymbol
 from sympy.tensor.array.expressions.conv_array_to_matrix import convert_array_to_matrix
 from sympy.tensor.array.expressions.conv_indexed_to_array import convert_indexed_to_array, _convert_indexed_to_array
 from sympy.testing.pytest import raises
@@ -15,7 +15,6 @@
 
 A, B = symbols("A B", cls=IndexedBase)
 i, j, k, l, m, n = symbols("i j k l m n")
-d0, d1, d2, d3 = symbols("d0:4")
 
 I = Identity(k)
 
@@ -168,32 +167,3 @@ def test_arrayexpr_convert_indexed_to_array_out_of_bounds():
     raises(ValueError, lambda: convert_indexed_to_array(expr))
     expr = Sum(M[i, j]*N[j,m], (j, 1, k-1))
     raises(ValueError, lambda: convert_indexed_to_array(expr))
-
-
-def test_arrayexpr_convert_indexed_to_array_broadcast():
-    A = ArraySymbol("A", (3, 3))
-    B = ArraySymbol("B", (3, 3))
-
-    expr = A[i, j] + B[k, l]
-    O2 = OneArray(3, 3)
-    expected = ArrayAdd(ArrayTensorProduct(A, O2), ArrayTensorProduct(O2, B))
-    assert convert_indexed_to_array(expr) == expected
-    assert convert_indexed_to_array(expr, [i, j, k, l]) == expected
-    assert convert_indexed_to_array(expr, [l, k, i, j]) == ArrayAdd(PermuteDims(ArrayTensorProduct(O2, A), [1, 0, 2, 3]), PermuteDims(ArrayTensorProduct(B, O2), [1, 0, 2, 3]))
-
-    expr = A[i, j] + B[j, k]
-    O1 = OneArray(3)
-    assert convert_indexed_to_array(expr, [i, j, k]) == ArrayAdd(ArrayTensorProduct(A, O1), ArrayTensorProduct(O1, B))
-
-    C = ArraySymbol("C", (d0, d1))
-    D = ArraySymbol("D", (d3, d1))
-
-    expr = C[i, j] + D[k, j]
-    assert convert_indexed_to_array(expr, [i, j, k]) == ArrayAdd(ArrayTensorProduct(C, OneArray(d3)), PermuteDims(ArrayTensorProduct(OneArray(d0), D), [0, 2, 1]))
-
-    X = ArraySymbol("X", (5, 3))
-
-    expr = X[i, n] - X[j, n]
-    assert convert_indexed_to_array(expr, [i, j, n]) == ArrayAdd(ArrayTensorProduct(-1, OneArray(5), X), PermuteDims(ArrayTensorProduct(X, OneArray(5)), [0, 2, 1]))
-
-    raises(ValueError, lambda: convert_indexed_to_array(C[i, j] + D[i, j]))
diff --git a/sympy/utilities/tests/test_iterables.py b/sympy/utilities/tests/test_iterables.py
index 0c391d507c..c2c22e6295 100644
--- a/sympy/utilities/tests/test_iterables.py
+++ b/sympy/utilities/tests/test_iterables.py
@@ -560,11 +560,6 @@ def test_derangements():
         'cccabba', 'cccabab', 'cccaabb', 'ccacbba', 'ccacbab',
         'ccacabb', 'cbccbaa', 'cbccaba', 'cbccaab', 'bcccbaa',
         'bcccaba', 'bcccaab']
-    assert [''.join(i) for i in D('books')] == ['kbsoo', 'ksboo',
-        'sbkoo', 'skboo', 'oksbo', 'oskbo', 'okbso', 'obkso', 'oskob',
-        'oksob', 'osbok', 'obsok']
-    assert list(generate_derangements([[3], [2], [2], [1]])) == [
-        [[2], [1], [3], [2]], [[2], [3], [1], [2]]]
 
 
 def test_necklaces():
@@ -799,8 +794,6 @@ def test_has_dups():
     assert has_dups(set()) is False
     assert has_dups(list(range(3))) is False
     assert has_dups([1, 2, 1]) is True
-    assert has_dups([[1], [1]]) is True
-    assert has_dups([[1], [2]]) is False
 
 
 def test__partition():