Time-dependent models and gate set tomography

This tutorial demonstrates how time dependence can be added to models in pyGSTi and, since gate set tomography (GST) just optimizes model parameters, how to run time-dependent GST.

**Notice: this topic describes "beta level" functionality in pyGSTi!** It may contain bugs and holes in its implementation, which will be addressed in future releases.

In [ ]:
import pygsti
from pygsti.modelpacks import smq1Q_XYI
import numpy as np

Time dependent models

To make a time-dependent Model, you create a time dependent gate or operation and add this to any of the models in pyGSTi. (Expert note: this isn't quite true - currently, only models with sim_type="map" support time-dependent evaluation of circuit outcomes, so we're currently limited to using this simulation type.) Here's an example of how to make a custom idle operation that depolarizes its input state more and more over time:

In [ ]:
class MyTimeDependentIdle(pygsti.obj.DenseOperator):
    """And idle that depolarizes over time with a parameterized rate"""
    def __init__(self, initial_depol_rate):
        #initialize with no noise
        super(MyTimeDependentIdle,self).__init__(np.identity(4,'d'), "densitymx") # this is *super*-operator, so "densitymx"
    def num_params(self): 
        return 1 # we have two parameters
    def to_vector(self):
        return np.array([self.depol_rate],'d') #our parameter vector
    def from_vector(self, v, close=False, nodirty=False):
        #initialize from parameter vector v
        self.depol_rate = v[0]
        if not nodirty: self.dirty = True # mark that paramvec (self.to_vector()) may have changed
    def set_time(self,t):
        a = 1.0-min(self.depol_rate*t,1.0)
        # .base is a member of DenseOperator and is a numpy array that is 
        # the dense Pauli transfer matrix of this operator
        # Technical note: use [:,:] here b/c we don't want to change id of self.base
        self.base[:,:] = np.array([[1,   0,   0,   0],
                                  [0,   a,   0,   0],
                                  [0,   0,   a,   0],
                                  [0,   0,   0,   a]],'d')
    def transform(self, S):
        # Update self with inverse(S) * self * S (used in gauge optimization)
        raise NotImplementedError("MyTimeDependentIdle cannot be transformed!")

The key piece to note in the above class is the set_time method, which will be called sometime after from_vector and takes over responsiblility (from from_vector) for setting the object's .base member to the process matrix based on the parameters (in from_vector's v and the time given to set_time).

Here's an example of how to see what a MyTimeDependentIdle(1.0) gate looks like at the time 0.1:

In [ ]:
t = 0.1
Gi_at_t = MyTimeDependentIdle(1.0)
Gi_matrix_at_t = Gi_at_t.todense()

We can add a MyTimeDependentIdle gate to a model just like any other operator (in pyGSTi all operators are considered potentially time-dependent, and so the base class of our idle gate is DenseOperator just as it would be if we were creating a custom time-independent gate):

In [ ]:
mdl = smq1Q_XYI.target_model(sim_type="map")
mdl['Gi'] = MyTimeDependentIdle(1.0)

There you have it - mdl is a time-dependent model, where Gi depolarizes with strength equal to the current time. To compute the probability of a circuit, GiGi for example, we just call the usual probs function but specify a time argument:

In [ ]:
mdl.probs( ('Gi','Gi'), time=0.1)

The zero probability is equal to 0.5 * (1 + 0.9**2) = 0.905, where the 0.9 comes from the Gi gate depolarization rate of 0.1 at time 0.1. Note that this is the same as what you'd get using the Gi_matrix_at_t above (since our "t" was 0.1):

In [ ]:
E = mdl['Mdefault']['0']
rho = mdl['rho0']
print(np.dot(E.T, np.dot(Gi_matrix_at_t, np.dot(Gi_matrix_at_t, rho))))

Time-dependent (or "time aware") circuits

Circuit objects may include time information: labels within a circuit (e.g. "Gi") may contain a relative time giving the duration of the operation being labeled. By default, all labels have zero duration, meaning all the operations within the circuit are interpreted as occurring at the same time. The below example gives the Gi gate a duration of 0.1, so that in the circuit simulation the first Gi occurs at time 0.1 and the second at 0.2:

In [ ]:
Gi_with_duration = pygsti.obj.Label('Gi',time=0.1)
mdl.probs( (Gi_with_duration,Gi_with_duration), time=0.1)

This is the same as the following:

In [ ]:
Gi_matrix_at_t1 = Gi_at_t.todense()
Gi_matrix_at_t2 = Gi_at_t.todense()
print(np.dot(E.T, np.dot(Gi_matrix_at_t2, np.dot(Gi_matrix_at_t1, rho))))

You can also use the following "!"-shorthand (exclamation point followed by time) notation to specify label durations:

In [ ]:
mdl.probs( (('Gi','!0.1'),('Gi','!0.1')), time=0.1)

Time dependent data

When DataSet objects contain timestamped data, these timestamps indicate at what absolute time the relevant circuit began executing when it produced certain data. These time values correspond to those given to the time argument of probs above.

At first, we don't bother with "time-aware" circuits, and just create a list of two sample circuits. We then use the times argument of generate_fake_data to construct a DataSet with 100 samples of data taken at each of three times: 0, 0.1, and 0.2 (arbitrary time units). By setting sample_error="none" we can see the underlying outcome probabilities in the data (and how the depolarization caused by Gi increases with time):

In [ ]:
circuits = pygsti.construction.circuit_list([ ('Gi',), ('Gi','Gi')]) # just pick some circuits

ds = pygsti.construction.generate_fake_data(mdl, circuits, n_samples=100,
                                            sample_error='none', seed=1234, times=[0,0.1,0.2])

A DataSet with timestamps displays 3 parallel arrays for each circuit: "Outcome Label Indices", "Time stamps", and "Repetitions". Each index corresponds to a bin of some number (given by "Repetitions") of X-outcomes (X given by "Outcome Label Indices") occuring at some time (given by "Time stamps"). We see that for each of the two circuits there are bins of 0- and 1-outcomes at each of times 0, 0.1, and 0.2. Summing the bin counts (outcome repetitions) at each time, for a given circuit, gives 100.

We can also add a duration of 0.05 time units to each "Gi" gate. This makes the depolarization of the length-2 sequence a bit worse because the second application of "Gi" occurs at a time 0.05 units after the start of the circuit, at which point the noise on the gate as increased:

In [ ]:
circuits = pygsti.construction.circuit_list([ (('Gi','!0.05'),), (('Gi','!0.05'),('Gi','!0.05'))])

ds = pygsti.construction.generate_fake_data(mdl, circuits, n_samples=100,
                                            sample_error='none', seed=1234, times=[0,0.1,0.2])

Time-dependent gate set tomography (TD-GST)

To run gate set tomography, we'll need more sequences than the two in the example above. We'll generate some timestamped data for the standard set of GST sequences for a 1-qubit $X(\pi/2)$, $Y(\pi/2)$, $I$ gate set. In particular, we create a data-generating model that has a MyTimeDependentIdle "Gi" gate with a depolarization "acceleration" rate of 1.0, and we generate 10 counts at each of 10 equally spaced times between 0 and 0.3.

In [ ]:
prep_fiducials, meas_fiducials = smq1Q_XYI.prep_fiducials(), smq1Q_XYI.meas_fiducials()
germs = smq1Q_XYI.germs()
maxLengths = [1,2]

mdl_datagen = smq1Q_XYI.target_model(sim_type="map").depolarize(op_noise=0.01, spam_noise=0.001)
mdl_datagen['Gi'] = MyTimeDependentIdle(1.0)

edesign = pygsti.protocols.StandardGSTDesign(smq1Q_XYI.target_model(), prep_fiducials,
                                             meas_fiducials, germs, maxLengths)
#listOfExperiments = pygsti.construction.make_lsgst_experiment_list(
#    mdl_datagen, prep_fiducials, meas_fiducials, germs, maxLengths)

#Data for initial non-sparse mode
ds = pygsti.construction.generate_fake_data(mdl_datagen, edesign.all_circuits_needing_data, n_samples=10,
                                            sample_error="binomial", seed=1234, times=np.linspace(0,0.3,10))

We can run GST on this timestamped data similar to any other data, using do_long_sequence_gst. The key difference is that the advanced option "timeDependent" is set to True. This tells do_long_sequence_gst to take the timestamps in the given DataSet seriously, and perform time-dependent circuit simulations rather than aggregating the counts across all times (the behavior when the default, "timeDependent": False, is used).

Running time-dependent GST with 10 timesteps requires 10 times the number of circuit simulations (each circuit needs to be simulated 10 times). This, coupled with the fact that this the time-dependent simulation routines are less optimized in pyGSTi, means this running time-dependent GST is significantly slower than normal GST. Note also that we set gauge_opt_params=False. This disables gauge optimization, and this is necessary since it won't work because our MyTimeDependentIdle operation doesn't implement transform (the action of a gauge transformation).

The cell below will take around 5 minutes to run.

In [ ]:
target_model = smq1Q_XYI.target_model("TP", sim_type="map") # TP-constraints on the non-Gi gates
target_model['Gi'] = MyTimeDependentIdle(0.0)
target_model.set_simtype('map', max_cache_size=0)

builders = pygsti.protocols.GSTObjFnBuilders([pygsti.objects.TimeDependentPoissonPicLogLFunction.builder()],[])
gst = pygsti.protocols.GateSetTomography(target_model, gaugeopt_suite=None,
                                         objfn_builders=builders, optimizer={'tol': 1e-4} )
data = pygsti.protocols.ProtocolData(edesign, ds)
results = gst.run(data)

# results = pygsti.do_long_sequence_gst(ds, target_model, prep_fiducials, meas_fiducials,
#                                       germs, maxLengths, gauge_opt_params=False, verbosity=3,
#                                       advanced_options={'timeDependent': True, # enable time-dependent circuit simulation
#                                                        'starting_point': 'target', # so we don't try to do LGST
#                                                        'tolerance': 1e-4}) # to speed up

We can extract the (non-gauge-optimizeed) best-fit model from results, and see what depolarization "acceleration" was found. We find that the value is reasonably close to the value of 1.0 that we used to generate the data.

In [ ]:
final_mdl = results.estimates['GateSetTomography'].models['final iteration estimate']
print("Time-dependent idle parameters = ",final_mdl['Gi'].to_vector())