Updated README figures and examples, and fixed plotting bugs for rect…

…angle functions.

README.md

@@ -1,4 +1,5 @@
-# :hammer_and_wrench: `pyQSP`: a Python Package for Quantum Signal Processing
+# :hammer_and_wrench: `pyQSP`
+## A Python Package for Quantum Signal Processing
 
 ![test workflow](https://github.com/ichuang/pyqsp/actions/workflows/run_tests.yml/badge.svg)
 
@@ -124,26 +125,24 @@ The method specified above works well up to above `polydeg=100`, and is limited
 
 ### A mélange of methods for phase finding and plotting
 
-This package includes various tools for plotting aspects of the computed QSP unitaries, many of which can be run from the command line. As an example, in the chart below the dashed line shows the target ideal polynomial QSP *response function* approximating a scaled version of $1/a$ over a sub-interval of $[-1,1]$. The dark line shows the real part of an actual response function, i.e., the matrix element $P_x(a)$, achieved by a QSP circuit with computed phases, while the blue line shows the imaginary part of the QSP response, with `n = 20` (the length of the QSP phase list less one).
+This package includes various tools for plotting aspects of the computed QSP unitaries, many of which can be run from the command line. As an example, in the chart below the dashed line shows the target ideal polynomial QSP *response function* approximating a scaled version of $1/a$ over a sub-interval of $[-1,1]$. The blue line shows the imaginary part of an actual response function, i.e., the matrix element $P_x(a)$, achieved by a QSP circuit with computed phases, while the red line in the bottom plot shows the error between achieved and desired functions on a logarithmic plot. Here `n = 147` (the length of the QSP phase list less one).
 
 <p align="center">
-    <img src="https://github.com/ichuang/pyqsp/blob/master/docs/ex_inversion.png" alt="QSP response function for the inverse function 1/a" width="75%"/>
+    <img src="https://github.com/ichuang/pyqsp/blob/master/docs/ex_qsp_inversion.png" alt="QSP response function for the inverse function 1/a" width="75%"/>
 </p>
 
-This was generated by running `pyqsp --plot --tolerance=0.01 --seqargs 3 invert`, which also spits out the the following verbose text:
+This was generated by running `pyqsp --plot --tolerance=0.01 --seqargs 6 --method sym_qsp invert`, which also spits out the the following verbose text:
 
 ```python
-    b=30, j0=14
-    [PolyOneOverX] minimum [-3.5325637] is at [-0.20530335]: normalizing
-    [PolyOneOverX] bounding to 0.5
-    [pyqsp.PolyOneOverX] pcoefs=[ 0.00000000e+00  4.24568213e+00  0.00000000e+00 -6.14813187e+01
-      0.00000000e+00  5.70160728e+02  0.00000000e+00 -3.77110116e+03
-      0.00000000e+00  1.86774294e+04  0.00000000e+00 -7.07245446e+04
-      0.00000000e+00  2.06037512e+05  0.00000000e+00 -4.61025085e+05
-      0.00000000e+00  7.86778785e+05  0.00000000e+00 -1.01130011e+06
-      0.00000000e+00  9.59305607e+05  0.00000000e+00 -6.49556764e+05
-      0.00000000e+00  2.96436030e+05  0.00000000e+00 -8.15921088e+04
-      0.00000000e+00  1.02215704e+04]
+b=147, j0=35
+[PolyOneOverX] minimum [-7.75401864] is at [-0.09251427]: normalizing
+[pyqsp.PolyOneOverX] pcoefs=[
+    0.00000000e+00  2.21344679e-01  0.00000000e+00 -1.99904501e-01
+    0.00000000e+00  1.78896005e-01  0.00000000e+00 -1.58587792e-01
+    0.00000000e+00  1.39221019e-01  0.00000000e+00 -1.21000963e-01
+    ...
+    0.00000000e+00  6.32566663e-85  0.00000000e+00 -4.30307535e-87
+    0.00000000e+00  1.45866961e-89]
 ```
 
 The sign function is also often useful to implement, e.g. for oblivious amplitude amplification. Running instead

pyqsp/main.py

@@ -269,7 +269,11 @@ def float_list(value):
             sys.exit(0)
         if isinstance(coefs, str):
             coefs = list(map(float, coefs.split(",")))
-        print(f"[pyqsp] polynomial coefficients={coefs}")
+
+        basis = "Chebyshev" if is_sym_qsp else "Monomial"
+        print(f"[pyqsp] {basis} polynomial coefficients={coefs}")
+        print(f"[CHECK] Coefficients expected in {basis} basis.")
+
         if is_sym_qsp:
             (phiset, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(
             coefs, **qspp_args)
@@ -508,9 +512,12 @@ def float_list(value):
         if args.plot:
             response.PlotQSPResponse(
                 phiset,
+                # target=lambda x: scale *
+                # (1 - (np.sign(x + 1 / args.seqargs[2]) -
+                #      np.sign(x - 1 / args.seqargs[2])) / 2),
                 target=lambda x: scale *
-                (1 - (np.sign(x + 1 / args.seqargs[2]) -
-                     np.sign(x - 1 / args.seqargs[2])) / 2),
+                (1 - (np.sign(x + 3 / (4 * args.seqargs[2])) -
+                     np.sign(x - 3 / (4 * args.seqargs[2]))) / 2),
                 signal_operator="Wx",
                 title="Rect Function",
                 **plot_args)
@@ -564,11 +571,12 @@ def float_list(value):
             print(
                 f'Known polynomial generators: {",".join(polynomial_generators.keys())}')
             return
+        # If one uses the poly argument, rather than a specific generator, then the name fed specifies the generator, and arguments fed are not seqargs but polyargs (which are in many cases the same format).
         pg = polynomial_generators[args.polyname]()
         if not args.polyargs:
             print(pg.help())
             return
-        pcoefs = pg.generate(*args.polyargs, chebyshev_basis=is_sym_qsp)
+        pcoefs, scale = pg.generate(*args.polyargs, return_scale=True, chebyshev_basis=is_sym_qsp)
         print(f"[pyqsp] polynomial coefs = {pcoefs}")
         if is_sym_qsp:
             (phiset, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(

pyqsp/poly.py

@@ -263,7 +263,7 @@ def generate(
                 pcoefs = g.coef
             else:
                 pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef)
-            print(f"[pyqsp.PolyOneOverX] pcoefs={pcoefs}")
+            # print(f"[pyqsp.PolyOneOverX] pcoefs={pcoefs}")
             if ensure_bounded and return_scale:
                 return pcoefs, scale
             else:
@@ -483,7 +483,7 @@ def erf_delta(x):
         pcoefs = the_poly.coef
         # force even coefficients to be zero, since the polynomial must be odd
         pcoefs[0::2] = 0
-        if ensure_bounded and return_scale:
+        if (ensure_bounded and return_scale) or chebyshev_basis:
             return pcoefs, scale
         else:
             return TargetPolynomial(pcoefs, target=lambda x: np.sign(x))
@@ -639,6 +639,10 @@ def rect(x):
             return 1 + (erf_delta(x - 3 / (4 * kappa)) +
                         erf_delta(-x - 3 / (4 * kappa))) / 2
 
+        # target=lambda x: scale *
+        # (1 - (np.sign(x + 1 / args.seqargs[2])
+        # - np.sign(x - 1 / args.seqargs[2])) / 2)
+
         if ensure_bounded and return_scale:
             the_poly, scale = self.taylor_series(
                 rect,
@@ -770,7 +774,7 @@ def gibbs(x):
         pcoefs = the_poly.coef
         # force odd coefficients to be zero, since the polynomial must be even
         pcoefs[1::2] = 0
-        if ensure_bounded and return_scale:
+        if (ensure_bounded and return_scale) or chebyshev_basis:
             return pcoefs, scale
         else:
             return TargetPolynomial(pcoefs, target=lambda x: gibbs(x))

pyqsp/sym_qsp_opt.py

@@ -246,7 +246,7 @@ def newton_solver(coef, parity, **kwargs):
     if 'maxiter' in kwargs:
         maxiter = kwargs['maxiter']
     else:
-        maxiter = 1e5
+        maxiter = 1e3
 
     # Currently deprecated, but real and imaginary parts can be alternately targeted by flipping overall sign of coef.
     # # targetPre = True
@@ -261,6 +261,8 @@ def newton_solver(coef, parity, **kwargs):
     qsp_seq_opt = SymmetricQSPProtocol(reduced_phases=reduced_phases, parity=parity)
     curr_iter = 0
 
+    print(f"[sym_qsp] Iterative optimization to err {crit:.3e} or max_iter {int(maxiter)}.")
+
     # Start the main loop
     while True:
         # Recover evaluated differences and Jacobian.
@@ -269,7 +271,8 @@ def newton_solver(coef, parity, **kwargs):
         err = np.linalg.norm(res, ord=1) # Take the one norm error.
         curr_iter = curr_iter + 1
 
-        print("iter: %s --- err: %.4g"%(str(curr_iter).zfill(3), err))
+        # Format string to show running error computation
+        print(f"iter: {curr_iter:03} --- err: {err:.3e}")
 
         # Invert Jacobian at evaluated point to determine direction of next step.
         lin_sol = np.linalg.solve(DFval, res)
@@ -278,10 +281,10 @@ def newton_solver(coef, parity, **kwargs):
 
         # Break conditions (beyond maxiter or within err.)
         if curr_iter >= maxiter:
-            print("Max iteration reached.\n")
+            print("[sym_qsp] Max iteration reached.")
             break
         if err < crit:
-            print("Stop criteria satisfied.\n")
+            print("[sym_qsp] Stop criteria satisfied.")
             break
 
     return (qsp_seq_opt.reduced_phases, err, curr_iter, qsp_seq_opt)