diff --git a/CMakeLists.txt b/CMakeLists.txt index 80ccd8585..805ce3665 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,6 +81,7 @@ set(SOURCE_FILES_NO_MAIN src/stim/search/hyper/graph.cc src/stim/search/hyper/node.cc src/stim/search/hyper/search_state.cc + src/stim/simulators/dem_sampler.cc src/stim/simulators/detection_simulator.cc src/stim/simulators/error_analyzer.cc src/stim/simulators/error_matcher.cc @@ -136,6 +137,7 @@ set(TEST_FILES src/stim/search/hyper/graph.test.cc src/stim/search/hyper/node.test.cc src/stim/search/hyper/search_state.test.cc + src/stim/simulators/dem_sampler.test.cc src/stim/simulators/detection_simulator.test.cc src/stim/simulators/error_analyzer.test.cc src/stim/simulators/error_matcher.test.cc @@ -166,6 +168,7 @@ set(BENCHMARK_FILES src/stim/mem/sparse_xor_vec.perf.cc src/stim/probability_util.perf.cc src/stim/search/graphlike/algo.perf.cc + src/stim/simulators/dem_sampler.perf.cc src/stim/simulators/error_analyzer.perf.cc src/stim/simulators/frame_simulator.perf.cc src/stim/simulators/tableau_simulator.perf.cc diff --git a/doc/python_api_reference_vDev.md b/doc/python_api_reference_vDev.md index 00118828f..6baf74b8e 100644 --- a/doc/python_api_reference_vDev.md +++ b/doc/python_api_reference_vDev.md @@ -20,7 +20,6 @@ API references for stable versions are kept on the [stim github wiki](https://gi - [`stim.Circuit.__str__`](#stim.Circuit.__str__) - [`stim.Circuit.append`](#stim.Circuit.append) - [`stim.Circuit.append_from_stim_program_text`](#stim.Circuit.append_from_stim_program_text) - - [`stim.Circuit.append_operation`](#stim.Circuit.append_operation) - [`stim.Circuit.approx_equals`](#stim.Circuit.approx_equals) - [`stim.Circuit.clear`](#stim.Circuit.clear) - [`stim.Circuit.compile_detector_sampler`](#stim.Circuit.compile_detector_sampler) @@ -30,7 +29,6 @@ API references for stable versions are kept on the [stim github wiki](https://gi - [`stim.Circuit.detector_error_model`](#stim.Circuit.detector_error_model) - [`stim.Circuit.explain_detector_error_model_errors`](#stim.Circuit.explain_detector_error_model_errors) - [`stim.Circuit.flattened`](#stim.Circuit.flattened) - - [`stim.Circuit.flattened_operations`](#stim.Circuit.flattened_operations) - [`stim.Circuit.from_file`](#stim.Circuit.from_file) - [`stim.Circuit.generated`](#stim.Circuit.generated) - [`stim.Circuit.get_detector_coordinates`](#stim.Circuit.get_detector_coordinates) @@ -79,6 +77,9 @@ API references for stable versions are kept on the [stim github wiki](https://gi - [`stim.CircuitTargetsInsideInstruction.target_range_end`](#stim.CircuitTargetsInsideInstruction.target_range_end) - [`stim.CircuitTargetsInsideInstruction.target_range_start`](#stim.CircuitTargetsInsideInstruction.target_range_start) - [`stim.CircuitTargetsInsideInstruction.targets_in_range`](#stim.CircuitTargetsInsideInstruction.targets_in_range) +- [`stim.CompiledDemSampler`](#stim.CompiledDemSampler) + - [`stim.CompiledDemSampler.sample`](#stim.CompiledDemSampler.sample) + - [`stim.CompiledDemSampler.sample_write`](#stim.CompiledDemSampler.sample_write) - [`stim.CompiledDetectorSampler`](#stim.CompiledDetectorSampler) - [`stim.CompiledDetectorSampler.__init__`](#stim.CompiledDetectorSampler.__init__) - [`stim.CompiledDetectorSampler.__repr__`](#stim.CompiledDetectorSampler.__repr__) @@ -144,12 +145,15 @@ API references for stable versions are kept on the [stim github wiki](https://gi - [`stim.DetectorErrorModel.append`](#stim.DetectorErrorModel.append) - [`stim.DetectorErrorModel.approx_equals`](#stim.DetectorErrorModel.approx_equals) - [`stim.DetectorErrorModel.clear`](#stim.DetectorErrorModel.clear) + - [`stim.DetectorErrorModel.compile_sampler`](#stim.DetectorErrorModel.compile_sampler) - [`stim.DetectorErrorModel.copy`](#stim.DetectorErrorModel.copy) + - [`stim.DetectorErrorModel.flattened`](#stim.DetectorErrorModel.flattened) - [`stim.DetectorErrorModel.from_file`](#stim.DetectorErrorModel.from_file) - [`stim.DetectorErrorModel.get_detector_coordinates`](#stim.DetectorErrorModel.get_detector_coordinates) - [`stim.DetectorErrorModel.num_detectors`](#stim.DetectorErrorModel.num_detectors) - [`stim.DetectorErrorModel.num_errors`](#stim.DetectorErrorModel.num_errors) - [`stim.DetectorErrorModel.num_observables`](#stim.DetectorErrorModel.num_observables) + - [`stim.DetectorErrorModel.rounded`](#stim.DetectorErrorModel.rounded) - [`stim.DetectorErrorModel.shortest_graphlike_error`](#stim.DetectorErrorModel.shortest_graphlike_error) - [`stim.DetectorErrorModel.to_file`](#stim.DetectorErrorModel.to_file) - [`stim.ExplainedError`](#stim.ExplainedError) @@ -198,7 +202,6 @@ API references for stable versions are kept on the [stim github wiki](https://gi - [`stim.PauliString.__truediv__`](#stim.PauliString.__truediv__) - [`stim.PauliString.commutes`](#stim.PauliString.commutes) - [`stim.PauliString.copy`](#stim.PauliString.copy) - - [`stim.PauliString.extended_product`](#stim.PauliString.extended_product) - [`stim.PauliString.random`](#stim.PauliString.random) - [`stim.PauliString.sign`](#stim.PauliString.sign) - [`stim.Tableau`](#stim.Tableau) @@ -781,21 +784,6 @@ def append_from_stim_program_text( """ ``` - -```python -# stim.Circuit.append_operation - -# (in class stim.Circuit) -def append_operation( - self, - name: object, - targets: object = (), - arg: object = None, -) -> None: - """[DEPRECATED] use stim.Circuit.append instead - """ -``` - ```python # stim.Circuit.approx_equals @@ -1246,38 +1234,6 @@ def flattened( """ ``` - -```python -# stim.Circuit.flattened_operations - -# (in class stim.Circuit) -def flattened_operations( - self, -) -> list: - """[DEPRECATED] - - Returns a list of tuples encoding the contents of the circuit. - Instead of this method, use `for instruction in circuit` or, to - avoid REPEAT blocks, `for instruction in circuit.flattened()`. - - Examples: - >>> import stim - >>> stim.Circuit(''' - ... H 0 - ... X_ERROR(0.125) 1 - ... M 0 !1 - ... ''').flattened_operations() - [('H', [0], 0), ('X_ERROR', [1], 0.125), ('M', [0, ('inv', 1)], 0)] - - >>> stim.Circuit(''' - ... REPEAT 2 { - ... H 6 - ... } - ... ''').flattened_operations() - [('H', [6], 0), ('H', [6], 0)] - """ -``` - ```python # stim.Circuit.from_file @@ -1285,9 +1241,14 @@ def flattened_operations( # (in class stim.Circuit) @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.Circuit: - """Args: + + """Reads a stim circuit from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + + Args: file: A file path or open file object to read from. Returns: @@ -1777,6 +1738,8 @@ def to_file( """Writes the stim circuit to a file. + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + Args: file: A file path or an open file to write to. @@ -2391,6 +2354,261 @@ def targets_in_range( """ ``` + +```python +# stim.CompiledDemSampler + +# (at top-level in the stim module) +class CompiledDemSampler: + """A helper class for efficiently sampler from a detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ +``` + + +```python +# stim.CompiledDemSampler.sample + +# (in class stim.CompiledDemSampler) +def sample( + self, + shots: int, + *, + bit_packed: bool = False, + return_errors: bool = False, + recorded_errors_to_replay: Optional[np.ndarray] = None, +) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + + """Samples the detector error model's error mechanisms to produce sample data. + + Args: + shots: The number of times to sample from the model. + bit_packed: Defaults to false. + False: the returned numpy arrays have dtype=np.bool8. + True: the returned numpy arrays have dtype=np.uint8 and pack 8 bits into each byte. + + Setting this to True is equivalent to running np.packbits(data, endian='little', axis=1) + on each output value, but has the performance benefit of the data never being expanded + into an unpacked form. + return_errors: Defaults to False. + False: the first entry of the returned tuple is None. + True: the first entry of the returned tuple is a numpy array recording which errors were sampled. + recorded_errors_to_replay: Defaults to None, meaning sample errors randomly. + If not None, this is expected to be a 2d numpy array specifying which errors to apply (e.g. one + returned from a previous call to the sample method). The array must have + dtype=np.bool8 and shape=(num_shots, num_errors) or + dtype=np.uint8 and shape=(num_shots, math.ceil(num_errors / 8)). + + Returns: + A tuple (detector_data, obs_data, error_data). + + Assuming bit_packed is False and return_errors is True: + If error_data[s, k] is True, then the error with index k fired in the shot with index s. + If detector_data[s, k] is True, then the detector with index k ended up flipped in the shot with index s. + If obs_data[s, k] is True, then the observable with index k ended up flipped in the shot with index s. + + The dtype and shape of the data depends on the arguments: + if bit_packed: + detector_data.shape == (num_shots, num_detectors) + detector_data.dtype == np.bool8 + obs_data.shape == (num_shots, num_observables) + obs_data.dtype == np.bool8 + if return_errors: + error_data.shape = (num_shots, num_errors) + error_data.dtype = np.bool8 + else: + error_data is None + else: + detector_data.shape == (num_shots, math.ceil(num_detectors / 8)) + detector_data.dtype == np.uint8 + obs_data.shape == (num_shots, math.ceil(num_observables / 8)) + obs_data.dtype == np.uint8 + if return_errors: + error_data.shape = (num_shots, math.ceil(num_errors / 8)) + error_data.dtype = np.uint8 + else: + error_data is None + + Note that bit packing is done using little endian order on the last axis + (i.e. like `np.packbits(data, endian='little', axis=1)`). + + Examples: + >>> import stim + >>> import numpy as np + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + + >>> # Taking samples. + >>> det_data, obs_data, err_data_not_requested = sampler.sample(shots=4) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data_not_requested is None + True + + >>> # Recording errors. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + + >>> # Bit packing. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True, bit_packed=True) + >>> det_data + array([[6], + [6], + [6], + [6]], dtype=uint8) + >>> obs_data + array([[1], + [1], + [1], + [1]], dtype=uint8) + >>> err_data + array([[2], + [2], + [2], + [2]], dtype=uint8) + + >>> # Recording and replaying errors. + >>> noisy_dem = stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... error(0.25) D1 + ... ''') + >>> noisy_sampler = noisy_dem.compile_sampler() + >>> det_data, obs_data, err_data = noisy_sampler.sample(shots=100, return_errors=True) + >>> replay_det_data, replay_obs_data, _ = noisy_sampler.sample(shots=100, recorded_errors_to_replay=err_data) + >>> np.array_equal(det_data, replay_det_data) + True + >>> np.array_equal(obs_data, replay_obs_data) + True + """ +``` + + +```python +# stim.CompiledDemSampler.sample_write + +# (in class stim.CompiledDemSampler) +def sample_write( + self, + shots: int, + *, + det_out_file: Union[None, str, pathlib.Path], + det_out_format: str = "01", + obs_out_file: Union[None, str, pathlib.Path], + obs_out_format: str = "01", + err_out_file: Union[None, str, pathlib.Path] = None, + err_out_format: str = "01", + replay_err_in_file: Union[None, str, pathlib.Path] = None, + replay_err_in_format: str = "01", +) -> None: + + """Samples the detector error model and writes the results to disk. + + Args: + shots: The number of times to sample from the model. + det_out_file: Where to write detection event data. + If None: detection event data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + det_out_format: The format to write the detection event data in (e.g. "01" or "b8"). + obs_out_file: Where to write observable flip data. + If None: observable flip data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + obs_out_format: The format to write the observable flip data in (e.g. "01" or "b8"). + err_out_file: Where to write errors-that-occurred data. + If None: errors-that-occurred data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + err_out_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + replay_err_in_file: If this is specified, errors are replayed from data instead of generated randomly. + If None: errors are generated randomly according to the probabilities in the detector error model. + If str or pathlib.Path: the file at the given path is opened and errors-to-apply data is read from there. + NOT IMPLEMENTED: io.IOBase + replay_err_in_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + + Returns: + Nothing. Results are written to disk. + + Examples: + >>> import stim + >>> import tempfile + >>> import pathlib + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(0) D1 + ... error(0) D0 + ... error(1) D1 D2 L0 + ... error(0) D0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> with tempfile.TemporaryDirectory() as d: + ... d = pathlib.Path(d) + ... sampler.sample_write( + ... shots=1, + ... det_out_file=d / 'dets.01', + ... det_out_format='01', + ... obs_out_file=d / 'obs.01', + ... obs_out_format='01', + ... err_out_file=d / 'err.hits', + ... err_out_format='hits', + ... ) + ... with open(d / 'dets.01') as f: + ... assert f.read() == "011\n" + ... with open(d / 'obs.01') as f: + ... assert f.read() == "1\n" + ... with open(d / 'err.hits') as f: + ... assert f.read() == "3\n" + """ +``` + ```python # stim.CompiledDetectorSampler @@ -4015,6 +4233,68 @@ def clear( """ ``` + +```python +# stim.DetectorErrorModel.compile_sampler + +# (in class stim.DetectorErrorModel) +def compile_sampler( + self, + *, + seed: object = None, +) -> stim::DemSampler: + """Returns a CompiledDemSampler, which can quickly batch sample from detector error models. + + Args: + seed: PARTIALLY determines simulation results by deterministically seeding the random number generator. + Must be None or an integer in range(2**64). + + Defaults to None. When set to None, a prng seeded by system entropy is used. + + When set to an integer, making the exact same series calls on the exact same machine with the exact + same version of Stim will produce the exact same simulation results. + + CAUTION: simulation results *WILL NOT* be consistent between versions of Stim. This restriction is + present to make it possible to have future optimizations to the random sampling, and is enforced by + introducing intentional differences in the seeding strategy from version to version. + + CAUTION: simulation results *MAY NOT* be consistent across machines that differ in the width of + supported SIMD instructions. For example, using the same seed on a machine that supports AVX + instructions and one that only supports SSE instructions may produce different simulation results. + + CAUTION: simulation results *MAY NOT* be consistent if you vary how many shots are taken. For + example, taking 10 shots and then 90 shots will give different results from taking 100 shots in one + call. + + Returns: + A seeded stim.CompiledDemSampler for the given detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ +``` + ```python # stim.DetectorErrorModel.copy @@ -4037,6 +4317,43 @@ def copy( """ ``` + +```python +# stim.DetectorErrorModel.flattened + +# (in class stim.DetectorErrorModel) +def flattened( + self, +) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model without repeat blocks or detector_shift instructions. + + Returns: + A `stim.DetectorErrorModel` with the same errors in the same order, + but with loops flattened into repeated instructions and with + all coordinate/index shifts inlined. + + Examples: + >>> import stim + >>> stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... REPEAT 5 { + ... error(0.25) D0 D1 + ... shift_detectors 1 + ... } + ... error(0.125) D0 L0 + ... ''').flattened() + stim.DetectorErrorModel(''' + error(0.125) D0 + error(0.25) D0 D1 + error(0.25) D1 D2 + error(0.25) D2 D3 + error(0.25) D3 D4 + error(0.25) D4 D5 + error(0.125) D5 L0 + ''') + """ +``` + ```python # stim.DetectorErrorModel.from_file @@ -4044,9 +4361,14 @@ def copy( # (in class stim.DetectorErrorModel) @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.DetectorErrorModel: - """Args: + + """Reads a detector error model from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md + + Args: file: A file path or open file object to read from. Returns: @@ -4214,6 +4536,49 @@ def num_observables( """ ``` + +```python +# stim.DetectorErrorModel.rounded + +# (in class stim.DetectorErrorModel) +def rounded( + self, + arg0: int, +) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model but with rounded error probabilities. + + Args: + digits: The number of digits to round to. + + Returns: + A `stim.DetectorErrorModel` with the same instructions in the same order, + but with the parens arguments of error instructions rounded to the given + precision. + + Instructions whose error probability was rounded to zero are still + included in the output. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0.019499) D0 + ... error(0.000001) D0 D1 + ... ''') + + >>> dem.rounded(2) + stim.DetectorErrorModel(''' + error(0.02) D0 + error(0) D0 D1 + ''') + + >>> dem.rounded(3) + stim.DetectorErrorModel(''' + error(0.019) D0 + error(0) D0 D1 + ''') + """ +``` + ```python # stim.DetectorErrorModel.shortest_graphlike_error @@ -4309,7 +4674,9 @@ def to_file( file: Union[io.TextIOBase, str, pathlib.Path], ) -> None: - """Writes the stim circuit to a file. + """Writes the detector error model to a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md Args: file: A file path or an open file to write to. @@ -5341,19 +5708,6 @@ def copy( """ ``` - -```python -# stim.PauliString.extended_product - -# (in class stim.PauliString) -def extended_product( - self, - other: stim.PauliString, -) -> Tuple[complex, stim.PauliString]: - """[DEPRECATED] Use multiplication (__mul__ or *) instead. - """ -``` - ```python # stim.PauliString.random diff --git a/doc/stim.pyi b/doc/stim.pyi index 865b84ec1..a0f2c21d1 100644 --- a/doc/stim.pyi +++ b/doc/stim.pyi @@ -781,9 +781,14 @@ class Circuit: """ @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.Circuit: - """Args: + + """Reads a stim circuit from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + + Args: file: A file path or open file object to read from. Returns: @@ -1196,6 +1201,8 @@ class Circuit: """Writes the stim circuit to a file. + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + Args: file: A file path or an open file to write to. @@ -1556,6 +1563,240 @@ class CircuitTargetsInsideInstruction: Includes coordinate data with the targets. """ +class CompiledDemSampler: + """A helper class for efficiently sampler from a detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ + def sample( + self, + shots: int, + *, + bit_packed: bool = False, + return_errors: bool = False, + recorded_errors_to_replay: Optional[np.ndarray] = None, + ) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + + """Samples the detector error model's error mechanisms to produce sample data. + + Args: + shots: The number of times to sample from the model. + bit_packed: Defaults to false. + False: the returned numpy arrays have dtype=np.bool8. + True: the returned numpy arrays have dtype=np.uint8 and pack 8 bits into each byte. + + Setting this to True is equivalent to running np.packbits(data, endian='little', axis=1) + on each output value, but has the performance benefit of the data never being expanded + into an unpacked form. + return_errors: Defaults to False. + False: the first entry of the returned tuple is None. + True: the first entry of the returned tuple is a numpy array recording which errors were sampled. + recorded_errors_to_replay: Defaults to None, meaning sample errors randomly. + If not None, this is expected to be a 2d numpy array specifying which errors to apply (e.g. one + returned from a previous call to the sample method). The array must have + dtype=np.bool8 and shape=(num_shots, num_errors) or + dtype=np.uint8 and shape=(num_shots, math.ceil(num_errors / 8)). + + Returns: + A tuple (detector_data, obs_data, error_data). + + Assuming bit_packed is False and return_errors is True: + If error_data[s, k] is True, then the error with index k fired in the shot with index s. + If detector_data[s, k] is True, then the detector with index k ended up flipped in the shot with index s. + If obs_data[s, k] is True, then the observable with index k ended up flipped in the shot with index s. + + The dtype and shape of the data depends on the arguments: + if bit_packed: + detector_data.shape == (num_shots, num_detectors) + detector_data.dtype == np.bool8 + obs_data.shape == (num_shots, num_observables) + obs_data.dtype == np.bool8 + if return_errors: + error_data.shape = (num_shots, num_errors) + error_data.dtype = np.bool8 + else: + error_data is None + else: + detector_data.shape == (num_shots, math.ceil(num_detectors / 8)) + detector_data.dtype == np.uint8 + obs_data.shape == (num_shots, math.ceil(num_observables / 8)) + obs_data.dtype == np.uint8 + if return_errors: + error_data.shape = (num_shots, math.ceil(num_errors / 8)) + error_data.dtype = np.uint8 + else: + error_data is None + + Note that bit packing is done using little endian order on the last axis + (i.e. like `np.packbits(data, endian='little', axis=1)`). + + Examples: + >>> import stim + >>> import numpy as np + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + + >>> # Taking samples. + >>> det_data, obs_data, err_data_not_requested = sampler.sample(shots=4) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data_not_requested is None + True + + >>> # Recording errors. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + + >>> # Bit packing. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True, bit_packed=True) + >>> det_data + array([[6], + [6], + [6], + [6]], dtype=uint8) + >>> obs_data + array([[1], + [1], + [1], + [1]], dtype=uint8) + >>> err_data + array([[2], + [2], + [2], + [2]], dtype=uint8) + + >>> # Recording and replaying errors. + >>> noisy_dem = stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... error(0.25) D1 + ... ''') + >>> noisy_sampler = noisy_dem.compile_sampler() + >>> det_data, obs_data, err_data = noisy_sampler.sample(shots=100, return_errors=True) + >>> replay_det_data, replay_obs_data, _ = noisy_sampler.sample(shots=100, recorded_errors_to_replay=err_data) + >>> np.array_equal(det_data, replay_det_data) + True + >>> np.array_equal(obs_data, replay_obs_data) + True + """ + def sample_write( + self, + shots: int, + *, + det_out_file: Union[None, str, pathlib.Path], + det_out_format: str = "01", + obs_out_file: Union[None, str, pathlib.Path], + obs_out_format: str = "01", + err_out_file: Union[None, str, pathlib.Path] = None, + err_out_format: str = "01", + replay_err_in_file: Union[None, str, pathlib.Path] = None, + replay_err_in_format: str = "01", + ) -> None: + + """Samples the detector error model and writes the results to disk. + + Args: + shots: The number of times to sample from the model. + det_out_file: Where to write detection event data. + If None: detection event data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + det_out_format: The format to write the detection event data in (e.g. "01" or "b8"). + obs_out_file: Where to write observable flip data. + If None: observable flip data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + obs_out_format: The format to write the observable flip data in (e.g. "01" or "b8"). + err_out_file: Where to write errors-that-occurred data. + If None: errors-that-occurred data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + err_out_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + replay_err_in_file: If this is specified, errors are replayed from data instead of generated randomly. + If None: errors are generated randomly according to the probabilities in the detector error model. + If str or pathlib.Path: the file at the given path is opened and errors-to-apply data is read from there. + NOT IMPLEMENTED: io.IOBase + replay_err_in_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + + Returns: + Nothing. Results are written to disk. + + Examples: + >>> import stim + >>> import tempfile + >>> import pathlib + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(0) D1 + ... error(0) D0 + ... error(1) D1 D2 L0 + ... error(0) D0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> with tempfile.TemporaryDirectory() as d: + ... d = pathlib.Path(d) + ... sampler.sample_write( + ... shots=1, + ... det_out_file=d / 'dets.01', + ... det_out_format='01', + ... obs_out_file=d / 'obs.01', + ... obs_out_format='01', + ... err_out_file=d / 'err.hits', + ... err_out_format='hits', + ... ) + ... with open(d / 'dets.01') as f: + ... assert f.read() == "011\n" + ... with open(d / 'obs.01') as f: + ... assert f.read() == "1\n" + ... with open(d / 'err.hits') as f: + ... assert f.read() == "3\n" + """ class CompiledDetectorSampler: """An analyzed stabilizer circuit whose detection events can be sampled quickly. """ @@ -2725,6 +2966,61 @@ class DetectorErrorModel: >>> model stim.DetectorErrorModel() """ + def compile_sampler( + self, + *, + seed: object = None, + ) -> stim::DemSampler: + """Returns a CompiledDemSampler, which can quickly batch sample from detector error models. + + Args: + seed: PARTIALLY determines simulation results by deterministically seeding the random number generator. + Must be None or an integer in range(2**64). + + Defaults to None. When set to None, a prng seeded by system entropy is used. + + When set to an integer, making the exact same series calls on the exact same machine with the exact + same version of Stim will produce the exact same simulation results. + + CAUTION: simulation results *WILL NOT* be consistent between versions of Stim. This restriction is + present to make it possible to have future optimizations to the random sampling, and is enforced by + introducing intentional differences in the seeding strategy from version to version. + + CAUTION: simulation results *MAY NOT* be consistent across machines that differ in the width of + supported SIMD instructions. For example, using the same seed on a machine that supports AVX + instructions and one that only supports SSE instructions may produce different simulation results. + + CAUTION: simulation results *MAY NOT* be consistent if you vary how many shots are taken. For + example, taking 10 shots and then 90 shots will give different results from taking 100 shots in one + call. + + Returns: + A seeded stim.CompiledDemSampler for the given detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ def copy( self, ) -> stim.DetectorErrorModel: @@ -2740,11 +3036,46 @@ class DetectorErrorModel: >>> c2 == c1 True """ + def flattened( + self, + ) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model without repeat blocks or detector_shift instructions. + + Returns: + A `stim.DetectorErrorModel` with the same errors in the same order, + but with loops flattened into repeated instructions and with + all coordinate/index shifts inlined. + + Examples: + >>> import stim + >>> stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... REPEAT 5 { + ... error(0.25) D0 D1 + ... shift_detectors 1 + ... } + ... error(0.125) D0 L0 + ... ''').flattened() + stim.DetectorErrorModel(''' + error(0.125) D0 + error(0.25) D0 D1 + error(0.25) D1 D2 + error(0.25) D2 D3 + error(0.25) D3 D4 + error(0.25) D4 D5 + error(0.125) D5 L0 + ''') + """ @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.DetectorErrorModel: - """Args: + + """Reads a detector error model from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md + + Args: file: A file path or open file object to read from. Returns: @@ -2882,6 +3213,42 @@ class DetectorErrorModel: ... ''').num_observables 400 """ + def rounded( + self, + arg0: int, + ) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model but with rounded error probabilities. + + Args: + digits: The number of digits to round to. + + Returns: + A `stim.DetectorErrorModel` with the same instructions in the same order, + but with the parens arguments of error instructions rounded to the given + precision. + + Instructions whose error probability was rounded to zero are still + included in the output. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0.019499) D0 + ... error(0.000001) D0 D1 + ... ''') + + >>> dem.rounded(2) + stim.DetectorErrorModel(''' + error(0.02) D0 + error(0) D0 D1 + ''') + + >>> dem.rounded(3) + stim.DetectorErrorModel(''' + error(0.019) D0 + error(0) D0 D1 + ''') + """ def shortest_graphlike_error( self, ignore_ungraphlike_errors: bool = False, @@ -2965,7 +3332,9 @@ class DetectorErrorModel: file: Union[io.TextIOBase, str, pathlib.Path], ) -> None: - """Writes the stim circuit to a file. + """Writes the detector error model to a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md Args: file: A file path or an open file to write to. diff --git a/doc/usage_command_line.md b/doc/usage_command_line.md index 9b0b04d2b..1812e6dd0 100644 --- a/doc/usage_command_line.md +++ b/doc/usage_command_line.md @@ -48,6 +48,18 @@ - [--seed](#--seed) - [--shots](#--shots) - [--skip_reference_sample](#--skip_reference_sample) +- **(mode)** [stim sample_dem](#sample_dem) + - [--err_out](#--err_out) + - [--err_out_format](#--err_out_format) + - [--in](#--in) + - [--obs_out](#--obs_out) + - [--obs_out_format](#--obs_out_format) + - [--out](#--out) + - [--out_format](#--out_format) + - [--replay_err_in](#--replay_err_in) + - [--replay_err_in_format](#--replay_err_in_format) + - [--seed](#--seed) + - [--shots](#--shots) ## Modes @@ -437,6 +449,56 @@ Flags used with this mode: - [--shots](#--shots) - [--skip_reference_sample](#--skip_reference_sample) + +### stim sample_dem + +*Samples detection events and observable flips from a detector error model.* + +stdin (or --in): The detector error model to sample from, specified using the [detector error model file format](https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md). + +stdout (or --out): The detection event data is written here. + +- Example: + + ```bash + echo "error(0) D0" > example.dem + echo "error(0.5) D1 L0" >> example.dem + echo "error(1) D2 D3" >> example.dem + stim sample_dem \ + --shots 5 \ + --in example.dem \ + --out dets.01 \ + --out_format 01 \ + --obs_out obs_flips.01 \ + --obs_out_format 01 \ + --seed 0 + cat dets.01 + # 0111 + # 0011 + # 0011 + # 0111 + # 0111 + cat obs_flips.01 + # 1 + # 0 + # 0 + # 1 + # 1 + ``` + +Flags used with this mode: +- [--err_out](#--err_out) +- [--err_out_format](#--err_out_format) +- [--in](#--in) +- [--obs_out](#--obs_out) +- [--obs_out_format](#--obs_out_format) +- [--out](#--out) +- [--out_format](#--out_format) +- [--replay_err_in](#--replay_err_in) +- [--replay_err_in_format](#--replay_err_in_format) +- [--seed](#--seed) +- [--shots](#--shots) + ## Flags - **`--after_clifford_depolarization`** @@ -600,6 +662,16 @@ Flags used with this mode: (e.g. must be larger than 2 or must be odd or etc). +- **`--err_out`** + Specifies a file to write a record of which errors occurred. + + This data can then be analyzed, modified, and later given to for example a --replay_err_in argument. + + +- **`--err_out_format`** + The format to use when writing error data (e.g. b8 or 01). + + - **`--fold_loops`** Allows the output error model to contain `repeat` blocks. @@ -671,6 +743,17 @@ Flags used with this mode: See `stim help formats` for a list of supported formats. +- **`--replay_err_in`** + Specifies a file to read error data to replay from. + + When replaying error information, errors are no longer sampled randomly but instead driven by the file data. + For example, this file data could come from a previous run that wrote error data using --err_out. + + +- **`--replay_err_in_format`** + The format to use when reading error data to replay. (e.g. b8 or 01). + + - **`--rounds`** The number of times the circuit's measurement qubits are measured. diff --git a/glue/python/generate_api_reference.py b/glue/python/generate_api_reference.py index 4615b2f29..55ab26cad 100644 --- a/glue/python/generate_api_reference.py +++ b/glue/python/generate_api_reference.py @@ -74,7 +74,11 @@ def main(): else: version = "v" + version is_dev = False - objects = list(generate_documentation(obj=stim, full_name="stim", level=0)) + objects = [ + obj + for obj in generate_documentation(obj=stim, full_name="stim", level=0) + if all('[DEPRECATED]' not in line for line in obj.lines) + ] print(f"# Stim {version} API Reference") print() diff --git a/glue/python/src/stim/__init__.pyi b/glue/python/src/stim/__init__.pyi index 865b84ec1..a0f2c21d1 100644 --- a/glue/python/src/stim/__init__.pyi +++ b/glue/python/src/stim/__init__.pyi @@ -781,9 +781,14 @@ class Circuit: """ @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.Circuit: - """Args: + + """Reads a stim circuit from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + + Args: file: A file path or open file object to read from. Returns: @@ -1196,6 +1201,8 @@ class Circuit: """Writes the stim circuit to a file. + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + Args: file: A file path or an open file to write to. @@ -1556,6 +1563,240 @@ class CircuitTargetsInsideInstruction: Includes coordinate data with the targets. """ +class CompiledDemSampler: + """A helper class for efficiently sampler from a detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ + def sample( + self, + shots: int, + *, + bit_packed: bool = False, + return_errors: bool = False, + recorded_errors_to_replay: Optional[np.ndarray] = None, + ) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + + """Samples the detector error model's error mechanisms to produce sample data. + + Args: + shots: The number of times to sample from the model. + bit_packed: Defaults to false. + False: the returned numpy arrays have dtype=np.bool8. + True: the returned numpy arrays have dtype=np.uint8 and pack 8 bits into each byte. + + Setting this to True is equivalent to running np.packbits(data, endian='little', axis=1) + on each output value, but has the performance benefit of the data never being expanded + into an unpacked form. + return_errors: Defaults to False. + False: the first entry of the returned tuple is None. + True: the first entry of the returned tuple is a numpy array recording which errors were sampled. + recorded_errors_to_replay: Defaults to None, meaning sample errors randomly. + If not None, this is expected to be a 2d numpy array specifying which errors to apply (e.g. one + returned from a previous call to the sample method). The array must have + dtype=np.bool8 and shape=(num_shots, num_errors) or + dtype=np.uint8 and shape=(num_shots, math.ceil(num_errors / 8)). + + Returns: + A tuple (detector_data, obs_data, error_data). + + Assuming bit_packed is False and return_errors is True: + If error_data[s, k] is True, then the error with index k fired in the shot with index s. + If detector_data[s, k] is True, then the detector with index k ended up flipped in the shot with index s. + If obs_data[s, k] is True, then the observable with index k ended up flipped in the shot with index s. + + The dtype and shape of the data depends on the arguments: + if bit_packed: + detector_data.shape == (num_shots, num_detectors) + detector_data.dtype == np.bool8 + obs_data.shape == (num_shots, num_observables) + obs_data.dtype == np.bool8 + if return_errors: + error_data.shape = (num_shots, num_errors) + error_data.dtype = np.bool8 + else: + error_data is None + else: + detector_data.shape == (num_shots, math.ceil(num_detectors / 8)) + detector_data.dtype == np.uint8 + obs_data.shape == (num_shots, math.ceil(num_observables / 8)) + obs_data.dtype == np.uint8 + if return_errors: + error_data.shape = (num_shots, math.ceil(num_errors / 8)) + error_data.dtype = np.uint8 + else: + error_data is None + + Note that bit packing is done using little endian order on the last axis + (i.e. like `np.packbits(data, endian='little', axis=1)`). + + Examples: + >>> import stim + >>> import numpy as np + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + + >>> # Taking samples. + >>> det_data, obs_data, err_data_not_requested = sampler.sample(shots=4) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data_not_requested is None + True + + >>> # Recording errors. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + + >>> # Bit packing. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True, bit_packed=True) + >>> det_data + array([[6], + [6], + [6], + [6]], dtype=uint8) + >>> obs_data + array([[1], + [1], + [1], + [1]], dtype=uint8) + >>> err_data + array([[2], + [2], + [2], + [2]], dtype=uint8) + + >>> # Recording and replaying errors. + >>> noisy_dem = stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... error(0.25) D1 + ... ''') + >>> noisy_sampler = noisy_dem.compile_sampler() + >>> det_data, obs_data, err_data = noisy_sampler.sample(shots=100, return_errors=True) + >>> replay_det_data, replay_obs_data, _ = noisy_sampler.sample(shots=100, recorded_errors_to_replay=err_data) + >>> np.array_equal(det_data, replay_det_data) + True + >>> np.array_equal(obs_data, replay_obs_data) + True + """ + def sample_write( + self, + shots: int, + *, + det_out_file: Union[None, str, pathlib.Path], + det_out_format: str = "01", + obs_out_file: Union[None, str, pathlib.Path], + obs_out_format: str = "01", + err_out_file: Union[None, str, pathlib.Path] = None, + err_out_format: str = "01", + replay_err_in_file: Union[None, str, pathlib.Path] = None, + replay_err_in_format: str = "01", + ) -> None: + + """Samples the detector error model and writes the results to disk. + + Args: + shots: The number of times to sample from the model. + det_out_file: Where to write detection event data. + If None: detection event data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + det_out_format: The format to write the detection event data in (e.g. "01" or "b8"). + obs_out_file: Where to write observable flip data. + If None: observable flip data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + obs_out_format: The format to write the observable flip data in (e.g. "01" or "b8"). + err_out_file: Where to write errors-that-occurred data. + If None: errors-that-occurred data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + err_out_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + replay_err_in_file: If this is specified, errors are replayed from data instead of generated randomly. + If None: errors are generated randomly according to the probabilities in the detector error model. + If str or pathlib.Path: the file at the given path is opened and errors-to-apply data is read from there. + NOT IMPLEMENTED: io.IOBase + replay_err_in_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + + Returns: + Nothing. Results are written to disk. + + Examples: + >>> import stim + >>> import tempfile + >>> import pathlib + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(0) D1 + ... error(0) D0 + ... error(1) D1 D2 L0 + ... error(0) D0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> with tempfile.TemporaryDirectory() as d: + ... d = pathlib.Path(d) + ... sampler.sample_write( + ... shots=1, + ... det_out_file=d / 'dets.01', + ... det_out_format='01', + ... obs_out_file=d / 'obs.01', + ... obs_out_format='01', + ... err_out_file=d / 'err.hits', + ... err_out_format='hits', + ... ) + ... with open(d / 'dets.01') as f: + ... assert f.read() == "011\n" + ... with open(d / 'obs.01') as f: + ... assert f.read() == "1\n" + ... with open(d / 'err.hits') as f: + ... assert f.read() == "3\n" + """ class CompiledDetectorSampler: """An analyzed stabilizer circuit whose detection events can be sampled quickly. """ @@ -2725,6 +2966,61 @@ class DetectorErrorModel: >>> model stim.DetectorErrorModel() """ + def compile_sampler( + self, + *, + seed: object = None, + ) -> stim::DemSampler: + """Returns a CompiledDemSampler, which can quickly batch sample from detector error models. + + Args: + seed: PARTIALLY determines simulation results by deterministically seeding the random number generator. + Must be None or an integer in range(2**64). + + Defaults to None. When set to None, a prng seeded by system entropy is used. + + When set to an integer, making the exact same series calls on the exact same machine with the exact + same version of Stim will produce the exact same simulation results. + + CAUTION: simulation results *WILL NOT* be consistent between versions of Stim. This restriction is + present to make it possible to have future optimizations to the random sampling, and is enforced by + introducing intentional differences in the seeding strategy from version to version. + + CAUTION: simulation results *MAY NOT* be consistent across machines that differ in the width of + supported SIMD instructions. For example, using the same seed on a machine that supports AVX + instructions and one that only supports SSE instructions may produce different simulation results. + + CAUTION: simulation results *MAY NOT* be consistent if you vary how many shots are taken. For + example, taking 10 shots and then 90 shots will give different results from taking 100 shots in one + call. + + Returns: + A seeded stim.CompiledDemSampler for the given detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + """ def copy( self, ) -> stim.DetectorErrorModel: @@ -2740,11 +3036,46 @@ class DetectorErrorModel: >>> c2 == c1 True """ + def flattened( + self, + ) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model without repeat blocks or detector_shift instructions. + + Returns: + A `stim.DetectorErrorModel` with the same errors in the same order, + but with loops flattened into repeated instructions and with + all coordinate/index shifts inlined. + + Examples: + >>> import stim + >>> stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... REPEAT 5 { + ... error(0.25) D0 D1 + ... shift_detectors 1 + ... } + ... error(0.125) D0 L0 + ... ''').flattened() + stim.DetectorErrorModel(''' + error(0.125) D0 + error(0.25) D0 D1 + error(0.25) D1 D2 + error(0.25) D2 D3 + error(0.25) D3 D4 + error(0.25) D4 D5 + error(0.125) D5 L0 + ''') + """ @staticmethod def from_file( - file: object, + file: Union[io.TextIOBase, str, pathlib.Path], ) -> stim.DetectorErrorModel: - """Args: + + """Reads a detector error model from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md + + Args: file: A file path or open file object to read from. Returns: @@ -2882,6 +3213,42 @@ class DetectorErrorModel: ... ''').num_observables 400 """ + def rounded( + self, + arg0: int, + ) -> stim.DetectorErrorModel: + """Creates an equivalent detector error model but with rounded error probabilities. + + Args: + digits: The number of digits to round to. + + Returns: + A `stim.DetectorErrorModel` with the same instructions in the same order, + but with the parens arguments of error instructions rounded to the given + precision. + + Instructions whose error probability was rounded to zero are still + included in the output. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0.019499) D0 + ... error(0.000001) D0 D1 + ... ''') + + >>> dem.rounded(2) + stim.DetectorErrorModel(''' + error(0.02) D0 + error(0) D0 D1 + ''') + + >>> dem.rounded(3) + stim.DetectorErrorModel(''' + error(0.019) D0 + error(0) D0 D1 + ''') + """ def shortest_graphlike_error( self, ignore_ungraphlike_errors: bool = False, @@ -2965,7 +3332,9 @@ class DetectorErrorModel: file: Union[io.TextIOBase, str, pathlib.Path], ) -> None: - """Writes the stim circuit to a file. + """Writes the detector error model to a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md Args: file: A file path or an open file to write to. diff --git a/src/stim/circuit/circuit.pybind.cc b/src/stim/circuit/circuit.pybind.cc index 2de8dd1ed..777832c2d 100644 --- a/src/stim/circuit/circuit.pybind.cc +++ b/src/stim/circuit/circuit.pybind.cc @@ -947,6 +947,11 @@ pybind11::class_ pybind_circuit(pybind11::module &m) { }, pybind11::arg("file"), clean_doc_string(u8R"DOC( + @signature def from_file(file: Union[io.TextIOBase, str, pathlib.Path]) -> stim.Circuit: + Reads a stim circuit from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + Args: file: A file path or open file object to read from. @@ -1020,6 +1025,8 @@ pybind11::class_ pybind_circuit(pybind11::module &m) { @signature def to_file(self, file: Union[io.TextIOBase, str, pathlib.Path]) -> None: Writes the stim circuit to a file. + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_stim_circuit.md + Args: file: A file path or an open file to write to. diff --git a/src/stim/dem/detector_error_model.cc b/src/stim/dem/detector_error_model.cc index 1bc3171f5..c3b27e363 100644 --- a/src/stim/dem/detector_error_model.cc +++ b/src/stim/dem/detector_error_model.cc @@ -580,9 +580,9 @@ void DetectorErrorModel::clear() { blocks.clear(); } -DetectorErrorModel DetectorErrorModel::rounded(uint8_t sig_figs) const { +DetectorErrorModel DetectorErrorModel::rounded(uint8_t digits) const { double scale = 1; - for (size_t k = 0; k < sig_figs; k++) { + for (size_t k = 0; k < digits; k++) { scale *= 10; } @@ -591,13 +591,15 @@ DetectorErrorModel DetectorErrorModel::rounded(uint8_t sig_figs) const { if (e.type == DEM_REPEAT_BLOCK) { auto reps = e.target_data[0].data; auto &block = blocks[e.target_data[1].data]; - result.append_repeat_block(reps, block.rounded(sig_figs)); - } else { + result.append_repeat_block(reps, block.rounded(digits)); + } else if (e.type == DEM_ERROR) { std::vector rounded_args; for (auto a : e.arg_data) { rounded_args.push_back(round(a * scale) / scale); } - result.append_dem_instruction({rounded_args, e.target_data, e.type}); + result.append_dem_instruction({rounded_args, e.target_data, DEM_ERROR}); + } else { + result.append_dem_instruction(e); } } return result; @@ -615,6 +617,68 @@ uint64_t DetectorErrorModel::total_detector_shift() const { return result; } +void flattened_helper( + const DetectorErrorModel &body, + std::vector &cur_coordinate_shift, + uint64_t &cur_detector_shift, + DetectorErrorModel &out) { + for (const auto &op : body.instructions) { + if (op.type == DEM_SHIFT_DETECTORS) { + while (cur_coordinate_shift.size() < op.arg_data.size()) { + cur_coordinate_shift.push_back(0); + } + for (size_t k = 0; k < op.arg_data.size(); k++) { + cur_coordinate_shift[k] += op.arg_data[k]; + } + if (!op.target_data.empty()) { + cur_detector_shift += op.target_data[0].data; + } + } else if (op.type == DEM_REPEAT_BLOCK) { + const auto &loop_body = body.blocks[op.target_data[1].data]; + auto reps = op.target_data[0].data; + for (uint64_t k = 0; k < reps; k++) { + flattened_helper(loop_body, cur_coordinate_shift, cur_detector_shift, out); + } + } else if (op.type == DEM_LOGICAL_OBSERVABLE) { + out.append_dem_instruction(DemInstruction{{}, op.target_data, DEM_LOGICAL_OBSERVABLE}); + } else if (op.type == DEM_DETECTOR) { + while (cur_coordinate_shift.size() < op.arg_data.size()) { + cur_coordinate_shift.push_back(0); + } + + std::vector shifted_coords; + for (size_t k = 0; k < op.arg_data.size(); k++) { + shifted_coords.push_back(op.arg_data[k] + cur_coordinate_shift[k]); + } + std::vector shifted_detectors; + for (DemTarget t : op.target_data) { + t.shift_if_detector_id(cur_detector_shift); + shifted_detectors.push_back(t); + } + + out.append_dem_instruction(DemInstruction{shifted_coords, shifted_detectors, DEM_DETECTOR}); + } else if (op.type == DEM_ERROR) { + std::vector shifted_detectors; + for (DemTarget t : op.target_data) { + t.shift_if_detector_id(cur_detector_shift); + shifted_detectors.push_back(t); + } + + out.append_dem_instruction(DemInstruction{op.arg_data, shifted_detectors, DEM_ERROR}); + } else { + throw std::invalid_argument("Unrecognized instruction type: " + op.str()); + } + } +} + +DetectorErrorModel DetectorErrorModel::flattened() const { + DetectorErrorModel result; + std::vector shift; + uint64_t det_shift = 0; + flattened_helper(*this, shift, det_shift, result); + return result; +} + uint64_t DetectorErrorModel::count_detectors() const { uint64_t offset = 1; uint64_t max_num = 0; diff --git a/src/stim/dem/detector_error_model.h b/src/stim/dem/detector_error_model.h index 1ee0332c6..175423af0 100644 --- a/src/stim/dem/detector_error_model.h +++ b/src/stim/dem/detector_error_model.h @@ -184,7 +184,11 @@ struct DetectorErrorModel { /// Gets a python-style slice of the error model's instructions. DetectorErrorModel py_get_slice(int64_t start, int64_t step, int64_t slice_length) const; - DetectorErrorModel rounded(uint8_t sig_figs) const; + /// Rounds error probabilities to a given number of digits. + DetectorErrorModel rounded(uint8_t digits) const; + + /// Returns an equivalent detector error model with no repeat blocks or detector_shift instructions. + DetectorErrorModel flattened() const; }; void print_detector_error_model(std::ostream &out, const DetectorErrorModel &v, size_t indent); diff --git a/src/stim/dem/detector_error_model.pybind.cc b/src/stim/dem/detector_error_model.pybind.cc index 91c72cd93..e321ba2e4 100644 --- a/src/stim/dem/detector_error_model.pybind.cc +++ b/src/stim/dem/detector_error_model.pybind.cc @@ -23,6 +23,7 @@ #include "stim/io/raii_file.h" #include "stim/py/base.pybind.h" #include "stim/search/search.h" +#include "stim/simulators/dem_sampler.h" using namespace stim; using namespace stim_pybind; @@ -821,6 +822,11 @@ void pybind_detector_error_model(pybind11::module &m) { }, pybind11::arg("file"), clean_doc_string(u8R"DOC( + @signature def from_file(file: Union[io.TextIOBase, str, pathlib.Path]) -> stim.DetectorErrorModel: + Reads a detector error model from a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md + Args: file: A file path or open file object to read from. @@ -892,7 +898,9 @@ void pybind_detector_error_model(pybind11::module &m) { pybind11::arg("file"), clean_doc_string(u8R"DOC( @signature def to_file(self, file: Union[io.TextIOBase, str, pathlib.Path]) -> None: - Writes the stim circuit to a file. + Writes the detector error model to a file. + + The file format is defined at https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md Args: file: A file path or an open file to write to. @@ -920,4 +928,135 @@ void pybind_detector_error_model(pybind11::module &m) { 'error(0.25) D2 D3\n' )DOC") .data()); + + c.def( + "compile_sampler", + [](const DetectorErrorModel &self, const pybind11::object &seed) { + return DemSampler(self, *make_py_seeded_rng(seed), 1024); + }, + pybind11::kw_only(), + pybind11::arg("seed") = pybind11::none(), + clean_doc_string(u8R"DOC( + Returns a CompiledDemSampler, which can quickly batch sample from detector error models. + + Args: + seed: PARTIALLY determines simulation results by deterministically seeding the random number generator. + Must be None or an integer in range(2**64). + + Defaults to None. When set to None, a prng seeded by system entropy is used. + + When set to an integer, making the exact same series calls on the exact same machine with the exact + same version of Stim will produce the exact same simulation results. + + CAUTION: simulation results *WILL NOT* be consistent between versions of Stim. This restriction is + present to make it possible to have future optimizations to the random sampling, and is enforced by + introducing intentional differences in the seeding strategy from version to version. + + CAUTION: simulation results *MAY NOT* be consistent across machines that differ in the width of + supported SIMD instructions. For example, using the same seed on a machine that supports AVX + instructions and one that only supports SSE instructions may produce different simulation results. + + CAUTION: simulation results *MAY NOT* be consistent if you vary how many shots are taken. For + example, taking 10 shots and then 90 shots will give different results from taking 100 shots in one + call. + + Returns: + A seeded stim.CompiledDemSampler for the given detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + )DOC") + .data()); + + c.def( + "flattened", + &DetectorErrorModel::flattened, + clean_doc_string(u8R"DOC( + Creates an equivalent detector error model without repeat blocks or detector_shift instructions. + + Returns: + A `stim.DetectorErrorModel` with the same errors in the same order, + but with loops flattened into repeated instructions and with + all coordinate/index shifts inlined. + + Examples: + >>> import stim + >>> stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... REPEAT 5 { + ... error(0.25) D0 D1 + ... shift_detectors 1 + ... } + ... error(0.125) D0 L0 + ... ''').flattened() + stim.DetectorErrorModel(''' + error(0.125) D0 + error(0.25) D0 D1 + error(0.25) D1 D2 + error(0.25) D2 D3 + error(0.25) D3 D4 + error(0.25) D4 D5 + error(0.125) D5 L0 + ''') + )DOC") + .data()); + + c.def( + "rounded", + &DetectorErrorModel::rounded, + clean_doc_string(u8R"DOC( + Creates an equivalent detector error model but with rounded error probabilities. + + Args: + digits: The number of digits to round to. + + Returns: + A `stim.DetectorErrorModel` with the same instructions in the same order, + but with the parens arguments of error instructions rounded to the given + precision. + + Instructions whose error probability was rounded to zero are still + included in the output. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0.019499) D0 + ... error(0.000001) D0 D1 + ... ''') + + >>> dem.rounded(2) + stim.DetectorErrorModel(''' + error(0.02) D0 + error(0) D0 D1 + ''') + + >>> dem.rounded(3) + stim.DetectorErrorModel(''' + error(0.019) D0 + error(0) D0 D1 + ''') + )DOC") + .data()); } diff --git a/src/stim/dem/detector_error_model.test.cc b/src/stim/dem/detector_error_model.test.cc index 155461fd3..553da1258 100644 --- a/src/stim/dem/detector_error_model.test.cc +++ b/src/stim/dem/detector_error_model.test.cc @@ -751,8 +751,8 @@ TEST(detector_error_model, rounded) { repeat 2 { error(0) D1 D2 L3 } - detector(0,0) D0 - shift_detectors(5,0) 3 + detector(0.0200000334,0.12345) D0 + shift_detectors(5.0300004,0.12345) 3 )DEM")); ASSERT_EQ(dem.rounded(1), DetectorErrorModel(R"DEM( @@ -760,8 +760,8 @@ TEST(detector_error_model, rounded) { repeat 2 { error(0.1) D1 D2 L3 } - detector(0,0.1) D0 - shift_detectors(5,0.1) 3 + detector(0.0200000334,0.12345) D0 + shift_detectors(5.0300004,0.12345) 3 )DEM")); ASSERT_EQ(dem.rounded(2), DetectorErrorModel(R"DEM( @@ -769,8 +769,8 @@ TEST(detector_error_model, rounded) { repeat 2 { error(0.12) D1 D2 L3 } - detector(0.02,0.12) D0 - shift_detectors(5.03,0.12) 3 + detector(0.0200000334,0.12345) D0 + shift_detectors(5.0300004,0.12345) 3 )DEM")); ASSERT_EQ(dem.rounded(3), DetectorErrorModel(R"DEM( @@ -778,8 +778,8 @@ TEST(detector_error_model, rounded) { repeat 2 { error(0.123) D1 D2 L3 } - detector(0.020,0.123) D0 - shift_detectors(5.030,0.123) 3 + detector(0.0200000334,0.12345) D0 + shift_detectors(5.0300004,0.12345) 3 )DEM")); } @@ -802,3 +802,74 @@ TEST(detector_error_model, surface_code_coords_dont_infinite_loop) { ASSERT_EQ(coords1.size(), n); ASSERT_EQ(n, 168); } + +TEST(detector_error_model, flattened) { + ASSERT_EQ(DetectorErrorModel().flattened(), DetectorErrorModel()); + + ASSERT_EQ( + DetectorErrorModel(R"DEM( + error(0.125) D0 D1 L0 + )DEM") + .flattened(), + DetectorErrorModel(R"DEM( + error(0.125) D0 D1 L0 + )DEM")); + + ASSERT_EQ( + DetectorErrorModel(R"DEM( + error(0.125) D0 D1 L0 + shift_detectors 5 + )DEM") + .flattened(), + DetectorErrorModel(R"DEM( + error(0.125) D0 D1 L0 + )DEM")); + + ASSERT_EQ( + DetectorErrorModel(R"DEM( + shift_detectors 5 + error(0.125) D0 D1 L0 + )DEM") + .flattened(), + DetectorErrorModel(R"DEM( + error(0.125) D5 D6 L0 + )DEM")); + + ASSERT_EQ( + DetectorErrorModel(R"DEM( + detector(10, 20) D0 + detector(10, 20, 30, 40) D1 + logical_observable L0 + shift_detectors(1, 2, 3) 5 + detector(10, 20) D0 + detector(10, 20, 30, 40) D1 + logical_observable L1 + )DEM") + .flattened(), + DetectorErrorModel(R"DEM( + detector(10, 20) D0 + detector(10, 20, 30, 40) D1 + logical_observable L0 + detector(11, 22) D5 + detector(11, 22, 33, 40) D6 + logical_observable L1 + )DEM")); + + ASSERT_EQ( + DetectorErrorModel(R"DEM( + repeat 5 { + error(0.125) D0 + shift_detectors(3) 2 + } + detector(10, 20, 30, 40) D0 + )DEM") + .flattened(), + DetectorErrorModel(R"DEM( + error(0.125) D0 + error(0.125) D2 + error(0.125) D4 + error(0.125) D6 + error(0.125) D8 + detector(25, 20, 30, 40) D10 + )DEM")); +} diff --git a/src/stim/dem/detector_error_model_pybind_test.py b/src/stim/dem/detector_error_model_pybind_test.py index c575a7372..41e1afd09 100644 --- a/src/stim/dem/detector_error_model_pybind_test.py +++ b/src/stim/dem/detector_error_model_pybind_test.py @@ -399,3 +399,53 @@ def test_dem_to_file(): c.to_file(object()) with pytest.raises(ValueError, match="how to write"): c.to_file(123) + + +def test_flattened(): + dem = stim.DetectorErrorModel(""" + shift_detectors 5 + repeat 2 { + error(0.125) D0 D1 + } + """) + assert dem.flattened() == stim.DetectorErrorModel(""" + error(0.125) D5 D6 + error(0.125) D5 D6 + """) + + +def test_rounded(): + dem = stim.DetectorErrorModel(""" + error(0.1248) D0 D1 + """) + assert dem.rounded(1) == stim.DetectorErrorModel(""" + error(0.1) D0 D1 + """) + assert dem.rounded(2) == stim.DetectorErrorModel(""" + error(0.12) D0 D1 + """) + assert dem.rounded(3) == stim.DetectorErrorModel(""" + error(0.125) D0 D1 + """) + assert dem.rounded(4) == stim.DetectorErrorModel(""" + error(0.1248) D0 D1 + """) + assert dem.rounded(5) == stim.DetectorErrorModel(""" + error(0.1248) D0 D1 + """) + + dem = stim.DetectorErrorModel(""" + error(0.01248) D0 D1 + """) + assert dem.rounded(1) == stim.DetectorErrorModel(""" + error(0) D0 D1 + """) + assert dem.rounded(2) == stim.DetectorErrorModel(""" + error(0.01) D0 D1 + """) + assert dem.rounded(3) == stim.DetectorErrorModel(""" + error(0.012) D0 D1 + """) + assert dem.rounded(4) == stim.DetectorErrorModel(""" + error(0.0125) D0 D1 + """) diff --git a/src/stim/help.cc b/src/stim/help.cc index e02171f34..eb1475113 100644 --- a/src/stim/help.cc +++ b/src/stim/help.cc @@ -129,6 +129,56 @@ stdout: The sample data. {"--out_format", "--seed", "--in", "--out", "--skip_reference_sample", "--shots"}, }; + modes["sample_dem"] = CommandLineSingleModeData{ + "Samples detection events and observable flips from a detector error model.", + R"PARAGRAPH( +stdin (or --in): The detector error model to sample from, specified using the [detector error model file format](https://github.com/quantumlib/Stim/blob/main/doc/file_format_dem_detector_error_model.md). + +stdout (or --out): The detection event data is written here. + +- Example: + + ```bash + echo "error(0) D0" > example.dem + echo "error(0.5) D1 L0" >> example.dem + echo "error(1) D2 D3" >> example.dem + stim sample_dem \ + --shots 5 \ + --in example.dem \ + --out dets.01 \ + --out_format 01 \ + --obs_out obs_flips.01 \ + --obs_out_format 01 \ + --seed 0 + cat dets.01 + # 0111 + # 0011 + # 0011 + # 0111 + # 0111 + cat obs_flips.01 + # 1 + # 0 + # 0 + # 1 + # 1 + ``` +)PARAGRAPH", + { + "--in", + "--out", + "--out_format", + "--obs_out", + "--obs_out_format", + "--seed", + "--shots", + "--err_out", + "--err_out_format", + "--replay_err_in", + "--replay_err_in_format", + }, + }; + modes["explain_errors"] = CommandLineSingleModeData{ "Describes how detector error model errors correspond to circuit errors.", R"PARAGRAPH( @@ -367,6 +417,23 @@ only in shots where the corresponding sweep data has the bit at index 5 set to T flags["--sweep_format"] = R"PARAGRAPH(Specifies the format sweep data is stored in (e.g. b8 or 01). )PARAGRAPH"; + flags["--err_out"] = R"PARAGRAPH(Specifies a file to write a record of which errors occurred. + +This data can then be analyzed, modified, and later given to for example a --replay_err_in argument. +)PARAGRAPH"; + + flags["--err_out_format"] = R"PARAGRAPH(The format to use when writing error data (e.g. b8 or 01). +)PARAGRAPH"; + + flags["--replay_err_in"] = R"PARAGRAPH(Specifies a file to read error data to replay from. + +When replaying error information, errors are no longer sampled randomly but instead driven by the file data. +For example, this file data could come from a previous run that wrote error data using --err_out. +)PARAGRAPH"; + + flags["--replay_err_in_format"] = R"PARAGRAPH(The format to use when reading error data to replay. (e.g. b8 or 01). +)PARAGRAPH"; + flags["--obs_out"] = R"PARAGRAPH(Specifies a file to write observable flip data to. When sampling detection event data, this is an alternative to --append_observables which has the benefit diff --git a/src/stim/io/measure_record_reader.cc b/src/stim/io/measure_record_reader.cc index aa4b3822b..325e90544 100644 --- a/src/stim/io/measure_record_reader.cc +++ b/src/stim/io/measure_record_reader.cc @@ -112,6 +112,15 @@ void MeasureRecordReader::move_obs_in_shots_to_mask_assuming_sorted(SparseShot & } } +size_t MeasureRecordReader::read_into_table_with_major_shot_index(simd_bit_table &out_table) { + size_t read_shots = 0; + size_t max_shots = out_table.num_major_bits_padded(); + while (read_shots < max_shots && start_and_read_entire_record(out_table[read_shots])) { + read_shots++; + } + return read_shots; +} + /// 01 format MeasureRecordReaderFormat01::MeasureRecordReaderFormat01( @@ -144,6 +153,25 @@ bool MeasureRecordReaderFormat01::expects_empty_serialized_data_for_each_shot() return false; } +size_t MeasureRecordReaderFormat01::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t max_shots = out_table.num_minor_bits_padded(); + size_t read_shots = 0; + while (read_shots < max_shots) { + bool more = start_and_read_entire_record_helper( + [&](size_t k) { + out_table[k][read_shots] &= 0; + }, + [&](size_t k) { + out_table[k][read_shots] |= 1; + }); + if (!more) { + break; + } + read_shots++; + } + return read_shots; +} + /// B8 format MeasureRecordReaderFormatB8::MeasureRecordReaderFormatB8( @@ -168,6 +196,29 @@ bool MeasureRecordReaderFormatB8::start_and_read_entire_record(simd_bits_range_r return true; } +size_t MeasureRecordReaderFormatB8::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t max_shots = out_table.num_minor_bits_padded(); + size_t n = bits_per_record(); + if (n == 0) { + return 0; // Ambiguous when the data ends. Stop as early as possible. + } + for (size_t read_shots = 0; read_shots < max_shots; read_shots++) { + for (size_t bit = 0; bit < n; bit += 8) { + int c = getc(in); + if (c == EOF) { + if (bit == 0) { + return read_shots; + } + throw std::invalid_argument("b8 data ended in middle of record."); + } + for (size_t b = 0; b < 8 && bit + b < n; b++) { + out_table[bit + b][read_shots] = ((c >> b) & 1) != 0; + } + } + } + return max_shots; +} + bool MeasureRecordReaderFormatB8::start_and_read_entire_record(SparseShot &cleared_out) { size_t n = bits_per_record(); size_t nb = (n + 7) >> 3; @@ -236,6 +287,22 @@ bool MeasureRecordReaderFormatHits::expects_empty_serialized_data_for_each_shot( return false; } +size_t MeasureRecordReaderFormatHits::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t max_shots = out_table.num_minor_bits_padded(); + size_t read_shots = 0; + out_table.clear(); + while (read_shots < max_shots) { + bool more = start_and_read_entire_record_helper([&](size_t bit_index) { + out_table[bit_index][read_shots] |= 1; + }); + if (!more) { + break; + } + read_shots++; + } + return read_shots; +} + /// R8 format MeasureRecordReaderFormatR8::MeasureRecordReaderFormatR8( @@ -262,6 +329,22 @@ bool MeasureRecordReaderFormatR8::expects_empty_serialized_data_for_each_shot() return false; } +size_t MeasureRecordReaderFormatR8::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t max_shots = out_table.num_minor_bits_padded(); + size_t read_shots = 0; + out_table.clear(); + while (read_shots < max_shots) { + bool more = start_and_read_entire_record_helper([&](size_t bit_index) { + out_table[bit_index][read_shots] |= 1; + }); + if (!more) { + break; + } + read_shots++; + } + return read_shots; +} + /// DETS format bool MeasureRecordReaderFormatDets::start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) { @@ -291,16 +374,39 @@ bool MeasureRecordReaderFormatDets::expects_empty_serialized_data_for_each_shot( return false; } +size_t MeasureRecordReaderFormatDets::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t max_shots = out_table.num_minor_bits_padded(); + size_t read_shots = 0; + out_table.clear(); + while (read_shots < max_shots) { + bool more = start_and_read_entire_record_helper([&](size_t bit_index) { + out_table[bit_index][read_shots] |= 1; + }); + if (!more) { + break; + } + read_shots++; + } + return read_shots; +} + +/// PTB64 format + MeasureRecordReaderFormatPTB64::MeasureRecordReaderFormatPTB64( FILE *in, size_t num_measurements, size_t num_detectors, size_t num_observables) : MeasureRecordReader(num_measurements, num_detectors, num_observables), in(in), - buf((bits_per_record() + 63) / 64 * 64 * 64), + buf(0), num_unread_shots_in_buf(0) { } bool MeasureRecordReaderFormatPTB64::load_cache() { size_t n = bits_per_record(); + size_t expected_buf_bits = (n + 63) / 64 * 64 * 64; + if (buf.num_bits_padded() < expected_buf_bits) { + buf = simd_bits(expected_buf_bits); + } + size_t nb = bits_per_record() * (64 / 8); size_t nr = fread(buf.u8, 1, nb, in); if (nr == 0) { @@ -371,3 +477,83 @@ bool MeasureRecordReaderFormatPTB64::start_and_read_entire_record(SparseShot &cl bool MeasureRecordReaderFormatPTB64::expects_empty_serialized_data_for_each_shot() const { return bits_per_record() == 0; } + +size_t MeasureRecordReaderFormatPTB64::read_into_table_with_minor_shot_index(simd_bit_table &out_table) { + size_t n = bits_per_record(); + if (n == 0) { + return 0; // Ambiguous when the data ends. Stop as early as possible. + } + size_t max_shots = out_table.num_minor_bits_padded(); + assert(max_shots % 64 == 0); + for (size_t shots_read = 0; shots_read < max_shots; shots_read += 64) { + for (size_t bit = 0; bit < n; bit++) { + size_t read = fread(&out_table[bit].u64[shots_read >> 6], 1, sizeof(uint64_t), in); + if (read != sizeof(uint64_t)) { + if (read == 0 && bit == 0) { + // End of file at a shot boundary. + return shots_read; + } else { + // Fragmented file. + throw std::invalid_argument("File ended in the middle of a ptb64 record."); + } + } + } + } + return max_shots; +} + +size_t MeasureRecordReaderFormatPTB64::read_into_table_with_major_shot_index(simd_bit_table &out_table) { + size_t n = bits_per_record(); + if (n == 0) { + return 0; // Ambiguous when the data ends. Stop as early as possible. + } + uint64_t buffer[64]; + size_t max_shots = out_table.num_minor_bits_padded(); + assert(max_shots % 64 == 0); + for (size_t shot = 0; shot < max_shots; shot += 64) { + for (size_t bit = 0; bit < n; bit += 64) { + for (size_t b = 0; b < 64; b++) { + if (bit + b >= n) { + buffer[b] = 0; + } else { + size_t read = fread(&buffer[b], 1, sizeof(uint64_t), in); + if (read != sizeof(uint64_t)) { + if (read == 0 && bit == 0 && b == 0) { + // End of file at a shot boundary. + return shot; + } else { + // Fragmented file. + throw std::invalid_argument("File ended in the middle of a ptb64 record."); + } + } + } + } + inplace_transpose_64x64(buffer); + for (size_t s = 0; s < 64; s++) { + out_table[shot + s].u64[bit >> 6] = buffer[s]; + } + } + } + return max_shots; +} + +size_t stim::read_file_data_into_shot_table( + FILE *in, + size_t max_shots, + size_t num_bits_per_shot, + SampleFormat format, + char dets_char, + simd_bit_table &out_table, + bool shots_is_major_index_of_out_table) { + auto reader = MeasureRecordReader::make( + in, + format, + dets_char == 'M' ? num_bits_per_shot : 0, + dets_char == 'D' ? num_bits_per_shot : 0, + dets_char == 'L' ? num_bits_per_shot : 0); + if (shots_is_major_index_of_out_table) { + return reader->read_into_table_with_major_shot_index(out_table); + } else { + return reader->read_into_table_with_minor_shot_index(out_table); + } +} diff --git a/src/stim/io/measure_record_reader.h b/src/stim/io/measure_record_reader.h index 572b982a2..7c2084927 100644 --- a/src/stim/io/measure_record_reader.h +++ b/src/stim/io/measure_record_reader.h @@ -123,6 +123,28 @@ struct MeasureRecordReader { /// std::invalid_argument: A record was only partially read. virtual bool start_and_read_entire_record(SparseShot &cleared_out) = 0; + /// Reads many records into a shot table. + /// + /// Args: + /// out_table: The table to write shots into. + /// Must have num_minor_bits >= bits_per_shot. + /// num_major_bits is max read shots. + /// + /// Returns: + /// The number of shots that were read. + virtual size_t read_into_table_with_major_shot_index(simd_bit_table &out_table); + + /// Reads many records into a shot table. + /// + /// Args: + /// out_table: The table to write shots into. + /// Must have num_major_bits >= bits_per_shot. + /// num_minor_bits is max read shots. + /// + /// Returns: + /// The number of shots that were read. + virtual size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) = 0; + protected: void move_obs_in_shots_to_mask_assuming_sorted(SparseShot &shot); }; @@ -139,6 +161,8 @@ struct MeasureRecordReaderFormatPTB64 : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_major_shot_index(simd_bit_table &out_table) override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; private: bool load_cache(); @@ -152,6 +176,7 @@ struct MeasureRecordReaderFormat01 : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; private: template @@ -199,6 +224,7 @@ struct MeasureRecordReaderFormatB8 : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; }; struct MeasureRecordReaderFormatHits : MeasureRecordReader { @@ -209,6 +235,7 @@ struct MeasureRecordReaderFormatHits : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; private: template @@ -246,6 +273,7 @@ struct MeasureRecordReaderFormatR8 : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; private: template @@ -290,6 +318,7 @@ struct MeasureRecordReaderFormatDets : MeasureRecordReader { bool start_and_read_entire_record(simd_bits_range_ref dirty_out_buffer) override; bool start_and_read_entire_record(SparseShot &cleared_out) override; bool expects_empty_serialized_data_for_each_shot() const override; + size_t read_into_table_with_minor_shot_index(simd_bit_table &out_table) override; private: template @@ -352,6 +381,15 @@ struct MeasureRecordReaderFormatDets : MeasureRecordReader { } }; +size_t read_file_data_into_shot_table( + FILE *in, + size_t max_shots, + size_t num_bits_per_shot, + SampleFormat format, + char dets_char, + simd_bit_table &out_table, + bool shots_is_major_index_of_out_table); + } // namespace stim #endif diff --git a/src/stim/io/measure_record_reader.test.cc b/src/stim/io/measure_record_reader.test.cc index 1a2129419..b5894cea0 100644 --- a/src/stim/io/measure_record_reader.test.cc +++ b/src/stim/io/measure_record_reader.test.cc @@ -698,3 +698,39 @@ TEST(MeasureRecordReader, start_and_read_entire_record_ptb64_sparse) { SparseShot discard; ASSERT_FALSE(reader->start_and_read_entire_record(discard)); } + +TEST(MeasureRecordReader, read_file_data_into_shot_table_vs_write_table) { + for (const auto &format_data : format_name_to_enum_map) { + SampleFormat format = format_data.second.id; + size_t num_shots = 500; + if (format == SAMPLE_FORMAT_PTB64) { + num_shots = 512 + 64; + } + size_t bits_per_shot = 1000; + + simd_bit_table expected(num_shots, bits_per_shot); + for (size_t shot = 0; shot < num_shots; shot++) { + expected[shot].randomize(bits_per_shot, SHARED_TEST_RNG()); + } + simd_bit_table expected_transposed = expected.transposed(); + + RaiiTempNamedFile tmp; + FILE *f = fopen(tmp.path.c_str(), "w"); + write_table_data(f, num_shots, bits_per_shot, simd_bits(0), expected_transposed, format, 'M', 'M', 0); + fclose(f); + + f = fopen(tmp.path.c_str(), "r"); + simd_bit_table output(num_shots, bits_per_shot); + read_file_data_into_shot_table(f, num_shots, bits_per_shot, format, 'M', output, true); + ASSERT_EQ(getc(f), EOF) << format_data.second.name << ", not transposed"; + fclose(f); + ASSERT_EQ(output, expected) << format_data.second.name << ", not transposed"; + + f = fopen(tmp.path.c_str(), "r"); + simd_bit_table output_transposed(bits_per_shot, num_shots); + read_file_data_into_shot_table(f, num_shots, bits_per_shot, format, 'M', output_transposed, false); + ASSERT_EQ(getc(f), EOF) << format_data.second.name << ", yes transposed"; + fclose(f); + ASSERT_EQ(output_transposed, expected_transposed) << format_data.second.name << ", yes transposed"; + } +} diff --git a/src/stim/io/raii_file.cc b/src/stim/io/raii_file.cc index ca558a189..4ba69cba2 100644 --- a/src/stim/io/raii_file.cc +++ b/src/stim/io/raii_file.cc @@ -18,7 +18,15 @@ using namespace stim; -RaiiFile::RaiiFile(const char *path, const char *mode) { +RaiiFile::RaiiFile(FILE *claim_ownership) : f(claim_ownership), responsible_for_closing(true) { +} + +RaiiFile::RaiiFile(RaiiFile &&other) noexcept : f(other.f), responsible_for_closing(other.responsible_for_closing) { + other.responsible_for_closing = false; + other.f = nullptr; +} + +RaiiFile::RaiiFile(const char *path, const char *mode) : f(nullptr), responsible_for_closing(true) { if (path == nullptr) { f = nullptr; return; @@ -40,8 +48,13 @@ RaiiFile::RaiiFile(const char *path, const char *mode) { } RaiiFile::~RaiiFile() { - if (f != nullptr) { + done(); +} + +void RaiiFile::done() { + if (f != nullptr && responsible_for_closing) { fclose(f); f = nullptr; + responsible_for_closing = false; } } diff --git a/src/stim/io/raii_file.h b/src/stim/io/raii_file.h index 443f83563..20cac3fca 100644 --- a/src/stim/io/raii_file.h +++ b/src/stim/io/raii_file.h @@ -21,10 +21,13 @@ namespace stim { struct RaiiFile { FILE* f; + bool responsible_for_closing; RaiiFile(const char* path, const char* mode); + RaiiFile(FILE* claim_ownership); RaiiFile(const RaiiFile& other) = delete; - RaiiFile(RaiiFile&& other) = delete; + RaiiFile(RaiiFile&& other) noexcept; ~RaiiFile(); + void done(); }; } // namespace stim diff --git a/src/stim/io/read_write.pybind.h b/src/stim/io/read_write.pybind.h index 7f919a7e4..47e20f2bb 100644 --- a/src/stim/io/read_write.pybind.h +++ b/src/stim/io/read_write.pybind.h @@ -24,6 +24,7 @@ namespace stim_pybind { stim::simd_bit_table numpy_array_to_transposed_simd_table( const pybind11::object &data, size_t expected_bits_per_shot, size_t *num_shots_out); + pybind11::object transposed_simd_bit_table_to_numpy( const stim::simd_bit_table &table, size_t bits_per_shot, size_t num_shots, bool bit_pack_result); diff --git a/src/stim/main_namespaced.cc b/src/stim/main_namespaced.cc index b632e331c..9e6130f9b 100644 --- a/src/stim/main_namespaced.cc +++ b/src/stim/main_namespaced.cc @@ -17,8 +17,10 @@ #include "stim/arg_parse.h" #include "stim/gen/circuit_gen_main.h" #include "stim/help.h" +#include "stim/io/raii_file.h" #include "stim/io/stim_data_formats.h" #include "stim/probability_util.h" +#include "stim/simulators/dem_sampler.h" #include "stim/simulators/detection_simulator.h" #include "stim/simulators/error_analyzer.h" #include "stim/simulators/error_matcher.h" @@ -82,40 +84,39 @@ int main_mode_detect(int argc, const char **argv) { find_argument("--shots", argc, argv) ? (uint64_t)find_int64_argument("--shots", 1, 0, INT64_MAX, argc, argv) : find_argument("--detect", argc, argv) ? (uint64_t)find_int64_argument("--detect", 1, 0, INT64_MAX, argc, argv) : 1; - if (num_shots == 0) { - return EXIT_SUCCESS; - } if (out_format.id == SAMPLE_FORMAT_DETS && !append_observables) { prepend_observables = true; } - FILE *in = find_open_file_argument("--in", stdin, "r", argc, argv); - FILE *out = find_open_file_argument("--out", stdout, "w", argc, argv); - FILE *obs_out = find_open_file_argument("--obs_out", stdout, "w", argc, argv); - if (obs_out == stdout) { - obs_out = nullptr; + RaiiFile in(find_open_file_argument("--in", stdin, "r", argc, argv)); + RaiiFile out(find_open_file_argument("--out", stdout, "w", argc, argv)); + RaiiFile obs_out(find_open_file_argument("--obs_out", stdout, "w", argc, argv)); + if (obs_out.f == stdout) { + obs_out.f = nullptr; } - auto circuit = Circuit::from_file(in); - if (in != stdin) { - fclose(in); + if (out.f == stdout) { + out.responsible_for_closing = false; + } + if (in.f == stdin) { + out.responsible_for_closing = false; + } + if (num_shots == 0) { + return EXIT_SUCCESS; } + + auto circuit = Circuit::from_file(in.f); + in.done(); auto rng = optionally_seeded_rng(argc, argv); detector_samples_out( circuit, num_shots, prepend_observables, append_observables, - out, + out.f, out_format.id, rng, - obs_out, + obs_out.f, obs_out_format.id); - if (obs_out != nullptr) { - fclose(obs_out); - } - if (out != stdout) { - fclose(out); - } return EXIT_SUCCESS; } @@ -291,6 +292,73 @@ int main_mode_repl(int argc, const char **argv) { return EXIT_SUCCESS; } +int main_mode_sample_dem(int argc, const char **argv) { + check_for_unknown_arguments( + { + "--seed", + "--shots", + "--out_format", + "--out", + "--in", + "--obs_out", + "--obs_out_format", + "--err_out", + "--err_out_format", + "--replay_err_in", + "--replay_err_in_format", + }, + {}, + "sample_dem", + argc, + argv); + const auto &out_format = find_enum_argument("--out_format", "01", format_name_to_enum_map, argc, argv); + const auto &obs_out_format = find_enum_argument("--obs_out_format", "01", format_name_to_enum_map, argc, argv); + const auto &err_out_format = find_enum_argument("--err_out_format", "01", format_name_to_enum_map, argc, argv); + const auto &err_in_format = find_enum_argument("--replay_err_in_format", "01", format_name_to_enum_map, argc, argv); + uint64_t num_shots = find_int64_argument("--shots", 1, 0, INT64_MAX, argc, argv); + + RaiiFile in(find_open_file_argument("--in", stdin, "r", argc, argv)); + RaiiFile out(find_open_file_argument("--out", stdout, "w", argc, argv)); + RaiiFile obs_out(find_open_file_argument("--obs_out", stdout, "w", argc, argv)); + RaiiFile err_out(find_open_file_argument("--err_out", stdout, "w", argc, argv)); + RaiiFile err_in(find_open_file_argument("--replay_err_in", stdin, "r", argc, argv)); + if (obs_out.f == stdout) { + obs_out.f = nullptr; + } + if (err_out.f == stdout) { + err_out.f = nullptr; + } + if (err_in.f == stdin) { + err_in.f = nullptr; + } + if (out.f == stdout) { + out.responsible_for_closing = false; + } + if (in.f == stdin) { + out.responsible_for_closing = false; + } + if (num_shots == 0) { + return EXIT_SUCCESS; + } + + auto dem = DetectorErrorModel::from_file(in.f); + in.done(); + + DemSampler sampler(std::move(dem), optionally_seeded_rng(argc, argv), 1024); + sampler.sample_write( + num_shots, + out.f, + out_format.id, + obs_out.f, + obs_out_format.id, + err_out.f, + err_out_format.id, + err_in.f, + err_in_format.id); + + return EXIT_SUCCESS; +} + int stim::main(int argc, const char **argv) { try { const char *mode = argc > 1 ? argv[1] : ""; @@ -298,7 +366,10 @@ int stim::main(int argc, const char **argv) { mode = ""; } auto is_mode = [&](const char *name) { - return find_argument(name, argc, argv) != nullptr || strcmp(mode, name + 2) == 0; + if (name[0] == '-') { + return find_argument(name, argc, argv) != nullptr || strcmp(mode, name + 2) == 0; + } + return strcmp(mode, name) == 0; }; if (is_mode("--help")) { @@ -307,6 +378,7 @@ int stim::main(int argc, const char **argv) { bool mode_repl = is_mode("--repl"); bool mode_sample = is_mode("--sample"); + bool mode_sample_dem = is_mode("sample_dem"); bool mode_detect = is_mode("--detect"); bool mode_analyze_errors = is_mode("--analyze_errors"); bool mode_gen = is_mode("--gen"); @@ -318,7 +390,7 @@ int stim::main(int argc, const char **argv) { mode_analyze_errors = true; } int modes_picked = - (mode_repl + mode_sample + mode_detect + mode_analyze_errors + mode_gen + mode_convert + + (mode_repl + mode_sample + mode_sample_dem + mode_detect + mode_analyze_errors + mode_gen + mode_convert + mode_explain_errors); if (modes_picked != 1) { std::cerr << "\033[31m"; @@ -353,6 +425,9 @@ int stim::main(int argc, const char **argv) { if (mode_explain_errors) { return main_mode_explain_errors(argc, argv); } + if (mode_sample_dem) { + return main_mode_sample_dem(argc, argv); + } throw std::out_of_range("Mode not handled."); } catch (const std::invalid_argument &ex) { diff --git a/src/stim/main_namespaced.test.cc b/src/stim/main_namespaced.test.cc index b83d7eaa7..5da043a4b 100644 --- a/src/stim/main_namespaced.test.cc +++ b/src/stim/main_namespaced.test.cc @@ -27,12 +27,7 @@ std::string execute(std::vector flags, const char *std_in_content) // Setup input. RaiiTempNamedFile raii_temp_file; if (std_in_content != nullptr) { - FILE *tmp_in = fdopen(raii_temp_file.descriptor, "w"); - if (tmp_in == nullptr) { - throw std::runtime_error("Failed to open temporary stdin file."); - } - fprintf(tmp_in, "%s", std_in_content); - fclose(tmp_in); + raii_temp_file.write_contents(std_in_content); flags.push_back("--in"); flags.push_back(raii_temp_file.path.data()); } @@ -829,8 +824,7 @@ TEST(main, detection_event_simulator_counts_measurements_correctly) { TEST(main, m2d) { RaiiTempNamedFile tmp; - FILE *f = fopen(tmp.path.data(), "w"); - fprintf(f, "%s", R"CIRCUIT( + tmp.write_contents(R"CIRCUIT( X 0 M 0 1 DETECTOR rec[-2] @@ -944,9 +938,7 @@ TEST(main, explain_errors) { ASSERT_EQ(execute({"explain_errors"}, ""), ""); RaiiTempNamedFile tmp; - FILE *f = fopen(tmp.path.data(), "w"); - fprintf(f, "error(1) D0\n"); - fclose(f); + tmp.write_contents("error(1) D0\n"); ASSERT_EQ( trim(execute({"explain_errors", "--dem_filter", tmp.path.data()}, R"input( @@ -969,3 +961,37 @@ ExplainedError { } )output")); } + +TEST(main, sample_dem) { + ASSERT_EQ(execute({"sample_dem"}, ""), "\n"); + + RaiiTempNamedFile obs_out; + + ASSERT_EQ( + trim(execute( + { + "sample_dem", + "--obs_out", + obs_out.path.data(), + "--out_format", + "01", + "--obs_out_format", + "01", + "--shots", + "5", + "--seed", + "0", + }, + R"input( + error(0) D0 + error(1) D1 L2 + )input")), + trim(R"output( +01 +01 +01 +01 +01 + )output")); + ASSERT_EQ(obs_out.read_contents(), "001\n001\n001\n001\n001\n"); +} diff --git a/src/stim/py/compiled_measurement_sampler.pybind.cc b/src/stim/py/compiled_measurement_sampler.pybind.cc index 4e0400867..fe61a224f 100644 --- a/src/stim/py/compiled_measurement_sampler.pybind.cc +++ b/src/stim/py/compiled_measurement_sampler.pybind.cc @@ -87,19 +87,20 @@ std::string CompiledMeasurementSampler::repr() const { return result.str(); } -pybind11::class_ pybind_compiled_measurement_sampler_class(pybind11::module &m) { +pybind11::class_ stim_pybind::pybind_compiled_measurement_sampler_class( + pybind11::module &m) { return pybind11::class_( m, "CompiledMeasurementSampler", "An analyzed stabilizer circuit whose measurements can be sampled quickly."); } -CompiledMeasurementSampler py_init_compiled_sampler( +CompiledMeasurementSampler stim_pybind::py_init_compiled_sampler( const Circuit &circuit, bool skip_reference_sample, const pybind11::object &seed) { simd_bits ref_sample = skip_reference_sample ? simd_bits(circuit.count_measurements()) : TableauSimulator::reference_sample_circuit(circuit); return CompiledMeasurementSampler(ref_sample, circuit, skip_reference_sample, make_py_seeded_rng(seed)); } -void pybind_compiled_measurement_sampler_methods(pybind11::class_ &c) { +void stim_pybind::pybind_compiled_measurement_sampler_methods(pybind11::class_ &c) { c.def( pybind11::init(&py_init_compiled_sampler), pybind11::arg("circuit"), diff --git a/src/stim/py/compiled_measurement_sampler.pybind.h b/src/stim/py/compiled_measurement_sampler.pybind.h index c98e80849..5bbf0f5df 100644 --- a/src/stim/py/compiled_measurement_sampler.pybind.h +++ b/src/stim/py/compiled_measurement_sampler.pybind.h @@ -22,6 +22,8 @@ #include "stim/circuit/circuit.h" #include "stim/mem/simd_bits.h" +namespace stim_pybind { + struct CompiledMeasurementSampler { const stim::simd_bits ref_sample; const stim::Circuit circuit; @@ -46,4 +48,6 @@ void pybind_compiled_measurement_sampler_methods(pybind11::class_ + +#include "stim/io/measure_record_reader.h" +#include "stim/io/measure_record_writer.h" +#include "stim/probability_util.h" + +using namespace stim; + +DemSampler::DemSampler(DetectorErrorModel init_model, std::mt19937_64 rng, size_t min_stripes) + : model(std::move(init_model)), + num_detectors(model.count_detectors()), + num_observables(model.count_observables()), + num_errors(model.count_errors()), + rng(rng), + det_buffer((size_t)num_detectors, min_stripes), + obs_buffer((size_t)num_observables, min_stripes), + err_buffer((size_t)num_errors, min_stripes), + num_stripes(det_buffer.num_minor_bits_padded()) { +} + +void DemSampler::resample(bool replay_errors) { + det_buffer.clear(); + obs_buffer.clear(); + if (!replay_errors) { + err_buffer.clear(); + } + size_t error_index = 0; + model.iter_flatten_error_instructions([&](const DemInstruction &op) { + simd_bits_range_ref err_row = err_buffer[error_index]; + if (!replay_errors) { + biased_randomize_bits((float)op.arg_data[0], err_row.u64, err_row.u64 + err_row.num_u64_padded(), rng); + } + for (const auto &t : op.target_data) { + if (t.is_relative_detector_id()) { + det_buffer[(size_t)t.raw_id()] ^= err_row; + } else if (t.is_observable_id()) { + obs_buffer[(size_t)t.raw_id()] ^= err_row; + } + } + error_index++; + }); +} + +void DemSampler::sample_write( + size_t num_shots, + FILE *det_out, + SampleFormat det_out_format, + FILE *obs_out, + SampleFormat obs_out_format, + FILE *err_out, + SampleFormat err_out_format, + FILE *err_in, + SampleFormat err_in_format) { + for (size_t k = 0; k < num_shots; k += num_stripes) { + size_t shots_left = std::min(num_stripes, num_shots - k); + + if (err_in != nullptr) { + size_t errors_read = read_file_data_into_shot_table( + err_in, shots_left, (size_t)num_errors, err_in_format, 'M', err_buffer, false); + if (errors_read != shots_left) { + throw std::invalid_argument("Expected more error data for the requested number of shots."); + } + } + resample(err_in != nullptr); + + if (err_out != nullptr) { + write_table_data( + err_out, shots_left, (size_t)num_errors, simd_bits(0), err_buffer, err_out_format, 'M', 'M', false); + } + + if (obs_out != nullptr) { + write_table_data( + obs_out, + shots_left, + (size_t)num_observables, + simd_bits(0), + obs_buffer, + obs_out_format, + 'L', + 'L', + false); + } + + if (det_out != nullptr) { + write_table_data( + det_out, shots_left, (size_t)num_detectors, simd_bits(0), det_buffer, det_out_format, 'D', 'D', false); + } + } +} diff --git a/src/stim/simulators/dem_sampler.h b/src/stim/simulators/dem_sampler.h new file mode 100644 index 000000000..4c46272fb --- /dev/null +++ b/src/stim/simulators/dem_sampler.h @@ -0,0 +1,74 @@ +/* + * Copyright 2021 Google LLC + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef _STIM_SIMULATORS_DEM_SAMPLER_H +#define _STIM_SIMULATORS_DEM_SAMPLER_H + +#include + +#include "stim/dem/detector_error_model.h" +#include "stim/io/stim_data_formats.h" +#include "stim/mem/simd_bit_table.h" + +namespace stim { + +/// Performs high performance bulk sampling of a detector error model. +struct DemSampler { + DetectorErrorModel model; + uint64_t num_detectors; + uint64_t num_observables; + uint64_t num_errors; + std::mt19937_64 rng; + // TODO: allow these buffers to be streamed instead of entirely stored in memory. + simd_bit_table det_buffer; + simd_bit_table obs_buffer; + simd_bit_table err_buffer; + size_t num_stripes; + + /// Compiles a sampler for the given detector error model. + DemSampler(DetectorErrorModel model, std::mt19937_64 rng, size_t min_stripes); + + /// Clears the buffers and refills them with sampled shot data. + void resample(bool replay_errors); + + /// Samples from the dem, writing results to files. + /// + /// Args: + /// num_shots: The number of samples to take. + /// det_out: Where to write detection event data. Set to nullptr to not write detection event data. + /// det_out_format: The format to write detection event data in. + /// obs_out: Where to write observable data. Set to nullptr to not write observable data. + /// obs_out_format: The format to write observable data in. + /// err_out: Where to write recorded error data. Set to nullptr to not write recorded error data. + /// err_out_format: The format to write error data in. + /// replay_err_in: If this argument is given a non-null file, error data will be read from that file + /// and replayed (instead of generating new errors randomly). + /// replay_err_in_format: The format to read recorded error data to replay in. + void sample_write( + size_t num_shots, + FILE *det_out, + SampleFormat det_out_format, + FILE *obs_out, + SampleFormat obs_out_format, + FILE *err_out, + SampleFormat err_out_format, + FILE *replay_err_in, + SampleFormat replay_err_in_format); +}; + +} // namespace stim + +#endif diff --git a/src/stim/simulators/dem_sampler.perf.cc b/src/stim/simulators/dem_sampler.perf.cc new file mode 100644 index 000000000..25348a794 --- /dev/null +++ b/src/stim/simulators/dem_sampler.perf.cc @@ -0,0 +1,41 @@ +// Copyright 2021 Google LLC +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "dem_sampler.h" + +#include "stim/benchmark_util.perf.h" +#include "stim/gen/gen_surface_code.h" +#include "stim/simulators/error_analyzer.h" + +using namespace stim; + +BENCHMARK(DemSampler_surface_code_rotated_memory_z_distance11_100rounds_1024stripes) { + auto params = CircuitGenParameters(100, 11, "rotated_memory_z"); + params.before_measure_flip_probability = 0.001; + params.after_reset_flip_probability = 0.001; + params.after_clifford_depolarization = 0.001; + auto circuit = generate_surface_code_circuit(params).circuit; + auto dem = ErrorAnalyzer::circuit_to_detector_error_model(circuit, true, true, false, false, false, false); + std::mt19937_64 rng(0); + DemSampler sampler(dem, std::mt19937_64(0), 1024); + size_t count = 0; + benchmark_go([&]() { + sampler.resample(false); + count += sampler.det_buffer[0].popcnt(); + count += sampler.obs_buffer[0].popcnt(); + }).goal_millis(35); + if (count == 0) { + std::cerr << "Data dependence."; + } +} diff --git a/src/stim/simulators/dem_sampler.pybind.cc b/src/stim/simulators/dem_sampler.pybind.cc new file mode 100644 index 000000000..10c8e8c03 --- /dev/null +++ b/src/stim/simulators/dem_sampler.pybind.cc @@ -0,0 +1,347 @@ +// Copyright 2021 Google LLC +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "stim/simulators/dem_sampler.pybind.h" + +#include "stim/io/raii_file.h" +#include "stim/io/read_write.pybind.h" +#include "stim/py/base.pybind.h" + +using namespace stim; +using namespace stim_pybind; + +RaiiFile optional_py_path_to_raii_file(const pybind11::object &obj, const char *mode) { + try { + auto path = pybind11::cast(obj); + return RaiiFile(path.data(), mode); + } catch (pybind11::cast_error &ex) { + } + + auto py_path = pybind11::module::import("pathlib").attr("Path"); + if (pybind11::isinstance(obj, py_path)) { + auto path = pybind11::cast(pybind11::str(obj)); + return RaiiFile(path.data(), mode); + } + + return RaiiFile(nullptr); +} + +pybind11::object dem_sampler_py_sample( + DemSampler &self, size_t shots, bool bit_packed, bool return_errors, pybind11::object &recorded_errors_to_replay) { + bool replay = !recorded_errors_to_replay.is_none(); + if (replay && min_bits_to_num_bits_padded(shots) != self.num_stripes) { + DemSampler perfect_size(self.model, std::move(self.rng), shots); + auto result = dem_sampler_py_sample(perfect_size, shots, bit_packed, return_errors, recorded_errors_to_replay); + self.rng = std::move(perfect_size.rng); + return result; + } + + if (replay) { + size_t out_shots; + simd_bit_table converted = + numpy_array_to_transposed_simd_table(recorded_errors_to_replay, self.num_errors, &out_shots); + if (out_shots != shots) { + throw std::invalid_argument("recorded_errors_to_replay.shape[0] != shots"); + } + assert(converted.num_minor_bits_padded() == self.err_buffer.num_minor_bits_padded()); + assert(converted.num_major_bits_padded() == self.err_buffer.num_major_bits_padded()); + self.err_buffer = std::move(converted); + } + + self.resample(replay); + + pybind11::object err_out = pybind11::none(); + if (return_errors) { + err_out = transposed_simd_bit_table_to_numpy(self.err_buffer, self.num_errors, shots, bit_packed); + } + pybind11::object det_out = + transposed_simd_bit_table_to_numpy(self.det_buffer, self.num_detectors, shots, bit_packed); + pybind11::object obs_out = + transposed_simd_bit_table_to_numpy(self.obs_buffer, self.num_observables, shots, bit_packed); + return pybind11::make_tuple(det_out, obs_out, err_out); +} + +pybind11::class_ stim_pybind::pybind_dem_sampler(pybind11::module &m) { + return pybind11::class_( + m, + "CompiledDemSampler", + clean_doc_string(u8R"DOC( + A helper class for efficiently sampler from a detector error model. + + Examples: + >>> import stim + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + )DOC") + .data()); +} + +void stim_pybind::pybind_dem_sampler_after_types_all_defined( + pybind11::module &m, pybind11::class_ &c) { + c.def( + "sample", + &dem_sampler_py_sample, + pybind11::arg("shots"), + pybind11::kw_only(), + pybind11::arg("bit_packed") = false, + pybind11::arg("return_errors") = false, + pybind11::arg("recorded_errors_to_replay") = pybind11::none(), + clean_doc_string(u8R"DOC( + @signature def sample(self, shots: int, *, bit_packed: bool = False, return_errors: bool = False, recorded_errors_to_replay: Optional[np.ndarray] = None) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + Samples the detector error model's error mechanisms to produce sample data. + + Args: + shots: The number of times to sample from the model. + bit_packed: Defaults to false. + False: the returned numpy arrays have dtype=np.bool8. + True: the returned numpy arrays have dtype=np.uint8 and pack 8 bits into each byte. + + Setting this to True is equivalent to running np.packbits(data, endian='little', axis=1) + on each output value, but has the performance benefit of the data never being expanded + into an unpacked form. + return_errors: Defaults to False. + False: the first entry of the returned tuple is None. + True: the first entry of the returned tuple is a numpy array recording which errors were sampled. + recorded_errors_to_replay: Defaults to None, meaning sample errors randomly. + If not None, this is expected to be a 2d numpy array specifying which errors to apply (e.g. one + returned from a previous call to the sample method). The array must have + dtype=np.bool8 and shape=(num_shots, num_errors) or + dtype=np.uint8 and shape=(num_shots, math.ceil(num_errors / 8)). + + Returns: + A tuple (detector_data, obs_data, error_data). + + Assuming bit_packed is False and return_errors is True: + If error_data[s, k] is True, then the error with index k fired in the shot with index s. + If detector_data[s, k] is True, then the detector with index k ended up flipped in the shot with index s. + If obs_data[s, k] is True, then the observable with index k ended up flipped in the shot with index s. + + The dtype and shape of the data depends on the arguments: + if bit_packed: + detector_data.shape == (num_shots, num_detectors) + detector_data.dtype == np.bool8 + obs_data.shape == (num_shots, num_observables) + obs_data.dtype == np.bool8 + if return_errors: + error_data.shape = (num_shots, num_errors) + error_data.dtype = np.bool8 + else: + error_data is None + else: + detector_data.shape == (num_shots, math.ceil(num_detectors / 8)) + detector_data.dtype == np.uint8 + obs_data.shape == (num_shots, math.ceil(num_observables / 8)) + obs_data.dtype == np.uint8 + if return_errors: + error_data.shape = (num_shots, math.ceil(num_errors / 8)) + error_data.dtype = np.uint8 + else: + error_data is None + + Note that bit packing is done using little endian order on the last axis + (i.e. like `np.packbits(data, endian='little', axis=1)`). + + Examples: + >>> import stim + >>> import numpy as np + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(1) D1 D2 L0 + ... ''') + >>> sampler = dem.compile_sampler() + + >>> # Taking samples. + >>> det_data, obs_data, err_data_not_requested = sampler.sample(shots=4) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data_not_requested is None + True + + >>> # Recording errors. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True) + >>> det_data + array([[False, True, True], + [False, True, True], + [False, True, True], + [False, True, True]]) + >>> obs_data + array([[ True], + [ True], + [ True], + [ True]]) + >>> err_data + array([[False, True], + [False, True], + [False, True], + [False, True]]) + + >>> # Bit packing. + >>> det_data, obs_data, err_data = sampler.sample(shots=4, return_errors=True, bit_packed=True) + >>> det_data + array([[6], + [6], + [6], + [6]], dtype=uint8) + >>> obs_data + array([[1], + [1], + [1], + [1]], dtype=uint8) + >>> err_data + array([[2], + [2], + [2], + [2]], dtype=uint8) + + >>> # Recording and replaying errors. + >>> noisy_dem = stim.DetectorErrorModel(''' + ... error(0.125) D0 + ... error(0.25) D1 + ... ''') + >>> noisy_sampler = noisy_dem.compile_sampler() + >>> det_data, obs_data, err_data = noisy_sampler.sample(shots=100, return_errors=True) + >>> replay_det_data, replay_obs_data, _ = noisy_sampler.sample(shots=100, recorded_errors_to_replay=err_data) + >>> np.array_equal(det_data, replay_det_data) + True + >>> np.array_equal(obs_data, replay_obs_data) + True + )DOC") + .data()); + + c.def( + "sample_write", + [](DemSampler &self, + size_t shots, + pybind11::object &det_out_file, + const std::string &det_out_format, + pybind11::object &obs_out_file, + const std::string &obs_out_format, + pybind11::object &err_out_file, + const std::string &err_out_format, + pybind11::object &replay_err_in_file, + const std::string &replay_err_in_format) { + RaiiFile fd = optional_py_path_to_raii_file(det_out_file, "w"); + RaiiFile fo = optional_py_path_to_raii_file(obs_out_file, "w"); + RaiiFile feo = optional_py_path_to_raii_file(err_out_file, "w"); + RaiiFile fei = optional_py_path_to_raii_file(replay_err_in_file, "r"); + self.sample_write( + shots, + fd.f, + format_to_enum(det_out_format), + fo.f, + format_to_enum(obs_out_format), + feo.f, + format_to_enum(err_out_format), + fei.f, + format_to_enum(replay_err_in_format)); + }, + pybind11::arg("shots"), + pybind11::kw_only(), + pybind11::arg("det_out_file"), + pybind11::arg("det_out_format") = "01", + pybind11::arg("obs_out_file"), + pybind11::arg("obs_out_format") = "01", + pybind11::arg("err_out_file") = pybind11::none(), + pybind11::arg("err_out_format") = "01", + pybind11::arg("replay_err_in_file") = pybind11::none(), + pybind11::arg("replay_err_in_format") = "01", + clean_doc_string(u8R"DOC( + @signature def sample_write(self, shots: int, *, det_out_file: Union[None, str, pathlib.Path], det_out_format: str = "01", obs_out_file: Union[None, str, pathlib.Path], obs_out_format: str = "01", err_out_file: Union[None, str, pathlib.Path] = None, err_out_format: str = "01", replay_err_in_file: Union[None, str, pathlib.Path] = None, replay_err_in_format: str = "01") -> None: + Samples the detector error model and writes the results to disk. + + Args: + shots: The number of times to sample from the model. + det_out_file: Where to write detection event data. + If None: detection event data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + det_out_format: The format to write the detection event data in (e.g. "01" or "b8"). + obs_out_file: Where to write observable flip data. + If None: observable flip data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + obs_out_format: The format to write the observable flip data in (e.g. "01" or "b8"). + err_out_file: Where to write errors-that-occurred data. + If None: errors-that-occurred data is not written. + If str or pathlib.Path: opens and overwrites the file at the given path. + NOT IMPLEMENTED: io.IOBase + err_out_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + replay_err_in_file: If this is specified, errors are replayed from data instead of generated randomly. + If None: errors are generated randomly according to the probabilities in the detector error model. + If str or pathlib.Path: the file at the given path is opened and errors-to-apply data is read from there. + NOT IMPLEMENTED: io.IOBase + replay_err_in_format: The format to write the errors-that-occurred data in (e.g. "01" or "b8"). + + Returns: + Nothing. Results are written to disk. + + Examples: + >>> import stim + >>> import tempfile + >>> import pathlib + >>> dem = stim.DetectorErrorModel(''' + ... error(0) D0 + ... error(0) D1 + ... error(0) D0 + ... error(1) D1 D2 L0 + ... error(0) D0 + ... ''') + >>> sampler = dem.compile_sampler() + >>> with tempfile.TemporaryDirectory() as d: + ... d = pathlib.Path(d) + ... sampler.sample_write( + ... shots=1, + ... det_out_file=d / 'dets.01', + ... det_out_format='01', + ... obs_out_file=d / 'obs.01', + ... obs_out_format='01', + ... err_out_file=d / 'err.hits', + ... err_out_format='hits', + ... ) + ... with open(d / 'dets.01') as f: + ... assert f.read() == "011\n" + ... with open(d / 'obs.01') as f: + ... assert f.read() == "1\n" + ... with open(d / 'err.hits') as f: + ... assert f.read() == "3\n" + )DOC") + .data()); +} diff --git a/src/stim/simulators/dem_sampler.pybind.h b/src/stim/simulators/dem_sampler.pybind.h new file mode 100644 index 000000000..627c57d86 --- /dev/null +++ b/src/stim/simulators/dem_sampler.pybind.h @@ -0,0 +1,29 @@ +// Copyright 2021 Google LLC +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef _STIM_SIMULATORS_DEM_SAMPLER_PYBIND_H +#define _STIM_SIMULATORS_DEM_SAMPLER_PYBIND_H + +#include + +#include "stim/simulators/dem_sampler.h" + +namespace stim_pybind { + +pybind11::class_ pybind_dem_sampler(pybind11::module &m); +void pybind_dem_sampler_after_types_all_defined(pybind11::module &m, pybind11::class_ &c); + +} // namespace stim_pybind + +#endif diff --git a/src/stim/simulators/dem_sampler.test.cc b/src/stim/simulators/dem_sampler.test.cc new file mode 100644 index 000000000..267e76891 --- /dev/null +++ b/src/stim/simulators/dem_sampler.test.cc @@ -0,0 +1,100 @@ +// Copyright 2021 Google LLC +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "stim/simulators/dem_sampler.h" + +#include "gtest/gtest.h" + +#include "stim/test_util.test.h" + +using namespace stim; + +TEST(DemSampler, basic_sizing) { + std::mt19937_64 irrelevant_rng(0); + DemSampler sampler(DetectorErrorModel(R"DEM()DEM"), irrelevant_rng, 700); + ASSERT_EQ(sampler.det_buffer.num_major_bits_padded(), 0); + ASSERT_EQ(sampler.obs_buffer.num_major_bits_padded(), 0); + ASSERT_GE(sampler.det_buffer.num_minor_bits_padded(), 700); + ASSERT_GE(sampler.obs_buffer.num_minor_bits_padded(), 700); + sampler.resample(false); + ASSERT_FALSE(sampler.obs_buffer.data.not_zero()); + ASSERT_FALSE(sampler.det_buffer.data.not_zero()); + + sampler = DemSampler( + DetectorErrorModel(R"DEM( + logical_observable L2000 + detector D1000 + )DEM"), + irrelevant_rng, + 200); + ASSERT_GE(sampler.det_buffer.num_major_bits_padded(), 1000); + ASSERT_GE(sampler.obs_buffer.num_major_bits_padded(), 2000); + ASSERT_GE(sampler.det_buffer.num_minor_bits_padded(), 200); + ASSERT_GE(sampler.obs_buffer.num_minor_bits_padded(), 200); + sampler.resample(false); + ASSERT_FALSE(sampler.obs_buffer.data.not_zero()); + ASSERT_FALSE(sampler.det_buffer.data.not_zero()); +} + +TEST(DemSampler, resample_basic_probabilities) { + DemSampler sampler( + DetectorErrorModel(R"DEM( + error(0) D0 + error(0.25) D1 L0 + error(0.5) D2 + error(0.75) D3 + error(1) D4 ^ D5 + )DEM"), + SHARED_TEST_RNG(), + 1000); + for (size_t k = 0; k < 2; k++) { + sampler.resample(false); + ASSERT_EQ(sampler.det_buffer[0].popcnt(), 0); + ASSERT_GT(sampler.det_buffer[1].popcnt(), 0); + ASSERT_LT(sampler.det_buffer[1].popcnt(), 500); + ASSERT_GT(sampler.det_buffer[2].popcnt(), 250); + ASSERT_LT(sampler.det_buffer[2].popcnt(), 750); + ASSERT_GT(sampler.det_buffer[3].popcnt(), 500); + ASSERT_LT(sampler.det_buffer[3].popcnt(), 1000); + ASSERT_EQ(sampler.det_buffer[4].popcnt(), sampler.det_buffer[4].num_bits_padded()); + + ASSERT_EQ(sampler.det_buffer[1], sampler.obs_buffer[0]); + ASSERT_EQ(sampler.det_buffer[4], sampler.det_buffer[5]); + } +} + +TEST(DemSampler, resample_combinations) { + DemSampler sampler( + DetectorErrorModel(R"DEM( + error(0.1) D0 D1 + error(0.2) D1 D2 + error(0.3) D2 D0 + )DEM"), + SHARED_TEST_RNG(), + 1000); + for (size_t k = 0; k < 2; k++) { + sampler.resample(false); + ASSERT_GT(sampler.det_buffer[0].popcnt(), 340 - 100); + ASSERT_LT(sampler.det_buffer[0].popcnt(), 340 + 100); + ASSERT_GT(sampler.det_buffer[1].popcnt(), 260 - 100); + ASSERT_LT(sampler.det_buffer[1].popcnt(), 260 + 100); + ASSERT_GT(sampler.det_buffer[2].popcnt(), 380 - 100); + ASSERT_LT(sampler.det_buffer[2].popcnt(), 380 + 100); + + simd_bits total = sampler.det_buffer[0]; + total ^= sampler.det_buffer[1]; + total ^= sampler.det_buffer[2]; + ASSERT_FALSE(total.not_zero()); + } +} diff --git a/src/stim/simulators/dem_sampler_pybind_test.py b/src/stim/simulators/dem_sampler_pybind_test.py new file mode 100644 index 000000000..5f1f109b6 --- /dev/null +++ b/src/stim/simulators/dem_sampler_pybind_test.py @@ -0,0 +1,73 @@ +import numpy as np +import pathlib +import pytest +import stim +import tempfile + + +@pytest.mark.parametrize("bit_packed", [False, True]) +def test_dem_sampler_sample(bit_packed: bool): + noisy_dem = stim.DetectorErrorModel(""" + error(0.125) D0 + error(0.25) D1 + """) + noisy_sampler = noisy_dem.compile_sampler() + det_data, obs_data, err_data = noisy_sampler.sample(shots=100, return_errors=True, bit_packed=bit_packed) + replay_det_data, replay_obs_data, _ = noisy_sampler.sample(shots=100, recorded_errors_to_replay=err_data, bit_packed=bit_packed) + np.testing.assert_array_equal(det_data, replay_det_data) + np.testing.assert_array_equal(obs_data, replay_obs_data) + + +def test_dem_sampler_sampler_write(): + dem = stim.DetectorErrorModel(''' + error(0) D0 + error(0) D1 + error(0) D0 + error(1) D1 D2 L0 + error(0) D0 + ''') + sampler = dem.compile_sampler() + with tempfile.TemporaryDirectory() as d: + d = pathlib.Path(d) + sampler.sample_write( + shots=1, + det_out_file=d / 'dets.01', + det_out_format='01', + obs_out_file=d / 'obs.01', + obs_out_format='01', + err_out_file=d / 'err.hits', + err_out_format='hits', + ) + with open(d / 'dets.01') as f: + assert f.read() == "011\n" + with open(d / 'obs.01') as f: + assert f.read() == "1\n" + with open(d / 'err.hits') as f: + assert f.read() == "3\n" + + sampler = stim.DetectorErrorModel(''' + error(1) D0 # this should be overridden by the replay. + error(1) D1 + error(1) D0 + error(1) D1 D2 L0 + error(1) D0 + ''').compile_sampler() + sampler.sample_write( + shots=1, + det_out_file=d / 'dets.01', + det_out_format='01', + obs_out_file=d / 'obs.01', + obs_out_format='01', + err_out_file=d / 'err2.01', + err_out_format='01', + replay_err_in_file=d / 'err.hits', + replay_err_in_format='hits', + ) + with open(d / 'dets.01') as f: + assert f.read() == "011\n" + with open(d / 'obs.01') as f: + assert f.read() == "1\n" + with open(d / 'err.hits') as f: + assert f.read() == "3\n" + with open(d / 'err2.01') as f: + assert f.read() == "00010\n" diff --git a/src/stim/simulators/matched_error.pybind.cc b/src/stim/simulators/matched_error.pybind.cc index 02bb0282b..c6b3b822e 100644 --- a/src/stim/simulators/matched_error.pybind.cc +++ b/src/stim/simulators/matched_error.pybind.cc @@ -632,7 +632,7 @@ void pybind_MatchedError(pybind11::module &m) { c.def("__str__", &ExplainedError::str); } -void pybind_matched_error(pybind11::module &m) { +void stim_pybind::pybind_matched_error(pybind11::module &m) { pybind_CircuitErrorLocationStackFrame(m); pybind_GateTargetWithCoords(m); pybind_DemTargetWithCoords(m); diff --git a/src/stim/simulators/matched_error.pybind.h b/src/stim/simulators/matched_error.pybind.h index 2b9dfd404..beac30088 100644 --- a/src/stim/simulators/matched_error.pybind.h +++ b/src/stim/simulators/matched_error.pybind.h @@ -17,6 +17,10 @@ #include +namespace stim_pybind { + void pybind_matched_error(pybind11::module &m); +} + #endif diff --git a/src/stim/simulators/tableau_simulator.pybind.cc b/src/stim/simulators/tableau_simulator.pybind.cc index 0b5acc611..c3bb05683 100644 --- a/src/stim/simulators/tableau_simulator.pybind.cc +++ b/src/stim/simulators/tableau_simulator.pybind.cc @@ -116,7 +116,7 @@ TempViewableData args_to_target_pairs(PyTableauSimulator &self, const pybind11:: return result; } -void pybind_tableau_simulator(pybind11::module &m) { +void stim_pybind::pybind_tableau_simulator(pybind11::module &m) { auto c = pybind11::class_( m, "TableauSimulator", diff --git a/src/stim/simulators/tableau_simulator.pybind.h b/src/stim/simulators/tableau_simulator.pybind.h index feefc07dd..4d3d88612 100644 --- a/src/stim/simulators/tableau_simulator.pybind.h +++ b/src/stim/simulators/tableau_simulator.pybind.h @@ -19,6 +19,8 @@ #include "stim/simulators/tableau_simulator.h" +namespace stim_pybind { + struct PyTableauSimulator : stim::TableauSimulator { std::shared_ptr rng_reference; explicit PyTableauSimulator(std::shared_ptr rng); @@ -26,4 +28,6 @@ struct PyTableauSimulator : stim::TableauSimulator { void pybind_tableau_simulator(pybind11::module &m); +} // namespace stim_pybind + #endif diff --git a/src/stim/stabilizers/pauli_string.pybind.cc b/src/stim/stabilizers/pauli_string.pybind.cc index fb2b9f866..2ea08e1a9 100644 --- a/src/stim/stabilizers/pauli_string.pybind.cc +++ b/src/stim/stabilizers/pauli_string.pybind.cc @@ -215,7 +215,7 @@ PyPauliString PyPauliString::from_text(const char *text) { return value; } -void pybind_pauli_string(pybind11::module &m) { +void stim_pybind::pybind_pauli_string(pybind11::module &m) { auto c = pybind11::class_( m, "PauliString", diff --git a/src/stim/stabilizers/pauli_string.pybind.h b/src/stim/stabilizers/pauli_string.pybind.h index ce1f06e62..7b3ed5229 100644 --- a/src/stim/stabilizers/pauli_string.pybind.h +++ b/src/stim/stabilizers/pauli_string.pybind.h @@ -20,6 +20,8 @@ #include "stim/stabilizers/pauli_string.h" +namespace stim_pybind { + struct PyPauliString { stim::PauliString value; bool imag; @@ -52,4 +54,6 @@ struct PyPauliString { void pybind_pauli_string(pybind11::module &m); +} // namespace stim_pybind + #endif diff --git a/src/stim/stabilizers/tableau.pybind.cc b/src/stim/stabilizers/tableau.pybind.cc index 2df83a7c1..fb60b652e 100644 --- a/src/stim/stabilizers/tableau.pybind.cc +++ b/src/stim/stabilizers/tableau.pybind.cc @@ -25,7 +25,7 @@ using namespace stim; using namespace stim_pybind; -void pybind_tableau(pybind11::module &m) { +void stim_pybind::pybind_tableau(pybind11::module &m) { auto c = pybind11::class_( m, "Tableau", diff --git a/src/stim/stabilizers/tableau.pybind.h b/src/stim/stabilizers/tableau.pybind.h index 8d4465949..f66a594e6 100644 --- a/src/stim/stabilizers/tableau.pybind.h +++ b/src/stim/stabilizers/tableau.pybind.h @@ -17,6 +17,8 @@ #include +namespace stim_pybind { void pybind_tableau(pybind11::module &m); +} #endif \ No newline at end of file diff --git a/src/stim/stabilizers/tableau_iter.pybind.cc b/src/stim/stabilizers/tableau_iter.pybind.cc index 2844fdf01..0bd5893a8 100644 --- a/src/stim/stabilizers/tableau_iter.pybind.cc +++ b/src/stim/stabilizers/tableau_iter.pybind.cc @@ -19,7 +19,7 @@ using namespace stim; using namespace stim_pybind; -pybind11::class_ pybind_tableau_iter(pybind11::module &m) { +pybind11::class_ stim_pybind::pybind_tableau_iter(pybind11::module &m) { auto c = pybind11::class_( m, "TableauIterator", @@ -39,7 +39,8 @@ pybind11::class_ pybind_tableau_iter(pybind11::module &m) { return c; } -void pybind_tableau_iter_after_types_all_defined(pybind11::module &m, pybind11::class_ &c) { +void stim_pybind::pybind_tableau_iter_after_types_all_defined( + pybind11::module &m, pybind11::class_ &c) { c.def( "__iter__", [](TableauIterator &self) -> TableauIterator { diff --git a/src/stim/stabilizers/tableau_iter.pybind.h b/src/stim/stabilizers/tableau_iter.pybind.h index 0610735c4..b07e73424 100644 --- a/src/stim/stabilizers/tableau_iter.pybind.h +++ b/src/stim/stabilizers/tableau_iter.pybind.h @@ -19,7 +19,9 @@ #include "stim/stabilizers/tableau_iter.h" +namespace stim_pybind { pybind11::class_ pybind_tableau_iter(pybind11::module &m); void pybind_tableau_iter_after_types_all_defined(pybind11::module &m, pybind11::class_ &c); +} // namespace stim_pybind #endif \ No newline at end of file diff --git a/src/stim/test_util.test.cc b/src/stim/test_util.test.cc index dad9bcc7b..d5932c52f 100644 --- a/src/stim/test_util.test.cc +++ b/src/stim/test_util.test.cc @@ -57,3 +57,31 @@ RaiiTempNamedFile::~RaiiTempNamedFile() { path = ""; } } + +std::string RaiiTempNamedFile::read_contents() { + FILE *f = fopen(path.c_str(), "r"); + if (f == nullptr) { + throw std::runtime_error("Failed to open temp named file " + path); + } + std::string result; + while (true) { + int c = getc(f); + if (c == EOF) { + break; + } + result.push_back(c); + } + fclose(f); + return result; +} + +void RaiiTempNamedFile::write_contents(const std::string &contents) { + FILE *f = fopen(path.c_str(), "w"); + if (f == nullptr) { + throw std::runtime_error("Failed to open temp named file " + path); + } + for (char c : contents) { + putc(c, f); + } + fclose(f); +} diff --git a/src/stim/test_util.test.h b/src/stim/test_util.test.h index 3e4b631da..e3c8dd563 100644 --- a/src/stim/test_util.test.h +++ b/src/stim/test_util.test.h @@ -32,6 +32,8 @@ struct RaiiTempNamedFile { std::string path; RaiiTempNamedFile(); ~RaiiTempNamedFile(); + std::string read_contents(); + void write_contents(const std::string &contents); }; #endif