{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use QuTiP operators in graphs\n",
    "**Incorporate QuTiP objects and programming syntax directly into graphs**\n",
    "\n",
    "It is possible to use Boulder Opal with [QuTiP](https://qutip.org/) operators.\n",
    "This is possible in [graph](https://docs.q-ctrl.com/boulder-opal/toolkit/design/calculate-with-graphs/how-to-represent-quantum-systems-using-graphs) execution, where every operator can be replaced by a [QuTiP `Qobj`](https://qutip.org/docs/latest/apidoc/classes.html#qobj).\n",
    "\n",
    "## Example: QuTiP operators in a graph-based simulation of a single qubit with leakage\n",
    "\n",
    "In this example we to simulate the dynamics of a single qubit subject to leakage to an additional level. The resulting qutrit is treated as an oscillator (truncated to three levels) with an anharmonicity of $\\chi$, described by the Hamiltonian:\n",
    "$$\n",
    "H(t) = \\frac{\\chi}{2} (a^\\dagger)^2 a^2 + \\frac{\\Omega(t)}{2} a + \\frac{\\Omega^*(t)}{2} a^\\dagger,\n",
    "$$\n",
    "where $a = |0\\rangle \\langle 1| + \\sqrt{2} |1\\rangle \\langle 2|$ is the lowering operator and $\\Omega(t)$ is a time-dependent Rabi rate.\n",
    "\n",
    "In the example code below we illustrate how the noise-free infidelity $\\mathcal{I}_0$ can be extracted from the coherent simulation results. To do this, we choose as target an X gate between the states $|0\\rangle$ and $|1\\rangle$. Notice that this target is not unitary in the total Hilbert space, but is still a valid target because it is a partial isometry—in other words, it is unitary in the subspace whose basis is $\\{ |0\\rangle, |1\\rangle \\}$.\n",
    "\n",
    "Note the use of QuTiP to define operators where convenient in defining the relevant graph nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import qctrlopencontrols as oc\n",
    "import qctrlvisualizer as qv\n",
    "import boulderopal as bo\n",
    "\n",
    "plt.style.use(qv.get_qctrl_style())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your task (action_id=\"1829175\") has completed.\n",
      "Noise-free infidelity at end: 8.7e-02\n",
      "Time evolution operator at end:\n",
      "[[ 0.057-0.145j  0.235-0.926j -0.088+0.235j]\n",
      " [ 0.235-0.926j -0.067+0.187j -0.193+0.102j]\n",
      " [-0.088+0.235j -0.193+0.102j -0.825+0.456j]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAFUCAYAAABP8bodAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAAsTAAALEwEAmpwYAABquUlEQVR4nO3dd3hT9dvH8Xd22qTpoqWlbGRP2RsFRZYTwb1R3HtvH7f+3AsVN4gg4GSogAzZQ9l7j+6VJs0+5/kjpVCJzJYk7f26yJX05CS521OaT77raFRVVRFCCCGEEFFJG+4ChBBCCCHEyZMwJ4QQQggRxSTMCSGEEEJEMQlzQgghhBBRTMKcEEIIIUQUkzAnhBBCCBHF9OEuIFy27s3FoNeFuwwhhBBCiGPy+QM0rZcS8r4aG+YMeh0N05PCXYYQQgghxDHtyiz4z/ukm1UIIYQQIopJmBNCCCGEiGIS5oQQQgghopiEOSGEEEKIKCZhTgghhBAiikmYE0IIIYSIYhLmhBBCCCGiWMSuMzfvjwUsXbCczL2ZdOzekWtGX/Gf+86ZMY9Z0+bg83jp0LU9I6+/FIMhYr81IYQQQohKE7Etc/EJ8Zx3wbl079vtqPttXLOJWb/O5q5Hb+O5t58iLyef6VNnnqYqhRBCCCHCK2Kbrzp0aQfA3p178RYU/+d+S/9aTvd+3UivmwbAoIvO5auPxnPhZcNOS53/ZdkkN7PeLUWrA40WtDoNWh1o9WAwazCYNRhjwWjWYIjVEBOnISZeS2y8BnO8hth4LZYkDXEpWuJqadEZNGH9foQQQggRmSI2zB2vzH1ZtO3YpvzrjPp1KCkuwVnixBJnqbDvwjmLWTh3MQBDrroUqvB0XvZshb2r/ZX2fJbEYLCz1daSkKElMUNH4sHrulpSGukwx0VsQ6sQQgghqkjUhzmvx0tMrLn865iYGADcbs8RYa5X/x706t8DOPo5zipD18vNNO1tQA2AEgAloKIooPjB51bxuVW8pSpeF3hLVdwOFVexgqtIpbRYobRIxVmg4MhTKclTcBaqOAsDZG0J/OdrWpM11Gqko1YjHamNddRupie9uY7aTfUYY6VlTwghhKiOoj7MGU1G3C5P+ddulxsAs9kUrpIAsKVosaVUTkuZElApLVSx5yoUZykU7Q9QuF+hsOy6YE+AvF0BHPkqjnw/u1Yc2SKYVF9Legs9Ga30ZLTVU7etntpNdej0EvKEEEKIaBb1YS69bhr79xygY7cOAOzbc4C4+LgjWuWimVanwVpLg7WWljotQ++jKCr2LIXcnQHydgbI2R5sxcva7Cdne4CCPQoFe7ys/91b/hi9CdJb6mnQQU+DjgYadNJTp6VexucJIYQQUSRiw1wgEEAJKCiKgqoq+Lw+tDotOp2uwn5de3dm3Cff0blnR+IT4vntpz/o1qdLmKoOH61WQ0IdHQl1dDTtVfG+gE8ld2eAzI1+9q/3s2+tn33r/OTvUtj7j5+9//j568tgi6bBDHXb6mnc1UDjHgaadDcQX1sX4hWFEEIIEQk0qqqq4S4ilOlTZzLjh98rbBt88UC69+3Gi4++yhOvPEJSrUQA5syYy6xf5+Dz+mjfpR2X3TDimOvM7cosoGEVToCIBq5ihX1r/ez+28/uVT52r/KTu+PIMXnJDbU06W6gaU8jzfoaSGmsQ6OR1jshhBDidDlabonYMFfVJMyF5ixQ2L3Kx/alPnYs9bFzuR+Po+KvSGKGlmZ9jTTrY6DFWUaS6knLnRBCCFGVJMyFIGHu+AT8KgfW+9m22MfWv3xsWeDFWVDxV6Z2Ux0tBxhpNcBI094GzFZZIkUIIYSoTBLmQpAwd3IUReXAhgBb5nvZMt/L5vk+3CWHfoV0Bjijp4G2g0y0HWQk9YyIHZYphBBCRA0JcyFImKscAZ/KzuU+Ns7xsmG2l90r/Rz+G1W7qY62g4y0HWLijB4GtDoZayeEEEKcKAlzIUiYqxqOfIUNs72snelh/R9eXEWHfr2stTS0G2LizAtMND/LiMEkwU4IIYQ4HhLmQpAwV/UCfpXtS3ysne5h9XQvudsPzZQ1x2loM8hIp4vNtD7XiMEswU4IIYT4LxLmQpAwd3qpqsqBjQH++dnDPz972Lf20FkqzDYN7Yea6DzcRIuzjeiNEuyEEEKIw0mYC0HCXHjl7vCz6kcPK6d62Lv6ULCzJGroeLGJbpebadzdIOvZCSGEEEiYC0nCXOTI3upn5VQPK6a4ydx4qCu2ViMtXS8z0+0ys8yKFUIIUaNJmAtBwlxk2r/ez9IJbpZNclOcqZRvb9LdQM9rzHS82IQ5TtaxE0IIUbNImAtBwlxkUwIqWxb4WDrBzd8/efA4g7+mJouGjheZ6HGNmTN6SjesEEKImkHCXAgS5qKH26Hw948eFn3jZtsiX/n22k119L4hhu5XmrEmS2udEEKI6kvCXAgS5qJT9lY/i8e7WTLeTXFWsBtWb4KOF5nofUOMtNYJIYSoliTMhSBhLroF/CprZ3hZ8IWLjbO85WedSG+p46zRsXS9zCTniBVCCFFtSJgLQcJc9ZG3K8DCr1ws+tqNPSfYWme2aehxtZl+o2Ko3VRmwgohhIhuEuZCkDBX/fi9Kn//7GHexy62Lzk0tq7VOUYG3BFDywFG6YIVQohqSlVVVJeLgN0evBQXE3A4UZylKKVOlNLS4G2XC8XjQfV4UL3e4G2vFzUQgEAA1R8ov13u4HuHBjQ6Hej1aHR6NHodGp0OQ3o66Y89UqXf39FyizRZiGpDb9TQ5VIzXS41s3e1j7mfuFg+yc2GWV42zPKS3kJH/9tj6Xq5GWOMhDohhIgWAYcTf04OvpwcfNk5BPLz8RcUlF8CBYX4i4rA7z/mc1VJfSWOsLzuQdIyJ6o1R77CX1+4mPuJq3zdOmuyhj43xXDWrbHYUmRcnRBChJuqKPhzcvHu34/vwAG8+w8Erw8cwJ+djeIsPa7n0ZjN6Gxx6Gw2dDYbWqsVndWCNtaCNjYWrSUWbUwMGpMJrcmExmgsv9YY9KDTodEGW9vQaQENoFI+MFsFNXCw5c6P6vej+gNoDAZiWraosp8PSDdrSBLmaha/V2XVDx5mf1DKnr+Dn9z0JuhxlZlz7oqVM0wIIcRpoKoqvswsPDt24N29G8/uPXh378G7dy+q2/2fj9OYzRhSU9Cnpgava9VCn5SELjkZfVIi+qRkdIkJaI3G0/jdnF4S5kKQMFczqarKtkU+Zr1byprpXiA4FKL9+SYG3htLoy6GMFcohBDVg+r349m1C/fmLXi2bcezYwee7TtQSkO3sukSEzHWrYuhTjrGjDoYMjIw1knHkJaO1hZX48c8S5gLQcKcyNzkZ9Z7pSz7zo0/mOto3s/AoActNO8n69UJIcTxUlUVf3Y2rvUbcG/ajGvTJjzbtqN6PEfsq0tMxNSkMaaGDTE2qI+xQX1M9eujs9nCUHn0kDAXgoQ5cVBxVoA5H7qY/5kLtz3436FhZz2DHrTQdrARrVZCnRBCHE4NBPDs3Ilr7Tpc69bjWrsOf17eEfsZ6tTB3LwZpqZnYG7SGFOTJuiT5L33ZEiYC0HCnPi30iKFeZ+6mPNBKY784H+LOq10DH7YQseLTGh1EuqEEDWTqqp4d+2m9O9/gpfVq1EcFWdwauPiiGndCnPLFsS0aI65WXN08dLaVlkkzIUgYU78F2+pyl9fuZj1TimF+4MzYNOaB0Nd5+ES6oQQNYO/sBDn8hU4l6+gdNUqAoVFFe43pKUR07ZN8NKmDcb69dBoZYWAqiJhLgQJc+JY/F6VJePdzHzDSf7uYKhLPUPHkIdj6TzCjE4voU4IUX2oioJ70yacS5bhXL4c9+YtFe7XJSdhOfNMYs/sQOyZHTCkpYWp0ppJwlwIEubE8Qr4VJZ862bG/5zk7wqGutpNdQx9zEKn4SYZUyeEiFqKx0Ppqr9xLFyEY/ESAoWF5fdpDAZiOrTH0qULli6dMNavLxPDwkjCXAgS5sSJCvhUln7nZsbrTvJ2BkNdeksdw56w0OF8CXVCiOiguFw4liylZN58nMuWV1jfTV+7NtYe3bF07UJsh/ZozeYwVioOJ2EuBAlz4mQFfCqLx7uZ8ZqTgr3BUFe3nZ4LnrLQ5jw5/6sQIvJUCHBLl1VYMsTUtCnWXj2w9uqJqXFj+RsWoSTMhSBhTpwqn0dl0dcuZrxeWn6qsCbdDVz4rIWmvarvKuRCiOigBgI4V6zEPms2joWLKrTAmVu3wnZWP6x9emNITQ1jleJ4SZgLQcKcqCw+t8r8sS5m/s9ZvqRJ64FGLnzaQr32ckYJIcTpo6oqni1bKP5jFiV/zq0wA9XcqiVxZ51FXF8JcNFIwlwIEuZEZXPZFWa/X8qs91x4HMH/Vp0vNXHB01ZSGunCXJ0QojrzFxdj/30WxTNm4t21q3y7sV49bOcOIG5Af4zp6eErUJwyCXMhSJgTVaUkV+G3N5zMG+vC7wGdAfqOimHwQxbiUmQNJiFE5VADAUpXrqJo+gwcixaD3w+ALj4e24D+2M4dgKlZMxkDV01ImAtBwpyoavl7AvzygpNl37lRVTDHaRh4XywD7ojFGCt/XIUQJ8dfWEjxjJkU/TINf3Z2cKNWi6VrF+IHnYe1R3c0BhniUd1ImAtBwpw4Xfat9fHDM042/OEFID5dy4VPW+h2pVmWMxFCHBdVVXGtXUfRz79QMn9BeSucIT2N+CGDsQ0ciCGlVpirFFVJwlwIEubE6bZpnpcfnnKw5+/gH+G67fQMf8lKi34y81UIEZri9VIyew6FU37As2NHcKNWi7V7NxIuOJ/Yzp3kFFo1hIS5ECTMiXBQFJXlkzz89Kyj/LyvbQcZufh5K+kt9GGuTggRKfz5+RT9/AtFv0wjUFQEgC4xgfghQ0gYNhRDbZmNWtNImAtBwpwIJ69LZfb7pfz2Zikeh4pWB31vjmHYYxYsSfIpW4iayrNzJwWTvsc++8/yrlTTGU1IHH4JcWefhdYoLfk1lYS5ECTMiUhgz1H45QUHC79yoypgSdQw7AkLfW6KQaeX8XRC1ASqquJas4aC7ybhXLosuFGrxdqzB4nDLyGmXVuZkSokzIUiYU5Ekn3r/Hz/SAlb5vsASG+h49JXrLQaYApzZUKIqqIqCo6FiyiYMBH3pk0AaEwm4gedR+KISzHWkXXhxCES5kKQMCcijaqqrP7Vy9QnHeTuCADQfqiR4S/HyaLDQlQjaiBAydx55I+fUL7Ar85mI+GiC0m46AL0CQlhrU9EJglzIUiYE5HK51H586NSpr8aHE+nN8E5d8cy6AELJot0tQgRrVS/H/us2eR/OwHfvv0A6FNSSLpsJPFDBqE1m8NcoYhkURnmnA4n346dyKa1W7DEWbhg5BA69+x0xH4+n58p435gzYq1BAIBGjdtxGU3XEpCUsJRn1/CnIh0xVkBfnjaydIJwZNjJ9TRcskLVjpfapLxM0JEETUQwP77H+SPG48vMwsAQ3o6SVdeTvzAc2WBX3FcojLMffHBN6iqylWjLmPf7v2MeWMs9z99N+l10yrsN+vXOSxftJI7HrmVmBgzEz7/Ho/Hw8333HDU55cwJ6LFjqU+Jj5UUr4+XdPeBi5/M446LWUpEyEimRoIYJ89h/yvx+E7cAAInis16aorsA3oj0YnwyfE8TtabonINRA8bg+rl69h2PBBmMwmmjRvTNuOrVm2cMUR++bnFtCybQts8XEYjAY6du9A1r6sMFQtRNVo3M3AI3MTufr9OKzJGrb+5ePFngVMebwEd4kS7vKEEP+iKgr22XPYdePNZL3yGr4DBzDUzSD9icdo+PmnwdY4CXKiEkXkR/ucrFy0Oi2p6YcWRcyoV4dtm7YfsW+Pft2YPO4HiguLiYmNYcWilbRq3+J0litEldNqNfS6LoYO55v4+QUnC8a6mPWei+Xfe7jkRStdRkjXqxDhpqoqzmXLyRv7GZ7twbM1GNLTSb72amznDJAAJ6pMRIY5j8eLOabiQFBzrBm323PEvilptUhMSuDJu59Dq9VSp146I669JOTzLpyzmIVzFwMw5KpLQbpZRZSxJGm54s04el5j5rv7S9i1ws8XN9lZ9LWBy9+II615RP6XFqLac63fQO6nY3GtWQsEJzYkX3M18YMGotHL/0tRtSLyN8xkMuJ2uStsc7s8mM1Hrrk16asp+P1+XvnoeYwmE7OnzeGj1z/lwefuPWLfXv170Kt/DyDY9yxEtGpwpoGHZiey+Bs3PzztYPM8Hy/0KODce2MZ/JAFY4y00glxOnh27yZv7Oc4Fi4CQGuLI/nKK0i48AK0JlknUpweETlmLjUtBSWgkJOVW75t/54DpP1r8gPA/t0H6NanKxarBYNBT99z+7B7xx4cJY7TWbIQp93BrtdnVyXT6zozAR/MfL2U/+uaz9qZR7ZiCyEqj7+gkKy33mHXTbfgWLgIjdlM8tVX0XjcNySNHCFBTpxWERnmTGYT7Tu3ZdqUmXjcHnZs2cnaVevo2qvzEfvWb1yPZX8tx1XqIuAPsGD2QuITbVjjrGGoXIjTz5qs5er3bTz4RyIZbfTk71L4cEQxn1xdTNGBQLjLE6JaUdxu8seNZ8c111H8y68AxJ8/lMbjvqLWjdejs1rCXKGoiSJ2aRKnw8n4Tyeyed0WLHGxXDByKJ17dmLb5h189PonvDH2leB+JU4mf/MDm9ZtIRDwk143nYuvvICGTRoc9fllaRJRHQX8Kn9+5OLXF514nCrmOA0XPmOh76gYtDrpehXiZKmKgn3WbPLGfo4/Lw8AS/fupNwyClPDo7/fCFEZonKduaomYU5UZwV7A0x8sIQ1070ANOik56p346jXThYnFeJEuTZsIOeDj3BvDJ4/1XTGGaTceguWjmeGuTJRk0iYC0HCnKgJ/vnFw8QHSyg6oKDVQf/bYxj2hFVOCybEcfDl5pH36Vjss2YDoEtKImXUjdgGnotGG5GjlEQ1JmEuBAlzoqZwlyj8/IKTuWNcqAokN9Ry5dtxtBogA7SFCEXxein8fjL54yegut1oDAYSRwwn+cor0MbGhrs8UUNJmAtBwpyoaXat9DH+rhL2rQ2eFqzr5SYufSmOuBRpYRDiIMfSZeS8/wG+/cHTb1n79CZl9C0Y66SHuTJR00mYC0HCnKiJAj6VWe+VMu1lJz43WJI0jHjFStfLzXIGCVGjeQ9kkvvhRzgWBReWNzaoT+pdd2Dp2DHMlQkRJGEuBAlzoibL3eFn/D0lbJ7rA6DVuUaufDuO5PpyuiFRsyheLwXfTaRg/ARUnw9NTAy1rruGxEsuljM3iIgiYS4ECXOiplNVlSXj3Ux+3EFpoYrJquGi58qWMdFKK52o/pwrV5H9zrv49u0HwHbOAFJG34w+OTnMlQlxJAlzIUiYEyKoODvApAcdrPoxeNaIJt0NXP2+nOdVVF/+/HxyPhxDyZ9zgWCXau177ya2ffvwFibEUUiYC0HCnBAV/f2Tm+8ecGDPVtCbYNjjFs65OxadXlrpRPWgBgIU/fIreZ99juIsRWMykXzN1SSNGI7GIGswisgmYS4ECXNCHMlZqDDlcQeLx7kBaNBRzzUf2shoLa10Irp5duwk6823cG/YCIClezdS77oDY7rMUhXRQcJcCBLmhPhvG2Z7GH9XCQV7FXQGGPywhUEPxKIzSCudiC6K10v+N+Mp+G4iBALokpOofdedWPv0lhncIqpImAtBwpwQR+eyK/z4jJP5Y10A1G2r59oxckowET1K//6HrLfeLp/gkHDB+dQadRM6qyXMlQlx4iTMhSBhTojjs2WBl2/usJO3U0Grh8EPxTLoQQt6o7RqiMgUcDjJ/fgTiqdNB8DYoAG1H7iX2DZtwlyZECdPwlwIEuaEOH4ep8qPzziY+7G00onI5li0mOy33sGfn4/GYCDp6itJvvwymeAgop6EuRAkzAlx4rb85eWb26WVTkQef2EhOe9/WL7ciLlVS9IefABTwwbhLUyISiJhLgQJc0KcnH+30tVrr+e6j2XGqwgPVVUpmTuPnHfeI2C3ozGbSbnpBhIuuhCNTs5oIqoPCXMhSJgT4tRsWeDl69vs5O8Ozngd9riFc++VdenE6eMvKCT7nXdxLPgLgNiOZ1L7gftkuRFRLUmYC0HCnBCnzl2iMPVJBws+D65L17CznuvG2OTsEaJKqapKyZ9zyXn3/WBrXEwMqbeNJn7oEFluRFRbEuZCkDAnROXZMMvDN3eUUHRAwWCGC5+1cvZtco5XUfn8hYVkv31Ya1ynjqQ9cD+GtNphrkyIqiVhLgQJc0JUrtIihUkPO1g6IdhK16yPgWvH2EiuL+OWROUomb+A7LfeIVBcjDY2lpRbRxM/dLC0xokaQcJcCBLmhKga//ziYfzddhx5KuY4DSNetdLjarO84YqTFrDbyX7vA0pmzwGCY+PSHnoQQ+3UMFcmxOkjYS4ECXNCVB17rsKEe0r45xcPAG0HG7n6fRu2VG2YKxPRxrFkKVlvvEkgvyA4U/WWUSRccD4arfwuiZpFwlwIEuaEqFqqqrLsOzcTH3LgKlax1tJw1bs2OpxvCndpIgooLhc5H44pP4uDuXUr0h95CGPdumGuTIjwkDAXgoQ5IU6Pgn0Bvr7Nzua5PgB6XmPm0lesxNikZUWEVrpuHVkvv4YvMxONwUDyDdeRNOJSWTdO1GgS5kKQMCfE6aMoKnPHuPjxGQc+NyQ30HLdxzaa9jKGuzQRQVSfj7yvvqbgu0mgKJiaNCb9sUcxNW4U7tKECDsJcyFImBPi9Mvc5OeLm+3s/cePRgPn3hvL+U/K6cAEeHbuIvPlV/Bs2w4aDUmXjST5+mvRGiXwCwES5kKSMCdEePi9KtNfdTLzf6WoCtRtp+eGsTbqtJSFhmsiVVEo+vEncj/+FNXnw5CeRtqjjxDbtk24SxMiokiYC0HCnBDhtX2Jjy9vKSZvp4LeBBc/b+Ws0bLQcE3iz88n87X/Ubp8BQDxgweResdtaGNjw1yZEJFHwlwIEuaECD93icL3jzhY9E1woeGW/YMLDSeky0D36q5kwV9kv/EWAbsdnc1G7QfuI65P73CXJUTEkjAXgoQ5ISLHP794GHenHWeBiiVRw1XvxXHmheZwlyWqgOJykfPBRxRPnwFAbJfOpD/8IPrk5DBXJkRkkzAXgoQ5ISJLcVaAr28vYcMfXiC4hMmI16yYrbKESXXh3ryZAy++jG/ffjQGAymjbybhogtlAWAhjoOEuRAkzAkReVRVZe7HLqY+6cDvgZTGOm4Ya6NRF0O4SxOnQA0EKJj4PXlffAmBAMbGjajzxGOYGsmSIzWVqqoUOVw4Sj0oikqNDCJlNIBWq8EaayLBGvOfpz6UMBeChDkhIteBjX4+v9HO/nV+tDoY8qiFQQ/GotPL5Iho48vNJfPlV3H9sxqAxEsuptYto2TJkRouu6AEDZBoi0Wv09boczerqoo/oFBoL0UFaifFhdzvaLlF2raFEBGnTks9j8xN5Jy7YlAC8OuLTt4aXET+7kC4SxMnoGTBX+waNRrXP6vRJSaQ8fKLpN55uwQ5gdvjo1aiFYNeV6ODHIBGo8Gg11Er0Yrb4zup55AwJ4SISAaThuEvxXH3zwnEp2nZvsTHCz0LWP69O9yliWNQ3G6y3nqbA888h1JSgqVrFxp++gnWbl3DXZqIECqgreEh7t+0Gs1JdzdLmBNCRLSWZxt5cnES7YcacdtVPr/Rzpe32HHZlXCXJkJwb9/B7tvupPiXaWgMBlLvuI2Ml19En5QY7tKEqLYkzAkhIp61lpbRE+K58p04DDGwdIKbF3sVsGPZyXVJiMqnqiqFP/zIntvvxLt7N8b69aj/wbskDr+kxnejCVHVTijMvfLE/5j3xwJKnaVVVY8QQoSk0Wjoc2MMjy9Iol57Pfm7FN4YWMiM150ogRo5jytiBIrtHHjqGXLe+wDV5yN+6BAafPQB5jPOCHdpQtQIJxTmWndoxexpf/LkXc/yxQffsHn9lqqqSwghQkprrueh2YcmR/z8f07eHlpEwT6ZHBEOpatXs+vm0TgWLUZrsVDnmadIe+A+tDEx4S5NiBrjhJcmUVWVDWs2sXT+Mtb+vR5bvI3ufbvQrU9XkmpFz5gIWZpEiOi3YbaHr24pwZ6jEFt25oiOcuaI00INBMj/Zjz548aDomBu3Yo6TzyOIa12uEsTUSCa34O3btzGuE8m8NxbTx33Y5bMX0aDJvVJz0g76n5H+7kc7T79cVdSRqPR0Lp9S1q3b4nT4WThnMXM+PF3Zv74B81aN+XsQX1p1a7liT6tEEKcsFYDTDyxxMA3t9lZ95uXT6+20/sGLyNeicMYK+O0qoovN5fMF17GtXYtaDQkXXUlta6/Fo1Ozqkraqb5sxYye9qf2IvtpGekccnVF3FG88bl9+/ZuZc9O/cy8rrhVfL6Jz0BYue2Xfw8cRp//DqH+AQbgy8eSK3UZD579yumjPuhMmsUQoj/ZEvRcvv38Yx83YreBH994eblvgXsWyuTI6qCY9Fidt08GtfateiSk6j7+quk3HSDBDlRY61c8jdTxv3AwAsG8MjzD9CoaUM+ev0TCvIKy/dp17ENa1etr7IaTqhlrqS4hGV/rWDJgmXk5eTT9szW3HTXdbRo27x8n669OvP+q2MYfvXFp1SY0+Hk27ET2bR2C5Y4CxeMHELnnp1C7rt31z6mjPuRvbv2YTIZGXjBOZx1Xt9Ten0hRPTQaDScfWssTXsZ+Ox6O1mbA7x6diGXvGDlrNH/fXoccfwUr5e8T8ZSODX4Yd3StQtpjzyEPjF6htcIURX+nDGPbn260OvsHgCMuPYSNq7ZxF+zF3LBZcMAaNryDDxuN3t27qV+o3qVXsMJhbmn7vk/UmrXonu/bnTr0xlrnPWIfdLqptGg8akXOumrqej0el764Dn27d7PmDfGklE/g/S6FfubHSUOPnztEy656kI6dG1PwO+nqKD4lF9fCBF96rY18NiCJL5/tIS/vnAz6SEHG+d4ufZDG9ZashLTyfLu28eB51/Cs3Ur6HSk3HwTiZcOR6OVn6moPLfF5YTldT8qST3px/r9fvbu2seAIWdV2N6iTXN2bt1V/rVOr6Nlu5asWbmuSsLcCf1PvPOx23ji1UcYMOSskEEOICbGzN2P33FKRXncHlYvX8Ow4YMwmU00ad6Yth1bs2zhiiP2nTNjHi3bNadLr04YDHrMMWbSMmQArhA1lTFWw1Xv2rh5nI2YBA1rZ3h5oWcBm+d7w11aVLLPnsOu0bfj2boVQ3oa9d99m6SRIyTICQE4S5woikJcfMXzqcbFW7EXl1TY1q5jG9asXFsldZxQy9z0KTMZdc8NxFoqTjl3udx8+tbn3P347ZVSVE5WLlqdltT0Q2k5o14dtm3afsS+u7btpk69dN587l1ys/No2KQ+I64bHnJm7cI5i1k4dzEAQ666FKJ0Jo0Q4tg6XmimYUcDn99oZ/sSH+8MK2LQQ7EMfcyCTi/drseiuN3kvP8BxdNnAhB3Vj9q338fOqslzJWJ6upUWsiiQav2Lfjqo3EU5BVW+uofJxTmtm3aTsDvP2K73+tj+5YdlVaUx+PFHFNxeQFzrBm323PEvkWFxezbvY87HrmVOnXT+em7X/jyw2+4/+m7j9i3V/8e9Oof7NPelVlQafUKISJTUj0d981IYPorTma8VsqM10rZPM/HjZ/bSK4vA/b/i2fnLg783wt4d+9GYzSSeuftxA8dImMPhfgXS5wFrVZLyb9a4UqKHdj+1VpXkFeIwWggzha6Z/NUHFc7+d5d+9i7ax8A+/dmln+9d9c+du/Yw8I/l5CQGF9pRZlMRtyuiifTdrs8mM2mI/Y1GPS069SWBo3rYzAaGHzxeezcugtXqavS6hFCRC+dXsP5T1q5d1oCCXW07Fjq48WeBaz6yX3sB9cwqqpSNH0Guw87JVeDD98nYdhQCXJChKDX66nXsC6b1lU8icKm9Vto1LRhhW1rV62jZdvmGIyGyq/jeHZ6/em3ym9/+NrHR9xvMBi49NpTm716uNS0FJSAQk5WLqlpKQDs33OAtLpHLraXUb9OxT8y8vdGCBFCsz5GnliYxDd32FkzPbgmXZ8bvVz6ShzGGPnDoZSWkvXWO5TMngOAbdB51L7rDjmTgxDHcPbgfnwz5lsaNKlP46aN+GvOIooLi+k9oGeF/dasWke/gX2qpIbjCnPPvvkEqgrPPfAiDz57L1bboTETOr2eOJsVbSUOhjWZTbTv3JZpU2Zy5U0j2b/nAGtXrQvZddqtT1c+e/dL+g3sQ3pGGjN//IPGzRoREyt/gIQQFVlrabn1u3jmfuxi6hMOFnzuZvsSHzd9GU+dlie8hnq14d62jQP/9wK+ffvRmM3Uvvdu4geeG+6yhIgKnbqfidNRym8//YG9yE563XRue/BmkmodGpdfWFDE/j0HaNOhVZXUcMKn8zpdnA4n4z+dyOZ1W7DExXLByKF07tmJbZt38NHrn/DG2FfK910wayG//fwHXo+PJs0aMfL64SQmH31wYTSfSkQIcer2rvHx2fV2srcGMMTAyFfj6HW9uUZ1J6qqStHPv5D74RhUnw9T48akP/0Epvr1w12aqOai+T34ZE7nNX/WQv5Z9s8xV/uostN5/bN8DW3PbI1Or+Of5WuOum+HLu2O9XTHzWK1cMt9Nx6x/YzmjSsEOYA+5/Sizzm9Ku21hRDVX712Bh6dn8ikBx0sHu9m/N0lbJrn5ap34oiJr/7LbgQcTrL+9waO+QsAiD9/KKm334bWdOTYZCHEqVm7ah3tOrWtsuc/Zpj7/L2vePG9Z4mLj+Pz97466r7vfv1GpRUmhBBVzWzVcu0YGy3ONvLtvSWsnOJh10ofo76Ip2Hnyh+kHCncmzcHu1Uzs9DGxlL7/nux9T873GUJUW3d8fDoKn3+iO1mrWrR3MQrhKh8Odv8jL3Bzt5//Gj1cNGzVgbcFYNWW326XVVVpWjqD+R8/Cn4/ZiaNqXO009gzMgId2mihonm9+D83ALWrFzL2YP6Vfpzn2w3a/XvSxBCiOOQeoaeh2Yl0v+OGBQ/TH3SwYcjirHnKuEurVIESko48PRz5HzwEfj9JFx0IfXfe1uCnBAnKDklqUqC3Kk4rjFzx6syx8wJIcTpZjBpGPFKHM37Gvn6Vjvrf/fyUs8CbvjMRvO+xnCXd9JcGzdy4P9exJ+djdZiIe2hB4jrWzVLJAghTr/jGjN3vGTMnBCiOmg3xMQTi5P4/EY72xYFTwU2+JFYhj5qQauLnm5XVVUp/H4yuZ9+BoEA5ubNSX/6CYzp6eEuTQhRiWTMnBBC/IeAX2Xay05mvl6KqkLTXgZu/NxGQp3IPxVYoNhO5quv41yyBIDE4ZeQcssoNIbqO7FDRA95Dw5NxswJIUQl0+k1XPCUlbt/TsBWW8vWhT5e6FHA2plHnic6krjWrWfX6FtxLlmC1mqlzvPPkXrHbRLkhKimInadOSGEiBQtzjLyxKIkvrzFzsbZXj4cUcw5d8Vw4bNW9MbI6XZVFYWCSd+TN/ZzUBTMLVtQ56knMaTVDndpQogqJOvMCSHEcbClarlzajx/vF3Kz//nZNZ7LrYt9nHTF/HUahj+bld/cTFZL7+Kc9lyABJHDCdl1E3SGidEDSBj5oQQ4gRtX+Lj8xuLKdirEBOv4er34+h4kTls9ZSuWUvmCy/hz8tDa4sj/eGHsPbsEbZ6hDgWeQ8OTcbMCSHEadKku4HH/0qi/TAjrmKVT6+xM+G+Enzu0/vZWFUU8sd/y977H8Sfl4e5dSsafjxGgpwQNcwxu1n/be+uffw5cz5ZB7IASKtTm7MH9aNew7qVXpwQQkQqS5KW0d/GM+8TF1MedzB/rIsdS33c9KWNtGYn/Kf1hPkLC8l8+VVKV6wEIOnykdS68QY0+qp/bSFEZDmhlrnlC1fy+tNvYS+y07p9S1q3b0lJcQn/e+Ztli9cUVU1CiFERNJoNJw1OpaHZieS0kTHvrV+XulbyNIJrip93dK//2HXzbdSumIlOpuNjJdeIOWWmyXICXEabN24jWfue/6EHrNk/jIy92dVUUUn2DL36+TpDL10MOddcE6F7b//PItfJ8+gS6/OlVqcEEJEg/odDDw2P5Fv7y1hxfcevrylhM3zfFz2RhwmS+XNdlUDAfLHfUv+N+NAUYhp25b0Jx/HkFKr0l5DCHFitm3azuzpc9m7ay/FhXauuvlyuvftWmGfPTv3smfnXkZeN7xKajihljmH3UnHbu2P2H5mt/aU2B2VVpQQQkSbGJuWGz+zcfX7cRhiYPF4N6/0K2D/en+lPL8/P599Dz9K/ldfg6qSfPVV1HvzdQlyQoSZx+0hvW4aw6++GIMx9Ozxdh3bsHbV+iqr4YTCXNNWZ7B14/Yjtm/duJ0zWjSptKKEECIaaTQael0XwyNzk0hvoSNrc4BXzyrgry9cnMrCAc4VK9l1y62U/v0PusQE6r76MrVuvB6NLvxLoghR07Xu0IoLRg7lzK7t0WhCt8Q3bXkGHrebPTv3VkkNx7Vo8EGt2rXgl0nT2LNzLw2bNABg1/bdrF6+lsGXnFclBQohRLTJaKXnkblJTHqohEXfuBl/dwmb53u58p04YmzH/xlaDQTI+/IrCr79DlSV2A4dSH/iUfTJyVVYvRDhsbn/uWF53eZz/qjy19DpdbRs15I1K9dRv1G9Sn/+41o0+N8W/bmERX8uqbBt8tdT6XtOr8qrTAghopjJouGaD20062tkwr0lrJjsYfcqP6O+slG/w7EX8vXl5JD54su41q4DrZbk668l+corpDVOiCjVrmMbfvv5D4ZdOrjSn/uYYU7O6iCEECev2+VmGnbSM/Y6O/vW+nl9QCGXvGjlrNEx/9kl41i8hMxXX0Oxl6BPTib9yceIbX/keGUhqpPT0UIWTq3at+Crj8ZRkFdIUq3ESn1uWTRYCCGqWO2meh6ek0jfUTH4vTDpIQcfX1mMs1CpsJ/q85Hz4Rj2P/EUir0ES9cuNPh0jAQ5IaqBgrxCDEYDcTZrpT/3CS9KVOosZcPqjRTkFxHwV5ylNfhiGTcnhBChGMwarngrjub9DHxzRwmrf/Wyd00BN30RT+OuBrz7D5D5wou4N28BnY5aN91A0sgRaLTymVuI6mDtqnW0bNv8P2e8nooTCnM7t+1izBtj0ev1OEocJCTGYy+yo9frSUpJkjAnhBDH0PEiM/XaG/jshmJ2r/TzxnmFXHHT39Te9CFKaSn62rWp89TjxLRqFe5ShRDHweP2kJudB4CqqhTmF7Jv935iLbEVulPXrFpHv4F9qqSGEwpzP074hS49OjL8mot56JbHueux2zGajHz5wTf06NetSgoUQojqJqWRjgd/T+SnJ/NRp39Cyqo/UQBz117UfeIBdHFx4S5RCHGc9uzcy7svfVj+9fSpvzF96m907d2Fa0ZfAUBhQRH79xygTYeq+ZB2QmHuwN5Mrhp1GRqNBq1Wg9/np1ZqMhdePoyvPhxHl16dqqRIIYSobgL7d9Mx70W8qbvwKwbm77mK3fsGcuMwE816h7s6IcTxatryDN775s2j7rN21XqaNGuExWqpkhpOaDCGXn9oSnycLY6C/AIATCYTxYX2yq1MCCGqIVVVKfp1OrtvuwPvrl0Y69Uj5aV3cDYZQnGmyttDi5j2ihMlcPKLDAshIsvaVeto16ltlT3/CbXM1W1Yl9079pKankrTlk34dfIMSoodLF+4kjr106uqRiGEqBYCDifZb75Fydx5ANgGnUftu+5AGxPDfTNUfn3JyW//K+XXF51sWeDlxs9sxKfJunJCRLs7Hh5dpc9/Qi1z5186hPhEGwBDLx2CNc7K5K+nUlpayhU3jqiSAoUQojpwbdzI7tG3UjJ3HpqYGNIff5T0hx9EGxMDgE6v4cKnrdz1YwJxKRq2zPfxQo8CNszyhLlyIcThkmolcdZ5fcNdRgUa9VROGBjFdmUW0DA9KdxlCCGqOVVRKPhuEnlffAmBAKamTanz9BMYMzL+8zHF2QG+GGVn81wfAAPvi+WCpyzoDKEXGRYi2sh7cGhH+7kc7b4TXmcOIDc7j+wD2QCkZaRRK1XOEyiEEP/mz88n85XXKF25CoDEEcNJGXUTGsPR15mKr63j7h8T+O3NUn55wcnvb5WybaGXG7+IJ7m+dLsKISo6oTDnLHEyfuxE1v29vvw0NKqq0ubMVlw16nIscVUzS0MIIaKNY+kysl59nUBREbqEBNIeeQhrt67H/XitTsPghyw07WXg8xvt7Fjm58VeBVz7gY0OF5iqsHIhRLQ5oW7WT9/+nNysPC6/cQQNmtQHYPf2PUz8cjK1atfi5ntuqLJCK5s08QohqoLi9ZI39nMKJ08BILZTR9IffRh98sn3YDjyFb6+1c7amV4A+t0Sw/AXrRjM0u0qopO8B4d2st2sJzQBYuPazVxx00gaN2uETqdDp9PRuFkjLr9hBJvWbD7xqoUQohrx7tnLnjvvCQY5nY5ao26i7qsvn1KQA7Ama7ltUjyXvmJFZ4B5n7h4bUAh2Vv9x36wEKLaO6EwZ42zYjQZj9huMBqwxMVWWlFCCBFNVFWlaPoMdt16O55t2zCkp1P/3bdJvvLySju3qkajYcAdsTw0K5GUxjr2rfHzcp9Clk5wVcrzCyGi1wn9lRl80blMGfcjRQVF5duKCor4YcLPDLpoYGXXJoQQES9QUkLm8y+S/b83Ud1ubOcMoMEnHxHTskWVvF6DjgYeW5BI50tNeJwqX95Swpe32HE7lCp5PSFE5DvmBIiXHnutfLIDQH5uAc/c/wIJifEAFBUWYzAYcNgd9Dyre9VVKoQQEaZ0zVoyX3oZf04umpgYat9zF/EDz63y142xabnxcxvN+7mZ9HAJSye42bncx6ivbNRrd/SZskKI6ueYYa5Dl/anow4hhIgaqt9P/tfjyP92AigK5ubNSX/ysaOuHVfZNBoNva+PoXE3A2OvKyZzY4DXzi5k+EtW+t0SU+FDuBCiepNFg4UQ4gR4D2SS+dLLuDdsBI2GpCsvp9Z116LRn9SynZVTk0tl8qMlLPjcDUD7oUau+dCGJalyxusJUdnkPTi00zKb9aDN67cy748FzP/jL7Zu3HYyTyGEEFFFVVWKf/+D3bfcinvDRvQpKdR743VSbroxrEEOwBij4cp3bIz62kZMvIbV07y82KuArQu9Ya1LiOpo68ZtPHPf8yf0mCXzl5G5P6uKKjrBRYOLCor49J0v2LtzH/FlY+aKC4up36geN997Q/k2IYSoTgIlJWS//S4lf84FwNq3D2n334vOZgtvYf/S6WIzDc408NkNxexa4eetIUUMfczC4Idi0eqk21WIqvD7z7NYvWItOZk56A16GjZpwPkjh1KnXnr5Pnt27mXPzr2MvG54ldRwQi1zk7/5Aa1Wy9P/e5zn33ma5995mqf/9zharZbJ3/xQJQUKIUQ4la5eza6bR1Py51w0ZjNpDz1AnWeeirggd1Cthjoe/D2R8+6PBRV+fdHJ28OKKNwfCHdpQlRLWzdtp885vbjv6bu567Hb0Oq0vP/qGJwOZ/k+7Tq2Ye2q9VVWwwmFuc3rtjDyuuEVzsVaKzWZS6+5mE3rtlRqYU6Hk0/f/pwHbnqUp+99nhWLVh51f7/fzwuPvMJTdz9XqXUIIWom1ecj99PP2Hv/Q/hzcjG3bEHDT8cQP3hQxE8u0Bk0XPSclbt+TMCWqmXrXz5e7FnAmumecJcmRLVzx8Oj6d63K3XqpVOnXh2uvfUqHHYHO7bsKt+nacsz8Ljd7Nm5t0pqqJyBHlXwd23SV1PR6fW89MFz7Nu9nzFvjCWjfgbpddNC7j972p9Y46x4XPLHSghxajx79pD54it4tm4FrZbkq68k+Zqrwz427kS17G/kicVJfDXazoZZXj66rJizbo3hkuflVGAi8nQaPzEsr7vyqssq9fncbg+qqhJriSnfptPraNmuJWtWrqN+o3qV+npwgi1zzVo3ZfI3P1CYX1i+rSCvkCnjfqRZ66aVVpTH7WH18jUMGz4Ik9lEk+aNaduxNcsWrgi5f15OPssXreTc8wdUWg1CiJpHVVUKf/yJ3aNvx7N1K4a0NOq9/Qa1brg+6oLcQbZULXdMieeSF4KnAps7xsWrZxeSuUlOBSZEVZjyzQ/UbZBBo6YNK2xv17ENa1aurZLXPKG/TpdeczGfvPU5zz7wIvEJZRMgioqpUzedS6+5uNKKysnKRavTkpqeWr4to14dtm3aHnL/yd/8wPkjhmA0Hn2xzIVzFrNw7mIAhlx1Kci0aCFEGX9BAVmv/Q/nsuUA2M47l9Q770BnsYS5slOn1Wo4955YmvUx8NmNdvav8/NKvwJGvhZHz2vNEd9tLGqGym4hC4ep439i+5ad3PfUXWj/dSq/Vu1b8NVH4yjIKySpVmKlvu4JhTmL1cKDz97L1o3byM7MAaB2ndq0aNOsUovyeLyYY8wVtpljzbjdR3ahrl6xBkVRaN+53TGXSenVvwe9+vcAguu1CCEEQMlfC8l+4y0CxcVobXGk3Xcvcf36hrusStego4HHFyTy3QMOlk5wM+7OEjbM9nLVu3HEJsiadEKciinjfmTVkr+56/HbK8wtOKggrxCD0UCczVrpr33cYU5RFB4a/TiPvvggLdo2p0Xb5pVezEEmkxG3y11hm9vlwWw2VdjmcXv46btfufXBm6usFiFE9RVwOsl5/0Psv/0OQGynjqQ9/BCGlFphrqzqmOO0XP+JjZb9jUy4r4RVP3jYtcLHjZ/ZaNLDGO7yhIhKk7/5gVVL/+Hux24nrU7tkPusXbWOlm2bYzhGL+LJOO4wp9VqSUpOJOCv+untqWkpKAGFnKxcUtNSANi/5wBp/5r8kJudR35eAW+/8D4AAb8fV6mbx+98hgeeuYfkFOlGFUKEVrp6NZmvvI4/OxuN0UjKLaNIuOhCNNqa0ULV7XIzjbro+fwmO7tX+nljUBFDHg2uSafTS7erEMdr0pdTWL5wBaPuvZFYSwz2IjsAJrMJ02GNUGtWraPfwD5VUsMJdbMOumggP0/8lWtvuwprXOU3Ex5kMpto37kt06bM5MqbRrJ/zwHWrlrH/U/fXWG/9LppPP/20+Vf79i6i++/nsojz9+PtQqaMYUQ0U/xesn7/EsKv58MqoqpWVPSH3sEU4MG4S7ttEttoufB3xP55QUnf7xdyrSXnGz608sNY20k19eFuzwhosKC2QsBeP+VjypsH3zxQIZcMgiAwoIi9u85QJsOraqkhhMKc7On/0l+bgFP3f0cCUkJGE0Vm+Qfe+mhSits5PXDGf/pRB6/4xkscbFcdv1w0uumsW3zDj56/RPeGPsKOp0OW8KhhTst1li0Gk2FbUIIcZB72zYyX34V785dUb3kSGXSGzVc/H9WWvY38uXNdrYv9vFirwKueieOTpeYj/0EQtRw733z5jH3WbtqPU2aNcJirZoJVSf0F6xDl/ZoNKCqVVJLBRarhVvuu/GI7Wc0b8wbY18J+ZimLc/g+XefqerShBBRRg0EKJjwHXlffQOBAIaMOqQ/9ggxrarmU3I0anGWkScXJ/H17XbWzvAy9jo76//wMvJ1K2Zrzeh6FqKqrF21jnad2lbZ8x9XmPN6vPw44RfWrFpLwK/QrHVTRlx7cZV2tQohRGXw7tlL5quv4d64CYCECy8g5ZZRaGNijvHImsdaS8ttE+OZ/5mLKY85WDzOzbbFwckRDTtV/qBtIWqKOx4eXaXPf1wft6ZPncnSBctp3b4VnXqcyZb1W5j4xZQqLUwIIU6FqigUTv2BXaNvw71xE/qUFOq+/iq177lLgtxRaDQa+o2K5dH5SWS00ZO7PcDr5xQy839OlMBp6JYRIsIl1UrirPMia+mi42qZW71iLVeOuoxOPc4EoHPPjrz1/HsoinLEonhCCBFu3sxMsl5/A9c/qwGwDTyX1DtvR2eV3oTjVaelnkf+TOTHZx3M+cDFT8852TDby/Wf2EiqJ5MjRM2VnJLE2YP6hbuMCo4riRXmF9GkeaPyrxs2aYBOq6W4sLjKChNCiBOlqipFv/zKrlGjcf2zGl1CAnWee4b0Rx+WIHcSDGYNI16J486p8dhStWz9y8cLPQpYMcV97AcLIU6b4wpziqKg+9dsL61ORyCgVElRQghxonw5Oex7+FGy33oH1eUirl9fGn7+KXF9eoe7tKjX+lwTTy5Jou1gI65ilc+ut/PlLXZcdnkPECISHPds1q/HjEd/WKDz+XxM+HwSRuOh5UlG339T5VYnhBDHoKoq9pm/kfPhRyjOUrS2OGrfcze2s88Kd2nVSlxKcHLEgs/dTH6shKUT3Gxb7OWGT+TMEeLEaQj+35XzAh+iqion+9M4rjDXtXfnI7Z16dnpJF9SCCEqhy8nh+w338a5bDkA1p49qH3/veiT5OwvVUGj0dD3phia9Tbw+U129q4Onjli0IOxDH3Ugs4gb8zi+Oh1Wjy+AGZjzV3j8d88vgB63cnNQ9Co6ulYNS7y7MosoGG6/MEXIhqpqkrxjJnkfjQm2BoXF0ftu+4gbkB/+aR/mvi9Kr++6OT3t0pRVWjQSc8Nn9qo3VTenMWxOV0eCuwuUhKtmAy6Gv3/VlVVPL4AuYUOkmwxWGJMIfc7Wm6RMCeEiCq+7Byy3zqsNa5XT2rfd4+0xoXJ1r+8fHmLnYK9CsZYuPTlOHrfYK7Rb87i+DhdHopKXPgDCjUyiJTREGypTIj77yAHEuZCkjAnRHRRVZXiX6eR+/GnKKXSGhdJXMUK3z1YwrLvPAC0Oc/I1R/EEV9bljARorJImAtBwpwQ0cN7IJPsN96k9O9/ALD27kXte++W1rgIs3yymwn3leAqUrEma7jqPRsdzv/vlgYhxPGTMBeChDkhIp+qKBT9+BO5Yz9HdbvRJSSQevedxPXrK61xEapwf4Cvb7Oz6U8fAD2uNjPiVSsxNllgXohTIWEuBAlzQkQ2z+7dZP3vTdzrNwAQ1/9sUu+8HX1CQngLE8ekKCrzPnHxw1MOfG5IbqDlujE2mvaWJUyEOFkS5kKQMCdEZFJ9Pgq+m0T+uPGoPh+65CRq33sPcb16hrs0cYIyN/n54mY7e//xo9HAgLtiueApCwaztKoKcaIkzIUgYU6IyOPatJms/72Bd8dOAOKHDiFl9M1yKq4o5veqzHjNycz/laIEIL2ljus/sVG/gyHcpQkRVSTMhSBhTojIobhc5H35FYVTfgBFwZCeTu0H7sPS8cxwlyYqyc7lPr68xU7OtgBaPQx9zMJ598ei00srnRDHQ8JcCBLmhIgMzuXLyX7rXXxZWaDVkjj8YmrdcD1aszncpYlK5i1V+eEZB3PHuIDgQsPXjbGR3kIWGhbiWCTMhSBhTojw8hcVkfvhGOyzZgNgOqMJaQ/ch7l58zBXJqraxj+9fHO7ncJ9CnoTXPi0lf53xKDVSSudEP9FwlwIEuaECA9VVbH/MYvcD8cQsNvRGI0kX3cNSSMuRaOXFpqawlWs8P2jDhaPcwPQpLuBa8fEkdpEfgeECEXCXAgS5oQ4/bz79pH99ruUrvobgNgzO1D7/nsxZmSEuTIRLmtneBh3Vwn27ODpwC7+Pyt9b45Bq5VWOiEOJ2EuBAlzQpw+Ryw3YrORcust2M4bKIv/Chz5ChMfKmHF98HTgTXrY+DqD2ykNJLTgQlxkIS5ECTMCXF6lK5dR/abb+PdvRsA23nnknLraPTx8WGuTESav38Kng6sJFfFZNFw0f9Z6DtKWumEAAlzIUmYE6Jq+YuLyftkLMUzZgJgqJtB7XvvkeVGxFE58spa6SZLK50Qh5MwF4KEOSGqhqqq2Gf+Ru7HnxKw20GvJ/mKy0i66kq0Rjmdkzg+h7fSGWPhwmetnDVaWulEzSVhLgQJc0JUPs+u3WS/9Q6utWsBiO3QgdR778JUv36YKxPRyJGnMPHhQ2PpmnQ3cPUHcaQ1kxmvouaRMBeChDkhKo/icpH/zTgKvp8CgQC6xARSb7uVuAH9ZYKDOGWrf/Uw4b4SirOC69Kd/4SFAXfJ2SNEzSJhLgQJc0KcOlVVcSz4i5wPPsKfmwsaDfHDhpAy6iZ0cXHhLk9UI85ChSmPOVg8PrguXYOOeq7+wEbdNtJKJ2oGCXMhSJgT4tR49x8g5733cS5bDoCpWVNq33sPMS3kDA6i6qz/w8P4u0so3Keg1cPA+2IZ8rAFg1la6UT1JmEuBAlzQpwcxe2mYMJ3FHw3CdXnQ2uxUGvUjSQMG4pGJzMORdVzlyj8+KyTeZ8Ez/Ga1kzH1e/H0aSHTLAR1ZeEuRAkzAlxYlRVxbFwUbBLNTsbANvAc0m55Wb0SYlhrk7URNsWeRl3ZwnZWwNoNND35hgufMZCjE0b7tKEqHQS5kKQMCfE8fPu20fO+x8e6lJt0pjUe+4itk2bMFcmajqfW2X6a05+f6sUxQ8JdbRc/kYc7YeZwl2aEJVKwlwIEuaEOLbgLNXxFEyeAn5/sEv1xutJuOB86VIVEWXfOj/j77Kza4UfgA4XmLjsdSsJdeT3VFQPEuZCkDAnxH9TVZWS2XPI/fhT/Pn5ANgGnUfKzTehT5QuVRGZlIDKvE9d/PScE49DxWzTcNGzFvrcJIsNi+gnYS4ECXNChObeuo2c997HtW49AOYWLUi96w5iWrYIc2VCHJ+CfQG+u7+EtTO8ADTsrOfKd+Ko184Q5sqEOHkS5kKQMCdERf7CQvI+/5Li6TNAVdElJpAy6iZs5w1Eo5UB5SK6qKrK3z97mPSQg+JMBa0Ozr4thmFPWDBb5fdZRB8JcyFImBMiSPX5KPzhR/K/GYfiLAWdjsSLLyL52mvQWS3hLk+IU+KyK/zygpO5H7tQFUjM0DLy9Tg6nC8TJER0kTAXgoQ5UdOpqopzyVJyPhqDb99+ACxdu5By+61yLlVR7ez+28eEe0vYvSo4QaLtICMjX4+jVkOZICGig4S5ECTMiZrMs3MnOR+OoXTlKgCM9eqRcttorN27hbkyIarOwQkSPz/vxG1XMZhh0EMWzr0nFoNJJkiIyCZhLgQJc6Im8hcWkvflVxRPmwGKgtZiIfnaq0m86EI0BhkcLmqG4uwAU59wsGyiB4DUJjoue8NKqwHS9Soil4S5ECTMiZpE8XopnPoDBeO/DY6L02pJOH8Yyddfiz4+PtzlCREWm+d7+e7+ErI2BwA480ITw1+yklxful5F5InKMOd0OPl27EQ2rd2CJc7CBSOH0LlnpyP2mzVtDssWrKAgvxCL1UKfc3pyztD+x3x+CXOiJlBVlZI/55I39nN8WVlA2bi4W0djatggzNUJEX5+r8rs90uZ8VopHqeKIQbOuz/Y9WqMka5XETmOllv0p7mW4zbpq6no9Hpe+uA59u3ez5g3xpJRP4P0umkVd1ThmluvpE69dPJy8vng1Y9JTEqkU48zw1O4EBGidN06cj/6GPfGTQAYGzYk9dZbsHTtEubKhIgceqOG8+630PUyM1OfdLBisodfX3SyeLyLES/H0W6oEY1GQp2IbBHZMudxe3jk1id5/OWHSE1PBeDrMeOJT4znwsuGHfWxk7+eigqMuPaSo+4nLXOiuvLu30/up5/hmL8AAF1iIrVuuJ74wefJKbiEOIYtC7xMfLCEAxuCXa8tBxgZ8YqV9BYR2/Yhaoij5ZaIXDkxJysXrU5bHuQAMurVIWtf1lEfp6oq27fsJD0j7aj7CVEd+QsLyX73fXZefxOO+QvQmEwkX3MVjb/5koRhQyTICXEcmvUx8vjCJEa8ZiUmXsPG2V5e6F7ApEdKcBYq4S5PiJAi8qOGx+PFHGOusM0ca8bt9hz1cdOn/oaiKHTr2zXk/QvnLGbh3MUADLnqUpCWOVENKC4XhVOmUvDdJJTSUtBosJ03kFo3Xo8hJSXc5QkRdXR6Df1vi6XLpWZ+edHJX1+4+PNDF8u+c3P+UxZ6Xx+DTi9dryJyRGSYM5mMuF3uCtvcLg9m839PG5/3xwKW/bWCe5+6E4Mh9LfVq38PevXvAQSbK4WIZmogQPGM38j/6mv8+fkAWLp1JeXmUZgaNwpzdUJEv7gULVe+HUffm8xMetjB1r98fHefg/mfuhj+kixlIiJHRIa51LQUlIBCTlYuqWnBloX9ew6Q9u/JD2UWz1vKrF/mcM+Td5KYlHAaKxXi9FNVFceCv8j77Au8e/cCYG7ejJRbbib2zA7hLU6IaqhuWwP3TU/g7589TH3CwYENAd67qJjWA40Mf1HG04nwi8gxcyazifad2zJtykw8bg87tuxk7ap1dO3V+Yh9ly9cyS/fT+eOR26lVmpyGKoV4vQp/fsf9txxNwee/T+8e/diyKhD+lNPUP+D9yTICVGFNBoNHS8088yKZC7+Pwtmm4b1vwfH0317bwn2XBlPJ8InImezQnCdufGfTmTzui1Y4mK5YORQOvfsxLbNO/jo9U94Y+wrADxz3wsUFRah1x/6ZNSlVycuv2HEUZ9fZrOKaOLesoXcz76gdPkKIDhDNfnaa0gYOhiNXloFhDjd7LkK014KjqdTAmCyahh4XywD7ojFZJHxdKLyReWiwVVNwpyIBp7du8n74qvyZUa0lliSLhtJ4vBL0MbEhLk6IUTmJj9Tn3Cw7ncvAPFpWoY9YaHH1WaZJCEqlYS5ECTMiUjmy8oi76uvsf8xGxQFjdFI4iUXkXTZZejibeEuTwjxL5vne5n6pIM9f/sBSG+h48JnrbQbIosOi8ohYS4ECXMiEvlyc8kf9y3F02dAIAA6HQlDh5B89ZXoa9UKd3lCiKNQFJWVUz389JyD/F3BMXSNuxm46FkLTXsbw1ydiHYS5kKQMCciib+ggPxvJ1D8yzRUnw+0WmwD+pN83bUY66SHuzwhxAnweVTmj3Ux83UnjvzgW2zrgUYufMZCvXaGMFcnopWEuRAkzIlI4C8spGDi9xT99DOqJ7godtxZ/Ui+7hpMDRqEuTohxKlwlyjMet/FrHdL8TiCb7WdLzUx9DELac1k4pI4MRLmQpAwJ8LJX1RE4aTvKfzxZ1R3cIFsa69eJF9/LeYmjcNcnRCiMpXkKvz2hpN5n7rwe0Gjha6XmRn6aCwpjSXUieMjYS4ECXMiHEKFOEv37tS67mrMzZuHuTohRFUq2BtgxutOFn3jRvGDVg89rjIz5BELSfXk3Mni6CTMhSBhTpxO/oJCCr6fHOxOLQ9x3Ui+9hpiWkiIE6Imyd0ZYPorTpZ+50ZVQGeAnteYOe8BC8n1JdSJ0CTMhSBhTpwO/rw8Cr6bRNG06eVj4ixdu5B83bXEtGwR5uqEEOGUtcXPtFecrJzsQVUPtdQNetBCrYYS6kRFEuZCkDAnqpIvK5uCiRMpnj4zODsVsPbsQfI1V0l3qhCigqzNfma87mT59x5UJRjqul9h5rz7Y0k9Q8bUiSAJcyFImBNVwbNnDwUTvsM+a05wnTiNBmuf3sEQ16RJuMsTQkSw7K1+ZrxeyvJJbpRAcKJEx4tNDHoglrptZUmTmk7CXAgS5kRlcm/ZSv63E3As+Itgf4mWuP5nk3zF5ZgaNQx3eUKIKJKz3c8fb5eyeLybQLBhnzYDjQx6MJYmPWTx4ZpKwlwIEubEqVJVldK//6Hgu4mUrlgJgMZgwDZoIEmXXSaL/QohTknh/gCz3ivlry9ceEuD25p0N3DuPbG0HWJEq5XThNUkEuZCkDAnTpYaCOBY8BcFEyfh3rwFAI3ZTML5w0gaMVxOuyWEqFQluQpzPipl3qcuXEXBt+zaTXWcc1cs3a4wYzBLqKsJJMyFIGFOnCjF7cb++x8UfD8Z3/4DAOgSEki8+CISLjwfnc0W5gqFENWZu0Rh4ddu5nxQSsHe4Llfbala+t0SQ58bY4hL0Ya5QlGVJMyFIGFOHC9/YSFFP/1M0Y8/E7DbATCkp5M4cgTxgwaiNZnCXKEQoiYJ+FRW/uDhj3dK2bfGD4DeFDyrRP/bY8loLTNgqyMJcyFImBPH4tm9m8LJU7H//kf58iLm5s1JHHkpcX37oNHJOlBCiPBRVZXN83zM+bCUdTO9HHw3b36WgbNvjaXtICNanXTBVhcS5kKQMCdCURWF0hUrKZgyldLlK4IbNRqsPbqTOHIEMW3boNHIH0chRGTJ2ebnzzEuFo9z43EG39aT6mvpe1MMva6NwVpLumCjnYS5ECTMicMdHA9XOPUHvHv2AqAxmbANPIek4cMx1q8X5gqFEOLYSosUFn3jZv5YF7k7AkCwC7bTJWb63RxDw856+UAapSTMhSBhTgB4D2RS9NPPFM+YieJwAKCvVYuEiy8kYegQmdQghIhKiqKyYZaXeZ+4WP/7oS7YjDZ6et9gputIM7EJ0loXTSTMhSBhruZSFYXSlaso/OEnnEuXcvCvnLlVSxKHX0Jcn95o9DKAWAhRPeTuDLDgMxeLx7lw5Af/3hlioPMlZnpdZ6Zxd4O01kUBCXMhSJireQJ2O8W//UHRL7/i27cPCC7yG9f/LBIvulDOmSqEqNZ8HpU10zws+MLF5rm+8u2pZ+jocZWZ7leaSagjE7silYS5ECTM1QyqquLeuJGiX36l5M95qF4vAPqUFBIuPJ/4IYPRJySEt0ghhDjNcrb7WfS1myXfuinOCq5Zp9FCy/5Gul9lpv0QE8ZYaa2LJBLmQpAwV70FHA7ss+dQPG06nm3bgxs1GixdOhN//jCs3bvJ0iJCiBov4FfZONvL4vFu1kzz4A9+3sVk1XDmBSa6XmameT+DLHESASTMhSBhrvpRVRXXuvUUT5tOybz5qB4PALr4eOIHn0f8sGFyvlQhhPgPjnyFFZPdLP3Oza4V/vLtttpaOl9qotMlZhp1kdmw4SJhLgQJc9WHPy+P4j9mYf/t9/JlRQBiz+xA/NAhWHv3Qms0hrFCIYSILjnb/Cz73sOyiW5ytwfKtyfV09LxomCwa9BJgt3pJGEuBAlz0U3xenEuWkzxb7/jXL4ClOCYD11SEvGDBhI/eDDGjDphrlIIIaKbqqrsWuFn5RQ3q370ULhfKb8vqb6W9sNMdDjfRJPuBnR6CXZVScJcCBLmoo+qqrjXb6D4j1mUzJuHYi8J3qHXY+3RnfjzBmLp2kWWFRFCiCqgKCo7l/lZ+YObVT94KM48FOysyRraDg4GuxZnGWXyRBWQMBeChLno4d27D/usWdhnzcGXmVm+3XRGE+IHnUfcgP7o4+PDWKEQQtQsiqKya7mff37x8M+vngpdsQYzNOtrpM15RtqeZyK5QfRNNlNVFXcggN3rxe7xUuL1Yvf6KPF5KfH6cPh8OLzesmsfyTFmHunSqUprkjAXgoS5yObLzqFk7jzsf/6JZ8vW8u365GTizhmA7ZwBmJs0DmOFQgghIBh8MjcF+OcXD2ume9i90l/h/vQWOloNMNKiv5GmvYyYLKe31U5VVRw+H4VuD0Wewy9eijweisuu7V4vxR4vxV4Pdo8Xr6Ic+8nLNIiLY+oFQ6rwu5AwF5KEucjjz8+nZMFflMz5E9e69eXbtbGxWPv0xnbuAGLbt5clRYQQIoLZcxTW/+5h3W9eNszx4rYfihl6IzTuZqDF2Uaa9zPS4Ew9OsOJhbuD4azA7aHA7abQ7aHQ46bA7aHQ7abQ4ynb5ikPcIGTiDpGrZY4oxGb0YjNaCi/HWc0YDUasRoMWI0G4gxGkswmOtVOPeHXOBES5kKQMBcZfNk5lCz4C8f8+bjWbyg/tZbGZMLavRtxZ5+FpVtXtCZTeAsVQghxwgI+le1LfWyc42XTn152r/RzeOowWTQ07m6gcW8dqd0VYpsGsAeCIe1gWDsY0vLLgluBx4P/BFrNAGL1ehLNJhJMBy9GEkwm4stux5tMxBuNxJfdthmNmHW6iJqte7TcIiPFxWmlqire3XtwLFqEY8FC3Js3l9+nMRiwdOlM3NlnYe3RHW1sbPgKFUIIcdIUVaXY46XQ46akqYfYeh4aXeYmtsjNrl2l7M9xkV/qoVTvYVGcj4DFD7sJXo5DrF5PktlMktlE4sFrk6k8sCWazRW+NlXzHh0Jc6LKqYEArnXrcCxajGPhYnwHDpTfpzGbsXTtQly/vli7dZUAJ4QQEcgTCATHmbkrjjcLdmW6KfJ4KfR4KCrr5izyeFGO1vGXVHYpo1E1GEsN6AqM6EsMGEqMGOxG9A4DiUYTddNiaNwolmatLLQ500JCkqHKv+doImFOVAl/QSHO5ctxLluOc/kKFIej/D5dfDyW7t2w9uyBpUtntGZzGCsVouqoqopfVVHKLn5FQVVBQUVVVZSy28F/wTe+f7/9acqvNaABLRq0GtBoNGjRoNGAVqNBq9GgO+w6krqHRGRQVZVSv58Sr48S76EZmnavh+KyWZuHJgF4KS4LZcUeD+5A4Ngv8C82o5HEshazJLOZRHPwOslkIsFsItlsLt9uMxrRajQ4CxR2LvexfamPHdt97Frhw1sK+QQvy/ECXlIa68hooyejjZ66rfVktNGR3FCHVlszf+9lzJyoFKrPh2vDBkpX/o1z+XLcm7dUuN9QNwNrz55Ye/UkplVLmcQgTgtFVXH7/bj8AVx+f4WLOxDA7Q/gDviD134/nkAAdyCAJxDAW36t4FUCeAIKvrJtPkXBryj4Drv4lWBYC5SFNr+iHBHMTqeDoU6n0aDTatFrNOi0GvQaLXqtFr1WU3atxVB2v16nQ6/RYNBqMei0GLW68vsNWi1GnRaDVld2X/CxxrL9DBWutRh1uvLHHO1+vVYbxp9S5AsoSvnvpcvvp9QX/P0tLfs9LvX5cfp9OH1+nD4fpT4fDp8/uGTG4ctnlC2hcTITAQD0Wm35OLPgWLPg7YNdmYkmc1mXppFEs5kEkwlDJRzbgF8la3OA3at87P7bz56/fexb4y8/h+zhjLFQu6metOY60prpqd1MR+2memo11GK2Rv/vmUyACEHC3KlRFQXPzp2Urvqb0pWrKF2zFtXtLr9fYzQS26E9lq5dsXTrgjEjI4zVimjjCwQqvAGVlN12ln3t9B1683L4fJT6/WXbfLj8fpyHveGF2+GtZTqNBq32YItaWQtbWdub5vDbZY9Vy6/LWu3KWvAUNdjKoqCiKCoKweAaOKwVMFpoD4bHw8Ni+W1teeD79+0KgVRTMZweHmC1h/3cdRpt+fE4vHUz2JijKTsGFY/JwR+lWnYUgl8Hf9aqGvy5K6gEFLX8GPgVJdgSq6r4Agp+VTn0AaDsw0HwWsFb9uHAU/bhwqsE8PgDFT5YVCazTldhRma80Uhc2cB/W9nXthATAiz6yDl1l9+rkr01wP51fvav87NvnZ99a/3Ys/97UkRcioaURjpqNdJRq6GOxAwdCRlaEuvqSKyjJSY+8luzJcyFIGHuxKiBAJ7tOyhdvRrX6jWUrl2HUlJSYR9jo4ZYOnYktlNHYs/sIDNQazBVVXH6/ZSUddfYy7p0Sry+CrcPLsD5768r8w3MpNMRq9cTU37Rld8263WYdXpMOh0mvQ6zTodJp8OsD14bdTpM2uD1wXBh1OnKW6UM/woceo022Pp1MFCEqbtTVQ8Fi0MthSoB9VALol9Vyrf7FaVCa+O/Wx19geD1wRDiUxS8AQWfEqhw21sWUHyBQIXrw+8/+DwHb0dT8AwHDWDW6zGX/V7G6g3E6HXEGgzE6PXE6vVYDHosBgMWg6Hs6+CSGVaDgbjDblsNBgzVuFfEWaCQvTVA1hY/2VsCZG72k7MtQP7uQMiWvMOZLBriUjXYUrTEpWiJSw1eWxKDQS82QUtsgobYBA0miwZj7KHL6TqNmYS5ECTMHV3A4cC9cSOu9RtxrV+Pe9MmFGdphX30qSnEduhAbKeOWDqeiT45OUzViqriUxRKysbQ2MtCmf2wgGb3BMfalHh9FHs8weuysHay3TkQbM2yHv4mZDRiNejL35AsFS76Cm9isYbgG1yswUCMTodOuvEimv+wwOgtC4cHA2PwEsCvqEds/3cg9SkKAVUhoBwKsAFFxa8qh8YsqmUtmap6qIWTg62ZAOphLXHB+zTl7XTBljrKWvIOjlPUlo1l1Je1Bh7eEnuoBfFQS6JBW7Hb2agLfkAw6Q59mDj4dUzZh4pIbzGKdIqiUpypkLszQN7OYLgr3K9QdCBA4T6Fwv0KHufJ/73SGyG9lZ7HF1RtppClScRRKV4vnu3bcW/eUn7x7t4N/3ozNqSnEdOuHbHt2xHTvh2GtDT5IxMFDp6W5vBQdjBwlYe0wwJacdlK6HaPF+cpdFPG6PXYjIayBTeDXTlxZV8ffttqOLQgZ3ABzmCLg/xu1QwHg06MvBuJKqLVakjMCHatNut95P2qquK2q5TkKthzlOB1rkJJrkppkUJpoYqrWKG0KHjtdqr4SsFTquJ1qvi9oIR5RIf896lh/EVFeHbswLN9J54d2/Fs34Fn5y74d7eWXo+5WVNiWrciplUrzK1aYUipFZaaRcXzBJZ3SR5sHTvs9qHZaYfCWon3xE5LczitRlMeuoLjZ4Kzzmymw8balH1tMx5+qd7dOUKI6kOj0RATryEmXkvqGSf2WFVV8XuC4/jCScJcNaQqCv78fLx79uLds6fsshfP7j0ECgqOfIBWi7FhQ8zNm2Ju1gxzs2aYmp6B1mg8/cVXQweD2MEB/IcG8ftx+A6dqPmIAf9eLyU+X3lwO5Vuy4OnpSkf5Pzv8GUyEm80lYUyA/HG4Gw1i8GAVlrIhBAiJI1Gg8EMBnN4/05GbJhzOpx8O3Yim9ZuwRJn4YKRQ+jcs9MR+6mqys8Tf2XRvKUA9OzXjQsuG1atu2hURSFQXIw/Nxdfbh7+7Gx8BzLxHsjElxm8qN7Qoz01MTGYGjfC1Lhx8NKkEeYmTdDGxJzm7yLyqGpw3M3hy1WUL2VRtiyA2x/cXnpwaQCfn1L/oRmUB2dVlpZtP7itMgZ5m3Q6rAZDhfMExh1xzsBDLWOHt5aZ9RH7X10IIcQpiti/8JO+mopOr+elD55j3+79jHljLBn1M0ivm1Zhv4V/LmbNynU8+uKDaIAPXv2Y5JRkeg/oGZ7CT5CqKKgeD0ppKQGHA8XhLLt2ECi2Eyguwl9URKCwiEBREf78fPx5+ag+31GfVxcfj7F+fYz165VdgrcNtWujOc0Dwg/OrCsfgFy+gGpwdl1AVcsGLVecaRcov11xtl1wsHOgwow7byB438Gp/t7y2XTBKf8HZ88dvn5YMLQFbx+8VNXMuoNBLLZsEL/FYCgf4F9+22AILhdgMJYP/C9fQsBgwCjdlkIIIUKIyDDncXtYvXwNj7/8ECaziSbNG9O2Y2uWLVzBhZcNq7DvsgUr6D/4LBKTEgDoP7gfi+YuCXuYm/rLNCbsz6RsaXcOzZJSKVucCNTgKvDHpDVAckrw0qRpcJtOh0avR6PXg16PxmBAYzh4bYCywKaqZStU5RWg5h7sYlUPrV+lHlo/icNuH9x+cB9FPbg9eF1hNtjB9a7KHnhwravAYftGC4NWW75chblsqYpDy1gEl7Uw64JLARxcFiCm7LZFHwxrlsO+Pji7UhZGFUIIUVUiMszlZOWi1WlJTU8t35ZRrw7bNm0/Yt/M/Vlk1K9zaL/6GWTuzw75vAvnLGbh3MUADLnqUqjCpUkKnKXsiIursucPyecPXnCd3tc9DhqoMG3/4NR9fdlCnvrDFvTUH7ZO1+HT+g3/WiT04ErywdXqyxYaLbt9cIX5g6vPm3QVlwAw6rSH1hc7bF0xCV1CCCGiTUSGOY/Hizmm4vk6zbFm3G7Pkfu6PZhjzRX287g9qKp6xLi5Xv170Kt/DyC4XktVOr9fHzpnZYFWg0ajK7vWgk5b3nqmMRmDrWuc4Pi+ELuHeo7D10b69/0VVpsPcVtDxTWVNOWrox9aOV0D5Yui/vsckcH1lw59XZ3HMAohhBDhFJFhzmQy4na5K2xzuzyYzUeeUcBkNlXY1+1yYzKbwh4eaqenUTs97dg7CiGEEEKcgojsU0pNS0EJKORk5ZZv27/nAGl1jwxH6Rlp7N9zoMJ+6Rm1T0udQgghhBDhFpFhzmQ20b5zW6ZNmYnH7WHHlp2sXbWOrr06H7Fv196d+XPmPIoKiiguLGbOjLl069M1DFULIYQQQpx+EXtuVqfDyfhPJ7J53RYscbFcMHIonXt2YtvmHXz0+ie8MfYVIDiT8qfvfmXxvCUA9OjXnQsvP/Y6c3JuViGEEEJEi6PllogNc1VNwpwQQgghosXRcktEdrMKIYQQQojjI2FOCCGEECKKSZgTQgghhIhiEuaEEEIIIaKYhDkhhBBCiCgmYU4IIYQQIopF5Om8TgefP1Dl52cFcNgdWG3WKn8dcfzkmEQeOSaRSY5L5JFjEplOx3Hx+QP/eV+NDXNN66Wcltd57aMvefj/7j8tryWOjxyTyCPHJDLJcYk8ckwiU7iPi3SzCiGEEEJEMQlzQgghhBBRTMJcFet1Vo9wlyD+RY5J5JFjEpnkuEQeOSaRKdzHpcaem1UIIYQQojqQljkhhBBCiCgmYU4IIYQQIorV2KVJqprT4eTbsRPZtHYLljgLF4wcQueencJdVrU2748FLF2wnMy9mXTs3pFrRl9Rft/m9VuY9NVUCvMLadikPlffcgVJtZIA8Pn8TPpyMv8sW43BZOScoWfTf/BZYfouqpeDP9vN67dS6iylVmoy548cSuv2LQE5LuHy1Ufj2LJ+K16Pl7gEG+cMPZueZ3UH5JhEgpysXF5+/HU6dGnHdbddDcCKRSv5edJ0nCVOmrdpxlU3X4bFagHk/aYqvfPiB+zavhutNtj2lZAYz1OvPwZE1jGRMFdFJn01FZ1ez0sfPMe+3fsZ88ZYMupnkF43LdylVVvxCfGcd8G5bFq7Ga/XV77dUeJg7DtfcuVNI2lzZmumTZnBF+9/zQPP3gvAjKkzycnK5bm3n8JeVMK7L39IWkZtWrVrGabvpPpQAgESkxO454k7SExOYMPqjXzx/tc89tJDmMxGOS5hMvD8c7hy1OUYDHqyDmTz7ksfUrdBBkm1EuWYRIDvv5pC/Ub1yr/O3JfFd19M5tYHRlGvYV0mfD6JSV9O4YY7rwXk/aaqjbj2kvIPOwdF2jGRbtYq4HF7WL18DcOGD8JkNtGkeWPadmzNsoUrwl1atdahSzvad26LxRpbYfvq5WtJz0jjzG4dMBgNDL74PPbvOUDWgWwAlv61gkEXDSTWEktaRm16ntWdpfOXh+NbqHZMZhNDLhlEckoSWq2WNme2Jjklib279spxCaP0umkYDMHP8hqNBg2Ql5MvxyQCrFz8NzGxMTRv3bR82/JFK2lzZivOaNEEk9nE0OGDWb1iLW6XW95vwiTSjomEuSqQk5WLVqclNT21fFtGvTpk7csKY1U1V+b+LDLq1yn/2mQ2USu1Fln7syh1lmIvsle4P6N+HTL3y7GqCvbiEnKycknLSJPjEmYTv5zM/Tc9wgsPv4ItwUbr9i3lmISZy+Vm2tSZXHzVhRW2Z/3ruKTUroVOryMnK1feb06DXyZN49HbnuLN/3uXrRu3AZF3TKSbtQp4PF7MMeYK28yxZtxuT5gqqtk8bs8R58wzx5pxuzx4yo5JzGHHKybGXL5dVJ6AP8BXH42jW+/OpNWpLcclzC67/lJGXHsJO7fuYuvG7ej1ejkmYTZt8gx69OtKYlJChe0et5eYmJgK22Jigz97jVYr7zdV6MLLhpGWURudXs+qJX/z8Zuf8cgLD0TcMZGWuSpgMhlxu9wVtrldHsxmU5gqqtlMZlOI4+HGHGPCVHZMDr/f7XaXbxeVQ1EUvh4zHr1Ox4hrhwNyXCKBVqulSfPGFBUWsWD2QjkmYbRv9342r9/C2YP6HXGfyRz6PcVkNsn7TRVreEYDzDFmDAY93fp0oXHThmxYvTHijom0zFWB1LQUlIBCTlYuqWkpAOzfc4A0GYwaFukZaSz969C4Ho/bQ15OPmkZacRaYrEl2Ni/5wAt2jYHgscqPUOOVWVRVZVvx06kxF7CrQ/ejE6vA+S4RBIloJCXky/HJIy2btxGQW4hT9/7PBD82auKwqv736Bluxbs33OgfN+8nHz8Pj+paSloNBp5vzmdNBpUVSUtIy2ijom0zFUBk9lE+85tmTZlJh63hx1bdrJ21Tq69uoc7tKqtUAggM/rQ1EUVFXB5/URCARo17ktmfuy+Gf5anxeHzN//J2Meumk1akNQNfenfntpz8odZaSdSCbRX8uoVvfLmH+bqqPiV9OJvtANqPvH4XRaCzfLsclPEqKS1i5+G88bg+KorBxzSZWLv6b5q2byjEJo15n9+CZNx7n0Rce4NEXHqB3/5607tCK2x8eTZeenVj393q2bd6Bx+1h2pQZtO/cFnOMWd5vqlCp08XGNZvK30uWL1zJ9k07aNWuRcQdEzmdVxVxOpyM/3Qim9dtwRIXywUjh8q6P1Vs+tSZzPjh9wrbBl88kCGXDGLTui18//VUCvMKaNCkAVffcgXJKSHWzjIaOGdYf1k7q5IU5BXwzH0voDfoy9dpArj8hhF06dVJjksYlNgdfP7ul+zfewBVUUmslUi/gX3odXbw3JJyTCLD9Kkzyc3O+9c6c9NwlpTSvE1Trrr58gprmsn7TeUrsTsY879Pyc7MQavVUDs9laHDB5e3TEfSMZEwJ4QQQggRxaSbVQghhBAiikmYE0IIIYSIYhLmhBBCCCGimIQ5IYQQQogoJmFOCCGEECKKSZgTQgghhIhiEuaEEEIIIaKYhDkhhBBCiCgmYU4IIYQQIopJmBNCCCGEiGIS5oQQQgghopiEOSGEEEKIKCZhTgghhBAiikmYE0IIIYSIYhLmhBBCCCGimIQ5IYQQQogoJmFOCCGEECKKSZgTQgghhIhiEuaEEDXKkvnLeGDUo2F7/VJnKY/f8TS52XmV8nw+n5+n7/0/9uzYWynPJ4SIPhpVVdVwFyGEEJXhrmvuP+r9XXt34bIbhuNxeYiLjztNVVX044SfcZQ4ufqWKyrtOef9voA1K9dx12O3VdpzCiGihz7cBQghRGV58b1ny2+v+2cDEz6bVGGbwWjAaDRiNBpPf3GA1+Nl0dyljL7/pkp93s49O/LjhJ/J3JdFet20Sn1uIUTkkzAnhKg2bAm28tsxsTFHbINgN+v3X0/ljbGvADB96kz+WbaGAUPPZvrUmTjsTs7s1p7LbxzBorlL+eOX2Xi9Xrr17sJFV5yPVhscneL3+5k2eQYrFq3C6SwlPSONYZcOpmW7Fv9Z3/rVG9FooHGzRuXbtm7cxrsvfcidj97KL5Omc2BfJmkZtbnixpHUa1gXAFepi++/nsrGNZtxu93EJ9joN7APZw/qB4DFaqFR00asXLyKYSOGVMJPUggRTSTMCSFqvPy8AtasWsfo+0dRXFjM2He/xF5kx5Zg446HR5Odmc3n739N42YN6dClPQDjP/mOvJw8rrv9ahKSEli/egMfv/kZDz53L3UbZIR8ne2bd1CvYT00Gs0R9/08aRoXXjYMW4KNKeN+5KuPxvHEK4+g0Wj4dfIMDuzNZPQDN2GLjyM/twCH3VHh8Q2a1Gfrpu2V/8MRQkQ8mQAhhKjxVEXh6psvp069dFq2a0Grdi3Yu2sfl984grSM2rTv3I7GTRuxZcM2AHKz81i55G9uuPM6zmjRhFqpyfQ7tw+t2rdk4Z+L//N1CvIKiU+0hbxv2PDBNGvVlLQ6tRl80UCyD+RQVFhc/rh6DevSsEkDkmol0bTlGZzZrUOFx8cn2CjIK6ycH4gQIqpIy5wQosZLTE4s75YFiLPFkZKWgl5/6E9kXHxceWvYvl37UFWVFx99tcLz+P1+mrVq+p+v4/P5sBmsIe+rU79O+e2DXcMOu4PEpAR6D+jJ5+9+xZ6de2nRpjltzmxF05ZnVHi8wWjA5/Ud53cshKhOJMwJIWo8nU5XcYPmyG0aQCmb/K+oKhqNhoeeuw+dvmIHh8Fg+M/XsVotlDpd/1HDoec52A2rKsHXa92+Jc+9/SQbVm9i8/qtjHljLGd2bV9hRmypsxSrzXL0b1QIUS1JN6sQQpygeg0yUFUVe7GdlNopFS4JSQn/+bi6DTLI2p99Uq9pjbPStXdnrhl9BVeOuoxlf63A5/OX35+5L4t6Deqe1HMLIaKbhDkhhDhBqempdO7ZkXGffMffy1aTl5PPnh17mT3tT/5ZvuY/H9eyXXOyDmTjLHGe0OtNmzKD1SvWkpOVS9b+bFavWEtyShIGw6HOle2bdxx1Jq0QovqSblYhhDgJV998Bb/9/Ac/ffcLRQXFxFpjadC4Pk1bnfGfj6lTrw4NmtRn5ZK/6Xtu7+N+Lb1ez6+Tp5OfW4DBoKdhkwaMvn9U+f07t+7CVeqmQ9d2p/Q9CSGik5wBQgghTqMNazYy5ZsfeeLVR8rXrDtVn737FXUbZnDeBedUyvMJIaKLtMwJIcRp1KpdS3LOyaWooIikWkmn/Hw+n5+M+umcPahvJVQnhIhG0jInhBBCCBHFZAKEEEIIIUQUkzAnhBBCCBHFJMwJIYQQQkQxCXNCCCGEEFFMwpwQQgghRBSTMCeEEEIIEcX+H5y0vmDG9tPRAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define system parameters.\n",
    "chi = 2 * np.pi * 3 * 1e6  # Hz\n",
    "omega_max = 2 * np.pi * 1e6  # Hz\n",
    "total_rotation = np.pi\n",
    "\n",
    "# Define pulse using pulses from Q-CTRL Open Controls.\n",
    "pulse = oc.new_primitive_control(\n",
    "    rabi_rotation=total_rotation,\n",
    "    azimuthal_angle=0.0,\n",
    "    maximum_rabi_rate=omega_max,\n",
    "    name=\"primitive\",\n",
    ")\n",
    "\n",
    "sample_times = np.linspace(0, pulse.duration, 100)\n",
    "target_operator = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])\n",
    "\n",
    "try:\n",
    "    import qutip as qt\n",
    "\n",
    "    # Define matrices for the Hamiltonian operators.\n",
    "    a = qt.destroy(3)\n",
    "    n = qt.num(3)\n",
    "    ad2a2 = n * n - n\n",
    "\n",
    "    # Define the target.\n",
    "    target_operation = qt.Qobj(target_operator)\n",
    "\n",
    "    # Define the initial state.\n",
    "    initial_state = qt.basis(3)\n",
    "\n",
    "    graph = bo.Graph()\n",
    "\n",
    "except ModuleNotFoundError:\n",
    "    graph = bo.Graph()\n",
    "\n",
    "    # Define matrices for the Hamiltonian operators.\n",
    "    a = graph.annihilation_operator(3)\n",
    "    n = graph.number_operator(3)\n",
    "    ad2a2 = n @ n - n\n",
    "\n",
    "    # Define the target.\n",
    "    target_operation = target_operator\n",
    "\n",
    "    # Define the initial state.\n",
    "    initial_state = graph.fock_state(3, 0)[:, None]\n",
    "\n",
    "    pass\n",
    "\n",
    "# Define the anharmonic term.\n",
    "anharmonic_drift = 0.5 * chi * ad2a2\n",
    "\n",
    "# Define Rabi drive term.\n",
    "rabi_rate = graph.pwc(\n",
    "    durations=pulse.durations,\n",
    "    values=pulse.rabi_rates * np.exp(1j * pulse.azimuthal_angles),\n",
    ")\n",
    "rabi_drive = graph.hermitian_part(rabi_rate * a)\n",
    "\n",
    "hamiltonian = anharmonic_drift + rabi_drive\n",
    "\n",
    "# Calculate the time-evolution operators.\n",
    "unitaries = graph.time_evolution_operators_pwc(\n",
    "    hamiltonian=hamiltonian, sample_times=sample_times, name=\"time_evolution_operators\"\n",
    ")\n",
    "\n",
    "# Calculate the infidelity.\n",
    "infidelity = graph.unitary_infidelity(\n",
    "    unitary_operator=unitaries[-1], target=target_operation, name=\"infidelity\"\n",
    ")\n",
    "\n",
    "evolved_states = unitaries @ initial_state\n",
    "evolved_states.name = \"evolved_states\"\n",
    "\n",
    "# Run simulation.\n",
    "graph_result = bo.execute_graph(\n",
    "    graph=graph,\n",
    "    output_node_names=[\"evolved_states\", \"infidelity\", \"time_evolution_operators\"],\n",
    ")\n",
    "\n",
    "# Extract and print final infidelity.\n",
    "print(\n",
    "    f\"Noise-free infidelity at end: {graph_result['output']['infidelity']['value']:.1e}\"\n",
    ")\n",
    "\n",
    "# Extract and print final time evolution operator.\n",
    "print(\"Time evolution operator at end:\")\n",
    "print(np.round(graph_result[\"output\"][\"time_evolution_operators\"][\"value\"][-1], 3))\n",
    "\n",
    "# Extract and plot state populations.\n",
    "state_vectors = graph_result[\"output\"][\"evolved_states\"][\"value\"]\n",
    "qv.plot_population_dynamics(\n",
    "    sample_times,\n",
    "    {rf\"$|{state}\\rangle$\": np.abs(state_vectors[:, state]) ** 2 for state in range(3)},\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
