{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2c935a4a-aa7d-499a-b358-3ef7fac417c2",
   "metadata": {},
   "source": [
    "# Obtaining the deuteron energy descibed on two qubits (Dumitrescu et al., PRL 120, 210501 (2018))\n",
    "\n",
    "Our starting point is the 2‑qubit Hamiltonian\n",
    "\n",
    "$H_2 = 5.906709\\, I + 0.218291\\, Z_0 - 6.125\\, Z_1\n",
    "      - 2.143304\\left(X_0X_1 + Y_0Y_1\\right).$\n",
    "      \n",
    "The trial state is \n",
    "$\n",
    "|\\psi(\\theta)\\rangle =\n",
    "e^{\\,i\\frac{\\theta}{2}\\left(X_0 Y_1 - X_1 Y_0\\right)} |00\\rangle,$\n",
    "\n",
    "1. Define the Hamiltonian and the two qubit circuit\n",
    "2. For a given theta, prepare the ansatz |\\psi(\\theta)\\rangle\n",
    "3. Compute the expectation values for $Z_0$, $Z_1$, $X_0X_1$, $Y_0Y_1$ with a finite number of measurements.  \n",
    "4. Deduce the expectation value of $E(\\theta)=\\langle H_2\\rangle$\n",
    "5. Then reproduce the Fig. 2 of (Dumitrescu et al., PRL 120, 210501 (2018)) where the different Pauli strings are shown as a function of \\(theta\\), show also \\(E(\\theta)\\) as a function of \\(theta\\), and deduce the ground state energy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4766d432-a17d-4e44-8f69-1734e11c0ff6",
   "metadata": {},
   "source": [
    "## 0. Imports\n",
    "\n",
    "- `QuantumCircuit` to build circuits\n",
    "- `StatevectorSampler` for **shot-based** counts (like running the circuit many times)\n",
    "- `Statevector` for exact (ideal) checks (no shot noise)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "09303518",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit\n",
    "from qiskit.quantum_info import Statevector\n",
    "from qiskit.primitives import StatevectorSampler\n",
    "from qiskit.visualization import plot_histogram\n",
    "import numpy as np\n",
    "\n",
    "# for vizualization \n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae79692c-9c5b-4cf1-bcd1-049d2d064d72",
   "metadata": {},
   "source": [
    "## 1. Circuit for the trial state vector\n",
    "\n",
    "The trial state is \n",
    "\\[\n",
    "|\\psi(\\theta)\\rangle =\n",
    "e^{\\,i\\frac{\\theta}{2}\\left(X_0 Y_1 - X_1 Y_0\\right)} |00\\rangle,\n",
    "\\]\n",
    "It could be prepared by \n",
    "1) performing a X gate on the first qubit\n",
    "2) a Y-rotation with angle \\( theta \\) on the first qubit\n",
    "3) a CNOT where 1 control 0\n",
    "\n",
    "**Task: For a fixed $\\theta$, make the circuit and perform measurements by counts** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47e471f4-5345-405d-91cf-c474a4257ecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here Make the circuit\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a86c0e-08a2-468f-8c8f-7782c1289d2d",
   "metadata": {},
   "source": [
    "## 3. Computing expectation values of observables for one single angle value \n",
    "\n",
    "To measure on one qubit:\n",
    "\n",
    "- **Z**: do nothing, measure in computational basis\n",
    "- **X**: apply **H**, then measure Z\n",
    "- **Y**: apply **Sdg** then **H**, then measure Z\n",
    "\n",
    "For two-qubit Pauli strings (like $X_0X_1$), apply the rule on each qubit.\n",
    "\n",
    "### From counts to expectation values\n",
    "\n",
    "Eigenvalues:\n",
    "- bit `0` → +1\n",
    "- bit `1` → −1\n",
    "\n",
    "Qiskit bitstrings are shown as `'b1b0'` (highest classical bit on the left). \n",
    "Tips: You need to convert counts to expectation values, you can use a similar \n",
    "routine used in the two-qubit tutorial. Examples of routines are also given below.  \n",
    "\n",
    "**Task: For a fixed $\\theta$, make the circuits and computes the different expectation values** \n",
    "\n",
    "In this intermediate version, the example of $\\langle Z_0 \\rangle$ is given as an example. \n",
    "The participant should do similarly $\\langle Z_1 \\rangle$, $\\langle X_0 X_1 \\rangle$ and $\\langle Y_0 Y_1 \\rangle$. \n",
    "\n",
    "The routines from counts to expectations are given. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "63288326-c79b-4d9c-9e3a-b62473ee9ec2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "counts_z : {'01': 743, '10': 257}\n",
      "\n",
      "Expectation values (from shots):\n",
      "<Z0>   = -0.486\n"
     ]
    }
   ],
   "source": [
    "sampler = StatevectorSampler()\n",
    "Nshot = 1000\n",
    "theta = np.pi/3\n",
    "\n",
    "# -----------------------\n",
    "# 1) Define circuits independently\n",
    "# -----------------------\n",
    "\n",
    "# Z0 (and Z1 measurement) (computational basis)\n",
    "qc_z = QuantumCircuit(2, 2, name=\"Measure_Z0Z1\")\n",
    "qc_z.ry(theta, 1)\n",
    "qc_z.x(0)\n",
    "qc_z.cx(1, 0)\n",
    "qc_z.measure([0, 1], [0, 1])\n",
    "\n",
    "# -----------------------\n",
    "# 2) Run shots and get counts\n",
    "# -----------------------\n",
    "#sampler = StatevectorSampler()\n",
    "\n",
    "counts_z  = sampler.run([qc_z],  shots=Nshot).result()[0].data.c.get_counts()\n",
    "\n",
    "print(\"counts_z :\", counts_z)\n",
    "\n",
    "# -----------------------\n",
    "# 3) Convert counts -> expectation values for Z0, Z1, X0X1, Y0Y1\n",
    "# -----------------------\n",
    "def bit_to_ev(b: int) -> int:\n",
    "    # Z-eigenvalue mapping\n",
    "    return 1 if b == 0 else -1\n",
    "\n",
    "def exp_Z0(counts: dict, shots: int) -> float:\n",
    "    # Qiskit bitstrings are 'b1b0' (bit 1 left, bit 0 right)\n",
    "    return sum(bit_to_ev(int(s[1])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_Z1(counts: dict, shots: int) -> float:\n",
    "    return sum(bit_to_ev(int(s[0])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_two_qubit_product(counts: dict, shots: int) -> float:\n",
    "    # After basis rotations, Z-measurement corresponds to X or Y measurement.\n",
    "    # The product eigenvalue is (ev on qubit0) * (ev on qubit1).\n",
    "    total = 0.0\n",
    "    for s, c in counts.items():\n",
    "        b1 = int(s[0])\n",
    "        b0 = int(s[1])\n",
    "        total += (bit_to_ev(b0) * bit_to_ev(b1)) * c\n",
    "    return total / shots\n",
    "\n",
    "Z0 = exp_Z0(counts_z, Nshot)\n",
    "\n",
    "print(\"\\nExpectation values (from shots):\")\n",
    "print(\"<Z0>   =\", Z0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f647fec-3bbe-4311-8e78-9d8aeca8e207",
   "metadata": {},
   "source": [
    "## 3. Computing expectation values of observables for a set of angles\n",
    "\n",
    "Tips: \n",
    "\n",
    "### Making loops on $\\theta$\n",
    "\n",
    "- You should define a list of angle values and make a loop on the angles, for instance:  \n",
    " \n",
    "``` \n",
    "n_theta = 41  # number of theta points\n",
    "theta_list = np.linspace(-np.pi, np.pi, n_theta)\n",
    "```\n",
    "\n",
    "### Example of useful routines for the $\\langle Z_0 \\rangle $ expectation value\n",
    "\n",
    "- Prior to the loop, you can define a routine that builds the circuit for a given $\\theta$ value. \n",
    "An example of such a circuit to measure $Z_0$ (or $Z_1$) is \n",
    "\n",
    "```\n",
    "def circuit_z(theta: float) -> QuantumCircuit:\n",
    "    qc = QuantumCircuit(2, 2, name=\"Measure_Z0Z1\")\n",
    "    qc.ry(theta, 1)\n",
    "    qc.x(0)\n",
    "    qc.cx(1, 0)\n",
    "    qc.measure([0, 1], [0, 1])\n",
    "    return qc\n",
    "```\n",
    "\n",
    "The routine can be called in a loop on $\\theta$ value for instance using the syntax\n",
    "  \n",
    "```\n",
    "counts_z  = sampler.run([circuit_z(theta)],  shots=Nshot).result()[0].data.c.get_counts()\n",
    "```\n",
    "\n",
    "### Storing\n",
    "\n",
    "To store results as a function of $\\theta$ this could be done by defining a list: \n",
    "\n",
    "```\n",
    "Z0_list = []\n",
    "```\n",
    "\n",
    "Results could be stored in the $\\theta$ loop using \n",
    "\n",
    "```\n",
    "    Z0_list.append(exp_Z0(counts_z, Nshot))\n",
    "```\n",
    "\n",
    "with `exp_Z0` the routine to transform bits into qubits and compute expectation values (given below).  \n",
    "After the loop, you do \n",
    "```\n",
    "Z0_list = np.array(Z0_list)\n",
    "```\n",
    "that are the $\\langle Z_0 \\rangle$ results.\n",
    "\n",
    "- As an intermediate step, just do the code that shows the expectation value of $\\langle Z_0 \\rangle$ and display this expectation value as a function of $\\theta$ and display it\n",
    "\n",
    "To display the result, use matplotib standard routines. For instance, a minimal script is:\n",
    "\n",
    "```\n",
    "fig, ax = plt.subplots(figsize=(6.0, 4.0))\n",
    "ax.plot(theta_list, Z_0_list, marker = 'o', label=r\"$\\langle Z_0\\rangle$\")\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "- Then, do the same for $\\langle Z_1 \\rangle$, $\\langle X_0X_1\\rangle $ and $\\langle Y_0Y_1\\rangle $.\n",
    "- Note that, while for $\\langle Z_1 \\rangle$, the same circuit can be used,for $\\langle X_0X_1\\rangle $ and $\\langle Y_0Y_1 \\rangle$ additional circuit should be built with proper rotations before measurements.   \n",
    "\n",
    "### Task: \n",
    "In this intermediate version, the participant should do a similar work to get the $\\langle Z_1 \\rangle$, $\\langle X_0 X_1 \\rangle$ and $\\langle Y_0 Y_1 \\rangle$ expectation value, store them as a function of $\\theta$ and display them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cbbd0d97-b74c-4fb3-9e3b-3bc89749a707",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAGACAYAAABWTZ3rAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZgtJREFUeJzt3Qd0FNX3B/CbhFRSKAkJvQQChBbpSegBFRFEQUCpFhQEewHUH9W/iKCggoAggoBIlSLSkd5LQgklkARCIKRAei//c1+YdXezabubbfP9nLNnd2dmZyeb8m7eu+8+q4KCggICAAAAAAXr/x4CAAAAAAIkAAAAAA3QgwQAAACgBgESAAAAgBoESAAAAABqECABAAAAqEGABAAAAKAGARIAAACAmkrqG0B/4uPjac+ePdSgQQNydHTERwsAAGAgGRkZFBkZSc888wy5u7uX+/UIkCrQzp07acyYMRX5FgAAAFCClStX0ujRo0mWAVJiYiJNmDCB/vjjD4qIiBA9NqagZs2a4n7ZsmXUtm1bkoNBgwbR5s2bSS7w9Vqu1NRU6t69Ox0+fJicnZ1JDuT28yzHr1lOX++FCxdo7NixirZYdgHS3r17ady4ceTk5KRVYDV58mQ6ePAg2draUrVq1WjOnDkUEBBQ5NgTJ06IYxMSEignJ4d69uwpjq1SpUqx53dwcBD3Pj4+sgmQeChRLl8rw9druZKTk8W9n58fubq6khzI7edZjl+znL7e1NRUlbZYdknas2bNEtHw4MGDy/W6vLw86tu3L127do2Cg4Pp6tWrNHToUAoKCqLz58+rHMvPefuQIUPEcXz89evXxev5PKUJfZBCefnyWBOYe/LkBF8vWBK5/TzL8WuWy9ebl18g2l5dWBUUFJh1y52bm0uVKlWi6dOn04wZM8o8xLZq1SqRH3Tq1Cnq1KmTYnuzZs2oVq1aoldJwr1FMTExIpiS8Ov8/f3FeUaNGqXxPY4cOSK66L1GL6AGTVvStP6+9GzLmsV+M89EPKLYlEyq4eJAHRtWIxtrq3J+GgCgzx4kNzc3SkpKkk0PEoCpyitHG7n7ygOasSOUIm9coZhVH4hh8m7duslviI2DI21s2LCBXFxcqGPHjirbuado8eLFFBsbSzVq1KCHDx+KD/edd95ROY5fx69fv359sQGSnZ2d4nFMUiaNX3OBFo9oWyRIkr6ZD5IyFdtqujmUGFABQMWyt7enadOmiXsAMJ7ytJF8LLe1BcW0xbIaYtMWD5M1bNiQrKxUI1Bvb2/iTrVLly6J53zPz3m7Mmtra9FTFRISUux7KH9TCp7cpm67SkkZOeKcyt9M5W+8ckDF+wHA8Dgw4p5pBEgAxlOeNpJ7maZtv6oSHOkSIJl9D5K24uLiqFGjRkW2S13pvF/5XlMXO2+7ceNGud43NiWL2szYS7Y2VlTF0ZYep+cU+WYy3sahG0fNfXy9MNwGAACykpdfINrA4tpI9vGGENpy4R5FJ2ZSRFwqpefk6+39ZduDZGw5eQUUl5pNuSUkb/Mejpp53BUAAEBOzkQ8KtJzpC4tO4/2hsbS1fvJeg2OZN2DxFU1pWm8yqRtHh4eiuOUt6sfKx1XHstG+JFvnWq09WI0zd1Teg9UbHLJPyAAAACWJjalbG3foLa16blWNelhYhp9vu2/yVS6km0PEtc24RLk6pP4wsPDRV5S69atxfM2bdqI57xdWX5+vng97y8rqyeJZb18a1HtKo7Utl7VMr1u+bFwuhWr23RFAAAAcxtpKYvB7epSUHNPGtqpoWhj9TX/WxYBEtcq4llpyrimEfcAnT17VmU7T+/nqfk8g415enqK6YEHDhxQOY5fl5KSIs5TFtI3jLPupamJPE2xLN/My9HJ9OyCozRzR6hI8BZfU34BnbydQNuCo8W9XOosAQCA5dsWHE3/23q5TJ0O3JYyblu5jZX26UoWARJP0edS41wNWzJy5Ejq3LkzTZo0idLT08U2nt7PdZTmzZun8np+ztt//vln8ZyP56ra/PoRI0aU6Rq83ByKTPEv6Ztp9eQ2c0AL6t3cU+QqrTgeQb3mHaIv/rpMgd8cpFeWnaL3/wwW913mHMSMNwAAMGtZuXn05dbLom3LyMknH0/nYttI9U4Hxm0st7Xc5pLcA6SZM2eK4bIlS5aI588995x4ztP4JdwLxEuCKM9Es7GxoV27dollQPj4Fi1aiLXc9u/fT+3atVN5j/bt24sepHXr1onj+PgmTZqI1/N5SjP7xVZ0bFIvjTWNivtmSgHVqIAGtHx0e1r1ekfy9qhMCWnZtPb0XYpRy0tCWQAAADAXeRpGQaIepdPgxSdpzam74piJPRvTrve70ZIS2sji2lVuc7ntlXUlbVNfKI+DLV6qpLS1b8pSJTQzJ486/N9+SsnM1XgOqyc/NPyDgSrcAABgLoUfqzrZijaOe4348fdD/ahn08JUF21XmyhPG6yJbGexGXr1ZF4gkNfAKW4dHP5G+3tXL/E8F+8mFhscqZcFKO1cAAAAhrZbQ6VrxjUBWUN3J1r7ZmeqVcWx3G2kZNGiReKWkZGh07UiQDIAXkxXH6snl3XKY1mPAwAAMIXCj5LMnHzydNUtf0jqjJB6kGSbgyQn3K1YtuOwdhQAAJhf4UdTKo6MAMmMlLUswI5L9yknT78VRQEAAOQ0CoIAyYyUVhZA8sfpKBrz2xlKejKmCwAAYGxVnWz1OlpS0RAgmZmSygLwVMjlo9pTZTsbOn4rgV78+TiFx6Ua7VoBAACkJbO+33eTylP40diQpG2mQVIfX69ipzxuGh9Ab646R+HxafTizydo8fC21KlR9XJPkQQAANBVcFQivb36HD1MziInOxtKz84TwVBBGQo/GhMCJDNV0pTH5jVdaeuEQPEDeeFuIo349TQ521eiZKUSARyl8w+ipiJbAAAA+rDp/D36/K/LlJ2bT01qONOyUe3pekxykTpIXibYJqFQZAWSphg2aNCg1DpIFYGLbo1acUbjjAApPi+uEikAAEBZqRdybFuvCs3ZfUMskcX6+HrS/KF+4p91Tcfrc1RDuQ4SLyqvbaFIBEgVSNcqnrriH8DAOQfFMiSaoPI2AABURGVsOxtryn4ym/q9oCb0QVATsjbw0JmubTCStC0YR+fFBUfqlbcBAAC0rYz9QK2tkYKjt7s1pI/6+Bg8ONIHBEgWzNxqTgAAgGVVxt4e8kAcZ44QIFmwslfeNo2aEwAAYD7OmFll7PJCgCTzytuV7W1MpuYEAACYj1gLH6VAgCTTytuStKw8WnGscJYBAABAWdnZWFv0KAUCJAMYNGgQ+fr6immHplJ5m3uWBvrVFo//759rolYFAABAWVyJTqIZO66aZGVsbmu5zeW2VxeY5m/B0/yVaao5wZMKvv7nGi07GiF6m5aOaEe9fT2Nep0AAGDadl+JoQ/XB1NGTh55utqLCtnFVcY2Zq09TPOHclXefsGvtrjn51ZWVjSlb3Ma1LaOCKAm/HHBbJPpAACgYhUUFNCif2/RuDXnRXDUtYk77f2wu1gHVNP6oOZeiBhLjcgc16b4ZlArSkzPpgPXY+mNVWdp/Vv+5FvL1diXBgAAJjLq0LqOG3259Qr9dTFa7B8T0IC+7NecKtlYl7o+qLlCgARka2NNi4a3pZG/nqazkY9p9G9naMNb/hSTnGlRP+wAAKBdZWxbGyvKySsQ7cD0AS1oZOf6ZV4f1FwhQALBwdaGlo/uQEOXnqTrMSnUe/5hleJeWNwWAEA+lbEL1LZzcMQm9mxcJDiyVJjFBgpujrb0WkAD8Vi98ikvWcK/NPzLAwAA8qyMveFclNlWxi4vBEigwD/0Cw6EafxEpF8H/uWRyy8HAICcWHpl7PJCgGThdZD0+cuBxW0BACyXpVTGXqSnOkjIQTKAzZs3G70Okpx+OQAAoPwS03MsojL2hAkTxE2qgyTrACkzM5OmT59Of/31F9nb25OdnR1NnTqVBgwYUOLrevToQffu3SNnZ2eV7enp6RQWFqZS4HHMmDF07NixIse2aNGC1q5dS5YAi9sCAMizvtGK45H09c7QEo+zelLfSC7rd1pEgDRy5Ei6cuUKHT9+nNzd3WnHjh304osv0tatW+n5558v8bXLly8XgZKy7777jlasWFGk10fTsZa4uC0nZBeXZeTubCebXw4AAEuXnJlDn228RLuvxojnbetXoYt3EsVjTZWxeX1PuZR8MfscpMOHD9OmTZtEDxIHR6x///7Uu3dvev/990VkXJwvvviCmjdvXmT7L7/8QuPGjSO5Kcvitlm5+fQwGUNsAADmhCfXnLydQNuCo8U9P+f11Pr/dEwER1znaMaAFrR5XIDG9Tu9LKAytux6kDZs2CDug4KCVLbz8z179tC5c+eoQ4cOGl/bp0+fItv+/fdfMew2atQokiNpcVv1ImFerg7ECzdHJ2bSG6vO0aZx/lTZ3ux/fAAAZFn40dXBltKzcyk3v4BqV3EUxYL96lYR+yy1MnZ5mX0LFxwcTK6uroreI4m3t7e4DwkJKTZA0mTp0qX06quvkpubW5F9nGs0bdo0SkhIEHlOzz77LE2aNEnjseasuF+OB0kZNHDRcbr2IJk+WB8sFrflpUoAAMC8Cj/y0BprWduV1rzRiao42Vl8ZWzZBUhxcXEiQFInbeP95TkXJ3qfPHmyyD4XFxexuCvnN/G5OfB6+eWXadu2bXTq1CmxvzipqamUnJyseM6J5HwzZZp+OepUdaKlI9vTK8tO0b7QhzRnz3Wx2C0AAJhn4ceE1GxycbAlS5CVlSVuym2vrHOQ9Om3334jPz8/jVPyf/rpJ/rxxx8VgVebNm1o/vz5FBoaWmp9o+7du4teJuk2e/ZsMlft6lelbwe1Fo+XHg6njeeijH1JAACggdwKP86ePVulreW2V9YBEg+tKffOSKRtHh4eZToPJ3Nzcvb48ePL/N7+/v7inmfPlZZInpSUpLhNmTKFzNnAp2rTu70ai8ef/3XZYn65AAAsidxq202ZMkWlreW2V9YBEvf4cDDEeUHKwsPDFT09ZXHgwAF69OgRDR06tMi+vLw8jUN1NjY24j4/P7/Ec3PtJO55km6mPrxWFh/29qHnWnmJBQzfXn2OIuLSisyQAAAA45FbbTt7e3uVtla9bqHsAqQhQ4YoAhxl/LxRo0bUvn17RZATGxtbYnI2F4N0dHQssi8qKorq168vzqGMZ8ix8iSBWwpOzv7uZT9qXceNHqfnUO/5h0Vu0vt/Bov7LnMOYmFbAAATqG1XHJ5iw/tR285CAyQu3Dh48GBRByk+Pl5s27lzJ+3bt48WLFggEqvZO++8QzVr1qQTJ04UOcfDhw9p+/btJdY+ysjIoP/973+KIOn+/fv0ySefUN26dWnixIkkR452NvRqx3risXqPEReb5JkTPIMCAACMM9nmnR6FM7rVybHwo+xmsbHVq1eLACkwMFB0sdna2tKWLVtEwUiJp6cnValSReOMN66a3bVrV/Lx8dF4/lq1atGyZcto48aN1Lp1a5GvxMuRcDHKmTNnFikxIBccFP1wIEzjPg6X+FeOZ1BwyQD8AgIAGN6t2MKZXHaVrCk7N1+l8CMHR3Iq/FheVgUllZoGnUgL5Smv6WZJONeIh9NKs25sZ9nX0wAAMLTHadkU8M1BysjJo9Wvd6RKNtayKvx4Qcc22CJ6kEzdoEGDRG6TtMKwpZDbDAkAAHOy9vQdERz51nSlLk3cFSknlm7RokXixqkxukCAZACbN2+2yB4kuc2QAAAwF5k5ebTyxB3x+K1ujWQTHDGpM0LqQZJtkjYYf4ZESb92mCEBAGB4XHIlPjWLark5UL/WyDPSBgIk0BqPX3OSHysuSPqiX3OLH+cGADAl+fkFtOxohHj8epeGZMsrjUO54VMDnfAMiMUj2ooZEcqk3tz7ibqNAQMAQPkcuhkrZq+52FeioR3q4uPTEnKQQC9BEk/l5yVHpBkSdx6l0eTNl2n+vjDq17oW1a5StAAnAADoH6+TyV7tVM9iFqI1BvQggV7wMJq/d3V6wa+2uB/Sri51bFBNzKCYsf0qPmUAAAMIiUqk0xGPqJK1FY0JbIDPXAcIkKDCliL56sWW4pd0b+hD2hf6EJ80AEAFW3a0sPdoQJtaVNMNPfe6QIBkoDpIvr6+oi6DnPh4utCbXRuJx9O3X6X07FxjXxIAgMWKepRO/1wuXN5J+tsrR4sWLRJtLre9ukCAZKA6SKGhoRZVJLKs3gtqLPKPohMzil2WBAAAdLfieATxsphdm7iTb62iy2rJxYQJE0Sby22vLhAgQYVysqtEMwa0EI9/PRpBN2JS8IkDAOhZUnoOrT8bJR6PlXHvkT4hQIIK19vXk5729aTc/AL6cutlUaMDAAD0Z+2ZO5SenUfNvFxEDxLoDgESGMS0AS3Iyc6GzkY+pk3n7+FTBwDQk6zcPFp5PFLReySnZUUqEgIkMAjOQ/qwt494/PWua/QoLRufPACAHmwPvk+xKVnk6WpP/dvUwmeqJwiQwGC4Jgd3/yam59DX/4TSydsJYr0gvs/DsBsAQLkVFBQopva/FtiQ7CqhWdcXVNI2AJ5q6OjoqFhhWK54PaD/e7ElDVp8kjadjxY35UVteV03rsoNAABlc/hmHN18mEqV7WzolY718LFR4TR/vmVk6LbUFUJNA5DzNH91cSlZGrfHJGXS+DUXaPeVwhoeAABQPO515973r/+5Jp4P6VCX3ByxrIg+p/mjBwkM+gs9Y0eoxn08r43TCnk/r+vGS5cAAEBR/I8k/618kJSp2Pb3pQfUqWE19MLrEXqQwGB4MVvlX2hNQRLv5+MAAEBzcMS97ep/S+NTstALr2cIkMBgYlMy9XocAIAce+E1VZKTtvF+THrRDwRIYDA1XBz0ehwAgJygF96wECCBwXRsWE3MVisuu4i3834+DgAAVKEX3rAQIIHBcOI1T+VnVsV0EfN+JGgDABR1Oza1TB8LeuH1AwGSgeog+fr6iroMcsd1jhaPaEtebkWH0XxrumIGBgCAhmKQ8/fdpB8P3irxs0EvfCFua7nN5bZXF1YF/MlDhbhw4QK1a9eOzp8/T23btsWnrISTCHk8nbuM+Ufw4w0hlFdAtHGcP3VogCE2AACWmZNHH28MoZ2XCmvE9fH1pP2hD8Vj5cZb6pXnf0BRcFc/bbBF9CBlZmbS5MmTqWnTptS6dWtq3749bd++vdTXHTp0iKpVq0Z+fn5Fbrdv3y5y/I4dO8S5+T34vfg9+b2h/HgYzd+7Or3gV5sGPlWHhnQorAA7d88NETABAMjdw+RMGrr0pAiObG2s6NtBrWnZqPYae+H5OYIj/bKIQpEjR46kK1eu0PHjx8nd3V0EMi+++CJt3bqVnn/++RJfO2DAAFq5cmWp7yGd86+//qL+/ftTfHw8de3alcLDw2nDhg16/Grk6b2gxrT5wj3Rq3Q0LJ66+XgY+5IAAIzSq845RE52NvT26vMUk5xJVZ1safGIdtS5UXVxLPcQcUFd5eN5cgvyN/XL7AOkw4cP06ZNm+jPP/8UwRHjAKZ37970/vvvU79+/cjKSreqzNyjwefq06ePODfj95o2bRq98sordOTIEerWrZtevh65qunmSCM61acVxyNo3t4b1LWJu87fNwAAc62MLWlcw5l+Hd2e6levrLEXHiqO2Q+xSb03QUFBKtv5OffunDt3Tuf34HNERERofA+2fv16nd8DiN7p6S3+a7p0L4n2PhljBwCQY2VsyTs9vIsER2AYZh8gBQcHk6urq6L3SOLt7S3uQ0JCSnx9WFgYDR48mDp06EA+Pj40bNiwIkEVv4fyOSUeHh7k4uJS6ntA2bg729NrgQ3E4+/33qT8fOQiAYA8K2Mzqyd5maiMbRxmHyDFxcWJAEmdtI33F8fW1pays7PFUNnZs2dFYOTs7Ez+/v60e/dulfdQPqf6+5T0Hiw1NZWSk5MVt6wszSvaA9FbXb3JxaES3XiYQjsu3cdHAgAWC5Wx9YvbVuW2ltteWQdIuggMDBSBUatWrRTBzs8//0yenp70wQcf6O19unfvTm5uborb7Nmz9XZuS+PmZEtvdW0kHi/YH0a5efnGviQAgAqBytj6xW2rclvLba+sAyQeWuNIUZ20jYfBysPOzk7US7hx4wYlJCQo3kP5nOrvU9p7cCJ5UlKS4jZlypRyXZPcvNalIVWrbEcR8WliZhsAgCXC+pT6xW2rclvLba9RAiROgB47dqzIy2nUqPA//hkzZtDOnTvJkLhmEQcpUjCjfH2sTZs2xb720aNHYohNnY2NjbjPz89XvIfyOSU81T8lJaXE92A8bMe9U9LN3t6+zF+fHDnbVxKJiezHA7coKzfP2JcEAKB3PDW/hkvx7QEqY5cPt63KbS23vQYPkDgp+amnnhKzt/gipECCA4WJEyfS33//TYYyZMgQcX/gwAGV7fycAzcu7Mjy8vIoNjZW5ZiXXnpJlAhQxsdxUnb9+vUVPUOcwN2wYUON76F8DaA/IzrXJ09Xe4pOzKB1p+/iowUAi8NT9dvUraJxn1TkBOtTGo9WARJXkH7rrbfo4cOHdPHiRapatarYPnDgQJHcPHfuXH1fZ7F69OghZqFNnz5d9Ogw7sXat28fLViwQFFL55133qGaNWvSiRMnVF7/7bff0v37hcnAHOh9+eWXFBkZSd98843iGD7HDz/8QHv37lX0kHGPFfeY8XvrOs4JRTnY2tC7vZqIxwv/vU0Z2ehFAgDLEp+aRcfCCtstLgapDJWxzbRQ5M2bN2nXrl2K58oF/XgJjvT0dDKk1atXiwCJk665i41np23ZskVR1JFx4nWVKlVUZqJ9/fXX9Ntvv9HTTz9N1tbW9PjxY2rSpIn42p599lmV9+Bz8TmnTp0qxjl5iRGurM1BElSMIe3r0tIjtynqUQatOhlJ47qrllkAADBny46GU0ZOHrWu40ZbxgfQ2cjHqIxtQrRarJaHrpTzcTipmReFk9SrV4/u3sWwCBar1d3m8/fEQo1VnGzpyGc9ydVB9b8sAABzlJCaRV3m/CsCpBVj2lOvZp7GviSLc8EYi9VyYjYPs+Xk5BTZN3PmTPL19dXmtBZr0KBB4jNZtGiRsS/F7Ax8qjZ5e1SmxPQcmrH9Km0LjqaTtxNQOA0AzNovSr1HPZvWMPblWJRFixaJNpfbXoP3IJ0+fVrk/nCGOCdBcy0hzsPh5O0HDx7Q0aNHtYrWLA16kPTj/3aG0rKjESrbaro5iORFXrQRAMCcoPfIgnuQOnXqJOoLcITGydA8XX7btm1Ut25dsR3BEehznaLlasERi0nKFOsX8X4AAHPC//Ch98hCk7RZx44dRTCUkZEhkpt5Jpujo6N+rw5kraR1ingbTw3g/X18vcR0WQAAc+g9+v1kpHj8flATlUlOYFp0rqTNQVGtWrVUgiOe7QWgK6xTBACW2HuUnp1HrWq7Ua9myD0yZRWy1MhXX31VEacFmcE6RQBgqb1HH/RG75FFDrFJS3EAVCSsUwQAlgS9RzIIkGrUqEHjxo1T2cZrkl27do2uXLlCr732mr6uD2S+ThHPVuOE7OKmWvJ+Pg4AwJQ9SstG7pEcAqSgoCCaNm2axn08xZ9ntsF/uBYD52hNmDBB3KBsOPGap/LzbDVOY9QUJHVt4o4EbQAwi6rZnHvUsrYrBTVH7lFF10HiG08iM3gdpNL4+/vTyZMnSe5QB0k/eCo/z1Z7kJSp2ObqUImSM3PJvpI17XyvKzWuoduqzQAAFdl71GXOQREgLR/Vnnr7omq2ObTBWk/zLw4Xirx3756+TwsyxsUgeSo/z2rjxG3OTWpfvyq9vuosHQ2Lp482BNPm8QFka1Mhcw4AAHSC3iPzpFWA9PrrrxfZxh1RXDDy+PHj1LdvX31cG4DKcJu/d3WVT2Tu4Db09PzDdOleEi08eIs+7OODTwwATK73aNWJJzPXgnxQ98jSA6S1a9eK2kfqM9s8PT1p/Pjx9Nlnn+nr+gCK5eXmQLMGtqT3/wymhf/eop7NapBf3Sr4xADAJArdcq/3imOFuUctarkg90gOARIvMXLx4kX9Xw1AOb3gV5v2X4ulHSH36aP1wSIfydEOZSgAwLTyJu8nZtKeqzFYP9KMaJW0sXr1av1fCYCWZr3Qgjxd7Sk8Po2+2XUNnyMAGDU44pm3ysERS0zPwfqRcgiQWrZsWeL+p59+WtvrASi3Kk52Ih+JrTp5h47cjMOnCAAmt34k4/18HFjIENvMmTPLdVIuFgn/QR2kitfNx4NG+den30/eoU83hdCeD7qJwAkAwBTXj1SfdAKmVwepTAHS9OnTy3VSrE6savPmzVrVYIDymdK3OR0LixdDbV9uvULDO9VXlAXgats8Ew4AoKJg/UjTIBVlluogVWiA1KZNm3IlZT/11FNaXxCAtjg5+/uhfvTSz8fp70sPxE15SRKuys01lQAAKgLWj5RhDtKnn35arpOW93gAfYlJyiBNw/u8nhsnTnICJQBAReCeai4/Uhzuw8b6kRYWIL366qvlOumZM2e0vR4AnRMkNUGCJABUNB7GH9a+rsZ90gA/92RjuN88aL3UCFfO5vVNbt++TVlZWSr7tm7dSgsWLNDH9QGUGRIkAcDYLkcniXsnOxtRIFLCPUsY5pdBgBQTE0P9+/cXARInZCuvd4sEbTAWJEgCgDGFx6XSgeuxZGVFtH1iF4pLycJEEbnVQZo8eTJ17tyZQkNDqVmzZhQRESFu0jps3333nf6v1Myn+XP1cZ52CBUHCZIAYEwrn6y5FtSsBjWu4Sym8nO1f77HsJrhcFvLbS63vQbvQQoODha9R7z+mr29PdWvX19s5/uNGzfSgAED6IMPPtDpwiwJpvkbLkGSEyA5IbugmBwA7ubm4wAA9CkpPYc2nrsnHr8e2BAfrgVM89eqB6lSpUoiOGK5ubkq+xwdHSkuzrCVjDMzM0WvVtOmTal169bUvn172r59e6mvu3TpkvgQOdJs1aoVNW/enAYPHiy2qxszZgw1btyY/Pz8VG7Dhw+voK8Kyov/Q+Mxfqap4hEHTUiQBICKsP7cXcrIyaNmXi4oAmkhtAqQ8vLy6NGjR+Kxl5cX/fHHH4p9f//9NyUmJpIhjRw5krZt2yaG+Di4mTZtGr300kviWkoyZMgQkWR+4sQJunz5Mp07d45ycnKoY8eOdPbs2SLHL1++XPSeKd/Wrl1bgV8ZlBfXOVo8oq3GqbaOttbk7+2ODxUA9Co3L59Wnbij6D1CLq6MA6Ru3bpRly5dKDo6msaOHUsjRowQPTfcozJw4MBylwXQxeHDh2nTpk2i2re7e2HjxwnkvXv3pvfff18lgVyT2bNnU5UqVcTjypUr07x588SsvJ9++skg1w8VEyQdm9SL1o3tTD8M86O1b3aixh6VKSMnn5YfDcdHDgB6tefqQ4pOzKDqle1ogF8tfLpyDpCmTp1Kv/32G1WrVo1efvll+uWXX8jV1ZVsbW3FvhkzZpChbNiwQdwHBQWpbOfn4eHholeoONzbpF71u06dOuL+8ePHFXK9YLjhNilBMrCxO33yTFOxfcWxCEpIVS1LAQCgixXHI8T98M71ycG2MP0EZBog8TT/Tp06iXwj9uabb9KxY8fEsBQHSBwoGQoPc3FwJvUeSby9vcV9SEhIsa+1syu6mOmNGzfEfa9evYrs4+G07t27U8uWLcXaap9//jklJRXWvADT9kwLL2pV243SsvNoyeHbxr4cALAQwVGJdP7OY7K1saIRnesZ+3LA2AEST51TLw5pLJwQzgGSOmlbeRPGf/75Z/Lx8aG3335bZbuLi4sICHfs2EFXrlwRPWg8tBcQEEApKSklnjM1NZWSk5MVN1P57OSEcwI+ftpHPP795B16mFz8itsAAGX125Peo/5tapW51AhUDG5bldtabnsNHiDdu3dPzOh65513xHR/S7Fz507asmWLmJbv5OSkso9zkn788UdF4MUL+M6fP1/UgiqtvhH3Orm5uSlunPcEhtfdx4M6NKhKWbn59NPBMHwLAEAnXFJk55NFsTG13/i4bVVua7ntNXiA1KRJE7p+/bqoL/Dee++JPJ6FCxcafPYa46E1jhTVSds8PDzKdJ4jR46IXqNdu3aJIbSy8Pf3F/c8e660RHIeipNuU6ZMKdP5oSJ6kQpzkf48E0VRj9LxEQOA1lafiqTc/AJRW61lbTd8kkbGbatyW8ttr8EDJO414hlfb7zxhggOeJr/nTt3FHWBDh48SIbC78nBUEJCgsp2TtCWenpKs2/fPho1apQoC9ChQweNZQ00DdVJtaDy8/NLPL+zs7PoeZJuXFwTjKNzo+rUtYm7+KO2YD96kQBAOxnZebT29F3x+I0uKAxpCrhtVW5rue01eIBkba36Mi6w+OKLL4rp/+vWraM+ffqQoXAtI3bgwAGV7fy8UaNGomikFOTExsYWeT3nFEk9RxxssQcPHohq4JKoqChRJZzPoUyaIacpqALTJfUi/XXxHt2KLTl/DABAk78uRlNieg7VreZIvZt74kOyQNbaFmZk3KvC6661aNGCunbtKqpXv/XWW3Tq1CkylB49eojq11wHKT4+XpFLxL1CCxYsUBTs4nypmjVriqKQyiUC+LXce8S9YmvWrBG39evXF6mmnZGRQf/73/8UQdL9+/fpk08+obp169LEiRMN9vWC7vzqVqE+vp6UX0A0H71IAFBOXF9Pmto/JqAh1lmzUFqtxfbvv/+KmWw8JMWVpzk44lldXBNJmvpvSKtXrxYBUmBgoOhi4zIDnGzNBSMlnp6eoiCk8oy3cePGUXZ2tsa6TdL6cqxWrVq0bNkysc4cF8TkX4709HRRjHLmzJlFSgyA6eMZbfuvPRQJlu/0SKIWtZA/AABlczQsnm7FppKzfSUa0r6wdh5YHquC0kpNFzPExgHH6NGjRR4SJ21DUdJCedw7xXWTwLS8t+4ibQ+5L1be/nUMhkkBoGzG/HaGDt2Io9cCG9C0/i3wsVloG6zVEFuzZs3EMiPffPMNgqMy4N42XhC3tHIAYFgf9G4iusYPXI+lC3dROR0ASsd5ixwccfbGawFIzjZF3NZym8ttry60CpB4aE09URuKx3WVuF7ShAkT8DGZkEYezjSobW3xeN6ewgrqAACa5OUX0MnbCTRt21XxvHezGlSvumq9PDAN3NZym8ttr8FzkHh2GIAleC+oiZiNcuJ2Ai0/Ek4ervaiGi7XNeHeJQCA3Vce0IwdofQg6b8K/OfvJortvDg2WCatAiQAS1GnqhMFervToZtx9NU/1xTba7o50LT+vvjjByBzHASNX3OB1JN1H6dli+2LR7TF3wkLhXEyILn/8ePgSNMSAvzHj/cDgHyH1bjnSNNMJmkb7+fjwPIgQAKS+x8/TfDHDwDORDxSGVbT9HeC9/NxYHkqJEDiddoATB3++AFASWJTMvV6HJiXCgmQXn311Yo4LYBe4Y8fAJSEJ2zo8ziQQZJ2WloazZ07V6x3FhMTU2SNMl6GA/7DtRi4wjhPPcRUf9OBP34AUBKezerl6kAxyZp7iHieq5db4axXMK06SHzjJcIMHiCNHz9erLvGS3vwlH/lmkhcmJvrJMF/uBYDKmmbHv6jxrPVOCFbU4ol/vgByBuX+ujVzIP+OBNVZJ9UBIRnu6IkiGmROiOkStoGDZAOHjxIV65coTp1NK9B06dPH60vCMBQ+I8a/3Hj2Wr8x049SOLn+OMHIF+ZOXl08HrhLFdXh0qUnJmr2Mc9RygFYtm0CpA4MCouOGL79u3T5ZoADIaLvHEdE/UicGxA65qobwIgY2tO3RHDa7XcHGj/x90pJCpJ5C6imKw8aBUgcQ/RiRMnKCAgQOP+kSNH0urVq3W9NgCDBUl9fL3ErDb+43fzYQot+vc2nY58RDl5+WRrg2oYAHKTlpVLiw/dVlTcd7KrRP7e1Y19WWDqARLnHPFMNT8/P2ratClVrlxZZT8nbwOY23Cb9McvOzefNpy7Rw+Ts2j3lRjq36aWsS8PAAzst+MRlJCWTQ2qO9GgdsWPmIDl0ipAmjVrlri/e/euxv1WvMwxgJmyq2RNwzvVowX7w2jliUgESAAyk5SeQ0uPhIvHH/bxQS+yTGk1dtCmTRvKz88v9ta6dWv9XymAAb3aqR7Z2ljR+TuP6dK9RHz2ADKy7Gg4pWTmUlNPF+rfGj3IcqVVgDRlypQS98+ZM0fb67HYOki+vr6iLgOYB07C7NeqcJVu7kUCAHmIT82iFccjxOOPnvYha2uMiJgbbmu5zeW2VxdWBVy4SAdZWVn06NEjqlatGtnb2+t0MZZGqsFw/vx51EEyQ8FRiTRw0XGys7GmE1N6kbszfr4BLN2sv0Pp12MR1LqOG22bEIiUERm3wVpPzwkODqbevXuTi4uLmPLP9zy7LSQkRNtTApgUv7pVxC07L5/WndacbwcAluNBUgatPnVHPP746aYIjmTOWtvgqEuXLnTp0iXq1asXDRs2TNzzc96OIAksxWuBDcQ9/9HkKf8AYLl+OnhLzGLt2KAadWvibuzLAXMMkL744gtRxjs6Opp2795Na9euFff37t2jiRMn0ueff67/KwUwgr4ta5KHiz3FpmTRrisx+B4AWKi7Cem04WzhkiKfPIPeI9AyQAoNDRWJ2La2tirb+fnXX38t9gNY0pR/tgrJ2gAWa8GBm5SbX0DdfDyw+CxoHyDZ2NgUu49rICkvXgtg7jDlH8Cy3YpNoa0Xo8Xjj/v4GPtywERoFcnUrl2bli1bpnHfr7/+KvbDfzDN3/yn/D//pBYKpvwDWJ7v992k/AKip309qU3dKsa+HDDnaf779++nZ599VhSM7NSpk5jiz1P9z5w5IxK09+zZI5K25Q7T/C0HpvwDWJa8/AKx/uKFu49p7p4bYtueD7pRUy8XY18amPM0f57ev2XLFhEULVmyROQd8T0/5+2GDo4yMzNp8uTJYl04ruLdvn172r59e5lem5iYSOPGjSMfHx9q0aIFde3aVSzEqwlv79atmziOj3/77bfF68HyYco/gOXYfeUBdZlzkF5ZdkoRHDnYWlNEfKqxLw1MiNbJQgMGDKCIiAi6fv06HTt2TNyHh4dT//79ydBGjhxJ27Zto+PHj4tSA9OmTaOXXnqJ/v777xJfl5eXR3379qVr166J0gVXr16loUOHUlBQkIg4lfFz3j5kyBBxHB/PXzO/ns8Dlg9T/gEsIzgav+YCPUjKVNmemZMvtvN+AKZzNjX3pAQEBIh7ycKFCw326R4+fJg2bdpE06dPJ3f3wroVHKRxL9f7779PJY0grlmzhk6dOkXffvstOTk5iW1cpqB+/fr06aefqhz7ySefUIMGDcR+xsfzTD5+PZc5AMuHKf8A5j+sNmNHKJWUV8L7+TgA6/IMY+Xm5orHd+/eLfHGw22GsmHDBnHPvTvK+Dn3aJ07d67E13IF8I4dOxZ57aFDhyg2NlY8f/jwoQjE1N+DX8evX79+vR6/IjCHKf8rn6zVBADmg3OO1HuOlHFYxPv5OIBKZf0IOO+Ge1AOHDgg7nk6vyngoS5XV1dF75HE29tb3HPSeIcOHYp9bcOGDYt8Lfxa7nni4TruieJ7fi6dU8LlDPizQOVweU35X/TvLbpwN5Eu3Uuk1nUw4wXAXMSmZOr1OLBsZQ6QeLqcp6eneFy3bl2aOXOmxuM4kODhLkOJi4sTAZI6aRvvL+m1jRo1KvW10n1x73PjRmGSX3FSU1MpOTlZ8ZwX9cXCvuY95f+vi9EiuXNwuzpiW8eG1cgGq34DmDT+XdXncWBasrKyxE257TVIgMR5OhIeaho9enSxxx45ckSni7I03bt3V3nOSeSGDCJBv3w8C6cBHw2LFzdW082BpvX3pWdb1sTHDWCi+B8Z/l0tbpiNxxK83Ar/4QHzM3v2bJoxY4Zxk7QXLFigcXtYWJhIfF66dCkZCg+tKffOSKRtHh4eOr9WGr4r7tiS3oNx/lJSUpLiNmXKlFK/LjBNPMPl293Xi2yPScrEDBgAE8e9vFP6NtO4T0q04H900BtsnrhtVW5rue01eIDUo0cPjdtTUlJo8eLFYtq9ofj5+YkgJSEhQWU7J2gzLmZZ0msjIyOLzHTj13JeEtdUks7Bz6VzSvLz88XrS3oP5uzsLIbipBuG1yxvBoy0DTNgAExbdl7hb6v6iDj3HC0e0Ra9wGbM3t5epa3lttfgAVJxU+e5UiXXIuK6QobCdYkYJ48r4+ecX8RFIxnXKpJmpSm/loOrs2fPqmw/ePCgGBarUaOGeM65V1wgUv09+HUcFErXAJYNM2AAzN/qk5Hi/sM+PrRubGf6YZifuD82qReCI9AuB4mn73NvCUtLS6OjR48WCZT4+b179zQORVUU7s0aPHiwyOnhCt48HLZz507at28f/fXXX4oZau+88w4tX75cXDfXbWLc08XDgZMmTRKv4dpG3APGBTBXr16t8j7z5s0TVbZ//vlnca709HRRvbtz5840YsQIg329YDyYAQNg3kKiEinkXhLZ2VjTKx3rkbuzvbEvCSwhQPrtt99E8pMUcGgaZuMAiae+T506lQyJgxkOkAIDA0UXm62trVjyRLmqN/cCValSRWUmmo2NDe3atUsESDzcxq/jdeV4rTlev0UZ90RxDxIfywvh5eTkiM+Ak9f5PGD5MAMGwLytOXVH3D/XygvBEehvsdo7d+4o8nXGjh0remPUcYDBdYFq1Spc+VzusFit5eUg8fpNnJBdUMIMGO6qR5IngGlJTM+mTl8foKzcfNo83p/a1cdMNUt3QcfFasvcg8TLb/CN8SKt6lPXoeQaUo6OjjRhwgRxA/PEQQ/PcOH1mjgYUg+S+DlmwACYpo3n7ongyLemK7WtV9XYlwMViEd5+JaRkWH4JG1el6ysNZOAaPPmzRQaGorgyAJwnSOe6cI9Req6+3ggyRPABOXnF9Ca04XDayP965vMShBQMbgjgttcbnt1UeYeJE3i4+Pp9u3bKpUrGQ+/ffbZZzpdGIApB0l9fL3ErDZO3I5LyaKvdl6j0xEJ9Dgtm6pWtjP2JQKAkqO34ulOQjq5OFSiF/yQAgIVGCDxLLUxY8bQtm3btHk5gEUMt/l7VxePOS+Plx65ej+Z1p6+QxN7NTH25QGAhqn9vDSQk51O/QIgI1oNsfEsNe6i5KnxXGuI6wbxjatocyIUltEAOeHfhTe7NhSPV528Q1m5eca+JAB4IupROh24XlgDb0TnwjxagLLQKpTm8t0nTpwQiccuLi4qCdvPPvssDR8+XJvTApgtXsB2zq4bFJOcSduD79PL7esa+5IAgIjWnblLPFc7sHF18vbQrbIyyItWPUhc64iDI5abm6uyr2rVqhQXF6efqwMwE7Y21jQmsIF4/OuxiGKrzQOA4XBv7vqzUeLxSPQegSECJF62IzOzcDVkLqz477//KvZxvYGYmBhtTgtg1l7pUI+c7GzoekwKHQ2LN/blAMjerssxlJCWTV6uDtS7uafsPw8wQIDUsWNH6tu3r+gpevXVV+m5556jAQMG0MCBA8WaZf369dPmtBZdB8nX11fUZQDL5eZkS0OeDK0tPxZh7MsBkL3VTypnv9qpHlWy0aq5AzPEbS23udz2GjwH6fPPP6fTp0+TnZ0dvfnmm2KdtpUrV4rp/kOHDqW5c+fqdFGWhmsxaFPFE8zPG10a0u8nI+nIzTi6EZNCTb1cjH1JALJ09X4Snb/zmCpZW9GwDsgJlJMJT4oyS5W0taVVSM3LiXAg5ObmJvKRvvrqK7FILfcorVixQiRuA8hR3WpO9GxLL/F4+dFwY18OAMl93TX+fazhWrSwK0CFBEjFVcresWMHNW7cmHbv3q3NaQEswhtdGon7bcH3RSFJADCspIwc2nrxvniM5GwwaID0559/atweFBQkxv54xXsAuWpXvyq1rVeFsvPyafXJwv9iAcBwtly4Rxk5eeTj6UwdG2JRWtCOXrPWnJyc6Jlnniky9R9AbsZ2baTo5s/IRuFIAEPhEhtScvZI/wZYdw0qPkl71apV4sZu3bpFvXr10viDGR0dLab+A8jZ0y28qG41R4p6lEGbLtxDNz+AgZy4nUDhcWlU2c6GXnyqNj53MEwPEgdAJd1sbW2pZ8+etHr1au2vCMBC1mp7PbBw+ZEVxyLEauIAUHHy8gvo5O0EmrP7ung+8Kna5GyPdddAe1YFWpT8feqpp+jixYs6vK08SFMMedYfVx6Xph6CPKRl5ZL/7AOUnJlLy0a1pz6+KFQHUBF2X3lAM3aE0oOk/yZFuDvb0VcDW9KzLWviQ5eZRYsWiVtGRgZFRkaKAtbalNrRKgfp6NGj2rxM1nWQQkNDERzJTGX7SvRqp8LFMZdhyj9AhQVH49dcUAmOWEJqttjO+0FeJkyYINpcbnt1oVWAlJCQQD/++CMtWbJEZfv8+fPFRQFAodEB9UWhujMRj+jSvUR8LAB6HlbjniNNwyDSNt7PxwEYJEBauHAhTZ8+XXRdKUtJSRFT/U+cOKHNaQEsTk03R+rfppZ4PHvXddoWHC3yJPAHG0B3/I+Hes+RMg6LeD8fB1BeWmWw7du3TwyztWjRQmX71KlTxbpsn376qcoCtgBy1tzLhf4iEoER31hNNwea1t8X+REAOihrIVYUbAWD9SBZWVkVCY4k7du3p6SkJK0uBsDScP4D9xypi0nKRH4EgI5quDjo9TgAnQOkx48fi2n9muTl5YkcJQC5Q34EQMXiKtkezvbF7rd60luLatpgsACpY8eO9MYbb4hASdmjR4/ozTffpM6dO2t1MZZq0KBB5OvrK6YdgnwgPwKgYllbEbm72BUbHDEeyua6ZCAfixYtEm0ut70Gz0GaPXu2CILWrVtHDRs2FJWzOTiKiIggFxcXOnXqlE4XZWl4qqE2NRjAvCE/AqBi/XM5hq49SBEzRas62VFcapZinxfy/GRrwpOag1ItQoMGSN7e3qLwEidl7927Vyw94u7uTsOGDROz2+rXL6z9YgiJiYk0efJkOnjwoKjkzcHanDlzKCAgoNTXHj9+nH755Rc6ffq0eC2vIceBH39dHPgp69GjB8XGxpKdnep/K3369KG5c+fq/esC84f8CICKk5qVS7P+Liwr807PxvR+UBPRa8v/mPDvHg+roecIdKF1HfZ69erRypUryZg436lv374iaAkODhaL5XIJAi41cOzYsVIjRz7ulVdeEVXBudJ1XFwc9e/fX7yOt6kHev/884+oig1QFvwHmvMfOCG7uCosbo62yI8A0MKPB8IoJjlTrHn4Tg9vEQz5e1fHZwnGzUFSFhUVpVh2JD8/nwxpzZo1Yjjv22+/FcERmzhxoghsuNRAaaytrem7774TwRHz8PCgWbNmidyqX3/9tcKvHywb/8Hm/AdWXAZEUkYO7UKlX4ByuRGTQr8eixCPZwxoQQ62NvgEwXQCpA0bNlCTJk1Ej8rzzz8vto0YMYImTZqkz+sr9Ro454mTxtV7hg4dOiSGxErCgRAPySmrU6eOYh+ArngdqMUj2op8CGXcs9SlceF/ux+uD6Z/b5T8swoAhXgG9f+2XRGzRJ/29aRezbDGIZhQgMRJxzw0xT017777LlWuXFlsnzlzJoWEhIheGUPgYTXOFeK6TOo5UvxLdOnSpRJfb29fdHrojRs3xH2vXr2K7OOlVAIDA0V2POcqca5TVtZ/SYEAxQVJxyb1onVjO9MPw/zEPT9f9Xoner51TcrJK6Bxq8/T6XCUxwAozV8Xo0WukYOtNU190kMLUBG0nsW2detWka/Djhw5Iu4bN25Mf/75J/Xu3Zs+/vhjqmicM9SoUaMi211dXRX7y2vx4sXUrVs3euGFF1S2V6lShTw9PenAgQMisOKveciQIbR7927av38/2dgU38WbmppKycnJiuf8ek3BGViu4vIj5g/1o/TsPDp4PZbeWHVOBE+t6rgZ5RoBTB0PSX/9zzXx+N1eTahO1cLUCgDGHRbKnRbc9hq8B4nfVAqO1HEgwcnT5cVDYtwTVJbb9etFKxPrw88//0w3b94U5Qs4P0kZB4Sff/45OTg4iGvo3r27mO3G171x48YSz8vHurm5KW4cYAIwWxtr+nl4W+rUsJqYlTP6tzN0KzYFHw6ABt/tvUHxqdnUyKMyje1a9J9jkLfZs2ertLXc9hq8BykzM1PcOFhQl5aWJmoilRcPW3FgUha1a9cW91xaQLlnRiJt46Tr8uQzzZs3T/QQ1apVuLhoafz9/RXlArjEQXEOHz5Mfn5+iufoPQJlnGC6fHR7Gr78NF26lyTu17/lLxbZxJRlgEJXopNozak74vGsF1qSXSWd5xiBhZkyZQp99NFHKmk4ugRJWgVIPATFi9JyQKFcAPHu3btiFhkPsZVXjRo1SgwyNOGgg4MTzjdSzkMKDw8Xz1u3bl2m8/z+++/09ddfi94gLl+gLjs7W/SaqSd0S8Nqpc3ec3Z2Vgz7AWji4mBLK1/rSEOXnqSw2FQK+v6wSEKVYHFbkLP8/AL6YusV4l+J/m1qUWBjd2NfEpgg9fQVbnt1oVUIzsnJkZGR1KFDB3EBoaGhVL16dZEwzcNfhhpC4hwg7i06e/asynYuGslRIwddykGOpjXilixZIsoE/Pvvv4rgiItgvv3224pjTpw4oXH5lHPnzol7/hwAdFWtsh290aWwQKlycMSwuC3I2Z9noygkKpGc7SvRl/2aG/tyQCa0CpBq1qwpah99+eWXopeGA4umTZuKnBwOVpQDk4o0cuRIEbhwaYH09HRFkjUvecK9W8o4Z4qn8HNgpzwr7cMPPxTryvHQGtdV4tuWLVsUs9kkYWFh4twSDgR51h73YvGMPgBdcVD0w4EwjfukcGnGjtAiwROAJeKf85O3E+iP03fo//4prJj9YR8f8nQtmtoBYFKVtDkBasaMGeJmLDzEtWvXLhEgcaAiLTXCs8rUq2h7eXmJXi6pKCT3PEljlcpjlhLlcUs+14IFC0SOFC+Cx0uS5OTkiIXwpk2bhpwiMPjitqgYDJZs95UH4p8B5d8HXm/NyxWzf8EMAiSWkpJCly9fpujoaJE43apVK1G40ZB41tzSpUtLPW7VqlUqzzkniHOXyoK/pvfff1/cACoKFrcFKAyOxq+5UGR5ntz8Apr4x0VRMoNriwFUNK2nAfzf//0f1a1bl7p27SqSq/men3OyM6jiniaepce9TwDFweK2IHc8rMY9RyX964phZigNt7Xc5nLba/AeJM694QCJlxjhZT54WIsToE+fPi2G3HhWF+cnwX+Vx5Vn+wFos7gtz9PkJUv4OABLhGFm0IcJEyaI24ULF0pdtF7vAdIvv/wiVrbnNc/Ucf7PmDFjECABaLm4LQ8vcDCkHiTxc97PxwFYIgwzg9kPsfHUfk3BEeMaSIbOQwKw9MVtWf3qTvRMCy+jXBeAIWCYGUyJVj1IvCYZV8yWFqlVxgUVpUrXEp42/9JLL2l/lQAyC5L6+HqJ4Qb+j9q+kjW9u+4i3UlIp/N3HlP7BhhiA8vEw8ceLvYUl6J5EXAMM4PJ9yDNmjWLhg4dSidPnlQsBscz2rig4ptvvlmkBtFXX32ln6sFkNniti/41RYB06C2dcT2ZUfDjX1pABX6c1+/muYFaKWBZQwzg0n3IPXo0UPccw0idTx1vrTFWwGgfLjCNlcT3hv6kO4kpFH96kV7bwHMXXBUIp2781g8dne2EwvTSnjYmYMjTPEHkx9iGzduXJmO5YCJk7oBQHtNPF2oR1MPOnQjjlYci6AZL7TExwkWhduK/9tZWDGbe0y/HdxaMczMuUk8/IYJCmDyAVKTJk1EBemy4rXR5IxrMXAFb2nqIYA2xnZtJAKkDefuiSUXqjjZ4YMEi7Hn6kM6G/lY5Nx98oyPYpgZQJs6SHzLyMggg+cgaRpaU6Y+xHbkyBGSex0kXtAXwRHoIsC7OjWv6UoZOXm09vRdfJhgMbJz8+mbXdcU/wjUdCtcEgpAG9zWcpvLba/BAySuml0cXqcMRSIB9M/Kyore7NJQPF51IlI0KgCWgBekjUxIF3lH43p4G/tyALQPkDgy48VbNW3v3Lkz3bp1S5vTAkAp+repRZ6u9hSbkkU7Qu7j8wKzl5SRQz8cCBOPP+jtQ872Oi0RCmDcAKlRo0YUFhZGPXv2pDt37ohtc+fOpfbt21N8fDzVqFFDf1cIAAp2laxpdEADxZT/si64DGCqfv73Fj1Oz6HGNZxpWIe6xr4cAAWtQvWVK1eKNdj+/vtv6t69O7m7u9PFixfptddeo/nz59PSpUu1OS0AlMHwjvXppwO36HpMCp24nUCBjd3xuYFZinqUTr8djxSPP3+uGVWy0Xr9dAC90+qnUcoMDwkJofv379PVq1fJ39+ffvjhB7HMyCeffKLv6wSAJ9ycbGlIexSOBPM3d88Nys7LFxMQejbFyANYQIDEPUVcLPJ///ufWEIkKiqKhg8fTgEBAWLGWt++ffV/pQCg8HqXhmRlRWLaf9jDFHwyYJZFIbeH3Bc/x58/11xMQgAw+wApMjKSrly5QuvWraM///xTDLGNHz+e/vrrL5o6dSodPnxY/1dq5nWQfH19RV0GAH3gStrP+BYuXLv8aAQ+VDArnDv39c7Caf0vPVWHWtZ2M/YlgQVZtGiRaHO57dWFVYEWWZ5169al8+fPa0zGzs/PpxYtWtC1a4U//HJ24cIFateunfis2rZta+zLAQtzLvIRDV5yUiRuH5/USyzyCWAOdl+JoXFrzouikIc+7YG6R2CSbbBWPUgff/xxsTPVrK2tafbs2dqcFgDKoV39quRXt4qoh7T6VOFsUgBTh6KQYC60CpA++OCDEvcPHDhQ2+sBgDLinA2uOszWnLpDmTl5+OzAZOXlF9DJ2wk0ZcslURSyemVbFIUEy5jm36tXL3HPvUOdOnXSeAx3ZT1+/Fj84b59+7b+rhIANHqmhSfVqepI9x5n0Lw9N6hVHTcs7AkmZ/eVBzRjRyg9SMpUbMvJK6BjYXH0bMuaRr02AJ0DJJ6ptnz5crFQ7d27qutA1atXT9xv27ZNLDVS0lIkAKA/XDemc6PqtOn8PVp+7L9k7ZpuDjStvy8aHzCJ4Gj8mguknuyakpkrti8e0RY/p2DeAZKzs7MoCin1FLEHDx5QrVq1KDw8XDyvU6ewNoutrW3FXC0AFGl8Np+/V+RTiUnKROMDJjGsxj1HmmYC8Tae2M/7+/h6kY01pvmDmeYgKdeoiIiIELfmzZsrgiMoHqb5gzEaH8b7+TgAYzgT8UhlWE0d/2Tyfj4OwNSm+etU1x2Fvcpm8+bNYiHfCRMm6PJxA6hA4wOmLjYlU6/HAZQFt7Xc5nLbqwuzX/gmMTGRxo0bRz4+PqL+Euc/nThxosxrynl5eZGfn1+Rm7ScirIVK1ZQ69atxY2j02+//VbUfQIwBjQ+YOpquDjo9TgAk8xBys7OFonaynUlNW1jnKhtCHl5eWJZEzs7OwoODiYnJydauHAhBQUF0bFjxxS5UiXh4Gr69Oll6rLjNeZ4KZUOHTqIoUVefy4+Pl4ESgCGhsYHTF3HhtXI09WeHiZnadzPiRtebg7iOACzDZC4u6pBgwYq2zgwUt9mSGvWrKFTp06JGwdHbOLEiSJI+vTTT+ngwYN6eZ+UlBT6/PPPacyYMSI4Yo0aNaIPP/yQvvzySxFk8XMAQ+JGhWercUJ2cVlGvB+NDxgLJ177N6pOW4PvF9knZbXybEskaINZB0ienp4iECgNB00//vgjGcKGDRvIxcWFOnbsqLKde5AWL15MsbGxxVb8Lo/du3dTcnKyOK/6+0yZMkWMc3JABmBI3Khw48JTpbmx0RQkvdWtERofMBouXnrsVoJ47OZoS0kZOYp93HOEUhRgEQES5+pMmzatTMeuWrWKDIGH1Ro25FXNVaeHent7i0Dt0qVL1Lt37xLPcebMGerXr58oWcB5RwEBATRp0iSR06T8PtJ51d+HhYSE6PGrAig7LrLHdWTUi/DZ2VhTdl4+/XH6Lg3tUJec7Mr8qw6gN39djKb41CzRk3nw4x4UHJUocud4eJh7NtFzBKaszH81T548SeUZjjOEuLg4jUNbrq6uiv0lcXBwEPlSv/zyixgq5B6n119/XeQuHT16VCRrK59HOm953yc1NVX0QEns7e3FDUBfQRLXkeFZbVLj09C9Mg1YeIzCYlPpy7+u0HdD2mDWKRgUl5f45UhhGZg3ujQkRzsb8veuju8CVJisrCxxU257DTKLjYOJsnJ0dCz3hRw6dEj8AS/L7fr166QPw4YNo7179yryqHg4jme2cfL35MmTSV+4wKabm5vihsV8oUJyPbyr0wt+tcU9D1/89MpTxLX3tlyMpg3novChg0HtvRpDEfFpYmjtlY6Fqy0AVCRuW5XbWqm4tbZMpt+dp82vW7euTMfWrl1b3Lu7u6v0zEikbR4eHuW+Dj5n48aN6fjx4yrblM9b3vc5fPiwojeKofcIDKFTo+r0yTNN6dvdN2jqtqvUqnYV8q2l2gsKUBE4xWHJ4cL1OEf516fK9ibT1IAFmzJlCn300Ucq6TG6BEkm81PLvTfco1MeHHRwIMO/jMp5SDwFn59zvaKS8NBYtWrVyMbGRmU7P1eubyQFN3zep556SuV9WJs2bUpdpkV9eA7AEMZ186azEY/o3xtxNOGPC7R9YiC5OGApIKhYJ8MTKOReEtlXsqbRAcab6QzyYq+WvsJtr2wLRQ4ZMkT04pw9e1ZlO0/v56hReQYb12xKSCicTSHhKfunT59W2cbnCwsLU0znZ1xriWfLHThwQOVYfs7BlK7lzAEqirW1FX0/xI9quTmI4Y7Jmy8XqVsGoG9LDxf+8/hy+zrk7ox8SzBPZh0gjRw5kjp37ixmnaWnp4ttPL2f14mbN2+eyrH9+/cXi+lGRkaqbJ8xY4aoxs04uYvrKHEwNWvWLMUxHBzx2CbnJ50/f15s4/eYP3++qIWEGkhgyqpWtqOFw9uSrY0V7bz8gH4/ecfYlwQWLPR+Mh2+GSfy397qqjrzF8CcmMwQmza492bXrl0iQOJhMFtbWzFktn///iJVtLlMQfXq1VUSyJcsWUK///47BQYGinM9evRIDKHxDLZOnToVWduFX8vFInn4LicnRwRHqH8E5qBtvao0pW9zmvl3KH21M5Ra1XajrNx8TLkGvVt6pDD36LlWNale9cICvgDmyKoA/e0V5sKFCyJQ416ntm3bVtwbAZQB/6pzUcndV2PEf/f5SiNtXKcGRftAV1GP0qnHvENiiv/f73ahlrXd8KGC2bbBZj3EZi44R4ln6fF6bgDGwj2fT7fwFI+VgyPGy5WI4OnKA+NcHFiEX49FiOCoaxN3BEdgNNzWcpura34wAiQD4KVIuHgmD9MBGAs3XHP33NC4T4qXuCI3HwdQXo/SsunPs3fF47e7IfcIjIfbWm5zue3VBQIkAJngStvKy5Go47CI9/NxAOW16kQkZebkU8varhTYGBWzwfwhQAKQCV6GRJ/HAUjSs3Np1cnCGcLjuntjWRuwCAiQAGSC12jT53EAkvVnoygxPYfqVXOivi1r4oMBi4AACUAmePV0nq32X815Vbyd9/NxAGWVk5dPy49GiMdjuzUS6wICWAIESAAywQ0XT+VnVsXkIPF+NHBQFpzMf/J2As3cEUrRiRlUvbItvdyuDj48sBgIkABk5NmWNWnxiLbk5VZ0GK2Kky31aPrf8jwAxeFyEF3mHKRXlp2i1acKK7Nn5RbQoRux+NDAYiBAMgDUQQJTC5KOTepF68Z2ph+G+dGq1ztQTVd7kUOyGsuQQBmCI66ZpT4jMi0rF7W0wKLqIKGSdgVCJW0wFxvORdFnmy6JXqQjn/UkVwdbY18SmOiwGvccFVcugoduuXeSA3AM1YKxoZI2AOjspadqU+MazqIX6ZcnK7EDqEMtLZATDLEBAFWysaZPn2mqWC4CtZBAE9TSAjlBgAQAwtO+nvRUvSqUkZNHPx24hU8FikAtLZATBEgAoFjMdtKzzcTjdWfu0p2ENHwyoIJrZFV1Kj4/DbW0wJIgQAIAhc6NqlN3Hw/KzS+g7/bexCcDKuJSsig7N1/jpyLV1kItLbAUCJAAQIWUi7Q95D5dvZ+ETwcUM9g+XB9Madl5VLeqI3m52qt8Mjx7jWtscRkJAEtQydgXIAdci8HR0ZEmTJggbgCmrGVtNxrQppYIkObuuUErX+to7EsCE7D40C06GZ5ATnY2tOr1jlS/emUxq40Ttzk3iYffMLUfTKUOEt8yMjJ0Og/qIFUg1EECcxUZn0a9vz8shtr+fKuzGHoD+Tp/5xENWXpK9CLNe7kNDcaSImAGUAcJAPSugXtleqVjPfF4zu7rVFDAK7WBHCWl59B764JFcDTQrxYNalvb2JcEYBDIQQIAjd7t1ZgcbW3o4t1E2hv6EJ+SDHFgPHnLJbEYbf3qTjRrYEsx2xFADhAgAYBGNVwd6PUuDcRjzkXiHgSQlz/O3KVdV2KokrUV/TjsKXLBEjQgIwiQAKBYb3f3Fuuz3YpNpbm7r9O24Gg6eTsBwZIM3IhJoZk7QsXjz55tSm3qVjH2JQEYFGaxAUCxeNHaXs1q0JYL0bTkyH9rtNV0cxD1bjCl2zJlZOfRu+suUFZuPnXz8aA3uzQy9iUBGBx6kAw0zd/X11dMOwQwJ7uvPKC/LkQX2R6TlEnj11wQ+8H88fAp9wxKPYQz/r5KNx+mkruzPX33chuytkbeEZgPbmu5zeW2VxeY5l+BMM0fzL3R7DLnID1IytS43+pJccBjk3qh/o0Z4yB3xo5Qjd/n1W90pK5NPIxyXQDGboPNfogtMTGRJk+eTAcPHiRbW1uqVq0azZkzhwICAkp9bYMGDcjJyYns7OyKnDMqKooePnxI7u7uYluPHj0oNja2yLF9+vShuXPn6vmrAjA+LgJYXHDEOGWb9/Nx/t6ok2SuwRH3BBaXfp+WlWvgKwIwHWYdIOXl5VHfvn1F0BIcHCyCnYULF1JQUBAdO3ZMRI6l+eeff0SgpOzdd9+l27dvK4Kjko4FsFRcIVmfx4Hp9RByz1FxwRH3EPL+Pr5e6CEEWTLrHKQ1a9bQqVOn6NtvvxXBEZs4cSLVr1+fPv3001Jfv2DBAvLwUO0+Tk9Pp9WrV9P48eMr7LoBzAEvH6HP48B8ewgB5MisA6QNGzaQi4sLdeyoulYU9yAdOnRIDImVZODAgVS5cmWVbevXryc3Nzfq169fhVwzgLngtbV4tlpJ6bkezvbiODA/6CEEsOAAiYfVGjZsWKSyq7e3t6gAe+nSpXKfc+nSpfTWW2+RtXXRj2b+/PkUGBgosuM7d+4scp2ysrJ0+hoATBUvPMpT+ZlVCcM0j9KyDXpdoB/oIQSw4AApLi6OXF1di2yXtvH+8ggJCRFZ72+88UaRfVWqVCFPT086cOAAXb16VQRH33//PT377LMiF6okqamplJycrLghqAJzwXWOFo9oK2arKfN0tacaLvb0KD2b3l59jjJzSv4dANPDPX9ujrbF7uegmHsQ0UMI5oLbVuW2ltteiwiQeEiMe4LKcrt+/XqFXAP3Hr344ovk5eVVZN/WrVvp888/JwcHB3EN3bt3p6lTp4rr3rhxY4nn5WN52E66zZ49u0KuH6CigiSeyr9ubGf6YZifuD8xOYj+fKuzaGAv3E2kSZsvYUFbM/PP5QeUlJGjcZ/UY8g9iNyTCGAOuG1Vbmu57bWIWWw8bLVu3boyHVu7duFq0jzLjKNEddI29QTskqSlpdHatWtp27ZtZX6Nv7+/uD9+/DgNGzas2OMOHz5Mfn5+iuf29vZlfg8AU8CNpPpU/kYezrR4eFsateIMbQu+T409nOndoCZGu0You3+vx9KH64PF465N3CksNlUU/5RwjyEqpYO5mTJlCn300UcqaTi6BEkmEyDVqFGjxCBDEw46ODjhfCPlPKTw8HDxvHXr1mU+FwdnHHhxvSN12dnZoquOaywps7GxEff5+fklntvZ2VnjUCCAuQto7E4zX2hJn/91mb7bd1METf1a1zT2ZUEJToUn0Lg15yk3v4Be8KtF84f4iRlrPFuNE7c5N4mH1dBzBOaGOx+UOyC47bWIITZtDBkyRPQWnT17VmU7F43kqJGDLuUgJyEhocThtXHjxmncd+LECZGUre7cuXPivkOHDjp8FQDm7dVO9ej1wIbi8ccbg+nSvURjXxIUg783b646J9ZYC2pWg+Y9WUZE6iF8wa+2uEdwBGDmAdLIkSNF4DJp0iRRv4gtXryYIiIiaN68eSrH9u/fn+rUqUORkZFFzsOJ2aGhoTRq1Khi3yssLEycW8J5UDNnzhS9WK+88opevy4Ac/NFv+bUs6kHZebk09jfz1H04wyVtb14thsYV9jDFBq94gylZuWSf6PqtGh4W7K1MesmAKBCmcwQmzZ4iGvXrl0iQOJARVpqZP/+/UWqaHPidfXq1cnR0VFj7xEHOTxTTRM+FxeV5GE4XgQvNzeXcnJyxEJ406ZNQ04RyB73OPz4ylM0aPEJschp97n/iiEcCc+GQk6L4XBAqjxkxp//iF9P0+P0HGpTx42WjW5PDraFKQIAoBkWq61AWKwW5GbNqTv05dYrRbZLGYJcMoBnxYFhF5+1sbKivIIC8vF0pvVv+VPVyqprSgJYogs6LlaL/lUA0FuvxaJ/b2ncJ/UlccON4baKX3xWfQkRDo7YG4ENERwBlBECJAPgoTguY8DDcwCWCmt7mfbis2zBgTAEqGDxFi1aJNpcbnt1gQDJADZv3iySwCdMmGCItwMwCqztZdoBKsPisyAHEyZMEG0ut726QIAEAHqBtb2MCwEqgH4hQAIAveDigjxbqqSFKXgNN6ztVTHcne31GsgCyB0CJADQ21R/nsrPiguSnO0rUW4pleeh/B6nZdPiQ5oT5CVYfBagfBAggcF88cUXtGrVKvH41q1bYjkYrl/10ksviWVm+NazZ0+xvWnTpqL6eXGef/55sc4OmBaews9T+XktL2UezvbkYGtNt+PS6LNNWNhWn65EJ9HzPx2jY7cSyO5J4Uf1ABWLzwLIrFAklL1QnLHXVuJK5xs3bqRLly6J53v27KE//vhDpQo5r6nXu3dvsrOzE0U5+b44b7/9Ns2fP18RcIFpBUl9fL2K/PzxGmCjnyxsW796Zfqoj4+xL9Xsf583nIsSdaeyc/OpfnUnWjqyHUXGpxWpg4TFZwHKDwGSAfBUQ67gzZn1hpjJpqlQnLErGXMgw70+Dg4Oisrm6lMwv/32W7GO3nfffVdqUa/nnnuOPvzwQ4qJiRHnAtMire2lLLCxO/3fiy1p0ubL9OOBMKpfzYkGtatjtGs0F5p+n71c7amJpwsdDYsXz3ldte+H+pGboy0183LVGKBifTWQ0zT/RYsWUUZGhk7nQSVtC6ukLRWKKzChSsbcM9S8eXPatm2bGDrThBf+DQgIoF69eonlY3iYrTS8Fh4v+TJr1qwKuGqoKHN2X6fFh26TrY0VrX6jE3VupBpIQem/z8q4J25iz8Zi0VkA+A8qaVsgDijSs3PLfUvJzKFp269q/GMqbZu+PVQcV95z8zVp6++//6ZatWoVGxylpaXR8OHDxTp63NNUluCIvfnmm7R8+XKd/0sAw/r06abUr1VNyskroLdXn6fbcan4FmhZ+LFaZTuagOAIoEJgiM0EZeTkke/UPXo/L/+hjUnOpFbT95b7taEznyEnO+1+XL7//nt65513it3/3nvvUVhYmOg58vT0VGznYIkTte/cuUMDBgygjh07qryOgy5/f39avXo1vfXWW1pdGxge93R8N6QNRSdmUHBUIr2+8ixtGhdAt2JTMSRUzsKPj9KyxXHqw5kAoDsESFChLl68SNevX6eBAwdq3M+J2ytWrKCPPvqInnnmGcX2K1euiCG5LVu2UFZWFnXr1o1OnTpVpHeJk7U5F2ns2LFl7nkC4+OV5JePbk8DFx2nOwnpFPDNAdGjZCo5c6YAhR8BjAsBkglytLURPTblxf9JjvntbKnHrXytQ7mL9fE1aWPhwoU0ZswYsrW1LbIvKipKBDicnzV79myVfTzDjXOSmL29PVWvXp1OnDhBgYGBKsc9/fTTlJmZKZK7g4KCtLpGMF5hw9cDG9LMv0NVgiMWk5Qpcm+MkTNnKlCZHMC4UAfJBHFPCA9nlffWtYlHiZWMpUJxfFx5z61t7wwPg92/f7/I9vz8fBoxYoQYQlOf0s89RiEhIeK1yufhtXXUcf7R48ePVYbmwHxybJYdDde4TwqXOAeHj5Mj/iemhkvx1bFR+BGgYiFAkkklY2MViuOyBjxUxkGMMu4xOnLkCP3000/k4+NTpNeJA6fKlSsrtjk7O1N8fOGUZmV//vkndejQgVq2bFmBXwUYI8emQOaLq+bk5ZOzg+ZOfhR+BKh4CJAMgOv9+Pr6iroMxqpkzM+NMVzBNYpeeOEF+v333xXbTp8+TdOnTxeVs1977TWV4zloOnr0KFWpUkUUl5TwY57lpm7JkiUifwnMD3JsipefX0CfbrpE4XFpogK5u7OdSfw+A5gDbmu5zVWvtVdeyEEygM2bNxusDlJJlYyNVSiOk6h5Gv/7778vhs/4cW5uLiUlJYkgiXE9o+joaDp79qwoFMn7Hjx4oDgHD9ONGjWqSI0LLhGgnNwNlpdjs+bUHfL2cKaWtd1MtlK8vs3ff5N2hNynStZWtGJ0B+rUqLpFf70A+iQVZZbqIGkLAZKMKhkbC6+3VqNGDdEz1LVrV7EOW2l4SRLuZeLeIWmqP0/pV+89+uCDDzB7zUxxI885cZyQXVKW0dnIx2KtsR5NPURBxPjULJOrFK9Pm87fo58OFv6OfP1SKwpo7C4em8rvM4BcYIgNDNaLtHTp0jIf37p1a1FVe/HixfTVV1+JvCQbm/9m0qWkpND+/ftp5MiRFXTFYOycOb5Nfd6XBvrVIu4sOXQjjgYvOUnj1lwokrskzXrjytPm7OTtBJqypXC9wgk9vWlI+7rGviQA2UIPEhgEr8P2xRdf0KNHjzTmEmkyceLEYvdxcUieBSet7QbmScqZK21x1Q/7+NDPh27R+rP3NJ6n4ElAxefh4WVzHH7iQplvrz4nSh4837omfdxHc+V5ADAMBEhgENbW1jR58mQ6efIk9evXT+fzXbt2jb788ku9XBsYV1ly5upXr0wD/eoUGyCpz3oz9eEo9Rwqb4/KoqJ4cmYuta1Xhea93AZrqwEYGQIkMBhOztYXLg8A8sqZs5RZbzwMqN5jxgv3cs9R3WqOtGxUe1FpHACMCzlIFjbNH0Dus94i49N0Wly5ooOj8RpyqKRK4lxZvLpz8cUhAcBw0/wtIkC6d++eWHLCVNfi4mn+XAWapx0CgG6z3kr7LZ+/P4xeW3mW7iSkqQxpcQL0tuBocW+M6tz8ntxzVNI7/3IkXLaVwwH0hdtabnO57ZV1gLR27VqxZldkZKRWr+faO1yLh6s5N2vWjJ577jmR36LJjh07qH379mKGVdOmTUVODa8DBgCmMeutb0svsrOxFjPe+sw/Qgv236TtwdHUZc5BemXZKXr/z2Bxz8/1OeOtLAFYaZXDSeaVwwFMjVkHSFxMkKeOHzp0SLGwaXkkJydT9+7dRc8TR5scGPGSFbxy/N27d4sERy+++CJNmzZN1Og5fvy4WEJDvXghABivUvziEe1o1wddqUtjd8rOzacF+8PovT+Dy10WoDw9TnyOkgIwHu67+TCF1p6+YxE5VAByYVVgqoP1ZcALnkozpHjF+FWrVpUr92DGjBk0c+ZMsaq8tDAq9wjx8hgvvfQSrVixQmzjc3p7e4teo127dqmsA/bKK6/Q4cOHRVClTqrief78eYNW0gawdKVV0ubf2e3B9+nDDcFUXGxj9SSwOjapl8prNSVRF1eIUsopUn8LPhtv69Wshpi+f/fRf8vmlGbd2M4mPwsPwBzo2gabdQ8SB0Z809aGDRuoefPmKqvGc12dwMBA2rRpkyIAO3fuHEVERFBQUJDK66Xn69ev1/oaAED7WW8v+NUW9+p1j7hXuIarQ7HBkXJZAB6Cy83LLzGJWlOPU0k5RdK2g9djRXDEw37dfdzJtZjFZ8U1PwnEONgDAOOT7TR/7im6fv26KGCojnuL/vnnHwoPD6fGjRtTcHCwYrsyDw8PcnFxoZCQEINdNwCQXoeqPtwQQp9tvkR1qzpSdKLmZU+kbZ9uvERHw+IpLiWLbsellppTJM7f24fe7NqQKttXUgRgyudkUnjHvVTmWOQSwBLJNkDiis7cQ+Tq6lpkn7QtLi5OBEh8r7xd/Vhpf3FSU1NFvpPE3t5e3ADA+GUBpBpE4fGlD4OlZOXS2tOq+YmlaeDuJIKj8lQOB4Dy48XQ+abc9lpEgMSJ1j179izTsZxMzTPOzAUngivjRG9eiBUAjLcYrpSDdOTTnhSbmkVrT92hnw/dLvW8z/h6UhcfD0pMy6bv9t0sd6BWlsrhAFB+s2fPFrnF+mIyARIXdVq3bl2Zjq1du7bO78frgXH+knLPjkTaxkNozN29cDXt4o5t0KBBie/FSdy8or0EvUcAhisLwENaUtK0piEt20rWVLuKI3Vt4lGmAGlMYEOR98Q5SH+cuVtqAKYpp6gslcMBoHymTJlCH330keI5p8eod1CYZYBUo0YNUY/IUDgZm2elcZ6ROt7GuUWNGjUSz6XgRv3Y+Ph4sap8mzZtSnwvZ2dnjcNzAFCxyjOkVdYeJyngKWsAhp4hAMNQT1/htlcXZj2LrTyys7MpISFBZduQIUNE/aMHD/6bmcLjl1zjiEuUSzPkOnToQA0bNqQDBw6ovF56zucp7j2V7wHA8DgI4qn8PH3+h2F+4p6fq+f7lFaIUlPAU1pdJuQUARiPrm2wyfQgVbT+/fvTkSNHRP6SNCTGXXGrV6+mTz75RNRQsrGxEflBfK88jslThn/44QdRKHLnzp1iNXoOtviYwYMHF9uFhwAJwDSUdUhLmyRq5BQBmCZd22Cz70EaN26cGALbvn27eM6P+RYbG6tyHBd/rF69Ojk6Oiq28bAX5wfl5uaKHCiuicRT9jmQqlevXpEAa8uWLTR16lSx1Ii/v7/YxgEW/EduC/Li67XcHqdfh7eiuO3fintNPU7lqctkLuT28yzHr1luX69sK2mbOg60uHepuErblogDTR62lAt8vZaLJ2C4ubmJJY3kkkMot59nOX7Ncvp6j+jYBstmiM0YpIVsb968qXOymLnIyMgQ5d3lAl+v5ZJqqPBMGPz+Wi78DluumzcLy3Bou6g8epAqEOc18RpxAAAAYBwrV66k0aNHl/t1CJAqEJcB2LNnj0gKV859AgAAgIrvHYyMjKRnnnlGUc+wPBAgAQAAAFjaLDYAAAAAfUOABAAAAKAGs9gAANRw9ZP58+fTw4cPydbWlmJiYsRzXoIIAOQBOUgAABpWBd+/f79iOSEuEMvT/aWCtABg+RAgAQAoefToEdWtW5eWLFlCI0eOFNt4Jgyvx8iF57p27YrPC0AGkIMEAKBk9+7dlJ6eTu3atVNs41IdvFTRhg0b8FkByAQCJAOaN2+e+C+UF8PlBXClG6/rBgCm4fLly+K+Tp06Ktv5OQ+zAYA82lokaRsIL2rLyZ5OTk70zTffUOXKlWnBggU0c+ZM8vb2NtRlAEAZhtgY/44q4+VG1BfBBgDLbWvRg2QgQUFBIp8hMTGRPv74Y3JwcKB+/frRsGHDqEOHDoa6DAAohZ2dnbjn/ziV8XNpHwBYfluLAEkH3OWu3H2n6cZrwLBatWqJ/AVeUdja2po2btwovpEAYFpq1Kgh7tPS0lS283P+PQYA06XPthZDbDo4deoU5ebmlniM8vov3PU3btw4evz4MR0+fBgJnwAmyM/PT9xHR0dTs2bNFNv5ec+ePY14ZQBQFvpqaxEg6UA9ibMkx44do/Pnz9Pzzz9Pu3btoqZNm4qic5cuXaLWrVvrchkAoEe9e/cW+QsXLlxQBEhRUVEi/2jw4MH4rAFMmD7bWgyxGcC6devEasJDhw6lqlWr0vXr1yk8PJzeffdd8YcXAEyHo6MjffDBB7Rq1SrFNh4q79KlC3Xu3Nmo1wYAhmtrESAZQGhoKNWvX19k1LOAgACRn8R/iDl5DABMy6xZs0RC53vvvUfTpk2jsLAw2rZtW5HEbQCw3LYWlbQBAAAA1KAHCQAAAEANAiQAAAAANQiQAAAAANQgQAIAAABQgwAJAAAAQA0CJAAAAAA1CJAAAAAA1CBAAgAAAFCDAAkAAABADQIkAAAAADWV1DcAAABRfn4+ff/991SpUiU6d+4cBQUF0WuvvYaPBkAmsBYbAICagoICGj58uAiK3njjDTp69Cg999xzFBMTQ5UrV8bnBSADGGIDAFAze/ZsysnJEcER416k1NRUioiIwGcFIBMYYgMAUBIeHk5fffUVXb58WbHt4cOH4t7KygqfFYBMoAcJAEDJ119/TX369CFvb2/FtrNnz4p7Dw8PfFYAMoEcJACAJzIyMsjd3Z0CAgKoRYsWis9l8+bNlJeXR/fv38dnBSATGGIDAHji+PHjlJ6eTgsXLqSmTZsqgqalS5fSSy+9hM8JQEYwxAYA8MSVK1fELDUpOGL79u2jzMxMGjVqFD4nABlBgAQA8ERcXFyRPKPly5eL4bann34anxOAjGCIDQDgCXt7e5WZardv36Z//vmH9u7dixlsADKDHiQAgCc6duxIDx48EDWQuFjke++9R1OmTKFevXrhMwKQGcxiAwB4goOicePGkYuLC1lbW1OzZs3o9ddfx+cDIEMIkAAAAADUYIgNAAAAQA0CJAAAAAA1CJAAAAAA1CBAAgAAAFCDAAkAAABADQIkAAAAADUIkAAAAADUIEACAAAAUIMACQAAAEANAiQAAAAANQiQAAAAANQgQAIAAABQgwAJAAAAgFT9P1+/+rNao8iMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### Here run the code and complete it to get all expectation values.   \n",
    "\n",
    "Nshot = 1000\n",
    "n_theta = 41  # number of theta points\n",
    "theta_list = np.linspace(-np.pi, np.pi, n_theta)\n",
    "\n",
    "sampler = StatevectorSampler()\n",
    "\n",
    "# -----------------------\n",
    "# Helpers: counts -> expectation values\n",
    "# -----------------------\n",
    "def bit_to_ev(b: int) -> int:\n",
    "    # Z-eigenvalue mapping\n",
    "    return 1 if b == 0 else -1\n",
    "\n",
    "def exp_Z0(counts: dict, shots: int) -> float:\n",
    "    # Qiskit bitstrings are 'b1b0' (bit 1 left, bit 0 right)\n",
    "    return sum(bit_to_ev(int(s[1])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_Z1(counts: dict, shots: int) -> float:\n",
    "    return sum(bit_to_ev(int(s[0])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_two_qubit_product(counts: dict, shots: int) -> float:\n",
    "    # After basis rotations, Z-measurement corresponds to X or Y measurement.\n",
    "    total = 0.0\n",
    "    for s, c in counts.items():\n",
    "        b1 = int(s[0])\n",
    "        b0 = int(s[1])\n",
    "        total += (bit_to_ev(b0) * bit_to_ev(b1)) * c\n",
    "    return total / shots\n",
    "\n",
    "# -----------------------\n",
    "# Circuit builders (independent, no compose)\n",
    "# -----------------------\n",
    "def circuit_z(theta: float) -> QuantumCircuit:\n",
    "    qc = QuantumCircuit(2, 2, name=\"Measure_Z0Z1\")\n",
    "    qc.ry(theta, 1)\n",
    "    qc.x(0)\n",
    "    qc.cx(1, 0)\n",
    "    qc.measure([0, 1], [0, 1])\n",
    "    return qc\n",
    "\n",
    "# -----------------------\n",
    "# Loop over theta\n",
    "# -----------------------\n",
    "Z0_list = []\n",
    "\n",
    "for theta in theta_list:\n",
    "    counts_z  = sampler.run([circuit_z(theta)],  shots=Nshot).result()[0].data.c.get_counts()\n",
    "\n",
    "    Z0_list.append(exp_Z0(counts_z, Nshot))\n",
    "\n",
    "Z0_list = np.array(Z0_list)\n",
    "\n",
    "# -----------------------\n",
    "# APS-like plotting style (no manual colors)\n",
    "# -----------------------\n",
    "plt.rcParams.update({\n",
    "    \"font.family\": \"serif\",\n",
    "    \"font.size\": 12,\n",
    "    \"mathtext.fontset\": \"stix\",\n",
    "    \"axes.linewidth\": 1.0,\n",
    "    \"xtick.direction\": \"in\",\n",
    "    \"ytick.direction\": \"in\",\n",
    "    \"xtick.minor.visible\": True,\n",
    "    \"ytick.minor.visible\": True,\n",
    "    \"xtick.top\": True,\n",
    "    \"ytick.right\": True,\n",
    "    \"legend.frameon\": False,\n",
    "})\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6.0, 4.0))\n",
    "ax.plot(theta_list, Z0_list, marker ='o',label=r\"$\\langle Z_0\\rangle$\")\n",
    "\n",
    "ax.set_xlabel(r\"$\\theta$\")\n",
    "ax.set_ylabel(\"Expectation value\")\n",
    "ax.set_xlim(-np.pi, np.pi)\n",
    "ax.set_ylim(-1.05, 1.05)\n",
    "\n",
    "# Nice tick labels at -pi, 0, +pi\n",
    "ax.set_xticks([-np.pi, 0, np.pi])\n",
    "ax.set_xticklabels([r\"$-\\pi$\", r\"$0$\", r\"$\\pi$\"])\n",
    "\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "970f568b-3ee4-417a-b34d-115b78f04765",
   "metadata": {},
   "source": [
    "## 4. Expectation value of $H_2$ from measurements as a function of $\\theta$\n",
    "\n",
    "## Task \n",
    "From the different expectation values above, compute the expectation value of $\\langle H_2 \\rangle$ \n",
    "as a list and display the result as a function of $\\theta$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e1fe9799-f4b7-4102-a9ea-b71d8f609454",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here show the result \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c1c50a-56b6-4cf1-8a8f-e817f0d02529",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}