Added command line support for symmetric QSP methods with symbolic fu…

…nction targets.

pyqsp/angle_sequence.py

@@ -218,7 +218,7 @@ def QuantumSignalProcessingPhases(
         # We could allow for the zero polynomial, but choose not to.
         if (is_even and is_odd) or (not is_even and not is_odd):
             raise AngleFindingError(
-                "Polynomial must have definite parity and be non-zero: {}".format(str(pcoefs)))
+                "Polynomial must have definite parity and be non-zero: {}".format(str(coefs)))
 
         # Set parity, and reduce Chebyshev coefs accordingly.
         parity = 0 if is_even else 1
@@ -231,7 +231,7 @@ def QuantumSignalProcessingPhases(
         Note that this method currently ignores choice of signal,
         measurement, and eps, suc, tolerance. These can eventually be matched with the known internal parameters for the sym_qsp method, but many of them are instroduced for numerical stability in root finding methods, which are deprecated here.
 
-        TODO: we also need to provide an internal check for good
+        TODO: we can also provide an internal check for good
         agreement with the requested polynomial.
 
         adat = np.linspace(-1., 1., 100)
@@ -248,7 +248,7 @@ def QuantumSignalProcessingPhases(
                 "The angle finding program failed on given instance, with an error of {}. Please relax the error budget and / or the success probability.".format(max_err))
         """
 
-        # For now, return reduced and full phases
+        # For now, return reduced, full phases, and parity together.
         return (qsp_seq_opt.full_phases, phases, parity)
 
     else:

pyqsp/main.py

@@ -10,14 +10,11 @@
 from pyqsp import angle_sequence, response
 from pyqsp.phases import phase_generators
 
-# TODO: create a function within pyqsp.poly that computes interpolants of desired functions up to a given degree; these are rescaled by a computed value if they need to be (with an announced flag) or otherwise left alone and fed to later methods. We also need to check these for parity, which can be done at compute time, again cleaning things up-stream.
-
 from pyqsp.poly import (StringPolynomial, TargetPolynomial,
                         polynomial_generators, PolyTaylorSeries)
 
 # -----------------------------------------------------------------------------
 
-
 class VAction(argparse.Action):
     def __call__(self, parser, args, values, option_string=None):
         curval = getattr(args, self.dest, 0) or 0
@@ -44,6 +41,7 @@ def CommandLine(args=None, arglist=None):
     angles      - generate QSP phase angles for the specified --seqname and --seqargs
     poly        - generate QSP phase angles for the specified --polyname and --polyargs, e.g. sign and threshold polynomials
     polyfunc    - generate QSP phase angles for the specified --func and --polydeg using tensorflow + keras optimization method (--tf)
+    sym_qsp_func    - generate QSP phase angles for the specified --func and --polydeg using symmetric QSP iterative method. Note that the desired polynomial is automatically scaled, and appears as the imaginary part of the top-left unitary matrix element in the standard basis.
     response    - generate QSP polynomial response functions for the QSP phase angles specified by --phiset
 
 Examples:
@@ -512,23 +510,25 @@ def float_list(value):
         if args.plot:
             response.PlotQSPResponse(phiset, signal_operator="Wx", **plot_args)
 
-    # TODO: Given some polymorphism in the options below, we should be able to call something like `pyqsp --plot --func "np.cos(3*x)" --polydeg 6 sym_qsp_func` and have this thing spit out a plot and a series of phases.
-    #
-    # New option added to support general Chebyshev interpolation
+    # New --sym_qsp_func option added to support general Chebyshev interpolation, in parallel with usage of --polyfunc
     elif args.cmd == "sym_qsp_func":
         if (not args.func) or (not args.polydeg):
             print(f"Must specify --func and --polydeg")
             return
 
-        # TODO: Currently test with hardcoded polynomial, as the use of StringPolynomial is restricted only to the tensorflow case; we will need a secondary method, most likely in poly.py, which returns the proper interpolating polynomial, given an evaluatable expression.
+        # Try to cast string as function; if no immediate error thrown at 0.5, generate anonymous function for later use.
+        try:
+            ret = eval(args.func, globals(), {'x': 0.5})
+        except Exception as err:
+            raise ValueError(
+                f"Invalid function specifciation, failed to evaluate at x=0.5, err={err}")
 
-        # poly = Chebyshev([0, 0, 0, 0, 0, 1])
+        # Generate anonymous function evaluated at x.
+        func = lambda x: eval(args.func, globals(), {'x': x})
 
-        # NOTE TEMPORARY HARDCODED FUNC HERE AND IN CALL BELOW.
-        func = lambda x: np.cos(3*x)
+        # Note that the polynomial is renormalized once according to its maximum magnitude (both up and down), and then again according to specified scale in a multiplicative way. We ought to be able to specify this, or at least return the metadata.
 
-        # TODO: Note that the polynomial is renormalized once according to its maximum magnitude, and then again according to specified scale in a multiplicative way.
-        # TODO: note that QuantumSignalProcessingPhases will determine if the parity is not definite, allowing us to use this method raw.
+        # Note that QuantumSignalProcessingPhases will determine if the parity is not definite, allowing us to use below method raw.
         base_poly = PolyTaylorSeries()
         poly, scale = base_poly.taylor_series(
             func=func,
@@ -537,7 +537,7 @@ def float_list(value):
             return_scale=True,
             max_scale=0.9,
             chebyshev_basis=True,
-            cheb_samples=50)
+            cheb_samples=4*args.polydeg) # Set larger to prevent aliasing
 
         # Modify choice of method globally.
         qspp_args['method'] = 'sym_qsp'
@@ -548,12 +548,13 @@ def float_list(value):
             chebyshev_basis=True,
             **qspp_args)
 
-        # TODO: determine which measurement, signal operator, re and im parts, and target functions are being used so that everything matches in the Chebyshev picture.
+        # Finally, if plotting called for, plot while passing 'sym_qsp' flag.
         if args.plot:
             response.PlotQSPResponse(
                 phiset,
                 target=poly,
                 signal_operator="Wx",
+                sym_qsp=True,
                 **plot_args)
 
     elif args.cmd == "polyfunc":

pyqsp/response.py

@@ -9,7 +9,8 @@ def ComputeQSPResponse(
         adat,
         phiset,
         signal_operator="Wx",
-        measurement=None):
+        measurement=None,
+        sym_qsp=False):
     """
     Compute QSP response.
 
@@ -69,12 +70,21 @@ def qsp_op(phi): return H @ np.array(
     for phi in phiset:
         pmats.append(qsp_op(phi))
 
-    for a in adat:
-        W = sig_op(a)
-        U = pmats[0]
-        for pm in pmats[1:]:
-            U = U @ W @ pm
-        pdat.append((p_state.T @ U @ p_state)[0, 0])
+    if not sym_qsp:
+        for a in adat:
+            W = sig_op(a)
+            U = pmats[0]
+            for pm in pmats[1:]:
+                U = U @ W @ pm
+            pdat.append((p_state.T @ U @ p_state)[0, 0])
+    else:
+        # For symmetric QSP, we plot matrix element directly, rather than the amplitude according to a measurement.
+        for a in adat:
+            W = sig_op(a)
+            U = pmats[0]
+            for pm in pmats[1:]:
+                U = U @ W @ pm
+            pdat.append(U[0, 0])
 
     pdat = np.array(pdat, dtype=np.complex128)
 
@@ -99,7 +109,8 @@ def PlotQSPResponse(
         plot_positive_only=False,
         plot_real_only=False,
         plot_tight_y=False,
-        show_qsp_model_plot=False):
+        show_qsp_model_plot=False,
+        sym_qsp=False):
     """
     Plot QSP response.
 
@@ -140,7 +151,8 @@ def PlotQSPResponse(
     qspr = ComputeQSPResponse(adat,
                               phiset,
                               signal_operator=signal_operator,
-                              measurement=measurement)
+                              measurement=measurement,
+                              sym_qsp=sym_qsp)
     pdat = qspr['pdat']
     plt.figure(figsize=[8, 5])
 
@@ -174,7 +186,13 @@ def PlotQSPResponse(
         ymin = np.min(np.real(pdat))
 
     # format plot
-    plt.ylabel("response")
+
+    # Modify labels depending on plotting matrix element or amplitude
+    if not sym_qsp:
+        plt.ylabel("Response")
+    else:
+        plt.ylabel("Matrix element")
+
     plt.xlabel("a")
     plt.legend(loc="upper right")
 

pyqsp/sym_qsp_plotting.py

@@ -15,13 +15,13 @@ def main():
     # Call existing methods to compute Jacobi-Anger expression for cosine.
     freq = 16
     pg = poly.PolyCosineTX()
-    
+
     # Note new argument to use Chebyshev basis.
     pcoefs = pg.generate(tau=freq, epsilon=1e-12, chebyshev_basis=True)
-    
+
     # Generate anonymous function (approx to cosine) using pcoefs.
     cos_fun = lambda x: np.polynomial.chebyshev.chebval(x, pcoefs)
-    
+
     # Initialize definite parity coefficients, and slice out nontrivial ones.
     parity = 0
     coef = pcoefs[parity::2]
@@ -30,7 +30,7 @@ def main():
     true_fun = lambda x: 0.5*np.cos(freq*x)
 
     # Optimize to the desired function using Newton solver.
-    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_Solver(coef, parity, crit=1e-12)    
+    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_solver(coef, parity, crit=1e-12)
     print("phases: %s\nerror: %s\niter: %s\n"%(str(phases), str(err), str(total_iter)))
 
     # Plot achieved versus desired function over samples.
@@ -90,13 +90,13 @@ def main():
     # Call existing methods to compute Jacobi-Anger expression for sine.
     freq = 16
     pg = poly.PolySineTX()
-    
+
     # Note new argument to use Chebyshev basis.
     pcoefs = pg.generate(tau=freq, epsilon=1e-12, chebyshev_basis=True)
-    
+
     # Generate anonymous function (approx to cosine) using pcoefs.
     sin_fun = lambda x: np.polynomial.chebyshev.chebval(x, pcoefs)
-    
+
     # Initialize definite parity coefficients, and slice out nontrivial ones.
     full_coef = pcoefs
     parity = 1
@@ -106,7 +106,7 @@ def main():
     true_fun = lambda x: 0.5*np.sin(freq*x)
 
     # Optimize to the desired function using Newton solver.
-    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_Solver(coef, parity, crit=1e-12)    
+    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_solver(coef, parity, crit=1e-12)
     print("phases: %s\nerror: %s\niter: %s\n"%(str(phases), str(err), str(total_iter)))
 
     # Plot achieved versus desired function over samples.
@@ -165,24 +165,24 @@ def main():
 
     # Generate inverse polynomial approximation
     pg = poly.PolyOneOverX()
-    
+
     # Underlying parameters of inverse approximation.
     # We use return_scale=True for ease of plotting correct desired function.
     kappa=5
     epsilon=0.01
     pcoefs, scale = pg.generate(kappa=kappa, epsilon=epsilon, chebyshev_basis=True, return_scale=True)
-    
+
     # Generate anonymous approximation and ideal function for scaled reciprocal.
     inv_fun = lambda x: np.polynomial.chebyshev.chebval(x, pcoefs)
     ideal_fun = lambda x: scale*np.reciprocal(x)
-    
+
     # Using odd parity and instantiating desired coefficeints.
     parity = 1
     coef = pcoefs[parity::2]
 
     # Optimize to the desired function using Newton solver.
     crit = 1e-12
-    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_Solver(coef, parity, crit=crit)    
+    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_solver(coef, parity, crit=crit)
     print("phase len: %s\nerror: %s\niter: %s\n"%(str(len(phases)), str(err), str(total_iter)))
 
     # Generate samples for plotting.
@@ -194,7 +194,7 @@ def main():
 
     # Grab im part of QSP unitary top-left matrix element.
     im_vals = np.array(qsp_seq_opt.gen_response_im(samples))
-    
+
     # Generate plotted values.
     approx_vals = np.array(list(map(inv_fun, samples)))
     # NOTE: For some reason this map casts to have an additional dimension.
@@ -211,11 +211,11 @@ def main():
     # Plot difference between two on log-plot
     diff = np.abs(im_vals - approx_vals)
     approx_diff = np.abs(ideal_vals - approx_vals)
-    
+
     # Plot QSP output polynomial versus desired polynomial, and error bound.
     axs[1].plot(samples, diff, 'r', label="Approx vs QSP")
     axs[1].plot(samples, [crit]*len(samples), 'y', label="QSP error limit")
-    
+
     # Plot approximation versus ideal function, and error bound.
     axs[1].plot(samples, approx_diff, 'g', label="Approx vs ideal")
     axs[1].plot(samples, [epsilon]*len(samples), 'b', label="Approx error limit")
@@ -251,20 +251,20 @@ def main():
     # Call existing methods to compute approximation to rect.
     freq = 16
     pg = poly.PolySign()
-    
+
     # TODO: note that definition of PolySign has been changed to return bare pcoefs and not TargetPolynomial
     pcoefs, scale = pg.generate(
             degree=161,
             delta=25,
-            chebyshev_basis=True, 
+            chebyshev_basis=True,
             cheb_samples=250,
             return_scale=True)
     # Cast from TargetPolynomial class bare Chebyshev coefficients if not using return_scale.
     # pcoefs = pcoefs.coef
-    
+
     # Generate anonymous function (approx to cosine) using pcoefs.
     sign_fun = lambda x: np.polynomial.chebyshev.chebval(x, pcoefs)
-    
+
     # Initialize definite parity coefficients, and slice out nontrivial ones.
     parity = 1
     bare_coefs = pcoefs
@@ -274,7 +274,7 @@ def main():
     true_fun = lambda x: scale * scipy.special.erf(x * 20)
 
     # Optimize to the desired function using Newton solver.
-    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_Solver(coef, parity, crit=1e-12)    
+    (phases, err, total_iter, qsp_seq_opt) = sym_qsp_opt.newton_solver(coef, parity, crit=1e-12)
     print("phases: %s\nerror: %s\niter: %s\n"%(str(phases), str(err), str(total_iter)))
 
     # Plot achieved versus desired function over samples.
@@ -324,4 +324,4 @@ def main():
     plt.show()
 
 if __name__ == '__main__':
-    main()
+    main()

pyqsp/test/test_main.py

@@ -18,6 +18,7 @@
     "--plot-positive-only --plot-real-only --seqargs 20,0.3 efilter",
     "--plot-positive-only --plot-real-only --seqargs=20,3.5 gibbs",
     "--plot-real-only --seqargs=20,0.6,15 relu",
+    "--return-angles --func np.cos(3*x) --polydeg 41 sym_qsp_func",
 ]
 
 """