clean up 10

.gitignore

@@ -505,3 +505,4 @@ results/**.npy
 /python/energy_flow_proper/11_new_ho_comparison/taurus/
 /python/energy_flow_proper/11_new_ho_comparison/ho_data_local/
 /python/energy_flow_proper/11_new_ho_comparison/local_ho_data/
+/python/energy_flow_proper/10_antizeno_engine/taurus/

python/energy_flow_proper/10_antizeno_engine/10_first_anti_zeno.py

@@ -21,202 +21,14 @@ from hops.util.logging_setup import logging_setup
 import logging
 logging_setup(logging.INFO, show_stocproc=False)
 
-def anti_zeno_engine(
+from anti_zeno_engine import *
-    Δ=5,
-    γ=0.1 / 2,
-    ω_c=1,
-    ω_0=6,
-    ε=1e-3,
-    ε_couple=1 / 2,
-    n=4,
-    cycle_scale=1,
-    switch_cycles=3,
-    detune=0,
-    δ=[10] * 2,
-    T_h=5,
-    therm_initial_state=False,
-    s=[1, 1],
-    dt=0.1,
-    T_c=0,
-    ε_init=0,
-    δ_init=1,
-    L_mat=qt.sigmax().full(),
-    therm_methods=["fft", "fft"],
-    terms=5,
-    sp_tol=1e-5,
-    init_time_steps=None,
-):
-    # τ_bath = 1 / ω_c
-    τ_mod = 2 * np.pi / Δ
-
-    smallest_exponent = (
-        QubitModelMutliBath(δ=[2, 2], ω_c=[ω_c, ω_c])
-        .bcf_coefficients()[1][0]
-        .real.min()
-    )
-
-    τ_off = int(np.around(-np.log(ε) / smallest_exponent / τ_mod)) * τ_mod
-    τ_bath = -np.log(ε_couple) / smallest_exponent
-
-    cycles = int(np.around(τ_bath / τ_mod * cycle_scale))
-    τ_c = cycles * τ_mod
-
-    H_mat = 1 / 2 * (qt.sigmaz().full())
-    Δ_switch = τ_mod * switch_cycles
-    t_max = n * (τ_c + τ_off)
-
-    H = ω_0 * ConstantMatrix(H_mat) + γ * Δ * (
-        Harmonic(H_mat, Δ)
-        @ (
-            Periodic(
-                SmoothStep(np.eye(2), 0, Δ_switch, 2)
-                - SmoothStep(np.eye(2), τ_c - Δ_switch, τ_c, 2),
-                τ_c + τ_off,
-            )
-        )
-    )
-
-    L = Periodic(
-        SmoothStep(L_mat, 0, Δ_switch, 2) - SmoothStep(L_mat, τ_c - Δ_switch, τ_c, 2),
-        τ_c + τ_off,
-    )
-
-    ts = [
-        np.concatenate(
-            [
-                ut.linspace_with_strobe(
-                    max((τ_c + τ_off) * i - Δ_switch, 0),
-                    (τ_c + τ_off) * (i) + τ_c + Δ_switch,
-                    int((τ_c  + 2 * Δ_switch) / dt),
-                    Δ,
-                ),
-                np.linspace((τ_c + τ_off) * (i) + τ_c, (τ_c + τ_off) * (i + 1), 10),
-            ]
-        )
-        for i in range(n)
-    ]
-
-    initializing_period = 0
-    if ε_init > 0:
-        initializing_period = (
-            int(-(np.log(ε_init) / (smallest_exponent * τ_mod) + 1)) * τ_mod
-        )
-        t_max += initializing_period + τ_off
-        H = Piecewise(
-            [
-                ω_0 * ConstantMatrix(H_mat)
-                + Harmonic(H_mat, Δ)
-                * γ
-                * Δ
-                @ (
-                    ConstantMatrix(np.eye(2))
-                    - SmoothStep(
-                        np.eye(2),
-                        initializing_period - Δ_switch,
-                        initializing_period,
-                        2,
-                    )
-                ),
-                ω_0 * ConstantMatrix(H_mat),
-                Shift(
-                    H,
-                    initializing_period + τ_off,
-                ),
-            ],
-            [0, initializing_period, initializing_period + τ_off, np.inf],
-        )
-
-        L = Piecewise(
-            [
-                (
-                    SmoothStep(
-                        L_mat,
-                        0,
-                        Δ_switch,
-                        2,
-                    )
-                    - SmoothStep(
-                        L_mat,
-                        initializing_period - Δ_switch,
-                        initializing_period,
-                        2,
-                    )
-                )
-                * np.sqrt(δ_init),
-                ConstantMatrix(np.zeros((2, 2))),
-                Shift(
-                    L,
-                    initializing_period + τ_off,
-                ),
-            ],
-            [0, initializing_period, initializing_period + τ_off, np.inf],
-        )
-
-        ts = [times + initializing_period + τ_off for times in ts]
-        ts.insert(
-            0,
-            np.concatenate(
-                [
-                    ut.linspace_with_strobe(
-                        0,
-                        initializing_period + Δ_switch,
-                        init_time_steps or int((initializing_period + Δ_switch) / (dt)),
-                        Δ,
-                    ),
-                    np.linspace(initializing_period, initializing_period + τ_off, 10),
-                ]
-            ),
-        )
-
-    model = QubitModelMutliBath(
-        δ=δ,
-        ω_c=[ω_c, ω_c],
-        ω_s=[ω_0 - Δ - ω_c + detune, ω_0 + Δ - ω_c + detune],
-        t=np.unique(np.concatenate(ts)),
-        ψ_0=(
-            np.sqrt(1 / (np.exp(ω_0 / T_h) + 1)) * qt.basis([2], [0])
-            + np.sqrt(1 / (np.exp(-ω_0 / T_h) + 1)) * qt.basis([2], [1])
-            if therm_initial_state
-            else qt.basis([2], [1])
-        ),
-        description=f"Building the anti-zento-engine",
-        k_max=4,
-        bcf_terms=[terms, terms],
-        truncation_scheme="simplex",
-        driving_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
-        thermal_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
-        T=[T_c, T_h],
-        L=[L] * 2,
-        H=H,
-        therm_methods=therm_methods,
-        s=s,
-    )
-
-    return (
-        model,
-        Δ,
-        (
-            τ_mod,
-            τ_c,
-            τ_bath,
-            cycles,
-            model.ω_s,
-            ω_0,
-            τ_mod,
-            τ_off,
-            n,
-            Δ_switch,
-            initializing_period,
-        ),
-    )
 
 (
     model,
-    Δ,
+    params
-    (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s, τ_off, n, Δ_switch, τ_init),
 ) = anti_zeno_engine(
     Δ=11,
-    ε=.1,#.1,
+    ε=1,#.1,
     ω_c=1,
     ε_couple=0.7,
     n=10,
@@ -235,264 +47,15 @@ def anti_zeno_engine(
     init_time_steps=10,
 )
 model.k_max = 4
-# model, params = anti_zeno_engine(ε=1/2, ε_couple=1e-4, n=1, detune=.5, δ=[.1,.1])
 
-fs.tex_value(
+params.cycles
-    model.bcf_scales[0] * model.bcf(0)(0).real,
-    prec=1,
-    prefix="α(0)=",
-    save="omega_alpha",
-), fs.tex_value(
-    ω_0,
-    prec=0,
-    prefix="ω_0=",
-    save="omega_zero",
-), fs.tex_value(
-    5 / 10,
-    prec=1,
-    prefix="γ=",
-    save="gamma",
-),fs.tex_value(
-    Δ,
-    prec=0,
-    prefix="Δ=",
-    save="delta",
-),fs.tex_value(
-    model.T[0],
-    prec=0,
-    prefix="T_c=",
-    save="tc",
-),fs.tex_value(
-    model.T[1],
-    prec=0,
-    prefix="T_h=",
-    save="th",
-)
 
-cycles
+plot_az_coupling_diagram(model, params)
 
-ts = np.linspace(0,50,1000)
+plot_az_sd_overview(model, params)
-fig, ax = pu.plot_complex(ts, model.bcf(0)(ts))
-pu.plot_complex(ts, model.thermal_correlations(0)(ts))
 
-ts = np.linspace(0,10,1000)
+plot_total_power(model, params)
-proc = model.thermal_process(1)
-import hops
-z=hops.core.utility.uni_to_gauss(np.random.rand(proc.get_num_y() * 2))
-proc.new_process(z)
-pu.plot_complex(ts, proc(ts) * model.bcf_scales[1])
 
-G_h = model.thermal_spectral_density(1)
+plot_excited_state(model, params)
-G_c = model.spectral_density(0)
-def chi(v, ω_0):
-    return G_h(ω_0 + v) - G_c(ω_0 - v)
 
-vs = np.linspace(0.1, 10, 100)
+plot_power_output(model, params)
-plt.plot(vs, chi(vs, ω_0))
-plt.plot(vs, G_h(vs))
-
-aux.integrate(model, 10_000)
-
-#_, ax = pu.plot_energy_overview(model, markersize=1, ensemble_args=dict(gc_sleep=0.05))
-
-fig, ax = plt.subplots()
-with aux.get_data(model) as data:
-    pu.plot_with_σ(model.t, model.total_power(data, gc_sleep=0), ax=ax)
-    #pu.plot_with_σ(model.t, model.total_energy_from_power(data), ax=ax)
-
-#ax.legend()
-
-with aux.get_data(model) as data:
-    pu.plot_with_σ(
-        model.t,
-        abs(EnsembleValue(
-            (data.rho_t_accum.mean[:, 0, 0], data.rho_t_accum.ensemble_std[:, 0, 0])
-        ) - .5),
-    )
-    plt.yscale("log")
-    ρ_ee = data.rho_t_accum.mean[:, 0, 0].real
-
-plt.plot(model.t[101:], ut.smoothen(model.t[101:], ρ_ee[101:], frac=.5, it=0))
-
-with aux.get_data(model) as data:
-  _, ax = plt.subplots()
-  pu.plot_with_σ(model.t, model.bath_energy(data), bath=0, ax=ax)
-  pu.plot_with_σ(model.t, model.bath_energy(data), bath=1, ax=ax)
-  pu.plot_with_σ(model.t, model.bath_energy(data).sum_baths(), ax=ax)
-  pu.plot_with_σ(model.t, model.total_energy_from_power(data), ax=ax)
-
-  pu.plot_with_σ(model.t, model.total_energy(data), ax=ax)
-  #pu.plot_with_σ(model.t, model.interaction_energy(data).for_bath(1), ax=ax)
-  #pu.plot_with_σ(model.t, model.system_energy(data), ax=ax)
-  #pu.plot_with_σ(model.t, model.system_energy(data) + model.bath_energy(data).sum_baths(), ax=ax)
-  #ax.plot(model.t, ut.smoothen(model.t, model.total_energy(data).value, frac=.1))
-  #ax.plot(model.t, np.gradient(model.total_energy(data).value))
-
-ts_begin =  τ_init + τ_off + ((τ_c + τ_off) * np.arange(0, n))
-ts_end = τ_c + ts_begin
-ind_begin = np.searchsorted(model.t, ts_begin)
-ind_end = np.searchsorted(model.t, ts_end)
-ts_begin = model.t[ind_begin]
-ts_end = model.t[ind_end]
-
-# #plt.plot(model.t, model.L[0].operator_norm(model.t))
-
-with aux.get_data(model) as data:
-    tot_power = (model.total_power(data).sum_baths()) 
-    #model.total_energy(data).value
-
-fig, ax = plt.subplots()
-
-for t in ts_begin:
-    plt.axvline(t, linestyle="dashed", color="orangered", linewidth=1)
-
-for t in ts_end:
-    plt.axvline(t, linestyle="dotted", color="lightblue", linewidth=1)
-
-pu.plot_with_σ(
-    model.t[ind_begin[0]:], tot_power.slice(slice(ind_begin[0], None, 1)), ax=ax, linewidth=0.5
-)
-ax.plot(
-    model.t[ind_begin[0]:],
-    ut.smoothen(model.t[ind_begin[0] :], tot_power.value[ind_begin[0] :], frac=0.05),
-)
-
-powers = []
-for begin, end in zip(ind_begin, ind_end):
-    powers.append(
-        (tot_power.slice(slice(begin, end))).mean
-    )
-
-steady_index = 0
-power = sum(powers[steady_index:]) * (1/(len(powers[steady_index:])))
-
-ax.set_ylabel(r"$\langle{P}\rangle$")
-ax.set_xlabel(r"$\tau$")
-
-
-fig.suptitle(fr"$\omega_0={ω_0},\,\omega_c={1},\,N={tot_power.N},\,\delta={.5}\,\Delta={Δ},\,\lambda={.2},\, T_h={model.T[0]},\, T_c={model.T[1]},\,n={cycles}$")
-fs.export_fig("anti_zeno_maybe_note_quite_steady", tikz=False)
-fs.tex_value(
-    (power.value),
-    err=power.σ,
-    prefix="P=",
-    save="power_with_cool",
-)
-
-mean_norm = np.zeros_like(model.t)
-with aux.get_data(model) as data:
-    for i in range(min(data.samples, 2)):
-        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[40:], mean_norm[40:] / mean_norm.max())
-plt.yscale("log")
-mean_norm.max()
-
-model, _ = anti_zeno_engine(Δ=5, γ=0.1 / 2, ω_c=1, ω_0=2, ε=1e-1, ε_couple=1 / 3, n=6, cycle_scale=1, switch_cycles=3)
-
-_, ax = pu.plot_energy_overview(model, markersize=1)
-
-model, params = anti_zeno_engine(ω_0=6)
-
-_, ax = pu.plot_energy_overview(model, markersize=1)
-
-model, params = anti_zeno_engine(ω_0=6, ε=1e-1)
-
-_, ax = pu.plot_energy_overview(model, markersize=1)
-
-model, params = anti_zeno_engine(ε=1e-1, ε_couple=1/4)
-
-model, params = anti_zeno_engine(ε=1e-1, ε_couple=1e-3, n=3)
-
-model, Δ, params = anti_zeno_engine(ε=1/2, ε_couple=1e-10, n=1, detune=5, ω_0 = 10)
-
-model, Δ, params = anti_zeno_engine(ε=1/5, ε_couple=1/4, n=5, detune=5, ω_0 = 10)
-#model, params = anti_zeno_engine(ε=1/2, ε_couple=1e-4, n=1, detune=.5, δ=[.1,.1])
-params
-
-model, Δ, params = anti_zeno_engine(
-    ε=1 / 10,
-    ε_couple=1 / 100,
-    n=1,
-    detune=4,
-    ω_0=10,
-    δ=[2] * 2,
-    γ=0.05 / 2,
-    therm_initial_state=True,
-)
-
-# model, Δ, (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s) = anti_zeno_engine(Δ = 12, ε=1/2, ω_c=1, ε_couple=.9, n=10, detune=-1, ω_0 = 15, T_h=20, δ=[6]*2, γ=.2, therm_initial_state=False)
-model, Δ, (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s) = anti_zeno_engine(
-    Δ=1.5,
-    ε=.4,
-    ω_c=0.2,
-    ε_couple=0.5,
-    n=2,
-    detune=-0.2,
-    ω_0=2,
-    T_h=6,
-    δ=[0.1] * 2,
-    γ=0.5/10,
-    switch_cycles=4,
-    therm_initial_state=False,
-    s=[1.2]*2
-)
-model.k_max = 4
-# model, params = anti_zeno_engine(ε=1/2, ε_couple=1e-4, n=1, detune=.5, δ=[.1,.1])
-
-model, Δ, (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s, τ_off, n, Δ_switch) = anti_zeno_engine(
-      Δ=1,
-      ε=.7 ,
-      ω_c=0.1,
-      ε_couple=0.9,
-      n=7,
-      detune=-0.08,
-      ω_0=2,
-      T_h=4,
-      δ=[1] * 2,
-      γ=0.5 / 10,
-      switch_cycles=1,
-      therm_initial_state=False,
-      s=[1] * 2,
-  )
-
-model, Δ, (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s, τ_off, n, Δ_switch) = anti_zeno_engine(
-    Δ=1,
-    ε=.94 ,
-    ω_c=0.1,
-    ε_couple=0.9,
-    n=7,
-    detune=-0.08,
-    ω_0=2,
-    T_h=4,
-    δ=[1] * 2,
-    γ=0.5 / 10,
-    switch_cycles=1,
-    therm_initial_state=False,
-    s=[1] * 2,
-)
-
-model, Δ, (τ_mod, τ_c, τ_bath, cycles, model.ω_s, ω_0, τ_s, τ_off, n, Δ_switch) = anti_zeno_engine(
-    Δ=1.3,
-    ε=.5,
-    ω_c=0.2,
-    ε_couple=0.7,
-    n=4,
-    detune=-0.19,
-    ω_0=2,
-    T_c=3,
-    T_h=11,
-    δ=[1] * 2,
-    γ=0.5 / 10,
-    switch_cycles=1,
-    therm_initial_state=False,
-    s=[1] * 2,
-    ε_init=.6,
-    δ_init=2
-)
-model.k_max = 4

python/energy_flow_proper/10_antizeno_engine/10_first_anti_zeno_cc.py

@@ -0,0 +1,61 @@
+import figsaver as fs
+from hiro_models.one_qubit_model import QubitModelMutliBath, 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
+from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix, Piecewise, Shift
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib
+import plot_utils as pu
+
+import ray
+ray.shutdown()
+ray.init()
+
+from hops.util.logging_setup import logging_setup
+import logging
+logging_setup(logging.INFO, show_stocproc=False)
+
+from anti_zeno_engine import *
+
+(
+    model,
+    params
+) = anti_zeno_engine(
+    Δ=11,
+    ε=1,#.1,
+    ω_c=1,
+    ε_couple=0.7,
+    n=10,
+    detune=.5,
+    ω_0=20,
+    T_c=1e3,
+    T_h=1e4,
+    δ=[3.2*.01/4, 1*.01/4],
+    γ=.2,
+    switch_cycles=1,
+    therm_initial_state=False,
+    ε_init=.00001/2,
+    terms=6,
+    dt=0.01,
+    sp_tol=1e-3,
+    init_time_steps=10,
+)
+model.k_max = 4
+
+params.cycles
+
+plot_az_coupling_diagram(model, params)
+
+plot_az_sd_overview(model, params)
+
+plot_total_power(model, params)
+
+plot_excited_state(model, params)
+
+plot_power_output(model, params)

python/energy_flow_proper/10_antizeno_engine/anti_zeno_engine.org

@@ -37,194 +37,7 @@ Init ray and silence stocproc.
 
 * Model Definition
 #+begin_src jupyter-python :results none
-  def anti_zeno_engine(
+  from anti_zeno_engine import *
-      Δ=5,
-      γ=0.1 / 2,
-      ω_c=1,
-      ω_0=6,
-      ε=1e-3,
-      ε_couple=1 / 2,
-      n=4,
-      cycle_scale=1,
-      switch_cycles=3,
-      detune=0,
-      δ=[10] * 2,
-      T_h=5,
-      therm_initial_state=False,
-      s=[1, 1],
-      dt=0.1,
-      T_c=0,
-      ε_init=0,
-      δ_init=1,
-      L_mat=qt.sigmax().full(),
-      therm_methods=["fft", "fft"],
-      terms=5,
-      sp_tol=1e-5,
-      init_time_steps=None,
-  ):
-      # τ_bath = 1 / ω_c
-      τ_mod = 2 * np.pi / Δ
-
-      smallest_exponent = (
-          QubitModelMutliBath(δ=[2, 2], ω_c=[ω_c, ω_c])
-          .bcf_coefficients()[1][0]
-          .real.min()
-      )
-
-      τ_off = int(np.around(-np.log(ε) / smallest_exponent / τ_mod)) * τ_mod
-      τ_bath = -np.log(ε_couple) / smallest_exponent
-
-      cycles = int(np.around(τ_bath / τ_mod * cycle_scale))
-      τ_c = cycles * τ_mod
-
-      H_mat = 1 / 2 * (qt.sigmaz().full())
-      Δ_switch = τ_mod * switch_cycles
-      t_max = n * (τ_c + τ_off)
-
-      H = ω_0 * ConstantMatrix(H_mat) + γ * Δ * (
-          Harmonic(H_mat, Δ)
-          @ (
-              Periodic(
-                  SmoothStep(np.eye(2), 0, Δ_switch, 2)
-                  - SmoothStep(np.eye(2), τ_c - Δ_switch, τ_c, 2),
-                  τ_c + τ_off,
-              )
-          )
-      )
-
-      L = Periodic(
-          SmoothStep(L_mat, 0, Δ_switch, 2) - SmoothStep(L_mat, τ_c - Δ_switch, τ_c, 2),
-          τ_c + τ_off,
-      )
-
-      ts = [
-          np.concatenate(
-              [
-                  ut.linspace_with_strobe(
-                      max((τ_c + τ_off) * i - Δ_switch, 0),
-                      (τ_c + τ_off) * (i) + τ_c + Δ_switch,
-                      int((τ_c  + 2 * Δ_switch) / dt),
-                      Δ,
-                  ),
-                  np.linspace((τ_c + τ_off) * (i) + τ_c, (τ_c + τ_off) * (i + 1), 10),
-              ]
-          )
-          for i in range(n)
-      ]
-
-      initializing_period = 0
-      if ε_init > 0:
-          initializing_period = (
-              int(-(np.log(ε_init) / (smallest_exponent * τ_mod) + 1)) * τ_mod
-          )
-          t_max += initializing_period + τ_off
-          H = Piecewise(
-              [
-                  ω_0 * ConstantMatrix(H_mat)
-                  + Harmonic(H_mat, Δ)
-                  ,* γ
-                  ,* Δ
-                  @ (
-                      ConstantMatrix(np.eye(2))
-                      - SmoothStep(
-                          np.eye(2),
-                          initializing_period - Δ_switch,
-                          initializing_period,
-                          2,
-                      )
-                  ),
-                  ω_0 * ConstantMatrix(H_mat),
-                  Shift(
-                      H,
-                      initializing_period + τ_off,
-                  ),
-              ],
-              [0, initializing_period, initializing_period + τ_off, np.inf],
-          )
-
-          L = Piecewise(
-              [
-                  (
-                      SmoothStep(
-                          L_mat,
-                          0,
-                          Δ_switch,
-                          2,
-                      )
-                      - SmoothStep(
-                          L_mat,
-                          initializing_period - Δ_switch,
-                          initializing_period,
-                          2,
-                      )
-                  )
-                  ,* np.sqrt(δ_init),
-                  ConstantMatrix(np.zeros((2, 2))),
-                  Shift(
-                      L,
-                      initializing_period + τ_off,
-                  ),
-              ],
-              [0, initializing_period, initializing_period + τ_off, np.inf],
-          )
-
-          ts = [times + initializing_period + τ_off for times in ts]
-          ts.insert(
-              0,
-              np.concatenate(
-                  [
-                      ut.linspace_with_strobe(
-                          0,
-                          initializing_period + Δ_switch,
-                          init_time_steps or int((initializing_period + Δ_switch) / (dt)),
-                          Δ,
-                      ),
-                      np.linspace(initializing_period, initializing_period + τ_off, 10),
-                  ]
-              ),
-          )
-
-      model = QubitModelMutliBath(
-          δ=δ,
-          ω_c=[ω_c, ω_c],
-          ω_s=[ω_0 - Δ - ω_c + detune, ω_0 + Δ - ω_c + detune],
-          t=np.unique(np.concatenate(ts)),
-          ψ_0=(
-              np.sqrt(1 / (np.exp(ω_0 / T_h) + 1)) * qt.basis([2], [0])
-              + np.sqrt(1 / (np.exp(-ω_0 / T_h) + 1)) * qt.basis([2], [1])
-              if therm_initial_state
-              else qt.basis([2], [1])
-          ),
-          description=f"Building the anti-zento-engine",
-          k_max=4,
-          bcf_terms=[terms, terms],
-          truncation_scheme="simplex",
-          driving_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
-          thermal_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
-          T=[T_c, T_h],
-          L=[L] * 2,
-          H=H,
-          therm_methods=therm_methods,
-          s=s,
-      )
-
-      return (
-          model,
-          Δ,
-          (
-              τ_mod,
-              τ_c,
-              τ_bath,
-              cycles,
-              model.ω_s,
-              ω_0,
-              τ_mod,
-              τ_off,
-              n,
-              Δ_switch,
-              initializing_period,
-          ),
-      )
 #+end_src
 
 * Model Definition
@@ -269,7 +82,7 @@ Init ray and silence stocproc.
   ), fs.tex_value(
       5 / 10,
       prec=1,
-      prefix="γ=",
+      prefix="λ=",
       save="gamma",
   ),fs.tex_value(
       Δ,

python/energy_flow_proper/10_antizeno_engine/anti_zeno_engine.py

@@ -0,0 +1,371 @@
+from hiro_models.one_qubit_model import QubitModelMutliBath, StocProcTolerances
+import numpy as np
+import qutip as qt
+import utilities as ut
+from hops.util.dynamic_matrix import (
+    SmoothStep,
+    Periodic,
+    Harmonic,
+    ConstantMatrix,
+    Piecewise,
+    Shift,
+)
+import numpy as np
+from types import SimpleNamespace
+from plot_utils import wrap_plot, plot_with_σ
+import figsaver as fs
+import matplotlib.pyplot as plt
+import hiro_models.model_auxiliary as aux
+from hopsflow.util import EnsembleValue
+
+
+def anti_zeno_engine(
+    Δ=5,
+    γ=0.1 / 2,
+    ω_c=1,
+    ω_0=6,
+    ε=1e-3,
+    ε_couple=1 / 2,
+    n=4,
+    cycle_scale=1,
+    switch_cycles=3,
+    detune=0,
+    δ=[10] * 2,
+    T_h=5,
+    therm_initial_state=False,
+    s=[1, 1],
+    dt=0.1,
+    T_c=0,
+    ε_init=0,
+    δ_init=1,
+    L_mat=qt.sigmax().full(),
+    therm_methods=["fft", "fft"],
+    terms=5,
+    sp_tol=1e-5,
+    init_time_steps=None,
+):
+    # τ_bath = 1 / ω_c
+    τ_mod = 2 * np.pi / Δ
+
+    smallest_exponent = (
+        QubitModelMutliBath(δ=[2, 2], ω_c=[ω_c, ω_c])
+        .bcf_coefficients()[1][0]
+        .real.min()
+    )
+
+    τ_off = int(np.around(-np.log(ε) / smallest_exponent / τ_mod)) * τ_mod
+    τ_bath = -np.log(ε_couple) / smallest_exponent
+
+    cycles = int(np.around(τ_bath / τ_mod * cycle_scale))
+    τ_c = cycles * τ_mod
+
+    H_mat = 1 / 2 * (qt.sigmaz().full())
+    Δ_switch = τ_mod * switch_cycles
+    t_max = n * (τ_c + τ_off)
+
+    H = ω_0 * ConstantMatrix(H_mat) + γ * Δ * (
+        Harmonic(H_mat, Δ)
+        @ (
+            Periodic(
+                SmoothStep(np.eye(2), 0, Δ_switch, 2)
+                - SmoothStep(np.eye(2), τ_c - Δ_switch, τ_c, 2),
+                τ_c + τ_off,
+            )
+        )
+    )
+
+    L = Periodic(
+        SmoothStep(L_mat, 0, Δ_switch, 2) - SmoothStep(L_mat, τ_c - Δ_switch, τ_c, 2),
+        τ_c + τ_off,
+    )
+
+    ts = [
+        np.concatenate(
+            [
+                ut.linspace_with_strobe(
+                    max((τ_c + τ_off) * i - Δ_switch, 0),
+                    (τ_c + τ_off) * (i) + τ_c + Δ_switch,
+                    int((τ_c + 2 * Δ_switch) / dt),
+                    Δ,
+                ),
+                np.linspace((τ_c + τ_off) * (i) + τ_c, (τ_c + τ_off) * (i + 1), 10),
+            ]
+        )
+        for i in range(n)
+    ]
+
+    initializing_period = 0
+    if ε_init > 0:
+        initializing_period = (
+            int(-(np.log(ε_init) / (smallest_exponent * τ_mod) + 1)) * τ_mod
+        )
+        t_max += initializing_period + τ_off
+        H = Piecewise(
+            [
+                ω_0 * ConstantMatrix(H_mat)
+                + Harmonic(H_mat, Δ)
+                * γ
+                * Δ
+                @ (
+                    ConstantMatrix(np.eye(2))
+                    - SmoothStep(
+                        np.eye(2),
+                        initializing_period - Δ_switch,
+                        initializing_period,
+                        2,
+                    )
+                ),
+                ω_0 * ConstantMatrix(H_mat),
+                Shift(
+                    H,
+                    initializing_period + τ_off,
+                ),
+            ],
+            [0, initializing_period, initializing_period + τ_off, np.inf],
+        )
+
+        L = Piecewise(
+            [
+                (
+                    SmoothStep(
+                        L_mat,
+                        0,
+                        Δ_switch,
+                        2,
+                    )
+                    - SmoothStep(
+                        L_mat,
+                        initializing_period - Δ_switch,
+                        initializing_period,
+                        2,
+                    )
+                )
+                * np.sqrt(δ_init),
+                ConstantMatrix(np.zeros((2, 2))),
+                Shift(
+                    L,
+                    initializing_period + τ_off,
+                ),
+            ],
+            [0, initializing_period, initializing_period + τ_off, np.inf],
+        )
+
+        ts = [times + initializing_period + τ_off for times in ts]
+        ts.insert(
+            0,
+            np.concatenate(
+                [
+                    ut.linspace_with_strobe(
+                        0,
+                        initializing_period + Δ_switch,
+                        init_time_steps or int((initializing_period + Δ_switch) / (dt)),
+                        Δ,
+                    ),
+                    np.linspace(initializing_period, initializing_period + τ_off, 10),
+                ]
+            ),
+        )
+
+    model = QubitModelMutliBath(
+        δ=δ,
+        ω_c=[ω_c, ω_c],
+        ω_s=[ω_0 - Δ - ω_c + detune, ω_0 + Δ - ω_c + detune],
+        t=np.unique(np.concatenate(ts)),
+        ψ_0=(
+            np.sqrt(1 / (np.exp(ω_0 / T_h) + 1)) * qt.basis([2], [0])
+            + np.sqrt(1 / (np.exp(-ω_0 / T_h) + 1)) * qt.basis([2], [1])
+            if therm_initial_state
+            else qt.basis([2], [1])
+        ),
+        description=f"Building the anti-zento-engine",
+        k_max=4,
+        bcf_terms=[terms, terms],
+        truncation_scheme="simplex",
+        driving_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
+        thermal_process_tolerances=[StocProcTolerances(sp_tol, sp_tol)] * 2,
+        T=[T_c, T_h],
+        L=[L] * 2,
+        H=H,
+        therm_methods=therm_methods,
+        s=s,
+    )
+
+    model_params = SimpleNamespace(
+        Δ=Δ,
+        τ_mod=τ_mod,
+        τ_c=τ_c,
+        τ_bath=τ_bath,
+        cycles=cycles,
+        ω_0=ω_0,
+        τ_off=τ_off,
+        ω_c=ω_c,
+        n=n,
+        Δ_switch=Δ_switch,
+        initializing_period=initializing_period,
+        τ_init=initializing_period,
+        τ_s=τ_mod,
+        λ=γ,
+        δ=detune,
+    )
+
+    return (model, model_params)
+
+
+@wrap_plot
+def plot_az_coupling_diagram(model, params, ax=None):
+    def thermal_bcf(t):
+        return model.bcf(1)(t) + 2 * (model.thermal_correlations(1)(t).real)
+
+    # plt.plot(model.t, np.abs(thermal_bcf(model.t))/np.abs(thermal_bcf(0)), alpha=.5, zorder=5, label="BCF (cold bath)")
+    # plt.plot(model.t, abs(model.bcf(0)(model.t)) / abs(model.bcf(0)(model.t)).max(), linewidth=1)
+
+    t = np.linspace(0, model.t.max(), 10000)
+    ax.plot(t, model.L[0].operator_norm(t), label=r"$||L||$")
+    ax.plot(
+        t,
+        (model.H.operator_norm(t) - params.ω_0 / 2) / params.ω_0,
+        label=r"$||H_{\mathrm{S}}||$",
+    )
+    ax.plot(
+        t,
+        np.exp(-t * np.min(np.array(model.bcf_coefficients()[1][0]).real)),
+        label="Slowest Decaying BCF Term",
+    )
+    ax.set_xlabel(r"$\tau$")
+    ax.legend()
+
+
+@wrap_plot
+def plot_az_sd_overview(model, params, ax=None):
+    ωs = np.linspace(0.01 + model.ω_s[0] - 1, 2 * params.ω_0, 10000)
+
+    def total_sd(n, ω):
+        return model.bcf_scales[n] * (
+            model.spectral_density(n)(ω)
+            * (1 / (np.expm1(ωs / model.T[n])) + np.heaviside(ω, 0))
+        )
+
+    ax.plot(
+        ωs,
+        total_sd(0, ωs),
+        label="Cold Bath",
+    )
+    ax.plot(
+        ωs,
+        total_sd(1, ωs),
+        label="Hot Bath",
+    )
+    ax.axvline(model.H.operator_norm(0) * 2 + params.Δ)
+    ax.axvline(model.H.operator_norm(0) * 2 - params.Δ)
+    for mult in [1, 2, 3]:
+        ax.plot(
+            ωs,
+            total_sd(1, ωs).max()
+            * np.sinc((ωs - params.ω_0 - params.Δ) * params.τ_c * mult),
+            label=rf"Filter $n={params.cycles * mult}$",
+        )
+
+    ax.set_xlabel(r"$\omega$")
+    ax.legend()
+
+
+@wrap_plot
+def plot_total_power(model, params, ax=None):
+    with aux.get_data(model) as data:
+        mask = model.t > params.τ_init
+        plot_with_σ(
+            model.t[mask], model.total_power(data, gc_sleep=0).slice(mask), ax=ax
+        )
+
+
+@wrap_plot
+def plot_excited_state(model, params, ax=None):
+    with aux.get_data(model) as data:
+        mask = model.t > params.τ_init
+
+        plot_with_σ(
+            model.t[mask],
+            abs(
+                EnsembleValue(
+                    (
+                        data.rho_t_accum.mean[mask, 0, 0],
+                        data.rho_t_accum.ensemble_std[mask, 0, 0],
+                    )
+                )
+            ),
+            ax=ax,
+        )
+
+
+def plot_power_output(model, params):
+    fig, ax = plt.subplots()
+    ts_begin = (
+        params.τ_init
+        + params.τ_off
+        + ((params.τ_c + params.τ_off) * np.arange(0, params.n))
+    )
+    ts_end = params.τ_c + ts_begin
+    ind_begin = np.searchsorted(model.t, ts_begin)
+    ind_end = np.searchsorted(model.t, ts_end)
+    ts_begin = model.t[ind_begin]
+    ts_end = model.t[ind_end]
+
+    # #plt.plot(model.t, model.L[0].operator_norm(model.t))
+
+    with aux.get_data(model) as data:
+        tot_power = model.total_power(data, gc_sleep=0).sum_baths()
+
+    for t in ts_begin:
+        plt.axvline(t, linestyle="dashed", color="orangered", linewidth=1)
+
+    for t in ts_end:
+        plt.axvline(t, linestyle="dotted", color="lightblue", linewidth=1)
+
+    plot_with_σ(
+        model.t[ind_begin[0] :],
+        tot_power.slice(slice(ind_begin[0], None, 1)),
+        ax=ax,
+        linewidth=0.5,
+    )
+    # ax.plot(
+    #     model.t[ind_begin[0] :],
+    #     ut.smoothen(
+    #         model.t[ind_begin[0] :], tot_power.value[ind_begin[0] :], frac=0.05
+    #     ),
+    # )
+
+    powers = []
+    for begin, end in zip(ind_begin, ind_end):
+        part_power = (tot_power.slice(slice(begin, end))).mean
+        powers.append(part_power)
+
+        plot_with_σ(
+            model.t[[begin, end]],
+            EnsembleValue(
+                (
+                    np.array([part_power.value] * 2),
+                    np.array([part_power.σ] * 2),
+                )
+            ),
+            color="black",
+            ax=ax,
+        )
+
+    steady_index = 0
+    power = sum(powers[steady_index:]) * (1 / (len(powers[steady_index:])))
+
+    ax.set_ylabel(r"$\langle{P}\rangle$")
+    ax.set_xlabel(r"$\tau$")
+
+    fig.suptitle(
+        fr"$\omega_0={params.ω_0},\,\omega_c={params.ω_c},\,N={tot_power.N},\,\delta={params.δ},\,\Delta={params.Δ},\,\lambda={params.λ},\, T_h={model.T[0]},\, T_c={model.T[1]},\,n={params.cycles}$"
+    )
+    ax.set_title(
+        fs.tex_value(
+            (power.value),
+            err=power.σ,
+            prefix=r"\bar{P}=",
+        )
+    )
+
+    return fig, ax, power

python/energy_flow_proper/10_antizeno_engine/continously_couple.org

@@ -0,0 +1,130 @@
+#+PROPERTY: header-args :session 10_anti_zeno_engine_cc :kernel python :pandoc no :async yes :tangle 10_first_anti_zeno_cc.py
+Here we try to reproduce the anti zeno engine from the paper.
+
+* Boilerplate
+#+begin_src jupyter-python :results none
+  import figsaver as fs
+  from hiro_models.one_qubit_model import QubitModelMutliBath, 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
+  from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix, Piecewise, Shift
+  import matplotlib.pyplot as plt
+  import numpy as np
+  import matplotlib
+  import plot_utils as pu
+#+end_src
+
+Init ray and silence stocproc.
+#+begin_src jupyter-python
+  import ray
+  ray.shutdown()
+  ray.init()
+#+end_src
+
+#+RESULTS:
+: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '141.30.17.225', 'raylet_ip_address': '141.30.17.225', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-12_20-07-54_888036_22466/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-12_20-07-54_888036_22466/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-12_20-07-54_888036_22466', 'metrics_export_port': 60638, 'gcs_address': '141.30.17.225:50826', 'address': '141.30.17.225:50826', 'node_id': 'a3660dfc275e05eee6674ca8f0ba5c4fc2754aa9c49703a883bebc9e'})
+
+#+begin_src jupyter-python :results none
+  from hops.util.logging_setup import logging_setup
+  import logging
+  logging_setup(logging.INFO, show_stocproc=False)
+#+end_src
+
+* Model Definition
+#+begin_src jupyter-python :results none
+  from anti_zeno_engine import *
+#+end_src
+
+#+begin_src jupyter-python
+  (
+      model,
+      params
+  ) = anti_zeno_engine(
+      Δ=11,
+      ε=1,#.1,
+      ω_c=1,
+      ε_couple=0.7,
+      n=10,
+      detune=.5,
+      ω_0=20,
+      T_c=1e3,
+      T_h=1e4,
+      δ=[3.2*.01/4, 1*.01/4],
+      γ=.2,
+      switch_cycles=1,
+      therm_initial_state=False,
+      ε_init=.00001/2,
+      terms=6,
+      dt=0.01,
+      sp_tol=1e-3,
+      init_time_steps=10,
+  )
+  model.k_max = 4
+#+end_src
+
+#+RESULTS:
+
+Let's test the assumptions of the paper.
+
+#+begin_src jupyter-python
+  params.cycles
+#+end_src
+
+#+RESULTS:
+: 9
+
+#+begin_src jupyter-python
+  plot_az_coupling_diagram(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:xlabel= | $\tau$ | > |
+[[file:./.ob-jupyter/a941e4ac84ce34fd566734c1aba25daec397b27d.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_az_sd_overview(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:xlabel= | $\omega$ | > |
+[[file:./.ob-jupyter/997e7a1fda9ba22feba8ddba7eae9bde652f463f.svg]]
+:END:
+
+
+#+begin_src jupyter-python
+  plot_total_power(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:> |
+[[file:./.ob-jupyter/a3b4954409608a3f24b4cfc3ffd2462e0584aad2.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_excited_state(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:> |
+[[file:./.ob-jupyter/e80d88373632fdca29311847c2cd654fec89c2bf.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_power_output(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:title= | (center : \(\bar{P}=-0.00044\pm 0.00020\)) | xlabel= | $\tau$ | ylabel= | $\langle{P}\rangle$ | > | EnsembleValue | (((100000 -0.0004421216901391656 0.00020060870097189892))) |
+[[file:./.ob-jupyter/5f54c22a0a8a8ef9cef5cb03b167e0152b4504dc.svg]]
+:END:

python/energy_flow_proper/10_antizeno_engine/nil

@@ -1,19 +1,56 @@
 def thermal_bcf(t):
     return model.bcf(1)(t) + 2 * (model.thermal_correlations(1)(t).real)
-plt.plot(model.t, np.abs(thermal_bcf(model.t))/np.abs(thermal_bcf(0)), alpha=.5)
+plt.clf()
-plt.plot(model.t, abs(model.bcf(0)(model.t)) / abs(model.bcf(0)(model.t)).max())
-plt.plot(model.t, model.L[0].operator_norm(model.t))
-plt.plot(model.t, model.H.operator_norm(model.t) - ω_0/2)
-plt.plot(model.t, np.exp(- model.t * np.min(np.array(model.bcf_coefficients()[1][0]).real)))
-#plt.xlim(0,10)
 
-ωs = np.linspace(0.01, 4*ω_0, 10000)
+#plt.plot(model.t, np.abs(thermal_bcf(model.t))/np.abs(thermal_bcf(0)), alpha=.5, zorder=5, label="BCF (cold bath)")
-mx = (model.spectral_density(1)(ωs) + model.thermal_spectral_density(1)(ωs)).max()
+#plt.plot(model.t, abs(model.bcf(0)(model.t)) / abs(model.bcf(0)(model.t)).max(), linewidth=1)
-plt.plot(ωs, model.spectral_density(0)(ωs)/mx + model.thermal_spectral_density(0)(ωs)/mx)
+plt.plot(model.t, model.L[0].operator_norm(model.t), linewidth=1, label=r"$||L||$")
-plt.plot(ωs, model.spectral_density(1)(ωs)/mx + model.thermal_spectral_density(1)(ωs)/mx)
+plt.plot(model.t, (model.H.operator_norm(model.t) - ω_0/2)/ω_0, linewidth=.5, label=r"$||H_{\mathrm{S}}||$")
+plt.plot(model.t, np.exp(- model.t * np.min(np.array(model.bcf_coefficients()[1][0]).real)), label="Slowest Decaying BCF Term")
+plt.xlabel(r"$\tau$")
+plt.legend()
+#fs.export_fig("modulation_setup")
+
+ωs = np.linspace(0.01, 4 * ω_0, 10000)
+
+def total_sd(n, ω):
+    return model.bcf_scales[n] * (
+        model.spectral_density(n)(ω)
+        * (1/(np.expm1(ω / model.T[n])) + np.heaviside(ω, 0))
+    )
+
+
+mx =  1#(total_sd(1, ωs)).max()
+
+plt.clf()
+plt.plot(
+    ωs,
+    total_sd(0, ωs) / mx,
+    label="Cold Bath",
+)
+plt.plot(
+    ωs,
+    total_sd(1, ωs) / mx,
+    label="Hot Bath",
+)
 plt.axvline(model.H.operator_norm(0) * 2 + Δ)
 plt.axvline(model.H.operator_norm(0) * 2 - Δ)
-plt.plot(ωs, np.sinc((ωs - ω_0 - Δ) * τ_s * cycles))
+plt.plot(ωs, np.sinc((ωs - ω_0 - Δ) * τ_s * cycles), label=rf"Filter $n={cycles}$")
-plt.plot(ωs, np.sinc((ωs - ω_0 - Δ) * τ_s * cycles * 2), color="orange", linewidth=.5)
+plt.plot(
-plt.plot(ωs, np.sinc((ωs - ω_0 - Δ) * τ_s * cycles * 10), color="yellow", linewidth=.4)
+    ωs,
-#plt.xlim(2, 4)
+    np.sinc((ωs - ω_0 - Δ) * τ_s * cycles * 2),
+    color="orange",
+    linewidth=0.5,
+    label=rf"Filter $n={cycles*2}$",
+)
+plt.plot(
+    ωs,
+    np.sinc((ωs - ω_0 - Δ) * τ_s * cycles * 10),
+    color="yellow",
+    linewidth=0.2,
+    label=rf"Filter $n={cycles*10}$",
+)
+plt.xlim(8, 50)
+plt.xlabel(r"$\omega$")
+plt.legend()
+fs.export_fig("sd_setup")

python/energy_flow_proper/10_antizeno_engine/with_gap.org

@@ -0,0 +1,148 @@
+#+PROPERTY: header-args :session 10_anti_zeno_engine_gap :kernel python :pandoc no :async yes :tangle 10_first_anti_zeno_cc.py
+Here we try to reproduce the anti zeno engine from the paper.
+
+* Boilerplate
+#+begin_src jupyter-python :results none
+  import figsaver as fs
+  from hiro_models.one_qubit_model import QubitModelMutliBath, 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
+  from hops.util.dynamic_matrix import SmoothStep, Periodic, Harmonic, ConstantMatrix, Piecewise, Shift
+  import matplotlib.pyplot as plt
+  import numpy as np
+  import matplotlib
+  import plot_utils as pu
+#+end_src
+
+Init ray and silence stocproc.
+#+begin_src jupyter-python
+  import ray
+  ray.shutdown()
+  ray.init(address="141.30.17.16:6379")
+#+end_src
+
+#+RESULTS:
+: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '141.30.17.225', 'raylet_ip_address': '141.30.17.225', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-03_14-53-16_736667_825118/sockets/plasma_store.7', 'raylet_socket_name': '/tmp/ray/session_2022-08-03_14-53-16_736667_825118/sockets/raylet.3', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-03_14-53-16_736667_825118', 'metrics_export_port': 59318, 'gcs_address': '141.30.17.16:6379', 'address': '141.30.17.16:6379', 'node_id': '4261dafb07982d866fc5f9805325751682bf9a59ad6e8cfeafab3a7a'})
+
+#+begin_src jupyter-python :results none
+  from hops.util.logging_setup import logging_setup
+  import logging
+  logging_setup(logging.INFO, show_stocproc=False)
+#+end_src
+
+* Model Definition
+#+begin_src jupyter-python :results none
+  from anti_zeno_engine import *
+#+end_src
+
+#+begin_src jupyter-python
+  (
+      model,
+      params
+  ) = anti_zeno_engine(
+      Δ=11,
+      ε=.1,
+      ω_c=1,
+      ε_couple=0.7,
+      n=10,
+      detune=.5,
+      ω_0=20,
+      T_c=1e3,
+      T_h=1e4,
+      δ=[3.2*.01/4, 1*.01/4],
+      γ=.2,
+      switch_cycles=1,
+      therm_initial_state=False,
+      ε_init=.00001/2,
+      terms=7,
+      dt=0.01,
+      sp_tol=1e-3,
+      init_time_steps=10,
+  )
+  model.k_max = 3
+#+end_src
+
+#+RESULTS:
+
+#+begin_src jupyter-python
+  aux.integrate(model, 1)
+#+end_src
+
+#+RESULTS:
+#+begin_example
+  [INFO    hops.core.integration     2872983] Choosing the nonlinear integrator.
+  [INFO    hops.core.integration     2872983] Using 21 integrators.
+  [INFO    hops.core.integration     2872983] Some 1 trajectories have to be integrated.
+  [INFO    hops.core.integration     2872983] Using 680 hierarchy states.
+    0% 0/1 [00:00<?, ?it/s][INFO    hops.core.signal_delay    2872983] caught sig 'SIGINT'
+  [INFO    hops.core.signal_delay    2873543] caught sig 'SIGINT'
+  [INFO    hops.core.signal_delay    2873545] caught sig 'SIGINT'
+  [INFO    hops.core.signal_delay    2873542] caught sig 'SIGINT'
+  [INFO    hops.core.signal_delay    2873537] caught sig 'SIGINT'
+  [INFO    hops.core.signal_delay    2872983] caught sig 'SIGINT'
+#+end_example
+
+Let's test the assumptions of the paper.
+
+#+begin_src jupyter-python
+  params.cycles
+#+end_src
+
+#+RESULTS:
+: 9
+
+#+begin_src jupyter-python
+  plot_az_coupling_diagram(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:xlabel= | $\tau$ | > |
+[[file:./.ob-jupyter/8f2785b27b379ec3b747417c016c2570ce97e1de.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_az_sd_overview(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:xlabel= | $\omega$ | > |
+[[file:./.ob-jupyter/faf0641b41b72d81d9b03797ef6791f9292383b1.svg]]
+:END:
+
+
+#+begin_src jupyter-python
+  plot_total_power(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:> |
+[[file:./.ob-jupyter/591e7e4c4bbaf60e276e6f0a64c69b9ea09cbfbf.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_excited_state(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:> |
+[[file:./.ob-jupyter/3f26d886f102ce62e82565c14d154a8fc08b0668.svg]]
+:END:
+
+#+begin_src jupyter-python
+  plot_power_output(model, params)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <Figure | size | 520x320 | with | 1 | Axes> | <AxesSubplot:title= | (center : \(\bar{P}=-0.0016\pm 0.0006\)) | xlabel= | $\tau$ | ylabel= | $\langle{P}\rangle$ | > | EnsembleValue | (((10000 -0.0016086506160326898 0.0006489645812092305))) |
+[[file:./.ob-jupyter/5a037227671990870089ffb27396fa2c66fd3628.svg]]
+:END: