diff --git a/pygsti/algorithms/dfe.py b/pygsti/algorithms/dfe.py index a8e4d0641..3a8267634 100644 --- a/pygsti/algorithms/dfe.py +++ b/pygsti/algorithms/dfe.py @@ -13,17 +13,25 @@ from pygsti.tools import internalgates as _gates from pygsti.algorithms import randomcircuit as _rc -def sample_dfe_circuit(pspec, circuit, clifford_compilations, seed=None): +def sample_dfe_circuit(pspec, circuit, clifford_compilations, create_same_measurement_reference=False, seed=None, + fixed_identity_locations=None): """ - ... ... ... ... + create_same_measurement_reference: If true returns another DFE circuit with the same measurement but a trivial circuit + in the middle """ qubit_labels = circuit.line_labels n = len(qubit_labels) rand_state = _np.random.RandomState(seed) # Ok if seed is None - - rand_pauli, rand_sign, pauli_circuit = _rc._sample_random_pauli(n = n, pspec = pspec, qubit_labels=qubit_labels, - absolute_compilation = clifford_compilations, - circuit = True, include_identity = False) + + if fixed_identity_locations is None: + + rand_pauli, rand_sign, pauli_circuit = _rc._sample_random_pauli(n = n, pspec = pspec, qubit_labels=qubit_labels, + absolute_compilation = clifford_compilations, + circuit = True, include_identity = False) + + else: + rand_pauli, rand_sign, pauli_circuit = _rc._sample_random_pauli_with_identities_fixed(n = n, is_identity = fixed_identity_locations, pspec = pspec, qubit_labels=qubit_labels, + absolute_compilation = clifford_compilations, circuit = True) s_inputstate, p_inputstate, s_init_layer, p_init_layer, prep_circuit = _rc._sample_stabilizer(rand_pauli, rand_sign, clifford_compilations, qubit_labels) @@ -66,4 +74,17 @@ def sample_dfe_circuit(pspec, circuit, clifford_compilations, seed=None): outcircuit = full_circuit - return outcircuit, measurement, sign \ No newline at end of file + return outcircuit, measurement, sign + + + # def create_dfe_reference_circuit_from_measurement_and_sign(measurement, sign, qubit_labels, pspec, clifford_compilations): + + # pauli_circuit = _circuit_from_pauli_and_sign(measurement, sign, pspec = pspec, absolute_compilation = clifford_compilations) + # s_inputstate, p_inputstate, s_init_layer, p_init_layer, prep_circuit = _rc._sample_stabilizer((measurement, sign, clifford_compilations, qubit_labels) + + # s_pc, p_pc = _symp.symplectic_rep_of_clifford_circuit(pauli_circuit, pspec=pspec.subset(gate_names_to_include='all', qubit_labels_to_keep=qubit_labels)) + + # # Turn the stabilizer into an all Z and I stabilizer + # s_stab, p_stab, stab_circuit = _rc._stabilizer_to_all_zs(measurement, qubit_labels, clifford_compilations) + + # return \ No newline at end of file diff --git a/pygsti/algorithms/randomcircuit.py b/pygsti/algorithms/randomcircuit.py index de8bf9ad9..e54d39651 100644 --- a/pygsti/algorithms/randomcircuit.py +++ b/pygsti/algorithms/randomcircuit.py @@ -908,6 +908,95 @@ def create_random_circuit(pspec, length, qubit_labels=None, sampler='Qeliminatio circuit.done_editing() return circuit + +def create_random_circuit_with_fixed_feature_values(pspec, width, depth, xi, qubit_labels=None, rand_state=None): + """ + todo + """ + if rand_state is None: + rand_state = _np.random.RandomState() + if qubit_labels is not None: + assert(isinstance(qubit_labels, list) or isinstance(qubit_labels, tuple)), "SubsetQs must be a list or a tuple!" + qubits = list(qubit_labels[:]) # copy this list + assert(len(qubits) == width) + else: + qubits = list(pspec.qubit_labels[:width]) # copy this list + + # Initialize an empty circuit, to populate with sampled layers. + circuit = _cir.Circuit(layer_labels=[], line_labels=qubits, editable=True) + + #edgelists = [] + possible_twoQgate_locations = [] + for i in range(depth): + # Prep the sampling variables. + #sampled_layer = [] + edgelist = pspec.compute_2Q_connectivity().edges() + edgelist = [e for e in edgelist if all([q in qubits for q in e])] + selectededges = [] + #print(edgelist) + # Go through until all qubits have been assigned a gate. + while len(edgelist) > 0: + + edge = edgelist[rand_state.randint(0, len(edgelist))] + selectededges.append(edge) + # Delete all edges containing these qubits. + edgelist = [e for e in edgelist if not any([q in e for q in edge])] + + #print(edgelist) + + possible_twoQgate_locations += [(i,e) for e in selectededges] + + target_num_2q_gates = xi * width * depth / 2 + + # We can only have an integer number of two-qubit gates, so we round the target up or down to the nearest integer. + # We do so randomly so that the expected number of two-qubit gates in the circuit is the target value. + ceiling_target_num_2q_gates = int(_np.ceil(target_num_2q_gates)) + floor_target_num_2q_gates= int(_np.floor(target_num_2q_gates)) + prob_ceiling = target_num_2q_gates - floor_target_num_2q_gates + num_2q_gates = rand_state.choice([ceiling_target_num_2q_gates, floor_target_num_2q_gates], p=[prob_ceiling, 1 - prob_ceiling]) + + #print(possible_twoQgate_locations) + #print(num_2q_gates) + selected_2q_gate_locations_indices = rand_state.choice([i for i in range(len(possible_twoQgate_locations))], size=num_2q_gates, replace=False) + selected_2q_gate_locations = [possible_twoQgate_locations[i] for i in selected_2q_gate_locations_indices] + + #print('selected', selected_2q_gate_locations) + selected_2q_gates_by_layer = [[] for i in range(depth)] + for (d, edge) in selected_2q_gate_locations: + #print(d, edge) + selected_2q_gates_by_layer[d].append(edge) + + for i in range(depth): + #print(i) + sampled_layer = [] + selected_edges_for_layer = selected_2q_gates_by_layer[i] + used_qubits = [q for edge in selected_edges_for_layer for q in edge] + unusedqubits = [q for q in qubits if q not in used_qubits] + #print(selected_edges_for_layer) + #print('used', used_qubits) + #print('unused', unusedqubits) + + ops_on_qubits = pspec.compute_ops_on_qubits() + for edge in selected_edges_for_layer: + possibleops = ops_on_qubits[edge] + #print(possibleops) + gate_label = possibleops[rand_state.randint(0, len(possibleops))] + sampled_layer.append(gate_label) + #print(selectededges) + #print(sampled_layer) + + for q in unusedqubits: + possibleops = ops_on_qubits[(q,)] + gate_label = possibleops[rand_state.randint(0, len(possibleops))] + sampled_layer.append(gate_label) + + circuit.insert_layer_inplace(sampled_layer, i) + + + circuit.done_editing() + return circuit + + #### Commented out as this code has not been tested since a much older version of pyGSTi and it is probably #### not being used. # def sample_simultaneous_random_circuit(pspec, length, structure='1Q', sampler='Qelimination', samplerargs=[], @@ -3204,10 +3293,6 @@ def _sample_random_pauli(n,pspec = None, absolute_compilation = None, qubit_labe # - qubit_labels: # - circuit: Boolean that determines if a list of single-qubit Paulis or a compiled circuit is returned. - if circuit is True: - if qubit_labels is not None: qubits = qubit_labels[:] # copy this list - else: qubits = pspec.qubit_labels[:] - pauli_list = ['I','X','Y','Z'] if include_identity is False: @@ -3225,17 +3310,56 @@ def _sample_random_pauli(n,pspec = None, absolute_compilation = None, qubit_labe if circuit is False: return pauli, sign else: - pauli_layer_std_lbls = [_lbl.Label(pauli_list[rand_ints[q]], (qubits[q],)) for q in range(n)] - # Converts the layer to a circuit, and changes to the native model. - pauli_circuit = _cir.Circuit(layer_labels=pauli_layer_std_lbls, line_labels=qubits).parallelize() - pauli_circuit = pauli_circuit.copy(editable=True) - pauli_circuit.change_gate_library(absolute_compilation) - if pauli_circuit.depth == 0: - pauli_circuit.insert_layer_inplace([_lbl.Label(())], 0) - pauli_circuit.done_editing() + pauli_circuit = _circuit_from_pauli_and_sign(pauli, sign, pspec=pspec, absolute_compilation=absolute_compilation, qubit_labels=qubit_labels) + return pauli, sign, pauli_circuit + + +def _sample_random_pauli_with_identities_fixed(n, is_identity, pspec = None, absolute_compilation = None, qubit_labels = None, circuit = False): + # Samples a random Pauli along with a +-1 phase. Returns the Pauli as a list or as a circuit depending + # upon the value of "circuit" + # - n: Number of qubits + # - pspec: Processor spec + # - absolute_compilation: compilation rules + # - qubit_labels: + # - circuit: Boolean that determines if a list of single-qubit Paulis or a compiled circuit is returned. + + non_identity_pauli_list = ['X','Y','Z'] + rand_ints = _np.random.randint(0, 3, n) + pauli = [] + for i in range(n): + if is_identity[i]: + pauli.append('I') + else: + pauli.append(non_identity_pauli_list[rand_ints[i]]) + if set(pauli) != set('I'): sign = _np.random.choice([-1,1]) + else: sign = 1 + + if circuit is False: + return pauli, sign + else: + pauli_circuit = _circuit_from_pauli_and_sign(pauli, sign, pspec=pspec, absolute_compilation=absolute_compilation, qubit_labels=qubit_labels) return pauli, sign, pauli_circuit + +def _circuit_from_pauli_and_sign(pauli, sign, pspec = None, absolute_compilation = None, qubit_labels = None): + + pauli_list = ['I','X','Y','Z'] + if qubit_labels is not None: qubits = qubit_labels[:] # copy this list + else: qubits = pspec.qubit_labels[:] + n = len(pauli) + + pauli_layer_std_lbls = [_lbl.Label(pauli[q], (qubits[q],)) for q in range(n)] + # Converts the layer to a circuit, and changes to the native model. + pauli_circuit = _cir.Circuit(layer_labels=pauli_layer_std_lbls, line_labels=qubits).parallelize() + pauli_circuit = pauli_circuit.copy(editable=True) + pauli_circuit.change_gate_library(absolute_compilation) + if pauli_circuit.depth == 0: + pauli_circuit.insert_layer_inplace([_lbl.Label(())], 0) + pauli_circuit.done_editing() + + return pauli_circuit + def _select_neg_evecs(pauli, sign): # Selects the entries in an n-qubit that will be turned be given a -1 1Q eigenstates diff --git a/pygsti/tools/__init__.py b/pygsti/tools/__init__.py index 4fcbf81d4..7e10e7220 100644 --- a/pygsti/tools/__init__.py +++ b/pygsti/tools/__init__.py @@ -34,3 +34,5 @@ from .slicetools import * from .symplectic import * from .typeddict import TypedDict + +from .vbtools import * diff --git a/pygsti/tools/vbtools.py b/pygsti/tools/vbtools.py new file mode 100644 index 000000000..99b967770 --- /dev/null +++ b/pygsti/tools/vbtools.py @@ -0,0 +1,215 @@ +""" +Utility functions for featuremetric and volumetric benchmarking +""" +#*************************************************************************************************** +# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights +# in this software. +# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +# in compliance with the License. You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. +#*************************************************************************************************** + +import numpy as _np +from scipy.stats.qmc import Sobol as _Sobol + +def max_depth_log2(error_per_layer, min_fidelity): + """ + A function for choose a maximum depth for benchmarking + experiments based on a guess of the error per circuit + layer. + + Parameters + ---------- + error_per_layer : float + The error per circuit layer. + + min_fidelity : float + The desired minimum fidelity of the circuits + + Returns + ------- + The smallest integer k such that the fidelity of + a depth 2^k will be less than or equal to `min_fidelity` + and a depth 2^{k-1} circuit will have a fidelity of greater + than `min_fidelity`, if each layer in the circuit has an + error per layer of `error_per_layer` and the noise is uniform + depolarization. If that k is less than 0, the function returns + -1. + + This function uses the approximation that the fidelities of + depolarizing channels multiple, which is a good approximation + except for very few qubit (narrow) circuits. + """ + return max([-1,int(_np.ceil(_np.log2(_np.log(min_fidelity) / _np.log(1 - error_per_layer))))]) + + +def max_depths_log2(widths, error_per_layer, min_fidelity, include_all_widths=True): + """ + A function for choose a maximum depth for benchmarking + experiments, versus circuit width, based on a guess of + the error per circuit layer.i + + Parameters + ---------- + widths : list or array + The widths at which to compute a maximum depth + + error_per_layer : dict + A dictionary containing the error per circuit layer + (values) versus circuit width (keys) + + min_fidelity : float + The desired minimum fidelity of the circuits. + + Returns + ------- + dict + + A dictionary where the keys are the integers w in widths, + and the values are the smallest integer k such that the fidelity of + a depth 2^k will be less than or equal to `min_fidelity` + and a depth 2^{k-1} width-w circuit will have a fidelity of greater + than `min_fidelity`, if each layer in the circuit has an + error per layer of `error_per_layer` and the noise is uniform + depolarization. + + This function uses the approximation that the fidelities of + depolarizing channels multiple, which is a good approximation + except for very few qubit (narrow) circuits. + """ + + if include_all_widths: + out = {w: max([0, max_depth_log2(error_per_layer[w], min_fidelity)]) for w in widths} + else: + out = {} + for w in widths: + mdlog2 = max_depth_log2(error_per_layer[w], min_fidelity) + if mdlog2 >= 0: + out[w] = mdlog2 + return out + + +def sample_circuit_shapes_logspace_uniform(max_depths_log2, num_samples, min_depth_log2=0): + """ + Samples circuit shapes (widths and depths) so that they are uniformly + spaced in log depth space. + + Parameters + ---------- + max_depths_log2 : dict + The maximum depth (value) at each width (key) + + num_samples: float + The number of circuit shapes to sample + + Returns + ------- + np.array + The sampled circuit widths + + np.array + The sampled circuit depths + """ + widths = list(max_depths_log2.keys()) + widths_weights = _np.array([max_depths_log2[w] + 1 for w in widths]) + widths_probablities = widths_weights / _np.sum(widths_weights) + sampled_widths = [] + sampled_depths = [] + for i in range(num_samples): + w = _np.random.choice(widths, p=widths_probablities) + log2_sampled_depth = -1 + while log2_sampled_depth < min_depth_log2: + log2_sampled_depth = max_depths_log2[w] * _np.random.rand() + d = int(_np.round(2**log2_sampled_depth)) + + sampled_widths.append(w) + sampled_depths.append(d) + + return sampled_widths, sampled_depths + +def quasi_uniform_feature_vectors(num_samples, num_features, max_values, min_values, logspace, integer_valued): + """ + todo + """ + + sobol = _Sobol(num_features) + rounded_base2_num_samples = int(_np.ceil(_np.log2(num_samples))) + samples = sobol.random_base2(rounded_base2_num_samples) + + rescaled_max_values = _np.zeros((num_features,), float) + rescaled_min_values = _np.zeros((num_features,), float) + for j in range(num_features): + if logspace[j]: + rescaled_max_values[j] = _np.log2(max_values[j]) + rescaled_min_values[j] = _np.log2(min_values[j]) + else: + rescaled_max_values[j] = max_values[j] + rescaled_min_values[j] = min_values[j] + + rescaled_samples = (rescaled_max_values - rescaled_min_values) * samples + rescaled_min_values + + sampled_feature_vectors = [] + + for j in range(num_features): + sampled_feature_value = rescaled_samples[:,j].copy() + if logspace[j]: + sampled_feature_value = 2**sampled_feature_value + if integer_valued[j]: + sampled_feature_value = _np.round(sampled_feature_value) + sampled_feature_value = sampled_feature_value.astype(int) + + sampled_feature_vectors.append(sampled_feature_value) + + return sampled_feature_vectors + + +def sample_circuit_shapes_logspace_quasi_uniform(max_depths_log2, num_samples, min_depth_log2=0): + """ + Samples circuit shapes (widths and depths) so that they are uniformly + spaced in log depth space. + + Parameters + ---------- + max_depths_log2 : dict + The maximum depth (value) at each width (key) + + num_samples: float + The number of circuit shapes to sample + + Returns + ------- + np.array + The sampled circuit widths + + np.array + The sampled circuit depths + """ + widths = list(max_depths_log2.keys()) + widths_weights = _np.array([max_depths_log2[w] + 1 for w in widths]) + widths_probablities = widths_weights / _np.sum(widths_weights) + sampled_widths = [] + sampled_depths = [] + + for i in range(num_samples): + w = _np.random.choice(widths, p=widths_probablities) + sampled_widths.append(w) + + num_each_width = {w: sampled_widths.count(w) for w in widths} + for w in widths: + sobol = _Sobol(1) + rounded_up_log2_num_samples = int(_np.ceil(_np.log2(num_each_width[w]))) + unscaled_ds = sobol.random_base2(rounded_up_log2_num_samples)[:,0] + sampled_depths_for_width[w].append(d) + + return sampled_widths, sampled_depths + + +def sample_random_subset(w, qubits): + """ + Select a uniformly random set of w qubits from `qubits`. Returned + list is in the same order as `qubits`. + """ + indices = list(np.random.choice(np.arange(0,len(qubits)), w, replace=False)) + indices.sort() + return [qubits[i] for i in indices]