{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Particle Systems",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "1e24c9aa62954c4099b28c89d46d55ca": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "VBoxModel",
          "state": {
            "_view_name": "VBoxView",
            "_dom_classes": [
              "widget-interact"
            ],
            "_model_name": "VBoxModel",
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.5.0",
            "box_style": "",
            "layout": "IPY_MODEL_ebd017bc9dca41099c4fa560e9f68922",
            "_model_module": "@jupyter-widgets/controls",
            "children": [
              "IPY_MODEL_d34c46ec2e8f4170af80954555f8a2b7",
              "IPY_MODEL_1dbff5ccf90247b69e19cb78cbe86e25"
            ]
          }
        },
        "ebd017bc9dca41099c4fa560e9f68922": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "d34c46ec2e8f4170af80954555f8a2b7": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "IntSliderModel",
          "state": {
            "_view_name": "IntSliderView",
            "style": "IPY_MODEL_2d5a3a03f7cb421cbb0cae446cbdfdb1",
            "_dom_classes": [],
            "description": "t",
            "step": 1,
            "_model_name": "IntSliderModel",
            "orientation": "horizontal",
            "max": 19,
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "value": 13,
            "_view_count": null,
            "disabled": false,
            "_view_module_version": "1.5.0",
            "min": 0,
            "continuous_update": true,
            "readout_format": "d",
            "description_tooltip": null,
            "readout": true,
            "_model_module": "@jupyter-widgets/controls",
            "layout": "IPY_MODEL_f17ea6f4ce5245009f234130c9273294"
          }
        },
        "1dbff5ccf90247b69e19cb78cbe86e25": {
          "model_module": "@jupyter-widgets/output",
          "model_name": "OutputModel",
          "state": {
            "_view_name": "OutputView",
            "msg_id": "",
            "_dom_classes": [],
            "_model_name": "OutputModel",
            "outputs": [
              {
                "output_type": "display_data",
                "metadata": {
                  "tags": [],
                  "needs_background": "light"
                },
                "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de3Add5Xnv+deSbYcKYBlM84mku2BBLDBcZAqgWXX4A0TFA+VjJMUm4zsEYkZ5bEh9rI7VaG8RaWgXLMzU7Nr5+EkBofRJJphGPJYF/E4DGBKVZ6EjZSHQdkETPyQAya2HEDGSvS4Z/+498p9r/rX/ev3r7vPp8pl3e6+fc/t++vT53deP2JmCIIgCOmnkLQAgiAIQjiIQhcEQcgIotAFQRAygih0QRCEjCAKXRAEISOIQhcEQcgIrgqdiB4hojeJ6KeK/URE9xLRISI6SEQfDV9MQRAEwQ0dC/3vAHQ77L8awMWVf30AHgwuliAIguAVV4XOzIMATjscci2Av+cyzwF4NxFdEJaAgiAIgh4NIZzjQgCjltfHK9t+VX8gEfWhbMXjvPPO6/zgBz8Ywsenh8nTZ3D26JtAyVKdWyAsWPpeNC1sSU6wOmzlJAI0q4qpWMS7Vi+LRjhhDpOnz+DtN8ZQmpxGoakB8y9sM2o8CeEyPDx8ipkX2+0LQ6Frw8y7AOwCgK6uLh4aGorz4xNn37JbMVE6WbuxBDSXFqN76OFkhLLBVk4vHSJmgOZTizFx7BSaOxZh5bYetPesAQCMDgxiZOuA7T7BO6MDg3ix70HMTL5T3jAJFH89D5d97Xa5rhmFiI6q9oWh0N8A0G55fVFlm1DHxLFTnrYnRRjyTBw9Ofv/i33nwiov9j2ImbPvzO4b2rgdYwdexeqdfYE/M4+MbB2YvZ5VZs6+g5GtA6LQc0gYaYt7APxZJdvlYwB+y8xz3C0C0NyxyNP2sBgdGMS+ZbfiycL12LfsVowODDoeH7Y8VQVjp3zAwOGH9rnKJNiTFiPBDa9jVLBHJ23xHwE8C+ADRHSciDYR0W1EdFvlkL0AXgdwCMDXAdwRmbQpZ+W2HhQXzKvZVlwwDyu39UT2mdUp+cTRkwDzrMXsdMPYyRmUiWOn1EqGy5am4J2kjIQw8TNGBXt0slxuYuYLmLmRmS9i5t3M/BAzP1TZz8z8X5j5fcz8EWbOl2PcA+09a3DZrtvRvHQxQITmpYtx2a5ofZ1OU3ItOUOiuWORo5JJm0VpCkkYCWHjZ4wK9sQaFBXKyjJO36bblFwVpKzKaPV5+8WqYIY2brcNsKbJojSJ6u+U5kBzVtxGJiAKPeM0dyyaDVDWb5/NkLAEKasBzPaeNTi4eXdgZU7FQs0sZOzAqzj80L4apZ42i9I04jYSwsZpjArekF4uGcdpSu401R0dGMTk2Higzy4umIfO/rtqlM3qnX3oenRLrG4nwWyy4DYyBUpqCbo85qEnhcqt8mThevtiISKl1VRPsWU+OjZ+Cif2DmPi6ElQsQCeKaF56eLUTf2F5JD6BH2IaJiZu2z3iULPL/uW3Wo/1V1aLgpSVYZ2PbZFbjZBSAgnhS4ulxzjNNVV+S8b21pFmQuCoYhCzzFOaZQqZX/pjk0JSStkASkgihbJcsk5qgyJLKTDCWbhllUlBEd86IIgxIJTzKb7iDnN6UxHfOiCICSOFBBFjyh0IXeIHzcZstB3xnREoQu5QhpBJYcUEEWPKHQhV0gjqORIojld3pAsl4SRCrl4ET9usqS974zpiIWeIDL9jx/x4wpZRhR6gsj0P37EjytkGVHoCSLT//gx3Y8rGThCEMSHniDSBzoZTPXjSiWlEBSx0BNEpv+CFXHBCUERCz1BpF+KYEVccEJQRKEnjKnTf1PIU1pnUi64PF3jrCMulxSQ10BZ3tI6k3DB5e0aZx1R6IaT5xsubz7lJDJw8naNs460zzWcPLccdVrzdH3p8fgFSilOLhXlNQZm15YVF4xZSPvcFJPnQJlUdQbHbYbneC1zNiPMAhIUNZw856qv3NZTk5cNSFpnPXbWN3Auc4oKBJ4p1bxn5uw7GO69F4D9Na7HerxY6mYjFrrh5DlX3fSqzqSxs76Hb74PL9xy/+y2emVehWdKs0VL1musonp8UpZ6XhMDvCI+9BQgaWWCHar4ihfqYzFu50widlNfQQuUjZq8PtydfOii0AUhpTgGND3QvHTxrJFgpzxrSCAgnefEADucFLr40IXMkZcZTePCFkyNjQc+j13PmOHee23dNUnEbvKcGOAV8aELmSIvefujA4OYGZ8I7XzW3PP2njXo7L/LmNiNZDvpIwpdyATVoNnQhu25KJQZ2TqA0uR0qOe0WrwmBaTznBjgFXG5JERe3AJx4Or3hbnTc7/jIIrvU2/xmtJnSJrY6SMKPQGk73W42JWv12Pi9DzIOAjLf17FdIvXlIeL6YjLJQGkf0a4uFmrpior3XFgl4Otzhj3juT3ZwcthU5E3UT0GhEdIqK7bfZ3ENF+InqRiA4S0brwRc0OErUPFyfr22RlpTMOVEHeyRCt84ljpzCydSBzgeM84qrQiagI4AEAVwNYAeAmIlpRd9j/APBtZr4MwI0AdoYtaJaQqH24qIJmXY9tQfeRh41U5oDeOFBZ8VQMcXKd4WygvKEzKi4HcIiZX2fmSQDfAnBt3TEM4PzK3+8C8MvwRMweErUPF5MyMrygMw5UVjzPlFBocg+BVR9sXY9tmfNZ9UTt9pPy/ejRCYpeCGDU8vo4gCvqjrkHwPeI6IsAzgPwabsTEVEfgD4A6Ojo8CprZpCoffikMWimMw5UzdkK581DcX4TSlbXC6FsWlmoKulqRWX1s1QVplG5/SQRIB5cS/+J6AYA3cz8hcrrjQCuYOY7Lcd8qXKuvyWijwPYDeDDzGzfGQhS+m8yklJpDqMDgxi++T7w1IzjccUF8zyV68ddTi/l++ERtB/6GwDaLa8vqmyzsgnAtwGAmZ8FMB+AOIRTSF4qLU1ntlBq4w5lx0QrTn51O1993G4/SQSIBx2F/jyAi4loORE1oRz03FN3zDEAVwIAEX0IZYUerA2ckAiSUpk89Q9VlPQacPFMCdRYrNmmUtJxxx0kESAeXH3ozDxNRHcCeAZAEcAjzDxCRF8FMMTMewD8NwBfJ6L/irIX7/OcVBtHIRBiSSWPTqGUCiJCQ1srpk6fcXWXxRl3kMVK4kGrUpSZ9wLYW7ftK5a/XwHwiXBFE5IgiyskRR0TCPv8QR6epclpzGuZj8+e6vd9jiiQRIB4kNJ/oQbTLSmvyjPq7Ioozq96qFKxAC4xmjsWYcm6Thx+cJ/t+02dTfmdEUiQXh8p/RdqMDmn20/ANuqYQBTnVwUsO/vvwvrS4+g+8jBW7+wr/0Y2pHk2VY8E6b0hCl2YQ3vPGnQfeXhWeZigzAF/yjPqmEAU59d9qOahQE2C9N4Ql4uQGvwoz6hjAlGdX8c9kQe/tATpvSEWupAa/KS+RW3FhnV+v2Xxps6mwkLSHb0hCl1IHF1l5kd5Rh0TCOP84idWkwe3Upi4lv5HhZT+x4+J2QJ2qw0VF8xTKkUTv0NQpCy+lvrfeMm6TpzYO5yp3zwITqX/otBzglfFGReizIAnC9fbN8uy6cGSdUwdpyYRtJeLkAFMzRbwE/TKWhtW8ROfw9RxmhZEoecEk7IFrAqZCvaLqVGBbBV2Fv3N4ic+h0njNI2IQs8JpliB9QpZ1UmQZ0q2CjuLFpzJxVxxY8o4TSui0HOCKVagqvEUFQsAkW0LWKvCzqoF5zf9MGvuJ1PGaVoRhZ4TTLEClUuqlRjrS4+DFa1iq+9rXNhiu1+13STCVr5ZdD+ZMk7TilSK5ggTlmlzq6x022/vcVdvN4Uomng5uZ+S/p2DYMI4TStioQux4jaldts/efqM7XlV200hCt9/Vt1Pgn9EoQux4jaldtvvFDQz2Z8chfKVAKJQj7hchNhxm1I77Vf1a1+yrtPoVeWjaOJleu96IX7EQhdShcqCP7F32Oh0RjtXUqGpAdNn3vY9o5AAolCPlP4LmSAN5fMv3bELR3Z9r5xjX6ByJHfmnMxZLnHPYg+epJDSfyH1uPnHk/Ynu8k3OjCIY/37zxVSlbhGmQNmzSjCJIvplaYiCl0wHh2FkGRBio58qoKqerKYoZLF6l5TEYUuGI9KIby8efesVTyydQAdvWtD9SfrZs3oKCxdRZ3FDBVJr4wPyXIRjEd140+NjWNqbLx8zNGTONa/P5ASn/XzHj1Z9nFbqladsmZ0FJYqy8VKVjNUol4GUDiHWOiC8eje+PVWsZe89Jfu2IWhjdvPKR6bFgTW8+t0jLTKbecSosYimtpaM5+hIv1Z4kMsdMF4Vm7rwQu33I/S5LTrsVWr2Eup/ejAIA4/tA/QSPiaOHZqzrl5Zu4b6xVWHhZ0VpHn7x43krYoGM/owCCGb74PPDXjemx1pSMvKyGpjlWdH4Dt8VQsgEssCkuIFKe0RbHQBeMZ2Tpgr8wJNVa11Sr2EojTDc5Vzz+0cYft/mrHSEFICvGhC5ETtMeKUuEyfPV80dk2BwI6eteivWeNr5x3k/vMCNlBLHQhUsJoG6vMknBYSNqpz0l91eJ577/A3eXCwLH+/Wj7xAdtzw0ClqzrtH1rFK1zBcEOsdCFSAmjqMRPloSqz8nYgVfPZbNUioBO/eCglhzWXuMdvWtrm7BXFH79+qf7lt2KoQ3bpbBGiAUJigqR4qfHil3fDyB4lsTowCCGNm7XymZxonnpYqVFX5011FvlthjUZ0ZIDxIUFRLDa1GJnXtiaON2LL+tW+le0WVk60BgZQ6yz3CpUvX365T6S2GNEDbichEixau7xFYRMnD4wX14elFvoGBiKKXmLg+EqpJ2+6y0FNZIMDddiIUuRIrXohInRTg5Nu4pmFjvumla2ILJSquAKLAqaadS/+ali1ORpy7B3PSh5UMnom4AOwAUAXyDmf+nzTGfA3APyjbMy8z8p07nFB+6YIdOkY9bAc/owCAObt49R3lTYxFEVFtxSkDLh9px5tXjtuX+uliVtOrz09bv3EtxlhAfgfqhE1ERwAMArgawAsBNRLSi7piLAXwZwCeYeSWALYGlFnLJym09tdkjNvBMSdmmdnRgEC/ccr+tJc5TMyi2Ntdkviy/rRsTR96co8ypqagtc1Nba40yf7HvwTmf39jWmiplDkiXxDSi43K5HMAhZn4dAIjoWwCuBfCK5Zg/B/AAM78FAMz8ZtiCCvmgvWcNxg68qt1bxZpKCAAHN+927PkydfoMPnuqf/b1vmW32gYvedK9zUAVqytIFQxtaJmfKmUOmN0lUVZAskcnKHohgFHL6+OVbVYuAXAJER0goucqLpo5EFEfEQ0R0dDJk3q9M4T8sXpnH7oe3YLGtlat460Wo5uPvF4ZhWVtVh8sWbJqTe2SKCsgqQkry6UBwMUAPgXgJgBfJ6J31x/EzLuYuYuZuxYvXhzSRwtZpL1nDT57qh9dj22ZdZFQ0X64aluMhFllVM3esM2R90nVWrTDBKvWK6YuQi0rIKnRcbm8AaDd8vqiyjYrxwH8mJmnABwmop+hrOCfD0VKIVZMms6296yZ/Wy7Yp16i7GxrXV20Ys5cFkZjB14FUd3f1/tmqlr+qVL9Vq5yZgmrNffFLI0CwobHQv9eQAXE9FyImoCcCOAPXXHPIWydQ4iWoSyC+b1EOXMBSbk/Jo8na2xGFHOdqlaZlX5Lt2xyfEcE0dP4vCD+5x7qzNcA7P1VJW2qVZtlsjSLChsXC10Zp4mojsBPINy2uIjzDxCRF8FMMTMeyr7riKiVwDMAPgLZh6LUvCsYUrOr9N01gSlVJXB6Vq9cOuDKP3efUFmR6pKnStpkjMl5aHNSxdjybpOjGwdwNDGHYnParJO1mZBYaJVWMTMewHsrdv2FcvfDOBLlX+CD0xRpGmYzqqu1XDvvQCA0tnJcD6Iz+WXq3rAVFMWTXgY5wVZAUmNNOcyBD9NrKIgDcUkymsFAFS21maCWugWGttaMX1mAvzOXDcNNRbRcP4CW7+9SddMyA6BCouEeAjbL6jyx7v56b2kqvn1+ft5n86izAAARqjKHACmxsZtlTlQLlZSBWFNmtUI+UB6uRhCmH5BlT9+7MCrONa/39E1oDud9evz13nfbJbN0ZO2/mu7RZlNRIJ0QtyIy8UgwkoXVLlNVME9P64Bv64Zt/dp9RE3jKa2VsxMTM55GEt2ixAF0g89JdTn/L50xy4M995bVsIFQqG5CaWzk747FqoyNfy4BvwGT93ep9NHvAafOeOhQeXq1Ka2VhSamzB1+kwqg3Qm1R4I/hEfuqG8dMcuHH5w3zklXOJyKp5Gbrhqqh+40lLjPdXtKj+52/u8PlyW32bbZSIeLA+TybFxlCYm0fXoZnQfeThVytDk2gPBG6LQDeXIru857reWOtcrzyXrOm0Dm8v6rgqtN4dT8NRJQbgFXT09XAqEww89g+J589yPDUDB7vw2M4O0lp+nqZTehOI7kxGFbihOhSxVJo6etFWex/r3o6N37ZxqxdU7+wJVMVpvppGtA7af0d6zxlVBFJqbZrc31bWVtVP4SkoMMIee1WKFigX7IiWFmyeNmS1pqD0AZCahgwRFDeWphhtclToVC5h/UVsseeOqPip2DwSnPPHignm1yr5i6dYvEFHNckkbbotvmEgaag+A9MgZNZKHnkKW9V3legzPlJRKLyzrqmqVD23Yrj0tb1zYYn8ygu16oUBlMegN2/EkXYehDdsxfeZtdD22BSCPTVXixEY0p8U3TMXUNrn1pGUmkSSi0A1l9c4+LLpyle/3Ny1scfU36uyfneIqsLuZlCrYw2Rwamwcwzffp344JA0Bi/7TKmWgGTDXD11PWhqKSVMudyRt0VBGBwbx1rOv+X7/5Ftn8MIt9892FZw4ehJDG7dj7MCrWL2zT6vARyeFkAqE0YHBmqKgsBZi5qkZEGzcNCbAwFvPvubqFkuL9Whim9x6pCmXO6LQDcVzPnY9JZ7bIpaBww/tQ9snPqjVDExHGfFMafZBAKDm73rsCnDcmBwbjzyLxS8638POepScb39IUy53RKEbSmSWHZ+7IWw/1+JeUa0pWY/VtaBScsUF87Cq0qvca8AzyiyWKLGzHk1pk5xW0jCTSBLxoRtKlH5Bp6XSQJj1pXtJIZw4dsrxIVT1ybb3rEH3kYexnp9Ak+aaoWmCigVHP3Sacr6F9CEK3VA85WN7pDpVtY1eViz4qltg5uw7s4G/5qWL1Qs3Myu7IDYvXWzb3CssX7sxUNkF5eQKkEwNIUokD91gRgcG8fLm3eo1Mn1gzR1/kq5zPM4u5xyAp+ZZhaYGFFubMTU2PtscrKmtFVO/Owuemgn2ZWLCacWi2X11laOqHH1VLjVQm4svCCokDz2ltPesQUPLfN/vLy6Yh+W3dyvT0aprc9ZTXavTijVgal3X0/Hzz5sHZp59IFWV4uTYeDTK3KlPut9TNjWA5jfa7qPGIjr77ypfC802AE4zrzTlrgtmIkFRw/E8FbepvFShSgNTWd9VWaq+cMeVgwCU3p7SamFQT2NbK6Z/83vv7y2FO9ucnUkogrJV+by4UWoyNWwsdZPWbxXSh1johuOlsKZ56WJ0PboF6/kJrY5/qoISlfVdH0h1C9z6UebFBfNw6Y5Nvt7rhca2VhQUlncVBpxnEiWurDVq/yBRXZ9qYFhVBSv+dMEvotANZnRgEDPjE1rHVvtZeLXsZrNOSo/Pvl+3FDzswK21UZdTBWYYlCYmsfTmKx2P0YpdKCYFOgUvUvkohI0o9ATQbQE6snVgbnGQgjCtOqdScLuOi+pafw0q5+96bAv++FQ/gHLgMGoLfebsOzj+7QOR+N11S+fT0kNFSA/iQ48ZL4UlXpR02FadXQGHnezH+vf7XjGovkte3MvPTY2No2VFO868MhreSYm0O/9J5aMQNqLQY0an5L6KbqUmgFisOpXsTml9KuwsUd12B8WW+Si9M1Xj36bGIhrPX4DJ02ccA7X1hKrM4f3BKpWPQpiIyyVmvGRE6PqoG9taY1EKTmuVevGlN9YtauF2fgCzbpn1/ASuGf8HdH7zizUuoc5vfhF/fKofXY9ujtz/roTiebAKggqx0GNGZXXbWXZuKW5A2TK9tNIjJWqUsldSJA9u3j2n+nO2sEhj8WSn89e7MZxcQlH7322h8vqmYm0LSSIWesx4DYRVs05UgcfG8xfEpkScZG/vWVO2kB/bUmM5f/SRO/HZU/01WTR+zq9D4A6VPqmmi67e2We7X9bBDB+5pvaIhR4zfgJhI1sHlIHHydNnohDTFh3Zg/iEgwYJk8jftj7QrNQso2dpCyDdFYMjHSvVSC+XFOBUkZm39RSdUPZJqeuzEjZ+snXkd/NP3tcWderlIha6YdgtfqDMdpEgXA2qVgbv+fgHcOoHByP73ImjJ8tKpvKbTZ9529X1I9Wg/pGOlWrEh24QNWt4WhYaXrKuc24WSYhBOB1/ZBp8lqqCqN8f+lW0H0yo+c10KkylGtQ/UmGrRix0g1DleZ/YO4zLdt0eSQGKjj8yTT5LOx/+0MYd0X6oR3eOVIMGQ9YWVSMKPSHsXCtOU8moClB0Cp28FEOZiMplpSqIomIBXCov2BFaCqSHLpiCM1Jhq0YUegKoLN6mhS22q/hEvRyd2/a0+yxVFl1H71oc69+vXMhjaMN213PrVMlSsYDO/rtE4YSIVNjaIz70BFBZvNNvT9oev2RdZ2Sy6Pgj0+6zVPnWV+/sm7O9o3ctXt68W0uZFxfMw7K+q1yrZLnEonyEWNBS6ETUTUSvEdEhIrrb4bjriYiJyDalRiijsmxLioUUTuwdjkwWnWKeLHQFtGsTXL995bYeHOvf7xzUrBR41T8UnNoNpOXBJ6QfV5cLERUBPADgjwAcB/A8Ee1h5lfqjmsFsBnAj6MQNEt4aboFROva0C0WcjsmC+hUmnY9umXO966+fuGW++e0O6bGYqoefEK60fGhXw7gEDO/DgBE9C0A1wJ4pe64rwH4KwB/EaqEGWN0YBDTZ96es724YB6KzU2x+9ABPX9kHnyWbg9OpyZo1e3WfjaNba24dMemzF+3uLFLKJBrXEZHoV8IwNpj9DiAK6wHENFHAbQz89NEpFToRNQHoA8AOjo6vEubclQVhE1trVhVabClk46VpQFt0ndxmznNjE/M5t/byVz/0Kt+t6GNOxL/blkhTSm0SeBa+k9ENwDoZuYvVF5vBHAFM99ZeV0A8EMAn2fmI0T0IwD/nZkd6/rzWPqvU7LspuBeumMXDj+0ryb3uZqZkbYBbfeAS/K76JTsN7W1YmZics4xxZb5uOyh2wDAtocLkN7fKSr8PMzzXvYPOJf+6yj0jwO4h5k/U3n9ZQBg5r+svH4XgF8AqHaJWgLgNIBrnJR6HhW6W08WtwE9OjBYWZTY/v1pG9Am3pw1TbW8QgRqKDguLJ3G3yksrAq8aWELpn53tuZa6TzwlPcQEdaXHo9CbONwUug6WS7PA7iYiJYTUROAGwHsqe5k5t8y8yJmXsbMywA8BxdlnlecfOHVqaNTSb1T18W05IRbiSu/3UvbgmrWS/PSxd4/iNlRmQPp/J3CoL6txeTY+JxrVS1WcyLtKbRR46rQmXkawJ0AngHw/wB8m5lHiOirRHRN1AJmiZXbelBoUoct3Aa044o+KRzQcdycqv44br1oVKmajW2tgeRJ4+8UBrq96t0eeFlIoY0SrTx0Zt7LzJcw8/uYeVtl21eYeY/NsZ8S69ye9p41KLY2Ox7jS2mnrOti1WKe9TNbCPvmdGpb4ISqGOnSHZuUi424UVwwD+e9/wI81XADnqTr8FTDDXjpjl3+TpYydGcmbg881e8icYkyUvofM1MuC1JYB3R90GjJus45peppW/psTuCREVqfE6/9cdxQpWoe+eYPPbfjpWIBzcveW/M+ninh8IP7AEC52lGShJmBpFN7ofswz0MKrV+k9D9mnCwQ64C2cxUc69+Pjt61NdaJ09JnJmI79eZzwcIgytzOtVJc0GR7fNPCFt+f89azr3l+H8+UcOaVUdt9R3Z9z5csUeLXVaXCzlVSaGoou7DE0g4NsdBjZHRgEDM2RUXAuVx0tw6HJ/YOpzpLIqpAqOp6oWDvH/G7gFEU65Ymsqi1C2F32MxLtXHSiEKPCVWOs6qaMO0dDlWoOkr6tZirKK9LyV51u7m+PH9OEBQPnSSJYvyJqyR6xOUSEyrLrqFlvu0gz2p6lsoyDrrkp+q6qJpm+b2OyvcFVMpVV4YpK0NldfxlHVHoMeHV4gmSnmWKUrBDZRn7tZirqK6XXXvbIJk0qs9ZfutnXNvoKikxXt68O3S/dRAkPTCdiEKPCa8Wj9/0rKSVgtvDJCrLz0vP8yDBN63PATynNk6NjePlzbt9pVhGgaQHphPX0v+oyFvpf1x9S5Isp9f5ji/dsWs2Vc/K8tu7U5Wt44TqN/BNCsvaTWq6ljWClv4LIRCXxZNkMFWniEe1WEeUi3jETdjXutElYGyaiy3pWWKekSyXGIkjyq8q4IgjmJWH9Ul18LqAiRvVtr12Y8fEdrJpX1Q8zYiFnjGSDGZlaX3SIFav3W8QhNLkNF7evNtWvqEN222V53DvvYlZxHl4aJuKKPSMkWQwK+n1ScNyPQR1Gdj9Bstv7/bXwbHC1Ng4vruoF6MDg7XyKeCZUmJujrQ8tLOIBEVTjImBJx2ZopA7zKBzlIHl0YFBDG3Y7vv9TksV2uFH5qC/j2kLl2SNQAtcRIUo9GCMDgxi+Ob7anpKU2MRnd/8Yi5vmjCVcBiLKDgpRVWmTyR4zJAJSxmbaGxkBclySYgosw9e3rx7zgIBPDUzx9eaF5R+26MnPf8GQV0Gbi6b1Tv70PXYFq1zBcWrm8Nvu+F6qguFdD26GQAwtHHH7PU3LSsnS4hCj4ioU7emFFNu1fas49QrPqzFLXT9/DpKsb1nDZp8LpbR2NaqFXT1E5sIM6Bpdw8M33wfXrjlfklpjAhR6BERlqUjqLFaetNn3p67GlTdIs1AsMUtdF0Gukpx1Y5NoI6oBeIAABJGSURBVAZvt2BxwTxc9LlPoNBs3xa4it9geFgBzdGBQQz33jvnHuCpGZQmp2u2yX0RHpKHHhFOLoAwaGprte9aGHCJtHpM9YXW+3qnxsZBjcXydTl9xjEXPMjiFjro1gJUzz/UuwOYcY9lNS9dbL/Iic1xfoO3K7f12PrQvVj61d/GS1tgSWkMB7HQI8LJBRDG9HLVjk1zLNJCUwNW7dgU+NxVTKz4c8q/5qkZFFvmY33pcceFnqNOn/PismnvWYP104+jZUW7/ckKhOW3d2M9P4HuIw/jxN5hR2UeNAU0jLRXPz3jJaUxHMRCj4iV23owtHH73L6wXA5oBrVy41gwwLSKP1VPeStWSy8Ma9MPOr+N3czH7T2Ay5qzAZfws8of5BxOMlJjEURU43aRLo7hIWmLEfIkXafc1/XYFiNcF06Ekb4XJjpNr+rdDWG5jMJ0PQVJDQwzPXN0YBAvb949G0ivXzXLLyoZqVhAZ/9dAGTloiA4pS2KhR4hzUsXKxVQGvpaJNkXxg43P6udpRdG/5yw+6UEmfmENeuwq2OYHBvHC7fcDyBYH5gl6zpt8+yX9V01e17Tx35aER96hDjdZLpBoCRzdk1b5MDpQRJli4OwM5aCpAaG1dphZOvAnDoGoNw3JmjGSR46apqKWOgR0t6zBgc377bNRtGxckcHBvHCLffP+hsnjp4MxYLSxbSFfVXWadQl5WE3mwo68wlj1uEke9CME2nOlRxioUfMqh2bfFu5BzfvnpOzW5qcxsEYq0GrFX/VzJEkp8pJNR7zm5utml2ZMPNxnO0EcKmNDgyCFOurhuWqk0pTNWKhR0wQK1fVgEm3MVOUJJWfnsTK8X781jp+9yRnPiu39czxoQPl1FfdB0v9GKjmyNvln4fZUdO0/u8mIVkuBjJ7ozhkdKznJ2KUqBY/WRqmFijp4lX+JJcC1CVIlottCqlNZS5wLrtF9/d2utZpuK5RI1kuKUIn17ox5GpQr3jN0siCVeV1ZpAGP3KQ2Y5t8ZDCNuQSe1LmTmMlDdc1ScSHbhhuVXbUWMSlIVaD+sHrTZXHvjZZX+TBiwL18p3dxkrWr2tQRKEbhlsl4LIv/BFGtg4kGhDyelPl0aoyIfAZJU6tLax4/c5ubZCXrOvM9HUNiij0mNCNzCuVZaWs+1j//sR7q3hVVnm0qpJcCjAOVGNg+W3dgb6z05iYOHoSx/r3o6N3bWava1DEhx4DXnzIqoyKJes6Mdx775wMgiR6q1Q/y5pj79TONameKmHiJ6ibREZOXESVqWM3VqzMnH0HJ/YO5yYA6hVR6DHgJYhod6MsWdeJI9/4V2U70qRcFzMTk7N/T42NKx9SJqTpBSELQd0oiOKBVTNWArQ/ziuSthgDQZtcfXdRr+NKREmkbOUpfSyt3zXtqaJpve5RE3hNUSLqJqLXiOgQEd1ts/9LRPQKER0koh8Q0dKgQmeJoD5kJ2WelOsijEBnWir+0hjUNbGXvVfs/PQAMH3m7VR9jzhxVehEVATwAICrAawAcBMRrag77EUAXcy8CsB3APx12IKmmSgzHpIKCEW9kLJJpDGoa0KqaNAHdjWwXL8KV9W9Z+JYSRodC/1yAIeY+XVmngTwLQDXWg9g5v3MfLby8jkAF4UrZroJmvGgWlauqa010UZZUS+kbApBvmtSs5CkZxVhPbDbe9ag2DJ/znZTx0rS6ARFLwQwanl9HMAVDsdvAvAvdjuIqA9AHwB0dHRoipgNnAJIbr7OVTs21XRdBMJfbs4rQQOdSSscL/j9rkkGU5PuZR/maldpGitJE2qWCxFtANAF4JN2+5l5F4BdQDkoGuZnm4yTwta56esj/1Qs1PStTspKj2MhZVPw812TXMIv6VTRMJVw2sZKkui4XN4AYF3B9qLKthqI6NMAtgK4hpm9rRCbYdymnrquh/aeNbNT/2r6Yv250hJkBLJfSQkka1kmXdgUZtwhD2MlLHQs9OcBXExEy1FW5DcC+FPrAUR0GYCHAXQz85uhS5li3Kw0Lze9m/JPU6502nPTdUjaskyysCnMGUIexkpYuCp0Zp4mojsBPAOgCOARZh4hoq8CGGLmPQD+BkALgH8mIgA4xszXRCh3anBT2F5ueqdzJTm990Lac6O9kLTbI0nCVsJZrroNEy0fOjPvBbC3bttXLH9/OmS5MoObwvZy0zudKw2Bo7xVXObdshQlHD/SnCti3Px/XnydqnMtWdcZ+bJfYZCmVMWwMGkJv6RJU4wnrUgvl4jRsdJ0LRlVn5eol/0KizTMIoRoyNvsLCmkl0tChOVLVvW78LrsVxxIb478Ir99eATu5SKES5hl7yrrlkvlB7XbFDfOabCkn+UXp4UrhPAQhZ4AYfqSVT7yxoUtrg+NuPupJJ0bLSSH0wpHXseb+OLViEJPgDB9ySqrlwDXh0YSQUqnIKHcqNll5baeOcvTAQAYnsZbmpq6JYEo9AQIs4pOZfVOnj5je7z1oWFSkFJu1BygCNd5GW95zJTygij0BPDiS9axWu2sXp2HhkltYeVGzS7Vh7UKL+PNJCPEREShJ4CuLzmI1arz0DApSCk3anaxe1hX8TreTDJCTEQUegLopiwGsVp1HhomBSnlRs0uTg9lr+PNJCPERKSwKGa8FFgEtVp1Cpasx1QfNEMbd3jKjQ8jpz7PfU+yjrJlxdLFWj3lq2OrcWHLnGA/FQvo6F0rmVIVxEKPGS9Wd5xWq1/3Tpgr05gyWxDCxa9VXT+2psbGMVm3vi7PlHCsf78EzytIpWjMPFm4HrC75kRYX3q8ZlO9NQ+Ub4SO3rU4sXc41IZPfiv5pAJQ0MHPLE41tuzI03hzqhQVl0vMeGmX69S7JeyeGH7dOxLMFHTw03nRyxiS8VZGXC4x43X6WZ+SeGLvcCTpfX7dO2kIZkrBkhqTr42nMcSMpxpuwEt37IpOoBQgCj1m3HzFbjdYVBaxXz+n6VkHUrCkxvRrs2Rdp6fjeaaEww/uy7VSFx+6Qah85laFH6XP2m+2ismrEImPX43p18aLD90KFQv4k+nvRCCRGYgP3QB0lJ7OMnJRpvf5XWHG5JVpxMevxvRr41cOu7UB8oK4XGJAd2qrskas24Ok940ODOLpRb14kq7Dk3Qdvruo14jpdZR+3DT4+JPC9GvjVw4q5let5febx4hu7rlqINZv97Os2ejAIF645f6aPN6psXEM33xfoko9aj+u6T7+JDH92tjJV+3Y2Lx0MRZducr2fcv6ropYMnMRhR4DulNb1VQxjCnkyNYBlCan5557aibRBlhRN+XyMqMxOeMjCkwv5rKTr+vRLVjPT6D7yMP4j9+/B8tv7541eKhYwPLbu7F6Z1/CkieHBEVjQDf4FGWQSlnQBNgWNcWFl0KrKNEJSAvesSvdnzx9xrjgeZqQJegSRndqG+UU2MkfmaTP1BQ/rrTvDQfrLOfpRb0Yvvm+uaX7BqZIZgVR6DGgO7WNcgq8clsPCk1zk5qosZioz1TlJ504ejJWt4fpGR9poD4eMjk2Dp6aUR4vD8zwkbTFmNBN7YsqBbB6zoObd88GRhvbWnHpjk2JTntr2hscPVkOelU8MGG1NdDBS0sGwR6nvucq5IEZLuJDF2pIskgoyUIX8aEHxzFOo8CUIqY0IT70lJFUtkXSpeBJuj1Mz/hIA15nM0HiQ3nLSNJFFLphJKlUX9682zYwOLRhO55e1IvvLuqN9AZKOkDa3rMGK7f1lN0vx05hZOuAKAoP2MVDCk0NaGxrBYjQ2NaKpsrfQR6YSRseJiM+dMPQKf+PgtGBQUzVLR5gxVqQVL2Bxg68Gmpf9qRXLfKympQwF7t2z1G47JK6R9KAKHTDSMrt4DXbYObsOzj80L5QA5hxKQQVoiiCE0dfH8lIUiMK3TCSyrbwdTPUxb+q7pmRrQO+FXGSjb50FYXJ3SXzgGQkqREfumF4KS4KMzAU5s2QVp+mjg9f/LfJY3oPmiQRhW4YutkWYSsW2wKfAKSxaERHUUhFafJIRpIayUNPKVHkbNe7Epas65wNenrNLwaQaI8Yv7i5U0zpPSPkF1ngIoNEERhy8l8/1XCD566PQdw4Sfmp3Xz44r8VTEbL5UJE3UT0GhEdIqK7bfbPI6J/quz/MREtC1vQPOLkI487Z9trj+mgRSOm+qnFf6uPFP/Ej6tCJ6IigAcAXA1gBYCbiGhF3WGbALzFzO8H8L8B/FXYguYNN6UWt2JZvbMPy2/v1jo2qE/TZD+1+G/1MPmhnGVcfehE9HEA9zDzZyqvvwwAzPyXlmOeqRzzLBE1ADgBYDE7nFx86M7o+MiTcEu4ul5C8CWLnzr9mL4AdZpx8qHrKPQbAHQz8xcqrzcCuIKZ77Qc89PKMccrr39ROeZU3bn6AFSXE/kAgNf8fSUsApCmKgLP8r4P7+5U7fsFfjMcWCI1jrI6yQUA0yhNHsXvfhJEgKU4/yMNKDRpnjtNYyFNsgIB5E1g/Obm2gJYysyL7XbEGhRl5l0AdgU9DxENqZ5QJpImedMkK5AuedMkK5AuedMkKxCdvDpB0TcAtFteX1TZZntMxeXyLgBjYQgoCIIg6KGj0J8HcDERLSeiJgA3AthTd8weAL2Vv28A8EMn/7kgCIIQPq4uF2aeJqI7ATwDoAjgEWYeIaKvAhhi5j0AdgN4lIgOATiNstKPksBum5hJk7xpkhVIl7xpkhVIl7xpkhWISN7EKkUFQRCEcJFeLoIgCBlBFLogCEJGSIVCJ6KFRPSvRPTzyv/vsTlmNRE9S0QjRHSQiP5zAnKmpkWChqxfIqJXKtfyB0S0NAk5K7I4ymo57noiYiJKNH1NR14i+lzl+o4Q0T/ELaNFDrdx0EFE+4noxcpYWJeEnBVZHiGiNyt1L3b7iYjurXyXg0T00bhlrJPHTd6eipw/IaJ/I6JLA38oMxv/D8BfA7i78vfdAP7K5phLAFxc+fvfAfgVgHfHKGMRwC8A/CGAJgAvA1hRd8wdAB6q/H0jgH9K6HrqyLoWwILK37ebLGvluFYAgwCeA9CVhKweru3FAF4E8J7K6/caLOsuALdX/l4B4EiC13YNgI8C+Kli/zoA/wKAAHwMwI+TklVT3n9vGQNXhyFvKix0ANcC6K/83Q/gT+oPYOafMfPPK3//EsCbAGyrqSLicgCHmPl1Zp4E8C2U5bZi/R7fAXAlEVGMMlZxlZWZ9zPz2crL51CuP0gCnesKAF9DuYfQ23EKZ4OOvH8O4AFmfgsAmPnNmGWsoiMrAzi/8ve7APwyRvlqBWEeRDmLTsW1AP6eyzwH4N1EdEE80s3FTV5m/rfqGEBI91haFPofMPOvKn+fAPAHTgcT0eUoWxy/iFowCxcCGLW8Pl7ZZnsMM08D+C2AtlikU8hRwU5WK5tQtnySwFXWytS6nZmfjlMwBTrX9hIAlxDRASJ6joj0up6Fj46s9wDYQETHAewF8MV4RPOF13FtEqHcY8b0Qyei7wNYYrNrq/UFMzMRKXMtK0/kRwH0MrO3Bt7CHIhoA4AuAJ9MWhY7iKgA4H8B+HzConihAWW3y6dQtsoGiegjzPybRKWy5yYAf8fMf1tp1PcoEX1Y7q3wIKK1KCv0/xD0XMYodGb+tGofEf2aiC5g5l9VFLbtFJWIzgfwNICtlSlXnHhpkXA84RYJOrKCiD6N8gP1k8z8Tv3+mHCTtRXAhwH8qOK9WgJgDxFdw8xJtPPUubbHUfaXTgE4TEQ/Q1nBPx+PiLPoyLoJQDcAcLmb6nyUG0sl5SZyQmtcmwQRrQLwDQBXM3NgXZAWl4u1tUAvgP9Tf0ClLcGTKPvQvhOjbFXS1CLBVVYiugzAwwCuSdDHC7jIysy/ZeZFzLyMmZeh7ItMSpkDeuPgKZStcxDRIpRdMK/HKWQFHVmPAbgSAIjoQwDmA5jbF9cM9gD4s0q2y8cA/NbiqjUOIuoA8ASAjcz8s1BOmmQU2EO0uA3ADwD8HMD3ASysbO8C8I3K3xsATAF4yfJvdcxyrgPwM5R991sr276KsoIByjfDPwM4BOD/AvjDBK+pm6zfB/Bry7XcY6qsdcf+CAlmuWheW0LZTfQKgJ8AuNFgWVcAOIByBsxLAK5KUNZ/RDl7bQrlWc4mALcBuM1yXR+ofJefGDAO3OT9BoC3LPfYUNDPlNJ/QRCEjJAWl4sgCILggih0QRCEjCAKXRAEISOIQhcEQcgIotAFQRAygih0QRCEjCAKXRAEISP8f6aEsBYlYVUkAAAAAElFTkSuQmCC\n",
                "text/plain": "<Figure size 432x288 with 1 Axes>"
              }
            ],
            "_view_module": "@jupyter-widgets/output",
            "_model_module_version": "1.0.0",
            "_view_count": null,
            "_view_module_version": "1.0.0",
            "layout": "IPY_MODEL_b974db90bbbf44c482c1dd12f61013a3",
            "_model_module": "@jupyter-widgets/output"
          }
        },
        "2d5a3a03f7cb421cbb0cae446cbdfdb1": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "SliderStyleModel",
          "state": {
            "_view_name": "StyleView",
            "handle_color": null,
            "_model_name": "SliderStyleModel",
            "description_width": "",
            "_view_module": "@jupyter-widgets/base",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.2.0",
            "_model_module": "@jupyter-widgets/controls"
          }
        },
        "f17ea6f4ce5245009f234130c9273294": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "b974db90bbbf44c482c1dd12f61013a3": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "43db062d32ac4b39b5f554d12c335620": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "VBoxModel",
          "state": {
            "_view_name": "VBoxView",
            "_dom_classes": [
              "widget-interact"
            ],
            "_model_name": "VBoxModel",
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.5.0",
            "box_style": "",
            "layout": "IPY_MODEL_a4dc0e31bfea45bcb9f064488a7c40ed",
            "_model_module": "@jupyter-widgets/controls",
            "children": [
              "IPY_MODEL_7b0b70baaaaa4da68b7e71469ee49731",
              "IPY_MODEL_a6fb1ecdfa004558849b54cd293b6645"
            ]
          }
        },
        "a4dc0e31bfea45bcb9f064488a7c40ed": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "7b0b70baaaaa4da68b7e71469ee49731": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "IntSliderModel",
          "state": {
            "_view_name": "IntSliderView",
            "style": "IPY_MODEL_197250af0ee648579e97f82dd8006c14",
            "_dom_classes": [],
            "description": "t",
            "step": 1,
            "_model_name": "IntSliderModel",
            "orientation": "horizontal",
            "max": 19,
            "_view_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "value": 10,
            "_view_count": null,
            "disabled": false,
            "_view_module_version": "1.5.0",
            "min": 0,
            "continuous_update": true,
            "readout_format": "d",
            "description_tooltip": null,
            "readout": true,
            "_model_module": "@jupyter-widgets/controls",
            "layout": "IPY_MODEL_0071e8e86e8949f088f9d82f1e5fa3e6"
          }
        },
        "a6fb1ecdfa004558849b54cd293b6645": {
          "model_module": "@jupyter-widgets/output",
          "model_name": "OutputModel",
          "state": {
            "_view_name": "OutputView",
            "msg_id": "",
            "_dom_classes": [],
            "_model_name": "OutputModel",
            "outputs": [
              {
                "output_type": "display_data",
                "metadata": {
                  "tags": [],
                  "needs_background": "light"
                },
                "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29f5RU53km+Lz1S10FarW6seOsEC17MBOrLeImiElWezZk5RZyI2GOLTTONITYzmCEMgd7NqtjCTubtYVQtDtjMeeAZMYjB6POZITsZVqADxA7HO9h7ACmEuSWTwiR1RjFiaXG7RbqanV197d/VN3qe299P++vqur+nnN0RFfduve7v97v/d73eZ+XGGOwsLCwsGh9pBo9AAsLCwuLaGANuoWFhcU8gTXoFhYWFvME1qBbWFhYzBNYg25hYWExT2ANuoWFhcU8gdKgE9FzRPRzIvqR4Hsiov9ERJeJ6CIRrYp+mBYWFhYWKuh46H8G4F7J9x8B8P7qf9sAPBN+WBYWFhYWplAadMbY9wBck2zyUQDfYBX8AEAHEf1qVAO0sLCwsNBDJoJ93ALgp66/r1Y/+5l/QyLahooXj0WLFv3Gr/3ar0Vw+IWJ0rUJvPX6OGanZpDKpdF2UxsmRieAWVflb4pwU3cH8p0F7u9/OTLm3Z4AJFw4nMqlceMt7bUx/vzlf8Ls1Ax3u3ff8R7pvkS/ndtJ/fWQ/eam996MfGeBf60k11Z7PJr7mQ/4px++LvzuPb9xS+BtFyJ++MMfvskYexfvuygMujYYYwcAHACA1atXs/Pnzyd5+HmD4cEijm87gumpcuWDKQBvot4YzwLtsx14+PwjdfvYd9tTGJ8d837YCBWIKSDzz1n0f3kjegZ6sSf1GH+7MvDo+SekuxL+1gHnesh+42zLvVazAF0h3P/lB9Az0BtsPIIxNSuGB4s4vesUxq+MoX1ZB9bu7hOeux/7bnsK4yNj9V8QsOFzmzz7EW3b3t0a1yluENGI6LsoDPrrAG51/b20+pmFAkFfkNO7TmF6ouz9UGCMx69wXiLJ543A9EQZL219EQDQvqyD/zIv6/D8zbt2ot+64T9v2W+cbUXXis0wHN92BAC4901nPAAwPjJWMWIBDGWcGB4s4uTOo5gcLdV9Nz4yJj139z5O7zolvg4M+Pb2/+65l8v7V+DiwaLnGc8Usli7uy/cCS0AREFbHALwe1W2y28C+CVjrC7cYuGF42WPj4wBbO4FGR4sKn9rYozbOvPYd9tT2JN6DPtue6q2f7+BbDQc47i8fwUyhaznu0whi+X9K2rn8fSSx/HSJ79Zd+2W968AZeWPtP+8ZUbCCYPIrpUzGfHu29rdfXXnwgUh0HMQJ4YHi3jpk9/kGnMH0xNlnN51SrqP2jMuQfn6lOf8Lx4sYuXWXrR3dwBU8cz7D2xsikmu2aFDW/yvAL4P4F8S0VUi+jQRbSei7dVNjgN4FcBlAP8ZwI7YRjuPwPOyVS+IA6GBId+f2RTKb03VGYsTO46gfH2q/ucKYxg3pifKuHz8EvoPbPS8zCu39uLiwWLtPEqjE2DlWe5v29rbhPvneXk9A73ILOIbXVZd9qgMszMZOUZ4eLCIfbc9haEth5HNZ9HWlQcIaOvKI51Le3/MyVvoPgdx4vSuU3XXmAeZc8FdSWrAuZe1VdeVMZzedarhk1wrQBlyYYz9ruJ7BuDhyEa0QBA0FDI8WMTU9XfqPs8Usrjlt27FldM/AZthoDQhc0OmznBPT5Rx4dmzdUakrSuPX/nQr2LkO6+anUjEGB8Zw9CWw2jrzCPfWcD4lTEUD5wDm1EH+FWeoMjLExmdyWsV79T5zUtbXxSOwzHCV8+MeK5vaXQCmUIWGw5V4sT+UJEw3DMyhj2pxzwhmDAxbFPorgJlq5cwYT3H+XDujW6IZ6Ej0aToQoPsBdSNFfv3537IHeS7Cnj3h96Dke++WjMkbIZxvfDKl/yPR77bWGNeA4Nnqa9jzFVo7+4QGgKde+H8dmjLYXG+YmSMO1k6xr5noLf2nwNhshDwrKqunhnxxJXjNnA68X9VXFs3h8ADpUi4grUGXQxb+h8TVDFy3jJe9YLIlrBuYx4Ek6OlxrBcEkA6l5ZeV9170TPQK71GlCbj5LROnH16oozigXOBQ3RBsHZ3nzQEpxPX1s4hcMBm+ReymZL5zQjroccEWYzc7aWZLKFFD3NpdEJ/YA3gmsvQ1pWXJt6iQPbGnPS6mtyL9m6B10nylYRo5aUTzgHE+47LwDnjcrNc8l0F9O1dr+0h867r1PV3Qt3vZkvmNxuoUS3o5jsPfU/qMb7hJODRWTmfWoSnlzxuZrw5WPXQmjpKGKUpkrBGEGQKWe3EGaUJbJaBUubjbe+OJubMDXsRsGr7Glw+fknMtT60SXps4fPi7EJwj9q7O/Dwa63DzRaFDXWQKWQt2wUAEf2QMbaa950NucQEkScR1MMYHixicnyy7vN0Ll1hUeiMqbsD6/Zv9LBI8l2FhhlzkDgh6UemkMX9Bx/Ao7NP4P6DDxgv5aOiA/YM9KL/wEbPNc93FrD0rm5hmGLV9jVKIyR7LjKFLHq33WkcoguCEzuO4MnMLuyhx/BkZhdO7DgS6f6d69ferX4P2rrylrpoCGvQY4JOXNaht/k54jyIaGTZG3O4Z+99SgPnPnbPQC8efu0RPDr7BLKLcyanZYx8V6GeqgdUnjzJPNLWlUe+q8B9mf1GgdJUO5Ys7htlzHm6NF37d2l0opa4TJGXO5rOpbH0rm7l/oTx5mpy8PLxS7Fzs0/sOIILz5ytTfBshuHCM2djMeoPv/YINjy/Sfrc3v7gHbXn9OHXHrHGXAM2hh4TVHFZ/9LTzVq4emakRtWjNKF3253CWOnktRL3WMv7V+DHL/yoFqLJ5Pm3Os4kU6aQBQPDDEfPJH9zAdnFOXGJtyKM4GeLOFBVJkZxvqL8CI9eOTM1o8XM8N/Dts48ym9N1a6dU3Cja8SDUByLB84JP1+3f6PymKZQ5Q9eeeHlWI47n2ENeowQGR1AbBS+/ZkjKL8997njJWUXZT2fO3CW6v5jndhxBKVrc/H2ydESl+YWhlomgxOzHtpymPt96doE+vaur4unhg0jONdBqAcSIqmmmizCJi6dsQ8PFrlGzr3CkDkK/nJ9XYqjaPxxhuR6BnqFz8jkaAnDg0XrmRvAGvQGQfSS84w2AJRL03UJRJHxGx4sCvnQQ5sP49TOY2BgmLxWQltnpXqR50UHQb6r6nlXq/vaOvkslvZlHcpVTBis3d0nnSxMPVidZJ4wucwqfHOdc3OOI5wcJAU3/qImN3Q43KLxOyEt3lh1rqFqO5lTYXnnZrAGvUEw9oxnGfoPbNR6gU7vOiWNT7uZMpOjJVA2hXxXITSDJp1LY3J8sraf8ZExpHNpUDblif/74/lxvLCyyUIW7jJZUbmRKWRrEgW87XS9ZNVxKM0vuDm582ilslVy31UrhWVr38utFO7ddicAr2HmhYSGthzG0ObDHkaRzrVeu7sPQ5v5XrrlnZvBGvQGgedBykBp0jZ+pi8BK8+CgYXmhGdvzNX9fmZqxuO1R12yLvP+RNdLFO4a2nwYp3ed4o5Pdk3dBmzpXd3CsIxjeIPUGgByiqfOfVNVIb/+/Z/Wfd599/uwbv/GOsPMPV51MnEbbVU9BlC5T6d2HuM6FJZ3bgbLcmkQePQ3GRwvSQdBmiVMjpbERkHjKWnv7qhpn/hRujYRmK0gYwIFVayUGU3RPoQ01GoC122cHn7tkTqhNAdOXFgE0XEoTdp0PxGmrr8jPLZoZfCLy9ek34vgGG1dzaK+veuFrDATNthChzXoCcL/YAJAbvEN8h+lCKseWlPzklQPtoivHhTt3R1Y9Zk1UnqZ8+LFwb2XGWyR9yeSs9UdD4/eyKUVErC8fwV33JQSWPTquEUQ0V3vP1hppCH6Pt+lnsSdxDjv2qgMb5DQh7Nq4sH/uYeKWq2PyOQzGNp8GENbDjedvHCzwhr0hCAyTqo4evutN3mWvKoHW1f2VAsEPPzaI7h8/JLQO3PzoYPo08igkhhWNZ4QTXg8tUo//PvuGejFyq29Xs+bARcPFrmrBhkzxL1v3iTv59g75+wwPvzywv0HNnI9XB5EXHzhJMeAPfSYdIISwQmB6T4TNX76oU0ol8pzK8YmlBduVtgYekIQGSdV2b3z8uvEIt3b89DeXdGWzncWaiyX9mUdKF+fksYvhfusGnwHuqwVXXaEymuUJZZ518ak7Jxn4C4fvyRVUgT0QhPOvkUJw1t+69baOTrPhj+ZKApbiZKLbvCuqyqnY0pddIy2M053jFxUE+FA5xraZCkf1qAnBJk3KUt21Yyqpm62UApWUqzDM3RuL0q0T6cbkt8wm9D/ZOwP0XHznYU5nrlEbMx/zYWGwrcPkQepEw9WGRr3vkWTtEiTXkU97Bnolbd7q4I3WemKhGmBgJVbvc9BuTR3nqKaCAc6xtomS/mwIZeEIEuq9R/YyI2BupNCoiQbgLo2bKZhD9Ey3nnZlvevqD9+qvJimsY2TTo18ZbrDjWyZrQktsd/zYWGgqGuQ9LpXae02/a5P5cZGvd1HR4sBiroUhk7lWStTEq4Z6BXKFtrBFZdzVRh2p1LZaxtf1ExrIceMzzVhQJP0F0hyAtF7LvtKS3JW3cLN9NiHVkp/cWDxfrjc8L0OsUrJp2awsiv8l56ndWLbAWhKlaSbeOeIJ1jBAGlSFo9qfK0VVLCRvURmqsj2T0fHix6wjFtXXnc/uAd9Xz+6rGiUs2cr7DyuTFCJLVq+mCqpFU9CCHPy4O0o44AjzLx8UX7ozTVmBwyaF2LqpytXwdEFFpyG1uhZEDV6OvE/1XbfGXJl0Px/XVkZGXXycml8MZmKm+rI+sruqZtXd7iJNG+rRH3Qiafaz30GMGN2TJzDWsTrykSrRJ34tTQ8IjKxB2Ikm8OMwWQV1JqXYvqkn8d6o3ryq29Fd1ygUFTrSB0irtk2wwPFkM39HComc6xTuw4MicMliKkMuLOScBcPoaXv/CsikbGlEl73ne6qxYCSSUnnPxSI3qqtiqsQY8RQRtB+8F7IdK5NGYZE5bUm8LvmQWVAVAl1GQhAZ2QjW6FrbOc94dPVIqFQXq9miAqup0zAf7t13/oTaLOMsxO6a+6eRWyvAlJd6XW1pXHPXvv8/xexH4SiXL5x+dcM9s0Wg2bFI0RURXa8JKW65/7GO7/+scj08c+tfOY/jJb0WtSBVnyTTXZ+a+FaEXQvqzDOBkHmPd6lTWE4BWCRUm3kzFiTOFosewhftGabn/Q3OIbuM+gwzF3qoUBaHPbHaG3JHuqtiqshx4jdJJouhAt46PwToYHi9oeuRPP5Cn7mZxbGE/YfS1klEuRBygzqiYKkE5DCAeO1PHLh/4Gd2z5kCex53iUIvVJoLGtAAFwtVhEoRgRdCYsneIrN9qXdUS22p3vsEnRmNGouJ/JcU2W059784uBjsEbnypBqQsZOyhoAw0dPJnZJTZKAgZIvquAcqnMPW8A3CR69//2Plw5/ZPkjb0gwS5L6joNWWR5CpNEu3NtRBNJq/VUjQKypKg16PMQpuwaXRZNvquAz775hUjHGedkF2TSMBnTHnrMfFDVhtGyBhUiyd9jn/pWZLr1JvA/M6bn7b/mus+b+7gndhzhrggXYp9Ra9BbBFEZOB0PyJ280vaYIqZEJgGTa2o6AUg9dAGCepTDg0UM/d6LQBSFP0FQdQiC6ubr0BhF2/OMuYiauhBgaYsNQhhjEiaLrxNXdJdf37y8U8ugt2K5tUkDDV29HAe92+70xNBVCJo/qRUiNcqYAzVjGpT95H4m1+7uw0uf/KZYRI7g6SzFvcYuaqrFHCzLJSaYanVHmcXXNbxOw4WR76qZEguh3No08bZu/0asemiNXJbBAUffRBemWuTNCPcz2TPQi7b2NuG2q7avAVDx5GViYzYhWg9r0GOCykD7KW1RdqnXpZgB1c4zMscvAkpkqyAIzXTd/o14dPYJbHh+k5xK6dM3MUHLGy6Xx+3A3cCcBy1p6RZcMcYNG3KJCSr9Cn94RcSKUGl3OBBWRCoUCWVoJgZBEmyhMDRTd2hHlDQMIsYFBOg/22RYtX2NtpImgLmqVwXm+4oxCKyHHhNk3p5IEoC3dJc1a3DAC+9cPFjE2t19eJQ9gQ2HNgnVHIWdbjheVaMQtNWcKVSqk6Kx+YuHZPIHqjHz9sfritRK4CUuZc+WjjHPLMrO+xVjEFiWS0yQMSaGthwWesw6Ykd+6PKteV4uwOc+NxODIG4+eVCI7rEs3i1jzfD2x5N4aCXI7lEYkbINz29asAbdslwaAFHFIVAJo4iMtk5izm+YdePvMsZHM4seNWuVYJAuVG7WjP8+Tl1/p25/jeCdR4VMIYvl/Su4TVAA4J699xkpOwKoORvN9Hw2E6xBjxF+AyoreXZitcKKOFfbMjflSxZb1U0amVD7GoG4BbPcMInVB+1CJcyjzCO0d3dgef+KOvmDoc2HcWrnMfTtXV/n9KjyPFZGVw1r0BOEiH5GafIsw2WJuZM7j2otv8MqLzaTxx6lJo4MprUAsoYZa3f3CZtMCPMo8wiy8y+NTnD7ozZraK2VoJUUJaJ7iejviOgyEX2e8/0yIvorIioS0UUi6o9+qK0PoUc3yzwiSLLEnCzmGIXyYlIJSOdY/gQgD0GSlUFgWgsgU2bsGejF/QcfEH4/3zxyNzKLskrxLd515dJtqfIMyp4PizkoPXQiSgPYB6APwFUA54hoiDH2imuzLwB4gTH2DBHdDuA4gNtiGG9LQzd0EDQEEoUXY1otGRSm3rDJNQm6wjCN1auUGWXfi7xXShFuvPWmljX4lE0h25bVqijlXdds3heqkihAWtRDJ+SyBsBlxtirAEBEfwHgowDcBp0BaK/++yYA/xjlIOcLoggdiLQ0hPRDQySVgIxr4uBNFLy4LQ9BYvWqiUb0vch7ZbMMy/tXGEkKNAscfSCdxhWA97rqtL6TJZQbHRZsFuiEXG4B8FPX31ern7nxJwA2E9FVVLzzf8fbERFtI6LzRHT+jTfeCDDc1kYUoYO+veuRzqU9n6VzafTtXR/JGKNqyqFCXBOHKDbtxG1ly3bT5hYyqMJJMq76hWejNeZRTfYybHh+Ez735hfRM9Cr9az4r6tuTsGdUHaHBYe2HPY0F1moiCop+rsA/owx9h+I6LcAHCKiDzLGPNk7xtgBAAeACg89omO3FMIySkwaMASBqN3d1PV3sCf1mPR4Jl5TXMwV2YQg8vDaOvMgEErXJtDWmUc2n0Xp2kTga6sTTpIWz0T8Zrglj/0sqajgvkaqNoGUImTyGQxtmWt9pzuRywrzLjx7Fkvv6l7QnrqOh/46gFtdfy+tfubGpwG8AACMse8DaAOwJIoBWtTD384rygfYv4rIdxUwy1hN80WUJDVNpkbpDbvR1pmXfs/z8CZHS5UwVvXf5VIZGw5tCnxtee38/ElAnVZ9UYDS5FklqISxgsJ9n3sGerFyq/i6sVnv8/TSJ7+JfKd6FVFLKIuMP4uuZ2urQsegnwPwfiJ6LxHlAHwCwJBvmysA7gYAIvoAKgZ94cVU5gncE0Z2ca7Om+MxFEwZInExV0ghfahDGdRVuXSHVZ5e8ji+suTL2EOPCROCfglZXQG1MGAzrG6CVQljAQA0+306OPapb3mM+isvvKw/xvIsypPlulAiZVJo68rXPR+yVVyji80aDWXIhTE2TUR/COAEgDSA5xhjw0T0JQDnGWNDAP53AP+ZiD6HyoLx91mjNAUsIoVurDtITDyOgiaZsVL1GnVDZRj8YRUdVodfQhaoePNBNcZN4UxUSrEvAjZ84wGjDkUzUzOecJZpSf/02+W65uNsehaTo6W6gqLaPeRYmIWuwKjFQ2eMHWeMrWCM/QvG2O7qZ39cNeZgjL3CGLuLMfbrjLEPMcZOxjnoVoYu97pZoJskTSqZqoLoeO7iLZ0xibZx7t/Q5sPGhUH+cFLPQC8+++YXEklaOhi/MqZcHWQKWQxtOawlNubfNxA87CGK6/vDdz0DvRXNdN+wFoJmvwpWnCtBmPb6NNmvSWckk4Sqbls20bklLfKlM14VRU4koKVDrRPBofS5PXLns6tnRiKnKapE3mrPwchYbdt8VwGT45OBE6b5rgKyi3OxcejbuvLILb6h9uwu718hbUY9X2F7ijYJVL0UgzS9NemDGaRpsvM7nUmgWRr56oxXxHKRnZ9Jt3o3MoUsVm7txd/+lx/WiW1RNoVMLo3y29HJADjHc+uoOMh3FYR8/KDn54CyqcCTQWZRFtOG18A2ieZ8Zw16ctDpdm6qW6HSv3AbLpnKYxRVps2kxRFH4Ylut3oAdSsvkehaVKA0gc0yz7kODxZxcufRung2zxCe2HGkocVM+a4CyqWy8erH77UvBC/dyuc2CXQ6z5hm6U06I4m4z6J9mBrFZpG5jbLhthuy+5ddnEP6hjQmr5W410q3ejIo2CzDo7NPeD7rGejF6V2n6gy6vyK30cYcqCSzV21fo92tyMHkaKl2flYewHYsShQ6VDXTJKJxZyTNfQQR6WqWxGiUDbfdkN0/Nstwz977sOHQJgDwJBUBNT9eC1QpyuFBdI11JtnigXNah2/v7qjQCAVj436smVBt68zj4sGi15jT3HGzi/QonlHc51aGNegJwsO9BiLJ0ssKdHQ8Y9ExgxjFuIqFeJCxhYKsFHTYR8794xmp6YkyTu48KpwEVfx4HbQv60DvZ+40usY6k6xWy7fqMe7Ze18dvZCyKazavqZuXJRNIZNXG+JMIQsCcas/nZDVrEFsfiFz0W3IJWG4uddRxHllUgCiuC0v3upHUF65aCxRQhVS0ZEVcF/7fKeX3eHfH48RwgOPe+1MglFwzZ1esbUG4BrXWEcQTnZOzvcrt3prBnj3eOld3Z5Ec/mtKZSvTwl2Co/BFoWkxq+M4fSuU9zOTZSqPMd+LGQuuk2KJoRGqMMFZbUAzZXg9EMnESw7b136oWNsglIVayBx28EgME0Eqp49nRi6KaNEhzHjfpZk91TWzcjfGWohMF9kSVEbckkASTaNcCNMeX2S4RNTqFYPqvM2UfaLorNQ+7KOyIw5UFkJmDxLKu2fdfs3IqOIUZvGpnXCHs42w4NFridP2VTlc8Glc+5r3I1PWgk25JIAkmoawUPQ8vqkwic8qDxKnZCK7LxNlP3CxmNrRilGhH2WhgeLyLapeeAm10KH0dW+rEO4WqJcCmxqVhiqcneGWsgG3A9r0BOA6MF2Wms1K4c2qZfFX+RTfmuqFjPlNai4eXkn95ou71+hdTwdY6PyDnWRIkpEqyXoxGNS/WoSm9YJVTn3lgc2JU6C2mbRYliDHjNUYRXHsMTBoW2Fri5+gyISdXIaVFw9M4KR777K3eby8UtYp3FMkeZ79sYcJq+VapNKWENMaeIm80yw6qE1ngRo+Tp/XEETgTypXwC1pKUDf7hN9Wz5V3i1atywkxtF02pxvsIa9BjhGCtdRBmGkTFBgMaEUngwiVFPT5QrnGmB16zrparCSftue8pYLdBf9u5P1hnDpYPjnqRECd+gXZWEBpbNJST910e3cIu3wgsrL7CQGSw6sAY9RgRJqPGMUhBPWxS3P7nzKKZL05FXUQYdr2moQEqvS5Gyq5KDoDH2OkNdNbxuyp5zfFEjaB3kOwtYelc3d9xANBOyLMkpYzOFyQmFzUk0Q1K+mWENeoyQPry+Ja0DvwcStIxddGwZVzrulQFv/zrxbDdknGnn87CTlDDp6tZl4RhT/7FE8WEdOCEm3n6jym3Ink+Z4Qwj8aB7v9O5dH24ynLylLCXKEbItLl5lXW8pXPQMnbTpWlU1XWylQEPPHokZVPcUu90Lo10m54PEqYEXETZXN6/os6YAxBWmIZtMxd3GbvoGWnrynPVKZ3zNJUfcGPt7r66StO6/XR3IHtjrv6LWdtiTgVr0GOEyDDcf/ABrNu/UYtDK2TIKAywaYuzqGKTspWBrJzefR3u//rH8UfX/y9seH5TXW9TP7VO1hwi6CTFG9PKrb0ofu28h//90tYX8dInvynkhK/d3SfUONFFnGXsoufznr33eT7z11HwVki6cfyegV7c//WPCzVhnP1MXuPnMBZyWb8ObKVozAjDNBkeLIpbbWlUbA4PFvWW/RE2opAlvcJUmcr2q2rmEAW+suTL2olS93F12rfp7isO6DyfomuvIyER9PjNXKncaFj53AYiTLzz9K5TfEYH6SWHZHouHjDg4sEilt7VHTo2u3Z3n3ASCeNdyX4bxmPUgWmPTKe+YO3uvgpTJCCrI4nKXJ3nU3TteZK9buhMFqLjyzRohgeL3M5PzUbJbQRsyKWJITRiTC/ZNzxYxNT1d+q/4IQBoorX9gz0CsMgsrCOSu1QJyREaaqFZzL5TJ2EbRCYUk8dOEUzpTffBmXMX7O4y9hNetsGkUUOK3chkm8AgGOf+paHbjk5WsJLn/xm0/fnTQLWoDcxhC+SRrLNeaH8nmW+qxCax61C3971RjowOi+/ThUom2XYcGgTyqVy5bwVhkTHqIXVcim/XQabNmvL5oQV4jTmJsY2iK5PGE36WiPuqgLjhkObatdDpLzIyrM2YQpr0JsaYQSyRIYouzgnnBCiSoyqxLFqhpQew5OZXRjafFj58l8+fkl5XFFTD54hGR4s1iU0hzYfxtNLHtfSVo8LSYRZTI1tEJE3IbVxZEzqSasmG1n4yiZMbQy9aSCLNwZJqsq4whsObYqs2lAEUWxUty2eM1bev3lwxi/T1Xbj5M6j3IbGpdEJDG05jKtnRrBu/0ZjnnwYJKVRIjO2osIsno7/0JbDwmdSdt1kNQLKyUZQv+Ecc6HDGvQmgKoYJ9LimGUdDVVSNAlhuF9QmYFwG0JREtj/skuTnAy48OxZLL2r20wPXWJsZIhawzuoWiUAj1cM1Btd3cIx2XWTFbLJHBEhSQCV2gVbRWpDLk2BOHpgqsI1bo1sp/pRJ0EWFrrLYjejoUZh47Ts2/D8Jk+8WXbe7pi5Egw1o+MON8jg9BNVIZ1LV3jYLqeALbEAACAASURBVI57VNdfJz6uU6Mgev50n9WegV6s3GoekhEW46VIulK6/+sftywXWA89UfBamclobWFigrpeeFBpgaDQ0snunqvC9Hh5rH4bHg0OqD/vq2dGcOHZs0YetLthRm2forG7GhqbcLajvv46Oiv+ayRMknPknU3K/mV5D5HhFnn2bIaJ5TK6O6wxr8Ia9IQgih3XPE8NXRdT6IRrkm6+IVuK+0MP+257ShjqWN6/Qjg+/3kPDxaNjTkAgAF//uGv4ed/809q2deqRy/iT4tCKlFff12D675GwqItqpd3buvMc8NVvGc1iFaMMyausJmgJkNXB38hwIZcEoI0dszADSckERMMI7QUBJ4QBqrccfCZE7IxXHj2rHZoQhZ7VWHkO69qa3iPj4wZM0Kivv5BOOPCEIzvmk1PlEEgbeaV6Jj5roJSCZPX/JmLalGc5aBXYD30hKB8QSX603FCp51blPAk7BSsDlXyTubFuo8TtuuQLpzJyR/ScMeX/aGgqK+/rMJShLqQkiS5W7o2gQ2HNmkl1EVj6du7XnkeJuyipNo5tgKslktCUAn7N0qjQtQwIY4qRdNjybRsHPgnBX9ZeNJ4lD3BPU/KppAibwejTCGLlVt7cfFgMdLrH0Y/KOrnNOhYTFrjAQAIUhmC+QSZlos16AlB+oBWPaJG9UpMqlVdEMGlEzuO4MIzZ6X7dQwgUCkL1277lgJgVsQphZP4pJRYs90PlcZ60tiTekw4gcY10Yvgfi5V13QhiXZZca4mgH9ZW1MIdC1v42aYAGLjncRLGiRe7ChAypKazpK7fH1KbcwJXuaLYrIwgWNwTLoUjV8Zi/T6yyZn0Xc6hpPSlKgxB+qLmWQOkU2MVmANeoLwx1Z5sq9xxgOjpMgF8eqDxovX7d841+JNRvFU2FG/F9cM2h9R5ipUfWR53109M+IJ+YiUK6Mw5kGeGfdv8p0V0bU6lk2EaqGtDhtySRBGcUGXJxnVQxqVxnTQuLvu72QvvuwcpEk0qhT+eDzShEr6a0PIpDxCXVGHMGTXBuDroIi05KPQOncjyDMj+g3XqGPhhF1kIRct2iIR3UtEf0dEl4no84JtHiSiV4homIj+PMyA5yuMlPsCSI4CcgXBqChyQStb/ZQ+nsytqtJRVgkq6oIDoMaKObHjyNz+owABGU67PB5uuOkGI4ErU8jur1DTXNSftap1HpXqY5BnRvQbkWyDFefSCLkQURrAPgB9AK4COEdEQ4yxV1zbvB/AowDuYoz9gojeHdeAWxlBHjh3CEa1ZD2x44gn1uwPqURFkQszMThxUVF4IJvPSgtteLkI5/vbH7xDGhMfHxmLNGYOoNKSrTwLyqa4Yl9uTI6W8Lk3v6jcZdAkter+mnjoUdNWgzwzpu+LFefS89DXALjMGHuVMTYF4C8AfNS3zb8FsI8x9gsAYIz9PNphzg/ImkbL4EiOyjxXUTWk2wsKI8ercx4mL5TI+xLRDd0vd89Ab+1c3BW3Fw8WQbnka+VmpmbQ1t6m1HpR3WcgXGMI2f0Vfde77c5IngkVgjwzssKkJMbcitB5+m8B8FPX31ern7mxAsAKIjpDRD8gont5OyKibUR0nojOv/HGG8FG3MLgVuQRsGzte6ViSZQm5ZJVVg3p1iQRVTGadLDRFcCS7cd4tcLg2Z/oerCpCHmIBiiNTiiTsjrsl6DhLMern54oc6tvRfdet1l5WARxJkS/6du7PpExtyKiYrlkALwfwFoASwF8j4juYIx53lrG2AEAB4BKUjSiY7cMegZ660WiGDDy3VelxoDNMOWSVWYg3Z4OjyIXhP2SyWdq2+e7CrXqP939iMIDbV15TJemubkG9/6aLl6qIZ2r02kqSGiCpxPkGEuRprkbSdBWg0g2q35jDXg9dAz66wBudf29tPqZG1cB/DVjrAzgJ0R0CRUDfy6SUc4jXD5+qf7FVxgCShNuXHoTP5FX9VxFokk6DaVNBKJ4zINyqWy8H1FZ+D1776vti3e+0xNlnNp5zKh4JxEohqIbEtDNc6i4481YDh9k4kiqRmK+QCfkcg7A+4novUSUA/AJAEO+bY6g4p2DiJagEoJ5NcJxzhsE8SzZDJNqWI+PjKH81hQo67udBKzavkb5Qph4hTKjbdJ2TBb+cbTaRTHp0uhEcxlzBfJdBe2QgE5owh9nF12LplvFWMQOpYfOGJsmoj8EcAJAGsBzjLFhIvoSgPOMsaHqd/cQ0SsAZgD8H4yx0TgH3qoI0tKsrSuv1OSemZpBvquA7OJcpOwIP+NCVthj2nZM5X0l2f4tFqQIG77xgFCnRlaxKwtN6NJfw7I+kpKEsIgOtrAoIXiKWQxblaVzaax/7mO1l2kPiTvuBNGD4RY8OWP0j1XSZEDVrs208MNPw2w1+O+bA971TufSyN6Yw+S1ktJ4yvRWHEQh8JWUaJvueOzkUkHowiKLcPAskQGu/rkMM1MzNZbD8GBR+ltdmpubjXJ61yms3No7l7RzG21evF+g3e6EUYRjMwgBDA8WcfFgsWWNOeC9b27wPOyZqZlKDkSDqiilv0bE+oijLWJQhKFyLjRYDz0BqCRJddHe3YGp6+/IGxy7thV5wzLvS7cknqfdriqpN1EWjOqaNQP8qyYdD9v5He8eJuE9y8boyCy42yjG6TFHJVkxX2A99AYjquTU+MiYljF3tt1Dj+HpJY/XeTJBEptu8Ers61YhPmQKWSzvX6HtaZkY83QujVUPrdGiBTYC/vPUjW2L7oVpV6QgEI7R1ZbOXdQVp8ecdFetVoZVW0wAouRevquAcqmsr+8SAKXRCRz71LcAzCUkZS+IKhFJ2RTKb03VJpbxkbGKBvn0LCBoG+b2zHkTycmdR+vEuUzyDDNTM3jlhZdBJnGshOGmEapyDQ5khj8OOp87Tt3WmUc6l/bKEUvuSZw0yaS7arUyrIeeAHQr3vJdhViO74/lCmOwKcLy/hVCeiSlCRn/S17dv8iYg1ATeBJNJJOjJfz5h7+GJzO7sIcew9BmeZci0T4a1aVIF+7zz+TnfKns4lwd5VSHt25S3auCP049OVrCLGOVZ7K6ClDdk7g85qgkKxYCrEFPADqc60dnn8Bn3/xCbGNwv2xrd/chnUvXbcNmGC4eLOKW37qVqzvCZhjKb5utJtyTh8yjGvnOqy3FLQ8ChwZ6fNsRT+iMzTL0/sFqoxBK1IlC3uqJlWeRXZyrqS6qQlpxecxJhJjmC2zIJSHoLpHzXYVYPM2sz8OZFSTDpyfKSikCbfg6yazd3Vfxvuc52rryeGf8HY/6ImVT0rDT5eOXjBJ8JlW5OtCJU8tCRc75xQVbMaoH66E3EYYHiyhdiydsUH67jBM7KsU9p3edkku9RuUoVzvJOF5jz0BvbGGlZkGmkMXtD96BFHlXOCkiXD0zIi3MMgmhRJ0o1FFDdDxl3j30n29SiDLsNB9gDXqTwFlCx8m7vvDMWTyZ2ZUoHdDPXe7bu16qLNnKcPpuXj5+iZtnkGmx5zsLRiGUKCSM3dCNU/cM9CK7OFf3e3+tRBJG1vLT62ENeoPgf+hP7jwaK9vFgTJOHYOj5dcy98dDu+9+X/QHbQDYLJMmf0VI59JgYEaFPFEnCk3i1LLVQVJGdniwiJe2vtg0xU/NAhtDbwB4crXNgq4PvAujP34j0pWC32v0d3OfLy+gc56mGjSzjBm3VQsiR6uCbpxaRiOMOrbPg/P+WFGyeliD3gAY9RZNGKOvRNt4xJ0s8+txLO9f4ek438pwn6dp8peVZwO1gguaKAyriyKSPl67uw9DW/jnHaWRVb0/C5mfbkMuDcBC8iDa2ts8PUTdS/ELz56dF8YcBNz/9Y97Gi+YJn+dphRuxMG1jiIkIgvPRB3b50H2/ix0fro16A1AvpP/smcX5+YdC6Q0OoF9tz2Foc2H6413M9PODXMJfg+Xl/xN59LC/TpGMW6udVSiW+76CadwDEimCEgmTrbQ+ek25NIAMJElY6zW/Wc+oZlyBLrIFLKY1i2iqnaN4rVH84eYXnnh5bp4uVutMm5jJGtCsu+2p0LH4+OI7fshCvksdGMOWLXFhkBXbc8iAhhqz4eBzKhwNecx14/V1BAFjYPrqlg2u4FcyProMrVFa9AbgFaShs13FVC6NtF8PTwbCckk4Zd0PbHjCIoHzgmvXRAJ2DDyuaKJJaqxWcQPK5/bZJD1B2029O1dj0dnn8D9Bx9omTHHiUwhiw2HNglj4e6QxokdR3DhmbPSiTDIxB4mDu5JaCqQdPLeVn2GhzXoDYDJSxUXKJsCUurMn2MkZGXfcYGnu95ImLI5igfOKffJE0FTIWzZv6oJt4Mk6X+m7Btr/PmwBr1B0H2pokJmUdYjhXr/1z+ODd94QHl8f5XnZ9/8AjY8vynQZNTe3THXiKI6DpHRbu/uQG7xDcbHiA0ET5xWh82hE6IKEsaKihoo2z5p+p/JqsOW/IthWS4NRtSd7TOFLFZu7cXl45e0EkaqAhjeS++wMUxyAe547DrX56IG1U2XY2DwNOLwsznynQUwMAxtOYzTu05h7e4+YbGQG0E8dFlhT9j9AJWV0T1770s0yWiy6kiiGrVVYQ16g6DqvxkE2UVZfOSrlZDAOvXmAOb6Q/KgMhK1ykANJ5P3Yg4PFnFq57GW4adPjpbw9JLHUbo2UaMhAgAYKiqZ1XE7HuOyte/FyHdele4ziIceFTUwCYqhLky6EtmWdGJYg94AmDANTJBfssj4ZRR5aTp0up6BXlw9M4ILz55VGmH/izk8WKy0rvOpEjY7HK368ZExr3qi7/ynJ8r4xeVrWPXQGiXLJQii4qzHyX03oRaarDpsSzoxrEFvAOLScgnioYT10tbt34ild3XXVhuUIjBfOzrKplC+PoU9qcfQ1pkHgZq+XVwUGB8Zw+Xjl8BmGdq68ii/NeWZwOZzmTpPgO74tooeP+/ZMnkOowo5zUdYg54A/J5KXPFhShGGB4uBlt5RFrb4mw2X35qqGXCRquC8hCsXMDlaAmVTNV7/fC+GiTPO3UyhomaDLSyKGXGFV0RIosLPpLCllYqoksB8LtZxT+TCEBwBj84+wf2tLefXgy0saiCSlspNQuDfhGIWS6KqMd3OIsF8Tdz5qYQiiOLcomfqpa0vWjqiAaxBjxmNeIHjPqYJyyCWRFUTsWBIozjLv/18NFA6josszi16ptgMsxxzA1iDHjPaOpOvdmxf1oETO47gycwu7KHH8GRmV61BdFT754F3rq0kcxAEpiHL+WqgpE6EhhywbOJf6G3lTGANeow4seNI8klAAm5e3unREGEzDBeeORuZUV+7u68iHeBD+a2pOkPlb4YQRTl/s2jGU5oCrUCaxUBFWT4vrF7t7qjTTOdBNfHP11BV1LAGPSYMDxYr/OwkQcCq7Wtw5fRPuF/raIvooGegF23tbXWfuzu/A3MGw2lLtuHQJnzuzS+G0rDJFLLo27u+KXRe2AyrxIwDxPQbbaCiLp8P29jCmfhFlbOWY64Ha9BjwuldpxKP9W44tAnr9m8UFrGwGRaZmFHpGp9H7hgqmcEIGobJdxVqy/Z79t7XPKEcBmOj3mgDFbZzkd+7BxC641LPQC9X1TOdS2Pq+jtWiEsDloceExrhgTkvo1RDhNUXeQRpFqCq1pMZDIe2Z9JIObsoi8+++YXa3874TPYRK0wmb8KcbEAIuO+boyUzea2kdQ/DlM+Liob6D2wMTcnkaeRMjk/WQpeqAqWFDuuhx4RGeGBuDREVHOPK86SHNh/G00sel3pCqiW2rNWZU/ykG3qhbAof+erGus9N9tFUYMDFg8VQnqb/vpVGJypGTzN8EkaxMaq+pCK4+5VmF+fAyrN1x3pp64s4seOIldD1wRr0mNAodsf0RBn//Dc/w6qH1iiV/MavjAnpZqXRCalRkHV+B+SGQSv0Uh26I/Uraus2df0d6Tn6kcrp6cAHAoGbLOZheqKMkzuPBjZIKpqgysCGiXknKY4lozNeeOasldD1QevpI6J7iejviOgyEX1est3HiYgREbeKaSGB1xAiqUTe5GgJS+/qxuend2PD8+LuOu3LOqQvocooiDq/A/IJzV0CLkyEsbmqSlmPzjoWUYrQfff7hGOenZrFozO78Sh7AtlFEU+4DEiR/mQxOVoKbJB0jKdsG9WELENUeuw6MNlns7CHGgmlQSeiNIB9AD4C4HYAv0tEt3O2uxHATgB/HfUgWxnl0pwXFTWFURZucCrshMnZasMG1QszfmUsEL3NMRiy/Trb+cW8/NvwIPRQZxmu/n8jwt+5J4+PfHUj0rm0cNsgCKMeaWKQdAydahvZhCxDWEYLD6JnzHSl22j2UKOh46GvAXCZMfYqY2wKwF8A+Chnuy8D+FMAkxGOr6URd9m/7AVyCliEOiqs8kKrXphsIRuO3iZwWN1FSEE8PtmLKzOqbIbVDEbPQC/WP/exporDqwySY/hUdMk41QfDePc8yBhRKjqjH41mDzUaOgb9FgA/df19tfpZDUS0CsCtjLFjsh0R0TYiOk9E59944w3jwbYakvAWZEU20xNlMa+3asRqnrQgrlwuTQdKgDkvqYj9QS5rFMTjC/Pijo+M4dinvlUzGA+/9kigzkFB0d7dIbxvMmkAj+EDPHTJfFehEtKLwMDqwLluGw5tAgAMbTkcODHJa3LifsZEdEYeomAPtTJC0xaJKAXgPwL4fdW2jLEDAA4AFbXFsMdudsQplQtUWqLd/uAd3kYLPrAZhkwhq9aOFoQ9RJ+rJivV6sTNYw8ihypqzKGLmakZnNp5rHaMIJ2DgsDJC4hUOJ2VFVBPy+NeU1euISl4um0R6jo1AfqUwuHBolAbf3xkDF9Z8uVaqDK7OIe2rrw0dHnxYBFL7+pesJRGHYP+OoBbXX8vrX7m4EYAHwRwmioJofcAGCKiDYyx+a+PK8Ha3X146ZPfrKNdRYXJ0RIuH78k3YbShP4DG6Xa5c5LKPo9z9jpxN5l8P/eVJPd2fbUzmOBm2W4fydrxSeEy5jpwslJyCY8kW64LrskSF2BLuomIk6nppM7j2ofX7XScxvv8vUplK8D3Xe/D69//6fc67fQe4vqGPRzAN5PRO9FxZB/AsC/cb5kjP0SwBLnbyI6DeCPFroxr8HwhXe8rT30mNb2KiPEZpjHWDov+9CWw2hf1oHy9SmhYXEaTl88WDTuDiNbnQSJ74qMlKcwKsBqaE/qsVp/UP95Kps8B3Dq0zdktIqhxkfGamNzzlWn9ZpppyBT6OSFJkdL2oVAQcKSqj6tCzkxqoyhM8amAfwhgBMAfgzgBcbYMBF9iYg2xD3AVsbpXafAps28c8fQ6STqsouyypJz9354ySeZd9t/YCPW7d8YKAEmEvACgJVbzbxxle6IE88VXjPZNarur3jgnEcKN7Moi1Q6+jKNmclp/Y01JBP8k2PcRT9BjKVTCOQof+6hOSZLo6mO8w1aTyxj7DhjbAVj7F8wxnZXP/tjxtgQZ9u11juvwPjhp8oLuSf1mFbBTKYtK/USdV52Edq7O2pGNwi9TSTgBUAZJvJD10iJDN6q7WuUBT9shqF8fWpu/2+Xm6aBtZ+3L5tc4y76CWos3cqfwNxEtbx/hXYxlg4Wem9Rq+USI4yTor4elEgBkDj4Mu+6vbs+dqn9UkekNaIS8NKFrpGSJVfdjaybEZSuNtcWTNBu3r5sQtUJy4QBNxldzSW0d1dCeLo5jemJMi4fv4TMDRmUy1PqH2hgobess6X/MWLt7j5u4QplU+i++301qhylCZRL1RtvVbRGVAEqqLDUfqkj0BoBxM09eJ/LipfynXyKH+98RKuJnoHeyiTVpO3r7j/4AB6dfUIYNtK9d2GLflRFZD0DvVi5tdfz7K7avgaPssr17tu73qwQaGTMszIKA/eqcqHCeugxgsfEaOvK456993kevOHBYjDVQIeH7PLqZC/v2t192seJgi0gihf7P5cl8gBgcry+Vi2dSxstrWv69CHYiZlCFkRA+e1oi8W6735f7TrzPGBTXXFATAGVMWB0EqrDg0VcPFisa55y4ZmztVXhyq292tdamXjWxEIPtTiwBj0BZBfnULo2IX65JGGAtq48pjnFPTVUl7o6FLGegV4jmp87pBGECicyfOW3yx4GhyhGfnLnUbwzNsl94bM35oykf4Po02cWZVFYssizbwA49qlveeLr6Vw6VLy9a8WSSvVn9Tgrt/bilRdenuNf5/U93jAGW5arcE8UomfR2V8mn9G71hQN/5/nJC1UWIMeI0Qv0NUzI/jxCz/SMqz37L0PQEWbhcsHFxSVuF/sts48CITStQnkOwugbMrLjRfwqZ1lfixUOBeDQ2QgZAUkk9dKRuMKkhTs/6o4Hus3mmHi8+7CsPGRscrfrmBoaXQCQ1sO4+qZEazbX6+Po1voozLYOrkK1XWcnijrF3tFVMtFzRpHawCsQY8RohdIdzma7yp4DIpoKe73yvycardhLI1OIJ1LI9uVrzVD4HGw3UtYHc9NNH7VpOXIE5h6au3LOozGJUtQ8zzs7OKc8NiixGSYytU6+PMnDLjw7Nm6KkidQh9dg62TUI27+jkIHKlnwDa9sEnRGCH0ZjRsl9M7092XM5PPVDRAqpS1lVt78e3PHMHQ5sMejvaFZ89KDcvM1AwIVEscqrjmQalwfXvXa6kZshlmlKx0JhuTcYlEyLKLc/j1T/9GnbRx+fqUVIRM1IItVk0YVlmpuROWOlRUt8Hmwfmce418jKe1u/uaMrHshOgWOoixxkiqrF69mp0/P7/p6jVVPEM4ySUAddIBlE3h/q9/nPudKTY8v0nLoxGdh46GiHv1QCmFJ65Bf6M04f6DD6BnoFc6rloYxBf/PrnzaF0oJ1PIIpvPco+X7yoguzhXtx/easmRC47UU5fAr9Ejgkw/xhm38xyc2HGkbgXJ3UaiH1RDAGmEsNB9plsZRPRDxhi354T10GNEkK5F+a5CjW53cufROoPNyrM1rYywGjG61YNBqXD+UFDvtjvl18MlNMWjv2UK2Zoxl41ref8KbmUpAOQW31B32OmJsjA0VBqdqNvPyZ1HpaEe92onTsjUNB2475NOYdLl45e4YZuhzXNqiuv2b5SqfNZgYMyj0qVf6A0ubAw9RvQM9OLqmRE9b6YKdzGOKCk4OVrC5DVFswwN74gXmpCxJExYLryE5cWDRazc2ovLxy8JVy7O5zrHFG0ji62HrZiUJf14xT+6mjxBwWZYfQ7AtdLhXbMgzaMB78TYt3d9qBXiqofWVJ4D130LI7RWG+MC1nEBrEGPHaZl7o4etmrZKE1OEbBq+9wLA0DKYnGgYo2YLGVFRvXy8Ut4+LVH8GRmFzf84vY4dY7J22ZoC59r7zB+eBNlZlFFRiFMqIQXow6k4miIWcYqCegqNXZ5/4ravXdriutAlfR0JsaHX3sklAFet38j1nE+DxuuWsg6LoANucSK4cGi8cvMZhiGNh/G/7P4/xQyLfJdBan41arta7Bu/8ZaxeSGQ5u0QiZRCjupEpaiWHoUvGRZ8k9Eccu2ZevCEaIesPmugnYISlQtHCVYeRal0QlQinDz8k5cPFj0hImOfepb+MqSL2u1ENQJEzr3UCTtoII0TBQyTLXQi4ushx4TVDrjKpTfLgMpgDKpOsXGDzz4QX44p+qZL72r21OosnZ3n1QT3UGUwk4yCtzwYFHMfeeUvpsUNQ0PFrnCZk5sXRT+Kl2bqPP2RUnEvr3rAeiFoJzPAlUCG4LNMK607MzUDGY05Gx1WTOOdENQCmPvtjvrjunn0FsEg2W5xISgDBctSB58yqQqFXiu2KafpSBCGDaLH6KOPG1d1SIn3lKdgA2HNmkZVd75iI6Z7yrgAw9+sI5rr3OOUTWLiDuWHgTuGLvo2vGQzqWx/rmPATAPkWQXZfGRr26MzYgn3b2pEZCxXKyHHhNijZvKei5w9Nd1dVnC6oi4IeooJKv+dBpXA3JZBNH5iLzL7OIcLh+/JG3kITpHWXMQHePu/KYZ4fbWTaSVZ6ZmanF057fuojZZ0rv8dllaCBUWNilqEQuiEh2KCuMjY8pkaxA2iwxOublu4swJt+h4i7wXN2jISGf1EkT+wMTrbRSCsn9kcr7rIF7tUZpivR4LPSlqDXpMaCZj7uDYp74FQM54MGWzqKBrKFRSA35wGSWK0nVROEnnfIPIH5h4vXGiTrvHBxn7RwQeQ0omP6Ezjiiw0JOiluUSE7QKLxKG0+k+SYg8prauvFhqQBGukjFKROyTsDrhQbz/Zln+3//1j0tbGsrYPwDqmCfu6zY8WMTTSx6vk59wag6ce9zWlUeK4q20yi7KzvsqURWsQY8JLEhwkCoTQZx6IGELN0whMqT37L2P24iixoARIN9VEIZIZJWQOlWSMqh0UEy/SwqOwNvDrz2CDc/X01cByLsMVRPVvOvmhJR4v3XXHDw6+wRyi28QSwxX73d7d0coR6hcmg7dlKXVYUMuMUFZyelDdlEW+ar2NjQ9GbeuSbPGa03j8id3HpUmyrKLc4FDRmHCSbyEMWVTKF+f8mi7u/fPbdeGyr2eLk1XWs7FDIdiCczdC7+ejbSV4bIO4XVThZR0ZXfdzKY9qRBsoFm24FUXrUGPCaYc3fLbZZTfdqo61S+6n7rnN5wgAjgGQ1QsYwJTKp+uIR0eLCrjuI0KY/ivb1tnHuW35jxbXpLU+f/xzxzBtKvZR9Qdj0Ro68rXXXcnUa0TL1eFpFT3Qkt2lyoc/aHNh5HvKhjH8v2IotNWK8Ma9Jgg8s6iAi9c4KfY8ZQanYYZQcFje7z0yW9W6Imcrkz+38omAh16X5JhDDd10mEttXd3YMOhTVyjKDImjVg1OWEtHnQmRUqTMiSlkp8YHxnDvtuewvL+FdxiLwCe1VhpdAKUSYVOnjZL7qIRl4DQZQAAHTVJREFUsAY9JjgvgqjTUBik2zJSqhyvU1FYCqID3jLbKT0HKi8xr7uODu1P9SKKGnpEcV5++Mfr3ENVlyX/OQRpfRcUlCawWaa8JjqrRzbL6gq8/Ndc6rS4uiaZiNOx6VnkuwqYHCsFfm+aIXfRKNikaIzoGeiNhb448w6/+bJjhBy2weRoCeVSGRsObfIkHsNAy/updtdxJ6h0dGKkCcZqMg4AVxrXnwxTda9XQRYflsnW+s8hKW/RkRbmJZn910FHryXfOZec9D9X7knbn2jOLjKTi+ahdG0C9x98wFh6GgBAC5u6aA16jFAxNgKDgWukVEYzrJEDDLwf5g2h6ND+RIyYDc/PTUg6EwPPAMm6D6nGxT29GVY31nQujanr73iur9swxgIJY0d0HQDMGWIB3Cwt7uRWnbQBeNhKUeQHnERsIF15tnATooA16LEizuW2E7t+esnjNQMi1Bi/MhaJkQPMmna4jaKM9udus5fNZyuJW4Gh0pkYolCNVE1cztgcg5PvKmCWsUpc3aVyWPpFfDTR9u6OOo/cDVUx1MOvPSI0lm6WlqyVYlBZA0qJrbS7IYczWehCNkktBFiDHiO0l9sBvfha7LpqQET7UTVUNoHfc5KyZlwrCd3uQqXRCUyXpoVhIh0+eBSqkbKJy4nluw1OdnGuLpE3MzVT3+w5IugURelcB0c50Q/39ZStMhxJidr+NFhUmUIWy37nvdzvuu9+H3dy0uGnB9Udmk+wBj1GiIxPvqvgiTtG5sVz9hOkobIKbkP2uTe/iFUPrRFOJtxlvsv7/vELPzKaaHQqPoMUAfnhmbgwp+EtCm8YX0dS6IJLoFsUpboOw4NFlN+aqh9aNuWpBJ0cn5Qex623TqB6q5KCp7l5/4GN+MXla9x9iT7v27teer2yi3NGhWLzFdagxwiR8enbu94Td4x6mcgLWURh5ERYt3/jXDUhB/5lvnPegLioRWQgdSo+w5b5u4/lVFjeuPQm6UrK5Do6oRLTxF9bV96TT1BBdR1O7zrFrd7M3DDHotLpXTszNVMLNZVGJ5DOpD3P4IZvbMJn3/yCJzxk6mD0DPSirUPs/ZffnsLVMyPScS4EWNpijNCtkoyas55bfAM+9+YX645x7FPf8rzA6Vw6siWqw4Hfk3qMu1Lgxfdl4R6ZgVQVKkWpGqmrssi7h+lcGrOM1WnTu2PEgITamvIWh02OlowqIVXXQWQ8y9en8PSSxysdiQKsHmemZrDI9wz6aY/5zgJ3Mpfdd2mHpGqSduld3QvaS7cGPWboVElGzVkXvaizM7PSv1X87hM7jqB44BzYDAOlCb3b7vRwzQF5RaBfvlcWplBNNKqxRqUaKco9nNp5rO74vK5Qzj5Uja55OvTZfLbO6JlWQvqvg5OAHr8yBkqJJZ6jbNbsL3ITJe9VqyjRJFBDNUlrDbpFwyF6sYOA5+Wc2nms7uVlM6ymvuhvROH3RE/sOOIpEGEzDBeeOYsLz5z1dL5Zu7uv0qTZbyc4L5vI+PNK1t3Q8ZqjKj4STTql0Ym6sv/+Axu53XJ0J3T/eGXNroNAVCwVB9zP4MmdR9WVnwSs3Npb0yXiTYyqWD6wsKtEAWvQmw6ZfEbfoHPad7kTWm6IPJvS6ISw36XbGyweOCcchr86VLQ/pxRcppktK1l3oKLjBWlGIYKuJk9YDRHeikLUsSlo3iO0PjvB25XoyhjynQVMjk8Kw0qAokuVAwZcPn4JSwX3LpvPaskBLOQqUcAa9KZBILVEIo+QV1tXHvfsvS/aBhVVWprSm3PFMNu7xUbQ+dytme0YB5En7ffYZHx7IFgzChFM8htRe4dhWwLqXjcd+Ht1rpMcJ8xqSHTvdK6/pS1qGnQiuhfAXgBpAF9jjD3p+/7fA/gDANMA3gDwKcaYTTlD/2EP5D25EmaO8JboRWrrCq5id3zbkUq8VSX3Wg2r6BpBt2a2CDxvW9RYuNaZKGKKJuANh0xdf4d7LcN4h+7nJN9ZAAPD5LUS8p0FZPIZTF4rGRlLk+uW7yoguzjnUZF0J89VhlKVr8h3KWLfVbQv6zC+R7r6NQsFSoNORGkA+wD0AbgK4BwRDTHGXnFtVgSwmjE2QUQPAXgKwL+OY8CtBJOlf1jvjpVncXLnUWEM+fYH70Dxa+cDqdhNT5SRWZT1SMCKMH5lrM4IypJvqvMWlZ37jZPb6Kha0QWR/3WfT76zUKcIGMY79D8nbuNXGp2oyB+4NMN1YHLd+vauVwpxhTGUfXvX1zGs/HCunyjMJPqN5Z57ocNDXwPgMmPsVcbYFIC/APBR9waMsb9ijDlP4Q8ALI12mK0Jk+rMKGJ/jtfIK/O/eLCI3j9YPaeNYVjTMj1RxqqH1iiLYZzzcHPOZZ697LyHB4vil5tByEWX8a+DSCD4f1ManUCKSCpRYALV6ixIRa+sXF/VtclfLxDWYPYM9GL9cx8TPjtuqV5daYmw13y+QifkcguAn7r+vgrgX0m2/zSAb/O+IKJtALYBwLJlyzSH2LoQLv19yUGlFKkhRBOJE94YHiwqPSY/2pd1YN3+jTWa4okdRyriTAIv2f9bEZVR5NU6RlQ4Hl9M1w0Z/3rfbU9F0uyZx7UOCp3VmekKTrhKkVy3OCGjZwobtQgm83xXoSHn0AqItFKUiDYDWA3g/+Z9zxg7wBhbzRhb/a53vSvKQzclhN5nVfzfVAFPheziHAB1DFlUISgCz1B7qkMVXirX6yJg1fY1RkZUNh4/RF5mlM2eo0qC6qzOTFdwUVXLmkKm6KlT5ets51TopnPpumNMjk8u+N6hIuh46K8DuNX199LqZx4Q0YcB7ALw24wxQXuShQWu181JTDkeotvwfGXJl42TmNPvVJrkirwzSlEljGFgiChNWLm1Ils7tOWwx9vVLd7Rqdw0YWSEWWqr4utR/cYEqtVZUNkCQF0tG2W8XCdn5B+XE0riHbNnoBcndx7FjO89YOVZD0U17mYnrQRiiv6VRJQBcAnA3agY8nMA/g1jbNi1TS+AFwHcyxj7e50Dr169mp0/fz7ouFsG2oaK4JEJldIYU0Amz09SOkU+ot+KKhCdMfhDKCu39nK54lHGL7nnKmBkOF3oozI6gPp8gvwmyLh4LJc4jVTU5yWScHaHeUyPKZKScJ6DuO9LM4KIfsgYW839TmXQqzvoB/A0KrTF5xhju4noSwDOM8aGiOgvAdwB4GfVn1xhjG2Q7XOhGHQ/dB56B8ODRXz7M0eMmwa0d1eKP5wyfT/auvKYLk3X6Y5QNlWbJLKLssi0CQy/a7w6cgDOuYg8KZmWu+zYQRHEq5uPnqDJs6gDmfF1nBXTY8q2B/gyAo3KEyQFmUHX4qEzxo4DOO777I9d//5wqBEuIJgUi/QM9OKlrS8aH8NhtYiogpPXSrVGx55qP19netlEMn5lTCgHAMCon6hpLNq/fRAaoqkxjkobppkQdW5AJzRlekzZ+xK1NMJ8gJXPTQDuRNHpXaewcmuvVjIRkOttUFZ8+0z6Xk6OlYz56e3LOoRyAP7PVfRNUSxaZ/xRdWKazxAlKqOWVI5Dq16WSI1TErpVYQ16zBBxwpf3r6hVxp3edYprgKRGiYAUyTnhvL6XvC5BpiJNzksq+p3/c5VXJjIEvdvuVBqIqDoxNRui6P/q7Ec04UXNhNFhsQQ5poix1CgmTzNDK4YeBxZKDF0YH+Yk/fJdBU/VnjS2LEoauuAkSB1OL6WpFusOqrTnVlZ8MrOLux9KEz4/vbv2t+g8/CXnBELp2oQnbKIKp6jitq0Y+44yWSm69k7JvOi6x4nhwaJH3TOMBlEr3t+wCB1DtwgOWcWeH6XRCf3YssIeu/teAt6CjqDGfNVDazyx8d5td3pi6O7P3RA1f5gcn6y91JOjJW6Juyx2PTxYFMoKOM2no1JdDANToxOluJjoGXKumei6x41yae78TBt3uDEfcxthYEMuMcM0nqcTW1Yh31XweHOmwl+ijuyXj1/y/L30ru5aMVPlh5Umv5ePX/KECnhL8eyN9U2VeaESUejBMdY8Y+7WBWl0OCZIjD/KZKXOM5T0NQl6X6IKQ81nWIMeM0RVkjLIYss6yC7OebwWU0Mg0l7xd6E5vu0Iytfnmgyns2lc+d5rXOPlj4NOXuMXTfGOwdufaJJy64LEXeGpgyDGK8pkn+4zlOQ1CXJfbPJbD9agxwyed7pq+xrpS+YWuHL/Vhf+FyOqrL97PyJ9Ex2vWzYm1TGc/QlDCbOsNpk1AwsiiPGKMtnnf4Z0mU9xoq2T3+xZ9DkgfhZO7jxqvXYXrEFPAH7vdN3+jeg/sLGi1ueD/8V1/1Zb54XB83CLDATv+CL4x2Xi0TlNMtxYu7uvjnbp77Yka2ShY6xVqotJGIIgk4qu5oku3M/Q/QcfaCgzZHiwiHd+wW8lRxKvRfS8TY6WrNfugk2KNghOMsckYWbUPYeTAOT1aXQ37/VD1jzAtAMOL+mVIoJbIsxNwxweLEobWegUaMnOO4lk6fBgEVPX62WNdMXF4kj26Wq8xIFa3kMQ0itdEzfBSKoVYKvD0habHH6D7+7n6Plb8LDLyqCVMro+fRn/b3nMlVnGhBOEeyyqEnAZ3dNhZOhOhv7thB2HIpQU4HX+Aeqpqa2GMDRBlcSDSP6iJqWrQdUFIH1u5wMsbbGBCPMC8Gh3Fw8W65bf6yDmYzva67zjqmR0VWEBZx9+71fYJNq1bFbFlmV0T7dyn4zSyDMEMoMSJjHov1cipUx/wrqVEJYGKru+/lXL8GARJ3ce9V5H1/Pd1lXlz3O0hhZypag16DEi7AugYki4jWlbp7hnqHPcq2dGPN69zLiFCQvodKtX6X7IGjSoULd60FyEhjEEutTQONgkSRXXhOXHC2WdXcwkQK9h+nRpWqgEupArRW1SNEac2nksFA9a1vHIT+F6Z4yfaHIf98KzZz2/EeWg/C+YKURUuanr7ygTtc7LGIbpEaThdlhDoGuoo/Yek6TzhaWBiu7p/Qcf8DxrOvfP6cAVZfJ4PsB66DFheLAolJ41efm5XnSK6h54repP/yaCpsFhXwrnt+7yboBfESjyLMMk73Sur1t2IAqvVidpZzpp6HjeUVaVqhC20YfuPdVNto+PjNlKUR9sUjQmyBJAusm34cGilIUSFdq7O2JZrkettx32uA6cSQuIju3BCxNQNoW29rZAOikqPRdPjoCHGBKDSTX6GNpyWDtM5pejWAiwSdEGQFU4ooOegd46LzdqOMbVMRBDWw7j9K5T2owRmZFqVKWmrPWfIy4GREtdjJoOqMqfqGLMUYZ2/N2UMvlMbN2UTu86pW3MAeDCs2ex9K5u66VXYT30mCDyEtu68kad4oVqghykc2kva0VB83J7qjwDke8q4AMPfrCWSNWh4rlffpFwlozfHhVUE49p56ikedsyFUmdhHZUXnMSXrkbJs+7g/neocgP66E3AKLCl3v23me0H5MCnvXPfUxbKted+Nx321Ncb680OuFRUxSxaByVyKtnRjysA5Veepzqh6rYqu7qIQrFxiATgixeLVvhuOWNo0CSMXrAvGANWNgdivywBj0mRLUE160Obe/u0JbK9XtYUbwQ0xNlYQ9TxyPneeyNquzTTfCFNWhBJ4Tl/Stw4dmzdQlrt7593dhj8FTDhs24za9HSzVHwz8BmVRDO1jIvHM/rEGPEVFk4P0TQ63/pytR6nQhkiUDoyzjF0Hokc8yPDr7RGU5zcH4lbHEwxq6vV3DGrQgE8LwYBEXDxa9oQcCVm7t5U7aorFHgTDMFv9k5s4FiVZpsoI1XsP0hc4798Ma9BaAf2JwMxwoTTWOuSz26BhVP0R6I1FCVSzU1plPvBGF7gpKWAyTopossAxBJgQuD5tV9OjXCca+vH9FLakd5YRo0tRc6zw48E9wIkfIVPtoIcIa9IQRxQPJ9dJU7eg4HpVORV5Y+IuFeMaBUM+rTyIUo7OCEoUA2AzTmnSCeLg6k4B77HF2Zoq7HsB0W7+onYqVtdBgDXqC4L14Q1sOY2jz4TltCk3Oskk1pMijkjWJcMIzNy/vxMh3XtU8Qy/88VGRcRjaotZ+SRr+2O90qVw3aepMOkE8XNNJIIo4v8xgBw0dmoTyTOLgzdJasBlhDXqCEC2lAS+DZHxkDC9tfRFDv/ciMFthqvRuu9NTQKFdbSphPciaRDgNlp0XpQ4KSqQoQcczDjraL0lCFvv1Q3Ufgni4ppNAmDh/nMZRN8GpmuD8E075+lRDVnStAGvQE4SJx+lOMLIZVqMPOkZd6f0QsGq7uIpO1mAZDHh6yeMVfWqR0VaEeJxzFXl/fg+Ysqm6RG+jkl0mqx+dScfUwzWdBMIkLuOkJfIS+iqWix+8CUcES1+0Bj1RhGWTFA+cqxlopffDgAvPnMWFZ856qiPrJEkFCFud2taZF3p/fr56aXQC6Vwa2a58bBWIJtA1DHFOOiaTQJjEZdzVvGGZXiaTa76zEPg48wXWoCeIIBxbN9zetN/7EXrbqIZwPvlNYJbpiXhFAAIJvT93sZKDmakZLFp8g1EVbVyQMVtuuLmtKSYdN8IkLsMKbgWFLjnAZGKZHJ/UYh7NZ1iDniA8L14AT93f4Nft/Yg43g7iFvjyo3RtQtpSjIdGL5lV3XHYLMPktZJxaXoSCOoJh/Hug8Ikbi9c1fLuT3l2wcfRrR56wugZ6BXqhWcX5wACsovqvwOA3m13CvfbbNVy7cs6jMfUyHPw6IoDytxBWN3xpJpUq9AzEG1Dah2ohMfcEGmoi+5Po52CRsN66AljeLCIl7a+yA195LsK+KO3/gQAcGLHkVopPY/l4keN/tcE3iNlU1xFQxl0vcK4CkuCNMUImjxsNtpd0priJnF7UTip2ZhRzQJr0BNEreu5KNbteqDX7d9opPPcM9CLq2dGlBWjiWCWVWR/r03Mya1KErFOtavjoZkwHqIyhEE9uyC/S1rwqhGQTbymcXvRhJN0qKgVYA16glB5gaIHWtcrXbd/I5be1V2nuChCvquA8mQZ028HS9KKwGZYjSVTGp1AppDFqofW1PV/pGwKKaKaHK/KQAc1hDrXzyRW6/+dKRqlEx8X/Nd3ef8Kz73239co4vZR68/PF9gYeoIw6XoOVF6Up5c8jqHNhz09I1/65Dfx9JLHufHXnoFePPzaI3iUPYHPT++Wjuezb34h1rJ/B9MTZfz4hR/VxWrb2tvqtNVlPVeDGELdnpuiWO2q7WvmGlP7erAG9QhFk0Arhgt41/fCs2elMfKwcXsn/+BUGG84tAkPv/bIgjfmgPXQE4Vu13NArrPCyrM1D1jm1coSbQ5jJiqlRRVKoxO4embEUz0qVF8cGePSz4JQ7HS9epnHt666TVTx+yiZJY0Wq5JVP/sh0qIxQbPlH5oN1qAnCNGLzPNOTJJ00xNlvPR7L2Joc8VjaevK4/YH76hIsArghGLCcuNN4G8XJptMeC9pEENomoCTGYWokodRhQsabdyGB4tGzkAUK5CFkH8IA62QCxHdS0R/R0SXiejznO9vIKL/Vv3+r4notqgHOh9gstQ0jaey2Tm3aHK0hAvP1C973XDCCLwxrXrIFWYwhY8r7x0kPOEUEX0T4IdegizVmzW8UQuNzT4ROFxgQv+LGlKdHx5Iv5eu7JiiCaRV8w9RQ+mhE1EawD4AfQCuAjhHREOMsVdcm30awC8YY8uJ6BMA/hTAv45jwK0OXS8vzlCI36vljWkdxH03ZWhfehOmrr8jZLX4l90AaisL2bayscrQiMKZpNDI5KoxzZOFWzWoJpBGT9DNAh0PfQ2Ay4yxVxljUwD+AsBHfdt8FMDB6r9fBHA3EUlcNQsVRN5rZlEW6Vw68H558XrZGEwxfmWs0jdVcPf9L17PQK9wNRDFS9qIwpmk0MjVh+mkEXjFV4VsApkvE3QUIMbkpGUiegDAvYyxP6j+vQXAv2KM/aFrmx9Vt7la/fsfqtu86dvXNgDbqn/+SwB/F3DcSwC8qdyqeRBovIuQ77wRN96SRio3g9mpt/DW62+jdM39OYPQdtaBgc3+EuMjb6N0TXesv4p3/3oKKe1cywxmp/4JP3/5Zty0rID8u3SOvwj5zpvQ3k2glGpb1XibHJGONeR104FwvO/Bu+9II5Xzfz6L2WkCpaIe0y14z2+IvhvDL3/yNkop0VibFGGehW7G2Lt4XySaFGWMHQBwIOx+iOg8Y2x1BENKBK003lYaK9Ba422lsQKtNd5WGisQ33h1Qi6vA7jV9ffS6mfcbYgoA+AmAKNRDNDCwsLCQg86Bv0cgPcT0XuJKAfgEwCGfNsMAdha/fcDAL7LVLEcCwsLC4tIoQy5MMamiegPAZwAkAbwHGNsmIi+BOA8Y2wIwH8BcIiILgO4horRjxOhwzYJo5XG20pjBVprvK00VqC1xttKYwViGq8yKWphYWFh0RqwWi4WFhYW8wTWoFtYWFjME7SEQSeiTiI6RUR/X/3/zZxtPkRE3yeiYSK6SESJV6q2kkSCxlj/PRG9Ur2W3yGi7kaMszoW6Vhd232ciBgRNZS+pjNeInqwen2HiejPkx6jaxyq52AZEf0VERWrz0J/I8ZZHctzRPTzat0L73siov9UPZeLRLQq6TH6xqMa70B1nC8T0f8gol8PfVDGWNP/B+ApAJ+v/vvzAP6Us80KAO+v/vt/AvAzAB0JjjEN4B8AvA9ADsDfArjdt80OAM9W//0JAP+tQddTZ6y/A6BQ/fdDzTzW6nY3AvgegB8AWN2IsRpc2/cDKAK4ufr3u5t4rAcAPFT99+0AXmvgtf1fAawC8CPB9/0Avo1Krd1vAvjrRo1Vc7z/s+sZ+EgU420JDx1eaYGDAOpa+TDGLjHG/r76738E8HMA3GqqmNBKEgnKsTLG/oox5nR5/gEq9QeNgM51BYAvo6IhNJnk4DjQGe+/BbCPMfYLAGCM/TzhMTrQGSsD0F79900A/jHB8XkHwtj3UGHRifBRAN9gFfwAQAcR/Woyo6uHaryMsf/hPAOI6B1rFYP+K4yxn1X//U8AfkW2MRGtQcXj+Ie4B+bCLQB+6vr7avUz7jaMsWkAvwTQlcjoBOOogjdWNz6NiufTCCjHWl1a38oYO5bkwATQubYrAKwgojNE9AMiujex0XmhM9Y/AbCZiK4COA7g3yUztEAwfa6bCZG8Y02jh05EfwngPZyvdrn/YIwxIhJyLasz8iEAWxljs9GOcuGBiDYDWA3gtxs9Fh6IKAXgPwL4/QYPxQQZVMIua1Hxyr5HRHcwxppRA/Z3AfwZY+w/ENFvoVJv8kH7bkUHIvodVAz6/xJ2X01j0BljHxZ9R0T/TES/yhj7WdVgc5eoRNQO4BiAXdUlV5IwkUi42mCJBJ2xgog+jMqE+tuMsXcSGpsfqrHeCOCDAE5Xo1fvATBERBsYY+cTG+UcdK7tVVTipWUAPyGiS6gY+HPJDLEGnbF+GsC9AMAY+z4RtaEiLNWoMJEMWs91M4GIVgL4GoCPMMZC24JWCbm4pQW2Avjv/g2qsgT/LyoxtBcTHJuDVpJIUI6ViHoBfBXAhgbGeAHFWBljv2SMLWGM3cYYuw2VWGSjjDmg9xwcQcU7BxEtQSUE82qSg6xCZ6xXANwNAET0AQBtAN5IdJT6GALwe1W2y28C+KUrVNt0IKJlAL4FYAtj7FIkO21kFtggW9wF4DsA/h7AXwLorH6+GsDXqv/eDKAM4G9c/30o4XH2A7iESux+V/WzL6FiYIDKy3AYwGUAZwG8r4HXVDXWvwTwz65rOdSsY/VtexoNZLloXltCJUz0CoCXAXyiicd6O4AzqDBg/gbAPQ0c639Fhb1WRmWV82kA2wFsd13XfdVzebkJngPVeL8G4Beud+x82GPa0n8LCwuLeYJWCblYWFhYWChgDbqFhYXFPIE16BYWFhbzBNagW1hYWMwTWINuYWFhMU9gDbqFhYXFPIE16BYWFhbzBP8/OmuQzuF4j7QAAAAASUVORK5CYII=\n",
                "text/plain": "<Figure size 432x288 with 1 Axes>"
              }
            ],
            "_view_module": "@jupyter-widgets/output",
            "_model_module_version": "1.0.0",
            "_view_count": null,
            "_view_module_version": "1.0.0",
            "layout": "IPY_MODEL_242932679844407388684f5af4adcdda",
            "_model_module": "@jupyter-widgets/output"
          }
        },
        "197250af0ee648579e97f82dd8006c14": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "SliderStyleModel",
          "state": {
            "_view_name": "StyleView",
            "handle_color": null,
            "_model_name": "SliderStyleModel",
            "description_width": "",
            "_view_module": "@jupyter-widgets/base",
            "_model_module_version": "1.5.0",
            "_view_count": null,
            "_view_module_version": "1.2.0",
            "_model_module": "@jupyter-widgets/controls"
          }
        },
        "0071e8e86e8949f088f9d82f1e5fa3e6": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        },
        "242932679844407388684f5af4adcdda": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "state": {
            "_view_name": "LayoutView",
            "grid_template_rows": null,
            "right": null,
            "justify_content": null,
            "_view_module": "@jupyter-widgets/base",
            "overflow": null,
            "_model_module_version": "1.2.0",
            "_view_count": null,
            "flex_flow": null,
            "width": null,
            "min_width": null,
            "border": null,
            "align_items": null,
            "bottom": null,
            "_model_module": "@jupyter-widgets/base",
            "top": null,
            "grid_column": null,
            "overflow_y": null,
            "overflow_x": null,
            "grid_auto_flow": null,
            "grid_area": null,
            "grid_template_columns": null,
            "flex": null,
            "_model_name": "LayoutModel",
            "justify_items": null,
            "grid_row": null,
            "max_height": null,
            "align_content": null,
            "visibility": null,
            "align_self": null,
            "height": null,
            "min_height": null,
            "padding": null,
            "grid_auto_rows": null,
            "grid_gap": null,
            "max_width": null,
            "order": null,
            "_view_module_version": "1.2.0",
            "grid_template_areas": null,
            "object_position": null,
            "object_fit": null,
            "grid_auto_columns": null,
            "margin": null,
            "display": null,
            "left": null
          }
        }
      }
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C9UDWHVz7Cnf"
      },
      "source": [
        "# Particle system simulation and fitting\n",
        "\n",
        "This tours shows how to use [PyTorch](https://pytorch.org/) and [Keops](https://www.kernel-operations.io/keops/index.html) to compute the evolution of a particles system with a simple interaction energy. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aJNoYic02Cg_"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "import torch.optim as optim\n",
        "from time import time\n",
        "import progressbar"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u1ofi5ng2nPB"
      },
      "source": [
        "Check if cuda is available (be sure to go to the parameter of the notebook to activate GPU)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HLmCgwm-2kkL",
        "outputId": "4a6dfd9a-177d-4436-d542-0466bbe466f5"
      },
      "source": [
        "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
        "print(device)"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "cuda\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Oh1WY5ttkR9g"
      },
      "source": [
        "Draw random particles."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CGKSujUO2w8X",
        "outputId": "d33c7c97-742d-42b7-bcb3-6cc9f35c7f1c"
      },
      "source": [
        "n = 10000 # number of points\n",
        "m = 10100 # number of points on the second cloud\n",
        "d = 2 # dimension\n",
        "X = torch.rand(n,d)\n",
        "Y = torch.rand(m,d)\n",
        "print( X.is_cuda )\n",
        "X = X.to(device); # put it on gpu\n",
        "Y = Y.to(device); # put it on gpu\n",
        "print( X.is_cuda )"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "False\n",
            "True\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e5f6IgJp04gn"
      },
      "source": [
        "Handling boundary conditions."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        },
        "id": "gBG00zg003qq",
        "outputId": "baf4071f-5cbf-4be1-a14e-a89bf565bd4f"
      },
      "source": [
        "boundary = 'no' # no boundary condition\n",
        "boundary = 'per' # periodic \n",
        "if boundary=='no': # change this for usual BC\n",
        "  print('No boundary.')\n",
        "  def bc_pos(X): return X\n",
        "  def bc_diff(D): return D\n",
        "else:\n",
        "  print('Periodic boundary.')\n",
        "  def bc_pos(X): return torch.remainder(X,1.0)\n",
        "  def bc_diff(D): return torch.remainder(D-.5,1.0)-.5\n",
        "\n",
        "t = torch.tensor(np.linspace(-1.5,1.5,1000))\n",
        "plt.plot( t, bc_pos( t ) )\n",
        "plt.plot( t, bc_diff( t ), '--' )"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Periodic boundary.\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f863a6a8890>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5gUVdaH30MGSZKzRMkoMCQxKwiooCKKZEWCiK6L68qq6yquruu67q4roqAgShIQERREUFzDJwxDzjkNYQgDzBBmmHC/P6obGpiZrp6u6qquvu/zzDMdau49NVV16tS55/6uKKXQaDQajfcp4LQBGo1Go4kM2uFrNBpNjKAdvkaj0cQI2uFrNBpNjKAdvkaj0cQIhZw2IDcqVKigateu7bQZGo1GE1WsXLnymFKqYk7fudbh165dm4SEBKfN0Gg0mqhCRPbm9p1O6Wg0Gk2MoB2+RqPRxAja4Ws0Gk2MoB2+RqPRxAja4Ws0Gk2MYInDF5GJInJERDbk8r2IyLsiskNE1olIKyv61Wg0Go15rIrwPwG65PF9V6CB72coMM6ifjUajUZjEkscvlLqJyA5j016AJ8qg2VAWRGpakXfXmZ/8lm+XnfQaTM0ASilmL0ykSOpaU6boglgW1Iq329OctoM1xOpHH51YH/A+0TfZ5cgIkNFJEFEEo4ePRoh09xJSloGAybG88yMNU6boglgws+7+MOstcxfe8hpUzQ+jqSk0fej5bz4ZY4ZZU0Arhq0VUqNV0rFKaXiKlbMcWZwTKCU4o+z1rH72Bkys/UCNW5h+a7j/P3brQBkZWc7bI0GIDMrm5HTV3M0NV1fKyaIlMM/ANQMeF/D95kmBz7+ZTffbjxMtTLFnDZF4+NIahojp6+mcqmiTpuiCeDt77YRvztZXysmiZTDnwcM8FXrtAdOKaX0M3EOJOxJ5s2FW+jcpDIPxtUM/gca28nMyubp6atJTcvgv31aAqBXBnWexZuS+OB/O+nTrha3NqoE6IMSDKvKMqcDvwENRSRRRAaLyHARGe7bZAGwC9gBTABGWNGv1zh2Op0np62i+tXF+Uev6xCnDdIA8M/F21i2K5m/3tecRlVKO22OBth7/AyjZq6hefUyvHxPE6fNiRosUctUSj0S5HsFPGlFX14lK1vx9PTVnDybwZcj2lKmeGHE5/GVUoho9+8EizclMe7HnTzStiYPtq7B2fOZgI4lnSQtI4snpqyigAjv921FscIFEfRTlxlcNWgby7yzeCv/t/M4r93XjCbVjChSfDG+PpGdYd/xs4yauYZm1Uvzl3ubAhePicY5/vLVRjYdSuFfD19HzXIlABDRN2EzaIfvAr7fnMTYpTt5OK4mD+m8vStIy8jiiakrEWBc39YUK1zwku/1TdgZZibs5/OE/Yy8rT63N6rstDlRh2sXQIkV9ief5fefr6FJ1dK82qPpJd9dSOk4YFes88q8jWw8mMLHA+MuRJEQeEz0UYk0Gw+e4s9zN9Cxfnl+3+naS74TBKXvwkHREb6D+KNIBXzQ78ooUuMMsxL2M2PFfkbcWo87Guso0g2cOpfBiKmruLpEEf7TuyUFC+jUWn7QEb6DjPl6ExsOpDBhQBy1ype44nv/KW1ELvoEjwSbDqbw0twNdKhbnlGXRZGB6GAyciileG7WWg6cOMfnw9pToeSVcyF0Dt8cOsJ3iDmrEpm2fB/Db6lHpyY5R5E6pRNZUtIyGDF1JWWKF+bdR1pSqOCVl4culoo8E37exXebkvhTt8a0vqZcjtvoKh1z6AjfAbYcTuGFL9fTrk45/tA59yhSEzn8chb7T5xjxtD2VNQzal2BX86iW/MqPNaxttPmRD06wo8wqWkZPDFlFaWKFea/fXKOIv34a+915GI/H/1syFmM7tKINrVzjiIhsFRWHxS78ctZXFOuBH/v2SLPuSgietDWDDrCjyBKKf44ex37ks8y7fF2VCql9T/cQPzuZN78dgtdmlbh8ZvqOG2OBkPO4qlphpzFZ4PbUqpYYadN8gQ6wo8gH/+ym4UbDvPHuxrSrm5503+nSwDt40hqGiOnraLm1cV5q1feUSQEjKvoQ2Irb3+3jeW7k3njfvNyFvqQBEdH+BEiUBRt6M11nTZHw0VRtJS0DCY/1pbSJqJIPWZrP4GiaA+0qmHqb/Rgujm0w48AV4iimTw7dTRpL35RtLd7XUfjqqGJoulDYg9hiaLpgxIU7fBtJidRNLNo3Rb7uFwUzSx6IN0+chJFM4sg2t+bQDt8m/nX4m38387jvPVgiwuiaBpn2Xf8LM9eJoqmcR6/KNrEQZfKWWisQw/a2sgPW5J4b+mOfIui6ZSO9aRlZDFi2kogZ1G0YFyY/azjSUvxi6I9eVu9fImiiehSWTPoCN8mDFG0tTmKommc49X5G9lwIIWPBuQvitSDg9az6WAKf567gRvqlWdUp4b5akMfFnPoCN8G0jKyGDF1FdlKMa5faLnIQHQ0aS2zVyYyPX4/T9xajztzkbMwiw4mrcEvZ1G2hCFnEY4omj4kwbFqicMuIrJVRHaIyOgcvq8lIktFZLWIrBORblb061bGfL2J9QdO8c9e13FN+avy3Y5O6VjH5kMpvPjlejrULc+zeYiiBUOvPGYdSin+MHMtiSfOMbZPqxxF0cxipHQsNM6jhO3wRaQgMBboCjQBHhGRy+upXgJmKqVaAr2B98Pt1634RdGG3VKXzk2rOG2OBiOKfGJK3qJooaJ9S/j4RdFGd21EXB5yFhrrsCLCbwvsUErtUkqdB2YAPS7bRgH+EpUywEEL+nUdgaJoz3XOXy4ykAu6LWG3FLsEiqK916eVdaJoOpwMi0BRtME3hi9nISI69WkCKxx+dWB/wPtE32eBvAL0E5FEYAHwVE4NichQEUkQkYSjR49aYFrkCEUUTRM5AkXR2taxJorUWZ3wCEUUzSz6kJgjUl7pEeATpVQNoBvwmYhc0bdSarxSKk4pFVexYsUImRY+Sime/8IQRXvvkZaWiaJdzOHryCU/2CmKpo9I/ggURXu/XytLRdH0ZRIcKxz+ASCwyLyG77NABgMzAZRSvwHFgAoW9O0KJv66hwXrQxdFM4s+j0PnaGp6SKJooaCjyfyTH1E0U+gVr0xhhcNfATQQkToiUgRjUHbeZdvsA+4AEJHGGA4/unI2uZCwJ5m/LdisRdFcRKAo2rh+rU2JooWKjiZDxy+K9khb86JoGmsJ2+ErpTKBkcAiYDNGNc5GERkjIt19mz0LDBGRtcB0YJDyQJ4iv6JoZtG6LfnjncXb+G3Xcf56X/OQRdHMoAcIQ2ff8bOMuiBnEaIomglEh/imsGSmrVJqAcZgbOBnLwe83gR0tKIvt5CVrfjdDEMUbc6INiGJomnsY8mmJN7/cSe924QmihYKOqUTGmkZWTwxdSUFRPIlZ2EGPZBuDi2tkE/+tXgbv+44zls9W9C0Whlb+rhwDuvIxRT+KLJptdK80t1eOQv91GWeV+ZtZONB+0XR9FNXcHTtYD7wi6I9FFeDh9qELooWKvpEDk64omihoKNJ88xK2M+MFfkXRTOLoG/CZtARfogEiqKN6dHM1r60YzFPoCharfL2S+tq3xKcTQdTeClMUTSNtegIPwTSM60RRTPLBfE07V3yxEpRNDMIoo9JEKwURTOD6DFbU+gIPwTGzDdE0cb3bx2WKJrGOjYfSuGlueGLommsI1AUbcbQ9mGJommsRUf4JvlydSJTIyyKdqEsMyK9RR9GFLmK0sWsE0UzhehxlbxwQhTNeOrSxyQYOsI3wdbDqfxpjnWiaKGiT+Qr8Yui7Us+y/Qh7a0TRTOBHlrJHb8oWtdm1oiimUWndMyhI/wgpPqkdZ0QRdODtrnz8S/Wi6KFhPYuV+AXRatVrgRvPWitnIXGGnSEnwd+UbS9yWeZ9ng7y0TRzHJxxStNICv2JPO3hVu4q2lly0XRzKD92JX45SxS0zL4bHBbS0XRzKDLMs2hI/w88IuiPWeTKJomdI6mpvPkVEMUzQ45C7No33Ip/1y8jWW7knn9PotF0TSWoiP8XPCLonVqUplhTomiaS2dS/BHkafOZfDJo21tEUUzgx4gvJTFm5IY96MhitbTJjmLoOjHLlPoCD8Hjp1OZ+S01VS/ujhvOxhF+tEVIQYXRdGa0aSac1Gk9i0XsVsUzSwX56zoayUvdIR/GX5RtBNnzzNnxA2OiqJpv3KRQFG0XnH2y1kEQ/uVyIiiaaxFR/iX8e8lhijaaz2a2SaKZhbRo7aAIWcRKVE0M+gbsYFfFO2dh66zVRTNDBdXh3PUDNejHX4AS7cc4b8/RE4UTRMcfxQJ9ouihUKs+xW/KNqIW+txR2P75Sw01qBTOj72J5/lmc/XREQUzSyCnmn76vxNERVFM4NIbGvp+EXROtQtzyiXyFnoa8UcOsLHEEV7clrkRNFCJVadyxcrE5kevy9iomhmieWUjl8UrUzxCMtZBOFiSidGLxaT6AgfQxRtXaL7RNFiuRpky+EUXpy7nvZ1y7lSFC0WK6f8chb7faJokZSz0FiDJbdnEekiIltFZIeIjM5lm4dEZJOIbBSRaVb0awVOiKKFSqw5l5S0DJ6Y4oAomlli9Eb80c+GnMWfujaiTYRE0cyi6xvMEXaELyIFgbFAJyARWCEi83zr2Pq3aQD8CeiolDohIpXC7dcK/KJobR0SRQtGLPqVy0XRIi1nYZZYyxzE707mzW+3RFwUTWMtVoRObYEdSqldSqnzwAygx2XbDAHGKqVOACiljljQb1icTs/kiSkrKVm0MO+5MYokNkvN/KJoz3dp6Iwomgli7UZ8JDWNkdNWuVoULRavlfxghZerDuwPeJ/o+yyQa4FrReRXEVkmIl1yakhEhopIgogkHD161ALTckYpxfOzDVG09/q0pFJpd0aRfmLlHA4URRtyk0NyFiZwo8OzC7+cRUpaBuP6tYq4KJpZLq4dEStXS/6IVFhbCGgA3Ao8AkwQkbKXb6SUGq+UilNKxVWsWNE2Yyb9uodv1h/iubsa0t7FomgSQ7GkW0TRzBIr1SBaFM1bWOHwDwCBs5Rq+D4LJBGYp5TKUErtBrZh3AAizsq9ybzhtChaiHjdufjlLE6dy+D9vq0dE0Uzi8vvRZZxURStpnOiaCHi8UslbKxw+CuABiJSR0SKAL2BeZdtMxcjukdEKmCkeHZZ0HdIHDudzpNT3SOKFhSXm2cV7yzeyv/tdF4ULRS87lcuFUVzXs5CYw1hO3ylVCYwElgEbAZmKqU2isgYEenu22wRcFxENgFLgeeUUsfD7TsUAkXR3u/bylFRNLNcVAB01Axb+X5zEmOXukcUzQxevw/75SwEd8lZ5IXbYze3YMnEK6XUAmDBZZ+9HPBaAaN8P47gF0V7q2cLx0XRNAb7k8/ye5+chRtE0czidWmFV+cbomgfD4xzXBTNLLE03hUO7qtFtIFoFUVzfcopDPxRpAI+6BcdUWQgXq0Gmb0ykenx0SuK5uUbsRV4Xloh8YQhitbYRaJooeLFk9gvijbBRaJoZvHqbXjzoRRe/HK9q0TRzHKhDt+jN2Kr8HSEn56ZxYipq8jOVozr6z5RtGB41bH4RdGG31KPTi4SRQsFr92EDTkL94miaazF0xH+a18bomgf9m9N7QruEUUzixejlkBRtD90jq4o0o/XMm1eEEWLhQIHK/DsbXzu6gNMWbaPYTfX5S6XiqLFGqluF0UzjXjoFnxRFG10F/eJopnFazdhu4jWKy5PtiUFiKLd5T5RNLN4SR9EKcUfZxuiaO/1aeVaUTSzeOGYwEVRtC5Nq/D4TdEviuaFw2LnREvPOfzT6ZkMn7KSq4oWcq0oWqh44ST++JfdLNzgblE0s3glmjx6/BjvT/mcmlcX561e7hRFM8uFFa88cCf+dPYcXvtiGdnZ1u+L53L4p9MyKVO8MG/c39z1omjB8EptccKeZN6MAlG00Ihux5KZcZ7y/63PJyg299ntejmLWGHN4ml0WP9Xzl77ZwoUaG95+55z+FXKFGPOEzdEdbRyOdEctRxPSkRN7s3k4sVo2fwJROKcNilsvHBm/TTlr9zuu2k1rlzSYWvC52KBQ/RyYNdG6v76LEmFq/Foz+7B/yAfRH++Iwe84uyjfTeyMjM5NLEvzbO30DFrBSUOLnfaJEsQie4c/pJNSYzceh3HivgF0aJ4ZzxC2tnTpE3tBwhXd3mBYpmnbenHkw7fa0Tr5Rg/6Vmapa9h3fV/gWJlid498Q4H9u7kxZnLqFu9EqXvGQPXdsUbzywG0XojXjdhKPWydrHnxn9Q4ZvHYNWntvSjHX4UEI0n8ZrvZ9DhwCfEX303be9/KvofVwIQolNLJ+3cGc592otxvMG4Pq0o0qIn9JkBhYo4bVrYSBTndL78bSPlj69iWfVBtLjtYd+n9uyI53L4XiJaU1P7k8+S+vMH7CxYlxZDxl/8Ihq9ZC5E42S4tROG0y5rJ2tu/ICa5aNvIqIX2XjwFKO/2ceNtd5n/KM32x4Y6Qg/Koge55KWYchZPM1zFB04h2IlfAOCBQqDeON0i8b7cPzcsbRLnsdv1QZw/Z2PGB8uHw9vXwsZ55w1zgIuzLSNomvl1MnjrJ70e6oUV/y9T0cKFgqIv20KjnSE72Ki0K/w7eQ32HfgWt4ecCs1agXo5Dy33TmjLEaIroeVXRuW03z1K2wseh1tHv3nxS8yzsDppOjamVyItpuwys5m1/j+PHx+Ga3v6UeFkn45Cx3hxyzRNtN2xVfvc1/i27xbf3XUiqJ5jZS0DF74ZjcJBZpT+bGpFCocmK+PMi9pgmi5VpZPe5WWZ39l5bXP0LjN7Re/EIH7P4RGd9vSr3b4UUA0nMO7Ny6n2aq/sKlIczoOHHPlBgufhzXTI2+YDYhEh5aOys7mjzPXsuJUaYoOmE2FKrmtBRENe5M3F1M67mfTbwuJ2/4uq0reTLtHXrr0SxG4rjdUsUfK3RKHLyJdRGSriOwQkdF5bNdTRJR4YfZNBIiWmbapp5IpNHsQZ6QElR6bdlkU6WPDF5AYH3njbCIaIsnl01/j/u3P82Kn2jnLWURbHsQDHDl1hhKLRnGoQBUaDJmMFMjBBe/9Pzi5z5b+w3b4IlIQGAt0BZoAj4hIkxy2KwX8DvDG7JsI4mbnopRi9cdPUzX7MEfu+oAKVWrltXHkDItxNi9fRNy2f1OxdHEevaVRzhtVaAjNe4FE1zoROeGvaHPzrPTMrGyemrGOYVnPkdHrU0qVyUVTalJXWD3VFhusiPDbAjuUUruUUueBGUCPHLZ7Dfg7kGZBnzFBNARgH/+ym2ePdGVpszdo0qFrHltGwc6EgJurQY4d3k+FhcM4XKAy9XOLIgEadoGeH0Hh6NacihamfzGb5buPM+z+ztRt0ibI1vacX1Y4/OrA/oD3ib7PLiAirYCaSqlv8mpIRIaKSIKIJBw9etQC07yBW53LuvVr+MfCjbRs0pA7Hxxu4i/cuR+h4uYbcWbGeZIm9qWkOsP5npMpXba80yZFBLfPu1qzeBr9Nw3hnfrreKBVjeB/YNOTiu2DtiJSAHgHeDbYtkqp8UqpOKVUXMWKFe02zfW4eRWf40mJVPnifv5VYhL/6HVd8ElixctC4ehauzY3RHCtZ5n09Y9USt/N+uv/Qt1m7fLeePmH8Ho1SEuJjHE24uZr5cCuTdT99Vl2FKxHt76/M/EX9p1gVtThHwACh/9r+D7zUwpoBvzocwpVgHki0l0plWBB/57FrZGkXxStvkqlYffnKFPchLTuyBX2GxZBXOhXWLIpideXZ3Cg1Qxeub9j8D/IyjBq8V25N94g7dwZzk3tR0mEEv2mUay4iRnONl74VkT4K4AGIlJHRIoAvYF5/i+VUqeUUhWUUrWVUrWBZYB29iHgtqglftIfaJa+hvXX/Zl6LTo4bU7EcWP11IFdW9g882WaV72K0feb1FF3a0SRH/yDti67ea2dMJz6WTvZc/M7VKuTy+D55Tz0GTR70BZ7wo7wlVKZIjISWAQUBCYqpTaKyBggQSk1L+8WNLnjvgty7Q8z6HBgEivK3k2bB8w8nvqY/zuo0gLaDLbPuAjipmoQI4rsw0AO80CPP1CscIhVNy7aFy8xK2E/8w83Jq1+TW65vbf5P2x8j202WSKtoJRaACy47LOXc9n2Viv6jCXcErXsTz7Lmz+d4KnC7Wk1dHzwPwhk2yLIzrLHsAjjtsB43YRhtPWJol1fu0EIf+myHQmDC3vijkuFzfuO8NLcDbSuczs3DgoylnI52xfD1bWhQijH0hx6pq2LcZO0Qtr5DEZMXcUG6lJ9eIAommm85VxccEgAWDH3Pdomz79UFM0slZtA60FQ0AvyyMZvNxyXlJPHuWrSrQwu+j3vPtKSggVCPPen94a19sxK1+JpLsZNLnLth4/zQNI5qvf+D9fkW1rXDZejd9i2N5GGq1+/UhTNLHVvNX40lqGys9k5vj/Nsg/TvUvnAFG0UBuK0rJMTfSz4qtxtDs+l0Y1K9K5WdX8NSLiGX8v4vwCKClpGQybtYNnCv2Zyo9NyVnOIhhKXfyJcvwD6U7vyrJpYy6IojVq1zmfrbi7SkdjE25YAGXXpniarXqZjUWa0+axf+W/odLVjFp8j+CkX1HZ2Xzw2TT2JZ9leL/eectZ5EX8eHi1LJw7Ya2BMcqm3xbSZvt/chZFCxn31uFrbMapqCX1VDKFZw3kjJSgcm6iaGZ5fIl1hjmM07fh5dNf448H3qFxu/E5i6KFitNhsQVczOE7sy9HUtOYu2QpfQpUy10UzSwur8PX2ISTq/gopZgway5ls08EF0WLNcS5sky/KNrqq27knnt7hdma07cu63Bypm1mVjZPT1/Np+dvI+3x/+UuimaWvrOhZX9rjLsMHeG7GCczOhN/3cO7Oypz9Z3f8miH68Nv8MvhULkp3PBU+G3FKIGiaPWGfBpeFHkJ0R/hO8mvn7zAVXtK8HrPx2hUvUL4Dda9Jfw2ckE7/Cgg0lHL5vjF7Pj2Ozo3eZhBd1xnTaN7fsUrEaUTZZmZmZkkTexLPXWalJ4zrBFFc8EYkVU4VZa5ZvE0btk/jmJVutOutQlRNDNsnAvl60GV5ta0F4B2+C7GievxeFIi5RcM5cnCRSh130sWDxx7KJKM8K68s2QHR8+048E2D9OuuUnphGBUaQ4dRkIhLY+cHw7s2nxBFO26IR9Y1/CXw6HtEO3wY5VI+ZZAUbQDD8yjRuky1jUebSt/50Gkq6eWbtjH+z/upHeb/rS7r4V1Dddqb/x4gItlmZE5xww5i76UCkUULSR0HX7MEemT+FJRtBssbt1D6QMiN5B+cPcWmsy+hUEVtvBK96bWNp553pBGzs62tl0niPCs9PlT3qV+1k523xSCKJpZdJWOxm6Wxf9Gu8RPiC/bLTRRNLOUrwelqljfrodJO3eGs1P6UJx0hjzQLXRRtGCs/hTerAlnj1nbrseZlbCf53Y257NmE7nujhBE0ULBpjuXTum4mQgNRO1PPsuwhancW+oVXhr6uD2d9P/SnnYdQCQykeTaCU/QLmsnazqO4/q6je3ryAOptkg9P+7cEM/Eueu4oV4z+jwQoiiaaezbG+3wXUwkTuL0tLO88+nnZKtqDHl0CMVKWJ2L9CZ2+8gVc9+nXfJX/FZ1AB069bGpl4szPTTBSTl5nKJfDOCDQsJVD68MXRTNLIO+hqvsWfFPO/wowE7nsmbCE/z9xHyW3f1dGKJoJpg1CCo2gltH29dHhLB7AZQth1NYvWoZJYq1oM1j+RBFM4unyjLt1dJR2dnsnDCA5tlJbO82g2tK27hcZ/VWtjWtHb6LsbsaJGGeIYq2rFpfbm4bZ2tfHFoLBbxxuhk6cPZ4lpS0DJ6YsoozRQdw34i24clZmEWndIKyfNoY2p/5hWXXjqJ9u7vs7WzNdEMLv4b116QetI0KrL8gd29aQdOVL7OpSHPiHvu35e3niAcci52o7GwSxg2lzIn1vNenFZXKlra3w2ot4dY/QRHvpPHsuBFvjP+BuO3/YfVVN9HukT9b3v4VLHgONnxhS9OWOHwR6SIiW0Vkh4hc8cwuIqNEZJOIrBOR70XkGiv69Tp26YOkppyg0KwBnJUSVHp0amSiSFctGxI+dty7lk//K7efmsPzjY9bI4oWjGotjRRbMZtvLBHArsWCjqSmMfS7dKYW6UX9cEXRXEDY1otIQWAs0BVoAjwiIk0u22w1EKeUagHMBt4Kt99YwI6MjlKK0fN2MDvjBpI6j6NC1Qjdez2UL7aDzcu/o7VPFK19nxxXB7We82cg5aBnlp60msyM8zw/9ReOp0O7wW9Tygo5CzPYWAZmxe2qLbBDKbVLKXUemAH0CNxAKbVUKXXW93YZYJHoRGxg5aGf/NNWvtlwhOKdXqDJDd0sbDkIlZtBubqR689GRMTSY3I8KZHyC4eRVKCSxaJoQVg/G95pDKmHI9OfjdihpbNi0rO8cvAJ3upWi0ZVIvkUZN/TsBVnVnVgf8D7RN9nuTEYWGhBv57H6mqQLfGL6fJDFwbXO8WwmyPsfB+aDLeHuyiEO7BSJSIrW/HLZ2MorVJJf+ATa0TRQib6U21WXyurl0ynw8FPSarQju4dLJ7hHAwbH4YjWjYhIv2AOCBH/U8RGQoMBahVS+uv+7HCuRxPSqTcgqGcL1CUpx/s7IrVtDTwzuKtjDt2Lx/e+QCdrBJFM4sHzwErZEgO7NpMvV8MUbQWQz60wKoQefx7KGrPE4UVEf4BoGbA+xq+zy5BRO4EXgS6K6XSc2pIKTVeKRWnlIqrWNGeiQfRxMWBqPBOYr8oWhmVSvr9kyhztQWa3aEyvQ8sejHy/dqAcVzCdywrf/qamUsT6BV3DZ3u7BJ2e/nGA9VTVqV0/KJogE2iaCao0ABKVbalaSsc/gqggYjUEZEiQG9gXuAGItIS+BDD2R+xoM+YwKr4K/6T52iWvoZ1LV6yQRTNJMd3wKlEZ/p2IQd3b6HBD0N5v+REXu0R4ZTBBbwX4YfLP+av5HhGYXbZIYpmlhUfw+6fbGk6bIevlMoERgKLgM3ATKXURhEZIyLdfZv9AwAPydUAAB0QSURBVCgJzBKRNSIyL5fmNDkQTtSydPNhkvduYEXZbrTp+YxlNuWP6I8kIfwiCr8oGiiq9RlrvSiaWaq3hs6vQzELZbAdJpzjMntlIh+vSuV/HT7hertE0czww2uweb4tTVuSw1dKLQAWXPbZywGv77Sin5gjzABsf/JZnpm5jurl/8yc4W2ssSm/REpxLAII4VXprJ0wPDKiaMGo3MT48QDhjknt2rCcq7/6E52v+SOjOje0yKowcHFZpsZm8nPs09POsvXDAVRThxnXvzXFihW33rCQ0OkDgOXzJtAueR6/VR3A9baJopkk7RQc227o4nuG0C+WlJPHKfLFQFoU2MUb9zelUEGn3aLWw49JLiyAko+TeM2EJ7gzfTFv3FDAXlE0s9RsYyxi7gGMh5XQj8mWwyk8taIsX5TsY68ommmDvoH34iD1oNOWhE1+Z6Wr7Gx2TBhIlewkjnX9kApVagb/o4ig9fA1Jlkx7wNDFK1KX9p37ue0OQbd/+u0BY6SmnKCpz9LgGJXc9Owf0VIzsIkHkm15QdDFO3nyIiimcXGclnt8F3MheMewvW4e1OCIYpWtDlxgyMkihZjhDoPUmVns/3DAfwt9TCZA76hUim3LBrunTRbfsoyE3Ycour2qawqeXNkRNHMMuxnKGKP/LJO6biYUC/H0+mZ7JvzEmelOJUGRUoUzSRTesJ8p6uELEIkpKB4+fS/0urMT2Q06Eq7+pXssyvfRH+EH+pM26Op6Yz4fCNPX/U29Yd84i5RtDLVofjVtjTtor3U5IaZy1EpxfOz1zHi7DAO3DudCtVcJkiacgjOHHXaioizxS+KVqIj7fq+4rQ5l+LJmbbBt8nMOM/8j17lbNo53ux/m0NyFnnw21jYtsiWprXDdzGhrOLz3VdTWbp+N0/d1YLrWjs0uSovPFWWae4mfDwpkXILh5FUoCJ1IymKZpYabeDe/0AJlzm8fHAxpRP8yCRMepbHTo3lw/YnaFzVhdLQv/wbttojN+ayM1CTH7bEL+b21U/zbqWvIy+KZhrvRZN5kZWteH3Oco6oMqQ/MNkZOYtglK8HrQd5auJVMNYsmU77g58SX+5eOt49wGlz8kDX4cccZqIWvyhaUoGKtHn0Hy4XRfNIhG+iLPOdxVuZs7cYm+6ZR71Ii6KZ5WwyHFgFGWlOWxI2ZsoyD+zaTF0nRdHMYuM1rB2+iwl22LMyMzk4sZ+zomhmqXuLLWt0upE138+gzi/P0a9VRR5q42LV1x1LYMJtkHKF1qHnSDufyclpjwFQot9UZ0TRQsGm9Kcuy4wCcjv2v336Ijemr2ZFi1dp45Qomlnuet1pCywjrxvxwT1bqfPzKEoVqcw990SJbIEHxlaCLXH46tebWH1mAH/rXJmWdRyUszCFrsOPSfKqLV669QgvbWvMi9cMpZvjomixheRSlpl27gxnPutDSRTF+k6lWImSkTcuJPIx0cO15D4rfeFPvzE9PpkRt95Gy9scUsAMhZHxUKCwLU3rlE4UcuDQQUbNWEXpqvW4/fE3nTbHHJPvhS8ed9oKW1n70RM0yNrBzo5vU71uFET3dq387SJ2b1zObd/fy6uVfmJUp2udNsccxcrYNvFKR/iuxl+WefGCTE87y9mPuzMmuwrN+37unLRuqKSdgkJOC7hZg1GWeamTXPDrSm48tpjfqvWnQ6e+zhgWw+R070o5eZzCsweSKiW5u89IF4iimeSnt40Kqqb3W950lPwHYpOcBuvXTHiCBpnbqdHxYWpXcPnA0yW4uXooPLYcTmHUt0cYXfkD2jz2jtPmmKdmW+j5sW2rKzlJoCja0S4fuEgUzQQJE40BdRvQEX4U4A9aEi4RRevvqE35wxupg8A5ZKmnkvl60tuUKnobr/S/zV1yFsEoW8v48QCXhxPx08bQ7szPLGvwe9q3d3D5yPxi06WiHb6LuXASK9izOYEm0SyK5ur5AaEhGIO2Kjub7eMH8kzaL3Tufo+LRNFMcvooHN1srHxVJJqeFq8kcFb6ij3JTN6syL66M+37vBzkL92IrsOPaU6nZ/LP+SvYL9XcJ4pmlgZ3QZ2bnbbCUpZNf51WZ35iRf2naeFGOYtg7PnZGEw/ud9pSyzjWGoaT05dxaayt9J05DT3yVmYxsUzbUWki4hsFZEdIjI6h++Lisjnvu+Xi0htK/r1Ov6o5c2FW1hw8hpO9F/iPlE0s9z2J7jhKaetsAaBEkkJxG37F6tLdKR93784bVGYRH+qTYCCZFF81sN0TV/IuH6tKV3MntJG23GzHr6IFATGAp2ARGCFiMxTSm0K2GwwcEIpVV9EegN/Bx4Ot+9YYGDBRVQ8fZLCnV+ifT0Xz6SNIQqp84zJfMe9omhm8VCaDWBUoVm0z15NwVY93SmKZpan14DYc05Z0WpbYIdSapdS6jwwA+hx2TY9gMm+17OBO8Tdoi+uoNTRVbxUaAo3lj7K0JvrO21OeEzsCtM8co8vWJQXs4eT9sAn7pazMIsH6vArHV7Kk4XmEV/uHtrcF+VPkgULgU1BhBWDttWBwCRgItAut22UUpkicgooDxwL3EhEhgJDAWrV8kb1QL45c4y6P47kdImq1Bv6WfRGkX6yzhs/0U7ybp7tfC0ZWQ2oXzfaZYU9EnOd2EPz+D+SXLoRLYaMd9qa8Pn+NaMO/3rrF7l3lRdRSo1XSsUppeIqVqzotDnOkZ0FXwxGzhyj1IBplCrrgSjSCw902xbBf1vTOm057aPe2QO12sMjn0OZGk5bEh67fkSkAOUGzXC/KJoZ1s+EXf+zpWkrIvwDQOCshhq+z3LaJlFECgFlgOMW9O1NDq6Bvf8Hd/8Tql7ntDXWEc2pgxN7Yc5QqNwE6t7qtDXWUKoKNIzCGvXLaT0IGneHEuWctsQi3F2WuQJoICJ1RKQI0BuYd9k284CBvtcPAj+oYILisUyN1vBkPLRy8wINoRLFEX5GGswcYNywHvoUCntDIoLUw7BlgSF7EY1smAO7fzZee8bZ+3FpWaZSKhMYCSwCNgMzlVIbRWSMiHT3bfYxUF5EdgCjgCtKNzXAyX2wca7xulwdZ22xmibdoWFXp63IH9+OhkNr4P5xUM6tK4rlg8QVMOMR47yLNpI2wtwR8PPb0f3kmBM2LgdqyUxbpdQCYMFln70c8DoN6GVFX54lM92IIo/vMiYoeS1iidYafKWgYkO46VlodLfT1thDtDnMtFPweX9DVfKBCd4YHwqkUHEoaM8cAi2t4BYWvQAHV8PDU73n7AGyMo3fBaPolFPKcCbtn3DaEpuIQkepFHz1JJzYA4O+hpKVnLbIep5cZlvTrqrSiVnWzYIVH8ENT0Pje5y2xh4+uRumPOC0FeZJS4GP7jAqczxPFEX4W76GzfOh06twTRTKWTiMdvhOk3oY5j8NtW6AO6J9in4wosSxKAXzRhrVUlEuKpYn0ZgKaXQPPPQZdBjptCX28d1LED/Blqa1w3eaUlXg3neh16ToSneEio0DUZazbBxs+grueBlq3+i0NfZRsz0M+iY6BqJPHzFKY0WMAoBovFmZZcsC2GdPWkc7fKdQCpJ3G69b9DIcv6eJkgt03zJY/GdoeDd0/J3T1tjLVeWNG1rRUk5bkjdZmTD7MZjYxSiRjQlcWpapySfLP4T328PhDU5boglk6wIoUxPue9/bUSRAykFYNxPOJjttSd4sfd2Qcr79JSgcZWsO5Acbzzvt8J1gfzx89yLUux0qN3XamsjQohc06+m0FcG581UY8gMUL+u0JfZzaB3MGWJUvLiVrQvhl3eMSYgtY2itYDfX4WtC4MwxmDXI0C+5b5z3o0g/cY85bUHexE8w0huVGnuzLDZPXDq2cmIvfDkMqrSArv9w2prIUaK8bWk27fAjSXYWfPG44fQfXxwbUaSftBTjdzEX6pRv+w4W/AHiBsM9UbQIebi4PdgoUR6a3g8dn4mNVI6fwd/Z1rR2+JFEKSOF0/R+b4mimWFGH1DZ8OiC4NtGkhN7jbRGleZw1+tOW+MMbgzwM89D0ZJw73+ctsRT6Bx+pFDKKLu863VoPTD49l7EbWWZmekwa6D3RNFM49IIf810+PBmY45KLLLgOfjl37Y0rR1+JDi5Dz68yZBO0LiHZeOMY+I1UTSz1GoHQ/8HlRo5bclFkjbC17+HqypACQ+sA5Efdv8EB1fZ0rRO6dhNZjrMHGikDoq6MH8dUVwW4bcfARWuhUbdnLbEGYqVgWrXO23FRQJF0Xp+7O2JiMGw6WlYR/h2s+gF42593/vGsmWxipsGCI/vNGrPCxWJXWcPcCoREiYas1idJlAUrdckKFXZaYscxL5rJYZvoRHggijaU9D4XqetcZaW/Y0qJadJS4FpD0HxckY1hJtuRJHmyBYjfVK5mfOqk2mnjNTnna9oUTTArqdh7fDtZNNcnyjaK05b4jwtHnLagouiaMm7YeC7se3swV1jtsXLwuDFULCI05Y4T5kacJU9N2Dt8O3koU8hPTW2c5F+Th8xHK6Tj+p+UbROY6B2R+fscBtOVk+dPgo/vAadXzNy9xroN9u2psPK4YtIORFZLCLbfb+vzmGb60XkNxHZKCLrROThcPp0PUrB//5h6JQUKBhbk6vyYs5Q+Lyfc/3vjzdE0RrdY6w7oOFiiO+Qw8/Ogi8eg3Wfw8n9ztgQY4Q7aDsa+F4p1QD4npzXqj0LDFBKNQW6AP8WEe96wfjxsPSvsPFLpy1xIQ5GkmWvgRYPQ4+xOpXjx+n/w9LXjRLEu9+BKs2ctcVNzHsafrBnEmC4uYYewK2+15OBH4HnAzdQSm0LeH1QRI4AFYGTYfbtPvbHw6IX4dqu0M6ry+LlE6ecS3bWxVTSfe87Y4NbqdkOnloFpatFvu9ti+DnfxqD+bEkimaGAyuNAMUGwo3wKyulDvleHwbyTNCKSFugCLAzl++HikiCiCQcPXo0TNMijF8UrXQ1YyJPAV3xegVO5IqXvgGf3QcZ5yLft9spcpVRKhzpGcbZWcaqTlWaQ7cYEkUzjYNlmSKyBMhpdY4XA98opZSI5HpFi0hV4DNgoFIqO6dtlFLjgfEAcXFxLpulE4QlrwSIol0xlKFxoiRk2yL4+W1o2S8GZRNMcHIfbJxryFaXqR65fgsUhAFfQVaGPi654lBZplLqzty+E5EkEamqlDrkc+g5zuAQkdLAN8CLSin7lmR3krteh2YPxJ4omlniHoOs9Mj1d2KvMVBcuTl0ezty/UYTybuMgewacZFz+Fu/hQadnEkjRQuCa2fazgP8SmADga8u30BEigBfAp8qpeyrN3KKw+uNZdeKlTEWNNHkTON7IrcAygVRtGx4aLKOIoMRqVTb2hkw/WFYMzUy/UUr5RtA2Vq2NB2uw38T6CQi24E7fe8RkTgR+ci3zUPAzcAgEVnj+3GRgEcYnNwHk++Fr59x2hL3c3KfEVFGgpQDhnTCfeNiW84iKBFMsyVthPnPQO2b4Lo+kes3Guk1Cbq9ZUvTYVXpKKWOA3fk8HkC8Ljv9RRgSjj9uBK/KFp2Ftz8nNPWuJ9vnoUzR2Hoj/b3Va4uPBkfW4tmhIXNEX5aihZFcwm6lCS/aFG00LE7dXBks1EWm3leO3szRKpU9uvfa1G0UPhyOHz7J1ua1rfa/LB+thZFCxmbnUt6qhFFpp0yZtJqxxKcmu3g2a32V5W1Gw51btKiaGY5shlK2nP+aoefH6peB60GwB1/cdqSKMOmCF8p+GokJO+EgfO1szdLoaJQKqeKa4s4d9KQFqnZxvjRhIA7q3Rii8zzhnOp0AC6/xcKFnbaoujBzvTB8g8MZdI7XobaN9rXj9c4sdeYmJa82/q2Tx+F9zvAr+9a37bXEbEt/akjfLMoBV8OhYJF4f4PnNchiTbaDbdntuu5k7D0b9CwG3TU1VIhcSoR/vd3uKYjlKtjXbt+UbRzyVDvNuvajRn0AijOEz/eEES78xXt7PODXRd+8bIweBGUqqqPS76xOJr0i6L1GGvIJ2hCo0pzKFrKlqa1wzfD/hUXRdFu+J3T1kQnx7YbpaxWqSJmZ8H2xdCwC1RqbE2bsYYdN8it3waIojkohx3NdLcvDaZz+ME4c8yYtalF0cJj8cswd7h17f34pjFrc8+v1rUZq1iZL05PMap/tCiaK9ERfjCO7zSiyd7TtChauFjlV7Yvhp/eguv76ZWrwsKGCL/FQ9DsQR0YhcPsx6BwCejxnuVNa4cfjFrt4Hdr9USesBEs8fgn98GcIYYo2t1aFC0saraFPx2AQhac29++YOSer39EO/twObnPthy+PjK5sX0J/PJv43FXO/vwsaLULDvbWHMgO0uLollBgYJQtGT4UgdrZ8CysXB0izV2aXRZZkQ5uQ/mPA6lqkG7YdqxuIUCBeDGUSAFtJyFFZzYA8vHQ9yjxtyS/OAXRbvmRrj9z5aaF7tY9DScA9rhX05m+sUo8uHPtLO3ihuehvOn8//3/lmbje+xzqZYJzXJiMzr354/h39BFK00PDhRi6JZhY3lxTqlczmLXjTWlOwxVkeRVlKrHdS/QljVHEc2w79bwKYrllvQhEO4jmXLN8ZTwoNaFM1SarSBaq1saVrfkgM5ug0SPoYOI6FJd6et8RaH18P5s4bjDwW/KFqhoka5n8Z68ps9uP4RwzlVqG+pOTHPXa/b1rR2+IFUvBYGL4GqLZy2xHv87y1j8tWTIaxweUEUbRcMnGev0FdM4o/wQ/T4+1cY4ynVW2tnH2XolA5A+mnY9aPxukZrLYpmGyE6Fi2KZi/+lE4oFSGnj8LM/jB3hDHOpbGez/sZ44g2oCN8pWDeU7B5Hjy92ra1JGOe/OSLM84a6w101HIWtlAjDl45ZX777Cz4YjCcOwF9ZxllnRrrOXPctv9tWA5fRMoBnwO1gT3AQ0qpE7lsWxrYBMxVSo0Mp19LiR8PG+cYUaR29vYSam3xTc8atfdaFM0dLH0Ddv9Pi6LZjYurdEYD3yulGgDf+97nxmvAT2H2Zy0XRNG6QMffO22NxzF5EmdnwezBsHOp8V7P2rSPE3uNp9vDG4Jvu28Z/Py2FkWLFDZNvAr3auoBTPa9ngzcl9NGItIaqAx8F2Z/1pF2KkAU7QPtWOzmplFw77+Db/fjm7BhNpzab79Nsc7ZY7DqU0g5EHzbGm2gy9+1KFpEcK8efmWl1CHf68MYTv0SRKQA8E+gH3BnXo2JyFBgKECtWjanV4qWhht/b5zIWhTNfqpeF3ybQFG0VgPst0ljkFc0mZFmKGCWrATtLVQ71eROnZuck1YQkSVATvVwLwa+UUopEcnJyhHAAqVUogTJTSmlxgPjAeLi4mxaAJWLszbbDrGtC81lJK6EtJO5T766IIrWTEeREcNEJPntaNi6EJ5cblwzGvu5Na/MeHgEdfhKqVyjchFJEpGqSqlDIlIVOJLDZh2Am0RkBFASKCIip5VS9u1VXmxfYsiP9p9jVCloIsOysXBoLTy1MufvV072iaJ9CkVKRNa2mCeX2GrtDFg5yVg6Ujt7TxBuSmceMBB40/f7irnvSqm+/tciMgiIc8zZn9xviKKVrQmVmjhiQkyT12Pq7S/B9X20nEUkKVAIipQCyaEEMGmTFkVzimm9ITsT+s22vOlwRyrfBDqJyHaM/PybACISJyIfhWucpWSmG4O0Oop0iFzSB9u+g+TdRimadvaRpWoLeCERru186edpKcbkKi2K5gznT4cnNJgHYR1JpdRx4IqkrFIqAXg8h88/AT4Jp8984xdFe+gz7Vgc47II/8hm4yZc5xboM8MZkzRXorKhYiNoP0KLojmBFWtH5EJs1CJmZ0PWeS2K5iSXD9j7RdGKlDRXrqmxnhN7jfGsAwHjKkoZ+freU/XykY7h3rLM6KBAAWMl+Oxspy2JXW55/uJjql/OInknDNCiaI6Rdgo2fAFNH4DqGJOrFv8Fek0y5qdoHERH+KGTnmoMgPhnEurJVc5RoQFUa2m8XjsdNn5pyFnUuclZuzSAMkTRZg2C00nGAtoa52jQCRp0Dr5dPvBuhK8UzP8dbF8EHZ502hrN3v+DM0ehSQ9ocp8xMNh2qNNWxTb+NFugKNrgxboE02lueMq2pr3r8OMnGI+rd/xFR5FuYOUnsPVbqHeHsXC2nrXpHn78m7EAeff39FoQHsebOY79K2DRCz5RtGectkYDRuVH+imY1FXrqLuFgkWgWFnD2bfsB636O22RBmBKT/hYp3TMs2ysFkVzG0e3Gr/bPK511N1CxYYweq+xAJA+Ju6hSEnbmhZlU71nuMTFxamEhIT8/XFWhqEAeHVtS23ShMHunwxphQ4jtb69RmMjIrJSKZWjbow3I/yChbWzdxt1bjZ+NBqNY+h8h0aj0cQI2uFrNBpNjKAdvkaj0cQI2uFrNBpNjKAdvkaj0cQI2uFrNBpNjKAdvkaj0cQI2uFrNBpNjODambYichTYG0YTFYBjFpnjJF7ZD9D74la8si9e2Q8Ib1+uUUpVzOkL1zr8cBGRhNymF0cTXtkP0PviVryyL17ZD7BvX3RKR6PRaGIE7fA1Go0mRvCywx/vtAEW4ZX9AL0vbsUr++KV/QCb9sWzOXyNRqPRXIqXI3yNRqPRBKAdvkaj0cQInnH4ItJLRDaKSLaI5FrOJCJ7RGS9iKwRkXwuqWUfIexHFxHZKiI7RGR0JG00i4iUE5HFIrLd9/vqXLbL8h2PNSIyL9J25kWw/7OIFBWRz33fLxeR2pG3Mjgm9mOQiBwNOA6PO2FnMERkoogcEZENuXwvIvKubz/XiUirSNtoFhP7cquInAo4Ji+H3alSyhM/QGOgIfAjEJfHdnuACk7bG85+AAWBnUBdoAiwFmjitO052PkWMNr3ejTw91y2O+20rfn9PwMjgA98r3sDnzttdz73YxDwntO2mtiXm4FWwIZcvu8GLAQEaA8sd9rmMPblVuBrK/v0TISvlNqslNrqtB3hYnI/2gI7lFK7lFLngRlAD/utC5kewGTf68nAfQ7akh/M/J8D93E2cIeI6xbtjZbzJShKqZ+A5Dw26QF8qgyWAWVFpGpkrAsNE/tiOZ5x+CGggO9EZKWIDHXamHxSHdgf8D7R95nbqKyUOuR7fRionMt2xUQkQUSWiYibbgpm/s8XtlFKZQKngPIRsc48Zs+Xnr40yGwRqRkZ0ywnWq4Ns3QQkbUislBEmobbWFQtYi4iS4AqOXz1olLqK5PN3KiUOiAilYDFIrLFd6eNGBbthyvIa18C3yillIjkVgN8je+Y1AV+EJH1SqmdVtuqyZP5wHSlVLqIDMN4arndYZtinVUY18ZpEekGzAUahNNgVDl8pdSdFrRxwPf7iIh8ifG4G1GHb8F+HAACI7Aavs8iTl77IiJJIlJVKXXI91h9JJc2/Mdkl4j8CLTEyDk7jZn/s3+bRBEpBJQBjkfGPNME3Q+lVKDNH2GMv0Qjrrk2wkUplRLweoGIvC8iFZRS+RaIi6mUjohcJSKl/K+BzkCOI+QuZwXQQETqiEgRjMFCV1W3+JgHDPS9Hghc8fQiIleLSFHf6wpAR2BTxCzMGzP/58B9fBD4QflG3FxE0P24LM/dHdgcQfusZB4wwFet0x44FZBWjCpEpIp/PEhE2mL46/CCCadHqi0c8b4fI1+XDiQBi3yfVwMW+F7XxahQWAtsxEihOG57qPvhe98N2IYRCbtuP3w2lge+B7YDS4Byvs/jgI98r28A1vuOyXpgsNN2X7YPV/yfgTFAd9/rYsAsYAcQD9R12uZ87sfffNfEWmAp0Mhpm3PZj+nAISDDd50MBoYDw33fCzDWt5/ryaNiz+kfE/syMuCYLANuCLdPLa2g0Wg0MUJMpXQ0Go0mltEOX6PRaGIE7fA1Go0mRtAOX6PRaGIE7fA1Go0mRtAOX6PRaGIE7fA1Go0mRvh/ubOhOfplaZYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mgiD8_ci3mzG"
      },
      "source": [
        "Example of how to compute a pairwise distance matrix\n",
        "$$\n",
        "  D_{i,j} = \\|x_i-y_j\\|^2\n",
        "$$\n",
        "efficiently using the `None` keyword."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_abJlmdr29Ql"
      },
      "source": [
        "def distmat_square(X,Y):\n",
        "  return torch.sum( bc_diff(X[:,None,:] - Y[None,:,:])**2, axis=2 )"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "id": "wWJaya0J2vbv",
        "outputId": "fed0fb35-825d-41c6-ffd1-9f0beccae609"
      },
      "source": [
        "plt.imshow( distmat_square(t[:,None],t[:,None]) )"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f863a196290>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2da4wk13Xff6equntmKHLJlZbMarkWd03GgWDA0QN6QEEgiDYiMYL5hZYoGQqlMNgvki07DizS+WAHcAAJMCwzcEBnYcaQBMEriRZCQhGs2JT0IYDNSJQEyxZNa/kSd8XHklw+Z3q6u+rkQ93uqe6uW9Mz04/bXecP9HR3VXXVqVO3fnPuPadvi6piMpnqq2jRBphMpsXKIGAy1VwGAZOp5jIImEw1l0HAZKq5DAImU801EwiIyHtF5CEROSsit83iGCaTaTqSadcJiEgM/BPwS8A54DvAh1T1R1M9kMlkmopmEQm8DTirqo+oagc4A9w4g+OYTKYpKJnBPo8BTxTenwPePrqRiJwCTgHExG/Z4LLCyvyPiPQ3LnmARu51JKiARvm2GuWvVYCR10QKkRJFShJlJFFKM0ppRj1akj+aAgpsq9DWBp0soaMx3Sx/ZJmQZQKpgAqSARlIBqLDryVTGHsuezB4rbjXVSr6SHw+ElQk90H/2fkK55PcZwx8lvsHJFLiKCOOsoF/mpLSirq0SEkkoqsZ25q4R8P5J6KXRWRZhGYCmez4I8P5YPS1840qZJr70L0e+IcdP+ng/S4tUcS5abjdDPkn2lmev3Y+8rQfHWlDcZSRFHyUt6EuLdd0txW2tcG2Jnk7ymJ6WUSaRXkbco9BG3I+GXqtOuSn/L2//fja0MtcfFZVj4y6aRYQmEiqeho4DXCZHNa3y/U7K0WQOAaJkEYCcYwkCSQJ0mxAI0GbjfzRislaCVkrJm3FpGtCrxXRWxPSNfLndfLX60q6kaHrKclGj42Nba7Y2OLKjZc5vn6R42vPc13rKU4mz3M8iWhrytnuGg93r+TR7SM80b6CJ7cOcWHzEl58dZ32ZhPdTIi2IuItIdkU4jYkbYi3lGRbSdpK3Fbidkq8nRFt94g6KdLpQaeLbHeh10O7/ecepCmaZmiagmbjQBDJfRPJwDcSx9BIkEYDGs5HrQbaTNBmQtqKcx+tRfTWItKW0FuTgX96a5Cua/7YyJD1lOZGh8suaXN4fZOjGy/xM+vPc6J1gZPNZziZvMLhqMnTaYcfd6/gkU7uo/Nbl/P01qU89+oGr2y26G42ka2YaDMi2SL3U98/bYgHPsqItzPi7ZRou4d0esh2D+l0obvjH7o9tNfLfZOmaKaQpeMNrO+jOEbiCBoN99xEkhhc+6HZIGsmrh3F9NZyH6WtiF7LtaH1fluCdEPprSu6nhJf0mN9Y5tD622u3HiZ16+/xBvWn+Vk8wI/27jAiUYGwKPdiIe7R3ikc4THt17HT7cu45nNS3lhc52tzSbZZkK0GRO3hXiz7x+cX8jb0VZG1MlI2inRdopsp0TbuW+k04VOF+2l0O2gaQbd7lgb+mu9+/Gye3EWEDgPHC+8v9otm1yqaJoiMWi3h5CDrvgshc2H+zQxkLmlxU/gniNSoAdsVhrxPMeTmGsbbSJ5yrtVG8hIClbIyHPRLui7fNS6snME0JT8fIog0Dzc0CxCSAf/EMvOdvToOyp6bfQTuY86NHlp7HMFJa9wVdwELno3eQXo0iQDekSUn/GoPePW9bcuWpl/OkWJx0HQ91Ha/0wXpYHQQWkO7TNip8UM3xC+NgQ9YlJgy3vmADkIchhcKN1CRNl0x1fpX6XisfrHjkjI222/tfVtHvZo051jA6Gbu2LQhsqtnAUEvgNcJyInyG/+m4EP73kvAYHgZNIBDASlMhAcGATAAASpzhIE5Zo6BFS1JyKfAL5B3u7+p6r+Q+WH+v200bC3biCIxEBgINhfRCCyKwjwgGAmYwKq+nXg65N/Ih8DyPsvNQdBpgYCA8H+IgI3IFgJghItbGCwKBkM4hAkCFKNgGctIjAQlNg5OxDsXKXpRAQ+44KAAIA0kvwGDwQEIsPHjyWjGBFk+oz3XCwiMBD4tTsIVIUtJgCBCgnRxBGBT2FAQCRPdUFQICjXTtYAZgyCkdcGgvqAALePahBEIDrxYKFPwUBAkmSnsRsIDAQGAmAHBJWDhb6IgPJ7YlTBQIAkGW7sAYFAtcyFBgIDwXxB4B8snCwi8CkYCEizMdbIQwJBuaygyECwBAVFbrDQp2AgQCMpbeTLAgIrKKqQgWDx6cOKowcCAdBmw9vIDQQF6yx9iIEA9pM+9CkICGgk+Zc58P+3MxAUrLP0IQaCvXcNfAoCAoigrZ1mt4wgmHtBEeU+AgMBI1sXrVw1EPT2EBH4FAYEIiFrJd6LuCwgmHtBEeU+AgMBI1sXrVwFEBQLiirHCArpQ5+CgIAKZC4SWDYQjFYW5rL0oYEgrIKiuKI/EAYEIkgL3YFlA0G5DAQGgvmmD3cbLPQpCAggQromFJvasoHACooMBCF0DapBUK4gIKAR9Fp9Ny4vCMo1x4IiSx9SVxBMFhGUKxwIrPVDltUFgWUNKmQgmG36UEbb0I7CgIDkE4Hu9F0MBFZHUC4DwUEqC8sVBAQYRALFQQwDgXUNymUgOEhB0biCgIBGkK7DsHNWDwRzLyiyiAADQbGgqFxhQGDQHYBVA8HCZyii3EdgIGBk66KVqwCCYtagCgRBQIAo/2GQYWesDgjKZelDA8G8C4rKFQgE8l+92XEK1A0ENh+BRwaCqaQPt0abSkHBQEDXUzct+uqCoKqgyNKHFTIQTCUi8CkICESRkmz06MHKg6BcBgIDweJAEAQEkihjY2ObTTAQWPrQLwPBTEAQCARSDl+S3yJ1BYGlDw0EO1bMFwRBQKAZpRxZf2XQZ64LCCx9aCAIAQSBQKDH8fXhC1oXENi3Dw0EiwZBEBBoSY/ja8+PLa8LCGxiEgPBIkEQDASua5X3gesCAosIDASzB0G5goBAU+BkMh4J9FUXEJTLCooMBNMBgU9BQECB40kEjIOgHyobCKyOwEBwMBD4FAQEtlVoa8rxJKYMBH0ZCKyOwEBwEBA8Ubo0CAi0tcHZ7hrXNtoGAu+Zg01MYiA4OAjGFQQEOlnCw90rgWcMBARUUGQRAXUAwa4QEJHjwOeBq9xRT6vqHSJyGPgScA3wGPABVb0oIgLcAdxA3p4/qqrfqzpGR2Me3T6Sn548xcmk4wWBFRTNuaDIIgJWHQSTRAI94LdU9XsicinwgIj8FfBR4D5V/bSI3AbcBnwKeB9wnXu8HbjTPXvVzWKeaF9RWFINgr7qAoJyWfrQQDAdEOwKAVV9EnjSvX5ZRB4EjgE3Au92m30O+DY5BG4EPq+qCvytiFwuIkfdfkrVzWKe3Do0stRAYCAwEMwDBHsaExCRa4A3AfcDVxVu7KfIuwuQA6I4DHnOLRuCgIicAk4BJEcOcWHzkpIjGgj2UlBkdQQeGQjYHGm6RU0MARF5DfAXwG+o6ktSmMdcVVXKa1+9UtXTwGmA1slj+uKr66XbZVo9WGh1BGDpQwPBbiA48PRiItIgB8AXVfWrbvHT/TBfRI6yE4+eB44XPn61W+ZXKrQ3mxUbVIOgboOF5bL0oYGgGgQ+TZIdEOAu4EFV/cPCqnuBW4BPu+d7Css/ISJnyAcEX6waDwBABd1MaFduZBFBUCCwiIBlA4FPk0QC7wI+AvxQRH7glv0O+c3/ZRG5FXgc+IBb93Xy9OBZ8jb7sd0OIBlEWxEZBwNBX3UFgU1MYiDYsWJyEEySHfi/FXu4vmR7BT6+236HlEG8lTf03UBgdQSB1RFQ7iMwEDCyddHKxYCgXEFUDEoGyeaOQ3ePCCxrYF9DNhDsCQSlbcS5wbtmjhKFuA1FhxoIbGKSUQsMBAeNCDwu8K6ZpzJIBne8gSCXRQQGAqYMAs/pe9fMUZJBvFV22XcHwW51BH3VBQTlsoIiA8GwTUOn7l0zR4lCsj3qlr1EBJY+tPShgaAKBFJRLRQGBDIlaRdvnOmDoC8Dgc1QVEcQ9KpOuWLd/KQQt0duGgPB8oPA6ggICQTe061YNzdJpsTtlPFLujsIsrHBMANBUAVF1jUgFBB4T7Vi3fykEG/3Oy17B0E0kh6zgqLACoosIiAEEHhPs2Ld3CSZEm332DFn7yAYl6UPbT4CA0ERBN5TrFg3P6kSdfonaCBYFAgsfejRioDAe3oV6+YnVaTTKzhkfiCoW/qwqqDIvmtQoRUAgffUKtbNT6rQ6Q5us1zTA4EVFAUWEdhgIYsAgfe0KtbNT6rIdnfoNss1zYjAJiax9GG9QeA9Je+aeUoVesMNfxEgqFvXYFSWPlxtEHhPx7tmnlLQbnfMRbMAgaUPLX1YVxB4T8W7Zp5ykUCZi6YBgtGCIpu8tHqwsBgRWNagQksGAu9peNfMU6pjjX3aILCCor11DWyGotUDgfcUqlfPRwqQpuONPJIZjxFYHYEVFBkIgoBA3tiz8kae6UxBYOnDvc1QZF0Dj5YABF7TvWvmLE1zw32NfLYRgYFgLwVF9jVkj4IHgcds/xnNWVocxAgPBHUbLCyXTUyy1CDwmTzRVvPQ2GhmeCDoy0BgEcEqgSAcCMBCQGDzEewNBPYDJ6sHgrAgAAsBgaUP9wYCSx8uIQg029lg1Myqc1iYJgGBpQ+DyhpY+tCjUECQDX92yMQq++cqkfKLWAWCGacPDQRWR7BKIMDTKwgDAvm/drwXcaGDhQYC+4GT1QCBr3I4DAggeXifhQkCKygKbD4Cyn0EBgJGti5a6TVrl/VzkYhAHBf6LuGBwOoIAksfWh0B0wJBEBBAQJJkZBBj+UBQ96xBLpuYJFQQeM3ZZf18JILEubuXHQQWEYBFBMsFgmAgQCPZuXiEBQIrKNobCKygKEwQeM2osnFuEkEajeGLxwFBMOU6AisoqgbBwmcootxHYCAYPY8xE6rsK0pEYuC7wHlVfb+InADOAK8FHgA+oqodEWkBnwfeAjwHfFBVH9tl59BojF88DggCqyOw9GHfNgOB//ATbNPXJ4EHgcvc+88An1XVMyLyJ8CtwJ3u+aKqXisiN7vtPli5Z9cdgJKLR1gRwbgMBFZZuNwgmAgCInI18G+B/wr8RxER4D3Ah90mnwN+jxwCN7rXAHcDfywiolrxhWYRtNUovS1gCiCwiCCoiMDqCDxaEAgmjQT+CPht4FL3/rXAC6ra/9nzc8Ax9/oY8ASAqvZE5EW3/bPFHYrIKeAUwFrzENrciQT6z1MFAcwUBFZQFFjWgHIfQX1B4D3kbhuIyPuBZ1T1ARF590R7nUCqeho4DXDZpce0DwFYXhBY+jAwEFj6kElAMEkk8C7gl0XkBmCNfEzgDuByEUlcNHA1cN5tfx44DpwTkQQ4RD5A6JWKkLbiIYeuKgj6MhBY+jAUEOwKAVW9HbgdwEUC/0lVf1VEvgLcRJ4huAW4x33kXvf+b9z6b1aOBzgbs1buSgPBqNWrB4K51xFQ7iMwEMDB6gQ+BZwRkd8Hvg/c5ZbfBXxBRM6S3wE377YjjYR0bcfIZQSBFRTtJyKwnzybLwg8h6hcOyJV/Tbwbff6EeBtJdu0gV/Z034FemvjzayvmYHACorqVVBkXYPy3Vcde26KIG0N34owJxBY+jCAiMDqCGYOgtL6DbfrquPOS3kkUN53WcauwbgMBFZZGAAIfLutOua8pFERAlBHEFj6EKygaNYg8Oyy6nhzk0C67l4MtHogsIKiwNKHtRss9Oyu6ljzUh4J9N+tNghsYpLAQFCrwULPrqqOMy9pBOm6z4X1A4F1DcAigumDwLubqmPMTZE6CICBYFh1BYFNTDJ9EHh3UbX/uSmCdKPfaCAYEMz4a8hWRxBYHQHlPoLVAIH341X7npckUmQ9LTRyCAYEVkcQTNfA0ocVmgAE3o/u+RMzUBxlNDc6dGiGCQIwEAQEAksferRPEAQDgcsuafMS1BYElj6cvKDIvnRUoX2AIAgINKOUw+v55a8zCGywMLCsAeU+gtUCQSAQ6HF0Y+e0DASWPgwGBCucPhxsuusWc1BTUn5mfbjBGwgsIrD0YV+zBUEQEGhFXU60LowtrxMIbD6C/UQENjHJNEAQBgRIOdksT/3UCQQ2H8HeQGB1BHsFgWd11WfnpUQiTib+coagQWC/a2C/a+C2Dh0E3lVVn5uXuppxOGpWGho0CKygKJj0oYFg7woCAtua8HTayUOWZQUBGAgWPkZgBUU+EOSDqWdL1gQEgR93rwAuGgjGznhHuxUUWfoQrI5g7xFBIBBo8EjnSvfOQGDpwyUCwbLWERQUBAS6Wcyj20cKSwwENjHJEoFgWesIhs5owepmEee3Lh9ZaiCwysIlKigKPCIoH0wtns2C1csint66tGTNkoPA5iOo13wEAUcEL4+f8kBBQCDLIp57dcOzdslBADOOCCxrYPMRTAYCn4KAgGbCK5utii2WHARWR2AgcFsvEgQ+BQEBMqG72aycAmnpQWCVhVZQ1LdtESDQaHTjgQKBAMhWTJcVB8GMIwKbmMQKiqojgnIFAQHJINqMyGD1QQAzjggMBEGlDyn3EcwXBBJ6dkAySLbyH0gwEFj6cOVAEEL6UMc2GCgICKAQb+WnaCAYtaDMZosIlg4EwaQPxxUEBCSDpA39UzQQjFpQZrNNTLJ0BUWU+wjmBYJyBQOBeGvYTQaCUQvKbLaJSWxiknEL9gqCICCA9iOBmoLA0oc2MYnbehEgCAICkkG8Xf4/sjYgsIKiYMYIVhcE5ZoIAiJyOfCnwM+7Y/574CHgS8A1wGPAB1T1oogIcAdwA7nPP6qq36vcf6Yk7bLTqhkIYKEgqFvWoJ4FReOaNBK4A/hLVb1JRJrABvA7wH2q+mkRuQ24DfgU8D7gOvd4O3Cne/ZLIW73L2D/NAwEiygoqvsYQa45FhQtJH04rF0hICKHgH8NfDS3QTtAR0RuBN7tNvsc8G1yCNwIfF5VFfhbEblcRI6q6pPeY6gSb2funYFg9hGB1REElT6k3EcwHxBMEgmcAC4AfyYivwA8AHwSuKpwYz8FXOVeHwOeKHz+nFs2BAEROQWcAlhrHiLeLn7FwUCwaBD0ZSBYnToCnyaBQAK8Gfg1Vb1fRO4gD/0L9qhK+ZCrV6p6GjgNcGjj9Rpt90a2MBAYCGoGgll2DTL/7TkJBM4B51T1fvf+bnIIPN0P80XkKDujKOeB44XPX+2W+aWKdHolSYyIHC0GgmmBYK8FRXUfI1jFnzwb1a4QUNWnROQJEfk5VX0IuB74kXvcAnzaPd/jPnIv8AkROUM+IPhi1XgAAJkiLhIoA0HeDmsMginXEey1oKivuoBg4TMUUe4jOBgIfJo0O/BrwBddZuAR4GPO7i+LyK3A48AH3LZfJ08PniX368d23bsq0ukO3paBIFeNQWB1BMF0DZY5fVimiSCgqj8A3lqy6vqSbRX4+CT7LXwI3IXty0Aw/zECm48gsPQh5T6C/YHApyAqBkHRbneXi1hcYiCYFQh26xrULX1YVVC0bOlDn8KAgCr0ehNcxOISA8GiIoK6DxbmWj4Q+BQIBIBub/DSQLB4EFhBUWAgmFL6sExBQEBV0V5v5AYxEISePuyrriBYtvShT0FAAPKLCaM3iIEg9PRh3boGy5w+9CkMCCiQjg9iGAiYDASWPgxmsDDk9KFPYUAA0EwpG800ELCQiGBcBoLln5ikXMFAgCxFiQ0EFhEMtgwdBEsXEXgUDgTAQMAUQAAzBYEVFC15QVGJwoIAGAgIHwSWPgwsfUi5j2AyEIQHATAQsPwg6MtAEFIdQfnxw4QAGAgwEBgIis9Fu2B/dQTlCgcCIuMX0UAQHAisoGhvIJh7QVFFROCrFwoDAiIgJY0cDARMAQRWUFSvgqKKiKBMQUBAAIlj/yCGgeDgILD0YTBdg0WlD30KAgKIIHF+wgaC5ekajMtAEDIIfAoGAjQaCPnsQgYCA8GqgGC3gqJ51hH4FAwEJI5QDASrAILdCoqsjgAWkT70KRAIAI0mQsdAwGqAwGYxDgwEFVYEAgFBkhjFQNBXHUBgEQHMEwQ+BQMBmg1n7HxBYL9rsFgQ9FVXEMyzjsCnYCCgDgKLAIH9rsFsv4ZsdQRh1BH4FAQENMojgeFb0LoGfVkdwajVyw+CRXwN2acgIIAIWTMhYvQWNBD0tRpdAwPBoicmKVMgEABtxQMXGgjKLTQQjFq9/CCYW0QQlR0nVxAQ0EjIWvHgkhoIVhsENjHJIiYm8Y8KBAEBBHpr8eBUDATF53ELVwEElj5cRPqwXEFAQCMhXctPwEBQtLTeIOjLQDBbEAQCAUhbw83MQGAgKMpAcEAQVLAgDAgI9FrjzSwUEFhB0XRBYBOT7A0Esyko2lEQECCCdA3KmlkoILCCIpuYZKkLioKPBCLorfv/34QCgp1P1RQEVlAUTNdgz+lDCTwSUClGAgaCoEGA7wqNWlBms4FgsenDcgUBASLorfXfGAjAQFAFgrqlD6sKivbeNRjXRBAQkd8E/oOz+IfAx4CjwBngtcADwEdUtSMiLeDzwFuA54APqupjVfvXCNKNUWcYCOoMAisomkVEUK5dISAix4BfB96oqlsi8mXgZuAG4LOqekZE/gS4FbjTPV9U1WtF5GbgM8AHKw8SKb310dvRQAD1BoFNTDK99OF26fcUck3aHUiAdRHpAhvAk8B7gA+79Z8Dfo8cAje61wB3A38sIqI66sGCIkXXU3qDC2QgKMpAYJWFB00fvnQQCKjqeRH5A+AnwBbwf8jD/xdUtec2Owccc6+PAU+4z/ZE5EXyLsOzxf2KyCngFEDjyCHiS3ruhA0EZaozCCx9OL30YZkm6Q5cQf7f/QTwAvAV4L17OkqJVPU0cBrgNf/8n+n6xjZbsLQgOGhBUcZFjhoIBhotKLLJS6sHC4sRwW5jBKOapDvwi8CjqnoBQES+CrwLuFxEEhcNXA2cd9ufB44D50QkAQ6RDxD6jYgyDq3nTWCZQXDQgqLIIgIrKDpA12C/EcEkEPgJ8A4R2SC/R68Hvgt8C7iJPENwC3CP2/5e9/5v3PpvVo4HAM0o5cqNlwfvlxkE43bWrGsw46nKrI5gfwVFEf5bcJIxgftF5G7ge+743ycP4/83cEZEft8tu8t95C7gCyJylvwq3bzbMZpRj9evvzS0zECwxCCYcWWhpQ/3NkNRJE+RT9lTromyA6r6u8Dvjix+BHhbybZt4Fcm2W9fLenxhvVnx5YbCJYYBPiu0KgFZTbbtw+nXVBUVWIcRMVgS7qcbF4oXWcgMBD4ZYOFe6sjKFcgEICfbZRDAAwEBoIqWUQwKQh8CgICACcaGWAgKJeBwOYjOHhBkU9BQGDb+cxAUCUDgaUPDxoR/LR0aSAQaPBoN+JEIzMQeM8clh4Elj4MIGswrkAgkPBw9whwwUCw6iCwiUkWEhGkJa2xryAg0MkSHukcce8MBCsPAnxXaNSCMpsNBNNJH+4oEAjEPL71usISA4GBoKi9g8AKiiYZI8gVBAR6WcRPty4bWWogMBAUtZ+IwOoIJgFBEBBIs4hnNi8tWWMgMBAUNX0QWNYgEAhkmfDC5rpnrYHAQFCURQT7BYFPQUCATNjabFZssPogsB84mS4IrKBoBAQjhysqGAhkmwlbUJHbXH0Q2A+c2A+czAwEGniKUDKINvMvO1YPYqw+CMbtrBkIrI5ghhFBuYKAABnEbUHFQGAgmH5EMK66gqBcQUBAFOLN3NBUDQQGAosI5gmCMCCQQdKGvqHFiMBf7WQgWHkQ4LtCoxaU2WwFReMgKFcQECCDeKv/Zjgi2PJ9BjAQGAiGtZ+IoE7pw3IFAQFRSNplgxgGAgPB4kHQ1/KDoFxhQCCDuA3DTcdAML7EQGAgOAgIyhUEBFAlaY/eahSeDQRgBUXTBkH9CorKFQQE8oFBzfsFY00nf+4PFlpB0aidNQOBFRTtDwT+ccFQIKBEnf6PJ/sLHSx9WFxSYxBY+nDvIAgdAigk7ZQe7AoCKygqLqkxCPBdoVELymyuKQg8CgQCSrSdksBEIChGBFZHYCBYFAiWKX0oo58tKAgISKbIdj4L2qQgsMHC4hIDwaxAsCoFRRW9gTAgkEcC3cFNNQQC7TuverDQLwOBgaCo/UQEqzUxyaiCgQDd3tBNtQOCsluNwbNVFhaXGAgWBYJl6Br4FAwEpNNFYV8gsK5BvuSgdQQZFzlqIKjQcncNfAoEAkCnOzA2EikHQUUdQREEVkewfxBEyx4RzPhryMtcR+BTIBBQtJcOnXSkWj5GMMFgoaUPR+2sWdfA6ghKQeBTMBCg23E3SwEEpREBGAgMBIvvGiwfCHwKBgKaZu4mGQGBRQQGgkBBsGzpQ5+CgIACdLvu5igBAWWDhVAFApuYpLjEQDC7iGB5QOBTEBDYiQT2AYKKOgJLHxaXGAgWBYJQ0oc+hQEBGBi5dxBUpw+toKi4xEBQ54jAp2AggBYv4nRBYBFBcYmBYFogWLr5CDwKCAKjF7ECBCVZgzgaqSOQ8oigLfDiAesINjW/PGMg0OGLKNlOt2UqINAMzaLpgECmA4J05EctVIWXKYDArRcVN/gzAQgynR0IfHUE0fCS/IDiPjQMgm1RXnJtqD8fQYQSkxHJM0NZg75/UrfvsToC7ftoFxBovjgCRLUUbPsFgeiEtJilRORl4KFF27EHvQ54dtFGTKhlshWWy95lshXgDap6ZHRhKJHAQ6r61kUbMalE5LvLYu8y2QrLZe8y2Vql8YjRZDLVSgYBk6nmCgUCpxdtwB61TPYuk62wXPYuk61eBTEwaDKZFqdQIgGTybQgGQRMpppr4RAQkfeKyEMiclZEbgvAnuMi8i0R+ZGI/IOIfNItPywifyUiP3bPV7jlIiL/zdn/dyLy5gXYHIvI90Xka+79CRG539n0JRFpuuUt9/6sW3/NAmy9XETuFpF/FJEHReSdofpWRH7TtYG/F5E/F5G1kH27Xy0UAiISA/8deB/wRuBDIvLGRdpEXsj1W6r6RuAdwMedTbcB96nqdcB97j3ktl/nHqeAO+dvMjrTHKQAAAKnSURBVJ8EHiy8/wzwWVW9FrgI3OqW3wpcdMs/67abt+4A/lJV/wXwC+R2B+dbETkG/DrwVlX9efKCxpsJ27f7k6ou7AG8E/hG4f3twO2LtKnExnuAXyKvaDzqlh0lL3AC+B/AhwrbD7abk31Xk9847wG+Rl45+iyQjPoY+AbwTvc6cdvJHG09BDw6eswQfQscA54ADjtffQ34N6H69iCPRXcH+o7u65xbFoRcSPcm4H7gKlV90q16CrjKvV70OfwR8Nvs/Ozsa4EXVLU/i0TRnoGtbv2Lbvt56QT5lzL+zHVf/lRELiFA36rqeeAPgJ8AT5L76gHC9e2+tWgIBCsReQ3wF8BvqOpLxXWa437huVUReT/wjKo+sGhbJlQCvBm4U1XfBLzKTugPBOXbK4AbycH1euAS4L0LNWpGWjQEzgPHC++vdssWKhFpkAPgi6r6Vbf4aRE56tYfBZ5xyxd5Du8CfllEHgPOkHcJ7gAuF5H+90KK9gxsdesPAc/NyVbI/3OeU9X73fu7yaEQom9/EXhUVS+oahf4Krm/Q/XtvrVoCHwHuM6NuDbJB17uXaRBIiLAXcCDqvqHhVX3Are417eQjxX0l/87N5L9DuDFQmg7U6nq7ap6tapeQ+67b6rqrwLfAm7y2No/h5vc9nP7r6uqTwFPiMjPuUXXAz8iQN+SdwPeISIbrk30bQ3StwfSogclgBuAfwIeBv5zAPb8K/Jw9O+AH7jHDeT9u/uAHwN/DRx22wt5huNh4Ifko8mLsPvdwNfc65PA/wPOAl8BWm75mnt/1q0/uQA7/yXwXeff/wVcEapvgf8C/CPw98AXgFbIvt3vw8qGTaaaa9HdAZPJtGAZBEymmssgYDLVXAYBk6nmMgiYTDWXQcBkqrkMAiZTzfX/AU7lHVsaJ/6cAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2G4vkyh8z7_0"
      },
      "source": [
        "A more memory efficient way (specially in high dimension $d$) way to compute the distance by expanding the squares."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3F7QiKf0zeo2"
      },
      "source": [
        "def distmat_square2(X, Y):\n",
        "  X_sq = (X ** 2).sum(axis=-1)\n",
        "  Y_sq = (Y ** 2).sum(axis=-1)\n",
        "  cross_term = X.matmul(Y.T)\n",
        "  return X_sq[:, None] + Y_sq[None, :] - 2 * cross_term"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LiWlq3PtP1e9",
        "outputId": "8a30b20a-d157-4586-be8e-2823228e9437"
      },
      "source": [
        "t0 = time()\n",
        "distmat_square(X, Y)\n",
        "print(time() - t0)\n",
        "\n",
        "t0 = time()\n",
        "distmat_square2(X, Y)\n",
        "print(time() - t0)"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.027776718139648438\n",
            "0.04250526428222656\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dBoeQ8xT3VrZ"
      },
      "source": [
        "# Maximum Mean Discrepencies"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h-m3t5it57Ds"
      },
      "source": [
        "When training ML model with point clouds, it is important to be able to compare to pair of clouds. A simple loss function is to use an MMD norm\n",
        "$$\n",
        "    \\text{MMD}(X,Y) \\triangleq\n",
        "    \\frac{1}{n^2} \\sum_{i,i'} k(x_i,x_{i'}) \n",
        "    +\n",
        "    \\frac{1}{m^2} \\sum_{j,j'} k(y_j,y_{j'})\n",
        "    -2\n",
        "    \\frac{1}{nm} \\sum_{i,j} k(x_i,y_{j}).\n",
        "$$\n",
        "We use here the energy distance kernel, so that MMD$(X,Y)$ is a dual Sobolev norm between the discrete measure $\\frac{1}{n}\\sum_i \\delta_{x_i}$ and $\\frac{1}{m}\\sum_i \\delta_{y_j}$\n",
        "$$\n",
        "  k(x,y) = -\\|x-y\\|.\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_oQaky5t55zK"
      },
      "source": [
        "def kernel(X,Y):\n",
        "  return -torch.sqrt( distmat_square(X,Y) )\n",
        "\n",
        "def MMD(X,Y):\n",
        "  n = X.shape[0]\n",
        "  m = Y.shape[0]\n",
        "  a = torch.sum( kernel(X,X) )/n**2 + \\\n",
        "      torch.sum( kernel(Y,Y) )/m**2 - \\\n",
        "      2*torch.sum( kernel(X,Y) )/(n*m)\n",
        "  return a.item()"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qxnErW9v6oMr",
        "outputId": "f49df341-a7a7-4b79-a811-caad4073a63e"
      },
      "source": [
        "print( MMD(X,X) ) # should be 0\n",
        "print( MMD(X,Y) ) # should be >0"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.0\n",
            "6.395578384399414e-05\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CihPkBXe3blG"
      },
      "source": [
        "# Gradient flow with Pytorch"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m16xCozj4I1s"
      },
      "source": [
        "We compute a speed vector field to advances the particlesusing an interaction kernel of the form\n",
        "$$\n",
        "  v(X)_i = \\frac{1}{n} \\sum_j \\psi(\\|x_i-x_j\\|^2) (x_j-x_i)\n",
        "$$\n",
        "Note that the evolution\n",
        "$$\n",
        "  \\dot X = -v(X)\n",
        "$$\n",
        "is the Wasserstein gradient flow of the energy\n",
        "$$\n",
        "  \\sum_{i,j} \\phi(\\|x_i-y_j\\|^2 ).\n",
        "$$\n",
        "when defining $\\psi(r)=4\\phi'(r)$.\n",
        "\n",
        "If $\\phi$ is decreasing (resp. increasing), the flow is repulsive (resp. attractive)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P37QnQCe4H-e"
      },
      "source": [
        "sigma = .1;\n",
        "def psi(r):\n",
        "  return torch.exp( -r/(2*sigma**2) )"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "neZBQfE53EqM"
      },
      "source": [
        "def Speed(X):\n",
        "  return 2/X.shape[0] * 1/sigma**2 * torch.sum( psi(distmat_square(X,X))[:,:,None] * bc_diff( X[:,None,:] - X[None,:,:] ), axis=1 )"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nsipu8UblUt5"
      },
      "source": [
        "Discretize the evolution $\\dot X = -v(X)$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "v69IX-Ov0L9B",
        "outputId": "0e29574d-7e9e-4899-fb6d-ad1e75c78a1e"
      },
      "source": [
        "if boundary=='no':\n",
        "  tau = 1/500 # time step\n",
        "else:\n",
        "  tau = 1/200\n",
        "niter = 200\n",
        "save_per = 10 # periodicity of saving\n",
        "Zsvg = torch.zeros((n,2,niter//save_per)) # to store all the intermediate time\n",
        "Z = X\n",
        "for it in progressbar.progressbar(range(niter)):\n",
        "  if np.mod(it,save_per)==0:\n",
        "    Zsvg[:,:,it//save_per] = Z.clone().detach() # for later display\n",
        "  Z = bc_pos( Z - tau*Speed(Z) )"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "100% (200 of 200) |######################| Elapsed Time: 0:00:15 Time:  0:00:15\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S8vh1YvTlaB0"
      },
      "source": [
        "Display the evolution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 301,
          "referenced_widgets": [
            "1e24c9aa62954c4099b28c89d46d55ca",
            "ebd017bc9dca41099c4fa560e9f68922",
            "d34c46ec2e8f4170af80954555f8a2b7",
            "1dbff5ccf90247b69e19cb78cbe86e25",
            "2d5a3a03f7cb421cbb0cae446cbdfdb1",
            "f17ea6f4ce5245009f234130c9273294",
            "b974db90bbbf44c482c1dd12f61013a3"
          ]
        },
        "id": "Lwvb9cOCgeaA",
        "outputId": "15abf4e5-35de-4afa-f49b-cc187cdc4c63"
      },
      "source": [
        "import ipywidgets as widgets\n",
        "@widgets.interact(t=(0,niter//save_per-1))\n",
        "def display_frame(t=0):\n",
        "  s = t/(niter//save_per-1)\n",
        "  plt.scatter(Zsvg[:,0,t], Zsvg[:,1,t], color=[s,0,1-s])\n",
        "  plt.axis('equal')\n",
        "  plt.axis([0,1,0,1])"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "1e24c9aa62954c4099b28c89d46d55ca",
              "version_minor": 0,
              "version_major": 2
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=0, description='t', max=19), Output()), _dom_classes=('widget-interact',…"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1-lalpcK3ind"
      },
      "source": [
        "# Computing Gradient with Auto-diff\n",
        "\n",
        "Instead of computing \"by hand\" the gradient of the interaction energy, one can directly rely on Pytorch auto-diff functionality. This simplifies coding and it reduces bugs in code."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "akQKlOnh3xpu",
        "outputId": "1de5abb3-5ee7-4775-8e04-e0bec140ef32"
      },
      "source": [
        "X.requires_grad = True\n",
        "L = -1/X.shape[0] * torch.sum( psi(distmat_square(X,X)), axis=(0,1) )\n",
        "[g] = torch.autograd.grad(L, [X])\n",
        "# compare with the \"by hand\" computation\n",
        "print( 'Difference \"hand\" vs. pytorch\" : ' + str( torch.norm( g-Speed(X) ).item()  / torch.norm( g ).item() ) )"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Difference \"hand\" vs. pytorch\" : 1.7333064950325788e-07\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2DVdojtPpSyd"
      },
      "source": [
        "# Computations using Keops\n",
        "\n",
        "In order to speed up (by a factor 10 to 100) the evaluation of the kernel, it is possible to use [Keops](https://www.kernel-operations.io/keops/index.html), which allows one to define \"lazy tensor\". This both reduces the memory footprint (by evaluating the kernel on the fly) and accelerate the computation by a careful mapping of the tensor entries evaluations on the tiles of the GPU."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D5LIYwHWpVSc"
      },
      "source": [
        "!pip install pykeops[colab] > install.log"
      ],
      "execution_count": 16,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "C8iPf7bTpdd9",
        "outputId": "24ab29e7-60c0-48e3-8bf2-00b55f23be64"
      },
      "source": [
        "import pykeops\n",
        "import pykeops.torch as keops\n",
        "pykeops.clean_pykeops()          # just in case old build files are still present\n",
        "pykeops.test_torch_bindings()    # perform the compilation"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cleaning /root/.cache/pykeops-1.5-cpython-37/...\n",
            "[pyKeOps] Initializing build folder for dtype=float32 and lang=torch in /root/.cache/pykeops-1.5-cpython-37 ... done.\n",
            "[pyKeOps] Compiling libKeOpstorch180bebcc11 in /root/.cache/pykeops-1.5-cpython-37:\n",
            "       formula: Sum_Reduction(SqNorm2(x - y),1)\n",
            "       aliases: x = Vi(0,3); y = Vj(1,3); \n",
            "       dtype  : float32\n",
            "... \n",
            "[pyKeOps] Compiling pybind11 template libKeOps_template_574e4b20be in /root/.cache/pykeops-1.5-cpython-37 ... done.\n",
            "Done.\n",
            "\n",
            "pyKeOps with torch bindings is working!\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YXtRegSD-A8O"
      },
      "source": [
        "Compute the gradient of a kernel interaction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jN_iunlV4mRZ",
        "outputId": "e1426001-4794-4600-a990-b9131364eaaf"
      },
      "source": [
        "X.requires_grad = True\n",
        "D = keops.Vi(X) - keops.Vj(X)\n",
        "if boundary=='per':\n",
        "  D1 = (D-.5).mod(1.0)-.5 # for periodic BC\n",
        "else: \n",
        "  D1 = D\n",
        "D2 = ( D1 ** 2 ).sum( dim=2 )\n",
        "K = ( -D2 / (2*sigma**2) ).exp()\n",
        "L = -1/X.shape[0] * (K.sum(dim=1)**1).sum()\n",
        "[g] = torch.autograd.grad(L, [X])\n",
        "# compare with the \"by hand\" computation\n",
        "print( 'Difference \"hand\" vs. \"keops+pytorch\" : ' + str( torch.norm( g-Speed(X) ).item()  / torch.norm( g ).item() ) )"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[pyKeOps] Compiling libKeOpstorchea57037698 in /root/.cache/pykeops-1.5-cpython-37:\n",
            "       formula: Sum_Reduction(Exp((Minus(Sum(Square((Mod(((Var(0,2,0) - Var(1,2,1)) - Var(2,1,2)), Var(3,1,2), IntCst(0)) - Var(4,1,2))))) / Var(5,1,2))),0)\n",
            "       aliases: Var(0,2,0); Var(1,2,1); Var(2,1,2); Var(3,1,2); Var(4,1,2); Var(5,1,2); \n",
            "       dtype  : float32\n",
            "... \n",
            "Done.\n",
            "[pyKeOps] Compiling libKeOpstorch17a5f3fca7 in /root/.cache/pykeops-1.5-cpython-37:\n",
            "       formula: Grad_WithSavedForward(Sum_Reduction(Exp((Minus(Sum(Square((Mod(((Var(0,2,0) - Var(1,2,1)) - Var(2,1,2)), Var(3,1,2), IntCst(0)) - Var(4,1,2))))) / Var(5,1,2))),0), Var(0,2,0), Var(6,1,0), Var(7,1,0))\n",
            "       aliases: Var(0,2,0); Var(1,2,1); Var(2,1,2); Var(3,1,2); Var(4,1,2); Var(5,1,2); Var(6,1,0); Var(7,1,0); \n",
            "       dtype  : float32\n",
            "... \n",
            "Done.\n",
            "[pyKeOps] Compiling libKeOpstorchfdfa2e7b49 in /root/.cache/pykeops-1.5-cpython-37:\n",
            "       formula: Grad_WithSavedForward(Sum_Reduction(Exp((Minus(Sum(Square((Mod(((Var(0,2,0) - Var(1,2,1)) - Var(2,1,2)), Var(3,1,2), IntCst(0)) - Var(4,1,2))))) / Var(5,1,2))),0), Var(1,2,1), Var(6,1,0), Var(7,1,0))\n",
            "       aliases: Var(0,2,0); Var(1,2,1); Var(2,1,2); Var(3,1,2); Var(4,1,2); Var(5,1,2); Var(6,1,0); Var(7,1,0); \n",
            "       dtype  : float32\n",
            "... \n",
            "Done.\n",
            "Difference \"hand\" vs. \"keops+pytorch\" : 2.946057011364191e-07\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7x3ZP-176ONd",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c3b14b63-f2bc-4d3a-c9a1-bf6da9f680b5"
      },
      "source": [
        "Zsvg = torch.zeros((n,2,niter//save_per)) # to store all the intermediate time\n",
        "Z = X\n",
        "Z.requires_grad = True\n",
        "for it in progressbar.progressbar(range(niter)):\n",
        "  if np.mod(it,save_per)==0:\n",
        "    Zsvg[:,:,it//save_per] = Z.clone().detach() # for later display\n",
        "  D = keops.Vi(Z) - keops.Vj(Z)\n",
        "  if boundary=='per':\n",
        "    D1 = (D-.5).mod(1.0)-.5 # for periodic BC\n",
        "  else: \n",
        "    D1 = D\n",
        "  D2 = ( D1 ** 2 ).sum( dim=2 )\n",
        "  K = ( -D2 / (2*sigma**2) ).exp()\n",
        "  L = 1/X.shape[0] * (K.sum(dim=1)**1).sum() # There is a bug, I needed to add **1 here !!\n",
        "  [g] = torch.autograd.grad(L, [Z])\n",
        "  Z = bc_pos( Z + tau*g )"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "100% (200 of 200) |######################| Elapsed Time: 0:00:01 Time:  0:00:01\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 301,
          "referenced_widgets": [
            "43db062d32ac4b39b5f554d12c335620",
            "a4dc0e31bfea45bcb9f064488a7c40ed",
            "7b0b70baaaaa4da68b7e71469ee49731",
            "a6fb1ecdfa004558849b54cd293b6645",
            "197250af0ee648579e97f82dd8006c14",
            "0071e8e86e8949f088f9d82f1e5fa3e6",
            "242932679844407388684f5af4adcdda"
          ]
        },
        "id": "2PtiN6706_GG",
        "outputId": "33f30928-e7ae-4c60-ad4a-be1868bc6091"
      },
      "source": [
        "import ipywidgets as widgets\n",
        "@widgets.interact(t=(0,niter//save_per-1))\n",
        "def display_frame(t=0):\n",
        "  s = t/(niter//save_per-1)\n",
        "  plt.scatter(Zsvg[:,0,t], Zsvg[:,1,t], color=[s,0,1-s])\n",
        "  plt.axis('equal')\n",
        "  plt.axis([0,1,0,1])"
      ],
      "execution_count": 20,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "43db062d32ac4b39b5f554d12c335620",
              "version_minor": 0,
              "version_major": 2
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=0, description='t', max=19), Output()), _dom_classes=('widget-interact',…"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XrEJTQa08tqg"
      },
      "source": [
        ""
      ],
      "execution_count": 20,
      "outputs": []
    }
  ]
}