diff --git a/docs/source/changes.rst b/docs/source/changes.rst index 242300ab4..6e4e28713 100644 --- a/docs/source/changes.rst +++ b/docs/source/changes.rst @@ -4,6 +4,14 @@ Changelog v2.4 (Development) ------------------ +- Operator estimation data structures introduced in v2.2 have changed. ExperimentSetting's had + two members: in_operator and out_operator. out_operator is unchanged. in_operator has been + renamed to in_state and its data type is now ``TensorProductState`` instead of ``PauliTerm``. + It was always an abuse of notation to interpret pauli operators as defining initial states. + Analogous to the Pauli helper functions sI, sX, sY, and sZ, TensorProductStates are constructed + by multiplying together terms generated by the helper functions plusX, minusX, plusY, minusY, + plusZ, and minusZ. This functionality enables process tomography and process DFE (gh-770). + v2.3 (January 28, 2019) ----------------------- diff --git a/pyquil/gate_matrices.py b/pyquil/gate_matrices.py index 396a19785..89250e002 100644 --- a/pyquil/gate_matrices.py +++ b/pyquil/gate_matrices.py @@ -287,6 +287,39 @@ def bitphase_flip_operators(p): 'bitphase_flip': bitphase_flip_operators, } +SIC0 = np.array([1, 0]) +SIC1 = np.array([1, np.sqrt(2)]) / np.sqrt(3) +SIC2 = np.array([1, np.exp(-np.pi * 2j / 3) * np.sqrt(2)]) / np.sqrt(3) +SIC3 = np.array([1, np.exp(np.pi * 2j / 3) * np.sqrt(2)]) / np.sqrt(3) +""" +The symmetric informationally complete POVMs for a qubit. + +These can reduce the number of experiments to perform quantum process tomography. +For more information, please see http://info.phys.unm.edu/~caves/reports/infopovm.pdf +""" + +STATES = { + 'X': [ + np.array([1, 1]) / np.sqrt(2), + np.array([1, -1]) / np.sqrt(2), + ], + 'Y': [ + np.array([1, 1j]) / np.sqrt(2), + np.array([1, -1j]) / np.sqrt(2), + ], + 'Z': [ + np.array([1, 0]), + np.array([0, 1]), + ], + 'SIC': [ + SIC0, + SIC1, + SIC2, + SIC3, + ] +} + __all__ = list(QUANTUM_GATES.keys()) + ['relaxation_operators', 'dephasing_operators', 'depolarizing_operators', 'phase_flip_operators', - 'bit_flip_operators', 'bitphase_flip_operators'] + 'bit_flip_operators', 'bitphase_flip_operators', + 'STATES', 'SIC0', 'SIC1', 'SIC2', 'SIC3'] diff --git a/pyquil/operator_estimation.py b/pyquil/operator_estimation.py index 6f2bd1e9a..2aa67ee57 100644 --- a/pyquil/operator_estimation.py +++ b/pyquil/operator_estimation.py @@ -1,10 +1,15 @@ +import functools import itertools import json import logging +import re import sys +import warnings from json import JSONEncoder from math import pi -from typing import List, Union, Iterable, Dict +from operator import mul +from typing import List, Union, Iterable, Dict, Tuple + import networkx as nx import numpy as np from networkx.algorithms.approximation.clique import clique_removal @@ -23,6 +28,129 @@ @dataclass(frozen=True) +class _OneQState: + """ + A description of a named one-qubit quantum state. + + This can be used to generate pre-rotations for quantum process tomography. For example, + X0_14 will generate the +1 eigenstate of the X operator on qubit 14. X1_14 will generate the + -1 eigenstate. SIC0_14 will generate the 0th SIC-basis state on qubit 14. + """ + label: str + index: int + qubit: int + + def __str__(self): + return f'{self.label}{self.index}_{self.qubit}' + + @classmethod + def from_str(cls, s): + ma = re.match(r'\s*(\w+)(\d+)_(\d+)\s*', s) + if ma is None: + raise ValueError(f"Couldn't parse '{s}'") + return _OneQState( + label=ma.group(1), + index=int(ma.group(2)), + qubit=int(ma.group(3)), + ) + + +@dataclass(frozen=True) +class TensorProductState: + """ + A description of a multi-qubit quantum state that is a tensor product of many _OneQStates + states. + """ + states: Tuple[_OneQState] + + def __init__(self, states=None): + if states is None: + states = tuple() + object.__setattr__(self, 'states', tuple(states)) + + def __mul__(self, other): + return TensorProductState(self.states + other.states) + + def __str__(self): + return ' * '.join(str(s) for s in self.states) + + def __repr__(self): + return f'TensorProductState[{self}]' + + def __getitem__(self, qubit): + """Return the _OneQState at the given qubit.""" + for oneq_state in self.states: + if oneq_state.qubit == qubit: + return oneq_state + raise IndexError() + + def __len__(self): + return len(self.states) + + def states_as_set(self): + return frozenset(self.states) + + def __eq__(self, other): + if not isinstance(other, TensorProductState): + return False + + return self.states_as_set() == other.states_as_set() + + def __hash__(self): + return hash(self.states_as_set()) + + @classmethod + def from_str(cls, s): + if s == '': + return TensorProductState() + return TensorProductState(tuple(_OneQState.from_str(x) for x in s.split('*'))) + + +def SIC0(q): + return TensorProductState((_OneQState('SIC', 0, q),)) + + +def SIC1(q): + return TensorProductState((_OneQState('SIC', 1, q),)) + + +def SIC2(q): + return TensorProductState((_OneQState('SIC', 2, q),)) + + +def SIC3(q): + return TensorProductState((_OneQState('SIC', 3, q),)) + + +def plusX(q): + return TensorProductState((_OneQState('X', 0, q),)) + + +def minusX(q): + return TensorProductState((_OneQState('X', 1, q),)) + + +def plusY(q): + return TensorProductState((_OneQState('Y', 0, q),)) + + +def minusY(q): + return TensorProductState((_OneQState('Y', 1, q),)) + + +def plusZ(q): + return TensorProductState((_OneQState('Z', 0, q),)) + + +def minusZ(q): + return TensorProductState((_OneQState('Z', 1, q),)) + + +def zeros_state(qubits: Iterable[int]): + return TensorProductState(_OneQState('Z', 0, q) for q in qubits) + + +@dataclass(frozen=True, init=False) class ExperimentSetting: """ Input and output settings for a tomography-like experiment. @@ -38,11 +166,45 @@ class ExperimentSetting: number of these :py:class:`ExperimentSetting` objects will be created and grouped into a :py:class:`TomographyExperiment`. """ - in_operator: PauliTerm + in_state: TensorProductState out_operator: PauliTerm + def __init__(self, in_state: TensorProductState, out_operator: PauliTerm): + # For backwards compatibility, handle in_state specified by PauliTerm. + if isinstance(in_state, PauliTerm): + warnings.warn("Please specify in_state as a TensorProductState", + DeprecationWarning, stacklevel=2) + + if is_identity(in_state): + in_state = TensorProductState() + else: + in_state = TensorProductState([ + _OneQState(label=pauli_label, index=0, qubit=qubit) + for qubit, pauli_label in in_state._ops.items() + ]) + + object.__setattr__(self, 'in_state', in_state) + object.__setattr__(self, 'out_operator', out_operator) + + @property + def in_operator(self): + warnings.warn("ExperimentSetting.in_operator is deprecated in favor of in_state", + stacklevel=2) + + # Backwards compat + pt = sI() + for oneq_state in self.in_state.states: + if oneq_state.label not in ['X', 'Y', 'Z']: + raise ValueError(f"Can't shim {oneq_state.label} into a pauli term. Use in_state.") + if oneq_state.index != 0: + raise ValueError(f"Can't shim {oneq_state} into a pauli term. Use in_state.") + + pt *= PauliTerm(op=oneq_state.label, index=oneq_state.qubit) + + return pt + def __str__(self): - return f'{self.in_operator.compact_str()}→{self.out_operator.compact_str()}' + return f'{self.in_state}→{self.out_operator.compact_str()}' def __repr__(self): return f'ExperimentSetting[{self}]' @@ -54,7 +216,7 @@ def serializable(self): def from_str(cls, s: str): """The opposite of str(expt)""" instr, outstr = s.split('→') - return ExperimentSetting(in_operator=PauliTerm.from_compact_str(instr), + return ExperimentSetting(in_state=TensorProductState.from_str(instr), out_operator=PauliTerm.from_compact_str(outstr)) @@ -239,115 +401,83 @@ def read_json(fn): return json.load(f, object_hook=_operator_object_hook) -def _local_pauli_eig_prep(op: str, idx: int): - """ - Generate gate sequence to prepare a the +1 eigenstate of a Pauli operator, assuming - we are starting from the |00..00> ground state. - - :param op: A string representation of the Pauli operator whose +1 eigenstate we'd like to - prepare. - :param idx: The index of the qubit that the preparation is acting on - :return: A program which will prepare the requested state. - """ - if op == 'X': - return Program(RY(pi / 2, idx)) - elif op == 'Y': - return Program(RX(-pi / 2, idx)) - elif op == 'Z': +def _one_q_sic_prep(index, qubit): + """Prepare the index-th SIC basis state.""" + if index == 0: return Program() - raise ValueError(f'Unknown operation {op}') + theta = 2 * np.arccos(1 / np.sqrt(3)) + zx_plane_rotation = Program([ + RX(-pi / 2, qubit), + RZ(theta - pi, qubit), + RX(-pi / 2, qubit), + ]) + if index == 1: + return zx_plane_rotation -def _local_pauli_eig_meas(op, idx): - """ - Generate gate sequence to measure in the eigenbasis of a Pauli operator, assuming - we are only able to measure in the Z eigenbasis. - - """ - if op == 'X': - return Program(RY(-pi / 2, idx)) - elif op == 'Y': - return Program(RX(pi / 2, idx)) - elif op == 'Z': - return Program() - raise ValueError(f'Unknown operation {op}') + elif index == 2: + return zx_plane_rotation + RZ(-2 * pi / 3, qubit) + elif index == 3: + return zx_plane_rotation + RZ(2 * pi / 3, qubit) -def _ops_commute(op_code1: str, op_code2: str): - """ - Given two 1q op strings (I, X, Y, or Z), determine whether they commute + raise ValueError(f'Bad SIC index: {index}') - :param op_code1: First operation - :param op_code2: Second operation - :return: Boolean specifying whether the two operations commute - """ - if op_code1 not in ['X', 'Y', 'Z', 'I']: - raise ValueError(f"Unknown op_code {op_code1}") - if op_code2 not in ['X', 'Y', 'Z', 'I']: - raise ValueError(f"Unknown op_code {op_code2}") - - if op_code1 == op_code2: - # Same op - return True - elif op_code1 == 'I' or op_code2 == 'I': - # I commutes with everything - return True - else: - # Otherwise, they do not commute. - return False +def _one_q_pauli_prep(label, index, qubit): + """Prepare the index-th eigenstate of the pauli operator given by label.""" + if index not in [0, 1]: + raise ValueError(f'Bad Pauli index: {index}') -def _all_qubits_diagonal_in_tpb(op1: PauliTerm, op2: PauliTerm): - """ - Compare all qubits between two PauliTerms to see if they are all diagonal in an - overall shared tensor product basis. More concretely, test if ``op1`` and - ``op2`` are diagonal in each others' "natural" tensor product basis. + if label == 'X': + if index == 0: + return Program(RY(pi / 2, qubit)) + else: + return Program(RY(-pi / 2, qubit)) - Given some PauliTerm, the 'natural' tensor product basis (tpb) to - diagonalize this term is the one which diagonalizes each Pauli operator in the - product term-by-term. + elif label == 'Y': + if index == 0: + return Program(RX(-pi / 2, qubit)) + else: + return Program(RX(pi / 2, qubit)) - For example, X(1) * Z(0) would be diagonal in the 'natural' tensor product basis - {(|0> +/- |1>)/Sqrt[2]} * {|0>, |1>}, whereas Z(1) * X(0) would be diagonal - in the 'natural' tpb {|0>, |1>} * {(|0> +/- |1>)/Sqrt[2]}. The two operators - commute but are not diagonal in each others 'natural' tpb (in fact, they are - anti-diagonal in each others 'natural' tpb). This function tests whether two - operators given as PauliTerms are both diagonal in each others 'natural' tpb. + elif label == 'Z': + if index == 0: + return Program() + else: + return Program(RX(pi, qubit)) - Note that for the given example of X(1) * Z(0) and Z(1) * X(0), we can construct - the following basis which simultaneously diagonalizes both operators: + raise ValueError(f'Bad Pauli label: {label}') - -- |0>' = |0> (|+>) + |1> (|->) - -- |1>' = |0> (|+>) - |1> (|->) - -- |2>' = |0> (|->) + |1> (|+>) - -- |3>' = |0> (-|->) + |1> (|+>) - In this basis, X Z looks like diag(1, -1, 1, -1), and Z X looks like diag(1, 1, -1, -1). - Notice however that this basis cannot be constructed with single-qubit operations, as each - of the basis vectors are entangled states. +def _one_q_state_prep(oneq_state: _OneQState): + """Prepare a one qubit state. - :param op1: PauliTerm to check diagonality of in the natural tpb of ``op2`` - :param op2: PauliTerm to check diagonality of in the natural tpb of ``op1`` - :return: Boolean of diagonality in each others natural tpb + Either SIC[0-3], X[0-1], Y[0-1], or Z[0-1]. """ - all_qubits = set(op1.get_qubits()) & set(op2.get_qubits()) - return all(_ops_commute(op1[q], op2[q]) for q in all_qubits) + label = oneq_state.label + if label == 'SIC': + return _one_q_sic_prep(oneq_state.index, oneq_state.qubit) + elif label in ['X', 'Y', 'Z']: + return _one_q_pauli_prep(label, oneq_state.index, oneq_state.qubit) + else: + raise ValueError(f"Bad state label: {label}") -def _expt_settings_diagonal_in_tpb(es1: ExperimentSetting, es2: ExperimentSetting): +def _local_pauli_eig_meas(op, idx): """ - Extends the concept of being diagonal in the same tpb (see :py:func:_all_qubits_diagonal_in_tpb) - to ExperimentSettings, by determining if the pairs of in_operators and out_operators are - separately diagonal in the same tpb + Generate gate sequence to measure in the eigenbasis of a Pauli operator, assuming + we are only able to measure in the Z eigenbasis. - :param es1: ExperimentSetting to check diagonality of in the natural tpb of ``es2`` - :param es2: ExperimentSetting to check diagonality of in the natural tpb of ``es1`` - :return: Boolean of diagonality in each others natural tpb """ - in_dtpb = _all_qubits_diagonal_in_tpb(es1.in_operator, es2.in_operator) - out_dtpb = _all_qubits_diagonal_in_tpb(es1.out_operator, es2.out_operator) - return in_dtpb and out_dtpb + if op == 'X': + return Program(RY(-pi / 2, idx)) + elif op == 'Y': + return Program(RX(pi / 2, idx)) + elif op == 'Z': + return Program() + raise ValueError(f'Unknown operation {op}') def construct_tpb_graph(experiments: TomographyExperiment): @@ -371,7 +501,9 @@ def construct_tpb_graph(experiments: TomographyExperiment): if expt1 == expt2: continue - if _expt_settings_diagonal_in_tpb(expt1, expt2): + max_weight_in = _max_weight_state([expt1.in_state, expt2.in_state]) + max_weight_out = _max_weight_operator([expt1.out_operator, expt2.out_operator]) + if max_weight_in is not None and max_weight_out is not None: g.add_edge(expt1, expt2) return g @@ -400,34 +532,46 @@ def group_experiments_clique_removal(experiments: TomographyExperiment) -> Tomog return TomographyExperiment(new_cliqs, program=experiments.program, qubits=experiments.qubits) -def _validate_all_diagonal_in_tpb(ops: Iterable[PauliTerm]) -> Dict[int, str]: - """Each non-identity qubit should result in the same op_str among all operations. Return - said mapping. +def _max_weight_operator(ops: Iterable[PauliTerm]) -> Union[None, PauliTerm]: + """Construct a PauliTerm operator by taking the non-identity single-qubit operator at each + qubit position. + + This function will return ``None`` if the input operators do not share a natural tensor + product basis. + + For example, the max_weight_operator of ["XI", "IZ"] is "XZ". Asking for the max weight + operator of something like ["XI", "ZI"] will return None. """ mapping = dict() # type: Dict[int, str] for op in ops: for idx, op_str in op: if idx in mapping: - assert mapping[idx] == op_str, 'Improper grouping of operators' + if mapping[idx] != op_str: + return None else: mapping[idx] = op_str - return mapping + op = functools.reduce(mul, (PauliTerm(op, q) for q, op in mapping.items()), sI()) + return op -def _get_diagonalizing_basis(ops: Iterable[PauliTerm]) -> PauliTerm: - """ - Find the Pauli Term with the most non-identity terms +def _max_weight_state(states: Iterable[TensorProductState]) -> Union[None, TensorProductState]: + """Construct a TensorProductState by taking the single-qubit state at each + qubit position. + + This function will return ``None`` if the input states are not compatible - :param ops: Iterable of PauliTerms to check - :return: The highest weight PauliTerm from the input iterable + For example, the max_weight_state of ["(+X, q0)", "(-Z, q1)"] is "(+X, q0; -Z q1)". Asking for + the max weight state of something like ["(+X, q0)", "(+Z, q0)"] will return None. """ - # obtain qubit: operation mapping - dict_pt = _validate_all_diagonal_in_tpb(ops) - # convert this mapping to PauliTerm - pt = sI() - for q, op in dict_pt.items(): - pt *= PauliTerm(op, q) - return pt + mapping = dict() # type: Dict[int, _OneQState] + for state in states: + for oneq_state in state.states: + if oneq_state.qubit in mapping: + if mapping[oneq_state.qubit] != oneq_state: + return None + else: + mapping[oneq_state.qubit] = oneq_state + return TensorProductState(list(mapping.values())) def _max_tpb_overlap(tomo_expt: TomographyExperiment): @@ -452,26 +596,25 @@ def _max_tpb_overlap(tomo_expt: TomographyExperiment): found_tpb = False # loop through dict items for es, es_list in diagonal_sets.items(): - # update the dict value if es is diagonal in the same tpb as expt_setting - if _expt_settings_diagonal_in_tpb(es, expt_setting): - # shared tpb was found + trial_es_list = es_list + [expt_setting] + diag_in_term = _max_weight_state(expst.in_state for expst in trial_es_list) + diag_out_term = _max_weight_operator(expst.out_operator for expst in trial_es_list) + # max_weight_xxx returns None if the set of xxx's don't share a TPB, so the following + # conditional is True if expt_setting can be inserted into the current es_list. + if diag_in_term is not None and diag_out_term is not None: found_tpb = True - # determine the updated list of ExperimentSettings - updated_es_list = es_list + [expt_setting] - # obtain the diagonalizing bases for both the updated in and out sets - diag_in_term = _get_diagonalizing_basis([expst.in_operator for expst in updated_es_list]) - diag_out_term = _get_diagonalizing_basis([expst.out_operator for expst in updated_es_list]) - assert len(diag_in_term) >= len(es.in_operator), \ - "Highest weight in-PauliTerm can't be smaller than the given in-PauliTerm" + assert len(diag_in_term) >= len(es.in_state), \ + "Highest weight in-state can't be smaller than the given in-state" assert len(diag_out_term) >= len(es.out_operator), \ "Highest weight out-PauliTerm can't be smaller than the given out-PauliTerm" + # update the diagonalizing basis (key of dict) if necessary - if len(diag_in_term) > len(es.in_operator) or len(diag_out_term) > len(es.out_operator): + if len(diag_in_term) > len(es.in_state) or len(diag_out_term) > len(es.out_operator): del diagonal_sets[es] new_es = ExperimentSetting(diag_in_term, diag_out_term) - diagonal_sets[new_es] = updated_es_list + diagonal_sets[new_es] = trial_es_list else: - diagonal_sets[es] = updated_es_list + diagonal_sets[es] = trial_es_list break if not found_tpb: @@ -497,10 +640,51 @@ def group_experiments_greedy(tomo_expt: TomographyExperiment): return grouped_tomo_expt -def group_experiments(experiments: TomographyExperiment, method: str = 'greedy') -> TomographyExperiment: +def group_experiments(experiments: TomographyExperiment, + method: str = 'greedy') -> TomographyExperiment: """ Group experiments that are diagonal in a shared tensor product basis (TPB) to minimize number - of QPU runs, using a specified method (greedy method by default) + of QPU runs. + + Background + ---------- + + Given some PauliTerm operator, the 'natural' tensor product basis to + diagonalize this term is the one which diagonalizes each Pauli operator in the + product term-by-term. + + For example, X(1) * Z(0) would be diagonal in the 'natural' tensor product basis + {(|0> +/- |1>)/Sqrt[2]} * {|0>, |1>}, whereas Z(1) * X(0) would be diagonal + in the 'natural' tpb {|0>, |1>} * {(|0> +/- |1>)/Sqrt[2]}. The two operators + commute but are not diagonal in each others 'natural' tpb (in fact, they are + anti-diagonal in each others 'natural' tpb). This function tests whether two + operators given as PauliTerms are both diagonal in each others 'natural' tpb. + + Note that for the given example of X(1) * Z(0) and Z(1) * X(0), we can construct + the following basis which simultaneously diagonalizes both operators: + + -- |0>' = |0> (|+>) + |1> (|->) + -- |1>' = |0> (|+>) - |1> (|->) + -- |2>' = |0> (|->) + |1> (|+>) + -- |3>' = |0> (-|->) + |1> (|+>) + + In this basis, X Z looks like diag(1, -1, 1, -1), and Z X looks like diag(1, 1, -1, -1). + Notice however that this basis cannot be constructed with single-qubit operations, as each + of the basis vectors are entangled states. + + + Methods + ------- + + The "greedy" method will keep a running set of 'buckets' into which grouped ExperimentSettings + will be placed. Each new ExperimentSetting considered is assigned to the first applicable + bucket and a new bucket is created if there are no applicable buckets. + + The "clique-removal" method maps the term grouping problem onto Max Clique graph problem. + This method constructs a NetworkX graph where an edge exists between two settings that + share an nTPB and then uses networkx's algorithm for clique removal. This method can give + you marginally better groupings in certain circumstances, but constructing the + graph is pretty slow so "greedy" is the default. :param experiments: a tomography experiment :param method: method used for grouping; the allowed methods are one of @@ -565,21 +749,21 @@ def measure_observables(qc: QuantumComputer, tomo_experiment: TomographyExperime # estimate. log.info(f"Collecting bitstrings for the {len(settings)} settings: {settings}") - # 1.1 Prepare a state according to setting.in_operator + # 1.1 Prepare a state according to the amalgam of all setting.in_state total_prog = Program() if active_reset: total_prog += RESET() - in_mapping = _validate_all_diagonal_in_tpb(setting.in_operator for setting in settings) - for idx, op_str in in_mapping.items(): - total_prog += _local_pauli_eig_prep(op_str, idx) + max_weight_in_state = _max_weight_state(setting.in_state for setting in settings) + for oneq_state in max_weight_in_state.states: + total_prog += _one_q_state_prep(oneq_state) # 1.2 Add in the program total_prog += tomo_experiment.program # 1.3 Measure the state according to setting.out_operator - out_mapping = _validate_all_diagonal_in_tpb(setting.out_operator for setting in settings) - for idx, op_str in out_mapping.items(): - total_prog += _local_pauli_eig_meas(op_str, idx) + max_weight_out_op = _max_weight_operator(setting.out_operator for setting in settings) + for qubit, op_str in max_weight_out_op: + total_prog += _local_pauli_eig_meas(op_str, qubit) # 2. Run the experiment bitstrings = qc.run_and_measure(total_prog, n_shots) @@ -596,8 +780,6 @@ def measure_observables(qc: QuantumComputer, tomo_experiment: TomographyExperime # either the first column, second column, or both and multiplying along the row. for setting in settings: # 3.2 Get the term's coefficient so we can multiply it in later. - assert setting.in_operator.coefficient == 1, 'in_operator should specify a state and ' \ - 'therefore cannot have a coefficient' coeff = complex(setting.out_operator.coefficient) if not np.isclose(coeff.imag, 0): raise ValueError(f"{setting}'s out_operator has a complex coefficient.") diff --git a/pyquil/tests/test_operator_estimation.py b/pyquil/tests/test_operator_estimation.py index f3eb232b7..0548201fc 100644 --- a/pyquil/tests/test_operator_estimation.py +++ b/pyquil/tests/test_operator_estimation.py @@ -1,21 +1,32 @@ +import functools import itertools import random from math import pi - -import numpy as np -import functools from operator import mul +import numpy as np import pytest +from pyquil import Program, get_qc from pyquil.api import WavefunctionSimulator +from pyquil.gates import * from pyquil.operator_estimation import ExperimentSetting, TomographyExperiment, to_json, read_json, \ - _all_qubits_diagonal_in_tpb, group_experiments, ExperimentResult, measure_observables, \ - _get_diagonalizing_basis, _max_tpb_overlap, group_experiments_greedy, \ - _expt_settings_diagonal_in_tpb + group_experiments, ExperimentResult, measure_observables, SIC0, \ + SIC1, SIC2, SIC3, plusX, minusX, plusY, minusY, plusZ, minusZ, \ + _max_tpb_overlap, _max_weight_operator, _max_weight_state, \ + TensorProductState, zeros_state from pyquil.paulis import sI, sX, sY, sZ, PauliSum, PauliTerm -from pyquil import Program, get_qc -from pyquil.gates import * + + +def _generate_random_states(n_qubits, n_terms): + oneq_states = [SIC0, SIC1, SIC2, SIC3, plusX, minusX, plusY, minusY, plusZ, minusZ] + all_s_inds = np.random.randint(len(oneq_states), size=(n_terms, n_qubits)) + states = [] + for s_inds in all_s_inds: + state = functools.reduce(mul, (oneq_states[pi](i) for i, pi in enumerate(s_inds)), + TensorProductState([])) + states += [state] + return states def _generate_random_paulis(n_qubits, n_terms): @@ -29,19 +40,19 @@ def _generate_random_paulis(n_qubits, n_terms): return operators -def test_experiment(): - in_ops = _generate_random_paulis(n_qubits=4, n_terms=7) +def test_experiment_setting(): + in_states = _generate_random_states(n_qubits=4, n_terms=7) out_ops = _generate_random_paulis(n_qubits=4, n_terms=7) - for iop, oop in zip(in_ops, out_ops): - expt = ExperimentSetting(iop, oop) + for ist, oop in zip(in_states, out_ops): + expt = ExperimentSetting(ist, oop) assert str(expt) == expt.serializable() expt2 = ExperimentSetting.from_str(str(expt)) assert expt == expt2 - assert expt2.in_operator == iop + assert expt2.in_state == ist assert expt2.out_operator == oop -def test_experiment_no_in(): +def test_setting_no_in_back_compat(): out_ops = _generate_random_paulis(n_qubits=4, n_terms=7) for oop in out_ops: expt = ExperimentSetting(sI(), oop) @@ -51,7 +62,17 @@ def test_experiment_no_in(): assert expt2.out_operator == oop -def test_experiment_suite(): +def test_setting_no_in(): + out_ops = _generate_random_paulis(n_qubits=4, n_terms=7) + for oop in out_ops: + expt = ExperimentSetting(zeros_state(oop.get_qubits()), oop) + expt2 = ExperimentSetting.from_str(str(expt)) + assert expt == expt2 + assert expt2.in_operator == functools.reduce(mul, [sZ(q) for q in oop.get_qubits()], sI()) + assert expt2.out_operator == oop + + +def test_tomo_experiment(): expts = [ ExperimentSetting(sI(), sX(0) * sY(1)), ExperimentSetting(sZ(0), sZ(0)), @@ -72,7 +93,7 @@ def test_experiment_suite(): assert prog_str == 'X 0; Y 1' -def test_experiment_suite_pre_grouped(): +def test_tomo_experiment_pre_grouped(): expts = [ [ExperimentSetting(sI(), sX(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sX(1))], [ExperimentSetting(sI(), sZ(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sZ(1))], @@ -91,13 +112,13 @@ def test_experiment_suite_pre_grouped(): assert prog_str == 'X 0; Y 1' -def test_experiment_suite_empty(): +def test_tomo_experiment_empty(): suite = TomographyExperiment([], program=Program(X(0)), qubits=[0]) assert len(suite) == 0 assert str(suite.program) == 'X 0\n' -def test_suite_deser(tmpdir): +def test_experiment_deser(tmpdir): expts = [ [ExperimentSetting(sI(), sX(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sX(1))], [ExperimentSetting(sI(), sZ(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sZ(1))], @@ -113,43 +134,40 @@ def test_suite_deser(tmpdir): assert suite == suite2 -def test_all_ops_belong_to_tpb(): - expts = [ - [ExperimentSetting(sI(), sX(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sX(1))], - [ExperimentSetting(sI(), sZ(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sZ(1))], - ] - for group in expts: - for e1, e2 in itertools.combinations(group, 2): - assert _all_qubits_diagonal_in_tpb(e1.in_operator, e2.in_operator) - assert _all_qubits_diagonal_in_tpb(e1.out_operator, e2.out_operator) - - assert _all_qubits_diagonal_in_tpb(sZ(0), sZ(0) * sZ(1)) - assert _all_qubits_diagonal_in_tpb(sX(5), sZ(4)) - assert not _all_qubits_diagonal_in_tpb(sX(0), sY(0) * sZ(2)) - # this last example illustrates that a pair of commuting operators - # need not be diagonal in the same tpb - assert not _all_qubits_diagonal_in_tpb(sX(1) * sZ(0), sZ(1) * sX(0)) +@pytest.fixture(params=['clique-removal', 'greedy']) +def grouping_method(request): + return request.param -def test_group_experiments(): +def test_group_experiments(grouping_method): expts = [ # cf above, I removed the inner nesting. Still grouped visually ExperimentSetting(sI(), sX(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sX(1)), ExperimentSetting(sI(), sZ(0) * sI(1)), ExperimentSetting(sI(), sI(0) * sZ(1)), ] suite = TomographyExperiment(expts, Program(), qubits=[0, 1]) - grouped_suite = group_experiments(suite) + grouped_suite = group_experiments(suite, method=grouping_method) assert len(suite) == 4 assert len(grouped_suite) == 2 -def test_experiment_result(): +def test_experiment_result_compat(): er = ExperimentResult( setting=ExperimentSetting(sX(0), sZ(0)), expectation=0.9, stddev=0.05, total_counts=100, ) - assert str(er) == '(1+0j)*X0→(1+0j)*Z0: 0.9 +- 0.05' + assert str(er) == 'X0_0→(1+0j)*Z0: 0.9 +- 0.05' + + +def test_experiment_result(): + er = ExperimentResult( + setting=ExperimentSetting(plusX(0), sZ(0)), + expectation=0.9, + stddev=0.05, + total_counts=100, + ) + assert str(er) == 'X0_0→(1+0j)*Z0: 0.9 +- 0.05' def test_measure_observables(forest): @@ -238,14 +256,59 @@ def test_no_complex_coeffs(forest): res = list(measure_observables(qc, suite)) -def test_get_diagonalizing_basis_1(): - pauli_terms = [sZ(0), sX(1) * sZ(0), sY(2) * sX(1)] - assert _get_diagonalizing_basis(pauli_terms) == sY(2) * sX(1) * sZ(0) +def test_max_weight_operator_1(): + pauli_terms = [sZ(0), + sX(1) * sZ(0), + sY(2) * sX(1)] + assert _max_weight_operator(pauli_terms) == sY(2) * sX(1) * sZ(0) + + +def test_max_weight_operator_2(): + pauli_terms = [sZ(0), + sX(1) * sZ(0), + sY(2) * sX(1), + sZ(5) * sI(3)] + assert _max_weight_operator(pauli_terms) == sZ(5) * sY(2) * sX(1) * sZ(0) + +def test_max_weight_operator_3(): + pauli_terms = [sZ(0) * sX(5), + sX(1) * sZ(0), + sY(2) * sX(1), + sZ(5) * sI(3)] + assert _max_weight_operator(pauli_terms) is None -def test_get_diagonalizing_basis_2(): - pauli_terms = [sZ(0), sX(1) * sZ(0), sY(2) * sX(1), sZ(5) * sI(3)] - assert _get_diagonalizing_basis(pauli_terms) == sZ(5) * sY(2) * sX(1) * sZ(0) + +def test_max_weight_state_1(): + states = [plusX(0) * plusZ(1), + plusX(0), + plusZ(1), + ] + assert _max_weight_state(states) == states[0] + + +def test_max_weight_state_2(): + states = [plusX(1) * plusZ(0), + plusX(0), + plusZ(1), + ] + assert _max_weight_state(states) is None + + +def test_max_weight_state_3(): + states = [plusX(0) * minusZ(1), + plusX(0), + minusZ(1), + ] + assert _max_weight_state(states) == states[0] + + +def test_max_weight_state_4(): + states = [plusX(1) * minusZ(0), + plusX(0), + minusZ(1), + ] + assert _max_weight_state(states) is None def test_max_tpb_overlap_1(): @@ -254,9 +317,12 @@ def test_max_tpb_overlap_1(): tomo_expt_program = Program(H(0), H(1), H(2)) tomo_expt_qubits = [0, 1, 2] tomo_expt = TomographyExperiment(tomo_expt_settings, tomo_expt_program, tomo_expt_qubits) - expected_dict = {ExperimentSetting(sX(0) * sZ(1) * sX(2), sZ(0) * sY(1) * sY(2)): - [ExperimentSetting(sZ(1) * sX(0), sY(2) * sY(1)), - ExperimentSetting(sX(2) * sZ(1), sY(2) * sZ(0))]} + expected_dict = { + ExperimentSetting(plusX(0) * plusZ(1) * plusX(2), sZ(0) * sY(1) * sY(2)): [ + ExperimentSetting(plusZ(1) * plusX(0), sY(2) * sY(1)), + ExperimentSetting(plusX(2) * plusZ(1), sY(2) * sZ(0)) + ] + } assert expected_dict == _max_tpb_overlap(tomo_expt) @@ -288,29 +354,58 @@ def test_group_experiments_greedy(): PauliTerm.from_compact_str('(1+0j)*Z4X8Y5X3Y7Y1'))], [ExperimentSetting(sZ(7), sY(1))]], program=Program(H(0), H(1), H(2)), qubits=[0, 1, 2]) - grouped_tomo_expt = group_experiments_greedy(ungrouped_tomo_expt) + grouped_tomo_expt = group_experiments(ungrouped_tomo_expt, method='greedy') expected_grouped_tomo_expt = TomographyExperiment( - [[ExperimentSetting(PauliTerm.from_compact_str('(1+0j)*Z7Y8Z1Y4Z2Y5Y0X6'), - PauliTerm.from_compact_str('(1+0j)*Z4X8Y5X3Y7Y1')), - ExperimentSetting(sZ(7), sY(1))]], + [[ + ExperimentSetting(TensorProductState.from_str('Z0_7 * Y0_8 * Z0_1 * Y0_4 * ' + 'Z0_2 * Y0_5 * Y0_0 * X0_6'), + PauliTerm.from_compact_str('(1+0j)*Z4X8Y5X3Y7Y1')), + ExperimentSetting(plusZ(7), sY(1)) + ]], program=Program(H(0), H(1), H(2)), qubits=[0, 1, 2]) assert grouped_tomo_expt == expected_grouped_tomo_expt def test_expt_settings_diagonal_in_tpb(): - expt_setting1 = ExperimentSetting(sZ(1) * sX(0), sY(1) * sZ(0)) - expt_setting2 = ExperimentSetting(sY(2) * sZ(1), sZ(2) * sY(1)) + def _expt_settings_diagonal_in_tpb(es1: ExperimentSetting, es2: ExperimentSetting): + """ + Extends the concept of being diagonal in the same tpb to ExperimentSettings, by + determining if the pairs of in_states and out_operators are separately diagonal in the same + tpb + """ + max_weight_in = _max_weight_state([es1.in_state, es2.in_state]) + max_weight_out = _max_weight_operator([es1.out_operator, es2.out_operator]) + return max_weight_in is not None and max_weight_out is not None + + expt_setting1 = ExperimentSetting(plusZ(1) * plusX(0), sY(1) * sZ(0)) + expt_setting2 = ExperimentSetting(plusY(2) * plusZ(1), sZ(2) * sY(1)) assert _expt_settings_diagonal_in_tpb(expt_setting1, expt_setting2) - expt_setting3 = ExperimentSetting(sX(2) * sZ(1), sZ(2) * sY(1)) - expt_setting4 = ExperimentSetting(sY(2) * sZ(1), sX(2) * sY(1)) + expt_setting3 = ExperimentSetting(plusX(2) * plusZ(1), sZ(2) * sY(1)) + expt_setting4 = ExperimentSetting(plusY(2) * plusZ(1), sX(2) * sY(1)) assert not _expt_settings_diagonal_in_tpb(expt_setting2, expt_setting3) assert not _expt_settings_diagonal_in_tpb(expt_setting2, expt_setting4) def test_identity(forest): qc = get_qc('2q-qvm') - suite = TomographyExperiment([ExperimentSetting(sI(), 0.123 * sI(0))], + suite = TomographyExperiment([ExperimentSetting(plusZ(0), 0.123 * sI(0))], program=Program(X(0)), qubits=[0]) result = list(measure_observables(qc, suite))[0] assert result.expectation == 0.123 + + +def test_sic_process_tomo(forest): + qc = get_qc('2q-qvm') + process = Program(X(0)) + settings = [] + for in_state in [SIC0, SIC1, SIC2, SIC3]: + for out_op in [sI, sX, sY, sZ]: + settings += [ExperimentSetting( + in_state=in_state(q=0), + out_operator=out_op(q=0) + )] + + experiment = TomographyExperiment(settings=settings, program=process, qubits=[0]) + results = list(measure_observables(qc, experiment)) + assert len(results) == 4 * 4 diff --git a/pyquil/tests/test_unitary_tools.py b/pyquil/tests/test_unitary_tools.py index f947a6c5e..9b0b5a640 100644 --- a/pyquil/tests/test_unitary_tools.py +++ b/pyquil/tests/test_unitary_tools.py @@ -4,9 +4,10 @@ from pyquil import Program from pyquil import gate_matrices as mat from pyquil.gates import * +from pyquil.operator_estimation import plusX, minusZ from pyquil.paulis import sX, sY, sZ from pyquil.unitary_tools import qubit_adjacent_lifted_gate, program_unitary, lifted_gate_matrix, \ - lifted_gate, lifted_pauli + lifted_gate, lifted_pauli, lifted_state_operator def test_random_gates(): @@ -326,3 +327,40 @@ def test_lifted_pauli(): true_matrix[0, 0] = 2 true_matrix[-1, -1] = -2 np.testing.assert_allclose(trial_matrix, true_matrix) + + +def test_lifted_state_operator(): + xz_state = plusX(5) * minusZ(6) + + plus = np.array([1, 1]) / np.sqrt(2) + plus = plus[:, np.newaxis] + proj_plus = plus @ plus.conj().T + assert proj_plus.shape == (2, 2) + + one = np.array([0, 1]) + one = one[:, np.newaxis] + proj_one = one @ one.conj().T + assert proj_one.shape == (2, 2) + + np.testing.assert_allclose( + np.kron(proj_one, proj_plus), + lifted_state_operator(xz_state, qubits=[5, 6]), + ) + + +def test_lifted_state_operator_backwards_qubits(): + xz_state = plusX(5) * minusZ(6) + plus = np.array([1, 1]) / np.sqrt(2) + plus = plus[:, np.newaxis] + proj_plus = plus @ plus.conj().T + assert proj_plus.shape == (2, 2) + + one = np.array([0, 1]) + one = one[:, np.newaxis] + proj_one = one @ one.conj().T + assert proj_one.shape == (2, 2) + + np.testing.assert_allclose( + np.kron(proj_plus, proj_one), + lifted_state_operator(xz_state, qubits=[6, 5]), + ) diff --git a/pyquil/unitary_tools.py b/pyquil/unitary_tools.py index 1dc19f254..5b8c68d73 100644 --- a/pyquil/unitary_tools.py +++ b/pyquil/unitary_tools.py @@ -17,7 +17,8 @@ import numpy as np -from pyquil.gate_matrices import SWAP, QUANTUM_GATES +from pyquil.gate_matrices import SWAP, QUANTUM_GATES, STATES +from pyquil.operator_estimation import TensorProductState from pyquil.paulis import PauliSum, PauliTerm from pyquil.quilbase import Gate @@ -55,7 +56,11 @@ def qubit_adjacent_lifted_gate(i, matrix, n_qubits): In general, this takes a k-qubit gate (2D matrix 2^k x 2^k) and lifts it to the complete Hilbert space of dim 2^num_qubits, as defined by - the rightward tensor product (1) in arXiv:1608.03355. + the right-to-left tensor product (1) in arXiv:1608.03355. + + Developer note: Quil and the QVM like qubits to be ordered such that qubit 0 is on the right. + Therefore, in ``qubit_adjacent_lifted_gate``, ``lifted_pauli``, and ``lifted_state_operator``, + we build up the lifted matrix by performing the kronecker product from right to left. Note that while the qubits are addressed in decreasing order, starting with num_qubit - 1 on the left and ending with qubit 0 on the @@ -311,6 +316,10 @@ def lifted_pauli(pauli_sum: Union[PauliSum, PauliTerm], qubits: List[int]): [ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [ 0.+2.j, 0.+0.j, 0.+0.j, 0.+0.j]] + Developer note: Quil and the QVM like qubits to be ordered such that qubit 0 is on the right. + Therefore, in ``qubit_adjacent_lifted_gate``, ``lifted_pauli``, and ``lifted_state_operator``, + we build up the lifted matrix by performing the kronecker product from right to left. + :param pauli_sum: Pauli representation of an operator :param qubits: list of qubits in the order they will be represented in the resultant matrix. :returns: matrix representation of the pauli_sum operator @@ -345,3 +354,24 @@ def tensor_up(pauli_sum: Union[PauliSum, PauliTerm], qubits: List[int]): :returns: matrix representation of the pauli_sum operator """ return lifted_pauli(pauli_sum=pauli_sum, qubits=qubits) + + +def lifted_state_operator(state: TensorProductState, qubits: List[int]): + """Take a TensorProductState along with a list of qubits and return a matrix + corresponding to the tensored-up representation of the states' density operator form. + + Developer note: Quil and the QVM like qubits to be ordered such that qubit 0 is on the right. + Therefore, in ``qubit_adjacent_lifted_gate``, ``lifted_pauli``, and ``lifted_state_operator``, + we build up the lifted matrix by using the *left* kronecker product. + + :param state: The state + :param qubits: list of qubits in the order they will be represented in the resultant matrix. + """ + mat = 1.0 + for qubit in qubits: + oneq_state = state[qubit] + assert oneq_state.qubit == qubit + state_vector = STATES[oneq_state.label][oneq_state.index][:, np.newaxis] + state_matrix = state_vector @ state_vector.conj().T + mat = np.kron(state_matrix, mat) + return mat