added fitting and plotting helpers to su2rbsims. Made SU2.angles2irre…

…pchars use the more efficient method of computing characters. Removed an old, unused function for computing characters.

pygsti/adhoc/su2rbsims.py

@@ -10,6 +10,8 @@
 from typing import List, Tuple, Callable
 from tqdm import tqdm
 import functools
+from matplotlib import pyplot as plt
+from matplotlib import axes as plt_axes
 
 
 def get_M():
@@ -343,9 +345,9 @@ def __init__(self, su2rep, N, lengths, statepreps, povm, seed):
     def _compute_chars(self):
         # This is a separate function since it might take a while.
         for j in range(self.num_lens):
-            for k in range(self.N):
-                angles = self.all_circuits[j][k]
-                self.chars[j,k,:] = self.angles2irrepchars(angles[0,:])
+            circuits = np.array(self.all_circuits[j])
+            firstgate_angles = circuits[:,0,:]
+            self.chars[j,:,:] = self.angles2irrepchars(firstgate_angles)
         return
 
 
@@ -434,4 +436,95 @@ def mean_charweighted_empirical_distribution(distributions, wchars):
             survival_probs[ :, j] = P[diag_statepreps]
             synthetic_probs[:, j] = Q[diag_povmeffects]
 
-        return synthetic_probs, survival_probs
+        return synthetic_probs, survival_probs
+
+
+class Analysis:
+
+    plot_colors = default_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+
+    @staticmethod
+    def exp_decay_logvar(x,a,b):
+        return a * np.exp(-b * x)
+
+    @staticmethod
+    def exp_decay(x,a,b):
+        return a * b ** x
+
+    @staticmethod
+    def fit_exp(x, y, fitlog):
+        assert x.size == y.size
+        assert x.ndim == y.ndim == 1
+
+        if fitlog:
+            model = Analysis.exp_decay_logvar
+            p0 = np.array([1, 0.05])
+            bounds = ((-np.inf, -np.inf), (np.inf, np.inf))
+        else:
+            p0 = np.array([1, np.exp(-0.05)])
+            model = Analysis.exp_decay
+            bounds = ((-np.inf, 0), (np.inf, 1))
+
+        from scipy.optimize import curve_fit
+        p, pcov = curve_fit(model, x, y, p0=p0, bounds=bounds, maxfev=6000)
+        stddevs = np.sqrt(np.diag(pcov))
+        if fitlog:
+            stddevs[1] = np.exp(-p[1])*abs(2*np.sinh(stddevs[1]))
+        return p, stddevs
+
+    @staticmethod
+    def fit_exp_batch(x, ys, fitlog):
+        num_series = ys.shape[0]
+        ps = np.zeros((num_series, 2))
+        sigmas = np.zeros((num_series, 2))
+        for i,y in enumerate(ys):
+            p, sigma = Analysis.fit_exp(x, y, fitlog)
+            ps[i,:] = p
+            sigmas[i,:] = sigma
+        return ps, sigmas
+
+    @staticmethod
+    def fit_and_get_rates(x, ys, fitlog, F=None):
+        assert np.all(x > 0)
+        if F is None:
+            F = get_F()
+        num_series = ys.shape[0]
+        assert F.shape == (num_series, num_series)
+        ps, sigmas = Analysis.fit_exp_batch(x, ys, fitlog)
+        f_mu = ps[:,1]
+        rates = la.solve(F, f_mu)
+        return rates, ps, sigmas
+
+    @staticmethod
+    def fit_and_plot(x, pkm, ax : plt_axes.Axes, fitlog):
+        ps, sigmas = Analysis.fit_exp_batch(x, pkm, fitlog)
+        f_mu = ps[:,1]
+        f_sig = sigmas[:,1]
+        default_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+        model = Analysis.exp_decay_logvar if fitlog else Analysis.exp_decay
+        num_series = f_mu.size
+        for i in range(num_series):
+            y = pkm[i,:]
+            textstr = f'{i}: f = {f_mu[i]:.3f}±({f_sig[i]:.3f})'
+            ax.scatter(x, y, label=textstr, color=default_colors[i % len(default_colors)], s=20)
+            ax.plot(x, model(x, ps[i,0], ps[i,1]), color=default_colors[i % len(default_colors)], linestyle='-')
+        ax.legend()
+        return ps, sigmas
+
+    @staticmethod
+    def fit_and_plot_with_rates(x, ys, ax : plt_axes.Axes, fitlog, F=None):
+        if F is None:
+            F = get_F()
+        rates, ps, sigmas = Analysis.fit_and_get_rates(x, ys, fitlog)
+        f_mu = ps[:,1]
+        f_sig = sigmas[:,1]
+        rates = la.solve(F, f_mu)
+        default_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+        model = Analysis.exp_decay_logvar if fitlog else Analysis.exp_decay
+        for i,r in enumerate(rates):
+            y = ys[i,:]
+            textstr = f'{i}: f = {f_mu[i]:.3f}±({f_sig[i]:.3f}) | rate = {r:.3f}'
+            ax.scatter(x, y, label=textstr, color=default_colors[i % len(default_colors)], s=20)
+            ax.plot(x, model(x, ps[i,0], ps[i,1]), color=default_colors[i % len(default_colors)], linestyle='-')
+        ax.legend()
+        return rates, ps, sigmas

pygsti/tools/su2tools.py

@@ -251,10 +251,10 @@ def character_from_unitary(cls, U, j=1/2):
     @classmethod
     def angles2irrepchars(cls, angles):
         angles = np.atleast_1d(angles)
-        assert angles.ndim == 1
-        assert angles.size == 3
-        U = SU2.unitaries_from_angles(angles[0], angles[1], angles[2])[0]
-        out = np.array([ SU2.character_from_unitary(U, j) for j in cls.irrep_labels ])
+        angles = np.atleast_2d(angles).T
+        assert angles.shape == (3, angles.size // 3)
+        out = cls.characters_from_euler_angles(angles).T
+        assert out.shape == (angles.size // 3, cls.irrep_labels.size)
         return out
 
     @classmethod
@@ -266,14 +266,6 @@ def char(_k):
         out = np.array([ char(_k) for _k in cls.irrep_labels])
         return out
 
-    @staticmethod
-    def characters_from_angles(alphas, betas, gammas, j):
-        U2x2s = SU2.unitaries_from_angles(alphas, betas, gammas)
-        characters = np.zeros(len(U2x2s))
-        for i,U in enumerate(U2x2s):
-            characters[i] = SU2.character_from_unitary(U, j)
-        return characters
-
 
 def clebsh_gordan_matrix_spin72():
     """
@@ -505,27 +497,3 @@ def all_characters_from_unitary(cls, U):
             idx += b_sz
         out = np.array(out)
         return out
-
-
-def monte_carlo_projectors(N:int, seed=0) -> np.ndarray:
-    """
-    This function computes an array "Pis" of shape (8, 64, 64), where Pis[a] is
-    an N-sample Monte-Carlo estimate for the projector onto the a-th irrep of
-    the spin-7/2 superoperator representation of SU2. 
-
-    We compute five intermediate arrays on the way to getting "Pis."
-
-            su2s[b] = Euler angles for the b-th SU2 element
-            Us[b]   = 8-by-8 unitary representation of su2s[b]
-            Gs[b]   = 64-by-64 superoperator representation of su2s[b]
-
-            chars[a,b] = a-th irrep character for su2s[b]
-        wchars[a,b] = [dim of a-th irrep] / N * chars[a, b]
-    """
-    su2s = np.row_stack(SU2.random_euler_angles(N, True, seed))
-    Us = Spin72.unitaries_from_angles(*su2s)
-    Gs = [unitary_to_std_process_mx(U) for U in Us]
-    chars  = Spin72.new_characters_from_euler_angles(su2s)
-    wchars = Spin72.irrep_block_sizes[:,np.newaxis] * chars / N
-    Pis = np.tensordot(wchars, Gs, axes=1)
-    return Pis