Skip to content

Commit 2f61126

Browse files
authored
Merge pull request #448 from sandialabs/feature-faster-layout-creation
Faster Layout and MDC Store Creation
2 parents 74648a6 + 76ff5bf commit 2f61126

File tree

18 files changed

+699
-346
lines changed

18 files changed

+699
-346
lines changed

pygsti/algorithms/core.py

Lines changed: 30 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,8 @@
3333
from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
3434
from pygsti.optimize.customlm import CustomLMOptimizer as _CustomLMOptimizer
3535
from pygsti.optimize.customlm import Optimizer as _Optimizer
36+
from pygsti import forwardsims as _fwdsims
37+
from pygsti import layouts as _layouts
3638

3739
_dummy_profiler = _DummyProfiler()
3840

@@ -400,7 +402,7 @@ def _construct_ab(prep_fiducials, effect_fiducials, model, dataset, op_label_ali
400402
for j, rhostr in enumerate(prep_fiducials):
401403
opLabelString = rhostr + estr # LEXICOGRAPHICAL VS MATRIX ORDER
402404
dsStr = opLabelString.replace_layers_with_aliases(op_label_aliases)
403-
expd_circuit_outcomes = opLabelString.expand_instruments_and_separate_povm(model)
405+
expd_circuit_outcomes = model.expand_instruments_and_separate_povm(opLabelString)
404406
assert(len(expd_circuit_outcomes) == 1), "No instruments are allowed in LGST fiducials!"
405407
unique_key = next(iter(expd_circuit_outcomes.keys()))
406408
outcomes = expd_circuit_outcomes[unique_key]
@@ -429,7 +431,7 @@ def _construct_x_matrix(prep_fiducials, effect_fiducials, model, op_label_tuple,
429431
for j, rhostr in enumerate(prep_fiducials):
430432
opLabelString = rhostr + _circuits.Circuit(op_label_tuple, line_labels=rhostr.line_labels) + estr
431433
dsStr = opLabelString.replace_layers_with_aliases(op_label_aliases)
432-
expd_circuit_outcomes = opLabelString.expand_instruments_and_separate_povm(model)
434+
expd_circuit_outcomes = model.expand_instruments_and_separate_povm(opLabelString)
433435
dsRow_fractions = dataset[dsStr].fractions
434436
assert(len(expd_circuit_outcomes) == nVariants)
435437

@@ -677,16 +679,10 @@ def run_gst_fit(mdc_store, optimizer, objective_function_builder, verbosity=0):
677679
if _np.linalg.norm(mdc_store.model.to_vector() - v_cmp) > 1e-6:
678680
raise ValueError("MPI ERROR: *different* MC2GST start models"
679681
" given to different processors!") # pragma: no cover
680-
681-
#MEM from ..baseobjs.profiler import Profiler
682-
#MEM debug_prof = Profiler(comm)
683-
#MEM debug_prof.print_memory("run_gst_fit1", True)
684-
682+
685683
if objective_function_builder is not None:
686684
objective_function_builder = _objfns.ObjectiveFunctionBuilder.cast(objective_function_builder)
687-
#MEM debug_prof.print_memory("run_gst_fit2", True)
688685
objective = objective_function_builder.build_from_store(mdc_store, printer) # (objective is *also* a store)
689-
#MEM debug_prof.print_memory("run_gst_fit3", True)
690686
else:
691687
assert(isinstance(mdc_store, _objfns.ObjectiveFunction)), \
692688
"When `objective_function_builder` is None, `mdc_store` must be an objective fn!"
@@ -705,14 +701,8 @@ def run_gst_fit(mdc_store, optimizer, objective_function_builder, verbosity=0):
705701

706702
printer.log("Completed in %.1fs" % (_time.time() - tStart), 1)
707703

708-
#if target_model is not None:
709-
# target_vec = target_model.to_vector()
710-
# targetErrVec = _objective_func(target_vec)
711-
# return minErrVec, soln_gs, targetErrVec
712704
profiler.add_time("do_mc2gst: total time", tStart)
713-
#TODO: evTree.permute_computation_to_original(minErrVec) #Doesn't work b/c minErrVec is flattened
714-
# but maybe best to just remove minErrVec from return value since this isn't very useful
715-
# anyway?
705+
716706
return opt_result, objective
717707

718708

@@ -888,10 +878,30 @@ def _max_array_types(artypes_list): # get the maximum number of each array type
888878
#The ModelDatasetCircuitsStore
889879
printer.log('Precomputing CircuitOutcomeProbabilityArray layouts for each iteration.', 2)
890880
precomp_layouts = []
881+
882+
#pre-compute a dictionary caching completed circuits for layout construction performance.
883+
unique_circuits = list({ckt for circuit_list in circuit_lists for ckt in circuit_list})
884+
if isinstance(mdl.sim, (_fwdsims.MatrixForwardSimulator, _fwdsims.MapForwardSimulator)):
885+
precomp_layout_circuit_cache = mdl.sim.create_copa_layout_circuit_cache(unique_circuits, mdl, dataset=dataset)
886+
else:
887+
precomp_layout_circuit_cache = None
888+
891889
for i, circuit_list in enumerate(circuit_lists):
892890
printer.log(f'Layout for iteration {i}', 2)
893-
precomp_layouts.append(mdl.sim.create_layout(circuit_list, dataset, resource_alloc, array_types, verbosity= printer - 1))
894-
891+
precomp_layouts.append(mdl.sim.create_layout(circuit_list, dataset, resource_alloc, array_types, verbosity= printer - 1,
892+
layout_creation_circuit_cache = precomp_layout_circuit_cache))
893+
894+
#precompute a cache of possible outcome counts for each circuits to accelerate MDC store creation
895+
if isinstance(mdl, _models.model.OpModel):
896+
if precomp_layout_circuit_cache is not None: #then grab the split circuits from there.
897+
expanded_circuit_outcome_list = mdl.bulk_expand_instruments_and_separate_povm(unique_circuits,
898+
completed_circuits= precomp_layout_circuit_cache['completed_circuits'].values())
899+
else:
900+
expanded_circuit_outcome_list = mdl.bulk_expand_instruments_and_separate_povm(unique_circuits)
901+
outcome_count_by_circuit_cache = {ckt: len(outcome_tup) for ckt,outcome_tup in zip(unique_circuits, expanded_circuit_outcome_list)}
902+
else:
903+
outcome_count_by_circuit_cache = {ckt: mdl.compute_num_outcomes(ckt) for ckt in unique_circuits}
904+
895905
with printer.progress_logging(1):
896906
for i in range(starting_index, len(circuit_lists)):
897907
circuitsToEstimate = circuit_lists[i]
@@ -908,7 +918,8 @@ def _max_array_types(artypes_list): # get the maximum number of each array type
908918
mdl.basis = start_model.basis # set basis in case of CPTP constraints (needed?)
909919
initial_mdc_store = _objfns.ModelDatasetCircuitsStore(mdl, dataset, circuitsToEstimate, resource_alloc,
910920
array_types=array_types, verbosity=printer - 1,
911-
precomp_layout = precomp_layouts[i])
921+
precomp_layout = precomp_layouts[i],
922+
outcome_count_by_circuit=outcome_count_by_circuit_cache)
912923
mdc_store = initial_mdc_store
913924

914925
for j, obj_fn_builder in enumerate(iteration_objfn_builders):

pygsti/circuits/__init__.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
1111
#***************************************************************************************************
1212

13-
from .circuit import Circuit
13+
from .circuit import Circuit, SeparatePOVMCircuit
1414
from .circuitlist import CircuitList
1515
from .circuitstructure import CircuitPlaquette, FiducialPairPlaquette, \
1616
GermFiducialPairPlaquette, PlaquetteGridCircuitStructure

pygsti/circuits/circuit.py

Lines changed: 32 additions & 115 deletions
Original file line numberDiff line numberDiff line change
@@ -1095,13 +1095,12 @@ def _proc_lines_arg(self, lines):
10951095
def _proc_key_arg(self, key):
10961096
""" Pre-process the key argument used by many methods """
10971097
if isinstance(key, tuple):
1098-
if len(key) != 2: return IndexError("Index must be of the form <layerIndex>,<lineIndex>")
1099-
layers = key[0]
1100-
lines = key[1]
1098+
if len(key) != 2:
1099+
return IndexError("Index must be of the form <layerIndex>,<lineIndex>")
1100+
else:
1101+
return key[0], key[1]
11011102
else:
1102-
layers = key
1103-
lines = None
1104-
return layers, lines
1103+
return key, None
11051104

11061105
def _layer_components(self, ilayer):
11071106
""" Get the components of the `ilayer`-th layer as a list/tuple. """
@@ -1191,22 +1190,41 @@ def extract_labels(self, layers=None, lines=None, strict=True):
11911190
`layers` is a single integer and as a `Circuit` otherwise.
11921191
Note: if you want a `Circuit` when only selecting one layer,
11931192
set `layers` to a slice or tuple containing just a single index.
1193+
Note that the returned circuit doesn't retain any original
1194+
metadata, such as the compilable layer indices or occurence id.
11941195
"""
11951196
nonint_layers = not isinstance(layers, int)
11961197

11971198
#Shortcut for common case when lines == None and when we're only taking a layer slice/index
1198-
if lines is None:
1199-
assert(layers is not None)
1200-
if nonint_layers is False: return self.layertup[layers]
1201-
if isinstance(layers, slice) and strict is True: # if strict=False, then need to recompute line labels
1202-
return Circuit._fastinit(self._labels[layers], self._line_labels, not self._static)
1199+
if lines is None and layers is not None:
1200+
if self._static:
1201+
if not nonint_layers:
1202+
return self._labels[layers]
1203+
if isinstance(layers, slice) and strict is True: # if strict=False, then need to recompute line labels
1204+
#can speed this up a measurably by manually computing the new hashable tuple value and hash
1205+
if not self._line_labels in (('*',), ()):
1206+
new_hashable_tup = self._labels[layers] + ('@',) + self._line_labels
1207+
else:
1208+
new_hashable_tup = self._labels[layers]
1209+
ret = Circuit.__new__(Circuit)
1210+
return ret._copy_init(self._labels[layers], self._line_labels, not self._static, hashable_tup= new_hashable_tup, precomp_hash=hash(new_hashable_tup))
1211+
else:
1212+
if not nonint_layers:
1213+
return self.layertup[layers]
1214+
if isinstance(layers, slice) and strict is True: # if strict=False, then need to recompute line labels
1215+
return Circuit._fastinit(self._labels[layers], self._line_labels, not self._static)
1216+
#otherwise assert both are not None:
1217+
12031218

12041219
layers = self._proc_layers_arg(layers)
12051220
lines = self._proc_lines_arg(lines)
12061221
if len(layers) == 0 or len(lines) == 0:
1207-
return Circuit._fastinit(() if self._static else [],
1208-
tuple(lines) if self._static else lines,
1209-
not self._static) if nonint_layers else None # zero-area region
1222+
if self._static:
1223+
return Circuit._fastinit((), tuple(lines), False) # zero-area region
1224+
else:
1225+
return Circuit._fastinit(() if self._static else [],
1226+
tuple(lines) if self._static else lines,
1227+
not self._static) # zero-area region
12101228

12111229
ret = []
12121230
if self._static:
@@ -4427,107 +4445,6 @@ def done_editing(self):
44274445
self._hashable_tup = self.tup
44284446
self._hash = hash(self._hashable_tup)
44294447

4430-
def expand_instruments_and_separate_povm(self, model, observed_outcomes=None):
4431-
"""
4432-
Creates a dictionary of :class:`SeparatePOVMCircuit` objects from expanding the instruments of this circuit.
4433-
4434-
Each key of the returned dictionary replaces the instruments in this circuit with a selection
4435-
of their members. (The size of the resulting dictionary is the product of the sizes of
4436-
each instrument appearing in this circuit when `observed_outcomes is None`). Keys are stored
4437-
as :class:`SeparatePOVMCircuit` objects so it's easy to keep track of which POVM outcomes (effects)
4438-
correspond to observed data. This function is, for the most part, used internally to process
4439-
a circuit before computing its outcome probabilities.
4440-
4441-
Parameters
4442-
----------
4443-
model : Model
4444-
The model used to provide necessary details regarding the expansion, including:
4445-
4446-
- default SPAM layers
4447-
- definitions of instrument-containing layers
4448-
- expansions of individual instruments and POVMs
4449-
4450-
Returns
4451-
-------
4452-
OrderedDict
4453-
A dict whose keys are :class:`SeparatePOVMCircuit` objects and whose
4454-
values are tuples of the outcome labels corresponding to this circuit,
4455-
one per POVM effect held in the key.
4456-
"""
4457-
complete_circuit = model.complete_circuit(self)
4458-
expanded_circuit_outcomes = _collections.OrderedDict()
4459-
povm_lbl = complete_circuit[-1] # "complete" circuits always end with a POVM label
4460-
circuit_without_povm = complete_circuit[0:len(complete_circuit) - 1]
4461-
4462-
def create_tree(lst):
4463-
subs = _collections.OrderedDict()
4464-
for el in lst:
4465-
if len(el) > 0:
4466-
if el[0] not in subs: subs[el[0]] = []
4467-
subs[el[0]].append(el[1:])
4468-
return _collections.OrderedDict([(k, create_tree(sub_lst)) for k, sub_lst in subs.items()])
4469-
4470-
def add_expanded_circuit_outcomes(circuit, running_outcomes, ootree, start):
4471-
"""
4472-
"""
4473-
cir = circuit if start == 0 else circuit[start:] # for performance, avoid uneeded slicing
4474-
for k, layer_label in enumerate(cir, start=start):
4475-
components = layer_label.components
4476-
#instrument_inds = _np.nonzero([model._is_primitive_instrument_layer_lbl(component)
4477-
# for component in components])[0] # SLOWER than statement below
4478-
instrument_inds = _np.array([i for i, component in enumerate(components)
4479-
if model._is_primitive_instrument_layer_lbl(component)])
4480-
if instrument_inds.size > 0:
4481-
# This layer contains at least one instrument => recurse with instrument(s) replaced with
4482-
# all combinations of their members.
4483-
component_lookup = {i: comp for i, comp in enumerate(components)}
4484-
instrument_members = [model._member_labels_for_instrument(components[i])
4485-
for i in instrument_inds] # also components of outcome labels
4486-
for selected_instrmt_members in _itertools.product(*instrument_members):
4487-
expanded_layer_lbl = component_lookup.copy()
4488-
expanded_layer_lbl.update({i: components[i] + "_" + sel
4489-
for i, sel in zip(instrument_inds, selected_instrmt_members)})
4490-
expanded_layer_lbl = _Label([expanded_layer_lbl[i] for i in range(len(components))])
4491-
4492-
if ootree is not None:
4493-
new_ootree = ootree
4494-
for sel in selected_instrmt_members:
4495-
new_ootree = new_ootree.get(sel, {})
4496-
if len(new_ootree) == 0: continue # no observed outcomes along this outcome-tree path
4497-
else:
4498-
new_ootree = None
4499-
4500-
add_expanded_circuit_outcomes(circuit[0:k] + Circuit((expanded_layer_lbl,)) + circuit[k + 1:],
4501-
running_outcomes + selected_instrmt_members, new_ootree, k + 1)
4502-
break
4503-
4504-
else: # no more instruments to process: `cir` contains no instruments => add an expanded circuit
4505-
assert(circuit not in expanded_circuit_outcomes) # shouldn't be possible to generate duplicates...
4506-
elabels = model._effect_labels_for_povm(povm_lbl) if (observed_outcomes is None) \
4507-
else tuple(ootree.keys())
4508-
outcomes = tuple((running_outcomes + (elabel,) for elabel in elabels))
4509-
expanded_circuit_outcomes[SeparatePOVMCircuit(circuit, povm_lbl, elabels)] = outcomes
4510-
4511-
ootree = create_tree(observed_outcomes) if observed_outcomes is not None else None # tree of observed outcomes
4512-
# e.g. [('0','00'), ('0','01'), ('1','10')] ==> {'0': {'00': {}, '01': {}}, '1': {'10': {}}}
4513-
4514-
if model._has_instruments():
4515-
add_expanded_circuit_outcomes(circuit_without_povm, (), ootree, start=0)
4516-
else:
4517-
# It may be helpful to cache the set of elabels for a POVM (maybe within the model?) because
4518-
# currently the call to _effect_labels_for_povm may be a bottleneck. It's needed, even when we have
4519-
# observed outcomes, because there may be some observed outcomes that aren't modeled (e.g. leakage states)
4520-
if observed_outcomes is None:
4521-
elabels = model._effect_labels_for_povm(povm_lbl)
4522-
else:
4523-
possible_lbls = set(model._effect_labels_for_povm(povm_lbl))
4524-
elabels = tuple([oo for oo in ootree.keys() if oo in possible_lbls])
4525-
outcomes = tuple(((elabel,) for elabel in elabels))
4526-
expanded_circuit_outcomes[SeparatePOVMCircuit(circuit_without_povm, povm_lbl, elabels)] = outcomes
4527-
4528-
return expanded_circuit_outcomes
4529-
4530-
45314448
class CompressedCircuit(object):
45324449
"""
45334450
A "compressed" Circuit that requires less disk space.

0 commit comments

Comments
 (0)