Skip to content

WIP: crosstalk-free GST#641

Draft
nkoskelo wants to merge 149 commits intodevelopfrom
nick_xfgst
Draft

WIP: crosstalk-free GST#641
nkoskelo wants to merge 149 commits intodevelopfrom
nick_xfgst

Conversation

@nkoskelo
Copy link
Contributor

@nkoskelo nkoskelo commented Aug 19, 2025

TODO: Add description.

dhothem and others added 30 commits January 20, 2025 15:49
…c parent class of CrosstalkFreeExperimentDesign.
… just raise an more informative error). Add a copy in PermutationOperator.pp_braiding_operators.
…ys compute model violation w.r.t. logl objective function, rather than whatever might have been used as the final objective function.
…s. These are hacky changes. A proper change would be to use try-catch if encountering linalg errors.
…functions only use idles with explicit target qubits, rather than implicit global idles
… explicitmodel.py. Update convert_members_inplace so it can handle failures of conversion when dealing with embedded ops.
if isinstance(x, _Label): return x
# # do this manually when desired, as it "boxes" a circuit being inserted
#elif isinstance(x,Circuit): return x.to_circuit_label()
elif isinstance(x,Circuit): return x.to_label()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Neither Nick nor I remember making this change. The commented-out branch was broken. The uncommenting seems valid based on the to_label member function of Circuit objects.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tensor circuit construction should retain knowledge of original lane structure.

#self._times = None # for FUTURE expansion
self.auxinfo = {} # for FUTURE expansion / user metadata
self.saved_auxinfo = {}
self.saved_auxinfo["lanes"] = {tuple(line_labels): self._labels}
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We are saving off the lane structure here so that it can be used in the future if the error models are cross-talk free.

rileyjmurray added a commit that referenced this pull request Jan 28, 2026
The changes in this PR were originally made for crosstalk-free GST work
(see #641). We've factored out those changes into this PR to reduce the
size of #641 (and get these incidental changes in sooner rather than
later).

Changes are overwhelmingly type annotations.

---------

Co-authored-by: Nick Koskelo <koskelo2@illinois.edu>
Co-authored-by: nkoskelo <129830924+nkoskelo@users.noreply.github.com>
Copy link
Contributor

@rileyjmurray rileyjmurray left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nkoskelo now that we've got the xfgst-incidentals merged in we'll want to create a new (smaller scope) crosstalk-free GST PR. The new PR should include everything that's useful in XFGST modeling except the things related to a new forward simulator. I've left comments on the code that matches that description. Please do a pass of your own and let me know if you agree with what I've highlighted and/or if you think I missed something.

Once we agree on the plan here we can close this PR and create a new branch starting from develop!

Comment on lines +373 to +417
def create_explicit(self, composedops_as_views=False, spam_type='full TP', op_type='CPTPLND'):
from pygsti.models import ExplicitOpModel
if isinstance(composedops_as_views, str):
# Old positional argument behavior.
spam_type = composedops_as_views
op_type = composedops_as_views
composedops_as_views = False

state = self.__getstate__()
state['povms'] = state['povm_blks']['layers']
state['preps'] = state['prep_blks']['layers']
opdict = _OrderedMemberDict(None, spam_type, None, dict())
opdict.update(state['_opcaches']['complete-layers'])
state['operations'] = opdict

for v in state['operations'].values():
v.unlink_parent(force=True)
for v in state['povms'].values():
v.unlink_parent(force=True)
for v in state['preps'].values():
v.unlink_parent(force=True)

eom = ExplicitOpModel(self.state_space)
eom.preps.update(state['preps'])
eom.povms.update(state['povms'])
eom.operations.update(state['operations'])

eom.convert_members_inplace(
to_type=spam_type, categories_to_convert=('preps', 'povms'),
allow_smaller_pp_basis=True, flatten_structure=True
)
eom.convert_members_inplace(
to_type=op_type, categories_to_convert=('operations',),
allow_smaller_pp_basis=True
)
if composedops_as_views:
from pygsti.modelmembers.operations import ComposedOp
allop_keys = eom.operations.keys()
for k in allop_keys:
curr_op = eom[k]
if isinstance(curr_op, ComposedOp) and all([subk in allop_keys for subk in k]):
view_op = ComposedOp([eom[subk] for subk in k])
eom[k] = view_op
eom._rebuild_paramvec()
return eom
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Worth merging.

Comment on lines +361 to +410
def convert_members_inplace(self, to_type,
categories_to_convert='all', labels_to_convert='all',
ideal_model=None, flatten_structure=False, set_default_gauge_group=False,
cptp_truncation_tol= 1e-7, spam_cp_penalty = 1e-7, allow_smaller_pp_basis=False
):
"""
TODO: docstring -- like set_all_parameterizations but doesn't set default gauge group by default

allow_smaller_pp_basis : bool

Setting allow_smaller_pp_basis=True allows dimension mismatches between
this ExplicitOpModel's operations and the dimensions we'd expect for
operations based on the properties of self.basis.

We can ONLY handle mismatches when self.basis.name indicates a Pauli product basis
or a tensor product thereof. We handle a failed conversion by trying a second time,
passing in the string literal `pp` instead of `self.basis` to _op.convert(...).

"""
if isinstance(categories_to_convert, str): categories_to_convert = (categories_to_convert,)
assert all(c in ['all', 'ops', 'operations', 'instruments', 'preps', 'povms'] for c in categories_to_convert)
fallback_basis = '' if not allow_smaller_pp_basis else self.basis.name.replace('pp','').replace('*','') + 'pp'
ideal_model = ModelView.cast(ideal_model)
roster = Roster(labels_to_convert)
if any([c in categories_to_convert for c in ('all', 'ops', 'operations')]):
for lbl, gate in self.operations.items():
if labels_to_convert == 'all' or lbl in labels_to_convert:
ideal = ideal_model.operations.get(lbl, None) if (ideal_model is not None) else None
self.operations[lbl] = _op.convert(gate, to_type, self.basis, ideal, flatten_structure, cptp_truncation_tol)
op_items = [(k,v) for (k,v) in self.operations.items() if k in roster]
for lbl, gate in op_items:
ideal = ideal_model.operations[lbl]
try:
op = _op.convert(gate, to_type, self.basis, ideal, flatten_structure, cptp_truncation_tol)
except ValueError as e:
if not fallback_basis == 'pp':
raise e
op = _op.convert(gate, to_type, 'pp', ideal, flatten_structure, cptp_truncation_tol)
self.operations[lbl] = op
if any([c in categories_to_convert for c in ('all', 'instruments')]):
for lbl, inst in self.instruments.items():
if labels_to_convert == 'all' or lbl in labels_to_convert:
ideal = ideal_model.instruments.get(lbl, None) if (ideal_model is not None) else None
self.instruments[lbl] = _instrument.convert(inst, to_type, self.basis, ideal, flatten_structure)
inst_items = [(k,v) for (k,v) in self.instruments.items() if k in roster]
for lbl, inst in inst_items:
ideal = ideal_model.instruments[lbl]
self.instruments[lbl] = _instrument.convert(inst, to_type, self.basis, ideal, flatten_structure)
if any([c in categories_to_convert for c in ('all', 'preps')]):
for lbl, prep in self.preps.items():
if labels_to_convert == 'all' or lbl in labels_to_convert:
ideal = ideal_model.preps.get(lbl, None) if (ideal_model is not None) else None
self.preps[lbl] = _state.convert(prep, to_type, self.basis, ideal, flatten_structure, cp_penalty=spam_cp_penalty)
prep_items = [(k,v) for (k,v) in self.preps.items() if k in roster]
for lbl, prep in prep_items:
ideal = ideal_model.preps[lbl]
self.preps[lbl] = _state.convert(prep, to_type, self.basis, ideal, flatten_structure, cp_penalty=spam_cp_penalty)
if any([c in categories_to_convert for c in ('all', 'povms')]):
for lbl, povm in self.povms.items():
if labels_to_convert == 'all' or lbl in labels_to_convert:
ideal = ideal_model.povms.get(lbl, None) if (ideal_model is not None) else None
self.povms[lbl] = _povm.convert (povm, to_type, self.basis, ideal, flatten_structure, cp_penalty=spam_cp_penalty)
povm_items = [(k,v) for (k,v) in self.povms.items() if k in roster]
for lbl, povm in povm_items:
ideal = ideal_model.povms[lbl]
self.povms[lbl] = _povm.convert(povm, to_type, self.basis, ideal, flatten_structure, cp_penalty=spam_cp_penalty)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe worth upstreaming. Would need the Roster and ModelView classes defined at the top of the file if so.

Comment on lines +204 to +211

def tensor_circuits(self, other_circuitlist, new_name=None):
assert len(self) == len(other_circuitlist)
circuits = []
for c1,c2 in zip(self._circuits, other_circuitlist._circuits):
circuits.append(c1.tensor_circuit(c2))
out = CircuitList(circuits, name=new_name)
return out
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe worth merging.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TODO: determine if we want to keep these new functions.

Can see these functions in action in test/unit/objects/test_circuit_splitting.py.

Comment on lines 710 to +910
@@ -739,6 +751,9 @@ def _compute_product_cache(self, layout_atom_tree, resource_alloc):
eval_tree = layout_atom_tree
cacheSize = len(eval_tree)
prodCache = _np.zeros((cacheSize, dim, dim), 'd')
# ^ This assumes assignments prodCache[i] = <2d numpy array>.
# It would be better for this to be a dict (mapping _most likely_
# to ndarrays) if we don't need slicing or other axis indexing.
scaleCache = _np.zeros(cacheSize, 'd')

for iDest, iRight, iLeft in eval_tree:
@@ -752,7 +767,10 @@ def _compute_product_cache(self, layout_atom_tree, resource_alloc):
else:
gate = self.model.circuit_layer_operator(opLabel, 'op').to_dense("minimal")
nG = max(_nla.norm(gate), 1.0)
# ^ This indicates a need to compute norms of the operation matrices. Can't do this
# with scipy.linalg if gate is represented implicitly.
prodCache[iDest] = gate / nG
# ^ Indicates a need to overload division by scalars.
scaleCache[iDest] = _np.log(nG)
continue

@@ -765,10 +783,12 @@ def _compute_product_cache(self, layout_atom_tree, resource_alloc):
scaleCache[iDest] = scaleCache[iLeft] + scaleCache[iRight]

if prodCache[iDest].max() < _PSMALL and prodCache[iDest].min() > -_PSMALL:
nL, nR = max(_nla.norm(L), _np.exp(-scaleCache[iLeft]),
1e-300), max(_nla.norm(R), _np.exp(-scaleCache[iRight]), 1e-300)
nL = max(_nla.norm(L), _np.exp(-scaleCache[iLeft]), 1e-300)
nR = max(_nla.norm(R), _np.exp(-scaleCache[iRight]), 1e-300)
sL, sR = L / nL, R / nR
prodCache[iDest] = _np.dot(sL, sR); scaleCache[iDest] += _np.log(nL) + _np.log(nR)
# ^ Again, shows the need to overload division by scalars.
prodCache[iDest] = sL @ sR
scaleCache[iDest] += _np.log(nL) + _np.log(nR)

nanOrInfCacheIndices = (~_np.isfinite(prodCache)).nonzero()[0] # may be duplicates (a list, not a set)
# since all scaled gates start with norm <= 1, products should all have norm <= 1
@@ -839,6 +859,8 @@ def _compute_dproduct_cache(self, layout_atom_tree, prod_cache, scale_cache,

tSerialStart = _time.time()
dProdCache = _np.zeros((cacheSize,) + deriv_shape)
# ^ I think that deriv_shape will be a tuple of length > 2.
# (Based on how swapaxes is used in the loop below ...)
wrtIndices = _slct.indices(wrt_slice) if (wrt_slice is not None) else None

for iDest, iRight, iLeft in eval_tree:
@@ -852,6 +874,9 @@ def _compute_dproduct_cache(self, layout_atom_tree, prod_cache, scale_cache,
#doperation = self.dproduct( (opLabel,) , wrt_filter=wrtIndices)
doperation = self._doperation(opLabel, wrt_filter=wrtIndices)
dProdCache[iDest] = doperation / _np.exp(scale_cache[iDest])
# ^ Need a way to track tensor product structure in whatever's
# being returned by self._doperation (presumably it's a tensor ...)

continue

tm = _time.time()
@@ -862,18 +887,31 @@ def _compute_dproduct_cache(self, layout_atom_tree, prod_cache, scale_cache,
# since then matrixOf(circuit[i]) = matrixOf(circuit[iLeft]) * matrixOf(circuit[iRight])
L, R = prod_cache[iLeft], prod_cache[iRight]
dL, dR = dProdCache[iLeft], dProdCache[iRight]
dProdCache[iDest] = _np.dot(dL, R) + \
_np.swapaxes(_np.dot(L, dR), 0, 1) # dot(dS, T) + dot(S, dT)
term1 = _np.dot(dL, R)
term2 = _np.swapaxes(_np.dot(L, dR), 0, 1)
# ^ From the numpy docs on .dot :
#
# If a is an N-D array and b is an M-D array (where M>=2),
# it is a sum product over the last axis of a and the second-to-last axis of b:
#
# dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])
#
dProdCache[iDest] = term1 + term2 # dot(dS, T) + dot(S, dT)
# ^ We need addition of tensor-product-structured "doperators."

profiler.add_time("compute_dproduct_cache: dots", tm)
profiler.add_count("compute_dproduct_cache: dots")

scale = scale_cache[iDest] - (scale_cache[iLeft] + scale_cache[iRight])
if abs(scale) > 1e-8: # _np.isclose(scale,0) is SLOW!
dProdCache[iDest] /= _np.exp(scale)
if dProdCache[iDest].max() < _DSMALL and dProdCache[iDest].min() > -_DSMALL:
# ^ Need the tensor-product-structured "doperators" to have .max() and .min()
# methods.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes might be worth keeping.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

todo: figure out if we want to merge these changes

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we'll want to merge this

Comment on lines +1281 to +1292
if len(inds) == 1:
if add:
a[inds] += rhs # all integers or slices behave nicely
a[inds] += rhs
else:
a[inds] = rhs # all integers or slices behave nicely
a[inds] = rhs
elif all([isinstance(i, (int, slice)) for i in inds]):
assert len(inds) == rhs.shape[1]
for ind, rhs_vec in zip(inds, rhs.T):
if add:
a[ind] += rhs_vec # all integers or slices behave nicely
else:
a[ind] = rhs_vec # all integers or slices behave nicely
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we'll want to keep this

Comment on lines +15 to +17
pygsti.optimize.customsolve.CUSTOM_SOLVE_THRESHOLD = 10
wcomm = MPI.COMM_WORLD
print(f'Running with CUSTOM_SOLVE_THRESHOLD = {pygsti.optimize.customsolve.CUSTOM_SOLVE_THRESHOLD}')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be merged

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably want to keep this

@rileyjmurray
Copy link
Contributor

@nkoskelo please see my most recent review. I'd like for you to create a new "Crosstalk-free GST infrastructure" PR based on the files indicated there. As soon as that PR is open we can/should close this one.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants