{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear cross-entropy benchmarking using Piquasso\n",
    "\n",
    "\n",
    "Linear cross-entropy benchmarking, or LXEB, is a statistical test used to quantify whether samples produced by a photonic experiment are correlated with the ideal output distribution predicted by quantum theory. In photonic quantum advantage proposals, the experiment typically consists of preparing an input state $\\rho$, evolving it through an $m$-mode interferometer $U$, and measuring the output photon-number pattern <cite data-footcite=\"lxeb\"></cite>.\n",
    "\n",
    "We denote by $p_{\\rho,U}(\\mathbf n)$ the ideal probability of observing an output occupation pattern $\\mathbf n$. If the input state has fixed total photon number $n$, then the possible outcomes are occupation vectors $\\mathbf n \\in \\mathcal S_{m,n}$, where $\\mathcal S_{m,n}$ is the set of all $m$-mode photon-number configurations with total photon number $n$. If the physical input state does not have a sharply fixed photon number, one can instead define the distribution after postselecting on runs with total detected photon number $n$.\n",
    "\n",
    "The actual experimental device samples from another distribution, which we denote by $q$. The purpose of LXEB is to test whether $q$ is correlated with the ideal distribution $p_{\\rho,U}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of linear cross-entropy\n",
    "\n",
    "The linear cross-entropy between the ideal distribution $p_{\\rho,U}$ and the experimental distribution $q$ is\n",
    "$$\n",
    "\\mathrm{LXE}(p_{\\rho,U},q)\n",
    "=\n",
    "\\mathbb E_{\\mathbf n \\sim q}\\, p_{\\rho,U}(\\mathbf n)\n",
    "=\n",
    "\\sum_{\\mathbf n \\in \\mathcal S_{m,n}}\n",
    "q(\\mathbf n)\\,p_{\\rho,U}(\\mathbf n).\n",
    "$$\n",
    "Operationally, this quantity asks the following question:\n",
    "\n",
    "> If we draw samples from the experimental device, how large are their ideal probabilities?\n",
    "\n",
    "Thus, LXEB rewards experimental samples that land on outcomes that are assigned high probability by the ideal quantum distribution. If $q$ is uncorrelated with $p_{\\rho,U}$, then the LXE score is expected to be comparatively small. If $q = p_{\\rho,U}$, meaning that the experiment perfectly samples from the ideal distribution, then\n",
    "\n",
    "$$\n",
    "\\mathrm{LXE}(p_{\\rho,U},p_{\\rho,U})\n",
    "=\n",
    "\\sum_{\\mathbf n \\in \\mathcal S_{m,n}}\n",
    "p_{\\rho,U}(\\mathbf n)^2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical LXEB score\n",
    "\n",
    "In an experiment, we do not know the full distribution $q$ directly. Instead, we collect a finite set of samples\n",
    "$$\n",
    "\\mathcal X = \\{\\mathbf n_1,\\ldots,\\mathbf n_s\\},\n",
    "\\qquad\n",
    "\\mathbf n_i \\sim q.\n",
    "$$\n",
    "The empirical LXEB score is\n",
    "$$\n",
    "\\widehat{\\mathrm{LXE}}(p_{\\rho,U},\\mathcal X)\n",
    "=\n",
    "\\frac{1}{s}\n",
    "\\sum_{i=1}^s\n",
    "p_{\\rho,U}(\\mathbf n_i).\n",
    "$$\n",
    "This estimator is unbiased, since\n",
    "$$\n",
    "\\mathbb E_{\\mathcal X \\sim q}\n",
    "\\left[\\widehat{\\mathrm{LXE}}(p_{\\rho,U},\\mathcal X)\\right]\n",
    "=\n",
    "\\mathrm{LXE}(p_{\\rho,U},q).\n",
    "$$\n",
    "Therefore, the empirical LXEB score estimates the correlation between the experimental samples and the ideal quantum probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LXEB pass criterion\n",
    "\n",
    "Given a tolerance parameter $\\Delta \\in [0,1)$, one can say that the sample set $\\mathcal X$ passes the LXEB test if\n",
    "$$\n",
    "\\frac{\n",
    "\\left|\n",
    "\\widehat{\\mathrm{LXE}}(p_{\\rho,U},\\mathcal X)\n",
    "-\n",
    "\\mathrm{LXE}_{\\mathrm{ref}}^{(n)}(\\rho)\n",
    "\\right|\n",
    "}{\n",
    "\\mathrm{LXE}_{\\mathrm{ref}}^{(n)}(\\rho)\n",
    "}\n",
    "<\n",
    "\\Delta.\n",
    "$$\n",
    "In words, the empirical LXEB score must approximate the Haar-averaged ideal reference value $\\mathrm{LXE}_{\\mathrm{ref}}^{(n)}(\\rho)$ up to relative error $\\Delta$.\n",
    "This criterion tests whether the experimental data have the same average ideal-probability weight as samples from the target distribution. A successful LXEB test therefore provides evidence that the device is producing samples correlated with the intended photonic quantum distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LXEB experiment for boson sampling with Piquasso\n",
    "\n",
    "In this notebook, we simulate a small boson sampling experiment and evaluate its\n",
    "linear cross-entropy benchmarking score.\n",
    "\n",
    "The experiment has three parts:\n",
    "\n",
    "1. Prepare an input Fock state with a fixed number of photons.\n",
    "2. Apply a Haar-random interferometer and sample output photon-number patterns.\n",
    "3. Compute the empirical LXEB score\n",
    "\n",
    "$$\n",
    "\\widehat{\\mathrm{LXE}}\n",
    "=\n",
    "\\frac{1}{s}\n",
    "\\sum_{i=1}^s p_U(\\mathbf n_i),\n",
    "$$\n",
    "\n",
    "where the samples $\\mathbf n_i$ are produced by the simulated experiment and\n",
    "$p_U(\\mathbf n_i)$ is the ideal probability of observing that sample.\n",
    "\n",
    "We compare the empirical score with the Haar-averaged reference value $\\mathrm{LXE}_{\\mathrm{ref}}^{(n)}(\\rho)$ computed\n",
    "using `piquasso.lxeb.lxe_ref_boson_sampling` in the following script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated value: 2.3201076223045534e-12\n",
      "Reference value: 2.3130728523936438e-12\n",
      "Relative error: 0.0030413092711843275\n"
     ]
    }
   ],
   "source": [
    "import scipy\n",
    "import piquasso as pq\n",
    "\n",
    "from piquasso.lxeb import lxe_ref_boson_sampling\n",
    "\n",
    "m = 32\n",
    "n = 16\n",
    "\n",
    "config = pq.Config(seed_sequence=42)\n",
    "shots = 2000\n",
    "\n",
    "U = scipy.stats.unitary_group.rvs(m, random_state=config.seed_sequence)\n",
    "\n",
    "instructions = [pq.StateVector([1] * n + [0] * (m - n)), pq.Interferometer(U)]\n",
    "\n",
    "state = pq.simulate(pq.Program(instructions), number_of_modes=m).state\n",
    "\n",
    "samples = pq.simulate(\n",
    "    pq.Program(instructions + [pq.ParticleNumberMeasurement()]),\n",
    "    number_of_modes=m,\n",
    "    shots=shots,\n",
    ").samples\n",
    "\n",
    "sum_ = 0.0\n",
    "\n",
    "for sample in samples:\n",
    "    sum_ += state.get_particle_detection_probability(sample)\n",
    "\n",
    "lxe_estimate = sum_ / shots\n",
    "print(\"Estimated value:\", lxe_estimate)\n",
    "\n",
    "lxe_reference = lxe_ref_boson_sampling(n, m)\n",
    "print(\"Reference value:\", lxe_reference)\n",
    "\n",
    "relative_error = abs(lxe_reference - lxe_estimate) / lxe_reference\n",
    "print(\"Relative error:\", relative_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whenever there are losses in the experiment, the estimated LXE value dramatically differs from the reference value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lossy estimated value: 7.29166920647751e-15\n",
      "Relative error: 0.9968476266543304\n"
     ]
    }
   ],
   "source": [
    "samples = pq.simulate(\n",
    "    pq.Program(instructions + [pq.UniformLoss(0.85), pq.ParticleNumberMeasurement()]),\n",
    "    number_of_modes=m,\n",
    "    shots=shots,\n",
    ").samples\n",
    "\n",
    "sum_ = 0.0\n",
    "\n",
    "for sample in samples:\n",
    "    sum_ += state.get_particle_detection_probability(sample)\n",
    "\n",
    "lossy_lxe_estimate = sum_ / shots\n",
    "print(\"Lossy estimated value:\", lossy_lxe_estimate)\n",
    "\n",
    "relative_error = abs(lxe_reference - lossy_lxe_estimate) / lxe_reference\n",
    "print(\"Relative error:\", relative_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LXEB is useful because it reduces a complicated distributional comparison to a single scalar statistic. This is especially important in quantum advantage experiments, where the number of possible output patterns\n",
    "$$\n",
    "|\\mathcal S_{m,n}| = \\binom{m+n-1}{n}\n",
    "$$\n",
    "can be extremely large. However, several caveats are important:\n",
    "\n",
    "1. LXEB is not a full certification of the output distribution.\n",
    "   Passing the test shows correlation with the ideal probabilities, but it does not prove that $q = p_{\\rho,U}$.\n",
    "\n",
    "2. The ideal probabilities $p_{\\rho,U}(\\mathbf n_i)$ must be computed for the observed samples.\n",
    "   These probabilities may involve matrix permanents or related quantities, which can be computationally expensive.\n",
    "\n",
    "3. Finite-sample fluctuations matter.\n",
    "   Even if the experiment samples from the ideal distribution, the empirical score fluctuates around its expectation. The sample size $s$ must therefore be large enough for the desired relative precision, but numerical evidence suggests that this procedure remains sample-efficient <cite data-footcite=\"lxeb\"></cite>.\n",
    "\n",
    "In summary, LXEB provides a practical way to benchmark photonic quantum advantage experiments by checking whether experimentally observed samples are concentrated on outputs that the ideal quantum model predicts to be likely."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}