From c53315b949aeabd386260b48a9d1e68105095afd Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Mon, 20 Feb 2023 14:03:55 +0100
Subject: [PATCH] small changes in plots

---
 final_plots.py            | 122 ++++++++++++-
 plot.py                   | 352 ++++++++++++++++++++++----------------
 plot_pyqtgraph.py         |   4 +
 tikz/commutator_RG2.tex   |   2 +-
 tikz/contour.tex          |   4 +-
 tikz/floquet_matrices.tex |   4 +-
 tikz/ga.tex               |   5 +-
 tikz/preamble.tex         |  31 +++-
 tikz/pulse_charge.tex     |  14 +-
 9 files changed, 371 insertions(+), 167 deletions(-)

diff --git a/final_plots.py b/final_plots.py
index 94533c7..ccb8c0c 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -226,7 +226,8 @@ def export_interp(
         y_func = None,
         kxorder = 2,
         kyorder = 2,
-        spline_s = 2e-4,
+        spline_s = 1e-4,
+        spline_s_ac = 2e-4,
         special_selection = None,
         o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
         o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
@@ -344,7 +345,7 @@ def export_interp(
             selected &= special_selection(thedata, order, method, padding, **kwargs)
         return thedata.loc[selected]
 
-    def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
+    def gen_data(order, method, padding=False, s=spline_s, s_ac=spline_s_ac, extend_vdc=10, extend_vac=10, **kwargs):
         suffix = "_p" if padding else ""
         reduced_data = selection(order, method, padding, **kwargs)
         results = {}
@@ -361,7 +362,7 @@ def export_interp(
             if extend_vac:
                 ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
                 ac_mirrored.vac *= -1
-                ac_mirrored.ac_current_phase += np.pi
+                ac_mirrored.ac_current_abs *= -1
                 extended_data.append(ac_mirrored)
                 if extend_vdc:
                     acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
@@ -441,7 +442,7 @@ def export_interp(
             except:
                 settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
             try:
-                iac_tck = bisplrep(*xy_data_i[::-1], 2*i_data.ac_current_abs, s=s, kx=kxorder, ky=kyorder)
+                iac_tck = bisplrep(*xy_data_i[::-1], 2*i_data.ac_current_abs, s=s_ac, kx=kxorder, ky=kyorder)
                 results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
                 results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
             except:
@@ -615,6 +616,78 @@ def export_interp_relative_o3a(
             selected &= special_selection(merged, order, method, padding, **kwargs)
         return merged.loc[selected]
 
+    def gen_data_cmp_init(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
+        suffix = "_p" if padding else ""
+        reduced_data = selection(order, method, padding, voltage_branches[order], **kwargs)
+        reduced_data = reduced_data[(reduced_data.solver_flags_ref & i_flags_good == i_flags_good) & (reduced_data.solver_flags & g_flags_bad == 0)]
+        results = {}
+        if reduced_data.shape[0] < 100:
+            return results
+        if extend_vdc or extend_vac:
+            extended_data = [reduced_data]
+            if extend_vdc:
+                dc_mirrored = reduced_data[reduced_data.vdc<extend_vdc].copy()
+                dc_mirrored.vdc *= -1
+                dc_mirrored.dc_current *= -1
+                dc_mirrored.dc_current_ref *= -1
+                extended_data.append(dc_mirrored)
+                del dc_mirrored
+            if extend_vac:
+                ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
+                ac_mirrored.vac *= -1
+                ac_mirrored.ac_current_abs *= -1
+                ac_mirrored.ac_current_abs_ref *= -1
+                extended_data.append(ac_mirrored)
+                if extend_vdc:
+                    acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
+                    acdc_mirrored.vdc *= -1
+                    acdc_mirrored.dc_current *= -1
+                    acdc_mirrored.dc_current_ref *= -1
+                    extended_data.append(acdc_mirrored)
+                    del acdc_mirrored
+                del ac_mirrored
+            reduced_data = pd.concat(extended_data)
+        settings.logger.info(f"Starting with init cmp. {order}{suffix} {method}, using {reduced_data.shape[0]} data points")
+        xy_data = (x_func(reduced_data), y_func(reduced_data))
+        try:
+            idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current - reduced_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
+            results[f"idc_{method}_{order}{suffix}_spline_cmp_init"] = bisplev(y_arr, x_arr, idc_tck)
+            results[f"gdc_{method}_{order}{suffix}_ispline_cmp_init"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
+        except:
+            settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method} (cmp init)")
+        try:
+            results[f"idc_{method}_{order}{suffix}_cmp_init"] = griddata(
+                    xy_data,
+                    reduced_data.dc_current - reduced_data.dc_current_ref,
+                    mesh,
+                    method="cubic")
+        except:
+            settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method} (cmp init)")
+        try:
+            results[f"iac_{method}_{order}{suffix}_cmp_init"] = griddata(
+                    xy_data,
+                    2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref),
+                    mesh,
+                    method="cubic")
+        except:
+            settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method} (cmp init)")
+        try:
+            iac_tck = bisplrep(*xy_data[::-1], 2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref), s=s, kx=kxorder, ky=kyorder)
+            results[f"iac_{method}_{order}{suffix}_spline_cmp_init"] = bisplev(y_arr, x_arr, iac_tck)
+            results[f"gac_{method}_{order}{suffix}_ispline_cmp_init"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
+        except:
+            settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method} (cmp init)")
+        try:
+            sel = reduced_data.vac > 0
+            results[f"phase_{method}_{order}{suffix}_cmp_init"] = griddata(
+                    (xy_data[0][sel], xy_data[1][sel]),
+                    reduced_data.ac_current_phase[sel] - reduced_data.ac_current_phase_ref[sel],
+                    mesh,
+                    method="cubic")
+        except:
+            settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method} (cmp init)")
+        return results
+
     def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
         suffix = "_p" if padding else ""
         reduced_data = selection(order, method, padding, voltage_branches[order], **kwargs)
@@ -646,12 +719,19 @@ def export_interp_relative_o3a(
                 del ac_mirrored
             reduced_data = pd.concat(extended_data)
         g_data = reduced_data.loc[(reduced_data.solver_flags_ref | reduced_data.solver_flags) & g_flags_bad == 0]
+        if order == "o3a_vb7":
+            g_data = reduced_data
         i_data = reduced_data.loc[reduced_data.solver_flags_ref & reduced_data.solver_flags & i_flags_good == i_flags_good]
         settings.logger.info(f"Starting with {order}{suffix} {method}, using {g_data.shape[0]}/{i_data.shape[0]} data points")
         xy_data = (x_func(g_data), y_func(g_data))
         has_current = i_data.shape[0] >= 100
         if has_current:
             xy_data_i = (x_func(i_data), y_func(i_data))
+        try:
+            idc_tck_old = bisplrep(*xy_data[::-1], g_data.dc_current - g_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
+            results[f"gdc_{method}_{order}{suffix}_ispline_old"] = bisplev(y_arr, x_arr, idc_tck_old, dy=1)
+        except:
+            settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method} (old)")
         try:
             results[f"gdc_{method}_{order}{suffix}"] = griddata(
                     xy_data,
@@ -660,6 +740,22 @@ def export_interp_relative_o3a(
                     method="cubic")
         except:
             settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
+        try:
+            results[f"idc_{method}_{order}{suffix}_old"] = griddata(
+                    xy_data,
+                    g_data.dc_current - g_data.dc_current_ref,
+                    mesh,
+                    method="cubic")
+        except:
+            settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method} (old)")
+        try:
+            results[f"iac_{method}_{order}{suffix}_old"] = griddata(
+                    xy_data,
+                    2*(g_data.ac_current_abs - g_data.ac_current_abs_ref),
+                    mesh,
+                    method="cubic")
+        except:
+            settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method} (old)")
         if has_current:
             try:
                 results[f"idc_{method}_{order}{suffix}"] = griddata(
@@ -711,6 +807,7 @@ def export_interp_relative_o3a(
     np.savez(
             filename,
             **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
+            **gen_data_cmp_init("o3a", "mu"),
             **gen_data("o3", "mu"),
             **gen_data("o3a_vb7", "mu"),
             **gen_data("o3p", "mu", integral_method=-1),
@@ -755,7 +852,8 @@ def export_omega5_interp(
                   y_arr = yarr,
                   kyorder = 3,
                   special_selection = special_selection,
-                  spline_s = 2e-5,
+                  spline_s = 5e-6,
+                  spline_s_ac = 1e-5,
                   d = 1e9,
                   solver_tol_rel = 1e-8,
                   solver_tol_abs = 1e-10,
@@ -797,7 +895,7 @@ def export_omega5_interp_deviation(
                   y_arr = yarr,
                   kyorder = 3,
                   special_selection = special_selection,
-                  spline_s = 2e-5,
+                  spline_s = 1e-5,
                   d = 1e9,
                   solver_tol_rel = 1e-8,
                   solver_tol_abs = 1e-10,
@@ -1423,6 +1521,7 @@ def export_pulse_current_pgfplots(omega=1.5, pulse_duration=0.01):
                 solver_tol_abs = 1e-10,
                 has_fourier_coef = True,
                 bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+                        | DataManager.SOLVER_FLAGS["improved_initial_conditions"] \
                         | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
                         | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                         | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
@@ -1483,6 +1582,7 @@ def export_pulse_charge_pgfplots(omega=1.5):
                 solver_tol_rel = 1e-8,
                 solver_tol_abs = 1e-10,
                 bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+                        | DataManager.SOLVER_FLAGS["improved_initial_conditions"] \
                         | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
                         | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                         | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
@@ -1584,13 +1684,13 @@ def export_harmonic_modes(
         n_phase_crossings = 10,
         ):
     label = "-omega20"
+    #label = "-omega1"
     vac_omega = np.array([5, 10, 20, 40, 80])
     dm = DataManager()
     data = dm.list(
                 omega = omega,
                 d = 1e9,
                 vdc = 0,
-                method = "J",
                 voltage_branches = 0,
                 include_Ga = True,
                 integral_method = -15,
@@ -1601,7 +1701,9 @@ def export_harmonic_modes(
                         | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
                         | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
                         | DataManager.SOLVER_FLAGS["deleted"],
+                good_flags = DataManager.SOLVER_FLAGS["improved_initial_conditions"]
                 )
+    data = data.loc[data.method != "mu"]
     current = []
     t = np.linspace(-0.5, 0.5, 201)
     current = []
@@ -1611,8 +1713,12 @@ def export_harmonic_modes(
     phase_crossing_currents = []
     for vac in vac_omega:
         sel = np.abs(data.vac-vac*omega) < 1e-6
-        if sel.sum() != 1:
+        if sel.sum() == 0:
+            settings.logger.warning(f"No results found for vac={vac}")
             continue
+        if sel.sum() > 1:
+            settings.logger.warning(f"Unexpected number of results for vac={vac}: {sel.sum()}")
+            sel *= data.nmax >= data[sel].nmax.max()
         row = data[sel].iloc[-1]
         filename = os.path.join(row.dirname, row.basename)
         kondo, = KondoImport.read_from_h5(filename, row.hash)
diff --git a/plot.py b/plot.py
index c6b6287..07529c8 100644
--- a/plot.py
+++ b/plot.py
@@ -912,9 +912,9 @@ def check_convergence(dm, **parameters):
     c_j = drtol[j]
     c_mu = drtol[mu]
     c_mod = drtol[mod]
-    s_j = 0.05 * table.nmax[j]**2
-    s_mu = 0.08 * table.nmax[mu]**2
-    s_mod = 0.08 * table.nmax[mod]**2
+    s_j = 0.12 * table.nmax[j]**2
+    s_mu = 0.15 * table.nmax[mu]**2
+    s_mod = 0.15 * table.nmax[mod]**2
     lw_j = 0.3 * table.voltage_branches[j]
     lw_mu = 0.3 * table.voltage_branches[mu]
     lw_mod = 0.3 * table.voltage_branches[mod]
@@ -1061,57 +1061,71 @@ def check_convergence_fit(dm, **parameters):
         settings.logger.warning(
                 "check_convergence expects specification of all physical parameters")
     table = dm.list(**parameters)
-    mod = (table.solver_flags & DataManager.SOLVER_FLAGS["simplified_initial_conditions"]) != 0
+    flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] | DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+    mod_bad = table.solver_flags & flags == DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
+    mod_old = table.solver_flags & flags == 0
+    mod_good = table.solver_flags & flags == DataManager.SOLVER_FLAGS["improved_initial_conditions"]
     d_fit_max = 9e11
     d_fit_max_j_nopadding = 9e11
     d_fit_max_j_shifted_nopadding = 2e9
-    j = (~mod) & (table.method == "J")
-    mu = (~mod) & (table.method == "mu")
+    j = mod_old & (table.method == "J")
+    mu_good = mod_good & (table.method == "mu")
+    mu_old = mod_old & (table.method == "mu")
+    mu_bad = mod_bad & (table.method == "mu")
     logd_inv_arr = np.linspace(0, 1/np.log10(table.d.min()), 200)
     logd_inv = 1/np.log10(table.d)
     x = 1/np.log10(table.d)
-    x3 = 1/np.log10(table.d)**3
     x_j = x[j]
-    x_mu = x[mu]
-    x_mod = x[mod]
+    x_good = x[mu_good]
+    x_old = x[mu_old]
+    x_bad = x[mu_bad]
     drtol = table.d*table.solver_tol_rel
     norm = LogNorm(max(drtol.min(), 0.05), min(drtol.max(), 500))
     c_j = drtol[j]
-    c_mu = drtol[mu]
-    c_mod = drtol[mod]
+    c_good = drtol[mu_good]
+    c_old = drtol[mu_old]
+    c_bad = drtol[mu_bad]
     s_j = 0.05 * table.nmax[j]**2
-    s_mu = 0.08 * table.nmax[mu]**2
-    s_mod = 0.08 * table.nmax[mod]**2
+    s_good = 0.08 * table.nmax[mu_good]**2
+    s_old = 0.08 * table.nmax[mu_old]**2
+    s_bad = 0.08 * table.nmax[mu_bad]**2
     lw_j = 0.3 * table.voltage_branches[j]
-    lw_mu = 0.3 * table.voltage_branches[mu]
-    lw_mod = 0.3 * table.voltage_branches[mod]
+    lw_good = 0.3 * table.voltage_branches[mu_good]
+    lw_old = 0.3 * table.voltage_branches[mu_old]
+    lw_bad = 0.3 * table.voltage_branches[mu_bad]
 
     selections = {}
     for r in range(10):
         suffix = '_%d'%r if r else ''
         if r in table.resonant_dc_shift.values:
-            mu_shift = (~mod) \
+            mu_good_shift = mod_good \
+                    & (table.method == "mu") \
+                    & (table.d <= d_fit_max) \
+                    & (table.resonant_dc_shift == r)
+            mu_old_shift = mod_old \
                     & (table.method == "mu") \
                     & (table.d <= d_fit_max) \
                     & (table.resonant_dc_shift == r)
-            mu_mod_shift = mod \
+            mu_bad_shift = mod_bad \
                     & (table.method == "mu") \
                     & (table.d <= d_fit_max) \
                     & (table.resonant_dc_shift == r)
-            j_shift = (~mod) \
+            j_shift = mod_old \
                     & (table.method == "J") \
                     & (     ((table.padding > 0) & (table.d <= d_fit_max)) \
                             | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding)) ) \
                     & (table.resonant_dc_shift == r)
-            j_mod_shift = mod \
+            j_mod_shift = mod_bad \
                     & (table.method == "J") \
                     & (     ((table.padding > 0) & (table.d <= d_fit_max)) \
                             | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding)) ) \
                     & (table.resonant_dc_shift == r)
-            if mu_shift.any():
-                selections["mu" + suffix] = mu_shift
-            if mu_mod_shift.any():
-                selections["mu_mod" + suffix] = mu_mod_shift
+            if mu_good_shift.any():
+                selections["mu_good" + suffix] = mu_good_shift
+            if mu_old_shift.any():
+                selections["mu_old" + suffix] = mu_old_shift
+            if mu_bad_shift.any():
+                selections["mu_bad" + suffix] = mu_bad_shift
             if j_shift.any():
                 selections["J" + suffix] = j_shift
             if j_mod_shift.any():
@@ -1119,17 +1133,23 @@ def check_convergence_fit(dm, **parameters):
 
     fit_func = lambda logd_inv, a, b, c: a + b * logd_inv**c
 
-    x_mean = x3[~mod].mean()
-    x_max = x3[~mod].max()
-    x_shifted = x3[~mod] - x_mean
-    s = (~mod).sum()
-    s_xx = (x_shifted**2).sum()
+    x_bad_mean = (x[mod_bad]**2).mean()
+    x_bad_max = (x[mod_bad]**2).max()
+    x_bad_shifted = x[mod_bad]**2 - x_bad_mean
+    s_bad = (mod_bad).sum()
+    s_bad_xx = (x_bad_shifted**2).sum()
 
-    xmod_mean = x3[mod].mean()
-    xmod_max = x3[mod].max()
-    xmod_shifted = x3[mod] - xmod_mean
-    smod = mod.sum()
-    smod_xx = (xmod_shifted**2).sum()
+    x_old_mean = (x[mod_old]**3).mean()
+    x_old_max = (x[mod_old]**3).max()
+    x_old_shifted = x[mod_old]**3 - x_old_mean
+    s_old = (mod_old).sum()
+    s_old_xx = (x_old_shifted**2).sum()
+
+    x_good_mean = (x[mod_good]**4).mean()
+    x_good_max = (x[mod_good]**4).max()
+    x_good_shifted = x[mod_good]**4 - x_good_mean
+    s_good = (mod_good).sum()
+    s_good_xx = (x_good_shifted**2).sum()
 
     # Create figure
     fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, sharey=False)
@@ -1188,171 +1208,215 @@ def check_convergence_fit(dm, **parameters):
     ax1.scatter(
             x_j,
             table.dc_current[j],
-            marker = 'x',
+            marker = 'o',
             s = s_j,
             c = c_j,
             linewidths = lw_j,
             norm = norm)
     ax1.scatter(
-            x_mu,
-            table.dc_current[mu],
+            x_good,
+            table.dc_current[mu_good],
+            marker = '*',
+            s = s_good,
+            c = c_good,
+            linewidths = lw_good,
+            norm = norm)
+    ax1.scatter(
+            x_bad,
+            table.dc_current[mu_bad],
             marker = '+',
-            s = s_mu,
-            c = c_mu,
-            linewidths = lw_mu,
+            s = s_bad,
+            c = c_bad,
+            linewidths = lw_bad,
             norm = norm)
     ax1.scatter(
-            x_mod,
-            table.dc_current[mod],
-            marker = '*',
-            s = s_mod,
-            c = c_mod,
-            linewidths = lw_mod,
+            x_old,
+            table.dc_current[mu_old],
+            marker = 'x',
+            s = s_old,
+            c = c_old,
+            linewidths = lw_old,
             norm = norm)
 
-    bb = table.dc_current[~mod].mean()
-    aa = (x_shifted*table.dc_current[~mod]).sum()/s_xx
-    #ax1.plot(logd_inv_arr, aa*(logd_inv_arr**3 - x_mean) + bb, linewidth=0.5)
-
-    bbmod = table.dc_current[mod].mean()
-    aamod = (xmod_shifted*table.dc_current[mod]).sum()/smod_xx
-    #ax1.plot(logd_inv_arr, aamod*(logd_inv_arr**3 - xmod_mean) + bbmod, linewidth=0.5)
-
-    for name, selection in selections.items():
-        try:
-            (a, b, c), covar = curve_fit(
-                    fit_func,
-                    logd_inv[selection],
-                    table.dc_current[selection],
-                    (bb-aa*x_mean, aa, 3))
-            aerr, berr, cerr = covar.diagonal()**0.5
-            print(f"DC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
-            ax1.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
-            ax1.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
-        except:
-            pass
+    bb_good = table.dc_current[mod_good].mean()
+    aa_good = (x_good_shifted*table.dc_current[mod_good]).sum()/s_good_xx
+    ax1.plot(logd_inv_arr, aa_good*(logd_inv_arr**4 - x_good_mean) + bb_good, linewidth=0.5)
+
+    bb_old = table.dc_current[mod_old].mean()
+    aa_old = (x_old_shifted*table.dc_current[mod_old]).sum()/s_old_xx
+    ax1.plot(logd_inv_arr, aa_old*(logd_inv_arr**3 - x_old_mean) + bb_old, linewidth=0.5)
+
+    bb_bad = table.dc_current[mod_bad].mean()
+    aa_bad = (x_bad_shifted*table.dc_current[mod_bad]).sum()/s_bad_xx
+    ax1.plot(logd_inv_arr, aa_bad*(logd_inv_arr**2 - x_bad_mean) + bb_bad, linewidth=0.5)
+
+    #for name, selection in selections.items():
+    #    try:
+    #        (a, b, c), covar = curve_fit(
+    #                fit_func,
+    #                logd_inv[selection],
+    #                table.dc_current[selection],
+    #                (bb-aa*x_mean, aa, 3))
+    #        aerr, berr, cerr = covar.diagonal()**0.5
+    #        print(f"DC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
+    #        ax1.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
+    #        ax1.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
+    #    except:
+    #        pass
 
     ax2.set_title("DC conductance")
     ax2.scatter(
             x_j,
             table.dc_conductance[j],
-            marker = 'x',
+            marker = 'o',
             s = s_j,
             c = c_j,
             linewidths = lw_j,
             norm = norm)
     ax2.scatter(
-            x_mu,
-            table.dc_conductance[mu],
+            x_good,
+            table.dc_conductance[mu_good],
+            marker = '*',
+            s = s_good,
+            c = c_good,
+            linewidths = lw_good,
+            norm = norm)
+    ax2.scatter(
+            x_old,
+            table.dc_conductance[mu_old],
             marker = '+',
-            s = s_mu,
-            c = c_mu,
-            linewidths = lw_mu,
+            s = s_old,
+            c = c_old,
+            linewidths = lw_old,
             norm = norm)
     ax2.scatter(
-            x_mod,
-            table.dc_conductance[mod],
-            marker = '*',
-            s = s_mod,
-            c = c_mod,
-            linewidths = lw_mod,
+            x_bad,
+            table.dc_conductance[mu_bad],
+            marker = 'x',
+            s = s_bad,
+            c = c_bad,
+            linewidths = lw_bad,
             norm = norm)
 
     ax3.set_title("AC current")
     ax3.scatter(
             x_j,
             table.ac_current_abs[j],
-            marker = 'x',
+            marker = 'o',
             s = s_j,
             c = c_j,
             linewidths = lw_j,
             norm = norm)
     ax3.scatter(
-            x_mu,
-            table.ac_current_abs[mu],
+            x_good,
+            table.ac_current_abs[mu_good],
+            marker = '*',
+            s = s_good,
+            c = c_good,
+            linewidths = lw_good,
+            norm = norm)
+    ax3.scatter(
+            x_old,
+            table.ac_current_abs[mu_old],
             marker = '+',
-            s = s_mu,
-            c = c_mu,
-            linewidths = lw_mu,
+            s = s_old,
+            c = c_old,
+            linewidths = lw_old,
             norm = norm)
     ax3.scatter(
-            x_mod,
-            table.ac_current_abs[mod],
-            marker = '*',
-            s = s_mod,
-            c = c_mod,
-            linewidths = lw_mod,
+            x_bad,
+            table.ac_current_abs[mu_bad],
+            marker = 'x',
+            s = s_bad,
+            c = c_bad,
+            linewidths = lw_bad,
             norm = norm)
 
-    b = table.ac_current_abs[~mod].mean()
-    a = (x_shifted*table.ac_current_abs[~mod]).sum()/s_xx
-    #ax3.plot(logd_inv_arr, a*(logd_inv_arr**3 - x_mean) + b, linewidth=0.5)
-
-    bmod = table.ac_current_abs[mod].mean()
-    amod = (xmod_shifted*table.ac_current_abs[mod]).sum()/smod_xx
-    #ax3.plot(logd_inv_arr, amod*(logd_inv_arr**3 - xmod_mean) + bmod, linewidth=0.5)
-
-    for name, selection in selections.items():
-        try:
-            (a, b, c), covar = curve_fit(
-                    fit_func,
-                    logd_inv[selection],
-                    table.ac_current_abs[selection],
-                    (bb-aa*x_mean, aa, 3))
-            aerr, berr, cerr = covar.diagonal()**0.5
-            print(f"AC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
-            ax3.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
-            ax3.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
-        except:
-            pass
+    b = table.ac_current_abs[mu_good].mean()
+    a = (x_good_shifted*table.ac_current_abs[mu_good]).sum()/s_good_xx
+    ax3.plot(logd_inv_arr, a*(logd_inv_arr**4 - x_good_mean) + b, linewidth=0.5)
+
+    b = table.ac_current_abs[mu_old].mean()
+    a = (x_old_shifted*table.ac_current_abs[mu_old]).sum()/s_old_xx
+    ax3.plot(logd_inv_arr, a*(logd_inv_arr**3 - x_old_mean) + b, linewidth=0.5)
+
+    b = table.ac_current_abs[mu_bad].mean()
+    a = (x_bad_shifted*table.ac_current_abs[mu_bad]).sum()/s_bad_xx
+    ax3.plot(logd_inv_arr, a*(logd_inv_arr**2 - x_bad_mean) + b, linewidth=0.5)
+
+    #for name, selection in selections.items():
+    #    try:
+    #        (a, b, c), covar = curve_fit(
+    #                fit_func,
+    #                logd_inv[selection],
+    #                table.ac_current_abs[selection],
+    #                (bb-aa*x_mean, aa, 3))
+    #        aerr, berr, cerr = covar.diagonal()**0.5
+    #        print(f"AC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
+    #        ax3.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
+    #        ax3.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
+    #    except:
+    #        pass
 
     ax4.set_title("AC phase")
     ax4.scatter(
             x_j,
             table.ac_current_phase[j],
-            marker='x',
+            marker='o',
             s=s_j,
             c=c_j,
             linewidths=lw_j,
             norm=norm)
     ax4.scatter(
-            x_mu,
-            table.ac_current_phase[mu],
+            x_good,
+            table.ac_current_phase[mu_good],
+            marker='*',
+            s=s_good,
+            c=c_good,
+            linewidths=lw_good,
+            norm=norm)
+    ax4.scatter(
+            x_old,
+            table.ac_current_phase[mu_old],
             marker='+',
-            s=s_mu,
-            c=c_mu,
-            linewidths=lw_mu,
+            s=s_old,
+            c=c_old,
+            linewidths=lw_old,
             norm=norm)
     ax4.scatter(
-            x_mod,
-            table.ac_current_phase[mod],
-            marker='*',
-            s=s_mod,
-            c=c_mod,
-            linewidths=lw_mod,
+            x_bad,
+            table.ac_current_phase[mu_bad],
+            marker='x',
+            s=s_bad,
+            c=c_bad,
+            linewidths=lw_bad,
             norm=norm)
 
-    b = table.ac_current_phase[~mod].mean()
-    a = (x_shifted*table.ac_current_phase[~mod]).sum()/s_xx
-    #ax4.plot(logd_inv_arr, a*(logd_inv_arr**3 - x_mean) + b, linewidth=0.5)
-
-    bmod = table.ac_current_phase[mod].mean()
-    amod = (xmod_shifted*table.ac_current_phase[mod]).sum()/smod_xx
-    #ax4.plot(logd_inv_arr, amod*(logd_inv_arr**3 - xmod_mean) + bmod, linewidth=0.5)
-
-    for name, selection in selections.items():
-        try:
-            (a, b, c), covar = curve_fit(
-                    fit_func,
-                    logd_inv[selection],
-                    table.ac_current_phase[selection],
-                    (bb-aa*x_mean, aa, 3))
-            aerr, berr, cerr = covar.diagonal()**0.5
-            print(f"AC phase ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
-            ax4.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
-            ax4.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
-        except:
-            pass
+    b = table.ac_current_phase[mu_good].mean()
+    a = (x_good_shifted*table.ac_current_phase[mu_good]).sum()/s_good_xx
+    ax4.plot(logd_inv_arr, a*(logd_inv_arr**4 - x_good_mean) + b, linewidth=0.5)
+
+    b = table.ac_current_phase[mu_old].mean()
+    a = (x_old_shifted*table.ac_current_phase[mu_old]).sum()/s_old_xx
+    ax4.plot(logd_inv_arr, a*(logd_inv_arr**3 - x_old_mean) + b, linewidth=0.5)
+
+    b = table.ac_current_phase[mu_bad].mean()
+    a = (x_bad_shifted*table.ac_current_phase[mu_bad]).sum()/s_bad_xx
+    ax4.plot(logd_inv_arr, a*(logd_inv_arr**2 - x_bad_mean) + b, linewidth=0.5)
+
+    #for name, selection in selections.items():
+    #    try:
+    #        (a, b, c), covar = curve_fit(
+    #                fit_func,
+    #                logd_inv[selection],
+    #                table.ac_current_phase[selection],
+    #                (bb-aa*x_mean, aa, 3))
+    #        aerr, berr, cerr = covar.diagonal()**0.5
+    #        print(f"AC phase ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
+    #        ax4.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
+    #        ax4.plot((0,logd_inv_arr[-1]), (a,a), ':', label=name)
+    #    except:
+    #        pass
 
 
 def convergence_overview(dm, **parameters):
diff --git a/plot_pyqtgraph.py b/plot_pyqtgraph.py
index 258c62f..9fd52da 100644
--- a/plot_pyqtgraph.py
+++ b/plot_pyqtgraph.py
@@ -803,6 +803,10 @@ def parse():
             help="Minimal major version")
     parser.add_argument("--min_version_minor", metavar="int", type=int, default=-1,
             help="Minimal minor version")
+    parser.add_argument("--good_flags", metavar="int", type=int, default=0,
+            help="Required solver flags")
+    parser.add_argument("--bad_flags", metavar="int", type=int, default=0,
+            help="Forbidden solver flags")
     args = parser.parse_args()
 
     options = args.__dict__
diff --git a/tikz/commutator_RG2.tex b/tikz/commutator_RG2.tex
index 1d7cbe9..7cfb798 100644
--- a/tikz/commutator_RG2.tex
+++ b/tikz/commutator_RG2.tex
@@ -33,7 +33,7 @@
       ytick = {0, 0.01, 0.02},
       yticklabels = {$0$, $0.01$, $0.02$},
       scaled ticks = false,
-      ylabel = {$\|[\chi^{-1},ZG^2_{LR}]\|_\infty\,/\,J$},
+      ylabel = {$\|[\hat\chi^{-1},\hat{Z}\Gtwo_{LR}]\|_\infty\,/\,J$},
     ]
     \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.0288735933] {../figdata/colorbar_viridis.png};
   \end{axis}
diff --git a/tikz/contour.tex b/tikz/contour.tex
index 7cdc29e..9e002cc 100644
--- a/tikz/contour.tex
+++ b/tikz/contour.tex
@@ -34,11 +34,11 @@
 
     \legend{
       %$O(J^2)$,
-      %{$O(J^3)$, $G^a=0$},
+      %{$O(J^3)$, $\Ga=0$},
       %{$O(J^3)$, approx.~$I[X]$},
       %$O(J^3)$,
       leading order,
-      without $G^a$,
+      without $\Ga$,
       approx.~integral,
       full RG equations
     }
diff --git a/tikz/floquet_matrices.tex b/tikz/floquet_matrices.tex
index b6b1c3a..603ee16 100644
--- a/tikz/floquet_matrices.tex
+++ b/tikz/floquet_matrices.tex
@@ -70,7 +70,7 @@
     \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_z_mu_diff_p.png};
   \end{axis}
   %
-  \node[xshift=2mm, yshift=8mm] at (z_mu_small_p.north east) {$Z$};
+  \node[xshift=2mm, yshift=8mm] at (z_mu_small_p.north east) {$\hat{Z}$};
   %
   \begin{axis}[
       name = z_J_small,
@@ -217,7 +217,7 @@
     \addplot graphics[xmin=-9.5, xmax=9.5, ymin=-9.5, ymax=9.5] {../figdata/floquet_matrix_g201_mu_diff_p.png};
   \end{axis}
   %
-  \node[xshift=2mm, yshift=8mm] at (g201_mu_small_p.north east) {$G^2_{LR}$};
+  \node[xshift=2mm, yshift=8mm] at (g201_mu_small_p.north east) {$\Gtwo_{LR}$};
   %
   \begin{axis}[
       name = g201_J_small,
diff --git a/tikz/ga.tex b/tikz/ga.tex
index 99d19af..e90d31c 100644
--- a/tikz/ga.tex
+++ b/tikz/ga.tex
@@ -7,7 +7,7 @@
       width = 4cm,
       height = 4cm,
       scale only axis,
-      title = {absolute value $\displaystyle\|G^a\|_\infty$},
+      title = {absolute value $\displaystyle\|\Ga\|_\infty$},
       xmin = 0, xmax = 10,
       ymin = 0, ymax = 10,
       xlabel = {$\vdc~(\Omega)$},
@@ -48,7 +48,8 @@
       xshift = 28mm,
       width = 4cm,
       height = 4cm,
-      title = {relative to $J^2$: $\displaystyle\frac{\|G^a\|_\infty}{\|G^2\|^2_\infty}$},
+      title = {relative to $J^2$: $\displaystyle\frac{\|\Ga\|_\infty}{\|\Gtwo\|^2_\infty}$},
+      title style = {yshift = -2.72mm},
       scale only axis,
       xmin = 0, xmax = 10,
       ymin = 0, ymax = 10,
diff --git a/tikz/preamble.tex b/tikz/preamble.tex
index e3508d4..8959bd8 100644
--- a/tikz/preamble.tex
+++ b/tikz/preamble.tex
@@ -1,12 +1,37 @@
 \newcommand\tikzsetnextfilename[1]{}
-\newcommand\tkv{T_\mathrm{K}}
-\newcommand\tkrg{T_\mathrm{K}'}
+\newcommand\tkv{T_K}
+\newcommand\tkrg{T_K'}
 \newcommand\vdc{V_\mathrm{avg}}
 \newcommand\vac{V_\mathrm{osc}}
-\newcommand\gdc{G_\mathrm{avg}}
+\newcommand\gdc{G}
 \newcommand\gac{G_\mathrm{osc}}
 \newcommand\idc{I_\mathrm{avg}}
 \newcommand\iac{I_\mathrm{osc}}
+
+% Coupling vertex parametrization
+\newcommand\Ga{\hat{G}^a}
+\newcommand\Gb{\hat{G}^b}
+\newcommand\Gone{\hat{G}^1}
+\newcommand\Gtwo{\hat{G}^2}
+\newcommand\Gthree{\hat{G}^3}
+\newcommand\Gna{G^a}
+\newcommand\Gnb{G^b}
+\newcommand\Gnone{G^1}
+\newcommand\Gntwo{G^2}
+\newcommand\Gnthree{G^3}
+
+% Superoperator parametrization
+\newcommand\Lbar{\bar{L}}
+\newcommand\Labar{\Lbar^a}
+\newcommand\Lbbar{\Lbar^b}
+\newcommand\Lalgebra{L}
+\newcommand\La{\Lalgebra^a}
+\newcommand\Lb{\Lalgebra^b}
+\newcommand\Lone{\Lalgebra^1}
+\newcommand\Ltwo{\Lalgebra^2}
+\newcommand\Lthree{\Lalgebra^3}
+
+
 \usepackage{pgfplots}
 \pgfplotsset{compat=1.18}
 \pgfplotsset{surf shading/precision=pdf}
diff --git a/tikz/pulse_charge.tex b/tikz/pulse_charge.tex
index 2147a79..3cb5c3b 100644
--- a/tikz/pulse_charge.tex
+++ b/tikz/pulse_charge.tex
@@ -13,13 +13,13 @@
       mark size = 1pt,
     ]
     \node[anchor=north west] at (axis description cs:0.02,0.97) {(a)};
-    \draw[blue] (0.25, -0.02) -- ++(0, 5pt);
+    \draw[blue] (0.25, -0.02) -- ++(0, 4.5pt);
     \draw[blue] (0.25, 0.83) -- ++(0, -7.5pt);
-    \draw[green] (0.5, -0.02) -- ++(0, 5pt);
+    \draw[green] (0.5, -0.02) -- ++(0, 4.5pt);
     \draw[green] (0.5, 0.83) -- ++(0, -7.5pt);
-    \draw[orange] (0.75, -0.02) -- ++(0, 5pt);
+    \draw[orange] (0.75, -0.02) -- ++(0, 4.5pt);
     \draw[orange] (0.75, 0.83) -- ++(0, -7.5pt);
-    \draw[magenta] (1, -0.02) -- ++(0, 5pt);
+    \draw[magenta] (1, -0.02) -- ++(0, 4.5pt);
     \draw[magenta] (1, 0.83) -- ++(0, -7.5pt);
     %
     \addplot[no marks, cyan] table [x=phase, y=q] {../figdata/pulse_charge_phase_interp.dat};
@@ -42,7 +42,9 @@
       mark size = 1pt,
     ]
     \node[anchor=north west] at (axis description cs:0.02,0.97) {(b)};
-    \draw[cyan] (0.137701404, -0.02) -- ++(0, 5pt);
+    \draw[->,black] (0.034425351, -0.02) -- ++(0, 4.5pt);
+    \draw[->,black] (0.034425351, 0.83) -- ++(0, -7.5pt);
+    \draw[cyan] (0.137701404, -0.02) -- ++(0, 4.5pt);
     \draw[cyan] (0.137701404, 0.83) -- ++(0, -7.5pt);
     %
     \addplot[no marks, blue] table [x=t, y=q1] {../figdata/pulse_charge_duration_interp.dat};
@@ -130,5 +132,7 @@
     %
     \draw[cyan] (0.137701404, 0) -- ++(0, 5pt);
     \draw[cyan] (0.137701404, 3.3) -- ++(0, -5pt);
+    \draw[->,black] (0.034425351, 0) -- ++(0, 5pt);
+    \draw[->,black] (0.034425351, 3.3) -- ++(0, -5pt);
   \end{axis}
 \end{tikzpicture}%
-- 
GitLab