some updates on 8

python/energy_flow_proper/08_dynamic_one_bath/cutoff_dep.py

@@ -0,0 +1,247 @@
+import energy_shovel as es
+import numpy as np
+import qutip as qt
+import matplotlib.pyplot as plt
+import utilities as ut
+import figsaver as fs
+import hiro_models.model_auxiliary as aux
+import plot_utils as pu
+
+import ray
+ray.shutdown()
+ray.init(address="auto")
+
+from hops.util.logging_setup import logging_setup
+import logging
+logging_setup(logging.INFO, show_stocproc=False)
+
+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])
+
+for model in ω_models:
+      plt.plot(model.t, model.bcf_scale * model.full_thermal_bcf(model.t).real, label=model.ω_c)
+plt.legend()
+
+for model in ω_models:
+      vals = abs(model.full_thermal_bcf(model.t))
+      plt.plot(model.t, vals / vals.max(), 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.legend()
+
+bath_scales = [1/(model.bcf_coefficients()[1][0].real.min()) for model in ω_models]
+therm_scales = [Δ/(model.bcf_scale) for model in ω_models]
+plt.plot(ωs, bath_scales)
+#plt.plot(ωs, therm_scales)
+
+fig, ax = plt.subplots()
+for model, data in aux.model_data_iterator(ω_models):
+    _, _, bar = pu.plot_with_σ(
+        model.t,
+        model.interaction_energy(data).for_bath(0),
+        ax=ax,
+        label=fr"$\omega_c={model.ω_c:.2f}$",
+    )
+
+fig, ax = plt.subplots()
+import scipy
+for model, data in aux.model_data_iterator(ω_models):
+    _, _, bar = pu.plot_with_σ(
+        model.t,
+        model.total_energy_from_power(data) * (1/es.ergo(model.T)),
+        ax=ax,
+        label=fr"$\omega_c={model.ω_c:.2f}$",
+        strobe_frequency=Δ,
+        markersize=3,
+    )
+
+    pu.plot_with_σ(
+        model.t,
+        (model.total_energy_from_power(data) * (1/es.ergo(model.T))),
+        ax=ax,
+        # label=fr"$ω_c={model.ω_c:.2f}$",
+        # strobe_frequency=Δ,
+        linewidth=0.5,
+        alpha=0.2,
+        color=bar.lines[0].get_color(),
+    )
+
+    τ_off = 1 / (model.bcf_coefficients()[1][0].real.min())
+    a0 = abs(model.bcf(0))
+    τ_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,
+    )
+
+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.legend()
+#fs.export_fig("omegas_total")
+
+fig, ax = plt.subplots()
+for model, data in aux.model_data_iterator(ω_models):
+    fs.plot_with_σ(
+        model.t,
+        model.system_energy(data),
+        ax=ax,
+        label=fr"$ω_c={model.ω_c:.2f}$",
+        strobe_frequency=Δ,
+    )
+ax.legend()
+
+fig, ax = plt.subplots()
+for model,data in aux.model_data_iterator(ω_models):
+      fs.plot_with_σ(model.t, model.system_energy(data) + model.bath_energy(data), bath=0, ax=ax, label=fr"$ω_c={model.ω_c}$", strobe_frequency=Δ)
+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):
+    fs.plot_with_σ(
+        model.t,
+        (model.interaction_energy(data) - int_0) * (1 / abs(int_0).max.value),
+        ax=ax,
+        bath=0,
+        label=fr"$ω_c={model.ω_c:.2f}$",
+        linewidth=0.5,
+    )
+#ax.legend()
+
+fig, ax = plt.subplots()
+fig.set_size_inches(fs.get_figsize("poster", 0.49))
+with aux.get_data(ω_models[0]) as data:
+    first_flow = ω_models[0].bath_energy_flow(data)
+
+for model, data in aux.model_data_iterator(ω_models[1:]):
+    fs.plot_with_σ(
+        model.t,
+        (model.bath_energy_flow(data) - first_flow) * (1 / abs(first_flow.value.mean())),
+        ax=ax,
+        bath=0,
+        label=fr"$ω_c={model.ω_c:.2f}$",
+        linewidth=1,
+    )
+
+ax.set_xlabel(r"$\tau$")
+ax.set_ylabel(r"$(J-J_\mathrm{ref})/\overline{J_\mathrm{ref}}$")
+ax.plot(
+    model.t,
+    model.L.operator_norm(model.t) / 10,
+    linewidth=1,
+    label=r"$||H_I||$",
+    linestyle="--",
+    color="lightblue",
+)
+ax.legend()
+# ax.set_yscale("symlog", linthresh=0.01)
+
+fig, ax = plt.subplots()
+fig.set_size_inches(fs.get_figsize("poster", 0.49))
+max_flows = []
+max_errs = []
+min_flows = []
+min_errs = []
+mean_flows = []
+mean_errs = []
+max_diff = []
+max_diff_err = []
+fig.set_size_inches(fs.get_figsize("poster", .5))
+with aux.get_data(ω_models[0]) as data:
+    model = ω_models[0]
+
+    first_mean = abs(model.bath_energy_flow(data).for_bath(0).mean.value)
+    first_min = abs(model.bath_energy_flow(data).for_bath(0).min.value)
+    first_max = abs(model.bath_energy_flow(data).for_bath(0).max.value)
+    int_0 = abs(model.interaction_energy(data)).for_bath(0).mean.value
+
+
+for model, data in aux.model_data_iterator(ω_models):
+    minimum = model.bath_energy_flow(data).for_bath(0).min * (1 / first_min)
+
+    min_flows.append(minimum.value)
+    min_errs.append(minimum.σ)
+
+    maximum = model.bath_energy_flow(data).for_bath(0).max * (1 / first_max)
+    max_flows.append(maximum.value)
+    max_errs.append(maximum.σ)
+
+    mean = model.bath_energy_flow(data).mean * (1 / first_mean)
+    mean_flows.append(mean.value)
+    mean_errs.append(mean.σ)
+
+    diff = (abs(model.interaction_energy(data).for_bath(0)) * (1 / int_0)).mean
+    max_diff.append(diff.value)
+    max_diff_err.append(diff.σ)
+
+ax.errorbar(ωs, min_flows, min_errs, label="Minimum Flow")
+ax.errorbar(ωs, mean_flows, mean_errs, label="Mean Flow")
+ax.errorbar(ωs, max_flows, max_errs, label="Maximum Flow")
+ax.errorbar(ωs, max_diff, diff.σ, label="Mean Absolute Interaction")
+ax.set_xlabel(r"$\omega_c$")
+ax.legend()
+fs.export_fig("flow_interaction_overview")
+
+for model,data in aux.model_data_iterator(ω_models):
+      #fs.plot_with_σ(model.t, model.bath_energy(data), ax=ax, bath=0, label=fr"$ω_c={model.ω_c}$")
+      fig, ax = fs.plot_energy_overview(model, strobe_frequency=Δ)
+      #fs.plot_with_σ(model.t, model.total_energy(data), ax=ax, bath=0, label=fr"$ω_c={model.ω_c}$")
+      ax.plot(model.t, relative_entropy(data, strobe_indices[-1])[0])
+      fs.plot_with_σ(model.t, EnsembleValue(entropy(data)), ax=ax)
+ax.legend()
+
+model = ω_models[0]
+fig, ax = plt.subplots()
+fig.set_size_inches(fs.get_figsize(448.13094, .5))
+ax.set_xlabel(r"$\tau$")
+ax.set_ylabel("Energy")
+with aux.get_data(model) as data:
+    ax.plot(model.t, model.L.operator_norm(model.t) / 4, label=r"Coupling", linestyle="dotted")
+    #fs.plot_with_σ(model.t, model.bath_energy_flow(data).for_bath(0), label=r"Flow", ax=ax, alpha=.5)
+    fs.plot_with_σ(model.t, model.interaction_energy(data).for_bath(0), label=r"Interaction", ax=ax)
+    fs.plot_with_σ(model.t, model.total_energy(data), label=r"Total", ax=ax, strobe_frequency=Δ, marker="D", markersize=2)
+
+ax.legend()
+ax.set_xlim((0,10))
+fs.export_fig("energy_shovel_preview")
+
+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.legend()
+
+fig, ax = plt.subplots()
+for model,data in aux.model_data_iterator(ω_models):
+          ent = EnsembleValue(relative_entropy(data, strobe_indices[-1]))
+          ent = ent * (1/ent.value.max())
+          fs.plot_with_σ(model.t, ent, ax=ax, label=model.ω_c)
+ax.legend()
+
+final_states = []
+for model, data in aux.model_data_iterator(ω_models):
+    final_states.append(data.rho_t_accum.mean[strobe_indices[-1]])
+
+rel_entropies = [
+    relative_entropy_single(final_states[i], final_states[0])
+    for i in range(len(final_states))
+]
+plt.plot(rel_entropies)
+
+final_diag = []
+final_offdiag = []
+for model, data in aux.model_data_iterator(ω_models):
+    final_diag.append(abs(data.rho_t_accum.mean[strobe_indices[-1]][0,0] - 1/(1+np.exp(1/model.T))))
+    final_offdiag.append(abs(data.rho_t_accum.mean[strobe_indices[-1]][0,1]))
+
+
+plt.plot(ωs, final_offdiag, label=r"$|ρ_{01}|$")
+plt.plot(ωs, final_diag, label=r"$|ρ_{00}-1/(e^{βω}+1)|$")
+plt.legend()

python/energy_flow_proper/08_dynamic_one_bath/fric_v_nofric.py

@@ -42,7 +42,6 @@ with_friction_y, strobe_t, strobe_indices = es.energy_shovel(
 with_friction_y.description = r"With Friction $\sigma_y$"
 
 models = [with_friction, with_friction_y, without_friction]
-# print(es.models_table(models))
 
 aux.integrate_multi(models, 2000)
 

python/energy_flow_proper/08_dynamic_one_bath/mod_freq.py

@@ -0,0 +1,71 @@
+import energy_shovel as es
+import numpy as np
+import qutip as qt
+import matplotlib.pyplot as plt
+import utilities as ut
+import figsaver as fs
+import hiro_models.model_auxiliary as aux
+import plot_utils as pu
+
+import ray
+ray.shutdown()
+ray.init(address="auto")
+
+from hops.util.logging_setup import logging_setup
+import logging
+logging_setup(logging.INFO, show_stocproc=False)
+
+from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix
+
+
+Δs = np.linspace(1,10,20)#np.sort(np.concatenate((np.linspace(1, 5, 20), np.linspace(5, 7, int(20/5 * 2)))))
+Δ_models = []
+strobe_ts = []
+strobe_indices_s = []
+δs = [2, 1,.5, .1]
+for Δ in Δs:
+    for δ in δs:
+        proto, strobe_t, strobe_indices = es.energy_shovel(Δ, periods=20, k_max=5, modulate_system=False)
+        proto.δ = δ
+        strobe_ts.append(strobe_t)
+        strobe_indices_s.append(strobe_indices)
+        Δ_models.append(proto)
+
+aux.integrate_multi(Δ_models, 10_000)
+
+final_e = [[], []]
+final_e_error = [[], []]
+ensemble_arg = dict(overwrite_cache=False)
+fig, ax = plt.subplots()
+for (model, data), Δ, strobe_t, strobe_indices in zip(
+    aux.model_data_iterator(Δ_models),
+    np.array([[Δ]*len(δs) for Δ in Δs]).flatten(),
+    strobe_ts,
+    strobe_indices_s,
+):
+    energies = model.total_energy_from_power(data, **ensemble_arg)
+    idx = energies.value[strobe_indices].argmin()
+    energies = energies * (1 / strobe_t[idx])
+
+    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()
+    # fs.plot_with_σ(model.t, model.energy_change_from_interaction_power(data, **ensemble_arg).for_bath(0), ax=ax)
+    # fs.plot_with_σ(model.t, model.energy_change_from_system_power(data, **ensemble_arg), ax=ax)
+    # fs.plot_with_σ(model.t, model.total_energy_from_power(data, **ensemble_arg), ax=ax)
+    # fs.plot_with_σ(model.t,  model.interaction_power(data, **ensemble_arg).for_bath(0), ax=ax)
+    # fs.plot_with_σ(model.t, model.total_energy(data, **ensemble_arg), ax=ax)
+    # ax.plot(model.t, model.L.derivative()(model.t)[:,0,1])
+    # ax.plot(model.t, model.L(model.t)[:,0,1])
+    # print(strobe_t[1])
+
+plt.ylabel(r"$\frac{\Delta E}{T}$")
+ax.set_xlabel(r"$\Delta$")
+ax.set_ylabel(r"$P_\mathrm{max}$")
+for i, energy, σ in zip(itertools.count(0), final_e, final_e_error):
+    ax.errorbar(Δs, energy, σ, label=rf"$\alpha(0)$ = {δs[i]}")
+ax.legend()
+fs.export_fig("delta_dependence")

python/energy_flow_proper/08_dynamic_one_bath/mount_taurus.sh

@@ -0,0 +1,6 @@
+#!/usr/bin/env bash
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_08/project/python/energy_flow_proper/08_dynamic_one_bath/.data/ .data_taurus
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_08/project/python/energy_flow_proper/08_dynamic_one_bath/results results_taurus
+
+
+sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus  s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_08/project/python/energy_flow_proper/08_dynamic_one_bath/ taurus

python/energy_flow_proper/08_dynamic_one_bath/night.py

@@ -1,145 +1,55 @@
-# [[file:shovel_experiments.org::*Dependence on Cutoff][Dependence on Cutoff:1]]
-ωs = np.unique(np.sort(np.concatenate([[.5], np.linspace(1, 1.8, 6)])))
+import scipy
+import pathlib
+
+ωs = np.linspace(.5, 2.5, 6)
 ω_models = []
+Δ = 5
 
 for ω_c in ωs:
-     proto, strobe_t, strobe_indices = energy_shovel(Δ, periods=30, modulate_system=False)
-     proto.δ = 2 / ω_c
-     proto.ω_c = ω_c
-     proto.ω_s = 1 + Δ - ω_c
-     proto.k_max = 7
-     proto.therm_method = "fft"
+    proto, strobe_t, strobe_indices = es.energy_shovel(
+        Δ,
+        periods=3,
+        modulate_system=False,
+        δ=0.1,
+        ω_c = ω_c,
+    )
 
-     ω_models.append(proto)
+    proto.k_max = 5
+    proto.therm_method = "fft"
 
-fs.tex_value(
-    ω_models[1].bcf_scale * ω_models[1].bcf(0).real,
-    prec=1,
-    prefix="α(0)=",
-    save="omega_alpha",
-), fs.tex_value(
-    Δ,
-    prec=1,
-    prefix="Δ=",
-    save="omega_delta",
-), fs.tex_value(
-    ω_models[1].T,
-    prec=1,
-    prefix="T=",
-    save="omega_delta",
-)
-# Dependence on Cutoff:1 ends here
+    goal = 0.4
 
-# [[file:shovel_experiments.org::*Dependence on Cutoff][Dependence on Cutoff:2]]
-aux.integrate_multi(ω_models, 10000)
-# Dependence on Cutoff:2 ends here
 
-# [[file:shovel_experiments.org::*Dependence on Cutoff][Dependence on Cutoff:3]]
+    def cost(δ):
+        proto.δ = δ
+        aux.integrate(proto, 200)
+        with aux.get_data(proto) as data:
+            inter = proto.interaction_energy(data)
+
+        return (abs(inter.value).max() - goal) ** 2
+    cache_path = pathlib.Path(f"./.cache/{proto.hexhash}_delta_{goal}.npy")
+    if cache_path.exists():
+        with cache_path.open("rb") as cache:
+            δ = np.load(cache, allow_pickle=True)
+    else:
+        δ = scipy.optimize.minimize_scalar(cost,  bounds=(.1, 3), method="bounded", options=dict(xatol=1e-2, disp=3)).x
+        cache_path.parent.mkdir(parents=True, exist_ok=True)
+        with cache_path.open("wb") as cache:
+            np.save(cache, δ)
+
+    final, strobe_t, strobe_indices = es.energy_shovel(
+        Δ,
+        periods=40,
+        modulate_system=False,
+        δ=δ,
+        ω_c=ω_c
+    )
+    final.k_max = 5
+    final.therm_method = "fft"
+
+    ω_models.append(final.copy())
+
+aux.integrate_multi(ω_models, 10_000)
+
 from pprint import pprint
 print_diff_tbl(ω_models)
-# Dependence on Cutoff:3 ends here
-
-# [[file:shovel_experiments.org::*Energy stuff][Energy stuff:1]]
-fig, ax = plt.subplots()
-fig.set_size_inches(fs.get_figsize("poster", 0.49))
-for model, data in aux.model_data_iterator(ω_models):
-    _, _, bar = fs.plot_with_σ(
-        model.t,
-        model.total_energy_from_power(data) * (1/ergo(model.T)),
-        ax=ax,
-        label=fr"$ω_c={model.ω_c:.2f}$",
-        strobe_frequency=Δ,
-        markersize=3,
-    )
-
-    fs.plot_with_σ(
-        model.t[::5],
-        (model.total_energy_from_power(data) * (1/ergo(model.T))).slice(None, None, 5),
-        ax=ax,
-        # label=fr"$ω_c={model.ω_c:.2f}$",
-        # strobe_frequency=Δ,
-        linewidth=0.5,
-        alpha=0.2,
-        color=bar.lines[0].get_color(),
-    )
-
-    τ_off = 1 / (model.bcf_coefficients()[1][0].real.min())
-    a0 = abs(model.bcf(0))
-    τ_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,
-    )
-
-ax.set_xlabel(r"$\tau$")
-ax.set_ylabel(r"$(\langle H\rangle_\tau -\langle H\rangle_0)/\Delta E_\mathrm{max}$")
-# ax.axhline(, color="grey", linewidth=1, linestyle="dashed")
-# ax.plot(model.t, model.L.operator_norm(model.t), linewidth=1)
-ax.legend()
-fs.export_fig("omegas_total")
-# Energy stuff:1 ends here
-
-# [[file:shovel_experiments.org::*Modulation Frequency Dependence][Modulation Frequency Dependence:1]]
-from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix
-
-
-Δs = np.linspace(1,10,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 = [1,.5]
-for Δ in Δs:
-    for δ in δs:
-        proto, strobe_t, strobe_indices = energy_shovel(Δ, periods=10, k_max=5, modulate_system=False)
-        proto.δ = δ
-        strobe_ts.append(strobe_t)
-        strobe_indices_s.append(strobe_indices)
-        Δ_models.append(proto)
-# Modulation Frequency Dependence:1 ends here
-
-# [[file:shovel_experiments.org::*Modulation Frequency Dependence][Modulation Frequency Dependence:2]]
-aux.integrate_multi(Δ_models, 1000)
-# Modulation Frequency Dependence:2 ends here
-
-# [[file:shovel_experiments.org::*Modulation Frequency Dependence][Modulation Frequency Dependence:3]]
-final_e = [[], []]
-final_e_error = [[], []]
-ensemble_arg = dict(overwrite_cache=False)
-fig, ax = plt.subplots()
-fig.set_size_inches(fs.get_figsize("poster", 0.49))
-for (model, data), Δ, strobe_t, strobe_indices in zip(
-    aux.model_data_iterator(Δ_models),
-    np.array([[Δ]*len(δs) for Δ in Δs]).flatten(),
-    strobe_ts,
-    strobe_indices_s,
-):
-    energies = model.total_energy_from_power(data, **ensemble_arg)
-    idx = energies.value[strobe_indices].argmin()
-    energies = energies * (1 / strobe_t[idx])
-
-    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()
-    # fs.plot_with_σ(model.t, model.energy_change_from_interaction_power(data, **ensemble_arg).for_bath(0), ax=ax)
-    # fs.plot_with_σ(model.t, model.energy_change_from_system_power(data, **ensemble_arg), ax=ax)
-    # fs.plot_with_σ(model.t, model.total_energy_from_power(data, **ensemble_arg), ax=ax)
-    # fs.plot_with_σ(model.t,  model.interaction_power(data, **ensemble_arg).for_bath(0), ax=ax)
-    # fs.plot_with_σ(model.t, model.total_energy(data, **ensemble_arg), ax=ax)
-    # ax.plot(model.t, model.L.derivative()(model.t)[:,0,1])
-    # ax.plot(model.t, model.L(model.t)[:,0,1])
-    # print(strobe_t[1])
-
-plt.ylabel(r"$\frac{\Delta E}{T}$")
-ax.set_xlabel(r"$\Delta$")
-ax.set_ylabel(r"$P_\mathrm{max}$")
-for i, energy, σ in zip(itertools.count(0), final_e, final_e_error):
-    ax.errorbar(Δs, energy, σ, label=rf"$\alpha(0)$ = {δs[i]}")
-ax.legend()
-fs.export_fig("delta_dependence")
-# Modulation Frequency Dependence:3 ends here

python/energy_flow_proper/08_dynamic_one_bath/sys_v_no_sys.py

@@ -38,7 +38,12 @@ without_system.description = "No System Modulation"
 
 models = [without_system, with_system]
 
-aux.integrate_multi(models, 10_000)
+aux.integrate_multi(models, 10_00)
 
 f, _ = es.energy_friction_plot(models, Δ, strobe_indices)
 fs.export_fig("system_vs_no_system", f)
+
+print(es.models_table(models))
+
+f, ax, _ = pu.plot_ρ(models[0], 0, 1)
+pu.plot_ρ(models[0], 0, 1, strobe_frequency=Δ, ax=ax, markersize=5)