milked 08:mod_freq, have to do more detailed simulations now

python/energy_flow_proper/08_dynamic_one_bath/mod_freq.py

@@ -9,7 +9,7 @@ import plot_utils as pu
 
 import ray
 ray.shutdown()
-ray.init(address="auto")
+ray.init()
 
 from hops.util.logging_setup import logging_setup
 import logging
@@ -17,39 +17,57 @@ logging_setup(logging.INFO, show_stocproc=False)
 
 from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix
 
-
+Δs = [
-Δs = np.linspace(1,10,10)#np.sort(np.concatenate((np.linspace(1, 5, 20), np.linspace(5, 7, int(20/5 * 2)))))
+    1,
+    2,
+    3,
+    4,
+    5,
+    6,
+    7,
+    8,
+    9,
+    10,
+]  # np.sort(np.concatenate((np.linspace(1, 5, 20), np.linspace(5, 7, int(20/5 * 2)))))
 Δ_models = []
 strobe_ts = []
 strobe_indices_s = []
-δs = np.linspace(.1, 2, 10)
+δs = np.linspace(0.1, 1.5, 10)
 for Δ in Δs:
     for δ in δs:
-        proto, strobe_t, strobe_indices = es.energy_shovel(Δ, periods=20, k_max=5, modulate_system=False)
+        proto, strobe_t, strobe_indices = es.energy_shovel(
+            Δ,
+            periods=int(Δ / (2 * np.pi) * 20),
+            k_max=3,
+            bcf_terms=4,
+            modulate_system=False,
+            bcf_norm_method="unit_therm",
+        )
         proto.δ = δ
         strobe_ts.append(strobe_t)
         strobe_indices_s.append(strobe_indices)
         Δ_models.append(proto)
 
-aux.integrate_multi(Δ_models, 10_000)
+Δ_models[1].bcf(0) * Δ_models[-1].bcf_scale
 
-final_e = [[], []]
+aux.integrate_multi(Δ_models, 2000)
-final_e_error = [[], []]
+
+dim_models = len(δs)
+power = np.zeros((len(Δs), dim_models))
+
+# final_e_error = [[], []]
 ensemble_arg = dict(overwrite_cache=False)
 fig, ax = plt.subplots()
-for (model, data), Δ, strobe_t, strobe_indices in zip(
+for (model, data), Δ, strobe_t, strobe_indices, idx in zip(
     aux.model_data_iterator(Δ_models),
     np.array([[Δ]*len(δs) for Δ in Δs]).flatten(),
     strobe_ts,
     strobe_indices_s,
+    range(len(Δ_models))
 ):
     energies = model.total_energy_from_power(data, **ensemble_arg)
-    idx = energies.value[strobe_indices].argmin()
+    powers = energies.value[strobe_indices[1:]] * (1/strobe_t[1:])
-    energies = energies * (1 / strobe_t[idx])
+    power[int(idx/dim_models), idx % dim_models] = abs(powers[powers <= 0]).max()
-
-    energy_idx = δs.index(model.δ)
-    final_e[energy_idx].append(energies.value[strobe_indices[idx]])
-    final_e_error[energy_idx].append(energies.σ[strobe_indices[idx]])
 
     # fs.plot_energy_overview(model, ensemble_args=ensemble_arg)
     # fig, ax = plt.subplots()
@@ -69,3 +87,44 @@ for (model, data), Δ, strobe_t, strobe_indices in zip(
 #     ax.errorbar(Δs, energy, σ, label=rf"$\alpha(0)$ = {δs[i]}")
 # ax.legend()
 # fs.export_fig("delta_dependence")
+
+normed_power = power.copy()
+for i in range(dim_models):
+    normed_power[:, i] /= power[:, i].max()
+
+for i in range(dim_models):
+    plt.plot(Δs, power[:, i])
+
+f, (a, a2) = plt.subplots(ncols=2, sharey=True)
+im = a.imshow(
+    power,
+    interpolation="gaussian",
+    interpolation_stage="data",
+    aspect="auto",
+    origin="lower",
+    cmap="plasma",
+    extent=[δs.min(), δs.max(), min(Δs), max(Δs)],
+)
+f.colorbar(im, ax=a)
+
+im2 = a2.imshow(
+    normed_power,
+    interpolation="gaussian",
+    interpolation_stage="data",
+    aspect="auto",
+    origin="lower",
+    cmap="plasma",
+    extent=[δs.min(), δs.max(), min(Δs), max(Δs)],
+)
+f.colorbar(im2, ax=a2)
+a.set_ylabel(r"$\Delta$")
+for ax in (a, a2):
+    ax.set_xlabel(r"$\alpha_\beta(0)$")
+
+a.set_title(r"$\bar{P}$")
+a2.set_title(r"$\bar{P}/\bar{P}_{\mathrm{\max},\alpha}$")
+fs.export_fig("power_heatmap", f, tikz=False)
+
+for model, data in aux.model_data_iterator(Δ_models):
+    inter = model.interaction_energy(data)
+    print(abs(inter.value).max())