QCut to Three Parts

import QCut as ck
from QCut import cut
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.primitives import Estimator

#define initial circuit

circuit = QuantumCircuit(4)
circuit.h(0)
circuit.cx(0,1)
circuit.cx(1,2)
circuit.cx(2,3)
circuit.measure_all()

circuit.draw("mpl")

#insert cuts

cut_circuit = QuantumCircuit(4)
cut_circuit.h(0)
cut_circuit.cx(0,1)
cut_circuit.append(cut, [1])
cut_circuit.cx(1,2)
cut_circuit.append(cut, [2])
cut_circuit.cx(2,3)

cut_circuit.draw("mpl")

#find cut locations and separate into subcircuits

cut_locations, subcircuits, map_qubit = ck.get_locations_and_subcircuits(cut_circuit)

subcircuits[0].draw("mpl")

subcircuits[1].draw("mpl")

#insert qpd operations to get the 8^n experiment circuits, where n is the number of cuts made
#also returns coefficients used for reconstructing the expectation values
#and locations of identity basis measurements that get added with post processing

experiment_circuits, coefficients, id_meas = ck.get_experiment_circuits(subcircuits, cut_locations)

#define backend
backend = AerSimulator()

#run experiment circuits
#run_experiments() also post processes the results

results = ck.run_experiments(experiment_circuits, cut_locations, id_meas, backend=backend, mitigate=False)

#define observables to calculate expectation values for
observables = [0,1,2,3, [0,2], [0,1,3]]

#get the exact expectation values with ideal simulator
expectation_values = ck.estimate_expectation_values(results, coefficients, cut_locations, observables, map_qubit)

#transform observables into a PauliList for the Qiskit Estimator
paulilist_observables = ck.get_pauli_list(observables, circuit.num_qubits)

#get the exact expectation values with ideal simulator
estimator = Estimator(run_options={"shots": None}, approximation=True)
exact_expvals = (
    estimator.run([circuit] * len(paulilist_observables),  # noqa: PD011
                  list(paulilist_observables)).result().values
)

import numpy as np

#set numpy print options
np.set_printoptions(formatter={"float": lambda x: f"{x:0.6f}"})
expectation_values = np.array(expectation_values).astype(float)

#compare QCut to Estimator
print(f"QCut circuit knitting expectation values:{np.array(expectation_values)}")
print(f"Exact expectation values with ideal simulator :{np.array(exact_expvals)}")

QCut circuit knitting expectation values:[-0.011825 -0.011825 0.029135 0.006827 1.006896 0.006827]
Exact expectation values with ideal simulator :[0.000000 0.000000 0.000000 0.000000 1.000000 0.000000]