ham_sim_tutorial_L12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trotterized Hamiltonian simulation with `qiskit` and `qiskit-ibm-runtime`\n",
    "\n",
    "In this notebook, you will learn how to use a  **Multi-Product Formula (MPF)** to achieve a lower Trotter error on our observable compared to the one incurred by the deepest Trotter circuit that we will actually execute.\n",
    "You will do so by working through the steps of a **Qiskit pattern**:\n",
    "\n",
    "- **Step 1: Map to quantum circuits**\n",
    "    - Initialize our problem's Hamiltonian\n",
    "    - Generate a Trotter circuit\n",
    "- **Step 2: Optimize the problem**\n",
    "    - Transpile our circuits\n",
    "- **Step 3: Execute experiments**\n",
    "    - Use the `Estimator` primitive from qiskit runtime to compute our circuit's expectation values\n",
    "        - Run local simulation\n",
    "- **Step 4: Postprocess**\n",
    "    - Inspect our results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Map to quantum problem\n",
    "\n",
    "### 1a: Setting up our Hamiltonian\n",
    "\n",
    "We use the XXZ model on a line of 12 sites:\n",
    "\n",
    "$$\n",
    "\\hat{\\mathcal{H}}_{XXZ} = \\sum_{i=1}^{L-1} J_{i,(i+1)}\\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\\cdot Z_i Z_{(i+1)} \\right) \\, ,\n",
    "$$\n",
    "\n",
    "where $J_{i,(i+1)}$ is a random coefficient corresponding to edge $(i, i+1)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAB4sAAAA7CAYAAABi6GwJAAAn3ElEQVR4Ae2dCbwVVR3H/+/xWB6rICKyiAiyIwZJYMgSirmRoBYhIoghWCofl4zwk1qKwUeDLAgpRCWEHiYfRDLBREGKpcQQqFiULTZjkR0ETuc/Ou/euct7971ZzjK/+Xzwzpx758z/fP9f5907Z+acPCEXwgICIAACIAACIAACIAACIAACIAACIAACIAACIAACIAACIAACIAACIAACIBAnAnPy49RatBUEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQOALAugshgkgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgEEMC6CyOYdLRZBAAARAAARAAARAAARAAARAAARAAARAAARAAARAAARAAARAAARAAARAo0AnBmTNnaM+ePc6/gwcPEm8fPnyYTp8+TVWrVqXKlStTYWEhnXPOOXTBBRdQnTp1dAofsVhOAH5anmDDmwc/DU+g5eHDT8sTbHjz4KfhCbQ8fPhpeYINbx78NDyBlocPPy1PsOHNg5+GJ9Dy8OGn5Qk2vHnw0/AEWh6+6X4q6Sw+fvw4rVq1itasWUNr166ldevW0ebNm2nv3r1OB3GuzlSpUoUaNWpELVu2pHbt2lHbtm2pY8eO1KZNG8rLy8u1GnwOBDwE4KcHBzY0IwA/NUsIwvEQgJ8eHNjQjAD81CwhCMdDAH56cGBDMwLwU7OEIBwPAfjpwYENzQjAT80SgnA8BOCnBwc2NCMAPzVLCMLxELDVzzwhF09LQ9g4e/YsrVixgv70pz/Ru+++63QUnzx50nky2O3k5Q5fflq4QYMGdP755zvv5efnU40aNaigoICOHTtGvM+JEydo//79tHPnTtq1axdt376d1q9f73Q4/+tf/6JTp07ReeedR1deeSX16tWL+vbtSxdeeGEIrUKVthCAn7Zk0s52wE8782pLq+CnLZm0sx3w08682tIq+GlLJu1sB/y0M6+2tAp+2pJJO9sBP+3Mqy2tgp+2ZNLOdsBPO/NqS6ti4uecUDuLly1bRjNnzqR58+Y5nbvNmjWjnj17Uo8ePZx/QXfi8nDVH374IS1ZsoTee+8959+hQ4eoU6dOdPPNN9PgwYOdzmhbJEU7/BGAn/74Ye9wCcDPcPmidn8E4Kc/ftg7XALwM1y+qN0fAfjpjx/2DpcA/AyXL2r3RwB++uOHvcMlAD/D5Yva/RGAn/74Ye9wCcDPcPmidn8EYubnHOIni4Nc5BzDYtKkSaJ9+/b8xLK49NJLxRNPPCHkkNNBHianuuSTyOLNN98Uw4cPF3Xr1hXyCWXRr18/sXDhwpz2x4fsIwA/7cupTS2CnzZl0762wE/7cmpTi+CnTdm0ry3w076c2tQi+GlTNu1rC/y0L6c2tQh+2pRN+9oCP+3LqU0tgp82ZdO+tsTYz6LAOosZ4sSJE0X9+vWFnEtY3HrrrWLRokXa2MIdx0VFReKqq64Scj5jpxObt+Uj5NrEiEDCIwA/w2OLmv0TgJ/+GaKG8AjAz/DYomb/BOCnf4aoITwC8DM8tqjZPwH46Z8hagiPAPwMjy1q9k8AfvpniBrCIwA/w2OLmv0TgJ/+GaKG8AjAT+G/s/jMmTNi8uTJ4txzzxU1a9YUjz76qNi3b194WQug5n/84x9CzmXsdBp37txZrFy5MoBaUYWOBOCnjllBTC4B+OmSwKuOBOCnjllBTC4B+OmSwKuOBOCnjllBTC4B+OmSwKuOBOCnjllBTC4B+OmSwKuOBOCnjllBTC4B+OmSwKuOBOBncVb8dRZzp+vll18uKlasKB566CHtO4mLm/3lygcffCDkHMoiPz9fjBgxQhw4cCD1I9g2mAD8NDh5MQgdfsYgyQY3EX4anLwYhA4/Y5Bkg5sIPw1OXgxCh58xSLLBTYSfBicvBqHDzxgk2eAmwk+DkxeD0OFnDJJscBPhpyd55ess5qGbx40b53QSd+/eXaxdu9ZTq0kb3JYZM2Y4w2c3adJEyEmrTQofsWYgAD8zQEGRNgTgpzapQCAZCMDPDFBQpA0B+KlNKhBIBgLwMwMUFGlDAH5qkwoEkoEA/MwABUXaEICf2qQCgWQgAD8zQEGRNgTgpzapQCAZCMDPDFBEOYah3r9/v+jTp4/TUTx+/Hhr5vz99NNPxfXXXy8KCgoEtwuLmQTgp5l5i0vU8DMumTaznfDTzLzFJWr4GZdMm9lO+Glm3uISNfyMS6bNbCf8NDNvcYkafsYl02a2E36ambe4RA0/45JpM9sJP7PmrWxPFm/btk20bdtWNG7c2Mp5fvmOgmeffVZUqFBBjBw5Upw+fTorObyhHwH4qV9OEFGCAPxMsMCafgTgp345QUQJAvAzwQJr+hGAn/rlBBElCMDPBAus6UcAfuqXE0SUIAA/Eyywph8B+KlfThBRggD8TLDAmn4E4GeJOcm9s3jjxo2iUaNGon379mLHjh0l1mr6m3PnzhWFhYXilltuEZ9//rnpzYlF/PAzFmk2tpHw09jUxSJw+BmLNBvbSPhpbOpiETj8jEWajW0k/DQ2dbEIHH7GIs3GNhJ+Gpu6WAQOP2ORZmMbCT+NTV0sAoefpaY5t87inTt3iosvvlh07txZHDhwoNRabfjAkiVLRNWqVcWdd95pzVDbNuQlUxvg59lMWFCmCQH4CT81UTFjGPATfmYUQ5NC+Ak/NVExYxjwE35mFEOTQvgJPzVRMWMY8BN+ZhRDk0L4CT81UTFjGPATfmYUQ5NC+Ak/NVExYxjwMyc/S+8sPnr0qOjQoYNo1aqV4Hl947QsWLDAmZv5sccei1OzjWor/Kwo4Ke+ysJP+KmvnULAT/gJP/UkgO+feuYlOSqcP3H+TPZBt3X4CT91czI5HvgJP5N90G0dfsJP3ZxMjgd+ws9kH3Rbh5/wUzcnk+OBnzn7WXpn8YgRI0Tt2rXFli1bkhnHZn3KlCkiPz9fLF68ODZtNqmh8BN+6uwr/ISf8FNfAvj7rm9uODKcP3H+1NlQ+Ak/4ae+BPD3Xd/ccGQ4f+L8qbOh8BN+wk99CeDvu7654chw/sT5U2dD4WfOfpbcWfz666+LvLw88eqrr+qc79Bj47mLeb7muAzBHTrQgA4AP78ACT8DEirgauAn/AxYqUCrg5/wM1ChAq4MfsLPgJUKtDr4CT8DFSrgyuAn/AxYqUCrg5/wM1ChAq4MfsLPgJUKtDr4CT8DFSrgyuAn/AxYqUCrg59l8rMojz9OGZZTp05RmzZtqEuXLvT73/8+wyfiUyQ7ially5Y0ZMgQGj9+fHwarnFL4WciOfAzwUKXNfiZyAT8TLDQZQ1+JjIBPxMsdFmDn4lMwM8EC13W4GciE/AzwUKXNfiZyAT8TLDQZQ1+JjIBPxMsdFmDn4lMwM8EC13W4GciE/AzwUKXNfiZyAT8TLDQZQ1+JjKRo59zKFtX/YQJE0RhYaHYtm1bto/Eqvy5554TVapUie1w3LolG356MwI/vTxUb8FPbwbgp5eH6i346c0A/PTyUL0FP70ZgJ9eHqq34Kc3A/DTy0P1Fvz0ZgB+enmo3oKf3gzATy8P1Vvw05sB+OnloXoLfnozAD+9PFRvwU9vBuCnl4fqLfjpzUAOfhblJ/qXE2tnz54lCZNGjhxJjRs3TrwR47W7776b6tWrR5MmTYoxBT2aDj/T8wA/05moKoGf6eThZzoTVSXwM508/ExnoqoEfqaTh5/pTFSVwM908vAznYmqEviZTh5+pjNRVQI/08nDz3QmqkrgZzp5+JnORFUJ/EwnDz/TmagqgZ/p5OFnOhNVJfAznXwufmbsLH777bdJPlFM3/ve99JrjbBk9erV9Oijj9LUqVPpyJEjER45/VCVKlWioUOH0ssvv0yff/55+gdQkhMBZrd3796cPpvtQ7r4yfEtXLiQFixYkC3UyMrhZzCobfFz69atNHbsWBo1ahS9+OKLys9Z8BN+JhM4evQozZ49mx588EGaNWsWj3CS/Hbk6/AzGOS2nD+TafB34WHDhtHp06eTiyNdh5/B4LbJzw0bNjhTBPE0QfzvlVdeoTNnzgQDqoy1wM8yAsvycZv8dJu4a9cux8+nn36aNm/e7BZH+go/g8Ftg5887J57zkx9XbNmTTCgylgL/CwjsCwft8FPbhr/Ppo8eTI99NBD9LOf/Yz+85//ZGlxNMXwMxjOtvh5/PhxmjlzpvP7/Te/+Q199tlnwQAqZy3ws5zgUnYz3c/Srser6lOCnymilXPTdj8ZS2kOlxNdibvl5Kf3YeQvtmSnqOjatWumtyIre+GFF8S1117rDPv80ksviU6dOolPP/00suNnOtCWLVtEXl6eeOuttzK9jbIcCOzfv18UFBSI3r17C9mJJeSXjBz28n5EBz8XLVok+vTpwz0c4vHHH/cGqGgLfvoHb4Of69atE9WqVRONGjUSFStWdBzt2LGjOHz4sH9APmqAnz7gfbmrDX7Ki8eiRYsW4rrrrhO1atVy/Lz33nv9w/FZA/z0CVDuboOfyRRk55vo2bOn4+iJEyeS34p8HX76R26Tnz169HC85O+g/I/PpyoX+Omfvk1+Mo3nn39eXHHFFWL58uVC3tHvH5CPGuCnD3hf7mqDn/KmBc950z1/8uucOXP8QypnDfCznOCSdrPBT/kwhWjWrJmYMWOGOHbsmHPubN26tZg7d25SS6NfhZ/+mdvg5//+9z9x5ZVXiokTJ4qdO3eKhx9+2Pk9z36oXOCnf/qm+pnL9XjVfUrwE37yd8xs/UW5OOyfYPYaSvEz8zDUS5cupauvvlq2S82yfv16uv/++2natGnUpEkTGjx4MJ177rn0k5/8RE1AXx6VY2nevDm9//77SuMw/eD8hM7ixYvpzjvvpLp161K/fv3otddeI3kxNqemqfaTg+zWrRvJCyE5xRvVh+BnMKRN95PPm/z0/fbt2+mTTz6h73znO/TBBx/QU089FQygctYCP8sJLmU30/2cMmUK/f3vf3dGZOCnji699FLnb/2hQ4dSWhrtJvwMhrfpfiZT+MUvfkHyJsXkImXr8DMY9Db4uWTJEmrfvj3xnfLuPx5BROUCP4Ohb4Of8pIE3XTTTfSHP/yB/vKXv9DXvvY1kjdaBwOonLXAz3KCS9nNZD/Zy3nz5lFRURHxdSb+fcT/Vq5cSdWrVyd5w01Ka6PbhJ/BsDbZTyYwbtw4uuSSS2jQoEFUWFjonDtvv/12kp1ywQAqZy3ws5zgUnYz3U8eDaxChQrONfoLLriAfv7zn1OVKlXovvvuS2lptJvwMxjeJvpZ2vV4HfqU4Cf8LIlAaQ6XtG8Q75XmZ9ow1PKuIWeoKPlkcRDHL1cdPPQKf1niP0Tu8o1vfMO5oMwdICoXeZc0/e1vf1MZgtHH5h9rvPC48fyPhxV444036JZbbnFuCOAvyPPnz886bK4OfnL8/OWoYcOGvKrVAj/9pcN0Pw8ePOjcyNClSxcHBDvKPz75Qt2KFSv8wQlgb/jpD6LpfnLrR48eTTVq1HBA8MUQvhmM/eShUFQv8NNfBmzw0yXAQ1LyTTYDBw50i5S/wk9/KbDFTx7S98c//jFddtllxf/OO+88f3AC2Bt++oNoi5/PPPMMyaeJnaEq+beSLgv89JcJ0/3k6Xn4ZsVbb72V5NOadNFFFzn/+Iab66+/nqpWreoPkM+94ac/gKb7ya3fsWMH8U20blu4TI4U5lxz4nWVC/z0R9/NqanXP3m4aR66/5prrikGkZ+f70zR+Prrr9M///nP4nIVK/DTH3VT/SzterwufUrwE35mI1Caw9n2C7K8JD/TOov5yzSfMFq2bBlkDGWqiy/QyWEqPfvwl/pTp06RfFTbUx71BsclH9eO+rBWH4/vZGLn5JA7zp3offv2dZ445km3+Slu9w8YQ9DBTzcZfHedbgv8DD4jJvl5zjnnOE/qJ1PgO4batm3r3ICTXK5iHX4GT90kP7n1lStX9kDgJzd5bm0dLirDT09qAtkwzU9u9MmTJ50nOX71q18pfyIuOQnwM5lGMOum+bls2TL685//TK1ataKbb76ZVq1aFQyIAGqBnwFATKnCND/59/uYMWOc+Qzr16+f0hq1m/AzeP4m+cnXkTp06JAGQQ4/Td/+9rfTyqMugJ/BEzfJT249PxjDnW7uSIocP88Py7+RVC/wM/gMmOTnhg0bSE7N4zzxnkyCb7zhhW+6UbnAz+Dpm+JnSdfjdelTgp/wsyQCJTlc0n5BvVeSnwWpB+EnN3nhYZ9VLHz8PXv2pB2fv+TzwkMGqVyYC8fIPy6wlJ3A0aNHS9yJ/zDxwkOSTp8+naZOneo8wXvHHXfQd7/7XYc9v6/KTz62u7jDqrmvbrnKV/jpj75Nfrok+C5WPm8+8cQTbpGyV/jpD71tfvJw1PwD9I9//KM/MAHtDT/9gbTFT35qk+9G1uF7RnJG4GcyjbKv2+CnnFeMBgwYQB999BHJeQyJn+jgJ43ZV9UL/PSXARv8nDBhgnODb9OmTWnIkCHOzdWdOnVyOj9q1arlD5DPveGnP4A2+JlKgK/n8PfQa6+9NvWtyLfhpz/kNvg5dOhQmj17Nj355JPOVFJ8femuu+6iYcOG+YMTwN7w0x9E0/2sWbOmA4B/sycvPEIYL9u2bUsujnwdfvpDbrKf7nV499UloVOfEvx0s1K+Vxv9TCbhuuu+Jr8XxTr7uW/fvsyHkk9tehZ58UHITwr5FK+nPKqNd955xzm+vKvOc8jNmzc75bLT0FMe9YYcgkNUrFjRiYU54V90DOTda+JHP/qRw1yVn8m+yU44JxbZCZdcrHQdfkbnY+r/+7r56YrI53Q5Z5xgX1Uv8BN+soOHDx8WI0eOFPJHpnMOlXfNC/k0p2o9BfyEn3K+d/HAAw8Uuzh27FjH0RMnThSXqVqBn/Az2b0FCxYI+QPP8XPhwoXJbylZh5/wU96dLho0aCDkfMXO33l5M4Pzd14+0SnktENKvHQPCj/hp+uC+ypvSBfy5ht3U+kr/ISfLKAcZU/06tXL+bvesWNHIR+gUeqle3D4GW8/+e/3hRdeKOQIduLAgQOuFs7fer4eJof4Ly5TsQI/4+tntuvxOvUpwU/4WVJ/UTaHozqXsp9yOsBMhytKG4a6du3a8pxPJP8QOK9R/0dG6RxSdsh6Dn38+HFnW/WwVtzrzr3vHCf+lZ1B1rsWPNmm4vkreR62++67j5YuXUo8Sf03v/lN55Oq/EwJU7tN+Fl2J5P/P7bNT54TnJ86eumll7QYThV+wk8+aVavXp0mTZpES5Ysoa5du9LEiROpqKhI+fkUfsbbT/5e8eyzz5LsIFbuYqYA4Ge8/Ux14rrrrnOG/uMnPnjIdNUL/Iy3nwcPHqSNGzc6Q6nysL78d/7GG2+ke+65xxladdasWUoVhZ/x9jOTfLoMQc2xwU/4yR6sXLmS5A03zlQoPISqvNnbecqY31O5wM94+1lQUEC//e1viUeAbNOmjfNb6fHHHy8etY6nO1O5wM94+5nJPZ36lOAn/MzkqC5l7GedOnUyhpM2DDV3hPLC8wjWq1cv405hFjZq1MipPrUz0H38vF27dmEevtS6mYvLqNQP4wNlIsA3CHDnFl9g6NevHw0ePNi56JCfn7inwWWvys8yNUjBh+FneNBN9JPnOXrssceUzkGfnBH4mUwj2HXT/OShVr761a/Sm2++Sc2aNaM33niDBg0aFCyUMtYGP8sIrAwfN8HP0aNHOzfV8Ku78IU7Xn74wx/SZZddRjxMoKoFfoZH3gQ/M7W+cePGdNNNN9Hy5cszvR1pGfwMD7cJfvLvdr44V7duXQ+Ibt26OReWP/zwQ7r99ts970W5AT/Do22Cn6mt54tjugxBzbHBz9QMBbdtip/8d5ynXVuzZg3xTWBNmjShH/zgB84NN/Pnzw8OSDlqgp/lgJbjLqb42adPH2cKFP7dXrlyZerfvz89//zzdMkll9DXv/71HFsbzsfgZzhcuVZT/EwloFOfEvxMzU5w26b6GRwB/zWV5GdaZzGf8KtUqeLcra7iLiGem5h7tnft2uVp+datW51tFTElB8J3+bVv3z65COs+CPCE3nxxge9Y4zvQ+UvyNddcU/xkcWrVqv1MjUe3bfgZbEZM9vOXv/wlXX755cRPH+mywM9gM2Gyny4JnsewR48eJKcWcIuUvcLPYNGb5iffjLZp0ybnQp1LYvfu3c4qzxErh19zi5W8ws9gsZvmZ7bW84g7qTfYZvtsmOXwM1i6pvnJv99r1KhBO3fu9IDg0UN4qVatmqc86g34GSxx0/xMbT3P+c5zFfM1Lx0W+BlsFkz0k0dc6ty5s9NRzDS+//3vO3PB8lzwPP9m6o04wRIruTb4WTKfsr5rop/cRv47L6eRcpo7ffp0ksOkE19vUjXXpssdfrokgnk11c/k1uvUpwQ/kzPjf90GP/1TCK6GkvxMPLL55fH4TqGvfOUr9Ne//jW4CMpQkxwvmwYOHOgMO5y8G99lx0MS89AXqhbu1OS7/twfvqrisOW4fCcI/1Djocl4+LJXX33V6TBmB7Itqv1Mjssd3sJ9TX5PxTr8DJa6yX6+8MILzhf3IUOGFENhP/79738Xb0e9Aj+DJW6yn6kk+McmdxirXOBnsPRN9POpp54iOWex59+wYcMcMHwnvZxvJlhIZagNfpYBVg4fNdHPbM3iKVr4CQ+VC/wMlr6JfvLF4u7duzs3myfT2L59u7PJ76la4Gew5E30M5UAX3Pg4dJ1WOBnsFkw1U/+LcQXwZOXu+++2xl1b+/evcnFka7Dz2Bxm+pnMoXPPvuMxo0bR5MnT1Z+HoWfyZnxv26an5x/XtxXl4AufUrw081IMK+2+JlMw3XXfU1+L+z10vxM6yzmgHr37k083MmZM2fCji9j/Y888ogzJwLPU8vLkSNHaOrUqfTkk086w15k3CmCQp5fcf/+/c7QyBEczspD8MUE9ut3v/sd8Rdf9ox/rBUWFubcXtV+uoG6T8K5Q6S75ape4ad/8jb4OWXKFOf/Lx7C6sUXXyS+85PnM7zhhhucYc78UypfDfCzfNyS9zLdT57r6JVXXqEdO3YUN+vdd9+lY8eOFd+pXPxGxCvw0z9w0/30TyC8GuCnf7am+3n27FlnHkP+3szrvPD58+OPP6bkG8P8kyp7DfCz7MxS9zDdT24Pf9fk0RhmzpxZ3LwFCxbQ1VdfTVdddVVxWdQr8NM/cRv8dCnwtRwegppHZdBhgZ/+s2CDn3fddZczJc/x48eLgaxevZo6dOhArVq1Ki6LegV++idug58uBR6ytFevXsS+jhgxwi1W9go//aM32c+Srsfr0KcEP+EnEyipv6gkh/3TK7mGUv2Uvclpy+bNm4U8aQj5JEXae1EVrFixQshOQTF+/HghnzQWEydOjOrQWY8j51oScljXrO/jjWgI6OCnfPJe3HPPPXwrk2jevLmQQwcJOd9yNACyHAV+ZgETcbFKP2XHsOMke5n6r2nTpkJeYI6YRuJw8DPBQuWaSj/lXfNCTjMh5F2B4lvf+paQc22Ke++9V8jOYpVInGPDT+UpcAJQ6WcmAmPHjnXOpSdOnMj0dmRl8DMy1CUeSKWf8gZeIUdgcHxs0KCBc/58+umnhbwJp8SYo3gTfkZBufRjqPTTjU7ezCBat24t5FNH4v777xe33XabkBdJ3LeVvMJPJdjTDqqDnxzUtGnTxKBBg9LiU1UAP1WR9x5XtZ/8t3z06NFCdg6LX//612LMmDFiwIABQt4Q5g004i34GTHwLIdT7Sdfm+e/6/LhAzFjxowsUUZfDD+jZ57piCr8zOV6vOo+JfiZyZboy3T1MxeHw6RVip9FeXzwTP3NPXv2dJ7ifeuttzK9HVnZJ598Qk2aNKH8/IwPQUcWBz8J1aJFC5Kd1jR8+PDIjosDZSagi5+Zo4u+FH5Gz7ykI8JPLx346eWhekuln/yVQ35hc75fNG7cWDUK5/jwU4s0FAeh0s/iIDRagZ8aJUOGotrPXbt2OU8WN2zYUAsw8FOLNBQHodpPDoTvkuf53/n3u+q5iuFnsRparOjg54YNGxwvdTiHwk8ttCwOQgc/T5486Zw/69Wr50zBVxycghX4qQB6CYdU5Sf/dl+2bJkzZzGfN/lJVB0W+KlDFhIxqPIzEUH2NRV9SvAzez5UvKOznyp45ODnnKydxfxIsryDnbizuE+fPiri1+qYQ4cOdYZ74zk/ed5cLGoJwE8vf/jp5aF6C356MwA/vTxUb8FPbwbgp5eH6i346c0A/PTyUL0FP70ZgJ9eHqq34Kc3A/DTy0P1Fvz0ZgB+enmo3oKf3gzATy8P1Vvw05sB+OnloXoLfnozAD+9PFRvwU9vBnLwM3tnMVfVt29fZy4sntelSpUq3tpjtLV8+XLq1q0bvfzyyySHxI5Ry/VuKvz8Ij/wU09P4Sf81NPML6KCn/ATfupMAH7qnB2cP+En/NSZAPzUOTs4f8JP+KkzAfipc3Zw/oSf8FNnAvBT5+zg/FkmP0vuLN6+fTvJeTNIzutCzz33nM55Dy22I0eOUMeOHeniiy8mOYezNsNuhNZggyqGn0TwU19h4Sf81NdOIvgJP+GnzgTgp87ZwfkTfsJPnQnAT52zg/Mn/ISfOhOAnzpnB+dP+Ak/dSYAP3XODs6fZfJzDvGcxSUts2bNEnJeAjF79uySPmble2fOnBH9+/cX9evXF7t377ayjaY3Cn7CT50dhp/wE37qSQB/3/XMS3JUOH/i/Jnsg27r8BN+6uZkcjzwE34m+6DbOvyEn7o5mRwP/ISfyT7otg4/4aduTibHAz/hZ7IPuq3Dz5z9LCq1s5iTO2rUKCHn6RVvv/22brkONZ6RI0cKOfy2kOObh3ocVO6PAPyEn/4MCndv+Ak/wzXMX+3wE376MyjcveEn/AzXMH+1w0/46c+gcPeGn/AzXMP81Q4/4ac/g8LdG37Cz3AN81c7/ISf/gwKd2/4CT/DNcxf7fAzJz9z6yzmJ3Buu+02UaNGDfHOO+/4y4wBe589e1Y8+OCDokKFCmLu3LkGRBzvEOFnvPOve+vhp+4Zind88DPe+de99fBT9wzFOz74Ge/86956+Kl7huIdH/yMd/51bz381D1D8Y4PfsY7/7q3Hn7qnqF4xwc/c8p/bp3FXNWpU6fEgAEDnCeM+dFtW5eTJ0+KgQMHikqVKomZM2fa2kzr2gU/rUupVQ2Cn1al07rGwE/rUmpVg+CnVem0rjHw07qUWtUg+GlVOq1rDPy0LqVWNQh+WpVO6xoDP61LqVUNgp9WpdO6xsDPUlOae2cxV8U98A888IAzh/EjjzzidCCXegiDPrBlyxZxxRVXiJo1a4pFixYZFDlCZQLwEx7oTAB+6pwdxAY/4YDOBOCnztlBbPATDuhMAH7qnB3EBj/hgM4E4KfO2UFs8BMO6EwAfuqcHcQGP0t0oGydxW5V06ZNE9WqVRNdunQRmzZtcouNfp0zZ46oXbu2aNu2rVi7dq3RbYl78PAz7gbo3X74qXd+4h4d/Iy7AXq3H37qnZ+4Rwc/426A3u2Hn3rnJ+7Rwc+4G6B3++Gn3vmJe3TwM+4G6N1++Kl3fuIeHfzMaED5Oou5qvXr14sOHTqIwsJC8dOf/lScOHEi4xF0L/z444/FDTfcIIhIDB8+XBw7dkz3kBFfDgTgZw6Q8BFlBOCnMvQ4cA4E4GcOkPARZQTgpzL0OHAOBOBnDpDwEWUE4Kcy9DhwDgTgZw6Q8BFlBOCnMvQ4cA4E4GcOkPARZQTgpzL0OHAOBOBnGqTydxZzVTzO9/jx40X16tVF8+bNxYwZM8Tp06fTjqJjwd69ewUPpc2d3W3atBGLFy/WMUzE5IMA/PQBD7uGTgB+ho4YB/BBAH76gIddQycAP0NHjAP4IAA/fcDDrqETgJ+hI8YBfBCAnz7gYdfQCcDP0BHjAD4IwE8f8LBr6ATgZ+iIcQAfBOCnB56/zmK3qu3bt4s77rhDFBQUiJYtW4rp06eL48ePu29r9bpt2zbx8MMPOx3c9erVExMmTLBu7mWtgGsQDPzUIAkIISsB+JkVDd7QgAD81CAJCCErAfiZFQ3e0IAA/NQgCQghKwH4mRUN3tCAAPzUIAkIISsB+JkVDd7QgAD81CAJCCErAfiZFQ3e0IAA/HSSEExnsZvPjRs3iiFDhohKlSqJOnXqiFGjRmkx/y/fITB//nxx4403igoVKoj69euLZ555Rhw9etQNHa8xIAA/Y5Bkg5sIPw1OXgxCh58xSLLBTYSfBicvBqHDzxgk2eAmwk+DkxeD0OFnDJJscBPhp8HJi0Ho8DMGSTa4ifDT4OTFIPSY+xlsZ7Hry+7du8XYsWNF06ZNnbmAW7VqJcaMGSNWrVoV2TDVhw4dEvPmzXOeeOaO67y8PNG7d29RVFQkTp486YaK1xgSgJ8xTLpBTYafBiUrhqHCzxgm3aAmw0+DkhXDUOFnDJNuUJPhp0HJimGo8DOGSTeoyfDToGTFMFT4GcOkG9Rk+GlQsmIYakz9LMrjXFNIy9mzZ2nZsmX02muv0dy5c2nr1q1Uq1Yt6tatm/OvY8eO1K5dO2rQoIGvCOQ8ySR7/Wnt2rW0fPlyWrp0Ka1evZr4+F27dqX+/fs7/y666CJfx8HOdhGAn3bl07bWwE/bMmpXe+CnXfm0rTXw07aM2tUe+GlXPm1rDfy0LaN2tQd+2pVP21oDP23LqF3tgZ925dO21sBP2zJqV3ti5uecUDuLU9Xgztz33nuPlixZQu+//z7t3LnT+Yh88pdatGhBcnhoaty4Mcm5hJ1O5cqVK1PVqlWJXw8fPkzcKcyv8qlhkuOI0549e0jOQUwbNmwgOdQ0yTmTqXXr1tSjRw/q3r278+/8889PDQPbIJCRAPzMiAWFmhCAn5okAmFkJAA/M2JBoSYE4KcmiUAYGQnAz4xYUKgJAfipSSIQRkYC8DMjFhRqQgB+apIIhJGRAPzMiAWFmhCAn5okAmFkJGC5n9F2FqcS3rdvH3300Ue0bt062rRpE8nHu+m///2v0wnMHcJyuGiS8wo7HcHVq1enihUrUo0aNahmzZrUsGFDp3O5UaNGJIe5prZt21KbNm2cjuXU42AbBMpDAH6Whxr2iYoA/IyKNI5THgLwszzUsE9UBOBnVKRxnPIQgJ/loYZ9oiIAP6MijeOUhwD8LA817BMVAfgZFWkcpzwE4Gd5qGGfqAjAz6hI4zjlIWCZn2o7i8uTAOwDAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiDgm8CcfN9VoAIQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAHjCKCz2LiUIWAQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEQ8E8AncX+GaIGEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEDCOQIGMeI5xUSNgEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABEAABPwSW/x9ZSpLWk8rmogAAAABJRU5ErkJggg==",
      "text/plain": [
       "<PIL.PngImagePlugin.PngImageFile image mode=RGBA size=1931x59>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.transpiler import CouplingMap\n",
    "from rustworkx.visualization import graphviz_draw\n",
    "\n",
    "# L = length of our 1D spin chain\n",
    "L = 12\n",
    "\n",
    "# Generate the edge-list for this spin-chain\n",
    "elist = [(i-1,i) for i in range(1, L)]\n",
    "# Generate an edge-coloring so we can make hw-efficient circuits\n",
    "even_edges  = elist[::2] \n",
    "odd_edges = elist[1::2]\n",
    "\n",
    "elist = even_edges + odd_edges\n",
    "\n",
    "# Instantiate a CouplingMap object\n",
    "coupling_map = CouplingMap(elist)\n",
    "graphviz_draw(coupling_map.graph, method=\"circo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate the hamiltonian we will evolve by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIIIIIIIIIII', 'IIIIIIIIIIXX', 'IIIIIIIIIIYY', 'IIIIIIIIIIZZ', 'IIIIIIIIXXII', 'IIIIIIIIYYII', 'IIIIIIIIZZII', 'IIIIIIXXIIII', 'IIIIIIYYIIII', 'IIIIIIZZIIII', 'IIIIXXIIIIII', 'IIIIYYIIIIII', 'IIIIZZIIIIII', 'IIXXIIIIIIII', 'IIYYIIIIIIII', 'IIZZIIIIIIII', 'XXIIIIIIIIII', 'YYIIIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIIIYYI', 'IIIIIIIIIZZI', 'IIIIIIIXXIII', 'IIIIIIIYYIII', 'IIIIIIIZZIII', 'IIIIIXXIIIII', 'IIIIIYYIIIII', 'IIIIIZZIIIII', 'IIIXXIIIIIII', 'IIIYYIIIIIII', 'IIIZZIIIIIII', 'IXXIIIIIIIII', 'IYYIIIIIIIII', 'IZZIIIIIIIII'],\n",
      "              coeffs=[1.        +0.j, 0.52440675+0.j, 0.52440675+0.j, 1.0488135 +0.j,\n",
      " 0.60759468+0.j, 0.60759468+0.j, 1.21518937+0.j, 0.55138169+0.j,\n",
      " 0.55138169+0.j, 1.10276338+0.j, 0.52244159+0.j, 0.52244159+0.j,\n",
      " 1.04488318+0.j, 0.4618274 +0.j, 0.4618274 +0.j, 0.9236548 +0.j,\n",
      " 0.57294706+0.j, 0.57294706+0.j, 1.14589411+0.j, 0.46879361+0.j,\n",
      " 0.46879361+0.j, 0.93758721+0.j, 0.6958865 +0.j, 0.6958865 +0.j,\n",
      " 1.391773  +0.j, 0.73183138+0.j, 0.73183138+0.j, 1.46366276+0.j,\n",
      " 0.44172076+0.j, 0.44172076+0.j, 0.88344152+0.j, 0.64586252+0.j,\n",
      " 0.64586252+0.j, 1.29172504+0.j])\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from qiskit.quantum_info import SparsePauliOp, Pauli\n",
    "\n",
    "# Genereate random coefficients for our XXZ Hamiltonian\n",
    "np.random.seed(0)\n",
    "Js = np.random.rand(L-1) + 0.5*np.ones(L-1)\n",
    "Js_fixed = np.zeros(L-1)\n",
    "\n",
    "hamiltonian = SparsePauliOp(Pauli('I'*L))\n",
    "for i, edge in enumerate(elist):\n",
    "    Js_fixed[edge[0]] = Js[i]\n",
    "    hamiltonian += SparsePauliOp.from_sparse_list(\n",
    "    [(\"XX\", (edge), Js[i] /2),\n",
    "     (\"YY\", (edge), Js[i] /2),\n",
    "     (\"ZZ\", (edge), Js[i] )], \n",
    "    num_qubits=L)\n",
    "\n",
    "print(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the observable we want to estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIIIIZZIIIII'],\n",
      "              coeffs=[1.+0.j])\n"
     ]
    }
   ],
   "source": [
    "observable = SparsePauliOp.from_sparse_list(\n",
    "    [(\"ZZ\", (L//2 - 1, L//2), 1.0)], num_qubits=L\n",
    ")\n",
    "print(observable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1c: Construct a Trotter circuit\n",
    "We can use the `generate_time_evolution_circuit` function which is built on top of qiskit's `PauliEvolutionGate` to construct our trotter circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwcAAAC+CAYAAABgdXNvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYE0lEQVR4nO2deXiU1fm/72Qy2VcCJAFC2JewkwAiVrGKbFLsQl1qFbAUW5cQrNX+2m9V6gZIECvijtWKuFStJaKIigrKmrCEJWGH7DshezLz/v6YJBBJMu87JMyZ4bmvKxeTmfN5cs7cmZCTec85HpqmaQiCIAiCIAiCcNnj6ewOCIIgCIIgCIKgBjI5EARBEARBEAQBkMmBIAiCIAiCIAgNyORAEARBEARBEARAJgeCIAiCIAiCIDQgkwNBEARBEARBEACZHAiCIAiCIAiC0IBMDgRBEARBEARBAGRyIAiCIAiCIAhCAx02OfDw8ODjjz/W3X727NncdNNNF/U1T5w4gYeHB7t373a4xqOPPsrIkSMvqh+CIAiCIAiC4IoYnhzk5uaSkJBAv3798PX1JSIiggkTJrBq1SoqKys7oo/tysSJE/Hw8MDDwwNfX19iY2N54YUXnN0tw2zatKlpHOd/5ObmtpnTNI1nnnmGAQMG4OPjQ/fu3XniiSeatXn77bcZMWIE/v7+REVFMXfuXIqKipoe//DDD4mPjyc0NJSAgABGjhzJW2+91azGhx9+yA033EB4eHirE7bc3Fx++9vfEhkZSUBAAKNHj+Y///lPsza9evW6YIxPP/100+Pp6elce+21RERE4OvrS58+ffjb3/5GXV1dszrvv/8+gwYNwtfXl2HDhvHpp582e7yl59LDw4OlS5c2tXniiSe48sor8ff3JzQ09ILxvPHGG63Wyc/Pb2q3cuVKBg8ejJ+fHwMHDuTNN9+0W8fX17dZm/Lycu6991569OiBn58fsbGxvPjii4afX0EQBEEQhPPxMtL42LFjTJgwgdDQUJ588kmGDRuGj48P+/bt4+WXX6Z79+787Gc/66i+thvz5s1j0aJFVFZW8uabb3LPPfcQFhbGrbfe6uyuGSY9PZ3g4OCmz7t27dpm+4SEBDZs2MAzzzzDsGHDKC4upri4uOnxLVu2cMcdd7B8+XJmzJhBVlYWd999N/PmzePDDz8EoFOnTvz1r39l0KBBeHt7s27dOubMmUPXrl2ZPHkyABUVFVx11VX8+te/Zt68eS325Y477qC0tJRPPvmEzp07s2bNGn7961+zc+dORo0a1dRu0aJFzWoEBQU13Tabzdxxxx2MHj2a0NBQ9uzZw7x587BarTz55JMAfP/999x666089dRT3HjjjaxZs4abbrqJlJQUhg4dCkBOTk6zvq1fv5677rqLX/7yl0331dbWMmvWLMaPH89rr712wXhuvvlmpkyZ0uy+2bNnU11d3eRl1apV/OUvf+GVV15hzJgxbN++nXnz5hEWFsaMGTOacsHBwaSnpzd97uHh0azuwoUL+eqrr/j3v/9Nr1692LBhA3/84x/p1q1b02tQ7/MrCIIgCILQhGaAyZMnaz169NDKy8tbfNxqtTbdBrSPPvqo6fO9e/dq1157rebr66t16tRJmzdvnnb27Nmmx++8805t5syZ2qOPPqp17txZCwoK0ubPn6/V1NQ0tVm/fr02YcIELSQkROvUqZM2ffp07ciRI02PHz9+XAO01NTUVsdwzTXXaAkJCc3u69+/v3bLLbdomqZpjzzyiDZixAjtzTff1GJiYrTg4GDt5ptv1srKynT3o6amRrvnnnu0yMhIzcfHR+vZs6f25JNPNj1eUlKi3XXXXU3jvPbaa7Xdu3e32ueW+PrrrzVAKykp0Z05cOCA5uXlpR06dKjVNkuXLtX69OnT7L7nnntO6969e5u1R40apf3tb3+74P62nAQEBGhvvvlms/s6deqkvfLKK02fx8TEaMuXL2/za/+YxMRE7aqrrmr6/Ne//rU2ffr0Zm3GjRunzZ8/v9UaM2fO1H7605+2+Njq1au1kJAQu/3Iz8/XzGZzszGOHz9e+9Of/tSs3cKFC7UJEyYYqj9kyBBt0aJFze4bPXq09te//rXpcz3PryAIgiAIwvnovqyoqKiIDRs2cM899xAQENBimx//dbORiooKJk+eTFhYGDt27OD9999n48aN3Hvvvc3affnllxw8eJBNmzbxzjvv8OGHH/LYY481q7Nw4UJ27tzJl19+iaenJz//+c+xWq16h9Eifn5+1NbWNn1+9OhRPv74Y9atW8e6dev45ptvml3KYq8fzz33HJ988gnvvfce6enpvP322/Tq1aspP2vWLPLz81m/fj27du1i9OjRXHfddU1/wW9cO7Fp0ya7fR85ciRRUVFMmjSJLVu2tNn2f//7H3369GHdunX07t2bXr168bvf/a7ZOwfjx4/n9OnTfPrpp2iaRl5eHh988AHTpk1rsaamaXz55Zekp6dz9dVX2+3v+Vx55ZW8++67FBcXY7VaWbt2LdXV1UycOLFZu6effprw8HBGjRrF0qVLqa+vb7XmkSNH+Oyzz7jmmmua7vvhhx+4/vrrm7WbPHkyP/zwQ4s18vLySE5O5q677jI0nh/z5ptv4u/vz69+9aum+2pqai64RMjPz4/t27c3uxSqvLycmJgYoqOjmTlzJvv372+WufLKK/nkk0/IyspC0zS+/vprMjIyuOGGG5q10fP8CoIgCIIgNKF3FrF161YN0D788MNm94eHh2sBAQFaQECA9uc//7npfs575+Dll1/WwsLCmr3jkJycrHl6emq5ubmaptneOejUqZNWUVHR1GbVqlVaYGCgZrFYWuxTQUGBBmj79u3TNM34Owf19fXaW2+9pQHa888/r2ma7Z0Df3//Zu8UPPjgg9q4ceNarfnjftx3333aT3/602bvpDTy3XffacHBwVp1dXWz+/v27au99NJLmqZpWmZmpjZw4EBt27ZtrX7NQ4cOaS+++KK2c+dObcuWLdqcOXM0Ly8vbdeuXa1m5s+fr/n4+Gjjxo3Tvv32W+3rr7/WRo4cqV177bXN2r333ntaYGCg5uXlpQHajBkztNra2mZtSktLtYCAAM3Ly0vz8fHRXnvttRa/ZltOSkpKtBtuuEEDNC8vLy04OFj7/PPPm7VZtmyZ9vXXX2t79uzRVq1apYWGhmqJiYkX1Bo/frzm4+OjAdrvf//7Zt8zZrNZW7NmTbP2K1eu1Lp27dpinxcvXqyFhYVpVVVVLT6u952DwYMHa3/4wx+a3feXv/xFi4yM1Hbu3KlZrVZtx44dWkREhAZo2dnZmqZp2vfff6/961//0lJTU7VNmzZpN954oxYcHKydPn26qU51dbV2xx13ND133t7e2r/+9a9mX0vP8ysIgiAIgnA+Fz05OHr0qHb48GFt7NixzS7XOX9ykJiYqE2cOLFZrrS0VAO0b775RtM02+Tgx7+k7t69WwO0EydOaJqmaRkZGdott9yi9e7dWwsKCtICAgI0QEtOTtY0Tf/kwGw2awEBAZrZbNZ8fX21xMTEpl8mH3nkES02NrZZJikpSevdu3fT5/b6sWvXLq1Tp05a//79tfvuu6/ZL2TPP/+85unp2TShavzw9PRsNrlyhKuvvlq7/fbbW3183rx5GqClp6c33bdr1y4NaLrUaP/+/VpUVJS2ZMkSbc+ePdpnn32mDRs2TJs7d26zWhaLRTt8+LCWmpqqPfPMM1pISIj29ddfX/A123Jy7733amPHjtU2btyo7d69W3v00Ue1kJAQbe/eva2O4bXXXtO8vLwumFydOnVK279/v7ZmzRqte/fu2uLFi5seMzo5GDhwoHbvvfe22gc9k4Pvv/9eA7SdO3c2u7+ysrJpImcymbRu3bppf/7znzWgaaL8Y2pra7W+ffs2u2xr6dKl2oABA7RPPvlE27Nnj/bPf/5TCwwM1L744oumNo48v4IgCIIgXN7oXpDcr18/PDw8mi2SBOjTpw9guzSio5kxYwYxMTG88sordOvWDavVytChQ5tdEqSH3/zmN/z1r3/Fz8+PqKgoPD2bX11lNpubfe7h4dHs0iV7/Rg9ejTHjx9n/fr1bNy4kV//+tdcf/31fPDBB5SXlxMVFdXiJUMt7YBjhLFjx7J58+ZWH4+KisLLy4sBAwY03Td48GAATp06xcCBA3nqqaeYMGECDz74IADDhw8nICCAn/zkJzz++ONERUUB4OnpSb9+/QDbpU0HDx7kqaee0n3JytGjR3n++edJS0tjyJAhAIwYMYLvvvuOlStXXrDzTiPjxo2jvr6eEydOMHDgwKb7o6OjAYiNjcVisfD73/+eBx54AJPJRGRkJHl5ec3q5OXlERkZeUH97777jvT0dN59911d42iNV199lZEjRxIXF9fsfj8/P15//XVeeukl8vLyiIqK4uWXXyYoKIguXbq0WMtsNjNq1CiOHDkCQFVVFf/v//0/PvroI6ZPnw7YPO3evZtnnnmG66+/3uHnVxAEQRCEyxvdaw7Cw8OZNGkSzz//PBUVFYa+yODBg9mzZ0+z3JYtW/D09Gz2C96ePXuoqqpq+nzr1q0EBgYSHR1NUVER6enp/O1vf+O6665j8ODBlJSUGOpHIyEhIfTr14/u3btfMDGwh95+BAcHc/PNN/PKK6/w7rvv8p///Ifi4mJGjx5Nbm4uXl5e9OvXr9lH586dHRpPI7t372765b0lJkyYQH19PUePHm26LyMjA4CYmBgAKisrL3hOTCYTYFtf0BpWq5WamhrdfW3c9ralr9XWGpLdu3fj6enZ5q5MVquVurq6pjrjx4/nyy+/bNbmiy++YPz48RdkX3vtNeLi4hgxYoTusfyY8vJy3nvvvTbXLJjNZnr06IHJZGLt2rXceOONrX4vWiwW9u3b1+S2rq6Ourq6Np87R59fQRAEQRAubwxtZfrCCy8wYcIE4uPjefTRRxk+fDienp7s2LGDQ4cOXfBX0kZ+85vf8Mgjj3DnnXfy6KOPUlBQwH333cdvf/tbIiIimtrV1tZy11138be//Y0TJ07wyCOPcO+99+Lp6UlYWBjh4eG8/PLLREVFcerUKR5++OGLG70D6OlHUlISUVFRjBo1Ck9PT95//30iIyMJDQ3l+uuvZ/z48dx0000sWbKEAQMGkJ2dTXJyMj//+c+Jj48nKyuL6667jjfffJOxY8e22I9nn32W3r17M2TIEKqrq3n11Vf56quv2LBhQ1Ob559/no8++qjpF+Prr7+e0aNHM3fuXJ599lmsViv33HMPkyZNano3YcaMGcybN49Vq1YxefJkcnJyWLBgAWPHjqVbt24APPXUU8THx9O3b19qamr49NNPeeutt1i1alXT1y4uLubUqVNkZ2cDNL3jFBkZSWRkJIMGDaJfv37Mnz+fZ555hvDwcD7++GO++OIL1q1bB9gWEm/bto1rr72WoKAgfvjhBxITE7n99tsJCwsDbGcymM3mpm11d+7cyV/+8hduvvnmpneAEhISuOaaa1i2bBnTp09n7dq17Ny5k5dffrnZc1pWVsb777/PsmXLWnzOT5061TQui8XSdHZDv379CAwMbGr37rvvUl9fz+23335BjYyMDLZv3864ceMoKSkhKSmJtLQ0/vWvfzW1WbRoEVdccQX9+vWjtLSUpUuXcvLkSX73u98BtonnNddcw4MPPoifnx8xMTF88803vPnmmyQlJQHoen4FQRAEQRAuwOh1SNnZ2dq9996r9e7dWzObzVpgYKA2duxYbenSpc0WE+PgVqZ///vftfDwcC0wMFCbN29es2vLv/jiC23w4MGaj4+PNnz4cG3Tpk3Nvo6jW5meT+NWpuezfPlyLSYmRnc/Xn75ZW3kyJFaQECAFhwcrF133XVaSkpKU76srEy77777tG7dumlms1mLjo7WfvOb32inTp1qNo6WruFvZPHixVrfvn2bns+JEydqX3311QVjOb/fmqZpWVlZ2i9+8QstMDBQi4iI0GbPnq0VFRU1a/Pcc89psbGxmp+fnxYVFaX95je/0TIzM5se/+tf/6r169dP8/X11cLCwrTx48dra9eubVZj9erVGnDBxyOPPNLUJiMjQ/vFL36hde3aVfP399eGDx/ebOvNXbt2aePGjdNCQkI0X19fbfDgwdqTTz7Z7Hti7dq12ujRo7XAwEAtICBAi42N1Z588skLFhO/99572oABAzRvb29tyJAhTetDzuell17S/Pz8tNLS0haf8zvvvLPFMf3Y0/jx47XbbrutxRoHDhzQRo4cqfn5+WnBwcHazJkzL9hadsGCBVrPnj01b29vLSIiQps2bVqz7x9N07ScnBxt9uzZWrdu3TRfX19t4MCB2rJly5otgrf3/AqCIAiCIPwYD01r41oRQRAEQRAEQRAuG4xdcC8IgiAIgiAIgtsikwNBEARBEARBEACZHAiCIAiCIAiC0ICh3YqciaW2jpIDJ3W3D4u1bc1pNGPyNttvKDRD3KiLuFEXcaMu4kZdxI26iBv3wWUmByUHTrJuqv6tS29c/zSA4Uznkf0M9+1yR9yoi7hRF3GjLuJGXcSNuogb90EuKxIEQRAEQRAEAZDJgSAIgiAIgiAIDbT75GDjxo0kJSUxatQoPvzww/YuT89p4xg0e0qHZ/SwLR/mbYZFqWCV0yKUcrOjAH6/GR4TN4BabnYV2tw8miJuQC03KYXw+y3wSApYxI1SblIb3Pxd3ABqudldLG7ORyU3exvc/F8K1Fvbvbzb0u6Tg/T0dO6//37mzZuHr69ve5fH0+TJmWPZeIcGdmhGD+8ch9Ri+CwLjp1t19IuiUpu1hyDlAY3R8ratbRLopKbt4+ec5MhbtRycwxSimxu0kvbtbRLoqSbTDhQ0q6lXRKl3Bw552afuFHKzb+PnnOzt7hdS7s17b4guby8nMWLFzN06FCqqqpabLNixQpOntS/Oh3Ar6iG/oBmtRLStxs53+1rs/3y5csBDGeqwn109ykzeiLePa7Gw2rln4uWYLLW6s66E0q66TER7+ir8dCsrHx8CSaLuFHHzTV497gGsPLC40vxstTozroTSrrp/hO8o68FNF58com4UchNVver8I7+KaDxylNL8bJU6866E6q7ef3ppXjVixtl3HS7Eu+e1wOweslS/l3f8u+llwMxMTEkJCToauuhaVq7vwmWnJzM8OHDiY6ObreahbuPKLkKvrwOXjwEfxpmKOZWiBt1ETfqIm7URdyoi7hRF3HjPnTIVqbTp0/viLIXMGjuVExmL0rST2Gts+BpNmGts5C7Jc1upuxELtUFpfiEBVFTcpaClMMO9SHQDJ4ejo7AfRE36iJu1EXcqIu4URdxoy7ixjVx7d2KNA2f8CB8QgMpO57T9K+ejMnbi4LUI3j5+1CQeuTS9PdyQtyoi7hRF3GjLuJGXcSNuogbl8SlJwc1xWc59sG3BHTvQkjfblQXlRHSt5uuTFXhGdA0KvNKoP2vrLrsETfqIm7URdyoi7hRF3GjLuLGNXGZE5Jb4vh/twBQmpFpONNI/vZD7donwYa4URdxoy7iRl3EjbqIG3URN66JS79zIAiCIAiCIAhC++G0dw5q6yzsSbe/6eyIgZ3wNpsuQY8untWrV1NZWcnEiROpqqoiKyuL2NhYtm7dyqxZs3jjjTcYNWoU48aNw1JbR8kB/du5BvfrTtmRLN3tw2JjMHmbHRmGWyJu1EXcqIu4URdxoy4quQHxcz7iRh9OmxzsSS9m7G2f2G23fc3PGDO0ywX3R08eg2/nEApTDxMWG4Nvp2CK9h3DAw9yf9jfYi1HMkaYM2dO0+2cnBzi4+MB6N+/PwB33303ubm5AJQcOGlo+66rnr2HzQtW6m7vyHZf7YW4aRtx0xxxY0PctI24aY64sSFu7OMsP+LGPs587bSFy645OP35jqbbemd2jmQcJSoqqsX7IyMjO/TrqoC4URdxoy7iRl3EjbqIG3URN66LrDkQBEEQBEEQBAHooMnBxo0bSUpKIjY2lrS01g+6+DG3TO3Dryb1onOYb6ttuo4bTO+ZEwz1x5GMHtZnws1fw6YcqLc6VqPntHEMmj2lwzP2qLNCwlbbeD7Tv6lAM1Ry89l5bupc3E29FRY0uFnvBm42NLj5OgdqLY7VUMlN4jbbeJJPO1ZDJTdfZMEtbuRmYYObdaccq6GSmy+zz7mpcQc3221uPnHwj8MqufnKjdxYNPjTdvj11/CxG7jZlA23boKvcqDaDdw82PC6+ahj31TpmMlBeno6999/P/fff7+h3Le7chk/IoLCkupW2+RvO2i4P45k9PB5Jhw9C7lVcKLcsRqeJk/OHMvGOzSwQzP2OHYWthbYxvO5sfU0TSjlJss2lpwq29gcQRU3x8/CtkY3Dk4OVHLzWda5142ruzlZfu51s8ENXjefZcKRhtfN4TLHaqji5nQFbG90k+1YDZXcrD/PTfoZx2qo4iazErbn29x84Q5uss65OVjqWA1V3GRV2P6/OXYWNrqBm0+zbD/Lcqtgf4ljNVRxk1MJ2wttrxtH3eilQ9YclJeXs3jxYuLj48nIyGDo0KHNHl+xYgXb9uYBfZvdP/nK7qzffJroyABO51YAsHz5ciJDqvErqqE/EDqgB6EDo+32Yfny5QCGM1XhPrrGCJAfGY9f9E8BD55/7Bk8Nf3T0sbxaFYrIX27kfPdvjbbr3nnHXo23NaTMToWq4eJoGF3URHQjbxN61j47k7dWRXd5HWNwy/mOsCDVYtc303g0LmUB3Yn/9t1LHzP1d2Mxi/mesCDFx9fhqe1XndWbTfJLHxvh/1QA6q7eeUJV3fjid/QuVgCe1Dw7acsfG+77qyKbvK7jsQv5gbAk9eeXIbJWqc7q6Ib/yFzsQT1oHDzZyx8f6vurJJuuozEr9ckwMTqp9RyA8bGo3l44jdkDnVB0RRu/pyF7/+gKweKuuk8DP/eU9Aw8ebTy3jbxd34DJlNTVBPCrdsYOH73+vKNRITE0NCQoKuth6a1jHHziUnJzN8+HCio1uWvSOtwNBuRYW7jxhaNX7j+qcBDGcu1apxo+O5VLtHJKXBwqH2252PuBE34kbciBtxYxRx0/FuwLHxiBv3cmOUDtutaPr06R1V2m0ZNHcqJrMXJemnsNZZ8DSbsNZZyN3S+rqNxkzZiVyqC0rxCQuipsTBazWEVhE36iJu1EXcqIu4UZf2dFOQcvgS9tz9uVzcyG5FKqFp+IQH4RMaSNnxnKZ/9WRM3l4UpB7By9+HgtQjl6a/lxPiRl3EjbqIG3URN+oibtTlMnHjtHMORgzsxPY1P9PVDmynyDW+BaWHsNgYAIcyzqKm+Cy5W9LoMSmekL5lVBeVEdK3G5XZRXYzPuHBoGlU5pVAx1wp1iripu2MuOlYxI39jLMQN/YzzkLc2M84C3FjP+MsXNWNUZw2OfA2m1o8+bg1TN5mh64zU/HkudY4/t8tAJRm6N+SpjHTSP72Q+3aJz2Im7YzjYibjkHcqIu4URdxoy7iRl1c1Y1R5LIiQRAEQRAEQRAAJ75zYBRLbZ2ho7Qb33rSmwmLjcHkbXaob+7I6tWrqaysZOLEiVRVVZGVlUVsbCxbt25l1qxZrF27FpPJxKxZs9i7ew99fMN11xY3F4e4UReV3DRmxI8NcaMu4kZdxI26GHGzf/9+4uLidNd2mclByYGTHbpF1qXcHqsloiePwbdzCIWphwmLjcHL15uClMP4R4SR9fVuw5mLZc6cOU23c3JyiI+PB6B///4AzJ49u+nxgHIr634ubvRmLhZxcw5xYz/jLD/ixn5G3NgQN+dobzetZfQibs7hym66d+9uqLbLTA7cndOfnzs86fxZdFsz6rYyo/58S7v1LSoqqs3HIzrrXzviiogbdRE36iJu1EXcqEt7u2lPxI3ruomMjDRUT9YcXCTV9WBVe9G5bqwa1Fud3Yv2Q9yoi7hRl+p6sIgbJRE36iJu1EXcGKdDJgcbN24kKSmJa6+9lv/+97/tWrvntHEMmj3lkuXa4u0j8Iuv4K7Ntm8+R+g6bjC9Z07o8Iw9quph7mZYdxreOeZYDZXcrDkKv/zKNqYqF3dTXW/7Hlt32vY95wgquVl7zOZmzndQ6epuLPC7Bjf/dgM37x1veN18BxVu4GbeFpubNx08a0glN+8fP/e6Ka9zrIYqbmrOc/OGG7j54Dj86muY/a3ru6m1wO+3wP9Ow+tu4Oajk+7jps4K8xvdZLRr6Qvw0LT232x15cqVzJ8/n5deeokePXowc+bMZo+vWLGCkyeNvaXiV1RD/w159JoxnpozFRTtPUZtaXmr7Q/fEAFA/w15AHZzh2+IoCrcx1Cf0obeRXHnIWC1ELdrGQEVuYbHA9B75oQLtrr6MafGdaLntuKmz+1ljI6nPKAbqaMXoJm86FSYxtC013Vn1XQzl+LOQ8FqYfSu5QRWZBseD6jhpiIgipTRCWgmbzoV7mdo2mu6s2q6mUNx52E2NynPElieZXg8oIgb/whS4xKxKujGkfHsHzKboi7DQbMyctezBJfr365PNTeV/l1JiUvEavIhrOgAw/a9qjurvJuU5wg+e8rweEARN35dbG68fAkrOsiwfa/ozqrp5k6KuowATWNE6nOElOn//aaj3YCx8VT5dWZXXCJWLz9Ciw4y3MXdHIi9g8KuI21udj9PyJnjhscDirjx7URK3EIsZn9Ciw4xfN/LunKNxMTEkJCQoKtth6w5KC8vZ/HixdTX1zN27NgLHtfbufMp3H2EdRseRrNaCenbjZzv9rXZPjExEYB1G2wLXezlEhMTDS9y+T4PXs2AsjoTqxb9GU8P/dnG8YQO6EHowGi77W+79VY2b1sJoCtjdDwWDR5NhawKmHfFUMbPTdKdVdHND/nwSjqcqTOxatGfMLmwG2uDm8wK+N24IVzp4m625sPLjW4ee8Dl3TyWCqcr4K6xQ5igkJvGjJHxbC+Al9KhtNaTlx5b6PJuFqXCqQqYMzaWn8xxbTc7C+GFQ1Ba48mLjyzAy8D7/qq50TRYtBtOlsPssYO52sXd7GpwU1LjwUt/T1DKDRgbj6bBP3bD8XK4M34wE13cTUqDm+IaD176v/tc3s3ju+FYOfw2fhA/NeDGKB0yOXjooYdITk5m+PDhREfbf3KMcDJ52yXNtcWVEbaPpDQMTQzOpzQjk9Qlazs8Yw+TB/xj9MXVUMnN+K62j6Q0DP2Ccz6quPH0gEVu5OaKrrYPd3HzmBu5GdvF9uEubh51IzfxneH1q2xujPyCcz6quPHwgEdGXVwNldzEdYbX3MjN393IzejO8Kobufm/i3Sjlw7brWj69OkdVbqJQXOnYjJ7UZJ+CmudBU+zCWudhdwtaXYzZSdyqS4oxScsiJqSsx3e18sNcaMu4kZd2tNNQYqDFwwLLSJu1EXcqIu4cU1ce7ciTcMnPAif0EDKjuc0/asnY/L2oiD1CF7+PhSkOriSUGgdcaMu4kZdxI26iBt1ETfqIm5cEpc+56Cm+Cy5W9LoMSmekL5lVBeVEdK3G5XZRXYzPuHBoGlU5pXYLuS6xITFxjQdAKKH4H7duVHH9Wvn13cm4qbt+s5E3LRd35mIm7brOxNx03Z9ZyJu7H8NZyFu7H8NFXHpyUHjKvDSDP07avx45Xj+9kPt2ie9mLzNhhdyOvMkWqOIG3URN+oibtRF3KiLuFEXceOauPZlRYIgCIIgCIIgtBtOe+egts7CnvRiu+1GDOyEt9l0CXp08axevZrKykomTpxIVVUVWVlZxMbGsnXrVqZMmcL333/PqVOnGDt2LPGjRus+PjssNgaTtxlLbZ3hjGBD3KiLuFEXcaMu4kZdxI26iBt9OG1ysCe9mLG3fWK33fY1P2PM0C4tPhY9eQzBfaLI/mYP/pGdCOzRhcPvfIW1rvWjPVvLtAdz5sxpup2Tk0N8fDwA/fv3B2g6DC43N5eSAydZN/VhXXVvXP80nUf2cyjjLMSN/YyzEDf2M86iPd20ldGLuDmHuLGfcRbixn7GWYgb+xkVcek1B6c/39F0W+9MzZGMI0RFRbX6WGRkJIW57r3yXtyoi7hRF3GjLuJGXcSNuogb10TWHFwkJ8qh1ursXrQP5XW2EyvdhZPiRlnczc0JN3JzqhxqLc7uRfvQ6MYJG510CG7n5qy4UZGKevdyc1rcGKZDJgcbN24kKSmJG264gQ8++EB37papffjVpF50DvNttU3XcYPpPXOC4T45mmuL5w/AHd/AF1m2H3SO0HPaOAbNntLhGXuU1cLczfDbb+CFg47VUMnNCwdtY/kiC866uJvyunNuVrqBmxcPwW+/hQ1ZcKbWsRoqubmrwc0/DzhWQyU3Lx2C27+FDdlQ6uJuKurhdxfpBhx7njvCzSvpttfN59lQWuNYDVXcVDa6+RZWuIGb185zU+LibqrOc/PsfsfrqOLm9Qzbz7TPs6DYxd1U18O8zXD7N7D8ItzowUPT2n/+sXLlSubPn8+SJUuIjIxk7ty5zR5fsWIF2/bm8c6Ovs3u79bVnwfuGMYDz5w7QvvWMUeJDKnGr6iG/hvyAOg9c8IFW139mMM3RAA0ZezlDt8QQVW4j/5BAmlD51LceSgeljridi3DvzJfd7ZxPL1mjKfmTAVFe49RW9rynx8b++ZIRi8VAZGkjl6A1eRNp8I0hqa9bngsoJKbORR3HgaWeuJ2LSOgMs9+qAHl3PhHkBq3AKvJxy3c7B8yh6Iuw8BqYfSuJAIr7ByIcx6quan060JKXCJWL186Fe5naNprhscCHeOmMWNkPPuHzKaoy3Cbm5TlBJZn686q5qbKrzMpoxOxmP0IK9rPsH2OuQF9zzN0sJvYOynqOgKsFkalrCCoXP/WkMq58Q0nJW5hg5sDDNv3quGxNKKGmzso6joSrBZGpj5H8NnTurPquelEatxC6s3+hBUdZNi+VwyPpREV3ByIvYPCBjcjdq8kpOyE7qxqbqp9QkmJf4B6cwBhRYcYtu9l3VmAmJgYEhISdLXtkDUH5eXlLF68mPr6esLCwi54PCEhgR1pBbzzowXJk6/szvrNp4mODOB0bgUAiYmJjBnahcLdR1i34WFCB/QgVMchE4mJiQCs22BbGGIvl5iYaHhhSPoZWHUI+gSaue+mh/Hw0J9tHI9mtRLStxs53+2z2zdHMnrRNHjuAHyZDc/8cigD5iYZHotKbjIa3PQK9OL+mx5yeTf/PAAbs2HpL4cy0MXdHD4DLxyCmEATCTMfdHk3zx+0vUO1+OdDGKyQm8aMkfEcKbO9O9UzwMSCmX9yaTdge3d3QxY8fdMQYucYdwP6n2foWDfHyuCfByE6wETizIUu72blQfg8E566KZYhbuUm0eXdrDoI67PgiasGM8zF3Zw4a3t3qru/iQdm3u/ybl48BJ9mwuNXDWK4ATdG6ZDJwUMPPURycjLDhw8nOlr/aXGrPz5st01pRiapS9Ya7pOjubYYGALPjru4GieTt9lv1A4Ze3h4QMIQsGgwIMSxGiq5GRACy93Izf1DoF6zfc85gkpu+ruZm/tioc4Kgy/8O4guVHLTL9h93ADcG2tb2xLroBtw7HnuCDd93MzNPYOhxgJDxA2glps/DIYqCwzr5HgNVdz0CnIvN3cPsl2WN/wi3Oihw3Yrmj59ekeVdlsGzZ2KyexFSfoprHUWPM0mrHVtr6JpzJSdyKW6oBSfsCBqSs5eoh5fPogbdRE36iJu1EXcqIu4UZfLxY3TtjIdMbAT29f8TFc7sB0WceP6p3XXD4uNAdCdaWzvVDQNn/AgfEIDydt+iIixg8izd2x4Q8aUXUhB6hF63XgFmV+lXpr+NiBu2s6Imw5G3NjNOA1xYzfjNMSN3YzTEDd2M07DRd0YxWmTA2+zqdXDzVrC5G126LAIVQ+YaIma4rPkbkmjx6R4QvqWUV1URkjfbroyPuHBoGlU5pVc8v3HxE3bGXHTsYgbdRE36iJu1EXcqIurujGKSx+C5m40rtAvzdC/C8WPV/Xn25vBCg4hbtRF3KiLuFEXcaMu4kZdLhc3MjkQWmT16tVUVlYyceJEqqqqyMrKIjY2lq1btzJlyhQ2bdrE4cOHue666xg37iJX+wiGEDfqIm7URdyoi7hRF3GjLh3pxmUmB5baOkPHaDdel2Y0Y/I2G+6bOzJnzpym2zk5OcTHxwPQv39/AGbNmgVAbm6uuLnEiBt1ETfqIm7URdyoi7hRFyNujOIyk4OSAydZN/Vh3e0bF7gYzTjz2rfoyWMI7hNFxltfEDlhKEExEWS8vdFw+/qK6nbtV1RUVKuPRUZG2vb1FTe62oub9kfctJ0RNxcibsSNvYy4uRBx47pujOIyk4PLgdOf72jxdnu1FxxH3KiLuFEXcaMu4kZdxI26XC5uPJ3dAVfGosFX2XCq5ZOxXY5T5bYPi9qL6HUhbtTF2uDmpNrbPOvG3dx8neM+bjIrbG7qrc7uycXT6OaEuFEOqwabcuC4m7jJEjfKkl15adx0yORg48aNJCUlcdttt/HWW2+1a+2e08YxaPaUDs/oYcleeHgn3PsDlNQ4VqPruMH0njmhwzP2KKqGe36AzfmwtPXTvttEJTdL99nc3PMDFLu4m5Ia2/fY5nzb95wjqOTmmUY3W6HQwXdWVXFTWgP3brW5WewGbpalwUM74I8/uIGb2nM/09zBzfIGN/f8APlVjtVQxU1Z7bmfaU/tcayGSm6e3X/OTZ6LuzlbZ/vZvDkfnnQDNysa3Pzxe8itdKyGKm7K6869bh7f3a6lL8BD09p/s9WVK1cyf/58XnzxRXr16sWNN97Y7PEVK1Zw8qT+BSgAfkU19N+QR68Z46k5U0HR3mPUlrb+Z+HDN0QAGM5Uhfvo7tP+IXMo6jIMT0stcTuW4lddZHg8AL1nTrhgq6uW+uVIRi+Vfl1IHb0Ai9mP8IK9DNn/hu6smm5mU9RlOJ6WWkbvXIZ/VYHh8YAabqr8OpMyegEWs797uIm9k6KuI/Cw1BG38xnXduPbiZS4hVjM/nQq2MfQ/at1Z1V0cyD2Dgq7jsTDUsfoXcsIqMw3PB5Qw021bxi7Ridi8Q4kvGAfQ1zezW8p7DoKD0sdo1KeJbAix/B4QBE3PqGkxC2k3juQToX7GZr2mu6s+m5WEFiRbXg8oIabGp8QdsUtpN47iLCi/Qzb59puDg6+nYKI0XhY6hmZuoKg8izD4wFF3HgHkxL/AHXeQYQVHWDYvld1ZwFiYmJISEjQ1bZD1hyUl5ezePFievTowYkTJy54XG/nzqdw9xHWbXgYzWolpG83cr5r+8/biYmJAIYzRha65FXBCwcht8qbVU/+VXcOzo0ndEAPQgdG6+qXIxkjrD0Gh87AHycNp6tfku6cim7yq+CFQ5Bd4c2LT/xFdw7UdXOw1OYmwm3cmN3CzXvHYX8J/GHSMCJd3E1BNaw8aHPz0uP6FwmCmm7ePw5pJXD39cOI8ndtN0XV8PxByKow8/I/HtSdAzXdfHAc9pXA/OuG0G2ua7sprobnmtz8SXcO1HTz4QnYUwy/v24I3ee4tpuSavjnQThV4cUrix7QnQN13ewuht9fF0sPA26M0iGTg4ceeojk5GSGDx9OdHTbT45RTiZvuyQZPUT4wWOjISnN8RqlGZmkLlnb4Rk93NLn4vIquenqB4+OEjeNiJuOc/Pr3kBvx/Mqueni615uZvW2fTiKSm7CfeERN3Lzq962D0dRyU0nN3vd/KKX7cNRVHIT5gt/FzeG6bDdiqZPn95RpZsYNHcqJrMXJemnsNZZ8DSbsNZZyN3S+ndBY6bsRC7VBaX4hAVRU3KWgpTDHd7fywlxoy7iRl3EjbqIG3URN+oiblwT197KVNPwCQ/CJzSQvO2HiBg7iDx7x1I3ZEzZhRSkHqHXjVeQ+VXqpenveYTFxjTt8aunraMZpyFu7Gachrixm3Ea4sZuxmmIG7sZpyFu7Gachrixm1ERl54c1BSfJXdLGj0mxRPSt4zqojJC+najMrv1hcGNGZ/wYNA0KvNKoP3XZNvF5G02fK2ZIxlnIW7URdyoi7hRF3GjLuJGXcSNa+LSk4PGVeClGZmGM43k25vBCg4hbtRF3KiLuFEXcaMu4kZdxI1r4tKTA9VYvXo1lZWVTJw4kaqqKrKysoiNjWXr1q3MmjWLN954g/z8fG699VbKysqIi4tzdpcvG8SNuogbdRE36iJu1EXcqIu40YfTJge1dRb2pBfbbTdiYCe8zaZL0KOLZ86cOU23c3JyiI+PB6B///4A3H333U2PZ506TeHuI7rqhsXGYPI2Y6mto+SAvvMhGjOCDXGjLuJGXcSNuogbdRE36iJu9OG0ycGe9GLG3vaJ3Xbb1/yMMUO7XHB/9OQx+HYOoTD1MKEDo/EO8qcyt5iKnCKK9x1vsZYjGUeJiopq83Gf4hrWTdW3j/iN65+m88h+lBw4aTjjDMSNvowzEDf6Ms5A3OjLOANxoy/jDMSNvowzEDf6MirispcVnf58R9NtvbM0RzKCccSNuogbdRE36iJu1EXcqIu4cV08nd0BV6a6Hl46BDmVzu5J+7At3zae6npn9+TicTs3BW7kxgIvu5GbHQ1uqtzFTTpkVzi7J+3DzgJ48aB7uKmxwCvpkOVmbirFjXLsLIRVB6HCDdzUWuDVDPdxs+sSuemQycHGjRtJSkqiZ8+e7N+/X3fulql9+NWkXnQO8221Tddxg+k9c4Kh/jiS0cMTe+CVDPihAPKrHKvRc9o4Bs2e0uEZe+RWwf+l2Mbz5F7Haqjk5sm959zkuribvCr4e4ObJ/Y4VkMlN0/tgZcb3Dg6QVDFTUH1udeNO7h5eq9tcvBDAWS7uJvCBjevHoZ/7HashkpuFu+FlxrcZDr4i44qboqq4e+p7uNmyT6bm60FcKrcsRqquCmpsf1/89phWLTbsRoquVm6D148ZHvdHD/rWA1V3JTWwiOpDW46+NgHD01r/81jV65cyfz583nppZeIj49n3LhxzR5fsWIF2/bm8c6Ovs3u79bVnwfuGMYDz5w7RvvWMUeJDKnGr6iG/hvyAOg9c8IFW139mMM3RAAYzlSF++gbJHBw8O0URIzGVF9N3I4l+NaU6s42jqfXjPHUnKmgaO8xaktb/qnS2C9HMnqp9gkjNW4Bdd5BdMnbxeCDbxseC6jk5jcURMThWV9N/I6l+NaU6M6q5yaUlNELqPcJpkteCoMP/tvwWEAdN4cG/Yb8yDg862uI27kUv2r7GxM0opqbGu8QUuISqfMJpnN+KrEH3jI8FlDJza3kR47Bs76auJ3PuLibYFLjEqn1CaFz/m5iD7xpeCygnhuP+hridi3Dv6pQd1Y1N7XeQaTEJVLrE0rn/D3EHviX4bGAOm7SB95CXtRYPOprGb1rGQFVBbqzyrkxB5ISt5Ba31DC8/cwxOXd3Exe1Dg8LLWM3rWcgMo83Vn13AQ0uAkjvGAvQ/a/oTsLEBMTQ0JCgq62HbLmoLy8nMWLF1NfX8/JkycvmBwkJCSwI62Ad360IHnyld1Zv/k00ZEBnM61/WkkMTGRMUO7ULj7COs2PEzogB6EDoy224fExEQAwxkji0Mq6uG1DBgZ5svVv/i77hzQNB7NaiWkbzdyvttnt1+OZIywKQdez4AXp8Xh76V/+y5V3byeAcPCfJn4i//TnQM13XybY/vr9IvTRhPgNdrwWFRyU9nwuhkW5sPEX/xNdw7UdPNdru2v7aumjiLQPMrwWFRyU9XgZkiYL9e6g5s82yVsq6aOJNA8UndOZTeDQ3247hf/T3cO1HSzOc92Od6qqSMINCcZHotKbqob3AwK8ea6X/xFdw7UdPN9HrzQ4CbIDdy8fhgGBHtz/c8f0p0DNd38kAcrD8GqqcMNuTFKh0wOHnroIZKTkxk+fDjR0fZlN7L648N225RmZJK6ZK2h/jiS0UOAF9wfe3E1TiZvs9+oHTJ6mBgFKUXg7+B3hWpu7nMjN1dHwc4i27gcQSU3/m7m5ieRsKMQAh3ckU4lN35ecK87uYmwrQkRNzZUcnNVBGx3Eze+XnCPG7m5MsJ2iVSQm7j54+CLq6GSm/ERtkukHHWjlw7brWj69OkdVdptGTR3KiazFyXpp7DWWfA0m7DWWXRlyk7kUl1Qik9YEDUlDl5YJ7SKuFEXcaMu4kZdxI26iBt1uVzcOG0r0xEDO7F9zc90tQPbYRE3rn9ad/2w2BgAhzJOQ9PwCQ/CJzSQvO2HiBg7iDx7x4Y3ZEzZhRSkHqHXjVeQ+VUHr1T5EeKm7Yy46WDEjd2M0xA3djNOQ9zYzTgNcWM34zRc1I1RnDY58DabWjzcrDVM3maHDotQ9YCJlqgpPkvuljR6TIonpG8Z1UVlhPTtpivjEx4MmkZlXgm0/xrzNhE3bWfETccibtRF3KiLuFEXcaMururGKC57CJo70rhKvzQj03CmkXx7M1jBIcSNuogbdRE36iJu1EXcqMvl4kYmB0KLrF69msrKSiZOnEhVVRVZWVnExsaydetWpkyZwn/+8x/Cw8Pp1q0bEya0/97EQuuIG3URN+oibtRF3KiLuFGXjnTjMpMDS22doaO0G69LM5oxeXfwEnAXYc6cOU23c3JyiI+PB6B///4A3H333QDk5uaKm0uMuFEXcaMu4kZdxI26iBt1MeLGKC4zOSg5cJJ1Ux/W3b5xgYvRjDOvfYuePIbgPlGc+N/3dBrSm9D+3Ul74RPD7TWrtV37FRUV1epjkZGRtn19xY2u9uKm/RE3bWfEzYWIG3FjLyNuLkTcuK4bo7jM5OBy4PTnO5puV2QWNvu8PdoLjiNu1EXcqIu4URdxoy7iRl0uFzeezu6AK1NUDY/vho/1vyOmNB+egG0FUFzt7J5cPMXV8MRu25jcgY9Puo+bEjdz898GN0Xu4KYGntgD/znh7J60D5+csrkpcAM3pTXw5B744Lize9I+/K/BTX6Vs3ty8Zyptb1u3nMTN+vcyE1ZLTy1B9475uyetA/Jp21ucjvYTbtPDjZu3EhSUhLDhg1j69at7V2entPGMWj2lA7P6OHpffDxKXjuAGRVOFaj67jB9J5pbKGIIxl7nC6Hfx6Ao2dhceunfbeJam4+OmUb0+lyx2qo4iazAlY0uHnajdw8dwBOurib7Mrz3Ox1rIZKbpbsg49O2tyccHE3uZXw3H73cvPhSXjuIBwtc6yGKm7yqmzfY+7k5qOT8PwBOHzGsRqquMmvOvcz7Sk3cLN0H/znJPzzIBxycTeF1bDiIn+m6cVD09p3s9WVK1cyf/58li5dyoQJE7j66qsvaLNixQpOnjT253a/ohr6b8ij14zx1JypoGjvMWpLW//f6/ANEQCGM1XhPrr7lDFgFrlR4zDVVTFmx2K86/T/b9o4HoDeMydcsNVVS/1yJKOXGu8gUkcvoNYnlMicHxiQ8YHurJpufkVu1BWY6qttbmr1n0aomptacyCpcQuo8QkjMmcrAzLe151V0k3/X5HbzeYmfscSfGr1/6ajnpsAUkcnUuMbRkTuNgamv6c7q6Kbw/1/SU638Q1uluJTq/9/U9Xc1JkDSBm9gBrfTkTkbmdg+ru6s0q66fdzcrpfiam+lridS/GtKTU8HlDEjZc/KXELqPENJyJ3BwPT1+rOqujmSL+byO4+wW3cpMYtoNo3nK55Oxl06B3dWSXd9J1Jdo+r8KyvJW7XM/hVlxgeD6jhpt7Ll5TRiVT7dTbsBiAmJoaEhARdbdt9zUF5eTmLFy9m6NChrU4A9HbufAp3H2HdhofRrFZC+nYj57u2/4SamJgIYDhjZKFLvRXWZ8KOwgAW/WqR7hycG0/ogB6EDozW1S9HMkY4fMY2s54640q8PK/UnVPZzfYCf/7xy8d050BNN0fOwMEzMHXGeLw8x+vOqerms0zYVuDPP375qO4cqOnmaBnsL4VpPa7Ay/MK3TkV3Vg0WH+60c0junOgpptjZZBWCtN6jMPLc5zunNpu/PjHL/+uOweKujkLaSUwrcdYvDzH6s6p6uazTPghz4/H3cDNibOwpximRY/B7DlGd059N/+nOwfu5cYo7T45eOihh0hOTmb48OFER7f9xDjCyeRtlySjBy9PmNETDjv4Fi/YDtJIXaL/ryaOZvTQP8T24SjipuPc9AuxfTiKam5u7AkZbuKmb7Dtw1FUcmPycC83fYJtH44ibjrQTZDtw1FUczM9GtIdvGwF1HLTK8j24SjiRl03eumQ3YqmT5/eEWUvYNDcqZjMXpSkn8JaZ8HTbMJaZyF3S5rdTNmJXKoLSvEJC6Km5CwFKYcvSZ8vF8SNuogbdRE36iJu1EXcqIu4cU1ceytTTcMnPAif0EDyth8iYuwg8uwdS92QMWUXUpB6hF43XkHmV6mXpr/nERYb07THr562jmachrixm3Ea4sZuxmmIG7sZpyFu7Gachrixm3Ea4sZuRkVcenJQU3yW3C1p9JgUT0jfMqqLygjp243K7CK7GZ/wYNA0KvNKoH3XZOvC5G02fK2ZIxlnIW7URdyoi7hRF3GjLuJGXcSNa+LSk4PGVeClGZmGM43k25vBCg4hbtRF3KiLuFEXcaMu4kZdxI1r4tKTA9VYvXo1lZWVTJw4kaqqKrKysoiNjWXr1q3MmjWLTZs2kZeXx5VXXklZWRlxcXHO7vJlg7hRF3GjLuJGXcSNuogbdRE3+nDa5KC2zsKe9GK77UYM7IS32XQJenTxzJkzp+l2Tk4O8fHxAPTv3x+AKVPOHfCRdeo0hbuP6K4d3K87ZUeydLcPi43B5G3W3d7dETfqIm7URdyoi7hRF5XcgPg5H3GjD6dNDvakFzP2tk/sttu+5meMGdrlgvujJ4/Bt3MIhamHCYuNwbdTMGdP5lGVX0phassr2h3JOEpUVFSbj/sU17Bu6sO661317D1sXrBSd/sb1z/ttOvexE3biJvWETfipjXETeuIG3HTFs7yI27s48zXTlu47GVFpz/f0XS75IC+05YdyQjGETfqIm7URdyoi7hRF3GjLuLGdfHsiKIbN24kKSmJadOmsWzZMt25W6b24VeTetE5zLfVNl3HDab3zAmG+uNIRg/HyuBP222n1TlKz2njGDR7iv2GF5nRw6pDtvEcc/CQHXHTcW5edCM3J866l5uXGtwccQM3J8vhwe2wu9jxzUFUcvNKus3NYQcPQFLJzalGN0Xu4ebVBjcZbuDmtJu5ea3BjaMHh6nkJrMCHtwBqW7i5vUMm5uDJe1euhkdMjlIT0/n/vvvZ+rUqQwYMEB37ttduYwfEUFhSXWrbfK3HTTcH0cyenj+IGzKtb2ATpQ7VsPT5MmZY9l4hwZ2aMYex87C20ds41np4MYAKrlZecg2lkNnbGNzBFXcnDgLbx+1jed5B58uldz8s+F1c7DU8V+oVXFzqhzePtbwunEDN88fgK9z4VCp46eLq+Ims8L9XjdfN/xMO+TgL22quMmqgH8fdaPXzXlu9pc6VkMVNzmV5/1MO+BYDZXcrDwIX+fYfk/b5+Av1Kq4ya069zPthQ7ewKlDLisqLy9n8eLF1NfXc9ttt13w+IoVK9i2Nw/o2+z+yVd2Z/3m00RHBnA6twKA5cuXExlSjV9RDf2B0AE9CB0YbbcPy5cvBzCcqQr3sduukaN9bsSj+1V41NewbNHTmOurdGcbx6NZrYT07UbOd/vabL/mnXfo2XBbT8boWOq8/PEYdR/4debo9k0sXJusO6uim2N9puPR/Sd41tew/B+LMddX6s6q58YPRt1vc7PjGxa+u053Vkk3vafh0eNqTJYaVjzu2m7qvXzRRiWAX2eO7fyWhe/+T3dWRTfHe0/Fo8c1eFpq+OeTSzDXVejOKufG5It19P3g14Xju75j4Xv217g1oqSbXlPwjJ6Ip6WGlU8txbtO/1+k1HPjg3XU/eDfleMpm1n43n91Z1V0cyJmMp49r8XTUsuLTy9Ryg0YG0+9yYf6BjcnUrew8P2PdeVAVTeT8Ox5HZ6WWl5eshTvWv1/LVTNjcXTm7rR94N/BCdSv2fh+x/pyjUSExNDQkKCrrYemtYxJ0skJyczfPhwoqNblr0jrcDQguTC3UcMLQxpPKHOaMbIwhBNs70lGuEPod66YwCGx3MpFoiV1sIze+EfceDhoT8nbsSNUTeHy6CrL4Tq/xkPiBsQN0Y4UwtL3cQN2H6mdfGFMHFjqG/Q8W4On4HOCroB4+Mpq4XFe+FxN3IT7gudLmM3RumwBcnTp0/vqNLK4OEBA0Od3Yv2I9QbOvl27DfcpULcqIuHBwwIcXYv2g9xoy4hbuQGxI3K9HcjN8Hetl+mxY16XCo3LrtbkTsyaO5UTGYvStJPYa2z4Gk2Ya2zkLslzW6m7EQu1QWl+IQFUVPi4EX2QquIG3URN+oibtRF3KhLe7opSGnf7T8vdy4XN06bHIwY2Inta36mqx3YDopofAtKD2GxMQAOZZyGpuETHoRPaCB52w8RMXYQefaODW/ImLILKUg9Qq8bryDzq9RL098GxE3bGXHTwYgbuxmnIW7sZpyGuLGbcRrixm7GabioG6M4bXLgbTa1eLhZa5i8zQ4dFKHi4RKtUVN8ltwtafSYFE9I3zKqi8oI6duNyuwiuxmf8GDQNCrzShzfr8tBxE3bGXHTsYgbdRE36iJu1EXcqIurujGKXFakEMf/uwWA0oxMw5lG8u3NYAWHEDfqIm7URdyoi7hRF3GjLpeLG5eZHFhq6wydltf41pPRjMnbbLhv7sjq1auprKxk4sSJVFVVkZWVRWxsLFu3bmXWrFm88cYbhIeHM2PGDPbu3kMf33DdtcXNxSFu1EXcqIu4URdxoy7iRl2MuNm/fz9xcXG6a7vM5KDkwEklt8hyV+bMmdN0Oycnh/j4eAD69+8PwN133930eEC5lXU/FzeXCnGjLuJGXcSNuogbdRE36mLETffu3Q3VdpnJgbsTPXkMvp1DKEw9jH9kJzrFxnB6w048vEytzqrbyrQnUVFRbT4e0Vn/2hFXRNyoi7hRF3GjLuJGXdrbjZG/yttD3Lium8jISEP1ZHKgCKc/39F0u+TASbJ0rGR3JCMYR9yoi7hRF3GjLuJGXcSNulxObjw7oujGjRtJSkpi8ODBfPSRseOd7dFz2jgGzZ7S4Rk97CyE32+Bx3c7vvC867jB9J45ocMz9rBq8I9U2JhtG5cjqOZm/hbbmKzu4GY3fOEmblIa3CxyAzeaZnv9b8yG7QWO1VDJTWqDm0fdxM2Te2xutuY7VkMlN3uKbf/fPJICFjdw81SDmx/cwM3eYtvrxp3cfJEN3+c5VkMlN/sa3PzdTdws3mtzsyW3XUtfQIdMDtLT07n//vu59dZbqaioaNfaniZPzhzLxjs0sEMzevj3UUgpgk8z4ZiD58Dkbzt4STL2OFIGn2ZBfjW8fdSxGiq5efso7CqyjelomWM1VHFz7Cysz4SCatv3nCOo5ObfjW4y4bCLuzlebnOTXw1rjjlWQyU3bx+zuVmfCemljtVQxc1JN3Pz7yO2/2/WZ8LBUsdqqOLmdMV5btzgZ9pbR8+9btJKHKuhipusSvgsy/b/jVu8bhrcfJZpm8Q5gipusitt4yiohneOt3v5ZnTIZUXl5eUsXryYgQMHtvj4ihUrOHnS2LVWfkU19Ac0q5WQvt3I+W5fm+2XL18OYDhTFe6ju0+Z0ddi7nE1HpqF5xYtwWSt1Z1tHE/ogB6EDoy2237NO+/Qs+G2nozRsdSbfPAZ8Qesfl3I3PIFC9/dpDurpJseE/GOvgY0K88/vhiTxXXdWEzeeI/4Ixb/LmR+v5GF736tO6umm2vwjp4ImpWVjy/By1KjO6ucG09vzCP/QL1/VzJ/+IqF736pO6ukm+5X493zWtCsrHrS1d2Y8RrxR0wBXcna+hUL33N1Nz/Bu+dPQbPy8pNL8bJU686q6MY04g+YAiLI3v41C9/fqDuropus7lc1uIHXnl6CV706bsDYeKyeXniO+COmgAhytm9i4ftf6MqBom66XYl3z+sBeH3xUsz1VbqzyrnxMOEx8o+YAiLJ3v4tC9//XFeukZiYGBISEnS19dC0jjmJITk5meHDhxMdbf/J0UPh7iNK7lZ0phZeSYc/DTMUMzyeq569h80LVupu78hY6q1QWQ/B3oZi4uYSuamohxA3cVNWCy+LG2XdvJQOD4obQ30DcWMEcWOjo92A8fG4m5uzdfDiocvbjVE6bEHy9OnTO6p0E4PmTsVk9qIk/RTWOgueZhPWOgu5W9LsZspO5FJdUIpPWBA1JWcpSDnsUB9CvMHTw9ERqIWXp/GJQWuIm/bFy7P9fhio4CZY3LSIKm5M4uYCxE37Im7Uxd3cBJnFjVE6ZM3BJUPT8AkPwic0kLLjOU3/6smYvL0oSD2Cl78PBalHLk1/LyfEjbqIG3URN+oibtRF3KiLuHFJXHor05ris+RuSaPHpHhC+pZRXVRGSN9uVGYX2c34hAeDplGZV+L4VkMXQVhsTNNbanoI7tedG3Vcv3Z+fWcibtqu70zETdv1nYm4abu+MxE3bdd3JuLG/tdwFuLG/tdQEZeeHBz/7xYASjMyDWcayd9+qF37pBeTt9nwdXOudCqguFEXcaMu4kZdxI26iBt1ETeuiWtfViQIgiAIgiAIQrvhtHcOauss7Em3v+nsiIGd8Da37/HsHcXq1auprKxk4sSJVFVVkZWVRWxsLFu3bmXWrFm88cYbFBYWcsstt1BcWEQf33BddcNiYzB5m7HU1uk+brsxI9gQN+oibtRF3KiLuFEXcaMu4kYfTpsc7EkvZuxtn9htt33NzxgztMsF90dPHoNv5xAKUw/jH9mJkH7dKdp3jPrKGor2tHyqiiMZI8yZM6fpdk5ODvHx8QD0798fgLvvvrvp8f0nClj3c31bZDVudVVy4KTubbUc2e6rvRA3+jLOQNzoyzgDcaMv4wzEjb6MMxA3+jLOQNzoy6iIy645OP35jqbbJQdOkvVVaodkHCUqKqrNxyM6XzjhcRfEjbqIG3URN+oibtRF3KiLuHFdZM2BIAiCIAiCIAhAB0wONm7cSFJSEmPHjuXNN980lL1lah9+NakXncN8W23Tddxges+cYKiuIxk9fJYJt3wNm3Jsp9Y5Qs9p4xg0e0qHZ+xRZ4UF22zj2aB/U4FmqORmw3lu6lzcTb0VEhvcfOYGbr7Ici83CxvcfOoGbr7Mhls3wdc5UGtxrIZKbh7YDjd/DcmnHashbjrOzZ8a3Kw75VgNldxsOs9NjYu7sWjwYIOb/7mBm29yz7mpdgM3f95h+//mv/qWNThMu08O0tPTuf/++5kxYwanThn7zvp2Vy7jR0RQWFLdapv8bQcN98mRjB4+y4QjZyG3Ck6UO1bD0+TJmWPZeIcGdmjGHsfOwg/5tvGsz3Kshkpu1mfZxpJTZRubI6ji5vhZ2NrgxtHJgVJuMs+5OVrmWA1V3Jwshx8KbOP53E3cHC6zuTns4m5OV8C2fDh6Fj53h59p57k5dMaxGqq4yay0/UxzFzfJWefcHCx1rIYqbrIqYGuBG7k5fc7N/hLHaqjiJqcStjX8f7Mhu93Ktki7rzkoLy9n8eLFeHt707t37xbbrFixgm1784C+ze6ffGV31m8+TXRkAKdzKwBYvnw5kSHV+BXV0B8IHdCDUB2HTCxfvhzAcKYq3Mduu0byIsfgF/1TAJ5/7Bk8Nf3T0sbxaFYrIX27kfPdPrv9ciSjF6uHiaBhv6MiIIq8TetY+O5Ow2NRyk1EPH49rwNg1SLXdxMw7C7KA7qR920yC9/bYT/0o7Go5Ca3axx+MdcDHrz4+DN4Wut1Z1V0Ezh0LmcDu5P/7acsfG+74bGo5Cav6+gGN568/MQyTNY63VkV3fgNnUt9YA8Kvv2Uhe9tMzwW9dxMQvPw4LUnXd2NJ35D76I+sAeFm9ez8P2thseikpv8LiPx63UDmoeJ1U8949JuNA9PfIfOpS4wmsLNn7Hw/R8Mj0UpN51H4Nd7MpqHiTefXsbb1lrdWRXdeA+dS01gNIVbPmfh+9/rzgLExMSQkJCgq62HprX/sXPJyckMHz6c6OjWRe9IKzC0W1Hh7iO6V4ADTafaGc1cqpXjRsbT2C9HMkZJSoOFQ41lxI24ETfiRtyIG6OIG3ED4sYojrgxSofsVjR9+vSOKOv2DJo7FZPZi5L0U1jrLHiaTVjr2v6Ld2Om7EQu1QWl+IQFUVPi4HU0QquIG3URN+oibtRF3KiLuFGXy8WN7FakEpqGT3gQPqGBlB3PafpXT8bk7UVB6hG8/H0oSD1yafp7OSFu1EXcqIu4URdxoy7iRl0uEzdOO+dgxMBObF/zM13twHaSXONbUHoIi40BcCjjLGqKz5K7JY0ek+IJ6VtGdVEZIX276cr4hAeDplGZVwLtf6VYm4ibtjPipmMRN/YzzkLc2M84C3FjP+MsxI39jLNwVTdGcdrkwNtsavHk49YweZsdujZL1dPnWuL4f7cAUJqhf9uTxkwj+dsPtWuf9CBu2s40Im46BnGjLuJGXcSNuogbdXFVN0aRy4oEQRAEQRAEQQCc+M6BUSy1dZQc0H/qQ+NbT0YzJm+z4b65I6tXr6ayspKJEydSVVVFVlYWsbGxbN26lSlTprBx40a8vLzo1q0bV4wZK24uIeJGXcSNuogbdRE36iJu1MWImwkTjB0w5zKTg5IDJ91qi6yWiJ48huA+UWS89QWRE4bi5e/Dif+2vo9ta+01q4PHzp7HnDlzmm7n5OQQHx8PQP/+/QG49dZbAcjNzRU3BtqLm/ZH3LSdETc2xE1zxE3bGXFjQ9w0x1XdGMVlJgeXA6c/39Hi7fZq7yhRUVGtPhYZGUlhrtqr7tsDcaMu4kZdxI26iBt1ETfq4qpujCJrDi6SGovyi851Y9XAcvGTWWWosdjG5A6IG3XRxI2yiBt1ETfq4m5uasWNYdptcrBx40aSkpKIjY0lLS2N//3vf3z44YftVb6JntPGMWj2lA7P6GHNUfjlVzB3M1TXO1aj67jB9J5p7FowRzL2qKqHuzbD/zJh7THHaqjkZu0xm5u7NtvG5giquKmuh981uFlz1LEaKrl5153cWBpeN6fhbQf/aKaSm/eO29zM+Q4qXdxNjQV+t8Xm5t9u4OaD89xUuIGbeVvgk9Pw5mHHaqjk5sMTNjezv4PyOsdqqOKm1gK/b3jdvOEGbj46aXNzpxu4qbPC/O9tbl530I1ePDStff7uvXLlSubPn8+rr77KVVddxalTpwCYNm3aBW1XrFjByZP6F6AA+BXV0H9DHr1mjKfmTAVFe49RW1reavvDN0QAGM5Uhfvo7lPa0Lso7jwErBbidi0joEL/dV2N4wHoPXPCBVtdtdQvRzJ6KQ/oRkrcAvD0olNhGkPTXtedVd3N6F3LCazINjweUMNNRUAUKaMXoJnMdCrcz9C013Rn1XQzl+LOQ21uUp4lsDzL8HhAETf+EaTEJaKZvN3Czf4hcyjqMgw0KyN3PUtwuf7t+lRzU+nflZS4RKwmH8KKDjBs36u6s2q6mU1Rl+E2NynPEXz2lOHxgCJu/LrY3Hj5ElZ0kGH7XtGdVduNxojU5wgp0//7jWpuqvw6sytuIVYvX0KLDjLc1d3E3klR1xE2N7ufJ+TMccPjAUXc+HYiJW4hFrM/oUWHGL7vZd1ZgJiYGBISEnS1bbc1B+Xl5SxevJj4+HgyMjIwm83U1ta22FZv586ncPcR1m14GM1qJaRvN3K+29dm+8TERADDGSMLXbbkwquHoazOxAuL/ozJQ3e0aTyhA3oQOjBaV78cyejFosEjKZBdCb8bN5Qr5ybpzqro5vs8eDUDztSZWLXoTy7v5tFUyKyA340dwgQXd/NDPryS3uDmsQdc2o21wc3pCpg7dgg/cXE32/LhpXQorfXkxUcX4mXgvWUV3TyWCqcqYPbYWK6Z49puthfAi4ca3DyywG3c3Bk/mIku7mZnIbxwEEprPXjp7wku7UbTYNFuOFEOv40fzE9d3M2uQlh5CEprPHjp/+5zeTf/2A3Hy+H2+EFcZ8CNUdptcvDQQw+RnJxMbGws0dFtPyEXw8nkbZcko4cJkbaPpDQM/YJzPqUZmaQuWdvhGXuYPODxuIuroZKbKyNsH+7i5h+jL66GSm7Gd7V9uIMbTw9Y5EZuxnW1fSSlYeg/0fNRyc1jbuRmbBfbh7ixoZKb+M7w+k/cw42HBzwy6uJqqOQmrjO8fpX7uPn7RbrRS7vuVjR9+vT2LGeXQXOnYjJ7UZJ+CmudBU+zCWudhdwtaXYzZSdyqS4oxScsiJqSsxSkdPAFXJcZ4kZdxI26iBt1ETfqIm7URdy4Jq69W5Gm4RMehE9oIGXHc5r+1ZMxeXtRkHoEL38fClLdf/utS464URdxoy7iRl3EjbqIG3URNy6JS59zUFN8ltwtafSYFE9I3zKqi8oI6duNyuwiuxmf8GDQNCrzSpyyF2lYbEzTASB62jqacRbixn7GWYgb+xlnIW7sZ5yFuLGfcRbixn7GWYgb+xkVcenJQeMq8NIM/Ttq/HjleP72Q+3aJ72YvM2GT/lzJOMsxI26iBt1ETfqIm7URdyoi7hxTVz7siJBEARBEARBENoNp71zUFtnYU96sd12IwZ2wttsugQ9unhWr15NZWUlEydOpKqqiqysLGJjY9m6dSuzZs3iww8/pKKigilTppCfn09c3EVuDyToRtyoi7hRF3GjLuJGXcSNuogbfThtcrAnvZixt31it932NT9jzNAuF9wfPXkMvp1DKEw9TOjAaIJ7R5H7w37qK2so2tPyMbKOZIwwZ86cpts5OTnEx8cD0L9/fwBuu+22psfNZvNFfz1VETfqIm7URdyoi7hRF3GjLuLGdXHZNQenP9/RdLvkgL7TCB3JOEpUVFSbj0dGRnbo13cm4kZdxI26iBt1ETfqIm7URdy4LrLm4CI5VQ51Vmf3QmiJ0+JGWcSNumRWQK24URJxoy7iRl0yK6DW4uxeuBbtOjnYuHEjSUlJxMbGkpaWxgcffEBaWusHXfyYW6b24VeTetE5zLfVNl3HDab3zAmG+uVIRg8vHITffgsbsqC8rt3LuxwquVl1EG5vcHNW3Cjl5qVD59yU1bZ7eZdDJTevpMPt39jcnBE3Srl59Tw3peJGKTevZ9jcfJ4FpTXtXt7lUMnN6sMNbrKhWNzoxkPT2m/z2JUrVzJ//nxeffVVrrrqKgIDAykvL2fo0KHN2q1YsYJte/N4Z0ffZvd36+rPA3cM44Fnzh2jfeuYo0SGVONXVEP/DXkA9J454YKtrn7M4RsiAAxnqsJ99A0WSBs6l+LOQ/Gw1BG3axn+lfm6s+6Eym6w1BOXkkRARa7urDuhpps5FHceBlYLo3clEVhh50AcN0VFN/uHzKaoy3CwWhiVspyg8mzdWXdCSTexsynqanMzMuU5gstP6866E2q6uZOiriNsblL/SfDZU7qz7oSKbg7E3kFh15FgtTBi90pCyk7ozrobMTExJCQk6GrbrmsOysvLWbx4MfHx8WRkZLQ6OUhISGBHWgHv/GhB8uQru7N+82miIwM4nVsBQGJiImOGdqFw9xHWbXiY0AE9CB0YbbcviYmJAIYzRvanPVgCL6ZDnyAz99/0MB4euqNuhYpu0kvhhUPQK8iLBTf9Wdwo5CbjDKw8CL0CTSyY+aC4UcjN4TPw/EGICTSROPNP4kYhN0fKbG56+Jt4YGaiuFHIzbEyeO4A9Agw8cDMBeJGQTfd/U38aeb9l60bo7Tr5OChhx4iOTmZ2NhYoqPtS/4xqz8+bLdNaUYmqUvWGqrrSEYPg8NgxRXtXtZlUcnNwFBxcz4quRkQIm7ORyU3/cVNM1Ry0y8Ynh3X7mVdFpXc9AmGZ+V104S4cX3afbei6dOnt3dJQRAEQRAEQRAuAe265sAIRg9Bs9TWGdrWKiw2BjC2FVZYbAwm78t3X1tHETfqIm7URdyoi7hRF3GjLuLGfXDa5EAQBEEQBEEQBLWQcw4EQRAEQRAEQQBkciAIgiAIgiAIQgMyORAEQRAEQRAEAZDJgSAIgiAIgiAIDcjkQBAEQRAEQRAEQCYHgiAIgiAIgiA0IJMDQRAEQRAEQRAAmRwIgiAIgiAIgtCATA4EQRAEQRAEQQBkciAIgiAIgiAIQgMyORAEQRAEQRAEAZDJgSAIgiAIgiAIDfx/1QVZjGQc2mgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 980.827x207.356 with 1 Axes>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit import QuantumCircuit\n",
    "from qiskit.synthesis import LieTrotter\n",
    "from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit\n",
    "\n",
    "total_time = 0.6\n",
    "num_steps = 7\n",
    "\n",
    "# Generate an initial state\n",
    "initial_state = QuantumCircuit(L)\n",
    "for i in range(L):\n",
    "    if i%2:\n",
    "        initial_state.x(i)\n",
    "\n",
    "trotter_circuit = initial_state.compose(\n",
    "    generate_time_evolution_circuit(\n",
    "        hamiltonian,\n",
    "        time=total_time,\n",
    "        synthesis=LieTrotter(reps=num_steps),\n",
    "    ) \n",
    ")\n",
    "\n",
    "trotter_circuit.draw(\"mpl\", fold = -1, scale = 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Optimize the problem\n",
    "\n",
    "For the noiseless simulations we won't need to transpile, but for the hardware expeirments we will. Let's go ahead and transpile our PUBs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create PUBs (Primitives Unified Blocs) which specify our experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pub_circuits    = [trotter_circuit]\n",
    "pub_observables = [observable]\n",
    "pubs = list(zip(pub_circuits, pub_observables))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transpile PUBs (circuits and observables) to match the backend ISA (Instruction Set Architecture)\n",
    "By selecting `optimization_level=3`, the transpiler will choose a 1D chain of qubits which minimizes the noise affecting our circuit. Once we have converted our circuits into the format that the backend is prepared to accept, we will apply a complimentary transformation to our observables as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_ibm_runtime import QiskitRuntimeService\n",
    "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
    "\n",
    "# need to create a service that's not local\n",
    "service = QiskitRuntimeService(channel=\"ibm_quantum\", instance='ibm-q-internal/dev-services/middleware')\n",
    "backend = service.backend(\"ibm_torino\")\n",
    "\n",
    "# Transpile pubs (circuits and observables) to match ISA\n",
    "pass_manager = generate_preset_pass_manager(backend=backend, optimization_level=3)\n",
    "isa_circuits = [pass_manager.run(circuit) for circuit in pub_circuits]\n",
    "isa_observables = [obs.apply_layout(isa_circuit.layout) \n",
    "                                if isinstance(obs, SparsePauliOp) \n",
    "                                else [sub_obs.apply_layout(isa_circuit.layout) for sub_obs in obs]\n",
    "                                for (isa_circuit, obs) in zip(isa_circuits, pub_observables)]\n",
    "isa_pubs = zip(isa_circuits, isa_observables)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Execute quantum experiments using `AerSimulator` backend\n",
    "\n",
    "We can simulate our various Trotter circuits, including the back-propagated one using the ``Estimator`` primitive from ``qiskit-ibm-runtime`` and we will compare the results we get with a Trotter circuit of 100 steps. \n",
    "These simulations do not have any sources of noise apart from statistical uncertainty due to drawing finitely many shots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noiseless local simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Qiskit Runtime service and Estimator\n",
    "# V2 is the default version in qiskit-ibm-runtime>=0.28.0, no need to import EstimatorV2 as Estimator\n",
    "from qiskit_ibm_runtime import QiskitRuntimeService, Estimator\n",
    "\n",
    "num_shots = 2**20\n",
    "\n",
    "# Instantiate a local Estimator using a Simulator\n",
    "local_service = QiskitRuntimeService(channel=\"local\")\n",
    "simulator = local_service.backend(\"aer_simulator\")\n",
    "simulator_estimator = Estimator(simulator)\n",
    "\n",
    "# Simulate experiments, retrieve epxectation values\n",
    "noiseless_result = simulator_estimator.run(pubs, precision=1/np.sqrt(num_shots)).result()\n",
    "noiseless_expval = [pub_result.data.evs.tolist() for pub_result in noiseless_result]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expectation Value: -0.671061\n"
     ]
    }
   ],
   "source": [
    "print(f\"Expectation Value: {noiseless_expval[0]:3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}