master-thesis/python/energy_flow_proper/03_gaussian/naive_integration.org

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

#+begin_src jupyter-python :results none
  import scipy
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import integrate
  import sys
  sys.path.append("/home/hiro/python/energy_flow_proper/hops")
  sys.path.append("/home/hiro/python/energy_flow_proper")
  import utilities
  import quadpy
#+end_src


* System Parameters
#+begin_src jupyter-python :results none
  Ω = 1.5
  A = np.array([[0, Ω], [-Ω, 0]])
  η = 2
  ω_c = 1
  t_max = 16


  def α(t):
      return η / np.pi * (ω_c / (1 + 1j * ω_c * t)) ** 2


  def α_0_dot(t):
      return 2 * η / (1j * np.pi) * (ω_c / (1 + 1j * ω_c * t)) ** 3


  def K(τ):
      return np.array([[0, 0], [α(τ).imag, 0]])


  def L(t, s):
      return np.exp(-A * t) @ K(t - s) @ np.exp(A * s)
#+end_src

* Definition of the RHS
#+begin_src jupyter-python :results none
  def H_dot(t, _):
      return -integrate.trapezoid(
          [(L(t, ts[s]) @ np.array(H_to_date[s])).flatten() for s in range(len(ts))],
          ts,
          axis=0,
      )


  def K_dot(t, curr):
      return (A @ curr.reshape(2, 2)).flatten() - integrate.simpson(
          [(K(t - ts[s]) @ np.array(H_to_date[s])).flatten() for s in range(len(ts))]
          + [(K(0) @ curr.reshape(2, 2)).flatten()],
          ts + [t],
          axis=0,
      )
#+end_src

* Scipy Hack
Hacking together an integro-differential solver.
#+begin_src jupyter-python :results none
  ts = [0]
  H_to_date = [np.eye(2)]
  r = integrate.ode(K_dot).set_integrator('vode')
  r.set_initial_value(np.eye(2, dtype="float").flatten(), 0)
  t1 = t_max
  dt = .001
  while r.successful() and r.t < t1:
      sol = r.integrate(r.t + dt)
      H_to_date.append(sol.reshape(2,2))
      ts.append(r.t + dt)
#+end_src

Dividing out the unitary dynamics.
#+begin_src jupyter-python
  #plt.plot(ts,  (np.array([scipy.linalg.expm(-A * np.array(t)) for t in ts]) @ np.array(H_to_date))[:,1,0])
  proper = np.array(H_to_date)
#+end_src

#+RESULTS:

\(\langle q\rangle\)
#+begin_src jupyter-python
  plt.plot(ts, proper @ np.array([1, 0]))
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd11c35820> | <matplotlib.lines.Line2D | at | 0x7efd11c35850> |
[[file:./.ob-jupyter/699a32087aa1ac7105aef99a1b81469bbeb8d28a.svg]]
:END:

This looks like a slightly dampened HO.

* IDESolver
With a proper solver it works like a charm.
See [[https://idesolver.readthedocs.io/en/latest/manual.html][the docs]].
#+begin_src jupyter-python
  from idesolver import IDESolver
  steps = 500
  t = np.linspace(0, t_max, steps)
  solver = IDESolver(
      x = t,
      y_0 = np.array([1, 0, 0, 1]),
      c = lambda x, y: (A @ y.reshape(2, 2)).flatten(),
      d = lambda x: 1,
      k = lambda x, s: -K(x-s)[1, 0],
      f = lambda y: np.array([0, 0, y[0], y[1]]),
      lower_bound = lambda x: 0,
      upper_bound = lambda x: x,
  )
  sol = solver.solve()
#+end_src

#+RESULTS:
#+begin_example
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 0
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 1
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 2
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 3
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 4
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 5
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 6
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 7
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 8
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 9
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 10
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 11
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 12
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 13
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 14
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 15
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 16
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 17
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 18
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 19
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 20
    warnings.warn(
  /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 24
    warnings.warn(
#+end_example


Reshape \(G\) into a time series of matrices.
#+begin_src jupyter-python
  G = np.einsum("ijk->kij", sol.reshape((2,2,steps)))
#+end_src

#+RESULTS:

And plot the time development of the mean values of \(p\) and \(p\).
#+begin_src jupyter-python
  plt.plot(t, G @ np.array([1, 0]))
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd11708f70> | <matplotlib.lines.Line2D | at | 0x7efd11708490> |
[[file:./.ob-jupyter/fad23f518c037f04e6108539af9a0652901fe1ff.svg]]
:END:
Looks OK.


#+begin_src jupyter-python
  plt.plot(
      t,
      (np.array([scipy.linalg.expm(-A * np.array(τ)) for τ in t]) @ G).reshape(steps, 4),
  )
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd11773f10> | <matplotlib.lines.Line2D | at | 0x7efd11773f40> | <matplotlib.lines.Line2D | at | 0x7efd11799f70> | <matplotlib.lines.Line2D | at | 0x7efd11799e50> |
[[file:./.ob-jupyter/e731bf5509fa205c540e9062c3a15c878b8a4a21.svg]]
:END:
We see that the derivative at the begining is really 0 if we divide
out the unitary dynamics.


The scipy integration arrives at pretty much the same picture :).
#+begin_src jupyter-python
  plt.plot(
      ts,
      (np.array([scipy.linalg.expm(-A * np.array(t)) for t in ts]) @ np.array(H_to_date)).reshape(len(ts), 4),
      label="scipy"
  )
  plt.plot(
      t,
      (np.array([scipy.linalg.expm(-A * np.array(τ)) for τ in t]) @ G).reshape(steps, 4),
      linestyle="--",
      label="idesolve"
  )
  plt.legend()
  plt.xlabel('t')
  plt.ylabel('$e^{-At}G_{ij}$')
#+end_src

#+RESULTS:
:RESULTS:
: Text(0, 0.5, '$e^{-At}G_{ij}$')
[[file:./.ob-jupyter/18f4e12520e5e460dfc0d1a5b4ba87d7d087be39.svg]]
:END:

The initial slip is removed in idesolve due to the iteration. You have
to look at the off diagonal curves to see that.

#+begin_src jupyter-python
  plt.plot(t, np.linalg.det(G))
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0f57e190> |
[[file:./.ob-jupyter/0d736728dd6cac1d188363bc857527c2c70fcc4e.svg]]
:END:

We see that the "volume" bleeds out of the system.

* Interpolation
#+begin_src jupyter-python
  import fcSpline
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  class FCSWrap:
      def __init__(self, t_min, t_max, G):
          self._G = G
          self._t_min = t_min
          self._t_max = t_max

          self._splines = np.array(
              list(map(lambda g: fcSpline.FCS(self._t_min, self._t_max, g), G))
          )

      def __call__(self, t):
          t = np.asarray(t).astype('float64')

          res = np.array(list(map(lambda g: g(t), self._splines)))
          if t.size == 1:
              return res.reshape((2,2))

          return res.reshape((2,2,t.size)).swapaxes(0,2).swapaxes(1,2)

      def __getitem__(self, key):
          return self._splines.reshape(2,2)[key]

#+end_src

#+RESULTS:

#+begin_src jupyter-python
  G_inter = FCSWrap(t.min(), t.max(), sol)
#+end_src

#+RESULTS:

* Calculate the Convolutions
We use quadpy for complex integration.
#+begin_src jupyter-python
    import quadpy
#+end_src

#+RESULTS:
We will need \(G_{12}\) often.
#+begin_src jupyter-python
  G12 = G_inter[0, 1]
#+end_src

#+RESULTS:

Calculate \(G_{12}\ast \dot{\alpha_0}\).
#+begin_src jupyter-python :results none
  G12_star_α_0_dot = np.array(
      [
          quadpy.quad(lambda τ: G12(τ) * α_0_dot(t_now - τ), 0, t_now)[0]
          for t_now in t
      ]
  )
  G12_star_α_0_dot = fcSpline.FCS(t.min(), t.max(), G12_star_α_0_dot)
#+end_src

Calculate \(\tau_r\alpha\ast G_{12}\).
#+begin_src jupyter-python :results none
  τ_r_alpha_star_G12 = np.array(
      [
          [
              quadpy.quad(lambda τ: G_inter[0, 1](t_now - τ) * α_0_dot(τ - r), 0, t_now)[
                  0
              ]
              for r in t
          ]
          for t_now in t
      ]
  )
  τ_r_alpha_star_G12_real = scipy.interpolate.interp2d(t, t, τ_r_alpha_star_G12.real)
  τ_r_alpha_star_G12_imag = scipy.interpolate.interp2d(t, t, τ_r_alpha_star_G12.imag)

  def τ_r_alpha_star_G12_inter(t, r):
      return τ_r_alpha_star_G12_real(t, r) + 1j * τ_r_alpha_star_G12_imag(t, r)
#+end_src


#+begin_src jupyter-python
  plt.plot(t, G12_star_α_0_dot(t).imag)
  plt.plot(t, G12_star_α_0_dot(t).real)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0de08520> |
[[file:./.ob-jupyter/71d347c687e062a6feb23ab0f42f6b4f328155b6.svg]]
:END:

#+begin_src jupyter-python
  plt.imshow(τ_r_alpha_star_G12.real)
#+end_src

#+RESULTS:
:RESULTS:
: <matplotlib.image.AxesImage at 0x7efd0d9abbe0>
[[file:./.ob-jupyter/0dcad0dda95f3d0c88222bfc06e3ab453b898aa5.svg]]
:END:

#+begin_src jupyter-python
    plt.imshow(τ_r_alpha_star_G12.imag)
#+end_src

#+RESULTS:
:RESULTS:
: <matplotlib.image.AxesImage at 0x7efd0dd1f820>
[[file:./.ob-jupyter/1c3acd3be3596cb399b651d18fe9b3777dcdf185.svg]]
:END:

And the final integral.
#+begin_src jupyter-python
  flow_half = np.array(
      [
          quadpy.quad(
              lambda r: G12_star_α_0_dot(t_now - r) * τ_r_alpha_star_G12_inter(t_now, r)[0],
              0,
              t_now,
          )[0]
          for t_now in t
      ]
  )
  flow_half = -1/2 * flow_half.imag
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  plt.plot(t, flow_half)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0deab580> |
[[file:./.ob-jupyter/3fd6e2e90774a32ecd7037d93121032e347c999f.svg]]
:END:

And now the \(G\) part of the flow.
#+begin_src jupyter-python
  n = 1
  flow_G_half = np.array(
      [
          quadpy.quad(
              lambda s: (
                  (
                      G_inter[0, 0](t_now) * G_inter[0, 1](s)
                      - G_inter[0, 0](s) * G_inter[0, 1](t_now)
                  )
                  ,* α_0_dot(t_now - s).real
                  + (2 * n + 1)
                  ,* (
                      G_inter[0, 0](t_now) * G_inter[0, 0](s)
                      + G_inter[0, 1](t_now) * G_inter[0, 1](s)
                  )
                  ,* α_0_dot(t_now - s).imag
              ),
              0,
              t_now,
          )[0]
          for t_now in t
      ]
  )
  flow_G_half = -1 / 2 * flow_G_half
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  plt.plot(t, flow_G_half)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0de3fe80> |
[[file:./.ob-jupyter/df69413ad6fb4f37ec6ad22639acc3cebdf7ee0b.svg]]
:END:

#+begin_src jupyter-python
  flow = flow_G_half + flow_half
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  plt.plot(t, flow)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0e01be20> |
[[file:./.ob-jupyter/88b32148be473f0d22726d145bbacf7a6acdf794.svg]]
:END:

#+end_src
#+begin_src jupyter-python
  scipy.integrate.trapz(flow, t)
#+end_src

#+RESULTS:
: 1.5215449446799811

There remains the zero point energy.

#+begin_src jupyter-python
  def integrand(r, s, l, t_now):
      return -1 / 2 * (G12(t_now - l) * G12(s - r) * α(l - r) * α_0_dot(t_now - s)).imag


  flow_half_alt = np.array(
      [
          scipy.integrate.nquad(
              integrand, [lambda s, _, __: [0, s], [0, t_now], [0, t_now]], args=[t_now]
          )
          for t_now in np.linspace(0, t.max(), 100)
      ]
  )
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  flow_half_alt_int = fcSpline.FCS(0, t.max(), flow_half_alt[:, 0])
#+end_src

#+RESULTS:


#+begin_src jupyter-python
  plt.plot(t, flow_half_alt_int(t))
  plt.plot(t, flow_half)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0e3b1eb0> |
[[file:./.ob-jupyter/1d9aaf9144c9a0f5e43b0c97e9b8d1966188f3a8.svg]]
:END:

#+begin_src jupyter-python
  flow_alt = flow_G_half + flow_half_alt_int(t)
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  plt.plot(t, flow_alt)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7efd0bea4dc0> |
[[file:./.ob-jupyter/8a927282d95f0cc9b714cba88fa99ed00ecf5f4f.svg]]
:END:

#+begin_src jupyter-python
  scipy.integrate.trapz(flow_alt, t)
#+end_src

#+RESULTS:
: 1.5885826683889948

* Analytical Solution
** Exponential Fit
First we need an exponential fit for our BCF.
#+begin_src jupyter-python
  W, G_raw = utilities.fit_α(α, 3, 80, 10_000)
  τ = np.linspace(0, t_max, 1000)
#+end_src

#+RESULTS:

Looks quite good.
#+begin_src jupyter-python
  fig, ax = utilities.plot_complex(τ, α(τ), label="exact")
  utilities.plot_complex(
      τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--"
  )
#+end_src

#+RESULTS:
:RESULTS:
| hline | <AxesSubplot:> |
[[file:./.ob-jupyter/afa5b41a1ab07e555bddc54be560b06c11bf3906.svg]]
:END:

** Calculate the Magic Numbers
We begin with the $\varphi_n$ and $G_n$ from the original ~G~.
#+begin_src jupyter-python
  φ = -np.angle(G_raw)
  φ
#+end_src

#+RESULTS:
: array([-1.18710245,  0.64368323,  2.11154332])

#+begin_src jupyter-python
  G = np.abs(G_raw)
  G
#+end_src

#+RESULTS:
: array([0.51238635, 0.62180167, 0.10107935])

Now we calculate the real and imag parts of the ~W~ parameters and
call them $\gamma_n$ and $\delta_n$.

#+begin_src jupyter-python
  γ, δ = W.real, W.imag
#+end_src

#+RESULTS:

Now the \(s_n, c_n\).
#+begin_src jupyter-python
  s, c = np.sin(φ), np.cos(φ)
#+end_src

#+RESULTS:

Now we calculate the roots of $f_n(z)=-G_n ((z+\gamma_n) s_n + \delta_n
c_n)$.
Normally we should be vary of one of the \(\deltas\) being zero, but
this is not the case.
#+begin_src jupyter-python
  roots_f = -(γ + δ * c/s)
  roots_f
#+end_src

#+RESULTS:
: array([-1.19725058, -2.0384181 , -0.23027627])

Now the \(z_k\) the roots of \(\delta_k^2 + (\gamma_k + z)^2\). *We
don't include the conjugates.*
#+begin_src jupyter-python
  roots_z = -W
  roots_z
#+end_src

#+RESULTS:
: array([-2.29230292-2.71252483j, -1.07996581-0.719112j  ,
:        -0.28445344-0.09022829j])

** Construct the Polynomials
#+begin_src jupyter-python
  from numpy.polynomial import Polynomial
#+end_src

#+RESULTS:

We later need \(f_0(z) = \prod_k (z-z_k) (z-z_k^{\ast})\).
#+begin_src jupyter-python
  f_0 = utilities.poly_real(Polynomial.fromroots(np.concatenate((roots_z, roots_z.conj()))))
  f_0
#+end_src

#+RESULTS:

Another polynomial is simply \(p_1 = (z^2 + \Omega^2)\prod_k (z-z_k)
(z-z_k^{\ast})\) and we can construct it from its roots.
#+begin_src jupyter-python
  p_1 = Polynomial([Ω**2, 0, 1]) * f_0
  p_1
#+end_src

#+RESULTS:

The next ones are given through \(q_n=\Omega f_n(z) \prod_{k\neq
n}(z-z_k) (z-z_k^{\ast})\)
#+begin_src jupyter-python
  q = [
      -G_c
      ,* Ω * s_c
      ,* utilities.poly_real(Polynomial.fromroots(
          np.concatenate(
              (
                  [root_f],
                  utilities.except_element(roots_z, i),
                  utilities.except_element(roots_z, i).conj(),
              )
          )
      ))
      for root_f, G_c, γ_c, δ_c, s_c, c_c, i in zip(roots_f, G, γ, δ, s, c, range(len(c)))
  ]
#+end_src

#+RESULTS:

With this we construct our master polynomial \(p = p_1 + \sum_n q_n\).
#+begin_src jupyter-python
  p = p_1 + sum(q)
  p
#+end_src

#+RESULTS:

And find its roots.
#+begin_src jupyter-python
  master_roots = p.roots()
  master_roots
#+end_src

#+RESULTS:
: array([-2.28139877-2.68284887j, -2.28139877+2.68284887j,
:        -0.93979297-0.62025112j, -0.93979297+0.62025112j,
:        -0.24204514-0.10390896j, -0.24204514+0.10390896j,
:        -0.19348529-1.49063362j, -0.19348529+1.49063362j])

Let's see if they're all unique. This should make things easier.
#+begin_src jupyter-python
  np.unique(master_roots).size == master_roots.size
#+end_src

#+RESULTS:
: True

Very nice!

** Calculate the Residuals
These are the prefactors for the diagonal.
#+begin_src jupyter-python
  R_l = f_0(master_roots) / p.deriv()(master_roots)
  R_l
#+end_src

#+RESULTS:
: array([ 0.00251589-0.00080785j,  0.00251589+0.00080785j,
:         0.06323421+0.02800137j,  0.06323421-0.02800137j,
:         0.02732669-0.00211665j,  0.02732669+0.00211665j,
:        -0.09307679+0.36145133j, -0.09307679-0.36145133j])

And the laplace transform of \(\alpha\).
#+begin_src jupyter-python
  def α_tilde(z):
      return (
          -G[None, :]
          ,* ((z[:, None] + γ[None, :]) * s[None, :] + δ[None, :] * c[None, :])
          / (δ[None, :] ** 2 + (γ[None, :] + z[:, None]) ** 2)
      ).sum(axis=1)
#+end_src

#+RESULTS:

And these are for the most compliciated element.
#+begin_src jupyter-python
  R_l_21 = (Ω + α_tilde(master_roots))* f_0(master_roots) / p.deriv()(master_roots)
  R_l_21
#+end_src

#+RESULTS:
: array([-0.00325014-0.02160514j, -0.00325014+0.02160514j,
:         0.00074814-0.05845205j,  0.00074814+0.05845205j,
:        -0.00094159-0.00084894j, -0.00094159+0.00084894j,
:         0.00344359+0.56219924j,  0.00344359-0.56219924j])


Now we can calculate \(G\).
#+begin_src jupyter-python
  def G_12_ex(t):
      t = np.asarray(t)
      t_shape = t.shape

      if len(t.shape) == 0:
          t = t.reshape((1,))
      return Ω * (R_l[None, :] * np.exp(t[:, None] * master_roots[None, :])).real.sum(
          axis=1
      ).reshape(t_shape)


  def G_11_ex(t):
      t = np.asarray(t)
      t_shape = t.shape
      if len(t.shape) == 0:
          t = np.array([t])

      return (
          R_l[None, :]
          ,* master_roots[None, :]
          ,* np.exp(t[:, None] * master_roots[None, :])
      ).real.sum(axis=1).reshape(t_shape)


  def G_21_ex(t):
      t = np.asarray(t)
      return -(R_l_21[None, :] * np.exp(t[:, None] * master_roots[None, :])).real.sum(
          axis=1
      )


  def G_21_ex_alt(t):
      t = np.asarray(t)
      return (
          R_l[None, :]
          ,* master_roots[None, :] ** 2
          ,* np.exp(t[:, None] * master_roots[None, :])
      ).real.sum(axis=1) / Ω


  def G_ex(t):
      t = np.asarray(t)
      if t.size == 1:
          t = np.array([t])
      diag = G_11_ex(t)
      return (
          np.array([[diag, G_12_ex(t)], [G_21_ex(t), diag]]).swapaxes(0, 2).swapaxes(1, 2)
      )


  def G_ex_new(t):
      t = np.asarray(t)
      if t.size == 1:
          t = np.array([t])

      g_12 = R_l[None, :] * np.exp(t[:, None] * master_roots[None, :])
      diag = master_roots[None, :] * g_12
      g_21 = master_roots[None, :] * diag / Ω

      return (
          np.array([[diag, g_12 * Ω], [g_21, diag]])
          .real.sum(axis=3)
          .swapaxes(0, 2)
          .swapaxes(1, 2)
      )
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  #plt.plot(τ, G_ex(τ).reshape(len(τ), 4))
  plt.plot(τ, G_ex_new(τ).reshape(len(τ), 4))
  #plt.plot(ts, proper.reshape(len(ts), 4), linestyle="--")
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fcc1162c8b0> | <matplotlib.lines.Line2D | at | 0x7fcc1162cdc0> | <matplotlib.lines.Line2D | at | 0x7fcc10d36760> | <matplotlib.lines.Line2D | at | 0x7fcc10d36940> |
[[file:./.ob-jupyter/caa8f4b42f6de38f8733ac5169a063de1d5c5d53.svg]]
:END:

#+begin_src jupyter-python
  plt.plot(τ, G_21_ex_alt(τ) - G_21_ex(τ))
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fcc10516610> |
[[file:./.ob-jupyter/5cf02217f725f3e631bb65db0df160a7584ebf65.svg]]
:END:

* Correlation Function
We use the terminology form my notes.
#+begin_src jupyter-python
  def B_s(t):
      return Ω * (R_l * np.exp(t * master_roots)).real.sum()


  def A_s(t):
      return (
          R_l
          ,,* master_roots
          ,,* np.exp(t * master_roots)
      ).real.sum()


  A = G_11_ex
  B = G_12_ex
  n = 1
  q_s_0 = 1 + 2 * n
  p_s_0 = 1 + 2 * n
  qp_0 = 1j
#+end_src

#+RESULTS:

We'd like to calculate \(\langle q(t), q(s)\rangle\).

The convolution free part is easy.
#+begin_src jupyter-python
  def q_corr_conv_free(t, s):
      return (
          q_s_0 * A(t) * A(s)
          + p_s_0 * B(t) * B(s)
          + qp_0 * A(t) * B(s)
          + qp_0.conjugate() * B(t) * A(s)
      )
#+end_src

#+RESULTS:

Now the convolution part (we don't really need this actually).
#+begin_src jupyter-python
  def Λ_real(t, s):
      def integrand(l, r):
          return (B(t - l) * B(s - r) * α(l - r)).real

      return scipy.integrate.dblquad(
          integrand, 0, s, lambda _: 0, lambda _: t
      )[0]

  def Λ_real_mag(t, s):
      def integrand(l, r):
          return (B(t - l) * B(s - r) * α(l - r)).imag

      return scipy.integrate.dblquad(
          integrand, 0, s, lambda _: 0, lambda _: t
      )[0]
#+end_src

#+RESULTS:

With this we can define the correlation function.
#+begin_src jupyter-python
  def q_ts(t, s):
      return q_corr_conv_free(t, s) + Λ(t, s)
#+end_src

#+RESULTS:

For the energy flow we have to calculate the convolution with
\(\dot{\alpha_0}\).  We can reduce the triple integral to a double
integral by assuming the zero temp BCF to be a sum of exponentials.


* Energy Flow
This is "just" another convolution.
We have to integrate over a wedge.
#+begin_src jupyter-python
  from functools import cache, wraps

  @cache
  def flow(t):
        conv_free = quadpy.quad(
            lambda s: (q_corr_conv_free(t, s) * α_0_dot(t - s)).imag, 0, t
        )[0]

        def conv_integrand(l, r, s):
            return (B(t - l) * B(s - r) * α(l - r) * α_0_dot(t - s)).imag

        # def conv_integrand_quadpy(x):
        #     s, r, l = x
        #     if (s<r).any():
        #         print("ooo")
        #         return (B(t - l) * B(s - r) * α(l - r) * α_0_dot(t-s)).imag

        conv_part = scipy.integrate.tplquad(
            conv_integrand, 0, t, lambda s: 0, lambda s: s, lambda _, __: 0, lambda _, __: t
        )[0]

        return -1 / 2 * (conv_part + conv_free), -1/2 * conv_part, -1/2 * conv_free


    # conv_part = scheme.integrate(
    #     conv_integrand_quadpy,
    #     wedge,
    # )
    # r, l, s = generate_wedge_samples(t, int(1e6))
    # conv_part = t**3 / 2* conv_integrand(l, r, s).mean()
    # scheme = quadpy.w3.kubatko_yeager_maggi_9()
    # wedge_base = np.array([[0, 0, 0], [t, 0, 0], [t, t, 0]])
    # wedge_top = wedge_base.copy()
    # wedge_top[:, 2] = t
    # wedge = np.array([wedge_base, wedge_top])
    # return scheme.show(wedge, backend="mpl")
#+end_src

#+RESULTS:


#+begin_src jupyter-python
  τ = np.linspace(0, 16, 5000)
  from multiprocessing import Pool
  with Pool(8) as p:
      flow_τ = np.array(p.map(flow, τ))
#+end_src

#+RESULTS:
:RESULTS:
: Process ForkPoolWorker-1:
: Process ForkPoolWorker-2:
: Process ForkPoolWorker-4:
: Process ForkPoolWorker-3:
# [goto error]
#+begin_example
  ---------------------------------------------------------------------------
  KeyboardInterrupt                         Traceback (most recent call last)
  /tmp/ipykernel_15334/4272752106.py in <module>
        2 from multiprocessing import Pool
        3 with Pool(4) as p:
  ----> 4     flow_τ = np.array(p.map(flow, τ))

  ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in map(self, func, iterable, chunksize)
      362         in a list that is returned.
      363         '''
  --> 364         return self._map_async(func, iterable, mapstar, chunksize).get()
      365 
      366     def starmap(self, func, iterable, chunksize=None):

  ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in get(self, timeout)
      763 
      764     def get(self, timeout=None):
  --> 765         self.wait(timeout)
      766         if not self.ready():
      767             raise TimeoutError

  ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in wait(self, timeout)
      760 
      761     def wait(self, timeout=None):
  --> 762         self._event.wait(timeout)
      763 
      764     def get(self, timeout=None):

  ~/.conda/envs/hops/lib/python3.9/threading.py in wait(self, timeout)
      572             signaled = self._flag
      573             if not signaled:
  --> 574                 signaled = self._cond.wait(timeout)
      575             return signaled
      576 

  ~/.conda/envs/hops/lib/python3.9/threading.py in wait(self, timeout)
      310         try:    # restore state no matter what (e.g., KeyboardInterrupt)
      311             if timeout is None:
  --> 312                 waiter.acquire()
      313                 gotit = True
      314             else:

  KeyboardInterrupt: 
#+end_example
:END:

#+begin_src jupyter-python
  plt.plot(τ, flow_τ)
#+end_src

#+RESULTS:



#+begin_src jupyter-python
  scipy.integrate.trapz(flow_τ, τ)
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  with open("res.npy", 'wb') as f:
      np.save(f, τ)
      np.save(f, flow_τ)
#+end_src


#+begin_src jupyter-python
  with open("res.npy", 'rb') as f:
      τ = np.load(f)
      flow_τ = np.load(f)
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  with utilities.hiro_style():
      plt.plot(τ, flow_τ)
#+end_src

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


* Scratch
#+begin_src jupyter-python
  from hops.util import bcf
  α_nonzero = bcf.OhmicBCF_nonZeroTemp(1, 1, 5, 1/2)
  W, G_raw = utilities.fit_α(α_nonzero, 6, 150, 10_000)
  τ = np.linspace(0, t_max, 1000)
#+end_src

#+RESULTS:

Looks quite good.
#+begin_src jupyter-python
  fig, ax = utilities.plot_complex(τ, α_nonzero(τ), label="exact")
  utilities.plot_complex(
      τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--"
  )
#+end_src

#+RESULTS:
:RESULTS:
| hline | <AxesSubplot:> |
[[file:./.ob-jupyter/f0e3aad0ce345f9d66f2435a7f1dd6f85afa128b.svg]]
:END: