master-thesis/python/graveyard/richard_hops/energy_flow_gauss.org

#+PROPERTY: header-args :session rich_hops_eflow :kernel python :pandoc t :async yes

* Setup
** Jupyter
#+begin_src jupyter-python
  %load_ext autoreload
  %autoreload 2
  %load_ext jupyter_spaces
#+end_src

#+RESULTS:
: The autoreload extension is already loaded. To reload it, use:
:   %reload_ext autoreload
: The jupyter_spaces extension is already loaded. To reload it, use:
:   %reload_ext jupyter_spaces

** Matplotlib
#+begin_src jupyter-python
  import matplotlib
  import matplotlib.pyplot as plt

  #matplotlib.use("TkCairo", force=True)
  %gui tk
  %matplotlib inline
  plt.style.use('ggplot')
#+end_src

#+RESULTS:

** Richard (old) HOPS
#+begin_src jupyter-python
  import hierarchyLib
  import hierarchyData
  import numpy as np

  from stocproc.stocproc import StocProc_FFT, StocProc_KLE
  import bcf
  from dataclasses import dataclass
  import scipy
  import scipy.misc
  import scipy.signal
#+end_src

#+RESULTS:

** Auxiliary Definitions
#+begin_src jupyter-python
  σ1 = np.matrix([[0,1],[1,0]])
  σ2 = np.matrix([[0,-1j],[1j,0]])
  σ3 = np.matrix([[1,0],[0,-1]])
#+end_src

#+RESULTS:

* Model Setup
Basic parameters.
#+begin_src jupyter-python
  γ = 5 # coupling ratio
  ω_c = 0  # center of spect. dens
  δ = .1  # breadth BCF
  t_max = 10
  t_steps = 500
  k_max = 6
  seed = 100
  H_s = σ3 + np.eye(2)
  L =  1 / 2 * (σ1 - 1j * σ2) * γ
  ψ_0 = np.array([1, 0])
  W = ω_c * 1j + δ  # exponent BCF
  N = 100
#+end_src

#+RESULTS:

** BCF
#+begin_src jupyter-python
  @dataclass
  class CauchyBCF:
      δ: float
      ω_c: float

      def I(self, ω):
          return np.sqrt(self.δ) / (self.δ + (ω - self.ω_c) ** 2 / self.δ)

      def __call__(self, τ):
          return np.sqrt(self.δ) * np.exp(-1j * self.ω_c * τ - np.abs(τ) * self.δ)

      def __bfkey__(self):
          return self.δ, self.ω_c

  @dataclass
  class GaussBCF:
      σ: float
      ω_c: float

      def I(self, ω):
          return (
              np.exp(-(((ω - self.ω_c) / self.σ) ** 2) / 2)
              ,* 1
              / (np.sqrt(2 * np.pi) * self.σ) * np.pi
          )

      def __call__(self, τ):
          return (
              np.exp(-(((τ - self.ω_c) / self.σ) ** 2) / 2)
              ,* 1
              / (np.sqrt(2 * np.pi) * self.σ) * np.pi
          )
      #np.exp(1j * self.ω_c * τ - self.σ**2 * τ**2/2)

      def __bfkey__(self):
          return self.σ, self.ω_c


  α = GaussBCF(δ, ω_c)
#+end_src

#+RESULTS:

*** Plot
#+begin_src jupyter-python
  %%space plot
  t = np.linspace(0, t_max, 1000)
  ω = np.linspace(ω_c - 10, ω_c + 10, 1000)
  fig, axs = plt.subplots(2)
  axs[0].plot(t, np.real(α(t)))
  axs[0].plot(t, np.imag(α(t)))
  axs[1].plot(ω, α.I(ω))
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7f6798323760> |
| <matplotlib.lines.Line2D | at | 0x7f6799a6ea60> |
| <matplotlib.lines.Line2D | at | 0x7f6799a6eca0> |
[[file:./.ob-jupyter/e79e0652b99bba8df94a7f60c64af072947ade03.png]]
:END:

** Hops setup
#+begin_src jupyter-python
  HierachyParam = hierarchyData.HiP(
      k_max=k_max,
      # g_scale=None,
      # sample_method='random',
      seed=seed,
      nonlinear=False,
      # normalized=False,
      # terminator=False,
      result_type=hierarchyData.RESULT_TYPE_ALL,
      # accum_only=None,
      # rand_skip=None
  )
#+end_src

#+RESULTS:

Integration.
#+begin_src jupyter-python
  IntegrationParam = hierarchyData.IntP(
      t_max=t_max,
      t_steps=t_steps,
      # integrator_name='zvode',
      # atol=1e-8,
      # rtol=1e-8,
      # order=5,
      # nsteps=5000,
      # method='bdf',
      # t_steps_skip=1
  )
#+end_src

#+RESULTS:

And now the system.
#+begin_src jupyter-python
  SystemParam = hierarchyData.SysP(
      H_sys=H_s,
      L=L,
      psi0=ψ_0,  # excited qubit
      g=np.array([np.sqrt(δ)]),
      w=np.array([W]),
      H_dynamic=[],
      bcf_scale=1,  # some coupling strength (scaling of the fit parameters 'g_i')
      gw_hash=None,  # this is used to load g,w from some database
      len_gw=1,
  )
#+end_src

#+RESULTS:

The quantum noise.
#+begin_src jupyter-python
  Eta = StocProc_KLE(
      α,
      t_max,
      seed=seed,
      tol=1e-3
  )
#+end_src

#+RESULTS:
:RESULTS:
#+begin_example
  stocproc.method_kle - INFO - check 33 grid points
  stocproc.method_kle - INFO - check 65 grid points
  alpha_k is real
  alpha_k is real
  stocproc.method_kle - INFO - check 129 grid points
  alpha_k is real
  stocproc.method_kle - INFO - check 257 grid points
  alpha_k is real
  stocproc.method_kle - INFO - check 513 grid points
  alpha_k is real
#+end_example
# [goto error]
#+begin_example

  KeyboardInterruptTraceback (most recent call last)
  <ipython-input-777-b4dd44af7ffd> in <module>
  ----> 1 Eta = StocProc_KLE(
        2     α,
        3     t_max,
        4     seed=seed,
        5     tol=1e-3

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/stocproc/stocproc.py in __init__(self, alpha, t_max, tol, ng_fac, meth, diff_method, dm_random_samples, seed, align_eig_vec, scale)
      286         key = alpha, t_max, tol
      287
  --> 288         sqrt_lambda_ui_fine, t = method_kle.auto_ng(
      289             corr=alpha,
      290             t_max=t_max,

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/stocproc/method_kle.py in auto_ng(corr, t_max, ngfac, meth, tol, diff_method, dm_random_samples, ret_eigvals, relative_difference)
      543                     ui_super_fine.reshape(1, -1)
      544                 )
  --> 545             md = np.max(np.abs(diff) / abs_alpha_res)
      546             time_calc_diff += time.time() - t0
      547

  KeyboardInterrupt:
#+end_example
:END:
* Actual Hops
Generate the key for binary caching.
#+begin_src jupyter-python
  hi_key = hierarchyData.HIMetaKey_type(
      HiP=HierachyParam,
      IntP=IntegrationParam,
      SysP=SystemParam,
      Eta=Eta,
      EtaTherm=None,
  )
#+end_src

#+RESULTS:

Initialize Hierarchy.
#+begin_src jupyter-python
  myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=N, desc="run a test case")
#+end_src

#+RESULTS:
: init Hi class, use 14 equation
: /home/hiro/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py:1058: UserWarning: sum_k_max is not implemented! DO SO BEFORE NEXT USAGE (use simplex).HierarchyParametersType does not yet know about sum_k_max
:   warnings.warn(

Run the integration.
#+begin_src jupyter-python :results none
  myHierarchy.integrate_simple(data_name="energy_flow.data")
#+end_src


Get the samples.
#+begin_src jupyter-python
  # to access the data the 'hi_key' is used to find the data in the hdf5 file
  with hierarchyData.HIMetaData(hid_name="energy_flow.data", hid_path=".") as metaData:
      with metaData.get_HIData(hi_key, read_only=True) as data:
          smp = data.get_samples()
          print("{} samples found in database".format(smp))
          τ = data.get_time()
          rho_τ = data.get_rho_t()
          s_proc = np.array(data.stoc_proc)
          states = np.array(data.aux_states).copy()
          ψ_1 = np.array(data.aux_states)[:, :, 0:2]
          ψ_0 = np.array(data.stoc_traj)
          y = np.array(data.y)
#+end_src

#+RESULTS:
:RESULTS:
# [goto error]
#+begin_example

  KeyErrorTraceback (most recent call last)
  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData_fname(self, key)
      779         try:
  --> 780             hdf5_name = self.db[hashed_key][0]
      781         except KeyError:

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in __getitem__(self, key)
      243         if item is None:
  --> 244             raise KeyError(key)
      245         return self.decode(item[0])

  KeyError: '47177f42579f772ba59f8489393910e4b2e9b1f2567e082f6b944d00382a9df7793b33c105854cabecd012bc9f3fb59d'

  During handling of the above exception, another exception occurred:

  PicklingErrorTraceback (most recent call last)
  <ipython-input-781-12d991d6efe0> in <module>
        1 # to access the data the 'hi_key' is used to find the data in the hdf5 file
        2 with hierarchyData.HIMetaData(hid_name="energy_flow.data", hid_path=".") as metaData:
  ----> 3     with metaData.get_HIData(hi_key, read_only=True) as data:
        4         smp = data.get_samples()
        5         print("{} samples found in database".format(smp))

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData(self, key, read_only)
      787
      788     def get_HIData(self, key, read_only=False):
  --> 789         hdf5_name = self.get_HIData_fname(key)
      790
      791         if key.HiP.result_type == RESULT_TYPE_ZEROTH_ORDER_ONLY:

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData_fname(self, key)
      781         except KeyError:
      782             hdf5_name = self._new_rand_file_name(pre=self.name + "_", end=".h5")
  --> 783             self.db[hashed_key] = (hdf5_name, key)
      784             self.db.commit()
      785

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in __setitem__(self, key, value)
      250
      251         ADD_ITEM = 'REPLACE INTO "%s" (key, value) VALUES (?,?)' % self.tablename
  --> 252         self.conn.execute(ADD_ITEM, (key, self.encode(value)))
      253         if self.autocommit:
      254             self.commit()

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in encode(obj)
       95 def encode(obj):
       96     """Serialize an object using pickle to a binary format accepted by SQLite."""
  ---> 97     return sqlite3.Binary(dumps(obj, protocol=PICKLE_PROTOCOL))
       98
       99

  PicklingError: Can't pickle <class '__main__.GaussBCF'>: it's not the same object as __main__.GaussBCF
#+end_example
:END:

Calculate energy.
#+begin_src jupyter-python
  energy = np.array([np.trace(ρ * H_s).real/np.trace(ρ).real for ρ in rho_τ])
  plt.plot(τ, energy)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7f67981b3e80> |
[[file:./.ob-jupyter/892241bc3127d1dc4581b7f31e8b392425c9716c.png]]
:END:

#+begin_src jupyter-python
  %%space plot
  plt.plot(τ, np.trace(rho_τ.T).real)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7f6799c2caf0> |
[[file:./.ob-jupyter/d79cca7773c99c7279d79f331d77f301acbdb71d.png]]
:END:

* Energy Flow
:PROPERTIES:
:ID:       eefb1594-e399-4d24-9dd7-a57addd42e65
:END:
#+begin_src jupyter-python
  ψ_1.shape
#+end_src

#+RESULTS:
| 160 | 500 | 2 |

Let's look at the norm.
#+begin_src jupyter-python
  plt.plot(τ, (ψ_1[0].conj() * ψ_1[0]).sum(axis=1).real)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7f679a203700> |
[[file:./.ob-jupyter/5783eff51f0815225fe39255d41152e4b447963f.png]]
:END:

And try to calculate the energy flow.
#+begin_src jupyter-python
  def flow_for_traj(ψ_0, ψ_1):
      a = np.array((L @ ψ_0.T).T)

      return np.array(2 * (1j * -W * np.sum(a.conj() * ψ_1, axis=1)).real).flatten()


  def flow_for_traj_alt(ψ_0, y):
      Eta.new_process(y)
      eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-8)
      a = np.array((L @ ψ_0.T).T)

      return -(2j * eta_dot.conj() *
          np.array((np.sum(ψ_0.conj() * a, axis=1))).flatten()
      ).real
#+end_src

#+RESULTS:

Now we calculate the average over all trajectories.
#+begin_src jupyter-python
  j = np.zeros_like(τ)
  for i in range(0, N):
      j += flow_for_traj(ψ_0[i], ψ_1[i])
  j /= N
#+end_src

#+RESULTS:
:RESULTS:
# [goto error]
#+begin_example

  ValueErrorTraceback (most recent call last)
  <ipython-input-787-b9081128ed40> in <module>
        1 j = np.zeros_like(τ)
        2 for i in range(0, N):
  ----> 3     j += flow_for_traj(ψ_0[i], ψ_1[i])
        4 j /= N

  <ipython-input-786-c57f86a6b31b> in flow_for_traj(ψ_0, ψ_1)
        1 def flow_for_traj(ψ_0, ψ_1):
  ----> 2     a = np.array((L @ ψ_0.T).T)
        3
        4     return np.array(2 * (1j * -W * np.sum(a.conj() * ψ_1, axis=1)).real).flatten()
        5

  ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)
#+end_example
:END:

And do the same with the alternative implementation.
#+begin_src jupyter-python
  ja = np.zeros_like(τ)
  for i in range(0, N):
      ja += flow_for_traj_alt(ψ_0[i], y[i])
  ja /= N
#+end_src

#+RESULTS:
:RESULTS:
# [goto error]
#+begin_example

  ValueErrorTraceback (most recent call last)
  <ipython-input-788-5dc6ccf09941> in <module>
        1 ja = np.zeros_like(τ)
        2 for i in range(0, N):
  ----> 3     ja += flow_for_traj_alt(ψ_0[i], y[i])
        4 ja /= N

  <ipython-input-786-c57f86a6b31b> in flow_for_traj_alt(ψ_0, y)
        8     Eta.new_process(y)
        9     eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-8)
  ---> 10     a = np.array((L @ ψ_0.T).T)
       11
       12     return -(2j * eta_dot.conj() *

  ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)
#+end_example
:END:

And plot it :)
#+begin_src jupyter-python
  %matplotlib inline
  plt.plot(τ, j)
  plt.plot(τ, ja)
  plt.show()
#+end_src

#+RESULTS:
[[file:./.ob-jupyter/b159f74a1bffc1460ac7f285f76ab1b81f31bd07.png]]

\Let's calculate the integrated energy.
#+begin_src jupyter-python
  E_t = np.array([0] + [scipy.integrate.simpson(j[0:n], τ[0:n]) for n in range(1, len(τ))])
  E_t[-1]
#+end_src

#+RESULTS:
: 0.0

#+begin_src jupyter-python
  E_t = np.array([0] + [scipy.integrate.simpson(ja[0:n], τ[0:n]) for n in range(1, len(τ))])
  E_t[-1]
#+end_src

#+RESULTS:
: 0.0

With this we can retrieve the energy of the interaction Hamiltonian.
#+begin_src jupyter-python
  E_I = 2 - energy - E_t
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  %%space plot
  plt.rcParams['figure.figsize'] = [10, 8]
  #plt.plot(τ, j, label="$J$", linestyle='--')
  plt.plot(τ, E_t, label=r"$\langle H_{\mathrm{B}}\rangle$")
  plt.plot(τ, E_I, label=r"$\langle H_{\mathrm{I}}\rangle$")
  plt.plot(τ, energy, label=r"$\langle H_{\mathrm{S}}\rangle$")

  plt.xlabel("τ")
  plt.legend()
  plt.show()
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7f6798fb5eb0> |
| <matplotlib.lines.Line2D | at | 0x7f6798f8cb50> |
| <matplotlib.lines.Line2D | at | 0x7f6798f75d60> |
: Text(0.5, 0, 'τ')
: <matplotlib.legend.Legend at 0x7f6798f8ceb0>
[[file:./.ob-jupyter/ca75f76563aca062310d5779c0c5df539b728d3f.png]]
:END:
#+RESULTS:

* Derivatives
#+begin_src jupyter-python
  Eta.new_process(y[0])
  #plt.plot(τ, Eta(τ).real)
  eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-3)
  plt.plot(τ, eta_dot)
#+end_src

#+RESULTS:
:RESULTS:
: /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
:   return array(a, dtype, copy=False, order=order)
| <matplotlib.lines.Line2D | at | 0x7f678d445f40> |
[[file:./.ob-jupyter/aaa6bdcc61452758e092eccc2d62ef8cdfe5f258.png]]
:END: