Fidelity guard hits 49qubit circuit (Distance 5 surface QEC in T+ state)

[sergei-mironov]

Could you please advise me on running a 49qubit logical T+ circuit in the tensornetwork mode? The code below emits the QUnit fidelity estimate is effectively 0! warning and it is poor indeed if the fidelity guard is disabled.

A similar program for the distance 3 code shows correct results and nice timings in this mode, so I am very interested in making running a larger version too. Please let me know what can I try.

Package versions are: pennylane 0.43.1, pennylane_catalyst 0.13.0, pennylane_qrack 0.25.0 (compiled from source), isTensorNetwork flag is enabled.

---

The full runnable program:

import json
import time
from contextlib import contextmanager
from signal import ITIMER_REAL, setitimer
from functools import reduce

import pennylane as qml
from catalyst import qjit, measure as catalyst_measure

num_shots = 200
nattempts = 1
timeout = 3000.0
dry_run = False


@contextmanager
def with_alarm(timeout: float):
    timer_set = False
    try:
        if timeout is not None and 0 < timeout < float("inf"):
            setitimer(ITIMER_REAL, timeout)
            timer_set = True
        yield
    finally:
        if timer_set:
            setitimer(ITIMER_REAL, 0)


def dump_logical_observables(measurement_attempts) -> None:
    """ For each observable, calculate its value from measurements (by xor-ing them) and dump to the
    configured JSON file.  """
    logical_values_attempts = []
    logical_measurements = {'q1_X_0': ('c_(0, 0, 0)_0', 'c_(1, 0, 0)_0', 'c_(2, 0, 0)_0', 'c_(3, 0, 0)_0', 'c_(4, 0, 0)_0', 'c_(1, 2, 1)_0', 'c_(2, 1, 1)_0', 'c_(3, 2, 1)_0', 'c_(4, 1, 1)_0', 'c_(0, 1, 1)_0', 'c_(5, 2, 1)_0')}

    for attempt_shots in measurement_attempts:
        if not attempt_shots:
            continue
        nshots = len(attempt_shots)
        logical_values = {}

        for lol, meas_labels in logical_measurements.items():
            lo_shots = []
            for i in range(nshots):
                shot_dict = attempt_shots[i]
                samples = [int(shot_dict[l]) for l in meas_labels]
                lo_value = reduce(lambda a,b: a ^ b, samples)
                lo_shots.append(int(lo_value))
            logical_values[lol] = lo_shots

        logical_values_attempts.append(logical_values)

    with open("logical_measurements.json", "w") as f:
        json.dump(logical_values_attempts, f)

if not dry_run:
    @qjit
    @qml.qnode(qml.device("qrack.simulator",wires=49,isTensorNetwork=True,isStabilizerHybrid=True,isOpenCL=True,isPaged=False,isCpuGpuHybrid=False,), shots=num_shots)
    def pl_main():

        # PennyLane Program Generated from Eka
        # Initialize Wires from Eka quantum channels Ids
        measurements = {}
        qml.registers({'20cd9525-88d4-4f79-b0d6-69896659d801': 1, '4c114191-2117-4026-98d6-3c79c90aad94': 1, '5fe8fcfe-24d7-4b0e-9160-c1b359ac4b3e': 1, '17cc13ef-3f1b-4359-bce9-1e22f1e18d2b': 1, '7583a83f-e428-4e35-a0ea-c2d94ddb6c1b': 1, 'd98a1861-4cc6-40f2-9aa7-36d8c7c05971': 1, '2a268fb7-af6f-4710-9a78-2f00ec6ccca8': 1, 'bc194ed6-3130-4f08-a66e-9ec86012437a': 1, 'e33d8911-19e9-4d2f-b9bd-76ae59655d82': 1, '4a002fb1-3b77-4441-b328-d86c037afb04': 1, '44e83a60-a55f-43ea-a182-d5e6ced3d3ec': 1, '8dd6e03b-13fe-43eb-8b1c-19aecded6131': 1, 'aed16e51-d9a9-4137-acf3-7ea4e264d21d': 1, 'b88ff4d7-6a7c-4e3c-9081-958b8e7b0052': 1, '7c7692d8-1462-4ac0-b905-93f9a503d5e9': 1, '6788af4c-aa90-4930-8977-038f899ba40f': 1, '10b37912-72a6-447e-b6bb-906fce07f686': 1, 'c927451d-b188-416a-b3a7-6e94a4e4517e': 1, '1c315269-f3a7-4d53-9f2f-cd9df110048d': 1, '808e0aef-21a5-41ba-b765-2ce0edec2d73': 1, 'c3c8ec05-3c5c-4c08-9a98-45c19d8c6679': 1, 'e1f462bf-407b-4a9e-8960-dfe06263ca40': 1, '647804f4-829b-4531-a47a-a9301181b966': 1, 'd02b3986-77c1-4227-b083-8baef94f5200': 1, 'dd77ba9f-072d-4922-8ae9-b0f6555fb09b': 1, 'a429cce7-3c33-4233-86d0-2855d24bc043': 1, 'a8736677-f3b7-4d9d-be5e-480cd4f915e6': 1, 'def4773d-0891-4ae9-8eb4-ad4e3c025748': 1, '79f9cc85-b749-4d4f-bc21-907bf8a4e32c': 1, '3779c449-f260-45ec-b793-b6482da82ba7': 1, '435ee73a-d126-4dc6-987a-8dff11f177a9': 1, '4c66ae56-f120-446d-96f9-7900fe71ce37': 1, 'fe065b07-a59a-4019-bcbe-f2d10531c282': 1, 'c68e755f-dcc0-4884-9937-d95dc0427d5c': 1, '5f646221-9634-4b07-8365-5a8c9fa3d19d': 1, '8658056b-df30-442d-84fc-4fbc08c62bac': 1, '91f930ec-81e2-48c3-8c36-35ada1ccc276': 1, '865330de-b339-423c-8fbd-ebb4ea6c00eb': 1, 'faed7ac9-ef9e-4622-ad7a-6986adbaa5a6': 1, '735f0add-4d65-4a89-9f17-ea88b17aa5af': 1, '6bb6ace0-998a-4b57-bf2e-25a14abaaa34': 1, '89bcd509-8b8e-4b7f-81b7-d458cef8d7b9': 1, 'aaf26ee1-f4d9-4ba0-8d9c-c51b1d14b422': 1, 'c13fcf73-ef5f-4fac-bfca-e9266837afdf': 1, '97e70477-6c86-401b-9b43-165aa528087c': 1, 'fd103e3c-7411-41c3-ba53-a718de15ea20': 1, 'd1e36d0d-3114-49ea-b8c5-9c5536fc28e3': 1, '0a84f1f2-d755-47d6-9570-0ed03f0108b3': 1, '12a5a18e-6568-43b5-b622-8da3d17320b7': 1})
        
        # Start of circuit: final circuit
        # Start of circuit: inject t into block q1 and measure syndromes
        # Start of circuit: inject t into block q1
        # Start of circuit: resource state reset
        catalyst_measure(21, reset=True)
        qml.Hadamard([21])
        qml.T([21])
        #  End of circuit: resource state reset
        # Start of circuit: reset four quadrants
        catalyst_measure(8, reset=True)
        qml.Hadamard([8])
        catalyst_measure(1, reset=True)
        qml.Hadamard([1])
        catalyst_measure(10, reset=True)
        qml.Hadamard([10])
        catalyst_measure(3, reset=True)
        qml.Hadamard([3])
        catalyst_measure(0, reset=True)
        qml.Hadamard([0])
        catalyst_measure(4, reset=True)
        qml.Hadamard([4])
        catalyst_measure(32, reset=True)
        qml.Hadamard([32])
        catalyst_measure(39, reset=True)
        qml.Hadamard([39])
        catalyst_measure(41, reset=True)
        qml.Hadamard([41])
        catalyst_measure(43, reset=True)
        qml.Hadamard([43])
        catalyst_measure(45, reset=True)
        qml.Hadamard([45])
        catalyst_measure(30, reset=True)
        qml.Hadamard([30])
        catalyst_measure(27, reset=True)
        catalyst_measure(19, reset=True)
        catalyst_measure(28, reset=True)
        catalyst_measure(37, reset=True)
        catalyst_measure(7, reset=True)
        catalyst_measure(17, reset=True)
        catalyst_measure(25, reset=True)
        catalyst_measure(34, reset=True)
        catalyst_measure(12, reset=True)
        catalyst_measure(14, reset=True)
        catalyst_measure(6, reset=True)
        catalyst_measure(23, reset=True)
        #  End of circuit: reset four quadrants
        #  End of circuit: inject t into block q1
        # Start of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure q1 syndromes - cycle 0
        catalyst_measure(9, reset=True)
        catalyst_measure(13, reset=True)
        catalyst_measure(22, reset=True)
        catalyst_measure(26, reset=True)
        catalyst_measure(29, reset=True)
        catalyst_measure(33, reset=True)
        catalyst_measure(42, reset=True)
        catalyst_measure(46, reset=True)
        catalyst_measure(11, reset=True)
        catalyst_measure(15, reset=True)
        catalyst_measure(20, reset=True)
        catalyst_measure(24, reset=True)
        catalyst_measure(31, reset=True)
        catalyst_measure(35, reset=True)
        catalyst_measure(40, reset=True)
        catalyst_measure(44, reset=True)
        catalyst_measure(2, reset=True)
        catalyst_measure(5, reset=True)
        catalyst_measure(47, reset=True)
        catalyst_measure(48, reset=True)
        catalyst_measure(18, reset=True)
        catalyst_measure(38, reset=True)
        catalyst_measure(16, reset=True)
        catalyst_measure(36, reset=True)
        qml.Hadamard([9])
        qml.Hadamard([13])
        qml.Hadamard([22])
        qml.Hadamard([26])
        qml.Hadamard([29])
        qml.Hadamard([33])
        qml.Hadamard([42])
        qml.Hadamard([46])
        qml.Hadamard([11])
        qml.Hadamard([15])
        qml.Hadamard([20])
        qml.Hadamard([24])
        qml.Hadamard([31])
        qml.Hadamard([35])
        qml.Hadamard([40])
        qml.Hadamard([44])
        qml.Hadamard([2])
        qml.Hadamard([5])
        qml.Hadamard([47])
        qml.Hadamard([48])
        qml.Hadamard([18])
        qml.Hadamard([38])
        qml.Hadamard([16])
        qml.Hadamard([36])
        qml.CZ([9, 7])
        qml.CZ([13, 10])
        qml.CZ([22, 19])
        qml.CZ([26, 23])
        qml.CZ([29, 27])
        qml.CZ([33, 30])
        qml.CZ([42, 39])
        qml.CZ([46, 43])
        qml.CNOT([11, 8])
        qml.CNOT([15, 12])
        qml.CNOT([20, 17])
        qml.CNOT([24, 21])
        qml.CNOT([31, 28])
        qml.CNOT([35, 32])
        qml.CNOT([40, 37])
        qml.CNOT([44, 41])
        qml.CNOT([2, 0])
        qml.CNOT([5, 3])
        qml.CZ([16, 14])
        qml.CZ([36, 34])
        qml.CZ([9, 0])
        qml.CZ([13, 3])
        qml.CZ([22, 8])
        qml.CZ([26, 12])
        qml.CZ([29, 17])
        qml.CZ([33, 21])
        qml.CZ([42, 28])
        qml.CZ([46, 32])
        qml.CNOT([11, 10])
        qml.CNOT([15, 14])
        qml.CNOT([20, 19])
        qml.CNOT([24, 23])
        qml.CNOT([31, 30])
        qml.CNOT([35, 34])
        qml.CNOT([40, 39])
        qml.CNOT([44, 43])
        qml.CNOT([2, 1])
        qml.CNOT([5, 4])
        qml.CZ([16, 6])
        qml.CZ([36, 25])
        qml.CZ([9, 8])
        qml.CZ([13, 12])
        qml.CZ([22, 21])
        qml.CZ([26, 25])
        qml.CZ([29, 28])
        qml.CZ([33, 32])
        qml.CZ([42, 41])
        qml.CZ([46, 45])
        qml.CNOT([11, 1])
        qml.CNOT([15, 4])
        qml.CNOT([20, 7])
        qml.CNOT([24, 10])
        qml.CNOT([31, 19])
        qml.CNOT([35, 23])
        qml.CNOT([40, 27])
        qml.CNOT([44, 30])
        qml.CNOT([47, 39])
        qml.CNOT([48, 43])
        qml.CZ([18, 17])
        qml.CZ([38, 37])
        qml.CZ([9, 1])
        qml.CZ([13, 4])
        qml.CZ([22, 10])
        qml.CZ([26, 14])
        qml.CZ([29, 19])
        qml.CZ([33, 23])
        qml.CZ([42, 30])
        qml.CZ([46, 34])
        qml.CNOT([11, 3])
        qml.CNOT([15, 6])
        qml.CNOT([20, 8])
        qml.CNOT([24, 12])
        qml.CNOT([31, 21])
        qml.CNOT([35, 25])
        qml.CNOT([40, 28])
        qml.CNOT([44, 32])
        qml.CNOT([47, 41])
        qml.CNOT([48, 45])
        qml.CZ([18, 7])
        qml.CZ([38, 27])
        qml.Hadamard([9])
        qml.Hadamard([13])
        qml.Hadamard([22])
        qml.Hadamard([26])
        qml.Hadamard([29])
        qml.Hadamard([33])
        qml.Hadamard([42])
        qml.Hadamard([46])
        qml.Hadamard([11])
        qml.Hadamard([15])
        qml.Hadamard([20])
        qml.Hadamard([24])
        qml.Hadamard([31])
        qml.Hadamard([35])
        qml.Hadamard([40])
        qml.Hadamard([44])
        qml.Hadamard([2])
        qml.Hadamard([5])
        qml.Hadamard([47])
        qml.Hadamard([48])
        qml.Hadamard([18])
        qml.Hadamard([38])
        qml.Hadamard([16])
        qml.Hadamard([36])
        measurements["c_(1, 1, 1)_0"] = catalyst_measure(9)
        measurements["c_(1, 3, 1)_0"] = catalyst_measure(13)
        measurements["c_(2, 2, 1)_0"] = catalyst_measure(22)
        measurements["c_(2, 4, 1)_0"] = catalyst_measure(26)
        measurements["c_(3, 1, 1)_0"] = catalyst_measure(29)
        measurements["c_(3, 3, 1)_0"] = catalyst_measure(33)
        measurements["c_(4, 2, 1)_0"] = catalyst_measure(42)
        measurements["c_(4, 4, 1)_0"] = catalyst_measure(46)
        measurements["c_(1, 2, 1)_0"] = catalyst_measure(11)
        measurements["c_(1, 4, 1)_0"] = catalyst_measure(15)
        measurements["c_(2, 1, 1)_0"] = catalyst_measure(20)
        measurements["c_(2, 3, 1)_0"] = catalyst_measure(24)
        measurements["c_(3, 2, 1)_0"] = catalyst_measure(31)
        measurements["c_(3, 4, 1)_0"] = catalyst_measure(35)
        measurements["c_(4, 1, 1)_0"] = catalyst_measure(40)
        measurements["c_(4, 3, 1)_0"] = catalyst_measure(44)
        measurements["c_(0, 1, 1)_0"] = catalyst_measure(2)
        measurements["c_(0, 3, 1)_0"] = catalyst_measure(5)
        measurements["c_(5, 2, 1)_0"] = catalyst_measure(47)
        measurements["c_(5, 4, 1)_0"] = catalyst_measure(48)
        measurements["c_(2, 0, 1)_0"] = catalyst_measure(18)
        measurements["c_(4, 0, 1)_0"] = catalyst_measure(38)
        measurements["c_(1, 5, 1)_0"] = catalyst_measure(16)
        measurements["c_(3, 5, 1)_0"] = catalyst_measure(36)
        #  End of circuit: measure q1 syndromes - cycle 0
        #  End of circuit: measure q1 syndromes 1 time(s)
        #  End of circuit: inject t into block q1 and measure syndromes
        # Start of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure q1 syndromes - cycle 0
        catalyst_measure(9, reset=True)
        catalyst_measure(13, reset=True)
        catalyst_measure(22, reset=True)
        catalyst_measure(26, reset=True)
        catalyst_measure(29, reset=True)
        catalyst_measure(33, reset=True)
        catalyst_measure(42, reset=True)
        catalyst_measure(46, reset=True)
        catalyst_measure(11, reset=True)
        catalyst_measure(15, reset=True)
        catalyst_measure(20, reset=True)
        catalyst_measure(24, reset=True)
        catalyst_measure(31, reset=True)
        catalyst_measure(35, reset=True)
        catalyst_measure(40, reset=True)
        catalyst_measure(44, reset=True)
        catalyst_measure(2, reset=True)
        catalyst_measure(5, reset=True)
        catalyst_measure(47, reset=True)
        catalyst_measure(48, reset=True)
        catalyst_measure(18, reset=True)
        catalyst_measure(38, reset=True)
        catalyst_measure(16, reset=True)
        catalyst_measure(36, reset=True)
        qml.Hadamard([9])
        qml.Hadamard([13])
        qml.Hadamard([22])
        qml.Hadamard([26])
        qml.Hadamard([29])
        qml.Hadamard([33])
        qml.Hadamard([42])
        qml.Hadamard([46])
        qml.Hadamard([11])
        qml.Hadamard([15])
        qml.Hadamard([20])
        qml.Hadamard([24])
        qml.Hadamard([31])
        qml.Hadamard([35])
        qml.Hadamard([40])
        qml.Hadamard([44])
        qml.Hadamard([2])
        qml.Hadamard([5])
        qml.Hadamard([47])
        qml.Hadamard([48])
        qml.Hadamard([18])
        qml.Hadamard([38])
        qml.Hadamard([16])
        qml.Hadamard([36])
        qml.CZ([9, 7])
        qml.CZ([13, 10])
        qml.CZ([22, 19])
        qml.CZ([26, 23])
        qml.CZ([29, 27])
        qml.CZ([33, 30])
        qml.CZ([42, 39])
        qml.CZ([46, 43])
        qml.CNOT([11, 8])
        qml.CNOT([15, 12])
        qml.CNOT([20, 17])
        qml.CNOT([24, 21])
        qml.CNOT([31, 28])
        qml.CNOT([35, 32])
        qml.CNOT([40, 37])
        qml.CNOT([44, 41])
        qml.CNOT([2, 0])
        qml.CNOT([5, 3])
        qml.CZ([16, 14])
        qml.CZ([36, 34])
        qml.CZ([9, 0])
        qml.CZ([13, 3])
        qml.CZ([22, 8])
        qml.CZ([26, 12])
        qml.CZ([29, 17])
        qml.CZ([33, 21])
        qml.CZ([42, 28])
        qml.CZ([46, 32])
        qml.CNOT([11, 10])
        qml.CNOT([15, 14])
        qml.CNOT([20, 19])
        qml.CNOT([24, 23])
        qml.CNOT([31, 30])
        qml.CNOT([35, 34])
        qml.CNOT([40, 39])
        qml.CNOT([44, 43])
        qml.CNOT([2, 1])
        qml.CNOT([5, 4])
        qml.CZ([16, 6])
        qml.CZ([36, 25])
        qml.CZ([9, 8])
        qml.CZ([13, 12])
        qml.CZ([22, 21])
        qml.CZ([26, 25])
        qml.CZ([29, 28])
        qml.CZ([33, 32])
        qml.CZ([42, 41])
        qml.CZ([46, 45])
        qml.CNOT([11, 1])
        qml.CNOT([15, 4])
        qml.CNOT([20, 7])
        qml.CNOT([24, 10])
        qml.CNOT([31, 19])
        qml.CNOT([35, 23])
        qml.CNOT([40, 27])
        qml.CNOT([44, 30])
        qml.CNOT([47, 39])
        qml.CNOT([48, 43])
        qml.CZ([18, 17])
        qml.CZ([38, 37])
        qml.CZ([9, 1])
        qml.CZ([13, 4])
        qml.CZ([22, 10])
        qml.CZ([26, 14])
        qml.CZ([29, 19])
        qml.CZ([33, 23])
        qml.CZ([42, 30])
        qml.CZ([46, 34])
        qml.CNOT([11, 3])
        qml.CNOT([15, 6])
        qml.CNOT([20, 8])
        qml.CNOT([24, 12])
        qml.CNOT([31, 21])
        qml.CNOT([35, 25])
        qml.CNOT([40, 28])
        qml.CNOT([44, 32])
        qml.CNOT([47, 41])
        qml.CNOT([48, 45])
        qml.CZ([18, 7])
        qml.CZ([38, 27])
        qml.Hadamard([9])
        qml.Hadamard([13])
        qml.Hadamard([22])
        qml.Hadamard([26])
        qml.Hadamard([29])
        qml.Hadamard([33])
        qml.Hadamard([42])
        qml.Hadamard([46])
        qml.Hadamard([11])
        qml.Hadamard([15])
        qml.Hadamard([20])
        qml.Hadamard([24])
        qml.Hadamard([31])
        qml.Hadamard([35])
        qml.Hadamard([40])
        qml.Hadamard([44])
        qml.Hadamard([2])
        qml.Hadamard([5])
        qml.Hadamard([47])
        qml.Hadamard([48])
        qml.Hadamard([18])
        qml.Hadamard([38])
        qml.Hadamard([16])
        qml.Hadamard([36])
        measurements["c_(1, 1, 1)_1"] = catalyst_measure(9)
        measurements["c_(1, 3, 1)_1"] = catalyst_measure(13)
        measurements["c_(2, 2, 1)_1"] = catalyst_measure(22)
        measurements["c_(2, 4, 1)_1"] = catalyst_measure(26)
        measurements["c_(3, 1, 1)_1"] = catalyst_measure(29)
        measurements["c_(3, 3, 1)_1"] = catalyst_measure(33)
        measurements["c_(4, 2, 1)_1"] = catalyst_measure(42)
        measurements["c_(4, 4, 1)_1"] = catalyst_measure(46)
        measurements["c_(1, 2, 1)_1"] = catalyst_measure(11)
        measurements["c_(1, 4, 1)_1"] = catalyst_measure(15)
        measurements["c_(2, 1, 1)_1"] = catalyst_measure(20)
        measurements["c_(2, 3, 1)_1"] = catalyst_measure(24)
        measurements["c_(3, 2, 1)_1"] = catalyst_measure(31)
        measurements["c_(3, 4, 1)_1"] = catalyst_measure(35)
        measurements["c_(4, 1, 1)_1"] = catalyst_measure(40)
        measurements["c_(4, 3, 1)_1"] = catalyst_measure(44)
        measurements["c_(0, 1, 1)_1"] = catalyst_measure(2)
        measurements["c_(0, 3, 1)_1"] = catalyst_measure(5)
        measurements["c_(5, 2, 1)_1"] = catalyst_measure(47)
        measurements["c_(5, 4, 1)_1"] = catalyst_measure(48)
        measurements["c_(2, 0, 1)_1"] = catalyst_measure(18)
        measurements["c_(4, 0, 1)_1"] = catalyst_measure(38)
        measurements["c_(1, 5, 1)_1"] = catalyst_measure(16)
        measurements["c_(3, 5, 1)_1"] = catalyst_measure(36)
        #  End of circuit: measure q1 syndromes - cycle 0
        #  End of circuit: measure q1 syndromes 1 time(s)
        # Start of circuit: measure logical x of q1
        qml.Hadamard([0])
        qml.Hadamard([1])
        qml.Hadamard([3])
        qml.Hadamard([4])
        qml.Hadamard([6])
        qml.Hadamard([7])
        qml.Hadamard([8])
        qml.Hadamard([10])
        qml.Hadamard([12])
        qml.Hadamard([14])
        qml.Hadamard([17])
        qml.Hadamard([19])
        qml.Hadamard([21])
        qml.Hadamard([23])
        qml.Hadamard([25])
        qml.Hadamard([27])
        qml.Hadamard([28])
        qml.Hadamard([30])
        qml.Hadamard([32])
        qml.Hadamard([34])
        qml.Hadamard([37])
        qml.Hadamard([39])
        qml.Hadamard([41])
        qml.Hadamard([43])
        qml.Hadamard([45])
        measurements["c_(0, 0, 0)_0"] = catalyst_measure(0)
        measurements["c_(0, 1, 0)_0"] = catalyst_measure(1)
        measurements["c_(0, 2, 0)_0"] = catalyst_measure(3)
        measurements["c_(0, 3, 0)_0"] = catalyst_measure(4)
        measurements["c_(0, 4, 0)_0"] = catalyst_measure(6)
        measurements["c_(1, 0, 0)_0"] = catalyst_measure(7)
        measurements["c_(1, 1, 0)_0"] = catalyst_measure(8)
        measurements["c_(1, 2, 0)_0"] = catalyst_measure(10)
        measurements["c_(1, 3, 0)_0"] = catalyst_measure(12)
        measurements["c_(1, 4, 0)_0"] = catalyst_measure(14)
        measurements["c_(2, 0, 0)_0"] = catalyst_measure(17)
        measurements["c_(2, 1, 0)_0"] = catalyst_measure(19)
        measurements["c_(2, 2, 0)_0"] = catalyst_measure(21)
        measurements["c_(2, 3, 0)_0"] = catalyst_measure(23)
        measurements["c_(2, 4, 0)_0"] = catalyst_measure(25)
        measurements["c_(3, 0, 0)_0"] = catalyst_measure(27)
        measurements["c_(3, 1, 0)_0"] = catalyst_measure(28)
        measurements["c_(3, 2, 0)_0"] = catalyst_measure(30)
        measurements["c_(3, 3, 0)_0"] = catalyst_measure(32)
        measurements["c_(3, 4, 0)_0"] = catalyst_measure(34)
        measurements["c_(4, 0, 0)_0"] = catalyst_measure(37)
        measurements["c_(4, 1, 0)_0"] = catalyst_measure(39)
        measurements["c_(4, 2, 0)_0"] = catalyst_measure(41)
        measurements["c_(4, 3, 0)_0"] = catalyst_measure(43)
        measurements["c_(4, 4, 0)_0"] = catalyst_measure(45)
        #  End of circuit: measure logical x of q1
        #  End of circuit: final circuit
        return measurements


else:
    pl_main = None

if not dry_run:
    with with_alarm(timeout or float("inf")):
        pl_main.compile()

times = []
measurements = []

print("Running", end="")
for _ in range(nattempts):
    begin = time.time()
    shot_results = []
    if not dry_run:
        with with_alarm(timeout or float("inf")):
            raw = pl_main()
            per_key_lists = {
                k: v.astype(int).tolist() for k, v in raw.items()
            }
            nshots = (
                len(next(iter(per_key_lists.values())))
                if per_key_lists else 0
            )
            for i in range(nshots):
                shot_results.append({
                    k: vals[i] for k, vals in per_key_lists.items()
                })
    end = time.time()
    times.append(end - begin)
    measurements.append(shot_results)

    with open("running_times.json", "w") as f:
        json.dump(times, f)
    with open("measurements.json", "w") as f:
        json.dump(measurements, f)

    print(f", t={times[-1]}", end="", flush=True)
print(" (dry-run)" if dry_run else "")

if not dry_run:
    dump_logical_observables(measurements)

[WrathfulSpatula]

My guess as to what's happening is this (although the emergent behavior of Qrack can sometimes be hard to properly intuit, without significant consideration and effort): seeing as it looks like this only has a single T gate, for any non-Clifford gates (if I have that right), it should be a candidate for near-Clifford simulation, but Qrack needs to stay conservative, as a JIT state machine, that any further gates applied might take it over any automatically detected limit for system memory footprint. Hence, it's probably detecting that your system might support something like roughly 30 qubits of state vector (probably), but you're requesting 49 plus 1 T gate (which adds an ancillary qubit), and it can't be sure ahead of time that your circuit will stay near-Clifford, so it proactively chooses "automatic circuit elision" (ACE) mode over near-Clifford.

There's a way to force it to basically always use near-Clifford mode, if your circuit is Clifford+RZ (or Clifford+T):

1. Set environment variable "export QRACK_MAX_PAGING_QB=-1. (Qrack reads the -1 as wrap-around to maximum integer, so it effectively removes any limit.)
2. Set environment variable "export QRACK_MAX_CPU_QB=-1. (Same idea.)
3. @qml.qnode(qml.device("qrack.simulator",wires=49, isTensorNetwork=False, isSchmidtDecompose=False, isStabilizerHybrid=True), shots=num_shots) (isTensorNetwork optimizes the circuit representation before dispatching to the appropriate state machine, ultimately, and it might work whether True or False in this case, but that's based on some of the most recent updates, and I haven't tested it for your case, yet.)

With these settings, I've had your script running for about the past ten minutes: there's no longer any QUnit fidelity error, but what the memory and processor utilization seem to say me is that the unitary portion including T-gate "gadgetization" runs in something like a fraction of a second (which is always the easy part, for our near-Clifford algorithm), then multi-threaded marginal probability calculation for measurement sampling is a "monster," for 49 qubits. I can't promise you that effectively 49+1 qubits worth of marginal probability calculation for measurement sampling will complete in reasonable time, but, based on my own test, just now, I can tell that the working memory footprint will be trivial, and it'll fix the behavior you didn't want where Qrack tries ACE instead, for safety's sake.

[WrathfulSpatula]

On top of that, I look at your circuit again, and I realize it has multiple rounds of measurements, which is extremely costly. So, Qrack's near-Clifford is different from IBM's. I came up with the idea for it a few years ago based on Appendix A of 10.1103/PRXQuantum.3.020361. Using their "reverse phase-injection gadget," you can handle Clifford+RZ in a way that is exact, with post-selection on the gadget qubits, just like they describe. For amplitude or marginal probability calculation, you end up needing to run an auxiliary state-vector simulator just over the gadget ancillae, so your memory is effectively limited by power-of-2 growth over non-Clifford phase gate count. This is great in a case like when you want to calculate only a very small number of amplitudes or basis-dimension marginal probabilities for large near-Clifford circuits. Similarly, if the overall count of nonzero stabilizer amplitudes is small, you can basically get away with something like a small constant factor on a single traversal of a sparse set amplitudes that would fit in a laptop's memory at once anyway (to my recall, at least). But the problem with sampling is that, while this will all fit in memory perfectly well, even if the stabilizer amplitude count isn't sparse, so long as the ancillary qubit count is small, Qrack still doesn't have a better way to sample (for now) than to traverse marginal probabilities by basis dimension until it exceeds the probability threshold for whatever random [0.0, 1.0) number it pulled via random number generator for a sample. It even batches this when you call the multi-shot measurement methods, so it only needs to traverse the marginals calculation once for multiple shots, worth of random floats, but, if the circuit isn't sparse, that's still, like, ~2^49 marginals to traverse.

We can't offer IBM's (approximate, stablizer-rank-based) method, which will be appropriate for different use cases, but neither do we seek to compete with it: they're in a different niche, with a different algorithm, and they might work better for different use cases. I'm saying, keep their method in mind, alternatively. But, with Qrack, you'll fit an exact calculation in memory, but sampling and mid-circuit measurement is extremely costly.

[WrathfulSpatula]

I just hit your timeout condition, in fact. But memory usage was trivial, and near-Clifford was used instead of ACE.

[WrathfulSpatula]

Sorry to spam, but I'll just add one more thing: I don't suspect that ACE is an appropriate alternative for your case, but it will at least complete, quickly, with approximate simulation. The exception you saw it raise was that ACE tracks an internal first-principles, worst-case fidelity estimator based on the entangling gates it replaces with semi-classical gate shadows, to artificially truncate long-range entanglement. However, this estimator can be extremely pessimistic when it comes to something like RCS or even QFT-based algorithms. The exception is raised by default exactly so that people simply don't stumble into and misunderstand what I'm about show you, when intentionally trying to use ACE approximation. But how you would typically use ACE, in general, is to disable the fidelity estimator guard and tell Qrack to split your circuit into about 2 or 4 "patches" without long-range entanglement between them:

1. export QRACK_DISABLE_QUNIT_FIDELITY_GUARD=1
2. export QRACK_MAX_PAGING_QB=25 (for ~2 patches on 49 qubits)
3. @qml.qnode(qml.device("qrack.simulator",wires=49), shots=num_shots)
4. Run as normal. (But validating the approximation is on you. That's exactly why we include the "guard," so nothing is accidentally misinterpreted.)

[WrathfulSpatula]

Actually, you asked about isTensorNetwork, so let me just take the opportunity, for future reference of other users, to round-out an explanation of what that does and all the major alternatives we offer that might suit any specific circuit case.

isTensorNetwork caches and optimizes a circuit diagram, including with light-cone-based optimizations when local measurements or observables are requested. We call it isTensorNetwork because it's working on the circuit diagram rather than Schrödinger picture, not entirely unlike a tensor network tree. However, what it's not doing is replacing the role of software like quimb to just automatically pick optimal contraction paths in a conventional tensor network tree. (Use quimb, for that! Similarly, if one wants IBM's approximate stabilizer-ranked-based near-Clifford method, use Aer! We differentiate and work together, as qrack.) However, qrack is a JIT state machine. After isTensorNetwork does its local circuit simplification and light-cone tricks, that's the interface contract Qrack by design and, therefore, downstream users are asking from it. So it ultimately must dispatch to a Schrödinger picture simulator, finally (at least from the simulator options available to it in Qrack). Light-cone optimization and local circuit simplification can do a lot of work in specific cases, but, if you want conventional tensor network trees, look to other (complementary) software in the ecosystem like quimb.

Qrack does its best to try pick the right tool for you, transparently, from its own (differentiated) toolbox. That only goes so far, though. Near-Clifford, ACE, quantum binary decision diagrams (QBDD), or sparse simulation all have trade-offs that require a degree of thought and intentionality from the user, to first understand what they're actually trying simulate, what observables they want back, and which approximations are acceptable.

For circuits that produce intermediate states that are regular in ways amenable to QBDD (see arXiv:2302.04687), keep that method in mind. It also tends to work well for at least a few shallow layers of high-entropy intermediate states, even for RCS.

There are many other little rounding-tolerance knobs to tune via environment variable, to trade accuracy for performance, from the Qrack README, but, if none of the above are actually appropriate to your case (and the most obvious but semi-trivial such example is fully-connected RCS, in my own experience), last line of defense might be sparse (ACE) simulation:

1. export QRACK_SPARSE_MAX_ALLOC_MB=1 (Start here, with just 1 MB in footprint, and work upwards if necessary. I know that's tiny, but, besides the method being CPU-based and potentially a bit slow, it's surprising how much fidelity you often get in the end while retaining even only the top 512 KB to 1 MB of most significant amplitude components).
2. @qml.qnode(qml.device("qrack.simulator",wires=49, isOpenCL=False, isSparse=True), shots=num_shots)
3. Validate on small-scale circuits you can compare against ideal, like you also have to do with general ACE, to build a model of how it will extrapolate, because that's just the territory the field is in, mostly, with approximating NISQ-scale circuits.

But I'll let you digest this all, now. Sorry to bombard you with all this information at once, but the issue will be useful as a reference for other users.