Repository: sympy/sympy
Pull request: Add 'evaluate' option to MatAdd and MatMul
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #18764

#### Brief description of what is fixed or changed
'evaluate' option is added to MatAdd and MatMul, which canonicalizes the object.

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- matrices
    - Added `evaluate` option to ``MatAdd`` and ``MatMul``
<!-- END RELEASE NOTES -->

Problem statement:
Introduce 'evaluate' keyword argument in MatAdd and MatMul
``MatAdd`` and ``MatMul`` are behaving differently from ``Add`` and ``Mul``.

Here is an example:
```python
>>> from sympy import Add, MatAdd, Symbol, MatrixSymbol
>>> x = Symbol('x')
>>> Add(x, x, evaluate=True)
2*x
>>> Add(x, x, evaluate=False)
x + x

>> A = MatrixSymbol('A', 2,2)
>>> MatAdd(A, A)
A + A
>>> MatAdd(A, A, evaluate=True)
A + A
```

I believe it would be better to add ``evaluate`` option which canonicallizes the object, so that

```python
>>> MatAdd(A, A, evaluate=True)
2*A
```

Hints from issue discussion:
I had seen https://github.com/sympy/sympy/issues/16484#issuecomment-477878619 that things were not originally designed to use any automatic evaluation.
But I think that things can look different at this moment, and I also didn't like the inconsistency there before. But this may need some agreements because this is more of a design issue.
After `evaluate` is introduced in matadd, matmul, it should also be added in matpow or hadamard stuff.
@sylee957 
I added it to matpow and hadamard.

If sympy's matrices are designed that way, then setting `evaluate`'s default value as `False` honours this.
This is similar to `Derivative`: since `Derivative` should not be automatically evaluated, its default `evaluate` is `False`.

`evaluate` is used everywhere in sympy core. Since this change does not break backward compatibility, I strongly suggest that it is necessary.
Handling evaluate consistently is important. I think it should be done in general across all of sympy as in #17372.

I don't see why we should change the default to `evaluate=True` for MatAdd etc though. You can use `a+b` or `MatAdd(a, b).doit()` if you want it to evaluate.

Commit message: Improve code quality

Patch to review:
diff --git a/sympy/core/parameters.py b/sympy/core/parameters.py
index 3c5ad8d6be5f..4db5fe371417 100644
--- a/sympy/core/parameters.py
+++ b/sympy/core/parameters.py
@@ -65,7 +65,11 @@ def __setattr__(self, name, value):
             clear_cache()
         return super(_global_parameters, self).__setattr__(name, value)
 
-global_parameters = _global_parameters(evaluate=True, distribute=True)
+global_parameters = _global_parameters(\
+    evaluate=True,\
+    distribute=True,\
+    evaluate_matrix=False\
+)
 
 @contextmanager
 def evaluate(x):
@@ -126,3 +130,30 @@ def distribute(x):
         yield
     finally:
         global_parameters.distribute = old
+
+
+@contextmanager
+def evaluate_matrix(x):
+    """ Control automatic evaluation of MatExpr
+
+    Classes in matrices module are designed to prefer unevaluated result.
+
+    Examples
+    ========
+
+    >>> from sympy import MatrixSymbol, MatAdd
+    >>> from sympy.core.parameters import evaluate_matrix
+    >>> A = MatrixSymbol('A', 2,2)
+    >>> print(MatAdd(A,-A))
+    -A + A
+    >>> with evaluate_matrix(True):
+    ...     print(MatAdd(A,-A))
+    0
+    """
+    old = global_parameters.evaluate_matrix
+
+    try:
+        global_parameters.evaluate_matrix = x
+        yield
+    finally:
+        global_parameters.evaluate_matrix = old
diff --git a/sympy/core/tests/test_parameters.py b/sympy/core/tests/test_parameters.py
index 84885533e3ff..f152aadfea3f 100644
--- a/sympy/core/tests/test_parameters.py
+++ b/sympy/core/tests/test_parameters.py
@@ -1,7 +1,7 @@
 from sympy.abc import x, y
-from sympy.core.parameters import evaluate
+from sympy.core.parameters import evaluate, evaluate_matrix
 from sympy.core import Mul, Add, Pow, S
-from sympy import sqrt, oo
+from sympy import sqrt, oo, MatrixSymbol, MatAdd
 
 def test_add():
     with evaluate(False):
@@ -82,3 +82,10 @@ def test_nested():
         expr = (x + x) + (y + y)
         assert expr.args == ((x + x), (y + y))
         assert expr.args[0].args == (x, x)
+
+def test_MatAdd():
+    A = MatrixSymbol('A', 2,2)
+    assert MatAdd(A,A).args == (A, A)
+    with evaluate_matrix(True):
+        expr = MatAdd(A,A)
+        assert expr.args == (2, A)
diff --git a/sympy/matrices/expressions/hadamard.py b/sympy/matrices/expressions/hadamard.py
index 06ab396a498d..8e68f4fad0af 100644
--- a/sympy/matrices/expressions/hadamard.py
+++ b/sympy/matrices/expressions/hadamard.py
@@ -1,6 +1,7 @@
 from __future__ import print_function, division
 
 from sympy.core import Mul, sympify
+from sympy.core.parameters import global_parameters
 from sympy.matrices.expressions.matexpr import (
     MatrixExpr, ShapeError, OneMatrix, ZeroMatrix
 )
@@ -60,13 +61,20 @@ class HadamardProduct(MatrixExpr):
     """
     is_HadamardProduct = True
 
-    def __new__(cls, *args, **kwargs):
+    def __new__(cls, *args, evaluate=None, **kwargs):
         args = list(map(sympify, args))
         check = kwargs.get('check', True)
         if check:
             validate(*args)
 
-        return super(HadamardProduct, cls).__new__(cls, *args)
+        obj = super(HadamardProduct, cls).__new__(cls, *args)
+
+        if evaluate is None:
+            evaluate = global_parameters.evaluate_matrix
+
+        if evaluate:
+            obj = obj.doit(deep=False)
+        return obj
 
     @property
     def shape(self):
@@ -80,7 +88,7 @@ def _eval_transpose(self):
         return HadamardProduct(*list(map(transpose, self.args)))
 
     def doit(self, **ignored):
-        expr = self.func(*[i.doit(**ignored) for i in self.args])
+        expr = self.func(*[i.doit(**ignored) for i in self.args], evaluate=False)
         # Check for explicit matrices:
         from sympy import MatrixBase
         from sympy.matrices.immutable import ImmutableMatrix
@@ -279,7 +287,7 @@ def absorb(x):
             else:
                 new_arg.append(HadamardPower(base, exp))
 
-        x = HadamardProduct(*new_arg)
+        x = HadamardProduct(*new_arg, evaluate=False)
 
     # Commutativity
     fun = condition(
diff --git a/sympy/matrices/expressions/matadd.py b/sympy/matrices/expressions/matadd.py
index 84f826683b9c..fef9311063f6 100644
--- a/sympy/matrices/expressions/matadd.py
+++ b/sympy/matrices/expressions/matadd.py
@@ -4,6 +4,7 @@
 from operator import add
 
 from sympy.core import Add, Basic, sympify
+from sympy.core.parameters import global_parameters
 from sympy.functions import adjoint
 from sympy.matrices.matrices import MatrixBase
 from sympy.matrices.expressions.transpose import transpose
@@ -33,7 +34,7 @@ class MatAdd(MatrixExpr, Add):
 
     identity = GenericZeroMatrix()
 
-    def __new__(cls, *args, **kwargs):
+    def __new__(cls, *args, evaluate=None, **kwargs):
         if not args:
             return cls.identity
 
@@ -44,10 +45,20 @@ def __new__(cls, *args, **kwargs):
         check = kwargs.get('check', False)
 
         obj = Basic.__new__(cls, *args)
+
         if check:
             if all(not isinstance(i, MatrixExpr) for i in args):
                 return Add.fromiter(args)
             validate(*args)
+
+        if evaluate is None:
+            evaluate = global_parameters.evaluate_matrix
+
+        if evaluate:
+            if all(not isinstance(i, MatrixExpr) for i in args):
+                return Add(*args, evaluate=True)
+            obj = canonicalize(obj)
+
         return obj
 
     @property
diff --git a/sympy/matrices/expressions/matmul.py b/sympy/matrices/expressions/matmul.py
index 9da77056932d..5bdefadf0411 100644
--- a/sympy/matrices/expressions/matmul.py
+++ b/sympy/matrices/expressions/matmul.py
@@ -2,6 +2,7 @@
 
 from sympy import Number
 from sympy.core import Mul, Basic, sympify, S
+from sympy.core.parameters import global_parameters
 from sympy.functions import adjoint
 from sympy.strategies import (rm_id, unpack, typed, flatten, exhaust,
         do_one, new)
@@ -34,7 +35,7 @@ class MatMul(MatrixExpr, Mul):
 
     identity = GenericIdentity()
 
-    def __new__(cls, *args, **kwargs):
+    def __new__(cls, *args, evaluate=None, **kwargs):
         check = kwargs.get('check', True)
 
         if not args:
@@ -46,13 +47,22 @@ def __new__(cls, *args, **kwargs):
         args = list(map(sympify, args))
         obj = Basic.__new__(cls, *args)
         factor, matrices = obj.as_coeff_matrices()
+
         if check:
             validate(*matrices)
+
         if not matrices:
             # Should it be
             #
             # return Basic.__neq__(cls, factor, GenericIdentity()) ?
             return factor
+
+        if evaluate is None:
+            evaluate = global_parameters.evaluate_matrix
+
+        if evaluate:
+            return canonicalize(obj)
+
         return obj
 
     @property
@@ -111,7 +121,7 @@ def as_coeff_matrices(self):
 
     def as_coeff_mmul(self):
         coeff, matrices = self.as_coeff_matrices()
-        return coeff, MatMul(*matrices)
+        return coeff, MatMul(*matrices, evaluate=False)
 
     def _eval_transpose(self):
         """Transposition of matrix multiplication.
@@ -355,7 +365,7 @@ def combine_permutations(mul):
         else:
             result.append(B)
 
-    return MatMul(*result)
+    return MatMul(*result, evaluate=False)
 
 rules = (
     any_zeros, remove_ids, combine_powers, unpack, rm_id(lambda x: x == 1),
diff --git a/sympy/matrices/expressions/matpow.py b/sympy/matrices/expressions/matpow.py
index 9d3cc529f161..e1498373cc4c 100644
--- a/sympy/matrices/expressions/matpow.py
+++ b/sympy/matrices/expressions/matpow.py
@@ -2,6 +2,7 @@
 
 from .matexpr import MatrixExpr, ShapeError, Identity, ZeroMatrix
 from sympy.core import S
+from sympy.core.parameters import global_parameters
 from sympy.core.sympify import _sympify
 from sympy.matrices import MatrixBase
 
@@ -10,12 +11,21 @@
 
 class MatPow(MatrixExpr):
 
-    def __new__(cls, base, exp):
+    def __new__(cls, base, exp, evaluate=None, **options):
         base = _sympify(base)
         if not base.is_Matrix:
             raise TypeError("Function parameter should be a matrix")
         exp = _sympify(exp)
-        return super(MatPow, cls).__new__(cls, base, exp)
+
+        obj = super(MatPow, cls).__new__(cls, base, exp)
+
+        if evaluate is None:
+            evaluate = global_parameters.evaluate_matrix
+
+        if evaluate:
+            obj = obj.doit(deep=False)
+
+        return obj
 
     @property
     def base(self):
@@ -84,7 +94,7 @@ def doit(self, **kwargs):
             return Inverse(base).doit(**kwargs)
         elif exp is S.One:
             return base
-        return MatPow(base, exp)
+        return MatPow(base, exp, evaluate=False)
 
     def _eval_transpose(self):
         base, exp = self.args
diff --git a/sympy/matrices/expressions/tests/test_matadd.py b/sympy/matrices/expressions/tests/test_matadd.py
index c7ce61761e30..5a2c9bdf8e9e 100644
--- a/sympy/matrices/expressions/tests/test_matadd.py
+++ b/sympy/matrices/expressions/tests/test_matadd.py
@@ -7,6 +7,9 @@
 X = MatrixSymbol('X', 2, 2)
 Y = MatrixSymbol('Y', 2, 2)
 
+def test_evaluate():
+    assert MatAdd(X, X, evaluate=True) == MatAdd(X, X).doit()
+
 def test_sort_key():
     assert MatAdd(Y, X).doit().args == (X, Y)
 
diff --git a/sympy/matrices/expressions/tests/test_matmul.py b/sympy/matrices/expressions/tests/test_matmul.py
index 7d93e2f4e88d..e0d7b5333460 100644
--- a/sympy/matrices/expressions/tests/test_matmul.py
+++ b/sympy/matrices/expressions/tests/test_matmul.py
@@ -19,6 +19,8 @@
 D = MatrixSymbol('D', n, n)
 E = MatrixSymbol('E', m, n)
 
+def test_evaluate():
+    assert MatMul(C, C, evaluate=True) == MatMul(C, C).doit()
 
 def test_adjoint():
     assert adjoint(A*B) == Adjoint(B)*Adjoint(A)


Code review comments:
diff_hunk: @@ -65,7 +65,11 @@ def __setattr__(self, name, value):
             clear_cache()
         return super(_global_parameters, self).__setattr__(name, value)
 
-global_parameters = _global_parameters(evaluate=True, distribute=True)
+global_parameters = _global_parameters(\
+    evaluate=True,\
+    distribute=True,\
+    evaluate_matrix=False\
path: sympy/core/parameters.py
text: @user1:
These global flags are problematic and already don't work properly. We should not add new global flags.

In any case the argument for having evaluate does not apply to evaluate_matrix because evaluate is used for *disabling* automatic evaluation. If you want to make the matrix expressions evaluate you just have to call `.doit()`.

@user2:
I agree with Oscar, but as another note, the backslashes here are not needed and should be omitted. 

@author:
You're right. I will remove `global_parameters.evaluate_matrix`.