One other way we can use the pygsti.report.reportables module described in the ModelAnalysisMetrics tutorial is to procedurally generate error bars for any quantity you want.
First, let's simulate a noisy GST experiment
import pygsti
from pygsti.modelpacks import smq1Q_XY
from pygsti.report import reportables as rptbl, modelfunction as modelfn
target_model = smq1Q_XY.target_model()
L=128
edesign = smq1Q_XY.create_gst_experiment_design(L)
noisy_model = target_model.randomize_with_unitary(.1)
noisy_model = noisy_model.depolarize(.05)
N=64
dataset = pygsti.data.simulate_data(noisy_model,edesign,N)
gst_proto = pygsti.protocols.StandardGST(modes=['full TP','CPTPLND','Target'],verbosity=2)
data = pygsti.protocols.ProtocolData(edesign,dataset)
results = gst_proto.run(data)
Now let's compute error bars on the CPTP estimate, and then get a 95% confidence interval "view" from the ConfidenceRegionFactory.
crfact = results.estimates['CPTPLND'].add_confidence_region_factory('stdgaugeopt', 'final')
crfact.compute_hessian(comm=None, mem_limit=3.0*(1024.0)**3) #optionally use multiple processors & set memlimit
crfact.project_hessian('intrinsic error')
crf_view = results.estimates['CPTPLND'].confidence_region_factories['stdgaugeopt','final'].view(95)
Finally, we can construct pygsti.report.ModelFunction objects that take a function which computes some observable from a model and the extracted view from above to compute error bars on that quantity of interest.
One common thing to check is error bars on the process matrices. The ModelFunction in this case only needs to return the operation:
final_model = results.estimates['CPTPLND'].models['stdgaugeopt'].copy()
def get_op(model, lbl):
return model[lbl]
get_op_modelfn = modelfn.modelfn_factory(get_op)
rptbl.evaluate(get_op_modelfn(final_model, ("Gxpi2", 0)), crf_view)
rptbl.evaluate(get_op_modelfn(final_model, ("Gypi2", 0)), crf_view)
But we can also create model functions that perform more complicated actions, such as computing other reportables.
# Note that when creating ModelFunctions in this way, the model where you want the quantity evaluated must be the first argument
def ddist(model, ideal_model, lbl, basis):
return rptbl.half_diamond_norm(model[lbl], ideal_model[lbl], basis)
ddist_modelfn = modelfn.modelfn_factory(ddist)
rptbl.evaluate(ddist_modelfn(final_model, target_model, ("Gxpi2", 0), 'pp'), crf_view)
rptbl.evaluate(ddist_modelfn(final_model, target_model, ("Gypi2", 0), 'pp'), crf_view)