add data for 05

hops_data/mount_workhorse.sh

@@ -0,0 +1,3 @@
+#! /usr/bin/env bash
+
+sudo mount.nfs 141.30.17.16:/data_local/hiro/hops $(dirname "$0")

hops_data/readme.md

@@ -0,0 +1,2 @@
+This is the directory for hops simulation data.
+If there's nothing in it, it means that you have to mount it first.

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.org

@@ -220,9 +220,10 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
   from hopsflow import gaussflow_two as gf
   from hopsflow import util as util
 
-  α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=8, resolution=0.1)
+  α_0 = util.BCF(params.IntP.t_max, hops_bcf, num_terms=12, resolution=0.1/2)
 #+end_src
 
+
 #+begin_src jupyter-python :results none
   α = util.BCF(
       params.IntP.t_max,
@@ -232,13 +233,13 @@ 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=7,
-      resolution=0.1,
+      num_terms=9,
+      resolution=0.1/2,
   )
 #+end_src
 
 #+begin_src jupyter-python
-  t_test = np.linspace(0, α.t_max, 500)
+  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")
@@ -246,8 +247,8 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 
 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/4869651a72aa97ebc614977ef2209cc6d04cf1fd.svg]]
-[[file:./.ob-jupyter/9a7489ea460c1f7af4a4b8d3a5a6ad6c08879e88.svg]]
+[[file:./.ob-jupyter/eb3f1978b8663a27a45d24124456861cbb22b4c1.svg]]
+[[file:./.ob-jupyter/fd3c2dabd7d087e246a11747045d85fe0947099b.svg]]
 :END:
 
 #+begin_src jupyter-python :results none
@@ -269,7 +270,69 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-: 3135
+: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:
 
 * System Energy
 #+begin_src jupyter-python
@@ -285,8 +348,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-: 100% 3134/3134 [01:06<00:00, 46.88it/s]
-: INFO:root:Writing cache to: results/system_energy_sandwhich_operator_3135_1000_bdddf4a64ad8bc388abd582bceb8d039f8815ee5a66beea77716ab355baa3949.npy
+: INFO:root:Loading cache from: results/system_energy_sandwhich_operator_10000_1000_d0ed047dab9d3ae1be74b1e5b17755f0c1a41bd611212b504a35555740e6bf85.npy
 
 
 #+begin_src jupyter-python
@@ -295,7 +357,7 @@ This will be tangled into the [[file:config.py][config file]] that can be used w
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/379210c06234a6a5a5d340133be180e808247f50.svg]]
+[[file:./.ob-jupyter/5bd904b9184abb5bc91b85d68a42bebe5ba1f5ae.svg]]
 
 Now with Gaussflow.
 #+begin_src jupyter-python
@@ -304,35 +366,34 @@ Now with Gaussflow.
 #+end_src
 
 #+RESULTS:
-: 32a8d025-ef3b-4ce6-80d2-1489d98eb193
 
 
 #+begin_src jupyter-python
   with ut.hiro_style():
       ax.plot(result.τ, energy_gf, color="green", linestyle="dashdot")
-      #plt.xlim(30,40)
+ #plt.xlim(30,40)
       #plt.ylim(-.005,.005)
   fig
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/5222fd8000ff9df3cc5ea289ac13ae17e062c8d7.svg]]
+[[file:./.ob-jupyter/9975e4a825c9842e39abf68b76523672581c7c2d.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)
+      # plt.ylim(0,.01)D
       plt.legend()
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/9c1bdc31919fb6a74990e6570510914cfa14e8d7.svg]]
+[[file:./.ob-jupyter/2f82b9edaf51d90e93a9cbc50fda894b9a8d96fb.svg]]
 
 * Thermal State
 #+begin_src jupyter-python
-  v, w = np.linalg.eig(system.H_sys)
+  v, w = np.linalg.eig(params.SysP.H_sys)
   v = np.sort(v.real)
   v
 #+end_src
@@ -364,37 +425,17 @@ Now with Gaussflow.
 #+end_src
 
 #+begin_src jupyter-python
-  thermal_e(v, 1 / .3)
+  thermal_e(v, 1.0 / (params.SysP.__non_key__["T"][0]))
 #+end_src
 
 #+RESULTS:
-: 1.1140772133416454
+: 1.0591063695586198
 
 #+begin_src jupyter-python
-  ρ_therm = sc.linalg.expm(-system.H_sys/(system.__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)
 #+end_src
 
 #+RESULTS:
-: 0.16441343359804375
-
-* Propagator
-#+begin_src jupyter-python
-  _, q_1_traj, σ_q_1_traj = util.operator_expectation_ensemble(
-      iter(result.ψ),
-      p_1.full(),
-      result.N,
-      params.HiP.nonlinear,
-      save="q_1"
-  )
-  plt.errorbar(result.τ, q_1_traj.real, yerr=σ_q_1_traj.real, ecolor="yellow")
-#+end_src
-
-#+RESULTS:
-:RESULTS:
-: 100% 1021/1021 [00:02<00:00, 428.10it/s]
-: INFO:root:Writing cache to: results/q_1_sandwhich_operator_1022_None_3422a7089748f5a2a9b27570e4c85e3ea7e48865e48733fd74e8438c6fa0c512.npy
-: <ErrorbarContainer object of 3 artists>
-[[file:./.ob-jupyter/4bbc638bea9e95d5b2c8ef3e17a14a718e177eef.svg]]
-:END:
+: 0.07726133069400465

python/energy_flow_proper/05_gaussian_two_baths/comparison_with_hops.py

@@ -0,0 +1,126 @@
+params = make_config(
+    max_HO_level=6,
+    bcf_terms=4,
+    t_max=40,
+    k_fac=1.6,
+    sp_tol=1e-5,
+)
+
+params = make_config(
+    max_HO_level=8,
+    bcf_terms=4,
+    t_max=30,
+    k_fac=1.3,
+    sp_tol=1e-5,
+    bcf_scale=[.3, .3],
+    T=.2
+)
+
+params = make_config(
+    max_HO_level=5,
+    bcf_terms=4,
+    t_max=50,
+    k_fac=1.3,
+    sp_tol=1e-5,
+    bcf_scale=[.3, .3],
+    T=.2
+)
+
+params = make_config(
+    max_HO_level=8,
+    bcf_terms=4,
+    t_max=50,
+    k_fac=1.4,
+    sp_tol=1e-5,
+    bcf_scale=[.3, .3],
+    T=.2
+)
+
+import numpy as np
+import scipy as sc
+import matplotlib.pyplot as plt
+from hopsflow import util, hopsflow
+import utilities as ut
+from hops.core.hierarchy_data import HIMetaData
+import logging
+logger = logging.getLogger()
+logger.setLevel(logging.INFO)
+
+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)
+
+α = util.BCF(
+    params.IntP.t_max,
+    hops.util.bcf.OhmicBCF_nonZeroTemp(
+        s,
+        1,
+        wc,
+        beta=1.0 / (params.SysP.__non_key__["T"][0]),
+    ),
+    num_terms=9,
+    resolution=0.1/2,
+)
+
+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)
+
+class result:
+    hd = HIMetaData("data", ".").get_HIData(params, read_only=True)
+    N = hd.samples
+    seeds = hd.rng_seed
+    τ = hd.get_time()
+    ψ_1 = hd.aux_states
+    ψ = hd.stoc_traj
+    ρ = hd.get_rho_t()
+
+result.N
+
+e_sys = util.operator_expectation_ensemble(
+    iter(result.ψ),
+    params.SysP.H_sys,
+    result.N,
+    params.HiP.nonlinear,
+    every=1000,
+    real=True,
+    save="system_energy"
+)
+
+with ut.hiro_style():
+    fig, ax = ut.plot_convergence(result.τ, e_sys, transform=np.real, linestyle="dashdot")
+
+C = gf.CorrelationMatrix(gf_params, gf.initial_correlation_pure_osci(0, 0), αs=[α, None])
+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 ut.hiro_style():
+    ut.plot_diff_vs_sigma(result.τ, e_sys, energy_gf)
+    # plt.xlim(35,40)
+    # plt.ylim(0,.01)D
+    plt.legend()
+
+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))
+
+def trace_norm(ρ: np.ndarray):
+      return np.trace(sc.linalg.sqrtm(ρ @ ρ.conj().T)).real
+
+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)