optimization2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: -7.573357607487814e-05\n",
      "            Iterations: 3\n",
      "            Function evaluations: 12\n",
      "            Gradient evaluations: 3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHFCAYAAAANG6v4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABpM0lEQVR4nO3dd1hT1x8G8DfsjQoKDkTcu65qXVXcq63bWkeddWAVN05w4kIcta5abd1VW0e1rjpqrbRurVqrVVyIW3Ayz++P8yMYGQYFTsb7eZ48XJLLzZuQ5Jt77xkaIYQAERERGRUL1QGIiIgo41jAiYiIjBALOBERkRFiASciIjJCLOBERERGiAWciIjICLGAExERGSEWcCIiIiPEAk5ERGSEWMBTERYWhnbt2iFv3rywsbGBp6cn2rZtiyNHjrzTdr/++musWLEixfXh4eHQaDSp3va2smKb+lixYgU0Gg3Cw8PTXS8oKAgajQb3799P9fayZcuibt26mR/QROzatQuNGjVCvnz5YGtri3z58qFu3bqYNm2aznoajQYDBgxQlPLd/ffff7C1tdV573Xr1g0ajUZ7sbGxQZEiRTBs2DBER0e/1f2sWbMGc+bMeaesab2/M4NGo0FQUNAb15szZw5at24NHx8faDSadN9D+/fvR8OGDZEnTx44OTmhfPnymDdvHhISElKsu3fvXlSvXh0ODg5wd3dHt27dcPfu3bd+PKo+n1R69OgRcuTIgc2bN2faNlnAXzN//nzUrFkTN2/exIwZM7B3717MmjULt27dQq1atfDVV1+99bbTeoPnzZsXR44cQfPmzd8hedZvkwzDokWL0KRJE7i4uOCrr77Crl27MH36dJQqVQobN25UHS9TDRs2DA0bNkT16tV1rre3t8eRI0dw5MgRbN26Fb6+vggJCUHbtm3f6n4MvYDra9GiRbh27Rrq1auH3Llzp7ne3r170aBBA8THx2Pp0qXYvHkz6tati0GDBmHIkCE66x48eBBNmzaFh4cHtmzZgrlz52Lv3r2oX78+YmJisvohmYycOXNi8ODBGD58OGJjYzNno4K0fv/9d2FhYSFatGgh4uLidG6Li4sTLVq0EBYWFuL3339/q+2XKVNG1KlTJxOSGq7ly5cLAOLq1avprhcYGCgAiHv37qV6uzk8V2+rYMGC4sMPP0z1toSEBJ3fAQg/P7/siJXpzp8/LwCInTt36lz/+eefC0dHxxTr+/r6CgDiypUrGb6v5s2bC29v77eNKoTI2tcsABEYGPjG9V79/6eXp1OnTsLW1lY8ffpU5/pGjRoJFxcXnevef/99Ubp0aZ3PxMOHDwsA4uuvv9b/Qbzi6tWrAoBYvnz5W/39u3j+/Hm232eSyMhIYWVlJVavXp0p2+Me+CuCg4Oh0WiwcOFCWFlZ6dxmZWWFr7/+GhqNRucwZdKh4JMnT6J169ZwcXGBq6srOnfujHv37mnXK1SoEM6dO4eDBw9qD/0VKlQIQOqHk5K2e+bMGbRr1w6urq7IlSsXhgwZgvj4eFy8eBFNmjSBs7MzChUqhBkzZujkTW2brx52fP3y6iHvY8eO4eOPP0auXLlgZ2eHihUr4ocffkjxfIWFhaFmzZqws7NDvnz5MGrUKMTFxb3FM6+f2NhYTJ48GSVLloStrS1y586N7t276zzP06ZNg4WFBbZt26bzt926dYODgwPOnj0LAHj58iWGDh2KChUqaJ/b6tWrY8uWLSnu9/Hjx+jZsydy5coFJycnNG/eHFeuXEn1sObvv/+O+vXrw9nZGQ4ODqhRowa2b9+us07SaYb9+/ejX79+cHd3h5ubG1q3bo2IiIg3Pg8PHjxA3rx5U73NwiL1t/TKlStRqlQpODg44L333sPPP/+sc/vly5fRvXt3FCtWDA4ODsifPz8++ugj7fOV5MCBA9BoNFi1ahWGDBkCT09P2Nvbo06dOjh58mSK+9X3tZSahQsXwtPTEw0bNtRr/SpVqgAA7ty5o70uMTERM2bM0L5m8uTJg65du+LmzZvaderWrYvt27fj2rVrOu+JJBMmTEC1atWQK1cuuLi4oFKlSli2bBnEK/NApff+BoDo6GgMGzYMPj4+sLGxQf78+eHv749nz57pPIbo6Gj07t0bbm5ucHJyQpMmTfDvv//q9fiBtP//r7O2toaNjQ3s7e11rs+RIwfs7Oy0v9+6dQtHjx5Fly5ddD4Ta9SogeLFi+Onn356431FRESgffv2cHZ2hqurKzp06IDIyMhU19X39fL777+jevXqsLOzQ/78+TFu3Dh88803KT7LChUqhBYtWuDHH39ExYoVYWdnhwkTJgAAIiMj0adPHxQoUAA2Njbw8fHBhAkTEB8fr3Nf+nzuAMC+fftQt25duLm5wd7eHgULFkSbNm3w/Plz7ToeHh5o2LAhFi1a9MbnTS+Z8jXABMTHxwsHBwdRrVq1dNerWrWqcHBwEPHx8UKI5D1Jb29vMXz4cLFr1y4xe/Zs4ejoKCpWrChiY2OFEEKcOHFCFC5cWFSsWFEcOXJEHDlyRJw4cUIIkfq30aTtlihRQkyaNEns2bNHjBgxQgAQAwYMECVLlhTz5s0Te/bsEd27dxcAxKZNm7R/n9o2k+436bJv3z6RP39+4enpKaKiooQQQuzbt0/Y2NiI2rVri/Xr14udO3eKbt26pdjWuXPnhIODgyhdurRYu3at2LJli2jcuLEoWLBghvbAIyMjRVxcXIrL63sPCQkJokmTJsLR0VFMmDBB7NmzR3zzzTcif/78onTp0tpv1YmJiaJZs2YiZ86cIjw8XAghxLfffisAiG+++Ua7vcePH4tu3bqJlStXin379omdO3eKYcOGCQsLC/Hdd9/p3G+tWrWEnZ2dmDZtmti9e7eYMGGCKFasWIq9ogMHDghra2tRuXJlsX79erF582bRqFEjodFoxLp167TrJR2lKFy4sPjyyy/Frl27xDfffCNy5swpfH19033ehBCiQYMGwsrKSgQGBopTp05pX4upASAKFSokqlatKn744QexY8cOUbduXWFlZSX+++8/7XoHDx4UQ4cOFRs3bhQHDx4UP/30k2jZsqWwt7cX//zzj3a9/fv3CwDCy8tLfPLJJ2Lbtm1i1apVomjRosLFxUVnm/q+ltJSuHBh0b59+xTXp7UH3rZtW2FlZSXu3Lmjve6LL77Qvmd27twpFi1aJHLnzi28vLy0R3/OnTsnatasKTw9PXXeH0m6desmli1bJvbs2SP27NkjJk2aJOzt7cWECRO066T3/n727JmoUKGCcHd3F7NnzxZ79+4Vc+fOFa6urqJevXoiMTFRCCFfu76+vsLW1lZMmTJF7N69WwQGBorChQvrvQf+qvT2wMPCwoStra3w8/MTt27dEo8ePRLff/+9sLa2FrNmzdKut3PnTgFAbN++PdXnO2/evOlmeP78uShVqpRwdXUV8+fPF7t27RIDBw7Ufk68+jrQ9/Vy+vRpYWdnJ8qXLy/WrVsntm7dKpo1ayYKFSqU4rPH29tb5M2bVxQuXFh8++23Yv/+/eKvv/4St2/fFl5eXsLb21ssXrxY7N27V0yaNEnY2tqKbt26af9e38+dq1evCjs7O9GwYUOxefNmceDAAbF69WrRpUsX8ejRI53nZPr06cLCwiLF9W+DBfz/IiMjBQDx6aefprtehw4dBADth0RSIRo8eLDOeqtXrxYAxKpVq7TXpfWGSq+Ah4SE6KxboUIFAUD8+OOP2uvi4uJE7ty5RevWrdPd5qvi4+PFJ598IpycnMTx48e115csWVJUrFgxxSmEFi1aiLx582oP0XXo0EHY29uLyMhInW2WLFkyQwU8vcurz9XatWtTfEkRQoijR4+mOJR3//59UaBAAVG1alVx4sQJ4eDgIDp37pxunvj4eBEXFyd69uwpKlasqL1++/btAoBYuHChzvrBwcEpPlQ/+OADkSdPHvHkyROd7ZYtW1YUKFBA+0GdVMD79++vs80ZM2YIAOL27dvpZr18+bIoW7as9nmyt7cX9evXF1999ZX2C2MSAMLDw0NER0drr4uMjBQWFhYiODg43ecjNjZWFCtWTOe1nVTAK1WqpH08QggRHh4urK2tRa9evbTX6ftaSs2dO3cEADFt2rQUtyUV8KQve/fv3xcLFy4UFhYWYvTo0dr1Lly4kOrz/OeffwoAOuvqewg9ISFBxMXFiYkTJwo3Nzed5yCt93dwcLCwsLAQR48e1bl+48aNAoDYsWOHEEKIX375RQAQc+fO1VlvypQpmV7AhZCHwfPly6d9HVlaWooZM2borJP0OfbqF5okX3zxhbCxsUk3w8KFCwUAsWXLFp3re/funeLzSd/XS7t27YSjo6PO6beEhARRunTpVAu4paWluHjxos42+/TpI5ycnMS1a9d0rp81a5YAIM6dOyeE0P9zJ+l/eerUqXSfDyGE2LNnjwAgfvnllzeu+yY8hJ5B4v+HzV49xAYAnTp10vm9ffv2sLKywv79+9/p/lq0aKHze6lSpaDRaNC0aVPtdVZWVihatCiuXbum93YHDBiA7du3Y8OGDahUqRIAeRj1n3/+0T6W+Ph47aVZs2a4ffs2Ll68CEC2YK1fvz48PDy027S0tESHDh0y9Pj27t2Lo0ePprgUKVJEZ72ff/4ZOXLkwEcffaSTq0KFCvD09MSBAwe067q5uWH9+vU4ceIEatSogYIFC6Z6yGrDhg2oWbMmnJycYGVlBWtrayxbtgwXLlzQrnPw4EEA8v/5qo4dO+r8/uzZM/z5559o27YtnJycdJ6TLl264ObNm9rnLsnHH3+s83v58uUB4I3/xyJFiuD06dM4ePAgJkyYgAYNGuDo0aMYMGAAqlevjpcvX+qs7+vrC2dnZ+3vHh4eyJMnj879xMfHY+rUqShdujRsbGxgZWUFGxsbXLp0Sef5SPLZZ5/pvAe8vb1Ro0YN7es9I6+l1CSdSsiTJ0+qtz979gzW1tawtraGu7s7+vXrhw4dOmDKlCnadZKydOvWTedvq1atilKlSuHXX39N8/5ftW/fPjRo0ACurq6wtLSEtbU1xo8fjwcPHujVEvvnn39G2bJlUaFCBZ3noXHjxtBoNNrXblLe1z9LPvvsM71yZsTx48fRqlUrVK5cGdu2bcO+ffswatQojB07FpMmTUqx/uufd2+6Psn+/fvh7Oyc4rX++mPKyOvl4MGDqFevHtzd3bV/b2FhkeI9mqR8+fIoXry4znU///wzfH19kS9fPp37SvpcTXrf6/u5U6FCBdjY2OCLL77Ad999hytXrqT5nCS9pm/dupXuc6cPqzevYh7c3d3h4OCAq1evprteeHg4HBwckCtXLp3rPT09dX63srKCm5sbHjx48E65Xr8fGxsbODg46JynSrpe3y40kydPxqJFi7Bs2TI0adJEe33SucNhw4Zh2LBhqf5tUrevBw8epHjMQMrn4U3ee+89nTdiktcf3507d/D48WPY2NikmytJtWrVUKZMGZw+fRr9+vWDo6Ojzu0//vgj2rdvj3bt2mH48OHw9PSElZUVFi5ciG+//Va73oMHD2BlZZXi//DqFxdAdhERQqR6bjpfvnzabb3Kzc1N53dbW1sAwIsXL1J9jK+ysLDAhx9+iA8//BCALGg9e/bE+vXr8e2336J///5p3k/Sfb16P0OGDMGCBQswcuRI1KlTBzlz5oSFhQV69eqVap60/venT58GkLHXUmqS7vP110ESe3t7/PbbbwDkucyQkBCsXbsW5cuXR0BAAIDk5zut/4k+X3j/+usvNGrUCHXr1sXSpUu150s3b96MKVOm6PW/unPnDi5fvgxra+tUb3/1PZX0ufGqjL6n9OHn5wcPDw/89NNPsLS0BCC/6FlYWCAoKAidOnVC4cKFtVlS+xx7+PBhivfF6x48eJDivQKkfEwZ/exJbZupXQek/v+/c+cOtm3b9sb/ib6fO0WKFMHevXsxY8YM+Pn54dmzZyhcuDAGDhyIQYMG6fxN0mtan9fOm7CA/5+lpSV8fX2xc+dO3Lx5EwUKFEixzs2bN3H8+HE0bdpU+6JPEhkZifz582t/j4+Px4MHD1L98FRpxYoVGDduHIKCgtCjRw+d25IK6ahRo9C6detU/75EiRIAZFFIrSFKWo1T3lVSQ6+dO3emevure5gAEBgYiLNnz6Jy5coYP348WrRogcKFC2tvX7VqFXx8fLB+/XqdvYjXu8W4ubkhPj4+xYfV648zqeDdvn07RbakvcnUvqhkFkdHR4waNQrr16/H33//neG/X7VqFbp27YqpU6fqXH///n3kyJEjxfpp/e+TXu8ZeS2lJunvHz58mOrtFhYW2kZrANCwYUNUrlwZEyZMQKdOneDl5aXNcvv27RTv54iICL3+H+vWrYO1tTV+/vlnnS8TGenL6+7uDnt7e50vhq/fDiS/1l7/3MiK99SpU6fQsWPHFJ9j77//PhITE3HhwgUULlwYZcuWBQCcPXsWzZo101n37Nmz2tvT4ubmhr/++ivF9a8/pox+9rzaUDGtbSZJ7SiBu7s7ypcvr3PE5lVJX7oz8rlTu3Zt1K5dGwkJCTh27Bjmz58Pf39/eHh44NNPP9Wul/SazozPAx5Cf8WoUaMghED//v1TDGaQkJCAfv36QQiBUaNGpfjb1atX6/z+ww8/ID4+Xmcghdf3erLbzp070bt3b/To0QOBgYEpbi9RogSKFSuG06dPo0qVKqlekl6wvr6++PXXX3XeSAkJCVi/fn2WZG/RogUePHiAhISEVHO9Wgz27NmD4OBgjB07Fnv27NG2fH2172XSACCvvrkjIyNTtEKvU6cOAKR4XOvWrdP53dHREdWqVcOPP/6o8z9OTEzEqlWrUKBAgRSH8d5Wal8SAGgPdSd9+GSERqPRHgFIsn379jQP861du1anFfa1a9fwxx9/aF/vGXktpcbb2xv29vb477//9Mpva2uLBQsW4OXLl5g8eTIAoF69egDkl5NXHT16FBcuXED9+vV1/j6196ZGo4GVlZVOoXvx4gVWrlyZaobUttGiRQv8999/cHNzS/V5SGqt7uvrCyDlZ8maNWv0eQoyJF++fDh27FiKz7mkAXOSvvDkz58fVatWxapVq3TWDQsLw8WLF9Mstkl8fX3x5MkTbN26Vef61x9TRl4vderUwb59+3SO4CQmJmLDhg16P/4WLVrg77//RpEiRVK9r6T3UEY+d5JYWlqiWrVqWLBgAQDgxIkTOrcnHV4vXbq03nnTwj3wV9SsWRNz5syBv78/atWqhQEDBqBgwYK4fv06FixYgD///BNz5sxBjRo1Uvztjz/+CCsrKzRs2BDnzp3DuHHj8N577+mclylXrhzWrVuH9evXo3DhwrCzs0O5cuWy5bFdvXoV7dq1Q+HChdG9e3eEhYXp3F6xYkXY2tpi8eLFaNq0KRo3boxu3bohf/78ePjwIS5cuIATJ05o3yRjx47F1q1bUa9ePYwfPx4ODg5YsGBBim4xmeXTTz/F6tWr0axZMwwaNAhVq1aFtbU1bt68if379+OTTz5Bq1atcPv2bXTu3Bl16tRBYGAgLCwssH79enz44YcYMWKEdrCOpK4l/fv3R9u2bXHjxg1MmjQJefPmxaVLl7T326RJE9SsWRNDhw5FdHQ0KleujCNHjuD7778HoNttJzg4GA0bNoSvry+GDRsGGxsbfP311/j777+xdu3aN54v1FeZMmVQv359NG3aFEWKFMHLly/x559/IiQkBB4eHujZs2eGt9miRQusWLECJUuWRPny5XH8+HHMnDkz1SNRAHD37l20atUKvXv3RlRUFAIDA2FnZ6fz5Vbf11JqbGxsUL169RSv0/TUqVMHzZo1w/LlyxEQEIASJUrgiy++wPz582FhYYGmTZsiPDwc48aNg5eXFwYPHqz923LlyuHHH3/EwoULUblyZe0efvPmzTF79mx89tln+OKLL/DgwQPMmjUrxZedpG2k9v729/fHpk2b8OGHH2Lw4MEoX748EhMTcf36dezevRtDhw5FtWrV0KhRI+3r9NmzZ6hSpQoOHz6c6peFtBw7dkzbjSo6OhpCCO3gPu+//z68vb0BAIMHD8bAgQPx0UcfoU+fPnBwcMCvv/6KkJAQNGjQAO+99552m9OnT0fDhg3Rrl079O/fH3fv3kVAQADKli2L7t27p5una9euCA0NRdeuXTFlyhQUK1YMO3bswK5du1Ksq+/rZcyYMdi2bRvq16+PMWPGwN7eHosWLdJ+9ujTlW7ixInYs2cPatSogYEDB6JEiRJ4+fIlwsPDsWPHDixatAgFChTQ+3Nn0aJF2LdvH5o3b46CBQvi5cuX2iMuDRo00LnvsLAwuLm5Zc5n/zs3gzNBR44cEW3bthUeHh7CyspK5MmTR7Ru3Vr88ccfKdZNak19/Phx8dFHHwknJyfh7OwsOnbsqNOdRQjZUrdRo0bC2dlZ2/VMiPRbob8+0ElaXWjq1KkjypQpo/399W0mtR5O6/Jqy83Tp0+L9u3bizx58ghra2vh6ekp6tWrJxYtWqRzn4cPHxYffPCBsLW1FZ6enmL48OFiyZIlWTaQS1xcnJg1a5Z47733hJ2dnXBychIlS5YUffr0EZcuXRLx8fGiTp06wsPDI0VL7pkzZwoA4qefftJeN23aNFGoUCFha2srSpUqJZYuXarN9aqHDx+K7t27ixw5cggHBwfRsGFDERYWlmqL4UOHDol69eoJR0dHYW9vLz744AOxbds2nXWSWqG/3io56X+0f//+dJ+7xYsXi9atW4vChQsLBwcHYWNjI4oUKSL69u0rbty4obMu0hjIxdvbW3z++efa3x89eiR69uwp8uTJIxwcHEStWrXEoUOHRJ06dXT+D0kZV65cKQYOHChy584tbG1tRe3atcWxY8dS3I++r6XULFu2TFhaWoqIiAid69N6DwghxNmzZ4WFhYXo3r27EEK2Tp4+fbooXry4sLa2Fu7u7qJz584pnqeHDx+Ktm3bihw5cgiNRqPzGvj2229FiRIlhK2trShcuLAIDg4Wy5YtS/E6T+v9LYQQT58+FWPHjhUlSpQQNjY2wtXVVZQrV04MHjxYpyfH48ePRY8ePXRea//884/erdA///zzNN/jr/dI2bRpk6hVq5Zwd3cXjo6OokyZMmLSpEkpBncRQojdu3eLDz74QNjZ2YlcuXKJrl27pvh8S8vNmzdFmzZttJ+Nbdq0EX/88UeqmfR9vRw6dEhUq1ZN57Nn+vTpAoB4/Pixdj1vb2/RvHnzVHPdu3dPDBw4UPj4+Ahra2uRK1cuUblyZTFmzBid5+BNnztCyJrRqlUr4e3tLWxtbYWbm5uoU6eO2Lp1q859JiYmCm9vb/Hll1/q9dy9CQv4O3pTISLTlNS95vDhw6qjZKukAr5hw4Ysv68XL16I3Llzp9qVjOh1DRs2FMWKFVMdI1179+4VFhYW4sKFC5myPR5CJ3qDtWvX4tatWyhXrhwsLCwQFhaGmTNn4sMPP0z1dApljqRRs4KCgjBgwIAUPQnIfA0ZMgQVK1aEl5cXHj58iNWrV2PPnj1YtmyZ6mjpmjx5Mnr06IGSJUtmyvZYwInewNnZGevWrcPkyZPx7Nkz5M2bF926ddM2lqKs88UXX+Dx48e4cuVKtrUXIcOXkJCA8ePHIzIyEhqNBqVLl8bKlSvRuXNn1dHS9OjRI9SpU0eni+e70gjxSlNSIiIiMgrsRkZERGSEWMCJiIiMEAs4ERGRETKrRmyJiYmIiIiAs7Nzpg2qQURElJmEEHjy5Any5cuX7sA0ZlXAIyIi4OXlpToGERHRG924cSPN0RABMyvgSWPp3rhxAy4uLorTEBERpRQdHQ0vL6905wsAzKyAJx02d3FxYQEnIiKD9qZTvWzERkREZIRYwImIiIwQCzgREZERYgEnIiIyQizgRERERogFnIiIyAixgBMRERkhFnAiIiIjxAJORERkhFjAiYiIjBALOBERkRFiASciIjJCZjWZCZmWhATg5k3VKYjUc3MDnJxUp6DsxgJORikxEahfHzh4UHUSIvWcnYHjx4FixVQnoezEAk5GadUqWbw1GsDWVnUaInXi44EnT4ChQ4GtW1WnoezEAk5G5+lTICBALgcHAyNHqs1DpNI//wDlygHbtgG7dwONGqlORNmFjdjI6AQHA7dvA0WKAP7+qtMQqVWyJDBggFwePBiIi1Obh7IPCzgZlStXgJAQuRwSwsPnRAAwfjzg7g6cPw8sWqQ6DWUXFnAyKsOHAzExQIMGwMcfq05DZBhy5gQmT5bL48cD9++rzUPZgwWcjMa+fcCPPwKWlkBoqGzARkRSr15A+fLA48dAYKDqNJQdWMDJKMTHJ5/v7tcPKFtWaRwig2NpCcydK5cXLQLOnlWbh7IeCzgZhW++kR9IOXMCQUGq0xAZprp1gTZt5DgJ/v6AEKoTUVZiASeD9+gRMHasXJ44UY46RUSpmzlTNu7ctw/YskV1GspKLOBk8CZMAB48AMqUAfr2VZ2GyLD5+ADDhsnloUOBly/V5qGswwJOBu38eeCrr+TynDmAFYceInqjgAAgXz7Z7XLOHNVpKKuwgJPBEgIYMkROWvLJJ7LrGBG9mZMTMH26XJ48WQ58RKaHBZwM1o4dwK5dgLU1MGuW6jRExuWzz4Bq1YBnz4DRo1WnoazAAk4GKTZWDgsJyJ9Fi6rNQ2RsLCySu5WtWAEcPao0DmUBFnAySPPnA5cuAR4ewJgxqtMQGadq1YCuXeXywIHsVmZqWMDJ4Ny5I7uLAXLiEhcXtXmIjFlwMODoCISFAWvWqE5DmYkFnAzO2LFAdDRQuTLw+eeq0xAZt3z5ko9ijRwpz4mTaWABJ4Ny8iSwbJlcnjtXnscjonczeLDsH37rVnLrdDJ+/HgkgyEEMGiQ/NmxI1CzpupERKbBzi65J8fMmUB4uNI4lElYwMlgbNgAHDoE2NtzL4Eos7VqBfj6ypHZRoxQnYYyAws4GYTnz+Vc34AcRcrLS20eIlOj0chR2Sws5JflgwdVJ6J3xQJOBmHWLOD6daBgweRxnIkoc5UvD/TpI5cHDZKjHJLxYgEn5W7cAKZNk8szZgAODmrzEJmyiROBHDmA06eTG4yScWIBJ+UCAoAXL4BatYD27VWnITJt7u5AUJBcHjsWePxYZRp6FyzgpNThw3JwCY1GdhvTaFQnIjJ9/fsDJUsC9+4BkyapTkNviwWclElMlOfhAKBnT6BSJbV5iMyFtXXyNKPz5gEXLyqNQ2+JBZyU+e474PhxOVTq5Mmq0xCZl8aNgRYtgPh4OW0vGR8WcFIiOhoYNUoujxsnJy0houwVEiL3xnfsAH75RXUayigWcFJi6lQ5aUmxYnKWJCLKfsWLJ7//Bg8G4uLU5qGMYQGnbHf5MhAaKpdnzwZsbNTmITJn48YBuXPL8+ALFqhOQxnBAk7ZbtgwIDZWnoNr3lx1GiLz5uoqj4gBsnvZvXtK41AGsIBTttq7F9iyBbC0lHvh7DZGpF737kDFikBUlNwjJ+PAAk7ZJj4e8PeXy35+QKlSSuMQ0f9ZWspxGABg6VI5ShsZPhZwyjaLFwPnzgFubskjQRGRYahdW46EmJgov2gLoToRvYnRFvDg4GBoNBr4J+3SkUF78CD50NykSUDOnGrzEFFKM2bIucMPHAB+/FF1GnoToyzgR48exZIlS1C+fHnVUUhPQUHAo0dAuXJA796q0xBRary9k+cKHzZMzlFAhsvoCvjTp0/RqVMnLF26FDm5G2cUzp0DFi6Uy3PmAFZWSuMQUTpGjAAKFADCw2U3TzJcRlfA/fz80Lx5czRo0EB1FNKDEHKAiIQEoFUroF491YmIKD2OjsD06XJ56lTg1i21eShtRlXA161bhxMnTiA4OFiv9WNiYhAdHa1zoey1bRuwZ48crGXWLNVpiEgfHTsCNWoAz58nD3lMhsdoCviNGzcwaNAgrFq1CnZ2dnr9TXBwMFxdXbUXLy+vLE5Jr4qJSZ4kYehQoHBhtXmISD9J0/sCwMqVQFiY2jyUOo0QxtFZYPPmzWjVqhUsLS211yUkJECj0cDCwgIxMTE6twFyDzwmJkb7e3R0NLy8vBAVFQUXF5dsy26uZswARo4E8uaVwzQ6O6tOREQZ0aMHsHw5ULUqcOQIYGE0u3zGLTo6Gq6urm+sVUZTwJ88eYJr167pXNe9e3eULFkSI0eORNmyZd+4DX2fFHp3kZFyopKnT+W0oV27qk5ERBnF97Ea+tYqo/k+5ezsjLJly+pcHB0d4ebmplfxpuw1Zox801etCnTurDoNEb0NT09g7Fi5HBAg39NkOIymgJPxOHZMHnYDZLcxHnYjMl7+/kCRIsDt24Ce7YcpmxjNIfTMwEPoWU8IoFYt4I8/5J73ypWqExHRu9qyBWjZErC1Bc6fZ4PUrGZyh9DJOKxbJ4u3gwMwbZrqNESUGT7+GGjQQPYsGT5cdRpKwgJOmebZs+RhGEePBvLnV5uHiDKHRiOn/7W0lGOk79unOhEBLOCUiWbOBG7elOMpJ/X/JiLTULYs0K+fXPb3l9MDk1os4JQprl9PHn5x1izA3l5tHiLKfEFBcibBs2eBb75RnYZYwClTjBgBvHwJ1KkDtGmjOg0RZQU3N2DiRLk8dqycYZDUYQGnd3boELB+vewuNmeOPF9GRKapb1+gTBngwQNgwgTVacwbCzi9k4QEYNAgudy7N1ChgtI4RJTFrKzkF3UA+Oor2a2M1GABp3eyfDlw8iTg6gpMmqQ6DRFlhwYNgE8+kV/ghwyR4z9Q9mMBp7cWFSWHTAWAwEAgd261eYgo+8yaBVhbA7t2ATt2qE5jnljA6a1NngzcvQuUKAH4+alOQ0TZqWhRYPBguTx4MBAbqzaPOWIBp7fy77/J8wWHhgI2NmrzEFH2GzMG8PAALl0C5s9Xncb8sIDTWxk6FIiLA5o1A5o2VZ2GiFRwcUme4GTiRODOHbV5zA0LOGXYrl3Azz/L1qizZ6tOQ0Qqff45ULkyEB2dPPUoZQ8WcMqQuLjk815ffinPfxOR+bKwSD6dtmyZ7JVC2YMFnDJk4ULgwgXA3R0YP151GiIyBDVrAh07yu5kgwaxW1l2YQEnvd2/L7uLAcCUKUCOHErjEJEBmT5dzoFw6BCwYYPqNOaBBZz0Nn488Pgx8N57QM+eqtMQkSHx8gICAuTy8OHA8+dq85gDFnDSy5kzwOLFcnnuXDkvMBHRq4YNAwoWlLMTzpqlOo3pYwGnNxJCzv+bmAi0bStnHCMiep2DAzBjhlyeNg24cUNtHlPHAk5vtHkzsH8/YGsLzJypOg0RGbL27YFatYAXL5IPqVPWYAGndL18KQdtAeR5rUKFlMYhIgOn0cjTbBoNsGYNcPiw6kSmiwWc0hUaCly9CuTPz2/TRKSfSpWSG7oOGiRPv1HmYwGnNEVEyO5igOwi4uioNg8RGY/Jk+VQq8ePA999pzqNaWIBpzSNHg08ewZ88AHw2Weq0xCRMfHwAMaNk8ujRsmhVilzsYBTqv76K/lbc9L5LCKijBg4EChWTE5yMnWq6jSmhwWcUkhMlG88QE5UULWq2jxEZJxsbJInPAoNBS5fVpvH1LCAUwpr1gB//gk4OSVPFUhE9DaaNwcaNwZiY+VAL5R5WMBJx9OnwMiRcnnMGCBvXrV5iMi4aTRy79vSEtiyBdi7V3Ui08ECTjqmT5etz3185OhrRETvqlQpwM9PLvv7A/HxSuOYDBZw0goPTx5pLSQEsLNTGoeITEhQEODmBpw7lzyvAr0bFnDSGj4ciIkB6tUDWrZUnYaITEnOnMCkSXJ53DjgwQO1eUwBCzgBAA4cADZuBCwsgDlz2G2MiDJf795AuXLAo0dyj5zeDQs4ISEh+Xx3377yDUZElNmsrOQOAgAsXAj8/bfSOEaPBZywbBlw+jSQIwcwYYLqNERkyurVA1q1kjsOgwfL6Yrp7bCAm7nHj2V3MUAWb3d3pXGIyAzMmiUHedm7F9i2TXUa48UCbuYmTgTu35fdPPr1U52GiMxB4cLJ0xQPGSIbz1LGsYCbsX/+AebPl8uhoYC1tdo8RGQ+Ro2SA0X995+cb4EyjgXcjA0ZIgdU+OgjOdQhEVF2cXYGpk2Ty5MmAZGRavMYIxZwM7VjB/DLL3KvOyREdRoiMkedO8vJkp4+ldMXU8awgJuh2Fi59w0AgwbJ6f6IiLJb0rgTALBiBXDsmMo0xocF3AwtWABcvAjkyQOMHas6DRGZs+rV5Z64EHKHgt3K9McCbmbu3k3u6z11KuDqqjYPEdG0aYCDA/DHH8C6darTGA8WcDMzbhwQFQVUqgR066Y6DRERkD9/8jnwESOAZ8/U5jEWLOBm5NQpYOlSuTxnjpyfl4jIEAwZAnh7AzdvAjNmqE5jHFjAzYQQcrxzIYAOHYDatVUnIiJKZm8vR2gDZAG/fl1tHmPAAm4mNm0CDh6Uc3zz2y0RGaI2bYA6dYCXL+WhdEofC7gZePECGDZMLo8cCRQsqDYPEVFqNBp5es/CAli/Hjh0SHUiw8YCbgZCQoBr14ACBfitlogMW4UKct5wQHYrS0hQGsegsYCbuFu3gOBguTxjhuyqQURkyCZNkl1cT54Eli9XncZwsYCbuIAA4PlzoGZN4NNPVachInqz3LmBwEC5PGaM7PpKKbGAm7AjR4BVq+R5pblz5U8iImPg5weUKCEHn5o8WXUaw8QCbqISE+X5IwDo3h2oXFltHiKijLCxkdMcA3IH5N9/1eYxRCzgJmrlSuDoUTll35QpqtMQEWVc06ZAs2ZAXBwwdKjqNIaHBdwEPXkCjBoll8eOBTw91eYhInpbs2cDVlbAzz8Du3apTmNYjKaABwcH4/3334ezszPy5MmDli1b4uLFi6pjGaTgYOD2baBIkeTD6ERExqhECeDLL+Xy4MFyb5wkoyngBw8ehJ+fH8LCwrBnzx7Ex8ejUaNGeMZR73VcuSL7fQPym6utrdo8RETvavx4wN0duHABWLhQdRrDoRHCOGdfvXfvHvLkyYODBw/iww8/1OtvoqOj4erqiqioKLi4uGRxQjVatwZ++glo2FAebmLLcyIyBUuWAH36ADlyAJcuyYJuqvStVUazB/66qP93DMyVK1ea68TExCA6OlrnYsr27ZPF29JStt5k8SYiU9GzJ/Dee8Djx3KPnIy0gAshMGTIENSqVQtly5ZNc73g4GC4urpqL15eXtmYMnvFx8vZxgCgXz+gTBmlcYiIMpWlpexOBgCLFwNnzqjNYwiM8hC6n58ftm/fjt9//x0FChRIc72YmBjExMRof4+OjoaXl5dJHkJfuBDo3x/IlUseXkrnwAQRkdFq1w7YuBHw9QV+/dU0jzSa7CH0L7/8Elu3bsX+/fvTLd4AYGtrCxcXF52LKXr4EBg3Ti5PnMjiTUSma+ZM2Th3/35g82bVadQymgIuhMCAAQPw448/Yt++ffDx8VEdyWBMmAA8eCAPm/fpozoNEVHWKVQIGD5cLg8dKucON1dGU8D9/PywatUqrFmzBs7OzoiMjERkZCRevHihOppS588DCxbI5blz5YAHRESmLCAAyJ8fuHo1ebhVc2Q058A1aZzoWL58Obp166bXNkytG5kQQJMmwO7dwCef8HASEZmP1auBzp0BR0c5Tnq+fKoTZR6TOwcuhEj1om/xNkXbt8vibWMDzJqlOg0RUfb57DPggw+AZ8+A0aNVp1HDaAo46YqNlcMKAvJn0aJq8xARZaekaZIB4LvvgL/+UptHBRZwIzVvHnD5spyoZMwY1WmIiLJf1arA55/L5YED5TTK5oQF3AjduSO7iwFy4hJnZ7V5iIhUCQ4GnJyAP/8E1qxRnSZ7sYAbobFj5ZShVaoAXbuqTkNEpE7evMlHIUeOBJ4+VZsnO7GAG5kTJ4Bly+Ty3LmABf+DRGTm/P0BHx8gIgKYPl11muzDj38jIoSc31sI2QKzRg3ViYiI1LOzS55GeeZMIDxcaZxswwJuRH74Afj9d8DBwby+ZRIRvUnLlkC9ekBMTPJIbaaOBdxIPH+e/KIMCADeMAw8EZFZ0WiAOXPkacWNG4EDB1Qnynos4EZi1izgxg2gYEFg2DDVaYiIDE+5ckDfvnLZ3x9ISFAaJ8uxgBuBGzeAadPk8syZgL292jxERIZqwgQgRw7g9OnkBr+migXcCIwcCbx4AdSuLefCJSKi1Lm7yyIOyO5ljx8rjZOlWMAN3O+/A2vXJg8baIqT1xMRZaZ+/YBSpYD795MHvTJFLOAGLDFRnscBgF69gIoVlcYhIjIK1tbJ04zOnw/884/aPFmFBdyAffcdcPw44OICTJ6sOg0RkfFo3Bj46CMgPh4YMkR1mqzBAm6goqOBUaPk8vjxQJ48avMQERmbkBC5N/7LL8COHarTZD4WcAM1ZYqctKRYMeDLL1WnISIyPsWKydErAbkXHhurNk9mYwE3QJcuJZ+/CQ0FbGzU5iEiMlZjx8ojmBcvAgsWqE6TuVjADdCwYUBcHNCkCdCsmeo0RETGy9UVmDpVLk+YANy9qzZPZmIBNzB79gBbtwJWVsDs2ew2RkT0rrp1AypVAqKigHHjVKfJPCzgBiQ+PrnbmJ+f7MdIRETvxtJSjpMOAEuXAqdOqUyTeVjADciiRcD584CbGxAYqDoNEZHpqF0b6NBBTsfs7y9/GjsWcAPx4IHsLgbIPt85c6rNQ0RkambMkHOHHzwIbNqkOs27YwE3EIGBwKNHQPnyQO/eqtMQEZmeggXl3BKAbCz84oXaPO+KBdwA/P23PHwOyPM0lpZK4xARmawRI4ACBYBr1+RAL8aMBVyxpPMxCQlA69aAr6/qREREpsvBQR5KB4DgYODWLbV53gULuGJbtwK//grY2sq5vomIKGt9+ilQsybw/DkQEKA6zdtjAVcoJiZ5kP2hQ4HChdXmISIyB69Oz7xqFXDkiOpEb4cFXKE5c4ArV4C8eZMnLiEioqxXuTLQvbtcHjRITt9sbFjAFYmMTJ4idPp0wMlJbR4iInMzZQrg7AwcPQqsXKk6TcaxgCsyejTw9ClQtSrQqZPqNERE5sfTM3lo1YAA4MkTtXkyigVcgaNHgeXL5fLcuYAF/wtEREoMHAgUKSKPigYHq06TMSwd2UyI5Plpu3QBPvhAbR4iInNmaysnjgJkv/ArV9TmyQgW8Gy2dq1s8ejoaHzf9oiITNFHHwENGwKxsXKENmPBAp6Nnj1LHsZv9Gggf361eYiISHYnCw2Vo2D+9JMcm8MYsIBnoxkzgJs3gUKFkvt/ExGRemXKAP37y2V/fzm9s6FjAc8m164lD983a5acEYeIiAxHUBCQK5ecn2LpUtVp3owFPJuMGAG8fAnUrSvHPCciIsOSKxcwcaJcHjcOePhQbZ43YQHPBr/9Bvzwg+wuNmeOPN9CRESGp08feTj9wQNgwgTVadLHAp7FEhLk+RQA+OIL4L33lMYhIqJ0WFnJ8TkAYMEC4Px5tXnSwwKexZYvB06eBFxdkw/NEBGR4apfH2jZUu6ADR4sx+8wRBku4Dt37sTvv/+u/X3BggWoUKECPvvsMzx69ChTwxm7qCjZXQyQjSNy51Yah4iI9DRrFmBjA+zeDWzfrjpN6jJcwIcPH47o6GgAwNmzZzF06FA0a9YMV65cwRD2jdIxaRJw7x5QsiTg56c6DRER6atIEbn3DcifsbFq86QmwwX86tWrKF26NABg06ZNaNGiBaZOnYqvv/4av/zyS6YHNFb//pt8HiU0FLC2VpuHiIgyZswYOeHJ5cvAvHmq06SU4QJuY2OD58+fAwD27t2LRo0aAQBy5cql3TMnYOhQORBA8+ZAkyaq0xARUUY5OycPeT1xInDnjto8r7PK6B/UqlULQ4YMQc2aNfHXX39h/fr1AIB///0XBQoUyPSAxmjnTuDnn2VrxpAQ1WlMW0JCAuLi4lTHIFLOxsYGFpzaMNN17Spbox87JvfIv/lGdaJkGS7gX331Ffr374+NGzdi4cKFyP//Ab1/+eUXNOGuJuLiks+bDBwIlCihNo+pEkIgMjISjx8/Vh2FyCBYWFjAx8cHNjY2qqOYFAsLeTq0Zk3g22/lcKuVKqlOJWmEMNQG8pkvOjoarq6uiIqKgouLS5bcx9y5st937tzyPHiOHFlyN2bv9u3bePz4MfLkyQMHBwdoODoOmbHExERERETA2toaBQsW5PshC3TqBKxZA9SqJQfnysqnWN9apdceeHR0tHYjbzrPnVWF0Rjcvy+7iwHAlCks3lklISFBW7zd3NxUxyEyCLlz50ZERATi4+NhzVazmW76dGDzZuD33+XImh06qE6kZyO2nDlz4u7duwCAHDlyIGfOnCkuSdebs/HjgcePgQoVgB49VKcxXUnnvB0cHBQnITIcSYfOExISFCcxTQUKAAEBcnn4cOD/bbmV0msPfN++fciVK5d2mYdnUjpzBli8WC7PnSvnlaWsxdchUTK+H7LesGGyEdv168DMmUBgoNo8ehXwOnXqaJfr1q2bVVmMlhDyvHdiItCuHfDhh6oTERFRZrO3l4W7Qwd5SL1HD8DLS12eDPc5GDduXKqHaKKiotCxY8dMCWVsfvoJ2L9fzvGdNOc3Eb2doKAgVKhQ4Z22ER4eDo1Gg1OnTmVKptSsWLECObK4oYtGo8HmzZsBZM9jojdr1w6oXRt48QIYOVJtlgwX8O+//x41a9bEf//9p73uwIEDKFeuHMLDwzMzm1F4+VIO2gLI8yKFCimNQwbs7t276NOnDwoWLAhbW1t4enqicePGOHLkiHadQoUKYc6cOWluY9OmTahWrRpcXV3h7OyMMmXKYGjSC/D/Xrx4gcDAQJQoUQK2trZwd3dH27Ztce7cOZ31goKCoNFooNFoYGlpCS8vL/Tq1Qv37t1L93HcuHEDPXv2RL58+WBjYwNvb28MGjQIDx48yPBz8mqBSjJs2DD8+uuvGd7Wq7y8vHD79m2ULVv2nbZD9DqNRp4m1WiAtWtlozZVMlzAz5w5g0KFCqFChQpYunQphg8fjkaNGqFbt246k5yYi9BQIDwcyJ9f/bcxMmxt2rTB6dOn8d133+Hff//F1q1bUbduXTx8+FCvv9+7dy8+/fRTtG3bFn/99ReOHz+OKVOmIPaVQZpjYmLQoEEDfPvtt5g0aRL+/fdf7NixAwkJCahWrRrCwsJ0tlmmTBncvn0b169fx8KFC7Ft2zZ07do1zQxXrlxBlSpV8O+//2Lt2rW4fPkyFi1ahF9//RXVq1fX+7Gkx8nJ6Z17F1haWsLT0xNWVhke6oLojSpWBHr1ksuDBsnTp0qItzR69Gih0WiEtbW12Lt379tuJsMWLFggChUqJGxtbUWlSpXEb7/9pvffRkVFCQAiKioqU7LcuiWEo6MQgBCrVmXKJkkPL168EOfPnxcvXrxQHUVvjx49EgDEgQMH0l3P29tbhIaGpnrboEGDRN26ddP9+2nTpgmNRiNOnTqlc31CQoKoUqWKKF26tEhMTBRCCBEYGCjee+89nfUmT54sLCwsxPPnz1PdfpMmTUSBAgVS3H779m3h4OAg+vbtq/NYJk6cKDp27CgcHR1F3rx5xbx583RuB6C9eHt7p5rr888/F5988omYMmWKyJMnj3B1dRVBQUEiLi5ODBs2TOTMmVPkz59fLFu2TPs3V69eFQDEyZMntdt49b6SLvv37xdCCBETEyOGDx8u8uXLJxwcHETVqlW1tyVZvny58PLyEvb29qJly5Zi1qxZwtXVNdXnKcmNGzdEhw4dRM6cOYWDg4OoXLmyCAsL096+detWUalSJWFrayt8fHy0jysJAPHTTz+l+pheZ4zvC2N2544QLi7y8/+Vl16m0LdWvdW4e/Pnz0doaCg6duyIwoULY+DAgTh9+vS7f5t4g/Xr18Pf3x9jxozByZMnUbt2bTRt2hTXr1/P8vtOzahRwLNnQPXqwGefKYlAkI0Inz1Tc9F3GCQnJyc4OTlh8+bNiImJeavH6enpiXPnzuHvv/9Oc501a9agYcOGeO+993Sut7CwwODBg3H+/Pl036v29vZITExEfHx8itsePnyIXbt2oX///rC3t0+RrVOnTli/fj3EK0/KzJkzUb58eZw4cQKjRo3C4MGDsWfPHgDA0aNHAQDLly/H7du3tb+nZt++fYiIiMBvv/2G2bNnIygoCC1atEDOnDnx559/om/fvujbty9u3LiR6t/PnTsXt2/f1l4GDRqEPHnyoGTJkgCA7t274/Dhw1i3bh3OnDmDdu3aoUmTJrh06RIA4M8//0SPHj3Qv39/nDp1Cr6+vpg8eXKaeQHg6dOnqFOnDiIiIrB161acPn0aI0aMQOL/d9d27dqFzp07Y+DAgTh//jwWL16MFStWYMqUKelulwxDnjyy6zAgp41WMhVIRr8ZNGnSRLi5uYkNGzYIIYR4/vy56Nu3r7CzsxPTp09/q28b+qpatarON3whhChZsqQICAjQ6+8zcw/84EH5zQsQ4q+/3nlzlAGv72k8fZr8v8juy9On+ufeuHGjyJkzp7CzsxM1atQQo0aNEqdPn9ZZJ7098KdPn4pmzZpp91Y7dOggli1bJl6+fKldx87OTgwaNCjVvz9x4oQAINavXy+ESLmne+HCBVG0aFFRtWrVVP8+LCxMZ4/wdbNnzxYAxJ07d7SPpUmTJjrrdOjQQTRt2lT7e2rbS20P3NvbWyQkJGivK1GihKhdu7b29/j4eOHo6CjWrl0rhEh/b3XTpk3C1tZWHDp0SAghxOXLl4VGoxG3bt3SWa9+/fpi1KhRQgghOnbsmOpjSW8PfPHixcLZ2Vk8ePAg1dtr164tpk6dqnPdypUrRd68ebW/v/r8cA/c8MTECFGsmPwsGDEi87abZXvg8fHxOHPmDNq2bQtAfmNfuHAhNm7ciNDQ0Ez6WpFSbGwsjh8/rp39LEmjRo3wxx9/pPo3MTExiI6O1rlkhvv3gaSedR98ALz/fqZslkxcmzZttHtjjRs3xoEDB1CpUiWsWLFCr793dHTE9u3bcfnyZYwdOxZOTk4YOnQoqlatqp0hMD3i/3vGr46VffbsWTg5OcHe3h6lS5eGl5cXVq9e/VaPL2n7r/ZHrl69us461atXx4ULFzK87TJlyuhM1OHh4YFy5cppf7e0tISbm5t2wKm0nDx5El27dsWCBQtQq1YtAMCJEycghEDx4sW1R0qcnJxw8OBBbWPdCxcupPpY0nPq1ClUrFhRO4bG644fP46JEyfq3Gfv3r1x+/Ztvf6fpJ6NDTB7tlwODQXe8PLLdBlu4ZF0+Ot1zZs3x9mzZ985UFru37+PhIQEeHh46Fzv4eGByMjIVP8mODgYEyZMyPQsrw7S8vJlpm+eMsjBAXj6VN19Z4SdnR0aNmyIhg0bYvz48ejVqxcCAwPRrVs3vbdRpEgRFClSBL169cKYMWNQvHhxrF+/Ht27d0exYsVw/vz5VP/un3/+AQAUL15ce12JEiWwdetWWFpaIl++fLC1tU3zfosWLQqNRoPz58+jZcuWqW4/Z86ccHd3Tzf/2ww48vrQoBqNJtXrEtNpTRQZGYmPP/4YPXv2RM+ePbXXJyYmwtLSEsePH4flayMwOTk5AUj+cpIRr59meF1iYiImTJiA1q1bp7jNzs4uw/dHaiQVbXt72ZU4O2Xq3HNveuNmhtff/EKIND8QRo0ahaioKO0lrfNjGZUzp5yoBABOnQJ27MiUzdJb0mgAR0c1l3cd/Kp06dJ49uzZW/99oUKF4ODgoN1Gx44dsXfv3hTnuRMTExEaGooqVaqgdOnS2uttbGxQtGhR+Pj4pFu8AcDNzQ0NGzbE119/jRcvXujcFhkZidWrV6NDhw4678fXW72HhYVpzzsDsjBnx9CfL1++xCeffIKSJUtidtIu0/9VrFgRCQkJuHv3LooWLapz8fT0BCD/T6k9lvSUL18ep06dSrNlfqVKlXDx4sUU91m0aFFOC2okoqNlWyhAjsqW3VOBZHgPPCEhAaGhofjhhx9w/fp1nS4sADKlG0lq3N3dYWlpmWJv++7duyn2ypPY2tq+8UPpbRUrJvt9z5wppw9t0EAeTiFKzYMHD9CuXTv06NED5cuXh7OzM44dO4YZM2bgk08+0Vn31q1bKQbrKFiwIObNm4fnz5+jWbNm8Pb2xuPHjzFv3jzExcWhYcOGAIDBgwdjy5Yt+OijjxASEoJq1arhzp07mDp1Ki5duoTDhw+/0+P46quvUKNGDTRu3BiTJ0+Gj48Pzp07h+HDhyN//vwpGmAdPnwYM2bMQMuWLbFnzx5s2LAB27dv195eqFAh/Prrr6hZsyZsbW2zbD6FPn364MaNG/j11191+rnnypULxYsXR6dOndC1a1eEhISgYsWKuH//Pvbt24dy5cqhWbNmGDhwIGrUqKF9LLt378bOnTvTvc+OHTti6tSpaNmyJYKDg5E3b16cPHkS+fLlQ/Xq1TF+/Hi0aNECXl5eaNeuHSwsLHDmzBmcPXv2jQ3kyDBMniz3wIsXBwYMUBAgoyfXx40bJ/LmzStmzpwp7OzsxKRJk0TPnj2Fm5ubmDt37tucr9db1apVRb9+/XSuK1WqlJJGbHJ7QuTJIxswhIRkyiZJD8bYWOfly5ciICBAVKpUSbi6ugoHBwdRokQJMXbsWJ0uWa93rUq6LF++XOzbt0+0adNGeHl5CRsbG+Hh4SGaNGmibYyV5OnTp2LMmDGiSJEiwsrKSgAQRYsWFTdu3NBZL7VuZPoIDw8X3bp1E56ensLa2lp4eXmJL7/8Uty/f19nPW9vbzFhwgTRvn174eDgIDw8PMScOXN01tm6dasoWrSosLKyemM3slfVqVMnRWO9VxsAvt7gK63nNamrWGxsrBg/frwoVKiQsLa2Fp6enqJVq1bizJkz2u0vW7ZMFChQQNjb24uPPvpIr25k4eHhok2bNsLFxUU4ODiIKlWqiD///FN7+86dO0WNGjWEvb29cHFxEVWrVhVLlizR3g42YjNY//4rhLW1/Pz/+efM3ba+tSrDBbxw4cLi5/+ndXJyEpcvXxZCCDF37lzRsWPHt4iqv3Xr1glra2uxbNkycf78eeHv7y8cHR1FeHi4Xn+f2QVcCNn/DxDC1VX2C6Ssxw+qjNmxY4ewtbUV8+fPz9b7Ta9FPWU+vi+y18cfy8/+Jk2E+P/QCpkmy1qhR0ZGalt/Ojk5ISoqCgDQokULnUNjWaFDhw6YM2cOJk6ciAoVKuC3337Djh074O3tnaX3m55u3YBKlYCoKGDcOGUxiNLUtGlT/PLLL3j48CHu37+vOg6R0du9G9i6FbCykq3QVU0El+ECXqBAAdy+fRuAbJW6e/duAHJQhqw63/yq/v37Izw8HDExMTh+/Dg+VDz1l4WFHBcXAJYulY3aiAyNr68vxo8fny0NTYlMWVycbPcEyPPepUqpy5LhAt6qVSvtRAODBg3CuHHjUKxYMXTt2hU9evTI9IDGoFYt4NNP5dAegwbpPzoXkSkLDw+Hv7+/6hhEmWrRIuD8ecDdPXkkNlU0QrxbuQkLC8Mff/yBokWL4uOPP86sXFkiOjoarq6uiIqKgksmt/e/fh0oWVJOMffDD3LKOcoaL1++xNWrV+Hj48P+skT/x/dF1nvwQPZAevQIWLgQ6Ns3a+5H31r1zp0NP/jgAwwZMsTgi3dWK1gweTay4cNlISciItMRGCiLd/nyQO/eqtO8YwF3cXHBlStXMiuL0Rs+HPDyAq5dA0JCVKcxfemNukVkbt7xYCq9wdmzcq8bAObM0R2RUxW9B3K5efMmChQooHMdXzC6HByAGTOAjh2B4GDZQv21p4wygY2NDSwsLBAREYHcuXPDxsbmrYbnJDIVQgjcu3cv1SFm6d0JAfj7y3m/27QBfH1VJ5L0LuBly5bF/Pnz0aVLl6zMY/Q6dAC++go4fBgICABWrVKdyPRYWFjAx8cHt2/fRkREhOo4RAZBo9GgQIECKcZzp3e3ZQuwbx9gaytH3zQUehfwqVOnws/PD5s3b8aSJUvg5uaGzp07Z3pjMGOn0chuZe+/D6xeDfj5yfnCKXPZ2NigYMGCiI+Pz5axtIkMnbW1NYt3FoiJAYYOlcvDhgE+PmrzvCpDrdCvXr2Knj174vz581iyZInRNVzLylbor+vVC1i2TBbysDDZX5yIiIzL9OnyaGrevHISq/9PUJel9K1Vb9WN7KuvvsLgwYNRqlQpWFnp7sSfOHEi42mzSXYW8Dt3ZHeDJ0+AFSuAzz/P0rsjIqJMdvu2nKjk6VPg+++B7DqDrG+tyvBsZNeuXcOmTZuQK1cufPLJJykKOEkeHnJo1REj5Le31q0BZ2fVqYiISF+jR8viXa0a0KmT6jQpZaj6Ll26FEOHDkWDBg3w999/I3fu3FmVyyQMHAgsWQJcvgxMnSpbphMRkeE7elQePQVkuyZDPA2qd6QmTZpg5MiR+Oqrr/Djjz+yeOvB1lYOdA/In//9pzYPERG9WdKw2ADQtavcAzdEehfwhIQEnDlzBl27ds3KPCanRQugUSMgNlYO9EJERIZt7VrgyBHA0dGwj5zqXcD37NmTYiAXejONBggNlaP2/PQT8P95YIiIyAA9eybbLgHyHHi+fGrzpMcAj+qbntKlgf795bK/PxAfrzQOERGlYfp04NYtoFAhYMgQ1WnSxwKeTYKCgFy5gL//lg3biIjIsISHJ4+0FhICGPqkbizg2SRXLmDSJLk8bhzw8KHaPEREpGvECODlSznWeatWqtO8GQt4NvriC6BsWVm8J0xQnYaIiJL89huwYYPsLjZnjmy/ZOhYwLORlZV8YQDAggXA+fNK4xAREYCEhORuY198Ief7NgYs4Nmsfn2gZUv5gvH3l/0NiYhInW+/BU6dAnLkACZOVJ1GfyzgCsyaBdjYAHv2AD//rDoNEZH5evwYGDNGLgcFAcY0RhkLuAJFiiR3TxgyRE5XR0RE2W/SJODePaBkyeTuvsaCBVyR0aMBT085Tvq8earTEBGZn4sXkz9/Q0MBa2u1eTKKBVwRZ+fkIfomTZLTjxIRUfYZOlQOrNW8OdCkieo0GccCrlDXrkCVKnLO8KRzMERElPV++QXYvl32DkqadMrYsIArZGGRfPjm22+B48fV5iEiMgdxccDgwXJ50CCgeHG1ed4WC7hi1avLieKTpq9jtzIioqy1YIE8/507txwZ01ixgBuAadMABwfg8GFg/XrVaYiITNe9e7K7GABMmQK4uiqN805YwA1AgQJAQIBcHjECeP5cbR4iIlM1fjwQFQVUqAD06KE6zbthATcQw4YBBQsCN24kz4ZDRESZ5/Tp5Nkg584FLC3V5nlXLOAGwt5ejtAGyPlor19Xm4eIyJQIIYevTkwE2rcHPvxQdaJ3xwJuQNq2lS+qFy+AkSNVpyEiMh0//ggcOCDn+J4xQ3WazMECbkA0GnlYR6MB1q0Dfv9ddSIiIuP38qU8TQkAw4cD3t5q82QWFnADU6EC0KuXXB40SB7uISKitzd7NhAeDuTPb1pHN1nADdDkyYCLC3DiBLBiheo0RETG69YtYOpUuTxjBuDoqDZPZmIBN0B58gCBgXJ51CggOlptHiIiYzVqFPDsGVCjBtCxo+o0mYsF3EANGCCH97t7V+6RExFRxoSFAStXyuWk9kWmhAXcQNnYyOntAGDOHODSJaVxiIiMSmKi7DYGAN26yYmjTA0LuAFr1kxOcRcXl9yCkoiI3mz1auDPPwEnp+Rz4KaGBdzAzZ4tp7vbuhXYvVt1GiIiw/f0aXJr87Fjgbx51ebJKizgBq5UKXk+HJDT38XFqc1DRGTogoOB27eBIkWSD6ObIhZwIzB+PODuDpw/DyxapDoNEZHhunoVCAmRyyEhgK2t2jxZiQXcCOTMmdwSPTAQePBAbR4iIkM1fDgQEwPUrw98/LHqNFmLBdxI9OoFlC8PPHqU3EeciIiS7d8PbNoEWFjI3jum1m3sdSzgRsLSUr4gAWDhQuDsWaVxiIgMSnx88vnufv2AsmWVxskWLOBGxNcXaNMmuX+jEKoTEREZhm++Ac6ckaccJ0xQnSZ7sIAbmZkzZaOMffuALVtUpyEiUu/RI9ldDAAmTgTc3NTmyS4s4EbGxyd5UJehQ2VjDSIiczZxomzcW7o00Lev6jTZhwXcCAUEAPnyAVeuJJ8XJyIyRxcuAF99JZfnzJEDX5kLFnAj5OQETJsmlydPlgMWEBGZGyHkAFfx8bLLWMOGqhNlLxZwI9WpE1CtmhwycPRo1WmIiLLfjh3Arl2AtXXy4C3mhAXcSFlYyOnxAGDFCuDoUaVxiIiyVWys3PsG5M+iRdXmUYEF3IhVqwZ07SqXBw1itzIiMh/z58tplj08gDFjVKdRwygKeHh4OHr27AkfHx/Y29ujSJEiCAwMRGxsrOpoygUHA46OwJEjwNq1qtMQEWW9u3dly3NAThXq4qI2jypGUcD/+ecfJCYmYvHixTh37hxCQ0OxaNEijObJX+TLl3wOfMQI4NkztXmIiLLa2LFAdDRQuTLQrZvqNOpohDDOA68zZ87EwoULceXKFb3/Jjo6Gq6uroiKioKLCX1le/lSTjsaHg6MG5f8zZSIyNScPCkLtxDA778DNWuqTpT59K1VRrEHnpqoqCjkypVLdQyDYGeX3AJz5kxZyImITI0Qye19OnY0zeKdEUZZwP/77z/Mnz8ffd8w5E5MTAyio6N1LqaqVSs5VvrLl/JQOhGRqdmwATh0CLC3B6ZPV51GPaUFPCgoCBqNJt3LsWPHdP4mIiICTZo0Qbt27dCrV690tx8cHAxXV1ftxcvLKysfjlIajRyFyMJCvsh/+011IiKizPPihZzrGwBGjgRM+ONcb0rPgd+/fx/3799Pd51ChQrBzs4OgCzevr6+qFatGlasWAELi/S/f8TExCDmlcHCo6Oj4eXlZXLnwF/Vrx+waBFQoQJw7JichpSIyNhNmgSMHy8L9z//AA4OqhNlHX3PgRtNI7Zbt27B19cXlStXxqpVq2D5FpXJVBuxverePaB4ceDxY2DJEqB3b9WJiIjezY0bQIkSci983TqgQwfVibKWSTVii4iIQN26deHl5YVZs2bh3r17iIyMRGRkpOpoBid3biAoSC6PGSMLORGRMQsIkMW7Vi2gfXvVaQyHURTw3bt34/Lly9i3bx8KFCiAvHnzai+UUv/+QMmScm980iTVaYiI3t4ffwBr1sh2PnPnyp8kGUUB79atG4QQqV4oJWvr5GlG580DLl5UGoeI6K0kJspuYwDQowdQqZLaPIbGKAo4ZVzjxkDz5nKavaFDVachIsq477+XjXGdnYEpU1SnMTws4CZs9mw5uf327cAvv6hOQ0Skv+hoee4bkK3PPTzU5jFELOAmrHjx5MNPgwcDcXFq8xAR6WvqVODOHaBYMWDgQNVpDBMLuIkbN062TL94EViwQHUaIqI3++8/IDRULs+eDdjYqM1jqFjATZyrq/wmC8juZffuKY1DRPRGw4YBsbFAo0ayLQ+ljgXcDHTvLkdmi4qS55KIiAzV3r3A5s1yFMnQUHYbSw8LuBmwtJT9JwE5Otvp02rzEBGlJj4e8PeXy35+QOnSSuMYPBZwM/Hhh3IEo8RE+QZhF3oiMjSLFwPnzgFubskjSlLaWMDNyIwZcu7wAweAH39UnYaIKNnDh8mn+CZNAnLmVJvHGLCAmxFv7+S5wocNk3OHExEZgqAgWcTLluUkTPpiATczI0YA+fMD4eGyewYRkWrnzgFffy2X58yRA1DRm7GAmxlHR3koHZDdy27dUpuHiMybEHKgqYQEoFUroH591YmMBwu4GerYEahRA3j2DBg1SnUaIjJn27YBe/bIwVpmzVKdxriwgJuhpGn5AGDlSiAsTG0eIjJPMTHAkCFyeehQoHBhtXmMDQu4mapSRQ7wAsjx0hMT1eYhIvMzb54cNtXTk0cD3wYLuBmbOhVwcgL++gtYvVp1GiIyJ5GRsrsYAEybJqcMpYxhATdjnp7A2LFyeeRI4OlTtXmIyHyMGQM8eQK8/z7QpYvqNMaJBdzM+fsDRYoAt28DwcGq0xCROTh2DFi+XC7PnQtYsBK9FT5tZs7WFggJkcshIcCVK2rzEJFpEyJ5OOfOnYHq1VUnMl4s4ISPPwYaNJAtQocPV52GiEzZ+vXA4cOAg4M8901vjwWcoNHIafssLOQY6fv3q05ERKbo+fPknYRRo+SokPT2WMAJgBx/uF8/uezvL6f1IyLKTDNmADdvynkZhg5Vncb4sYCT1oQJcgagM2eAb75RnYaITMn168D06XJ51izA3l5tHlPAAk5abm7AxIlyeexY4NEjtXmIyHSMGCFnQKxTB2jTRnUa08ACTjr69gXKlAEePJB75ERE7+rQIdl4TaORs41pNKoTmQYWcNJhZSUbtAHAggXAhQtq8xCRcUtIkMM1A3Ke7woVlMYxKSzglELDhrJrWXy8nOZPCNWJiMhYrVgBnDwJuLoCkyerTmNaWMApVSEhgLU1sGsXsGOH6jREZIyiooDRo+VyYCCQO7faPKaGBZxSVbSo3PsG5M/YWLV5iMj4TJ4M3L0LlCgB+PmpTmN6WMApTWPGAB4ewKVLwPz5qtMQkTG5dEmOcw4As2cDNjZq85giFnBKk4uLnHIUkN3L7t5Vm4eIjMfQoUBcHNC0KdCsmeo0pokFnNLVrRtQuTIQHZ089SgRUXp27QK2bZO9WmbPVp3GdLGAU7osLJIPg33zjWxNSkSUlri45PYzX34JlCypNo8pYwGnN6pZE+jYUXYnGzSI3cqIKG0LF8rxI9zdgfHjVacxbSzgpJfp0+XYxYcOARs2qE5DRIbo/n3ZXQyQLdBz5FAax+SxgJNevLyAkSPl8vDhwIsXavMQkeEJDAQePwbKlwd69VKdxvSxgJPehg+Xhfz6dTmbEBFRkrNngUWL5PLcuYClpdo85oAFnPTm4ADMnCmXg4OBGzfU5iEiw5DUPiYxEWjbFqhbV3Ui88ACThnSvj1Qq5Y8hB4QoDoNERmCzZuB/fsBW9vkL/mU9VjAKUM0Gnl4TKMB1qwB/vhDdSIiUunlSzloCwAMGwYUKqQ0jllhAacMq1QJ6NFDLicdNiMi8zRnDnD1KpAvH4/KZTcWcHorU6YAzs7AsWPA99+rTkNEKkREJE8ROn064OSkNo+5YQGnt+LhkTxIQ0CAHGqViMzL6NHAs2fABx8An32mOo35YQGntzZwIFCsGHDnTvKkJ0RkHv76C/juO7k8d64cdpmyF59yems2NskTFYSGAv/9pzYPEWWPpG5jANC1K1C1qto85ooFnN5J8+ZAo0ZAbKxsgUpEpm/NGiAsDHB0lGNCkBos4PRONBq5921pKfuC7t2rOhERZaWnT4ERI+TymDGy9TmpwQJO76x0acDPTy77+wPx8UrjEFEWmj5dtj738UmeNpTUYAGnTBEUBLi5AefOAYsXq05DRFkhPDx5pLWQEMDOTmkcs8cCTpkiZ05g0iS5PH488PCh2jxElPmGDwdiYgBfX6BlS9VpiAWcMk3v3kDZsrJ4BwWpTkNEmengQWDjRtldbM4c2f6F1GIBp0xjZSXf2ADw9dfycDoRGb+EhORuY336yPm+ST0WcMpU9esDrVrJN/zgwbK/KBEZt2XLgNOngRw5gIkTVaehJCzglOlmzZKDvOzZA2zbpjoNEb2Lx49ldzEAmDABcHdXGodeYXQFPCYmBhUqVIBGo8GpU6dUx6FUFC6cPL3gkCGy0QsRGaeJE4H794FSpYB+/VSnoVcZXQEfMWIE8nHkAIM3ahTg6SmHV503T3UaInobFy8C8+fL5dBQwNpabR7SZVQF/JdffsHu3bsxa9Ys1VHoDZydgWnT5PKkSUBkpNo8RJRxQ4bIgZlatAAaN1adhl5nNAX8zp076N27N1auXAkHBwe9/iYmJgbR0dE6F8o+XboA778PPHmSfA6NiIzDjh3yYm0tB20hw2MUBVwIgW7duqFv376oUqWK3n8XHBwMV1dX7cXLyysLU9LrLCzkNIMAsHw5cOyY2jxEpJ/YWLn3DcjuY8WLq81DqVNawIOCgqDRaNK9HDt2DPPnz0d0dDRGjRqVoe2PGjUKUVFR2suNGzey6JFQWqpXBzp3lt3J/P3ZrYzIGCxYIM9/58kDjB2rOg2lRSOEuo/U+/fv4/79++muU6hQIXz66afYtm0bNK8M/ZOQkABLS0t06tQJ3yXNKv8G0dHRcHV1RVRUFFxcXN4pO+nv1i35Df75c2DtWuDTT1UnIqK03LsHFCsGREUBS5cCvXqpTmR+9K1VSgu4vq5fv65z/joiIgKNGzfGxo0bUa1aNRQoUECv7bCAqzN5MjBuHFCggPxmr2czBiLKZn37ygmJKlYEjh6VUwVT9tK3VlllY6a3VrBgQZ3fnZycAABFihTRu3iTWkOHAt98A1y7BsyYwbHSiQzRqVPAkiVyee5cFm9DZxSN2Mj42dvLEdoAOZ/w9etq8xCRrlfbqXToANSurToRvYlRFvBChQpBCIEKFSqojkIZ0KYNUKcO8PIlMGKE6jRE9KpNm+SMY3Z28igZGT6jLOBknDSa5GkI168HDh1SnYiIAODFCznXNyC/XL921pIMFAs4ZasKFeS84YDsX5qQoDQOEQGYPRsID5eNTHl0zHiwgFO2mzwZcHUFTp4EVqxQnYbIvN26BUydKpdnzAAcHdXmIf2xgFO2y50bCAyUy6NHy/6mRKRGQIAco6FmTY7RYGxYwEkJPz+gRAng7l25R05E2S8sDFi1SrZLmTtX/iTjwQJOStjYyPNugPzguHRJbR4ic5OYKNuhAEC3bkDlykrj0FtgASdlmjUDmjYF4uLkQC9ElH1WrQL++gtwcko+B07GhQWclJo9G7CyArZtA3btUp2GyDw8eSLPfQNyiGNPT7V56O2wgJNSJUsCX34plwcPlnvjRJS1goOB27eBIkWSD6OT8WEBJ+XGjwfc3YELF4CFC1WnITJtV64AISFyefZswNZWbR56eyzgpFyOHMkt0QMDgTfMMEtE72DYMCA2FmjQAPjoI9Vp6F2wgJNB6NULKF8eePw4uY84EWWuffuAn36Ss4yFhrLbmLFjASeDYGkpu5MBwKJFwNmzavMQmZr4eDnbGAD06weULas0DmUCFnAyGHXrAm3bJvdPFUJ1IiLTsXSp/GKcKxcwYYLqNJQZWMDJoMycKRvV7N8PbN6sOg2RaXj4UHYXA4CJE2URJ+PHAk4GpVAh2cgGkIO7vHypNA6RSZgwAXjwAChTBujTR3Uayiws4GRwAgKAfPmAq1fl/OFE9PbOnwcWLJDLc+bIgZPINLCAk8FxcgKmT5fLkycDERFq8xAZKyHkAEkJCcAnn8iuY2Q6WMDJIH32GfDBB8CzZ3LKUSLKuO3bgd275eRBs2apTkOZjQWcDJKFRXK3su++k5MuEJH+YmPl3jcgfxYtqjYPZT4WcDJYVasCXbvKZXYrI8qYefOAy5cBDw9gzBjVaSgrsICTQQsOBhwdgbAwYM0a1WmIjMOdO8CkSXI5OBhwdlabh7IGCzgZtHz5kvceRowAnj5Vm4fIGIwdC0RHA5UrA59/rjoNZRUWcDJ4gwcDPj6yNXpS63QiSt2JE8CyZXJ53jzZnoRME/+1ZPDs7JKnP5w5EwgPVxqHyGAJkdxe5LPPgBo1VCeirMQCTkahZUvA1xeIiQGGD1edhsgw/fAD8PvvgL09MG2a6jSU1VjAyShoNHIUKQsLYONG4OBB1YmIDMvz57KdCCBHM/TyUpuHsh4LOBmN8uWTx3EeNEiOLkVE0qxZwPXrQMGCyfMJkGljASejMnEikCMHcPp0ckMdInN340byIfOZMwEHB7V5KHuwgJNRcXdPnst4zBjg8WOlcYgMwsiRwIsXQO3aQLt2qtNQdtEIYT7jW0VHR8PV1RVRUVFwcXFRHYfeUlwc8N57wIULgJ8fG7WRefv7b6BFC9lO5PhxoGJF1YnoXelbq1jAySjt2gU0aaI6BZHh6NULWLpUdQrKDPrWKs4MS0apcWP5gbV6NcdIJ/LxAaZMUZ2Cshv3wImIiAyIvrWKjdiIiIiMEAs4ERGREWIBJyIiMkIs4EREREaIBZyIiMgIsYATEREZIRZwIiIiI8QCTkREZIRYwImIiIwQCzgREZERYgEnIiIyQizgRERERogFnIiIyAixgBMRERkhs5oPPGnm1OjoaMVJiIiIUpdUo94027dZFfAnT54AALy8vBQnISIiSt+TJ0/g6uqa5u0a8aYSb0ISExMREREBZ2dnaDSad9pWdHQ0vLy8cOPGjXQnXDdlfA74HJj74wf4HJj74wcy/zkQQuDJkyfIly8fLCzSPtNtVnvgFhYWKFCgQKZu08XFxWxftEn4HPA5MPfHD/A5MPfHD2Tuc5DenncSNmIjIiIyQizgRERERogF/C3Z2toiMDAQtra2qqMow+eAz4G5P36Az4G5P35A3XNgVo3YiIiITAX3wImIiIwQCzgREZERYgEnIiIyQizgRERERogFPJNs374d1apVg729Pdzd3dG6dWvVkZSIiYlBhQoVoNFocOrUKdVxsk14eDh69uwJHx8f2Nvbo0iRIggMDERsbKzqaFnq66+/ho+PD+zs7FC5cmUcOnRIdaRsERwcjPfffx/Ozs7IkycPWrZsiYsXL6qOpVRwcDA0Gg38/f1VR8lWt27dQufOneHm5gYHBwdUqFABx48fz5b7ZgHPBJs2bUKXLl3QvXt3nD59GocPH8Znn32mOpYSI0aMQL58+VTHyHb//PMPEhMTsXjxYpw7dw6hoaFYtGgRRo8erTpallm/fj38/f0xZswYnDx5ErVr10bTpk1x/fp11dGy3MGDB+Hn54ewsDDs2bMH8fHxaNSoEZ49e6Y6mhJHjx7FkiVLUL58edVRstWjR49Qs2ZNWFtb45dffsH58+cREhKCHDlyZE8AQe8kLi5O5M+fX3zzzTeqoyi3Y8cOUbJkSXHu3DkBQJw8eVJ1JKVmzJghfHx8VMfIMlWrVhV9+/bVua5kyZIiICBAUSJ17t69KwCIgwcPqo6S7Z48eSKKFSsm9uzZI+rUqSMGDRqkOlK2GTlypKhVq5ay++ce+Ds6ceIEbt26BQsLC1SsWBF58+ZF06ZNce7cOdXRstWdO3fQu3dvrFy5Eg4ODqrjGISoqCjkypVLdYwsERsbi+PHj6NRo0Y61zdq1Ah//PGHolTqREVFAYDJ/r/T4+fnh+bNm6NBgwaqo2S7rVu3okqVKmjXrh3y5MmDihUrYunSpdl2/yzg7+jKlSsAgKCgIIwdOxY///wzcubMiTp16uDhw4eK02UPIQS6deuGvn37okqVKqrjGIT//vsP8+fPR9++fVVHyRL3799HQkICPDw8dK738PBAZGSkolRqCCEwZMgQ1KpVC2XLllUdJ1utW7cOJ06cQHBwsOooSly5cgULFy5EsWLFsGvXLvTt2xcDBw7E999/ny33zwKehqCgIGg0mnQvx44dQ2JiIgBgzJgxaNOmDSpXrozly5dDo9Fgw4YNih/Fu9H3OZg/fz6io6MxatQo1ZEznb7PwasiIiLQpEkTtGvXDr169VKUPHu8Pi2vEOKdp+o1NgMGDMCZM2ewdu1a1VGy1Y0bNzBo0CCsWrUKdnZ2quMokZiYiEqVKmHq1KmoWLEi+vTpg969e2PhwoXZcv9mNZ1oRgwYMACffvppuusUKlQIT548AQCULl1ae72trS0KFy5s9I159H0OJk+ejLCwsBTjAFepUgWdOnXCd999l5Uxs5S+z0GSiIgI+Pr6onr16liyZEkWp1PH3d0dlpaWKfa27969m2Kv3JR9+eWX2Lp1K3777bdMn6rY0B0/fhx3795F5cqVtdclJCTgt99+w1dffYWYmBhYWloqTJj18ubNq/PZDwClSpXCpk2bsuX+WcDT4O7uDnd39zeuV7lyZdja2uLixYuoVasWACAuLg7h4eHw9vbO6phZSt/nYN68eZg8ebL294iICDRu3Bjr169HtWrVsjJiltP3OQBkdxJfX1/tURgLC9M9wGVjY4PKlStjz549aNWqlfb6PXv24JNPPlGYLHsIIfDll1/ip59+woEDB+Dj46M6UrarX78+zp49q3Nd9+7dUbJkSYwcOdLkizcA1KxZM0X3wX///TfbPvtZwN+Ri4sL+vbti8DAQHh5ecHb2xszZ84EALRr105xuuxRsGBBnd+dnJwAAEWKFDGbvZKIiAjUrVsXBQsWxKxZs3Dv3j3tbZ6engqTZZ0hQ4agS5cuqFKlivaIw/Xr1032vP+r/Pz8sGbNGmzZsgXOzs7aIxGurq6wt7dXnC57ODs7pzjn7+joCDc3N7NpCzB48GDUqFEDU6dORfv27fHXX39hyZIl2Xf0TVn7dxMSGxsrhg4dKvLkySOcnZ1FgwYNxN9//606ljJXr141u25ky5cvFwBSvZiyBQsWCG9vb2FjYyMqVapkNt2o0vpfL1++XHU0pcytG5kQQmzbtk2ULVtW2NraipIlS4olS5Zk231zOlEiIiIjZLon6YiIiEwYCzgREZERYgEnIiIyQizgRERERogFnIiIyAixgBMRERkhFnAiIiIjxAJORNniwIED0Gg0ePz4seooRCaBBZzIzCQkJKBGjRpo06aNzvVRUVHw8vLC2LFjs+R+a9Sogdu3b8PV1TVLtk9kbjgSG5EZunTpEipUqIAlS5agU6dOAICuXbvi9OnTOHr0KGxsbBQnJKI34R44kRkqVqwYgoOD8eWXXyIiIgJbtmzBunXr8N1336VZvEeOHInixYvDwcEBhQsXxrhx4xAXFwdAzs7VoEEDNGnSBEn7BI8fP0bBggUxZswYACkPoV+7dg0fffQRcubMCUdHR5QpUwY7duzI+gdPZCI4GxmRmUqaDrNr1644e/Ysxo8fjwoVKqS5vrOzM1asWIF8+fLh7Nmz6N27N5ydnTFixAhoNBp89913KFeuHObNm4dBgwahb9++8PDwQFBQUKrb8/PzQ2xsLH777Tc4Ojri/Pnz2pnsiOjNeAidyIz9888/KFWqFMqVK4cTJ07Aykr/7/QzZ87E+vXrcezYMe11GzZsQJcuXTBkyBDMnTsXJ0+eRPHixQHIPXBfX188evQIOXLkQPny5dGmTRsEBgZm+uMiMgc8hE5kxr799ls4ODjg6tWruHnzJgCgb9++cHJy0l6SbNy4EbVq1YKnpyecnJwwbtw4XL9+XWd77dq1Q+vWrREcHIyQkBBt8U7NwIEDMXnyZNSsWROBgYE4c+ZM1jxIIhPFAk5kpo4cOYLQ0FBs2bIF1atXR8+ePSGEwMSJE3Hq1CntBQDCwsLw6aefomnTpvj5559x8uRJjBkzBrGxsTrbfP78OY4fPw5LS0tcunQp3fvv1asXrly5gi5duuDs2bOoUqUK5s+fn1UPl8j0ZNvM40RkMJ4/fy6KFSsm/Pz8hBBCXLt2TTg7O4uFCxemuv6sWbNE4cKFda7r2bOncHV11bmub9++omTJkmL37t3CyspK/Prrr9rb9u/fLwCIR48epXofAQEBoly5cm//oIjMDPfAicxQQEAAEhMTMX36dABAwYIFERISguHDhyM8PDzF+kWLFsX169exbt06/Pfff5g3bx5++uknnXW2b9+Ob7/9FqtXr0bDhg0REBCAzz//HI8ePUo1g7+/P3bt2oWrV6/ixIkT2LdvH0qVKpXpj5XIZKn+BkFE2evAgQPC0tJSHDp0KMVtjRo1EvXq1ROJiYkpbhs+fLhwc3MTTk5OokOHDiI0NFS7B3737l3h4eEhpk6dql0/Li5OVK1aVbRv314IkXIPfMCAAaJIkSLC1tZW5M6dW3Tp0kXcv38/8x8wkYliK3QiIiIjxEPoRERERogFnIiIyAixgBMRERkhFnAiIiIjxAJORERkhFjAiYiIjBALOBERkRFiASciIjJCLOBERERGiAWciIjICLGAExERGSEWcCIiIiP0P8abFj/MIYAaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "from honeycomb_analysis import youngs_modulus_E1, youngs_modulus_E2, shear_modulus_G12\n",
    "\n",
    "# Fixed length of the cell (height of the hexagon)\n",
    "l = 6.0\n",
    "\n",
    "# Define weight coefficients\n",
    "lambda_1 = 0.33\n",
    "lambda_2 = 0.33\n",
    "lambda_3 = 0.34\n",
    "\n",
    "# Ensure the sum of weights is 1\n",
    "assert lambda_1 + lambda_2 + lambda_3 == 1.0, \"The weights should sum to 1.\"\n",
    "\n",
    "def objective(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    try:\n",
    "        E1 = youngs_modulus_E1(beta_l, l, alpha_l, theta_deg)\n",
    "        E2 = youngs_modulus_E2(beta_l, l, alpha_l, theta_deg)\n",
    "        G12 = shear_modulus_G12(beta_l, l, alpha_l, theta_deg)\n",
    "        \n",
    "        if np.any(np.isnan([E1, E2, G12])) or np.any(np.isinf([E1, E2, G12])):\n",
    "            return np.inf  # Return a large value if any modulus is invalid\n",
    "        \n",
    "        # Normalize the values to avoid extreme numbers\n",
    "        E1 /= 1e9\n",
    "        E2 /= 1e9\n",
    "        G12 /= 1e9\n",
    "        \n",
    "        return -(lambda_1 * E1 + lambda_2 * E2 + lambda_3 * G12)  # Negative for maximization since the optimizer is a minimizer\n",
    "    except Exception as e:\n",
    "        print(f\"Error in objective function: {e}\")\n",
    "        return np.inf\n",
    "\n",
    "def constraint_alpha_lower(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    alpha = alpha_l / l\n",
    "    return alpha - 0.5\n",
    "\n",
    "def constraint_alpha_upper(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    alpha = alpha_l / l\n",
    "    return 1 - alpha\n",
    "\n",
    "def constraint_beta_lower(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    beta = beta_l / l\n",
    "    return beta - 0.004\n",
    "\n",
    "def constraint_beta_upper(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    beta = beta_l / l\n",
    "    return 0.008 - beta\n",
    "\n",
    "def constraint_theta_lower(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    return theta_deg + 30\n",
    "\n",
    "def constraint_theta_upper(H):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    return 30 - theta_deg\n",
    "\n",
    "def calculate_relative_density(alpha_l, beta_l, theta_deg):\n",
    "    t = beta_l\n",
    "    return (2 / np.sqrt(3)) * (t / l) * (1 - (t / (2 * np.sqrt(3) * l)))\n",
    "\n",
    "def constraint_relative_density(H, rho_ori):\n",
    "    alpha_l, beta_l, theta_deg = H\n",
    "    rho_star = calculate_relative_density(alpha_l, beta_l, theta_deg)\n",
    "    return rho_ori - rho_star\n",
    "\n",
    "def run_optimization(initial_guess, rho_ori):\n",
    "    try:\n",
    "        bounds = [(3.0, 12.0), (0.024, 0.048), (-30, 30)]\n",
    "        constraints = [\n",
    "            {'type': 'ineq', 'fun': constraint_alpha_lower},\n",
    "            {'type': 'ineq', 'fun': constraint_alpha_upper},\n",
    "            {'type': 'ineq', 'fun': constraint_beta_lower},\n",
    "            {'type': 'ineq', 'fun': constraint_beta_upper},\n",
    "            {'type': 'ineq', 'fun': constraint_theta_lower},\n",
    "            {'type': 'ineq', 'fun': constraint_theta_upper},\n",
    "            {'type': 'ineq', 'fun': lambda H: constraint_relative_density(H, rho_ori)},\n",
    "        ]\n",
    "\n",
    "        solution_slsqp = minimize(objective, initial_guess, method='SLSQP', bounds=bounds, constraints=constraints, options={'disp': True})\n",
    "\n",
    "        return solution_slsqp\n",
    "    except Exception as e:\n",
    "        print(f\"Error during optimization: {e}\")\n",
    "        return None\n",
    "\n",
    "def plot_hexagon(ax, alpha_l, beta_l, theta_deg, color, linestyle, label):\n",
    "    theta_rad = np.radians(theta_deg)\n",
    "    h = alpha_l\n",
    "    t = beta_l\n",
    "\n",
    "    # Regular hexagon vertices calculation\n",
    "    x = l * np.cos(np.linspace(0, 2 * np.pi, 7))\n",
    "    y = l * np.sin(np.linspace(0, 2 * np.pi, 7))\n",
    "\n",
    "    # Rotate the hexagon 180 degrees\n",
    "    x = -x\n",
    "    y = -y\n",
    "\n",
    "    ax.plot(x, y, color=color, linestyle=linestyle, label=label)\n",
    "\n",
    "def plot_results(solution_slsqp):\n",
    "    fig, ax = plt.subplots()\n",
    "    \n",
    "    if solution_slsqp.success:\n",
    "        plot_hexagon(ax, *solution_slsqp.x, 'b', '-', 'SLSQP Optimized cell')\n",
    "        \n",
    "    ax.set_aspect('equal')\n",
    "    ax.legend()\n",
    "    plt.title(\"Optimized Hexagon Shape (Rotated 180 degrees)\")\n",
    "    plt.xlabel(\"X-axis\")\n",
    "    plt.ylabel(\"Y-axis\")\n",
    "    plt.show()\n",
    "\n",
    "# Initial guesses for optimization\n",
    "initial_guess = [6.0, 0.04, 30]  # Example initial guess\n",
    "rho_ori = 0.0153  # Original relative density\n",
    "\n",
    "solution_slsqp = run_optimization(initial_guess, rho_ori)\n",
    "plot_results(solution_slsqp)\n"
   ]
  },
  {
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}