{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 3D Time/Frequency/Phase Representation of Audio for Speech Recognition.\n", "\n", "I recently participated in a speech recognition kaggle competition. Although I didn't come close to the top of the leaderboard (238th place with 87% accuracy vs 91% accuracy for the winners) I learned quite a bit about handling audio data and had a lot of fun. One of the more novel things I tried during the competition was to spatially encode the phase information in the audio and pass the results into a 3D CNN.\n", "\n", "A common pre-processing step in speech recognition is to turn the 1D audio into a 2D spectrogram. The spectrogram the volume of the audio as a function of time and at a particular frequency. Spectrograms are a great way of summarizing the important information in an audio clip in a way that makes it accessible visually. Here is a spectrogram of an utterance of the word \"marvin\".\n", "\n", "![marvin_specgram](../images/marvin_specgram.jpg)\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although the plot of the audio trace is interesting it is nothing compared to the richness of the visual information available from the spectrogram. Since spectrograms are not only beautiful but are also intuitively interpretable it is no wonder that they would be a popular choice for preprocessing audio data before digesting with a neural network.\n", "\n", "However not all information in the input audio makes it into the output spectrogram. To find the volume at a particular time and frequency we first decompose the audio into a sum of windowed sinusoids using a short time Fourier transform (STFT). The STFT gives us a set of complex numbers with the magnitude of those numbers represnting a magnitude of the associated sinusoid and the angle of the complex number representing its phase. To get the spectrogram you want just the volume information which is the magnitude of the STFT and you throw away the phase information entirely.\n", "\n", "If you look at the phase information corresponding to the same audio sample as above it is easy to see why the phase is usually just thrown away.\n", "\n", "![phase_angle_example](../images/marvin_phase_angle.jpg)\n", "\n", "You might be able to see some weak signs of structure in the raw phase data but you could be forgiven for assuming that it was all just structureless random noise. The fact that we can do good word classification without incorporating the phase information seems to back this idea up. However if the phase information really were just random noise which carried no additional information then any 2 sound clips which share the same spectrogram should sound the same to our ear. Lets listen to the example sound clip." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Using TensorFlow backend.\n" ] } ], "source": [ "import os\n", "import numpy as np\n", "import scipy.signal\n", "\n", "import IPython.display as ipd\n", "import tf_speech_recognition as tsr #local python library written for the competition" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "example_file = os.path.join(tsr.data_dir, \"raw\", \"train\", \"audio\", 'marvin/fda46b78_nohash_0.wav')\n", "audio = tsr.read_wav(example_file)\n", "ipd.Audio(audio, rate=16000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nice and understandable and human sounding. Now lets see what happens to the audio when we replace the phase information with random phases. We simply generate white noise in the time domain and then apply the same STFT to both noise and voice audio. But at reconstruction time we take the phase information from the noise and the magnitude information from the noise. The resulting audio sample will have exactly the same spectrogram as the original audio but truely randomized phases. Will it sound the same?" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stft = scipy.signal.stft(audio)[2]\n", "stft_mag = np.abs(stft)\n", "\n", "noise_stft = scipy.signal.stft(np.random.normal(size=audio.shape))[2]\n", "random_phase = noise_stft/(np.abs(noise_stft)+1e-5)\n", "phase_randomized_audio = scipy.signal.istft(stft_mag*random_phase)[1]\n", "ipd.Audio(phase_randomized_audio, rate=16000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sound clip is still understandable but now the speaker has seemingly been assimilated by the borg. Instead of sounding like a single coherent voice it suddenly sounds like a cacophony of many nearly coherent voices but all with tiny random offsets. Although the second clip is certainly still understandable it no longer sounds much like a normal human voice. That may be enough to show that that the phases are not random but it isn't enough to show that the phase information adds anything above and beyond the regular spectrogram when it comes to speech recognition. After all we can still make out what the word is supposed to be just fine even with totally random phases.\n", "\n", "A good way to address whether or not the phase information adds anything to the classification problem is to try and incorporate it into a machine learning system and see if we can get improved classification performance using phase information than without it.\n", "\n", "The first problem we have to solve is how to best process the audio to make the phase information available to our models. We can of course simply feed the raw audio in to our neural network but then we lose all the nice aspects of using spectrogram inputs. We can simply directly feed the phase in radians in to the network as an extra channel along with the magnitude but as we saw above the phase information as a raw phase angle doesn't have any apparent image like structure, it just looks like a random mess. \n", "\n", "An interesting idea is to try and map the audio into a 3D time/frequency/phase space just in the same way that audio is commonly mapped to a 2D time/frequency space. Just like frequency and time there is a natural sense in which one phase is closer to another, $\\pi$ radians is close to $7/8 \\pi$ radians and far from 0 radians. We can somewhat naturally then encode the phase information as an index along the third dimension of a 3D tensor in the same way that the STFT encodes frequency by an index along the second dimension of a 2D tensor. We can either then treat the resulting phase information as extra channels in a 2D input image to a 2D CNN or double down on the spatial representation and feed the result as a 1 channel 3D \"image\" into a 3D CNN. The 3D CNN idea is too juicy to ignore and so that is what I decided to do. I also experimented with passing the phase information as channels into a 2D CNN and the results seemed comparable but the idea just isn't as much fun. Unfortunately the 3D CNN is much more memory greedy and so takes longer to train and requires that smaller batch sizes be used.\n", "\n", "No implementation that I attempted \n", "\n", "## Base TF Phased Spectrogram Code.\n", "\n", "Here is some base code for processing the audio inputs and generating a \"phased\" spectrogram in which the magntidue of the spectrogram is modulated by how closely aligned it is to a target phase angle, one phase angle per slice in a 3D output tensor.\n", "\n", "The amount which the spectrogram gets attenuated is determined by ak \"phase_kernel\" which takes as an input the euclidean distance between the location on the unit circle of the target phase and the central phase of the corresponding index in the 3D tensor. \n", "\n", "The reason for taking the linear distance in the 2D plane is to handle the fact that 0.9$\\pi$ radians and $3.1\\pi$ radians are equally close to $\\pi$ radians since any difference of $2\\pi$ radians is the same as a difference of 0. Ideally we would take a \"distance\" between phase angles directly in radians since that would cause a minimum of distortion but I can't think of a way to express that in TensorFlow code that will be totally free of jump discontinuities. Just so long as our phase kernel drops off fast enough the 2D euclidean distance should be a servicable approximation to the phase difference in radians. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import time, copy\n", "\n", "import keras\n", "import keras.backend as K\n", "import tensorflow as tf\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "input_length = 16000\n", "target_length = 16384\n", "\n", "sound_in = keras.layers.Input(shape=(target_length,), name=\"audio\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def make_spectrogram_function(\n", " frame_length, \n", " frame_step,\n", " magnitude_activation=lambda x: tf.log(1+x),\n", " fft_length=None\n", "):\n", " \"helper function for making spectrogram Lambda layers in keras\"\n", " def _spectrogram(x):\n", " stft = tf.contrib.signal.stft(\n", " x,\n", " frame_length=frame_length,\n", " frame_step=frame_step,\n", " fft_length=fft_length,\n", " pad_end=True,\n", " )\n", " mag = tf.abs(stft)\n", " return magnitude_activation(mag) \n", " return _spectrogram\n", "\n", "def make_phased_spectrogram_function(\n", " frame_length, \n", " frame_step, \n", " phase_kernel,\n", " fft_length=None,\n", " phase_mag_damp=1e-5,\n", " phase_resolution=8,\n", " phase_only=False,\n", " magnitude_activation=lambda x: tf.log(1.0+x),\n", "):\n", " \"helper function for making phased spectrogram Lambda layers in keras\"\n", " def _spectrogram(x):\n", " stft = tf.contrib.signal.stft(\n", " x,\n", " frame_length=frame_length,\n", " frame_step=frame_step,\n", " fft_length=fft_length,\n", " pad_end=True,\n", " )\n", " mag = tf.abs(stft)\n", " real_part = tf.real(stft)\n", " imag_part = tf.imag(stft)\n", " xpos = real_part/(mag+phase_mag_damp)\n", " ypos = imag_part/(mag+phase_mag_damp)\n", " \n", " #generate phase_resolution equally spaced angles\n", " phi_knots = np.linspace(0, 2*np.pi, phase_resolution+1)[:-1]\n", " x_knots = np.cos(phi_knots)\n", " y_knots = np.sin(phi_knots)\n", " \n", " mag_act = magnitude_activation(mag)\n", " \n", " if phase_only:\n", " activations = []\n", " else:\n", " activations = [mag_act]\n", " for i in range(phase_resolution):\n", " dist_sq = (xpos-x_knots[i])**2 + (ypos-y_knots[i])**2\n", " phase_act = phase_kernel(dist_sq)\n", " if phase_only:\n", " act = mag_act*phase_act\n", " activations.append(act)\n", " \n", " activations = tf.stack(activations, axis=3)\n", " \n", " return activations \n", " return _spectrogram" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Picking a phase kernel\n", "\n", "To encode the phase information spatially we need to pick both a number of phase slices to use and how rapidly the activation of a particular phase slice should fall off as the phase of the input spectrogram becomes more distant. Since we are using the 2D euclidean distance as an approximation it makes sense to simply plot up the value of a proposed kernel in 2D. We want to pick enough points on the unit circle so that we can get a good approximation to the phase angle but not so many that we will run out of memory resources when building our network. We also need to make sure that the kernel function we choose drops off fast enough so that the euclidean distance approximation is still roughly valid (we don't want $-\\pi$ and $\\pi$ to start both simultaneously having high non-zero activations) but we don't want it to drop off so fast that a phase half way in between two slices in our phase grid wont significantly activate its two neighboring slices. The ever ubiquitous Gaussian kernel seems like as good as any to use and I tuned the kernel width and the number of phase slices by tweaking the look of the plot below. The larger the number of phase slices the better resolution in phase but also the more redundant our representation and the less well the 2D euclidean distance between adjacent phases approximates the phase difference in radians." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEzCAYAAACc82y3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsvXu8JFV1t/+s7jPDOCACgoSAXEyQqCReIBjj71VQUTRGjLdgYgRFiUaTGKJvJBrkhxpR4y2RJE4EUZOISrwMCiLogCaCAYxRwCCICiMEuQsODDOn1/vHrupe1bV3V1V3nT7dc9bDpz+neu9du3Z3D717fWtdRFVxHMdxnOWks9wLcBzHcRzfjBzHcZxlxzcjx3EcZ9nxzchxHMdZdnwzchzHcZYd34wcx3GcZaeVzUhETheRn4rIFYl+EZG/FZFrReQ7IvI403e0iFyTPY5uYz2O4zjOfNGWZXQGcMSI/mcC+2eP44B/ABCRXYC3AI8HDgHeIiI7t7Qmx3EcZ05oZTNS1a8Bt48YciTwMQ1cAuwkInsAzwDOV9XbVfUO4HxGb2qO4zjONsi07hntCdxgnm/M2lLtjuM4zgpiYUrXkUibjmgvTyByHEHio0v3oLXs2N7qnNki9q+incGjz5DmcxUnqzi/cvrEgHHOS7/IsZYw3rUaMM57PyY/u/emW1V1t7bnfcZh2+ttty82Pu/y72w+T1VXvCI0rc1oI/BQ83wv4Mas/dCh9gtjE6jqOmAdwI6yiz5enroU63QsbX5BSH0jXDoNrhubt+J8Sb2uTmQuOzbSX5grdl27vtS68jkSY6PrTa2rP1diXVVzmWPtROaqXEvVa4z3a6y9iW7T4N/ql7/91h83mLk2t92+yH+et3fj87p7XLPrEixn7piWTLceeGnmVfcbwF2qehNwHvB0Edk5c1x4etbmOI4zVyjQG+M/J9CKZSQinyBYOLuKyEaCh9wqAFX9R+Ac4FnAtcAm4GVZ3+0i8lbg0myqk1V1lCOEs9SsZGso+ut8DGsIBmtMWSUVr0Fi1khqrbH1NLBmNLnGyaysgrUT+6eQ+kwi7VHLyTITEZPKovrmMi6tbEaq+uKKfgVek+g7HTi9jXU4S8NhvR9zLFewG5u4hbWcxoFs6Oyz3MuaWw697xpe9vP/ZLfePdzS2YGPbP94LnzAw5d7Wc6EBMvIS/KMy7TuGTlzymG9H3M8l7OGcGN2dzZxPJdDj/iG1MAaakRq3nEsojGsocJcnQoLx44ZsoYOvff7vO7ur7GGrQDs3ruH1919EdLpcOH2B4xeQ5U1ZNvHvQ/UwMrqWyt2qRUWWaW1lDivymLX6fk/jMRlt/GZCePWmV2O5Yr+RpSzhkWOJZpsw6ngZT//Zn8jylnDVo65++JlWpHTFoqyqM0fTsAtI2cku7GpUbszmt169zRqd+YLl+nGxzcjZ6QEcgtr2T2y8dzC2skvW+W0ME1pzoxJOzuMdsEeWkR07C2dHdg9svHc0n1gtYNClaOBOS44JcTmrTFXXRmukfRWmD8xNvKeRh0YZkSay1Fg0TejsXGZzhnJaRzIfXQLbffR5TQOXKYVzTcf2eE3uG/oN+B9ssAZD/yNZVqR0yY9tPHDCbhl5IxkQ2cf6OHedC1x4drgNfeyey7pe9OdseMTuHDtARVnOrOOgt8DmgDfjFYyNWOKNnT2YQMVm0+bMUVjeORVxhHViDOKy3sRz7mqWJvE2Hz+i7Y/gIt2+JXyGjoJ+a+uh5xtbzBXUmbrS36dkf1NpLfCtezYqOddYq4msXBT1n7cl258fDNyHMdpAUX9ntEE+GbkOI7TBgqLvheNjW9GzvhMU5pLnD9pUGttac62p9L6xNL5VCU0tccVHm7arZEOKD9OjI16yJl1RSW7lDzYrS+9RSW3hDxY13OuMkXQlAkZGJxx8c3IcRynFYTFWfM3nyN8M1ppTPprch6soehcFdaQvV4qjqhJqYcmyU27ESeLmINCaq5ITFHSKaFmfxhTnl8j19WENRS1diTSNtReuv4w0TikxNgpokDPZbqx8c3IcRynJdwyGh8PenUcx3GWHbeMVgIzdqO3kogM02Z11uS16tYgqpLm7HVjclxp3poOCimnhFgcUFV/TI4rrSE2V1mGS8ltsfaCnFaQ/0x7P80QcaLOEImxUySkA5qBhcwpvhk5juO0RG8Wbl7NKb4ZOY7jtIBbRpPhm5FTTdvlw2uWDW9TmkuWCh8jjqhSmrPHKcnPymSx68bS+aTigboxyS9xrVwui8lxEE0tVJgrIsNpwRuvgfSWOI553lk09s9xBvYARVj02/Bj45uR4zhOS7hMNz6+GW2rtOG0sMzJTwukHBFi8TyF7ry/whqyc42bVSFmzVRZQ7a9wtkhbe2UraTi/GUHhGhy1eFrCKWxRWeH8mssWC0V1k6jsWM6LUxzb3CZbjJ8M3Icx2kFYTGqITp18M3IcRynBUJuOt+MxsU3I2f5aJLupzigojsiUVVJc3ZMkxQ/qbpBEektWQo8kg6oEBsUk8NSkl0seWlMpks5JXTK19BEf2WckSkQXCW9xZ0SUs4Q5aGzcqvGZbrx8c3IcRynBVRdppsE34wcx3FaoueW0di0shmJyBHAB4Au8GFVPWWo/33AYdnTtcBDVHWnrG8R+G7Wd72qPqeNNa1Y5iErd92YoqpUPmZMrTiiGOPEEVVJc+a8ZIqfblnqS8UO9SW7QlsDb7nIeTE5LqyB0thC9u3CtdJtQFSSK3rxlftte7U0l/j3t0z7QfCmc8toXCbejESkC5wKHA5sBC4VkfWqelU+RlX/zIz/Y+CxZop7VfUxk67DcRxneXGZbhLaeOcOAa5V1etU9X7gTODIEeNfDHyihes6juPMDLk3XdOHE2hDptsTuME83wg8PjZQRPYB9gO+aprXiMhlwFbgFFX9XAtrcpaLSYNboVG6n0rPu/45KeltjKDWCmkOjDxXIc2FsbHg0bI3XVFai0uFMW+5qLSWlPwotRe84iIBsJqQ8WIBrinpLRpYm/hoq/qjY6fE4qy49c0hbWxGsXc/Ve/wKOAsVV00bXur6o0i8jDgqyLyXVX9QekiIscBxwGsYe2ka3Ycx2kVz003GW1sRhuBh5rnewE3JsYeBbzGNqjqjdnf60TkQsL9pNJmpKrrgHUAO8ouXtzXMktOCymaJEIdDDDnj1GXyLankp/GrldV3jtVHjzmrFBlDZn2ytihiAUUrtHAQaHK2om0pyyfuGUUH0tsbMpKilk+Vf1E+p25oo1t/FJgfxHZT0RWEzac9cODROQAYGfgYtO2s4hslx3vCjwRuGr4XMdxnHmgp53GDycwsWWkqltF5LXAeQTX7tNV9UoRORm4TFXzjenFwJmqaq2aRwAfEpEeYWM8xXrhOY7jzAvu2j0ZrcQZqeo5wDlDbScOPT8pct43gF9tYw3O0tFmjSKoEVO0FOl+Upm4Yw4KVaXCCxJaRRxRhTRXaE+MjUpvC/EM4FEHhSYynZXR+vIg8f5YnFHEaaFwnOivkuGiKYdmTKZTxB0YJsAzMDiO47SEu2qPj29GjuM4LaCKB71OgG9G88qUi+ctCZUecol4oso0RBHPuVQm7pg8WJEOKFWYLprip0qaS4y18l9l7NBCc2+5XmGu0WOT/Z1Yf8VxIrYonjqIaP+oc5YX8dx0E+CbkeM4TgsobhlNgm9GjuM4LeHedOPjm5GTZKys3EtUMK9yDam1VmXijh0nPOQGWawjhfGG2qs85KJBrRFprtgf95DrFebN2iLSXTiv3NYrSHrmuMqbLirj2eOIFFjHm64idVA0qzfl/mmjCL3Z0AvnEt/GHcdxWmKRTuPHvCIivykizzbPHywinxCR74rI32QVHWrjltG8McXUP0tGTceFWjWK6qb7qXJaSMylMcsnlQ4o5uxQw4Ehv0aVU0La2rHnlfsLlk/E2kk6M8Qsn4gVlU4nFDmuikMyY6osn6p6RyPblwCFlZZR4RTgK8AXsufvBp4FXAC8GrgLeGvdyVbUO+c4jrN0CItjPOaYRwCXAYjIKuAFwJ+p6vOBNwG/12Qyt4wcx3FaYAVaRjsAP8uODwG2Z2AlfQvYu8lkvhk5BZa8lHhxgDl/zP+Jq9L91M3Ebcc0cHCI1R2CGil+YnFECRkul+fScUaDJfQizg4xB4Vkf3SsbauIQ6qQ6VLf1YX2cRwYZDYS+c+5pdOUnwCPBr4OPBO4QlV/mvXtDGxqMplvRo7jOC2gKivNMvoE8NcicijhXtFbTN/jgGuaTLai3jnHcZylZFE7jR9ViMgRInK1iFwrIm+M9O8tIhtE5L9E5Dsi8qwleXFlTgLeCWxHcGZ4r+l7NPDpJpO5ZbQSmAcPOsNAWqvwoCvMn0j3U5XiJ3Ic9aCDvjyXlOYiUmC6/HfZcy4mzdk5CtLaQlkus3OkpLeBjJfqt+2ROKPIecn+mOdcjdRBueTWJKt3gW1IKcvco08FDicUMr1URNYPldp5M/ApVf0HEXkkoYLCvku9tqxi99sTfc9tOp9vRhEO0+s5livYjU3cwlpO40A2SKN7cY7j1OTZ37mc479yDnvcdQc3PWhn3nP4M/nCow9a7mU1RmEpctMdAlyrqtcBiMiZwJEUi5AqsGN2/CDSlbZnGt+MhjhMr+d4LmcNiwDsziaO53JQlm9DaiMp6qjpp+m0EAZl51ckQq2qUZRaW8wxIhmzFLFWqsYWftGnrKSytVNwWoiMjVlD9jhlDRXihGKJUCOWT8oyijk2pMdmfyXRH3NK6Gih/9nfvpy3nv1p1m7ZAsCed93B2z7/abQL6x93UOn8wvVSTgvLZhnJuLnpdhWRy8zzdaq6LjveE7jB9G0EHj90/knAl0XkjwkebU8bZxFNySp7n0AonLo3Qa6zqKrW3mN8MxriWK7ob0Q5a1jkWK5gQzNPRcdxKnj9+ef2N6KctVu28Przzi1uRnNAcO0eaye8VVUPTvTFJhzehV8MnKGq7xGRJwAfF5EDVbU3zmIa8G7gNcC5wGeAzZNM5pvRELslvBFT7Y7jjM8ed97RqH3WWYL0PhuBh5rne1GW4Y4FjgBQ1YtFZA2wK/BTlpYXAG9R1eh9o6b4ZjTELaxl98jGcwtrl2E1EzBnTgvjXyP8cKysURSLJ7L99jiV7icSO1SU7CpS/KQSncbiiKL9g2ulZbhc0iPaH5XxEmOjMl2jdEFabh+S227aaWf2jGw8N+68E73VxgCYA5luiRKlXgrsLyL7EeJ6jqKc2eB64KnAGSLyCGANcEvbC4mwA3BxW5PNwDfWbHEaB3Ifxfx+99HlNA5cphU5zrbLu494JptWrSq0bVq1ir951jOXaUWT0aPT+DEKVd0KvBY4D/gewWvuShE5WUSekw37c+CVIvLfhNifY1R1GlHAZwNPamsyt4yG2CB7g+LedI4zBc7O7gu9/kvn8ot33sGNO+3M3/zWEaw/aL7uFwFZ2fH2zTJVPYfgrm3bTjTHVwFPbP3C1fwd8DER6RHWd/vwgNwLsA6+GUXYIHsvv7PCEnvQ1VvDEtUrqpuVe1wPulh7UnqLxBSl0v1EYoeK6X4i3nKpOKNKzztKx0lpznjZxaS1omRXT8YrXnd0f1SOI+FN19VS/+d/43F87jcfZ8YqoEMZvs0P/ao4oxhTShe0wuoZ5RLdSRSzL1hql5HwzchxHKcFwj2jFXXn4+WUPfvGxjcjx3GcllhJiVJV9Yw25/PNaAVTGew6DkvlQTdpVu5YcGuNsdEA15iH3dC66Eakt1RBu3xsJPu2Pa6S5uwcxbaIjFfhQWfbi9Lc4Idwr1/Ib9BfeL+6ZfmuILclxg5kuPjYfrv955tMDTS9bN4TxBk5tLQZicgRwAcI+uCHVfWUof5jCAFSP8maPqiqH876jibkVgJ4m6p+tI01OY7jTJcVJ9MhIg8hBN0eQHApt6iqHlt3rok3o5qJ/AA+qaqvHTp3F8KNr4MJPywuz86dz4i3WWDS+KIG57dar6iqlHjqunUToTZKlJpYd6zMdsVxMZ4nEZ8UdVAYneInZg3Z9oI1FLGC6jgw9C2fiIVjj9P9ZcumaAElLJ/cejJWlEQcGEQi1tQQ0/YDWoLcdDOLiBwAXEIwQrYHbgV2yZ7fQSg7Xps2tvF+Ij9VvR/IE/nV4RnA+ap6e7YBnU8WSew4jjNP5K7dTR9zzLuB/wR2J/wceCbwAOAVhMJ6v9NksjY2o1givz0j456f1do4S0Ty9BZ1z3Ucx5l5etpp/Jhjfh34ewY56TqqulVVTyfEIL2/yWRt3DOqk8jvbOATqrpZRF4FfBR4Ss1zw0VEjgOOA1gzb6l56jIFTWEsp4U650xar2jcNcSuG3FQqHRaIOGgEJEaYyXDQ3skdqjQxsjzig4O9rjswJB0OlgoxxHFJLlkHNJCWXJLxhllYwvSW1KmK0tvdOMyXH4sBZ8ULR8X2uJjp8kSpQOaZXYAblfVnojcRciHl3MZcGL8tDhtbMuVifxU9TZVzXfPfwIOqnuumWOdqh6sqgevKmUqdxzHcabMj4BfyI6vBl5o+p4N3NlksjY2o34iv6y+xVHAejtARPYwT59DyLEEId/S00VkZxHZGXh61uY4jjN39JDGjznmfILjGoSS4y/LyqNfCfwpcHqTySaW6VR1q4jkify6wOl5Ij/gMlVdD/xJltRvKyF/0THZubeLyFsJGxrAyapaym/kzBHjSo0xGS3l2bdU3nSdSFvMw62qOB9Weov3x1P4VKQDSkpvo9MBVc1lpbmopBeR7mx7VI6DggyXH4tp63Qj0hsgnVCGp2Pmssf52JQ0t3wy3YqLMzqBrKCeqn5KRO4FfhdYSwj1+acmk7USZ1Qjkd8JhIXHzj2dhjuo4zjOLDLnDgmNyG69bDbPzyb4B4yFZ2DYFliq2KJxy4r3z29QVrzi+sl6RZGxTUqJVyZCjTgzFMts27GmPb+GdVqwVkfEmSHmtGDbk+XBC1ZOpC3irNBbpaVzymO11EbBSsqOrQW0MCguGrOCOh3Tb6ydbnfQnltBdmw3YvlYC6iTsJIG/VOwlnTFOTC0im9GjuM4LaBs+0GvIvJV4I9U9X+y41Goqj617ty+GTmO47TECrCM7AvsMDprd6M3wzcjZ3ymkRS131aR72VcB4ZCAs7y2Ggi1Fo1iiL9yRLk6TZ7nKpnFE/hU24D43RQIc0VxyacEhbKTgmyUJbbADqZDGflOHu80LHnZWOtjGf7M8mtm5DeOsvkzLASHBhU9TBzfGibc/tm5DiO0xLb+ma0lPhm5DiO0wIrIQODiDypyXhV/Vrdsb4ZzQJLnAZoSeoWlS6Se4/VkO6qyor3m2qkYq5br6jCgw6MTJZIB5R/z6RqFEWzdi+kxpaPk/1VsUMRL72qOKKkNLeqLNPFpDkYSHKdhbLEBrCwUJbkrBy3YCW7Tlmys2PjMl25bZipeNEZtnUHBuBCBveJxByn8LLjjuM4U0VXhEx3mDneiZAQ9QpCtYabCRm8Xww8CnhNk4l9M3Icx2mBFeLAcFF+LCJnAF9W1VcMDfuYiJwGPI8GQbC+Ga1kYl5rkwa6ps6bNEN3qqx4bG1V3nQxDzpzXFlKPOFtF5XWEqXGYwGuyezYsbLjyazd5fPjxfUictzQcV+eM22dVYuD46y/uzBoK3jIFY7DmFUJD7pV3cVS+6qOuZaR2xYyea6TkOlSx9NgW9+MhjgSeFGi75MEa6k2Kyd3heM4zhKSOzA0fcwxHeCXE3370+B+EbhlNL9MmgJoXNqMLSpYJRWxRU1KnFti6YBS9YokYvnYFD/92kjx68YcGNJOCxGLKmFF9R0cKpwW7JiCtVNISVSOHbJOCwVnhVURBwXTn1tE1lHBWjirjMWUt682/auN5WMdFPL2hU7ZWoK4ZWRjjjosn2Wk8725NOWLwDtE5FbgM6q6KCJd4PnA24AvNJnMNyPHcZyWWAHedJY/IdSj+ySwVUTuAHYm7Cv/nvXXxjcjx3GcFtCV4U3XR1VvBf6PiBwO/AawB3ATcLGqXtB0Pt+MtlGmGltUa6zRiirWVukw0aSseER6Szkz9L9H7OVjqX9ScUZWRoumAyr32+O0A0N6/uH2XrecaTsaZ1SR4gcG8pyV5haMA8PCQu6UYKQ3I81ZSS4fs11366C/Y8cO2nMZbrtIGwzkuVWRtnC8PPWMViqqej6h0N5E+GbkOI7TErN6z2gpMye0hW9GjuM4rTDT3nEXUp0tAQZZFWp5wonIccCrgQPIqr5aVNUzMMw8S5wCqA3GKaQ39diiqqzc/bkGhwXJzlxDu+G4IK0V5qXcnzwu/k31h+uO7s/lu6SMV4gTioy1/5fHUvyksm5n7TFpDmD1wtbs76Btu649HshsuUy3prtlcH4nPnZ1JxxbGc561q2ScGzluLxtGJsyaBrMqmVEMXNCK4jISwkZGD4KPJpQsXsV8BzgFuBfmsznm5HjOE4LzHIGBps5oUVeB7wDeCvwCuDvVfVbIrIzwRK7rclkvhmtNGpmXZgKqVipSa3GmDNDldMCDBwXItYQmEwIsXNKx1VxRuXjVCLU/nrtWqocJ4y1Ez1OWUM2uWk/jqjsiAADi8haQ2sWjOVj2zOLyFpG21nLqLPVHIcxhQwNUraMbFs34cwwVTR41M0LIrIrwQvuwcDZqnq7iKwB7lfVOm/i/sDXgF72WA2gqneIyNuBtwMfrLsez8DgOI7TEj2k8WPaSODdwEZgPUFe2zfr/jzwpppT3Qt0VFWB/wUeZvruAX6xybp8M3Icx2kBJdwzavpYBk4AXgucDDyeok1/NvDsmvN8l0E6oK8DfykiTxCRXwdOAv6nyaJcpps3liMNUNvlxevWLko5Q8RS8FSl+0nGGUXaE84O/e+NComscJyQ1og4NqRkvF4szshKb51ye2osWSlvKcQTmRQ/NpFpJtlZaW67VUZOy50SjDS3xjgi2PYHZPLcaiPHPaAg2ZnzMpkuJs3ZY9tWSA1EWWGaTuzRTHvTWV4BnKyq78jS91iuBX6p5jzrGFhDfwVcQMi8AHA38Nwmi/LNyHEcpyXm5J7RnsAlib77ge3rTKKqnzTH14rIo4AnAGuBb2QZGmrjMp3jOE5LzIlM9xPgwETfo4EfVk0gIqtF5H2ZJAeAqv5cVS9Q1fVNNyJoyTISkSOADxACpT6sqqcM9R9PMA23EvzPX66qP876FgnaI8D1qvqcNtbkLBFtpgBKefE1kQVjMlxhrvAnlQ4oWla8IAOWr6UV59vz6qQD6v8kTMh4/fimVJxRTDYsSHNmrkySEyPH2VLhVqbrZ9q23nRmbB4bVCXN2eMHdO4fnG+kubXdQXsuv+VedbbNHq+Wwfkd4hm8+20R6a5tVGc6zsjyaeBEEfkWAwtJReThwJ8T5LeRqOr9IvKHwGfbWtTEllGmOZ4KPBN4JPBiEXnk0LD/Ag5W1V8DzgLeZfruVdXHZA/fiBzHmVvmpJ7RSQTngq8B12RtnyYYBdcAp8RPK/FfwK+2tag2ZLpDgGtV9TpVvZ9Q3e9IO0BVN6jqpuzpJcBeLVzXcRxnplBt/pj+GvVe4FDgGOAbBMeDS4HjgMOz7/E6/DnwehF5toxbDtrQhky3J3CDeb6R4C6Y4ljgXPN8jYhcRpDwTlHVz7WwphXLVLJ15xg5rf9vsc3rV6UAanpe3evZU2IyWyr4tBOR/CJyW/m8dBvYdEGjPegKc6XS/WQynBSK5FkPOpP6Jxtrg0+3W9haOq6S5mAgz9m2tUayW9OxnnVlb7o1Mhi7up8OqFdqG6YzBXnOMicyHaq6CHw8e4zLp4EHEWKTtorITynmv1NV3afuZG1sRrF3P7rfi8hLgIOBJ5vmvVX1RhF5GPBVEfmuqv4gcu5xhJ2bNaydfNWO4zgtoiybQ8Jy8RXqJV+tRRub0UZCtb+cvYAbhweJyNMIkb1PVtXNebuq3pj9vU5ELgQeC5Q2I1VdR3ZjbUfZZT4cKOeQFqztGhepKDHeaK5IjaHKRKn2/PKYlFNCzJEglQ6oHztUx4qKWUZRZ4jUXGYNucUjkTYGFlHKaaFQ3juvQWStIeusUBE7ZJ0V8nZrDa3t9r8GWCMmVimzjGxb0YEhS9Bq44wKDgxla6jb3nfmSObhi0lEVhMCX18M7E0527aqauXeoKrHtLmuNjajS4H9RWQ/gsvgUcDv2QEi8ljgQ8ARqvpT074zsElVN2d5kp5I0bnBcRzHaZd3A68h3C75DLB59PDpMPFmpKpbReS1wHkE1+7TVfVKETkZuExV1xNe/A7Ap7Nf3rkL9yOAD4lIj/C79BRVvWrSNTmO40ydJXLtrgqdyca8iOAlp8B/q+rvDY8xvAB4i6q+vfXFTkArcUaqeg5wzlDbieb4aYnzvkGLroEzzxzUMFoyxnEkqEoBVGf+aNbuqtRB8Xk1GmdE9HhQ+8j0V2XtHrO/mFookpXbyHCS9XetA0NBmjNxRln7KlsevFMuG/6AQo2iuGSXy3NWmis4MEjZmcHKdKsjMl0ha7eV6YxYNvUM3i3rdCZ05nDCLZFLRWS9/dEuIvsTZLcnZhmzH1Ix7Q7Axe2udHI8A4PjOE5LLEEGhsrQGeCVwKmqekdYw+BWSIKzgUZlyKeB56ZzHMdpiSWIG6oTOvNwABH5D4KUd5KqfmnEnH8HfCy7PXIOcPvwAFW9bpJFj4NvRiuBVKbvqpigtrN1D9GK516TEuORttgP0yrpLSnjxbzhahXiK58fld4K/fHieYPUQkamM8e5F531oLMync3Q3U8HZAvq2QzemSRXSOuTiB2KSW8xaQ5g+06Q8laxGO1fTSTOyMh0sQzd0/Cmy0tIjMGuWaxlzrrMexjqhc4sEArdHUrwZv66iByoqncmrpdLdCcBb0mMGc7mveT4ZuQ4jtMGSvzXTTW3qurBib46oTMbgUtUdQvwQxG5mrA5XZqY8+VMeHcruy+Vp4DbJWu+HbgKOLeGVFjCNyPHcZyWWAKZrjJ0BvgcIWbojCxE5uHAKJnts8B9Nt6zLiLSIZQT/zNCmfFNwB0EC24nQvmI+0Xk/cAJWRXYWvhm5CwdTVIDjZOpu6qtTn8hdVC+lkgbxqs2wmtEAAAgAElEQVStSsazx4WxozN8W2JphJLnFMbmkp6V6czQTjno1XrQ2YzXuRfdgsS96fIAWCvTFY+tN93mwl+IS3MwkPJSQa+5JLfKyHQd80PfvoZpBbv2aflyNUNnzgOeLiJXAYvAG1T1tth8IrIA3Ab8DsGRoSknAK8jbEgfV9UfDc2/D/AHhAQHd2fjauGbkeM4TissTTqgGqEzChyfParm2ioiNwPxZH7VvJJg8bw/Mf+PgbeJyM+BP8U3o22M5Sg1Xrj+mDWMUlRZTFWlwivXYGOKwp+kBWOIphGKzh8/jtdGSpwXs3wqkqom+6XswGCdFkSsA0MWZyTxOCPrwJCnBrJOC9YaWejXHdpq+m1y03L76kL/lpHHxfPNGjPLaHXBMhrQLX/802Me8gHBPxPqy51TNTDC7oTSEVV8KxtbG9+MHMdx2mB+iuv9CPg9EbmUkHH7Joa2UVU9PXHu9wj3rS6quMaLCTWTauObkeM4TlvMh2V0avZ3T+CgSL8Cqc3oZOAsETmAYGFdSXBgUIJX3aOA3ydUZnhBk0X5ZuQ0I1bDCNqtY9SESa8bcx6oExtUd2wqa3fMCaIyvsnWLUqUFc/aCzKdkeFiDgxWslsoyHC9clunHGdUkOYKpcK3ltoLsUfJrNyL2d+eGWv6s2/8VeY9skExXfPvcvoC91xYRvuNe6Kqfk5Efgt4B/BhytuvAP8N/Laqnjt8/ih8M3Icx2mLObCMMieDSc4/DzhPRPYiWEK7EDah24ErVfWGUeen8M3IcRynLeZgM8oRkV8j5Kh7MPAhVf1fEfll4GZVvbvqfFXdSAi4bQXfjJzpsBQF9Qy1PeFGXVfK/RrN2h05p7Se4t/ktQrXtW1a0V8xb8IRMfes69oUQRL3vMslOZu1u+DV1pfTyhJb+Xhr4S+M8qxbLPyFgTQX5sjPF9NvpLnlyo4/fgaGqSIi2xHu9zyP8C9FCTFH/0uoJ/d94I3TXpdn7XYcx1lZvB14GiE4dXeKP23OBZ6xHItyy2gbQCa8iT+VUuMTEq1hlCIWn5Q4v5FFFZu/MJcdM/R36DjmwJDK7KCRtrQVVXZg6ESsoE4izmghcmyzG1gHhnyOVYU4o/hxbgXZ5Ke232ZKyOOIYtZQmEuyfmMZGcu7s4y/sZcgHdBS8GLgzar6r1m9JMsPgX2nvyTfjBzHcdpjPjajBxPihWJ0gO2muJbChR3HcZw2UGn+mD4/BJ6Q6DsEuHqKa+njlpEzG8RilupIaC1KjJXfC1VxRtE548lRK2kg+cWcFQrdUpbpYm2l4+xnvpXmbHnv3EGhWOZ7cGwdFDrZeV3jAFGQ5qQc92TT+sScFVLSXDHOaLq/tyOllGaRjwF/KSI/Aj6TtamIHEbIxn1S0wmzkhJrhttV9fq6c/hm5DiO0wbKvMh07wIeDXycELgK8O+EzeRMVf27OpOIyI7AB4DfJS3t1S7S55uR4zhOKyyb7NYIVV0EjhKRUwmecw8hlJX4kqpW5ZyznAo8HzgN+C7QuD6SxTcjZ8VRLa2ZJ5NKa3XGVp0fTRcUiUOCSp1IIt50MWnOtncTY3PJrVOQ3mz5715prI096iQku/w4VbIpjyOqI811pp2eZ0YtIxHZRVVvt22q+nXg65GxR6nqmTWmfQahdtKplSNr4A4MjuM4baFjPKbD+ZmsNhIROYYg39VBaNHZwTcjx3GctpjdzWg/4EsiskNqgIgcR5Dczqs555nAb7ewNsBlOsdpzNi3BeqkBpqUviNiPMVPFVY6i/eX5ypIbInzY+W/C154heO8346V0nEdaa47zcKUs50O6BnA+cAXReQIVb3XdorInwDvJ9Q3elHNOb8MvF9EHkgo1Hf78ABV/WrdBbbySYnIESJytYhcKyKlnEYisp2IfDLr/6aI7Gv6TsjarxaRZUlD4TiO0waizR/TQFUvBZ4FPBZYn+WnC2sW+b+EjehTwAtUdUt8lhKfJ1hcx2TnXpA9zjd/azOxZZSlkzgVOJyQwfVSEVmvqleZYccCd6jqL4vIUcA7gd8VkUcSqgY+CvhF4AIReXjm7ZHk4dzBP+s5nMaBbJC9J30JjjPfjPGNZn+/dxJaUW4lJfsrNKa4NRR3WiiMkfxvfSujljX0mZ8h77idg+IF5dphRh0YAFT1GyLy2wQr5jMi8lzgL4G3EBKnHqOqo03jIoe1ub42ZLpDgGtV9ToAETkTOBKwm9GRDAKpzgI+KCEh2pEEv/bNwA9F5NpsvourLro7mziey0HxDclxnNF85mfI629B7p3h3WIKqOpF2Sa0nlCl9ZcI94mOU22WWa+hG3glbWxGewK2mNJG4PGpMaq6VUTuIuRH2hO4ZOjcPeteeA2LHMsVbMA3I8dx0sg7bp/KRjSrGRhE5GHm6Q+AvyBIc18ATgH2swmTc+Oi5ty7ENIL7UKIV7pk2I28Dm1sRjFbOlaKNjamzrlhguDpcRwUbezd2FRjiY7jrGh+srV6TBvMrgPDtcS/W58N/FakvVbmBBF5G/DnwGoG3+ebReRvVPWvmiywjc1oI/BQ83wv4MbEmI0isgA8iOB5UedcAFR1HbAO4GDjHnQLaydcvuM42zx7LsDGJd6QZjsd0MvanlBEXke453Qa4Z7T/wK/ALyEkPvuFlX927rztbEZXQrsLyL7AT8hOCT83tCY9cDRhHtBLwC+qqoqIuuBfxWR9xIcGPYH/rPuhe+jy2kc2MJLcBxnW0ZP2AVW8D0jVf3oEkz7KuADqvpnpu1q4CIRuQf4I2B6m1F2D+i1hECpLnC6ql4pIicDl6nqesLO+fHMQeF2woZFNu5TBGeHrcBrqjzpcm5mrXvTOQ6MJQ3Zr+ReItCpp53R/RUBUouR/kUTTRLrB1jUvN9m9R59rZ59RWrTEGXXe96OYcQ7bl9aC2ll7XX7Al9M9H0ReHWTyVoJelXVcwjugrbtRHN8H/DCxLlvJ5TBrc332ZmXyFPHWKnjOCuW5+2IPm9HLt/j2suX6hKz6sCwRNwGHEiIKRrmUVl/bTwDg+M0ZOwvnNR5bX6BZXOpWs+o+pZTbg2l+2PWzqBtMXF+3EoatPUKx3m/HVvO1tCxnsiFzA9mDY3CZlpgZW1GnwXeKiK3EUJ0tmQ+AS8ETgYaSYOem85xHKctZjc33VJwAvBtwqazSURuBu4F/gX4b4JzQ23cMnIcx2mBaab3mQVU9W4ReRLBNfz/EOKMbgcuAs5tGkTrm5Gz4rBfGLH/Wwr9Tf530sRx1diq87XcJgUZzi54tCSXS3ZWbiscU25fTIzNJTkr7VkHhUJ7drxFB+Erq40Qt2gcFHL5zjolWLGtl73enpXmrDyYkuymwezGGS0J2YbzhewxES7TOY7jtMWcyHQi8lgR+YyI3CoiW0XkcVn7X4vIEcuxJt+MHMdxWmJWs3YX1ijy/xFiPn8F+FeK+0CPED80dVymc2YDKzX1suM6/6c2k6VHkl8uOWOko2qJYtYnTSScCpkudayR16ARmS3WVjrO5LKtvYG0ttgdfG/lkttiRNoDuN9Icmuy7zvrbWeluS2mfVUmylkPuvsLGmX+wgZy3Cr71hZkw1phi+0xH/eMTiHEhT6XEBv6WtP3LeCly7Eo34wcx3HaYH4cGB4HPC/LgjO84luB3ZZhTb4ZOY7jtMZ8bEb3QTKp5x7AXVNcSx/fjLYBtGekoM4YqWGslNSgoNk0EfMaK+Iyi9JdfmzOt/mIrYymdet/F6TBwTnSQFqTyHGszV4h6eVn5b9chusZ6axXltEKXnG9wRu6NXLc64ppG7x5vezf2pbe4Gtki5hjk9krl+xWSdf0D8auksHYLZmk143pj+FJ+FOQ5gaSXaw0+tSYj83o34HXicjnTVu+8mOB2qXC28Q3I8dxnJaYE5nur4D/IASmnkXYiI7OElYfBPx66kQR+SENtlxVfVj1qIBvRs50KKRlyX4dj+t8EDmvkYWTsGwG1ox1OqB0XHBESLyEmLVTXEPkOGntxLwSKuZNGBW548JixFqy/TCwgrYYa2hLp+zAYGOH0scL2V9jLRnngo6uSrygDGM5DV7Q4IVZZ4eu+UzdXbiMqv53Fqz6buBNhP8JXgt8HXiyql494vSLKP7reyqwO2Fzuzk7fiKhnMRXmqzLNyPHcZy2mA/LCFX9FvBUEVlDyJxwp6pWVipV1WPy46zg6eOB31TVjab9oQRvvYubrMl/ODiO47TBGDFGyy3rqep9qnojsEZEDhKR7Rqc/gbgLXYjyua8ATiJUNq8Nm4ZOc3oDeQ2NZJNPwNLrWLFba4n+7953J9VEfUw7ShQPj3qtGCPrcxnrhV1Vkg4MAzmMs4SNhOOdc7olR0Y1Dgl9LLPr9ezqXqMNGe8Q/LjQpuR7DZnjgvbdQb1ge7rDeQ265SQt3fNG97pmLpDtgZRNqabMDMWs/6encv02/P813YZEXkzsL2qnpA9fxIhnc/2wE9E5Kmqek2NqfYieObF2Azs2WRd/lk5juO0RZ30P8OP6fMS4Drz/F0EZ4bnEu77vLXmPFcBb8ikvj4i8gCC1XRVk0W5ZeQ4jtMW83HPaE/gGgAR2Y3gPfdUVb1QRFZTv1T4/yVUdL1eRM5h4MDwLOBBwDObLMo3o3nAeqLJMhizhSJmVZ5qEd1rmFxW6ibm6nuPpfSyCgqphcIf+7alpso98qq98eLHsdQ/SRmvV24vSG8ReS/Zb73hsve2INNpOeZoMRFntGVx8Jmt6gSZ7X7TtrpjZbzQvrkXjxe6z8YRaZDybIog60HXrSiCZ9P69PpZva3XnPWmix8vNcLy3wOqySKwOjt+EkFq+4/s+S0Eh4ZKVPUrIvJY4M2EEhJ7ADcBXwbepqr/02RRLtM5juO0xRLIdCJyhIhcLSLXisgbR4x7gYioiBxcMeWVwEtEZAfg5cBFqrol63so8NPqVQVU9Xuq+vuq+kuqujb7+5KmGxG4ZeQsJYWsBxXWhnGMoFvhBREzbVKxQ6POGW7PrYVepA1joVQ5GtjjCsvJzmuRmOWUcIYoGBX91yDR/twysg4MWxeNM4P5nPL4oq3GmrnfWEGbF8PEC8YaspZRV6xTQgOTIVuOTaq6prOlf9zLrKT7jXm5uuDMUL7WVCykJfCOE5EucCpwOLARuFRE1qvqVUPjHgj8CfDNGtOeDHwe+H1gC/AM0/csQrLUJmvsAI8EHgxcpqo/b3J+jltGjuM4bdG+ZXQIcK2qXqeq9wNnAkdGxr2V4IiQ8m4bLFH1POARwIuAR6nqRab7a8A7K1eVISKvIQS4foeQRuiArP1zIvIndecB34wcx3Hao/3NaE/gBvN8I0Mu09l9m4eqau1qq6r6Q1X9N1X9wVD7h1T1kjpziMgrgQ8AnyNsbFaS+Drw/LrrAZfpVgYpB4h+jE5CQsuls87S/GZpJUFrTH4rvNyyU4LEpDkG3wspGS6ShaZwfenZa8TGjj5OSW/Sjx0qxxPBUKacrpb6i0lTs7pCi4O5FjtxB4Zu9iLu7wzaFsQeh6+PhZ5J62PesMJx9iI75kV2Io4mljVipDnrZCHBGWIVNrmqjVmyaxjtGNE2Y8p0u4rIZeb5OlVdl08ZGd+/SiaRvQ84pulFRWRnYH9gzXCfqn6txhTHA+9R1b/I5ETL/xDcu2vjm5HjOE5bjLcZ3aqqKaeDjQSngpy9gBvN8wcCBwIXZj/ofgFYLyLPUVW7wfXJ4oJOp2zNWOqEr+9HSPsT4+fATjXm6OMyneM4ThuMI9FVb16XAvuLyH5ZDNBRwPr+JVXvUtVdVXVfVd0XuARIbkQZfwUcChzNIEnqKwilJX4APLvmK74V2DfRdwDwk5rzABNaRiKyC/DJbEE/Al6kqncMjXkM8A/AjgT/9rer6iezvjOAJzMo5nSMqn57kjXNNE3idbY1Cq+95lgrSxn9I1rbqMpbLlbWHKBjPdwySS+WMdssoSjjja5nVEgBVHE8bn9U0rPSnfGWy2W8xUUTZ9SNe9blGbpt1m7rTbeQybidRRMvZN6Ee2PaW4Ke+f+hl/1G7omV5qzHXljPatPWsTJd4c3J2qYUb9S2N52qbhWR1xIskC5wuqpeKSInEzzX1o+eIcrzCR51ZwIfB76ZJU79iIh8GjgCOLfGPGcDJ4rIhcCP8yWLyK7AnxHuJdVmUsvojcBXVHV/QrrwmA/8JuClqvoowot8v4hY8+0NqvqY7LHtbkSO42z7LEGckaqeo6oPz2J43p61nRjbiFT10AqrCGBv4EpVXSS4dm9v+k4Hfrd6VUAIdt0MXAFcQHg1fwt8j2B4nFxzHmDyzehI4KPZ8UcJuY0KqOr386R7WXbYn7JMNdYdx3GWkjnJ2n0bsEN2fAPwaNO3K/CAOpOo6m3AwcA7gFUEiW8B+CDwBFVtVL58UgeG3VX1pmxhN4nIQ0YNFpFDCGkorDvh20XkRDLLSlU3T7gmZwKmUoK8ryuZe6STFtorlF43Ml6uCSYcCguX7acOstKclMZa77XCbV4r32XrSQW3Wqkxl9aSWb1ryngAmstvXTOBkeTy/p7xoLOedVtNezfzrNuc8pATLbXduzje51hYT/Yb+X5TwnyN3N8/Xp3JdPeZF7664FI4OOw0kApbYXk2l6ZcAjyWIMX9G/DWLGh2K/DnhHtHI8m85w4EblTVt1I/uWqSys1IRC4geGgM86YmFxKRPQj65NGq/W+jEwgBU6uBdYT6F1HTLivkdBzAGtY2ubTjOM7SU1N2mwHeSZDqAN4G/DLhe7dL2KheXWMOBS4DfouQi25iKjcjVX1aqk9EbhaRPTKraA8SOY1EZEdCdtc324Cq3KoCNovIR4DXj1jHOsKGxY6yy3x85MtAHosiqdihNonUNiqEdUxa2yiVKLXKipo0wap1YIhYMCZLTdra0bK1k6xtVGX5ZD/6pWutLHMtY/lIbhFZa8j6Y2QOCmqsv0VjKtpwkdwZoWusls1bB18Z/dihMb+Bbblzm/Q0L1e+nUkBtMWsK085tFpMHaWCA0N5Pd0pWEhCtW/OLJDdU7osO74beH5WVG87Vf1ZzTl6InIDxftNEzHpPaP1BPdAsr+fHx6QuSN+FviYqn56qG+P7K8Q7jddMeF6HMdxlo8lcGCYBqq6ue5GZPgQ8LrsO35iJr1ndArwKRE5FrgeeCFAljX2Var6CkJg1ZOAB4vIMdl5uQv3v2T1NAT4NvCqCdfjOI7jVJCpVc8iyHXDGRg0uw9UxQOBXwKuE5EvEcpHFAIcVPUtddc00WaUeVM8NdJ+GSGIClX9Z+CfE+c/ZZLrO8vAuLWNTDqXfgbvVFbvXP6ryt5t15OS4yJOCQUnjVhqoESKn/y8VLqgmLNCzFHBrqtwHGuDQcySuUdvnTAKx9kYtTKtee1szaVGU4rcnG+dGbZkaYDExnhZB4bFhVJbFYU6SpSdFmCQoXuLlqU5e2zbOoUM4WVJrjMlt7V5qGckIk8kxAilMiQo9RwS/tIcvzwxT+3NyDMwOI7jtMV8yHTvJyQp+HVgjap2hh617vZGzhtrnhzPTec4jtMWc2AZkZWPUNXLl3shFt+M5g3t607Tu6YtfNdGBu9+wI5tisQ3pVL4RFIDaTS9tk3xk/Cwi6UGshmxIzFHGvF0C2PNcT4m4UGnES+8QtZv+5YvlpZV8JYUE1OUe9bZddus7LknoG41nmxSmHhwmL2n5YTM9bAecrnnnC3e1zNuiVtNZvA81mmhM3gTqmQ6mwJo2pm6+yxfEGtTrge2a2OizPnst8n8AoCTVPXHIvJk4Jos0UEtXKZzHMdpi/mQ6f5/4I2ZE8PYZCUovkHIQfcK4KWEDQnglcTTwyVxy8hxHKclZtUyEpGPDTXtDvxQRC4Gbh/qU1U9mmreTShv8URCdvH7Td8FeD0jZySxQnu9uKQz1bW0kRooJybpJQq5RVMDJYrvyWI5oLiw1KJTaxibkN6iQa8JyW+QmTxxfuQ8Gwhb0D9ymc5Ic2ouZqXAfNqU02SVF53Nfp7LdDbQdWtnsPCtJnP4luwFrzLy8EKnLNNZD7lVhRxNA7rTluxmdDMiyGjD/1p/BjwqMrbuqzgSeL2qXhwprnc9xTpMlfhm5DiO0xKzahlltY7aZgfSNYvW0DAhhW9Gy8Uc1Dbqx9XUWV/2K9bGrzRKDRSLObKOE1LwdigfV6UAipQiB1Bzjf6P60QmVe2n+LEOFCSOI8lPEzFH/XQ/nXh/3yiwFeOtz4F5b6SfvDQ+tm8RFZwWbL+1ksJcWys+PC0kky1bQwBbM2snaRkZZ4bNWYLU7bqDdD8LvbKDwqqE08K0YopKzFBGhVFk9YbuUdX7JpzqauDpBElumCcD320ymTswOI7jtMWMOjCISFdEThKRO4GbgZ+JyL8N1ZZryqmEdEBvYpB4dScReRmheuypTSZzy8hxHKcFhNmV6Qip1k4ELiQ4GzwM+B3CfaOXjTOhqv6TiPwSwTsvr7ZwPkGHeJeq/kuT+Xwz2kbRQn2fpapLNGFqoG78nEp5MOKgUChFHok5koLcZp0ZbExR4ZSASY+TS2CxMt+QcCRISXMRea9jS/JYzSJX1ipSAAF0Onm9InNdo8PlElav4MBAlEHZ78FirWQ3UEirZbrFiAPDopF07zelzVdnb4RtszFHC9kbVkgBZNMUGXNj6pLd7G5GrwT+SVX/MG8QkT8EPigif6iq96dPTaOqbxSRfwAOBx5CKNx3vqpe13Qu34wcx3FaQib1BF06Hka5RM8ngX8A9gGuaTqhiDwJ+Jaq/hj48FDfDsDjVPVrdefze0aO4zhtMM79ountXTsQJDnL3dnfB4455wbgkYm+A7L+2rhlNK/E4oWmQZupgWI1v4G+613KQy4m3yWzdmftFemCQnNZ0itWpivHDhXSGEUL5tUvO24/0o6RB3uZ9FaQASV+nAcHxTzsgH7W7pQCWnwXyx5wg+gjiLlIFjMtGUkue42LJkX4FpMCaFV3sdS+qmOzcg8mjsl0sXLow8fTYIbvGQHsKSIPM8+7pv1OO7CmzDZKm9+O4j+WSnwzchzHaYvZ3ozOSrR/LtIW9ecXkX0Jkl/OwZkkZ3kAoaTE9U0W55vRSiaWdDWSjSGaxLTO9AmroXbMUazGUViEvUi6zR4nMjDELCYxVknBGSIWD2R/+1lrpD920G2tnaJpU34NxaSoZYeOgk+K/fgi0xdijvptUmoLJw4OY9+rvYIDQ680rlewGiMZGMznuNg1jhHGYsqdFazl1O2UraBujdiiZYs5mj3G8piLcDShRlEuMv4dRQtJs+dbgdc0mdg3I8dxnJaY1b1PVT/a0lRnENzDBfgqYcO5amjMZuD7qjqc824kvhk5juO0xYxuRm2Rec79GEBEDiN40909+qx6+GY0C8Skpjann2bMUR0Hh1wetPWuI3E1SXkwVo8oUuMIjMxm3wMj7xTSF2UympXmos4MvcS6ImmCOltNt5Xxos4O5bawrvxamEZzGDkuXNfGGZHHGWHaEpJdROG0X7aap/ix0tyCOTafby7JFaU5mwh1cLwlO+4aOTQu01VLc01Kok+Mzq5ltBSo6kVtzueu3Y7jOG0xu67drSMiq0XkLSLyPyKySUQWhx5bq2cZ4JaR4zhOC8x4OqCl4N2Ee0bnAp8h3CsaG9+MnPFpuxx5Tj/gpkaNo6qs3VX9Mc86Kw/GPOsK8T5lbzswcmhBPjT9hTny2J+4/NeX3kx/L9IfxmT9mDZznEt2KfkqKtlFYq3CYSZrWlnUvs1d257HUpnYo8WyjAfQyeS5jpHmuhHpTQrxRJT6C69rWrvE7GZgWApeALxFVd/exmS+GTmO47TECrOMdgAubmsyv2fkOI7TBuPcL5rvzetsQgXZVnDLaFtg0tRAqfNz+aVQZrtBAGwv4bXWj7ysiIS18o/NUm3lwW63NDbmWVcIXq3wrEsFvfY966xsZcd2rLdcfl3TZm7nFjzgsrEdWyrcrjcb3EtIUTFprSjNYcb2Zxis23ymhZjnXFqLFBgMT/I2e7HBFaxMt5gd90wgc8f094xkJ9m/i475N9OxgcZ9mY5S2/DxtJl2lfNl5u+Aj0n4n+gcoBRX1CR790SbkYjsQsj8ui/wI+BFqnpHZNwig6p/16vqc7L2/YAzgV2AbwF/MG4qc8dxnGVnvi2dpuQS3UmErAwxqvKt9JnUMnoj8BVVPUVE3pg9/4vIuHtV9TGR9ncC71PVM0XkH4FjCSnNnXll3HLqhZ/k/UYzIOHMELth3MSBIebMUHBKKP/iLlhZEQsnnJefQ3SsjQPqWz7Ya0UcGIhjz+tErajyOTGnhtBedkYoJEq1H4lG2rpxiyv/fHWxbC3BUIqlfk0m+xoilk+hLT522qywe0Yvp8Xtd9LN6Ejg0Oz4o4Q0EbHNqISE/9ueAvyeOf8kfDNyHGceUVaUN52qntHmfJM6MOyuqjcBZH8fkhi3RkQuE5FLROS5WduDgTtVNf+NuBHYc8L1OI7jLBuizR/zjoh0RORAEXmyiGw/7jyVlpGIXAD8QqTrTQ2us7eq3pjV0viqiHyXcqEnGGHyichxwHEAa1jb4NJzxLgSV5NLjJMaKJLJuzwmz39T/fumXzdonEzeEM/m3e3Gx1aVJY84MxScLezN++xaxUzd5vUWsnITGRt3UJBIDaJO4bzs5j9l6Q6Gs3JnzhCJUjN96c20dQr9NtN21l9wYCiPtf4v9rh4kWyskeaKzh9mjk7/wv02K+P139uEXGpZov+N0mwDm0sTROQ1hPtFuxJe/a8D3xKRzwFfVdW/rTtX5Wakqk8bsZCbRWQPVb1JRPYAfpqY48bs73UiciHwWODfgJ1EZCGzjvYCbhyxjnXAOoAdZZcV9pE7jjPrrLQMDCLySuADwOnAl4FPme6vA88Ham9Gk8p06wn1Lcj+fn54gIjsLCLbZce7Ak8ErtLw83gDIYo3eb7jOM5coII5SjwAABfeSURBVDreY345HniPqh4HfHao738IpcdrM6kDwynAp0TkWEJVvxcCiMjBwKtU9RXAI4APZb7oHeAUVc3rX/wFcKaIvA34L+C0CdfjxArmjXN+jTkqY44K3lQVqYPsdXM3qkgm79J1Y9eKraGON12suJ7VmiRIgZKQ+Qo/i3PVEttkYprMe5vLdIV0P+a8/L21nm6FQnxmbKfflohD0uI4GPqYCnFGWb/Nym0G9A9taXUjL2o3cmzlRxu/1o28NvN+aOHfQu7VaNuIs5JMlemzH3Beou/nwE5NJptoM1LV24CnRtovA16RHX8D+NXE+dcBh0yyBsdxnFlhhe19txJiTGMcAPykyWSegWEFk//KbbXG0ZInT6UY5Z8dSuq6desdJcYWfvXHHByqSpSbl1B0UIhZQfF4n7w5FQ9ELANDwUgzVlJe/tucTcQaAhtnZMdKaayak6w1FLOYCtaS/edhs0/k85n3Szsxy4n6TC1R6nQuMyOcDZyY+QH8OGvT7HbMnwGfazKZ56ZzHMdpiRXm2v1mQtmIK4ALCFvx3wLfI/wcO7nJZL4ZOY7jtIESzMimjzklu01zMPAOYBXwA4La9kHgCap6V5P5XKabVaYQc1S9hogzRCLmqB87lFprxJkhmjwVBglUU2XJY2mCxnRgkEXjVBA53b6aaCJVK7LZm/PZggtiWqr2UX4Re7FFKfVbaS4p2UVkuKKzghbGhTaix4M4I9MfiSmyCWIL/RH5TqwzhHVKsYvslPutlNxvT5kV0X+CU/p/aH73lrFQ1buBt2aPifDNKMJhej3HcgW7sYlbWMtpHMgG2Xu5l+U42yS//a3Lef2XzuUX77yDG3famb/5rSNYf9BBy72ssZhz2a2SLEznOOAyVY3WMhKR3wQOAv5RVbfUndtluiEO0+s5nsvZnU10gN3ZxPFczmF6/XIvzXG2OX77W5fz1//2afa68w46wF533sFff+osnnP55cu9tPFYgjgjETlCRK4WkWuzhNTD/ceLyFUi8h0R+YqI7LMkry1wHCGH6PdHjLmakJXh5U0mdstoiGO5gjWFYAhYwyLHcgUbmCPraNIaR22wVJ51hWtkMSepekf5dQvafGJdlZJedq1I7BEMZ4uOSJyL8S+eQWxQnEF/XJqLt8dluMoUPzbddy7DWektkjy9ILcVZDpb9ylvK8YOveFL57J2S/HH89otW3jDF7/E2b928KDRfkwxWdOyTKo2tG8ZiUgXOBU4nJC/81IRWW9iNSHEaB6sqptE5NXAu4DfbXclfV4IrMvuF0VR1dtEZB0hCfaH6k7sltEQu7GpUbvjOOOzx52l8mcj22caHfMxmkOAa1X1uqzW25mEagmDy6puUNX8C+oSQmq1peLXgItqjPtaNrY2vhkNcUsiCWuq3XGc8blpp50btc8yITedNn4Au2ZVDfLHcWbaPYEbzPOq6gbHAue2/uIGrCFkV6ji58ADmkzsMt0Qp3Egx3N5Qaq7jy6nceDyLWqJPesqM3lXlSWHQVbmOmXJ+9GUDcqSV3rWJYrv5deQsvwU2svBsKms3v2xkUBYAI0V9U7Jg/ZF5N5yW+MjI2eMKIgn+YFpMyf2S4kT7e+ZJ52KrNwDb7r4XDFvuYJXXAfe+5Rn8bazP8UDjFS3adUq3vO0Zw4VIIzMm/r3tYwyXVJrHc2tqnpwoi/2aqL2lIi8hOBq/eSxVlGP/wX2JyRCHcX+2dja+GY0xAbZGxT3pnOcKfCFXwtec8d/5Rz2uOsObnrQzrzn8GfyhUfPqzdd6+50G4GHmufR6gYi8jRCWZ8nq+rmthdh2AC8SkTOUNXo1isiHeBVwFebTOybUYQNsvd8OStUMWfODNF6RykrKa9tZFP8lI2O6jgkc9wo9qhgDZUtPUk4LRTudOcWgPm/sdJKSlk2mYHYiaQICudl1l+kblFp3txStEZnIWlq9tdaRgUHhvJx7LM551cO4ouPNJtPBzpbSCdCzY2/QtsM+FTXuwfUlEuB/UVkP0Kut6MYVMcGQEQeS3AUOEJVo2V8WuRdwLcICbL/aPh6IvIQ4O8J+UiPbjKxb0aO4zitUM9Vu9GMqltF5LWE7Nhd4HRVvVJETibE+qwH3g3sAHw6k8avV9XntLqQwXq+JyIvBT4K3CAilzLIS7cPQSbsAX+gqt9rMrdvRo7jOC2xFAaaqp4DnDPUdqI5ThZAXQpU9dMi8l+EekZPBR6Xdd1AKAP0PlW9tum8vhk5BSZ2ZoikCIIW6h2l6M8xYVZvO0cs9sgeGxkvKdnlNXcKObmLLghmQeGPvWEfkexSmbZj2b6LTgvmsFt2NCmkBorEFBXbtDQ2nQ4ocpxIAVSsV1RaVpGoTDcjzgzzXSyvNtlm80dtzumu3Y7jOM6y45aR4zhOG+iQ043TCN+M5o1cBhg33miWPOtSclzdrN4QL1FeldXbSnax+CMbexTzrEvIeNatV7NrWGmuWrIzbVay65Z6Cx5usbCpnpX5tCydFlL12HUXXpqU22IyXqcs3Q1fo0qmi7U38qaj3L8srBCZbinwzchxHKctfC8aG9+MnCSVZcljVlZFvSOo4czQxOrL1xCLPTLricYeDV83x1hDdK11FkmUmiC3fFLWULw9VYUogn2bbXNmLhSyMkQMuZ6xrAqxRZHPrOhoUBFnVDiOrKFgAdmCRaY9t3wSltOoc5abJQh6XTH4ZuQ4jtMWvhmNjW9GjuM4baCMm5vOwTej+aWN5KmxsuLTpCq2qCqRKhSkwMHgSOxQJPYozFWxhojjQ0FuK5QSL38T1XFgyNuTcUgauWO/kHB2yKUxO9ZIcnl7oZ5RwdGgnHS1SnrTWEnw0nnZQaHf/hsuLXGoLSLppdTjZZLsBHWZbgJ8M3Icx2kL34zGxjcjx3GctvDNaGwm2oxEZBfgk8C+wI+AF6nqHUNjDgPeZ5p+BThKVT8nImcQam/clfUdo6rfnmRNTvtUpggKg/IBg7YKz7qCV11MLkvIjwXPvNq1jxLpggrS2dD1Ie5Z14ucA2g3cV5/7GjJLh2HFGGrGdu17+3oekV9aS3SFl5D+bykN13f6y3huRfxliv0J+KIomXFjRxaVXZ82Tzr/J7RREx6s+CNwFdUdX/gK9nzAllJ3Meo6mOApwCbgC+bIW/I+30jchxnnhmz0qvD5JvRkYRU4mR/n1sx/gXAuaZeu+M4zraDavOHA0x+z2h3Vb0JQFVvygorjeIo4L1DbW8XkRPJLKslrlK4bTKpZ12DFEG1JLsY45Qob5LVuzD/GOmC7PWqMnwXo0z7h4XUQblkF5HrIOFNZ9P2WI+/vN2+31YStN9neXPBg872R7ztCt6F5rRc0uuW2+wUaWku4nmXkuZiHnlVMlyNf37Tlex8c5mEys1IRC4AfiHS9aYmFxKRPQjV/84zzScQ6qSvBtYBfwGcnDj/OOA4gDWsbXJpx3GcpUfxzWgCKjejUYWbRORmEdkjs4r2AEaVvH0R8FlV3WLmvik73CwiHwFeP2Id6wgbFjvKLv6JzyptJGKNJVJNWEn98t5VFmHBcoq392sfxdZSIBEPZGOhMouoyqkhXC9zYCjc77cJWjvldcXikOyxaat0cCiUFTeX6NeCMv3RUuLxOKV4ip9UCqDyeclEqf1zZqSGkcUdGMZm0ntG6xnUOT8a+PyIsS8GPmEbsg0MCd8kzwWumHA9juM4y4Y7MIzPpJvRKcDhInINcHj2HBE5WEQ+nA8SkX2BhwIXDZ3/LyLyXeC7wK7A2yZcj+M4jjOHTOTAoKq3EWqgD7dfBrzCPP8RsGdk3FMmub4TYYr1jpa8RHmdUuR1ax+lXouWnR2SGb7HkOwKTg0VDh9WaiyUMM/nLaT4SRUG6k8QP444OEhBqjRriJQdLzg7VDklWPkvf22xeCHiL8H+G45fIz7XsuKWzth4BgbHcZw2UIr3Jp1G+GbkOI7TCu7aPQm+GTnLR5N0QZaKuKrx0gVBZYbv/rCUy1REsjPrS6UDGgywYw25B1w3MZdNTxSLSarwtiu8hxH5rhgDZKSxbD2SkNsk5vCXkuZiGboTkl7snBRTl+98Mxob34wcx3HawjejsfHNyHEcpw38ntFE+Ga0rTLl4nu5Z11lVu+a85VIedZVpPDpS34xDztoluG7dlCsmayg3I0urpf8nLL3NCnNdSMyXMrbLv+ytMG4NuA0Jt91y01hDfn0ca/K2FQpdbAya7ch7nkXHztddMj10GmCb0aO4zht4TLd2Phm5FTTdiLVmrWPkolUY1ZQRVLVqFMDNEuqOk4cUiE2qBMdK30HhXh/7C58KiYpT4SaTLqaW3cJYyhmMUnPzmWcKLJraU9K54R5y+2FazWwoixRA3IWLCOX6SbCNyPHcZy2cMtobHwzchzHaQvfjMbGN6OVQBvODNOkbrogqJ/hO5lOKFIHadI4pNR1bdZuKct3kuiPphGyMUmRGuIF2UsXy/1WJqzyZbDXt5Jc3i7xf1/SKUt2MQeI4fY8zU8yRijyfsxGOiAPep0E34wcx3HaQBnhXelU4ZuR4zhOW7hlNDa+Ga00ZqlEedVcFemCoEGG70hBvsL59noTxiFBInVQSirM5blOxKUM421nzynoXeXXIIWy5Oa9y9ZTeLvNwqPynZHmCl58eSxTxGtu+LxCe2QuCsvN46rsYHO8qIVxpbHLiW9GY+ObkeM4Tiuou3ZPgG9GzvjMkJVUGYeU6I9aSTGnBqiOQ7J+AvlrsHE5Kau0v4aU40R2Sh0Hh5i124lYNhFrCYbe+syyse8RNtuC5k4JibpDdu157FgNK0qqHBj6llONL/5Jy4c2QUE9A8PYTPOjchzHcZwobhk5juO0hct0Y+Ob0Upm0hLlhbmWIKlqg4SqjeKQKiS7qFMD1IhDsowRk2SJ5cepcHCw7RqJPUrNJQlniPywEKfUi1yrEFtUfin2vDqS3iAxrF1W4RmlASkWq4e0ijswjI1vRk4lh/V+zLFcwW5s4hbWchoHsqGzz3Iva245dNP3edk9l7Bb7x5u6ezAGTs+gQvXHrDcy3ImRdXjjCbANyNnJIf1fszxXM6a7Cfm7mzieC6HHr4hjcGhm77P6352IWvYCsDuvXv407s2AHDhDo9YzqU5beCW0dj4ZuSMjD06liv6G1HOGhY5livYwGSb0dgedgnPuv5pTeKQIpKdJuKBKkuYWxIxSS+755L+RpSzRrdyzM8uLlpHsfRGqWJAEflOIrFHYDzvInLb8LFUjO3HA0U88Erz9j3k6kt6lmQJ+v5aRnZPDXXLaGx8M3JGshubGrU7o9mtd0+jdmee8Nx0k+Cu3c5IbmFto3ZnNLd0dmjU7swReT2jpg8HcMvIqeA0DizcMwK4jy6ncWD8hEnLi6eokuwSnnm1g2LDoJH9lSXMLYkA2Y9sfwivu/trBanuPhb4yPaPLwafxqTTZOZxu574kP5Ui2UPymTQbN+dLi4JDiRQe4GyjGfXVfQiTEXIxrJyj9bh2vynNhEe9Do2E32EIvJCEblSRHoicvCIcUeIyNUicq2IvNG07yci3xSRa0TkkyKyepL1OO2zobMP7+UgbmYtPeBm1vJeDnLnhTG5cM3+vP+BT+Lmzg7h/ezswPsf+GQufMDDl3tpzoQo4T5o04cTmNQyugJ4HvCh1AAR6QKnAocDG4FLRWS9ql4FvBN4n6qeKSL/CBwL/MOEa3ImIfKLfENnn/GcFdouVz5q3nGcGqBRgtXqEuaRYj2RmKQLt/tlLtzul81rkDC+KiZpcbF4TmzdFdbMoKZ3Oa3P8Hl9iymVlDW6lkTy0l5kLpPOJxpHlIrLipGysqaJqltGEzCRZaSq31PVqyuGHQJcq6rXqer9wJnAkRK+EZ4CnJWN+yjw3EnW4ziOs5y4ZTQ+07hntCdwg3m+EXg88GDgTlXdatr3nMJ6HMdxlga3jMamcjMSkQuAX4h0vUlVP1/jGjH7WUe0p9ZxHHBc9nTzBXrWFTWuvdTsCty63ItgqdbR/EdbO+uYPIXLtv25NMfXUWRJ0l3czR3nXaBn7TrGqbPwniw7lZuRqj5twmtsBB5qnu8F3Ej4AHYSkYXMOsrbU+tYB6wDEJHLVDXpMDEtfB2+Dl/HfK5jKeZV1SOWYt6VwjQcIi8F9s8851YDRwHrNdwF3gC8IBt3NFDH0nIcx3G2MSZ17f4dEdkIPAH4ooicl7X/ooicA5BZPa8FzgO+B3xKVa/MpvgL4HgRuZZwD+m0SdbjOI7jzCcTOTCo6meBz0babwSeZZ6fA5wTGXcdwduuKevGOGcp8HUU8XUU8XUU8XU4SUQ9l5LjOI6zzMxKEg3HcRxnBTOzm9GspBoSkV1E5PxsnvNFZOfImMNE5NvmcZ+IPDfrO0NEfmj6HrNU68jGLZprrTft03w/HiMiF2ef33dE5HdN39jvR+qzNv3bZa/t2uy17mv6TsjarxaRZzR/5Y3WcbyIXJW99q+IyD6mL/r5LNE6jhGRW8z1XmH6js4+w2tE5OglXsf7zBq+LyJ3mr4234/TReSnIhIN+5DA32br/I6IPM70tfZ+OGOiqjP5AB5BiAe4EDg4MaYL/AB4GLAa+G/gkVnfp4CjsuN/BF495jreBbwxO34j8M6K8bsAtwNrs+dnAC9o4f2otQ7gnkT71N4P4OHA/tnxLwI3ATtN8n6M+qzNmD8C/jE7Pgr4ZHb8yGz8dsB+2TzdMV9/nXUcZj7/V+frGPX5LNE6jgE+mPg3el32d+fseOelWsfQ+D8GTm/7/cjmehLwOOCKRP+zgHMJMY6/AXyz7ffDH+M/ZtYy0tlJNXRkdn7deV4AnKuqbRf8abqOPtN+P1T1+6p6TXZ8I/BTYLcxr5cT/axHrO0s4KnZaz8SOFNVN6vqD4FrGc9xptY6VHWD+fwvIcTQtU2d9yPFM4DzVfV2Vb0DOB8YN0am6TpeDHxizGuNRFW/RvghmOJI4GMauIQQ57gH7b4fzpjM7GZUk1iqoT1pN9XQ7qp6E0D29yEV44+i/D/b2zNZ4H0ist0Sr2ONiFwmIpfkUiHL+H6IyCGEX8w/MM3jvB+pzzo6JnutdxFee51z69J0rmMJv8ZzYp/PUq7j+dl7fZaI5MHny/J+ZHLlfsBXTXNb70cdUmtt8/1wxmRZ6xnJ7KQaSq6jxhrsPHsAv0qIqco5AfhfwhfyOkJs1clLuI69VfVGEXkY8FUR+S7ws8i4ab0fHweOVu0n7ar9fgxPF2kbfg2t/HtoYR1hoMhLgIOBJ5vm0uejqj+Ind/COs4GPqGqm0XkVQSr8Sk1z21zHTlHAWepqk341Nb7UYdp/PtwxmRZNyOdnVRDyXWIyM0isoeq3pR9uf50xHpeBHxWVbeYuW/KDjeLyEeA1y/lOjJZDFW9TkQuBB4L/BtTfj9EZEfgi8CbM0kkn7v2+zFE6rOOjdkoIgvAgwiyTZ1z61JrLhF5GmHzfrKqbs7bE5/POF++letQ1dvM038ilGzJzz106NwLx1hDrXUYjgJeM7TGtt6POqTW2ub74YzJvMt000g1tD47v848JT08+8L+f+3dP0vDQBjH8e+56ODin0UnFURHhQ6ie0GHouDgIIJ28T10cdMX4Sg4OAgdBEHrKDgpoqDW0UE3R3Gowz2Fs1Rb08RD+X0gpL3kuCfPlVyaHk39d5tF/DOgMonDOddXv+3lnBsE5oDb386H9cUh/v78QcO2pPlo2tffxLYMVOzYy8CK87PtRoFx4KLNdn8ch3NuGv+Mr0KtVnsJypv2T4ZxDAVvC/h/QAH/zT1v8fQBeT5/m081DotlAj854DwoSzMf7SgDazarbgZ4tYujNPMhScWeQfHVAizhr1jegGfg2MqHgaNgvwXgHn81VQrKx/AnnCpwAHQnjGMAOAUebN1v5TlgN9hvBHgCuhrqV4Br/El3D+jNKg5g1tq6snUxRj6AVeAduAyWqU7z0ayv8bf4Cva6x46tasc6FtQtWb07YL7Dz2arOE7sM1s/9nKr/skojm3gxto7AyaDuhuWpyqwnmUc9n4L2Gmol3Y+9vEzN9/x544isAls2naHf9Dno7WXC+qmlg8tyRb9A4OIiET312/TiYjIP6DBSEREotNgJCIi0WkwEhGR6DQYiYhIdBqMREQkOg1GIiISnQYjERGJ7gPqQqYoUYIH0wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "phase_resolution = 10\n", "phase_sigma = 0.8*2*np.pi/phase_resolution\n", "\n", "xc, yc = np.meshgrid(np.linspace(-1.0, 1.0, 100), np.linspace(-1.0, 1.0, 100))\n", "target_x, target_y = 1, 0\n", "theta_grid = np.linspace(0, 2*np.pi, phase_resolution+1)[:-1]\n", "radii_squared = (xc-target_x)**2 + (yc-target_y)**2\n", "phase_similarity = np.exp(-0.5*radii_squared/phase_sigma**2)\n", "\n", "fig, ax = plt.subplots(figsize=(8, 5))\n", "\n", "kernel_value_img = ax.imshow(phase_similarity, extent=(-1, 1, -1, 1))\n", "cbar = plt.colorbar(kernel_value_img)\n", "cbar.set_label(\"Phase Kernel \\n Centered on 0 radians\", fontsize=16)\n", "ax.plot(np.cos(theta_grid), np.sin(theta_grid), linestyle=\"none\", marker=\"o\", color=\"r\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now to load in some helpers for data loading and batching from a python helper module I wrote for the speech competition." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/tim/anaconda3/lib/python3.6/site-packages/scipy/io/wavfile.py:273: WavFileWarning: Chunk (non-data) not understood, skipping it.\n", " WavFileWarning)\n" ] } ], "source": [ "command_words = [\"yes\", \"no\", \"up\", \"down\", \"left\", \"right\", \"on\", \"off\", \"stop\", \"go\"]\n", "command_encoder = tsr.preprocessing.LimitedVocabularyEncoder(\n", " vocabulary=command_words,\n", " encoding_type=\"onehot\",\n", ")\n", "\n", "preprocessor = tsr.SoundPreprocessor(\n", " input_transforms={\"audio\":lambda x, y: x},\n", " output_transforms={\n", " \"command_probs\":lambda x, y: command_encoder.encode(x[\"word\"]),\n", " },\n", " input_length=input_length,\n", " target_length=target_length,\n", " downsampling_rate=1,\n", " t_shift_max=0,\n", " silence_fraction=0.1,\n", " additive_bg_fraction=0.3,\n", " bg_mag_range=(0.01, 0.1),\n", " stretch_fraction=0,\n", " stretch_range=None,\n", " stretch_center=0.5,\n", " repitching_fraction=0,\n", " repitching_range=None,\n", " relabel_fraction=0.0,\n", " bg_noise_samples=tsr.get_background_noise_samples(),\n", " mean_subtract=True,\n", " norm_type=\"rms\",\n", " norm_damp=1e-4,\n", ")\n", "\n", "load_batch = tsr.data_loading.simple_sound_batch_loader\n", "\n", "train_labels = tsr.get_labels_df(\"train\")\n", "\n", "train_loading_kwargs = dict(\n", " labels_df=train_labels,\n", " preprocessor=preprocessor,\n", " augment=True,\n", " data_directory=os.path.join(tsr.data_dir, \"raw\", \"train\", \"audio\"),\n", " forcing_kwargs={},\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Phased Spectrogram Visualizations\n", "\n", "Because these phased spectrograms are 3D tensors they are significantly more difficult to visualize than regular spectrograms which you can just visualize as images. The easiest way to visualize the 3D tensor is to simply take slices and look at them one at a time. But it is difficult to get a sense of the global structure that way. Chain the slice visualizations together and you have a movie. \n", "\n", "With a little work we can make a nice animation in matplotlib which takes a 3D array and animates a movie of slices through it." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "plt.rc('animation', html='html5')\n", "\n", "class SlicingAnimation():\n", " \n", " def __init__(\n", " self,\n", " data,\n", " animation_axis,\n", " reduction_axis,\n", " axis_labels,\n", " extents,\n", " vmin=None,\n", " vmax=None,\n", " reducer=np.mean,\n", " figsize=(10, 6),\n", " ):\n", " assert reduction_axis != animation_axis\n", " left_over_axis ,= list(set(range(3)) - set([animation_axis, reduction_axis]))\n", " \n", " #permute the input data so we won't have to do any reindexing later\n", " axis_order = [animation_axis, reduction_axis, left_over_axis]\n", " self.axis_order = axis_order\n", " self.axis_labels = [axis_labels[ao] for ao in axis_order]\n", " self.extents = [extents[ao] for ao in axis_order]\n", " data = data.transpose(axis_order)\n", " \n", " #extract a reduced 2D form of the input data to aid with orienting ourselves\n", " reduction_axis = axis_order[reduction_axis]\n", " self.reduced_approx = reducer(data, axis=1)\n", " \n", " self.data = data\n", " \n", " if vmin is None:\n", " vmin = np.min(data)\n", " self.vmin = vmin\n", " if vmax is None:\n", " vmax = np.max(data)\n", " self.vmax = vmax\n", " \n", " fig, axes = plt.subplots(2, 1, figsize=figsize)\n", " self.fig = fig\n", " self.axes = axes\n", " \n", " def plot_init(self):\n", " im_kwargs = {\n", " \"vmin\":self.vmin,\n", " \"vmax\":self.vmax,\n", " \"aspect\":\"auto\",\n", " \"origin\":\"lower\",\n", " }\n", " self.axes[0].set_xlabel(self.axis_labels[0])\n", " self.axes[0].set_ylabel(self.axis_labels[2])\n", " self.axes[1].set_xlabel(self.axis_labels[2])\n", " self.axes[1].set_ylabel(self.axis_labels[1])\n", " \n", " extent = []\n", " extent.extend(self.extents[0])\n", " extent.extend(self.extents[1])\n", " static_im = self.axes[0].imshow(\n", " self.reduced_approx.transpose(),\n", " extent=extent,\n", " **im_kwargs\n", " )\n", " extent = []\n", " extent.extend(self.extents[2])\n", " extent.extend(self.extents[1])\n", " changing_im = self.axes[1].imshow(\n", " self.data[0],\n", " extent=extent,\n", " **im_kwargs,\n", " )\n", " \n", " self.temporal_indicator = self.axes[0].plot([], [], alpha=0.8, c=\"r\")[0]\n", " self.static_im = static_im\n", " self.changing_im = changing_im\n", " return [static_im, changing_im, self.temporal_indicator]\n", " \n", " def update(self, frame):\n", " self.changing_im.set_data(self.data[frame])\n", " delta = self.extents[0][1]-self.extents[0][0]\n", " ind_coord = self.extents[0][0] + frame/len(self.data) * delta\n", " span = self.axes[0].get_ylim()\n", " self.temporal_indicator.set_data([ind_coord, ind_coord], span)\n", " return [self.static_im, self.changing_im, self.temporal_indicator]\n", " \n", " def get_anim(self):\n", " ani = FuncAnimation(\n", " self.fig, \n", " self.update, \n", " frames=len(self.data),\n", " init_func=self.plot_init, \n", " blit=True,\n", " )\n", " return ani" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have everything we need to visualize the phase structure in a way which hopefully will give us more insight than simply displaying a heat map of phase angle." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "time_resolution = 256\n", "size_ratio = int(target_length/time_resolution)\n", "phased_spectrogram = keras.layers.Lambda(\n", " make_phased_spectrogram_function(\n", " frame_length=4*size_ratio, \n", " frame_step=size_ratio,\n", " phase_resolution=phase_resolution,\n", " phase_mag_damp=1e-5,\n", " phase_kernel = lambda x: tf.exp(-0.5*x/phase_sigma**2),\n", " phase_only=True,\n", " )\n", ")(sound_in)\n", "\n", "phased_specgram_model = keras.models.Model(sound_in, phased_spectrogram)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "During training it will be desirable to add extra noise and/or stretch and repitch the input sounds but for the moment all that noise is just going to make it more difficult to understand the structure in our samples. So we want to turn augmentation off for these visualizations." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "train_loading_kwargs_no_aug = copy.copy(train_loading_kwargs)\n", "train_loading_kwargs_no_aug[\"augment\"] = False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I've picked an example of the word \"sheila\" and cropped the time dimension down so that we don't spend much time looking at the silence before or after. If you animate the time axis and display frequency/phase slices for each time step you get the following. In order to save on filesize of this notebook I have saved the animation out to file instead of displaying it inline. The %%capture magic in the cell below is to suppress the usual inline output so that we don't see two versions of the same file. However you can embed the animation directly in the notebook without saving it out by removing the %%capture line at the top of the cell and ending the cell simply with phase_anim.get_anim(). " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": false }, "outputs": [], "source": [ "%%capture\n", "word = \"sheila\"\n", "example_file = train_labels.index[train_labels.word.values == word][4]\n", "example_batch = load_batch([example_file], **train_loading_kwargs_no_aug)\n", "spec_out = phased_specgram_model.predict(example_batch[0])[0]\n", "\n", "print(\"word=\", word)\n", "\n", "t_keep = 0.7\n", "cur_dat = spec_out[:int(len(spec_out)*t_keep)]\n", "vmax = np.percentile(cur_dat, 99)\n", "extents = [\n", " (0, 1000*t_keep), \n", " (0, 8000),\n", " (0, 2*np.pi*(1-1.0/phase_resolution)), \n", "]\n", "phase_anim = SlicingAnimation(\n", " data=cur_dat,\n", " animation_axis=0,\n", " reduction_axis=2,\n", " axis_labels=[\"time\", \"frequency\", \"phase\"],\n", " extents=extents,\n", " vmin=0,\n", " vmax=vmax,\n", " reducer=np.max,\n", ")\n", "mpl_anim = phase_anim.get_anim()\n", "mpl_anim.save(\"../files/videos/sheila_phase_animation.mp4\")" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%HTML\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Something you might notice is that the phases tend to form little groups which move somewhat coherently over a small band of frequencies. This behavior is particularly pronounced in the regions of the word which correspond to the vowel sounds. I think these coherent phase patterns are due to resonances of the vocal tract of the speaker (also apparently sometimes called formants ). \n", "\n", "## 3D Convolutional Neural Net\n", "\n", "Even though this representation of the phase information certainly makes it more accessible than a raw phase heatmap it still is difficult to grasp completely. Fortunately for us we don't really need to be able to clearly see the structure in our data in order for a neural network to exploit it. Now that we have transformed the phase information into a spatial dimension we can process it like volumetric data by passing it in to a 3D CNN. \n", "\n", "The phase dimension is different than the time and frequency dimensions in that it wraps around on itself. The 0'th slice is close to the n-1'th slice. In order to give the neural network access to the adjacent phase slices even at the edges of the phase grid we add Lambda layers which concatenate the last theta layer to the beginning of the phase stack and the first theta layer to the end. If we then use convolutions with padding=\"valid\" then we will recover the same result we would have ahieved if the convolutional kernel wrapped around the theta dimension.\n", "\n", "I have tinkered with a number of fairly different architectures. Most of the architectures I have tried peter out around an F1 score averaged over words of around 0.93-0.94. This is the generally same sort of performance that I get from 2D CNN models in which the phase information is thrown away. That seems to indicate that either the phase information isn't adding much or that the way I have coded up the network isn't able to take advantage of it. However all of my 2D CNN models benefit from a greatly increased batch size which is something that I can't do as easily for the 3D CNN models since I run out of memory for even modestly sized batches. The fact that the 3D models do so well with a maximum batch size of 32 might be an indication that they could perform better than the 2D models if only I could train them effectively." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def concat_theta_layer(activation):\n", " #n_theta = activation.shape[3]\n", " return K.concatenate([activation[:, :, :, :1], activation, activation[:, :, :, -1:]], axis=3)\n", "\n", "def theta_concat_output_shape(input_shape):\n", " nb, nr, nz, nt, nf = input_shape\n", " return (nb, nr, nz, nt+2, nf)\n", "\n", "from keras.layers import Conv1D, Conv2D, Conv3D\n", "from keras.layers import MaxPooling1D, MaxPooling2D, MaxPooling3D\n", "from keras.layers import LeakyReLU, concatenate\n", "\n", "def make_conv(\n", " x,\n", " n_filters,\n", " kernel_size,\n", " conv_strides,\n", " leak,\n", " pool_size,\n", " pool_strides,\n", " conv_kwargs=None,\n", "):\n", " n_dim = len(x.shape)\n", " if n_dim == 3:\n", " conv = Conv1D\n", " pool = MaxPooling1D\n", " elif n_dim == 4:\n", " conv = Conv2D\n", " pool = MaxPooling2D\n", " elif n_dim == 5:\n", " conv = Conv3D\n", " pool = MaxPooling3D\n", " else:\n", " raise ValueError()\n", " \n", " default_conv_kwargs = {\n", " \"padding\":\"same\",\n", " }\n", " \n", " if conv_kwargs is None:\n", " conv_kwargs={}\n", " default_conv_kwargs.update(conv_kwargs)\n", " x = conv(\n", " n_filters, \n", " kernel_size=kernel_size,\n", " strides=conv_strides,\n", " **default_conv_kwargs\n", " )(x)\n", " if not pool_size is None:\n", " x = pool(\n", " pool_size=pool_size,\n", " strides=pool_strides,\n", " padding=\"same\",\n", " )(x)\n", " else:\n", " assert pool_strides is None\n", " x = LeakyReLU(alpha=leak)(x)\n", " return x" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "x = phased_spectrogram\n", "\n", "#add a channels dimension to the spectrogram data\n", "x = keras.layers.Lambda(lambda x:tf.expand_dims(x, axis=4))(x)\n", "\n", "x=keras.layers.Lambda(concat_theta_layer, output_shape=theta_concat_output_shape)(x)\n", "x = make_conv(x, 16, (3, 3, 3), (2, 2, 1), 0.2, None, None, conv_kwargs={\"padding\":\"valid\"})\n", "x=keras.layers.Lambda(concat_theta_layer, output_shape=theta_concat_output_shape)(x)\n", "x = make_conv(x, 32, (3, 3, 3), (2, 2, 1), 0.2, None, None, conv_kwargs={\"padding\":\"valid\"})\n", "x=keras.layers.Lambda(concat_theta_layer, output_shape=theta_concat_output_shape)(x)\n", "x = make_conv(x, 64, (3, 3, 3), 1, 0.2, None, None, conv_kwargs={\"padding\":\"valid\"})\n", "x=keras.layers.Lambda(concat_theta_layer, output_shape=theta_concat_output_shape)(x)\n", "x = make_conv(x, 128, (3, 3, 3), 1, 0.2, None, None, conv_kwargs={\"padding\":\"valid\"})\n", "\n", "#leave the concatenated theta dimensions padding=\"same\"\n", "x=keras.layers.Lambda(concat_theta_layer, output_shape=theta_concat_output_shape)(x)\n", "x = make_conv(x, 256, (3, 3, 3), 1, 0.2, None, None)\n", "\n", "#reduce the size of the phase dimension by max pooling\n", "x = make_conv(x, 128, 1, 1, 0.2, (1, 3, 4), (1, 3, 4))\n", "\n", "#stack the phase and frequency dimensions together as extra channels\n", "reduction_shape = (int(x.shape[1]), int(x.shape[2])*int(x.shape[3])*int(x.shape[4]))\n", "x = keras.layers.Reshape(reduction_shape)(x)\n", "#now the convolutions will be 1D\n", "\n", "#collapse the channels down to a manageable number again\n", "#x = make_conv(x, 128, 1, 1, 0.2, None, None)\n", "x = make_conv(x, 32, 3, 1, 0.2, None, None)\n", "timewise_features = x\n", "\n", "dense_predictor = keras.layers.Dense(12)(keras.layers.Flatten()(timewise_features))\n", "\n", "time_indep_features = keras.layers.GlobalMaxPooling1D()(timewise_features)\n", "max_pooled_predictor = keras.layers.Dense(12)(time_indep_features)\n", "\n", "\n", "x = keras.layers.Add()([dense_predictor, max_pooled_predictor])\n", "x = keras.layers.Activation(\"softmax\", name=\"command_probs\")(x)\n", "\n", "command_probs = x" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "model = keras.models.Model(sound_in, command_probs)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "__________________________________________________________________________________________________\n", "Layer (type) Output Shape Param # Connected to \n", "==================================================================================================\n", "audio (InputLayer) (None, 16384) 0 \n", "__________________________________________________________________________________________________\n", "lambda_1 (Lambda) (None, 256, 129, 10) 0 audio[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_14 (Lambda) (None, 256, 129, 10, 0 lambda_1[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_15 (Lambda) (None, 256, 129, 12, 0 lambda_14[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_13 (Conv3D) (None, 127, 64, 10, 448 lambda_15[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_15 (LeakyReLU) (None, 127, 64, 10, 0 conv3d_13[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_16 (Lambda) (None, 127, 64, 12, 0 leaky_re_lu_15[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_14 (Conv3D) (None, 63, 31, 10, 3 13856 lambda_16[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_16 (LeakyReLU) (None, 63, 31, 10, 3 0 conv3d_14[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_17 (Lambda) (None, 63, 31, 12, 3 0 leaky_re_lu_16[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_15 (Conv3D) (None, 61, 29, 10, 6 55360 lambda_17[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_17 (LeakyReLU) (None, 61, 29, 10, 6 0 conv3d_15[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_18 (Lambda) (None, 61, 29, 12, 6 0 leaky_re_lu_17[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_16 (Conv3D) (None, 59, 27, 10, 1 221312 lambda_18[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_18 (LeakyReLU) (None, 59, 27, 10, 1 0 conv3d_16[0][0] \n", "__________________________________________________________________________________________________\n", "lambda_19 (Lambda) (None, 59, 27, 12, 1 0 leaky_re_lu_18[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_17 (Conv3D) (None, 59, 27, 12, 2 884992 lambda_19[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_19 (LeakyReLU) (None, 59, 27, 12, 2 0 conv3d_17[0][0] \n", "__________________________________________________________________________________________________\n", "conv3d_18 (Conv3D) (None, 59, 27, 12, 1 32896 leaky_re_lu_19[0][0] \n", "__________________________________________________________________________________________________\n", "max_pooling3d_3 (MaxPooling3D) (None, 59, 9, 3, 128 0 conv3d_18[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_20 (LeakyReLU) (None, 59, 9, 3, 128 0 max_pooling3d_3[0][0] \n", "__________________________________________________________________________________________________\n", "reshape_3 (Reshape) (None, 59, 3456) 0 leaky_re_lu_20[0][0] \n", "__________________________________________________________________________________________________\n", "conv1d_3 (Conv1D) (None, 59, 32) 331808 reshape_3[0][0] \n", "__________________________________________________________________________________________________\n", "leaky_re_lu_21 (LeakyReLU) (None, 59, 32) 0 conv1d_3[0][0] \n", "__________________________________________________________________________________________________\n", "flatten_3 (Flatten) (None, 1888) 0 leaky_re_lu_21[0][0] \n", "__________________________________________________________________________________________________\n", "global_max_pooling1d_2 (GlobalM (None, 32) 0 leaky_re_lu_21[0][0] \n", "__________________________________________________________________________________________________\n", "dense_4 (Dense) (None, 12) 22668 flatten_3[0][0] \n", "__________________________________________________________________________________________________\n", "dense_5 (Dense) (None, 12) 396 global_max_pooling1d_2[0][0] \n", "__________________________________________________________________________________________________\n", "add_1 (Add) (None, 12) 0 dense_4[0][0] \n", " dense_5[0][0] \n", "__________________________________________________________________________________________________\n", "command_probs (Activation) (None, 12) 0 add_1[0][0] \n", "==================================================================================================\n", "Total params: 1,563,736\n", "Trainable params: 1,563,736\n", "Non-trainable params: 0\n", "__________________________________________________________________________________________________\n" ] } ], "source": [ "model.summary()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def damped_log_loss(y_true, y_pred):\n", " damp_param = 0.0001\n", " clipped = (y_pred+damp_param)/(1.0+2.0*damp_param)\n", " costs = -1.0*(y_true*K.log(clipped) + (1.0-y_true)*K.log(1.0-clipped))\n", " return K.mean(costs)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "optimizer = \"adadelta\"\n", "model.compile(loss=damped_log_loss, optimizer=optimizer)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "dev_labels = tsr.get_labels_df(\"dev\")\n", "\n", "dev_loading_kwargs = dict(\n", " labels_df=dev_labels,\n", " preprocessor=preprocessor,\n", " augment=False,\n", " data_directory=os.path.join(tsr.data_dir, \"raw\", \"train\", \"audio\"),\n", " forcing_kwargs={},\n", ")" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "beginning training\n", "\n", " epoch 1\n", "0.28770 0.21902 0.17688 0.17231 0.14983 0.14786 0.13679 0.12461 0.11239 0.10515 0.09558 0.09674 0.09450 0.08528 0.07854 0.07652 0.07930 0.07123 0.06363 0.06660 0.05967 0.06002 0.05993 0.05732 0.05041 0.05183 0.05647 0.04738 0.04740 0.04375 0.04720 0.04538 \n", "\n", "17.87 minutes per epoch\n", "{'f1': [0.90243902439024382, 0.79078694817658346, 0.83089770354906045, 0.8125, 0.82823529411764707, 0.87688984881209509, 0.80266075388026625, 0.780287474332649, 0.90109890109890112, 0.76039603960396041, 0.9440276976341605]}\n", "{'f1': 0.83911088050868798}\n", "\n", " epoch 2\n", "0.04075 0.03771 0.04208 0.04154 0.04150 0.03635 0.03854 0.03568 0.03651 0.03064 0.03522 0.03248 0.03725 0.03486 0.03090 0.03479 0.03496 0.03141 0.02840 0.02984 0.03260 0.03188 0.03131 0.02726 0.02604 0.03027 0.02663 0.02936 0.02565 0.02692 0.02270 0.02711 \n", "\n", "18.20 minutes per epoch\n", "{'f1': [0.9158878504672896, 0.85211267605633811, 0.89494163424124507, 0.90258449304174948, 0.92693110647181631, 0.92213114754098358, 0.86680761099365744, 0.87475915221579958, 0.93775933609958517, 0.82560706401766004, 0.96501838017312946]}\n", "{'f1': 0.89859458648356849}\n", "\n", " epoch 3\n", "0.02822 0.02019 0.02083 0.02347 0.02319 0.02334 0.02501 0.02180 0.02407 0.02153 0.02536 0.02193 0.02637 0.02251 0.01991 0.02181 0.02035 0.02115 0.02259 0.02057 0.02520 0.02007 0.01940 0.02008 0.02340 0.02031 0.02051 0.01626 0.01762 0.02146 0.02298 0.02135 \n", "\n", "18.21 minutes per epoch\n", "{'f1': [0.90605427974947816, 0.88336520076481828, 0.90909090909090895, 0.88888888888888884, 0.9291666666666667, 0.92217898832684819, 0.9158110882956878, 0.8927875243664718, 0.93043478260869561, 0.8565488565488566, 0.96390608573764758]}\n", "{'f1': 0.90893029736772446}\n", "\n", " epoch 4\n", "0.01966 0.01882 0.01449 0.01423 0.01529 0.01621 0.01603 0.01693 0.01711 0.01669 0.01668 0.01837 0.01734 0.01734 0.01965 0.01532 0.01606 0.01958 0.01873 0.01839 0.01948 0.01473 0.01518 0.01338 0.01723 0.01506 0.01498 0.01726 0.01800 0.01742 0.01393 0.01679 \n", "\n", "17.77 minutes per epoch\n", "{'f1': [0.93885601577909283, 0.87050359712230208, 0.91015625, 0.91840607210626191, 0.91552062868369355, 0.91368421052631565, 0.92585170340681366, 0.90980392156862744, 0.94468085106382982, 0.91089108910891081, 0.97195817490494296]}\n", "{'f1': 0.92093750129734453}\n", "\n", " epoch 5\n", "0.02020 0.01273 0.01337 0.01409 0.01150 0.01457 0.01182 0.01302 0.01603 0.01285 0.01167 0.01661 0.01403 0.01350 0.01562 0.01359 0.01410 0.01394 0.01615 0.01301 0.01015 0.01165 0.01295 0.01611 0.01344 0.01204 0.01179 0.01152 0.01013 0.01240 0.01199 0.01491 \n", "\n", "17.73 minutes per epoch\n", "{'f1': [0.95256916996047425, 0.9020332717190388, 0.92187499999999989, 0.94716242661448136, 0.9347368421052632, 0.93581780538302273, 0.93787575150300595, 0.90522243713733075, 0.94042553191489364, 0.8701298701298702, 0.97260114259064945]}\n", "{'f1': 0.92913174991436653}\n", "\n", " epoch 6\n", "0.01529 0.00947 0.00957 0.00838 0.01003 0.00867 0.01384 0.01028 0.01080 0.00992 0.01239 0.01203 0.00943 0.01124 0.01363 0.01309 0.01269 0.01199 0.01346 0.01104 0.01096 0.00937 0.00944 0.01342 0.01545 0.00954 0.01273 0.01262 0.01176 0.01625 0.01039 0.01155 \n", "\n", "18.16 minutes per epoch\n", "{'f1': [0.96484375, 0.9050279329608939, 0.90945674044265601, 0.9477756286266924, 0.94845360824742264, 0.94093686354378814, 0.92957746478873238, 0.88848263254113335, 0.94514767932489441, 0.89539748953974907, 0.97408232672686756]}\n", "{'f1': 0.93174382879480255}\n", "\n", " epoch 7\n", "0.01171 0.00804 0.00551 0.00946 0.00902 0.00850 0.01049 0.01027 0.00999 0.00933 0.00904 0.00998 0.00891 0.01084 0.00777 0.00784 0.00846 0.00896 0.01018 0.00904 0.00884 0.01197 0.01120 0.01091 0.00652 0.00757 0.01041 0.01069 0.00941 0.00887 0.01258 0.01043 \n", "\n", "18.00 minutes per epoch\n", "{'f1': [0.94339622641509435, 0.9147005444646098, 0.90671641791044788, 0.9390962671905696, 0.94583333333333341, 0.95334685598377267, 0.93861386138613856, 0.905811623246493, 0.96265560165975106, 0.90041493775933601, 0.97617924528301891]}\n", "{'f1': 0.93516044678477861}\n", "\n", " epoch 8\n", "0.01114 0.00815 0.00735 0.00808 0.00703 0.00617 0.00875 0.00645 0.00729 0.00718 0.01046 0.00805 0.00744 0.00752 0.01164 0.00732 0.00752 0.00934 0.00899 0.00813 0.00665 0.00680 0.01227 0.00780 0.00832 0.00650 0.00744 0.00580 0.00848 0.00924 0.00831 0.00947 \n", "\n", "17.62 minutes per epoch\n", "{'f1': [0.95890410958904115, 0.91009174311926611, 0.90873786407766988, 0.9390962671905696, 0.94142259414225937, 0.94758064516129037, 0.92307692307692313, 0.90097087378640772, 0.94646680942184147, 0.90000000000000002, 0.97464064508589465]}\n", "{'f1': 0.93190804315010578}\n", "\n", " epoch 9\n", "0.00764 0.00580 0.00500 0.00893 0.00416 0.00964 0.00774 0.00561 0.00545 0.00792 0.00763 0.00629 0.00492 0.00574 0.00727 0.00963 0.00668 0.00738 0.00624 0.00738 0.00787 0.00650 0.00553 0.00620 0.00870 0.00745 0.00864 0.00674 0.00742 0.00553 0.00728 0.00706 \n", "\n", "18.14 minutes per epoch\n", "{'f1': [0.95057034220532322, 0.91461100569259968, 0.90076335877862601, 0.9242424242424242, 0.93495934959349591, 0.95445544554455441, 0.93253968253968256, 0.90000000000000002, 0.94824016563146996, 0.9117043121149897, 0.97541952257149611]}\n", "{'f1': 0.93159141899224207}\n", "\n", " epoch 10\n", "0.00692 0.00445 0.00445 0.00542 0.00458 0.00855 0.00604 0.00634 0.00525 0.00531 0.00630 0.00482 0.00437 0.00738 0.00532 0.00378 0.00556 0.00740 0.00651 0.00509 0.00541 0.00647 0.00514 0.00821 0.00591 0.01020 0.00697 0.00691 0.00586 0.00931 0.00678 0.00710 \n", "\n", "17.99 minutes per epoch\n", "{'f1': [0.96153846153846156, 0.91078066914498146, 0.91461100569259968, 0.94497153700189762, 0.96280991735537202, 0.94331983805668007, 0.92720306513409956, 0.89200000000000002, 0.94042553191489364, 0.91093117408906876, 0.97597173144876337]}\n", "{'f1': 0.93496026648880171}\n" ] } ], "source": [ "import sklearn.metrics\n", "\n", "train_losses = []\n", "model_quality_results = []\n", "\n", "def do_training_iteration(): \n", " train_hashes = train_batcher.make_batch(batch_size)\n", " xi, yi, aug_dict = load_batch(\n", " train_hashes, \n", " **train_loading_kwargs,\n", " )\n", " ctime = time.time()\n", " yres = model.train_on_batch(xi, yi)\n", " train_losses.append(yres)\n", "\n", "batch_size = 32\n", "batches_per_epoch = len(train_labels)//batch_size + 1\n", "\n", "train_batcher = tsr.data_loading.BalancedBatcher([train_labels.index,])\n", "\n", "print(\"beginning training\")\n", "\n", "for epoch_idx in range(10):\n", " stime = time.time()\n", " \n", " print(\"\\n epoch\", epoch_idx+1)\n", " for iter_idx in range(batches_per_epoch):\n", " do_training_iteration()\n", " if (iter_idx % 50) == 0:\n", " print(\"{:4.5f} \".format(np.mean(train_losses[-50:])), end=\"\")\n", " print(\"\\n\")\n", " etime = time.time()\n", " print(\"{:4.2f} minutes per epoch\".format((etime-stime)/60.0))\n", " mq_res = tsr.evaluate_command_model_quality(\n", " eval_labels=dev_labels,\n", " command_model=model,\n", " metrics={\"f1\":lambda x, y: sklearn.metrics.f1_score(x, y>0.5)},\n", " batch_loader=load_batch,\n", " batch_size=batch_size,\n", " loading_kwargs=dev_loading_kwargs,\n", " )\n", " print(mq_res)\n", " avg_mq = {metric_name:np.mean(mq_res[metric_name]) for metric_name in mq_res}\n", " print(avg_mq)\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Final Thoughts\n", "\n", "Wrapping the phase convolutions around in the phase dimension so that the spots close to 0 phase interact with those close to 2$\\pi$ seemed like a good idea but it adds significant computational overhead since the entire tensor needs to get copied between the input and output of each layer. Also because the inputs can only ever have just one phase the phase dimension adds a huge amount of redundancy. If we have 10 phase slices then at most only around 1/10th of our numbers will be non-zero. \n", "\n", "Although the use of the 3D CNN is fun the time/frequency/phase structure of the speech data doesn't look enough like volumetric data to justify a using 3D CNN. The memory and computation constraints all go through the roof without much apparent benefit. I can train up a 1D CNN model with quality comparable to my 3D CNN in less than a 2% of the training time. It should be noted however that likely one of the places where incorporating the phase information would have one of the biggest advantages is in telling the difference between noise and speech. In speech the phases aren't random at all but in pure noise even noise with strange power spectra like brown or blue noise the phase information is usually still totally random. There is no pure noise in the held out dev set and so the way we are evaluating the 3D CNN does not include in a meaningful way the ability of the CNN to tell the difference between noise and speech.\n", "\n", "Phase information really may not be important for distinguishing words from each other but it clearly is still important to how we perceive sounds as humans. Because of the large amount of redundancy involved, transforming phase into position in a 3D tensor is probably not the right way to encode the phase information. Feeding a raw phase angle into a 2D convolution doesn't seem like a good way to do it either however since the value of that convolution would be dramatically changed by shifting a phase from -$\\pi$ to $\\pi$ when such a phase change represents absolutely no change in the corresponding efffective phase information. That problem could possibly be somewhat alleviated by a data augmentation strategy where we randomly add a global phase shift of between 0 and 2$\\pi$ and then normalize back to the -$\\pi$ to $\\pi$ range but there has to be a more elegant way." ] } ], "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.6" }, "nikola": { "author": "Tim Anderton", "category": "", "date": "2018-1-22", "description": "", "link": "", "slug": "audio-3DCNN", "tags": "audio, neural networks, speech recognition", "title": "3D CNN for audio data", "type": "text" } }, "nbformat": 4, "nbformat_minor": 2 }