added an updated system for statreprep and fixed some bugs along the way

examples/gradient_test.py

@@ -2,23 +2,32 @@
 import numpy as np
 from qsearch.gates import *
 
-gate = SingleQutritStep()
+gate = ProductGate(KroneckerGate(U3Gate(), U3Gate(), U3Gate()), KroneckerGate(IdentityGate(), ProductGate(CNOTGate(), KroneckerGate(XZXZGate(), U3Gate()))), KroneckerGate(ProductGate(CNOTGate(), KroneckerGate(XZXZGate(), U3Gate())), IdentityGate()), KroneckerGate(IdentityGate(), ProductGate(CNOTGate(), KroneckerGate(XZXZGate(), U3Gate()))))
 totaldiff = [0] * gate.num_inputs
-eps = 5e-11
+eps = 5e-5
+
+
+def func(v):
+    return qsearch.utils.distance_with_initial_state(np.array([1,0,0,0,0,0,0,0]),np.array([0.5,0,0.5,0,0.5,0,0,0.5]), np.eye(8),gate.matrix(v))
+
+def jac_func(v):
+    return qsearch.utils.distance_with_initial_state_jac(np.array([1,0,0,0,0,0,0,0]),np.array([0.5,0,0.5,0,0.5,0,0,0.5]), np.eye(8),*gate.mat_jac(v))[1]
 
 for _ in range(100):
     v = np.random.rand(gate.num_inputs)
-    M, Js = gate.mat_jac(v)
+    M = func(v)
+    Js = jac_func(v)
     for i in range(gate.num_inputs):
         v2 = np.copy(v)
         v2[i] = v[i] + eps
-        U1 = gate.matrix(v2)
+        U1 = func(v2)
         v2[i] = v[i] - eps
-        U2 = gate.matrix(v2)
+        U2 = func(v2)
 
         FD = (U1 - U2) / (2*eps)
 
         diffs = np.abs(np.sum(FD - Js[i]))
         totaldiff[i] += diffs
 
+print(np.array(totaldiff) < 2*eps)
 print(totaldiff)

examples/stateprep.py

@@ -8,14 +8,21 @@
 p = qsearch.Project("stateprep-example")
 
 # configure the project with the stateprep defaults instead of the standard synthesis defaults
-p.configure(**qsearch.defaults.stateprep_defaults)
-p["compiler_class"] = leap_compiler.LeapCompiler
-p["solver"] = qsearch.solvers.LeastSquares_Jac_Solver()
+
+stateprep_options = qsearch.Options(smart_defaults=qsearch.defaults.stateprep_smart_defaults)
+
 
 # add states that are converted using generate_stateprep_target_matrix
-p.add_compilation("basic_state_test", qsearch.utils.generate_stateprep_target_matrix([0.5,0,0,0.5j,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.5j,0,0,0,0,0,0,0,-0.5,0,0,0,0,0]))
+
+# It may seem strange that we have to pass an identity.  Qsearch is still doing unitary synthesis;
+# stateprep is achieved by modifying the way we compare unitaries such that they are compared by
+# the state produced by acting on an initial state.
+# In this case, we are synthesizing a circuit designed to act on the |000> state to match the state
+# resulting by performing the Identity on the desired state.
+# Note that you can specify the intiial state for the circuit as initial_state (the default state is
+# the zero state that matches target_state in number of qudits).
+toffoli_magic_state = np.array([0.5,0,0.5,0,0.5,0,0,0.5], dtype='complex128')
+p.add_compilation("toffoli_magic_state", np.eye(8,dtype='complex128'), target_state=toffoli_magic_state, options=stateprep_options)
 
 p.run()
 
-# run post-processing to improve circuits that were generated with LEAP
-p.post_process(reoptimizing_compiler.ReoptimizingCompiler(), solver=multistart_solvers.MultiStart_Solver(12), parallelizer=parallelizers.ProcessPoolParallelizer, depth=7)

qsearch/compiler.py

@@ -50,10 +50,11 @@ def __init__(self, options=Options()):
         Args:
             options: See class level documentation for the options SearchCompiler uses
         """
-        self.options = options
+        self.options = Options()
         self.options.set_defaults(**standard_defaults)
         self.options.set_defaults(verbosity=1, logfile=None, stdout_enabled=True)
         self.options.set_smart_defaults(**standard_smart_defaults)
+        self.options = self.options.updated(options)
 
     def compile(self, options=Options()):
         """

qsearch/defaults.py

@@ -9,15 +9,13 @@
     stateprep_defaults : A dictionary containing defaults for stateprep synthesis.
 """
 
-from . import utils, gatesets, solvers, backends, parallelizers, heuristics, logging, checkpoints, assemblers
+from . import utils, gatesets, solvers, backends, parallelizers, heuristics, logging, checkpoints, assemblers, compiler
 from functools import partial
+import numpy as np
 
 
 
 def default_eval_func(options):
-    if options.error_func == utils.matrix_residuals:
-        return utils.matrix_distance_squared
-    else:
         return options.error_func
 
 def default_heuristic(options):
@@ -47,6 +45,37 @@ def default_logger(options):
 def default_checkpoint(options):
     return checkpoints.FileCheckpoint(options=options)
 
+def stateprep_error_func(options):
+    return partial(utils.distance_with_initial_state,options.target_state,options.initial_state)
+
+def stateprep_error_jac(options):
+    return partial(utils.distance_with_initial_state_jac,options.target_state,options.initial_state)
+
+def stateprep_error_resi(options):
+    return partial(utils.residuals_with_initial_state,options.target_state,options.initial_state)
+
+def stateprep_error_resi_jac(options):
+    return partial(utils.residuals_with_initial_state_jac,options.target_state,options.initial_state)
+
+def stateprep_initial_state(options):
+    v = np.zeros(options.target_state.shape,dtype='complex128')
+    v[0] = 1
+    return v
+
+def stateprep_target(options):
+    return np.eye(options.target_state.shape[0], dtype='complex128')
+
+def stateprep_default_solver(options):
+    opt = options.copy()
+    opt.make_required("error_jac", "error_func")
+    if "error_func" in opt and "error_jac" not in opt:
+        return solvers.COBYLA_Solver()
+    else:
+        return solvers.BFGS_Jac_Solver()
+
+def default_compiler(options):
+    return compiler.SearchCompiler # this gets around some pesky import loops
+
 def identity(U):
     return U
 
@@ -68,6 +97,9 @@ def identity(U):
         "write_location" : None,
         "unitary_preprocessor": utils.nearest_unitary,
         "timeout" : float('inf'),
+        "blas_threads" : None,
+        "verbosity" : 1,
+        "stdout_enabled" : True,
         }
 standard_smart_defaults = {
         "eval_func":default_eval_func,
@@ -77,12 +109,16 @@ def identity(U):
         "heuristic":default_heuristic,
         "logger" :default_logger,
         "checkpoint":default_checkpoint,
+        "compiler_class" : default_compiler,
         }
 
-stateprep_defaults = {
-        "error_residuals" : partial(utils.matrix_residuals_slice, (0, slice(None))),
-        "error_residuals_jac" : partial(utils.matrix_residuals_slice_jac, (0, slice(None))),
-        "eval_func" : partial(utils.eval_func_from_residuals, partial(utils.matrix_residuals_slice, (0, slice(None)))),
-        "unitary_preprocessor": identity
+stateprep_smart_defaults = {
+        "error_residuals" : stateprep_error_resi,
+        "error_residuals_jac" : stateprep_error_resi_jac,
+        "error_func" : stateprep_error_func,
+        "error_jac" : stateprep_error_jac,
+        "initial_state" : stateprep_initial_state,
+        "target" : stateprep_target,
+        "solver" : stateprep_default_solver,
         }
 

qsearch/options.py

@@ -58,10 +58,11 @@ def __getattr__(self, name):
             # to generate an error rather than recurse infinitely when there are interdependent smart default functions
             smartfunc = self.smart_defaults[name]
             self.required.add(name)
-            retval = smartfunc(self)
-            self.smart_defaults[name] = smartfunc
-            self.cache[name] = retval
-            self.required.remove(name)
+            try:
+                retval = smartfunc(self)
+                self.cache[name] = retval
+            finally:
+                self.required.remove(name)
             return retval
         elif name in self.defaults:
             return self.defaults[name]
@@ -75,6 +76,8 @@ def __setattr__(self, name, value):
     def __contains__(self, name):
         if name in self.__dict__:
             return True
+        elif name in self.required:
+            return False
         elif name in self.defaults:
             return True
         elif name in self.smart_defaults:
@@ -84,22 +87,22 @@ def __contains__(self, name):
 
     def empty_copy(self):
         """Create an Options object with the same defaults but without any specific values."""
-        newOptions = Options(**self.defaults)
+        newOptions = Options(self.defaults)
         newOptions.smart_defaults = self.smart_defaults
         newOptions.required = self.required.copy()
         return newOptions
 
     def copy(self):
         """Create a full copy of an Options object."""
-        newOptions = Options(**self.defaults)
+        newOptions = Options(self.defaults)
         newOptions.smart_defaults.update(self.smart_defaults)
         newOptions._update_dict(self.__dict__)
         newOptions.required = self.required.copy()
         newOptions.cache = self.cache.copy()
         return newOptions
 
     def __copy__(self):
-        newOptions = Options(**self.defaults)
+        newOptions = Options(self.defaults)
         newOptions.smart_defaults.update(self.smart_defaults)
         newOptions._update_dict(self.__dict__)
         newOptions.required = self.required.copy()

qsearch/project.py

@@ -49,14 +49,14 @@ def __init__(self, path, use_mpi=False):
                     os.mkdir(path)
                 with open(self.projfile, "rb") as projfile:
                     self._compilations, self.options = pickle.load(projfile)
-                    self.set_defaults()
                     self.logger = logging.Logger(self.options.stdout_enabled, os.path.join(path, "{}-project-log.txt".format(self.name)), self.options.verbosity)
                     self.logger.logprint("Successfully loaded project {}".format(self.name))
                     self.status(logger=self.logger)
             except IOError:
                 self._compilations = dict()
                 self.options = Options()
                 self.set_defaults()
+                self.set_smart_defaults()
                 self.logger = logging.Logger(True, os.path.join(path, "{}-project-log.txt".format(self.name)), verbosity=1)
 
     def _save(self):
@@ -91,7 +91,6 @@ def add_compilation(self, name, U, options=None, handle_existing=None, **extraar
             elif s == Project_Status.PROGRESS or s == Project_Status.COMPLETE:
                 warn("A compilation with name {} already exists.  To change it, remove it and then add it again.".format(name), RuntimeWarning, stacklevel=2)
                 return
-
         compopt = Options(statefile=self._checkpoint_path(name))
         compopt.update(options, **extraargs)
         compopt.target = U
@@ -196,10 +195,13 @@ def __exit__(self, exc_typ, exc_val, exc_tb):
         else:
             sys.exit(0)
 
-    def set_defaults(self):
-        """Updates the Project Options with the standard defaults from defaults.py"""
-        self.options.set_defaults(verbosity=1,stdout_enabled=True,blas_threads=None,compiler_class=SearchCompiler,**standard_defaults)
-        self.options.set_smart_defaults(**standard_smart_defaults)
+    def set_defaults(self, defaults=standard_defaults):
+        """Updates the Project Options with the standard defaults from defaults.py, or a provided dictionary."""
+        self.options.set_defaults(**defaults)
+
+    def set_smart_defaults(self, smart_defaults=standard_smart_defaults):
+        """Updates the Project Options with the standard smart_defaults from defaults.py, or a provided dictionary"""
+        self.options.set_smart_defaults(**smart_defaults)
 
     def run(self):
         """Runs all of the compilations in the Project."""

qsearch/solvers.py

@@ -18,16 +18,15 @@
 def default_solver(options, x0=None):
     """Runs a complex list of tests to determine the best Solver for a specific situation."""
     options = options.copy()
-    # re-route the default behavior for error_func and error_jac because the default functions for those parameters often rely on the return valye from default_solver
-    options.set_defaults(logger=Logger(), U=np.array([]), error_func=None, error_jac=None, target=None)
-    options.remove_smart_defaults("error_func", "error_jac")
+    options.make_required("error_func", "error_residuals", "error_jac", "error_residuals_jac")
+
 
     # Choosse the best default solver for the given gateset
     ls_failed = False
 
     # check if Rust works on the layers
     gateset = options.gateset
-    qudits = 0 if options.target is None else int(np.log(options.target.shape[0]) // np.log(gateset.d))
+    qudits = 0 if "target" not in options else int(np.log(options.target.shape[0]) // np.log(gateset.d))
 
     rs_failed = True
     if native_from_object is not None:
@@ -41,18 +40,17 @@ def default_solver(options, x0=None):
             rs_failed = False
 
     # Check to see if the gateset and error func are explicitly supported by LeastSquares
-    error_func = options.error_func
-    error_jac = options.error_jac
     logger = options.logger
 
     if type(gateset).__module__ != QubitCNOTLinear.__module__:
         ls_failed = True
-    elif error_func is not None and (error_func.__module__ != utils.matrix_distance_squared.__module__ or (error_func.__name__ != utils.matrix_distance_squared.__name__ and error_func.__name__ != utils.matrix_residuals.__name__)):
+
+    if "error_func" in options and "error_residuals" not in options:
         ls_failed = True
 
     if not ls_failed:
         # since all provided gatesets support jacobians, this is the only check we need
-        if rs_failed or options.error_residuals not in (utils.matrix_residuals, matrix_residuals) or options.error_residuals_jac not in (utils.matrix_residuals_jac, matrix_residuals_jac):
+        if rs_failed or "error_residuals" not in options or "error_residuals_jac" not in options:
             logger.logprint("Smart default chose LeastSquares_Jac_Solver", verbosity=3)
             return LeastSquares_Jac_Solver()
         else:
@@ -71,7 +69,7 @@ def default_solver(options, x0=None):
         except:
             jac_failed = True
 
-    if error_jac is None and error_func not in [None, utils.matrix_distance_squared, utils.matrix_residuals]:
+    if "error_func" in options and not "error_jac" in options: 
         jac_failed = True
 
     if jac_failed:
@@ -162,7 +160,7 @@ def solve_for_unitary(self, circuit, options, x0=None):
         def eval_func(v):
             M, jacs = circuit.mat_jac(v)
             return options.error_jac(options.target, M, jacs)
-        result = sp.optimize.minimize(eval_func, np.random.rand(circuit.num_inputs)*2*np.pi if x0 is None else x0, method='BFGS', jac=True)
+        result = sp.optimize.minimize(eval_func, np.random.rand(circuit.num_inputs)*2*np.pi if x0 is None else x0, method='BFGS', jac=True,options={"gtol":options.threshold*0.1})
         xopt = result.x
         return (circuit.matrix(xopt), xopt)
 

qsearch/utils.py

@@ -108,6 +108,28 @@ def matrix_residuals_blacklist(badrows, badcols, A, B, I):
     Im = np.reshape(Im, (1,-1))
     return np.append(Re, Im)
 
+def distance_with_initial_state(stateA,stateB,A,B):
+    return np.abs(1-np.abs(np.vdot(A.dot(stateA),B.dot(stateB))))
+
+def distance_with_initial_state_jac(stateA,stateB,A,B,J):
+    si = A.dot(stateA)
+    s = np.vdot(si,B.dot(stateB))
+    dist = 1-np.abs(s)
+    vu = np.array([np.vdot(si, K.dot(stateB)) for K in J])
+    jacs = -(np.real(s)*np.real(vu) + np.imag(s)*np.imag(vu))/np.abs(s)
+    return (dist, jacs)
+
+def residuals_with_initial_state(stateA, stateB,A,B,I):
+    v = 1-np.conj(A.dot(stateA)) * B.dot(stateB)
+    Re, Im = np.real(v), np.imag(v)
+    return np.append(Re,Im)
+
+def residuals_with_initial_state_jac(stateA,stateB, U, M, J):
+    vc = np.conj(U.dot(stateA))
+    vu = [-vc * K.dot(stateB) for K in J]
+    vu = np.array([np.append(np.real(v),np.imag(v)) for v in vu])
+    return vu.T
+
 def matrix_residuals_blacklist_jac(slices, A, B, J):
     JU = np.array([np.append(np.reshape(np.real(np.delete(np.delete(K, badrows, 0), badcols, 1)), (1,-1)), np.reshape(np.imag(np.delete(np.delete(K, badrows, 0), badcols, 1)), (1,-1))) for K in J])
     return JU.T