update project 06

.gitignore

@@ -493,3 +493,9 @@ get_neq_from_truncation_scheme
 
 data
 /python/energy_flow_proper/06_two_qubit_first_experiments/.zz_xx_test/
+
+results
+results/**.json
+results/**.npy
+/python/energy_flow_proper/05_gaussian_two_baths/data_conv/
+/python/energy_flow_proper/05_gaussian_two_baths/.git.b/

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.org

@@ -25,7 +25,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
   s = 1
   Ω = 1
   Λ = 1
-  γ = 0.1
+  γ = 1
 
 
   hops_bcf = hops.util.bcf.OhmicBCF_zeroTemp(
@@ -35,17 +35,22 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
   )
 
 
+  max_HO_level_set = 10
   def make_config(
       max_HO_level: int,
       bcf_terms: int = 4,
       t_max: float = 10,
       k_fac: float = 1.4,
       inf_tol: float = 0.1,
+          k_max: int = 7,
       sp_tol: float = 1e-3,
       bcf_scale: list[float] = [0.5, 0.5],
       truncation_scheme: str = "power",
       T: float = 0.3,
   ):
+      global max_HO_level_set
+      max_HO_level_set = max_HO_level
+
       # The BCF fit
       g, w = get_ohm_g_w(bcf_terms, s, wc)
 
@@ -94,7 +99,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
                   system, nonlinear=True, influence_tolerance=inf_tol
               )
               if truncation_scheme == "bath_memory"
-              else TruncationScheme_Simplex(7),
+              else TruncationScheme_Simplex(k_max),
               save_therm_rng_seed=True,
           ),
           Eta=[
@@ -203,7 +208,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 
 - better
 *** 8 HO Levels
-#+begin_src jupyter-python :results none :async no :tangle config.py
+#+begin_src jupyter-python :results none :async no
   params = make_config(
       max_HO_level=8,
       bcf_terms=4,
@@ -215,7 +220,22 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
       truncation_scheme="bath_memory"
   )
 #+end_src
+
 *** 9 HO Levels, Higher Temp
+#+begin_src jupyter-python :results none :async no
+  params = make_config(
+      max_HO_level=9,
+      bcf_terms=5,
+      t_max=50,
+      inf_tol=4,
+      sp_tol=1e-5,
+      bcf_scale=[0.4, 0.4],
+      T=1,
+      truncation_scheme="bath_memory",
+  )
+#+end_src
+
+*** 9 Levels, Slightly Lower Temp, Less Terms
 #+begin_src jupyter-python :results none :async no
   params = make_config(
       max_HO_level=9,
@@ -223,12 +243,69 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
       t_max=50,
       inf_tol=.3,
       sp_tol=1e-5,
-      bcf_scale=[.4, .4],
-      T=1,
-      truncation_scheme="bath_memory"
+      bcf_scale=[0.3, 0.3],
+      T=.8,
+      truncation_scheme="bath_memory",
   )
 #+end_src
 
+*** 6 Levels
+#+begin_src jupyter-python :results none :async no
+  params = make_config(
+      max_HO_level=6,
+      bcf_terms=4,
+      t_max=50,
+      inf_tol=0.1,
+      sp_tol=1e-5,
+      bcf_scale=[0.25, 0.25],
+      T=0.7,
+      truncation_scheme="bath_memory",
+  )
+#+end_src
+
+*** 8 Levels
+#+begin_src jupyter-python :results none :async no
+  params = make_config(
+      max_HO_level=8,
+      bcf_terms=4,
+      t_max=50,
+      inf_tol=0.1,
+      sp_tol=1e-5,
+      bcf_scale=[0.25, 0.25],
+      T=0.6,
+      truncation_scheme="bath_memory",
+  )
+#+end_src
+
+*** 9 Levels
+#+begin_src jupyter-python :results none :async no
+  params = make_config(
+      max_HO_level=9,
+      bcf_terms=4,
+      t_max=50,
+      inf_tol=.1,
+      sp_tol=1e-5,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      truncation_scheme="bath_memory",
+  )
+#+end_src
+
+*** 9 Levels
+#+begin_src jupyter-python :results none :async no :tangle config.py
+  γ = .5
+  params = make_config(
+     max_HO_level=9,
+     bcf_terms=5,
+     t_max=50,
+     inf_tol=.05,
+     sp_tol=1e-5,
+     bcf_scale=[0.2, 0.2],
+     T=.6,
+     k_max=5,
+     truncation_scheme="simplex",
+  )
+#+end_src
 
 ** Setup Boilerplate
 #+begin_src jupyter-python :results none
@@ -237,22 +314,24 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
   import matplotlib.pyplot as plt
   from hopsflow import util, hopsflow
   import utilities as ut
+  import figsaver as fs
+
   from hops.core.hierarchy_data import HIMetaData
   import logging
   logger = logging.getLogger()
   logger.setLevel(logging.INFO)
 #+end_src
 
-** Convergence Testing
+** Convergence Testing 9
 #+begin_src jupyter-python :results none
   with_power = make_config(
       max_HO_level=9,
       bcf_terms=4,
       t_max=50,
-      inf_tol=.3,
       sp_tol=1e-5,
-      bcf_scale=[.4, .4],
-      T=1,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      k_fac=1.5,
   )
 #+end_src
 
@@ -262,9 +341,9 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
       bcf_terms=4,
       t_max=50,
       sp_tol=1e-5,
-      bcf_scale=[.4, .4],
-      T=1,
-      inf_tol=1,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      inf_tol=.1,
       truncation_scheme="bath_memory"
   )
 #+end_src
@@ -274,10 +353,10 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
       max_HO_level=9,
       bcf_terms=4,
       t_max=50,
-      inf_tol=.3,
       sp_tol=1e-5,
-      bcf_scale=[.4, .4],
-      T=1,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      k_max=5,
       truncation_scheme="simplex",
   )
 #+end_src
@@ -288,16 +367,30 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
   supervisors = [HOPSSupervisor(config, number_of_samples=1, data_name="data_conv") for config in configs]
 #+end_src
 
+#+begin_src jupyter-python
+  for conf in configs:
+      print(len(conf.indices))
+#+end_src
+
+#+RESULTS:
+: 1287
+: 798
+: 900
+
 #+begin_src jupyter-python
   for sup in supervisors:
+      print(len(sup.params.indices))
       sup.integrate_single_process()
 #+end_src
 
 #+RESULTS:
 :RESULTS:
 : INFO:hops.core.integration:Some 1 trajectories have to be integrated.
-:   0% 0/1 [00:00<?, ?it/s]INFO:hops.core.integration:Some 1 trajectories have to be integrated.
-:   0% 0/1 [00:00<?, ?it/s]INFO:hops.core.integration:Some 1 trajectories have to be integrated.
+: 1287
+:   0% 0/1 [00:00<?, ?it/s]798
+: INFO:hops.core.integration:Some 1 trajectories have to be integrated.
+:   0% 0/1 [00:00<?, ?it/s]900
+: INFO:hops.core.integration:Some 1 trajectories have to be integrated.
 :   0% 0/1 [00:00<?, ?it/s]
 :END:
 
@@ -311,12 +404,17 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-: 6435
-: 543
-: 676
+: 1287
+: 798
+: 900
 
 #+begin_src jupyter-python :results none
-  diffs = [np.linalg.norm(trajectories[0] - trajectories[1], axis=1), np.linalg.norm(trajectories[0] - trajectories[2], axis=1)]
+  diffs = [
+      np.linalg.norm(trajectories[0] - trajectories[1], axis=1)
+      / np.linalg.norm(trajectories[0], axis=1),
+      np.linalg.norm(trajectories[0] - trajectories[2], axis=1)
+      / np.linalg.norm(trajectories[0], axis=1),
+  ]
 #+end_src
 
 #+begin_src jupyter-python
@@ -327,17 +425,120 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/b06c30144eebd8d3d0380b8f8acd3fcc169fc5b9.png]]
-** Init Gaussflow
+[[file:./.ob-jupyter/30b2331fd5aa64028e68709753ab52aa8af2ab56.svg]]
+
+** Convergence Testing 10
 #+begin_src jupyter-python :results none
+  with_power = make_config(
+      max_HO_level=10,
+      bcf_terms=4,
+      t_max=50,
+      sp_tol=1e-5,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      k_fac=1.8,
+  )
+#+end_src
+
+#+begin_src jupyter-python :results none
+  with_inftol = make_config(
+      max_HO_level=10,
+      bcf_terms=4,
+      t_max=50,
+      sp_tol=1e-5,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      inf_tol=.05,
+      truncation_scheme="bath_memory"
+  )
+#+end_src
+
+#+begin_src jupyter-python :results none
+  with_simplex = make_config(
+      max_HO_level=10,
+      bcf_terms=5,
+      t_max=50,
+      sp_tol=1e-5,
+      bcf_scale=[0.2, 0.2],
+      T=.6,
+      k_max=6,
+      truncation_scheme="simplex",
+  )
+#+end_src
+
+#+begin_src jupyter-python :results none
+  configs = [with_simplex, with_inftol, with_power]
+  from hops.core.integration import HOPSSupervisor
+  supervisors = [HOPSSupervisor(config, number_of_samples=1, data_name="data_conv") for config in configs]
+#+end_src
+
+#+begin_src jupyter-python
+  for conf in configs:
+      print(len(conf.indices))
+#+end_src
+
+#+RESULTS:
+: 8008
+: 2142
+: 2209
+
+#+begin_src jupyter-python
+  for sup in supervisors:
+      print(len(sup.params.indices))
+      sup.integrate_single_process()
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: 8008
+: INFO:hops.core.integration:Some 1 trajectories have to be integrated.
+:   0% 0/1 [00:00<?, ?it/s]
+:END:
+
+#+begin_src jupyter-python
+  trajectories = []
+  for supervisor in supervisors:
+        with supervisor.get_data(True) as data:
+            τ = data.get_time()
+            print(len(data.params.indices))
+            trajectories.append(data.get_stoc_traj(0))
+#+end_src
+
+#+RESULTS:
+: 890aec2e-817c-4312-9a77-26cd35e5db0e
+
+#+begin_src jupyter-python :results none
+  diffs = [
+      np.linalg.norm(trajectories[0] - trajectories[1], axis=1)
+      / np.linalg.norm(trajectories[0], axis=1),
+      np.linalg.norm(trajectories[0] - trajectories[2], axis=1)
+      / np.linalg.norm(trajectories[0], axis=1),
+  ]
+#+end_src
+
+#+begin_src jupyter-python
+  with ut.hiro_style():
+      plt.plot(τ, diffs[0])
+      plt.plot(τ, diffs[1], label="power")
+      plt.legend()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/e2bb2f09a3c7beb51b771c7ac321bc80d9ef3d46.svg]]
+
+** Init Gaussflow
+#+begin_src jupyter-python
   from hopsflow import gaussflow_two as gf
   from hopsflow import util as util
 
-  α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=12, resolution=0.1/2)
+  α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=12, resolution=0.01/5)
 #+end_src
 
+#+RESULTS:
+: INFO:root:Looking up bcf fit at .cache/bcf_fit/2328985928938661957.npy.
+: INFO:root:Loading bcf fit from .cache/bcf_fit/2328985928938661957.npy.
 
-#+begin_src jupyter-python :results none
+#+begin_src jupyter-python
   α = util.BCF(
       params.IntP.t_max,
       hops.util.bcf.OhmicBCF_nonZeroTemp(
@@ -346,33 +547,43 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
           wc,
           beta=1.0 / (params.SysP.__non_key__["T"][0]),
       ),
-      num_terms=9,
-      resolution=0.1/2,
+      num_terms=12,
+      resolution=0.01/5,
   )
 #+end_src
 
+#+RESULTS:
+: /home/hiro/src/hops/hops/util/bcf.py:213: UserWarning: this implementation uses mpmath to evaluate the zeta_function! for a better performance consider the 'OhmEnv' package
+:   warnings.warn(
+: INFO:root:Looking up bcf fit at .cache/bcf_fit/-3203860338031140913.npy.
+: INFO:root:Loading bcf fit from .cache/bcf_fit/-3203860338031140913.npy.
+
 #+begin_src jupyter-python
   t_test = np.linspace(0, α.t_max, 4000)
-  with ut.hiro_style():
-        fig, ax = ut.plot_complex(t_test, α(t_test) - α.approx(t_test), label="finite")
-        fig, ax = ut.plot_complex(t_test, α_0(t_test) - α_0.approx(t_test), label="zero")
+
+  with fs.hiro_style():
+        fig, ax = fs.plot_complex(t_test, α(t_test) - α.approx(t_test), label="$α$")
+        fs.plot_complex(t_test, α_0(t_test) - α_0.approx(t_test), label="$α_0$", ax=ax)
+        fig.set_size_inches(fs.get_figsize(239, 1, .8))
+
+        ax.set_xlabel("$τ$")
+        fs.export_fig("bcf_fits")
+
+  # fs.tex_value(len(α.factors), prefix="K=", save="bcf_terms")
 #+end_src
 
 #+RESULTS:
-:RESULTS:
-[[file:./.ob-jupyter/7860e5186a807c5b41a3d8b01e841284e959ff4b.svg]]
-[[file:./.ob-jupyter/23c7abbd150839cddfec7ce4f096451f25b6fa04.svg]]
-:END:
+[[file:./.ob-jupyter/a8f04e9a2da787db38bb667c3ff290a1e782348e.svg]]
 
 #+begin_src jupyter-python :results none
-  gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2)
+  gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2, root_tol=1e-5)
 #+end_src
 
 * Load Hops Data
  #+begin_src jupyter-python
    class result:
        hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
-       N = hd.samples
+       N = hd.samples - 200
        seeds = hd.rng_seed
        τ = hd.get_time()
        ψ_1 = hd.aux_states
@@ -380,10 +591,22 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
        ρ = hd.get_rho_t()
 
    result.N
+   fs.tex_value(result.N, prefix="N=", save="samples")
 #+end_src
 
 #+RESULTS:
-: 10000
+: \(N=5191\)
+
+#+begin_src jupyter-python
+  fs.tex_value(
+      params.SysP.bcf_scale[0], prec=1, save="bcf_scale", prefix="η="
+  ), fs.tex_value(wc, prec=0, save="cutoff_freq", prefix="ω_c="), fs.tex_value(
+      params.SysP.__non_key__["T"][0], prec=1, save="temp", prefix="T="
+  ), fs.tex_value(Ω, prec=0, save="omega", prefix="Ω="), fs.tex_value(max_HO_level_set, save="max_HO_level")
+#+end_src
+
+#+RESULTS:
+| \(η=0.2\) | \(ω_c=2\) | \(T=0.6\) | \(Ω=1\) | \(9\) |
 
 * System Energy
 #+begin_src jupyter-python
@@ -392,55 +615,60 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
       params.SysP.H_sys,
       result.N,
       params.HiP.nonlinear,
-      every=1000,
+      every=int(result.N // 5),
       real=True,
-      save="system_energy"
+      save="system_energys"
   )
 #+end_src
 
 #+RESULTS:
-: WARNING:root:Loading cache from: results/system_energy_sandwhich_operator_10000_1000_a34265c8f80fbc84d6a3b074cffc30d9abe302353134db3cf684ee0c0fded8ea.npy
-
+: WARNING:root:Loading cache from: results/system_energys_sandwhich_operator_5191_1038_7b17e8304ee2a0f3b39c1715ca57d4999edd9d0fe23c629cbda401088cd94d0b.npy
 
 #+begin_src jupyter-python
-  with ut.hiro_style():
-      fig, ax = ut.plot_convergence(result.τ, e_sys, transform=np.real, linestyle="dashdot")
+  with fs.hiro_style():
+      fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real, linestyle="dashdot")
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/222be3d5cd5d5b527ea800e08bf73b2ee0de89ef.svg]]
+[[file:./.ob-jupyter/d35e70d489882edecd6565bf17fed76bdb242b60.svg]]
 
 Now with Gaussflow.
 #+begin_src jupyter-python
-  C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=[α, None])
+  C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=[α, α_0])
   energy_gf = C.system_energy(result.τ)
 #+end_src
 
 #+RESULTS:
 
-
 #+begin_src jupyter-python
-  with ut.hiro_style():
-      ax.plot(result.τ, energy_gf, color="green", linestyle="dashdot")
- #plt.xlim(30,40)
-      #plt.ylim(-.005,.005)
-  fig
+  with fs.hiro_style():
+      fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real)
+      fig.set_size_inches(fs.get_figsize(239, 1, .8))
+
+      ax.plot(result.τ, energy_gf, color="green", label="analytic", linestyle="--", zorder=1000)
+      ax.set_ylabel(r"$\langle H_S\rangle$")
+      ax.set_xlabel(r"$τ$")
+      ax.legend()
+      fs.export_fig("system_energy")
+      # plt.xlim(0,30)
+      # plt.ylim(1.5, 1.6)
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/78cf68864d651dd1405f0dfa80f7c831a7bde996.svg]]
+[[file:./.ob-jupyter/c574ce2d97056f5772a61da13eca6d4b8ac5ab7a.svg]]
 
 
 #+begin_src jupyter-python
-  with ut.hiro_style():
-      ut.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
-      # plt.xlim(35,40)
-      # plt.ylim(0,.01)D
-      plt.legend()
+  with fs.hiro_style():
+      fig, ax = fs.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
+      fig.set_size_inches(fs.get_figsize(239, 1, .8))
+      ax.set_xlabel(r"$τ$")
+      ax.legend()
+      fs.export_fig("e_sys_convergence", fig)
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/4e49054f9b64e79c49a9c4af01b2daf9ffa0124d.svg]]
+[[file:./.ob-jupyter/ba0a30288d48d6da5e22938d414e4e6bffdec659.svg]]
 
 * Thermal State
 #+begin_src jupyter-python
@@ -451,25 +679,29 @@ Now with Gaussflow.
 
 #+RESULTS:
 #+begin_example
-  array([ 1.04772256,  2.04772255,  2.14316767,  3.04772255,  3.14316761,
-          3.23861273,  4.04760557,  4.14204639,  4.18286322,  4.18458163,
-          4.2394244 ,  4.33415536,  5.04772633,  5.14275084,  5.21406448,
-          5.22036624,  5.24633389,  5.33647238,  5.430074  ,  6.0492956 ,
-          6.15007353,  6.250944  ,  6.25844588,  6.25964261,  6.34972826,
-          6.43985241,  6.52787895,  7.05033123,  7.14899946,  7.24380084,
-          7.29230796,  7.29395284,  7.31504986,  7.35756411,  7.46316661,
-          7.55516337,  7.63134403,  8.05550357,  8.17553157,  8.30665127,
-          8.33818705,  8.35038921,  8.50017485,  8.64458643,  8.73500426,
-          9.07062116,  9.24040077,  9.36698584,  9.41225713,  9.48373566,
-          9.70479266,  9.83379515, 10.10351883, 10.14569718, 10.37338888,
-         10.71719495, 10.71787655, 10.95409824, 11.17729002, 11.64575746,
-         11.99096368, 12.35910821, 13.11566094, 13.80670653])
+  array([ 1.20710676,  2.20710635,  2.62131951,  3.20707431,  3.62124166,
+          4.03548025,  4.2070741 ,  4.62090764,  5.03449736,  5.1680822 ,
+          5.26132568,  5.32641331,  5.44895674,  5.66761406,  6.06336422,
+          6.2064203 ,  6.36059303,  6.3621128 ,  6.45468408,  6.63682254,
+          6.85898472,  7.09723329,  7.22842636,  7.44603437,  7.47571157,
+          7.53606078,  7.79286457,  7.89128096,  8.24644979,  8.25842562,
+          8.27423355,  8.59867018,  8.60800633,  8.68840482,  8.76148041,
+          9.07963262,  9.17511679,  9.23759914,  9.32094813,  9.39108179,
+          9.56706803,  9.73835633,  9.7611353 ,  9.87036682,  9.88238241,
+         10.19179189, 10.31273478, 10.58379054, 10.58415872, 10.66166785,
+         10.9022923 , 11.1194365 , 11.16339017, 11.31589812, 11.38936437,
+         11.52631398, 11.79148772, 12.17989834, 12.18134628, 12.20571458,
+         12.26147546, 12.53524288, 12.86610389, 12.89302842, 13.06538653,
+         13.20646914, 13.8943837 , 14.06670366, 14.11916572, 14.85727309,
+         15.15284881, 15.6347988 , 15.67239516, 15.97586999, 16.09266967,
+         18.44274693, 18.59751137, 18.61468424, 18.87454828, 22.72078481,
+         22.84247174])
 #+end_example
 
 
 #+begin_src jupyter-python :results none
   def thermal_e(v: np.ndarray, β: float):
-        return np.sum(v * np.exp(-β * v)) / np.sum(np.exp(-β * v))
+        return np.sum(v[None, :] * np.exp(-β * v[None, :]), axis=1) / np.sum(np.exp(-β * v[None,:]), axis=1)
 
   def trace_norm(ρ: np.ndarray):
         return np.trace(sc.linalg.sqrtm(ρ @ ρ.conj().T)).real
@@ -480,7 +712,7 @@ Now with Gaussflow.
 #+end_src
 
 #+RESULTS:
-: 1.0591063695586198
+: array([1.59542423])
 
 #+begin_src jupyter-python
   ρ_therm = sc.linalg.expm(-params.SysP.H_sys/(params.SysP.__non_key__["T"][0]))
@@ -489,7 +721,25 @@ Now with Gaussflow.
 #+end_src
 
 #+RESULTS:
-: 0.07726133069400465
+: 0.20088856709757003
+
+#+begin_src jupyter-python
+  e_final = e_sys[-1][1][-1]
+  e_final
+#+end_src
+
+#+RESULTS:
+: 1.4292409539116322
+
+#+begin_src jupyter-python
+  import scipy.optimize
+  closest_temp = scipy.optimize.minimize(lambda β: np.abs(thermal_e(v, β) - e_final)**2, 1.0 / (params.SysP.__non_key__["T"][0]))
+  1/closest_temp.x, (params.SysP.__non_key__["T"][0])
+#+end_src
+
+#+RESULTS:
+| array | ((0.48000139)) | 0.6 |
+Interestingly in the middle.
 
 * Heat Flow
 #+begin_src jupyter-python :results none
@@ -501,7 +751,6 @@ Now with Gaussflow.
       L=params.SysP.L,
       G=params.SysP.g,
       W=params.SysP.w,
-      bcf_scale=[1, 1],
       fock_hops=True
   )
   run = hf.HOPSRun(np.array(result.ψ[0, :]), np.array(result.ψ_1[0, :]), sys)
@@ -526,8 +775,8 @@ Now with Gaussflow.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f6f9892ae80> |
-[[file:./.ob-jupyter/1d3f0e09cc8a1373f485a5feded2f90d6a589871.svg]]
+| <matplotlib.lines.Line2D | at | 0x7f060bd82700> |
+[[file:./.ob-jupyter/494f4a872592e6928904bfe575a46b0f7a14a40a.svg]]
 :END:
 
 #+begin_src jupyter-python
@@ -535,50 +784,93 @@ Now with Gaussflow.
       iter(result.ψ),
       iter(result.ψ_1),
       sys,
-      500,
+      result.N,
       (iter(result.seeds), therm_params),
-      every=int(result.N / 5),
-      save="heat_flow"
+      every=int(result.N / 20),
+      save="heat_flow1"
   )
 #+end_src
 
 #+RESULTS:
-: 100% 499/499 [01:15<00:00,  6.60it/s]
-: INFO:root:Writing cache to: results/heat_flow__heat_flow_ensemble_body_500_2000_682ce0d1b2d8b028cb7c6a02e6faff13ee5fb01582be1253db6bbe9db4e15312.npy
+: WARNING:root:Loading cache from: results/heat_flow1__heat_flow_ensemble_body_5191_259_8cdbcb067351c07e5d518185d3ebcf81eae0e54c15f4682f5c460f7b867a50d3.npy
 
 #+begin_src jupyter-python
   with ut.hiro_style():
       _, ax = ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=0)
-      ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1, ax=ax)
+      _, ax = ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1)
+      #ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1, ax=ax)
       ax.legend()
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/bac8b6fe9327df5899f5d317aa4af42cfb092b43.svg]]
+:RESULTS:
+[[file:./.ob-jupyter/66e584bd18b5548fe7c1baf5fc243539e23e6841.svg]]
+[[file:./.ob-jupyter/eaa35bd6cdea4956b1360c0e9bda4212a45cc8a1.svg]]
+:END:
 
 
-#+begin_src jupyter-python :results none
-  exact_flow = [C.flow(result.τ, 0), C.flow(result.τ, 1)]
-#+end_src
-
 #+begin_src jupyter-python
-  plt.plot(result.τ, exact_flow[0])
-  plt.plot(result.τ, -full_flow[-1][1][0])
+  exact_flow = [C.flow(result.τ, 0, False), C.flow(result.τ, 1, False)]
+#+end_src
+
+#+RESULTS:
+
+#+begin_src jupyter-python
+  for bath in range(2):
+      with fs.hiro_style():
+          fig, ax = ut.plot_convergence(
+              result.τ, full_flow[::5], transform=lambda y: -y, bath=bath
+          )
+          fig.set_size_inches(fs.get_figsize(239, 1, .8))
+          ax.set_ylabel("$-J$")
+          ax.set_xlabel("$τ$")
+          ax.plot(
+              result.τ,
+              exact_flow[bath],
+              label="analytic",
+              color="green",
+              linestyle="--",
+              zorder=1000,
+          )
+          ax.set_title(f"Bath {bath + 1}, $T={params.SysP.__non_key__['T'][bath]}$")
+          ax.legend(loc="upper right")
+          fs.export_fig(f"flow_{bath}")
 #+end_src
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f6f70916670> |
-[[file:./.ob-jupyter/755f01900a0ee80481360bbf33eea3c62c2c717d.svg]]
+[[file:./.ob-jupyter/0bdb5ce9fa4f4929667ce0d8398dbc3b6c8441a2.svg]]
+[[file:./.ob-jupyter/834cd5ce6f522f761ae1e5214a3d71c6a18c0918.svg]]
 :END:
 
+
 #+begin_src jupyter-python
-  with ut.hiro_style():
-        ut.plot_diff_vs_sigma(result.τ, [(n, flow[1], tol[1]) for (n, flow, tol) in full_flow], -exact_flow[1])
-        # plt.xlim(35,40)
-        # plt.ylim(0,.01)D
-        plt.legend()
+  for bath in range(2):
+      with fs.hiro_style():
+          fig, ax = fs.plot_diff_vs_sigma(
+              result.τ,
+              full_flow[::5],
+              -exact_flow[bath],
+              bath=bath,
+              ylabel="$J$"
+          )
+          fig.set_size_inches(fs.get_figsize(239, 1, .8))
+
+          ax.legend()
+          ax.set_xlabel("$τ$")
+          ax.legend()
+          ax.set_title(f"Bath {bath + 1}")
+          fs.export_fig(f"flow_conv_{bath}")
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/9620d5702095804e22f9572a2266d3e4a3699180.svg]]
+:RESULTS:
+[[file:./.ob-jupyter/c2efc39704e07f80cb3bb4b1f8655e2659e81fa1.svg]]
+[[file:./.ob-jupyter/dca19c821f5a6f25aae2d7168bbcb19f4820faa0.svg]]
+:END:
+
+* Thoughts
+- had to take many hierarchy states
+- initial slip appears sensitive to HO hilbert space truncation
+- maybe more investigation needed
+- initial slip -> both baths similar, then steady flow

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.py

@@ -36,15 +36,70 @@ params = make_config(
     T=.2
 )
 
+params = make_config(
+    max_HO_level=8,
+    bcf_terms=4,
+    t_max=50,
+    inf_tol=.2,
+    sp_tol=1e-6,
+    bcf_scale=[.3, .3],
+    T=.2,
+    truncation_scheme="bath_memory"
+)
+
+params = make_config(
+    max_HO_level=9,
+    bcf_terms=5,
+    t_max=50,
+    inf_tol=4,
+    sp_tol=1e-5,
+    bcf_scale=[0.4, 0.4],
+    T=1,
+    truncation_scheme="bath_memory",
+)
+
 params = make_config(
     max_HO_level=9,
     bcf_terms=4,
     t_max=50,
     inf_tol=.3,
     sp_tol=1e-5,
-    bcf_scale=[.4, .4],
-    T=1,
-    truncation_scheme="bath_memory"
+    bcf_scale=[0.3, 0.3],
+    T=.8,
+    truncation_scheme="bath_memory",
+)
+
+params = make_config(
+    max_HO_level=6,
+    bcf_terms=4,
+    t_max=50,
+    inf_tol=0.1,
+    sp_tol=1e-5,
+    bcf_scale=[0.25, 0.25],
+    T=0.7,
+    truncation_scheme="bath_memory",
+)
+
+params = make_config(
+    max_HO_level=8,
+    bcf_terms=4,
+    t_max=50,
+    inf_tol=0.1,
+    sp_tol=1e-5,
+    bcf_scale=[0.25, 0.25],
+    T=0.6,
+    truncation_scheme="bath_memory",
+)
+
+params = make_config(
+    max_HO_level=9,
+    bcf_terms=4,
+    t_max=50,
+    inf_tol=.1,
+    sp_tol=1e-5,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    truncation_scheme="bath_memory",
 )
 
 import numpy as np
@@ -52,6 +107,8 @@ import scipy as sc
 import matplotlib.pyplot as plt
 from hopsflow import util, hopsflow
 import utilities as ut
+import figsaver as fs
+
 from hops.core.hierarchy_data import HIMetaData
 import logging
 logger = logging.getLogger()
@@ -61,10 +118,10 @@ with_power = make_config(
     max_HO_level=9,
     bcf_terms=4,
     t_max=50,
-    inf_tol=.3,
     sp_tol=1e-5,
-    bcf_scale=[.4, .4],
-    T=1,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    k_fac=1.5,
 )
 
 with_inftol = make_config(
@@ -72,9 +129,9 @@ with_inftol = make_config(
     bcf_terms=4,
     t_max=50,
     sp_tol=1e-5,
-    bcf_scale=[.4, .4],
-    T=1,
-    inf_tol=1,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    inf_tol=.1,
     truncation_scheme="bath_memory"
 )
 
@@ -82,10 +139,10 @@ with_simplex = make_config(
     max_HO_level=9,
     bcf_terms=4,
     t_max=50,
-    inf_tol=.3,
     sp_tol=1e-5,
-    bcf_scale=[.4, .4],
-    T=1,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    k_max=5,
     truncation_scheme="simplex",
 )
 
@@ -93,7 +150,11 @@ configs = [with_simplex, with_inftol, with_power]
 from hops.core.integration import HOPSSupervisor
 supervisors = [HOPSSupervisor(config, number_of_samples=1, data_name="data_conv") for config in configs]
 
+for conf in configs:
+    print(len(conf.indices))
+
 for sup in supervisors:
+    print(len(sup.params.indices))
     sup.integrate_single_process()
 
 trajectories = []
@@ -103,7 +164,74 @@ for supervisor in supervisors:
           print(len(data.params.indices))
           trajectories.append(data.get_stoc_traj(0))
 
-diffs = [np.linalg.norm(trajectories[0] - trajectories[1], axis=1), np.linalg.norm(trajectories[0] - trajectories[2], axis=1)]
+diffs = [
+    np.linalg.norm(trajectories[0] - trajectories[1], axis=1)
+    / np.linalg.norm(trajectories[0], axis=1),
+    np.linalg.norm(trajectories[0] - trajectories[2], axis=1)
+    / np.linalg.norm(trajectories[0], axis=1),
+]
+
+with ut.hiro_style():
+    plt.plot(τ, diffs[0])
+    plt.plot(τ, diffs[1], label="power")
+    plt.legend()
+
+with_power = make_config(
+    max_HO_level=10,
+    bcf_terms=4,
+    t_max=50,
+    sp_tol=1e-5,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    k_fac=1.8,
+)
+
+with_inftol = make_config(
+    max_HO_level=10,
+    bcf_terms=4,
+    t_max=50,
+    sp_tol=1e-5,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    inf_tol=.05,
+    truncation_scheme="bath_memory"
+)
+
+with_simplex = make_config(
+    max_HO_level=10,
+    bcf_terms=5,
+    t_max=50,
+    sp_tol=1e-5,
+    bcf_scale=[0.2, 0.2],
+    T=.6,
+    k_max=6,
+    truncation_scheme="simplex",
+)
+
+configs = [with_simplex, with_inftol, with_power]
+from hops.core.integration import HOPSSupervisor
+supervisors = [HOPSSupervisor(config, number_of_samples=1, data_name="data_conv") for config in configs]
+
+for conf in configs:
+    print(len(conf.indices))
+
+for sup in supervisors:
+    print(len(sup.params.indices))
+    sup.integrate_single_process()
+
+trajectories = []
+for supervisor in supervisors:
+      with supervisor.get_data(True) as data:
+          τ = data.get_time()
+          print(len(data.params.indices))
+          trajectories.append(data.get_stoc_traj(0))
+
+diffs = [
+    np.linalg.norm(trajectories[0] - trajectories[1], axis=1)
+    / np.linalg.norm(trajectories[0], axis=1),
+    np.linalg.norm(trajectories[0] - trajectories[2], axis=1)
+    / np.linalg.norm(trajectories[0], axis=1),
+]
 
 with ut.hiro_style():
     plt.plot(τ, diffs[0])
@@ -113,7 +241,7 @@ with ut.hiro_style():
 from hopsflow import gaussflow_two as gf
 from hopsflow import util as util
 
-α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=12, resolution=0.1/2)
+α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=12, resolution=0.01/5)
 
 α = util.BCF(
     params.IntP.t_max,
@@ -123,20 +251,27 @@ from hopsflow import util as util
         wc,
         beta=1.0 / (params.SysP.__non_key__["T"][0]),
     ),
-    num_terms=9,
-    resolution=0.1/2,
+    num_terms=12,
+    resolution=0.01/5,
 )
 
 t_test = np.linspace(0, α.t_max, 4000)
-with ut.hiro_style():
-      fig, ax = ut.plot_complex(t_test, α(t_test) - α.approx(t_test), label="finite")
-      fig, ax = ut.plot_complex(t_test, α_0(t_test) - α_0.approx(t_test), label="zero")
 
-gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2)
+with fs.hiro_style():
+      fig, ax = fs.plot_complex(t_test, α(t_test) - α.approx(t_test), label="$α$")
+      fs.plot_complex(t_test, α_0(t_test) - α_0.approx(t_test), label="$α_0$", ax=ax)
+      fig.set_size_inches(fs.get_figsize(239, 1, .8))
+
+      ax.set_xlabel("$τ$")
+      fs.export_fig("bcf_fits")
+
+# fs.tex_value(len(α.factors), prefix="K=", save="bcf_terms")
+
+gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2, root_tol=1e-5)
 
 class result:
     hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
-    N = hd.samples
+    N = hd.samples - 200
     seeds = hd.rng_seed
     τ = hd.get_time()
     ψ_1 = hd.aux_states
@@ -144,41 +279,55 @@ class result:
     ρ = hd.get_rho_t()
 
 result.N
+fs.tex_value(result.N, prefix="N=", save="samples")
+
+fs.tex_value(
+    params.SysP.bcf_scale[0], prec=1, save="bcf_scale", prefix="η="
+), fs.tex_value(wc, prec=0, save="cutoff_freq", prefix="ω_c="), fs.tex_value(
+    params.SysP.__non_key__["T"][0], prec=1, save="temp", prefix="T="
+), fs.tex_value(Ω, prec=0, save="omega", prefix="Ω="), fs.tex_value(max_HO_level_set, save="max_HO_level")
 
 e_sys = util.operator_expectation_ensemble(
     iter(result.ψ),
     params.SysP.H_sys,
     result.N,
     params.HiP.nonlinear,
-    every=1000,
+    every=int(result.N // 5),
     real=True,
-    save="system_energy"
+    save="system_energys"
 )
 
-with ut.hiro_style():
-    fig, ax = ut.plot_convergence(result.τ, e_sys, transform=np.real, linestyle="dashdot")
+with fs.hiro_style():
+    fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real, linestyle="dashdot")
 
-C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=[α, None])
+C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=[α, α_0])
 energy_gf = C.system_energy(result.τ)
 
-with ut.hiro_style():
-     ax.plot(result.τ, energy_gf, color="green", linestyle="dashdot")
-#plt.xlim(30,40)
-     #plt.ylim(-.005,.005)
- fig
+with fs.hiro_style():
+    fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real)
+    fig.set_size_inches(fs.get_figsize(239, 1, .8))
 
-with ut.hiro_style():
-    ut.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
-    # plt.xlim(35,40)
-    # plt.ylim(0,.01)D
-    plt.legend()
+    ax.plot(result.τ, energy_gf, color="green", label="analytic", linestyle="--", zorder=1000)
+    ax.set_ylabel(r"$\langle H_S\rangle$")
+    ax.set_xlabel(r"$τ$")
+    ax.legend()
+    fs.export_fig("system_energy")
+    # plt.xlim(0,30)
+    # plt.ylim(1.5, 1.6)
+
+with fs.hiro_style():
+    fig, ax = fs.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
+    fig.set_size_inches(fs.get_figsize(239, 1, .8))
+    ax.set_xlabel(r"$τ$")
+    ax.legend()
+    fs.export_fig("e_sys_convergence", fig)
 
 v, w = np.linalg.eig(params.SysP.H_sys)
 v = np.sort(v.real)
 v
 
 def thermal_e(v: np.ndarray, β: float):
-      return np.sum(v * np.exp(-β * v)) / np.sum(np.exp(-β * v))
+      return np.sum(v[None, :] * np.exp(-β * v[None, :]), axis=1) / np.sum(np.exp(-β * v[None,:]), axis=1)
 
 def trace_norm(ρ: np.ndarray):
       return np.trace(sc.linalg.sqrtm(ρ @ ρ.conj().T)).real
@@ -189,13 +338,19 @@ thermal_e(v, 1.0 / (params.SysP.__non_key__["T"][0]))
 ρ_therm = ρ_therm / np.trace(ρ_therm)
 trace_norm(result.ρ[-1] - ρ_therm)
 
+e_final = e_sys[-1][1][-1]
+e_final
+
+import scipy.optimize
+closest_temp = scipy.optimize.minimize(lambda β: np.abs(thermal_e(v, β) - e_final)**2, 1.0 / (params.SysP.__non_key__["T"][0]))
+1/closest_temp.x, (params.SysP.__non_key__["T"][0])
+
 from hopsflow import hopsflow as hf
 
 sys = hf.SystemParams(
     L=params.SysP.L,
     G=params.SysP.g,
     W=params.SysP.w,
-    bcf_scale=[1, 1],
     fock_hops=True
 )
 run = hf.HOPSRun(np.array(result.ψ[0, :]), np.array(result.ψ_1[0, :]), sys)
@@ -214,24 +369,53 @@ full_flow = hopsflow.heat_flow_ensemble(
     iter(result.ψ),
     iter(result.ψ_1),
     sys,
-    500,
+    result.N,
     (iter(result.seeds), therm_params),
-    every=int(result.N / 5),
-    save="heat_flow"
+    every=int(result.N / 20),
+    save="heat_flow1"
 )
 
 with ut.hiro_style():
     _, ax = ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=0)
-    ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1, ax=ax)
+    _, ax = ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1)
+    #ut.plot_convergence(result.τ, full_flow, transform=lambda y: -y, bath=1, ax=ax)
     ax.legend()
 
-exact_flow = [C.flow(result.τ, 0), C.flow(result.τ, 1)]
+exact_flow = [C.flow(result.τ, 0, False), C.flow(result.τ, 1, False)]
 
-plt.plot(result.τ, exact_flow[0])
-plt.plot(result.τ, -full_flow[-1][1][0])
+for bath in range(2):
+    with fs.hiro_style():
+        fig, ax = ut.plot_convergence(
+            result.τ, full_flow[::5], transform=lambda y: -y, bath=bath
+        )
+        fig.set_size_inches(fs.get_figsize(239, 1, .8))
+        ax.set_ylabel("$-J$")
+        ax.set_xlabel("$τ$")
+        ax.plot(
+            result.τ,
+            exact_flow[bath],
+            label="analytic",
+            color="green",
+            linestyle="--",
+            zorder=1000,
+        )
+        ax.set_title(f"Bath {bath + 1}, $T={params.SysP.__non_key__['T'][bath]}$")
+        ax.legend(loc="upper right")
+        fs.export_fig(f"flow_{bath}")
 
-with ut.hiro_style():
-      ut.plot_diff_vs_sigma(result.τ, [(n, flow[1], tol[1]) for (n, flow, tol) in full_flow], -exact_flow[1])
-      # plt.xlim(35,40)
-      # plt.ylim(0,.01)D
-      plt.legend()
+for bath in range(2):
+    with fs.hiro_style():
+        fig, ax = fs.plot_diff_vs_sigma(
+            result.τ,
+            full_flow[::5],
+            -exact_flow[bath],
+            bath=bath,
+            ylabel="$J$"
+        )
+        fig.set_size_inches(fs.get_figsize(239, 1, .8))
+
+        ax.legend()
+        ax.set_xlabel("$τ$")
+        ax.legend()
+        ax.set_title(f"Bath {bath + 1}")
+        fs.export_fig(f"flow_conv_{bath}")

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_old_gaussflow.org

@@ -68,7 +68,7 @@ coordinate operators in a coherent state.
 #+end_src
 
  #+RESULTS:
- [[file:./.ob-jupyter/bb9441635e6497907c9917b75c02e8f904a0f1f9.svg]]
+ [[file:./.ob-jupyter/566145787562dead5fc0ba613e3241da213aa9b7.svg]]
 
 * Correlations
 And now the correlation matrix.
@@ -88,7 +88,7 @@ And now the correlation matrix.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/0fd2993ff0130aea65d99d3b4db885cc9733d554.svg]]
+[[file:./.ob-jupyter/f15a52ab49778444ed143cb46bcadb18b61fe955.svg]]
 
 #+begin_src jupyter-python
   flow = [C.flow(t_points, i) for i in range(2)]
@@ -103,7 +103,7 @@ And now the correlation matrix.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/de4d62a0fbb03cc3ddd053d161e3c87dbbd25576.svg]]
+[[file:./.ob-jupyter/965f06261af587cbf7ff303fdf8daf1aae7ac173.svg]]
 
 #+begin_src jupyter-python :results none
   from hopsflow import gaussflow as gf_old
@@ -126,7 +126,7 @@ And now the correlation matrix.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/007b510b5fdb124f2a2be6bd827739a91dc346e4.svg]]
+[[file:./.ob-jupyter/87c368143a6db2eaacde956a6d787cf32ffc3d3c.svg]]
 
 
 Checks out :P.

python/energy_flow_proper/05_gaussian_two_baths/config.py

@@ -20,7 +20,7 @@ wc = 2
 s = 1
 Ω = 1
 Λ = 1
-γ = 0.1
+γ = 1
 
 
 hops_bcf = hops.util.bcf.OhmicBCF_zeroTemp(
@@ -30,17 +30,22 @@ hops_bcf = hops.util.bcf.OhmicBCF_zeroTemp(
 )
 
 
+max_HO_level_set = 10
 def make_config(
     max_HO_level: int,
     bcf_terms: int = 4,
     t_max: float = 10,
     k_fac: float = 1.4,
     inf_tol: float = 0.1,
+        k_max: int = 7,
     sp_tol: float = 1e-3,
     bcf_scale: list[float] = [0.5, 0.5],
     truncation_scheme: str = "power",
     T: float = 0.3,
 ):
+    global max_HO_level_set
+    max_HO_level_set = max_HO_level
+
     # The BCF fit
     g, w = get_ohm_g_w(bcf_terms, s, wc)
 
@@ -89,7 +94,7 @@ def make_config(
                 system, nonlinear=True, influence_tolerance=inf_tol
             )
             if truncation_scheme == "bath_memory"
-            else TruncationScheme_Simplex(7),
+            else TruncationScheme_Simplex(k_max),
             save_therm_rng_seed=True,
         ),
         Eta=[
@@ -137,13 +142,15 @@ def make_config(
 
     return params
 
+γ = .5
 params = make_config(
-    max_HO_level=8,
-    bcf_terms=4,
-    t_max=50,
-    inf_tol=.2,
-    sp_tol=1e-6,
-    bcf_scale=[.3, .3],
-    T=.2,
-    truncation_scheme="bath_memory"
+   max_HO_level=9,
+   bcf_terms=5,
+   t_max=50,
+   inf_tol=.05,
+   sp_tol=1e-5,
+   bcf_scale=[0.2, 0.2],
+   T=.6,
+   k_max=5,
+   truncation_scheme="simplex",
 )

python/energy_flow_proper/05_gaussian_two_baths/flake.nix

@@ -14,7 +14,7 @@
       shellOverride = (oldAttrs: {
         shellHook = ''
                     export PYTHONPATH=/home/hiro/src/hopsflow/:$PYTHONPATH
-#                    export PYTHONPATH=/home/hiro/src/hops/:$PYTHONPATH
+                    export PYTHONPATH=/home/hiro/src/hops/:$PYTHONPATH
 #                    export PYTHONPATH=/home/hiro/src/stocproc/:$PYTHONPATH
                     '';
       });

python/energy_flow_proper/05_gaussian_two_baths/poetry.lock

@@ -384,11 +384,11 @@ typer = "^0.4.0"
 type = "git"
 url = "git@gitlab.hrz.tu-chemnitz.de:s8896854--tu-dresden.de/hops.git"
 reference = "BathMemoryTruncation"
-resolved_reference = "22ee9b9365203d32f0a755123b678f54237ecfd4"
+resolved_reference = "06853a819161a70f2903494633affe131369b721"
 
 [[package]]
 name = "hopsflow"
-version = "1.0.5"
+version = "1.0.6"
 description = "Calculating open system bath energy changes with HOPS and analytically."
 category = "main"
 optional = false
@@ -407,7 +407,7 @@ tqdm = "^4.62.3"
 type = "git"
 url = "https://github.com/vale981/hopsflow"
 reference = "main"
-resolved_reference = "a70bac6c4ccebbd1d9a7e3f90b9cfe125329fb6e"
+resolved_reference = "df67a6b95e214ecf14c33a5482bbd645884c5a85"
 
 [[package]]
 name = "humanfriendly"

python/energy_flow_proper/05_gaussian_two_baths/utilities.py

@@ -120,14 +120,14 @@ def plot_convergence(
     transform=lambda y: y,
     slice=None,
     linestyle="-",
-    bath: int = 0,
+    bath=None,
 ):
     label = label + ", " if (len(label) > 0) else ""
     slice = (0, -1) if not slice else slice
     for n, val, _ in y[slice[0] : slice[1]]:
         plt.plot(
             x,
-            transform(val[bath]),
+            transform(val[bath] if bath is not None else val),
             label=f"{label}n={n}",
             alpha=n / y[-1][0],
             linestyle="--",
@@ -135,8 +135,8 @@ def plot_convergence(
 
     ax.errorbar(
         x,
-        transform(y[-1][1][bath]),
-        yerr=y[-1][2][bath],
+        transform(y[-1][1][bath] if bath is not None else y[-1][1]),
+        yerr=y[-1][2][bath] if bath is not None else y[-1][2],
         ecolor="yellow",
         label=f"{label}n={y[-1][0]}",
         color="red",
@@ -160,7 +160,7 @@ def plot_diff_vs_sigma(
     ax.fill_between(
         x,
         0,
-        y[-1][2],
+        y[-1][2] / np.mean(np.abs(reference)),
         color=ecolor,
         label=fr"{label}$\sigma$",
     )
@@ -170,7 +170,7 @@ def plot_diff_vs_sigma(
         within = (diff < y[-1][2]).sum() / y[-1][2].size
         ax.plot(
             x,
-            diff,
+            diff / np.mean(np.abs(reference)),
             label=fr"{label}n={n}  $\Delta<\sigma = {within * 100}\%$",
             alpha=n / y[-1][0],
         )

python/energy_flow_proper/05_gaussian_two_baths/values/bcf_terms.tex

@@ -0,0 +1 @@
+\(K=12\)

python/energy_flow_proper/05_gaussian_two_baths/values/max_HO_level.tex

@@ -0,0 +1 @@
+\(9\)