Repository: sympy/sympy

Issue:
error in calculating between an UGate operator and well-defined Unitary operators.
Hi,
I have a problem in calculating between an UGate operator and well-defined Unitary operators in  sympy-1.1.1-py3.6.

I tested two codes as the following:

```
>>> from sympy import *
>>> from sympy.physics.quantum import *
>>> from sympy.physics.quantum.qubit import Qubit
>>> from sympy.physics.quantum.gate import UGate,H
>>> u1 = UGate(0, ImmutableMatrix([[1/sqrt(2),1/sqrt(2)],[I/sqrt(2),-I/sqrt(2)]]))
>>> represent(H(0)*u1,nqubits=1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/represent.py", line 207, in represent
    result = represent(expr.args[-1], **options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/represent.py", line 145, in represent
    return expr._represent(**options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/qexpr.py", line 393, in _represent
    result = self._format_represent(result, format)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/qexpr.py", line 398, in _format_represent
    return to_sympy(result)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/matrixutils.py", line 98, in to_sympy
    raise TypeError('Expected sympy/numpy/scipy.sparse matrix, got: %r' % m)
TypeError: Expected sympy/numpy/scipy.sparse matrix, got: None
```

or

```
>>> from sympy import *
>>> from sympy.physics.quantum import *
>>> from sympy.physics.quantum.gate import X,UGate
>>> v,w,f,g = symbols('v w f g')
>>> u2 = UGate(0,ImmutableMatrix([[v,w],[f,g]]))
>>> u3=X(0)+u2
>>> represent(u3,nqubits=2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>  
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/represent.py", line 165, in represent
    result = represent(expr.args[0], **options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/represent.py", line 145, in represent
    return expr._represent(**options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/qexpr.py", line 387, in _represent
    result = self._represent_default_basis(**options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/gate.py", line 256, in _represent_default_basis
    return self._represent_ZGate(None, **options)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/gate.py", line 278, in _represent_ZGate
    controls, targets, target_matrix, nqubits, format
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/gate.py", line 1094, in represent_zbasis
    return matrix_tensor_product(*product)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/matrixutils.py", line 247, in matrix_tensor_product
    return _sympy_tensor_product(*product)
  File "/Users/guest/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/sympy/physics/quantum/matrixutils.py", line 194, in _sympy_tensor_product
    'Sequence of Matrices expected, got: %s' % repr(matrices)
TypeError: Sequence of Matrices expected, got: (Matrix([
[1, 0],
[0, 1]]), Matrix([
[v, w],
[f, g]]))
```

According to some examination, I found an issue in physics/quantum/matrixutils.py.
Perhaps matrixutils.py had better to use DenseMatrix class instead of Matrix class. 
The test code would work well without any errors when I have replaced the 'Matrix' to 'DenseMatrix'.
I propose to fix physics/quantum/matrixutils.py.

Regards,

Discussion:

Patch:
diff --git a/sympy/physics/quantum/matrixutils.py b/sympy/physics/quantum/matrixutils.py
index 2bf954a029..4e1ea978de 100644
--- a/sympy/physics/quantum/matrixutils.py
+++ b/sympy/physics/quantum/matrixutils.py
@@ -2,7 +2,7 @@
 
 from __future__ import print_function, division
 
-from sympy import Matrix, I, Expr, Integer
+from sympy import MatrixBase, I, Expr, Integer
 from sympy.core.compatibility import range
 from sympy.matrices import eye, zeros
 from sympy.external import import_module
@@ -54,12 +54,12 @@ def sympy_to_numpy(m, **options):
     if not np:
         raise ImportError
     dtype = options.get('dtype', 'complex')
-    if isinstance(m, Matrix):
+    if isinstance(m, MatrixBase):
         return np.matrix(m.tolist(), dtype=dtype)
     elif isinstance(m, Expr):
         if m.is_Number or m.is_NumberSymbol or m == I:
             return complex(m)
-    raise TypeError('Expected Matrix or complex scalar, got: %r' % m)
+    raise TypeError('Expected MatrixBase or complex scalar, got: %r' % m)
 
 
 def sympy_to_scipy_sparse(m, **options):
@@ -67,27 +67,27 @@ def sympy_to_scipy_sparse(m, **options):
     if not np or not sparse:
         raise ImportError
     dtype = options.get('dtype', 'complex')
-    if isinstance(m, Matrix):
+    if isinstance(m, MatrixBase):
         return sparse.csr_matrix(np.matrix(m.tolist(), dtype=dtype))
     elif isinstance(m, Expr):
         if m.is_Number or m.is_NumberSymbol or m == I:
             return complex(m)
-    raise TypeError('Expected Matrix or complex scalar, got: %r' % m)
+    raise TypeError('Expected MatrixBase or complex scalar, got: %r' % m)
 
 
 def scipy_sparse_to_sympy(m, **options):
     """Convert a scipy.sparse matrix to a sympy matrix."""
-    return Matrix(m.todense())
+    return MatrixBase(m.todense())
 
 
 def numpy_to_sympy(m, **options):
     """Convert a numpy matrix to a sympy matrix."""
-    return Matrix(m)
+    return MatrixBase(m)
 
 
 def to_sympy(m, **options):
     """Convert a numpy/scipy.sparse matrix to a sympy matrix."""
-    if isinstance(m, Matrix):
+    if isinstance(m, MatrixBase):
         return m
     elif isinstance(m, numpy_ndarray):
         return numpy_to_sympy(m)
@@ -101,7 +101,7 @@ def to_sympy(m, **options):
 def to_numpy(m, **options):
     """Convert a sympy/scipy.sparse matrix to a numpy matrix."""
     dtype = options.get('dtype', 'complex')
-    if isinstance(m, (Matrix, Expr)):
+    if isinstance(m, (MatrixBase, Expr)):
         return sympy_to_numpy(m, dtype=dtype)
     elif isinstance(m, numpy_ndarray):
         return m
@@ -113,7 +113,7 @@ def to_numpy(m, **options):
 def to_scipy_sparse(m, **options):
     """Convert a sympy/numpy matrix to a scipy.sparse matrix."""
     dtype = options.get('dtype', 'complex')
-    if isinstance(m, (Matrix, Expr)):
+    if isinstance(m, (MatrixBase, Expr)):
         return sympy_to_scipy_sparse(m, dtype=dtype)
     elif isinstance(m, numpy_ndarray):
         if not sparse:
@@ -126,7 +126,7 @@ def to_scipy_sparse(m, **options):
 
 def flatten_scalar(e):
     """Flatten a 1x1 matrix to a scalar, return larger matrices unchanged."""
-    if isinstance(e, Matrix):
+    if isinstance(e, MatrixBase):
         if e.shape == (1, 1):
             e = e[0]
     if isinstance(e, (numpy_ndarray, scipy_sparse_matrix)):
@@ -137,7 +137,7 @@ def flatten_scalar(e):
 
 def matrix_dagger(e):
     """Return the dagger of a sympy/numpy/scipy.sparse matrix."""
-    if isinstance(e, Matrix):
+    if isinstance(e, MatrixBase):
         return e.H
     elif isinstance(e, (numpy_ndarray, scipy_sparse_matrix)):
         return e.conjugate().transpose()
@@ -153,13 +153,13 @@ def _sympy_tensor_product(*matrices):
     Parameters
     ==========
 
-    matrices : tuple of Matrix instances
+    matrices : tuple of MatrixBase instances
         The matrices to take the tensor product of.
 
     Returns
     =======
 
-    matrix : Matrix
+    matrix : MatrixBase
         The tensor product matrix.
 
     Examples
@@ -189,7 +189,7 @@ def _sympy_tensor_product(*matrices):
     [1] http://en.wikipedia.org/wiki/Kronecker_product
     """
     # Make sure we have a sequence of Matrices
-    if not all(isinstance(m, Matrix) for m in matrices):
+    if not all(isinstance(m, MatrixBase) for m in matrices):
         raise TypeError(
             'Sequence of Matrices expected, got: %s' % repr(matrices)
         )
@@ -243,7 +243,7 @@ def _scipy_sparse_tensor_product(*product):
 
 def matrix_tensor_product(*product):
     """Compute the matrix tensor product of sympy/numpy/scipy.sparse matrices."""
-    if isinstance(product[0], Matrix):
+    if isinstance(product[0], MatrixBase):
         return _sympy_tensor_product(*product)
     elif isinstance(product[0], numpy_ndarray):
         return _numpy_tensor_product(*product)
@@ -336,7 +336,7 @@ def _scipy_sparse_matrix_to_zero(e):
 
 def matrix_to_zero(e):
     """Convert a zero matrix to the scalar zero."""
-    if isinstance(e, Matrix):
+    if isinstance(e, MatrixBase):
         if zeros(*e.shape) == e:
             e = Integer(0)
     elif isinstance(e, numpy_ndarray):


Test patch:
diff --git a/sympy/physics/quantum/tests/test_gate.py b/sympy/physics/quantum/tests/test_gate.py
index 92228e1db1..82acb7c550 100644
--- a/sympy/physics/quantum/tests/test_gate.py
+++ b/sympy/physics/quantum/tests/test_gate.py
@@ -1,6 +1,6 @@
 from sympy import exp, symbols, sqrt, I, pi, Mul, Integer, Wild
 from sympy.core.compatibility import range
-from sympy.matrices import Matrix
+from sympy.matrices import Matrix, ImmutableMatrix
 
 from sympy.physics.quantum.gate import (XGate, YGate, ZGate, random_circuit,
         CNOT, IdentityGate, H, X, Y, S, T, Z, SwapGate, gate_simp, gate_sort,
@@ -133,6 +133,19 @@ def test_UGate_CGate_combo():
         assert u2Rep*qubit_to_matrix(IntQubit(i, 2)) == \
             qubit_to_matrix(qapply(u2*IntQubit(i, 2)))
 
+def test_UGate_OneQubitGate_combo():
+    v, w, f, g = symbols('v w f g')
+    uMat1 = ImmutableMatrix([[v, w], [f, g]])
+    cMat1 = Matrix([[v, w + 1, 0, 0], [f + 1, g, 0, 0], [0, 0, v, w + 1], [0, 0, f + 1, g]])
+    u1 = X(0) + UGate(0, uMat1)
+    assert represent(u1, nqubits=2) == cMat1
+
+    uMat2 = ImmutableMatrix([[1/sqrt(2), 1/sqrt(2)], [I/sqrt(2), -I/sqrt(2)]])
+    cMat2_1 = Matrix([[1/2 + I/2, 1/2 - I/2], [1/2 - I/2, 1/2 + I/2]])
+    cMat2_2 = Matrix([[1, 0], [0, I]])
+    u2 = UGate(0, uMat2)
+    assert represent(H(0)*u2, nqubits=1) == cMat2_1
+    assert represent(u2*H(0), nqubits=1) == cMat2_2
 
 def test_represent_hadamard():
     """Test the representation of the hadamard gate."""