Merge branch 'master' into initialshape

pyproject.toml

@@ -84,6 +84,7 @@ docs-clean = "make -C doc clean"
 test-docs = "make -C doc doctest"
 
 
+
 [tool.poetry.group.cuda11]
 optional = true
 
@@ -93,6 +94,7 @@ cuquantum-python-cu11 = "^23.3.0"
 qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 
+
 [tool.poetry.group.cuda12]
 optional = true
 

src/qibo/models/dbi/double_bracket.py

@@ -105,29 +105,60 @@ def __init__(
     def __call__(
         self, step: float, mode: DoubleBracketGeneratorType = None, d: np.array = None
     ):
-        """Performs one double bracket rotation."""
+        r"""We use convention that $H' = U^\dagger H U$ where $U=e^{-sW}$ with $W=[D,H]$
+        (or depending on `mode` an approximation, see `eval_dbr_unitary`).
+        If $s>0$ then for $D = \Delta(H)$ the GWW DBR will give a $\sigma$-decrease,
+        see https://arxiv.org/abs/2206.11772."""
+
+        operator = self.eval_dbr_unitary(step, mode, d)
+        operator_dagger = self.backend.cast(
+            np.array(np.matrix(self.backend.to_numpy(operator)).getH())
+        )
+        self.h.matrix = operator_dagger @ self.h.matrix @ operator
+        return operator
+
+    def eval_dbr_unitary(
+        self,
+        step: float,
+        mode: DoubleBracketGeneratorType = None,
+        d: np.array = None,
+    ):
+        """In call we are working in the convention that $H' = U^\\dagger H
+        U$ where $U=e^{-sW}$ with $W=[D,H]$ or an approximation of that by a group commutator.
+        That is handy because if we switch from the DBI in the Heisenberg picture for the
+        Hamiltonian, we get that the transformation of the state is $|\\psi'\rangle = U |\\psi\rangle$
+        so that $\\langle H\rangle_{\\psi'} = \\langle H' \rangle_\\psi$ (i.e. when writing the unitary
+        acting on the state dagger notation is avoided).
+        The group commutator must approximate $U=e^{-s[D,H]}$. This is achieved by setting $r = \\sqrt{s}$ so that
+        $$V = e^{-irH}e^{irD}e^{irH}e^{-irD}$$
+        because
+        $$e^{-irH}De^{irH} = D+ir[D,H]+O(r^2)$$
+        so
+        $$V\approx e^{irD +i^2 r^2[D,H] + O(r^2) -irD} \approx U\\ .$$
+        See the app in https://arxiv.org/abs/2206.11772 for a derivation.
+        """
         if mode is None:
             mode = self.mode
 
         if mode is DoubleBracketGeneratorType.canonical:
             operator = self.backend.calculate_matrix_exp(
-                1.0j * step,
+                -1.0j * step,
                 self.commutator(self.diagonal_h_matrix, self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.single_commutator:
             if d is None:
                 d = self.diagonal_h_matrix
             operator = self.backend.calculate_matrix_exp(
-                1.0j * step,
+                -1.0j * step,
                 self.commutator(self.backend.cast(d), self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.group_commutator:
             if d is None:
                 d = self.diagonal_h_matrix
             operator = (
-                self.h.exp(-step)
+                self.h.exp(step)
                 @ self.backend.calculate_matrix_exp(-step, d)
-                @ self.h.exp(step)
+                @ self.h.exp(-step)
                 @ self.backend.calculate_matrix_exp(step, d)
             )
         elif mode is DoubleBracketGeneratorType.group_commutator_third_order:
@@ -141,11 +172,20 @@ def __call__(
                 @ self.h.exp(-step * (3 - np.sqrt(5)) / 2)
                 @ self.backend.calculate_matrix_exp(-step, d)
             )
-        operator_dagger = self.backend.cast(
-            np.array(np.matrix(self.backend.to_numpy(operator)).getH())
-        )
+            operator = (
+                self.backend.calculate_matrix_exp(step, d)
+                @ self.h.exp(step * (3 - np.sqrt(5)) / 2)
+                @ self.backend.calculate_matrix_exp(-step * (np.sqrt(5) + 1) / 2, d)
+                @ self.h.exp(-step)
+                @ self.backend.calculate_matrix_exp(step * (np.sqrt(5) - 1) / 2, d)
+                @ self.h.exp(step * (np.sqrt(5) - 1) / 2)
+            )
+        else:
+            raise NotImplementedError(
+                f"The mode {mode} is not supported"
+            )  # pragma: no cover
 
-        self.h.matrix = operator @ self.h.matrix @ operator_dagger
+        return operator
 
     @staticmethod
     def commutator(a, b):

src/qibo/models/variational.py

@@ -96,7 +96,6 @@ def minimize(
             loss = lambda p, c, h: dtype(loss_func(p, c, h))
         elif method != "sgd":
             loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(loss_func(p, c, h))
-
         result, parameters, extra = self.optimizers.optimize(
             loss,
             initial_state,

src/qibo/quantum_info/entropies.py

@@ -461,7 +461,7 @@ def von_neumann_entropy(
 
     if purity(state) == 1.0:
         if return_spectrum:
-            return 0.0, backend.cast([1.0], dtype=float)
+            return 0.0, backend.cast([0.0], dtype=float)
 
         return 0.0
 

src/qibo/transpiler/decompositions.py

@@ -466,3 +466,23 @@ def _u3_to_gpi2(t, p, l):
     gates.ECR, [gates.S(0), gates.SX(1), gates.CNOT(0, 1), gates.X(0)]
 )
 standard_decompositions.add(gates.CCZ, [gates.H(2), gates.TOFFOLI(0, 1, 2), gates.H(2)])
+standard_decompositions.add(
+    gates.TOFFOLI,
+    [
+        gates.H(2),
+        gates.CNOT(1, 2),
+        gates.TDG(2),
+        gates.CNOT(0, 2),
+        gates.T(2),
+        gates.CNOT(1, 2),
+        gates.T(1),
+        gates.TDG(2),
+        gates.CNOT(0, 2),
+        gates.CNOT(0, 1),
+        gates.T(2),
+        gates.T(0),
+        gates.TDG(1),
+        gates.H(2),
+        gates.CNOT(0, 1),
+    ],
+)

src/qibo/transpiler/pipeline.py

@@ -225,7 +225,7 @@ def __call__(self, circuit):
         physical (keys) to logical (values) qubit. If `int_qubit_name` is `True`
         each key `i` correspond to the `i-th` qubit in the graph.
         """
-        final_layout = self.initial_layout
+        self.initial_layout = None
         for transpiler_pass in self.passes:
             if isinstance(transpiler_pass, Optimizer):
                 transpiler_pass.connectivity = self.connectivity
@@ -234,6 +234,9 @@ def __call__(self, circuit):
                 transpiler_pass.connectivity = self.connectivity
                 if self.initial_layout == None:
                     self.initial_layout = transpiler_pass(circuit)
+                    final_layout = (
+                        self.initial_layout
+                    )  # This way the final layout will be the same as the initial layout if no router is used
                 else:
                     raise_error(
                         TranspilerPipelineError,

tests/test_callbacks.py

@@ -97,7 +97,7 @@ def test_entropy_in_circuit(backend, density_matrix, base):
     values = [backend.to_numpy(x) for x in entropy]
     backend.assert_allclose(values, target, atol=PRECISION_TOL)
 
-    target_spectrum = [1.0] + list([0, 0, np.log(2), np.log(2)] / np.log(base))
+    target_spectrum = [0.0] + list([0, 0, np.log(2), np.log(2)] / np.log(base))
     entropy_spectrum = np.ravel(np.concatenate(entropy.spectrum)).tolist()
     backend.assert_allclose(entropy_spectrum, target_spectrum, atol=PRECISION_TOL)
 

tests/test_cirq.py

@@ -4,8 +4,9 @@
 import numpy as np
 import pytest
 
-from qibo import gates, models
+from qibo import Circuit, gates, matrices
 from qibo.backends import NumpyBackend
+from qibo.models import QFT
 from qibo.quantum_info import random_statevector, random_unitary
 
 numpy_backend = NumpyBackend()
@@ -54,7 +55,7 @@ def assert_gates_equivalent(
     if ndevices is not None:
         accelerators = {"/GPU:0": ndevices}
 
-    c = models.Circuit(nqubits, accelerators)
+    c = Circuit(nqubits, accelerators)
     c.add(qibo_gate)
     assert c.depth == target_depth
     if accelerators and not backend.supports_multigpu:
@@ -307,7 +308,7 @@ def test_unitary_matrix_gate_controlled_by(backend, nqubits, ntargets, ndevices)
 
 @pytest.mark.parametrize("nqubits", [5, 6, 7, 11, 12])
 def test_qft(backend, accelerators, nqubits):
-    c = models.QFT(nqubits, accelerators=accelerators)
+    c = QFT(nqubits, accelerators=accelerators)
     initial_state = random_statevector(2**nqubits, backend=numpy_backend)
     final_state = backend.execute_circuit(c, np.copy(initial_state)).state()
     final_state = backend.cast(final_state, dtype=final_state.dtype)

tests/test_gates_gates.py

@@ -6,6 +6,7 @@
 from qibo import Circuit, gates, matrices
 from qibo.parameter import Parameter
 from qibo.quantum_info import random_hermitian, random_statevector, random_unitary
+from qibo.transpiler.decompositions import standard_decompositions
 
 
 def apply_gates(backend, gatelist, nqubits=None, initial_state=None):
@@ -1205,6 +1206,13 @@ def test_toffoli(backend, applyx):
     assert not gates.TOFFOLI(0, 1, 2).clifford
     assert gates.TOFFOLI(0, 1, 2).unitary
 
+    # test decomposition
+    decomposition = Circuit(3)
+    decomposition.add(standard_decompositions(gates.TOFFOLI(0, 1, 2)))
+    decomposition = decomposition.unitary(backend)
+
+    backend.assert_allclose(decomposition, backend.cast(matrices.TOFFOLI), atol=1e-10)
+
 
 def test_ccz(backend):
     nqubits = 3