router mishandling mid-circuit measurements

[richrines1]

Description of the issue

cirq.RouteCQC will sometimes commute gates through mid-circuit measurements in a way that changes the expected outcome

How to reproduce the issue

q0, q1, q2 = cirq.LineQubit.range(3)

circuit = cirq.Circuit(
    cirq.CZ(q0, q1),
    cirq.CZ(q1, q2),
    cirq.CZ(q0, q2),
    cirq.measure(q0, q1, q2, key="key"),
    cirq.H.on_each(q0, q1, q2),
    cirq.measure_each(q0, q21, q2),
)

graph = nx.Graph([[q0, q1], [q1, q2]])
router = cirq.RouteCQC(graph)

print(circuit)
print(router.route_circuit(circuit)[0])

prints

0: ───@───────@───M('key')───H───M───
      │       │   │
1: ───@───@───┼───M──────────H───M───
          │   │   │
2: ───────@───@───M──────────H───M───
0: ───×───────@───────────────@───M──────────H───M('q(2)')───
      │       │               │   │
1: ───×───@───@───H───M───×───@───M('key')───H───M('q(0)')───
          │               │       │
2: ───────@───────────────×───────M──────────────────────────

where one of the hadamards and one of the terminal measurement have incorrectly been commuted through the mid-circuit measurement

Cirq version

1.3.0.dev20230830191034

[shef4]

I can look into this issue and try to reproduce it.

[shef4]

was able to reproduce the issue, gonna start looking through RouteCQC class and test.

[shef4]

It seems like the routing algorithm is sometimes reordering gates in a quantum circuit in a way that they are incorrectly placed before or after measurements, leading to discrepancies between the expected and actual outcomes of the circuit. New to opensource contribution so please let me know if I am misunderstanding the issue.

I am going to start looking at which stage of gates we start seeing a mismatch between the actual and expected circuit.

[shef4]

After some investigation seems like the issue is using measurement gates in the circuit. Cirq has a test case to check behavior with measurements.

def test_circuit_with_measurement_gates():
    device = cirq.testing.construct_ring_device(3)
    device_graph = device.metadata.nx_graph
    q = cirq.LineQubit.range(3)
    circuit = cirq.Circuit(cirq.MeasurementGate(2).on(q[0], q[2]), cirq.MeasurementGate(3).on(*q))
    hard_coded_mapper = cirq.HardCodedInitialMapper({q[i]: q[i] for i in range(3)})
    router = cirq.RouteCQC(device_graph)
    routed_circuit = router(circuit, initial_mapper=hard_coded_mapper)
    cirq.testing.assert_same_circuits(routed_circuit, circuit)

Going to compare this example with the test case and see if there needs to be a new test case for this situation.

I tried out multiple versions of the code based on the test case implementation. Outside of modifying the cirq.RouteCQC function, changing the nx network graph device and adding initial_mapper=hard_coded_mapper seem to affect the result.

Debug Output
Constructed graph from cirq.testing.construct_ring_device(3).metadata.nx_graph with initial_mapperr=hard_coded_mapper parameter:

import cirq 
import networkx as nx

device = cirq.testing.construct_ring_device(3)
device_graph = device.metadata.nx_graph

q = cirq.LineQubit.range(3)

circuit = cirq.Circuit(cirq.CZ(q[0], q[1]), cirq.CZ(q[1], q[2]), cirq.CZ(q[0], q[2]),cirq.MeasurementGate(3).on(*q),cirq.H.on_each(*q), cirq.measure_each(*q))

hard_coded_mapper = cirq.HardCodedInitialMapper({q[i]: q[i] for i in range(3)})

router = cirq.RouteCQC(device_graph)
routed_circuit = router(circuit, initial_mapper=hard_coded_mapper)

print(circuit)
print(routed_circuit)

result:

0: ───@───────@───M('')───H───M───
      │       │   │
1: ───@───@───┼───M───────H───M───
          │   │   │
2: ───────@───@───M───────H───M───
              ┌──┐
0: ───@─────────@────────M('')───H───M───
      │         │        │
1: ───@───@────H┼────M───M───────────────
          │     │        │
2: ───────@─────@────────M───────H───M───

[shef4]

Hey @richrines1 was able to get a 2nd opinion, the issue seems to be that the algorithm was originally designed with 2 qubit circuits in mind. So there is a small bug in one of the helper functions which is used to decomposes the single and two quit circuit operations into timesteps.

[shef4]

Hey @tanujkhatta, step 3 in the route_circuit function for RouteCQC class was causing the issue. Thanks for the hint!

# 3. Get two_qubit_ops and single-qubit operations.
two_qubit_ops, single_qubit_ops=self._get_one_and_two_qubit_ops_as_timesteps(circuit)

by changing this if statement if protocols.num_qubits(op) == 2 and not protocols.is_measurement(op): to if protocols.num_qubits(op) >= 2: the output looks fixed.

0: ───@───────@───M('')───H───M───
      │       │   │
1: ───@───@───┼───M───────H───M───
          │   │   │
2: ───────@───@───M───────H───M───

Would this fix work or is it necessary that measurements are excluded from the two_qubit_ops?

I can also update the variable names from two_qubit_... to multi_qubit_... for _get_one_and_two_qubit_ops_as_timesteps function?

 @classmethod
    def _get_one_and_two_qubit_ops_as_timesteps(
        cls, circuit: 'cirq.AbstractCircuit'
    ) -> Tuple[List[List['cirq.Operation']], List[List['cirq.Operation']]]:
        """Gets the single and two qubit operations of the circuit factored into timesteps.

        The i'th entry in the nested two-qubit and single-qubit ops correspond to the two-qubit
        gates and single-qubit gates of the i'th timesteps respectively. When constructing the
        output routed circuit, single-qubit operations are inserted before two-qubit operations.
        """
        two_qubit_circuit = circuits.Circuit()
        single_qubit_ops: List[List[cirq.Operation]] = []
        for moment in circuit:
            for op in moment:
                timestep = two_qubit_circuit.earliest_available_moment(op)
                single_qubit_ops.extend([] for _ in range(timestep + 1 - len(single_qubit_ops)))
                two_qubit_circuit.append(
                    circuits.Moment() for _ in range(timestep + 1 - len(two_qubit_circuit))
                )
                if protocols.num_qubits(op) >= 2:# and not protocols.is_measurement(op):
                    two_qubit_circuit[timestep] = two_qubit_circuit[timestep].with_operation(op)
                else:
                    single_qubit_ops[timestep].append(op)
        two_qubit_ops = [list(m) for m in two_qubit_circuit]
        return two_qubit_ops, single_qubit_ops

[shef4]

@richrines1 tried it out on the original example and looks like the expected result, given the network graph = nx.Graph([[q0, q1], [q1, q2]]).

import cirq
import networkx as nx
q0, q1, q2 = cirq.LineQubit.range(3)

circuit = cirq.Circuit(
    cirq.CZ(q0, q1),
    cirq.CZ(q1, q2),
    cirq.CZ(q0, q2),
    cirq.measure(q0, q1, q2, key="key"),
    cirq.H.on_each(q0, q1, q2),
    cirq.measure_each(q0, q1, q2),
)

graph = nx.Graph([[q0, q1], [q1, q2]])

hard_coded_mapper = cirq.HardCodedInitialMapper({q[i]: q[i] for i in range(3)})
router = cirq.RouteCQC(graph)
rcirc, initial_map, swap_map = router.route_circuit(circuit)
print(circuit)
print(router.route_circuit(circuit, initial_mapper=hard_coded_mapper)[0])

result:

0: ───@───────@───M('key')───H───M───
      │       │   │
1: ───@───@───┼───M──────────H───M───
          │   │   │
2: ───────@───@───M──────────H───M───
0: ───×───────@───────@───M──────────H───M('q(2)')───
      │       │       │   │
1: ───×───@───@───×───@───M('key')───H───M('q(0)')───
          │       │       │
2: ───────@───────×───────M──────────H───M('q(1)')───

[NoureldinYosri]

@shef4 per cirqcync, we will go about this in stages

1. Errorout if there are intermediate measurements on 3+ qubits.
2. Break intermediate measurements on 3+ qubits into single qubit measurements if -and only if- the result isn't stored (i.e. key is None)
3. Check if it's possible to create a ClassicallyControlledOperation controlled by the result of a multiqubit measurement
   • if it's possible then we stop here until there is a need to support this case because then we will know how this should be handled.
   • if it's not possible then break all intermediate multiqubit measurements into single measurements even those that have a key (because they are not using it anyway) while logging a warning to the user that the key is being dropped.

[shef4]

@NoureldinYosri Thanks, will start work working on 1st stage.

[shef4]

@NoureldinYosri Started working on Stage 2 of implementation and was wondering if there is a preference between methods 1 & 2.

Method 1: Filter measurement gates outside of the current algorithm - possible helper function

        two_qubit_circuit = circuits.Circuit()
        single_qubit_ops: List[List[cirq.Operation]] = []
        filtered_circuit = circuits.Circuit()

        if any(
            protocols.num_qubits(op) > 2 and protocols.is_measurement(op)
            for op in itertools.chain(*circuit.moments[:-1])
        ):
            # There is at least one non-terminal measurement on 3+ qubits
            for moment in circuit:
                for op in moment:
                    timestep = filtered_circuit.earliest_available_moment(op)
                    if protocols.num_qubits(op) < 3 or not protocols.is_measurement(op):
                        filtered_circuit[timestep].append(op)   
                    elif op.gate.key == "":
                        filtered_circuit[timestep] += [cirq.MeasurementGate(1).on(qubit) for qubit in multi_qubit_measurment.qubits]
                    else:
                        raise ValueError('Non-terminal measurements on three or more qubits when result is stored are not supported')      
            circuit = filtered_circuit

        for moment in circuit:
            for op in moment:
                timestep = two_qubit_circuit.earliest_available_moment(op)
                single_qubit_ops.extend([] for _ in range(timestep + 1 - len(single_qubit_ops)))
                two_qubit_circuit.append(
                    circuits.Moment() for _ in range(timestep + 1 - len(two_qubit_circuit))
                )
                if protocols.num_qubits(op) == 2:
                    two_qubit_circuit[timestep] = two_qubit_circuit[timestep].with_operation(op)
                else:
                    single_qubit_ops[timestep].append(op)
        two_qubit_ops = [list(m) for m in two_qubit_circuit]
        return two_qubit_ops, single_qubit_ops

Method 2: Incorporate filtering measurements into algorithm

        two_qubit_circuit = circuits.Circuit()
        single_qubit_ops: List[List[cirq.Operation]] = []

        for moment in circuit:
            for op in moment:
                timestep = two_qubit_circuit.earliest_available_moment(op)
                single_qubit_ops.extend([] for _ in range(timestep + 1 - len(single_qubit_ops)))
                two_qubit_circuit.append(
                    circuits.Moment() for _ in range(timestep + 1 - len(two_qubit_circuit))
                )
                if protocols.num_qubits(op) > 2 and protocols.is_measurement(op):
                    if op.gate.key == "":
                        single_qubit_ops[timestep] += [cirq.MeasurementGate(1).on(qubit) for qubit in op.gate.qubits]
                    else:
                        raise ValueError('Non-terminal measurements on three or more qubits when result is stored are not supported')
                elif protocols.num_qubits(op) == 2:
                    two_qubit_circuit[timestep] = two_qubit_circuit[timestep].with_operation(op)
                else:
                    single_qubit_ops[timestep].append(op)
        two_qubit_ops = [list(m) for m in two_qubit_circuit]
        return two_qubit_ops, single_qubit_ops

[NoureldinYosri]

@shef4 lets go with option 2

[shef4]

@NoureldinYosri Thanks, I was able to get option 2 working and my tests are passing but test_circuit_with_measurement_gates() is now failing.

#Diagram of expected circuit:
0: ───M('')───M('')───
      │       │
1: ───┼───────M───────
      │       │
2: ───M───────M───────

#Diagram of actual circuit:
0: ───M('')───M───
      │
1: ───┼───────M───
      │
2: ───M───────M───

Since it's a terminal measurement on 3+ qubits, I'm planning on adding it to the two_qubit_circuit list.

        ...
        current_moment = 0

        for moment in circuit:
            current_moment += 1
            for op in moment:
                ...
                if protocols.num_qubits(op) > 2 and protocols.is_measurement(op):
                    if len(circuit.moments) == current_moment:
                        two_qubit_circuit[timestep] = two_qubit_circuit[timestep].with_operation(op)
                    elif op.gate.key == "":
                        ...
                    else:
                        ...
                elif protocols.num_qubits(op) == 2:
                    ...
                else:
                    ...
        ...
        return two_qubit_ops, single_qubit_ops

I could also update the test and leave implementation as is.

[NoureldinYosri]

@shef4 multiqubit terminal states are currently handled correctly. Look at how the algorithm handles them currently and preserve that behaviour, don't change how they are handled. your work should only affect intermediate measurements.

[shef4]

@NoureldinYosri I think I see what you're saying. The original algorithm would place measurements in single_qubits_ops. I'll maintain the same input-to-output path for functions that are dependent on it.

if protocols.num_qubits(op) == 2 and not protocols.is_measurement(op):
   two_qubit_circuit[timestep] = two_qubit_circuit[timestep].with_operation(op)
else:
   single_qubit_ops[timestep].append(op)

[shef4]

Part 3: From looking at the documentation it seems to be possible to create a ClassicallyControlledOperation controlled by the result of a multiqubit measurement. Will stop here until there is a need to support this case.

cirq.ClassicallyControlledOperation(
    sub_operation: 'cirq.Operation',
    conditions: Sequence[Union[str, 'cirq.MeasurementKey', 'cirq.Condition', sympy.Basic]]
)

[shef4]

Hi @NoureldinYosri, I was wondering if it's fine to close this issue unless I misunderstood the goal for part 3?