From dc195c1aa19edd5fc1c4c18101bf19841f13983e Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 20 Jun 2024 15:44:15 +0400
Subject: [PATCH 01/16] fix `controlled_by` method

---
 src/qibo/backends/numpy.py | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index f280a05640..da705ff67e 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 from scipy import sparse
+from scipy.linalg import block_diag
 
 from qibo import __version__
 from qibo.backends import einsum_utils
@@ -110,13 +111,25 @@ def matrix(self, gate):
         if callable(_matrix):
             _matrix = _matrix(2 ** len(gate.target_qubits))
 
+        num_controls = len(gate.control_qubits)
+        if num_controls > 0:
+            _matrix = block_diag(
+                np.eye(2 ** len(gate.qubits) - len(_matrix)), self.to_numpy(_matrix)
+            )
+
         return self.cast(_matrix, dtype=_matrix.dtype)
 
     def matrix_parametrized(self, gate):
         """Convert a parametrized gate to its matrix representation in the computational basis."""
         name = gate.__class__.__name__
-        matrix = getattr(self.matrices, name)(*gate.parameters)
-        return self.cast(matrix, dtype=matrix.dtype)
+        _matrix = getattr(self.matrices, name)(*gate.parameters)
+
+        num_controls = len(gate.control_qubits)
+        if num_controls > 0:
+            _matrix = block_diag(
+                np.eye(2 ** len(gate.qubits) - 2), self.to_numpy(_matrix)
+            )
+        return self.cast(_matrix, dtype=_matrix.dtype)
 
     def matrix_fused(self, fgate):
         rank = len(fgate.target_qubits)

From 79b5bcfecb849595628006c9e5754ef088d671f1 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 20 Jun 2024 15:51:59 +0400
Subject: [PATCH 02/16] fix dimensions

---
 src/qibo/backends/numpy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index da705ff67e..d05c10bd76 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -127,7 +127,7 @@ def matrix_parametrized(self, gate):
         num_controls = len(gate.control_qubits)
         if num_controls > 0:
             _matrix = block_diag(
-                np.eye(2 ** len(gate.qubits) - 2), self.to_numpy(_matrix)
+                np.eye(2 ** len(gate.qubits) - len(_matrix)), self.to_numpy(_matrix)
             )
         return self.cast(_matrix, dtype=_matrix.dtype)
 

From 6559920195b254591c9c08a22e72f661387dab82 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 20 Jun 2024 16:34:48 +0400
Subject: [PATCH 03/16] fix `apply_gate`

---
 src/qibo/backends/numpy.py | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index d05c10bd76..fa4555be2f 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -186,9 +186,17 @@ def control_matrix(self, gate):
         return self.np.concatenate([part1, part2], axis=1)
 
     def apply_gate(self, gate, state, nqubits):
+        from qibo.gates.abstract import ParametrizedGate
+
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
-        matrix = gate.matrix(self)
+        # matrix = gate.matrix(self)
+        matrix = (
+            gate.__class__(*gate.qubits, *gate.parameters)
+            if isinstance(gate, ParametrizedGate)
+            else gate.__class__(*gate.qubits)
+        )
+        matrix = matrix.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
             ncontrol = len(gate.control_qubits)
@@ -213,9 +221,17 @@ def apply_gate(self, gate, state, nqubits):
         return self.np.reshape(state, (2**nqubits,))
 
     def apply_gate_density_matrix(self, gate, state, nqubits):
+        from qibo.gates.abstract import ParametrizedGate
+
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
-        matrix = gate.matrix(self)
+        # matrix = gate.matrix(self)
+        matrix = (
+            gate.__class__(*gate.qubits, *gate.parameters)
+            if isinstance(gate, ParametrizedGate)
+            else gate.__class__(*gate.qubits)
+        )
+        matrix = matrix.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
             matrixc = self.np.conj(matrix)

From 2c8ff415f94b2b1448ce4f8da720dc6493083396 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Fri, 21 Jun 2024 11:02:20 +0400
Subject: [PATCH 04/16] fix apply gate and remove `control_matrix` method

---
 src/qibo/backends/abstract.py |  5 ----
 src/qibo/backends/numpy.py    | 54 ++++++++++++-----------------------
 tests/test_backends.py        | 32 ---------------------
 tests/test_gates_gates.py     | 17 -----------
 4 files changed, 18 insertions(+), 90 deletions(-)

diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index 4a725c22e0..8e4175b8e2 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -133,11 +133,6 @@ def matrix_fused(self, gate):  # pragma: no cover
         """Fuse matrices of multiple gates."""
         raise_error(NotImplementedError)
 
-    @abc.abstractmethod
-    def control_matrix(self, gate):  # pragma: no cover
-        """ "Calculate full matrix representation of a controlled gate."""
-        raise_error(NotImplementedError)
-
     @abc.abstractmethod
     def apply_gate(self, gate, state, nqubits):  # pragma: no cover
         """Apply a gate to state vector."""
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index fa4555be2f..67a647b0b9 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -161,41 +161,20 @@ def matrix_fused(self, fgate):
 
         return self.cast(matrix.toarray())
 
-    def control_matrix(self, gate):
-        if len(gate.control_qubits) > 1:
-            raise_error(
-                NotImplementedError,
-                "Cannot calculate controlled "
-                "unitary for more than two "
-                "control qubits.",
-            )
-        matrix = gate.matrix(self)
-        shape = matrix.shape
-        if shape != (2, 2):
-            raise_error(
-                ValueError,
-                "Cannot use ``control_unitary`` method on "
-                + f"gate matrix of shape {shape}.",
-            )
-        zeros = self.np.zeros((2, 2), dtype=self.dtype)
-        zeros = self.cast(zeros, dtype=zeros.dtype)
-        identity = self.np.eye(2, dtype=self.dtype)
-        identity = self.cast(identity, dtype=identity.dtype)
-        part1 = self.np.concatenate([identity, zeros], axis=0)
-        part2 = self.np.concatenate([zeros, matrix], axis=0)
-        return self.np.concatenate([part1, part2], axis=1)
-
     def apply_gate(self, gate, state, nqubits):
         from qibo.gates.abstract import ParametrizedGate
+        from qibo.gates.gates import Unitary
 
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
-        # matrix = gate.matrix(self)
-        matrix = (
-            gate.__class__(*gate.qubits, *gate.parameters)
-            if isinstance(gate, ParametrizedGate)
-            else gate.__class__(*gate.qubits)
-        )
+        if isinstance(gate, Unitary):
+            matrix = gate.__class__(gate.init_args[0], *(gate.init_args[1:]))
+        else:
+            matrix = (
+                gate.__class__(*gate.init_args, *gate.parameters)
+                if isinstance(gate, ParametrizedGate)
+                else gate.__class__(*gate.init_args)
+            )
         matrix = matrix.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
@@ -222,15 +201,18 @@ def apply_gate(self, gate, state, nqubits):
 
     def apply_gate_density_matrix(self, gate, state, nqubits):
         from qibo.gates.abstract import ParametrizedGate
+        from qibo.gates.gates import Unitary
 
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
-        # matrix = gate.matrix(self)
-        matrix = (
-            gate.__class__(*gate.qubits, *gate.parameters)
-            if isinstance(gate, ParametrizedGate)
-            else gate.__class__(*gate.qubits)
-        )
+        if isinstance(gate, Unitary):
+            matrix = gate.__class__(gate.parameters[0], *gate.target_qubits)
+        else:
+            matrix = (
+                gate.__class__(*gate.target_qubits, *gate.parameters)
+                if isinstance(gate, ParametrizedGate)
+                else gate.__class__(*gate.target_qubits)
+            )
         matrix = matrix.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
diff --git a/tests/test_backends.py b/tests/test_backends.py
index 876e931b2f..c59654f9dc 100644
--- a/tests/test_backends.py
+++ b/tests/test_backends.py
@@ -89,38 +89,6 @@ def test_matrix_rotations(backend, gate, target_matrix):
     backend.assert_allclose(gate.matrix(backend), target_matrix(theta))
 
 
-def test_control_matrix(backend):
-    theta = 0.1234
-    rotation = np.array(
-        [
-            [np.cos(theta / 2.0), -np.sin(theta / 2.0)],
-            [np.sin(theta / 2.0), np.cos(theta / 2.0)],
-        ]
-    )
-    target_matrix = np.eye(4, dtype=rotation.dtype)
-    target_matrix[2:, 2:] = rotation
-    gate = gates.RY(0, theta).controlled_by(1)
-    backend.assert_allclose(gate.matrix(backend), target_matrix)
-
-    gate = gates.RY(0, theta).controlled_by(1, 2)
-    with pytest.raises(NotImplementedError):
-        matrix = backend.control_matrix(gate)
-
-
-def test_control_matrix_unitary(backend):
-    u = np.random.random((2, 2))
-    gate = gates.Unitary(u, 0).controlled_by(1)
-    matrix = backend.control_matrix(gate)
-    target_matrix = np.eye(4, dtype=np.complex128)
-    target_matrix[2:, 2:] = u
-    backend.assert_allclose(matrix, target_matrix)
-
-    u = np.random.random((16, 16))
-    gate = gates.Unitary(u, 0, 1, 2, 3).controlled_by(4)
-    with pytest.raises(ValueError):
-        matrix = backend.control_matrix(gate)
-
-
 def test_plus_density_matrix(backend):
     matrix = backend.plus_density_matrix(4)
     target_matrix = np.ones((16, 16)) / 16
diff --git a/tests/test_gates_gates.py b/tests/test_gates_gates.py
index 372bf0b202..a966ddad7a 100644
--- a/tests/test_gates_gates.py
+++ b/tests/test_gates_gates.py
@@ -1574,23 +1574,6 @@ def test_controlled_unitary(backend):
     backend.assert_allclose(final_state, target_state)
 
 
-def test_controlled_unitary_matrix(backend):
-    nqubits = 2
-    initial_state = random_statevector(2**nqubits, backend=backend)
-
-    matrix = np.random.random((2, 2))
-    gate = gates.Unitary(matrix, 1).controlled_by(0)
-
-    target_state = apply_gates(backend, [gate], nqubits, initial_state)
-
-    u = backend.control_matrix(gate)
-    u = backend.cast(u, dtype=u.dtype)
-
-    final_state = np.dot(u, initial_state)
-
-    backend.assert_allclose(final_state, target_state)
-
-
 ###############################################################################
 
 ################################# Test dagger #################################

From e1b2e80ffe4715d457efce2a4cbb126f20da9a40 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Fri, 21 Jun 2024 11:19:51 +0400
Subject: [PATCH 05/16] fix `apply_gate_density_matrix`

---
 src/qibo/backends/numpy.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 67a647b0b9..1c99aefa58 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -206,12 +206,12 @@ def apply_gate_density_matrix(self, gate, state, nqubits):
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
         if isinstance(gate, Unitary):
-            matrix = gate.__class__(gate.parameters[0], *gate.target_qubits)
+            matrix = gate.__class__(gate.init_args[0], *(gate.init_args[1:]))
         else:
             matrix = (
-                gate.__class__(*gate.target_qubits, *gate.parameters)
+                gate.__class__(*gate.init_args, *gate.parameters)
                 if isinstance(gate, ParametrizedGate)
-                else gate.__class__(*gate.target_qubits)
+                else gate.__class__(*gate.init_args)
             )
         matrix = matrix.matrix(self)
         if gate.is_controlled_by:

From 4e50a32d196f4661b08dd16674b2a5b29ec7e5f5 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Fri, 21 Jun 2024 13:21:05 +0400
Subject: [PATCH 06/16] test solution

---
 src/qibo/backends/numpy.py | 33 +++++++++++++++++++++++----------
 1 file changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 1c99aefa58..80b3b839d9 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -168,7 +168,13 @@ def apply_gate(self, gate, state, nqubits):
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
         if isinstance(gate, Unitary):
-            matrix = gate.__class__(gate.init_args[0], *(gate.init_args[1:]))
+            matrix = gate.__class__(
+                gate.init_args[0],
+                *(gate.init_args[1:]),
+                trainable=gate.init_kwargs["trainable"],
+                name=gate.init_kwargs["name"],
+                check_unitary=gate.init_kwargs["check_unitary"],
+            )
         else:
             matrix = (
                 gate.__class__(*gate.init_args, *gate.parameters)
@@ -205,15 +211,22 @@ def apply_gate_density_matrix(self, gate, state, nqubits):
 
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
-        if isinstance(gate, Unitary):
-            matrix = gate.__class__(gate.init_args[0], *(gate.init_args[1:]))
-        else:
-            matrix = (
-                gate.__class__(*gate.init_args, *gate.parameters)
-                if isinstance(gate, ParametrizedGate)
-                else gate.__class__(*gate.init_args)
-            )
-        matrix = matrix.matrix(self)
+        # if isinstance(gate, Unitary):
+        #     matrix = gate.__class__(
+        #         gate.init_args[0],
+        #         *(gate.init_args[1:]),
+        #         trainable=gate.init_kwargs["trainable"],
+        #         name=gate.init_kwargs["name"],
+        #         check_unitary=gate.init_kwargs["check_unitary"],
+        #     )
+        # else:
+        #     matrix = (
+        #         gate.__class__(*gate.init_args, *gate.parameters)
+        #         if isinstance(gate, ParametrizedGate)
+        #         else gate.__class__(*gate.init_args)
+        #     )
+        # matrix = matrix.matrix(self)
+        matrix = gate.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
             matrixc = self.np.conj(matrix)

From 7acfbe4e42ddd347cf78b494e736d79fe4659b5d Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Fri, 21 Jun 2024 13:21:19 +0400
Subject: [PATCH 07/16] remove comments

---
 src/qibo/backends/numpy.py | 15 ---------------
 1 file changed, 15 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 80b3b839d9..e16b734641 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -211,21 +211,6 @@ def apply_gate_density_matrix(self, gate, state, nqubits):
 
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
-        # if isinstance(gate, Unitary):
-        #     matrix = gate.__class__(
-        #         gate.init_args[0],
-        #         *(gate.init_args[1:]),
-        #         trainable=gate.init_kwargs["trainable"],
-        #         name=gate.init_kwargs["name"],
-        #         check_unitary=gate.init_kwargs["check_unitary"],
-        #     )
-        # else:
-        #     matrix = (
-        #         gate.__class__(*gate.init_args, *gate.parameters)
-        #         if isinstance(gate, ParametrizedGate)
-        #         else gate.__class__(*gate.init_args)
-        #     )
-        # matrix = matrix.matrix(self)
         matrix = gate.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))

From 7810d38c236a48646379e7cdecae0d0d5117863d Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Fri, 21 Jun 2024 13:45:28 +0400
Subject: [PATCH 08/16] remove extra control definition

---
 src/qibo/backends/numpy.py |  7 ++----
 src/qibo/gates/gates.py    | 51 --------------------------------------
 tests/test_gates_gates.py  | 36 +++++++++------------------
 3 files changed, 14 insertions(+), 80 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index e16b734641..64a456c117 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -162,8 +162,8 @@ def matrix_fused(self, fgate):
         return self.cast(matrix.toarray())
 
     def apply_gate(self, gate, state, nqubits):
-        from qibo.gates.abstract import ParametrizedGate
-        from qibo.gates.gates import Unitary
+        from qibo.gates.abstract import ParametrizedGate  # pylint: disable=C0415
+        from qibo.gates.gates import Unitary  # pylint: disable=C0415
 
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
@@ -206,9 +206,6 @@ def apply_gate(self, gate, state, nqubits):
         return self.np.reshape(state, (2**nqubits,))
 
     def apply_gate_density_matrix(self, gate, state, nqubits):
-        from qibo.gates.abstract import ParametrizedGate
-        from qibo.gates.gates import Unitary
-
         state = self.cast(state)
         state = self.np.reshape(state, 2 * nqubits * (2,))
         matrix = gate.matrix(self)
diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py
index 5b995b9e69..324989bfcb 100644
--- a/src/qibo/gates/gates.py
+++ b/src/qibo/gates/gates.py
@@ -71,17 +71,6 @@ def clifford(self):
     def qasm_label(self):
         return "x"
 
-    @Gate.check_controls
-    def controlled_by(self, *q):
-        """Fall back to CNOT and Toffoli if there is one or two controls."""
-        if len(q) == 1:
-            gate = CNOT(q[0], self.target_qubits[0])
-        elif len(q) == 2:
-            gate = TOFFOLI(q[0], q[1], self.target_qubits[0])
-        else:
-            gate = super().controlled_by(*q)
-        return gate
-
     def decompose(self, *free, use_toffolis=True):
         """Decomposes multi-control ``X`` gate to one-qubit, ``CNOT`` and ``TOFFOLI`` gates.
 
@@ -187,15 +176,6 @@ def clifford(self):
     def qasm_label(self):
         return "y"
 
-    @Gate.check_controls
-    def controlled_by(self, *q):
-        """Fall back to CY if there is only one control."""
-        if len(q) == 1:
-            gate = CY(q[0], self.target_qubits[0])
-        else:
-            gate = super().controlled_by(*q)
-        return gate
-
     def basis_rotation(self):
         from qibo import matrices  # pylint: disable=C0415
 
@@ -234,15 +214,6 @@ def clifford(self):
     def qasm_label(self):
         return "z"
 
-    @Gate.check_controls
-    def controlled_by(self, *q):
-        """Fall back to CZ if there is only one control."""
-        if len(q) == 1:
-            gate = CZ(q[0], self.target_qubits[0])
-        else:
-            gate = super().controlled_by(*q)
-        return gate
-
     def basis_rotation(self):
         return None
 
@@ -571,17 +542,6 @@ def _dagger(self) -> "Gate":
             self.target_qubits[0], -self.parameters[0]
         )  # pylint: disable=E1130
 
-    @Gate.check_controls
-    def controlled_by(self, *q):
-        """Fall back to CRn if there is only one control."""
-        if len(q) == 1:
-            gate = self._controlled_gate(  # pylint: disable=E1102
-                q[0], self.target_qubits[0], **self.init_kwargs
-            )
-        else:
-            gate = super().controlled_by(*q)
-        return gate
-
 
 class RX(_Rn_):
     """Rotation around the X-axis of the Bloch sphere.
@@ -847,17 +807,6 @@ def __init__(self, q, trainable=True):
 
         self.init_kwargs = {"trainable": trainable}
 
-    @Gate.check_controls
-    def controlled_by(self, *q):
-        """Fall back to CUn if there is only one control."""
-        if len(q) == 1:
-            gate = self._controlled_gate(  # pylint: disable=E1102
-                q[0], self.target_qubits[0], **self.init_kwargs
-            )
-        else:
-            gate = super().controlled_by(*q)
-        return gate
-
 
 class U1(_Un_):
     """First general unitary gate.
diff --git a/tests/test_gates_gates.py b/tests/test_gates_gates.py
index a966ddad7a..34b215129b 100644
--- a/tests/test_gates_gates.py
+++ b/tests/test_gates_gates.py
@@ -521,8 +521,7 @@ def test_cnot(backend, applyx):
 
 @pytest.mark.parametrize("seed_observable", list(range(1, 10 + 1)))
 @pytest.mark.parametrize("seed_state", list(range(1, 10 + 1)))
-@pytest.mark.parametrize("controlled_by", [False, True])
-def test_cy(backend, controlled_by, seed_state, seed_observable):
+def test_cy(backend, seed_state, seed_observable):
     nqubits = 2
     initial_state = random_statevector(2**nqubits, seed=seed_state, backend=backend)
     matrix = np.array(
@@ -544,15 +543,10 @@ def test_cy(backend, controlled_by, seed_state, seed_observable):
         initial_state=initial_state,
     )
 
-    if controlled_by:
-        gate = gates.Y(1).controlled_by(0)
-    else:
-        gate = gates.CY(0, 1)
+    gate = gates.CY(0, 1)
 
     final_state = apply_gates(backend, [gate], initial_state=initial_state)
 
-    assert gate.name == "cy"
-
     backend.assert_allclose(final_state, target_state)
 
     # testing random expectation value due to global phase difference
@@ -566,15 +560,15 @@ def test_cy(backend, controlled_by, seed_state, seed_observable):
         @ backend.cast(target_state),
     )
 
-    assert gates.CY(0, 1).qasm_label == "cy"
-    assert gates.CY(0, 1).clifford
-    assert gates.CY(0, 1).unitary
+    assert gate.name == "cy"
+    assert gate.qasm_label == "cy"
+    assert gate.clifford
+    assert gate.unitary
 
 
 @pytest.mark.parametrize("seed_observable", list(range(1, 10 + 1)))
 @pytest.mark.parametrize("seed_state", list(range(1, 10 + 1)))
-@pytest.mark.parametrize("controlled_by", [False, True])
-def test_cz(backend, controlled_by, seed_state, seed_observable):
+def test_cz(backend, seed_state, seed_observable):
     nqubits = 2
     initial_state = random_statevector(2**nqubits, seed=seed_state, backend=backend)
     matrix = np.eye(4)
@@ -590,15 +584,10 @@ def test_cz(backend, controlled_by, seed_state, seed_observable):
         initial_state=initial_state,
     )
 
-    if controlled_by:
-        gate = gates.Z(1).controlled_by(0)
-    else:
-        gate = gates.CZ(0, 1)
+    gate = gates.CZ(0, 1)
 
     final_state = apply_gates(backend, [gate], initial_state=initial_state)
 
-    assert gate.name == "cz"
-
     backend.assert_allclose(final_state, target_state)
 
     # testing random expectation value due to global phase difference
@@ -612,9 +601,10 @@ def test_cz(backend, controlled_by, seed_state, seed_observable):
         @ backend.cast(target_state),
     )
 
-    assert gates.CZ(0, 1).qasm_label == "cz"
-    assert gates.CZ(0, 1).clifford
-    assert gates.CZ(0, 1).unitary
+    assert gate.name == "cz"
+    assert gate.qasm_label == "cz"
+    assert gate.clifford
+    assert gate.unitary
 
 
 def test_csx(backend):
@@ -1457,8 +1447,6 @@ def test_controlled_u1(backend):
     target_state = np.zeros_like(final_state)
     target_state[1] = np.exp(1j * theta)
     backend.assert_allclose(final_state, target_state)
-    gate = gates.U1(0, theta).controlled_by(1)
-    assert gate.__class__.__name__ == "CU1"
 
 
 def test_controlled_u2(backend):

From ba1ed764b191317547033678464a471bc0ba0ce4 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 09:08:10 +0400
Subject: [PATCH 09/16] new solution

---
 src/qibo/backends/numpy.py | 37 +++++++------------------------------
 1 file changed, 7 insertions(+), 30 deletions(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 64a456c117..1fa121b944 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -110,25 +110,12 @@ def matrix(self, gate):
         _matrix = getattr(self.matrices, name)
         if callable(_matrix):
             _matrix = _matrix(2 ** len(gate.target_qubits))
-
-        num_controls = len(gate.control_qubits)
-        if num_controls > 0:
-            _matrix = block_diag(
-                np.eye(2 ** len(gate.qubits) - len(_matrix)), self.to_numpy(_matrix)
-            )
-
         return self.cast(_matrix, dtype=_matrix.dtype)
 
     def matrix_parametrized(self, gate):
         """Convert a parametrized gate to its matrix representation in the computational basis."""
         name = gate.__class__.__name__
         _matrix = getattr(self.matrices, name)(*gate.parameters)
-
-        num_controls = len(gate.control_qubits)
-        if num_controls > 0:
-            _matrix = block_diag(
-                np.eye(2 ** len(gate.qubits) - len(_matrix)), self.to_numpy(_matrix)
-            )
         return self.cast(_matrix, dtype=_matrix.dtype)
 
     def matrix_fused(self, fgate):
@@ -140,6 +127,13 @@ def matrix_fused(self, fgate):
             # small tensor calculations
             # explicit to_numpy see https://github.com/qiboteam/qibo/issues/928
             gmatrix = self.to_numpy(gate.matrix(self))
+            # add controls if controls were instantiated using
+            # the ``Gate.controlled_by`` method
+            num_controls = len(gate.control_qubits)
+            if num_controls > 0:
+                gmatrix = block_diag(
+                    np.eye(2 ** len(gate.qubits) - len(gmatrix)), gmatrix
+                )
             # Kronecker product with identity is needed to make the
             # original matrix have shape (2**rank x 2**rank)
             eye = np.eye(2 ** (rank - len(gate.qubits)))
@@ -162,25 +156,8 @@ def matrix_fused(self, fgate):
         return self.cast(matrix.toarray())
 
     def apply_gate(self, gate, state, nqubits):
-        from qibo.gates.abstract import ParametrizedGate  # pylint: disable=C0415
-        from qibo.gates.gates import Unitary  # pylint: disable=C0415
-
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
-        if isinstance(gate, Unitary):
-            matrix = gate.__class__(
-                gate.init_args[0],
-                *(gate.init_args[1:]),
-                trainable=gate.init_kwargs["trainable"],
-                name=gate.init_kwargs["name"],
-                check_unitary=gate.init_kwargs["check_unitary"],
-            )
-        else:
-            matrix = (
-                gate.__class__(*gate.init_args, *gate.parameters)
-                if isinstance(gate, ParametrizedGate)
-                else gate.__class__(*gate.init_args)
-            )
         matrix = matrix.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))

From 02d438afeac631e6bd7a49d251e3352cc1c4ba7e Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 09:14:20 +0400
Subject: [PATCH 10/16] fix issue

---
 src/qibo/backends/numpy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 1fa121b944..38ecf991b3 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -158,7 +158,7 @@ def matrix_fused(self, fgate):
     def apply_gate(self, gate, state, nqubits):
         state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
-        matrix = matrix.matrix(self)
+        matrix = gate.matrix(self)
         if gate.is_controlled_by:
             matrix = self.np.reshape(matrix, 2 * len(gate.target_qubits) * (2,))
             ncontrol = len(gate.control_qubits)

From 448952ec93459885a9b0e7bca873631dd6155687 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 09:27:36 +0400
Subject: [PATCH 11/16] revert one change

---
 src/qibo/gates/gates.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py
index 324989bfcb..e6d2299038 100644
--- a/src/qibo/gates/gates.py
+++ b/src/qibo/gates/gates.py
@@ -71,6 +71,17 @@ def clifford(self):
     def qasm_label(self):
         return "x"
 
+    @Gate.check_controls
+    def controlled_by(self, *q):
+        """Fall back to CNOT and Toffoli if there is one or two controls."""
+        if len(q) == 1:
+            gate = CNOT(q[0], self.target_qubits[0])
+        elif len(q) == 2:
+            gate = TOFFOLI(q[0], q[1], self.target_qubits[0])
+        else:
+            gate = super().controlled_by(*q)
+        return gate
+
     def decompose(self, *free, use_toffolis=True):
         """Decomposes multi-control ``X`` gate to one-qubit, ``CNOT`` and ``TOFFOLI`` gates.
 

From 0f98825f1e4e3677d1adefef5f37c8b95e47f8fa Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 09:45:33 +0400
Subject: [PATCH 12/16] revert changes

---
 src/qibo/gates/gates.py | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py
index e6d2299038..5b995b9e69 100644
--- a/src/qibo/gates/gates.py
+++ b/src/qibo/gates/gates.py
@@ -187,6 +187,15 @@ def clifford(self):
     def qasm_label(self):
         return "y"
 
+    @Gate.check_controls
+    def controlled_by(self, *q):
+        """Fall back to CY if there is only one control."""
+        if len(q) == 1:
+            gate = CY(q[0], self.target_qubits[0])
+        else:
+            gate = super().controlled_by(*q)
+        return gate
+
     def basis_rotation(self):
         from qibo import matrices  # pylint: disable=C0415
 
@@ -225,6 +234,15 @@ def clifford(self):
     def qasm_label(self):
         return "z"
 
+    @Gate.check_controls
+    def controlled_by(self, *q):
+        """Fall back to CZ if there is only one control."""
+        if len(q) == 1:
+            gate = CZ(q[0], self.target_qubits[0])
+        else:
+            gate = super().controlled_by(*q)
+        return gate
+
     def basis_rotation(self):
         return None
 
@@ -553,6 +571,17 @@ def _dagger(self) -> "Gate":
             self.target_qubits[0], -self.parameters[0]
         )  # pylint: disable=E1130
 
+    @Gate.check_controls
+    def controlled_by(self, *q):
+        """Fall back to CRn if there is only one control."""
+        if len(q) == 1:
+            gate = self._controlled_gate(  # pylint: disable=E1102
+                q[0], self.target_qubits[0], **self.init_kwargs
+            )
+        else:
+            gate = super().controlled_by(*q)
+        return gate
+
 
 class RX(_Rn_):
     """Rotation around the X-axis of the Bloch sphere.
@@ -818,6 +847,17 @@ def __init__(self, q, trainable=True):
 
         self.init_kwargs = {"trainable": trainable}
 
+    @Gate.check_controls
+    def controlled_by(self, *q):
+        """Fall back to CUn if there is only one control."""
+        if len(q) == 1:
+            gate = self._controlled_gate(  # pylint: disable=E1102
+                q[0], self.target_qubits[0], **self.init_kwargs
+            )
+        else:
+            gate = super().controlled_by(*q)
+        return gate
+
 
 class U1(_Un_):
     """First general unitary gate.

From 533ba0e19e02d3748c4f2a61dcfcc877b4438904 Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 10:19:28 +0400
Subject: [PATCH 13/16] update documentation

---
 src/qibo/gates/abstract.py | 59 +++++++++++++++++++++++++++++---------
 1 file changed, 45 insertions(+), 14 deletions(-)

diff --git a/src/qibo/gates/abstract.py b/src/qibo/gates/abstract.py
index 84f66c963f..a4eac1571c 100644
--- a/src/qibo/gates/abstract.py
+++ b/src/qibo/gates/abstract.py
@@ -138,8 +138,7 @@ def control_qubits(self) -> Tuple[int]:
 
     @property
     def qubits(self) -> Tuple[int]:
-        """Tuple with ids of all qubits (control and target) that the gate
-        acts."""
+        """Tuple with ids of all qubits (control and target) that the gate acts."""
         return self.control_qubits + self.target_qubits
 
     @property
@@ -197,8 +196,7 @@ def _set_targets_and_controls(
 
     @staticmethod
     def _find_repeated(qubits: Sequence[int]) -> int:
-        """Finds the first qubit id that is repeated in a sequence of qubit
-        ids."""
+        """Finds the first qubit id that is repeated in a sequence of qubit ids."""
         temp_set = set()
         for qubit in qubits:
             if qubit in temp_set:
@@ -213,14 +211,14 @@ def _check_control_target_overlap(self):
         if common:
             raise_error(
                 ValueError,
-                f"{set(self._target_qubits) & set(self._control_qubits)} qubits are both targets and controls "
+                f"{set(self._target_qubits) & set(self._control_qubits)}"
+                + "qubits are both targets and controls "
                 + f"for gate {self.__class__.__name__}.",
             )
 
     @property
     def parameters(self):
-        """Returns a tuple containing the current value of gate's
-        parameters."""
+        """Returns a tuple containing the current value of gate's parameters."""
         return self._parameters
 
     def commutes(self, gate: "Gate") -> bool:
@@ -230,7 +228,7 @@ def commutes(self, gate: "Gate") -> bool:
             gate: Gate to check if it commutes with the current gate.
 
         Returns:
-            ``True`` if the gates commute, otherwise ``False``.
+            bool: ``True`` if the gates commute, ``False`` otherwise.
         """
         if isinstance(gate, SpecialGate):  # pragma: no cover
             return False
@@ -297,8 +295,7 @@ def dagger(self) -> "Gate":
         action of dagger will be lost.
 
         Returns:
-            A :class:`qibo.gates.Gate` object representing the dagger of
-            the original gate.
+            :class:`qibo.gates.Gate`: object representing the dagger of the original gate.
         """
         new_gate = self._dagger()
         new_gate.is_controlled_by = self.is_controlled_by
@@ -322,12 +319,24 @@ def wrapper(self, *args):
     def controlled_by(self, *qubits: int) -> "Gate":
         """Controls the gate on (arbitrarily many) qubits.
 
+        To see how this method affects the underlying matrix representation of a gate,
+        please see the documentation of :meth:`qibo.gates.Gate.matrix`.
+
+        .. note::
+            Some gate classes default to another gate class depending on the number of controls
+            present. For instance, an :math:`1`-controlled :class:`qibo.gates.X` gate
+            will default to a :class:`qibo.gates.CNOT` gate, while a :math:`1`-controlled
+            :class:`qibo.gates.X` gate defaults to a :class:`qibo.gates.TOFFOLI` gate.
+            Other gates affected by this method are: :class:`qibo.gates.Y`, :class:`qibo.gates.Z`,
+            :class:`qibo.gates.RX`, :class:`qibo.gates.RY`, :class:`qibo.gates.RZ`,
+            :class:`qibo.gates.U1`, :class:`qibo.gates.U2`, and :class:`qibo.gates.U3`.
+
         Args:
             *qubits (int): Ids of the qubits that the gate will be controlled on.
 
         Returns:
-            A :class:`qibo.gates.Gate` object in with the corresponding
-            gate being controlled in the given qubits.
+            :class:`qibo.gates.Gate`: object in with the corresponding
+                gate being controlled in the given qubits.
         """
         if qubits:
             self.is_controlled_by = True
@@ -343,7 +352,7 @@ def decompose(self, *free) -> List["Gate"]:
             free: Ids of free qubits to use for the gate decomposition.
 
         Returns:
-            List with gates that have the same effect as applying the original gate.
+            list: gates that have the same effect as applying the original gate.
         """
         # TODO: Implement this method for all gates not supported by OpenQASM.
         # Currently this is implemented only for multi-controlled X gates.
@@ -354,6 +363,28 @@ def decompose(self, *free) -> List["Gate"]:
     def matrix(self, backend=None):
         """Returns the matrix representation of the gate.
 
+        If gate has controlled qubits inserted by :meth:`qibo.gates.Gate.controlled_by`,
+        then :meth:`qibo.gates.Gate.matrix` returns the matrix of the original gate.
+
+        .. code-block:: python
+
+            from qibo import gates
+
+            gate = gates.SWAP(3, 4).controlled_by(0, 1, 2)
+            print(gate.matrix())
+
+        To return the full matrix that takes the control qubits into account,
+        one should use :meth:`qibo.models.Circuit.unitary`, e.g.
+
+        .. code-block:: python
+
+            from qibo import Circuit, gates
+
+            nqubits = 5
+            circuit = Circuit(nqubits)
+            circuit.add(gates.SWAP(3, 4).controlled_by(0, 1, 2))
+            print(circuit.unitary())
+
         Args:
             backend (:class:`qibo.backends.abstract.Backend`, optional): backend
                 to be used in the execution. If ``None``, it uses
@@ -374,7 +405,7 @@ def generator_eigenvalue(self):
         """This function returns the eigenvalues of the gate's generator.
 
         Returns:
-            (float) eigenvalue of the generator.
+            float: eigenvalue of the generator.
         """
 
         raise_error(

From 5717858a52f36fd43d9df26cbeb653f4af80b2ef Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 10:41:39 +0400
Subject: [PATCH 14/16] update documentation

---
 src/qibo/models/circuit.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/qibo/models/circuit.py b/src/qibo/models/circuit.py
index 67ec3a670b..930626cdf0 100644
--- a/src/qibo/models/circuit.py
+++ b/src/qibo/models/circuit.py
@@ -737,13 +737,17 @@ def gate_names(self) -> collections.Counter:
     def gates_of_type(self, gate: Union[str, type]) -> List[Tuple[int, gates.Gate]]:
         """Finds all gate objects of specific type or name.
 
+        This method can be affected by how :meth:`qibo.gates.Gate.controlled_by`
+        behaves with certain gates. To see how :meth:`qibo.gates.Gate.controlled_by`
+        affects gates, we refer to the documentation of :meth:`qibo.gates.Gate.controlled_by`.
+
         Args:
-            gate (str, type): The QASM name of a gate or the corresponding gate class.
+            gate (str or type): The name of a gate or the corresponding gate class.
 
         Returns:
-            List with all gates that are in the circuit and have the same type
-            with the given ``gate``. The list contains tuples ``(i, g)`` where
-            ``i`` is the index of the gate ``g`` in the circuit's gate queue.
+            list: gates that are in the circuit and have the same type as ``gate``.
+                The list contains tuples ``(k, g)`` where ``k`` is the index of the gate
+                ``g`` in the circuit's gate queue.
         """
         if isinstance(gate, str):
             return [(i, g) for i, g in enumerate(self.queue) if g.name == gate]

From b7a821249902dcc6ebce1676971d3fd520c7a4cb Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Mon, 24 Jun 2024 15:23:12 +0000
Subject: [PATCH 15/16] Update src/qibo/gates/abstract.py

Co-authored-by: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
---
 src/qibo/gates/abstract.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibo/gates/abstract.py b/src/qibo/gates/abstract.py
index a4eac1571c..f222c027ed 100644
--- a/src/qibo/gates/abstract.py
+++ b/src/qibo/gates/abstract.py
@@ -325,7 +325,7 @@ def controlled_by(self, *qubits: int) -> "Gate":
         .. note::
             Some gate classes default to another gate class depending on the number of controls
             present. For instance, an :math:`1`-controlled :class:`qibo.gates.X` gate
-            will default to a :class:`qibo.gates.CNOT` gate, while a :math:`1`-controlled
+            will default to a :class:`qibo.gates.CNOT` gate, while a :math:`2`-controlled
             :class:`qibo.gates.X` gate defaults to a :class:`qibo.gates.TOFFOLI` gate.
             Other gates affected by this method are: :class:`qibo.gates.Y`, :class:`qibo.gates.Z`,
             :class:`qibo.gates.RX`, :class:`qibo.gates.RY`, :class:`qibo.gates.RZ`,

From 3637416a622505fc134266aae570cfda3bd48e0a Mon Sep 17 00:00:00 2001
From: Renato Mello <renato.msf@gmail.com>
Date: Thu, 27 Jun 2024 08:59:22 +0400
Subject: [PATCH 16/16] add test suggested by Stavros

---
 tests/test_models_circuit_features.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/tests/test_models_circuit_features.py b/tests/test_models_circuit_features.py
index f20b983c97..b8a1ade5e7 100644
--- a/tests/test_models_circuit_features.py
+++ b/tests/test_models_circuit_features.py
@@ -58,6 +58,28 @@ def test_circuit_unitary_and_inverse_with_noise_channel(backend):
         circuit.invert()
 
 
+def test_circuit_unitary_non_trivial(backend):
+    target = np.array(
+        [
+            [1, 0, 0, 0, 0, 0, 0, 0],
+            [0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 1],
+        ],
+        dtype=complex,
+    )
+    target = backend.cast(target, dtype=target.dtype)
+    nqubits = 3
+    circuit = Circuit(nqubits)
+    circuit.add(gates.SWAP(0, 2).controlled_by(1))
+    unitary = circuit.unitary(backend)
+    backend.assert_allclose(unitary, target)
+
+
 @pytest.mark.parametrize("compile", [False, True])
 def test_circuit_vs_gate_execution(backend, compile):
     """Check consistency between executing circuit and stand alone gates."""