{"cells":[{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os\nimport seaborn as sns\nimport librosa\nimport librosa.display\nfrom IPython import display\nfrom sklearn.utils.class_weight import compute_class_weight\nfrom sklearn import metrics","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.core.display import HTML\nHTML(\"\"\"\n<style>\n.output_png {\n    display: table-cell;\n    text-align: center;\n    vertical-align: middle;\n}\n</style>\n\"\"\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore dataset","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Load csv file with information about the training data","execution_count":null},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"birdcall_meta = pd.read_csv('/kaggle/input/birdsong-recognition/train.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see how many training files we have and how many columns with information about the audio files","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Dataset has %d rows and %d columns' % birdcall_meta.shape, end=\"\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.set_option('display.max_columns', 35)\nbirdcall_meta.sample(5, random_state = 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Explore the distributions of the number of files and file duration","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We want to see how many distinct bird species there are in the data set and the distribution of observations and audio clip duration over species.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print('There are %d unique bird species in the dataset' % birdcall_meta['ebird_code'].nunique(), end=\"\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the bar chart below we can see that most species have between 90 and 100 observations","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"species_count = birdcall_meta.groupby(['species']).size().reset_index()\nspecies_count['Number of audio files interval'] = pd.cut(species_count[0], np.arange(0,110,10))\nspecies_count_bins = species_count.groupby(['Number of audio files interval']).size()\nspecies_count_bins.plot(kind=\"barh\", title=\"Count of species by number of audio files\", color='green');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And below we see that most species have between 5 and 120 minutes of recordings","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"species_duration = birdcall_meta.groupby(['species']).sum()['duration'].reset_index()\nspecies_duration['duration_mins'] = np.round(species_duration['duration']/60)\nspecies_duration['Duration interval'] =  pd.cut(species_duration['duration_mins'], 10)\nspecies_duration_bins = species_duration.groupby(['Duration interval']).size()\nspecies_duration_bins.plot(kind=\"barh\", title=\"Count of species by total duration of recordings\", color='yellow');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Because we will split the audio files in chunks of 3 seconds for example, we want to see which are the species with longest and shortest total duration of recorded audio files","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"species_duration_top = (\n    species_duration\n      .sort_values('duration', ascending=False)\n      .head(10)[['species', 'duration_mins']]\n      .set_index('species')\n)\nax = species_duration_top.plot(kind=\"barh\", title=\"Top 10 species by total duration of recordings\", color='darkblue');\nax.invert_yaxis()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"species_duration_bottom = (\n    species_duration\n      .sort_values('duration', ascending=True)\n      .head(10)[['species', 'duration_mins']]\n      .set_index('species')\n)\nax = species_duration_bottom.plot(kind=\"barh\", title=\"Bottom 10 species by total duration of recordings\", color='lightblue');\nax.invert_yaxis()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Look at pitch, speed, time of day, month of year, elevation and country","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Below I plot a few more interesting features: pitch, speed, time of day, month of year, elevation and country","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pitch_count =  birdcall_meta.groupby(['pitch']).size()\npitch_count.name = 'Pitch distribution'\npitch_count.plot.pie(y='Pitch distribution', figsize=(6, 6));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"speed_count =  birdcall_meta.groupby(['speed']).size()\nspeed_count.name = 'Speed distribution'\nspeed_count.plot.pie(y='Speed distribution', figsize=(6, 6));","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I define a function to extract the time of day from various formats in the dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_hour_of_day(time):\n    time = time.lower()\n    hour = time[:time.find(':')]\n    if hour.isnumeric():\n        hour = int(hour)\n    else:\n        hour = np.nan\n        \n    if ('pm' in time) & (hour !=12):\n        hour = hour+12    \n    if ('am' in time) & (hour ==12):\n        hour = 0    \n    return hour","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta['hour_of_day'] = list(map(extract_hour_of_day, birdcall_meta['time']))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Extract month of year from date string","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta['month_of_year'] = birdcall_meta['date'].str[5:7]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pivot time and month to create a dataframe for a nice heatmap","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"time_count = pd.pivot_table(birdcall_meta, values='rating', index=['hour_of_day'],\n                    columns=['month_of_year'], aggfunc='count')\ndel time_count['00']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks like most data was collected in the morning around May and June","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(time_count);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I create a function to extract the elevation values from various formats","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_elevation(elevation):\n    elevation = elevation.replace('m', '')\n    elevation = elevation.replace('~', '')\n    elevation = elevation.replace(',', '').strip()\n    if elevation.isnumeric():\n        elevation = float(elevation)\n    else:\n        elevation = np.nan\n    return elevation","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta['elevation_clean'] = list(map(extract_elevation, birdcall_meta['elevation']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(birdcall_meta['elevation_clean'], kde=False);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"country_count = birdcall_meta.groupby(['country']).size().sort_values(ascending=False).head(10)\ncountry_count.name = 'count'\nax = country_count.plot(kind=\"barh\", title=\"Count of recordings by country\", color='darkgreen');\nax.invert_yaxis()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Melspectrogram","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Load an example file and plot melspectrogram","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ex_file = ('/kaggle/input/birdsong-recognition/train_audio'+ '/' + \n           birdcall_meta['ebird_code']+ '/' + \n           birdcall_meta['filename']).iloc[4423] #4423\nx, sr = librosa.load(ex_file)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display.Audio(data=x, rate=sr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(14, 5))\nlibrosa.display.waveplot(x, sr=sr);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Fourier transform converts an audio signal from a time and amplitude domain to a frequency and amplitude domain.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"![ft-birdcall.gif](attachment:ft-birdcall.gif)","attachments":{"ft-birdcall.gif":{"image/gif":"R0lGODlhkAHnAAAAACwAAAAAkAHnAIcAAAAAADMAAGYAAJkAAMwAAP8AKwAAKzMAK2YAK5kAK8wAK/8AVQAAVTMAVWYAVZkAVcwAVf8AgAAAgDMAgGYAgJkAgMwAgP8AqgAAqjMAqmYAqpkAqswAqv8A1QAA1TMA1WYA1ZkA1cwA1f8A/wAA/zMA/2YA/5kA/8wA//8zAAAzADMzAGYzAJkzAMwzAP8zKwAzKzMzK2YzK5kzK8wzK/8zVQAzVTMzVWYzVZkzVcwzVf8zgAAzgDMzgGYzgJkzgMwzgP8zqgAzqjMzqmYzqpkzqswzqv8z1QAz1TMz1WYz1Zkz1cwz1f8z/wAz/zMz/2Yz/5kz/8wz//9mAABmADNmAGZmAJlmAMxmAP9mKwBmKzNmK2ZmK5lmK8xmK/9mVQBmVTNmVWZmVZlmVcxmVf9mgABmgDNmgGZmgJlmgMxmgP9mqgBmqjNmqmZmqplmqsxmqv9m1QBm1TNm1WZm1Zlm1cxm1f9m/wBm/zNm/2Zm/5lm/8xm//+ZAACZADOZAGaZAJmZAMyZAP+ZKwCZKzOZK2aZK5mZK8yZK/+ZVQCZVTOZVWaZVZmZVcyZVf+ZgACZgDOZgGaZgJmZgMyZgP+ZqgCZqjOZqmaZqpmZqsyZqv+Z1QCZ1TOZ1WaZ1ZmZ1cyZ1f+Z/wCZ/zOZ/2aZ/5mZ/8yZ///MAADMADPMAGbMAJnMAMzMAP/MKwDMKzPMK2bMK5nMK8zMK//MVQDMVTPMVWbMVZnMVczMVf/MgADMgDPMgGbMgJnMgMzMgP/MqgDMqjPMqmbMqpnMqszMqv/M1QDM1TPM1WbM1ZnM1czM1f/M/wDM/zPM/2bM/5nM/8zM////AAD/ADP/AGb/AJn/AMz/AP//KwD/KzP/K2b/K5n/K8z/K///VQD/VTP/VWb/VZn/Vcz/Vf//gAD/gDP/gGb/gJn/gMz/gP//qgD/qjP/qmb/qpn/qsz/qv//1QD/1TP/1Wb/1Zn/1cz/1f///wD//zP//2b//5n//8z///8AAAAAAAAAAAAAAAAI6AD3CRxIsKDBgwgTKlzIsKHDhxAjSpxIsaLFixgzatzIsaPHjyBDihxJsqTJkyhTqlzJsqXLlzBjypxJs6bNmzhz6tzJs6fPn0CDCh1KtKjRo0iTKl3KtKnTp1CjSp1KtarVq1izat3KtatXm8QyfR1LtizDMDjMql1LVhkAAMrYyp17VcwBAGLo6t3blBgAGDHg8h1MeKjYNHnFFl7MeKckNI0jS66JeLLlyywfY97MeWTlzqBDZ5wEWbTp0xAl5UXNuvXBNKVdy249afXs26bRxMbNm3Pt3sA36w5OfPLv4sgXD08+OJP/beDHmc/NFBhHjLjBl0tXOymwGFD7osUoHn27V2VovCuLtmyZsmXXiWs3r5WYmBgx0qxnr+z9+/HElbeTMnBgR19FzgFwwxvrtbcMMf65J56Bvc2nkyYOIOCAJgdK1B1emkRTj3vu9RfKe6CQCGBwAuKkjIY4OLBhhwyhF5gPxETD34kkRtifhPEFB1tPDzjwwzI/OJAWjQfZh59+IpboH4oouscjMfAVpxpPysgY1zIycsikQJngsMKCDbanTIo+9liilesFCdyQO/2AQBoDpWEkjcp0F4MYmSzDn5QksnmiMjwiOmV7K0L33E31yAieQMTIuIx5ykgSQAxoaCLo/6L+JfqmMsVQ6WacFPJmoU0Y5lBQkWI2pMwPaVwamn0A5Nfgjz6KyuYmVUaIJZueNtpbiyrJaZCePxSkJ54NYaihA6lOpkkYAdwQqI6nHnrog6CO6h+b/okXjXy7uUTMDQeFWRCGDjQUDQ4P4JCDA65O1mcMBgAqaIlY+uhfqXAW2qOoofpnbEjQEAPNWsi6FB2YOJxLUDRFVrusjOzJCEdjNnKKpYO89qjoMmyifDCvKQfbYzQBmBTWJJlozNWqMsWQSasHxRgrQvPOuM8bCDRLmDJi5Krfp1LymCIobYIadY/kojywMguDtPM+yhBDs81YRVxQNJqkoQbYH+UAAP+0BRFtdEKaaGgxmNTuhavO/wZ7cpuG/ti31Srv/aCbUGdZ0s6KRhNWJoF2hfO7MsqIwBsnxRhKJuwShOGSCSXJ9j4xfiwX5gpue+qpigp8uo+l+jdswW0aYJImgXbN6+IOayW2QPA6kMa9d5ZEt4E3TCIQmAhYjJCGtvKOAOcNqRyVMtSB+K+43wo+rsGvr/lmyi1v7/1/Jg2jydOoU59J7lY9fryMY1j8Bqwk8WwQXPeCvbnyXHPskCYxcgAxmoI0/HxnUP1xmamCVYyWlSh7pUJYr6ZGvsN5KmDAyl49uva1quwuSRUjSJLqppH1bOxtBnEAANSQEGaNLQcI+Jn/QiJnqaQgTWlpahPCSAQ4XvnIgSWbUqKylzCVwaB85xtf104ELDZFQxnD6GBU3JcJ5hnkXii0SBXHkEKhHSRJb0iDJNrlxYG8wXcNiVEOkKQko2BuU5ngVrhAxaaqBWuIATMYuRDWMh55qmQHMAniwHWiQ6UuFBD6FONq9hSxFelzA4GXDCkCr4J0yQHNM4j99nGDAb4Pk5qUEUPoBp7hBaVPN1gBDjJRD/as7pB5jCWJ8qg6E5Wsj7U0WY+y9pGdBS5gDRycJlBEMPfg7mFLeRwclMQ/EeIrI9GQFEGepZBKEkQZ7NokQeiWydc8UyBY9EnIxBAiNSXQYNoLnNPg/2SqwqVTYK0bnKJKBT73pOiIJTGf36xGrjqWaJjveSLu0PaTiMUIkgSpFCgxMgY0DqRIogNajLopEBwAgHIGCdokC+Iu3jHzQvfJT474UzJ7jgqdUvtetyb4oOwFroc7jFDMSkK7ZQzTPQTzVIqIsdOUDbOY/WEcQXmCM3hRlCD3QuhEQKGhgVzyqAQJXc82BKgrOiCiB7HmQDK0UZhgLgba2tVJnzbBOQ6RgvI0qd+KyNZRfStl+CSJLxFFjBORq4ELfOuafOQ1Rg4FWUldCLyaSRGPgfObBdHHNBE7kGUUyVb6EIxAzqjUgcChaM5yAEblpYlhgqQe1SMnyU7nQP9ZvuePPVRp6qzmUrQqEZ0r61Eg8zlMO5KVVKCgZyJ3ilOo+TATNBMKzjqaEI1mRA1KQpmXGAKKNhYkE84lyDDGEzfoGcRnkMvXKGOEWY7YqF8peiKvvFVL8sKOpXMkXG/V+sBcrtaH7pmpXG9qyz9GMIFO08ShoAa1BgbzPQ1jHIR80iJttnBPGZERvXwXjTT44AchQgj9RDi5hNhloQbh5tiWy5AR4mCoDimgSF3pWvYCcVF/26fL3FmoqM0zQn2DMcFOFteRdJauJjtUMCG0iSTl4A04FZ8y7DslodZjJ6uSKnOfpxFiFAlfyKXhGJBZkHttNjyWK66MlBXJ6CIltYwIgZcnL4K5M0FJl0XcY3kV2KM8/pFU/kGtO2E8uPBBcGp9VP9B+WpntQxCaL8oS9KT9eu93I7PrlbzzyJBzBIBRRMBk2IIcS+ijDRoYpn4gkNDEfyupkbS0wm5F4fcktksCsRzDAmnFqszTATyDc3yNDQRV7a3HbL41Sll82rJelJlzLZ+tQ0YjIdVKj094F4Ze1DKAhawC56oGCVatE0spCft7kOxoXaoRuim3WVWeMNevBdWD0JZgqBhNfTaaBVNbckaTgQNK8ALKFqJ61n/MaawvXWLVRtbZQBr3+E7WCIXiE42rUCQfCYrsFKHKAdkDLkPAHJPFcVrd/pWTceciYCs/BAMPYAjaizIJCYtEOQuqbkkhJvDDxKDA6ScIKH4aEL/lmndGk3CAAAYw36M2e+e0xG2ZW0trOcsLtKyc99X6/UyZFcSYgzzRLTUcYr0tCGU3SsHhdtp0k8bKp7ysGvA7Q9M6CSQaGK4IRzGyGBfyGmnyuhsgV3IJQ9CN15qCG3YbYh9DKArHb1n4LNm8bDEl6hFef3v6q1SivMNu3XSOWopTtTBDxeo3dq0P/R0T5FywCOqX7C9dvW33/I4cdcJ1SVbGkilEBARO427IjJS6iUnienIQfVieV+s0YhhLFEnBHmEfa4Y+G46OPectzAu8R2FaEt55lp7Oxy44VkacOP359cjcTqvFU+iM0bcnpqA6F4X3j1FOQ3OJ7s1cP16/5InksaZlT1w/CGiVavWfB960lBXDYLqKoNZDMaDf9pmENBlbSLnHZ5CUjlGQYmCSy81PrCGaw/4gO+VV12nTrZUb/ukCbzkEYOEW6L3dylyLw6QPUVyJNvjTv/2Z7tFV3+HSBfnHgLGaBkhBukxENHgewOBbdUkShgRd9dmSUUSadsUYQaBbRRSf+GRdgVxAzHCbgJBTQQIAJ0iR+ekYtjDMjNWdGtCMGpWMNUncIUyY1tISxMYRFSjdLzSgR2BOHYkeinicfoRTKBgbE8HTBZHOCjmIwv3HsYkYFQmEtRhPObjPxBhdrf3P+6GEEC4EDxIJgDoMAB0dgZWEMtwF//7QFFmR4QDMQl4QiihADt6tDppqGslsoVv2IDu8Wa95ngxNipA5B7/ljDpZ09syBH6lHlPxyuYBgeJtgxR9kdQ8y03pVNeV1d1dVqs1XV7tYoC9oge8X77UCYOgH0K8YgwdGUT4UI9+HILwT9JAwDZwnEDoWoI4Vjx4hxc01gTdRCZgAa0xjfhAj5slSj31Xjp1XwMCF9o+E7n11rWdyq3uBE3hjCDtwzIhnnEIn45tTdBtFqHJ2wLFzi3Qzvs1xHnVo5tFBYQcUYGGBGPxokZlXsPoTg3wHdwsUmXRFC5JwbDwDWVKBCeeCr/FYve82futUAo9ibJaFNGZ0u+gk7/Qrd48nhW7KUJmSNXtcN1dkVPmwc4PKUMJ0hxbVZIXmdoVFlfVnKMKINXbXZ6HCGNT8U1kHGRKhcvFbE5R/WIyGU00HgQDXMfxRMY2MFVUchYCJEkrycGPhA87giP+tgjreN4aNaAW5h4j7eA4PJKpXgqd8WPhvaTtTSQTXY+l+dWQ+Yx4rMJPnBsPDdkDqlStDSLqfVau8hswDUMwfduq4FcBug1+5AuGeaDFNF/C8FUZxeEBeEkeIEdVXWXDvBEJEduflllB8Au1TKTzJeYXeNz8rgMmXBRpLic/WZnLoaBQMdH7JWYxHJ/vcRn5heVy5AGsGI1XlltZCUqd8hE/2sljD2UW/3VhVUzjFNCMwNmEdIYI4phEJPQHwGYUTo4ERoikghBkpbUHf0yZmRyMfeyYB/ZJLSJg/5TVQaCCZNwhqWJUmAoikNGhe3hTue3Uu3Vc4enMgOXVgwEgXwYRKSJNTITbFBHIp4iI2vUnFkXjH9HV9VHcTMWmfSlk10JCgdJn2BXM0fGmv3zAImIHtOooKcWQ5SElwqRJKZGPYERBvsHOQqWiBKaPJ2WRWLwJZ4IDUKZfIXnJic6Xu8xCWsjcEGHVoUTMPjWTrmGmEgHKkgpSKhlS+7kcD9AcegXfle1nXGoKH9EevuWjLPoZ2u1iovJJotjlg3xG1TqEP8ciQYGckbeyRBSqHeGSD0htZ8PQQxvcDYPYZunNoBOFQbOgWYUVE8gSoskEp2A+ZyQd3RCVorpFHmy+l6iuFNJORLmE4eL+UdwACvNN4zwUis9FWRl+jQ4eZAmcniKsnBeCa36SD1S1BBDMpwLsT6xmYkRChHeqpsXIyMQElKYehKbw1EClBCTcABiAJ1qomxqKHT90aInkgYemnheWCIt6pNvCnSkSHSI11YKpAmPEhJzxaclgmmEZq+JRC9HQq2IdlOH968taK/jJWzQijIQsnxfuT5D5RztShGToA+a4HJcSndyAxE5EAMB0C80uBHA12VoiRBoEEfLcAO1Ql//BbOnKIowiRINbDqvjZq0eqhLNomi/pSBsUWUPeQjlJkRiPNsCvkedoJ1IKgyUDOVAOtDiLZEzWh+waJT44JjjYpIxqSQpyJtCSEJ9gJmElEPB1Avb7kQHhli6WEAB7CpJXFQA+F6CgGAJBIN1AEKQMsyAFdLiIcGeIFAq1N9RqliTxs1s/iv/DguqGhSeUpbNpo6DYRFMGU7SFIvhdR1bDJjEBSnIiiY8SlELNgf6AN1UIO2atKfCVFFhmgRGHIAHPKfqaZZcvdVOiOoLoG84uqNBSEGDPImIvIGMUCvJ7Y6Z/Ue/IoDZYqhyKc6TXsyPVlWaVVP4oOmwSoSBYl5/4XmcUDGdXASCuhJthFSWyC7lfSoivSFaHCQBgUCON2rPcIGnWAzcnQ7EV2SL2mgDxypEGN5ENSDAyoQA8Lbji0Re0HzoAThL6ojKDgAjwmXbxQUDUkjBuLFhWzLYuTrfPd4nYyLXrN0OiU6ZOcbEsOwlE15WmFiNZ11Pq2Do/YEh6SFYniIflbjdcbmO8RElZ8Hn9WKKLqLEGlwAz9wpBghTWVnPOAKoTlLEM6BH5JAWOTIEnGDL4uIEGKwNGzlh/Cxwwp0sE+7DJB7A0GGtWcqNWdKj9/TTkD3Xmp2SLhWtWp3U1jpPXY4ZGOABmmQBjVzIsv0AL4YkPm6Jsl4X/8TyV8mJaSBU3viN615LDWKBjappxGNKBDUE66ZxUWqB7mcAjZu8xIjVMDNWzv42h/sQR2X5x//tmxukjTa616Vq5hRw7YuRrR49kp1/LkkMayDeTI+tgxhVDtlkwZYIi0FQpVCSkR2JWz/1VLGh8lSmWlPOZFDFjhNFEw6tQxP/BrPkbf013bP5TD/eS/gYSO9yRBOtsUsQQxw4KRm/Ly2+ibhIQbau7gsJcK5wh5/M5+JN766prZ1DEQS1G9FxHt6SqZaVyqCammW1kNpwCCakFT125zEZEepu1N6BT6XdVWV5nCCLLoq6tKNgxDSWBGPeEn8LHJcQ3VpgB848JL/DhE0LTsTGhw4DphA78E1OqO4fmiwy5A0MVDC+aiTc1pr+ci0U1LJ2zmByCysGJt1yvUAtAJQXlk2IeKnCZO6nMvEi3lO7oTR/XF1niKoDIJ8L3WaJiWpA5GRHAFDr+fALZcrovoQAcoTOHChVWOGERKKPyIoHAide3oo85IrhUmLMzaLQuYrEog+pkhB4esf0ABAswNQfoTDvrPIP3JTmWBppCukrpOVPxSkl4dXqbsMvWgi9MJ5tHvLEDm240kP8LqwF6GpC0EM6XEAB6C7stlhqqoTWPp4fsxwZ/pEZdJgbZI0BjAiKSZnP7nHwByCa625jhewwbIJ0BkGgpSx/7wiaJYmgsFiaT5mKMl4JSd8ky61uraDvbDSeZLifAIb37wCLHgtEHpts0zoVGIAAwGAyDLHSn3Scd3VE6skoq51piRyxpUWUElzAKjDmPpm3+u12Uu7uXlMkSajsEj0w4hGgjqVyOs9rWkgaEZ9ikokzOmXjFRS0ldnS41cIKnreG5bOOnMn8l9EYTrVN0hr9jRqQMhFpOQpaQ8rjrhL/rm0K4VIfyBA05XJgCgAuKFtfyGx5BpzPpIzNxMsGncHupoQWo7LA6HddBpaak9TAv3BvSCA4TWROapoicsfSdcrEXyBojkH48kJYPXOhmkvgKcEANOkDU0CTcAVoE9YY0GAQ1xIQZUTEY16xLNfadZ6JWHNiU/8tliEI71uk92JKJESXCOO6cR7U5yqgyijkRrdSLwAmRk7R6KnCgf7XCHXTBXGcBNnIF2ZYc6jJD4Qqx5OJ+Koh+/7RHLgAPHiR+SelmAS8qUDtygc1U+MdRw9svK93VAAgAHQN0faqt7g9jvi3i8uoAVaOY7lOb/wFYlf1bIk8Agp7XRXjdCw5SMFZejy8iKWreYPzCEPlLIPIl8TEwwQV4QboDtCIItLWczQdPXCDFAwUUQlzV/NGHnaUzMtDjfRM1DN7DlTb0zNWNOC+3cs9YmOMmchKNbzulHDt8R+hTobPKUt/4jijw+O46mymi24yuw1lwMNGpPwYRp0Jo9dE2fC08QuwMRCRIDb6BQXaWElroPYepRH9cTUq6GyIevSlc4gQEArgRA6Nw3MwZn1rlSvtqrGxpnHvpaAfcGM88Rc6Xr91IE0Nm+yUjWcaYh7cvEDdjatA7qO/ZPDgdk7CSo1dxeCze2BMMgil739owGBrAC6yoQ/zRHWEGj8de4DNCQFwWOE9yuofn4VnX2MmN/Upgznmr6gKLSkzYZU5knfbOKF4K5k+/xjrLeH08HsZW2i/6x0XE46HsEIZGp2VuJ1rwG/LTzdHIIslCXMIPHJk3fiUNeI9SRoO2CAK9HNMGJEQOEAMDL9fkrLCH6im3sJssAA2/xRG7CHgB4OY7rqo+KUlBHmCLLryTMeAChbBmoZZnE7EOYUOFChg0dIsykaeAygcrSPHAQSlMaicQISsyURtlIZTkc/AhFTBmoYspCVdS0MRNFUMRejqT4sqWyTcp+OMixMc1QjiVPVhRIkKDKl8uIUVwm8iGagw8bKssUQ2smq/+aHCDQpNCrAzhWzTJ8cyANGmJn3b6FuxDHTKQDK0K9W7EYwbouR/ZUFi0AAAB4kQamiENTxLt8JfrNWRdUX8NQ+RqONgkADrsUIefVVDXu6ITDJrd8qQljDmUbcZJctjH2smU/WRt+OlloUc8pXz5diRfjgzSZBC5OA+oHApQ2W3vm6/cjVmUPJ6F5i1UMgBuTorn96UDT97FjSFtd5sDByH1iqp+HHzeMpqR4e1e8ObFv07wUBwOoL6/XYuOMt/4ksy/AzvpKSpmWlkGDu3oSdHDATLCLL66IaDpOvZmKo2mgvZRBQ6LJLnIAlJpeeukjtgriSMTGIKMvFAceUIz/tslgRPExpHqaTRngBJLKIarMgoYYMbQqEi6TgMJBvRy+y5AhjMLapzplRKuyS4UUk0zBHQW0LEH+loGBsMsgI6gyaCYRKQ37dlyTItQqzItF/PKKZrsYoqnoMr6aMshLt4ahz6mXfsCoJ95U8iwq+gSCQz0ghbQTlDE02bE46XgSCNJEx3pjIFGj+umByVJSqsGJLkuDHuswbEjJGABA473zwlNvDCoNRSiNkxZqqzhgu5zPPuCU1Y9Np1zFqU3CAIxUr+D6A1SMmOhKkK/oIqvswKaWwQGAGMA9UKA3uDyWoUw88suoH1oTqcWJQBJpIE0QeGCmVl8KacyN6Bsy/8To1LhxUntHetJHp0zNSVQidWXoyIVIjMEA97xsjb52E/IKgV8T0qetjT8e7YZEnYXWs0urLbMuAwh7TL86BURMGRwsWvnMA/XjL6mCAwuDMMroFGiZ61BuSBOmKL0xjWhCeljIazfacSw1RgLuKd4SLYpFjyjKzSXbpnOKr0oRgMOlsUf6reUmGVoaIWUmuWEFHIZhum/1sLyqOnb7digGaI9OcEGf9aNo2gEFvCk/O3EqSIw+84T78Jg3J+jWAACNrD+CCiU8IcZ2Ck8iEFd1ybOBaTIJJb9SizHzgGH2C7jYAQs9tqirZSlz/l6ahOKFjswEhxi6g6Z0poX94f+sti50viGVEz+uspvxCp7PaWvMXMG8Bq0oGtqK0/Tnu5yrcMHJo7mV2vUZF4j66k2b/cYpiXk0VMjCNsptWlM/kbjtOGlQSU9YQieKjEVOfnkMYDDyAxEh7UEeGVFUjKeQSRhgK/WoXt/GskGHmIwrIdyH4fYUtBAhpSmCwkt+BDKtZf2MfOHyTNKIgYaozKSF9EOX+5YxLXSxkEiDY9qGRkIJ9bQNawaDmEWMIxBGiYd8RQHMTQzUm7LB4UqpWYzCSoKjSNUkMt+a2EM0Ya4TopBpy0nDaJaxjxuQ8GMxqNlrxEe/CqExhjMkDLcc5K1qNSZ0i4NGaJahLTJpbk//gYnfmvZTP1qVDlEuAcVqCPiXAXVKTgR5Q8Ii9UShUeQNnwRVDluSOiIRpWapUgpw+NKSWL5kbgyhRxq0kis3fmwsI3sLCJVUOsO572YyBGLL7vI9QPkRMuzD3g8FAqhMlEoMbbrWHxNHQ/HppCJpQCLKMrGhTGAkDbZcWQWFRDuredEBbxgkvVQymfq4xjKgSo1JVLWMU27SJcIST+aWsqoaKs2OF1tSDIrXS2CZJI7wqQ4aJsG0GKzpZdC6E+bumUFlTMs4ZmqfIcPHR2RSBIEhSdqeyli+71XLiJmpJOE0kSvVJGwlWHzcTVyTGvWoYUwgelj3BjZAmL1Bk5oY/8M8+Zcoo5KFbDFrkBYPWis0BIB5DK3SGxDggDlmqDo3aMuxVKg99c0Id/ZZI83AhcY6OSeaKu2k5cQgyEH+sRgejcbQgKa0cH5sGGlQQ3jmFZt8BUc3ZEPncVZTNo4s7H8+rCt9QJGGKFGwmpAJyUsq5YAHvoR3JvKfBs9DjLzFoI1YJU3sgAWNGAALj5tr0F5+eDPZssxcACjsJO9DpujoaU72IeyWJPJWj661kAWJad/4F56gjCQNk4jig6Iz1JU8yTEFRE3X+CnZbLamipwykEYQ2JrV0DNremGdaOOTFQCcDLVwUcZXymIoEFosQ8Uc0+I6Wcbe6TE2ADgAAP+4tb3PDM9/RvwMmaCx3QtZZlmguK3U1geZZURDEzHoa7s0gYM05AAHP9AEOgUVL/3Ihoo4gkMD06BA6EwkYA8SEkHggAOMyNNUzi1VTXPgQ2iSeGy3hE8aCOPe95pls/P62FzjU9ED5RFdMN4JSD0T4d+SOF1vbdYffxYnyP5mMwCQBG0+Q5vQTKu11RtnGt4Ah+LRCyYVcU5uYOSXUBJHXcVr4ERU4mZTNYV/+ush/VwDChrnQCR8cc6o8ALk+GBMoVMt8j4AKrKPVefMozkXZVwlScbB0JBrvAEAvBPa1+Qnv5vO1OSAFq2ztoaH8Gxge3EbqQqXeVprmWj1PCL/lJlQtykm6t1QYzNBk4rxWffaom6qGKfHOKcmMZpgcnZiYLtI5A1hPRZptwJp1GrCw12t9D5IGxeV7YnAcJVhoCKVyC+L4T4UuRSUH0ch3Kkv3WJe1hoj5G7zGWRauAoLTEP4LovkKyQSEWOolqLFSekYFDvds56AU5QF1s+6WMscVEDkgwcYWrLPihtLTpmGMXB7vTHQG98ifRXnRTTXZsm0lAus6TMt6ylr3A4aogFjjd7wbZ25oaeDyOnXZAK37d3Sv3NlPo1LYuAEgVNMQIQTSGU8xgWZRFMa5hq/XKpBaWibZ1rS1AdgXUbrxIrUTJKDN6xsKZjliGucx16w/64ctc2r+0PGyukx1XDM6BYIu8tluQDde7/qwyGkxtTzyERDyHD4N9KZbh+JhpAYHinOUMBnL9AKmjdFuNEbJrEyuI0tUZmALus20kSeCQ8nrknVoKcNkzHQMyIhXIYkYnAAttgdtcq4gUNuIOYgwjWlJB4TGnHeXuxR27+nRvXf+/sz2mwnrdfch/lIgpRojOG0hHvXKRv4mnRyMjfOmXpNf1AUhGMPa1He8JV4FsVEI3Czr8u4RyZRqjS6cUu7NDnfOxYQqpuEwKNrAQ5vihlkWopoMjohu4FoMDWRgq2IKasW4o/aOhMLy4RQUxNiiAbtYxzGqTw0G6pfo5Huiv8cAvqIKNkxDqk2u1CgkpuMnmiqHBCva2kVIomJGzk4A0I2hZkEUEAtg9g2AfS/9gir4aO3a8mLp9Av51sjIcOBergyVoMrLZORJqS+yYofcxmv47KyUICQ7+ubv5qiH/q5FRmJtyOsRPkJjhCUIYmy+esUjCgVlEqbqrHBiGiitiOJlCCgybm2SAuNI0TC6mkeOlqBlLqhR8oyxWMQqHhAcxFBofOv1zi3F7IM6YKK2ECDL4yBazOfBdMoMilB53GaBeqa/HCrkaitVkEqdcmB9TuO+vgXIjEOWmoYSXHCoMEac4I4CKI4k8K2SMuKlEtEFLoBXUqpKkucYwoQghDmMiH7k6qbuVWjn0Hxu2NTKQuLkGm5gZkQQZUKnt2iiD4xQ6ZBlDJxDLu4xRcCjp1QM6lzABzIIjg7jQHaqdaQkpcokWcBkpr4CJHIgQcAMQD6pnRKjkRkL21ZRmJKxxtwN3XrpADZI5vLDAGzxIy0FlULKZcyx3Rsje2YFm2hDTFzoUwcM9oghgN4udIxDnoCucbgORhrldiAk6HwwRtTih8StlTJseTwMbyAPaCYv3iyJ7A5RgGMhlCEASKLSNcanwPEi4uikAlcNABYoz/plntqoXhxPvXZsH+7phDMRJn7RHUJtfn/QKF2BJU1CZ5cbJGdGJEBAoV/DBBayqgYIbQrmbNjU71iQxiu8pSV0MUFabuptJtQjAFGY8zziLnj04QV2BFpekT6eYmN7ErDwErE07QR7MBpWYESyb6qOas6uUyL0AoxaEo0s7oaGZqdwJTZUanaUA8Q04T14xTogBTR80eELDbDGhqXGL06kzqXfEOFUxrAicxhwAEVMK3IrBL86o1IEQMmpAmea8JPjJCuBB25XMDoe6GkARRdMjORKEfR2c7XM88dCgAY4CXUwqCRcDKCqks9wZNm2xF9upH/DIo2ua5QCA8nKgoY842VGLseJAuMm4zJ6jxlWEzqVAjoZCMK/yUNAyA+P6KIfcCVCgOpVGoQIuHKAOvOxWu1ZQEUNFiB9AyM0Jy3/tgHUNgOmLS7d4GhScmhQKkIFlMQyXqS/wRQpOGf5Om40akXxGRAGKmptkEgSiwgRbkb58RQiGjRvKtSmMOd8hyIffivLARTCClRGNA+6ROTd9M+MfgP7vAOEH0wDBTLu7CwcplOAcSf2EAQYKuZHruJDPqIWgS73VQPoGib3Fg9zjIpbmlFA6IXUIgdKZqNqcuPCc1ShciMUOu9Sm0IDVWmt9Il7sGhwNgOULgV5Sw80CSzcCSMcQTRHCq+s2qgUIsBlXNIOXuO7mkWx2IVKY2XExmPpNlNjP/IiODYCMDcCHj5SVh8Jo6wDX56IGXYFLuQrdHT1KsIRd4LwIhcAdDBSAM7EBxwt3cblEUqUQCY0cpQvLL5RFsbMonw0u1TSeDBHl06gDB4TSRUos9QCrtsQE4yjLl8EKqBFyJZjcfYzQ47kVJBTfChp5doO4R5AD7DuDp8g3WsVrvxk4XSVE4NGp8BzWjgSsPYiT4BAFItjHeTPtpInmnRGInQvsvQETGBsVQNgBiYT+qcJzxJiRZLviBRUKKjnaq7iVTxAU6Jwzwcrx3Ei+yCkTp7Nt3Il1iKinu9WLEIgwPojkrd1qozDANbmATJhAB4RoEY1VvhVhyipi/EFcT/MMUZmbfLAJQNg4E6zVKC00+7UEARTda4nI2aEIj9k1qC8M8nMTTNgxYFGpGboJeaoo8YAREdVReLrVqFqIdhCDUlo1ANjaGC2h6kgQrzwYEITMdyCYVbscg2MSm1lYqUhNFl0Ik66cCalVzq/KhTCanpiC1IxELaSZTkQ5H/pI+YSApaChdbBbuekqLWg7fgoNbJfQtdUgEcaMqYXLltHcNJBE3YEIjHTKGSLdW7SFVRLJXTpLkB4drWyIoYkF7n3YdVzFXE5J1Scpt99Cz7jBTRE0Q8PQ5GMRAdnJyPeJUHKY7YecPkgMHIZV+4eMpdCisxSK734tRweSagAzzi/0MKwtiEUEvJ2BCDNCEMGTAOEaw6TPQtcJQBisxW2qU6RPujZFXN2jRKFfMXXFRegnoq1qnDNlSxKjoREHErsKHSBD6LLbkB0wo1qmUo62WZCuwL8RyQUBMyADGImRkycrRgIBqTmW0N5emOFMbQdjwvselTPNlDM/qXt7mulUkUOHmmqDoww+IJxCo4mxqKTLgJrp2E2RXihyCG6AwwPQ4hzT0Qr+McLCMI+CmMcjHL8QBR7RXJqpkmZaCEIg4DQJ7cNEywHjWsBcHIamk2U0Ig/jilSemJunS9UDkknDgljEgxqoEbvOWZPT4Po1MBAxDb6g1P7E23vA2eFVTZwf+gUXMpkV+FV0KysqQZCPSUSllWCGIwkZ0o5bc5jHgCEqawiVdslWIIiYgQCmAjCU+0GikFWAWll8qKiuS4ZjjbLmaOC+AbL9/T0B1JwCw8RyyktdgoGsL4k5ioMPAVCP0CrpkCQHa2iilKsDJ5Zt9S0pVMpad5OKJgDOZ80Aq52yaODG8DipAoCqG5ieYl6IudzGTK3k8sCA8EMEm4C1MUSZuksORZnklYxI9+iEvaWWKNouhgTpszZliuGmc2q3ruDGdLXLfaM014A0BdDDyrGWq0ZJleRkH+Spa5oQrzN8Iozae44sMDS8TowOVpapkGH3IOk+zij1x8mKfZCYb/7USL3FV1Zegy3mE52QjFzV9K1QeEuOt9yOu9xuu+1mu/zuu/JhmnNpRtDV68hAq9ohwLM0nCWLp+Rpy6WBaVbVHXJGy4uKR5rsuejRYynie8VTy8kCVaeg48QQpES41MIbFEcUFfjdOooN7Lps4AUGlxsedsMUvEiNdp7NrygZBbEYOYlm236N2nKm0pldLHqF+PNWU9caafXCA9kyFJ8lN++ht+gq6x+c0gHu6phIHxfGSKUNOTHI8yzbi2mgjLEIMAkwRg6u6zQENOkY57ikuz0iPENauqK2UR/bnaxMVpYx/ACM55EbbM2b/3plAAKMfEhtsd+sK5QkuXkuAz/6kcrYhtBH8LC3su4R2/3tHuPrtVdZsU4jUjPUpQmnC2KBreUhveZVg7plYnScHwyAyAc3VC1k0DFZiWuWC6C8TEdeEOTQjsGYcojVYiPA2qAzaTGlJoe2HDQLkTkIzuTd4zxDyyxR0TGCHyqRRkOf3tHbdi3OnljwwMbc4Yy95yYCEGONA8pFm4W5yOuZRWWNxH1MzdON3bzPlslrAMfUJIKeHNUPDoNBdAc0Vm1hRHd93RreavBduIWnY3QmeatIu7wwBw+y6bFWwVBaqr3qJDBfUmMpTE3g3S4VCPBxhlL5Z0wlHw1DUzK45RgG6ghHrgVW8XKzw4fk4vCIKKzuqbxzDh1eyJFo4aqZwspcmqxVoECmElC1v3Pa78wu5o5JFEitoaCJceR2d3o7w6JW7motPmrwfF0/z9FocJ6wAuParbVcR8FYtou4O8R3DTdtR69WmitbRs8DQoYjGg1XnvpXpADsa4l6rJ34Kqj6cpT0JOFDkjH5o8bocfmMlqdn9/LwFrVUe+DIcJX5tFYopHIXrgtceIJTg+rMMKDonpWrdS2oVXbVNpQIMtDpmQd49HIQWvC06cNwKZW5imeSSsB43W0S3MIwTt5M2DDR2FJujewwahp56QOo4Q7p6vHnPNo2r/DE1oy6//lnqGogeg793hFffzyrPFA46XmaXcIZvkm0dO7iEQ2voQavWMewztmwQc0JiOf3sBTLuIGBJhL5jLcI79vomK+xSe1lvETNyRaLuZz3uUqQebtw/zEegDWObGx9mQWPxvYZ2eGBsURbu1TjXg8BGiXkGv8RjLJ5yQZmwP0ljUv1hi2GbQos3IkCUQjyfo4KieRe2PhLqZcH2mgQZzrTVZ1frfZ0ys4wh4we/5ZnIoL3ewr5oRy0kQMX6UwSNZxfvqr1asyAQ46M10L2OxvqgzIT1vqnOuBUTtH8D2yn71v9iHg/08Fl7bROVfJKpzr0Ednm8AhlT3B4h9AgcSLGjw/yDChAoXMmzo8CHEiBInUoxYD9RAYpkyadKkjNgyZZuUgQpZDJSyYspChlK2kljLkDJBqVwG0mVIkssmeazo8yfQoEKHEi1q9CjSpAJDTuKoCeVKklFLQlXZcqVHmihDUoWpTJPMNJmWKS1r9izatGrXsnUYbdjGjjSVhfIKkirOrVGxRp3JUtmbnm0HEy5s+DBiszbjesT5NeTNkTZbjowJilhXyMvSCE7s+TPo0KIJf2RMdyVe1C6j3pWZNypnZaNn065t+3ZDfcuGaeL41aXKkV5H4r1rN2Zs3MqXM29+WBlcTV6v5jQZNSbq6mLJOu/u/Tv4oMo2ZopMFzNJr//nS35s6TQ8/Pjy5w+sN973adcp9XoNmZw+gAEKeBs0GpV3WmXAUcXVZk8N+CCEESJW2jDEQLMMWFDdBJYyTRUjIYghimjWfcTgdNIyejWI0YgtuviiT+RZKFNWy7yRyYcw6rgjjwkZaOIyMjXFYo9FGqljaZlEBYd0Rzr5ZIv1bAQXMftEAyWWWUI4nmxaevklmGGKOSaZZZp5Zm1dPqQmmm26uZAyacAhED3RfKXJG2k88ACbBK2mCRxp5LDnD28aeihBgyLww6AOPJCAAwkk8AACc9Zz4VdvDDoppAkgQGkOiIrqpjKTPhqpqTn8oMkPhe7zKaWSOvqAAzmoAYcEMtFUOf8qr2fu6QAcJemTEKUPqMoqpb0q++YPkzZEEkFwPHrlstWWqScCE2E7rLXdhvnDngbZyWo9CEmbALXeqptlsw4Ayuien07qqkHY0rMuvk8qs+epk8JKawKamDttvgUbCWsOOQTm0b0LYWswxDw2OpG0CKQbMcYigvvARGpE2nDG9H00cmOTVGRgn0ChUS6I7U70BsHdHagkaD8W1ZRtAOi8cyZiAEBRJjvrnDJFygCQSYh6JkAxpCAzZ/TOYnwG9c5EUwSA1LeJEUNCarKpW0E+DwQA1wMR3SXYfipkdXwTS6T0xctN8rNBXabsEkN3F4R3QdzuM/dAMdAt0NmE183/NtuHxXDD2EgfLbgYQR8tUBo6Mx745X8DIBsxPAt0tM6ZzL25QJIDwPjW+3yyguRIB7gvpdxBpKe73YVxAEE4SBIDDvugoXPZygg+9D6TG12l6bLFMMkNOg/kMwDELC4Q788PbjTSwuuMxueZCI5G66UPz7XRl26uc9ahTf65bM77zD3W+4httJoAmFz6z0ZzL3bnUguupP18F73ipWEf05sb/MoWIH4x7QF+Y04AABCGG5hMBQDgniTyt7llxA9w9NtH5/bxwQDqbDekG5psUvdBgRggfccrHuMq174YLCMMABgDDPcRNJMd8Ged613zRBONAepQg7IBnAGlhjXv/w2ueK6T38+CKMLNBY16ZZtcAItXpQCKzYAKDNADRSTBNIhBYPErXuQEtw8piq2K+6icF3s2QAAMY4pFJEgVb/BFABQQfwJsXOeG9TvqKTFrAdTjHzUnGjdCUZFrNCQA9Xg/P46tkAIBXOoIuA96AMCEfoTGCblXvEmGbChGk50+SGc03knic/eLXxenR7YYoAEad1RkF8dmND+R7nMr8KJAhrE5JMYAkh8cIfdiULbphSaX08vkBdlXj14SBIlQ+9wVGddC9u1DEz9DovNCWDyBFS9xTwKUOQ/DyG7SjYNsil/QkIYDxnWuldG0nkDEgLtc+jJztzQgAGyJA7odgP9rN7jc+jrZOeXRDYt9zCJomDlKNH4OeXQT3P8CJ7TMdQ6j3QMh3cQGSrINkn8AuNc6w5QDB/xgV6Lhpw01mlGfHcCGKSSbDQU2N8gZEKFSGykfBSJF6+0sBnbbmQquFMAPuhF4gxxhlV4YGhMNZBhXslDpBFKgjNDMbB5xqdnKkxGt7opv4xFIVGWjD3Lq0EymAlY62UJVs7FpPC490ODE+kS7orVLxPAbb7h1RroSw6oGEWtf0Uq4u3LOSrt64hNL+aKNnSoHbPXOIIWySsnuaDWeTQ8oOiLa0YqWSWp4w2lTi1rU/kANrVJVwmA7KH79KgFpiOtouODPn5CNszrsWgYCOMUpeUXKUbM6brEcxalHMZdfs3Lur2gr3WJ1xrfWTQttnzvd7XKXVtDt7ne1G93LXre8SdEHpMC7p0lFKl7R9e5z29up5LLXVNMtrhpkZ979JuW9xn0ADuA7qIQxqlU/AIKqGAWENPxATqdlEoQBlQk4MEkTabBve8nL3w3TB1ztzYEtOSwho11slxAD7p7e0KNMRLCXafiiQagpqgGKs3gRe8APcBuhoEkiToRr4kEoeCg0iEE2bmye69aXL74VKQa4I0gMWtlNlwpmEjA2EydvIIYtWvR05QKAhkUMoN+Rkm5GSwY+fXe85TJ5kyCDgygdxRwh6F3yyYh88T6Eecndnil0INXl2OQ84viFAZKP3McbLtrHN7WRbjWusaAfZACu0bKiTgRmOd+EtEF6D39d0nKkBTQJNBDjd1USMgyDFjkaA7nPOosq3UZqx1AHyHSuE0Pv8gyAABQvoZFFVCcLEoNf03pAJkbIJJ7MK004TUrFDhGx8fjsaVO72ta+Nrazre1tc7vb3v42uMMt7nGTu9zmPje6063udbO73e5+N7xbFBAAOw=="}},"execution_count":null},{"metadata":{},"cell_type":"markdown","source":"To retain the time information we can use apply the Fourier transform over widows of the audio signal. The result can be visualized in a spectrogram, in which one dimension represents time, the other frequency and the colors represent the amplitude.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Further the melspectrogram shows the frequency in the mel scale (a non-linear transformation of the Hz scale).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(14, 5))\ns = librosa.feature.melspectrogram(x, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\ns_db = librosa.power_to_db(s, ref=np.max)\nlibrosa.display.specshow(s_db, sr=sr);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Remove silence","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"I am following the procedure in this paper: http://ceur-ws.org/Vol-1609/16090547.pdf","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"norm_s = (s-s.min())/(s.max()-s.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.ndimage.morphology import binary_erosion,binary_dilation","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For each image keep the pixels which are 3 times higher that column and row median","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"column_medians = np.median(norm_s, axis=0)\nrow_medians = np.median(norm_s, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"filtered_spectrogram = np.greater(norm_s, column_medians*3)&np.greater(norm_s.T, row_medians*3).T*1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.specshow(filtered_spectrogram);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Apply binary erosion","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"eroded_spectrogram = binary_erosion(filtered_spectrogram)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.specshow(eroded_spectrogram);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Apply dilation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"dilated_idx = binary_dilation(eroded_spectrogram.sum(axis=0)>0,  iterations=3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(dilated_idx,'ro')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dilated_idx.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(np.round(np.interp(np.arange(x.shape[0]), np.arange(dilated_idx.shape[0])*x.shape[0]/dilated_idx.shape[0], dilated_idx)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(np.round(np.interp(np.arange(x.shape[0]), np.arange(dilated_idx.shape[0])*x.shape[0]/dilated_idx.shape[0], dilated_idx)),'ro')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(14, 5))\ns = librosa.feature.melspectrogram(x, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\ns_db = librosa.power_to_db(s[:,dilated_idx], ref=np.max)\nlibrosa.display.specshow(s_db, sr=sr);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Classification based on melspectrogram","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"To showcase the general ideas in this notebook I only pick a few species out to the total 264 species.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(0)\nsample_classes = 3\nsample_species = list(np.random.choice(birdcall_meta['ebird_code'].unique(), sample_classes, replace=False))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta_samp = birdcall_meta[(birdcall_meta['ebird_code'].isin(sample_species))]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"species_duration_samp =  birdcall_meta_samp.groupby(['species']).sum()['duration']\nspecies_duration_samp.plot.pie(y='Duration distribution', figsize=(6, 6));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta_samp['path'] = '/kaggle/input/birdsong-recognition/train_audio'+ '/' +  \\\n                            birdcall_meta_samp['ebird_code'] + '/' + \\\n                            birdcall_meta_samp['filename']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta_samp['chunks'] = np.floor(birdcall_meta_samp['duration']/3).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_meta_samp = birdcall_meta_samp[birdcall_meta_samp['chunks']>0]\nbirdcall_meta_samp = birdcall_meta_samp[birdcall_meta_samp['duration']<120]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\nle = preprocessing.LabelEncoder()\nbirdcall_meta_samp['class_code'] = le.fit_transform(birdcall_meta_samp['ebird_code'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nbirdcall_train, birdcall_test = train_test_split(birdcall_meta_samp, test_size=0.2, random_state=0, stratify=birdcall_meta_samp[['ebird_code']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"birdcall_train[['path','chunks','duration','class_code']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_size = birdcall_train.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_size","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Split audio files in 3s pieces","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sec_split = 3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"classes_size = birdcall_train['ebird_code'].nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"classes_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"obs_train = birdcall_train['chunks'].sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"obs_train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = np.zeros((obs_train, 128, 130))\nY_train = np.zeros((obs_train, classes_size))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler, MinMaxScaler\nscaler = StandardScaler()\n#minmaxscaler = MinMaxScaler()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"i=0\nfor r in birdcall_train[['path','class_code']].iterrows():\n    x, sr = librosa.load(r[1]['path'])\n    S = librosa.feature.melspectrogram(x, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\n    norm_S = (S-S.min())/(S.max()-S.min())\n    column_medians = np.median(norm_S, axis=0)\n    row_medians = np.median(norm_S, axis=1)\n    eroded_spectrogram = binary_erosion(np.greater(norm_S, column_medians*3)&np.greater(norm_S.T, row_medians*3).T*1)\n    dilated_idx = binary_dilation(eroded_spectrogram.sum(axis=0)>0,  iterations=3)\n    x =x[np.round(np.interp(np.arange(x.shape[0]),\n                            np.arange(dilated_idx.shape[0])*x.shape[0]/dilated_idx.shape[0],\n                            dilated_idx)).astype(bool)]\n    x=x[:int(np.floor(x.shape[0]/sr/sec_split)*sec_split*sr)]\n    if x.shape[0]>0:\n        for n in np.array_split(x, np.floor(x.shape[0]/sr/sec_split)):        \n            print('Loading train data [%.2f%%]\\r'% np.round(i/obs_train*100, 2), end=\"\")\n            S = librosa.feature.melspectrogram(n, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\n            S_DB_sc = scaler.fit_transform(librosa.power_to_db(S))\n    #         S_DB_mm = minmaxscaler.fit_transform(S_DB_sc)\n            X_train[i, :, :] = S_DB_sc\n            Y_train[i, r[1]['class_code']] = 1\n            i += 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"i","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = X_train[:i, :, :]\nY_train = Y_train[:i, :]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"obs_test = birdcall_test['chunks'].sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"obs_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = np.zeros((obs_test, 128, 130))\nY_test = np.zeros((obs_test, classes_size))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"j=0\nfor r in birdcall_test[['path','class_code']].iterrows():\n    x, sr = librosa.load(r[1]['path'])\n    S = librosa.feature.melspectrogram(x, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\n    norm_S = (S-S.min())/(S.max()-S.min())\n    column_medians = np.median(norm_S, axis=0)\n    row_medians = np.median(norm_S, axis=1)\n    eroded_spectrogram = binary_erosion(np.greater(norm_S, column_medians*3)&np.greater(norm_S.T, row_medians*3).T*1)\n    dilated_idx = binary_dilation(eroded_spectrogram.sum(axis=0)>0,  iterations=3)\n    x =x[np.round(np.interp(np.arange(x.shape[0]),\n                            np.arange(dilated_idx.shape[0])*x.shape[0]/dilated_idx.shape[0],\n                            dilated_idx)).astype(bool)]\n    x=x[:int(np.floor(x.shape[0]/sr/sec_split)*sec_split*sr)]\n    if x.shape[0]>0:\n        for n in np.array_split(x, np.floor(x.shape[0]/sr/sec_split)):        \n            print('Loading test data [%.2f%%]\\r'% np.round(j/obs_test*100, 2), end=\"\")\n            S = librosa.feature.melspectrogram(n, sr=sr, n_fft=1028, hop_length=512, n_mels=128)\n            S_DB_sc = scaler.fit_transform(librosa.power_to_db(S))\n    #         S_DB_mm = minmaxscaler.fit_transform(S_DB_sc)\n            X_test[j, :, :] = S_DB_sc\n            Y_test[j, r[1]['class_code']] = 1\n            j += 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"j","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = X_test[:j, :, :]\nY_test = Y_test[:j, :]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = X_train.reshape(X_train.shape[0], 128, 130, 1)\nX_test = X_test.reshape(X_test.shape[0], 128, 130, 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Train model","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Conv2D, MaxPool2D, Flatten, Dense, Dropout\nfrom keras import backend as K","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"K.clear_session()\nmodel = Sequential()\nmodel.add(Conv2D(16, (4,4), strides=(1, 1), input_shape = (128, 130, 1), padding='same', activation = 'relu'))\n\nmodel.add(MaxPool2D((4,4)))\nmodel.add(Flatten())\n\nmodel.add(Dense(32))\n\nmodel.add(Dense(classes_size, activation = 'softmax'))\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class_weights = compute_class_weight(class_weight='balanced',\n                                     classes=np.arange(classes_size),\n                                     y=np.argmax(Y_train, axis=1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class_weights_dict = {}\nfor c in np.arange(classes_size):\n    class_weights_dict[c] = class_weights[c]","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"model.compile('Adam', loss = 'categorical_crossentropy',\n              metrics = ['categorical_crossentropy'])\nmodel.fit(x = X_train, y = Y_train, \n          batch_size = 64, \n          epochs = 20, \n          validation_split=0.2,\n          class_weight=class_weights_dict)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Y_pred_test = model.predict_classes(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(metrics.confusion_matrix(np.argmax(Y_test, axis=1), Y_pred_test))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(metrics.classification_report(np.argmax(Y_test, axis=1), Y_pred_test, digits=3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The model is performing ok on a few classes with a simple CNN and data pre-processing. Extending to more classes will probably involve using a more sophisticated CNN and more data manipulation like superimposing noise or multiple birds calls.","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}