Source code for pennylane.devices.default_clifford
# Copyright 2018-2024 Xanadu Quantum Technologies Inc.# 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."""This module contains the Clifford simulator using ``stim``."""importconcurrent.futuresfromcollections.abcimportSequencefromdataclassesimportreplacefromfunctoolsimportpartialfromtypingimportUnionimportnumpyasnpimportpennylaneasqmlfrompennylane.measurementsimport(ClassicalShadowMP,CountsMP,DensityMatrixMP,ExpectationMP,MutualInfoMP,ProbabilityMP,PurityMP,SampleMP,ShadowExpvalMP,StateMP,VarianceMP,VnEntropyMP,)frompennylane.ops.qubit.observablesimportBasisStateProjectorfrompennylane.tapeimportQuantumScript,QuantumScriptBatchfrompennylane.transformsimportconvert_to_numpy_parametersfrompennylane.transforms.coreimportTransformProgramfrompennylane.typingimportResult,ResultBatchfrom.importDevicefrom.default_qubitimportaccepted_sample_measurementfrom.execution_configimportDefaultExecutionConfig,ExecutionConfigfrom.modifiersimportsimulator_tracking,single_tape_supportfrom.preprocessimport(decompose,null_postprocessing,validate_adjoint_trainable_params,validate_device_wires,validate_measurements,validate_multiprocessing_workers,validate_observables,)has_stim=Truetry:importstimexcept(ModuleNotFoundError,ImportError)asimport_error:# pragma: no coverhas_stim=False# Updated observable list_OBSERVABLES_MAP={"PauliX","PauliY","PauliZ","Hermitian","Identity","Projector","LinearCombination","Sum","SProd","Prod",}# Clifford gates and error channels_OPERATIONS_MAP={"Identity":"I","PauliX":"X","PauliY":"Y","PauliZ":"Z","Hadamard":"H","S":"S","Adjoint(S)":"S_DAG","SX":"SX","Adjoint(SX)":"SX_DAG","CNOT":"CNOT","SWAP":"SWAP","ISWAP":"ISWAP","Adjoint(ISWAP)":"ISWAP_DAG","CY":"CY","CZ":"CZ","GlobalPhase":None,"BasisState":None,"StatePrep":None,"Snapshot":None,"Barrier":None,"PauliError":"CORRELATED_ERROR","BitFlip":"X_ERROR","PhaseFlip":"Z_ERROR","DepolarizingChannel":"DEPOLARIZE1",}
[docs]defoperation_stopping_condition(op:qml.operation.Operator)->bool:"""Specifies whether an operation is accepted by ``DefaultClifford``."""returnop.namein_OPERATIONS_MAP
[docs]defobservable_stopping_condition(obs:qml.operation.Operator)->bool:"""Specifies whether an observable is accepted by ``DefaultClifford``."""returnobs.namein_OBSERVABLES_MAP
@qml.transformdef_validate_channels(tape,name="device"):"""Validates the channels for a circuit."""ifnottape.shotsandany(isinstance(op,qml.operation.Channel)foropintape.operations):raiseqml.DeviceError(f"Channel not supported on {name} without finite shots.")return(tape,),null_postprocessingdef_pl_op_to_stim(op):"""Convert PennyLane operation to a Stim operation"""try:stim_op=_OPERATIONS_MAP[op.name]stim_tg=map(str,op.wires)exceptKeyErrorase:raiseqml.DeviceError(f"Operator {op} not supported with default.clifford and does not provide a decomposition.")frome# Check if the operation is noisyifisinstance(op,qml.operation.Channel):stim_op+=f"({op.parameters[-1]})"# get the probabilityifop.name=="PauliError":stim_tg=[pauli+wireforpauli,wireinzip(op.parameters[0],stim_tg)]returnstim_op," ".join(stim_tg)def_pl_obs_to_linear_comb(meas_obs):"""Convert a PennyLane observable to a linear combination of Pauli strings"""meas_rep=meas_obs.pauli_rep# Use manual decomposition for enabling Hermitian and partial Projector supportifisinstance(meas_obs,(qml.Hermitian,BasisStateProjector)):meas_rep=qml.pauli_decompose(meas_obs.matrix(),wire_order=meas_obs.wires,pauli=True)# A Pauli decomposition for the observable must existifmeas_repisNone:raiseNotImplementedError(f"default.clifford doesn't support expectation value calculation with {type(meas_obs).__name__} at the moment.")coeffs,paulis=np.array(list(meas_rep.values())),[]meas_op_wires=list(meas_obs.wires)forpwinmeas_rep:p_wire,p_word=pw.keys(),pw.values()ifnotp_word:# empty pauli word correspond to identityr_wire,r_word=meas_op_wires[:1],"I"else:# reorder the wires based on original meas_op# reorder the pauli terms based on above.r_wire=sorted(p_wire,key=meas_op_wires.index)r_word="".join(map(pw.get,r_wire))paulis.append((r_word,r_wire))returncoeffs,paulis# pylint:disable = too-many-instance-attributes
[docs]@simulator_tracking@single_tape_supportclassDefaultClifford(Device):r"""A PennyLane device for fast simulation of Clifford circuits using `stim <https://github.com/quantumlib/stim/>`_. Args: wires (int, Iterable[Number, str]): Number of wires present on the device, or iterable that contains unique labels for the wires as numbers (i.e., ``[-1, 0, 2]``) or strings (``['aux_wire', 'q1', 'q2']``). Default ``None`` if not specified. shots (int, Sequence[int], Sequence[Union[int, Sequence[int]]]): The default number of shots to use in executions involving this device. check_clifford (bool): Check if all the gate operations in the circuits to be executed are Clifford. Default is ``True``. tableau (bool): Determines what should be returned when the device's state is computed with :func:`qml.state <pennylane.state>`. When ``True``, the device returns the final evolved Tableau. Alternatively, one may make it ``False`` to obtain the evolved state vector. Note that the latter might not be computationally feasible for larger qubit numbers. seed (Union[str, None, int, array_like[int], SeedSequence, BitGenerator, Generator]): A seed-like parameter matching that of ``seed`` for ``numpy.random.default_rng``, or a request to seed from numpy's global random number generator. The default, ``seed="global"`` pulls a seed from numpy's global generator. ``seed=None`` will pull a seed from the OS entropy. max_workers (int): A ``ProcessPoolExecutor`` executes tapes asynchronously using a pool of at most ``max_workers`` processes. If ``max_workers`` is ``None``, only the current process executes tapes. If you experience any issue, try setting ``max_workers`` to ``None``. **Example:** .. code-block:: python dev = qml.device("default.clifford", tableau=True) @qml.qnode(dev) def circuit(): qml.CNOT(wires=[0, 1]) qml.X(1) qml.ISWAP(wires=[0, 1]) qml.Hadamard(wires=[0]) return qml.state() >>> circuit() array([[0, 1, 1, 0, 0], [1, 0, 1, 1, 1], [0, 0, 0, 1, 0], [1, 0, 0, 1, 1]]) The devices execution pipeline can be investigated more closely with the following: .. code-block:: python num_qscripts = 5 qscripts = [ qml.tape.QuantumScript( [qml.Hadamard(wires=[0]), qml.CNOT(wires=[0, 1])], [qml.expval(qml.Z(0))] ) ] * num_qscripts >>> dev = DefaultClifford() >>> program, execution_config = dev.preprocess() >>> new_batch, post_processing_fn = program(qscripts) >>> results = dev.execute(new_batch, execution_config=execution_config) >>> post_processing_fn(results) (array(0), array(0), array(0), array(0), array(0)) .. details:: :title: Clifford Tableau :href: clifford-tableau-theory The device's internal state is represented by the following ``Tableau`` described in the `Sec. III, Aaronson & Gottesman (2004) <https://arxiv.org/abs/quant-ph/0406196>`_: .. math:: \begin{bmatrix} x_{11} & \cdots & x_{1n} & & z_{11} & \cdots & z_{1n} & &r_{1}\\ \vdots & \ddots & \vdots & & \vdots & \ddots & \vdots & &\vdots\\ x_{n1} & \cdots & x_{nn} & & z_{n1} & \cdots & z_{nn} & &r_{n}\\ & & & & & & & & \\ x_{\left( n+1\right) 1} & \cdots & x_{\left( n+1\right) n} & & z_{\left( n+1\right) 1} & \cdots & z_{\left( n+1\right) n} & & r_{n+1}\\ \vdots & \ddots & \vdots & & \vdots & \ddots & \vdots & & \vdots\\ x_{\left( 2n\right) 1} & \cdots & x_{\left( 2n\right) n} & & z_{\left( 2n\right) 1} & \cdots & z_{\left( 2n\right) n} & & r_{2n} \end{bmatrix} The tableau's first `n` rows represent a destabilizer generator, while the remaining `n` rows represent the stabilizer generators. The Pauli representation for all of these generators are described using the :mod:`binary vector <pennylane.pauli.binary_to_pauli>` made from the binary variables :math:`x_{ij},\ z_{ij}`, :math:`\forall i\in\left\{1,\ldots,2n\right\}, j\in\left\{1,\ldots,n\right\}` and they together form the complete Pauli group. Finally, the last column of the tableau, with binary variables :math:`r_{i},\ \forall i\in\left\{1,\ldots,2n\right\}`, denotes whether the phase is negative (:math:`r_i = 1`) or not, for each generator. Maintaining and working with this tableau representation instead of the complete state vector makes the calculations of increasingly large Clifford circuits more efficient on this device. .. details:: :title: Probabilities for Basis States :href: clifford-probabilities As the ``default.clifford`` device supports executing quantum circuits with a large number of qubits, the ability to compute the ``analytical`` probabilities for ``all`` computational basis states at once becomes computationally expensive and challenging as the system size increases. While we don't manually restrict users from doing so for any circuit, one can expect the underlying computation to reach its limit with ``20-24`` qubits on a typical consumer grade machine. As long as number of qubits are below this limit, one can simply use the :func:`qml.probs <pennylane.probs>` function with its usual arguments to compute probabilities for the complete computational basis states. We test this for a circuit that prepares the ``n``-qubit Greenberger-Horne-Zeilinger (GHZ) state. This means that the probabilities for the basis states :math:`|0\rangle^{\otimes n}` and :math:`|1\rangle^{\otimes n}` should be :math:`0.5`, and :math:`0.0` for the rest. .. code-block:: python import pennylane as qml import numpy as np dev = qml.device("default.clifford") num_wires = 3 @qml.qnode(dev) def circuit(): qml.Hadamard(wires=[0]) for idx in range(num_wires): qml.CNOT(wires=[idx, idx+1]) return qml.probs() >>> circuit() tensor([0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5], requires_grad=True) Once above the limit (or even otherwise), one can obtain the probability of a single target basis state by computing the expectation value of the corresponding projector using :mod:`qml.expval <pennylane.expval>` and :mod:`qml.Projector <pennylane.Projector>`. .. code-block:: python num_wires = 4 @qml.qnode(dev) def circuit(state): qml.Hadamard(wires=[0]) for idx in range(num_wires): qml.CNOT(wires=[idx, idx+1]) return qml.expval(qml.Projector(state, wires=range(num_wires))) >>> basis_states = np.array([[0, 0, 0, 0], [0, 1, 0, 1], [1, 0, 1, 0]]) >>> circuit(basis_states[0]) tensor(0.5, requires_grad=True) >>> circuit(basis_states[1]) tensor(0.0, requires_grad=True) >>> circuit(basis_states[2]) tensor(0.0, requires_grad=True) .. details:: :title: Error Channels :href: clifford-errors This device supports the finite-shot execution of quantum circuits with the following error channels that add Pauli noise, allowing for one to perform any sampling-based measurements. - *Multi-qubit Pauli errors:* :mod:`qml.PauliError <pennylane.PauliError>` - *Single-qubit depolarization errors:* :mod:`qml.DepolarizingChannel <pennylane.DepolarizingChannel>` - *Single-qubit flip errors:* :mod:`qml.BitFlip <pennylane.BitFlip>` and :mod:`qml.PhaseFlip <pennylane.PhaseFlip>` .. code-block:: python import pennylane as qml import numpy as np dev = qml.device("default.clifford", shots=1024, seed=42) num_wires = 3 @qml.qnode(dev) def circuit(): qml.Hadamard(wires=[0]) for idx in range(num_wires): qml.CNOT(wires=[idx, idx+1]) qml.BitFlip(0.2, wires=[1]) return qml.counts() >>> circuit() {'0000': 417, '0100': 95, '1011': 104, '1111': 408} .. details:: :title: Tracking :href: clifford-tracking ``DefaultClifford`` tracks: * ``executions``: the number of unique circuits that would be required on quantum hardware * ``shots``: the number of shots * ``resources``: the :class:`~.resource.Resources` for the executed circuit. * ``simulations``: the number of simulations performed. One simulation can cover multiple QPU executions, such as for non-commuting measurements and batched parameters. * ``batches``: The number of times :meth:`~.execute` is called. * ``results``: The results of each call of :meth:`~.execute`. .. details:: :title: Accelerate calculations with multiprocessing :href: clifford-multiprocessing See the details in :class:`~pennylane.devices.DefaultQubit`'s "Accelerate calculations with multiprocessing" section. Additional information regarding multiprocessing can be found in the `multiprocessing docs page <https://docs.python.org/3/library/multiprocessing.html#contexts-and-start-methods>`_. """@propertydefname(self):"""The name of the device."""return"default.clifford"# pylint:disable = too-many-argumentsdef__init__(self,wires=None,shots=None,check_clifford=True,tableau=True,seed="global",max_workers=None,)->None:ifnothas_stim:raiseImportError("This feature requires stim, a fast stabilizer circuit simulator. ""It can be installed with:\n\npip install stim")super().__init__(wires=wires,shots=shots)self._max_workers=max_workersself._check_clifford=check_cliffordself._tableau=tableauseed=np.random.randint(0,high=10000000)ifseed=="global"elseseedself._rng=np.random.default_rng(seed)self._debugger=Nonedef_setup_execution_config(self,execution_config:ExecutionConfig)->ExecutionConfig:"""This is a private helper for ``preprocess`` that sets up the execution config. Args: execution_config (ExecutionConfig) Returns: ExecutionConfig: a preprocessed execution config """updated_values={}ifexecution_config.gradient_method=="best":# pragma: no coverupdated_values["gradient_method"]=Noneupdated_values["use_device_jacobian_product"]=Falseifexecution_config.grad_on_executionisNone:updated_values["grad_on_execution"]=Falseupdated_values["device_options"]=dict(execution_config.device_options)# copyif"max_workers"notinupdated_values["device_options"]:updated_values["device_options"]["max_workers"]=self._max_workersif"rng"notinupdated_values["device_options"]:updated_values["device_options"]["rng"]=self._rngif"tableau"notinupdated_values["device_options"]:updated_values["device_options"]["tableau"]=self._tableaureturnreplace(execution_config,**updated_values)
[docs]defpreprocess(self,execution_config:ExecutionConfig=DefaultExecutionConfig,)->tuple[TransformProgram,ExecutionConfig]:"""This function defines the device transform program to be applied and an updated device configuration. Args: execution_config (Union[ExecutionConfig, Sequence[ExecutionConfig]]): A data structure describing the parameters needed to fully describe the execution. Returns: TransformProgram, ExecutionConfig: A transform program that when called returns QuantumTapes that the device can natively execute as well as a postprocessing function to be called after execution, and a configuration with unset specifications filled in. This device currently does not intrinsically support parameter broadcasting. """config=self._setup_execution_config(execution_config)transform_program=TransformProgram()transform_program.add_transform(validate_device_wires,self.wires,name=self.name)transform_program.add_transform(qml.defer_measurements,allow_postselect=False)# Perform circuit decomposition to the supported Clifford gate setifself._check_clifford:transform_program.add_transform(decompose,stopping_condition=operation_stopping_condition,name=self.name)transform_program.add_transform(_validate_channels,name=self.name)transform_program.add_transform(validate_measurements,sample_measurements=accepted_sample_measurement,name=self.name)transform_program.add_transform(validate_observables,stopping_condition=observable_stopping_condition,name=self.name)# Validate multi processingmax_workers=config.device_options.get("max_workers",self._max_workers)ifmax_workers:transform_program.add_transform(validate_multiprocessing_workers,max_workers,self)# Validate derivativestransform_program.add_transform(validate_adjoint_trainable_params)ifconfig.gradient_methodisnotNone:config.gradient_method=Nonereturntransform_program,config
[docs]defsimulate(self,circuit:qml.tape.QuantumScript,seed=None,debugger=None,)->Result:"""Simulate a single quantum script. Args: circuit (QuantumTape): The single circuit to simulate debugger (_Debugger): The debugger to use Returns: tuple(TensorLike): The results of the simulation This function assumes that all operations are Clifford. >>> qs = qml.tape.QuantumScript([qml.Hadamard(wires=0)], [qml.expval(qml.Z(0)), qml.state()]) >>> qml.devices.DefaultClifford().simulate(qs) (array(0), array([[0, 1, 0], [1, 0, 0]])) """# Account for custom labelled wirescircuit=circuit.map_to_standard_wires()# Build a stim circuit, tableau and simulatorstim_circuit=stim.Circuit()tableau_simulator=stim.TableauSimulator()ifself.wiresisnotNone:tableau_simulator.set_num_qubits(len(self.wires))# Account for state preparation operationprep=Noneiflen(circuit)>0andisinstance(circuit[0],qml.operation.StatePrepBase):prep=circuit[0]use_prep_ops=bool(prep)# TODO: Add a method to prepare directly from a Tableauifuse_prep_ops:stim_tableau=stim.Tableau.from_state_vector(qml.math.reshape(prep.state_vector(wire_order=list(circuit.op_wires)),(1,-1))[0],endian="big",)stim_circuit+=stim_tableau.to_circuit()# Iterate over the gates --> manage them manually or apply them to circuitglobal_phase_ops=[]foropincircuit.operations[use_prep_ops:]:gate,wires=_pl_op_to_stim(op)ifgateisnotNone:# Note: This is a lot faster than doing `stim_ct.append(gate, wires)`stim_circuit.append_from_stim_program_text(f"{gate}{wires}")else:ifisinstance(op,qml.GlobalPhase):global_phase_ops.append(op)ifisinstance(op,qml.Snapshot):self._apply_snapshot(circuit,stim_circuit,op,global_phase_ops,debugger)tableau_simulator.do_circuit(stim_circuit)global_phase=qml.GlobalPhase(qml.math.sum(op.data[0]foropinglobal_phase_ops))# Perform measurements based on whether shots are providedifcircuit.shots:meas_results=self.measure_statistical(circuit,stim_circuit,seed=seed)else:meas_results=self.measure_analytical(circuit,stim_circuit,tableau_simulator,global_phase)returnmeas_results[0]iflen(meas_results)==1elsetuple(meas_results)
@propertydef_analytical_measurement_map(self):"""Maps measurement type to the desired analytic measurement method."""return{DensityMatrixMP:self._measure_density_matrix,StateMP:self._measure_state,# kwargs -> circuit, global_phaseExpectationMP:self._measure_expectation,# kwargs -> circuit, stim_circuitVarianceMP:self._measure_variance,VnEntropyMP:self._measure_vn_entropy,# kwargs -> circuitMutualInfoMP:self._measure_mutual_info,# kwargs -> circuitPurityMP:self._measure_purity,# kwargs -> circuitProbabilityMP:self._measure_probability,# kwargs -> circuit, stim_circuit}@propertydef_statistical_measurement_map(self):"""Maps measurement type to the desired sampling measurement method."""return{ExpectationMP:self._sample_expectation,VarianceMP:self._sample_variance,ClassicalShadowMP:self._sample_classical_shadow,ShadowExpvalMP:self._sample_expval_shadow,}def_apply_snapshot(self,circuit:qml.tape.QuantumScript,stim_circuit,operation:qml.Snapshot,global_phase_ops:Sequence[qml.GlobalPhase],debugger=None,):"""Apply a snapshot operation to the stim circuit."""ifdebuggerisnotNoneanddebugger.active:meas=operation.hyperparameters["measurement"]orStateMP()measurement_func=self._analytical_measurement_map.get(type(meas),None)ifmeasurement_funcisNone:# pragma: no coverraiseqml.DeviceError(f"Snapshots of {type(meas)} are not yet supported on default.clifford.")# Build a temporary simulator for obtaining statesnap_sim=stim.TableauSimulator()ifself.wiresisnotNone:snap_sim.set_num_qubits(len(self.wires))snap_sim.do_circuit(stim_circuit)global_phase=qml.GlobalPhase(qml.math.sum(op.data[0]foropinglobal_phase_ops))snap_result=measurement_func(meas,snap_sim,circuit=circuit,stim_circuit=stim_circuit,global_phase=global_phase,)# Add to the debugger snapshotdebugger.snapshots[operation.tagorlen(debugger.snapshots)]=snap_result@staticmethoddef_measure_observable_sample(meas_obs,stim_circuit,shots,sample_seed):"""Compute sample output from a stim circuit for a given Pauli observable"""meas_dict={"X":"MX","Y":"MY","Z":"MZ","_":"M"}ifisinstance(meas_obs,BasisStateProjector):stim_circ=stim_circuit.copy()stim_circ.append_from_stim_program_text("M "+" ".join(map(str,meas_obs.wires)))sampler=stim_circ.compile_sampler(seed=sample_seed)return[qml.math.array(sampler.sample(shots=shots),dtype=int)],qml.math.array([1.0])coeffs,paulis=_pl_obs_to_linear_comb(meas_obs)samples=[]forpauli,wireinpaulis:stim_circ=stim_circuit.copy()forop,wrinzip(pauli,wire):ifop!="I":stim_circ.append(meas_dict[op],wr)sampler=stim_circ.compile_sampler(seed=sample_seed)samples.append(qml.math.array(sampler.sample(shots=shots),dtype=int))returnsamples,qml.math.array(coeffs)# pylint:disable=protected-access
[docs]defmeasure_statistical(self,circuit,stim_circuit,seed=None):"""Given a circuit, compute samples and return the statistical measurement results."""# Compute samples via circuits from tableaunum_shots=circuit.shots.total_shotssample_seed=seedifisinstance(seed,int)elseself._rng.integers(2**31-1,size=1)[0]results=[]formeasincircuit.measurements:measurement_func=self._statistical_measurement_map.get(type(meas),None)ifmeasurement_funcisnotNone:res=measurement_func(meas,stim_circuit,shots=num_shots,seed=sample_seed)else:# Decide wire ordermeas_wires=meas.wiresifmeas.wireselserange(stim_circuit.num_qubits)wire_order={wire:idxforidx,wireinenumerate(meas.wires)}# Decide measurement opmeas_op=meas.obsorqml.prod(*[qml.Z(idx)foridxinmeas_wires])samples=self._measure_observable_sample(meas_op,stim_circuit,num_shots,sample_seed)[0]# Check if the rotation was permissibleiflen(samples)>1:raiseqml.QuantumFunctionError(f"Observable {meas_op.name} is not supported for rotating probabilities on {self.name}.")# Process the result from samplesres=meas.process_samples(samples=np.array(samples),wire_order=wire_order)# Post-processing for special casesifisinstance(meas,CountsMP):res=res[0]elifisinstance(meas,SampleMP):res=np.squeeze(res)results.append(res)returnresults
[docs]defmeasure_analytical(self,circuit,stim_circuit,tableau_simulator,global_phase):"""Given a circuit, compute tableau and return the analytical measurement results."""results=[]formeasincircuit.measurements:# signature: measurement_func(meas, tableaus_simulator, **kwargs) -> meas_resultmeasurement_func=self._analytical_measurement_map.get(type(meas),None)ifmeasurement_funcisNone:# pragma: no coverraiseNotImplementedError(f"default.clifford doesn't support the {type(meas)} measurement analytically at the moment.")results.append(measurement_func(meas,tableau_simulator,circuit=circuit,stim_circuit=stim_circuit,global_phase=global_phase,))returnresults
@staticmethoddef_measure_density_matrix(meas,tableau_simulator,**_):"""Measure the density matrix from the state of simulator device."""wires=list(meas.wires)state_vector=qml.math.array(tableau_simulator.state_vector(endian="big"))returnqml.math.reduce_dm(qml.math.einsum("i, j->ij",state_vector,state_vector),wires)def_measure_state(self,_,tableau_simulator,**kwargs):"""Measure the state of the simualtor device."""wires=kwargs.get("circuit").wiresglobal_phase=kwargs.get("global_phase",qml.GlobalPhase(0.0))ifself._tableau:# Stack according to Sec. III, arXiv:0406196 (2008)tableau=tableau_simulator.current_inverse_tableau().inverse()x2x,x2z,z2x,z2z,x_signs,z_signs=tableau.to_numpy()pl_tableau=np.vstack((np.hstack((x2x,x2z,x_signs.reshape(-1,1))),np.hstack((z2x,z2z,z_signs.reshape(-1,1))),)).astype(int)ifpl_tableau.shape==(0,1)andlen(wires):returnnp.array([[1.0,0.0,0.0],[0.0,1.0,0.0]])returnpl_tableaustate=qml.math.array(tableau_simulator.state_vector(endian="big"))ifstate.shape==(1,)andlen(wires):# following is faster than using np.eye(length=1, size, index)state=qml.math.zeros(2**len(wires),dtype=complex)state[0]=1.0+0.0jreturnstate*qml.matrix(global_phase)[0][0]def_measure_expectation(self,meas,tableau_simulator,**kwargs):"""Measure the expectation value with respect to the state of simulator device."""meas_obs=meas.obsifisinstance(meas_obs,BasisStateProjector):kwargs["prob_states"]=qml.math.array([meas_obs.data[0]])returnself._measure_probability(qml.probs(wires=meas_obs.wires),tableau_simulator,**kwargs).squeeze()# Get the observable for the expectation value measurementcoeffs,paulis=_pl_obs_to_linear_comb(meas_obs)expecs=qml.math.zeros_like(coeffs)foridx,(pauli,wire)inenumerate(paulis):pauli_term=["I"]*max(np.max(list(wire))+1,tableau_simulator.num_qubits)forop,wrinzip(pauli,wire):pauli_term[wr]=opstim_pauli=stim.PauliString("".join(pauli_term))expecs[idx]=tableau_simulator.peek_observable_expectation(stim_pauli)returnqml.math.dot(coeffs,expecs)def_measure_variance(self,meas,tableau_simulator,**_):"""Measure the variance with respect to the state of simulator device."""meas_obs1=meas.obs.simplify()meas_obs2=(meas_obs1**2).simplify()# use the naive formula for variance, i.e., Var(Q) = ⟨𝑄^2⟩−⟨𝑄⟩^2return(self._measure_expectation(ExpectationMP(meas_obs2),tableau_simulator)-self._measure_expectation(ExpectationMP(meas_obs1),tableau_simulator)**2)def_measure_vn_entropy(self,meas,tableau_simulator,**kwargs):"""Measure the Von Neumann entropy with respect to the state of simulator device."""wires=kwargs.get("circuit").wirestableau=tableau_simulator.current_inverse_tableau().inverse()z_stabs=qml.math.array([tableau.z_output(wire)forwireinrange(len(wires))])returnself._measure_stabilizer_entropy(z_stabs,list(meas.wires),meas.log_base)def_measure_mutual_info(self,meas,tableau_simulator,**kwargs):"""Measure the mutual information between the subsystems of simulator device."""wires=kwargs.get("circuit").wirestableau=tableau_simulator.current_inverse_tableau().inverse()z_stabs=qml.math.array([tableau.z_output(wire)forwireinrange(len(wires))])indices0,indices1=getattr(meas,"_wires")returnself._measure_stabilizer_entropy(z_stabs,list(indices0),meas.log_base)+self._measure_stabilizer_entropy(z_stabs,list(indices1),meas.log_base)def_measure_purity(self,meas,tableau_simulator,**kwargs):r"""Measure the purity of the state of simulator device. Computes the state purity using the monotonically decreasing second-order Rényi entropy form given in `Sci Rep 13, 4601 (2023) <https://www.nature.com/articles/s41598-023-31273-9>`_. We utilize the fact that Rényi entropies are equal for all Rényi indices ``n`` for the stabilizer states. Args: stabilizer (TensorLike): stabilizer set for the system wires (Iterable): wires describing the subsystem log_base (int): base for the logarithm. Returns: (float): entanglement entropy of the subsystem """wires=kwargs.get("circuit").wiresifwires==meas.wires:returnqml.math.array(1.0)tableau=tableau_simulator.current_inverse_tableau().inverse()z_stabs=qml.math.array([tableau.z_output(wire)forwireinrange(len(wires))])return2**(-self._measure_stabilizer_entropy(z_stabs,list(meas.wires),log_base=2))# pylint: disable=protected-access@staticmethoddef_measure_stabilizer_entropy(stabilizer,wires,log_base=None):r"""Computes the Rényi entanglement entropy using stabilizer information. Computes the Rényi entanglement entropy :math:`S_A` for a subsytem described by :math:`A`, :math:`S_A = \text{rank}(\text{proj}_A {\mathcal{S}}) - |A|`, where :math:`\mathcal{S}` is the stabilizer group for the system using the theory described in Appendix A.1.d of `arXiv:1901.08092 <https://arxiv.org/abs/1901.08092>`_. Args: stabilizer (TensorLike): stabilizer set for the system wires (Iterable): wires describing the subsystem log_base (int): base for the logarithm. Returns: (float): entanglement entropy of the subsystem """# Get the number of qubits for the systemnum_qubits=qml.math.shape(stabilizer)[0]# Von Neumann entropy of a stabilizer state is zeroiflen(wires)==num_qubits:return0.0# Build a binary matrix desribing the stabilizers using the Pauli wordspauli_dict={0:"I",1:"X",2:"Y",3:"Z"}terms=[qml.pauli.PauliWord({idx:pauli_dict[ele]foridx,eleinenumerate(row)})forrowinstabilizer]binary_mat=qml.pauli.utils._binary_matrix_from_pws(terms,num_qubits)# Partition the binary matrix to represent the subsystempartition_mat=qml.math.hstack((binary_mat[:,num_qubits:][:,wires],binary_mat[:,:num_qubits][:,wires],))# Use the reduced row echelon form for finding rank efficiently# tapering always come in handy :)rank=qml.math.sum(qml.math.any(qml.qchem.tapering._reduced_row_echelon(partition_mat),axis=1))# Compute the entropyentropy=qml.math.log(2)*(rank-len(wires))# Determine wether to use any log baseiflog_baseisNone:returnentropyreturnentropy/qml.math.log(log_base)# pylint: disable=too-many-branches, too-many-statementsdef_measure_probability(self,meas,_,**kwargs):r"""Measure the probability of each computational basis state. Computes the probability for each of the computational basis state vector iteratively according to the follow pseudocode. 1. First, a complete basis set is built based on measured wires' length `l` by transforming integers :math:`[0, 2^l)` to their corresponding binary vector form, if the selective target states for computing probabilities have not been specified in the ``kwargs``. 2. Second, We then build a `stim.TableauSimulator` based on the input circuit. If an observable `obs` is given, an additional diagonalizing circuit is appended to the input circuit for rotating the computational basis based on the `diagonalizing_gates` method of the observable. 3. Finally, for every basis state, we iterate over each measured qubit `q_i` and peek if it can be collapsed to the state :math:`|0\rangle` / :math:`\rangle`1` corresponding to the bit `0` / `1` in the basis state vector at `i`th index. 4. If the qubit can be collapsed to the correct state, we do the post-selection and continue. If not, we identify it as an `unattainable` state and assign them with a zero probability. We do so for all the other basis states with the same `i`th index and keep this information stored in a visit-array. 5. Alternatively, if the qubit is in a superposition state, then it can collapse to either of the states :math:`|0\rangle` / :math:`\rangle`1`. We identify this as a `branching` scenario. We half the current probability and post-select based on the `i`th index of the basis state we are iterating. """circuit=kwargs.get("circuit")# Set the target statestgt_states=kwargs.get("prob_states",None)# Obtain the measurement wires for getting the basis statesmobs_wires=meas.obs.wiresifmeas.obselsemeas.wiresmeas_wires=mobs_wiresifmobs_wireselsecircuit.wires# Build the complete computational basis,# this will be expensive for larger circuits (> 24 qubits onwards)iftgt_statesisNone:num_wires=len(meas_wires)basis_vec=np.arange(2**num_wires)[:,np.newaxis]tgt_states=(((basis_vec&(1<<np.arange(num_wires)[::-1])))>0).astype(int)# Rotate the circuit basis to computational basisdiagonalizing_cit=kwargs.get("stim_circuit").copy()diagonalizing_ops=[]ifnotmeas.obselsemeas.obs.diagonalizing_gates()fordiag_opindiagonalizing_ops:# Check if it is Cliffordifdiag_op.namenotin_OPERATIONS_MAP:# pragma: no coverraiseValueError(f"Currently, we only support observables whose diagonalizing gates are Clifford, got {diag_op}")# Add to the circuit to rotate the basisstim_op=_pl_op_to_stim(diag_op)ifstim_op[0]isnotNone:diagonalizing_cit.append_from_stim_program_text(f"{stim_op[0]}{stim_op[1]}")# Build the Tableau simulator from the diagonalized circuitcircuit_simulator=stim.TableauSimulator()circuit_simulator.do_circuit(diagonalizing_cit)ifnotself._tableau:state=self._measure_state(meas,circuit_simulator,circuit=circuit)returnmeas.process_state(state,wire_order=circuit.wires)iflen(meas_wires)>=tgt_states.shape[1]:meas_wires=meas_wires[:tgt_states.shape[1]]else:# pragma: no covercgc_states=[]forstateintgt_states:iflist(state[meas_wires])notincgc_states:cgc_states.append(list(state[meas_wires]))tgt_states=np.array(cgc_states)# Maintain a representaiton of basis states to build a visit-arraytgt_integs=np.array([int("".join(map(str,tgt_state)),2)fortgt_stateintgt_states])# Iterate over the required basis states and for each of them compute the probability# If an impossible branch occur keep a note of it via a visit-array# Worst case complexity O(B * M * N), where N is #measured_qubits and B is #basis_states.visited_probs=[]prob_res=np.ones(tgt_states.shape[0])fortgt_index,(tgt_integ,tgt_state)inenumerate(zip(tgt_integs,tgt_states)):iftgt_integinvisited_probs:continueprob_sim=circuit_simulator.copy()foridx,wireinenumerate(meas_wires):expectation=prob_sim.peek_z(wire)# (Eig --> Res) | -1 --> 1 | 1 --> 0 | 0 --> 0 / 1 |outcome=int(0.5*(1-expectation))ifnotexpectation:prob_res[tgt_index]/=2.0else:nope_idx=(np.where(np.squeeze(np.all(tgt_states[:,:idx]==tgt_state[:idx],axis=-1))&tgt_states[:,idx]!=outcome)[0]ifidxelsenp.where(tgt_states[:,idx]!=outcome)[0])nope_idx=np.setdiff1d(nope_idx,visited_probs)prob_res[nope_idx]=0.0visited_probs.extend(tgt_integs[nope_idx])iftgt_state[idx]!=outcome:prob_res[tgt_index]=0.0breakprob_sim.postselect_z(wire,desired_value=tgt_state[idx])visited_probs.append(tgt_integ)returnprob_resdef_sample_expectation(self,meas,stim_circuit,shots,seed):"""Measure the expectation value with respect to samples from simulator device."""# Get the observable for the expectation value measurementmeas_op=meas.obssamples,coeffs=self._measure_observable_sample(meas_op,stim_circuit,shots,seed)ifisinstance(meas_op,BasisStateProjector):matches=np.where((samples[0]==meas_op.data[0]).all(axis=1))[0]returnlen(matches)/shotsexpecs=[qml.math.mean(qml.math.power([-1]*shots,qml.math.sum(sample,axis=1)))forsampleinsamples]returnqml.math.dot(coeffs,expecs)def_sample_variance(self,meas,stim_circuit,shots,seed):"""Measure the variance with respect to samples from simulator device."""# Get the observable for the expectation value measurementmeas_obs1=meas.obs.simplify()meas_obs2=(meas_obs1**2).simplify()# use the naive formula for variance, i.e., Var(Q) = ⟨𝑄^2⟩−⟨𝑄⟩^2return(self._sample_expectation(qml.expval(meas_obs2),stim_circuit,shots,seed)-self._sample_expectation(qml.expval(meas_obs1),stim_circuit,shots,seed)**2)@staticmethoddef_measure_single_sample(stim_ct,meas_ops,meas_idx,meas_wire):"""Sample a single qubit Pauli measurement from a stim circuit"""stim_sm=stim.TableauSimulator()stim_sm.do_circuit(stim_ct)res=[0]*meas_idx+meas_ops+[0]*(meas_wire-meas_idx-1)res=[int(r)forrinres]returnstim_sm.measure_observable(stim.PauliString(res))def_sample_classical_shadow(self,meas,stim_circuit,shots,seed):"""Measures classical shadows from the state of simulator device"""meas_seed=meas.seedorseedmeas_wire=stim_circuit.num_qubitsbits=[]recipes=np.random.RandomState(meas_seed).randint(3,size=(shots,meas_wire))# Random Pauli basis to be used for measurementsforrecipeinrecipes:bits.append([self._measure_single_sample(stim_circuit,[int(rec)+1],idx,meas_wire)foridx,recinenumerate(recipe)])returnnp.asarray(bits,dtype=int),np.asarray(recipes,dtype=int)def_sample_expval_shadow(self,meas,stim_circuit,shots,seed):"""Measures expectation value of a Pauli observable using classical shadows from the state of simulator device."""bits,recipes=self._sample_classical_shadow(meas,stim_circuit,shots,seed)# TODO: Benchmark scaling for larger number of circuits for this existing functionalitywires_map=list(range(stim_circuit.num_qubits))shadow=qml.shadows.ClassicalShadow(bits,recipes,wire_map=wires_map)returnshadow.expval(meas.H,meas.k)