1. pyGSTi
  2. jupyter_notebooks
  3. Examples

How to generate error bars for 2Q-GST

In [1]:
import pygsti
import pickle
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 file, run the 2Q-GST example)
with open("example_files/easy_2q_results.pkl","rb") as f:
    results = pickle.load(f)

#Set a memory limit
print("Memory limit was = ", results.estimates['default'].parameters.get('memLimit',"none given"))
results.estimates['default'].parameters['memLimit'] = 2.5*(1024.0)**3 # 2.5GB
print("Memory limit is now = ", results.estimates['default'].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['default'].add_confidence_region_factory('go0', '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.48GB
 mem(1 subtrees, 1,1 param-grps, 1 proc-grps) in 0s = 6773.17GB (6773.17GB fc)
Created evaluation tree with 1 subtrees.  Will divide 1 procs into 1 (subtree-processing)
 groups of ~1 procs each, to distribute over (1616,1616) params (taken as 1616,4 param groups of ~1,404 params).
 Memory estimate = 2.08GB (cache=1317, wrtLen1=1, wrtLen2=404, subsPerProc=1).
rank 0: 29.3325s: block 0/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 58.4175s: block 1/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 87.3489s: block 2/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 116.28s: block 3/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 146.016s: block 4/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 175.125s: block 5/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 204.208s: block 6/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 233.866s: block 7/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 262.991s: block 8/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 292.132s: block 9/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 321.038s: block 10/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 349.834s: block 11/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 378.618s: block 12/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 407.43s: block 13/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 436.11s: block 14/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 464.901s: block 15/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 493.786s: block 16/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 522.456s: block 17/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 551.266s: block 18/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 579.995s: block 19/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 608.664s: block 20/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 637.449s: block 21/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 666.115s: block 22/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 694.955s: block 23/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 723.691s: block 24/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 752.808s: block 25/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 781.978s: block 26/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 811.188s: block 27/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 840.21s: block 28/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 869.157s: block 29/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 898.371s: block 30/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 927.436s: block 31/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 957s: block 32/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 986.078s: block 33/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1015.71s: block 34/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1044.65s: block 35/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1073.76s: block 36/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1102.93s: block 37/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1132.23s: block 38/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1161.42s: block 39/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1191.09s: block 40/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1220.73s: block 41/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1250.07s: block 42/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1279.96s: block 43/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1309.22s: block 44/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1339.77s: block 45/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1368.74s: block 46/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1397.34s: block 47/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1426.26s: block 48/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1455.08s: block 49/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1483.91s: block 50/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1512.55s: block 51/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1541.71s: block 52/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1571.52s: block 53/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1601.31s: block 54/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1631.02s: block 55/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1660.61s: block 56/4043, sub-tree 0/1, sub-tree-len = 1317
rank 0: 1690.12s: block 57/4043, sub-tree 0/1, sub-tree-len = 1317

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-2-7208f561b6ac> in <module>()
      7 # initialize a factory for the 'go0' gauge optimization within the 'default' estimate
      8 crfact = results.estimates['default'].add_confidence_region_factory('go0', 'final')
----> 9 crfact.compute_hessian(comm=comm) #optionally use multiple processors
     10 crfact.project_hessian('intrinsic error')
     11 

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

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

~/research/pyGSTi/packages/pygsti/objects/matrixforwardsim.py in bulk_hprobs_by_block(self, evalTree, wrtSlicesList, bReturnDProbs12, comm)
   2705             #Fill arrays
   2706             self._fill_result_tuple((None, dprobs1, dprobs2, hprobs), evalTree,
-> 2707                                     slice(None), slice(None), calc_and_fill)
   2708             hProdCache = hGs = dProdCache2 = dGs2 =  None # free mem
   2709             if bReturnDProbs12:

~/research/pyGSTi/packages/pygsti/objects/forwardsim.py in _fill_result_tuple(self, result_tup, evalTree, param_slice1, param_slice2, calc_and_fill_fn)
    524             #          all of the raw operation sequences which need to be computed
    525             #          for the current spamTuple (this list has the SAME length as fInds).
--> 526             calc_and_fill_fn(spamTuple,fInds,gInds,pslc1,pslc2,False) #TODO: remove SumInto == True cases
    527 
    528         return

~/research/pyGSTi/packages/pygsti/objects/matrixforwardsim.py in calc_and_fill(spamTuple, fInds, gInds, pslc1, pslc2, sumInto)
   2652             if deriv2MxToFill is not None:
   2653                 _fas(deriv2MxToFill, [fInds,pslc2], self._dprobs_from_rhoE( 
-> 2654                     spamTuple, rho, E, Gs[gInds], dGs2[gInds], scaleVals[gInds], wrtSlice2), add=sumInto)
   2655 
   2656             _fas(mxToFill, [fInds,pslc1,pslc2], self._hprobs_from_rhoE( 

~/research/pyGSTi/packages/pygsti/objects/matrixforwardsim.py in _dprobs_from_rhoE(self, spamTuple, rho, E, Gs, dGs, scaleVals, wrtSlice)
   1677         # dp_dOps      = squeeze( dot( E, dot( dGs, rho ) ), axis=(0,3))
   1678         old_err2 = _np.seterr(invalid='ignore', over='ignore')
-> 1679         dp_dOps = _np.squeeze( _np.dot( E, _np.dot( dGs, rho ) ), axis=(0,3) ) * scaleVals[:,None]
   1680         _np.seterr(**old_err2)
   1681            # may overflow, but OK ; shape == (len(circuit_list), nDerivCols)

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 file
with open("example_files/easy_2q_results_withCI.pkl","wb") as f:
    pickle.dump(results, f)
In [ ]: