jcheong0428/tds-matrices.ipynb

## tds-matrices.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "TDS-matrices",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyODZLVioJOcFmhXnoHETH9R",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/jcheong0428/ea4b62f6e2bb56e2f04102aa580004c5/tds-matrices.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Measuring similarity between two correlation matrices. \n",
        "\n",
        "\n",
        "Author: [Jin Hyun Cheong](https://jinhyuncheong.com)"
      ],
      "metadata": {
        "id": "S3P3M72s9k9o"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import seaborn as sns\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scipy.stats as stats\n",
        "from scipy.spatial.distance import squareform\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# n is the number of subjects or items that will determine the size of the correlation matrix. \n",
        "n = 20 \n",
        "data_length = 50\n",
        "nr = .58\n",
        "\n",
        "# Random data\n",
        "np.random.seed(0)\n",
        "m1_u = np.random.normal(loc=0, scale=1, size=(data_length, n))\n",
        "m2_u = nr*np.random.normal(loc=0, scale=1, size=(data_length, n)) + (1-nr)*m1_u\n",
        "\n",
        "m1 = pd.DataFrame(m1_u).corr()\n",
        "m2 = pd.DataFrame(m2_u).corr()\n",
        "\n",
        "f,axes = plt.subplots(1,2, figsize=(10,5))\n",
        "sns.set_style(\"white\")\n",
        "for ix, m in enumerate([m1,m2]):\n",
        "  sns.heatmap(m, cmap=\"RdBu_r\", center=0, vmin=-1, vmax=1, ax=axes[ix], square=True, cbar_kws={\"shrink\": .5}, xticklabels=True)\n",
        "  axes[ix].set(title=f\"m{ix+1}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 270
        },
        "id": "0JK_Fk_r-FPh",
        "outputId": "9ec0d1cd-ee3c-4c01-f428-ec9d48d4a10c"
      },
      "execution_count": 71,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAD9CAYAAABz/0lcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhTZdoG8DtLS/fSNQUsIAgWFUQFAXGKFMpi6QCFWsbBEQGZcRyQTcfiguIAOqKizufC4D6OIoh0oIOCVayDLKIO8CnIhxYpS1toS/c2TXK+P7yM1Cb0yUnSnAP377pyXTR9OOdNmt59c3LO8xoURVFARERERGLGQA+AiIiISG84gSIiIiLyECdQRERERB7iBIqIiIjIQ5xAEREREXmIEygiIiIiD3ECRUREROQhTqDI6YEHHsDo0aORkpKC9evXB3o4REQiRUVFuOOOOzB48GBce+21mDFjBr7//vtAD4vOc5xAkVNKSgoeeughXHbZZYEeChGRWE1NDdLS0vD+++9j+/bt6Nu3L/74xz8Gelh0nuME6jyXlpaG1atXIzMzE/3798eiRYtw+vRpzJw5E1dddRWmTZuGqqoqAMBvf/tbDBkyBB06dAjwqImI5PnVr18/ZGdno2PHjggKCsK0adNQVFSEysrKQD8EOo9xAnUB2LJlC1555RV88MEH+Pjjj3H77bdj/vz52LlzJxwOB954441AD5GIyCU1+bVnzx4kJCQgJiYmACOmC4U50AMg/5s6dSri4+MBAAMGDEBsbKzzY7r09HTs2LEjkMMjInLL0/wqKSnBww8/jHvvvbfdx0oXFh6BugD8FD4A0KFDhxZfh4SEoL6+PhDDIiJqkyf5VVFRgenTp+Pmm2/GuHHj2nWcdOHhBIqIiHSvqqoK06dPR1paGu64445AD4cuAPwIj5ysVisURYGiKLDZbGhqakJQUBCMRs6ziUi7amtrMWPGDFx99dVYuHBhoIdDFwj+ZSSnGTNmoF+/fvjqq6/wwAMPoF+/fvj8888DPSwionPaunUr9u/fj/Xr1+Oqq65y3k6cOBHoodF5zKAoihLoQRARERHpCY9AEREREXmIEygiIiIiD3ECRUREROQhTqCIiIiIPMQJFBEREZGHvOoDVVhYiKVLl8LhcCA7OxuzZs06Z/0fDN1F211Zf0BUZ3DYRHWwN4vKrOZQUV2jTXbhYuxR2RIptm7XiOoAoNEoW+g3tLlGVKeYgmV1Ztl+g09+Lapr2vepqM6YPlNUBwCGhipR3XFEi+pO18teX1dGNIjqHCGy/QJASFi4uJbU8yTDpPn1P0c3iurerYwT1U0yHxLVeZIjilEW/abqElGdUfi797mhq6juokhZLr1/uFxUd2s3h6gOAOyRiaI6U02ZqM4RKvu9dwSHieoqGuyiupgQk6iuqkm2PQBIPCFra2O67AbxNvVM9REou92OJUuWYPXq1cjPz8emTZtw+PBhX46NiMhvmGFE5A3VE6h9+/ahW7duSE5ORnBwMDIyMlBQUODLsRER+Q0zjIi8oXoCVVpaiqSkJOfXFosFpaWlPhkUEZG/McOIyBs8iZyIiIjIQ6onUBaLBSUlP59gWFpaCovF4pNBERH5GzOMiLyhegLVt29fHDlyBMXFxbBarcjPz0daWpovx0ZE5DfMMCLyhuo2BmazGQ8++CBmzpwJu92OSZMmoVevXr4cGxGR3zDDiMgbXvWBGjZsGIYNG+arsRARtStmGBGpxZPIiYiIiDxkUBRF1lbbhdzcXGzbtg1xcXHYtGlTm/WNDbKOzXPD+ojq/vbfVaK62p7Xi+rC9m8W1VVdPlZUF22vFtX90CzvOt3zzD5RneKQdd41mINEdfZI4cm1imy/to5dRHVGa71svwCKm2SPpatJ1qXdYLOK6g7ZY2XbM4jKAACXJUXJi0kVT/PLXrxftN07u2aK6pZVfyOqiy6XdSI32GWvVwA4k3iFqC6y8bSozh4u66pubJL97kkzsVtQnaiu1iz/fTpU0SiquzpatlLBGYPssTTaZX+Ku5TsEdV50pleymqQfWgVGSZb1UPvvDoClZWVhdWrV/tqLERE7Yb5RUTe8GoCNXDgQERHy9f3IiLSCuYXEXmD50AREREReYgTKCIiIiIPcQJFRERE5CFOoIiIiIg85NUEav78+ZgyZQqKioqQmpqKtWvX+mpcRER+xfwiIm941Yn8ySef9NU4iIjaFfOLiLzBj/CIiIiIPORVJ3JPNdXJutCaDu8Q1f2p/yxR3bNlhaI6OOyiMiVY1mXV0CTrkuuIiBfVAUADZN22o0pkHculrF2uFNUV18i68yaEmbwZjkshJlmrb5vwFd8k7Aws3C0qGmWvLwDoER8prqX28c6+E6K6kRd3FNUtirpMVLey/oCoztAsW+kBAJrNsgw7cFrWlbuioVlUN6Lxv6I6xCeLyhzf7xXVmS2y7QGArbxEtu8rx4jqjMKfi+nMcVEdGmQrXBxLuEpU18lYK9svgD1Vsr8/13WXdabXO9Uf4Z08eRL33HMPysvLYTAYcNNNN+HWW2/15diIiPyGGUZE3lA9gTKZTLj33ntx+eWXo7a2FpMmTcLQoUNxySWX+HJ8RER+wQwjIm+oPgcqMTERl19+OQAgIiICPXr0QGlpqc8GRkTkT8wwIvKGT04iP3bsGA4cOIArr5SdJ0NEpCXMMCLylNcTqLq6OsyZMweLFi1CRESEL8ZERNRumGFEpIZXE6jm5mbMmTMHmZmZGDVqlK/GRETULphhRKSW6gmUoii477770KNHD9x2222+HBMRkd8xw4jIG6onUF988QXy8vKwc+dOjB8/HuPHj8cnn3ziy7EREfkNM4yIvKG6jcGAAQPw7bff+nIsRETthhlGRN5o307k1RWiugZjiKguvOG0qG52Yqqo7rFaWcff8DNHRHVKkKzbb324RVQHABFlsjEeDJH1sonuIDsI2fm0rONvU9drRHVWYZfvEMg6mwOAuVjW5djRsbOorjRY9nNJ6CB7LP8uknf8nXhFJ3Gtt3Jzc7Ft2zbExcVh06ZNrb6vKAqWLl2KTz75BCEhIXj00Uedl/9fSOzfbJMVhsk6kTdbLhXVzQ3rI6p7ukb2OwoAp6yylQA6ndglqnMk9BDVnQxKENV1bjopqttrk22vb5Q8R5TPN4rqjH1vENUdM8lWmrA6ZDnSPVy29MGROtn2ooLlq0LEKbLVRIJjksTb9IVAZZjqj/CampowefJk/PrXv0ZGRgaeeeYZrwdDRO0vKysLq1evdvv9wsJCHDlyBFu2bMEjjzyChx56qP0G50fMMKLzQ6AyTPVHeMHBwXjttdcQHh6O5uZm3HzzzUhNTUX//v19MjAiah8DBw7EsWPH3H6/oKAAEyZMgMFgQP/+/VFdXY2ysjIkJia24yh9jxlGdH4IVIapPgJlMBgQHh4OALDZbLDZbDAYhKuqEpFulJaWIinp50PySUlJ50XHbmYY0YXBXxmm+ggUANjtdmRlZeHo0aO4+eab2cWXKID+YOju8v7hbz+GNWvWOL/OyclBTk5O+wxK45hhRNrgLr8A7WaYVxMok8mEvLw8VFdX484778ShQ4fQu3dvX42NiDwQbHR99MTbsLFYLCgpKXF+XVJSAotFfuGDljHDiLTBXX4B2s0wn6yFFxUVhUGDBuHTTz/1xeaISIVgo8HlzVtpaWnYsGEDFEXBf//7X0RGRur+/KdfYoYRBZa7/NJyhqk+AlVRUQGz2YyoqCg0Njbis88+w+233+71gIhIHbVBM3/+fOzevRuVlZVITU3F7NmzYbP9eNn3b37zGwwbNgyffPIJ0tPTERoaimXLlvly2AHDDCPSDm8mSoHKMNUTqLKyMtx7772w2+1QFAVjxozB8OHDfTIoIvKc2gB68sknz/l9g8GAxYsXq9q2ljHDiLTDmwlUoDJM9QQqJSUFGzZs8OVYiMgLQbyCzCPMMCLt0GN+eXUSuaesZlln7rC9/5ZtsMfVojJph/E/R8g6/uae/l9RXVKQVVTnycT7eLTsBNdeBlnX61pzlKiuPH+dqM448ypRXWWjXVQHAF3NdaK6o/Gy/j3J1hOiOmkX5lPJg0V1F8fIXv9qmfSXP7pi6ybrsh9UIssbQ3ODqE7aYfyuSPkVhAvL9ovqSrrIXttNdoeornOorOu1LfQiUV1MnazDuGKWd9tuHDJFVBfiaBLVXVRXJqpzhMeJ6gw22WPuHBEuqguC7GcHABu/k9VOihFv0kmP+eX1SeR2ux0TJkzA73//e1+Mh6gF6eSJ/HcS+fmOGUYUeP48idxfvD4C9frrr6Nnz56orZWv80VEvqfloNEyZhhR4Okxv7w6AlVSUoJt27Zh8uTJvhoPEamkt3dvWsAMI9KGC+4I1LJly3D33Xejro4fsxAFmpaDRquYYUTaoMf8Un0E6uOPP0ZsbCyuuOIKX46HiFTS27u3QGOGEWnHBXUE6ssvv8RHH32EwsJCNDU1oba2FgsXLsSKFSt8OT4iEtJy0GgRM4xIO/SYX6onUAsWLMCCBQsAALt27cLLL7/M4CEKID0GUCAxw4i0Q4/51a59oIjIf0xmnyxtSUTU7vSYXz6ZQA0aNAiDBg3yxaaISCWjHjvRaQQzjCiw9JhfBkVRlPba2anqep9uL9peLaoz1leK6o6HJIvqlsfLTjr9w7H/iup6xnQQ1QFAkPA1ZrTKripSjLI5tCMoRFRnqq+Q1dWcEtUBQHOirPt6UOm3ojp7RLxsxybZc9PUIVpUZ7XLf9XiIsPEtT/5oLfrzvyjD33p8baotcYGWefwepvs5xws/INxRti1v9GD19eKxL6iuvvLZasuWKylorr3ymTd+CdFlojq7LFdRXV7KuV/nJuFz+PgBOERE4Os7rQtSFQXK+zmvueE7G/AddZvRHUAgA6y7uami2WrhJzNXX4B2s0wr45ApaWlITw8HEajESaTCevXr/fVuIgAyCdPBJiC5MtV0I+YYUTaoMf88vojvNdeew2xsbG+GAsRecEUrL9zCLSAGUYUeHrML55ETnSeMAXr7x0cERGgz/zyegI1Y8YMGAwG5OTkICcnxxdjIiIV9HgIXAuYYUSBp8f88moC9dZbb8FisaC8vBy33XYbevTogYEDB/pqbETkAaMOD4EHGjOMSBv0mF9ejdhisQAA4uLikJ6ejn379vlkUETkOVOwyeWN3GOGEWmDu/zScoapnkDV19ejtrbW+e/t27ejV69ePhsYEXnGaDS4vJFrzDAi7XCXX1rOMNUf4ZWXl+POO+8EANjtdowbNw6pqak+GxgReUbL79S0iBlGpB16zC/VE6jk5GT861//8uVYiMgLRh0GUCAxw4i0Q4/5xTYGROcJU5D+TsIkIgL0mV9eTaCqq6tx//3349ChQzAYDFi2bBmuuuoqt/WxR3eItmvvcrmoztAkXK4kSLZ8QFKQVVQnXaLlhYv6i+qWVstb6TcKV2sICZK13D9e0yyqiw2RLScR21gjqmva/YGozjRatrwOACjlx0V1Xxpkyz90DJH9eiQFyX4o/l7qSY+HwAPNkwwzVcuWF4k0yn4O/22IFNVdXS1bxqKky2BRHSBfouUvcbJlq5YJMywrQfY7ejRIdi5arFn2O9ot2iGqA4A44cpaQSUHRHU/RF0q3LNsjGV1sry5zlgsqjvT5RpRHQBEH/pYVnixeJNOeswvryZQS5cuxa9+9Ss888wzsFqtaGxs9NW4iMhDxiAeUPYUM4xIG/SYX6qPmdXU1ODzzz/H5MmTAQDBwcGIiory2cCIyDOmYLPLG7nGDCPSDnf5peUMUz2BOnbsGGJjY5Gbm4sJEybgvvvuQ319vS/HRkQeMAaZXd4kCgsLMXr0aKSnp2PVqlWtvr9+/XoMHjwY48ePx/jx47F27VpfD7/dMcOItMNdfkkyLFD5pXoCZbPZ8M033+A3v/kNNmzYgNDQUJcDJ6L2YQwOcnlri91ux5IlS7B69Wrk5+dj06ZNOHz4cKu6G2+8EXl5ecjLy0N2drY/HkK7YoYRaYe7/GorwwKZX6onUElJSUhKSsKVV14JABgzZgy++UZ+MjQR+ZbRaHR5a8u+ffvQrVs3JCcnIzg4GBkZGSgoKGiHEQcWM4xIO9zlV1sZFsj8Uj2BSkhIQFJSEr7//nsAwI4dO9CzZ0+fDYyIPGMMNru8taW0tBRJSUnOry0WC0pLS1vVbdmyBZmZmZgzZw5Onjzp07EHAjOMSDvc5VdbGRbI/PLq7KwHHngACxcuRHNzM5KTk7F8+XKfDIqIPOfuXIE1a9ZgzZo1zq9zcnKQk5Pj0baHDx+OcePGITg4GG+//Tb+/Oc/4/XXX/dqvFrADCPShnOd6+Rthvkrv7yaQPXp0wfr16/3ehBE5D13V6u0FTYWiwUlJT/3OCotLXUusvuTmJgY57+zs7Px+OOPezlabWCGEWnDua62O1eGBTK/9Nf6k4hcUnsFS9++fXHkyBEUFxfDarUiPz8faWlpLWrKysqc//7oo4/4URcR+ZTaq/ACmV+qj0B9//33mDdvnvPr4uJizJkzB9OmTXP7f2zdZB1PjzbIOpJ2jZB1GG80BIvqpIs+94xpuwaQdxi/L+oy2QYBLCjdL6qLCJbNjWNDZM91hLBLrN1oabsIQPDQ8aI646nWV1O483Xn60V1PcJkL3u7sOt7ZaOwS3uof/uZqG1EZzab8eCDD2LmzJmw2+2YNGkSevXqhaeffhpXXHEFRowYgTfeeAMfffQRTCYToqOjz4uPujzNMGNDlWi7zYmyLtoVFbKVFBwJPUR1TXZ5t+1kW1nbRZB3GF8kzLB7Tsnyq4NJll+HymWNTxPD274a9ScGRfaLbw+T/SGIF+ZNdZPs51fVJMubmvjeorpa4X4BIPJi9yuNeEuP+WVQFOGr5RzsdjtSU1PxzjvvoEuXLm7rmmrOiLYnnkCFyX7wvp5ASVmFf4EDOYGSLi8inUCZbLJAM9W0PsnPFekfLQD4Oky2ZILFxxOoRpvsdejJBCo6XPbm4GxnVi1yeX/HWcs83taFRpJh9iOyJZykE6hPimUTqLSOsrpiU7yoDpBPoGrCO4nqAjWBOlEjW37LkwlUUojsF99YI3sOG6M6i+p8PYHqHCHLmyoPJlCdDNWiuuBY2WM+m7v8ArSbYT55S7xjxw4kJyefc/JERP6lx6UQtIIZRhRYeswvn4w4Pz8f48aN88WmiEglScsCco0ZRhRYeswvr08it1qt+OijjzBmzBhfjIeIVDIYjS5vdG7MMKLAc5dfWs4wr6d8hYWFuPzyyxEfL//8nYh8z2CWnetHLTHDiAJPj/nl9QQqPz8fGRkZvhgLEXnBECQ/UZZ+xgwjCjw95pdXx8bq6+vx2WefYdSoUb4aDxGpZQ52fSO3mGFEGuEuvzScYV4dgQoLC8OuXbt8NRYi8oLBrL93cIHGDCPSBj3ml/5Oeycil/QYQEREgD7zy6sJ1Kuvvoq1a9fCYDCgd+/eWL58OTp06OC2vtHo/ntn63nmS1FddZisK2pk2QFR3fFoWefWxFBZU8lGYSNGaXNMAHjC0ldU93SV7Dk0fP2xqM6YnCKqa46Xtci3x3QT1Tk6yvu8xltlDeFMwo6p0YqsKagSJntdN4uqvKDhQ91a5UmGfW7oKtrm1U01oroRjftEdSeCrhPVdRbmEgC893+yRq1ZCcdFddIGmX9NkOXXM5W7RXXx328V1RlTBovqAMAgzBFrx2RRXeihQlFdUGmxqC7h6tGiOsjiC8ebo2SFAA4q4aK6fuItnkWH+aX6HKjS0lK8/vrrePfdd7Fp0ybY7Xbk5+f7cmxE5AFDUJDLG7nGDCPSDnf5peUM8+oIlN1uR2NjI8xmMxobG5GYmOircRGRp4zyIxD0I2YYkUboML9UT6AsFgumT5+O4cOHo0OHDhg6dCiuv162mCsR+Z4ezyEIJGYYkXboMb9Uf4RXVVWFgoICFBQU4NNPP0VDQwPy8vJ8OTYi8oDBHOzyRq4xw4i0w11+aTnDVE+gPvvsM1x00UWIjY1FUFAQRo0aha+++sqXYyMiD+jt/IFAY4YRaccFdQ5U586dsXfvXjQ0NCAkJAQ7duzAFVdc4cuxEZEnNPxOTYuYYUQaosP8Uj2BuvLKKzF69GhMnDgRZrMZffr0QU5Oji/HRkQe0OM5BIHEDCPSDj3ml1dX4c2ZMwdz5szx1ViIyBtG9sX1FDOMSCN0mF/6GzERuaToMICIiAB95pdBURR5q2cvWc+UiepMJd+K6qSH/A6EXyaq6xVSL6qTzpSbgmRdW2uFnW8BINbYJKq7K/pqUd2zJbJO5NVhFlGd3SF7OdU1yx5zsMmr9a5dsjSeENU5jsg6RRt6yJ5rT/qcBCXIul6fzV7suhu0KVnW/ZnO7VhFraiuwSb7Hbi4+Zhsx8K8sXW8SLY9AEElstUZjkb2EtVJf08TDHWiujkx14rqnjv8jqju+6g+ojoA6H5su6ju1MW/EtXFmny7BoGxvlJUVxki62kWFiTP2JAq2WvW3En2ujmbu/wCtJthXv11eu211zBu3DhkZGTg1Vdf9dGQiEgNxWB0eSP3mGFE2uAuv7ScYapHdujQIaxduxZr165FXl4etm3bhh9++MGXYyMiTxjNrm/kEjOMSEPc5ZeGM0z1BOq7775Dv379EBoaCrPZjIEDB2LLli2+HBsReUAxml3eyDVmGJF2uMsvLWeY6glU79698cUXX6CyshINDQ0oLCxESUmJL8dGRJ4wmlzfyCVmGJGGuMsvDWeY6qldz549MXPmTMyYMQOhoaFISUmB0ajdzyqJzndafqemRcwwIu3QY355NeLs7GxkZ2cDAJ588klYLLIrtYjID0z6C6BAY4YRaYQO88urt1vl5eUAgBMnTmDLli3IzMz0yaCISAUvTsAsLCzE6NGjkZ6ejlWrVrX6vtVqxdy5c5Geno7s7GwcOya8BF/jmGFEGuHFSeSByi+vpnyzZ8/GmTNnYDabsXjxYkRFRflkUETkObWHwO12O5YsWYJXXnkFFosFkydPRlpaGi655BJnzdq1axEVFYWtW7ciPz8fK1aswMqVK3019IBhhhFpgx7zy6sJ1D//+U+vB0BEPmIwqPpv+/btQ7du3ZCcnAwAyMjIQEFBQYsA+uijj/CnP/0JADB69GgsWbIEiqLAoHKfWsEMI9IIHeZXu37oqJhkqy37elHB6A6yTyprzbJ3n+GwiuqO18g60MaGyK8yMHwt6xwu7TA+O2m4qO7pqi9FddWGEFFdjPAxVzTaRXUAEBks22Z1RBdRXVSy7Oe8pzFaVHd1tE1Up5bad3ClpaVISkpyfm2xWLBv375WNZ06dQIAmM1mREZGorKyErGxseoHrDPvHy4X1f2uVwdRnePbvaK6/d3HiOpi6uSvr4tiZZ3uY82y19Sh8kZRXfz3W0V10g7jf7zkJlFd6pf/EdUBQNKlw0R1ZmHH+Uaj7PUQ2lwjqrNHyjqMx9Sekm0vKEFUBwCHjbJzBFPEW/yZHvOrzRHn5uZi27ZtiIuLw6ZNmwAAZ86cwbx583D8+HF06dIFK1euRHS07I8IEfmH4uZy3zVr1mDNmjXOr3NycpCTk9Newwo4ZhiR9rnLL0C7GdbmoZmsrCysXr26xX2rVq3CkCFDsGXLFgwZMsTlSVtE1L7sDsXlLScnB+vXr3fefhk8FoulRf+j0tLSVlejWSwWnDx5EgBgs9lQU1ODmJgY/z8oH2CGEWmfu/xqK8MCmV9tTqAGDhzY6p1ZQUEBJkyYAACYMGECPvzwQ68HQkTesSuub23p27cvjhw5guLiYlitVuTn5yMtLa1FTVpaGt577z0AwAcffIDBgwfr5vwnZhiR9rnLr7YyLJD5pepDx/LyciQm/vg5bEJCgvNSYCIKHLtDdk7GL5nNZjz44IOYOXMm7HY7Jk2ahF69euHpp5/GFVdcgREjRmDy5Mm4++67kZ6ejujoaDz11FM+Hn37YoYRaYse88vrk8gNBoNu3okSnc8kR5vcGTZsGIYNa3ny7F133eX8d4cOHfDMM8+o34GGMcOIAk+P+aWqkWZcXBzKysoAAGVlZRfUlThEWmVXFJc3ao0ZRqQt7vJLyxmmagKVlpaGDRs2AAA2bNiAESNG+HRQROQ5u8P1jVpjhhFpi7v80nKGtTmBmj9/PqZMmYKioiKkpqZi7dq1mDVrFrZv345Ro0bhs88+w6xZs9pjrER0DoqiuLxd6JhhRNrnLr+0nGFtngP15JNPurz/tdde8/lgiEg9b84hOJ8xw4i0T4/5ZVDacXrX2NAgqjNVl7RdBMAeldR2EYDgYlkX7fL8daK6iDuWi+pqmmRdtCOEHbQBILj8O1FdVfTFsn0rsg7Cd0VfLaob980uUd2AzhGiOgCIrz8hqjPYZI9FSto5H4rsGHNtdDfxvjtGhIlrf/L9adedjHvER3q8LWqtubRIVFcTJuvW3PHkV6I6m+VSUZ1ilnW8BoAvTsm67HeLlq0sYBNeQdW5uVRU94NJ1m1797EqUV3h1deL6gBgZf0BUV1Q2SFRnT26k3jfEg1Bst/nEEeTqM5UfVK878pIWYYlRPkuvwDtZlibH+Hl5uZiyJAhGDdunPO+zZs3IyMjAykpKdi/f79fB0gXNunkifR3/kB7YYYRad95eQ6Uqy6+vXv3xrPPPouBAwf6bWBE5Bm9XcHSXphhRNqnx6vw2jwHauDAgTh27FiL+3r27Om3ARGROlp+pxZIzDAi7dNjfnndSJOItEHL79SIiM5Fj/nFCRTReaJZj5exEBFBn/nFCRTReaLZocNj4ERE0Gd+cQJFdJ5w6PAQOBERoM/8anMCNX/+fOzevRuVlZVITU3F7Nmz0bFjRzzyyCOoqKjA73//e/Tp0wcvvfRSe4yXiNzQ4yHw9sAMI9I+PeaX6k7k6enpPh8MEanXLGxmeKFhhhFpnx7zq10/wgs++bWozh4WI6orrrGJ6i7qeo2ozjjzKlGdqb5CVBfb6L6zagt1gD1S1r24OV52+bVd2AW92iDrNCztML7pskGiujHHPxDVAUBJSBdRXeeGg6K6r0N7ieq6Rsk6kZUx93IAAB8ASURBVDfZZJ/dR5hUrd0t1qzH64B1xB4p64596JSsI/7V5bIVF5Sjso7XjUOmiOoAoNku61IdJ2xubhB+/GKwyl6j3Y9tF9UlXTpMVJcl7C4OAHPD+ojqnjn1H9kGhSsVmGpPieoqQkNFdXZFtsLFukPyacCsAeJSj+kxv9p85nJzc7Ft2zbExcVh06ZNAIDHHnsMH3/8MYKCgtC1a1csX74cUVFRfh/s+Uo6eboQSSdPpM93cO2BGUakfXrML1WdyIcOHYpNmzZh48aN6N69O1588UW/DZCIZJrtDpe3Cx0zjEj73OWXljOszQnUwIEDER0d3eK+66+/Hmbzjwev+vfvj5IS2aFoIvKfZofi8nahY4YRaZ+7/NJyhnl9DtS7776LsWPH+mIsROQFPV7FogXMMKLA02N+eTWBev7552EymfDrX//aV+MhIpX02Ecl0JhhRNqgx/xSPYFav349tm3bhldffRUGg8GXYyIiFbR8roAWMcOItEOP+aVqAlVYWIjVq1fjH//4B0KFl1QSkX9p+VwBrWGGEWmLHvNLVSfyVatWwWq14rbbbgMAXHnllViyZInfB0tE7ll1+A6uPTDDiLRPj/mlqhN5dna2XwZDROpZhQ09LzTMMCLt02N+cTFhovOEHgOIiAjQZ36p6kS+cuVKFBQUwGg0Ii4uDsuXL4fF0nY37aZ9n8pGNeYOUVmC8Am3Ci+PrGyULX/SsV7Wcr9pt2y5kuCh40V1AGCP6Saqq2uWPTcxIbJ2/wM6R4jqpEu0/KnLaFHdc9++KaoDAFus7LlJsZ0W1W0qChPVDe8e3XYRgO/OyJbPAIC+obIlds7mjwA6c+YM5s2bh+PHj6NLly5YuXJlq55KANCnTx/07t0bANCpUye88MILPh+LWr7KMFNNmWh/V0eHi+ocCWNEddL9hjjkr6/BCbJlhYJKZEugSJffsnZMFtVVhXYW1ZltsmyPrfw/UR0gX6JlTsL1oroZP3wlqkuO7iGqMwn/ntUI/57d0l/2XANAjXApnljxFn/mrwmUPzNMVSfymTNnYuPGjcjLy8MNN9yA//mf/5E+FiLyE6vN4fLmjVWrVmHIkCHYsmULhgwZglWrVrmsCwkJQV5eHvLy8jQ1eQKYYUR64C6/tJxhqjqRR0T8fDSioaGBlwATaUCTzeHy5o2CggJMmDABADBhwgR8+OGHvhhqu2KGEWmfu/zScoapPgfqqaeewoYNGxAZGYnXX3/dZwMiInXsfrgMuLy8HImJiQCAhIQElJeXu6xrampCVlYWzGYzZs2ahZEjR/p8LL7GDCPSDn/kF+DfDFM9gZo3bx7mzZuHF198Ef/4xz8wZ84ctZsiIh9wdxnwmjVrsGbNGufXOTk5yMnJcX49bdo0nD7d+rywuXPntvjaYDC4PVLz8ccfw2KxoLi4GLfeeit69+6Nrl27qnkY7YYZRqQd52pjoNUM8/oqvMzMTMyaNYvhQxRg7s4VyJnSMmx+6dVXX3X7vbi4OJSVlSExMRFlZWWIjXV9euhPJ2AnJyfj2muvxTfffKP5CdRPmGFEgXeuc520mmGySzF+4ciRI85/FxQUoEcP2dUDROQ/Vpvd5c0baWlp2LBhAwBgw4YNGDFiRKuaqqoqWK1WAEBFRQW+/PJLXHLJJV7t19+YYUTa4i6/tJxhqjqRFxYWoqioCAaDAV26dMHDDz/s6WMiIh/zx2XAs2bNwty5c7Fu3Tp07twZK1euBADs378fb7/9NpYuXYrvvvsOixcvhsFggKIouP322zU1gWKGEWmfv9oY+DPD2Imc6Dzh7dUqrsTExOC1115rdX/fvn3Rt29fAMDVV1+NjRs3+nzfvsIMI9I+f+QX4N8MYydyovOEHjv5EhEB+swvVZ3If/Lyyy/jsccew44dO9yemHU2Y/pM2ais9bI6o6xbcwhsorqu5jpRXXNib1GdabSs667x1GFRHQA4Osou9Qw2yU5vqxB2q+3aXCKqKwnpIqqTdhj/46W/FdUBwOwT+0R1KY2y57t/UqKoLhTNorogo397DekxgNqDrzLMESrrOF+ldBDVdWxuENUdM8WL6i6qk3UsBwBHeJyo7oeoS0V18WGy9+KhhwpFdbE9rhXVNRplz7U9upOoDgCgyH6PpB3GX+p2laju6QLZYtbKENnR0y5Vh0R1ZyJTRHUA0Cj7U6qKHvNLVSdyADh58iS2b9+Ozp3lbeCJyH8Uh+LydqFjhhFpn7v80nKGqepEDgDLly/H3XffzQ6+RBphtztc3i50zDAi7XOXX1rOMFXnQH344YdITExESor80B8R+Zddh4fAA4UZRqQteswvjydQDQ0NePHFF/Hyyy/7YzxEpJJDuEr7hY4ZRqQ9eswvjxtpHj16FMeOHcP48eORlpaGkpISZGVl4dSpU/4YHxEJ2W0OlzdqiRlGpD3u8kvLGebxEahLL70UO3bscH6dlpaGdevWia7CIyL/0fK5AlrCDCPSHj3mV5tHoObPn48pU6agqKgIqampWLt2bXuMi4g85LA5XN4udMwwIu1zl19azjBVncjP9tFHH/lsMESknh7fwbUHZhiR9ukxv9iJnOg8Ybfp7yRMIiJAn/mlqhP5s88+i3feecd5zsD8+fMxbNiwNndmaKgSDeqoIuv4e1GErH+L+Yf/yvYb319U17X0W1GdUn5cVPd15+tFdQAQb/XtLD0y2CSqMzQ0iuo6NxwU1dliu4nqpN3FAeDZzv1EdUvOfC2qiw+SXWMRdGKvqO6SsBhR3Y+iPKj9kZYbzgWSrzLMERwm2l9jvaxds6lGlg/W0ItFddLu4gBw2hYkrJTlTXWTrC6otFhUZxB2Ig9trhHVecJUK7uYIDm6h6hO2mH8rhEPiuoWV04S1SUaZdkeAauoDgAi7bXCynDxNn+ix/xqcwKVlZWFqVOn4s9//nOL+6dNm4YZM2b4bWBE5Bk9HgJvD8wwIu3TY361OYEaOHAgjh071h5jISIvaPlky0BihhFpnx7zy+M+UD958803kZmZidzcXFRVyT6aIyL/0dsyCIHGDCPSDj0u5aJqAvWb3/wGW7duRV5eHhITE/Hoo4/6elxE5CG9NaELJGYYkbbosZGmqglUfHw8TCYTjEYjsrOzsX//fl+Pi4g85LBZXd6oNWYYkba4yy8tZ5iqNgZlZWVITEwE8OOinL169fLpoIjIc45m7QaN1jDDiLRFj/nV5gRq/vz52L17NyorK5GamorZs2dj9+7dOHjwx8vVu3TpgiVLZJdpEpH/aPmdWiAxw4i0T4/5paoTeXZ2tl8GQ0Tq6TGA2gMzjEj79Jhf7EROdJ5QHPZAD4GISBU95le7TqCOQ9ZhvKtJ1l3Wqsg6O5s6dhbVJVtPiOrsEfGiui8NXUV1PcLkPwaTUdZ9PbpW1uW4OqKLeN8SX4fKziVJsZ2W1TUeFu9b2mH8wY6Xi+qWVn8jqutwSvZc114h65QOALHiyp/ZdfgOTk8qGmQB36Vkj2yDwk7R3eNlv/MGm6wDOgDEhoaI6srqZN2hq5pkz03C1aNFdYb6SlGdPTJRVNekyK+XqggNFdWZ7LLnRhkiO9op7TD+cIwsv3773ReiuuuOyy+gOBh9pahONsKW9Jhfbb6qcnNzMWTIEIwbN67F/W+88QbGjBmDjIwM/PWvf/XbAIlIxtFsdXm70DHDiLTPXX5pOcNULeWyc+dOFBQU4F//+heCg4NRXl7u10ESUdv0eA5Be2CGEWmfHvNL1VIub731FmbNmoXg4GAAQFycfBFLIvIPh6050EPQJGYYkfbpMb9UNdI8cuQI9uzZg+zsbEydOhX79u3z9biIyEP+aEK3efNmZGRkICUl5ZzNJgsLCzF69Gikp6dj1apVXu2zPTDDiLTFX400/ZlhqiZQdrsdVVVVeOedd3DPPfdg7ty5UBTZCXVE5B/2ZqvLmzd69+6NZ599FgMHDnS/X7sdS5YswerVq5Gfn49Nmzbh8GH5yf+BwAwj0hZ3+aXlDFN1FZ7FYkF6ejoMBgP69esHo9GIyspKxMaquXaIiHzBH+cQ9OzZs82affv2oVu3bkhOTgYAZGRkoKCgAJdcconPx+MrzDAibfHXOVD+zDBVR6BGjhyJXbt2AQCKiorQ3NyMmBhZSwEi8g/FYXd587fS0lIkJSU5v7ZYLCgtLfX7fr3BDCPSFnf5peUMU7WUy6RJk7Bo0SKMGzcOQUFBePTRR2EwyHqVEJF/NH3xd5f3r1mzBmvWrHF+nZOTg5ycHOfX06ZNw+nTrftyzZ07FyNHjvT9QNsZM4xI+9zlF6DdDFO1lAsArFixwueDISLf+2XY/NKrr77q1fYtFgtKSkqcX5eWlsJisXi1TV9ihhHpm1YzrF07kZ+ul3XKvShE9llok7AT7JlgWZh3OrFLVIcusj6rHUNkT6/wYQAAopVGUZ3jiOyqoqhk2XOtmIJFdV2jZHWbisJEdf2TZJ2GASA+SPaJtLTD+H1Rl4nqnjkj6zzdsex/RXUAgMhr5bUB1rdvXxw5cgTFxcWwWCzIz8/HE088Eehh+VxMiKxzuK3bNaK6kkbZEa/GOoeornNEuKgOAP57ok5Ud52xWFRXE99btmNZfKEyRPZ7H1N7SlQXEtpRtmMAdkX2c65pFHamrzokqksUdqaXdhh/s6fsddipTN6JvJPwd0Bv1GZYm3/hc3NzsW3bNsTFxWHTpk0AfjwsVlRUBACoqalBZGQk8vLyvHwIRKQ1W7duxSOPPIKKigr8/ve/R58+ffDSSy+htLQU999/P/7+97/DbDbjwQcfxMyZM2G32zFp0iT06iVb0qc9MMOILlz+zDBVnchXrlzp/Pejjz6KiIgIlQ+NiLQsPT0d6enpre63WCz4+99/Pmdh2LBhGDZsWHsOTYwZRnTh8meGtfmZx8CBAxEd7XoRYEVRsHnz5lZrTBERaQUzjIj8QVUbg5/s2bMHcXFx6N69u4+GQ0TUfphhRKSWVxOoTZs28Z0bEekWM4yI1FI9gbLZbNi6dStuvPFGX46HiKhdMMOIyBuqJ1CfffYZevTo0aJ7JxGRXjDDiMgbbU6g5s+fjylTpqCoqAipqalYu3YtAODf//43MjIy/D5AIiJvMMOIyB9UdyJ/9NFHfT4YIiJfY4YRkT8YFEXxoA+2d6wVJ0R1/9ck68nSJVLW6TvUIOsYW2mTfaIZLux4bXfIntpKYUdbAOgUJtu3qaZMVPdFo+vLu39pgFm2OGxFRFdRXbBJvu5YKJpFdUElB0R1tlPHZTu+TNYTZE7HAaK6pxoOyvYLIDQkRFxL7ePkGVn37tggWedwg1W2vXJDpKguRrYIAAAg6OiXorozXWTdrGutsscsresWHSSqM0O2PXPlUVEdAKw4JPu7ckv/zqK6ELMs6yIgWxXCfFzWObworr+obkViX1EdAKysl2VsSGioeJt61uZf49zcXAwZMqTFlSoHDhzATTfdhPHjxyMrKwv79smWDSHylHTyROQOM4yI/KHNCVRWVhZWr17d4r7HH38cd955J/Ly8nDXXXfh8ccf99sAiYi8wQwjIn9Q1YncYDCgru7Hw881NTVITJQv+EpE1J6YYUTkD7IPe39h0aJFmDFjBh577DE4HA68/fbbvh4XEZHfMMOIyFuq+kC99dZbyM3NxSeffILc3Fzcd999vh4XEZHfMMOIyFuqJlDvvfceRo0aBQAYO3YsT8AkIl1hhhGRt1RNoBITE7F7924AwM6dO7kQJxHpCjOMiLzV5jlQ8+fPx+7du1FZWYnU1FTMnj0bjzzyCJYtWwabzYYOHTpgyZIl7TFWIiKPMcOIyB9UdyJfv369zwdDRORrzDAi8od27UTeWC/rvPt9tawzt7TD696SWlHdxTGy7qnJUbKWvx4024bBICs2CrcZUndKVKcEyTpe15tl3eGDTbJPhb870ySqA4Ag4YO+xH5SVFcdmSyq61j2v6K6pi79RHXzQlNEdQDwgnJEXEvtw/7NNlFdffdBorr9ZfWiugHRsmayG4/JunIDwISIElGdUiH7nXJcfJWo7mBjuKjuMoNsJYXDRouoLi5UfsG5NLdrhF3VzcL8SrBXiur+rzlKVNcpQvaYw4R/RwFgblgfUd2Fkl+qOpEfPHgQOTk5yMzMxB/+8AfU1somKOSadPJ0IZJOnojcYYYRkT+o6kR+3333YcGCBdi4cSNGjhzZ6vtERFrBDCMif1DVifzIkSMYOHAgAGDo0KHYsmWLf0ZHROQlZhgR+YOqNga9evVCQUEBAOD999/HyZOyz8mJiLSAGUZE3lI1gVq6dCn++c9/IisrC3V1dQgOlp1UTUSkBcwwIvKWqrXwevbsiZdffhkAUFRUhG3btvlyTEREfsUMIyJvqToCVV5eDgBwOBx4/vnnMWXKFJ8OiojIn5hhROQtVZ3I6+vr8c9//hMAkJ6ejkmTJvl9oEREajDDiMgfVHciv/XWW30+GCIiX2OGEZE/qPoIj4iIiOhC1q5LuRARERGdD3gEioiIiMhDnEAREREReYgTKCIiIiIPcQJFRERE5CFOoIiIiIg8xAkUERERkYc4gSIiIiLykOmhhx56qD13+N1332HdunXYvHkzCgsLcfDgQXTs2BGxsbGqt/ftt98iNja2xYrqhYWF6Natm/Prffv2obS0FBaLBYcPH0ZeXh6qq6vRvXv3c27/nnvuQXp6epvj2LNnDz744APU1dWha9euzvv37t2LiIgIBAcHo7GxEc899xxeeeUVHDhwAP369UOHDh0AAK+//joSExMRGRl5zv1YrVZs3LgRFRUVSE5OxsaNG7FmzRocO3YMffr0gclkctYWFxdj3bp1+Pe//43t27fj2LFjuPjii7nyPJFKgcovQF2GMb+I/KddG2muWrUK+fn5yMjIgMViAQCUlpY675s1a1ab23j33Xed61a9/vrrePPNN9GzZ08cPHgQixYtwsiRIwEAEydOxHvvvQcA+Nvf/obCwkLYbDYMHToUe/fuxaBBg/DZZ5/h+uuvxx133AEA+MMf/tBqf7t27cKgQYMAAC+88ILz/smTJ2PdunUAgHfeeQdvvvkm0tPT8Z///AdpaWnOx5KRkYG8vDyYzWY88MADCAkJwejRo7Fz504cPHgQf/vb3wAA11xzDUJDQ9G1a1dkZGRg7NixLkN5wYIFsNvtaGxsRGRkJOrr65Geno6dO3dCURQ89thjzudm27ZtGDBgAAoLC9GnTx9ERUVh69atWLx4sfMxEZGML/IL+DnDpPkFyDKM+UXUzpR2NGrUKMVqtba6v6mpSUlPTxdtY9iwYc5/jxs3TqmtrVUURVGKi4uViRMnKq+++qqiKIoyfvz4FnU2m02pr69XrrrqKqWmpkZRFEVpaGhQxo0b56ybMGGCsmDBAmXnzp3Krl27lJ07dypDhw5Vdu3apezatavFOM7eflZWllJeXq4oiqLU1dW12OaYMWNabP9sv/71r1tsz263K59++qmSm5urDBo0SJk+fbqyfv1653h/eiyKoijNzc3KkCFDFJvNpiiKojgcjhb7/ekxK4qi1NfXK1OnTlUURVGOHz/eYux6dfr0aZ9ur6Kiwqfbo/OPL/JLUX7OMGl+/VTbVoYxv/SFGaZ/7XoOlMFgQFlZWav7T506BYPB4Pw6MzPT7e306dPOOofDgfDwcADARRddhDfeeAOFhYVYvnw5lLMOrJlMJphMJuc7pIiICABASEgIjMafn4J3330XV1xxBV544QVERkZi0KBB6NChA6699lpce+21LcbscDhQVVWFyspKKIrifLcVFhbW4jB0r1698O677wIAUlJSsH//fgBAUVERzOaf13I2GAwwGo24/vrrsWzZMnz66ae4+eab8emnnzrflQKAoiiwWq2oq6tDQ0MDampqAPx4aNxms7UYo91ud36vrq4OANC5c+dWdTU1NVixYgXGjBmDa6+9FoMGDcLYsWOxYsUKVFdXt/p5uTJz5kznv2tra/HEE0/g7rvvxsaNG1vUnf2J8alTp7B48WI8/PDDqKysxLPPPovMzEzcddddLV4nZ86caXGrrKxEdnY2qqqqcObMGWddYWFhi8e0aNEiZGZmYsGCBS1eNytWrEBFRQUAYP/+/RgxYgRuuukmDB8+HLt3724x3okTJ+K5557D0aNHz/n49+/fj1tuuQULFy7EyZMncdttt+Gaa67BpEmT8M033zjr6urq8PTTTyMjIwPXXHMNBg8ejJtuugnr168/5/Yp8KT5BcgyTJpfgCzDmF/azC8gcBnG/PIvc9slvrNo0SJMmzYN3bp1Q6dOnQAAJ06cwNGjR/HAAw8468rLy/HSSy8hKiqqxf9XFAVTpkxxfh0XF4cDBw6gT58+AIDw8HC8+OKLWLRoEQ4dOuSsCwoKQkNDA0JDQ1v8oGtqalpMoIxGI6ZNm4YxY8Zg2bJliI+Pd/4S/1JtbS2ysrKgKIozWBMTE1FXV9ci/JYuXYqlS5fi+eefR0xMDKZMmYKkpCR06tQJS5cubfHYzhYUFIQRI0ZgxIgRaGhocN4/efJkjB07Fg6HA/PmzcNdd92F5ORk7N27FxkZGS3qJk2ahCuvvBJ79uzB7bffDgCoqKhAdHR0i33NnTsXgwYNwhtvvIGEhAQAP4bDe++9h7lz5+Lll18GAHz99dcunwtFUXDw4EHn17m5uejWrRtGjx6NdevWYcuWLXjiiScQHByMvXv3Ouvuvfde3HDDDWhoaMDvfvc7ZGZmYtWqVfjwww+xePFiPP/88wCAwYMHo3Pnzi32WVpaiokTJ8JgMKCgoAAA8NRTTyE1NRUA8OijjyIhIQEvvPACtm7digcffBDPPfccAOCTTz7BwoULAQB//etf8dRTT6Ffv34oKirCggULWrxGqqqqUFNTg9/97neIj4/HuHHjMHbsWOdHOD95+OGHMXv2bNTU1GDKlCnIzc3FK6+8gh07duDhhx/GmjVrAAALFy5Eeno6XnrpJWzevBn19fXIyMjA888/jyNHjmD+/Pkun2MKPGl+AbIMk+YXIMsw5pc28wsIXIYxv/ysvQ952e125auvvlLef/995f3331e++uor56Han+Tm5iqff/65y/8/f/58579PnjyplJWVuazbs2eP899NTU0ua8rLy5WDBw+6HevHH3+sPPHEE26/70p9fb1y9OjRVvfX1NQoBw4cUPbv36+cOnWq1fe///578T5KSkqUkpISRVEUpaqqStm8ebOyd+/eVnWHDh1SNm/erBw+fPic2xs1apToeykpKcott9yiTJ06tdWtb9++zrqzD+0riqI899xzSk5OjlJRUdHiY4CzD8Wf/dHsL7fx0ksvKdOnT2/xsxo+fHirsZ697V+O4eyvx4wZozQ3NyuKoijZ2dkt6s7+GOGX2/z888+VxYsXK9ddd50ydepU5e233xY9lrO/l5mZ2eJ7WVlZiqL8+HsxevToVo+JtEWSX4oiyzBpfimKugxjfmkjvxQlcBnG/PKvdj0CBfz4Lql///7nrFm2bJnb7z3xxBPOfyclJbmtu+aaa5z/dnfVRmxs7Dmvnrnhhhtwww03nGOkrYWGhiI5ObnV/REREUhJSXH7/y6++GLxPs5+9xAVFYUxY8a4rOvVqxd69erV5va6dOmCv//975g4cSLi4+MBAKdPn8b69eud77QBoGfPnliyZInLq36GDRvm/LfVaoXD4XC+M77jjjtgsVgwdepU1NfXO+scDofz3+PHj2+xvbO/N336dNx4441YtmwZOnXqhNmzZ7f6yAT48V3/K6+8AkVRUFtb63x3/cvt3XzzzZg1axZuv/12/OpXv8Jf/vIXjBo1Cjt37jznz2jAgAEYMGAAHnjgAWzfvh2bN29GTk4OAKBDhw74z3/+g5qaGhgMBnz44YcYOXIkdu/e3eIoZ1hYGPbs2YMBAwagoKAAHTt2BPDj74XSftdzkEqS/AJkGSbNL0BdhjG/tJFfgDYyjPnle+0+gSLteeqpp7Bq1SpMnTrV+bl6XFwc0tLS8PTTTzvr/vSnP7UKhp+c/RHG8OHDsXPnTlx33XXO+7KyshAfH4+//OUvzvtGjBiBuro6hIeHY968ec77f/jhh1aBnJSUhGeeeQYFBQWYPn06GhsbW43hpptucp4rMXHiRFRWViI2NhanTp1yfkwCALfccgt69+6Nt956C0eOHIHdbscPP/yAkSNH4o9//GOLbboKW5PJhNTUVOehduDHQ+CPP/44DAYDVq9ejbfeegv33nsvLBYLHnnkEWfdQw89hPvvvx8//PADLrnkEucf2oqKCvz2t791+dwSkXt6yC8gMBnG/PKzQB7+Iu1bt26d5uoaGhqUb7/9tt336+9terJvImqbFvNLUbSRYcwv73ECRef0y8/DWaeNfRNR25gj2qs7n/AjPEJmZqbb75196eyFVhfofRNR25gj2qu7UHACReK2ERdaXaD3TURtY45or+5CwQkU4YYbbkBdXV2LkxR/cvaSCRdaXaD3TURtY45or+5C0a5r4RERERGdD9p1KRciIiKi8wEnUEREREQe4gSKiIiIyEOcQBERERF5iBMoIiIiIg/9PyCWxBJSCZQhAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 720x360 with 4 Axes>"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def upper(df):\n",
        "    '''Returns the upper triangle of a correlation matrix.\n",
        "\n",
        "    You can use scipy.spatial.distance.squareform to recreate matrix from upper triangle.\n",
        "\n",
        "    Args:\n",
        "      df: pandas or numpy correlation matrix\n",
        "\n",
        "    Returns:\n",
        "      list of values from upper triangle\n",
        "    '''\n",
        "    try:\n",
        "        assert(type(df)==np.ndarray)\n",
        "    except:\n",
        "        if type(df)==pd.DataFrame:\n",
        "            df = df.values\n",
        "        else:\n",
        "            raise TypeError('Must be np.ndarray or pd.DataFrame')\n",
        "    mask = np.triu_indices(df.shape[0], k=1)\n",
        "    return df[mask]\n",
        "\n",
        "# Compare that the upper triangle is properly extracted. \n",
        "print(\"Extracted upper triangle from m1\")\n",
        "print(np.round(upper(m1)[:25],3))\n",
        "print(\"Top 3 rows from m1\")\n",
        "display(m1.iloc[:3, :].round(3))\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 275
        },
        "id": "iLfFHclR-FNI",
        "outputId": "3a81f1a6-fae6-4551-c407-b7ebf84689d9"
      },
      "execution_count": 77,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Extracted upper triangle from m1\n",
            "[ 0.002  0.023  0.04  -0.024 -0.063  0.228  0.036 -0.035  0.056  0.034\n",
            "  0.003  0.181  0.367  0.022  0.079 -0.126 -0.222  0.075  0.011  0.235\n",
            " -0.042  0.276 -0.065  0.069 -0.153]\n",
            "Top 3 rows from m1\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "\n",
              "  <div id=\"df-74026d0e-f275-40ed-a748-c24fde4fba9b\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>0</th>\n",
              "      <th>1</th>\n",
              "      <th>2</th>\n",
              "      <th>3</th>\n",
              "      <th>4</th>\n",
              "      <th>5</th>\n",
              "      <th>6</th>\n",
              "      <th>7</th>\n",
              "      <th>8</th>\n",
              "      <th>9</th>\n",
              "      <th>10</th>\n",
              "      <th>11</th>\n",
              "      <th>12</th>\n",
              "      <th>13</th>\n",
              "      <th>14</th>\n",
              "      <th>15</th>\n",
              "      <th>16</th>\n",
              "      <th>17</th>\n",
              "      <th>18</th>\n",
              "      <th>19</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>1.000</td>\n",
              "      <td>0.002</td>\n",
              "      <td>0.023</td>\n",
              "      <td>0.040</td>\n",
              "      <td>-0.024</td>\n",
              "      <td>-0.063</td>\n",
              "      <td>0.228</td>\n",
              "      <td>0.036</td>\n",
              "      <td>-0.035</td>\n",
              "      <td>0.056</td>\n",
              "      <td>0.034</td>\n",
              "      <td>0.003</td>\n",
              "      <td>0.181</td>\n",
              "      <td>0.367</td>\n",
              "      <td>0.022</td>\n",
              "      <td>0.079</td>\n",
              "      <td>-0.126</td>\n",
              "      <td>-0.222</td>\n",
              "      <td>0.075</td>\n",
              "      <td>0.011</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.002</td>\n",
              "      <td>1.000</td>\n",
              "      <td>0.235</td>\n",
              "      <td>-0.042</td>\n",
              "      <td>0.276</td>\n",
              "      <td>-0.065</td>\n",
              "      <td>0.069</td>\n",
              "      <td>-0.153</td>\n",
              "      <td>0.203</td>\n",
              "      <td>0.216</td>\n",
              "      <td>0.196</td>\n",
              "      <td>0.105</td>\n",
              "      <td>0.117</td>\n",
              "      <td>-0.007</td>\n",
              "      <td>0.043</td>\n",
              "      <td>-0.146</td>\n",
              "      <td>0.074</td>\n",
              "      <td>0.098</td>\n",
              "      <td>-0.184</td>\n",
              "      <td>-0.179</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>0.023</td>\n",
              "      <td>0.235</td>\n",
              "      <td>1.000</td>\n",
              "      <td>0.090</td>\n",
              "      <td>0.104</td>\n",
              "      <td>0.067</td>\n",
              "      <td>0.113</td>\n",
              "      <td>0.031</td>\n",
              "      <td>-0.047</td>\n",
              "      <td>0.197</td>\n",
              "      <td>0.198</td>\n",
              "      <td>0.002</td>\n",
              "      <td>-0.151</td>\n",
              "      <td>-0.045</td>\n",
              "      <td>-0.043</td>\n",
              "      <td>0.001</td>\n",
              "      <td>-0.022</td>\n",
              "      <td>-0.062</td>\n",
              "      <td>-0.041</td>\n",
              "      <td>-0.138</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-74026d0e-f275-40ed-a748-c24fde4fba9b')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-74026d0e-f275-40ed-a748-c24fde4fba9b button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-74026d0e-f275-40ed-a748-c24fde4fba9b');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ],
            "text/plain": [
              "      0      1      2      3      4   ...     15     16     17     18     19\n",
              "0  1.000  0.002  0.023  0.040 -0.024  ...  0.079 -0.126 -0.222  0.075  0.011\n",
              "1  0.002  1.000  0.235 -0.042  0.276  ... -0.146  0.074  0.098 -0.184 -0.179\n",
              "2  0.023  0.235  1.000  0.090  0.104  ...  0.001 -0.022 -0.062 -0.041 -0.138\n",
              "\n",
              "[3 rows x 20 columns]"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# You can go back to the original matrix using scipy's squareform. \n",
        "# But note that the diagonal will be 0s. \n",
        "recovered_m1 = pd.DataFrame(squareform(upper(m1))).round(2)\n",
        "\n",
        "f,axes = plt.subplots(1,2, figsize=(10,5))\n",
        "sns.set_style(\"white\")\n",
        "for ix, m in enumerate([m1,recovered_m1]):\n",
        "  sns.heatmap(m, cmap=\"RdBu_r\", center=0, vmin=-1, vmax=1, ax=axes[ix], square=True, cbar_kws={\"shrink\": .5})\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        },
        "id": "WVUh2IqdC61J",
        "outputId": "1df17711-8a76-479e-f301-da45fe52c862"
      },
      "execution_count": 73,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAADqCAYAAABtLHt/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3RU9bk38O9cMiTkSgKZCASQAILIpdYc5cUTSiBcDAgNxGjrWYpSaF9LhKi0gXP0HFyAqFWq73ktlHJAWxFRTGqip2gwxCU3bzUq0FZN5CKZQAi5kWQu2e8fHvOWMkOe2XvPzN7k+1lr1iLJ42+eJJOvv9mz59kWRVEUEBEREZGYNdINEBEREZkNN1BEREREQeIGioiIiChI3EARERERBYkbKCIiIqIgcQNFREREFCS7lv+4qqoKa9euRVdXF/Lz87FkyZLL1v/UMky07sYLR0V1li6vqA4+j6jMbY8R1XV4ZZMfko8fENV5h35fVAcAHdY+oroYT4uoTrE5ZHV22f06Tn8uquusfldUZ81ZLKoDAEt7k6juFBJFdWcvyB5fE+LaRXVd0bL7BYDovrHiWlIvmAzraJf9nqX0zi+PML8AoNMny7Ckr4UZNuR7orp2a7Sorq+3VVQnzSVpzgHyDOv4uFJUZ5u1VFQXsfyK7xTVAUBXn3hRXW/JL9VHoHw+H9asWYMtW7agvLwcZWVl+OKLL/TsjYgoZJhhRKSF6g1UdXU1hg4divT0dDgcDuTm5qKiokLP3oiIQoYZRkRaqN5AuVwupKWldX/sdDrhcrl0aYqIKNSYYUSkBU8iJyIiIgqS6g2U0+lEXV1d98culwtOp1OXpoiIQo0ZRkRaqN5AjRs3DrW1tThx4gTcbjfKy8uRnZ2tZ29ERCHDDCMiLVSPMbDb7Xj44YexePFi+Hw+LFiwACNHjtSzNyKikGGGEZEWmuZATZkyBVOmTNGrFyKisGKGEZFaPImciIiIKEgWRVFkI2n9KC4uRmVlJVJSUlBWVtZjvXSS7/K+Y0R1/+fPm0V1rRk3i+r6fvqmqK5p7GxRXaKvWVT3tUc+tTXjfLWoTunqEtVZ7FGiOl+88ORaRXa/3qRBojqr+4LsfgGc6JR9L0NssintFq9bVPdXX7JsPYuoDABwbVqCvJhUCVV+SUV9uV9U1zJskqhOml8A0HKdLMPivfpmWEazbMq34u4Q1VkcssnmvoS0nou+W1N65YqkdFGdzSPLsOMdsvwaGtUmqrN4ZD9DaX4FY+xVvSO/NB2BysvLw5YtW/TqhYgobJhfRKSFpg1UZmYmEhPl1/ciIjIK5hcRacFzoIiIiIiCxA0UERERUZC4gSIiIiIKEjdQREREREHStIEqKirC7bffjpqaGmRlZWHXrl169UVEFFLMLyLSQtMk8qeeekqvPoiIwor5RURa8CU8IiIioiBpmkQerM422QRo2xcHRHU/n7hEVPdsfZWoDl0+UZniiBHVWTplE2O74vqL6gCgHbJptQl1sonlUu5BE0R1J1q8oroBfW1a2vEr2iYb9e0VPuI7fbJC4d3iXIfs8QUAw/vHi2spPKT5Za85LKrzZPwvUZ2t9YyoTnoVAABQomQTvPXOsDZFll9JrsjkFxC5DDN6fgHyDOst+aX6JbzTp09j5cqVaGhogMViwW233Ya77rpLz96IiEKGGUZEWqjeQNlsNvzyl7/E2LFj0draigULFmDy5MkYMWKEnv0REYUEM4yItFB9DlRqairGjh0LAIiLi8Pw4cPhcrl0a4yIKJSYYUSkhS4nkZ88eRJHjx7FhAny15mJiIyCGUZEwdK8gWpra0NhYSFWrVqFuLg4PXoiIgobZhgRqaFpA+XxeFBYWIi5c+dixowZevVERBQWzDAiUkv1BkpRFKxevRrDhw/HokWL9OyJiCjkmGFEpIXqDdSHH36I0tJSHDx4EPPmzcO8efOwb98+PXsjIgoZZhgRaaF6jMENN9yAv/zlL3r2QkQUNswwItIivJPIm8+J6tqtsim5se1nRXXLUrNEdRtaj8ru93ytqE6Jkk0svxDrFNUBQFy9rMdj0bJZNol9ZAchB579RFTXOeT7ojq3cEpuNGRTgQHAfuLPorqupIGiOpdD9nsZ0Ef2vbxR0yqqA4AfXneVuFar4uJiVFZWIiUlBWVlZZd8XVEUrF27Fvv27UN0dDQee+yx7rf/9yaRyi9f3ABRnSeIJI9prBXVSTOsrW+qqC6h/oio7mjMKFFdvDC/Bp/7TFQHyKeWS3/eDkWWYfavPxTVKf0Giepc0bKc6++QT7CXZlg48wuIXIapfgmvs7MTCxcuxK233orc3Fw888wzmpshovDLy8vDli1bAn69qqoKtbW12LNnDx599FH8+7//e/iaCyFmGNGVIVIZpvolPIfDge3btyM2NhYejwc/+tGPkJWVhYkTJ+rSGBGFR2ZmJk6ePBnw6xUVFZg/fz4sFgsmTpyI5uZm1NfXIzVVdtTBqJhhRFeGSGWY6g2UxWJBbGwsAMDr9cLr9cJiCeKqhESkq59ahvn9/NSXNmDnzp3dHxcUFKCgoEC8rsvlQlpaWvfHaWlpcLlcpt9AMcOIjCNQfgHGzTDVGygA8Pl8yMvLw/Hjx/GjH/2IU3yJIijQVdWDDZvehBlGZAyB8gswboZpGqRps9lQWlqKffv2obq6Gn/961/16ouIguSwWvzetHI6nairq+v+uK6uDk6n/I0PRsYMIzKGQPll5AzT5Vp4CQkJuPHGG/Huu+/qsRwRqWCzWPzetMrOzkZJSQkURcGf//xnxMfHm/7lu3/EDCOKrED5ZeQMU/0S3rlz52C325GQkICOjg7s378fP/nJTzQ3RETqqH2mVlRUhMOHD6OxsRFZWVlYtmwZvN5v33p9xx13YMqUKdi3bx9ycnIQExODdevW6dl2xDDDiIxDy5GmSGWY6g1UfX09fvnLX8Ln80FRFMyaNQtTp07VpSkiCp7aAHrqqacu+3WLxYJHHnlE1dpGxgwjMg4tG6hIZZjqDdTo0aNRUlKiZy9EpEEU30EWFGYYkXGYMb80vQsvWG67bKpt30/ekC04/HpRmXTC+C/ixojqis/KptqmRblFdcFsvE8lyib0jrTIJsa22hNEdQ3lr4jqrIu/J6pr7PCJ6obY20R1AHC8v2x+T7r7G1HdVd8cEtWdSb9JVHd1P9njXy09TrakwDx651eGbGq/dOJ1VDA5Ep0uqpNmmE342DuZNFpUp3d+nS3ZIaoDANvS8aK6pk5Zhg22yTLseKrs8SDNr7RTB0V10vwCQpthZswvzSeR+3w+zJ8/H0uXLtWjHyJSKVTvYLnSMcOIIi+U78ILFc1HoJ5//nlkZGSgtVV+nS8i0p+Rg8bImGFEkWfG/NJ0BKqurg6VlZVYuHChXv0QkUoOq/8bBcYMIzKGQPll5AzTdARq3bp1eOihh9DWJj9PhYhCQ495Kb0NM4zIGMyYX6r3du+88w6Sk5Nx3XXX6dkPEalktvMHIo0ZRmQcveocqI8++gh79+5FVVUVOjs70draigcffBBPPvmknv0RkZCRg8aImGFExmHG/FK9gXrggQfwwAMPAAAOHTqErVu3MniIIsiMARRJzDAi4zBjfoV1DhQRhY7NbuCzLYmILsOM+aXLBurGG2/EjTfeqMdSRKSSzWGLdAumxQwjiiwz5pdFURThnFvtzjRf0HW9RF+zqM56oVFUJ53Ou76/7KTTn578s6guo18fUR0gnzZsdcveVaRYZXvorqhoUZ3twjlZXcsZUZ0nVTZ5HQCiXH8R1fni+ssWtMl+Np19EkV1bp/8Ty0lvq+49jtvjfE/yTjn6IdBr0WXamjRN7/ivZHJLwBw9pX9z+rYuU5RnTTDotAlqrN62kV1ii1KVOezyTPW3nZWVGdraxDVSTPM6PkFyDNMz/wCjJthmo5AZWdnIzY2FlarFTabDbt379arLyIKki3KfM/gIo0ZRmQMZswvzS/hbd++HcnJyXr0QkQaWE14CNwImGFEkWfG/OJJ5ERXCKvNfO9iISICzJlfmjdQ9957LywWCwoKClBQUKBHT0SkghkPgRsBM4wo8syYX5o2UDt27IDT6URDQwMWLVqE4cOHIzMzU6/eiCgIViNfNMqgmGFExmDG/NLUsdPpBACkpKQgJycH1dXVujRFRMGzOWx+bxQYM4zIGALll5EzTPUG6sKFC2htbe3+93vvvYeRI0fq1hgRBccWZfN7I/+YYUTGESi/jJxhql/Ca2howH333QcA8Pl8mDNnDrKysnRrjIiCYzPhIfBIYoYRGYcZ80v1Bio9PR1//OMf9eyFiDQw8qFuI2KGERmHGfMrrGMMko8fENX5Bo0V1Vk6hdO2o2JEdWlRblGddML4bwZPFNWtbT4iqgOADuEw6+ioWFHdqRaPqC452ier62gR1XUe/pOoDvgTbDMXiyqVhlOiuo8sQ0R1SdGyP4+0KNkvJdTv0jXyoe4rQdLXV0Z+AcDRc7Jn+6OTZRO827yyvwE3ZH8EDp3zKyVGNgEdAPoJr+IgzTDbTNmEeKPnFwDhb08dM+aXpmNmzc3NKCwsxKxZszB79mx8/PHHevVFBADizRN9O4jO340CY4YRGUOg/DJyhmk6ArV27Vr88z//M5555hm43W50dHTo1RcRBcliNd85BJHGDCMyBjPml+qOW1pa8P7772PhwoUAAIfDgYSEBN0aI6Lg2Bx2vzfyjxlGZByB8svIGaZ6A3Xy5EkkJyejuLgY8+fPx+rVq3Hhgr5XKyciOWuU3e9NoqqqCjNnzkROTg42b958ydd3796Nm266CfPmzcO8efOwa9cuvdsPO2YYkXEEyi9JhkUqv1RvoLxeL44cOYI77rgDJSUliImJ8ds4EYWH1RHl99YTn8+HNWvWYMuWLSgvL0dZWRm++OKLS+puueUWlJaWorS0FPn5+aH4FsKKGUZkHIHyq6cMi2R+qd5ApaWlIS0tDRMmTAAAzJo1C0eOyN9NRkT6skXZ/d56Ul1djaFDhyI9PR0OhwO5ubmoqKgIQ8eRxQwjMo5A+dVThkUyv1S/uDhgwACkpaXhq6++wvDhw3HgwAFkZGTo2RsRBSHQoe6dO3di586d3R//40VzXS4X0tLSuj92Op1+L2myZ88evP/++7j66qtRXFyMq666Ssfuw48ZRmQcl3up7nIZFsn80nR21r/927/hwQcfhMfjQXp6OtavX6+5ISJSxxrgZMt/3DCpMXXqVMyZMwcOhwMvvfQSfvGLX+D555/XtKYRMMOIjCFQfgHaMyxU+aVpAzVmzBjs3r1bcxNEpJ3k5Tp/nE4n6urquj92uVzdF9n9Tr9+/br/nZ+fjyeeeEJdkwbDDCMyBjPmV1jfH+gd+n1R3fF22eCsIXGyCb0dFoeoziocs5rRr+caILgJ46sTrhXVPeD6VFQXJ7yuUHK07GcdJxxm5rM6ey4C4Jg8T1SH1jOwtjeJSj8feLOobnhf2cPeJxzQ29ghnNIeE9o/N4tN3SmN48aNQ21tLU6cOAGn04ny8nL86le/uqimvr4eqampAIC9e/f2ype6vEO+J6o73tHzifuAPL/aIVvPJg0wyDNMOmE81i677xMtXlGdNL+ShPkVGyX/2/DFDRDVSTPMeubSE5r9keZXRqzs8eDtkv3upPkFAP2F2amGGfNL9U/jq6++wooVK7o/PnHiBAoLC3H33Xfr0VevIt089UbSzRNd/hyCy7Hb7Xj44YexePFi+Hw+LFiwACNHjsSvf/1rXHfddZg2bRpeeOEF7N27FzabDYmJiVfES13MMCLjMGN+WRRFkV8IJwCfz4esrCy8/PLLGDRoUMC6zpbzovXER6D6yq5vpPcRKCm38BBGMBsovY9ASa/PJj0CZfPKJjnbWlyiumA2UJ/3vUZU59T5CFSHV/Y4DOYIVGKs7OjE3zu/eZXfzyctWRf0Wr2NJMM6m8+J1hIfgRLmVyiOQElJMyxSR6Aswm85IYjLgdg8sllgttYzojrxEXRhfqXpfARKml+A/AhUfF/98gswbobpcjzuwIEDSE9Pv+zmiYhCS+0zOGKGEUWaGfNLl47Ly8sxZ84cPZYiIpUu9y4WujxmGFFkmTG/NF+9z+12Y+/evZg1a5Ye/RCRSlZ7lN8bXR4zjCjyAuWXkTNM85avqqoKY8eORf/+/fXoh4jUMnDQGBkzjMgATJhfmjdQ5eXlyM3N1aMXItLAEiV7swRdjBlGFHlmzC9NL+FduHAB+/fvx4wZM/Tqh4hUstij/N4oMGYYkTEEyi8jZ5imI1B9+/bFoUOH9OqFiLSwyt+qTd9ihhEZhAnzy3ynvRORX0Z+pkZEdDlmzC9NG6ht27Zh165dsFgsGDVqFNavX48+ffoErO+wBv7a38s4/5Gorrmv7NIK8fVHRXWnEkeJ6lJjZDvlDuEgRulwTAD4lXOcqO7XTbKfoeXzd0R11vTRojpPf9mIfF+/oaK6riT5nNf+btlAOOnAwURFNhRU6St7XHtEVRrYzXcOQaQFk2Ht1mjRmhnNl14J3p/zMeNFdUlnZJeEOpkk+xsFgNRo2d+AG/oOyEyPl/0vx+puE9XhyD5RmW2wLNsBwNN/hKjOrXOGSfNLOi812dopqpPmFxDiDDNhfqk+B8rlcuH555/Hq6++irKyMvh8PpSXl+vZGxEFwRIV5fdG/jHDiIwjUH4ZOcM0HYHy+Xzo6OiA3W5HR0dH98X6iCj8LCZ8BhdpzDAiYzBjfqneQDmdTtxzzz2YOnUq+vTpg8mTJ+Pmm2VXkyYi/ZnxHIJIYoYRGYcZ80v1S3hNTU2oqKhARUUF3n33XbS3t6O0tFTP3ogoGPYo/zfyixlGZCCB8svAGaZ6A7V//34MHjwYycnJiIqKwowZM/Dxxx/r2RsRBcES5fB7I/+YYUTGESi/jJxhqjdQAwcOxCeffIL29nYoioIDBw4gI0P2DiwiCgGrzf+N/GKGERlIoPwycIapPgdqwoQJmDlzJn74wx/CbrdjzJgxKCgo0LM3IgqCGc8hiCRmGJFxmDG/NL0Lr7CwEIWFhXr1QkRaWDkXN1jMMCKDMGF+ma9jIvJLMWEAEREB5syvsHYc42kR1SldsomsCXWyib9HY68V1Y20tIrq4Jb92KKjYkV1cQ75qWjSCeP3J14vqnu2TjaJvKmvU1Tn6/SJ6to8st+xw6bpetd+JbaeEtV11coeX9bhsp+1LZjX8vsOkdd234FxzxW4EvT1yvJBccsm2Ce5IpRfAOCRZZhD5wyTThjvcsju13b1RFHd+dirRHUA0OWRTQ5vdcumr+udYbrnV8b3xfdtswi/l16SX5p+s9u3b8ecOXOQm5uLbdu26dQSEamhWO1+bxQYM4zIGALll5EzTPUG6q9//St27dqFXbt2obS0FJWVlfj666/17I2IgmG1+7+RX8wwIgMJlF8GzjDVG6gvv/wS48ePR0xMDOx2OzIzM7Fnzx49eyOiIJjt2VukMcOIjKNXHYEaNWoUPvzwQzQ2NqK9vR1VVVWoq6vTszciCobJZqhEGjOMyEB60xyojIwMLF68GPfeey9iYmIwevRoWK36n/BLRDKK9ARPAsAMIzISM+aXpmNj+fn5yM/PBwA89dRTcDpl79QiohCwGfdQt1Exw4gMwoT5pWnL19DQAAD45ptvsGfPHsydO1eXpohIBQ0nYFZVVWHmzJnIycnB5s2bL/m62+3G8uXLkZOTg/z8fJw8eVLv7iOCGUZkEBpOIo9Ufmna8i1btgznz5+H3W7HI488goSEBF2aIqLgqT3Z0ufzYc2aNfiv//ovOJ1OLFy4ENnZ2RgxYkR3za5du5CQkIC33noL5eXlePLJJ7Fx40a9Wo8YZhiRMZgxvzRtoF588UXNDRCRTlSebFldXY2hQ4ciPT0dAJCbm4uKioqLAmjv3r34+c9/DgCYOXMm1qxZA0VRYLFYtPcdQcwwIoMwYX6F9UVHxeYQ1el9UcHEPrJXKlvtsmefsXCL6k61eER1ydHyB47lc9nkcOmE8WVpU0V10gnozZZoUV0/4fd8rkM22RwA4h2yNZvjBonqEtJlv+cPOhJFddcnyiYXqxXoGdzOnTuxc+fO7o8LCgouumiuy+VCWlpa98dOpxPV1RdPMXa5XLjqqm+nOdvtdsTHx6OxsRHJycl6fguGptj7iOosDtnfgFS8zvkFALEWWTZJMyxJmmFH9onKpBPGfYkDRXVxnnZRHQA0IzIZpnt+DZHlzQft8sdNKDPsckegLpdhkcyvHjdQxcXFqKysREpKCsrKygAA58+fx4oVK3Dq1CkMGjQIGzduRGKi7H8iRBQaXfD/bOofN0y9DTOMyPgC5Rdg3Azr8alNXl4etmzZctHnNm/ejEmTJmHPnj2YNGmS35O2iCi8fIri99YTp9N50fwjl8t1ybvRnE4nTp8+DQDwer1oaWlBv3799P0GQoQZRmR8gfKrpwyLZH71uIHKzMy85JlZRUUF5s+fDwCYP38+3n77bc2NEJE2vi7F760n48aNQ21tLU6cOAG3243y8nJkZ2dfVJOdnY3XXnsNAPCnP/0JN910k2nOf2KGERlfoPzqKcMimV+qzoFqaGhAamoqAGDAgAHdbwUmosgR7JX8stvtePjhh7F48WL4fD4sWLAAI0eOxK9//Wtcd911mDZtGhYuXIiHHnoIOTk5SExMxNNPP61v82HGDCMyFjPml+aTyC0Wi2meiRJdyXwqAwgApkyZgilTplz0ufvvv7/733369MEzzzyj/g4MjBlGFHlmzC9VgzRTUlJQX18PAKivr+9V78QhMiq150D1RswwImNRew5UJKnaQGVnZ6OkpAQAUFJSgmnTpunaFBEFz9fl/0aXYoYRGUug/DJyhvW4gSoqKsLtt9+OmpoaZGVlYdeuXViyZAnee+89zJgxA/v378eSJUvC0SsRXYbZnr2FCzOMyPjMeASqx3OgnnrqKb+f3759u+7NEJF6knfc9UbMMCLjM2N+hXcSuXCSry9edkV0X0Jaz0UABp6QTdFuKH9FVNf1s/WiuuRo2QTaOOEEWgCwpo8W1TX1lf0MpRPG70+8XlQ358ghUd0NA+NEdUM8dT0X/Q9Le4e4VkI6Of8Gi0tU12ofKr5v2T1fTMtJmNQz6eNBmkvSnBt86hNR3dmSHaI6APDdt0FUlxIje/0kNkp2Noht8ChR3fnYq0R10gnjSlSMqA4ADtU2i+r0zrCI5Zddll8A0GyTZZiaMwrNmF89PuqLi4sxadIkzJkzp/tzb775JnJzczF69Gh8+umnIW2QiGTMdvg7XJhhRMZnxpfwVE0iHzVqFJ599llkZmaGrDEiCo6i+L/1dswwIuMLlF9GzrAeX8LLzMzEyZMnL/pcRkZGyBoiInWM/G6VSGKGERmfGfMrrOdAEVHoGPlQNxHR5Zgxv7iBIrpCeMx4FiYREcyZX9xAEV0hPF0mPAZORARz5hc3UERXCDM+gyMiAsyZXz1uoIqKinD48GE0NjYiKysLy5YtQ1JSEh599FGcO3cOS5cuxZgxY/C73/0uHP0SUQBmfAYXDswwIuMzY36pnkSek5OjezNEpJ4Zn8GFAzOMyPjMmF9hfQnPcfpzUZ2vbz9R3YkWr6hu8JDvi+qsi78nqrNdOCeqS+5oEdX5rLKJxADg6S97+7WvUzYFvdkSLaqTThgvu/ZGUd2sU38S1dVFDxLVAcDA9mOius9jRorqhiTIJvl2emXPnOJsqq7dLeYx4aUQzESaX10xiaI6aX6lD5ogqrMtHS+qAwB721lRXT93m6jOFzdAVOfpP0JU1+WRPZabIcsv6XRxAJg+LEFUZ288Lqr7Rphhg3XOr2GJsvzq8MpzI9YeugwzY371uIEqLi5GZWUlUlJSUFZWBgDYsGED3nnnHURFRWHIkCFYv349EhJkDzoiCo0uE74NOByYYUTGZ8b8UjWJfPLkySgrK8Prr7+OYcOGYdOmTSFrkIhkPL4uv7fejhlGZHyB8svIGdbjBiozMxOJiRcfkr755ptht3978GrixImoq5Nf8JWIQsPTpfi99XbMMCLjC5RfRs4wzedAvfrqq5g9e7YevRCRBmY8CdMImGFEkWfG/NK0gXruuedgs9lw66236tUPEalkxrcBRxozjMgYzJhfqjdQu3fvRmVlJbZt2waLxaJnT0SkghmfwUUSM4zIOMyYX6o2UFVVVdiyZQt+//vfIyYmRu+eiEgFMz6DixRmGJGxmDG/VE0i37x5M9xuNxYtWgQAmDBhAtasWRPyZokoMLdwHlVvwwwjMj4z5peqSeT5+fkhaYaI1PMZ+N0qkcQMIzI+M+ZXWCeRd1a/Kyuc9TNR2QDhjtUtfG21sUM2vTvpwhlRXedh2bRtx+R5ojoA8PUbKqpr88h+Nv2ibaK6GwbGieqkE8Z/PmimqA4A/u9f/iCq8ybLfjajvbIpzGU1fUV1U4fJJk9/eb5TVAcA42JkE5b/nhmfwZlJx8eVssJbfi4qGyCcbyMcyo0m4dUHACCprUFUp3eGuYX51eqWTWnXO78A+YRxb78horrBJz6Srad3fn0ly68pQ2X5BcgzrLfkl6pJ5Bs3bkRFRQWsVitSUlKwfv16OJ3yy5EQSUk3TxSaADp//jxWrFiBU6dOYdCgQdi4ceMlM5UAYMyYMRg1ahQA4KqrrsJvfvMb3XtRixlGZHyh2kCFMsNUTSJfvHgxXn/9dZSWluIHP/gB/vM//1P6vRBRiLi9XX5vWmzevBmTJk3Cnj17MGnSJGzevNlvXXR0NEpLS1FaWmqozRPADCMyg0D5ZeQMUzWJPC7u/x8ObW9v51uAiQyg09vl96ZFRUUF5s+fDwCYP38+3n77bT1aDStmGJHxBcovI2eY6nOgnn76aZSUlCA+Ph7PP/+8bg0RkTqBnqnt3LkTO3fu7P64oKAABQUFojUbGhqQmpoKABgwYAAaGvyfO9PZ2Ym8vDzY7XYsWbIE06dPD7L78GOGERnH5Y40GTXDVG+gVqxYgRUrVmDTpk34/e9/j8LCQrVLEZEO3AFOSu4pbO6++26cPXvpianLly+/6GOLxRLwSM0777wDpxMee9kAABVeSURBVNOJEydO4K677sKoUaMwZIjsJNtIYYYRGUeg/AKMm2Ga34U3d+5cLFmyhOFDFGFqzxXYtm1bwK+lpKSgvr4eqampqK+vR3Jyst+6707ATk9Pxz/90z/hyJEjht9AfYcZRhR5Ws51ilSG9XgOlD+1tbXd/66oqMDw4cPVLENEOnJ7fX5vWmRnZ6OkpAQAUFJSgmnTpl1S09TUBLfbDQA4d+4cPvroI4wYMULT/YYaM4zIWALll5EzTNUk8qqqKtTU1MBisWDQoEH4j//4j2C/JyLSWSgG0S1ZsgTLly/HK6+8goEDB2Ljxo0AgE8//RQvvfQS1q5diy+//BKPPPIILBYLFEXBT37yE0NtoJhhRMYXqkGaocwwTiInukJofbeKP/369cP27dsv+fy4ceMwbtw4AMD111+P119/Xff71gszjMj4QpFfQGgzLKyTyK05i2WF7gvCBWXTTqMhm2o7xN4mqvOkjhLV2Wami+p8AKLOfCGq7UqS7dIdNtmrs+eE09eHeOpEdXXRg0R1wQzI/N/X/FhUt+ybalHd6A7Zz3piWqqoLgYeUV2UNbRvlTfjJF8zsc1aKiv06JtfDkWWX4NtsvwC9M8w6xWSXwDwjTDDpBPG3enXi+r+dt4tqru2UzaJXJpfsVb5S2ShzDAz5lePj9Li4mJMmjQJc+bMueRrW7duxTXXXINz586FpLneQrp56o2kmycKzSDNKwEzjMj4QjVIM5RUTSIHgNOnT+O9997DwIEDQ9IYEQXH5+3ye+vtmGFExhcov4ycYaomkQPA+vXr8dBDD3GCL5FBdPm6/N56O2YYkfEFyi8jZ5iqc6DefvttpKamYvTo0Xr3Q0Qq+byheRfLlYgZRmQsZsyvoDdQ7e3t2LRpE7Zu3RqKfohIJZ+Bn6kZCTOMyHjMmF9BD9I8fvw4Tp48iXnz5iE7Oxt1dXXIy8vDmTNnQtEfEQl1ebv83uhizDAi4wmUX0bOsKCPQF1zzTU4cOBA98fZ2dl45ZVXAo5HJ6Lw6FLMdwg8EphhRMZjxvzq8QhUUVERbr/9dtTU1CArKwu7du0KR19EFCSzPXsLF2YYkfFdkUeg/E3x/Xt79+7VrRkiUs+M5xCEAzOMyPjMmF9hnURORKFjxnexEBEB5syvHjdQxcXFqKysREpKCsrKygAAzz77LF5++eXucwaKioowZcqUHu/M0t4kauq4cunMFn8Gx8nmt9i//rPsfvtPFNUNcf1FVKc0nBLVfT7wZlEdAPR367tLj3fYRHWW9g5R3cD2Y6I6b/JQUZ308iwA8OzA8aK6Nec/F9X1j5K9xyLqm09EdSP69hPVfSshiNpvGflQdyTplWHS/Pq6S5Zf6fHC/Kr5UFR3PPX7ojogchlm9PwCgME6Z5j0Ei0jkxyiujPukaI6vfMLCCbDekd+9biBysvLw5133olf/OIXF33+7rvvxr333huyxogoOGY8BB4OzDAi4zNjfvW4gcrMzMTJkyfD0QsRaWDkSx5EEjOMyPjMmF9Bz4H6zh/+8AfMnTsXxcXFaGqSHdomotAx22UQIo0ZRmQcZryUi6oN1B133IG33noLpaWlSE1NxWOPPaZ3X0QUpK4uxe+NLsUMIzKWQPll5AxTtYHq378/bDYbrFYr8vPz8emnn+rdFxEFqcvr9nujSzHDiIwlUH4ZOcNUjTGor69HamoqgG8vyjlypOxdAUQUOl0e4waN0TDDiIzFjPnV4waqqKgIhw8fRmNjI7KysrBs2TIcPnwYx459+1bPQYMGYc2aNSFvlIguz8jP1CKJGUZkfGbML1WTyPPz80PSDBGpZ8YACgdmGJHxmTG/OImc6ArhM2EAEREB5syvsG6gTkE2oXeIrUVU51ZkU1FtSQNFdenub0R1vrj+orqPLENEdcP7yn8NNqtsenFiq2yCcHPcIPF9S3weIzuXZLT3rKyu4wvxfUsnjD+cNFZUt7b5iKiuzxnZz7r1OtmkdABIFlf+f2Z8Bmcm0vwaGtUmqnMLr7hg7yf7G5XmF6B/hmXERonqhPEVsfwC9M+waztlddIJ40kO2Xu/2oSXRpHmFyDPsN6SXz3+JoqLizFp0iTMmTPnos+/8MILmDVrFnJzc/H444+HrEEikunyuP3eejtmGJHxBcovI2eYqku5HDx4EBUVFfjjH/8Ih8OBhoaGkDZJRD0z4zO4cGCGERmfGfOrxyNQmZmZSEy8+FDzjh07sGTJEjgc3178MCUlJTTdEZGY0tXl99bbMcOIjC9Qfhk5w1QN0qytrcUHH3yA/Px83Hnnnaiurta7LyIKUiiG0L355pvIzc3F6NGjLztssqqqCjNnzkROTg42b96s6T7DgRlGZCyhGqQZygxTtYHy+XxoamrCyy+/jJUrV2L58uVQFOOOWyfqDXwet9+bFqNGjcKzzz6LzMzMwPfr82HNmjXYsmULysvLUVZWhi++kJ/8HwnMMCJjCZRfRs4wVe/CczqdyMnJgcViwfjx42G1WtHY2IjkZDXn3hORHkJxDkFGRkaPNdXV1Rg6dCjS09MBALm5uaioqMCIESN070cvzDAiYwnVOVChzDBVG6jp06fj0KFDuOmmm1BTUwOPx4N+/WQjBYgoNDre/43fz+/cuRM7d+7s/rigoAAFBQW63a/L5UJaWlr3x06n0/AviTHDiIwlUH4Bxs0wVZdyWbBgAVatWoU5c+YgKioKjz32GCwW4YAPIgqrnsLm7rvvxtmzl86qWb58OaZPnx7K1sKCGUZkbkbNMFWXcgGAJ598UvdmiCj8tm3bpum/dzqdqKur6/7Y5XLB6XRq7Eo/zDCiK1ukMiysk8jPXvCK6gZHy14L7fTJTvo875CF+VXfHBLVYZBsknVStOzHK/w2AACJSoeorqtW9hJKQrrsZ63YHKK6IQmyurKavqK6iWmpojoA6B8le0+EdML46oRrRXXPnP9AVJdU/5moDgAQ/0/y2ggbN24camtrceLECTidTpSXl+NXv/pVpNvSnTy/ZH+jnfYEUV1TtOxKCmmnDorqAOieYd4uWYglWztFdeL8GiL7nUjzCwCGJQoz7Ct9M0yaX9IJ47F24RHVa6fI6hBEhpkovwD1GdbjX0dxcTEqKyuRkpKCsrIyAN8eFqupqQEAtLS0ID4+HqWlpRq/BSIymrfeeguPPvoozp07h6VLl2LMmDH43e9+B5fLhX/913/Fb3/7W9jtdjz88MNYvHgxfD4fFixYgJEjZZelCAdmGFHvFcoMUzWJfOPGjd3/fuyxxxAXF6fyWyMiI8vJyUFOTs4ln3c6nfjtb3/b/fGUKVMwZYr8mWw4McOIeq9QZpiqSeTfURQFb7755iXXmCIiMgpmGBGFgqpBmt/54IMPkJKSgmHDhunUDhFR+DDDiEgtTRuosrIyPnMjItNihhGRWqo3UF6vF2+99RZuueUWPfshIgoLZhgRaaF6A7V//34MHz78oumdRERmwQwjIi163EAVFRXh9ttvR01NDbKysrBr1y4AwBtvvIHc3NyQN0hEpAUzjIhCQfUk8scee0z3ZoiI9MYMI6JQsCiKEsQcbG3c574R1f2tUzaTZVC8bEpujMUnqmv0yl7RjBVOjPUJp/M2dsj6A4Cr+sru29ZSL6r7sMP/27v/0Q12l6juXNwQUZ3DJpuSGwOPqA4AouqOiuq8Z07JFhRO6C1MukFU93T7Mdn9AoiJjhbXUni4G+t6LgLwtw7ZhOpI5RcQuQwT51frGVHdh+2yae7S/AKAxvihojphhCHWKvvZRJ3+XFSnd351OWJl6wUhOiZG9zWNqMdHc3FxMSZNmnTRO1WOHj2K2267DfPmzUNeXp7hr7xORL0XM4yIQqHHDVReXh62bNly0eeeeOIJ3HfffSgtLcX999+PJ554ImQNEhFpwQwjolBQNYncYrGgra0NwLfXkUpNlV/wlYgonJhhRBQKshfh/8GqVatw7733YsOGDejq6sJLL72kd19ERCHDDCMirVTNgdqxYweKi4uxb98+FBcXY/Xq1Xr3RUQUMswwItJK1Qbqtddew4wZMwAAs2fP5gmYRGQqzDAi0krVBio1NRWHDx8GABw8eJAX4iQiU2GGEZFWPZ4DVVRUhMOHD6OxsRFZWVlYtmwZHn30Uaxbtw5erxd9+vTBmjVrwtErEVHQmGFEFAqqJ5Hv3r1b92aIiPTGDCOiUAjrJPKOC22iuq+aZZNbo+2yUbCf1LWK6q7uJ5uemp7gENVJJ9VaLMJCAFZhaXSbbJKvEiWbeH3BLpsO77DJXhX+8nynqC5K+g0DGOE7Laprjk8X1SXVfyaq6xw0XlS3Ima0qA4AfqPUimspPKT59WWTLL9ioiKTX4A8w6R/fVHSsBPSO79abbL8AoBoe2QyLFL55R44TlQXDE4i/x/+pvgeO3YMBQUFmDt3Ln7605+itVX2B05EFG7MMCIKBVWTyFevXo0HHngAr7/+OqZPn37J14mIjIIZRkShoGoSeW1tLTIzMwEAkydPxp49e0LTHRGRRswwIgoFVWMMRo4ciYqKCgDAf//3f+P0adlrt0RERsAMIyKtVG2g1q5dixdffBF5eXloa2uDwyE7IZGIyAiYYUSklapr4WVkZGDr1q0AgJqaGlRWVurZExFRSDHDiEgrVUegGhoaAABdXV147rnncPvtt+vaFBFRKDHDiEgrVZPIL1y4gBdffBEAkJOTgwULFoS8USIiNZhhRBQKqieR33XXXbo3Q0SkN2YYEYWCqpfwiIiIiHqzsF7KhYiIiOhKwCNQREREREHiBoqIiIgoSNxAEREREQWJGygiIiKiIHEDRURERBQkbqCIiIiIghTRDVRVVRVmzpyJnJwcbN68WdNap0+fxr/8y7/glltuQW5uLrZv365Ljz6fD/Pnz8fSpUs1r9Xc3IzCwkLMmjULs2fPxscff6x5zW3btiE3Nxdz5sxBUVEROjs7g/rvi4uLMWnSJMyZM6f7c+fPn8eiRYswY8YMLFq0CE1NTZrX3LBhA2bNmoW5c+fivvvuQ3Nzs6b1vrN161Zcc801OHfunOb1XnjhBcyaNQu5ubl4/PHHNa139OhR3HbbbZg3bx7y8vJQXV0tXo/MobflF6B/hmnNL0D/DOtt+RVoTWaYgBIhXq9XmTZtmnL8+HGls7NTmTt3rvK3v/1N9Xoul0v57LPPFEVRlJaWFmXGjBma1vvO1q1blaKiImXJkiWa11q5cqXy8ssvK4qiKJ2dnUpTU5Om9erq6pSpU6cq7e3tiqIoSmFhofLqq68Gtcbhw4eVzz77TMnNze3+3IYNG5RNmzYpiqIomzZtUh5//HHNa7777ruKx+NRFEVRHn/88aDW9LeeoijKN998o9xzzz3KD37wA6WhoUHTegcOHFDuuusupbOzU1EURTl79qym9RYtWqRUVlYqiqIolZWVyp133ilej4yvN+aXouibYXrkl6Lon2G9Lb8CrckM61nEjkBVV1dj6NChSE9Ph8PhQG5uLioqKlSvl5qairFjxwIA4uLiMHz4cLhcLk091tXVobKyEgsXLtS0DgC0tLTg/fff717L4XAgISFB87o+nw8dHR3wer3o6OhAampqUP99ZmYmEhMTL/pcRUUF5s+fDwCYP38+3n77bc1r3nzzzbDbv71y0MSJE1FXV6dpPQBYv349HnroIVgsFs397dixA0uWLIHD4QAApKSkaFrPYrGgra0NwLe/+2B/L2RsvS2/gNBkmNb8AvTPsN6WX4HWZIb1rMdr4YWKy+VCWlpa98dOp1O3Q4QnT57E0aNHMWHCBE3rrFu3Dg899FD3g0hrT8nJySguLsaxY8cwduxYrF69Gn379lW9ptPpxD333IOpU6eiT58+mDx5Mm6++WbNvTY0NHT/sQwYMKD7yvV6efXVVzF79mxNa7z99ttITU3F6NGjdemptrYWH3zwAZ5++mn06dMHK1euxPjx41Wvt2rVKtx7773YsGEDurq68NJLL+nSJxlDb8uv7/rSM8NClV9AaDOsN+QXwAyTuOJOIm9ra0NhYSFWrVqFuLg41eu88847SE5OxnXXXadLX16vF0eOHMEdd9yBkpISxMTEaD5voqmpCRUVFaioqMC7776L9vZ2lJaW6tLvdywWS9DPkC7nueeeg81mw6233qp6jfb2dmzatAn333+/bn35fD40NTXh5ZdfxsqVK7F8+XIoGq5ytGPHDhQXF2Pfvn0oLi7G6tWrdeuVrlxGzS9A/wwLR34B+mZYb8kvgBkmEbENlNPpvOgwqMvlgtPp1LSmx+NBYWEh5s6dixkzZmha66OPPsLevXuRnZ2NoqIiHDx4EA8++KDq9dLS0pCWltb9rHLWrFk4cuSIph7379+PwYMHIzk5GVFRUZgxY4YuJ6anpKSgvr4eAFBfX4/k5GTNawLA7t27UVlZiSeffFJToB0/fhwnT57EvHnzkJ2djbq6OuTl5eHMmTOq13Q6ncjJyYHFYsH48eNhtVrR2Nioer3XXnut+zE4e/ZsnoB5helt+QXon2Ghyi8gNBnWm/ILYIZJRGwDNW7cONTW1uLEiRNwu90oLy9Hdna26vUURcHq1asxfPhwLFq0SHN/DzzwAKqqqrB371489dRTuOmmm/Dkk0+qXm/AgAFIS0vDV199BQA4cOAAMjIyNPU4cOBAfPLJJ2hvb4eiKLqsCQDZ2dkoKSkBAJSUlGDatGma16yqqsKWLVvw3HPPISYmRtNa11xzDQ4cOIC9e/di7969SEtLw+7duzFgwADVa06fPh2HDh0CANTU1MDj8aBfv36q10tNTcXhw4cBAAcPHsSwYcNUr0XG09vyC9A/w0KVX4D+Gdbb8gtghklYFK3H+TTYt28f1q1bB5/PhwULFuBnP/uZ6rU++OAD/PjHP8aoUaNgtX67LywqKsKUKVM093no0CFs3boVmzZt0rTO0aNHsXr1ang8HqSnp2P9+vV+Ty4MxjPPPIM33ngDdrsdY8aMwdq1a7tPJJQoKirC4cOH0djYiJSUFCxbtgzTp0/H8uXLcfr0aQwcOBAbN25EUlKSpjU3b94Mt9vdvc6ECROwZs0a1evl5+d3fz07OxuvvPKK+Fmmv/XmzZuHVatW4dixY4iKisLKlSsxadIk1etdffXVWLduHbxeL/r06YNHHnlE15dTKPJ6W34B+meY1vwC9M+w3pZfgdZkhvUsohsoIiIiIjO64k4iJyIiIgo1bqCIiIiIgsQNFBEREVGQuIEiIiIiChI3UERERERB4gaKiIiIKEjcQBEREREFiRsoIiIioiD9PwRBs4QqW9T4AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 720x360 with 4 Axes>"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Now lets measure the similarity \n",
        "print(stats.spearmanr(upper(m1), upper(m2)))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OS_ArNqZDlDq",
        "outputId": "f1b8bf6d-6a0e-43b4-f3c3-de5301ca45ef"
      },
      "execution_count": 78,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "SpearmanrResult(correlation=0.14991317735875465, pvalue=0.038974535666595055)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "\"\"\"Nonparametric permutation testing Monte Carlo\"\"\"\n",
        "np.random.seed(0)\n",
        "rhos = []\n",
        "n_iter = 5000\n",
        "true_rho, _ = stats.spearmanr(upper(m1), upper(m2))\n",
        "# matrix permutation, shuffle the groups\n",
        "m_ids = list(m1.columns)\n",
        "m2_v = upper(m2)\n",
        "for iter in range(n_iter):\n",
        "  np.random.shuffle(m_ids) # shuffle list \n",
        "  r, _ = stats.spearmanr(upper(m1.loc[m_ids, m_ids]), m2_v)  \n",
        "  rhos.append(r)\n",
        "perm_p = ((np.sum(np.abs(true_rho) <= np.abs(rhos)))+1)/(n_iter+1) # two-tailed test\n",
        "\n",
        "\n",
        "f,ax = plt.subplots()\n",
        "plt.hist(rhos,bins=20)\n",
        "ax.axvline(true_rho,  color = 'r', linestyle='--')\n",
        "ax.set(title=f\"Permuted p: {perm_p:.3f}\", ylabel=\"counts\", xlabel=\"rho\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "id": "IwysAEWcDk7K",
        "outputId": "5cc43df1-38b0-4258-f071-371e2494e97d"
      },
      "execution_count": 85,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de1xVdb7/8dcOvKACCuFGkS4amaOldjQkj9hAeMMLmGZppnTRJpOwsoc6lT7GvHVs5NE0WpzpFFZnphPpptQUxQs2mjVHGSu19HhMJfamEEXHVMD1+2NP66dHxC269trq+/l4rMdmL9ba+/21R3z29/td67sdhmEYiIiIANfZHUBERAKHioKIiJhUFERExKSiICIiJhUFERExqSiIiIhJRUHEAlu2bCExMdHuGCIXTUVBbJWUlMQdd9xB165dufvuu5kyZQr/+Mc/7I7FlClTWLBggd0xzuudd96hZ8+e3HnnnUydOpVTp06d99jNmzfTr18/OnfuzOjRoykpKTnnmMOHD9OjRw8efPBBc19xcTEZGRncdddd9OjRg8zMTMrKyixpjwQOFQWx3RtvvMG2bdtYunQpX3/9NYsWLbqo8w3D4PTp0xalCzwbN24kJyeHd955h3Xr1nHw4EFee+21Wo89dOgQTz31FE8//TRffPEFnTp1YtKkSeccN3/+fNq1a3fWviNHjnD//fezdu1a1q1bR9OmTZk6daolbZLAoaIgAcPpdNKrVy92794NeD+pPvDAA3Tr1o3BgwezZcsW89jRo0ezYMECHnjgATp37syBAwdo374977//Pn369KFr165kZ2ezf/9+HnjgAe68806efvpp8xP1kiVLzvpUDNC+fXu+//57PvjgAz755BPeeustunbtyhNPPAGAx+Nh4sSJ9OjRg6SkJBYvXmyee+LECaZMmUL37t0ZMGAAX331VZ1tbd++PYsXLyY5OZn4+HjmzZvnc2FzuVwMGzaMuLg4wsPDefLJJ1m6dGmtx65evZq4uDj69+9Po0aNmDhxIrt27eJ//ud/zGO2bt3K7t27GTp06Fnn9u7dm/79+9OsWTNCQkJ46KGH2Lp1q08Z5coVbHcAkV+UlpZSVFRESkoKHo+H8ePH88orr9CrVy82b95MZmYmn376KREREQDk5+fz7//+79x88838slrLZ599xpIlSygtLSU9PZ1t27bxb//2bzRv3pwRI0awfPly0tPT68wxYsQItm3bhtPpND9Vnz59mt/85jckJSXx6quv4vF4GDt2LDfffDO9evXi9ddfZ//+/axevZqff/6Zxx9//ILtXb16NR999BHHjx8nIyODtm3bMnz4cH744QcGDx7Mxx9/TOvWrc85b/fu3SQnJ5vP27dvz08//URFRQUtWrQ459j27dubz5s0acINN9zAnj17aNeuHTU1NcycOZOZM2fy3Xff1Zn3yy+/JC4u7oLtkiubegpiuwkTJtCtWzdGjhxJ9+7deeKJJ8jPzycxMZHevXtz3XXX0bNnTzp16sSGDRvM89LT04mLiyM4OJgGDRoA8Nhjj9GsWTPi4uK49dZb6dmzJ7GxsYSGhpKYmMiOHTvqlfGrr74yh2IaNmxIbGws999/PytWrADg008/5YknnqB58+a0atWK0aNHX/A1H3/8cZo3b07r1q15+OGHWbZsGQCtW7fmb3/7W60FAeD48eM0a9bMfB4aGgpQ61zM8ePHzd//olmzZuax7777LnfccQedOnWqM+uuXbtYuHAhzz///AXbJVc29RTEdn/84x+5++67z9r3ww8/sHLlStatW2fuq66uJj4+3nzeqlWrc17r+uuvN39u1KjROc9/+umnemUsKSmhrKyMbt26mftqamrM52VlZWflOd8f9DOdeXxMTIzPk7hNmjTh2LFj5vNffm7atOkFjwVv8WjatCkej4fFixezZMmSOt/v+++/5/HHH2fatGlntV+uTioKEpBatWrFkCFDePnll897jMPhqPfrh4SEcOLECfP5jz/+WOdrt2rVijZt2lBQUFDr60VFRVFaWmoOr5SWll4ww5nH//DDD7Rs2dKn7HFxcXz77bcMGDAA8H6Kv/76688ZOvrl2DPnG44fP87+/fu55ZZb+Oqrr/jxxx9JTU0FvPMiJ0+epGfPnhQVFREUFERJSQkZGRk8+eSTpKWl+ZRPrmwaPpKANHjwYNatW8fGjRupqanh5MmTbNmyBbfbfVle/7bbbmP37t3s3LmTkydP8oc//OGs30dGRnLw4EHz+R133EHTpk3JycnhxIkT1NTU8N1337F9+3YA+vfvT05ODkeOHMHtdvPuu+9eMMNbb73FkSNHKC0tZfHixeYf+QsZMmQIeXl57Nmzh8rKShYtWnTeeZKUlBR2797NqlWrOHnyJH/84x9p37497dq1IzExkbVr1+JyuXC5XGRmZtKhQwdcLhdBQUF4PB7GjBnDqFGjzpmUl6uXioIEpFatWrFw4ULefPNNEhIS6N27N2+99dZlu/T05ptvZsKECYwdO5Y+ffrwL//yL2f9ftiwYezZs4du3brx5JNPEhQUxBtvvMGuXbtITk6mR48evPDCC+bQzFNPPUXr1q1JTk7mkUceYciQIRfMkJyczNChQ0lLS+Oee+5h2LBhgLfX0LVrV3744Ydaz0tMTOSxxx7j4Ycf5p577iEmJobMzEzz96mpqXz88ccARERE8Ic//IEFCxbQvXt3tm/fzu9//3sAGjZsSFRUlLmFhoYSHBxMVFQUAB9++CEHDhzg9ddfp2vXruYmVzeHvmRHxP/at29PQUEBN954o91RRM6inoKIiJhUFERExKThIxERMamnICIipiv6PoX4+HhiYmLsjiEiF+OX5TRuvdXeHNewkpKSs9YSO9MVXRRiYmIueDemiASYe+7xPur/Xdv838UPz6ThIxERMV3RPQURuQK98ILdCaQOKgoi4l/33mt3AqmDho9ExL+Ki72bBCT1FETEv7KyvI/r19saQ2qnnoKIiJhUFERExKSiICIiJhUFERExWTbRvHfvXiZNmmQ+P3DgAJmZmaSlpTFp0iRKSkqIiYkhOzub8PBwDMNg1qxZbNiwgcaNGzN37lw6duxoVTy5RpyoqqFxgyC/nyt1mD3b7gRSB8uKQtu2bcnPzwe8X3CemJhISkoKOTk5JCQkMG7cOHJycsjJyWHy5MkUFRWxb98+CgoK+Pvf/86MGTP48MMPrYon14jGDYK4acryep27b27qZU4jANx9t90JpA5+GT7avHkzsbGxxMTEUFhYaH4BeFpaGmvWrAEw9zscDrp06UJlZSVlZWX+iCci/rRpk3eTgOSX+xSWL1/OwIEDASgvL6dly5YAREVFUV5eDoDH4yE6Oto8Jzo6Go/HYx4rIleJadO8j7pPISBZ3lM4deoUa9eupV+/fuf8zuFw4HA4rI4gIiI+srwoFBUV0bFjR66//noAIiMjzWGhsrIyIiIiAHA6nbjdbvM8t9uN0+m0Op6IiJzB8qKwfPlyUlP//4RdUlISLpcLAJfLRXJy8ln7DcOguLiY0NBQDR2JiPiZpUXh+PHjbNq0iT59+pj7xo0bx1//+lf69OnDpk2bGDduHAC9e/cmNjaWlJQUXnzxRaZPn25lNBERqYWlE81NmjQ55yvfWrRoQW5u7jnHOhwOFQKRa0F2tt0JpA5aJVVE/KtLF7sTSB20zIWI+NeaNd5NApJ6CiLnoSUyLPLyy95HfQNbQFJREDkPLZEh1yINH4mIiElFQURETCoKIiJi0pyCiPjXm2/anUDqoKIgIv7Vvr3dCaQOGj4SEf/65BPvJgFJPQURC1zqfQpX9X0Or77qfRw0yN4cUisVBRELXMo9DqD7HMQ+Gj4SERGTioKIiJhUFCTgnaiqsTuCyDVDcwoS8LQG0VXm3XftTiB1UFEQEf+KjbU7gdRBw0ci4l8ffODdJCCppyAi/rVokfdxxAh7c0it1FMQERGTpUWhsrKSzMxM+vXrR//+/dm2bRuHDx8mIyODPn36kJGRwZEjRwAwDIOXX36ZlJQUBg0axDfffGNlNBERqYWlRWHWrFn06tWLlStXkp+fT7t27cjJySEhIYGCggISEhLIyckBoKioiH379lFQUMDMmTOZMWOGldFERKQWlhWFo0eP8uWXXzJs2DAAGjZsSFhYGIWFhaSlpQGQlpbGmn9+gfcv+x0OB126dKGyspKysjKr4omISC0sm2g+ePAgERERTJ06lV27dtGxY0d++9vfUl5eTsuWLQGIioqivLwcAI/HQ3R0tHl+dHQ0Ho/HPFZErhJ5eXYnkDpY1lOorq5mx44dPPjgg7hcLkJCQsyhol84HA4cDodVEUQkEF1/vXeTgGRZUYiOjiY6OprOnTsD0K9fP3bs2EFkZKQ5LFRWVkZERAQATqcTt9ttnu92u3E6nVbFExG7vPOOd5OAZFlRiIqKIjo6mr179wKwefNm2rVrR1JSEi6XCwCXy0VycjKAud8wDIqLiwkNDdXQkcjVSEUhoFl689qLL77Ic889R1VVFbGxscyZM4fTp0+TlZVFXl4erVu3Jjs7G4DevXuzYcMGUlJSCAkJYfbs2VZGExGRWlhaFDp06MCSJUvO2Z+bm3vOPofDwfTp062MIyIiF6A7mkVExKSiICIiJi2IJyL+tWKF3QmkDioKIuJfTZrYnUDqoOEjEfGvhQu9mwQkFQUR8a//+i/vJgFJRUFEREwqCiIiYlJREBERk4qCiIiYdEmqiPjX+vV2J5A6qKcgIiImFQUR8a/5872bBCQVBRHxr2XLvJsEJBUFERExqSiIiIhJRUFEREy6JFVE/CskxO4EUgcVBRHxr08/tTuB1EHDRyIiYrK0p5CUlETTpk257rrrCAoKYsmSJRw+fJhJkyZRUlJCTEwM2dnZhIeHYxgGs2bNYsOGDTRu3Ji5c+fSsWNHK+OJiB1mzvQ+vviivTmkVpb3FHJzc8nPz2fJkiUA5OTkkJCQQEFBAQkJCeTk5ABQVFTEvn37KCgoYObMmcyYMcPqaCJih8JC7yYBye/DR4WFhaSlpQGQlpbGmjVrztrvcDjo0qULlZWVlJWV+TueiMg1zfKi8OijjzJ06FA++OADAMrLy2nZsiUAUVFRlJeXA+DxeIiOjjbPi46OxuPxWB1PRETOYOmcwp///GecTifl5eVkZGTQtm3bs37vcDhwOBxWRhARkYtgaU/B6XQCEBkZSUpKCtu3bycyMtIcFiorKyMiIsI81u12m+e63W7zfBG5ikRGejcJSJYVhePHj3Ps2DHz57/+9a/ExcWRlJSEy+UCwOVykZycDGDuNwyD4uJiQkNDzWEmEbmKfPSRd5OAZNnwUXl5ORMmTACgpqaGgQMHkpiYyO23305WVhZ5eXm0bt2a7OxsAHr37s2GDRtISUkhJCSE2bNnWxVNRETOw7KiEBsby8cff3zO/hYtWpCbm3vOfofDwfTp062KIzY7UVVD4wZBdse4YlzKv1fA/1tPnep9nDPH3hxSKy1zIX7RuEEQN01ZXq9z981NvcxpAt9V/e+1ebPdCaQOWuZCRERMKgoiImJSURAREZPmFETEv9q0sTuB1EFFQUT867337E4gddDwkYiImFQURMS/srK8mwQkDR+JiH8VF9udQOqgnoKIiJhUFERExKSiICIiJs0piIh/3Xqr3QmkDioKIuJfOTl2J5A6aPhIRERMKgoi4l/jxnk3CUgaPhIR//ruO7sTSB3UUxAREZNPRSE3N5djx45hGAbTpk0jPT2dzz77zOpsIiLiZz4VhY8++ohmzZrx2WefUVlZySuvvMKrr75qdTYREfEzn4qCYRgAbNiwgSFDhhAXF2fuu5CamhrS0tIYP348AAcOHGD48OGkpKSQlZXFqVOnADh16hRZWVmkpKQwfPhwDh48WJ/2iEig69LFu0lA8qkodOrUiUceeYSioiL+9V//lWPHjnHddb5NRyxevJh27dqZz+fPn8/YsWNZvXo1YWFh5OXlAfDhhx8SFhbG6tWrGTt2LPPnz69Hc0Qk4GVnezcJSD79ZZ81axbPPvsseXl5hISEUFVVxezZsy94ntvtZv369QwbNgzw9jg+//xz+vbtC0B6ejqFhYUArF27lvT0dAD69u3L5s2bfe6NiIjI5eFTUcjIyKBjx46EhYUB0KJFC+bMmXPB82bPns3kyZPNXkVFRQVhYWEEB3uvhI2Ojsbj8QDg8Xho1aoVAMHBwYSGhlJRUXHxLRKRwPbQQ95NAlKd9ymcPHmSn3/+mYqKCo4cOWJ+cj927Jj5x/x81q1bR0REBJ06dWLLli2XL7GIXNk0XxjQ6iwKf/nLX8jNzaWsrIyhQ4eaRaFZs2Y8dIFKv3XrVtauXUtRUREnT57k2LFjzJo1i8rKSqqrqwkODsbtduN0OgFwOp2UlpYSHR1NdXU1R48epUWLFpepmSIi4os6i8KYMWMYM2YM7777LqNHj76oF3722Wd59tlnAdiyZQv/8R//wauvvkpmZiarVq0iNTWVpUuXkpSUBEBSUhJLly6la9eurFq1ih49euBwOOrZLBERqQ+flrkYPXo0W7dupaSkhJqaGnN/WlraRb/h5MmTmTRpEtnZ2XTo0IHhw4cDMGzYMCZPnkxKSgrh4eEsWLDgol9bREQujU9FYfLkyRw4cIDbbruNoKAgABwOh89FIT4+nvj4eABiY2PNy1DP1KhRI1577TVfc4vIlSohwe4EUgefisLXX3/NihUrNJwjIpfOhysXxT4+XZIaFxfHjz/+aHUWERGxmU89hYqKClJTU7njjjto0KCBuf+NN96wLJiIXKXuu8/7+NFH9uaQWvlUFCZOnGh1DhG5VpSX251A6uBTUbjrrrusziEiIgHAp6LQtWtXc5K5qqqK6upqQkJC2Lp1q6XhRETEv3wqCtu2bTN/NgyDwsJCiouLLQslIiL2uOiv43Q4HNx777365jURqZ/kZO8mAcmnnkJBQYH58+nTp/n6669p1KiRZaFE5Cr24ot2J5A6+FQU1q1bZ/4cFBRETEwMCxcutCyUiIjYw6ei4Mt3J4iI+KR/f+/jp5/am0Nq5dOcgtvtZsKECSQkJJCQkMDEiRNxu91WZxORq9HPP3s3CUg+FYWpU6eSlJTExo0b2bhxI7/+9a+ZOnWq1dlERMTPfCoKhw4d4r777iM4OJjg4GCGDh3KoUOHrM4mIvVwoqrmwgdZcK5cHXyaU2jevDn5+fkMHDgQgGXLltG8eXNLg4lI/TRuEMRNU5bX69x9c1Mvcxq50vhUFGbPns3MmTOZM2cODoeDrl27MnfuXKuzicjV6J8fLiUw+VQUXnvtNebNm0d4eDgAhw8fZt68eboqSUQu3nPP2Z1A6uDTnMK3335rFgTwDift3LnTslAiImIPn4rC6dOnOXLkiPn88OHDZ31Xs4iIz+65x7tJQPJp+OiRRx5hxIgR9OvXD4CVK1fyxBNPWBpMRET8z6eikJaWRqdOnfj8888BeP3117nlllvqPOfkyZOMGjWKU6dOUVNTQ9++fcnMzOTAgQM888wzHD58mI4dO/LKK6/QsGFDTp06xfPPP88333xD8+bNWbBgAW3atLn0FoqIiM98KgoAt9xyywULwZkaNmxIbm4uTZs2paqqipEjR5KYmMjbb7/N2LFjSU1N5aWXXiIvL4+RI0fy4YcfEhYWxurVq1m+fDnz588nOzu7Xo0SEZH6ueils33lcDho2rQpANXV1VRXV+NwOPj888/p27cvAOnp6RQWFgKwdu1a0tPTAejbty+bN2/GMAyr4omISC187inUR01NDUOHDmX//v2MHDmS2NhYwsLCCA72vm10dDQejwcAj8dDq1atvKGCgwkNDaWiooKIiAgrI4qIv91/v90JpA6WFoWgoCDy8/OprKxkwoQJ7N2718q3E5ErwZNP2p1A6mDZ8NGZwsLCiI+Pp7i4mMrKSqqrqwHv6qtOpxMAp9NJaWkp4B1uOnr0KC1atPBHPPGB1sSRy+b4ce8mAcmynsKhQ4cIDg4mLCyMEydOsGnTJh5//HHi4+NZtWoVqampLF26lKSkJACSkpJYunQpXbt2ZdWqVfTo0QOHw2FVPLlIl7KeDmhNHTnDgAHex/XrbY0htbOsKJSVlTFlyhRqamowDIN+/frx61//mltuuYVJkyaRnZ1Nhw4dGD58OADDhg1j8uTJpKSkEB4ezoIFC6yKJiIi52FZUbjttttwuVzn7I+NjSUvL++c/Y0aNeK1116zKo6IiPjAL3MKIiJyZVBREBERk6WXpIqInGPsWLsTSB1UFETEv1QUApqGj0TEv376ybtJQFJPQUT8a9gw76PuUwhI6imIiIhJRUFEREwqCiIiYlJREBERkyaaRcS/fvMbuxNIHVQURMS/RoywO4HUQcNHIuJfBw54NwlI6imIiH+NHu191H0KAUk9BRERMakoiIiISUVBRERMKgoiYjpRVWPLuRI4NNEsIqbGDYK4acryep27b26qbwc++2y9Xl/8Q0VBRPxr0CC7E0gdLBs+Ki0tZfTo0QwYMIDU1FRyc3MBOHz4MBkZGfTp04eMjAyOHDkCgGEYvPzyy6SkpDBo0CC++eYbq6KJiJ2+/da7SUCyrCgEBQUxZcoUVqxYwQcffMB//ud/smfPHnJyckhISKCgoICEhARycnIAKCoqYt++fRQUFDBz5kxmzJhhVTQRsdP48d5NApJlRaFly5Z07NgRgGbNmtG2bVs8Hg+FhYWkpaUBkJaWxpo1awDM/Q6Hgy5dulBZWUlZWZlV8UREpBZ+ufro4MGD7Ny5k86dO1NeXk7Lli0BiIqKory8HACPx0N0dLR5TnR0NB6Pxx/xRETknywvCv/4xz/IzMxk2rRpNGvW7KzfORwOHA6H1RFERMRHlhaFqqoqMjMzGTRoEH369AEgMjLSHBYqKysjIiICAKfTidvtNs91u904nU4r44mIyP9hWVEwDIPf/va3tG3bloyMDHN/UlISLpcLAJfLRXJy8ln7DcOguLiY0NBQc5hJRK4iL7zg3SQgWXafwn//93+Tn5/PrbfeypAhQwB45plnGDduHFlZWeTl5dG6dWuys7MB6N27Nxs2bCAlJYWQkBBmz55tVTQRsdO999qdQOpgWVHo1q0b357nWuRf7lk4k8PhYPr06VbFEZFAUVzsfezSxd4cUivd0Swi/pWV5X3U9ykEJC2IJyIiJhUFERExqSiIiIhJRUFEREyaaL6GnKiqoXGDILtjyLVOl5sHNBWFa4hfvkBF5ELuvtvuBFIHDR+JiH9t2uTdJCCppyAi/jVtmvdR9ykEJPUURETEpKIgIiImFQURETGpKIiIiEkTzSLiX/9cLl8Ck4qCiPiXlswOaBo+EhH/WrPGu0lAUk9BRPzr5Ze9j/oGtoCknoKIiJhUFETksjhRVWPr+XJ5aPhIRC4LXxdc/MvecgAe+D/HatHFwKCegoiImCzrKUydOpX169cTGRnJsmXLADh8+DCTJk2ipKSEmJgYsrOzCQ8PxzAMZs2axYYNG2jcuDFz586lY8eOVkUTERtN6/uU3RGkDpb1FIYOHcqf/vSns/bl5OSQkJBAQUEBCQkJ5OTkAFBUVMS+ffsoKChg5syZzJgxw6pYImKzvZFt2BvZxu4Ych6WFYXu3bsTHh5+1r7CwkLS0tIASEtLY80/r1X+Zb/D4aBLly5UVlZSVlZmVTQRsVHyni0k79lidww5D79ONJeXl9OyZUsAoqKiKC/3Tjh5PB6io6PN46Kjo/F4POaxInL1ePyLpQAU3hJvcxKpjW0TzQ6HA4fDYdfbi4hILfxaFCIjI81hobKyMiIiIgBwOp243W7zOLfbjdPp9Gc0ERHBz0UhKSkJl8sFgMvlIjk5+az9hmFQXFxMaGioho5ERGxg2ZzCM888wxdffEFFRQWJiYlMnDiRcePGkZWVRV5eHq1btyb7n0vo9u7dmw0bNpCSkkJISAizZ8+2KpaIiNTBsqLw+9//vtb9ubm55+xzOBxMnz7dqihXlRNVNTRuEGR3DJF6mzTwWbsjSB20zMUVxtelBGqjZQQkEJSGRdkdQeqgZS5ExK8G7ixi4M4iu2PIeainICJ+9dC2FQAs65BocxKpjXoKIiJiUlEQERGTioKIiJhUFERExKSJZhHxq9+kTa11/6Xcg6P7dy4fFQUR8auKJuG17tc9OIFBw0c20BeUy7Vs2FdrGPbVGrtjyHmop2ADfSKSa9kvBSHv9nttTiK1UU9BRERMKgoiImJSURAREZOKgoiImDTRLCJ+NXb4DLsjSB3UUxARvzrRoDEnGjS+vK95CZd56xLxs6mnICJ+9dBW7+XY7915+S6v1mXel496CvWkTxci9TNw10YG7tpodww5D/UU6kmfTETkahRQPYWioiL69u1LSkoKOTk5dscREbnmBExRqKmp4Xe/+x1/+tOfWL58OcuWLWPPnj12xxIRuaYETFHYvn07N954I7GxsTRs2JDU1FQKCwstez/NCYgI2Hfl0qX+DbLqb5jDMAzDkle+SCtXrmTjxo3MmjULAJfLxfbt23nppZfOe058fDwxMTH+iigiclUoKSlhy5Yttf7uip5oPl+jRESkfgJm+MjpdOJ2u83nHo8Hp9NpYyIRkWtPwBSF22+/nX379nHgwAFOnTrF8uXLSUpKsjuWiMg1JWCGj4KDg3nppZd47LHHqKmp4b777iMuLs7uWCIi15SAmWgWERH7BczwkYiI2E9FQURETCoKl9Hhw4fJyMigT58+ZGRkcOTIkXOO2blzJyNGjCA1NZVBgwaxYsUKG5LWny9tBHj00Ufp1q0b48eP93PC+rvQMiunTp0iKyuLlJQUhg8fzsGDB21IWX8Xat+XX35Jeno6v/rVr1i5cqUNCS/dhdr49ttvM2DAAAYNGsSYMWMoKSmxIeWluVAb//znPzNo0CCGDBnCgw8+ePErQxhy2cybN8948803DcMwjDfffNN45ZVXzjlm7969xv/+7/8ahmEYbrfb6Nmzp3HkyBF/xrwkvrTRMAxj06ZNRmFhoTFu3Dh/xqu36upqIzk52di/f79x8uRJY2qhi00AAATUSURBVNCgQcbu3bvPOua9994zXnzxRcMwDGPZsmXG008/bUfUevGlfQcOHDB27txpTJ482fj0009tSlp/vrRx8+bNxvHjxw3DMIz333//ivpvaBi+tfHo0aPmz2vWrDEeeeSRi3oP9RQuo8LCQtLS0gBIS0tjzZo15xxz8803c9NNNwHeezMiIiI4dOiQP2NeEl/aCJCQkEDTpk39Ge2S+LLMytq1a0lPTwegb9++bN68GeMKuU7Dl/a1adOG2267jeuuuzL/LPjSxh49ehASEgJAly5dzro36krgSxubNWtm/vzzzz/jcDgu6j0C5pLUq0F5eTktW7YEICoqivLy8jqP3759O1VVVdxwww3+iHdZXGwbrxQej4fo6GjzudPpZPv27ecc06pVK8B7CXVoaCgVFRVERET4NWt9+NK+K93FtjEvL4/ExER/RLtsfG3j+++/z9tvv01VVRW5ubkX9R4qChdp7Nix/PTTT+fsz8rKOuu5w+Gos0KXlZUxefJk5s2bF3CfzC5XG0UCVX5+Pl9//TXvvfee3VEsMWrUKEaNGsUnn3zCokWLmDdvns/nqihcpHfeeee8v4uMjKSsrIyWLVtSVlZ23k+Qx44dY/z48UyaNIkuXbpYlLT+LkcbrzS+LLPidDopLS0lOjqa6upqjh49SosWLfwdtV6uhWVkfG3jpk2beOONN3jvvfdo2LChPyNesov975iamsqMGTMu6j0C6yPqFS4pKQmXywV4V3lNTk4+55hTp04xYcIEhgwZQr9+/fwd8ZL50sYrkS/LrCQlJbF06VIAVq1aRY8ePa6YntK1sIyML23csWMHL730EosWLSIyMtKmpPXnSxv37dtn/rx+/XpuvPHGi3uTyzEjLl6HDh0yHn74YSMlJcUYM2aMUVFRYRiGYWzfvt2YNm2aYRiG4XK5jF/96lfG4MGDzW3Hjh12xr4ovrTRMAzjwQcfNOLj443bb7/d6NWrl1FUVGRXZJ+tX7/e6NOnj5GcnGwsXLjQMAzDyM7ONtasWWMYhmGcOHHCmDhxonHvvfca9913n7F//3474160C7Xv73//u9GrVy+jc+fOxl133WUMGDDAzrj1cqE2jhkzxkhISDD/3xs/frydcevlQm2cOXOmMWDAAGPw4MHGQw89ZHz33XcX9fpa5kJEREwaPhIREZOKgoiImFQURETEpKIgIiImFQURETGpKIhcBgcPHmTgwIF2xxC5ZCoKIpeBYRicPn3a7hgil0z3KYjU08GDB3n00Ufp3LkzBQUFREREcPfdd7Nt2zacTicLFy6kcePG7Ny5k+nTp/Pzzz9zww03MHv2bMLDw+2OL1Ir9RRELsH333/PyJEjWbZsGW63m1GjRrF8+XJCQ0NZtWoVAM8//zzPPfccn3zyCbfeeiuvv/66zalFzk9FQeQStG7d2lzUsE2bNnTo0AGAjh07UlJSwtGjRzl69Ch33XUXAOnp6fztb3+zLa/IhagoiFyCJk2amD+fueJmUFAQNTU1dkQSuSQqCiIWCg0NJSwszOwd5Ofn0717d5tTiZyfvk9BxGLz5s0zJ5pjY2OZM2eO3ZFEzktXH4mIiEnDRyIiYlJREBERk4qCiIiYVBRERMSkoiAiIiYVBRERMakoiIiI6f8BYJC5JqDGRd8AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Real life example comparing city distances and their average yearly temperature differences. "
      ],
      "metadata": {
        "id": "AvNyohP09FQt"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\"\"\"Real data of including latitude, longitude, average temperature, and city name. \n",
        "\n",
        "References:\n",
        "  https://keisan.casio.com/exec/system/1224587128\n",
        "  https://en.wikipedia.org/wiki/List_of_cities_by_average_temperature\n",
        "\"\"\"\n",
        "seoul = (37.566535, 126.977969, 12.5, \"Seoul\")\n",
        "tokyo = (35.676192, 139.650311, 15.4, \"Tokyo\")\n",
        "sf = (37.77493, -122.419415, 14.6, \"San Francisco\")\n",
        "la = (34.052234, -118.243685, 18.6, \"Los Angeles\")\n",
        "chicago = (41.878114, -87.629798, 9.8, \"Chicago\")\n",
        "beijing = (39.904211, 116.407395, 12.9, \"Beijing\")\n",
        "cairo = (30.04442, 31.235712, 21.4, \"Cairo\")\n",
        "moscow = (55.755826, \t37.6173, 5.8, \"Moscow\")\n",
        "sydney = (-33.86882, 151.209296, 17.7, \"Sydney\")\n",
        "london =(51.507218, -0.127586, 11.3, \"London\")\n",
        "paris = (48.856614, 2.352222, 12.3, \"Paris\")\n",
        "berlin = (52.520007, 13.404954, 10.3, \"Berlin\")\n",
        "rome = (41.902784, 12.496366, 15.2, \"Rome\")\n",
        "lima = (-12.046373, -77.042754, 19.2, \"Lima\")\n",
        "mexicocity = (19.432608, -99.133208, 17.5, \"Mexico City\")\n",
        "capetown = (-33.924868, 18.424055, 16.2, \"Cape Town\")\n",
        "newdelhi = (28.613939, 77.209021, 25, \"New Delhi\")\n",
        "auckland = (-36.850883, 174.764488, 15.2, \"Auckland\")\n",
        "saopaulo = (-23.555771, -46.639557, 19.2, \"Sao Paulo\")\n",
        "madrid = (40.416775, -3.70379, 15, \"Madrid\")\n",
        "bangkok = (13.756331, 100.501765, 28.6, \"Bangkok\")\n",
        "phnompenh = (11.556374,104.92821, 28.3, \"Phnom Penh\")\n",
        "mumbai = (19.075984, 72.877656, 27.1, \"Mumbai\")\n",
        "singapore = (1.352083, 103.819836, 27, \"Singapore\")\n",
        "vladivostok = (43.133248, 131.911297, 4.9, \"Vladivostok\")\n",
        "city_list = [seoul, tokyo, sf, la, chicago, beijing, cairo, moscow, sydney, london, paris, \n",
        "             berlin, rome, lima, mexicocity, capetown, newdelhi, auckland, saopaulo, madrid,\n",
        "             bangkok, phnompenh, mumbai, singapore, vladivostok]"
      ],
      "metadata": {
        "id": "VL2ymxJgEc9O"
      },
      "execution_count": 91,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Measure similarity \n",
        "n = len(city_list)\n",
        "city_names = [city[-1] for city in city_list]\n",
        "m1 = pd.DataFrame(np.zeros((n,n)), index = city_names, columns = city_names)\n",
        "m2 = pd.DataFrame(np.zeros((n,n)), index = city_names, columns = city_names)\n",
        "for i, city1 in enumerate(city_list):\n",
        "  for j, city2 in enumerate(city_list):\n",
        "    # m1.iloc[i,j] = distance(city1[0], city1[1], city2[0], city2[1])\n",
        "    m1.iloc[i,j] = np.abs(city1[0]-city2[0]) \n",
        "    m2.iloc[i,j] = np.abs(city1[2]-city2[2])\n",
        "\n",
        "stats.spearmanr(upper(m1), upper(m2))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7OMZvcQdEc0h",
        "outputId": "21776e13-24c8-4ff8-c3a8-6315b22526ee"
      },
      "execution_count": 92,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "SpearmanrResult(correlation=0.21213613537496678, pvalue=0.00021460054002915316)"
            ]
          },
          "metadata": {},
          "execution_count": 92
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Permute for significance\n",
        "rhos = []\n",
        "n_iter = 5000\n",
        "true_rho, _ = stats.spearmanr(upper(m1), upper(m2))\n",
        "# matrix permutation, shuffle the groups\n",
        "m_ids = list(m1.columns)\n",
        "m2_v = upper(m2)\n",
        "for iter in range(n_iter):\n",
        "  np.random.shuffle(m_ids) # shuffle list \n",
        "  r, _ = stats.spearmanr(upper(m1.loc[m_ids, m_ids]), m2_v)  \n",
        "  rhos.append(r)\n",
        "perm_p = ((np.sum(np.abs(true_rho) <= np.abs(rhos)))+1)/(n_iter+1) # two-tailed test\n",
        "\n",
        "f,ax = plt.subplots()\n",
        "plt.hist(rhos,bins=20)\n",
        "ax.axvline(true_rho,  color = 'r', linestyle='--')\n",
        "ax.set(title=f\"Permuted p: {perm_p:.3f}\", ylabel=\"counts\", xlabel=\"rho\")\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "id": "pf3K6t8oNgHs",
        "outputId": "df288457-90cf-47ab-b94f-bb4e8d31f0f5"
      },
      "execution_count": 93,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de1xUdeL/8dcoXkhBxWDwQu2mhK6a6GqIPtQWQk0lwTS7mZdau5ikln21rfSxpmlr6VqbymP7Gla73zaVodQUxQu2krWrrHmp9NvXVdEZCkE0QwXP74+p88sFYwY9DDO+n4/HeQxz5hx42yN98/mcm80wDAMRERGgnq8DiIhI3aFSEBERk0pBRERMKgURETGpFERExKRSEBERk0pBxAI7d+6kX79+vo4h4jWVgvhUQkICt9xyC926daN3795Mnz6d7777ztexmD59OgsXLvR1jMt666236NOnD927d2fGjBmcP3/+stvm5eUxaNAgunbtyujRoykoKDA/mz9/PgMGDKBbt24MGjQIh8NR5fdwOBzExMTw/vvvX/U/i9QtKgXxuaVLl7J7924yMzPZu3cvS5Ys8Wp/wzC4ePGiRenqnu3bt5Oens5bb73Fli1bOHbsGIsXL65y25MnT/LEE0/w5JNP8umnn9K5c2emTJlifh4cHMySJUv45z//yfz585kzZw67du265HucOnWKpUuXEh0dbemfS+oGlYLUGXa7nb59+3Lw4EEA8vPzueeee+jRowd33nknO3fuNLcdPXo0Cxcu5J577qFr164cPXqUmJgY3n33XfM330WLFnHkyBHuueceunfvzpNPPmn+Rr169WruvffeS35+TEwM//73v3nvvff48MMPefPNN+nWrRuPPvooAC6Xi0mTJtGrVy8SEhJYsWKFuW9ZWRnTp0+nZ8+eDB48mM8///xn/6wxMTGsWLGCxMRE4uLimD9/vsfF5nA4GDFiBNHR0TRr1ozHH3+czMzMKrfduHEj0dHR3HHHHTRq1IhJkybxxRdf8L//+78ApKWl0a5dO+rVq0fXrl359a9/TX5+/iXf45VXXmH06NG0aNHCo3zi31QKUmecOHGC3NxcOnbsiMvl4pFHHuGxxx7j008/5b/+679IS0vj5MmT5vZZWVnMnj2bXbt20bp1awA+/vhjVq9ezd/+9jf+/Oc/8/zzz/OHP/yBbdu2cfDgQdauXVttjlGjRpGcnMxDDz3E7t27Wbp0KRcvXuSxxx4jJiaG3NxcMjIyyMjIYPv27QC8/vrrHDlyhI0bN/Lmm29edhrmpzZu3MiqVavIzMxk8+bNrFq1CoDjx4/To0cPjh8/XuV+Bw8epEOHDub7mJgYvv32W4qLi6vcNiYmxnx/3XXXccMNN3Do0KFK25aVlbF3717at29vrtuzZw979+6tVKASuFQK4nMTJ06kR48e3HffffTs2ZNHH32UrKws+vXrR//+/alXrx59+vShc+fObNu2zdwvNTWV6OhogoKCaNCgAQAPP/wwTZs2JTo6mptvvpk+ffoQFRVFSEgI/fr1Y//+/TXK+Pnnn5tTMQ0bNiQqKoq7776bdevWAfDRRx/x6KOP0rx5c1q1asXo0aOr/Z6//e1vad68Oa1bt+bBBx9kzZo1ALRu3Zp//OMfZtH9p7Nnz9K0aVPzfUhICECVx2LOnj1rfv6jpk2bVrntzJkziYmJoW/fvgBUVFQwa9YsXnjhBerV0z8V14ogXwcQ+dOf/kTv3r0vWXf8+HHWr1/Pli1bzHXl5eXExcWZ71u1alXpe11//fXm140aNar0/ttvv61RxoKCAgoLC+nRo4e5rqKiwnxfWFh4SZ7L/YP+Uz/dvk2bNhQWFnqU5brrruPMmTPm+x+/btKkSbXbgrs8/nPb+fPnc/DgQVasWIHNZgPgL3/5CzExMcTGxnqUSwKDSkHqpFatWjFs2DBefPHFy27z4z9eNREcHExZWZn5/ptvvvnZ792qVSvatm1LdnZ2ld8vPDycEydOmAdjT5w4UW2Gn25//PhxIiIiPMoeHR3Nl19+yeDBgwH44osvuP7666uc84+Ojr7keMPZs2c5cuTIJVNEixcvZvv27bz99tuXjEDy8vL47LPPyM3NBdwHnPfv38+BAwd44YUXPMoq/kdjQqmT7rzzTrZs2cL27dupqKjg3Llz7Ny5E6fTeVW+f4cOHTh48CAHDhzg3LlzvPbaa5d83rJlS44dO2a+v+WWW2jSpAnp6emUlZVRUVHBV199xZ49ewC44447SE9P59SpUzidTt5+++1qM7z55pucOnWKEydOsGLFCvMf+eoMGzaMlStXcujQIUpLS1myZAmpqalVbpuUlMTBgwfZsGED586d409/+hMxMTG0a9cOgGXLlrFmzRqWL19eqVTmzZvHunXrcDgcOBwOOnfuzBNPPHHJ2UsSeFQKUie1atWKN954g2XLlhEfH0///v158803r9qpp7/85S+ZOHEiY8eOZcCAAfz617++5PMRI0Zw6NAhevToweOPP079+vVZunQpX3zxBYmJifTq1YvnnnvOnJp54oknaN26NYmJiYwfP55hw4ZVmyExMZHhw4eTkpLCbbfdxogRIwD3qKFbt26XPdDcr18/Hn74YR588EFuu+022rRpQ1pamvn5kCFD+OCDDwAICwvjtddeY+HChfTs2ZM9e/bw6quvmtu++uqrHD9+3Dxjq1u3bixduhSA0NBQwsPDzaVBgwY0bdq00jEKCSw2PWRHpPbFxMSQnZ3NjTfe6OsoIpfQSEFEREwqBRERMWn6SERETBopiIiIya+vU4iLi6NNmza+jiFSN331lfv15pt9m0PqnIKCgkvuJfZTfl0Kbdq0YfXq1b6OIVI33Xab+1V/R+Q/DB8+/LKfafpIRERMfj1SEJGf8dxzvk4gfkilIBKobr/d1wnED2n6SCRQ5ee7FxEvaKQgEqgmT3a/bt3q0xjiXzRSEBERk6Wl8NZbbzFkyBCGDh3K1KlTOXfuHEePHmXkyJEkJSUxefJk85m558+fZ/LkySQlJTFy5MhLblssIiK1w7JScLlcrFixglWrVrFmzRoqKipYu3YtCxYsYOzYsWzcuJHQ0FBWrlwJwPvvv09oaCgbN25k7NixLFiwwKpoIiJyGZaOFCoqKigrK6O8vJyysjLCw8P55JNPGDhwIOB+xm5OTg4AmzdvNh8UMnDgQPLy8tBtmUREapdlB5rtdjvjx4/nN7/5DY0aNaJPnz506tSJ0NBQgoLcPzYyMhKXywW4RxY/PrM2KCiIkJAQiouLCQsLsyriNafsQgWNG9Sv9X3FR+bO9XUC8UOWlcKpU6fIyckhJyeHkJAQnnzySbZv327VjxMPNG5Qn19MX1ujfQ/PG3KV04jlevf2dQLxQ5ZNH+3YsYO2bdsSFhZGgwYNGDBgALt27aK0tJTy8nIAnE4ndrsdcI8sfnzYeXl5OadPn67yQeQi4qEdO9yLiBcsK4XWrVvzr3/9i++//x7DMMjLy6N9+/bExcWxYcMGADIzM0lISAAgISGBzMxMADZs2ECvXr2w2WxWxRMJfM8+615EvGBZKXTt2pWBAweSmppKcnIyFy9eZNSoUUybNo3ly5eTlJRESUkJI0eOBNwPSi8pKSEpKYnly5fz9NNPWxVNREQuw9IrmtPS0khLS7tkXVRUlHka6k81atSIxYsXWxlHRESqoSuaRUTEpFIQERGTbognEqgWLfJ1AvFDKgWRQBUb6+sE4oc0fSQSqDZtci8iXtBIQSRQvfii+1VPYBMvaKTgZ8ouVPg6gogEMI0U/IzuXyQiVtJIQURETCoFERExafpIJFAtW+brBOKHVAoigSomxtcJxA9p+kgkUH34oXsR8YJGCiKB6pVX3K/Jyb7NIX5FIwURETGpFMQjV3rRnC66E/EPlk0fff3110yZMsV8f/ToUdLS0khJSWHKlCkUFBTQpk0bFi1aRLNmzTAMgzlz5rBt2zYaN27MvHnz6NSpk1XxxEtXctEc6MI5EX9h2UjhpptuIisri6ysLFavXk1wcDBJSUmkp6cTHx9PdnY28fHxpKenA5Cbm8vhw4fJzs5m9uzZzJo1y6poIiJyGbUyfZSXl0dUVBRt2rQhJyeHlJQUAFJSUtj0w10cf1xvs9mIjY2ltLSUwsLC2ognEpjeftu9iHihVkph7dq1DB06FICioiIiIiIACA8Pp6ioCACXy0VkZKS5T2RkJC6XqzbiiQSmqCj3IuIFy0vh/PnzbN68mUGDBlX6zGazYbPZrI4gcm167z33IuIFy0shNzeXTp06cf311wPQsmVLc1qosLCQsLAwAOx2O06n09zP6XRit9utjicSuJYscS8iXrC8FNauXcuQIf//zJOEhAQcDgcADoeDxMTES9YbhkF+fj4hISHmNJOIiNQOS0vh7Nmz7NixgwEDBpjrJkyYwN///ncGDBjAjh07mDBhAgD9+/cnKiqKpKQknn/+eWbOnGllNBERqYKlt7m47rrr2Llz5yXrWrRoQUZGRqVtbTabikBExMd0RbOIiJh0QzyRQLVypa8TiB9SKYgEqh/O+BPxhqaPRALVW2+5FxEvqBREApVKQWpApSAiIiaVgoiImFQKIiJiUimIiIhJp6SKBKp163ydQPyQSkEkUF13na8TiB/S9JFIoHrjDfci4gWVgkig+tvf3IuIF1QKIiJiUimIiIhJpSC1ouxChU/2FRHv6OwjqRWNG9TnF9PX1mjfw/OGVL+RiFwVlo4USktLSUtLY9CgQdxxxx3s3r2bkpISxo0bx4ABAxg3bhynTp0CwDAMXnzxRZKSkkhOTmbfvn1WRhMJfFu3uhcRL1haCnPmzKFv376sX7+erKws2rVrR3p6OvHx8WRnZxMfH096ejoAubm5HD58mOzsbGbPns2sWbOsjCYiIlWwrBROnz7NZ599xogRIwBo2LAhoaGh5OTkkJKSAkBKSgqbNm0CMNfbbDZiY2MpLS2lsLDQqngigW/BAvci4gXLSuHYsWOEhYUxY8YMUlJS+N3vfsfZs2cpKioiIiICgPDwcIqKigBwuVxERkaa+0dGRuJyuayKJxL41qxxLyJesKwUysvL2b9/P/feey8Oh4Pg4GBzquhHNpsNm81mVQQREfGSZaUQGRlJZGQkXbt2BWDQoEHs37+fli1bmtNChYWFhIWFAWC323E6neb+TqcTu91uVTwREamCZaUQHh5OZGQkX3/9NQB5eXm0a9eOhIQEHA4HAA6Hg8TERABzvWEY5OfnExISYk4ziYhI7bD0OoXnn3+ep59+mgsXLhAVFcVLL73ExYsXmTx5MitXrqR169YsWrQIgP79+7Nt2zaSkpIIDg5m7ty5VkYTCXzBwb5OIH7I0lLo2LEjq1evrrQ+IyOj0jqbzcbMmTOtjCNybfnoI18nED+k21yIiIhJpSASqGbPdi8iXlApiASqnBz3IuIFlYKIiJhUCiIiYlIp+ICeDyAidZWep+ADeraA1IqWLX2dQPyQSkEkUK1a5esE4oc0fSQiIiaVgkigmjHDvYh4QdNHIoEqL8/XCcQPaaQgIiImlYKIiJhUCiIiYtIxBZFA1batrxOIH1IpiASqd97xdQLxQ5o+EhERk6UjhYSEBJo0aUK9evWoX78+q1evpqSkhClTplBQUECbNm1YtGgRzZo1wzAM5syZw7Zt22jcuDHz5s2jU6dOVsYTCWyTJ7tff3jkrYgnLB8pZGRkkJWVZT6WMz09nfj4eLKzs4mPjyc9PR2A3NxcDh8+THZ2NrNnz2bWrFlWRxMJbPn57kXEC7U+fZSTk0NKSgoAKSkpbNq06ZL1NpuN2NhYSktLKSwsrO14IiLXNMtL4aGHHmL48OG89957ABQVFREREQFAeHg4RUVFALhcLiIjI839IiMjcblcVscTEZGfsPSYwl//+lfsdjtFRUWMGzeOm2666ZLPbTYbNpvNyggiIuIFS0vBbrcD0LJlS5KSktizZw8tW7aksLCQiIgICgsLCQsLM7d1Op3mvk6n09xfRGrg5pt9nUD8kGXTR2fPnuXMmTPm13//+9+Jjo4mISEBh8MBgMPhIDExEcBcbxgG+fn5hISEmNNMIlID6enuRcQLlo0UioqKmDhxIgAVFRUMHTqUfv360aVLFyZPnszKlStp3bo1i344Xa5///5s27aNpKQkgoODmTt3rlXRRETkMiwrhaioKD744INK61u0aEFGRkal9TabjZkzZ1oVR+TaM2GC+1WjBfGCbnMhEqi++srXCcQP6TYXIiJi8qgUMjIyOHPmDIZh8Oyzz5KamsrHH39sdTYREallHpXCqlWraNq0KR9//DGlpaW8/PLLvPLKK1ZnExGRWubRMQXDMADYtm0bw4YNIzo62lwnInVUbKyvE4gf8qgUOnfuzPjx4zl27BhPPfUUZ86coV49HY4QqdN0d1SpAY9KYc6cORw4cICoqCiCg4MpLi7WdQQiIgHIo1/3x40bR6dOnQgNDQXc1xq89NJLlgYTkSv0wAPuRcQLPztSOHfuHN9//z3FxcWcOnXKPI5w5swZ3cFUpK47dszXCcQP/Wwp/M///A8ZGRkUFhYyfPhwsxSaNm3KA/oNREQk4PxsKYwZM4YxY8bw9ttvM3r06NrKJCIiPuLRgebRo0eza9cuCgoKqKioMNf/+AQ1EREJDB6VwrRp0zh69CgdOnSgfv36gPsGdioFkTosPt7XCcQPeVQKe/fuZd26dXpKmog/0RmCUgMenZIaHR3NN998Y3UWERHxMY9GCsXFxQwZMoRbbrmFBg0amOuXLl1qWTARuUJ33eV+XbXKtznEr3hUCpMmTbI6h4hcbUVFvk4gfsijUrj11ltr/AMqKiq46667sNvtLFu2jKNHjzJ16lRKSkro1KkTL7/8Mg0bNuT8+fM888wz7Nu3j+bNm7Nw4ULatm1b458rIiLe8+iYQrdu3ejevTvdu3enS5cudOzYke7du3v0A1asWEG7du3M9wsWLGDs2LFs3LiR0NBQVq5cCcD7779PaGgoGzduZOzYsSxYsKAGfxwREbkSHpXC7t272bVrF7t27WLPnj289tpr3HfffdXu53Q62bp1KyNGjADct+D+5JNPGDhwIACpqank5OQAsHnzZlJTUwEYOHAgeXl5uj23AFB2oaL6jSzYV+Ra5PUzmm02G7fffjuvv/46Tz/99M9uO3fuXKZNm8Z3330HuA9Yh4aGEhTk/rGRkZHmPZRcLhetWrVyhwoKIiQkhOLiYsLCwryNKAGmcYP6/GL62hrte3jekKucxo8kJvo6gfghj0ohOzvb/PrixYvs3buXRo0a/ew+W7ZsISwsjM6dO7Nz584rSyki3nv+eV8nED/kUSls2bLF/Lp+/fq0adOGN95442f32bVrF5s3byY3N5dz585x5swZ5syZQ2lpKeXl5QQFBeF0OrHb7QDY7XZOnDhBZGQk5eXlnD59mhYtWlzBH01ERLzlUSnU5NkJTz31FE899RQAO3fu5L//+7955ZVXSEtLY8OGDQwZMoTMzEwSEhIASEhIIDMzk27durFhwwZ69eqlK6hFrsQdd7hfP/rItznEr3h0oNnpdDJx4kTi4+OJj49n0qRJOJ3OGv3AadOmsXz5cpKSkigpKWHkyJEAjBgxgpKSEpKSkli+fHm1xytEpBrff+9eRLzg0UhhxowZDB06lD/+8Y8AfPDBB8yYMYPly5d79EPi4uKIi4sDICoqyjwN9acaNWrE4sWLPc0tIiIW8GikcPLkSe666y6CgoIICgpi+PDhnDx50upsIiJSyzwqhebNm5OVlUVFRQUVFRVkZWXRvHlzq7OJiEgt82j6aO7cucyePZuXXnoJm81Gt27dmDdvntXZRORKDB3q6wTihzwqhcWLFzN//nyaNWsGQElJCfPnz6/RWUkiUkt0sobUgEfTR19++aVZCOCeTjpw4IBloURExDc8KoWLFy9y6tQp831JScklz2oWkTrottvci4gXPJo+Gj9+PKNGjWLQoEEArF+/nkcffdTSYCIiUvs8KoWUlBQ6d+7MJ598AsDrr79O+/btLQ0mIiK1z+O7pLZv315FICIS4Dw6piAiItcGr5+nICJ+4u67fZ1A/JBKQSRQPf64rxOIH9L0kUigOnvWvYh4QSMFkUA1eLD7detWn8YQ/6KRgoiImFQKIiJiUimIiIjJsmMK586d4/777+f8+fNUVFQwcOBA0tLSOHr0KFOnTqWkpIROnTrx8ssv07BhQ86fP88zzzzDvn37aN68OQsXLqRt27ZWxRMRkSpYNlJo2LAhGRkZfPDBBzgcDrZv305+fj4LFixg7NixbNy4kdDQUPPRnO+//z6hoaFs3LiRsWPHsmDBAquiiVwbxo51LyJesKwUbDYbTZo0AaC8vJzy8nJsNhuffPIJAwcOBCA1NZWcnBwANm/eTGpqKgADBw4kLy8PwzCsiicS+FQKUgOWHlOoqKhg2LBh9O7dm969exMVFUVoaChBQe5Zq8jISFwuFwAul4tWrVoBEBQUREhICMXFxVbGEwls337rXkS8YOl1CvXr1ycrK4vS0lImTpzI119/beWPE5GfGjHC/arrFMQLtXL2UWhoKHFxceTn51NaWkp5eTkATqcTu90OgN1u58SJE4B7uun06dO0aNGiNuKJiMgPLCuFkydPUlpaCkBZWRk7duygXbt2xMXFsWHDBgAyMzNJSEgAICEhgczMTAA2bNhAr169sNlsVsUTEZEqWDZ9VFhYyPTp06moqMAwDAYNGsRvfvMb2rdvz5QpU1i0aBEdO3Zk5MiRAIwYMYJp06aRlJREs2bNWLhwoVXR5BpSdqGCxg3q1/q+Iv7KslLo0KEDDoej0vqoqCjzNNSfatSoEYsXL7YqjlyjGjeozy+mr63RvofnDbnKaUTqPt0QTyRQPfaYrxOIH1IpiASqUaN8nUD8kO59JBKojh51LyJe0EihhnQQUuq80aPdr7pOQbygUqghHcAUkUCk6SMRETGpFERExKRSEBERk44piASqp57ydQLxQyoFkUCVnOzrBOKHNH0kEqi+/NK9iHhBIwWRQPXII+5XXacgXtBIQURETCoFERExqRRERMSkUhAREZMONIsEquee83UC8UOWlcKJEyd45plnKCoqwmazcffddzNmzBhKSkqYMmUKBQUFtGnThkWLFtGsWTMMw2DOnDls27aNxo0bM2/ePDp16mRVPJHAd/vtvk4gfsiy6aP69eszffp01q1bx3vvvcdf/vIXDh06RHp6OvHx8WRnZxMfH096ejoAubm5HD58mOzsbGbPns2sWbOsiiZybcjPdy8iXrCsFCIiIszf9Js2bcpNN92Ey+UiJyeHlJQUAFJSUti0aROAud5msxEbG0tpaSmFhYVWxROpVtmFCp/se9VMnuxeRLxQK8cUjh07xoEDB+jatStFRUVEREQAEB4eTlFREQAul4vIyEhzn8jISFwul7mtSG3TMzPkWmT52UffffcdaWlpPPvsszRt2vSSz2w2GzabzeoIIiLiIUtL4cKFC6SlpZGcnMyAAQMAaNmypTktVFhYSFhYGAB2ux2n02nu63Q6sdvtVsYTEZH/YFkpGIbB7373O2666SbGjRtnrk9ISMDhcADgcDhITEy8ZL1hGOTn5xMSEqKpIxGRWmbZMYV//vOfZGVlcfPNNzNs2DAApk6dyoQJE5g8eTIrV66kdevWLFq0CID+/fuzbds2kpKSCA4OZu7cuVZFE7k26O+Q1IBlpdCjRw++vMxtezMyMiqts9lszJw506o4Itee3r19nUD8kG5zIRKoduxwLyJe0G0uRALVs8+6X/U8BfGCRgoiImJSKYiIiEmlICIiJpWCiIiYdKBZJFD9cA2QiDdUCiKBKjbW1wnED2n6SCRQbdrkXkS8oJGCSKB68UX3q57AJl7QSEHEAlf6kJ068ZAeuSZppCBigSt5QA/oIT3iOxopiIiISaUgIiImTR+JBKply3ydQPyQSkEkUMXE+DqB+CFNH4kEqg8/dC8iXrCsFGbMmEF8fDxDhw4115WUlDBu3DgGDBjAuHHjOHXqFOB+nvOLL75IUlISycnJ7Nu3z6pYIteOV15xLyJesKwUhg8fzp///OdL1qWnpxMfH092djbx8fGkp6cDkJuby+HDh8nOzmb27NnMmjXLqlgiIvIzLCuFnj170qxZs0vW5eTkkJKSAkBKSgqbfrgE/8f1NpuN2NhYSktLKSwstCqaiIhcRq0eUygqKiIiIgKA8PBwioqKAHC5XERGRprbRUZG4nK5ajOaiIjgwwPNNpsNm83mqx8vIiJVqNVTUlu2bElhYSEREREUFhYSFhYGgN1ux+l0mts5nU7sdnttRhMJPG+/7esE4odqdaSQkJCAw+EAwOFwkJiYeMl6wzDIz88nJCTEnGYSkRqKinIvIl6wbKQwdepUPv30U4qLi+nXrx+TJk1iwoQJTJ48mZUrV9K6dWsW/fBkqP79+7Nt2zaSkpIIDg5m7ty5VsUSuXa89577ddQo3+YQv2JZKbz66qtVrs/IyKi0zmazMXPmTKuiiFyblixxv6oUxAu6ollEREwqBRERMakUROqgK3nymp7aJldCd0kVqYOu5MltemqbXAmVgkigWrnS1wnED6kURAJM2YUKGjeoD9dfX/N95ZqlUhAJMD9OPY343H3DyZVdbvd4X009iQ40iwSoEZ9vMotBxFPXbCnoDA0Rkcqu2emjKzm7AzTMlsB0JccUdDwiMFyzpSAilelUWLlmp49ERKQyjRREAtTYkbN8HUH8kEpBJECVNWjs6wjihzR9JBKgHti1lgd21fxkCrk2qRREAtTQL7Yz9Ivtvo4hfkalICJ1gu4MWzfUqWMKubm5zJkzh4sXLzJy5EgmTJjg60gi4qErvU5Bp8PWDXWmFCoqKvj973/P8uXLsdvtjBgxgoSEBNq3b+/raCLiAV9eEKqL7q6eOlMKe/bs4cYbbyQqKgqAIUOGkJOTo1IQkWpdSSF9MXtQjX9uIJaRzTAMw9chANavX8/27duZM2cOAA6Hgz179vDCCy9cdp+4uDjatGlTWxFFRKIe77wAAAb8SURBVAJCQUEBO3furPKzOjNSqInL/aFERKRm6szZR3a7HafTab53uVzY7XYfJhIRufbUmVLo0qULhw8f5ujRo5w/f561a9eSkJDg61giIteUOjN9FBQUxAsvvMDDDz9MRUUFd911F9HR0b6OJSJyTakzB5pFRMT36sz0kYiI+J5KQURETAFRCiUlJYwbN44BAwYwbtw4Tp06VWmbAwcOMGrUKIYMGUJycjLr1q3zQdLKPMkO8NBDD9GjRw8eeeSRWk5YWW5uLgMHDiQpKYn09PRKn58/f57JkyeTlJTEyJEjOXbsmA9SVq267J999hmpqan86le/Yv369T5I+POqy798+XIGDx5McnIyY8aMoaCgwAcpq1Zd9r/+9a8kJyczbNgw7r33Xg4dOuSDlFWrLvuPNmzYQExMDJ9//nktprvKjAAwf/58Y9myZYZhGMayZcuMl19+udI2X3/9tfF///d/hmEYhtPpNPr06WOcOnWqNmNWyZPshmEYO3bsMHJycowJEybUZrxKysvLjcTEROPIkSPGuXPnjOTkZOPgwYOXbPPOO+8Yzz//vGEYhrFmzRrjySef9EXUSjzJfvToUePAgQPGtGnTjI8++shHSavmSf68vDzj7NmzhmEYxrvvvutX/+1Pnz5tfr1p0yZj/PjxtR2zSp5kNwx3/vvuu88YOXKksWfPHh8kvToCYqSQk5NDSkoKACkpKWzatKnSNr/85S/5xS9+AbiviQgLC+PkyZO1GbNKnmQHiI+Pp0mTJrUZrUo/vR1Jw4YNzduR/NTmzZtJTU0FYODAgeTl5WHUgfMZPMnetm1bOnToQL16de+vhif5e/XqRXBwMACxsbGXXPvjS55kb9q0qfn1999/j81mq+2YVfIkO8Af//hHfvvb39KoUSMfpLx66t7/+TVQVFREREQEAOHh4RQVFf3s9nv27OHChQvccMMNtRHvZ3mb3ddcLheRkZHme7vdjsvlqrRNq1atAPepxiEhIRQXF9dqzqp4kr0u8zb/ypUr6devX21Eq5an2d99911uv/12/vCHP/Dcc8/VZsTL8iT7vn37cDqd3HbbbbWc7uqrM9cpVGfs2LF8++23ldZPnjz5kvc2m+1nf8MoLCxk2rRpzJ8/v9Z+G7xa2UU8lZWVxd69e3nnnXd8HcUr999/P/fffz8ffvghS5YsYf78+b6OVK2LFy8yb948XnrpJV9HuSr8phTeeuuty37WsmVLCgsLiYiIoLCwkLCwsCq3O3PmDI888ghTpkwhNjbWoqSVXY3sdYUntyOx2+2cOHGCyMhIysvLOX36NC1atKjtqJX4+61UPM2/Y8cOli5dyjvvvEPDhg1rM+JlefvffsiQIcyaNasWklWvuuzfffcdX331FQ8++CAA33zzDY899hhLliyhS5cutZ73SgXE9FFCQgIOhwNw3101MTGx0jbnz59n4sSJDBs2jEGDan6r3KvNk+x1iSe3I0lISCAzMxNwn43Rq1evOjEC8vdbqXiSf//+/bzwwgssWbKEli1b+ihpZZ5kP3z4sPn11q1bufHGG2s5ZdWqyx4SEsLOnTvZvHkzmzdvJjY21m8LAQiMs49OnjxpPPjgg0ZSUpIxZswYo7i42DAMw9izZ4/x7LPPGoZhGA6Hw/jVr35l3Hnnneayf/9+X8Y2DMOz7IZhGPfee68RFxdndOnSxejbt6+Rm5vrq8jG1q1bjQEDBhiJiYnGG2+8YRiGYSxatMjYtGmTYRiGUVZWZkyaNMm4/fbbjbvuuss4cuSIz7L+p+qy/+tf/zL69u1rdO3a1bj11luNwYMH+zJuJdXlHzNmjBEfH2/+P/7II4/4Mu4lqss+e/ZsY/Dgwcadd95pPPDAA8ZXX33ly7iXqC77Tz3wwAN+ffaRbnMhIiKmgJg+EhGRq0OlICIiJpWCiIiYVAoiImJSKYiIiEmlIHIVHDt2jKFDh/o6hsgVUymIXAWGYXDx4kVfxxC5YrpOQaSGjh07xkMPPUTXrl3Jzs4mLCyM3r17s3v3bux2O2+88QaNGzfmwIEDzJw5k++//54bbriBuXPn0qxZM1/HF6mSRgoiV+Df//439913H2vWrMHpdHL//fezdu1aQkJC2LBhAwDPPPMMTz/9NB9++CE333wzr7/+uo9Ti1yeSkHkCrRu3dq8uWLbtm3p2LEjAJ06daKgoIDTp09z+vRpbr31VgBSU1P5xz/+4bO8ItVRKYhcgeuuu878+qd3JK1fvz4VFRW+iCRyRVQKIhYKCQkhNDTUHB1kZWXRs2dPH6cSuTy/eZ6CiL+aP3++eaA5KioqYB7GIoFJZx+JiIhJ00ciImJSKYiIiEmlICIiJpWCiIiYVAoiImJSKYiIiEmlICIipv8Hea+qeHl4v6kAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {}
        }
      ]
    }
  ]
}