diff --git a/python/energy_flow_proper/08_dynamic_one_bath/coupling_modulation.org b/python/energy_flow_proper/08_dynamic_one_bath/coupling_modulation.org index 88c1d70..447e25d 100644 --- a/python/energy_flow_proper/08_dynamic_one_bath/coupling_modulation.org +++ b/python/energy_flow_proper/08_dynamic_one_bath/coupling_modulation.org @@ -13,8 +13,6 @@ Here we try the same thing as in [[file:coupling_and_system_cycle.org][here]] bu import scipy import utilities as ut import matplotlib - matplotlib.rcParams.update(fs.MPL_RC_POSTER) - import itertools #+end_src @@ -170,8 +168,6 @@ Init ray and silence stocproc. strobe_t, strobe_indices = ut.strobe_times(model.t, Δ) return model, strobe_t, strobe_indices - - #+end_src #+RESULTS: @@ -2138,7 +2134,7 @@ effect persists. Ok, it was the exponents in the bcf expansion :P. [[file:./.ob-jupyter/c9f9872a8dfae5d029676840920eeaaff011fc40.svg]] :END: -* TODO Model: Testing the Inital Slip +* Testing the Inital Slip #+begin_src jupyter-python from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix, Piecewise diff --git a/python/energy_flow_proper/08_dynamic_one_bath/energy_shovel.py b/python/energy_flow_proper/08_dynamic_one_bath/energy_shovel.py new file mode 100644 index 0000000..622c015 --- /dev/null +++ b/python/energy_flow_proper/08_dynamic_one_bath/energy_shovel.py @@ -0,0 +1,66 @@ +import figsaver as fs +from hiro_models.one_qubit_model import QubitModel, StocProcTolerances +import hiro_models.model_auxiliary as aux +from hops.util.utilities import relative_entropy, relative_entropy_single, entropy, trace_distance +from hopsflow.util import EnsembleValue +import numpy as np +import qutip as qt +import scipy +import utilities as ut +import matplotlib +import itertools + +from hops.util.dynamic_matrix import ( + SmoothStep, + Periodic, + Harmonic, + ConstantMatrix, + ScaleTime, +) + +def energy_shovel( + Δ=2, + γ=1 / 10, + periods=5, + L_op=(1 / 2 * qt.sigmax()).full(), + H_op=1 / 2 * (qt.sigmaz() + 1).full(), + H_time_op=qt.sigmay().full() / 2, + dt=.01, + detune=0, + ω_c=1, + modulate_system=True, + k_max=7, +): + t_max = 2*np.pi/Δ * periods + L = ConstantMatrix(L_op) - Harmonic( + L_op, Δ, np.pi / 2 + ) + H = ConstantMatrix(H_op) + if modulate_system: + H = H + γ * Δ * Harmonic(H_time_op, Δ) + + model = QubitModel( + δ=1, + ω_c=ω_c, + ω_s=1 + Δ - ω_c + detune, + t=ut.linspace_with_strobe(0, t_max, int(t_max/dt + 1), Δ), + ψ_0=(0.5 * 0 * qt.basis([2], [0]) + qt.basis([2], [1])), + description=f"Testing the time dependent coupling with smooth step.", + k_max=k_max, + bcf_terms=7, + truncation_scheme="simplex", + driving_process_tolerance=StocProcTolerances(1e-5, 1e-5), + thermal_process_tolerance=StocProcTolerances(1e-5, 1e-5), + T=5, + L=L, + H=H, + therm_method="fft", + ) + + strobe_t, strobe_indices = ut.strobe_times(model.t, Δ) + + return model, strobe_t, strobe_indices + +def ergo(T, ω=1): + """The ergotropy of a qubit coupled to a bath.""" + return (T * np.log(1 + np.exp(-ω / T))) diff --git a/python/energy_flow_proper/08_dynamic_one_bath/flake.lock b/python/energy_flow_proper/08_dynamic_one_bath/flake.lock index b0e3ac6..5b1c002 100644 --- a/python/energy_flow_proper/08_dynamic_one_bath/flake.lock +++ b/python/energy_flow_proper/08_dynamic_one_bath/flake.lock @@ -2,11 +2,11 @@ "nodes": { "flake-utils": { "locked": { - "lastModified": 1656065134, - "narHash": "sha256-oc6E6ByIw3oJaIyc67maaFcnjYOz1mMcOtHxbEf9NwQ=", + "lastModified": 1659877975, + "narHash": "sha256-zllb8aq3YO3h8B/U0/J1WBgAL8EX5yWf5pMj3G0NAmc=", "owner": "numtide", "repo": "flake-utils", - "rev": "bee6a7250dd1b01844a2de7e02e4df7d8a0a206c", + "rev": "c0e246b9b83f637f4681389ecabcb2681b4f3af0", "type": "github" }, "original": { @@ -17,11 +17,11 @@ }, "flake-utils_2": { "locked": { - "lastModified": 1656065134, - "narHash": "sha256-oc6E6ByIw3oJaIyc67maaFcnjYOz1mMcOtHxbEf9NwQ=", + "lastModified": 1659877975, + "narHash": "sha256-zllb8aq3YO3h8B/U0/J1WBgAL8EX5yWf5pMj3G0NAmc=", "owner": "numtide", "repo": "flake-utils", - "rev": "bee6a7250dd1b01844a2de7e02e4df7d8a0a206c", + "rev": "c0e246b9b83f637f4681389ecabcb2681b4f3af0", "type": "github" }, "original": { @@ -32,11 +32,11 @@ }, "nixpkgs": { "locked": { - "lastModified": 1656753965, - "narHash": "sha256-BCrB3l0qpJokOnIVc3g2lHiGhnjUi0MoXiw6t1o8H1E=", + "lastModified": 1660485612, + "narHash": "sha256-sSLW1KaB1adKTJn9+Ja3h3AaS7QCZyhUKiSUStcLg80=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "0ea7a8f1b939d74e5df8af9a8f7342097cdf69eb", + "rev": "6512b21eabb4d52e87ea2edcf31a288e67b2e4f8", "type": "github" }, "original": { @@ -47,11 +47,11 @@ }, "nixpkgs_2": { "locked": { - "lastModified": 1656239181, - "narHash": "sha256-wW1xRFBn376yGloXZ4QzBE4hjipMawpV18Lshd9QSPw=", + "lastModified": 1660396586, + "narHash": "sha256-ePuWn7z/J5p2lO7YokOG1o01M0pDDVL3VrStaPpS5Ig=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "f2537a505d45c31fe5d9c27ea9829b6f4c4e6ac5", + "rev": "e105167e98817ba9fe079c6c3c544c6ef188e276", "type": "github" }, "original": { @@ -62,11 +62,11 @@ }, "nixpkgs_3": { "locked": { - "lastModified": 1656492412, - "narHash": "sha256-CdTRqPFyE/t/nRfnfN3fxm+OOKjeJdaDbXm7E7k+Fxc=", + "lastModified": 1660546381, + "narHash": "sha256-rEzCjeWVGhK5AyHxm1zet0lF6+AVSW3JuU5LAU2SMYU=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "72f17cca89ec773e5845f5fa1bffcf0b7a4e2094", + "rev": "eb642f80f9aecc19312909e08601a3c2020b5ce2", "type": "github" }, "original": { @@ -81,11 +81,11 @@ "nixpkgs": "nixpkgs_3" }, "locked": { - "lastModified": 1656405165, - "narHash": "sha256-p2c9QAnKsiUu0a3iHAkpHkO9jAkMWCbG1HCU7/TrYc0=", + "lastModified": 1660542485, + "narHash": "sha256-XfklMwJMLB7bLI5ZnQTrNaK7KyBnElLGoWOL++XO3zk=", "owner": "nix-community", "repo": "poetry2nix", - "rev": "abc47c71a4920e654e7b2e4261e3e6399bbe2be6", + "rev": "8b6239cf2ded121f8f3d570d2563d69f05d4208f", "type": "github" }, "original": { @@ -107,11 +107,11 @@ "poetry2nix": "poetry2nix" }, "locked": { - "lastModified": 1656499586, - "narHash": "sha256-tx99W8VI8xxsnPRK0zZs5V2xFo5/hfhxhrTxCQi/mQ8=", + "lastModified": 1660547406, + "narHash": "sha256-ZRLGxqHmX5xtCBPCt1pb2E3Z70IEmMT3TA+JruosF4Q=", "owner": "vale981", "repo": "hiro-flake-utils", - "rev": "3c96df6681a2475d131a5f30e7cb02c223b7c567", + "rev": "dd1e89c5fabbca9a3dbec91cab3e3e3867da10fc", "type": "github" }, "original": { diff --git a/python/energy_flow_proper/08_dynamic_one_bath/shovel_experiments.org b/python/energy_flow_proper/08_dynamic_one_bath/shovel_experiments.org new file mode 100644 index 0000000..f39cac2 --- /dev/null +++ b/python/energy_flow_proper/08_dynamic_one_bath/shovel_experiments.org @@ -0,0 +1,1344 @@ +#+PROPERTY: header-args :session 08_shovel :kernel python :pandoc no :async yes :tangle no + +#+begin_src jupyter-python :results none + 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 +#+end_src + +Init ray and silence stocproc. +#+begin_src jupyter-python :tangle night.py + import ray + ray.shutdown() + ray.init(address="auto") +#+end_src + + +#+begin_src jupyter-python :results none :tangle night.py + from hops.util.logging_setup import logging_setup + import logging + logging_setup(logging.INFO, show_stocproc=False) +#+end_src + +* Shovel +** Fricionless vs Friction +- we modulate only the system, once with σx once with σz +#+begin_src jupyter-python :results none + proto, strobe_t, strobe_indices = energy_shovel(Δ, periods=30, modulate_system=False) + + for ω_c in ωs: + + proto.δ = 2 / ω_c + proto.ω_c = ω_c + proto.ω_s = 1 + Δ - ω_c + proto.k_max = 7 + proto.therm_method = "fft" + + ω_models.append(proto) + +#+end_src +** No Bath mod vs Bath mod +** No System mod vs System Mod + +* Only Coupling Modulation +** Dependence on Cutoff +- found that with low temperature only a small cutoff produces a net + energy gain + +#+begin_src jupyter-python :tangle night.py + ωs = np.unique(np.sort(np.concatenate([[.5], np.linspace(1, 1.8, 6)]))) + ω_models = [] + + 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" + + ω_models.append(proto) + + 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", + ) +#+end_src + +#+RESULTS: +| \(α(0)=1.0\) | \(Δ=5.0\) | \(T=5.0\) | + +#+begin_src jupyter-python :tangle night.py + aux.integrate_multi(ω_models, 10000) +#+end_src + +#+RESULTS: +#+begin_example + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. + 0it [00:00, ?it/s] + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. +0it [00:00, ?it/s] + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. +0it [00:00, ?it/s] + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. +0it [00:00, ?it/s] + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. +0it [00:00, ?it/s] + [INFO hops.core.integration 746376] Choosing the nonlinear integrator. + [INFO hops.core.integration 746376] Using 21 integrators. + [INFO hops.core.integration 746376] Some 5292 trajectories have to be integrated. + [INFO hops.core.integration 746376] Using 3432 hierarchy states. + 0% 0/5292 [00:00 | +[[file:./.ob-jupyter/649456f5f15190bf9b35f20208fe8b1f2cb9df73.svg]] +:END: + +#+begin_src jupyter-python + for model in ω_models: + plt.plot(model.t, model.bcf_scale * model.bcf(model.t).real, label=model.ω_c) + plt.legend() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/edb55c566d04f86a831b264320b519f71f35b6ae.svg]] +:END: + +#+begin_src jupyter-python + for model in ω_models: + plt.plot(model.t, model.bcf_scale * (model.bcf(model.t) + 2* model.thermal_correlations(model.t)).real, linewidth=1) + plt.xlim(-.1, 6) +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/31ca683a025fb7cf6bed55cd9d8acf77746c9cd6.svg]] + +The time scale of the bath and the modulation is also interesting. + +#+begin_src jupyter-python + 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) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/4dbbb5d116c58927a7b40c32ba1a75696f8327ef.svg]] +:END: + +*** Energy stuff +#+begin_src jupyter-python :tangle night.py + 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") +#+end_src + +#+RESULTS: +:RESULTS: +[[file:./.ob-jupyter/8f1bed759f6ce8711ccd62e466af732f4e417984.svg]] +[[file:./.ob-jupyter/910357c65d8a55952256eec9476985cab4225257.svg]] +:END: + + #+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/6776ea0f3456a859f2ee41eed4245eb4bc985e7e.svg]] +:END: + + + +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/185d7cb94397008b5746bb2b06f25509780904d6.svg]] +:END: + +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/dada909ebf98f84e3ab6ff83510d6e169fe46d1e.svg]] + +- pretty similar +- substantial :P + + +#+begin_src jupyter-python + 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) +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/01514106235220e693d4580a46e48e54efb24075.svg]] +:END: + + +#+begin_src jupyter-python + 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") +#+end_src + +#+RESULTS: +:RESULTS: +[[file:./.ob-jupyter/b92b3c59f650ebb74ff6b23dae048cbabc060b24.svg]] +[[file:./.ob-jupyter/c151c4405d329b0b8ffa802707e2ba269b5b2a8b.svg]] +:END: + + +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +:RESULTS: +: /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. +: warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) +: /home/hiro/src/hops/hops/util/utilities.py:114: RuntimeWarning: invalid value encountered in sqrt +: np.sqrt(np.trace(_square_mat(Δρ @ (logm2(ρ) + 1 - logm2(ρ_ref))).real)) +: /home/hiro/src/hops/hops/util/utilities.py:43: RuntimeWarning: invalid value encountered in sqrt +: np.sqrt(np.trace(_square_mat(Δρ @ (logm2(ρ) + 1))).real) +: +[[file:./.ob-jupyter/84436d134701d5e37b89e467ef94173d0e8be9fa.svg]] +[[file:./.ob-jupyter/c0a9f52f580f8ed08d60f31357e2503145284592.svg]] +[[file:./.ob-jupyter/d24c21c2d3996121f9c075fc4bd540b88cb7abec.svg]] +[[file:./.ob-jupyter/e313ce79770980ea22acc37621e497d3636497f1.svg]] +[[file:./.ob-jupyter/fd8a45f1cb65721b26d2ecf991c22904520f9b53.svg]] +[[file:./.ob-jupyter/554cb98d566f2d180ea3d0f67c11199c892810f4.svg]] +:END: + +#+begin_src jupyter-python + 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") +#+end_src + +#+RESULTS: +:RESULTS: +[[file:./.ob-jupyter/21dd4fb6de7275b03ea6e97caaff4557b94ab7a7.svg]] +[[file:./.ob-jupyter/5b5e7c45fdb0dce47d26fc005c76f10e7f3a80f2.svg]] +:END: + +*** State stuff +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +:RESULTS: +: /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. +: warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) +: /home/hiro/src/hops/hops/util/utilities.py:114: RuntimeWarning: invalid value encountered in sqrt +: np.sqrt(np.trace(_square_mat(Δρ @ (logm2(ρ) + 1 - logm2(ρ_ref))).real)) +: +[[file:./.ob-jupyter/db58e138f71c5d6e18f644e3090e6f3fb6dd5863.svg]] +:END: + +Whole different steady state -> phase transition? + + +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/36bad68d4442eeef3c646d2dae5e6bd4c667cd36.svg]] +:END: + +#+begin_src jupyter-python + 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) +#+end_src + +#+RESULTS: +: d48f5002-7709-471d-af92-f93f1222654a + +#+begin_src jupyter-python + 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() +#+end_src + +#+RESULTS: +: f2137693-bd52-4382-bf39-dbe2238b0924 + +*** TODO Find out what happens when we use the same peak heights + +** Modulation Frequency Dependence +#+begin_src jupyter-python :tangle night.py + 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) +#+end_src + +#+RESULTS: + +#+begin_src jupyter-python :tangle night.py + aux.integrate_multi(Δ_models, 1000) +#+end_src + +#+RESULTS: +#+begin_example + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. + 100% 1000/1000 [04:59<00:00, 3.34it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [04:58<00:00, 3.35it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:52<00:00, 5.80it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:43<00:00, 6.12it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:24<00:00, 6.92it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:13<00:00, 7.47it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:11<00:00, 7.59it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:00<00:00, 8.28it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:59<00:00, 8.35it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:49<00:00, 9.12it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:01<00:00, 8.21it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:46<00:00, 9.38it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:03<00:00, 8.12it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:52<00:00, 8.88it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:07<00:00, 7.84it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:54<00:00, 8.75it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:06<00:00, 7.92it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:54<00:00, 8.72it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [02:01<00:00, 8.24it/s] + [INFO hops.core.integration 27050] Choosing the nonlinear integrator. + [INFO hops.core.integration 27050] Using 21 integrators. + [INFO hops.core.integration 27050] Some 1000 trajectories have to be integrated. + [INFO hops.core.integration 27050] Using 792 hierarchy states. +100% 1000/1000 [01:54<00:00, 8.71it/s] +#+end_example + +#+begin_src jupyter-python :tangle night.py + 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") +#+end_src + +#+RESULTS: +:RESULTS: +: /tmp/ipykernel_27050/88422595.py:14: RuntimeWarning: divide by zero encountered in double_scalars +: energies = energies * (1 / strobe_t[idx]) +: /home/hiro/src/hopsflow/hopsflow/util.py:233: RuntimeWarning: invalid value encountered in multiply +: [(N, val * other, np.abs(σ * other)) for N, val, σ in self._value] +: /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/numpy/core/_methods.py:44: RuntimeWarning: invalid value encountered in reduce +: return umr_minimum(a, axis, None, out, keepdims, initial, where) +: /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/numpy/core/_methods.py:40: RuntimeWarning: invalid value encountered in reduce +: return umr_maximum(a, axis, None, out, keepdims, initial, where) +[[file:./.ob-jupyter/186eb54734f796b02b1ee7f8f37c7392300fb50f.svg]] +[[file:./.ob-jupyter/111ffb226d35d3a9765a8019b3ff0eb4ee7a1311.svg]] +:END: + +Looks very much like a speed limit. Likely due tho lower temperature +suppression. + +** Modulation Tuning +#+begin_src jupyter-python + from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix + + + tunings = np.sort(np.array([0, *np.linspace(-.5, .5, 8), *np.linspace(-.1, .1, 4)])) + tune_models = [] + Δ = 2 + for tune in tunings: + proto, strobe_t, strobe_indices = energy_shovel(Δ, periods=2, modulate_system=False, detune=tune) + + tune_models.append(proto) +#+end_src + +#+RESULTS: + +#+begin_src jupyter-python +aux.integrate_multi(tune_models, 10000) +#+end_src + +#+RESULTS: +:RESULTS: +#+begin_example + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 100% 5000/5000 [03:45<00:00, 22.22it/s] + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 100% 5000/5000 [04:25<00:00, 18.81it/s] + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 100% 5000/5000 [03:46<00:00, 22.05it/s] + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 100% 5000/5000 [04:03<00:00, 20.56it/s] + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 100% 5000/5000 [04:02<00:00, 20.59it/s] + [INFO hops.core.integration 370330] Choosing the nonlinear integrator. + [INFO hops.core.integration 370330] Using 22 integrators. + [INFO hops.core.integration 370330] Some 5000 trajectories have to be integrated. + [INFO hops.core.integration 370330] Using 792 hierarchy states. + 73% 3631/5000 [02:35<00:58, 23.39it/s][INFO hops.core.signal_delay 370330] caught sig 'SIGINT' + [INFO hops.core.signal_delay 370330] caught 1 signal(s) + [INFO hops.core.signal_delay 370330] emit signal 'SIGINT' + [INFO hops.core.signal_delay 370330] caught sig 'SIGINT' + 73% 3637/5000 [02:35<00:58, 23.37it/s] + 2022-06-08 18:39:16,138 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,144 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,147 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,151 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,153 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,158 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,163 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,165 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,174 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,180 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,183 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,186 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,189 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,193 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,197 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,199 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,204 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,213 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,215 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,217 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,223 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + 2022-06-08 18:39:16,228 ERROR worker.py:92 -- Unhandled error (suppress with 'RAY_IGNORE_UNHANDLED_ERRORS=1'): The worker died unexpectedly while executing this task. Check python-core-worker-*.log files for more information. + [INFO hops.core.signal_delay 370330] caught 1 signal(s) + [INFO hops.core.signal_delay 370330] emit signal 'SIGINT' +#+end_example +# [goto error] +#+begin_example + --------------------------------------------------------------------------- + KeyboardInterrupt Traceback (most recent call last) + Input In [303], in () + ----> 1 aux.integrate_multi(tune_models, 10000) + + File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:82, in integrate_multi(models, *args, **kwargs) +  75 """Integrate the hops equations for the ``models``. +  76 Like :any:`integrate` just for many models. +  77 +  78 A call to :any:`ray.init` may be required. +  79 """ +  81 for model in models: + ---> 82 integrate(model, *args, **kwargs) + + File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:106, in integrate(model, n, data_path, clear_pd) +  96 # with model_db(data_path) as db: +  97 # if hash in db and "data" db[hash] +  99 supervisor = HOPSSupervisor( +  100 model.hops_config, +  101 n, +  102 data_path=data_path, +  103 data_name=hash, +  104 ) + --> 106 supervisor.integrate(clear_pd) +  108 with supervisor.get_data(True) as data: +  109 with model_db(data_path) as db: + + File ~/src/hops/hops/core/integration.py:1288, in HOPSSupervisor.integrate(self, clear_pd) +  1285 break +  1287 integration.update() + -> 1288 data.new_samples( +  1289 idx=index, +  1290 incomplete=incomplete, +  1291 psi0=psi0, +  1292 aux_states=aux_states, +  1293 stoc_proc=stoc_proc, +  1294 result_type=self.params.HiP.result_type, +  1295 normed=self._normed_average, +  1296 rng_seed=seed, +  1297 ) + + File ~/src/hops/hops/core/signal_delay.py:87, in sig_delay.__exit__(self, exc_type, exc_val, exc_tb) +  84 if len(self.sigh.sigs_caught) > 0 and self.handler is not None: +  85 self.handler(self.sigh.sigs_caught) + ---> 87 self._restore() + + File ~/src/hops/hops/core/signal_delay.py:68, in sig_delay._restore(self) +  66 for i, s in enumerate(self.sigs): +  67 signal.signal(s, self.old_handlers[i]) + ---> 68 self.sigh.emit() + + File ~/src/hops/hops/core/signal_delay.py:42, in SigHandler.emit(self) +  40 for s in self.sigs_caught: +  41 log.info("emit signal '{}'".format(SIG_MAP[s])) + ---> 42 os.kill(os.getpid(), s) + + KeyboardInterrupt: +#+end_example +:END: + +#+begin_src jupyter-python + final_e = [] + for (model, data), tune in zip(aux.model_data_iterator(tune_models), tunings): + final_e.append((model.system_energy(data) + model.bath_energy(data)).value[0,-1]) + + plt.ylabel(r"$\Delta E$") + plt.xlabel(r"$\delta$") + plt.plot(tunings, final_e, linestyle="none", marker="o") +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/2b1e3a706d93b9972fd199062fa6ba37d0b4e695.svg]] +:END: + +There is an optimum. Slightly away from zero. + +#+begin_src jupyter-python + plt.plot( + tune_models[0].t, + np.abs(tune_models[0].bcf(tune_models[0].t) + + 2 * tune_models[0].thermal_correlations(tune_models[0].t).real), + ) + plt.plot( + tune_models[0].t, + np.abs(tune_models[0].bcf(tune_models[0].t) + + 0 * tune_models[0].thermal_correlations(tune_models[0].t).real), + ) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/51323bae0cfdcfa0608b37c4d9e9d7e067dc17b5.svg]] +:END: + +#+begin_src jupyter-python + ωs = np.linspace(1e-3, 10, 1000) + def nonzero_density(ωs, model): + return model.spectral_density(ωs) * (1 + 1/np.expm1(ωs / model.T)) + for tuning, model in zip(tunings, tune_models): + freqs = np.array([Δ + 1, 2*Δ + 1]) + line = plt.plot(tuning, nonzero_density(freqs[0], model), label=tuning, markersize=2, marker="o") + plt.plot(tuning, nonzero_density(freqs[1], model), color=line[0].get_color(), markersize=4, marker="*") + #plt.plot(ωs, nonzero_density(ωs, model)) + + plt.legend() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/cb820f528d2bd53259cbbc0e4fe1dc2ca37430a9.svg]] +:END: + +- maybe not only first harmonic important ;P + + +*** HOLD Linear With Temperature :HOLD: + +* Misc Projects +** Decay of Auxiliary States +#+begin_src jupyter-python + from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix + + L_op = (1 / 2 * qt.sigmax()).full() + Δ = 2 + γ = 1 / Δ * 1 / 10 + t_max = 60 + L = (ConstantMatrix(L_op) - Harmonic(L_op, Δ, np.pi/2)) @ (ConstantMatrix(np.eye(2)) - SmoothStep(np.eye(2), 5, 6, 2)) + H = 1 / 2 * (ConstantMatrix(qt.sigmaz() + 1)) + γ * Δ * Harmonic(qt.sigmay().full(), Δ) + + model = QubitModel( + δ=1, + ω_c=1, + ω_s=2, + t=ut.linspace_with_strobe(0, t_max, 1000, 2*Δ), + ψ_0=(0.5 * 0 * qt.basis([2], [0]) + qt.basis([2], [1])), + description=f"Testing the time dependent coupling with smooth step.", + k_max=5, + bcf_terms=7, + truncation_scheme="simplex", + driving_process_tolerance=StocProcTolerances(1e-5, 1e-5), + thermal_process_tolerance=StocProcTolerances(1e-4, 1e-4), + T=5, + L=L, + H=H, + therm_method="fft", + ) +#+end_src + +#+RESULTS: + +#+begin_src jupyter-python + aux.integrate(model, 50) +#+end_src + +#+RESULTS: +: [INFO hops.core.integration 234181] Choosing the nonlinear integrator. +: [INFO hops.core.integration 234181] Using 21 integrators. +: [INFO hops.core.integration 234181] Some 40 trajectories have to be integrated. +: [INFO hops.core.integration 234181] Using 792 hierarchy states. +: 100% 40/40 [00:24<00:00, 1.61it/s] + +#+begin_src jupyter-python + mean_norm = np.zeros_like(model.t) + with aux.get_data(model) as data: + for i in range(data.samples): + aux_state = data.aux_states[i, :] + # plt.plot(model.t, np.linalg.norm(aux_state, axis=1)) + + mean_norm += np.linalg.norm(aux_state, axis=1) + + mean_norm /= data.samples + + plt.plot(model.t, mean_norm) + plt.plot(model.t, model.L.operator_norm(model.t)) + #plt.plot(model.t, abs(model.bcf(model.t - 4.9))) + τ = 1/model.bcf_coefficients()[1][0].real.min() + slowest_decay = lambda t, a: a * np.exp(-t / τ) + shift = scipy.optimize.curve_fit(slowest_decay, model.t[-100:], mean_norm[-100:])[0] + plt.plot(model.t, slowest_decay(model.t, shift)) + plt.yscale("log") + plt.ylim((slowest_decay(model.t[-1], shift)/2, 1.1)) +#+end_src + +#+RESULTS: +:RESULTS: +| array | ((0.00106639)) | 1.1 | +[[file:./.ob-jupyter/f2f7e727122b13fed5b494964c654071d04788ac.svg]] +:END: + +They do decay. Slower than the bcf though. +#+begin_src jupyter-python +1/model.bcf_coefficients()[1][0].real.min() +#+end_src + +#+RESULTS: +: 23.371125928961806 + +We find exponential decay with the smallest exponent of the BCF +expansion. Without the expansion in exponentials we would probably see +the real BCF. + +*** TODO Decay weighted by G +** Testing the Inital Slip +#+begin_src jupyter-python + from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix, Piecewise + + L_op = (1 / 2 * qt.sigmax()).full() + + H = 1 / 2 * (ConstantMatrix(qt.sigmaz() + 1)) + + t_max = 2 + models = [ + QubitModel( + δ=0.2, + ω_c=1, + t=np.linspace(0, t_max, 500), + ψ_0=qt.basis([2], [0]), + description=f"Testing the time dependent coupling with smooth step.", + k_max=4, + bcf_terms=7, + truncation_scheme="simplex", + driving_process_tolerance=StocProcTolerances(1e-6, 1e-6), + thermal_process_tolerance=StocProcTolerances(1e-3, 1e-3), + T=0, + L=SmoothStep(L_op, 0, 2, s) if s is not None else Piecewise([ConstantMatrix(np.zeros((2,2))), ConstantMatrix(L_op)], [0, 1, np.inf]), + H=H, + ) + for s in (None, 0, 1, 2, 3, 4) + ] +#+end_src + +#+RESULTS: + +#+begin_src jupyter-python + for model in models: + plt.plot(model.t, model.L.operator_norm(model.t)) +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/deb60dc0fb8f4c9b582144d5ee00570e39c3f4f7.svg]] + +#+begin_src jupyter-python + fs.tex_value( + models[0].bcf_scale * models[0].bcf(0).real, + prec=1, + prefix="α(0)=", + save="init_slip_alpha", + ), + fs.tex_value(models[0].ω_c, prec=1, prefix="ω_c=", save="initial_slip_cuttoff") +#+end_src + +#+RESULTS: +: \(ω_c=1.0\) + +#+begin_src jupyter-python + aux.integrate_multi(models, 100) +#+end_src + +#+RESULTS: +:RESULTS: +: [INFO hops.core.integration 76514] Choosing the nonlinear integrator. +# [goto error] +#+begin_example + --------------------------------------------------------------------------- + RaySystemError Traceback (most recent call last) + Input In [10], in () + ----> 1 aux.integrate_multi(models, 100) + + File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:84, in integrate_multi(models, *args, **kwargs) +  77 """Integrate the hops equations for the ``models``. +  78 Like :any:`integrate` just for many models. +  79 +  80 A call to :any:`ray.init` may be required. +  81 """ +  83 for model in models: + ---> 84 integrate(model, *args, **kwargs) + + File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:108, in integrate(model, n, data_path, clear_pd) +  98 # with model_db(data_path) as db: +  99 # if hash in db and "data" db[hash] +  101 supervisor = HOPSSupervisor( +  102 model.hops_config, +  103 n, +  104 data_path=data_path, +  105 data_name=hash, +  106 ) + --> 108 supervisor.integrate(clear_pd) +  110 with supervisor.get_data(True) as data: +  111 with model_db(data_path) as db: + + File ~/src/hops/hops/core/integration.py:1238, in HOPSSupervisor.integrate(self, clear_pd) +  1235 with self.get_data_and_maybe_clear(clear_pd) as data: +  1236 t = data.get_time() + -> 1238 num_integrators = int(ray.available_resources().get("CPU", 0)) +  1240 if num_integrators == 0: +  1241 raise RuntimeError("No cpu available for integration!") + + File /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/ray/_private/client_mode_hook.py:105, in client_mode_hook..wrapper(*args, **kwargs) +  103 if func.__name__ != "init" or is_client_mode_enabled_by_default: +  104 return getattr(ray, func.__name__)(*args, **kwargs) + --> 105 return func(*args, **kwargs) + + File /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/ray/state.py:913, in available_resources() +  899 @DeveloperAPI +  900 @client_mode_hook(auto_init=False) +  901 def available_resources(): +  902 """Get the current available cluster resources. +  903 +  904  This is different from `cluster_resources` in that this will return idle +  (...) +  911  resource in the cluster. +  912  """ + --> 913 return state.available_resources() + + File /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/ray/state.py:725, in GlobalState.available_resources(self) +  713 def available_resources(self): +  714 """Get the current available cluster resources. +  715 +  716  This is different from `cluster_resources` in that this will return +  (...) +  723  resource in the cluster. +  724  """ + --> 725 self._check_connected() +  727 available_resources_by_id = self._available_resources_per_node() +  729 # Calculate total available resources. + + File /nix/store/l77qg6h691cglkdhnl4bn06f7flrwkr0-python3-3.9.13-env/lib/python3.9/site-packages/ray/state.py:48, in GlobalState._check_connected(self) +  46 # _really_init_global_state should have set self.global_state_accessor +  47 if self.global_state_accessor is None: + ---> 48 raise ray.exceptions.RaySystemError( +  49 "Ray has not been started yet. You can start Ray with 'ray.init()'." +  50 ) + + RaySystemError: System error: Ray has not been started yet. You can start Ray with 'ray.init()'. +#+end_example +:END: + + + +#+begin_src jupyter-python + %matplotlib inline + for model in models: + _, ax = fs.plot_energy_overview(model) + ax.plot(model.t, model.L.operator_norm(model.t) * np.sqrt(model.δ)) + ax.legend() +#+end_src + +#+RESULTS: +:RESULTS: +# [goto error] +#+begin_example + --------------------------------------------------------------------------- + RuntimeError Traceback (most recent call last) + Input In [27], in () +  1 get_ipython().run_line_magic('matplotlib', 'inline') +  2 for model in models: + ----> 3 _, ax = fs.plot_energy_overview(model) +  4 ax.plot(model.t, model.L.operator_norm(model.t) * np.sqrt(model.δ)) +  5 ax.legend() + + File ~/Documents/Projects/UNI/master/masterarb/python/energy_flow_proper/08_dynamic_one_bath/figsaver.py:549, in plot_energy_overview(model, ensemble_args, **kwargs) +  546 ax.set_ylabel("Energy") +  547 ax.set_xlabel(r"$\tau$") + --> 549 with aux.get_data(model) as data: +  550 system_energy = model.system_energy(data, **ensemble_args) +  551 bath_energy = model.bath_energy(data, **ensemble_args) + + File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:146, in get_data(model, data_path, read_only, **kwargs) +  135 return HIData( +  136 path, +  137 hi_key=model.hops_config, +  (...) +  142 **kwargs, +  143 ) +  145 else: + --> 146 raise RuntimeError(f"No data found for model with hash '{hexhash}'.") + + RuntimeError: No data found for model with hash '07424b6686738dcd1a83299021af44e081e78959f6d428cd5c06fc449984e44a'. +#+end_example +[[file:./.ob-jupyter/5a64809253057fe3638b890f14e8d512e2a730cb.svg]] +:END: + +#+begin_src jupyter-python + fig, ax = plt.subplots() + for model, data in aux.model_data_iterator(models): + flow = model.bath_energy_flow(data) + fs.plot_with_σ( + model.t, + flow, + ax=ax, + label=model.L.order if isinstance(model.L, SmoothStep) else "sudden", + bath=0, + ) + ax.legend() +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/40eea3e28d1225c6d1c991bc0462a806001a265c.svg]] +:END: + +#+begin_src jupyter-python + from hopsflow.util import integrate_array + import scipy + + fig, ax = plt.subplots() + fig.set_size_inches(fs.get_figsize("poster", .5)) + pure_dephasing = [] + pure_dephasing_int = [] + ax.tick_params(axis='y', which='major', pad=10) + + for model in models: + L = model.L + L_fun = L.factor if isinstance(L, SmoothStep) else lambda t: np.heaviside(t - 1, 0) + # plt.plot(model.t, L_deriv(model.t)) + # plt.plot(model.t, L_fun(model.t)) + def α_dot(t): + return -2j * model.ω_c / (1 + 1j * model.ω_c * t) ** 3 + + def α(t): + return 1 / (1 + 1j * model.ω_c * t) ** 2 + + def integrand(s, t): + return np.prod(L_fun(np.array([s, t]))) * α_dot(t - s).imag / 4 + + def integrand_int(s, t): + return np.prod(L_fun(np.array([s, t]))) * α(t - s).imag / 4 + + res = np.array([scipy.integrate.quad(integrand, 0, t, (t,))[0] for t in model.t]) + pure_dephasing.append(res) + pure_dephasing_int.append( + np.array([scipy.integrate.quad(integrand_int, 0, t, (t,))[0] for t in model.t]) + ) + with aux.get_data(model) as data: + _, _, (line, _) = fs.plot_with_σ( + model.t, + abs(-1 * model.bath_energy_flow(data)), + bath=0, + ax=ax, + label=(fr"$s={model.L.order}$") + if isinstance(model.L, SmoothStep) + else "sudden", + ) + ax.plot( + model.t, + abs(-res * model.bcf_scale * model.bcf._c1 * 2), + linestyle="--", + color=line[0].get_color(), + ) + + # fs.plot_with_σ( + # model.t, + # -1 * model.interaction_energy(data), + # linestyle="-.", + # bath=0, + # ax=ax, + # color=line[0].get_color(), + # ) + ax.set_xlim(0.1, 2) + ax.set_ylim(1e-8, 0.3) + ax.set_xscale("log") + ax.set_yscale("log") + + + ax.legend(fontsize="xx-small", loc="upper left") + ax.set_xlabel(r"$\tau$") + ax.set_ylabel(r"$-J$") + + inset = ax.inset_axes([.65, .16, .3, .3]) + for model in models: + inset.plot(model.t, model.L.operator_norm(model.t), linewidth=1) + + inset.set_xlabel(r"$\tau$", labelpad=-20) + inset.set_ylabel(r"$||L||$", labelpad=-20) + fs.export_fig("initial_slip_modcoup") +#+end_src + +#+RESULTS: +:RESULTS: +[[file:./.ob-jupyter/15aeb702b04e742aa50223046bc5c5cfd07b4748.svg]] +[[file:./.ob-jupyter/d982608f086b49fd947ee298a198b0f0a3c1134c.svg]] +:END: + +We see good agreement for short times. + +#+begin_src jupyter-python + fig, ax = plt.subplots() + with aux.get_data(models[-1]) as data: + _, _, (line, _) = fs.plot_with_σ( + model.t, + abs(-1 * model.bath_energy_flow(data)), + bath=0, + ax=ax, + label=model.L.order if isinstance(model.L, SmoothStep) else "sudden", + ) + ax.plot( + model.t, + abs(-pure_dephasing[-1] * model.bcf_scale * model.bcf._c1 * 2), + linestyle="--", + color=line[0].get_color(), + ) + + + + _, _, (line, _) = fs.plot_with_σ( + model.t, + -1 * model.interaction_energy(data), + linestyle="-", + bath=0, + ax=ax, + ) + + ax.plot( + model.t, + abs(-pure_dephasing_int[-1] * model.bcf_scale * model.bcf._c1 * 2), + linestyle="--", + color=line[0].get_color(), + ) + + fs.plot_with_σ( + model.t, + abs(model.system_energy(data) - model.system_energy(data).value[0]), + linestyle="-.", + ax=ax, + ) + ax.plot(model.t, model.L.operator_norm(model.t)) + ax.set_xlim(.7, 2.5) + ax.set_ylim(1e-6) + ax.set_yscale("log") + ax.set_xscale("log") +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/52fc0759fd1eaa929ba61ed90ca9c8ed70b3dbf1.svg]] + +System energy dynamics start only after they diverge. Significant +interaction energy only after that point. Interaction energy is +overestimated (and always wrong) and flow underestimated. -> System energy change. +