old working code for two level arrays

src/qutip_qip/decompose/_utility.py

@@ -29,3 +29,64 @@ def check_gate(gate, num_qubits):
         raise ValueError("Input is not unitary.")
     if gate.dims != [[2] * num_qubits] * 2:
         raise ValueError(f"Input is not a unitary on {num_qubits} qubits.")
+
+
+def _binary_sequence(num_qubits):
+    """ Defines the binary sequence list for basis vectors of a n-qubit gate.
+    The string at index `i` is also the row/column index for a basis vector in
+    a numpy array.
+    """
+    old_sequence = ['0', '1']
+    full_binary_sequence = []
+
+    if num_qubits == 1:
+        full_binary_sequence = old_sequence
+    else:
+        for x in range(num_qubits-1):
+            full_binary_sequence = []
+            zero_append_sequence = ['0' + x for x in old_sequence]
+            full_binary_sequence.extend(zero_append_sequence)
+            one_append_sequence = ['1' + x for x in old_sequence]
+            full_binary_sequence.extend(one_append_sequence)
+            old_sequence = full_binary_sequence
+
+    return(full_binary_sequence)
+
+
+def _gray_code_sequence(num_qubits, output_form=None):
+    """ Finds the sequence of gray codes for basis vectors by using logical
+    XOR operator of Python.
+    For print( _gray_code_sequence(2)), the output is [0, 1, 3, 2] i.e. in
+    terms of the binary sequence, the output is ['00', '01', '11', '10'].
+    https://docs.python.org/3/library/operator.html#operator.xor'
+
+    Parameters
+    ----------
+    num_qubits
+        Number of qubits in the circuit
+    output_form : :"index_values" or None
+        The format of output list. If a string "index_values" is provided then
+    the function's output is in terms of array indices of the binary sequence.
+    The default is a list of binary strings.
+
+    Returns
+    --------
+    list
+        List of the gray code sequence in terms of array indices or binary
+        sequence positions.
+    """
+    gray_code_sequence_as_array_indices = []
+
+    for x in range(2**num_qubits):
+        gray_code_at_x = x ^ x // 2  # floor operator to shift bits by 1
+        # when the shift is done, the new spot is filled with a new value.
+        gray_code_sequence_as_array_indices.append(gray_code_at_x)
+        if output_form == "index_values":
+            output = gray_code_sequence_as_array_indices
+        else:
+            gray_code_as_binary = []
+            binary_sequence_list = _binary_sequence(num_qubits)
+            for i in gray_code_sequence_as_array_indices:
+                gray_code_as_binary.append(binary_sequence_list[i])
+            output = gray_code_as_binary
+    return(output)

src/qutip_qip/decompose/decompose_general_qubit_gate.py

@@ -0,0 +1,94 @@
+import numpy as np
+import cmath
+from qutip import Qobj
+from qutip_qip.decompose._utility import (
+    check_gate,
+)
+
+
+def _decompose_to_two_level_arrays(input_gate, num_qubits, expand=True):
+    """Decompose a general qubit gate to two-level arrays.
+
+    Parameters
+    -----------
+    input_gate : :class:`qutip.Qobj`
+        The gate matrix to be decomposed.
+    num_qubits : int
+        Number of qubits being acted upon by the input_gate
+    expand : True
+        Default parameter to return the output as full two-level Qobj. If
+    `expand = False` then the function returns a tuple of index information
+    and a 2 x 2 Qobj for each gate. The Qobj are returned in reversed order.
+    """
+    check_gate(input_gate, num_qubits)
+    input_array = input_gate.full()
+
+    # Calculate the two level numpy arrays
+    array_list = []
+    index_list = []
+    for i in range(2**num_qubits):
+        for j in range(i+1, 2**num_qubits):
+            new_index = [i, j]
+            index_list.append(new_index)
+
+    for i in range(len(index_list)-1):
+        index_1, index_2 = index_list[i]
+
+        # Values of single qubit U forming the two level unitary
+        a = input_array[index_1][index_1]
+        a_star = np.conj(a)
+        b = input_array[index_2][index_1]
+        b_star = np.conj(b)
+        norm_constant = cmath.sqrt(
+            np.absolute(a*a_star)+np.absolute(b*b_star))
+
+        # Create identity array and then replace with above values for
+        # index_1 and index_2
+        U_two_level = np.identity(2**num_qubits, dtype=complex)
+        U_two_level[index_1][index_1] = a_star/norm_constant
+        U_two_level[index_2][index_1] = b/norm_constant
+        U_two_level[index_1][index_2] = b_star/norm_constant
+        U_two_level[index_2][index_2] = -a/norm_constant
+
+        # Change input by multiplying by above two-level
+        input_array = np.dot(U_two_level, input_array)
+
+        # U dagger to calculate the gates
+        U__two_level_dagger = np.transpose(np.conjugate(U_two_level))
+        array_list.append(U__two_level_dagger)
+
+    # for U6 - multiply input array by U5 and take dagger
+    U_last_dagger = input_array
+    array_list.append(U_last_dagger)
+
+    if expand is True:
+        array_list_with_qobj = []
+        for i in reversed(range(len(index_list))):
+            U_two_level_array = array_list[i]
+            array_list_with_qobj.append(Qobj(
+                U_two_level_array, dims=[[2] * num_qubits] * 2))
+        return(array_list_with_qobj)
+    else:
+        compact_U_information = []
+        for i in reversed(range(len(index_list))):
+            U_non_trivial = np.full([2, 2], None, dtype=complex)
+            index_info = []
+            U_index_together = []
+
+            # create index list
+            index_1, index_2 = index_list[i]
+            index_info = [index_1, index_2]
+            U_index_together.append(index_info)
+
+            # create 2 x 2 arrays
+            U_two_level = array_list[i]
+            U_non_trivial[0][0] = U_two_level[index_1][index_1]
+            U_non_trivial[1][0] = U_two_level[index_2][index_1]
+            U_non_trivial[0][1] = U_two_level[index_1][index_2]
+            U_non_trivial[1][1] = U_two_level[index_2][index_2]
+            U_index_together.append(
+                Qobj(U_non_trivial, dims=[[2] * 1] * 2))
+
+            compact_U_information.append(U_index_together)
+
+        return(compact_U_information)

tests/decomposition_functions/test_general_decomposition.py

@@ -0,0 +1,34 @@
+
+
+import numpy as np
+import pytest
+
+from qutip import (
+    Qobj, average_gate_fidelity, rand_unitary
+)
+
+
+from qutip_qip.decompose.decompose_general_qubit_gate import (
+    _decompose_to_two_level_arrays)
+
+
+@pytest.mark.parametrize("num_qubits", [2, 3, 4, 5, 6])
+def test_two_level_full_output(num_qubits):
+    """ Check if product of full two level array output is equal to the input.
+    """
+    input_gate = rand_unitary(2**num_qubits, dims=[[2] * num_qubits] * 2)
+    array_decompose = _decompose_to_two_level_arrays(
+        input_gate, num_qubits, expand=True)
+
+    product_of_U = array_decompose[-1]
+
+    for i in reversed(range(len(array_decompose)-1)):
+        product_of_U = product_of_U
+        product_of_U_calculated = np.dot(product_of_U, array_decompose[i])
+        product_of_U = product_of_U_calculated
+
+    product_of_U = Qobj(product_of_U, dims=[[2] * num_qubits] * 2)
+    fidelity_of_input_output = average_gate_fidelity(
+        product_of_U, input_gate
+    )
+    assert np.isclose(fidelity_of_input_output, 1.0)