Fig5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy \n",
    "from sklearn.model_selection import ParameterGrid\n",
    "from sklearn.manifold import Isomap\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "\n",
    "import librosa\n",
    "from librosa import cqt\n",
    "from librosa.core import amplitude_to_db\n",
    "from librosa.display import specshow\n",
    "\n",
    "import os\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir= '/Users/sripathisridhar/googleDrive/projects/helix/helix2019_data/SOL'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_paths= sorted(glob.glob(os.path.join(data_dir, '**', '*.wav')))\n",
    "\n",
    "file_names= []\n",
    "for file_path in file_paths:\n",
    "    file_names.append(os.path.basename(file_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "hop_size= 512\n",
    "q= 24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1212/1212 [05:17<00:00,  4.49it/s]\n"
     ]
    }
   ],
   "source": [
    "features_dict= {}\n",
    "feature_key= ''\n",
    "\n",
    "for file_path in tqdm(file_paths, disable=False):\n",
    "\n",
    "    # Read audio files\n",
    "    waveform, sample_rate= librosa.load(file_path, sr=None)\n",
    "\n",
    "    # Compute CQTs\n",
    "    cqt_complex= cqt(y= waveform,\n",
    "                    sr= sample_rate,\n",
    "                    hop_length= hop_size,\n",
    "                    bins_per_octave= q,\n",
    "                    n_bins= q*7,\n",
    "                    sparsity=1e-6,\n",
    "                    )\n",
    "    scalogram= np.abs(cqt_complex)**2\n",
    "\n",
    "    # Find frame with maximum RMS value\n",
    "    rms= librosa.feature.rms(y= waveform,\n",
    "                              hop_length= hop_size,\n",
    "                             )\n",
    "    rms_argmax= np.argmax(rms)\n",
    "    feature= scalogram[:,rms_argmax]\n",
    "\n",
    "    # Stack in dict\n",
    "    file_name= os.path.basename(file_path)\n",
    "    feature_key= f'{file_name}'\n",
    "    features_dict[feature_key]= feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = {\n",
    "    'Q': [24],\n",
    "    'k': [3],\n",
    "    'comp': ['none','cuberoot'],\n",
    "    'instr': ['all'],\n",
    "    'dyn': ['all']\n",
    "}\n",
    "\n",
    "settings = list(ParameterGrid(grid))\n",
    "\n",
    "for setting in settings:\n",
    "    \n",
    "    if setting[\"instr\"] == 'all':\n",
    "        setting['instr'] = ''\n",
    "    \n",
    "    if setting['dyn'] == 'all':\n",
    "        setting['dyn'] = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fig 5a: Linear loudness mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_str= []\n",
    "q= 0\n",
    "CQT_OCTAVES= 7\n",
    "\n",
    "features_keys= list(features_dict.keys())\n",
    "\n",
    "setting= settings[0]\n",
    "    \n",
    "q= setting['Q']\n",
    "# Batch process and store in a folder\n",
    "batch_str= [setting['instr'], setting['dyn']]\n",
    "#     print(batch_str)\n",
    "\n",
    "batch_features= []\n",
    "for feature_key in features_keys:\n",
    "    # Get features that match setting\n",
    "\n",
    "    if all(x in feature_key for x in batch_str):\n",
    "        batch_features.append(features_dict[feature_key])\n",
    "\n",
    "batch_features= np.stack(batch_features, axis=1)\n",
    "# print(batch_features.shape)\n",
    "\n",
    "# Isomap parameters\n",
    "hop_size= 512\n",
    "compression= setting['comp']\n",
    "if compression== 'none':\n",
    "    features= batch_features\n",
    "elif compression== 'cuberoot':\n",
    "    features= np.power(batch_features, 1/3.0)\n",
    "n_neighbors= setting['k']\n",
    "n_dimensions= 3\n",
    "n_octaves= 3    \n",
    "\n",
    "# Prune feature matrix\n",
    "bin_low = np.where((np.std(features, axis=1) / np.std(features)) > 0.1)[0][0] + q\n",
    "bin_high = bin_low + n_octaves*q \n",
    "X = features[bin_low:bin_high, :]\n",
    "\n",
    "# Z-score Standardization- improves contrast in correlation matrix\n",
    "mus = np.mean(X, axis=1)\n",
    "sigmas = np.std(X, axis=1)\n",
    "X_std = (X - mus[:, np.newaxis]) / (1e-6 + sigmas[:, np.newaxis]) # 1e-6 to avoid runtime division by zero\n",
    "\n",
    "# Pearson correlation matrix\n",
    "rho_std = np.dot(X_std, X_std.T) / X_std.shape[1]\n",
    "\n",
    "# Isomap embedding\n",
    "isomap= Isomap(n_components= n_dimensions, n_neighbors= n_neighbors)\n",
    "coords = isomap.fit_transform(rho_std)\n",
    "\n",
    "# Get note value\n",
    "freqs= librosa.cqt_frequencies(q*CQT_OCTAVES, fmin=librosa.note_to_hz('C1'), bins_per_octave=q) #librosa CQT default fmin is C1\n",
    "chroma_list= librosa.core.hz_to_note(freqs[bin_low:bin_high])\n",
    "\n",
    "notes=[]\n",
    "reps= q//12\n",
    "for chroma in chroma_list:\n",
    "    for i in range(reps):\n",
    "        notes.append(chroma)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAACjCAYAAACZvbxxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de3hU9bnvv++6zC0zuZGQBDIkASEJECYRGERAQa1oWytq7TlHqj229mBbbU/VHd1td8/edrfHprh3W3x65NTn0dbT3b2t99qn21sVb0AAIWKAyC1kgEAScp1kbmut9/wxa4ZJSEIus5JJWJ9/8jxrZtb8Zuab9/f+3t/7vj9iZpiYGIkw2QMwmf6YIjMxHFNkJoZjiszEcEyRmRiOKTITwzFFZmI4pshMDMcUmYnhmCIzMRxTZCaGI032AC6G+sQ2N4BqAF4AtQBqxPuu9k3uqExGA6XyBrkusDoATgAygAgAPwCPKbSpQ0paMvWVZ6PWK1+6HUFbBnrSBagSEBWaE1HLdv9kjtFk5KScyHSBRa2XrMgBtR0fd+yF3JUHkWWU5xfJDovNO5p7dn6l54IpN/M515S0hA1thRd8ltKckyn9WVJOZIh+gU4A8v5GH0JKBFeWXwbqdUE5l44DZ04ox8+dbryN1i7Xn9/AzN0DbxJ+8LgbQDUzVsuZ1jKlxyGyKsoAqgBs7PxKj2eqCU0XWKL7UAVgY0NboSeVhZaKIvMys/zRocMonV2AnHRX9KolDEmUIktmz/MvmT3vIWb2EZEAYAERlQFgAATgaOiBYw7oPwYRZMkZgpQWRvBMJnShTdUpN/4PGA4zLBaaEp8l5UQWVpRd79c3XL6y7DLJYbVGLzI0RKRWAH9CwuqSmTUAh2KvJSICMPfVI2/8rCAtz0UEiUDRWzAQ6CMsEJbBJbtkRKebqYYXgOxrVHD0MwXlHglSJstBle9uPJELADXri1pTzqKllMiIyLFmUen7bz76yB2SKJ5fURL8cPYtFzcOv6Lk6FL5aPjB48UY5LOpIQm1DUchkaQsylxca8iHMJZaf49W1dmhyauvs+KN93pxWYUEi43SAGwCsPH1E7meVBNaygRjiSgHQNX79Q3PSaLoAbAVUcd2KwCPePOdo/niahENd8RhBrSwBE9WZWSWfZb/5m03vUtEWUn7ABNAMMg1dXvCwYoqS6QroqB8mYRPd0egh6ESp86UIiUsGRHNAZDDzB/ql3wYn49RA2AjdP+FGREwNMVvPwjggyxrds22s+/6iGgJERUz897xfQLj6O7wxleTe9+wN5V6rJ8Hwv8lqPLdRJS2eLmMI/Uq5i+Oh3hSzg2YdJHpTjuY+eNk3dPyeIkv/OBxD/Qfhwi1INRk/FtmP2vIzJ8QkYuI1gA4yMxtB9ZHLggRLHxdnpTpRxdYHQDnkSNBeYXXWTV7tuVaAJ5GPg4AmwQZcmdQxak+wCqQ5pKFA5Mx1uGY1Ig/EV0OoI2ZmyZtEOfHUr429wtzn6h84Vki4YIdhskQWneHdwuivpb80fYeXLnSBX1MW7d3H69RNN5/qk/NOF6voGhR3F50Aai4d1FHyvhlk2LJ9FXgaujWYzLGMBBmPrjzmrb793R+6CpyzJdyrfnAgB2GV69RLrByX/qbZOSP6dXHAEmk2DUZgHd9Uavv/9ZnvaQBdxFBcAgSyp1ZyLHY0oOa8lp3h/eL6Vm1KSG0CRGZuuOn8R8nFFZ2raks+uD9fSdeZea+iXj/keKSM5Yuy1ojNfZ+huZAE2TBgvnOxbIkSF5dYBcEQl+9RvEYKLRa/X3kBaU2HDoUQFmZPaJfhwYsBCBYScKNM+dAIgEiEWnMFQDquju8nlQQmuEi0wVWB8DZGwjL2z9tuvzNJ75xhyQKHwJIKZFB/1GL0xbIABDWQmjwf6I0B5oam1XaPN+1zJljK5T1505EIDS+gMnMkOSGhoCC6PRdkzjeErtL1gUGABCiM0XKBGknIoRRDcDZ2tEr7/usGdctv0ySRCEll9qI/nh+6OEPi2CNLEq/vOe6vA0PrczdUAxA3tf+FkJqIPZ8Q1dzuhWKh3Py8yzP3P+94xsSrFMNAH+mbGGRCNt39KCjU5mQsY2GiRCZF4B82HcOqzxFsWsp8wUkojv3F8To9Ou1ObbCiCfrWhzzxyMe8anLKNKzan3pWbX3p2fVrqiorPvm759ti8f/dOfe4w8phyOKxjabgAMH4v8Aho9tpEyET1bb0ROompFhlxOupcwXMBBdUINNMTUANhKR0yamyQGlJ2KXXIlT10TRSEQlzNEYxr2LOnwP/LbgzuUL1NcLCuS09nOKXH+gT1m00DEZYxuUibBkNYcaWwPz5+TE/gNjYYGU+AJGiu7cewBsLU6rqN3T/p+vADDS6R8UZm4GkJ947V9/deaTH/2D78szc+WtnZ3qpz5f+IU33+q8IhWcfmCC4mSfX1V6058fv+t6JKZQX/HDlPgCxgoRLQbQNFia0QS8txsAmNmXcG3podbZLW//NbB53Xpb8V9eCvhuus3x/VRIATJ8uiSiAgB14hU//LPR7zWRMPOnRLQSwPZJeG/fDWsLb1Ya747H7L79rdm7ALw5M190CQJJK9dYq442RK4DCismW2gT4ZO5mTkl/a/xsjxno/3RZSe2ZlpmV0K30N/9yPjpU2m82/3Lf1zxu+aWPkfBTIcMoOqLN2V8E9wlkBD9TXNmivLJJjUt0Kc9AuA7Ro9pOFImC2Oq8esrFfdX5z/zfEeo6euIWpNNAOp+faXinoC3r76sON1x3OePx+wUsc8Cgpzo/XiWytLenZHPKY13b1Ea796p/52I8fXDUJER0VwAx418j0mkGoDTIWVL54KNwMSm2ngByDlZVrS1BwEAFs1JYPRzsAkUqSiZ7T7W1LMJCf8IEy00oy1ZDjO3Gvwek4UXgJxnL4VAIo73bEdL4LOJiv/VAogsmJuBw43RdYczMkshCGHWEIvGRgiCtii3XDzTGhhsl2LCMMwn0zfBpzPxfcUsqxtZVjf8kbbIvnMvnCD6r5UA6vRMXTxUEbpgY33zfut4fLf4dlNWhkU+1xGKCKoc6Dwy8wvdXY3fBjAPQK27b+VqmW2VDnsvFEWDJAnAJATCDQthENFCAD5m7jHkDSYZ3fcatPD4e9vlcwCWAFDvWdDQmmEp+niw541HaPqUVw3A+/xfGpv++Mqxp199s8kH4DAzB/XnbAGw6dSZXlkUCfm5Dujvv1UqfnrC9jSNFJl3uq4qY+hC62ehEleXRCTdPOe5Pzikmbfm2iokWUiLPRQBsHXzfmtSfuiNN1Ree+2yefeVzM6af1VVyTvQ45C6EOtCYdXZcLRLXlKeHRe4VPz0hIU1DBEZEYkAqph5d9JvPsV4qCK0k5m9vt5tmONcm/hQ7eb91hXjvb+646fugNi9v+bPL6cXLorQqooSrah3WY9dTa9IEFr1+7Vn1q7x5r8LoGYiBQYY55MtAfCJQfeeatQSUZVVTDdk77ZHan10z4znM6QiP5raNRQ7eoUz9oaMpee+/GgmcLcuqPuvKaGlzLwnGe85WoxaXcrMHDbo3lONGgB+Zi2+6kMS925PpO25UaUI0rMBfzeDSYNKEZxI23PjgKdO2kIs6SIjIhuAULLvO1XRnXuPX2l+EQnpQ+NcXcbpspwGk4ZwkGG1Ra8xaei0nE7G7ZOCEdPlYgBJqzyaDmzeb/UR3V5jxHTVh+B7zLi9tVlDXqEIANAYCCD4XrLfa6wYMV0KevsAkwlgV6e/V2VGIMCQZEBjhsqMXZ3+3skeW4ykWjIiykS0JMtkguhVeGFtawhdYQ1yBrD3YBihHCCkYeFkjy1GUiyZ0ni3W2m8e8ub/7b+g9DRr903GZuwUwCjHO/atmaOBDMAX5qG490qWk5pKZV5PG6RxQJ+ADY57NIiQaBJ2YS9hKnpaeOAI4sUAJhRLEYC3RzY/kzo6QHP48na6kuGJasG4Gxu6ZPzcuxACjf+mGQM2Vp5bEW3b8czoTv8bRxfveaXiovfezIsDWgo0wkg04gxXIxk+GReAPLJ5j4s9+TErqVkNdJkQUQWDOgylEyOfKC2HPkg8IvEHZbHGD4iWk1EH+tF1C0AZgHoMGocQ5EMS1YLIDLAEKeUT5ACuABMeC0AgA8BXEFEMqIBYNckjCEpIosWxDIMiWhPE9IBGJKNou8TZwJoH/iYnmr0LoCrjHjvkTLu6VIqftr38lPXrnQ5Lf8AYD70bISJ3oRNcVwATiXzhrE29Pt+9bM1+46fOFrunjWoD8zMGhG9D2AtDBL6xUhKFobeAqqOmdXxD2n6QZTczenENvS7jxyTmaEsnz+3B8N0pCSiNAD3MPOvkjWOkZKsiL9gCmxCiXfBBgAiSLjIip6ZewEcJKJlEzLCBCa90+J05xa3371ixn/83S1ufwl0V+Iln3O8rkS8bxkA6B2+R7KibwdwiogqmHn/OMcwYsySOAO5xe13A6izifm3IaFaSL8+HuKNl+n8RsKIVvR6m4MeIlowzjGMmHGLjIhmYJCVjQkAfVojUGzGSFagOt7iisHQmAf2LRsWZm5EtNaneJzjGBFjFll3h9fd3eHd8pc/l23raFv2fb2Jrkl/vADkAeWQ4w5U6869B8DWvmC4vrm980WMsg09MzcAcOltJAxlTD5ZYldmh0OQBYEWANiYKu0jU4h42VzCtaQEqnVB3X/NBloOQEtsvjIMQSKyM3MAAJh5PxEtI6IQMxs2G43V8e+3uoF5ROBQ6PWR7EL0uzYiUJ0HYN/FnqTu+Kn7s+cf+F5EUZepO376AfSKJmbeTUSrVv7vtNY8r3Q/EiqvXr62KykGY6zTZXx1kxBnM/crB6CvIj1hrf15jUO7oKdeJ2F1mUgBgObhnhDr2zsnP/Nr/kC4EvoCRL+Om99KbyIRuzWF+7Uz2PB2RlJcoLFaMsOmgenGSz6nj+jW7wOQmPlksu77ZH2WW1W5+pYHbF+4/uu2JU/WZ9UM07u/GoBTEoVYQ5aBM0917uWirW2vKudUiiAhuafPjdWS9WvgC3O/8mK0ApiZrJs9WZ/lBlAXDmBTZp5QAt3y6NcHo19cTSdx5vESkTxjiYi2T9TBHh8XYxJZYlfmQIDroU8DptM/OPpuSDITBqsBOEO9LNvT4oHY4UIjtQAi3b1BuByW2LXEmacWQESQCFmlItrrFSiaGGntzbYtfbF859IXy7csfbF8zFPnmCP+uqDu/8JNtAzAHjaiFH16kUyReQHIAT/Dmtb/pJIhnl8DYOPx0x2usqLcwRYg8QYukp1k2yw5cnSfS3LMs5Yj4WCMpS+We/bcenDUhiQZEf8zGNAo18RwagFEwkGG1R6/NqRPrPfn9Zxo7nzRapHOH++o9+3VV5Hx1vK9lqyDUoakBHzhpLScSsbe5SkAS3GRFY5JUolaHoYLNLLQiHjFD323rqRfDNWfRBfa/QCw9MXynZYcyP6GAJgZemnAmH20cVsyc5qceGKHRHS18gvhAD6FbpmSeDJcLYCINV9GuDWeNT7m6IGZhTFFuXdRh48W02ZE/8+TXZleA2CjnCE5e48EZOvM8UUPzCwMkwvQnXsPgK1a8Hz0YCxOP5C8zNhlZi+y4SGiKgD7OZoxMW4e2ZnuPvK+8gurE4vcVdK7AGoeW9E9pAiIyA5gLjPXj+Z9iGg5M+8az1iTJbIKJLSRNOmPZ16rOxhu+CeB7JUWec6HAGrqjuaO2X96ZGe6G0Bd80HVRQQpv0yMd1AcSmh6e9UTeobsiCGiEgD+8TSYTtZ0eQzA3CTda1rhmdfqBlBnled9VdW6qqBH5/XrYyWap0aQBtkmGgrHaAWm0wigeAyvi5MUkemDdyTjXtMQXRCSzFCB5CQuRvPU+k9ChiQoJCN6YDr+xnOxfcNR44i4DhCTFuof8U/ZBAVTZMYTz8dPYMyCUJ/Y5v72gV/eYlUdQk8zkFFAEDQRxNQH4xIUevWSujFhisx4kp2xUp0Vnum4v34L3O2LUBRagOWt67Xv1P/qpeFWl+OkAUDZWF+c7CZ4ZO4A9KfuaK7PM6/VA6Ba0/rWAXgH41tdegHIGaEcVLSvwb0HbwaixmLIpnfjbfjCzCoRjdkgJVNkLYjmTJ1N4j2nBbqg7ifakIx4Yi2AqrM97XJ+enbs2sWm3yIAJ8b5vmMmmdPlSQCFSbzfdETTG6SMhxoAfl9HqzI7MxcY2fSbwcyd43zfiN4daNQkzZLpjT2m+6Fd4+UsooUfY+5/Lt53tU99Ypunuattc6bduQj69Cved7WhCaOVP7qqJ+IP/27NUzfPg15o8v49r4zoPZN67I25vTQ8uhXzMPO4W9AT0WoA3cx80ZNfxvu7rHnqZjeAup7jHS5XSVZiapFnJEJL2nTZ0Fbo/s2z2Q81tBXubGgr3NLQVmgW+w5AT8NO1ndejIk7sLY63BNyqkFlTJXwSZkudUHV5eaJsfrCKgAbG9oKPZN9yPp0RF/ppSXjmMfrH3v4gpPu3njk5z79fXIAzF384JXrbDPscmZ5buJLRxxQTtZ/VTUApyyTFA5zbABmc2LjKMEIrJj6yrPuwAtPP/Huz374tPrKs1v0vmZxdIHVIbqf6gWwKdzb90lu+WVfJCIvoguG2syynHdsuWljDignS2ReALK7WISvMZ7JYhb7GkcWLtLkJtYo79S5jv9ROrtgMWIFvf2FVg3Aqamq3NV0Gl1Np2UlEEqr+tpt65m5lpmP6s8bV0A5WSKrBRBJcwro64svJFJ2L226oz6xzY2A7TUwMlu7uuWcdBcwyOzCzN6uk81yz+kWpBcWIGPOLDhysi4wDrpzHy800f+OyOkHkhfCiJdUof8RyGax7wSjPrEtOgUKWuaJ1jayyjISIktxARHRoor/9qWm3PLLqkRZumgnAF1QY6omT4ol0537aLFv7/l0XdPp70/pugPu2Ut++8D8q/bUlq47sKV03YGxrsAFDH34RDUAZ9PpLgqFI/CUzEl8LFL72dFDRLQKwKn8JWUPiLJkeCeApB8PTUT5iB6qagosAV1QdWqky6lGOmWLoygea2p4Z+GIvysisgK4FsDHzHxm4OOBX7+z1yIIlbt8n2L5FRkAMUBA87lOPnLmrH9b/cF1/+sPL8QLT4ZbXSaLpFcrMfMZIloOwBRZf6oBOAUpXY4EmoCxt9uaj+jBExfsEV/94FPu/zkjt+xGVwagScCZAnQKzahvOcozHBkNaxaVXr/2B//c73fRBWVouy+zJG7i8AKQB+y8jWUFbgMQHCLbpfoPnR3iujQXeiNhfHT4ENJsDizNvypsE4TrxZuN3XoaCqPyycx0nwtJavLiEHhbVUX+St12PH+iAZkF83AsLRv3nGo6aPTe5nAYZck6iSjbyBaRU5D4CpwEi6ypoYggWpPtZNf6W32Xt/s7pT53Oe47cxKICvmDJL7HqDHKkh0FMM+ge09JdOfeA2CrxVGyLxLed9ya33zC7j5ZXXnXGyNaZW54O8PtfdT+d6sedzyz4e2MLQM7IX768pb/lO0uv2x3KqJsBVIklGSIyMzs2MFpeGehz+4+WeMoaimSs3rmCrJayRq+wyod9lxTN6xvtvrR+7yhdvmws1C+XbTSIubzLTcpypr2Y598ZHVmLgn5O+NnXwLwbHv8nkldhJmO/8RTzQwnUbS3Pwkg1tgKUdnmmde6YLC07NJ1B9z5N/Zs62uJWB15GgSJQASZGU7WuBrASwC2612tO4ho83irvpOJkYUkAb003qQ/XiLIrClg1tDX1oBg13GEcNgaDB/6JyLKJiL7gATQaltBgxWkwn9Sg31m9CEtwvLZncp1AN6LtU3XMzRSaiYxUmSHAJQbeP8px9UPPuWWLH02JdgJ0eoCQNC0CCzO2ZDCBSRQWiWiZ1fOBXA5ES0loqV9nbvXttTm0YnXFIQ6GYJEiPQy2uuZ86+Q3hrQX+MyAEcm4/MNhWHTJTMrSchnnzZc/eBTbgB1dmeHs+vUOaTlVgAgyJZMaKEgpJaqiChbPmTmYwNfW7ruwLtddf6yvGs+lGy5EYQ6Nfh9ArLLXdRVf92z+Fy/p2cw82cT9LFGhFl3OXFUA3AqoU45I0+FxeaHIIZgtxcitEdUELYMuQo89cm33iCtQm195yc4t/sKdB4oBgVvw8mXH490frzxzon9GKPHaMdfIyLRPAsTgB7xD3WfgzOvCETR8+ZV0dXbqxX+x6m2h58uf/wc1jyFLQxttSr3CRFHq9Zee+5w9m3qT+QzTYfVgHtx5OTP0NWmQHL5IUjqlMjZM9qSHYTpl8U4f3zgeZ8+Imb0PN3gW/UNtfC1XAB1DN5EECrFSNqS8F5XZfrCjA25N2RtC2u7jzvy5Yi94BwsWV0QJBUYsGOQqtVihoqMmfsAmCvMKHp2adxJ7xcoXfCNy9f7m7pcBJLD3SF0f9YOR346bFKOLLCQsaK8abWdBA3Dp+UUI9rqKaUwVGSvn8h1/3ir66HXT+TufP1E7pbXT+ReshVMekDUE/J3DRootaRbq8JdQamroQ1qUEHGghmQHRYIig0aaUJXWlvWc9lFwi22dFUG7Yu9ft/vr0+Mq80YT7M6o0h6PlkMXVB1DfsirgUeSSKieP7U+qLWSzYNaLAaSCKism8t/3dbjuNWLaJK6fOi7QcYDMXaCbZ3YOmOSvySHwWiFmyr5fGSC9JzUrXu1UjHvxqAs6BIlJqbNMwqEi/54woHC5QSUREA9z0tP/73+WfLb+06cxbyAQdeWfl7+HI/A1u7wT7Gt9X/Hpt3poSzn4iRIvMCkNOzBBw/GEbBHAF6L4Up9QUli+sfe9i94r6v/cTidFRc/9jDH3WfOvPLnU/8fhaApg+uDYfTu/GCABA552BP84f4+789hsdW/yAiqEuk9Yfm04LCeLfUQdODdKf/knP846up8stl1O9WgEu0gilW38iattGa7rzcf6b1Xk1RP77m0QcamfmEhfBbQV909im9yLXmoz10Bne89VO5+6lN9KePPGjrE4DhsyoKMI4eG0ZipMjitXoWG6Fwrhg58qkSwKVZwaQHYsNSZ+NJWNOdUmbRbLsoS9UAIBHmEQEd4Ta0hc9gjmMefIFjmGkhWORG5NnK+P/VOVuhO/uWx0sG82ln41ITme7cx2v1MmcIW9/4U/DzNxS3jan90BTHC0DOLJqNrBI3ZIcdSPCtFMZRZuBkoBFzHNE0vLmOMuzr/BRWUUS6VaAz3WnZP3hjVs0QAgOQuilWhkb8daHFnfz1zwD0O/IQUW4qLrUNpBbAkPWNYcY3bYQdzFGfihlIt8zAa2f/hJL0jdF5lJkwRRdNE753ycx1ZSulNf+yI/M3T9Zn7XyyPmvLMCfOTheGLfNf9ZZc+1brm9cLQtqhoIbAyRDj1VYNaZYKtAb3A4ivTKfkosmwONlQPFmf5WbmukM7FFfZFf3jZ0k85SzluFh9IxEt77i9+wyAupdOiRnb204KRBL6lLPItpZDFuwRAFs377cOFh/LBeBk5olqJTUqJkNkWwBsUiIsH92roNQrA3qA8d5FHVNuKkgWROTtuL37TgCbuiKQ/3H/x8i1r4DKQHPfdm1W2pVdADyb91vjwgw/eNwNoHqbb8e6NYXedwQShvXZJovJSPXxApAlmVC8WEJDbQTMfMnGzwCAiLIBdCDW2VoGqstLcXmWhtl2RpnL1tHQ9fyqQQRWB2BTmmxfJJCwCUCdfj2lmAyRxeNnVgdhrkfCoR2Koqp8ycXPEriMmQ8j4bspsKfhhnwVXy9RIt8tXfzH13xfzRrwmmqcb3ADpHBPuMkQWT8nWLZSZP4yqeeXX/d/aGbSDrtAOKd3PozhBSC39LUh05oRu5aSM8KEiyx2tDESel1JMnmO7FH/BGDdWNt4T1USpkpkPucatA9Y5nMuHzM3AFiQ8NJaAJFjnU24LKs4di0ld1Qm3PEfDn2ZvhbAR5fK2ZlE5GUematwVcF3Vq3J+9bX8xxli7Md4QMrS/bdqrJiL85wJ/aEG2pHYNJIKZEB8Y3eqwHsGuP5jFOKkYrs11cqbgB1x7o/cs1Nv1ICEGkOvBf6+zXznrfJ2kLoYZFUExiQgsW9+tbIu0S0moj2M3PXZI9ptOgHpvaLiQ1WtEtEWQBGelJINQBntrVIag81AWDZLszFM7tn+7/7kbQiWWM3gpQTWQxm/oCIVt7wiK2n8hZLrDtzLS5y1vZkEzupF+dXflUANnrmtXoGEdp8ACOt9PY6RE32zszF6b7dOBdSQfIquU9NPUd/ICldEvfwDtfJ3Pnih73tWrwFOIA6/QzuVGVUoYWRbmpnydqBm/JCmO9U0RU+jtW5c3BTXghZsnYgKaM2kJQWGYDq2YtFe08Lyz2tGpDCsaAEvABkjUPoC+6NXbsgtEBEmQBG7Ap8LicEmQCRgNZgCwrsMyETcM2M8JgPO50oUl1kXgByfpkIJQic2q+i67SWkrGgBGoBRILhQ4nXBgstLAAw4kpvq4BKQc97tYhWAIBAADNuffUaJZUte8qLLB4Bz3ILmF0hQrQgcujtSBMRVaRonWE0oMqKAho+m3U0+V8hRRAGPpsZCCuCgNS27Knr+OtccD6AM0fwl10rPPDKD4OdAK4gIgXA3gFNRwxluIyKuqO5vvmFby1XtI5/FihtMYB3MWB1OdqpEgC6gpKWmxYGACzNXgrmqMh6QhIhtS176sXJBqI7+f1+0MTVpd7nqwrRf5h6Zu4Gol0JB77u5Wu7xr0qTTiPaODBGJ6Eg6+WAdgLoJKZ9wy8h35u0a7RWLID6yNbBPC3nTZFsIgawqoAf1BCu0pam0a/+dLfpJTNYEl1SwZdUEN+gboF26VPnQuJKL30Lmuk7C7rGxgQRtjwdoZnMKHp5zn2E+RQR7poqvpwX2u7S40okiXNDnt25mClfoJ+bvcwwx71f3eNynRnV0DOiCbKAhqATg09SPG6iZQX2UjRf7R6AFj7f5zPtNcrLtlJkqtIBPRVaViVHl36YrkfCWJytC8ABolrrXnqZs/797zi0/cW4/1vl2264+r0WTMl0WJB96kzsR4MicfJWACEhxonEWUAGPURggtfl30H1kcqVMajGuPGAAOdKv6qgMa0PQ8AAAPXSURBVH78pb9JKRs3BKaRyBLJXCCWA5DCPYz2g1FXzZIryZ1yzp2IGoC4mJiUV4glJwCZNUbgrF9Wg4or3BncTESbET2NbXfM8lz/2MPvAijF+TgY0H/1WAFg3y1uv3tVzp8fuGn22fmSkLYTQM1LPqdPf+2YWm0ufF32Abh7LK+dTKalyKAXblhcJGeXRz/i6RMOzd8SJBDJFI0FyMzsinS33ySSRQaiU5A93wnJJkkAiofou3qxw8rEDYU9swDUZVg8zpDWIktCSaVVxJ2bV0Reum/hm8vnZ6x959dXKjXf/Si1LVCySPUQxli5IC9LyHZw2gKHkDbfDsc8GxzzbEi7zC6lLbAqzqLMiLMoE67iTEg2CRgmZUZ37j0AtqrhyKfQ03HeeOTnPn2qjCB2yCxlymGtA1YR8vJ8a4aqdd/lkvMWQd+50De9pz3T0pK9fG2Xb8PbGR4kOPMqi05ELVC/aS5i6/irFM64GaM4RjF2HhERLR2weqxA1L/rd8TNHJcEkYDTwQZhjnM5MPZzlaYk01JkQFRoSPgBl75Y7gZwgZhYCv8YwI8xwtXlRRCYWbnF7a9F1OeTASDdIkAgAkMb9OzJ6c60FdlA9tx60Lf0xfJ+1g1AzZ5bD8bENC6Lomf0xgLC/fy27rCmhbRGIdtakviSlMxiNYKUD8amMkS0dMHa+hYA1YHOPWut6RXbBMHy84Z3FvpucfvdAKrbQh+sK3Ku2TXDsfvL89JXWTEggHspOP+XjCVLJnq79OrFt/zoc+FOtRgsCBBkWRAspQDuKF13wNPgW+h7qCJUw93+mVm2vsUtgcDbxU7uEAWKZ7FeCgIDTJGNmlg/fgDOUM8smZjAmgL9FJu4Q/9QRagGQJ1TKnCdCx2UZtnXek71an4MKNC9FJiuIQwjiSclKgEHAAGh3sOwps2PPR5z6KPPi52hFPXZUj0XzhBMSzZ6vNBXjYLcA4DBrICEeGQk5tDr1eDFEMkSe+ySWVEmYlqy0VMLIKKpIshG0LTemLFCtB0s9yG6uqwFELFLM2ARXbHXXjIrykRMkY2eGgD+UF+65phRgr7AH+HIKwDJYYhpvZo1r+Ul/QDVYdtFXUqYIhslsX78StjWSoKEjMIqWGf0wpbfAkt2pyBYIgsBQHfuL6gGv9ScfsCMk42Zyrve2ILoHuTAbIyt+35//bTfKhoNpiUbO+Z0OEJMSzYO9APq+21TDTiGxgSmyEwmAHO6NDEcU2QmhmOKzMRwTJGZGI4pMhPDMUVmYjimyEwMxxSZieGYIjMxHFNkJoZjiszEcEyRmRjO/weTbOBWl2zbdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 396x198 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAACnCAYAAACB6MdcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO1de5QdRZn/Vc/NhDAkMxOGJEDIAyNPiUjkJQRRkUWUI+ruYkQEFBBZ36virhr3ILCi+MCDiyAKyOKKiIARFI4iBFCQhxAhhFdCQghkMsxkkpDJzNx7a/+o+qqrvqq6fSfEuXNz6nfOnDt9u7uqum/3V7/6nkJKiYSEsY6s0QNISKgH6UFNaAqkBzWhKZAe1ISmQHpQE5oC6UFNaAqURrOzrq52OXPmNABboxIT+lNGtkdy7lb0LmP9xrYBKbZuzELS+e7/7jjcPkS54n5fKuVd6nNi12DakOFt3UiNfX5btVBe3gsAeHqj2t5SeYV29UgpdwmdM6oP6syZ0/DA/ZcBsuruqOrtLHP/t/e16KHSuULvr5T944vOrbL+OexxaIhyoJ9a2wDkuNbwmGPXT33pbZllzv/OOKgP/UCK/n73+8k75+MZHqp5DaYNvZ9v2/cjuC/QVhD62HUf/D8AwDvvUm080XedPqCyMnbqqD6oBoIxjsz6PmPH0La5IaytjLeV+efyB4Ofw8eVZQDCDwTYA2MeIDous9rgY+YvCr1AfDihF8lu1+53oxJLYv16tX/CBLWf+hou5+dEriG2337ojJwM7HP2x8ZdrZr7u2LdZADAE33fAwDs33my3v5ZrJXEUROaA42RqHwqtL83U6j+zkiXsrtdawr2ztXdsWnUgLaJRsjMl5hsWqS2jSQdv4Pa3qCn4ExATtjRPTcmybkEtdvmUo/aIFpBErSrS22PHx9u0wKnD+Z7uj/sOFkq5WMq1/c78LZ+cfh9+PrKJQCA5X2LAAD3Hfk5AMDsXRRn3e2m6JCTRE1oDoyuRBXQxDy2vwZHjW23uNwNwmqcC+xCbho4toi/cYliJEtLtC2XCwPItHSs+NJT8H6I1+pjSZKjtZX1bx1P7WopLCtMooow7zTbLaVcT1HnYpLf68f7W4wk3bPzBADAIXe+S3ekzx3/U8SQJGpCU2B0JaqEevPYGy1s1QbjqDFu6L1iFV/15HHSmETV43FWwfQ/9cu36ZLGuZJMdnSabaHVQp60MW1EpDGpk1pKuTqou1vta2tTbehPMbjFvTauJmqBz/NJKtP4wPhnaJaISVACbet+N5z9cwDAp26fDQC4rvsCnDt7IQDgDe1a58slew00ZjHFHhgzrdhTP+3jD1cLmwpDF0sPBO9Pn0sLHjmpXX1fHXLPs6Zcs5jiYy65U7C5hMH84fTogfcjuwsTerFsaiAGBtRXr76q+tWLJTFEDxmcccKart2+rHbpAaWXbJi1xfWqWeap30TkwaX99yzbA4B6QAHg5ClfwfnL3uKOKfbQB5Cm/oSmwChLVMlUUFyVUfYtMUzNkaue1EdIjWQkAp3LpIvcaaLapimWxqPbElmW0wF+bkDaANYUbCxmFU8J7o2LSyf6ntRaA5vNmOSuu6rPtp3Uvk0b9DkyH7PVpiOlxkUWWkySEjz1VcBSx+nK0x+4HQBw9oOqr8X9lwAAbjn4XADA/H1WAJjvXn+w9zCSRE1oCoyyRBVBFZStFuG8UjJ+mSvlIwsly3QpOZ+lbrX0o3ONpHOOivBL3abU/NLjlRYv5RK0ECRhyWgwMGAU+tSuGNistvkiLmLaxbjW/J6ZRRSTmKS2ikleazHlmUr1uBatVr4ki/vPAwAc1f4ZAMDxi9+qzqtWzW8bVRPWQJKoCU2Bxij8SZUUUguRhOAKbOFyItNkTKVitQ/JnDBsLgpLohNnE5ZUHucq4w0XHR7W/WoZQ252xvFE5tKHxsSuyUg/Mr/2rFP7aTwdHTlffXWTc918FjC8jwwAtmpMsHvDV/tc4jIOK8plJd0BoHUcAGDwq78EAJx9w+sAANeuVav7i/ZRKqgTpqtrcbQTEWeYtOpP2G7QWDc/ek1sycYdVca53FSw1bfHZW03P4SlMElrjzu2+PzSSEjaZ6SBlmUknUTAQMC5MR8HSRjNO0Wvcs6QM2eacZvVPfVDiEkjkoIkAcePD3DCiFmWQPd6yNIvM/3svQ8qPSlJ0lOmfgUA8O8PH6KO4wYS20+YuyqiGEmiJjQFxoQJ1cCWptzyxM2MnJMFuKnHhYibclc93hbgOTmLoS3uWLX+0kiOihsKgkolcA1akuzAxkwSZvfdWRvlnPu2sPHQMLnbH9cV207j3FoVk6T8mioVXPu2RwAAP3peSf0H+i8DANx+6BcBAEce/Jw6tnyw26b9u3iOO4mjJmxnaAxHJbCVtRjc4q+yuVMK4zlmFR6SipwL8Xgezm9lvho3K/Eyk6ScGxPfZbZ3SIuTVcN2eDGg+KdZ2TOnFLTknN3rh6QeSWdPV2rxSx5fxnXSZtZwfR4MD28p4eZVSrqSJD20/RMAgHf8Yb7bFvORIIgAR42FtYSQJGpCU6Ahtn7P1kveSwMDUa8cz+JiLEQ+PNs9d9Wj0ZA0YlxNtpRyC5CRcsPsUohXMadnwtBwLg1JYmrpT20b6W8k3Ob8XGqauDDzAzD9cS0J55+2ZEeY9xtOOjio+1W60vIFNwAArr9jT9zS900AwHs7vwwAOHGGdgyn+1J2pbW3krddBWNhLTWQJGpCU6Ahtn7DQ7kv48SJxcF7nF8RzAq+7FltohyVRsUknLC9iiIhKCQNPE8jGl8rPPu86OlRbXR2Om0bnm10shanpX7oPpAfAPmnTpnijMPTAsDyfGK6aO/eaklK51540+sBAOc9ex6ueaPSk55876HqGJKk48Y5n2aW4FoT697wWbAePWpjF1ME66K8xROBmwFb2I9MCE0xsX4jU7OrQmMGBzqXqZzMFE2m3nGt/guzs4pnlztNcvvjD7ut6qLFCS046djJk51rMJcUingVzJ2Qq5/0Q7bkxD8BAC57ugMA8OM13wAALJyzECcd+4w+Z547Vm7yNovcIf+42AI5Tf0J2wsaFIXqTkX2FOSpofhigVBhksN2eIiY6owU3qwdPEIOHOY8Lcm4C2BM5USwJCD1YyQo3QbmZhijJM79oOslx2miE9QWd+IhVKvwTMkRNd2Pn3Ul6Zm7fQ0AsPCxQ4DhA9WxRG3GMdVSKPQFjCpxdVRSTyVsb2gsR+Wqp1IpHnriRU4q+Kquqi9VmEKfJJwXimLzUh5wx9tkJlNuFhXDQ74qibcVk4LWuEkJbxZ6pDYzbbqhKGY83KWOXx+APxx7HwDgtpeUc/Zlq5Uk/cR0JUnPnKPzWQ0P5Sqz8a4BwjPDcrUiN4TA/81imVucZguPSEgYA2iMU4q0uaDLWSTnppJLO7jn8u1SKR6Kwk2UNKyAOssLwIuFrXDJpdVGaGvznV54aApvk7sw2llO+PXHQlECs4l3Dfrzh0+rz0V9Knzk0zOU0/N3nzzEaRtZlquWOCclMPNskHcWZWSpgbokqhDiRCGEFELsw76fJIRYLYS4tJ52EhK2FvVK1AUA7tWfX7e+/waAxfV3x8KlOUSWS1DJ+CsPn46FGttOGEYaR0x3BO4GaGkmDPdjelRPOtInZdcbv4NnhvXcC7lJl0vnatV3QOb6W8nGWXbHiXI5l6Q//BUA4LtXzwIALOq7CABwQqcKaT5+1wGnDVvP6mljCB7/ZpzUDnuRkfu/LfSoQoidABwJ4GMAPmh9Pw/AVAB3FPaSkPAaUY9EfS+A30spnxZCvKIf0L8B+A6ADwM4pv7uWLg0wdbB8Wx1XCfKTg1myjOclA4Kc0UDFhIR5KjsWPPZq/PPU84putKBzXHuGdNKhHSRNUO6c8R4qCiVzP9LfqdSGJ27THFSCsT7/N8OU8dyR3N7DLEx0uxFIUPDzFXQTh/EeLaXZr0G6nlQFwC4RP//C719OIDbpJSrRUFxASHEWQDOAoAZM6bU0V1Cgo+aD6oQYjKAtwM4QAghoXLDSQD3AJgvhDgHwE4AWoUQm6SUX+ZtSCmvAHAFAMybt7crEEMShLugRTipx1GJy7a25pI05rjCxxjaz1wEucOK0W92uUU8jN6wknNDT39Kdntyhmll4cu2ZDP28czfZ485El6z7pQbcNif1gIAVvb9HgDw8Fs/BQCY+67n9bGHmDEr6LZsJ2juDOTptclixXTDlL+2Us41ApGwoloo4qj/DOBaKeVMKeUsKeUeAFYAuFxKOUNKOQvAFwD8LPSQJiRsKxRN/QsAXMS+u1F/P4LVPkPMmlEp+wFv7FQv9Jcn+qqW/NcvYmP23NwsfsfTThpQGLL2OIrmtEfJ57cEkkraVc5LcGaNt15XRb6frnVFdydW9l0DAJjZeRwAYO7v3+WM3ZxbpCtVg3WvIXZtQUtimAPXo0eteYSU8m2B737Atq8GcHVhTwkJrwGN9Z4i1FuoDMg5EGkFyJ5ue1nxZBXUfcG2zVG9IEMtGeREHYZsAvSY7R05V6yV4AKAyRrtSR9bz8p9DjhHpfHp7cXHKJ/S05c9AUDx0j8f+VkAwOwpfW5bvG0OO/lGRILmqT3ZNfKQGFuK8lmQJ9cIoDFOKbF4bj7Nhr7jzre8Op4V7cinbZ47NFpQwoLnfW8eoi3OOZ6JtVz2MxHGrpurdCzaETPl2hTD3r90g3pxaeE0s/M4HPzH48yYnLYYTaiJWCQrxVmZFwfhax4eCqu9gOQ4nbD9oDEStYYKyisQwVQ5sSx7BFGt5s4okaIWHF4Mz9AQAD3VM2phJCxJdhoHN11Wq77ZEwXjYtO7k/maX682oV49/yEAwP+uVE7ad/V/HwBw6f7KVW+/SZujWQyj6jpmypUlyxDDnWEof6uVm0G14ZqvnUzgNA7KtUVF3GogSdSEpkBj3PwIAVc9r5IJy+EpSyzKkbdlm+r4wis2LpIgJLXLlTzTs+GkWsJyTkrhLJnrhC1LJUu9E3HVi4XKwA5YdDmoYNzvplWqX5KkR7erhdPHHzrcjMfjpEXGE8aD1SwF9xporMP6HB5ew4MMbTMs8Vor42ARkkRNaAqMMkfVmVJqFMQ1/IUKgOljqrvtpraZKsdTMdkFwARzK+OqHfokZ2dbbcS5MSnnY5dGyvJx1urYzths98eNGkaiwrk2NeaK0+/AF1ShsWvuVTH3t67/FgDg3R1fAgC8bwYLq66UvcouxmhgFV5zrpHuKY/Rt1FUUYXaorh/INccaAkqdci3UfnVQJKoCU2Bhqz6o2XIAWC9CigTK1erYzt1db1p09T3ZVbRhGddsXigWYEyXmmOXa8V4P06q57O9CyqVYBy5htzYovbFkl+ciyxa0PpbTG02b3uWI0qLvHNBWT5L6SPvfIeJUk/t1QF4tHq/uy/Huq2YetZQ/ppat8+NnPnC5Mxprc3z8iyRfNnroPlGgTWlmxtNZLZaFJG4OaXJGpCU6Cx4dIaRtIM5Do/OZ0q1e1Y++SYtLBgOBflGN24Ue0gXtmhEi84IcncrBvTMvAEESQtbA7HQ0+KUg1Z6DvrJgDAqXeo+3Gb5qTf209J0lOP1Kl2qge747R1szEnHH4s3Q9eiaWtzbc0BUJe7P0m6YZ2LBfDw3ntWR6mXgeSRE1oCoyyRA2HopCEExs3GiuFpFV+zO2vVla/mFsfWZd0Vr0qVR8hrmhxWi+ko8g1jVuuqtVcgnNOyttgksXef/8KdR9uW6+8LY/Xq/tPPag5qZaknn7X1oFGbPteHQCmUTGSdXgo/y6WYZtAfWlpbKxPbW2mhDshJUlL2O7QGD2qrlCXrXlZfUtSYPasPN9mIN2Ps61bDFU+obc+e+ZZdcxL3c4oqm/cX/3D7NOOXpFVl/ZS6nDpxHwRUClDkNDhoddsluC2/tPnPAAAWDK4Go+tV07Pi3SV5sNmr3HGyvWsIR8EL+Q7klqIJ1wTVmiMoNDpFpZpuqS1IcRJtSbF+AJobYGcMMHXbtTjtaWRJGpCU6AhHNVUlevXq8pdVKixtNLgFKYojPiQCiulj+hVOlm5Rq8825RUlO1q9emlf3RWtozf0gqZxr7DeLoiBw7/i+mLuR8s23/34IMAgFXr78AbO04FABxHVZo5N2R6Vs83IMs8TuqNg1uuhlwHbpRK4GmYSJJ6DuR0fyZNys8F856KhfnUwKg+qGJoCNkLqyDWKYKNdpXrEzvpz6qfcdrbjmWCswwAvKSPaNcqownj3XOoDfrHcsfjCyGj6KZMKDp7n/lReRaQShm2t78N83LpBVjLotsBAA9dqo5btV7l9JjRcSzmjp/uthsyZ9rgXvvWPY2mna9EDA52xkK7AAbge+XTA0r3hxlAIKuF7oW1kKb+hKbA6E79W7ZALH0W6FROCJWDdBZj600vnOqprch+0bPOUIvqXLVootgkzymbteVkSmHOMSa+idoa4uoaty1UM3+Kj0z9z1+t1HOHLL4SAPDXoz4DAHjT7f/kzyjuyHMDgwa5CEZzRQFeGfSYe5+BPfWz4sXixRdVf1RTQDuYmFnNnuZjhoeknkrYXjC6EnXHHVGdd0BeFnEbQJCLni2VSSXC3NSijtMBkLOL1M4wUSkUc9lr2ylXB5HLYlcXAODLe/0ZAHBtn+Km3RuUOmr58WcAAGadptqoWO1GwQMDCVZWFu6c40XpcjVZACZLNw+u3EOVQ/cqy4Scwo0EjRgLaiBJ1ISmwKhKVFkqebmaAOSvi21eNS5otbeNgwlJh4kT4w6/zJTK9zscssIMCUxyeVmTeTjJwOZ8ddunc+Frlc3GYSXbSZJOmaTMoXv8RmX1rFiOMF5oTsBEam97PLhSNqEfRdkMg6o+apNmEpLOJgTHNz8HYUt8L1NO+JQRHpKQ0Hg0NlzaWwWy//k+ANnKlQAA0aOU+HKq5qNaiS/b2jydZkz6GPD91rnRlSo3AzLOKltbjeR44JwVAIBbXlTOMFfoOk5n6TpOE8cJ977YutBIP05Ytv1JsJyi68rIErgmYzbe0J+XkORclRkgYlm0gTKieavqUPwniZrQFGhsSp+tgOhWUgmrVM5POed16pMCxCzJUpStmnNWR9JyLsja5G5/3iq3WjXBej96Rq2Ir117HgDgS7MXAgAufPQgdU4rTxdkSfPYOGJhHKH0QUbqMg1BgSuj2KRCdERfXx6Kwuu3FplDrXCbGDeOVme0Dy08IiFhDKAxEjXi5hZ6K7OXlSugeEl9YqMu/z13L/WpuanN6wqz9/HEatxOn2U+F2QZqMM1R/P9az74a1z5pLLTX7v2AgDAKVNVKfH37q4dxVeuUud2qlAYUzHabjuSzc/TeUbSB4lASIwXTkJSkVXjM9K7q8t3c+RlziNc2nGO5imcSPoX+S8gSdSEJsGY4Ki13L3Es2rFXPmbkj4tr1dcqfr6OeqAUKhKzC+giBtbx8UyPUfTP7I+fvj4dFy0QnHShXMUJ/3asqOc7qQO1zCeWLEEGfC9x/j3tdIHFfFsfl+ImxprE3yHcT/DNms7lFI0puutI8YvSdSEpkBDQlG8YDbLQyhbo0ItxJKnAADVlWqV37KvDp+ePlWdFPC7pE9vdcmzIMcqyNVTuc8LKlTbfzzmHgDAdc8rn9efrT0f5+rV/Rn7qmQanl2eOy6Tz6vhbFZIDLm90rXxNJSx9EGhUBQeFs38d+WOyj/Y8Fu7urThpnq+otCUCguJYSmQ5PgdIDZvgg3T/zaqM7UNEY5CtV3SRI9yqq6u0HFV05SCuXLs0eqYWFRq3oOv9uDT1DDLNEfn2pmeY+513JNe77/sGfWD3dx7PgDgI1O/ivOX6kJjEWdvY4akByKwEPH65aV/zI/NnEFsqsDG7GUv1G6RpOKzY6UAqBeJ/KTppRKW9z8ALwKAHmD7pdTxXSNx7yPUNfXzor1CiJlCiEeEEI8KIZ4QQpxdd48JCVuBeiUqL9r7EoDDpZSDulbq40KI30gp19RuJjz1Z0ufBACIlWtQfVE5cGR7q6kek5n6KRKaYjAwgKy/Xx2jJYSk8o98ejQDYOoae+ov88WJvmU/URlMLr1CqaBu7v0mAODEyarc1smzBn2VEl+IxKJC+fewpC/PyMJd9bRJ0w2vcdVP3vWTiZXugynqpiVgSH1kh6kAubQUbsiMrZ4ymVKYI8s2cZwOFe2VUg5JKXU2Voyvp52EhNeCeiSqV7RXSvmwEGIPALcCmAPgi8XSFIhxVCxR8fcDD/Zh3BQdJ36WirosjEqlNvR21tNjeC699YZ7xdzcaHR1cFTafmqRCmKjrHqUC+qTD+tMz7bUiphbCYYrU19WWIvkHJDGyjgqMpZbn6SW7bBMoNh7XmiYZeozKJV8F8BqhKPGTKqW66QYHHLHiGLUIwkXQBXrBfKivZBSviClnAv1oJ4qhJgaOlkIcZYQ4iEhxEM9Pevr6C4hwUdNiRor2iuE+KKUUgKAlHKNEOJxAPMB/Iq34Rbt3UtCVpHdcCsAoP93vc6xHWfuBzlR8aUqN7fRmPi24TnqGzlpkpfJ2CsAFrtgi0tyyX37O+4FAJy+7G4AudPz4+84BwCwzwmrnTZqhUt7GUq4a1xrvrIvckLxMrQQD6dQmgk75nmg6N6R87rHwxnvtMykph+jnmP8le5/xCwuBrd4db48V8UaKJKosaK984UQEwBACNEJxWGfKuwtIWErUcRRY0V7LwcwrKWsAHCxlPLvxd1pjtqv9Hbr+xTP221/5WhSOfyQuJ40sm1Wm3Ree7uXJS6aNY7rUfWnzDIvwcKqzWqbh4/s/dv36HH4mZ69pBncvBjSecLiilsGAZYVmpsfo4V3bSU/5YnSySF4QoqortjW53JdNNeTshKTEm4QYOh6+TXUQs0HtUbR3h8EDk9I+IdhVC1Tw8/1Ys37r8faDcpEN2+B4j1i1mwAKqhtpIV2jW7OTpAQqZjiOXTEBpplePZdim6f8Vd1i+7rvxQAcNkbvgoAmLEj08kGqvMZ02VMj2qFNAMAdHZmUzWks9NIVMecCeSSlpkujTueFW5i7pGGV32QdvAZx95PUpmB63c9q5YtnfkskPKjJmxvGFWJunm4hEfXduF1k5TjcPVj/6J2WMF+XkVo+ifGUUNhzUX8lrdBfVlJJu7v3hkAcF+/cno+ov2TAIAzHj3a6VeyFazN96Iugrx/mgF0Mg2Tgsgqg27a5aEmnPfRzEIZosfv4CeciOiGPS5tj5Py1nI9KU/1Uw2Hlds+B0UZGUNIEjWhKTCqErVj5zLefXI3RJfiqFVyt7Pf5IIwlWjYykhCUTxnY4UfHKB0pPetLePXvf8NALjqjSp85LApr4T745VOrO2o/pZzVFopUyoiY6EayqUhn0F0dWnZpl3yiHdS8jbtqgf4vJFzZc8ngfdl1arinNT0Qf0T+LjL5bpCj2JIEjWhKTC6KX12mQx5zgLIWCqd0Dn1hq3UE4pSML4fv6xsFsv6fon3T/4PAMCH7z/CacPTKMQSDpfL0RAQz3G6wnShQzbPY/yWc1SSlpRIl7ikxUuLwmmi4TWW/6rRJlDqnnEszz+3/fNxW5nAvSoseI161KaCPfXzKYYfq28QqaDWvaocO5b1/RIAsE/nv+KIqa4nvTe1FU2jgak/lJHFOYdgT7lsCuVe+SbNPBV9aGWZoKtVj3J428y9kDtYS8v44S18WNQE/967B/b/IrAvgjT1JzQFGiNReWaMzPo+Y8fUu60hs8zsixUto+nymAfVVP/i+j8BABa/5fMAgMMXH5+fw6ZlRKZzb1oNuMblGVlqUx8n8zVX83SrDDFCFzc2C7A2linGPo+HfvBriBVIs8JeuKO2KA8554KbRXlb1u9ifrM64vkJSaImNAUaI1Fjb1KtN4wXmKVXjDtLVKtAC3v/9JucvbAKAND/X3cCyCXp7h3KpWGXNrds+YgQiXN34F0Dk3SsLQH4AYAkKdt1Vj3taOJlIdFN1VLXeRkImSR1MhMGwmOcbba49NqqlOPXXweSRE1oCoyNVX/oDXsNb58HfW71RhV7P/nanwMAek/5EABg4hUfVsfxWHmr/6gGIQJZKvlcnOcFNfvDM4nMMlOO03BSXQcAPHth5P4Eq5HEQJxUZ0pxXB15bimu2mI5sYz6aoMOtGQFe0eKJFETmgJjQ6La2BYSVOPiuar6yLfX/AYA0LvxUQDA8LdU2cbsA/NVlzwTM+ow80WcURx+V5j5Osz3snvvV19sHADalSSqvO0o3ZFuk5eB5LAkHeek3phZ9hcvU0oIXI/MslQbzkouhtZs4t3bOn7rJFETmgLbnR7VfjuXrldOxSRJJ0/UlQI/fwqAPIAwVD2l0L0wZg60Z4RKWHLE3A3N9wMqZYJ8dQuElqiedOYJOCJt2xw1Ws3PmHK5hSpvI6QnDl2bB7suAb9nhJRxOmF7wdjQo3LdqH1MVG8abvqut/4BN69WukWe6Xm/Dso8F3aCcKqiFIXAZKSbfUFtd3ToTzd9UOhcj6NRuMmrumZWhw7zntwOObnDvW6uhWD1nnhVPFuPmn+pj41ljSbYnDYWTsM1Ctp+T/zW5uVRvh/5PWwkiZrQFBhdiSqg3kC+UOW8s9a+WscCuHr5eFzXrdLs/Nt0lWbn+8vn1zU82zZfxFGFTsRmqrTo8GFTU6DGuZ4kJU8o3Wb1YKtaSiS02oCkkWRcWVpSm/XnZXqWkVW41XfUVZHn4Y+0JZ3+wcaDQiSJmtAUGF2JKuHoFw1s3hnjPhGOunBfpSt9VQuW67ovwMlTFCc9cbr21eR8l4NZV2SW+ZYo7pdKx3bp9EFaktqhKF7QHOtWbOx3xiW13d4JeaZ+yM+U8zwWTmMcm22JZwIg9UF0G3g4SY1qhZIl5DBg1aZDySs8xH7TGmh69dS3V6l4+3JZ5bH69IyF+M6zR6pjYkV7I7Cnq2i0K5Vd1wUiqjuraFVPqW+d69EGXWiYPk0e1061ELMXIoKrlHjBW+6sElAXmZeLYrMotxV3JIlkMJFZ5hkazCKOXZsXq2/fc/pf8t8/qacSthM0Vj1Vh6KXo+XG2wAAw0+qRQxJ0lJJTcFtta6oYOoPHsciRXOFv6jdpr2PL8CodHts8WpjejEAAAfsSURBVGKfF2u/yFkmVOq9yK0w1lZIKtI1DbMx87bscyOhOSkUJWG7wdhQTznH1DavnfnvymH4Jy+pPG2DX1OuetnXTlMHhJyvi9Qg9ahJyGEj4l4XVGeRUeAZlVFbLHlaHTNrmmryzW8OjtnmiuZ/FkZjjuVSmeWqkqVSHkbCUJiVxnZ7NFn64Iwjxkkln4HshRtXl3FH9wCSRE1oCjRWPRXiqGzf0uNuAQB8f6laEV/18vkAgI/tqrLqjdtXOflWbBVUzNwa5aj+/miGFr6ftpmzSLbsKYjndVmDHXSs/QyVPV527exeayS8Q9rtsn0EHqoiqixTXii7oeduGJnFbH7O2vAkKG/DXBvy4/l9JkmaivYmbC8Y8xz1t6uV9Lnq5fMAAKdPU5L0RyuVIzHLOq9evZGa6gL7R5pV0NNBrl6LyjPdAICWI/dWY513UGG/vO1oFRiCqTTi6kqd/RGp5+mKuaui1YaXt4uv3Ll51uh3rT749XJTbg0kiZrQFBhzHPXF910PAPjovWp1f2f/9wAAF+y1EADwnuksq17IHFsvRzWWEp8rjpSjtjy6RP2zQvHSavdGtLxphjpm12mR/sPjsvmu5yits1IbZ5gOZbqtTtN98ADFeq6Jr8wF27ba8bJF8zAeGb9ffoZpf4wxJIma0BQYGxzV4jd/71bhwHf2q2Isb2//HADgS48fGWkz4DdQL0eN+Rxg5BwVvdpF72Wlhchm74LKUUeEx8z75xw1pEclKfiqrnvQp8OQJ010xxErtR66Bs4vufOKfTzn4qHQG/sauJtfiKNuaze/QHXpA4UQf9GVpZcIIU6qp52EhK1FvRKVV5feDOAjUspnhBC7AXhYCHG7lLJ2DUniqAw/nbcYAHDR6sewvG8RAGDRwecCAA6Y0sPaYPw2xFG5pSOWzCLmbmgjwlGpInZl0SPq1J2VPjM7Yj91XNfkYu41Ao4K8g/Q+U+rr99TbXe4oSoe77Tb5qt743QdtvU7fLTK7jsf85BrqQqOY6T6bQuFEjVSXfppKeUz+v81ALoB7FLYW0LCVqIeiRqsLk07hRCHAGgF8FzoZCHEWQDOAoAZM6a6HkH6bVvarz6X9y3Cnp0nAACO+8s71TH8bSvwBQi+etvA1u+lndyo/FKHu5Umd8JcJdkqB87Nx13kZ1kHR6X6TmJQh1B36iDCmD9oqB/ORbl0jIWP2OPQ/Zs0mKRdkGyW4uOwr/E1cNR6HtQFAC7R/1N16YcBQAixK4BrAZwqeb5zDbdo794S1SqWv+dGAMDCR5Qy//qeCwEAn5mxEPu1s2l6BDk0zXn/APVUy133AgCG7lwBANi4Rn3ffowy7cqZuzltiHI5Pg2a/sPjytZo0+umTaaEppypVV2UcTpWJpL3lWXxhQ9zcPFUT9YDbnJHxTJwR83AVop27u1fFHlhoeaDWqu6NICJAG4F8BUp5f2FPSUkvAYUSVSqLv1x+kIIcTdUyfOvA/iZlNIrfR6FVk/d9ZKis9f3KAeTk7r+EwBw8fKjiqf6ImzNFFPLhEpZop9XWfXWPKEkS0eniseqfui96rhQtGjMYaNGvwAgevvUP/2bgHZd6oi5F8YKV3hTro2Qyx0fs32cbY4lW3UkAzenE1x6OlGo3rjCX4/kkAUAbmLf3QjgGgBHAThNCPGo/juwuLuEhK1DTYm6ratLL3u0F4d2/QKPrL8SAHD5AcrB5OhdlbQKcpWYOiq2vY2Q3XArAODJq5Qo2TysuOG8k3RwX7sKfamLQcdCL1pIia/azJYu09/rHAGzdjeRqTFEM+OF1FNFGUm4Q0uIs3sDYLNCPXltt+I327a/bkLCPwijakLdXOnBI+uvxEEdZwAAPvq3o9WOUDAbdwEr2iaIrPa+WrD2y55NAIDnNqiqIwdOVYaH6pnaCMelQogPx0It+PYWvYLu0+HT+81Rfe0xI94PjTPmZkew72mtMPXQuZZE9UJi7Pbtxrimxb4HW/u7IEnUhCaBkLLerPTboDMh1gFYWXBYF4CegmPGGppxzMDYG/dMKWXQwjmqD2o9EEI8JKV8c6PHMRI045iB5hp3mvoTmgLpQU1oCozFB/WKRg9gK9CMYwaaaNxjjqMmJIQwFiVqQoKHhj6ozRriEhj3TCHEI9rn4QkhxNmNHiMHH7P1/SQhxGohxKWNGls9aOjUL4S4HsBuAO6UUn5dCLEXAGmHuADYtzDEZZQRGHcr1L0c1BERjwN4i45+GBPgY7a+vwQqOqNXSvnJRo2vCA2TqM0a4hIZ95CUclAfMh5jjFKFxqy/nwdgKoA7GjS0utHIG2pCXAC8om+aQVGISwMRHLcQYg8hxBIALwC4aCxJUwTGLITIAHwHwBcaO7T60MgHdQFUaAuQh7gAcEJcTo+FuDQQwXFLKV+QUs4FMAfAqUKIqQ0aXwihMZ8D4DYp5eqGjWoEaAhH1SEuqwGsg3JxbNGfM6FCXO4CcOGIogdGAbXGLa0bKYT4KdRD0PDx1xjzPVCRGlUAO0HNXv8jpfxyg4ZaG1LKUf+Dikq9nH13N1TUwB8BfLYR43qN456gtzsBPA3ggEaPt2jM1vZpAC5t9Fhr/TVq6m/WEJfYuC8H8IAQ4jGoh+BiKeXfR3twEcTGvCBw7JhFskwlNAXGlBolISGG9KAmNAXSg5rQFEgPakJTID2oCU2B9KAmNAXSg5rQFEgPakJT4P8BA774WC+MPXQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "curr_fig= plt.figure(figsize=(5.5, 2.75))\n",
    "ax= curr_fig.add_subplot(121)\n",
    "ax.axis('off')\n",
    "\n",
    "import colorcet as cc\n",
    "subsampled_color_ids = np.floor(np.linspace(0, 256, q, endpoint=False)).astype('int')\n",
    "color_list= [cc.cyclic_mygbm_30_95_c78[i] for i in subsampled_color_ids]\n",
    "\n",
    "# Plot embedding with color\n",
    "for i in range(coords.shape[0]):\n",
    "    plt.scatter(coords[i, 0], coords[i, 1], color= color_list[i%q], s=30.0)\n",
    "\n",
    "plt.plot(coords[:, 0], coords[:, 1], color='black', linewidth=0.2)\n",
    "\n",
    "# Plot Pearson correlation matrix\n",
    "rho_frequencies = freqs[bin_low:bin_high]\n",
    "\n",
    "freq_ticklabels = ['A2', 'A3', 'A4']\n",
    "freq_ticks = librosa.core.note_to_hz(freq_ticklabels)\n",
    "\n",
    "tick_bins = []\n",
    "tick_labels= []\n",
    "for i,freq_tick in enumerate(freq_ticks):\n",
    "    tick_bin = np.argmin(np.abs(rho_frequencies-freq_tick))\n",
    "    tick_bins.append(tick_bin)\n",
    "    tick_labels.append(freq_ticklabels[i])\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "plt.imshow(np.abs(rho_std), cmap='magma_r')\n",
    "plt.xticks(tick_bins)\n",
    "plt.gca().set_xticklabels(freq_ticklabels)\n",
    "# plt.xlabel('Log-frequency (octaves)')\n",
    "plt.yticks(tick_bins)\n",
    "plt.gca().set_yticklabels(freq_ticklabels)\n",
    "# plt.ylabel('Log-frequency (octaves)')\n",
    "plt.gca().invert_yaxis()\n",
    "\n",
    "plt.clim(0, 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fig 5b: Cube-root loudness mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_str= []\n",
    "q= 0\n",
    "CQT_OCTAVES= 7\n",
    "\n",
    "features_keys= list(features_dict.keys())\n",
    "\n",
    "setting= settings[1]\n",
    "    \n",
    "q= setting['Q']\n",
    "# Batch process and store in a folder\n",
    "batch_str= [setting['instr'], setting['dyn']]\n",
    "#     print(batch_str)\n",
    "\n",
    "batch_features= []\n",
    "for feature_key in features_keys:\n",
    "    # Get features that match setting\n",
    "\n",
    "    if all(x in feature_key for x in batch_str):\n",
    "        batch_features.append(features_dict[feature_key])\n",
    "\n",
    "batch_features= np.stack(batch_features, axis=1)\n",
    "# print(batch_features.shape)\n",
    "\n",
    "# Isomap parameters\n",
    "hop_size= 512\n",
    "compression= setting['comp']\n",
    "if compression== 'none':\n",
    "    features= batch_features\n",
    "elif compression== 'cuberoot':\n",
    "    features= np.power(batch_features, 1/3.0)\n",
    "n_neighbors= setting['k']\n",
    "n_dimensions= 3\n",
    "n_octaves= 3    \n",
    "\n",
    "# Prune feature matrix\n",
    "bin_low = np.where((np.std(features, axis=1) / np.std(features)) > 0.1)[0][0] + q\n",
    "bin_high = bin_low + n_octaves*q \n",
    "X = features[bin_low:bin_high, :]\n",
    "\n",
    "# Z-score Standardization- improves contrast in correlation matrix\n",
    "mus = np.mean(X, axis=1)\n",
    "sigmas = np.std(X, axis=1)\n",
    "X_std = (X - mus[:, np.newaxis]) / (1e-6 + sigmas[:, np.newaxis]) # 1e-6 to avoid runtime division by zero\n",
    "\n",
    "# Pearson correlation matrix\n",
    "rho_std = np.dot(X_std, X_std.T) / X_std.shape[1]\n",
    "\n",
    "# Isomap embedding\n",
    "isomap= Isomap(n_components= n_dimensions, n_neighbors= n_neighbors)\n",
    "coords = isomap.fit_transform(rho_std)\n",
    "\n",
    "# Get note value\n",
    "freqs= librosa.cqt_frequencies(q*CQT_OCTAVES, fmin=librosa.note_to_hz('C1'), bins_per_octave=q) #librosa CQT default fmin is C1\n",
    "chroma_list= librosa.core.hz_to_note(freqs[bin_low:bin_high])\n",
    "\n",
    "notes=[]\n",
    "reps= q//12\n",
    "for chroma in chroma_list:\n",
    "    for i in range(reps):\n",
    "        notes.append(chroma)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAACjCAYAAACZvbxxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deXQcV5X/v7e6Xu+txdqXtuVFXiXLsuX2ghPHgST8IJCYDAGOISFMGAMHQ8gi+A3DMMyQARSHGeLMHAwhMAk/ZiaZkI2ZSVhiOwmO1fEmL/JuS5Zsy5ZtbS31Ut11f390dbsltaRWL2rJ7s85OTkudVW9rv7We/fde999xMzIkCGVSOluQIbrn4zIMqScjMgypJyMyDKknIzIMqScjMgypJyMyDKknIzIMqScjMgypJyMyDKknIzIMqQcOd0NmIzc/sNv2gHUA3AAcAJo+P23ftSW3lZNXSgTIB+MJrAmAFYAAoACwAWgJiO0+MgMl8OpxzWBQfu/VTueIQ4yw+VwHK5LV4Tf7QFJwXdQknXCUpC3Is3tmrJkRBYBEWVXferOtrw5FbXWwrxQTwbV71eunmo9RUR12qEAgGZm9kaev/aRZ4bZcjuefPCGH2IzNpkGES0F4L/tB/VdGMMmIyIBYAEAg3a6Wn3Pw725MxY0RjvvRhfaDd+TEVExgDkA9jOzCwBu/+E3azDK7JKZFQAHIq6hA/i53gtnbESQZYMFptzCSFtu08R9o8nHDduTEZEEYAWAK8x8PNHrrX3kmUYERQmvqxuenk5klcwGSZJzx5MP3tD23A3Rk930zF2DbKXT/3Hw3wEQgEZm9ifpNk4AtQCEwZoDvdmGnvMn/XpL9sEkXX/Kct33ZJrAmgBYVSUg+lq6/YZpJpcxz7z4nQdfjWorxeOM1Yz+YbbcsTeevedi83tXmPnAaOdfz9wIfrJ6AFZ3h0v0t/cia8402ZhntmAEv1eEM3YjgiLbCKBJOz4imnFfA2ArgsLcCqCm4/DObQDaiegWIjKMdo3rlet+uGRmR9+pq8KQb4Gt2Bo6LKDZT1EYzRk7qgGvCW3YZ5j5KhHtALCSiM4x89nxf5Opy3UtMiLKnfuF2vb85WW1kiyJiD8pCPY20XAAEP2Xr0Lpd0OSdbAWFwhJpxtJlDHBQbvkPSKaRUSO6s+8dkGSDYOG5P3P3X5dujqu2+GSiCoBzChcZX9IkiUXgsICrvmvGkY41em6eFkhIuTMKIO1uACujk7/pcPHW4hoOREZE2kXM5+eeev3u/s7mw8HfP2DhuQl9/1+1CF5qjLle7Jlv10wyEhXfWrD/s8ctwNoZeYTAHDTM3cN83uNZPQf+PfXXqv88NrPmXKzzQCEpNMpWWXFrqyy4kcBnAdQrdlWjKDX3xU6N1aPf1bZ8k0AjANXTgidMMOQVXZd+9Sm9OxSE1h4RhfoDygDpz0enVlXe+SxM6fGez0iWgbgwm0/qNdhjNklERGARQAsADDvjgdcRYtWv4MYPP5L7vt92Kfm678EZeAKLAULAMC5/7nbrzuf2lTvycJGuue8D6pHFbZqCwA8hHH2CES0AsBpZu7UDo16vmZjHdLOJWN2wXO9F07bwCzbiitAkm603insU9NbCqHTZ6Hvwj6/OX/+7vG0eaowJUR295+yhw1Dr3ywp83bqXzA3+0XAKAvEDCW6oHRZ45RIaIPIDj0dcXTPmbmtY88MxeAzKyi++xRZJfPhaSTR2pLA4AN0F4QnTAq1uIlrtN/+us/Et1dzMwd8bRjsjLpRaYJLDwk+nrUWle7el/xt8WG3L8oP2CpNFXhmrsBGH3mOMj7zyo7D/1459sA3mDmvgSb6gRQSySJnOkL0N12FNmlcxRJFsPasv+529uW3Pf7QXYiETX0nd/TRkSLiCiXmY8k2J5Jw6QXGYB6/wBbe1sDwR4rSxLTFslY+bh8R1uP8TsAPo7hdlDUmWOk959VFj1HLy9d+LWVG3R6XQ2AREUW7p2ISOTY5yvdZ494AHoq2oc1d0U0n9phIiohotUA3uOpbDRrTFoXBhHJRLSsc59/3cBFVeTO12HaAhnWMgnQhsQ9nzgS1cuuHY9GPQArB1TRc/QysubmyTq9LilZr0M9/kS0NatszqKDL/24lIisY5w+CGa+AGAfgHVEZE60bekmbT3ZtxqzhtlZP1rZ1w5gPoAsBBMDmwpq5W0A5mKEIVETVKxGvgOA6Dl5Fdnz80ESAXHYcCMR1eP/k6+0EdFNRHSAmXtivRYzu4loG4JRgrNzbzksYcjzOrZt4ZRw3qalJ9MEFo4P9l5UN7YfCByatVp3J4DzzNzIzLu1vK0GBIfAWJ2po+Hsb+/1W8qyQgILXW9EGy4ZMPM7ABYRUe44z2Nmfi+r+O5Kr+vYIQyJp85b1zwlnLfpGi7rAVg7TwXE+cMBgCDKF+tMn/yx+fahb/srH+yJOiRqx8fFuT+cegbAgLDqkyHYccHMOwHMJaL88Z5bPP/xe3T6PNNA9x6hmWhTanFLuoZLx4XmgMgqJlimhXU+4rClCSohT7jmPM1Z8/OPVyFG73+yYeZGLTQlMfOlcZy6RtbnC0lngbfvEIxZ1UDwea1JTUuTS1pEdvRPytkZdXKtKZtidj0kgVUA3nvnwVd9SGPohpnfJ6KlJUvutxct3vB5xBYglwBA0pkgyVnwey9BNhSGj092JrSRFORmZvytKZuSZWfFct9ZANqZ2ZeK64+Xms+92Wkrc7zl6788LEBORDYiWkBEy4iojojq3N17DJ7eQ3D3HoAa6IfP3QbmAACoaf0iMZLSnuynh3PDM0hVZWdlne7dE7sDr7/6bffAtxqzhgWtf7iiN+nDFhGZABQy865kXzsB6s15lSZPd4sYuHIVAAmAbcrA1ScAfBdAW2Tgfd665j8AmAVABJRueF3H4elp8ptylr6bpvaPi5SJTBNYMO05wOJYo3/p135m3SDraSeAAU1QcQ1bI4WZQn+PzMyYsankXO5q29cT/T5JxgFAGHMqIo/JyMNMZj4W5fNhR69O5Ahj9hKl7+Lrvp6Ol14CXp6A5iZGKofLegDWgMLimNOPuQ5ZlvWU8IwoIsw0aKjRjkdmZmwcaPE4suusd0p6qUk7Pllw4pqpEGJEm1Tzh4Vn2JKk32or/OiC3o5XlhNRUWqbmjipHC4dAETLwQDmrZAhJc/xGfTaqyz6WlRAgpAErAEP/p4+RN+sfmbO3+gskjUwoArSEWTrqNkQ6WJQgBwx2KSa0Aa1nwibATxERC8yc3vqmpsYqezJnN2XVCUrn0ICA5Iwg1QD7Og6GhBdRwOwlEmwlEowZEtCn0WLARj9vf4V3g5F+F0BmOzhdRtJ8+onA20WOcz3N970ay2u+WsAK4hoXtIbmiRS2ZM1XGpV75+7XAZifFtHg4hkAMuWf9fUWrxKrpXka+4PnR6KsNJOZj677LcL3gGwEOPIzEgHIwXIxwszdxLRKQCCiJYy897EW5dcUtaTfbmqO6+/W3VgyNv6pUVd43pbiUhPRKsALAWwu/Qm8Ygkj+r+SGYYakrAzPutxpvmzZ++90szS/6jefGsi1tqZndOGhs0JenXWubAQmaOOdNzaMC8bZ//J7/5srsQgA/AXmYO+4TGM7sM/X2UzIwpT83sTntA7WvyKqdtRv18ud/T6LcYV/ZLpH8ZwV7dCaCh6VRBWp5BqkS2BsCfY82FigiYW31uFpeOq36dHgMFs6SqzTf3XbfiSBY1szu3ANjo9h4UBjEHREb0e3bBpF+o6nTZEiLWG6RDaEkfLoloBoCz40y2q/cNsLX9QEB0t6soWyzJJQt0JtlAUyIAPAlwABA6XR4Cah+ICFbTKniUI5KqeoA0B9STKjItCG0fzwppIjKdafTf2n1OFWXVEgordQheZnLNCCc5TgCKRHowrkXOzIZl8ConQv9M2/NM9uxyGYA90f4QGWIC4Dx3PPBP3/9EXzEAT8Vy3VskUSUm+YxwEtMAYAOR3spqr0AwpikRCTACoc+k7XkmrScjoiwACjO7h/4tIsS0MeBnx+n9/i953bzvyZ1Z7cy8lyS64WaEyUSzs2oI8s8Dge6DAJ4D0INJ8jyTZvgT0Rpmjhqw/enh3C0ANl5qDYj+Hsb0hTroZFIAbP3Soq5NQPR07FQEzK9nNHOllpn3ai6Mepd75zqrafU2TPXZJRHNBjCgLYAYxk8P5zaqKjtaDgYwq2bQCO380qKu627FdDohomXMvCfi33MQrCYZ15rSZJDwcKmVxSweSWAazrbmgN++QBd5LGNzTQynEEwTShvJsMlCw9uI+DzcoPjgFnqaFDbCjYTmSqIxP5hC4hLZm60F9jdbC7b8++68PT/6TdZX32jJLx7t81+v68mXdFiKBENMGaYm47bJ3mwtCHvnD+9WxKI6EfYm3zGjc5hoNGN0JTO/l4wGZxidoTaZdmw5gN3pWo0eT09WD8B67kxAlM/SAWN7kxcjKMoM6aMdQFm6bh6PyBwAhGeAYbGFh/qo3mStF7Mw80D8TcwQC+vtLvt6u2vLmoI3frXe7tqy3u4KZ2Fok7JRTZpUEo/InACUnHxCz9VwYsRIM8UqROzckSE1aIJqArBRJnMVtJT0SKEhWBkyLcQTVmoAsCF7mmRtPR4QeUWjzhRtkatuMkRHqxWbSJHiESt218zubABQP6vklXU1szvT4pSNyxmrGf/1BxqVdYtXiG0AGoYa/URUheDSrpiLjEwUkynfTBNY1A3DYhXaersrXB60y7cXufqlAAAF6v42vXcGAOuAt0mY9NUKkTThKT9xBcg1QW36cAXVjZKYmJVMgSUr7DS0ziyCZTU3LPvtgtFKTqWS4KouX79Q/W4Ic348C1/C5UHpmktMuSwrknYtYRCz4FVOC6OYk2Ni6Xfr7a47X26zTsj3TUmOPxEtAJC0SoGRSY2IEMa3GrNq4hBaPfvZ2n/CLQyleohsOd2rmYITqe4WCEsB+juPAoAwZttXRn5Iq+AzUumo8OonBgu3v11xB867OwMugySZBLMHQi6Hyi6AQCpzNYI2W81ECC1VOf65SY6VJWXLZiIi19GBWwbOeIRlngm+iwoC/YHQ9dKVu+b09XcqwpwHvaUQloL5MOfPUwauHDtORA4iqpmzZtd0RFlrWnnzXjsRTX+l3Va86/KnPnfB/d8vd3q2tyjs+vU0g6PKalr1B7NhiWIxroSqDsCnnAUYMLBEmMAkxqT1ZCE77eRh/83/9m6u883WgrZoztk4cQQUFm37AjBYB7lN1tJKkqPt9DZ0HcCBLe7/AABjqf7PcpY8D4CwzDXB1TwA82yjIhmkdMVRG3yuC/dbixYDmk1GRC5bydJvMXMbERnd3e//m85QYJMkvcyqglBZA7/30mYAX2fm97Vr/TcR1R3u+dvdAFAz2x3s4Rg5Rv1ccnsPQVFOIwcLQveakBcrbpEdu1xu/5fnpj127HJ5hS+gNgNYD8Ds87AoWiQvAHDPm60FUaMA48U3wLvPHw4snbFcJ2tZswCgeHr5zwhu3hD6Hgzgwl1/zJKgDa/eHlW4zqpLKz9t+OziTabFbT3y4wDuhdYzWhaYlL4DAx5Dif7JRNsZD03P32EsWfrFm61Fi/8SUWaXzOyZt665AoCsqj5Ikj50qgygIrJStpasEPYrNZ0qaKuZ3VljYul3KnN1nq6Scv0GuPkEbGLehCUoxCWyY5fL7QCaCot1NgBynz9QB0A63xpA6YxhUYBE64oZ7bW6tz/zr6bPaLVXwzMwYxb9IzMPndWWXfizstmYJ9lYZVlYCXnVsozgpg71ez5xZNOy3y4IF3shIqd5tvGJA/efqKQvU2saQi8F5/f8bCdGf05OALWSpB8rc9iKIQWWm04VtK23u+4E0NSlyrZs/Ry539+i9CiH3NmiakISFOLtyYI2EgXP96osAYCrR0XpjPBzSLg71vYxWtW2L/ACEe1EDLNLZj5395+yKzD8u4XbE63OLN1PHQBuBrAjkTaPByJaDCCWTVdjLWuQhShVvF9us7att7tqPIGOJwDMtMgVzuae7z23rW9VLsApN/zjFVmwKo2R4HEzDBLBpw7rABLqjrWSTysAbNd6l/GsuA5P6WNtDzP7iGgPEa3WSm+mFG1os8Syf8CxbQvb5q1rHlZqK0phYhuA1mjXeLnN2kb0qSeuBc+fANHmEiKqZeZ9CXyVMYlXZE4AtSXlOnH+nB/TZ8no9/sAbQED4sgXi1xoMtCj7p3rkP903Ol/Kc7ha9wFTQCAmV1EdCxaJkMKqMMIi26iEa3gShTMAIatsRgJZr5ARAGiUf2dCROvyBoAbLBaJesZl1/IEimlJv2A84rvbQBF0N60WI3+yFpm3gEW7ccCSx/6heVTRNSIYA82Ll75YE/b3X/KHvbmx1LMmJmvEJGhcEX5h+Z9cdldkecnq7astsucmoLKjzTel5KZL2lCW3F3ed95DHlmyfCjxZ3jrxn/9e/v9K5bvtqwDUDD/IJzhfH0AKGFJu4+Fu3HA6hcJgPB3ie80GQiuemZu+zuS/0HwWw2FVkHhXqSITSttseuZE8yxuqRRuuhl03busgsV7yXb/iAEUPCW4kKLW4Xxrz89jYAm+bfFfnFqDDOyzkutwdXMmkCA9LrIK03FVrMAx0u4bkyAGOeOeHZcmgvTMXtWlV9zzeO5c5Y2I44eukxiFu00y0bvqSofcarXqeYZnAASfQQTIrqyQe2KWcAKDMWDdJ8OheaOAAIc7EVAY8f3qtuIAHRawILVn/s6liWO2PhJwE0accnCw4h2YRelwdvILQbY3Je9LSKjIgEEa3z9PPf55frJtPi3nC5TUtZFvwDCrxX3YmIvh6A1d3TKfTmLCB1tSkSWTDiBKBY5dno958OHUvKi542kRFRIYCVAHY8+83+ZkSpPJjGhSaDVrRbyrMUxeVzt75y9N/ivJ7DdemsADNMOWGLYrLV+gh/ZyFlw6d2Je1FT0bs0kdEBmb2xnoCEc0HIGt7DgEANEFNipqu7zz4atvQfcut07Mb9v/DjlwiKmPmc7Fei4jkeR/+wrm8WTWKbDSnutZH3DZZyGELoN4m5jvODbx8tsy8/uG0zi7DFyDKAZDPzCeHzm5CQXNoP5TXww13zb9iB9DCzOcTunGaIKIaAJ2xtJ+ICgDMX/WVf24VRst+xLA/eTw8Wu21A6hv639nnd1y0zYADZsPGqKtHIvZ/6dFI04nI7M5WWUK6ph5d6TIIpfOARDuAVaa9yges4WWPrS++2TCN00jRLQEwMXRVs1rhYL1zHwQuDa7RIQPKokCawJg7XDvEcWmZWEBDxXaOEWWtKWMqSxMHM4B6zwfwJWLqli6RoCIvo5JMizGCzPvJ6IlRIShQgv9OAgWAgwPq1H3wkwO13Lthu8gF/f9mJmJyEtERmb2JNLAlNfxB4CL51RULRfA5DN240YTWm3lHTeXzFy36gEADsXt2Zs7076t60zb69FKaKUIBwChqP2QJVPoWLKe8z4E48cJbRmUSpGFg9QG48RtYDqR3PaD+ss97RcOefv6TayqwtPdu3TZFz/9KSJ6D8l3tI6EE0DtZU+zKDbVhY4l5TkzM8/O+kj+Q4t6n5Ylw3LtmlHtvdFIpQvjRihsV59dXmLy9rlEwKcgZ0aZrOW8TWRt1gYALobq1xI6oz7nR6u99runv/TYo9XexkervVs0W25UHq322u+a/uJzlzz7B6V9x3JuJCkTmRYcrwGw1ePmQ9B8X0lMyZ4MOACIrNIiWAqmhY5NqEmw+aChrdX11gqf2vciInyMkb1NaHJgEcX3YHxiqZdIthLCmcdxOZGTNVz6iUg/9GCMS+emMuPOW0sF/9XykRwAn43c62AI9QBGEstokwMHAEE0qK7cuF+iZPVkbQDKk3StqcRkMQloFIEBmliGfCQWsYxrN7uRSJbIrgLIG+XvaavDkEp+/60fRd2ISzs+IWgZxGO5GOIVS/AlurYaLK6XKCnO2GOXy+1//B/35pxcaWHdKsN2AA1aKlDwJhOTaXpDQkR1GLIt0FBCNlnHwG5bsblORkTEYayZ4hfmHpzZ4zv9jzNtH56FOGeXCYsstHLp4D6fDQS5eok+/AVCQsuILHUQkYOZxxy+7q/cM8ultD8+XrFoceZzsaxFGIlkGP4hj/N4jcoMCaKt5oopMeG5k3UWAF+Iw0lsTURgQHJssjUIza6udYpCO54htVQhtiV1AGCcwCjEIBIWmQGGa9egEY5nSBXSGLPKSNI2+UpYCFlStkogZGVL6L6qIhBgEAhZUnasXz5DHGgrnsaz2iltZdYTFplM8rsFUoGyYFYW4NfhcotAgVSgyCRH3QInQ2I8tdpvf2q1f8vXFr21+8cr+h98arV/zBCPtsntuIdKiig8kgjJGNIadCS7snU5yk1LSnDqkC6gI/l6i1FOCjRBBWvDSsYqnaT/KwBNMQhtLoBjcdyyDEDMWcAjkfhwmesMOyRLSvTOXY2unW6PWqMdz5Bc6gFYfQG3kINRvFhjiTIzD3XGxkIxgI4xPzXWzRO9ABAW2iYA+NnP6Zaf/fxSHvOEpbrcSDgAiLb+PZhl+0DoWEoD8slYgJzUfLLeLof9N7+e82WLRVd18ULddpNR+mGmR0sqzsueU7V5hgoRWacNkzxHL2luht4uhx1A09xK4/qaGvPCAwf6/wpAk3Y8QxIY8Hc90e0758kxlMcckNc2u03ImZooyfRl1QOwzpplFK2tXtisOrm725+2zdWvR/7v+4Wl0/QzFmNIQP5rO+XRRotKACdG+XvKSeZw6QAgdDqCyxXAB1bb0NjoEitX2hwAfpXE+9yQEFExgKt/t3d2C8YXrhsrDWik+2UjuMV0wiSzJwunkyxbasWePf0oKpGUfac6bL9+Lf9Xxy6Xb9GC6Rnio5KZj0/g/WYBOJOMCyVtD/KQTQZteVbL2QHlknJFvtThDxhMJFfVDM/OyBAbRFQL4OR4AtWa76z+ZO/b6+Zk3bwNQMMYw+rQeyYtmzlpPVmkvwyAM6vMfaS7KxAotevkk0f9QOqKjFzXaIY7xyGwpj7l0sZcvX0RtJz+WKIDqSCpLoxIf9mFy+WN1bV6ufHPXvT1qlBVhiTRdbPuMkT3vX3DVobnvGBLZk9dE1kzJEbqAVivelvFdEsdkOb0q1RmSjgBKItrBVgFDu5TAIZqVUpm+Fse2OJveWDK22eawIbtFKIdj4vXbvXbX7vVv+W1W/2NP1i6/fnPz/5RPO4HBwChsoIIf1raXvCUr7s0mSWlqFSHi+dVePolKc83pwjaj3EdCK0egFVlVfhUH5CgSfDarddik35VcehI/vT66Q+/pR0fD85z/Qf8JebqyGMxOW2fWu23P7Gi+1++UfXus0+t9m9JxhCbMpFpxn0NgK1lxeYrs4uK1bNvz4LMRuD6sc8cAETjlUYc6WnGnqu7cbD7oAhwYFCP0XyHYm++Q9nSfIfSqP1/pB8uXNfiRK8Tc7McMuJ4Tlc9rU96A30DRp1tXKuoQrZcl/fsF0vMVdVIki2X0sTCefntbfPy2zfl9lZ13L5wudTR4ce5jv7Qn68H+8x5pKfZP9s6GzW5S7BsWh0qbZXKm+ffaNE235q/90O9UYfUSKERkYGIqg527Vh3one3ONH7PgpNFZBIAuJ4Tt/bN6c43zh7vE5bIBSAVweEQWcJ3TvhziCVtTDCdPf6DgOY+6mPzRRbf30UX//LRZAkYgBGf8sDdrnil1PSpfHPR//p2dtLbv9cobHQDK3umFFndH2k7KOPaptv5bxz+Y2fFRpKbQWGYjnfUIzW/hPCo7ptnd6OzYvons3apbwATlbnrt2GYFpO3IuFtdlo/9/sLm/F+I18zZYbtB9awp3BhKRI/7+XTz2r+FWX2SQrH7ttOv7z9TNAMFNzAaaobUZE4nsHv2tZmL2oGkN6jNDskpm7bytaP6MmZ4UsJD3O9B9DuXkWqrOXy7cWfqyCmd/X/jvAzANIzmLhqlBNtDhw9imXFKsoiDyWcAB+QkT22/9t3f+dJ/b+HwBH5lRksawjHD3VDUxt22w1gJ05L9jacl6wbcp5wbZC+//QXtkJQCkwlGCebTFMOjMwwg/38bfkqIuFteNjQkRlABKpYNnQ6T7hzjfOTuqK+KR5/Eeit8th9/u53ul03XJTXXmB8BQVSSzwg6eb8NXPL4TNKgDAKVf8ckVKG5JEtEqLrbFsHKvZXlH3GF/4pkiqmUBEH2DmPydwvmHxtLtu/ct5//URRPj9xhMpiHrdVIosMtS0d1+/WLrEooIlyeiaA1XR4alnm3HfJ+co+bnGrXLFLyflGs2hztZ/OPj3z/346OZuZo45s0ET2iCHbQoENh/AZWa+nMA1VgJwxhNQH41UG/6Dt3UmSIAKxdCpGrhU+sRHKpTtOzu8TUeuPveDp1PckjiIcLYG69763bW3ldx233eq/7ZqPNfRBJWyl0hb8JHHzEcTvEZcGRtjkWqbLFzS0+dVwcwAAars6gTgrCi3bv3Ffx6/74f/csBDRHNS3JZ4qAdg7VP6xBXvFey5ulusyFtpwuSzIZcBSDSYXY3YFwqPi1T3ZOH6XYsXm7G/aQC1SywK65QXQ8PjH98hCUEj2kVE1QnMjJIKEdEfbn1rLUMVVtmGgUA/lkyrBdHkir9qpQp4PPsojIA5GeXUo5FqkYX3nTSbdUKvJ/+Vq/7+vGlyeLbCzCoRtXz5Eysqvnrvygfe+tcHV6xdOnMbgAbdym+nxH/me+TMMBtJ/+TMUHGYQgAzPzzzlsKqmcWUbbTyqcvnyDqQB6tsBSZBTr0WZqoH4Hhs0W/OOvLufDiR6xHRdKSwxu2EzC4R8YPetf7Y9m3be16K/Exg1+P2d5taDi9fWG50exTRdKLDf1NtRZ9EVDNUaIFXnx8mEN1dnxvzAQWe3mEHUM8q1rDLOp+7s3QIyAKA4lcDrk+//pXPvn7qD1cAdHofPq0AaGKGtc/XJ052taK2sBqejhyFAzoXInxhE01EfNPa4+sUfcoVpdwy34VxuDqGQkQrmTmhCtejkXI/WVausy0r17kpK9e5IivXuWn7jt63iKh6yMfqV1VPNzoPt4scmy/9dBEAAAivSURBVAkrq+zytt2nbV6f/5uRH9IENixEox0fEU1gTQA2koQlZHMZpbIL4mz/WTgv7Bf7Lx22Pf/Rn3yUmRuZ+TTCExYW+y4dxtLiKoCYRa7rCNIoMI3wZOr8wHGUW+Yn5GuciIUmExJWioSZu4hoFhFZmDkUyHToJEnMn1GA3UeCC5azLAb5H3+1/WOPr/27FwD0A8DLf/2Nx0qm5diMQsi5VjNMer0w6vVWo5DrdaPP3sI/jEfx4sC501CZUZ49Ew7LEiD4HOoiPh8MfF/Yh5UltQAAIpBsUjzmGAWWwjwzR7fvkuhwn0KJuTJ0LBE7cRESrNM/FhMuMo29CJaWCiXjOQHUFuRaREGuJfQZZfnC8tce/+V25Z+/8VHlq/eu+iIUcScUg+zpMqGr14vu/gG4u7pFl6v/lluIlo10s21f/8ktZr1RqMwQOh3qps+DJElgjx5qcD+RoXaW80TXmVq7rVQYZEO4PYjRFhvq+kBw8rOh+96+hHpBIip9bNFvzs601tTOz16VcDFkClYcDiR7B+GhpEVk2pYqJ4horrY4YsSN6f3vfR+7j5w71Haxx2QvyhaQFRiN/SjRlwBqeBvp7aNVcgw8vWM7gHmICDwzA+zVh84fFDp5/vBLP12YV3lfZe5MDG1PjF8x3HM2Xm6E3WIXpabSuDNTNcO8DMCFNYWffBjABxHlWY33ugiKf18c542LdPVkYOYOIppJRIKZ2wK7Hh+2Mb1u5bfbArse31K3oMx0tKVT7DrUhgUVBci2GAFLL9A3bcQHPGiCUCo142LxAAKyGYBghgImlXttRwC8S7ndz0s5vfWBp886VFadvz79mx1vLnqpamh7QjPQGHAAECf7TqDSNgfugAd7ru4WA/6BdR+nO0uH7jA3UkSAiGYBKALQFrmR1mu3+oc9qziN/nhrZIyLlM8uR705BX1kzDximanArscbodkbqqri1Lmr6HF5oSo6d0eL/Prz29795cu7dr8TYd9FThCuve2MAVwsfhmKfiFCIv7q2raISYEVgNh5+pB/qX1un1Hoa3RfXRvX0NZ9b98WABvfv+IUy/PCppICYGvui1kNAEq1Y22Hb/fphrb1lOvIQMOxx+5/98rvdzHzxXjaMBIh98eZvqabrSK3scA4/R/inZXGStp6MiDoI3v4i1X+93/38V/VVuUtgPbjD8kvCzt0JUlCpT0fCP5gv9D91bc3fYJIBlBJRGFj7vXvPPLILdULrCa9PjQ8ChDMKO5w6e763NBAfHhoO9nZDntuoWwU+kQXXTR0ei59dpo+z4IhQxozt0HzSRGR/Y8XX9lcaCi1FRvL5V5/F1z+PjHDXGneuux3H1r4png1zvtHJdL94VPdosBYswDAX7x2qz9u90cspFVk/pYH7A3fXv4/7zov2lSVZUmiWgAb/C0P1EQIbUR7DQA4WGP+SOR1+198dvaJ8x3C6w+OBBQsMigGvL5b1t09eIIQmhQAgEkYMKegHEgwUS/nBVvbJ0vuuf/Fm166DaPMLpm5rfkOpQKAfMHThiJDOeZYsxO+/yjUy2CrqlwU84wFyJdYdKuw+kEpXcWUVpFB60VW1BbIO3dfwhpH0bClW5pdFtVeG+miRr1orK6wL8HwDNNhE4RokwIk6NXXtgA6l/OCLZYfzgmgtsRoT/nWOQK8plyGaO5vR3XOMjADNgmi3c8pLSKdbpE5AAghJCyozMbbuzqwxlEkJIkGvcWaoMbzpo3a+yXw2VhZAiDWfQtScf+o5OsgSQAMuuA2WESAxMHjyb5XJOmuUB2un5GXa8SqZYXY/t4F/5Uu795ELqqFmYZlmEYLP2nG/fDPxmn0a0jMHIjlg1oa0LD7JzvfDABMBJUImGu7FnAhCh5P9r0iSevsUsvtH5Y1+rEH/vD5N7afe3+0Pb4nG6Gg+96LB9eWWIreK7EWfn8cLo8JofkOZQszNhIN8hcqRNi68E2RMpssrT2ZZtwPe4v/d1v7awBytWzPSY8msCYAG5WAv7rEWvgAgCbt+GSigWjwQhXt3yktIp3WnmwsiKgcQGksewelE98jZ7YA2HjBdUl4Az5UZJcDml9M/+TMSZVWPhGp4ENJt+E/KszcTkQ9RLQOwDvMgxcETiIcbsUjTnSdwc32sBtuUiU3hkh1Kng00m34j4lWMultAGuIKCfd7YmG2+/Z/d6FPf6bygdpKu3JjZOFSS8yAGDmADNvBzBHCxZPGohI+uhL97271r6qj4jSvYPvpGRS22TRIKJKAAZmPjQJ2kIIZkRs9z58ugTxB9Sva6acyACAiIoAzAGwM9W5UGO0Yy2ARmYea3vmG5opKTIgvPf26sZffrll2fyyhxBjyCmJ918N4GCiG47eCExZkQGAa8ff2Z2H2w8vmFlgLMy1DioBkEqhEZEDwULBV1N1j+uJSe3CGAuTQdSvXTrTeOBkh2i90I3q2UXCaBBJrY0aWuUErae89xff/R2ApozAYmdK92SRCY2KP4CjLZ3wKgEMeHyHb/3KL+4HcGy0BatjLa8bmtB46PwZf5bR7Jo+rWhxgrHNG4op3ZMhIqFRyDpUzykGgkPmNgRz1+cRkTXi86dCPVCU7NlaABsCrz4fGUgPJzQev9QGs94gT59WZEGaqkhPVaaEn2wURiwax8wqMx8JFZpDsFbENCJaTkTLX921Z/P5q13hYjD+QECcu9Jl/Z/d+xuIqI6I6t4+sX/dnrPHxJ6zxyBLOszKLwUmqSd/MjOle7LxJDRqro6ToX8HXn2+oq3zinj/xGkQATpJQlFOtrhtSdWs0E4cgad3JFxeM8MUFxkQV0JjCKe9IK/WXpA3moAmLKHwemaqD5eJMGZ91hQlNN5wTOnZZaLEW7wlw/i4oUWWYWK4kYfLDBNERmQZUk5GZBlSTkZkGVJORmQZUk5GZBlSTkZkGVJORmQZUk5GZBlSTkZkGVJORmQZUk5GZBlSzv8H9AAt69418ecAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 396x198 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAACnCAYAAACB6MdcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19eZgcV3Xv71Z3z6KZ0Wik0TZaZrRvliV5AWyMbRyWF5ZnSCDYkAQTshDyAoQXHo4hCXtiAjgQTGIT8xkbsImxwQTMEmNsYWJsy1qsXbKkGWk00oxmlTRbL3XzxznnVtWp7pmx4bXcdp3vm6+nu6tu3aqu+t2z/o6x1iKRRJ7v4p3rCSSSyFQkuVETqQhJbtREKkKSGzWRipDkRk2kIiS5UROpCEmX82DNMxts28LZgEfPh03Rq8nlaAPfB9Jp/i5F3/l+8F341fAzxmNYL3jmTIG3yWbpVVxwVVX06pnIvGRfk8/TB7k8IP/z+KipiZ6MMdF9CwWeHx/LABgZ5XnwWFV8uWtro2PJ/OQc5Pp4xn1nZNxCPrqvnOswHcsfoXmYKpqfmVYF1E3j4/CLXHfnmuTr4ReC8wdC1yAVzK0qE52zvPrRc4gdo+DDpuk3Pbv/NACgb5y26cufkLPptdbORhEp643atnA2nvjBJ2Dr6gDAvXonT9IGo6Owzc30XUMDAMCMj9F3w8P0fniE3ldX03aN0+m1mm8kz4M5c4a27ThKr1l6EPzWRbSNuuks3zimv49eT/bA9A7Qlw30I/trV0dPJsUPVC19b87SxXc3ZzoNb9tO2qazl7ZpmUljbVpP28hDJze53MD8aqtr3I1pxsZ4jv3Recj1eIyONbZtCABQtZBuqNSFrfAv3hQ5nru28iDLTSjX+GQ3H2uQtm9sAGrpmtkFC2hbARe5mcfHo3PXxxgehm1qAgD88pU/BQDcfpiu3e0nP8UnU+hACSnrjWozGfjzW9wN4R3vpC8yhHR2wQL3w8PSRTUDdLHkB7Hz5tIrXwj3Y/ONY7pPwvTS+GieBQAoLFhI38nNJMIoIEgq+5mOLtgltI+/dk10Hz5OcJMphJtGP5R38BD8vcfp/4uW0y7r19E28qMKwmfopnLnLmPnssGDynO1gmj8YJqHnwQADG2mh7PxMn5w3/4G2i2Xg63nayXoL2OKjI5G5mW6T9EYy5fQa2Mj7PRG+q6XvnMrXyN/ztfQrWxynTy+PnV1MPwg3LSXzvv+AbpBr5v3UQDA7Sc/jlKS6KiJVISUFVFFjDzBY7Qs2Bm0JMDzAiTl5dvyk4s5syPbmuGz9F4t48b3A4SS5UeQlJ9oQSM3FiM88rwEL5gLy2gcRwo5B15yZQkUtOrpofdDZ+Gd3woA8JcvjRxfdHLL6ouTNF2PMIrKdXCS42220FKfO0ArzvQLGWEvWEFjnx4KxuDrIegfVpOAkB3Av4tD0pkz3WHd9W6eHZmjuz4ypsxPjs/nuPvqn+COQ3RN7x/4RwDA1U3XAwCuW0rX8vaTKCkJoiZSEVJWRDW5HLwTXQGSNM0AANh60qvM8Fmnx5gheiLtvHm0reiR8mSLfsdoYEZpPztjBuz8FvpMkEn0SHnqGWHMIBlMpod1skULeexMYCxppJAxGfkFMczhw/Ta1cPzng27hBAVY6yTsv5qfR6TPRtuFRHDIxugtLtWbHia+34CABjfTqtA1YW80rzxt+hyyGpVE6C1rEpGDDM5Bzb8zCAbTXPm0CuiYuvq3QpSaiVzRp8YbKwXi077iR3TcW/fJwAA/7f17wAAN27bwOPTuaH6n1BKpoSoxpg3GWOsMWa1+ny6MabTGPPlqYyTSCLPVaaKqNcCeJRf/z70+ScBbJ7y0axPlqW4pxhJnfU/Ogqwr83OJH3GskfA+TjZleOxLuh0RxmzpgaG0c/pj0rMsc7I93ZOCZdYeB9BOxFGfq+TfYBsMdsl5AKz8+YGriptCfu56PzkeHKuoitmMjDdfJ6PPA0AyB6jbQVJcT7ppM5yZ3Fus9GRYDxB5b6+yLzsIp6zrFLiQZF55bKBrSDf8TV2aM3XQ45x/M3fBgDcuINWt3v7/gG/O+tvAAB/uIyPLyvbwAAmk0lvVGNMPYDLALwSwH+Cb1RjzIUA5gL4MYCLJj0ST8zOnx+ctCzbx7vo/dBZ2Fby07nlQAyAunoAQOrR/3bbAgAWkrvKP28OT9iLGQvy6u3ZS5vw8ow57NdsjvqYnZsoNEc5vkhq52765yj5HO06MpjEfYaxceevdcsg//D6vXNHyY0jD11DA8z9DwMAeh6hbeZcRXOzb34t7SPGluzDN5szRuvq3HIshp7XSdfbX0ZGk3NLyfUSQzET3LhurvzARq4RAl+010kg8IWd9Dv+ayct93++8G/xL0/wgswuLSfyW08gU1n6rwbwY2vtAQB9xpgLjTEegM8D+OvJdjbG/KkxZosxZsup3qEpHC6RROIylaX/WgBf5P/v5veXAHjAWttpjCm5IwBYa28FcCsAXHjRKmtrp8Hbtw8AYPoZ8kd52WyZ4yJTbqkX40rcHeO8BHPECI3szDYSjg2WKVEPII58CRoIaosLKhtdxpBKB2jMkrrvh7RNHxsT8sVLVvOYi+l9NqRuiCHExlIsisXbumWTnfneUQoU+I/vx9gBWrbnXMXL9mXn0baCpDrgwCIOeu9oO8wpjmbxObkInaCwnLeEpdPxoIZGUDGu5PNfXPEzAMDn9pLq9sDgZwEA71tMhtMH1x+HbX4Fn/+ZyBxLnUNYJrxRjTEzAVwFYL0xxoK8iBbALwC8whjzXgD1AKqMMWettddPesREEnkOMhmivgXAndbaP5MPjDGPALjFWvsOfn8dgIumdJNaS2HBQUbHXnaLLKMn3F93XqDHCSqy4WOG+ClkY8VvZdePoICgVjp0Sv3ifiJEtYtaeIw2etVGjkgYTflpL7QTKhVOEfpWbSK9tnAehUWNCyZwDkLtNJgR1qNF1xMjTwILopNKjF8Qbs8vAQDDO8ZQt4FdWn/8u7QNu+2c+ylDK0rMjZYVnTIHnOCw53x2Zc2ilcQhpha5HqnQNRDU89XqcIL03dsP0zwfGKSw6Otm/D8AwE0/p5XGX3hpYIDJ+cpc5dpNIJPpqNcC+K767F7+PJFEyiYTIqq19pVFPvuSen87gNuncjAzPILUU9tc6BTz2AUloTrfBzg7yONXQUM0ktVtJWVNkFSedAkH9vUFySWik7awI1sscnmSJRlGkEMkl6V5ArDtFNeznBpYdTEFIOwKQgrntmLnfUSXEySV88+r9Dlx7Zzg2OH2AwCA7B5CzfqX1MOsY0Qajbq6BEGdRS9WuYzNeiDyBdjlvPowChttuQuyShAD0eQV63nB/9Pod9h+1fcBADftJUT/1qlPAwD+cC4lmFy3lNMOWzj4MjwceHK0+0u51opJEkJNpCKkvEkpBZ98pcsYJZawH090okIe3r6DAAB7gj0Cq3jb89bS5wqlHGoJsgwOwvQSGtuF8+krdmgb5Vi3yrJ3iJKpckhaOEhjpf8X5ZAWNp5PY0kyRlU0IOG8E7W1gTNc/JKiP3NgwX2/dT8AYPRJ0tlrN1EgxL/udwLr2oV/JSgQtbrdMU6TR8Hr4pWgeaYLQzsUlrxbOW0Z2+msUeQ36bRDUo910n87SHP81inSSd8++yMAgK9tWUljchg78EMHvlJthxidnFNEEkRNpCKkvIiaSZOvdAYlo4ie6fyqJ7phnyEkMLNYJ5XMckEjL6qbmm7ORj9DFrYZOgPbTDqv6L6mhHXrUJHfp57eRe/3H4XfSciYXkOWsr+YfK9aV3MpgjKvUHZ+kBTDVr54Mn70EH28j6zxsZPkT61Zy+i4vo1eR0cCBFVzdUjKiTUShhX902+Z5+bjkFTmI3qurAbijRBdUdICJSxcOw2PXfljAMC/7Kex7un7DADg3fNJJ33PCkJylxCkVgKTz0c9MgiiZ0hPfhsmiJpIRUh5S1Gm1cJfvz7+xVZC1NzBfqSXUzKzf9Wl9J0uhBORSBRb9s7POqsJ/tJltI0ghejAOgIjKCjbHSH9q3BoAOkNhEiFN1BMXeuZoqu5FEFO7PalNGN8LEAQ0aN5FcjvolVg+DDhY8NFbH2/83/TGOHIkCuToboriZPbBp67IOlJSS9kDwcn9cDzgrlrf6VcB104KYjLieXesaP4ykH6He7pI+v+rbNuAADc+isusxGdNBf1KLi0wNppAcpy3oRReQoTSYKoiVSEnJNSFK+jnf4RP2cv65cZD2Y6V33KU6aiOqntW+n9aX4apehtPiNJGIFV9hT86Oeph39BH2+n4sfxE6STVW2YBdvaEtnWZU+5BGrW3zg+LwnOzq+ayzmr2nyXkp1zOwgVsxx6n34JIwvrpAXOtpKK03B0y+mTnL8gieUuVXE+66TKC2HT6UBPDJelA4EPWooN5ZwYtX/8corf392Rxt3sJ72Grfv3rmA/qeikIwFyRiRsU4ivW5W1eM8cxGSSIGoiFSHlRVRrgUIe5sgxet9PlqLhIjisWgZ/OvnndFTHyXHO/+QsJtPGSLqYy0iqawKds8QYznI/Qh6Gs7to+/rzOE7/5tcEvkXxLnCysSvM43lK0rFOKEYm48bIbiE99tgeQqplb2CegXe/lfaRMg/R4cSjMHI2KIlx1vQhepWE6rmchxtKHAfgSkBMKu3q6V3ZtFwHmavkJ6jE6psO0PcPDd2Id8whJL3jMfJ9+4sWR8Zw10tEfx7KwBJEdT7n0Xiiupby1kyNZ+EdOQL0kEvFSl29hEebm2M3l3eAnOGSyGLPcpJvM98oM9nVJftpJ35IUtu20z/H6AbNHSX1oW4lhT+9ZZxoUl0TqoGPZvZbcd2IaiLbicEiSdr79sN/hJKrBzvoRli4nM7BW7sKAOBLmJOXeih3ka2uCcg0uK5JWF7sAk6wkeCBys4Hp9DZdDowPJWj3d1M/GB0v+UuAMBXdtHD99DQTQCAqxr/Cte00gPrblBxi4nKI3X+LgCTjZwLAIBvUBFRX8QQnkiSpT+RipDyLv3ZLEzHcWAJo8FCfpWklExVsNSKW6aDkogxwMvFGlITJL0uMJRCLhblDnJJGPuJ4qdwiNBAXGH2Lb9Nn4cdz2JoCOpIalwqikYiVhlC/iO7ceinhDZLX8WI+YE/ou9Gosjilm0xRMKJxGzgeEyz4zOSQiFpLCAg4cl8PggVK/eTC7gco+vyqadIffq341Q+csNSSnr+xKOr3Pm7oIGoC1KNWipV0pXfTIcZ5bRLWSWeOULfzS1KNxUdbtItEknkeSDlRVTPA+pqg6pPLjtxksvC28c66VGu7hzmJ3iehEWZ3UTQIR9N6IXnuc+8H/0cQJD0nO8jFMys53KXxWSI+ILANhRqFWe4fKf0PBErSPKFbwAABrYSpnX3NWDZRkKQ1AWUqOG7VDwOqYq+a6PuIqe79fc795MvKYpSMKgKGGUeRumBsD6Q4oLAqqj+/9WXbwEA/KCTrssPuXzkPQv+FgDw3vPI6LVNLwsYYCSBxLH7yZj8O5goWuvVDQjQ3q5eER1zAkkQNZGKkPIi6rRp8DdtiLt+Qs5pM8DWfQ+jynnsDpE0v5DuRV8ovcv3nV7nHyVEGztE76ddwJ4CYRVRJdASEAgXCOoSD6fvSbIIu1iGdtDOuzsJpVfO7UPmQxQSLUg4U5BTW+gqcdlwqp45eSpIsJHVR5AzlE4IhKxtmXcoDKuTqiUc+7Vj5MJ6cugWAMDruXzk5s3MmbWIi/HGx4JQtkZUrZtalQDEiS6mutqhcoy6slSYPCQJoiZSEVLepBSjkpVF7+MED69vEOghf6HzkzYqbs9SvKTMvOw9vRv2IHkKxjsJSWpamZFDiHQ1b5VmRvE8lCp1EQRJ/fwRGmIn6XF9/aQHLp1F85+5MhsghTtnNabjzWIfrPhVhdmkeaZjeHYriIRohVFGpcjFCxbzAKciHryOigbvaSfUFyS9uJFqN9+wkBPKpWQ9VJyoy0gc+o8H4V4Ajq06xhQ+POzSO52nQnsjJpAEUROpCClvZMryk6aiOGb3MwAAe3IQZgX5CQuXvCS6r37qRCdTnFS2vdsVx1VfwN6Fq19FmyidNCjfGIqMaRsaQskUUWQwZ0jvzT9K7H29T9MlXPpqRro/eRONUVcPG2KOBkLow5xb4fIZIJSqx3qpv3BhjBLdOkp4Jm+I+S+rIscyvb2OMulj22jcb/eSn/QTK8hP+pEfS0kzh6FlrBCbt9OJVWGiOwfexXkWmE/K66AVx29d5NgMg7myxyJJnE7khSJlTkohi1x8pdLJww5xUkbLTOdjjfnfQp4BIEBYc4L8rc5bAErTAwCzgnn4dcRHJwzLsbRFi5BueNt9NI1jNOfBI6Q7Np/PPts1hEqFcAmyVR4JXSYshW/yvVj44bJmIdaQngUpVdut9ezRKOfpw+88iH/ZT8f9Xj8xPb+tmZKe39omCH5BdCzVhMIWoTiKNMQA4B2mZBmvbyCyry/Rx1mzguQbtbK4vIEJJEHURCpCyouoviWrtZuzp87w07+cnjp//XnxdDFJc5MAUT4agxckxVkmU1u12LXa8TXBgiCaKispluzrIjxcvnF2K6HA6T5Ckrnn03vvL64GABREZwwlCrsEaNbbNJeo40VldPKZ8Tnc8kZ8ri5FT+UvoMDXUOLnXCaNHtJpv3GkBt/rp5LmN80k1qW7vk+6ut+6JnIdYmXkIbqkGHWmLhgUVm1FHwQpsBweDtoklfKBTyAJoiZSEVJeRB0bI4IJ8cFxrF3oH21VVcmny2VTcVmHPUkI4ldzLumFVNBnmxpj/P56zFjRnYjosAUf5hs/AACMbCf0y43QPvMvoWhKah3rv+Jf1H5XIChiE8IHIdcV/6rq/+T2E4SvqopHrwT1eFuHrKyb3vkeOsZt7YRwm4c+i3dxH6ffX8JoKLmsQmAsuqnQYaai0S0zNOS8DM4HKv0PhqKct3Y1d4CREh3x7HgmQF8pbZeGeE9uw2RS/gz/8SzsSm4Ps6SNXvVyDwQXT9Wc++1kJOR66Iap3kBOZH8lt62pq0OsylQHB7RLR4kZHsb4frqYne00/sq3kwPGv+aNNGSpEGJ4kXIuLXbkS52XpOjpOizHCRW6OX1VraC2dbwGB6nu6ztH6RbePESUtpc3vh+33cLjruFrpMKxOqEkFqYu5CNhbgCB0casND6z0sjD57L52a1lGxpCyeXszOKknMJTJRv2OUmW/kQqQsqLqOk07JxZrgZeJ6cAgOY79TYTZ78kqfg5+rx6PdfPL+NlTJbTYnVSpRBWwn1OraAenWO/6sWJw+Q6aV3DzvhV50eH1NxP2n02NubQzrFQz2CDS9xgKimnmNrjxk1FUU7cT4+8kxgAv7iP0Ok/B24EALyx6cMAgLcsTgNzGOUkhKm5VEsEU1yiS1VVYKwxQsYYWepVSYxyq4UbSjherB9QFfB4Z/GmIGFJEDWRipDyJqXU1sJfty5ktJR2UzjXzmFyShdOkU6UfhmnoF1KIVZfSiFMXDd0opDKIZjUszM6jD1B6W/7ds/Gmk3kZqn62O/RkJoDqjpaguHCjjKPs2dg2PktQQyfmwe7TUVnc2zZcdecZhw03czN1Umvtx2ife4foLp7adv43Tt46BULg9VGp9PJ3CWGoFxOzhDKZt010owsUkTorrFiG4yIuOF+uQcAcPopuv6NlzfGt1WSIGoiFSFlZkqxgPXjOlmIazT1KKWiCaufP0RPXXo167ULWCfSnTvcIxd69nQamSC4WKz3kE56+jFyQR3uoGOsWN6L6k3RhmwuGYPFIasgiJSPiIU/Ph6U3Ei9vCRfiw6oCvOcyDxTaXc+4g76j7e2AwDuPELbSPcRqbt/9zJGw1XMUzp3bpCwo70r4iZUq4QOTJiBgYCRRZBUthUdXZJgxAUlSSyiW/ecAp4mRpR8N23TeBWXul+8DpNJgqiJVISUGVENIp31REKWuu0kPTF/lPyYmSvJke9f9jL6XvPii59VEjyAmPM9EpIEXMvF8b2EfoKky5aQY37aP7wJvpSPiHjRRA0ncgxBECFTaKh3aXPaKe9S4hTna8xzkB1xZSPmOCXf3H6YxvrJ0OcABN1H7vwPTv9bwr7SUJJyEG4V24BeTKzYsCoyX2fpd/c4RhYXJFDhaDdnGVN8tJJM8+B/I/sU6f3VV3Ja4csvjhxnIpkSouqmvcaYVmPMVmPMdmPMbmPMe6YyTiKJPFeZKqLqpr0nAFxirR3nXqm7jDHft9Z2TWk0SX9jFEz99GH6uLMfuaOEepm1nPI2j/U8zVYsVqWwOXO0I9JHUFn35l5ipxvbS2h9ZC8do62F9L/680LpfpMkTHgdFE0xXewrldKQBRShsQ0N8YQaQX2dKKzKNgTxMTaGHe+i8OI/7yXLWJD0tY3U3fMPlvAxFjJPPyeBhMPEpcKtLlVvOod6NcWQ6NTNs4JrGOqxCoR0VGHedumX7Cu950EAQLZ9BFUbWe9fR5yq0CVAE8ikiBpq2vtuANcAgLU2a61l3w6qpzJOIon8OjIVRHVNe40xfcaYC621TxljFgH4IYDlAD40JTQ1oBS609EnN/cEdecbbjeYfiHrWldTSbMQeDkLVZBUfKDuTMLx8aCgDIDjMh15mpD02CEisViwgKI70//+KjoEx8BtTU2ge+aKo6DpoDnbduYr5cRpIRGL+HIFWTX9jc5BYMvaMH2P6e3HrQdJ17yj+5MAgj5Ot9/LPsnVTG4xXfkiwx4Pd7zRyPXQHbxda3ehD2LSCzu/xRVguiQUKYmR47mUSfbVMpL2PERjz7miBvZ3mL2by2rAaGzrdaOvuEwFCa8FNesFgqa9sNYes9aeD7pR32mMmVts53B36d5TSXfpRJ6bTIiopZr2GmM+ZK21AGCt7TLG7ALwCgDf0WNEuktfuNIil4X34KMAAL+LEM0fo+dy+ktr4K3iJGpdliv6piY4cP1F+fuxMae35v75RwCA7n305HYNkR634XzSn6rXiJ7FKKB9s6H/XecW6Tp4lv2GS1g35DKS4MT9gOFajeUoGcULIHFw5gm94x1k4d/WfgaPDlGjRKHZ+dMV3D92AaXTCc9+jBAupKM6j4RQVrJHwBw/zvvkgjkjRB8klJKjI0ADU4OqvAmx8r1dFG3yt1LR4zh3xZ77Grb6X7YWVvcQkPnUFMmeUzIZokrT3lZrbZu1dhGAI6DO0rUAYIxpAumw+yc9WiKJPEeZTEe9FsCN6rN7AdwCIMcoawB8zlq7c7KDGd+HGR2FPUU6auEU6UrVl3ES8qteQSS6QJGSi3C0BoC4U1Vc3wwPu94APfsJSff1EtptXES6V817LwMQJGxLj1AX2/b9IMYun51hNOCCRCxkgjXJg5WsoVABofbxBoWK7KkQRD3G6j132/7eURrr0aEv4bLG9wEAvvLNDB+P/ckNUfZonUvqpJAPSCvEUpeIGCOp63Q4j6N/irzOjI8FLNjjyiMgeuZuoq4c3kHfS1ds/4+oK7bp6wv0WCmvUfkSE8mEW0zQtPdLRTafXMbG4O3bD1tDh02v5GWrhZ3Inhcs9fLjitHk3EX8I0vITqoeeb/83Y/jxNPRG/S8ueQ0n7lOGpOx8i7cR3IzhmqGdOKKbWDDQ5oGyyurGe6BCqsPstS5h46PJy2Hemhe3/wgfX53O+37w0GqFr266Xq8aTG7vUS1qImGOWMPsFckzVGMJlYBPKlvkmQYTnoWuncR506bVh9vX7STVCG7sx0AkH2GDNX6yzj9cm0bXQYJGuRywZzl+uuq3AkkcSslUhFS/irUkXHgAuKwd7U7siSl0yWbcwVjiHHF6CguHUar04fS2M9IevFSMkpm/A2zrjQyn5Vb1k3xY3hegFRinLSy20knqQgq6/fZLMwII+rIWPS7ve20CzcEvuUIqSBiOAmr3ne/mnU1SK6FO9f1x0pDfJUELWhYVeWWWkFSXT6ia6ccOg9zok0uGzj2hVF7C7HbjO0hpK19aZS9255WtVS1tYFLS9Qkue5JG/REXihSXkRNpWCbpgOi70mlYqT1ozI8NLO0oJLof09we8ouej86WoMLWsn9NH1TKnKcUhz2Ot3Q5LKQRmvQij9kV7WvDoPmcgAbYIJg4tL65McJfX7WQ/s+yoV4Yjhd08bzba4HuGlxjJVZJbJoXTVIbMkGerYgmaTqSRWsY9MeiV4X0bszdUg9+RSd3lZiRBk/wEh6MafqbaTAgwuy8Apnw0WJwkwoCCpzTxA1kReKlBdRqzKwixYGSFpdpIxERD7TM8xFC8H8PkKp7ElCwHkXZ5F6yyU0vqS6aR20hPUZ6KW5AIWra6PbuLKRKAeWFut5MMITNcbocppQ6Nu9FH7dM0R9nS5vfD8A4MEPkBcgfQWz+S1ri7mUYuJHvSSxflPZbJD8wdfD12mWRZhRAAQehtppzro/s4XGql/J7rLfJ/ZCnarnwrNhjliH1ApBNZ9WEUkQNZGKkPIiqjGEZs7np/Q8wCX1uqc6H/Wbeo9R2psEDaTUNjODLPjUrJrgKdahOV3yMaKe7HCnj1J9kwrKzytphooBzwwPw5yghJV7P0Zz+8Upsu73DNwMAFjbdC0A4LeayfeYWsGoKIyGtbVBeFehvi6FEdKKWNJOPu/KRmLXO8QMAyBYrXil8bY+DQDw9xzD8A66VnVLad8Ulwb5ck01b6paNU0uF/icRWcu9vuXkARRE6kIKW+5tOcRQmi0CnUFiel8Yokzsub3cPn0EL2v3cCW66uprMFW1wQWcgku0YBTlFGHkcVK/9D6hiBCJkgporlVpWxDoYQ5eQr+fgqNfpoMZWwf/BoA4H2Lien5pk8yJ9V05old3EYbStQrHZq/IOlY1CfrylqkjERS6Lh3LKY3wOemvY68gjmmHPO1IKv4ORkFR++j9I3OZxrRuoHDvR9+B5+++p3CvnAgXjadyQTXWfpr8fG8/34Ck0mCqIlUhJS5uE+JStkz+XysAM9+n2hfCr2EJKNHaZ+6daz7zGPLXvx1mUzIqudX3VHEla9wIkdVEYIGneRRrJ8VQqjcS+jo8fvtHzuFz++hyM/2QSKH2DjjjwAAr5jNaPTHpR8AABn+SURBVDiTV4OaTDB3AFY8Dqm4Pu96rso82Nr2xK8sLcWnM8VOKHklYDlUP7v0RH2cfKWFre0AgKEThHiL1g4hs4F1UpnPKFv5vGrFOGYFPUVHlVKV0Ll4O6n7tu3qx2SSIGoiFSHnBlE1OVhYR5UYMWcYnX6CkHVkiJ7ueVdyxOht3H1PLNqwPlfML8vjR96KXiU+v3DfgBLZSC7TRxBMVoCjXMzWRdlctx6cj2/2UPmIkEN8/aOUe+CtYhSazxEiQfhpKu2tmDWsCgY97oYnmVh2FumjriQmxJ6tS2Ecq7aUj/+ICNeO7iTdddlrWO/+wB/DynnLNRR9VvzJktIor3Itw9lv3BXFeROepgJJb82C+HkqSRA1kYqQc9NnSkT7IvN5FL5BHfEGd9MzlB2nJ3fuy7j7yDomSUsp/VPGSqepYCb8nYwv4hkUFZ09FJmj6h0Q6lcKAB//EKHFXX3Ug/7AwL/Hyke889sAhBKUhepHokpiDQvSZ6oCBOMXR/so1j2fo7+YUUmzWBfyAZKq3Ib8Z74NAOg/QOjXf5rQePVv02rmXUFUO/7oSNzLoAnVnI0g3oeQTgoA02phjjBhL+dReC+n/ra+ZKZNIOXnnioWcgxl8WdP0k3VO0DLz/LL2Vi4lipFfTEOdLJKeJl0P65FRIpxpyJ8o4QDD0o9kcQJcUdJokkv3TiPnKJjHRiisrGVTW/BzTdwa6FNlKrnL2I2P3GfSca7Pgd5tX7MledClUwvLpyrdvbc6LmHRdGry5Lfs4cersN9ZJCuW8iNIv6EsvL98Hwy0QSiWGBBu6Nc8IQTUU52u4faLllEr5K6qAMvRSRZ+hOpCClz4rRPyORqiQiF/NuogcRIu8Wxo5SQ0bqEjJJUG9c1qVS8APGYqU6W91w2HlDQiboa1XWydiEfGE1SVyRtIIdoyb/v/YRs3zlKbpmHh8gFdWXjBwAAV8yaDtPCLjap5hTkLtWktsgq4R0hBjyH4CxOfZAVRjPzyTUwQRJ46kFuNLyXkH4sS9dlfSsZZNM3KBeYcNRW17il3JNWl+LKE/VFM9cIknZQLZUZOgPbwvxV3KbIiW49WUQSRE2kIqT8XVFyWRjRr1i6nqQn+GBfEy5cQWHHug9fSbtwahp0dWWen36j9LswWrlwXia6b4lWkxEXC+tNxrGXEOpIpejXD9PY/8lMz9fMJhfUt25mZGn2YOe10f9iHIlu6lLxSoR4QwWDwixtjzDT9Aoympz7Sa0eSiunc2IXWv4JQjfXaPi1fL7XvZn2rY+iszPCCvmg+4okRLO+rdMdZR9vFznzzX42oJYthF3KrX0cU4u0zZz8NkwQNZGKkLIX95mxcfj3UacT6RV1dozr7hd3o2E9W5eiA0k4UTvxS+mZnhdHLAk/SsII7xLjk3Kh1lysMO8j19Fcd/TTmA8M6u4jrJu1MH9V88yg/FgCCsqqd/NRLjaHXkOngTNcWt3G+l1LlDmpFD+phDS9Bzc7bq8Bbt7rGg2vJ4QrcDJOjLtf5jU8HDD/STI688c6pulnuGnvXnLP2SG25Ne30aktnB+sZJnob1qU719JgqiJVISckxDqWAeh09AJekpXvZJJCt75asdxFIQKle9TUEfe6zBjKhX4J7Vea6JIFpMw4glJBDNIf73/YQBA1+BmAAHT8/23MUosIg4qn0vAMa026NckYd5QqDgssbIWPjb6hxw5R2FjtM+VE19Z++pYhd0n0P00fTZnNaGc95e/Q98Jb1VOpTKKv5k5sUwuB3/RosgczcjZyPG8g6T/ZreRrzSzlH7HwsteStsPnw1+F10Ko1kNi0iCqIlUhJQXUUdGYbbtQe4s6Xm1dYSsXhMXkWUy8cRkzS2l/Jvu+1LFb0DcqixEUciNIaUYo6N45kN7AQCPnCBEEyRtmXE5AGDDTEaFZtapmdwiUl4tyBQj09AeCjm+lLPwelFXC1tTfEWJCY+d2radjr//GACgfzswfRZds/Qi1jM1LyqLDSflRD7PBOei2sOn9lIii3+U/d6zOUI1qz42RbeS5aJlO1ORBFETqQgpbynKaB75Xd2Y8UaOqmxkah+xjmuq48nUGlklrS4fTTSJ6aMI6bEaORXDsmO141Q5nD6D9z1OPr6fDH0KAHDremJ6/uObOb1NYuxCE6S744VRXKJo4yV6nqo8ArtgQTBmKR1UxI+i8djdlELXsZes84VtedR+/A20TaPqsicicXrF1Cfbh6+teBNSP3oYAJDfSzppagX3Dnjr6+mUVX6B66SNuHfB27MPk0mCqIlUhJQ3za86hdSSGUALd8cQJJ1WBIUEKRT9i4vtl7TcTZEolup9KmNJhIwR9Ud/SYh655FMrPvIFfM586dxDc85pFcDcXqi8PzEgxHuLQoEyC7bqutgPa8kSZwry/4OsWqPb6NYfNeBaFfs6lX1QYJ4KTR2RXesh8o5hFDR27GDjtdBkcPCYTqeIKlpmxfdR8T1XghRS3LehMflO+IrnkgSRE2kIqS8Vv/0ethXXxrEkFU5A/wi5dKCPhJ/1hGrfLRQD8aLewzEmpdeSTqnlMnM7uD4/bd7Px3E7r9F4/qrz48e33kbFD1REd3UeTLkVVNVyhhSIBei6YkxX6s80OHNdA4HnyG9f/U66Yr9Vhqqri5OWaRogEQcouociPExYDdFnPKH2LpfxiUv11zN56R8sfp3DHkSvC7KW8BxisDZ1hZMJuVnSqmqimfYu8zzIKnCbVGlXDzupuYtdUvDdDq4adUSK1nxQnN+z/to6f/aIbpxZLl/W/MN+IMlfDNNZ9YScelIYol2ORXjKlBtGJ0aI/syJbsLJSomF+t5wc0k5/LV7wEAzu6i94eO0PxWr6cbtOYl3HJd5hs6nuMEEBVIwqLykClDNvUANZDzD51C9gTT2F/MIdwlfHNpxmttsEq4uLomaKYm9W3LOIG6WbXzLCLJ0p9IRUj5ETWTmZKj1y15qq7IiW6IFkYep8ALghH6mnYyBGw/Pdn3HyNDQBtOd90J2AWcxiYhXVWK4VBK1AiFsNbzAh4tqXMSroJpjHaCQqIK5KJBDFOVCdINed/Bpwip2rsIhZa2kkFS/VFagp0zv5hhKtdMlukw0x4Q1GXx94W9tEQPHyig/jxC/8LVr0FESlT8OndaVVBDZQboOkiKoEP0Iq5FLVNC1CJNezcaYx7jhr1PG2PeNpVxEknkucpUEVU37R0B8IfW2oPGmBYATxljfmKtHZxoEOsZesJ0knP4idaPjiSSQBXbibgCPk429v0AMc5EHfo3XU+ocHcXPdlPcvnI29lwunoRz6ehLjBwxGibbBVQCB/Wwl1wIhbKDTVxA+JFib4P3PsQAOD0k4SoJ3rIkb9mA7cUX8OMKBK6LRb4kOsh+rXWTVlSv3qcPuck7cIozbvh5Q2Qxm8FVcTnbApx7Mt1khDrvv1uXi7ZmpNhXEWt4rwtJpMiaommvQestQf5/y4APQBmT3q0RBJ5jjIVRC3atFe+NMa8BNSe7NDkQxlCFV3SG0ZWnXPh6tPVe0mKFoe7azbmA9wG0nRx+cgQPd0PdNHpPjl0CwDg4sY/AwDc+Vn2BqzgJ37+3MANNa24HhVreal16FAfAIcyWicUpFW8A0ZYWPIFjO2l/090EQotWUNzrf4r0hVdcZ+gpV6tCvkgjCnuL5UcYkY5dHuUmRIP08KY2cCpi799VcBIowMOKoQrvbqMOPN72AXYusD9VmaQe4OJ9a/L2ovIVHTUok17AcAYMx/AnQDeZW2xgnJEm/b2TqgZJJJISZkQUSdq2gugAdQG/SPW2l+VGiPStPeiVTaKmjr9LfS/YkOOdf+QRJKzTMjg2jaOuWTnWz5MT+xPughJfjZEHfF+q/GDAIDXtbBHYWa0O5+tqoqHEXWShSQ7u3kXIt/TZ/yt7r6ieVoVAUXhrkcBAP27Uzh8klBt/TrSG2vWcPPcWpWyp8k1XNlHVZH0PU4JfJBSF4VNL3eYrltmHfs1ORncZqriSd/awc/X33S202s/oaYrkZ43NxR4YZ08zXMfC4VXS8hkiFqyaS+A7wK4w1ob6yidSCK/aZlMRy3VtPfrABYCmGWMuY4/v85au33C0ayduOwghEYx696F/1ifyUeR1flKT55yCHH7MUK0x4f+FQDwppnXAwDu+zJHipaxxdzEuqkwPVfXBLqvtvZ1r1YRTceTTpfu9qF5/8VjwOfSs40QcFd3M1Y307lMu551Ut1tpBR7d7HPJALGyc+FndQGXUqDapZxSPvNdCxfc54CwWoY6hsLAF4HJVCLXWDnECr7a9fwMYeDFYaRVBivJfl6IpnwRv2NN+1NJJHnKOVN8/P9wKIFgicspKM5C1kz8UlqnkSbmFrHHOcEhzw96t//OHBXO0WcHh/6DADgpY1/DgB4bQufbqPoooKaPA9BpzBrtYheCUr1Qg2XapdKNxRhpPO/Tyr+6W10rru6ial6zex+zFnF/kqx6kv1E52o3EXKVJ7cQvsepPLpLJer167gAr3F5FkoFCvr0Z0U+XfwVP6EUA1hZlN0eyAWifOOd/K5KIKQIpLE+hOpCCk/pU8YVWLRnnSg+7k+U6qPkaTqcSKubacIjX+G9L3bD8/H9/oJSd/WfAMA4Fs38lO/jPWqZkJc1EiPJkkdVCgVOm6MUVUhfXBOoS1LrQqOb59QaWALnevWDkJS0UsXfmYDbBNTGoluGukbi9LFfmEWb7lWgqQ7OD9gPZNJvJE0PL9ETyv4fuBVYLXbdFF3aUmVlGvqr1oZnUfIWyLeB4+J0wRZ7RT4URNETaQi5Nxw+E/U5dlZwBzpED+pEEKcZK76bgoefPSfCIX+q5+oFLcM/qOz7q9tYwu9hQvPGjmKU6ciNMVKtLWXQYtY7HoF8EW/KwToKtuKPvsAURoVuI/r6Cjt47piL8sH81S0RE7keKUIxnj+3rFjQdYSf1W1gel42ii31OrO1UXKmc0JyjwTG8Hw7+HoL4XMTraX7KlQdxmXl8sZaWKPhDu3lJJz274nrX4E34//qOJIFha5TrpRC8fJxfJf/XRjbxn8KgDgohl/gnuup5s2tYZuYruoJXI8p7zrmzDseC+1pIrohy1XRBUQjlCVXpdrp6U/x4G6uRvphs28jqsIpB6rri5epcA3v+PPkmOJ60luMrmhevuAYf6fQ8T+qhW0rzxIJXhjnYswl4UZYH5YUSOaqELVX7qMJ1KiWkDceSEj2p+vMvoTppREXihybhFVJMQLanSlKC9bGCR31Kf+kUJyj58iZ/GWwc8CICQFgFfPnA+vgRAVaVEpxNGs2p77jFKCBtK81k8HiCQhU81tqubutgujsqwGW7nx10lCpfwZwsGaVkJcbxYjaDHmQm1UapeXQ3Z+eXIr7XaIlmp/JAfvIkZQRkGnTuiEGl2ftYfYYkz/YJBgI/z7Wn3KRUtSvKPt9J7D2XbObNjm5qLHdXwCE0iCqIlUhJQZUQ09gdqICT1hrvZdkIP1KwzTU3dHL2UTHhq8H0DAqve9DxKCpNr6YUQHq2Wk0qzVpdBRpFh7H51eqIoOXfJYuFSFkcI/RC603HF6X7uey1teuSl6DNf4lo9RUw1ki7AVItABnU4qLCsHKSw6vpNWoqrFNfDXsMtIJ65oo0m3z5Qmb72nIfysftsS2kajsaxSopMepd/D1ew3zwqSrKU3gPQK4ODNRJIgaiIVIee2z5QOP4ZDrLzdvR+mp37bAKHQoQFC0mVNVMz20mZ6KlML6Hszd0agiwlyCo+A01UVm4m8hpJDnM4p6OZKUxjJ5BwKKr1PkjS27YS/mxzsw7tp3NpWGsvMJneMr90yOvk5laaUdCAoyVGJIubQYXrlUHK+k7whmTnscZjX4NDXrVYyV51Y0kUo6EmK3hnW02fVuybApRgJzSlmyRZuV9FZG9kV1dDgCiGFh8olUOtCzSKSIGoiFSHnRkcV0YkdQMCoN0pP2YcOE+dR+8ADAIAblv4dAOBTn2XeqDoeo5n8ebYqA1vPOqBDzmjSSYx9RPd3ymUDK9/pjRM00kXIV8regsJTHej+Fe0zZxMnhbzr1QAAX1Bayjt0575wGYmgXlrpgOIV2UOImt/HrHpLOSz6+ivptbY2SPLOxh35kbE6aAXASQ45Lyd/Z2HTxmBj5fOUfb0O4mOVMiC7jMKifttS2m446IQjPatMD/N5yQo4gSSImkhFSPmTUrLZoDw2xthXwObrSW95rI9QUZC0rel1AIBNTbzvDAmHckmGEEWkQ+l1rs+UQiVdViKl2GFvxCQIGvNvPkIpdIV28pX2bMlgWiOhrNcs5cnRHqiR1u3FjpVKB3NmHdXbu4fecijZserNpmOYOcxpWqQkWlqny3l7R4hPSkp3BEnRxHolc8CGvQWuC4roxsLSN8Sc/nOEqZH3Zc+CrZ0WeCi44bELZYfph0pIgqiJVISUF1ELBWofLpYhI4rp46D3eBZ/t5us6M1DnwAA/MMq0kk//CVGv0Ucv5enUFHHIJMp3ZtKpER3ElfeUV0dY3I2ukCPuVUlSWP4EUK2zg7SEZubzqLxr19C20pvphqV/CFjyq+QiqK49bw4l+p26o2a20fol17M6Pd2Ynr2dcvzsLiVhM9pfzvt28tJ6IzG/uWX0udhggqdQ7CP+0mN8Qq3YTntu3o1fa/0bnP2DAx3WbHMZB3rFD6BJIiaSEVIeRE1l6Oe7Yph+N9voGk83uth89AXAQCXN74fAHDJLLYWa7m0QXe90/1Oq6oC617rpKW6S7Pozn4Agvi/e88IfhdRMp7eSYhy/DgTRKwkpKual3I6mEvMVn7LkhLSlb3NlBIoPtlcN+m9mdWE0mYh00zqQrzwMUQnPcAFeOJZGeeS82W8SjXFOfsBsti9h35J353m9L4RRuXlTFIxk5mnw93+AHi9fKyqjOv251ZDyfg6M1TyUrjTmXSLRBJ5HkjZdVT0D8Ge4EReRqsvdhJK7B74Jt49n7qPfPVrrJMu4CdWSjJUrN2hU7i8WWc06aplRbOoH1cb2teJeCbYYzHI5SPbj5Kv8eLVFNWp+dvX83yqgmyoSbwOwYEVIYXvw99P455+io7beCmhkf8O6gjtdNcYUU1A5OGoMY9wMR2Xk5uNpFe6roA6bh+y0vM7KSPN5mibzOXkty5c/nLaRunw3jNcqiLl00sWud/OqKRzI6g7gZS3fc9YHoV9PfjCHZTYsG+QJrp7gFrkrGt6B17KmWColsYHJVw48rk2GvL50I9WwsVUjEUwcgwEN4Dwk37357TpAD1U3X1kCKyYRQbCtDZVyZpOB2FPXeekRd2g3i5yQaHrFMY76Ph1PL5ZODOyqwtLFqv3AuB1dDgeKDvABuA8Nu5mzYgc16k+7ILyOM0PQ2eBeg4aNLPxNn9u5HxdBYBUpbL7TEKoqK2NhU5duFxUowkkWfoTqQgpL6KOW4weyuPmk1THfnTwpwCAd83j5f6GLnhL+OmexdDqDBFeRjU66aRn1SgNQJwnSYQ3jS3zQICk7JwefJTrm4ZpHkvXM6vee66gaWgGE6B0sEA3fVBtM+2+DjqVAwPILGIX0e+9inZtaoqOpQMTgnDimO/qdix9Zj43iLiUGunqchG30jAq2v1ULWqq0jCvuoA2Wbs2uq2rzZLwK7P4NRNaFy66yM3PnGU32ElKYJEqYJ/ZvSeSBFETqQgpK6J2DFbjL36wFEcH7wIALJ5BHEerZ7BulAonrKi2PIo7P5YiGBaneyonvUZhdyzFFOhbeM9QiBBdnEDhUdpv43xCWkmjC46pCuXCrC9KN9VOfGdMnAl4UQEg3TItSAl0qYol9G1B56Fowgf6h2AaGOXrVOshlbLnCROfuK/knGZOD8Kdcg6SWMLNzTDCSTlzuOJXwq9yjqeHg8JN5baLsOeUkARRE6kIKSui9udP4M7uT+MP5hJn/s1vINSqXcbJEfPnAA2SqMBPXWYSXqJn0Uo7ZuVrPqV8wLLneOyPEGLMeDPXwG9aGxnShUV1mUdYdDmwXg04KcQMcLLGGk6RW70qYC/RLdXl+GLti/P8OFvfPawrzmyEv5pKc2JMMBIeFjQ+wbojl//YC4iJz1+yJMYfYE7Q9XGs3mzdFy7cFJmv1PKb/n6XGK51Uo89BhNJgqiJVIQYayfnT/+NHcyYUwA6ynbAiaUZwOSe5uefVOK8pzrnVmtt0aYlZb1Rn09ijNlirb3oXM/j2Uolzvs3Medk6U+kIiS5UROpCHkx36i3nusJPEepxHn/2nN+0eqoiVSWvJgRNZEKkhfNjVqJHbKLzLnVGLPVGLOd5/2ecz1HLXrOoc+nG2M6jTFffk7jvliWfmPMtwG0AHjIWvv3xpiVAGy4QzaANZN1yC6nFJlzFeg3G+dmyrsAXMqNk58Xoucc+vyLoMbO/dba//Nsx31RIGoldsguMeestVaImqrxPPv9is2ZP78QwFwAP32uYz+vTvT/o7gO2QD6+MI5eXYdsssmRedsjFlkjHkawDEANz6f0BRF5myozvzzAP761xn4xXKj/lodss+RFJ2ztfaYtfZ8AMsBvNMYM/ccza+YFJvzewE8YK3t/HUGfsHrqNwhuxPAKVDdXopfW0Edsh8G8JnnU/PhieZsQz+YMeZroJvgnM99gjn/AtTk2QdQD1q5vmKtvf5ZHcBa+4L+A/CnAG5Rnz0C4HIAPwPwgXM9x2c551p+3wTgAID153q+k8059P46AF9+LuO/GJb+a0Et28MiHbIvB3Adu3u2G2M2xvY+N1JqzrcAeNwYswN0E3zOWruz3JMrIaXmfG2RbZ+1vOCX/kReGPJiQNREXgCS3KiJVIQkN2oiFSHJjZpIRUhyoyZSEZLcqIlUhCQ3aiIVIcmNmkhFyP8AFfhj38UqOjMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "curr_fig= plt.figure(figsize=(5.5, 2.75))\n",
    "ax= curr_fig.add_subplot(121)\n",
    "ax.axis('off')\n",
    "\n",
    "import colorcet as cc\n",
    "subsampled_color_ids = np.floor(np.linspace(0, 256, q, endpoint=False)).astype('int')\n",
    "color_list= [cc.cyclic_mygbm_30_95_c78[i] for i in subsampled_color_ids]\n",
    "\n",
    "# Plot embedding with color\n",
    "for i in range(coords.shape[0]):\n",
    "    plt.scatter(coords[i, 0], coords[i, 1], color= color_list[i%q], s=30.0)\n",
    "\n",
    "plt.plot(coords[:, 0], coords[:, 1], color='black', linewidth=0.2)\n",
    "\n",
    "# Plot Pearson correlation matrix\n",
    "rho_frequencies = freqs[bin_low:bin_high]\n",
    "\n",
    "freq_ticklabels = ['A2', 'A3', 'A4']\n",
    "freq_ticks = librosa.core.note_to_hz(freq_ticklabels)\n",
    "\n",
    "tick_bins = []\n",
    "tick_labels= []\n",
    "for i,freq_tick in enumerate(freq_ticks):\n",
    "    tick_bin = np.argmin(np.abs(rho_frequencies-freq_tick))\n",
    "    tick_bins.append(tick_bin)\n",
    "    tick_labels.append(freq_ticklabels[i])\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "plt.imshow(np.abs(rho_std), cmap='magma_r')\n",
    "plt.xticks(tick_bins)\n",
    "plt.gca().set_xticklabels(freq_ticklabels)\n",
    "# plt.xlabel('Log-frequency (octaves)')\n",
    "plt.yticks(tick_bins)\n",
    "plt.gca().set_yticklabels(freq_ticklabels)\n",
    "# plt.ylabel('Log-frequency (octaves)')\n",
    "plt.gca().invert_yaxis()\n",
    "\n",
    "plt.clim(0, 1)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}