Repository: qutip/qutip-qip (language: python)

Issue:
Align the signatures of gate functions with qutip-v5
In qutip-v5 the gate functions no longer take parameters `N` , `targets` and `control`, which permute the qubits if needed and generate gate operators in a larger Hilbert space. Although many of the gates are moved to `qutip.core.gates`, these additional parameters are removed in qutip-v5. Also, the signature of `expand_operator` is slightly changed to be more general.

To avoid confusion and code duplication, we should adapt the functions in `qutip-qip` accordingly.

Proposals:
- Import the gate functions from qutip when v5 is used. Otherwise still use the one defined in `qutip-qip`.
- Align the signatures of gate functions and `expand_operator` with qutip-v5, and issue deprecation warnings when needed.

Optional:
- In principle, we can try to keep the compatibility and define a wrapper that expands the obtained `Qobj`. This will need some complications but should be feasible.
- Do we also want to add some deprecation warning for this to the gates in qutip-v4 in the module `qutip.qip`. Some people are still using it. It would make the migration from 4 to 5 easier.

@hodgestar Do you think we should implement the optional ones?


Discussion / hints:

Fix patch:
diff --git a/doc/source/apidoc/qutip_qip.operations.rst b/doc/source/apidoc/qutip_qip.operations.rst
index 917bbb97..d239fd4c 100644
--- a/doc/source/apidoc/qutip_qip.operations.rst
+++ b/doc/source/apidoc/qutip_qip.operations.rst
@@ -64,9 +64,6 @@ qutip\_qip.operations
       cz_gate
       expand_operator
       fredkin
-      gate_expand_1toN
-      gate_expand_2toN
-      gate_expand_3toN
       gate_sequence_product
       globalphase
       hadamard_transform
diff --git a/src/qutip_qip/algorithms/qft.py b/src/qutip_qip/algorithms/qft.py
index b9c2a640..f3830a3b 100644
--- a/src/qutip_qip/algorithms/qft.py
+++ b/src/qutip_qip/algorithms/qft.py
@@ -4,7 +4,7 @@
 
 
 import numpy as np
-from ..operations import Gate, snot, cphase, swap
+from ..operations import Gate, snot, cphase, swap, expand_operator
 from ..circuit import QubitCircuit
 from qutip import Qobj
 from ..decompose import decompose_one_qubit_gate
@@ -69,13 +69,22 @@ def qft_steps(N=1, swapping=True):
         for i in range(N):
             for j in range(i):
                 U_step_list.append(
-                    cphase(np.pi / (2 ** (i - j)), N, control=i, target=j)
+                    expand_operator(
+                        cphase(np.pi / (2 ** (i - j))),
+                        dims=[2] * N,
+                        targets=[i, j],
+                    )
                 )
-            U_step_list.append(snot(N, i))
+            U_step_list.append(
+                expand_operator(snot(), dims=[2] * N, targets=i)
+            )
         if swapping:
             for i in range(N // 2):
-                U_step_list.append(swap(N, [N - i - 1, i]))
-
+                U_step_list.append(
+                    expand_operator(
+                        swap(), dims=[2] * N, targets=[N - i - 1, i]
+                    )
+                )
     return U_step_list
 
 
diff --git a/src/qutip_qip/circuit/circuit.py b/src/qutip_qip/circuit/circuit.py
index 9ef5af7e..a414103f 100644
--- a/src/qutip_qip/circuit/circuit.py
+++ b/src/qutip_qip/circuit/circuit.py
@@ -923,7 +923,7 @@ def propagators(self, expand=True, ignore_measurement=False):
                 if expand:
                     all_targets = gate.get_all_qubits()
                     qobj = expand_operator(
-                        qobj, N=self.N, targets=all_targets, dims=self.dims
+                        qobj, dims=self.dims, targets=all_targets
                     )
             else:
                 if expand:
diff --git a/src/qutip_qip/circuit/circuitsimulator.py b/src/qutip_qip/circuit/circuitsimulator.py
index d756fcd8..d39035ed 100644
--- a/src/qutip_qip/circuit/circuitsimulator.py
+++ b/src/qutip_qip/circuit/circuitsimulator.py
@@ -104,9 +104,11 @@ def _mult_sublists(tensor_list, overall_inds, U, inds):
     ind_map = {ind: pos for ind, pos in zip(revised_inds, sorted_positions)}
 
     U_sublist = expand_operator(
-        U_sublist, N, [ind_map[ind] for ind in inds_sublist]
+        U_sublist, dims=[2] * N, targets=[ind_map[ind] for ind in inds_sublist]
+    )
+    U = expand_operator(
+        U, dims=[2] * N, targets=[ind_map[ind] for ind in inds]
     )
-    U = expand_operator(U, N, [ind_map[ind] for ind in inds])
 
     U_sublist = U * U_sublist
     inds_sublist = revised_inds
@@ -125,7 +127,9 @@ def _expand_overall(tensor_list, overall_inds):
 
     U_overall = tensor(tensor_list)
     overall_inds = _flatten(overall_inds)
-    U_overall = expand_operator(U_overall, len(overall_inds), overall_inds)
+    U_overall = expand_operator(
+        U_overall, dims=[2] * len(overall_inds), targets=overall_inds
+    )
     overall_inds = sorted(overall_inds)
     return U_overall, overall_inds
 
@@ -188,7 +192,9 @@ def _gate_sequence_product(U_list, ind_list):
                 tensor_list, overall_inds
             )
             U_left, rem_inds = _gate_sequence_product(U_list[i:], ind_list[i:])
-            U_left = expand_operator(U_left, num_qubits, rem_inds)
+            U_left = expand_operator(
+                U_left, dims=[2] * num_qubits, targets=rem_inds
+            )
             return U_left * U_overall, [
                 sorted_inds[ind] for ind in overall_inds
             ]
@@ -292,6 +298,7 @@ def __init__(
         """
 
         self.qc = qc
+        self.dims = qc.dims
         self.mode = mode
         self.precompute_unitary = precompute_unitary
 
@@ -468,7 +475,7 @@ def _compute_unitary(self, U_list, inds_list):
 
         if len(overall_inds) != self.qc.N:
             U_overall = expand_operator(
-                U_overall, N=self.qc.N, targets=overall_inds
+                U_overall, dims=self.qc.dims, targets=overall_inds
             )
         return U_overall
 
@@ -573,7 +580,9 @@ def _decimal_to_binary(decimal, length):
             if apply_gate:
                 if self.precompute_unitary:
                     U = expand_operator(
-                        U, self.qc.N, operation.get_all_qubits()
+                        U,
+                        dims=self.qc.dims,
+                        targets=operation.get_all_qubits(),
                     )
                 self._evolve_state(U)
         else:
diff --git a/src/qutip_qip/device/cavityqed.py b/src/qutip_qip/device/cavityqed.py
index 4614d451..180d1a82 100644
--- a/src/qutip_qip/device/cavityqed.py
+++ b/src/qutip_qip/device/cavityqed.py
@@ -17,7 +17,6 @@
 from ..operations import Gate
 from .processor import Processor, Model
 from .modelprocessor import ModelProcessor, _to_array
-from ..operations import expand_operator
 from ..pulse import Pulse
 from ..compiler import GateCompiler, CavityQEDCompiler
 
diff --git a/src/qutip_qip/device/optpulseprocessor.py b/src/qutip_qip/device/optpulseprocessor.py
index e911aa11..b454fb9f 100644
--- a/src/qutip_qip/device/optpulseprocessor.py
+++ b/src/qutip_qip/device/optpulseprocessor.py
@@ -181,16 +181,12 @@ def load_circuit(
             for label in control_labels:
                 qobj, targets = self.model.get_control(label)
                 full_ctrls_hams.append(
-                    expand_operator(
-                        qobj, len(self.dims), targets=targets, dims=self.dims
-                    )
+                    expand_operator(qobj, dims=self.dims, targets=targets)
                 )
 
             full_drift_ham = sum(
                 [
-                    expand_operator(
-                        qobj, len(self.dims), targets=targets, dims=self.dims
-                    )
+                    expand_operator(qobj, dims=self.dims, targets=targets)
                     for (qobj, targets) in self.model.get_all_drift()
                 ],
                 Qobj(
diff --git a/src/qutip_qip/gates.py b/src/qutip_qip/gates.py
deleted file mode 100644
index 6696f3e9..00000000
--- a/src/qutip_qip/gates.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import warnings
-
-from .operations.gates import *
-
-warnings.warn(
-    "Importation from qutip_qip.gates is deprecated."
-    "Please use e.g.\n from qutip_qip.operations import cnot\n",
-    DeprecationWarning,
-    stacklevel=2,
-)
diff --git a/src/qutip_qip/operations/gateclass.py b/src/qutip_qip/operations/gateclass.py
index 0ca9136d..b5ea7d08 100644
--- a/src/qutip_qip/operations/gateclass.py
+++ b/src/qutip_qip/operations/gateclass.py
@@ -391,9 +391,8 @@ def get_qobj(self, num_qubits=None, dims=None):
             num_qubits = max(all_targets) + 1
         return expand_operator(
             self.get_compact_qobj(),
-            N=num_qubits,
-            targets=all_targets,
             dims=dims,
+            targets=all_targets,
         )
 
 
diff --git a/src/qutip_qip/operations/gates.py b/src/qutip_qip/operations/gates.py
index 902ad2e7..e26755cf 100644
--- a/src/qutip_qip/operations/gates.py
+++ b/src/qutip_qip/operations/gates.py
@@ -48,9 +48,6 @@
     "controlled_gate",
     "globalphase",
     "hadamard_transform",
-    "gate_expand_1toN",
-    "gate_expand_2toN",
-    "gate_expand_3toN",
     "qubit_clifford_group",
     "expand_operator",
     "gate_sequence_product",
@@ -60,6 +57,15 @@
 #
 # Single Qubit Gates
 #
+def _deprecation_warnings_gate_expansion():
+    warnings.warn(
+        "The expansion of output gate matrix is no longer included "
+        "in the gate functions. "
+        "To expand the output `Qobj` or permute the qubits, "
+        "please use expand_operator.",
+        DeprecationWarning,
+        stacklevel=2,
+    )
 
 
 def x_gate(N=None, target=0):
@@ -73,7 +79,8 @@ def x_gate(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(x_gate(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(x_gate(), dims=[2] * N, targets=target)
     return sigmax()
 
 
@@ -88,7 +95,8 @@ def y_gate(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(y_gate(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(y_gate(), dims=[2] * N, targets=target)
     return sigmay()
 
 
@@ -105,7 +113,10 @@ def cy_gate(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(cy_gate(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            cy_gate(), dims=[2] * N, targets=(control, target)
+        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1j], [0, 0, 1j, 0]],
         dims=[[2, 2], [2, 2]],
@@ -123,7 +134,8 @@ def z_gate(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(z_gate(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(z_gate(), dims=[2] * N, targets=target)
     return sigmaz()
 
 
@@ -140,7 +152,10 @@ def cz_gate(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(cz_gate(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            cz_gate(), dims=[2] * N, targets=(control, target)
+        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
         dims=[[2, 2], [2, 2]],
@@ -158,7 +173,8 @@ def s_gate(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(s_gate(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(s_gate(), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, 1j]])
 
 
@@ -175,7 +191,10 @@ def cs_gate(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(cs_gate(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            cs_gate(), dims=[2] * N, targets=(control, target)
+        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1j]],
         dims=[[2, 2], [2, 2]],
@@ -192,7 +211,8 @@ def t_gate(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(t_gate(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(t_gate(), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, np.exp(1j * np.pi / 4)]])
 
 
@@ -209,7 +229,10 @@ def ct_gate(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(ct_gate(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            ct_gate(), dims=[2] * N, targets=(control, target)
+        )
     return Qobj(
         [
             [1, 0, 0, 0],
@@ -231,7 +254,8 @@ def rx(phi, N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(rx(phi), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(rx(phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [np.cos(phi / 2), -1j * np.sin(phi / 2)],
@@ -250,7 +274,8 @@ def ry(phi, N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(ry(phi), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(ry(phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [np.cos(phi / 2), -np.sin(phi / 2)],
@@ -269,7 +294,8 @@ def rz(phi, N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(rz(phi), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(rz(phi), dims=[2] * N, targets=target)
     return Qobj([[np.exp(-1j * phi / 2), 0], [0, np.exp(1j * phi / 2)]])
 
 
@@ -283,7 +309,8 @@ def sqrtnot(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(sqrtnot(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(sqrtnot(), dims=[2] * N, targets=target)
     return Qobj([[0.5 + 0.5j, 0.5 - 0.5j], [0.5 - 0.5j, 0.5 + 0.5j]])
 
 
@@ -306,7 +333,8 @@ def snot(N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(snot(), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(snot(), dims=[2] * N, targets=target)
     return 1 / np.sqrt(2.0) * Qobj([[1, 1], [1, -1]])
 
 
@@ -335,7 +363,8 @@ def phasegate(theta, N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(phasegate(theta), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(phasegate(theta), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, np.exp(1.0j * theta)]], dims=[[2], [2]])
 
 
@@ -361,7 +390,8 @@ def qrot(theta, phi, N=None, target=0):
         a rabi pulse.
     """
     if N is not None:
-        return expand_operator(qrot(theta, phi), N=N, targets=target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(qrot(theta, phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [
@@ -403,8 +433,9 @@ def qasmu_gate(args, N=None, target=0):
 
     theta, phi, gamma = args
     if N is not None:
+        _deprecation_warnings_gate_expansion()
         return expand_operator(
-            qasmu_gate([theta, phi, gamma]), N=N, targets=target
+            qasmu_gate([theta, phi, gamma]), dims=[2] * N, targets=target
         )
     return Qobj(rz(phi) * ry(theta) * rz(gamma))
 
@@ -437,6 +468,8 @@ def cphase(theta, N=2, control=0, target=1):
     U : qobj
         Quantum object representation of controlled phase gate.
     """
+    if N != 2 or control != 0 or target != 1:
+        _deprecation_warnings_gate_expansion()
 
     if N < 1 or target < 0 or control < 0:
         raise ValueError("Minimum value: N=1, control=0 and target=0")
@@ -481,7 +514,8 @@ def cnot(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(cnot(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(cnot(), dims=[2] * N, targets=(control, target))
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]],
         dims=[[2, 2], [2, 2]],
@@ -513,7 +547,8 @@ def csign(N=None, control=0, target=1):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(csign(), N, control, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(csign(), N, (control, target))
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
         dims=[[2, 2], [2, 2]],
@@ -545,7 +580,8 @@ def berkeley(N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(berkeley(), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(berkeley(), N, targets=targets)
     return Qobj(
         [
             [np.cos(np.pi / 8), 0, 0, 1.0j * np.sin(np.pi / 8)],
@@ -582,7 +618,8 @@ def swapalpha(alpha, N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(swapalpha(alpha), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(swapalpha(alpha), N, targets=targets)
     return Qobj(
         [
             [1, 0, 0, 0],
@@ -628,7 +665,8 @@ def swap(N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(swap(), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(swap(), dims=[2] * N, targets=targets)
     return Qobj(
         [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]],
         dims=[[2, 2], [2, 2]],
@@ -658,7 +696,8 @@ def iswap(N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(iswap(), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(iswap(), dims=[2] * N, targets=targets)
     return Qobj(
         [[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]],
         dims=[[2, 2], [2, 2]],
@@ -678,7 +717,8 @@ def sqrtswap(N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(sqrtswap(), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(sqrtswap(), dims=[2] * N, targets=targets)
     return Qobj(
         np.array(
             [
@@ -720,7 +760,8 @@ def sqrtiswap(N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return gate_expand_2toN(sqrtiswap(), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(sqrtiswap(), N, targets=targets)
     return Qobj(
         np.array(
             [
@@ -756,7 +797,10 @@ def molmer_sorensen(theta, N=None, targets=[0, 1]):
         N = 2
 
     if N is not None:
-        return expand_operator(molmer_sorensen(theta), N, targets=targets)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            molmer_sorensen(theta), dims=[2] * N, targets=targets
+        )
     return Qobj(
         [
             [np.cos(theta / 2.0), 0, 0, -1.0j * np.sin(theta / 2.0)],
@@ -801,8 +845,9 @@ def fredkin(N=None, control=0, targets=[1, 2]):
         N = 3
 
     if N is not None:
-        return gate_expand_3toN(
-            fredkin(), N, [control, targets[0]], targets[1]
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            fredkin(), dims=[2] * N, targets=(control,) + tuple(targets)
         )
     return Qobj(
         [
@@ -848,7 +893,10 @@ def toffoli(N=None, controls=[0, 1], target=2):
         N = 3
 
     if N is not None:
-        return gate_expand_3toN(toffoli(), N, controls, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(
+            toffoli(), dims=[2] * N, targets=tuple(controls) + (target,)
+        )
     return Qobj(
         [
             [1, 0, 0, 0, 0, 0, 0, 0],
@@ -879,7 +927,8 @@ def rotation(op, phi, N=None, target=0):
 
     """
     if N is not None:
-        return gate_expand_1toN(rotation(op, phi), N, target)
+        _deprecation_warnings_gate_expansion()
+        return expand_operator(rotation(op, phi), N, target)
     return (-1j * op * phi / 2).expm()
 
 
@@ -1066,7 +1115,7 @@ def qubit_clifford_group(N=None, target=0):
 
         # Finally, we optionally expand the gate.
         if N is not None:
-            yield gate_expand_1toN(op, N, target)
+            yield expand_operator(op, N, target)
         else:
             yield op
 
@@ -1076,264 +1125,7 @@ def qubit_clifford_group(N=None, target=0):
 #
 
 
-def gate_expand_1toN(U, N, target):
-    """
-    Create a Qobj representing a one-qubit gate that act on a system with N
-    qubits.
-
-    Parameters
-    ----------
-    U : Qobj
-        The one-qubit gate
-
-    N : integer
-        The number of qubits in the target space.
-
-    target : integer
-        The index of the target qubit.
-
-    Returns
-    -------
-    gate : qobj
-        Quantum object representation of N-qubit gate.
-
-    """
-
-    if N < 1:
-        raise ValueError("integer N must be larger or equal to 1")
-
-    if target >= N:
-        raise ValueError("target must be integer < integer N")
-
-    return tensor(
-        [identity(2)] * (target) + [U] + [identity(2)] * (N - target - 1)
-    )
-
-
-def gate_expand_2toN(U, N, control=None, target=None, targets=None):
-    """
-    Create a Qobj representing a two-qubit gate that act on a system with N
-    qubits.
-
-    Parameters
-    ----------
-    U : Qobj
-        The two-qubit gate
-
-    N : integer
-        The number of qubits in the target space.
-
-    control : integer
-        The index of the control qubit.
-
-    target : integer
-        The index of the target qubit.
-
-    targets : list
-        List of target qubits.
-
-    Returns
-    -------
-    gate : qobj
-        Quantum object representation of N-qubit gate.
-
-    """
-
-    if targets is not None:
-        control, target = targets
-
-    if control is None or target is None:
-        raise ValueError("Specify value of control and target")
-
-    if N < 2:
-        raise ValueError("integer N must be larger or equal to 2")
-
-    if control >= N or target >= N:
-        raise ValueError("control and not target must be integer < integer N")
-
-    if control == target:
-        raise ValueError("target and not control cannot be equal")
-
-    p = list(range(N))
-
-    if target == 0 and control == 1:
-        p[control], p[target] = p[target], p[control]
-
-    elif target == 0:
-        p[1], p[target] = p[target], p[1]
-        p[1], p[control] = p[control], p[1]
-
-    else:
-        p[1], p[target] = p[target], p[1]
-        p[0], p[control] = p[control], p[0]
-
-    return tensor([U] + [identity(2)] * (N - 2)).permute(p)
-
-
-def gate_expand_3toN(U, N, controls=[0, 1], target=2):
-    """
-    Create a Qobj representing a three-qubit gate that act on a system with N
-    qubits.
-
-    Parameters
-    ----------
-    U : Qobj
-        The three-qubit gate
-
-    N : integer
-        The number of qubits in the target space.
-
-    controls : list
-        The list of the control qubits.
-
-    target : integer
-        The index of the target qubit.
-
-    Returns
-    -------
-    gate : qobj
-        Quantum object representation of N-qubit gate.
-
-    """
-
-    if N < 3:
-        raise ValueError("integer N must be larger or equal to 3")
-
-    if controls[0] >= N or controls[1] >= N or target >= N:
-        raise ValueError(
-            "control and not target is None." " Must be integer < integer N"
-        )
-
-    if (
-        controls[0] == target
-        or controls[1] == target
-        or controls[0] == controls[1]
-    ):
-
-        raise ValueError(
-            "controls[0], controls[1], and target" " cannot be equal"
-        )
-
-    p = list(range(N))
-    p1 = list(range(N))
-    p2 = list(range(N))
-
-    if controls[0] <= 2 and controls[1] <= 2 and target <= 2:
-        p[controls[0]] = 0
-        p[controls[1]] = 1
-        p[target] = 2
-
-    #
-    # N > 3 cases
-    #
-
-    elif controls[0] == 0 and controls[1] == 1:
-        p[2], p[target] = p[target], p[2]
-
-    elif controls[0] == 0 and target == 2:
-        p[1], p[controls[1]] = p[controls[1]], p[1]
-
-    elif controls[1] == 1 and target == 2:
-        p[0], p[controls[0]] = p[controls[0]], p[0]
-
-    elif controls[0] == 1 and controls[1] == 0:
-        p[controls[1]], p[controls[0]] = p[controls[0]], p[controls[1]]
-        p2[2], p2[target] = p2[target], p2[2]
-        p = [p2[p[k]] for k in range(N)]
-
-    elif controls[0] == 2 and target == 0:
-        p[target], p[controls[0]] = p[controls[0]], p[target]
-        p1[1], p1[controls[1]] = p1[controls[1]], p1[1]
-        p = [p1[p[k]] for k in range(N)]
-
-    elif controls[1] == 2 and target == 1:
-        p[target], p[controls[1]] = p[controls[1]], p[target]
-        p1[0], p1[controls[0]] = p1[controls[0]], p1[0]
-        p = [p1[p[k]] for k in range(N)]
-
-    elif controls[0] == 1 and controls[1] == 2:
-        #  controls[0] -> controls[1] -> target -> outside
-        p[0], p[1] = p[1], p[0]
-        p[0], p[2] = p[2], p[0]
-        p[0], p[target] = p[target], p[0]
-
-    elif controls[0] == 2 and target == 1:
-        #  controls[0] -> target -> controls[1] -> outside
-        p[0], p[2] = p[2], p[0]
-        p[0], p[1] = p[1], p[0]
-        p[0], p[controls[1]] = p[controls[1]], p[0]
-
-    elif controls[1] == 0 and controls[0] == 2:
-        #  controls[1] -> controls[0] -> target -> outside
-        p[1], p[0] = p[0], p[1]
-        p[1], p[2] = p[2], p[1]
-        p[1], p[target] = p[target], p[1]
-
-    elif controls[1] == 2 and target == 0:
-        #  controls[1] -> target -> controls[0] -> outside
-        p[1], p[2] = p[2], p[1]
-        p[1], p[0] = p[0], p[1]
-        p[1], p[controls[0]] = p[controls[0]], p[1]
-
-    elif target == 1 and controls[1] == 0:
-        #  target -> controls[1] -> controls[0] -> outside
-        p[2], p[1] = p[1], p[2]
-        p[2], p[0] = p[0], p[2]
-        p[2], p[controls[0]] = p[controls[0]], p[2]
-
-    elif target == 0 and controls[0] == 1:
-        #  target -> controls[0] -> controls[1] -> outside
-        p[2], p[0] = p[0], p[2]
-        p[2], p[1] = p[1], p[2]
-        p[2], p[controls[1]] = p[controls[1]], p[2]
-
-    elif controls[0] == 0 and controls[1] == 2:
-        #  controls[0] -> self, controls[1] -> target -> outside
-        p[1], p[2] = p[2], p[1]
-        p[1], p[target] = p[target], p[1]
-
-    elif controls[1] == 1 and controls[0] == 2:
-        #  controls[1] -> self, controls[0] -> target -> outside
-        p[0], p[2] = p[2], p[0]
-        p[0], p[target] = p[target], p[0]
-
-    elif target == 2 and controls[0] == 1:
-        #  target -> self, controls[0] -> controls[1] -> outside
-        p[0], p[1] = p[1], p[0]
-        p[0], p[controls[1]] = p[controls[1]], p[0]
-
-    #
-    # N > 4 cases
-    #
-
-    elif controls[0] == 1 and controls[1] > 2 and target > 2:
-        #  controls[0] -> controls[1] -> outside, target -> outside
-        p[0], p[1] = p[1], p[0]
-        p[0], p[controls[1]] = p[controls[1]], p[0]
-        p[2], p[target] = p[target], p[2]
-
-    elif controls[0] == 2 and controls[1] > 2 and target > 2:
-        #  controls[0] -> target -> outside, controls[1] -> outside
-        p[0], p[2] = p[2], p[0]
-        p[0], p[target] = p[target], p[0]
-        p[1], p[controls[1]] = p[controls[1]], p[1]
-
-    elif controls[1] == 2 and controls[0] > 2 and target > 2:
-        #  controls[1] -> target -> outside, controls[0] -> outside
-        p[1], p[2] = p[2], p[1]
-        p[1], p[target] = p[target], p[1]
-        p[0], p[controls[0]] = p[controls[0]], p[0]
-
-    else:
-        p[0], p[controls[0]] = p[controls[0]], p[0]
-        p1[1], p1[controls[1]] = p1[controls[1]], p1[1]
-        p2[2], p2[target] = p2[target], p2[2]
-        p = [p[p1[p2[k]]] for k in range(N)]
-
-    return tensor([U] + [identity(2)] * (N - 3)).permute(p)
-
-
-def _check_qubits_oper(oper, dims=None, targets=None):
+def _check_oper_dims(oper, dims=None, targets=None):
     """
     Check if the given operator is valid.
 
@@ -1343,10 +1135,9 @@ def _check_qubits_oper(oper, dims=None, targets=None):
         The quantum object to be checked.
     dims : list, optional
         A list of integer for the dimension of each composite system.
-        e.g ``[2, 2, 2, 2, 2]`` for 5 qubits system. If None, qubits system
-        will be the default.
+        e.g ``[2, 2, 2, 2, 2]`` for 5 qubits system.
     targets : int or list of int, optional
-        The indices of qubits that are acted on.
+        The indices of subspace that are acted on.
     """
     # if operator matches N
     if not isinstance(oper, Qobj) or oper.dims[0] != oper.dims[1]:
@@ -1371,19 +1162,19 @@ def _targets_to_list(targets, oper=None, N=None):
     Parameters
     ----------
     targets : int or list of int
-        The indices of qubits that are acted on.
+        The indices of subspace that are acted on.
     oper : :class:`qutip.Qobj`, optional
-        An operator acts on qubits, the type of the :class:`qutip.Qobj`
+        An operator, the type of the :class:`qutip.Qobj`
         has to be an operator
         and the dimension matches the tensored qubit Hilbert space
         e.g. dims = ``[[2, 2, 2], [2, 2, 2]]``
     N : int, optional
-        The number of qubits in the system.
+        The number of subspace in the system.
     """
     # if targets is a list of integer
     if targets is None:
         targets = list(range(len(oper.dims[0])))
-    if not isinstance(targets, Iterable):
+    if not hasattr(targets, "__iter__"):
         targets = [targets]
     if not all([isinstance(t, numbers.Integral) for t in targets]):
         raise TypeError("targets should be an integer or a list of integer")
@@ -1403,26 +1194,28 @@ def _targets_to_list(targets, oper=None, N=None):
     return targets
 
 
-def expand_operator(oper, N, targets, dims=None, cyclic_permutation=False):
+def expand_operator(
+    oper, N=None, targets=None, dims=None, cyclic_permutation=False
+):
     """
-    Expand a qubits operator to one that acts on a N-qubit system.
+    Expand an operator to one that acts on a system with desired dimensions.
 
     Parameters
     ----------
     oper : :class:`qutip.Qobj`
-        An operator acts on qubits, the type of the :class:`qutip.Qobj`
-        has to be an operator
-        and the dimension matches the tensored qubit Hilbert space
-        e.g. dims = ``[[2, 2, 2], [2, 2, 2]]``
-    N : int
-        The number of qubits in the system.
-    targets : int or list of int
-        The indices of qubits that are acted on.
-    dims : list, optional
+        An operator that act on the subsystem, has to be an operator and the
+        dimension matches the tensored dims Hilbert space
+        e.g. oper.dims = ``[[2, 3], [2, 3]]``
+    dims : list
         A list of integer for the dimension of each composite system.
-        E.g ``[2, 2, 2, 2, 2]`` for 5 qubits system. If None, qubits system
-        will be the default option.
+        E.g ``[2, 3, 2, 3, 4]``.
+    targets : int or list of int
+        The indices of subspace that are acted on.
+        Permutation can also be realized by changing the orders of the indices.
+    N : int
+         Deprecated. Number of qubits. Please use `dims`.
     cyclic_permutation : boolean, optional
+        Deprecated.
         Expand for all cyclic permutation of the targets.
         E.g. if ``N=3`` and `oper` is a 2-qubit operator,
         the result will be a list of three operators,
@@ -1431,20 +1224,70 @@ def expand_operator(oper, N, targets, dims=None, cyclic_permutation=False):
     Returns
     -------
     expanded_oper : :class:`qutip.Qobj`
-        The expanded qubits operator acting on a system with N qubits.
+        The expanded operator acting on a system with desired dimension.
 
-    Notes
-    -----
-    This is equivalent to gate_expand_1toN, gate_expand_2toN,
-    gate_expand_3toN in ``qutip_qip.gate.py``, but works for any dimension.
+    Examples
+    --------
+    >>> from qutip_qip.operations import expand_operator, x_gate, cnot
+    >>> import qutip
+    >>> expand_operator(x_gate(), dims=[2,3], targets=[0]) # doctest: +NORMALIZE_WHITESPACE
+    Quantum object: dims = [[2, 3], [2, 3]], shape = (6, 6), type = oper, isherm = True
+    Qobj data =
+    [[0. 0. 0. 1. 0. 0.]
+     [0. 0. 0. 0. 1. 0.]
+     [0. 0. 0. 0. 0. 1.]
+     [1. 0. 0. 0. 0. 0.]
+     [0. 1. 0. 0. 0. 0.]
+     [0. 0. 1. 0. 0. 0.]]
+    >>> expand_operator(cnot(), dims=[2,2,2], targets=[1, 2]) # doctest: +NORMALIZE_WHITESPACE
+    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8,
+    8), type = oper, isherm = True
+    Qobj data =
+    [[1. 0. 0. 0. 0. 0. 0. 0.]
+     [0. 1. 0. 0. 0. 0. 0. 0.]
+     [0. 0. 0. 1. 0. 0. 0. 0.]
+     [0. 0. 1. 0. 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. 1.]
+     [0. 0. 0. 0. 0. 0. 1. 0.]]
+    >>> expand_operator(cnot(), dims=[2, 2, 2], targets=[2, 0]) # doctest: +NORMALIZE_WHITESPACE
+    Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8,
+    8), type = oper, isherm = True
+    Qobj data =
+    [[1. 0. 0. 0. 0. 0. 0. 0.]
+     [0. 0. 0. 0. 0. 1. 0. 0.]
+     [0. 0. 1. 0. 0. 0. 0. 0.]
+     [0. 0. 0. 0. 0. 0. 0. 1.]
+     [0. 0. 0. 0. 1. 0. 0. 0.]
+     [0. 1. 0. 0. 0. 0. 0. 0.]
+     [0. 0. 0. 0. 0. 0. 1. 0.]
+     [0. 0. 0. 1. 0. 0. 0. 0.]]
     """
+    if N is not None:
+        warnings.warn(
+            "The function expand_operator has been generalized to "
+            "arbitrary subsystems instead of only qubit systems."
+            "Please use the new signature e.g.\n"
+            "expand_operator(oper, dims=[2, 3, 2, 2], targets=2)",
+            DeprecationWarning,
+        )
+    if dims is not None and N is None:
+        if not isinstance(dims, Iterable):
+            f"dims needs to be an interable {not type(dims)}."
+        N = len(dims)  # backward compatibility
     if dims is None:
         dims = [2] * N
     targets = _targets_to_list(targets, oper=oper, N=N)
-    _check_qubits_oper(oper, dims=dims, targets=targets)
+    _check_oper_dims(oper, dims=dims, targets=targets)
 
     # Call expand_operator for all cyclic permutation of the targets.
     if cyclic_permutation:
+        warnings.warn(
+            "cyclic_permutation is deprecated, "
+            "please use loop through different targets manually.",
+            DeprecationWarning,
+        )
         oper_list = []
         for i in range(N):
             new_targets = np.mod(np.array(targets) + i, N)
@@ -1453,7 +1296,7 @@ def expand_operator(oper, N, targets, dims=None, cyclic_permutation=False):
             )
         return oper_list
 
-    # Generate the correct order for qubits permutation,
+    # Generate the correct order for permutation,
     # eg. if N = 5, targets = [3,0], the order is [1,2,3,0,4].
     # If the operator is cnot,
     # this order means that the 3rd qubit controls the 0th qubit.
diff --git a/src/qutip_qip/operations/measurement.py b/src/qutip_qip/operations/measurement.py
index 140e7b23..8ce28ae6 100644
--- a/src/qutip_qip/operations/measurement.py
+++ b/src/qutip_qip/operations/measurement.py
@@ -86,7 +86,7 @@ def measurement_comp_basis(self, state):
             raise ValueError("target is not valid")
 
         measurement_ops = [
-            expand_operator(op, N=n, targets=self.targets)
+            expand_operator(op, dims=[2] * n, targets=self.targets)
             for op in measurement_ops
         ]
         return measurement_statistics(state, measurement_ops)
diff --git a/src/qutip_qip/pulse.py b/src/qutip_qip/pulse.py
index 55fd9e5a..e9e9ea7e 100644
--- a/src/qutip_qip/pulse.py
+++ b/src/qutip_qip/pulse.py
@@ -53,7 +53,7 @@ def get_qobj(self, dims):
         else:
             qobj = self.qobj
             targets = self.targets
-        return expand_operator(qobj, len(dims), targets, dims)
+        return expand_operator(qobj, dims=dims, targets=targets)
 
     def _get_qobjevo_helper(self, spline_kind, dims):
         """


Test patch:
diff --git a/tests/decomposition_functions/test_single_qubit_gate_decompositions.py b/tests/decomposition_functions/test_single_qubit_gate_decompositions.py
index 333938c2..da72b8d5 100644
--- a/tests/decomposition_functions/test_single_qubit_gate_decompositions.py
+++ b/tests/decomposition_functions/test_single_qubit_gate_decompositions.py
@@ -15,11 +15,11 @@
 from qutip_qip.operations.gates import snot, sqrtnot
 
 # Fidelity closer to 1 means the two states are similar to each other
-H = snot(1, 0)
+H = snot()
 sigmax = sigmax()
 sigmay = sigmay()
 sigmaz = sigmaz()
-SQRTNOT = sqrtnot(N=1, target=0)
+SQRTNOT = sqrtnot()
 T = Qobj([[1, 0], [0, cmath.rect(1, np.pi / 4)]])
 S = Qobj([[1, 0], [0, 1j]])
 target = 0
diff --git a/tests/test_gates.py b/tests/test_gates.py
index 41603f2c..e1f7284f 100644
--- a/tests/test_gates.py
+++ b/tests/test_gates.py
@@ -10,7 +10,7 @@
 from qutip_qip.operations import (
     X, Y, Z, RX, RY, RZ, H, SQRTNOT, S, T, QASMU, CNOT, CPHASE, ISWAP, SWAP,
     CZ, SQRTSWAP, SQRTISWAP, SWAPALPHA, SWAPALPHA, MS, TOFFOLI, FREDKIN,
-    BERKELEY)
+    BERKELEY, expand_operator)
 
 
 def _permutation_id(permutation):
@@ -32,7 +32,7 @@ def _make_random_three_qubit_gate():
     def gate(N=None, controls=None, target=None):
         if N is None:
             return operation
-        return gates.gate_expand_3toN(operation, N, controls, target)
+        return expand_operator(operation, dims=[2]*N, targets=controls + [target])
     return gate
 
 
@@ -104,9 +104,8 @@ def test_swap(self):
                              tuple(itertools.permutations([0, 1, 2])),
                              ids=_permutation_id)
     def test_toffoli(self, permutation):
-        test = gates.toffoli(N=3,
-                             controls=permutation[:2],
-                             target=permutation[2])
+        test = expand_operator(
+            gates.toffoli(), dims=[2] * 3, targets=permutation)
         base = (qutip.tensor(1 - qutip.basis([2, 2], [1, 1]).proj(),
                              qutip.qeye(2))
                 + qutip.tensor(qutip.basis([2, 2], [1, 1]).proj(),
@@ -178,6 +177,7 @@ def test_gate_normalises_pauli_group(self, gate):
         assert len(pauli_gates) == 0
 
 
+@pytest.mark.filterwarnings('ignore::DeprecationWarning')
 class TestGateExpansion:
     """
     Test that gates act correctly when supplied with controls and targets, i.e.
@@ -267,8 +267,10 @@ class Test_expand_operator:
         ids=_permutation_id)
     def test_permutation_without_expansion(self, permutation):
         base = qutip.tensor([qutip.rand_unitary(2) for _ in permutation])
-        test = gates.expand_operator(base,
-                                     N=len(permutation), targets=permutation)
+        test = gates.expand_operator(
+            base,
+            dims=[2] * len(permutation),
+            targets=permutation)
         expected = base.permute(_apply_permutation(permutation))
         np.testing.assert_allclose(test.full(), expected.full(), atol=1e-15)
 
@@ -283,12 +285,14 @@ def test_general_qubit_expansion(self, n_targets):
         for targets in itertools.permutations(range(n_qubits), n_targets):
             expected = _tensor_with_entanglement([qutip.qeye(2)] * n_qubits,
                                                  operation, targets)
-            test = gates.expand_operator(operation, n_qubits, targets)
+            test = gates.expand_operator(
+                operation, dims=[2] * n_qubits, targets=targets)
             np.testing.assert_allclose(test.full(), expected.full(),
                                        atol=1e-15)
 
     def test_cnot_explicit(self):
-        test = gates.expand_operator(gates.cnot(), 3, [2, 0]).full()
+        test = gates.expand_operator(
+            gates.cnot(), dims=[2] * 3, targets=[2, 0]).full()
         expected = np.array([[1, 0, 0, 0, 0, 0, 0, 0],
                              [0, 0, 0, 0, 0, 1, 0, 0],
                              [0, 0, 1, 0, 0, 0, 0, 0],
@@ -316,17 +320,6 @@ def test_hadamard_explicit(self):
         expected = expected / np.sqrt(8)
         np.testing.assert_allclose(test, expected)
 
-    def test_cyclic_permutation(self):
-        operators = [qutip.sigmax(), qutip.sigmaz()]
-        test = gates.expand_operator(qutip.tensor(*operators), N=3,
-                                     targets=[0, 1], cyclic_permutation=True)
-        base_expected = qutip.tensor(*operators, qutip.qeye(2))
-        expected = [base_expected.permute(x)
-                    for x in [[0, 1, 2], [1, 2, 0], [2, 0, 1]]]
-        assert len(expected) == len(test)
-        for element in expected:
-            assert element in test
-
     @pytest.mark.parametrize('dimensions', [
         pytest.param([3, 4, 5], id="standard"),
         pytest.param([3, 3, 4, 4, 2], id="standard"),
@@ -340,8 +333,8 @@ def test_non_qubit_systems(self, dimensions):
                          for n, dimension in enumerate(dimensions)]
             expected = qutip.tensor(*operators)
             base_test = qutip.tensor(*[operators[x] for x in targets])
-            test = gates.expand_operator(base_test, N=n_qubits,
-                                         targets=targets, dims=dimensions)
+            test = gates.expand_operator(base_test, dims=dimensions,
+                                         targets=targets)
             assert test.dims == expected.dims
             np.testing.assert_allclose(test.full(), expected.full())