Repository: sympy/sympy

Issue:
RandomMatrixSymbol doesn't respect pspace for substitution
```python3
from sympy import *
from sympy.stats import GaussianSymplecticEnsemble as GSE, density
from sympy.stats.random_matrix import RandomMatrixPSpace
from sympy.stats.rv import RandomMatrixSymbol

X = MatrixSymbol('X', 2, 2)

G = GSE('U', 2)
h_pspace = RandomMatrixPSpace('P', model=density(G))
H = RandomMatrixSymbol('H', 2, 2, pspace=h_pspace)
H2 = RandomMatrixSymbol('H', 2, 2, pspace=None)
(H*2).xreplace({H2: X})
```
I have tackled to this problem in #19715, but I think that they should be different objects when different `pspace` are used and I think that this is a side effect of inheriting redefined `_hashable_content` from `MatrixSymbol` which only takes account of name and shape.

Discussion:

Patch:
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 1443b296d8..87a40b0b57 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -3,7 +3,7 @@
 from functools import wraps, reduce
 import collections
 
-from sympy.core import S, Symbol, Tuple, Integer, Basic, Expr, Mul, Add
+from sympy.core import S, Symbol, Integer, Basic, Expr, Mul, Add
 from sympy.core.decorators import call_highest_priority
 from sympy.core.compatibility import SYMPY_INTS, default_sort_key
 from sympy.core.sympify import SympifyError, _sympify
@@ -776,25 +776,14 @@ def __new__(cls, name, n, m):
         obj = Basic.__new__(cls, name, n, m)
         return obj
 
-    def _hashable_content(self):
-        return (self.name, self.shape)
-
     @property
     def shape(self):
-        return self.args[1:3]
+        return self.args[1], self.args[2]
 
     @property
     def name(self):
         return self.args[0].name
 
-    def _eval_subs(self, old, new):
-        # only do substitutions in shape
-        shape = Tuple(*self.shape)._subs(old, new)
-        return MatrixSymbol(self.args[0], *shape)
-
-    def __call__(self, *args):
-        raise TypeError("%s object is not callable" % self.__class__)
-
     def _entry(self, i, j, **kwargs):
         return MatrixElement(self, i, j)
 
@@ -802,13 +791,6 @@ def _entry(self, i, j, **kwargs):
     def free_symbols(self):
         return {self}
 
-    def doit(self, **hints):
-        if hints.get('deep', True):
-            return type(self)(self.args[0], self.args[1].doit(**hints),
-                    self.args[2].doit(**hints))
-        else:
-            return self
-
     def _eval_simplify(self, **kwargs):
         return self
 


Test patch:
diff --git a/sympy/stats/tests/test_random_matrix.py b/sympy/stats/tests/test_random_matrix.py
index 5cfd7dcfcc..87f9bca09a 100644
--- a/sympy/stats/tests/test_random_matrix.py
+++ b/sympy/stats/tests/test_random_matrix.py
@@ -1,5 +1,5 @@
 from sympy import (sqrt, exp, Trace, pi, S, Integral, MatrixSymbol, Lambda,
-                    Dummy, Product, Abs, IndexedBase, Matrix, I, Rational)
+                   Dummy, Product, Abs, IndexedBase, Matrix, I, Rational)
 from sympy.stats import (GaussianUnitaryEnsemble as GUE, density,
                          GaussianOrthogonalEnsemble as GOE,
                          GaussianSymplecticEnsemble as GSE,
@@ -12,7 +12,7 @@
                          Normal, Beta)
 from sympy.stats.joint_rv_types import JointDistributionHandmade
 from sympy.stats.rv import RandomMatrixSymbol
-from sympy.stats.random_matrix_models import GaussianEnsemble
+from sympy.stats.random_matrix_models import GaussianEnsemble, RandomMatrixPSpace
 from sympy.testing.pytest import raises
 
 def test_GaussianEnsemble():
@@ -110,3 +110,15 @@ def test_issue_19841():
     G1 = GUE('U', 2)
     G2 = G1.xreplace({2: 2})
     assert G1.args == G2.args
+
+    X = MatrixSymbol('X', 2, 2)
+    G = GSE('U', 2)
+    h_pspace = RandomMatrixPSpace('P', model=density(G))
+    H = RandomMatrixSymbol('H', 2, 2, pspace=h_pspace)
+    H2 = RandomMatrixSymbol('H', 2, 2, pspace=None)
+    assert H.doit() == H
+
+    assert (2*H).xreplace({H: X}) == 2*X
+    assert (2*H).xreplace({H2: X}) == 2*H
+    assert (2*H2).xreplace({H: X}) == 2*H2
+    assert (2*H2).xreplace({H2: X}) == 2*X