02_predict_with_pretrained_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predict with Pre-trained Model\n",
    "  \n",
    "This notebook shows how to load a pre-trained model and to use this for predicting pitch-class estimates of an unknown audio file.\n",
    "\n",
    "&copy; Christof Weiss and Geoffroy Peeters, Télécom Paris 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "basepath = os.path.abspath(os.path.dirname(os.path.dirname('.')))\n",
    "sys.path.append(basepath)\n",
    "import numpy as np, os, scipy, scipy.spatial, matplotlib.pyplot as plt, IPython.display as ipd\n",
    "from numba import jit\n",
    "import librosa\n",
    "import libfmp.b, libfmp.c3, libfmp.c5\n",
    "import pandas as pd, pickle, re\n",
    "from numba import jit\n",
    "import torch.utils.data\n",
    "import torch.nn as nn\n",
    "import libdl.data_preprocessing\n",
    "from libdl.data_loaders import dataset_context, dataset_context_segm\n",
    "from libdl.nn_models import basic_cnn_segm_sigmoid, basic_cnn_segm_blank_logsoftmax\n",
    "from libdl.data_preprocessing import compute_hopsize_cqt, compute_hcqt, compute_efficient_hcqt, compute_annotation_array_nooverlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Specify and load model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir_models = os.path.join(basepath, 'models_pretrained')\n",
    "\n",
    "# Set model parameters\n",
    "num_octaves_inp = 6\n",
    "num_output_bins, min_pitch = 12, 60\n",
    "model_params = {'n_chan_input': 6,\n",
    "                'n_chan_layers': [20,20,10,1],\n",
    "                'n_ch_out': 2,\n",
    "                'n_bins_in': num_octaves_inp*12*3,\n",
    "                'n_bins_out': num_output_bins,\n",
    "                'a_lrelu': 0.3,\n",
    "                'p_dropout': 0.2\n",
    "                }\n",
    "mp = model_params\n",
    "\n",
    "\n",
    "### Select model according to paper's experiments (uncomment as desired) ##############################\n",
    "\n",
    "\n",
    "# Experiment 1 (Table 3) - Loss and Model variants ####################################################\n",
    "\n",
    "# Separable CTC (SCTC) loss:\n",
    "# fn_model, model = 'exp136b_traintest_schubert_sctcthreecomp_pitchclass.pt', basic_cnn_segm_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=3, n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "# Non-Epsilon MCTC (MCTC:NE) loss:\n",
    "# fn_model, model = 'exp136d_traintest_schubert_mctcnethreecomp_pitchclass.pt', basic_cnn_segm_blank_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=3, n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "# MCTC with epsilon (MCTC:WE) loss:\n",
    "# fn_model, model = 'exp136e_traintest_schubert_mctcwe_pitchclass.pt', basic_cnn_segm_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=mp['n_ch_out'], n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "# Strongly-aligned training (BCE loss):\n",
    "# fn_model, model = 'exp136h_traintest_schubert_aligned_pitchclass.pt', basic_cnn_segm_sigmoid(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "\n",
    "# Experiment 2 (Figure 4) - Cross-dataset experiment ###################################################\n",
    "\n",
    "# Train MusicNet & MAESTRO, test others (aligned / MCTC):\n",
    "fn_model, model = 'exp131b_trainmaestromunet_testmix_mctcwe_pitchclass_basiccnn_normtargl_SGD.pt', basic_cnn_segm_blank_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=mp['n_ch_out'], n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "# fn_model, model = 'exp131e_trainmaestromunet_testmix_aligned_pitchclass_basiccnn_SGD.pt', basic_cnn_segm_sigmoid(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'],     n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "# Test MusicNet, train others (aligned / MCTC):\n",
    "# fn_model, model = 'exp137a_trainmix_testmusicnet_mctcwe_pitchclass_basiccnn.pt', basic_cnn_segm_blank_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=mp['n_ch_out'], n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "# fn_model, model = 'exp137b_trainmix_testmusicnet_aligned_pitchclass_basiccnn_segmmodel.pt', basic_cnn_segm_sigmoid(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'],     n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "# Test MAESTRO, train others (aligned / MCTC):\n",
    "# fn_model, model = 'exp138a_trainmix_testmaestro_mctcwe_pitchclass_basiccnn.pt', basic_cnn_segm_blank_logsoftmax(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'], n_ch_out=mp['n_ch_out'], n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "# fn_model, model = 'exp138b_trainmix_testmaestro_aligned_pitchclass_basiccnn_segmmodel.pt', basic_cnn_segm_sigmoid(n_chan_input=mp['n_chan_input'], n_chan_layers=mp['n_chan_layers'],     n_bins_in=mp['n_bins_in'], n_bins_out=mp['n_bins_out'], a_lrelu=mp['a_lrelu'], p_dropout=mp['p_dropout'])\n",
    "\n",
    "\n",
    "########################################################################################################\n",
    "\n",
    "path_trained_model = os.path.join(dir_models, fn_model)\n",
    "\n",
    "model.load_state_dict(torch.load(path_trained_model, map_location=torch.device('cpu')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Load Test Audio and compute HCQT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 22050\n",
    "audio_folder = os.path.join(basepath, 'data', 'Schubert_Winterreise', 'audio_wav')\n",
    "\n",
    "# fn_audio = 'Schubert_D911-02_HU33.wav'\n",
    "# fn_audio = 'Schubert_D911-02_SC06.wav'\n",
    "# fn_audio = 'Schubert_D911-23_HU33.wav'\n",
    "fn_audio = 'Schubert_D911-23_SC06.wav'\n",
    "\n",
    "path_audio = os.path.join(audio_folder, fn_audio)\n",
    "f_audio, fs_load = librosa.load(path_audio, sr=fs)\n",
    "\n",
    "bins_per_semitone = 3\n",
    "num_octaves = 6\n",
    "n_bins = bins_per_semitone*12*num_octaves\n",
    "num_harmonics = 5\n",
    "num_subharmonics = 1\n",
    "\n",
    "f_hcqt, fs_hcqt, hopsize_cqt = compute_efficient_hcqt(f_audio, fs=22050, fmin=librosa.note_to_hz('C1'), fs_hcqt_target=10, \\\n",
    "                                                    bins_per_octave=bins_per_semitone*12, num_octaves=num_octaves, \\\n",
    "                                                    num_harmonics=num_harmonics, num_subharmonics=num_subharmonics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Predict Pitch Classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1920, 12)\n"
     ]
    }
   ],
   "source": [
    "# Set test parameters\n",
    "test_params = {'batch_size': 1,\n",
    "              'shuffle': False,\n",
    "              'num_workers': 1\n",
    "              }\n",
    "device = 'cpu'\n",
    "\n",
    "test_dataset_params = {'context': 75,\n",
    "                       'seglength': 10,\n",
    "                       'stride': 10,\n",
    "                       'compression': 10\n",
    "                      }\n",
    "half_context = test_dataset_params['context']//2\n",
    "\n",
    "inputs = np.transpose(f_hcqt, (2, 1, 0))\n",
    "targets = np.zeros(inputs.shape[1:]) # need dummy targets to use dataset object\n",
    "inputs_context = torch.from_numpy(np.pad(inputs, ((0, 0), (half_context, half_context+1), (0, 0))))\n",
    "targets_context = torch.from_numpy(np.pad(targets, ((half_context, half_context+1), (0, 0))))\n",
    "\n",
    "\n",
    "test_set = dataset_context_segm(inputs_context, targets_context, test_dataset_params)\n",
    "test_generator = torch.utils.data.DataLoader(test_set, **test_params)\n",
    "\n",
    "pred_tot = np.zeros((0, num_output_bins))\n",
    "\n",
    "\n",
    "### Use this variant for MCTC loss ###############################################################\n",
    "for test_batch, test_labels in test_generator:\n",
    "    # Model computations\n",
    "    y_pred = model(test_batch)\n",
    "    pred_log = torch.squeeze(y_pred.to('cpu')).detach().numpy()\n",
    "    pred_tot = np.append(pred_tot, pred_log[1, :, 1:], axis=0)\n",
    "predictions = np.exp(pred_tot)\n",
    "\n",
    "\n",
    "### Use this variant for aligned training (BCE loss) ############################################    \n",
    "# for test_batch, test_labels in test_generator:\n",
    "#     # Model computations\n",
    "#     y_pred = model(test_batch)\n",
    "#     # pred_log = torch.squeeze(torch.squeeze(y_pred.to('cpu'),2),1).detach().numpy()\n",
    "#     pred_log = torch.squeeze(y_pred.to('cpu')).detach().numpy()\n",
    "#     pred_tot = np.append(pred_tot, pred_log, axis=0)\n",
    "# predictions = pred_tot\n",
    "\n",
    "print(predictions.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Plot Predictions and Annotations (Ground truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAD0CAYAAACGjNCJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5weZb338c+PzSbZ9N47BNLQhIQgSpGiFEFAEEQ8FroclANylOPhUfA8HETFBxQRAgJSpEgRRKrSIUAIREqIJARCGiSb3svmev6479/MNbP33rvAbnZJvu/XK6/dnfnNXP2aayezc1sIARERERERKdihuTMgIiIiItKSaIEsIiIiIhLRAllEREREJKIFsoiIiIhIRAtkEREREZGIFsgiIiIiIhEtkEWkSZjZg2b2rY9w3BAzC2bWqhHycIOZ/d+Pe56WwszeNbMDi9//2Myu/YjnecPMPt+omRMR2YZ87AuQiGy/zOxdoDdQA6wBHgC+F0JYHUI4JIr7NnByCGGv5sjntiiE8L8NiTOzG4B5IYTzo2NHN1W+RES2BbqDLCIf1+EhhA7AbsDuwPn1xAvQGHfIRUSkaWiBLCKNIoQwH3gQGANgZk+Y2clmNhK4CtjTzFab2fLi/iozu9TM5pjZCjN7xsyqolOeYGbvmVm1mf13ubTNbC8ze87MlpvZ3OId63xMVzO738wWm9my4vcDov3fNrPZZrbKzN4xsxOK23cysyeLeaw2s9vryIM/GnKqmS0ws4Vm9oNo/wVmdqeZ3WxmK4Fvm9kOZnaemb1tZkvM7A4z6xYd82/F+lmSr4Pi+W4uVwdmdipwAvDDYt3/tRgbP6rRxswuK+Z5QfH7NsV9nzezeWb2AzNbVCzTd6I0DzWz6cU6m29m55ZrJxGRTwotkEWkUZjZQOBQ4JV4ewjhTeB0YHIIoUMIoUtx16+A8cBngW7AD4Et0aF7AbsABwA/KS60S6U7iMLC/LdAT2AsMK1E6A7A9cBgYBCwDriieI72wG+AQ0IIHYt58nP8D/AI0BUYUEynnP2A4cAXgfN8IVp0BHAn0AW4Bfg+cCSwL9APWAb8rpinUcDvgX8r7uteTL/BdRBCmFRM5xfFuj+8xOH/DXymeMyngYlk/xegD9AZ6A+cBPzOzLoW9/0BOK1YZ2OAx+qpGxGRTwQtkEXk4/pL8a7wM8CTQL3PxprZDsCJwFkhhPkhhJoQwnMhhA1R2IUhhHUhhH8C/6SweCvlBODvIYRbQwibQghLQgi1FsjF7XeFENaGEFYBF1FYmLotwBgzqwohLAwhvFHcvonCorpfCGF9COGZeop3YQhhTQjhNQoL8uOjfZNDCH8JIWwJIawDTgP+O4Qwr1j2C4Bjio9fHAPcH0J4qrjv/5D9BeJD10GZY38WQlgUQlgMXEhhUe42FfdvCiE8AKym8IuL7xtlZp1CCMtCCC83ME0RkRZNC2QR+biODCF0CSEMDiGcUVz41acH0BZ4u0zM+9H3a4EOAMVHBfzfIGBgPeeheFw7M7u6+MjCSuApoIuZVYQQ1gDHUbjTvdDM/mZmI4qH/hAw4MXi2x9OrCepudH3cyjc/S21DwoL73uKj0UsB96k8AePvYvHJfHFPC6pI80G1UEd+hXzWVeel4QQNkc/J20BHE3hfw3mFB9D2fMj5kFEpEXRAllEtoaQ+7kaWA/s+KFPVHhUwP+9R2ER2ZDz/IDCnc89QgidgH2K26143odDCF8A+gIzgGuK298PIZwSQuhH4Y7vlWa2U5l0BkbfDwIWxNnPxc6l8FhHl+hf2+Lz3Avjc5lZOwqPWZRSrg7yaeYtoLBQryvPdQohTAkhHAH0Av4C3NGQ40REWjotkEVka/gAGGBmrQFCCFuA64Bfm1k/M6swsz39j8M+pFuAA83sWDNrZWbdzWxsibiOFJ47Xl78Q7if+g4z621mXy4+i7yBwmMENcV9X43+mG8ZhQVnTZn8/J/i3erRwHeAkn/UV3QVcJGZDS6m1dPMjijuuxM4rPjHd62Bn1H3nF2uDj4AhpXJw63A+cW0ewA/AW4uE08xr63N7AQz6xxC2ASspHy9iIh8YmiBLCJbw2PAG8D7ZlZd3HYu8BowBVgKXMJHmJOKd5EPpXCHeCmFP64r9bzyZUAVhbvXzwMPRft2KB6/oHiOfYEzivt2B14ws9XAfRSem36nTJaeBGYB/wB+FUJ4pEzs5cVzPmJmq4r52qNYrjeAfwf+ROFu8jJgXqmT1FMHf6DwnPByM/tLicP/L/AS8CqF9ni5uK0h/g14t/jIyunANxp4nIhIi2Yh1Pe/byIiUh8zGwK8A1TmntkVEZFPGN1BFhERERGJaIEsIiIiIhLRIxYiIiIiIhHdQRYRERERibRq7gzkdenSJfTp04fVq1cDsHlz+rcuvq1Nm8KboLZsST9UqnXr1gD4HfEePXoAsGFD+sFc7du3B2CHHXbIfC2l3D6Xv/tuZvUeE+c5n0apu/n5c/rx8bF+nO/zY0ql5ftWrVpVK8a3+VdXU1NT6/uKiorM1512qv1aWE/T87d+/fpkX1VVFQBr1qwBoEOHwucObNq0KYnxc+frae3atbXS8OO8H8TlWrx4MQDLli3LnKeysjL5fuPGjQB06VL4FOSOHTtm8hCn6/2wU6dOACxdujSJadu2bSYNP0+ptvW0li9fnomFtN9WVxde+OD1Fcd43t55551aZc63t48jPw+kdeX58Jju3dNX7eb7ycqVKzNfoXaf8HqN67dVq8JU43Xo43PRokW1yrNu3brMz3Gened1yJAhQHacexre3p73uA28Xrp2LXxi8ooVK2qVvVevXpl9Xk7/OY55//3CZ5p4+8fj1vPas2dPAD744INaac2dW/g8kN69ewNpvfhcF8d7Xfp54z7h27xt82MwjvEx41/j9vI+7nXnacRjL99Ocb04P7en7/kp1Vf9q+cvrkMvj2/r06dPJm1I55fOnTsD6biK89WuXbvMV+/7kM4B+TEcz39eDj+38/aP473feh79/JDWnde5X5s8X5DWtdeV/xy3k3/v86jH+vkgrTtvUx8zcd35WPG+We5/lr1+vJzeHyGtTx8Xnqb3b6jdJ1zc3v699wnPq58XYN68eZlYj4nHjH/vdeD9xvMFad15H58/f36t/Pn3HuPzcqnrsPffUmuI/HgsFeNt53NsqfWK1323bt0y5Ynb1GO833nf8LqI94UQqkMIPWtlJsfMmvORg4dDCAdvzQRb3AK5T58+TJo0ieeffx5IOyLAU089BcCwYYVXesYLrsGDC++598550kknATBz5swkZvfddwfSDhhPInk+uZWaKOpakMadM7/P+cQI6WSYH1jxMT75eIwfHy8a/DifQH1SiCeB/C8HXpdxfp588kkAHn/88Uye4wWgD1q/CPnXe++9N4nxevCJ1PM3Y8aMJGb06NEAvPjiiwDsvffeQDo5QTrZxhcNgFdeeSX53uth4cKFAAwdOjSTT4Brr70WgNtvL7yONn+RBZgzp/BBYvvttx8A++5b+ARiX0DF6T777LMAfPGLXwTgpptuSmJGjRoFpO1/4IEHAtmFv09cX/7yl4G07jwW4K233gLghhtuAGDkyJEAfP7zn09iPG8nnHACkG1Ln8j9Auq/JIwZMyaJ6dev8GFpRx55JABLlhQ+pO1b3/pWEuN9wS8i//jHPwB4+OGHkxiv6/zFMa5fX9y99tprQDo+f/Ob3yQxfrzHeBk+9alPJTFed75Aueaaa4D0lwSAadMKn7B82223Aek8UeqXr6OPPjpTnm98I31L2fe//30AHnzwwUya/jPAGWcU3gT3q1/9CoCdd94ZyF6kfdFx+umnA3DZZZcB2Xo+++yzATjnnHMyMT7XxfF+UfT8eF+FdKz2798fSMdgPI/5nOoLDO8bvjgHmDx5MpC2u6cRj718Oz3ySOFtdvHi1xeO/guMj+V4oe2/HHudef7iC7kv6nyb15PPPwCzZs0C4KCDDgLgb3/7W+YrpH1pwoQJABx++OHJvvfeew+AESMKH6Doc2V8w8DH2P3330/sF7/4RfK9t4GPNR+nXt+Q1p2PlYkTJwKw2267JTFe1z6P//Of/wSy46pv375A2l6+4Bk/fnwSM3Bg4bNmnnvuOQCuvPLKTL4A/vznPwNpv/MxU+qmj7eF58t/2QU44ojCK7x97DzxxBNA2r8h/eUvvn5D6YW/9w2fG3y8Afznf/4nkPYbj/G1AMDw4cMz27zfeH1Beg3yPv7jH/8YyN6c8+/32afw2ULXX389kP1lymP8lw2/NsV16Hn1ucjLF/O+8OijjwLpmImv+T72jz32WACefvppIL1pA2n7eL/2MTdgwIAkxvvk+vXr40/SbKl61B/SuFrcAllEREREWp6G/O96U4h/8d5atEAWERERkXo15FHSbYUWyCIiIiJSLy2QG4mZ1VD46FIDaoAzQwjPNWWaIiIiItK4zKzZHrFoDk19B3ldCGEsgJkdBFwM7Fv+EBERERFpaXQHuWl0ApbVGyUiIiIiLY7uIDeeKjObBrQF+gL7lwoys1OBUyH7qiERERERaX5mpjvIjSh+xGJP4EYzGxNyLxcOIUwCJgGMGDFCn30tIiIi0sJogdwEQgiTzawH0BNYVF+8iIiIiLQcesSiCZjZCKACWLK10hQRERGRxqE7yI3Hn0GGwqvevhVCqCl3gIiIiIi0LHrNWyMKIVQ05flFREREZOvQHWQRERERkcj2tEC23Aslmt24cePCk08+yZIlhUeVV65cmezr27cvAG+88QYAN954Y7LvySefBGDz5s0AXHLJJQAMGTIkiRkzZgyQPmReWVkJwKZNm5KY/Lb27dsDsHHjxiSmoqIi89XrcPHixUlM165dAdiyZUvmvBs2bEhiqqurM/nx8sX58fi2bdtmYtesWZPE+L6amsLTK+vWrQOgTZs2tfLsx7/55pu1yjV06FAAOnbsmMmH1y3AvvsWPufF69nrx8sX14cPJI+N+5rvW7hwIVD+9X5eDj8+zrPXj6fvX+P/BvJy+LZWrQq/F3rbACxbln1Ft9eB5z1Oy8vsdfree+8lMZ7GihUrAJgyZQoAn/rUp5KY22+/HYBzzz0XSPt4ly5dkhhPw9vW04r7hpfDlZq4ZsyYAcDOO+8MZOvO68rr1c8X91Gvo2eeeSZzzKWXXprEtG7dGkjbyfvx5z73uSTmwAMPBKCqqiqT5mWXXZbEHHnkkQDcdtttmdgvfvGLSYyP+fvvvx+AQw89FMiO81dffRWAESNGAGnbnnbaaUmMj71BgwZljvnsZz+bxHh5vF59fMW8PmbPng1Az549gWz/87J62+XnBIBnn30WgPHjxwNpHznqqKOSmHbt2mXy423pfTVOw9N8/fXXARg7dmwS89prrwEwePBgAB577DEgW4cLFiwAYI899siUx7dD2jdHjhwJpP3mxRdfTGJ8vvAye569fwOsWrUqc/zq1atrlcuvB96fly5dCqT1DXDAAQcAMG1a4am+Un1+zpw5AMybNw+Az3/+88k+n1P79OkDpHUYz6O+bf369cR8HovL6vXj7eVjGODBBx8EYOLEiZm8xjE+L3i+1q5dC6TzNKRz/dSpU4G0LvfZZ58kxvut16eP15jPc378I488AmTHg5enQ4cOmZ9jXvb8tTHm27xPeV3G+fI681g/bzyu/Hrrc62PZe8bcbn69+8PpH0rzrvXj6fp4zIe756uf/V98Xzs85XPLT169ACy/c/L42WNr0H5tDw/pWLysaXq2evV8+U/x9eJmTNnAvDpT396aghhQp0JFVVWVob4OrU1VVdXNyiPjUl3kEVERESkXtvTHWQtkEVERESkXvojPRERERGRiO4gi4iIiIgU6aOmRURERERy9IiFiIiIiEhEd5BFRERERIr0SXoiIiIiIjm6gywiIiIiEtECWUREREQkokcsRERERESK9Jo3EREREZEc3UEWEREREYnoDrKIiIiISEQLZBERERGRou3tPcgWQmjuPGSMHz8+PPfcc2zYsAGA1q1bJ/tatcqu5+PfZLZs2QKAl8djly5dmsRs2rQJgFdeeQWAP/3pTwC8/fbbScxrr70GwJFHHgnAiSeeCMAVV1yRxPTp0weAPffcE4Bdd90VgJ///OdJzDPPPAPAoEGDAOjSpQsA1157bRLzyCOPADBt2jQAzjzzTADmzJmTxOyxxx4APProowAMGTIEgG7duiUxy5YtA6BXr15A+oxQ586dk5i1a9dm8rF69WoA2rZtm8R4XW/cuBGANm3aAGmdQlrnvq2mpgaAioqKWjG+b/PmzZl8Qdo+CxcuBKBDhw6ZfMb58X1+vLdjHOP7vB/E+alLXC5Pt7KyMnP8qlWrkhivM68Xz0ecH+ft43XhxwJ0794dSOsl369L5X/dunWZ/EFaZv+6fv36WsfPnDkTgJ133rlWTMeOHTN5LMXz6DHvvPMOAAMGDKiVln/1ssZl8Hby83kdLliwIInp168fkPY/T2unnXZKYvL14OdduXJlrXK9/PLLAAwdOhQgmVMgHedLliwB4Oyzz878DGn/2H///TNp//CHP0xi9tprr0xed9llFyDt+5DWuY89z3O5Pup12K5du2Sb9+28UnO459WP93oHWLFiBZC2wdSpUwEYN25cEuN5zZ87/jmfnzVr1gDZPupz26hRozL56tmzZxLj9eD92Ns/zrPX2fLlywHo1KkTkO3PPo58zHbt2hXItrvPlT73xmPvnHPOyaTr5bnuuuuSmPfffx9I6+fYY48lz/OUHzPeNwCuv/56AA499FAABg4cCGTb2+va8zNlyhQAevfuncR4Gg8//DCQjie/fgEMHz4cSOvFr0lxH83X/fTp04H0ehOXy+verx2l5g/vG/nrMsDjjz8OwFVXXQWk49P7JaTX2EMOOQSAk08+GYD/+q//SmJuvPFGAGbPnp2J9T4X5/mmm24C4H//93+BtO8DPP/885m8+3kmT56cbPN28eu394O4r/u87n3L+0Z8TZs1axaQXtf9+LgOPd7rrn379kB6rYR0HHjbel+N+8+DDz4IwIEHHgik7e/rBIC33noLgNGjR08NIUygHlVVVWHYsGH1hTWJ6dOnNyiPjUl3kEVERESkXtvTHWQtkEVERESkrO3tNW/bz68CIiIiIvKR+SJ5a/9rQL4ONrN/mdksMzuvxP7OZvZXM/unmb1hZt+p75xNfgfZzI4C7gZGhhBmNHV6IiIiItL4WuIjFmZWAfwO+AIwD5hiZveFEKZHYf8OTA8hHG5mPYF/mdktIYSNdZ13a5T0eOAZ4GtbIS0RERERaQIt9A7yRGBWCGF2ccF7G3BELiYAHa1wsg7AUmAzZTTpHWQz6wB8DtgPuA+4oCnTExEREZHG18yveethZi9FP08KIUwqft8fmBvtmwfskTv+Cgrr0AVAR+C4EELpVwMVNfUjFkcCD4UQ3jKzpWa2Wwjh5SZOU0REREQaWTP+kV51mde8lcpU/v2XBwHTgP2BHYFHzezpEMLKWkcWNfWvAsdTuNVN8evxpYLM7FQze8nMXlq8eHETZ0lEREREPqwW+ojFPGBg9PMACneKY98B7g4Fs4B3gBHlTtpkd5DNrDuFlfoYMwtABRDM7Ich9/b54m3ySVD4oJCmypOIiIiIfHjN/IhFOVOA4WY2FJhP4W/evp6LeQ84AHjazHoDuwCzy520KR+xOAa4MYRwmm8wsyeBvYCnmzBdEREREWlkLfE9yCGEzWZ2JvAwhZux14UQ3jCz04v7rwL+B7jBzF6j8EjGj0II1eXO25QL5OOBn+e23UVhVa8FsoiIiMgnSAu9g0wI4QHggdy2q6LvFwBf/DDnbLIFcgjh8yW2/aap0hMRERGRptMS7yA3FX3UtIiIiIiU1dBPtdtWaIEsIiIiIvVqqY9YNAXLvVCi2U2YMCFMmTKl5D7/zWXz5sKHn7Rqla7vt2wpvO/ZG6+mpiazPbZxY+GTBdu0aQPA+vXrk31+7vbt22diq6qqkph169YBsGLFCgCWLl0KwJAhQ5KYu+66K3P89ddfD0Dr1q1r5aNv374AdOzYEYCjjjoqiamsrARgxozCp3QfeeSRADz++ONJzF577QXAoEGDAPA2jcvl5/Y685hSvw1u2LABSOsnlu8vfvymTZuSbV6Hnn7nzp2B7MDydvHjfV98/nzePA3PH0C7du0ysf51/vz5SUynTp2AtJ06dOgApHUSp19RUVEy7TjP+fyVivU68PPF5cqXtdxv5GvWrMmU08sQl8P7Ybdu3Wod7+PA8xGXIV9mz48fE8d4Gt5/477lY8Pz6jFe7+XKWKpePI++z/MXb/PzlRrffh4fX56fuI96H/I6LNcW3pZuzpw5yfc+5n2b/xyfJz9f+fwRl+uUU04B4LLLLgNg5syZAHTp0iWJ8fHt43L16tVAWu9xWsuXLwfg4YcfrlWGXXbZBYC5cwvv1f/ud79bq8w+7+TrMK7v/Dzh/SYe5/59Q+46eRuUuiZ5/r0N27ZtWysm33+8fteuXZvEeDk8P6X6X358xvONlzWfR6+nOF3nc7e3V1ye/Px59913JzF77FH4nANvi1/84hcAXHnllUnMiBEjMudZtGgRkL0WeV79PKXk6z4uj/O6a6wFkterjx2/DkI63+XnrThtrzPv/wsWFN7qNXLkyCQm3xdKyc+R3k7xtdrrLt9v4n7t5fE0fV4sNefmrVq1Kvk+vi7VxdP3fhTP2fly+TXDzxu/SnfYsGEAtGrVamqZdwzHeQvjx4+vN39N4cknn2xQHhuT7iCLiIiISFkt+DVvTUILZBERERGpl55BFhERERGJaIEsIiIiIhLRIxYiIiIiIkV6zZuIiIiISI4WyCIiIiIiET1iISIiIiIS0R1kEREREZEivQdZRERERCRHd5BFRERERCJaIIuIiIiIFOkRCxERERGRHN1BFhERERGJbE93kC2E0Nx5yNh1113Dfffdx8UXXwzAzJkzk309evQA4KyzzgLgqaeeSvZNmTIFgH/+858AHHTQQQAsXLgwifn6178OwIIFCwDYY489ABg5cmQS442/efNmALp16wbApk2bkpiKiopMjIt/3rJlCwDt27cHYOnSpZnzAWzYsAGAtm3bAuBtsX79+iSmsrIyk6afN+b7PI+ejzjPVVVVmW0rVqwAYN68eUnMgAEDAHj77bcBWL58OQB33nlnEnPmmWcCMHfu3Mx5Bw0alMT4cTfddBMADz30EADvv/9+EjNkyJBMGXbeeWcAJk+enGzzunOHHHIIAI8//niybdGiRUBaL56/XXbZJYnxPrT77rsD0K5dOwBuvvnmJObWW28FoHXr1pmvq1evTmImTpwIwGuvvZaJ8ToAWLx4cSbP++yzDwDPPvtssq1NmzYAjBo1CoD58+cDcNJJJyUxDzzwAADf/OY3AejSpQsAs2bNSmJOPfVUAM444wwARo8eXev4d955B4BevXoB0Lt37yRmwoQJAMyYMQOA66+/HkjbFtJ69Db9yU9+AkCrVunv1l4PPtZ8XMX5GTZsGJD2t3HjxgHwm9/8Jok57rjjAHjmmWcyaY4dOzaJufzyywH40pe+BKRjOp4n/vGPf2TK9/zzzwPQuXPnJMb7wsqVKzPlietn8ODBANTU1ADw4osvAtk2vuiiizLn8T4+fPjwJMbLvGrVqkxacf14ux511FFA2u+8LgAGDhyYSctjP/OZzyQxPnf4+Fq2bBmQ9jVI++Jhhx0GwAsvvADAOeeck8QcfPDBANxxxx1AOn/1798/iTnmmGMAOP/884F0PvPxDuk87HedPD+dOnVKYnzs+riYPXt2rRiflz2me/fumfPG3/s86uf1uQ7gb3/7G5D2x7hv+Vzi9ezn83JB2hf8PE8//TQA11xzTRLjx/u49rz7GI7T8jJecMEFAAwdOjSJ2bhxIwDf+973Mmn79QLSvtihQ4dMmeN+7N97vTz66KOZ8wPMmTMHgLPPPhuAxx57rNZ5fJ7zuvdrmacdl/GUU04B0rr3ORfS+vQy33DDDZmyAHTs2DGTfnV1NZC9hvz1r38F0nF++OGHA/DrX/86iVm7di0Ap59+OpC26e9///sk5rvf/S6QXht9XRHP/ZMmTQKgX79+QDpmYl73vp6YNm0akB1706dPB9I5c8cdd8zkC9Lri7e374uv/XvvvTcAH3zwAZDOCfF6Lt+ffd2z2267JTE+969cuXJqCGFCrULldOvWLRxwwAH1hTWJO++8s0F5bEy6gywiIiIi9dIjFiIiIiIike3pEQstkEVERESkXrqDLCIiIiJSpNe8NSIz6w38P+AzwDJgI/CLEMI9TZmuiIiIiDQu3UFuBFaoxb8AfwwhfL24bTDw5aZKU0RERESahhbIjWN/YGMI4SrfEEKYA/y2CdMUERERkUamRyxKMLP+wOA4PoTwVN1HADAaeLmB5z8VOBXS9wyKiIiISMuhO8gRM7sEOA6YDtQUNwegvgVy/jy/A/aicFd593hfCGESMAkKHxTyYc4rIiIiIk1Pd5CzjgR2CSFsqDcy6w3gaP8hhPDvZtYDeOlDnkdEREREmllLvYNsZgcDlwMVwLUhhJ+XiPk8cBlQCVSHEPYtd86G/Cowu3iyD+sxoK2ZfTfa1q6uYBERERFpmcys2f7Vk68K4HfAIcAo4HgzG5WL6QJcCXw5hDAa+Gp95a3zDrKZ/ZbCoxRrgWlm9g8guYscQvh+uROHEIKZHQn8PzP7IbAYWAP8qL5MiYiIiEjL0kIfsZgIzAohzAYws9uAIyg8Guy+DtwdQngPIISwqL6TlnvEwh+FmArc91FyHEJYCHztoxwrIiIiIi1HMz5i0cPM4kd0JxX/fg2gPzA32jcP2CN3/M5ApZk9AXQELg8h3FguwToXyCGEPwKYWXtgfQihpvhzBdCm/rKIiIiIyLagmV/zVh1CmFDHvlKr9vwLH1oB44EDgCpgspk9H0J4q64EG1LSfxRP5qqAvzfgOBERERHZRrTEZ5Ap3DEeGP08AFhQIuahEMKaEEI1hTexfbpsWUMo/1Y1M5sWQhhb37bGMn78+PDcc8+xZMkSAKqq0ruhn10AACAASURBVLX5W28VFvoDBxbqoU2b9EZ2q1aFm+Ft27YFoHXr1gBs3rw5iamoqABg06ZNAKxatQqAWbNmJTFdu3YFYO3atQD07NkTgIcffjiJ+dznPgfAffcVnjxZsKDQDmvWrEli/vrXvwLp8zqDBg2qVR7//vjjjwdg//33B2Dy5MlJzO67F96IN3164VGaKVOmAPCd73wnifH0vT569OgBQGVl+reVXvZOnToBMG3aNAAGDBiQxPhxnmc/ZsuWLbXOk+83pTqwH+expWI2btyYifXzQ9qm+TTiNl26dCkA3bp1y+Q9zp+3i8cOHjy4Vn78nOV+O/b+4vXsX1euXJnEeP3myxrnx/fl04rr2XlMvi5LHe91GafhZe7SpQtQun7LTT41NTWZn5cvXw7AihUrauXD++iYMWOAdJwB9O3bN5OWl+drX0ufwHr99dcBePfddzNp+viAtPw77bQTAIcddhgAjz76aBLz3nvvAXDHHXcAaZv06dMniZk4cSIA559/PgAvvVT4n7sDDjggifE5xOvg7bffBmDu3PR/8vbYo/C/eFdeeSUAn/50Yb4dN25cEuPtkp9bOnTokMRceumlQFq/fr54jttzzz2BtK+2a1f4m+drrrkmifF9PvY977169UpivO69LbxPxP3Aj8uP97jdndfv+vXrAVi0KH20z/PhdXDuuecCaX8EeOKJJwCYM2dOJn9x/XzlK1/J5MPreejQoUmMl9HnyGHDhgHwr3/9K4nZcccdgXTOHD16dLLv8ccfB+Bb3/pWJu8PPfRQEvPCCy8AcOGFFwJwzz33AOn8Dmnf9PHVvXt3AKZOnZrEeBm9H3qdjRgxIonxfuLXND9+/PjxtWK8fr1cw4cPT2Ly7ZvvG5D2G2/fBx98EID+/fsnMV7X69atA2DIkCGZ/EE6F8RzNMBdd92VfH/GGWcAaX/ZsGFDJn+Qjj0vx+LFiwE466yzkpirr74agOrq6swxfj5Ir2/XXnstkNbvsmXLkhi/dvj84/tKfSbDyJEjgfRaMH/+/Fr7HnnkEQC+/OXaHxrs9eLj2vtznGdvH0/Dx4HP5QCrV68G0vbJz/NAsobKX5Pat2+fxHzwwQd+nqll7s4mevXqFb761Xr/tq1JXHnllXXm0cxaAW9RuDs8H5gCfD2E8EYUMxK4AjgIaA28CHwthPB6XWk25DVva8xstxDCy8VExgPrGlYkEREREdkWtMQ/0gshbDazM4GHKbzm7boQwhtmdnpx/1UhhDfN7CHgVWALhVfB1bk4hoYtkP8D+LOZ+e3qvugP70RERES2Gw183KFZhBAeAB7Ibbsq9/MvgV829JwNWSC/CowAdqHwIPQMGvbssoiIiIhsI1riHeSm0pAF8uQQwm5AcivazF4GdmuyXImIiIhIi9JS7yA3hXIfFNKHwrvlqsxsHOlrNDqhT8QTERER2a5ogVxwEPBtCq/L+HW0fRXw4ybMk4iIiIi0IM38HuStrr4PCvmjmR0dQrirrjgRERER2fbpDjJgZt8IIdwMDDGzc/L7Qwi/LnGYiIiIiGyDdAe5wN8m3aFMjIiIiIhs41rya96aQrlHLK4ufr1w62VHRERERFoiLZAjZjYMuBz4DBCAycDZIYTZTZw3EREREWkhtqdHLBpS0j8Bd1D4BL1+wJ+BW5syUyIiIiLSsvhjFlv7X3NoyAeFWAjhpujnm4ufeS0iIiIi2wG95q22x83sPOA2Co9YHAf8zcy6AYQQljZh/kRERESkBdienkG2EEL5ALN3yuwOIYRhjZmh3XbbLTzzzDNs2rQJgHbt0g/tq6ioAGDLli2ZnwHWrl3r+c0c5+cBaNWqlWc6E1uqwT2Ncr8t+bkrKysz542/96+evzVr1iQx1dXVAAwePBiA1q1b1ypX/H18vlWrViXb2rdvn4nx8sT58XJ4uTxm9erVScySJUsAuOeeewD4yle+AsCee+6ZxIwdOxaAESNGZPZVVVUlMbvvvjsAbdu2BeCtt94CoF+/fklMjx49MuXztDt27JjE+Dk9rzU1NbXKlW87b+ONGzcm2/x7by9Pw88HsG7dOgA6depU8ryQ1p1/zfenOI3NmzcDaZvG+fG28OM3bNgAZPu6H++xK1asANK+Bmm7+/Fe3/G2ZcuWAWl9ez8EaNOmTeZrvlwA69evz9TBSy+9BMBuu6WfNP+DH/wAgJkzZwJw9tlnA7DvvvsmMd6WnobX2W9/+9sk5oQTTgDg7bffzuRj6NChScwLL7wAwHPPPQfAeeedB6TtB/Dwww8D0Lt3bwB23XVXINvXu3btCkCXLl0yacXtnh9H+TaBtP8uXVq4T9C5c+daMc7L7v1u5cqVyb7XX38dgH322SeTtqcZp5XPXzme9zg/+TbwffFc6fvyc1tch14Ob1sfw3Fa3bp1K5nnuB9723kbeP7ivtqhQ+FlSr/85S8BOPPMwn9ier0D9OzZE0jH3OLFi2vlwfNWas72ff61XP3m59H4PF4vnkapa4m3vY95b9u4vf14bxcva/fu3ZMYT9/rys8Tzzded76vVJ7z89YTTzwBwKhRo2rluW/fvkDa1+M+4XOrp1VqfeHpex5LXYc9H14H/nN8nZkxYwaQtq/3o/h66PXr54nntnx+8uOi1Jztc6XPr/H5PCaf91Lrgnz/iXlMflzEaweft/J9NR7DixYtAtLx4MfH12G/PvTt23dqCGFCrczk9OvXL5x66qn1hTWJCy+8sEF5bEz13kEOIQytL0ZEREREtl16xEJEREREJGd7esRCC2QRERERqZfuIIuIiIiIRHQHGTCz3eraBxBCeLkhCZhZDfBatOnIEMK7DcqdiIiIiDQ7fdR06tIy+wKwfwPTWBdCGNvwLImIiIhIS6NHLIAQwn5bMyMiIiIi0nLpDnKRmXUHvg6MKG56E/jTh/xwkCozm1b8/p0QwlEl0jkVOBVg4MCBH+LUIiIiItLUtrfXvNVZUjMbCbwOjAfeAmYCuwOvm9mIuo4rYV0IYWzxX63FMUAIYVIIYUIIYYJ/oIGIiIiItBz+HPLW/tccyt1B/h/grBDCHfFGMzsauAg4uikzJiIiIiIthx6xKNg1hHBMfmMI4S4z+98mzJOIiIiItDDb0yMW5RbIaz7iPhERERHZhug1b6leZnZOie0G9GxoAiGEDh86VyIiIiLSougOcsE1QMc69l3bBHkRERERkRZqe7qDbCGE5s5DxoQJE8JLL71ETU0NABUVFck+z2tDGmjLli1A9rcd35Y/Pq6Dus7tx8bxK1asAEjyumjRoiTm7rvvBuC0004DYMOGDQDEb+lYsmQJAN27dwdg1apVANxyyy1JTLt27QAYNGgQAAMGDACga9euSUzv3r2BtK68DHFZ8nXneY7rp66yl9q+efNmoHSdrly5EoBZs2YB8O1vfzvzM6T18MEHHwDQpUsXAHbaaackpnXr1gCcfPLJmWNeffXVJObqq68G4OijC38zOm7cuFrnGT16NABvv/02AKNGjapVdi+P10urVoXfHdevX5/EtGnTBkjr+c033wRg5MiR1KWyshLI9h9P19P0tOIY5/Vaal/+PPFY8X3Lli0DoKqqKhMLad/K95eGjK9yY2b16tUAdOhQ938eleo369aty+Q1P84g7fcbN24E0j4SW7Om8ASY14e3W6k69DS8vuI+4fvyc1GpucC3+c+bNm1KYtq2bZspa/6YmPcFj/W6hLTtOnfunDkmPo/H+PHe/+qKj9OMt+fb1OsgjvFzf9wLptdHufm9rvzEvH3y5ynVV0tdH/L5KZW2H+f9z+sgHnsNkU+/rmsTQHV1NZCOp7hv+Rj2/Hg+4phy4zCvrj4fn9PzWmo8+Nzv4zQfC2mZG3Idzte3j+V4Xz7Wxz+kc4n3cU/D6xTS60qp8ZjPc7mx63OR10H79u1rlTOefwHWrl0LwAUXXJBsW7BgQeY83/jGNwB4/fXXk5iJEycCcOihhwIwb948ALp165bEzJkzB4COHQv3Oe+4o/C+hX322adWWscdd9zUEMKEUmWPDRo0KPzoRz+qL6xJnHnmmWXzaGYHA5cDFcC1IYSf1xG3O/A8cFwI4c5yaZb7qOnflDswhPD9cvtFREREZNvQUt+DbGYVwO+ALwDzgClmdl8IYXqJuEuAhxty3nKPWEz9iHkVERERkW1MC33EYiIwK4QwG8DMbgOOAKbn4r4H3EXhMz3qVe6jpv/40fIpIiIiItuaZlwg9zCzl6KfJ4UQJhW/7w/MjfbNA/aIDzaz/sBRwP583AWymd1X7sAQwpcbkoCIiIiIfLI18yMW1WWeQS61as//gd1lwI9CCDUNXeSXe8RiTwor8luBF+rIQL3MrAZ4Ldp0W10PT4uIiIhIy9RCH7GYBwyMfh4ALMjFTABuK+a/B3ComW0OIfylrpOWWyD3ofDA8/HA14G/AbeGEN74kBlfF0IY+yGPEREREZEWpCX+kR4wBRhuZkOB+cDXKKxbEyGEof69md0A3F9ucQzln0GuAR4CHjKzNhQWyk+Y2c9CCL/9qKUQERERkU+elngHOYSw2czOpPB2igrguhDCG2Z2enH/VR/lvOXuIFNcGH+JwuJ4CPAb4O4PmUaVmU2Lfr44hHB7Lp1TgVMhfd+viIiIiLQMLfmjpkMIDwAP5LaVXBiHEL7dkHOW+yO9PwJjgAeBC0MIr9cVW496H7Eo/iXiJCh8UMhHTEdEREREmkgLfcSiSZS7g/xvwBpgZ+D70W8NBoQQQqcmzpuIiIiItBAt9Q5yUyj3DPL282uCiIiIiNSppX6SXlMp+wxyI8k/g/xQCOG8rZCuiIiIiDQS3UFuRCGEiqZOQ0RERESalhbIIiIiIiJFesRCRERERCRne7qDbCG0rLeqTZgwIUyZMgXP1+bNm5N9FRWFpzV8X01NTbLPf6vZsmULAJWVlZnY+Htv4Pz5Yn4ej1m/fn2yz8/tMf417jhLly7NHLdu3ToA7rnnnlppeDnOO6/waPbJJ5+cxJx11lkA7LTTTgC0alX4naaqqqpW2Z3no1S5fJt/9TyUO0+pbX6cH1MqNp/GihUrkn3t27cH0vr1cjWkP8Z5Xr58eeb4Uu3u2rRpk4ktZdOmTZmvcT17W3oaXvY1a9bUKpfzGM8nQOfOnTNpeExc/94nfJv3n3bt2iUx5SaqfF/3rxs3bkxi8nWVj415Xr3uVq9enex76aWXANhrr70yMfF58mPWY9auXZvEeF3PnTsXgIsuugiAY489NokZM2YMAL17986UZ8mSJUlMz549M+fJxwJ07Ngxk5/WrVvXKvOH4fOUlzluy7rGYzy3eT48xvtqXIf58/jx3o4AGzZsyJTHx0qpebRUv3P5PpGfNwCqq6uBtD97PkrlOT9HPvXUU0nMrrvuCqR9rEOHDrXy07Zt28zx/nXBgvTTZPv06QPAypUrgXSsxPNFPo24XrxPeluUmkvy9VGqbfPjx9OPryGeN9+WnwdLnacUb4MePXpk8h6n5XXnabg4z57HxYsXZ87nP0Nav57Gm2++CcDo0aOTGO9LPj94Ob1fAsyYMQOAkSNHZvL34osvJjEjRozInKdXr16Z8wNceeWVAJx++umZ8sTjIX+dKiU/nvLjI1bqWu98W7wuic9XTql52a8ZXj/+Nc5zvnzev+N08/mKy+XHt2nTZmoIYUJ9+dxxxx3DJZdcUm95msJXv/rVBuWxMekOsoiIiIjUa3u6g6wFsoiIiIiU1ZI/Sa8paIEsIiIiIvXSH+mJiIiIiER0B1lEREREpEiveRMRERERydEdZBERERGRiBbIIiIiIiIRPWIhIiIiIlKk17yJiIiIiOToDrKIiIiISER3kEVEREREIlogi4iIiIgU6T3IIiIiIiI5uoPcjBYtWsQVV1zB3XffDUCrVmkWP/jgAwAqKysB+MIXvpDse+uttwAYPnw4AKeccgoA1113XRKzdu1aAHr06AHA/PnzAXjiiSeSmJ49ewLwzW9+E4D3338fgJkzZyYxnTt3BmDTpk0A/OUvfwHgs5/9bBLzwAMPADBy5MhMnhcuXJjErFmzBoAjjjgCgB/84AeZfAK88sormdibb74ZgGOOOaZWfrwcixYtyqQJcPjhhwNw7bXXAjB06FAAxo8fn8QccMABALzwwgsAfOYznwFg5cqVSYzXz+bNmzP10qlTpySmb9++AMyePRuAbt26AVBRUZHE1NTUEPPYdu3aJdu8rjzN1157DYA2bdrUKnvHjh0BGDJkCJC2TVwPt9xyCwBz584FoHXr1knM3//+dyDtGyNGjADS9o/P8+KLLwLw/e9/H4Af//jHSUyXLl0yZT7uuOMAuPrqq5OY9u3bA3DJJZcAcNlllwFw8MEHJzGPPfYYAF/60pcAWLFiBQA33XRTEuN1HkIAYMaMGck+L1vc3wAOPPDA5Pvq6mogbeeLL74YyLa3T4bTp08H4Je//CUA7777bhLjddSnTx8g7RveFgDdu3cH4OijjwZg1KhRAIwbNy6J2XHHHQGYNm1aJu1nn302iXnzzTcz5TvvvPMAuPfee5MYHz+eL6+fLVu2JDG+7bvf/S4Al19+OQC77LJLEtO1a9fMcevWrQNg1apVSYyPJx/D++67LwAbNmxIYvbcc08Afv/73wNpW+69995JjJf5jDPOAOCGG24A0nEPMHr0aCDt/95//vCHPyQxp59+OgB33HEHkI7zgw46KInxtvO+3rZt20zaABdddBGQ1o+3qccC9OrVK7Pv+eefB2Ds2LFJzOrVq4H0D3u8z/k8DfDqq69m8uV9xWPjeO/znmY8Pn0O8Lnlvffey6QdH+9tG5s3bx4Av/vd7wCYOnUqAJ/61KeSmMmTJwPw5JNPZo6Nr1NeL95ePqf59QvgzjvvBNJ55nvf+14m73EZfd70OcTPD7D//vsDaXv5vrj/+fzpfeknP/kJkO3HPv/5dffD8DkP0jnA0/Dx6XMepGO4Iby9fOzF7eb9w+dhHxelrjM+Dg455BAAHn/88STmsMMOA2Dx4sUA/Otf/6qVTx8r119/PZCuL3wNAWl7+fwzbNgwIL32Ajz00EMA/OxnPwPS8fTGG28kMd5PPH3Pu1//AL7yla8Q8/P6+Ae4/fbbATj00EMB+POf/wykbQTw9NNP82HpDrKIiIiISJFe8yYiIiIikqMFsoiIiIhIRI9YNAIzqwFeAyqBzcAfgctCCFvKHigiIiIiLY7uIDeOdSGEsQBm1gv4E9AZ+GkTpikiIiIijWx7e83bVilpCGERcCpwpm1Pv36IiIiIbCP8D/W29r8G5OtgM/uXmc0ys/NK7D/BzF4t/nvOzD5d3zm32jPIIYTZZrYD0Av4oL54EREREWk5WuI9TjOrAH4HfAGYB0wxs/tCCNOjsHeAfUMIy8zsEGASsEe5827tP9IrWbNmdiqFO8wl308pIiIiIs2nBT9iMRGYFUKYDWBmtwFHAMkCOYTwXBT/PDCgvpNutQWymQ0DaoBF+X0hhEkUVvMMGjQobK08iYiIiEjDNOMd5B5m9lL086Ti2hGgPzA32jeP8neHTwIerC/BrbJANrOewFXAFcE/wkpEREREPjGa8Q5ydQhhQh37Sq3aS641zWw/CgvkvepLsCkXyFVmNo30NW83Ab9uwvREREREpIm0xGeQKdwxHhj9PABYkA8ys08B1wKHhBCW1HfSJlsghxAq6o8SERERkZauBX/U9BRguJkNBeYDXwO+HgeY2SDgbuDfQghvNeSk+iQ9EREREalXS/wjvRDCZjM7E3gYqACuCyG8YWanF/dfBfwE6A5cWVzkby7zyAagBbKIiIiINEALvYNMCOEB4IHctqui708GTv4w57SW9jdz48ePDy+88ELy85Yt6SdTV1ZWAuB5jhsqXw7fF2/333w2bNgAQOvWrWvF5Bu/1Hk8T4sXLwagV69eAKxcuTKJWbp0KQC9e/cGYO3atQB069Ytidm4cWPJcsVl9jyvWLEic3xDOmlNTU3yfatW2d+FStXhh+HH+9f4t0rP//r164G0nuO0Nm3aBMDMmTMBGDZsGJDWU3weP+5Pf/oTAOPGjUtiBgwovKmlf//+ALRp0+ZjlauxeN63xm/b5dJat24dkLZ/HON1lO/jpdoyH7N58+Ykxvuv9zc/xrfH8u0S91H/3vPqfaFDhw5JjKfrXz3W+xpAVVVVrXPn0/Yyel5L5dnjfd+qVasyacbfz5o1C4A+ffrUWWbPl8f6uAA48cQTAbj44osBuP3224F0fACcd17h3fc777wzAPfeey8AY8eOTWJGjhyZKbt/9XkIao99r4t4jvO5yceT7/N5DaBLly4AVFRUZNLysR2nke8jcR/Lp+8/x+fJ96X8/BOnle+rcUx+Torz4emV6rf5vHqfKDVm8n281FybP97rpVT/8zz7datTp061YvLni+XrpdQ1Lb/N5422bdsmMT4e27dvnzkmHmf5PpUvb7zPx7D/HB/TkOvTu+++C8CgQYOAtP38Whmn4dfofB2UyqsfU6ouva/n+1G8z8eO7/OfIR1PrtTcnS+713u8vujevXsmTS97fP58/kuVx9uudevWU+u7mwowZsyYcMcdd9QX1iRGjx7doDw2Jt1BFhEREZF6NffNp61JC2QRERERqZcWyCIiIiIikZb4R3pNRQtkERERESmrBb/mrUlogSwiIiIi9dIdZBERERGRiO4gi4iIiIhEtEAWERERESkyMz1iISIiIiIS0x1kEREREZGIFsgiIiIiIkV6xEJEREREJEd3kEVEREREIrqDLCIiIiIS0R1kEREREZGi7e2jpi2E0Nx5yNh5553DFVdcwZgxYwDo2LFjsm/u3LkAdOjQAYBp06Yl+55++mkATj75ZACGDRsGwKZNm5KYyspKALZs2QJAmzZtasV4469evTqT/po1a5IYP27Dhg0AtG/fHoD33nsvibnlllsAOP744wH47W9/C8B5552XxDzxxBMA9O3bF4CBAwcCcO+99yYxJ554IgAbN27MpFlRUZHE9OzZE4DNmzeXLGe8La9U+3s99+7dG8j+l8q8efMAGDx4cOb4devWJTGzZ88G4NZbbwXgjjvuAGDp0qW18rN48WIgrWc/L8CvfvUrAPr37w+k9fTcc88lMZMnTwbglFNOyZy3Vav0d79u3bplyuHt7e0Iabt7efzntWvXJjErV67MnM/7iNdXnMcLLrgAgD333BOAl19+OYkZNGgQAJ/97GeBtG29zwIsWbIEgOHDhwMwZ86czM9xXj2PVVVVyT5v+7ifxGUHaNeuHfWpqanJpPXOO+8A2fr1Mvs2r5eFCxcmMd6u3gd69eoFwGmnnZbE7LPPPkA6hr0MXbt2TWJWrVoFpO0zdOhQIFv2c889F4CnnnoKgP322w9Ix0d8ngkTJgBpfe+9995JjM8z8TiC7FzgY7+u+o7zWm6ureuiE6eVby8/xtso3uZ9ytuie/fuSUw+r56vuJw+nn2MeMyiRYtqncfHQ3V1NZDWCUCXLl0AWL58OZD21Tg/F154IZC2iZfHx3Z8Tm8Tn/O8fADXX389AJdeeikAU6ZMAbJz7o9+9CMg7Xd9+vRJ9nk/GzFiBAArVqwAoG3btkmM59vrt3PnzgCsX78+ifE5yOcLHxetW7dOYnybn8f3xf3H29Lr2dOIx57XuW/r1KlTrbT8eG9DTyOe1z0ffh6Pjfusj+d+/fpl9sVzih/v16l8/4F0LvP5zssVt+Xrr78OpNeg++67D4CTTjopibnyyisBmD59OgA//elPAbj//vuTmPfffx9Ir63nn38+kB0z/n1+nMfXoquuugqA/fffH4Abb7wRgG9+85tJjOffz+d9I24Lb6/8NSnm9elt4nPbsmXLkhjvq/k2ja/DPuY8Pz5mSq0FKioqpoYQJtTakTN27Njw97//vb6wJtGzZ88G5bEx6Q6yiIiIiNRre7qDrAWyiIiIiJSl17yJiIiIiOToDnIjMrM+wGXA7sAG4F3gP0IIbzV12iIiIiLSOLRAbiRWqMl7gD+GEL5W3DYW6A1ogSwiIiLyCaBHLBrXfsCmEMJVviGEMK1MvIiIiIi0QLqD3HjGAFPrCzKzU4FTIX39k4iIiIi0HLqDvJWFECYBk6DwHuRmzo6IiIiI5OgOcuN5AzimidMQERERkSa0vX2SXlPfK38MaGNmp/gGM9vdzPZt4nRFREREpBHtsMMOzfKvOTTpHeQQQjCzo4DLzOw8YD3F17w1ZboiIiIi0ri2pzvITf4McghhAXBsU6cjIiIiIk1Dr3kTEREREcnRHWQRERERkcj2tEC2EFrWW9XMbDGwBqhu7rxIk+qB2nhbpzbePqidt31q423b4BBCz/qCzOwhCn2hOVSHEA7emgm2uAUygJm9FEKY0Nz5kKajNt72qY23D2rnbZ/aWLZH28/T1iIiIiIiDaAFsoiIiIhIpKUukCc1dwakyamNt31q4+2D2nnbpzaW7U6LfAZZRERERKS5tNQ7yCIiIiIizUILZBERERGRSLMukM1soJk9bmZvmtkbZnZWcfsFZjbfzKYV/x3anPmUj8fM2prZi2b2z2I7X1jc3s3MHjWzmcWvXZs7r/LRlGljjeVtjJlVmNkrZnZ/8WeN421MiTbWOJbtTrM+g2xmfYG+IYSXzawjMBU4EjgWWB1C+FWzZU4ajRU+eqd9CGG1mVUCzwBnAV8BloYQfm5m5wFdQwg/as68ykdTpo0PRmN5m2Jm5wATgE4hhMPM7BdoHG9TSrTxBWgcy3amWe8ghxAWhhBeLn6/CngT6N+ceZLGFwpWF3+sLP4LwBHAH4vb/0jhlyP5BCrTxrINMbMBwJeAa6PNGsfbkDraWGS702KeQTazIcA44IXipjPN7FUzu07/ZffJV/wvu2nAIuDREMILQO8QwkIo/LIE9GrOPMrHU0cbg8bytuQy5CS59QAABTVJREFU4IfAlmibxvG2pVQbg8axbGdaxALZzDoAdwH/EUJYCfwe2BEYCywELm3G7EkjCCHUhBDGAgOAiWY2prnzJI2rjjbWWN5GmNlhwKIQwtTmzos0jTJtrHEs251mXyAXn1e8C7glhHA3QAjhg+LFdgtwDTCxOfMojSeEsBx4gsKzqR8Un0P359EXNWPWpJHEbayxvE35HPBlM3sXuA3Y38xuRuN4W1KyjTWOZXvU3G+xMOAPwJshhF9H2/tGYUcBr2/tvEnjMbOeZtal+H0VcCAwA7gP+FYx7FvAvc2TQ/m46mpjjeVtRwjhv0IIA0IIQ4Cvwf9v795CvCjDOI5/f9mFQRZkGRLUhVYWVAYhhGUHouuKiqKMsqCgk0bdGEFCkUUUhdGN2UEiEDvSUVBQUiODUslCKCQ6UdZFSQdIny7mXRpDzTa31d3vB4adeeedd172zyzPPvvsvKyoqqvxOR4xdvcZ+xxrNDp4mO8/HZgJbGy1iwBzgSuTTKX7J58twI3DMz3tIxOBZ5OMofulbElVvZ5kLbAkyfXAF8BlwzlJ/Se7+4wX+yyPePPxOR7pHvI51mjjUtOSJElSz7DXIEuSJEn7EwNkSZIkqccAWZIkSeoxQJYkSZJ6DJAlSZKkHgNkSQeMJOOTfNS2b5N81TteM0T3PD3JwqEYezCSPJPk0j2cvyXJdf/nnCRppBnu9yBL0l6rqh/olrslyb3Atqp6eIhvOxe4b4jvsS8tAlYDTw/3RCTpQGUGWdKIkGRb+3pukpVJliTZnGR+kquSvJ9kY5JJrd9RSV5Msq5t03cx5jjg1Kpa347P6WWsP2znSXJXG2NDknm9669pbeuTLG5txyVZ3tqXJzm2tT+T5PEka5J8PpAlTmdBkk1J3gAm9Maf39o3JHkYoKp+AbYkcTlgSRokM8iSRqLTgJOAH4HPgYVVNS3J7cCtwGzgMeDRqnq3BanvtGv6zmDnZXXvBG6uqtVJDgV+S3IhcDwwDQjwWpIZwA/A3cD0qtqa5Ig2xgLguap6Nsks4HHgonZuInAWMIVuCeeldEv7ngicAhwNbAIWtfEuBqZUVQ0s9d18AJwNvD+Yb54kjXYGyJJGonVV9Q1Aks+AZa19I3Be278AODnJwDWHJRlXVT/3xpkIfN87Xg08kuR54KWq+rIFyBcCH7Y+h9IFzKcBS6tqK0BV/djOnwlc0vYXAw/1xn+lqnYAm5Ic3dpmAC9U1Xbg6yQrWvtPwG/AwpZZfr03znd0QbYkaRAssZA0Ev3e29/RO97BX4mBg4Azq2pq2475W3AM8CswduCgquYDNwCHAO8lmUKXNX6gN87kqnqqtddezLXfpz/v7KbPwFz+oMtav0iXgX67d3psm7skaRAMkCWNVsuAWwYOkkzdRZ9PgMm9PpOqamNVPUhXxjCFrjRjViu5IMkxSSYAy4HLk4xv7QMlFmuAK9r+VcC7/zDPVcAVScYkmUjLgLf7HV5Vb9KVjPTnfwI7l4ZIkv4FSywkjVa3AU8k2UD3s3AVcFO/Q1V9muTwXunF7CTnAdvpaoHfqqrfk5wErG3lGtuAq6vq4yT3AyuTbKcrwbi23XdRkrvoyjf+6ZVsLwPn05WHbAZWtvZxwKtJxtJlm+f0rpkOzEOSNCip2pu/AErS6JRkDvBzVe0370LekySnA3dU1czhnoskHagssZCkPXuSnWuD93dHAvcM9yQk6UBmBlmSJEnqMYMsSZIk9RggS5IkST0GyJIkSVKPAbIkSZLUY4AsSZIk9fwJZL6c2N0GXLUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x252 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAD0CAYAAACGjNCJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5hkZXXv8e9vRiJEUFQu4TaKSERDFJkB9cgJihrBJKLRKKBGNAY5kXg7iRqTc4STRNAYRJRIJoQIXiAmoCGKEKNRgoAygwgOKCGIMIDAgBdQIGFY54+9GzZFd3VPT1dXdfX38zz1dNe+vO+7a9fuWbVm1btTVUiSJElqLBn2ACRJkqRRYoAsSZIkdRggS5IkSR0GyJIkSVKHAbIkSZLUYYAsSZIkdRggS5pTSb6Q5LWz2O/xSSrJw+ZgDB9L8mcb2864S7IsyZ1Jlg57LJI0SgyQJW2wJNcmuasNrm5O8ndJNgeoqgOq6pR2u0OTnD/c0Y6n2XygaM/b8yeeV9V1VbV5Va0fzCglaWEyQJY0W79RVZsDewJ7AX8y5PFIkjQnDJAlbZSqugH4ArA7QJKvJHlDkicDJwLPajPNP2rXb5bkL5N8P8mPk5yfZLNOk69Kcl2SdUn+uF/fSfZJckGSHyW5Psmhk2zz6CSfS3Jrkh+2v+/YWX9okmuS3JHke0le1S5/YpKvtmNcl+Tv+4zjH5L8oN32vCS/1Fn3sSQnJPl828fXk+zSWV9JDk/yH+34TkiSdt2SJH/Svla3JDk1yaPaXc9rf/6ofX2flWSXJF9Ocls75k8m2bJt6+PAMuCf2+3f0ZuFTrJ9krOS3J7k6iS/2xnnkUk+3Y7hjiRrkqzorH9nkhvadd9N8rx+506SRpkBsqSNkmQn4EXAN7vLq+pK4HDgwva/8bdsV30AWA78D+AxwDuA+zq77gM8CXge8H/bQHuyfpfRBOYfBrYG9gAunWTTJcDfAY+jCRDvAj7StvEI4HjggKraoh3TRBt/CvwL8Ghgx7afqXwB2BXYBrgE+GTP+oOBo9q2rgb+vGf9r9Nk4Z8GvAJ4Ybv80PbxXOAJwOYTYwd+pf25Zfv6XggEOBrYHngysBNwJEBVvQa4jjbzX1Xvn+Q4TgPWtvu/HHhvT6D7YuB0YEvgLB54HZ8EHAHs1b6OLwSunaR9SVoQDJAlzdZn26zw+cBXgfdOt0OSJcDrgbdU1Q1Vtb6qLqiqezqbHVVVd1XVt4Bv0QSNk3kV8K9VdVpV/XdV3VZVDwmQ2+VnVNXPquoOmuB0384m9wG7J9msqm6qqjXt8v+mCaq3r6q7q2rKWuqqOrmq7miP40jgaZ1ML8CZVfWNqrqXJnjeo6eJY6rqR1V1HfBvnfWvAo6tqmuq6k7gj4CDpqo7rqqrq+qLVXVPVd0KHNtzrFNqP+jsA7yzPd5LgZOA13Q2O7+qzm5rlj/OA+dmPfBw4ClJNqmqa6vqP2fSrySNIgNkSbP1kqrasqoeV1W/V1V3zWCfrYBNgX7B0w86v/+MJmtKWxYw8VhGkx2dNghL8vNJ/rotU/gJTWnClkmWVtVPgVfSZLpvassgdmt3fQdNRvYbbTnB66dof2mSY5L8Z9v+tZ1j7XtMM1i/PfD9zrrvAw8Dtp1iLNskOb0tdfgJ8ImecfSzPXB7+yGi298Ofca5aZKHVdXVwFtpPhzc0o5h+xn2K0kjxwBZ0iBVz/N1wN3ALpNs27+hpixg4nEdcP0M2/nfNCUbz6iqR/JAaULads+tqhcA2wHfAf6mXf6DqvrdqtoeeCPwV0meOEn7hwAHAs8HHgU8vtv+RrqRJos9YRlwL3AzD31toSmvKOCp7bG+umcck+3T7esxSbbo6e+GmQy0qj5VVfu04y3gfTPZT5JGkQGypEG6Gdgxyc8BVNV9wMnAse0Xwpa2Xy57+Cza/iTw/CSvSPKwJI9N0lu6ALAFTd3xj5I8BnjPxIok2yZ5cVuLfA9wJ025AEl+q/Nlvh/SBH2TTYe2RbvvbcDPM4NSkw1wGvC2JDunmUbvvcDft6Uat9KUhzyhZyx30hzrDsAf9rR3c8/296uq64ELgKOTbJrkqcDv8NB66odI8qQk+7Xn8W6a19up4yQtWAbIkgbpy8Aa4AdJ1rXL/gC4HLgYuJ0m07jBf4vaLPKLaDLEt9N8uW6yeuXjgM1ostcXAed01i1p97+xbWNf4PfadXsBX09yJ80X0t5SVd+bpP1TaUoRbgCuaPuYKyfT1PqeB3yPJvj8fYCq+hlNPfXX0szi8UyaLwLuCfwY+DxwZk97RwN/0m7/B5P0dzBNBvxG4DPAe6rqizMY58OBY2he4x/QfFnx3TM/TEkaLanq9z9ukiRJ0uJiBlmSJEnqMECWJEnSgpXk5PZmSt+eYn2SHN/eAOmyJHtO16YBsiRJkhayjwH791l/AM3NnHYFDgM+Ol2DBsiSJElasKrqPJovWk/lQODUalxEMxf+dv3anPRuTMOUZN6/Nbh8+fL57nIkrV69eiDtDvv1HdRxaXxsyHt0XN5Pvcc8Lsc1Exv7N2khv1b9jn0hH9dMzPa8j8LrMsh/R1evXr2uqraebrthxGcda2hm8ZmwsqpWbsD+O9DMnT9hbbvspql2GLlZLIZxAkbtNRiWZC7ua/BQw359B3VcGh8b8h4dl/dT7zGPy3HNxMb+TVrIr1W/Y1/IxzUTsz3vo/C6DPLf0SSrq2rFDLYb5j/m044xyeOBz1XV7pOs+zxwdFWd3z7/EvCOqpry08/IZZAlSZI0epYsGU5l7n333bexTawFduo835FmvvcpWYMsSZKkaSUZymMOnAX8djubxTOBH1fVlOUVYAZZkiRJ00gytAzy+vX971yf5DTgOcBWSdYC7wE2AaiqE4Gzae68ejXwM+B10/U50AA5yXqaW8oGWA8cUVUXDLJPSZIkzb1RqMeeTFUdPM36At60IW0OOoN8V1XtAZDkhcDRwL4D7lOSJElzbFQD5EGYzxKLRwI/nMf+JEmSNEeGVWIxDIMOkDdLcimwKbAdsN9kGyU5jObOJpIkSRoxc/iFuQVhPkssngWcmmT36pnQr53seWW7nZMSS5IkjRgzyANQVRcm2QrYGrhlvvqVJEnSxjODPABJdgOWArfNV5+SJEmaGwbIc2eiBhmaqd5eW1X9J7OTJEnSSBnmPMjDMNAAuaqWDrJ9SZIkzQ8zyJIkSVKHAfIQLV++nFWrVg17GItSz+QiY2Ncj0vDMa7vp3E9rkEY19dqXI9rY/m6NCyxkCRJknqYQZYkSZI6zCBLkiRJHWaQJUmSpJa3mpYkSZJ6WGIhSZIkdZhBliRJklpO8yZJkiT1MIMsSZIkdRggS5IkSR2WWEiSJEktp3mTJEmSephBliRJkjrMIEuSJEkdBsiSJElSy3mQF6HF9Imon6oa9hCkeTWba39crpOZHPu4HGuvjf2bP66vy7ib7Xn3fD9gMcVLBsiSJEmalhlkSZIkqbXYpnlbPB8FJEmSNGsTQfJ8P2Ywrv2TfDfJ1UneNcn6RyX55yTfSrImyeuma3PgGeQkLwXOBJ5cVd8ZdH+SJEmae6NYYpFkKXAC8AJgLXBxkrOq6orOZm8Crqiq30iyNfDdJJ+sqv+aqt35ONKDgfOBg+ahL0mSJA3AiGaQ9waurqpr2oD3dODAnm0K2CJNY5sDtwP39mt0oBnkJJsDzwaeC5wFHDnI/iRJkjT3hjzN21ZJVnWer6yqle3vOwDXd9atBZ7Rs/9HaOLQG4EtgFdW1X39Ohx0icVLgHOq6qoktyfZs6ouGXCfkiRJmmND/JLeuqpaMcW6yQbVOzffC4FLgf2AXYAvJvn3qvrJVB0O+qPAwTSpbtqfB0+2UZLDkqxKsurWW28d8JAkSZK0oUa0xGItsFPn+Y40meKu1wFnVuNq4HvAbv0aHVgGOcljaSL13ZMUsBSoJO+onlm32zT5SoAVK1Y4I7ckSdIIGeE76V0M7JpkZ+AGmu+8HdKzzXXA84B/T7It8CTgmn6NDrLE4uXAqVX1xokFSb4K7AP8+wD7lSRJ0hwbxXmQq+reJEcA59IkY0+uqjVJDm/Xnwj8KfCxJJfTlGS8s6rW9Wt3kAHywcAxPcvOoInqDZAlSZIWkBHNIFNVZwNn9yw7sfP7jcCvbkibAwuQq+o5kyw7flD9SZIkaXBGMYM8KN5qWpIkSX3N9K5248IAWZIkSdMa1RKLQTBABnom1ZC0SCzma99j12Ljed94ZpAlSZKk1ghP8zYQBsiSJEmalhlkSZIkqcMAWZIkSeqwxEKSJElqOc2bJEmS1MMMsiRJktRhBlmSJEnqMECWJEmSWs6DLEmSJPUwgyxJkiR1mEGWJEmSWk7zJkmSJPUwQJYkSZI6LLEYotWrV8/7J5Sqmtf+RtVcve69r+dC+sQ5MfaFNOaFbCFee+P23hjGe34hnveucXgPTHYOBnVcnu+NNwqvoSUWkiRJUg8DZEmSJKnDEgtJkiSpwwyyJEmS1PJOenMoybbAB4FnAj8E/gt4f1V9ZpD9SpIkaW6ZQZ4DaV7FzwKnVNUh7bLHAS8eVJ+SJEkaDAPkubEf8F9VdeLEgqr6PvDhAfYpSZKkOWaJxSSS7AA8rrt9VZ03zW6/BFwyw/YPAw6bybaSJEmaf2aQO5K8D3glcAWwvl1cwHQBcm87JwD70GSV9+quq6qVwMp2u+HPhi1JkqQHMYP8YC8BnlRV92xg22uAl008qao3JdkKWLWB7UiSJGnIRjWDnGR/4EPAUuCkqjpmkm2eAxwHbAKsq6p9+7U5k48C17SNbagvA5sm+V+dZT8/i3YkSZI0RBO3mh7GY5pxLQVOAA4AngIcnOQpPdtsCfwV8OKq+iXgt6Y73ikzyEk+TFNK8TPg0iRfAu7PIlfVm/s1XFWV5CXAB5O8A7gV+CnwzukGJUmSpNEyoiUWewNXV9U1AElOBw6kKQ2ecAhwZlVdB1BVt0zXaL8Si4lSiNXAWbMZcVXdBBw0m30lSZI0OoZYYrFVkm6J7sr2+2sAOwDXd9atBZ7Rs/8vApsk+QqwBfChqjq1X4dTBshVdQpAkkcAd1fV+vb5UuDh0x+LJEmSxsGQp3lbV1Urplg3WdTeO+HDw4DlwPOAzYALk1xUVVdN1eFMjvRLbWMTNgP+dQb7SZIkaUyMYg0yTcZ4p87zHYEbJ9nmnKr6aVWto5mJ7Wn9Gp3JLBabVtWdE0+q6s4kA/uy3fLly1m1an4nupjJfxlUPfBhZFS/xbmxuse4ENodpIU45oVooV9L4/Q+mc9j8byPpkEdl+d7fIzoubwY2DXJzsANNKW9h/Rs80/AR5I8DPg5mhKMD/ZrdCYB8k+T7FlVlwAkWQ7ctYGDlyRJ0gI2il/Sq6p7kxwBnEszzdvJVbUmyeHt+hOr6sok5wCXAffRTAX37X7tziRAfivwD0km0tXb4RfvJEmSFo0ZljsMRVWdDZzds+zEnud/AfzFTNucSYB8GbAb8CSaQujvMLPaZUmSJI2JUcwgD8pMAuQLq2pP4P5UdJJLgD0HNipJkiSNlFHNIA9CvxuF/ALN3HKbJXk6D0yj8Ui8I54kSdKiYoDceCFwKM10Gcd2lt8BvHuAY5IkSdIIGfI8yPNuuhuFnJLkZVV1xjyOSZIkSSPGDDKQ5NVV9Qng8Une3ru+qo6dZDdJkiSNITPIjUe0Pzefj4FIkiRpNI3yNG+D0K/E4q/bn0fN33AkSZI0igyQO5I8AfgQ8EyggAuBt1XVNQMemyRJkkbEYiqxmMmRfgr4NM0d9LYH/gE4bZCDkiRJ0miZKLOY78cwzORGIamqj3eef6K957UkSZIWAad5e6h/S/Iu4HSaEotXAp9P8hiAqrp9gOOTJEnSCLAG+cFe2f58Y8/y19MEzE+Y0xENQVUNdPuFYq7e+OP6+mju+V5ZnDzvi4vne3wYIHdU1c7zMRBJkiSNJkssJEmSpB5mkCVJkqQOM8iSJElShxlkIMme/Xasqktm0kGS9cDlnUUvqaprZzQ6SZIkDZ23mn7AX/ZZV8B+M+zjrqraY+ZDkiRJ0qixxAKoqufO50AkSZI0uswgt5I8FjgE2K1ddCXwqQ28OchmSS5tf/9eVb10kn4OAw4DWLZs2QY0LUmSpEFbbNO8TXmkSZ4MfBtYDlwF/AewF/DtJLtNtd8k7qqqPdrHQ4JjgKpaWVUrqmrF1ltvvQFNS5IkaT5M1CHP92MY+mWQ/xR4S1V9urswycuAPwdeNsiBSZIkaXRYYtH45ap6ee/CqjojyXsHOCZJkiSNmMVUYtEvQP7pLNdJkiRpjDjN2wO2SfL2SZYHmHGhcFVtvsGjkiRJ0kgxg9z4G2CLKdadNICxSJIkaUSZQQaq6qj5HIiGr6oG0u6wL6hBHZckSYvJsP89n0qS/YEPAUuBk6rqmCm22wu4CHhlVf1jvzb73Wr6+H47VtWbpx2xJEmSFrxRnQc5yVLgBOAFwFrg4iRnVdUVk2z3PuDcmbTbr8Ri9SzHKkmSpDEzohnkvYGrq+oagCSnAwcCV/Rs9/vAGTT39JhWvxKLU2Y3TkmSJI2bIQbIWyVZ1Xm+sqpWtr/vAFzfWbcWeEZ35yQ7AC8F9mNjA+QkZ/XbsapePJMOJEmStLANucRiXVWtmGLdZFF775ePjgPeWVXrZxrk9yuxeBZNRH4a8PUpBjCtJOuByzuLTp+qeFqSJEmjaURLLNYCO3We7wjc2LPNCuD0dvxbAS9Kcm9VfXaqRvsFyL9AU/B8MHAI8HngtKpas4EDv6uq9tjAfSRJkjRCRvFLesDFwK5JdgZuAA6iiVvvV1U7T/ye5GPA5/oFx9C/Bnk9cA5wTpKH0wTKX0ny/6rqw7M9CkmSJC08o5hBrqp7kxxBMzvFUuDkqlqT5PB2/YmzabdfBpk2MP41muD48cDxwJkb2MdmSS7tPD+6qv6+p5/DgMMAli1btoHNS5IkaZBG+VbTVXU2cHbPskkD46o6dCZt9vuS3inA7sAXgKOq6tszHumDTVti0X4TcSXAihUrvKuDJEnSiBnREouB6JdBfg3wU+AXgTd3PjUEqKp65IDHJkmSpBExqhnkQehXg7x4PiZIkiRpSqN6J71B6VuDPEd6a5DPqap3zUO/kiRJmiNmkOdQVS0ddB+SJEkaLANkSZIkqWWJhSRJktTDDLI0h6qcuU+SpIXODLIkSZLUYQZZkiRJao3ynfQGwQBZkiRJ07LEQpIkSeowgyxJkiS1nOZNkiRJ6mEGWZIkSeowQJYkSZI6LLGQJEmSWk7zJkmSJPUwgyxJkiR1mEGWJEmSOgyQJUmSpJbzIEuSJEk9zCAP0erVq0fyBFTV/b+P4vjmQvcYx8mGnK+J12Bcz7HUy/f84jTZ33vfA6NrVP59NoMsSZIktZzmTZIkSephgCxJkiR1WGIxB5KsBy4HNgHuBU4Bjquq+wbVpyRJkgbDDPLcuKuq9gBIsg3wKeBRwHsG2KckSZLm2GKb5m1ejrSqbgEOA47IYvr4IUmSNCYmvqg3348ZjGv/JN9NcnWSd02y/lVJLmsfFyR52nRtzlsNclVdk2QJsA1w83z1K0mSpI03ijnOJEuBE4AXAGuBi5OcVVVXdDb7HrBvVf0wyQHASuAZ/dqd7y/pTfrKJjmMJsMsSZKkETPCJRZ7A1dX1TUASU4HDgTuD5Cr6oLO9hcBO07X6LwFyEmeAKwHbuldV1UraaJ5kozGbNiSJEm63xAzyFslWdV5vrKNHQF2AK7vrFtL/+zw7wBfmK7DeQmQk2wNnAh8pEbldjCSJEmasSFmkNdV1Yop1k0WtU8aayZ5Lk2AvM90HQ4yQN4syaU8MM3bx4FjB9ifJEmSBmQUa5BpMsY7dZ7vCNzYu1GSpwInAQdU1W3TNTqwALmqlg6qbUmSJM2fEb7V9MXArkl2Bm4ADgIO6W6QZBlwJvCaqrpqJo16Jz1JkiRNaxS/pFdV9yY5AjgXWAqcXFVrkhzerj8R+L/AY4G/aoP8e/uUbAAGyJIkSZqBEc0gU1VnA2f3LDux8/sbgDdsSJsjFyAvX76cVatWTb/hEPk9w4VlNufLc6zFxve8fA+onxGe5m0gRi5AliRJ0ugZ1QzyIBggS5IkaVoGyJIkSVKHJRaSJElSa4SneRsIA2RJkiRNywBZkiRJ6rDEQpIkSeowgyxJkiS1nAdZkiRJ6mEGWZIkSeowQJYkSZJallhIkiRJPcwgS5IkSR1mkCVJkqQOM8iSJElSy1tNS3Ns2BdUVQ21f01t2O+NjeV7a3Y87xqG2b7vPN8PsMRCkiRJ6ljoH243hAGyJEmS+nKaN0mSJKmHGeQ5lOQXgOOAvYB7gGuBt1bVVYPuW5IkSXPDAHmOpHklPwOcUlUHtcv2ALYFDJAlSZIWAEss5tZzgf+uqhMnFlTVpQPuU5IkSXPMDPLc2R1YPd1GSQ4DDgNYtmzZgIckSZKkDWUGeZ5V1UpgJcCKFSuccFCSJGnEmEGeO2uAlw+4D0mSJA3QYruT3qBz5V8GHp7kdycWJNkryb4D7leSJElzaMmSJUN5DMNAM8hVVUleChyX5F3A3bTTvA2yX0mSJM2txZRBHngNclXdCLxi0P1IkiRpMEZ5mrck+wMfApYCJ1XVMT3r065/EfAz4NCquqRfmyPxJT1JkiSNtlHMICdZCpwAvABYC1yc5KyquqKz2QHAru3jGcBH259TMkCWJEnStEYxQAb2Bq6uqmsAkpwOHAh0A+QDgVOrqoCLkmyZZLuqummqRkcuQF69evW6JD8F1g17LBqorZinczyiF/RiMG/neFh8bwGL4Dz3WoTnfdGd465FcL4fN5ONVq9efe6SJUu2GvRgprBpklWd5yvbKYIBdgCu76xby0Ozw5NtswOwcALkqto6yaqqWjHssWhwPMfjz3O8OHiex5/nWABVtf+wxzCFyT7B9N5TYybbPMhoVltLkiRJ01sL7NR5viNw4yy2eRADZEmSJC1UFwO7Jtk5yc8BBwFn9WxzFvDbaTwT+HG/+mMYwRKL1srpN9EC5zkef57jxcHzPP48xxpZVXVvkiOAc2mmeTu5qtYkObxdfyJwNs0Ub1fTTPP2uunaTfOFPkmSJElgiYUkSZL0IAbIkiRJUsdQA+QkOyX5tyRXJlmT5C3t8iOT3JDk0vbxomGOUxsnyaZJvpHkW+15Pqpd/pgkX0zyH+3PRw97rJqdPufYa3nMJFma5JtJPtc+9zoeM5OcY69jLTpDrUFOsh2wXVVdkmQLYDXwEuAVwJ1V9YGhDU5zpr0H+iOq6s4kmwDnA28BfhO4vaqOSfIu4NFV9c5hjlWz0+cc74/X8lhJ8nZgBfDIqvr1JO/H63isTHKOj8TrWIvMUDPIVXVTVV3S/n4HcCXNnU00RqpxZ/t0k/ZRNLd+PKVdfgrNhyMtQH3OscZIkh2BXwNO6iz2Oh4jU5xjadEZmRrkJI8Hng58vV10RJLLkpzsf9ktfO1/2V0K3AJ8saq+Dmw7MQ9h+3ObYY5RG2eKcwxey+PkOOAdwH2dZV7H42Wycwxex1pkRiJATrI5cAbw1qr6CfBRYBdgD5r7ZP/lEIenOVBV66tqD5q71+ydZPdhj0lza4pz7LU8JpL8OnBLVa0e9lg0GH3OsdexFp2hB8htveIZwCer6kyAqrq5/cf2PuBvgL2HOUbNnar6EfAVmtrUm9s69Il69FuGODTNke459loeK88GXpzkWuB0YL8kn8DreJxMeo69jrUYDXsWiwB/C1xZVcd2lm/X2eylwLfne2yaO0m2TrJl+/tmwPOB79Dc+vG17WavBf5pOCPUxprqHHstj4+q+qOq2rGqHk9zK9cvV9Wr8ToeG1OdY69jLUbDvtX0s4HXAJe3tYsA7wYOTrIHzZd8rgXeOJzhaY5sB5ySZCnNh7JPV9XnklwIfDrJ7wDXAb81zEFqo0x1jj/utTz2jsHreNy93+tYi423mpYkSZI6hl6DLEmSJI0SA2RJkiSpwwBZkiRJ6jBAliRJkjoMkCVJkqQOA2RJC0aSxya5tH38IMkNnecXDKjPpyc5aRBtz0aSjyV5eZ/1RyR53XyOSZLGzbDnQZakGauq22hud0uSI4E7q+oDA+723cCfDbiPuXQy8DXg74Y9EElaqMwgSxoLSe5sfz4nyVeTfDrJVUmOSfKqJN9IcnmSXdrttk5yRpKL28ezJ2lzC+CpVfWt9vm+nYz1N9v1JPnDto3LkhzV2f+322XfSvLxdtnjknypXf6lJMva5R9LcnySC5JcM5ElTuMjSa5I8nlgm077x7TLL0vyAYCq+hlwbRJvByxJs2QGWdI4ehrwZOB24BrgpKraO8lbgN8H3gp8CPhgVZ3fBqnntvt0reDBt9X9A+BNVfW1JJsDdyf5VWBXYG8gwFlJfgW4Dfhj4NlVtS7JY9o2PgKcWlWnJHk9cDzwknbddsA+wG40t3D+R5pb+z4J+GVgW+AK4OS2vZcCu1VVTdzqu7UK+J/AN2bz4knSYmeALGkcXVxVNwEk+U/gX9rllwPPbX9/PvCUJBP7PDLJFlV1R6ed7YBbO8+/Bhyb5JPAmVW1tg2QfxX4ZrvN5jQB89OAf6yqdQBVdXu7/lnAb7a/fxx4f6f9z1bVfcAVSbZtl/0KcFpVrQduTPLldvlPgLuBk9rM8uc67dxCE2RLkmbBEgtJ4+iezu/3dZ7fxwOJgSXAs6pqj/axQ09wDHAXsOnEk6o6BngDsBlwUZLdaLLGR3faeWJV/W27vGYw1u423XFnim0mxnIvTdb6DJoM9Dmd1Zu2Y5ckzYIBsqTF6l+AIyaeJNljkm2uBJ7Y2WaXqrq8qt5HU8awG01pxuvbkguS7JBkG+BLwCuSPLZdPlFicQFwUPv7q4DzpxnnecBBSZYm2Y42A97296iqOpumZKQ7/l/kwaUhkqQNYImFpMXqzcAJSS6j+Vt4HnB4d4Oq+k6SR3VKL96a5LnAeppa4C9U1T1Jngxc2JZr3Am8uqrWJPlz4KtJ1tOUYBza9ntykj+kKd+Ybkq2zwD70ZSHXAV8tV2+BfBPSTalyTa/rbPPs4GjkCTNSqpm8j+AkrQ4JXkbcNmV500AAABZSURBVEdVjcxcyP0keTrw9qp6zbDHIkkLlSUWktTfR3lwbfCo2wr4P8MehCQtZGaQJUmSpA4zyJIkSVKHAbIkSZLUYYAsSZIkdRggS5IkSR0GyJIkSVLH/wcVK+wVRwzlMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x252 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "start_sec = 25\n",
    "show_sec = 50\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [1, 0.05]}, figsize=(10, 3.5))\n",
    "im = libfmp.b.plot_matrix(predictions.T[:, int(start_sec*fs_hcqt):int(show_sec*fs_hcqt)], Fs=fs_hcqt, ax=ax, cmap='gray_r', ylabel='MIDI pitch')\n",
    "ax[0].set_title('Pitch-class predictions')\n",
    "ax[0].set_yticks([0, 2, 4, 5, 7, 9, 11])\n",
    "ax[0].set_yticklabels(['C', 'D', 'E', 'F', 'G', 'A', 'B'], fontsize=10)\n",
    "ax[0].set_xticklabels(np.arange(start_sec-5, show_sec+5, 5))\n",
    "ax[1].set_ylim([0, 1])\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "annot_folder = os.path.join(basepath, 'data', 'Schubert_Winterreise', 'ann_audio_note')\n",
    "fn_annot = os.path.join(annot_folder, fn_audio[:-4]+'.csv')\n",
    "\n",
    "df = pd.read_csv(fn_annot, sep=';', skiprows=1, header=None)\n",
    "note_events = df.to_numpy()[:, :3]\n",
    "        \n",
    "f_annot_pitchclass = compute_annotation_array_nooverlap(note_events.copy(), f_hcqt, fs_hcqt, annot_type='pitch_class', shorten=1.0)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [1, 0.05]}, figsize=(10, 3.5))\n",
    "cfig, cax, cim = libfmp.b.plot_matrix(f_annot_pitchclass[:, int(start_sec*fs_hcqt):int(show_sec*fs_hcqt)], ax=ax, Fs=fs_hcqt, cmap='gray_r', ylabel='MIDI pitch')\n",
    "ax[0].set_title('Pitch-class annotations')\n",
    "ax[0].set_yticks([0, 2, 4, 5, 7, 9, 11])\n",
    "ax[0].set_yticklabels(['C', 'D', 'E', 'F', 'G', 'A', 'B'], fontsize=10)\n",
    "ax[0].set_xticklabels(np.arange(start_sec-5, show_sec+5, 5))\n",
    "ax[1].set_ylim([0, 1])\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}