Fig2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir= '/Users/sripathisridhar/googleDrive/projects/helix/helix2019_data/SOL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data management"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "## Feature extraction"
   ]
  },
  {
   "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:09<00:00,  4.37it/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": "markdown",
   "metadata": {},
   "source": [
    "with h5py.File('fig2_SOL_temp.h5', 'w') as h5_file:\n",
    "    for key in features_dict.keys():\n",
    "        h5_file[key] = features_dict[key]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Settings list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = {\n",
    "    'Q': [24],\n",
    "    'k': [3],\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": [
    "## Isomap embedding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1/1 [00:00<00:00, 60.87it/s]\n"
     ]
    }
   ],
   "source": [
    "batch_str = []\n",
    "CQT_OCTAVES = 7\n",
    "\n",
    "features_keys = list(features_dict.keys())\n",
    "\n",
    "for setting in settings:\n",
    "    \n",
    "    q = setting['Q']\n",
    "    # Batch process and store in a folder\n",
    "    batch_str = [setting['instr'], setting['dyn']]\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",
    "    \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": "markdown",
   "metadata": {},
   "source": [
    "## Isomap and Pearson correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAACjCAYAAACZvbxxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deXwU15Xvf6eW3lv7goSEhFgECCFA0Nisxg428e4k48TPcSZO/ELixLETbNmTmWTeZF4SBzsfJybJmMQztgd7MslM7Njj5NkkMRhvqDGGBgSITYBA+65Wb1Vd5/3R1aK1YJDULbWgvp+PPqDqulW3q36699xzzr2XmBkGBolEmOgKGFz+GCIzSDiGyAwSjiEyg4RjiMwg4RgiM0g4hsgMEo4hMoOEY4jMIOEYIjNIOIbIDBKOITKDhCNNdAWSlYVf2FYIoAqAC4AbwKZ9/359/cTWanJCRhbGUHSBeQA4AMgAFABeABWjEdqajc8OEezbP7nvihGs0V0OTxUAh6YGZF/7Ufjaj8nB3gZHqK/1u0Q0omemC8wDYAMiItsAwKMfvyIwusthCHmbVyj+DlkQTbBmzAJYQ1jxyaq/3QWg4iJCI/1fBoD5d3zzYVtGnkOy2GXRZAERyYi0kFUAHkjsN0kODJHFQERzAKTNuvHn++3Zc+cj0lUCJEIyOxXJ7HyHmfeO5JqrHnpmhhr0yyFfN8KdzdHDMmvha4jIzMzB4cpd//ijQ7rYbY/9eFJ2sVeMyErXHhry0mq3z6snIhHAYgAigFpmPrLwC9vOAbgVQ22yTSO9ryBK1Sabc6HJ5pRjDithJbgTQCkRmWKOH2fmLl1gsTbhIgB3X//4oxWTUWhXhE2mC2yAXRRWevbbM1ffiojA9jLzLmbuBADduK8AsAXAPgCHAZwGUKUPCkbCJkQEqui/KwC8omx+nJn3M/OHzPwhgD0AsoloSUvNsSe76xudga6eqDBju9hJxxUhMuiGPABZCTTA3+2R1VCLvWDBM+uYeTczq4ML6ELbBKAIwFwAC6Eb7SMRmj6KjArWrf9bMXh0yRGOMfOHOWWzilML8yQSBbQfPwUtHAYiQnON5stPNFeKyFxKoEH2d38EkAhragXM9pmX8tKqADhCfa1ysOccWFNH3KKM0n3hBqCYnQ5klBShs+4sNDWs6McnHZe9yIgop7fljTMAKdbUxZDNudGPLvjSiMhKREu8TZ61vrZamUiAaE5BoOs0fO3HZG/z/rVEtET/qSSihUSUQUQUe50xuC/6u1gSCBkzpimddWcCPeeafjqGRzFhXLYiI6IUIloBINWZs/4J2ZKnQXcrYJAhTxFmEZGLiFwAZgHY65hSsd2WVarItkxIZiesGTNgy5ylOHIXbI/aUsy8B8BBAGkAKmPF13Z835M9jXVOf1eLzJoGXKJtpRv3/V0sEW1JKy4sq/75vxcRkTXuDyvBXHYefyIyA1gCoIeZD8QY/dGRGgMI+rv33li/9/N9MUWPM3NH7LXG6vlfs/HZagAuJdCHQFcLWNNgzy6AKJvdb//kvmWj+G4EYA2A3czcd7Hzk4XLpiUjIpGIlgEoB/A+Mx/QP6oC4GBmOdh3HP6e/eTv2S8xh7+MyMty6z8dg685aJTZb7SPILTkBqDIFjucU6bDmVcCf1eL2nr0wzNElDPS76gPDnYg0mI6Rlp+oph0frKKGa0DDGlNC2w6eGpaFgALgA+ZWRlUxBXy18vhYAvMjjkQJDsQ+d6z+BKacV1Qo/XMbwJwN/SWkIgUR3ah15Fd+O3D2JJCRDMANDFz3Uguysw7iWjVnA1L2rKXTr0fMYOKd+57Nen8aJOqu9QF1t99BUK1ihpuD2jsX1LX+Lmjg88noswp857YYkutvFUy5w5whgLYUrt93seKZ7SZGINGlIf0w/MwzOiSiPIQcZP0MnNN9PjFPP6rnr2tsLu2rcZemGqRbPKArjzZhDbZRLYZwIaQUi8r4QaY5VmQxAwFwBbPiex+wehe9CUAOmdfU+PFBeyq2u3zLvgyRmuPxYwoh5T7ONcFEaUAKAOgrv7O/S1mp2PvcNeICm3Vs7dtBrCh53iHbMtzQLKboJ+35Z37Xk2qmOhks8lcQeWkzFBgtyyDJGYAMU5KfZS4CBHH6S5mPqwLaYhd9XEC0+l34Oq/X6qPbFTlmLmHmT8A4Ok6fe6JrtPnnGFF+bhruADIKTMz4GvshepTouclncN2UtlkvuC+IwLZFllMswd3fW4img4gH4CHmb2x5XRBjeivO9jbsDIc1MM6JEBPrpA1NXAN0Q2VFypXcecj1wiyWQYA1jSIJjNsGXmy7hq5KMwcuv7xR4uYWeptaAYJApx5OcBQAbkRiWnKKTMz0X20HfbCFEWyyknnsJ00IiOiotz0v/tTbvq3PouIG4IAKGq4o+9My1ffAhBi5vcuVP5CAfJB9ygEkAcAs2/ZsteWNacM51skICLoHbpvbFjWbHx2B4DSaDk1FEBvU50a8nadIqIKAAeZOXyRr+smokUpU6fIIa8PHSfPIL24UCGBYgU0YFCROjtT6TrcGtDC2uaLXHvcmRTdJREVpNrvmJWb/q1/QaTOxKxwX6Ca23te+Gqvb/srzHzuQuWHC5AD8JSuPVRIRNOIaJne0oSjLg1rWvE/YpjANi6eiTEgIC6ZLEpKXklv1qzFDwOoBbBId/pW6rllA9ANfof+PTWTw4a0ogKl4+TpwLnd+7dEz9ON+wFmgHNGRtmhn1VPTTaHbdIb/vroK2dBSct9iIhDDoRqoWl9sJoXKkTCAKN/OErXHtocLQsASqABaqhNVYPNrzTWPPQQMzcMVy5Oo8th45VEJCHi04s6iGvW/agqEwMHDZr+8yKA7/357zZlIeJkPnGhe+sO27WI+AoDF6vreJDUIiOiXABTmfmjihmt1QBc/uABSGIOZKk/Bun2nMj+WO956dpD1czsCnqPgLUQZEs+JHM2ALhrt88bsec93uiZtmUVn7/9H8xOx6dsOZmSbDFHP1YAbNn22I8f0M8tApDLzBe0vXShXQvgHWYOJbj6FyVpbTIiygJQqOdaAYDbF/QslsUpUozALikzwd9z4Ag0ZbHZUSrpzthLLjseMLMG4MD1jz9azMySr60DfU2tsGVnwGS3DTD4mfk0EXUQ0bUA3lv561tzMKjVZOZ6InoLwHVE9PYwDupxJSltMiLKADAjRmA41XzvH2Ux1ytLuZdsIxFRFhGt8HXs3GxNW9wrSPaR2lfjjZuIFHt2JtKKCxDy+tB9tlEd3Goxcy+A7Wnzsm8Mdvr3Y5CtuerZ2wr1aMZbANboXfOEkXQiI6JUAKXMXB1zbHFP3x/3ydKUBRjk7/KcyB5iIxGRSc/AyGLm99rqfvEhRucrG28GDBocuVmKIyer96N//d0O3XToh5l5/reuvlb1K/a+c8Nn0Oot5HZEhCaO27cYRFLZZETkBFDBzO/GHFsEoJGZmy6hPCHiOxIB7NEf8qTiQuEkIipFJJ3IvfLXtxbo59wLwB7qCSDUFYRjWmr0Mu537nu139bUW7JrALw1Ec8kaURGRDYAS5n57ZhjiwE0xApscIAcwCbPiex6IpoJIBeRfH3f+NZ+fCAiS+rcrJtm37voWXO61Y4YH56/pQ8kECxZtmFDS7q7ZPXC38w+KpiEAc9vz6cOJ7RFn7C++o5Cb79YFK13T5Z51fa24Dv/Hf38YwQ2YBaPorZ83mm780sAPmDm4+P7LcYXZg6seva2a/zNXnvPiQ7ZOT0dJESSca05dvSe7FQlh8krWaQhtiYzK/OeKqnzHvLXOCtsFl10iwDcXfny3IpECm1CRKYLzAPAoXFI7grtWbwi6/XPEQm7ANQTUSWAs8zcPKhof1xQ0wLwhzyyKGTYS/J+d63nRPYr4/w1JgqXNdchmzM1dNe2wTkjA6JJBIA+Z0n6c/ufeO/dnqPtZ3Hf0ILWIvO3zPmyxXvQJzvL7cBAGy5hQfWJMvyrADjC7Jfbgx8gy7xGIhLSALyeb73lRgwvMEAPCqvhVvhDHtjMLlhMs5IyKJxA3AAUQRKQOicLvXWdCIfCCoDn3rnv1Qd6jra/CuDqC5R1CbIg22Za4DvV76dN+PObKJG5AuEWuSvkQZZ5NfT5F9QZ2ltekf7US7cX9JouUM6thtuUkHoOdsuyaLmk8XeNE+cnmRAhdXam0nO8I9B9tP1nQKRLBXBaHygMxg1AEa0itCCDwwyMw/ObEJF1BN21wXCLmmm+KioUdIX2wipOJYs4xY4LpMW0df/610HllN9mXpjs/q6EMThmSURbUmZmlB144r0CIrIAgB7HtRBR5qDi/QK1z7DAdyKgYhye37iPLolo0QzH1y3laY//EZEhOfUqRyAL6bCI50NFr9Q7lg0qlwlgxoKSlkYMM7ocx6+QlOjum2sQGQAF9GNrALwbm/VR+fLc/gFX9x7vSZLpe8f+z5ljCa3beIlMfwirANQwc7tu/L8e5kB5j1JD6ab+FC0FwJZX6h2xma7piDhod41LZScpg4WmO2BXMPPOjzl/+celSMWlXokQ2cPlwQG+rLN97/zst3XrpiHyV9UfsL2j0FvYGthZk2VeFR1S96cZv1LvqAcAIkoDMFfPGjW4CIOzMPQYcD4z77/A+aUAOpi5NVF1irtNpgusP3erT23ZANBHD87rPDY4I+APZ52mtuDOa4loQLgnRmCpAOYZArt09JjldgBXE5GFmdsABIlo6gXOr0VkMnPCSISfrN+X1RU8AYV9coF9ZfR4bBcoAphypOfx9wB8OPgi+sSK+Yluyi9HmJmJaAeAtUT0PjPXEtHy0h8VWx2l1gcx1Nt/hoiKmflUIuqTCJG5AMjMjL5wC6bargaG98UsA9BvY8VGAILh1r3Z5rXbWgJvvZyA+l0R6ELbDuAaIvpg4W9m13sP+WsQmZ86wNuvpwatAHAqEXVJhMjcABb5w22yU8qPHhvgi9Gb7qboqCc2AqBqXrlXPbZ4RfZrd95R6N0d7ToNRk5si8ZhfNFearH3HfcL9plWYKi3v6ZwnenaysesdyCmpfvDdd1jfv6J8JNtAuDtUc6oTrkQGGZxEwDTmflkTJkqAA5V65O7lH3IMi+XMIkXfUsmmJkX/mb2sb5j/rtJjKx1qzthgZge5ra/pDiLb5T/wBoPyE27/a+pY15AOe4ie/KAuR5AhVdpeJlI6Dfm9eMAsBRDbTAXALlL2Yss88rosSstXJQwBJNQ5Sy3oa/WD0uBGb66/pCShvM9TFX6XNHSdUyL++qOCQmQP3nAXE/0N5sGTx0jomwA3mEmOLi7lZrFKXJZbH2utHBRInERkeCYHxFaNHMDkQksUW+/S5BJ5jCDmaORmLj8oScyrETDHJvDzIcGH+xRap4Ka16fSUi/YsNFCcYNPdZJAsGUKyPQENIAbI1J8XEDUJzFInpP9+c1xuUPfdxil/rE1mEdgm81X5WXalowH4PSow2jP270xywtBSaoPWFF7Va9AL43+BzZRorq7w+cx+UPPSEefyISZjhvWX970X99EoDLp7bs/UvDN/9wtPuVN4Y5dyoA86CBgEGciY1ZdlX3nura7f2Xju3d+5i5K3qObuRXtXnUtRll4g5Boh/HY3SZEJEty37EtTjzgTftcq4dgHy27x21wL6qFwMHAJEKEK2Mzek3SDz6YoFuAMuGiwfrq1XO4xFujHEhEmL459mWbbRJOXYAcnuwFhnm0liXRKzXfyEi/jGD8eUIgNLcZVL2rdtSNgsiDfCLMXOQBm5iMSYSYpOZxbSy6DoPwXAnbNLQVWn09RokfQ6hwTjCzN25y6SZS//R+kLPCe2rGN4vFoqX0BIismC4qwaA4lUa4JD647KDRyqViOzCYTABFN0kbxBkpGoqR3uzwX6x/QAWxONece0un16uFgKoKnJkVaSbKdzsP0u5VpeEoV7/aQDqL2XNVoP4c/tfUwszy6VPes9qgm2KgL5GDfY8AYjpbZg5HK8JwXFryXSBeQBsMInWUodMYpZV0CTCPsR4/fWw0jRmPh2vexuMmCqTk0jxMiwZAoKd/X6x2AgAAPiIyD60+MiIZ0tWBcChaAFZJBPCHJJlwazm28V3v/m+FDvdajGAj+J4X4OR40JMAyM7CKFehslJsREAILIJxlKM0SEbT5G5AMhNvhpMtS/EuT4P8u0LJAw09p0AtMt1hvckwg1gkT1fkPsaNDinieg4pGoZ86StsX6xR3c5C2q3Kw8/Vp1SpJfZ9PiynhH7zeJp+LsBKBqHIZAIjVWIJClhjQ89XB7c/HB5sPqmghee2zg/0BbHexqMjk0AvOZUQQl2a4BuM/uatH+KnvBYdUohAE9WiXi7v5v7R5/68RERT5Ft0jjsBQkqADCgMrOvyafdAWBDd6jONdW+8lYAHj1F22CC0FurCgBbwgHUANhiyRAq/vx5b3bMaVUAHJlFgtxZf+n7Qg1HXD3+lVmf+8SdJb/4VE+ocbkomHYr4SJJ0XA3ALnBtwv5tqsAfTbSkwfMSbXW/JWKvlyqzMxnZiyXbv2bp2zrEDFxihBZwAaNh8PIm9s/0HQ/vqxnRKtTxtWF8VH7b3v2tP3n/URZS5l598PlwWoAsqL5YBJSoqcZeWJJBDM3EtGyx6pT+JZ/sr7QVhe2Z00Xo+vVAgBiNlgcVVZGorMw3ACUjkhoKXrMyBNLIh6rTim8/UeWjQB2W1IoJejtX44qqo2o2Ea/D/vYqxmBIptRnaDI7mXRbfI2AfBqrKpCxK9n5IklEVHjfsoc8XZvu5YLQBCG9m2tip8PQvd1jmZ0Gc/uMpOZTxDREujhoicPmOsfLg9W9KlNTwIohj4MHpyJYTBhVAFwpOYJcuOhMByZA7pGINIo/Nd/fM3//MdtkHExEpGFQbHhop8ctDQgsiKzEadMPlzQV2ukmD5NT7+O7XVGvDdnLHHpLoloCoAmIhou5boUkdQSg+TDjfM7rgAARBM0NYhmjKF7HExcWrKytBuXfGXuqzfUez9anWYu+ODp5WrTN9+XopWz8yTayvgKY8D+TAAUk5V8e3+vrH/r6cC+mPOGazwumTG3ZE8vVwvXFTz2EoANGsILnHLOlwB49IC5QRKjt1IVALYoAT6ohXnL6d3qyu2bg3HN8Rt1S/batZG0nnSh506/Oc0hnhfs4JnJRjpPEqML7QG6inIBmAA0Apg/6LQxvcNRiUwXmAeAozlwRC51LAUTcMZ7PrCPyFy/QgBnx1JBg3GjBUClvi5GXDeWGG13OWB3WkEgCAAc56sWdbhOuZRNHgwmnkQmkI5WZOeHvrpNSASYIleLHfoaXeUkYsFaKfeZmvTN337e8fwzNembn6lJj4tdPVqRDRn6MkMLhtV26EPfmNGlwSTgmZp01ye/YnkZwNdNVsyHntoTD6GNVmQDNppi1hQidP/p3KYv6J9V/Whp84cPlu140BhlJj+6kN4mggkA6f1PdAD3/a9ttj/yTE169Whbt1GJ7Na3pP6hb661ZM+h7vd/C6DCF1lj2QNgQ0fwVOV059WfheHOmAxUATCDIt7+GK+YDOCe1Gz6NGKmzY1UaKN2YehCewDIBdGUpcxc//qce5+APiBgMIiEcdlWxWDMuACQ1UHw9w4xowmRXfeAUW6TE68sDAIAq5Rahpidy2IqZuSPJTduAIo9jdDXNUBkDEAY5O8f8fuMl8hCRGT2q9010O00Ol8zI38s+dkEwGt1kuKLtGQMIADgvwEog3wEI36fcRHZtlXHun6x6JUXVjjTFpaaOWwmNtYZm0R8tayzHkCFINCW3g4+HPDyfwCYDWAjAC8Dqn7qqN7nmHP8D92gFALwHOje7RRIlOY5FysaoG1tfO9UgXPlnwFsMtwZkwciuhnAPmY+CwC3PWhdXL5GfqCgVJwHPR9QF+UlE48sjCoADqtok4JaEESQRUApkjr3fX3gpF6DyUEKgO7oL689HRBeezrwpbFEBOLRXboAyDMdZbHHZKecWnaB8w2Sm9j0eQBjDzmNWWRhxqHBVWAGvGpf3VivbTC+vHk6u/DTX7He9caprA/ePJ29+c3T2XHxb465uzynAlMlQOCII485Mr2lPXzRogZJhC4oT16hkIKIX2wRgLuv+oTpnrFee8wtWQg0r14FujUgxJF/61VAFO0lY722wbhSBcBBdN7x2tulOT57v/WzY71wPGwytwpS2jRCcxhoDQMqSOlTu2ricG2D8aM/syZKU70ml1ZIw20zPSLiIbL+YPlU22w0+I8qALx/bXzhhThc22D8GJJZo2msCiJN/Dr+scFym5TiPuurfQVARXXba81jvbbBuLIVQLg/BwNQtDB8SMZ1/IloyatrlWZ32+tPuLJung7dgaeL0SAJiRr9QT87dv4xKK/7jIUBhH71f70bfv9r/5h7pLivheHKvCWXWfNkmPI+zXx+XSt9XoBBclIFwHGuLixPnS4CAPl9LCy71rQuHhePu8g2zH76tkb/CYcsmCU1sht03HYbM0gYLgCyojBE3al17mRYLnPJn4iHryzuIsuyFFT0qV1ytmUaGv0nooeNdJ/kxg1AOVcXRkGkJUM4zJAkygbgGavQErF0lJsB1SalwK/2RI8Z6T7JzSYAXl8vsz1Fl0TEVBcQh14o7iLLFHhrGoXDUwWNZ1imoCtw1kj3SXKK7WbkiNJrghBZiyzgY5gs/fmAY+6F4iqyQzcoheki3ii1lwjtwVNUai9iUk9zpsDrjdFlclLbVlAIwHPmuHb3nLmyCAB1R1QUl/ZPoh1zLxTvliyymK05W+5Q2kAEckhOAeG2L8X5PgbxowqAY//ekHS1ywoBkTm0+gJNcemF4i2yAaGJkBbCbMd8qbZ3/6o438cgfrgAyKrCsJoFhBoElM6UIRP1QZ9De0NR65h6oXiLrD80MT+lEjU9e0BEilftORjv9RUM4saAcFJfD2P6FLMy1WZ67oai1gfGKjAg/iLrj2MKJCBVzlA6Qq19Mx3zHkMkdcQg+dh0pk715RdIqq9Pg8VCKuI8UIt7WEnP+a9CpBl2/+89N771XtufXyGiq5n5g7jezCAuTJsufefN6tw893uhta4Vpu2CQJtKs87GbaCWkO2hB9wgstSnA0A7gFxmNpb2TBL0kWXVlp/23rzhIefrG+5q27bjz4H/ifd9Ei4yACCiFQDexwX2vDYYf6KuCwCOl/7VK6+4xqKkpFFfRqa4IJ6tGJD4zSKifITIFoSNRJQ/Tvc0+HiqADiOHVHkmXNkdHdpckamaEcCYszjIjJm9uv/bUFkzx6DiccFQN71ThDzK2RYrAQkKMY8Xi0ZcL4169N3LTGYWNytzWElPVPA8VoVs+dKQIJizOMmMn3uXjMirdnghW8Nxp+tZ+pUqWh6JLdH9/BLiGTIxpXxbMnAzMcBzADARDSu9zYYyLvbAw8VlUhqw9kwpk3v95OrAMY8BW4wE/GiawH4YDhnJwwisphMtOBMnSpnZglIz+gX2aS3yQAAzNx2be51M56/6sV/avu0t6nh9r6mo58MPqc7cQ3GhyWVy0zb91SHwpVXmWKPJ8QmGxc/WSxdd/YWBtXggZ2tO1OLHcWY6ZgFZqDVa+4OM5XPe1M2UoISiL63gnj1anPlNx5xPld5ldkCfcsbRMJJFfH2kyVil7iLUWUSzc4p1lyc9dVjlnMWAMBuVp09AdlY9jOBEBFdsybFdf/Xcj9TPt82u+ZI73ZCqIPB/ctCxVtgwMSIzEUEoTxtAf77zH9hf6cHC9IrYBI1AcY8gIRy801p6x99JP/f0tIka3uHKq9ekV5uFgUvgIqUdHfCepCJEJmbGUuIIHy68DN46dSLEElCkXWBBmMeQFzp6XT1Jyu0tSt7169Pmy6IZJs2zSzt+cgLs3l8Fo6eCJFtAnAPM1KJCHcV3Y0X67ai02LuK3HMN+YBxAldYB7oq5F7PL7Fs2ZaaH6ZTezoVCHLA/fBSmRdxn10mfY7Zz0RyjXG82oYzUFVar42854XvrHvs/9cts1kHe/6XMb0739VdyqA1lZFWuZyUG9vWDt+LIAF5bboeQmfSTbuo8sLQUQWAP/rb0vuPfbTyp/dCT0fDcCmtN85jRHnCOnpdFUDcHV2qXjxxTZs+EoOmAH3bm941coUDYNGlIm0yZJGZADwm5W/LW0JtOxZlbPaVOIoGfAQDKGNjJ5O1+aGhtCGP/+lW77l5nSkporY+U6PumZ1youCQF7E/BEnUmDAxNhkF+ST+Td+ozPUaXqz8Q2ZWcMM50xjR5MREjX2Dx3yXXfqdJDXrE5R0tNFecfbPeqqVSm9gkDfS7SoBpNUIgPgSjely5+Ysg7VbbvgC/tQnrbAWOLgEoka+/sP+JyKokkFBWa1uNgcfnNb1/FVK53bJJEeH2+BARMTu/w43ACULHMWFmdUojvUjfdb31M11gzXxqVRtau615mWJkqhEGNBuU36YFevWFlp/+uUqXu+PhECA5JPZFsRWdeY86x5KLIXK3aTJfBa73Ozeh484g5trNsc2lhnxDiHgYjE//dG53VWqyC1tii4apkDez7yorTUKmVlyhOajJA0Iuu6s7cQwBuI1IkA8LSUfC0rn7E0f8G1u5s8SzsDXRsAeAyhDYSIbADWdXSEa3NyZKWy0oGDNX4UTDUjI12a8MVukkZk0P06naEOWR/xkuT0yTPTi20as5xty8TJrjNyfU+DsdZZDESUBeAz66+Zarr3JldXXq5FqD3i19LTROTmykmx2E0yicwFQLaIVrzX+i66Qp0QTGGBCEJpRgnSzSnwqQG0BzrlfS01qye6sskAEc0EcNd3H1x46H+eW/e8BOvd5w7mihZKQ+GUVBWa+BIS7AO7FJJpdOkGsMgqWuWVOatwoOsA2ptD2vzCIhBByLVnI8eWhT3NB1RPy6G6ZXTLEmb+cKIrPV6op+4dMGl6+e2vf4jIKmK/+N5DC38GwFHf0CeHFQFzp80U0IcwAK9U/NyE+xeTSWSbANwNPRRSnlautIQa/DvqPxBWFSw1S4IkE5GyZMoC74dNnu8CmE5EtzTdv9eTbkl9BDHORdNPpk/4g40nusA8AByKosm/3Hq48ol/WNq3YknufKn4OYXn69wAAAa/SURBVC108ouu3Z42OT3VhDkz0qLFksb1kzTdpe7Rr0BkJRk3gC05pvz5C3Pmlf3+6J9eq+uu36t/VvHNv37vADO/dlPJddq/HfjtYb8S2ICYPbIvw4FBFQBHzdFO+alnD+Lzd8wQVyzJtQKoIqK8/3zt5NkFc9OV2SWpsWUm3OCPklRhpY+DiGYDcDLznuix0Ma6zV2B7g3P7v+NvHzqElydXwkiUgBsMf1k+mUTIWjee9eHz//ueGXpjFTcdF0BiAjMjB0fNNVcf/ebn1LqvuhHTMYFYsJxydBdThqRAQARORGZu7mHmb2hjXXV0LuEDxr2YH/rYSzMLsPc3OL96TN7dyKmCxW/sWbCH/alEGt7Nbf6966/Z9uBdavyP/PDxypXSKIgA0Bjiw/H63rUirKMX2fMf+n+weWgf+dkEBgwyUQWhYiWAOgKfvvkg4h0kTIAeEN92NO6Tw1ntrDDYqEcR5o0I3tq/191sgstant1dgcdNUe75ENHu8KfWJXvLZnmvB7AGxqzw72vVU5PNSulJalJ01JdjEkpMgAgooJPFK1c+tqnnv9XkYT+boIyO8KH+w6K3qBPLsnMx4m2cwCgzs8v+VXqxvVfn9hanye86wdDWp569nz3WF3PveeafFJxoQNXLcqGJAkKgC0/2Ox5fnqh4zt3rC8qtFqkaiRRS3UxJq3IAICI5LWFy2/95bofri9Jm7YAgFsoOrOSBCwMKEF4zp2AQAKYGdWnDp381u83fxHAaWY+AwDhn7895EWPR2unC8wDwMHMsvvQWfV9z+mQOaurdcXSnKKy2WkQxciYLMb2uoOZjyW6bolgUossChGVARCC3z7ZJeS2vA5roJwIFHOK0hPoezbjkZu2ACgGUDgja6rzDxt+8GhpbpFNIBpgLMcKLfzq1qFCvO2ei34+XEslXvX39QDQ/MZ3njl8qvXLx860SV5/CEvmFmD5gmkKpzUch9lXCoqM+huafThxqkddND/zV2llLyZNKzxSLguRAcAv1/1w/oy0oveXF1ZYrEUdMgQGEcAMhWioeLxPbfvFzuOer5xoPSel25zIdWbAYbaqqVb7v8/753u+DPQLaNhRm3jbPfUX/NzRtR72njegt1TnWnuUv+4+EfrBczt+ePJcR+dT37rp4dWLiktKp2XBYj6/xSQLioqseokRRrWnFZmpFsyantINQvlk6RqHI5mcsWPivgV3bWBmi7txn5zWY8ec6VmAOcTQhMOwBW4e3A1aZfOSG+a6JMwFuny9aPF2oTfgk6pPHVpVRrQBAD98x42fK8nNSZk9NU90WM0AIDPD2dDe+eSnb//C99/f9I9VTqvV4Q+FZGZGQFHko+caUzrU5pcdTiGFCKLGjPysFPm2NfPwtzctzhOv+vsfhnf9YB5iBiw6GmmyULffibM957B0cT7MZNPQnvqKtOT7k1ZgwGUkMgAuIpKX5S9Cu78T7pozAEAOk40rnr9+uJfkRmQ9DjnN5kSazQkASuW00jfv3fqjLQAQfnXrl8NhTTzW2ARvIAACgQhShtMxD0D6ufZOV4ZTkW1mE4gAp9WCu1ZfLVqmtKdBVgav9h3rge+PboQUVT5wollVVY1AEKflpmHV7IWR9jLiLJ8X5+c07lxOIusXTaY1HZnWdABQTnad2U1E0ZfbzMyn9f8PCGMhEgfUMHDpJLcoCovmFOTHtjgKgB3M/G741a1/QWSVooGfa+IJQJk75LjugZeu/ofwTStKv/DFmyvvyUm3z1s8J3+n1Szb9PoMW2Yyc9nYZHooaVj7KRrL1BdJnqYXaVA378gH8DYAMyI5bAOM/1HYZBoAhhj6PTKbbgDBBkBW1LCy/3iT//vPvnX3H9+rbQHQyMwDWtfYEeeQe+kDhslK0sQux4oupAGxT8QIDACYuYmZ3czsBsCveN7ZuKuuRmrsbo+ORAfszamPIodcMzq6jPn8JUTW9gIAEWHTHc3H7PT6zto3t+85WfPGB8deBrD49XePvK7ff4hodCENvdckFxhwGbVkoyH887erAbhOdzShpbcLQmS1QfhCgYNrf/bgF2NOpeHKI9LF0st/99AjeRlpnyKCpGmR55mV4lRKpuRsEW+757KJoY6Wy8kmGw1uAIuKMqbIRRlToseiNteeCxcbSPjVrcUY+iyTJtVmorlsustR0r9Nj/77aNOVB+xPFHOtSW+0x4MrWmS672yoHTTy0FK8xHpZckXbZPHkYuGnKxlDZAYJ54ruLg3GB0NkBgnHEJlBwjFEZpBwDJEZJBxDZAYJxxCZQcIxRGaQcAyRGSSc/w+4eqoMHxWGVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 396x198 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAACnCAYAAACB6MdcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19eZwcV3Xud6tnRiONpNG+L6NltFqStVuy5d1hMzFhs42JwRAIxDxCHg7kEcAPMBAITjDgACYBY0LAeSEYjB0bY2FjW7LQvu/SjGak0S6NpFm7u+7745zvVtXtWRTZGblFnd9Pv1bXVN26VV313bN85xxjrUUqqbzWJbjYE0gllfOR9EFNpSgkfVBTKQpJH9RUikLSBzWVopD0QU2lKKSkJ082ZECFrRoxECjvpVuMfORz8tmeBcpKAQC2rEz2yOXlbzbUffXT6LGleglG3jkbGJhQ92lu0WPyiX1tnz7JiamLzrS1yfdcHmhpl/9n9Dy9y+WzhOeT7TaT0XnqNfDc1gItrfL/rJ6/vFTH6p04LwIdS8c2oXXXbAMdn+PyXvFYbj8n58qdle+ZUr2mfmVAn968UL0Gvcd6Xjd33i/eW56jtETuCeB+H2RzyX31A6VlyWvh79aeddeZO3gWAHC8WZ6Dg62H9GAct9YORQfSow9q1YiBWP29uxFOmyIbAnm4zMmT8nngEOy4UQCAcPx42XbqlHy26E08e04++UAPGybfe8lF217lMC3NMvz6jfK3xib5HDEYAJCfOzc5saw8lMH+/Tqf07Bba+VvFfpSza6W8YcM1vPrD9K/Uo45fiw5z2wWZstuuZZ6uYZg+mj5Pmu67JvNyhj64NpBMrZp0we8rQ22bz8dtzlxr6DHgud7eRsA4OSz8r3vaPl7r2vHIZx3mRyrDyTvme1VrmOekPlt3CJj6UvPc9hRw2COyTXw9zFHj8s+pfrg6kMYjh6l92lo4lrMwYPuNzr9qeUAgH9ZPxEA8Ikdn5cxkNebXig9+qDaigrkr1gMc1ofvuN6sfpD5ZctATLJKQU75MfGWfmh7Hh9kCdMkDH0B6RkXl4FW3NY9p0pN8LOnimfvQVJeX7bTx8C/QHNzhoAQG7DIZTcfLmcZ/w4GZgvQkXfxHlNs744ihy2okK2/3w5Tq+QF2TAm4fL36ZM1GPlYbKDB8lnn76JefF+2L79om36UjtEVSTP/+dqAMDG38pDPmuJoFXpXdfI/AcMcC+A1Tk6BOVYRHi9x/kdR+UUM0fKcWPHwg6Xawj2ycuM9qzen7HyvZ9cAx9+NzYlCNwLed8KuQ//WCsP6FenfRYA8Ikd96IzSXXUVIpCehRRYQAEAczp0/JVl/FwwAD9fhZ2wED5/9EjAAA7TJfaMSPk+1BZttwy1re/fG9SZMtkYIbLeOFwVQv0PEZ1MDtQz8Hl+rgsfSgRlCq5phrhhCrZlvOQgdLUpOeVTyJpsFmW4LAtFyHpvOnJY21Sh6Z6x2snGpm2VpgzjToP1RGpk760SaaxX/TIOdfLNZRcJWqV1WsxbW0RknL1UT3SqVUHVUfsLatG5rqZifuBTMapOuGU6sR1O/TnqqArjDl8JHFfNnxgC/6tRq6PSPpX4wVJ3z93HwDgEzvQqaSImkpRSI8iqmlpRbBjB+xAMUDCKtUzz4lehcAg2CM6qTkkb2Q4Td5gVMoxRBRalf7+dlIVQjUWTKOikaffmbNyvqC2ToZcs0e2L5gs36dUO/3VzZ16rBoeTs9UpAk2bgUAtC/fCwAoWzoW4cI5yfNTF9S5WyKbGnMO8RRFTXsWaFXvQ7nsi6fXAgD2/lLGnHCDGDPm9htk7jpfO2RIdG4aR/Qm6H0x53QV2l4jx1w+VT4HDdQxxCAKGg7BUkfXbU4X54oTt+7jnwNkfl/f0R8/OiJI+nFF0q+sUqOypUr2/ckP0ZmkiJpKUUjP6qjWAm1t7g2n68K5YM40wvYXnZOI4Nw/ijZEMOcdUD3KVsoYtl8/h3YO9ej75Nh0LQ2QsQMiKS3YIPb+Eo2IHIp+4SjxPgSHxcNA/a7szeIKCidWRX5ZXj71S7rlFJWcq8n5NfNuux2h1vY2UeDyTbLvxDepDnqNutrOclUKEtcKGzqdlJ4LJ3RtXTFLdiUKK/LynocjRsC0K7L71ryKuw96n47cK8j/o+3iZvzRkS/iT4f/LQDgXVWqG585o+cL0Z307INa3gt2SnV0M3WCbjk9fsL9MPTPOWOC/jp9QLk9HCvuIzOUzvWsOx2XJ+ek1x8ms0v9m2PHyH4jRiT2NzEDyqkl/KH0R8+s2yDH7K6X7VNlHuGoEXryEHbgAD2WD2Lyup14LwG/24o+ME8+DwA497zco743yNIbLluQGMupKhyL52rLAjl9cfVhNgcb9Lp1Ge9PtUoeUAcceu2mqanw3vCTRpXno/3Carkf366X5f4r0z6Ljz+nRiX9s2pEn4+c19JvjHmLMcYaY6Z52/sbY+qNMd867zOmksoFyPki6u0AXtTPuFf2CwB+d74ns8bAlpUhOKTukJPqzKaiPnZM5FA/p8sCkUGXCYdw/ZPGlVPmwzDapqFaU6/nY/SkepJ855JIVYGoEASRYeOFLoMd6vAu01t3jSBbfqSqAnUHZIyysmh8L5TrVB6uKFyCuayeVZfXz16Azck+fe8QQye8XJZpBktcVIuIRwNSVxw7ZKi739wnXLQweT6uUjTqiOwaVEAQODR2K11Dg/sbAJz67AsAgL9fJc78b9d/AQDw4TGfAQD86fRamKOqWmhUMaSqken+MewWUY0xfQFcBeD9AG6LbZ8PYDiAX3d7llRSeYVyPoh6C4CnrLW7jDEn9AFdD+B+AO8GcOP5nsxYC9PeDrNLHLy2XlDBzKiS7zaMQob65gbPvyR/O6xBgqkSL89ffZV8J3IoshjAOceD7bvkbzskhGwmaRx61MjkxEpJ0tD3trTMoSxRKLNqHQAgt1oQs+TGGYnzOiQlwg0eHAUDhonezaCEPacufurT/NR5hM+I/rv+VwMw9/US27eT5B5Rn7fDVL9kvP6ohD0xSMKyUNUxqK+H2SEuM8xQVx9XCerwp7wwLUlDNOoQCzd7+3IeD64VJP3KftFJPzlBXFBf3LRA5zc+0qNJbFFh4KUrOZ8H9XYAD+j/f6rflwB40lpbb8ie6USMMR8E8EEAGDe2Q2JMKql0K10+qMaYQQCuBzDLGGMBZCBcsRcALDPG/AWAvgDKjDHnrLV/449hrX0IwEMAsGDKKBu8vDpyYi8UVArHCEoitM4pbujuoF43ik7oQckT0BHPkGZjI0yDOP+dtX2lEkyUWBEhmDdZOuCbziF4Sqzt3E5BsLBSds78ieh3+WliVzo0cGikCJfLRa416tutbTp+c3Lup2VVyC/fLrvVyvb5d7QANyyTfanHenN1nwxD792n81KWVe9yhIvnyf9jDDMg0kXtYGVt0YLnfVPEtb3K3YrBe7flQxLCvX+r/C6PHLkPAPC5akHSu+frPI4qqaekBJZBG7onNdwaHHQ0v06lOx317QB+ZK0db62tstaOBbAfwHetteOstVUA7gHwSEcPaSqpvFrS3dJ/O4CveNt+ptvP29p30paF3dUAM1fpd5Pk06FFSQbB2vXyt90H5XOxONCt+jwdCZqWqa/v7NqHcKdYpMHlEqK1owWxqYOafD5xjPM0kNiSzaJ9oyDliZ2CQsOvDxJjOSRV6mA4NApaABIedbxSkk+I5NQRFeHyzwmSblsuTvrpy0QftzcvkzAqYs5442EL7wPvIXXHvXIPzIQRkT/Z80G7kC5J51zNvECJuy+xeTy0S1a2R46IdX/n8E8DAP52taxe5ujoxJjo1y8KGjAQwjmH3deW6PJBtdZe18G2b3jfHwbwcLdnSiWVVyA9G5nq3Qtm1niEk9WPyTeL/sT6I8AJ9ZeOUlIxCbnU/fhWUidUCl9QH+k5wSKxbklJY3jR8O1W/6DT1RQxgo2bZfvmWgdcI96kFvA8pc8xIsOUi75J8koi5Gs9P6mLFsl122fWAABa9sqYM14v88hcI2QW29oWISRJ355/15Fl6Oekb3SxksWDIPI+KI3Q0Qn9kG4bqZNyTcHBerdf7f96GQDwza3iOXlQ/aR3q5/0g1M0S+PoseTY1HNLokctYEbH3hrZUOYbC4WSklJSKQrp2VSUvhXIL10S+UoV2YIXVgEAwt1HYRaI3hrOnS3HEEkYzfFIKcF+8ZHazeorXTAF+cvlWIcy9N8RwRld4Tw4v61i2Z58vg2D3z9J5yGRIBfbpo9RETVo0KgP/Zcul6nZzdHNmfv8TNT7fc/I9qorlWBy9zvknJxXRUV0/URFxuGPqWeDFvR+oSzaCUKsoV5qBw6MCCM+L4AeA6IfI3S8T4rOtn8lHlQk9UnPf79e0N9FykiG8SiNCILofjCKptRFl4fWhaSImkpRSM8Sp8MQpqUZpk50H1Ov5OjdElUJxgyI6Hqqk7qkOdUnqTe5LMgmRQH1JITVkyI9ztNnibCWUa99NXLMy8Kmaj8oqD3opkrYaUL9o1Vth0s0y1n77UkGkLOQFY1MS0s051USaWp9UdC39ZggaNUVMkbJjaJPhswwpZVuItJzgU+WVLycnC+cqno/dXpN0UG23ZGe4Vv/mu7DNBOmj6BSjt31YYnGPbBtKL5zUHRSkp7vninXYhqr5BhG95jmw3sdty2aTiXOF04Ur0ywpYscFJUUUVMpCul54nQ+B6yQtI2WrWLl9r5KOIzhG66NdDIiWUmS2QON4mC/Wvnj5Nj8okWyX0tzZM3TdFe/KZHU5eTvFL3u5Cr5++A3qlV8y40RGjPVRCMzpk78u+EEIQS7SBr1YdYb6FPhCj+0/k7mumOThJBnLJbVIPgr1UlzSb+u06Hr6gDyEjSqxdUh2C2Rn7BafdFES0bXlBOAbHvkoeDfmIDXrHov9cgK9QkrP/fBHTLf7xz8Aj40Wqz7r6yeqdcrESfq5s4XfTrJiIsnP9Lad75V/g5DB6I76dkHtT0Lc6Ae+WOyXJcOUxrZYC5TWZi8LhkkgyhBGcfFCY42fYD0AQ2rtFAFl/tstoAoTWMh2CM/bn6FLPUt+2WswTdqJutlGogAXC68y7OisaBZAY5B35IkJVtSC59dhSYlO586Itcy+2bN1Vo6S6eljm7SCzXESBecrewPc0jpdM4oUarijKn6PbkoOmINs3mbzzlyiXP0c84n9J7qPFbcJUvwI2qgfu+QLPcfGfsZ3D3tWHIMGkZc2hncoHjkcNPailCpkJmNWhjkkKpRZa8CzS+VVF4L0rOIerYZdvl6lMyWpSW8dikAIB8vW6PuHUdi3lYj+xwQNCqZI8fm52qoji4XLmuIBQWIQsyrUhfW0Zfl/Ry2SJfpW1+fGCNBwmA4kUYTM2eJLCrOBaVqR/NLx7F6m4R9r7pW1AW87+2yL8nNrepyY/kbLtd0gQ0YGKHfXqURTvDyuujqqVX3HMnRsbI9xjDfiTW49LoH6zVpWPT+7UIKf+ykIOkHRsly/497rnQuLuf2csadqjx039F5z0xbqnLnzkaGqCJ8fqusFplFqr50ISmiplIU0rOIGhiYfr1gJyjBhJXwiCRlZYDqOsFeebvDZkHYzER9+7X2lHNWkzBBhGlqisKLrMiyV4ym9hpR7IctkMvOzJF5OGOGkbxsO6yPDKTCsW4WEZUOd3WX5Z+XSimH9ldi4QxB0tKlokdbko7pUNdzuGAGdVTW5IqJnVwln+r+cQYJ7wMRlARnripnzzqjzo2r973xXqEyfn+dINpjJ/8OAPCWQUKEu3MCw9P1AKmJNB45Z4/I4ow+rnCxmg2GqUenxcWWmaaVbGxSn+1IUkRNpSikZ0OoQwYifM/bC/Q7uLBoBsE6cV2FO0UnChaJIzu/WAjLzrrPJ9OCC/RSAGaTOvLXiyO7dI7WCnizVhWh24roxLFOnIAdqukjTBk+wfpUJcljFFHzTwmhZc2zom/OW3IEJZ98q+xLjwBRWpHT1dlSJ75zcalrJ9ixC6GmoFD3deFHRWeHsP1ien5sXgCcTmj7iguJIeUvvST39ms1Ehb92nTRST+2UmoFkGBi+/WL5sykSi/1282LercjoSuh/UwjsGGn/G26XJMlacgnhXcgKaKmUhTSwyFUK7qN50w39aLLmboGWC26GyzTNJWxoxNj0IqM6ikl9Ztg1TrYekG/bI3WCp2lpGOtl1pAL2SKBlGisrIgDMvQqLN6H34SAHB6jYyxt0HOMe9KQe+y102NrHuGVTlJ6sQ8H6l66oBnqk44ZbLzALBanqP9seivvzoRnag7NjfB6srxxFtltXqsTpD9B4clfeSeKgmLvm+eJAGa4+JZMKdYu8tIACM+Z947BmIoJLzQ/xzX5RcpCd6rWEPyeVeSImoqRSE9H0LN5RI6KQCYPeIDzO84iszCKgBw5dOjWv7qV43TxoBIZ9Tv+a2H0LJT3vKKBYJC4W23yBg+ra8yRnIGolBjEERUN0dc1jRoffuPvCjbX6gTxF86WpA088X3yznjKcDU76hHMolO9TlXK6tG02+0FiwQkUuM0gkd6YQI62qYJn/KqJxRxo3xT7skIe/pRskuumuEpI/83WamNI9LHOuqR/frG52PCOoV73B+VXpwiPiH5T6EM6a564VX6bvTGrQxSRE1laKQnkVUWImOsP79CUW4dtGhMnPGwI73kviInEQdlz6hJXY2SYMEu1fQqL2hDX1mavxZi1U4FCA6xml0iOm9PJcNordc3/5gtVDesr+VohZHGgUxFo8UiuLIBYqO9Ovmc5HuRR2UFQmJPiTJaJ1SW6XzZRx9SKwOAvU6j0juipTxXNT/9RoPvOtJ/Hy/EFuebvwaAOB1lfcAAN4yVkneLo2FnU8ULbla9ekb6aIkmbAvAK+lLbk6gH0SRiupJl5JkIQjnbMr8dSFpIiaSlFIDyOqEV/pSkmJDms0yetyidyEC+dGliB1IBXnt+RI1Fl3SdSpdZXoRr2vHYn8La+T8ejDo6+P5XdYNIwUNacHqw8wCFwhCVq3rb8UZtHa1YIQi28UJDWfuUs+vbI0dsjQCEGpk7HSM9RvqgmJrtQmj2U/gNoap0e7MbzVISpEUZo4luf+wroR+MFhid3fr37Sj26WjilO7yXdTsvz+B4Vc/pU9H/q7kRMrk70bGzTLjZaPiheudslQjISxRWtJGVPpXKJSM8iamurVE5uV/1yjiaiTVIScjzmSz8ckZPJdDXatKxei9G2ip5VPl+5ANXjojGor/VKJuTxLXe6KaNcROAggPn5bwAA51YI+h2ql/HnzZHzliwSFlXo8TDjXgpnNdMiJ7ooslpN+XDH8u+KNOGoUYUlGX2r20un3vOG/wAAfGO7oOMPDt/nrPs/maBFKRgJY7ye5SV1aJdOwt+jqTXSSR1xW/oeBGrVo0K9ItMmRXNHLNqWbXc6r2NiMQK3bgu6k559UJtaEa7eheAyeZjyV14BIOauMUF08/gjU1x4U5ecOiXyThXjK79Uxzp3NnI3eTnlbrvrGBgjw8SltBQnn5F9t9bLUr9gtiyTZV9/H4BYpqjnYHcPVq41esl0eFbcc5WoaTQxI+GIZov6DnEgeon9Fo507Osy/oM9QvRgpecPj/kMvrllbuI8ri4J8/td9z2PhM2WPO3ZiGRCFYfq1H6Zs5mkzdOoonC/eB1bGsKuuYbcu7bnDqA7SZf+VIpCehZRK8oRLJ7mmjo4IomrO29cdqdh3v5//RYAEB4QBAv11QoW6BIzRbJFnbsml4uWei7poV4mkZvCv+d0OfuN1GJtWXEULS2ylC65QRCjdJ7MOXQpL6zKrGMz1MqxS0qivqWKMuFodb357YS4imiqBl1cpqWloE6VC//qmL++TmoEPLRHxiRV7xNan/SuyUcjl9YJbwWJ1fwCAFumWbO6NJN8jYqKSF3gPSTCL9FGFX7NVd5quhVPnoxI3hrYaH9GXH3nDqaVUlK5RKSHm/b2QX7+vIiowLfQtRYPIiezvvW5HaLrtOwXZ3TfBZpqMV/b1nhVPVBREVHMfCc0kczTxUhgblsr59q8cTgWvktdS++WavChO486xWmIUb/j3+OV6ti6Uo0mUgadk5xJfDTAYtUE3diazOeMJj022CcEkkf2y/kfO/lFABHp+Uu/k5XGjrw60s3pcitNIphDXBo+NCrj9Dtu8yuy0ACkjs5MUzalIwLHyEO5F8UQ27xSKZF36Gr47+hUUkRNpSik50kp2fakTgokvrPxl90rVmzbSflbxTzVb+aoW8jV7tQ3lqTkXC5CaL8eqVe9LlgrJI3sCiHFHNomVu7s+UdgZmo1PCIJaX4kgTDVmWQVj5Rh2tpcTVdHLPEq8Tlrn+kaXg2oOApx362v+yUA4IHt4jn46TFB0tuGSrOxOyfoNWswI+6sj0LJFehIXJU/oiPvbVub82CEM4Qs5DwqfoCGv4t6OFzIO5uF/cmzMuw+bTR8na44VyzscD5xSRE1laKQnidOt7cndVIg8d1uFXRr2ypWZvl1ogvl36xhUTa09dtwu7TcM9HfiNSuslySEpf7vfjvtr4kzuyZV8obnvnEbS786nRDWt0+2pBwovMKtK5WOHpUhC6e/9ZPo3HCczhieWnky1Sk/t5uQT2GRenM/+fnhNASTtAq2zqkiYeeQ71Xno/adVphjX1Prwx274sqstjkquQQk14Xz0fqQt/Pb8L+Z+X6Jr1dz389+115lWI6kBRRUykK6Xmanw0LdNRglTR4DXcfRu6IvJnlNyiJd6b6Sf1Oej46OZ9ovnOd9PHnAIifFABqdwg6zVwmiFK6QNOns9nCt5zfS5LhUOcrVT0znFAl8+5f6dCHCW8+gjricGvMBxw/lwmwfNkzAIBH9osO+OOjopOy0vMHqtW/PPFqOcSvzx+vPUWdWGtPUXd2IVPq45pm7soHXTY9VlOrb2IMOrZdJz/6ZLWoBbvL5FvzmPTHev1XqsfGRdtSP2oql4j0PM0v5it15IvdEglp3tgSdU/Wcj/Ocidiks5G3yMjQkRcIBqftfoVyVp/L0i6eaPQ6mbO0PSRj79TzhnvAk19kd6E5qQlTiQ1h5RgMWxIdCwA2DCqZMd4txfNsqoyu3borqpeRJ37/l45/6PHBUnvGCbW/QNrtHMg02mySQKOo86VlkXIznukfQ9AdGSPA9IQPV8pwjDyXLDri0cndELy0AZJHz/3rPim+75nBuwopTMyfdvnGnQh54WofndpY8x4Y8w6Y8wGY8xWY8yHzmecVFK5ULnQ7tINAJZYa9u0qe8WY8wvrbXd5BSIjhr85kUAQH6PlJhhV5CKpZWwUz3Kn/MlZhLbja/Pcb9cLtJJV0rXEfYvPbRDUGHWLEHS8oVavIHpy/Fe937CmZ43qE0yfVgn1el5Gk83iIjC8T5NIh752dN/971VfKXf3jECjx7/EgDg1iGfAgDcOUGjWYOSJYYiP7KOFbPsXUlMIpjfwZsWO6t7T1KqJFeY8l4FhSYK0kdUNw1//Gxic8U1Ms9wxrSoonV5sneAfTXSpTvqLm2tbbfWaqwQvc5nnFRSeSVyQd2lrbVrjTFjATwBYDKAv+4MTZNNe4cBpWVoXSEWYsNOQZxxy5Q99M43FvSQt4OiWvQACkv4kPTLz5KM05OyL9YAAGrWyhhjZ4h+VfpV4ZRaWsi0VPtrKkY+54pSUEcO9khs3e7Ry5yjNf7pK2WcnL5IIFoVfC6pn0ZSxmLBos99Y5voct+q+3xB9xF3nc0eStMnHVh3De5PJGqzhJCXNhNovydXRZvdW1gQ7lADoB0L/ZXOFS7+lfShWv+UHDv3Zu0+eMNi+fvJk66iNT0CBZG5LuRCu0uvtdbWAZhtjBkF4DFjzH9Ya4/4B8eb9s4fM8TaBx51tK6xi0URz8zVnCmv6klCSKfzVQLrLR/3/xhNu2T5q6+TH2byPLnhZUvUKU4DgMs7AwHOYCuLmt+yit0R/dEmCkE4HKoPNR86vjh8oeJGo4ZbrcfKpyHy6xtEFfrhPtmfy/1Hxn4GH56mZOpeUvqdpA/3wPrZuiyrTsJL7z6Obe9qvpI4o2oCayi47bwGVV1s/8pkWyLANQXOPynVo9sOye8w92YZI7hJ6te66n69ernwszOEPQOwK7mg7tLGmL+21loAsNYeMsZsAbAMwH90e8ZUUrkA6Q5R2V36z7nBGPM8pAX6amttizFmIESH/cfuTpZvAc5sthh8ozY1uEsqMIexEKPt7VXk8HOnSAMsT+Yj0ag4tLIUu4/L0jmvWpbpko++QYbSJceRrBl+pXERqzdlzqihoS0v7SSltU2oku3OuEgem8ioZJYlXTmsDFgn7h9zVAyvB4U/jCdOJw2nr788BRiwGAnxQ7lekzc3H2aYlpdHIWXey3p1qbGOAO8t74/L6VLXV1NTRJRRtSD/Y8kp2/6CkGNmXC/3KXjHNYlzUTVARUW0KvLTN/K6kAvtLv1dAFlFWQPga9bazd2eLZVULlAutLv0NzrYvfuT9TGoXFgCM1M7mTjUjIUrPad4VNlYSSc0BFjr/zcrAQCtaxVRTw/H7LGiKvdfolQ0rb/vkNS13U62VI83xnWtErWmqMu9px7LlBmnhyYNJoRhdH36GeyQ+qBP/KUg1uMHZd8nTn8VAPCmAZ8AALxnIhFuaISG1EmpV2djwQlEK4xD9uFSvypR7Y+6sk96ZuVrIitdX0x+LClxZJvcD34NAGitlXs340atHn3tdDk2TvpO3A8LlKpurm6xgnpVXUjqVkqlKKRnU1GGD4L96LsR+g5w51oqibU5Z2JeMvzndEe1xomk9bvk7Zx3w3EEH5AuJ6G2hSyo5Rkj8wIx1KblfPo0wilTEud116A6oBuDeh8teW3baLJZoFET4nbVyJSfEb3xvl2ik61q/DYA4IPafeQ7P9ShlUhsMyWR7klXDnV0LwztNzV2q0dLC4zOw45XRz5daRyL4dEjmvdv9O+0AwKD3NNSW7XmZVmlJt6svaveuETGILHc9RZQRI3RDF2KkKsYo5eQVpxO5VKRi5OK4vv+4t/pFy1JpjjwbcdKSR/J7RIkbTwoVuXoifKWluDYLOYAABxFSURBVFw2CqFfAcUVhdDvzveqxRs8/60dMiRCpAK6XzK0anxkJTo0NcOcFKf3P39O5vNvB0RvXNX4dQDA4soPAwDePFp15GFaxIHzz7YXVNR2qM/tGa8ioZunrgSZDOyYUYl9DFeSs8kUGCKvq1ZdIyT29l9sRnO9nGfiTUoCulr9pNRJvc6CLpExNl+36rCytk+u7kJSRE2lKKRnU1GsFcuSOinftsT3GD0Nka5lSA07IW9jizbDGzJbkCxzu/jv8pOro3QIWrFeFMsJkfa4Rn9oKbe3RkhBvyXn6kKVimS8Nq+WvTlzFjgs0ZufHhCUeU6R9NrKjwEAHr9NzlvxHi1uMTai1QEe0hOx/ehdjikg2padRJhTGsIcNszppKwi6Py6rsdAa+JaXQXBn0lxi22rhmDa5ept+fDt8qnJe52l6lh6VuJ9CnLeCqcSbN6O7iRF1FSKQnpYR4WgRSfx+sQ2leDXmsqwU6zp1r2CHL3VoC+ZKtGUkL7SpnOx8eWtNiRQkyjMBEGSIegz7aqDnI9kvv7LT7V+/+t/N+IX9YJuyxslEe/6yr8CANw2ThCt4ibRHcMZ0+T0sb6hcvEmQvBOOhW6Y/JnEseylTkyJY4o4wjSOR2f5Xb0ksw2KQyR0/Txw5tl7Gmzj6H8Kp0ruQYu5VuRlN1bGBnzEgRNW1ukr3J1WCdIGh7S1JcuJEXUVIpCehZRDQQJCE4u2S3yDTodTD9zG8X3eE4LGfebJf654B2SzJYfw76q8dRk6pOkoiXpZCTwumgTu7O0Rr7IAn9lYvxou2nVgmZH1H+6Q1hXX98JPKNI6vyk/6BINkLpfoO0CBnR0uMGxLcVpIC4sjtK1NbIEYYMTlyTLYv5VHktXt1YRqiaHhMkXb1NOADs5Vr2+XciZPoK/branaUgcZCRsEov/aYJEboqkp5+XO7ZgFuiLjCdSYqoqRSFXBwdtbOCA2GI4HcSu89vFSRtrhE9s99MQYjMdHn7QrJyqBt6bCsABX0AHIIokvq9SJ1Om8tFehSn7qVe06PA1OJ1fyaRm2/uFH33mcb7cFPlxwEAt4xRhtMY4bDa/hpj1+JpjqfKc8b7unodA53Oyh5VtPZHa9dtxu9jnoMCngRXh1pB4ex/SfpIQ53MZ9Fs2d7nSiU4N7fEond6R0zvxNzd34n4zisg+9tevRCsFu4qddLKpcrFuEx7onYhF2fp96ptRG6qErT/rgYAcGKnPFTDr9cw3nv+GACQr/BogBw6buxkksu2o62RbMGHjOQILlvtXogVHbDxyXpnFupeWS6/s1t+5EeOSNvGO4d/Gg9/Wx3a0zSvyiMIW7/SM9/fuFHHB5LfWb3koLzI4Viv5ipdcrGgh6srQCK03sP2pyU3a91KeXEX3qzX9Gdvk/2oMvTpHdWY8o04LvV8gRke9SoVBrU1aF8h9MayGyU7IpwlmbTGr1vbgaRLfypFIRdn6ffE/Ku82S3rzuDsUUGdEa9T6pk2NWPL8kRlacSWTRdSLEmGIIEoQzOfDKE6qh4nEm+my6Uzo+fxCNr/8jZBh/93QJbmZxrvBwC8f6TUgvpQ9RnYy7Qi4MiRyfNTOKbncnLpG7k87GBZBYL6+sTcExVZEEM2Glk8R4zwEX7nVwCAhjWygjQ0ijG34AZRCTLXSPXokG4romJpWVSb38vOdfVaD+r89L6FY4UAk1m/wW0vu16rhCuSOnJ1GkJN5VKRHtZRDVBSEul9iihtG0VX3LVjKC67Stwt4XtvTezjah35bbcpLpRXXohQvlHlf6dBFg+P+jVMj6n7abP4yX56QNBgeaNk4NBw+t739bgZ8xxh228T5EK7XkU8d65YXVdzSnXio1IDwY5SV45XRdCFen1idb9+btz63wtyrWoQJL2xWlaF4M/fKGMTfWl0xsKiriKNGo+u3wDJJ2VJonSgqTB2pyYUXrMgul5WkBmsDenYELgLSRE1laKQHq49JWK+LcXaW3eJvrl3l+hh02YdQ+lMeduZOm2DZDjUSWASX+M1P/3Ev05RmGMaUvU8+iHgej+tv0Oqrjy0R3QyPyz6jnHqrmES4KBBUbDAr1XgBtfznFXU9Kr7mWPHpMcTIiRlW0xHfvb7E3i1V4PnXnJ1YE81ybHLxgplcsCCqE5W4tNLvcbp0xHKMgGQFVqYxqMNhjOr1+oY+vtUax+wESOi1p50tXG1jKfLdCIpoqZSFNKzqSiBge1Vjn1Pydu56YRYwzfNEl9k2T+8D2Fn6Of5HH0Csav5Gdd3fCvbaxLrEMxrGmvLypzvMdgmucx05D9yRJD0AxoW/e6D6uiepfVcmfobBJEuSrT3PR68hrLktQX7auS4fQ3AdPXB+unRRCW/ygh9paoXtzyxH+s2yn1edL3o2cFn3ivXzzQbBhUYDqUlv0P6KZjTja6fF/x0Hr027pv7zTYAQMm1U2W3GfJpjh9zSEp/NsURuLuQFFFTKQrpWUStO47s//5nNJwTi/GaieJ767dQaxIlOqYkQ5XOcu6EbucofB0gsrOmvQa7BaHcWEToX64SnfQXdTLek5rSfOdw8ZN+cLLoVeEiTcQbOjwxb2TbHdoU1PAnKrpmuerPPay+yGNiWZvJoxB6Vr6bM61+EsoZsXpeakA1Py0Wfd3eAZi/UHTS0uuEThiycwtTQah/qs/WhVpJ3Rs72umgLnmxQYtYvCShY5uT85cslJUlrJZeAu736tsvqtPlJXea42m6dCqXiPQoorY0l2LT+uFYvERbiX/5LgCKpBC9kwlgTo/sTK/ztsdToa2ncxbU9KfE0rTj50QY4mHtvbqy8UEAwBu1OMTDj+guMxbIMYqkjpQRm7f1deROisCRdmhqtQzPKK26XT0pQlK/8rWfAqKf2XVCzVuzSVatuVMbUPY57T7Ie+cn1THFmr2hzrIfgaCiaWmJkivpg14nuvuR5XJNQ+eod+L6ZXKM31Wmdx+3YjmdlPPZtLvD+xKXFFFTKQrpUUQtK8lj7OBGlE0WfTL0LffYvn60piD1g2+nj1JhGKtG7R3jd1TxUmKO3fYTAMDy/aMdki6pvBsAcMtYQZLwigWJ0zn9kynI8bqk/A91UoqL/KgOy/kMVtpfRazPamvXenWwR4jadqv4So9skDFnjtHa+Uv7w/IYohxXDr9HF6+FSBvrnE1d2GyREmPZfRJNHDxVh5yqq4DfsZq2RZwRpucP9tfIn5q8+9OBpIiaSlFIjyJqaVUlhn//zQjjPZCAhG/UEZX9EjnUq7RrMXmONki+ayaunwbJfRxKu32SKSCfWilRpYcP34d/vVxKP9668qrkHHmeds+70EE6tTuv382ZeqxGvVhqh5a7jXFujccD9aNGrf8pVvemtaIrz5oj96vXF0UvtdlspIvS9+qlY9PKp3/TkaEZMSorA7RH14lnBHUHXaM+0Vu1o6JX9tJPWUdbW8Sl3SRzzj0vumnJTTPQnfSwwz+A7d27oNZpfEkucCXxAWGuuSNutBccm/ge3+a1y+H2TTf8AgDww/3ipH/4sJCe3zvi07h+Ql1iX7cs8kHkckqDzNEM9Xs2G70YDOmapMrBl87Px0pUNIyz/QEEKyVEmf29uLJYXn7OIjGEes3XEHRM3XDzcLlZXnh6cJJQzjz/QKtu2y370L5NHtpB1yndb6a+XJ09oB4omNOnENTqPW2S8Usu1wougwagO0mX/lSKQi4KKcUhqeeKQhA4soNrOcOlnsgWc2UBKFQVYv8vWC4pavD8pFaQ9BsHPg8A+Og4aexw/56rCgjaBQ2GKRzbr08fCy64VBBWIiSxhNeiq0S88bATLv06fnaNIGlBo+G/fgeAWFM3utzipGRXH0uR229+4RtC62RpPvGbFgxarHN8xx/Jp1/h2q/QSOOSYeHaOoRrpG6AuWa2HDtxonz697QDSRE1laKQi0Pz83TU+HfXjpEuHC/L0RkCvhEVc+q7oAERSlH62Tf9HgDwg30y5k+PSdvGe6oESW8fr9Q1G0aurdDTSSn83klai+3bP6LCcT5a2yqqI6VZr0Qy6qNMRSkpgXl8OQCg+WUZq26vrAKXXSOJeCWLtKMMj3HzjkKXrhaB3g93X0jKZgLeqnUyxEZB7ZZauabBrx8AM1XrY/H6iKR+zQAKv7s6/ZUIFmhSn9YecAmDmRRRU7lE5OLqqJQ4UZl6G9HW/+4I1R2/Y/EAgNN3tYrIvx8QlP7pMbHubxsqDXC/vFMbBPPNzucK33IvzOonuUUnDdy5/Yosxm9iRrcV00aIzjFHfNtquYat25Lup+CedwFIWvcJKYtcTE7fH6UeEz+EqZ6K/Hq5T0fXyPdhV2hO/i03wHrVVZxO6rf4ZM8uFdb+z182M5oHEwX/G11RzgtRO2jae7kxZqU27N1kjLn1fMZJJZULlfNFVL9pbzOAO621u7Vz31pjzNPW2i7LshmriNeZ7zOepsyDSGfzw60+2ZmIF0PCx6+RjniPiDsQj538OwBRSvM7x5Ful/TnmlyucI7uIjxPBcWvUdXWCjtKK0izZikR1fheCVYTEcTLaOWS1t8dRMM+CavOnq+Nhhepn5S6qEeLdEKrvHefgnpV7lK0iAW2yg3KNcpYI25S1Jyl+m97e1SRhejP2lasTq2IG2yUEKsdIp4NFsgw7e2Rj1V1Y9IKM+vXozvpFlE7adq7y1q7W/9/CMBRAEO7PVsqqVygnA+idti0l380xiyCxCL3djeQNZ2ESeM7+cjg+1opREGOwcrUtbWuT+d3dwtCPN34NQDAWwb9DQDgoSfV6tS33UsbTNINSYqJE6Jj4nRibmcyXmVlRBSm1c2xvKIZUShTUqLbVwrS1W2rxISlWi3wk3fqIR5lsSSZkOj8pt59SfyNnpSnpKJ0drOct+xqRdBl0nfVklJYXh5F10grZK1VFw5WfXyvovQg1Ttp2Z85U6CTO3/u4VeHOH07pFkvEDXtlRMZMxLAjwDcZW3HVXCNMR80xqwxxqw5fqz7bMNUUulIukTUrpr2AugHaYP+t9balzsbI95desG8KTauo/qFv0zs/51GlfwSPtTVtDjXo3ccxG+PiP5EJH1d5T0AgDsn6pyIpLQ2vcIQiXH5t6xHaCGi+2hJH2W8Zj3LEdFj4dEPzb8+DgA4s1J0yIMHZV7Vi06hZGGVXC6tbj89mvNhV2kiKOfX1ORQLbNB9cedQq7OHRC0Lpur1Q011TtR2ggQfdSrIO24B6vXJbdrQbhQO7E4/3IuFxWFYz3a/ZLUicpYPdhOpDtEZdPe8dbaKmvtWAD7IZ2kfw7gEWtt2lE6lf9x6U5H7axp7w8BjAEw2BjzXt3+Xmvthq4G605H7ehvTvxiBSzSoElmZr+gxH01Ddh2SgjQLqX5p4Io+Su0UzOtT78jdBzNu4uW+BEg6n3Uu+K+WJv0/TrrW/c9+Vu5lk11gkLzp4qel/n0HQhdZ2gP5ejl8OqmJs4PwA6MfJThhhoAQPNW8XZULBJvRHjjlbIDUdvF8aNCEcb15pIVI7N2nY65X6bFqFO1FEKLR9cA0dmpI2delBq49og4icJbbkJ30uWv8Wo37U0llQuVixOZ6iSq1JXQB2kOK4KqZb/pXrE2/71W/HbbTv0EMwaKvXfdcE0lVp9egXTUz8n/m++n9PtMMebtIVly8slibcET0us+t0XmvrdByuRUD5FrqphCTmusc1+B7q4d8lq9WHtbEsmCbdthdkuaSlud7Nt7qvpJx6j3g/P07ALHQSgtdeWPMhukcyIa5fcIpqqveKSmi5d4/NyYBKxC3Vf9ur1Frw/9roMdSM8/qDGnvs+aRxBE//faoFPxDp8R7aJ5uyxTH1kn7tsXG4Wq96mJn8V9X5WbGC7V+qRsfOAvwV7VlYQK4gcUuLSzIgsb3fq1AuLUOR6jDPqgQTI5m38jS/vBfTIvx8r/v1rpmeXNW1ujB9DPQeL9ISWQ7jGGX0kwX7cTLWtliS1fqnlNb3uDfHpZCgU1obi9dx9kNguhx66X4ICZp0v9dGluTAPMnZeDUiXZsxemUbNb52t7ypLkb9yVpKSUVIpCeh5RwzBy0+imxHeHFOpQVie42S5K+0OPSArEo3WCJC82PgAAuKryowCAd44/Cav56AzR+ehckGfV2fIa39eryOL+HqsiWCB6TLBcjIfcTnFsnzkiS++YKWIYls8Xgyf0m1+YoEDlcCqGyxD1mua2ivGVWSsrj82FKF+sqSaTRievl8IVxaW96N/VqZ/ZvMUt9S6sSvcTxbWr92oGcOzKfi671qlLXpOPriRF1FSKQi6uMdVVQh6NhN9KtPb07wQpfqL5YS82iuPh6sq/BAA8/lZR1PvdNQPhSCUomw7G72I+Tq+KZZC6MSgMADR7ehXJKqTXnTwBoy6a7GbRUVsa5AyDZwgKlrxd3GWWTcxoiMUNNr9yDA2cMBaqRWT0BVukmp6t1ezd6lEIF85LXGcUyvXqonr1ZDPaBh0HjgITBEHz8+cl5+oI2v2Tt8OrdxBWTYzuHVc21beD/fvRnaSImkpRyMVJRfF01ITwDVVk+s43xe3xkzrR54ik1ElvHStvcp/pWlt+yCCgvBMiCSmDnoPfTyc2ba2Rde3NK9Jn2eZbcdiFMtUFtWYz8tvFldYuajb6ztbxB6t1Ha+IEhNHAEGETJ3p0azwbDcKKoVtiqxLtBHw4IGFZG/PhebOsVfGcDWw9Brt1PHAgMrEsU6fJcFktybuqdswnCqOf1flMNsenZ/5/ZqOfT6SImoqRSGvOR2V3UeC/eKk/kG9oNKaxu8BiJD0lzeLBT3gTTrmdPHN2VEjI8vcr5Xq+03972yP3tgI6zUD7qyKoBN1tBM18lsbcHar7MtGw+Zt18r56BTXREaXbEfdzWs5KRtZ4blX4phQ/ZpN6wQVK2bK3/MLF0Tz5TheCrMjm/P8W5RAvU9Wp5IF2o9gSnVhK0nq4lxp9kvKiT2kCY3VwgByuioQ+U21QgyR2142Bd1JiqipFIW85nTU9bevBgD86pCUeVlzWpB0wYAPAABuHyMoVHm1vtkLkt3xbElJLB3Z0+t8ndQvD+SqW/cu8L1SXIovP9VCz7y4Qk6xRVaCpl159K0W/TUzVfS0UDuKuPNqgpwjxxTUmSqLup8Q9dWvHGi/q9a94g3pPUFDpuPUJxsnTLPAhZ/Oc0JWJVYEDM/IMZlxGhkbon2ySkqiqB5r9msfKbA9elbv+Qg9hlEuprD0r4wqWdP324/p22nnvlQuEbkosX6fBBH3WX5xi7xlj52U2P29k6U4xL0Pi+4VThN9xqXvqoVcUEfVO2dH53M8AkWpePkgP4rlOlL7vUA5ppKRWYF52NJS4GN3yJR8Eo5fxdCbu2mP6XWsNM3I1Faxrlt/Kw7lXvPFBxve9mb5ZP/SDojnjl6oxOWgTuYc/l6yiIJlUlUvvEw/Y94C66XVGO0WYw+JlW+v1QrcumpEvlnVz880Itghq0A4Tdqeh1OmJq67K0kRNZWikIuro+rbvvqGpwAAzxzp71KamYh38ygl107Q7iMs+eNb4b5eGv9bJ5WtHZIpWylR8rLAb5ocK/jFrwEA2a2Cxse2CXIMnqpRnZlVkZ7oz8OPQPE7rf14gTX6NnfWJs7Xa6aytyZpT1KiJZFfkdj2KgdCJYqTN1GrFroSl4OpEskLx2gXbHaXJuL37Ydgq9Q0JQMKrepBqNbu3yNGJI5lmo8ra9TeDjtKU16GDUdcgtpadCcpoqZSFPKaKOlz/3ax+n524vP4XLXopJ/aqpWeyVrizi7pzhvU57YChWktnZSfdPzPeHcW71g/etWyQrgFOzYLx3PKNO15/4X3ydD5XGSpF5T/8VK925Neh0Ssf4PodedWiI+3zxzlm96mHaGpC7pjTOK7aWuN/r+7RuamRdAyc8RXHF6lKTp+ceJYWjMRPazXHlivV52ULCpH2NbEPXo0YmWNXGkjRVlX2ufAIXQnF+VBXX61VAJZfpQP6JcBAG8b/H9w03ClfnXWTKIT8SmDADrNdvVVgYIQYwcVW8zDjwEAzq0WY6L2gDYanikPaPkcbRQRd3V1utR7D4TfQGKFkJTz6+vQvFeO5QMaVOsy7TVGcyFc77tpbITZq0vrWVmWM9PUaBw7Mnl+GmJaIj3zsuT948QZhIfO6Pm1UsswrTfiZawG29XI6idBFxJubP9KR93kHAPN1rDjPMpgB5Iu/akUhfQoopowhGlqwpd3yJLyXOPfAwA+OUGW+/t2LC04ptNwayffu1r6O1UFOqrG4jXaPfmcLG3P7xPi8BsWipFT9nVZ6q1vOMXH95b6gk+GMlWtyK0SBKxb1Qfjlum47/ljGb6bpd5lkLLm6doNCDeKK8v8kdL9tLlFAbKzjpQiX26lXGP2aB6931Al57/6iuSxrO7Hev/bamT75eKCYoML09bqanBlNm6UfR3qqkurC0kRNZWikB5F1IbtLfjS4o14rvHrAIBrKz8GALh+WIvbx1Wf9o0Y/r2b7/Ft3emoXVW+Dr73KACgfY8YBdsOij43e7A4uMunaKDBq6jiJJ5e4RqvdayjBpu18e1WQdITW2W/kVPPIjNZznu+OindQxlWMDnWiKBKK+sNZmWYpE7qaqnWSviXBBOuDeUz+8KOIV1Pw5+k6u3YqefRQo5jNEw8KGoHDwBobHSkH1up+vyAAYmxupIUUVMpCulRRD3Udgif3f15fJ4uqFVzAMQIDCh0XXVUaa/D7zGdtauOKYljO6mxalpacPhpQY66E2LlLlogLpTSf/gz2Znh19DXGRkODmNzQuJvjnyiOmn4e3FBNa4WdBw4V3P2//LdCOlF6EQn9Sv0uQomSo4JFk5BOENI1AU6qR/oWCfoGJ4VfTzzVk2VGTcm6tXFSikrtNrJZlkFzBRF/isWJq6RrjfTcNjps+E0nY/eu8zLq9CdpIiaSlGIsbb7VNVX7WTGHAPQfbzsf06GADh+Ec9/IfKHNOfx1toOC0L36IN6scUYs8Zau6D7PV87ks5ZJF36UykKSR/UVIpC/tAe1Icu9gQuQNI54w9MR02leOUPDVFTKVK5ZB/UYuw22MGcxxtj1hljNui8P3Sx59iR+POObe9vjKk3xnzrFZ/jUl36jTGPAhgFYLm19l5jzBQANt5tEMD07roN9qR0MOcyyG/Upo3ptgBYqk3oXjPizzu2/QFIo7yT1tqPvJJzXJKIWozdBjuZc7u1ll15e+E1+Ht1NG/dPh/AcAC/fjXO85q78FdJXLdBACf0pjn573Qb7EHpcM7GmLHGmE0A6gB85bWGpuhg3saYAMD9AO55tU5yqT6or6jb4EWSDudsra2z1s4GMBnAe4wxwzs5/mJJR/P+CwBPWmvrX62TXHI6qnYbrAdwDEK+yujneEi3wecAfOm11Mitqznb2A9kjPk+5AF4Tcy9i3m/AGmaFwLoC1m9/sla+zcXfDJr7SX1D8AHAXzX2/Y8gKsBPAvgYxd7jv/NOffW7wMB7AIw62LP93zmHfv+XgDfeqXnuhSX/tsh7S/jwm6DVwN4r7p7NhhjLu/x2XUsnc35uwBWGWM2Qh6Ar1lrN/f05LqQzuZ9ewf7viK55Jb+VC5NuRQRNZVLUNIHNZWikPRBTaUoJH1QUykKSR/UVIpC0gc1laKQ9EFNpSgkfVBTKQr5/0t2XBrnqyh3AAAAAElFTkSuQmCC\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": [
    "## Shortest path distance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x185a00da0>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAACZCAYAAAA1tzjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO192XLj2JbdBjhPoigOIimSoqTMqlt1u213ux3hV/+E/9Ef4PBPuMN9r7u7ppREUSQ1UJTEeQb8sNY+TCArdS8jOuSMaJwXplIicHCAs7CHtde2XNeVYATjPYf9/3sCwfj3N4KHLhjvPoKHLhjvPoKHLhjvPoKHLhjvPsL7/HE+HncrsQMREVlv8bwuHXyuHAufW3jDK9fB72UhIiIbdykiIq6zwcEsfC9kRSRkRUVEJCYxERGJWDhWNITPiI1j6g6J8uewjXNEwlv8Psy/i+EvXc5FcBhxV/h5vQpxriEz7zU/lw7nr5+yxu/dOQ7FuYUljnNzCWO2zblanjnG+GmJy2vCnO2QI6EYfxfhlcV4O6JYD1lgzWSJNXOWnBPnr/dg5ei94DWYe7DlNeAeaJwiJnGJWjhGjGvMj8/mje9GI/gMRTnXOL4nvF6Jc64h/n8MP//jP/765LpuUX5n7PXQNdJp+R//+b+LiMjdJCUiIlfThIiItKY46c0Ek7ydz/B7+19FRGSwvBIRkcVqwDnipmViVTkKN0VE5MI5ExGRIi+knsIxa0ncqDgXpJ7EIubj+DwpDUVEJFnAzYmd4eF1hnhgrDgWaH6N3z92MpjzMCvdOc7Vnev8cY72FMe+te5wvZt/wbxtHLtkfxARkbJzLCIizVRSRESqSZyrkcQ6XKRxnJCFa6hnx7jug4VkLvA34WOsoXXOe3RaxedP17iOa6zZ/BOup3eTxZx4D25mWMvLCW7npxGu83aJc93YP4uIiCuYQ9P5UZpxgMdZBtd9EMF1nyZxjg+ZCa6nMMLvGysREYn+gHNbaayDfNfAsXOH+Dw9xTWF/tuNfGUEr9dgvPvYC+lcR2S8BDKMN/jqaI3ndrLGThmtsHuHFnbIZNsXEZHl+gXH4G6Lh7FjEqGcpB38Ox3GMVNhmz/jmCm+kmIh/hzBTs4m8PpJHHkRzj4viIiI8yeg1LKN378+AI3uJ2kREbmbR+WOCHcHYJa7GXZ6x3oQEZFH55NnDZI2dnTazXBuEfx/2OLcXc7d4c84XogonUgAMeLZjYRymK/VxHzlrIbPmy7m33nl/HGM4SMQ7WGK6+jN8f37Ba7hYY5z9pa4mLtQW0REBvPfRESklvgHzNWKSYqv9N188ZnkWiejOGcshbULZfD3BuEyQGc3g7VUhLNuvgpwZgRIF4x3H3sh3XoTkjZ32f0CO/yWCNGZYke0HSBb14Ut9zL9RUREouGciIjkE9gROasuIiJFpyTFMGyT0wymU+BmOk3imDXacFEatydHQNH0IZAu1sBcFr/h5/j2UUREJgSpzRq/f6X9+Uy0fl6H5BnAIy9L7PCBM8Xf2j0REXFc7Pha+D9i3g6uoxLFOhQTQJkmNrzUEphzI4WFqRdgb9pEumwTJwwXIztE/hWoajuYgwzx3e0jrmc+wLoMOf+nz+YvIvJirgHrM7BhAw43QMxKAnM/d5qcW1yaWHI5I5IdhPFZpy1XrmKNkxd0+M6P8IWPuG+SxgHciwsREbEuL/E5xvffGgHSBePdx15It9ra0pphl90jgiA3Y+yQ1vZJRETazp9EROR5+quI7BAuFSuJiMiBXRERkdwWuzxrJ+Qwhh2bBSDJIT2pQ7XdYtjxGdpwB0UgX/wcXwh9hNcXt4Gyiyug0+Mj7K8pke6WKH3Da7ieiPSI0N01PL1r688iIjJadEREpBj/A+ZCuzNHr/sgiv16GLU41y0/ce5Dzjl5ABjS8Ei4hHPbH4oiH+H56c53PmH+2wG+M7zC7en1Dzh/oIveg655y+DclxsgfHf7TziXneTcS5wbvpeNWnIY1TVe83r4qW+VQ4af8kBXOYIN6x7hfqoNZ/8Ge1E+0ZYb88F4Y+z10G1dS3iPZII5ynSLf0wtvEYWG3yGbNycYvIHXIzgwstbhBjyUfw+HwtJIY4bV0tg8YoxHPMkjVfdcR4PRCSG38fKuE36cCUieJ0sWvh5MaSTw9fQysFDPdvge3oN07Ursy2OObFxjuUGr5VS/EfMgWGcegyLnmUMsMx7UYo5nPvSM+diEcdLnOLaLDpHdh0bwWkNxI6Gzb9FdiGe7RQ3fDHHZpnTaZtuQp75z/U6NhteA9Y+7OJhO7FxDQ07LyIiVYagTpIiJ3RodL66oXMneGii5ziGcXKaJyLyOw7DBN/Xh82dMLb4xgher8F497EX0i0dS67H2IV3c8DwJ/uSP+O1unWwg8ppNV7/KCK711IhhVPmiW75mEgxCrQ5T+N9UUhi15RLQJ1Ug1Cf9Bq1VvRZRETmv2B3PXVhzQ8ZSvg0xutII/adOT41AHw9nUnXxmv0YfMT5hNB0LfpnIuISDWBeWugWl9LlTjmfBzH9Z4eYK6FEgzpFL4u4T8QKULc3xcIi9iRnjg/34vIzgFShJ5Ncc7WM17ptzPA6s0Mc/Dfg5bF0MgGAfhz+7+IiEjDxjqpg3aWxvfOkitpHgCJqxWgYywLxI5+YKD6e5hB7gWQzT0D4tu/wmySFpwUecZxNr8CrbdjOkNvjADpgvHuYz9HwnGlPQMa9RhSuF/8s4iIrLf4/+PU34qISN35Hp9xoM0hDe8SbaE8EeMoujVooSmibA7HSitaNLHj11fYlXJFhLtmmOMeCHc7hN31uoItdzsPc974Woe2UnsGhLgJXUt/+TPnAdf/dEsbLomJalqrzlTcURT203F8yevFXEsVrw0X/gCDWy5gC8kvQCO5Zijmsm9s0Kcu1mhMhJ6sNBzlRbgOzafbOd4EHetWRETu10jR1SJ/LyIizRAQrsq3Sh3mmZwmuc7p6e4twjW2cziXzVSce1bnJ9bDur7m/HFOtwXH0RlgLrMWfr2cMhf7xgiQLhjvPvZGulYINtxghcjreAaX+SAFZEsJdnjaYkJf0y0Mh3DzSYppokx4a1JFmiKKpsnEyHDXZHCs0AG2unqtr31s4e4ICKfI8LzCOdtEBmW+dIhwLbslIkDpTBS2y5HLzwjQJhdjKMQgMhAup59EugMm9CO0l3ZzJqQzTaQ/a1hk2V7Lcwfz77wiJDIkwo3ppSrC6XV0pjjntY2113tQieDtcm6BKFBPY5FPiHANItwpPdVKcSjpM6xx+AODvjnOk6m4L4K+inCfEJZZXQHhV69Yp4c7XMOM1/DWCJAuGO8+9kK6hcylO/1HERGZL+H1KcJVYthtZw68nXoKiNFgukURo8QY3BGDkfnEQo4y2DXZOnZkpILdYuexVTc/AR3UhuvfA9k6RLjWVG0f7CFNC10xdrQkr6wdgnennnYp/qPUHOzoWhTHOiXVp5HEfGuJNa8HKJFNANlyBcw5USGnr8a4VoPIUUM8Un7CObe/ACHmlzjeoJeWGyKczn/IOOJ4rfPHHNSGa9k41oRpvkYECfymVRYRkQui7EVaKWBYiCZTWydlEAjS55951d+TmnQIu9l4qT/D1pVLIJxzxcD7L5jLSxdzHpFW1eK9mG3/Mo4FSBeMdx97Id3GWciKCJeI4d1fiH4UEZH6FjumSoSr0Os7YQI8yzRRiUiRT2DHZLNzSRWxtU0UvAwE2F4i9rPqkS41wO/7TGc9LnGuhyXO9YhDG3qSxq+WFlDpcQ4SQioKFCq6dSmEAMVHZMTmSTYoMPuRZzorl8Z8MwdMwRVpE5Wx078gYBqEgJenduigt/O0b4kSt4wfjohwo5UX4dr2Da8T89c1P3CwThkyjpVelSGTOs3PVJS28gFt5WxMJIs1dA+Y3vrCS+1w/kC41SfMRe3QnkYKlhopwOdkQwryGyNAumC8+9gL6Rx3I9EIdnQ6Bjvia15fzuf1ZWnD5UgxP8jgM5lbS6RAjv+vQKTIGDtz1cN3x30cezDGLntYkMxI8qIiXJ/JyNEG3x9b9LRY3xAlherY/k5ERCrWoVSSsB81l1plpqFKRD6mx3dUwtyiR0zc1xjXqsIWMvSkFebgPiFup8n7xQhLPV7gWobriLyuvQj3yhqIyRqINBTYYiPnjmuIt8u5w/nHMYcasyXNFOZ+RlSuZjCHYg3HiX/EutkfSiIfmGnwe6kjfEde8B3nBUi/HNLeZCxxRC91xLzwcI37Pg2QLhjf4tgL6UQsKaX+RkREygK74oLId5rW3ebznJjjU8pMtoBPMp0kVIiJfQK0kJ9h/yi9XCP19yPYQV3G4W5nXgKp0quu6dU920CGhzFoStlkU0RETsL/SUR28ayzTFQqRDhluDSVfJlDxD5Xxs9Ko1Kqj1Wjl1rBhfjpSet7XP+kiyXuPIBdcjPVYpqoib/p/IdbUpqIcC0XFKUUWSINh3lhXxFQ3VcEdHoELzVXAeLFP5LG/x096g+n4nzE/fuCmkQC6fpnZH0m1zjH/SPsxyt6qVrQ9LoCsrV4LdP1Vv7SCJAuGO8+9kI62wpJxiIZ0wH6aI0qyz1FOQaui//oskxuShtgw2h7htmB5GQpkSV2lebvxi9Ak+4rdlVvrl4ejtExEXovgVQp8sM5PMdUDDv7IARkyztAjFwCx8lFRY6imlOlzUlbLk0vNVrg9RXoWecZuS8wt3rDXGoL17C8AVq9djHnx1f8/RWLgW6I0u2pSJukxGsH9uDQAkrOnAHnDbu5uWUsMYE5NNOYU5XofMq3SPMQuel8QzlxtOGIcO538FDdi4sdwv3awvxpkyqnb/SJBUt9Ml1MySlQUwmkQ3raLZZszty/zKfb66EL2REpsM6zEMEFHbGm1IQaot5Qg6amtnwItSpKh2XPxSI/f/KKY75MsLjPTNwPVpjm8xLHGLAW4HmNcwxs0GxGC3yut3g9HSdBqyq4SLoXwpjLjlblSjG28cz36IBB32M8jCGtSa3g9Sh1OE7SBS3J7ZPa84SHTesZnun06OtUTQKzYSYbuWGlXMf5v/juGg+uspX1YWskcaxaCmutdQ1VEi9rdBgKJ7ju2Hc+AuYFA8CfOw18nTpMay0ucazVCOfQh+1msjMHRHYpuR7Zyq9rXLemFmcuqv7eGsHrNRjvPvZCuogkDG27zCqoCpPKdQaB6yRgnhxi9ylV/GGBnRKyYl8c12LwVhFuwJDIYAl0eCLCPTOkMFjiFdA3VU94na432OnFNBCu6YAqrwhXS3upPrXEWmpJnLvM+R6cMJDa9KFFhcFfBk3lBVveeSQJ4YHhjyHnTpPgaRn2XIOi9NNmLgMSSNdbHOMkAWpSg6nE87SGRPDdBqvjzhnGKWfwqQiX+IGBaiVgKj2JCGcImJdtQyBVAmy/g9f/lG+XT3Te2kQ4ddqumaNrb7H2rzaQUlOL+TgclLdGgHTBePexF9LFrLA0GBrRUEMl7q3zrJKIeVTFz84WuzTBwKwGFR+JBFvHElWgfaT912dqxZ/eemTw9552Q1+QspksECIpphHOqQkQ7jSBcxbiGlrAcZophnPSEykfeREu9h2+Y10wvMCCFLlESs3tIRzhvGBSGsAePTM1RxtOU3RPvIb+wuG1AFke7J68rnDMahTU/lMHwd+TJNbmlA5DXR2GFN4i9SM4DIdlpuQuWGGmDsNHOgxaRMPAr/y2s+OUIq8IdzOEDTeho9fy2XCdCcNSdHrUaXM3mJvW1laYDm3L/5KvjQDpgvHuYy+kC1uWlLAJpcj0VoleX0mpP3m661UWskxmnmPEWJr4MNdAb9LYfR2mWPpL/Kxuudamtrfw7noC1YCnKeg3edKrFOEapGtrekgViUoxpZoDIQoHU0mXvHQqq0ovtcrotY+evWYUdE3m/KiP69HwTodo3Z6FeX34u1sSMG+sXSFQMQYvtU6EqxHhTlJesoSqVFVJUTqsEuGaDFg3Gagm2cBfJmj5rmF1tZCXO8yzN9Z5Y+3HpFd9be0V4TYO/kAD7vktyxxjfJ28MQKkC8a7j/1supDIxzR2bIXUpOMsdp8W6cbPaF9QpyMThncT7eDvskNShJ5xnPEiZlBv6iuGnm2AUFrQPbNwjMUGNt1xGsTRpoOi4lMKNtbSOM7HNOyNDGlV9aRXXyR3vpRIw+elfgRKGBuuTXrVFRBO01pTBkmfGI9ThOsygN0zSIG531rwFh82QOdy+I/ywcW5Tkm+rJHwqnooF/ROTw6Rkjuqc43/CDvMagBdNA7nfAcigP2vQCMTuGbgd/4TEPKpm5bWK2y4ywlJmFMtSMeaX5J0ceXiGHcOin/yYVTypC18/yKEN0Iuie8fJ4hjuO2/OwKkC8a7j/2QznYM9bl4iM+DCuNamnI5ho2w/Qm7a/Irds7LE3bniNSeR6LbcB2WB8ayrifw1p7o6d3MsdM7IXipjywXVBKjSYAncOxG2psA/0AafIIZjyrVng6bQNnIh4xYp4y/+ejlmrhftXF9z7e0RUleHDGepfHH1hTXoLZQe4rvXbq04bawQxvhv8OnVZIzIpwqKNWpgnnCWGezCERXGr8huf7A4pk6bTh/IbQqeDI1N/8N379v403QHmXkigXdVxMg1C0VVMeUqPjNhT14t0a2pBxB7PPCbYqISCrCAm7S+3OsRyozmvHWCJAuGO8+9kK6qO1ImRq06Qqe6EiDCGe01vAyX5Ji/joACj2RYq6F0BqLe1nbJpb1tMB3+gvszEf7nn+DXedHOC3k1pyk6vw2ksxJ0haKKZG0hv+PnmssrizuKeJwFtUv3XvYe4aa9Agb7ZkSFRpfVCVS9bRNLJHX8ODgOE8ukLMaAlI0SJhopKJS02wOEe6UEQDNjhiE+8D5+guh67DlDMWcdqhqFCtFXuU2OvRUb2ZxuTUF3JjvzRJrNbLx6Ue4potzNajEmeCTU2MB01FEya+BgE4wvsGxX5wu5EiqSITTMsEK7AS3DftBbaDhHRDhnghxR4QYrXFKzTaM1pY8Uiu3y0KUB+rMqad3EIbtorJdtRiOqfEsVV9X+asKsyNKMbepDRc98cXiygWDcNIm+bIFW3Xcwd/2n0kgJVIreXG08RcDYaf3VkCpHoWFVHH0xC1xbvh+JSlSjStbxBsJ0EyDvkWMl0qGi0G42zY/4WGqp61yFa/33nugGsW9uS29GdbsboFzd22gpdKq8lG8TWoO5l+hkJDm2rVovhrHuQoxrL3KbLw1AqQLxruP/ZAu5kjsb2AXKF3b+VfkPReXeNL7t0QG0povJ/BSH2n7DElv1lzk62or/RV2m/acGK7h8alcQt2Ft1ZPYbepDacU87MUkKFBanyxDMRI/ghUseK4TIMYGov77WYXpW9hh762SFIcAME7M2xt5cP1TLkg0FOLgVQF884COfJIgBAXFuZ+doC5aPF5PbGWc8bhaowbml4NynD5gUKEPjEbE4e7xdqr1MPiF1zD/TXJr2Pci2sWc2ufifbEkes5bLcrG/G34QprfhxDzPN7l/YzbTiNIZ6Ty5dieeM51zzLvHC2ROh/Y+zHHI7ZYuVYL6kV39cMKfS8uiJtuuRqsD5QFXREDr0mvl/dmfRDeMUttlgIfdiaLha7mtQL9z5sxmHQMI4+bOdk+35gOISKl3ICJq6mtpxPj7JqYZHGPTwUvWc8bEq+7M01nYVjKnlRQwt3DggAmh46svFAN1yYBLU0A7+GAqZOw8yEcMzD9oFBXwaqv66cRKfHBK7xsD3f4uHqGG0X3IO2cRqw0TvzmbSph7Jxcf2q0KC0qprvYTtNKunAG4Y6LlCtKs+2Co1AyyQY3+DYrxosFhany5BCBzvktYfdpBXfBuGIEB3alY9zGrcboNM9K7bGVl/mTGud2gicllzUH6geSjVJio8/JJImLZ3acIpw4XM6CtoMhP2qlIC5Ux5ayLDHBiDPDCdMvK/TLhG6M1FHAed8pb6v0quKbNtUcxFkbqSxDg2DFKqJwprU3EgOG74QzhlJBqrv60M46zOEFhFZt3EPhrcUvqY2irZt0tZTO10+3IxW6FKmJGE2LSTss6x5qScSvnlrCtFLzo0TsTNN9g4jE8Suce3fGAHSBePdx35IN1/L4hciHN3xmwGe7BZDCtdMHCvCXU+wQ3psezQmvfl51cIE7KgchmFHHLk4Vo5dAFUt4MivEUc6VTaDY8fyLKIpaqUW62h/wTm0U5/qcix/w/f6N6kvamqvp7tqLZHfoVVZSGdt2dSkZAHhzi2EM6ppfL9hDG/8nToNJ3kgZLa+kugf6JR9VHo5C2hULl/TWj5tOP896Lwwec/QyDVTcmrD3cxwvZdsTDffvhiEu4jAudIWBXVS49VhUBtOEa5whp9D/LvIH1gVVwTKmsKlN0aAdMF497F3Q7rZCzUsaDdoWuuFuhwsm5QhC1BeBDtkKEhpTTdAm2QYIZecXZeSA6/yhPZEMWHxZ5efbJXkK6LJVmkT+YtolGL+ZyCFIty6hR2vNal3w7T0SDzofEFJwk7vbGFvKnHUYp/UE4YUylQerVJqVNUv1cOuqVK8r/An0kwYsoF7ygS+v0cDFcz9YR0lYD4MvSjdoQ3XJcJ1WVusqk/9OZCuFP9Rsi6+m6ZS6gGbsGTofBrlJ9YDJ0nxD2fZD0M7zWRV0JhrroVLb4wA6YLx7mM/1aalyP0T3t09X8ujDuNYbXp57RXiV7cW1NenK+zWo2hTRERygjhWwT2SYgwepKa1ijE/wnnR4vCE1fdNX+tKXxxu+wke2vqBWiEdBn4Zi7udJaUz91GSJl7VAEW4LNUuDx2cqxLGMUoJrw2nqk9aRFNTTZQG53xOtacPJXG1CNrvpd7Qy76EHawE0lGba/3Ce2Ao/95CblWPv1Vq/BIB4EQEb5esFCUbwrEOaTcfKHBRC1pVtrSTTjyLdbFz9FKbsAWVAmaz07Wi8lsjQLpgvPvYr/Xm2jYyA7vWm1QyJ0LcUBqhayMdZLm0gVhmp/ZbPgykzMfCpkRQ9dVU6sFQfVgmuKMmeXsefEExv1KNYpYHsnim96Iavzsvz1+AospPty4QWtuma4bhMKZ6drguLVTaST0wJZfz6orEPtDu/MBY3HnjS204nf+Nql+yDWibiXoi9DVjib2FNxZ6MyHFnGQD7ZFxGAWi5plOPLWLpouOaisfRtjuPemd/1GNqlUXPuUnFgEZtSpS4rcq+PzGCJAuGO8+9uwjEZIWO6LcL/C8trT1Jm24Gwvkv62D/9d4UI7tx0tx7Jg80S0XFTmOY5edMtNQZIlglYTRg9rO4xOhkqTILib0m1cMZnWF7z92cM6HiYrYqAq70uNFeqSVd13YIopwdQuF2ycCG65Oarm22ixzzsdx1bVjjwbanVpEY8RsfqfXll8bzmQa6GUPrr0IfUM7usX53/mo8Z8Exxms2V8iinzqsUP7+TN5DbVBtfA8yxac2rUo38DBNZdqbDif4qjzG2OHV94OlG+NAOmC8e5j7y6IbXqp9zP1UoFGtxZlBlzaBkSKJuNx2ShtoIRKiynN2TH9XRtaTEy5Lm09bsoE1YbTnCoJmIoQy0sg3EuH+d8REEKLZzQf/HkHmlvBzr1zYf+YHqkWzqXxN5WkOKSeXZnkRS3c1hbv6lnHLkjAJCr7uwla19ciV23P/FdXQLhRV71UasORmqT6fDr/uxnWp+3i+5+XN4qInNKGKyX4dmGv2lpq1ydMu+hk2AEoX6Xc2BltOCKcnNGGWxPh6LVu+liHyRPmPJ8HvcGC8Q2OPW06y8Th7sm2aNvokzrbICamMgMVwQ45ZgVHlraQ2m95iiceRTdyTLp2/oi9UkmJVwVz04XmhJ7TpVcmQRFO+WSGbUHGi1LjTRdBenkt6cqjQ/vHphChhXPU094Mg0qhaW+wMudcYL8yE4drECE+81JFficWd3nzGduFXWjapJdTouLWdAJiHM7EEnGuLtXje1vE4aphJb0iQqAsnSLt50JM9aDXJjJQKcJLjR/guoy2snIRtT36JWJ+GvtcdnD/pgMyXAZA5en6Lz9S+z10W9eERO7JqZ+skd5S4uUZg7411gI0mIPXHvfHDIcckVN/lJgb9cvMKR+2slaY+RqCaMW63qyW92HToOntTAPXWGwlkHZocF8JFvDV7RpHp2YjndX0JexPEl5K0iHJBmaDHJO8qGpPPvXLLwSl+Up1Pz2aFgaDtnf+qhZw5SMf6PyvLTgMLw4233EIWi4fLKSilHSgJoE2Yil8FoqqHcPxOzgnNYlSqlrVJ99jk8gY1ymXTMlRAWBLvT3HwedWP92gTVMwvsGxH9K5jnQExupoBRJmLQri5blgZ9dJz941tsXuOqRLXuJr6YgNNtIHC0lqU7ePgGij76uvU20Z1FJ6NsmLJGDq61QRThWTlLz4QALpDUkHSj44c//WNKKrmgYgLufNanvV3aOjkDwA2sSZ3w6fkGL+Pavtterer37pC+ssLxcyaHO+dBhupt7E/TVBRl+nbQtIPxWQEJoWAu7H1BU5Y7a+alpuepW18gxcH2Vmkj7hW+UMa2eVSE1SpYNfSI3n63XzCecc/8Yalz7+/pU6fC0GrCebAOmC8Q2O/dJgspL5FraAFnI0WRepVUO7aievQqfaQkc5GCjJPJvQFWwJnZDM+IHB3jJtOU0LdbDL1m3suvE9dURegDL3rKm9Y1rogQVJinA9F3bogMHTU5fhnPjBrnGeQWaGb0iFNyk4VtuH80ChUI3IUIWTY0IiWpOq6pdadc+g6vIa6/B8mzAIdz3RkIhXl++GBFjVtHulgtKFC4SrRjXgznahX7FDVTvwIIuFSeQ3EimzUq6AY8hFnfNVSjxCIg41hqcQKjBNTLRVw+sa6KoB96BNUzC+ybFf6013KdkwW25u4ZYXE9gxR6TIFOkplbTym8TLDHXttM1mpETEKCZ29oQ2BOmwR0MHCLVhrGB2j+88vaquCBBCyQf3c2/gWhHuwaKeiEvFzijOV0nan5EuvTacX4vYpODy+gWiMds0fbXqXgmYbZZcKoH0NWOasvQWfuVRzKVDiv+zA4Q+sX7k/DOcv5d04KeCqf7zwRGVO9kbI1yIGlPwPGQAAB6sSURBVPVOtZNtlmkapYYW5rtWJa1HtmciwqlqqjbV07nPNtq+5usjQLpgvPvYr2OOFZELokU5rgrgvraPrLbXLi6lGj6jjO+Gq0QM6thJIbvT92UhjSonrT55dUWUQKotg7SxbWuqBFLYkbdLnFM76dQczLkRPfTM+TztGi0ORTittjcpuKZXMUnR2K0A6d0q00Mah1P1JxIwl7/Cnhq0vowlXrHivvVZWk5ETGpuQI24hvUfRGSXUtSW70rn17fLGVNalZyP7FrFeewSbedi5qvUpM29NzIwolLnLRWfrqda3ohvvrJNk9KqVDX1rREgXTDefeyHdBL5rJkwntdkCE96jD2/QhbLAW18KoEytcROSAjsjTDf/fZqI+Lg39qjYXPnteGeX7w23AM76TwsLP4/zv24xC69D8HbK2+b+AyRWh732kCl2Nqks4osEcyU6aVqCk4VnkiFd4+IdEQ4qwfbTe4ostujHdrF8aYP7JgzSnHOuwbJu/4YQAftj/EoSM1peWPFwhyKTCnq/JVWVWR2p0Dl0YOiF+FCVb5VjlmaWciZ+To9b/H8ok+twCHm+0LlVCVNPDIT8UClrdEKn6rHN7HZpPiNESBdMN597Il0ITnQkrWIFnQ4/IQ9dRDFrkuRHNijHbZaYSe57D+QWMF+CS+3EnKAVOqlLuC8yvOA1HiSF+/mXgKp5lTv5zhWzwbqLF3qIVtJzon9JEg6yDIPnI1sJE2U0BK7UJbiO1nCSQ4o4Wb56UM4646T7SCute0gc8GGjDJ41liiz9Ne2HJPxae7Db7zYMPLDrn4m0MHyJSJ4TaZNWd5YCbsX3NqKRPQtIhGchrAY+HS5a1p977pYBH1rTJk29MHZnde2IHyjh72Pddcuxe9bPl2Mfp2eFu9NQKkC8a7j/00hy3bRPCVrl1lbKjK2NBxDiiTqdLOWMG2UZsmMsYutcPUprU3YkXwN6sBEG82gh0xXGgsKMJP9ZiEnyzodlmix5xqTUBTqtDDLlOAp5LgnOmxVtJTKeUx33SNuchTxgyb7LOluVR/HG5I2+UJdpgzAEpvBlvPNWhuUhFD24+/rlx5XVFUyEYsz3Wp2yvwtqtRoGSVpZnVhFdxVD3uI9qjmk81Hrfq8WkHR9WzG0zEGeB+rZ4xn/EIazVk7FDn++ybtxbRP2+p50chofEG3u9fMwKkC8a7jz3V1S0j16VoUVYJKQrDmH6pLBPMrfH7CO2oCTvNzO5Y6DKdS2qO76gK+EDJi1ONflO2S/tskUjaob6vUuXrLqnmjMdpXwlF5Top2kqLPy6OJd1glL7JGNYFvdSvKJhrpkFecG71uP1Ucy0X3FHNfb3CJivDGhm6QOiGMOMQgodcYzyxntyRLz+ff5GcvnSJeWyfqKLwGgzDhaWZ24eFrJ7o8fe8xUsDIrMSRxXhbn0q7B2WOY63mHvIxve0ZHMg/1u+NvZrMmyL1AjtGmoo+ZuYNL0V4PElq4dCbG05xCvgdQBD9eE5I/kNXg8qyNxnWEFVMO9ovPao0aEacXrhZeqKnIRhQWt6SF9HWs+gJkCJJkC64Ui46W159NXXqT5sLe11DwNaDXF92FTnbieM7d0wmuLquk/y6CI0osFffdiqVJE+YXvTKskTquWiD1umxpSiJu/1daoOg08LRQkT8wdbJiM6NhNcf880Rcba6Zq/LLnmC7KV6TD0F6C45eIgex4qeZf6fKip+/0RvF6D8e5j7+YlCu2FQwZTq97XqaGYn4HlaDMIHI3B4A6/MFQSJ9qMFvL6ClQYUgFquKb7vga0j+g4qM6vGq9HLqhQzRBeJw22S1ed3CbraLUGo5LHq0FTXOGLrFhNbw3tF1X3mri/xGtkc4PX6naI15PSrLRGoGs0XrxU89sJUKnl7pr0moo53/x3JEyt2PISSQ+aRDhF6TJjJKTIa6spffUvP+H7isbPo6S8ENmUOKqkg8FS54u1fl6zZShbZakI+UEMxI8GW9bnbSx6hdVz8kbkJEC6YLz72LP1piOVY6BMqo6dHq752o+XacT+Blto+xt1OTrYnctXhj2edw3qNL3VmmmTYU3g4xzdGRPhbF+ZctmeiTbQ2YFXl0MN7ibl5rWIRgt/DC3+u+pOG4423BcFNFQwX3/CsSYdnGtGpYM+7dDWxOv0KMK1iHBXLsIVdw5sobr1N3LGtpVGV8S03sR3NIFvlA5Uhf17puYafKuYNYfDsP2FSgF3OM7zjVcpYLCMyTPb2l9NvUHfZ7aZutrADnymNvQrW2XVw38vIiJpB7brWQxrqWlRbW7y1giQLhjvPvZs07SVVIMIV/cFIBs+XRGSA5dsGTR8wG6bKJrN1I6LSH+pdZ3YA89LL8JpfWfMxTGqIex09fIqPi9P2x4VGPhNVbV5ns6ZfRpOa194qSaASoTbsG2TItwjG+xN2FDv3qfkufNSGdahkqcqCGh97YkUdl6qn0iqenwkkhq7uanGnip4spme0vq1DrjNqvtHzpF1tKpn11+G5XmFtVY9lLsZm9dtga6qrTxb4U2luih1wRstGWGEgIHrIxb2V6hA+tYIkC4Y7z72S/hHXQmfeds0GeUk1RXRMkEWQqs+rqbB1EN9NOmhkLHh7kmX6S8YyyLChV1Ms24DVU+Sfn3fDX+ml0rlpHSDCFfzaRITIdx6YxeHU/XLFs5pEK7NxnNEuDvabkNWsj8uvQh3x0ZvGrhWxNCgqSJFNRX9Yv6K0OUDtp43mspx7/xpw5l2odSzW7eZmmJvD11zbTWlaDxYWqIycl3qoWjnn76NwpzxEohfjlO9agvbVxvqpUg+0GvQ1puVRKBPF4xvcOzXGyweEusjI/aMelu/snW3FhGzR4MqJ7Wf4d3cM8ugMbjHhaqxW6Y53dUU9sQje01kSO3RTENN43C+YmLVhlPlJKMrol5ejXYnEc4j9aBo0WbU/hd458MWEZlljm02d+uRmjRce1tv3pjUHOagGi+HgnOeC5vqpXeN6VS90z9/o/yk2naqi6JpLV8cTjWJtZneLddcCaOqKagpuMHClZeVetWIQz64KAqfMoNUS/wD5uY0OV8c6ywDhNPWm2eMJapMSJmlm2+NAOmC8e5j795gWgitfapU33f1CQj3pXISc3xLnGq09lLMR6utPJDgeRtqiYhIYQs78SwGpNJcqkpVaByuzlxkncU0ptcW24/LR9puNWYbtBPNZwqYu7wk0eLam2HoGfVO/P+dtt5kQYq2T79ew5ZVGy4vONep4I1wqghHQaF6YmMQzj9/LXe0PjKHqoVLGhmg+uWqy8gAVeNvn7R4Gwv1uPQV/jDL8LSZy8DGfLubf8KxmP+uJFDIfe5DOO2k0+TbJUkV9qZqCqrA0HHQZDgY3+DYD+mi0V2fKl+vUbXhjJiNj9Kz6/eKHaLFKENnLo82Mg11FtIUo9TGTSlbBN/1I1yVPRp2qutexU7Ta0sp5tqJRu24q77pQjPpktEy0PlTwXyuMUTsdG2bPiKB9J49au/YP1UZL1U54jWw36tPbqOenH8xf1PQ7Vcc9WVHtF36+MGruq4CQrrmam8qwml/sye7K6MtbLmwje8cxtmndnvO+WpDZ1x3jYydBktMtd9rmdmSZA5zijX+jfXpxHVFqCDkMphIVJYVGabzDdNEW4DofGvxZ7KCSQCYOdTIsOaSdXGDyjEsesFXtVVmr/sypVZLFCLULsrRGqk9NSoEaNW9v56hgzCAdJja6sxk2sN8tZbh4XdqGURE7hkKuWdz5JHLLtJ2S0REjl3crKqFazlmU5NjbpjjmApjs9d9emqqtiLHrEutkLW8YNiBFXRCUWl3QelVhpiWvjUf00lTE+aJr/6nDXb8I+f6umobh+E4DVpVTPDwxdiGKhrCMdinTqKs7ovaOGaKKlxhhkqiZTYvWQQV/sH4BseerTcXsv4Tdsi8w/biPRj7bbbwviIzWJuaqOT/3QrIMCQtKWbh78pWTgpsAqzuuLZpOmN1U4Wv00qe2hya+NbmvDS4/WLS9p/RXkBrPNXp+bzq/s6XIrpieGEnJu1tytKzQHmK2UDG2hZUqGYcKFVmDckxUVrDOqpeVdOmIPXF10MiJF/Kn+HwbK4QGpldYV1uO0D0O0OjwlpeTnDuT+wKqK9TbZXVH0MmNhJKSy3zX0VE5MJBWu4gpI4OEPoi7dXpUxWrGgkf2rbJyMV+xDU4V0GbpmB8g2O/1pvLrUxbQKPBE3Z6h0FT3W2KcJ2pV/J/aCPxXXaQBspHAQW5WEgKVHzS1t0q8V8jteeYCKf0bINwxmFgPYNPTFqu6TBQ3844Pe2d06NtK5Ua3/GFF7ob7OyOhYT9AUMg1S3sxVJMDW7YQmXqi1Ti6jBok14GrmvamC6+00fxNdLTMM62z5pUIlz/nmvOVFx3rkFfNnbWZoBEuI4gQP00wWeUTfQOYieGAKuKDaoekGPi/oitCw5ZU5vVxnQ5IlwTCKdNS+Ihkg2ov/fWCJAuGO8+9my96crdA9NaPnuiTYTTSi1VTtLRZDqolMRWyseBbrnoTnXolO54QenlRaCDSdwrPVsJo5q491PMr71qkqsOdqu/9WZ7umu9qYUo/tab9xbsKq00ywlQppxk9RMbIp+mvEVAZer7KsKZtk1samJ/KO3o5Qz6+psJLwdUMriD3el/q3TMW4Vt05faSAY23GCKQHUmgXXKhvFGyLsnhiJfT/v0Uehla+vQExJhCyew6UxjOhZexRys8eISCDd/Dto0BeMbHHu23gxJi17qHeNYLdKdb1i5f72EDeQIdsxFDJ5WlnRm3VG7Nk1bE7s6zVIbLkfl9TOtSSW9vOlreeRPa2kQVVtv/kaEe9y/9eazhZjemUPyYgzXfUhdFBN/i+M6zonSqh6vSgca+P2i9eZJeZfWojacql+qNtxwjJMo2aBNhFM9vi7t5tsFDNGWjfrfwQTrUUgBnYsWbN2iA3Q6isQMwqmavCpAqRK7tt48LOJ64g1c9/wX2nZrrPHsFxbwPGKOQyp0vjUCpAvGu4/9kG4bMuqXGqlXL1VtOEU4rbKvscr+kJVpJVbbF6K7zjlqwxUq2LExmmyRc5bWKcJpEQ0R7vfSWiI7qs/gjrp2I2296e8zIdJlZqVH1fiNxbiaA9Kl6omod3pI78603qSnfUJP25Rmlr3F534CqVx3TMW9enyqR6z08pel2su04eZeL7WzxLnaTMFNN0DpQgqxN9VEKQveNnlmSfIx+wvyhDbWK6o+SpWtN2vUhj5HpiW2YiqOXuvklXp79Kj/mjZNAdIF493H3g3pbk2/LXpMc+y2uYXdehaB3aDIoGWB2Yg395iPYWfl0nNjwyXOcGy7TP7POTMNfoRTXRGNa1HfVxFu1CY9iYihBNJbX2O3znQttw4LiCzM59xC/E014bTwpO5rO660bLXhyiUmvstAoUjDl7yvayyO3X8+Pe4a0XW8CKf0clVMMgjn81JV4uFpBRuuEQXxMuli/RqhI84Rx1EF/KOYaxBOG9PlqeKZzdOGY6ZhOwLyyxVif9MWfhy+eN8ibUYzZkHHnGB8i2NvpOtyt91RUGVswZY7DWFHK8LVTdEJdpRGtkssntFGwsncSmLH2IGhc9LKy9T19Ssn0YbzF9GsW14V9rsBW1fOlDHiVUwyVB+nb2w4bVd5QirSgSmpE14Pr4PsirLaPuwAZNSftPhcGS/awZEIp30a1q25iRv2fN0PVXTnmYpJinC3M6y5CgeNNmDPlGKw4eoOsgypEHuEpbwxuEJM255uvixzLJIeVsX9s89wL+wxzrn+RG1iEkPH9FJXDv5+STbRwgk65gTjGxz7ea+OSHeBHdK14DF+tGBn+csCT3wdaLL0jnJH1Bw+ou1zHNrZcHV6qaTEfyHXRYQTUyYI22bSJlGUCHdnSu68ZEYlYGo+dWNv5MQtcd7YuSpek43AhvPLdGXZdvwozx5ntOGMvp2WZmoHR+1Ry9JMIy3WiRpZMT/Cddn554VF5z0jq4HrV4Q7DgPhSltFaUBaMuwtD9QsQ56e9lF0ZcocDypEuBrOHSqlOF8gskOW0LRNVfWBt3POK73VLu3OyV/RG2yvh27pOPJkYTJNVxvRebn/pz52r4olJrNklrJMNnzMYGkxLVJmU7ePIEK6JTwIXygnXbHBxhUfthtWqbMNZNukibz0pIe5kg8QFhmTXnVh1czrVDU4lHSgYthVcx1/ofUmxRSlpmLSbERnulfjyR9RDPJucGBEH9ucb5epOA36vqxwLm1Id7cGVasWgZ6ImgTaKkvrGLRSS+tpjfOmGz89k1wFJzOy/xd4nTrX3pDIinQpbUSna9xh6OyVhFFd6+kmqPAPxjc49mtIJ1spudgRxTie9AKN1CKpMEWmUQp8raYLbHd0hB0RrnyGcCIipZypMDMI9wh0sO4RPJV7iknf45VAprWRn38yGh0Rfnrp2v0VdvVTCKktU4uRinzWSE/njx2epeOTZ+I7nedrKO+7DtWGU+Ukztm9Y7CZbY+mfQpGj/C9/iIufRrlOt++Khyw3cEzRd60IV02gtDRsYP1UoQrxG3OHa/jVFivhWLjDOvkqHOXPlwZevn6DmhoZ0b8mUHfJxx7TkR7mv3+GpMvKk+k0k+oIfjWCJAuGO8+9m5ecpbGE6+hhHNtRMf26Joo1jRK/HtSebRlZZVhkRI1fgt5ccuwTUzrcUWLT0AmdRhmaszeU/WSxuzVVFUlvUHUqxXsz6FNO9SBzXiawuSbaUtqVBmqkqR4lqXtxuBvtsK2lWfY8XaRZV1sLmxsOG2m90jEoJ6d1tFqlVnPVGxF5WaKPX+jhNElkPyWTUy0YuuQlCSt1DpL4RgnDFxrKOSc9P40G8lU2RTwUJ23IsM6xbBJa8nPoFMtfqEOClNxSv+abthEmGusenY9Fiq9UuJfiRIT2stvjQDpgvHuYz91dcuSUvz3bSC1G7ShrdoM9gltnkPacLR93CN6rOWyWPdsdWRsOLrrxoaDrTJ8wQ5XG06bpPWXbKLGEMMTS/c2QrqOAzQ6JrW8aGwgx8xf01mHh9jxd33MO0z7KETkC6+14R5DAzZL79gudKvN9O6otUcb7sFndz4sLHlaeG04rf993cCGS4Vh4xa3cJW1HriU8CLccYxeKq/hiorpWi5ocaquy3TjZiNWjO1N71kS+sKA9Hhnc4qIjDchzpFrvKCywWd1yyIifRve+nwD+/utESBdMN597Nlk2JKDiDeBf0gPKZ2k0lCOSKA2nCKc2nCfIZyIiHV/L9YTYkPyiF3i9oFw6we2PHoBor1MsNOfqYL5zOJi1Vp7oX3x5FL9yEUgM2ljFx+wclgDv7noVrJR7NgDBn3jB0CNzJh0I/ZbiA3olTn4tGIITGkbci2iWT3h+idEuGcinCbvn6n29LISGSypLcLY54sDAsPGwQVlBF7qIRvrHUa1GSDnz3uQpQ2X4TWkeE5tRBJiW1QdrrOQUJpNgQdYy9EYyDbQfhKc75QJfFVdH3CNBw6u/4UtpkYrxBIXm8CmC8Y3OPbTHLZFjllEo1HuEqkxRiKhSokHbc5LxSH3ED9/jnAiItbzs8gjkE4RbkPhk/kzkExjW4OldxcOGCsa0M4YUP1paWNbnrOHezJMqjyLaEqfUbOLzDiohxerknS6wM+DPrzV8QtQw9ng2FaU/TDC8BDXj0TlAW0jHyoPVtqnQee8lcGGJZFsGbqmfXQUbYqISJkJ/FKCZYK0RfUeaBxO7VFNMR6zkZ/aZYPFlxTy6CNVqojkL/xbVV3X+ZLF/8UaP4Vw/0YOIgwTKndafwWOBUgXjHcfe9p0u5zkQcTbTJjsdLFyDODlYMsZhPOL2Rg7bmAyDhtWyWjG4eWZOmszr7eqKp6mHTkFZ+5tfLHAiL02F6a83Wf2KPOr0ZWk2VUndsj27px/LMtMBJPtGmdb0460bMa+2B5+8kC7c8g5m9bvJCPonJl1eFwujbf6vGqJiMhhFASHjODtcBiOcd60RaMuf97w73Hdao/GsryuFyI+aUct6tVtXaDs1rEk/MRsjRJGl6oFHeGaYs2mvoyDrvFgg1jidIn4nG3hOmMRULrWpM7/3giQLhjvPvaM07m77oeU60ofM35Vo3rkCcmLVTbp9SPcK5tG0Y6T14lp0Lt+ppzYkD3CaGe8mp5h7LbD3Tdkh74XtmjPsgClEsXuLdOGS4Rw3HLMKzlWOJhK9ljlxigEQ/JllJ5lesv85drbLDlKGneItPvxxNtU2N8YeUjlzqE2RpaxabmZiiB2WabcWE7gdZd8cmOlmCqYs4iGudTDvBbR4HYejrxqmColpm+KkOVKakYN6JVqKNN+Xvn6sTEuqW3TR9SDVoRbMS6XS6EIKE0y73j2k3xtBEgXjHcfe/cGaxSBVFpiF//B2xvM/Z6cOC2E1nxqTw01eKiqDL4dLExU/LENO/BJJSt8PVMN3Xzq5cdFeRmNOBBCKfOqXp4MsckwGSNGrut08UV/WvnIHmH0eOMxIHI+Qv5ZD+eeUlLs4ZbayiN8PrIIXcsclWfWntI2dHG8W/efd3JjDrsh0jDOxjj/z/SJRXYqmA1yFDUvHCPFPPQ9bME0W5tHe1hruwuUjbLwZ7CMyZVP2u1lpUqjWKvWFMeeuNSDtiFRcT/5k4iIbB1FW4gqnlgo7M46LGiX/ylfG/u9XkNbydZJVaphsrsuyl+pSdWKrS4Dvy+4eesejrN8sWVEh6HDhejzFXXra33kf9hsArVW3+/qM1h7SydBpUq/VE5K7B42U8vAJixXqiuCv9V6hgdK/I/4yt+6mEPbNNOjpouZM4mYgpBCn42FQ1ZYTigVu6sRxvyzrBH2P2yqK2Lm3yAJgalGnfP8Gtc/ZCNnrYrTmpGnVci84rQyzjQTpuPUtpDWmlq4b32+LmMRnCsexubUh61he503eSNGHLxeg/HuY7/gcNSR6He+9uPfQSdDlZOMrojWd1KnQ2sDCFIy6muDurhJ4F/Rtdfkcpt9MLQKX5u7JQSoogh3mtEaWyoOEeGah9huUaoJGIT4wJrU7453FfdtIJFjVAK89HIVk1aJ/xFrA5asfrqaYA7aTE9fp1cCxLjfQAUzGcK61d0fTcuCBhGulvSmt85Zi1rNAuFM1b0qeDaOOGfSk66BqoMH3KPHqSobeEXHB0sRh5mxm7G3mfBtCPO9W4Iav1xjDSMhnLMW+zsREUm5QLzzMN50NXLk05GgGiwY3+DYs02TLdZHVtbUVOfXqw0nn7xaa9q2SYOnaoDrLnxdRg3d52am9Bkcqk17qOPgWBkm8E8TbDJMG65JIqkiXE1bvlNTLUzepV/yX+qVXT2qasN9pereL/GvCuYz1nve0t7sM1DdsmDLajFNOoL10vbjjeiBnFI5qW4QGt/VoK86PIdVL5FUqf7G7iQqv3Sp8jSkU2NsY9b90ql5XjiydgF1/mbCT0u8qUaqbZeE3amN6ZoObDmtra2ba8CxDyJBYU4wvsGxn7p6PLILiZzAWzUhES25+xkBT1Ui6t94NS9Ga4YaFruArxZ5XI219aYXLerazM001MBeUar5mYYS6N0Vy0C4xA9UTGJNrnVG8gGVA6xPrV0zYbpx4x6R+BVocm80iSP8ZJumNZBC6VRXC6CSJsLVJtLk/amD4KmGdRrpkEHoOgmiOv8D0sTyF7RBm16Edv6FBAEW0ajd+TBk0zxW3/dJjNBGwvdUin9aL2RkYY1unP8jIiKTJeY9XwL5c+k/iohIJfK3nD/ugZYrUITd6NvVidKZcFCYE4xvcOzZpiliEM7qstq+gx2iiLHuAOHG90QMIpy2Hdd0y8Ny193lkbQZ7UbzaOFYJy48y0ociHVMmrZqw1XiXl2RfIEIV8Px7BpZCGlt6sCUXJdz7jzLhobO9JZli6/e+fZ8nXNUrUrbTanXd8duNMM11uUgAtu37MK7L1PnTvtMlOOO0Seucv4F6qJoy6NIjUh9nDHzFRFZ3bLbEMsa+0TlO+32s/QqG9yz2PxujfV5sh9k7CKNNZxj3qs1bLtEDPc3H6U+H1XkK8k4P3FMLejWa6gwPZqJBU2Gg/ENjj1turjYP6OfgiKc8xPsC229+XTjjYJfsUxQkU0Tyopuo9XGFKa8CDIGFRexnwYRSjU56rTh6rR51EstV1k2eIE9ZNSffgDKSIQGiCp2Us9u+fNERl1tJsz5+qQp/KrrvaW388+LhesfM/V0FME5T7bIzNTjKV6LN5bYSC7ljDZcidp2aZjLYmeYafgOqcWdNAV7fLVwzCcW0ag2XItzvqSwfX+OOauH2rWgSTxadGW5RsDUdTGfVKKJucSRYbhwmC3xNaZrUnbjkLHPujZILjGeF//LvcEs13X/4h+ZP7asvojc/NVfCMa/53Hqum7x936x10MXjGD8W4zApgvGu4/goQvGu4/goQvGu4/goQvGu4/goQvGu4/goQvGu4/goQvGu4/goQvGu4/goQvGu4//B6lMgKlHIAeHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kng = isomap.nbrs_.kneighbors_graph\n",
    "A = kng().toarray()\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "specshow(isomap.dist_matrix_, cmap='magma_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trumpet scalogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x183b566a0>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAACZCAYAAAA1tzjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfTklEQVR4nO2dW4wdx5nff19V97nM7cwMyRE5pEhR1IWmLpZs+bbZ7MZXrWAbG8fZReBsNkiARZKHPC2wGyBAsHnP2z4ZQYJkHwJsYicx4sV6s/LaVmzrYutqSStRoiiK5JAcDufCM2fOpbvry0N19zmHN81IZFu260cQp6dPd3VVnX9X9fd9VdWiqgQCVWJ+3hkI/OoRRBeonCC6QOUE0QUqJ4guUDnRTg7evbuld9yx91blJfBLxLPPHl9R1T3X+m5Horvjjr08/czXb06uAr/URPbTp673XeheA5UTRBeonCC6QOUE0QUqJ4guUDlBdIHKCaLbCc79vHPwS0EQ3U4wobpuBqEWr8VgMNwOrdtNJ4juWtRqw+3Qut10Qo0GKmdHsdefO90uNJu+++t2odWCNIXIF0NWL0GaodPTyPIyxDE6PQ353xrHMD2Nee11dGEPurAAgPnLx+HQImy0oRbDOxd8t3pgAVbWIY4gsj4P9ZrfjixsbKJ334EsXYB2x7eKtQjEoAcX0X2LPm/dLtJu+7z2+0h7E7IMpib9Z7eHu+sI0ulAMkAXboN6/dplH8W5d2+J223fchfpbWwgWYrOtDBvvQWX1iCK0Nv3+3wmCdpsgFOk00HrdWjUodlEmxNX5+E98IsluqLAtdqwC4yGRdD5XcPt228f+0EKgQG42w94wRbsvw1d2OPTr9eRZgPaHfTwITiwD7l4CRp1/wNMTkKng+7di6ytoYuLYCMvpjRFd83DIEF37x7mrdlEm03o96FeR7tdL0DnkCyFTgemp/0N8m5lH2U7Xf+VabZaFLNi3KFDyOysF/peP3pIOx1fRud8PsH/fRORnUzMeeSRezWMMrk2snIR3X3NkTy/2LTb/tOYHYkvsp9+VlUfueZ37zdPH4jKHulit81ggKyv+9ZrtNUD35U557vrZhO1EWZlBbe4eM2Kl5MnkVNnkNkl9MAiwI3rJE39NYpWaGPj6jzcJOTsGb9hI7TVQrpbaL3hv1u9BIMEcQ538CDS2YT2JpKXn3YH+gPoDcAItKbRI3dc3ZLvkPctup+74OC9Fb5WG+tyx8i7Mh0RmJubu25Sevgwevjw9q8dRePd3i0SHIDuPzD+90g3XZSv6Ou0Xof5XdzqSanBeg1Uzvtr6fp9zOvH0b23+VajsKbabW+JdbtIvw+DBNYv+yZacp0bQQ8fQut1ZHUVcc43/ysruGPHkLNnkBOnYN8CbvfuPLdXtBB5t2reeAOdm4X2Jnr4MHLqFHroEHL2DDo7hzn+Brp/0T+w12IQg1xYhvMrsHc37uGHkJMnoTUDa+vokSPI6Xd81zrTQk6fhtaM75amp71BsLUFo63fxgY0GmNWp3n+BW9BNxtIv49OT/k0895Bzp4pW6KiDLp7D7JyETmzhLv7Lt/NX97w56QpsrKCHjq0rZ9HTp2CZhPZ2ID1y7iHP+zzXrTg7bYvzxXdu6xc9N2/jZBz52FmGo0imJy8fu+wA95XSycdb/qXGSmsqclJiCw6Owsrq95aPL8Ka23Y7MB6G5LUW4Dr6wBoqwXJAHfkiE974zLsnvNinZxEksR/rl4aZqDoVjtbyNo6TE/5tPbv9/vjms/LpXW02URW15ATp7yl5px3hWz1yuR0fpd3D6QpevtBZPki5s03c5fLxaHgez3/3DNKqzUmODl50j/7TE4g3R7EMdLZgo3Lw3NG3SJJgk7m+d+9ByYnkH7Pi26r659BNzaGZdsOziGnTsPqOnR7Xkij1yzKU6vB2prP9+ol7xoxBtShsy20UQdj/O95EwjWa+CWcEutV7pdpLvlraDli+jiXkgz5PJlEMEdPox5/kXvfDTG+4RmWpiTb8PKGsy30F1zSGcLnZ7yLdH+vej8PLK66rujVgs5f977yvYfQJaWvF8sSVBrfbr1OnL6LLSmfQtzYBE5s+S706VlWMh9eOsbMNvyd36awlwL3b3Ll2N5BdLMH1eLYXICen3fKkaRz197E5IEJifQqSnfGuxbHPOZmRdeRPfs8uc59S1Ov+9bvZVVn+auOWhvgjG4e+7GnD3ru+G1DWg2vNUIaGsa2Wijcy1fR7Oz0O9BvZFbwFPo1DQkCdJue0OhVvPb8/NXGVmyvOz9gXmdXdPpfIt5T6KT1Uvwze/SfbHN6pkJRJRGM2HQj6jVn2N1bZLUGZwK3eQ1jCjt5CJWhsHzTA01k6FsYuU0ImBQBs5iZQ0R6KWW2DgaUUbmBAAjz1GzjshmbPTqiEA3jRAURaiZZTIVJqKTbAxqWFE6acREdJrECYkzzNdX6aQxiRNis85c/TidJMKpMHCG2dqArTSil1mm4pReZrGiNGxKpobUCZFRmlFKzWbMTW0xf6hHugVp1/DO0hy97DSDzJRlBZDcLvTnXmC9X2MyTonNO1we1GjYFCtK4gyJM0TGH+8UYuPrrmYdsXFM1AakmUVEMaLUaylbvRqRzUgzS+IMmRMaUYaI0k8j6lEKgDWOWpxRr6dEtQwbKy6F+qwj6wqN/YI0DHbvFEx79wr33uFv2POr/rEkzWCiDvv24B768K0Xnc7vQr76WZpf2OTAGyfRY/f4Z7LcW7/r+RdLg8E9eRzzwH6SJ/yMNDsXoakikzHmnttgsweZQrMGgxQiA90Ud3HT/1ATMTLtwzIl8/6hV89cQmYauLOXkJoFI8juKej2Ye8ceuIC0mqSvHyOaF8T10nQgSO6ax69tIr2Un9Oq4GuriPzE2i7j7lzN1zaRDsDiAy6lSCNCJmbQDtd9HI/z1cd6jX04x9Fbz9IbWMDc+4csy+9gbtwGbcxwM7Vydb6SCSYhUlIMmTXJCzswr14CnNgDhoN3GvnMXsm4cBueHsZzRyIQJJBbH2+9uTPypsp8siHoLPlw2j1OkxP+ha9NQVLF6HTR7cGiBE0U9K3Vokf2QfTTR/CSzNoNtAjd/hn66Ulb/wB7uBBqNe5VeNrwjNd4JZwo2e6m+enGwz8c9co/f5NSz7wy8N76l7Ny6/Q/69PsfJmk7PrM2wmMQD9zNCwLzBwhm7+PGMEMpXymcvpcF/ihFQhFiUy/plnKzNEAjXjSJxg/KNceV6RJkCSP+cVOAWHT88I9PPvMxWsKJkKmULNKAMnSJ5WLEqiUu6rG8Xl6UVC2c00jCNRn2YkEBslEkfdKA2bkangVGindizvV5ajKEPifL4c/u43Amn+fSy+vqz4HbFRBF8uK0qcl6f43uTPgjYvi8n3Sf5ZM45eZoc/vNHyXEGJjcMaxaA04xRrHI0opVZLqdUzojgjTSzqBGMdUc3XSmOPo/4bB3BffnTb+nlv1uvJM2S5m8qgHF1YoTmZkA4sW1s1LnWaeQULF3oN5msDTm81cQp1qyQOGtaxUE8Y5BXlqwfmc3Fupb6CrCixUVSFmnEoMBMnGFFW+nViUVYHERORIxZfkT1nmIoyNpKISZuxnkRMRhmd1JKplNvFg/1U5GinlgnrcArztZSeMwycoWZceexU5EgVUvX7J2yGFaVVGzBRS+gmcW7w1Gnn1zKS+8RRpmMvzKnIn7fSj5mKvLW8mVqmoozYKN3M0MsMsXGYvA4yFWLjShEtNLuoCg7B5OXoZZaGzdhKIwbO0s/8N5upZeBiFpveJ2kEGjZFBJpxQuYMg8wyGSckzjBZG5RiqzUzjFVqu5V0w9Fv29LwsLFiJw36sQoMCfflR2l8Gfbh/4+yC7j9vSQa+MATAY1r7N9prDbEXgOVc+tFNzrJZYfIysXxv0dDYFBOmin356EcctO/JB8TZt54wzuVi3SL48E7SUfp9338FbyB1G6P52ckbHRd0vT63+XI0tK183wl1ytbwTXqeTS/RVk+COy4ezWPf5/j/+ESf7W0hzNbcLaTkalXb6aKiNBNMxL1zyoOyNThUAyCFYNDUVUcSkaGxZbfJ6Q4lAh/nMWiOFIcDkcmPl2DkJKhONyI09nhMPm9lJJiMGSkCIYUb01bYhwZhuGDdSL9/PyMWOv+mpIQaVymH2t97LyYGlYjIix1reHyjqYrPTJJy/y4EY+X5k9hEREJA2Tkvjf5v+L4ohwmdy5bIl+HGhFd0V4UdQiQkZXbJn9WrhPjeNnXPQ6LJSHFIHl+DFZM/gxpiMUQG6FmDNZAzXgjDMAKWCPEAnN14e6plI/uWeWe73x1Wxra+TPd4gLTk2eZjzPasSVpWqYiYSLyBsLACesDi+aW5ObAERnhcuILGOdmnAjUjMnNBiF1iuQWXqZKppA5xRopBW1F6DtHbAxWoJt6MfdcRiymrGCHEothy6U0jGXLpcQYBup/jLpY0lzk/gcTtlxSCn3CxCTqSFxGLJaeG2CxWDFofmMZhLpYrBEa1lCzQub8tTvpBD2XDW8kzXzZxeCAhrEYhE2X+DRE6Lms3E7Ukajzx6ti8oqx+XVFYDq2Y5awFSFx3rodZEqmSpLXaVG301GUNwwQG3+8FW/sZA5qVnJjz4ssNlCzIAiTkdLLpHx+E7xFvafu+FCrzd6FkYEMN1t07tgxbvvmMf7x8ePw9llYXsf93lfGjjFPPu3jh90e7kevY+5bpP83pzB1wc5GIILZNQF37fOe8cjChTVf8tYEpA7WNsmW2ogVZLbhoxL1GNpdWJiF1hTuqeOYg/MMnloivnMKYguquAsd7B2zZGcuYw/MkLy6it1Tx60NwIDd0yS70MW0YrSfIXVLstQnPtiEfoo5sgc9v4FuJRAZBm/1sC3BTEW4zRQzFYFTzFwDIoMcmPex3Y02nF/DXdwkO9/DDRRTE5IN/1M17m2inQR77wJMN0m+d4r47hbULOnxNey+CWRhBrKM5PllooOTaC9F+xm6lWFaMTJVQzd6mC887KMRb52FuWlfh+0taNZheR291CG71MNMWJJzCb1LhtaX9qAbXXQr8ZGQLI8zL8yhr5/15Rik0OnB3nmYbML0pB/smabIydM+8tHu+GN2z6J337nj4U4hIhG4JVQTkfhFoJitf6XRcC3e7ZhrzfwfNUx2wjYMjnfl3QyRDxA3RXSyvIz51ncw3/g25oUXIU0xP3oK88Mf3/jE3EKUlYtDC/P7/284meSKY8tzTpzA/M0T/tonT/r9g8H4ecUsptEfoxh+1Gz6H3pjw6dRWKl5Wubx74O1ZdnAj/Itts1rr40PZXr6J5gfPTWcf5sfd70yyIkTQ2uy3UbOnvXpvPQzzMuvbGspC/P8C2V5ZeWiH+C6tIScP++t9OsxejOlKebZ55ATJ/zfoxbwFTeCnLruEsI7DnfuuHv9yb//HVb+bImn397HhV7MuZ7hUh96mX+AHYZ+8ucYK2ylykytGOYDqkov8/9VlcgMH5D7ma/wqdiOPRiLDM9N8zxH+QN2ZLwh4sNgSmwMSW5wjH6qegOmKHKxnajDipDp8MHfjhgmhYFj841iPwwfyK/EGsGOR+lI8+OK4wujQdUbVYm6cn/deKPDyvB6kSk+oZf6OjEyDBEWwUanPo3IFPmGuhU2E29gZA6s8fnwYUGhYb11asU7e4sQoRVfjkyhYYeGRWx8OHEqgqZVPjzX5uN/ehh3/30+jzdtEGeS0PvOKR5/4zA/WolY6WVcGnRIcxO/z4CB9BnIgIykdGdkJEjia0BxOPWuDiXD4bCZj936gI4/PsrqOBLAux0KV4NTv89IXG6P4nCYPO7rMu8+yTJ/XKbjd2/hntA8P+V+iVHnH7JFLE4TBIOILd0YIgbBErnh8G/N08pIIAPBXuUCGXefZGguNJxPs0Ayf66IRXJXiuC3rcSl28aowYkrPwF//aI+8rJZia/6O8vLZcRi8/g5gFUvC0tcunHqWi+/8zmxxGJp2Rr7JiIynebjr74FuehuxM5EF8c0vnCQ+1+9TM/NcrFvuZxMMVtTagYSB5cTf0c5vDHaSZXYwOVB7sjNW5fICHUrpM7fuanLg+95OoXJD5StUPEZG+966aVD18Boy1Lcvf3MUbeGTpqNtXZFOkULZhAGzrt2VKFhTZlubIStLMvv+GHTVbh/IjN0A2XqW5lu6sh0mPfCZ1Y3Pt26NahCL3NlvnuZo27MsPXNr13UWXFNyVufyVjKltt/R+lHG+St/iDzt3bxORNHpYsqNv742PiWUqFs0WojLSr432UqEvpOx1pvK8JcTbl/psfh1mX0N359WzLaucvki5/n2Bfh2E5PDPzSs90HtV8t6zXwgSCILlA5QXSBynn/UxD7fczjT+DeXEb7GaZVRz7zEXjtLT+97oufx/zVd3GPfnbsNPP8C3Bxzc8ompxA5+f8FMBT56Ae4x79NObtt3F33eX9UXENefEVH256+xy0u+jnPon87G9xj34W8+xz0PELv+iDx5AXXvGTVM4s4x77DOaVV6HZwB09ivnhj9HDB5G3z+Du/xDm6edgcQFOnIZmHfezM5j7F/1k73YHPXavnyW/soreexfyzhl0723IyVN+muL5FVhe9yG8hXm/ckG+4oGsrAzXwXvyaT+1Mc1gdQO95zBy4h1ob+G+8hjm8R/4qYdz0+jRu/1iNp0tdPcu5PJl3KFD3ofXbGK+8W2Yn4G5GV+hSxf9350uHNqPGoPmE9dHkZMnkbPn/DTM1XW4uA57d/l0uj3/OxSrMvQHfupoZJEfPIM+cr+fqG7t2NRG88Mfo8fuHVuq7UbszE937HZ96muf5aXvzvPkSotLfWGlByu9jF7m6LuMhrU0rCl9UjXjrceZmh3zI3VTR89l+cgTKa22nvNujZm4RuJcaVmKQJI7TQf5CJaaWBQlFkuiWTmipSaWQe5ry/LAeZKPdInF0teUWIYjTHrqA/oZPj99GWA1IiYqR8GAtx69BTk04SKkTHvo1xPvbBDBjdRv4ShKcrdFTeJrbsdE5eAAixDlI3NiKdw10M1SbP73qD9Rc1+lFSmPNwJ1a9hM/CCEwrLO8oEFNbHUjS09BKMjVES8PzRxSsMa6tYQGe99sOJ9d1bg6Izjdz5+gol/9xi6uHjz/HS60eXEE9P89blZfrqSsp72WZON0i+XmD5GLaR+iFCmCZo5UvpI35Q+Locb87EJtjR9vH8vI0obpDqy5EMuuMKvJmJRl5X+u1E/W5GWOoeIQZ0b84mp+v0iBlV3lf9u9JjRbe83G/fTGYnySw7TVrJhua5MN/+u9M+NUOSnSHt4nXEfocl9buhw+FNRHy73f5Z1oH5olU2v/qkzTVHNsMSIM2WaUt5kQz+hJUZSQ5w1MGqIqBMRMeWm2G0ngZjuRszka2/6hSJvwI5EJ9MN7vzkZR7YmEGkzuUkpptO0qp5f08vEwZO6WXe95MpdBI/6aadFN723D8kQjMaDpUphtlYIR8i5crhUXHuVysiD97bLnQzh47404oWKHNKzUrp++plbsz7X7QIxV1tEPouI8oFVUQHijt/4DI/6iz31RXnFteMZeinA+9zy/LxgqPRi8j4ctTzEFrf+VbYTyJy5f5Utbx2ca0i0uOHVUEzGrZi4CMIhZ+uiDoMnJ/zMXBeitN5l1j0LKPXKOq38IO6K647FQv9TEtfpfcXwmxNODKVcs/0ZXb/0YO4hx96Vx3t7JluahLzJ/+MR/8Etj/3J/CrwnYnZwfrNVA5QXSBygmiC1ROEF2gcoLoApUTRBeonCC6QOW8/5eXnDqFPP0S2VuraKpEd7TQz33Kr8odx7i778I89wLuU58YP+/8eeT4CbAWveOAX2W83cZ8/8dQi3ysNl89Xc6fR+fny2Vk5cVXoZegHzmGnF/GfeJjyOol5PwFOHkW95lfx/zkOdwnP4b57hO4xz6HefVVdHYWXVzEPP593Of+nl/dfGIS88xz6Ecf9ItC1+vwwnE4uAD7FqCzhXvow+VLWmR52S8AubWFObuELuxBXnoVNjpw2xxYizt6D0xMjC8qnaaYN9/E7duHefZF2Oigj9yHnDrjY6+Pfc7P++gP4NAi7tgxzBtv+AUK09SnVcQ7ncP8n//rY6b1GiSpX0p3724fF80yn48HHxj/sZzLVypY8YtvOwenluDAbbiHH8K8/Aru/vswL/3ML1V7uY0u7IZeH3ntLdyvPTK+onyO+V9/Ca1J3MMPXPP7qzSzk9jr0em9+sdH/iVPXhRe62ywaTZpyyqp9shISbWHYMdCTqMhL7miYXV5+KkI/xTnqTqMicbCVm4kzFQMHweuCoMVYaQrw0w+/eGxwjAM5q4Ig2n5d5Ffv5jXaAisKM9oqGysYmW8rGVZcDcMgRXbwki4beR6AMbEN7zetdKPTG3sNwHIXH8Y4pO4DC+OplmE/ozE5afBf1qJqcsUkzrLEbOPP7gr4/M//oK/3s2KvTbilI/Mr7OezNOMZummLQZukanYEBsYZNDL/PBo8KGVrSwrZ7P7CwoOqIulZopYrI4N706d0tfh0ggxhgzFImRoufRBXzMydWVQfnyG/zCwP9CkDOinuHIm/2iIakBKlC/NEBORkJJJitWIRBIf/NZh/LiYJ+Bn7ttyqQyAhLRcBsOM3GgGKdMvjiuWh0jzQD/k80jy/BWfMnJdg9Aww5/OqRKJIS2FraQjAwsG+fIaU9LI0yC/pk+7WPZjOEDC5kt7SBn+m7C2XDFgNDQ3EVn2Txpun3B86uj21kvZkeiiO3dx/+N/nwdOnhx/LdHamn/P1uwsstW56tVAO2IwGH/J73WQ8+fLN/ddk/fyvrDtpv1B4MryFS8iuQ6ytPSugfgbIcvL25jJ/5vbSuu9LXR95Xuw5ubQvC/XbfTpN2QbggPeXRTvUXDbSvuDwJXlu9FrO+F9CQ64KW/KKQjWa6BygugClRNEF6icILpA5QTRBSoniC5QOUF0gcoJogtUThBdoHKC6AKVE0QXqJwgukDlBNEFKufmiS5Nx18ynK9c/vPkeu8SuyHbWV5/Y2M8rXZ7WN403f4S/aN1NLq9nXyOHrfd4693fsXsbPzPqq8Y8/wLuLuOwPQ08vU/x13q0j3pSHqG5nxG44t3wvom7vQ6/OE/wfzZ/8T9oy9BrTYcEv2Nb/uFhgH2zsHsjF9a7Hsvod0U+Tf/1I/hajYxTz+L+9THME88CYcW0R+8hA4y5EufRJ58Cfd7X8EcPw4n3oFLbdxXfwv5L99Cv/YY8h//O/oHv+uHVC/M4f7OJ5Gv/zl85dNos4msr8F3nkL+7oNw6iwcvp3sfz+F2dVE/+EXME88BXcfQptNv2TWmXO4o/dg3nwTfvoqfOQo/PRvyd7ZwN4573/IgwvooQPo7QfLqjMv/cy/BqCzhXvmDV/GP/p95D9/A+0M4F9/DfPNv0DPX0YeOox75GHM68fR1gx66BDmR0/5vBfj2v70v2EWWz7xzEEvgfvugG4Ptvp+6a/pSdzRo/6YwQDzzE9hYTe8cgLmpiDNcC+8g1mYwv3+P8B86zu43/4tzF/8NcQR1CLo9mFhHve9l/1xX30U89rruI9+BPP0T3Cf+BjpH/8n0g5Ek1B7YB73td++oYx2NFz9zolF/cMj/4JzXcOZTkYvU9aTActyiS1pk2iPOhPUaOJw5ao+fRkwo5PE+WjfWCyb2mPTtMvjipeu9aWHI2NGZ/MMGrZMB4MhJfWrQ+Uvj4s0RjDU8+sl9POXrcX0ZYtY62SS+s98xXFLTCJ9Yq2XQ9d7slm+pE5x5WpRsUyUQ9sBzBWrMBUvpxvoFhkJFr+akpWYmjaJqefHuXK1+VQSUu37lalkuDJVJA2/yhUZhpg6E778mPK6lqgcjdyTDi5fAarIm9WoHLZvialrA0uEVUudGm1plytsgX8JX6aJH3auE2UZivot01ZDIn0svlwGQ6y1/DMmkYQ5prl3psnvHtzk175+J+a+f3VzhqtvpsJ3l1Le1gtc1JP03Sb9rE2adcncIB+DP9pkD+chGBONzU9wmqKaXrUkVzGHYHSewHAOQ3aNa/jreEaXC3MwOk9hZN5AscI4YmBkaa/y+2Kuwshwdr/jir/La7vy5RTDtIV8hTe0XJvKHzf8+/oM5zyYYT5G9l17uZor68GMpGO5qv7UlXm5sqyjeRz7TqI8Tb9MmrUNYjtJq7afevsTnOvWYfkGryNlh6I7+OAU/+OZz/ou4Z//W7+z28WcOeO7pshCZHEffRjqdcyJE7i7775xolcOu97YQPr9dx2pao4fx91zz3W/L7uhdxnGfc1zz5/3M76aze0df/aMn63l1HfD73K94s08eviwf9vM1pafRdXtIt2t7a1oubbmR1lHEdTr7zqcXM6eGZ9GkE8L2O7QfPP8C9deBqzbLWfKmedexH3q4zhzY1MhvJAucEsIL6QLfKAIogtUThBdoHKC6AKVE0QXqJwgukDlBNEFKieILlA5QXSBygmiC1ROEF2gcoLoApUTRBeonCC6QOUE0QUqJ4guUDlBdIHKCaILVE4QXaBygugClRNEF6icILpA5QTRBSoniC5QOUF0gcoJogtUThBdoHKC6AKVE0QXqJwgukDlBNEFKieILlA5QXSBygmiC1ROEF2gcoLoApUTRBeonCC6QOUE0QUqJ4guUDlBdIHKCaILVE4QXaBygugClRNEF6icILpA5QTRBSoniC5QOUF0gcoJogtUThBdoHKC6AKVE0QXqJwgukDlBNEFKkdUdfsHi1wETt267AR+iTikqnuu9cWORBcI3AxC9xqonCC6QOUE0QUqJ4guUDlBdIHKCaILVE4QXaBygugClRNEF6ic/w8Q1XFU+ezEoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wav_name= 'TpC-ord-C4-pp.wav'\n",
    "tpc_dir= 'TpC-ord'\n",
    "wav_path= os.path.join(data_dir, tpc_dir, wav_name)\n",
    "\n",
    "wav_raw, sr= librosa.core.load(wav_path, sr=None)\n",
    "# wav_trim,_= librosa.effects.trim(y=wav_raw)\n",
    "wav_trim= np.roll(wav_raw, shift=0)\n",
    "\n",
    "# Compute spectrogram\n",
    "\n",
    "scalogram= librosa.core.cqt(y=wav_trim, \n",
    "                            bins_per_octave=24, \n",
    "                            hop_length=512, \n",
    "                            n_bins=24*7, \n",
    "                            sr=sr, \n",
    "                            sparsity=1e-7)\n",
    "scalogram= np.abs(scalogram)**2\n",
    "\n",
    "# Second plot with log magnitude scaling\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "# plt.imshow(scalogram, cmap='magma_r', origin='lower')\n",
    "librosa.display.specshow(librosa.amplitude_to_db(scalogram), cmap='magma_r')\n",
    "# plt.ylim((50,167))\n",
    "# plt.xlim((0,500))\n",
    "# plt.clim(0,1)\n",
    "# plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harp scalogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 500)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAACZCAYAAAA1tzjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAFNElEQVR4nO3du29bZRzG8eccX3Krk+biJqGXpKpaisrUlg0JVVBV3CYWFlYkmJC6scBfwA6iAiY2hgomOpQiIQpEDYghQrRRUApKE+K2ceKkTs7LUKkwJCWv3PMg7O9nsnIs+6foq3OO7XNJQggCnNL/egB0HqKDHdHBjuhgR3SwK8Y8eWRkIExOjuU1C9rI1NQvSyGE6nbLoqKbnBzT1e/efzRToa0VC2fmdlrG5hV2RAc7ooMd0cGO6GBHdLAjOtgRHeyIDnZEBzuigx3RwY7oYEd0sIuLbmFZyexsTqOgU8RFlyZKlms5jYJOERdddVDZqZM5jYJOwT4d7IgOdkQHO6KDHdHBjuhgR3SwIzrYER3siA52RAc7ooMd0cEuLrq1Rk5joJNErum4/D9ax+YVdnHR9fbmNAY6CWs62BEd7IgOdkQHO6KDHdHBjuhgR3Swi7pNk24tK52ZUSiVlDSbCvuqCkPDOY2GdhUVXdgK0uy8kq6ytHFPSe2OdGBVYXRMKpfzmhFtJiq6ZHxY2fPP5TULOgT7dLAjOtgRHeyIDnZEBzuigx3RwY7oYEd0sCM62BEd7IgOdnHRBc7wR+uijjJp3qipef6CthpS2iWVqgUVHqtIh6pST5c0PKjQ368wNMShTthR3EGcQWosprpT61WSBvUtbGhP7bbKjabUVVTSv6ikv0/JviGFsVGFkRGpGPcWaH9RRZQmBzTw3kvae3VaWr8n7alKhw8qO3ZMEpfXwe7wQQJ2cdu+YkFhpKrw4tmcxkEnYE0HO6KDHdHBjuhgR3SwIzrYER3siA52RAc7ooMd0cGO6GBHdLAjOtgRHeyIDnZEBzuig11cdM3NnMZAJ2FNB7u46Eqcw4rWxUW3VFN66XI+k6BjRK26Zm8WdemdDZ354RMVnjku3a1LWZAOH5Bm5x/c2CSZm1OYmMhlYPz/RUW3vLmgD38taeKrPv30WV0DpU1thURnX/le314c1ukrF1SoFPTbNz06dPGgtt79WMXxXoU3XlU6/aOUpsq+nJbOv6bk1i3p8yvSC08rWVmRZm4oe/mcJCm5fl2SFI4cUXptWtnkhDQ4qPTatFQqKXvyxKP/T8Amcict1R/Nuj6dHdVmJnUVStpbCqp+MaD5tW51X62qmGaqN8safesj1RfKGupeVyJJC0vKzj2r7OKUipe/llbWtDm/ouLPM1KlT2rce/AuSWP9/ondkjS/oPT3RWn/PunGTam7rLSrLNVXlR1/XOrpefjIjcb2z/nH35OlRYWR6v3Hy39yk72cJSHi8l/7u/eHt4++rkM9TfWXmiqlmcYrdQ2PrkqSSnuCkrJUGCioeKAilQvSUyeUHT3694usrkqbm1J3t5JajSs8tali4cxUCOH0tstiXmj8iR69+UFF2amT//rcbKcFfX0PHoaxsZi3R5vgezrYxUWXZbtaywEPExfd7XpOY6CTxEW3teOeGrBrcdEV2QVE6/jtFXZx0Q1UchoDnYTtJeyIDnZEBzuigx3RwY7oYEd0sCM62BEd7IgOdkQHu7jo1tdzGgOdJC66xkZOY6CTxEWXJDmNgU7CPh3s4qLb25/TGOgkcdHd5cQctC4uujU+SKB1cdFV/uW6IcAuRJ4Nxok5aF1cdF1c6Aat4ysT2BEd7IgOdkQHO6KDHdHBjuhgR3SwIzrYER3siA52RAc7ooMd0cGO6GBHdLCLugtikiSLkubyGwdtZCKEUN1uQVR0wKPA5hV2RAc7ooMd0cGO6GBHdLAjOtgRHeyIDnZ/AfVeB+OMIVGQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wav_name= 'Hp-ord-D3-ff.wav'\n",
    "hp_dir= 'Hp-ord'\n",
    "wav_path= os.path.join(data_dir, hp_dir, wav_name)\n",
    "\n",
    "wav_raw, sr= librosa.core.load(wav_path, sr=None)\n",
    "# wav_trim,_= librosa.effects.trim(y=wav_raw)\n",
    "# wav_trim= np.roll(wav_raw, shift=0)\n",
    "\n",
    "# Compute spectrogram\n",
    "\n",
    "scalogram= librosa.core.cqt(y=wav_trim, \n",
    "                            bins_per_octave=24, \n",
    "                            hop_length=512, \n",
    "                            n_bins=24*7, \n",
    "                            sr=sr, \n",
    "                            sparsity=1e-7)\n",
    "scalogram= np.abs(scalogram)**2\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "librosa.display.specshow(librosa.amplitude_to_db(scalogram), cmap='magma_r')\n",
    "plt.ylim((50,167))\n",
    "plt.xlim((0,500))\n",
    "# plt.clim(0,1)\n",
    "# plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trumpet data matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x18629aa90>"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAACZCAYAAAA1tzjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAYbklEQVR4nO2deXxVRZbH671skIVASICwBVBBdgyCgCBhF0QR24GWbUBc8CPabsPQoo7tMjq0S8uILSOIgOintQVElB2CrIlByAIEEAgQEsi+kZDtvfmnJTn1O+bWvb53CaG+/9VJvftu3v1V1alTp+o63G630GjsxHmtb0Bz46FFp7EdLTqN7WjRaWxHi05jO75mKjf1a+yObBR6teznUw11/G8OI+WCY0VQx89JZ8xBXUKhjqjGaxecuEzKZdXYZlo1uwI2R1gwKZenl0GdKhdeS75+Yx8X1OHuvfR4ASkXVeHP3CII77OgNICUfRwYWQjt2gRsrrQcUq6uYv4X6R4qXA6oE+DE7/Nx0P+5UUAV1KlmnkNSQU6O2+2OgD8Ik6KLbBQqVkZPvVpuFVoCdVqve4iUNwzYhnUa04cevfNe/LLCQjBtGPMTKScX+kOdeQ/8Aja/qQNI+fR/pECdvNLGYEspDCHlbsz/e/vO8WA7PHw9KW+/1BTqPDEI73PDoQ6kHOKHDW/s/lFgK330U1Iuyg2AOskX6fNPL8PfrhPTEEL9K0i58005UOdyAV4rcs2ys2D8F6ZE53I7RHGl39Vy34+GQh25rYQHVECd+DzaWqOZ79ox7gDYUouoMPZn4Y+0Y2dbsJVtyyDlTiH4I3HCl++rcPpnzJ0i0R/eTMrdluyHOkeTsRM4dZk+jieizyt9X15WICm3/faP+H13biblWYdilK5d+vgyUv5w9y1QJ9BH6VJXMSW6KrdD5FbUPLAPxmBr3XcplZTndsFufE4iFasjLxfqtArEIfCFDf1puVkzqJM58QuwtehGxZmTij3ByuhYsM34OYaU0zLCoM6eXj+CbVrv06QcumIW1Lnl0c/A1mfHfaS8cWA51NnaEb/vjdGXwSYTFkCv9catu6HOS6lDwJaUGknK8z9GhbnuuhNsz9QhRFOi83G4RVO/yqvl6Ajsap2iFSlfqUYf4KeYDaScVBACdSZH54OtauEaUvZ9azbUOZTZAmzj1o4k5Yp7/wF18iqM51Q9h2LjCPsJG4ePH+3vN965HepENML7vF0qt2YaXvdQ7KWdfrR8bjz+f4H+9FHfHYnuC8eA7aNJef3gWKgzHgelOjElOocQwlmr4+rw/SSoM1cqr+m/A+okFzUi5TR0lcT0UHTaP95Au/a5b+Hn7h5nPCSFNisF29weJ5mad5FSCXbsot13k8G2e8gmUv45H3vW/5xwru6bFEL03n4f2ph6a/pT3y+9DB9r/zD6Iw/YNAwvVIGu0Oahu0j5AOMvoldbN6ZE5+9TLTqE1rQQx+IvoY77STqRaOKLPd3gcPrQu4XgP1KcgS3agSM1kLyrOdjKhv5Ayo19sU7vT8YYXvt0On6uD1Nv4KALpOzY1xrqOBtjz7pt0BZS7hiKM/+bNj4Its5NaL0gnyCoEy5N3n4YthfqVDHL8Pdto5OwgLH7sJJJTImuotop0otqhsKdS3HGN+tJWg6pNRz/yh2v0ImEa/hdUCd5xLdgmxvbnZS5VIXbXggEW8JC+oO3DcehZf1d6OM4JZHfcz8O+Wn3fAW2knIqzsG7x0Kd88wQnyaFTKKjsEfm8PelPd0YZoZ7auw/6WecOJJ0blIMNue2PaRc5UJn7ds70H2oC1Oia+xfJXq0z7paHjIcH7AsBB8m9vPlPPojTU7A74pqk4fXDsOeRibxHXSqo3dOIOUd0kxOCCH252LPekXqpNttxYnLxgwMh7xwPzdUS9dihuVHDD8lhHN/HNiO5dN76Mx87mwR9ZtH7hvN1EKSRq4j5dhsnExVon7rxJTorlT4iqPna6b6g2aPM/xMSaUf2AorjZ32oA4o1s+lGebJEmx1w1vgUC3Tt1Mm2MpdbcCWW07vvWMbnEhEMZMg3+b0c1yv3SoMe5WIr6eR8om7v4E6my5gqOXxYWlgk2kqxdvcC1dCHce8GWDr9Qh9Vuf+hn4fx7t13JIp0fn5uERksPH0vDadmhaAbchw+RoY7zu2D1vU8WIqsoQ8vJfXX8XYndwQ/Rph0+SCrkgMWB7CSiJhGJ2df3UWVxEeqMLGKMtp7yXs2Y8XYmMMWISzeBmntLrxP8s7QJ358/Bzx5bT7n7cgZFYicPnzd/8kynRlVb5isM5NWL4phvGjBJy6E0u6NkI6rR/jQ4tjrMYvO7SG3uV17fdb3iP3ArIvv+i97lghLXEVVlMQvArErKtPzMkugaiH+vIoEFsLoCLET8hNg/cSso/F2Bvv2ApDdFELT7OXAnvKawJ9Svf7Y7P/Pkj+Lm6MCW6IL9K0S8y+2q5Zzm21iDfcFIOZlYkfpCEkcCEFCa2DQdbT4V77N8eh87Tl6PoPQ5APyyRGQJ7bqe+YN+x6GdyyAHqyLVToM7RUevAtuEC7d3n//UI1HHdgz3ykB40TNTmDI4S656j68/3xf071OFIvUifQ0GFQgjBAHPLYC6HKK0ltM6b/gB1bpbKX/bdCXWOFtGvPVmIM9w/j8fAswotZkWC7anO1F+LfwRDEV+k4RD4nlTmfB6OlEv0QeEdCZGQi8I/VkCH/ZUvY6Oedg9eK/CTh0m5B/N9V5heWsa5BZ/V0L/QyWL3T84wn8SVjLowJbriSl8Rm1nzg3KzJJneYRieaN2YDrldQnAIDhjeHmyrJAFPPYgBzrWvYkucEHcTKbcMxulyRCMUnQzXO3HiebBvuuG1RkVlgK2pHx0CL15B0XG4315Byo752IsFSvFSeVImhBDny3Bi1tyfuiKPJao1vLowJbrmweVi+sCasHzsYHTaY/bQIOuWTOzqn0kx9gEKZiwH26kSuR9Fxo/8zeSGq+Rexvjin48Z31OzYIybhRTi7FXO8ljUFf2gF1/GcFPkFEUnXWLJlx1Iec58rFMppW5xGSy9/TGQPzya/p6Le2Gjzrpibsg1JbrcywFi9YGaB384H78sRirLyy+qNFmJLvMrCp/zG9MFbPJc9UAu5sBxKwsynG82kam3tj8NlsZno1/72TvYQ87AyyMfrAbTnJnmJ0YT4kaAzXEel+a+uJ/2fnHZKNbsCkxMqAtzIROHEBG1kvg6NzHu/qN7oGMfO5gGZ+XeUQghciatAlv4V9MNv2/HiyjyE8V0/XBkpDV/UfwvZrCIp1Apg9pdJOW00nZQp9JlLJR/9sN1693ZeK0P0ox76a2ZVOTcGq67Hbo0/SJoDmNpNU7w8iuw196KkbKrmMscbusQE/9a439xKS0yK/fikLg/m/7gMcznki5gFsZww28T4lwpzoR3X6J93WOJOAFS4asVrcA26Sms1/KbqaT8LHex9z83/L7EAmzUJ4qZOGkVHRYd5zHpIboZZqyoENSY9mKzD2NMlWO+p1KbLmdUi7hXatYf+8UafybUDwOxbYKMs/4GdLtgWIdjyqhTYCveSLNTuBWCLRdxiUvehx7ZWM13cW6LJWXXyBis9Ow0tEmMaYW9dqdg9Ee33kV7xIP5ODF7dhSd3MjZzUJgPp8QQmxOo8kKM9k7NYcp0QnB5+3XRRdmEXlyAv5zMlkZwWBbcjN1yNNL0L9YMgFM4skkOvxwQd6yasZBLqP/6wtHmXQghs+kyP6Mn7GOnFEihBDrL1CxLDqLiQKDme87NnotKecw8dM1ezuRcsYVbPicXzttNG3EiSMw3MSlYNWFueBwl1CS1u1MYYKXPWgmCLc/QMVp/5FJxjyST2fLRx1JUGfpvv5ge1oqs6sICvfEwf0G0x9MM/xcShH2RidL1LJKZLpuodOZruXo2DuPnyBldzOcTHHdyaebqHuUUoCNc5HCPdbG3DLY8QIyNK04w6SLl9JY2tNdrP2Qd7ZEZ7/SRVcnI7JRKk18jVMenJsxFcc1BmdzMlzvtD0LxfPWLzOpoRh7+/ujLoItzJ82tANSHqAQQgzYZZxkIQLQr417moY++sWqpV4G+lApRgbavCJRXOUrdmbVCC2pAH2OVMdBUl6962XD6zpjMZeNS1a8Scpsne3PJHpmYNBVbsFfLsAsl8nGOZwisRAF9lMeCgoIwVgel3Xd6dBhUv7TRLzPAWBRY08ODX73U/zcjCXSvV/GTsRkZpNwmDm16fbbu7jj4peY+oK/9cDAaEo+/c4fyw9DncWdu4JtlEIOmDMpGWyuXnTVtmzOMqjT+GPjTA0Ox0Xssf4x/igpLzqDyZ97co1n0M6D6Ay6+uLeubel4PN8hUA3+31M4//0Weo3z1TcRebrM+yg2+1mvRZzE4ncAuFcWbM5xjXjAcOP9ArFqXq4P+3+Q/Jugzp5FThJUGHLHIwLJuTThz7/mDWBcbhbYRhlaBQdvk+UYPwrZSQuqfWQsmg4gXF0DDL+reQlvMQ89On2ZGPPOvdWqcFU4aqF8DUnI1O1sy/6iI/frcnxekxhGW74XuNxiwseyOuJqly8gkPugWyFiPllJv4VhHsNVJDjdC8xdV7tjL0Kt1CvwuR449jZRmn75BlmoSipBNPJ5EmKJzC97zW73Jwj6fzqO7C5JjE7+iW4RWsVHhp7GmyX1tEZmLwFUgghPmeyTORIP7e1r/0GTDtXoUsI02NYxWmcid07lM78WzfCsEqHIEwaTRxB43nfpuPk8VKZOa/O3LESPYLFgnhz/sKCBRgyeR19aICLpa06Q4XBLf/4vPkw2J6Xklg50am4toeycQkIB048CSB01Uyo8xCTISOz/LZYsB0pRIG9c9r4mYyeQ7s2FddICCHW30FzH7k4rUNlm14tTAeHa+PciFm6rrE0U+LhW9DRViEuF3ueE8XWwi8yXLhAZTbHLZJzrDhMA7FynFCVg3kosF+4ZTDp3BfnmTSo8t5C2mCeUcxQui+OPk/VMDBzEMBVTImu4FgRSQffmIFd9GIpUaHjD/9m5iuuclcEzvgCnNhryvy91y6wPZGktl7oKVoGWJsEyYxqhcmtXUPRz1w4MJGUjxficNc8wNqpcHJMUyWeaYQp0ZVUOcWPWTWOemqJ2tEEVpBTxYUQoldWFilzI+Kui/jAn/DUTSnyx010SsCO3MypVOUv072pE+Jwlu04hWvLkwfSwO9RF2b3ft0DfTEVPppHJTJHIZ5phNfjdE+2x54npYTmvezKYxZMvchjbfGe/i/duDd05GSjLR2D0a4+XOIQ5Q8R6JpsK6MTlcKSTwyvw+E8ehRsrm7dqCEfRxLBHEgkb9LmgvYcnovTWSAaE4dFZCBjlOCWnFQ3CBsR3sjaUs7Skfgwd13CoWylFCrkdruNa4OuiV8Gt6HRPCAwIcRn0qQkMR+H2/eZiZmqyMzgddGp5l/JrDiND8VaMjfy9vsYsFaZ9CdgGEvsqzyIRkFnplwAefbhKLA9+gUN4JpdXqqLWGk+90sF9tp2YW7tNbWIZP1yGb+eoqXZk/ZM4LrX2n2PaIX+4pULOILIGb9fnEH5rMnBJuSaYryv1yrR4bRniyrHLB4O+X85kIuSUQnZ1MaU6PIqHOLrczWL3jEKn7E6TJr9R36FO+zvTDEVy7ILMZau/eBPmLvMDT6rpV1rwX7WGtBLtzAHIJ609rs8nWztczkV9N6b+NmcZRLVslIsfs54e11tFjEbyT01THLIAhNCiB/LD0mWGC/eAW6NnPob9YxIysPddhxy0gE3nMNsOZQ5XJxBPjXVE5jz6VqGCfEncz9hqJ/a3k1PMb0jrrNWnqJpo7JTLQT6PEIIcayCGuNy8SxfJZjwiMpDH9vG+DAgIYR4ZwRN0Hz+CIpu/Wi6wUYO+tqJ1ycSSyfhWqjsaKvi/J6e2SGYQOmwFbgRKOYmutl6VutYqLO9Am2ZBfLBgRZFxwhMJQVLNaidnGc85Vh0kjZGcwnmnsXrolM5UYiDO3h6TxbNINl/BVPFk3MGGV67fzj6JfkZuLMtPtSDi/ISr03CPXovpTIVFRjT2jjW2sSp1mvKFEyl2T5NV1tLxKiNKdFlpJQQR72KaWBH8ml279fZ1pZNqpnfsVMIFcvBK7gzSgWuB+FXLehhNfImaiGEmBhv7f+7JdhzglZJHuBmyyosPdSRlF+wdBWKKdEVVVWJrVk1Jxe18MEdW3uq5Eg781CktHPHsjVQZdYh46GMOT3Bq/ycj7l6VrPNJidYczHsZt6LNDjpidihud1gPr6ib2jNUkmvpngLrnQqsjdvxWn/JSk2e6IYfZ5Ndi+YKjC6ldqBkHLSgd0JB57ENe0aJ3G26xks3jPIp5splYeG4cbmUid9eBcELi8JgXs+VXAswURL9+PWEi1lhuy+W6nevizaGOth+2GP73eswmflnm0tS6guvD6R2JXJPCh5ixyXKm6Rh19rCbZlj3vs8kqsulj/hs6RYd+Tcr4TF/yz3Zidck5lHmjy+XlddNweTMDiXgSOGR2NA6of98YsE28EQesTvULpVkK3G7dFVrjwcB751NTvLmDcNUc+ht4A74vOZoZuND7URx7+hBBijjdu5nfCbW8U36GP7H7UOP//vZS+1KDY0OXUtL1l+NqgHCZ/ry7qrei4kyKnSS+IY2E2NssManF9vFv56QEnwHamGLOn1z+qcDGLo0m3pjRM5XDgWxCzyjqAbY3AN/L8iinRpSeXiHmdalraQouL8ipwb3TxFPV2KHXR/zmCyftzu/GUAXk7Y1IeLgVajdPJhw+p4uvzxm//zcyFiqorxPYi4xepeYJJCSqn0V2/OKTUeyGEcK+mGTkvpTI7gplZ572RNKvlsMD3y3o3zcIc5uJ0Tj8xIBCdzbp48Wb0Qf77F+/1kNcLC+7ENa/YPLrCsoc7TZE5v2VEJP1cs9x62pP/C3OpTb1CxGKF3eS12ZFvbQvidYPFkwES8zB7+XiVnAuodmKoysHh9QmvTyQO5CrsrL5OmN0mFmzx5fjyueQc6tlzQ+l3uZi97EylKeye9GqHhNGDE3fneX6lQZV6O3utjwzDuLNolINv05BfRbkpA3/mv5/HdHHXrbdavzkD7mhq/AZJDvkl0Q/E/35fu8GJzpmKvpKnHiYXsuEO/3muI/Vj48vsmXzVhdX0/62ZNF1d7TCKumlwovto8iWw5VfQ4U3lRSW/h8ERNELv6+BOPLk+6G28W9Q0DU50p4sxtpVcSM9A+bOX70EegjzRO1wrrv0eiWuMc+1GUnZNxEwUlYD13CgM4+y/jEPgwRzjY/c15qm3ouPenrglk8ajllucgPVvjvNCt9tc/FFjnd93VNjhRLCpnOOhAncUaUqFtReayMxgJgS//91+GlVMiS7rSIlY3KtmaEpnjot7C5MQLDEwHDdJRDLv2NJcf5gSXaWbvmaxfaBx+PLDnug/XZCC8W8xy2JKGSWa6xJTomvTI1j8Jd7cW4y3ZWLGQ6HL2svRNA0Dr08k+oRhKo4QnE3jTWq7RUIIcYp558p7Z+xZw/W66F49Ya5n1HiHrRl0xMl0M+ee2US9DZlohEffbbEuZ5RxJYYpLenaa3I1Hj8rJzgYoUVXj3ktGl/T9Mpxe0eO7s2oRMLK8HX1ZtGiq8ckMmnndrMg1fN+nhZdPWZMa7Vj1pb1oTu2rB65axdadNcI7t1g8qTrMcXF9jB/750u5Q206K4RyfmeGzqtnhx1rdCiM8HSPngywCMWh7LxbW/cn/7G/c8twB3p/4jFa81SfFlvQ0SLzgTDW6q982uDdP7H+AP1Z89pfaDBiY47NpZLZbKC6gZwl/qbr25IGpzoTpV476UnqlzLk8uvBxqc6KZ0NH790MaBW8E2dr+1ZSKNeRqc6LpsNt4Gs/MSnpln7dxPjRXqrejkvaNCeC71xpunTWmMMSW6nKPFxFH3lIPOMSQC3+qsaRiYewtipRDbMmumZt7czHK9Rdk16pgSXaCvQ0Q3Nzc7XM68h+tGDoxqTIquZfdg8YzBkf4yhZXXx1GrDZ36lIni9YnE9XZ2WkNlx0W6c8/aG9s8Q72dvWo8i1m36Fe8saSnRXeD8PwRayPOpkx63Ox4D9yLFp2mTj4863n3SHv5GtvRotPYjhadxna06DS2o0WnsR0tOo3taNFpbEeLTmM7WnQa29Gi09iOFp3GdrToNLajRaexHS06je1o0WlsR4tOYztadBrb0aLT2I4WncZ2tOg0tqNFp7EdLTqN7WjRaWxHi05jO1p0GtvRotPYjhadxna06DS2o0WnsR0tOo3taNFpbEeLTmM7WnQa29Gi09iOFp3GdrToNLajRaexHS06je1o0WlsR4tOYztadBrb0aLT2I4WncZ2tOg0tqNFp7EdLTqN7WjRaWxHi05jOw63261e2eHIFkKc9d7taBoQUW63O4L7gynRaTSeQA+vGtvRotPYjhadxna06DS2o0WnsR0tOo3taNFpbEeLTmM7WnQa2/l/Vp0M03q6iTMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "batch_str=['TpC','ff']\n",
    "\n",
    "data_keys=[]\n",
    "features_keys= features_dict.keys()\n",
    "for key in features_dict.keys():\n",
    "    if all(x in key for x in batch_str):\n",
    "        data_keys.append(key)\n",
    "\n",
    "# list(sorted(tpc_keys, key= lambda x: x.strip('#')))\n",
    "data_keys=sorted(data_keys)\n",
    "\n",
    "data_mat= []\n",
    "for key in data_keys:\n",
    "    data_mat.append(features_dict[key])\n",
    "\n",
    "data_mat= np.stack(data_mat, axis=1)\n",
    "\n",
    "from librosa import cqt_frequencies\n",
    "from librosa.core import note_to_hz\n",
    "freqs=cqt_frequencies(n_bins=24*3,bins_per_octave=24,fmin=note_to_hz('C1'))\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "specshow(librosa.amplitude_to_db(data_mat), cmap='magma_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harp data matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x18629a9e8>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAACZCAYAAAA1tzjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19eWBV1bX+d87NnEAYAmRgSKCMAgIyyGhAnHDAqYqgPpRif621r1ZtfdVnbautz9Lnz2rtsz4GqaBSUREtgoDMiMog8zwTY4CEkDm5Oef9sdY+99wz7hNCBHu/f25y777n7HvOOnuv4VtrKbquI4YYmhLqtz2BGP71EBO6GJocMaGLockRE7oYmhwxoYuhyREXZHBGRrqem5t5vubyr4twmF4rq+i1ebPzd66qaiA5yX9cdTW9JkmMdcDGjXtP6brexumzQEKXm5uJDZ+/2qBJxOAOpfg0vW7aBgDQxuaft3OpO3dC69XLf9zu3TSXHj0adJ640Ogjrp816IhNAOXYUXrdcxDA+b0RslCOHYVSXQMA0Lp2bbTj6iG+DXu+ptexjXZoG7RevYAa+g1ITHQfJylsyoEDAAC9SxfpOTS50Kn79gGnigEA2tAhruP0Dh2jXi8E6B06IpAr/U9/p9dH7vEel55Ox3/wrgbNKwjURUuhXdd4Um0VNnXrNv85NNrZAahbvoK65SvPMVrXrtA+3w/t8/2NeWrUPTYddY9Nb9RjWlH/xAyo+/bRgyMBtXUK1NYp53VOQaFdNxbqoqVQFy09P8fv1BFaJ++FonFXurqw7xDlxHHo/z6pUU8LAPWsg8c3+pEjCD17P/DGe/SPzPaazgJXWmqsZo0KTaPXKv7xqaly3yuravy5MKqfeMd3TKMKnTboMt8xek57Q19z3DpLS+k14E3S/eX9nKGuWQd9eD/5L5RW0uv5EDgAUGmjUt76JwBAn/J9qa9pd9wof45ps4FH7/Wfyt690odsVKFTTp0EAOgZjpZyZNyyL2jcZAeha+ANSv6fKQAA9YuNALwfAHXVWhozanigc2gjhtFNAIBH8/zHT74NAKAUFEDPzg50riCQFTaBkomvAwBazv0337FKs0QpPVb76HMAwLrN7X3HNq7QFRYB8Bc67ZbGN89U3va06/P9B7dt3fDz9Myi8wT4jp6cHOwcrDPKWsgz+68AANy3OV9qfGIz+W1Byc3wFDrD3ZPVHAAw9llWKa50/06jCp3e1lvYDHisZmKl0r+kC6//aILUIbW7b5E7NxruewIA7fqr5AdXVNBry5aBzhF+m1Zi9UlvoVM30I6REvJfi9SdOwGQyyTl1SnSc9GuGu35ud6KHmB9wk0AgMJb5/ges1GETvnLm3RiWZNfKL4OK4DYFqumb6Eh5z69JoesmuGGUHYafd9vIBsSNZrif9CKiPEw57JPAQCTNnoLFACo6z+HNvxy/+MzMofV0x8L3Mc0jtC1JqtJZu9Xt3wFrUtn33GGjsYuGK3fpVJzUZfQBdWu9r+gQSFWYcBbZwwqbMqJ4/S9HNKHPnudHrUh93t/T/g5J8/iSIbXWJ5v3SPTkZeaJT+5Eyflxr3IK5wwOn4503VokzmHhXUjKzwGUoP5ufQBvf3nsmI1zSV/ZLC5qAq0ywYE+44EhLCpCz4GAFx6SWGg74ff/ZK+37eP79jE/q2AL/yPqc7/iA9eLzUHtX8uADldt1GETuP93At6K9JrlL+8GczzXlUdbDISF0k/+A39kR/s0OdD4MwIb6YwWMpVZOnKGit6OECcJC8Lwz66xHdY2T/pGjV7yN8NBgTzBDTOSldWRq/N3NkRYsvRK+ogoYEY0Dye3qBWnjGX+28PND7qu8+Ru0F53N/dEBShKWTy1fz5EwBAvI9tpBw6RHOqlxc6Wf0sdTrt7UGsdFk0jtAJYRMUnTj3wyq/8Hc0RkFYgA7e9ihhe2kuAEAP82V6+G77uS26U0OgxMk/MqfvpNhri+51CP3WR0FDxFmeOKE/AO8brm7fAa03rVjxPyaLWkb0ah+ejqR7afXS+rs7upW/vkXH9PAeqAsXA9m0mATZBc5J6GxWmoewNRiSoZ2S9XUAgNZ3tgPgfMPORdgElFakY8rc4NZvU6BfrI6ykLqBhScBVl+DkCKSru8CvVMH33FSrqrMDOBQAf3dVEJns9IkVjoAUDduAmC/uErxacPvExTpPUjMVk0jcRgxvkGH8caf/g7djzFigngolcnXBGOneEBdvIz+aNE8cp5CMjz0TH+CrTY2X5794nesQZcBgyxv1tb6fq9xlybJlU5sC1Y4CdyakYsAACNWX+d5TPVJ0rEKPlghNQcBZQYFqGX0PDUnHXB5YByPvZqsSu14CfDQROk5qes30PecqF/tMuj1RFHkvaoAAfySEoP50qj6GqtBSlGRz8DGDoPJ6kwe5EErVpwkR+kIn3Ha7ykmeseXk6WPDQBKVgsActulnj9YajUxxhcyeSGAwAHePEPhcqr6n+lIzqcbref5x4FFmFAfM8SICTcq4uOl59K4YbC3yDHrt2yrb30AwNvVIgS4X4saqXOH7iKfW9BtrG4tMV5CHgtpQwkCGD+Kvm9S+h2Pzz4x7bbr6Q0JZq9wngOAcoSY4XqnTq7jzWFCET7Thlj3xnNAQoL00EYVug/mkaf7pke8x518j1aA1h66qrKD3CHLCunH3OAwRl24GACg3XgNkCS/epoRlyNhGIRCkXMGyB0QzBLl2AnvgZnRaoVSQMq5zKoBmIgWHkJnOOe7dQN20YMGCaFTV6yWc6JLPCgCjSJ0wjCoCId8RgIoLcXawyycHsNEGKtj2ir3QV0iVlhDLVMZKy3Kt3WcowUBSAPauj2eN9jqO1O27qG5SQod9hyjV69zdDBdqxp5lom2+QggIXQq+wxlHsbGiUiwUp0e70+BVhevQnaKvK+rQ7L7BXLMavLw6znOR1CiPFgqQjnW27aFnufvbrCdo01qMKU9IeBtSZHY2tRIZoKS1xaApB5bUy/lzDeEramt14JqE1nc5eZXrSnCL7+iH73M8n2lqAh627ZR77VPkQuDRW21ZvhRxTu2df+MIeakLlwMtCRHuC0hZfsOoIKZwvtIH4XGt7WPP8Ehes6VwcanSuSmim2vrAza0IG+w9W57wMA/vj3PDzyuP/hlw+n6z9m7TU+IxtZ6B54np5nDXBdaRKfvQNvF77v+JlV4ABgbxnpXJfDLljqmnUAgPo1+4COLoLlw0QOZByoKvQeziG3KEOBtzl11nz6zMPz74gBwfh+QbK7lIoKdwu8pMRYEbWJNwMAtj2yQuq4u8qIGTNGYmwwoautc35fKJEVEqtSaio+OUD6l0zYf3i7U8bfel/LzWAFX73nSujvr5Y4GkMktKgBk+FSk4LRlto1jHqvdesmN7ABTl7lzBkggd0bVr+oiWwq1I7ZL4SkVIOQIu83CCZ0CS65Vrx06yUV/sdQVdTp8jpdTpdS429hnYmlXzyN6r590JNpbmb9y3bq2e8CACo/KwEAJL0yxZFQ6pZnoffoahxDu/dW37mLFUh950NotzvZ3w2HunETDi6n6975mh10Pg+3jDGnHj2gvvMhAED3mtNgOtakUSfwd4utJVxISE5C/ad07r1nc6XnHnx7NYW6DJYDW1m736aVo7s5tl1CN9j8FA3KKHE8tFJUhLoXFtLh/0B+qOSH8+kcpnH1h8/QeP5f69rVSAlUl5O1q48xCR2zYISgJJk4Byr7uMxWlxtBUzlwGPXHz0ad2zjO9h2ASu9aDZyCuaXIDEBs0Z4mAqT69H2RN3k3UddTAkzdJ/vw1Ca67rPepfdC6axvesRiFw/9BNdOdd6R6p+aYRATxGr7YPf9trClWSUJHSE61iXp52ulA6JCXVaTvkNuRJisAmlGzyXOlqLeti0S7ieHqljSldVEWzcr7sqvJrtOTxszyv5m2MNF4KAyuKZI7juB0DgXBkjR6Qgl3CJ0maMVqXCbuGZxndPt5xCGACcVxV8FXL+YflfcJfSAaRKB/yOVCa5z2L42A1aKbY/2J6FdNs71eGIFP/YrefUmuNCZ8lJFQFvASPioqoJSSiuCVf6VUyeh7GRHpYMSb01KqdtP5wuFw4bAi+3Bc8sykw88EmMEX0/dHtmiaqYRny3hxegEFq/wkaOwCzw0EQpHYZzWA+Fw1j+lB0zz8B0aq2gvIDuZmMZGJMMJTCMXzN7Canc9tmvnU7b3EtLkfHpDWsuTbYMLnckaFEq1sCL1EaxkJye70tKV4wWeFqP65OSo/0PP0HIvy0CR8bs54kxEd4zPidZdjZhySiqUd8kXKZNrauiG63ZDc6lqoM5+F9oEengO/ZxyHTpf675LmDGgG0UuhErhKPh83rqnZgAAHuwXD8DZ2Zvymp3zJ5s5Nm4y1aeR2WQDC506j3Quxyxxs2/OwiYWKwlSvPO7RB6lELDqH1N9kqRXIj9eH2RnE9f+O42zrk4APClXZyYR163FnAgT2Eo0NaIdL84B7rnWc/5mHH+WVvT2Nye5W4Bdc6C+T66gyZ+Rk3fVp5wA5CN0c74kleOB6e6rrLieQlf7Q+dVmOY22LSbCFLsoUXxyPun/wO2dg7d52E/9B0aXOj0QfbEl8q5uwAASSOGRd60UNeFZVX/5AyEnnF3CQhhE/pNUld7LC/8GrmVxSoIwCgFZLVsAXhSrtRQtDiom7dAX0e/B3ddDQA4dh9tt58VZOGOdrT96xJ5IRmd2MmbYHrQhMdeBMhPnQEq6b1RrTlScL9cwlDnVHZVORhrAtbdoW2SXd8WBIz6wyUGDV9JomtWWy8R2gTQNkWeXhVc6ByevoS28v6uuJwUKBIsB+M8D9vPF9fJnhea8AKtcJrwXUmixS8G0/fEG4kJUDq2or8PkWXbcRpZs3+4sgC3bCe9R+ZWbNiUAwAYUF6AVGY3qeuYK8fxTO3GawwheGrHchqzj5y3frkfY++iuajHaPvXTELnlns7ofPXUA7QrzWucVfORiurjlzT8fkAgC4HF/r/UADdplJUpMmyweJ6t5E+oT7hGijbdp7T+fSpd9jeq5hCOkuzH1wiPRfA7tvSevWyWZ/CCTusTXvEj8pzPb7V6l1wgqIpI8ebLP7e9miDuPmhZvTw6q1b+c5bKSiAzjw9rcZO/3JzYnf88E7b3A1CqomYKvycwnXlOIfX5hn3wtOYsaBxwmC9vxdoeGBeGkMU3KueTn6ppLspvKQNugzPr6ZV4Xe/tG8xVj3RF1Yd8Cd0YcesmQ9tuDsP3upi+e/Z/P2KHONGC2EwU8yFsKpPkV/OURkXTuy//IO+370dlL1cBbMb6XZeRXpEPTrZkJmnccIwP/wGHUuiUFBwQ+I9oo9rt5hYj6fPRI/Zus01dVApLoYesLaHgDhmQi4XXiw5G/lM3KmSUlihHDhMr2uZPj6WfYF/ehvqgyxE5hDQItIZDfIA+8iatayBeoK4cUaoyiP9Us+iJCHha4z6zBT/VHZykUWOVzs+HBwxqT9LD4RyvAT6BFINlEo24Bx0OyE8+MZ+XTwRkOkSpCpVsCOXl6NyGUl0kskjoWe3ixqmrdwGuAid3kGOGiS84NhDq4A28eZI6YI4DkqbSkf8Rz4X3PmSlXaTvigiDPuvpQB81+NUaOOnr3XGy086KN8DnMNJX+7NxpitXIdNCB0Lm/bMLJu7Rzi06/62wvtCiyKFby4BAKiXd6djOuRhGKmMr82DMpf8dEo/2p4VEREx091bNocflFffpvn+8E7jPc1kFLolUkUhQDw7mNClpSHlOtpCzHqB1bhQR17irlN50JrNLNXK6ZsBACEmq8RPBFRW8Ot3Ufa5EYqqrUVcGi11yrCeNCfzcT8i6zM1mXQfpTUJypHyKsfiP1ZCqMjAKq6NB3o4k0W3LGuNAU9Gv1f5Q3LjpF2d6Xo91J07oY8ZSuedQ64TpTrCSRMqhXXn0KfegfB/0PHjb7DUJqmpMVZnK8NFXbXWpt7Uf01WtlpSYqySZteYVBLS17y9SpBpA2+vMqWytH6XBkqLE9j6zGn0zqe/P9tKkx8zk7hoGiJbujK2LPqLdXU4tpuc1l0caET6IHJUZ13f1jgWAEz8r09Rf5o86Z7PJyed5Hc5DK23M3mnc3ax7b1wFa/IHkq2OU6r3EzpR1HhLI+eDwnX0wOmedDUjWOL++FAzYqK8Yp5eVTrrPkpCXvin02+0wDM7WBCV68ZSSRrXlQwYjr9gHU/oK3t8pUUo1Nnv+vKwjgz6XXDEVv7MDt02d3RZ3p/Y4UayN52bGTLzOw+sOpPqalIb0ZblFOpCSfGCQDctXF0JNrg8Lm69jP6rBO5PtK72Itzi5L2ZueyQFwiHVXduxdaG7YmLfqsmbhqGCKmfAPxO5wiLYbOKNoflJXTGAdGtXj41Y2boAuuo2yNYguS72EDrkHfDip0VZXQxvBe/+J6aB/QTampjzbP9ZH9jb+NPgOsaB4saAWxWMe1jvZ2mRNLmv+clGRs2iM1tZSWtCUFrWuiHCeWhOOTqvNl5ZsT/28jbRfaq39Cchf6fVq3boYFDYuVZ34gvOra6aX22Kb13FI8j9q6BgubgExtaS8EEjq9tNp4Ukc8kw6AtrTLj+yOHpiYaNTCKFlDT22zLlRNad2p72EAm/9WxRuwpyeqJ+TKZqXd35e+J/lbBLwc1EKZdqXCI1LhEhVVqHydrNSkl2nlNhfZ0RXaatVisjDN1p5rCmFVlWG1nmuPCaPVgmTptUVDSQ++bn2AyqOSCCR0SkJkZdJGDTdCOvVzSOiEwo68HCCddJFQPAtYGp1qUp9DUD8iV4cTS0Q7EW3aV310GACQeA2gTGcfVS82ZthKUxcuhj6ceP/mMqcAoC5d4dltx+qPUgoLbXpo9VLaghNutHeI0bJIiVePHUd9lTM5Vd3ylUGA0BzcRcoWCrsZQsd+QnXVOmjXOBfvVRcvg3Y55zo4UPKtq6a+ns/hxGBh5pAybzH0SqJ6bSohL8N1iNTN08a7x509qxJYEGylq6mPuknh/yRP/dptNMGrXuWn4qW5RlZ7M0ty+196rMaTc8lCfb4XHesXOyMOSDWLXRD8f/k3ZO0mv/MhipeTlaUtpS23FTdAqa+oQ/0airkmTBocdb6S1w+jperu6Pzmr6QPtcuMCKuVOZw0mbeTWfMjFdN5u1x7M13sYavGIXV6tIVplFioqjaiGuEielCXrCQBG/fZWPvNFE7pFu7uDv1ECdRtTMcSK7IpKdy2RcfZTaWz9xJZNCWLrvayVR1xDa9sH/yWdqonX30bmsmVAiBC7JixAMpt+fT3AdbBh7pOOTIV/yERKMlxUTdOqDzJIdo6xUqn1YRdg9B902sMB3NZXZrtHPrIaPPcqHx06BBa9ueYoTUT69RJhDbQ1matSpnSth5afxYGhx4VGZdxEcUaEgZ18xajyDaE7pJGOlD48FnEsWdf70W6Y1LIvQij0Ndq/rAQ8e3oUieMyQUAXHHKtV+b4fPStxx0zWVdMSMNV6xhYeNYNjKZamZaWa21RdQ16wwh/XQP7RjjZ9Nqet2Cj6Gzz653Eq32yuhsI1ZePXcrACBcSSv66u0d8OenyaXz0Wl5KlkwQyIuWvGPa0n/92rPZM4wXeRv1oeQlctPnWULvXbYEWiHyYf08EBBAo2wKoQFJ1aS+pdJx1Ofus8w+63QM9pAZ1dOqGAevcefJU4dbmxpBlHTJHSG3sUPAurCQKLlsnBqYVz/LCCVdazmtAoNnEaWraMuyUIel64a5xHjUv0z9aB0aQdY4rn1T1KMefR/doqcU1QfFSzos+XGMQxDhUNWiijrCqB1QnSOqjb+WmNBmH4i35ivWBF1nROhTNXcd+BL/kviBzGCCZ0WfWlD7enCt27HP5rbNGVNbAFtrHPJm6S7+yO8iJ6YVm9Odj2V3pwEQ001BcslfH5Kq+ii21G5DzKJKw7+KYOQ2u9Sw2mKUSTI2iq+iU4OVBbuUIc0Yz5O1CvhfDb0N/bq6907o/xpWlnTBpIBMP/jXADAlrfi8Hv24Oh53LzvTc4kNmWHCSNF+YIbxXWIRI9GTeVICLuGtOGXG/dYuMb0kYPwwXZy/t+zgHymeju6D9fNXYjxT9OCEaQboqLr8gkVA/t10dfdQN7s+DFdjeiBsQKxf8mrpoYfrEWoxZMdeub+iJEgju9k+rMSrnzDzlAfp6WgviOToh3m8I9AVMlXS76Fupzn61HNXSkqwurbSE90KnnmRCQF2FiwGBL681Sdav2HGRjxMunStjh3TQ22jSPlv88N5Ei/4nekK6/aO8ym7zlZ5yI0pqQnQy+hld6zKqelqn1caPRGXdcds7oDbq9xiB/F2erNIjc8SNRBOXLEu7qQpW6GmUokOs8op4hHJpycyqFDkVAcK+FmYfOKjoS30xaveuRblHxFOkwr0/GNDDOJ1gF627YYMrjA8TP1vUV4cBmtVNa2H06Wq2A1j2j/vjupoqAAZ2qZ/PowERrum72CPis6aRM6J1fQ6RW0gGT8+lKgRQvH80QdI0ALheBhMAf3g5dT05qHqqzcCP1e+ZXQXKLfLWdAqapydYyqaz/zLu4swT/V6h1cIR5FvZ0Q/6doXtqhceT+6fJ4Nuq0gBXkAVSuOY0kl7J3el4ervghJ/mI81WwPm5i5nihopyEttU5dBdywznx6YyEnA6kTCPD/zt6WAvESJDJqPdsD15f78mSiOtGNCKztircQjUfkRX76GrSU2ZpWrCqAB60p6IK0tHyRg3HhNzl8sdk1Nd4J6xbE5NaJvBjWedSpUGA1YeMnHLvceeAgKQpPTrYW0sTFNuloQf0bG8wGaxxT/3+2yOxws+ZlLmGtr+EF6bYEn8qp5JOl3ZVG9cgdMnE142Oflayot6rm12HMTuErX0namux6RniB760h4Tt3s4cCzUJnJX9UffIdKzdQFt672xa3TPmkUKvzn7XKKIjWB9DVl5v/L7bprCaIKzrXj2Nebo5ttMGphhpjdbiktvHvo/eS9lQYZ/az7Yzh9ChVKzwS1a+vgXJT9E1Tn4lQgIQKZICTuXAHHNTXBBI6MKHivFGHbGEp8DubK3ZTU9H4jhvzpxyhDn9zL5IPENbjQ4YCr3A6u10rGt/5c5O3ng4y2h7r38dTSjF/OWA1blZbCJ/8opgLmE/4I8kPD/7BdV967OM9KK6x6Yj4WE2BLaTtaZyqY2XF3fF8q/phnY+Tg/hi3yO8NGzwGGiasWx+0UYLKcLU5HGJTFQxB0FM0hovSIpZzdUobyU1uecbC7zkEUPeJv0SHkPdTc50rVuHJN2KFoonOApNTVR7iQAVARoMD0E6350GAAw7DVmz5hqrhS9T6t6hkSl20BCp2kKWiW4J98mjadJ6G6FdsRxLBZiVA5pVvTKqHJhFuWL7caKumAoWUrjN5CivbwoKSJ0lux1pXs2lM2s3/Aqo5+O3BTBdEl8NKJMi224zzJ6NfSvCS2hieqawmdZTatgZRhIUOm9U9XRq2fogXHQZpE1eeD3xDzOo77A6HizClzOeqtwCrPBU//EDOqm7YBms+6Dcj/332KnsL6Kfme7+aYc28NEaFC4Y5FZL7ZWFdVGDHMsDStCfaW1zJrOzbXNJ7mZz7ZtQiChS+iQilvfZl8N7IFqvQPdkG0/3IbcjrSFpc2wc7UMOhH74sw6j76CVgRBZbr6Hk7iHXMl1AVk2s88QNO+gd0pX5XkRI5tyYfQ8kcaRW8guHZqRB/ShXxUu9c27jSYKUP3Rn5L+SdkPKXwlv+j/rNx9izFm4+djQ5fKVt2IHQTheeS1m+K+ky/ZzyUeezrG9gz6rNwmW5knVnTBuufnIG0GdFdbVSFfpfZqDK4fA7bqtM2aSVAaJNvM9SRa6eSs1t3IOKKDjsyCKbTJSZGPSnKVi6HIHS6neStrK6PQ8lJUpTtgS5v35nSj5NM+H9Bsow7csygVbVMIMEsP0LL/EcnI0wI5ZQ9CUdw+4Thc2Y1rXTNpwLJP+AURFH4cMmnNvM/bgj7w0zJyIkWil6rN++FUAw6C38i61PajhPQt9BWnfWepWNQfb1rNQIzSVJhBV9cl7isZChCNx1B1rknrStgo2MneNVhCdKtMrAhEYUqSxiFE14G3tiwH1j32HT86UO6cI9zvnNcH84BNVmeMwvyo75Xes8spL9MN05Ps4t5wc1ESW+ZSfrUnC2kH/4IMJi56k7eahz8TU4ZVKHfkPA4kStt1vRDE6BYEr6DFMwG7AKlP3iXJ39ORDn0HvRbnXyjwjFevZbo/4JMazv3Zo4Te1Q8DNLdMpjQKdFmumFNcnBfZactDh61OzY5WqEUfQOFtzK9JTkdq54mK2zJF52wtThaZ9RuvNr4W/TayniOQmxi1T17JhktznD5MF5dzHy17Hlc+Ia3hfu4VAVwhaG8yzRqUw4ciIR5hBD19PA5mspsmAv0AMDRX5L+1eHllAhjWKKxXxSspApT6qRx/S2hyyiaFTvEky5ltjXg2AtkzxKyrrv/2L05dBAE99NxjQt15CV47Dba4//4Z9aHktky6mBvYqsuZK5dThtjddH56S0+Tltxu6QaPNk3WiE1aEZDBuF0KY1rY9lqQiENOiu7Qs80o+qnJKyi90LSC6ZmdanyF1B/by3wKAmdIURiBTYRLiPH5qhNOGyL+5ZV8bX6+DNgKgmdup8sYmHwqPM/8syvUBeQf0/QrQxL1XQuQ4hYsPU2pvRGXgjMq+i+KUTV+t7HLHSlpeiYx4vKRlJrnEKFQfJqAwvd0cWk2lYsPIiVZ1hAKvhi81Pl5KzVjtFKpOZl27bp9j8iBSl79Ei8MICCz+II5f9L/rCU3r3Q42FaGa1bTXJKrbGKaRnRHmplxjtI7hqt+JrLz+utHUiVrzFTxVJJoHJPGIJ3q3Ww6KUeT7+6boMtA+uSx9nQ2XHcdOLonaRyWSESdnCBxKfsBlntLrqm4iY6ERrqXqGHXWTqm1cwdRNb9Sbi5ZpvaF7CQaXu3ms8rFYmjTpvobHbaTvJSoZ3N62o+Uqj43xy/ikV5Rg1kJ/2O9hfx0K3euTHWHOKbs8PepMjeNinZO3te6hvpM6cQwWoYa1pi9dkDz0AABbmSURBVBEricrJLcqZEmPLsAqFcAzTP0xjEo3b2qUbxWCMi7aPV8MePaCsJwW45r+5Jt0LU4Bro7daoSSnXdEyYinu5hp7bO1VTp3hWGoLgGNFA8PHufnvUJeuoPcsfrmkV6ZEOIACL7wBAFAmXY2EwZlRv+uTYZQ3e/XbPYyV7fBmelC7OQT1nVi++ezYNljgcapj5Swg+r6pvbKi5uKF4NurqC+cmIjuzTnbyXLROqWX4qszJHTp7WkJb6UzrfuNBcZ2UPcVCWLItKAM+ZA95yyYSS/T06u/8IbRw1UZyP5A07SiOsIA0C+jnAk9o00kSZthvvB6GVnH4XJaZRIAW3kIoXvqZmPB4lqocyg2aJBar7/KVV+rOVKLpBaCj2gpZxEO26joag7/v323bevt2pId46Ytu+stdExtD3fzHmVqcSBKygpV4dK+WHKcqE8PHOQuPw9NRHJvl4w5U1gwSIX3c4q93j+RVww1N+r9jh/eiZ8a5Sfo6Xh2wxI+Y8QxGz+JTX3zl12UaPXyroDg4a/lDDFzjqvFz2YOfalX9LafR4y7nQyV1OvYAHE6uchZWL/BNQeg5Q+72Y+fbqJezWAyajfuR8s3KalnKvQzHPC3WLjK1wW2B0CsLgZb2ITcj/jpZW4bgEjW/rTZPCeTEPMCsvcJUjfyBm/B/jLaWM2dEt1cXNrvZzsmV/khmNCJ5hwMo3ggxyGjJnRL9OZ+1TvsdV+zMTLGKT/z97Po2Ja6wtrQIREjId6ef+qkmwGUB2s154W/ThsxzLgJuke/CeF/RA/3JiSOPLyciEElqhOUv0+/IYUvjzKgC2rfIVKrtXa9Y8FqsRomuldKcCJSapWcUuD6LaCuWMeotuwGMzVEsXIchXrz+vxc3GepaiCDYEJnSV8TLgyRxyBQ+/B0m89Hn78CAKAM9s5LPbia/GxO/FPD1/Tv0W4K5S9vQr/buRhz+M01UB/h0FUNrSj6Vyy8TlbYzp22h6F6OW0vSW39S3hFHyxyi8VDmGL1A1dUIaEHRTBs9ZkLC7FnMsVV05I4oecICfLkzflSUzBCeL9hdeOUnU/3vY8jtZSv4ygPOplyW7lWzbpRFLtLiSOB3HYmYB8Oxjltrwmib5elumTSLT1sW039N7yFOCSamNnBXR+gLclQ2CVKVtXsr0SCy0qlPnWfQaVWRPGb7nTjnLZSp9VXxJS1IYMM35u6kAyVyjWkjyU9e7tN/zIiNWb/nhXJia71i5UjxzBtO1nj2Sl0g4+U04p1v8dWb0aXX3EOhzlB2qN/2perScCGPBPZHYSePPwesnaVvuQfbffcPt/zOyGY0FXXYPtYorD0Xnozmr3I4aVlLBispzhZa6EW7qc6uZN0i0wA+hXRKYTaqMhqVPcYF4z5Y/QqmjSmvafVpAgOmdCxsr37gVlJqVFsZnGj4sl1tH83uRj6HjtuY2gIOAocW6Xa0CG27UtAGzII8epKmjrz4dokkfCFl+2C6iF0gmY2fzoJ7W1vmAwtlwx/pbAQiSK7zSFzTngLxMOaM/4oNGs5WwkEY5kUlWGDSrpTb/OEOvnXJlMvJwXVSTiyJzQ3PtNDlimJQHWzZli6hnQcqyvIiW4d9bkINUmGnIz4rQeBVFiOPTfRKi2T9BMFvnZfjv4QAxdf7Tps2jjSJ5MvIT/gyVV0BeN6ZLg+aLUPT8fH62iVPVZFmqJM6ye9bVv0u/Fs1PwAu05noH3rQMImECwxZ2B3fcPnr0becHgajIlKdK92gluLJMCUoXQ5L/0iZ2LXPtcwVtnkmWg2ixyr0vFOEbI7zq6C9my9eTTQVQoKvCthspogCKIC/W48C+VmJrxa6PiHr58XsUjFcSR6aKibtyC8gKIHe9aRu6fXbC5HFiCfxQ9eqQCNl5hjzYTaS0+hY0GVDP9ifE7Q2zpz3pW/vhVp6sGrn1LF1nRFFdwo8OYcTV0iZqg/9zqUOE4mXkAr3YjV/ml1emamjVZlpn4JnXT/Lyl0lZFIW/7yd9pgzAPO231Sop2jpg/176io9e8Hld1JHZhzp0vGc127BTmhKnhuB9CAbDAztF7uK4YX69UY88wsAIjy9bhliim5GRHFXxSVeYVaHymDOrsKXdqz44zvebYb5y3kpbl5GNyKtpi/7iMrzTmD13oAFcoWcrLqoi4K6zvmveSGQeQTS76EPAHL5rZy1bEy37U3PFGOktPWyXfmFOFx4jMa4znfFfuYnDn5Nihf8m+QEboAcWszAlKbrCc9t5JT9SX+viMBs8dbbDHHV9P02+cnR0JrFl6X7HZi5Nnqq1ARdun2aIZg1nA5/dPPbUJFBW2/7ZkCJGLE5q1XhJLEpez8z3m2Q1ujK1Gocieb6icdkmmYoCFqy0R/gR/ULNqClcJCaMedmwU6QcZ6dkIwoTtrSl8rK4sIXdC+qQxrWp4jTBdNbF+Vyymks/cUPe3tm5k4dCe4dkfAvr4CD91zGOU76WHo3MKhrBYLm/5XqlusJZIVO2trHnJTaDu0rkHKuk3Qb3fW93J6ldnfPM7J605CV2IaL4on/o1qKesOglW5nVQRpwJhVoe2cugQ1GZcNtZhvJdObFQSaPR26ukmPc1DRyi4+U1kv39u9dQERLBe2bgJB58iHfLFHbSC1HJJ9TEmXdOpPK3YOvV9XADRoQ+F8OXhF/dC/DLHX8iEAmvEpM/8JRi7ztkK1W517yKYkJti20AMn2RtLdT9rDez/zAq3sqGjZLCTYMdji/b2wsgQ0Yp+Mb1c5uwmVrVGykLEudpsHPYK1O/tDIJ8gXe7ccFTPrX9VyHNzsb5TWk5LIcIqT5Nysumfg6aqrpC8uO0KwmTXU4L/vy5G35aFw16bT7dz12grrjle43ISEBmkQ1euUyOwHCOLVTCwbxPafeD4Wnfc9nfL+42DN86IYG63TK9r1Q2Iq0PgGX/MrkQ7L2WjUZI+JHK7vpaa5bvh/x41kXE0+O6YL0GkFb2z0sa2dqhI/IPct77vY8TB1NnvO3VtLS79SPMOo3eDSwE7DGiJWrLmuQwIbLfG6CZUdxYvYaZcFE+Yy3lxndD+v3kLrh9Hjqc4l4iUcjeRt6kX27dyuOo4TDDfrNwYSuLDIh/XgxlGasKWRF++vMTkTBGNb78tNomvhbNxFbpEU83eTi2jzc9Uh0UN0cnFdSSH/q/QR3SFxACreeanLMWrsvQofCv7J5vKQj00PYBEq20erVmqMXWrdurvQl9YuNEbeSxYMfCmgAKqI5S7jettMIoylcWG3cWOH+cUS8vcOZVso+StN7O39EhISeS6KFzkqYlUUgoSs5AUzOWgEAyErNwh/2k1NTLOHIpkmYuWZHZpDx0XEm64O1tUZFpTmHaEtbUU0hm+yU/pi0gnw/YjsILyfBVEcMQ6i7qEJOq19C+8P0v/ni8w0XKXsP/qoU2km6s6Mzue6vpUSsGeru3caqJwwX5QhlcuHAcWhX0m9u9SatDuYnXXmLroNNZ9Q0qJx7u+jHJKSi4mV8x2Q7i4PnXva7FXaXB7NLlA1bjU474Y20+gnXk7mfl9qdS4NZcjQAGKuhMfatD6D/gPRPRVjQnTrhdJWlrYCw3D/bCO2qfPqbqfJ6tr+3IJDQna7RcTaRVqX9JaURdmlrDmOxsEVVUWIYIaXaWsMH9MoIWsVe2UHM3+xkQD8VXfgw9AAJnw5A9GJSjvG23JmfNIeSWspGvsgTbjJIpi3iySrVPaoQffXgHvRZRkJnOHlXfh4Z4NFiyslAAaIfwrFjZkR9dvaLGrTob6kyz9du1c4DsJogOuuHSkY60Ip2lrgR7hZneCNdqzgPWpZA3bZTCE3g3mVckBuJiRj1XMvo44trUBtJydTycuk9CTdaYEPiJ93pxr28Jx36US7ZNTnaffrra45jaGvS0z7+mhgdL4rV5XuR8hAdf0dP3XPHWNAyW0NbF60DOvrZSsgfVb+Py3ylxtsKQZrDRNpw8idlJq0AELmgTmGrGi0iAEbpVJGc3dHdPFIKCqBwRppXQR9rxn5leQJahJwbefZqbW+Iomxht8W1+bYb7JTbES4hUQk5qAzCH4hK0s3Lv1bQfBa5X7S7I433rKEuUccFHU33JoDPNpDQpScouHIiXYhvZmRGiidbwkvvFx9AWKcna8tZjjXuYTqN6YYY9dXE6wtvQOGq7DZuWUEB0IYtJfZVhb+h7Tmhh3tXGfP82ohGuHHONxkABv8mkgdh6Gh5OdHzhYmJksx6bVwcUMgpmB5FpJygZzqHwaxxVxrMV8bhJh/9gD7rMNU+3Kknm76GdDX9boqP69om1BfSA23O07U5qzltU+8ZrPulQCCha5MXb/QzGL1sDnCLs19q3S0aUl6l1W/e4BUAAO1qf++1cucY16C5NvuTSNdl1rXiD/JNHmcnYzqh52NcbdMj4G+mZQmjx8lCszJQ9NRU6GPcKVOiT5jVb5Z9czI0n7bpgClkdWl31zEdFt5pey/5LqbqO6x0++bRCv7ME6Tm7KhLxMbxlq0UgL6S9FH1KG3V1QvJmk3IlwoQ2hBsezX1qUpKDrtSf1KfvNLQPTql8uoisfx6sTTMjUBEPoTak8Zrkn2pZOLBgLuLwAvqoqWeySnxlsIDxhbVT7LDj+hwKNsRSOQG57nHmzt2p10obRfdx5SaVCDXnrNctp52qY//RqL4RTHN4b9ebJibt+HOYdXDQ2OiAKXEuVd5ssKcR+l4TmZAhP+H/EvxfSRM9pfmQvk+BUMVUbLV6caZ2bQruYhPAKHTdhd65nxaQ35G6wGHZnAGps2GMpEE+f2f0dxv2iA3H2U2kW296o8kvEhzeimLah4/8b950HPsD/7JQnpiHtrP2W063dP/gjub2wvBhK70rFEBKS7BnatrVv5z2ro3t7W5Ljp5m9vaHCoRtnkdbWOD+/hnWVbvqkLKIY5ybCFXBByEzqAlpaaiaBG5bdpe4b7iCbeGcop03PqaekcHrO17vIqWHaDRXv78O57PQfw0ckVsC9NvkGUnHn2PDL4OHsWUlBnE0lG5Sn55HRyb95Wwy6Q6TPeyXnMnHcggIF09DH0sWTJJez92HabO/wjaOFpdElPcV7rTz1AaXcu5LHQHTrg261AXLYX2c9InB1RzHmteb8exZqT8bCRQSEq/0jq63L8ZZirP+qMk/Dd5rHQnHiYdq/2PaezKDzOR/7jvdAwB3nWUhGnYvn2uW2ZuWjxmFdOKpcDd+BFYPZLuyamaeKw/RdvkH936jgFQOvJOUUMG2XNX7wccVq/sluRrLX2NOvvM+WOwestWBBK6qlLV6BuvNvcIEQ3oZViMyXnucceozHwA6Ouur+iHTxle/IrD9F665r/Sad26QZtLGVWeW5kJFWH/G5zzErN9E2kVOFu3VerYAqc4hOeloz1/cBRuHElE1a2l/o3kVnF5tj8cnY7Lkih/Re/goexnsL9NGBnqUcdhWdcx5Z1Z4F1f/dB3Ll4IJHS19apRozeU3cyVo1/z/5cj8ddUfiKU5b66WKF1dTfBzT0MVu4m8/+mZk7V7+yoPkFn97ptZtpOGwfGrm0+orMPx4+La70v5cz+KwAAUx4j5X1RAd1wd9I5YeQHZPV3nfJP90Gsq/ZqTvOurinA2RRWa7xoZ0nsVGZrvuzUFudrZPE8tEurcBoljUBClxiqh36E/VO3ultqB3a1wiUcOtLz/VtxR07gnoMAwLDIdp+lVeIGyWSYGi5n7yV0+mpaqZR125EWL9eekg4up9+sYT/z+l+QsO2sOiV3fPb+Z17hPkQ4sW+bQw/CrcMfR7d0fyKqoO+LRnYZPZx/i7X4ZXb/JhS6pG4toXTi+rYe4aCeS25pUGqaE8wlqNRt2+mQ/rtfFFoOpRvg1e/i738lBXrK9pdwRRpZdUskji1CVneNXA4g33XczD/TjfrdE6QPpcX7tyowwyhH4TQHiz/xH+tVhGeu5P/cLUyrnhd61p46CgC1ayi2G8eOhdBv5Uu9OiGwy0TaC91ANrEVwvelnDiOqpmUKdY8Tt6VAQDqQO7N6pFSeHnbCI+sf0v/bVv0Ny2fRQ/CvM+7YLLXF/qQ7vbEbvLqi8w22abIQQrU6F26IH7oQd/jG41l/NpqNc6tNBBc6CQdsTL0oEDnTUmFwitcrQR5M+q7Hg5SAVFa4emebdA51b2dpsDpl2jL2X6MQmSvHznrKXQihCR8jV7FDp0QKEsLMEJ3XtAFNclnV0q4lbj/QbuGu6FxJeM8QikuRvER7pYdTOYC5Xo+8XwVdr3sr7M0/x4JZssi8umlKXLEOMMA8djqBdTNWyJtCFq4qzOO55EJwIvdyEcF0vr6u6aC4KIROj07G0lcO2VzSUCpCwDtxmvQadF033FCr+nNte8erfYhHVgho34UnAT6898Be5GdSydKG/wMvIC4aIQOyclGa8rvj3B3TDth19VUAb3nErnuy8lPk8YsRcVmIuToNu8H2n7MLQfc4JRk9F3Atyp0jokhXmCXyehfcvBb8jwZ6cFMfKdQkB9kemLFQAgsdNbOLecChQsk1qwlmrcIQLtCsFR9CuZY0eYfd/sPiqHJEMwYLi+H/sE66B+sa5ST68P7Qx/eHzWlKmpKA0wlHLZnmcVw0SBQ1aa+Ldrom4v/4Tvu7L0z0WIcl6wPWLVJBrae9zFccGi0qk0JOXLhoaW7OyH3GOlRA9zbvjcY5e8eBgCkyO6yHiXNGgtKYWGjluFqSmy4ghzVogft+UYwnS4UMjzxXsVTUkP1qNMa2Y1dWgokkVvi66MkPNJxiYDuBmsLeKnvfLwW+uTb/AcKSCR0n3dwsvy0nXRdffewoG2kXBDsFyfEA/X+3vrrHixH0QeN3I47Pd1YsbI7uRNDHREwJNcQ61Wv878uFxw4EnFUOyk3/hyFTSDwY+ZUut425u5b0GYcM3GDz8kdvD2KtkEXFCQSj6Igs8KZmoOcF/D1nHWZt/NXlKtYMJ7izDdtkI8DO+G8re0yzs/vEvRWAcv9S6D6JzORcjsRBbyqyxtoILOn1y3cldLlc+UE+VOLaiTq9kng4olIXOhopK3HjHC1iiMvUHmzDh69VgXUzVwQ0oXy74prvVdp0enx3tGmlqXngJjQXcBo/ugQYJpk+hcAFJOuqxQVBdJLZR39bk2IgyImdBcwtF69kDYjQLmALsRg+ez7XzSZ+6MhiAmdgOhX0Qi97p0glPET/4/SKBurUqkZgrM39OcHG437dj5wHk2jiwy1tVRu9YuN/mMbAD0zE3pmJk6WpeJk2bkVCPeDNv7a83r8c0VspRNg98Hh3+xHxw8d+mI0EvpNspf5/1dDTOgsSEw8v0QCrzIP/yq4qLZX7bczof125nk9R7v5ThWJ3aHOW2g0DYlBDhfVSqdVkXoc6EkpKzsvPjQBr4I/MTjjohI6cy1dWagnTvg3oIuhSXFRba8Ngb5s87c9hRgs+FaFTnt6JrSnz6+OJiqHxnDh4FvdXuPGcJvyb3MSMTQ5vlWhc2q7HsN3H995nS6GCw8xoYuhyRETun9hqF9slIo1q8tXQV2+qtHOe1H56WJoXOid/Nt6AgAK5NtxyiAmdP/CkCV6anfL1YCRRWx7vZAhUcj7YkRM6C5knM9MsG8R381fFcMFjZjQfZdQURFpN3UBI2ZICHC3Zq8mwhc8AvRc/TYREzqBi1nYLjLEttcYmhwxoYuhyRETuhiaHDGhi6HJERO6GJocgWoOK4pyEsCR8zedGL5D6KTrumM7oEBCF0MMjYHY9hpDkyMmdDE0OWJCF0OTIyZ0MTQ5YkIXQ5MjJnQxNDliQhdDkyMmdDE0OWJCF0OT4/8AXb7nlUPn/HIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "batch_str=['Hp','ff']\n",
    "\n",
    "data_keys=[]\n",
    "features_keys= features_dict.keys()\n",
    "for key in features_dict.keys():\n",
    "    if all(x in key for x in batch_str):\n",
    "        data_keys.append(key)\n",
    "\n",
    "# list(sorted(tpc_keys, key= lambda x: x.strip('#')))\n",
    "data_keys=sorted(data_keys)\n",
    "\n",
    "data_mat= []\n",
    "for key in data_keys:\n",
    "    data_mat.append(features_dict[key])\n",
    "\n",
    "data_mat= np.stack(data_mat, axis=1)\n",
    "\n",
    "from librosa import cqt_frequencies\n",
    "from librosa.core import note_to_hz\n",
    "freqs=cqt_frequencies(n_bins=24*3,bins_per_octave=24,fmin=note_to_hz('C1'))\n",
    "\n",
    "plt.figure(figsize=(2.5,2.5))\n",
    "specshow(librosa.amplitude_to_db(data_mat), cmap='magma_r')"
   ]
  }
 ],
 "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
}