{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "prep-00",
   "metadata": {},
   "source": [
    "# Alias-sampling state preparation\n",
    "\n",
    "**Download Notebook** - {nb-download}`alias_prep.ipynb`\n",
    "\n",
    "Alias sampling loads a classical probability distribution into a quantum index register. It prepares the index **together with workspace**:\n",
    "\n",
    "$$\n",
    "P|0\\rangle=\\sum_j\\sqrt{\\widetilde p_j}|j\\rangle|g_j\\rangle,\n",
    "\\qquad \\Pr(\\text{index}=j)=\\widetilde p_j.\n",
    "$$\n",
    "\n",
    "The probabilities $\\widetilde p_j$ approximate the requested weights at a chosen binary precision. The index need not be a pure state on its own: the workspace can be entangled with it.\n",
    "\n",
    "Run from a source checkout with the development dependencies installed. The repository default is little endian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "prep-01",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from math import ceil, log2\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from guppylang import guppy\n",
    "from guppylang.std.quantum import qubit, discard_array, discard\n",
    "from guppylang.std.debug import state_output\n",
    "\n",
    "from guppyalgos.algorithms.state_preparation.alias_sampling import alias_samp_prep\n",
    "from guppyalgos.algorithms.state_preparation.alias_sampling.alias_table import discretize_distribution\n",
    "from guppyalgos.primitives.subroutines.fanout import fanout_basic\n",
    "from guppyalgos.utils import qarray"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prep-02",
   "metadata": {},
   "source": [
    "## 1. Choose the distribution and precision\n",
    "\n",
    "Use $p=(1/2,1/4,1/4)$. Two index qubits address the three entries, leaving address 3 unused. Two probability bits represent these weights exactly:\n",
    "\n",
    "$$\n",
    "n_{\\mathrm{index}}=\\lceil\\log_2L\\rceil=2,\n",
    "\\qquad n_{\\mathrm{prob}}=\\lceil\\log_2(1/\\epsilon)\\rceil=2,\n",
    "\\qquad \\epsilon=1/4.\n",
    "$$\n",
    "\n",
    "`discretize_distribution` gives the probabilities the finite-precision circuit will load. Smaller `precision` values use more probability bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "prep-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_dist = np.array([0.5, 0.25, 0.25])\n",
    "precision = 0.25\n",
    "n_index_qubits = ceil(log2(len(prob_dist)))\n",
    "n_prob_bits = ceil(log2(1 / precision))\n",
    "loaded_probabilities = discretize_distribution(prob_dist, n_prob_bits)\n",
    "np.testing.assert_allclose(loaded_probabilities, prob_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prep-04",
   "metadata": {},
   "source": [
    "## 2. Build and run PREPARE\n",
    "\n",
    "The alias QROM loads an alternative index and a keep threshold. A uniform comparison register chooses whether to keep the original index or swap in the alternative.\n",
    "\n",
    "| Register | Purpose |\n",
    "| :-- | :-- |\n",
    "| `index_qreg` | The sampled table index. |\n",
    "| `alternative_qreg` | The alternative index loaded by QROM. |\n",
    "| `keep_qreg` | The fixed-point keep threshold. |\n",
    "| `comparison_qreg` | The uniform value compared with that threshold. |\n",
    "| `comparison_result` | The comparison flag controlling the swap. |\n",
    "\n",
    "All registers start at zero. `fanout_op` chooses how QROM writes the selected data bits; here it uses basic CNOT fanout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "prep-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "alias_prepare = alias_samp_prep(prob_dist, precision, fanout_op=fanout_basic)\n",
    "\n",
    "\n",
    "@guppy\n",
    "def main() -> None:\n",
    "    index_qreg = qarray(n_index_qubits)\n",
    "    alternative_qreg = qarray(n_index_qubits)\n",
    "    keep_qreg = qarray(n_prob_bits)\n",
    "    comparison_qreg = qarray(n_prob_bits)\n",
    "    comparison_result = qubit()\n",
    "    alias_prepare(\n",
    "        index_qreg, alternative_qreg, keep_qreg,\n",
    "        comparison_qreg, comparison_result, False,\n",
    "    )\n",
    "    state_output(\"result_state\", index_qreg)\n",
    "\n",
    "    # This diagnostic is finished; no coherent operation follows.\n",
    "    discard_array(index_qreg)\n",
    "    discard_array(alternative_qreg)\n",
    "    discard_array(keep_qreg)\n",
    "    discard_array(comparison_qreg)\n",
    "    discard(comparison_result)\n",
    "\n",
    "\n",
    "total_qubits = (\n",
    "    3 * n_index_qubits + 2 * n_prob_bits + 5\n",
    "    + int(np.ceil((n_index_qubits - 3) / 2))\n",
    ")\n",
    "result = main.emulator(total_qubits).with_seed(42).run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prep-06",
   "metadata": {},
   "source": [
    "## 3. Check the index probabilities\n",
    "\n",
    "Trace out the workspace before comparing with the requested distribution:\n",
    "\n",
    "$$\n",
    "\\Pr(j)=\\sum_g|\\langle j,g|P|0\\rangle|^2.\n",
    "$$\n",
    "\n",
    "The simulator represents the reduced index state as a weighted mixture. Sum each component's squared **magnitude**, weighted by its probability; squaring complex amplitudes directly would be incorrect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "prep-07",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAE3CAYAAADmP0YBAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMWZJREFUeJzt3QeUVOX9//EvndCb9KY06QoRUEBRiv4QjCAC0kEIRVGEYMRQREEMEcGggtI1iRKJIgIiSFOkSZMiRQxIkd5FWOr/fJ7875yZZXeZXWb37rLv1zlzdubObXN3ducz3+e5z01z9erVqwYAAABfpPVnswAAABDCGAAAgI8IYwAAAD4ijAEAAPiIMAYAAOAjwhgAAICPCGMAAAA+IowBAACk9jB24cIFO3r0qMV3/Nlff/3Vzp8/n2j7BQAAcFOHsStXrli/fv0sV65cVrJkSStatKh9+umn111uxYoVVqNGDcufP78VL17cunbt6oIZAABASuNrGBs1apRNmTLFli9fbqdPn7Y///nP1qpVK9u6dWusy2zatMnq169vDz30kJ08edIOHjxodevWte3btyfpvgMAAERCGj+vTalqWMuWLW3kyJGBaaVLl7amTZva6NGjY1zmkUcesSNHjrjqGAAAQErnW2Xs8OHD9vPPP1vt2rVDptepU8dWr14d4zKXLl2yBQsWWIsWLVz/smPHjrmmTgAAgJQqvV8bVnVL8uXLFzJdj9VsGdsy6rB//PhxK1WqlGum1ON27drZmDFjLEuWLDEuFxUV5W4eBTitI2/evJYmTZqIvi4AAABR4ejMmTNWuHBhS5s2bfILY95OqdoV7OLFi5YuXboYl/GC04QJE2zJkiVWoUIF27Fjh9WrV8+yZs0aa9PmiBEjbOjQoRF/DQAAANezd+9ed5JisgtjRYoUcT/VAT/YoUOHAs9Fd8stt1imTJnsiSeecEFMypYta+3bt7eZM2fGGsYGDBhgffv2DTw+deqUOwtTBydHjhwRfFUAAAD/o5MTixUrZtmzZ7e4+BbGFILuuOMO1wdMZ1B6VbKFCxda7969A/OdOHHCjUNWoEABVzHTmZMq+QXT42zZssW6LQU43WLaB8IYAABITNfrEuVbGJNBgwa5IHbXXXfZ3Xffba+//rqb3rNnz8A8/fv3t5UrV9rmzZvd48GDB7thLRTK1Pl/1apVNnXq1FirYgAAAMmZr+OMNW/e3P7xj3/YtGnTrFmzZq6ct3TpUtcc6cmTJ48VLFgw8FghTE2SWq5x48ZunLKJEyda9+7dfXoVAAAAKXScMb8o9OXMmdP1HaOZEgAA+Jk3fG2mBADE7vLly+4McwDJk/qyp0+f/oaHySKMAUAypOvt7tu3z41TBCD50hinhQoVsowZMyZ4HYQxAEiGFTEFMf2TVx9aBqcGkh99UdJoDxqQfteuXVamTJk4B3aNC2EMAJIZNU3qH72C2O9+9zu/dwdALPT3mSFDBnd5RwWzzJkzW4o7mxIAEDsqYkDyl9BqWMg6IrInAAAASBDCGAAASWjFihXuOoXqG4j4+dvf/maPP/54xOf1G33GACCFKPnCnCTd3u7XHg5rPoWKEiVKxDlPy5Yt7Y033jC/rF+/3po2bWo7duxwJ0b4uZ6oqCjbv38/Z8omwKlTp1yH+UjP6zfCGADghsda0mXrPB9++KE9//zztnfv3sC0rFmz+nqUvQB05cqVZLEeJMzzzz9vzz77bMTn9RvNlACAG6ZmN++WO3fua6Y98sgj7mfx4sXdtYhHjBhxzYC2el6XumvRooWVKlXK/v73v7vpq1evtgYNGli5cuXsD3/4g3366adu3uCqx4EDB6xHjx5Wvnx5u/POO+25555zo5/L7t273fZF69CyzzzzzDWv4eDBg4H9vfXWW61+/fr2n//8J/B8uOvx9qdXr15WsWJFq169urt+cvQA980337jL+t1+++3udf34448hz2/dutU1s5UuXdqtY+TIkSHr0PEbMmSIe65SpUrusoCHDx8O65gEN5cuXrw4zv0I5i0zf/58e/DBB91wDo899pg7NtHn+eqrr+y+++6zkiVL2nfffeeeW7ZsmT388MPuNd177702YcKEayqEmkfHWfM0bNjQvv7668Bz7777rjuuwb+T9u3bW4UKFdz7atSoUYHm3+jzyowZM9x2b7vtNvf7XbBgQcjzL730knXo0MHefPNNq1Onjjuuffv2tXPnzlliIowBABLdxx9/7Kpn+qAdOnSouyaxPviCqeLUu3dvFwj0Qd65c2cXuBTE9KGvYNS2bVvr2rWrm9f70FVzVO3atd1ZbR999JFNnjzZdu7c6ZoT9UFfrFgxmzRpkptXwUP7oX2ITkOJ6DndtP2OHTu6m/eBHe56jh8/7oLBTz/9ZO+9954LHKoSfvHFFyHz9evXz1Vv9Lq0nwqhnhMnTrjQoKbQTz75xAYPHuyChn56FM5UhVTQU8ioWbNm4Jhe75gEV/ni2o/ovGUU/BRS/v3vf7shHRTMLl26FDJPnz593P4qdFatWtUWLVrktq9rUc+ZM8cGDBhgw4cPd/vvUchTSNL8Ct0vvPBCyPsketOjwqreB9qPt956y4VR3Y9p3s8//9y9f1q1amWzZ8+2Bx54wIXQVatWBeY5efKk/fOf/7SNGzfaO++849ap966+PCQmmikBAImuYMGCgfuqjilI9OzZ030YB1P1RpUOj/qZaVl9MGqoD1UqNKaTwoNn/Pjxli9fPjePRx+oefLkcX28qlWr5oKWFC5c2LJlyxZrc6sqOh5V5zZs2GBTp051FRo9H856xo0b5wLKzJkzA+PEaR+iV8YU0lTVkldeecXuuOMOO3TokBUoUMDefvtt17SrEKXtVqlSxQUFHTMFFG1706ZN1qhRIxfaRJWtTp06hX1MwtmP2Oj3ogAmH3zwgQuqCoStW7cOzKNte/smf/nLX6x///4uTHvVxVdffdW9HgU7b542bdq4/fDcf//9se6HjoGCkt4XotcRWxOygnO3bt3sqaeeco9VTVuzZo0NGzbMBTWPXotCtI67aH/nzZtnL7/8siUWwhgAINGp6er111+3bdu2uaYyVU+OHj3qfmbKlCkwX3BI8D5sVWUKHnNNFZ9gy5cvtx9++ME1h6myE3xTNSj6OuOiZtIpU6bYnj17XNPUmTNnAh/08XmtCiHRB+yNPh5V8Hq94KNjovvff/+93XPPPYFAIPXq1XP7pJMH9JrUPKjKnZor/+///s9VeryLUcfnmMS1H7FRE54nV65cVrlyZVdNCg5jwdtQ1UzNzV610NuX8+fPu+qV3gc6PmvXrnWVuriOWzBVxlRBVcjS69d7JfiYeRTQNm/eHBLivWOqimMwNesGr0PHQccjMRHGAACJ6r///a9ralRVZNCgQa46ow9mfZAqSASHsegjmKvC5PVB8wTPL/pAV4XI62MWTNsKl5rznn76aRs7dqyrsCjYqAlt6dKl8Xi1/wse0fcxJjGFhuAmxOiVN2+dek50/FQNU0VKgeKJJ55wlUZVmeJzTOLaj9hEvw6jHnv7FdPvUsdEgUgVL/UZi2l9UVFRbrvhHDvP+++/75o81QSsUKbXrWbFWrVqhcynpsyYfi96HH2/E3I8bhR9xgAAiUr9xNRkpr4/agJTM6U6y4dDnbhV0YheLQumTvJqTlRzZvBJA7p5w094H7Bxfaiqv5L6q3kdwrW8mkSDhbMer/nrRj7A1YSn1xRMzYuqEKr/nEcVKTW/qV+WKk6q/CiQhHNMbkTw70CBWicbBO9XdApmOiliy5Yt1+xP0aJF3evy5vE6+4dDyzVp0sQ166piqH5zwU2cHl2ySOtWxTH6MS1btqz5jTAGAEhU+rBV+NIHnygkqK9QONT8pCqaKiCi5sPonal1xpw6bqvz/9mzZwPzqcqlZkYpUqSI+7l9+/Y491OduY8dO+aC1PTp010n8mDhrEf9unTh6IEDB7qqiyoy6lQefFbg9egsSAUX9bvSvuj4vfjii676pWArqjJ+++23ruqjeXRxeVXzVGUK55jcCPXz0okK2rbO6FQFM7iJMiYKijrDUcdVVTItq+bUQYMGBeZRVU+d5tVHS69JTdrRmxY9eh3q/+UF5t9++8397rx+fdHpWGjd3vtwyZIlrr9bbGfEJiXCGAAgUakvjwKKqhZqInvooYfCHhld/XdU8dEHpoJGjRo1XPUquKlM1TOd/agP2Zw5c7o+TOoLpAqVN76ZOtzrg1t9nRSoYvoAVhBQINO8Wo+a/KLvZzjrUYVIzWbqFK6mRoUnhTpVscKlkwfUbPraa6+5161qoio7akL1qAO9QpFer/ZXgVVNlupjFc4xuRHesBbaN21XYVPbuF7AVLOv+oTpuGh+7f+D//9EAFFY1BmY7dq1c+vWEBT58+ePcX16HaoAavgMNWUrhKnqp5NDYqLfm4at0O9O69b7SCcMJIdR+tNcTeyG0GRISVtvTp326nV2BIDkQs1Mqqzowzd6H6qUQBUKVU2Cz0wUVU/0f1fhRE1bqtwo0Hid81XZ0QdqTH2GVF3S/26FOfUJ6tKli1tX9M7dv/76q6u6xPa/XfugztjqXB+9L1rwOrQ9hQVVX9RpPnogCGc9ojMgNU/wa1K1TJ3Wg4+P9vmXX35xzYrp04d259a8Ci/RTwjwaP9UZYrt7M7Yjkl898OrJunsRv3+9LwqUXr9wb+HmNYbnZbLnj37NX3PgvdD76G8efOGnLzhnfwRvfql46xwpubI682rfdfQIVp39P5hek/pd6/nPKos6n0QfEZwuH+v4eYNwhhhDEAyk9LDWKTpGoOqaOisNp0JqIqGKmQ66xFJK3oYg0UkjHEkAQDJmiosv//9712FR5UODdqpEdKBmwVhDACQrKnTum5qilKVIbjZCklL43jpagJUxSKLMAYASBGu10EciU993+LqC4aE4WxKAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAASUaXNpo2bdpNt62YaNuTJk1K9P273naQ/BHGAAARs3XrVncRcF2zUdch1MjkwRYtWmTffPNNkhzxhGxr/Pjx9o9//CMi29e2Fy9eHOvzkToW19sOkj/fxxnbvXu3u6TFoUOH3EVUn3zyyTgv/6FvANHfdBrzZNiwYUmwtwDgo5dyJvH2TsVr9smTJ1uvXr3cpYv0/1yDgz788MPuAtxdu3Z183Tv3j3k+oHJjS6urWtn6kLViS25HwukkjC2efNmq127tjVq1Mhq1arlvpF88MEHLuXH9gZdsWKFbdmyxV193RPXRVoBAEnjxRdftD//+c82dOjQwLQRI0a4L92eH3/80Q0cWqdOHfd47Nix7iLc+j/+9ddfu+v8abT96tWr2xdffGHz5893F8Du3Lmz3XbbbSHrrVKligt7Hn2x1wWm9aU+JtOnT7cFCxa4+7rguK5v+dhjjwVG9FdFbO3ata4g4IXHIUOGWLFixezAgQPu82nPnj1WsmRJa9++vbtWZrDvvvvOPvroI8uSJYv7XLue2I5F4cKF3X7qYuTNmzd3+xnf7cS1v1p+woQJrojhXQBd1/x87bXXrH///lauXLnr7jtuomZK/dHqTfbxxx9bv3793DeSDRs22Pvvvx/nciVKlLBOnToFbrpoLADAPwpBJ06ccJcrCpYuXTorVapUrE1zX375pT399NM2fPhwV5Hav3+/u+RO27Zt3QXC9f/+hx9+cJ8VuhyS57PPPrP169eHbGvp0qVxNtdpXfrir5tG89dnkEKeRxd6zps3rwst3nxZs2Z1AU3Bb/v27Va2bFlXSKhYsaJt27YtsOy8efPsnnvucdfPzJ49u6t6zZkzJ85jFtOxUBVxwIABLpzq4tIqWKgIEZ/tXG9/FXR133vtly5dsjZt2rjfH0EslVXGlPiV/N95553ANP0BPPDAA/b555/H+s1G9CZSKVx/9HXr1rXGjRsn0V4DAGKSNm1aa9mypQ0ePNhVY/S/XGHGq7zERQFo4cKF7nqHV69edZUcBQiFLa1XYU3VIn02qMKTUF7A8mhdCmiDBg1ygVHBR9tWKPQqY1K/fv1ASPJ069bNLadigqig0KdPHxcgRU21wZW8cKkKuGTJksC1H3/++WebOnWqC6jhbkf7Htf+6piqala1alVXjVMVTSFYQQ+pLIzpj/XixYvujR9Mj+Pq0Kg3kZL77bff7t48Kmc3bdo0zg6XUVFR7uY5ffp0hF4FACC4mVCBZsaMGe7svrNnz7omOHVBqVChQqwH6r777guEDzUZKhipz5n+34ue02fDvn37buhgK+ipiqRK05EjR1w1T11i9AU/uHoXTJ8zarFR4OnRo4dbh25a5uDBg24erUvVO71+T8GCBd3riq969eqFXIRbn3c7duwIezvh7K8ohOr3otaly5cvuyZhNd0ilYWxc+fOuZ8qAQdT2dV7LiYvvfRSyDcttafrG4NKrLFVyNS3ILgPQ1Io+ULc5enkbnfmNpZixbPTMYDIUIhQANBNQefbb791zWjqSqL+UbH53e9+F/JYISymaQoNN0ItLqrAqXP+HXfcYRkzZrQPP/zQzpw5E+syR48edT81f6FChQLT1WyqPlty+PDhQIUvmCpsCkLxEf11q5nXe93hbCec/fWouVPhV1VHr98aUlkYy5Ejh/sZ3AdAjh8/fk2fg2DRS941a9Z0nStXrlwZaxhTqVYl2+DKmJYBACQOhSd1Ixk4cKDr/6Uz5qN3eL8R6viu7i7B1OdJX+hjoiqdzsbXSQKq3on6Xan5LpjXmd+joCLly5e3Fi1axLhu7/NEZ48GV9jUxFi8ePEEvb6Ebiec/RWFZTVx6nekoPzCCy/YmDFjIravSCEd+PWmUujSmZHB1E+gUqVK8VrXb7/9Fue3D/3RKvwF3wAAkaXms+gBadmyZa6z/C233BLRbZUuXdqWL18eEkjCGWsruAqmVpPo1FSnooBH+/3ggw+6Mw+DiwcKfjrpTPSZon5l6n/lfRbpjEW99kgKZzvh7K/89a9/dU2e6jum29tvv+3OXEUqC2P61tSqVSs3Lo2+nciqVavcTU2OwX/c+mblnfGhM2iCjRs3zpVl6cQPAP5SKFDFRmFATYE6o0/DSegMea//V6SotUMd/NXU9vjjj1vDhg2v6YMcTF1iNBCtPndat27tlps5c+Y1XWUeeeQRd6KA5lNHeFWh9DmkL/Xqv6Vqk16fzkgMDm0a4FZVt2rVqrnhMnRTv7dIC2c719vfNWvWuCE71K9P1Uo1Uaoypv5jx44di/g+4/rSXI1vg3YE6Y3RoEEDl9j1ZtIZJDrV9s033wzMoz8GNUGqYqZ2c/0RKc2rM6hOAvBGe9bZNuFSM6WqcjptOLGqZPQZ8xF9xpDCaawtjVyvoRbiGgQ7OdKHuYZWUGdz9VlSn97gflCqXqnTvNdHSdUY/R8OPstRncnVD+quu+4KTJs1a5ZrilNfKI+2oRO+tH4tv2nTJtf8pk7wMW1LVq9e7YZ8UAi5//77XVhUMAs+I1HPK+ipUKB+yaqW6aNSy6ozvV6XushEbxLVZ5lGCVDfLO2POs2riODtT3ThHAttU+tVoIrPduLaX/WbU4Xw0UcfDcyv5VUhU8jTWZaIzN9ruHnD1zDmvQGU8r0R+KM3UaoMrT+44LHEfvrpJ3e2iMZh0R9mfM8AIYxdHx34Af+k5DAGpDbnIxDG0ieHs280Hk1s9I0lOpXBYzsNGQAAICXhQuEAAAA+IowBAAD4iDAGAADgI8IYAACAjwhjAJBM+XyyO4AwaDiVG+X72ZQAgFAae0qX5dGwPhpRPfolegAkjy9LuuKE/k41qLGudZpQhDEASGZ0ceiiRYvavn37bPfu3X7vDoA4aPBdDUh8I1eZIIwBQDKULVs2K1OmjF28eNHvXQEQxxcnjZd6o9VrwhgAJON/9LoBuLnRgR8AAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAUnsYO3PmjO3evdsuXrwYr+XOnTtn27Zts8OHDyfavgEAANy0YezSpUvWvXt3y5cvn9WoUcMKFChg//znP8NevkuXLla+fHl79dVXE3U/AQAAbsow9tprr9mnn35qmzZtctWt119/3Tp27GgbN2687rKTJk2yffv2WaVKlZJkXwEAAG66MPbuu+9a165drWzZsoFKV+nSpW3ChAlxLrd161YbNGiQffDBB5YuXbok2lsAAICbKIwdPHjQVbZq1aoVMv3uu++2tWvXxrrc+fPnrXXr1jZy5EgrWbJkEuwpAADATRjGjh075n6qv1gwPT569Gisy/Xt29cqVKhg7dq1C3tbUVFRdvr06ZAbAABAqg5jXvPihQsXrglO6dOnj3GZBQsW2IcffmjPPvusO4tSN81/4sQJdz82I0aMsJw5cwZuxYoVi/CrAQAASJiYU08SKFKkiKVJk8YOHDhwTfNl0aJFY1xGnfx1xmWnTp0C0zQkxqFDh2zVqlW2ZcuWGPuQDRgwwFXUPKqMEcgAAECqroxlz57dqlevbvPmzQtMU5Xsq6++snr16oWEs127drn7bdu2DVTEvNvtt99uHTp0cPdj68yfKVMmy5EjR8gNAAAgxYaxcuXKuWEpole14mvo0KFuXLFRo0bZ8uXLrU2bNpY5c2br0aNHYJ6BAwda06ZNb2g7AAAAN1UYU+d5DT+hpr4mTZq4scLiO3q+NG7c2GbOnGlffvml9ezZ0wWxZcuWWZ48eQLzFCpUyG677bZY13Hrrbe6pksAAICUKM3Vq1evJmRBLbZkyRKbPHmy/ec//3HNjgppGiusYsWKlpypz5g68p86dSrRmixLvjDHUrLdmdtYivXSKb/3AAAACzdvJLjPmDrf33///W7gVfXr6tevn7311ltuRHyNHTZ9+nQX2AAAAJBIZ1MqbC1evNhVxz755BMrUaKEq4zprEddc1LPjR8//kY2AQAAcFNLUBjbs2ePTZ061d3Uib958+Y2d+7ckLMg1QesatWqhDEAAIBIhzF1mldz5HPPPef6ieXOnfuaecqUKWP33ntvQlYPAACQaiSoz9iwYcPs+++/t969e18TxGbMmBG4HzyGGAAAACIUxl588cVYn3v88ccTskoAAIBUKaIj8Kvjvk7hBAAAQCL0GWvQoEGM9+XKlSu2Y8eOkE78AAAAiGAY+/3vf+9+Lly4MHDfkyFDBnftyNatW8dnlQAAAKlavMKYrkcp+fLlsz/96U+JtU8AAACpRoL6jBHEAAAAkrgy5vUF0/Uor9cvTPMAAAAggmGsTp06Md4HAABAEoQxDfQa030AAAAkk3HGAAAAkMh9xsJBnzEAAIBE7DMGAAAAH/uMAQAAIDLoMwYAAOAjxhkDAADwEeOMAQAA+IhxxgAAAFLKhcKjO3r0qO3YscPdL1u2rLuAOAAAABK5A//p06etU6dOVrBgQatdu7a76X6XLl3s119/TcgqAQAAUqUEhbGePXvad999Z7Nnz3bVMd10f9WqVdajR4/I7yUAAMBNKkHNlDNnzrTVq1dbxYoVA9MeeughK1asmNWsWTOS+wcAAHBTS1BlLHfu3FagQIFrpufPn989BwAAgEQMY82bN7cBAwbY+fPnA9N0/8UXX7RmzZolZJUAAACpUtjNlE2aNAkJXgsXLrRPPvnEKlWqZFevXrUtW7bY8ePHrX79+om1rwAAAKk3jJUuXTrksUJYsGrVqiVoB7Zt22YTJkywQ4cOWeXKle2pp56ybNmyxbnM3Llz7csvv3RnbqrfWufOnWkeBQAAKVKaqypr+WTdunVWt25da9mypdWqVcsmTpxoly5dspUrV1qmTJliXEbDZ2hoDVXg0qdPb1OnTrX9+/e7sztvueWWsLar5XPmzGmnTp2yHDlyWGIo+cIcS8l2Z25jKdZLp/zeAwAALNy84WsYe/DBB12gmjPnf8Hl2LFj7ozMN954I9YhMg4fPuxOFPCcPXvWcuXKZZMmTbIOHTqEtV3C2PURxgAAuDHh5o0Ej8B/5MgRW7Bgge3Zs8dVs4INHDjwustHRUXZokWL7L333gtMy5s3r6t4qRkytjAWHMRk9+7dbvslS5ZM6EsBAADwTYLC2LfffmsPP/ywG3V/+/btVrVqVXdZpHPnzrm+Y+GEMS/EFS9ePGS6Hi9dujTOZTds2GDDhg1ziXPz5s2uqfLee++NM/jp5tFyAAAAKXZoi/79+9vgwYNd53svHKnfls64vP/++8NahxeOsmTJEjJdnfeDh8yISaFChax169ZuiI0SJUrYmDFj3FUAYjNixAhXJvRuagoFAABIsWFs06ZN1rVr1/+tIG1aF6w02Otbb71l//rXv8Jah0KRnDhxImS6+o2pD1hcNOBsixYtXFPm4sWLXRAbPXp0rPNrTDS113q3vXv3hrWPAAAAyTKMaUgJryOagtHPP//s7mta9HAVm6JFi7oAt3HjxpDpelylSpWw9yVz5syuv5iaPWOjMzO1b8E3AACAFBvGgqnD/XPPPeeuV/nHP/4x7PHG0qRJY23btnVnQZ48edJNU18xDVHRrl27wHzjx4+3vn37uvuqwE2fPj1kPWvWrHG3uPqMAQAA3FRhbOzYsYH7I0eOtCtXrlinTp1c81/w2ZHXM3z4cHd2ZPny5e2BBx6wxo0bu0sq6b5HQWv+/Pnufrp06eyLL76wUqVKuXnvueceu++++6xnz55u/DEAAICUxtdxxkSbVzVMI/BrVP9bb7015Pm1a9e6yyw1bNgwMO3AgQOuOTNr1qwuyGlIjPhgnLHrY5wxAACS+Thjoo7zGtJCypYta/ny5Yv3OtRcWaNGjVifr169eoxnU+oGAACQKpsplfTULKlxxmrXru1uuq+mQnXuBwAAQCKGMfXRUtPi7NmzXXVMN91ftWpVrCPnAwAAIELNlDpzcvXq1VaxYsXAtIceesgNplqzZs2ErBIAACBVSlBlTOODaXyx6HRmpJ4DAABAIoYxXYZIo9oHX7ZI9zUsRbNmzRKySgAAgFQp7GZKXXcyOHgtXLjQPvnkEzcchYan2LJlixuCQoPAAgAAIMJhrHTp0iGPFcKChTvyPgAAABIQxsaMGRPurAAAAEiqa1MCAADAhzC2YMECa9CggRUpUsQKFy7s7msaAAAAEjmMTZs2zV2oW6Pu66zKv/zlL+6+pr3//vsJWSUAAECqlKBBX1999VWbNGmSdejQITDtqaeechfzHj58eMh0AAAARLgytmvXLnv00Uevma4xxvQcAAAAErEypn5iS5cutaZNm4ZMX7x4sRUtWjQhqwRuaiVfmGMp2e7MbSzFeumUJTe8H3yUDN8PQILC2LPPPmtt27a1Xr16WY0aNdw0XST8nXfesVdeeYWjCgAAkJhhrE+fPpYnTx4bOXKkvfHGG25a2bJl7e2336a/GAAAQGKHsWXLlrnQpdvly5fdtHTp0iVkVQAAAKlagjrw16tXL3BfIYwgBgAAkIRhrESJEvbjjz8mcJMAAAC4oTA2ZMgQ69Spk61Zs8Z+++03u3TpUsgNAAAAidhnrGPHju7nXXfdFePzV69eTchqAQAAUp0EhTGNJwYAAAAfwthnn31ms2bNctWvRx55JMaR+AEAAJAIfcYmT57swtfKlSvdIK+6/NGUKVPiswoAAAAkNIyNGTPGxo4da1u2bHE3PR49enR8VgEAAICEhrGdO3da586dA4+ffPJJ++mnn+KzCgAAACQ0jJ07d86yZs0aeJwtWzY3tAUAAACSqAP/sGHDrjtt4MCBCdwdAACA1CVeYaxIkSI2fvz4606Lbxg7duyYHTlyxEqWLGmZM2cOa5m9e/da+vTprVChQvHaFgAAQIoNY/v27Yvoxi9evGhdu3a1jz76yPLnz2+nTp1yJwV06dIl1mXefPNNGzVqlLt//vx5y5Ejh7377rtWv379iO4bAABAsr0cUqQMHz7c5s+fb9u3b3eVLlXYunXrZuvXr49x/suXL9uuXbtsxYoVtmfPHjt48KA1b97cDbFx+PDhJN9/AACAFB3GJkyY4Cpjap6UNm3aWLly5WzixIkxzp8uXTpXOVPTqKRNm9b69u1rZ86csbVr1ybpvgMAAKToMHbgwAH75ZdfrEaNGiHTa9WqZevWrQt7PRs3bnQ/S5QoEfF9BAAASJbXpowEddqXvHnzhkzX46NHj4a1jpMnT9rTTz/tLstUoUKFWOeLiopyN8/p06cTvN8AAAA3RWUsQ4YM7mdwSPIee8/F5ezZs9akSRM37tm0adPinHfEiBGWM2fOwK1YsWI3uPcAAAApPIyp31eaNGlcc2UwNV1eLywpiDVu3Nj1Ffvqq68sV65ccc4/YMAAd6amd9PJAgAAAKk6jGn0fvUXmzt3bmCahqpYuHBhyDAVCk4629KjEf8ffvhh10SpeaM3c8YkU6ZMbgiM4BsAAECq7jMmL7/8sgtW5cuXt7vvvtudKZk9e3br3r17YJ6hQ4faypUrbfPmzXbp0iVr2rSpbd261f71r3+5cc+8sc+KFy9uefLk8fHVAAAApLAw1qhRI1cZGzt2rM2aNcsqV65sy5Ytc/26PApZqoLJr7/+6jr+a9T9fv36haxryJAhbrwxAACAlMTXMCYNGzZ0t9gMHjw4cF99wzZs2JBEewYAAHCTD/oKAACQ2hHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAgNQcxjZu3GhPPfWUtWjRwoYOHWonT5687jI7duywP/3pT9akSRPbtGlTkuwnAADATRfGVq9ebTVr1rQrV65Y06ZN7csvv7Q6derYuXPnYl1mxIgRbt506dLZnDlz7NixY0m6zwAAADdNGBswYIA1atTIxo0bZx07drS5c+fa7t27bdKkSbEuo/m2bdtmvXv3TtJ9BQAAuKnC2Pnz523p0qXWvHnzwLRcuXJZ/fr1bd68ebEuV7hwYUuTJk0S7SUAAEDiSm8+2bNnj12+fNmKFi0aMr1YsWK2ePHiiG4rKirK3TynT5+O6PoBAABSXGXswoUL7meWLFlCpuux91ykqJ9Zzpw5AzcFPgAAgFQdxtQkKcePHw+Zrg753nOR7Jt26tSpwG3v3r0RXT8AAECKC2NFihSxvHnz2vfffx8yfcOGDVa1atWIbitTpkyWI0eOkBsAAECqDmPqhN++fXubOHFiYHiKBQsW2Lp166xDhw6B+caOHevGIQMAALgZ+daBX1555RVXGStXrpyVLVvW1q9fby+//LLde++9gXn0/MqVKwOPFy5caKNHj3ZnY3pNkKqwtWnTxt0AAABSEl/DWLZs2WzRokVuFP5Dhw5ZxYoV3dAVwZ555pmQSpmCW48ePdz9Pn36BKYrzAEAAKQ0voYxT5UqVcJ+TkNhRB8OAwAAIKXy/dqUAAAAqRlhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEAAAAfEcYAAAB8RBgDAADwEWEMAADAR4QxAAAAHxHGAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAFJ7GNu7d6+tWbPGTp8+najLAAAAJDe+hrHz58/bY489ZuXKlbP27dtbwYIFbezYsRFfBgAAILlK7+fGhw4daqtXr7affvrJChUqZDNnzrRmzZpZjRo1rGbNmhFbBgAAILnytTI2ZcoU69q1qwtV8uijj1qlSpXc9EguAwAAkFz5Vhn75Zdf7NChQ1a9evWQ6apwrV+/PmLLSFRUlLt5Tp065X4mZn+zK1G/WUp2Os1VS7GSYT9C3g8+4v0Q+UPK/wcgLF7OuHr1avIMY8ePH3c/8+bNGzJdj73nIrGMjBgxwjVvRlesWLEE7XtqkNNSsNdS9N4nSyn6iPJ+iDjeD0D8nDlzxnLmzJn8wliGDBkCHfKDnTt3zjJmzBixZWTAgAHWt2/fwOMrV6648KYQlyZNmht6HTdrkldQ1RmrOXLk8Ht34DPeD+D9AP4/JIwqYgpihQsXjnM+38KYPuzTpk1r+/fvD5mux8WLF4/YMpIpUyZ3C5YrV64b2v/UQEGMMAbeD+D/A/i8SLi4KmK+d+DPkiWL3XPPPTZr1qzAtLNnz9pXX31lDRs2DEzbuXNnoD9YuMsAAACkFL6eTTls2DA3NIWaERWwdGZk/vz57Y9//GNgntdee82NJxafZQAAAFIKX8PYfffdZ4sXL7aff/7Z3nzzTatYsaItW7bMsmXLFpinTJkyVq1atXgtgxujJt0hQ4Zc07SL1In3A3g/gP8PiSvN1eudbwkAAICb+9qUAAAAqRVhDAAAwEeEMQAAAB8RxhDi5MmTtmbNGtu3bx9HBvbbb7/Zt99+a3v27OFowL0PNm3a5N4XwOHDh93QU8eOHeNgEMYQKa+//rq7AHuHDh2sbNmy1rp1a7tw4QIHOBXSNWD79OljpUuXtgYNGth7773n9y7BRzNmzLDy5ctb3bp1rU2bNlagQAH761//yu8kldqwYYPVq1fP7rzzTnvyySfdgOytWrVyV8NBwlAZg7N06VJ7/vnn7bPPPrMffvjBtm3b5oYQ0ThvSH127dplJUuWtC1btliJEiX83h34TEMJ6X+DfqoypnD24osv2ty5c/3eNfhg9+7d7su7rn6zbt062759uy1cuNBGjRrF7yOBGNoCTseOHW3Hjh22YsWKwBHp37+/ffzxx+4PD6nX7bffbi1atHADLgPBY0A+/vjj9uqrr3JQYDVr1nRjgo4bN46jkQC+XZsSyYva/e+9996QaTVq1HDffk6cOGG5c+f2bd8AJL9m7L1797pmbKROly9fdl/ez58/7y5JqH7GU6ZM8Xu3UizCGJzjx49b3rx5Q46G91jPEcYAyJUrV6xLly5WvHhx168UqZP6h73wwgv266+/2o8//mjPPPOMq5YiYQhjcDJkyOC+4QTzOmNmzJiRowTAdMGWbt262dq1a23JkiWWJUsWjkoqpUsQ6lKEor6EalnRZ8jo0aP93rUUiQ78cNRJW50xg+mxQlrBggU5SkAqpyDWvXt3mz17ti1atMj1JQS8z4+WLVvaF198wQFJIMIYnIYNG9r8+fMtKioqcER09pROX1YgA5C6g1iPHj3c/wQFsQoVKvi9S/DR2bNnr5m2c+fOa7q6IHw0U8Lp1auXG0uqefPm7p+umiDUKVNDXiD10fhyq1evDjRXq7O2miRy5sxplStX9nv3kMT69etnkydPtrfeesud0OM1T2m8MfoJpT46u1oneN11112uD+Hnn3/uKqazZs3ye9dSLIa2QEizpAZy1NhSGvy1d+/e7nRlpD5HjhyxZs2aXTO9SpUq9s477/iyT/BP27ZtXb+g6Bo3buzGG0PqoiswjB8/3r7++msXxjRIuJqwCeYJRxgDAADwEX3GAAAAfEQYAwAA8BFhDAAAwEeEMQAAAB8RxgAAAHxEGAMAAPARYQwAAMBHhDEACMOuXbvcSOMAEGmEMQA3NY0c/9FHH7nrK94IXRpMV6UAgEgjjAG4qX3zzTf2xBNP2OXLl/3eFQCIEWEMQKqydetWmz9/vl28eNHWrFljc+bMsQMHDlwzn665pyA3b948O3jwYKzr++9//2ufffaZLV++3M6fPx+Yfvz4cVeRi76s5l27dm2EXxWAlCy93zsAAElJYWjcuHF2yy23WM6cOe3ChQu2bt06N71BgwaBCyE/9NBDtn37dqtevbp9//33Vr58+WvCWo8ePezTTz+1mjVrukB34sQJmzlzprugeu7cuW3q1Kn27rvv2sKFCy1t2rQ2efJke+aZZ2z9+vX80gEEEMYApDp79uxxIUmBS3r27GlDhgwJhLExY8a4ebZs2WL58uVz1a0777zTMmXKFFjHm2++6aphO3fudKFOBgwYYF26dHEVtzRp0rgwVrlyZRs5cqQ99thjLohp3WXKlPHplQNIjghjAFKdW2+9NRDEpF69evbxxx8HHv/73/+2zp07uyAmBQsWtA4dOtj06dMD80yZMsVVwBYsWOBODtAtV65crgny9OnTliNHDrecqmEKYtOmTbNGjRpZ165dk/jVAkjuCGMAUp08efKEPFbFK7i/l6piJUuWvCbABdu9e7dlyZLFZsyYETK9VatWdu7cORfGpGnTpq6q9t1337n+ZwAQHWEMAKLJmzev6/8VLPpjha3GjRvb4MGD4zx+qrJt3LjRhbtXXnnFJk6cyPEGEIKzKQEgmjp16riO+B41QaqjfjA1c6rpUVWwYPv37w/c37t3r3Xv3t31GVMF7YMPPrhmPQBAZQwAohk4cKBVq1bNWrRo4fp5zZ49243AnzVr1sA8w4cPt7p161qNGjVcp/2MGTPaihUr7OjRo645Umdbqp9ZrVq17Omnn3Yd+l9++WXr1q2bm1aoUCGOOwCHyhiAm5qaB9WPS0NLSIUKFezBBx8Mmado0aIueHlKlSrl+ngVL17cDXvxhz/8wTU3qv+Xp0CBAm6Iil69etnmzZvdMBhNmjSxuXPnuudXrlzpApc6+iuISf/+/a1du3ZcVglAiDRXb/QaIQAAAEgwKmMAAAA+IowBAAD4iDAGAADgI8IYAACAjwhjAAAAPiKMAQAA+IgwBgAA4CPCGAAAgI8IYwAAAD4ijAEAAPiIMAYAAOAjwhgAAID55/8BOgrLz5btrEgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mixture = result.partial_state_dicts()[0][\"result_state\"].state_distribution()\n",
    "observed = sum(component.probability * np.abs(component.state)**2 for component in mixture)\n",
    "expected = np.zeros(2**n_index_qubits)\n",
    "expected[:len(prob_dist)] = loaded_probabilities\n",
    "np.testing.assert_allclose(observed, expected, atol=1e-10)\n",
    "\n",
    "indices = np.arange(len(expected))\n",
    "fig, axis = plt.subplots(figsize=(6, 3), layout=\"constrained\")\n",
    "axis.bar(indices - 0.18, expected, width=0.36, label=\"Target at chosen precision\")\n",
    "axis.bar(indices + 0.18, observed, width=0.36, label=\"Simulated index\")\n",
    "axis.set(xlabel=\"Index\", ylabel=\"Probability\", xticks=indices, ylim=(0, 0.6))\n",
    "axis.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prep-08",
   "metadata": {},
   "source": [
    "## Use alias preparation inside an algorithm\n",
    "\n",
    "Pass `True` as the final argument to apply UNPREPARE, using the **same registers**. PREPARE followed immediately by UNPREPARE restores the initial state.\n",
    "\n",
    "After an intervening SELECT, the workspace need not return to zero. Retain it for subsequent coherent operations and include it in any preparation-state reflection. The {doc}`block-encoding demo </examples/block_encoding/block_encoding_demo>` shows the complete PREPARE \u2192 SELECT \u2192 UNPREPARE composition and how to change the alias-QROM fanout."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}