Last active
December 11, 2024 20:52
-
-
Save splch/c252f8ef58374cf2ef41b2863fcf4bef to your computer and use it in GitHub Desktop.
use a gaussian mixture model to fit results
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Mean absolute error: 0.001643681639981239\n", | |
| "Median absolute error: 0.0010639625051655334\n", | |
| "GMM parameters: {'weights': array([0.16367951, 0.0841679 , 0.0800853 , 0.09506292, 0.12197406,\n", | |
| " 0.13252617, 0.13591838, 0.18658577]), 'means': array([[ 35.63736163],\n", | |
| " [180.88951189],\n", | |
| " [120.83791113],\n", | |
| " [233.48248211],\n", | |
| " [ 13.7775845 ],\n", | |
| " [100.23318107],\n", | |
| " [159.81852624],\n", | |
| " [ 51.14557047]]), 'covariances': array([[[ 9.16729036]],\n", | |
| "\n", | |
| " [[158.38423532]],\n", | |
| "\n", | |
| " [[156.72131461]],\n", | |
| "\n", | |
| " [[147.97689595]],\n", | |
| "\n", | |
| " [[ 86.23842402]],\n", | |
| "\n", | |
| " [[117.68017745]],\n", | |
| "\n", | |
| " [[136.27107456]],\n", | |
| "\n", | |
| " [[ 80.85641489]]]), 'precisions': array([[[0.10908349]],\n", | |
| "\n", | |
| " [[0.00631376]],\n", | |
| "\n", | |
| " [[0.00638075]],\n", | |
| "\n", | |
| " [[0.00675781]],\n", | |
| "\n", | |
| " [[0.01159576]],\n", | |
| "\n", | |
| " [[0.00849761]],\n", | |
| "\n", | |
| " [[0.00733831]],\n", | |
| "\n", | |
| " [[0.0123676 ]]])}\n", | |
| "Space saved: 98.16%\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAJOCAYAAACqS2TfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAACWIklEQVR4nOzdd3gU5drH8d8u6UAIJRCCQEKRXgQEgyiokQRQCSL9SBFBgSAYRV8QCUXFBiJdREHPIYIooiIiIbQjRDooVUAQBEIRIdJS5/2DkzVLNskmZLMp38915cruzDMz98w+O7P3PM/MmAzDMAQAAAAAAPKc2dkBAAAAAABQVJF0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdANAAWQymTR+/Hhnh3Hb/v3vf6tu3bpydXWVj49Pvi23f//+CggIyNW048ePl8lkytuAbrF+/XqZTCatX7/eocvJiYCAAPXv39/ZYaAAa9eundq1a+eQeefXPs/Wd69du3Zq2LChw5ctScePH5fJZNLChQvzZXkACgaSbgAF0tGjR/XMM8+oRo0a8vDwkLe3t+699169//77un79urPDgx0OHjyo/v37q2bNmvrwww81b968bKfZtGmTunTpokqVKsnd3V0BAQF65plndOLEiXyIGPZKSEjQjBkz1KZNG5UtW1Zubm7y9/fXY489ps8++0wpKSmWsmlJhslk0muvvWZzfn369JHJZFKpUqWshrdr104mk0m1a9e2OV10dLRl3l988UW2caedUEn7c3V1VUBAgJ577jldunTJ/g1QgG3evFnjx4936voEBARYtrHZbJaPj48aNWqkwYMHa8uWLXm2nKioKE2bNi3P5peXCnJsAPKfi7MDAIBbfffdd+rWrZvc3d3Vt29fNWzYUImJifrxxx81atQo7du3z64ErjC7fv26XFwK9y56/fr1Sk1N1fvvv69atWplW37GjBkaMWKEatSooeHDh6ty5co6cOCA5s+fryVLlmjlypVq3bq1Xcv+8MMPlZqamqu4x44dq//7v//L1bTFwfnz59WhQwft2LFDISEhGjt2rMqVK6e4uDitWbNGvXv31pEjR/Tqq69aTefh4aHPPvtMY8eOtRp+9epVff311/Lw8LC5PA8PDx05ckRbt25Vy5YtrcYtWrRIHh4eunHjRo7WYc6cOSpVqpSuXr2qmJgYzZgxQzt37tSPP/6Yo/kURJs3b9aECRPUv3//fO1dcqumTZvqhRdekCT9/fffOnDggJYuXaoPP/xQzz//vKZOnWpVPjf7vKioKO3du1cjR460e5r7779f169fl5ubW46WlVOZxVa9enVdv35drq6uDl0+gIKlcP+iA1DkHDt2TD179lT16tW1du1aVa5c2TJu2LBhOnLkiL777jsnRug4qampSkxMlIeHR6YJSGFy7tw5SbLrh/+mTZs0cuRItWnTRqtWrZKXl5dl3JAhQ3TvvffqiSee0L59+1S2bNlM53P16lWVLFnytn7Quri4FPoTHo705JNPateuXfryyy/1+OOPW40bPXq0tm/frkOHDmWYrmPHjlq2bJn27NmjJk2aWIZ//fXXSkxMVGhoqNauXZthupo1ayo5OVmfffaZVdJ948YNffXVV+rUqZO+/PLLHK3DE088oQoVKkiSnnnmGfXs2VNLliyxmdg7SlpdLaqqVKmif/3rX1bD3nrrLfXu3VvvvfeeateurSFDhljGOXqfd+PGDbm5uclsNjt1/2oymYrE/h1AztC9HECB8vbbb+vKlSv66KOPrBLuNLVq1dKIESMs75OTkzVp0iTVrFnT0h15zJgxSkhIsJouICBAjzzyiNavX68WLVrI09NTjRo1slzXt2zZMjVq1EgeHh5q3ry5du3aZTV9//79VapUKf32228KCQlRyZIl5e/vr4kTJ8owDKuy7777rlq3bq3y5cvL09NTzZs3t9n11WQyKTw8XIsWLVKDBg3k7u6uVatWWcalv77x77//1siRIxUQECB3d3dVrFhRDz/8sHbu3Gk1z6VLl6p58+by9PRUhQoV9K9//UunTp2yuS6nTp1SWFiYSpUqJV9fX7344otW3YKzMnv2bEvM/v7+GjZsmFV31oCAAEVGRkqSfH19s71ec9KkSTKZTPrkk0+sEm7pZtL19ttv68yZM/rggw8yrMfRo0fVsWNHlS5dWn369LGMu/Wa7j///FNPPvmkvL295ePjo379+mnPnj0Zrq+0dU132me1fPlyNWzYUO7u7mrQoIHl80rz+++/a+jQoapTp448PT1Vvnx5devWTcePH89mi9pm7/wWLlwok8mkTZs2KSIiQr6+vipZsqS6dOmi8+fPW5U1DEOvvfaa7rjjDnl5eemBBx7Qvn377IonNjZWP/zwgwYPHpwh4U7TokULy+eQXlBQkAIDAxUVFWU1fNGiRQoNDVW5cuUyXW6vXr20ZMkSq94L3377ra5du6bu3bvbFXtW7rvvPkk3L2tJb8uWLQoNDVWZMmXk5eWltm3batOmTRmmP3XqlAYOHCh/f3+5u7srMDBQQ4YMUWJioqR/Pp8NGzZo6NChqlixou644w7L9N9//73uu+8+lSxZUqVLl1anTp0yfCY///yz+vfvb7nkxs/PT0899ZT+/PNPS5nx48dr1KhRkqTAwEBLF+/09eU///mPZR9Rrlw59ezZUydPnsywTvPmzVPNmjXl6empli1b6r///W8Ot2pGnp6e+ve//61y5crp9ddft9p35nSf165dO3333Xf6/fffLeuZ9p1Pu2578eLFGjt2rKpUqSIvLy/Fx8dneT+FHTt2qHXr1vL09FRgYKDmzp1rNT7tc7z1+3frPLOKLbNruteuXWupAz4+PurcubMOHDhgVSZt33TkyBFLT4YyZcpowIABunbtmlXZ6OhotWnTRj4+PipVqpTq1KmjMWPGZPHpAHAkTuUDKFC+/fZb1ahRw+5uxE8//bQ++eQTPfHEE3rhhRe0ZcsWTZ48WQcOHNBXX31lVfbIkSPq3bu3nnnmGf3rX//Su+++q0cffVRz587VmDFjNHToUEnS5MmT1b17dx06dEhm8z/nJlNSUhQaGqp77rlHb7/9tlatWqXIyEglJydr4sSJlnLvv/++HnvsMfXp00eJiYlavHixunXrphUrVqhTp05WMa1du1aff/65wsPDVaFChUxv/vXss8/qiy++UHh4uOrXr68///xTP/74ow4cOKBmzZpJuvmDcMCAAbr77rs1efJknT17Vu+//742bdqkXbt2WbU4p6SkKCQkRK1atdK7776rNWvWaMqUKapZs6ZV65Mt48eP14QJExQcHKwhQ4bo0KFDmjNnjrZt26ZNmzbJ1dVV06ZN06effqqvvvrK0pW3cePGNud37do1xcTE6L777lNgYKDNMj169NDgwYO1YsUKq67fycnJCgkJUZs2bfTuu+9mSNjTpKam6tFHH9XWrVs1ZMgQ1a1bV19//bX69euX5bqm9+OPP2rZsmUaOnSoSpcurenTp6tr1646ceKEypcvL0natm2bNm/erJ49e+qOO+7Q8ePHNWfOHLVr10779+/PNL7M5HR+w4cPV9myZRUZGanjx49r2rRpCg8P15IlSyxlxo0bp9dee00dO3ZUx44dtXPnTrVv396SIGbl22+/laQMLZj26tWrl/7zn//ozTfflMlk0oULF7R69Wr9+9//znACI73evXtr/PjxWr9+vR588EFJN7vvPvTQQ6pYsWKuYkkvLYlK34ti7dq16tChg5o3b67IyEiZzWYtWLBADz74oP773/9aWsRPnz6tli1b6tKlSxo8eLDq1q2rU6dO6YsvvtC1a9esujEPHTpUvr6+GjdunK5evSrp5s0G+/Xrp5CQEL311lu6du2a5syZozZt2mjXrl2WfUJ0dLR+++03DRgwQH5+fpbLbPbt26effvpJJpNJjz/+uH799Vd99tlneu+99yyt+b6+vpKk119/Xa+++qq6d++up59+WufPn9eMGTN0//33W+0jPvroIz3zzDNq3bq1Ro4cqd9++02PPfaYypUrp6pVq97Wti5VqpS6dOmijz76SPv371eDBg1slstun/fKK6/o8uXL+uOPP/Tee+9Z5p3epEmT5ObmphdffFEJCQlZdin/66+/1LFjR3Xv3l29evXS559/riFDhsjNzU1PPfVUjtbRntjSW7NmjTp06KAaNWpo/Pjxun79umbMmKF7771XO3fuzHBc6N69uwIDAzV58mTt3LlT8+fPV8WKFfXWW29Jkvbt26dHHnlEjRs31sSJE+Xu7q4jR47YPGEEIJ8YAFBAXL582ZBkdO7c2a7yu3fvNiQZTz/9tNXwF1980ZBkrF271jKsevXqhiRj8+bNlmE//PCDIcnw9PQ0fv/9d8vwDz74wJBkrFu3zjKsX79+hiRj+PDhlmGpqalGp06dDDc3N+P8+fOW4deuXbOKJzEx0WjYsKHx4IMPWg2XZJjNZmPfvn0Z1k2SERkZaXlfpkwZY9iwYZlui8TERKNixYpGw4YNjevXr1uGr1ixwpBkjBs3LsO6TJw40Woed911l9G8efNMl2EYhnHu3DnDzc3NaN++vZGSkmIZPnPmTEOS8fHHH1uGRUZGGpKsto0taZ/jiBEjsizXuHFjo1y5chnW4//+7/8ylO3Xr59RvXp1y/svv/zSkGRMmzbNMiwlJcV48MEHDUnGggULMsSdniTDzc3NOHLkiGXYnj17DEnGjBkzLMNu/ewNwzBiY2MNScann35qGbZu3boMdcwWe+e3YMECQ5IRHBxspKamWoY///zzRokSJYxLly4ZhvHP59epUyercmPGjDEkGf369csyni5duhiSLPNLc/36deP8+fOWv7/++ssy7tixY4Yk45133jH27t1rSDL++9//GoZhGLNmzTJKlSplXL161ejXr59RsmRJq/m2bdvWaNCggWEYhtGiRQtj4MCBhmEYxl9//WW4ubkZn3zyiWVbLl26NMvYDeOfz/bQoUPG+fPnjePHjxsff/yx4enpafj6+hpXr141DOPmd7t27dpGSEiI1Xa6du2aERgYaDz88MOWYX379jXMZrOxbdu2DMtLmzbt82nTpo2RnJxsGf/3338bPj4+xqBBg6ymi4uLM8qUKWM13FZd+OyzzwxJxsaNGy3D3nnnHUOScezYMauyx48fN0qUKGG8/vrrVsN/+eUXw8XFxTI8bV/StGlTIyEhwVJu3rx5hiSjbdu2GeK4VfXq1Y1OnTplOv69994zJBlff/21ZVhO93mGYRidOnWy+p6nSasTNWrUyLDdbH332rZta0gypkyZYhmWkJBgNG3a1KhYsaKRmJhoGMY/n+Ot29bWPDOLLe37kH6fk7acP//80zJsz549htlsNvr27WsZllZ/n3rqKat5dunSxShfvrzlfdr2zW7fCyD/0L0cQIERHx8vSSpdurRd5VeuXClJioiIsBqedvOeW6/9rl+/voKCgizvW7VqJUl68MEHVa1atQzDf/vttwzLDA8Pt7xO63KcmJioNWvWWIZ7enpaXv/111+6fPmy7rvvvgxdwSWpbdu2ql+/fjZrevO66C1btuj06dM2x2/fvl3nzp3T0KFDra4X7NSpk+rWrWvzOvhnn33W6v19991nc53TW7NmjRITEzVy5EirXgCDBg2St7d3rq63//vvvyVl/7mXLl3aUkfSy65lXpJWrVolV1dXDRo0yDLMbDZr2LBhdscZHBysmjVrWt43btxY3t7eVtss/WeflJSkP//8U7Vq1ZKPj4/Nzz87OZ3f4MGDrbrG33fffUpJSdHvv/8u6Z/Pb/jw4Vbl7L0RVdr2v7XVbu7cufL19bX8tWnTxub0DRo0UOPGjfXZZ59Jutla3blzZ7t6APTu3VvLli1TYmKivvjiC5UoUUJdunSxK+5b1alTR76+vgoICNBTTz2lWrVq6fvvv7fEsXv3bh0+fFi9e/fWn3/+qQsXLujChQu6evWqHnroIW3cuFGpqalKTU3V8uXL9eijj6pFixYZlnPrZQqDBg1SiRIlLO+jo6N16dIl9erVy7KMCxcuqESJEmrVqpXWrVtnKZu+Lty4cUMXLlzQPffcI0l21a1ly5YpNTVV3bt3t1qWn5+fateubVlW2r7k2WeftWoZ7t+/v8qUKWPP5s1WWv1J++7bkt0+zx79+vWz2m5ZcXFx0TPPPGN57+bmpmeeeUbnzp3Tjh07ch1Dds6cOaPdu3erf//+VpdYNG7cWA8//LDlOJeerX33n3/+afl+pvVY+Prrr3N9Q0kAeYukG0CB4e3tLSnrH2Lp/f777zKbzRnujO3n5ycfHx9LopEmfWItyfID8tbukmnD//rrL6vhZrNZNWrUsBp25513SpLVNX4rVqzQPffcIw8PD5UrV06+vr6aM2eOLl++nGEdMutOfau3335be/fuVdWqVdWyZUuNHz/eKtlLW9c6depkmLZu3boZtoWHh4ely2masmXLZljnW2W2HDc3N9WoUSPDcuyRlmxn97n//fffGRJzFxcXq2tjM/P777+rcuXKGZI7e+6qnubW+iNl3GbXr1/XuHHjVLVqVbm7u6tChQry9fXVpUuXbH7+2cnp/G6NMa27dFqMaZ/PrY/g8vX1zfIGdWnStv+VK1eshnft2lXR0dGKjo7O9DKCNL1799bSpUt15MgRbd68Wb179852uZLUs2dPXb58Wd9//70WLVqkRx55xOaJmsTERMXFxVn93Xqvgi+//FLR0dGKiorSPffco3PnzlklZ4cPH5Z0M2lLfzLB19dX8+fPV0JCgi5fvqzz588rPj7e7mc83/p9T1vOgw8+mGE5q1evttyMUJIuXryoESNGqFKlSvL09JSvr69lfvbUrcOHD8swDNWuXTvDsg4cOGBZVmZ1xNXVNcP+L7fS6k9WJ9qy2+fZw979qyT5+/tnuLGdrf17Xstq312vXj3LyZ70svue9+jRQ/fee6+efvppVapUST179tTnn39OAg44Edd0AygwvL295e/vr7179+ZoultbkzKTvoXJnuHGLTdIs8d///tfPfbYY7r//vs1e/ZsVa5cWa6urlqwYEGGG0hJsrsVpnv37rrvvvv01VdfafXq1XrnnXf01ltvadmyZerQoUOO48xsnZ2hVq1acnFx0c8//5xpmYSEBB06dChDa6K7u7tVi7sj2VNPhg8frgULFmjkyJEKCgpSmTJlZDKZ1LNnz1z94M3p/PKyLttSt25dSdLevXt17733WoZXrVrVcvKqbNmyunDhQqbz6NWrl0aPHq1BgwapfPnyat++vV3Lrly5stq1a6cpU6Zo06ZNmd6xfPPmzXrggQeshh07dszqutj777/fcr3zo48+qkaNGqlPnz7asWOHzGazZdu+8847atq0qc3llCpVShcvXrQr9jS3ft/TlvPvf/9bfn5+Gcqnv4t+9+7dtXnzZo0aNUpNmzZVqVKllJqaqtDQULvqVmpqqkwmk77//nub9SSra47zWto+PquTXnmxz7N3/2qvzI419t6AMq9k9z339PTUxo0btW7dOn333XdatWqVlixZogcffFCrV68uUPt/oLgg6QZQoDzyyCOaN2+eYmNjrbqC21K9enWlpqbq8OHDqlevnmX42bNndenSJVWvXj1PY0tNTdVvv/1maf2QpF9//VWSLD/ov/zyS3l4eOiHH36Qu7u7pdyCBQtue/mVK1fW0KFDNXToUJ07d07NmjXT66+/rg4dOljW9dChQ5YbTaU5dOhQnm2L9MtJ3+qVmJioY8eOKTg4OMfzLFmypB544AGtXbtWv//+u81YP//8cyUkJOiRRx7Jddzr1q3TtWvXrFq7jxw5kqv5ZeaLL75Qv379NGXKFMuwGzduWN3Z3ZnzS9u2hw8ftvr8zp8/n20vB+nm9/PNN9/UokWLrJLunKhWrZruvfderV+/XkOGDMnR49l69+6tp59+Wj4+PurYsaPNMk2aNFF0dLTVMFsJbZpSpUopMjJSAwYM0Oeff66ePXtaLiPw9vbOsk77+vrK29s7xycK06Qtp2LFilku56+//lJMTIwmTJigcePGWYantZSnl1liWLNmTRmGocDAQKt92K3S15H0+5KkpCQdO3bM6nFvuXHlyhV99dVXqlq1qtV+25as9nmS/Sdc7XH69OkMj3G7df+e1qJ86/fPVg8fe2NLv0+91cGDB1WhQoVcPVrObDbroYce0kMPPaSpU6fqjTfe0CuvvKJ169blaj8N4PbQvRxAgfLSSy+pZMmSevrpp3X27NkM448ePar3339fkiw/uqdNm2ZVZurUqZKU4U7heWHmzJmW14ZhaObMmXJ1ddVDDz0k6WYLhMlksmr5OH78uJYvX57rZaakpGToPlqxYkX5+/tbHo3WokULVaxYUXPnzrV6XNr333+vAwcO5Nm2CA4Olpubm6ZPn27VevrRRx/p8uXLuV7O2LFjZRiG+vfvr+vXr1uNO3bsmF566SVVrlzZ6prLnAgJCVFSUpI+/PBDy7DU1FTNmjUrV/PLTIkSJTK0Ks+YMSPXLWF5Pb/g4GC5urpqxowZVvO99TuUmXvvvVcPP/yw5s2bp6+//tpmGXta1V977TVFRkZq+PDhdi03zRNPPKHIyEjNnj070ztRly1bVsHBwVZ/2T0XuU+fPrrjjjssd39u3ry5atasqXfffTdDV3pJlsewmc1mhYWF6dtvv9X27dszlMtuW4SEhMjb21tvvPGGkpKSMl1OWsvkrfOz9bmlJWi3JoaPP/64SpQooQkTJmSYj2EYlkePtWjRQr6+vpo7d67VHe0XLlyY65M9aa5fv64nn3xSFy9e1CuvvJJly3F2+zzp5rrm5rINW5KTk60eSZiYmKgPPvhAvr6+at68uaR/TpJs3LjRKtZ58+ZlmJ+9sVWuXFlNmzbVJ598YrV99+7dq9WrV2d6cikrtnpgpPXYuPVxmgDyBy3dAAqUmjVrKioqSj169FC9evXUt29fNWzYUImJidq8ebOWLl2q/v37S7rZotWvXz/NmzdPly5dUtu2bbV161Z98sknCgsLy9DF9HZ5eHho1apV6tevn1q1aqXvv/9e3333ncaMGWO5PrpTp06aOnWqQkND1bt3b507d06zZs1SrVq1suw+nZW///5bd9xxh5544gk1adJEpUqV0po1a7Rt2zZLC6irq6veeustDRgwQG3btlWvXr0sjwwLCAjQ888/nyfbwNfXV6NHj9aECRMUGhqqxx57TIcOHdLs2bN199135/pRUvfff7/effddRUREqHHjxurfv78qV66sgwcP6sMPP1RqaqpWrlxp13XHtoSFhally5Z64YUXdOTIEdWtW1fffPON5cdpXrWYPfLII/r3v/+tMmXKqH79+oqNjdWaNWssjxRz9vzSnsc+efJkPfLII+rYsaN27dql77//3tLdOjv/+c9/FBoaqrCwMHXo0EHBwcEqW7as4uLitGbNGm3cuDHb7r9t27ZV27Ztcxx/mTJlsnzee265urpqxIgRGjVqlFatWqXQ0FDNnz9fHTp0UIMGDTRgwABVqVJFp06d0rp16+Tt7W15fNobb7yh1atXq23btho8eLDq1aunM2fOaOnSpfrxxx+tHtV3K29vb82ZM0dPPvmkmjVrpp49e8rX11cnTpzQd999p3vvvVczZ86Ut7e37r//fr399ttKSkpSlSpVtHr1ah07dizDPNMSxFdeeUU9e/aUq6urHn30UdWsWVOvvfaaRo8erePHjyssLEylS5fWsWPH9NVXX2nw4MF68cUX5erqqtdee03PPPOMHnzwQfXo0UPHjh3TggULcnRN96lTp/Sf//xH0s3W7f3792vp0qWKi4vTCy+8kOUJNHv2eWnrumTJEkVEROjuu+9WqVKl9Oijj9odY3r+/v566623dPz4cd15551asmSJdu/erXnz5snV1VXSzRsB3nPPPRo9erQuXryocuXKafHixUpOTs4wv5zE9s4776hDhw4KCgrSwIEDLY8My219nzhxojZu3KhOnTqpevXqOnfunGbPnq077rgj05scAnCwfL9fOgDY4ddffzUGDRpkBAQEGG5ubkbp0qWNe++915gxY4Zx48YNS7mkpCRjwoQJRmBgoOHq6mpUrVrVGD16tFUZw8j8ETaSMjyWJv0jjtKkPc7o6NGjRvv27Q0vLy+jUqVKRmRkpNWjswzDMD766COjdu3ahru7u1G3bl1jwYIFmT6GKrNH4ijd43MSEhKMUaNGGU2aNDFKly5tlCxZ0mjSpIkxe/bsDNMtWbLEuOuuuwx3d3ejXLlyRp8+fYw//vjDqoytRzMZhu1HZWVm5syZRt26dQ1XV1ejUqVKxpAhQ6weE5V+fjl5bM3GjRuNzp07GxUqVDBcXV2NatWqGYMGDTKOHz+eoWxm65E27tbH9Zw/f97o3bu3Ubp0aaNMmTJG//79jU2bNhmSjMWLF2eIO73MPqvq1atbPWbrr7/+MgYMGGBUqFDBKFWqlBESEmIcPHgwQzl7Hxlm7/zSHmV062OrbC0nJSXFmDBhglG5cmXD09PTaNeunbF3794M88zK9evXjWnTphlBQUGGt7e34eLiYvj5+RmPPPKIsWjRIqvHYtn6PtmS3SPDMpObR4bZqpOXL182ypQpY/VIrF27dhmPP/64Ub58ecPd3d2oXr260b17dyMmJsZq2t9//93o27ev4evra7i7uxs1atQwhg0bZnnkVmafT/p1CAkJMcqUKWN4eHgYNWvWNPr3729s377dUuaPP/4wunTpYvj4+BhlypQxunXrZpw+fTrDo7YMwzAmTZpkVKlSxTCbzRkecfXll18abdq0MUqWLGmULFnSqFu3rjFs2DDj0KFDVvOYPXu2ERgYaLi7uxstWrQwNm7caLRt29buR4ZJMiQZJpPJ8Pb2Nho0aGAMGjTI2LJli81pcrPPu3LlitG7d2/Dx8fHkGT5zmdVJzJ7ZFiDBg2M7du3G0FBQYaHh4dRvXp1Y+bMmRmmP3r0qBEcHGy4u7sblSpVMsaMGWNER0dnmGdmsdl6ZJhhGMaaNWuMe++91/D09DS8vb2NRx991Ni/f79Vmczq762PMouJiTE6d+5s+Pv7G25uboa/v7/Rq1cv49dff7W57QE4nskw8ujuKgBQhPXv319ffPGFza6mKLyWL1+uLl266Mcff8z1NcoAAABZ4ZpuAECxcOu14ikpKZoxY4a8vb3VrFkzJ0UFAACKOq7pBgAUC8OHD9f169cVFBSkhIQELVu2TJs3b9Ybb7yR548WAgAASEPSDQAoFh588EFNmTJFK1as0I0bN1SrVi3NmDFD4eHhzg4NAAAUYVzTDQAAAACAg3BNNwAAAAAADkLSDQAAAACAg3BNdy6lpqbq9OnTKl26tEwmk7PDAQAAAADkI8Mw9Pfff8vf319mc+bt2STduXT69GlVrVrV2WEAAAAAAJzo5MmTuuOOOzIdT9KdS6VLl5Z0cwN7e3s7ORrbkpKStHr1arVv316urq7ODgfFEHUQzkT9g7NRB+FM1D84U3Gpf/Hx8apataolN8wMSXcupXUp9/b2LtBJt5eXl7y9vYt0ZUfBRR2EM1H/4GzUQTgT9Q/OVNzqX3aXG3MjNQAAAAAAHISkGwAAAAAAByHpBgAAAADAQbimGwAAAHCAlJQUJSUlOWXZSUlJcnFx0Y0bN5SSkuKUGFB8FZX65+rqqhIlStz2fEi6AQAAgDxkGIbi4uJ06dIlp8bg5+enkydPZnuTJyCvFaX65+PjIz8/v9taD5JuAAAAIA+lJdwVK1aUl5eXU5KO1NRUXblyRaVKlZLZzBWlyF9Fof4ZhqFr167p3LlzkqTKlSvnel4k3QAAAEAeSUlJsSTc5cuXd1ocqampSkxMlIeHR6FNelB4FZX65+npKUk6d+6cKlasmOuu5oV3CwAAAAAFTNo13F5eXk6OBEBeSPsu3879GUi6AQAAgDxW2K9jBXBTXnyXSboBAAAAAHAQkm4AAAAAt+X48eMymUzavXu33dMsXLhQPj4+To8jP+Z1q4CAAE2bNs3y3mQyafny5Xm+HFvLQv7jRmoAAACAgw1cuC1fl/dh3+Y5nubkyZOKjIzUqlWrdOHCBVWuXFlhYWEaN25ctjeFq1q1qs6cOaMKFSrYvbwePXqoY8eOOY7zdrVr104bNmyQJLm5ualChQpq1qyZBgwYoMcff9xSLifrdPz4cQUGBmrXrl1q2rRptuW3bdumkiVL5nodbFm4cKFGjhyZ4VF1jlgWcoaWbgAAAKCY++2339SiRQsdPnxYn332mY4cOaK5c+cqJiZGQUFBunjxYqbTJiYmqkSJEvLz85OLi/1tep6enqpYsWJehJ9jgwYN0pkzZ3T06FF9+eWXql+/vnr27KnBgwdbyuRmnbKTmJgoSfL19c23m+3l57JgG0k3AAAAUMwNGzZMbm5uWr16tdq2batq1aqpQ4cOWrNmjU6dOqVXXnnFUjYgIECTJk1S37595e3trcGDB9vsiv3NN9+odu3a8vDw0AMPPKBPPvlEJpPJ0hJ7a/fy8ePHq2nTpvr3v/+tgIAAlSlTRj179tTff/9tKbNq1Sq1adNGPj4+Kl++vB555BEdPXo0x+vr5eUlPz8/3XHHHbrnnnv01ltv6YMPPtCHH36oNWvWSMrYvfyvv/5Snz595OvrK09PT9WuXVsLFiyQJAUGBkqS7rrrLplMJrVr106S1L9/f4WFhen111+Xv7+/6tSpY9mGt3b5PnPmjDp06CBPT0/VqFFDX3zxhWXc+vXrrbadJO3evVsmk0nHjx/X+vXrNWDAAF2+fFkmk0kmk0njx4+3uawTJ06oc+fOKlWqlLy9vdW9e3edPXs2R58DcqZAJN2zZs1SQECAPDw81KpVK23dujXL8kuXLlXdunXl4eGhRo0aaeXKlVbjx48fr7p166pkyZIqW7asgoODtWXLFqsyFy9eVJ8+feTt7S0fHx8NHDhQV65cyfN1AwAAAAqyixcv6ocfftDQoUMtzyVO4+fnpz59+mjJkiUyDMMy/N1331WTJk20a9cuvfrqqxnmeezYMT3xxBMKCwvTnj179Mwzz1gl7pk5evSoli9frhUrVmjFihXasGGD3nzzTcv4q1evKiIiQtu3b1dMTIzMZrO6dOmi1NTU29gCN/Xr109ly5bVsmXLbI5/9dVXtX//fn3//fc6cOCA5syZY+l6npa/rFmzRmfOnLGaR0xMjA4dOqTo6GitWLEi0+W/+uqr6tq1q/bs2aM+ffqoZ8+eOnDggF2xt27dWtOmTZO3t7fOnDmjM2fO6MUXX8xQLjU1VZ07d9bFixe1YcMGRUdH67ffflOPHj2symX3OSBnnH5N95IlSxQREaG5c+eqVatWmjZtmkJCQnTo0CGb3U02b96sXr16afLkyXrkkUcUFRWlsLAw7dy5Uw0bNpQk3XnnnZo5c6Zq1Kih69ev67333lP79u115MgR+fr6SpL69OmjM2fOKDo6WklJSRowYIAGDx6sqKiofF1/AAAAwJkOHz4swzBUr149m+Pr1aunv/76S+fPn7f8Pn/wwQf1wgsvWMocP37capoPPvhAderU0TvvvCNJqlOnjvbu3avXX389y1hSU1O1cOFClS5dWpL05JNPKiYmxjJd165drcp//PHH8vX11f79+y25QG6ZzWbdeeedGdYlzYkTJ3TXXXepRYsWkm62IKdJyzHKly8vPz8/q+lKliyp+fPny83NLcvld+vWTU8//bQkadKkSYqOjtaMGTM0e/bsbGN3c3NTmTJlZDKZMiw/vZiYGP3yyy86duyYqlatKkn69NNP1aBBA23btk133323pOw/B+SM01u6p06dqkGDBmnAgAGqX7++5s6dKy8vL3388cc2y7///vsKDQ3VqFGjVK9ePU2aNEnNmjXTzJkzLWV69+6t4OBg1ahRQw0aNNDUqVMVHx+vn3/+WZJ04MABrVq1SvPnz1erVq3Upk0bzZgxQ4sXL9bp06fzZb0BAACAgiR9S3Z20hLPzBw6dMiSwKVp2bJltvMNCAiwJHqSVLlyZZ07d87y/vDhw+rVq5dq1Kghb29vS+J74sQJu2PPimEYmT6XeciQIVq8eLGaNm2ql156SZs3b7Zrno0aNco24ZakoKCgDO/tbem214EDB1S1alVLwi1J9evXl4+Pj9WysvsckDNObelOTEzUjh07NHr0aMsws9ms4OBgxcbG2pwmNjZWERERVsNCQkIyvcV+YmKi5s2bpzJlyqhJkyaWefj4+FjtLIKDg2U2m7VlyxZ16dIlw3wSEhKUkJBgeR8fHy9JSkpKUlJSkn0rnM/S4iqo8aHoow7Cmah/cDbqYPGUlJQkwzCUmppq1eXZkP0JbV5IS6DTYslKjRo1ZDKZtH//fnXu3DnD+P3796ts2bIqX768ZV5eXl5W8017nbbehmFkWPatZdK/T4vV1dU1Q7zpyz766KOqVq2aPvjgA/n7+ys1NVWNGzfWjRs3Mswzq/W2tV1SUlJ0+PBhtWjRwua8QkJCdOzYMa1cuVJr1qzRQw89pKFDh+qdd97JdLmGYWTYVpnFYGva9NsnLca092m5SWbb09aybM3z1uXb8zlkJyf1r6BL2yZJSUkqUaKE1Th79+9OTbovXLiglJQUVapUyWp4pUqVdPDgQZvTxMXF2SwfFxdnNWzFihXq2bOnrl27psqVKys6OtpyzUVcXFyGrusuLi4qV65chvmkmTx5siZMmJBh+OrVqwv83QCjo6OdHQKKOeognIn6B2ejDhYvLi4u8vPz05UrVyx3qpak5KTkfI0j7aZX9tz8ytXVVQ888IBmz56tp556yuq67rNnzyoqKko9evSwzCs1NVU3btywNEJJstwb6erVq4qPj1dAQICio6OtymzatMkSk9ls1o0bN2QYhqVMQkKCUlJSrKZJS6bj4+N18eJFHTp0SFOnTrW0oqc11F2/fl3x8fEZ4rAlOTlZiYmJGcb/5z//0V9//aXQ0NBM5+Xu7q4uXbqoS5cuatGihSIjI/Xqq69aEuD4+Hir+SYlJSk5OTnDsmxtw//+978KCwuzvN+8ebMaNWqk+Ph4y2dy+PBh1a1bV5L0008/WbZ9fHy8UlJSMmy/W5dVrVo1nTx5Uvv379cdd9whSTp48KAuXbqk6tWrKz4+PtvPISeKws3XEhMTdf36dW3cuFHJydbf42vXrtk1D6df0+0oDzzwgHbv3q0LFy7oww8/VPfu3bVly5ZcP5Zg9OjRVi3s8fHxqlq1qtq3by9vb++8CjtPJSUlKTo6Wg8//LBcXV2dHQ6KIeognIn6B2ejDhZPN27c0MmTJ1WqVCl5eHhYhru45u/P7tKlS+vvv/9W6dKlM+0und7s2bPVpk0b9ejRQxMnTlRgYKD27dunl19+WVWqVNHbb79t+c1rNpvl4eFh9Ru4VKlSkm5ev+zt7a3hw4dr9uzZeuONN/TUU09p9+7dWrx4sSTJ29tb3t7e8vDwkMlksszH3d1dJUqUsJqvh4eHzGazvL29VapUKZUvX15RUVGqVauWTpw4ocjISEk3Hz+WViZ9HLa4uLgoOTlZ165dU3Jysv744w8tX75c06ZN07PPPqtOnTrZXKfIyEg1a9ZMDRo0UEJCgmJiYlSvXj15e3vLy8tLnp6e+vHHH1WnTh15eHioTJkycnV1lYuLS4ZYbG3Db775RkFBQWrTpo2ioqK0Y8cOffzxx/L29lbTpk1VtWpVTZkyRa+99pp+/fVXzZkzxxKnt7e36tWrpytXrmjbtm1q0qSJvLy85OXlZbWsxx57TI0aNdLQoUM1depUJScnKzw8XG3btlXbtm3t+hzsYRhGjupfQXbjxg15enrq/vvvt/pOS7L7JIRTk+4KFSqoRIkSVreol26eUcvsBgB+fn52lS9ZsqRq1aqlWrVq6Z577lHt2rX10UcfafTo0fLz88twTUJycrIuXryY6XLd3d3l7u6eYbirq2uBP5AWhhhRtFEH4UzUPzgbdbB4SUlJkclkktlsltn8z+2TTMrfxCMt0UmLJTt16tTR9u3bFRkZqZ49e1p+F4eFhSkyMlLly5fPMP/08017nbbeNWvW1BdffKEXXnhB06dPV1BQkF555RUNGTJEnp6eVtsn7X9azFbbLd0ws9msxYsX67nnnlPjxo1Vp04dTZ8+Xe3atbOMvzWOzMyfP99yc7Py5curefPmWrJkidVlprfOy93dXa+88oqOHz8uT09P3XfffVq8eLHMZrPc3Nw0ffp0TZw4UZGRkbrvvvssj/nK7DO4dfiECRP0+eefKzw8XJUrV9Znn31muTmcu7u7PvvsMw0ZMkRNmzbV3Xffrddee03dunWzxNemTRs9++yz6tWrl/78809FRkZaHhuWfllff/21hg8fbtluoaGhmjFjht2fgz3SupTbW/8KMrPZLJPJZHNfbu++3WTk5I4JDtCqVSu1bNlSM2bMkHTzA6pWrZrCw8P1f//3fxnK9+jRQ9euXdO3335rGda6dWs1btxYc+fOzXQ5NWvW1JNPPqnx48frwIEDql+/vrZv367mzZtLutlNPDQ0VH/88Yf8/f2zjTs+Pl5lypTR5cuXC3RL98qVK9WxY0cO9gXcwIXb9FH/u7MvWMhQB+FM1D84G3WweLpx44aOHTumwMDADK1i+SmtK7C3t3eBSXpef/11zZ07VydPnnR2KHCwglj/ciur77S9OaHTu5dHRESoX79+atGihVq2bKlp06bp6tWrGjBggCSpb9++qlKliiZPnixJGjFihNq2baspU6aoU6dOWrx4sbZv36558+ZJunnNxeuvv67HHntMlStX1oULFzRr1iydOnVK3bp1k3TzsQehoaEaNGiQ5s6dq6SkJIWHh6tnz552JdwAAAAAsjZ79mzdfffdKl++vDZt2qR33nlH4eHhzg4LyHdOT7p79Oih8+fPa9y4cYqLi1PTpk21atUqy83STpw4YXV2pHXr1oqKitLYsWM1ZswY1a5dW8uXL7d0vShRooQOHjyoTz75RBcuXFD58uV1991367///a8aNGhgmc+iRYsUHh6uhx56SGazWV27dtX06dPzd+UBAACAIurw4cN67bXXdPHiRVWrVk0vvPCC1VOLgOLC6Um3JIWHh2d61mv9+vUZhnXr1s3San0rDw8PLVu2LNtllitXTlFRUTmKEwAAAIB93nvvPb333nvODgNwusLdwR4AAAAAgAKMpBsAAAAAAAch6QYAAAAAwEFIugEAAAAAcBCSbgAAAAAAHISkGwAAAAAAByHpBgAAAIBCICAgQNOmTStw80qvf//+6tKli+V9u3btNHLkyDxfTtqywsLCHDLvvFQgntMNAAAAFGlRPfJ3eT0/y/EkcXFxmjx5sr777jv98ccfKlOmjGrVqqV//etf6tevn7y8vCTdTNZ+//13ffbZZ+rZs6fVPBo0aKD9+/drwYIF6t+/f67K32r8+PGaMGGCJMlsNsvf318dOnTQm2++qXLlyuV4PfNbQECARo4c6bDEM73026pEiRLy8fFR/fr19fjjj2vIkCFyd3e3lN22bZtKlixp13xzsg7vv/++UlJSchV/Zo4fP67AwEDt2rVLTZs2tVqWYRh5uixHoKUbAAAAKOZ+++033XXXXVq9erXeeOMN7dq1S7GxsXrppZe0YsUKrVmzxqp81apVtWDBAqthP/30k+Li4mwmcjktf6sGDRrozJkzOnHihBYsWKBVq1ZpyJAhuVhT+yQmJjps3o6WflutW7dO3bp10+TJk9W6dWv9/ffflnK+vr6WEyl5ISUlRampqSpTpox8fHzybL5Zyc9l3Q6SbgAAAKCYGzp0qFxcXLR9+3Z1795d9erVU40aNdS5c2d99913evTRR63K9+nTRxs2bNDJkyctwz7++GP16dNHLi4ZO9PmtPytXFxc5OfnpypVqig4OFjdunVTdHS0VZn58+erXr168vDwUN26dTV79myr8X/88Yd69eqlcuXKqWTJkmrRooW2bNki6WYLcdOmTTV//nwFBgbKw8NDknTp0iU9/fTT8vX1lbe3tx588EHt2bPHMs+jR4+qc+fOqlSpkkqVKqW7777b6gRFu3bt9Pvvv+v555+XyWSSyWSyjPvxxx913333ydPTU1WrVtVzzz2nq1evWsafO3dOjz76qDw9PRUYGKhFixZlu53Sbyt/f381atRIw4cP14YNG7R371699dZblnLpu5cbhqHx48erWrVqcnd3l7+/v5577rks12HhwoXy8fHRN998o/r168vd3V0nTpzI0L1ckpKTkxUeHq4yZcqoQoUKevXVV61aqE0mk5YvX241jY+PjxYuXChJCgwMlCTdddddMplMateunaSM3csTEhL03HPPqWLFivLw8FCbNm20bds2y/j169fLZDIpJiZGLVq0kJeXl1q3bq1Dhw7ZtW1zi6QbAAAAKMb+/PNPrV69WsOGDcu01Tl9sihJlSpVUkhIiD755BNJ0rVr17RkyRI99dRTNqfPafmsHD9+XD/88IPc3NwswxYtWqRx48bp9ddf14EDB/TGG2/o1VdftSzvypUratu2rU6dOqVvvvlGe/bs0UsvvaTU1FTLPI4cOaIvv/xSy5Yt0+7duyVJ3bp107lz5/T9999rx44datasmR566CFdvHjRMt+OHTsqJiZGu3btUmhoqB599FGdOHFCkrRs2TLdcccdmjhxos6cOaMzZ85Iupmsh4aGqmvXrvr555+1ZMkS/fjjjwoPD7fE079/f508eVLr1q3TF198odmzZ+vcuXM53l6SVLduXXXo0EHLli2zOf7LL7/Ue++9pw8++ECHDx/W8uXL1ahRoyzXQbr5Ob711luaP3++9u3bp4oVK9qc/yeffCIXFxdt3bpV77//vqZOnar58+fbHf/WrVslSWvWrNGZM2cyXY+XXnpJX375pT755BPt3LlTtWrVUkhIiOXzSvPKK69oypQp2r59u1xcXHJVD3OCa7oBAACAYuzIkSMyDEN16tSxGl6hQgXduHFDkjRs2DCrVlJJeuqpp/TCCy/olVde0RdffKGaNWtaXW97q5yWT++XX35RqVKllJKSYolp6tSplvGRkZGaMmWKHn/8cUk3W0b379+vDz74QP369VNUVJTOnz+vbdu2Wa4Dr1WrltUyEhMT9emnn8rX11fSzZborVu36ty5c5Zrod99910tX75cX3zxhQYPHqwmTZqoSZMmlnlMmjRJX331lb755huFh4erXLlyKlGihEqXLi0/Pz9LucmTJ6tPnz6Wa6Rr166t6dOnq23btpozZ45OnDih77//Xlu3btXdd98tSfroo49Ur149u7aXLXXr1tXq1attjjtx4oT8/PwUHBwsV1dXVatWTS1btpSkTNdBkpKSkjR79myrbWBL1apV9d5778lkMqlOnTr65Zdf9N5772nQoEF2xZ72mZQvXz5DDGmuXr2qOXPmaOHCherQoYMk6cMPP1R0dLQ++ugjjRo1ylL29ddfV9u2bSVJ//d//6dOnTrpxo0blh4OeY2WbgAAAAAZbN26Vbt371aDBg2UkJCQYXynTp105coVbdy4UR9//HG2rYU5LZ9enTp1tHv3bm3btk0vv/yyQkJCNHz4cEk3k62jR49q4MCBKlWqlOXvtdde09GjRyVJu3fv1l133ZXljdeqV69uSe4kac+ePbpy5YrKly9vNd9jx45Z5nvlyhW9+OKLqlevnnx8fFSqVCkdOHDA0tKdmT179mjhwoVW8w0JCVFqaqqOHTumAwcOyMXFRc2bN7dMU7du3du6ftkwjAw9FtJ069ZN169fV40aNTRo0CB99dVXSk5Oznaebm5uaty4cbbl7rnnHqtlBwUF6fDhw3l6w7WjR48qKSlJ9957r2WYq6urWrZsqQMHDliVTR9z5cqVJSnXvQjsQUs3AAAAUIzVqlVLJpMpw3WtNWrUkCR5enranM7FxUVPPvmkIiMjtWXLFn311VdZLien5dNzc3OztEy/+eab6tSpkyZMmKBJkybpypUrkm62arZq1cpquhIlSmS5Dund2rX+ypUrqly5stavX5+hbFry++KLLyo6OlrvvvuuatWqJU9PTz3xxBPZ3ojtypUreuaZZyzXTadXrVo1/frrr9nGm1MHDhywXBt9q6pVq+rQoUNas2aNoqOjNXToUL3zzjvasGGDXF1dM52np6dnpol8TphMpgx3IU9KSrrt+WYm/TqlxZ/+UoO8Rks3AAAAUIyVL19eDz/8sGbOnGl1Iy97PPXUU9qwYYM6d+6ssmXL5nn5zIwdO1bvvvuuTp8+rUqVKsnf31+//fabatWqZfWXlmQ2btxYu3fvznBtb1aaNWumuLg4ubi4ZJhvhQoVJEmbNm2y3DisUaNG8vPz0/Hjx63m4+bmlqFFt1mzZtq/f3+G+daqVUtubm6qW7eukpOTtWPHDss0hw4d0qVLl3K1vQ4ePKhVq1apa9eumZbx9PTUo48+qunTp2v9+vWKjY3VL7/8kuk65ETaDevS/PTTT6pdu7blpIivr6/VteKHDx/WtWvXLO/Trt/PKoaaNWvKzc1NmzZtsgxLSkrStm3bVL9+/VzHnhdIugEAAIBibvbs2UpOTlaLFi20ZMkSHThwQIcOHdJ//vMfHTx40JIc3apevXq6cOFChseBZSan5TMTFBSkxo0b64033pAkTZgwQZMnT9b06dP166+/6pdfftGCBQss13336tVLfn5+CgsL06ZNm/Tbb7/pyy+/VGxsbKbLCA4OVlBQkMLCwrR69WodP35cmzdv1iuvvKLt27dLunktdtqN1/bs2aPevXtnaDENCAjQxo0bderUKV24cEGS9PLLL2vz5s0KDw/X7t27dfjwYX399deWG6nVqVNHoaGheuaZZ7Rlyxbt2LFDTz/9tF0t9snJyYqLi9Pp06f1yy+/aMaMGWrbtq2aNm1qdV1zegsXLtRHH32kvXv36rffftN//vMfeXp6qnr16pmuQ06cOHFCEREROnTokD777DPNmDFDI0aMsIx/8MEHNXPmTO3atUvbt2/Xs88+a9UaXbFiRXl6emrVqlU6e/asLl++nGEZJUuW1JAhQzRq1CitWrVK+/fv16BBg3Tt2jUNHDgwxzHnJZJuAAAAoJirWbOmdu3apeDgYI0ePVpNmjRRixYtNGPGDL344ouaNGlSptOWL1/ermQwt+Uz8/zzz2v+/Pk6efKknn76ac2fP18LFixQo0aN1LZtWy1cuNDS0u3m5qbVq1erYsWK6tixoxo1aqQ333wz05MJ0s1uxytXrtT999+vAQMG6M4771TPnj31+++/q1KlSpJu3sytbNmyat26tR599FGFhISoWbNmVvOZOHGijh8/rpo1a1quGW/cuLE2bNigX3/9Vffdd5/uuusujRs3Tv7+/pbpFixYIH9/f7Vt21aPP/64Bg8enOndwdPbt2+fKleurGrVqqldu3b6/PPPNXr0aP33v/9VqVKlbE7j4+OjDz/8UPfee68aN26sNWvW6Ntvv1X58uUzXYec6Nu3r65fv66WLVtq2LBhGjFihAYPHmwZP2XKFFWtWlX33XefevfurRdffNHqGeIuLi6aPn26PvjgA/n7+6tz5842l/Pmm2+qa9euevLJJ9WsWTMdOXJEP/zww231qsgLJuPWzvOwS3x8vMqUKaPLly/L29vb2eHYlJSUpJUrV6pjx45ZXosB5xu4cJs+6n+3s8PIc9RBOBP1D85GHSyebty4oWPHjlk969kZUlNTFR8fL29vb5nNtLMhfxWl+pfVd9renLBwbwEAAAAAAAowkm6gEBq4cJsGLtzm7DAAAAAAZIOkGwAAAAAAByHpBgAAAADAQUi6AQAAAABwEJJuAAAAII/d+qxmAIVTXnyXXfIgDgAAAAC6+Txos9ms06dPy9fXV25ubjKZTPkeR2pqqhITE3Xjxo1C/8gmFD5Fof4ZhqHExESdP39eZrNZbm5uuZ4XSTcAAACQR8xmswIDA3XmzBmdPn3aaXEYhqHr16/L09PTKUk/ireiVP+8vLxUrVq12zp5QNINAAAA5CE3NzdVq1ZNycnJSklJcUoMSUlJ2rhxo+6//365uro6JQYUX0Wl/pUoUUIuLi63feKApBsAAADIYyaTSa6urk5LOEqUKKHk5GR5eHgU6qQHhRP1z1rh7GAPAAAAAEAhQNINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5SIJLuWbNmKSAgQB4eHmrVqpW2bt2aZfmlS5eqbt268vDwUKNGjbRy5UrLuKSkJL388stq1KiRSpYsKX9/f/Xt21enT5+2mkdAQIBMJpPV35tvvumQ9QMAAAAAFE9OT7qXLFmiiIgIRUZGaufOnWrSpIlCQkJ07tw5m+U3b96sXr16aeDAgdq1a5fCwsIUFhamvXv3SpKuXbumnTt36tVXX9XOnTu1bNkyHTp0SI899liGeU2cOFFnzpyx/A0fPtyh6woAAAAAKF6cnnRPnTpVgwYN0oABA1S/fn3NnTtXXl5e+vjjj22Wf//99xUaGqpRo0apXr16mjRpkpo1a6aZM2dKksqUKaPo6Gh1795dderU0T333KOZM2dqx44dOnHihNW8SpcuLT8/P8tfyZIlHb6+AAAAAIDiw6lJd2Jionbs2KHg4GDLMLPZrODgYMXGxtqcJjY21qq8JIWEhGRaXpIuX74sk8kkHx8fq+Fvvvmmypcvr7vuukvvvPOOkpOTc78yAAAAAADcwsWZC79w4YJSUlJUqVIlq+GVKlXSwYMHbU4TFxdns3xcXJzN8jdu3NDLL7+sXr16ydvb2zL8ueeeU7NmzVSuXDlt3rxZo0eP1pkzZzR16lSb80lISFBCQoLlfXx8vKSb15AnJSVlv7JOkBZXQY0P/3BRao4+JxelSir4ny11EM5E/YOzUQfhTNQ/OFNxqX/2rp9Tk25HS0pKUvfu3WUYhubMmWM1LiIiwvK6cePGcnNz0zPPPKPJkyfL3d09w7wmT56sCRMmZBi+evVqeXl55X3weSg6OtrZISAbHcvK6oaA9pSXcjaNM1EH4UzUPzgbdRDORP2DMxX1+nft2jW7yjk16a5QoYJKlCihs2fPWg0/e/as/Pz8bE7j5+dnV/m0hPv333/X2rVrrVq5bWnVqpWSk5N1/Phx1alTJ8P40aNHWyXq8fHxqlq1qtq3b5/tvJ0lKSlJ0dHRevjhh+Xq6urscJCF8EU7NbNPsxyVl5SjaZyBOghnov7B2aiDcCbqH5ypuNS/tN7P2XFq0u3m5qbmzZsrJiZGYWFhkqTU1FTFxMQoPDzc5jRBQUGKiYnRyJEjLcOio6MVFBRkeZ+WcB8+fFjr1q1T+fLls41l9+7dMpvNqlixos3x7u7uNlvAXV1dC3xFKgwxFnfJMufoM0r+3+0YCsvnSh2EM1H/4GzUQTgT9Q/OVNTrn73r5vTu5REREerXr59atGihli1batq0abp69aoGDBggSerbt6+qVKmiyZMnS5JGjBihtm3basqUKerUqZMWL16s7du3a968eZJuJtxPPPGEdu7cqRUrViglJcVyvXe5cuXk5uam2NhYbdmyRQ888IBKly6t2NhYPf/88/rXv/6lsmXLOmdDAAAAAACKHKcn3T169ND58+c1btw4xcXFqWnTplq1apXlZmknTpyQ2fzPTdZbt26tqKgojR07VmPGjFHt2rW1fPlyNWzYUJJ06tQpffPNN5Kkpk2bWi1r3bp1ateundzd3bV48WKNHz9eCQkJCgwM1PPPP2/VfRwAAAAAgNvl9KRbksLDwzPtTr5+/foMw7p166Zu3brZLB8QECDDMLJcXrNmzfTTTz/lOE4AAAAAAHLCqc/pBgAAAACgKCPpBgAAAADAQUi6AQAAAABwEJJuAAAAAAAchKQbAAAAAAAHIekGAAAAAMBBSLoBAAAAAHAQkm4AAAAAAByEpBsAAAAAAAch6QYAAAAAwEFIuoHCIqrHzT8AAAAAhQZJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0UZFE9nB0BAAAAgNtA0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOIiLswMAkHPDz47936sfnBoHAAAAgKzR0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3UBBE9XD2REAAAAAyCMk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOUiCS7lmzZikgIEAeHh5q1aqVtm7dmmX5pUuXqm7duvLw8FCjRo20cuVKy7ikpCS9/PLLatSokUqWLCl/f3/17dtXp0+ftprHxYsX1adPH3l7e8vHx0cDBw7UlStXHLJ+AAAAAIDiyelJ95IlSxQREaHIyEjt3LlTTZo0UUhIiM6dO2ez/ObNm9WrVy8NHDhQu3btUlhYmMLCwrR3715J0rVr17Rz5069+uqr2rlzp5YtW6ZDhw7pscces5pPnz59tG/fPkVHR2vFihXauHGjBg8e7PD1BQAAAAAUH05PuqdOnapBgwZpwIABql+/vubOnSsvLy99/PHHNsu///77Cg0N1ahRo1SvXj1NmjRJzZo108yZMyVJZcqUUXR0tLp37646deronnvu0cyZM7Vjxw6dOHFCknTgwAGtWrVK8+fPV6tWrdSmTRvNmDFDixcvztAiDgAAAABAbrk4c+GJiYnasWOHRo8ebRlmNpsVHBys2NhYm9PExsYqIiLCalhISIiWL1+e6XIuX74sk8kkHx8fyzx8fHzUokULS5ng4GCZzWZt2bJFXbp0yTCPhIQEJSQkWN7Hx8dLutmdPSkpKdt1dYa0uApqfPiHi1LTfU4u0v9e/3LqqhqlHy5JSUlKNbv+72XB/mypg3Am6h+cjToIZ6L+wZmKS/2zd/2cmnRfuHBBKSkpqlSpktXwSpUq6eDBgzaniYuLs1k+Li7OZvkbN27o5ZdfVq9eveTt7W2ZR8WKFa3Kubi4qFy5cpnOZ/LkyZowYUKG4atXr5aXl5ftFSwgoqOjnR0CstGxrP65N0HJ7lLa6zsH6WT64dLNcXcOkqR/xhVw1EE4E/UPzkYdhDNR/+BMRb3+Xbt2za5yTk26HS0pKUndu3eXYRiaM2fObc1r9OjRVi3s8fHxqlq1qtq3b29J5guapKQkRUdH6+GHH5arq6uzw0EWwhft1Mw+zW6+Wdpf6rZQkvTLtC5qNPKrf4ZLUreF+mXazd4YlnEFFHUQzkT9g7NRB+FM1D84U3Gpf2m9n7Pj1KS7QoUKKlGihM6ePWs1/OzZs/Lz87M5jZ+fn13l0xLu33//XWvXrrVKjP38/DLcqC05OVkXL17MdLnu7u5yd3fPMNzV1bXAV6TCEGNxlyxzus8oWfrfa3NqkvVwSXJ1lTk16X8vC8fnSh2EM1H/4GzUQTgT9Q/OVNTrn73r5tQbqbm5ual58+aKiYmxDEtNTVVMTIyCgoJsThMUFGRVXrrZbSF9+bSE+/Dhw1qzZo3Kly+fYR6XLl3Sjh07LMPWrl2r1NRUtWrVKi9WDQAAAAAA53cvj4iIUL9+/dSiRQu1bNlS06ZN09WrVzVgwABJUt++fVWlShVNnjxZkjRixAi1bdtWU6ZMUadOnbR48WJt375d8+bNk3Qz4X7iiSe0c+dOrVixQikpKZbrtMuVKyc3NzfVq1dPoaGhGjRokObOnaukpCSFh4erZ8+e8vf3d86GAAAAAAAUOU5Punv06KHz589r3LhxiouLU9OmTbVq1SrLzdJOnDghs/mfBvnWrVsrKipKY8eO1ZgxY1S7dm0tX75cDRs2lCSdOnVK33zzjSSpadOmVstat26d2rVrJ0latGiRwsPD9dBDD8lsNqtr166aPn2641cYAAAAAFBsOD3plqTw8HCFh4fbHLd+/foMw7p166Zu3brZLB8QECDDMLJdZrly5RQVFZWjOAEAAAAAyAmnXtMNAAAAAEBRRtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3UBRENXD2REAAAAAsIGkGwAAAAAAByHpBgAAAADAQUi6AQAAAABwEJJuAAAAAAAchKQbAAAAAAAHIekGAAAAAMBBSLoBAAAAAHAQkm4AAAAAAByEpBsAUKwNXLjN2SEAAIAijKQbAAAAAAAHIekGAAAAAMBBSLoBAAAAAHAQkm4AAAAAAByEpBsAAAAAAAch6QYAAAAAwEFIugEAAAAAcBCSbqCA2X3ykrNDAAAAAJBHSLoBAAAAAHAQkm6giBm4cJuzQwAAAADwPyTdAAAAAAA4CEk3AAAAAAAOkquke926dXkdBwAAAAAARU6uku7Q0FDVrFlTr732mk6ePJnXMQEAAAAAUCTkKuk+deqUwsPD9cUXX6hGjRoKCQnR559/rsTExLyODwAAAACAQitXSXeFChX0/PPPa/fu3dqyZYvuvPNODR06VP7+/nruuee0Z8+evI4TQCHB3dMBAACAf9z2jdSaNWum0aNHKzw8XFeuXNHHH3+s5s2b67777tO+ffvyIkYAAAAAAAqlXCfdSUlJ+uKLL9SxY0dVr15dP/zwg2bOnKmzZ8/qyJEjql69urp165aXsQIAAAAAUKi45Gai4cOH67PPPpNhGHryySf19ttvq2HDhpbxJUuW1Lvvvit/f/88CxQAAAAAgMImV0n3/v37NWPGDD3++ONyd3e3WaZChQo8WgwAAAAAUKzlqnt5ZGSkunXrliHhTk5O1saNGyVJLi4uatu27e1HCAAAAABAIZWrpPuBBx7QxYsXMwy/fPmyHnjggdsOCgAAAACAoiBXSbdhGDKZTBmG//nnnypZsuRtBwUAAAAAQFGQo2u6H3/8cUmSyWRS//79rbqXp6Sk6Oeff1br1q3zNkIAAAAAAAqpHCXdZcqUkXSzpbt06dLy9PS0jHNzc9M999yjQYMG5W2EAAAAAAAUUjlKuhcsWCBJCggI0IsvvkhXcgAAAAAAspCrR4ZFRkbmdRwAAAAAABQ5difdzZo1U0xMjMqWLau77rrL5o3U0uzcuTNPggMAAAAAoDCzO+nu3Lmz5cZpYWFhjooHAIBcG7hwmyTpo/53OzkSAACAm+xOutN3Kad7OQAAAAAA2cvVc7oBIC+ltU4CAAAARY3dLd1ly5bN8jru9C5evJjrgIDiaPjZsZJ+cHYYAAAAAPKY3Un3tGnTHBgGAAAAAABFj91Jd79+/RwZBwAAAAAARY7dSXd8fLy8vb0tr7OSVg4AAAAAgOIsR9d0nzlzRhUrVpSPj4/N67sNw5DJZFJKSkqeBgkAAAAAQGFkd9K9du1alStXTpK0bt06hwUEAAAAAEBRYXfS3bZtW5uvAQAAAACAbXYn3bf666+/9NFHH+nAgQOSpPr162vAgAGW1nAAAAAAAIo7c24m2rhxowICAjR9+nT99ddf+uuvvzR9+nQFBgZq48aNeR0jAAAAAACFUq5auocNG6YePXpozpw5KlGihCQpJSVFQ4cO1bBhw/TLL7/kaZAAAAAAABRGuWrpPnLkiF544QVLwi1JJUqUUEREhI4cOZJnwQGwz+6Tl5wdAgAAAAAbcpV0N2vWzHItd3oHDhxQkyZNbjsoAAAAAACKAru7l//888+W188995xGjBihI0eO6J577pEk/fTTT5o1a5befPPNvI8SAAAAAIBCyO6ku2nTpjKZTDIMwzLspZdeylCud+/e6tGjR95EBwAAAABAIWZ30n3s2DFHxgEAAAAAQJFjd9JdvXp1R8YBIBtpN0tr6tQoAAAAAORErh4Zlmb//v06ceKEEhMTrYY/9thjtxUUAAAAAABFQa6S7t9++01dunTRL7/8YnWdt8lkknTzmd0AAAAAABR3uXpk2IgRIxQYGKhz587Jy8tL+/bt08aNG9WiRQutX78+j0MEAAAAAKBwylVLd2xsrNauXasKFSrIbDbLbDarTZs2mjx5sp577jnt2rUrr+MEAAAAAKDQyVVLd0pKikqXLi1JqlChgk6fPi3p5s3WDh06lHfRAQAAAABQiOUq6W7YsKH27NkjSWrVqpXefvttbdq0SRMnTlSNGjVyNK9Zs2YpICBAHh4eatWqlbZu3Zpl+aVLl6pu3bry8PBQo0aNtHLlSqvxy5YtU/v27VW+fHmZTCbt3r07wzzatWsnk8lk9ffss8/mKG4AAAAAALKTq6R77NixSk1NlSRNnDhRx44d03333aeVK1dq+vTpds9nyZIlioiIUGRkpHbu3KkmTZooJCRE586ds1l+8+bN6tWrlwYOHKhdu3YpLCxMYWFh2rt3r6XM1atX1aZNG7311ltZLnvQoEE6c+aM5e/tt9+2O26gIBt+dqyzQwAAAADwP7m6pjskJMTyulatWjp48KAuXryosmXLWu5gbo+pU6dq0KBBGjBggCRp7ty5+u677/Txxx/r//7v/zKUf//99xUaGqpRo0ZJkiZNmqTo6GjNnDlTc+fOlSQ9+eSTkqTjx49nuWwvLy/5+fnZHSsAAAAAADmVq5bu9E6ePKmTJ0+qXLlyOUq4ExMTtWPHDgUHB/8TjNms4OBgxcbG2pwmNjbWqrx08wRAZuWzsmjRIlWoUEENGzbU6NGjde3atRzPAwAAAACArOSqpTs5OVkTJkzQ9OnTdeXKFUlSqVKlNHz4cEVGRsrV1TXbeVy4cEEpKSmqVKmS1fBKlSrp4MGDNqeJi4uzWT4uLi5H8ffu3VvVq1eXv7+/fv75Z7388ss6dOiQli1bluk0CQkJSkhIsLyPj4+XJCUlJSkpKSlHy88vaXEV1Pjwj1Szq+Vzyuq1dPPzTP86q2mcwUWpGepedvGknwa4HS66eelT+n0z9Q/OwnEYzkT9gzMVl/pn7/qZDMMwcjrzIUOGaNmyZZo4caKCgoIk3WyFHj9+vMLCwjRnzpxs53H69GlVqVJFmzdvtsxDkl566SVt2LBBW7ZsyTCNm5ubPvnkE/Xq1csybPbs2ZowYYLOnj1rVfb48eMKDAzUrl271LRp0yxjWbt2rR566CEdOXJENWvWtFlm/PjxmjBhQobhUVFR8vLyynL+AAAAAICi5dq1a+rdu7cuX74sb2/vTMvlqqU7KipKixcvVocOHSzDGjdurKpVq6pXr152Jd0VKlRQiRIlMiTLZ8+ezfRaaz8/vxyVt1erVq0kKcuke/To0YqIiLC8j4+PV9WqVdW+ffssN7AzJSUlKTo6Wg8//LBdvQ/gPL9M66JGI7/K9rUkNRr5ldXrrKZxhvBFOzWzTzNJ9tfB9NMAtyN80U5J0sw+zah/cDqOw3Am6h+cqbjUv7Tez9nJVdLt7u6ugICADMMDAwPl5uZm1zzc3NzUvHlzxcTEKCwsTJKUmpqqmJgYhYeH25wmKChIMTExGjlypGVYdHS0VUt5bqQ9Vqxy5cqZlnF3d5e7u3uG4a6urgW+IhWGGIs7c2qS5TPK6rV08/NM/zqraZwhWeYMy8+uDtqaBsiN5P/dqiR9faL+wdk4DsOZqH9wpqJe/+xdt1wl3eHh4Zo0aZIWLFhgSUQTEhL0+uuvZ5ow2xIREaF+/fqpRYsWatmypaZNm6arV69a7mbet29fValSRZMnT5YkjRgxQm3bttWUKVPUqVMnLV68WNu3b9e8efMs87x48aJOnDih06dPS5IOHTok6WYruZ+fn44ePaqoqCh17NhR5cuX188//6znn39e999/vxo3bpybzQEAAAAAgE12J92PP/641fs1a9bojjvuUJMmTSRJe/bsUWJioh566CG7F96jRw+dP39e48aNU1xcnJo2bapVq1ZZbpZ24sQJmc3/3GC9devWioqK0tixYzVmzBjVrl1by5cvV8OGDS1lvvnmG0vSLkk9e/aUJEVGRmr8+PFyc3PTmjVrLAl+1apV1bVrV40dy7ONAQAAAAB5y+6ku0yZMlbvu3btavW+atWquQogPDw809bx9evXZxjWrVs3devWLdP59e/fX/379890fNWqVbVhw4achgnkqYELt0mSPup/t5MjAQAAAOBIdifdCxYscGQcAAAAAAAUObm6pjvN+fPnLddM16lTR76+vnkSFAAAAAAARYE5+yIZXb16VU899ZQqV66s+++/X/fff7/8/f01cOBAXbt2La9jBAAAAACgUMpV0h0REaENGzbo22+/1aVLl3Tp0iV9/fXX2rBhg1544YW8jhEAAAAAgEIpV93Lv/zyS33xxRdq166dZVjHjh3l6emp7t27a86cOXkVHwAAAAAAhVauWrqvXbtmeaxXehUrVqR7OQAAAAAA/5OrpDsoKEiRkZG6ceOGZdj169c1YcIEBQUF5VlwAAAAAAAUZrnqXj5t2jSFhobqjjvuUJMmTSRJe/bskYeHh3744Yc8DRAAAAAAgMIqV0l3o0aNdPjwYS1atEgHDx6UJPXq1Ut9+vSRp6dnngYIAAAAAEBhleOkOykpSXXr1tWKFSs0aNAgR8QEAAAAAECRkONrul1dXa2u5QYAAAAAALbl6kZqw4YN01tvvaXk5OS8jgcAAAAAgCIjV9d0b9u2TTExMVq9erUaNWqkkiVLWo1ftmxZngQHAAAAAEBhlquk28fHR127ds3rWAAAAAAAKFJylHSnpqbqnXfe0a+//qrExEQ9+OCDGj9+PHcsBwAAAADAhhxd0/36669rzJgxKlWqlKpUqaLp06dr2LBhjooNAAAAAIBCLUdJ96effqrZs2frhx9+0PLly/Xtt99q0aJFSk1NdVR8AAAAAAAUWjlKuk+cOKGOHTta3gcHB8tkMun06dN5HhgAAAAAAIVdjpLu5ORkeXh4WA1zdXVVUlJSngYFAAAAAEBRkKMbqRmGof79+8vd3d0y7MaNG3r22WetHhvGI8MAAAAAAMhh0t2vX78Mw/71r3/lWTAAAAAAABQlOUq6FyxY4Kg4AAAAAAAocnJ0TTcAAAAAALAfSTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAuzg4AKI6Gnx37v1c/ODUOAAAAAI5FSzcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3UJRF9XB2BAAAAECxRtINAAAAAICDkHQDAAAAAOAgJN0AABQRAxdu08CF25wdBgAASIekGwAAAAAAByHpBpApWswAAACA20PSDQAAAACAg5B0AwAAAADgICTdAAAAAAA4CEk3AAAAAAAOQtINAAAAAICDkHQDAAAAAOAgJN0AAAAAADgISTcAAAAAAA5C0g0AAAAAgIOQdAMAAAAA4CAk3QAAAAAAOAhJNwAAAAAADkLSDQAAAACAg5B0A8gXAxdu08CF25wdBgAAAJCvSLoBAAAAAHAQkm4AAAAAAByEpBsAAAAAAAdxetI9a9YsBQQEyMPDQ61atdLWrVuzLL906VLVrVtXHh4eatSokVauXGk1ftmyZWrfvr3Kly8vk8mk3bt3Z5jHjRs3NGzYMJUvX16lSpVS165ddfbs2bxcLaDA4ZpqAAAAIP85NelesmSJIiIiFBkZqZ07d6pJkyYKCQnRuXPnbJbfvHmzevXqpYEDB2rXrl0KCwtTWFiY9u7daylz9epVtWnTRm+99Vamy33++ef17bffaunSpdqwYYNOnz6txx9/PM/XDwAAAABQvDk16Z46daoGDRqkAQMGqH79+po7d668vLz08ccf2yz//vvvKzQ0VKNGjVK9evU0adIkNWvWTDNnzrSUefLJJzVu3DgFBwfbnMfly5f10UcfaerUqXrwwQfVvHlzLViwQJs3b9ZPP/3kkPUEAAAAABRPTku6ExMTtWPHDqvk2Gw2Kzg4WLGxsTaniY2NzZBMh4SEZFrelh07digpKclqPnXr1lW1atVyNB8Ajkd3eAAAABR2Ls5a8IULF5SSkqJKlSpZDa9UqZIOHjxoc5q4uDib5ePi4uxeblxcnNzc3OTj45Oj+SQkJCghIcHyPj4+XpKUlJSkpKQku5efn9LiKqjxFWepZldJ/3w2qWZXu16nTWPv9JKL9L/XLkq1msYeLkrNcf1JP036/1ktP7Pl5Gb5KN7S1zN794FFqZ7l5nsOx+E4DGei/sGZikv9s3f9TIZhGA6OxabTp0+rSpUq2rx5s4KCgizDX3rpJW3YsEFbtmzJMI2bm5s++eQT9erVyzJs9uzZmjBhQoYboR0/flyBgYHatWuXmjZtahkeFRWlAQMGWCXQktSyZUs98MADmV4LPn78eE2YMCHD8KioKHl5edm1zgAAAACAouHatWvq3bu3Ll++LG9v70zLOa2lu0KFCipRokSGZPns2bPy8/OzOY2fn1+Oymc2j8TERF26dMmqtTu7+YwePVoRERGW9/Hx8apatarat2+f5QZ2pqSkJEVHR+vhhx+Wq6urs8NBOr9M6yJJajTyK8t7e16nTWPv9FraX+q2UJIUvminJGlmn2Z2xxm+aGeOyt86Tfo6+Pznv2S6/MyWk5vlo3hLX8/t3QcWpXqWm+85HKewH4eL0nejOCrs9Q+FW3Gpf2m9n7PjtKTbzc1NzZs3V0xMjMLCwiRJqampiomJUXh4uM1pgoKCFBMTo5EjR1qGRUdHW7WUZ6d58+ZydXVVTEyMunbtKkk6dOiQTpw4keV83N3d5e7unmG4q6trga9IhSHG4sacerMrStrnYk5Nsut12jT2Ti8lS/97nfy/WzjkpC4ky5zjumNrGldX1yyXn9lycrN8FG+26ll2+8CiVM9y8z2H4xXW43BR+m4UZ4W1/qFoKOr1z951c1rSLUkRERHq16+fWrRooZYtW2ratGm6evWqBgwYIEnq27evqlSposmTJ0uSRowYobZt22rKlCnq1KmTFi9erO3bt2vevHmWeV68eFEnTpzQ6dOnJd1MqKWbLdx+fn4qU6aMBg4cqIiICJUrV07e3t4aPny4goKCdM899+TzFgAAIHMDF27TR/3vdnYYAADgNjg16e7Ro4fOnz+vcePGKS4uTk2bNtWqVassN0s7ceKEzOZ/brDeunVrRUVFaezYsRozZoxq166t5cuXq2HDhpYy33zzjSVpl6SePXtKkiIjIzV+/HhJ0nvvvSez2ayuXbsqISFBISEhmj17dj6sMQAAAACgOHFq0i1J4eHhmXYnX79+fYZh3bp1U7du3TKdX//+/dW/f/8sl+nh4aFZs2Zp1qxZOQkVAAAAAIAccdpzugEAAAAAKOpIugEAAAAAcBCSbgAAAAAAHISkGwAAAAAAByHpBgAAAADAQUi6AQAAAABwEJJuAAAAAAAchKQbAAAAAAAHIekGABQLAxduc3YIAACgGHJxdgAA8sfws2P/9+oHp8YBAAAAFCe0dAMAABRx9PQAAOch6QYAAAAAwEFIugEAAAAAcBCSbgAAAAAAHISkGwAAAAAAByHpBgAAAADAQUi6AQAAAABwEJJuAAAAAAAchKQbAAAAAAAHIekGAACA3QYu3KaBC7c5OwwAKDRIuoH8EtXD2RE41fCzYzX87Nh/BhTz7QEAAIDigaQbgOMs7e/sCAAAAACnIukGAAAAAMBBSLoBOJ1Vt3MAAACgCCHpBpCnSKABAACAf5B0AwAAAADgICTdAAoW7moOAACAIoSkGzxvEwAAAAAchKQbwO2jdRoAAACwiaQbAAAAAAAHIekGAAAAAMBBSLoBOMXuk5ecHQIAAEC+4P5JxRtJN4ACi2d+AwAAoLAj6QaQO9w8DQAAAMgWSTcAAAAAAA5C0o0ih2tm7BDVg5ZqAAAAIB+QdAMAAAAA4CAk3QAAAAAAOAhJN4DCgS7xAAAAKIRIugHYh6QXuC3cbwIAgOKJpBsQP4YBAEUPxzYAKBhIugEAAAAAcBCSbgBAoUMLHgAAKCxIuoFiaPfJS9p98pKzwwAAFADOPonl7OUDgKORdAMAAAAA4CAk3Sg2OJOec8PPjr2t6X85dTmPIgEAAAAKJ5JuAIXzUWCFMWYAAAAUOyTdQBFWGK/bLowxAwAAAJkh6QZgjRZkAAAAIM+QdAMAAAAA4CAk3SiUBi7cxo3Rijm6oQNA7nAMBYD8RdINAAAAAICDkHQDuG20OgNA0UbrOADkHkk3AAAAAAAOQtINoPDjjusAAAAooEi6AVihqzgAFF90IQeAvEfSDQAAAACAg5B0AwAAAADgICTdAOyy++Qlup6jQKI7LAoS6iMA4FYk3QAAAAAAOAhJNwAAQA7wzGoAQE6QdAPIVbdxupoDAAoLTpQgL1GXkFMk3QAAAAAAOAhJdzHCWTkUC1E9nB0BCoN09WT42bHWw6lDAAAgD5F0AwCKBavk2l4FNAHnJCoAAIUHSTcAAPYooAk4AAAo2Ei6AUfiRzoAoIAZfnasVc8Pe14XNPT2AFCYuDg7AAAA8so/ScIPjl1QVA+p9xLHLgPIQze/G7n/XuTbdwsAiiCSbgAAgCLodhPtrOcrh8wbAIoiupej4KFLNgAbMutO6vQusOyzUJBQHwGgwCkQSfesWbMUEBAgDw8PtWrVSlu3bs2y/NKlS1W3bl15eHioUaNGWrlypdV4wzA0btw4Va5cWZ6engoODtbhw4etygQEBMhkMln9vfnmm3m+brBDVo/o4ccDAAAFG8dqAMiS05PuJUuWKCIiQpGRkdq5c6eaNGmikJAQnTt3zmb5zZs3q1evXho4cKB27dqlsLAwhYWFae/evZYyb7/9tqZPn665c+dqy5YtKlmypEJCQnTjxg2reU2cOFFnzpyx/A0fPtyh61oYOexGJTwLF0BRxL4NzlCA6tzut0KcHQIAFDhOT7qnTp2qQYMGacCAAapfv77mzp0rLy8vffzxxzbLv//++woNDdWoUaNUr149TZo0Sc2aNdPMmTMl3WzlnjZtmsaOHavOnTurcePG+vTTT3X69GktX77cal6lS5eWn5+f5a9kyZKOXl3cjgL0owKFBHWmyHJ6l3IAAAA7OfVGaomJidqxY4dGjx5tGWY2mxUcHKzY2Fib08TGxioiIsJqWEhIiCWhPnbsmOLi4hQcHGwZX6ZMGbVq1UqxsbHq2bOnZfibb76pSZMmqVq1aurdu7eef/55ubjY3iQJCQlKSEiwvI+Pj5ckJSUlKSkpKWcrnk/S4kr776JUm7G6KNWq3K3jbA0PX7RTM/s0y1lAS/tL3RZa5vy/hd7y+n/jsnmdm5gzG57duNtjO/5Us6ukf+JPNbva9TptGmdPn5/LTEq//W7OwPa4bF+jsEv/PbWnzt26D8yszmRVz+yrW7e8Tj99uv1ebvYzOd2fZbVvRN6x9xh06/88s7T/zf/pj6nZ1G1H7s8z/T5ZHfczl5ffjazKS8Xru+Gw+od8qbOFXXGpf/aun8kwDMPBsWTq9OnTqlKlijZv3qygoCDL8JdeekkbNmzQli1bMkzj5uamTz75RL169bIMmz17tiZMmKCzZ89q8+bNuvfee3X69GlVrlzZUqZ79+4ymUxasuTmI16mTp2qZs2aqVy5ctq8ebNGjx6tAQMGaOrUqTZjHT9+vCZMmJBheFRUlLy8vHK9DQAAAAAAhc+1a9fUu3dvXb58Wd7e3pmWK7aPDEvfWt64cWO5ubnpmWee0eTJk+Xu7p6h/OjRo62miY+PV9WqVdW+ffssN7AzJSUlKTo6Wg8//LBcXV0zbZ0OX7RTkjIdl5PhWUp/xjv9Wfpbz9jfWs7G61+mdZEkNRr5VZ7EnKv1sYed8f8yrYtdr9Omcfb09k6zZ3p3narVXw8//LAOzuqeq2VaPuNb6kmGcdnVHxR6Oa2nt+4DM6szWdUzu+pWZvu2W8blZj+T0/1ZVvtz5B17j5u31sE8Y2c9zcv9eVbHkCy/T5l9N9LJy+9GVuWl4vXdcFj9Q77U2cKuuNS/tN7P2XFq0l2hQgWVKFFCZ8+etRp+9uxZ+fn52ZzGz88vy/Jp/8+ePWvV0n327Fk1bdo001hatWql5ORkHT9+XHXq1Mkw3t3d3WYy7urqWuArUlqMyTLbjDX5f5f2ZzYuJ8MziOoh9V5imUqWaZLTgrvltY1yNl6bU5Psinngwm36qP/d2cZs9/rkmH3xm1OT7HqdNo2zp8/PZbqm3343Z2B7XLavUdjlpp6mvU6rg7bqTFb1zL66ldm+zXrckLPj5Oqas+ca53QfnNX+HHknp8dNV1dXPbtot+V4lCu3Hk9vzvif99nUbUfuz7P8PmX63fhHZvU5/THc3mkyk9VnltVyioLC8Fu1sMnN70bH/dbMfzn5zhT1+mfvujn1Rmpubm5q3ry5YmJiLMNSU1MVExNj1d08vaCgIKvykhQdHW0pHxgYKD8/P6sy8fHx2rJlS6bzlKTdu3fLbDarYsWKt7NKRYLD7lgOAAXE7pOXnB0CUHxxl38AxYzTu5dHRESoX79+atGihVq2bKlp06bp6tWrGjBggCSpb9++qlKliiZPnixJGjFihNq2baspU6aoU6dOWrx4sbZv36558+ZJkkwmk0aOHKnXXntNtWvXVmBgoF599VX5+/srLCxM0s2bsW3ZskUPPPCASpcurdjYWD3//PP617/+pbJlyzplO+A2pB24LS0AAFDAWbVaAgCAoszpSXePHj10/vx5jRs3TnFxcWratKlWrVqlSpUqSZJOnDghs/mfBvnWrVsrKipKY8eO1ZgxY1S7dm0tX75cDRs2tJR56aWXdPXqVQ0ePFiXLl1SmzZttGrVKnl4eEi62VV88eLFGj9+vBISEhQYGKjnn38+w13Rgcyk9QYoyt3RAKfjhBqKu2JycubmIwBzdukFABQmTk+6JSk8PFzh4eE2x61fvz7DsG7duqlbt26Zzs9kMmnixImaOHGizfHNmjXTTz/9lKtYYYdi8iMBAAAAALLj1Gu6UbhwrTcKNa4hREFF3QT+wXcBQBFE0o1Cw94bH93spgYAQCFE0gkARQ5JNwCg4Cjurb7Fff2BdDiJDqCoIOkupvK0qzg/ElGA8Cgo5AfqGQCgKBi4cBuXkOaDAnEjNQBwtLQkqalTowAA5ApPMwBQiNHSjSKNrmkACjP2YcUEvcUAoEgj6QayQHcbAI62++QluqsDAFCEkXQDAJyLVj77sJ0AC3qBAChMSLqRKxzsUOiRwABwFm5A6jD0UANQEJF0AwCQCbp9Iw3JXAHGCQwABRxJN4oNWucBFBXszwAAKDxIugEA+Y+WKQA5YO8NBwvqCSmehQwUbyTdAAAUZlwfXCgU1GQwv3CXfgDFGUk3kNf48QsAgPNxPAZQQJB0w27F/Sw9AGSH1jwUR3lZ5wvS94fu4HB2HXD28pF3XJwdAOAUaWe/ey9xbhwAkIf+OTn6g1PjAAAA/6ClGwCQP+jqaVHQWvNoTXGQpf0dMtuselQUpLpVkAw/O5YeewCchqQbQLG3+60QZ4cAOAYnOgCg2OAkasFF0g0AAIDihRNSAPIRSTeyRFcsACgiSDLyBcfNwsFR3fBpZQRgC0k3AAAAii1OlABwNJJu5LnCeBOX/DjgFsbtAtyWqB60rhZUt/m50JqXt4rb8YFH66Go4ppqZIakG8hj/JAAkKaoPr8YKAz4zgAoKEi6AQAoxniUUh6gR0fRwWeJ/JC+nt36Ot179s1FB0k3AABwCLpZFm+0NAP/uO0EuhieECpKxxCSbmRQ2M6q5fm1YcVwp4Z07Pz8i9KBAEDBRNLqXIXt9xAKIEf9puS3aqHj4uwAAABFSNoPgd5LnBsHci+qB59fNv5Jxn5wahwAChan7Bs47hYKtHQDKPRoDQLgTOyDAABZIekG8gLdfAAURezbgJt4BCIyU4DqBTfGLLhIuoE8QCsHUHzx/S+6uHcDgOxwDIA9SLpRLM+K2XvzteK2XYBcKUBn+eFAxfxz5ngAK8X8+1DcFZr9AfW0wCDpBgAAAADAQUi64XRZtTrTZQdOxTV8xU6eP4KwiErfyjNw4Ta7umHTVRtFFseJ4qGwf86FPf5CjqQbsFe6nVVx7JLvbCRD+c/eZArIqfT1Kqt65uz6V5D28+z/AOfiO4jbQdINAAByjZOQAIqqIrtvo9U735F0A0AmaF3PAgdsZKaw1w0nx88+J+cK0jYrskkagNvi4uwA4Bw3Dwo/ODsMAAAAoGBIO+nWe4lz48hH/5woIi9wJFq6ASAPOfsa1IKsILVG5UZhjz83bnedHdXq58j7DRTHzxmA87DPKR5IuuEUdNsFCh9OKKBI4ikFcBTqFoD/IekGAHvw44lrFXFbqD8ACqRifmy3wm8dhyHpRr6hZRsA8l5B6zlkTyxZJeC326PCrun5YYl8kP67yV3+i6eCtG+Gc5F0FyPs7AHHKI7drvkhgUKHJBsOwv6wECpC+wOHnngtQtvJ2bh7OQAAcBrunItiKapHsbpDNlDc0dKN21bQujYCaZxdLx15h+X8QHdIOEP6Opeb+mc1jR2tNI48hv1y6rJD5gsUVQ49btJqa+Hs30fFEUk3AOQCySiQz9L9YLY6IcT12QDgeOxnbwtJNwoczr6hICqMPToKe0s7APsUtn1TQeaUbUkyAzvcbt0sqPuJ4vJbhaQbDlVQv+BAcVQcDmqwD/tm4PbwHbKtQB1nctELpjCeYHcGevvlHEk3AOAftLgAgPOwDwaKJJJuwIE4W4rCoEC1TAAAkFuF/B4PBa2lPX0sBSmuwoikGwBuE92sAAB5rpAnkCgmqKd2IekGbsGZPDgKLcpA3itoLUMA4Gjs8wofkm7ATuzgipbMPs/i+Dmnb6kvjutfkOR1AlmUPs+itC7A7eIkbjq0shYsmXwexb3OknQDAAAABRldeIFCjaQbAAqowvjsSrr6Flz2fi7F4fOjniKnqDM559BjGCcgUMiQdANAHio0N1XjBwuKGRImFDZZJfqF5lgDQBJJNwAAAFB4FZeTqMVlPVEkkXQDQCGXvvteVl35aBkBCga6KsNRCvIlSZnFVpBjdpT8+v47Yz9D7wzbXJwdAAAUVf8cYH5wahwAgOLh5nGniBxz0lq2ey9xbhxAHqClGwCcLDdn+QtSy0BuzqTTygegOCuM+8D8OO4MPzs2X1pEC+P2zw/51gsni0sFCtLvm7xE0g0ATuDsHy/5cVC93YM3XXAB4KZc7Q+X9s901O0eg5w9PQq34nhsJ+kGgHxS3K9nAgA4UQG9ERnHxuItv3o3OBtJNwAUIcXhwAUA+EduWg3z8liRm1br3Cw//Xpm9jo3imOra0FWVH/HkHQDgBPk5qCS2Q+bonqAAgDkg3xqAc/psYpLjFCUkHQDgJPlJmkm0QYA5Ln0CXhUjxwn5OmPTcWl2zAcpyjVH5JuAChI0v3AufUHS1E6+KRHS0bhUBg/p8IYM4qfvOwqnZfLv7WlObPjUUE6NqWPmZZyFCQk3QAAAAAAOIiLswMAACBNWqtEU6dGAQAoanafvMSxxQE4btuHlm4AKEDoCgcAcASOL4DzkHQDAAAAAOAgBSLpnjVrlgICAuTh4aFWrVpp69atWZZfunSp6tatKw8PDzVq1EgrV660Gm8YhsaNG6fKlSvL09NTwcHBOnz4sFWZixcvqk+fPvL29paPj48GDhyoK1eu5Pm6AQAAAPntl1OXnR0CgP9xetK9ZMkSRUREKDIyUjt37lSTJk0UEhKic+fO2Sy/efNm9erVSwMHDtSuXbsUFhamsLAw7d2711Lm7bff1vTp0zV37lxt2bJFJUuWVEhIiG7cuGEp06dPH+3bt0/R0dFasWKFNm7cqMGDBzt8fQEAAAAAxYfTk+6pU6dq0KBBGjBggOrXr6+5c+fKy8tLH3/8sc3y77//vkJDQzVq1CjVq1dPkyZNUrNmzTRz5kxJN1u5p02bprFjx6pz585q3LixPv30U50+fVrLly+XJB04cECrVq3S/Pnz1apVK7Vp00YzZszQ4sWLdfr06fxadQAAAABAEefUpDsxMVE7duxQcHCwZZjZbFZwcLBiY2NtThMbG2tVXpJCQkIs5Y8dO6a4uDirMmXKlFGrVq0sZWJjY+Xj46MWLVpYygQHB8tsNmvLli15tn4AAAAAgOLNqY8Mu3DhglJSUlSpUiWr4ZUqVdLBgwdtThMXF2ezfFxcnGV82rCsylSsWNFqvIuLi8qVK2cpc6uEhAQlJCRY3l++fPM6mYsXLyopKSnL9XSWpKQkXbt2TX/++adcXV0Vnyj9+eefkpThtZT5OHte52b6tGmcPX1+rnNx3GZpdbC4b/PCtM4FZZvd7vQ52Qeyzahnjoo5/T7Q2dusuGzzgjS9s7dZbvaBxX2bOWOdi+I2y+s8pKD6+++/Jd3sbZ0lw4lOnTplSDI2b95sNXzUqFFGy5YtbU7j6upqREVFWQ2bNWuWUbFiRcMwDGPTpk2GJOP06dNWZbp162Z0797dMAzDeP31140777wzw7x9fX2N2bNn21xuZGSkIYk//vjjjz/++OOPP/74448//ix/J0+ezDLvdWpLd4UKFVSiRAmdPXvWavjZs2fl5+dncxo/P78sy6f9P3v2rCpXrmxVpmnTppYyt96oLTk5WRcvXsx0uaNHj1ZERITlfWpqqi5evKjy5cvLZDLZsbb5Lz4+XlWrVtXJkyfl7e3t7HBQDFEH4UzUPzgbdRDORP2DMxWX+mcYhv7++2/5+/tnWc6pSbebm5uaN2+umJgYhYWFSbqZzMbExCg8PNzmNEFBQYqJidHIkSMtw6KjoxUUFCRJCgwMlJ+fn2JiYixJdnx8vLZs2aIhQ4ZY5nHp0iXt2LFDzZs3lyStXbtWqampatWqlc3luru7y93d3WqYj49PLtc8f3l7exfpyo6CjzoIZ6L+wdmog3Am6h+cqTjUvzJlymRbxqlJtyRFRESoX79+atGihVq2bKlp06bp6tWrGjBggCSpb9++qlKliiZPnixJGjFihNq2baspU6aoU6dOWrx4sbZv36558+ZJkkwmk0aOHKnXXntNtWvXVmBgoF599VX5+/tbEvt69eopNDRUgwYN0ty5c5WUlKTw8HD17Nkz27MUAAAAAADYy+lJd48ePXT+/HmNGzdOcXFxatq0qVatWmW5EdqJEydkNv9zk/XWrVsrKipKY8eO1ZgxY1S7dm0tX75cDRs2tJR56aWXdPXqVQ0ePFiXLl1SmzZttGrVKnl4eFjKLFq0SOHh4XrooYdkNpvVtWtXTZ8+Pf9WHAAAAABQ5JkMI7tbraGwSkhI0OTJkzV69OgMXeOB/EAdhDNR/+Bs1EE4E/UPzkT9s0bSDQAAAACAg5izLwIAAAAAAHKDpBsAAAAAAAch6QYAAAAAwEFIuouwWbNmKSAgQB4eHmrVqpW2bt3q7JBQBI0fP14mk8nqr27dupbxN27c0LBhw1S+fHmVKlVKXbt21dmzZ50YMQq7jRs36tFHH5W/v79MJpOWL19uNd4wDI0bN06VK1eWp6engoODdfjwYasyFy9eVJ8+feTt7S0fHx8NHDhQV65cyce1QGGVXf3r379/hn1iaGioVRnqH3Jr8uTJuvvuu1W6dGlVrFhRYWFhOnTokFUZe467J06cUKdOneTl5aWKFStq1KhRSk5Ozs9VQSFkT/1r165dhn3gs88+a1WmONY/ku4iasmSJYqIiFBkZKR27typJk2aKCQkROfOnXN2aCiCGjRooDNnzlj+fvzxR8u4559/Xt9++62WLl2qDRs26PTp03r88cedGC0Ku6tXr6pJkyaaNWuWzfFvv/22pk+frrlz52rLli0qWbKkQkJCdOPGDUuZPn36aN++fYqOjtaKFSu0ceNGDR48OL9WAYVYdvVPkkJDQ632iZ999pnVeOofcmvDhg0aNmyYfvrpJ0VHRyspKUnt27fX1atXLWWyO+6mpKSoU6dOSkxM1ObNm/XJJ59o4cKFGjdunDNWCYWIPfVPkgYNGmS1D3z77bct44pt/TNQJLVs2dIYNmyY5X1KSorh7+9vTJ482YlRoSiKjIw0mjRpYnPcpUuXDFdXV2Pp0qWWYQcOHDAkGbGxsfkUIYoyScZXX31leZ+ammr4+fkZ77zzjmXYpUuXDHd3d+Ozzz4zDMMw9u/fb0gytm3bZinz/fffGyaTyTh16lS+xY7C79b6ZxiG0a9fP6Nz586ZTkP9Q146d+6cIcnYsGGDYRj2HXdXrlxpmM1mIy4uzlJmzpw5hre3t5GQkJC/K4BC7db6ZxiG0bZtW2PEiBGZTlNc6x8t3UVQYmKiduzYoeDgYMsws9ms4OBgxcbGOjEyFFWHDx+Wv7+/atSooT59+ujEiROSpB07digpKcmqLtatW1fVqlWjLsIhjh07pri4OKs6V6ZMGbVq1cpS52JjY+Xj46MWLVpYygQHB8tsNmvLli35HjOKnvXr16tixYqqU6eOhgwZoj///NMyjvqHvHT58mVJUrly5STZd9yNjY1Vo0aNVKlSJUuZkJAQxcfHa9++ffkYPQq7W+tfmkWLFqlChQpq2LChRo8erWvXrlnGFdf65+LsAJD3Lly4oJSUFKvKLEmVKlXSwYMHnRQViqpWrVpp4cKFqlOnjs6cOaMJEybovvvu0969exUXFyc3Nzf5+PhYTVOpUiXFxcU5J2AUaWn1ytb+L21cXFycKlasaDXexcVF5cqVo17itoWGhurxxx9XYGCgjh49qjFjxqhDhw6KjY1ViRIlqH/IM6mpqRo5cqTuvfdeNWzYUJLsOu7GxcXZ3EemjQPsYav+SVLv3r1VvXp1+fv76+eff9bLL7+sQ4cOadmyZZKKb/0j6QZwWzp06GB53bhxY7Vq1UrVq1fX559/Lk9PTydGBgD5r2fPnpbXjRo1UuPGjVWzZk2tX79eDz30kBMjQ1EzbNgw7d271+o+KkB+yaz+pb8/RaNGjVS5cmU99NBDOnr0qGrWrJnfYRYYdC8vgipUqKASJUpkuFPl2bNn5efn56SoUFz4+Pjozjvv1JEjR+Tn56fExERdunTJqgx1EY6SVq+y2v/5+flluKlkcnKyLl68SL1EnqtRo4YqVKigI0eOSKL+IW+Eh4drxYoVWrdune644w7LcHuOu35+fjb3kWnjgOxkVv9sadWqlSRZ7QOLY/0j6S6C3Nzc1Lx5c8XExFiGpaamKiYmRkFBQU6MDMXBlStXdPToUVWuXFnNmzeXq6urVV08dOiQTpw4QV2EQwQGBsrPz8+qzsXHx2vLli2WOhcUFKRLly5px44dljJr165Vamqq5ccBkFf++OMP/fnnn6pcubIk6h9uj2EYCg8P11dffaW1a9cqMDDQarw9x92goCD98ssvVid/oqOj5e3trfr16+fPiqBQyq7+2bJ7925JstoHFsv65+w7ucExFi9ebLi7uxsLFy409u/fbwwePNjw8fGxulMgkBdeeOEFY/369caxY8eMTZs2GcHBwUaFChWMc+fOGYZhGM8++6xRrVo1Y+3atcb27duNoKAgIygoyMlRozD7+++/jV27dhm7du0yJBlTp041du3aZfz++++GYRjGm2++afj4+Bhff/218fPPPxudO3c2AgMDjevXr1vmERoaatx1113Gli1bjB9//NGoXbu20atXL2etEgqRrOrf33//bbz44otGbGyscezYMWPNmjVGs2bNjNq1axs3btywzIP6h9waMuT/27v7oKiqNw7g3xVYhNgFUkBQYEVBxchEs8ARUlE0K6VCUgrTkgoJMzE1LYwZSytKrawcDbIxdSZJ7UUyTVBMXgU13RCVt3SVUUHlBykv5/cHs3e4srAruZL0/cww49577rnPOfecwWfPvZdXhL29vUhPTxc6nU76qa2tlcoY+73b0NAg7rvvPjF+/HhRWFgo0tLShJOTk1i8eHFnNInuIsbG36lTp0RiYqLIy8sTJSUlYseOHcLLy0sEBQVJdfxXxx+T7i7sk08+ER4eHkKpVIoRI0aIrKyszg6JuqCIiAjh6uoqlEql6N27t4iIiBCnTp2S9tfV1YmYmBjh6OgobG1tRVhYmNDpdJ0YMd3t9u3bJwC0+pkxY4YQovnPhr311lvCxcVFWFtbi7Fjx4qioiJZHZcuXRLTpk0TdnZ2Qq1Wi5kzZ4pr1651QmvobtPe+KutrRXjx48XTk5OwsrKSnh6eorZs2e3+sKb4486ytDYAyCSk5OlMqb83i0tLRUTJ04UNjY2omfPnmL+/Pmivr7+DreG7jbGxl95ebkICgoS9957r7C2thb9+/cXCxYsEFeuXJHV818cfwohhLhz6+pERERERERE/x18ppuIiIiIiIjITJh0ExEREREREZkJk24iIiIiIiIiM2HSTURERERERGQmTLqJiIiIiIiIzIRJNxEREREREZGZMOkmIiIiIiIiMhMm3URERERERERmwqSbiIi6tA0bNmD8+PGdHQbdBZYtW4YHHnigs8MAAKxbtw7u7u7o1q0bVq1a1Wp/aWkpFAoFCgsL26wjPT0dCoUC1dXVZouzJXP0340bN6DRaJCXl3db6yUiupOYdBMRdTEVFRWYNWsW3NzcoFQq4enpiblz5+LSpUu3VI8p/6n/t/v777/x1ltvISEhQdrWkcRAo9EYTHyoc/2TMapQKLB9+3bZtvj4eOzdu/f2BPcPXL16FbGxsVi4cCHOnj2L6OjoDtUTGBgInU4He3t7AEBKSgocHBxuqY5HHnkEr732mkllzdF/SqUS8fHxWLhw4W2tl4joTmLSTUTUhZw5cwbDhw9HcXExNm/ejFOnTuGLL77A3r17ERAQgMuXL3d2iHfUd999B7VajZEjR3Z2KLfNjRs3TCrX2NiIpqamTo3hbmNnZ4cePXp0dhgoLy9HfX09Jk2aBFdXV9ja2naoHqVSiV69ekGhUNzmCOWEEGhoaDBb/0VGRiIzMxPHjx+/7XUTEd0JTLqJiLqQOXPmQKlUYvfu3QgODoaHhwcmTpyIPXv24OzZs1iyZIlU1tBKn4ODA1JSUgAAffv2BQAMHToUCoUCjzzyiFTuq6++wuDBg2FtbQ1XV1fExsZK+8rLyzF58mTY2dlBrVZj6tSpuHDhgrRfv9L81VdfwcPDA3Z2doiJiUFjYyPef/999OrVC87Ozli+fLksturqarz44otwcnKCWq3GmDFjcOTIkXb7Y8uWLXj88cfbLfP8889jypQp+PDDD+Hq6ooePXpgzpw5qK+vB9C80ldWVoZ58+ZBoVDIEpjMzEyMGjUKNjY2cHd3R1xcHP73v/9J+3U6HSZNmgQbGxv07dsX3377batVc2Pt0vfX+vXr0bdvX3Tv3t1gO/SrmDt37oSvry+sra1RXl6O69evIz4+Hr1798Y999yDhx56COnp6a2O2759O7y9vdG9e3eEhoaioqLCaAzGYj9y5AhGjx4NlUoFtVqNYcOGyW4TNtZ/Go0G7777LmbNmgWVSgUPDw+sW7dO2t/WGM3NzcW4cePQs2dP2NvbIzg4GIcPH5bVCwBhYWFQKBTS55vvgmhqakJiYiL69OkDa2trPPDAA0hLS5P261faU1NTMXr0aNja2mLIkCE4dOiQwWuk194cSUlJgZ+fHwDAy8sLCoUCpaWlbdb1559/IjAwEN27d8d9992HjIwMaV/L28vT09Mxc+ZMXLlyRRrHy5YtAwCsXbtWuvYuLi54+umnATTPjYyMDKxevVo6prS0VKp3165dGDZsGKytrZGZmdmq/4zNLcC0OeLo6IiRI0diy5Yt7fYrEdG/FZNuIqIu4vLly/jll18QExMDGxsb2b5evXohMjISW7duhRDCpPpycnIAAHv27IFOp0NqaioA4PPPP8ecOXMQHR2NY8eOYefOnejfvz+A5iRl8uTJuHz5MjIyMvDrr7/izJkziIiIkNV9+vRp7Nq1C2lpadi8eTM2bNiASZMm4a+//kJGRgZWrlyJpUuXIjs7WzomPDwclZWV2LVrF/Lz8+Hv74+xY8e2u3qfmZmJ4cOHG23rvn37cPr0aezbtw9ff/01UlJSpC8fUlNT0adPHyQmJkKn00Gn00ltmDBhAp566ikcPXoUW7duRWZmpuwLiKioKJw7dw7p6enYtm0b1q1bh8rKStm5TWnXqVOnsG3bNqSmprZ7K3VtbS1WrlyJ9evX4/jx43B2dkZsbCwOHTqELVu24OjRowgPD8eECRNQXFwsO2758uXYuHEjDh48iOrqajzzzDOyug3FYCz2yMhI9OnTB7m5ucjPz8eiRYtgZWVlcv8BQFJSEoYPH46CggLExMTglVdeQVFREYC2x+i1a9cwY8YMZGZmIisrC97e3nj00Udx7do1AM1JOQAkJydDp9NJn2+2evVqJCUl4cMPP8TRo0cRGhqKJ554QtZ3ALBkyRLEx8ejsLAQPj4+mDZtGhoaGgzWaWyOREREYM+ePVL7dDod3N3dDdYFAAsWLMD8+fNRUFCAgIAAPP744wYfJQkMDMSqVaugVqulcRwfH4+8vDzExcUhMTERRUVFSEtLQ1BQkNT+gIAAzJ49WzqmZSyLFi3CihUroNVqcf/99xuMr725BZg2RwBgxIgROHDgQJv9QET0ryaIiKhLyMrKEgDE999/b3D/Rx99JACICxcuCCGEwbL29vYiOTlZCCFESUmJACAKCgpkZdzc3MSSJUsMnmP37t3CwsJClJeXS9uOHz8uAIicnBwhhBAJCQnC1tZWXL16VSoTGhoqNBqNaGxslLYNGDBAvPfee0IIIQ4cOCDUarX4+++/Zefr16+f+PLLLw3GUlVVJQCI/fv3y7YnJCSIIUOGSJ9nzJghPD09RUNDg7QtPDxcRERESJ89PT3Fxx9/LKvnhRdeENHR0bJtBw4cEN26dRN1dXVCq9UKACI3N1faX1xcLABIdZnSroSEBGFlZSUqKysNtlMvOTlZABCFhYXStrKyMmFhYSHOnj0rKzt27FixePFi2XFZWVnSfn3s2dnZbcZgSuwqlUqkpKQYjNdY/wnR3O/PPvustL+pqUk4OzuLzz//XAjR9hi9WWNjo1CpVOKHH36Qthka/zePDTc3N7F8+XJZmQcffFDExMTIzr9+/Xppv368a7Vag7GYMkcKCgoEAFFSUtJmm/TnXrFihbStvr5e9OnTR6xcuVIIIcS+ffsEAFFVVSWEaL7W9vb2snq2bdsm1Gq1bD62FBwcLObOnSvbpq93+/btsu23OrdMmSN6q1evFhqNps3+ICL6N7O8kwk+ERGZnzBxJbsjKisrce7cOYwdO9bgfq1WC3d3d9lqmK+vLxwcHKDVavHggw8CaL69V6VSSWVcXFxgYWGBbt26ybbpV7yOHDmCmpqaVs+L1tXV4fTp0wZjqaurA4A2b8duafDgwbCwsJA+u7q64tixY+0ec+TIERw9ehSbNm2Stgkh0NTUhJKSEpw8eRKWlpbw9/eX9vfv3x+Ojo6yOkxpl6enJ5ycnIy2Q6lUylYcjx07hsbGRvj4+MjKXb9+XXZOS0tL6doAwMCBA6VrNmLECIMxmBL766+/jhdffBHffPMNQkJCEB4ejn79+knHt9d/gwYNAgBZexQKBXr16mVwJbSlCxcuYOnSpUhPT0dlZSUaGxtRW1uL8vLydo9r6erVqzh37lyr9wGMHDmy1WMNLWN0dXUF0DxXBg4c2KpeU+eIqQICAqR/W1paYvjw4dBqtSYfP27cOHh6esLLywsTJkzAhAkTEBYWZtJz5KbcRdLe3CoqKjI6R/RsbGxQW1trSpOIiP51mHQTEXUR/fv3h0KhgFarRVhYWKv9Wq0Wjo6OUuKkUChaJegtn7U05Obb1jtKf4uxnkKhMLhN/yKwmpoauLq6yp5F1mvrbcw9evSAQqFAVVVVh+Ix9hKympoavPTSS4iLi2u1z8PDAydPnjR6XlPbdc899xitC2i+Pi2fOa+pqYGFhQXy8/NliQ/Q/NKwW3FzDKbEvmzZMkyfPh0//fQTdu3ahYSEBGzZsgVhYWFG+0+vI9dmxowZuHTpElavXg1PT09YW1sjICDAbC+Aaxmjvv/N9RK7202lUuHw4cNIT0/H7t278fbbb2PZsmXIzc01+qZzU8ZlR66fIZcvXzbpiycion8jJt1ERF1Ejx49MG7cOKxduxbz5s2TJcjnz5/Hpk2bEBUVJSUFTk5O0vPJAFBcXCxbSVIqlQCa34Ktp1KpoNFosHfvXowePbpVDIMGDUJFRQUqKiqklbwTJ06guroavr6+HW6bv78/zp8/D0tLS+mlV8YolUr4+vrixIkT//jvdCuVSlk/6GM6ceKE9Dz7zQYMGICGhgYUFBRg2LBhAJqfi275JUBH2nUrhg4disbGRlRWVmLUqFFtlmtoaEBeXp60ql1UVITq6mpptdkQU2P38fGBj48P5s2bh2nTpiE5ORlhYWFG+88UhsYoABw8eBBr167Fo48+CqD5z+hdvHhRVsbKyqrVcS2p1Wq4ubnh4MGDCA4OltWt76eOuN1zJCsrS3oGu6GhAfn5+a2ei9czNI6B5hXykJAQhISEICEhAQ4ODvjtt9/w5JNPtnnM7WDKHNH7448/MHToULPEQURkbnyRGhFRF/Lpp5/i+vXrCA0Nxf79+1FRUYG0tDSMGzcOvXv3lr0RfMyYMfj0009RUFCAvLw8vPzyy7JVKWdnZ9jY2CAtLQ0XLlzAlStXADSvXiYlJWHNmjUoLi7G4cOH8cknnwAAQkJC4Ofnh8jISBw+fBg5OTmIiopCcHCwSbeitiUkJAQBAQGYMmUKdu/ejdLSUvz+++9YsmSJ7G3YNwsNDUVmZmaHz6un0Wiwf/9+nD17VkreFi5ciN9//x2xsbEoLCxEcXExduzYISU8AwcOREhICKKjo5GTk4OCggJER0fLVqM72i5T+fj4IDIyElFRUUhNTUVJSQlycnLw3nvv4aeffpLKWVlZ4dVXX0V2djby8/Px/PPP4+GHH243uTQWe11dHWJjY5Geno6ysjIcPHgQubm5UiJvrP9M0dYY9fb2xjfffAOtVovs7GxERka2uktD/+XR+fPn27wbYsGCBVi5ciW2bt2KoqIiLFq0CIWFhZg7d67JMd7sds+Rzz77DN9//z3+/PNPzJkzB1VVVZg1a5bBshqNBjU1Ndi7dy8uXryI2tpa/Pjjj1izZg0KCwtRVlaGjRs3oqmpCQMGDJCOyc7ORmlpKS5evHhbV/BNmSN6Bw4c+MdfnhERdRYm3UREXYi3tzfy8vLg5eWFqVOnol+/foiOjsbo0aNx6NAh3HvvvVLZpKQkuLu7Y9SoUZg+fTri4+Nlz3FaWlpizZo1+PLLL+Hm5obJkycDaL51d9WqVVi7di0GDx6Mxx57THqbs0KhwI4dO+Do6IigoCCEhITAy8sLW7du/UftUigU+PnnnxEUFISZM2fCx8cHzzzzDMrKyuDi4tLmcS+88AJ+/vlnKRnrqMTERJSWlqJfv37SLa73338/MjIycPLkSYwaNQpDhw7F22+/DTc3N+m4jRs3wsXFBUFBQQgLC8Ps2bOhUqmk58w72q5bkZycjKioKMyfPx8DBgzAlClTkJubK7uF29bWFgsXLsT06dMxcuRI2NnZGb1mxmK3sLDApUuXEBUVBR8fH0ydOhUTJ07EO++8Y3L/GdPWGN2wYQOqqqrg7++P5557DnFxcXB2dpYdm5SUhF9//RXu7u5trqDGxcXh9ddfx/z58+Hn54e0tDTs3LkT3t7eJsd4s9s9R1asWIEVK1ZgyJAhyMzMxM6dO9GzZ0+DZQMDA/Hyyy8jIiICTk5OeP/99+Hg4IDU1FSMGTMGgwYNwhdffIHNmzdj8ODBAID4+HhYWFjA19cXTk5Ot/RcvCmMzREAOHToEK5cuSL9KTMioruNQpjzjTtERESdLDw8HP7+/li8eHFnh4K//voL7u7u2LNnT5svo7vTUlJS8Nprr6G6urqzQyEyOEciIiIwZMgQvPnmm50cHRFRx/CZbiIi6tI++OAD/PDDD51y7t9++w01NTXw8/ODTqfDG2+8AY1GIz2DS/RfZ2yO3LhxA35+fpg3b14nR0pE1HFMuomIqEvTaDR49dVXO+Xc9fX1ePPNN3HmzBmoVCoEBgZi06ZNrd7oTPRfZWyOKJVKLF26tJOjJCL6Z3h7OREREREREZGZ8EVqRERERERERGbCpJuIiIiIiIjITJh0ExEREREREZkJk24iIiIiIiIiM2HSTURERERERGQmTLqJiIiIiIiIzIRJNxEREREREZGZMOkmIiIiIiIiMhMm3URERERERERm8n98E4ZqQDR5TQAAAABJRU5ErkJggg==", | |
| "text/plain": [ | |
| "<Figure size 1000x600 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "import random\n", | |
| "import numpy as np\n", | |
| "from qiskit import transpile\n", | |
| "from qiskit.circuit.random import random_circuit\n", | |
| "from qiskit_ibm_runtime.fake_provider import FakeKawasaki as FakeBackend\n", | |
| "from qiskit_aer import AerSimulator\n", | |
| "from qiskit_aer.noise import NoiseModel\n", | |
| "from sklearn.mixture import GaussianMixture\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "\n", | |
| "# Build noise model from backend properties\n", | |
| "backend = FakeBackend()\n", | |
| "noise_model = NoiseModel.from_backend(backend)\n", | |
| "\n", | |
| "# Get coupling map from backend\n", | |
| "coupling_map = backend.configuration().coupling_map\n", | |
| "\n", | |
| "# Get basis gates from noise model\n", | |
| "basis_gates = noise_model.basis_gates\n", | |
| "\n", | |
| "\n", | |
| "# Function to generate and simulate a random quantum circuit\n", | |
| "def simulate_random_circuit(num_qubits, num_layers, shots):\n", | |
| " simulator = AerSimulator(\n", | |
| " noise_model=noise_model, basis_gates=basis_gates, coupling_map=coupling_map\n", | |
| " )\n", | |
| " circuit = random_circuit(num_qubits, num_layers, measure=True)\n", | |
| " compiled_circuit = transpile(circuit, backend)\n", | |
| " job = simulator.run(compiled_circuit, shots=shots)\n", | |
| " result = job.result()\n", | |
| " counts = result.get_counts()\n", | |
| " return counts\n", | |
| "\n", | |
| "\n", | |
| "# Function to convert counts to probabilities\n", | |
| "def counts_to_probabilities(counts):\n", | |
| " frequencies = np.array(list(counts.values()))\n", | |
| " probabilities = frequencies / frequencies.sum()\n", | |
| " bitstrings = list(counts.keys())\n", | |
| " int_values = np.array(\n", | |
| " [int(bs, 2) for bs in bitstrings]\n", | |
| " ) # Convert bitstrings to integers\n", | |
| " return int_values, probabilities\n", | |
| "\n", | |
| "\n", | |
| "# Function to fit a Gaussian Mixture Model to sampled data\n", | |
| "def fit_gmm(int_values, probabilities, num_samples, num_components):\n", | |
| " sampled_data = np.random.choice(\n", | |
| " int_values, size=num_samples, p=probabilities\n", | |
| " ).reshape(-1, 1)\n", | |
| " gmm = GaussianMixture(n_components=num_components)\n", | |
| " gmm.fit(sampled_data)\n", | |
| " return gmm\n", | |
| "\n", | |
| "\n", | |
| "# Function to evaluate the GMM probabilities for discrete outcomes\n", | |
| "def evaluate_gmm(gmm, int_values, num_qubits):\n", | |
| " full_range = np.arange(2**num_qubits).reshape(-1, 1)\n", | |
| " full_pdf = np.exp(gmm.score_samples(full_range))\n", | |
| " gmm_probabilities = (\n", | |
| " full_pdf / full_pdf.sum()\n", | |
| " ) # Normalize GMM density over discrete outcomes\n", | |
| " gmm_probs_for_observed = gmm_probabilities[int_values]\n", | |
| " return gmm_probabilities, gmm_probs_for_observed\n", | |
| "\n", | |
| "\n", | |
| "# Function to compute space savings\n", | |
| "def compute_space_savings(num_qubits, gmm_params, counts):\n", | |
| " # Space required to store original counts (bitstrings and frequencies)\n", | |
| " original_space = len(counts) * (\n", | |
| " num_qubits + np.ceil(np.log2(max(counts.values())))\n", | |
| " ) # Bits for each bitstring and count\n", | |
| "\n", | |
| " # Space required to store GMM parameters\n", | |
| " num_components = len(gmm_params[\"weights\"])\n", | |
| " gmm_space = num_components * (\n", | |
| " 1 + num_qubits + 1\n", | |
| " ) # Weights, means, and 1D variances (or covariances in full form)\n", | |
| "\n", | |
| " # Compute space savings\n", | |
| " space_saved = original_space - gmm_space\n", | |
| " savings_percentage = (space_saved / original_space) * 100\n", | |
| "\n", | |
| " return {\n", | |
| " \"original_space\": original_space,\n", | |
| " \"gmm_space\": gmm_space,\n", | |
| " \"space_saved\": space_saved,\n", | |
| " \"savings_percentage\": savings_percentage,\n", | |
| " }\n", | |
| "\n", | |
| "\n", | |
| "# Function to plot original and recreated distributions\n", | |
| "def plot_distributions(int_values, probabilities, gmm_probabilities, num_qubits):\n", | |
| " plt.figure(figsize=(10, 6))\n", | |
| " full_range = np.arange(2**num_qubits)\n", | |
| "\n", | |
| " # Original histogram\n", | |
| " plt.bar(\n", | |
| " int_values, probabilities, width=0.4, label=\"Original Distribution\", alpha=0.7\n", | |
| " )\n", | |
| "\n", | |
| " # Recreated GMM histogram\n", | |
| " plt.bar(\n", | |
| " full_range,\n", | |
| " gmm_probabilities,\n", | |
| " width=0.4,\n", | |
| " label=\"GMM Recreated Distribution\",\n", | |
| " alpha=0.7,\n", | |
| " )\n", | |
| "\n", | |
| " # Add labels and legend\n", | |
| " plt.xlabel(\"Outcome (Integer representation of bitstring)\")\n", | |
| " plt.ylabel(\"Probability\")\n", | |
| " plt.title(\"Comparison of Original and GMM-Recreated Distributions\")\n", | |
| " plt.legend()\n", | |
| " plt.grid(True)\n", | |
| " plt.tight_layout()\n", | |
| " plt.show()\n", | |
| "\n", | |
| "\n", | |
| "# Main pipeline\n", | |
| "def main(\n", | |
| " num_qubits=5, num_layers=5, num_components=3, shots=10000, num_samples_for_fit=10000\n", | |
| "):\n", | |
| " # Simulate a quantum circuit\n", | |
| " counts = simulate_random_circuit(num_qubits, num_layers, shots)\n", | |
| " if not counts:\n", | |
| " print(\"No results were returned from the quantum simulation.\")\n", | |
| " return\n", | |
| "\n", | |
| " # Convert counts to probabilities\n", | |
| " int_values, probabilities = counts_to_probabilities(counts)\n", | |
| "\n", | |
| " # Fit a GMM\n", | |
| " gmm = fit_gmm(int_values, probabilities, num_samples_for_fit, num_components)\n", | |
| "\n", | |
| " # Evaluate GMM probabilities for all outcomes\n", | |
| " gmm_probabilities, gmm_probs_for_observed = evaluate_gmm(\n", | |
| " gmm, int_values, num_qubits\n", | |
| " )\n", | |
| "\n", | |
| " # Compute errors\n", | |
| " errors = np.abs(probabilities - gmm_probs_for_observed)\n", | |
| " mean_error = np.mean(errors)\n", | |
| " median_error = np.median(errors)\n", | |
| " print(\"Mean absolute error:\", mean_error)\n", | |
| " print(\"Median absolute error:\", median_error)\n", | |
| "\n", | |
| " # Output GMM parameters\n", | |
| " gmm_params = {\n", | |
| " \"weights\": gmm.weights_,\n", | |
| " \"means\": gmm.means_,\n", | |
| " \"covariances\": gmm.covariances_,\n", | |
| " \"precisions\": gmm.precisions_,\n", | |
| " }\n", | |
| " print(\"GMM parameters:\", gmm_params)\n", | |
| "\n", | |
| " # Compute the space savings\n", | |
| " space_data = compute_space_savings(num_qubits, gmm_params, counts)\n", | |
| " print(f\"Space saved: {space_data[\"savings_percentage\"]:.2f}%\")\n", | |
| "\n", | |
| " # Plot distributions\n", | |
| " plot_distributions(int_values, probabilities, gmm_probabilities, num_qubits)\n", | |
| "\n", | |
| "\n", | |
| "if __name__ == \"__main__\":\n", | |
| " rnd = random.randint(2, 10)\n", | |
| " main(\n", | |
| " num_qubits=rnd,\n", | |
| " num_layers=rnd,\n", | |
| " num_components=rnd,\n", | |
| " shots=10000,\n", | |
| " num_samples_for_fit=10000,\n", | |
| " )" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "/Users/churchill/GitHub/qiskit-ionq/.venv/lib/python3.12/site-packages/sklearn/base.py:1389: ConvergenceWarning: Number of distinct clusters (1) found smaller than n_clusters (3). Possibly due to duplicate points in X.\n", | |
| " return fit_method(estimator, *args, **kwargs)\n", | |
| "/Users/churchill/GitHub/qiskit-ionq/.venv/lib/python3.12/site-packages/sklearn/base.py:1389: ConvergenceWarning: Number of distinct clusters (2) found smaller than n_clusters (7). Possibly due to duplicate points in X.\n", | |
| " return fit_method(estimator, *args, **kwargs)\n", | |
| "/Users/churchill/GitHub/qiskit-ionq/.venv/lib/python3.12/site-packages/sklearn/base.py:1389: ConvergenceWarning: Number of distinct clusters (4) found smaller than n_clusters (7). Possibly due to duplicate points in X.\n", | |
| " return fit_method(estimator, *args, **kwargs)\n" | |
| ] | |
| }, | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Job with 3 qubits and 7525 shots saved to database.\n", | |
| "Job with 6 qubits and 6045 shots saved to database.\n", | |
| "Job with 7 qubits and 15023 shots saved to database.\n", | |
| "Job with 7 qubits and 16680 shots saved to database.\n", | |
| "Job with 5 qubits and 11363 shots saved to database.\n", | |
| "\n", | |
| "Retrieved Job 4:\n", | |
| "{'num_qubits': 7, 'shots': 16680, 'gmm_weights': array([4.73321343e-01, 3.60911271e-02, 4.56115108e-01, 3.44724221e-02,\n", | |
| " 3.99360800e-19, 3.99360800e-19, 3.99360800e-19]), 'gmm_means': array([[ 2.],\n", | |
| " [34.],\n", | |
| " [ 0.],\n", | |
| " [32.],\n", | |
| " [ 0.],\n", | |
| " [ 0.],\n", | |
| " [ 0.]]), 'gmm_covariances': array([[1.e-06],\n", | |
| " [1.e-06],\n", | |
| " [1.e-06],\n", | |
| " [1.e-06],\n", | |
| " [1.e-06],\n", | |
| " [1.e-06],\n", | |
| " [1.e-06]])}\n", | |
| "Reconstructed Histogram for Job 4:\n", | |
| "[7607 0 7895 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 574 0 602 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", | |
| " 0 0]\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "import sqlite3\n", | |
| "import numpy as np\n", | |
| "from qiskit import transpile\n", | |
| "from qiskit.circuit.random import random_circuit\n", | |
| "from qiskit_aer import AerSimulator\n", | |
| "from sklearn.mixture import GaussianMixture\n", | |
| "\n", | |
| "\n", | |
| "# Database management\n", | |
| "class QuantumJobDB:\n", | |
| " def __init__(self, db_name=\"file::memory:?cache=shared\"):\n", | |
| " self.connection = sqlite3.connect(db_name)\n", | |
| " self._initialize_db()\n", | |
| "\n", | |
| " def _initialize_db(self):\n", | |
| " with self.connection:\n", | |
| " self.connection.execute(\n", | |
| " \"\"\"\n", | |
| " CREATE TABLE IF NOT EXISTS quantum_jobs (\n", | |
| " id INTEGER PRIMARY KEY AUTOINCREMENT,\n", | |
| " num_qubits INTEGER,\n", | |
| " shots INTEGER,\n", | |
| " gmm_weights TEXT,\n", | |
| " gmm_means TEXT,\n", | |
| " gmm_covariances TEXT\n", | |
| " )\n", | |
| " \"\"\"\n", | |
| " )\n", | |
| "\n", | |
| " def save_job(self, num_qubits, shots, gmm):\n", | |
| " weights = \",\".join(map(str, gmm.weights_))\n", | |
| " means = \",\".join(map(str, gmm.means_.flatten()))\n", | |
| " covariances = \",\".join(map(str, gmm.covariances_.flatten()))\n", | |
| "\n", | |
| " with self.connection:\n", | |
| " self.connection.execute(\n", | |
| " \"\"\"\n", | |
| " INSERT INTO quantum_jobs (num_qubits, shots, gmm_weights, gmm_means, gmm_covariances)\n", | |
| " VALUES (?, ?, ?, ?, ?)\n", | |
| " \"\"\",\n", | |
| " (num_qubits, shots, weights, means, covariances),\n", | |
| " )\n", | |
| "\n", | |
| " def get_job(self, job_id):\n", | |
| " with self.connection:\n", | |
| " job = self.connection.execute(\n", | |
| " \"\"\"\n", | |
| " SELECT num_qubits, shots, gmm_weights, gmm_means, gmm_covariances\n", | |
| " FROM quantum_jobs WHERE id = ?\n", | |
| " \"\"\",\n", | |
| " (job_id,),\n", | |
| " ).fetchone()\n", | |
| " if job:\n", | |
| " num_qubits, shots, weights, means, covariances = job\n", | |
| " return {\n", | |
| " \"num_qubits\": num_qubits,\n", | |
| " \"shots\": shots,\n", | |
| " \"gmm_weights\": np.array(list(map(float, weights.split(\",\")))),\n", | |
| " \"gmm_means\": np.array(list(map(float, means.split(\",\")))).reshape(\n", | |
| " -1, 1\n", | |
| " ),\n", | |
| " \"gmm_covariances\": np.array(\n", | |
| " list(map(float, covariances.split(\",\")))\n", | |
| " ).reshape(-1, 1),\n", | |
| " }\n", | |
| " return None\n", | |
| "\n", | |
| "\n", | |
| "# Quantum simulation and GMM fitting\n", | |
| "def simulate_random_circuit(num_qubits, num_layers, shots):\n", | |
| " simulator = AerSimulator()\n", | |
| " circuit = random_circuit(num_qubits, num_layers, measure=True)\n", | |
| " compiled_circuit = transpile(circuit, simulator)\n", | |
| " job = simulator.run(compiled_circuit, shots=shots)\n", | |
| " result = job.result()\n", | |
| " return result.get_counts()\n", | |
| "\n", | |
| "\n", | |
| "def counts_to_probabilities(counts):\n", | |
| " frequencies = np.array(list(counts.values()))\n", | |
| " probabilities = frequencies / frequencies.sum()\n", | |
| " bitstrings = list(counts.keys())\n", | |
| " int_values = np.array([int(bs, 2) for bs in bitstrings])\n", | |
| " return int_values, probabilities\n", | |
| "\n", | |
| "\n", | |
| "def fit_gmm(int_values, probabilities, num_samples, num_components):\n", | |
| " sampled_data = np.random.choice(\n", | |
| " int_values, size=num_samples, p=probabilities\n", | |
| " ).reshape(-1, 1)\n", | |
| " gmm = GaussianMixture(n_components=num_components)\n", | |
| " gmm.fit(sampled_data)\n", | |
| " return gmm\n", | |
| "\n", | |
| "\n", | |
| "def reconstruct_histogram_from_gmm(num_qubits, shots, gmm_params):\n", | |
| " full_range = np.arange(2**num_qubits).reshape(-1, 1)\n", | |
| " gmm = GaussianMixture(\n", | |
| " n_components=len(gmm_params[\"gmm_weights\"]), covariance_type=\"diag\"\n", | |
| " )\n", | |
| " gmm.weights_ = gmm_params[\"gmm_weights\"]\n", | |
| " gmm.means_ = gmm_params[\"gmm_means\"]\n", | |
| " gmm.covariances_ = gmm_params[\"gmm_covariances\"]\n", | |
| " gmm.precisions_cholesky_ = 1.0 / np.sqrt(gmm.covariances_)\n", | |
| "\n", | |
| " pdf = np.exp(gmm.score_samples(full_range))\n", | |
| " probabilities = pdf / pdf.sum()\n", | |
| " histogram = (probabilities * shots).astype(int)\n", | |
| " return histogram\n", | |
| "\n", | |
| "\n", | |
| "# Usage example\n", | |
| "import random\n", | |
| "\n", | |
| "\n", | |
| "def main():\n", | |
| " db = QuantumJobDB()\n", | |
| "\n", | |
| " # Generate multiple quantum jobs with random parameters\n", | |
| " num_jobs = 5\n", | |
| " for _ in range(num_jobs):\n", | |
| " # Randomize quantum job parameters\n", | |
| " num_qubits = random.randint(2, 7)\n", | |
| " num_layers = random.randint(2, 7)\n", | |
| " shots = random.randint(5000, 20000)\n", | |
| " # num_components = random.randint(2, 5)\n", | |
| "\n", | |
| " # Simulate quantum circuit\n", | |
| " counts = simulate_random_circuit(num_qubits, num_layers, shots)\n", | |
| " int_values, probabilities = counts_to_probabilities(counts)\n", | |
| "\n", | |
| " # Fit a GMM\n", | |
| " gmm = fit_gmm(\n", | |
| " int_values, probabilities, num_samples=shots, num_components=num_qubits\n", | |
| " )\n", | |
| "\n", | |
| " # Save the job to the database\n", | |
| " db.save_job(num_qubits, shots, gmm)\n", | |
| " print(f\"Job with {num_qubits} qubits and {shots} shots saved to database.\")\n", | |
| "\n", | |
| " # Fetch and reconstruct a random job from the database\n", | |
| " job_id = random.randint(1, num_jobs)\n", | |
| " job_data = db.get_job(job_id)\n", | |
| " if job_data:\n", | |
| " print(f\"\\nRetrieved Job {job_id}:\")\n", | |
| " print(job_data)\n", | |
| "\n", | |
| " # Reconstruct histogram from GMM\n", | |
| " histogram = reconstruct_histogram_from_gmm(\n", | |
| " num_qubits=job_data[\"num_qubits\"],\n", | |
| " shots=job_data[\"shots\"],\n", | |
| " gmm_params=job_data,\n", | |
| " )\n", | |
| " print(f\"Reconstructed Histogram for Job {job_id}:\")\n", | |
| " print(histogram)\n", | |
| " else:\n", | |
| " print(f\"No job found with ID {job_id}.\")\n", | |
| "\n", | |
| "\n", | |
| "if __name__ == \"__main__\":\n", | |
| " main()" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": ".venv", | |
| "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.7" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment