notebooks/textbook/CHSH_Inequality.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CHSH Inequality\n",
    "\n",
    "This tutorial shows how to run the CHSH (Clauser, Horne, Shimony and Holt) inequality experiment in Braket on local simulator and a QPU.\n",
    "\n",
    "For easy understanding, variables used here follow closely the reference: Scientific Background on the Nobel Prize in Physics 2022 [1]\n",
    "\n",
    "\n",
    "# Background \n",
    "\n",
    "The CHSH inequality is a generalization of Bell's inequality [2]. For a singlet state\n",
    "$$|\\psi ^{-}\\rangle = \\frac{1}{\\sqrt{2}}\\left(|+-\\rangle + |-+\\rangle \\right )$$\n",
    "the CHSH inequality is\n",
    "\n",
    "$$\n",
    "|E(A_1,B_1) + E(A_1,B_2) + E(A_2,B_1) - E(A_2,B_2)| \\leq 2\n",
    "$$\n",
    "for measurement settings $A_1, A_2$ on the first particle and settings $B_1,B_2$ on the second particle. This reduces to Bell's original inequality if $A_1=B_1$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References \n",
    "\n",
    "[1] The Nobel Committee for Physics, Scientific Background on the Nobel Prize in Physics 2022, https://www.nobelprize.org/uploads/2022/10/advanced-physicsprize2022.pdf\n",
    "\n",
    "[2] John F. Clauser, Michael A. Horne, Abner Shimony, and Richard A. Holt. Proposed Experiment to Test Local Hidden-Variable Theories. Phys. Rev. Lett. 23, 880 – Published 13 October 1969; Erratum Phys. Rev. Lett. 24, 549 (1970) https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.23.880"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run on a local simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from braket.devices import LocalSimulator\n",
    "from braket.tracking import Tracker\n",
    "\n",
    "from braket.experimental.algorithms.chsh_inequality import (\n",
    "    create_chsh_inequality_circuits,\n",
    "    run_chsh_inequality,\n",
    "    get_chsh_results,\n",
    ")\n",
    "\n",
    "tracker = Tracker().start()  # to keep track of Braket costs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CHSH Inequality experiment consists of four circuits acting on two qubits each. The four circuits are grouped together in the `create_chsh_inequality_circuits` function. The input arguments are the two qubits to act on and the two separate angles to rotate each qubits by. The default values for the angles are:\n",
    "$$A_1 = \\pi/2$$\n",
    "$$B_1 = \\pi/4$$\n",
    "$$A_2 = 0$$\n",
    "$$B_2 = 3\\pi/4$$\n",
    "which gives maximum violation of the CHSH inequality. \n",
    "\n",
    "Below, we print the four CHSH circuits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "circuit_a1b1: measurement setting for qubit0 at angle a1, qubit1 at angle b1\n",
      " T  : |0|1|2|   3    |Result Types|\n",
      "                                  \n",
      "q0 : -X-H-C-Rx(1.57)-Probability--\n",
      "          |          |            \n",
      "q1 : -X---X-Rx(0.79)-Probability--\n",
      "\n",
      "T  : |0|1|2|   3    |Result Types|\n",
      "\n",
      "circuit_a1b2: measurement setting for qubit0 at angle a1, qubit1 at angle b2\n",
      " T  : |0|1|2|   3    |Result Types|\n",
      "                                  \n",
      "q0 : -X-H-C-Rx(1.57)-Probability--\n",
      "          |          |            \n",
      "q1 : -X---X-Rx(2.36)-Probability--\n",
      "\n",
      "T  : |0|1|2|   3    |Result Types|\n",
      "\n",
      "circuit_a2b1: measurement setting for qubit0 at angle a2, qubit1 at angle b1\n",
      " T  : |0|1|2|   3    |Result Types|\n",
      "                                  \n",
      "q0 : -X-H-C----------Probability--\n",
      "          |          |            \n",
      "q1 : -X---X-Rx(0.79)-Probability--\n",
      "\n",
      "T  : |0|1|2|   3    |Result Types|\n",
      "\n",
      "circuit_a2b2: measurement setting for qubit0 at angle a2, qubit1 at angle b2\n",
      " T  : |0|1|2|   3    |Result Types|\n",
      "                                  \n",
      "q0 : -X-H-C----------Probability--\n",
      "          |          |            \n",
      "q1 : -X---X-Rx(2.36)-Probability--\n",
      "\n",
      "T  : |0|1|2|   3    |Result Types|\n"
     ]
    }
   ],
   "source": [
    "circuits = create_chsh_inequality_circuits(a1= np.pi / 2, a2= 0, b1= np.pi / 4,  b2=3 * np.pi / 4)\n",
    "\n",
    "print(\"\\ncircuit_a1b1: measurement setting for qubit0 at angle a1, qubit1 at angle b1\\n\",circuits[0])\n",
    "print(\"\\ncircuit_a1b2: measurement setting for qubit0 at angle a1, qubit1 at angle b2\\n\",circuits[1])\n",
    "print(\"\\ncircuit_a2b1: measurement setting for qubit0 at angle a2, qubit1 at angle b1\\n\",circuits[2])\n",
    "print(\"\\ncircuit_a2b2: measurement setting for qubit0 at angle a2, qubit1 at angle b2\\n\",circuits[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circuits can be run on the Braket local simulator with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "local_tasks = run_chsh_inequality(circuits, LocalSimulator(), shots=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results of the inequality experiment are called using the `get_chsh_results` function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_a1b1 = -0.7071067811865475, E_a1b2 = -0.7071067811865472, E_a2b1 = -0.7071067811865474, E_a2b2 = 0.7071067811865472\n",
      "\n",
      "CHSH inequality: 2.828427124746189 ≤ 2\n",
      "CHSH inequality is violated!\n",
      "Notice that the quantity may not be exactly as predicted by Quantum theory. This is may be due to finite shots or the effects of noise on the QPU.\n"
     ]
    }
   ],
   "source": [
    "chsh_value, results, E_a1b1, E_a1b2, E_a2b1, E_a2b2 = get_chsh_results(local_tasks, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the CHSH value is very close to $2\\sqrt{2}\\approx 2.828$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we plot the CHSH value for various values of the angles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "For all the iterations:\n",
      " Max CHSH_inequality: 2.8281587012240337 Corresponding theta: 2.3482611754105527\n"
     ]
    }
   ],
   "source": [
    "local_simulator = LocalSimulator()\n",
    "angles = np.linspace(0, 2 * np.pi, 100)\n",
    "\n",
    "chsh_inequality_lhs_max = 0\n",
    "chsh_inequality_lhs_max_theta = 0\n",
    "chsh_values = []\n",
    "\n",
    "for theta in angles:\n",
    "    circuits = create_chsh_inequality_circuits(a2=0, b1=theta, a1=2 * theta, b2=3 * theta )\n",
    "    local_tasks = run_chsh_inequality(circuits, local_simulator, shots=0)\n",
    "    chsh_value, results, E_a1b1, E_a1b2, E_a2b1, E_a2b2 = get_chsh_results(local_tasks, verbose=False)\n",
    "    if chsh_value > chsh_inequality_lhs_max:\n",
    "        chsh_inequality_lhs_max = np.abs(chsh_value)\n",
    "        chsh_inequality_lhs_max_theta = theta\n",
    "    chsh_values.append(chsh_value)\n",
    "\n",
    "print(\n",
    "    \"\\nFor all the iterations:\\n Max CHSH_inequality:\",\n",
    "    chsh_inequality_lhs_max,\n",
    "    \"Corresponding theta:\",\n",
    "    chsh_inequality_lhs_max_theta,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the CHSH Value against theta angle to determine which theta gives maximum violation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAHBCAYAAACFRxauAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/F0lEQVR4nO3df3RU9Z3/8ddNCCHBZAQn4AjEBBVYBNRKREQRj4JWWrFxNat2lbb63WWDxWqPBfwF1Rp119bWWtxWlFoXG7RBKlXXyM+lSKUoAgqI/BCVAKaWCSQhDJn7/YPJNJOZm0ySmblzZ56PcziHuZnMvOedz/183vOZz/2MYZqmKQAAAADKsDsAAAAAIFlQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAENDD7gBSgd/v1759+5SXlyfDMOwOBwAAAG2YpqnDhw/rtNNOU0aG9fwwxXEM7Nu3T4MGDbI7DAAAAHTgs88+08CBAy1/TnEcA3l5eZJOJDs/Pz/uz+fz+fTWW29p0qRJysrKivvzOQV5sUZuIiMv1shNZOTFGrmJjLxYS3Ru6urqNGjQoGDdZoXiOAZallLk5+cnrDjOzc1Vfn4+J1or5MUauYmMvFgjN5GRF2vkJjLyYs2u3HS0BJYL8gAAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACEiL4vjFF1+UYRhd+gcAAID0kRZfAnLo0CGZpml3GAAAAEhyKT9zvGPHDg0dOtTuMAAAAOAAKV8cr1ixQhMmTLA7DAAAADiAo4rjiooKlZSUKC8vT/369dO1116r7du3t/s7x44dC/u+7q48DoDUU+M9qh1eQzXeo3aHAgc51CSt2/WVaryNdocCIA4cVRyvWrVK5eXlWrdunaqrq+Xz+TRp0iTV19dHvP/nn3+uwsLCbj8OgNRTuX6vJjyxWr/8KFMTnlityvV7JUk13kat3VlL4YOIbeHlDZ9rznuZ+tfn/6pxjy4PthsAqcNRF+S9+eabIbcXLFigfv36acOGDRo/fry2bdumjIwMDRkyRJL01ltv6cYbb+z043SkqalJTU1Nwdt1dXWSJJ/PJ5/P1+nX1Vktz5GI53IS8mIt3XNT4z2qT//WoNNPyZXH1Us13qOaVbVZ/sB1un5TmlW1WX87clT/9dYO+U0pw5AenjJc158/0N7gbZLubeblDZ/rviUfhbSFi890674lH8nUiZ2MWtrN2OI+wXbVup2lm3RvM1bIi7VE5yba53FUcdyW1+uVJPXt21fHjh3Tq6++qq+++kqPP/64JKmhoUE5OTmdepxoVFRUaO7cuWHH33rrLeXm5kYbfrdVV1cn7LmchLxYS8fcvHPAUOWuDJkyZMhU2WC/3L0kv5kZcj+/Kf3n/34cUvjc++qH8u3dpJOz7Yg8OaRjmznUJM15LzOsLdx6lj9iu1n0+grVHlVYOxvbPz13SUrHNhMN8mItUblpaGiI6n6G6dA9zvx+v6655hodOnRIa9askSSZpqlLL71Uq1ev1pdffqkVK1bohhtu6PTjdCTSzPGgQYNUW1ur/Pz8rr+oKPl8PlVXV2vixIlh66nTGXmxlq65qfEe1YQnVgdniKUTs4CL/t8Y3fDrv4QcNyRF6gxf/O5ojSnum3azgunWZlr/fT/9W4P+9fm/ht3n52Wj9INFm6JqTxmGtPLu8WnRVlqkW5uJFnmxlujc1NXVye12y+v1tluvOXbmuLy8XFu2bAkpaA3D0KhRo7Rx40Zt3rxZU6ZM6dLjdCQ7O1vZ2eFTSVlZWQlt+Il+PqcgL9ZSPTc13kbtrq1Xsbu3PK4cfe71hhQs0omZPp/fUEXpyODSigxD+tFVw/TYm9tC7p9pGDqjf76qNtaE3LeidKTKSsKvZ0hFqd5mpBPrz9u2hQxDYW3hgsFuPTxluO599UOZMpRpGHqkdISO+RWxnX3hPaasrB4hbTIdpEOb6QryYi1RuYn2ORxZHE+fPl1Lly7V6tWrNXBg6HrA6667Tn/4wx/Uv3//Dmdx23scAM7StsCpKB2p8UMKIhY5Re5cjT3jFI0t7qNFr6/QDVdfpkJ3nk7OzdLsqi1qNs1g4SMpbH3y7KotGj+kIG2KnVRW420M+/s+/uZ2/ejrw/T4G9tD2oLHlaPrzx8o395NOuPcC3VG/3x5XDmq8TZGbGebvjikm59dl5ZvqgAnc9RuFaZpavr06Vq8eLGWL1+u4uLisPuMHz9eb731llwuV7ceB4BzRCpwZldtkXSiIMkMfBV86yJHkjyuXjrLZQY/+i4rKdSamZfppdsv1JqZl6mspFC7a+vDZgWbTVN7aqNbu4bkZvX3HTXg5LC20OLkbGlMcd9W7SgnrJ3dc9VQPfbGtrA2yS4oQPJz1MxxeXm5Fi5cqCVLligvL0/79++XJLlcruCFd5mZmTrvvPN01VVXdetxADhHewVsWUmhxg8p0J7aBhW5czuc7fW4ckLuU+zubTn7LIUv5UBya/v3au/v27YttKdtO2uvTdJOgOTmqOJ43rx5khT2jXfPP/+8pk6dGrz91FNPtbuuJNrHAeAMHRWwnSly2mqZFWy73MLjyom4lIOPzZOX1d/L6u/bWW3bWXttEkDyclRxHO3GGh0tuHboBh0ALLRXwMZCpNlnq6UcrEVOTu39vTr76UI04t0mAcSPo4pjALASjwKntbazgnxs7iwd/b268+mClXi3SQDxQXEMwJEirfWNR4FjpaOlHEgudv29rNoka9WB5OWo3SoAQDqxdnTco8t102/+onGPLlfl+r0JjyHSDgUtH5vXeBu1dmctOxPYqO3foL2/V6IlQ/sFYI2ZYwCOkkxrfSN9bM5Fevaz+hskwzKHZGq/ACJj5hiAoyTbvsMeV47GnnFKuxfpMYOcOB39DVr/veyQbO0XQDiKYwCO0rJ2tLVkWetL4WO/ZP8bJHP7BXACxTEAR0mmtaNtUfjYL9n/BsncfgGcwJpjAI6TDGtHI2FvW/s54W+QrO0XwAkUxwAcKZHbtnVGe4UP23fFR9u8OqH4TNb2C4DiGABiLlLhwy4W8WGVV4pPAF3FmmMASS0V9gxmF4v4SMW8pkJ7B5yOmWMASStVZlv5qun4SLW8pkp7B5yOmWMASSmVZgWTfQcFp0qlvKZSewecjuIYQFJK9v1qO4Ptu+IjlfKaSu0dcDqWVQBISi2zgq0LBqfOCkrsYhELkfLkhJ0popFq7R1wMopjAEnJCfvVdha7WHRde3lKhZ0pUrG9A05FcQwgaaXKrKAVq3Wm44cUpNxr7Y50yVOqt3fAKSiOASS1VJgVtJJquy3ESzrlKZXbO+AUXJAHADZJpd0W4ok8AUgkimMAsEkq7bYQT+QJQCKxrAIAbMQ60+iQJwCJQnEMICmk83ZmVutM0zUnVq87Xdfjpms7AOxCcQzAdmxnFi5dc5Kur9sK+QASjzXHAGzF1+aGS9ecpOvrtkI+AHtQHAOwFV+bGy5dc5Kur9sK+QDsQXEMwFZs0xUuXXOSrq/bCvkA7EFxDMBWbNMVLl1zkq6v2wr5AOzBBXkAbMc2XeGscpJKOxdEei20hVDkA0g8imMASSFdt+lqT9ucpNLOBe29FtpCKPIBJBbLKgDAAVJp54JUei0AUg/FMQA4QCrtXJBKrwVA6qE4BgAHSKWdC1LptQBIPRTHAOAAHe1cUONt1NqdtUm7NKF1fOzCACCZcUEeADiE1c4FyX6hnlV87MIAIBlRHAOAg7TducDq4rbxQwqSouDsKL5kiBEAWmNZBYCES/YlAE6S7Be3JXt8TsZ5BMQHM8cAEirZlwA4TcvFba0L0NYXtyX6S0PaPl9H8aFrOI+A+GHmGEDCsL9t7LV3cVvl+r0a9+hy3fSbv2jco8tVuX5vXGOJ9HxcfBd7nEdAfDFzDCBh2vuInWKp6yJd3JbotcjtPR8X38UW5xEQX8wcS6qoqFBJSYny8vLUr18/XXvttdq+fbvdYQEph/1t48fjytHYM04JFkftFVCxWKva9jE6WlvcNj50HecREF8Ux5JWrVql8vJyrVu3TtXV1fL5fJo0aZLq6+vtDg1IKXzEnjhWBdSmLw5ZLrWwKprbHo+0fIKCLXE4j4D4YlmFpDfffDPk9oIFC9SvXz9t2LBB48ePtykqIDXxEXtitBRQs6u2qNk0lWkYuueqoXrsjW0Rlz6s/vjLkAu8Hp4yXL0lvbzhc9235KPg8R9dNUyPvRn+GGtmXhb2fBRs8cN5BMQPxXEEXq9XktS3b9+IP29qalJTU1Pwdl1dnSTJ5/PJ5/PFPb6W50jEczkJebGWbLlx5/aQuzBfkr0xJVteYq30XI/GFvfR3q8aVNg3V5/+rSHi0od3d9WGrRe+b8lHuvNs6cl1H4Ucf/SNbWrzEGo2Te08UBf2fB5Xr5TLbTK1mWQ5j1okU26SCXmxlujcRPs8hmmabfu5tOb3+3XNNdfo0KFDWrNmTcT7zJkzR3Pnzg07vnDhQuXm8hEigOR0qEma816mTP1j/YMhU7ee5deCHZlh97/29Ga9+mn4cUNm2GPM+VqzTs6OT9wAEAsNDQ266aab5PV6lZ+fb3k/iuM2pk2bpjfeeENr1qzRwIEDI94n0szxoEGDVFtb226yY8Xn86m6uloTJ05UVlZW3J/PKciLNXITWTrmpe0yiYenDNfFZ7o14YnVIbPKGYZ059nH9eSHPcKO/3DSWfqvt3aEPMb150fuL1NNOraZaJGbyMiLtUTnpq6uTm63u8PimGUVrUyfPl1Lly7V6tWrLQtjScrOzlZ2dvgUSVZWVkIbfqKfzynIizVyE1k65eWmC4t12T+dGrZWte164Yem/JN6H9ikh6cM1/1LtoasIy4rKdS3vjYorde7plOb6SxyExl5sZao3ET7HBTHkkzT1B133KHFixdr5cqVKi4utjskAIgbjysnrKBte4GXO7eHXn99k64/f2DEYjrSYwBAKqA4llReXq6FCxdqyZIlysvL0/79+yVJLpdLOTl0/gDSQ+uCt/WFKxTCANIJ+xxLmjdvnrxeryZMmCCPxxP8V1lZaXdoAAAASCBmjnViWQUAAADAzDEAAAAQQHEMIG6svo4YQPxw3gHdw7IKAHFRuX5vyNcRV5SOVFlJod1hASmN8w7oPmaOAcRcjbcx7OuIZ1dtYSYLiCPOOyA2KI4BxNzu2vqQb1WTpGbT1J7aBnsCAtIA5x0QGxTHAGKu2N1bGUbosUzDUJE7156AgDTAeQfEBsUxgJjzuHJUUTpSmcaJkbrla4f5IgkgfjjvgNjggjwAcdH264gZoIH447wDuo/iGEDc8LXDQOJx3gHdw7IKAAAAIIDiGAAAAAigOAYAAAACKI4BAACAAIpjAAAAIIDiGAAAAAigOAYAAAACKI4BAACAAIpjAAAAIIDiGAAAAAigOAYAAAACKI4BdFuNt1Frd9aqxttodygALHCeAtHpYXcAAJytcv1ezaraLL8pZRhSRelIlZUU2h0WgFY4T4HoMXMMoMtqvI3BAVeS/KY0u2oLM1NAEuE8BTqH4hhAl+2urQ8OuC2aTVN7ahvsCQhAGM5ToHMojgF0WbG7tzKM0GOZhqEid649AQEIw3kKdA7FMYAu87hyVFE6UpnGiZE30zD0SOkIeVw5NkcGoAXnKdA5XJAHoFvKSgo1fkiB9tQ2qMidy4ALJCHOUyB6FMcAus3jymGwBZIc5ykQHZZVAAAAAAEUxwAAAEAAxTEAAAAQQHEMAAAABFAcAwAAAAEUxwAAAEAAxTEAAAAQQHEMAAAABFAcAwAAAAEUxwAAAEAAxTEAAAAQQHEMAAAABFAcA4hajbdRa3fWqsbbaHcoAGKAcxoI18PuAAA4Q+X6vZpVtVl+U8owpIrSkSorKbQ7LABdxDkNRMbMMYAO1Xgbg4OoJPlNaXbVFmabAIfinAasURwD6NDu2vrgINqi2TS1p7bBnoAAdAvnNGCN4hhAh4rdvZVhhB7LNAwVuXPtCQhAt3BOA9YojgF0yOPKUUXpSGUaJ0bTTMPQI6Uj5HHl2BwZgK7gnAascUEegKiUlRRq/JAC7altUJE7l0EUcDjOaSAyimMAUfO4chhAgRTCOQ2EY1kFAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAOIqMbbqLU7a1XjbbQ7FAAJxvmPdNbD7gAAJJ/K9Xs1q2qz/KaUYUgVpSNVVlJod1gAEoDzH+mOmWMAIWq8jcGBUZL8pjS7agszSEAa4PwHKI4BtLG7tj44MLZoNk3tqW2wJyAACcP5D1AcA2ij2N1bGUbosUzDUJE7156AACQM5z9AcQygDY8rRxWlI5VpnBghMw1Dj5SOkMeVY3NkAOKN8x/ggjwAEZSVFGr8kALtqW1QkTuXgRFII5z/SHcUxwAi8rhyGBSBNMX5j3TGsgoAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BgAAAAK6VRwfPXo0VnEAAAAAtut0cez3+/XQQw9pwIABOumkk7Rr1y5J0v3336/58+fHPEAAAAAgUTpdHD/88MNasGCBHn/8cfXs2TN4fMSIEXr22WdjGhwAAACQSJ0ujl944QX9+te/1s0336zMzMzg8XPOOUfbtm2LaXAAAABAInW6OP7iiy905plnhh33+/3y+XwxCQoAAACwQ6eL4+HDh+v//u//wo6/8sorOu+882ISFAAAAGCHHp39hQceeEC33nqrvvjiC/n9flVVVWn79u164YUXtHTp0njECAAAACREp2eOp0yZotdee01vv/22evfurQceeEBbt27Va6+9pokTJ8YjRgAAACAhOj1zLEmXXHKJqqurYx0LABvUeBu1u7Zexe7e8rhy7A4HQBKjv0A66FJxDCA1VK7fq1lVm+U3pQxDqigdqbKSQrvDApCE6C+QLjq9rCIjI0OZmZmW/wA4Q423MTjQSZLflGZXbVGNt9HewAAkHfoLpJNOzxwvXrw45LbP59P777+v3/72t5o7d27MAgMQX7tr64MDXYtm09Se2gY+LgUQgv4C6aTTxfGUKVPCjv3zP/+zzj77bFVWVup73/teTAIDEF/F7t7KMBQy4GUahorcufYFBSAp0V8gnXR6WYWVCy+8UMuWLYvVwwGIM48rRxWlI5VpGJJODHSPlI5gFghAGPoLpJOYXJDX2NioX/ziFxowYEAsHg5AgpSVFGr8kALtqW1QkTuXgQ6AJfoLpItOF8d9+vSREXjnKEmmaerw4cPKzc3Viy++GNPgAMSfx5XDIAcgKvQXSAedLo5/9rOfhRTHGRkZKigo0JgxY9SnT5+YBgcAAAAkUqeL46lTp8YhDAAAAMB+URXHmzZtivoBR40a1eVgAAAAADtFVRyfe+65MgxDpmm2ez/DMNTc3ByTwAAAAIBEi6o43r17d7zjAAAAAGwXVXF8+umnxzsOAAAAwHZd3uf4o48+0t69e3Xs2LGQ49dcc023gwIAAADs0OnieNeuXfrWt76lzZs3h6xDbtnejTXHAAAAcKpOf330jBkzVFxcrIMHDyo3N1cffvihVq9erdGjR2vlypVxCBEAAABIjE7PHL/zzjtavny53G63MjIylJGRoYsvvlgVFRX6/ve/r/fffz8ecQIAAABx1+mZ4+bmZuXl5UmS3G639u3bJ+nERXvbt2+PbXQAAABAAnV65njEiBH64IMPVFxcrDFjxujxxx9Xz5499etf/1qDBw+OR4wAAABAQnS6OL7vvvtUX18vSfrxj3+sb3zjG7rkkkt0yimnqLKyMuYBAgAAAIkSdXE8evRo3XbbbbrpppuUn58vSTrzzDO1bds2ffXVV+rTp09wxwoAAADAiaJec3zOOefonnvukcfj0S233BKyM0Xfvn0pjAEAAOB4URfH8+fP1/79+/X0009r7969uvzyy3XmmWfqkUce0RdffBHPGAHEwKEmad2ur1TjbbQ7FAAppMbbqLU7a1XjPWp3KEBMdGq3itzcXE2dOlUrV67Uxx9/rH/5l3/Rf//3f6uoqEiTJ09WVVVVvOIE0A0vb/hcc97L1L8+/1eNe3S5KtfvtTskACmgcv1ejXt0uW76zV804YnVeucAnyLD+Tq9lVuLM844Qw8//LD27Nmjl156SevWrdP1118fy9gAxECNt1H3LflIpk4MWn5Tml21hRlkAN1S423UrKrN8p/4olz5TalyVwYzyHC8LhfHkrRy5UpNnTpVU6dOVXNzs26//fZYxQUgRnbX1gcHrxbNpqk9tQ32BAQgJUTqW0wZ2vsVfQucrdNbuX3++edasGCBFixYoF27dumSSy7Rr371K11//fXKycmJR4wAuqHY3VsZhkIGsUzDUJE7176gADhepL7FkKnCvvQtcLaoZ44XLVqkq666SsXFxZo3b55uuOEGffzxx1q1apVuueUWRxfGTz/9tIqKitSrVy+NGTNG7777rt0hATHjceXo4SnDZejECJZpGHqkdIQ8LueeswDs53HlqKJ0pDIDu1VlGFLZYL88rl42RwZ0T9Qzx9/+9rc1efJkLV68WFdffbUyMrq1IiNpVFZW6q677tIzzzyjMWPG6Mknn9SVV16p7du3q1+/fnaHB8TE9ecPlG/vJp1x7oU6o38+hTGAmCgrKdT4IQXaU9ugAa6eev/Py+0OCei2qIvjzz//PCWLxZ/+9Ke6/fbb9Z3vfEeS9Mwzz+hPf/qTnnvuOc2cOdPm6IDYOTlbGlPcV1lZWXaHAiCFeFw58rhy5PP59L7dwQAxEHVxnIqF8bFjx7RhwwbNmjUreCwjI0NXXHGF3nnnHcvfa2pqUlNTU/B2XV2dJMnn88nn88Uv4ICW50jEczkJebFGbiIjL9bITWTkxRq5iYy8WEt0bqJ9HsM0TbPju6Wmffv2acCAAVq7dq3Gjh0bPH7PPfdo1apV+stf/hLx9+bMmaO5c+eGHV+4cKFyc7kQAQAAINk0NDTopptuktfrVX5+vuX9Or1bBaRZs2bprrvuCt6uq6vToEGDNGnSpHaTHSs+n0/V1dWaOHEiH5G3Ql6skZvIyIs1chMZebFGbiIjL9YSnZuWT/o7ktbFsdvtVmZmpg4cOBBy/MCBAzr11FMtfy87O1vZ2dlhx7OyshLa8BP9fE5BXqyRm8jIizVyExl5sUZuIiMv1hKVm2ifIzW2nOiinj176vzzz9eyZcuCx/x+v5YtWxayzAIAAADpIeqZ4z59+sgwOv7O9K+++qpbASXaXXfdpVtvvVWjR4/WBRdcoCeffFL19fXB3SsAAACQPqIujp988sng/03T1LRp0/TjH//Y8btYlJWV6csvv9QDDzyg/fv369xzz9Wbb76p/v372x0aAAAAEizq4vjWW28NuX3HHXfouuuu0+DBg2MeVKJNnz5d06dPtzsMAAAA2Cyt1xwDAAAArVEcAwAAAAEUxwAAAEBA1GuOW3/phXTiq5d/8pOfyOVyhRz/6U9/GpvIAAAAgASLujh+//33Q25fdNFF2rVrV8ixaLZ6AwAAAJJV1MXxihUr4hkHAAAAYLturzk+fvy4jhw5EotYAAAAAFtFXRy/9tprWrBgQcixn/zkJzrppJN08skna9KkSfr73/8e6/gAAACAhIm6OP7pT3+q+vr64O21a9fqgQce0P33369Fixbps88+00MPPRSXIAF0To23UWt31qrG22h3KADSFP0QnCrqNccffvhhyE4Ur7zyiiZOnKh7771XktSrVy/NmDGD3SoAm1Wu36tZVZvlN6UMQ6ooHanScz12hwUgjUTqh8pKCu0OC4hK1DPHhw8f1imnnBK8vWbNGl1++eXB22effbb27dsX2+gAdEqNtzE4IEmS35RmV21RjfeovYEBSBs13qMW/RAzyHCGqIvjAQMGaOvWrZKkI0eO6IMPPtBFF10U/Pnf/vY35ebmxj5CAFHbXVsfHJBaNJum9n7VYE9AANLOp39riNgP7amlH4IzRF0cX3/99brzzjv1u9/9TrfffrtOPfVUXXjhhcGf//Wvf9XQoUPjEiSA6BS7eyujzXbjmYahwr68cQWQGKefkhuxHypy0w/BGaIujh944AGVlJTo+9//vjZu3KgXX3xRmZmZwZ+/9NJL+uY3vxmXIAFEx+PKUUXpSGUGvpAn0zD0SOkIeVy9bI4MQLrwuHpZ9EM5NkcGRCfqC/JycnL0wgsvWP6cLwkBkkNZSaHGDynQntoGFblz5XHlyOfz2R0WgDQSqR8CnCLq4hiAc3hcOQxGAGxFPwSniro4Pu+882QYRof3e++997oVEAAAAGCXqIvja6+9Nvh/0zRVUVGhf//3f1ffvn3jERcAAACQcFEXxw8++GDI7SeeeEIzZszQ4MGDYx4UAAAAYIeod6sAAAAAUh3FMQAAABBAcQwAAAAERL3m+Be/+EXI7ePHj2vBggVyu90hx7///e/HJjIAAAAgwaIujn/2s5+F3D711FP1u9/9LuSYYRgUxwAAAHCsqIvj3bt3xzMOAAAAwHasOQYAAAACoi6Oly9fruHDh6uuri7sZ16vV2effbZWr14d0+AAAACARIq6OH7yySd1++23Kz8/P+xnLpdL//Zv/xa2LhkAAABwkqiL4w8++EBXXXWV5c8nTZqkDRs2xCQoAAAAwA5RF8cHDhxQVlaW5c979OihL7/8MiZBAQAAAHaIujgeMGCAtmzZYvnzTZs2yePxxCQoAAAAwA5RF8dXX3217r//fh09ejTsZ42NjXrwwQf1jW98I6bBAWhfjbdRa3fWqsbbaHcoABAV+i0ku6j3Ob7vvvtUVVWlIUOGaPr06Ro6dKgkadu2bXr66afV3Nyse++9N26BAghVuX6vZlVtlt+UMgyponSkykoK7Q4LACzRb8EJoi6O+/fvr7Vr12ratGmaNWuWTNOUdOJb8a688ko9/fTT6t+/f9wCBfAPNd7G4AAjSX5Tml21ReOHFMjjyrE3OACIgH4LThF1cSxJp59+ul5//XX9/e9/1yeffCLTNHXWWWepT58+8YoPQAS7a+uDA0yLZtPUntoGBhkASYl+C07RqeK4RZ8+fVRSUhLrWABEqdjdWxmGQgaaTMNQkTvXvqAAoB30W3AKvj4acCCPK0cVpSOVaRiSTgwwj5SOYPYFQNKi34JTdGnmGID9ykoKNX5IgfbUNqjIncsAAyDp0W/BCSiOAQfzuHIYXAA4Cv0Wkh3LKgAAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYcIAab6PW7qxVjbfR7lAAIObo45BMetgdAID2Va7fq1lVm+U3pQxDqigdqbKSQrvDAoCYoI9DsmHmGEhiNd7G4KAhSX5Tml21hdkVACmBPg7JiOIYSGK7a+uDg0aLZtPUntoGewICgBiij0MyojgGklixu7cyjNBjmYahIneuPQEBQAzRxyEZURwDSczjylFF6UhlGidGj0zD0COlI+Rx5dgcGQB0H30ckhEX5AFJrqykUOOHFGhPbYOK3LkMGgBSCn0ckg3FMeAAHlcOAwaAlEUfh2TCsgoAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BgAAAAIojgEAAIAAimMAAAAggOIYAAAACKA4BpJMjbdRa3fWqsbbaHcoAGAr+kPYoYfdAQD4h8r1ezWrarP8ppRhSBWlI1VWUmh3WACQcPSHsAszx0CSqPE2BgcCSfKb0uyqLcyYAEg79IewE8UxkCR219YHB4IWzaapPbUN9gQEADahP4SdKI6BJFHs7q0MI/RYpmGoyJ1rT0AAYBP6Q9iJ4hhIEh5XjipKRyrTODEiZBqGHikdIY8rx+bIACCx6A9hJy7IA5JIWUmhxg8p0J7aBhW5cxkIAKQt+kPYheIYSDIeVw6DAACI/hD2YFkFAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAAAABFMeATWq8jVq7s1Y13ka7QwEAR6H/RDz1sDsAIB1Vrt+rWVWb5TelDEOqKB2pspJCu8MCgKRH/4l4Y+YYSLAab2OwY5ckvynNrtrCDAgAdID+E4lAcQwk2O7a+mDH3qLZNLWntsGegADAIeg/kQgUx0CCFbt7K8MIPZZpGCpy59oTEAA4BP0nEoHiGEgwjytHFaUjlWmc6OEzDUOPlI6Qx5Vjc2QAkNzoP5EIXJAH2KCspFDjhxRoT22Dity5dOwAECX6T8QbxTFgE48rh04dALqA/hPxxLIKAAAAIIDiGAAAAAigOAYAAAACKI4BAACAAIpjAAAAIIDiGAAAAAigOAYAAAACKI4BAACAAIpjAAAAIIDiGAAAAAigOAbirMbbqLU7a1XjbbQ7FABIWfS1iJUedgdgt4qKClVVVWnbtm3KycnRRRddpMcee0xDhw61OzSkgMr1ezWrarP8ppRhSBWlI1VWUmh3WACQUuhrEUtpP3O8atUqlZeXa926daqurpbP59OkSZNUX19vd2iWarxHtcNrqMZ71O5Q0I4ab2Ows5YkvynNrtrCrAYAxBB9rXMlaz2T9jPHb775ZsjtBQsWqF+/ftqwYYPGjx9vU1TW/vHuOFO/2ro65N1xjbdRu2vrVezuLY8rx+ZIsbu2PthZt2g2Te2pbeDvAwAxQl+b/CLVJ+3VM3ZL++K4La/XK0nq27ev5X2amprU1NQUvF1XVydJ8vl88vl8cYutxns07N3xrKrNGlvcR2s+qdV9Sz4KfqT08JThuv78gXGLJRm15D6ef4POGOjKVoahkE47w5AGuHomPMZky02yIC/WyE1k5MWaXblJpr42knRvMy9v+DysPrn4TLdlPeNx9YpbLNH+DQzTNM2O75Ye/H6/rrnmGh06dEhr1qyxvN+cOXM0d+7csOMLFy5Ubm5u3OLb4TX0y48yw45PPatZv92RIVNG8JghU3O+1ixJ+vKooYJepk7OjltosPDOAUOVu078bQyZKhvs19j+nHIAEEv0tcnhUFNozXGoSZrzXmZYfXLrWX4t2BFez0wf3qyzXPH7uzU0NOimm26S1+tVfn6+5f1SsjieOXOmHnvssXbvs3XrVg0bNizk2LRp0/TGG29ozZo1GjjQetY10szxoEGDVFtb226yu6vGe1QTnlgd9u74p9eP0p2LNoXd/7aLT9dzf/40bWaTfT6fqqurNXHiRGVlZdkdTlCN96j2ftWgwr65cX1H3J5kzY3dyIs1chMZebFmd26Soa+NxO68JEqkGeJBfXL1r8//Ney+Py8bpR8s2hRWz6y8e3xc/3Z1dXVyu90dFscpuazi7rvv1tSpU9u9z+DBg0NuT58+XUuXLtXq1avbLYwlKTs7W9nZ4dOwWVlZcW34he4sVZSODLsid8wZ7vCPlKRgYSyd+Nn9S7bqsn86NeXXYMX779BZhe4sFbrz7A5DUvLlJlmQF2vkJjLyYs2u3CRTXxtJKreZGm9jsDCW/lFzVP3H2LD6JNMwdMFgd8R6Jt5/v2jzn5LFcUFBgQoKCqK6r2mauuOOO7R48WKtXLlSxcXFcY6ue8pKCjW2uI8Wvb5CN1x9WbAhVZSO1OyqLWo2TWUahr53cZF+/X+7Q36XCxQAAECsWV0U2XDMH1afPFI6Qh5XjmU9kwxSsjjujPLyci1cuFBLlixRXl6e9u/fL0lyuVzKyUnOItLj6qWzXGbIRw9lJYUaP6RAe2obVOQ+se752TW7w96ttfyMnS0AAEBXta4jit29I84QF7lzNfaMU0Lqk9Y1R6R6JhmkfXE8b948SdKECRNCjj///PMdLs1INh5XTkijs3q3xmbpAACgqyLVEVY1hxRenyS7tC+OU/B6xKC2s8keV47lZunjhxQ4quECAIDEs6oj1sy8TGtmXhZxhthp0r44TnVt362xWToAAOiq9uqIsWeckhK1RNp/fXS6aVkX1Frbtchrd9bytZsAACCsLuiojkgFzBynGY8rh7XIAACgQ1Z1QXvri1MBxXEaYi1yfLADCAAkH/rmrmmvLohUR6QSiuM0xVrk2GLWHQCSD31z13VUFzhtB4rOYM0xJKXHGqJ4sXp3zbptALAPfXP3pHNdQHEMSf9Yi5xpnDgTUnENUby09+4aAGAP+ubuSee6gGUVCLJaQ8R6rfa1981AAAB70Dd3TqSxPtXXFluhOEaItmuIWK/VsfZ2AAEA2IO+OXrtjfWpvLbYCsUxLLGDRfTS9d01ACQz+uaOMdaHoziGJXaw6Jx0fHcNAMmOvrl9jPXhuCAPltL5SlUAANIBY304imNYSucrVQEASAeM9eFYVoF2tbdei10sAABwFnal6BjFMToUab0Wu1gAAOAs7EoRHZZVoNP41iEAAJyFsTt6FMfoNL516EQns3ZnLZ0KADhUuvXjjN3RY1kFOi3dv3WIJSUA4Gzp2I+n+9jdGcwco9PS+cpWPpYCAGdL1348ncfuzmLmGF2SrrtYsFk6ADhbOvXjbcdjdqWIDsUxuiwdd7HgYykAcLZ06cetxmN2pegYyyoQM+nwURUfSwGAs6VDP54O43E8MXOMmEmXj6r4WAoAnC3V+/F0GY/jheIYMZMuH1VJbJYOAE6Xyv14Oo3H8cCyCsRMOnxUBQBAsmM87h5mjhFTqf5RFQAATsB43HUUx4g5q4+qnLjFmxNjBgB0nRP7fauYU3npSDxRHCMhnLjFmxNjBgB0nRP7fSfGnOxYc4y4c+KWMk6MGQDQdU7s950YsxNQHCPu2ttSJlk5MWYAQNc5sd93YsxOQHGMuGvZUqa1ZN9SxokxAwC6zon9vhNjdgKKY8Rde1vK1HgbtXZnbdJ9BMQ2OACQXpzQ77cdM50QsxNxQR4SItKWMsl+EQHb4ABAeknmft9qzEzmmJ2K4hgJ03pLGauLCMYPKbDlxGYbHACAFLnft3t7t47GTMaq2KI4hi2S6Xvfk30GGwBgn2QYI5JpzEwHrDmGLZLlIgK2wQEAWEmWMSJZxsx0QXEMW3R0EUGiLtRjGxwAgBU7x4jW4yAX3iUWyypgG6uLCBL5EVbLu/HWnR/vxgEAkn1jhNU4yIV3iZE2M8cvvviiDMPo0j/Ej8eVo7FnnNLhhXqxnEHm3TgAIBp2fMrZ3jjYdsxEfKTNzPGhQ4dkmmbHd4StOrrooLtXDPNuHADQGfH+lLPtuMbFd/ZLi+J4x44dGjp0qN1hIArtfYTV3Y4omq1wAABoq+0YEavtSCONa+OHFLDcz2ZpsaxixYoVmjBhgt1hIApWH2FJsvyY6R8fax0Neay2H3dx8R0AIBY6Gk8iLbeo8R4NOWZVYEtiuZ/NHDVzXFFRoaqqKm3btk05OTm66KKL9Nhjj3U4K3zs2DFlZWV1+3GQGJE+wlq7szZiR/T8mj16ds2u4LvuG4oNXS3ejQMA4qczn3I+PGW4Nh0w9IMnVoeMSYP65loW2Cz3s5ejiuNVq1apvLxcJSUlOn78uGbPnq1Jkybpo48+Uu/evSP+zueff67CwsJuP05rTU1NampqCt6uq6uTJPl8Pvl8vm68wui0PEcinssu7twechfmSzrxOge6ssM6IkMKFsbSiZ9V7srQt/bUhr0bn1W1WSvvHq+HpwzXfUs+CnZQD035J7lze6R0LqX0aDNdQV6skZvIyIu1dMqNO7dHxPHE5zseNv7c++pHkjLUMny1jEmL/t+YsHEtw5AGuHrK5/OFjYOpKNFtJtrnMUwHX6X25Zdfql+/flq1apXGjx8vSdq2bZsyMjI0ZMgQSdJzzz2nG2+8UTk51u+6Ij1Oe+bMmaO5c+eGHV+4cKFyc5mFjJd3Dhiq3JUhU4YMmZrg8WtFTWbY/a49vVmvfhp+fPrwZp3lMnWoSfryqKGCXqZOzk5E5ACAVNR2PNnhNfTLj8LHn0imD29W7VGFjGtlg/0a29+xZVnSa2ho0E033SSv16v8/HzL+zlq5rgtr9crSerbt6+kE8snXn31VX311Vd6/PHHJZ1IRHuFcaTH6cisWbN01113BW/X1dVp0KBBmjRpUrvJjhWfz6fq6mpNnDgxZLlIqrta0n94j2rvVw0q7HviTciqwMdULQyZuvGKEv3xuffC3o3fcPVl8rh6JTboJJGubaYj5MUauYmMvFgjNyfWFf9qa9txSZJMmfrH1rCtx6TW41q6jVGJbjMtn/R3xLHFsd/v15133qlx48ZpxIgTF2z17NlTP/rRj3TppZdK+seMcGcfpyPZ2dnKzg6fcszKykpoh5Do50sGhe4sFbrzgrcrSkdqdtUWNZtmYM2xX18rcoccb7mYofXvpat0bDPRIC/WyE1k5MVaOuem0J0VNv48NOWftGnTZi3anSm/qbAxqe24lo4S1WaifQ7bi+OZM2fqsccea/c+W7du1bBhw0KOlZeXa8uWLVqzZk3IccMwNGrUKG3cuFGbN2/WlClT2n1sq8eBM7S+aGGAq6fe//PysONczAAASJS24487t4d6H9ik/ygdry+8xxiTHMD24vjuu+/W1KlT273P4MGDQ25Pnz5dS5cu1erVqzVw4MCw+1933XX6wx/+oP79+7e7zKGjx4EztOw/6fP59H6E4wAAJFLr8aflIjCPq1fazxA7he3FcUFBgQoKCqK6r2mauuOOO7R48WKtXLlSxcXFEe83fvx4zZw5U9OnT+/W4wAAACC92F4cd0Z5ebkWLlyoJUuWKC8vT/v375ckuVyukIvuMjMzdd555+mqq67q1uMAAAAgvTjqG/LmzZsnr9erCRMmyOPxBP9VVlaG3fepp56ynJHuzOMAAAAgfThq5rgzWzK3d0Wig7d2BgAAQBw5auYYAAAAiCeKYwAAACCA4hgAAAAIoDgGAAAAAiiOAQAAgACKYwAAACCA4hgAAAAIoDgGAAAAAhz1JSDJquVLRerq6hLyfD6fTw0NDaqrq2v3y07SDXmxRm4iIy/WyE1k5MUauYmMvFhLdG5a6rSOvgyO4jgGDh8+LEkaNGiQzZEAAACgPYcPH5bL5bL8uWHyXcrd5vf7tW/fPuXl5ckwjLg/X11dnQYNGqTPPvtM+fn5cX8+pyAv1shNZOTFGrmJjLxYIzeRkRdric6NaZo6fPiwTjvtNGVkWK8sZuY4BjIyMjRw4MCEP29+fj4nWgTkxRq5iYy8WCM3kZEXa+QmMvJiLZG5aW/GuAUX5AEAAAABFMcAAABAAMWxA2VnZ+vBBx9Udna23aEkFfJijdxERl6skZvIyIs1chMZebGWrLnhgjwAAAAggJljAAAAIIDiGAAAAAigOAYAAAACKI4BAACAAIpjh3n66adVVFSkXr16acyYMXr33XftDgkAACCiiooKlZSUKC8vT/369dO1116r7du32x1WuyiOHaSyslJ33XWXHnzwQb333ns655xzdOWVV+rgwYN2h5ZUvvOd7+i+++4LO/7oo4/KMAzdeeediQ8qSbTkxomdVUfmzZunUaNGBb9paezYsXrjjTe69FiR2hDtJ/XMmTNHhmGE/Bs2bFiXHivV2kw8zqdU6HeKiorC2oxhGCovL+/0Y6Vam7GyatUqlZeXa926daqurpbP59OkSZNUX19vd2jWTDjGBRdcYJaXlwdvNzc3m6eddppZUVFhY1TJ5fjx46bb7Tb/8pe/hBx/9913zaKiInPUqFHmjBkz7AnOZq1zc+WVV5rPP/+8uWXLFnPjxo3m1VdfbRYWFppHjhyxO8wu++Mf/2j+6U9/Mj/++GNz+/bt5uzZs82srCxzy5YtnXqcSG0ondrPsGHDTEkR/z311FN2hxdTDz74oHn22WebNTU1wX9ffvllpx8nFdtMPM6nVOh3Dh48GNJeqqurTUnmihUrOvU4qdhmonXw4EFTkrlq1SrTNJOzz6E4doimpiYzMzPTXLx4ccjxW265xbzmmmvsCSrBojmBVq9ebXo8HtPv9wd/7/Dhw+ZZZ51lVldXm5deemlKdjpdzU2Ltp1VqujTp4/57LPPmqYZfQfcNk/p0H5a+/DDD01J5rJly8yamhpzz549ZkZGhvnyyy+bR48etTu8mHrwwQfNc845x/LntJlQsTifWkuFfmfGjBnmGWecEXx9tJmO7dixw5Rkbt682TTN5OxzWFbhELW1tWpublb//v1Djvfv31/79++3KarE+sMf/iBJWrZsmWpqarRnzx5lZGTo5Zdf1u233y5J+uMf/6hvfvObMgwj+Hvl5eWaPHmyrrjiClviToSu5qaF1+uVJPXt2zdxQcdRc3Ozfv/736u+vl5jx46VFF2OpPA8pUP7ae3AgQPq0aOHxo0bp1NPPVW1tbXy+/265JJLku5brGJhx44dOu200zR48GDdfPPN2rt3b/BntJkTYnk+teb0fufYsWN68cUX9d3vfjf4+mgz7fP7/brzzjs1btw4jRgxQlJy9jk9bHlWoAtan0DZ2dnasGFD2Am0ZMkS/exnPwv+zu9//3u99957Wr9+vV1hJ0RXctMiUmflVJs3b9bYsWN19OhRnXTSSVq8eLGGDx8uKbocSaF5Spf209rmzZs1ZMiQYE4++OAD9evXL+yNeSoYM2aMFixYoKFDh6qmpkZz587VJZdcoi1btigvLy/t20ysz6fWUqHfefXVV3Xo0CFNnTo1eCzd20xHysvLtWXLFq1ZsyZ4LBn7HIpjh3C73crMzNSBAwdCjh84cECnnnqqTVElVkcn0NatW7Vv3z5dfvnlkqTPPvtMM2bMUHV1tXr16mVb3InQ2dy0FqmzcqqhQ4dq48aN8nq9euWVV3Trrbdq1apVGj58eFQdcOs8pVP7aW3Tpk0aOXJk8PYHH3wQcjuVfP3rXw/+f9SoURozZoxOP/10LVq0SN/73vfSvs3E8nxqKxX6nfnz5+vrX/+6TjvttOCxdG8z7Zk+fbqWLl2q1atXa+DAgcHjydjnsKzCIXr27Knzzz9fy5YtCx7z+/1atmxZ8GOuVNfRCfTHP/5REydODHYwGzZs0MGDB/W1r31NPXr0UI8ePbRq1Sr94he/UI8ePdTc3Jzw1xAvnc1Ni5bOasWKFSGdlVP17NlTZ555ps4//3xVVFTonHPO0c9//nNJ0XXArfOUTu2ntU2bNmnUqFHB2x988EHI7VR28skna8iQIfrkk08k0WZieT61lgr9zqeffqq3335bt912W8jxdG8zkZimqenTp2vx4sVavny5iouLQ36ejH0OxbGD3HXXXfrNb36j3/72t9q6daumTZum+vp6fec737E7tITo6ARasmSJpkyZErx9+eWXa/Pmzdq4cWPw3+jRo3XzzTdr48aNyszMTGj88dTZ3HTUWaUKv9+vpqYmSdF1wK3zlE7tp4Xf79eHH34YkpedO3eqqKjIvqAS6MiRI9q5c6c8Ho8k2kxb3TmfpNTqd55//nn169dPkydPDjlOmwlXXl6uF198UQsXLlReXp7279+v/fv3q7GxMXn7HFsuA0SXPfXUU2ZhYaHZs2dP84ILLjDXrVtnd0gJ0dzcbObm5pqvvfZa8NjAgQPNn//856ZpmuaBAwfMrKysDrdhSsWrgLuSm2nTppkul8tcuXJlyLZEDQ0NCY8/VmbOnGmuWrXK3L17t7lp0yZz5syZpmEY5ltvvdVhjkwzujaUiu2ntY8//tiUZH766afBY1//+tfNk08+2VyzZo2NkcXH3Xffba5cudLcvXu3+ec//9m84oorTLfbbR48eDDt20w8zqdU6Xeam5vNwsJC80c/+lHY8XRuM1ZksXvH888/n7R9DsUxHKGjE+jZZ581x40b1+HjpFqnY5pdy017nZVTffe73zVPP/10s2fPnmZBQYF5+eWXm2+99ZZpmtEVfdG0oVRsP+msrKzM9Hg8Zs+ePc0BAwaYZWVl5ieffGKaJm0mHudTqvQ7//u//2tKMrdv3x5yPN3bTCoxTNM0EzRJDcTNNddco4svvlj33HOP3aEkHXITHfKEzqLNWCM3kZEXZ2DNMVLCxRdfrBtvvNHuMJISuYkOeUJn0WaskZvIyIszMHMMAAAABDBzDAAAAARQHAMAAAABFMcAAABAAMUxAAAAEEBxDAAAAARQHAMAOm3q1Km69tpr7Q4DAGKO4hgAUtA777yjzMxMTZ482e5QAMBRKI4BIAXNnz9fd9xxh1avXq19+/bZHQ4AOAbFMQCkmCNHjqiyslLTpk3T5MmTtWDBguDPVq5cKcMwtGzZMo0ePVq5ubm66KKLtH379pDHePjhh9WvXz/l5eXptttu08yZM3XuuedaPqff71dFRYWKi4uVk5Ojc845R6+88kqcXiEAxA/FMQCkmEWLFmnYsGEaOnSovv3tb+u5555T2y9Dvffee/XEE0/or3/9q3r06KHvfve7wZ/9z//8j37yk5/oscce04YNG1RYWKh58+a1+5wVFRV64YUX9Mwzz+jDDz/UD37wA33729/WqlWr4vIaASBe+PpoAEgx48aN0w033KAZM2bo+PHj8ng8evnllzVhwgStXLlSl112md5++21dfvnlkqTXX39dkydPVmNjo3r16qULL7xQo0eP1i9/+cvgY1588cU6cuSINm7cKOnEBXmHDh3Sq6++qqamJvXt21dvv/22xo4dG/yd2267TQ0NDVq4cGFCXz8AdAczxwCQQrZv3653331XN954oySpR48eKisr0/z580PuN2rUqOD/PR6PJOngwYPBx7jgggtC7t/2dmuffPKJGhoaNHHiRJ100knBfy+88IJ27twZk9cFAInSw+4AAACxM3/+fB0/flynnXZa8JhpmsrOzg6ZCc7Kygr+3zAMSSfWDXfFkSNHJEl/+tOfNGDAgJCfZWdnd+kxAcAuFMcAkCKOHz+uF154QU888YQmTZoU8rNrr71WL730koYNG9bh4wwdOlTr16/XLbfcEjy2fv16y/sPHz5c2dnZ2rt3ry699NKuvwAASAIUxwCQIpYuXaq///3v+t73vieXyxXys+uuu07z58/Xf/7nf3b4OHfccYduv/12jR49WhdddJEqKyu1adMmDR48OOL98/Ly9MMf/lA/+MEP5Pf7dfHFF8vr9erPf/6z8vPzdeutt8bk9QFAIlAcA0CKmD9/vq644oqwwlg6URw//vjj2rRpU4ePc/PNN2vXrl364Q9/qKNHj+qGG27Q1KlT9e6771r+zkMPPaSCggJVVFRo165dOvnkk/W1r31Ns2fP7tZrAoBEY7cKAECHJk6cqFNPPVW/+93v7A4FAOKKmWMAQIiGhgY988wzuvLKK5WZmamXXnpJb7/9tqqrq+0ODQDijpljAECIxsZGffOb39T777+vo0ePaujQobrvvvtUWlpqd2gAEHcUxwAAAEAAXwICAAAABFAcAwAAAAEUxwAAAEAAxTEAAAAQQHEMAAAABFAcAwAAAAEUxwAAAEAAxTEAAAAQQHEMAAAABPx/iQ4F/JKF4dEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "plt.plot(angles, chsh_values, \".\", label=\"CHSH Value\")\n",
    "\n",
    "plt.grid(which=\"major\", axis=\"both\")\n",
    "plt.yticks([-np.sqrt(2) * 2, -2, 0, 2, np.sqrt(2) * 2], [\"$-2\\sqrt{2}$\", \"$-2$\", \"$0$\", \"$2$\", \"$2\\sqrt{2}$\"])\n",
    "plt.xticks([i*np.pi/4 for i in range(9)], [\"$0$\", \"$\\pi/4$\", \"$\\pi/2$\", \"$3\\pi/4$\", \"$\\pi$\",\"$5\\pi/4$\",\"$3\\pi/2$\",\"$7\\pi/4$\",\"$2\\pi$\"])\n",
    "\n",
    "plt.xlabel(\"Angle\")\n",
    "plt.ylabel(\"CHSH Value\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the plot, we see that the maximum violation is for $\\theta = \\pi/4$ as expected. The red line shows the classical bounds of $\\text{CHSH} \\leq 2$ and the dotted black line shows the quantum bound at $\\text{CHSH} \\leq 2\\sqrt{2}$. \n",
    "\n",
    "\n",
    "Following sections how CHSH reduces to Bell's original inequality. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "For all the iterations:\n",
      "Max Bell_inequality: 1.4985175478267787 \n",
      "Corresponding theta: 1.0154642920694281 \n",
      "Corresponding E_a1b1, E_a1b2, E_a2b1, E_a2b2: 0.4440666126057739 -0.5272254676105024 -0.9999999999999998 -0.5272254676105024\n"
     ]
    }
   ],
   "source": [
    "# CHSH reduces to Bell's original inequality measurement configuration with a2=b1 and b2 middle of a1 and b1\n",
    "# thumb rule: left hand index finger within the right hand thumb and index finger. \n",
    "# left hand thumb aligned with righ hand index finger\n",
    "\n",
    "local_simulator = LocalSimulator()\n",
    "angles = np.linspace(0, 2 * np.pi , 100)\n",
    "\n",
    "bell_value_max = 0\n",
    "bell_value_max_theta = 0\n",
    "bell_value_max_E_a1b1=0\n",
    "bell_value_max_E_a1b2=0\n",
    "bell_value_max_E_a2b1=0\n",
    "bell_value_max_E_a2b2=0\n",
    "bell_values = []\n",
    "\n",
    "for theta in angles:\n",
    "    circuits = create_chsh_inequality_circuits(a2=0, b1=0, b2=theta, a1=2 * theta)\n",
    "    local_tasks = run_chsh_inequality(circuits, local_simulator, shots=0)\n",
    "    chsh_value, results, E_a1b1, E_a1b2, E_a2b1, E_a2b2 = get_chsh_results(local_tasks, verbose=False)\n",
    "    #since d=b, E_a2b1 should be -1 and abs(E_a2b1)=1, subtracting this value from both sides of the CSHS inequality gives Bell inequality\n",
    "    bell_value = np.abs(E_a1b1-E_a1b2)-E_a2b2\n",
    "    if bell_value > bell_value_max:\n",
    "        bell_value_max = np.abs(bell_value)\n",
    "        bell_value_max_theta = theta\n",
    "        bell_value_max_E_a1b1=E_a1b1\n",
    "        bell_value_max_E_a1b2=E_a1b2\n",
    "        bell_value_max_E_a2b1=E_a2b1\n",
    "        bell_value_max_E_a2b2=E_a2b2\n",
    "    bell_values.append(bell_value)\n",
    "\n",
    "print(\n",
    "    \"\\nFor all the iterations:\\nMax Bell_inequality:\",\n",
    "    bell_value_max,\n",
    "    \"\\nCorresponding theta:\",\n",
    "    bell_value_max_theta,\n",
    "    \"\\nCorresponding E_a1b1, E_a1b2, E_a2b1, E_a2b2:\",\n",
    "    bell_value_max_E_a1b1,bell_value_max_E_a1b2,bell_value_max_E_a2b1,bell_value_max_E_a2b2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHBCAYAAABg2dPCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEi0lEQVR4nO3deXQUVfr/8U9lISQYApggAQNJUEAQcAGRRZYBRtEBIo4yoCO44Bl+AgriV3B0gHE0biAujDrqCDouuIDiODpiVBARQSACAypLEJGwBCUxCyEk9fsD05K9O6nq6u56v87pc9LV3dVP36q69eTWrXsN0zRNAQAAAC4Q5nQAAAAAgL+Q/AIAAMA1SH4BAADgGiS/AAAAcA2SXwAAALgGyS8AAABcg+QXAAAArkHyCwAAANeIcDqAYFBWVqZ9+/YpNjZWhmE4HQ4AAAAqMU1TP//8s1q3bq2wsJrbd0l+vbBv3z4lJSU5HQYAAADq8P333+v000+v8XWSXy/ExsZKOlGYTZs2tf37SkpK9MEHH+i3v/2tIiMjbf8+nEC5O4Nydwbl7gzK3RmUuzP8Xe55eXlKSkry5G01Ifn1QnlXh6ZNm/ot+Y2JiVHTpk05SP2IcncG5e4Myt0ZlLszKHdnOFXudXVR5YY3AAAAuAbJLwAAAFyD5BcAAACuQfILAAAA1yD5BQAAgGsEVfK7cuVKDR8+XK1bt5ZhGHrrrbfq/Mzs2bNlGEaFR6dOnewPFgAAAAEnqJLfgoICde/eXQsWLPDpc126dFF2drbnsWrVKpsiRLDJzi3S6p05ys4t8mo5ACA4VFePZ+cepW5HcI3zO2zYMA0bNsznz0VERKhVq1Y2RIRgtnjdHs1cslllphRmSH8b2VlNJL2+fq/uenurZ3n6qK4a3bOt0+ECALxUXf2+6YChqXNXUrcjuJLf+tq+fbtat26txo0bq3fv3kpPT1fbtjXv8MXFxSouLvY8z8vLk3RisOaSkhLb4y3/Dn98l1tk5x7Vd4cL1e7UGCXGNVZ27lFPxShJZaZ019tbdWsXaf6arRWWz1yyWb1TmktShXXAGuzvzqDcnUG5W8+b+v3Pb22VFKZfFlWo28s/Q/1uPX/v795+j2Gapln32wKPYRhaunSp0tLSan3fe++9p/z8fHXs2FHZ2dmaM2eOfvjhB23ZsqXG6e9mz56tOXPmVFn+8ssvKyYmxorw4UefHzC0eFeYTBkyZGp0apniG0tPbA2v8t60dqV667uqy3+TWKaPs40K6+h9WlAeOgAQMnyp36szqXOpco6qyjqo34NTYWGhxo4dq9zc3Fpn5A355LeyI0eOqF27dpo3b55uuOGGat9TXctvUlKScnJy/Da98fLlyzV06FCmYayHk/+Dl6SBv1zmKhdmSK/d1EtX/eOLKstv7XJc8/8XUWG5IckwVOW9n9zWnxYCC7C/O4Nydwblbp3s3KNe1+8nJrs1Zcqo870n1++0CDeMv/f3vLw8xcfH15n8uqLbw8maNWumDh06aMeOHTW+JyoqSlFRUVWWR0ZG+rWy8vf3hYLK/bxu7JdSoVKTTiSxJWWG0kd11Z1LtqjUNBVuGLpn5FlqcmCT/jays+5+e5tn+Q39kvWPT7OqrOOH3GNqGx+r7NwiZeUUKCW+iRLjov34a0ML+7szKHdnUO6+q1zX7s3N9al+37Rps17LCleZKYUbhu4bdbaOlanadfyQe0yfZ/1U4XxCH+H689f+7u13uC75zc/P186dO/XHP/7R6VBgsezcoir9vJ79NEthlVptww1DyfEx6t3+VPXvkKDdOYVKjo9RfEyE/vOfTbry/NM16KxWnuWS9OyqrGrXUTnZpnIEAOtVV9f275DgU/3e5MAm/b9R/fVD7jElx8coMS5a2blF1a4jplFYlfPJnUu2qH+HBBo5QkBQDXWWn5+vzMxMZWZmSpKysrKUmZmpPXv2SJKeeOIJDR48uMJnpk+frhUrVmj37t1avXq1Lr/8coWHh2vMmDH+Dh82y8opqPofvKQb+6Uq3Dhxqav8v/3yyisxLlq9259apTI7eXliXLTSR3Wtsg5J1VaODKEDANaprmHjziVbJKnaurm2+j0xrnGFZTXV7wXHSqucT0pNU7tzCu38qfCToGr5/fLLLzVo0CDP82nTpkmSxo0bp4ULFyonJ0c7d+6s8Jm9e/dqzJgxOnz4sBISEtSvXz+tWbNGCQkJfo0d1qt8CSwlvkm1/8Ff1y9Z1/VL9rQA1Oe/9tE921ZoRUiMi9bqnTk1Vo60DACANapr2Civa6urm31V3TpqahEuvxpId7fgFlTJ78CBA1Xb/XmzZ8/W7NmzKyx79dVXbY4KTqipu0Hlfl6VWwEaorwVuFxNyXZ55QgAaLi66trKdXN9VF5HeYtwdecTursFv6BKfgGp5ktg/TskWNIK4K3aKsfyOGkZAADfVa4/a6tr7VJTizB9gYMfyS+CTm2XwMr/e/dXJVRTsk3LAADUT031p78aNk5W+XxS1/kHwSGobniDO1Wen738EtjJnOxuUPmmippaBrgRDgBqV1v9WdMNyv5U2/mn8rkKgYvkFwFt8bo96nv/Rxr7zBfqe/9HWrxuT4135wbKf921tQwAAGoW6PVnTeefld8eqnKuQuCi2wMCVqD07fUVN8IBQP0EQ/1Z+fwjSX3v/4h+wEGEll8ErLpaAALhElh1Ar1lGgACVbDUnyeffwK9tRpV0fKLgBUMLQA1CeSWaQAIZMFWfwbzucqtaPlFwAqWFoCaVNcyzQ0RAPCrmurEQL2yV51gP1e5ES2/CCiVx3YMthaA2jD8GQD8KpTqxNrOVYz5HnhIfhEwaqoI/Tlur10YGB0AfhWKdWJ156pQSvBDCd0eEBBCfWxcbogAgF+5oU4M9fNaMCP5RUAI9Yow0CbmAAAnuaFODPXzWjAj+UVACPWKkBsiAOBXbqgTQ/28Fszo8wtHVL4BoLwivHPJFpWaZkhWhKF08x4ANFSo14l1nde4Ec45JL/wu5puAAj1ilCq/oYIAHCrUK8TazqvcSOcs+j2AL+q6waAYBrbEQCAulQ+r3EjnPNIfuFX3ABQPSa/ABDKqON+xXnQeXR7gF8xDWRVXP4CEMqo4yriPOg8Wn7hV264w9cXXP4CEMqo46riPOg8Wn5hq+ruZnXDjW3equ3yl5vLBUBooI6rXk3nQUaA8A+SX9imtktdoX6Hr7e4/AUglFHH1azyeZDuIf5DtwfYgktd3uHyF4BQRh3nHc6Z/kXLL2zBpS7v0Q0EQCijjqsb50z/IvmFLbjU5Ru6gQAIZdRxteOc6V90e4AtuNQFAIB3OGf6Fy2/sASjOgAAUH+1nTMZBcJaJL9oMEZ1AACg4ao7ZzIKhPXo9oAG4Q5VAADswTnWHiS/aBDmKLdXdm6RVu/MoaIDEDSot6zDOdYedHtAg3CHqn241AUg2FBvWYtzrD1o+YXPTv6vnjtU7cGlLgDBhnrLenWdY2llrx9afuGTmv6rZ1QHazHgOYBgQ71lj5rOsbSy1x/JL7xW03/1/TskMKqDxbjUBSDYUG/Zp/I5tq7zMWpHtwd4jY73/kN3EgDBhnrLfzgfNwwtv/Aa/9X7F91JAAQb6i3/4HzcMLT8wmv8V+9/iXHR6t3+VMoYQNCg3rIf5+OGoeUXNWLKYgAAAlNN52OmQq4byS+qxZTFAAAEtsrnY0aA8A7dHlAFYzUCABBcOHd7j+QXVXAXKQAAwYVzt/dIflFlhpjyu0hPxl2kAAAErrrO3cwG9yuSX5dbvG6P+t7/kcY+84X63v+RFq/bw12kQYBKDIDTqIcCS23n7urO9W7GDW8uVtsMMYzqELi4oQGA06iHAlN1525mg6uKll8Xq6t/EGM1Bh5uaADgNOqhwFb53E1f4KpIfl3m5MtU9O0NPlRiAJxGPRRc6AtcFd0eXKS6y1Tpo7rqziVbVGqa9O0NAkxpCcBp1EPBpbwvcHXnerd2XyH5dYmaLlOtmjFIq2YMom9vkKitEgMAf6AeCj70Ba6I5NclartMRb/e4MLNiACcRj0UfCrPBldbXhDq25PkN0RVntuby1ShhSmmATiNeii41ZYXVM4hQg3JbwiqqQ8Pl6kAAIBUc/eVld8eCvl+wCS/Qa7yf2eM3QsAALxROS+QpL73f1RjP+BQaREO+eR35cqVeuihh7R+/XplZ2dr6dKlSktLczosS1TXwpvUIqbWPjxcpgIAAOVOzgtW78ypMYcIpRbhkB/nt6CgQN27d9eCBQucDsVr2blHtT3XUHbu0UrLfx2Lr6YW3iaNwhm7F0C1ahrP043jfAKoqqYxgWMahdU6sUnNdUv1+YzTQr7ld9iwYRo2bJjTYXjt19bccP1920rPf1aVW3lv7JdS7X9nhcfK6NsLoIrKdcjfRnZWE0mvr9+ru97eGhKtOQAapqZ+wAXHSn1uEa4pnwkEIZ/81kdxcbGKi4s9z/Py8iRJJSUlKikpse17s3OPVvnPauaSzWofH1Nl+TOfZlW5SzPMkNrENVKPtk3VO6W59vxYqLYtYpQY19jWuENFeRlRVv5FuVsvO/eovjtcqHannjj+q6tb7np7q27tIs1fs7VKndM7pbkkVVgHrMH+7gzK3XujzkmskkNk5x6tNueIDDO9zlvK6xY76xNvt69hmqZZ99tCg2EYXvX5nT17tubMmVNl+csvv6yYGPu6D2zPNfTE1vAqy9Paleqt76ou/01imT7ONmTKkCFTo1PL1Ps012xOANX4/IChxbvCKtQL8Y1F3QKgQayoWyZ1LtWZcfbVJYWFhRo7dqxyc3PVtGnTGt9H8luN6lp+k5KSlJOTU2thNlR27lENnLuyyn9Wr93US1f944sqyz+5rb8kVfjvDPVXUlKi5cuXa+jQoYqMjHQ6HJ9Vbu0LFsFe7oHE1zrk1i7HNf9/ERWWG5KMalp4Prmtf1DtV4EqGPf3YK1bThaM5R6IsnOPVmkR9jVvsXMfysvLU3x8fJ3Jb727PRw7dkxZWVlq3769IiJCq/dEVFSUoqKiqiyPjIy09aBpGx+p9FFdq/Sd6ZESX20fnLbxsb98Lta2mNzI7u1sh1CYnz0Yyz3Q7M3NrdIvr8yUSsqMKnXIPSPPUpMDm/S3kZ1199vbPMtv6Jesf3yaVWUdP+Qeo66xULDs76FQt5wsWMo9ULWNj6xQD5TnLZXzk/K8pfK+Y3cd4u229TlrLSws1OTJk7Vo0SJJ0rfffqvU1FRNnjxZbdq00YwZM3xdJU4yumdb9U5prtf+87GuunSQZ0dhjF7UxM3zs7udLzM59m5/aoU6JD4mQv/5zyZdef7pGnRWqwrjfD67KqvG2SBDZZxP1I26Bd6oKT+pKZ8JBD4PdTZz5kx99dVX+uSTT9S48a9N10OGDNHixYstDc4K+fn5yszMVGZmpiQpKytLmZmZ2rNnj7OB1SIxrrHOjDOrXBpIjItW7/anUumggtrmZ0foWrxuj/re/5HGPvOF+t7/kRav2+O5UzvcODFWUeXRXmqqQ05eXts6qvtOhC7qFnir5rql+nzGaT63/L711ltavHixLrzwQhnGr4PBdenSRTt37rQ0OCt8+eWXGjRokOf5tGnTJEnjxo3TwoULHYoKsE5trX0ITXbP5FjdOmgFdB/qFoQqn1t+Dx06pJYtW1ZZXlBQUCEZDhQDBw6UaZpVHiS+CBV1tfYh9NTVImfFVaLK66AV0H2oWxCqfG757dGjh959911NnjxZkjwJ77PPPqvevXtbGx0Ar9An3F2caJGjFdCdqFsQinxOfu+77z4NGzZMW7du1fHjx/Xoo49q69atWr16tVasWGFHjAC8cPL87AhtNc3CZOf2d+I7ERioWxBqfE5++/Xrp8zMTN1///3q2rWrPvjgA5133nn6/PPP1bVrVztiBADXqzzKghMtcrV9J6NAAAgW9Rqgt3379nrmmWesjgUAUI2axlp1okWuuu8MtbFgAYQ2n5PfuoYIa9uWCg8ArBLooywEenwAUJnPyW9ycnKtozqUlpY2KCAAwK9qG2UhEJLLQI8PACrzOfnduHFjheclJSXauHGj5s2bp3vvvdeywAAAgT/KQqDHBwCV+Zz8du/evcqyHj16qHXr1nrooYc0atQoSwIDAAT+KAuBHh8AVFavG96q07FjR61bt86q1QEAfhHoY60GenwAcDKfk9+8vLwKz03TVHZ2tmbPnq0zzzzTssAAAL8K9LFWAz0+ACjnc/LbrFmzKje8maappKQkvfrqq5YFBgBuFErj5YbSbwEQOnxOfj/++OMKz8PCwpSQkKAzzjhDERGW9aIAYBESkOARSuPlhtJvcQvqCriFz9nqgAED7IgDgA1IQIJHKI2XG0q/xS2oK+AmXiW/y5Yt83qFI0aMqHcwAKxDAhJcQmm83FD6LW5AXQG38Sr5TUtL82plhmEwyQUQIEhAgksojZcbSr/FDagr4DZh3ryprKzMqweJLxA4yhOQk5GABK7y8XLDf7mhOJjHyw2l3+IG1BVwG+5QA0IUkw8En1AaLzeUfkuoo66A29Qr+S0oKNCKFSu0Z88eHTt2rMJrU6ZMsSQwAA1HAhJ8Qmm83FD6LaGOugJu4nPyu3HjRl166aUqLCxUQUGBWrRooZycHMXExKhly5Ykv0CAIQEB4A3qCriFV31+TzZ16lQNHz5cP/30k6Kjo7VmzRp99913Ov/88/Xwww/bESMAhJzs3CKt3pmj7Nwip0PxOzf/dgDO87nlNzMzU08//bTCwsIUHh6u4uJipaam6sEHH9S4ceM0atQoO+IEgJDh5jFV3fzbAQQGn1t+IyMjFRZ24mMtW7bUnj17JElxcXH6/vvvrY0OAEJMTWOquqEV1M2/HUDg8Lnl99xzz9W6det05plnasCAAfrLX/6inJwcvfjiizr77LPtiBEAQoabx1R1828HEDi8bvktH8P3vvvuU2JioiTp3nvvVfPmzTVx4kQdOnRI//jHP+yJEgBChJvHVHXzbwcQOLxOftu0aaMZM2aoadOmGjRokKQT3R7ef/995eXlaf369erevbttgQJAKHDzBBBu/u0AAofX3R5uvvlmLVq0SA899JD69OmjG264QVdddZViYviPHQB84eYxVd382wEEBq9bfu+++27t2LFDGRkZSk1N1aRJk5SYmKgJEyboiy++sDNGAAg5iXHR6t3+VFcmf27+7QCc5/NoDwMHDtSiRYu0f/9+zZ07V9u2bVPv3r3VpUsXzZs3z44YAQAAAEv4nPyWO+WUU3TjjTdq1apVeuedd7R//37dfvvtVsYGAAAAWKreyW9hYaEWLlyoAQMGaMSIETr11FN17733WhkbAAAAYCmfx/ldvXq1/vnPf+r111/X8ePH9fvf/1733HOP+vfvb0d8AGySnVukrJwCpcQ3oe+lzShr71BO/kNZw828Tn4ffPBBPf/88/r222/Vo0cPPfTQQxozZoxiY2PtjA+ADZhi1n8oa+9QTv5DWcPtvO728NBDD+mSSy7RV199pS+++EI33XQTiS8QhJhi1n8oa+9QTv5DWQM+tPzu27dPkZGRdsYCwA+YYtZ/KGvvUE7+Q1kDPrT8kvgCoYEpZv2HsvYO5eQ/lDXQgNEeAAQnppj1H8raO5ST/1DWQD1GewAQ/Jhi1n8oa+9QTv5DWcPtSH4Bl0qMi+ak5yeUtXcoJ/+hrOFmPnd7GDBggF544QUVFXFnKAAAAIKLz8nvueeeq+nTp6tVq1aaMGGC1qxZY0dcAAAAgOV8Tn7nz5+vffv26fnnn9fBgwfVv39/de7cWQ8//LAOHDhgR4wAAACAJeo12kNERIRGjRqlt99+W3v37tXYsWN19913KykpSWlpafroo4+sjhMAAABosAYNdbZ27VrNmjVLc+fOVcuWLTVz5kzFx8frd7/7naZPn25VjAAQFLJzi7R6Zw6zZVmIMgVgNZ9Hezh48KBefPFFPf/889q+fbuGDx+uV155RRdffLGMX8YNHD9+vC655BI9/PDDlgcMAIFo8bo9nmljwwwpfVRXje7Z1umwghplCsAOPie/p59+utq3b6/rr79e48ePV0JCQpX3dOvWTT179rQkQAAIdNm5RZ4kTZLKTOnOJVvUv0MCw0nVE2UKwC4+J78ZGRm66KKLan1P06ZN9fHHH9c7KAAIJlk5BZ4krVypaWp3TiGJWj1RpgDs4nOf31mzZunIkSNVlufl5ek3v/mNFTEBQFBJiW+iMKPisnDDUHJ8jDMBhQDKFIBdfE5+V6xYoWPHjlVZfvToUX366aeWBAUAwSQxLlrpo7oq/Jf7HsINQ/eNOpsWygagTAHYxetuD5s2bZIkmaaprVu3av/+/Z7XSktL9f7776tNmzbWRwgAQWB0z7bq3yFBu3MKlRwfQ5JmAcoUgB28Tn7POeccGYYhwzCq7d4QHR2txx9/3NLgACCYJMZFk6BZjDIFYDWvk9+srCyZpqnU1FStXbu2wigPjRo1UsuWLRUeHm5LkAAAAIAVvE5+27VrJ0kqKyuzLRgAzsrOLVJWToFS4pvQ2gaEAI5poCqvkt9ly5Zp2LBhioyM1LJly2p974gRIywJDIB/MaEAEFo4poHqeZX8pqWlaf/+/WrZsqXS0tJqfJ9hGCotLbUqNgB+woQCQGjhmAZq5lXye3JXB7o9AKGHCQWA0MIxDdTM53F+AYQeJhQAQgvHNFAzr1p+H3vsMa9XOGXKlHoHA8AZ5RMK3Llki0pNkwkF6sBNRM5jG9SOYxqomVfJ7yOPPOLVygzDIPkFghQTCniHm4icxzbwDsc0UD2vkt+srCy74wAQAJhQoHbcROQ8toFvOKaBqujzCwBequ0mIvgH2wBAQ3k9ycXJ9u7dq2XLlmnPnj06duxYhdfmzZtnSWAAEGjKbyI6OfniJiL/YhsAaCifk9+MjAyNGDFCqamp+vrrr3X22Wdr9+7dMk1T5513nh0xAkBA4CYi57ENADSUz8nvzJkzNX36dM2ZM0exsbF688031bJlS1199dW65JJL7IgRAAIGNxE5j20AoCF8Tn63bdumV1555cSHIyJUVFSkU045RX/96181cuRITZw40fIgASCQcBOR89gGAOrL5xvemjRp4unnm5iYqJ07d3pey8nJsS4yAAAAwGI+t/xeeOGFWrVqlc466yxdeumluu2227R582YtWbJEF154oR0xAgAAAJbwOfmdN2+e8vPzJUlz5sxRfn6+Fi9erDPPPJORHgAAABDQfE5+U1NTPX83adJETz31lKUBAQAAAHZhkgsAAAC4hs8tv2FhYTIMo8bXS0tLGxQQAAAAYBefk9+lS5dWeF5SUqKNGzdq0aJFmjNnjmWBAUAgyM4tUlZOgVLimzC0VhBgewGoi8/J78iRI6ss+/3vf68uXbpo8eLFuuGGGywJDACctnjdHs1cslllphRmSOmjump0z7ZOh4UasL0AeMOyPr8XXnihMjIyrFodADgqO7fIk0hJUpkp3blki7Jzi5wNDNViewHwliXJb1FRkR577DG1adPGitUBgOOycgo8iVS5UtPU7pxCZwJCrdheALzlc7eH5s2bV7jhzTRN/fzzz4qJidG//vUvS4MDEBjc2I8yJb6JwgxVSKjCDUPJ8THOBYUauX17ufEYBerL5+T3kUceqZD8hoWFKSEhQb169VLz5s0tDQ6A89zajzIxLlrpo7rqziVbVGqaCjcM3TfqbBKLAOXm7eXWYxSoL5+T3/Hjx9sQBoBAVFM/yv4dElyRVIzu2Vb9OyRod06hkuNjXPGbg5kbt5fbj1GgPnxOfjdt2uT1e7t16+br6gEEkNr6UbrlxJoYF+2a3xoK3La9OEYB3/mc/J5zzjm1TnIhnegHbBgGE14AQc7t/SiBQMcxCvjO59EelixZopSUFP3973/Xxo0btXHjRv39739X+/bt9eabb2rXrl3KysrSrl277IgXgB+V96MM/+UfXjf1owSCAcco4DufW37vu+8+PfbYY7r00ks9y7p166akpCTdfffdWr9+vaUBAnCWG/tRAsGEYxTwjc/J7+bNm5WSklJleUpKirZu3WpJUAACi9v6UQLBhmMU8J7P3R7OOusspaen69ixY55lx44dU3p6us466yxLgwMAAACs5HPL71NPPaXhw4fr9NNP94zmsGnTJhmGoXfeecfyAAEAAACr+Jz8XnDBBdq1a5deeuklff3115Kk0aNHa+zYsWrSpInlAQIAAABW8Tn5laQmTZropptusjoWAHAM08OGLrYtgJPVK/l98cUX9fTTT2vXrl36/PPP1a5dOz3yyCNKTU3VyJEjrY4RAGzF9LChi20LoDKfb3h78sknNW3aNA0bNkw//fSTZyKL5s2ba/78+VbHBwC2qml62OzcImcDQ4OxbQFUx+fk9/HHH9czzzyjP//5z4qI+LXhuEePHtq8ebOlwQGA3WqbHhbBjW0LoDo+J79ZWVk699xzqyyPiopSQUGBJUEBgL+UTw97MqaHDQ1sWwDV8Tn5TUlJUWZmZpXl77//PuP8Agg6TA8buti2AKrj8w1v06ZN080336yjR4/KNE2tXbtWr7zyitLT0/Xss8/aESMA2IrpYUMX2xZAZT4nvzfeeKOio6N11113qbCwUGPHjlXr1q316KOP6g9/+IMdMQKA7ZgeNnSxbQGcrF5DnV199dW6+uqrVVhYqPz8fLVs2dLquAAAAADL1Sv5LRcTE6OYGG4cAAAAQHDw+Ya3AwcO6I9//KNat26tiIgIhYeHV3gAAAAAgcrnlt/x48drz549uvvuu5WYmCjDMOr+EAAAABAAfE5+V61apU8//VTnnHOODeEAAAAA9vG520NSUpJM06z7jQBCXnZukVbvzGG6WMCPOO6AhvG55Xf+/PmaMWOGnn76aSUnJ9sQEoBgsHjdHs1cslllphRmSOmjump0z7ZOhwWENI47oOF8Tn5Hjx6twsJCtW/fXjExMYqMjKzw+o8//mhZcAACU3ZukecELEllpnTnki3q3yGB8VQBm3DcAdaoV8svAHfLyinwnIDLlZqmducUBvxJODu3SFk5BUqJbxLwscJewbYvBPNxBwQSn5PfcePG2REHgCCSEt9EYYYqnIjDDUPJ8YE97jeXjFEuGPeFYD3ugEDj9Q1veXl5Xj0AhL7EuGilj+qq8F+GOgw3DN036uyAbn2q6ZIxNw25T7DuC8F43AGByOuW32bNmtU6pq9pmjIMQ6WlpZYEBiCwje7ZVv07JGh3TqGS42MC/gTMJWOUC+Z9IdiOOyAQeZ38fvzxx3bGASAIJcZFB83Jl0vGKBfs+0IwHXdAIPI6+R0wYICdcQCArcovGd+5ZItKTZNLxi7GvgC4m883vAFAsOKSMcqxLwDuRfILwFW4ZIxy7AuAO/k8vTEAAAAQrEh+AQAA4BokvwAAAHANr/r8jho1yusVLlmypN7BAAAAAHbyquU3Li7O60cgWrBggZKTk9W4cWP16tVLa9eudTokAAAAOMCrlt/nn3/e7jhss3jxYk2bNk1PPfWUevXqpfnz5+viiy/WN998o5YtWzodHgAAAPwo5Pv8zps3TxMmTNB1112nzp0766mnnlJMTIz++c9/Oh0aAAAA/Myrlt9zzz1XhmF4tcINGzY0KCArHTt2TOvXr9fMmTM9y8LCwjRkyBB9/vnnNX6uuLhYxcXFnud5eXmSpJKSEpWUlNgX8C/Kv8Mf34VfUe7OsLPcs3OP6rvDhWp3aowS4xpbvv5gxv5eMzv3G8rdGZS7M/xd7t5+j1fJb1paWkNicUxOTo5KS0t12mmnVVh+2mmn6euvv67xc+np6ZozZ06V5R988IFiYvw39/vy5cv99l34FeXuDKvL/fMDhhbvCpMpQ4ZMjU4tU+/TTEu/IxSwv1fkr/2GcncG5e4Mf5V7YWGhV+8zTNMM2bPBvn371KZNG61evVq9e/f2LP+///s/rVixQl988UW1n6uu5TcpKUk5OTlq2rSp7XGXlJRo+fLlGjp0qCIjI23/PpxAuTvDjnLPzj2qgXNXquyk2i3MkD65rT8twL9gf6/KH/sN5e4Myt0Z/i73vLw8xcfHKzc3t9Z8rV7TGx85ckRvvPGGdu7cqdtvv10tWrTQhg0bdNppp6lNmzb1Dtpq8fHxCg8P14EDByosP3DggFq1alXj56KiohQVFVVleWRkpF8PGn9/H06g3J1hZbnvzc2tkMBIUpkp/ZB7TG3jYy35jlDB/v4rf+43lLszKHdn+Kvcvf0On29427Rpkzp06KAHHnhADz/8sI4cOSLpxPi+J/etDQSNGjXS+eefr4yMDM+ysrIyZWRkVGgJBhBaUuKbKKzSbQrhhqHkeP91W0LwYb8B3MHn5HfatGkaP368tm/frsaNf70MdOmll2rlypWWBmeFadOm6ZlnntGiRYu0bds2TZw4UQUFBbruuuucDg2ATRLjopU+qqvCf7lRN9wwdN+os5UYF+1wZAhk7DeAO/jc7WHdunV6+umnqyxv06aN9u/fb0lQVho9erQOHTqkv/zlL9q/f7/OOeccvf/++1VuggMQWkb3bKv+HRK0O6dQyfExJDDwCvsNEPp8Tn6joqI8Q3+d7Ntvv1VCQoIlQVlt0qRJmjRpktNhAK6QnVukrJwCpcQ3cTxxSIyLdjwGBJ9A2W8C6VgCQonPye+IESP017/+Va+99pokyTAM7dmzR3fccYeuuOIKywMEEDwWr9ujmUs2q8w8cZd8+qiuGt2zrdNhAUGHYwmwj899fufOnav8/Hy1bNlSRUVFGjBggM444wzFxsbq3nvvtSNGAEEgO7fIc7KWTtwlf+eSLcrOLXI2MCDIcCwB9vK55TcuLk7Lly/XZ599pq+++kr5+fk677zzNGTIEDviAxAksnIKqgwTVWqa2p1TyCVbwAccS4C96jXOryT17dtXffv2tTIWAEGsfJiok0/aDBMF+I5jCbCXT90efv75Z61fv175+fmSpA0bNujaa6/VlVdeqZdeesmWAAEEB4aJAqzBsQTYy+uW35UrV+p3v/ud8vPz1bx5c73yyiv6/e9/rzZt2ig8PFxLlixRYWGhJkyYYGe8AAIYw0QB1uBYAuzjdcvvXXfdpSuvvFLff/+9br31Vo0ePVqTJk3Stm3btGXLFs2ZM0cLFiywM1YAQSAxLlq925/KyRpoII4lwB5eJ7+bNm3S7bffrjZt2uiOO+5QXl6eRo8e7Xn9D3/4g3bu3GlLkAAAAIAVvE5+8/Ly1KJFC0lSo0aNFBMTo9jYWM/rsbGxKiwstD5CAKhDdm6RVu/MYSgo2Ir9DAgNXvf5NQxDxi+d76t7DgBOYDIA+AP7GRA6vE5+TdPU4MGDFRFx4iOFhYUaPny4GjVqJEk6fvy4PRECQA1qmgygf4cE+knCMuxnQGjxOvmdNWtWhecjR46s8h6mNwbgT0wGAH9gPwNCS72TXwBwGpMBwB/Yz4DQ4tMkFwAQSJgMAP7AfgaElnpPbwwAgYDJAOAP7GdA6CD5BRD0EuOiSUZgO/YzIDTQ7QEAAACuQfILAAAA16hX8rtixQoNHz5cZ5xxhs444wyNGDFCn376qdWxAQAAAJbyOfn917/+pSFDhigmJkZTpkzRlClTFB0drcGDB+vll1+2I0YAAADAEj7f8HbvvffqwQcf1NSpUz3LpkyZonnz5umee+7R2LFjLQ0QAAAAsIrPLb+7du3S8OHDqywfMWKEsrKyLAkKAAAAsIPPyW9SUpIyMjKqLP/www+VlJRkSVAAAACAHXzu9nDbbbdpypQpyszMVJ8+fSRJn332mRYuXKhHH33U8gABoFx2bpGycgqUEt+E8VYRENgngeDjc/I7ceJEtWrVSnPnztVrr70mSTrrrLO0ePFijRw50vIAAUCSFq/bo5lLNqvMlMIMKX1UV43u2dbpsOBi7JNAcKrXDG+XX365Lr/8cqtjARDCGtJClp1b5EkyJKnMlO5cskX9OyTQ2gZHWLlP0noM+JfPyW9qaqrWrVunU089tcLyI0eO6LzzztOuXbssCw5AaGhoC1lWToEnyShXapranVNIsgBHWLVP0noM+J/PN7zt3r1bpaWlVZYXFxfrhx9+sCQoAKGjphay7Nwir9eREt9EYUbFZeGGoeT4GAsjBbxnxT5pxbEBwHdet/wuW7bM8/d///tfxcXFeZ6XlpYqIyNDycnJlgYHIPhZ0UKWGBet9FFddeeSLSo1TYUbhu4bdTatvnCMFfskVzQAZ3id/KalpUmSDMPQuHHjKrwWGRmp5ORkzZ0719LgAAS/8hayk0/y9Wm1Hd2zrfp3SNDunEIlx8eQHMBxDd0nrTo2APjG624PZWVlKisrU9u2bXXw4EHP87KyMhUXF+ubb77R7373OztjBRCEylvIwo0T14gb0mqbGBet3u1PJfFFwGjIPmnlsQHAez7f8MYsbgB8RastUD2ODcD/6jXUGQD4KjEumhM7UA2ODcC/fB7tAQAAAAhWJL8AAABwDZJfAAAAuEa9+vyWlZVpx44dnlEfTta/f39LAgMAAACs5nPyu2bNGo0dO1bfffedTLPi6NyGYVQ7+xsAAAAQCHxOfv/0pz+pR48eevfdd5WYmCjDMOr+EAD4IDu3SFk5BUqJb8Jd8AhK7MNA4PI5+d2+fbveeOMNnXHGGXbEA8DlFq/bo5lLNqvMlMIMKX1UV43u2dbpsACvsQ8Dgc3nG9569eqlHTt22BELAJfLzi3yJA3SiWlf71yyRdm5Rc4GBniJfRgIfD63/E6ePFm33Xab9u/fr65duyoyMrLC6926dbMsOADukpVT4EkaypWapnbnFHLpGEGBfRgIfD4nv1dccYUk6frrr/csMwxDpmlywxuABkmJb6IwQxWSh3DDUHJ8jHNBAT5gHwYCn8/Jb1ZWlh1xAIAS46KVPqqr7lyyRaWmqXDD0H2jzqbFDEGDfRgIfD4nv+3atbMjDgCQJI3u2Vb9OyRod06hkuNjSBoQdNiHgcBWr0kudu7cqfnz52vbtm2SpM6dO+uWW25R+/btLQ0OgDslxkWTMCCosQ8Dgcvn0R7++9//qnPnzlq7dq26deumbt266YsvvlCXLl20fPlyO2IEAAAALOFzy++MGTM0depU3X///VWW33HHHRo6dKhlwQEAAABW8rnld9u2bbrhhhuqLL/++uu1detWS4ICAAAA7OBz8puQkKDMzMwqyzMzM9WyZUsrYgIAAABs4XO3hwkTJuimm27Srl271KdPH0nSZ599pgceeEDTpk2zPEAAAADAKj4nv3fffbdiY2M1d+5czZw5U5LUunVrzZ49W1OmTLE8QAAAAMAqPie/hmFo6tSpmjp1qn7++WdJUmxsrOWBAXCH7Nyj2p5rKDv3qNrGR9b9ASCIsb8DzvO5z+/JYmNjSXwB1NvidXs0cO5KPbE1XAPnrtTidXucDgmwDfs7EBi8avk977zzlJGRoebNm+vcc8+VYRg1vnfDhg2WBQcgdGXnFmnmks0qM088LzOlO5dsUf8OCUwOgJDD/g4EDq+S35EjRyoqKsrzd23JLwB4IyunwJMIlCs1Te3OKSQZQMhhfwcCh1fJ76xZszx/z549265YALhISnwThRmqkBCEG4aS42OcCwqwCfs7EDh87vObmpqqw4cPV1l+5MgRpaamWhIUgNCXGBet9FFdFfbLhaQwQ7pv1Nm0giEksb8DgcPn0R52796t0tLSKsuLi4u1d+9eS4IC4A6je7ZV75Tmeu0/H+uqSwepbTw30CJ0sb8DgcHr5HfZsmWev//73/8qLi7O87y0tFQZGRlKSUmxNjoAIS8xrrHOjDOVGNfY6VAA27G/A87zOvlNS0uTdGKc33HjxlV4LTIyUsnJyZo7d66lwQEAAABW8jr5LSsrkySlpKRo3bp1io+Pty0oAAAAwA4+9/nNysqyIw4AAADAdj6P9jBlyhQ99thjVZY/8cQTuvXWW62ICQAAALCFz8nvm2++qb59+1ZZ3qdPH73xxhuWBAUAAADYwefk9/DhwxVGeijXtGlT5eTkWBIUAAAAYAefk98zzjhD77//fpXl7733HpNcAAAAIKD5fMPbtGnTNGnSJB06dEi/+c1vJEkZGRmaO3eu5s+fb3V8AAAAgGV8Tn6vv/56FRcX695779U999wjSUpOTtaTTz6pa6+91vIAAQAAAKv4nPxK0sSJEzVx4kQdOnRI0dHROuWUU6yOCwAAALBcvZLfcgkJCVbFAQAAANjO5+Q3JSVFhmHU+PquXbsaFBAAAABgF5+T38oTWZSUlGjjxo16//33dfvtt1sVFwAAAGA5n5PfW265pdrlCxYs0JdfftnggAAAAAC7+DzOb02GDRumN99806rVAQAAAJazLPl944031KJFC6tWBwAAAFjO524P5557boUb3kzT1P79+3Xo0CH9/e9/tzQ4AAAAwEo+J79paWkVnoeFhSkhIUEDBw5Up06drIoLAAAAsJzPye+sWbPsiAMAAACwnVfJb15entcrbNq0ab2DAQAAAOzkVfLbrFmzWie2OFlpaWmDAgIAAADs4lXy+/HHH3v+3r17t2bMmKHx48erd+/ekqTPP/9cixYtUnp6uj1RAgAAABbwKvkdMGCA5++//vWvmjdvnsaMGeNZNmLECHXt2lX/+Mc/NG7cOOujBAAAACzg8zi/n3/+uXr06FFleY8ePbR27VpLggIAAADs4HPym5SUpGeeeabK8meffVZJSUmWBAUAAADYweehzh555BFdccUVeu+999SrVy9J0tq1a7V9+3amNwYAAEBA87nl99JLL9W3336r4cOH68cff9SPP/6o4cOH69tvv9Wll15qR4wAAACAJXxu+ZVOdH247777rI4FAAAAsJXPLb+S9Omnn+qaa65Rnz599MMPP0iSXnzxRa1atcrS4AAAAAAr+Zz8vvnmm7r44osVHR2tDRs2qLi4WJKUm5tLazAAAAACms/J79/+9jc99dRTeuaZZxQZGelZ3rdvX23YsMHS4AAAAAAr+Zz8fvPNN+rfv3+V5XFxcTpy5IgVMQEAAAC28Dn5bdWqlXbs2FFl+apVq5SammpJUAAAAIAdfE5+J0yYoFtuuUVffPGFDMPQvn379NJLL2n69OmaOHGiHTECAAAAlvB5qLMZM2aorKxMgwcPVmFhofr376+oqChNnz5dkydPtiNGAAAAwBI+J7+GYejPf/6zbr/9du3YsUP5+fnq3LmzTjnlFBUVFSk6OtqOOAEAAIAGq9c4v5LUqFEjde7cWRdccIEiIyM1b948paSkWBkbAAAAYCmvk9/i4mLNnDlTPXr0UJ8+ffTWW29Jkp5//nmlpKTokUce0dSpU+2KEwAAAGgwr7s9/OUvf9HTTz+tIUOGaPXq1bryyit13XXXac2aNZo3b56uvPJKhYeH2xkrAAAA0CBeJ7+vv/66XnjhBY0YMUJbtmxRt27ddPz4cX311VcyDMPOGAEAAABLeN3tYe/evTr//PMlSWeffbaioqI0depUEl8AAAAEDa+T39LSUjVq1MjzPCIiQqeccootQQEAAAB28Lrbg2maGj9+vKKioiRJR48e1Z/+9Cc1adKkwvuWLFlibYQAAACARbxOfseNG1fh+TXXXGN5MAAAAICdvE5+n3/+eTvjAAAAAGxX70kuAAAAgGBD8gsAAADXIPkFAACAa5D8AgAAwDVIfgEAAOAaJL8AAABwDZJfAAAAuAbJLwAAAFwj5JPflStXavjw4WrdurUMw9Bbb73ldEgAAABwSMgnvwUFBerevbsWLFjgdCgAAABwmNfTGwerYcOGadiwYU6HAQAAgAAQ8i2/AAAAQLmQb/mtj+LiYhUXF3ue5+XlSZJKSkpUUlJi+/eXf4c/vgu/otydQbk7g3J3BuXuDMrdGf4ud2+/xzBN07Q5FsvNmDFDDzzwQK3v2bZtmzp16lRhmWEYWrp0qdLS0mr97OzZszVnzpwqy19++WXFxMT4HC8AAADsVVhYqLFjxyo3N1dNmzat8X1BmfweOnRIhw8frvU9qampatSoUYVl3ia/1bX8JiUlKScnp9bCtEpJSYmWL1+uoUOHKjIy0vbvwwmUuzMod2dQ7s6g3J1BuTvD3+Wel5en+Pj4OpPfoOz2kJCQoISEBNvWHxUVpaioqCrLIyMj/XrQ+Pv7cALl7gzK3RmUuzMod2dQ7s7wV7l7+x1Bmfz6Ij8/Xzt27PA8z8rKUmZmplq0aKG2bds6GBkAAAD8LeST3y+//FKDBg3yPJ82bZokady4cVq4cKFDUQEAAMAJIZ/8Dhw4UEHYrRkAAAA2YJxfAAAAuAbJLwAAAFyD5BcAAACuQfILAAAA1yD5BQAAgGuQ/AIAAMA1SH4BAADgGiS/AAAAcA2SXwAAALgGyS8AAABcg+QXAAAArkHyCwAAANcg+QUAAIBrkPwCAADANUh+AQAA4BokvwAAAHANkl8AAAC4BskvAAAAXIPkFwAAAK5B8gsAAADXIPkFAACAa5D8AgAAwDVIfgEAAOAaJL8AAABwDZJfAAAAuAbJLwAAAFyD5BcAAACuQfILAAAA1yD5BQAAgGuQ/AIAAMA1SH4BAADgGiS/AAAAcA2SXwAAALgGyS8AAABcg+QXAAAArkHyCwAAANcg+QUAAIBrkPwCAADANUh+AQAA4BokvwAAAHANkl8AAAC4BskvAAAAXIPkFwAAAK5B8gsAAADXIPkFAACAa5D8AgAAwDVIfgEAAOAaJL8AAABwDZJfAAAAuAbJLwAAAFyD5BcAAACuEVTJ78qVKzV8+HC1bt1ahmHorbfeqvMzs2fPlmEYFR6dOnWyP1gAAAAEnKBKfgsKCtS9e3ctWLDAp8916dJF2dnZnseqVatsihAAAACBLMLpAHwxbNgwDRs2zOfPRUREqFWrVjZEBAAAgGASVC2/9bV9+3a1bt1aqampuvrqq7Vnzx6nQwIAAIADgqrltz569eqlhQsXqmPHjsrOztacOXN00UUXacuWLYqNja32M8XFxSouLvY8z83NlST9+OOPKikpsT3mkpISFRYW6vDhw4qMjLT9+3AC5e4Myt0ZlLszKHdnUO7O8He5//zzz5Ik0zRrf6PpsDvuuMOUVOtj27ZtVT4nyVy6dKnP3/fTTz+ZTZs2NZ999tka3zNr1qw6Y+LBgwcPHjx48OAReI/vv/++1lzQ+CWRdMyhQ4d0+PDhWt+TmpqqRo0aVVhmGIaWLl2qtLQ0n7+zZ8+eGjJkiNLT06t9vXLLb1lZmX788UedeuqpMgzD5+/zVV5enpKSkvT999+radOmtn8fTqDcnUG5O4Nydwbl7gzK3Rn+LnfTNPXzzz+rdevWCguruWev490eEhISlJCQ4Lfvy8/P186dO/XHP/6xxvdERUUpKiqqwrJmzZrZHFlVTZs25SB1AOXuDMrdGZS7Myh3Z1DuzvBnucfFxdX5nqC64S0/P1+ZmZnKzMyUJGVlZSkzM7PCDWxPPPGEBg8e7Hk+ffp0rVixQrt379bq1at1+eWXKzw8XGPGjPF3+AAAAHCY4y2/vvjyyy81aNAgz/Np06ZJksaNG6eFCxdKknJycrRz507Pe/bu3asxY8bo8OHDSkhIUL9+/bRmzRq/tjYDAAAgMARV8jtw4MA67+CbPXu2Zs+e7Xn+6quv2hyV9aKiojRr1qwqXS9gL8rdGZS7Myh3Z1DuzqDcnRGo5e74DW8AAACAvwRVn18AAACgIUh+AQAA4BokvwAAAHANkl8AAAC4BslvgFmwYIGSk5PVuHFj9erVS2vXrnU6JAAAgGqlp6erZ8+eio2NVcuWLZWWlqZvvvnG6bBqRfIbQBYvXqxp06Zp1qxZ2rBhg7p3766LL75YBw8edDo017juuut01113SZKefPJJdevWzTMzTe/evfXee+85HGFgs7ISLN8WbAcEg/vvv1+GYejWW2+t1+fZ3703e/ZsGYZR4dGpU6d6rYtyb7gVK1bo5ptv1po1a7R8+XKVlJTot7/9rQoKCpwOrUYMdRZAevXqpZ49e+qJJ56QJJWVlSkpKUmTJ0/WjBkzHI4u9JWWlqpVq1Z69913dcEFF+idd95ReHi4zjzzTJmmqUWLFumhhx7Sxo0b1aVLF6fDDUiXXHKJ/vCHP6hnz546fvy47rzzTm3ZskVbt25VkyZNvF7PydviwIEDbAeLnHXWWfr666+rfe3xxx/XpEmT/BxRaFi3bp2uuuoqNW3aVIMGDdL8+fN9+jz7u29mz56tN954Qx9++KFnWUREhOLj431aD+Vuj0OHDqlly5ZasWKF+vfvH5j1jomAUFxcbIaHh5tLly6tsPzaa681R4wY4UxQIaRTp06mpGofjz/+uGmaprly5UozMTHRLCsrq3E9zZs3N5999ll/hR30Dh48aEoyV6xYYZqmd9vBNOveFmyH+vnf//5nSjIzMjLM7Oxsc/fu3WZYWJj5+uuvm0ePHnU6vKD0888/m2eeeaa5fPlyc8CAAeYtt9zieY393R6zZs0yu3fvXuPrlLuztm/fbkoyN2/ebJpmYNY7dHsIEDk5OSotLdVpp51WYflpp52m/fv3OxRV6HjzzTclSRkZGcrOztbu3bsVFham119/XRMmTJAkLVu2TMOHD5dhGFU+X1paqldffVUFBQXq3bu3X2MPZrm5uZKkFi1aSPJuO0g1bwu2Q8McOHBAERER6tu3r1q1aqWcnByVlZXpoosuCrgZmILFzTffrMsuu0xDhgyp8hr7u322b9+u1q1bKzU1VVdffbX27NnjeY1yd05ZWZluvfVW9e3bV2effbakwKx3gmp6Y6C+Tj74oqKitH79+ioH39tvv61HHnmkwuc2b96s3r176+jRozrllFO0dOlSde7c2YmfEHTqqgRr2g5S1W3BdrDG5s2b1aFDB09Zf/XVV2rZsmWVf7rhnVdffVUbNmzQunXrqn2d/d0evXr10sKFC9WxY0dlZ2drzpw5uuiii7RlyxbFxsZS7g66+eabtWXLFq1atcqzLBDrHZLfABEfH6/w8HAdOHCgwvIDBw6oVatWDkUVOuo6+LZt26Z9+/Zp8ODBFT7XsWNHZWZmKjc3V2+88YbGjRunFStWUCF6ob6VYHXbgu1gjU2bNqlr166e51999VWF5/De999/r1tuuUXLly9X48aNq30P+7s9hg0b5vm7W7du6tWrl9q1a6fXXntNN9xwA+XukEmTJunf//63Vq5cqdNPP92zPBDrHbo9BIhGjRrp/PPPV0ZGhmdZWVmZMjIyuOxigboOvmXLlmno0KFVTmKNGjXSGWecofPPP1/p6enq3r27Hn30Ub/FHazKK8GPP/7Y50qwum3BdrDGpk2b1K1bN8/zr776qsJzeG/9+vU6ePCgzjvvPEVERCgiIkIrVqzQY489poiICJWWlrK/+0mzZs3UoUMH7dixQxL1jL+ZpqlJkyZp6dKl+uijj5SSklLh9UCsd0h+A8i0adP0zDPPaNGiRdq2bZsmTpyogoICXXfddU6HFvTqOvjefvttjRw5ss71lJWVqbi42JYYQ4EVlaA324Lt4LuysjL973//q1DeO3fuVHJysnNBBbHBgwdr8+bNyszM9Dx69Oihq6++WpmZmQoPD2d/95P8/Hzt3LlTiYmJkqhn/O3mm2/Wv/71L7388suKjY3V/v37tX//fhUVFQVuvePIbXao0eOPP262bdvWbNSokXnBBReYa9ascTqkoFdaWmrGxMSY77zzjmfZ6aefbj766KOmaZrmgQMHzMjISPPQoUMVPjdjxgxzxYoVZlZWlrlp0yZzxowZpmEY5gcffODX+IPJxIkTzbi4OPOTTz4xs7OzPY/CwsI6t4NpVr8t2A7W+Pbbb01J5nfffedZNmzYMLNZs2bmqlWrHIwsdJw82gP7u31uu+0285NPPjGzsrLMzz77zBwyZIgZHx9vHjx4kHJ3gGoYWeP5558P2HqH5Bchr66D79lnnzX79u1b5XPXX3+92a5dO7NRo0ZmQkKCOXjwYCrCOjS0EqxuW7AdECxOTn7Z3+0zevRoMzEx0WzUqJHZpk0bc/To0eaOHTtM06Tc4R0muYDrjRgxQv369dP//d//OR2K67Et4Cbs786g3EGfX7hev379NGbMGKfDgNgWcBf2d2dQ7qDlFwAAAK5Byy8AAABcg+QXAAAArkHyCwAAANcg+QUAAIBrkPwCAADANUh+AQBVjB8/XmlpaU6HAQCWI/kFgCD0+eefKzw8XJdddpnToQBAUCH5BYAg9Nxzz2ny5MlauXKl9u3b53Q4ABA0SH4BIMjk5+dr8eLFmjhxoi677DItXLjQ89onn3wiwzCUkZGhHj16KCYmRn369NE333xTYR1/+9vf1LJlS8XGxurGG2/UjBkzdM4559T4nWVlZUpPT1dKSoqio6PVvXt3vfHGGzb9QgCwD8kvAASZ1157TZ06dVLHjh11zTXX6J///KcqT9b55z//WXPnztWXX36piIgIXX/99Z7XXnrpJd1777164IEHtH79erVt21ZPPvlkrd+Znp6uF154QU899ZT+97//aerUqbrmmmu0YsUKW34jANiF6Y0BIMj07dtXV111lW655RYdP35ciYmJev311zVw4EB98sknGjRokD788EMNHjxYkvSf//xHl112mYqKitS4cWNdeOGF6tGjh5544gnPOvv166f8/HxlZmZKOnHD25EjR/TWW2+puLhYLVq00IcffqjevXt7PnPjjTeqsLBQL7/8sl9/PwA0BC2/ABBEvvnmG61du1ZjxoyRJEVERGj06NF67rnnKryvW7dunr8TExMlSQcPHvSs44ILLqjw/srPT7Zjxw4VFhZq6NChOuWUUzyPF154QTt37rTkdwGAv0Q4HQAAwHvPPfecjh8/rtatW3uWmaapqKioCi25kZGRnr8Nw5B0ot9ufeTn50uS3n33XbVp06bCa1FRUfVaJwA4heQXAILE8ePH9cILL2ju3Ln67W9/W+G1tLQ0vfLKK+rUqVOd6+nYsaPWrVuna6+91rNs3bp1Nb6/c+fOioqK0p49ezRgwID6/wAACAAkvwAQJP7973/rp59+0g033KC4uLgKr11xxRV67rnn9NBDD9W5nsmTJ2vChAnq0aOH+vTpo8WLF2vTpk1KTU2t9v2xsbGaPn26pk6dqrKyMvXr10+5ubn67LPP1LRpU40bN86S3wcA/kDyCwBB4rnnntOQIUOqJL7SieT3wQcf1KZNm+pcz9VXX61du3Zp+vTpOnr0qK666iqNHz9ea9eurfEz99xzjxISEpSenq5du3apWbNmOu+883TnnXc26DcBgL8x2gMAQEOHDlWrVq304osvOh0KANiKll8AcJnCwkI99dRTuvjiixUeHq5XXnlFH374oZYvX+50aABgO1p+AcBlioqKNHz4cG3cuFFHjx5Vx44dddddd2nUqFFOhwYAtiP5BQAAgGswyQUAAABcg+QXAAAArkHyCwAAANcg+QUAAIBrkPwCAADANUh+AQAA4BokvwAAAHANkl8AAAC4BskvAAAAXOP/A67W0GYkpiRUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "plt.plot(angles, bell_values, \".\", label=\"Bell Inequality Value\")\n",
    "\n",
    "plt.grid(which=\"major\", axis=\"both\")\n",
    "plt.yticks([-1.5, -1, 0, 1, 1.5], [\"$-1.5$\", \"$-1$\", \"$0$\", \"$1$\", \"$1.5$\"])\n",
    "plt.xticks([i*np.pi/3 for i in range(7)], [\"$0$\", \"$\\pi/3$\", \"$2\\pi/3$\", \"$\\pi$\", \"$4\\pi/3$\",\"$5\\pi/3$\",\"$2\\pi$\"])\n",
    "plt.xlabel(\"Angle\")\n",
    "plt.ylabel(\"Reduction to Bell Inequality Value\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run on a QPU\n",
    "\n",
    "To run CHSH inequality on a QPU, we replace the LocalSimulator with an AwsDevice. \n",
    "To reduce the cost, we run the the experiment only the defualt angles which gives the maximum CHSH inequality value.\n",
    "\n",
    "The cost to run this experiment is \\\\$0.3 per task and \\\\$0.00145 per shot on the IQM Garnet device. Since we have four circuits of 1000 shots each, the total cost is \\$7 USD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Uncomment the following to run on QPU\n",
    "# from braket.aws import AwsDevice\n",
    "# device = AwsDevice(\"arn:aws:braket:eu-north-1::device/qpu/iqm/Garnet\")\n",
    "# circuits = create_chsh_inequality_circuits()\n",
    "# tasks = run_chsh_inequality(circuits, device, shots=1_000)\n",
    "# chsh_value, results, E_a1b1, E_a1b2, E_a2b1, E_a2b2 = get_chsh_results(tasks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that CHSH inequality is violated, so the device is demonstrating quantum behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated cost to run this example: 0.00 USD\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"Estimated cost to run this example: {tracker.qpu_tasks_cost() + tracker.simulator_tasks_cost() :.2f} USD\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Charges shown are estimates based on your Amazon Braket simulator and quantum processing unit (QPU) task usage. Estimated charges shown may differ from your actual charges. Estimated charges do not factor in any discounts or credits, and you may experience additional charges based on your use of other services such as Amazon Elastic Compute Cloud (Amazon EC2)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "5904cb9a2089448a2e1aeb5d493d227c9de33e591d7c07e4016fb81e71061a5d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}