r"""
.. role:: html(raw)
   :format: html

.. _ham_sim:

Hamiltonian simulation
======================

    "The problem of simulating the dynamics of quantum systems was the original motivation for quantum
    computers and remains one of their major potential applications." - Berry et al. [[1]_]

The simulation of atoms, molecules and other biochemical systems is another application uniquely
suited to quantum computation. For example, the ground state energy of large systems, the dynamical
behaviour of an ensemble of molecules, or complex molecular behaviour such as protein folding, are
often computationally hard or downright impossible to determine via classical computation or
experimentation [[2]_][[3]_].

In the discrete-variable qubit model, efficient methods of Hamiltonian simulation have been
discussed at-length, providing several implementations depending on properties of the Hamiltonian,
and resulting in a linear simulation time [[4]_][[5]_].
Efficient implementations of Hamiltonian simulation also exist in the CV formulation
[[6]_], with specific application to `Bose-Hubbard Hamiltonians
<https://en.wikipedia.org/wiki/Bose%E2%80%93Hubbard_model>`_ (describing a system of interacting
bosonic particles on a lattice of orthogonal position states [[7]_]). As
such, this method is ideally suited to photonic quantum computation.

Implementation
--------------

For a quick example, consider a lattice composed of two adjacent nodes:


.. raw:: html

    <br>

.. image:: /tutorials/images/graph.svg
    :align: center
    :width: 40%
    :target: javascript:void(0);


.. raw:: html

    <br>

This graph is represented by the :math:`2\times 2` adjacency matrix
:math:`A=\begin{bmatrix}0&1\\1&0\end{bmatrix}`. Here, each node in the graph represents a qumode, so
we can model the dynamics of Bosons on this structure via a 2-qumode CV circuit.

The `Bose-Hubbard Hamiltonian <https://en.wikipedia.org/wiki/Bose%E2%80%93Hubbard_model>`_ with
on-site interactions is given by

.. math::

    H = J\sum_{i}\sum_j A_{ij} \ad_i\a_j + \frac{1}{2}U\sum_i \hat{n}_i(\hat{n}_i-1)= J(\ad_1 \a_2 +
    \ad_2\a_1) + \frac{1}{2}U (  \hat{n}_1^2 - \hat{n}_1 + \hat{n}_2^2 - \hat{n}_2)

where :math:`J` represents the transfer integral or hopping term of the boson between nodes, and
:math:`U` is the on-site interaction potential. Here, :math:`\ad_1 \a_2` represents a boson
transitioning from node 1 to node 2, while :math:`\ad_2\a_1` represents a boson transitioning from
node 2 to node 1, and :math:`\hat{n}_i=\ad_i\a_i` is the number operator applied to mode :math:`i`.
Applying the Lie-product formula, we find that

.. math::
    e^{-iHt} = \left[\exp\left({-i\frac{ J t}{k}(\ad_1 \a_2 +
    \ad_2\a_1)}\right)\exp\left(-i\frac{Ut}{2k}\hat{n}_1^2\right)\exp\left(-i\frac{Ut}{2k}\hat{n}_2^2\right)
    \exp\left(i\frac{Ut}{2k}\hat{n}_1\right)\exp\left(i\frac{Ut}{2k}\hat{n}_2\right)\right]^k+\mathcal{O}\left(t^2/k\right),

where :math:`\mathcal{O}\left(t^2/k\right)` is the order of the error term, derived from the Lie
product formula. Comparing this to the form of various CV gates, we can write this as the product of
:ref:`symmetric beamsplitters <beamsplitter>` (:class:`~strawberryfields.ops.BSgate`), :ref:`Kerr gates <kerr>`
(:class:`~strawberryfields.ops.Kgate`), and :ref:`rotation gates <rotation>`
(:class:`~strawberryfields.ops.Rgate`):

.. math::
    e^{iHt} = \left[BS\left(\theta,\phi\right)\left(K(r)R(-r)\otimes
    K(r)R(-r)\right)\right]^k+\mathcal{O}\left(t^2/k\right).

where :math:`\theta=-Jt/k`, :math:`\phi=\pi/2`, and :math:`r=-Ut/2k`.

For the case :math:`k=2`, this can be drawn as the circuit diagram

|

.. image:: /tutorials/images/hamsim.svg
    :align: center
    :width: 70%
    :target: javascript:void(0);

|

For more complex CV decompositions, including those with interactions, see Kalajdzievski et al.
[[6]_] for more details.

Code
----

The Hamiltonian simulation circuit displayed above for the 2-node lattice, can be implemented using
Strawberry Fields:
"""
import numpy as np

# set the random seed
np.random.seed(42)

import strawberryfields as sf
from strawberryfields.ops import *

ham_simulation = sf.Program(2)

# set the Hamiltonian parameters
J = 1           # hopping transition
U = 1.5         # on-site interaction
k = 20          # Lie product decomposition terms
t = 1.086       # timestep
theta = -J*t/k
r = -U*t/(2*k)

with ham_simulation.context as q:
    # prepare the initial state
    Fock(2) | q[0]

    # Two node tight-binding
    # Hamiltonian simulation

    for i in range(k):
        BSgate(theta, np.pi/2) | (q[0], q[1])
        Kgate(r)  | q[0]
        Rgate(-r) | q[0]
        Kgate(r)  | q[1]
        Rgate(-r) | q[1]

######################################################################
# where, for this example, we have set ``J=1``, ``U=1.5``, ``k=20``, ``t=1.086``, ``theta =
# -J*t/k``, and ``r = -U*t/(2*k)``.
#
# Since are initial state is non-Gaussian, we will initialize the ``"fock"``
# backend, with a cutoff of 5:

eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})

######################################################################
# After constructing the circuit and running the engine,

results = eng.run(ham_simulation)

######################################################################
# the site occupation probabilities can be calculated via the
# ``fock_prob`` state method, which returns the probability of
# measuring some output Fock state:

state = results.state
print("P(|0, 2>) = ", state.fock_prob([0, 2]))
print("P(|1, 1>) = ", state.fock_prob([1, 1]))
print("P(|2, 0>) = ", state.fock_prob([2, 0]))

######################################################################
# As Hamiltonian simulation is particle preserving, these probabilities should add up to one:

result = [state.fock_prob([0, 2]), state.fock_prob([1, 1]), state.fock_prob([2, 0])]
print(np.sum(result))


######################################################################
# We can compare this result to the analytic matrix exponential :math:`e^{-iHt}`, where the matrix
# elements of :math:`H` can be computed in the Fock basis. Considering the diagonal interaction
# terms,
#
# .. math::
#     & \braketT{0,2}{H}{0,2} = \frac{1}{2}U\braketT{0}{(\hat{n}^2-\hat{n})}{0} +
#     \frac{1}{2}U\braketT{2}{(\hat{n}^2-\hat{n})}{2} = \frac{1}{2}U(2^2-2) = U\\[7pt] &
#     \braketT{1,1}{H}{1,1} = 0\\[7pt] & \braketT{2,0}{H}{2,0} = U
#
# as well as the off-diagonal hopping terms,
#
# .. math::
#     & \braketT{1,1}{H}{0,2} = J\braketT{1,1}{\left(\ad_1\a_2 + \a_1\ad_2\right)}{0,2} =
#     J(\sqrt{1}\sqrt{2} + \sqrt{0}\sqrt{3}) = J\sqrt{2}\\[7pt] & \braketT{1,1}{H}{2,0} = J\sqrt{2}
#
# and taking into account the Hermiticity of the system, we arrive at
#
# .. math::
#     H = \begin{bmatrix}U&J\sqrt{2}&0\\ J\sqrt{2} & 0 & J\sqrt{2}\\ 0 & J\sqrt{2} & U\end{bmatrix}
#
# which acts on the Fock basis :math:`\{\ket{0,2},\ket{1,1},\ket{2,0}\}`. Using the SciPy matrix
# exponential function ``scipy.linalg.expm``:

from scipy.linalg import expm

H = J*np.sqrt(2)*np.array([[0,1,0],[1,0,1],[0,1,0]]) + U*np.diag([1,0,1])
init_state = np.array([0,0,1])

theoretical_result = np.abs(np.dot(expm(-1j*t*H), init_state))**2
print(theoretical_result)

######################################################################
# which agrees, within the expected error margin, with our Strawberry Fields Hamiltonian simulation:

print(np.all(np.abs(theoretical_result - result) < 1e-2))

######################################################################
# References
# ----------
#
# .. [1] Dominic W. Berry, Andrew M. Childs, and Robin Kothari. Hamiltonian simulation with nearly optimal
#        dependence on all parameters. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science.
#        IEEE, Oct 2015. doi:10.1109/focs.2015.54.
#
# .. [2] Alán Aspuru-Guzik, Anthony D. Dutoi, Peter J. Love, and Martin Head-Gordon. Simulated quantum
#        computation of molecular energies. Science, 309(5741):1704–1707, 2005. doi:10.1126/science.1113479.
#
# .. [3] James D. Whitfield, Jacob Biamonte, and Alán Aspuru-Guzik. Simulation of electronic structure
#        Hamiltonians using quantum computers. Molecular Physics, 109(5):735–750, 2011.
#        doi:10.1080/00268976.2011.552441.
#
# .. [4] Andrew M. Childs and Nathan Wiebe. Hamiltonian simulation using linear combinations of unitary
#        operations. Quantum Information and Computation, 12(11-12):901–924, 2012. arXiv:1202.5822.
#
# .. [5] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders. Efficient quantum algorithms
#        for simulating sparse Hamiltonians. Communications in Mathematical Physics, 270(2):359–371, 2006.
#        doi:10.1007/s00220-006-0150-x.
#
# .. [6] Timjan Kalajdzievski, Christian Weedbrook, and Patrick Rebentrost. Continuous-variable gate
#        decomposition for the Bose-Hubbard model. 2018. arXiv:1801.06565.
#
# .. [7] Tomasz Sowiński, Omjyoti Dutta, Philipp Hauke, Luca Tagliacozzo, and Maciej Lewenstein. Dipolar
#        molecules in optical lattices. Physical Review Letters, 108:115301, 2012.
#        doi:10.1103/PhysRevLett.108.115301.