From 431ccbcb2f1cd650d17080cc3d0956c90aba9050 Mon Sep 17 00:00:00 2001 From: Everett Date: Wed, 7 May 2025 13:18:47 -0700 Subject: [PATCH 1/5] unify circuit architecture and add measurement support MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## Major Update Summary - Unified circuit architecture to support both unitary and measurement operations, and more flexible circuit class to compose any type of operation. - Unified measurement and post-selection framework to calculate state projection and measurement probability, using log2 scale to handle all probabilities. - Matrix representation of states and operators moves away form QuTiP dependency in favor of pure NumPy implementations. - Caching mechanisms for both measurement outcomes and random unitaries to support classical shadow collection via systematic forward/backward passes through the circuit. ## File-by-File Updates ### `__init__.py` - All imports now reference new class and function names - new state construction function: `bit_state`, - new circuit classes: `Measurement`, `Layer`, `Circuit` (removing `CliffordLayer`, `MeasurementLayer`, `CliffordCircuit`) - new gate and measurements: `SWAP`, `CZ`, `CX`, `measurement_layer` - Remove debug `print` statements ### `utils.py` - Unify projection trace and post-selection for stabilizer state. - Introduce `stabilizer_postselect` to post-select a stabilizer state based on a projection operator (given as another stabilizer state), and returns the base-2 logarithm `log2prob` of the likelihood for the post-selection to succeed. - All measurement and post-selection likelihood are given as log2 likelihood to avoid the underflow of the exponentially small probability for large system. - Remove `stabilizer_projection_trace` that computes the trace of a projection operator with a stabilizer state density matrix. Functionality achieved by `stabilizer_postselect`, from which the trace can be computed as `2**log2prob`. - Remove `stabilizer_postselection`, functionality replaced by `stabilizer_postselect`. - Remove `decompose`, which decomposes a Pauli operator to product of stabilizers and destabilizers with phase factor. This functionality has been achieved by `obj.transform_by(state.to_map().inverse())` where `state` is the `StabilizerState` that hosts the stabilizers and destabilizers, and `obj` can be any Pauli operator (or operator list) to be decomposed. The idea is to use the state-map correspondence to construct the state encoder, then use the inverse map to construct the state decoder, and use the decoder to transform any operator back to the logical/syndrome basis. - Add comments to explain the unified approach to realize `stabilizer_measure`, `stabilizer_postselect` and `stabilizer_project`. ### `paulialg.py` - Remove `import qutip as qt` , and all `to_qutip` functions in all classes. Replace them by `to_numpy` functions, that returns matrix representation of Pauli algebra objects as numpy array, such that installation of `qutip` package will no longer be required. - The `to_numpy` method for PauliList is now vectorized for efficiency, avoiding explicit Python loops. - `PauliList` and `Pauli` constructors now accept `**kwargs` to absorb extra arguments from subclasses. ### `stabilizer.py` - Bug fixed: The `StabilizerState.__init__(gs, r=0, **kwargs)` signature missed the `ps` argument. This results in the `.copy` method fails to make faithful copy for the stabilizer state (the phase factor will always be ignored under copy). The new signature `StabilizerState.__init__(*args, **kwargs)` handles the `.r` attribute initialization by `self.r=kwargs.pop(r,0)`. - `StabilizerState` now only accepts r as a keyword argument, not as a positional argument. - The `set_r` method has been removed; as `r` is set directly via the constructor or attribute assignment. - All code that previously used `.set_r(r)` now uses `r=r` in the constructor or direct assignment. - Default value for `r` is set to be 0 (indicating pure state) in all relevant classes and functions, instead of `None` in functions like `random_pauli_state` or `random_clifford_state`. - Error handling: Custom `Error` exceptions have been replaced with standard exceptions like `ValueError`. - The logic for stabilizer state expectation value and post-selection is more concise and robust, by using the `stabilizer_postselect` function uniformly for both trace evaluation and state projection. - The `get_prob` method is now achieved by post-selecting the target bit-string state `bit_state` (computational basis state) and return the successful probability of the post-selection. - Update the `.expect` method: - Stabilizer state expectation value (trace with a projection operator) is now realized by `state_postselect`. - Implement the fugacity setting to help re-weighting the Pauli observables expectation value for classical shadow purpos. ### `circuit.py` - Restructure all circuit classes (`CliffordGate`, `Measurement`, `Layer`, `Circuit`) to enable measurements. - All classes now has the `.unitary` property to indicate if the object is purely unitary or has measurement. - The `.take()` method has been replaced by `.append()` for adding operations in all container classes (`Layer`, `Circuit`). This is to be consistent with the standard notation in other popular packages like Qiskit or Cirq. - All `forward` and `backward` methods now always return the tuple (`obj`, `log2prob`), even for unitary processes (where `log2prob=0.0`). This unifies the interface and removes the need for branching on unitarity. - New class `Measurement` is introduce at the same level of `CliffordGate` to perform computational basis measurement on designated qubits, such that mid-circuit measurement can be implemented. - `CliffordLayer` and `MeasurementLayer` are unified into `Layer`, which can be either unitary or non-unitary, providing the ability to handle different type of operations together in a unified manner. - The functionality of `MeasurementLayer` is replaced by a new function `measurement_layer(N)`, which returns a `Circuit` object that contains a single layer of measurement of all the N qubits in computational basis simultaneously, and can be used by appending to other circuit object. - `Circuit` no longer takes qubit number `N` as input parameter. The `.N` attribute is also removed from the class, allowing circuit (as a collection of operations) to be defined without specifying the system size first. - Instead, `Circuit` class constructor now takes circuit objects (`CliffordGate`, `Measurement`, `Layer`, `Circuit`) directly to build the circuit. Objects passing to the constructor will be appended to the circuit in chronological order. - All references to `CliffordLayer`, `CliffordCircuit` are removed in favor of `Layer`, `Circuit`. - Measurement outcome caching - Measurement outcome will not be returned as output, but cached by the `Measurement` object during the forward pass, which is then reused for post-selection in the backward pass. This enables classical shadows to be collected by a forward and backward pass through the entire circuit systematically. - Measurement outcomes are now stored only at the `Measurement` object level, with `Layer` and `Circuit` aggregating via properties. - The measurement outcome property `out` is now a dynamic property at all levels, with setter and deleter for hierarchical management. - Introduce `.reset()` method to reset the cached measurement outcome in `Measurement` object. - Random unitary caching - Random Clifford gate will sample its unitary map in the forward pass and cache it in `CliffordGate` object to be used by the backward pass (implementing the inverse unitary map). This enables classical shadows to be collected by a forward and backward pass through the entire circuit systematically. - Introduce `.reset()` method for all classes to reset the cached Clifford maps in Clifford gates. --- pyclifford/__init__.py | 18 +- pyclifford/circuit.py | 976 ++++++++++++++++++++------------------- pyclifford/device.py | 10 +- pyclifford/paulialg.py | 155 ++++--- pyclifford/stabilizer.py | 196 +++++--- pyclifford/utils.py | 413 ++++++++--------- 6 files changed, 933 insertions(+), 835 deletions(-) diff --git a/pyclifford/__init__.py b/pyclifford/__init__.py index 2e76eee..ab83bbe 100644 --- a/pyclifford/__init__.py +++ b/pyclifford/__init__.py @@ -1,20 +1,14 @@ -# import paulialg -# import stabilizer -# import circuit -# import qchem from .paulialg import ( Pauli,PauliList,PauliMonomial,PauliPolynomial, pauli, paulis, pauli_identity, pauli_zero) from .stabilizer import( CliffordMap,StabilizerState, identity_map, random_pauli_map, random_clifford_map, clifford_rotation_map, - stabilizer_state, maximally_mixed_state, zero_state, one_state, ghz_state, - random_pauli_state, random_clifford_state,random_bit_state) + stabilizer_state, maximally_mixed_state, zero_state, one_state, bit_state, + ghz_state, random_pauli_state, random_clifford_state,random_bit_state) from .circuit import( - CNOT,C,X,Y,Z,H,S,CliffordGate,CliffordLayer,MeasureLayer,Circuit, - clifford_rotation_gate, - identity_circuit, brickwall_rcc, onsite_rcc, global_rcc, + CliffordGate,Measurement,Layer,Circuit, + CNOT,SWAP,CZ,CX,C,X,Y,Z,H,S,clifford_rotation_gate, + identity_circuit, brickwall_rcc, onsite_rcc, global_rcc, measurement_layer, diagonalize, SBRG) -from .device import ClassicalShadow -#from .qchem import qchem_hamiltonian - +from .device import ClassicalShadow \ No newline at end of file diff --git a/pyclifford/circuit.py b/pyclifford/circuit.py index 1626ba0..226b15a 100644 --- a/pyclifford/circuit.py +++ b/pyclifford/circuit.py @@ -1,12 +1,12 @@ import numpy -import numpy as np -from .utils import mask, condense, pauli_diagonalize1,stabilizer_measure -from .paulialg import Pauli, pauli, PauliMonomial, pauli_zero -from .stabilizer import (StabilizerState, CliffordMap, - zero_state, identity_map, clifford_rotation_map, random_clifford_map) +import warnings +from .utils import mask, condense, pauli_diagonalize1 +from .paulialg import Pauli, pauli, paulis, PauliMonomial, pauli_zero +from .stabilizer import (StabilizerState, CliffordMap, identity_map, + clifford_rotation_map, random_clifford_map) class CliffordGate(object): - '''Represents a Clifford gate. + '''Represents a Clifford unitary gate. Parameters: *qubits: int - the qubits that this gate acts on. @@ -20,6 +20,7 @@ class CliffordGate(object): described by the Clifford map, which is a table specifying how each single Pauli operator gets mapped to. (forward and backward maps must be inverse to each other). + unitary: bool - indicates that gate is unitary. Note: if either the geneator or Clifford maps are specified, the gate will represent the specific unitary transformation; otherwise, the gate @@ -27,10 +28,15 @@ class CliffordGate(object): def __init__(self, *qubits): self.qubits = qubits # the qubits this gate acts on self.n = len(self.qubits) # number of qubits it acts on + self.reset() + self.unitary = True # gate is unitary + + def reset(self): + # reset the gate to a random Clifford gate self.generator = None self.forward_map = None self.backward_map = None - + def __repr__(self): return '[{}]'.format(','.join(str(qubit) for qubit in self.qubits)) @@ -38,30 +44,40 @@ def set_generator(self, gen): if not isinstance(gen, Pauli): raise TypeError("Rotation generator must be a Pauli string") self.generator = gen + def set_forward_map(self,forward_map): if not isinstance(forward_map, CliffordMap): raise TypeError("Forward map must be a instance of CliffordMap") self.forward_map = forward_map + def set_backward_map(self,backward_map): if not isinstance(backward_map, CliffordMap): raise TypeError("Backward map must be a instance of CliffordMap") self.backward_map = backward_map - def copy(self): - gate = CliffordGate(*self.qubits) + new = CliffordGate(*self.qubits) if self.generator is not None: - gate.generator = self.generator.copy() + new.generator = self.generator.copy() if self.forward_map is not None: - gate.forward_map = self.forward_map.copy() + new.forward_map = self.forward_map.copy() if self.backward_map is not None: - gate.backward_map = self.backward_map.copy() - return gate + new.backward_map = self.backward_map.copy() + return new - def independent_from(self, other_gate): - return len(set(self.qubits) & set(other_gate.qubits))==0 + def independent_from(self, other): + return len(set(self.qubits) & set(other.qubits))==0 def forward(self, obj): + '''Apply the gate forward in time. (inplace update) + Forward transformation: O -> U O U^\dagger + + Input: + obj: Pauli, PauliList, StabilizerState - the object to be transformed + + Output: + obj: (same as input type) - the object after the gate is applied + log2prob: real - always 0.0 for unitary transformation''' if self.generator is not None: # if generator is given, use generator if self.n == obj.N: # global gate obj.rotate_by(self.generator) @@ -72,6 +88,7 @@ def forward(self, obj): if self.backward_map is None: # if both maps not given, treated as random gate clifford_map = random_clifford_map(self.n) + self.forward_map = clifford_map # record as forward map else: self.forward_map = self.backward_map.inverse() clifford_map = self.forward_map @@ -80,12 +97,20 @@ def forward(self, obj): if self.n == obj.N: # global gate obj.transform_by(clifford_map) else: # local gate - # print("Clifford gate acting: ",self.qubits) - # print("mask in CliffordGate: ",mask(self.qubits, obj.N)) obj.transform_by(clifford_map, mask(self.qubits, obj.N)) - return obj + log2prob = 0.0 # unitary gate is deterministic, log2prob is always 0.0 + return obj, log2prob def backward(self, obj): + '''Apply the gate backward in time. (inplace update) + Backward transformation: O -> U^\dagger O U + + Input: + obj: Pauli, PauliList, StabilizerState - the object to be transformed + + Output: + obj: (same as input type) - the object after the gate is applied + log2prob: real - always 0.0 for unitary transformation''' if self.generator is not None: # if generator is given, use generator if self.n == obj.N: # global gate obj.rotate_by(-self.generator) @@ -96,19 +121,21 @@ def backward(self, obj): if self.forward_map is None: # if both maps not given, treated as random gate clifford_map = random_clifford_map(self.n) + self.backward_map = clifford_map # record as backward map else: self.backward_map = self.forward_map.inverse() clifford_map = self.backward_map else: clifford_map = self.backward_map - if False and self.n == obj.N: # global gate + if self.n == obj.N: # global gate obj.transform_by(clifford_map) else: # local gate obj.transform_by(clifford_map, mask(self.qubits, obj.N)) - return obj + log2prob = 0.0 # unitary gate is deterministic, log2prob is always 0.0 + return obj, log2prob def compile(self): - '''construct forward and backward Clifford maps for this gate''' + '''Construct forward and backward Clifford maps for this gate''' if self.generator is not None: self.forward_map = clifford_rotation_map(self.generator) self.backward_map = clifford_rotation_map(-self.generator) @@ -123,286 +150,265 @@ def compile(self): self.backward_map = self.forward_map.inverse() return self +class Measurement(object): + '''Represents a computational basis measurement. + + Parameters: + *qubits: int - the qubits to be measured + + Data: + out: int (L) - array of measurement outcomes on corresponding qubits + None before measurement, populated after forward measurement. + unitary: bool - indicates that measurement is not unitary.''' + def __init__(self, *qubits): + self.qubits = qubits # the qubits to be measured + self.n = len(self.qubits) # number of qubits to be measured + self._out = None # container for measurement outcomes, will be populated after forward measurement + self.unitary = False # measurement is not unitary -class CliffordLayer(object): - '''Representes a layer of Clifford gates. + def reset(self): + self._out = None + + @property + def out(self): + if self._out is None: + raise ValueError("{} measurement outcome is not sampled yet.".format(repr(self))) + return self._out + + @out.setter + def out(self, out): + self._out = out - Parameters: - *gate: CliffordGate - the gates that this layer contains''' - def __init__(self, *gates): - self.gates = list(gates) # the gates this layer have - self.prev_layer = None # the previous layer - self.next_layer = None # the next layer - self.forward_map = None - self.backward_map = None - def __repr__(self): - return '|{}|'.format(''.join(repr(gate) for gate in self.gates)) + return '<{}>'.format(','.join(str(qubit) for qubit in self.qubits)) - def copy(self): - layer = CliffordLayer(*[gate.copy() for gate in self.gates]) - if self.forward_map is not None: - layer.forward_map = self.forward_map.copy() - if self.backward_map is not None: - layer.backward_map = self.backward_map.copy() - return layer - def independent_from(self, other_gate): - return all(gate.independent_from(other_gate) for gate in self.gates) + def copy(self): + new = Measurement(*self.qubits) + if self._out is not None: + new._out = self._out.copy() + return new - def take(self, gate): - if (self.prev_layer is None) or isinstance(self.prev_layer, MeasureLayer): - # if I have no previous layer - # or previous layer is measurement - self.gates.append(gate) # I will take the gate - else: # if I have a previous layer, check it - if self.prev_layer.independent_from(gate): # if independent (not overlapping) - self.prev_layer.take(gate) # previous layer take the gate - else: # if not independent - self.gates.append(gate) # I will have to keep the gate + def independent_from(self, other): + return len(set(self.qubits) & set(other.qubits))==0 def forward(self, obj): - if self.forward_map is None: - for gate in self.gates: - gate.forward(obj) - else: - obj.transform_by(self.forward_map) - return obj + '''Implements the measurement (non-deterministic sampling outcome). + (inplace update) + Forward transformation: rho -> M rho M^\dagger / Tr(M rho M^\dagger), + where M is sampled with probability Tr(M rho M^\dagger). + + Input: + obj: StabilizerState - the state to be measured + + Output: + obj: StabilizerState - the state after measurement + log2prob: real - log2 probability of the outcome''' + if not isinstance(obj, StabilizerState): # measurement only applicable to stabilizer state + raise NotImplementedError("the object {} is not a stabilizer state".format(repr(obj))) + # construct Z observables + obs = paulis(pauli({i: 'Z'}, obj.N) for i in self.qubits) + # perform measurement + self.out, log2prob = obj.measure(obs) + return obj, log2prob def backward(self, obj): - if self.backward_map is None: - for gate in self.gates: - gate.backward(obj) - else: - obj.transform_by(self.backward_map) - return obj - - def compile(self, N): - '''construct forward and backward Clifford maps for this layer''' - self.forward_map = identity_map(N) - self.backward_map = identity_map(N) - for gate in self.gates: - gate.compile() - self.forward_map.embed(gate.forward_map, mask(gate.qubits, N)) - self.backward_map.embed(gate.backward_map, mask(gate.qubits, N)) - return self - -class MeasureLayer(object): - ''' - Computational basis measurement - ''' - def __init__(self, *qubits,N): - self.qubits = qubits # qubits that are measured - self.result = None - self.log2prob = None - self.N = N - self.gs, self.ps = self.obs_gs_ps() - self.prev_layer = None # the previous layer - self.next_layer = None # the next layer - def __repr__(self): - return '|Mz[{}]|'.format(','.join(str(qubit) for qubit in self.qubits)) - def obs_gs_ps(self): - ps = np.zeros(len(self.qubits)).astype(int) - gs = np.zeros((len(self.qubits),2*self.N)).astype(int) - for i in range(len(self.qubits)): - gs[i,2*self.qubits[i]+1]=1 - return gs, ps - def forward(self,obj): - ''' - Measurement will in-place change the stabilizer state. - ''' - if not isinstance(obj,StabilizerState): - raise NotImplementedError("the object is not a stabilizer state") - tmp_gs_stb, tmp_ps_stb, tmp_r, tmp_out, tmp_log2prob = \ - stabilizer_measure(obj.gs,obj.ps,self.gs,self.ps,obj.r) - self.result = (-1)**tmp_out - self.log2prob = tmp_log2prob - return obj - def backward(self,obj,measure_result = None): - ''' - backward is post-selection + '''Postselect the measurement outcome (deterministic projection). + (inplace update) + Backward transformation: rho -> M^\dagger rho M / Tr(M^\dagger rho M) + where M is fixed by the measurement outcome generated in the forward pass. + Input: - - measure_result: list of integers, such as [1,-1,1] - ''' - if measure_result is not None: - if len(measure_result)!=len(self.qubits): - raise ValueError("classical bit does not match number of qubits.") - else: - for ii in range(1,len(measure_result)+1): - tmp = list(np.zeros(self.N).astype(int)) - tmp[self.qubits[-ii]]=3 - tmp_res = int((1-measure_result[-ii])/2) - prob = obj.postselect(pauli(tmp), tmp_res) - if prob == 0.0: - raise ValueError("Post-selection result is not possible, they are orthogonal states.") - return obj - else: - if self.result is None: - raise ValueError("post-selection result is not given.") - else: - for ii in range(1,len(self.result)+1): - tmp = list(np.zeros(self.N).astype(int)) - tmp[self.qubits[-ii]]=3 - tmp_res = int((1-self.result[-ii])/2) - prob = obj.postselect(pauli(tmp), tmp_res) - if prob == 0.0: - raise ValueError("Post-selection result is not possible, they are orthogonal states.") - return obj - -class CliffordCircuit(object): - '''Represents a circuit of Clifford gates. - - Examples: - # create a circuit - circ = CliffordCircuit() - # add a gate between qubits 0 and 1 - circ.gate(0,1) - # or take in a specific gate - g = pauli('-XX') - circ.take(clifford_rotation_gate(g))''' - def __init__(self,N): - self.N = N # number of qubits in the system - self.first_layer = CliffordLayer() - self.last_layer = self.first_layer + obj: StabilizerState - the state to be postselected + + Output: + obj: StabilizerState - the state after postselection + log2prob: real - log2 probability of successful postselection''' + if not isinstance(obj, StabilizerState): # postselection only applicable to stabilizer state + raise NotImplementedError("the object {} is not a stabilizer state".format(repr(obj))) + # construct Z observables + obs = paulis(pauli({i: 'Z'}, obj.N) for i in self.qubits) + log2prob = obj.postselect(obs, self.out) + if log2prob == -numpy.inf: # if postselection fails + warnings.warn("Impossible to postselect the recorded measurement outcomes on the current state in the backward pass. The resutlting state might be incorrect.") + return obj, log2prob + + +class Layer(object): + '''Representes a layer of Clifford gate or measurement operations. + + Parameters: + *ops: Operation - the operations (gates or measurements) in this layer + + Attributes: + prev_layer: Layer - the previous layer + next_layer: Layer - the next layer + forward_map: CliffordMap - the forward map of the layer (if unitary) + backward_map: CliffordMap - the backward map of the layer (if unitary) + unitary: bool - whether the layer is unitary''' + + def __init__(self, *ops): + self.ops = list(ops) # the operations in this layer + self.prev_layer = None # the previous layer + self.next_layer = None # the next layer self.forward_map = None self.backward_map = None - - def __repr__(self): - layout = '\n'.join(repr(layer) for layer in self.layers_backward()) - return 'CliffordCircuit(\n{})'.format(layout).replace('\n','\n ') + # the layer is unitary if all operations are unitary + self.unitary = all(op.unitary for op in self.ops) - def __getattr__(self, item): - if item == 'N': # if self.N not defined - # infer from gates (assuming last qubit is covered) - N = 0 - for layer in self.layers_forward(): - for gate in layer.gates: - N = max(N, max(gate.qubits)+1) - return N + def reset(self): + self.forward_map = None + self.backward_map = None + for op in self.ops: + op.reset() + + @property + def out(self): + outs = [op.out for op in self.ops if not op.unitary] + if outs: + return numpy.concatenate(outs) else: - return super().__getattribute__(item) + return numpy.array([], dtype=numpy.int_) + + @out.setter + def out(self, out): + if not self.unitary: + k = 0 + for op in self.ops: + if not op.unitary: + op.out = out[k:k+op.n] + k += op.n + def __repr__(self): + return '|{}|'.format(''.join(repr(op) for op in self.ops)) + + @property + def isempty(self): + return len(self.ops) == 0 + def copy(self): - circ = CliffordCircuit(self.N) - for i, layer in enumerate(self.layers_forward()): - new_layer = layer.copy() - if i == 0: - circ.first_layer = new_layer - circ.last_layer = new_layer - else: - circ.last_layer.next_layer = new_layer - new_layer.prev_layer = circ.last_layer - circ.last_layer = new_layer + new = Layer(*(op.copy() for op in self.ops)) if self.forward_map is not None: - circ.forward_map = self.forward_map.copy() + new.forward_map = self.forward_map.copy() if self.backward_map is not None: - circ.backward_map = self.backward_map.copy() - return circ - - def layers_backward(self): - # yield from last to first layers - layer = self.last_layer - while layer is not None: - yield layer - layer = layer.prev_layer + new.backward_map = self.backward_map.copy() + return new - def layers_forward(self): - # yield from first to last layers - layer = self.first_layer - while layer is not None: - yield layer - layer = layer.next_layer - - def take(self, gate): - if max(gate.qubits)>=self.N: - raise ValueError("The gate acting on unregistered qubits!") - if self.last_layer.independent_from(gate): # if last layer commute with the new gate - self.last_layer.take(gate) # the last layer takes the gate - else: # otherwise create a new layer to handle this - new_layer = CliffordLayer(gate) # a new layer with the new gate - # link to the layer structure - self.last_layer.next_layer = new_layer - new_layer.prev_layer = self.last_layer - self.last_layer = new_layer # new layer becomes the last - return self - - def gate(self, *qubits): - if max(qubits)>=self.N: - raise ValueError("The gate acting on unregistered qubits!") - return self.take(CliffordGate(*qubits)) # create a new gate - - def compose(self, other): - '''Compose the circuit with another circuit. - U = U_other U_self - - Parameters: - other: CliffordCircuit - another circuit to be combined. - - Note: composition will not update the compiled information. Need - compilation after circuit composition.''' - if not other.N == self.N: - raise ValueError("Qubit number does not match!") - for layer in other.layers_forward(): - for gate in layer.gates: - self.take(gate) + def independent_from(self, other): + return all(op.independent_from(other) for op in self.ops) + + def append(self, op): + '''append an operation into this layer.''' + if (self.prev_layer is not None) and self.prev_layer.independent_from(op): + # if the previous layer exists and is independent from the new operation + self.prev_layer.append(op) # append the operation to the previous layer + else: # otherwise, admit the new operation to the current layer + self.ops.append(op) + self.unitary = self.unitary and op.unitary + # Clear cached maps since layer had changed + self.forward_map = None + self.backward_map = None return self - + def forward(self, obj): - if self.forward_map is None: - for layer in self.layers_forward(): - layer.forward(obj) - else: - obj.transform_by(self.forward_map) - return obj - + '''Apply the layer forward. (inplace update)''' + log2prob = 0.0 + if self.unitary and self.forward_map is not None: + obj.transform_by(self.forward_map) + else: # otherwise, apply each operation forward (in parallel) + for op in self.ops: + obj, op_log2prob = op.forward(obj) + log2prob += op_log2prob + return obj, log2prob + def backward(self, obj): - if self.backward_map is None: - for layer in self.layers_backward(): - layer.backward(obj) - else: + '''Apply the layer backward. (inplace update)''' + log2prob = 0.0 + if self.unitary and self.backward_map is not None: + # if layer backward map has been compiled, use it obj.transform_by(self.backward_map) - return obj - - def compile(self, N=None): - '''Construct forward and backward Clifford maps for this circuit + else: # otherwise, apply each operation backward (in parallel) + for op in self.ops: + obj, op_log2prob = op.backward(obj) + log2prob += op_log2prob + return obj, log2prob + + def compile(self, N): + '''construct forward and backward Clifford maps for this layer - Note: The compilation creates a single Clifford map representing the - entire circuit, which allows it to run faster for deep circuits.''' - N = self.N if N is None else N + Input: + N: int - number of qubits in the system''' + if not self.unitary: + raise ValueError("only unitary layer can be compiled.") self.forward_map = identity_map(N) self.backward_map = identity_map(N) - for layer in self.layers_forward(): - layer.compile(N) - self.forward_map = self.forward_map.compose(layer.forward_map) - self.backward_map = self.backward_map.compose(layer.backward_map) + for op in self.ops: + op.compile() + self.forward_map.embed(op.forward_map, mask(op.qubits, N)) + self.backward_map.embed(op.backward_map, mask(op.qubits, N)) return self - - def povm(self, nsample): - '''Assuming computational basis measurement follows the circuit, this - will back evolve the computational basis state to generate prior POVM. - This returns a generator.''' - for _ in range(nsample): - zero = zero_state(self.N) - yield self.backward(zero) + class Circuit(object): + '''Represents a quantum process of Clifford circuit with gates and measurements. + + Attributes: + first_layer: Layer - the first layer of the circuit + last_layer: Layer - the last layer of the circuit + forward_map: CliffordMap - the forward map of the circuit (if unitary) + backward_map: CliffordMap - the backward map of the circuit (if unitary) + unitary: bool - whether the circuit is unitary + + Example: + >>> circ = Circuit() + >>> circ.gate(0) + >>> circ.measure(1) + >>> circ.gate(2) + >>> sigma, log2prob_fw = circ.forward(zero_state(N)) + >>> rho, log2prob_bw = circ.backward(sigma.copy()) + >>> print(sigma, rho, log2prob_fw, log2prob_bw) ''' - Reprsents a Clifford circuit with measurements - ''' - def __init__(self,N): - self.N = N # number of qubits in the system - self.first_layer = CliffordLayer() + def __init__(self, *objs): + self.first_layer = Layer() self.last_layer = self.first_layer - self.measure_result = [] - self.log2prob = 0.0 - self.unitary = True # if True, no measurement layer is inside + self.forward_map = None # forward map + self.backward_map = None # backward map + self.unitary = True # empty circuit is unitary + # append objects to the circuit + # each object can be CliffordGate, Measurement, Layer, or Circuit + for obj in objs: + self.append(obj) + + def reset(self): self.forward_map = None self.backward_map = None - self.num_of_measures = 0 + for layer in self.layers_forward(): + layer.reset() + + @property + def out(self): + outs = [layer.out for layer in self.layers_forward() if not layer.unitary] + if outs: + return numpy.concatenate(outs) + else: + return numpy.array([], dtype=numpy.int_) + + @out.setter + def out(self, out): + print('out:', out) # for debug + k = 0 + for layer in self.layers_forward(): + print(layer, layer.unitary) # for debug + if not layer.unitary: + n = sum(op.n for op in layer.ops if not op.unitary) + print('number measurement:', n) # for debug + print('out[k:k+n]:', out[k:k+n]) # for debug + layer.out = out[k:k+n] + k += n + def __repr__(self): layout = '\n'.join(repr(layer) for layer in self.layers_backward()) - c = 'CliffordCircuit(\n{})'.format(layout).replace('\n','\n ')+\ - '\n Unitary:{}'.format(self.unitary) + c = 'Circuit(\n{})'.format(layout).replace('\n','\n ') return c def layers_backward(self): @@ -411,142 +417,103 @@ def layers_backward(self): while layer is not None: yield layer layer = layer.prev_layer - + def layers_forward(self): # yield from first to last layers layer = self.first_layer while layer is not None: yield layer layer = layer.next_layer - - def take(self, gate): - ''' - it can take Clifford gate or measurement layers - ''' - if (max(gate.qubits)>=self.N): - raise ValueError("The gate acting on unregistered qubits!") - if isinstance(gate,CliffordGate): - if not isinstance(self.last_layer,MeasureLayer): - if self.last_layer.independent_from(gate): # if last layer commute with the new gate - self.last_layer.take(gate) # the last layer takes the gate - else: # otherwise create a new layer to handle this - new_layer = CliffordLayer(gate) # a new layer with the new gate - # link to the layer structure - self.last_layer.next_layer = new_layer - new_layer.prev_layer = self.last_layer - self.last_layer = new_layer # new layer becomes the last - else: # otherwise create a new layer to handle this - new_layer = CliffordLayer(gate) # a new layer with the new gate - # link to the layer structure - self.last_layer.next_layer = new_layer - new_layer.prev_layer = self.last_layer - self.last_layer = new_layer # new layer becomes the last - return self - elif isinstance(gate,MeasureLayer): - self.unitary = False - self.num_of_measures += len(gate.qubits) - self.last_layer.next_layer = gate - gate.prev_layer = self.last_layer - self.last_layer = gate + + def append_layer(self, layer): + # add a layer to the circuit (helper method not to be used directly) + # will not update unitary property and cached maps + if layer.isempty: # ignore empty layer return self - else: - raise NotImplementedError("Unknow input.") - def gate(self, *qubits): - if max(qubits)>=self.N: - raise ValueError("The gate acting on unregistered qubits!") - return self.take(CliffordGate(*qubits)) # create a new gate - def measure(self,*qubits): - if max(qubits)>=self.N: - raise ValueError("The gate acting on unregistered qubits!") - return self.take(MeasureLayer(*qubits,N = self.N)) # create measurement layer - def forward(self,obj): - if self.unitary: - if self.forward_map is None: - for layer in self.layers_forward(): - layer.forward(obj) + if self.last_layer.isempty: + # use the layer as the last layer + if self.last_layer is self.first_layer: + self.first_layer = layer + layer.prev_layer = self.last_layer.prev_layer + layer.next_layer = self.last_layer.next_layer + else: # otherwise, attach the layer to the last layer + self.last_layer.next_layer = layer + layer.prev_layer = self.last_layer + self.last_layer = layer + return self + + def append(self, op): + ''' Add an operation (gate or measurement) to the circuit. + This method tries to append the operation to the last layer if independent. + If not, a new Layer is created.''' + if isinstance(op, CliffordGate) or isinstance(op, Measurement): + # add an operation to the circuit + if self.last_layer.independent_from(op): + self.last_layer.append(op) else: - obj.transform_by(self.forward_map) - return obj + self.append_layer(Layer(op)) + elif isinstance(op, Layer): + # attach a layer to the circuit + self.append_layer(op) + elif isinstance(op, Circuit): + # compose a circuit to the circuit + for layer in op.layers_forward(): + self.append_layer(layer) else: + raise ValueError("Unsupported operation type {} to be appended to circuit.".format(type(op))) + # Update unitary property + self.unitary = self.unitary and op.unitary + # Clear cached maps since circuit had changed + self.forward_map = None + self.backward_map = None + return self + + def gate(self, *qubits): + '''Add a Clifford gate to the circuit''' + return self.append(CliffordGate(*qubits)) + + def measure(self, *qubits): + '''Add a measurement to the circuit''' + return self.append(Measurement(*qubits)) + + def forward(self, obj): + '''Apply the circuit forward to a quantum object. (inplace update)''' + log2prob = 0.0 + if self.forward_map is not None and self.unitary: + # if circuit forward map has been compiled, use it + obj.transform_by(self.forward_map) + else: # otherwise, apply each layer forward (in forward sequence) for layer in self.layers_forward(): - if not isinstance(layer,MeasureLayer): - layer.forward(obj) - else: - layer.forward(obj) - self.measure_result += layer.result.tolist() - self.log2prob += layer.log2prob - return obj - def backward(self,obj,measure_result = None): - if self.unitary: - if self.backward_map is None: - for layer in self.layers_backward(): - layer.backward(obj) - else: - obj.transform_by(self.backward_map) - return obj - else: - if measure_result is not None: - if len(measure_result)!=self.num_of_measures: - raise ValueError("classical bit does not match number of measurements.") - else: - pointer = 0 - for layer in self.layers_backward(): - if isinstance(layer,MeasureLayer): - new_pointer = pointer - len(layer.qubits) - if pointer == 0: - layer.backward(obj,measure_result = measure_result[new_pointer:]) - else: - layer.backward(obj,measure_result = measure_result[new_pointer:pointer]) - pointer = new_pointer - else: - layer.backward(obj) - return obj - else: - if len(self.measure_result)>0: - pointer = 0 - for layer in self.layers_backward(): - if isinstance(layer,MeasureLayer): - new_pointer = pointer - len(layer.qubits) - if pointer == 0: - layer.backward(obj,measure_result = self.measure_result[new_pointer:]) - else: - layer.backward(obj,measure_result = self.measure_result[new_pointer:pointer]) - pointer = new_pointer - else: - layer.backward(obj) - return obj - else: - raise ValueError("self measurement result is empty, and measurement result is not given!") - - def compile(self): - '''Construct forward and backward Clifford maps for this circuit + obj, layer_log2prob = layer.forward(obj) + log2prob += layer_log2prob + return obj, log2prob + + def backward(self, obj): + '''Apply the circuit backward to a quantum object. (inplace update)''' + log2prob = 0.0 + if self.backward_map is not None and self.unitary: + # if circuit backward map has been compiled, use it + obj.transform_by(self.backward_map) + else: # otherwise, apply each layer backward (in backward sequence) + for layer in self.layers_backward(): + obj, layer_log2prob = layer.backward(obj) + log2prob += layer_log2prob + return obj, log2prob + + def compile(self, N): + '''Compile the circuit into forward/backward maps where possible (unitary only). - Note: The compilation creates a single Clifford map representing the - entire circuit, which allows it to run faster for deep circuits.''' - N = self.N + Input: + N: int - number of qubits in the system''' if self.unitary: + for layer in self.layers_forward(): + layer.compile(N) self.forward_map = identity_map(N) self.backward_map = identity_map(N) for layer in self.layers_forward(): - layer.compile(N) self.forward_map = self.forward_map.compose(layer.forward_map) - self.backward_map = self.backward_map.compose(layer.backward_map) - return self - else: # has measurement, compile unitary layers - for layer in self.layers_forward(): - if not isinstance(layer, MeasureLayer): - layer.compile(N) - return self - def povm(self, nsample): - '''Assuming computational basis measurement follows the circuit, this - will back evolve the computational basis state to generate prior POVM. - This returns a generator.''' - if self.unitary: - for _ in range(nsample): - zero = zero_state(self.N) - yield self.backward(zero) - else: - raise NotImplementedError("The circuit has mid-circuit measurement. Therefore, it is unclear what is the POVM in the Heisenberg picture.") + self.backward_map = layer.backward_map.compose(self.backward_map) + return self # ---- gate constructors ---- def clifford_rotation_gate(generator, qubits=None): @@ -571,7 +538,7 @@ def identity_circuit(N = None): Parameters: N: int - number of qubits.''' - circ = CliffordCircuit(N) + circ = Circuit() if N is not None: circ.N = N # fix number of qubits explicitly return circ @@ -591,7 +558,7 @@ def brickwall_rcc(N, depth): def onsite_rcc(N): '''Construct random Clifford circuit of a layer of single-site gates. - (useful for implementing random Pauli measurements) + (useful for implementing random Pauli measurements) Parameters: N: int - number of qubits.''' @@ -602,7 +569,7 @@ def onsite_rcc(N): def global_rcc(N): '''Construct random Clifford circuit of a global Clifford gate. - (useful for implementing random Clifford measurements) + (useful for implementing random Clifford measurements) Parameters: N: int - number of qubits.''' @@ -610,176 +577,211 @@ def global_rcc(N): circ.gate(*range(N)) return circ -# ---- diagonalization ---- -def diagonalize(obj, i0 = 0, causal=False): - '''Diagonalize a Pauli operator or a stabilizer state (density matrix). - +def measurement_layer(N): + '''Construct a layer of computational basis measurements. + (to be combined with unitary circuit to form a measurement circuit) + Parameters: - obj: Pauli - the operator to be diagonalized, or - StabilizerState - the state to be diagonalized. - i0: int - index of the qubit to diagonalize to. - causal: bool - whether to preserve the causal structure by restricting - the action of Clifford transformation to the qubits at i0 and afterwards. - - Returns: - circ: CliffordCircuit - circuit that diagonalizes obj.''' - circ = identity_circuit(obj.N) - if isinstance(obj, (Pauli, PauliMonomial)): - if causal: - for g in pauli_diagonalize1(obj.g[2*i0:]): - circ.take(clifford_rotation_gate(Pauli(g), numpy.arange(i0,obj.N))) - else: - for g in pauli_diagonalize1(obj.g, i0): - circ.take(clifford_rotation_gate(Pauli(g))) - elif isinstance(obj, StabilizerState): - gate = CliffordGate(*numpy.arange(obj.N)) - gate.backward_map = obj.to_map() # set backward map to encoding map - # then the forward map automatically decodes (= diagonalize) the state - circ.take(gate) - else: - raise NotImplementedError('diagonalization is not implemented for {}.'.format(type(obj).__name__)) + N: int - number of qubits.''' + circ = identity_circuit(N) + circ.measure(*range(N)) return circ - -def SBRG(hmdl, max_rate=2., tol=1.e-8): - '''Approximately diagonalize a Hamiltonian by SBRG. - - Parameters: - hmdl: PauliPolynomial - model Hamiltonian. - - Returns: - heff: PauliPolynomial - MBL effective Hamiltonian. - circ: CliffordCircuit - Clifford circuit to diagonalize the Hamiltonian.''' - htmp = hmdl.copy() # copy of model Hamiltonian to workspace - N = htmp.N # system size - circ = identity_circuit(N) # create circuit - heff = pauli_zero(N) # create effective Hamiltonian - # SBRG iteration - for i0 in range(N): # pivot through every qubit - if len(htmp) == 0: - break - leading = numpy.argmax(numpy.abs(htmp.cs)) # get leading term index - # find circuit to diagonalize leading term to i0 - circ_i0 = diagonalize(htmp[leading], i0, causal=True) - circ.compose(circ_i0) # append it to total circuit - circ_i0.forward(htmp) # apply it to Hamiltonian - mask_commute = htmp.gs[:,2*i0] == 0 # mask diagonal terms - len_anti = sum(~mask_commute) # count number of offdiagonal terms - if len_anti != 0: # if offdiagonal terms exist - # split diagonal from offdiagonal terms - diag = htmp[mask_commute] - offdiag = htmp[~mask_commute] - # eleminate offdiagonal terms by perturbation theory - len_max = int(round(max_rate * len_anti)) # max perturbation terms to keep - prod = (offdiag @ offdiag).reduce(tol)[:len_max] - if len(prod) != 0: - htmp = diag + 0.5 * (htmp[leading].inverse() @ prod) - # mask terms that has become trivial on the remaining qubits - mask_trivial = numpy.all(htmp.gs[:,(2*i0+2):] == 0, -1) - heff += htmp[mask_trivial] # collect trivial terms to effective Hamiltonian - htmp = htmp[~mask_trivial] # retain with non-trivial terms - return heff, circ - -### gates #### +# ---- single qubit gates ---- def H(*qubits): if len(qubits)!=1: raise ValueError("Hadmand gate only acts on a single qubit.") gate = CliffordGate(*qubits) - f_map = CliffordMap(gs = np.array([[0,1],[1,0]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,0]]),ps = numpy.array([0,0])) gate.set_forward_map(f_map) return gate + def S(*qubits): if len(qubits)!=1: raise ValueError("S gate only acts on a single qubit.") gate = CliffordGate(*qubits) - f_map = CliffordMap(gs = np.array([[1,1],[0,1]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[1,1],[0,1]]),ps = numpy.array([0,0])) gate.set_forward_map(f_map) return gate + def X(*qubits): if len(qubits)!=1: raise ValueError("X gate only acts on a single qubit.") gate = CliffordGate(*qubits) - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([0,2])) gate.set_forward_map(f_map) return gate + def Y(*qubits): if len(qubits)!=1: raise ValueError("Y gate only acts on a single qubit.") gate = CliffordGate(*qubits) - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([2,2])) gate.set_forward_map(f_map) return gate + def Z(*qubits): if len(qubits)!=1: raise ValueError("Z gate only acts on a single qubit.") gate = CliffordGate(*qubits) - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([2,0])) gate.set_forward_map(f_map) return gate + def C(num, *qubits): # single qubit Clifford gate, num: 0-23 if len(qubits)!=1: raise ValueError("Single qubit Clifford gate acts on single qubit.") gate = CliffordGate(*qubits) if num == 0: - f_map = CliffordMap(gs = np.array([[1,0],[1,1]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[1,0],[1,1]]),ps = numpy.array([0,0])) elif num == 1: - f_map = CliffordMap(gs = np.array([[1,0],[1,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[1,1]]),ps = numpy.array([0,2])) elif num == 2: - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([0,0])) elif num == 3: - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([0,2])) elif num == 4: - f_map = CliffordMap(gs = np.array([[1,0],[1,1]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[1,0],[1,1]]),ps = numpy.array([2,0])) elif num == 5: - f_map = CliffordMap(gs = np.array([[1,0],[1,1]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[1,1]]),ps = numpy.array([2,2])) elif num == 6: - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([0,2])) elif num == 7: - f_map = CliffordMap(gs = np.array([[1,0],[0,1]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[1,0],[0,1]]),ps = numpy.array([2,2])) elif num == 8: - f_map = CliffordMap(gs = np.array([[1,1],[1,0]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[1,1],[1,0]]),ps = numpy.array([0,0])) elif num == 9: - f_map = CliffordMap(gs = np.array([[1,1],[1,0]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,1],[1,0]]),ps = numpy.array([0,2])) elif num == 10: - f_map = CliffordMap(gs = np.array([[1,1],[0,1]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[1,1],[0,1]]),ps = numpy.array([0,0])) elif num == 11: - f_map = CliffordMap(gs = np.array([[1,1],[0,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[1,1],[0,1]]),ps = numpy.array([0,2])) elif num == 12: - f_map = CliffordMap(gs = np.array([[1,1],[1,0]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[1,1],[1,0]]),ps = numpy.array([2,0])) elif num == 13: - f_map = CliffordMap(gs = np.array([[1,1],[1,0]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[1,1],[1,0]]),ps = numpy.array([2,2])) elif num == 14: - f_map = CliffordMap(gs = np.array([[1,1],[0,1]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[1,1],[0,1]]),ps = numpy.array([2,0])) elif num == 15: - f_map = CliffordMap(gs = np.array([[1,1],[0,1]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[1,1],[0,1]]),ps = numpy.array([2,2])) elif num == 16: - f_map = CliffordMap(gs = np.array([[0,1],[1,0]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,0]]),ps = numpy.array([0,0])) elif num == 17: - f_map = CliffordMap(gs = np.array([[0,1],[1,0]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,0]]),ps = numpy.array([0,2])) elif num == 18: - f_map = CliffordMap(gs = np.array([[0,1],[1,1]]),ps = np.array([0,0])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,1]]),ps = numpy.array([0,0])) elif num == 19: - f_map = CliffordMap(gs = np.array([[0,1],[1,1]]),ps = np.array([0,2])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,1]]),ps = numpy.array([0,2])) elif num == 20: - f_map = CliffordMap(gs = np.array([[0,1],[1,0]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,0]]),ps = numpy.array([2,0])) elif num == 21: - f_map = CliffordMap(gs = np.array([[0,1],[1,0]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,0]]),ps = numpy.array([2,2])) elif num == 22: - f_map = CliffordMap(gs = np.array([[0,1],[1,1]]),ps = np.array([2,0])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,1]]),ps = numpy.array([2,0])) elif num == 23: - f_map = CliffordMap(gs = np.array([[0,1],[1,1]]),ps = np.array([2,2])) + f_map = CliffordMap(gs = numpy.array([[0,1],[1,1]]),ps = numpy.array([2,2])) else: raise ValueError("There are only 24 single qubit Clifford gate. Input number exceed 0-23.") gate.set_forward_map(f_map) return gate + +# ---- two qubit gates ---- def CNOT(*qubits): if len(qubits)!=2: raise ValueError("CNOT gate acts on two qubit.") gate = CliffordGate(*qubits) if qubits[0] 0 + 0 - 0 = 0 (I) + # (1,0) -> 1 + 0 - 0 = 1 (X) + # (1,1) -> 1 + 3 - 2 = 2 (Y) + # (0,1) -> 0 + 3 - 0 = 3 (Z) + idx = self.g[2*i] + 3*self.g[2*i+1] - 2*self.g[2*i]*self.g[2*i+1] + matrices.append(sigma[idx]) + + # Compute tensor product + result = matrices[0] + for mat in matrices[1:]: + result = numpy.kron(result, mat) + + return (1j)**(self.p) * result class PauliList(object): '''Represents a list of Pauli operators. @@ -148,9 +168,10 @@ class PauliList(object): Parameters: gs: int (L, 2*N) - array of Pauli strings in binary repr. ps: int (L) - array of phase indicators (i powers).''' - def __init__(self, gs, ps = None): + def __init__(self, gs, ps=None, **kwargs): self.gs = gs self.ps = numpy.zeros(self.L, dtype=numpy.int_) if ps is None else ps + # kwargs ignored, in case subclass-specific arguments passed up def __repr__(self): return '\n'.join([repr(pauli) for pauli in self]) @@ -228,22 +249,47 @@ def transform_by(self, clifford_map, mask=None): def tokenize(self): return pauli_tokenize(self.gs, self.ps) - def to_qutip(self): - lists = [] - paulis = [qt.qeye(2),qt.sigmax(),qt.sigmay(),qt.sigmaz()] - for l in range(self.L): - tmp_list=[] - for i in range(self.N): - if (self.gs[l,2*i]==1)&(self.gs[l,2*i+1]==1): - tmp_list.append(paulis[2]) - elif (self.gs[l,2*i]==1)&(self.gs[l,2*i+1]==0): - tmp_list.append(paulis[1]) - elif (self.gs[l,2*i]==0)&(self.gs[l,2*i+1]==1): - tmp_list.append(paulis[3]) - else: - tmp_list.append(paulis[0]) - lists.append((1j)**(self.ps[l])*qt.tensor(tmp_list)) - return lists + + def to_numpy(self): + """Convert list of Pauli operators to numpy array representations in batch. + Returns a (L, 2^N, 2^N) array where L is the number of Pauli operators.""" + # Define all Pauli matrices as a single 4x2x2 array + sigma = numpy.array([ + [[1, 0], [0, 1]], # I (00) + [[0, 1], [1, 0]], # X (10) + [[0, -1j], [1j, 0]], # Y (11) + [[1, 0], [0, -1]] # Z (01) + ], dtype=complex) + + # Handle empty Pauli operator case + if self.N == 0: + return numpy.ones((self.L, 1, 1), dtype=complex) + + # For each qubit position, get the corresponding Pauli matrices for all operators + matrices = [] + for i in range(self.N): + # Map binary representation (x,z) to Pauli matrix index using formula: + # idx = x + 3z - 2xz for all operators at qubit i + idx = (self.gs[:,2*i] + 3*self.gs[:,2*i+1] - + 2*self.gs[:,2*i]*self.gs[:,2*i+1]) # shape: (L,) + # Select corresponding Pauli matrices for all operators + # sigma[idx] has shape (L,2,2) + matrices.append(sigma[idx]) + + # Compute tensor product for all operators simultaneously + result = matrices[0] # shape: (L,2,2) + for mat in matrices[1:]: + # Reshape for broadcasting: + # result: (L,m,m) -> (L,m,1,m,1) + # mat: (L,2,2) -> (L,1,2,1,2) + m = result.shape[1] + result = result.reshape(self.L, m, 1, m, 1) + mat = mat.reshape(self.L, 1, 2, 1, 2) + # Broadcast multiply and reshape back + result = (result * mat).reshape(self.L, m*2, m*2) + + # Apply phases + return (1j)**(self.ps[:,None,None]) * result class PauliMonomial(Pauli): '''Represent a Pauli operator with a coefficient. @@ -326,19 +372,10 @@ def as_polynomial(self): def inverse(self): return Pauli(self.g)/(self.c * 1j**self.p) - def to_qutip(self): - paulis = [qt.qeye(2), qt.sigmax(), qt.sigmay(), qt.sigmaz()] - tmp_list=[] - for i in range(self.g.shape[0]//2): - if (self.g[2*i]==1)&(self.g[2*i+1]==1): - tmp_list.append(paulis[2]) - elif (self.g[2*i]==1)&(self.g[2*i+1]==0): - tmp_list.append(paulis[1]) - elif (self.g[2*i]==0)&(self.g[2*i+1]==1): - tmp_list.append(paulis[3]) - else: - tmp_list.append(paulis[0]) - return self.c*(1j)**(self.p)*qt.tensor(tmp_list) + + def to_numpy(self): + """Convert Pauli monomial to numpy array representation.""" + return self.c * super().to_numpy() class PauliPolynomial(PauliList): '''Represent a linear combination of Pauli operators. @@ -346,7 +383,7 @@ class PauliPolynomial(PauliList): Parameters: gs: int (L, 2*N) - array of Pauli strings in binary repr. ps: int (L) - array of phase indicators (i powers). - cs: comlex - coefficients.''' + cs: comlex (L) - coefficients.''' def __init__(self, *args, **kwargs): super(PauliPolynomial, self).__init__(*args, **kwargs) self.cs = numpy.ones(self.ps.shape, dtype=numpy.complex_) # default coefficient @@ -425,22 +462,14 @@ def reduce(self, tol=1.e-10): cs = aggregate(self.cs * 1j**self.ps, inds, gs.shape[0]) mask = (numpy.abs(cs) > tol) return PauliPolynomial(gs[mask]).set_cs(cs[mask]) - def to_qutip(self): - paulis = [qt.qeye(2),qt.sigmax(),qt.sigmay(),qt.sigmaz()] - summation = 0 - for l in range(self.L): - tmp_list=[] - for i in range(self.N): - if (self.gs[l,2*i]==1)&(self.gs[l,2*i+1]==1): - tmp_list.append(paulis[2]) - elif (self.gs[l,2*i]==1)&(self.gs[l,2*i+1]==0): - tmp_list.append(paulis[1]) - elif (self.gs[l,2*i]==0)&(self.gs[l,2*i+1]==1): - tmp_list.append(paulis[3]) - else: - tmp_list.append(paulis[0]) - summation += self.cs[l]*(1j)**(self.ps[l])*qt.tensor(tmp_list) - return summation + + def to_numpy(self): + """Convert Pauli polynomial to numpy array representation. + Returns a (2^N, 2^N) array representing the sum of all terms.""" + # Get batch of Pauli matrices from parent class (shape: L x 2^N x 2^N) + matrices = super().to_numpy() + # Contract batch dimension with coefficients to get final matrix + return numpy.tensordot(self.cs, matrices, axes=(0,0)) # ---- constructors ---- def pauli(obj, N = None): @@ -513,4 +542,4 @@ def pauli_identity(N): def pauli_zero(N): '''Pauli polynomial of zero operator of N qubit''' - return 0 * pauli_identity(N) + return 0 * pauli_identity(N) \ No newline at end of file diff --git a/pyclifford/stabilizer.py b/pyclifford/stabilizer.py index 9da4eef..673fa47 100644 --- a/pyclifford/stabilizer.py +++ b/pyclifford/stabilizer.py @@ -1,11 +1,10 @@ import numpy -import qutip as qt from numba import njit from .utils import ( acq_mat, ps0, z2inv, pauli_combine, pauli_transform, binary_repr, random_pauli, random_clifford, map_to_state, state_to_map, clifford_rotate, - stabilizer_project, stabilizer_measure, stabilizer_expect, - stabilizer_entropy, mask, stabilizer_projection_trace,stabilizer_postselection) + stabilizer_measure, stabilizer_postselect, stabilizer_project, stabilizer_expect, + stabilizer_entropy, mask) from .paulialg import Pauli, PauliList, PauliPolynomial, pauli, paulis class CliffordMap(PauliList): @@ -24,6 +23,7 @@ class CliffordMap(PauliList): gs: int (2*N, 2*N) - strings of Pauli operators to be mapped to. ps: int (2*N) - phase indicators of Pauli operators to be mapped to.''' def __init__(self, *args, **kwargs): + # call superclass PauliList to handle arguments super(CliffordMap, self).__init__(*args, **kwargs) def __repr__(self): @@ -41,11 +41,11 @@ def __repr__(self): def copy(self): return CliffordMap(self.gs.copy(), self.ps.copy()) - def to_state(self, r=None): + def to_state(self, r=0): '''Interprete the Clifford map as a stabilizer state, such that the state is generated by the map from the zero state.''' gs, ps = map_to_state(self.gs, self.ps) - return StabilizerState(gs, ps).set_r(r) + return StabilizerState(gs, ps, r=r) def embed(self, small_map, mask): '''Embed a smaller map acting on a subsystem specified by qubit indices.''' @@ -85,11 +85,14 @@ class StabilizerState(PauliList): Parameters: gs: int (2*N, 2*N) - strings of Pauli operators in the stabilizer tableau. ps: int (2*N) - phase indicators (should only be 0 or 2). - r: int - number of logical qubits (log2 rank of density matrix)''' - # def __init__(self, *args, **kwargs): - def __init__(self, gs, r=0, **kwargs): - super(StabilizerState, self).__init__(gs, **kwargs) - self.r = r # pure state by default + r: int - number of logical qubits (log2 rank of density matrix) + (r can only be provided as a keyword argument)''' + def __init__(self, *args, **kwargs): + # extract r and remove it from kwargs, if present. + # otherwise, set r = 0 as pure state by default. + self.r = kwargs.pop('r', 0) + # call superclass PauliList to handle remaining arguments + super().__init__(*args, **kwargs) def __repr__(self): ''' will only show active stabilizers, @@ -105,81 +108,142 @@ def stabilizers(self): return self[self.r:self.N] def copy(self): - return StabilizerState(self.gs.copy(), self.ps.copy()).set_r(self.r) - - def set_r(self, r=None): - '''set log2 rank of the density matrix.''' - self.r = 0 if r is None else r - return self + return StabilizerState(self.gs.copy(), self.ps.copy(), r=self.r) def to_map(self): - '''Interprete the stabilizer code as its encoding Clifford map.''' + '''Interprete the stabilizer state as its encoding Clifford map.''' gs, ps = state_to_map(self.gs, self.ps) return CliffordMap(gs, ps) def measure(self, obs): - '''Perform measurement on the stabilizer state. + '''Perform Pauli observable measurement on the stabilizer state. + - sample measurement outcomes from: + p(out_k = 0|rho) = (1 + Tr(rho obs_k))/2 + p(out_k = 1|rho) = (1 - Tr(rho obs_k))/2 + - collapse the state to the posterior state conditioned on + the measurement outcomes. (in-place update) Parameters: obs: PauliList or StabilizerState (only active stabilizers measured) Returns: - out: int (L) - array of measurement outcomes of each observable. - log2prob: real - log2 probability of this set of outcomes.''' + out: int (L) - array of measurement outcomes of corresponding observables. + log2prob: real - log2 probability of sampling this set of outcomes.''' if isinstance(obs, StabilizerState): obs = obs.stabilizers self.gs, self.ps, self.r, out, log2prob = stabilizer_measure( self.gs, self.ps, obs.gs, obs.ps, self.r) return out, log2prob - def postselect(self, paulistring, postselect_res): - ''' - paulistring: is the measurement - postselect_res: is the post-selection result. 0 means (+1), 1 means (-1) - ''' - if self.r != 0: - raise ValueError("Currently, post-selection is only supported with pure states") - self.gs, self.ps, prob = stabilizer_postselection(self.gs, self.ps, paulistring.g, int(postselect_res*2)) - return prob - def expect(self, obs): + def postselect(self, obs, out=None): + '''Postselect the stabilizer state on a set of Pauli observables, + such that on the postselected state: + - default behavior: = +1, + - with target outcomes: = (-1)^{out_k}, + The postselected state is in-place updated. + + Parameters: + obs: PauliList or StabilizerState (only active stabilizers postselected) + out: int (L) - array of target measurement outcomes of corresponding observables. + default is None, meaning = +1. + + Returns: + log2prob: real - log2 probability for postselection to succeed.''' + if isinstance(obs, StabilizerState): + obs = obs.stabilizers + if out is None: + obs_ps = obs.ps + else: # encode target outcomes to operator phases + obs_ps = (obs.ps + 2*out)%4 # modify operator phases + self.gs, self.ps, self.r, log2prob = stabilizer_postselect( + self.gs, self.ps, obs.gs, obs_ps, self.r) + return log2prob + + def expect(self, obs, z=1.): '''Evaluate expectation values of observables on the statilizer state. Parameters: obs: observable, can be Pauli, PauliList, PauliPolynomial, StabilizerState z: fugacity of operator weight, it is not used when obs is StabilizerState + default is 1. Setting z = 3 for Pauli shadow tomography directly. Returns: out: output (depending on the type of obs) * Pauli: promote to PauliPolynomial * PauliPolynomial O: Tr(rho O z^|O|) - * StabilizerState sigma: Tr(rho sigma) + * StabilizerState sigma: Tr(rho sigma) = p(sigma|rho) * PauliList [O_i]: [Tr(rho O_i z^|O_i|)] ''' if isinstance(obs, Pauli): return self.expect(obs.as_polynomial()) # cast Pauli to PauliPolynomial elif isinstance(obs, PauliPolynomial): - xs = self.expect(PauliList(obs.gs, obs.ps)) # cast PauliPolynomial to PauliList - return numpy.sum(obs.cs * xs) + # cast PauliPolynomial to PauliList for expectation value calculation + xs = self.expect(PauliList(obs.gs, obs.ps)) + if z != 1.: # if fugacity not 1, multiply by fugacity to the power of operweight + zs = z ** obs.weight() + xs = xs * zs + return numpy.sum(obs.cs * xs) # combine expectation values by coefficients elif isinstance(obs, StabilizerState): - if self.r!=0: - raise NotImplementedError("Will be added in the next release!") - else: - _, _, _, trace = stabilizer_projection_trace(numpy.array(self.gs), numpy.array(self.ps), \ - numpy.array(obs.gs[obs.r:obs.N,:]), numpy.array(obs.ps[obs.r:obs.N]), 0) - return trace/2**obs.r - + rho = self.copy() # make a copy to avoid in-place update in postselection + log2prob = rho.postselect(obs) # use post-selection to calculate Tr(rho obs) + return 2. ** log2prob # convert log2 probability to probability + # WARNING: PauliList instance must be placed after StabilizerState instance + # otherwise StabilizerState will be shadowed by PauliList as subclass elif isinstance(obs, PauliList): xs = stabilizer_expect(self.gs, self.ps, obs.gs, obs.ps, self.r) + if z != 1.: # if fugacity not 1, multiply by fugacity to the power of operator weight + zs = z ** obs.weight() + xs = xs * zs return xs - def to_qutip(self): - ID = qt.tensor([qt.qeye(2) for i in range(self.N)]) - rho = ID - for i in range(self.r,self.N): - rho = rho*(ID+Pauli(self.gs[i],self.ps[i]).to_qutip())/2 - # rho = rho*(ID+pauli2pauli(state.gs[i],state.ps[i]))/2 - rho = rho/(2**self.r) - return rho - + else: + raise ValueError("Unsupported observable type: {}".format(type(obs))) + + def to_numpy(self): + """Convert stabilizer state to numpy density matrix representation. + Returns a (2^N, 2^N) array representing rho = 1/2^r prod_{a=1}^{N-r} (1+ Pauli[g_a,p_a])/2 + """ + # Define all Pauli matrices as a single 4x2x2 array + sigma = numpy.array([ + [[1, 0], [0, 1]], # I (00) + [[0, 1], [1, 0]], # X (10) + [[0, -1j], [1j, 0]], # Y (11) + [[1, 0], [0, -1]] # Z (01) + ], dtype=complex) + # Handle empty state case + if self.N == 0: + return numpy.ones((1, 1), dtype=complex) + # Only process active stabilizers + active_gs = self.gs[self.r:self.N] # Shape: (N-r, 2N) + active_ps = self.ps[self.r:self.N] # Shape: (N-r,) + # For each active stabilizer, get its Pauli matrices + matrices = [] + for i in range(self.N): + # Map binary representation (x,z) to Pauli matrix index for all active stabilizers + # This maps: (0,0)->0 (I), (1,0)->1 (X), (1,1)->2 (Y), (0,1)->3 (Z) + idx = (active_gs[:,2*i] + 3*active_gs[:,2*i+1] - + 2*active_gs[:,2*i]*active_gs[:,2*i+1]) # Shape: (N-r,) + # Select corresponding Pauli matrices for all active stabilizers + matrices.append(sigma[idx]) # Shape: (N-r, 2, 2) + # Compute tensor product for each active stabilizer simultaneously + result = matrices[0] # Shape: (N-r, 2, 2) + for mat in matrices[1:]: + # Reshape for broadcasting: + # result: (N-r,m,m) -> (N-r,m,1,m,1) + # mat: (N-r,2,2) -> (N-r,1,2,1,2) + m = result.shape[1] + result = result.reshape(self.N-self.r, m, 1, m, 1) + mat = mat.reshape(self.N-self.r, 1, 2, 1, 2) + # Broadcast multiply and reshape back + result = (result * mat).reshape(self.N-self.r, m*2, m*2) + # Apply phases to get Pauli operators + pauli_ops = (1j)**(active_ps[:,None,None]) * result # Shape: (N-r, 2^N, 2^N) + # Construct density matrix using the formula + rho = numpy.eye(2**self.N, dtype=complex) # Start with identity + for op in pauli_ops: + rho = rho @ (numpy.eye(2**self.N) + op) / 2 + # Apply normalization factor + return rho / (2**self.r) + def entropy(self, subsys): '''Entanglement entropy of the stabilizer state in a given region.''' if isinstance(subsys, (tuple, list)): @@ -199,15 +263,13 @@ def sample(self, L): C = numpy.random.randint(2, size=(L,self.N-self.r)) gs, ps = pauli_combine(C, self.gs[self.r:self.N], self.ps[self.r:self.N]) return PauliList(gs, ps) - def get_prob(self, readout): - ''' - Evaluate the probability of getting a bit string readout - ''' - if self.N != readout.shape[0]: - raise Error("readout is incompitile with system size!") - readout_state = identity_map(self.N).to_state() - readout_state.ps[:self.N]=readout - return self.expect(readout_state) + + def get_prob(self, out): + '''Evaluate the probability of getting a bit string readout. + Assume computational basis measurement.''' + if self.N != out.shape[0]: + raise ValueError("readout size {} is incompatible with system size {}!".format(out.shape[0], self.N)) + return self.expect(bit_state(self.N, out)) # !!! this function has exponential complexity. @property @@ -292,15 +354,28 @@ def one_state(N): ps = (2*numpy.ones(2*N)).astype(int) return StabilizerState(gs = gs,ps = ps) +def bit_state(N, bits): + # bits should be a binary array of length N. + # if bits is represented as integer or string, convert to array + if isinstance(bits, int): + bits = numpy.array(list(numpy.binary_repr(bits, width=N))).astype(int) + elif isinstance(bits, str): + if len(bits) != N: + raise ValueError("bitstring must be of length N") + bits = numpy.array(list(bits)).astype(int) + gs = zero_state(N).gs + ps = 2*bits + return StabilizerState(gs = gs, ps = ps) + def ghz_state(N): objs = [pauli({i:3,i+1:3},N) for i in range(N-1)] objs.append(pauli([1]*N)) return stabilizer_state(paulis(objs)) -def random_pauli_state(N, r=None): +def random_pauli_state(N, r=0): return random_pauli_map(N).to_state(r) -def random_clifford_state(N, r=None): +def random_clifford_state(N, r=0): return random_clifford_map(N).to_state(r) @njit @@ -311,6 +386,7 @@ def random_bit_state_gs_ps(N): gs[N+i,2*i]=1 ps = numpy.random.choice(numpy.array([0,2]),size = 2*N) return gs, ps + def random_bit_state(N): gs, ps = random_bit_state_gs_ps(N) return StabilizerState(gs = gs.astype(int), ps = ps) \ No newline at end of file diff --git a/pyclifford/utils.py b/pyclifford/utils.py index 4632eb7..73b9194 100644 --- a/pyclifford/utils.py +++ b/pyclifford/utils.py @@ -2,10 +2,10 @@ from numba import njit '''Conventions: -Binary represention of Pauli string. (arXiv:quant-ph/0406196) +Binary representation of Pauli string. (arXiv:quant-ph/0406196) A N-qubit Pauli string is represented as a (2*N)-dimensional binary - vector of the form - g = [x0,z0;x1,z1;...;x(N-1),z(N-1)] + array of the following form: + g = [x0,z0;x1,z1;...;x(N-1),z(N-1)] (flattened) which corresponds to sigma[g] = i^(x.z) prod_i X_i^xi Z_i^zi Commutation relation of Pauli strings: @@ -19,7 +19,7 @@ def front(g): '''Find the first nontrivial qubit in a Pauli string. Parameters: - g: int (2*N) - a Pauli string in binary repr. + g: int (2*N) - a Pauli string in binary representation. Returns: i: int - position of its first nontrivial qubit. @@ -40,10 +40,10 @@ def condense(g): a shorter string and the support. Parameters: - g: int (2*N) - a Pauli string in binary repr. + g: int (2*N) - a Pauli string in binary representation. Returns: - g_cond: int (2*n) - the condensed Pauli string in binary repr. + g_cond: int (2*n) - the condensed Pauli string in binary representation. qubits: int (n) - indices of supporting qubits.''' (N2,) = g.shape N = N2//2 @@ -59,7 +59,7 @@ def p0(g): '''Bare phase factor due to x.z for a Pauli string. Parameters: - g: int (2*N) - a Pauli string in binary repr. + g: int (2*N) - a Pauli string in binary representation. Returns: p0: int - bare phase factor x.z for the string.''' @@ -75,8 +75,8 @@ def acq(g1, g2): '''Calculate Pauli operator anticmuunation indicator. Parameters: - g1: int (2*N) - the first Pauli string in binary repr. - g2: int (2*N) - the second Pauli string in binary repr. + g1: int (2*N) - the first Pauli string in binary representation. + g2: int (2*N) - the second Pauli string in binary representation. Returns: acq: int - acq = 0 if g1, g2 commute, acq = 1 if g1, g2 anticommute.''' @@ -93,8 +93,8 @@ def ipow(g1, g2): '''Phase indicator for the product of two Pauli strings. Parameters: - g1: int (2*N) - the first Pauli string in binary repr. - g2: int (2*N) - the second Pauli string in binary repr. + g1: int (2*N) - the first Pauli string in binary representation. + g2: int (2*N) - the second Pauli string in binary representation. Returns: ipow: int - the phase indicator (power of i) when product @@ -118,7 +118,7 @@ def ps0(gs): '''Bare phase factor due to x.z for Pauli strings. Parameters: - gs: int (L,2*N) - array of Pauli strings in binary repr. + gs: int (L,2*N) - array of Pauli strings in binary representation. Returns: ps0: int (L) - bare phase factor x.z for all strings.''' @@ -135,7 +135,7 @@ def acq_mat(gs): '''Construct anticommutation indicator matrix for a set of Pauli strings. Parameters: - gs: int (L,2*N) - array of Pauli strings in binary repr. + gs: int (L,2*N) - array of Pauli strings in binary representation. Returns: mat: int (L,L) - anticommutation indicator matrix.''' @@ -186,7 +186,7 @@ def pauli_tokenize(gs, ps): '''Create a token of Pauli operators for learning tasks. Parameters: - gs: int (L, 2*N) - Pauli strings in binary repr. + gs: int (L, 2*N) - Pauli strings in binary representation. ps: int (L) - phase indicators. Returns: @@ -210,11 +210,11 @@ def pauli_combine(C, gs_in, ps_in): Parameters: C: int (L_out, L_in) - one-hot encoding of selected operators. - gs_in: int (L_in, 2*N) - input binary repr of Pauli strings. + gs_in: int (L_in, 2*N) - input binary representation of Pauli strings. ps_in: int (L_in) - phase indicators of input operators. Returns: - gs_out: int (L_out, 2*N) - output binary repr of Pauli strings. + gs_out: int (L_out, 2*N) - output binary representation of Pauli strings. ps_out: int (L_out) - phase indicators of output operators. ''' (L_out, L_in) = C.shape @@ -234,13 +234,13 @@ def pauli_transform(gs_in, ps_in, gs_map, ps_map): (right multiplication) Parameters: - gs_in: int (L, 2*N) - input binary repr of Pauli strings. + gs_in: int (L, 2*N) - input binary representation of Pauli strings. ps_in: int (L) - phase indicators of input operators. gs_map: int (2*N, 2*N) - operator map in binary representation. ps_map: int (2*N) - phase indicators associated to target operators. Returns: - gs_out: int (L, 2*N) - output binary repr of Pauli strings. + gs_out: int (L, 2*N) - output binary representation of Pauli strings. ps_out: int (L) - phase indicators of output operators.''' gs_out, ps_out = pauli_combine(gs_in, gs_map, ps_map) ps_out = (ps_in + ps0(gs_in) + ps_out)%4 @@ -252,9 +252,9 @@ def clifford_rotate(g, p, gs, ps): '''Apply Clifford rotation to Pauli operators. Parameters: - g: int (2*N) - Clifford rotation generator in binary repr. + g: int (2*N) - Clifford rotation generator in binary representation. p: int - phase indicator (p = 0, 2 only). - gs: int (L, 2*N) - input binary repr of Pauli strings. + gs: int (L, 2*N) - input binary representation of Pauli strings. ps: int (L) - phase indicators of input operators. Returns: gs, ps in-place modified.''' @@ -270,8 +270,8 @@ def clifford_rotate_signless(g, gs): '''Apply Clifford rotation to Pauli strings without signs. Parameters: - g: int (2*N) - Clifford rotation generator in binary repr. - gs: int (L, 2*N) - array of Pauli strings in binary repr. + g: int (2*N) - Clifford rotation generator in binary representation. + gs: int (L, 2*N) - array of Pauli strings in binary representation. Returns: gs in-place modified.''' (L, N2) = gs.shape @@ -307,7 +307,7 @@ def pauli_diagonalize1(g1, i0 = 0): to qubit i0 as Z. Parameters: - g1: int (2*N) - Pauli string in binary repr. + g1: int (2*N) - Pauli string in binary representation. i0: int - target qubit Returns: @@ -339,6 +339,7 @@ def pauli_diagonalize1(g1, i0 = 0): def pauli_diagonalize2(g1, g2, i0 = 0): '''Find a series of Clifford roations to diagonalize a pair of anticommuting Pauli strings to qubit i0 as Z and X (or Y). + Parameters: g1: int (2*N) - binary representation of stabilizer. g2: int (2*N) - binary representation of destabilizer. @@ -498,36 +499,129 @@ def state_to_map(gs_in, ps_in): return gs_out, ps_out # ---- stabilizer related ---- +''' Formulation: +A stabilizer state encoding r logical qubits on N physical qubits +is described by the density matrix of the following form: + rho = sum_{i=1}^{N-r} 2^{-r} (1 + S_i)/2. +where +* S_i are Pauli stabilizers (commuting with each other), + whose Pauli strings are stored in gs_stb[r:N] in binary representation, + and phase indicators are stored in ps_stb[r:N]. +* r is the log2 rank of the density matrix, i.e. the number of standby stabilizers. + +Structure of Stabilizer Tableau: +* gs_stb[0:r]: Pauli strings of standby stabilizers in binary representation +* gs_stb[r:N]: Pauli strings of active stabilizers in binary representation +* gs_stb[N:N+r]: Pauli strings of standby destabilizers in binary representation +* gs_stb[N+r:2*N]: Pauli strings of active destabilizers in binary representation +* ps_stb[0:r]: phase indicators of standby stabilizers +* ps_stb[r:N]: phase indicators of active stabilizers +* ps_stb[N:N+r]: phase indicators of standby destabilizers +* ps_stb[N+r:2*N]: phase indicators of active destabilizers +Symplectic structure: each stabilizer gs_stb[i] only anticommutes with + its destabilizer gs_stb[N+i]. + +Algorithm: +Stabilizer states can serve both as the prior density matrix +and as the set of commuting observables. +- As state: described by density matrix + rho = sum_{i=1}^{N-r} 2^{-r} (1 + S_i)/2 + where S_i are the active stabilizers. +- As observable: described by projective measurement operator + pi(o) = sum_{k=1}^{L} 2^{-N+L} (1 + (-)^o_k O_k)/2 + where O_k are the commuting observables, and o_k is the measurement outcome. +---- measure ---- +When measuring {O_k} on a stabilizer state of {S_i}, the algorithm is as follows: +for O_k in {O_k}: + if O_k anticommutes with: + case 1: any active stabilizer (the first of which being S_p) + case 2: any standy stabilizer or destabilizer (the first of which being S_p) + case 3: otherwise + then: + case 1: O_k is an error operator -> update + readout + case 2: O_k is a logical operator -> update + extend + readout + case 3: O_k is a trivial operator -> readout + if update: + - O_k replaces S_p to be the active stabilizer, + - S_p becomes its corresponding destabilizer, + - for any other stabilizers and destabilizers that anticommute with O_k, + update them to commute with O_k by multiplying O_k to them, + - keep track of the phase accumulated in the process, [1] + if extend: + - rank r is reduced by 1, + - include the new stabilizer S_p to active stabilizers, + - include the new destabilizer S_q to active destabilizers. + - sign of O_k is randomly sampled with half-to-half probability. [2] + - record measurement outcome, [3] + - update log2 probability. [4] + else: + - readout measurement outcome, [5] + - log2 probability = 0. [6] +--- postselect --- +Same as measure, but: +- [2] measurement outcome not sampled butcopied from observation value. +- [3] omitted. +- [6] if stabilizer readout not match observation value, log2 probability = -inf. +--- project --- +Same as measure, but lines [1-6] are omitted. +''' @njit -def stabilizer_project(gs_stb, gs_obs, r): - '''Project stabilizer tableau to a new stabilizer basis. +def stabilizer_measure(gs_stb, ps_stb, gs_obs, ps_obs, r): + '''Measure a set of commuting Pauli observables on a stabilizer state. + + Given the prior state rho = sum_{i=1}^{N-r} 2^{-r} (1 + S_i)/2, + specified by (gs_stb, ps_stb, r), and a set of commuting observables + {O_k}, specified by (gs_obs, ps_obs), sample measurement outcomes + out := {o_k} from the conditional probability distribution: + p(pi(o)|rho) := Tr(rho pi(o)) + where pi(o) = sum_{k=1}^{L} 2^{-N+L} (1 + (-)^o_k O_k)/2 is the + measurement operator. Upon measurement, the state is updated to: + rho' = pi(o) rho pi(o) / p(pi(o)|rho). + Returns rho', o, log2 p(pi(o)|rho). Parameters: gs_stb: int (2*N, 2*N) - Pauli strings in original stabilizer tableau. - gs_obs: int (L, 2*N) - Pauli strings of new stablizers to impose. + ps_stb: int (N) - phase indicators of (de)stabilizers. + gs_obs: int (L, 2*N) - strings of Pauli operators to be measured. + ps_obs: int (L) - phase indicators of Pauli operators to be measured. r: int - log2 rank of density matrix (num of standby stablizers). Returns: gs_stb: int (2*N, 2*N) - Pauli strings in updated stabilizer tableau. - r: int - updated log2 rank of density matrix.''' + ps_stb: int (N) - phase indicators of (de)stabilizers. + r: int - updated log2 rank of density matrix. + out: int (L) - measurment outcomes (0 or 1 binaries). + log2prob: real - log2 probability of this outcome.''' (L, Ng) = gs_obs.shape N = Ng//2 assert 0<=r<=N - for k in range(L): # loop over incoming projections gs_obs[k] + out = numpy.empty(L, dtype=numpy.int_) + ga = numpy.empty(2*N, dtype=numpy.int_) # workspace for stabilizer accumulation + pa = 0 # workspace for phase accumulation + log2prob = 0. + for k in range(L): # for each observable gs_obs[k] update = False extend = False - p = 0 # pointer + p = 0 # pointer to the first anticommuting operator + ga[:] = 0 + pa = 0 for j in range(2*N): if acq(gs_stb[j], gs_obs[k]): # find gs_stb[j] anticommute with gs_obs[k] if update: # if gs_stb[j] is not the first anticommuting operator - gs_stb[j] = (gs_stb[j] + gs_stb[p])%2 # update gs_stb[j] to commute with gs_obs[k] + # update gs_stb[j] to commute with gs_obs[k] + if j < N: # if gs_stb[j] is a stablizer, phase matters + ps_stb[j] = (ps_stb[j] + ps_stb[p] + ipow(gs_stb[j], gs_stb[p]))%4 + gs_stb[j] = (gs_stb[j] + gs_stb[p])%2 else: # if gs_stb[j] is the first anticommuting operator if j < N + r: # if gs_stb[j] is not an active destabilizer p = j # move pointer to j update = True if not r <= j < N: # if gs_stb[j] is a standby operator extend = True - # if gs_stb[j] is an active destablizer, gs_obs[k] alreay a combination of active stablizers, do nothing. + else: # gs_stb[j] an active destabilizer, meaning gs_obs[k] already a combination of active stabilizers + # collect corresponding stabilizer component to ga + pa = (pa + ps_stb[j-N] + ipow(ga, gs_stb[j-N]))%4 + ga = (ga + gs_stb[j-N])%2 if update: # now gs_stb[p] and gs_obs[k] anticommute q = (p+N)%(2*N) # get q as dual of p @@ -544,36 +638,49 @@ def stabilizer_project(gs_stb, gs_obs, r): s = (r+N)%(2*N) # get s as dual of r gs_stb[numpy.array([p,r])] = gs_stb[numpy.array([r,p])] # swap p,r gs_stb[numpy.array([q,s])] = gs_stb[numpy.array([s,q])] # swap q,s - return gs_stb, r + p = r + # as long as gs_obs[k] is not eigen, outcome will be half-to-half + ps_stb[p] = 2 * numpy.random.randint(2) + out[k] = ((ps_stb[p] - ps_obs[k])%4)//2 #0->0(+1 eigenvalue), 2->1(-1 eigenvalue) + log2prob -= 1. + else: # no update, gs_obs[k] is eigen, result is in pa + assert((ga == gs_obs[k]).all()) + out[k] = ((pa - ps_obs[k])%4)//2 + return gs_stb, ps_stb, r, out, log2prob @njit -def stabilizer_measure(gs_stb, ps_stb, gs_obs, ps_obs, r): - '''Measure Pauli operators on a stabilizer state. +def stabilizer_postselect(gs_stb, ps_stb, gs_obs, ps_obs, r): + '''Postselect stabilizer state given a set of Pauli observations. + + Given the prior state rho = sum_{i=1}^{N-r} 2^{-r} (1 + S_i)/2, + specified by (gs_stb, ps_stb, r), and a measurement operator + sigma = sum_{k=1}^{L} 2^{-N+L} (1 + O_k)/2. Postselect the state by + rho' = sigma rho sigma / p(sigma|rho), + where p(sigma|rho) = Tr(rho sigma) is the success probability. + Returns rho', log2 p(sigma|rho). Parameters: gs_stb: int (2*N, 2*N) - Pauli strings in original stabilizer tableau. ps_stb: int (N) - phase indicators of (de)stabilizers. - gs_obs: int (L, 2*N) - strings of Pauli operators to be measured. - ps_obs: int (L) - phase indicators of Pauli operators to be measured. + gs_obs: int (L, 2*N) - strings of Pauli operators to be postselected. + ps_obs: int (L) - phase indicators of Pauli operators to be postselected. r: int - log2 rank of density matrix (num of standby stablizers). Returns: gs_stb: int (2*N, 2*N) - Pauli strings in updated stabilizer tableau. ps_stb: int (N) - phase indicators of (de)stabilizers. r: int - updated log2 rank of density matrix. - out: int (L) - measurment outcomes (0 or 1 binaries). - log2prob: real - log2 probability of this outcome.''' + log2prob: real - log2 probability of successful postselection.''' (L, Ng) = gs_obs.shape N = Ng//2 assert 0<=r<=N - out = numpy.empty(L, dtype=numpy.int_) ga = numpy.empty(2*N, dtype=numpy.int_) # workspace for stabilizer accumulation pa = 0 # workspace for phase accumulation log2prob = 0. for k in range(L): # for each observable gs_obs[k] update = False extend = False - p = 0 # pointer + p = 0 # pointer to the first anticommuting operator ga[:] = 0 pa = 0 for j in range(2*N): @@ -589,7 +696,7 @@ def stabilizer_measure(gs_stb, ps_stb, gs_obs, ps_obs, r): update = True if not r <= j < N: # if gs_stb[j] is a standby operator extend = True - else: # gs_stb[j] anticommute with destabilizer, meaning gs_obs[k] already a combination of active stabilizers + else: # gs_stb[j] an active destabilizer, meaning gs_obs[k] already a combination of active stabilizers # collect corresponding stabilizer component to ga pa = (pa + ps_stb[j-N] + ipow(ga, gs_stb[j-N]))%4 ga = (ga + gs_stb[j-N])%2 @@ -610,14 +717,62 @@ def stabilizer_measure(gs_stb, ps_stb, gs_obs, ps_obs, r): gs_stb[numpy.array([p,r])] = gs_stb[numpy.array([r,p])] # swap p,r gs_stb[numpy.array([q,s])] = gs_stb[numpy.array([s,q])] # swap q,s p = r - # as long as gs_obs[k] is not eigen, outcome will be half-to-half - ps_stb[p] = 2 * numpy.random.randint(2) - out[k] = ((ps_stb[p] - ps_obs[k])%4)//2 #0->0(+1 eigenvalue), 2->1(-1 eigenvalue) + # stabilizer sign set by observation value via postselection + ps_stb[p] = ps_obs[k] log2prob -= 1. else: # no update, gs_obs[k] is eigen, result is in pa assert((ga == gs_obs[k]).all()) - out[k] = ((pa - ps_obs[k])%4)//2 - return gs_stb, ps_stb, r, out, log2prob + if (pa - ps_obs[k])%4 != 0: # if result not match observation + log2prob -= numpy.inf # log likelihood -inf + return gs_stb, ps_stb, r, log2prob + +@njit +def stabilizer_project(gs_stb, gs_obs, r): + '''Project stabilizer tableau to a new stabilizer basis. + + Parameters: + gs_stb: int (2*N, 2*N) - Pauli strings in original stabilizer tableau. + gs_obs: int (L, 2*N) - Pauli strings of new stablizers to impose. + r: int - log2 rank of density matrix (num of standby stablizers). + + Returns: + gs_stb: int (2*N, 2*N) - Pauli strings in updated stabilizer tableau. + r: int - updated log2 rank of density matrix.''' + (L, Ng) = gs_obs.shape + N = Ng//2 + assert 0<=r<=N + for k in range(L): # loop over incoming projections gs_obs[k] + update = False + extend = False + p = 0 # pointer to the first anticommuting operator + for j in range(2*N): + if acq(gs_stb[j], gs_obs[k]): # find gs_stb[j] anticommute with gs_obs[k] + if update: # if gs_stb[j] is not the first anticommuting operator + gs_stb[j] = (gs_stb[j] + gs_stb[p])%2 # update gs_stb[j] to commute with gs_obs[k] + else: # if gs_stb[j] is the first anticommuting operator + if j < N + r: # if gs_stb[j] is not an active destabilizer + p = j # move pointer to j + update = True + if not r <= j < N: # if gs_stb[j] is a standby operator + extend = True + # if gs_stb[j] is an active destablizer, gs_obs[k] alreay a combination of active stablizers, do nothing. + if update: + # now gs_stb[p] and gs_obs[k] anticommute + q = (p+N)%(2*N) # get q as dual of p + gs_stb[q] = gs_stb[p] # move gs_stb[p] to gs_stb[q] + gs_stb[p] = gs_obs[k] # add gs_obs[k] to gs_stb[p] + if extend: + r -= 1 # rank will reduce under extension + # bring new stabilizer from p to r + if p == r: + pass + elif q == r: + gs_stb[numpy.array([p,q])] = gs_stb[numpy.array([q,p])] # swap p,q + else: + s = (r+N)%(2*N) # get s as dual of r + gs_stb[numpy.array([p,r])] = gs_stb[numpy.array([r,p])] # swap p,r + gs_stb[numpy.array([q,s])] = gs_stb[numpy.array([s,q])] # swap q,s + return gs_stb, r @njit def stabilizer_expect(gs_stb, ps_stb, gs_obs, ps_obs, r): @@ -654,7 +809,7 @@ def stabilizer_expect(gs_stb, ps_stb, gs_obs, ps_obs, r): if trivial: xs[k] = (-1)**(((pa - ps_obs[k])%4)//2) return xs - + @njit def stabilizer_entropy(gs, mask): '''Entanglement entropy of the stabilizer state in a given region. @@ -814,167 +969,3 @@ def aggregate(data_in, inds, l): for i in range(data_in.shape[0]): data_out[inds[i]] += data_in[i] return data_out -# Aug 8th -@njit -def stabilizer_projection_trace(gs_stb, ps_stb, gs_obs, ps_obs, r): - '''Measure Pauli operators on a stabilizer state. - - Parameters: - gs_stb: int (2*N, 2*N) - Pauli strings in original stabilizer tableau. - ps_stb: int (N) - phase indicators of (de)stabilizers. - gs_obs: int (L, 2*N) - strings of Pauli operators to be measured. - ps_obs: int (L) - phase indicators of Pauli operators to be measured. - r: int - log2 rank of density matrix (num of standby stablizers). - - Returns: - gs_stb: int (2*N, 2*N) - Pauli strings in updated stabilizer tableau. - ps_stb: int (N) - phase indicators of (de)stabilizers. - r: int - updated log2 rank of density matrix. - trace: float - Tr(P * rho1 * P*)''' - (L, Ng) = gs_obs.shape - N = Ng//2 -# assert L==N # Current implementation is full state projection - assert 0<=r<=N - ga = numpy.empty(2*N, dtype=numpy.int_) # workspace for stabilizer accumulation - pa = 0 # workspace for phase accumulation - trace = 1 - for k in range(L): # for each observable gs_obs[k] - update = False - extend = False - p = 0 # pointer - ga[:] = 0 - pa = 0 - for j in range(2*N): - if acq(gs_stb[j], gs_obs[k]): # find gs_stb[j] anticommute with gs_obs[k] - if update: # if gs_stb[j] is not the first anticommuting operator - # update gs_stb[j] to commute with gs_obs[k] - if j < N: # if gs_stb[j] is a stablizer, phase matters - ps_stb[j] = (ps_stb[j] + ps_stb[p] + ipow(gs_stb[j], gs_stb[p]))%4 - gs_stb[j] = (gs_stb[j] + gs_stb[p])%2 - else: # if gs_stb[j] is the first anticommuting operator - if j < N + r: # if gs_stb[j] is not an active destabilizer - p = j # move pointer to j - update = True - if not r <= j < N: # if gs_stb[j] is a standby operator - extend = True - else: # gs_stb[j] anticommute with destabilizer, meaning gs_obs[k] already a combination of active stabilizers - # collect corresponding stabilizer component to ga - pa = (pa + ps_stb[j-N] + ipow(ga, gs_stb[j-N]))%4 - ga = (ga + gs_stb[j-N])%2 - if update: - # now gs_stb[p] and gs_obs[k] anticommute - q = (p+N)%(2*N) # get q as dual of p - gs_stb[q] = gs_stb[p] # move gs_stb[p] to gs_stb[q] - gs_stb[p] = gs_obs[k] # add gs_obs[k] to gs_stb[p] - if extend: - r -= 1 # rank will reduce under extension - # bring new stabilizer from p to r - if p == r: - pass - elif q == r: - gs_stb[numpy.array([p,q])] = gs_stb[numpy.array([q,p])] # swap p,q - else: - s = (r+N)%(2*N) # get s as dual of r - gs_stb[numpy.array([p,r])] = gs_stb[numpy.array([r,p])] # swap p,r - gs_stb[numpy.array([q,s])] = gs_stb[numpy.array([s,q])] # swap q,s - p = r - # the projection will change phase of stabilizer - ps_stb[p] = ps_obs[k] - trace = trace/2. - else: # no update, gs_obs[k] is eigen, result is in pa - assert((ga == gs_obs[k]).all()) - if not pa == ps_obs[k]: - trace = 0. - return gs_stb, ps_stb, r, trace - - - -############ Jan 26 Added by Hong-Ye Hu ############# -@njit -def decompose(g, gs, ps): - ''' Decompose a pauli string to phase*destabilizers*stabilizers - Parameters: - g: int(2*N) - the binary vector of a pauli string - gs: int(2*N,2*N) - the full tableau - ps: int(2*N) - phase of stabilizer and destabilizer - - Returns: - phase: int - phase in terms of imaginery power - b: int(N) - binary encoding of destabilizer decomposition - c: int(N) - binary encoding of stabilizer decomposition - ''' - phase = 0 - tmp_p = np.zeros_like(g) - N = gs.shape[0]//2 - # b = int(np.zeros(N)) # numbda does not support change data type - # c = int(np.zeros(N)) - b = np.array([0 for _ in range(N)]) - c = np.array([0 for _ in range(N)]) - for i in range(N): - if acq(g,gs[i]): #anti-commute - b[i] = 1 - phase = phase - ipow(tmp_p,gs[i+N]) + ps[i+N] - tmp_p = (tmp_p+gs[i+N])%2 - for i in range(N): - if acq(g,gs[i+N]): #anti-commute - c[i] = 1 - phase = phase - ipow(tmp_p,gs[i]) + ps[i] - tmp_p = (tmp_p+gs[i])%2 - return phase%4, tmp_p, b, c - -@njit -def stabilizer_postselection(gs_stb, ps_stb, gs_ob, ps_ob): - ''' - Stabilizer post-selection (pure state) - - Parameters: - gs_stb: int (2*N, 2*N) - Pauli strings in original stabilizer tableau. - ps_stb: int (N) - phase indicators of (de)stabilizers. - gs_ob: int (2*N) - strings of Pauli operators to be measured. - ps_obs: int (1) - phase indicators of Pauli operators to be measured. - Returns: - gs_stb: int (2*N, 2*N) - Pauli strings in updated stabilizer tableau. - ps_stb: int (N) - phase indicators of (de)stabilizers. - prob: float - ''' - (_, Ng) = gs_stb.shape - N = Ng//2 - ga = numpy.empty(2*N, dtype=numpy.int_) # workspace for stabilizer accumulation - pa = 0 # workspace for phase accumulation - prob = 1.0 - # one projection by gs_ob - update = False - p=0 - ga[:] = 0 - pa = 0 - for j in range(2*N): - if acq(gs_stb[j], gs_ob): # find gs_stb[j] anticommute with gs_obs[k] - if update: # if gs_stb[j] is not the first anticommuting operator - # update gs_stb[j] to commute with gs_obs[k] - if j < N: # if gs_stb[j] is a stablizer, phase matters - ps_stb[j] = (ps_stb[j] + ps_stb[p] + ipow(gs_stb[j], gs_stb[p]))%4 - gs_stb[j] = (gs_stb[j] + gs_stb[p])%2 - else: # if gs_stb[j] is the first anticommuting operator - if j < N: # if gs_stb[j] is not an active destabilizer - p = j # move pointer to j - update = True - - else: # gs_stb[j] anticommute with destabilizer, meaning gs_obs[k] already a combination of active stabilizers - # collect corresponding stabilizer component to ga - pa = (pa + ps_stb[j-N] + ipow(ga, gs_stb[j-N]))%4 - ga = (ga + gs_stb[j-N])%2 - if update: - # now gs_stb[p] and gs_obs[k] anticommute - q = (p+N)%(2*N) # get q as dual of p - gs_stb[q] = gs_stb[p] # move gs_stb[p] to gs_stb[q] - gs_stb[p] = gs_ob # add gs_obs[k] to gs_stb[p] - - # the projection will change phase of stabilizer - ps_stb[p] = ps_ob - prob = prob/2.0 - - else: # no update, gs_obs[k] is eigen, result is in pa - assert((ga == gs_ob).all()) - if not pa == ps_ob: - prob = 0. - return gs_stb, ps_stb, prob \ No newline at end of file From d6a8179a96915ddbc7ca944ad39209f7fd2b219a Mon Sep 17 00:00:00 2001 From: Everett Date: Fri, 9 May 2025 01:03:57 -0700 Subject: [PATCH 2/5] generalized stabilizer state --- dev/GeneralStabilizerState.ipynb | 2 +- dev/Untitled.ipynb | 33 -- dev/dev_generalized_stabilizer.ipynb | 380 ++++++++++++++++++ pyclifford/pyclifford.egg-info/PKG-INFO | 12 + pyclifford/pyclifford.egg-info/SOURCES.txt | 9 + .../pyclifford.egg-info/dependency_links.txt | 1 + pyclifford/pyclifford.egg-info/requires.txt | 2 + pyclifford/pyclifford.egg-info/top_level.txt | 1 + pyclifford/stabilizer.py | 2 +- pyclifford/utils.py | 50 ++- 10 files changed, 456 insertions(+), 36 deletions(-) delete mode 100644 dev/Untitled.ipynb create mode 100644 dev/dev_generalized_stabilizer.ipynb create mode 100644 pyclifford/pyclifford.egg-info/PKG-INFO create mode 100644 pyclifford/pyclifford.egg-info/SOURCES.txt create mode 100644 pyclifford/pyclifford.egg-info/dependency_links.txt create mode 100644 pyclifford/pyclifford.egg-info/requires.txt create mode 100644 pyclifford/pyclifford.egg-info/top_level.txt diff --git a/dev/GeneralStabilizerState.ipynb b/dev/GeneralStabilizerState.ipynb index b3f65ca..f5afcd9 100644 --- a/dev/GeneralStabilizerState.ipynb +++ b/dev/GeneralStabilizerState.ipynb @@ -178,7 +178,7 @@ "metadata": {}, "source": [ "From the result, we see pauli string $XY$ can be decomposed as\n", - "$$XY = i^1 d_0^0 d_1^1 g_0^1 g_1^1=i(XI)(-ZZ)(XX)$$" + "$$XY = i^1 d_0^0 d_1^1 g_0^1 g_1^1=i(ZI)(-ZZ)(XX)$$" ] }, { diff --git a/dev/Untitled.ipynb b/dev/Untitled.ipynb deleted file mode 100644 index 582e3b8..0000000 --- a/dev/Untitled.ipynb +++ /dev/null @@ -1,33 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "8334d599-4421-4ee6-81f6-58b3a2767b46", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.12" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/dev/dev_generalized_stabilizer.ipynb b/dev/dev_generalized_stabilizer.ipynb new file mode 100644 index 0000000..72c590a --- /dev/null +++ b/dev/dev_generalized_stabilizer.ipynb @@ -0,0 +1,380 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 104, + "id": "8334d599-4421-4ee6-81f6-58b3a2767b46", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "import numpy\n", + "import matplotlib.pyplot as plt\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "%matplotlib inline\n", + "# will autoupdate any of the packages imported:\n", + "import pyclifford as pc" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "id": "9033b484", + "metadata": {}, + "outputs": [], + "source": [ + "from pyclifford.utils import pauli_decompose, calculate_chi, pauli_combine\n", + "from pyclifford.paulialg import PauliList, Pauli, PauliPolynomial" + ] + }, + { + "cell_type": "code", + "execution_count": 106, + "id": "e7395e1b", + "metadata": {}, + "outputs": [], + "source": [ + "class PauliChannel(object):\n", + " '''Pauli channel.\n", + "\n", + " Parameters:\n", + " paulis: PauliList - a list of Pauli operators serving as operator basis.\n", + " phi: complex (L, L) - channel density matrix in the Pauli basis.'''\n", + " def __init__(self, paulis, phi):\n", + " self.paulis = paulis\n", + " self.phi = phi\n", + " self.N = self.paulis.N" + ] + }, + { + "cell_type": "code", + "execution_count": 120, + "id": "1d96ab45", + "metadata": {}, + "outputs": [], + "source": [ + "class GeneralizedStabilizerState(object):\n", + " '''Generalized stabilizer state.\n", + " rho = sum_{b,b'} chi_{b,b'} |b> is a basis state of destabilizer excitations,\n", + " |b> = prod_{i} D_i^{b_i} |0>,\n", + " with |0> being the state stabilzed by all stabilizers, i.e.\n", + " S_i |0> = |0>.\n", + " - The stabilizers and destabilizers are given by a StabilizerState \n", + " instance as the stabilizer frame.\n", + " - The basis of destabilizer excitations are represented as\n", + " binary array (tuple), e.g. (0,1,0,1) means the 2nd and 4th \n", + " destabilizers are excited.\n", + " - The basis dictionary keeps track of the mapping from unique\n", + " binary tuple to integer index of the basis.\n", + " - The density matrix is given by a complex array chi of size L x L,\n", + " where L is the number of excitation basis states.\n", + "\n", + " Parameters:\n", + " frame: StabilizerState - a stabilizer state serving as the frame.\n", + " basis: int (L, N) - a binary array encoding the basis of destabilizer excitations.\n", + " chi: complex (L, L) - density matrix in the excitation basis.'''\n", + " def __init__(self, frame, basis, chi):\n", + " self.frame = frame\n", + " self.basis = basis\n", + " self.chi = chi\n", + " self.N = self.frame.N\n", + "\n", + " def __repr__(self):\n", + " return f\"GeneralizedStabilizerState(\\nframe=\\n{self.frame},\\nbasis=\\n{self.basis},\\nchi=\\n{self.chi})\"\n", + " \n", + " def copy(self):\n", + " return GeneralizedStabilizerState(self.frame.copy(), self.basis.copy(), self.chi.copy())\n", + " \n", + " def rotate_by(self, generator, mask=None):\n", + " self.frame.rotate_by(generator, mask)\n", + " return self\n", + " \n", + " def transform_by(self, clifford_map, mask=None):\n", + " self.frame.transform_by(clifford_map, mask)\n", + " return self\n", + " \n", + " def evolve_by(self, pauli_channel):\n", + " # Evolve the state by a Pauli channel.\n", + " bs, cs, ps = pauli_decompose(pauli_channel.paulis.gs, pauli_channel.paulis.ps, \n", + " self.frame.gs, self.frame.ps)\n", + " # construct new basis, compute fusion map and fusion phase indicator\n", + " L_old = self.basis.shape[0]\n", + " L_add = bs.shape[0]\n", + " L_new = 0\n", + " basis_new = {}\n", + " fusion_map = numpy.zeros((L_old,L_add), dtype=numpy.int_)\n", + " fusion_p = numpy.zeros((L_old,L_add), dtype=numpy.int_)\n", + " for i in range(L_old):\n", + " for j in range(L_add):\n", + " b_new = tuple((self.basis[i] + bs[j])%2)\n", + " if b_new not in basis_new:\n", + " basis_new[b_new] = L_new\n", + " L_new += 1\n", + " k = basis_new[b_new]\n", + " fusion_map[i,j] = k\n", + " fusion_p[i,j] = ps[j] + 2*self.basis[i].dot(cs[j])\n", + " # perform fusion of state and channel density matrices\n", + " chi_new = calculate_chi(self.chi, pauli_channel.phi, fusion_map, fusion_p, L_new)\n", + " # in-place update basis and chi\n", + " self.basis = numpy.array(list(basis_new.keys()))\n", + " self.chi = chi_new\n", + " return self\n", + " \n", + " def represent(self, ops):\n", + " '''Represent a list of operators in the destabilizer excitation basis.\n", + " \n", + " Parameters:\n", + " ops: operator(s), can be Pauli, PauliList, PauliPolynomial'''\n", + " if isinstance(ops, Pauli):\n", + " # cast Pauli to PauliPolynomial\n", + " return self.represent(ops.as_polynomial())\n", + " if isinstance(ops, PauliPolynomial):\n", + " # cast PauliPolynomial to PauliList\n", + " mats = self.represent(PauliList(ops.gs,ops.ps))\n", + " return numpy.tensordot(ops.cs, mats, axes=(0,0))\n", + " if isinstance(ops, PauliList):\n", + " bs, cs, ps = pauli_decompose(ops.gs, ops.ps, self.frame.gs, self.frame.ps)\n", + " L = self.basis.shape[0]\n", + " K = ps.shape[0]\n", + " mats = numpy.zeros((K,L,L), dtype=numpy.complex_)\n", + " idx = {tuple(b): i for i, b in enumerate(self.basis)}\n", + " for k in range(K):\n", + " for i0 in range(L):\n", + " b0 = self.basis[i0]\n", + " b1 = tuple((b0 + bs[k])%2)\n", + " if b1 in idx:\n", + " i1 = idx[b1]\n", + " mats[k,i1,i0] = 1j**(ps[k] + 2*b0.dot(cs[k]))\n", + " return mats\n", + " \n", + " def expect(self, obs):\n", + " '''Evaluate expctationvalues of Pauli observables on the generalized\n", + " stabilizer state.\n", + " \n", + " Parameters:\n", + " obs: observable, can be Pauli, PauliList, PauliPolynomial\n", + "\n", + " Returns:\n", + " out: output (depending on the type of obs)\n", + " * Pauli: promote to PauliPolynomial\n", + " * PauliPolynomial O: Tr(rho O)\n", + " * PauliList [O_i]: [Tr(rho O_i)]''' \n", + " mats = self.represent(obs)\n", + " if mats.ndim == 2: # single matrix case\n", + " return numpy.trace(self.chi @ mats)\n", + " else: # batch case\n", + " return numpy.trace(numpy.matmul(self.chi, mats), axis1=1, axis2=2)\n", + " \n", + " def to_numpy(self):\n", + " '''Convert generalized stabilizer state to numpy density matrix representation.'''\n", + " destabilizers = self.frame[self.N:2*self.N]\n", + " gs, ps =pauli_combine(self.basis, destabilizers.gs, destabilizers.ps)\n", + " Ds = PauliList(gs,ps).to_numpy()\n", + " rho0 = self.frame.to_numpy()\n", + " L = Ds.shape[0]\n", + " rho = numpy.zeros_like(rho0)\n", + " for j1 in range(L):\n", + " for j2 in range(L):\n", + " rho += self.chi[j1, j2]* Ds[j1] @ rho0 @ Ds[j2]\n", + " return rho\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 130, + "id": "3e63172d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "GeneralizedStabilizerState(\n", + "frame=\n", + "StabilizerState(\n", + " +XX\n", + " +ZZ),\n", + "basis=\n", + "[[0 0]\n", + " [1 0]],\n", + "chi=\n", + "[[0.8+0.j 0.1+0.j]\n", + " [0.1+0.j 0.2+0.j]])" + ] + }, + "execution_count": 130, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state=GeneralizedStabilizerState(\n", + " pc.stabilizer_state(\"XX\",\"ZZ\"), \n", + " numpy.array([[0,0],[1,0]]), \n", + " numpy.array([[0.8,0.1],[0.1,0.2]],dtype=numpy.complex_))\n", + "state" + ] + }, + { + "cell_type": "code", + "execution_count": 131, + "id": "be822ce2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0.6+0.j, 0. +0.j, 0. +0.j, 0.3+0.j],\n", + " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", + " [0.3+0.j, 0. +0.j, 0. +0.j, 0.4+0.j]])" + ] + }, + "execution_count": 131, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state.to_numpy()" + ] + }, + { + "cell_type": "code", + "execution_count": 134, + "id": "9fb4ecae", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.3+0j)" + ] + }, + "execution_count": 134, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state.expect((pc.pauli(\"XX\")+1j*pc.pauli(\"XY\")+1j*pc.pauli(\"YX\")-pc.pauli(\"YY\"))/4)" + ] + }, + { + "cell_type": "code", + "execution_count": 135, + "id": "7155800c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.30000000000000004+0j)" + ] + }, + "execution_count": 135, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "numpy.trace(state.to_numpy() @ ((pc.pauli(\"XX\")+1j*pc.pauli(\"XY\")+1j*pc.pauli(\"YX\")-pc.pauli(\"YY\"))/4).to_numpy())" + ] + }, + { + "cell_type": "code", + "execution_count": 136, + "id": "9dcaff6b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "GeneralizedStabilizerState(\n", + "frame=\n", + "StabilizerState(\n", + " +XX\n", + " +ZZ),\n", + "basis=\n", + "[[0 0]\n", + " [1 0]],\n", + "chi=\n", + "[[0.5+0.j 0.1+0.j]\n", + " [0.1+0.j 0.5+0.j]])" + ] + }, + "execution_count": 136, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "channel = PauliChannel(pc.paulis(\"II\",\"ZI\"), numpy.array([[0.5,0.0],[0.0,0.5]],dtype=numpy.complex_))\n", + "state.evolve_by(channel)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 137, + "id": "3a6649ac", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0.6+0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", + " [0. +0.j, 0. +0.j, 0. +0.j, 0.4+0.j]])" + ] + }, + "execution_count": 137, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state.to_numpy()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c16191ca", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyclifford/pyclifford.egg-info/PKG-INFO b/pyclifford/pyclifford.egg-info/PKG-INFO new file mode 100644 index 0000000..2800d9d --- /dev/null +++ b/pyclifford/pyclifford.egg-info/PKG-INFO @@ -0,0 +1,12 @@ +Metadata-Version: 2.4 +Name: pyclifford +Version: 0.1.0 +Summary: A Python package for Clifford algebra and quantum computing +Author: PyClifford Team +Requires-Python: >=3.6 +Requires-Dist: numpy +Requires-Dist: matplotlib +Dynamic: author +Dynamic: requires-dist +Dynamic: requires-python +Dynamic: summary diff --git a/pyclifford/pyclifford.egg-info/SOURCES.txt b/pyclifford/pyclifford.egg-info/SOURCES.txt new file mode 100644 index 0000000..94df70c --- /dev/null +++ b/pyclifford/pyclifford.egg-info/SOURCES.txt @@ -0,0 +1,9 @@ +setup.py +pyclifford.egg-info/PKG-INFO +pyclifford.egg-info/SOURCES.txt +pyclifford.egg-info/dependency_links.txt +pyclifford.egg-info/requires.txt +pyclifford.egg-info/top_level.txt +tests/__init__.py +tests/test_paulialg.py +tests/test_stabilizer.py \ No newline at end of file diff --git a/pyclifford/pyclifford.egg-info/dependency_links.txt b/pyclifford/pyclifford.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/pyclifford/pyclifford.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/pyclifford/pyclifford.egg-info/requires.txt b/pyclifford/pyclifford.egg-info/requires.txt new file mode 100644 index 0000000..aa094d9 --- /dev/null +++ b/pyclifford/pyclifford.egg-info/requires.txt @@ -0,0 +1,2 @@ +numpy +matplotlib diff --git a/pyclifford/pyclifford.egg-info/top_level.txt b/pyclifford/pyclifford.egg-info/top_level.txt new file mode 100644 index 0000000..2b29f27 --- /dev/null +++ b/pyclifford/pyclifford.egg-info/top_level.txt @@ -0,0 +1 @@ +tests diff --git a/pyclifford/stabilizer.py b/pyclifford/stabilizer.py index 673fa47..8ac4e59 100644 --- a/pyclifford/stabilizer.py +++ b/pyclifford/stabilizer.py @@ -65,7 +65,7 @@ def inverse(self): '''Returns the inverse of this Clifford map, (such that it composes with its inverse results in identity map).''' gs_inv = z2inv(self.gs) - gs_iden, ps_mis = pauli_combine(gs_inv, self.gs, self.ps) + _, ps_mis = pauli_combine(gs_inv, self.gs, self.ps) ps_inv = (- ps_mis - ps0(gs_inv))%4 return CliffordMap(gs_inv, ps_inv) diff --git a/pyclifford/utils.py b/pyclifford/utils.py index 73b9194..c5683d8 100644 --- a/pyclifford/utils.py +++ b/pyclifford/utils.py @@ -202,7 +202,7 @@ def pauli_tokenize(gs, ps): ts[j,N] = 4 + x * (11 - 9 * x + 2 * x**2) // 2 return ts -# ---- combination and trasnformation ---- +# ---- combination, trasnformation, decomposition ---- @njit def pauli_combine(C, gs_in, ps_in): '''Combine Pauli operators by operator product. @@ -246,6 +246,40 @@ def pauli_transform(gs_in, ps_in, gs_map, ps_map): ps_out = (ps_in + ps0(gs_in) + ps_out)%4 return gs_out, ps_out +@njit +def pauli_decompose(gs_in, ps_in, gs_stb, ps_stb): + '''Decompose a Pauli operator into stabilizer and destabilizers. + + Parameters: + gs_in: int (L, 2*N) - Pauli strings in binary representation. + ps_in: int (L) - phase indicators of Pauli operators. + gs_stb: int (2*N, 2*N) - stabilizer tableaux in binary representation. + ps_stb: int (2*N) - phase indicators of (de)stabilizer. + + Returns: + bs_out: int (L, N) - binary encoding of destabilizer decomposition. + cs_out: int (L, N) - binary encoding of stabilizer decomposition. + ps_out: int (L) - phase indicators of decomposed operators.''' + (L, N2) = gs_in.shape + N = N2//2 + bs_out = numpy.zeros((L, N), dtype=numpy.int_) + cs_out = numpy.zeros((L, N), dtype=numpy.int_) + ps_out = ps_in.copy() + g_tmp = numpy.zeros(N2, dtype=numpy.int_) + for k in range(L): + g_tmp.fill(0) + for j in range(N): + if acq(gs_in[k], gs_stb[j]): + bs_out[k,j] = 1 + ps_out[k] = ps_out[k] - ps_stb[j+N] - ipow(g_tmp, gs_stb[j+N]) + g_tmp = (g_tmp + gs_stb[j+N])%2 + for j in range(N): + if acq(gs_in[k], gs_stb[j+N]): + cs_out[k,j] = 1 + ps_out[k] = ps_out[k] - ps_stb[j] - ipow(g_tmp, gs_stb[j]) + g_tmp = (g_tmp + gs_stb[j])%2 + return bs_out, cs_out, ps_out%4 + # ---- clifford rotation ---- @njit def clifford_rotate(g, p, gs, ps): @@ -969,3 +1003,17 @@ def aggregate(data_in, inds, l): for i in range(data_in.shape[0]): data_out[inds[i]] += data_in[i] return data_out + +# ---- generalized stabilizer utilities ---- +@njit +def calculate_chi(chi_old, phi, fusion_map, fusion_p, L_new): + L_old, L_add = fusion_map.shape + chi_new = numpy.zeros((L_new,L_new), dtype=numpy.complex_) + for i1 in range(L_old): + for j1 in range(L_add): + k1 = fusion_map[i1,j1] + for i2 in range(L_old): + for j2 in range(L_add): + k2 = fusion_map[i2,j2] + chi_new[k1,k2] += chi_old[i1,i2] * phi[j1,j2] * 1j**(fusion_p[i1,j1] + fusion_p[i2,j2]) + return chi_new \ No newline at end of file From b410dcea135296b1a9ba388d925ac7c541afcbf3 Mon Sep 17 00:00:00 2001 From: Everett Date: Mon, 12 May 2025 18:32:45 -0700 Subject: [PATCH 3/5] update --- dev/dev_generalized_stabilizer.ipynb | 4 +- dev/test_agent_benchmark.ipynb | 407 +++++++++++++++++++++++++++ 2 files changed, 409 insertions(+), 2 deletions(-) create mode 100644 dev/test_agent_benchmark.ipynb diff --git a/dev/dev_generalized_stabilizer.ipynb b/dev/dev_generalized_stabilizer.ipynb index 72c590a..50362f6 100644 --- a/dev/dev_generalized_stabilizer.ipynb +++ b/dev/dev_generalized_stabilizer.ipynb @@ -293,7 +293,7 @@ }, { "cell_type": "code", - "execution_count": 136, + "execution_count": null, "id": "9dcaff6b", "metadata": {}, "outputs": [ @@ -320,7 +320,7 @@ ], "source": [ "channel = PauliChannel(pc.paulis(\"II\",\"ZI\"), numpy.array([[0.5,0.0],[0.0,0.5]],dtype=numpy.complex_))\n", - "state.evolve_by(channel)\n" + "state.evolve_by(channel)" ] }, { diff --git a/dev/test_agent_benchmark.ipynb b/dev/test_agent_benchmark.ipynb new file mode 100644 index 0000000..b8b028f --- /dev/null +++ b/dev/test_agent_benchmark.ipynb @@ -0,0 +1,407 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# add the parent directory to the path, so that we can import the package\n", + "import sys\n", + "sys.path.insert(0, '..')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "import numpy\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "# will autoupdate any of the packages imported:\n", + "import pyclifford as pc" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Helper Functions" + ] + }, + { + "cell_type": "code", + "execution_count": 63, + "metadata": {}, + "outputs": [], + "source": [ + "def to_latex(obj):\n", + " if isinstance(obj, pc.PauliMonomial):\n", + " z = obj.c * (1j ** obj.p)\n", + " if z == 0.:\n", + " return \"0\"\n", + " elif z.real == 0.:\n", + " s = \"\"\n", + " elif z == 1.:\n", + " s = \"\"\n", + " elif z == -1.:\n", + " s = \"-\"\n", + " else:\n", + " s = str(z.real)\n", + " if z.imag > 0:\n", + " s += \"+\"\n", + " if z.imag != 1.:\n", + " s += str(z.imag)\n", + " s += \"i\"\n", + " elif z.imag < 0:\n", + " if z.imag != -1.:\n", + " s += str(z.imag)\n", + " else:\n", + " s += \"-\"\n", + " s += \"i\"\n", + " if len(s) > 0:\n", + " s += \" \"\n", + " for i in range(obj.N):\n", + " if numpy.array_equal(obj.g[2*i:2*i+2], [1,0]):\n", + " s += \"X\" + \"_{\" + str(i+1) + \"}\"\n", + " elif numpy.array_equal(obj.g[2*i:2*i+2], [1,1]):\n", + " s += \"Y\" + \"_{\" + str(i+1) + \"}\"\n", + " elif numpy.array_equal(obj.g[2*i:2*i+2], [0,1]):\n", + " s += \"Z\" + \"_{\" + str(i+1) + \"}\"\n", + " if s[-1] == \" \":\n", + " s += \"I\"\n", + " return s\n", + " elif isinstance(obj, pc.PauliPolynomial):\n", + " s = \"\"\n", + " for k, term in enumerate(obj):\n", + " s_term = to_latex(term)\n", + " if k != 0:\n", + " if s_term[0] == '-':\n", + " s_term = ' ' + s_term\n", + " else:\n", + " s_term = ' +' + s_term\n", + " s = s + s_term\n", + " return s\n", + " elif isinstance(obj, pc.Pauli):\n", + " return to_latex(obj.as_monomial())\n", + " elif isinstance(obj, pc.PauliList):\n", + " return [to_latex(op) for op in obj]\n", + " else:\n", + " return str(obj)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Benchmark Problems for Agents" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## I. Pauli Algebra" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### I.1. Product of Pauli Strings" + ] + }, + { + "cell_type": "code", + "execution_count": 64, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(+i Z_{1}) (-i X_{1}) = +i Y_{1}\n", + "(- Z_{1}) (-i I) = +i Z_{1}\n", + "(- Y_{1}Y_{2}) (+i X_{1}X_{2}) = +i Z_{1}Z_{2}\n", + "(+i Z_{1}Y_{2}) (Z_{1}X_{2}) = Z_{2}\n", + "(Z_{1}X_{2}Y_{3}) (-i Z_{1}Z_{2}X_{3}) = +i Y_{2}Z_{3}\n", + "(+i Z_{1}Y_{2}X_{3}) (-i X_{2}Y_{3}Y_{4}) = Z_{1}Z_{2}Z_{3}Y_{4}\n", + "(- X_{2}X_{4}) (- X_{1}Y_{3}X_{4}Z_{5}) = X_{1}X_{2}Y_{3}Z_{5}\n", + "(+i Y_{2}Z_{3}X_{4}X_{5}Z_{6}) (-i X_{1}X_{2}X_{3}Y_{4}Y_{5}Z_{6}) = - X_{1}Z_{2}Y_{3}Z_{4}Z_{5}\n", + "(-i X_{1}Y_{2}Z_{5}X_{6}Y_{7}Z_{8}) (-i Y_{1}Y_{2}Z_{3}X_{4}X_{5}X_{7}) = -i Z_{1}Z_{3}X_{4}Y_{5}X_{6}Z_{7}Z_{8}\n", + "(- X_{2}X_{3}X_{4}Z_{5}Z_{6}Y_{7}Y_{8}Z_{10}) (- Y_{3}X_{4}X_{5}Y_{6}X_{7}Z_{8}) = +i X_{2}Z_{3}Y_{5}X_{6}Z_{7}X_{8}Z_{10}\n", + "(-i X_{1}X_{2}Z_{4}Z_{5}Z_{6}Y_{8}Z_{9}Z_{10}X_{11}X_{12}Z_{14}) (X_{2}X_{3}Y_{4}X_{5}Z_{6}Z_{7}Z_{8}X_{9}Z_{10}Z_{11}Z_{12}X_{13}X_{14}) = X_{1}X_{3}X_{4}Y_{5}Z_{7}X_{8}Y_{9}Y_{11}Y_{12}X_{13}Y_{14}\n", + "(- X_{3}X_{4}Z_{5}X_{6}Z_{7}Y_{8}Y_{10}Z_{12}Z_{14}X_{17}Z_{18}) (-i Z_{1}X_{2}X_{3}X_{4}Y_{5}X_{6}Z_{8}Z_{9}X_{10}Z_{11}Y_{13}Z_{15}Y_{16}Y_{17}) = +i Z_{1}X_{2}X_{5}Z_{7}X_{8}Z_{9}Z_{10}Z_{11}Z_{12}Y_{13}Z_{14}Z_{15}Y_{16}Z_{17}Z_{18}\n", + "\n" + ] + } + ], + "source": [ + "def prod_pauli(Nmin, samples, seed=42):\n", + " rng = numpy.random.default_rng(seed)\n", + " txt = \"\"\n", + " for k in range(samples):\n", + " N = round(Nmin + 1.3**k - 1)\n", + " gs = rng.integers(0, 2, (2, 2*N))\n", + " ps = rng.integers(0, 4, (2,))\n", + " op0, op1 = pc.PauliList(gs, ps)\n", + " op = op0 @ op1\n", + " txt += \"({}) ({}) = {}\\n\".format(to_latex(op0), to_latex(op1), to_latex(op))\n", + " return txt\n", + "print(prod_pauli(1, 12))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### I.2. Krylov Space Construction" + ] + }, + { + "cell_type": "code", + "execution_count": 178, + "metadata": {}, + "outputs": [], + "source": [ + "def krylov_space(H, tol=1e-10):\n", + " basis = [pc.pauli_identity(H.N)]\n", + " def inner(A, B):\n", + " return (A @ B).trace()/2**H.N\n", + " while True:\n", + " # Generate next potential basis vector H^n\n", + " Knew = H @ basis[-1]\n", + " # Gram-Schmidt orthogonalization against existing basis\n", + " for K in basis:\n", + " proj = inner(K, Knew).item()\n", + " Knew = Knew - proj * K\n", + " # Check if new vector is linearly independent\n", + " norm = inner(Knew, Knew).real.item()\n", + " if norm < tol:\n", + " break\n", + " # Add normalized vector to basis\n", + " basis.append(Knew/numpy.sqrt(norm).item())\n", + " dim = len(basis)\n", + " for i, A in enumerate(basis):\n", + " for j, B in enumerate(basis):\n", + " for k, C in enumerate(basis):\n", + " ope = inner(C, A @ B).real.round(decimals=15)\n", + " if ope != 0.:\n", + " print(\"({},{},{}): {}\".format(i, j, k, ope))\n", + " return basis" + ] + }, + { + "cell_type": "code", + "execution_count": 179, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0,0,0): 1.0\n", + "(0,1,1): 1.0\n", + "(0,2,2): 1.0\n", + "(1,0,1): 1.0\n", + "(1,1,0): 1.0\n", + "(1,1,2): -1.0\n", + "(1,2,1): -1.0\n", + "(2,0,2): 1.0\n", + "(2,1,1): -1.0\n", + "(2,2,0): 1.0\n" + ] + }, + { + "data": { + "text/plain": [ + "[1 II, -0.71 ZZ -0.71 XX, -1 YY]" + ] + }, + "execution_count": 179, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H = - pc.pauli('XX') - pc.pauli('ZZ')\n", + "krylov_space(H)" + ] + }, + { + "cell_type": "code", + "execution_count": 180, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0,0,0): 1.0\n", + "(0,1,1): 1.0\n", + "(0,2,2): 1.0\n", + "(0,3,3): 1.0\n", + "(0,4,4): 1.0\n", + "(1,0,1): 1.0\n", + "(1,1,0): 1.0\n", + "(1,2,3): 0.408248290463863\n", + "(1,3,2): 0.408248290463863\n", + "(2,0,2): 1.0\n", + "(2,1,3): 0.408248290463863\n", + "(2,2,0): 1.0\n", + "(2,2,2): 1.632993161855453\n", + "(2,2,4): 1.0\n", + "(2,3,1): 0.408248290463863\n", + "(2,4,2): 1.0\n", + "(3,0,3): 1.0\n", + "(3,1,2): 0.408248290463863\n", + "(3,2,1): 0.408248290463863\n", + "(3,3,0): 1.0\n", + "(4,0,4): 1.0\n", + "(4,2,2): 1.0\n", + "(4,4,0): 1.0\n" + ] + }, + { + "data": { + "text/plain": [ + "[1 IIII,\n", + " -0.50 IIZZ -0.50 IZZI -0.50 ZZII -0.50 XXXX,\n", + " 0.41 IZIZ +0.41 ZIZI +0.41 ZZZZ -0.41 XXYY -0.41 XYYX -0.41 YYXX,\n", + " -0.50 ZIIZ +0.50 XYXY +0.50 YXYX -0.50 YYYY,\n", + " -1 YXXY]" + ] + }, + "execution_count": 180, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H = - pc.pauli('ZZII') - pc.pauli('IZZI') - pc.pauli('IIZZ') - pc.pauli('XXXX')\n", + "krylov_space(H)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## II. Random Clifford Circuit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### II.1. Random Clifford Unitary Forward" + ] + }, + { + "cell_type": "code", + "execution_count": 225, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "X_{1} &\\mapsto - Y_{1}Y_{2}X_{4}Z_{5}X_{6}Y_{7}Y_{8}Z_{10}\n", + "Z_{1} &\\mapsto Y_{3}Y_{4}Z_{5}Y_{6}Z_{7}Z_{8}Y_{9}Y_{10}\n", + "X_{2} &\\mapsto - Y_{1}Y_{2}Y_{3}Z_{4}Z_{5}Y_{6}Z_{7}Y_{8}X_{9}X_{10}\n", + "Z_{2} &\\mapsto - X_{1}Y_{3}Y_{4}Z_{5}Z_{6}X_{7}Z_{8}X_{9}X_{10}\n", + "X_{3} &\\mapsto - X_{1}X_{3}X_{5}X_{6}Z_{8}Y_{9}X_{10}\n", + "Z_{3} &\\mapsto Z_{2}X_{4}Y_{5}Y_{8}Z_{9}\n", + "X_{4} &\\mapsto Z_{1}Y_{3}Y_{4}X_{5}Y_{6}X_{7}Y_{8}X_{10}\n", + "Z_{4} &\\mapsto X_{2}Z_{3}Z_{4}Z_{7}Y_{8}X_{10}\n", + "X_{5} &\\mapsto X_{1}Z_{3}Z_{5}Z_{6}Y_{7}Y_{9}Z_{10}\n", + "Z_{5} &\\mapsto Z_{2}Z_{3}Z_{4}Y_{5}Z_{6}X_{7}X_{8}Z_{9}Z_{10}\n", + "X_{6} &\\mapsto Y_{1}Y_{2}Z_{3}Z_{4}X_{5}Z_{6}Z_{7}X_{8}X_{10}\n", + "Z_{6} &\\mapsto X_{1}X_{2}X_{4}Y_{5}Y_{9}Y_{10}\n", + "X_{7} &\\mapsto Z_{1}Z_{2}Z_{3}Z_{4}Z_{5}X_{6}Z_{7}Y_{8}Z_{9}Z_{10}\n", + "Z_{7} &\\mapsto X_{1}Z_{2}Z_{3}Z_{4}Y_{6}X_{9}Z_{10}\n", + "X_{8} &\\mapsto - Y_{1}Y_{2}Z_{3}Z_{4}X_{8}X_{9}\n", + "Z_{8} &\\mapsto Z_{1}X_{2}Z_{4}X_{6}Y_{7}X_{8}Y_{9}\n", + "X_{9} &\\mapsto - Z_{1}Z_{3}Z_{4}X_{5}Y_{6}Z_{9}\n", + "Z_{9} &\\mapsto - Y_{1}Y_{4}Z_{5}Z_{7}Y_{8}Z_{10}\n", + "X_{10} &\\mapsto - Z_{2}Y_{3}Y_{5}X_{6}Y_{7}Z_{8}Y_{9}X_{10}\n", + "Z_{10} &\\mapsto X_{1}X_{2}Y_{3}Y_{4}Z_{5}X_{6}X_{7}X_{9}Z_{10}\n" + ] + }, + { + "data": { + "text/plain": [ + "(['-i Z_{1}Z_{3}Z_{5}X_{8}Y_{9}X_{10}',\n", + " 'X_{1}Y_{2}Y_{3}Y_{4}Y_{5}Y_{7}X_{8}X_{9}',\n", + " '+i Y_{2}Y_{3}X_{4}Y_{5}Z_{6}Z_{7}Z_{8}X_{9}',\n", + " '- Z_{3}X_{4}Y_{5}Y_{6}Z_{7}Z_{8}Y_{9}Z_{10}',\n", + " 'Z_{1}Z_{2}X_{4}Y_{5}Y_{6}Y_{7}Y_{8}Y_{9}X_{10}',\n", + " '-i Z_{1}Y_{2}Y_{3}Y_{4}X_{5}X_{6}Y_{7}Z_{8}Z_{9}X_{10}'],\n", + " ['-i Z_{1}X_{2}Z_{3}Y_{4}Z_{5}Y_{6}Z_{7}Y_{9}Z_{10}',\n", + " 'X_{1}Z_{2}Z_{4}Y_{5}Y_{6}Y_{7}Y_{10}',\n", + " '+i X_{2}Z_{3}Z_{4}Z_{6}Y_{7}Y_{8}Y_{9}X_{10}',\n", + " '- Z_{1}Z_{3}X_{4}X_{6}X_{7}Z_{8}X_{9}X_{10}',\n", + " 'Y_{1}Z_{3}Y_{4}X_{6}Z_{7}X_{8}Z_{10}',\n", + " '+i X_{1}X_{4}X_{5}Z_{6}Y_{7}Z_{9}Z_{10}'])" + ] + }, + "execution_count": 225, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def random_clifford_forward(N, nop, seed=42):\n", + " numpy.random.seed(seed)\n", + " gs = numpy.random.randint(0, 2, (nop, 2*N))\n", + " ps = numpy.random.randint(0, 4, (nop,))\n", + " ops = pc.PauliList(gs, ps)\n", + " map = pc.random_clifford_map(N)\n", + " head = sum([[\"X_{\"+str(i+1)+\"} &\\mapsto\", \"Z_{\"+str(i+1)+\"} &\\mapsto\"] for i in range(N)],[])\n", + " for h, t in zip(head, to_latex(map)):\n", + " print(h+\" \"+t+\" \\\\\\\\\")\n", + " return to_latex(ops), to_latex(ops.transform_by(map))\n", + "\n", + "random_clifford_forward(10, 6)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 35aec8f91782a820e3c4e5846f3dc887adeb2f9a Mon Sep 17 00:00:00 2001 From: Everett Date: Thu, 15 May 2025 11:35:22 -0700 Subject: [PATCH 4/5] dev update --- dev/dev_generalized_stabilizer.ipynb | 2 +- dev/test_agent_benchmark.ipynb | 466 +++++++++++++++++++++++++-- 2 files changed, 438 insertions(+), 30 deletions(-) diff --git a/dev/dev_generalized_stabilizer.ipynb b/dev/dev_generalized_stabilizer.ipynb index 50362f6..1fb7657 100644 --- a/dev/dev_generalized_stabilizer.ipynb +++ b/dev/dev_generalized_stabilizer.ipynb @@ -293,7 +293,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 136, "id": "9dcaff6b", "metadata": {}, "outputs": [ diff --git a/dev/test_agent_benchmark.ipynb b/dev/test_agent_benchmark.ipynb index b8b028f..09c4e02 100644 --- a/dev/test_agent_benchmark.ipynb +++ b/dev/test_agent_benchmark.ipynb @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 261, "metadata": {}, "outputs": [ { @@ -309,33 +309,33 @@ }, { "cell_type": "code", - "execution_count": 225, + "execution_count": 266, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "X_{1} &\\mapsto - Y_{1}Y_{2}X_{4}Z_{5}X_{6}Y_{7}Y_{8}Z_{10}\n", - "Z_{1} &\\mapsto Y_{3}Y_{4}Z_{5}Y_{6}Z_{7}Z_{8}Y_{9}Y_{10}\n", - "X_{2} &\\mapsto - Y_{1}Y_{2}Y_{3}Z_{4}Z_{5}Y_{6}Z_{7}Y_{8}X_{9}X_{10}\n", - "Z_{2} &\\mapsto - X_{1}Y_{3}Y_{4}Z_{5}Z_{6}X_{7}Z_{8}X_{9}X_{10}\n", - "X_{3} &\\mapsto - X_{1}X_{3}X_{5}X_{6}Z_{8}Y_{9}X_{10}\n", - "Z_{3} &\\mapsto Z_{2}X_{4}Y_{5}Y_{8}Z_{9}\n", - "X_{4} &\\mapsto Z_{1}Y_{3}Y_{4}X_{5}Y_{6}X_{7}Y_{8}X_{10}\n", - "Z_{4} &\\mapsto X_{2}Z_{3}Z_{4}Z_{7}Y_{8}X_{10}\n", - "X_{5} &\\mapsto X_{1}Z_{3}Z_{5}Z_{6}Y_{7}Y_{9}Z_{10}\n", - "Z_{5} &\\mapsto Z_{2}Z_{3}Z_{4}Y_{5}Z_{6}X_{7}X_{8}Z_{9}Z_{10}\n", - "X_{6} &\\mapsto Y_{1}Y_{2}Z_{3}Z_{4}X_{5}Z_{6}Z_{7}X_{8}X_{10}\n", - "Z_{6} &\\mapsto X_{1}X_{2}X_{4}Y_{5}Y_{9}Y_{10}\n", - "X_{7} &\\mapsto Z_{1}Z_{2}Z_{3}Z_{4}Z_{5}X_{6}Z_{7}Y_{8}Z_{9}Z_{10}\n", - "Z_{7} &\\mapsto X_{1}Z_{2}Z_{3}Z_{4}Y_{6}X_{9}Z_{10}\n", - "X_{8} &\\mapsto - Y_{1}Y_{2}Z_{3}Z_{4}X_{8}X_{9}\n", - "Z_{8} &\\mapsto Z_{1}X_{2}Z_{4}X_{6}Y_{7}X_{8}Y_{9}\n", - "X_{9} &\\mapsto - Z_{1}Z_{3}Z_{4}X_{5}Y_{6}Z_{9}\n", - "Z_{9} &\\mapsto - Y_{1}Y_{4}Z_{5}Z_{7}Y_{8}Z_{10}\n", - "X_{10} &\\mapsto - Z_{2}Y_{3}Y_{5}X_{6}Y_{7}Z_{8}Y_{9}X_{10}\n", - "Z_{10} &\\mapsto X_{1}X_{2}Y_{3}Y_{4}Z_{5}X_{6}X_{7}X_{9}Z_{10}\n" + "X_{1} &\\mapsto - X_{1}Z_{2}Y_{3}Y_{4}X_{5}Z_{6}Y_{7}Y_{8}Y_{9}Y_{10} \\\\\n", + "Z_{1} &\\mapsto X_{1}Z_{3}Y_{4}Z_{5}Y_{6}Y_{7}Y_{8} \\\\\n", + "X_{2} &\\mapsto - X_{1}X_{2}Y_{3}Z_{5}X_{7}X_{9} \\\\\n", + "Z_{2} &\\mapsto - Y_{2}Z_{3}Y_{6}X_{7}X_{8}Z_{9}Y_{10} \\\\\n", + "X_{3} &\\mapsto - Z_{1}Y_{2}Z_{4}Z_{5}Z_{6}Y_{7}X_{8}Z_{9}Y_{10} \\\\\n", + "Z_{3} &\\mapsto Y_{1}X_{2}Z_{4}Z_{5}Z_{6}Y_{7}X_{8}Z_{10} \\\\\n", + "X_{4} &\\mapsto X_{1}Y_{2}Z_{3}Z_{4}X_{5}Y_{9}X_{10} \\\\\n", + "Z_{4} &\\mapsto X_{1}X_{3}Y_{6}Y_{7}Z_{8}X_{10} \\\\\n", + "X_{5} &\\mapsto Y_{1}Y_{3}Y_{5}X_{6}X_{7}Z_{8}X_{9} \\\\\n", + "Z_{5} &\\mapsto Z_{1}Z_{2}Z_{3}X_{4}Z_{5}X_{6}X_{7}Y_{8}Z_{9}Z_{10} \\\\\n", + "X_{6} &\\mapsto X_{1}X_{2}Y_{5}Z_{8}Y_{9}X_{10} \\\\\n", + "Z_{6} &\\mapsto X_{1}X_{2}Z_{4}Z_{5}X_{6}Y_{7}Z_{9}Z_{10} \\\\\n", + "X_{7} &\\mapsto X_{1}Z_{3}Y_{6}Z_{7}Z_{8} \\\\\n", + "Z_{7} &\\mapsto X_{1}Z_{2}X_{3}Y_{6}X_{7}Y_{8}X_{10} \\\\\n", + "X_{8} &\\mapsto - X_{1}Y_{2}Y_{3}X_{4}Y_{5}Y_{6}X_{7}X_{9} \\\\\n", + "Z_{8} &\\mapsto X_{1}X_{2}X_{7}X_{8}X_{9} \\\\\n", + "X_{9} &\\mapsto - X_{2}X_{3}X_{5}Z_{6}Y_{7}X_{8}Y_{9}X_{10} \\\\\n", + "Z_{9} &\\mapsto - Z_{3}Z_{4}X_{5}Z_{6}X_{8}X_{9}Y_{10} \\\\\n", + "X_{10} &\\mapsto - Y_{1}Y_{2}Y_{3}Y_{4}X_{5}Z_{8}Z_{9} \\\\\n", + "Z_{10} &\\mapsto X_{1}X_{3}X_{4}X_{5}X_{8}X_{9}Y_{10} \\\\\n" ] }, { @@ -347,15 +347,15 @@ " '- Z_{3}X_{4}Y_{5}Y_{6}Z_{7}Z_{8}Y_{9}Z_{10}',\n", " 'Z_{1}Z_{2}X_{4}Y_{5}Y_{6}Y_{7}Y_{8}Y_{9}X_{10}',\n", " '-i Z_{1}Y_{2}Y_{3}Y_{4}X_{5}X_{6}Y_{7}Z_{8}Z_{9}X_{10}'],\n", - " ['-i Z_{1}X_{2}Z_{3}Y_{4}Z_{5}Y_{6}Z_{7}Y_{9}Z_{10}',\n", - " 'X_{1}Z_{2}Z_{4}Y_{5}Y_{6}Y_{7}Y_{10}',\n", - " '+i X_{2}Z_{3}Z_{4}Z_{6}Y_{7}Y_{8}Y_{9}X_{10}',\n", - " '- Z_{1}Z_{3}X_{4}X_{6}X_{7}Z_{8}X_{9}X_{10}',\n", - " 'Y_{1}Z_{3}Y_{4}X_{6}Z_{7}X_{8}Z_{10}',\n", - " '+i X_{1}X_{4}X_{5}Z_{6}Y_{7}Z_{9}Z_{10}'])" + " ['-i Z_{1}Z_{2}Y_{3}Y_{6}Y_{7}Y_{8}Y_{9}Z_{10}',\n", + " '- X_{1}Y_{2}X_{3}X_{4}X_{5}Y_{6}X_{7}X_{8}Y_{10}',\n", + " '+i X_{1}Z_{2}Z_{3}X_{4}X_{5}Y_{6}Y_{8}X_{9}Y_{10}',\n", + " 'Z_{1}Z_{2}Z_{5}Y_{7}X_{8}Y_{9}Z_{10}',\n", + " '- Z_{1}Y_{3}Z_{4}X_{5}Z_{6}X_{7}Y_{8}Y_{9}',\n", + " '-i X_{1}Z_{2}X_{3}X_{5}Y_{7}Y_{8}Z_{9}X_{10}'])" ] }, - "execution_count": 225, + "execution_count": 266, "metadata": {}, "output_type": "execute_result" } @@ -375,6 +375,414 @@ "random_clifford_forward(10, 6)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### II.2. Random Clifford Unitary Backward" + ] + }, + { + "cell_type": "code", + "execution_count": 309, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "X_{1} &\\mapsto - Z_{1}Z_{3}X_{6}Y_{7}X_{8}X_{9}Y_{10} \\\\\n", + "Z_{1} &\\mapsto - Y_{1}Y_{2}Y_{3}Y_{5}X_{6}X_{7}Y_{10} \\\\\n", + "X_{2} &\\mapsto Y_{1}Y_{2}Y_{3}Y_{4}Z_{5}Y_{6}X_{7}X_{9}Y_{10} \\\\\n", + "Z_{2} &\\mapsto - Y_{1}X_{4}Y_{5}Z_{6}X_{7}Y_{8}Z_{9}X_{10} \\\\\n", + "X_{3} &\\mapsto - X_{1}Z_{2}X_{4}X_{5}Z_{6}Z_{8}Y_{9} \\\\\n", + "Z_{3} &\\mapsto - Y_{1}Y_{2}Y_{3}X_{5}X_{6}Z_{7}Y_{9} \\\\\n", + "X_{4} &\\mapsto Z_{1}Y_{2}Y_{3}Y_{5}Y_{7}Z_{9} \\\\\n", + "Z_{4} &\\mapsto Y_{1}Z_{2}Z_{3}X_{4}Z_{5}X_{8}X_{9}Z_{10} \\\\\n", + "X_{5} &\\mapsto - Y_{1}Y_{2}X_{3}X_{4}Z_{5}Z_{7}X_{8}Z_{10} \\\\\n", + "Z_{5} &\\mapsto X_{1}Y_{2}Y_{4}Y_{5}X_{7}Z_{9}X_{10} \\\\\n", + "X_{6} &\\mapsto X_{1}Y_{5}Z_{6}Y_{7}X_{8}Y_{9}X_{10} \\\\\n", + "Z_{6} &\\mapsto - Z_{1}X_{2}X_{3}Z_{4}Y_{5}X_{6}Z_{7}Y_{8}Y_{9}Y_{10} \\\\\n", + "X_{7} &\\mapsto X_{1}Y_{2}Z_{3}Z_{5}Z_{6}X_{7}X_{8}Z_{9}Y_{10} \\\\\n", + "Z_{7} &\\mapsto Z_{2}Y_{3}Y_{5}X_{6}Y_{7}Y_{8}X_{9}Y_{10} \\\\\n", + "X_{8} &\\mapsto X_{1}X_{2}Y_{4}Y_{5}Z_{6}Y_{7}Z_{8}Z_{9} \\\\\n", + "Z_{8} &\\mapsto X_{1}X_{2}Z_{3}Z_{5}Z_{7}Y_{9}X_{10} \\\\\n", + "X_{9} &\\mapsto - Z_{2}Y_{3}X_{4}Z_{7}Y_{8}Z_{9}Y_{10} \\\\\n", + "Z_{9} &\\mapsto Z_{1}Z_{2}Y_{3}Y_{4}Z_{5}X_{6}Y_{7}Z_{8} \\\\\n", + "X_{10} &\\mapsto Y_{1}Z_{2}X_{4}Z_{5}Y_{6}Z_{7}Z_{9}Y_{10} \\\\\n", + "Z_{10} &\\mapsto Y_{1}X_{4}Z_{6}Y_{7}Y_{8}Y_{9} \\\\\n" + ] + }, + { + "data": { + "text/plain": [ + "(['+i Y_{1}X_{2}Z_{4}Z_{5}Z_{6}Y_{7}Z_{8}Z_{9}Y_{10}',\n", + " 'X_{1}X_{2}Z_{6}Y_{7}Z_{8}X_{9}X_{10}',\n", + " 'Z_{1}Z_{2}Z_{3}Z_{4}Z_{7}X_{9}X_{10}',\n", + " '+i Y_{1}Z_{2}Z_{3}Y_{4}Y_{5}Y_{6}Y_{8}X_{9}',\n", + " 'Y_{2}Z_{5}X_{7}X_{8}Y_{9}',\n", + " 'Z_{2}Y_{3}X_{4}Z_{5}X_{6}Y_{7}'],\n", + " ['-i Y_{2}Y_{3}Y_{4}Z_{5}Z_{6}X_{7}Z_{8}X_{10}',\n", + " '- X_{1}Z_{2}X_{3}Z_{4}Y_{5}Z_{7}Z_{8}Y_{9}X_{10}',\n", + " '- X_{1}Z_{3}X_{4}Y_{5}Z_{6}Y_{7}Y_{8}X_{9}',\n", + " '-i Y_{1}Y_{2}X_{3}X_{4}Y_{5}Z_{8}Y_{9}Y_{10}',\n", + " 'X_{1}Y_{2}Z_{3}Z_{4}Z_{5}Y_{7}Z_{8}X_{9}Z_{10}',\n", + " '- Y_{1}X_{2}X_{3}X_{4}Z_{5}X_{6}X_{7}X_{8}Y_{9}X_{10}'])" + ] + }, + "execution_count": 309, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def random_clifford_backward(N, nop, seed=32):\n", + " numpy.random.seed(seed)\n", + " gs = numpy.random.randint(0, 2, (nop, 2*N))\n", + " ps = numpy.random.randint(0, 4, (nop,))\n", + " ops = pc.PauliList(gs, ps)\n", + " map = pc.random_clifford_map(N)\n", + " head = sum([[\"X_{\"+str(i+1)+\"} &\\mapsto\", \"Z_{\"+str(i+1)+\"} &\\mapsto\"] for i in range(N)],[])\n", + " for h, t in zip(head, to_latex(map)):\n", + " print(h+\" \"+t+\" \\\\\\\\\")\n", + " return to_latex(ops), to_latex(ops.transform_by(map.inverse()))\n", + "\n", + "random_clifford_backward(10, 6)" + ] + }, + { + "cell_type": "code", + "execution_count": 291, + "metadata": {}, + "outputs": [], + "source": [ + "from numba import njit" + ] + }, + { + "cell_type": "code", + "execution_count": 312, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 1, 0, 0])" + ] + }, + "execution_count": 312, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "@njit\n", + "def f(seed):\n", + " numpy.random.seed(seed)\n", + " return numpy.random.randint(0, 2, 4)\n", + "\n", + "f(42)" + ] + }, + { + "cell_type": "code", + "execution_count": 313, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "('MT19937',\n", + " array([ 42, 3107752595, 1895908407, 3900362577, 3030691166,\n", + " 4081230161, 2732361568, 1361238961, 3961642104, 867618704,\n", + " 2837705690, 3281374275, 3928479052, 3691474744, 3088217429,\n", + " 1769265762, 3769508895, 2731227933, 2930436685, 486258750,\n", + " 1452990090, 3321835500, 3520974945, 2343938241, 928051207,\n", + " 2811458012, 3391994544, 3688461242, 1372039449, 3706424981,\n", + " 1717012300, 1728812672, 1688496645, 1203107765, 1648758310,\n", + " 440890502, 1396092674, 626042708, 3853121610, 669844980,\n", + " 2992565612, 310741647, 3820958101, 3474052697, 305511342,\n", + " 2053450195, 705225224, 3836704087, 3293527636, 1140926340,\n", + " 2738734251, 574359520, 1493564308, 269614846, 427919468,\n", + " 2903547603, 2957214125, 181522756, 4137743374, 2557886044,\n", + " 3399018834, 1348953650, 1575066973, 3837612427, 705360616,\n", + " 4138204617, 1604205300, 1605197804, 590851525, 2371419134,\n", + " 2530821810, 4183626679, 2872056396, 3895467791, 1156426758,\n", + " 184917518, 2502875602, 2730245981, 3251099593, 2228829441,\n", + " 2591075711, 3048691618, 3030004338, 1726207619, 993866654,\n", + " 823585707, 936803789, 3180156728, 1191670842, 348221088,\n", + " 988038522, 3281236861, 1153842962, 4152167900, 98291801,\n", + " 816305276, 575746380, 1719541597, 2584648622, 1791391551,\n", + " 3234806234, 413529090, 219961136, 4180088407, 1135264652,\n", + " 3923811338, 2304598263, 762142228, 1980420688, 1225347938,\n", + " 3657621885, 3762382117, 1157119598, 2556627260, 2276905960,\n", + " 3857700293, 1903185298, 4258743924, 2078637161, 4160077183,\n", + " 3569294948, 2138906140, 1346725611, 1473959117, 2798330104,\n", + " 3785346335, 4103334026, 3448442764, 1142532843, 4278036691,\n", + " 3071994514, 3474299731, 1121195796, 1536841934, 2132070705,\n", + " 1064908919, 2840327803, 992870214, 2041326888, 2906112696,\n", + " 4182466030, 1031463950, 703166484, 854266995, 4157971695,\n", + " 4071962029, 2600094776, 2770410869, 3776335751, 2599879593,\n", + " 2451043853, 2223709058, 2098813464, 4008111478, 2959232195,\n", + " 3072496064, 2498909222, 4020139729, 785990520, 958060279,\n", + " 4183949075, 2392404465, 533774465, 4092066952, 3967420027,\n", + " 1726137853, 2907699474, 3158758391, 1460845905, 1323598137,\n", + " 2446717890, 3004885867, 3447263769, 1378488047, 3172418196,\n", + " 652839901, 1695052769, 226007057, 778836071, 1216725078,\n", + " 655651335, 1850195064, 427367795, 800074262, 2241880422,\n", + " 1713434925, 339981078, 1730571881, 672610244, 1952245009,\n", + " 2729177102, 3516932475, 4032720152, 3177283432, 411893652,\n", + " 2440235559, 3587427933, 43170267, 39225133, 3904203400,\n", + " 1935961247, 3843123487, 1625453782, 1337993374, 2095455879,\n", + " 3402219947, 634671126, 70868861, 3072823841, 851862432,\n", + " 1828056818, 2794213810, 1222863684, 2164539406, 4249334162,\n", + " 1380362252, 1512719097, 2773165233, 4063118969, 3041859837,\n", + " 529421431, 563872464, 2478730478, 3168749051, 4132953373,\n", + " 3922807735, 1124217574, 1970058502, 1744120743, 1906315107,\n", + " 1074758800, 1611130652, 2878846041, 886823888, 1175456250,\n", + " 1669874674, 2428820171, 1044308794, 3841962192, 138850094,\n", + " 1239727126, 1753711876, 2194286827, 872797664, 4276240980,\n", + " 690338888, 4087206238, 2279169960, 1117436170, 3344885072,\n", + " 3127829945, 315537090, 3802787206, 4157203318, 1637047079,\n", + " 3774106877, 3230158646, 1855823338, 1931415993, 667252379,\n", + " 4288528171, 1587598285, 1096793218, 1916566454, 101891899,\n", + " 2354644560, 3351208292, 1467125166, 2177732119, 4122299478,\n", + " 3904084887, 2653591155, 4201043109, 2867379343, 2660555187,\n", + " 3641744616, 4126452939, 326579197, 2697259239, 3365236848,\n", + " 3007834487, 4118919490, 3306741951, 2285455175, 1956645973,\n", + " 1879691841, 891565150, 1843460149, 2013381028, 819311674,\n", + " 123282948, 1436558519, 1154343666, 206804484, 1650349242,\n", + " 2142011886, 304163699, 2608574600, 2500624796, 2996744833,\n", + " 2344192475, 3152512202, 165571606, 691170269, 1806226529,\n", + " 568535825, 1243813863, 3068953841, 3843784723, 1540495237,\n", + " 4246006858, 1303595780, 3288680241, 864868851, 819595545,\n", + " 3230857496, 3574119395, 1545404573, 2970139338, 4292786727,\n", + " 1803072884, 1374565738, 1736333177, 1978645403, 3962597126,\n", + " 1068006206, 3458125500, 168085922, 1597587506, 2052497512,\n", + " 1323596727, 2421372441, 1468386547, 3574947527, 3363915938,\n", + " 860279252, 1309097460, 3065417722, 1490716202, 3476091722,\n", + " 1669402145, 895071221, 1432690175, 3353592973, 149850974,\n", + " 2789493615, 826939483, 666980418, 755367270, 3988951195,\n", + " 21783894, 1924727373, 1699517788, 1152431122, 2593798113,\n", + " 3522529522, 2797535609, 4018366956, 2350035889, 3010507270,\n", + " 2832621820, 627979167, 997422629, 365587204, 2302500352,\n", + " 1720920631, 689999548, 3713985947, 3267499624, 1971264680,\n", + " 1981530399, 1662926921, 1833821660, 1422522022, 3141447769,\n", + " 2727954526, 4172728772, 1787436028, 1902276939, 3145551277,\n", + " 4207627911, 2497093521, 4111966589, 3929089589, 2253454030,\n", + " 1069424637, 2165048659, 2848813944, 2435898022, 2546206777,\n", + " 3864777677, 3107311565, 3776562483, 1040285049, 3171631943,\n", + " 2404677828, 2522848682, 2930777301, 2831905121, 1436989598,\n", + " 602730315, 664177960, 3959954010, 3116042160, 2881899726,\n", + " 233404945, 4058465099, 1781994751, 485046222, 2776777695,\n", + " 432082123, 1989128370, 86344507, 2510576356, 2194076764,\n", + " 1742125237, 3715839140, 895100548, 147445686, 705462897,\n", + " 2245325113, 1052295404, 1956014786, 2916055958, 1829369612,\n", + " 2541711050, 1594343058, 3708804266, 150438233, 323857098,\n", + " 294681952, 783931535, 606075163, 2427042904, 121207604,\n", + " 3943199031, 1196785464, 1818211378, 1788241109, 3138862427,\n", + " 2037307093, 2306750301, 1644605749, 165986111, 542190743,\n", + " 486828112, 1757411662, 894543082, 4108143634, 1232805238,\n", + " 3801632949, 3863166865, 713767006, 2091486427, 3174776264,\n", + " 1157004409, 623072544, 1667151721, 3361539538, 696723008,\n", + " 3247069452, 682044344, 1382136166, 1385645682, 4219951151,\n", + " 2747881261, 2489355869, 786564174, 2040230554, 2967874556,\n", + " 1414286092, 2677969656, 1393412218, 2216095072, 935533444,\n", + " 3662643439, 3285199608, 3103672804, 522796956, 3952383595,\n", + " 1928659176, 3397717710, 4278554051, 1984736931, 3559102926,\n", + " 1878353094, 875578217, 2398931796, 2313634006, 1606027661,\n", + " 2790634022, 2334166559, 1857067101, 666458681, 1626872683,\n", + " 2155121857, 715449823, 1865157100, 2938814835, 4084911240,\n", + " 45488075, 3474982924, 1750873825, 2246019159, 125388929,\n", + " 1110287838, 652200437, 4212247716, 2702974687, 2963764270,\n", + " 208692058, 3170393729, 1378248367, 752591527, 591629541,\n", + " 2253399388, 2402291226, 3089656189, 3202324513, 3818308310,\n", + " 2828131601, 2690672008, 3676629884, 1007739430, 4072247562,\n", + " 3574795162, 518485611, 1889402182, 3687902739, 3410263649,\n", + " 2790674620, 779455241, 3573984673, 3053204735, 4089925351,\n", + " 789980683, 476440431, 3843536868, 2400661309, 3139919094,\n", + " 1643266656, 113318754, 428163528, 2386492935, 3807242009,\n", + " 574560611, 3174039857, 3774465602, 1164640969, 455942925,\n", + " 1374407495, 2562304709, 1024844203, 521375136, 417432138,\n", + " 1203241821, 2900988280, 2841030991, 2301700751, 369508560,\n", + " 2396447808, 1891459643, 4225682708, 3930667846, 1518293357,\n", + " 2697063889, 3113075061, 2411136298, 2836361984, 4105335811,\n", + " 914081338, 2675982621, 1816939127, 1596754123, 1464603632,\n", + " 1598478676, 1318403529, 4016663081, 2106416852, 2757323084,\n", + " 2042842122, 1175184796, 2212339255, 1334626864, 3994484893,\n", + " 3938045599, 2166620630, 3036360431, 397499085, 975931950,\n", + " 1868702836, 3530424696, 3532548823, 2770836469, 3537418693,\n", + " 3344319345, 3208552526, 1771170897, 4097379814, 3761572528,\n", + " 2794194423, 706836738, 2953105956, 3446096217, 220984542,\n", + " 309619699, 223913021, 3985142640, 1757616575, 2582763607,\n", + " 4018329835, 1393278443, 4121569718, 2087146446, 4282833425,\n", + " 807775617, 1396604749, 3571181413, 90301352, 2618014643,\n", + " 2783561793, 1329389532, 836540831, 26719530], dtype=uint32),\n", + " 624,\n", + " 0,\n", + " 0.0)" + ] + }, + "execution_count": 313, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "numpy.random.seed(42)\n", + "numpy.random.get_state()" + ] + }, + { + "cell_type": "code", + "execution_count": 320, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "('MT19937',\n", + " array([ 42, 3107752595, 1895908407, 3900362577, 3030691166,\n", + " 4081230161, 2732361568, 1361238961, 3961642104, 867618704,\n", + " 2837705690, 3281374275, 3928479052, 3691474744, 3088217429,\n", + " 1769265762, 3769508895, 2731227933, 2930436685, 486258750,\n", + " 1452990090, 3321835500, 3520974945, 2343938241, 928051207,\n", + " 2811458012, 3391994544, 3688461242, 1372039449, 3706424981,\n", + " 1717012300, 1728812672, 1688496645, 1203107765, 1648758310,\n", + " 440890502, 1396092674, 626042708, 3853121610, 669844980,\n", + " 2992565612, 310741647, 3820958101, 3474052697, 305511342,\n", + " 2053450195, 705225224, 3836704087, 3293527636, 1140926340,\n", + " 2738734251, 574359520, 1493564308, 269614846, 427919468,\n", + " 2903547603, 2957214125, 181522756, 4137743374, 2557886044,\n", + " 3399018834, 1348953650, 1575066973, 3837612427, 705360616,\n", + " 4138204617, 1604205300, 1605197804, 590851525, 2371419134,\n", + " 2530821810, 4183626679, 2872056396, 3895467791, 1156426758,\n", + " 184917518, 2502875602, 2730245981, 3251099593, 2228829441,\n", + " 2591075711, 3048691618, 3030004338, 1726207619, 993866654,\n", + " 823585707, 936803789, 3180156728, 1191670842, 348221088,\n", + " 988038522, 3281236861, 1153842962, 4152167900, 98291801,\n", + " 816305276, 575746380, 1719541597, 2584648622, 1791391551,\n", + " 3234806234, 413529090, 219961136, 4180088407, 1135264652,\n", + " 3923811338, 2304598263, 762142228, 1980420688, 1225347938,\n", + " 3657621885, 3762382117, 1157119598, 2556627260, 2276905960,\n", + " 3857700293, 1903185298, 4258743924, 2078637161, 4160077183,\n", + " 3569294948, 2138906140, 1346725611, 1473959117, 2798330104,\n", + " 3785346335, 4103334026, 3448442764, 1142532843, 4278036691,\n", + " 3071994514, 3474299731, 1121195796, 1536841934, 2132070705,\n", + " 1064908919, 2840327803, 992870214, 2041326888, 2906112696,\n", + " 4182466030, 1031463950, 703166484, 854266995, 4157971695,\n", + " 4071962029, 2600094776, 2770410869, 3776335751, 2599879593,\n", + " 2451043853, 2223709058, 2098813464, 4008111478, 2959232195,\n", + " 3072496064, 2498909222, 4020139729, 785990520, 958060279,\n", + " 4183949075, 2392404465, 533774465, 4092066952, 3967420027,\n", + " 1726137853, 2907699474, 3158758391, 1460845905, 1323598137,\n", + " 2446717890, 3004885867, 3447263769, 1378488047, 3172418196,\n", + " 652839901, 1695052769, 226007057, 778836071, 1216725078,\n", + " 655651335, 1850195064, 427367795, 800074262, 2241880422,\n", + " 1713434925, 339981078, 1730571881, 672610244, 1952245009,\n", + " 2729177102, 3516932475, 4032720152, 3177283432, 411893652,\n", + " 2440235559, 3587427933, 43170267, 39225133, 3904203400,\n", + " 1935961247, 3843123487, 1625453782, 1337993374, 2095455879,\n", + " 3402219947, 634671126, 70868861, 3072823841, 851862432,\n", + " 1828056818, 2794213810, 1222863684, 2164539406, 4249334162,\n", + " 1380362252, 1512719097, 2773165233, 4063118969, 3041859837,\n", + " 529421431, 563872464, 2478730478, 3168749051, 4132953373,\n", + " 3922807735, 1124217574, 1970058502, 1744120743, 1906315107,\n", + " 1074758800, 1611130652, 2878846041, 886823888, 1175456250,\n", + " 1669874674, 2428820171, 1044308794, 3841962192, 138850094,\n", + " 1239727126, 1753711876, 2194286827, 872797664, 4276240980,\n", + " 690338888, 4087206238, 2279169960, 1117436170, 3344885072,\n", + " 3127829945, 315537090, 3802787206, 4157203318, 1637047079,\n", + " 3774106877, 3230158646, 1855823338, 1931415993, 667252379,\n", + " 4288528171, 1587598285, 1096793218, 1916566454, 101891899,\n", + " 2354644560, 3351208292, 1467125166, 2177732119, 4122299478,\n", + " 3904084887, 2653591155, 4201043109, 2867379343, 2660555187,\n", + " 3641744616, 4126452939, 326579197, 2697259239, 3365236848,\n", + " 3007834487, 4118919490, 3306741951, 2285455175, 1956645973,\n", + " 1879691841, 891565150, 1843460149, 2013381028, 819311674,\n", + " 123282948, 1436558519, 1154343666, 206804484, 1650349242,\n", + " 2142011886, 304163699, 2608574600, 2500624796, 2996744833,\n", + " 2344192475, 3152512202, 165571606, 691170269, 1806226529,\n", + " 568535825, 1243813863, 3068953841, 3843784723, 1540495237,\n", + " 4246006858, 1303595780, 3288680241, 864868851, 819595545,\n", + " 3230857496, 3574119395, 1545404573, 2970139338, 4292786727,\n", + " 1803072884, 1374565738, 1736333177, 1978645403, 3962597126,\n", + " 1068006206, 3458125500, 168085922, 1597587506, 2052497512,\n", + " 1323596727, 2421372441, 1468386547, 3574947527, 3363915938,\n", + " 860279252, 1309097460, 3065417722, 1490716202, 3476091722,\n", + " 1669402145, 895071221, 1432690175, 3353592973, 149850974,\n", + " 2789493615, 826939483, 666980418, 755367270, 3988951195,\n", + " 21783894, 1924727373, 1699517788, 1152431122, 2593798113,\n", + " 3522529522, 2797535609, 4018366956, 2350035889, 3010507270,\n", + " 2832621820, 627979167, 997422629, 365587204, 2302500352,\n", + " 1720920631, 689999548, 3713985947, 3267499624, 1971264680,\n", + " 1981530399, 1662926921, 1833821660, 1422522022, 3141447769,\n", + " 2727954526, 4172728772, 1787436028, 1902276939, 3145551277,\n", + " 4207627911, 2497093521, 4111966589, 3929089589, 2253454030,\n", + " 1069424637, 2165048659, 2848813944, 2435898022, 2546206777,\n", + " 3864777677, 3107311565, 3776562483, 1040285049, 3171631943,\n", + " 2404677828, 2522848682, 2930777301, 2831905121, 1436989598,\n", + " 602730315, 664177960, 3959954010, 3116042160, 2881899726,\n", + " 233404945, 4058465099, 1781994751, 485046222, 2776777695,\n", + " 432082123, 1989128370, 86344507, 2510576356, 2194076764,\n", + " 1742125237, 3715839140, 895100548, 147445686, 705462897,\n", + " 2245325113, 1052295404, 1956014786, 2916055958, 1829369612,\n", + " 2541711050, 1594343058, 3708804266, 150438233, 323857098,\n", + " 294681952, 783931535, 606075163, 2427042904, 121207604,\n", + " 3943199031, 1196785464, 1818211378, 1788241109, 3138862427,\n", + " 2037307093, 2306750301, 1644605749, 165986111, 542190743,\n", + " 486828112, 1757411662, 894543082, 4108143634, 1232805238,\n", + " 3801632949, 3863166865, 713767006, 2091486427, 3174776264,\n", + " 1157004409, 623072544, 1667151721, 3361539538, 696723008,\n", + " 3247069452, 682044344, 1382136166, 1385645682, 4219951151,\n", + " 2747881261, 2489355869, 786564174, 2040230554, 2967874556,\n", + " 1414286092, 2677969656, 1393412218, 2216095072, 935533444,\n", + " 3662643439, 3285199608, 3103672804, 522796956, 3952383595,\n", + " 1928659176, 3397717710, 4278554051, 1984736931, 3559102926,\n", + " 1878353094, 875578217, 2398931796, 2313634006, 1606027661,\n", + " 2790634022, 2334166559, 1857067101, 666458681, 1626872683,\n", + " 2155121857, 715449823, 1865157100, 2938814835, 4084911240,\n", + " 45488075, 3474982924, 1750873825, 2246019159, 125388929,\n", + " 1110287838, 652200437, 4212247716, 2702974687, 2963764270,\n", + " 208692058, 3170393729, 1378248367, 752591527, 591629541,\n", + " 2253399388, 2402291226, 3089656189, 3202324513, 3818308310,\n", + " 2828131601, 2690672008, 3676629884, 1007739430, 4072247562,\n", + " 3574795162, 518485611, 1889402182, 3687902739, 3410263649,\n", + " 2790674620, 779455241, 3573984673, 3053204735, 4089925351,\n", + " 789980683, 476440431, 3843536868, 2400661309, 3139919094,\n", + " 1643266656, 113318754, 428163528, 2386492935, 3807242009,\n", + " 574560611, 3174039857, 3774465602, 1164640969, 455942925,\n", + " 1374407495, 2562304709, 1024844203, 521375136, 417432138,\n", + " 1203241821, 2900988280, 2841030991, 2301700751, 369508560,\n", + " 2396447808, 1891459643, 4225682708, 3930667846, 1518293357,\n", + " 2697063889, 3113075061, 2411136298, 2836361984, 4105335811,\n", + " 914081338, 2675982621, 1816939127, 1596754123, 1464603632,\n", + " 1598478676, 1318403529, 4016663081, 2106416852, 2757323084,\n", + " 2042842122, 1175184796, 2212339255, 1334626864, 3994484893,\n", + " 3938045599, 2166620630, 3036360431, 397499085, 975931950,\n", + " 1868702836, 3530424696, 3532548823, 2770836469, 3537418693,\n", + " 3344319345, 3208552526, 1771170897, 4097379814, 3761572528,\n", + " 2794194423, 706836738, 2953105956, 3446096217, 220984542,\n", + " 309619699, 223913021, 3985142640, 1757616575, 2582763607,\n", + " 4018329835, 1393278443, 4121569718, 2087146446, 4282833425,\n", + " 807775617, 1396604749, 3571181413, 90301352, 2618014643,\n", + " 2783561793, 1329389532, 836540831, 26719530], dtype=uint32),\n", + " 624,\n", + " 0,\n", + " 0.0)" + ] + }, + "execution_count": 320, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "numpy.random.seed(42)\n", + "numpy.random.randint(0, 1, 6)\n", + "numpy.random.get_state()" + ] + }, { "cell_type": "code", "execution_count": null, From 9c2dfb875ca202eff79203ec946d5fb0167df9d0 Mon Sep 17 00:00:00 2001 From: Everett Date: Sat, 17 May 2025 13:37:14 -0700 Subject: [PATCH 5/5] generalized stabilizer update --- dev/dev_generalized_stabilizer.ipynb | 216 ++++++++++++++++++++------- pyclifford/circuit.py | 16 +- pyclifford/paulialg.py | 34 ++++- pyclifford/stabilizer.py | 12 ++ pyclifford/utils.py | 25 ++-- 5 files changed, 221 insertions(+), 82 deletions(-) diff --git a/dev/dev_generalized_stabilizer.ipynb b/dev/dev_generalized_stabilizer.ipynb index 1fb7657..deaa518 100644 --- a/dev/dev_generalized_stabilizer.ipynb +++ b/dev/dev_generalized_stabilizer.ipynb @@ -2,20 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 104, + "execution_count": 1, "id": "8334d599-4421-4ee6-81f6-58b3a2767b46", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The autoreload extension is already loaded. To reload it, use:\n", - " %reload_ext autoreload\n" - ] - } - ], + "outputs": [], "source": [ + "import sys\n", + "sys.path.insert(0, '..')\n", "import numpy\n", "import matplotlib.pyplot as plt\n", "%load_ext autoreload\n", @@ -27,7 +20,7 @@ }, { "cell_type": "code", - "execution_count": 105, + "execution_count": 13, "id": "9033b484", "metadata": {}, "outputs": [], @@ -38,17 +31,21 @@ }, { "cell_type": "code", - "execution_count": 106, + "execution_count": 14, "id": "e7395e1b", "metadata": {}, "outputs": [], "source": [ "class PauliChannel(object):\n", " '''Pauli channel.\n", - "\n", + " C[rho] = sum_{m,n} phi_{mn} P_m rho P_n^H,\n", + " where P_m is the m-th element of paulis.\n", + " - The CPTP condition requires phi^H = phi, Tr phi = 1, and phi >=0,\n", + " as if phi is a density matrix.\n", + " \n", " Parameters:\n", - " paulis: PauliList - a list of Pauli operators serving as operator basis.\n", - " phi: complex (L, L) - channel density matrix in the Pauli basis.'''\n", + " paulis: PauliList - a list of Pauli operators {P_m} serving as operator basis.\n", + " phi: complex (L, L) - channel density matrix phi_{mn} in the Pauli basis.'''\n", " def __init__(self, paulis, phi):\n", " self.paulis = paulis\n", " self.phi = phi\n", @@ -57,7 +54,7 @@ }, { "cell_type": "code", - "execution_count": 120, + "execution_count": 15, "id": "1d96ab45", "metadata": {}, "outputs": [], @@ -81,13 +78,20 @@ "\n", " Parameters:\n", " frame: StabilizerState - a stabilizer state serving as the frame.\n", - " basis: int (L, N) - a binary array encoding the basis of destabilizer excitations.\n", + " basis: int (L, N-r) - a binary array encoding the basis of destabilizer excitations.\n", " chi: complex (L, L) - density matrix in the excitation basis.'''\n", " def __init__(self, frame, basis, chi):\n", " self.frame = frame\n", " self.basis = basis\n", " self.chi = chi\n", - " self.N = self.frame.N\n", + "\n", + " @property\n", + " def N(self):\n", + " return self.frame.N\n", + " \n", + " @property\n", + " def r(self):\n", + " return self.frame.r\n", "\n", " def __repr__(self):\n", " return f\"GeneralizedStabilizerState(\\nframe=\\n{self.frame},\\nbasis=\\n{self.basis},\\nchi=\\n{self.chi})\"\n", @@ -106,7 +110,7 @@ " def evolve_by(self, pauli_channel):\n", " # Evolve the state by a Pauli channel.\n", " bs, cs, ps = pauli_decompose(pauli_channel.paulis.gs, pauli_channel.paulis.ps, \n", - " self.frame.gs, self.frame.ps)\n", + " self.frame.gs, self.frame.ps, self.r)\n", " # construct new basis, compute fusion map and fusion phase indicator\n", " L_old = self.basis.shape[0]\n", " L_add = bs.shape[0]\n", @@ -143,10 +147,10 @@ " mats = self.represent(PauliList(ops.gs,ops.ps))\n", " return numpy.tensordot(ops.cs, mats, axes=(0,0))\n", " if isinstance(ops, PauliList):\n", - " bs, cs, ps = pauli_decompose(ops.gs, ops.ps, self.frame.gs, self.frame.ps)\n", + " bs, cs, ps = pauli_decompose(ops.gs, ops.ps, self.frame.gs, self.frame.ps, self.r)\n", " L = self.basis.shape[0]\n", " K = ps.shape[0]\n", - " mats = numpy.zeros((K,L,L), dtype=numpy.complex_)\n", + " mats = numpy.zeros((K,L,L), dtype=numpy.complex128)\n", " idx = {tuple(b): i for i, b in enumerate(self.basis)}\n", " for k in range(K):\n", " for i0 in range(L):\n", @@ -173,11 +177,11 @@ " if mats.ndim == 2: # single matrix case\n", " return numpy.trace(self.chi @ mats)\n", " else: # batch case\n", - " return numpy.trace(numpy.matmul(self.chi, mats), axis1=1, axis2=2)\n", + " return numpy.tensordot(self.chi, mats, axes=([1,0], [1,2]))\n", " \n", " def to_numpy(self):\n", " '''Convert generalized stabilizer state to numpy density matrix representation.'''\n", - " destabilizers = self.frame[self.N:2*self.N]\n", + " destabilizers = self.frame[self.N+self.r:2*self.N]\n", " gs, ps =pauli_combine(self.basis, destabilizers.gs, destabilizers.ps)\n", " Ds = PauliList(gs,ps).to_numpy()\n", " rho0 = self.frame.to_numpy()\n", @@ -190,9 +194,17 @@ " " ] }, + { + "cell_type": "markdown", + "id": "161714ca", + "metadata": {}, + "source": [ + "Test: create a generalized stabilizer state over three qubits, where the first two qubits is in a mixture of Bell state, and the third qubit is maximally mixed." + ] + }, { "cell_type": "code", - "execution_count": 130, + "execution_count": 25, "id": "3e63172d", "metadata": {}, "outputs": [ @@ -202,8 +214,8 @@ "GeneralizedStabilizerState(\n", "frame=\n", "StabilizerState(\n", - " +XX\n", - " +ZZ),\n", + " +XXI\n", + " +ZZI),\n", "basis=\n", "[[0 0]\n", " [1 0]],\n", @@ -212,88 +224,164 @@ " [0.1+0.j 0.2+0.j]])" ] }, - "execution_count": 130, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state=GeneralizedStabilizerState(\n", - " pc.stabilizer_state(\"XX\",\"ZZ\"), \n", + " pc.stabilizer_state(\"XXI\",\"ZZI\"), \n", " numpy.array([[0,0],[1,0]]), \n", - " numpy.array([[0.8,0.1],[0.1,0.2]],dtype=numpy.complex_))\n", + " numpy.array([[0.8,0.1],[0.1,0.2]],dtype=numpy.complex128))\n", "state" ] }, + { + "cell_type": "markdown", + "id": "7ecc9b19", + "metadata": {}, + "source": [ + "Density matrix of the state." + ] + }, { "cell_type": "code", - "execution_count": 131, + "execution_count": 24, "id": "be822ce2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[0.6+0.j, 0. +0.j, 0. +0.j, 0.3+0.j],\n", - " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", - " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", - " [0.3+0.j, 0. +0.j, 0. +0.j, 0.4+0.j]])" + "array([[0.4 , 0. , 0. , 0. , 0.05, 0. , 0. , 0. ],\n", + " [0. , 0.4 , 0. , 0. , 0. , 0.05, 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0.05, 0. , 0. , 0. , 0.1 , 0. , 0. , 0. ],\n", + " [0. , 0.05, 0. , 0. , 0. , 0.1 , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])" ] }, - "execution_count": 131, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "state.to_numpy()" + "state.to_numpy().real" + ] + }, + { + "cell_type": "markdown", + "id": "39e9fb12", + "metadata": {}, + "source": [ + "Represent a list of Pauli operators on the stabilizer frame basis." ] }, { "cell_type": "code", - "execution_count": 134, + "execution_count": 18, + "id": "3744ddcb", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[[ 1.+0.j, 0.+0.j],\n", + " [ 0.+0.j, -1.+0.j]],\n", + "\n", + " [[ 0.+0.j, 0.+1.j],\n", + " [-0.-1.j, 0.+0.j]],\n", + "\n", + " [[ 0.+0.j, 0.+1.j],\n", + " [-0.-1.j, 0.+0.j]],\n", + "\n", + " [[-1.+0.j, 0.+0.j],\n", + " [ 0.+0.j, 1.+0.j]]])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state.represent(pc.paulis(\"XXI\",\"XYI\",\"YXI\",\"YYI\"))" + ] + }, + { + "cell_type": "markdown", + "id": "4f06773a", + "metadata": {}, + "source": [ + "Compute the expectation value of a Pauli polynomial on the state." + ] + }, + { + "cell_type": "code", + "execution_count": 26, "id": "9fb4ecae", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(0.3+0j)" + "np.complex128(0.3+0j)" ] }, - "execution_count": 134, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "state.expect((pc.pauli(\"XX\")+1j*pc.pauli(\"XY\")+1j*pc.pauli(\"YX\")-pc.pauli(\"YY\"))/4)" + "obs = (pc.pauli(\"XXI\")+1j*pc.pauli(\"XYI\")+1j*pc.pauli(\"YXI\")-pc.pauli(\"YYI\"))/4\n", + "state.expect(obs)" + ] + }, + { + "cell_type": "markdown", + "id": "493a4c4f", + "metadata": {}, + "source": [ + "Verify that the result is correct by explicit calculation using matrix representations of state and operator." ] }, { "cell_type": "code", - "execution_count": 135, + "execution_count": 27, "id": "7155800c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(0.30000000000000004+0j)" + "np.complex128(0.30000000000000004+0j)" ] }, - "execution_count": 135, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "numpy.trace(state.to_numpy() @ ((pc.pauli(\"XX\")+1j*pc.pauli(\"XY\")+1j*pc.pauli(\"YX\")-pc.pauli(\"YY\"))/4).to_numpy())" + "numpy.trace(state.to_numpy() @ obs.to_numpy())" + ] + }, + { + "cell_type": "markdown", + "id": "b7231ecf", + "metadata": {}, + "source": [ + "Define a Pauli channel which maximally decohere the first qubit in the Z basis, then apply the channel to the state." ] }, { "cell_type": "code", - "execution_count": 136, + "execution_count": 20, "id": "9dcaff6b", "metadata": {}, "outputs": [ @@ -303,8 +391,8 @@ "GeneralizedStabilizerState(\n", "frame=\n", "StabilizerState(\n", - " +XX\n", - " +ZZ),\n", + " +XXI\n", + " +ZZI),\n", "basis=\n", "[[0 0]\n", " [1 0]],\n", @@ -313,38 +401,50 @@ " [0.1+0.j 0.5+0.j]])" ] }, - "execution_count": 136, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "channel = PauliChannel(pc.paulis(\"II\",\"ZI\"), numpy.array([[0.5,0.0],[0.0,0.5]],dtype=numpy.complex_))\n", + "channel = PauliChannel(pc.paulis(\"III\",\"ZII\"), numpy.array([[0.5,0.0],[0.0,0.5]],dtype=numpy.complex128))\n", "state.evolve_by(channel)" ] }, + { + "cell_type": "markdown", + "id": "3b71c863", + "metadata": {}, + "source": [ + "The density matrix shows that the decoherence indeed removes off-diagonal matrix elements." + ] + }, { "cell_type": "code", - "execution_count": 137, + "execution_count": 21, "id": "3a6649ac", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[0.6+0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", - " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", - " [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],\n", - " [0. +0.j, 0. +0.j, 0. +0.j, 0.4+0.j]])" + "array([[0.3, 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0.3, 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0.2, 0. ],\n", + " [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2]])" ] }, - "execution_count": 137, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "state.to_numpy()" + "state.to_numpy().real" ] }, { @@ -372,7 +472,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.13.2" } }, "nbformat": 4, diff --git a/pyclifford/circuit.py b/pyclifford/circuit.py index 226b15a..252b12a 100644 --- a/pyclifford/circuit.py +++ b/pyclifford/circuit.py @@ -69,8 +69,8 @@ def independent_from(self, other): return len(set(self.qubits) & set(other.qubits))==0 def forward(self, obj): - '''Apply the gate forward in time. (inplace update) - Forward transformation: O -> U O U^\dagger + '''Apply the gate forward in time. (inplace update) + Forward transformation: O -> U O U^H Input: obj: Pauli, PauliList, StabilizerState - the object to be transformed @@ -103,7 +103,7 @@ def forward(self, obj): def backward(self, obj): '''Apply the gate backward in time. (inplace update) - Backward transformation: O -> U^\dagger O U + Backward transformation: O -> U^H O U Input: obj: Pauli, PauliList, StabilizerState - the object to be transformed @@ -195,12 +195,12 @@ def independent_from(self, other): def forward(self, obj): '''Implements the measurement (non-deterministic sampling outcome). (inplace update) - Forward transformation: rho -> M rho M^\dagger / Tr(M rho M^\dagger), - where M is sampled with probability Tr(M rho M^\dagger). + Forward transformation: rho -> M rho M^H / Tr(M rho M^H), + where M is sampled with probability Tr(M rho M^H). Input: obj: StabilizerState - the state to be measured - + Output: obj: StabilizerState - the state after measurement log2prob: real - log2 probability of the outcome''' @@ -215,12 +215,12 @@ def forward(self, obj): def backward(self, obj): '''Postselect the measurement outcome (deterministic projection). (inplace update) - Backward transformation: rho -> M^\dagger rho M / Tr(M^\dagger rho M) + Backward transformation: rho -> M^H rho M / Tr(M^H rho M) where M is fixed by the measurement outcome generated in the forward pass. Input: obj: StabilizerState - the state to be postselected - + Output: obj: StabilizerState - the state after postselection log2prob: real - log2 probability of successful postselection''' diff --git a/pyclifford/paulialg.py b/pyclifford/paulialg.py index 1e7e3b1..9d4187b 100644 --- a/pyclifford/paulialg.py +++ b/pyclifford/paulialg.py @@ -48,6 +48,11 @@ def __repr__(self): def N(self): # number of qubits return self.g.shape[0]//2 + def expand(self, N): + if N is not None and N > self.N: + self.g = numpy.concatenate([self.g, numpy.zeros(2*(N-self.N), dtype=self.g.dtype)]) + return self + def __neg__(self): return type(self)(self.g, (self.p + 2) % 4) @@ -77,6 +82,10 @@ def __sub__(self, other): def __matmul__(self, other): if isinstance(other, Pauli): + if self.N != other.N: + N = max(self.N, other.N) + self.expand(N) + other.expand(N) p = (self.p + other.p + ipow(self.g, other.g)) % 4 g = (self.g + other.g) % 2 return Pauli(g, p) @@ -186,6 +195,11 @@ def L(self): @property def N(self): return self.gs.shape[1]//2 + + def expand(self, N): + if N is not None and N > self.N: + self.gs = numpy.concatenate([self.gs, numpy.zeros((self.L, 2*(N-self.N)), dtype=self.gs.dtype)], axis=1) + return self def __getitem__(self, item): if isinstance(item, (int, numpy.integer)): @@ -367,7 +381,7 @@ def as_polynomial(self): '''cast the Pauli monomial to a single-term Pauli polynomial''' gs = numpy.expand_dims(self.g, 0) ps = numpy.array([self.p], dtype=numpy.int_) - cs = numpy.array([self.c], dtype=numpy.complex_) + cs = numpy.array([self.c], dtype=numpy.complex128) return PauliPolynomial(gs, ps).set_cs(cs) def inverse(self): @@ -386,7 +400,7 @@ class PauliPolynomial(PauliList): cs: comlex (L) - coefficients.''' def __init__(self, *args, **kwargs): super(PauliPolynomial, self).__init__(*args, **kwargs) - self.cs = numpy.ones(self.ps.shape, dtype=numpy.complex_) # default coefficient + self.cs = numpy.ones(self.ps.shape, dtype=numpy.complex128) # default coefficient def __repr__(self): txt = '' @@ -420,6 +434,10 @@ def __add__(self, other): other = other.as_polynomial() else: # otherwise assuming other is a number other = other * pauli_identity(self.N) + if self.N != other.N: + N = max(self.N, other.N) + self.expand(N) + other.expand(N) gs = numpy.concatenate([self.gs, other.gs]) ps = numpy.concatenate([self.ps, other.ps]) cs = numpy.concatenate([self.cs, other.cs]) @@ -436,6 +454,10 @@ def __matmul__(self, other): other = other.as_polynomial() else: raise NotImplementedError('matmul is not implemented for between {} and {}'.format(type(self).__name__, type(other).__name__)) + if self.N != other.N: + N = max(self.N, other.N) + self.expand(N) + other.expand(N) gs, ps, cs = batch_dot(self.gs, self.ps, self.cs, other.gs, other.ps, other.cs) return PauliPolynomial(gs, ps).set_cs(cs) @@ -474,13 +496,13 @@ def to_numpy(self): # ---- constructors ---- def pauli(obj, N = None): if isinstance(obj, Pauli): - return obj + return obj.expand(N) elif isinstance(obj, (tuple, list, numpy.ndarray)): N = len(obj) inds = enumerate(obj) elif isinstance(obj, dict): if N is None: - raise ValueError('pauli(inds, N) must specify qubit number N when inds is dict.') + N = max(obj.keys()) + 1 if obj else 0 inds = obj.items() elif isinstance(obj, str): return pauli(list(obj)) @@ -490,6 +512,7 @@ def pauli(obj, N = None): h = 0 p = 0 for i, mu in inds: + i = int(i) assert i-h < N, 'qubit {} is out of bounds for system size {}.'.format(i, N) if mu == 0 or mu == 'I': continue @@ -532,6 +555,9 @@ def paulis(*objs, N = None): objs = objs[0] # otherwise construct data for Pauli operators objs = [pauli(obj, N = N) for obj in objs] + N = max(obj.N for obj in objs) + for obj in objs: + obj.expand(N) gs = numpy.stack([obj.g for obj in objs]) ps = numpy.array([obj.p for obj in objs]) return PauliList(gs, ps) diff --git a/pyclifford/stabilizer.py b/pyclifford/stabilizer.py index 9cd1764..148679f 100644 --- a/pyclifford/stabilizer.py +++ b/pyclifford/stabilizer.py @@ -38,6 +38,12 @@ def __repr__(self): lns2 = [dis.format(xz[i%2], i//2, pauli) for i, pauli in zip(range(self.L-10, self.L), self[-10:])] return 'CliffordMap(\n{}\n ...\n{})'.format('\n'.join(lns1),'\n'.join(lns2)).replace('\n','\n ') + def expand(self, N): + if N is not None and N > self.N: + return identity_map(N).embed(self, mask(range(self.N), N)) + else: + return self + def copy(self): return CliffordMap(self.gs.copy(), self.ps.copy()) @@ -106,6 +112,12 @@ def __repr__(self): @property def stabilizers(self): return self[self.r:self.N] + + def expand(self, N): + if N is not None and N > self.N: + return stabilizer_state(self.stabilizers.expand(N)) + else: + return self def copy(self): return StabilizerState(self.gs.copy(), self.ps.copy(), r=self.r) diff --git a/pyclifford/utils.py b/pyclifford/utils.py index 839c3b3..2dcb89b 100644 --- a/pyclifford/utils.py +++ b/pyclifford/utils.py @@ -169,7 +169,7 @@ def batch_dot(gs1, ps1, cs1, gs2, ps2, cs2): (L2, N2) = gs2.shape gs = numpy.empty((L1,L2,N2), dtype=numpy.int_) ps = numpy.empty((L1,L2), dtype=numpy.int_) - cs = numpy.empty((L1,L2), dtype=numpy.complex_) + cs = numpy.empty((L1,L2), dtype=numpy.complex128) for j1 in range(L1): for j2 in range(L2): ps[j1,j2] = (ps1[j1] + ps2[j2] + ipow(gs1[j1], gs2[j2]))%4 @@ -247,35 +247,36 @@ def pauli_transform(gs_in, ps_in, gs_map, ps_map): return gs_out, ps_out @njit -def pauli_decompose(gs_in, ps_in, gs_stb, ps_stb): - '''Decompose a Pauli operator into stabilizer and destabilizers. +def pauli_decompose(gs_in, ps_in, gs_stb, ps_stb, r): + '''Decompose Pauli operators into stabilizer and destabilizers. Parameters: gs_in: int (L, 2*N) - Pauli strings in binary representation. ps_in: int (L) - phase indicators of Pauli operators. gs_stb: int (2*N, 2*N) - stabilizer tableaux in binary representation. ps_stb: int (2*N) - phase indicators of (de)stabilizer. + r: int - number of standby stabilizer pairs (0:r to be excluded). Returns: - bs_out: int (L, N) - binary encoding of destabilizer decomposition. - cs_out: int (L, N) - binary encoding of stabilizer decomposition. + bs_out: int (L, N-r) - binary encoding of destabilizer decomposition. + cs_out: int (L, N-r) - binary encoding of stabilizer decomposition. ps_out: int (L) - phase indicators of decomposed operators.''' (L, N2) = gs_in.shape N = N2//2 - bs_out = numpy.zeros((L, N), dtype=numpy.int_) - cs_out = numpy.zeros((L, N), dtype=numpy.int_) + bs_out = numpy.zeros((L, N-r), dtype=numpy.int_) + cs_out = numpy.zeros((L, N-r), dtype=numpy.int_) ps_out = ps_in.copy() g_tmp = numpy.zeros(N2, dtype=numpy.int_) for k in range(L): g_tmp.fill(0) - for j in range(N): + for j in range(r, N): if acq(gs_in[k], gs_stb[j]): - bs_out[k,j] = 1 + bs_out[k,j-r] = 1 ps_out[k] = ps_out[k] - ps_stb[j+N] - ipow(g_tmp, gs_stb[j+N]) g_tmp = (g_tmp + gs_stb[j+N])%2 - for j in range(N): + for j in range(r, N): if acq(gs_in[k], gs_stb[j+N]): - cs_out[k,j] = 1 + cs_out[k,j-r] = 1 ps_out[k] = ps_out[k] - ps_stb[j] - ipow(g_tmp, gs_stb[j]) g_tmp = (g_tmp + gs_stb[j])%2 return bs_out, cs_out, ps_out%4 @@ -1008,7 +1009,7 @@ def aggregate(data_in, inds, l): @njit def calculate_chi(chi_old, phi, fusion_map, fusion_p, L_new): L_old, L_add = fusion_map.shape - chi_new = numpy.zeros((L_new,L_new), dtype=numpy.complex_) + chi_new = numpy.zeros((L_new,L_new), dtype=numpy.complex128) for i1 in range(L_old): for j1 in range(L_add): k1 = fusion_map[i1,j1]