Fig4.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:13<00:00,  4.22it/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": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = {\n",
    "    'Q': [24],\n",
    "    'k': [2,5],\n",
    "    'comp': ['log'],\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 4a: k=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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= 'log'\n",
    "features= amplitude_to_db(batch_features)\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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12a43fbe0>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAACjCAYAAACZvbxxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deXRU153nv/e9d2tRLZLQgiRUQkKAEFshASW8xWBPY3c8McSedk8Oxt057QTndEjmxLbsmcykT9zxGUeGnHTL3THTnsQx6Zx0Fmyc1ZkQjHEcJANGRuy7JCPQglSlqlJVveXOH/WqqJJKqFZVlXiffwRPr57uq/e99/62ex9hjEFDI5Nw2W6AxuxHE5lGxtFEppFxNJFpZBxNZBoZRxOZRsbRRKaRcTSRaWQcTWQaGUcTmUbG0USmkXE0kWlkHCGTF7/36ddsAFoBOAB0Amg7sPPJ3lVP/H7S8WNvbOzNZFs0sgfJVBWGKrAuAGYAFIAIwO1zFz7oHy/83cTjAOzH3tjY27Dh5CQBntm/VBNgHpNJkbUD2CaOu6nPNQwwBoBJHqfpElMsdSBEgPq3eb1F1Jkrdvn6atoQQ5gA7JrQ8pdMTpcOAHR85DqsVfU3/56uokqRdVF/V/I5qWege8P44LlyXl9uISCCrqAOnGCiCAquFcB2e/3gpFGu60KZJr4cJ5Mi6wx4nE20wEojjokcL11QZF0jgiNVsBGGQtFsWLmfD1Q7AAiMKQh4L0LxegGAUkP1HarAIke5JgBb7PWDdk1ouU0mvcu28dHBcWNRmaj+XwTg5nnxCwhOgVHHAbQhODqJhHDQmxbCaF0Jg2WFKI5fOen0/HqH13fEIsmDIXFGjnIaOUzGbDJCiKF0YfP6pQ9/6SHE6V2qRn9Mm8zQU/ZzAI6A2ANRHgBHDDDoloEQ0tl1oazlszb3pGu+2WvWRrgcIJMiWweggyX4B6byLu31g+0AtkGdZhXFg/HACUmSB/c0EMdzBULNUcQQpya07JMRkRFCCIAWxtihdF0zhk0mAnAHxN5mq/f0t/Vc+SOFdKUgcAWhj4gAdr3Za96erjZoJEembLKVCAoibajGvR3ALgRHuF0A7Kd6mi9XGR+unaNrEdzSOYwEDkPtOBTB0VAjy2TKuyxgjI2n+6Kq0GKNTJ2EkKYinZ3KzI8bgUMw8JWiSajtTHcbNBIn7SIjhNQBuJzu605DG4AtAMw80dMS/R2iV+rx/mnw4f2E7Dc+vdxXigl23o7jes1WmyEyMV1WMMb6M3DdKVGN+6iptECoWTHo3/9mveU/3zcw3tWNoNPgUH92PbPCb5vJNt7OpMXw3/jSczYArQG39w7P4PDp4jrbf//989/OiZHimRX+dp88um3Id5yWGVZCzxcCqlOw47j+tnIKslWYkLLIVIF1ATCPXvmEFs2fFw4f5ILQnlnh74DqAFzzHoaJVsJC5wFA547j+pasNi4GmUqdqQKLGYPMtNDSYZO1AjArskwJxwHRkfhcGCk6EUxB0YqCNRj2ncawfFosMSzJmlMwVeA4ntRZClUqrRHXBWbwOaVDZA4A1H1tEJaquaFjuRQ+CDsFAGiJYYnoCvR63rryX35Xc+ebkx7YVz4QMtqrVYFNEtJnbW47dLGFoFPIC5+1ud0KYXdXWIqXDBW7eInK4c82bDgZT5WKI3Rdz8AJgPDQW+dRQW8JP6ep6v9Sved0GP6dAETGGDg+fDlRPZ51VC8yyimw6mwr15U/cvms848nFCZHOQT/fKeUaYcgSkgyG6cu8aRlyP/Ov8yVP/5cnXKGlsknIPm74PUfg9/XRY3eE1udgY+fkmTnqmKX2VDfW0EFkQciRiN7/aDNXj/Ybq8f7FB/TryPTqj5YsLrUFDaAMk3Kjp7/nSZEOIoa1izUVHkLkxwkFThpUQ6RrLQSGFRrxeZ8M4JVKFFTQn/fOfWdlkJGM673qO15hbo+IIZmT4YUxwu6RSVFS8I4cETI6xCnVBXtPLTRE2WKIxhGRj+7PPALAqKleOJpIzxHukiGGQoUKC76oW3nIfOWEOJyN2NiNEaYGvAK1+yP/rebmLyfUO1udoAbFFk0Uw4nhJCREOhzW0otD3DGOu968vt/+oZ7LMool8wlddA0BnS9n2kLLLfP//t3o0vPWf3u9w7ANRCHWZzweifBgfP6egi63qcdx1AneUOCJw+bdO8GiIJTz03/Ge/84NzK8vXztl9Za7hPzUJ1BwudSoxcAoByEjgPCcqYwAAhTFYZVEZ94M5GcdzRA+BmCAxNzgQGAJ6iOMX4b3xvkycSpkO8616uoAHOBh0DRwnU+Ba6ROo79u06onf24+9sbF32cOHHvTeuPyGoFte7ZeMFzgqfuHE2+t6AUDQG1db5s4XGGMY678YqgFMy/eRlmDs75//di8h5GXG2OF0XG+G6ATQRAihC6334pzrXdRb7xZ5QlOe5lWBdQEwj0s36EjgbLNA+L99ea3nFwJPF/W7ZfnKmAS/DKqwgOgUu7lxhfDF+kUw8IvC1ylUlME+WfotD7IFEfV3ChhGC90wzakSAbj1cumgLI/MC4iXwOCHy/MOdEI1MODj6FynWWeuaG3YcLINsP5OGScWjs4RFB8aFZ/xdxH2XPj7AADGGAghaTF7bufVSm1Q69oIIVho/ZR41vlHn1caeTkN126VlYC5z/M+HZeHYDOtE5rKWsyUp1s5QlZVmXm+cc6I4hIPnh3w7dtTWdD0o8oCh2jgiyOvIRZw3M94kG8gov6OgYkgzD9cOHYMav6WgLwv8HPEAsNqmAwtMBs/BZmNwVhWD8Lrqfta14axwT/sGHcdtyhKIDSwTKzHC38fBSVV8A5flZAmsydtVRiEkDV5NpJBNfJvTmm+Kzu/+dHC+p99ynVWxxujvKyH/xi/1/nUkstHRgPnm+cV3A1CCIr1BGZKAADXxk/CL4/BTMvFUsOCXV/5QNgeOfJhQgxrx3F973S1cjdDH6wQIByggHF+uOf9DBbbcpEQsmu8t9rBmOyQAkOg+rmRze08s39pCxDtXQ6ePXKlbPHqp9PhXaZFZIQQPYAGxtjHKV8sy7Q7ji0YCVw7Zi++z0AIF/XA4xEaIWTBI/PferHO8uCjUKe4uQUcZObEJ54u2ExNMAjW0OmdX/lAaAEm23BIML9qrx+0QZBegD7wOCweDtXXOSZ4Rc/1Lp+pfMUyf399KyLq8VREALvO7F86ybAnhCwBcJ0xNhJvG6YiXSKrA+BijA2nfLEs8/Z9Urtf9m675O6iSwrvCB0WAex6+I/CBA81eiT8/pm//lXXjT3nnl7uExExMsnsnCJwXs5mbo78uAhg11c+ENLqyU5MHXmHz3333G++PH/VX/3a7R2yvQdGdABHACYC5JYrwQgh69JRE5iuKowSzHzlRaZw6PkCKnC6yGOTvCxVYF0AzIwxesF1sPm/LnjliX9Y+R97eIKlosL2do1IOOnc5yjWV3cvn9P0AIACRE+HaQ/zqOGKCOFuRNPn+ogc+H9HS5eu4jwD1UT0mBkt8MgcLz149Ce3TCkxQghJtLp5Imkr9Um1ITlEOA0VcSyWlxVMpzGZnnO+i8XWe4Qao87KAU8A4AKKs4mSIz5HSVFz65FV5yeOepiB7EIIa1X9V2UxYBjrP06LaoLxWAA8gK0x7itM8989doMp7N83vvRcPVIITWV0m4I8JRhcZtMGlx0A6AXX+1hUuB4llAcHgBBw/d7z8Mouap9zPwB8FcB2VVDZyuU6eKqjgtEERZbA8QIwTQxs40vP2UoW1naM9lwNfQ9NALZsfOm5hAsfbucQRkxU494+HPhkDyLKvGMY/Z297o+kqoIV4AgPPQcQAlxxd0PgdKi3NAO5k8PtBCDyVA9Z9IeOTRcDawVgpgUGIeAZB1JYgqiNZDF4+I9C7yby6C2Dy29dbv33O+d+YauJzikAQP0KlBu+i5yeL0CZoSZ0Wq7kcNsAbOGp3qyIfgqDKR6b0AGAGous8AyNQGcyAkl2Gm0km5ooG/Pt+yTb2/dJ7W/fJ3X88O5PXvUEzpeUGxetgJp4v+A+/7NxeWyswrgg1qLlrKLGuuyE474vjrs/hjo6TxMDCy605nkwWQ4dS6rTpBTCuOe1TTYArc6zQxsKF5fuB9B28Mm9uZ6zjAtCyGrG2BEgKDConqTCFNo9ekBqLl7vsQl4kydYOhwY+Ojr3U/+5pkVv/wIE4z7RIK4mYYQwgFYxRg7Ot25kcWozr5+WlhdmXQxatIiUwXWxRRm9vQ6qXl+UbgRs0FoE0QWXlh8cvR9NFrXYoFOBw5QZCZxh0cOSutKNowBsC99h+b0vSeSmQmV1Y9c6l1fXGd7F1nwLlsBmMUxP6UWPZB7FbGpQiL+HS74o5wBZYJOtTMYd3jkIBxz7hUwu+4dQLDwAcD2yA6XDKnYZA4ANODyg1r1oWO54k2lm5sFf4TAQIKe5DV/H+xFLeBIuOx8Nt57yqQisuAXrzBwQu5VxKaBSDsiXKEAAD4GhTGg0mCDkY/aFmG23PtEyPSnTE0qImsD4GYKk9T/54w3lW5CsTMAu/yyt9ul4A0ATsTe/irXmfHMTMre5fCxaztKVlXUQvWmZoPRD0Qb/rGOn3xAnJQmynWjH5j6vtL9majPp5pyzMc6sniY6r5S/cKzTTZEpgVjbz9Ssq+SIR0imy3VFxOZrfc146QjdznjPSMbhGywH67dt/7kA+K7yBMbLAb5ZfgDs9omC9shqsC6AJhPuI7SZdbmcHYjn4S28aXnbIOnzr9c1riwDgnUh+WCTXY7TCvhVd8Ruzjm1c7boVykzmJ+FBErxNXjGSUdIrsdpstwWim4HS6A/IvwtwIwE4KplsRlDM27nJrIztMJQBQVETzhQ8fyLcIf7igRzEhH0abL+GgD4B6XPaKRNwH5FeEPEc6/RhBvR8laWmm2E+48qnFvv+br/REl9COoRX/5ZPQjlAZkSCgNeM9rm2yNf+949p7XNnXc89qmdrXEKyFS8i7VRb2LGWPHk75IjhLLoyKENADoZ4y5stSslEjUuwzVDI5dHrVYaosiF9UkVDOYapysEMFE8e1CAQBPthuRLBEb48QbjriVsxB33Vyq02URgNEUr5GrxLJDOMaYHOP4bMUBgE6Y7BJ2FlIVmRnB4VNjdqK+tS/qWMJedaoi4xhjSorXyFVuB695OtJSM6h5lxpTohr39oDTH7XQOdGawVQNf623z3IOPrm3lxDSlu3c5WzldkiXzQiayDQyTqoi03r77YGWVtLIHOpbmLO6Cd5sNvzD9zb62JgNQOuv1/92/ehjY+8CaCv6qSWf8papYAEwlsoFUs1dzsqqWODmvakC6wJgPjbyEV1V3BTO3+WT0G5ujjO8vnBxybuIc/kiIWQ+AC9jbDDZv61Nl9OT91WxoUQ3gG28nl8OtSo2zoqKlFOHqYpsNk+XIWZNVSxu5iET6SiUMTaxDi0hkrLJVu9ptAFoXfSPNetX72l8F0DbkUdO5c3UkSDxblScy0R0lPCxGesoCY9kqsC6GGPbeD0XHnrV47OJ0CgdtdkK8rgqVvZL4HQzXz6ezHQZ3FrczyjREyAPbZQ4IQCgGvd2ALu8krcbav4un4x+qB0l4PSL6jZfiXSUlE2iZKbLYI2RyMAJeWujJIQqqO0Pkb/My30wDj65t/ee1zbZvf1jL89ZOTdcFTtTm+MkI7JOAE1MZpTcFFm+2Si3HWqie0cSIaeUszrJTJfBGiOJiYQnQH7aKBozSMIiU71Iu+SUfwzgKFQbZRZ7lxopklQI48gjp3oJIS8CGGCM3U4LSW5HUjb8UwnGht52Nlu5HQLNM0IqCfJZLbK/rPp0xehjY+2I2K4zy03KW1IRmQCEVyPPKkYfG7N9bcnTPwZgRLAzNQHY8umqhx7PbstmjtV7Gm1MYa2LXqjZsHpP436kkNVJKXc5i95xGYWs4AUOfOh94FB/mh+v2/pEFps1Y4SzOiLbJlj5ZUgxq6NVYUzg5AOiTVbIVkLAOQOj+GXf2yC8DL7QRWvmmR4OPH2pPfD0pbxLoRFCDAB8cZ7eCsAsuWUqmHkgxayOJrLJtAZkjgAE3c5u1FjmwVAxCmYcg1HPG6H26jwUWjmAgTjPdQCgklOGYA3nOpPO6mgim4zD7Re4ufpKLC9cCb7ADxAGiUkQ+ODrlZCfudoSAMNxntsJQGQKgxpwB1LI6mgvVcXN8moADr1ADX6JEzl5PuUFCYRKIASQFBkCSb1XZ5FE9vEIviIb074iO74/nMyHZhMR5dXbADgKjWIjAQSFQXT5KLzjOoUxQGYyeC5vd1lMiFBWR7whRa0cT9a71EayiKpRj+SBSTDRMotfdI7TU36J9/UMua9jMT4lK3IBT7hQbHDW52qPPHKqlzxKbvmK7Hi57UcyRFSNnhs7CwAQONASk+hb+g5tebLjkW8SghXDvpEfKYyFd1nU7azTcrVxktRItnlfoa3lH43Pbt5XGK5Neut+Z75+6dOWV+t21vWuIgv+AUCAMXZ9htuXTdISB014JNu8r9AGoEs/h4vaD149no/EW17NA8i7DfBW72m0rd7T2L7ohZofrt7T2J6NMvmky6+zsR98JkigvDrvRBaK3APYxulJQpF7QogR8Qdvb0ny5dcpbvGYS0xVXk1I1B6D+ZirDTs16v8T2fO1CsAn6WhEMiNZWrZ4zBMiUzF5N5JBHRAUkYEkvh6jCGnadDr58usE94PPUwwAxtV/56PIgm9SGRahKw37NXG/ICJdBRAJi0z1Iu3+G8ovEBGoy2PvEkB4agx398DTl2zH/uadnWf+7r2fB56+1L554QPzkH8iawPgltyyxJs4IEsDQtIbrhBCVgA4yxjzp7dJ2UE1dBcwxk6oye+uC6NXzHMMRbTYUCh+2N81vrysYXnhPzXmVWdavafRNto5tqPIYamFGm6KJ3Kfzs10Uon4uxDMbc0KkSFoEIe2SGoFYPZLfmrg9QAvUa5wjDfV3vhQfuXAzwC08V++Ny/Els7IfbKkEvF3AbCmqyE5QOQ7CYIGM1PAUwXcvH4Qs5cjVJoLNQwgv3IgX+OC8ZK2gtRURDaG2SUyC26KrBOAqICBL3QDHANHCCRZAvI8LpgNkhYZY0zC7EqwU8ZYQP13GwC3rMgSZxBBCGCfV4+jvefC5yKP44LTMSE+mDJagjwGavLb3u8Z+AUk3XXGoAh8sD+JwdEsb+KChBAOQKJvjSkGcCNdbdA2wZsC3c663s++9eTLgnV8LSFwAhBX1yzG4Z4zEvIrLlgB4FqCn6kC0J+uBmgjGYIxsV9s+j/PBp6+1DFxoYjqRdoB7OI5vnNg7Mae41cvrskX7xLJiUzPGEtL3hJIfWPivNxKKZJQTOzItY8tqytWRr04VP+dBeUxXqzKAXAwxg5lobkJk8wzSvdz1Uaym0nkuKpK1LfiSepbizXiQBNZRGVsBNN5j0cANGesReklbk9x875C2+Z9he13fafg9c37CtvTVSOYdAji+Q6r7dGXjc8+32ENV8e+1OLKFzslkk4ATYSQuDceZowxQoifEGJkjI1PdV42UQXSetfOgg2b9xXuxzTVy6FiVABmwUAogAYAWzbvK0w5L53USPZ8h9UGoMtUQqKqY9Xj+YZaVZLYi0PHfvra4E+e/fIb8t7dHfLe3e3y3t05c+8RgtkmGG8WK04zMqWyDfstSXa6DDZoFlTHRsXEIqpK+LoevPnFF5+VXznQIb9yoD0yjSTv3W0z6nRdCyrKN3t8/nAnyyGhtQIwM8ZiFStOhQMA9V5XYChJ717AyYosaMfMkupY3c663kf3fvFl3c66Ft3Ouu18XQ8AdFUWlkaN1BFCawVgblowX/jo4mUg9zqZAwAdH2AwlsUtmE4Aom9YgbE0LIu0BJ2TFVknANFcSjA2EA4m500UPA6C6xim9jgdACghBDVlJbh4bSB0Tq50sk4Aon+UQV8ct2DaJC9zC0aS9mLUZEXWBsBtKedE9xBLa4OySKQX5gBAK6xz0HMjvAIuUkSdUFc31ZSVYsDpgk8Uc6mThVZgxS2Yt+539n74rfG/tcznXkWai1GTEpnqRdoB7BJ97ESoQXnqXcaiE4A4r6gMfaPhl6NFiihqGZ1jcb3451PnfMiRThaqXvYNK1HbDNxKMIQQbqBTur73L1zb37rf2fLW/c7t6ap2TinirzZuLWPsw3Q0Jhuodlbre+eObfjUolX7cVMoXQDMfSMDVFIUsbakwg3AHkonqUZ+K9TtPne+9Zv/+9zrP/Ezxk5l4z5ikUh1KyGkCcCpdKaTwtdOg8gWAhhmjI2kp0kzhyqwLgDmwz2n6ZqaJeGUknpKKwDHT4/u73msecPXpstXPrJh2YPtz3xmc0WJpQlq7JBf9/WsjO6EEAHAcsbYsTjPX5epVFk6REYArMnH0Ux+5UA7gp4jPdxzGmtqlgDBKXAX/+V7w+sSCSFVAHSMsctTXuvQizYAXQeOXrLc21wXlQPNhtAIIQ0A+hljrjjOrQXgY4wlmkiPi5TTSnm+b6yDMUY7Lp/EXMuc0LFJXiJj7CqC5S+3IhjWaKgUPjzZF7pONsMalngEplKRKYEBt3nu8vxg3/ED549JK6vqYSsuDx2eyku8RAipu8XlHACo1WRAdXkhPj5/DcitsEZMCCGFCK7XyBgpl0+/eqLY9pV/M/2373YWLTaYyCEAbU8tG8kpL3OikT4eCLRZHnuy4vG1G3/8+hP/4xHcXMo/pavPGOsnhNwJ4NIUfya8O1BlqQWyoqD7wnVpef3cXAlrTEVjpsuWUrLJXj1RbAPQpcjM3HdGpjVLhbAdkitCUwXWBVVIg84xsftKr88vSWse+ubLZ0PeJSJeCjGVgf+tx/+q+b6VS7/a0rBwSfjcTVuD3qZqkyFCsBc/ueH9368f2PT9Xx4+kPEbncCtPEs1x9yqyMxx/qB0efF6+kwmw0+piixsOJ/uELF4jQCOJyKAXXVmoQ0THt4D8wdnXHjy3t3tALYpikKPXrwMo06HZTXVQeN+09bpNh2JvI4NQNefTp213NW4ONqwjxZa1D0Ld/xP/Re3NFT8y7fu+Ouo47U/yNxDDS5UrmeMdU/8Xai4AYD52mmZltVzIk+JGxmMc6Y6XYZrsepXCeg9JYMxUDGAv5i3jt+i05NQr24CsOWdK2X2WpMemPAgGkr70nJzsR7ytRHrnVcGhykBsLLWBoNOByRnK7UCMC+snCuc7ruKJdVVk3bIUb3IKOGKl87bTp4d/dXFnrGCBTWW8HchXf68Pd1CC41QW18ruLd0AXfw+Q7rSzGEE06eSwGApySRnX6SIlWRhe0QqieYv1wAAFHHmL/nvGyRRSYAwJxyjpZV8WYOeAHAJtycUpoAbDkzVG0PCc014pgkFGtxZy8ASJc/P3mkUB9U5HTFGKMnLg40uzy+JwR32R9aFi9ZgdRfVO8AQOcWFeLaiBPDrjGUWC3xiLV16eKighNnR2jvVQ9sVaaMPNTIEYrwoLoCsgTA557vsE4coRx+D6PXT8uoWjEzu3mnKrLQVtxRhnO5SVDoMhK+9vB1GWe6RAoZnxk1EUtFFS+UV0x608V2VWCRdk0TgC2uEYe9wLkMsX4XGhGuDY99/erQmEWWFUFhDEtrywWLSW+EwrkwCPfENiLxFFC4Q9nravDu8VO4a+likfL8dGJ1AKDLFhej+8wI5D4FtdVxiTNRbu5FNrkeLCzmq93yacKjuWa1EPnsM5p3TSmEoRr3dgRzY+EcGeXI+4h4g1zJXB4Ndio2r9VLK5t1giwzHD8WwPFjAXR3BeiHH/g3EELW/Po3IzsuXfKZRVGZVAelKKxVlBSzd1yiF3vGcPjjIXq4a8jyqz/0vEwIWT1ww3OHfWGFsHZpNVqW2WAJTssUnLI0VhtDdlQCROUr717aIB48cTqefGU4mb68oRg+v4xzl1yZeKjBzQkVBnZzlWVYzCRIy5l3xX+tbOTHcPP5ZLy4IeUQhiq0qGH/nStlMUe4Qsr/FsCWynkCrZwXPl0EsJ8xdnh0eG3t6KhMT532QZLCDgn1epUNBs91WC2UUoHD3FIjFtRYQu2vY4wdkQ+9+B6ARsSYFlVBpTQ18Zu29sp7d9uhTtcCz3Ve6B/Y9er+Xy396dyeaO80OsIf9V0sqS8Sz150er/43J9+/v3/SKVFk+gE0DRwTqGVDRzm6Aj0PGEyg2H7q+bFCAaTOzt2B7zPd1jD9xFqc856l7finStlk+wn1eiPcvOhemgNpX29rhFH2FuNuJQIYJc6Xcb8nVD7g+2xQgjIcFpHPvSi7WzP0HEd5QtqK4un/JuxbEla9zoPwMQYO5GOtoRssv5u2eJooQIHgBCCoT5ZdA2x8QWrhOXZCitlTGRTcWaoetIXPsHojymUGDZZ+HcTjP9bjSppRT70YjuAbR+e7KNL5peFpuhgeGTd16cdOQkhts2biht/+IOFn0EMRydR1n5Od8eDf2N4raqaaySEkCsnJJgKCUqreRHArqeWjWTEe5yOGRfZdCTrXWYD+dCLHWpbcPDYZaxaVBkSWie/7ust033eNeKwXe0PfHz9mmhuajJFxd6SERohZN33uov+6Ua/4rjRr6C6gYfBFK7F7Hxq2ci0bcoEObcrj/rlxuxxqqCy0hunIOxx3rOqFn8+3oOFthKxrMgUr1HfWlWpM+koEf70wRiaVhXQggJ+2vBGLFPkwdqhUQBzDv8u0LPALjQtbBZSDdmkjZwTWZ4RZdTfsaJG7DzRN/69n3f88H9VnI9n1HUAoKWlFCUlAj4+7sW4V6GNjcZ11uLYpsVljx+YEMq5ekXeOrea++71PuXf1jyo+xqA+5F6yCZt3NZVGKmi2ntR4RHHsurlZwcvFV244joO9c1z6s8udbqPJBzeIITAvtKElhaz+OFhz/l77jNsHnMpXROvoQa0zQDo1SsyTh4RqU4P69+/ZNryve4ibqqwUjZzyTlnk80GpMufbz9z0bkNDLShvjB0OOwJhw7cytHpl6+2njkhPuXxKILJzCEQYACD1O8JOMGjBADKq3gUlXK4cEJCWR1hwwobRQ4VJ4TQpsvM4GhYUEiv9Llx6KMBOJb632gAAAJBSURBVOxl4DgyKcpvLe7sdY04JsWsrMWdvf1D1Y6GZVQAgHGvAmMBBwBChY9IHlkRERHKCfgVMEEgCLCM5iCTRRNZZugE0DS/2kwry43489EBzC01igtrrZOM71s4OmGnQhUYAIiFlP+tR1Y2McaKxj2MXDwloaqOh1NUgBwtktRssswQTkHpdDzuWjNX1Ot4z4Nb39lPCDEDwanSNeJod404OtSfE+21qDQWY0xRGOMGAxJGh5SHOt4LXOi5IjGbXYDbyCAHrZ5cWvsZRrPJMsRUMT1CSNNDny6q/NHuhT/iCImyxVyK60EPc2+N+MxumbEvSQp73K8o3GhA4k51SZIswTtnIdk4zuO3iGHP5ZpNpoksC1y+0Py9U6fHn6ydrxfmzdMBAGQmiQPKoAIwDhGi6fH49/oDbMvlMzIVAwwLlwnQGYKFoZfc0qTC0FwTGKCJLCu4RhwdABy9fX4MDgZ3EnDLY+DNfjZ/gUA4Drh0XoJ7jEn9noCTGkhJzSIeAo3az67zgfmDWYngJ4pm+GeHTgBNtmo9tVUHdwUdkr1sZIwnp7pFMAWoWyhgwSJOqIzhTSJHba+p0Az/7DDpldQCqN9s4cTldh1WNOlgtgTf6qaWR8Xz+uqcRRNZFlDDFlFReY5w9yKGmHQ8942J5wKwZ2NRTrJoNlkOcasyqHxGE5lGxtGmS42Mo4lMI+NoItPIOJrINDKOJjKNjKOJTCPjaCLTyDiayDQyjiYyjYyjiUwj42gi08g4msg0Ms7/BwvapDEkHUTdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 396x198 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fig 4b: k=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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= 'log'\n",
    "features= amplitude_to_db(batch_features)\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x179c9fa58>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAACjCAYAAACZvbxxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19d3hT573/53t0joYteU/A2GDALGOWxZ5pIKPN6MhtmoYmTVty+0va+5TW6bpdt71NadJF2lvadGTc3GaUNOlIk2ZAIIDNNMNsM2ww4L20jqTv74/zSpFleUg6sk3iz/P4STh6dd5X0ud83+9+iZkxilEkEtJwL2AU732MkmwUCccoyUaRcIySbBQJxyjJRpFwjJJsFAnHKMlGkXCMkmwUCccoyUaRcIySbBQJxyjJRpFwyMO9gJGMjlZ7AYAKAHYAVQA2pqRX1Q3vqq490GiAPDIEwaoBWAEoAFQAXQDKRokWHUYlWd+owLsEg/ivVVx/MJ4bn2ga10tClmTVv2eJOyrJwkBEBgBT/v7Xqc9bLDQDACSJgq87Hf4jN37w+D19vN0HoF38dTCzN3yAIFhECfleJdr7WpIRkQRgIoD0kMs+ACeXLrG9BWAK3pVkgEYIf3tL+a8QQUcjIgVACoA0AOOJqNf3+8snM76Sk2ew+XyQJ0ySkZom6SYhRyquWUnmPXdvry1HLvpDn5KAiAjAeAA5IZcZwBlmbg0f/+LDv7M7O63bLp8sMY2dXkOL73xGTc29KgPwIg4d7UTTuEqxZhw+6EFOrgG5+QYAqCrJql8w2PtcS7gmSSYIFnHLCRCNiPIAFEAjUgB1zHxloPt/udQt7s9WgBQyeNlocvvvf2KdPzX3arhk25ySXjVoCXSiadwfANwT+PfpEypkmVBULP+xJKv+3sHe51rCtbpd9lDKW9rcSu35Ttulq45HPvI5+jEAfOX+Uun7X5l7tyRRqKQbkGCh93d6W5Uu70X4WSVy+KSnH11kWHTnM2BmpKbIKC42KURkj+UDtDT5cLHeh8klCpoafdi5zZVe8hHg1fPZvST02sLGa1pXu1ZJZm9udSknz3ZAkSWkpxoxb1amPJ+yiph5bwRJNwfAXd5z95aFbqliC80GMBYhjul/m/D6SlmyKBe7dyLNOBETbTeBiMjUZuW5c14iAGhr9+LAAYf3UoPn3F13Uzk0SXqCmf19LZqIkp54MavckkRIz5AwY5aCU8e98LgZrc3+aVPKlE9vejntkfB1v3o+u+xaJto1RzIiUp7ZtOLCjJL0OYvm5oRvXVXi/4OSrr3Dg/MXuxSPx99D0oXgKoBDzOwLXPhyqXtrg2PP1FkZ98kt7hPwswoDKWrelBM+AAYASlqqrM6dK3fNRfKXmbmOiGwA5griMgAXgGPilrMAGAE4yhcb35IkChoUJdMVAFBL5xhf8yT7rnvufxwppQsUQ2auhJyxhveEUXBNkYyISgEkfWDZmA3pqabr0Fsn2wgAx0+3LWvvVBUASE81omRiKkwmgwygiJn3DjTPjivfenJ62ifvVqSkpBxzmXLZudebn2Tvsti6bgBwNyJEAJi5E8BesU4CMBvAZwEQgEMAOgEclyTaCOCuCGt/qnyZ8eVZCxXDoUoVPh9wbL8Ki5WUoilyTFvySME1ofgT0TgAhQCOMHO72A6/B+BGMeSV3z5z4tHPf2NXMgCuefP2z02ekLoOvd0Pm+WiP/QrEYTbYeGGma7zELrRyfYt5yfabtzws6Np/W5ZRDQWwDjxzxPM3BbymgnANADKqrXmnH/7VPInFy41FZstVAngKQD/bHarqZ1ev8TMOLrHi8IpBkgGqPu3q3/5wec7fwxgf6jEvVYwoklGRFZoEqGemc8BPS1Lj8enVNe0eFUfO554/tSHfvenk2+Hj0Ef1mc/cy4D8E6obiXWUcTMRyKMT4VGHgZwkZnrB/nZFADTASi/fDLjK3ljDB8unCTLLVARmPjkIdWflil15Iw1rFXdvO7YAe/Ktmb/MaMJ3/z2fR0nBjPPSMCIJJnYbuZD80kd5JBFemrv2bT/cPP9VdWNssloQOnUdKSlGtXsDPPvssue+ffAuGj9aGLeEgDt7S3lSvh7UzP2jGHmSjHOCE3PMgDoAHCc4/giTzSNq/R62X7mpBeqj5E/heD2MxSiK7/7RdcXltxg/E1BsZwEQGFmtfaYz/mHjd13792qVjHz5VjnHSqMOJIR0WRoDtO9zOwOuV4I4KbvbZjz0MrF+YUL52SDiMDM6OhUcfp8x+GFt/ztXmg6EKBJFgq5tRPARQDtkQhBRCkAJrW3lDcighT8wX/Xr9n4SEMOgGYAHoQZC/HgRNO4TQDWA1COVnswYZKMpGRJBbD5XLcbV+p967s7WJk4PahCqwA231DU9BMAuQCamPm0HmtJBEaM4k9EWQCmAjjFzKfENRuAj0L7wc8DePxrD5RNh/hBxBikphjVeaVZ25h5Xz/3t0BzVUwionACAsAKAG+cPuP69sQJJqskkQIA5867lcarqq20NKkCwK8D0kxnBI2B6bMUZd9uj3f+IlPAkHkhd5xBaWv24+heFdPnyYGt1s7MZwGcJaIsIloIzaKtZmb+9dH0XpL8/hmtw+IGGXJJ5tv9g9APX+Nyq/T2wXPlZqN8YPmcCV+TF33zEoBbAYyBZpFtEZYbgPj0rb5ARPMAnATg2r+ntKqj0zebGWA/o6DAhBzNU1KVmrHnXmj6YUesn78vhGZmnKhRj77wdPfjT/6ma+er57ODUs7tYhzbr2LGfEVVjLR5bWFjDyNGPEizZi6Xcz730+QnFBP1+o6Gg2hDSjLf7h/YAWwDYAJAh05fhtOtonzaWLxz6IJ337GLnl/9ufIHtRdbnuxPgY5F3+oLRJQPwMbMJwGgo9Ue/FFDhgXDR0S0IFyaiTBUj/U8ctgU149JREsB7PznuayxCHmomFk9tFt1XTjtW/LYN7sOR3rvLw+lbbp4wne/V4U8saznFnv/jNYh97cNGcmEBDsFwFR3pQ3nL7dhQn4G/lV1Cl1ODxbPKsTckjEqgM2Ghd8Yki9CpPUsZubtgWsDJSsSUTlCXAnvxjl7j4+HaMK4mM3MVZFCTTcUNWVC0y/Phr/310fTKwHYO5r98DgZWeMMgZeq7p/ROuRB+KHUySoAGD2qFy+9fQz26QVo73ZhZnEe5k8bGxijQGQoDBEWAdgVeiElvaquo9Vehr7Trg9Cc6sE9L8KANZm1zEl0zwN0Cm5kZk9RNRFRBnMXBd+L2bUEVExEc1i5kNhb68CMCclU1LOHvIiS/PchUZEhhRDSTI7M9MvX6jEfbfMQ0qyGcyMvccuho4Zsi9CWLFnIyUWCkJFJAgzq2F5YvYrzoOK09eETEwLXNPlYWHmGrFt7ujj9TNElE1ESwDsDLGag4aEJEHxeVk1yBSMiAw1hrJaqaqqpt47t2QMUpLNAICGpk6MzU4JvN4jNJRIEFEygDRmvjjg4Ai4s/g33V+cufWpXyz2VjKdtiXJqapZSkOI6qHnw3KciKb19SIzN0KTqqtEVAFCuS8DsHlciaHq6HbvFgyT0g8MoU628cEbZ8+YmPP2jIm55oLcVAUA76mp986bOvZFSaIiiG3JsPAbun8Rtxd09dBpdjTevK3J/fafY3Gg/mKxtwBA9bnOSluOpURucp7xjrPONZxua/B6/C4lxVioi04WCuGe2MfMaj9jCMASaI7hprDX7Mw8LFslMIQkI6JlNyyakv63n6y7HuLHvr3i6X++/PaxvydyXkGwoGLe5jnoTZILu4xS+qwX66xRk+AXi72bAKxvctUqtZ3vwJ59N5hZdfpw7EDTVkNB8vK3oIN1GQqRJr6Imd8ZxNgyAJ3MXBtybSKALma+qteaosGQkIyIZgO4DCCbmQ+HXJ/XnwNVD9xe0BV0Sbh8V6H6W2BTpqoANr9YZ41KMb+9oKtgYZ5xj0WRck93bEeSIQ1GQzKyzBMBoOqLu5T/N5gsj1hARMUAnMx8aRBjJ0Jzy1SHXOvlehkqJFwnEwFkBlAEIDzAHO51TwTszKx0qDVw+upgU6YCMSjmAYnY7PJn13UdQL5lBsYkl4LZj8uOYwk3WJj5DLSil8GMrQVwkYiWiW0UQHBLHXIMheI/C8BhAAjVgQT5dPech4KIpAbnP862eCq9FsM4pBvnBV6KhRQVAKzHWs9IipQCi6wVOGWZi/2KZHH88eQnXtZv5X2iWuwKA0LoZVUAVlU8Y5u0cVvKvV/9k7X610fTN4mQ05AhoSQjounQskNnQfMvhWIigIQEdYnIJJTl8mS58OuZpoWdipQSUJpjsmJdvquLmt27FbfPh9q2sbjU5UOH24/L3b7GTHNR6YHm5w8DKBNO1IRAhNckYR0PZrz753tTT/m82A/Cp/1+lEJTHaqHkmgJ85MRkRlAMjM3EdEkZnaFDTH0lw8f45wpAGYAcAOoCtz/9oKuXs7VwSr94p4zl2a/cjzLtHQWAMXtA061eQGNsM//8JC57gtgENGrRqX4i2XFjeND56o+k62nxXwAmhUZ0XcWDqOZKibPk80XT/oURwfD42TFaKEhTelOmOIvkv92QMtodTNzQ9jr8/VSkokoF8AEaFbVUZ3umQTNs9/JzIfDrVSEhI8ChC0rbizodL59lGCwJJsXy0QUHKMn0e76QtL8BdcZHyyZrUzFABVNgRATAJzZ74VBAYpKZWAIQ0wJkWTCEqplZiaiPGbeHeu9yoobe8XtAj8YERUByAdwJZY5It37UG3OVQDzADgA7ArokS/WWesGIRErDGSzmE3T5attP0W67eOKUR7Tr9SINrj+6vnsgnUbkl87UqXa/H6WJYkGqmiqglb1pOQVS6h+Q0VRqTykISbdSSZCLrnMvFNsNb2Ue+GZ9gx0L0GCXqVtGbaPf6S161kXgHPMvKu/ewz23szeOQ73vnUpSWs/0eF49R+RHLWCUP1tMXYQyV7vFaTbPo4ux1YYlSLFINmWA6sB9HQMGyXUjLPJtwupGfx8Xy519+fIrQBgLZktyzV7vZhpVwaKlQZDTMmpkmK2krexzufILjAMWYgpEYr/AgABf8x0Zq6JMGYKtPytgRAsbWNmON2HlW7XHltextc/zsy7wrfgKFEBwMrsVxyu/XC6DylJpvmWorynbogjlboKzF7VdwVGeQzSbXcC8HtBhn1EtDDfctPNfvZWQ1O+7TajtI7Zn9roOqK0eWqBnsH1vmAHoChGQmGJATX71MD7IrpkQkNMAKqmLpAf37S++0dDGWLSRZL5XnqqAEDF6YYrS3du/PYh+5TiC0R0EUBfir0pgiEQCXaIvK5u1w4kmcohSWYZmq4UL+wAFKe7GhZTKbRk07gD2xv97PyUJCVbAMhEpFotS7tqG+74K4CT5ZlP/Feb54AN8MsSGSEbWHKxjEzTNLS4jwfuMdAagtufLVVC4RTgcKXqnVEu97n9CUIFpdyGC/QAERmZecDdRA/ELckEwar9zOsbWtpm26cU3wmg+raF89ait9siWlRBU7AhSSmQJDOgX/A5eG9BsLjvXX0mu+5q288+5fVdfUHcZzOAsk7HWy8CKHb7mhZlmMrlDNMCpCilSFXmI9c8F4qUFHqbgdawEZoxoQJAsk1SJ5fKnQ/d2b49Cmfr7wB8OtrPFyv02C4rAFj3nT6r2KdMBITI/8SKxff286QMdjvq8YVC30wNce9gqk/c9yYiW0lyXuf1qUl0X246An8AwMx7Or3HTzi8F7SHhmS0uf3wMxDiyhlwDUK5D25/ADabk6jscKX3JQCrRZyzXzCzE5q/zRTrZ40GepDMDkAhAkyKJhHONFxRFkwpLo40eDBfQgDCiiwDsNnvdxyBkAx6uAMC91Z9l7cgROrEc+9paWtX3z3piedsxryPQPte1gOoFpkbyDVfX+HlLkeXt1YFAB9Dre/0trt9jqd9rB4OrGGg4Prawsa6tYWND64tbFwg/lsnKru2QSPaYNSgfwFYE+tnjQZ66GRBHSFwoamj07+gZFJ4tmYAEwHU9vFaL4gf/UGiW3UPph+qzakH8Agz74n3XkRk/EzJnz9BJFmZvTIzB6qKgpaf5gaZXtri3vN9Veoqv2HsPHVskuRvdF/y/PPiy5995PDyuALYzOwlojeh5ZZt70/nemi3zbXtV+7/+OK/bN+2pNAuABsfXtCREGNAD0m2EZpPCQDQ2tWNNGuyBOB2oa+FI5UjNJ3rD+LH6pXBqgPGAdDri50zPf2mIgDK2OTZONu5M3C9hyL/Yp21bsPMZd+8rTB3jJePTDdINNvH3k99rOhrr7y82ht3qEdsvW8CWCKql3rhq5UpBQCqJy2Xlzef9c2DkLjiuu6Im2SGW++uA/AiWLMkT166jJKx+QCQhMimeCyZAOkAoiLmIJEHYLA9y/qEULglA8lVAFRFMiPLPAkNjqNAZEW+IsdcmGRT0uVzXYfh8nUpFoN1INfFoCFcMFsBlIva1XBUALBmTZBkZzvD0eofjOskZujljJ0O0ghreLeJr55FIekABsyjigXxtBcIwRxoMcV6CMdnijFXcfpa1Vb3he500/hwRd7u8nUrl51nMTNtOc507YdEBl2LaMTnepuIltzz38ltC2813o93owpLASimZII1S0JLHSMpPXFFPHqRrIoZcwHI9K6g0jN0YYVmdY1UBPx+db9Y7A2GnnItU6s2H7vltZq2V1q/EGJQN7svHrzsPDu3NH1F6PefkFDP/xxJu3DmgPdIZ6vfYkuXFABzrDJ8XV6oABSSAFs2oeOyX03JkxISatKFZE6PZyMzrwNgQYR+YQGIPq4xbU86SRzdIYo8gp7UL+yUezg+v4B/gIhWE9FbIpabPTfjhr9/u+yvH0MgXMbwInFFNBXFc2TL+aNehX1ASpakpBkldHv9ftZ+JyU1X1LrDnhdKXlSQkJNuoSVbHd8pu23r711C4DNDrcn6GoQ+looxkLbUoYdOmaJpg3CkNkBYCkRjQFQuK/5lZcR4utq8TT8GUDZLW/KibDu7ACUwhkyGuv8UD0MWSJlTJJ0DMBm1clHAGyur/at/NHCzoT4zXRJ9RFN6rzQpNQcZt7fx7iY0nsSUQsg2hMY+muHMIh7FAHwDDLvfi60WO7TEV5LWK3Dr4+mB2scmBmn9ngxxa4EWxaEzk1ESwZTrBIt9NLJ0qBluY503SkU4yDab0aDklU1wSyK/Ok/uWDLWfulgd4jCjucACqJqJCZz0e92tgRzMIgImVsiUE9f9TrLJwhB7bGoES/6T/NbXc+lvS/heXyJIi0Iz18Z3plYQQU34ipPSEYUXpVtHqeIFg1gPVeT7PdZJ16K4BqcT0iRGM9AzMfY60lVraobxgShGdhJKdKm9/+k+eWf5/Z1qO32lcrUwpKb1a2G4y4g5mD0Qo9fGd6SbLA05AKoFcDECDYa6xbp/mGC8HUI9VZB0vq7H5zuYhoJrR6x3OBa8y8l4hWEdE2vdPP+0J4Fsb9W7SuQUTUGDKsAoA1t8QgXz3lR+4U/Tpv651PZoHWiC0I30tPFfheemrTth9+c6fjhd8/0EcU4FqBHYDicdbBYMwIXIvoXyKiOQDaQgkWgu0AliVqkYPETgCLQ/6t5amZCf53ZZwuvjO9ScahW1AgDQjAepMizzTK8ucAVF/DRKvyOC+oYC+MlvGBa738W6SdUnKpL6NCNHk5LIg4LBBS9AiAUnEpmPoUAl18d3qTLNwt0N+ZkcMGIsoG0DjgwDBcOfm9Z+H3OoxJE4KpR9mQHU8oE601a9XKmrXqpqVZa24DcJoHOMOJmVsAtBPRhOg/gT5g5mYAqvBfboS2teueVpXo1lHBzNYQDHUPskgoBBDRzdIXiGgKgCu5U75VCmFdFpKx5ify+NsloruYWdnT+vbcR2c9c5dNSS0D0DLQPX++SFWPtP79xwTf/D+t8HabJfglwg4AGxPkM4uEYwAm/2hh566HdtvK3J34OhGWQEsb0sW6TDTJgmlACQo3xQqKRukWFqKP3+0w/SAA1KxVNwFIYmalsmUrZqctlM0Gy6CU5UB3oJlpa6zVTf+jGOGFRDKYMYMId7282pso52w4CJp+tuThBR1vE9F3Acjx+A/Dkeg2BXFntiYozWfQrgsimgpA5cgtzO0AlDa1GfPSl8JssACDl9QVAKwKOpTyzJtxrOMdMd/QqxSstSY9JqzhVGgnD+uGhJJMhJXKMHC4qT+kAWgbcFQCIOKS7tA2TGGoAqCmG7OgSFHXCdgBKBLcAPxQpB7dDYZcpRDN9BhaUqmuria9tss+JYMg1IOrbos5dJIOre3UkEL08XD04YIIoK/DuAYjqasAzPGyrBxv24Hrx/Q4T3VYVApmPkpEawG8pud99XbGJgI2aF2z40bJqpoC1dXwn4XlL9lLVtVsB7DxxFvTe0lVIpoBzdLqN/wz/VWlrmat2quqfPqrymAk9UYAdx1qq7TNSFkgMwNEADNUoqFpa9oH3imcb7jtq5UpKxHymeIxAPQKkM8Xnuw+pVWsQWC9gseBkJCr86jNZJ0qExlUAF2yreMGJa0jeLzg+e0/fKnt3NYTrHWcTihmpt/8wY9M+NmtTm/T0rnp5Z7hsC7DkxYmL5dvXvuQ+WlXBydnTTT06PkRK9FGzLE3QwDNZ8csa+37oQBsZb+0DeKgVGdL7dz8ufd9snDZ12YlejFiO95/9MA//g5gbqI7TvaDHlJm1oeM66xZUrKzzac42vxISpPiDi/FrfiLflxDUokcJ+wAlJ5nbpHi9xhNABRnyxnI5hTZmJyTjARbdkQ0HoCXmS+NgGTMHqqOyYoZAJTsSQa01geXFpchood1qbvJmyBUqe4rqmzMCrnELJncBAA+TxeUpCwgwZYdEWUAyGJxzE7gcqLmixbuLhyFcDkZ3t3n4jJE9CDZsLkYosRG1VnnUixj3/XZEXtkW5cKAOb0iXC2nNGuJ8iyExXbM/tK6hwJOPRXz5PQuWpfD5LZAHQmsDZSF5zcOoO7m7ddj5DyftnatVySfV0AVNlkAySDV3W2diMBlp1I914KLQMjHMO9ZQZxervvCoRvM5CajTiUfkAH6zJg/RFRJoAUjnCgVOi4wd43kIHqaNu7Milt/lb04W6IYp2LmXln+PXZ617r0YTuzL8e2tp1+eAWvXUl0o6v2cMhB8WGvKZb18loET536L+FcXKOmR193mAQ0NO6TIFOullIBqpVMlgUACUA7ipZVVMWC9GEoh3xfQefXNMjoY9orQnAQoQd7BUPRMfqY5EINgIR+nDVASgAENd553pslwGlNQXaIah6QO8UobGD9XsJIlwSRSJxg4gmQUtebO5vWj3m0glBI4S1btvWeG+oZ+xS5n7O/okSwt3QI1EiJquPtLPLL0TzHuHpz4u3tZLI0zIPEJoCRpB1mQgM5Slx0aAKgOpx1MKYFMzpi9XqG8uxnQZXCW3bjAlEZAVQxMzhp7C87zBSSaZlafqcXsmQBMRoRous03OxLEAo/kdF+ktUIC2kUB5FR+7h3C79Yr2B3nG6r0UPkum+KKHcl3ndV3s0qItB6c8bTOFtX2Dt6BhJOFCjwTIAb0cxfji3SxWAfNsbqQVrn7P+auXm5N/f9kbqptveSNWtDmPExi5Pbp1xBcDDzHwglveLgtqI7pRowMyHiGilKGGL+EDVrFWDbpDXr/zl3M/Knv32Fw/e4Ys0dgRCXfBfliIAu9gLqymdFADTANx12xupZXpMEDPJfI9tKwBQ8dYXf77S99i2rbeULn1FjwWFYDqAeE4XyYvkF4sRgVM9enVCFASrBmA93XVUmZW2YE6OKf/6mrVq2SBTfoDh3S69qZMNXwFgVZ2sJOdJQE9r/o/xThDTdikIVg1gfZLRPBPA579xw90v+x7b9gfxmh5QYm0BLk5EOaPTOgKNfBuFvy0cFQCsdY5aRZFMyDHlx+JuGe7tcj4AxecCQhJ0A9Z8R7wV77HqZEE/VlN3O7rdTglausw6ANU6Ei1W5A5UkhYtRAp2AYX0YweAc90nl1a3VSoGklGYNClweSRUZA0WXm83VyPYbr5XwU89tL4hMSNWkgVL3dZMLcfpxovYdfZogGxxF0GQdtReTKEMIpoMnTJpI2AXgEWBeYhooZ/9+8vSFqhjLD2E3KDdLS+v9hZ8vfTPX3l5tbfy5dXeTXr0jY0S6oVX1c3QrPde7eaFFI/Ye3awiJVkIQctEMrGTYK9cBoOX6rFrrNHlQ5Xd8z+JYFp0OoBY0GWKIrQHaKMzkBE9wJoZebdE61Tv4MYsxYEoaozjGM+jJCW7ENMNO/p5zxNAMpczf4e1vxfrmvXJTs3VsU/UECRCkHU8sKp2H7mEJYWz1Jfrak8/SG6eTGAamaOpfIlqrrIQJDb2Vq7fNrtT1bNXvdavYhJ6gKRQTEbgBnAPmipTSoQd56/UDs48Dvo1uQkCqgA5Jc+0FEL4BFmR9zt5sMRE8kMD6yo8z22rQzA9wB8EpovSZqWV6gevVTrvHHGwgpoe/kssfVVA9qPpXd2gyBYNQCr3+tSLOkTpwH4yOx1r5XFSzThnJwL7Xs6FJKNcCBt/NSPrtjw+AoAdszUiLXt0c9EO5/9QneNYpFTQq8NtT7nFXNmA7iaiAniTvUJuDIgnuKPPf6tv79Y/fZOZu4AglJgFjRd5tmBWl8Kx2dGH8W0vTB73WubAKx3tV9QZHMaZFMKoD2dmw8+uSYmaSCU+7nQXAsHw63cFRseL/B5XNWOlgabLW+CjJBii2iI9vXSPz87yTbv9izzuFBjQgWw+ZY35YRLstveSC3wefih9lO+FZ4OPpoxU674x20dF8LHNOxUH8lfrBRBSOlot9G4nbGGB1b0SJXZ8sA2ENFyCI83M7PvsW0t204dXOrxqQ/+36e/U5NiTvrGzb+q6OsowmJE1wHRDkDxqw5QUnbgWqzBdDO0tgoqgL3csyAgFBUGo9lqMJpl1dkJxWKLapsTEnLFl6Y98d0s87jrEVvdZlwQHv1qSYGVZFJMGTTVaKM1t72RGtTFAmPMmZINGlfmQDhpoyFaomKXh4ioDHjXpzYuPfuO+eOnTv/onJW35tjS96ydZr+NiJI9Gx8nZBwAABGDSURBVM4WeDac3eTZcLbSs+HsppsnXpcT5ZZaBUBNypoKR9MJsNZcK6pgOhHZiGgRNINjNzP3RzBAEDspIx/OtqbAtUERW5wSsgrAjkdr1tUg7DAuJK5BcTgqIFp8in/L6O0ZEBVeiKQzDhoJCSsxcxsRjSOiTO+mrRUArAaSFK/PB4lImVswBX/7/I+uu/OxRxzb6ytfmJMz02w1JisA5jy04PPrPBvOzjQ+OmGwX3SwituaW6p0Xan2WvPmdBFRUBq03dHZ6wjmtOdsdSKbdwq0c8ajSVIMOU+Kwdo5SgMSW7SsKgHwZuBBEoQaKiU/FNp5n01+mNII7jYGej8ounRlSljskpmPENEKP7NdIlIud7Rg4YQZgZcViST7Mx98DMxsPtx0XHF7PWjovqoAoL+e/tcjH6X1P4Z2KGs7tCNv2iNJl4NPrqmbve61MgAVJBnslozJ+2vf+PprnZf21QFBgvU4Yvqys+GTH83/8D0AjkdJrgDePXI5a6ziaL7oTc4a1+82J5Ig05l5RwzzJQJVAOY4r/qVjOky3G0Rd4Beh7NFGDMgEh0g3/u3w++cu3HGwjmywRBpoXYiUmZlT9MGX67G/LwyGVoe1l6RgpICrR9GYSAlpQ88If4AIJOIPgtg5841uzcUJhdZLQaLUu+sR4PzkpJjyk1+YdmWD6Q9Z3splg+17dHP1K3Y8HgZgAqDYrI7W6+eS84a9+W+lH6RLuSONdifIGz0e/kuksimdrGsJEXUBzc6r/o/ac6gZMShM+rSpqA/3LPwxpW3zFr615tnLjLJUs+yd9/Z8RUQPeaBIMlUAJuNj06IawsRccaMQzcd+X23zzHH6XUgx5yLsUljA0Oq0p6zLYhnjpC5pgG4GLCow15bAOACx3deekKQv1i5ef43LTe1n/YtT51s2GYw0o/CFfq8hcoHF/4gaS1CVI1orcuEk4yIqCA95/tnvvunFImk4EIND6yo82w4G9zKmFnZc7naa8+f3QmgLAqdrL+5p26//p2HZqaV3oXeIn9z2nM2XXQh4aaxM3NlyDUJwEpoVmp/beeHDUS0kJl391ctRUR2Zo6rDnUoSDYTQAOAyZEyRQXRKg43Hl+WZkqpKkgZ8196ECyAm8d+8KZfzN/0v5mmrHCRX5b2nE3PqMCCAMmEA3ohgG2iCfGIAxEVQPMw1fdFMiJKATCGmY/3vkMUcw0ByRYwc6VIv3H31SaSiMpZhxN0I2FV7uoP/3HRkx9INabOQ4h1qeccRJQOIAdayGkSgJ0joM9FnwhIMfH/fZHMDq1WNK7PkRCSBaTTyZbapX7275+aOek7xkcn1BHRMmhfftBKDMQduy5Xr7Lmlb0FYKOecUcguJ2tBvBWNDHRGOb5OLT6yupEzaEHhG+sLKSIt0+SxbtVAglwxoboWetbXe2zp2ZOuhtAtbi+CyEVQIG4I/t96yUlaQZEFoK4rhvEk7gNwAo97xsKIioF4AbQVyRjJGE2gINA0DnsCh8gujXpUuKYCI9/BQCrn/2Ky+cGQrzEQj+pp3d711cAsLra6xSzVreQsKa8Yu69wrOvK8Q9rwL4C7Q47UiHIURXnATt8LVwlAI4rMdkiSCZ1nCXJBSlFuDAlSNAiJdYFM7mC5GtFfH63AAFl5KwLARREV0bCHnFCyIyENFqAIeZ+cpI1sECiFDsHDh8LRyhRIwLiSBZMKGxMGUsxqeMxfb6Kq+f/aF7+25o22YVANWSOQWOxprAawltyivSsjtE+4CYQdqBZCugWZChxy+e06vFQYKQP1CZoNBhdXtgEkGyHr37My3p6vy8WZ3XP3fnO4FgrFC+z9Tt/vlzALqISLVkTIaj+WQij0kOQnQeMpF2km7UEO0HZjLzm+GhLkHiXB2WqTui0LPiyUzuPW8irUuEeIlNP5l4EcByAAeYuR3QrJdJN/y0JTl7+hcB2Lsba060X9j5y6tHn+tVepYICBP9JDMPuomfkIAWZu5TXyGiedCSHPXqDaILxOfdF/pghFqW87ZMKwBQ0Xmke5VtZvJbADbu+/Cx+J3iQ61GiA9az8yXiIhWj1/y4Vc++tRKCEJ+aMs9r7927u3XRAHDUKxnJYBdg2nrJFpAtXMfPdhCxhmgNRtOiN8vVoT6xuZtmVbg9/i/6qh1LbdOTdoK4CkA/wRg7T7lVJInW4JO63iJNuQkA4J98j3uL9W6rjqaDjU7W5OnZU5WAKh+5q61z3/i/m11u58borUQgOsAvNGX4i7GLAFwYrBFKqERgJEAYdE7mPmKkFjVrotuq5IhKwaLQYWW8SJ5GlWFZIKSLgMi/Lbvw8fiCr8NS5sC1k6+KPrn2a0bb5iwMrmh66pyqesKxlhzFYnI+n8f+uVqIjo7FJKAmZmI3k4tKPlIMGcf7+bsE5EMLQa5k6PrOHiJiGLtKJQI5IaE9SoAWH1Ov2K2aO3mfU4fOy94iAxA8qRgBZwulv6w9cJg5nPn1u+esvPiXmXRmHk411GPXZf2oTRrqpJlyZgDoImIJgy0NemB5V/6ba7q6v5tx6Vaa8qYibLEVC6z9Nn7b/zdjWLIm9FGCpi5TmRgDDvJRFp5qDpg97n8CnsZ7sseeDt8kMwSJU8xM4VU92KEHqoaFcZYc3fOzpmhbq3bhQJbPhaNmYfa9gveF0/98zw0X04WaX2+Eo0KxZScbErJkLubLsJPTI6uFpNL8rxydI3ndByhKB4gB26oMAuiYgwAXBfdR5pfb/ORRJAsEpKnWGAZb1KJyA0du14HMNz9yTYmKZau5eMWqDvq96DT06XOyp7WOTdnZgWAcgBOAOVhT1ciYAdBMVnTQCShrf4k3J4u+MfnGRFf9KEaIyMCQMzsJ6I8IvromR9ePJ59Y3pH8hSLqqQGda8uaH6/HvUGeliXw9o6yvjohDrPhrNlBslQsXL8IvuWk6+cP9V69uH/3PHjc/w4nyOt0UcugI8BSKQhUCUxlXe3XSGv2wFPVyvSc4tR1JFFiEMnYWa38E0NG4TLxUpEC6GF7I66LrlfmLdl2p8Q5mYShNLdET4s1mV/EIFmBzOfCbm2CMBkAH9j7exuPeczZBaX3TBh3k1bjGnZRtlqg+QDuuvO4Cfuh9QMNXnz9FeVmK0rkeLUkajWCf3MmwJgBoAFAH4DrdSwayh03F5rGWkkA4IN7KzMfCjk2lwARmhb/N7wgtsY5kiF5tn2Azj4/bt3zXZKnm0XLM2mImcWLWooUDvaG5ylqfNnRtFnrK+5hsydQVrDmUxoncjPQns4jdB8k8OSAj4iSQYARJQLYGJoNRERrYKWsjMfWnn9gWiD0sJflAtNutSEvhbaMRFA1fdqHnj22frfXGDmqLpnR5hzAYCqRAXQxZY8B9oDeDogNcW8JmgB/H4r9xOJEUsyINhB2g4tCO0jreX5fGZ+h4iSoH2xV5n5FNBvfaUELYfKCO2EjUGfBCx0mdA+GFFj7r0fm0KS9J3MyUXFgXW99tUfDUo6Lnv81l6faftnXqoTaxsDoBDaKX0Hw8JFBgAPAfhJH1kWQ4YRTTIg+GWtAFDJzN1ENBbaVnpCvJ4LYNI3Z37LuWHal19HSMl/h9rRdd/ue+95/fK/rkD7EaI+EURYtisBbI1FEq15+KECANVt5y/a0grH9uibMRDRBMFCa0ZVZu46vnnvnc37GtoBNHDY6cLztkwr8HX7vtb4WtuN6Ytsr5vyjN/Tw0KMB8PtwhgQzOxj5jehdQjKER50C4mO1CKP651lOcv/Y2fjTpvD61AudF/AnuYqpd5RZ3t26fPXM3NlLAQT92eIGtEYP0IFAKvRliy31zeg4+Jlpb3+su1qzalHiGi++JtORGNJa5dA4e8FoHgdKjrOtCidta22yetm38rMuyMRzNvpO+Q87/6sbUZSkSnP+CkA1SKMNGwYsd2vw8HMu4iojIiSmfkgEa0irSO1HwDsmfYSZpYPtx1CtjkH5Zl2QPt8cYdFhARtIKLiUKt3kBB9M9KAjLTANRnICxQwE7ROhqkAxgJIFjyj0q8sXiUZZQXMMJhkpBRnBD7TvEgTuerd3/G5/FbLRLPsueIBhqffWS+MeEkWClGgYSCtzG4HtKB1AFVEpM5KL0O+JT9wTbcESKH8pxBR2oCDeyKYxBmC4LpYg4OZG5j5ODPvE397U6dkvWUrSlNtE9KRNMbW672hIKKJPpe/3FJolrtPOGEeHzyxZ9j7115TJAMA1vqWNQOY8/UZ33D8beU/nm67o7MS2hPrQALCIiFzHwAwO8oIRI8kzijXNaj3knZkIEkmaa/jjNNrm5UUqcHwsGHEK/594Q+LnpiWbsyozDRlWPItY+QsU5YKjWQvQjsDIFH1lUZoR9q8M9j3COW/h4Woh3Up1jMXwBUAmZkfSE0r/Pf8v6B3vzNdwkOx4polWdsdnZtUv7p+V9NOxeP3YFXuahjIoGv7gb5ARPkA0phZtxTlGNexGNp5BVOh+cJaAtmt6B0uGjZcyySrhNZ/Ai6/C/tb9mNJ9hJAx0Yq/UHohVeGOlwk5paguXXOAxgDLdctYUXL8eKa08lCUAVAJSJYDBZMsU1BTftRL4ZI/2DtiMFpIqlxyCAc0qug9WwzMfOOkUww4NomWQ+lONucrTKz42fHf/r7IVzDO+hp4SYUIt66GFpIrX64t+vB4prdLoHIYaT051MmQNs+hqSbDhFlQTu4NaH9L0QZ3nxonSd79BMZ6bimSRYJQl9ZtvxLv61FGAFj6LM/2DlLoPWdjflszVCEW5S1zx555tLrtYXQlPt4Ts4bFrznSAYAk6+7qzRlTPEOa854C8JM+QQSbQm0FKSYwlcBhMcruy91qB2nWjw+l3fp2eePHtRjrUONayasFA3GzF71OUdLg8Xd1aaYrGnA0IRXdqZPKPjYmocfWoYY/GEhCMYr2441wtPuVvKWjgcZpPswjKGhePCeJBkAe1JGvtJ+8TQUczIkWQESHF65/ocV41SX+9ft9Q221HH5wYMV1jz8UK9sC2GRZkFrmtcjPTsQr3Re7oQxzYKcheMSvvZE471KsioAc1LGFCvtdSeQNn4qkPjwSoViNlmNSRa57Vw9SDYoAGzujq5HiOjHYWO9AJqgFQv32F6Xbv7QWx21rSVp07JlY6o5cHnYQ0Px4L1Kso0A7iIiqy2vSOloqPWm5E9MdCMXOwDFkpEGS49sC61d/GBuQEQZqSVZ22b8x8K7JFka8qNwEoVr2U/WJ4RyXwZgs8ForpJk5emT/3rq1kQp/QL9ZlsMBCKaAqCg7XjjC5Is9ToKJzReea3hPWldRgJpje/qmbk5EfcPZMAiQnC6P+VfZHQsAnB+BLU00BXvG5IBCJxelzBHbbTZFiI3bTYG2VXoWsX7jWQGAMuYeesIWMskADYeWUfhJATvK5IBwX77E5h5/zDNT9AKbi8y8zWrZ0WD96p12SeYuZWI0omoiJnPJXq+0BCRp8N9ILUk8432E81/He4ytaHE+06SBUBE86G18kzYuUehISJXY7eidnm8tgnpnbjGrcVo8Z50YQwGwnc1TwTUdQcRWVprrv6o82yrrbO2VYFEsE1Ij3Q67nse77vtMgzvAFiy5uGHziHGHHxR5V6M3ifculKKMyYZTHL4d3xNh4hiwftWkgEAM3sm37iys7Ph6hFoR+7YxX+rhTsiCCJKJaI5IQW588WWWwDtPKW9YX9HDCa5EnE4aN8reL9LMhQtt9/XfbXZ4mrvVMypNqgOp+Joabe5OzrDY44d0PK5ovGxBY+RxnskRBQL3teSTMCenJOpeJ0utNc3QHW5kTImR86ZPrkoTDKdjNaJK5T791SIKBa87yUZRMaGNS87rsPc+4Ig1DWZB6YXRiVZfBXeoxgE3rd+slDEU+E9ioExSrJRJByj2+UoEo5Rko0i4Rgl2SgSjlGSjSLhGCXZKBKOUZKNIuEYJdkoEo5Rko0i4Rgl2SgSjlGSjSLhGCXZKBKOUZKNIuH4/0EbfaTYujSGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 396x198 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)"
   ]
  }
 ],
 "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
}