Repository: sympy/sympy
Pull request: MatrixElement symbol property added; validate shape compatible with indices
#### Brief description of what is fixed or changed

fixes #22676

##### symbol attribute has been added
MatrixElement behaves more like Indexed in terms of argument access:

```python
>>> IndexedBase('x')[0].base
x
>>> MatrixSymbol('x',2,2)[0,0].symbol
x
```

##### invalid substitutions no longer succeed

In master,
```python
X = MatrixSymbol("X", 2, 2)
Y = MatrixSymbol("Y", 1, 2)
assert X[1, 1].subs(X, Y) == Y[1, 1]
>>> Y[1,1]
Traceback (most recent call last):
...
IndexError: indices out of bounds
```
in this PR that substitution no longer succeeds so you don't end up with an invalid expression

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->


#### Other comments

There is some divergence in how MatrixSymbol/IndexedBase and MatrixElement/Indexed behave in terms of free symbols and attributes. This should have a look by someone that works with matrix expressions to see if these changes make sense.

- [x] QUESTION: is there are reason to store Indexed.label as a Symbol but MatrixSymbol as a Str?
ANSWER: no (see #22391 and #22272)

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* matrices
    * replacement of MatrixSymbol in MatrixElement requires shape to be compatible with indices to succeed
    * MatrixElement now has `symbol` attributes (analogous to `base` attribute of `IndexedBase` )
    * `eval_hyper_sum` returns as unchanged arguments for which `is_hypergeometric` is False and those which produce singularities at `a` or `b + 1` where `and` and `b` are the lower and upper limits of evaluation, respectively.
* stats
    * removed deprecation warning about using `evaluate=False`
<!-- END RELEASE NOTES -->


Problem statement:
MatrixElement allows bad subs
```python
>>> m22 = MatrixSymbol("A",2,2)
>>> m11 = MatrixSymbol("B",1,1)
>>> m22.subs(m22,m11) == m11  # this is ok, just changing a symbol without context
True
>>> m22[1,1].subs(m22, m11)  # this is not ok -- it indicates a non-existent element in B
B[1, 1]
```
If you enter the element manually you will get the error:
```python
>>> B[1, 1]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
KeyError: (1, 1)
```
Such a substitution should either fail to change the object or else raise a KeyError.

Hints from issue discussion:

Commit message: watch for illegals at both limits

Patch to review:
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index 799667109e3c..f2c74da5eaaf 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -28,6 +28,7 @@
 from sympy.polys.partfrac import apart
 from sympy.polys.polyerrors import PolynomialError, PolificationFailed
 from sympy.polys.polytools import parallel_poly_from_expr, Poly, factor
+from sympy.polys.polyutils import illegal
 from sympy.polys.rationaltools import together
 from sympy.series.limitseq import limit_seq
 from sympy.series.order import O
@@ -1295,6 +1296,9 @@ def _eval_sum_hyper(f, i, a):
 def eval_sum_hyper(f, i_a_b):
     i, a, b = i_a_b
 
+    if f.is_hypergeometric(i) is False:
+        return
+
     if (b - a).is_Integer:
         # We are never going to do better than doing the sum in the obvious way
         return None
@@ -1307,10 +1311,15 @@ def eval_sum_hyper(f, i_a_b):
             if res is not None:
                 return Piecewise(res, (old_sum, True))
         else:
+            n_illegal = lambda x: sum(x.count(_) for _ in illegal)
+            had = n_illegal(f)
+            # check that no extra illegals are introduced
             res1 = _eval_sum_hyper(f, i, a)
+            if res1 is None or n_illegal(res1) > had:
+                return
             res2 = _eval_sum_hyper(f, i, b + 1)
-            if res1 is None or res2 is None:
-                return None
+            if res2 is None or n_illegal(res2) > had:
+                return
             (res1, cond1), (res2, cond2) = res1, res2
             cond = And(cond1, cond2)
             if cond == False:
diff --git a/sympy/concrete/tests/test_gosper.py b/sympy/concrete/tests/test_gosper.py
index 70cc488cb844..fc81b70769d7 100644
--- a/sympy/concrete/tests/test_gosper.py
+++ b/sympy/concrete/tests/test_gosper.py
@@ -8,8 +8,8 @@
 from sympy.functions.special.gamma_functions import gamma
 from sympy.polys.polytools import Poly
 from sympy.simplify.simplify import simplify
-from sympy.abc import a, b, j, k, m, n, r, x
 from sympy.concrete.gosper import gosper_normal, gosper_sum, gosper_term
+from sympy.abc import a, b, j, k, m, n, r, x
 
 
 def test_gosper_normal():
diff --git a/sympy/concrete/tests/test_sums_products.py b/sympy/concrete/tests/test_sums_products.py
index 10d7d14e0d80..93e924250596 100644
--- a/sympy/concrete/tests/test_sums_products.py
+++ b/sympy/concrete/tests/test_sums_products.py
@@ -29,7 +29,6 @@
 from sympy.simplify.combsimp import combsimp
 from sympy.simplify.simplify import simplify
 from sympy.tensor.indexed import (Idx, Indexed, IndexedBase)
-from sympy.abc import a, b, c, d, k, m, x, y, z
 from sympy.concrete.summations import (
     telescopic, _dummy_with_inherited_properties_concrete, eval_sum_residue)
 from sympy.concrete.expr_with_intlimits import ReorderError
@@ -38,6 +37,7 @@
 from sympy.matrices import (Matrix, SparseMatrix,
     ImmutableDenseMatrix, ImmutableSparseMatrix)
 from sympy.core.mod import Mod
+from sympy.abc import a, b, c, d, k, m, x, y, z
 
 n = Symbol('n', integer=True)
 f, g = symbols('f g', cls=Function)
@@ -1186,8 +1186,9 @@ def test_issue_14640():
 
 def test_issue_15943():
     s = Sum(binomial(n, k)*factorial(n - k), (k, 0, n)).doit().rewrite(gamma)
-    assert s == -E*(n + 1)*gamma(n + 1)*lowergamma(n + 1, 1)/gamma(n + 2
-        ) + E*gamma(n + 1)
+    assert s == Sum(gamma(n + 1)/gamma(k + 1), (k, 0, n))
+    s = s.doit()
+    assert s == (-E*(n + 1)*lowergamma(n + 1, 1)/factorial(n + 1) + E)*gamma(n + 1)
     assert s.simplify() == E*(factorial(n) - lowergamma(n + 1, 1))
 
 
@@ -1586,3 +1587,10 @@ def test_process_limits():
         raises(TypeError, lambda: D(x, x > 0))
         raises(ValueError, lambda: D(x, Interval(1, 3)))
         raises(NotImplementedError, lambda: D(x, (x, union)))
+
+
+def test_pr_22677():
+    b = Symbol('b', integer=True, positive=True)
+    assert Sum(1/x**2,(x, 0, b)).doit() == Sum(x**(-2), (x, 0, b))
+    assert Sum(1/(x - b)**2,(x, 0, b-1)).doit() == Sum(
+        (-b + x)**(-2), (x, 0, b - 1))
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 78db834466e6..ccbf3588a00a 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -274,8 +274,8 @@ def is_valid(idx):
             return isinstance(idx, (int, Integer, Symbol, Expr))
         return (is_valid(i) and is_valid(j) and
                 (self.rows is None or
-                (0 <= i) != False and (i < self.rows) != False) and
-                (0 <= j) != False and (j < self.cols) != False)
+                (i >= -self.rows) != False and (i < self.rows) != False) and
+                (j >= -self.cols) != False and (j < self.cols) != False)
 
     def __getitem__(self, key):
         if not isinstance(key, tuple) and isinstance(key, slice):
@@ -589,10 +589,14 @@ class MatrixElement(Expr):
     def __new__(cls, name, n, m):
         n, m = map(_sympify, (n, m))
         from sympy.matrices.matrices import MatrixBase
-        if isinstance(name, (MatrixBase,)):
+        if isinstance(name, MatrixBase):
             if n.is_Integer and m.is_Integer:
                 return name[n, m]
-        if isinstance(name, str):
+            name = _sympify(name)  # change mutable into immutable
+        elif isinstance(name, MatrixSymbol):
+            if not name.valid_index(n, m):
+                raise IndexError('indices out of range')
+        elif isinstance(name, str):
             name = Symbol(name)
         else:
             name = _sympify(name)
@@ -601,6 +605,10 @@ def __new__(cls, name, n, m):
         obj = Expr.__new__(cls, name, n, m)
         return obj
 
+    @property
+    def symbol(self):
+        return self.args[0]
+
     def doit(self, **kwargs):
         deep = kwargs.get('deep', True)
         if deep:
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index ff0c30e8537d..7f21fb1016bf 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -2,7 +2,7 @@
 from sympy.core.exprtools import gcd_terms
 from sympy.core.function import (diff, expand)
 from sympy.core.relational import Eq
-from sympy.core.symbol import (Dummy, Symbol)
+from sympy.core.symbol import (Dummy, Symbol, Str)
 from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.matrices.dense import zeros
 from sympy.polys.polytools import factor
@@ -48,10 +48,13 @@ def test_matrix_symbol_creation():
     raises(ValueError, lambda: MatrixSymbol('A', n, n))
 
 
-def test_shape():
+def test_matexpr_properties():
     assert A.shape == (n, m)
     assert (A*B).shape == (n, l)
     raises(ShapeError, lambda: B*A)
+    assert A[0, 1].indices == (0, 1)
+    assert A[0, 0].symbol == A
+    assert A[0, 0].symbol.name == 'A'
 
 
 def test_matexpr():
@@ -61,7 +64,7 @@ def test_matexpr():
     assert (A*B).shape == (n, l)
 
 
-def test_subs():
+def test_matexpr_subs():
     A = MatrixSymbol('A', n, m)
     B = MatrixSymbol('B', m, l)
     C = MatrixSymbol('C', m, l)
@@ -70,6 +73,30 @@ def test_subs():
     assert (A*B).subs(B, C) == A*C
     assert (A*B).subs(l, n).is_square
 
+    W = MatrixSymbol("W", 3, 3)
+    X = MatrixSymbol("X", 2, 2)
+    Y = MatrixSymbol("Y", 1, 2)
+    Z = MatrixSymbol("Z", n, 2)
+    # no restrictions on Symbol replacement
+    assert X.subs(X, Y) == Y
+    # it might be better to just change the name
+    y = Str('y')
+    assert X.subs(Str("X"), y).args == (y, 2, 2)
+    # it's ok to introduce a wider matrix
+    assert X[1, 1].subs(X, W) == W[1, 1]
+    # but for a given MatrixExpression, only change
+    # name if indexing on the new shape is valid.
+    # Here, X is 2,2; Y is 1,2 and Y[1, 1] is out
+    # of range so an error is raised
+    raises(IndexError, lambda: X[1, 1].subs(X, Y))
+    # here, [0, 1] is in range so the subs succeeds
+    assert X[0, 1].subs(X, Y) == Y[0, 1]
+    # and here the size of n will accept any index
+    # in the first position
+    assert W[2, 1].subs(W, Z) == Z[2, 1]
+    # but not in the second position
+    raises(IndexError, lambda: W[2, 2].subs(W, Z))
+
     A = SparseMatrix([[1, 2], [3, 4]])
     B = Matrix([[1, 2], [3, 4]])
     C, D = MatrixSymbol('C', 2, 2), MatrixSymbol('D', 2, 2)
@@ -187,11 +214,15 @@ def test_invariants():
         assert obj == obj.__class__(*obj.args)
 
 
-def test_indexing():
+def test_matexpr_indexing():
     A = MatrixSymbol('A', n, m)
     A[1, 2]
     A[l, k]
-    A[l+1, k+1]
+    A[l + 1, k + 1]
+    A = MatrixSymbol('A', 2, 1)
+    for i in range(-2, 2):
+        for j in range(-1, 1):
+            A[i, j]
 
 
 def test_single_indexing():
diff --git a/sympy/matrices/expressions/tests/test_trace.py b/sympy/matrices/expressions/tests/test_trace.py
index f00bd31a8181..3bd66bec2377 100644
--- a/sympy/matrices/expressions/tests/test_trace.py
+++ b/sympy/matrices/expressions/tests/test_trace.py
@@ -8,6 +8,8 @@
 )
 from sympy.matrices.expressions.special import OneMatrix
 from sympy.testing.pytest import raises
+from sympy.abc import i
+
 
 n = symbols('n', integer=True)
 A = MatrixSymbol('A', n, n)
@@ -95,8 +97,9 @@ def test_trace_constant_factor():
     assert trace(MatMul(2, X)) == 10
 
 
-def test_rewrite():
-    assert isinstance(trace(A).rewrite(Sum), Sum)
+def test_trace_rewrite():
+    assert trace(A).rewrite(Sum) == Sum(A[i, i], (i, 0, n - 1))
+    assert trace(eye(3)).rewrite(Sum) == 3
 
 
 def test_trace_normalize():
diff --git a/sympy/matrices/expressions/trace.py b/sympy/matrices/expressions/trace.py
index edecdd4b3924..8b462a5c4a5f 100644
--- a/sympy/matrices/expressions/trace.py
+++ b/sympy/matrices/expressions/trace.py
@@ -2,7 +2,7 @@
 from sympy.core.expr import Expr, ExprBuilder
 from sympy.core.singleton import S
 from sympy.core.sorting import default_sort_key
-from sympy.core.symbol import Dummy
+from sympy.core.symbol import uniquely_named_symbol
 from sympy.core.sympify import sympify
 from sympy.matrices.matrices import MatrixBase
 from sympy.matrices.common import NonSquareMatrixError
@@ -144,8 +144,9 @@ def get_arg_key(x):
 
     def _eval_rewrite_as_Sum(self, expr, **kwargs):
         from sympy.concrete.summations import Sum
-        i = Dummy('i')
-        return Sum(self.arg[i, i], (i, 0, self.arg.rows-1)).doit()
+        i = uniquely_named_symbol('i', expr)
+        s = Sum(self.arg[i, i], (i, 0, self.arg.rows - 1))
+        return s.doit()
 
 
 def trace(expr):
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 34384d4e2bc8..31691e1a75c9 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -39,11 +39,9 @@
 from sympy.sets.sets import FiniteSet, ProductSet, Intersection
 from sympy.solvers.solveset import solveset
 from sympy.external import import_module
-from sympy.utilities.misc import filldedent
 from sympy.utilities.decorator import doctest_depends_on
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import iterable
-import warnings
 
 
 x = Symbol('x')
@@ -841,11 +839,6 @@ def probability(condition, given_condition=None, numsamples=None,
     from sympy.stats.symbolic_probability import Probability
     if evaluate:
         return Probability(condition, given_condition).doit(**kwargs)
-    ### TODO: Remove the user warnings in the future releases
-    message = ("Since version 1.7, using `evaluate=False` returns `Probability` "
-              "object. If you want unevaluated Integral/Sum use "
-              "`P(condition, given_condition, evaluate=False).rewrite(Integral)`")
-    warnings.warn(filldedent(message))
     return Probability(condition, given_condition)
 
 
diff --git a/sympy/stats/tests/test_discrete_rv.py b/sympy/stats/tests/test_discrete_rv.py
index 0e3abb6d48c8..b011872abb98 100644
--- a/sympy/stats/tests/test_discrete_rv.py
+++ b/sympy/stats/tests/test_discrete_rv.py
@@ -1,7 +1,7 @@
 from sympy.concrete.summations import Sum
 from sympy.core.numbers import (I, Rational, oo, pi)
 from sympy.core.singleton import S
-from sympy.core.symbol import (Dummy, Symbol)
+from sympy.core.symbol import Symbol
 from sympy.functions.elementary.complexes import (im, re)
 from sympy.functions.elementary.exponential import log
 from sympy.functions.elementary.integers import floor
@@ -24,7 +24,7 @@
                                    FlorySchulz, Poisson, Geometric, Hermite, Logarithmic,
                                     NegativeBinomial, Skellam, YuleSimon, Zeta,
                                     DiscreteRV)
-from sympy.testing.pytest import slow, nocache_fail, raises, ignore_warnings
+from sympy.testing.pytest import slow, nocache_fail, raises
 from sympy.stats.symbolic_probability import Expectation
 
 x = Symbol('x')
@@ -43,14 +43,15 @@ def test_Poisson():
     l = 3
     x = Poisson('x', l)
     assert E(x) == l
+    assert E(2*x) == 2*l
     assert variance(x) == l
     assert density(x) == PoissonDistribution(l)
-    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
-        assert isinstance(E(x, evaluate=False), Expectation)
-        assert isinstance(E(2*x, evaluate=False), Expectation)
+    assert isinstance(E(x, evaluate=False), Expectation)
+    assert isinstance(E(2*x, evaluate=False), Expectation)
     # issue 8248
     assert x.pspace.compute_expectation(1) == 1
 
+
 def test_FlorySchulz():
     a = Symbol("a")
     z = Symbol("z")
@@ -107,8 +108,7 @@ def test_Logarithmic():
     assert E(x) == -p / ((1 - p) * log(1 - p))
     assert variance(x) == -1/log(2)**2 + 2/log(2)
     assert E(2*x**2 + 3*x + 4) == 4 + 7 / log(2)
-    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
-        assert isinstance(E(x, evaluate=False), Expectation)
+    assert isinstance(E(x, evaluate=False), Expectation)
 
 
 @nocache_fail
@@ -120,8 +120,7 @@ def test_negative_binomial():
     # This hangs when run with the cache disabled:
     assert variance(x) == p*r / (1-p)**2
     assert E(x**5 + 2*x + 3) == Rational(9207, 4)
-    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
-        assert isinstance(E(x, evaluate=False), Expectation)
+    assert isinstance(E(x, evaluate=False), Expectation)
 
 
 def test_skellam():
@@ -148,8 +147,7 @@ def test_yule_simon():
     x = YuleSimon('x', rho)
     assert simplify(E(x)) == rho / (rho - 1)
     assert simplify(variance(x)) == rho**2 / ((rho - 1)**2 * (rho - 2))
-    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
-        assert isinstance(E(x, evaluate=False), Expectation)
+    assert isinstance(E(x, evaluate=False), Expectation)
     # To test the cdf function
     assert cdf(x)(x) == Piecewise((-beta(floor(x), 4)*floor(x) + 1, x >= 1), (0, True))
 
@@ -291,17 +289,12 @@ def test_conditional():
 def test_product_spaces():
     X1 = Geometric('X1', S.Half)
     X2 = Geometric('X2', Rational(1, 3))
-    #assert str(P(X1 + X2 < 3, evaluate=False)) == """Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, """\
-    #    + """(-X2 + n + 3 >= 1) & (-X2 + n + 3 < oo)), (0, True)), (X2, 1, oo), (n, -oo, -1))"""
-    n = Dummy('n')
-    with ignore_warnings(UserWarning): ### TODO: Restore tests once warnings are removed
-        assert P(X1 + X2 < 3, evaluate=False).rewrite(Sum).dummy_eq(Sum(Piecewise((2**(-n)/4,
-         n + 2 >= 1), (0, True)), (n, -oo, -1))/3)
-    #assert str(P(X1 + X2 > 3)) == """Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, """ +\
-    #    """(-X2 + n + 3 >= 1) & (-X2 + n + 3 < oo)), (0, True)), (X2, 1, oo), (n, 1, oo))"""
-    assert P(X1 + X2 > 3).dummy_eq(Sum(Piecewise((2**(X2 - n - 2)*(Rational(2, 3))**(X2 - 1)/6,
-                                                 -X2 + n + 3 >= 1), (0, True)),
-                                       (X2, 1, oo), (n, 1, oo)))
-#    assert str(P(Eq(X1 + X2, 3))) == """Sum(Piecewise((2**(X2 - 2)*(2/3)**(X2 - 1)/6, """ +\
-#        """X2 <= 2), (0, True)), (X2, 1, oo))"""
+    assert str(P(X1 + X2 < 3).rewrite(Sum)) == (
+        "Sum(Piecewise((1/(4*2**n), n >= -1), (0, True)), (n, -oo, -1))/3")
+    assert str(P(X1 + X2 > 3).rewrite(Sum)) == (
+        "Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, "
+        "X2 - n <= 2), (0, True)), (X2, 1, oo), (n, 1, oo))")
+    assert str(P(X1 + X2 > 3).rewrite(Sum)) == (
+        "Sum(Piecewise((2**(X2 - n - 2)*(2/3)**(X2 - 1)/6, "
+        "X2 - n <= 2), (0, True)), (X2, 1, oo), (n, 1, oo))")
     assert P(Eq(X1 + X2, 3)) == Rational(1, 12)


Code review comments:
diff_hunk: @@ -589,10 +589,14 @@ class MatrixElement(Expr):
     def __new__(cls, name, n, m):
         n, m = map(_sympify, (n, m))
         from sympy.matrices.matrices import MatrixBase
-        if isinstance(name, (MatrixBase,)):
+        if isinstance(name, MatrixBase):
             if n.is_Integer and m.is_Integer:
                 return name[n, m]
-        if isinstance(name, str):
+            name = _sympify(name)  # change mutable into immutable
+        elif isinstance(name, MatrixSymbol):
+            if not name.valid_index(n, m):
path: sympy/matrices/expressions/matexpr.py
text: @user1:
I think that `valid_index` should always be checked, not just for MatrixSymbol e.g.:
```python
In [18]: n = Symbol('n', integer=True, positive=True)

In [19]: M = MatrixSymbol('M', n, n)

In [20]: M[3, m]
Out[20]: (M)[3, m]

In [21]: M[3, m].subs(M, Matrix([[1, 2], [3, 4]]))
Out[21]: 
⎛⎡1  2⎤⎞      
⎜⎢    ⎥⎟[3, m]
⎝⎣3  4⎦⎠ 
```

@author:
good point! To me this seems like a strange constructor in that you are creating an unevaluated indexed object rather than simply indexing an object. It works, but having to put these checks here seems odd to me.