diff --git a/pygsti/adhoc/su2rbsims.py b/pygsti/adhoc/su2rbsims.py index 96b8c0561..540cd65f8 100644 --- a/pygsti/adhoc/su2rbsims.py +++ b/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 \ No newline at end of file + 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 diff --git a/pygsti/tools/su2tools.py b/pygsti/tools/su2tools.py index c457171f1..41c6fa8b0 100644 --- a/pygsti/tools/su2tools.py +++ b/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