update cutoff dependence

python/energy_flow_proper/08_dynamic_one_bath/cutoff_dep.py

@@ -15,6 +15,8 @@ from hops.util.logging_setup import logging_setup
 import logging
 logging_setup(logging.INFO, show_stocproc=False)
 
+print(es.models_table(ω_models))
+
 plt.plot(ωs, [model.full_thermal_bcf(0).real * model.bcf_scale for model in ω_models])
 # plt.plot(ωs, [model.bcf_scale for model in ω_models])
 
@@ -23,14 +25,18 @@ for model in ω_models:
 plt.legend()
 
 for model in ω_models:
-      vals = abs(model.full_thermal_bcf(model.t))
+      vals = abs(model.bcf(model.t))
-      plt.plot(model.t, vals / vals.max(), linewidth=1)
+      plt.plot(model.t, model.bcf_scale * vals, linewidth=1)
       plt.xlim(-.1, 6)
 
 ωs = np.linspace(.01, 30, 1000)
 for model in ω_models:
-      plt.plot(ωs, model.full_thermal_spectral_density(ωs), linewidth=1, label=fr"$\omega_c={model.ω_c:.2f}$")
+      plt.plot(ωs, model.bcf_scale * model.spectral_density(ωs), linewidth=1, label=fr"$\omega_c={model.ω_c:.2f}$")
 plt.legend()
+plt.xlabel(r"$\omega$")
+plt.ylabel(r"$J(\omega)$")
+
+fs.export_fig("omega_sd", tikz=False, x_scaling=.49, y_scaling=.3)
 
 bath_scales = [1/(model.bcf_coefficients()[1][0].real.min()) for model in ω_models]
 therm_scales = [Δ/(model.bcf_scale) for model in ω_models]
@@ -38,6 +44,7 @@ plt.plot(ωs, bath_scales)
 #plt.plot(ωs, therm_scales)
 
 fig, ax = plt.subplots()
+ax.axhline(0, color="gray")
 for model, data in aux.model_data_iterator(ω_models):
     _, _, bar = pu.plot_with_σ(
         model.t,
@@ -45,13 +52,20 @@ for model, data in aux.model_data_iterator(ω_models):
         ax=ax,
         label=fr"$\omega_c={model.ω_c:.2f}$",
     )
+ax.set_xlim(0, model.t[strobe_indices][3])
+ax.set_ylabel(r"$\langle H_\mathrm{I}\rangle$")
+ax.set_xlabel(r"$\tau$")
+ax.legend()
+#fs.export_fig("omega_interactions", fig, tikz=False, x_scaling=.49, y_scaling=.3)
 
 fig, ax = plt.subplots()
 import scipy
+
 for model, data in aux.model_data_iterator(ω_models):
+    energy = model.total_energy_from_power(data) * (1 / es.ergo(model.T))
     _, _, bar = pu.plot_with_σ(
         model.t,
-        model.total_energy_from_power(data) * (1/es.ergo(model.T)),
+        energy,
         ax=ax,
         label=fr"$\omega_c={model.ω_c:.2f}$",
         strobe_frequency=Δ,
@@ -60,7 +74,7 @@ for model, data in aux.model_data_iterator(ω_models):
 
     pu.plot_with_σ(
         model.t,
-        (model.total_energy_from_power(data) * (1/es.ergo(model.T))),
+        energy,
         ax=ax,
         # label=fr"$ω_c={model.ω_c:.2f}$",
         # strobe_frequency=Δ,
@@ -69,27 +83,108 @@ for model, data in aux.model_data_iterator(ω_models):
         color=bar.lines[0].get_color(),
     )
 
-    τ_off = 1 / (model.bcf_coefficients()[1][0].real.min())
+    # τ_off = 1 / (model.bcf_coefficients()[1][0].real.min())
-    a0 = abs(model.bcf(0))
+    # # a0 = abs(model.bcf(0))
-    τ_off = scipy.optimize.fsolve(lambda τ: abs(model.bcf(τ)) - (a0/300), 10)
+    # # τ_off = scipy.optimize.fsolve(lambda τ: abs(model.bcf(τ)) - (a0/300), 10)
+    # ax.axvline(
+    #     τ_off,
+    #     color=bar.lines[0].get_color(),
+    #     linestyle="dotted",
+    #     linewidth=1,
+    #     zorder=-1,
+    # )
+
+    ts = model.t[strobe_indices]
+    e_s = energy.value[strobe_indices]
+
+
+    ind_min = e_s.argmin()
+    t_min = ts[ind_min]
+    e_min = e_s[ind_min]
+
     ax.axvline(
-        τ_off,
+        t_min,
         color=bar.lines[0].get_color(),
+        ymin=0,
+        ymax=abs((e_min + 0.41) / (0.45 + 0.1)),
         linestyle="dotted",
-        linewidth=1,
-        zorder=-1,
     )
+    ax.axhline(
+        e_min,
+        color=bar.lines[0].get_color(),
+        xmin=0,
+        xmax=t_min / model.t.max(),
+        linestyle="dotted",
+    )
+
 
 ax.set_xlabel(r"$\tau$")
 ax.set_ylabel(r"$(\langle H\rangle_\tau -\langle H\rangle_0)/\mathcal{W}_\mathrm{max}$")
 # ax.axhline(, color="grey", linewidth=1, linestyle="dashed")
 # ax.plot(model.t, model.L.operator_norm(model.t), linewidth=1)
+ax.set_xlim(0, model.t.max())
+ax.set_ylim(-0.41, 0.1)
 ax.legend()
-#fs.export_fig("omegas_total")
+fs.export_fig("omegas_total", tikz=False, y_scaling=0.4)
+
+import scipy
+
+opt_powers = []
+opt_energies = []
+memories = []
+for model, data in aux.model_data_iterator(ω_models):
+    energy = model.total_energy_from_power(data)
+    ts = model.t[strobe_indices]
+    e_s = energy.value[strobe_indices]
+
+    powers = np.divide(e_s, ts, where=ts > 0)
+    p_min = powers.min()
+    ind_min = powers.argmin()
+    # t_min = ts[ind_min]
+    # e_min = e_s[ind_min]
+
+    opt_powers.append(abs(p_min))
+    opt_energies.append(abs(e_s.min()) * (1 / es.ergo(model.T)))
+
+    # bcf = model.full_thermal_bcf
+    # α0 = bcf(0).real
+    # obj = lambda τ: (abs(bcf(τ))/α0 - 1e-1) ** 2
+    # τ_off = scipy.optimize.basinhopping(
+    #     obj, 1
+    # ).x
+    memories.append(3 / model.ω_c)
+    # memories.append(τ_off)
+
+fig, ax = plt.subplots()
+lines = ax.plot(memories, opt_powers, marker="o", markersize=3)
+ax.set_xlabel(r"$\tau_\mathrm{B}=3 / \omega_c$")
+ax.set_ylabel(r"$P_\mathrm{max}$", color=lines[0].get_color())
+
+ax2 = ax.twinx()
+
+ax2.plot([], [])
+lines = ax2.plot(memories, opt_energies, marker="o", markersize=3)
+ax2.set_ylabel(
+    r"$|\Delta \langle H\rangle|_\mathrm{max} / \mathcal{W}_\mathrm{max}$",
+    color=lines[0].get_color(),
+)
+
+# def f(t, s, h):
+#     t = t-s
+#     out = np.zeros_like(t)
+#     np.log(t, out=out, where=t>1)
+
+#     return out * h
+
+# opt,_ = scipy.optimize.curve_fit(f, memories, opt_powers)
+# smooth_mem = np.linspace(memories[0], memories[-1], 1000)
+# ax.plot(smooth_mem, f(smooth_mem, *opt)) #
+plt.axvline(1 * np.pi * 2/Δ, color="gray", linestyle="dotted")
+fs.export_fig("omega_energies_and_powers", fig, tikz=False, x_scaling=.5, y_scaling=.35)
 
 fig, ax = plt.subplots()
 for model, data in aux.model_data_iterator(ω_models):
-    fs.plot_with_σ(
+    pu.plot_with_σ(
         model.t,
         model.system_energy(data),
         ax=ax,
@@ -105,16 +200,16 @@ ax.legend()
 
 fig, ax = plt.subplots()
 int_0 = ω_models[0].interaction_energy(aux.get_data(ω_models[0])).for_bath(0)
-for model, data in aux.model_data_iterator(ω_models):
+for model, data in aux.model_data_iterator(ω_models[1:]):
-    fs.plot_with_σ(
+    pu.plot_with_σ(
         model.t,
-        (model.interaction_energy(data) - int_0) * (1 / abs(int_0).max.value),
+        (abs(model.interaction_energy(data)) - abs(int_0)) * (1 / abs(int_0).max.value),
         ax=ax,
         bath=0,
-        label=fr"$ω_c={model.ω_c:.2f}$",
+        label=fr"$\omega_c={model.ω_c:.2f}$",
         linewidth=0.5,
     )
-#ax.legend()
+ax.legend()
 
 fig, ax = plt.subplots()
 fig.set_size_inches(fs.get_figsize("poster", 0.49))
@@ -213,9 +308,11 @@ ax.legend()
 ax.set_xlim((0,10))
 fs.export_fig("energy_shovel_preview")
 
+from hops.util.utilities import relative_entropy
+
 fig, ax = plt.subplots()
 for model,data in aux.model_data_iterator(ω_models):
-          fs.plot_with_σ(model.t, EnsembleValue(relative_entropy(data, strobe_indices[-1])), ax=ax, label=model.ω_c, strobe_frequency=Δ)
+          ax.plot(model.t, relative_entropy(data, strobe_indices[-1])[0], ax=ax, label=model.ω_c, strobe_frequency=Δ)
 ax.legend()
 
 fig, ax = plt.subplots()