update the subproject according to richards tips

2025-03-06 02:21:38 -05:00 · 2022-03-16 11:27:53 +01:00 · 2022-03-16 11:27:53 +01:00 · 1bd4d0ed95
commit 1bd4d0ed95
parent 171e04b742
17 changed files with 336 additions and 326 deletions
--- a/python/energy_flow_proper/01_zero_temperature/notebook.org
+++ b/python/energy_flow_proper/01_zero_temperature/notebook.org
@ -70,7 +70,12 @@ hi 500 integrate
 #+end_src

 #+RESULTS:
- : z16zulpojyom5baq0jum
+ :RESULTS:
+Linux ArLeenUX 5.16.8-zen1 x86_64
+ 15:57:39  up 1 day  5:47,  1 user,  load average: 1.55, 1.58, 1.91
+ impure  ~/D/P/U/m/m/p/e/01_zero_temperature  hi 500 integrate
+Loading the configuration from config.py. This might take a while... -
+:END:

 And there we go. It is better to run the above command in a
 vterm-session.
@ -180,57 +185,30 @@ Let's export some infos about the model to TeX.
 #+end_src

 #+RESULTS:
-| \(η=0.8\) | \(ω_c=2\) |
+| \(η=0.5\) | \(ω_c=2\) |


 ** Load the Data
 #+begin_src jupyter-python :results none
  from hopsflow import hopsflow, util
-  from hops.core.hierarchyLib import HI
+  from hops.core.hierarchyData import HIMetaData
 #+end_src

 Now we read the trajectory data.
 #+begin_src jupyter-python
  class result:
-      with HI(params, 500).get_data(read_only=True) as hd:
-          N = hd.samples
-          τ = hd.get_time()
-          ρ = hd.get_rho_t()
-          ψ_1 = np.array(hd.aux_states)[0:N]
-          ψ = np.array(hd.stoc_traj)[0:N]
+      hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
+      N = 8000
+      τ = hd.get_time()
+      ψ_1 = hd.aux_states
+      ψ = hd.stoc_traj
+
+
  fs.tex_value(result.N, prefix="N=", save="samples")
 #+end_src

 #+RESULTS:
-:RESULTS:
-# [goto error]
-#+begin_example
-  [0;31m---------------------------------------------------------------------------[0m
-  [0;31mRuntimeError[0m                              Traceback (most recent call last)
-  [0;32m/tmp/ipykernel_86433/975144810.py[0m in [0;36m<module>[0;34m[0m
-  [0;32m----> 1[0;31m [0;32mclass[0m [0mresult[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m      2[0m     [0;32mwith[0m [0mHI[0m[0;34m([0m[0mparams[0m[0;34m,[0m [0;36m500[0m[0;34m)[0m[0;34m.[0m[0mget_data[0m[0;34m([0m[0mread_only[0m[0;34m=[0m[0;32mTrue[0m[0;34m)[0m [0;32mas[0m [0mhd[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      3[0m         [0mN[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0msamples[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      4[0m         [0mτ[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0mget_time[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      5[0m         [0mρ[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0mget_rho_t[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/tmp/ipykernel_86433/975144810.py[0m in [0;36mresult[0;34m()[0m
-  [1;32m      2[0m     [0;32mwith[0m [0mHI[0m[0;34m([0m[0mparams[0m[0;34m,[0m [0;36m500[0m[0;34m)[0m[0;34m.[0m[0mget_data[0m[0;34m([0m[0mread_only[0m[0;34m=[0m[0;32mTrue[0m[0;34m)[0m [0;32mas[0m [0mhd[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      3[0m         [0mN[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0msamples[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m----> 4[0;31m         [0mτ[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0mget_time[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m      5[0m         [0mρ[0m [0;34m=[0m [0mhd[0m[0;34m.[0m[0mget_rho_t[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      6[0m         [0mψ_1[0m [0;34m=[0m [0mnp[0m[0;34m.[0m[0marray[0m[0;34m([0m[0mhd[0m[0;34m.[0m[0maux_states[0m[0;34m)[0m[0;34m[[0m[0;36m0[0m[0;34m:[0m[0mN[0m[0;34m][0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/nix/store/xyw00clxkmn0hgfw4lr7plyn61clhxyy-python3-3.9.9-env/lib/python3.9/site-packages/hops/core/hierarchyData.py[0m in [0;36mget_time[0;34m(self)[0m
-  [1;32m    757[0m [0;34m[0m[0m
-  [1;32m    758[0m         [0;32mif[0m [0;32mnot[0m [0mself[0m[0;34m.[0m[0mtime_set[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m--> 759[0;31m             [0;32mraise[0m [0mRuntimeError[0m[0;34m([0m[0;34m"can not get time, time has not been set yet."[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    760[0m         [0;32mreturn[0m [0mself[0m[0;34m.[0m[0mtime[0m[0;34m[[0m[0;34m:[0m[0;34m][0m  [0;31m# type: ignore[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    761[0m [0;34m[0m[0m
-
-  [0;31mRuntimeError[0m: can not get time, time has not been set yet.
-#+end_example
-:END:
+: \(N=8000\)

 #+begin_src jupyter-python
  with fs.hiro_style():
@ -249,7 +227,7 @@ Now we read the trajectory data.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/307ffa8ed20e93410c6858e7aed624294982fbb9.svg]]
+[[file:./.ob-jupyter/e798cfe7da7c268545439bad0355a288a3a70905.svg]]


 #+begin_src jupyter-python
@ -269,7 +247,7 @@ Now we read the trajectory data.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/90244163620b237d5a8d493923f8dd87ba0cb01a.svg]]
+[[file:./.ob-jupyter/9ac035f227cbb51dc2039c1de623d9ef48959fe4.svg]]

 ** Calculate System Energy
 Simple sanity check.
@ -292,7 +270,10 @@ Simple sanity check.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/d47cf2752cb838d4024b053730d50295639dca3c.svg]]
+:RESULTS:
+: 100% 7999/7999 [00:07<00:00, 1127.70it/s]
+[[file:./.ob-jupyter/ad0141ab03421b5c7eaa983380f0c0d834dee689.svg]]
+:END:

 The energy bleeds out of the system. We don't reach the steady state
 yet. Also we don't loose all the energy.
@ -311,7 +292,7 @@ The begin and and values of the system energy expectation are.
 #+end_src

 #+RESULTS:
-| 0.5 | -0.44770384926040235 |
+| 0.5 | -0.47523035186448653 |

 And the total energy lost is:
 #+begin_src jupyter-python
@ -319,7 +300,7 @@ And the total energy lost is:
 #+end_src

 #+RESULTS:
-: 0.9477038492604024
+: 0.9752303518644865

 We do start in the state.
 #+begin_src jupyter-python
@ -350,12 +331,12 @@ Now we can apply our tooling to one trajectory for testing.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/72d61a32bd4c505790dd3f70698a0525f693e078.svg]]
+[[file:./.ob-jupyter/c2feea9c73a6bcb0ae3436ecef6c8d576fc40f70.svg]]

 And now for all trajectories.
 #+begin_src jupyter-python
  full_flow = hopsflow.heat_flow_ensemble(
-      result.ψ, result.ψ_1, hf_system, result.N, every=result.N // 4, save="results/flow_1.npy"
+      result.ψ, result.ψ_1, hf_system, result.N, every=result.N // 2, save="results/flow_1.npy"
  )

  with fs.hiro_style():
@ -368,30 +349,46 @@ And now for all trajectories.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/ab447d51128e3ac19bfd435aabb0a71bd9143a12.svg]]
+:RESULTS:
+: 100% 7999/7999 [00:30<00:00, 261.68it/s]
+[[file:./.ob-jupyter/f890fbe7a6ba040867524afc80b591d93a7ac335.svg]]
+:END:

 We can integrate the energy change in the bath:
 #+begin_src jupyter-python
-  e_bath = util.integrate_array(-full_flow[-1][1], result.τ)
+  import scipy.integrate
+
+  e_bath = np.array([0] + [
+      scipy.integrate.simpson(-full_flow[-1][1][:i], result.τ[:i])
+      for i in range(1, len(result.τ))
+  ])
  plt.plot(result.τ, e_bath)
+  σ_e_bath = np.sqrt(np.array([0] + [
+      scipy.integrate.simpson(full_flow[-1][2][:i]**2, result.τ[:i])
+      for i in range(1, len(result.τ))
+  ])).real
+  plt.errorbar(result.τ, e_bath, yerr=σ_e_bath, ecolor="yellow")
 #+end_src

 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f21fed6e550> |
-[[file:./.ob-jupyter/cf8ab2488d384d9cdbc0ac97168442808be05939.svg]]
+: <ErrorbarContainer object of 3 artists>
+[[file:./.ob-jupyter/1ce48a60f2383456c071f34d54a58326eb5da07c.svg]]
 :END:

 ** Calculate the Interaction Energy
 First we calculate it from energy conservation.
 #+begin_src jupyter-python
  e_int = (1/2 - e_sys - e_bath).real
-  with fs.hiro_style():
-      plt.plot(result.τ, e_int)
+  σ_e_int = np.sqrt(σ_e_sys ** 2 + σ_e_bath ** 2).real
+  plt.errorbar(result.τ, e_int, yerr=σ_e_int, ecolor="yellow")
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/1fad72f04d79f936bf6fccffdb3fb63e1f9e13e3.svg]]
+:RESULTS:
+: <ErrorbarContainer object of 3 artists>
+[[file:./.ob-jupyter/943a9e709a62c32f9b4f938b6ef327228131e226.svg]]
+:END:

 And then from first principles:
 #+begin_src jupyter-python
@ -402,22 +399,21 @@ And then from first principles:

 #+RESULTS:
 :RESULTS:
-: 100% 999/999 [00:06<00:00, 166.19it/s]
-[[file:./.ob-jupyter/b2e4348b02a6a8db2e23c0b3d40357350bc2e9a6.svg]]
+: 100% 7999/7999 [00:40<00:00, 195.63it/s]
+[[file:./.ob-jupyter/2b6983dd07ce6cb4fec48b2695276d6931797541.svg]]
 :END:

 And both together:
 #+begin_src jupyter-python
  with fs.hiro_style():
-       plt.errorbar(result.τ, e_int_ex, yerr=σ_e_int_ex, ecolor="yellow", label="direct")
-       plt.xlabel("$τ$")
-       plt.ylabel(r"$\langle H_I\rangle$")
-       plt.plot(result.τ, e_int, label="from energy conservation")
-       plt.legend()
-
-       plt.gcf().set_size_inches(fs.get_figsize(239, 1, .8))
-       fs.export_fig("interaction")
+      plt.errorbar(result.τ, e_int, yerr=σ_e_int, label="from energy conservation", ecolor="yellow")
+      plt.errorbar(result.τ, e_int_ex, yerr=σ_e_int_ex, label="direct", ecolor="pink")
+      plt.gcf().set_size_inches(fs.get_figsize(239, 1, .8))
+      plt.legend()
+      plt.ylabel(r"$\langle H_I\rangle$")
+      plt.xlabel(r"$τ$")
+      fs.export_fig("interaction")
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/8fa5f86b96629477c4c599359c0668f6750d4b88.svg]]
+[[file:./.ob-jupyter/f49262f572c16add9e7f28811933da7443f98bf5.svg]]
--- a/python/energy_flow_proper/01_zero_temperature/poetry.lock
+++ b/python/energy_flow_proper/01_zero_temperature/poetry.lock
@ -288,7 +288,7 @@ numpy = ">=1.14.5"

 [[package]]
 name = "hops"
-version = "1.0"
+version = "1.0.1"
 description = "HOPS implementation of the TU-Dresden Theoretical Quantum Optics Group"
 category = "main"
 optional = false
@ -311,11 +311,11 @@ typer = "^0.4.0"
 type = "git"
 url = "git@gitlab.hrz.tu-chemnitz.de:s8896854--tu-dresden.de/hops.git"
 reference = "main"
-resolved_reference = "d47c4c5f733b462df7a046a6293a5db21c141479"
+resolved_reference = "b18016185092da74b00581c3925e912b3f59ef18"

 [[package]]
 name = "hopsflow"
-version = "1.0.1"
+version = "1.0.2"
 description = "Calculating open system bath energy changes with HOPS and analytically."
 category = "main"
 optional = false
@ -334,7 +334,7 @@ tqdm = "^4.62.3"
 type = "git"
 url = "https://github.com/vale981/hopsflow"
 reference = "main"
-resolved_reference = "8146c0dbd2844ee0b41cd07cd8ee63ab17607dc1"
+resolved_reference = "af3e744c828d69cedfbba6c7430a19520d9ae393"

 [[package]]
 name = "ipykernel"
@ -1118,7 +1118,7 @@ tests = ["pytest", "typeguard", "pygments", "littleutils"]

 [[package]]
 name = "stocproc"
-version = "1.0.2"
+version = "1.0.3"
 description = "Generate continuous time stationary stochastic processes from a given auto correlation function."
 category = "main"
 optional = false
@ -1135,7 +1135,7 @@ scipy = "^1.7.3"
 type = "git"
 url = "https://github.com/vale981/stocproc"
 reference = "master"
-resolved_reference = "2c806e3bdbb571970dec678f21a558e54b232a50"
+resolved_reference = "88b0f979a5c2d9e1a0cf32c33b73227f90d92671"

 [[package]]
 name = "terminado"
--- a/python/energy_flow_proper/01_zero_temperature/pyproject.toml
+++ b/python/energy_flow_proper/01_zero_temperature/pyproject.toml
@ -11,7 +11,7 @@ numpy = "^1.20.3"
 scipy = "^1.7.3"
 stocproc = { git = "https://github.com/vale981/stocproc" }
 hops = { git = "git@gitlab.hrz.tu-chemnitz.de:s8896854--tu-dresden.de/hops.git", branch="main" }
-hopsflow = { git = "https://github.com/vale981/hopsflow", branch="main" }
+hopsflow = { git = "https://github.com/vale981/hopsflow", branch="project_02_fix_backport" }
 matplotlib = "^3.5.0"
 jupyter = "^1.0.0"

--- a/python/energy_flow_proper/01_zero_temperature/results/energy.npy
+++ b/python/energy_flow_proper/01_zero_temperature/results/energy.npy
--- a/python/energy_flow_proper/01_zero_temperature/results/flow.npy
+++ b/python/energy_flow_proper/01_zero_temperature/results/flow.npy
--- a/python/energy_flow_proper/01_zero_temperature/values/bcf_scale.tex
+++ b/python/energy_flow_proper/01_zero_temperature/values/bcf_scale.tex
@ -1 +1 @@
-\(η=0.8\)
+\(η=0.5\)
--- a/python/energy_flow_proper/01_zero_temperature/values/samples.tex
+++ b/python/energy_flow_proper/01_zero_temperature/values/samples.tex
@ -1 +1 @@
-\(N=1000\)
+\(N=8000\)
--- a/python/energy_flow_proper/02_finite_temperature/notebook.org
+++ b/python/energy_flow_proper/02_finite_temperature/notebook.org
@ -130,17 +130,14 @@ Simple sanity check.

  with fs.hiro_style():
      plt.gcf().set_size_inches(fs.get_figsize(239, 1, .8))
-      plt.errorbar(result.τ[::100], e_sys.real[::100], yerr=σ_e_sys.real[::100], ecolor="yellow")
+      plt.errorbar(result.τ, e_sys.real, yerr=σ_e_sys.real, ecolor="yellow")
      plt.ylabel(r"$\langle H_S\rangle$")
      plt.xlabel(r"$τ$")
      fs.export_fig("system_energy")
 #+end_src

 #+RESULTS:
-:RESULTS:
-: 100% 9999/9999 [00:09<00:00, 1009.10it/s]
-[[file:./.ob-jupyter/6ca9a74ecaeb4a206f93f97921ee4a4bb212254f.svg]]
-:END:
+[[file:./.ob-jupyter/e3ff0f83ee48dc3aa72be9fb48e2687e16f8bfed.svg]]

 ** Calculate the Heat Flow
 Now let's calculate the heatflow. In this simple case it is engouh to
@ -172,12 +169,12 @@ Now we can apply our tooling to one trajectory for testing.

 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f8956fadfd0> |
-[[file:./.ob-jupyter/0c5b760af32ab0310202f083aaa6c54692a8f8dc.svg]]
+| <matplotlib.lines.Line2D | at | 0x7fee23b5ebb0> |
+[[file:./.ob-jupyter/0775a369159dac940ed20b03310e57d003ec7c42.svg]]
 :END:

 And now for all trajectories.
-#+begin_src jupyter-python :results none
+#+begin_src jupyter-python
  full_flow = hopsflow.heat_flow_ensemble(
      iter(result.ψ),
      iter(result.ψ_1),
@ -190,84 +187,16 @@ And now for all trajectories.
  )

  with fs.hiro_style():
-      fig, ax = fs.plot_convergence(result.τ, full_flow, transform=lambda y: -y)
+      fig, ax = fs.plot_convergence(result.τ, [full_flow[1], full_flow[-1]], transform=lambda y: -y)
      fig.set_size_inches(fs.get_figsize(239, 1, .8))
      ax.legend()
      ax.set_xlabel("$τ$")
      ax.set_ylabel("$-J$")
-      #fs.export_fig("flow", fig)
+      fs.export_fig("flow", fig)
 #+end_src

 #+RESULTS:
-:RESULTS:
-:   1% 90/9999 [00:21<39:31,  4.18it/s]Process ForkPoolWorker-173:
-: Process ForkPoolWorker-172:
-# [goto error]
-#+begin_example
-  [0;31m---------------------------------------------------------------------------[0m
-  [0;31mIndexError[0m                                Traceback (most recent call last)
-  [0;32m/nix/store/dn4fwp0yx6nsa85cr20cwvdmg64xwmcy-python3-3.9.9/lib/python3.9/multiprocessing/pool.py[0m in [0;36mnext[0;34m(self, timeout)[0m
-  [1;32m    852[0m             [0;32mtry[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m--> 853[0;31m                 [0mitem[0m [0;34m=[0m [0mself[0m[0;34m.[0m[0m_items[0m[0;34m.[0m[0mpopleft[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    854[0m             [0;32mexcept[0m [0mIndexError[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;31mIndexError[0m: pop from an empty deque
-
-  During handling of the above exception, another exception occurred:
-
-  [0;31mKeyboardInterrupt[0m                         Traceback (most recent call last)
-  [0;32m/tmp/ipykernel_94001/1430602402.py[0m in [0;36m<module>[0;34m[0m
-  [0;32m----> 1[0;31m full_flow = hopsflow.heat_flow_ensemble(
-  [0m[1;32m      2[0m     [0miter[0m[0;34m([0m[0mresult[0m[0;34m.[0m[0mψ[0m[0;34m)[0m[0;34m,[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      3[0m     [0miter[0m[0;34m([0m[0mresult[0m[0;34m.[0m[0mψ_1[0m[0;34m)[0m[0;34m,[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      4[0m     [0mhf_system[0m[0;34m,[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m      5[0m     [0mresult[0m[0;34m.[0m[0mN[0m[0;34m,[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/nix/store/c9msmd5k6clygvhbasl6871wb2ldg58c-python3-3.9.9-env/lib/python3.9/site-packages/hopsflow/hopsflow.py[0m in [0;36mheat_flow_ensemble[0;34m(ψ_0s, ψ_1s, params, N, therm_args, only_therm, **kwargs)[0m
-  [1;32m    357[0m         [0;32mraise[0m [0mValueError[0m[0;34m([0m[0;34m"Can't calculate only thermal part if therm_args are None."[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    358[0m [0;34m[0m[0m
-  [0;32m--> 359[0;31m     return util.ensemble_mean(
-  [0m[1;32m    360[0m         [0miter[0m[0;34m([0m[0mzip[0m[0;34m([0m[0mψ_0s[0m[0;34m,[0m [0mψ_1s[0m[0;34m,[0m [0mtherm_args[0m[0;34m[[0m[0;36m0[0m[0;34m][0m[0;34m)[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    361[0m         [0;32mif[0m [0mtherm_args[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/nix/store/c9msmd5k6clygvhbasl6871wb2ldg58c-python3-3.9.9-env/lib/python3.9/site-packages/hopsflow/util.py[0m in [0;36mensemble_mean[0;34m(arg_iter, function, N, const_args, const_kwargs, n_proc, every, save)[0m
-  [1;32m    221[0m         )
-  [1;32m    222[0m [0;34m[0m[0m
-  [0;32m--> 223[0;31m         [0;32mfor[0m [0mres[0m [0;32min[0m [0mtqdm[0m[0;34m([0m[0mresult_iter[0m[0;34m,[0m [0mtotal[0m[0;34m=[0m[0;34m([0m[0mN[0m [0;34m-[0m [0;36m1[0m[0;34m)[0m [0;32mif[0m [0mN[0m [0;32melse[0m [0;32mNone[0m[0;34m)[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    224[0m             [0maggregate[0m[0;34m.[0m[0mupdate[0m[0;34m([0m[0mres[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    225[0m [0;34m[0m[0m
-
-  [0;32m/nix/store/c9msmd5k6clygvhbasl6871wb2ldg58c-python3-3.9.9-env/lib/python3.9/site-packages/tqdm/std.py[0m in [0;36m__iter__[0;34m(self)[0m
-  [1;32m   1178[0m [0;34m[0m[0m
-  [1;32m   1179[0m         [0;32mtry[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m-> 1180[0;31m             [0;32mfor[0m [0mobj[0m [0;32min[0m [0miterable[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m   1181[0m                 [0;32myield[0m [0mobj[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m   1182[0m                 [0;31m# Update and possibly print the progressbar.[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/nix/store/dn4fwp0yx6nsa85cr20cwvdmg64xwmcy-python3-3.9.9/lib/python3.9/multiprocessing/pool.py[0m in [0;36m<genexpr>[0;34m(.0)[0m
-  [1;32m    446[0m                     [0mresult[0m[0;34m.[0m[0m_set_length[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    447[0m                 ))
-  [0;32m--> 448[0;31m             [0;32mreturn[0m [0;34m([0m[0mitem[0m [0;32mfor[0m [0mchunk[0m [0;32min[0m [0mresult[0m [0;32mfor[0m [0mitem[0m [0;32min[0m [0mchunk[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    449[0m [0;34m[0m[0m
-  [1;32m    450[0m     def apply_async(self, func, args=(), kwds={}, callback=None,
-
-  [0;32m/nix/store/dn4fwp0yx6nsa85cr20cwvdmg64xwmcy-python3-3.9.9/lib/python3.9/multiprocessing/pool.py[0m in [0;36mnext[0;34m(self, timeout)[0m
-  [1;32m    856[0m                     [0mself[0m[0;34m.[0m[0m_pool[0m [0;34m=[0m [0;32mNone[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    857[0m                     [0;32mraise[0m [0mStopIteration[0m [0;32mfrom[0m [0;32mNone[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m--> 858[0;31m                 [0mself[0m[0;34m.[0m[0m_cond[0m[0;34m.[0m[0mwait[0m[0;34m([0m[0mtimeout[0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    859[0m                 [0;32mtry[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    860[0m                     [0mitem[0m [0;34m=[0m [0mself[0m[0;34m.[0m[0m_items[0m[0;34m.[0m[0mpopleft[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;32m/nix/store/dn4fwp0yx6nsa85cr20cwvdmg64xwmcy-python3-3.9.9/lib/python3.9/threading.py[0m in [0;36mwait[0;34m(self, timeout)[0m
-  [1;32m    310[0m         [0;32mtry[0m[0;34m:[0m    [0;31m# restore state no matter what (e.g., KeyboardInterrupt)[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    311[0m             [0;32mif[0m [0mtimeout[0m [0;32mis[0m [0;32mNone[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-  [0;32m--> 312[0;31m                 [0mwaiter[0m[0;34m.[0m[0macquire[0m[0;34m([0m[0;34m)[0m[0;34m[0m[0;34m[0m[0m
-  [0m[1;32m    313[0m                 [0mgotit[0m [0;34m=[0m [0;32mTrue[0m[0;34m[0m[0;34m[0m[0m
-  [1;32m    314[0m             [0;32melse[0m[0;34m:[0m[0;34m[0m[0;34m[0m[0m
-
-  [0;31mKeyboardInterrupt[0m:
-#+end_example
-:END:
+[[file:./.ob-jupyter/4baa612c7e7201f51f74886c2f74caa2a56383a8.svg]]


 We can integrate the energy change in the bath:
@ -289,7 +218,7 @@ We can integrate the energy change in the bath:
 #+RESULTS:
 :RESULTS:
 : <ErrorbarContainer object of 3 artists>
-[[file:./.ob-jupyter/a65a13edc562621973be18c8fbd4271bf3afa733.svg]]
+[[file:./.ob-jupyter/fd2493973fc2770c764cc73be94f2513def80b8e.svg]]
 :END:

 ** Calculate the Interaction Energy
@ -303,7 +232,7 @@ First we calculate it from energy conservation.
 #+RESULTS:
 :RESULTS:
 : <ErrorbarContainer object of 3 artists>
-[[file:./.ob-jupyter/5950975b59e103671be041250987da84adfdca83.svg]]
+[[file:./.ob-jupyter/64d9fac11c4d3356642e6e585b0ceb4b22d7e8ad.svg]]
 :END:

 And then from first principles:
@ -319,7 +248,6 @@ And then from first principles:
 #+end_src

 #+RESULTS:
-: 100% 499/499 [02:09<00:00,  3.86it/s]

 And both together:
 #+begin_src jupyter-python
@ -330,11 +258,11 @@ And both together:
      plt.legend()
      plt.ylabel(r"$\langle H_I\rangle$")
      plt.xlabel(r"$τ$")
-      #fs.export_fig("interaction")
+      fs.export_fig("interaction")
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/ecbcd0497b02327dd0a678f193edc9c3d9e7c553.svg]]
+[[file:./.ob-jupyter/48b4df7cf090fe5c307194ad23cdb51da1538247.svg]]

 Seems to work :P.
 #+begin_src jupyter-python
@ -344,10 +272,10 @@ Seems to work :P.

 #+RESULTS:
 :RESULTS:
-: /nix/store/maj240ris2x1a0r1wz26m0m785bc0xn9-python3-3.9.9-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
+: /nix/store/c9msmd5k6clygvhbasl6871wb2ldg58c-python3-3.9.9-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
 :   return np.asarray(x, float)
-| <matplotlib.lines.Line2D | at | 0x7f15e986fbb0> |
-[[file:./.ob-jupyter/8161e70494da51689f01e45a2beb0d57798a7436.svg]]
+| <matplotlib.lines.Line2D | at | 0x7fee23c017f0> |
+[[file:./.ob-jupyter/2bb493ebb9eb6bc714253b115f92b26457542e60.svg]]
 :END:

 * Close the Data File
@ -355,3 +283,8 @@ We need to release the hold on the file.
 #+begin_src jupyter-python :results none
 result.hd.close()
 #+end_src
+
+* Observations
+- convergence slower
+- error estimate is of course too small
+- energy change not as steep for smaller \(\omega_c\)
--- a/python/energy_flow_proper/04_gaussian_nonzero/comparison_with_hops.org
+++ b/python/energy_flow_proper/04_gaussian_nonzero/comparison_with_hops.org
@ -136,12 +136,10 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-: INFO:root:Loading cache from: results/system_energy_sandwhich_operator_10000_2000_8b3d580ffb3f80be6d0d9ff18a797be9c1efea107513b53e70df4b0dde8965b3.npy
-

 #+begin_src jupyter-python
  with fs.hiro_style():
-      fig, ax = fs.plot_convergence(result.τ, e_sys, transform=lambda y: y.real)
+      fig, ax = fs.plot_convergence(result.τ, [e_sys[2], e_sys[-1]], transform=lambda y: y.real)
      fig.set_size_inches(fs.get_figsize(239, 1, .8))
      ax.set_ylabel(r"$\langle H_S\rangle$")
      ax.set_xlabel(r"$τ$")
@ -150,7 +148,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/8034455491647f815b28b325d612c61b8a2f416b.svg]]
+[[file:./.ob-jupyter/99810f19c992a7cb6462154dcdca10cbe5cfc852.svg]]

 #+begin_src jupyter-python
  with fs.hiro_style():
@ -162,7 +160,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/45b4a504c688f20f85940b4001d83f53337c120a.svg]]
+[[file:./.ob-jupyter/3b6a8d1c8657b7f9d82cb7c101ab2769f6af0aba.svg]]

 * Heat Flow
 #+begin_src jupyter-python
@ -181,13 +179,13 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
      hf_system,
      result.N,
      (iter(result.seeds), hf_therm),
-      every=int(result.N / 5),
+      every=int(result.N / 2),
      save="heat_flow"
  )
 #+end_src

 #+RESULTS:
-: INFO:root:Loading cache from: results/heat_flow__heat_flow_ensemble_body_10000_2000_a787815f13e7fed6d3919548b33e06aea41908005c5f1c5886efc6103d44a0a5.npy
+

 #+begin_src jupyter-python
  with fs.hiro_style():
@ -201,7 +199,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/c94da92671284fb92f8f946bc9f57b822b953bf4.svg]]
+[[file:./.ob-jupyter/944b6da2b99f17761e2b016110ad63deb9c2559e.svg]]



@ -249,7 +247,7 @@ LGTM.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/6954e40a65e27a26a895bde9d7a919887f7bfb3c.svg]]
+[[file:./.ob-jupyter/77f10f7a79628b528c1510a2f5cc5d6be07bc6bc.svg]]

 #+begin_src jupyter-python :results none
  sys = gf.SystemParams(Ω=Ω, η=system.bcf_scale, α_0=α_0)
@ -264,7 +262,7 @@ LGTM.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/38ff6b1cf249e5618cc74b5e83c0e0bf58e55cf3.svg]]
+[[file:./.ob-jupyter/625a0c9d68d7d0759d58724a06837f0ad1201d98.svg]]

 * Compare
 #+begin_src jupyter-python
@ -279,7 +277,7 @@ LGTM.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/c9630ccb572710a1f0be1d22c0e74efe93dc5418.svg]]
+[[file:./.ob-jupyter/91f59ad05009d5cc07f4b905b664db9304d7d1a8.svg]]

 #+begin_src jupyter-python
  with fs.hiro_style():
@ -293,7 +291,7 @@ LGTM.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/9ee096ff0f380e5ef202106379f649e49f5e3a56.svg]]
+[[file:./.ob-jupyter/b717fd1b1e15f33ae40ed869ce4c2676874e8bad.svg]]

 * Cleanup
 #+begin_src jupyter-python :results none
--- a/python/energy_flow_proper/04_gaussian_nonzero/comparison_with_hops.py
+++ b/python/energy_flow_proper/04_gaussian_nonzero/comparison_with_hops.py
@ -35,7 +35,7 @@ e_sys = util.operator_expectation_ensemble(
 )

 with fs.hiro_style():
-    fig, ax = fs.plot_convergence(result.τ, e_sys, transform=lambda y: y.real)
+    fig, ax = fs.plot_convergence(result.τ, [e_sys[2], e_sys[-1]], transform=lambda y: y.real)
    fig.set_size_inches(fs.get_figsize(239, 1, .8))
    ax.set_ylabel(r"$\langle H_S\rangle$")
    ax.set_xlabel(r"$τ$")
@ -64,7 +64,7 @@ full_flow = hopsflow.heat_flow_ensemble(
    hf_system,
    result.N,
    (iter(result.seeds), hf_therm),
-    every=int(result.N / 5),
+    every=int(result.N / 2),
    save="heat_flow"
 )

--- a/python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.org
+++ b/python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.org
@ -573,10 +573,10 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/a8f04e9a2da787db38bb667c3ff290a1e782348e.svg]]
+[[file:./.ob-jupyter/236006fa4bf555c102b8c09a3365c8286a80c22e.svg]]

 #+begin_src jupyter-python :results none
-  gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2, root_tol=1e-5)
+  gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2, root_tol=1e-5/2)
 #+end_src

 * Load Hops Data
@ -602,11 +602,29 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
      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")
+  ), fs.tex_value(
+      Ω, prec=0, save="omega", prefix="Ω="
+  ), fs.tex_value(
+      max_HO_level_set, save="max_HO_level"
+  ), fs.tex_value(
+      γ, prec=1, save="gamma", prefix="γ="
+  )
 #+end_src

 #+RESULTS:
-| \(η=0.2\) | \(ω_c=2\) | \(T=0.6\) | \(Ω=1\) | \(9\) |
+| \(η=0.2\) | \(ω_c=2\) | \(T=0.6\) | \(Ω=1\) | \(9\) | \(γ=0.5\) |
+
+#+begin_src jupyter-python
+  with fs.hiro_style():
+      plt.plot(result.τ, np.linalg.norm(result.ψ[0], axis=1) - 1)
+      plt.xlabel("$τ$")
+      plt.gcf().set_size_inches(fs.get_figsize(398, 0.8, 0.3))
+      plt.ylabel("$||ψ|| - 1$")
+      fs.export_fig("norm")
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/ac8c3d83dd1b8577242c719db3e34a1d9cd3b72a.svg]]

 * System Energy
 #+begin_src jupyter-python
@ -630,7 +648,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/d35e70d489882edecd6565bf17fed76bdb242b60.svg]]
+[[file:./.ob-jupyter/20bbab9b2d03841eee868b45d588ddb14625a35f.svg]]

 Now with Gaussflow.
 #+begin_src jupyter-python
@ -642,7 +660,7 @@ Now with Gaussflow.

 #+begin_src jupyter-python
  with fs.hiro_style():
-      fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real)
+      fig, ax = fs.plot_convergence(result.τ, [e_sys[2], e_sys[-1]], 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)
@ -655,12 +673,12 @@ Now with Gaussflow.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/c574ce2d97056f5772a61da13eca6d4b8ac5ab7a.svg]]
+[[file:./.ob-jupyter/c1a26126332aecd8d2491637673e40e0780139ab.svg]]


 #+begin_src jupyter-python
  with fs.hiro_style():
-      fig, ax = fs.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
+      fig, ax = fs.plot_diff_vs_sigma(result.τ, [e_sys[2], e_sys[-1]], energy_gf)
      fig.set_size_inches(fs.get_figsize(239, 1, .8))
      ax.set_xlabel(r"$τ$")
      ax.legend()
@ -668,7 +686,7 @@ Now with Gaussflow.
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/ba0a30288d48d6da5e22938d414e4e6bffdec659.svg]]
+[[file:./.ob-jupyter/ce8435b75b0ea1e08dab50c9a1e048b7f3e78eb7.svg]]

 * Thermal State
 #+begin_src jupyter-python
@ -775,8 +793,8 @@ Interestingly in the middle.

 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f060bd82700> |
-[[file:./.ob-jupyter/494f4a872592e6928904bfe575a46b0f7a14a40a.svg]]
+| <matplotlib.lines.Line2D | at | 0x7f9f74830c10> |
+[[file:./.ob-jupyter/7ce8477df0437882a308254c616313e2876a2c21.svg]]
 :END:

 #+begin_src jupyter-python
@ -804,8 +822,8 @@ Interestingly in the middle.

 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/66e584bd18b5548fe7c1baf5fc243539e23e6841.svg]]
-[[file:./.ob-jupyter/eaa35bd6cdea4956b1360c0e9bda4212a45cc8a1.svg]]
+[[file:./.ob-jupyter/cf83f775e74276795856cc029aa19f48fd88480b.svg]]
+[[file:./.ob-jupyter/b274f2700432992dc2fbc623ac6c0687fc878eac.svg]]
 :END:


@ -816,57 +834,62 @@ Interestingly in the middle.
 #+RESULTS:

 #+begin_src jupyter-python
+  fig, ax = plt.subplots()
+  fig.set_size_inches(fs.get_figsize(239, 1, 0.8))
  for bath in range(2):
      with fs.hiro_style():
-          fig, ax = ut.plot_convergence(
-              result.τ, full_flow[::5], transform=lambda y: -y, bath=bath
+          fs.plot_convergence(
+              result.τ,
+              [full_flow[10], full_flow[-1]],
+              transform=lambda y: -y,
+              bath=bath,
+              ax=ax,
+              label=f"bath {bath + 1}"
          )
-          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",
+              label=f"analytic ({bath + 1})",
              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}")
+  fs.export_fig(f"flow")
 #+end_src

 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/0bdb5ce9fa4f4929667ce0d8398dbc3b6c8441a2.svg]]
-[[file:./.ob-jupyter/834cd5ce6f522f761ae1e5214a3d71c6a18c0918.svg]]
+[[file:./.ob-jupyter/15a8245249f94a0fcb1142488c0353a261f271ad.svg]]
+[[file:./.ob-jupyter/d6fada73d9df039fac276a4b1f1417ac9b5e7614.svg]]
 :END:

-
 #+begin_src jupyter-python
+  position = ["left", "right"]
  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))
+        with fs.hiro_style():
+            fig, ax = fs.plot_diff_vs_sigma(
+                result.τ,
+                [full_flow[10], full_flow[-1]],
+                -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}")
+            ax.legend()
+            ax.set_xlabel("$τ$")
+            ax.legend()
+            ax.set_title(f"Bath {bath + 1} ({position[bath]}, $T={params.SysP.__non_key__['T'][bath]}$)")
+            fs.export_fig(f"flow_conv_{bath}")
 #+end_src

 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/c2efc39704e07f80cb3bb4b1f8655e2659e81fa1.svg]]
-[[file:./.ob-jupyter/dca19c821f5a6f25aae2d7168bbcb19f4820faa0.svg]]
+[[file:./.ob-jupyter/06d4fce4b6ec92df105f7246128f6edb8c02342f.svg]]
+[[file:./.ob-jupyter/60916d79f1b29a8f642387cf43e83ad870bd34d2.svg]]
 :END:

 * Thoughts
--- a/python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.py
+++ b/python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.py
@ -267,7 +267,7 @@ with fs.hiro_style():

 # 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)
+gf_params = gf.SystemParams(Ω=Ω, Λ=Λ, η=params.SysP.bcf_scale, γ=γ, α_0=[α_0] * 2, root_tol=1e-5/2)

 class result:
    hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
@ -285,7 +285,20 @@ 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")
+), fs.tex_value(
+    Ω, prec=0, save="omega", prefix="Ω="
+), fs.tex_value(
+    max_HO_level_set, save="max_HO_level"
+), fs.tex_value(
+    γ, prec=1, save="gamma", prefix="γ="
+)
+
+with fs.hiro_style():
+    plt.plot(result.τ, np.linalg.norm(result.ψ[0], axis=1) - 1)
+    plt.xlabel("$τ$")
+    plt.gcf().set_size_inches(fs.get_figsize(398, 0.8, 0.3))
+    plt.ylabel("$||ψ|| - 1$")
+    fs.export_fig("norm")

 e_sys = util.operator_expectation_ensemble(
    iter(result.ψ),
@ -304,7 +317,7 @@ C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=
 energy_gf = C.system_energy(result.τ)

 with fs.hiro_style():
-    fig, ax = fs.plot_convergence(result.τ, e_sys, transform=np.real)
+    fig, ax = fs.plot_convergence(result.τ, [e_sys[2], e_sys[-1]], 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)
@ -316,7 +329,7 @@ with fs.hiro_style():
    # plt.ylim(1.5, 1.6)

 with fs.hiro_style():
-    fig, ax = fs.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
+    fig, ax = fs.plot_diff_vs_sigma(result.τ, [e_sys[2], e_sys[-1]], energy_gf)
    fig.set_size_inches(fs.get_figsize(239, 1, .8))
    ax.set_xlabel(r"$τ$")
    ax.legend()
@ -383,39 +396,45 @@ with ut.hiro_style():

 exact_flow = [C.flow(result.τ, 0, False), C.flow(result.τ, 1, False)]

+fig, ax = plt.subplots()
+fig.set_size_inches(fs.get_figsize(239, 1, 0.8))
 for bath in range(2):
    with fs.hiro_style():
-        fig, ax = ut.plot_convergence(
-            result.τ, full_flow[::5], transform=lambda y: -y, bath=bath
+        fs.plot_convergence(
+            result.τ,
+            [full_flow[10], full_flow[-1]],
+            transform=lambda y: -y,
+            bath=bath,
+            ax=ax,
+            label=f"bath {bath + 1}"
        )
-        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",
+            label=f"analytic ({bath + 1})",
            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}")
+fs.export_fig(f"flow")

+position = ["left", "right"]
 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))
+      with fs.hiro_style():
+          fig, ax = fs.plot_diff_vs_sigma(
+              result.τ,
+              [full_flow[10], full_flow[-1]],
+              -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}")
+          ax.legend()
+          ax.set_xlabel("$τ$")
+          ax.legend()
+          ax.set_title(f"Bath {bath + 1} ({position[bath]}, $T={params.SysP.__non_key__['T'][bath]}$)")
+          fs.export_fig(f"flow_conv_{bath}")
--- a/python/energy_flow_proper/05_gaussian_two_baths/longhopsidea.org
+++ b/python/energy_flow_proper/05_gaussian_two_baths/longhopsidea.org
--- a/python/energy_flow_proper/05_gaussian_two_baths/values/gamma.tex
+++ b/python/energy_flow_proper/05_gaussian_two_baths/values/gamma.tex
@ -0,0 +1 @@
+\(γ=0.5\)
--- a/python/energy_flow_proper/06_two_qubit_first_experiments/figsaver.py
+++ b/python/energy_flow_proper/06_two_qubit_first_experiments/figsaver.py
@ -0,0 +1 @@
+../figsaver.py
--- a/python/energy_flow_proper/06_two_qubit_first_experiments/test_convergence.org
+++ b/python/energy_flow_proper/06_two_qubit_first_experiments/test_convergence.org
@ -18,7 +18,7 @@ Use punishing parameters.
  models = [
      TwoQubitModel(
          γ=1,
-          δ=[0.6, 0.6],
+          δ=[0.4, 0.4],
          T=[0, 0],
          t_max=10,
          influence_tolerance=kf,
@ -26,10 +26,10 @@ Use punishing parameters.
          ω_c=[1, 1],
          truncation_scheme="bath_memory",
      )
-      for kf in [0.1, 0.08, 0.05, 0.03, 0.01] #[0.001, 0.002, 0.003, 0.005]
+      for kf in [0.1, 0.08, 0.05, 0.03, 0.01, 0.005, 0.001] #[0.001, 0.002, 0.003, 0.005]
  ]

-  for k_fac in np.linspace(1, 3, 10):
+  for k_fac in np.linspace(1, 2, 5):
      comparison = models[0].copy()
      comparison.k_fac = [k_fac] * 2
      comparison.truncation_scheme = "power"
@ -66,16 +66,13 @@ And now we get ourselves a single trajectory for each ~k_fac~.
 : 0it [00:00, ?it/s]
 : 0it [00:00, ?it/s]
 : 0it [00:00, ?it/s]
+:   0% 0/1 [00:00<?, ?it/s]
+:   0% 0/1 [00:00<?, ?it/s]
+: 0it [00:00, ?it/s]
+: 0it [00:00, ?it/s]
+: 0it [00:00, ?it/s]
+: 0it [00:00, ?it/s]
 : 0it [00:00, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
-:   0% 0/1 [00:00<?, ?it/s]
 : 0it [00:00, ?it/s]
 :END:

@ -89,15 +86,18 @@ And now we get ourselves a single trajectory for each ~k_fac~.


 #+begin_src jupyter-python
-  with ut.hiro_style():
-      fig, axs = ut.plot_trajectory_convergence(τ, trajectories)
+  import figsaver as fs
+  with fs.hiro_style():
+        fig, axs = ut.plot_trajectory_convergence(τ, trajectories)
+        plt.gcf().set_size_inches(fs.get_figsize(398, .8, .6))
+        axs.legend(prop={'size': 4})

-      #fig.set_size_inches((12,8))
+        #fig.set_size_inches((12,8))

-      plt.yscale('log')
-      fig.set_size_inches(10,8)
-      fig
+        plt.yscale('log')
+
+        fs.export_fig("convergence")
 #+end_src

 #+RESULTS:
-[[file:./.ob-jupyter/393ecbe7d268dc053e18395c0f4cc33172558a13.svg]]
+[[file:./.ob-jupyter/abba1e1105a25c94bd0474e60fdbf3cea50071a0.svg]]
--- a/python/energy_flow_proper/06_two_qubit_first_experiments/test_convergence.py
+++ b/python/energy_flow_proper/06_two_qubit_first_experiments/test_convergence.py
@ -0,0 +1,57 @@
+from two_qubit_model import TwoQubitModel
+import qutip as qt
+from hops.core.integration import HOPSSupervisor
+import utilities as ut
+from tqdm import tqdm
+import numpy as np
+
+models = [
+    TwoQubitModel(
+        γ=1,
+        δ=[0.4, 0.4],
+        T=[0, 0],
+        t_max=10,
+        influence_tolerance=kf,
+        bcf_terms=[4, 4],
+        ω_c=[1, 1],
+        truncation_scheme="bath_memory",
+    )
+    for kf in [0.1, 0.08, 0.05, 0.03, 0.01, 0.005, 0.001] #[0.001, 0.002, 0.003, 0.005]
+]
+
+for k_fac in np.linspace(1, 2, 5):
+    comparison = models[0].copy()
+    comparison.k_fac = [k_fac] * 2
+    comparison.truncation_scheme = "power"
+
+    models.append(comparison)
+
+comparison = models[0].copy()
+comparison.k_max = 9
+comparison.truncation_scheme = "simplex"
+models.append(comparison)
+
+configs = [model.hops_config for model in models]
+
+supervisors = [HOPSSupervisor(config, number_of_samples=1) for config in configs]
+
+for supervisor, model in zip(supervisors, models):
+    supervisor.integrate_single_process()
+
+trajectories = []
+for supervisor, model, config in zip(supervisors, models, configs):
+    with supervisor.get_data(True) as data:
+        τ = data.get_time()
+        trajectories.append((model, config, data.get_stoc_traj(0)))
+
+import figsaver as fs
+with fs.hiro_style():
+      fig, axs = ut.plot_trajectory_convergence(τ, trajectories)
+      plt.gcf().set_size_inches(fs.get_figsize(398, .8, .6))
+      axs.legend(prop={'size': 4})
+
+      #fig.set_size_inches((12,8))
+
+      plt.yscale('log')
+
+      fs.export_fig("convergence")