-
Notifications
You must be signed in to change notification settings - Fork 55
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Feature ddist wildcard integration #267
Changes from 8 commits
5ecd081
fd66810
7050f37
f56f2a2
22392c8
90e69d8
7d4497f
de5911e
f296ef6
8103406
1fcf5ce
1171298
bb00901
57de6d9
89fcd2b
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -43,6 +43,9 @@ | |
from pygsti.tools.legacytools import deprecate as _deprecated_fn | ||
|
||
|
||
|
||
|
||
|
||
#For results object: | ||
ROBUST_SUFFIX_LIST = [".robust", ".Robust", ".robust+", ".Robust+"] | ||
DEFAULT_BAD_FIT_THRESHOLD = 2.0 | ||
|
@@ -524,11 +527,28 @@ 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+': ... | ||
- 'robust+': ... | ||
- 'robust': ... | ||
- 'do nothing': ... | ||
|
||
wildcard_budget_includes_spam : bool, optional | ||
Include a SPAM budget within the wildcard budget used to process | ||
the `"wildcard"` action. | ||
|
||
wildcard_L1_weights : | ||
|
||
wildcard_primitive_op_labels: | ||
|
||
wildcard_methods: tuple, optional (default ('neldermead',)) | ||
|
||
wildcard_inadmissable_action: str, optional (default 'print') | ||
|
||
""" | ||
|
||
@classmethod | ||
|
@@ -556,7 +576,7 @@ def __init__(self, threshold=DEFAULT_BAD_FIT_THRESHOLD, actions=(), | |
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') | ||
valid_actions = ('wildcard', 'ddist_wildcard', '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) | ||
|
@@ -1834,6 +1854,31 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, | |
# 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 | ||
|
||
elif badfit_typ == 'ddist_wildcard': | ||
|
||
#If this estimate is the target model then skip adding the diamond distance wildcard. | ||
if base_estimate_label != 'Target': | ||
try: | ||
budget = _compute_wildcard_budget_ddist_model(base_estimate,objfn_cache, mdc_objfn, parameters, badfit_options, printer - 1) | ||
|
||
base_estimate.extra_parameters['ddist_wildcard' + "_unmodeled_error"] = budget | ||
base_estimate.extra_parameters['ddist_wildcard' + "_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 this.',3) | ||
continue | ||
|
||
elif badfit_typ == "do nothing": | ||
continue # go to next on-bad-fit directive | ||
|
@@ -1869,6 +1914,182 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, | |
results.estimates[base_estimate_label + '.' + badfit_typ].add_gaugeoptimized( | ||
gauge_opt_params.copy(), go_gs_final, gokey, comm, printer - 1) | ||
|
||
def _compute_wildcard_budget_ddist_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) | ||
|
||
# Not in cache, recompute | ||
ddists = _compute_ddists(estimate) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just point for discussion for tomorrow's dev meeting: I could imagine the user wanting to do this sort of scaling single-param wildcard with something other than the diamond distance. Is it worth plumbing that in now? Could be as simple as adding a kwarg to this function, and we can leave the rest of the plumbing for when we need it? |
||
|
||
critical_percircuit_budgets = _opt.wildcardopt._get_critical_circuit_budgets(mdc_objfn, redbox_threshold) | ||
|
||
#In Stefan's original code these were all passed in as a dictionary, instead pass these in as arguments. | ||
|
||
alpha = _bisect_alpha(0.1, 1e-3, mdc_objfn, two_dlogl_threshold, | ||
redbox_threshold, ddists, critical_percircuit_budgets) | ||
|
||
return _ddist_wildcard_model(alpha, ddists) | ||
|
||
def _bisect_alpha(guess, tol, mdc_objfn, two_dlogl_threshold, redbox_threshold, | ||
ddists, critical_percircuit_budgets): | ||
|
||
left = None | ||
right = None | ||
|
||
while left is None or right is None: | ||
print(f'Searching for interval [{left}, {right}] with guess {guess}') | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. We probably want to make this just debugger/info logging rather than print all the time. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Good point, just pushed a commit implementing that. |
||
# Test for feasibility | ||
guessval = _test_feasible_alphas([guess], mdc_objfn, two_dlogl_threshold, | ||
redbox_threshold, ddists, critical_percircuit_budgets) | ||
if guessval[0]: | ||
print('Guess value is feasible, ', end='') | ||
left = guess | ||
guess = left/2 | ||
else: | ||
print('Guess value is infeasible, ', end='') | ||
right = guess | ||
guess = 2*right | ||
print('Interval found!') | ||
|
||
# We now have an interval containing the crossover point | ||
# Perform bisection | ||
while abs(left - right) > tol: | ||
print(f'Performing bisection on interval [{left}, {right}]') | ||
test = left - (left - right)/2.0 | ||
|
||
testval = _test_feasible_alphas([test], mdc_objfn, two_dlogl_threshold, | ||
redbox_threshold, ddists, critical_percircuit_budgets) | ||
|
||
if testval[0]: | ||
# Feasible, so shift left down | ||
print(f'Test value is feasible, ', end='') | ||
left = test | ||
else: | ||
print(f'Test value is infeasible, ', end='') | ||
right = test | ||
|
||
print('Interval within tolerance!') | ||
print() | ||
|
||
return left # Return the feasible one | ||
|
||
def _test_feasible_alphas(alphas, mdc_objfn, two_dlogl_threshold, redbox_threshold, ddists, critical_percircuit_budgets): | ||
|
||
if not isinstance(alphas, list): | ||
alphas = [alphas] | ||
|
||
# Iterate through alphas and compute wildcard feasibilities | ||
feasible = [] | ||
for alpha in alphas: | ||
wcm = _ddist_wildcard_model(alpha, ddists) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just a note: we could avoid creating a new budget object on each iteration here if we wanted/needed the performance. I think this is fine for now - just an opportunity for later on if performance becomes an issue. |
||
|
||
## PYGSTI TRANSPLANT from pygsti.optimize.wildcardopt.optimize_wildcard_budget_barrier | ||
percircuit_budget_deriv = wcm.precompute_for_same_circuits(mdc_objfn.layout.circuits) | ||
|
||
global_percircuit_budget_deriv_cols = [] | ||
for i in range(percircuit_budget_deriv.shape[1]): | ||
global_percircuit_budget_deriv_cols.append( | ||
mdc_objfn.layout.allgather_local_array('c', percircuit_budget_deriv[:, i])) | ||
global_percircuit_budget_deriv = _np.column_stack(global_percircuit_budget_deriv_cols) | ||
|
||
initial_probs = mdc_objfn.probs.copy() | ||
current_probs = initial_probs.copy() | ||
probs_freqs_precomp = wcm.precompute_for_same_probs_freqs(initial_probs, mdc_objfn.freqs, mdc_objfn.layout) | ||
|
||
# Feasibility calculation, similar to penalty_vec function | ||
x = wcm.to_vector() | ||
wcm.update_probs(initial_probs, current_probs, mdc_objfn.freqs, mdc_objfn.layout, percircuit_budget_deriv, | ||
probs_freqs_precomp) | ||
f0 = _np.array([_opt.wildcardopt._agg_dlogl(current_probs, mdc_objfn, two_dlogl_threshold)]) | ||
fi = critical_percircuit_budgets - _np.dot(global_percircuit_budget_deriv, x) | ||
|
||
feasible.append(_np.all(_np.concatenate((f0, fi)) <= 0)) # All constraints must be negative to be feasible | ||
|
||
return feasible | ||
|
||
def _compute_ddists(estimate): | ||
final_model = estimate.models['final iteration estimate'] | ||
target_model = estimate.models['target'] | ||
gaugeopt_model = _alg.gaugeopt_to_target(final_model, target_model) | ||
|
||
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 | ||
|
||
def _ddist_wildcard_model(alpha, dd): | ||
keys = list(dd.keys()) | ||
wcb_values = _np.array([alpha * dd[key] for key in keys]) | ||
wcb_dict= {key:value for key,value in zip(keys,wcb_values)} | ||
wcm = _wild.PrimitiveOpsDiamondDistanceWildcardBudget(keys, start_budget=wcb_dict, ddists=dd, alpha=alpha) | ||
return wcm | ||
|
||
def _compute_robust_scaling(scale_typ, objfn_cache, mdc_objfn): | ||
""" | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe '1d_wildcard'?