How to generate error bars for 2Q-GST

In [1]:
import pygsti
import time

#If we were using MPI
# from mpi4py import MPI
# comm = MPI.COMM_WORLD
comm = None

#Load the 2-qubit results (if you don't have this directory, run the 2Q-GST example)
results = pygsti.io.load_results_from_dir("example_files/My2QExample", "StandardGST")

#Set a memory limit
print("Memory limit was = ", results.estimates['CPTP'].parameters.get('memLimit',"none given"))
results.estimates['CPTP'].parameters['memLimit'] = 2.5*(1024.0)**3 # 2.5GB
print("Memory limit is now = ", results.estimates['CPTP'].parameters['memLimit'])
Memory limit was =  3221225472
Memory limit is now =  2684354560.0
In [2]:
# error bars in reports require the presence of a fully-initialized
# "confidence region factory" within the relevant Estimate object.
# In most cases "fully-initialized" means that a Hessian has been 
# computed and projected onto the non-gauge space.
start = time.time()

# initialize a factory for the 'go0' gauge optimization within the 'default' estimate
crfact = results.estimates['CPTP'].add_confidence_region_factory('stdgaugeopt', 'final')
crfact.compute_hessian(comm=comm) #optionally use multiple processors
crfact.project_hessian('intrinsic error')

end = time.time()
print("Total time=%f hours" % ((end - start) / 3600.0))
Evaltree generation (deriv) w/mem limit = 2.47GB
bulk_evaltree: created initial tree (1082 strs) in 1s
bulk_evaltree: split tree (1 subtrees) in 0s
 mem(1 subtrees, 1,1 param-grps, 1 proc-grps) in 1s = 7853.61GB (7853.61GB fc)
Created evaluation tree with 1 subtrees.  Will divide 1 procs into 1 (subtree-processing)
 groups of ~1 procs each, to distribute over (1920,1920) params (taken as 1920,4 param groups of ~1,480 params).
 Memory estimate = 2.03GB (cache=1082, wrtLen1=1, wrtLen2=480, subsPerProc=1).
rank 0: 94.5588s: block 0/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 131.424s: block 1/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 168.311s: block 2/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 204.173s: block 3/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 239.922s: block 4/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 274.424s: block 5/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 309.599s: block 6/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 344.601s: block 7/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 380.99s: block 8/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 416.213s: block 9/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 451.373s: block 10/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 486.383s: block 11/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 522.758s: block 12/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 557.807s: block 13/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 592.904s: block 14/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 628.007s: block 15/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 664.472s: block 16/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 699.313s: block 17/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 734.498s: block 18/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 769.855s: block 19/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 806.698s: block 20/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 842s: block 21/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 877.752s: block 22/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 913.219s: block 23/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 949.306s: block 24/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 984.357s: block 25/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1019.19s: block 26/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1054.28s: block 27/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1090.41s: block 28/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1126.91s: block 29/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1162.49s: block 30/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1197.98s: block 31/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1234.56s: block 32/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1269.9s: block 33/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1305.11s: block 34/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1340.37s: block 35/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1376.68s: block 36/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1412.1s: block 37/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1447.11s: block 38/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1482.88s: block 39/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1519.37s: block 40/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1554.64s: block 41/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1590.87s: block 42/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1628.39s: block 43/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1664.2s: block 44/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1698.84s: block 45/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1733.53s: block 46/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1768.03s: block 47/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1803.87s: block 48/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1838.46s: block 49/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1873.87s: block 50/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 1910.03s: block 51/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 3837.1s: block 52/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 3872.04s: block 53/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 3908.06s: block 54/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 3942.97s: block 55/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 3979.41s: block 56/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4015.35s: block 57/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4050.74s: block 58/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4086.53s: block 59/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4123.42s: block 60/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4159.3s: block 61/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4194.81s: block 62/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4231.65s: block 63/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4270.02s: block 64/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4305.87s: block 65/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4341.45s: block 66/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4376.67s: block 67/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4413.29s: block 68/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4448.5s: block 69/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4483.95s: block 70/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4518.95s: block 71/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4555.33s: block 72/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4590.38s: block 73/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4625.86s: block 74/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4661.12s: block 75/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4697.37s: block 76/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4732.63s: block 77/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4767.74s: block 78/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4802.98s: block 79/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4839.52s: block 80/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4874.94s: block 81/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4909.99s: block 82/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4945.25s: block 83/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 4981.63s: block 84/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5016.7s: block 85/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5051.81s: block 86/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5086.85s: block 87/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5123.41s: block 88/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5158.31s: block 89/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5193.36s: block 90/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5228.24s: block 91/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5264.67s: block 92/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5300.09s: block 93/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5335.06s: block 94/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5370.76s: block 95/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5407.25s: block 96/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5442.45s: block 97/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5477.7s: block 98/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5512.84s: block 99/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5549.23s: block 100/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5584.49s: block 101/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5620.42s: block 102/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5655.66s: block 103/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5692.2s: block 104/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5727.45s: block 105/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5762.68s: block 106/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5798.16s: block 107/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5834.89s: block 108/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5870.06s: block 109/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5906.05s: block 110/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5941.48s: block 111/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 5977.28s: block 112/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6011.99s: block 113/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6047.07s: block 114/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6082.19s: block 115/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6118.64s: block 116/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6154.54s: block 117/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6189.9s: block 118/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6225.23s: block 119/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6261.84s: block 120/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6297.34s: block 121/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6332.55s: block 122/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6367.79s: block 123/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6404.29s: block 124/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6440.17s: block 125/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6475.38s: block 126/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6510.54s: block 127/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6547.22s: block 128/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6582.93s: block 129/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6618.17s: block 130/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6653.65s: block 131/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6690.67s: block 132/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6725.87s: block 133/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6761.19s: block 134/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6796.28s: block 135/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6832.49s: block 136/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6867.64s: block 137/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6903.53s: block 138/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6938.71s: block 139/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 6975.65s: block 140/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 7010.96s: block 141/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 7045.97s: block 142/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 7081.56s: block 143/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 7118.21s: block 144/4803, sub-tree 0/1, sub-tree-len = 1082
rank 0: 7154.8s: block 145/4803, sub-tree 0/1, sub-tree-len = 1082
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-2-dda12a7d9471> in <module>()
      7 # initialize a factory for the 'go0' gauge optimization within the 'default' estimate
      8 crfact = results.estimates['CPTP'].add_confidence_region_factory('stdgaugeopt', 'final')
----> 9 crfact.compute_hessian(comm=comm) #optionally use multiple processors
     10 crfact.project_hessian('intrinsic error')
     11 

~/pyGSTi/pygsti/objects/confidenceregionfactory.py in compute_hessian(self, comm, memLimit, approximate)
    238                                  minProbClip, probClipInterval, radius,
    239                                  comm=comm, memLimit=memLimit, verbosity=vb,
--> 240                                  opLabelAliases=aliases)
    241 
    242             nonMarkRadiusSq = max(2 * (_tools.logl_max(model, dataset)

~/pyGSTi/pygsti/tools/likelihoodfns.py in logl_hessian(model, dataset, circuit_list, minProbClip, probClipInterval, radius, poissonPicture, check, comm, memLimit, opLabelAliases, smartc, verbosity)
    853         k, kmax = 0, len(mySliceTupList)
    854         for (slice1, slice2, hprobs, dprobs12) in model.bulk_hprobs_by_block(
--> 855                 evalSubTree, mySliceTupList, True, blkComm):
    856             rank = comm.Get_rank() if (comm is not None) else 0
    857 

~/pyGSTi/pygsti/objects/matrixforwardsim.py in bulk_hprobs_by_block(self, evalTree, wrtSlicesList, bReturnDProbs12, comm)
   2702             hProdCache = self._compute_hproduct_cache(
   2703                 evalTree, prodCache, dProdCache1, dProdCache2,
-> 2704                 scaleCache, comm, wrtSlice1, wrtSlice2)
   2705             hGs = evalTree.final_view(hProdCache, axis=0)
   2706 

~/pyGSTi/pygsti/objects/matrixforwardsim.py in _compute_hproduct_cache(self, evalTree, prodCache, dProdCache1, dProdCache2, scaleCache, comm, wrtSlice1, wrtSlice2)
   1123 
   1124             dLdRa = _np.swapaxes(_np.dot(dL1, dR2), 1, 2)
-> 1125             dLdRb = _np.swapaxes(_np.dot(dL2, dR1), 1, 2)
   1126             dLdR_sym = dLdRa + _np.swapaxes(dLdRb, 0, 1)
   1127 

KeyboardInterrupt: 

Note above cell was executed for demonstration purposes, and was keyboard-interrupted intentionally since it would have taken forever on a single processor.

In [ ]:
#write results back to disk
results.write()
In [ ]: