Repository: sympy/sympy
Pull request: Fix for numpy arrays when the size 1 arrays were treated as scalars
#### References to other Issues or PRs
Revives and closes #14720 
Fixes #14706

#### Brief description of what is fixed or changed

--> Added a bug fix in sympify whereby the numpy arrays of size 1 won't be treated as scalars

#### Other comments

Example:
```
>>> from sympy import *
>>> from sympy.abc import *
>>> import numpy
>>> z = np.zeros((1,1))
```

Before the fix:
```
>>> f = x + z
x
```

After the fix:
```
>>> f = x + z
array([[x]], dtype=object)
```

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* core
  *  sympify no longer converts numpy arrays with size to scalars
<!-- END RELEASE NOTES -->

Problem statement:
Symbol violates the commutativity of addition when it sometimes converts one-element arrays to scalar
```
from sympy import symbols
import numpy as np

for i in [2, 1]:
    z = np.zeros(i)
    xz = symbols('x') + z
    print(z.shape)
    print(xz.shape)
    print()
```
prints:
```
(2,)
(2,)

(1,)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-165-14d0a2118753> in <module>()
      6     xz = symbols('x') + z
      7     print(z.shape)
----> 8     print(xz.shape)

AttributeError: 'Symbol' object has no attribute 'shape'
```

A scalar is  *definitely* not a numpy array with shape `(1,)`, `(1, 1)`, etc.
It's debatable whether a scalar is a numpy array with shape `()`, but I would simply return a numpy array in this case too.  Every library that has ever made the mistake of doing this dimension contraction had to eventually undo it (e.g., Theano).

The problem is that it's very easy to get the lone element from an array of size 1:
```
In [1]: o = np.ones([1 for _ in range(10)])

In [2]: o
Out[2]: array([[[[[[[[[[1.]]]]]]]]]])

In [3]: o.item()
Out[3]: 1.0
```
But it's very hard to go from an element and try to figure out what the shape it should have had and then give it the right shape.

Since this breaks my code, I can submit a patch within about a day if it will be accepted and if someone can point me to the right place in the code (I'm guessing tensor.py).

Hints from issue discussion:
Basic SymPy objects don't carry a shape attribute, `ImmutableDenseNDimArray` does however, have you had a look at: http://docs.sympy.org/latest/modules/tensor/array.html#module-sympy.tensor.array ?
Thanks for taking a look!  Unfortunately, I don't think `Array` helps me.

I have a function, say `f`, that transforms numpy arrays (say `A, B, C`) into one vector `v1`.  I want to pass in symbolic versions of those arrays to `f` to produce `v1`.  Let `v2` be a list of unique elements of `v1`.  From `v1` and `v2`, I want to find a matrix `M` that maps `v1` to `v2`.

Each time `A, B, C` are passed (this time they are regular numpy arrays of numbers rather than symbolic variables) I don't need to store `f(A, B, C)`, but rather `M @ f(A, B, C)`, which is shorter.  In my case, it is much shorter.

Any inconsistencies between numpy and sympy complicate `f` greatly.

How is `Array` better than `ndarray` having `dtype=np.object` and containing `Symbol`?  It seems to me that the `ndarray` will do the usual `numpy` broadcasting that `f` does with regular `ndarray`s?  That seems to me to be better.  This issue is about an inconsistency in sympy, which was apparently added as a feature (sympy.tensor.tensor.py, line 2258):

```
"""
    When all indices are contracted and components data are added to the tensor,
    accessing the data will return a scalar, no array object. In fact, arrays
    are dropped to scalars if they contain only one element.

    >>> A(i0, -i0)
    A(L_0, -L_0)
    >>> A(i0, -i0).data
    -18
"""
```

At least in my case, this is not a feature.
SymPy expressions aren't designed to interact directly with NumPy arrays. Generally if you want to work with NumPy arrays from SymPy expressions you should `lambdify` the expression first to create a compatible NumPy function. 
I don't need SymPy expressions to be "designed to interact with numpy arrays".  SymPy doesn't need to know anything about numpy arrays.  In fact, the whole reason this isn't working is because SymPy seems to know something about numpy arrays.

The numpy array is just like a rectangular list of objects.  For example, if I have:
```
xs = np.array(symbols('x:5'))
ys = np.array(symbols('y:3'))
```
Then I can add all pairs like this:
```
np.add.outer(xs, ys)
```
to give
```
array([[x0 + y0, x0 + y1, x0 + y2],
       [x1 + y0, x1 + y1, x1 + y2],
       [x2 + y0, x2 + y1, x2 + y2],
       [x3 + y0, x3 + y1, x3 + y2],
       [x4 + y0, x4 + y1, x4 + y2]], dtype=object)
```

This works with any object that supports addition (e.g., `Fraction`, `Decimal`, and `str`):
```
from fractions import Fraction
xs = np.array([Fraction(i) for i in range(1, 5)])
ys = np.array([Fraction(i) for i in range(6, 9)])
np.divide.outer(xs, ys)
```
gives:
```
array([[Fraction(1, 6), Fraction(1, 7), Fraction(1, 8)],
       [Fraction(1, 3), Fraction(2, 7), Fraction(1, 4)],
       [Fraction(1, 2), Fraction(3, 7), Fraction(3, 8)],
       [Fraction(2, 3), Fraction(4, 7), Fraction(1, 2)]], dtype=object)
```

The problem is that in some special cases when the shape of the numpy array contains one element (above), sympy is discarding shape information.
Playing with it a bit more, I see that SymPy is violating the commutativity of addition:
```
z = np.zeros((1, 1))
z + symbols('x')
```
gives the right answer: `array([[x]], dtype=object)`. But
```
symbols('x') + z
```
gives `x`!!

This suggests that `__add__` is handling addition with numpy arrays but `__radd__` is not.  Ideally, neither method would do anything with numpy arrays.  Just let numpy handle the addition in both cases.
Both SymPy and NumPy override the operators to try to "absorb" other objects into SymPy or NumPy (resp.) types. We could make `Symbol.__add__` return `NotImplemented` when it sees a NumPy array, but that could potentially break code that expects SymPy to convert NumPy arrays into SymPy NDimArray objects (the SymPy operators all just call `sympify`).

My best recommendation here is if you want to put SymPy objects (or any foreign objects, really) in NumPy arrays to always do it explicitly. Don't rely on NumPy's operators to create arrays automatically, because that will fail when the foreign objects have their own greedy operators. 

With that being said, perhaps making SymPy's NDimArray objects not automatically flatten would be a good change. Maybe @Upabjojr can comment more on why that behavior is there. 
> My best recommendation here is if you want to put SymPy objects (or any foreign objects, really) in NumPy arrays to always do it explicitly. 

How can I do the addition (`z + symbols(x)`) explicitly?  Using `np.broadcast_arrays` first?  It seems ugly.

> With that being said, perhaps making SymPy's NDimArray objects not automatically flatten would be a good change. Maybe @Upabjojr can comment more on why that behavior is there.

Thank you.  This is what I'm hoping for.
> How can I do the addition (z + symbols(x)) explicitly? Using np.broadcast_arrays first? It seems ugly.

Won't `z + numpy.array(x)` do what you want? I don't use NumPy that often, so I may be forgetting some of the details of broadcasting. 
> Won't z + numpy.array(x) do what you want?

Yup.

Nevertheless, I'd like to feel safe adding signals and arrays without worrying about order.  Therefore, I'm hoping that SymPy will keep addition commutative.

Also, since we can see that `ndarray + signal` does the right thing, I don't see why there should be a backwards compatibility question about changing `signal + ndarray` to match.
Addition is commutative in SymPy. If you mix NumPy and SymPy objects, you could have unexpected behaviours (such as your example).

Commit message: core: More test cases for sympify

Particularly related to the new changes, new tests were added to
increase the coverage.

Patch to review:
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 8ff2149dcf03..6933aae2b7d7 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -7,7 +7,8 @@
 import re as regex
 
 from .containers import Tuple
-from .sympify import converter, sympify, _sympify, SympifyError, _convert_numpy_types
+from .sympify import (SympifyError, converter, sympify, _convert_numpy_types, _sympify,
+                      _is_numpy_instance)
 from .singleton import S, Singleton
 from .expr import Expr, AtomicExpr
 from .evalf import pure_complex
@@ -1086,7 +1087,7 @@ def __new__(cls, num, dps=None, prec=None, precision=None):
             return num
         elif num is S.NaN:
             return num
-        elif type(num).__module__ == 'numpy': # support for numpy datatypes
+        elif _is_numpy_instance(num):  # support for numpy datatypes
             num = _convert_numpy_types(num)
         elif isinstance(num, mpmath.mpf):
             if precision is None:
diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py
index 29ec61491d45..e04b8112d1c3 100644
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -55,6 +55,16 @@ class CantSympify(object):
     pass
 
 
+def _is_numpy_instance(a):
+    """
+    Checks if an object is an instance of a type from the numpy module.
+    """
+    # This check avoids unnecessarily importing NumPy.  We check the whole
+    # __mro__ in case any base type is a numpy type.
+    return any(type_.__module__ == 'numpy'
+               for type_ in type(a).__mro__)
+
+
 def _convert_numpy_types(a, **sympify_args):
     """
     Converts a numpy datatype input to an appropriate SymPy type.
@@ -295,8 +305,7 @@ def sympify(a, locals=None, convert_xor=True, strict=False, rational=False,
         evaluate = global_parameters.evaluate
 
     # Support for basic numpy datatypes
-    # Note that this check exists to avoid importing NumPy when not necessary
-    if type(a).__module__ == 'numpy':
+    if _is_numpy_instance(a):
         import numpy as np
         if np.isscalar(a):
             return _convert_numpy_types(a, locals=locals,
@@ -324,15 +333,30 @@ def sympify(a, locals=None, convert_xor=True, strict=False, rational=False,
                 return Array(a.flat, a.shape)  # works with e.g. NumPy arrays
 
     if not isinstance(a, str):
-        for coerce in (float, int):
-            try:
-                coerced = coerce(a)
-            except (TypeError, ValueError):
-                continue
-            try:
-                return sympify(coerced)
-            except SympifyError:
-                continue
+        if _is_numpy_instance(a):
+            import numpy as np
+            assert not isinstance(a, np.number)
+            if isinstance(a, np.ndarray):
+                # Scalar arrays (those with zero dimensions) have sympify
+                # called on the scalar element.
+                if a.ndim == 0:
+                    try:
+                        return sympify(a.item(),
+                                       locals=locals,
+                                       convert_xor=convert_xor,
+                                       strict=strict,
+                                       rational=rational,
+                                       evaluate=evaluate)
+                    except SympifyError:
+                        pass
+        else:
+            # float and int can coerce size-one numpy arrays to their lone
+            # element.  See issue https://github.com/numpy/numpy/issues/10404.
+            for coerce in (float, int):
+                try:
+                    return sympify(coerce(a))
+                except (TypeError, ValueError, AttributeError, SympifyError):
+                    continue
 
     if strict:
         raise SympifyError(a)
diff --git a/sympy/core/tests/test_sympify.py b/sympy/core/tests/test_sympify.py
index e2d25e47b8a4..1e4e96958eed 100644
--- a/sympy/core/tests/test_sympify.py
+++ b/sympy/core/tests/test_sympify.py
@@ -632,6 +632,7 @@ def test_issue_13924():
     assert isinstance(a, ImmutableDenseNDimArray)
     assert a[0] == 1
 
+
 def test_numpy_sympify_args():
     # Issue 15098. Make sure sympify args work with numpy types (like numpy.str_)
     if not numpy:
@@ -695,3 +696,39 @@ def test_issue_16759():
 def test_issue_17811():
     a = Function('a')
     assert sympify('a(x)*5', evaluate=False) == Mul(a(x), 5, evaluate=False)
+
+
+def test_issue_14706():
+    if not numpy:
+        skip("numpy not installed.")
+
+    from sympy import symbols
+
+    x = symbols('x')
+
+    z1 = numpy.zeros((1,1), dtype=numpy.float)
+    z2 = numpy.zeros((2,2), dtype=numpy.float)
+    z3 = numpy.zeros((), dtype=numpy.float)
+
+    assert numpy.all(x + z1 == numpy.full((1, 1), x))
+    assert numpy.all(x + z2 == numpy.full((2, 2), x))
+    for z in [z3,
+              numpy.int(0),
+              numpy.float(0),
+              numpy.complex(0)]:
+        assert x + z == x
+        assert isinstance(x + z, Symbol)
+
+    assert x + numpy.array(x) == 2 * x
+    assert x + numpy.array([x]) == numpy.array([2*x], dtype=object)
+
+    assert sympify(numpy.array([1])) == ImmutableDenseNDimArray([1], 1)
+    assert sympify(numpy.array([[[1]]])) == ImmutableDenseNDimArray([1], (1, 1, 1))
+    assert sympify(z1) == ImmutableDenseNDimArray([0], (1, 1))
+    assert sympify(z2) == ImmutableDenseNDimArray([0, 0, 0, 0], (2, 2))
+    assert sympify(z3) == ImmutableDenseNDimArray([0], ())
+    assert sympify(z3, strict=True) == 0.0
+
+    raises(SympifyError, lambda: sympify(numpy.array([1]), strict=True))
+    raises(SympifyError, lambda: sympify(z1, strict=True))
+    raises(SympifyError, lambda: sympify(z2, strict=True))


Code review comments:
diff_hunk: @@ -695,3 +696,39 @@ def test_issue_16759():
 def test_issue_17811():
     a = Function('a')
     assert sympify('a(x)*5', evaluate=False) == Mul(a(x), 5, evaluate=False)
+
+
+def test_issue_14706():
+    if not numpy:
+        skip("numpy not installed.")
+
+    from sympy import symbols
+
+    x = symbols('x')
path: sympy/core/tests/test_sympify.py
text: @user1:
```suggestion
```

x is already imported at the top of the file