<jupyter_start><jupyter_text>1. Get OpenFermion representation of Hamiltonian<jupyter_code>from quchem.Misc_functions.conversion_scripts import Get_Openfermion_Hamiltonian

openFermion_H = Get_Openfermion_Hamiltonian(Hamilt_dictionary)
openFermion_H<jupyter_output><empty_output><jupyter_text>2. Get cliques defined by commutativity<jupyter_code>from quchem.Unitary_Partitioning.Graph import Clique_cover_Hamiltonian

commutativity_flag = 'AC' ## <- defines relationship between sets!!!
Graph_colouring_strategy='largest_first'


anti_commuting_sets = Clique_cover_Hamiltonian(openFermion_H, 
                                                     N_qubits, 
                                                     commutativity_flag, 
                                                     Graph_colouring_strategy)
anti_commuting_sets<jupyter_output><empty_output><jupyter_text>3. Example of X_sk operator  lexicographica order(maximises circuit reductions)<jupyter_code>key_larg, largest_AC_set = max(anti_commuting_sets.items(), key=lambda x:len(x[1])) # largest nonCon part found by dfs alg
largest_AC_set
from quchem.Misc_functions.Misc_functions import choose_Pn_index

Ps_index = choose_Pn_index(largest_AC_set)
B_Ps = largest_AC_set.pop(Ps_index) # REMOVE BETA_S_P_S from largest_AC_set
from quchem.Misc_functions.Misc_functions import lexicographical_sort_BASIS_MATCH
BASIS_MATCH_AC_set = lexicographical_sort_BASIS_MATCH(largest_AC_set)
BASIS_MATCH_AC_set.append(B_Ps) # add Ps back in last index!
BASIS_MATCH_AC_set
from quchem.Misc_functions.Misc_functions import lexicographical_sort_LADDER_CNOT_cancel
CNOT_LADDER_MATCH_AC_set = lexicographical_sort_LADDER_CNOT_cancel(largest_AC_set)
CNOT_LADDER_MATCH_AC_set.append(B_Ps) # add Ps back in last index!
CNOT_LADDER_MATCH_AC_set
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Xsk_op_list
from openfermion import QubitOperator

Ps_index=-1 # last term!
check_reduction_lin_alg = True
X_sk_theta_sk_list, full_normalised_set, Ps, gamma_l = Get_Xsk_op_list(BASIS_MATCH_AC_set,
                                                                        Ps_index, 
                                                                        N_qubits,
                                                                        check_reduction=check_reduction_lin_alg,
                                                                    )


angle_list=[]
op_list =[]
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]
    
    
    angle_list.append(theta_sk / 2 * const_X_sk)
    
    op_list.append(QubitOperator(pauliword_X_sk, -1j))
    
op_list
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import Optimized_LADDER_circuit

c_opt = Optimized_LADDER_circuit(op_list,
                                 angle_list, 
                                 check_reduction=True)
c_opt
from quchem.Misc_functions.Misc_functions import count_circuit_gates
print(count_circuit_gates(c_opt))
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import rewrite_Y_basis_change, ladder_or_stairs_circuit_array, circuit_array_to_circuit
unop_arr = ladder_or_stairs_circuit_array(op_list, ladder=True)
unop_arr = rewrite_Y_basis_change(unop_arr)
c_unop = circuit_array_to_circuit(unop_arr,
                                  angle_list,
                                 op_list)
c_unop
print(count_circuit_gates(c_unop))
IBM_opt, phase = optimized_cirq_circuit_IBM_compiler(c_unop, check_optimization=True)
print(count_circuit_gates(IBM_opt))
IBM_opt
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Xsk_op_list
from openfermion import QubitOperator

Ps_index=-1 # last term!
check_reduction_lin_alg = True
X_sk_theta_sk_list, full_normalised_set, Ps, gamma_l = Get_Xsk_op_list(CNOT_LADDER_MATCH_AC_set, # < --HERE
                                                                        Ps_index, 
                                                                        N_qubits,
                                                                        check_reduction=check_reduction_lin_alg,
                                                                    )


angle_list=[]
op_list =[]
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]
    
    
    angle_list.append(theta_sk / 2 * const_X_sk)
    
    op_list.append(QubitOperator(pauliword_X_sk, -1j))
    
    
c_opt_LADDER = Optimized_LADDER_circuit(op_list,
                                 angle_list, 
                                 check_reduction=True)
c_opt_LADDER
print(count_circuit_gates(c_opt))
op_list =[QubitOperator('Y0 Y1 Y2', -1j),
         QubitOperator('X0 X1 Y2', -1j)]
angle_list = [2, 3]
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import Optimized_LADDER_circuit
c_opt = Optimized_LADDER_circuit(op_list,
                                 angle_list, 
                                 check_reduction=True)
c_opt
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import Optimized_STAIRS_circuit
c_opt_STAIRS = Optimized_STAIRS_circuit(op_list,
                                 angle_list, 
                                 check_reduction=True)
c_opt_STAIRS
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import ladder_or_stairs_circuit_array, rewrite_Y_basis_change, circuit_array_to_circuit
non_opt_arr = ladder_or_stairs_circuit_array(op_list)
non_opt_arr = rewrite_Y_basis_change(non_opt_arr)
c_UNoptimized = circuit_array_to_circuit(non_opt_arr, angle_list, op_list)
print(np.allclose(c_UNoptimized.unitary(), c_opt.unitary()))
c_UNoptimized
from quchem.Misc_functions.Misc_functions import count_circuit_gates
print(count_circuit_gates(c_opt))
print(count_circuit_gates(c_UNoptimized))
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import ladder_or_stairs_circuit_array
arr1 = ladder_or_stairs_circuit_array(op_list)
arr1
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import cancel_common_change_basis
arr2 = cancel_common_change_basis(arr1)
arr2
Sd @ cirq.X._unitary_()
S @ H @ S.conj().T == -1J* S@S@H
X @ S.conj().T ==  -1j*(S.conj().T @ X)
Sd @ H @ Sd @ H 
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import full_exponentiated_PauliWord_circuit
from quchem.Qcircuit.Hamiltonian_term_measurement_functions import change_pauliword_to_Z_basis_then_measure
import cirq

circuit_list = []
for X_sk_Op, theta_sk in zip(op_list, angle_list):
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]

    full_exp_circ_obj = full_exponentiated_PauliWord_circuit(QubitOperator(pauliword_X_sk, -1j),
                                                             theta_sk / 2 * const_X_sk)

    circuit = cirq.Circuit(
        cirq.decompose_once((full_exp_circ_obj(*cirq.LineQubit.range(full_exp_circ_obj.num_qubits())))))

    circuit_list.append(circuit)

full_RS_circuit = cirq.Circuit(circuit_list)
full_RS_circuit
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import full_exponentiated_PauliWord_circuit
from quchem.Qcircuit.Hamiltonian_term_measurement_functions import change_pauliword_to_Z_basis_then_measure
import cirq

circuit_list = []
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]

    full_exp_circ_obj = full_exponentiated_PauliWord_circuit(QubitOperator(pauliword_X_sk, -1j),
                                                             theta_sk / 2 * const_X_sk)

    circuit = cirq.Circuit(
        cirq.decompose_once((full_exp_circ_obj(*cirq.LineQubit.range(full_exp_circ_obj.num_qubits())))))

    circuit_list.append(circuit)

full_RS_circuit = cirq.Circuit(circuit_list)
full_RS_circuit
np.allclose(full_RS_circuit.unitary(), c_opt.unitary())
cirq.testing.assert_allclose_up_to_global_phase(
    full_RS_circuit.unitary(),
    c_opt.unitary(),
    atol=1e-8,
    verbose=False)
dfsd
CNOT_LADDER_MATCH_AC_set
from quchem.Misc_functions.Misc_functions import count_circuit_gates
count_circuit_gates(c_opt)
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import ladder_or_stairs_circuit_array, circuit_array_to_circuit
non_opt_arr = ladder_or_stairs_circuit_array(CNOT_LADDER_MATCH_AC_set)
c_UNoptimized = circuit_array_to_circuit(non_opt_arr, [i for i in range(len(CNOT_LADDER_MATCH_AC_set))])
c_UNoptimized
opt_circuit = c_opt
angle_list=[i for i in range(len(CNOT_LADDER_MATCH_AC_set))]
P_Word_list=CNOT_LADDER_MATCH_AC_set
angle_list = [-1*angle_list[ind] if list(op.terms.values())[0]<0 else 1*angle_list[ind] for ind, op in enumerate(P_Word_list)]

N_qubits = max([lineQ.x for lineQ in list(opt_circuit.all_qubits())]) + 1
mat_list = [expm(qubit_operator_sparse(op, n_qubits = N_qubits)*angle_list[ind]/2).todense() for ind, op in enumerate(P_Word_list)]
lin_alg_mat = reduce(np.dot, mat_list[::-1])

if not np.allclose(lin_alg_mat, opt_circuit.unitary()):
    raise ValueError('circuit reduction incorrect!') 
np.allclose(c_UNoptimized.unitary(), c_opt.unitary())
count_circuit_gates(c_UNoptimized)
# from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import Optimized_STAIRS_circuit
# c_opt = Optimized_STAIRS_circuit(CNOT_LADDER_MATCH_AC_set, [i for i in range(len(CNOT_LADDER_MATCH_AC_set))])
# count_circuit_gates(c_opt)
# c_opt
count_circuit_gates
CNOT_LADDER_MATCH_AC_set
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import *

arr1 = ladder_or_stairs_circuit_array(CNOT_LADDER_MATCH_AC_set)
arr2 = NEW(arr1)
# arr3 = cancel_S_gates_through_controls(arr2)
# arr4 = test(arr3)
circuit_array_to_circuit(arr2, [i for i in range(len(CNOT_LADDER_MATCH_AC_set))])
def NEW(circuit_array):
    """
    Given a circuit array function will cancel any common change of basis gates
    """
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig in enumerate(qubit_operations):
            if sig in ['I', 'R' 'C']:
                continue
            for new_ind in range(ind+1, circuit_array.shape[1]):
                sig_new = qubit_operations[new_ind]
                if sig_new in ['C', 'X']:
                    continue
                elif sig == sig_new:
                    qubit_operations[ind]='I'
                    qubit_operations[new_ind]='I'
                    break
                elif sig_new == 'I':
                    continue
                elif sig!=sig_new:
                    break
    return circuit_array
arr2 = NEW(arr1)
for circuit_slice in arr2.T:
    non_I_indices = np.where(circuit_slice!='I')[0]

    if 'C' in circuit_slice:
        if 'X' not in circuit_slice:
#             print(circuit_slice)
            pass
    else:
        for qNo in non_I_indices:
            op = circuit_slice[qNo]
            if op not in ['H', 'R', 'S', 'Sdag']:
                print(op)
circuit_array_to_circuit(arr2, [i for i in range(len(CNOT_LADDER_MATCH_AC_set))])
def test(circuit_array):
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]-1):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig_control in enumerate(qubit_operations):
            if sig_control != 'C':
                continue
            print()
#             print(circuit_array[:,ind])
            print(qubit_row)
            next_row_ind = np.where(circuit_array[:,ind]=='X')[0][0]
            if next_row_ind.size>0:
                print(next_row_ind)
                qubit_operations_NEXT = circuit_array[next_row_ind, :]

                sig_target = qubit_operations_NEXT[ind]
                if sig_target != 'X':
                    continue

                for new_ind in range(ind+1, circuit_array.shape[1]):
                    sig_new_CONTROL = qubit_operations[new_ind]
                    sig_new_Target = qubit_operations_NEXT[new_ind]

                    if (sig_new_CONTROL in ['I', 'S', 'Sdag'] and sig_new_Target in ['I']):
                        continue
                    elif sig_new_CONTROL !='C':
                        break
                    elif sig_new_Target !='X':
                        break
                    else:
                        qubit_operations[ind]='I'
                        qubit_operations_NEXT[ind]='I'
                        qubit_operations[new_ind]='I'
                        qubit_operations_NEXT[new_ind]='I'
                        break
    return circuit_array
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Xsk_op_list
from openfermion import QubitOperator

Ps_index=-1 # last term!
check_reduction_lin_alg = True
X_sk_theta_sk_list, full_normalised_set, Ps, gamma_l = Get_Xsk_op_list(BASIS_MATCH_AC_set,
                                                                        Ps_index, 
                                                                        N_qubits,
                                                                        check_reduction=check_reduction_lin_alg,
                                                                    )


angle_list=[]
op_list =[]
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]
    
    
    angle_list.append(theta_sk / 2 * const_X_sk)
    
    op_list.append(QubitOperator(pauliword_X_sk, -1j))
    
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_LADDER
BASIS_MATCH_AC_circuit_LADDER_form = Circuit_from_list_PauliWords_LADDER(op_list, angle_list, check_reduction=True)
BASIS_MATCH_AC_circuit_LADDER_form
import cirq
from cirq.circuits import InsertStrategy

closest_together = cirq.Circuit(list(BASIS_MATCH_AC_circuit_LADDER_form.all_operations()), strategy=InsertStrategy.EARLIEST)
closest_together
# closest_together.moments
from quchem.Misc_functions.Misc_functions import count_circuit_gates
count_circuit_gates(BASIS_MATCH_AC_circuit_LADDER_form)
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_STAIRS
BASIS_MATCH_AC_circuit_STAIRS_form = Circuit_from_list_PauliWords_STAIRS(op_list, angle_list, check_reduction=True)
count_circuit_gates(BASIS_MATCH_AC_circuit_STAIRS_form)
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Xsk_op_list
from openfermion import QubitOperator

Ps_index=-1 # last term!
check_reduction_lin_alg = True
X_sk_theta_sk_list, full_normalised_set, Ps, gamma_l = Get_Xsk_op_list(CNOT_LADDER_MATCH_AC_set,
                                                                        Ps_index, 
                                                                        N_qubits,
                                                                        check_reduction=check_reduction_lin_alg,
                                                                    )


angle_list=[]
op_list =[]
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk = list(X_sk_Op.terms.keys())[0]
    const_X_sk = list(X_sk_Op.terms.values())[0]
    
    
    angle_list.append(theta_sk / 2 * const_X_sk)
    
    op_list.append(QubitOperator(pauliword_X_sk, -1j))
    
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_LADDER
CNOT_LADDER_MATCH_AC_circuit_LADDER_form = Circuit_from_list_PauliWords_LADDER(op_list, angle_list, check_reduction=True)
CNOT_LADDER_MATCH_AC_circuit_LADDER_form
count_circuit_gates(CNOT_LADDER_MATCH_AC_circuit_LADDER_form)
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_STAIRS
CNOT_LADDER_MATCH_AC_circuit_STAIRS_form = Circuit_from_list_PauliWords_STAIRS(op_list, angle_list, check_reduction=True)
count_circuit_gates(CNOT_LADDER_MATCH_AC_circuit_STAIRS_form)
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_STAIRS
circuit_STAIRS_form = Circuit_from_list_PauliWords_STAIRS(op_list, angle_list, check_reduction=True)
circuit_STAIRS_form
print(len(list(circuit_LADDER_form.all_operations())))
print(len(list(circuit_STAIRS_form.all_operations())))

dir(list(circuit_LADDER_form.all_operations())[0])
# type(list(circuit_LADDER_form.all_operations())[1].gate)
count_circuit_gates(circuit_LADDER_form)
circuit_LADDER_form
fadsdfsd
from openfermion.linalg import qubit_operator_sparse
from scipy.sparse.linalg import expm
R_sk_list = []
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk_MATRIX = qubit_operator_sparse(QubitOperator(list(X_sk_Op.terms.keys())[0], -1j),
                                                  n_qubits=N_qubits)
    const_X_sk = list(X_sk_Op.terms.values())[0]
    R_sk_list.append(expm(pauliword_X_sk_MATRIX * theta_sk / 2 * const_X_sk))

R_S_matrix = reduce(np.dot, R_sk_list[::-1])  # <- note reverse order!
R_S_matrix
np.allclose(circ.unitary(), R_S_matrix.todense())
 R_S_matrix.todense()
circ.unitary()
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import exp_pauliword_reduced_QC_STAIRS
def exp_pauliword_reduced_QC_STAIRS(P_A, theta_A, pre_P_A_term=None, post_P_A_term=None):
    """

    ## example
    from openfermion import QubitOperator
    import numpy as np

    op1 = QubitOperator('X0 X1 X2', -1j)
    op2 = QubitOperator('Z0 X1 X2', -1j)
    op3 = QubitOperator('X0 Y1 X2', -1j)
    
    OP_list = [op1, op2, op3]
    theta_A= np.pi

    circuit_Test = Pair_exp_pauliword_STAIRS(OP_list[1], 
                                         theta_A,
                                         pre_P_A_term=None,
                                         post_P_A_term=None)
    print(circuit_Test)
                                   
        0: ───────@─────────────────────@─────────
                  │                     │
        1: ───H───┼───@────────────@────┼─────H───
                  │   │            │    │
        2: ───H───X───X───Rz(2π)───X────X─────H───



    circuit_Test_REDUCED = Pair_exp_pauliword_STAIRS(OP_list[1], 
                                         theta_A,
                                         pre_P_A_term=OP_list[0],
                                         post_P_A_term=OP_list[2])
    print(circuit_Test_REDUCED)
            
        0: ───@─────────────────@─────
              │                 │
        1: ───┼────────────@────┼───H─
              │            │    │
        2: ───X───Rz(2π)───X────X─────
                              
                                   
    """

    input_list = [P_A, pre_P_A_term, post_P_A_term]
    active_list = [term for term in input_list if term is not None]
    
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, active_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = cirq.LineQubit.range(max_qubits)
    
    
    PauliStrs_A, coeff_A = tuple(*P_A.terms.items())
    qNos_A, sigmas_A = zip(*list(PauliStrs_A))

    circuit = cirq.Circuit()

    if pre_P_A_term is None:
        ## no cancellation on LHS for P_A

        # change of basis
        for qNo, sigma_str in zip(qNos_A, sigmas_A):
            gate = change_basis_initial(sigma_str, qubits[qNo])
            if gate: circuit.append(gate) 

        ## entangle initial for P_A
        target_q = qNos_A[-1]
        for ind, qNo in enumerate(qNos_A[:-1]):
            circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))

    else:
        ## check for cancellations on LHS for P_A
        PauliStrs_pre, coeff_pre = tuple(*pre_P_A_term.terms.items())
        qNos_pre, sigmas_pre = zip(*list(PauliStrs_pre))

        common_qubits = np.intersect1d(qNos_pre, qNos_A)
        
        reduction_possible_dict = {common_q: sigmas_pre[qNos_pre.index(common_q)]==sigmas_A[qNos_A.index(common_q)]
                                                            for common_q in common_qubits}
        
        for qNo, sigma_str in zip(qNos_A, sigmas_A):
            if reduction_possible_dict.get(qNo, False):
                pass
            else:
                gate = change_basis_initial(sigma_str, qubits[qNo])
                if gate: circuit.append(gate) 
                    
        ## entangle initial for P_A
        target_q = qNos_A[-1]
        
        pre_CNOT_red= reduction_possible_dict.get(target_q, True)
        
        if pre_CNOT_red:
            for ind, qNo in enumerate(qNos_A[:-1]):
                if reduction_possible_dict[qNo]:
                    continue
                else:
                    circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))
        else:
            for ind, qNo in enumerate(qNos_A[:-1]):
                circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))


    ### PERFORM ROTATION
    if coeff_A.imag==0:
        raise ValueError('not valid qubit rotation')
    elif coeff_A.imag < 0:
        circuit.append(cirq.rz(2 * theta_A * np.abs(coeff_A.imag)).on(qubits[qNos_A[-1]]))
    else:
        circuit.append(cirq.rz(-2 * theta_A * np.abs(coeff_A.imag)).on(qubits[qNos_A[-1]]))


    #### RHS

    if post_P_A_term is None:
        ## no cancellation on RHS for P_A

        ## final initial for P_A
        target_q = qNos_A[-1]
        for ind, qNo in enumerate(qNos_A[:-1][::-1]):
            circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))
            
        # change of basis final
        for qNo, sigma_str in zip(qNos_A, sigmas_A):
            gate = change_basis_final(sigma_str, qubits[qNo])
            if gate: circuit.append(gate) 

    else:
        ## check for cancellations on RHS for P_A
        PauliStrs_post, coeff_post = tuple(*post_P_A_term.terms.items())
        qNos_post, sigmas_post = zip(*list(PauliStrs_post))

        common_qubits = np.intersect1d(qNos_post, qNos_A)

        reduction_possible_dict = {common_q: sigmas_post[qNos_post.index(common_q)]==sigmas_A[qNos_A.index(common_q)]
                                                            for common_q in common_qubits}
        
        
        ## entangle final for P_A
        target_q = qNos_A[-1]
        
        post_CNOT_red= reduction_possible_dict.get(target_q, True)
        
        if post_CNOT_red:
            for ind, qNo in enumerate(qNos_A[:-1][::-1]):
                if reduction_possible_dict[qNo]:
                    continue
                else:
                    circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))
        else:
            for ind, qNo in enumerate(qNos_A[:-1][::-1]):
                circuit.append(cirq.CNOT.on(qubits[qNo], qubits[target_q]))

        # change of basis
        for qNo, sigma_str in zip(qNos_A, sigmas_A):
            if reduction_possible_dict.get(qNo, False):
                pass
            else:
                gate = change_basis_final(sigma_str, qubits[qNo])
                if gate: circuit.append(gate) 
                    

    return circuit
test_op = QubitOperator('X0 Y1 Z2 X3', -1j)
test_PRE= QubitOperator('X0 Y1 Z2 X3', -1j)
test_POST= QubitOperator('X0 Z1 Z2 Z3', -1j)

theta_A = np.pi/2

circuit_form = exp_pauliword_reduced_QC_STAIRS(test_op, theta_A, pre_P_A_term=test_PRE, post_P_A_term=test_POST)

# lin_alg_mat = expm(qubit_operator_sparse(test_op, n_qubits = 4)*theta_A)

# print(np.allclose(lin_alg_mat.todense(),circuit_form.unitary()))
circuit_form
test_op = QubitOperator('X0 Y1 Z2 X3', -1j)
theta_A = np.pi/2

circuit_form = exp_pauliword_reduced_QC_STAIRS(test_op, theta_A, pre_P_A_term=None, post_P_A_term=None)

lin_alg_mat = expm(qubit_operator_sparse(test_op, n_qubits = 4)*theta_A)

print(np.allclose(lin_alg_mat.todense(),circuit_form.unitary()))
circuit_form
op1 = QubitOperator('X0 Y1 X2', -1j)
op2 = QubitOperator('X0 X1 X2', -1j)
op3 = QubitOperator('Y0 Y1 X2', -1j)

OP_list = [op1, op2, op3]
theta_A= np.pi 


exp_pauliword_reduced_QC_LADDER(op1 , theta_A, pre_P_A_term=op2, post_P_A_term=op3)
cirq.inverse(A)
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords
op1 = QubitOperator('X0 X1 X2', -1j)
op2 = QubitOperator('X0 X1 Z2', -1j)
op3 = QubitOperator('Y0 X1 Z2', -1j)

# op1 = QubitOperator('X0 X1 X2', -1j)
# op2 = QubitOperator('Z0 X1 X2', -1j)
# op3 = QubitOperator('X0 Y1 Y2', -1j)

OP_list = [op1, op2, op3]
angle_list = [i+1 for i in range(len(OP_list))]
circuit_form = Circuit_from_list_PauliWords(OP_list, angle_list, check_reduction=True)
circuit_form
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_STAIRS
circuit_form = Circuit_from_list_PauliWords_LADDER(OP_list, angle_list)
circuit_form


from quchem.Qcircuit.Ansatz_quantum_circuit_functions import Circuit_from_list_PauliWords_LADDER
circuit_form = Circuit_from_list_PauliWords_STAIRS(OP_list, angle_list)
circuit_form
OP_list
exp_pauliword_reduced_QC_STAIRS(op1, 1, pre_P_A_term=None, post_P_A_term=None)
exp_pauliword_reduced_QC_STAIRS(op2, 1, pre_P_A_term=None, post_P_A_term=None)
import cirq

N_qbits=3
mat_list = [ expm(qubit_operator_sparse(op, n_qubits = N_qbits)*angle_list[ind]).todense() for ind, op in enumerate(OP_list)]

qubits = list(cirq.LineQubit.range(N_qbits))
circuit_mat = cirq.Circuit()

for mat in mat_list:
    gate = cirq.MatrixGate(mat).on(*qubits)
    circuit_mat.append(gate)

lin_alg = reduce(np.dot, mat_list[::-1])

np.allclose(lin_alg, circuit_mat.unitary())
np.allclose(circuit_form.unitary(), circuit_mat.unitary())
np.around(circuit_mat.unitary(), 10) == np.around(circuit_form.unitary(), 10)
from quchem.Unitary_Partitioning.Seq_Rot_circuit_functions import Build_R_SeqRot_Q_circuit
S_index=0
check_reduction_lin_alg = True
check_circuit = True
AC_set = anti_commuting_sets[key_larg]

Q_circuit_Rsl, Psl, gammal = Build_R_SeqRot_Q_circuit(AC_set,
                                         S_index,
                                         N_qubits,
                                         check_reduction_lin_alg=check_reduction_lin_alg,
                                        atol=1e-8,
                                        rtol=1e-05,
                                        check_circuit = check_circuit)
Q_circuit_Rsl
Q_circuit_Rsl.unitary()
## compare with linear algebra
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Xsk_op_list

S_index=0
check_reduction = True

X_sk_theta_sk_list, normalised_FULL_set, Ps, gamma_l = Get_Xsk_op_list(AC_set,
                S_index,
                N_qubits,
                check_reduction=check_reduction,
                atol=1e-8,
                rtol=1e-05,)

X_sk_theta_sk_list
from quchem.Unitary_Partitioning.Unitary_partitioning_Seq_Rot import Get_Rsl_matrix
R_Sl_matrix = Get_Rsl_matrix(X_sk_theta_sk_list, N_qubits)
## checking Rsl circuit unitary is the same as lin alg R_sl matrix
import numpy as np
np.allclose(Q_circuit_Rsl.unitary(), R_Sl_matrix.todense())
Ps
gamma_l<jupyter_output><empty_output><jupyter_text>4. Full Seq Rot circuit<jupyter_code>import cirq

qubits = list(cirq.LineQubit.range(N_qubits))

ansatz = cirq.Circuit([cirq.X.on(q) for q in qubits])
ansatz
from quchem.Unitary_Partitioning.Seq_Rot_circuit_functions import Full_SeqRot_Rl_Circuit

S_index=0
check_reduction_lin_alg = True
AC_set = anti_commuting_sets[key_larg]

full_circuit, Ps, gamma_l = Full_SeqRot_Rl_Circuit(ansatz, AC_set, S_index,N_qubits,
                                                   check_reduction_lin_alg=check_reduction_lin_alg)
full_circuit
Ps<jupyter_output><empty_output><jupyter_text>5. Circuit experiments (TODO)- Ansatz- Rsl- Ps measurement<jupyter_code># TODO
Seq_Rot_VQE_Experiment_UP_circuit_lin_alg
# TODO
Seq_Rot_VQE_Experiment_UP_circuit_sampling
qbits = list(cirq.LineQubit.range(2))

test1 = cirq.Circuit(
cirq.CNOT(qbits[0], qbits[1]),
cirq.H(qbits[1]),
cirq.CNOT(qbits[0], qbits[1])
    
)
print(test1)

# test2 = cirq.Circuit(
# cirq.I(qbits[0]),
# cirq.H(qbits[1]),
# # cirq.CNOT(qbits[0], qbits[1])
# cirq.X.on(qbits[1])
# )
# test2

test2 = cirq.Circuit(
cirq.CNOT(qbits[0], qbits[1]),
cirq.ry(np.pi/2).on(qbits[1]),
cirq.X.on(qbits[0]),
cirq.CNOT(qbits[0], qbits[1]),
cirq.X.on(qbits[0])  
)
test2
np.around(test2.unitary(), 3) == np.around(test1.unitary(), 3)
cirq.Circuit(cirq.decompose(test2))
test1.unitary() == test2.unitary()
np.around(test1.unitary(), 3)
np.around(test2.unitary(), 3)
cirq.rz(np.pi)._unitary_()
cirq.H._unitary_()
cirq.X._unitary_() @ cirq.H._unitary_() @ cirq.X._unitary_() 
cirq.X._unitary_() @ cirq.H._unitary_()
cirq.X._unitary_() @ cirq.ry(np.pi)._unitary_() == cirq.ry(-np.pi)._unitary_()
cirq.YPowGate(exponent=0.5)._unitary_()
cirq.H._unitary_()
cirq.X._unitary_() @ cirq.ry(np.pi/2)._unitary_()
np.allclose(cirq.H._unitary_(), cirq.X._unitary_() @ cirq.ry(np.pi/2)._unitary_())
R_sk_list = []
for X_sk_Op, theta_sk in X_sk_theta_sk_list:
    pauliword_X_sk_MATRIX = qubit_operator_sparse(QubitOperator(list(X_sk_Op.terms.keys())[0], -1j),
                                                  n_qubits=N_qubits)
    const_X_sk = list(X_sk_Op.terms.values())[0]
    R_sk_list.append(expm(pauliword_X_sk_MATRIX * theta_sk / 2 * const_X_sk))

R_S_matrix = reduce(np.dot, R_sk_list[::-1])  # <- note reverse order!
R_S_matrix
qbits = list(cirq.LineQubit.range(4))
theta= 2

P_mat = qubit_operator_sparse(QubitOperator('X0 Y1 Z2 Z3', -1j)
                              ,n_qubits=4)
lin_mat = expm(P_mat * theta / 2)

cccc = cirq.Circuit(

cirq.H.on(qbits[0]),
cirq.rx(np.pi/2).on(qbits[1]),

cirq.CNOT(qbits[0], qbits[1]),
cirq.CNOT(qbits[1], qbits[2]),
cirq.CNOT(qbits[2], qbits[3]),
cirq.rz(theta).on(qbits[3]),
cirq.CNOT(qbits[2], qbits[3]),
cirq.CNOT(qbits[1], qbits[3]),
cirq.CNOT(qbits[0], qbits[3]),  

cirq.H.on(qbits[0]),
cirq.rx(-np.pi/2).on(qbits[1]),
)
cccc
cccc.unitary()
lin_mat.todense()

qbits = list(cirq.LineQubit.range(3))

test1 = cirq.Circuit(
cirq.CNOT(qbits[1], qbits[2]),
cirq.CNOT(qbits[0], qbits[1]),
cirq.CNOT(qbits[1], qbits[2])
    
)
print(test1)


test2 = cirq.Circuit(
cirq.CNOT(qbits[0], qbits[1]),
cirq.CNOT(qbits[0], qbits[2]),
cirq.I.on(qbits[1]),  
)
test2
test1.unitary() == test2.unitary()
qbits = list(cirq.LineQubit.range(3))

test1 = cirq.Circuit(
cirq.CNOT(qbits[0], qbits[1]),
cirq.H(qbits[0]),
cirq.CNOT(qbits[0], qbits[1])
    
)
print(test1)


test2 = cirq.Circuit(
cirq.S(qbits[0]),
cirq.H(qbits[1]),
cirq.CNOT(qbits[1], qbits[0]),
cirq.S(qbits[0])**-1,
cirq.S(qbits[1]),
cirq.H(qbits[0]),
cirq.H(qbits[1]),
)
test2
np.allclose(test2.unitary(), test1.unitary())
test1 = cirq.Circuit(
cirq.rx(-np.pi/2).on(qbits[0]),
)

print(test1)
print()

test2 = cirq.Circuit(
cirq.H(qbits[0]),
cirq.S(qbits[0]),
# cirq.S(qbits[0])**-1,
# cirq.H(qbits[0]),
)
print(test2)

np.allclose(test2.unitary(), test1.unitary())
test2.unitary()
 test1.unitary()
(cirq.S**-1)._unitary_()
np.exp(-1j*np.pi/2)
(cirq.S)._unitary_() @ cirq.H._unitary_() 
(cirq.X**-0.5)._unitary_() 
cirq.rx(np.pi/2)._unitary_() 
from scipy.linalg import svd
H_S_dag = (cirq.S)._unitary_() @ cirq.H._unitary_() 
svd(test2.unitary())
from scipy.linalg import svd
svd(cirq.H._unitary_())a
qbits = list(cirq.LineQubit.range(3))

test1 = cirq.Circuit(
cirq.CNOT(qbits[0], qbits[1]),
cirq.H(qbits[0]),
cirq.S(qbits[0]),
cirq.CNOT(qbits[0], qbits[1])
    
)
print(test1)
print()
print()

test2 = cirq.Circuit(
cirq.S(qbits[0]),
cirq.H(qbits[1]),
cirq.CNOT(qbits[1], qbits[0]),
cirq.S(qbits[0])**-1,
cirq.S(qbits[1]),
cirq.H(qbits[0]),
cirq.H(qbits[1]),
cirq.S(qbits[0]),
)
test2

print(test2)

np.allclose(test2.unitary(), test1.unitary())
(cirq.S._unitary_()*-1) @ H @ S
S_dag @ H @ S
H @ S @ H
zero = np.array([[0],[1]])
H = cirq.H._unitary_()
S = cirq.S._unitary_()
S_dag = (cirq.S**-1)._unitary_()

H_zero = H @ zero

S_H_zero = S @ H_zero
S_H_zero
zero = np.array([[0],[1]])
Rx = cirq.rx(-np.pi/2)._unitary_()

r_zero = Rx @ zero
r_zero *1j*1j*1j
from quchem.Misc_functions.Misc_functions import lexicographical_sort_BASIS_MATCH
BASIS_MATCH_AC_set = lexicographical_sort_BASIS_MATCH(largest_AC_set)
BASIS_MATCH_AC_set.append(B_Ps) # add Ps back in last index!
BASIS_MATCH_AC_set
from openfermion import QubitOperator
from openfermion.utils import count_qubits
from functools import reduce
import numpy as np

from openfermion.linalg import qubit_operator_sparse
from scipy.sparse.linalg import expm
BASIS_MATCH_AC_set

op_list = [
    
    QubitOperator('Y0 Z1 Z2 Y3 X4 X5', 1),
    QubitOperator('Y0 Z1 Z2 Y3 X4 Z5', 1),
    QubitOperator('Y0 Z1 Z2 Y3 Z4 X5', 1),
    QubitOperator('X0 X1 X2 X3 X4 X5', 1),
    QubitOperator('X0 X1    Y3 X4 X5', 1),
]
def Change_Basis_op_array(list_P_ops):
    """
    maximises adjacent single qubit pauli terms in Pauliword list (allowing best change of basis cancellation)
    """
    fullOp = reduce(lambda Op1, Op2: Op1+Op2, list_P_ops)
    max_qubits = count_qubits(fullOp)

    P_word_array = np.zeros((max_qubits, len(list_P_ops)), dtype=str)
    for ind, op in enumerate(list_P_ops):

        P_dict =  dict(tuple(*op.terms.keys())) 
        arr = np.array([P_dict.get(qNo, 'I') for qNo in range(max_qubits)], dtype=str)
          
        P_word_array[:, ind] = arr
    
    
    return P_word_array
    
    
Change_Basis_op_array(op_list)
# def change_basis_between_terms_LEFT(sigma_list):
#     """
    
#     ## example
    
#     T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
#     change, newT = change_basis_between_terms_RIGHT(T)
    
#     print(change)
#     >> ['Y', 'Z', 'Y', 'I', 'X', 'I', 'I', 'I', 'I']
    
#     print(newT)
#     >> ['Z', 'Z', 'Z', 'Z', 'Z', 'Z', 'I', 'I', 'Z']
    
#     look to left of each term in T list. 
#     Expect change of basis between terms (note first term cannot have change of basis as nothing to LHS)
#     new sigma list is paulioperators after change of basis
#     """
    
#     sigma_list = np.asarray(sigma_list)
#     new_sigma_list = sigma_list.copy()
    
#     change_basis= [sigma_list[0]] # cannot simplify first term
#     new_sigma_list=[]
#     for ind in range(sigma_list.shape[0]):
#         s = sigma_list[ind]
#         new_sigma_list.append('Z' if s!='I' else 'I')

#         if s == 'I':
#             change_basis.append('I')
#             continue

#         for s_next in sigma_list[ind+1:]:
#             if s_next == 'I':
#                 continue
#             elif s_next == s:
#                 change_basis.append('I')
#                 break  
#             else:
#                 change_basis.append(s_next)
#                 break
# #     new_sigma_list=['Z' if s!='I' else 'I'for s in sigma_list]
#     return change_basis, new_sigma_list

# # T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
# T = ['I', 'Y','Z', 'Y', 'Y']
# change, newT = change_basis_between_terms(T)
# print(T)
# print(change)
# # print(newT)
def change_basis_between_terms_LEFT(sigma_list):
    """
    
    ## example
    
    T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
    change, newT = change_basis_between_terms_RIGHT(T)
    
    print(change)
    >> ['Y', 'Z', 'Y', 'I', 'X', 'I', None, None, 'I']
    
    print(newT)
    >> ['Z', 'Z', 'Z', 'Z', 'Z', 'Z', 'I', 'I', 'Z']
    
    look to left of each term in T list. 
    Expect change of basis between terms (note first term cannot have change of basis as nothing to LHS)
    new sigma list is paulioperators after change of basis
    """
    
    sigma_list = np.asarray(sigma_list)
    new_sigma_list = sigma_list.copy()
    
    change_basis= [sigma_list[0] if sigma_list[0]!='I' else 'N'] # cannot simplify first term
    new_sigma_list=['Z' if sigma_list[0]!='I' else 'I']
    for ind in range(1,sigma_list.shape[0]):
        s = sigma_list[ind]
        new_sigma_list.append('Z' if s!='I' else 'I')

        if s == 'I':
            change_basis.append('N')
            continue

        for prev_ind in range(ind-1, -1, -1):
            s_prev = sigma_list[prev_ind]
            if (s_prev == 'I' and prev_ind!=0):
                continue
            elif (s_prev == 'I' and prev_ind==0):
                change_basis.append(s)
                break
            elif s_prev == s:
                change_basis.append('I')
                break  
            else:
                if (s=='X' and s_prev=='Y'):
                    # XY = H HS = S gate only
                    change_basis.append('YX')
                elif (s=='Y' and s_prev=='X'):
                    # YX = S_dagger H H = S_dagger gate only
                    change_basis.append('XY')
                else:
                    change_basis.append(s)
                break
                
    return change_basis, new_sigma_list

T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
# T = ['I', 'X','Z', 'Y', 'Y']
T=['Y', 'Z', 'Z', 'Y', 'X', 'X']
change, newT = change_basis_between_terms_LEFT(T)
# print(T)
print(change)
# print(newT)
def change_basis_between_terms_RIGHT(sigma_list):
    """
    
    ## example
    
    T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
    change, newT = change_basis_between_terms_RIGHT(T)
    
    print(change)
    >> ['Y', 'Z', 'I', 'Y', 'I', 'I', None, None, 'X']
    
    print(newT)
    >> ['Z', 'Z', 'Z', 'Z', 'Z', 'Z', 'I', 'I', 'Z']
    
    look to RIGHT of each term in T list. 
    Expect change of basis between terms (note LAST term cannot have change of basis as nothing to RHS)
    new sigma list is paulioperators after change of basis
    """
    sigma_list = np.asarray(sigma_list)
    new_sigma_list = sigma_list.copy()
    
    change_basis= [] 
    new_sigma_list=[]
    for ind in range(sigma_list.shape[0]-1):
        s = sigma_list[ind]
        new_sigma_list.append('Z' if s!='I' else 'I')
        if s == 'I':
            change_basis.append('N')
            continue

        for ind_next in range(ind+1, sigma_list.shape[0]):
            s_next = sigma_list[ind_next]
            if (s_next == 'I' and ind_next==(sigma_list.shape[0]-1)):
                change_basis.append(s)
            elif s_next == 'I':
                continue
            elif s_next == s:
                change_basis.append('I')
                break  
            else:
                if (s=='X' and s_next=='Y'):
                    # XY = H H S = S gate only
                    change_basis.append('XY')
#                     change_basis.append('N') #included in LEFT optimization
                elif (s=='Y' and s_next=='X'):
                    # YX = S_dagger H H = S_dagger gate only
                    change_basis.append('YX')
#                     change_basis.append('N') #included in LEFT optimization
                else:
                    change_basis.append(s)
                break
                
    change_basis.append(sigma_list[-1] if sigma_list[-1]!='I' else 'N') # cannot simplify last term
    new_sigma_list.append('Z' if sigma_list[-1]!='I' else 'I')
    return change_basis, new_sigma_list

T = ['Y','Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
# T = ['Y', 'I','Y', 'X', 'Y']
# T = ['Z', 'Z', 'Z', 'X', 'I']
# T=['Y', 'Z', 'Z', 'Y', 'X', 'X']
change, newT = change_basis_between_terms_RIGHT(T)
# print(T)
print(change)
# print(newT)
newT
T =       ['Y', 'Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
LEFT = ['Y', 'Z', 'Y', 'I', 'YX', 'I', 'N', 'N', 'I']


T =      ['Y', 'Z', 'Y', 'Y', 'X', 'X', 'I', 'I', 'X']
RIGHT =     ['Y', 'Z', 'I', 'YX', 'I', 'I', 'N', 'N', 'X']

print(largest_AC_set)
Change_Basis_op_array(largest_AC_set)
import cirq
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import change_basis_initial, change_basis_final
change_basis_initial('YX', cirq.LineQubit(1))
L = np.array([['A', 'B'], ['C', 'D']], dtype=object)

L[0,:] = ['G', 'FX']
L
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import change_basis_initial, change_basis_final

def reduced_circuit(P_Word_list, angle_list):
    
    Pwords_array= Change_Basis_op_array(P_Word_list)

    Pwords_updated_array = Pwords_array.copy()
    LEFT_change_basis = np.zeros_like(Pwords_updated_array, dtype=object)
    RIGHT_change_basis=np.zeros_like(Pwords_updated_array, dtype=object)
    for qubitNo in range(Pwords_array.shape[0]):

        change_basis_LEFT, P_list_LEFT = change_basis_between_terms_LEFT(Pwords_array[qubitNo,:])
        LEFT_change_basis[qubitNo,:]= change_basis_LEFT
        
        change_basis_RIGHT, P_list_RIGHT = change_basis_between_terms_RIGHT(Pwords_array[qubitNo,:])
        RIGHT_change_basis[qubitNo,:]= change_basis_RIGHT
        
        if P_list_LEFT != P_list_RIGHT:
            raise ValueError('change of basis error')
        Pwords_updated_array[qubitNo,:] = np.asarray(P_list_LEFT)
    
    
    ## Q circuit
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, P_Word_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = list(cirq.LineQubit.range(max_qubits))
    circuit = cirq.Circuit()
#     print(Pwords_array)
#     print()
#     print(Pwords_updated_array)
#     print()
#     print(LEFT_change_basis)
#     print('*********')
    for col_ind in range(Pwords_array.shape[1]):
        change_B_left = LEFT_change_basis[:, col_ind]
        updated_P_word = Pwords_updated_array[:, col_ind]
        change_B_right = RIGHT_change_basis[:, col_ind]
#         print('###')
#         print(change_B_left)
#         print(updated_P_word)
#         print(change_B_right)
#         print(LEFT_change_basis[:, col_ind+1])
#         print()
#         print('old:', Pwords_array[:, col_ind])
        
        # change of basis initial!
        for qNo, P_str in enumerate(change_B_left):
            if P_str != 'N':
                gate = change_basis_initial(P_str, qubits[qNo])
                if gate is not None:circuit.append(gate)
        
        ## initial entangle
#         print(updated_P_word)
        active_locations = np.where(updated_P_word=='Z')[0]
        for i, qNo in enumerate(active_locations[:-1]):
            circuit.append(cirq.CNOT.on(qubits[qNo], qubits[active_locations[i+1]]))
        
        ### PERFORM ROTATION
        coeff_A = list(P_Word_list[col_ind].terms.values())[0]
        theta_A = angle_list[col_ind]
        if coeff_A.imag==0:
            raise ValueError('not valid qubit rotation')
        elif coeff_A.imag < 0:
            circuit.append(cirq.rz(2 * theta_A * np.abs(coeff_A.imag)).on(qubits[active_locations[-1]]))
        else:
            circuit.append(cirq.rz(-2 * theta_A * np.abs(coeff_A.imag)).on(qubits[active_locations[-1]]))
        
        ## FINAL entangler
        rev_active_locations = np.flip(active_locations)
        for i, qNo in enumerate(rev_active_locations[:-1]):
            circuit.append(cirq.CNOT.on(qubits[rev_active_locations[i+1]], qubits[qNo]))
            
        # change of basis final!
        for qNo, P_str in enumerate(change_B_right):
            if (P_str != 'N' and P_str != 'YX' and P_str != 'XY'):
                gate = change_basis_final(P_str, qubits[qNo])
                if gate is not None:circuit.append(gate)
    
    
#     N_qubits = max([lineQ.x for lineQ in list(circuit.all_qubits())]) + 1
#     mat_list = [expm(qubit_operator_sparse(op, n_qubits = N_qubits)*angle_list[ind]).todense() for ind, op in enumerate(P_Word_list)]
#     lin_alg_mat = reduce(np.dot, mat_list[::-1])
#     if not np.allclose(lin_alg_mat, circuit.unitary()):
#         raise ValueError('circuit reduction incorrect!') 

    
    return circuit
change = reduced_circuit([1j*t for t in op_list], [i for i in range(len(op_list))])
change
op_list
Change_Basis_op_array(op_list)

change_basis_between_terms_RIGHT(['X', 'Z', 'X', 'X', 'X'])
from quchem.Qcircuit.Ansatz_quantum_circuit_functions import change_basis_initial, change_basis_final

def reduced_circuit_CNOT_RED(P_Word_list, angle_list):
    
    Pwords_array= Change_Basis_op_array(P_Word_list)

    Pwords_updated_array = Pwords_array.copy()
    LEFT_change_basis = np.zeros_like(Pwords_updated_array, dtype=object)
    RIGHT_change_basis=np.zeros_like(Pwords_updated_array, dtype=object)
    for qubitNo in range(Pwords_array.shape[0]):

        change_basis_LEFT, P_list_LEFT = change_basis_between_terms_LEFT(Pwords_array[qubitNo,:])
        LEFT_change_basis[qubitNo,:]= change_basis_LEFT
        
        change_basis_RIGHT, P_list_RIGHT = change_basis_between_terms_RIGHT(Pwords_array[qubitNo,:])
        RIGHT_change_basis[qubitNo,:]= change_basis_RIGHT
        
        if P_list_LEFT != P_list_RIGHT:
            raise ValueError('change of basis error')
        Pwords_updated_array[qubitNo,:] = np.asarray(P_list_LEFT)
    
    
    ## Q circuit
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, P_Word_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = list(cirq.LineQubit.range(max_qubits))
    circuit = cirq.Circuit()
#     print(Pwords_array)
#     print()
#     print(Pwords_updated_array)
#     print()
#     print(LEFT_change_basis)
#     print('*********')
    for col_ind in range(Pwords_array.shape[1]):
        change_B_left = LEFT_change_basis[:, col_ind]
        updated_P_word = Pwords_updated_array[:, col_ind]
        change_B_right = RIGHT_change_basis[:, col_ind]
#         print('###')
        print(change_B_left)
#         print(updated_P_word)
        print(change_B_right)
#         print(LEFT_change_basis[:, col_ind+1])
#         print()
#         print('old:', Pwords_array[:, col_ind])
        
        # change of basis initial!
        for qNo, P_str in enumerate(change_B_left):
            if P_str != 'N':
                gate = change_basis_initial(P_str, qubits[qNo])
                if gate is not None:circuit.append(gate)
        
        ## initial entangle
        change_B_right_previous = RIGHT_change_basis[:, col_ind-1]
        active_locations = np.where(updated_P_word=='Z')[0]
        left_CNOT_cancel_possible=True
        for i, qNo in enumerate(active_locations[:-1]):
            if (col_ind>0 and left_CNOT_cancel_possible):
                q_next = active_locations[i+1]
                if ((change_B_left[qNo] == change_B_right_previous[qNo]) and (change_B_left[q_next] == change_B_right_previous[q_next])):
                    continue
                else:
                    circuit.append(cirq.CNOT.on(qubits[qNo], qubits[active_locations[i+1]]))
                    left_CNOT_cancel_possible=False
            else:
                circuit.append(cirq.CNOT.on(qubits[qNo], qubits[active_locations[i+1]]))
                
                
        ### PERFORM ROTATION
        coeff_A = list(P_Word_list[col_ind].terms.values())[0]
        theta_A = angle_list[col_ind]
        if coeff_A.imag==0:
            raise ValueError('not valid qubit rotation')
        elif coeff_A.imag < 0:
            circuit.append(cirq.rz(2 * theta_A * np.abs(coeff_A.imag)).on(qubits[active_locations[-1]]))
        else:
            circuit.append(cirq.rz(-2 * theta_A * np.abs(coeff_A.imag)).on(qubits[active_locations[-1]]))
        
        ## FINAL entangler
        right_CNOT_cancel_possible=True
        rev_active_locations = np.flip(active_locations)
        for i, qNo in enumerate(rev_active_locations[:-1]):
            q_next = rev_active_locations[i+1]
            if (col_ind<(Pwords_array.shape[1]-1) and right_CNOT_cancel_possible):
                change_B_left_next = LEFT_change_basis[:, col_ind+1]
                if ((change_B_right[qNo] == change_B_left_next[q_next]) and (change_B_right[q_next] == change_B_left_next[q_next])):
                        continue
                else:
                    circuit.append(cirq.CNOT.on(qubits[rev_active_locations[i+1]], qubits[qNo]))
                    right_CNOT_cancel_possible=False
            else:
                circuit.append(cirq.CNOT.on(qubits[rev_active_locations[i+1]], qubits[qNo]))
            
        # change of basis final!
        for qNo, P_str in enumerate(change_B_right):
            if (P_str != 'N' and P_str != 'YX' and P_str != 'XY'):
                gate = change_basis_final(P_str, qubits[qNo])
                if gate is not None:circuit.append(gate)
        
    return circuit
change = reduced_circuit_CNOT_RED([1j*t for t in op_list], [i for i in range(len(op_list))])
change
op_list
print([1j*t for t in largest_AC_set])
change = reduced_circuit([1j*t for t in largest_AC_set], [i for i in range(len(largest_AC_set))])
change
largest_AC_set
from quchem.Misc_functions.Misc_functions import choose_Pn_index, lexicographical_sort_BASIS_MATCH
key_larg, largest_AC_set = max(anti_commuting_sets.items(), key=lambda x:len(x[1])) # largest nonCon part found by dfs alg
Ps_index = choose_Pn_index(largest_AC_set)
B_Ps = largest_AC_set.pop(Ps_index) # REMOVE BETA_S_P_S from largest_AC_set

BASIS_MATCH_AC_set = lexicographical_sort_BASIS_MATCH(largest_AC_set)
# BASIS_MATCH_AC_set.append(B_Ps) # add Ps back in last index!

c = reduced_circuit([1j*t for t in BASIS_MATCH_AC_set], [i for i in range(len(BASIS_MATCH_AC_set))])
c
[1j*t for t in BASIS_MATCH_AC_set]
circuit_operators = list(c.all_operations())
circuit_operators
from cirq.circuits import InsertStrategy

closest_together = cirq.Circuit(list(BASIS_MATCH_AC_circuit_LADDER_form.all_operations()), strategy=InsertStrategy.EARLIEST)
closest_together

import networkx as nx
from networkx import grid_2d_graph

G = grid_2d_graph(5,5)
nx.draw(G)
import networkx as nx
import matplotlib.pyplot as plt

N_rows = 10
N_cols = 
G=nx.grid_2d_graph(N,N)
pos = dict( (n, n) for n in G.nodes() )
labels = dict( ((i, j), i + (N-1-j) * 10 ) for i, j in G.nodes() )
nx.draw_networkx(G, pos=pos, labels=labels)

plt.axis('off')
plt.show()
G.nodes()
A
BASIS_MATCH_AC_set
from openfermion.utils import count_qubits
from functools import reduce
def ladder_or_stairs_circuit_array(P_Word_list, ladder=True):
    """

    ## example
    from openfermion import QubitOperator
    import numpy as np

    op1 = QubitOperator('Y0 X1', -1j)
    op2 = QubitOperator('Y0 Z1', -1j)
    op3 = QubitOperator('Z0 X1', -1j)
    OP_list = [op1, op2, op3]


    circuit_arr =ladder_circuit_array(OP_list, ladder=True)
    print(circuit_arr)
    
     >> array([['H', 'S', 'C', 'I', 'C', 'Sdag', 'H', 'H', 'S', 'C', 'I', 'C','Sdag', 'H', 'I', 'C', 'I', 'C', 'I'],
               ['I', 'H', 'X', 'R', 'X', 'H',    'I', 'I', 'I', 'X', 'R', 'X', 'I',   'I', 'H', 'X', 'R', 'X', 'H']],
                dtype=object)
    """
    
    change_basis_initial_dict = {
        'X': 'H',
        'Y': 'Sdag',# S_dag then H (HSdag) 
        'Z': 'I'
    }
    
    change_basis_final_dict = {
        'X': 'H',
        'Y': 'S',#' H then S (HS)
        'Z': 'I'
    }
    
    
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, P_Word_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = cirq.LineQubit.range(max_qubits)
    
    empty_circuit_slice = np.array(['I' for _ in range(max_qubits)], dtype=object)
    circuit_array = np.empty((len(qubits),0), dtype=object)
    for ind, P_op in enumerate(P_Word_list):
        PauliStrs, coeff = tuple(*P_op.terms.items())
        qNos, sigmas = zip(*list(PauliStrs))
        
        # initial change basis
        circuit_slice = empty_circuit_slice.copy()
        Y_slice= empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_initial_dict[s]
            if s == 'Y':
                Y_slice[qNo]='H'
                
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        if 'H' in Y_slice:
            circuit_array = np.hstack((circuit_array, Y_slice.reshape(((max_qubits),1))))
            
            
        ### entangle_initial
        entangle_arr = np.empty((len(qubits),0), dtype=object)
        if ladder is True:
            # LADDER approach
            for ind, qNo in enumerate(qNos[:-1]):
                circuit_slice = empty_circuit_slice.copy()
                circuit_slice[qNo] = 'C'
                circuit_slice[qNos[ind+1]] = 'X'
                entangle_arr = np.hstack((entangle_arr, circuit_slice.reshape(((max_qubits),1))))
        else:
            # STAIRS approach
            target_qubit = qNos[-1]
            for ind, qNo in enumerate(qNos[:-1]):
                circuit_slice = empty_circuit_slice.copy()
                circuit_slice[qNo] = 'C'
                circuit_slice[target_qubit] = 'X'
                entangle_arr = np.hstack((entangle_arr, circuit_slice.reshape(((max_qubits),1))))
        circuit_array = np.hstack((circuit_array, entangle_arr))
        ####
        

        ## rotation
        circuit_slice = empty_circuit_slice.copy()
        circuit_slice[qNos[-1]]='R'
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        
        ## final entangle
        entangle_final = entangle_arr[:,::-1]
        circuit_array = np.hstack((circuit_array, entangle_final))
        
        
        # final change basis
        circuit_slice = empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_final_dict[s]
        
        if 'H' in Y_slice:
            circuit_array = np.hstack((circuit_array, Y_slice.reshape(((max_qubits),1))))
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))              
    return circuit_array
op1 = QubitOperator('X0 X1 ', -1j)
op2 = QubitOperator('X0 X1', -1j)
op3 = QubitOperator('Y0 Y1', -1j)

# op1 = QubitOperator('Y0', -1j)
# op2 = QubitOperator('X0', -1j)
# op3 = QubitOperator('X0', -1j)

# op1 = QubitOperator('Y0 X1', -1j)
# op2 = QubitOperator('Y0 Z1', -1j)
# op3 = QubitOperator('X1', -1j)

op_list = [op1, op2, op3]

out = ladder_or_stairs_circuit_array(op_list)
print(out)
out
def cancel_common_basis(circuit_array):
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig in enumerate(qubit_operations):
            if (sig == 'I') or (sig == 'R'):
                continue
            for new_ind in range(ind+1, circuit_array.shape[1]):
                sig_new = qubit_operations[new_ind]
                if sig_new == 'C':
                    break
                elif sig == sig_new:
                    qubit_operations[ind]='I'
                    qubit_operations[new_ind]='I'
                    break
                elif sig_new == 'I':
                    continue
                elif sig!=sig_new:
                    break
    return circuit_array
print(out)
print()
out2= cancel_common_basis(out)
print(out2)
[['I' 'I' 'I' 'I' 'I' 'I' 'I' 'H' 'C' 'I' 'I' 'I' 'C' 'H' 'H' 'S' 'C' 'I''I' 'I' 'C' 'Sdag' 'H']
 ['I' 'H' 'C' 'I' 'C' 'H' 'I' 'H' 'X' 'C' 'I' 'C' 'X' 'H' 'H' 'S' 'X' 'C' 'I' 'C' 'X' 'Sdag' 'H']
 ['H' 'S' 'X' 'R' 'X' 'Sdag' 'H' 'H' 'I' 'X' 'R' 'X' 'I' 'H' 'I' 'H' 'I''X' 'R' 'X' 'I' 'H' 'I']]
def cancel_S_gates_through_controls(circuit_array):
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig in enumerate(qubit_operations):
            if sig not in ['S', 'Sdag']:
                continue
            for new_ind in range(ind+1, circuit_array.shape[1]):
                sig_new = qubit_operations[new_ind]
                if sig_new in ['X', 'H', 'R']:
                    # cannot move through these gates
                    break
                elif sig_new in ['I','C']:
                    # S gates commute with controls
                    continue
                elif (sig == 'S') and (sig_new == 'Sdag'):
                    qubit_operations[ind]='I'
                    qubit_operations[new_ind]='I'
                    break
                elif (sig == 'Sdag') and (sig_new == 'S'):
                    qubit_operations[ind]='I'
                    qubit_operations[new_ind]='I'
                    break
                else:
                    raise ValueError(f'unkown opertion:{sig, sig_new}')
    return circuit_array
print(out2)
print()
out3= cancel_S_gates_through_controls(out2)
print(out3)
out3
def cancel_CNOT_Gates(circuit_array):
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]-1):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig_control in enumerate(qubit_operations):
            if sig_control != 'C':
                continue
            
            next_row_ind = np.where(circuit_array[:,ind]=='X')[0][0]
            qubit_operations_NEXT = circuit_array[next_row_ind, :]
            
            sig_target = qubit_operations_NEXT[ind]
            if sig_target != 'X':
                continue
            
            for new_ind in range(ind+1, circuit_array.shape[1]):
                sig_new_CONTROL = qubit_operations[new_ind]
                sig_new_Target = qubit_operations_NEXT[new_ind]
                
                if (sig_new_CONTROL in ['I', 'S', 'Sdag'] and sig_new_Target in ['I']):
                    continue
                elif sig_new_CONTROL !='C':
                    break
                elif sig_new_Target !='X':
                    break
                else:
                    qubit_operations[ind]='I'
                    qubit_operations_NEXT[ind]='I'
                    qubit_operations[new_ind]='I'
                    qubit_operations_NEXT[new_ind]='I'
                    break
    return circuit_array
out4 = cancel_CNOT_Gates(out3)
print(out4)
def circuit_array_to_circuit(circuit_array, angle_list):
    
    qubits = cirq.LineQubit.range(circuit_array.shape[1])
    circuit = cirq.Circuit()
    
    angle_ind=0
    for circuit_slice in circuit_array.T:
        non_I_indices = np.where(circuit_slice!='I')[0]

        if 'C' in circuit_slice:
            circuit.append(cirq.CNOT(qubits[non_I_indices[0]], qubits[non_I_indices[1]]))
        else:
            for qNo in non_I_indices:
                op = circuit_slice[qNo]
                if op == 'H':
                    circuit.append(cirq.H(qubits[qNo]))
                elif op == 'R':
                    circuit.append(cirq.rz(angle_list[angle_ind]).on(qubits[qNo]))
                    angle_ind+=1
                elif op == 'S':
                    circuit.append(cirq.S(qubits[qNo]))
                elif op == 'Sdag':
                    circuit.append((cirq.S**-1).on(qubits[qNo]))
                else:
                    print('unknown operation!')
    return circuit
    
    
circuit_array_to_circuit(out4, [i for i in range(len(op_list))])
op_list
test = cirq.Circuit()
def stairs_circuit_array(P_Word_list):
    """

    ## example
    from openfermion import QubitOperator
    import numpy as np

    op1 = QubitOperator('Y0 X1', -1j)
    op2 = QubitOperator('Y0 Z1', -1j)
    op3 = QubitOperator('Z0 X1', -1j)
    OP_list = [op1, op2, op3]


    circuit_arr =ladder_circuit_array(OP_list)
    print(circuit_arr)
    
     >> array([['H', 'S', 'C', 'I', 'C', 'Sdag', 'H', 'H', 'S', 'C', 'I', 'C','Sdag', 'H', 'I', 'C', 'I', 'C', 'I'],
               ['I', 'H', 'X', 'R', 'X', 'H',    'I', 'I', 'I', 'X', 'R', 'X', 'I',   'I', 'H', 'X', 'R', 'X', 'H']], dtype=object)
    """
    
    change_basis_initial_dict = {
        'X': 'H',
        'Y': 'S',#'H.S',
        'Z': 'I'
    }
    
    change_basis_final_dict = {
        'X': 'H',
        'Y': 'Sdag',#'Sdag.H',
        'Z': 'I'
    }
    
    
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, P_Word_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = cirq.LineQubit.range(max_qubits)
    
    empty_circuit_slice = np.array(['I' for _ in range(max_qubits)], dtype=object)
    circuit_array = np.empty((len(qubits),0), dtype=object)
    for ind, P_op in enumerate(P_Word_list):
        PauliStrs, coeff = tuple(*P_op.terms.items())
        qNos, sigmas = zip(*list(PauliStrs))
        
        # initial change basis
        circuit_slice = empty_circuit_slice.copy()
        Y_slice= empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_initial_dict[s]
            if s == 'Y':
                Y_slice[qNo]='H'
        if 'H' in Y_slice:
            circuit_array = np.hstack((circuit_array, Y_slice.reshape(((max_qubits),1))))
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))

        # entangle_initial
        target_qubit = qNos[-1]
        entangle_arr = np.empty((len(qubits),0), dtype=object)
        for ind, qNo in enumerate(qNos[:-1]):
            circuit_slice = empty_circuit_slice.copy()
            circuit_slice[qNo] = 'C'
            circuit_slice[target_qubit] = 'X'
            entangle_arr = np.hstack((entangle_arr, circuit_slice.reshape(((max_qubits),1))))
        circuit_array = np.hstack((circuit_array, entangle_arr))
        
        
        ## rotation
        circuit_slice = empty_circuit_slice.copy()
        circuit_slice[qNos[-1]]='R'
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        
        ## final entangle
        entangle_final = entangle_arr[:,::-1]
        circuit_array = np.hstack((circuit_array, entangle_final))
        
        
        # final change basis
        circuit_slice = empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_final_dict[s]
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        
        if 'H' in Y_slice:
            circuit_array = np.hstack((circuit_array, Y_slice.reshape(((max_qubits),1))))
            
    return circuit_array
op1 = QubitOperator('Y0 X1 X2', -1j)
op2 = QubitOperator('Y0 X1 Y2', -1j)
op3 = QubitOperator('Y0 Y1 Z2', -1j)

# op1 = QubitOperator('Y0', -1j)
# op2 = QubitOperator('X0', -1j)
# op3 = QubitOperator('X0', -1j)


op_list = [op1, op2, op3]

out = stairs_circuit_array(op_list)
print(out)
out
import cirq
import numpy as np

from openfermion.utils import count_qubits
from functools import reduce

def ladder_or_stairs_circuit_array(P_Word_list, ladder=True):
    """

    ## example
    from openfermion import QubitOperator
    import numpy as np

    op1 = QubitOperator('Y0 X1', -1j)
    op2 = QubitOperator('Y0 Z1', -1j)
    op3 = QubitOperator('Z0 X1', -1j)
    OP_list = [op1, op2, op3]


    circuit_arr =ladder_circuit_array(OP_list, ladder=True)
    print(circuit_arr)
    
     >> array([['H', 'S', 'C', 'I', 'C', 'Sdag', 'H', 'H', 'S', 'C', 'I', 'C','Sdag', 'H', 'I', 'C', 'I', 'C', 'I'],
               ['I', 'H', 'X', 'R', 'X', 'H',    'I', 'I', 'I', 'X', 'R', 'X', 'I',   'I', 'H', 'X', 'R', 'X', 'H']],
                dtype=object)
    """
    
    change_basis_initial_dict = {
        'X': 'H',
        'Y': 'SH',# S_dag then H (HSdag) 
        'Z': 'I'
    }
    
    change_basis_final_dict = {
        'X': 'H',
        'Y': 'HS',#' H then S (HS)
        'Z': 'I'
    }
    
    
    full_QubitOp = reduce(lambda Op1, Op2: Op1+Op2, P_Word_list)
    max_qubits = count_qubits(full_QubitOp)
    qubits = cirq.LineQubit.range(max_qubits)
    
    empty_circuit_slice = np.array(['I' for _ in range(max_qubits)], dtype=object)
    circuit_array = np.empty((len(qubits),0), dtype=object)
    for ind, P_op in enumerate(P_Word_list):
        PauliStrs, coeff = tuple(*P_op.terms.items())
        qNos, sigmas = zip(*list(PauliStrs))
        
        # initial change basis
        circuit_slice = empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_initial_dict[s]
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
            
            
        ### entangle_initial
        entangle_arr = np.empty((len(qubits),0), dtype=object)
        if ladder is True:
            # LADDER approach
            for ind, qNo in enumerate(qNos[:-1]):
                circuit_slice = empty_circuit_slice.copy()
                circuit_slice[qNo] = 'C'
                circuit_slice[qNos[ind+1]] = 'X'
                entangle_arr = np.hstack((entangle_arr, circuit_slice.reshape(((max_qubits),1))))
        else:
            # STAIRS approach
            target_qubit = qNos[-1]
            for ind, qNo in enumerate(qNos[:-1]):
                circuit_slice = empty_circuit_slice.copy()
                circuit_slice[qNo] = 'C'
                circuit_slice[target_qubit] = 'X'
                entangle_arr = np.hstack((entangle_arr, circuit_slice.reshape(((max_qubits),1))))
        circuit_array = np.hstack((circuit_array, entangle_arr))
        ####
        

        ## rotation
        circuit_slice = empty_circuit_slice.copy()
        circuit_slice[qNos[-1]]='R'
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        
        ## final entangle
        entangle_final = entangle_arr[:,::-1]
        circuit_array = np.hstack((circuit_array, entangle_final))
        
        
        # final change basis
        circuit_slice = empty_circuit_slice.copy()
        for qNo, s in zip(qNos, sigmas):
            circuit_slice[qNo] = change_basis_final_dict[s]
        circuit_array = np.hstack((circuit_array, circuit_slice.reshape(((max_qubits),1))))
        
    return circuit_array
op1 = QubitOperator('Y0 Y1 X2', -1j)
op2 = QubitOperator('Y0 X1 Y2', -1j)
op3 = QubitOperator('Y0 X1 Z2', -1j)

# op1 = QubitOperator('Y0', -1j)
# op2 = QubitOperator('X0', -1j)
# op3 = QubitOperator('X0', -1j)


op_list = [op1, op2, op3]
c_arr = ladder_or_stairs_circuit_array(op_list)
c_arr
def cancel_common_change_basis(circuit_array):
    """
    Given a circuit array function will cancel any common change of basis gates
    """
    circuit_array=circuit_array.copy()
    for qubit_row in range(circuit_array.shape[0]):
        qubit_operations = circuit_array[qubit_row, :]
        
        for ind, sig in enumerate(qubit_operations):
            if sig in ['I', 'R' 'C']:
                continue

            for new_ind in range(ind+1, circuit_array.shape[1]):
                sig_new = qubit_operations[new_ind]
                if sig_new in ['C', 'X']:
                    break
                elif sig == sig_new[::-1]:
                    qubit_operations[ind]='I'
                    qubit_operations[new_ind]='I'
                    break
                elif sig_new == 'I':
                    continue
                elif sig!=sig_new:
                    break
    return circuit_array
basis_arr = cancel_common_change_basis(c_arr)
basis_arr
def rewrite_Y_basis_change(circuit_array):
    """
    Re-write anywhere with HS or SH as seperate operations
    """
    
    print(circuit_array)
    print()
    
    circuit_array=circuit_array.copy()
    new_circuit_array = np.empty((circuit_array.shape[0],0), dtype=object)
    empty_circuit_slice = np.array(['I' for _ in range(circuit_array.shape[0])], dtype=object)
    
    for col_ind in range(circuit_array.shape[1]):
        term_arr = circuit_array[:, col_ind].copy()
        emplty_slice = empty_circuit_slice.copy()
        
        if 'S' in ''.join(term_arr):
            SH = np.where(term_arr=='SH')[0]
            HS = np.where(term_arr=='HS')[0]
            
            term_arr[SH]= 'Sdag'
            emplty_slice[SH]='H'
            
            term_arr[HS]= 'H'
            emplty_slice[HS]='S'
            
            
            new_slice = np.vstack((term_arr, empty_circuit_slice)).T
        else:
            new_slice = term_arr.reshape((circuit_array.shape[0],1))

        new_circuit_array = np.hstack((new_circuit_array, new_slice))
        
    return new_circuit_array
expanded_arr = rewrite_Y_basis_change(basis_arr)
expanded_arr
[['S', 'H', 'C', 'I', 'I', 'I', 'C', 'I', 'H', 'I', 'H', 'C', 'I','I', 'I', 'C', 'I', 'H', 'I', 'C', 'I', 'I', 'I', 'C', 'H', 'S'],
 ['S', 'H', 'X', 'C', 'I', 'C', 'X', 'H', 'S', 'H', 'S', 'X', 'C','I', 'C', 'X', 'I', 'S', 'I', 'X', 'C', 'I', 'C', 'X', 'H', 'S'],
 ['H', 'I', 'I', 'X', 'R', 'X', 'I', 'H', 'I', 'S', 'H', 'I', 'X','R', 'X', 'I', 'H', 'S', 'I', 'I', 'X', 'R', 'X', 'I', 'I', 'S']],
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import cancel_S_gates_through_controls
S_cancel = cancel_S_gates_through_controls(expanded_arr)
S_cancel
basis_arr2 = cancel_common_change_basis(S_cancel)
basis_arr2
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import cancel_CNOT_Gates
CNOT_cancel = cancel_CNOT_Gates(basis_arr2)
CNOT_cancel
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import circuit_array_to_circuit
circuit_array_to_circuit(CNOT_cancel, [i for i in range(len(op_list))], op_list)
from scipy.sparse.linalg import expm
from openfermion import qubit_operator_sparse
from functools import reduce
def Optimized_LADDER_circuit(P_Word_list, angle_list, check_reduction=False):

    circuit_array=ladder_or_stairs_circuit_array(P_Word_list, ladder=True)
    circuit_array_basis_cancelled=cancel_common_change_basis(circuit_array)
    circuit_array_S_cancelled=cancel_S_gates_through_controls(circuit_array_basis_cancelled)
    circuit_array_basis_cancelled=cancel_common_change_basis(circuit_array_S_cancelled)
    circuit_array_CNOT_cancel=cancel_CNOT_Gates(circuit_array_basis_cancelled)
    opt_circuit =circuit_array_to_circuit(circuit_array_CNOT_cancel, angle_list, P_Word_list)

    if check_reduction:
        N_qubits = max([lineQ.x for lineQ in list(opt_circuit.all_qubits())]) + 1
        mat_list = [expm(qubit_operator_sparse(op, n_qubits = N_qubits)*angle_list[ind]).todense() for ind, op in enumerate(P_Word_list)]
        lin_alg_mat = reduce(np.dot, mat_list[::-1])

        if not np.allclose(lin_alg_mat, opt_circuit.unitary()):
            raise ValueError('circuit reduction incorrect!') 

    return opt_circuit
CNOT_cancel
op_list
Optimized_LADDER_circuit(op_list, [i for i in range(len(op_list))], check_reduction=True)
x_pow_arr = YX_simplification(CNOT_cancel)
circuit_array_to_circuit(x_pow_arr, [i for i in range(len(op_list))], op_list)
from quchem.Qcircuit.Optimized_exp_pauliword_circuit_functions import *
c_arr = ladder_or_stairs_circuit_array(op_list)
circuit_array_basis_cancelled=cancel_common_change_basis(c_arr)
expanded_circuit_array = rewrite_Y_basis_change(circuit_array_basis_cancelled)
expanded_circuit_array
circuit_array_basis_cancelled
ff = np.insert(circuit_array_basis_cancelled,
          col_HS+1,
          I_arr,
          axis=1)
ff<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12381561190425457, 'fraction_numerical': 0.01105192130122878, 'mean_word_length': 2.7494540612516642, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 107, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '479686', 'license': 'MIT', 'n_tokens_mistral': 26562, 'n_tokens_neox': 25050, 'n_words': 4470, 'path': 'quchem_examples/3.2.3 notes - to delete.ipynb', 'repo_name': 'AlexisRalli/VQE-code', 'size': '70395'}