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.read_results_from_dir("example_files/My2QExample", "StandardGST")
# 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()
# Choose forward simulator and how computation should be split up. Here are a couple examples:
#results.estimates['CPTP'].models['stdgaugeopt'].sim = pygsti.forwardsims.MapForwardSimulator(num_atoms=100)
results.estimates['CPTP'].models['stdgaugeopt'].sim = pygsti.forwardsims.MatrixForwardSimulator(param_blk_sizes=(30,30))
# 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, mem_limit=3.0*(1024.0)**3) #optionally use multiple processors & set memlimit
crfact.project_hessian('intrinsic error')
end = time.time()
print("Total time=%f hours" % ((end - start) / 3600.0))
Layout creation w/mem limit = 3.00GB MatrixLayout: 1 processors divided into 1 x 1 x 1 (= 1) grid along circuit and parameter directions. 1 atoms, parameter block size limits (30, 30) *** Distributing 1 atoms to 1 atom-processing groups (1 cores) *** More atom-processors than hosts: each host gets ~1 atom-processors Atom-processors already occupy a single node, dividing atom-processor into 1 param-processors. *** Divided 1-host atom-processor (~1 procs) into 1 param-processing groups *** Param-processors already occupy a single node, dividing param-processor into 1 param2-processors. *** Divided 1-host param-processor (~1 procs) into 1 param2-processing groups *** Esimated memory required = 1.9GB rank 0: 8.10623e-06s: beginning atom 1/1, atom-size (#circuits) = 3884 rank 0: 27.3823s: block 1/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 56.8049s: block 2/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 86.7264s: block 3/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 116.186s: block 4/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 147.012s: block 5/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 176.884s: block 6/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 207.426s: block 7/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 237.48s: block 8/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 265.242s: block 9/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 293.074s: block 10/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 321.699s: block 11/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 349.566s: block 12/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 376.808s: block 13/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 404.092s: block 14/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 430.981s: block 15/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 458.278s: block 16/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 486.687s: block 17/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 513.101s: block 18/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 540.426s: block 19/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 568.835s: block 20/4096, atom 1/1, atom-size (#circuits) = 3884 rank 0: 596.736s: block 21/4096, atom 1/1, atom-size (#circuits) = 3884
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-2-b70a459f9a3e> in <module> 11 # initialize a factory for the 'go0' gauge optimization within the 'default' estimate 12 crfact = results.estimates['CPTP'].add_confidence_region_factory('stdgaugeopt', 'final') ---> 13 crfact.compute_hessian(comm=comm, mem_limit=3.0*(1024.0)**3) #optionally use multiple processors & set memlimit 14 crfact.project_hessian('intrinsic error') 15 ~/pyGSTi/pygsti/protocols/confidenceregionfactory.py in compute_hessian(self, comm, mem_limit, approximate) 311 hessian_fn = _tools.logl_approximate_hessian if approximate \ 312 else _tools.logl_hessian --> 313 hessian = hessian_fn(model, dataset, circuit_list, 314 minProbClip, probClipInterval, radius, 315 comm=comm, mem_limit=mem_limit, verbosity=vb, ~/pyGSTi/pygsti/tools/likelihoodfns.py in logl_hessian(model, dataset, circuits, min_prob_clip, prob_clip_interval, radius, poisson_picture, op_label_aliases, mdc_store, comm, mem_limit, verbosity) 421 regularization, {'prob_clip_interval': prob_clip_interval}, 422 op_label_aliases, comm, mem_limit, ('hessian',), (), mdc_store, verbosity) --> 423 local = -obj.hessian() # negative b/c objective is deltaLogL = max_logl - logL 424 return obj.layout.allgather_local_array('epp', local) 425 ~/pyGSTi/pygsti/objectivefns/objectivefns.py in hessian(self, paramvec) 4877 if self.ex != 0: raise NotImplementedError("Hessian is not implemented for penalty terms yet!") 4878 if paramvec is not None: self.model.from_vector(paramvec) -> 4879 return self._gather_hessian(self._construct_hessian(self.counts, self.total_counts, self.prob_clip_interval)) 4880 4881 def _hessian_from_block(self, hprobs, dprobs12, probs, counts, total_counts, freqs, resource_alloc): ~/pyGSTi/pygsti/objectivefns/objectivefns.py in _construct_hessian(self, counts, total_counts, prob_clip_interval) 1572 k, kmax = 0, len(slicetup_list) 1573 blk_rank = param2_resource_alloc.comm_rank -> 1574 for (slice1, slice2, hprobs, dprobs12) in self.model.sim._iter_atom_hprobs_by_rectangle( 1575 atom, slicetup_list, True, param2_resource_alloc): 1576 local_slice1 = _slct.shift(slice1, -global_param_slice.start) # indices into atom_hessian ~/pyGSTi/pygsti/forwardsims/distforwardsim.py in _iter_atom_hprobs_by_rectangle(self, atom, wrt_slices_list, return_dprobs_12, resource_alloc) 321 # Note: no need to index w/ [atom.element_slice,...] (compare with _iter_hprobs_by_rectangles) 322 # since these arrays are already sized to this particular atom (not to all the host's atoms) --> 323 self._bulk_fill_hprobs_dprobs_atom(hprobs, dprobs1, dprobs2, atom, 324 wrtSlice1, wrtSlice2, resource_alloc) 325 #Note: we give resource_alloc as our local `resource_alloc` above because all the arrays ~/pyGSTi/pygsti/forwardsims/distforwardsim.py in _bulk_fill_hprobs_dprobs_atom(self, array_to_fill, deriv1_array_to_fill, deriv2_array_to_fill, atom, param_slice1, param_slice2, resource_alloc) 257 param_slice2, resource_alloc) 258 --> 259 self._bulk_fill_hprobs_atom(array_to_fill, host_param_slice1, host_param_slice2, atom, 260 param_slice1, param_slice2, resource_alloc) 261 ~/pyGSTi/pygsti/forwardsims/matrixforwardsim.py in _bulk_fill_hprobs_atom(self, array_to_fill, dest_param_slice1, dest_param_slice2, layout_atom, param_slice1, param_slice2, resource_alloc) 1435 self._compute_dproduct_cache(layout_atom.tree, prodCache, scaleCache, 1436 resource_alloc, param_slice2) # computed on rank=0 only -> 1437 hProdCache = self._compute_hproduct_cache(layout_atom.tree, prodCache, dProdCache1, 1438 dProdCache2, scaleCache, resource_alloc, 1439 param_slice1, param_slice2) # computed on rank=0 only ~/pyGSTi/pygsti/forwardsims/matrixforwardsim.py in _compute_hproduct_cache(self, layout_atom_tree, prod_cache, d_prod_cache1, d_prod_cache2, scale_cache, resource_alloc, wrt_slice1, wrt_slice2) 991 dLdR_sym = dLdRa + _np.swapaxes(dLdRb, 0, 1) 992 --> 993 hProdCache[iDest] = _np.dot(hL, R) + dLdR_sym + _np.transpose(_np.dot(L, hR), (1, 2, 0, 3)) 994 995 scale = scale_cache[iDest] - (scale_cache[iLeft] + scale_cache[iRight]) <__array_function__ internals> in dot(*args, **kwargs) KeyboardInterrupt:
Note above cell was executed for demonstration purposes, and was keyboard-interrupted intentionally since it would have taken forever on a single processor.
#write results back to disk
results.write()