RB code updated so that Mohan's crosstalk detection code should work …

…again. Also some misc changes to the drift code that'll later be written over when merging in the new feature-newdrift branch

packages/pygsti/extras/drift/core.py

@@ -724,8 +724,7 @@ def estimate_probability_trajectories(results, modelSelector=(('per','per','avg'
 
             meandict = {o : _np.mean(results.timeseries[e][s][o])/results.number_of_counts for o in outcomes[:-1]}
             #outcomeIndlist = list(range(results.number_of_outcomes))
-            nullmodel = _mdl.ProbabilityTrajectoryModel(outcomes, modeltype='null', 
-                                                        hyperparameters={'basisfunctionInds':[0.,]}, 
+            nullmodel = _mdl.NullProbabilityTrajectoryModel(outcomes, hyperparameters=[0.,], 
                                                         parameters=meandict)
 
             results.add_reconstruction(ent, seq, nullmodel, modelSelector='null', estimator='null', auxDict={}, 
@@ -754,18 +753,19 @@ def estimate_probability_trajectories(results, modelSelector=(('per','per','avg'
                     null = False
 
                 # The hyper-parameters for the DCT model
-                hyperparameters = {"basisfunctionInds":freqs, 'starttime':results.timestamps[s][0], 
-                                   'timestep':results.meantimestepPerSeq[s], 'numsteps':results.number_of_timesteps[s]}
+                hyperparameters = freqs
                                    #starttime = hyperparameters['starttime']
-
+                modelparameters =  {'starttime':results.timestamps[s][0], 
+                                   'timestep':results.meantimestepPerSeq[s], 'numsteps':results.number_of_timesteps[s]}
                 # The parameters for the DCT filter without amplitude reduction model (a "raw" DCT filter). This
                 # is the basis for *all* the models, so we construct it and save it no matter what 
                 # estimator has been requested
                 parameters = _sig.amplitudes_at_frequencies(freqs, timeseries, transform=transform)
                 del parameters[outcomes[-1]]
                 # Creates the model, and records it into the results object with a standard 'key' to find it later.
-                ffmodel = _mdl.ProbabilityTrajectoryModel(outcomes, modeltype=transform, 
-                                                          hyperparameters=hyperparameters, parameters=parameters)
+                ffmodel = _mdl.DCTProbabilityTrajectoryModel(outcomes,
+                                                          hyperparameters=hyperparameters, parameters=parameters,
+                                                          modelparameters=modelparameters)
                 # Records this estimate.
                 results.add_reconstruction(ent, seq, ffmodel, modelSelector=modelSelector, estimator='FF-unbounded', 
                                            auxDict={'null':null}, overwrite=overwrite)

packages/pygsti/extras/drift/estimate.py

@@ -187,7 +187,7 @@ def objfunc(parameterslist):
         return trajectoryNegLogLikelihood(p, data, min_p, max_p)
 
     options = {'disp':False}
-    numparams = len(model.basisfunctionInds)*len(model.independent_outcomes)
+    numparams = len(model.hyperparameters)*len(model.independent_outcomes)
     if verbosity > 0:
         print("      - Performing MLE over {} parameters...".format(numparams),end='')
     if verbosity > 1:
@@ -229,7 +229,7 @@ def uniform_amplitude_compression(model, times, epsilon=0.001, stepsize=0.005, v
 
     """
     newmodel = model.copy()
-    if len(newmodel.basisfunctionInds) <= 1:
+    if len(newmodel.hyperparameters) <= 1:
         return model, False
 
     pt = newmodel.get_probabilities(times)

packages/pygsti/extras/drift/model.py

@@ -16,7 +16,7 @@ class ProbabilityTrajectoryModel(object):
     """
     todo
     """
-    def __init__(self, outcomes, modeltype, hyperparameters, parameters):
+    def __init__(self, outcomes, hyperparameters, parameters):
         """
         todo:
 
@@ -38,47 +38,14 @@ def __init__(self, outcomes, modeltype, hyperparameters, parameters):
         self.independent_outcomes = outcomes[:-1]
         self.constrained_outcome = outcomes[-1]
         self.numoutcomes = len(outcomes)
-        self.set_basisfunctions(modeltype, hyperparameters, parameters)
+        self.set_hyperparameters(hyperparameters, parameters)
+        #self.set_basisfunctions(modeltype, hyperparameters, parameters)
 
         return None        
 
-    def set_basisfunctions(self, modeltype, hyperparameters, parameters):
-        """
-        Todo:
-        """
-        self.modeltype = modeltype
-
-        if modeltype == 'null':
-
-            self.fullmodelsize = 1            
-            def basisfunction(i, times):
-                if i < self.fullmodelsize:
-                    return _np.array([1 for t in times]) #_np.cos(freqIns[i]*_np.pi*(t-starttime+0.5)/timedif)
-                else:
-                    raise ValueError("Invalid basis function index!")
-
-        elif modeltype == 'DCT':
-
-            try:
-                starttime = hyperparameters['starttime']
-                timestep = hyperparameters['timestep']
-                numsteps = hyperparameters['numsteps']
-                #timedif = endtime - starttime              
-            except:
-                raise ValueError("The hyperparameters are invalid for the model type! Need the start time and end time to creat the basis functions!")
-
-            self.fullmodelsize = _np.inf         
-            def basisfunction(i, times):
-                #return _np.array([_np.cos(i*_np.pi*(t-starttime+0.5)/timedif) for t in times])
-                return _np.array([_np.cos(i*_np.pi*((t-starttime)/timestep+0.5)/numsteps) for t in times])
-        else:
-            raise ValueError("Invalid model type!")
-
-        self.basisfunction = basisfunction
-        # This sets the hyperparameters and parameters so that they agree with the new model.
-        self.set_hyperparameters(hyperparameters, parameters)
+    def basisfunction(self):
 
-        return None
+        raise NotImplementedError("This should be defined in derived classes!")
 
     def get_basis_function(self, i, times):
         """
@@ -90,14 +57,16 @@ def set_hyperparameters(self, hyperparameters, parameters):
         """
         todo
         """
-        basisfunctionInds = list(hyperparameters['basisfunctionInds'])
-        basisfunctionInds.sort()
-        assert(max(basisfunctionInds) < self.fullmodelsize)
-        assert(0. in basisfunctionInds), "The constant function must be one of the basis functions!"
-        self.basisfunctionInds = basisfunctionInds
+        hyperparameters = list(hyperparameters)
+        hyperparameters.sort()
+        self.hyperparameters = hyperparameters
+        #assert(max(basisfunctionInds) < self.fullmodelsize)
+        #assert(0. in basisfunctionInds), "The constant function must be one of the basis functions!"
+        #self.basisfunctionInds = basisfunctionInds
         self.set_parameters(parameters)
-        self.modelsize_unconstrained = (self.numoutcomes-1)*len(basisfunctionInds)
-        self.modelsize_unconstrained_peroutcome= len(basisfunctionInds)
+        self.modelsize_unconstrained = (self.numoutcomes-1)*len(hyperparameters)
+        self.modelsize_unconstrained_peroutcome= len(hyperparameters)
+
         return None
 
     def set_parameters(self, parameters):
@@ -111,7 +80,7 @@ def set_parameters_from_list(self, parameterslist):
         """
         todo:
         """
-        numbasisfuncs = len(self.basisfunctionInds)
+        numbasisfuncs = len(self.hyperparameters)
         parameters = {o : parameterslist[oind:(1+oind)*numbasisfuncs] for oind,o in enumerate(self.independent_outcomes)}
         self.parameters = parameters
         return None
@@ -137,7 +106,7 @@ def get_probabilities(self, times):
         """ 
         p_for_constrained_outcome = _np.ones(len(times))
         for o in self.independent_outcomes:
-            p = _np.sum(_np.array([self.parameters[o][ind]*self.basisfunction(i,times) for ind, i in enumerate(self.basisfunctionInds)]), axis=0) 
+            p = _np.sum(_np.array([self.parameters[o][ind]*self.basisfunction(i,times) for ind, i in enumerate(self.hyperparameters)]), axis=0) 
             p_for_constrained_outcome = p_for_constrained_outcome - p
             probs = {o : p}    
         probs[self.constrained_outcome] = p_for_constrained_outcome 
@@ -147,11 +116,50 @@ def get_independent_probabilities(self, times):
         """
         todo
         """
-        return {o : _np.sum(_np.array([self.parameters[o][ind]*self.basisfunction(i,times) for ind, i in enumerate(self.basisfunctionInds)]), axis=0) for o in self.independent_outcomes}
+        return {o : _np.sum(_np.array([self.parameters[o][ind]*self.basisfunction(i,times) for ind, i in enumerate(self.hyperparameters)]), axis=0) for o in self.independent_outcomes}
 
     #def get_probabilities(times):
     #               
     #   return {o : [_np.sum(_np.array([self.parameters[o][0] + self.parameters[o][i]*_np.cos(freqIns[i]*_np.pi*(t-starttime+0.5)/timedif) for i in range(self.dof)])) for t in times]}    
 
     def copy(self):
-        return _copy.deepcopy(self)
+        return _copy.deepcopy(self)
+
+
+class NullProbabilityTrajectoryModel(ProbabilityTrajectoryModel):
+
+    def __init__(self, outcomes, hyperparameters, parameters):
+
+        self.fullmodelsize = 1
+        ProbabilityTrajectoryModel.__init__(self, outcomes, hyperparameters, parameters)          
+
+    def basisfunction(self, i, times):
+        """
+        
+        """
+        if i < self.fullmodelsize:
+            return _np.array([1 for t in times]) #_np.cos(freqIns[i]*_np.pi*(t-starttime+0.5)/timedif)
+        else:
+            raise ValueError("Invalid basis function index!")
+
+class DCTProbabilityTrajectoryModel(ProbabilityTrajectoryModel):
+
+    def __init__(self, outcomes, hyperparameters, parameters, modelparameters):
+
+        try:
+            self.starttime = modelparameters['starttime']
+            self.timestep = modelparameters['timestep']
+            self.numsteps = modelparameters['numsteps']
+            #timedif = endtime - starttime              
+        except:
+            raise ValueError("The hyperparameters are invalid for the model type! Need the start time and end time to creat the basis functions!")
+
+        self.fullmodelsize = _np.inf    
+
+        ProbabilityTrajectoryModel.__init__(self, outcomes, hyperparameters, parameters)            
+
+    def basisfunction(self, i, times):
+        """
+        Todo
+        """            
+        return _np.array([_np.cos(i*_np.pi*((t-self.starttime)/self.timestep+0.5)/self.numsteps) for t in times])

packages/pygsti/extras/rb/sample.py

@@ -786,6 +786,19 @@ def simultaneous_random_circuit(pspec, length, structure='1Q', sampler='Qelimina
 
     return circuit, idealout
 
+def _get_setting(l, circuitindex, substructure, lengths, circuits_per_length, structure):
+
+    lind = lengths.index(l)
+    settingDict = {}
+
+    for s in structure: 
+        if s in substructure:
+            settingDict[s] =  len(lengths) + lind*circuits_per_length + circuitindex
+        else:
+            settingDict[s] =  lind
+
+    return settingDict
+
 def simultaneous_random_circuits_experiment(pspec, lengths, circuits_per_length, structure='1Q', sampler='Qelimination', samplerargs=[], addlocal=False,
                                             lsargs=[], set_isolated=True, setcomplement_isolated=False,
                                             descriptor='A set of simultaneous random circuits', verbosity=1):
@@ -876,6 +889,8 @@ def simultaneous_random_circuits_experiment(pspec, lengths, circuits_per_length,
     experiment_dict['spec']['addlocal'] = addlocal
     experiment_dict['spec']['lsargs'] = lsargs
     experiment_dict['spec']['descriptor'] = descriptor
+    experiment_dict['spec']['createdby'] = 'extras.rb.sample.simultaneous_random_circuits_experiment'
+
 
     if isinstance(structure,str):
         assert(structure == '1Q'), "The only default `structure` option is the string '1Q'"
@@ -895,6 +910,7 @@ def simultaneous_random_circuits_experiment(pspec, lengths, circuits_per_length,
     experiment_dict['spec']['structure'] = structure
     experiment_dict['circuits'] = {}
     experiment_dict['probs'] = {}
+    experiment_dict['settings'] = {}
 
     for lnum, l in enumerate(lengths):
         if verbosity > 0:
@@ -907,23 +923,36 @@ def simultaneous_random_circuits_experiment(pspec, lengths, circuits_per_length,
             if (not set_isolated) and (not setcomplement_isolated):
                 experiment_dict['circuits'][l,j] = circuit
                 experiment_dict['probs'][l,j] = idealout
+                experiment_dict['settings'][l,j] = {s: len(lengths) + lnum*circuits_per_length + j for s in tuple(structure)}
             else:
                 experiment_dict['circuits'][l,j] = {}
                 experiment_dict['probs'][l,j] = {}
+                experiment_dict['settings'][l,j] = {}
                 experiment_dict['circuits'][l,j][tuple(structure)] = circuit
-                experiment_dict['probs'][l,j][tuple(structure)] = idealout
-
+                experiment_dict['probs'][l,j][tuple(structure)] = idealout           
+                experiment_dict['settings'][l,j][tuple(structure)] = _get_setting(l, j, structure, lengths, circuits_per_length, 
+                                                                                    structure)
             if set_isolated:
                 for subset_ind, subset in enumerate(structure):
                     subset_circuit = circuit.copy(editable=True)
-                    print(subset)
+                    #print(subset)
                     for q in circuit.line_labels:
                         if q not in subset:
-                            print(subset_circuit, q)
+                            #print(subset_circuit, q)
                             subset_circuit.replace_with_idling_line(q)
                     subset_circuit.done_editing()
                     experiment_dict['circuits'][l,j][(tuple(subset),)] = subset_circuit
                     experiment_dict['probs'][l,j][(tuple(subset),)] = idealout[subset_ind]
+                    # setting = {}
+                    # for s in structure: 
+                    #     if s in subset:
+                    #         setting[s] =  len(lengths) + lnum*circuits_per_length + j
+                    #     else:
+                    #         setting[s] =  lnum
+                    experiment_dict['settings'][l,j][(tuple(subset),)] = _get_setting(l, j, (tuple(subset),), lengths, circuits_per_length, 
+                                                                                    structure)
+                    # print(subset)
+                    # print(_get_setting(l, j, subset, lengths, circuits_per_length, structure))
 
             if setcomplement_isolated:
                 for subset_ind, subset in enumerate(structure):
@@ -941,6 +970,15 @@ def simultaneous_random_circuits_experiment(pspec, lengths, circuits_per_length,
                         experiment_dict['circuits'][l,j][subsetcomplement] = subsetcomplement_circuit
                         experiment_dict['probs'][l,j][subsetcomplement] = subsetcomplement_idealout
 
+                        # for s in structure: 
+                        #     if s in subsetcomplement:
+                        #         setting[s] =  len(lengths) + lnum*circuits_per_length + j
+                        #     else:
+                        #         setting[s] =  lnum
+                        experiment_dict['settings'][l,j][subsetcomplement] = _get_setting(l, j, subsetcomplement, lengths, circuits_per_length, 
+                                                                                    structure)
+
+
             if verbosity > 0: print(j+1,end=',')
         if verbosity >0: print('')
 
@@ -1636,7 +1674,7 @@ def simultaneous_direct_rb_circuit(pspec, length, structure='1Q', sampler='Qelim
     inversion_circuit.done_editing()
 
     if cliffordtwirl:
-        full_circuit = initial_circuit.copy()
+        full_circuit = initial_circuit.copy(editable=True)
         full_circuit.append_circuit(circuit)
         full_circuit.append_circuit(inversion_circuit)
     else:
@@ -1845,6 +1883,7 @@ def simultaneous_direct_rb_experiment(pspec, lengths, circuits_per_length, struc
     experiment_dict['spec']['compilerargs'] = compilerargs
     experiment_dict['spec']['partitioned'] = partitioned
     experiment_dict['spec']['descriptor'] = descriptor
+    experiment_dict['spec']['createdby'] = 'extras.rb.sample.simultaneous_direct_rb_experiment'
 
     if isinstance(structure,str):
         assert(structure == '1Q'), "The only default `structure` option is the string '1Q'"
@@ -1864,6 +1903,7 @@ def simultaneous_direct_rb_experiment(pspec, lengths, circuits_per_length, struc
     experiment_dict['spec']['structure'] = structure
     experiment_dict['circuits'] = {}
     experiment_dict['idealout'] = {}
+    experiment_dict['settings'] = {}
 
     for subsetQs in structure:
         subgraph = pspec.qubitgraph.subgraph(list(subsetQs))
@@ -1886,24 +1926,32 @@ def simultaneous_direct_rb_experiment(pspec, lengths, circuits_per_length, struc
             else:
                 experiment_dict['circuits'][l,j] = {}
                 experiment_dict['idealout'][l,j] = {}
+                experiment_dict['settings'][l,j] = {}
                 experiment_dict['circuits'][l,j][tuple(structure)] = circuit
                 experiment_dict['idealout'][l,j][tuple(structure)] = idealout 
+                experiment_dict['settings'][l,j][tuple(structure)] = _get_setting(l, j, structure, lengths, circuits_per_length, 
+                                                                                    structure)
 
             if set_isolated:
                 for subset_ind, subset in enumerate(structure):
-                    subset_circuit = circuit.copy()
+                    subset_circuit = circuit.copy(editable=True)
                     for q in circuit.line_labels:
                         if q not in subset:
                             subset_circuit.replace_with_idling_line(q)
+                    subset_circuit.done_editing()
                     experiment_dict['circuits'][l,j][(tuple(subset),)] = subset_circuit
                     experiment_dict['idealout'][l,j][(tuple(subset),)] = (idealout[subset_ind],)
+                    experiment_dict['settings'][l,j][(tuple(subset),)] = _get_setting(l, j, (tuple(subset),), lengths, circuits_per_length, 
+                                                                                    structure)
+
 
             if setcomplement_isolated:
                 for subset_ind, subset in enumerate(structure):
-                        subsetcomplement_circuit = circuit.copy()
+                        subsetcomplement_circuit = circuit.copy(editable=True)
                         for q in circuit.line_labels:
                             if q in subset:
                                 subsetcomplement_circuit.replace_with_idling_line(q)
+                        subsetcomplement_circuit.done_editing()
                         subsetcomplement = list(_copy.copy(structure))
                         subsetcomplement_idealout = list(_copy.copy(idealout))
                         del subsetcomplement[subset_ind]
@@ -1912,6 +1960,8 @@ def simultaneous_direct_rb_experiment(pspec, lengths, circuits_per_length, struc
                         subsetcomplement_idealout  = tuple(subsetcomplement_idealout)
                         experiment_dict['circuits'][l,j][subsetcomplement] = subsetcomplement_circuit
                         experiment_dict['idealout'][l,j][subsetcomplement] = subsetcomplement_idealout
+                        experiment_dict['settings'][l,j][subsetcomplement] = _get_setting(l, j, subsetcomplement, lengths, circuits_per_length, 
+                                                                                    structure)
 
             if verbosity > 0: print(j+1,end=',')
         if verbosity >0: print('')