Merge branch 'develop' into simplify-TimeIndependentMDCObjectiveFunction

Author: rileyjmurray

.github/CODEOWNERS

@@ -37,7 +37,6 @@ pygsti/algorithms/compilers.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeep
 pygsti/algorithms/mirroring.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers
 pygsti/algorithms/randomcircuit.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers
 pygsti/algorithms/rbfit.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers
-pygsti/extras/rb.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers # Should this just be deprecated and removed?
 pygsti/protocols/rb.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers
 pygsti/tools/rbtheory.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers
 pygsti/tools/rbtools.py @sandialabs/pygsti-rb @sandialabs/pygsti-gatekeepers

pygsti/algorithms/core.py

@@ -31,8 +31,7 @@
 from pygsti.modelmembers import states as _state
 from pygsti.circuits.circuitlist import CircuitList as _CircuitList
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
-from pygsti.optimize.customlm import CustomLMOptimizer as _CustomLMOptimizer
-from pygsti.optimize.customlm import Optimizer as _Optimizer
+from pygsti.optimize.simplerlm import Optimizer as _Optimizer, SimplerLMOptimizer as _SimplerLMOptimizer
 from pygsti import forwardsims as _fwdsims
 from pygsti import layouts as _layouts
 
@@ -619,7 +618,7 @@ def run_gst_fit_simple(dataset, start_model, circuits, optimizer, objective_func
     model : Model
         the best-fit model.
     """
-    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _CustomLMOptimizer.cast(optimizer)
+    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _SimplerLMOptimizer.cast(optimizer)
     objective_function_builder = _objfns.ObjectiveFunctionBuilder.cast(objective_function_builder)
     array_types = optimizer.array_types + \
         objective_function_builder.compute_array_types(optimizer.called_objective_methods, start_model.sim)
@@ -666,7 +665,7 @@ def run_gst_fit(mdc_store, optimizer, objective_function_builder, verbosity=0):
     objfn_store : MDCObjectiveFunction
         the objective function and store containing the best-fit model evaluated at the best-fit point.
     """
-    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _CustomLMOptimizer.cast(optimizer)
+    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _SimplerLMOptimizer.cast(optimizer)
     comm = mdc_store.resource_alloc.comm
     profiler = mdc_store.resource_alloc.profiler
     printer = VerbosityPrinter.create_printer(verbosity, comm)
@@ -843,7 +842,7 @@ def iterative_gst_generator(dataset, start_model, circuit_lists,
           (an "evaluated" model-dataset-circuits store).
     """
     resource_alloc = _ResourceAllocation.cast(resource_alloc)
-    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _CustomLMOptimizer.cast(optimizer)
+    optimizer = optimizer if isinstance(optimizer, _Optimizer) else _SimplerLMOptimizer.cast(optimizer)
     comm = resource_alloc.comm
     profiler = resource_alloc.profiler
     printer = VerbosityPrinter.create_printer(verbosity, comm)
@@ -1150,8 +1149,6 @@ def find_closest_unitary_opmx(operation_mx):
     # d = _np.sqrt(operation_mx.shape[0])
     # I = _np.identity(d)
 
-    #def getu_1q(basisVec):  # 1 qubit version
-    #    return _spl.expm( 1j * (basisVec[0]*_tools.sigmax + basisVec[1]*_tools.sigmay + basisVec[2]*_tools.sigmaz) )
     def _get_gate_mx_1q(basis_vec):  # 1 qubit version
         return _tools.single_qubit_gate(basis_vec[0],
                                         basis_vec[1],

pygsti/algorithms/fiducialselection.py

@@ -409,13 +409,6 @@ def final_result_test(final_fids, verb_printer):
         return prepFidList, measFidList    
 
 
-#def bool_list_to_ind_list(boolList):
-#    output = _np.array([])
-#    for i, boolVal in boolList:
-#        if boolVal == 1:
-#            output = _np.append(i)
-#    return output
-
 def xor(*args):
     """
     Implements logical xor function for arbitrary number of inputs.

pygsti/algorithms/gaugeopt.py

@@ -290,7 +290,7 @@ def gaugeopt_custom(model, objective_fn, gauge_group=None,
     gaugeGroupEl = gauge_group.compute_element(x0)  # re-used element for evals
 
     def _call_objective_fn(gauge_group_el_vec, oob_check=False):
-        # Note: oob_check can be True if oob_check_interval>=1 is given to the custom_leastsq below
+        # Note: oob_check can be True if oob_check_interval>=1 is given to the simplish_leastsq below
         gaugeGroupEl.from_vector(gauge_group_el_vec)
         return objective_fn(gaugeGroupEl, oob_check)
 
@@ -309,7 +309,7 @@ def _call_jacobian_fn(gauge_group_el_vec):
         assert(_call_jacobian_fn is not None), "Cannot use 'ls' method unless jacobian is available"
         ralloc = _baseobjs.ResourceAllocation(comm)  # FUTURE: plumb up a resource alloc object?
         test_f = _call_objective_fn(x0)
-        solnX, converged, msg, _, _, _, _, _ = _opt.custom_leastsq(
+        solnX, converged, msg, _, _, _, _ = _opt.simplish_leastsq(
             _call_objective_fn, _call_jacobian_fn, x0, f_norm2_tol=tol,
             jac_norm_tol=tol, rel_ftol=tol, rel_xtol=tol,
             max_iter=maxiter, resource_alloc=ralloc,

pygsti/algorithms/germselection.py

@@ -27,6 +27,7 @@
 from pygsti.baseobjs.statespace import ExplicitStateSpace as _ExplicitStateSpace
 from pygsti.baseobjs.statespace import QuditSpace as _QuditSpace
 from pygsti.models import ExplicitOpModel as _ExplicitOpModel
+from pygsti.forwardsims import MatrixForwardSimulator as _MatrixForwardSimulator
 
 FLOATSIZE = 8  # in bytes: TODO: a better way
 
@@ -57,10 +58,8 @@ def find_germs(target_model, randomize=True, randomization_strength=1e-2,
 
     Parameters
     ----------
-    target_model : Model or list of Model
-        The model you are aiming to implement, or a list of models that are
-        copies of the model you are trying to implement (either with or
-        without random unitary perturbations applied to the models).
+    target_model : Model
+        The model you are aiming to implement.
 
     randomize : bool, optional
         Whether or not to add random unitary perturbations to the model(s)
@@ -188,8 +187,14 @@ def find_germs(target_model, randomize=True, randomization_strength=1e-2,
         A list containing the germs making up the germ set.
     """
     printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, comm)
+    
+    if not isinstance(target_model.sim, _MatrixForwardSimulator):
+        target_model = target_model.copy()
+        target_model.sim = 'matrix'
+
     modelList = _setup_model_list(target_model, randomize,
                                   randomization_strength, num_gs_copies, seed)
+                                  
     gates = list(target_model.operations.keys())
     availableGermsList = []
     if candidate_germ_counts is None: candidate_germ_counts = {6: 'all upto'}
@@ -1351,6 +1356,10 @@ def test_germ_set_finitel(model, germs_to_test, length, weights=None,
         eigenvalues (from small to large) of the jacobian^T * jacobian
         matrix used to determine parameter amplification.
     """
+    if not isinstance(model.sim, _MatrixForwardSimulator):
+        model = model.copy()
+        model.sim = 'matrix'
+
     # Remove any SPAM vectors from model since we only want
     # to consider the set of *gate* parameters for amplification
     # and this makes sure our parameter counting is correct
@@ -3295,80 +3304,81 @@ def symmetric_low_rank_spectrum_update(update, orig_e, U, proj_U, force_rank_inc
     #return the new eigenvalues
     return new_evals, True
 
-#Note: This function won't work for our purposes because of the assumptions
-#about the rank of the update on the nullspace of the matrix we're updating,
-#but keeping this here commented for future reference.
-#Function for doing fast calculation of the updated inverse trace:
-#def riedel_style_inverse_trace(update, orig_e, U, proj_U, force_rank_increase=True):
-#    """
-#    input:
-#    
-#    update : ndarray
-#        symmetric low-rank update to perform.
-#        This is the first half the symmetric rank decomposition s.t.
-#        [email protected]= the full update matrix.
-#    
-#    orig_e : ndarray
-#        Spectrum of the original matrix. This is a 1-D array.
-#        
-#    proj_U : ndarray
-#        Projector onto the complement of the column space of the
-#        original matrix's eigenvectors.
-#        
-#    output:
-#    
-#    trace : float
-#        Value of the trace of the updated psuedoinverse matrix.
-#    
-#    updated_rank : int
-#        total rank of the updated matrix.
-#        
-#    rank_increase_flag : bool
-#        a flag that is returned to indicate is a candidate germ failed to amplify additional parameters. 
-#        This indicates things short circuited and so the scoring function should skip this germ.
-#    """
-#    
-#    #First we need to for the matrix P, whose column space
-#    #forms an orthonormal basis for the component of update
-#    #that is in the complement of U.
-#    
-#    proj_update= proj_U@update
-#    
-#    #Next take the RRQR decomposition of this matrix:
-#    q_update, r_update, _ = _sla.qr(proj_update, mode='economic', pivoting=True)
-#    
-#    #Construct P by taking the columns of q_update corresponding to non-zero values of r_A on the diagonal.
-#    nonzero_indices_update= _np.nonzero(_np.diag(r_update)>1e-10) #HARDCODED (threshold is hardcoded)
-#    
-#    #if the rank doesn't increase then we can't use the Riedel approach.
-#    #Abort early and return a flag to indicate the rank did not increase.
-#    if len(nonzero_indices_update[0])==0 and force_rank_increase:
-#        return None, None, False
-#    
-#    P= q_update[: , nonzero_indices_update[0]]
-#    
-#    updated_rank= len(orig_e)+ len(nonzero_indices_update[0])
-#    
-#    #Now form the matrix R_update which is given by P.T @ proj_update.
-#    R_update= P.T@proj_update
-#    
-#    #R_update gets concatenated with U.T@update to form
-#    #a block column matrixblock_column= np.concatenate([U.T@update, R_update], axis=0)    
-#    
-#    Uta= U.T@update
-#    
-#    try:
-#        RRRDinv= R_update@_np.linalg.inv(R_update.T@R_update) 
-#    except _np.linalg.LinAlgError as err:
-#        print('Numpy thinks this matrix is singular, condition number is: ', _np.linalg.cond(R_update.T@R_update))
-#        print((R_update.T@R_update).shape)
-#        raise err
-#    pinv_orig_e_mat= _np.diag(1/orig_e)
-#    
-#    trace= _np.sum(1/orig_e) + _np.trace( RRRDinv@(_np.eye(Uta.shape[1]) + Uta.T@pinv_orig_e_mat@Uta)@RRRDinv.T )
-#    
-#    return trace, updated_rank, True
+# Note: Th function below won't work for our purposes because of the assumptions
+# about the rank of the update on the nullspace of the matrix we're updating,
+# but keeping this here commented for future reference.
+'''
+def riedel_style_inverse_trace(update, orig_e, U, proj_U, force_rank_increase=True):
+    """
+    input:
+    
+    update : ndarray
+        symmetric low-rank update to perform.
+        This is the first half the symmetric rank decomposition s.t.
+        [email protected]= the full update matrix.
+    
+    orig_e : ndarray
+        Spectrum of the original matrix. This is a 1-D array.
+        
+    proj_U : ndarray
+        Projector onto the complement of the column space of the
+        original matrix's eigenvectors.
+        
+    output:
+    
+    trace : float
+        Value of the trace of the updated psuedoinverse matrix.
+    
+    updated_rank : int
+        total rank of the updated matrix.
+        
+    rank_increase_flag : bool
+        a flag that is returned to indicate is a candidate germ failed to amplify additional parameters. 
+        This indicates things short circuited and so the scoring function should skip this germ.
+    """
     
+    #First we need to for the matrix P, whose column space
+    #forms an orthonormal basis for the component of update
+    #that is in the complement of U.
+    
+    proj_update= proj_U@update
+    
+    #Next take the RRQR decomposition of this matrix:
+    q_update, r_update, _ = _sla.qr(proj_update, mode='economic', pivoting=True)
+    
+    #Construct P by taking the columns of q_update corresponding to non-zero values of r_A on the diagonal.
+    nonzero_indices_update= _np.nonzero(_np.diag(r_update)>1e-10) #HARDCODED (threshold is hardcoded)
+    
+    #if the rank doesn't increase then we can't use the Riedel approach.
+    #Abort early and return a flag to indicate the rank did not increase.
+    if len(nonzero_indices_update[0])==0 and force_rank_increase:
+        return None, None, False
+    
+    P= q_update[: , nonzero_indices_update[0]]
+    
+    updated_rank= len(orig_e)+ len(nonzero_indices_update[0])
+    
+    #Now form the matrix R_update which is given by P.T @ proj_update.
+    R_update= P.T@proj_update
+    
+    #R_update gets concatenated with U.T@update to form
+    #a block column matrixblock_column= np.concatenate([U.T@update, R_update], axis=0)    
+    
+    Uta= U.T@update
+    
+    try:
+        RRRDinv= R_update@_np.linalg.inv(R_update.T@R_update) 
+    except _np.linalg.LinAlgError as err:
+        print('Numpy thinks this matrix is singular, condition number is: ', _np.linalg.cond(R_update.T@R_update))
+        print((R_update.T@R_update).shape)
+        raise err
+    pinv_orig_e_mat= _np.diag(1/orig_e)
+    
+    trace= _np.sum(1/orig_e) + _np.trace( RRRDinv@(_np.eye(Uta.shape[1]) + Uta.T@pinv_orig_e_mat@Uta)@RRRDinv.T )
+    
+    return trace, updated_rank, True
+'''
+
 def minamide_style_inverse_trace(update, orig_e, U, proj_U, force_rank_increase=False):
     """
     This function performs a low-rank update to the components of