#@title Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Explanation: Copyright 2020 Google
End of explanation
"""


try:
    import recirq
except ImportError:
    print("Installing ReCirq...")
    !pip install git+https://github.com/quantumlib/recirq --quiet
    print("Installed ReCirq!")

"""
Explanation: Fermi-Hubbard experiment example
<table class="tfo-notebook-buttons" align="left">
  <td>
    <a target="_blank" href="https://quantumai.google/cirq/experiments/fermi_hubbard/experiment_example"><img src="https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png" />View on QuantumAI</a>
  </td>
  <td>
    <a target="_blank" href="https://colab.research.google.com/github/quantumlib/ReCirq/blob/master/docs/fermi_hubbard/experiment_example.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png" />Run in Google Colab</a>
  </td>
  <td>
    <a target="_blank" href="https://github.com/quantumlib/ReCirq/blob/master/docs/fermi_hubbard/experiment_example.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/github_logo_1x.png" />View source on GitHub</a>
  </td>
  <td>
    <a href="https://storage.googleapis.com/tensorflow_docs/ReCirq/docs/fermi_hubbard/experiment_example.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/download_icon_1x.png" />Download notebook</a>
  </td>
</table>

This notebook demonstrates how to define, execute, and plot the results of a single instance of the Fermi-Hubbard experiment. We show how to run the experiment using a Cirq simulator and a quantum processor through Google's Quantum Computing Service.
The Fermi-Hubbard model on a one-dimensional lattice of $L$ sites with open boundary conditions is defined by the Hamiltonian
$$
H = - J \sum_{j = 1}^{L - 1} \sum_{\nu = \uparrow, \downarrow} c_{j, \nu}^\dagger c_{j + 1, \nu} + \text{h.c.} + U \sum_{j = 1}^{L} n_{j, \uparrow} n_{j, \downarrow} + \sum_{j = 1}^{L} \sum_{\nu = \uparrow, \downarrow} \epsilon_{j, \nu} n_{j, \nu} 
$$
where $c_{j, \nu}$ ($c_{j, \nu}^\dagger$) are the fermionic annihilation (creation) operators associated to site number $j$ and spin state $\nu$, and $n_{j, \nu} = c_{j, \nu}^\dagger c_{j, \nu}$ are the number operators. The hopping term with coefficient $J$ describes particles tunneling between neighboring sites, the onsite interaction term with coefficient $U$ introduces an energy difference for doubly occupied sites, and the term $\epsilon_{j, \nu}$ represents spin-dependent local potentials. 
Our goal in this experiment is to compute the charge and spin densities which are defined as the sum and difference of the spin-up and spin-down particle densities, respectively
$$
\rho_{j}^{\pm} = \langle n_{j, \uparrow} \rangle \pm \langle n_{j, \downarrow} \rangle,
$$
after simulating the Fermi-Hubbard model for some evolution time. Here, the expectation is taken with respect to the extended Hamiltonian
$$
H' = H + V \sum_{j = 1}^{L - 1} \sum_{\nu = \uparrow, \downarrow}  n_{j, \nu} n_{j + 1, \nu}
$$
which has an additional interaction term between neighboring fermionic sites. This enables extended simulations beyond the standard Fermi-Hubbard model.
Setup
We first install ReCirq which contains code for running Fermi-Hubbard experiments.
End of explanation
"""


try:
    import ipywidgets
except ImportError:
    !pip install ipywidgets --quiet
    !jupyter nbextension enable --py widgetsnbextension --sys-prefix

"""
Explanation: To track the progress of simulating experiments, we use the tqdm package.
End of explanation
"""


import cirq

from recirq import fermi_hubbard
from recirq.fermi_hubbard import publication

# Hide numpy warnings
import warnings
warnings.filterwarnings("ignore")

"""
Explanation: We can now import Cirq and the fermi_hubbard module from ReCirq.
End of explanation
"""


"""Get all layouts for 8 sites on a 23-qubit subgrid of the Google Rainbow processor."""
layouts = publication.rainbow23_layouts(sites_count=8)
print(f"There are {len(layouts)} total qubit layouts.")

"""
Explanation: Experiment parameters
The first step is to decide on exact experiment parameters including problem Hamiltonian, initial state description, as well as a mapping from fermions to qubits on the device. Once we have this information, we can create circuits and run the experiment.
Qubit layout
We will simulate the Fermi-Hubbard model on $L = 8$ sites. Each site is represented by two qubits due to the two spin states, so we need a total of $16$ qubits to simulate the experiment.
The function rainbow23_layouts returns a set of $16$-qubit subgrids of the Google Rainbow processor. 

Note: We use multiple layouts to average results over different qubit assignments. One the quantum processor, this cancels some of the statistical errors which occur from calibration to calibration.
End of explanation
"""


"""Display an example layout."""
print(layouts[0].text_diagram())

"""
Explanation: We can see an example layout by printing out its text diagram.
End of explanation
"""


"""Get FermiHubbardParameters (problem descriptions) for each qubit layout with the above parameters."""
parameters = [
    publication.trapping_instance(
        layout, u=2, dt=0.3, up_particles=2, down_particles=2
    ) 
    for layout in layouts
]

"""
Explanation: The layout indicates the site index $j$ and spin state $\nu$, as well as which cirq.GridQubit on the Rainbow processor this combination of $(j, \nu)$ is encoded into. One can choose a different layout in the previous cell to see how the configurations vary.
Problem parameters
Let's use the Hamiltonian with uniform $J = 1$ and $U = 2$ on each site, initial state prepared as a ground state of a non-interacting Hamiltonian with trapping potential of a Gaussian shape, Trotter step size equal to 0.3, and two particles per chain. The problem parameters with this initial state can be prepared with the pre-defined function trapping_instance.
End of explanation
"""


"""Display the Hamiltonian for an example problem description."""
parameters_example = parameters[0]
print(parameters_example.hamiltonian)

"""
Explanation: Other configurations which support site-dependent $U$ and $J$ coefficients can be prepared by creating instances of the fermi_hubbard.FermiHubbardParameters data class explicitly.

The results are instances of the FermiHubbardParameters data class for each layout. This data class uniquely defines the configuration to run and contains information such as the Hamiltonian, initial state, layout, and time step. Below, we display these values for an example element of parameters.
End of explanation
"""


parameters_example.initial_state

"""
Explanation: We can also see the initial state:
End of explanation
"""


parameters_example.dt

"""
Explanation: And the time step:
End of explanation
"""


"""Create circuits from a problem description."""
initial, trotter, measurement = fermi_hubbard.create_circuits(parameters_example, trotter_steps=1)

"""
Explanation: Circuits
One can directly run an experiment from a FermiHubbardParameters instance (which we will do in the next section). However, it is illustrative to construct the circuits to see how the Fermi-Hubbard execution works.
Circuit creation
To create a circuit from a description of a problem, the function fermi_hubbard.create_circuits can be used. This function inputs a FermiHubbardParameters instance (i.e., a problem description) and number of Trotter steps. It returns circuits for constructing the initial state, simulating time-evolution via a number of Trotter steps, and measuring to compute observables.
End of explanation
"""


"""Display the total circuit to execute."""
circuit = initial + trotter + measurement
circuit

"""
Explanation: Below, we display the complete circuit to execute which is a sum of the three component circuits above.
End of explanation
"""


"""Convert the circuit to native hardware gates perfectly (without noise)."""
publication.ideal_sqrt_iswap_converter().convert(circuit)

"""
Explanation: Note: For a deeper explanation of these circuits and the gates used in them, see the Fermi-Hubbard experiment paper.

Circuit decomposition
The circuit above is constructed using gates which are not native to Google hardware, for example cirq.FSim or cirq.CZ with arbitrary exponent. To run these circuits on Google hardware, we have to convert them into native operations. For the Fermi-Hubbard experiment, a special converter called ConvertToNonUniformSqrtIswapGates is provided.
This converter has the ability to decompose gates to $\sqrt{\small \mbox{iSWAP}}$ both perfectly (i.e., without noise) and with unitary parameters deviating from the perfect ones and varying between qubit pairs. The function ideal_sqrt_iswap_converter creates an instance of the noiseless converter which decomposes $\sqrt{\small \mbox{iSWAP}}$ gates exactly as cirq.FSim(π/4, 0). The function google_sqrt_iswap_converter creates an instance of the noisy converter which approximates the average values on Rainbow processor (which are about cirq.FSim(π/4, π/24) on each two-qubit pair).
Below we show an example of the perfect decomposition into the $\sqrt{\small \mbox{iSWAP}}$ gate set.
End of explanation
"""


"""Set the number of Trotter steps to simulate."""
trotter_steps = range(10 + 1)

"""
Explanation: We will consider both ideal and noisy decompositions when executing the experiment below.
Cirq simulation
This section demonstrates how to simulate experiments using Cirq simulator. We will simulate the evolution from $0$ to $10$ Trotter steps. Physically, this corresponds to an evolution time of $t = 3 \hbar / J$.
End of explanation
"""


"""Get an ideal sampler to simulate experiments."""
ideal_sampler = fermi_hubbard.ConvertingSampler(
    cirq.Simulator(), publication.ideal_sqrt_iswap_converter().convert
)

"""
Explanation: Ideal
As mentioned above, we can use the ideal_sqrt_iswap_converter to convert circuits perfectly into the $\sqrt{\small \mbox{iSWAP}}$ gate set. The Fermi-Hubbard project provides ConvertingSampler that converts circuits before executing ands sampling from them. We get an ideal sampler below.
End of explanation
"""


"""Run the experiments on a perfect simulator for each qubit layout."""
from tqdm.notebook import tqdm

with tqdm(range(len(parameters) * len(trotter_steps))) as progress:
    experiments = [
        fermi_hubbard.run_experiment(
            params, 
            trotter_steps, 
            ideal_sampler, 
            post_run_func=lambda *_: progress.update()
        )
        for params in parameters
    ]

"""
Explanation: We can now run experiments using the run_experiment function. This function takes the parameters of a problem, a sampler, and a list of Trotter steps to simulate. Below, we provide the problem parameters defined on each $16$ qubit layout of the Rainbow processor and simulate the experiments using ten Trotter steps and the ideal_sampler.
End of explanation
"""


"""Post-process the experimental data for all qubit layouts."""
bundle = fermi_hubbard.InstanceBundle(experiments)
bundle.cache_exact_numerics()

"""
Explanation: The output of run_experiment is an instance of the ExperimentResult data class. A series of experiments for the same problem instance on different qubit layouts can be post-processed with the help of the InstanceBundle class. This class takes care of averaging results over qubits layouts, re-scaling the data by comparing against a reference run (perfect simulation in this case), and extracting various quantities.
End of explanation
"""


"""Show quantities which can be accessed from an InstanceBundle."""
for quantity_name in bundle.quantities:
    print(quantity_name)

"""
Explanation: A number of quantities of interest can be accessed from an InstanceBundle, as shown below.
End of explanation
"""


"""Example of getting a DataFrame from a quantity."""
charge_spin_density, _, _ = fermi_hubbard.quantity_data_frame(bundle, "charge_spin_density")
charge_spin_density.head()

"""
Explanation: Each quantity can be converted to a pandas DataFrame using the quantity_data_frame function. Our main goal in simulating the Fermi-Hubbard model was to compute the charge and spin densities
$$
\rho_{j}^{\pm} = \langle n_{j, \uparrow} \rangle \pm \langle n_{j, \downarrow} \rangle .
$$
We can get a DataFrame for the "charge_spin_density" quantity as follows.
End of explanation
"""


"""Plot the charge spin density."""
fermi_hubbard.plot_quantity(bundle, "charge_spin_density");

"""
Explanation: This data frame contains the value, standard error, and standard deviation of the "charge_spin_density" quantity at each site for each time (Trotter step). For convenience, this quantity (and others) can be plotted with the fermi_hubbard.plot_quantity helper function.
End of explanation
"""


"""Plot the charge spin spreading."""
fermi_hubbard.plot_quantity(bundle, "charge_spin_spreading");

"""
Explanation: This plotting function automatically adjusts the appearance of plots according to the data being plotted. We illustrate this by plotting the "charge_spin_spreading" below.
End of explanation
"""


"""Run the experiments on a noisy simulator for each qubit layout."""
parasitic_sampler = fermi_hubbard.ConvertingSampler(
    cirq.Simulator(), publication.google_sqrt_iswap_converter().convert
    )

with tqdm(range(len(parameters) * len(trotter_steps))) as progress:
    experiments = [
        fermi_hubbard.run_experiment(
            params, 
            trotter_steps, 
            parasitic_sampler, 
            post_run_func=lambda *_: progress.update()
        )
        for params in parameters
    ]

"""
Explanation: One can compare these plots to Figure 2 of the Fermi-Hubbard experiment paper.
Parasitic controlled-phase
We now run the same experiment but with the google_sqrt_iswap_converter. As mentioned, this decomposes $\sqrt{\small \mbox{iSWAP}}$ gates imperfectly as cirq.FSim(π/4, π/24) which is close to the average value of the parasitic controlled phase on the Rainbow processor.
End of explanation
"""


"""Post-process the experimental data for all qubit layouts."""
bundle = fermi_hubbard.InstanceBundle(experiments)
bundle.cache_exact_numerics()

"""
Explanation: As above, we can post-process the data using an InstanceBundle and plot quantities of interest using the plot_quantity helper function.
End of explanation
"""


fermi_hubbard.plot_quantity(bundle, "charge_spin_density");

"""
Explanation: We first plot the "charge_spin_density":
End of explanation
"""


"""Plot the charge spin spreading."""
fermi_hubbard.plot_quantity(bundle, "charge_spin_spreading", show_std_error=True);

"""
Explanation: And plot the "charge_spin_spreading" as well.
End of explanation
"""


"""Get an engine sampler."""
import os
import cirq_google

if "GOOGLE_CLOUD_PROJECT" in os.environ:
    engine_sampler = cirq_google.get_engine_sampler(
        processor_id="rainbow", gate_set_name="sqrt_iswap"
    )
else:
    # Use the simulator as a backup.
    engine_sampler = cirq.Simulator()

# Get a sampler for the Fermi-Hubbard experiment.
google_sampler = fermi_hubbard.ConvertingSampler(
    engine_sampler, publication.google_sqrt_iswap_converter().convert
)

"""
Explanation: One can compare these to the simulation with exact decompositions above to see the effect of the parasitic controlled phase.
Execution on Google's Quantum Computing Service
In order to run an experiment on Google's QCS, a QuantumEngine sampler is needed. To create an engine sampler, an environment variable GOOGLE_CLOUD_PROJECT must be present and set to a valid Google Cloud Platform project identifier.
End of explanation
"""


"""Run the experiments on Google's QCS and save the results."""
# Directory to save results in.
results_dir = "trapping"

with tqdm(range(len(layouts) * len(trotter_steps))) as progress:
    for index, params in enumerate(parameters):
        experiment = fermi_hubbard.run_experiment(
            params, 
            trotter_steps, 
            google_sampler, 
            post_run_func=lambda *_: progress.update()
        )
        
        fermi_hubbard.save_experiment(
            experiment, f"{results_dir}/trapping_{index + 1}.json"
        )

"""
Explanation: Now that we are running on a quantum computer, we follow good experimental practice and save the results on disk as soon as each experiment finishes using the fermi_hubbard.save_experiment function. Although rare, remote operation may fail for various reasons. More advanced execution workflow might include error handling, experiment pause and continuation, etc., which we omit here for simplicity.

Note: We do not include Floquet calibration, the calibration technique described in the Supplementary Information of the Fermi-Hubbard experiment paper, when executing the experiments below.
End of explanation
"""


"""Load experimental results."""
experiments = [
    fermi_hubbard.load_experiment(f"{results_dir}/trapping_{index + 1}.json") 
    for index in range(len(parameters))
]

"""
Explanation: We can now load the results using fermi_hubbard.load_experiment.
End of explanation
"""


"""Post-process the experimental data for all qubit layouts."""
bundle = fermi_hubbard.InstanceBundle(
    experiments,numerics_transform=publication.parasitic_cphase_compensation(0.138)
)
bundle.cache_exact_numerics()

"""
Explanation: When post-processing experimental data from hardware, we include effects due the the parasitic controlled phase as shown below. The value $\phi = 0.138$ was the approximate value of the parasitic controlled phase at the time when the experimental results in the paper were collected.
End of explanation
"""


"""Plot the charge spin density."""
fermi_hubbard.plot_quantity(bundle, "charge_spin_density", show_std_error=True);

"""Plot the charge spin spreading."""
fermi_hubbard.plot_quantity(bundle, "charge_spin_spreading", show_std_error=True);

"""
Explanation: We can now visualize these results using the same plotting functions from above. Here we show the standard deviation of results in the plots.
End of explanation
"""


