updates in project 05

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.org

@@ -204,7 +204,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
-  make_config(
+  params = make_config(
       max_HO_level=8,
       bcf_terms=4,
       t_max=50,
@@ -344,8 +344,8 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 
 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/eb3f1978b8663a27a45d24124456861cbb22b4c1.svg]]
-[[file:./.ob-jupyter/fd3c2dabd7d087e246a11747045d85fe0947099b.svg]]
+[[file:./.ob-jupyter/7860e5186a807c5b41a3d8b01e841284e959ff4b.svg]]
+[[file:./.ob-jupyter/23c7abbd150839cddfec7ce4f096451f25b6fa04.svg]]
 :END:
 
 #+begin_src jupyter-python :results none
@@ -367,69 +367,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-:RESULTS:
-: hi_key_bin could not be read from file data/_c/data_cb704de6c6755c80309dd26a16aadc98_1.h5
-# [goto error]
-#+begin_example
-  ---------------------------------------------------------------------------
-  OSError                                   Traceback (most recent call last)
-  Input In [322], in <module>
-  ----> 1 class result:
-        2     hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
-        3     N = hd.samples
-
-  Input In [322], in result()
-        1 class result:
-  ----> 2     hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
-        3     N = hd.samples
-        4     seeds = hd.rng_seed
-
-  File /nix/store/403l2yqlsk1rs7kdyjk3qhdakdwq6a28-python3-3.9.9-env/lib/python3.9/site-packages/hops/core/hierarchy_data.py:1240, in HIMetaData.get_HIData(self, key, read_only, overwrite_key)
-     1231 def get_HIData(
-     1232     self, key: HIParams, read_only: bool = False, overwrite_key: bool = False
-     1233 ) -> HIData:
-     1234     """Returns a :any:`HIData` instance initialized with ``key`` and an
-     1235     auto-generated file name.
-     1236
-     1237     For the arguments see :any:`hops.core.hierarchy_data.HIData`.
-     1238     """
-  -> 1240     hdf5_name, bin_key, hashed_key = self.get_HIData_fname(
-     1241         key, ret_bin_and_hash=True
-     1242     )
-     1244     assert isinstance(bin_key, bytes)
-     1245     return HIData(
-     1246         str(self.path / hdf5_name),
-     1247         read_only=read_only,
-     (...)
-     1251         overwrite_key=overwrite_key,
-     1252     )
-
-  File /nix/store/403l2yqlsk1rs7kdyjk3qhdakdwq6a28-python3-3.9.9-env/lib/python3.9/site-packages/hops/core/hierarchy_data.py:1217, in HIMetaData.get_HIData_fname(self, key, ret_bin_and_hash)
-     1213 with h5py.File(
-     1214     str(self.path / hdf5_name), "r", swmr=True, libver="latest"
-     1215 ) as h5File:
-     1216     hkb = h5File["hi_key_bin"]
-  -> 1217     hkb = bytearray(hkb[:])  # type: ignore
-     1218     if hkb == bin_key:  # this is the file we are looking for!
-     1219         break
-
-  File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()
-
-  File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()
-
-  File /nix/store/403l2yqlsk1rs7kdyjk3qhdakdwq6a28-python3-3.9.9-env/lib/python3.9/site-packages/h5py/_hl/dataset.py:710, in Dataset.__getitem__(self, args, new_dtype)
-      708 if self._fast_read_ok and (new_dtype is None):
-      709     try:
-  --> 710         return self._fast_reader.read(args)
-      711     except TypeError:
-      712         pass  # Fall back to Python read pathway below
-
-  File h5py/_selector.pyx:366, in h5py._selector.Reader.read()
-
-  OSError: [Errno 5] Can't read data (file read failed: time = Tue Feb 15 14:17:54 2022
-  , filename = 'data/_c/data_cb704de6c6755c80309dd26a16aadc98_1.h5', file descriptor = 55, errno = 5, error message = 'Input/output error', buf = 0x9d57510, total read size = 203999, bytes this sub-read = 203999, bytes actually read = 18446744073709551615, offset = 0)
-#+end_example
-:END:
+: 10000
 
 * System Energy
 #+begin_src jupyter-python
@@ -445,7 +383,7 @@ 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_1000_d0ed047dab9d3ae1be74b1e5b17755f0c1a41bd611212b504a35555740e6bf85.npy
+: WARNING:root:Loading cache from: results/system_energy_sandwhich_operator_10000_1000_a34265c8f80fbc84d6a3b074cffc30d9abe302353134db3cf684ee0c0fded8ea.npy
 
 
 #+begin_src jupyter-python
@@ -454,7 +392,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/5bd904b9184abb5bc91b85d68a42bebe5ba1f5ae.svg]]
+[[file:./.ob-jupyter/222be3d5cd5d5b527ea800e08bf73b2ee0de89ef.svg]]
 
 Now with Gaussflow.
 #+begin_src jupyter-python
@@ -474,7 +412,7 @@ Now with Gaussflow.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/9975e4a825c9842e39abf68b76523672581c7c2d.svg]]
+[[file:./.ob-jupyter/78cf68864d651dd1405f0dfa80f7c831a7bde996.svg]]
 
 
 #+begin_src jupyter-python
@@ -486,7 +424,7 @@ Now with Gaussflow.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/2f82b9edaf51d90e93a9cbc50fda894b9a8d96fb.svg]]
+[[file:./.ob-jupyter/4e49054f9b64e79c49a9c4af01b2daf9ffa0124d.svg]]
 
 * Thermal State
 #+begin_src jupyter-python
@@ -536,3 +474,95 @@ Now with Gaussflow.
 
 #+RESULTS:
 : 0.07726133069400465
+
+* Heat Flow
+#+begin_src jupyter-python :results none
+  from hopsflow import hopsflow as hf
+#+end_src
+
+#+begin_src jupyter-python
+  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)
+  ξ_1 = params.EtaTherm[0]
+  ξ_1.set_scale(params.SysP.bcf_scale[0])
+  ξ_1.calc_deriv = True
+  therm_params = hf.ThermalParams([ξ_1, None], result.τ)
+  therm_run = hf.ThermalRunParams(therm_params, result.seeds[0])
+#+end_src
+
+#+RESULTS:
+
+#+begin_src jupyter-python :results none
+  flow = hf.flow_trajectory(run, sys, therm_run)
+#+end_src
+
+
+#+begin_src jupyter-python
+  plt.plot(result.τ, flow[0])
+  plt.plot(result.τ, flow[1])
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f6f9892ae80> |
+[[file:./.ob-jupyter/1d3f0e09cc8a1373f485a5feded2f90d6a589871.svg]]
+:END:
+
+#+begin_src jupyter-python
+  full_flow = hopsflow.heat_flow_ensemble(
+      iter(result.ψ),
+      iter(result.ψ_1),
+      sys,
+      500,
+      (iter(result.seeds), therm_params),
+      every=int(result.N / 5),
+      save="heat_flow"
+  )
+#+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
+
+#+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.legend()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/bac8b6fe9327df5899f5d317aa4af42cfb092b43.svg]]
+
+
+#+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])
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f6f70916670> |
+[[file:./.ob-jupyter/755f01900a0ee80481360bbf33eea3c62c2c717d.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()
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/9620d5702095804e22f9572a2266d3e4a3699180.svg]]

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.py

@@ -175,3 +175,50 @@ thermal_e(v, 1.0 / (params.SysP.__non_key__["T"][0]))
 ρ_therm = sc.linalg.expm(-params.SysP.H_sys/(params.SysP.__non_key__["T"][0]))
 ρ_therm = ρ_therm / np.trace(ρ_therm)
 trace_norm(result.ρ[-1] - ρ_therm)
+
+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)
+ξ_1 = params.EtaTherm[0]
+ξ_1.set_scale(params.SysP.bcf_scale[0])
+ξ_1.calc_deriv = True
+therm_params = hf.ThermalParams([ξ_1, None], result.τ)
+therm_run = hf.ThermalRunParams(therm_params, result.seeds[0])
+
+flow = hf.flow_trajectory(run, sys, therm_run)
+
+plt.plot(result.τ, flow[0])
+plt.plot(result.τ, flow[1])
+
+full_flow = hopsflow.heat_flow_ensemble(
+    iter(result.ψ),
+    iter(result.ψ_1),
+    sys,
+    500,
+    (iter(result.seeds), therm_params),
+    every=int(result.N / 5),
+    save="heat_flow"
+)
+
+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.legend()
+
+exact_flow = [C.flow(result.τ, 0), C.flow(result.τ, 1)]
+
+plt.plot(result.τ, exact_flow[0])
+plt.plot(result.τ, -full_flow[-1][1][0])
+
+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()

python/energy_flow_proper/05_gaussian_two_baths/config.py

@@ -137,7 +137,7 @@ def make_config(
 
     return params
 
-make_config(
+params = make_config(
     max_HO_level=8,
     bcf_terms=4,
     t_max=50,

python/energy_flow_proper/05_gaussian_two_baths/notebook.org

@@ -56,7 +56,7 @@ Let us set up the system with some arbitrary parameters.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/507739453e96adcb21a75433edb50250d1bc49fe.svg]]
+[[file:./.ob-jupyter/3640b2dcf13a8267e9318699695747fb7e0dd296.svg]]
 
 #+begin_src jupyter-python :results none
   params = gf.SystemParams(Ω=1, Λ=1, η=[.5, .5], γ=.1, α_0=[α_0_1, α_0_1])
@@ -78,7 +78,7 @@ coordinate operators in a coherent state.
 #+end_src
 
  #+RESULTS:
- [[file:./.ob-jupyter/75baafbe4ae2c43b4bf91e8b7a3e91e335c64a36.svg]]
+ [[file:./.ob-jupyter/f0b1b1c19cf787ada6ba4019b7e064aeeabb1472.svg]]
 
 * Correlations
 And now the correlation matrix when both oscillators start in the excited state.
@@ -97,7 +97,7 @@ This allows us to calculate the system energy.
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/22efb914cba42853def8243dd5dc5d8b05bb31d9.svg]]
+[[file:./.ob-jupyter/b1ad7d91b497a32ea28f5ce9477a3eb40bb781c7.svg]]
 
 The fast oscillations come from the inter-oscillator coupling.  We
 arrive at some sort of steady state by virtue of the temperature
@@ -109,13 +109,13 @@ difference bytween the baths.
 
 #+begin_src jupyter-python
   with ut.hiro_style():
-      plt.plot(t_points, flow[0] + flow[1])
+      plt.plot(t_points, flow[0])
 
-      #plt.plot(t_points, flow[1])
+      plt.plot(t_points, flow[1])
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/821fb7d3798c55116d1f62bcdcd208233c420fca.svg]]
+[[file:./.ob-jupyter/94827d6dc906bee81ec076020f2479caeb368bce.svg]]
 
 Interestingly these flows don't seem to cancel out.
 
@@ -124,7 +124,7 @@ Interestingly these flows don't seem to cancel out.
 #+end_src
 
 #+RESULTS:
-: -2.467210888582958
+: nan
 
 * Notes
 - precision of nonzero BCF fit is important: zero and nonzero have to be conistent

python/energy_flow_proper/05_gaussian_two_baths/results/heat_flow__heat_flow_ensemble_body_100_2000_8a6864682469e106466d59946c340d6a1816e5ba420caf821b69560ebdf6fadd.json

@@ -0,0 +1 @@
+{"N": 100, "every": 2000, "const_args": ["<not serializable>", "<not serializable>", false], "const_kwargs": {}, "function_name": "_heat_flow_ensemble_body", "first_iterator_value": "<not serializable>"}

python/energy_flow_proper/05_gaussian_two_baths/results/heat_flow__heat_flow_ensemble_body_200_2000_96d9b3eba5de1698ef0b1611cba61509090e8e39e9d3d9258cb292498b9fdb69.json

@@ -0,0 +1 @@
+{"N": 200, "every": 2000, "const_args": ["<not serializable>", "<not serializable>", false], "const_kwargs": {}, "function_name": "_heat_flow_ensemble_body", "first_iterator_value": "<not serializable>"}

python/energy_flow_proper/05_gaussian_two_baths/results/heat_flow__heat_flow_ensemble_body_500_2000_682ce0d1b2d8b028cb7c6a02e6faff13ee5fb01582be1253db6bbe9db4e15312.json

@@ -0,0 +1 @@
+{"N": 500, "every": 2000, "const_args": ["<not serializable>", "<not serializable>", false], "const_kwargs": {}, "function_name": "_heat_flow_ensemble_body", "first_iterator_value": "<not serializable>"}

python/energy_flow_proper/05_gaussian_two_baths/utilities.py

@@ -113,18 +113,30 @@ def plot_complex(x, y, *args, ax=None, label="", **kwargs):
 
 @wrap_plot
 def plot_convergence(
-    x, y, ax=None, label="", transform=lambda y: y, slice=None, linestyle="-"
+    x,
+    y,
+    ax=None,
+    label="",
+    transform=lambda y: y,
+    slice=None,
+    linestyle="-",
+    bath: int = 0,
 ):
     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), label=f"{label}n={n}", alpha=n / y[-1][0], linestyle="--"
+            x,
+            transform(val[bath]),
+            label=f"{label}n={n}",
+            alpha=n / y[-1][0],
+            linestyle="--",
         )
+
     ax.errorbar(
         x,
-        transform(y[-1][1]),
-        yerr=y[-1][2],
+        transform(y[-1][1][bath]),
+        yerr=y[-1][2][bath],
         ecolor="yellow",
         label=f"{label}n={y[-1][0]}",
         color="red",