Repository: qiboteam/qibo

Issue:
trotter evolution and dense evolution does not match
We should check the trotterization evolution code, we are observing non convergence with this approach, while with dense matrices the evolution works well.

Discussion:

Patch:
diff --git a/src/qibo/core/hamiltonians.py b/src/qibo/core/hamiltonians.py
index 68769310f..f3a4056ec 100644
--- a/src/qibo/core/hamiltonians.py
+++ b/src/qibo/core/hamiltonians.py
@@ -254,7 +254,7 @@ class SymbolicHamiltonian(hamiltonians.SymbolicHamiltonian):
         self.form = sympy.expand(form)
         terms = []
         for f, c in self.form.as_coefficients_dict().items():
-            term = SymbolicTerm.from_factors(c, f, symbol_map)
+            term = SymbolicTerm(c, f, symbol_map)
             if term.target_qubits:
                 terms.append(term)
             else:
diff --git a/src/qibo/core/terms.py b/src/qibo/core/terms.py
index bf5d90168..32306801b 100644
--- a/src/qibo/core/terms.py
+++ b/src/qibo/core/terms.py
@@ -4,10 +4,12 @@ from qibo.config import raise_error
 
 
 class HamiltonianTerm:
-    """Single term in a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+    """Term of a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
 
-    This Hamiltonian is represented by a sum of such terms which are held in
-    its ``.terms`` attribute which is a list.
+    Symbolic Hamiltonians are represented by a list of
+    :class:`qibo.core.terms.HamiltonianTerm` objects storred in the
+    ``SymbolicHamiltonian.terms`` attribute. The mathematical expression of
+    the Hamiltonian is the sum of these terms.
 
     Args:
         matrix (np.ndarray): Full matrix corresponding to the term representation
@@ -38,11 +40,12 @@ class HamiltonianTerm:
 
     @property
     def matrix(self):
+        """Matrix representation of the term."""
         return self._matrix
 
     @property
     def gate(self):
-        """Qibo gate that implements the action of the term on states."""
+        """:class:`qibo.abstractions.gates.Unitary` gate that implements the action of the term on states."""
         if self._gate is None:
             self._gate = gates.Unitary(self.matrix, *self.target_qubits)
         return self._gate
@@ -52,12 +55,13 @@ class HamiltonianTerm:
         return K.qnp.expm(-1j * x * self.matrix)
 
     def expgate(self, x):
-        """Unitary gate implementing the matrix exponentiation of the term."""
+        """:class:`qibo.abstractions.gates.Unitary` gate implementing the action of exp(term) on states."""
         return gates.Unitary(self.exp(x), *self.target_qubits)
 
     def merge(self, term):
         """Creates a new term by merging the given term to the current one.
 
+        The resulting term corresponds to the sum of the two original terms.
         The target qubits of the given term should be a subset of the target
         qubits of the current term.
         """
@@ -98,7 +102,7 @@ class HamiltonianTerm:
 
 
 class SymbolicTerm(HamiltonianTerm):
-    """:class:`qibo.core.terms.HamiltonianTerm` constructed from ``sympy`` symbols.
+    """:class:`qibo.core.terms.HamiltonianTerm` constructed using ``sympy`` expression.
 
     Example:
         ::
@@ -113,59 +117,64 @@ class SymbolicTerm(HamiltonianTerm):
         coefficient (complex): Complex number coefficient of the underlying
             term in the Hamiltonian.
         factors (sympy.Expr): Sympy expression for the underlying term.
-        matrix_map (dict): Dictionary that maps symbols in the given ``factors``
+        symbol_map (dict): Dictionary that maps symbols in the given ``factors``
             expression to tuples of (target qubit id, matrix).
             This is required only if the expression is not created using Qibo
             symbols and to keep compatibility with older versions where Qibo
             symbols were not available.
     """
 
-    def __init__(self, coefficient, factors=[], matrix_map={}):
+    def __init__(self, coefficient, factors=1, symbol_map=None):
         self.coefficient = complex(coefficient)
-        self.factors = factors
-        self.matrix_map = matrix_map
         self._matrix = None
         self._gate = None
         self.hamiltonian = None
-        self.target_qubits = tuple(sorted(self.matrix_map.keys()))
 
-    @classmethod
-    def from_factors(cls, coefficient, factors, symbol_map=None):
-        if factors == 1:
-            return cls(coefficient)
-
-        _factors = []
-        _matrix_map = {}
-        for factor in factors.as_ordered_factors():
-            if isinstance(factor, sympy.Pow):
-                factor, pow = factor.args
-                assert isinstance(pow, sympy.Integer)
-                assert isinstance(factor, sympy.Symbol)
-                pow = int(pow)
-            else:
-                pow = 1
-
-            if symbol_map is not None and factor in symbol_map:
-                from qibo.symbols import Symbol
-                q, matrix = symbol_map.get(factor)
-                factor = Symbol(q, matrix, name=factor.name)
-
-            if isinstance(factor, sympy.Symbol):
-                if isinstance(factor.matrix, K.qnp.tensor_types):
-                    _factors.extend(pow * [factor])
-                    q = factor.target_qubit
-                    if q in _matrix_map:
-                        _matrix_map[q].extend(pow * [factor.matrix])
-                    else:
-                        _matrix_map[q] = pow * [factor.matrix]
+        # List of :class:`qibo.symbols.Symbol` that represent the term factors
+        self.factors = []
+        # Dictionary that maps target qubit ids to a list of matrices that act on each qubit
+        self.matrix_map = {}
+        if factors != 1:
+            for factor in factors.as_ordered_factors():
+                # check if factor has some power ``power`` so that the corresponding
+                # matrix is multiplied ``pow`` times
+                if isinstance(factor, sympy.Pow):
+                    factor, pow = factor.args
+                    assert isinstance(pow, sympy.Integer)
+                    assert isinstance(factor, sympy.Symbol)
+                    pow = int(pow)
                 else:
-                    coefficient *= factor.matrix
-            elif factor == sympy.I:
-                coefficient *= 1j
-            else: # pragma: no cover
-                raise_error(TypeError, "Cannot parse factor {}.".format(factor))
+                    pow = 1
+
+                # if the user is using ``symbol_map`` instead of qibo symbols,
+                # create the corresponding symbols
+                if symbol_map is not None and factor in symbol_map:
+                    from qibo.symbols import Symbol
+                    q, matrix = symbol_map.get(factor)
+                    factor = Symbol(q, matrix, name=factor.name)
+
+                if isinstance(factor, sympy.Symbol):
+                    if isinstance(factor.matrix, K.qnp.tensor_types):
+                        self.factors.extend(pow * [factor])
+                        q = factor.target_qubit
+                        # if pow > 1 the matrix should be multiplied multiple
+                        # when calculating the term's total matrix so we
+                        # repeat it in the corresponding list that will
+                        # be used during this calculation
+                        # see the ``SymbolicTerm.matrix`` property for the
+                        # full matrix calculation
+                        if q in self.matrix_map:
+                            self.matrix_map[q].extend(pow * [factor.matrix])
+                        else:
+                            self.matrix_map[q] = pow * [factor.matrix]
+                    else:
+                        self.coefficient *= factor.matrix
+                elif factor == sympy.I:
+                    self.coefficient *= 1j
+                else: # pragma: no cover
+                    raise_error(TypeError, "Cannot parse factor {}.".format(factor))
 
-        return cls(coefficient, _factors, _matrix_map)
+        self.target_qubits = tuple(sorted(self.matrix_map.keys()))
 
     @property
     def matrix(self):
@@ -178,6 +187,12 @@ class SymbolicTerm(HamiltonianTerm):
         """
         if self._matrix is None:
             def matrices_product(matrices):
+                """Product of matrices that act on the same tuple of qubits.
+
+                Args:
+                    matrices (list): List of matrices to multiply, as exists in
+                        the values of ``SymbolicTerm.matrix_map``.
+                """
                 if len(matrices) == 1:
                     return matrices[0]
                 matrix = K.np.copy(matrices[0])
@@ -191,11 +206,20 @@ class SymbolicTerm(HamiltonianTerm):
                 self._matrix = K.np.kron(self._matrix, matrix)
         return self._matrix
 
+    def copy(self):
+        """Creates a shallow copy of the term with the same attributes."""
+        new = self.__class__(self.coefficient)
+        new.factors = self.factors
+        new.matrix_map = self.matrix_map
+        new.target_qubits = self.target_qubits
+        return new
+
     def __mul__(self, x):
-        new = self.__class__(self.coefficient, self.factors, self.matrix_map)
-        new._matrix = self._matrix
-        new._gate = self._gate
+        """Multiplication of scalar to the Hamiltonian term."""
+        new = self.copy()
         new.coefficient *= x
+        if self._matrix is not None:
+            new._matrix = x * self._matrix
         return new
 
     def __call__(self, state, density_matrix=False):
@@ -211,7 +235,8 @@ class SymbolicTerm(HamiltonianTerm):
 class TermGroup(list):
     """Collection of multiple :class:`qibo.core.terms.HamiltonianTerm` objects.
 
-    Allows merging multiple terms to a single one for faster exponentiation.
+    Allows merging multiple terms to a single one for faster exponentiation
+    during Trotterized evolution.
 
     Args:
         term (:class:`qibo.core.terms.HamiltonianTerm`): Parent term of the group.
@@ -236,7 +261,9 @@ class TermGroup(list):
 
     @classmethod
     def from_terms(cls, terms):
-        """Groups a list of terms to multiple groups.
+        """Divides a list of terms to multiple :class:`qibo.core.terms.TermGroup`s.
+
+        Terms that target the same qubits are grouped to the same group.
 
         Args:
             terms (list): List of :class:`qibo.core.terms.HamiltonianTerm` objects.
@@ -245,6 +272,8 @@ class TermGroup(list):
             List of :class:`qibo.core.terms.TermGroup` objects that contain
             all the given terms.
         """
+        # split given terms according to their order
+        # order = len(term.target_qubits)
         orders = {}
         for term in terms:
             if len(term) in orders:
@@ -253,6 +282,8 @@ class TermGroup(list):
                 orders[len(term)] = [term]
 
         groups = []
+        # start creating groups with the higher order terms as parents and then
+        # append each term of lower order to the first compatible group
         for order in sorted(orders.keys())[::-1]:
             for child in orders[order]:
                 flag = True


Test patch:
diff --git a/src/qibo/tests/test_core_terms.py b/src/qibo/tests/test_core_terms.py
new file mode 100644
index 000000000..fbc9e677a
--- /dev/null
+++ b/src/qibo/tests/test_core_terms.py
@@ -0,0 +1,178 @@
+"""Tests methods defined in `qibo/core/terms.py`."""
+import pytest
+import numpy as np
+from qibo import K, matrices, models, gates
+from qibo.core import terms
+from qibo.tests.utils import random_state, random_density_matrix
+
+
+def test_hamiltonian_term_initialization(backend):
+    """Test initialization and matrix assignment of ``HamiltonianTerm``."""
+    matrix = np.random.random((2, 2))
+    term = terms.HamiltonianTerm(matrix, 0)
+    assert term.target_qubits == (0,)
+    K.assert_allclose(term.matrix, matrix)
+    matrix = np.random.random((4, 4))
+    term = terms.HamiltonianTerm(matrix, 2, 3)
+    assert term.target_qubits == (2, 3)
+    K.assert_allclose(term.matrix, matrix)
+
+
+def test_hamiltonian_term_gates(backend):
+    """Test gate application of ``HamiltonianTerm``."""
+    matrix = np.random.random((4, 4))
+    term = terms.HamiltonianTerm(matrix, 1, 2)
+    gate = term.gate
+    assert gate.target_qubits == (1, 2)
+    K.assert_allclose(gate.matrix, matrix)
+
+    initial_state = random_state(4)
+    final_state = term(np.copy(initial_state))
+    c = models.Circuit(4)
+    c.add(gates.Unitary(matrix, 1, 2))
+    target_state = c(np.copy(initial_state))
+    K.assert_allclose(final_state, target_state)
+
+
+def test_hamiltonian_term_exponentiation(backend):
+    """Test exp gate application of ``HamiltonianTerm``."""
+    from scipy.linalg import expm
+    matrix = np.random.random((2, 2))
+    term = terms.HamiltonianTerm(matrix, 1)
+    exp_matrix = expm(-0.5j * matrix)
+    K.assert_allclose(term.exp(0.5), exp_matrix)
+
+    initial_state = random_state(2)
+    final_state = term.expgate(0.5)(np.copy(initial_state))
+    exp_gate = gates.Unitary(exp_matrix, 1)
+    target_state = exp_gate(np.copy(initial_state))
+    K.assert_allclose(final_state, target_state)
+
+
+def test_hamiltonian_term_mul(backend):
+    """Test scalar multiplication of ``HamiltonianTerm``."""
+    matrix = np.random.random((4, 4))
+    term = terms.HamiltonianTerm(matrix, 0, 2)
+    mterm = 2 * term
+    assert mterm.target_qubits == term.target_qubits
+    K.assert_allclose(mterm.matrix, 2 * matrix)
+    mterm = term * 5
+    assert mterm.target_qubits == term.target_qubits
+    K.assert_allclose(mterm.matrix, 5 * matrix)
+
+
+def test_hamiltonian_term_merge(backend):
+    """Test ``HamiltonianTerm.merge``."""
+    matrix1 = np.random.random((2, 2))
+    matrix2 = np.random.random((4, 4))
+    term1 = terms.HamiltonianTerm(matrix1, 1)
+    term2 = terms.HamiltonianTerm(matrix2, 0, 1)
+    mterm = term2.merge(term1)
+    target_matrix = np.kron(np.eye(2), matrix1) + matrix2
+    assert mterm.target_qubits == (0, 1)
+    K.assert_allclose(mterm.matrix, target_matrix)
+    with pytest.raises(ValueError):
+        term1.merge(term2)
+
+
+@pytest.mark.parametrize("use_symbols", [True, False])
+def test_symbolic_term_creation(use_symbols):
+    """Test creating ``SymbolicTerm`` from sympy expression."""
+    if use_symbols:
+        from qibo.symbols import X, Y
+        expression = X(0) * Y(1) * X(1)
+        symbol_map = None
+    else:
+        import sympy
+        x0, x1, y1 = sympy.symbols("X0 X1 Y1", commutative=False)
+        expression = x0 * y1 * x1
+        symbol_map = {x0: (0, matrices.X), x1: (1, matrices.X),
+                      y1: (1, matrices.Y)}
+    term = terms.SymbolicTerm(2, expression, symbol_map)
+    assert term.target_qubits == (0, 1)
+    assert len(term.matrix_map) == 2
+    K.assert_allclose(term.matrix_map.get(0)[0], matrices.X)
+    K.assert_allclose(term.matrix_map.get(1)[0], matrices.Y)
+    K.assert_allclose(term.matrix_map.get(1)[1], matrices.X)
+
+
+def test_symbolic_term_matrix(backend):
+    """Test matrix calculation of ``SymbolicTerm``."""
+    from qibo.symbols import X, Y, Z
+    expression = X(0) * Y(1) * Z(2) * X(1)
+    term = terms.SymbolicTerm(2, expression)
+    assert term.target_qubits == (0, 1, 2)
+    target_matrix = np.kron(matrices.X, matrices.Y @ matrices.X)
+    target_matrix = 2 * np.kron(target_matrix, matrices.Z)
+    K.assert_allclose(term.matrix, target_matrix)
+
+
+def test_symbolic_term_mul(backend):
+    """Test multiplying scalar to ``SymbolicTerm``."""
+    from qibo.symbols import X, Y, Z
+    expression = Y(2) * Z(3) * X(2) * X(3)
+    term = terms.SymbolicTerm(1, expression)
+    assert term.target_qubits == (2, 3)
+    target_matrix = np.kron(matrices.Y @ matrices.X, matrices.Z @ matrices.X)
+    K.assert_allclose(term.matrix, target_matrix)
+    K.assert_allclose((2 * term).matrix, 2 * target_matrix)
+    K.assert_allclose((term * 3).matrix, 3 * target_matrix)
+
+
+@pytest.mark.parametrize("density_matrix", [False])
+def test_symbolic_term_call(backend, density_matrix):
+    """Test applying ``SymbolicTerm`` to state."""
+    from qibo.symbols import X, Y, Z
+    expression = Z(0) * X(1) * Y(2)
+    term = terms.SymbolicTerm(2, expression)
+    matrixlist = [np.kron(matrices.Z, np.eye(4)),
+                  np.kron(np.kron(np.eye(2), matrices.X), np.eye(2)),
+                  np.kron(np.eye(4), matrices.Y)]
+    initial_state = random_density_matrix(3) if density_matrix else random_state(3)
+    final_state = term(np.copy(initial_state), density_matrix=density_matrix)
+    target_state = 2 * np.copy(initial_state)
+    for matrix in matrixlist:
+        target_state = matrix @ target_state
+    K.assert_allclose(final_state, target_state)
+
+
+def test_symbolic_term_merge(backend):
+    """Test merging ``SymbolicTerm`` to ``HamiltonianTerm``."""
+    from qibo.symbols import X, Z
+    matrix = np.random.random((4, 4))
+    term1 = terms.HamiltonianTerm(matrix, 0, 1)
+    term2 = terms.SymbolicTerm(1, X(0) * Z(1))
+    term = term1.merge(term2)
+    target_matrix = matrix + np.kron(matrices.X, matrices.Z)
+    K.assert_allclose(term.matrix, target_matrix)
+
+
+def test_term_group_append():
+    """Test ``GroupTerm.can_append`` method."""
+    term1 = terms.HamiltonianTerm(np.random.random((8, 8)), 0, 1, 3)
+    term2 = terms.HamiltonianTerm(np.random.random((2, 2)), 0)
+    term3 = terms.HamiltonianTerm(np.random.random((2, 2)), 1)
+    term4 = terms.HamiltonianTerm(np.random.random((2, 2)), 2)
+    group = terms.TermGroup(term1)
+    assert group.can_append(term2)
+    assert group.can_append(term3)
+    assert not group.can_append(term4)
+    group.append(term2)
+    group.append(term3)
+    assert group.target_qubits == {0, 1, 3}
+
+
+def test_term_group_to_term(backend):
+    """Test ``GroupTerm.term`` property."""
+    from qibo.symbols import X, Z
+    matrix = np.random.random((8, 8))
+    term1 = terms.HamiltonianTerm(matrix, 0, 1, 3)
+    term2 = terms.SymbolicTerm(1, X(0) * Z(3))
+    term3 = terms.SymbolicTerm(2, X(1))
+    group = terms.TermGroup(term1)
+    group.append(term2)
+    group.append(term3)
+    matrix2 = np.kron(np.kron(matrices.X, np.eye(2)), matrices.Z)
+    matrix3 = np.kron(np.kron(np.eye(2), matrices.X), np.eye(2))
+    target_matrix = matrix + matrix2 + 2 * matrix3
+    K.assert_allclose(group.term.matrix, target_matrix)