diff --git a/pygsti/objectivefns/wildcardbudget.py b/pygsti/objectivefns/wildcardbudget.py index e7a74722b..e2206e7a1 100644 --- a/pygsti/objectivefns/wildcardbudget.py +++ b/pygsti/objectivefns/wildcardbudget.py @@ -67,7 +67,7 @@ def to_vector(self): def from_vector(self, w_vec): """ - Set the parameters of this wildcard budge. + Set the parameters of this wildcard budget. Parameters ---------- @@ -462,88 +462,83 @@ def update_probs(self, probs_in, probs_out, freqs, layout, precomp=None, probs_f return p_deriv if return_deriv else None -class PrimitiveOpsWildcardBudget(WildcardBudget): +class PrimitiveOpsWildcardBudgetBase(WildcardBudget): """ - A wildcard budget containing one parameter per "primitive operation". + A base class for wildcard budget objects that allocate wildcard error per "primitive operation". - A parameter's absolute value gives the amount of "slack", or - "wildcard budget" that is allocated per that particular primitive - operation. + The amount of wildcard error for a primitive operation gives the amount of "slack" + that is allocated per that particular operation. Primitive operations are the components of circuit layers, and so - the wilcard budget for a circuit is just the sum of the (abs vals of) - the parameters corresponding to each primitive operation in the circuit. + the wilcard budget for a circuit is just the sum of the wildcard errors + corresponding to each primitive operation in the circuit. Parameters ---------- - primitive_op_labels : iterable or dict + primitive_op_labels : iterable A list of primitive-operation labels, e.g. `Label('Gx',(0,))`, which give all the possible primitive ops (components of circuit - layers) that will appear in circuits. Each one of these operations - will be assigned it's own independent element in the wilcard-vector. - A dictionary can be given whose keys are Labels and whose values are - 0-based parameter indices. In the non-dictionary case, each label gets - it's own parameter. Dictionaries allow multiple labels to be associated - with the *same* wildcard budget parameter, - e.g. `{Label('Gx',(0,)): 0, Label('Gy',(0,)): 0}`. - If `'SPAM'` is included as a primitive op, this value correspond to a - uniform "SPAM budget" added to each circuit. + layers) that will appear in circuits. - start_budget : float or dict, optional - An initial value to set all the parameters to (if a float), or a - dictionary mapping primitive operation labels to initial values. + wildcard_vector : numpy.ndarray + An initial wildcard vector of the parameter values of this budget. + + idle_name : str, optional + The gate name to be used for the 1-qubit idle gate. If not `None`, then + circuit budgets are computed by considering layers of the circuit as being + "padded" with `1-qubit` idles gates on any empty lines. """ - def __init__(self, primitive_op_labels, start_budget=0.0, idle_name=None): - """ - Create a new PrimitiveOpsWildcardBudget. + def __init__(self, primitive_op_labels, wildcard_vector, idle_name=None): + #if isinstance(primitive_op_labels, dict): + # assert(set(primitive_op_labels.values()) == set(range(len(set(primitive_op_labels.values()))))) + # self.primOpLookup = primitive_op_labels + #else: + self.primitive_op_labels = primitive_op_labels + self.primitive_op_index = {lbl: i for i, lbl in enumerate(primitive_op_labels)} + self._idlename = idle_name + self.per_op_wildcard_vector = self._per_op_wildcard_error_from_vector(wildcard_vector) - Parameters - ---------- - primitive_op_labels : iterable or dict - A list of primitive-operation labels, e.g. `Label('Gx',(0,))`, - which give all the possible primitive ops (components of circuit - layers) that will appear in circuits. Each one of these operations - will be assigned it's own independent element in the wilcard-vector. - A dictionary can be given whose keys are Labels and whose values are - 0-based parameter indices. In the non-dictionary case, each label gets - it's own parameter. Dictionaries allow multiple labels to be associated - with the *same* wildcard budget parameter, - e.g. `{Label('Gx',(0,)): 0, Label('Gy',(0,)): 0}`. - If `'SPAM'` is included as a primitive op, this value correspond to a - uniform "SPAM budget" added to each circuit. + super(PrimitiveOpsWildcardBudgetBase, self).__init__(wildcard_vector) - start_budget : float or dict, optional - An initial value to set all the parameters to (if a float), or a - dictionary mapping primitive operation labels to initial values. + def _per_op_wildcard_error_from_vector(self, wildcard_vector): + """ + Returns an array of per-operation wildcard errors based on `wildcard_vector`, + with ordering corresponding to `self.primitive_op_labels`. + """ + raise NotImplementedError("Sub-classes need to implement this!") - idle_name : str, optional - The gate name to be used for the 1-qubit idle gate. If not `None`, then - circuit budgets are computed by considering layers of the circuit as being - "padded" with `1-qubit` idles gates on any empty lines. + def _per_op_wildcard_error_deriv_from_vector(self, wildcard_vector): + """ + Returns a MxN matrix, where M=(# of primitive ops) and N=(# of parameters), + such that its (i,j)-th element is the derivative of the wildcard error + for the i-th primitive op with respect to the j-th budget parameter. + """ + raise NotImplementedError("Sub-classes need to implement this!") + def from_vector(self, w_vec): """ - if isinstance(primitive_op_labels, dict): - assert(set(primitive_op_labels.values()) == set(range(len(set(primitive_op_labels.values()))))) - self.primOpLookup = primitive_op_labels - else: - self.primOpLookup = {lbl: i for i, lbl in enumerate(primitive_op_labels)} + Set the parameters of this wildcard budget. - if 'SPAM' in self.primOpLookup: - self.spam_index = self.primOpLookup['SPAM'] - else: - self.spam_index = None + Parameters + ---------- + w_vec : numpy array + A vector of parameter values. - self._idlename = idle_name + Returns + ------- + None + """ + super().from_vector(w_vec) + self.per_op_wildcard_vector = self._per_op_wildcard_error_from_vector(w_vec) - nParams = len(set(self.primOpLookup.values())) - if isinstance(start_budget, dict): - Wvec = _np.zeros(nParams, 'd') - for op, val in start_budget.items: - Wvec[self.primOpLookup[op]] = val - else: - Wvec = _np.array([start_budget] * nParams) - super(PrimitiveOpsWildcardBudget, self).__init__(Wvec) + @property + def num_primitive_ops(self): + return len(self.primitive_op_labels) + + @property + def wildcard_error_per_op(self): + return {lbl: val for lbl, val in zip(self.primitive_op_labels, self.per_op_wildcard_vector)} def circuit_budget(self, circuit): """ @@ -558,17 +553,18 @@ def circuit_budget(self, circuit): ------- float """ + error_per_op = self.wildcard_error_per_op + def budget_for_label(lbl): - if lbl in self.primOpLookup: # Note: includes len(lbl.components) == 0 case of (global) idle - return pos(Wvec[self.primOpLookup[lbl]]) - elif lbl.name in self.primOpLookup: - return pos(Wvec[self.primOpLookup[lbl.name]]) + if lbl in error_per_op: # Note: includes len(lbl.components) == 0 case of (global) idle + return pos(error_per_op[lbl]) + elif lbl.name in error_per_op: + return pos(error_per_op[lbl.name]) else: assert(not lbl.is_simple()), "Simple label %s must be a primitive op of this WEB!" % str(lbl) return sum([budget_for_label(component) for component in lbl.components]) - Wvec = self.wildcard_vector - budget = 0 if (self.spam_index is None) else pos(Wvec[self.spam_index]) + budget = error_per_op.get('SPAM', 0) layers = [circuit.layer_label(i) for i in range(circuit.depth)] if (self._idlename is None) \ else [circuit.layer_label_with_idles(i, idle_gate_name=self._idlename) for i in range(circuit.depth)] for layer in layers: @@ -616,19 +612,19 @@ def precompute_for_same_circuits(self, circuits): object """ def budget_deriv_for_label(lbl): - if lbl in self.primOpLookup: # Note: includes len(lbl.components) == 0 case of (global) idle - deriv = _np.zeros(len(self.wildcard_vector), 'd') - deriv[self.primOpLookup[lbl]] = 1.0 + if lbl in self.primitive_op_index: # Note: includes len(lbl.components) == 0 case of (global) idle + deriv = _np.zeros(self.num_primitive_ops, 'd') + deriv[self.primitive_op_index[lbl]] = 1.0 return deriv - elif lbl.name in self.primOpLookup: - deriv = _np.zeros(len(self.wildcard_vector), 'd') - deriv[self.primOpLookup[lbl.name]] = 1.0 + elif lbl.name in self.primitive_op_index: + deriv = _np.zeros(self.num_primitive_ops, 'd') + deriv[self.primitive_op_index[lbl.name]] = 1.0 return deriv else: assert(not lbl.is_simple()), "Simple label %s must be a primitive op of this WEB!" % str(lbl) return sum([budget_deriv_for_label(component) for component in lbl.components]) - circuit_budget_matrix = _np.zeros((len(circuits), len(self.wildcard_vector)), 'd') + circuit_budget_matrix = _np.zeros((len(circuits), self.num_primitive_ops), 'd') for i, circuit in enumerate(circuits): layers = [circuit.layer_label(i) for i in range(circuit.depth)] if (self._idlename is None) \ @@ -636,10 +632,11 @@ def budget_deriv_for_label(lbl): for layer in layers: circuit_budget_matrix[i, :] += budget_deriv_for_label(layer) - if self.spam_index is not None: - circuit_budget_matrix[:, self.spam_index] = 1.0 + if 'SPAM' in self.primitive_op_index: + circuit_budget_matrix[:, self.primitive_op_index['SPAM']] = 1.0 - return circuit_budget_matrix + dprimop_dwvec = self._per_op_wildcard_error_deriv_from_vector(self.wildcard_vector) + return _np.dot(circuit_budget_matrix, dprimop_dwvec) @property def description(self): @@ -655,11 +652,12 @@ def description(self): Keys are primitive op labels and values are (description_string, value) tuples. """ wildcardDict = {} - for lbl, index in self.primOpLookup.items(): + error_per_op = self.wildcard_error_per_op + for lbl, val in error_per_op.items(): if lbl == "SPAM": continue # treated separately below - wildcardDict[lbl] = ('budget per each instance %s' % str(lbl), pos(self.wildcard_vector[index])) - if self.spam_index is not None: - wildcardDict['SPAM'] = ('uniform per-circuit SPAM budget', pos(self.wildcard_vector[self.spam_index])) + wildcardDict[lbl] = ('budget per each instance %s' % str(lbl), pos(val)) + if 'SPAM' in error_per_op: + wildcardDict['SPAM'] = ('uniform per-circuit SPAM budget', pos(error_per_op['SPAM'])) return wildcardDict def budget_for(self, op_label): @@ -678,11 +676,10 @@ def budget_for(self, op_label): ------- float """ - return pos(self.wildcard_vector[self.primOpLookup[op_label]]) + return pos(self.per_op_wildcard_vector[self.primitive_op_index[op_label]]) def __str__(self): - wildcardDict = {lbl: pos(self.wildcard_vector[index]) for lbl, index in self.primOpLookup.items()} - return "Wildcard budget: " + str(wildcardDict) + return "Wildcard budget: " + str(self.wildcard_error_per_op) #For these helper functions, see Robin's notes @@ -914,3 +911,201 @@ def update_circuit_probs(probs, freqs, circuit_budget): #assert(_np.isclose(W, compTVD)), "TVD mismatch!" return updated_qvec + + +class PrimitiveOpsWildcardBudget(PrimitiveOpsWildcardBudgetBase): + """ + A wildcard budget containing one parameter per "primitive operation". + + A parameter's absolute value gives the amount of "slack", or + "wildcard budget" that is allocated per that particular primitive + operation. + + Primitive operations are the components of circuit layers, and so + the wilcard budget for a circuit is just the sum of the (abs vals of) + the parameters corresponding to each primitive operation in the circuit. + + Parameters + ---------- + primitive_op_labels : iterable or dict + A list of primitive-operation labels, e.g. `Label('Gx',(0,))`, + which give all the possible primitive ops (components of circuit + layers) that will appear in circuits. Each one of these operations + will be assigned it's own independent element in the wilcard-vector. + A dictionary can be given whose keys are Labels and whose values are + 0-based parameter indices. In the non-dictionary case, each label gets + it's own parameter. Dictionaries allow multiple labels to be associated + with the *same* wildcard budget parameter, + e.g. `{Label('Gx',(0,)): 0, Label('Gy',(0,)): 0}`. + If `'SPAM'` is included as a primitive op, this value correspond to a + uniform "SPAM budget" added to each circuit. + + start_budget : float or dict, optional + An initial value to set all the parameters to (if a float), or a + dictionary mapping primitive operation labels to initial values. + """ + + def __init__(self, primitive_op_labels, start_budget=0.0, idle_name=None): + """ + Create a new PrimitiveOpsWildcardBudget. + + Parameters + ---------- + primitive_op_labels : iterable or dict + A list of primitive-operation labels, e.g. `Label('Gx',(0,))`, + which give all the possible primitive ops (components of circuit + layers) that will appear in circuits. Each one of these operations + will be assigned it's own independent element in the wilcard-vector. + A dictionary can be given whose keys are Labels and whose values are + 0-based parameter indices. In the non-dictionary case, each label gets + it's own parameter. Dictionaries allow multiple labels to be associated + with the *same* wildcard budget parameter, + e.g. `{Label('Gx',(0,)): 0, Label('Gy',(0,)): 0}`. + If `'SPAM'` is included as a primitive op, this value correspond to a + uniform "SPAM budget" added to each circuit. + + start_budget : float or dict, optional + An initial value to set all the parameters to (if a float), or a + dictionary mapping primitive operation labels to initial values. + + idle_name : str, optional + The gate name to be used for the 1-qubit idle gate. If not `None`, then + circuit budgets are computed by considering layers of the circuit as being + "padded" with `1-qubit` idles gates on any empty lines. + + """ + if isinstance(primitive_op_labels, dict): + num_params = len(set(primitive_op_labels.values())) + assert(set(primitive_op_labels.values()) == set(range(num_params))) + + prim_op_lbls = list(primitive_op_labels.keys()) + self.primitive_op_param_index = primitive_op_labels + else: + num_params = len(primitive_op_labels) + prim_op_lbls = primitive_op_labels + self.primitive_op_param_index = {lbl: i for i, lbl in enumerate(primitive_op_labels)} + + self.trivial_param_mapping = all([self.primitive_op_param_index[lbl] == i + for i, lbl in enumerate(prim_op_lbls)]) + + #generate initial wildcard vector + if isinstance(start_budget, dict): + Wvec = _np.zeros(num_params, 'd') + for op, val in start_budget.items(): + Wvec[self.primitive_op_param_index[op]] = val + else: + Wvec = _np.array([start_budget] * num_params) + + super().__init__(prim_op_lbls, Wvec, idle_name) + + def _per_op_wildcard_error_from_vector(self, wildcard_vector): + """ + Returns an array of per-operation wildcard errors based on `wildcard_vector`, + with ordering corresponding to `self.primitive_op_labels`. + """ + if self.trivial_param_mapping: + return wildcard_vector + else: + return _np.array([self.primitive_op_param_index[lbl] for lbl in self.primitive_op_labels]) + + def _per_op_wildcard_error_deriv_from_vector(self, wildcard_vector): + """ + Returns a MxN matrix, where M=(# of primitive ops) and N=(# of parameters), + such that its (i,j)-th element is the derivative of the wildcard error + for the i-th primitive op with respect to the j-th budget parameter. + """ + if self.trivial_param_mapping: + return _np.identity(self.num_params) + else: + ret = _np.zeros((self.num_primitive_ops, self.num_params), 'd') + for i, lbl in enumerate(self.primitive_op_labels): + ret[i, self.primitive_op_param_index[lbl]] = 1.0 + return ret + + +class PrimitiveOpsSingleScaleWildcardBudget(PrimitiveOpsWildcardBudgetBase): + """ + A wildcard budget containing a single scaling parameter. + + This type of wildcard budget has a single parameter, and sets the wildcard + error of each primitive op to be this scale parameter multiplied by a fixed + reference value for the primitive op. + + Typically, the reference values are chosen to be a modeled metric of gate quality, + such as a gate's gauge-optimized diamond distance to its target gate. Then, + once a feasible wildcard error vector is found, the scaling parameter is + the fractional increase of the metric (e.g. the diamond distance) of each + primitive op needed to reconcile the model and data. + + Parameters + ---------- + primitive_op_labels : list + A list of primitive-operation labels, e.g. `Label('Gx',(0,))`, + which give all the possible primitive ops (components of circuit + layers) that will appear in circuits. + + reference_values : list, optional + A list of the reference values for each primitive op, in the same order as + `primitive_op_list`. + + alpha : float, optional + The initial value of the single scaling parameter that multiplies the reference + values of each op to give the wildcard error that op. + + idle_name : str, optional + The gate name to be used for the 1-qubit idle gate. If not `None`, then + circuit budgets are computed by considering layers of the circuit as being + "padded" with `1-qubit` idles gates on any empty lines. + """ + def __init__(self, primitive_op_labels, reference_values=None, alpha=0, idle_name=None, reference_name=None): + if reference_values is None: + reference_values = [1.0] * len(primitive_op_labels) + self.reference_values = _np.array(reference_values, 'd') + self.reference_name = reference_name + + Wvec = _np.array([alpha], 'd') + super().__init__(primitive_op_labels, Wvec, idle_name) + + @property + def alpha(self): + return self.wildcard_vector[0] + + @alpha.setter + def alpha(self, val): + self.from_vector(_np.array([val], 'd')) + + def _per_op_wildcard_error_from_vector(self, wildcard_vector): + """ + Returns an array of per-operation wildcard errors based on `wildcard_vector`, + with ordering corresponding to `self.primitive_op_labels`. + """ + alpha = wildcard_vector[0] + return self.reference_values * alpha + + def _per_op_wildcard_error_deriv_from_vector(self, wildcard_vector): + """ + Returns a MxN matrix, where M=(# of primitive ops) and N=(# of parameters), + such that its (i,j)-th element is the derivative of the wildcard error + for the i-th primitive op with respect to the j-th budget parameter. + """ + return self.reference_values.reshape((self.num_primitive_ops, 1)).copy() + + @property + def description(self): + """ + A dictionary of quantities describing this budget. + + Return the contents of this budget in a dictionary containing + (description, value) pairs for each element name. + + Returns + ------- + dict + Keys are primitive op labels and values are (description_string, value) tuples. + """ + wildcardDict = super().description + + #Add an entry for the alpha parameter. + wildcardDict['alpha'] = ('Fraction of diamond distance needed as wildcard error', self.alpha) + + return wildcardDict diff --git a/pygsti/optimize/wildcardopt.py b/pygsti/optimize/wildcardopt.py index 57051838e..977fc1f03 100644 --- a/pygsti/optimize/wildcardopt.py +++ b/pygsti/optimize/wildcardopt.py @@ -423,6 +423,75 @@ def _proxy_agg_dlogl_hessian(x, tvds, fn0s, percircuit_budget_deriv): return H +def _get_percircuit_budget_deriv(budget, layout): + """ Returns local_percircuit_budget_deriv, global_percircuit_budget_deriv """ + percircuit_budget_deriv = budget.precompute_for_same_circuits(layout.circuits) # for *local* circuits + + #Note: maybe we could do this gather in 1 call (?), but we play it safe and do it col-by-col + global_percircuit_budget_deriv_cols = [] + for i in range(percircuit_budget_deriv.shape[1]): + global_percircuit_budget_deriv_cols.append( + layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) + return percircuit_budget_deriv, _np.column_stack(global_percircuit_budget_deriv_cols) + + +def optimize_wildcard_bisect_alpha(budget, objfn, two_dlogl_threshold, redbox_threshold, printer, + guess=0.1, tol=1e-3): + printer.log("Beginning wildcard budget optimization using alpha bisection method.") + + layout = objfn.layout + critical_percircuit_budgets = _get_critical_circuit_budgets(objfn, redbox_threshold) # for *global* circuits + percircuit_budget_deriv, global_percircuit_budget_deriv = _get_percircuit_budget_deriv(budget, layout) + + initial_probs = objfn.probs.copy() + current_probs = initial_probs.copy() + probs_freqs_precomp = budget.precompute_for_same_probs_freqs(initial_probs, objfn.freqs, layout) + + def is_feasible(x): + budget.from_vector(x) + budget.update_probs(initial_probs, current_probs, objfn.freqs, objfn.layout, percircuit_budget_deriv, + probs_freqs_precomp) + f0 = _np.array([_agg_dlogl(current_probs, objfn, two_dlogl_threshold)]) + fi = critical_percircuit_budgets - _np.dot(global_percircuit_budget_deriv, x) + return _np.all(_np.concatenate((f0, fi)) <= 0) # All constraints must be negative to be feasible + + left = None + right = None + + while left is None or right is None: + printer.log(f'Searching for interval [{left}, {right}] with guess {guess}', 2) + # Test for feasibility + if is_feasible(_np.array([guess], 'd')): + printer.log('Guess value is feasible, ', 2) + left = guess + guess = left / 2 + else: + printer.log('Guess value is infeasible, ', 2) + right = guess + guess = 2 * right + printer.log('Interval found!', 2) + + # We now have an interval containing the crossover point + # Perform bisection + while abs(left - right) > tol: + printer.log(f'Performing bisection on interval [{left}, {right}]', 2) + test = left - (left - right) / 2.0 + + if is_feasible(_np.array([test], 'd')): + # Feasible, so shift left down + printer.log('Test value is feasible, ', 2) + left = test + else: + printer.log('Test value is infeasible, ', 2) + right = test + + printer.log('Interval within tolerance!', 2) + + budget.from_vector(_np.array([left], 'd')) # set budget to the feasible one + printer.log(f'Optimized value of alpha = {left}') + return + + def optimize_wildcard_budget_cvxopt(budget, L1weights, objfn, two_dlogl_threshold, redbox_threshold, printer, abs_tol=1e-5, rel_tol=1e-5, max_iters=50): """Uses CVXOPT to optimize the wildcard budget. Includes both aggregate and per-circuit constraints.""" @@ -439,14 +508,7 @@ def optimize_wildcard_budget_cvxopt(budget, L1weights, objfn, two_dlogl_threshol initial_probs = objfn.probs.copy() # *local* current_probs = initial_probs.copy() - percircuit_budget_deriv = budget.precompute_for_same_circuits(layout.circuits) # for *local* circuits - - #Note: maybe we could do this gather in 1 call (?), but we play it safe and do it col-by-col - global_percircuit_budget_deriv_cols = [] - for i in range(percircuit_budget_deriv.shape[1]): - global_percircuit_budget_deriv_cols.append( - layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) - global_percircuit_budget_deriv = _np.column_stack(global_percircuit_budget_deriv_cols) + percircuit_budget_deriv, global_percircuit_budget_deriv = _get_percircuit_budget_deriv(budget, layout) critical_percircuit_budgets = _get_critical_circuit_budgets(objfn, redbox_threshold) # for *global* circuits critical_percircuit_budgets.shape = (len(critical_percircuit_budgets), 1) @@ -524,14 +586,7 @@ def optimize_wildcard_budget_cvxopt_zeroreg(budget, L1weights, objfn, two_dlogl_ initial_probs = objfn.probs.copy() current_probs = initial_probs.copy() - percircuit_budget_deriv = budget.precompute_for_same_circuits(layout.circuits) - - #Note: maybe we could do this gather in 1 call (?), but we play it safe and do it col-by-col - global_percircuit_budget_deriv_cols = [] - for i in range(percircuit_budget_deriv.shape[1]): - global_percircuit_budget_deriv_cols.append( - layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) - global_percircuit_budget_deriv = _np.column_stack(global_percircuit_budget_deriv_cols) + percircuit_budget_deriv, global_percircuit_budget_deriv = _get_percircuit_budget_deriv(budget, layout) critical_percircuit_budgets = _get_critical_circuit_budgets(objfn, redbox_threshold) critical_percircuit_budgets.shape = (len(critical_percircuit_budgets), 1) @@ -615,15 +670,8 @@ def optimize_wildcard_budget_barrier(budget, L1weights, objfn, two_dlogl_thresho # for increasing values of t until 1/t <= epsilon (precision tolerance) printer.log("Beginning wildcard budget optimization using a barrier method.") layout = objfn.layout - percircuit_budget_deriv = budget.precompute_for_same_circuits(layout.circuits) # for *local* circuits critical_percircuit_budgets = _get_critical_circuit_budgets(objfn, redbox_threshold) # for *global* circuits - - #Note: maybe we could do this gather in 1 call (?), but we play it safe and do it col-by-col - global_percircuit_budget_deriv_cols = [] - for i in range(percircuit_budget_deriv.shape[1]): - global_percircuit_budget_deriv_cols.append( - layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) - global_percircuit_budget_deriv = _np.column_stack(global_percircuit_budget_deriv_cols) + percircuit_budget_deriv, global_percircuit_budget_deriv = _get_percircuit_budget_deriv(budget, layout) x0 = budget.to_vector() initial_probs = objfn.probs.copy() @@ -892,15 +940,7 @@ def optimize_wildcard_budget_cvxopt_smoothed(budget, L1weights, objfn, two_dlogl #initial_probs = objfn.probs.copy() #current_probs = initial_probs.copy() - percircuit_budget_deriv = budget.precompute_for_same_circuits(layout.circuits) - - #Note: maybe we could do this gather in 1 call (?), but we play it safe and do it col-by-col - global_percircuit_budget_deriv_cols = [] - for i in range(percircuit_budget_deriv.shape[1]): - global_percircuit_budget_deriv_cols.append( - layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) - global_percircuit_budget_deriv = _np.column_stack(global_percircuit_budget_deriv_cols) - + percircuit_budget_deriv, global_percircuit_budget_deriv = _get_percircuit_budget_deriv(budget, layout) critical_percircuit_budgets = _get_critical_circuit_budgets(objfn, redbox_threshold) critical_percircuit_budgets.shape = (len(critical_percircuit_budgets), 1) num_circuits = len(layout.circuits) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index f8b4676db..fb28857c3 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -524,11 +524,62 @@ class GSTBadFitOptions(_NicelySerializable): GST fit is considered satisfactory (and no "bad-fit" processing is needed). actions : tuple, optional - Actions to take when a GST fit is unsatisfactory. + Actions to take when a GST fit is unsatisfactory. Allowed actions include: + - 'wildcard': Find an admissable wildcard model... + - 'ddist_wildcard': Fits a single parameter wildcard model in which + the amount of wildcard error added to an operation is proportional + to the diamond distance between that operation and the target. + - 'robust': scale data according out "robust statistics v1" algorithm, + where we drastically scale down (reduce) the data due to especially + poorly fitting circuits. Namely, if a circuit's log-likelihood ratio + exceeds the 95% confidence region about its expected value (the # of + degrees of freedom in the circuits outcomes), then the data is scaled + by the `expected_value / actual_value`, so that the new value exactly + matches what would be expected. Ideally there are only a few of these + "outlier" circuits, which correspond errors in the measurement apparatus. + - 'Robust': same as 'robust', but re-optimize the final objective function + (usually the log-likelihood) after performing the scaling to get the + final estimate. + - 'robust+': scale data according out "robust statistics v2" algorithm, + which performs the v1 algorithm (see 'robust' above) and then further + rescales all the circuit data to achieve the desired chi2 distribution + of per-circuit goodness-of-fit values *without reordering* these values. + - 'Robust+': same as 'robust+', but re-optimize the final objective function + (usually the log-likelihood) after performing the scaling to get the + final estimate. + - 'do nothing': do not perform any additional actions. Used to help avoid + the need for special cases when working with multiple types of bad-fit actions. wildcard_budget_includes_spam : bool, optional Include a SPAM budget within the wildcard budget used to process the `"wildcard"` action. + + wildcard_L1_weights : np.array, optional + An array of weights affecting the L1 penalty term used to select a feasible + wildcard error vector `w_i` that minimizes `sum_i weight_i* |w_i|` (a weighted + L1 norm). Elements of this array must correspond to those of the wildcard budget + being optimized, typically the primitive operations of the estimated model - but + to get the order right you should specify `wildcard_primitive_op_labels` to be sure. + If `None`, then all weights are assumed to be 1. + + wildcard_primitive_op_labels: list, optional + The primitive operation labels used to construct the :class:`PrimitiveOpsWildcardBudget` + that is optimized. If `None`, equal to `model.primitive_op_labels + model.primitive_instrument_labels` + where `model` is the estimated model, with `'SPAM'` at the end if `wildcard_budget_includes_spam` + is True. When specified, should contain a subset of the default values. + + wildcard_methods: tuple, optional + A list of the methods to use to optimize the wildcard error vector. Default is `("neldermead",)`. + Options include `"neldermead"`, `"barrier"`, `"cvxopt"`, `"cvxopt_smoothed"`, `"cvxopt_small"`, + and `"cvxpy_noagg"`. So many methods exist because different convex solvers behave differently + (unfortunately). Leave as the default as a safe option, but `"barrier"` is pretty reliable and much + faster than `"neldermead"`, and is a good option so long as it runs. + + wildcard_inadmissable_action: {"print", "raise"}, optional + What to do when an inadmissable wildcard error vector is found. The default just prints this + information and continues, while `"raise"` raises a `ValueError`. Often you just want this information + printed so that when the wildcard analysis fails in this way it doesn't cause the rest of an analysis + to abort. """ @classmethod @@ -555,8 +606,8 @@ def __init__(self, threshold=DEFAULT_BAD_FIT_THRESHOLD, actions=(), wildcard_budget_includes_spam=True, wildcard_L1_weights=None, wildcard_primitive_op_labels=None, wildcard_initial_budget=None, wildcard_methods=('neldermead',), - wildcard_inadmissable_action='print'): - valid_actions = ('wildcard', 'Robust+', 'Robust', 'robust+', 'robust', 'do nothing') + wildcard_inadmissable_action='print', wildcard1d_reference='diamond distance'): + valid_actions = ('wildcard', 'wildcard1d', 'Robust+', 'Robust', 'robust+', 'robust', 'do nothing') if not all([(action in valid_actions) for action in actions]): raise ValueError("Invalid action in %s! Allowed actions are %s" % (str(actions), str(valid_actions))) self.threshold = float(threshold) @@ -567,6 +618,7 @@ def __init__(self, threshold=DEFAULT_BAD_FIT_THRESHOLD, actions=(), self.wildcard_initial_budget = wildcard_initial_budget self.wildcard_methods = wildcard_methods self.wildcard_inadmissable_action = wildcard_inadmissable_action # can be 'raise' or 'print' + self.wildcard1d_reference = wildcard1d_reference def _to_nice_serialization(self): state = super()._to_nice_serialization() @@ -577,7 +629,8 @@ def _to_nice_serialization(self): 'primitive_op_labels': self.wildcard_primitive_op_labels, 'initial_budget': self.wildcard_initial_budget, # serializable? 'methods': self.wildcard_methods, - 'indadmissable_action': self.wildcard_inadmissable_action}, + 'indadmissable_action': self.wildcard_inadmissable_action, + '1d_reference': self.wildcard1d_reference}, }) return state @@ -590,7 +643,8 @@ def _from_nice_serialization(cls, state): # memo holds already de-serialized ob wildcard.get('primitive_op_labels', None), wildcard.get('initial_budget', None), tuple(wildcard.get('methods', ['neldermead'])), - wildcard.get('inadmissable_action', 'print')) + wildcard.get('inadmissable_action', 'print'), + wildcard.get('1d_reference', 'diamond distance')) class GSTObjFnBuilders(_NicelySerializable): @@ -1835,6 +1889,32 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, # new_params['unmodeled_error'] = None continue # no need to add a new estimate - we just update the base estimate + elif badfit_typ == 'wildcard1d': + + #If this estimate is the target model then skip adding the diamond distance wildcard. + if base_estimate_label != 'Target': + try: + budget = _compute_wildcard_budget_1d_model(base_estimate, objfn_cache, mdc_objfn, parameters, + badfit_options, printer - 1) + + base_estimate.extra_parameters['wildcard1d' + "_unmodeled_error"] = budget + base_estimate.extra_parameters['wildcard1d' + "_unmodeled_active_constraints"] \ + = None + + base_estimate.extra_parameters["unmodeled_error"] = budget + base_estimate.extra_parameters["unmodeled_active_constraints"] = None + except NotImplementedError as e: + printer.warning("Failed to get wildcard budget - continuing anyway. Error was:\n" + str(e)) + new_params['unmodeled_error'] = None + #except AssertionError as e: + # printer.warning("Failed to get wildcard budget - continuing anyway. Error was:\n" + str(e)) + # new_params['unmodeled_error'] = None + continue # no need to add a new estimate - we just update the base estimate + + else: + printer.log('Diamond distance wildcard model is incompatible with the Target estimate, skipping.', 3) + continue + elif badfit_typ == "do nothing": continue # go to next on-bad-fit directive @@ -1870,6 +1950,101 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, gauge_opt_params.copy(), go_gs_final, gokey, comm, printer - 1) +def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, parameters, badfit_options, verbosity): + """ + Create a wildcard budget for a model estimate. This version of the function produces a wildcard estimate + using the model introduced by Tim and Stefan in the RCSGST paper. + TODO: docstring (update) + + Parameters + ---------- + model : Model + The model to add a wildcard budget to. + + ds : DataSet + The data the model predictions are being compared with. + + circuits_to_use : list + The circuits whose data are compared. + + parameters : dict + Various parameters of the estimate at hand. + + badfit_options : GSTBadFitOptions, optional + Options specifying what post-processing actions should be performed when + a fit is unsatisfactory. Contains detailed parameters for wildcard budget + creation. + + comm : mpi4py.MPI.Comm, optional + An MPI communicator used to run this computation in parallel. + + mem_limit : int, optional + A rough per-processor memory limit in bytes. + + verbosity : int, optional + Level of detail printed to stdout. + + Returns + ------- + PrimitiveOpsWildcardBudget + """ + printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, mdc_objfn.resource_alloc) + badfit_options = GSTBadFitOptions.cast(badfit_options) + model = mdc_objfn.model + ds = mdc_objfn.dataset + global_circuits_to_use = mdc_objfn.global_circuits + + printer.log("******************* Adding Wildcard Budget **************************") + + #cache construction code. + # Extract model, dataset, objfn, etc. + # Note: must evaluate mdc_objfn *before* passing to wildcard fn init so internal probs are init + mdc_objfn.fn(model.to_vector()) + + ## PYGSTI TRANSPLANT from pygsti.protocols.gst._compute_wildcard_budget + # Compute the various thresholds + ds_dof = ds.degrees_of_freedom(global_circuits_to_use) + nparams = model.num_modeltest_params # just use total number of params + percentile = 0.025 + nboxes = len(global_circuits_to_use) + + two_dlogl_threshold = _chi2.ppf(1 - percentile, max(ds_dof - nparams, 1)) + redbox_threshold = _chi2.ppf(1 - percentile / nboxes, 1) + + ref, reference_name = _compute_1d_reference_values_and_name(estimate, badfit_options) + primitive_ops = list(ref.keys()) + wcm = _wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[k] for k in primitive_ops], + reference_name=reference_name) + _opt.wildcardopt.optimize_wildcard_bisect_alpha(wcm, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, + guess=0.1, tol=1e-3) # results in optimized wcm + return wcm + + +def _compute_1d_reference_values_and_name(estimate, badfit_options): + final_model = estimate.models['final iteration estimate'] + target_model = estimate.models['target'] + gaugeopt_model = _alg.gaugeopt_to_target(final_model, target_model) + + if badfit_options.wildcard1d_reference == 'diamond distance': + dd = {} + for key, op in gaugeopt_model.operations.items(): + dd[key] = 0.5 * _tools.diamonddist(op.to_dense(), target_model.operations[key].to_dense()) + + spamdd = {} + for key, op in gaugeopt_model.preps.items(): + spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), + _tools.vec_to_stdmx(target_model.preps[key].to_dense(), 'pp')) + + for key in gaugeopt_model.povms.keys(): + spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + + dd['SPAM'] = sum(spamdd.values()) + return dd, 'diamond distance' + else: + raise ValueError("Invalid wildcard1d_reference value (%s) in bad-fit options!" + % str(badfit_options.wildcard1d_reference)) + + def _compute_robust_scaling(scale_typ, objfn_cache, mdc_objfn): """ Get the per-circuit data scaling ("weights") for a given type of robust-data-scaling. @@ -2044,7 +2219,7 @@ def _compute_wildcard_budget(objfn_cache, mdc_objfn, parameters, badfit_options, L1weights = _np.ones(budget.num_params) if badfit_options.wildcard_L1_weights: for op_label, weight in badfit_options.wildcard_L1_weights.items(): - L1weights[budget.primOpLookup[op_label]] = weight + L1weights[budget.primitive_op_param_index[op_label]] = weight printer.log("Using non-uniform L1 weights: " + str(list(L1weights))) # Note: must evaluate mdc_objfn *before* passing to wildcard fn init so internal probs are init @@ -2195,7 +2370,7 @@ def _evaluate_constraints(wv): # Note: active_constraints_list is typically stored in parameters['unmodeled_error active constraints'] # of the relevant Estimate object. primOp_labels = _collections.defaultdict(list) - for lbl, i in budget.primOpLookup.items(): primOp_labels[i].append(str(lbl)) + for lbl, i in budget.primitive_op_param_index.items(): primOp_labels[i].append(str(lbl)) for i, active_constraints in enumerate(active_constraints_list): if active_constraints: printer.log("** ACTIVE constraints for " + "--".join(primOp_labels[i]) + " **") diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index 451bb027e..11f6ca3c5 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -35,6 +35,7 @@ from pygsti.baseobjs.label import Label as _Lbl from pygsti.baseobjs.verbosityprinter import VerbosityPrinter as _VerbosityPrinter from pygsti.tools.legacytools import deprecate as _deprecated_fn +from pygsti.objectivefns.wildcardbudget import PrimitiveOpsSingleScaleWildcardBudget #maybe import these from drivers.longsequence so they stay synced? ROBUST_SUFFIX_LIST = [".robust", ".Robust", ".robust+", ".Robust+"] # ".wildcard" (not a separate estimate anymore) @@ -1212,6 +1213,13 @@ def construct_standard_report(results, title="auto", flags.add('ShowScaling') if est.parameters.get('unmodeled_error', None): flags.add('ShowUnmodeledError') + #check if the wildcard budget is an instance + #of the diamond distance model, in which case we + #will add an extra flag/plot to the report. + if (isinstance(est.parameters['unmodeled_error'], PrimitiveOpsSingleScaleWildcardBudget) + and est.parameters['unmodeled_error'].reference_name == 'diamond distance'): + flags.add('DiamondDistanceWildcard') + if combine_robust: flags.add('CombineRobust') diff --git a/pygsti/report/section/goodness.py b/pygsti/report/section/goodness.py index 0bec644a4..60c8b41e9 100644 --- a/pygsti/report/section/goodness.py +++ b/pygsti/report/section/goodness.py @@ -136,6 +136,10 @@ class GoodnessUnmodeledSection(_Section): def unmodeled_error_budget_table(workspace, switchboard=None, **kwargs): return workspace.WildcardBudgetTable(switchboard.wildcard_budget) + @_Section.figure_factory(1) + def unmodeled_error_ddist_bar_plot(workspace, switchboard=None, **kwargs): + return workspace.WildcardSingleScaleBarPlot(switchboard.wildcard_budget, reference_name='Diamond Distance') + @_Section.figure_factory(4) def final_model_fit_progress_bar_plot_ume(workspace, switchboard=None, max_lengths=None, comm=None, **kwargs): return workspace.FitComparisonBarPlot( diff --git a/pygsti/report/templates/standard_html_report/tabs/Goodness_unmodeled.html b/pygsti/report/templates/standard_html_report/tabs/Goodness_unmodeled.html index fac26cd9d..4f8fc494b 100644 --- a/pygsti/report/templates/standard_html_report/tabs/Goodness_unmodeled.html +++ b/pygsti/report/templates/standard_html_report/tabs/Goodness_unmodeled.html @@ -21,6 +21,14 @@