Jaccard index - 1000.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import sqlite3\n",
    "import pandas as pd\n",
    "import json\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "\n",
    "from math import factorial as fact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploration of the Jaccard Index\n",
    "\n",
    "In this part of the EDA, we will explore the importance of the Jaccard index as an indicator of similarity among songs. \n",
    "\n",
    "\n",
    "Before doing so, we have to connect to the database we have already created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new database by connecting to it\n",
    "conn = sqlite3.connect(\"spotifyDB.db\")\n",
    "cur = conn.cursor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>track</th>\n",
       "      <th>track_name</th>\n",
       "      <th>album</th>\n",
       "      <th>album_name</th>\n",
       "      <th>artist</th>\n",
       "      <th>artist_name</th>\n",
       "      <th>duration_ms</th>\n",
       "      <th>playlist_member</th>\n",
       "      <th>num_member</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4HPBYMNnpP9Yc8u51VUyaD</td>\n",
       "      <td>Mariconatico</td>\n",
       "      <td>5bZ7Jgk7YE35ULHWO3wgHP</td>\n",
       "      <td>Sincopa</td>\n",
       "      <td>07PdYoE4jVRF6Ut40GgVSP</td>\n",
       "      <td>Cartel De Santa</td>\n",
       "      <td>269333</td>\n",
       "      <td>1065,9003,10561,15769,21028,138066,147268,1567...</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3MCMHJ0p0tD0eRQC0Wg7Wu</td>\n",
       "      <td>Headshot!</td>\n",
       "      <td>3k4aYVzOplQp6MTye2KlQ0</td>\n",
       "      <td>Unbreakable</td>\n",
       "      <td>4Uc6biA3i00oRAWvIfjhlk</td>\n",
       "      <td>MyChildren MyBride</td>\n",
       "      <td>210813</td>\n",
       "      <td>3647,14323,35549,71836,74701,98691,113560,1189...</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2F950w5JasEsPPHfspO5yn</td>\n",
       "      <td>Know Your Onion - Spotify Sessions Curated by ...</td>\n",
       "      <td>5aeSTpjqwuarRBM1zPCbOY</td>\n",
       "      <td>Spotify Sessions</td>\n",
       "      <td>4LG4Bs1Gadht7TCrMytQUO</td>\n",
       "      <td>The Shins</td>\n",
       "      <td>152653</td>\n",
       "      <td>2136,5176,27542,55583,77394,85407,118872,18360...</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5GAU7SP15OkPwj2tjXlUJn</td>\n",
       "      <td>Goodfriend</td>\n",
       "      <td>2kyzu2gZ1SlkrzW6l4b7dD</td>\n",
       "      <td>Feathers &amp; Fishhooks</td>\n",
       "      <td>251UrhgNbMr15NLzQ2KyKq</td>\n",
       "      <td>Rayland Baxter</td>\n",
       "      <td>266746</td>\n",
       "      <td>3626,8788,17107,49060,52900,53173,64212,66550,...</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>7bTNGfXSJthLunBJzgOyBw</td>\n",
       "      <td>Harder</td>\n",
       "      <td>5nL2ZC8NBSjMZfgpOp9P0w</td>\n",
       "      <td>FLYTRAP</td>\n",
       "      <td>41yEdWozNYEzA2RfgYQHgr</td>\n",
       "      <td>CJ Fly</td>\n",
       "      <td>177842</td>\n",
       "      <td>17432,30221,51929,53295,74242,86976,109832,139...</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    track                                         track_name  \\\n",
       "0  4HPBYMNnpP9Yc8u51VUyaD                                       Mariconatico   \n",
       "1  3MCMHJ0p0tD0eRQC0Wg7Wu                                          Headshot!   \n",
       "2  2F950w5JasEsPPHfspO5yn  Know Your Onion - Spotify Sessions Curated by ...   \n",
       "3  5GAU7SP15OkPwj2tjXlUJn                                         Goodfriend   \n",
       "4  7bTNGfXSJthLunBJzgOyBw                                             Harder   \n",
       "\n",
       "                    album            album_name                  artist  \\\n",
       "0  5bZ7Jgk7YE35ULHWO3wgHP               Sincopa  07PdYoE4jVRF6Ut40GgVSP   \n",
       "1  3k4aYVzOplQp6MTye2KlQ0           Unbreakable  4Uc6biA3i00oRAWvIfjhlk   \n",
       "2  5aeSTpjqwuarRBM1zPCbOY      Spotify Sessions  4LG4Bs1Gadht7TCrMytQUO   \n",
       "3  2kyzu2gZ1SlkrzW6l4b7dD  Feathers & Fishhooks  251UrhgNbMr15NLzQ2KyKq   \n",
       "4  5nL2ZC8NBSjMZfgpOp9P0w               FLYTRAP  41yEdWozNYEzA2RfgYQHgr   \n",
       "\n",
       "          artist_name  duration_ms  \\\n",
       "0     Cartel De Santa       269333   \n",
       "1  MyChildren MyBride       210813   \n",
       "2           The Shins       152653   \n",
       "3      Rayland Baxter       266746   \n",
       "4              CJ Fly       177842   \n",
       "\n",
       "                                     playlist_member  num_member  \n",
       "0  1065,9003,10561,15769,21028,138066,147268,1567...          11  \n",
       "1  3647,14323,35549,71836,74701,98691,113560,1189...          13  \n",
       "2  2136,5176,27542,55583,77394,85407,118872,18360...          11  \n",
       "3  3626,8788,17107,49060,52900,53173,64212,66550,...          15  \n",
       "4  17432,30221,51929,53295,74242,86976,109832,139...          11  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read the table to ensure that everything is working\n",
    "pd.read_sql_query(\"select * from tracks WHERE num_member > 10 LIMIT 5;\", conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>track</th>\n",
       "      <th>track_name</th>\n",
       "      <th>album</th>\n",
       "      <th>album_name</th>\n",
       "      <th>artist</th>\n",
       "      <th>artist_name</th>\n",
       "      <th>duration_ms</th>\n",
       "      <th>playlist_member</th>\n",
       "      <th>num_member</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1Pm3fq1SC6lUlNVBGZi3Em</td>\n",
       "      <td>Concerning Hobbits (The Lord of the Rings)</td>\n",
       "      <td>3q8vR3PFV8kG1m1Iv8DpKq</td>\n",
       "      <td>Versus Hollywood</td>\n",
       "      <td>7zdmbPudNX4SQJXnYIuCTC</td>\n",
       "      <td>Daniel Tidwell</td>\n",
       "      <td>108532</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6BAwTzs5CWWuaaOpp6mTeL</td>\n",
       "      <td>End It Now!</td>\n",
       "      <td>4SKnp9LlgXiEeVVud74fTS</td>\n",
       "      <td>Arrows</td>\n",
       "      <td>4Xo1N4V7w6qX23OHAPcLCe</td>\n",
       "      <td>The Lonely Forest</td>\n",
       "      <td>261146</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3Ex4oMC0UeAQNAG78AuVfi</td>\n",
       "      <td>Pink City</td>\n",
       "      <td>4ukM0IdJpt2NL8F8kDebS7</td>\n",
       "      <td>A Promise</td>\n",
       "      <td>5JLqvjW3Nyom2OsRUyFsS9</td>\n",
       "      <td>Xiu Xiu</td>\n",
       "      <td>133373</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>17q3q8wMWgYB8HoqZNNMRT</td>\n",
       "      <td>So Much Love So Little Time</td>\n",
       "      <td>3jlchQb5BMfS7Cug5iceJn</td>\n",
       "      <td>Live A Little Love A Lot</td>\n",
       "      <td>5juac7bFYyLKmV0VGSyaKM</td>\n",
       "      <td>Moose</td>\n",
       "      <td>202533</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>6AAWJJ0zamfEVReLLDj4BC</td>\n",
       "      <td>Screaming</td>\n",
       "      <td>5xxnHsOdealeAwAjfJfWH6</td>\n",
       "      <td>. . . XYZ</td>\n",
       "      <td>5juac7bFYyLKmV0VGSyaKM</td>\n",
       "      <td>Moose</td>\n",
       "      <td>233733</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1KVU56xCG8SqaXRjf9n4Da</td>\n",
       "      <td>Crush The Camera</td>\n",
       "      <td>0ZvFi6ooJJOuMuMWfHAstk</td>\n",
       "      <td>10:1</td>\n",
       "      <td>2JSc53B5cQ31m0xTB7JFpG</td>\n",
       "      <td>Rogue Wave</td>\n",
       "      <td>180666</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7331X4zR3C2eaLh5iHN98k</td>\n",
       "      <td>Interference</td>\n",
       "      <td>5mmZbh9XqQaR0AEeevo5TS</td>\n",
       "      <td>Interference / Smile &amp; Gesture</td>\n",
       "      <td>3X03NfrYVLMXS2kz13t5WU</td>\n",
       "      <td>Lets Kill Janice</td>\n",
       "      <td>212661</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3O7HL3mzobCIxzErIMijxJ</td>\n",
       "      <td>Spread A Little Sunshine</td>\n",
       "      <td>3Ll91eivgGYSJUVnL1WV4h</td>\n",
       "      <td>Playground</td>\n",
       "      <td>3K6DKtqcjx8Vs5XEaTgZv7</td>\n",
       "      <td>The Truth</td>\n",
       "      <td>213466</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>5La66MFvn7AX8ZjXu2s8l8</td>\n",
       "      <td>Kalte Wut / Wenn Ich Einmal Reich Bin</td>\n",
       "      <td>6u6ISXvYPQkgzAQ170HUO0</td>\n",
       "      <td>Lucky Streik</td>\n",
       "      <td>04QGT0sr98LHIUjqPJ4RG4</td>\n",
       "      <td>Floh de Cologne</td>\n",
       "      <td>313066</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3vUTaQZGXUTqgYWNUZiNW8</td>\n",
       "      <td>Love's Been Good To Me</td>\n",
       "      <td>27MGNEcpVfyejgypNIkogw</td>\n",
       "      <td>One by One</td>\n",
       "      <td>62JorFOkIjXHHU7GMT9r77</td>\n",
       "      <td>Rod McKuen</td>\n",
       "      <td>183493</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    track                                  track_name  \\\n",
       "0  1Pm3fq1SC6lUlNVBGZi3Em  Concerning Hobbits (The Lord of the Rings)   \n",
       "1  6BAwTzs5CWWuaaOpp6mTeL                                 End It Now!   \n",
       "2  3Ex4oMC0UeAQNAG78AuVfi                                   Pink City   \n",
       "3  17q3q8wMWgYB8HoqZNNMRT                 So Much Love So Little Time   \n",
       "4  6AAWJJ0zamfEVReLLDj4BC                                   Screaming   \n",
       "5  1KVU56xCG8SqaXRjf9n4Da                            Crush The Camera   \n",
       "6  7331X4zR3C2eaLh5iHN98k                                Interference   \n",
       "7  3O7HL3mzobCIxzErIMijxJ                    Spread A Little Sunshine   \n",
       "8  5La66MFvn7AX8ZjXu2s8l8       Kalte Wut / Wenn Ich Einmal Reich Bin   \n",
       "9  3vUTaQZGXUTqgYWNUZiNW8                      Love's Been Good To Me   \n",
       "\n",
       "                    album                      album_name  \\\n",
       "0  3q8vR3PFV8kG1m1Iv8DpKq                Versus Hollywood   \n",
       "1  4SKnp9LlgXiEeVVud74fTS                          Arrows   \n",
       "2  4ukM0IdJpt2NL8F8kDebS7                       A Promise   \n",
       "3  3jlchQb5BMfS7Cug5iceJn        Live A Little Love A Lot   \n",
       "4  5xxnHsOdealeAwAjfJfWH6                       . . . XYZ   \n",
       "5  0ZvFi6ooJJOuMuMWfHAstk                            10:1   \n",
       "6  5mmZbh9XqQaR0AEeevo5TS  Interference / Smile & Gesture   \n",
       "7  3Ll91eivgGYSJUVnL1WV4h                      Playground   \n",
       "8  6u6ISXvYPQkgzAQ170HUO0                    Lucky Streik   \n",
       "9  27MGNEcpVfyejgypNIkogw                      One by One   \n",
       "\n",
       "                   artist        artist_name  duration_ms playlist_member  \\\n",
       "0  7zdmbPudNX4SQJXnYIuCTC     Daniel Tidwell       108532               2   \n",
       "1  4Xo1N4V7w6qX23OHAPcLCe  The Lonely Forest       261146               4   \n",
       "2  5JLqvjW3Nyom2OsRUyFsS9            Xiu Xiu       133373               4   \n",
       "3  5juac7bFYyLKmV0VGSyaKM              Moose       202533               4   \n",
       "4  5juac7bFYyLKmV0VGSyaKM              Moose       233733               4   \n",
       "5  2JSc53B5cQ31m0xTB7JFpG         Rogue Wave       180666               4   \n",
       "6  3X03NfrYVLMXS2kz13t5WU   Lets Kill Janice       212661               4   \n",
       "7  3K6DKtqcjx8Vs5XEaTgZv7          The Truth       213466               4   \n",
       "8  04QGT0sr98LHIUjqPJ4RG4    Floh de Cologne       313066               4   \n",
       "9  62JorFOkIjXHHU7GMT9r77         Rod McKuen       183493               4   \n",
       "\n",
       "   num_member  \n",
       "0           1  \n",
       "1           1  \n",
       "2           1  \n",
       "3           1  \n",
       "4           1  \n",
       "5           1  \n",
       "6           1  \n",
       "7           1  \n",
       "8           1  \n",
       "9           1  "
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read the table\n",
    "pd.read_sql_query(\"select * from playlists LIMIT 5;\", conn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, here we define functions that we will later use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function creates all the unique pairwise combination of tracks within a playlist\n",
    "def create_pair_from_playlist(playlist_id):\n",
    "    query = 'SELECT tracks FROM playlists WHERE playlist_id = {};'.format(playlist_id)\n",
    "    all_tracks = pd.read_sql_query(query, conn)\n",
    "    ids = all_tracks['tracks'].values[0].split(',')\n",
    "    track_pairs = [(track1, track2) for i, track1 in enumerate(ids) \n",
    "               for j, track2 in enumerate(ids) if i<j]\n",
    "    return track_pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function returns jaccard index between a pair of tracks (in tuple)\n",
    "def similarity_btw_pair(pair_tuple):\n",
    "    track_id1 = pair_tuple[0]\n",
    "    track_id2 = pair_tuple[1]\n",
    "    query = 'SELECT playlist_member, num_member FROM tracks WHERE track IN (\"{}\",\"{}\");'.format(track_id1, track_id2)\n",
    "    track_data = pd.read_sql_query(query, conn)\n",
    "    \n",
    "    mem1 = track_data['playlist_member'].values[0].split(',')\n",
    "    mem2 = track_data['playlist_member'].values[1].split(',')\n",
    "    num1 = track_data['num_member'].values[0]\n",
    "    num2 = track_data['num_member'].values[1]\n",
    "    # take intersect and union\n",
    "    intersect = len(set(mem1).intersection(mem2))\n",
    "    union = len(list(set(mem1) | set(mem1)))\n",
    "    \n",
    "    jaccard = intersect/union\n",
    "    \n",
    "    por1 = intersect/num1\n",
    "    por2 = intersect/num2\n",
    "    \n",
    "    return jaccard, intersect, union, por1, por2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jaccard index mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will pick 1000 different playlists, and then compute the mean Jaccard index from all the possible song combinations within each playlist. \n",
    "\n",
    "Then, we will use the average number of songs within a playlist to pick random songs from the database, and thus find the mean Jaccard index for 1000 sets of random songs (not bootstrapped)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we take 1000 random playlists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick 1000 random playlists\n",
    "N = 1000\n",
    "\n",
    "playlists = pd.read_sql_query('SELECT * FROM playlists ORDER BY RANDOM() LIMIT {};'.format(N), conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many songs per playlist?\n",
    "\n",
    "mean_tracks = playlists[['num_tracks']].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saving the jaccard means\n",
    "\n",
    "# ATTENTION-- this cell takes hours to run.\n",
    "\n",
    "jaccard_means = []\n",
    "\n",
    "for i, playlist in playlists.iterrows():\n",
    "    playlist_id = playlist.playlist_id\n",
    "    track_pairs = create_pair_from_playlist(playlist_id)\n",
    "    jaccard = []\n",
    "\n",
    "    for pair_tuple in track_pairs:\n",
    "        j,_,_,_,_ = similarity_btw_pair(pair_tuple)\n",
    "        jaccard.append(j)\n",
    "        \n",
    "    jaccard = np.array(jaccard)\n",
    "    jaccard_mean = np.mean(jaccard)\n",
    "    \n",
    "    jaccard_means.append(jaccard_mean)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEX5JREFUeJzt3X2MZXV9x/H3p/jQVG1Y3AER2A6Y1Yh9WHVCTImGFmsBDQ9NsBKjq5KuJpDWxKaukKhpY4LPsWmLWSsREkVsEaEBK1taJSbFuujyJCALXWVhu7uiEQ3GZuHbP+6Z9Lrcmbk79965s799v5Kbe87v/M45X87O/XDmN+fck6pCktSuX5t2AZKkyTLoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY17xrQLAFi7dm3Nzs5OuwxJOqTcfvvtP6qqmaX6rYqgn52dZdu2bdMuQ5IOKUl+MEw/h24kqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxq+LOWGkps5tvnNq+d172+qntWxoHz+glqXFLntEnOQG4CngB8BSwpao+leQo4BpgFtgJvLGqfpIkwKeAs4AngLdV1XcmU75W2jTPrCUtzzBn9PuB91TVS4FXARclORnYDNxSVeuBW7p5gDOB9d1rE3D52KuWJA1tyaCvqt3zZ+RV9TPgXuA44Bzgyq7blcC53fQ5wFXVcxtwZJJjx165JGkoBzVGn2QWeDnwLeCYqtoNvf8ZAEd33Y4DHu5bbVfXJkmagqGDPslzgWuBd1fV44t1HdBWA7a3Kcm2JNv27ds3bBmSpIM0VNAneSa9kP98VX25a94zPyTTve/t2ncBJ/Stfjzw6IHbrKotVTVXVXMzM0s+IEWStExLBn13Fc1ngXur6hN9i24ANnbTG4Hr+9rfmp5XAT+dH+KRJK28YW6YOhV4C3BXku1d2yXAZcCXklwI/BA4v1t2E71LK3fQu7zy7WOtWJJ0UJYM+qr6JoPH3QFOH9C/gItGrEuSNCbeGStJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatwwjxK8IsneJHf3tV2TZHv32jn/5Kkks0l+0bfs05MsXpK0tGEeJfg54O+Aq+YbqupP56eTfBz4aV//B6tqw7gKlCSNZphHCd6aZHbQsu7B4W8E/nC8ZUmSxmXUMfpXA3uq6oG+thOTfDfJN5K8esTtS5JGNMzQzWIuAK7um98NrKuqx5K8EvhKkpdV1eMHrphkE7AJYN26dSOWIUlayLLP6JM8A/gT4Jr5tqr6ZVU91k3fDjwIvHjQ+lW1parmqmpuZmZmuWVIkpYwytDNa4H7qmrXfEOSmSRHdNMnAeuBh0YrUZI0imEur7wa+E/gJUl2JbmwW/QmfnXYBuA1wJ1J7gD+GXhXVf14nAVLkg7OMFfdXLBA+9sGtF0LXDt6WZKkcfHOWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXGjfh+91LzZzTdOZb87L3v9VPar9nhGL0mNM+glqXEGvSQ1zqCXpMYN84SpK5LsTXJ3X9sHkzySZHv3Oqtv2fuS7Ehyf5I/nlThkqThDHNG/zngjAHtn6yqDd3rJoAkJ9N7xODLunX+Yf4ZspKk6Vgy6KvqVmDY576eA3yxqn5ZVf8N7ABOGaE+SdKIRhmjvzjJnd3Qzpqu7Tjg4b4+u7o2SdKULDfoLwdeBGwAdgMf79ozoG8N2kCSTUm2Jdm2b9++ZZYhSVrKsoK+qvZU1ZNV9RTwGf5/eGYXcEJf1+OBRxfYxpaqmququZmZmeWUIUkawrKCPsmxfbPnAfNX5NwAvCnJs5OcCKwH/mu0EiVJo1jyu26SXA2cBqxNsgv4AHBakg30hmV2Au8EqKp7knwJ+B6wH7ioqp6cTOmSpGEsGfRVdcGA5s8u0v9DwIdGKUqSND7eGStJjTPoJalxBr0kNc6gl6TGGfSS1DgfJXgImtaj7SQdmjyjl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4JYM+yRVJ9ia5u6/to0nuS3JnkuuSHNm1zyb5RZLt3evTkyxekrS0Yc7oPweccUDbVuC3q+p3ge8D7+tb9mBVbehe7xpPmZKk5Voy6KvqVuDHB7TdXFX7u9nbgOMnUJskaQzGMUb/DuCrffMnJvlukm8kefUYti9JGsFIX1Oc5FJgP/D5rmk3sK6qHkvySuArSV5WVY8PWHcTsAlg3bp1o5QhSVrEss/ok2wE3gC8uaoKoKp+WVWPddO3Aw8CLx60flVtqaq5qpqbmZlZbhmSpCUsK+iTnAG8Fzi7qp7oa59JckQ3fRKwHnhoHIVKkpZnyaGbJFcDpwFrk+wCPkDvKptnA1uTANzWXWHzGuCvk+wHngTeVVU/HrhhSdKKWDLoq+qCAc2fXaDvtcC1oxYlSRof74yVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxg0V9EmuSLI3yd19bUcl2Zrkge59TdeeJH+bZEeSO5O8YlLFS5KWNuwZ/eeAMw5o2wzcUlXrgVu6eYAz6T0UfD2wCbh89DIlScs1VNBX1a3AgQ/5Pge4spu+Eji3r/2q6rkNODLJseMoVpJ08JZ8OPgijqmq3QBVtTvJ0V37ccDDff12dW27+1dOsoneGT/r1q0boQypTbObb5zavnde9vqp7VvjN4k/xmZAWz2toWpLVc1V1dzMzMwEypAkwWhBv2d+SKZ739u17wJO6Ot3PPDoCPuRJI1glKC/AdjYTW8Eru9rf2t39c2rgJ/OD/FIklbeUGP0Sa4GTgPWJtkFfAC4DPhSkguBHwLnd91vAs4CdgBPAG8fc82SpIMwVNBX1QULLDp9QN8CLhqlKEnS+HhnrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcUM9eGSQJC8BrulrOgl4P3Ak8GfAvq79kqq6adkVSpJGsuygr6r7gQ0ASY4AHgGuo/fowE9W1cfGUqEkaSTjGro5HXiwqn4wpu1JksZkXEH/JuDqvvmLk9yZ5IokawatkGRTkm1Jtu3bt29QF0nSGIwc9EmeBZwN/FPXdDnwInrDOruBjw9ar6q2VNVcVc3NzMyMWoYkaQHjOKM/E/hOVe0BqKo9VfVkVT0FfAY4ZQz7kCQt0ziC/gL6hm2SHNu37Dzg7jHsQ5K0TMu+6gYgyW8AfwS8s6/5I0k2AAXsPGCZJGmFjRT0VfUE8PwD2t4yUkWSpLHyzlhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcSHfGHu5mN9847RIkaUme0UtS4zyjl/Q00/ptdedlr5/KflvnGb0kNc6gl6TGGfSS1DiDXpIaN/IfY5PsBH4GPAnsr6q5JEcB1wCz9J4y9caq+smo+5IkHbxxndH/QVVtqKq5bn4zcEtVrQdu6eYlSVMwqaGbc4Aru+krgXMntB9J0hLGEfQF3Jzk9iSburZjqmo3QPd+9Bj2I0lahnHcMHVqVT2a5Ghga5L7hlmp+5/CJoB169aNoQxJ0iAjn9FX1aPd+17gOuAUYE+SYwG6970D1ttSVXNVNTczMzNqGZKkBYwU9Emek+R589PA64C7gRuAjV23jcD1o+xHkrR8ow7dHANcl2R+W1+oqn9N8m3gS0kuBH4InD/ifiRJyzRS0FfVQ8DvDWh/DDh9lG1LksbDO2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcDweXtGr4UPLJ8Ixekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LhlB32SE5L8R5J7k9yT5C+69g8meSTJ9u511vjKlSQdrFHujN0PvKeqvtM9N/b2JFu7ZZ+sqo+NXp4kaVTLDvqq2g3s7qZ/luRe4LhxFSZJGo+xjNEnmQVeDnyra7o4yZ1JrkiyZhz7kCQtz8hBn+S5wLXAu6vqceBy4EXABnpn/B9fYL1NSbYl2bZv375Ry5AkLWCkoE/yTHoh//mq+jJAVe2pqier6ingM8Apg9atqi1VNVdVczMzM6OUIUlaxChX3QT4LHBvVX2ir/3Yvm7nAXcvvzxJ0qhGuermVOAtwF1JtndtlwAXJNkAFLATeOdIFUqSRjLKVTffBDJg0U3LL0eSNG7eGStJjTPoJalxPjNW0mFvWs+qhZV5Xm0TQT/NfyRJWu0cupGkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4yYW9EnOSHJ/kh1JNk9qP5KkxU0k6JMcAfw9cCZwMr3HC548iX1JkhY3qTP6U4AdVfVQVf0v8EXgnAntS5K0iEkF/XHAw33zu7o2SdIKm9SDRwY9NLx+pUOyCdjUzf48yf0TqmUha4EfrfA+D5Y1joc1js+hUOchVWM+PNJ2fmuYTpMK+l3ACX3zxwOP9neoqi3Algntf0lJtlXV3LT2PwxrHA9rHJ9DoU5rfLpJDd18G1if5MQkzwLeBNwwoX1JkhYxkTP6qtqf5GLga8ARwBVVdc8k9iVJWtzEHg5eVTcBN01q+2MwtWGjg2CN42GN43Mo1GmNB0hVLd1LknTI8isQJKlxh13QJ/nLJJVk7QLLNyZ5oHttXOHa/ibJnUm2J7k5yQsX6Pdk12d7khX9I/dB1DjN4/jRJPd1dV6X5MgF+u1Mclf337JtldY41a8SSXJ+knuSPJVkwatEpnwsh61xascyyVFJtnafh61J1izQbzKf7ao6bF70Lvn8GvADYO2A5UcBD3Xva7rpNStY32/2Tf858OkF+v18isdwyRpXwXF8HfCMbvrDwIcX6Ldz0M/BaqmR3oUMDwInAc8C7gBOXuE6Xwq8BPg6MLdIv2keyyVrnPaxBD4CbO6mNy/yMzmRz/bhdkb/SeCvOODmrT5/DGytqh9X1U+ArcAZK1VcVT3eN/scFq5zaoascdrH8eaq2t/N3kbvPo5VZcgap/5VIlV1b1Wt9M2MB2XIGqd9LM8BruymrwTOXcF9Hz5Bn+Rs4JGqumORblP/6oYkH0ryMPBm4P0LdPv1JNuS3JZkRX9gYKgap34c+7wD+OoCywq4Ocnt3Z3a07JQjavpOC5ltRzLhUz7WB5TVbsBuvejF+g3kc/2xC6vnIYk/wa8YMCiS4FL6P26vOgmBrSN9ax6sRqr6vqquhS4NMn7gIuBDwzou66qHk1yEvDvSe6qqgdXUY1TP45dn0uB/cDnF9jMqd1xPBrYmuS+qrp1FdU48ePY1bBknUOY+rFcahMD2lbsZ/IgNjORz3ZTQV9Vrx3UnuR3gBOBO5JA79fk7yQ5par+p6/rLuC0vvnj6Y37TbzGAb4A3MiAoK+qR7v3h5J8HXg5vfHH1VLj1I9j9wfgNwCnVzf4OWAb88dxb5Lr6P16P7ZwGkONS36VyDgcxL/3YtuY6rEcwsSP5WI1JtmT5Niq2p3kWGDvAtuYyGf7sBi6qaq7quroqpqtqll6/+ivOCDkofeH2tclWdP9Vfx1XduKSLK+b/Zs4L4BfdYkeXY3vRY4FfjeylQ4XI1M/zieAbwXOLuqnligz3OSPG9+uqvx7tVUI4fIV4lM+1gOadrH8gZg/uqzjcDTfguZ6Gd7pf7qvJpe9F0hAMwB/9i37B3Aju719hWu61p6H5A7gX8BjjuwRuD3gbvoXTVwF3DhaqtxFRzHHfTGY7d3r0937S8EbuqmT+qO4R3APfSGAFZVjd38WcD36Z3VrWiN3f7Po3di9EtgD/C1VXgsl6xx2scSeD5wC/BA935U174in23vjJWkxh0WQzeSdDgz6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatz/AaUMd1FzSggqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Distribution of the LOG of the jaccard index means\n",
    "\n",
    "plt.hist(np.log(jaccard_means));\n",
    "plt.title('Jaccard Index- means for indices within playlists')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For the random playlists:\n",
    "\n",
    "To find out how many random tracks we need to pick in order to explore the indices between random songs (instead of songs that exists in the same playlist), we need to compute the possible combinations of track from a 67 track playlist (67 an example mean size of a playlist), given by:\n",
    "$$ \\frac{n!}{(n-k)! k!}$$\n",
    "where n = 67 and k =2.\n",
    "\n",
    "In order to create the random pairs from which the Jaccard indices will be calculated, we will pick twice that number of tracks, split the list of them into two, and then pair the two resulting lists up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For the number of random tracks:\n",
    "\n",
    "mean_tracks_int = int(mean_tracks)\n",
    "\n",
    "# Multiply by 2 because we intend to split the list\n",
    "n_rand_tracks = int(2 *fact(mean_tracks_int)/(fact(mean_tracks_int-2)*fact(2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that splits a list into two\n",
    "\n",
    "def split_list(A):\n",
    "    half = len(A)//2\n",
    "    B = A[:half]\n",
    "    C = A[half:]\n",
    "    \n",
    "    return B, C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find random tracks, make sure there's no duplicates, split list and join as tuples (ONLY ONE)\n",
    "\n",
    "rand_tracks = pd.read_sql_query('SELECT track FROM tracks ORDER BY RANDOM() LIMIT {};'.format(n_rand_tracks), conn)\n",
    "rand_tracks = list(set(rand_tracks[['track']].values[:,0]))\n",
    "\n",
    "rand_list_a, rand_list_b = split_list(rand_tracks)\n",
    "rand_tracks_paired = list(zip(rand_list_a, rand_list_b))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute jaccard for one\n",
    "\n",
    "jaccard_temp=[]\n",
    "for pair in rand_tracks_paired:\n",
    "    j_rand,_,_,_,_ = similarity_btw_pair(pair)\n",
    "    jaccard_temp.append(j_rand)\n",
    "    \n",
    "jaccard_temp = np.array(jaccard_temp)\n",
    "mean_pls = np.mean(jaccard_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random tracks jaccard means\n",
    "\n",
    "jaccard_rand_means = []\n",
    "for i in range(N):\n",
    "    rand_tracks = pd.read_sql_query('SELECT track FROM tracks ORDER BY RANDOM() LIMIT {};'.format(n_rand_tracks), conn)\n",
    "    rand_tracks = list(set(rand_tracks[['track']].values[:,0]))\n",
    "\n",
    "    rand_list_a, rand_list_b = split_list(rand_tracks)\n",
    "    rand_tracks_paired = list(zip(rand_list_a, rand_list_b))\n",
    "    \n",
    "    jaccard_temp=[]\n",
    "    for pair in rand_tracks_paired:\n",
    "        j_rand,_,_,_,_ = similarity_btw_pair(pair)\n",
    "        jaccard_temp.append(j_rand)\n",
    "\n",
    "    jaccard_temp = np.array(jaccard_temp)\n",
    "    mean_rand_jaccard = np.mean(jaccard_temp)\n",
    "    \n",
    "    jaccard_rand_means.append(mean_rand_jaccard)\n",
    "    i += 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a16925e48>"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAG5CAYAAADLbpPTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu8VXWd//HXB1DQRJSLjoITlCSZINDREMzAW15S0TQ1SzSLmbIxtYtkN6ppfubPspz6WWYlOlbeUkkbZ7xApGiIgReERlTSo44eUVG8A5/fH3sd3MA5sLnsc9bxvJ6Px37stb7ru9b+7MUW3n7XLTITSZIklUOX9i5AkiRJbzGcSZIklYjhTJIkqUQMZ5IkSSViOJMkSSoRw5kkSVKJGM4kva1FxEkRcXsry/4xIpZGRNcN3PaiiNh/4yqUpFUZzqQOpAgDb0RE39Xa50ZERsTANq5nYPG53dryczeVzHwsM7fKzOXtXUtrinCZEfHl1dobI2LsJv6sj0XEzIh4JSKmt7B8eETcUyy/JyKGVy2LiPh+RCwuXudGRNSyrqRVGc6kjudR4PjmmYgYCmzRfuWUx4aOgHUAzwFnRcTWbfA5PwLOWX1BRGwOXA/8B7AtMAW4vmgHmAiMB3YHhgEfAf6pxnUlVTGcSR3PZcCJVfMTgEurO0RE94g4LyIei4inI+JnEbFFsWzbiLghIpoi4vliekDVutMj4rsRcUdEvBQR/736SF1rIuLQiJgTES9GxOMRMXm15XsXIzMvFMtPKtq3iIgfRMTfI2JJRNxeVe9VEfG/RfuMiHhf1fYuiYgLI+KPEfEyMC4i+kTE1KKGWcC711LvKiN/6/ruEfHJosbFEfG11bbVJSImRcTDxfIrI6J3sezYiHikOVxFxMHFd+pXy34F5gN3AmfU2H+DZOYtmXkl8GQLi8cC3YAfZebrmXkBEMC+xfIJwA8yszEznwB+AJxU47qSqhjOpI7nLmDriHhvMVJ0LJURiWrfB94DDAd2BvoD3yyWdQF+DbwT+EfgVeAnq63/ceBkYDtgc+BLNdb2MpXguA1wKPDZiBgPlfO7gP8E/h3oV9Q2t1jvPOD9wGigN/AVYEWx7D+BwUUtfwUub6HW7wE9gduBnwKvATsAnype66PF7x4RuwIXAp8EdgT6AAOq1juNysjRh4rlzxe1kJlXUAlXF0REH+CXwKczs2k96voGcEZz4FubIiS+0NprPT6z2vuA+3LVZ/7dV7Q3L7+3atm9qy1b27qSqhjOpI6pefTsAGAB8ETzguI8n88AZ2Tmc5n5EvBvwHEAmbk4M6/JzFeKZd+jEiiq/Toz/yczXwWupBKk1ikzp2fm/Zm5IjPvA35bte0TgFsy87eZ+WZRx9yI6EIlQH0hM5/IzOWZOTMzXy+2+avMfKmYnwzsHhG9qj72+sy8IzNXAG8CHwW+mZkvZ+YDVA6hrY/WvvvRwA2ZOaOo5Ru8FSChcgjva8XIUXOtR8db5+OdSmWkaDrwh8y8YX2Kysy5wH8DZ9XQ95zM3Ka11/p8bpWtgCWrtS2hEopbWr4E2Kr4Pa5rXUlVOuRJvJK4DJgBDGK1Q5pURqW2BO6pPh8b6AoQEVsC5wMHUTn/B6BnRHStOjH+f6u29wqVf1zXKSI+QOV8pd2ojDp1B64qFu8EPNzCan2BHi0tK0YGvwccU3yvFVXrNP9j/3jVKv2o/L1W3fb3Wmqv0tp337F6u5n5ckQsrur7TuDaiKgObMuB7YEnMvOFiLgKOJNKgNwQ3wRmRcT5G7j+xlgKrH7O29bAS60s3xpYmpkZEetaV1IVR86kDigz/07lwoBDgN+vtvhZKocq31c1WtIrM5tDxheBXYAPZObWwD5Fe7DxfgNMBXbKzF7Az6q2+zgtn//1LJXDkC0t+zhwBLA/0AsY2EKt1YfKmoBlVIJgs39cr2/Quqeqt1uE3D5Vyx8HDl5tlKpHcf4VxdWJn6IymnjBhhSQmQuo/HmfvbZ+EXF2VG4R0uJrQz4bmAcMi6rET+XE/3lVy3evWrb7asvWtq6kKoYzqeM6Bdg3M1+ubiwO7/0COD8itgOIiP4R8eGiS08q4e2F4vylb23CmnoCz2XmaxGxJ5Vw1exyYP+o3K6hW3Hi/vCi3l8BP4yIHSOia0TsFRHdi+29DiymMhr4b2v78GLk7/fA5IjYsjhPbMIm+m5XAx8pLmrYHPgOq/4d+jPgexHxToCI6BcRRxTTPaicF3g2lfPZ+kfE55pXLC5EmFxjHd8uttHq4cnM/LfiFiEtvlpbr9j3PaiMPnaJiB4RsVmxeDqVkcDTonLByeeL9tuK90uBM4vf2o5U/ifgkhrXlVTFcCZ1UJn5cGbObmXxWcBC4K6IeBG4hcpoGVRulbAFlRGru4CbNkU5xfvngO9ExEtUDsFdWVXvY1RG+r5I5ZYNc3lrpOVLwP3A3cWy71P5++lSKoclnwAeLOpdl89TORT5v1TCwa83/Gu9JTPnUTlv7DdURtGeBxqruvyYyqjhfxff/y7gA8Wy/wM0ZuaFxflonwD+NSIGF8t3Au6osY5HqRzWfsfGfaMWfZJKcL8Q+GAx/Yvic9+gcsHDicALVEYBxxftAD8H/kDlz/EB4MairZZ1JVWJVS+ekaTaRcQwYMZGnGTe6UXlNiZXZeZe7V2LpHLwggBJG6S4yvJjQGujd6pBZjYCBjNJKxnOJG2ox6gcOjy5vQuRpLcTD2tKkiSViBcESJIklUiHPqzZt2/fHDhwYHuXIUmStE733HPPs5m5zmfqduhwNnDgQGbP9lxkSZJUfhFR0xNLPKwpSZJUIoYzSZKkEjGcSZIklUiHPudMkqR6evPNN2lsbOS1115r71LUgfTo0YMBAwaw2WabrbtzCwxnkiS1orGxkZ49ezJw4EAior3LUQeQmSxevJjGxkYGDRq0QdvwsKYkSa147bXX6NOnj8FMNYsI+vTps1GjrYYzSZLWwmCm9bWxvxnDmSRJUol4zpkkSTUaOOnGTbq9ReccutblZ5xxBu985zs5/fTTAfjwhz/MTjvtxMUXXwzAF7/4Rfr3789xxx3HaaedxtVXX83cuXN58sknOeSQQwCYPHkyW221FV/60pfW2P7o0aOZOXPmJv1Otdpqq61YunTpeq+3aNEiPvKRj/DAAw8we/ZsLr30Ui644IJW+86cOZOPf/zjG1tum3LkTJKkkqoOTytWrODZZ59l3rx5K5fPnDmTMWPGsOOOO3L11VcDMHfuXP74xz/WtP32CmabSkNDQ6vBDCrh7De/+U0bVrRpGM4kSSqpMWPGrAxQ8+bNY7fddqNnz548//zzvP7668yfP58RI0awaNEidtttN9544w2++c1vcsUVVzB8+HCuuOIKAB588EHGjh3Lu971rlXCzFZbbQXA9OnTGTt2LEcffTRDhgzhhBNOIDPXWtsf/vAHPvCBDzBixAj2339/nn766TX6XHLJJRxxxBEcdNBB7LLLLnz7299eo8/SpUvZb7/9GDlyJEOHDuX6668H4Bvf+AY//vGPV/b72te+tkYQmz59Oh/5yEcA+NOf/sTw4cMZPnw4I0aM4KWXXmLSpEn8+c9/Zvjw4Zx//vnr3N9l4WFNSZJKascdd6Rbt2489thjzJw5k7322osnnniCO++8k169ejFs2DA233zzlf0333xzvvOd7zB79mx+8pOfAJXDmgsWLGDatGm89NJL7LLLLnz2s59d4x5cc+bMYd68eey4446MGTOGO+64g7333rvV2vbee2/uuusuIoKLL76Yc889lx/84Adr9Js1axYPPPAAW265JXvssQeHHnooDQ0NK5f36NGDa6+9lq233ppnn32WUaNGcfjhh3PKKadw1FFH8YUvfIEVK1bwu9/9jlmzZvHSSy+1WM95553HT3/6U8aMGcPSpUvp0aMH55xzDueddx433HDDeu339mY4kySpxJpHz2bOnMmZZ57JE088wcyZM+nVqxejR4+uaRuHHnoo3bt3p3v37my33XY8/fTTDBgwYJU+e+6558q24cOHs2jRorWGs8bGRo499lieeuop3njjjVbv6XXAAQfQp08fAI466ihuv/32VcJZZnL22WczY8YMunTpwhNPPMHTTz/NwIED6dOnD3PmzOHpp59mxIgR9OnTp9VwNmbMGM4880xOOOEEjjrqqDW+X0fiYU1Jkkqs+byz+++/n912241Ro0Zx5513rjzfrBbdu3dfOd21a1eWLVu2QX2q/cu//Auf//znuf/++/n5z3/e6n29Vr+txOrzl19+OU1NTdxzzz3MnTuX7bfffuW2Pv3pT3PJJZfw61//mk996lNrrWfSpElcfPHFvPrqq4waNYoFCxastX+ZGc4kSSqxMWPGcMMNN9C7d2+6du1K7969eeGFF7jzzjvZa6+91ujfs2fPVkeXNqUlS5bQv39/AKZMmdJqv5tvvpnnnnuOV199leuuu26NQLlkyRK22247NttsM6ZNm8bf//73lcuOPPJIbrrpJu6++24+/OEPr7Wehx9+mKFDh3LWWWfR0NDAggUL2mxfbGoe1pQkqUbruvVFPQwdOpRnn312ldtBDB06lKVLl9K3b981+o8bN45zzjmH4cOH89WvfrVudU2ePJljjjmG/v37M2rUKB599NEW++2999588pOfZOHChXz84x9f5ZAmwAknnMBhhx1GQ0MDw4cPZ8iQISuXbb755owbN45tttmGrl27rrWeH/3oR0ybNo2uXbuy6667cvDBB9OlSxe6devG7rvvzkknncQZZ5yx8V+8DcS6rsbYqI1HnAF8GkjgfuBkYAfgd0Bv4K/AJzPzjYjoDlwKvB9YDBybmYvWtv2GhoacPXt23eoHYHIvmLykvp8hSSql+fPn8973vre9y+iwLrnkklUuTlhfK1asYOTIkVx11VUMHjx4E1dXXy39diLinsxsaGWVlep2WDMi+gOnAQ2ZuRvQFTgO+D5wfmYOBp4HTilWOQV4PjN3Bs4v+kmSpE7owQcfZOedd2a//fbrcMFsY9X7sGY3YIuIeBPYEngK2BdoHpudAkwGLgSOKKYBrgZ+EhGR9RzakyRJdXPSSSdx0kknbdC6u+66K4888simLaiDqNvIWWY+AZwHPEYllC0B7gFeyMzmS0Aagf7FdH/g8WLdZUX/PvWqT5IkqYzqeVhzWyqjYYOAHYF3AAe30LV5ZKylR7ivMWoWERMjYnZEzG5qatpU5UqSJJVCPW+lsT/waGY2ZeabwO+B0cA2EdF8OHUA8GQx3QjsBFAs7wU8t/pGM/OizGzIzIZ+/frVsXxJkqS2V89w9hgwKiK2jMod5/YDHgSmAUcXfSYA1xfTU4t5iuW3eb6ZJEnqbOp2QUBm/iUirqZyu4xlwBzgIuBG4HcR8a9F2y+LVX4JXBYRC6mMmB1Xr9okSdogk3tt4u2t+1ZNXbt2ZejQoSvnr7vuOgYOHLhp66ij6dOnd8jnW1bb2FuCrK+6Xq2Zmd8CvrVa8yPAni30fQ04pp71SJLU0WyxxRbMnTu31eXLli2jW7f6/HNez22rdT6+SZKkDuaSSy7hmGOO4bDDDuPAAw8kM/nyl7/MbrvtxtChQ7niiiuAyqjVhz70IT72sY/xnve8h0mTJnH55Zez5557MnToUB5++OE1tj158mQmTpzIgQceyIknnsiiRYv44Ac/yMiRIxk5ciQzZ85cue2xY8dy9NFHM2TIEE444QSaz0a66aabGDJkCHvvvTe///3vV277ueeeY/z48QwbNoxRo0Zx3333rfzMCRMmcOCBBzJw4EB+//vf85WvfIWhQ4dy0EEH8eabb65R5wUXXMCuu+7KsGHDOO64ysG2WbNmMXr0aEaMGMHo0aP529/+tnJ/jR8/nsMOO4xBgwbxk5/8hB/+8IeMGDGCUaNG8dxzlVPcx44dy+mnn87o0aPZbbfdmDVr1hqf29TUxEc/+lH22GMP9thjD+64444N/nNsjXFYkqQSe/XVVxk+fDgAgwYN4tprrwXgzjvv5L777qN3795cc801zJ07l3vvvZdnn32WPfbYg3322QeAe++9l/nz59O7d2/e9a538elPf5pZs2bx4x//mH//93/nRz/60Rqfec8993D77bezxRZb8Morr3DzzTfTo0cPHnroIY4//nian84zZ84c5s2bx4477siYMWO44447aGho4DOf+Qy33XYbO++8M8cee+zK7X7rW99ixIgRXHfdddx2222ceOKJK0cFH374YaZNm8aDDz7IXnvtxTXXXMO5557LkUceyY033sj48eNXqfGcc87h0UcfpXv37rzwwgsADBkyhBkzZtCtWzduueUWzj77bK655hoAHnjgAebMmcNrr73GzjvvzPe//33mzJnDGWecwaWXXsrpp58OwMsvv8zMmTOZMWMGn/rUp3jggQdW+dwvfOELnHHGGey999489thjfPjDH2b+/Pkb94e8GsOZJEkl1tphzQMOOIDevXsDcPvtt3P88cfTtWtXtt9+ez70oQ9x9913s/XWW7PHHnuwww47APDud7+bAw88EKg8n3PatGktfubhhx/OFltsAcCbb77J5z//eebOnUvXrl35n//5n5X99txzTwYMGADA8OHDWbRoEVtttRWDBg1aeVf/T3ziE1x00UUr62wOS/vuuy+LFy9myZLKeXcHH3wwm222GUOHDmX58uUcdNBBK+tctGjRGjUOGzaME044gfHjx68MbkuWLGHChAk89NBDRMQqI27jxo2jZ8+e9OzZk169enHYYYet3H7zCB7A8ccfD8A+++zDiy++uDL4Nbvlllt48MEHV86/+OKLvPTSS/Ts2bPFfbkhDGeSJHVA73jHO1ZOr+3mBt27d1853aVLl5XzXbp0YdmyZS2uU73t888/n+233557772XFStW0KNHjxa33bVr15Xbq9ykYU0t1dnct7quzTbbbGV7a3XeeOONzJgxg6lTp/Ld736XefPm8Y1vfINx48Zx7bXXsmjRIsaOHbve+2H12lefX7FiBXfeeefK8FoPnnMmSVIHt88++3DFFVewfPlympqamDFjBnvuuca1dxtkyZIl7LDDDnTp0oXLLruM5cuXr7X/kCFDePTRR1eez/bb3/52lTovv/xyoHLOWt++fdl6663Xu6YVK1bw+OOPM27cOM4991xeeOEFli5dypIlS+jfv/LgoUsuuWS9twusPF/v9ttvp1evXvTqteoVugceeOAqV22u7WKNDeXImSRJtarh1hft4cgjj+TOO+9k9913JyI499xz+Yd/+AcWLFiw0dv+3Oc+x0c/+lGuuuoqxo0bt8qoWkt69OjBRRddxKGHHkrfvn3Ze++9V563NXnyZE4++WSGDRvGlltuyZQpUzaopuXLl/OJT3yCJUuWkJmcccYZbLPNNnzlK19hwoQJ/PCHP2TffffdoG1vu+22jB49mhdffJFf/epXayy/4IILOPXUUxk2bBjLli1jn3324Wc/+9kGfVZroiPf57WhoSGbT0qsm8m9SvsfoySpvubPn8973/ve9i5DbWTs2LGcd955NDQ0bPS2WvrtRMQ9mbnOjXtYU5IkqUQ8rClJkkTlPLgycORMkqS16Min/6h9bOxvxnAmSVIrevToweLFiw1oqllmsnjx4lVuObK+PKwpSVIrBgwYQGNjI01NTe1dijqQHj16rLw574YwnEmS1IrNNtuMQYMGtXcZ6mQ8rClJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUonULZxFxC4RMbfq9WJEnB4RvSPi5oh4qHjftugfEXFBRCyMiPsiYmS9apMkSSqruoWzzPxbZg7PzOHA+4FXgGuBScCtmTkYuLWYBzgYGFy8JgIX1qs2SZKksmqrw5r7AQ9n5t+BI4ApRfsUYHwxfQRwaVbcBWwTETu0UX2SJEml0Fbh7Djgt8X09pn5FEDxvl3R3h94vGqdxqJtFRExMSJmR8TspqamOpYsSZLU9uoeziJic+Bw4Kp1dW2hLddoyLwoMxsys6Ffv36bokRJkqTSaIuRs4OBv2bm08X8082HK4v3Z4r2RmCnqvUGAE+2QX2SJEml0Rbh7HjeOqQJMBWYUExPAK6vaj+xuGpzFLCk+fCnJElSZ9GtnhuPiC2BA4B/qmo+B7gyIk4BHgOOKdr/CBwCLKRyZefJ9axNkiSpjOoazjLzFaDPam2LqVy9uXrfBE6tZz2SJEll5xMCJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJVIXcNZRGwTEVdHxIKImB8Re0VE74i4OSIeKt63LfpGRFwQEQsj4r6IGFnP2iRJksqo3iNnPwZuyswhwO7AfGAScGtmDgZuLeYBDgYGF6+JwIV1rk2SJKl06hbOImJrYB/glwCZ+UZmvgAcAUwpuk0BxhfTRwCXZsVdwDYRsUO96pMkSSqjeo6cvQtoAn4dEXMi4uKIeAewfWY+BVC8b1f07w88XrV+Y9G2ioiYGBGzI2J2U1NTHcuXJElqe/UMZ92AkcCFmTkCeJm3DmG2JFpoyzUaMi/KzIbMbOjXr9+mqVSSJKkk6hnOGoHGzPxLMX81lbD2dPPhyuL9mar+O1WtPwB4so71SZIklU7dwllm/i/weETsUjTtBzwITAUmFG0TgOuL6anAicVVm6OAJc2HPyVJkjqLbnXe/r8Al0fE5sAjwMlUAuGVEXEK8BhwTNH3j8AhwELglaKvJElSp1LXcJaZc4GGFhbt10LfBE6tZz2SJEll5xMCJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJVIXcNZRCyKiPsjYm5EzC7aekfEzRHxUPG+bdEeEXFBRCyMiPsiYmQ9a5MkSSqjthg5G5eZwzOzoZifBNyamYOBW4t5gIOBwcVrInBhG9QmSZJUKu1xWPMIYEoxPQUYX9V+aVbcBWwTETu0Q32SJEntpt7hLIH/joh7ImJi0bZ9Zj4FULxvV7T3Bx6vWrexaFtFREyMiNkRMbupqamOpUuSJLW9bnXe/pjMfDIitgNujogFa+kbLbTlGg2ZFwEXATQ0NKyxXJIkqSOr68hZZj5ZvD8DXAvsCTzdfLiyeH+m6N4I7FS1+gDgyXrWJ0mSVDZ1C2cR8Y6I6Nk8DRwIPABMBSYU3SYA1xfTU4ETi6s2RwFLmg9/SpIkdRb1PKy5PXBtRDR/zm8y86aIuBu4MiJOAR4Djin6/xE4BFgIvAKcXMfaJEmSSqlu4SwzHwF2b6F9MbBfC+0JnFqveiRJkjoCnxAgSZJUIvW+WrNTGjjpxvYuYaMtOufQ9i5BkqROyZEzSZKkEjGcSZIklYjhTJIkqURqCmcRsVu9C5EkSVLtI2c/i4hZEfG5iNimrhVJkiR1YjWFs8zcGziByuOVZkfEbyLigLpWJkmS1AnVfM5ZZj4EfB04C/gQcEFELIiIo+pVnCRJUmdT6zlnwyLifGA+sC9wWGa+t5g+v471SZIkdSq13oT2J8AvgLMz89Xmxsx8MiK+XpfKJEmSOqFaw9khwKuZuRwgIroAPTLzlcy8rG7VSZIkdTK1nnN2C7BF1fyWRZskSZI2oVrDWY/MXNo8U0xvWZ+SJEmSOq9aw9nLETGyeSYi3g+8upb+kiRJ2gC1nnN2OnBVRDxZzO8AHFufkiRJkjqvmsJZZt4dEUOAXYAAFmTmm3WtTJIkqROqdeQMYA9gYLHOiIggMy+tS1WSJEmdVE3hLCIuA94NzAWWF80JGM4kSZI2oVpHzhqAXTMz61mMJElSZ1fr1ZoPAP9Qz0IkSZJU+8hZX+DBiJgFvN7cmJmH16UqSZKkTqrWcDa5nkVIkiSpotZbafwpIt4JDM7MWyJiS6BrfUuTJEnqfGo65ywiPgNcDfy8aOoPXFevoiRJkjqrWi8IOBUYA7wIkJkPAdvVqyhJkqTOqtZw9npmvtE8ExHdqNznTJIkSZtQreHsTxFxNrBFRBwAXAX8oX5lSZIkdU61hrNJQBNwP/BPwB+Br9erKEmSpM6q1qs1VwC/KF6SJEmqk1qfrfkoLZxjlpnv2uQVSZIkdWLr82zNZj2AY4Dem74cSZKkzq2mc84yc3HV64nM/BGwb51rkyRJ6nRqPaw5smq2C5WRtJ51qUiSJKkTq/Ww5g+qppcBi4CPbfJqJEmSOrlar9YcV+9CJEmSVPthzTPXtjwzf7iWdbsCs4EnMvMjETEI+B2VCwr+CnwyM9+IiO7ApcD7gcXAsZm5qKZvIUmS9DZR601oG4DPUnngeX/gn4FdqZx3tq5zz74AzK+a/z5wfmYOBp4HTinaTwGez8ydgfOLfpIkSZ1KreGsLzAyM7+YmV+kMro1IDO/nZnfbm2liBgAHApcXMwHlas8ry66TAHGF9NHFPMUy/cr+kuSJHUatYaH71rjAAAR8klEQVSzfwTeqJp/AxhYw3o/Ar4CrCjm+wAvZOayYr6RykgcxfvjAMXyJUX/VUTExIiYHRGzm5qaaixfkiSpY6g1nF0GzIqIyRHxLeAvVM4Pa1VEfAR4JjPvqW5uoWvWsOythsyLMrMhMxv69etXW/WSJEkdRK1Xa34vIv4T+GDRdHJmzlnHamOAwyPiECpPFdiaykjaNhHRrRgdGwA8WfRvBHYCGiOiG9ALeG69vo0kSVIHV+vIGcCWwIuZ+WMqAWrQ2jpn5lczc0BmDgSOA27LzBOAacDRRbcJwPXF9NRinmL5bZm5xsiZJEnS21lN4aw4lHkW8NWiaTPgPzbwM88CzoyIhVTOKftl0f5LoE/RfiYwaQO3L0mS1GHV+oSAI4ERVO5LRmY+GRE1P74pM6cD04vpR4A9W+jzGpUHqkuSJHVatR7WfKM4xJgAEfGO+pUkSZLUedUazq6MiJ9TOZn/M8AtwC/qV5YkSVLnVOvVmudFxAHAi8AuwDcz8+a6ViZJktQJrTOcFc/G/K/M3B8wkEmSJNXROg9rZuZy4JWI6NUG9UiSJHVqtV6t+Rpwf0TcDLzc3JiZp9WlKkmSpE6q1nB2Y/GSJElSHa01nEXEP2bmY5k5pa0KkiRJ6szWdc7Zdc0TEXFNnWuRJEnq9NYVzqJq+l31LESSJEnrDmfZyrQkSZLqYF0XBOweES9SGUHbopimmM/M3Lqu1UmSJHUyaw1nmdm1rQqRJElS7c/WlCRJUhswnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEDGeSJEklYjiTJEkqEcOZJElSiRjOJEmSSqRu4SwiekTErIi4NyLmRcS3i/ZBEfGXiHgoIq6IiM2L9u7F/MJi+cB61SZJklRW9Rw5ex3YNzN3B4YDB0XEKOD7wPmZORh4Hjil6H8K8Hxm7gycX/STJEnqVOoWzrJiaTG7WfFKYF/g6qJ9CjC+mD6imKdYvl9ERL3qkyRJKqO6nnMWEV0jYi7wDHAz8DDwQmYuK7o0Av2L6f7A4wDF8iVAnxa2OTEiZkfE7KampnqWL0mS1ObqGs4yc3lmDgcGAHsC722pW/He0ihZrtGQeVFmNmRmQ79+/TZdsZIkSSXQJldrZuYLwHRgFLBNRHQrFg0AniymG4GdAIrlvYDn2qI+SZKksqjn1Zr9ImKbYnoLYH9gPjANOLroNgG4vpieWsxTLL8tM9cYOZMkSXo767buLhtsB2BKRHSlEgKvzMwbIuJB4HcR8a/AHOCXRf9fApdFxEIqI2bH1bE2SZKkUqpbOMvM+4ARLbQ/QuX8s9XbXwOOqVc9Wj8DJ93Y3iVsEovOObS9S5Akab34hABJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJVK3cBYRO0XEtIiYHxHzIuILRXvviLg5Ih4q3rct2iMiLoiIhRFxX0SMrFdtkiRJZVXPkbNlwBcz873AKODUiNgVmATcmpmDgVuLeYCDgcHFayJwYR1rkyRJKqW6hbPMfCoz/1pMvwTMB/oDRwBTim5TgPHF9BHApVlxF7BNROxQr/okSZLKqE3OOYuIgcAI4C/A9pn5FFQCHLBd0a0/8HjVao1F2+rbmhgRsyNidlNTUz3LliRJanN1D2cRsRVwDXB6Zr64tq4ttOUaDZkXZWZDZjb069dvU5UpSZJUCnUNZxGxGZVgdnlm/r5ofrr5cGXx/kzR3gjsVLX6AODJetYnSZJUNvW8WjOAXwLzM/OHVYumAhOK6QnA9VXtJxZXbY4CljQf/pQkSeosutVx22OATwL3R8Tcou1s4Bzgyog4BXgMOKZY9kfgEGAh8Apwch1rkyRJKqW6hbPMvJ2WzyMD2K+F/gmcWq96JEmSOgKfECBJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpROoWziLiVxHxTEQ8UNXWOyJujoiHivdti/aIiAsiYmFE3BcRI+tVlyRJUpnVc+TsEuCg1domAbdm5mDg1mIe4GBgcPGaCFxYx7okSZJKq1u9NpyZMyJi4GrNRwBji+kpwHTgrKL90sxM4K6I2CYidsjMp+pVnzqHgZNubO8SNtqicw5t7xIkSW2orc852745cBXv2xXt/YHHq/o1Fm1riIiJETE7ImY3NTXVtVhJkqS2VpYLAqKFtmypY2ZelJkNmdnQr1+/OpclSZLUtto6nD0dETsAFO/PFO2NwE5V/QYAT7ZxbZIkSe2urcPZVGBCMT0BuL6q/cTiqs1RwBLPN5MkSZ1R3S4IiIjfUjn5v29ENALfAs4BroyIU4DHgGOK7n8EDgEWAq8AJ9erLkmSpDKr59Wax7eyaL8W+iZwar1qkSRJ6ijKckGAJEmSMJxJkiSViuFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSViOFMkiSpRAxnkiRJJWI4kyRJKhHDmSRJUokYziRJkkrEcCZJklQihjNJkqQSMZxJkiSVSLf2LkDS2g2cdGN7l7BJLDrn0PYuQZI6BEfOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI4UySJKlEvFpTUpt4O1x16hWnktqCI2eSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mSpBIxnEmSJJWI9zmTpBq9He7VBt6vTSo7w5kkdTJvh5BpwNTbmYc1JUmSSsSRM0mStMHeDiOxUK7R2FKFs4g4CPgx0BW4ODPPaeeSJEklZCDQ21lpwllEdAV+ChwANAJ3R8TUzHywfSuTJKk+3i4hU5tWmc452xNYmJmPZOYbwO+AI9q5JkmSpDZVmpEzoD/weNV8I/CB1TtFxERgYjG7NCL+Vue6+vLteLbOn/F21xdwH24c9+HGcf9tPPfhxnMfbpy67r/4fr22vIp31tKpTOEsWmjLNRoyLwIuqn85FRExOzMb2urz3o7chxvPfbhx3H8bz3248dyHG6cz7b8yHdZsBHaqmh8APNlOtUiSJLWLMoWzu4HBETEoIjYHjgOmtnNNkiRJbao0hzUzc1lEfB74Lyq30vhVZs5r57KgDQ+hvo25Dzee+3DjuP82nvtw47kPN06n2X+RucZpXZIkSWonZTqsKUmS1OkZziRJkkrEcFaIiIMi4m8RsTAiJrWwvHtEXFEs/0tEDGz7Ksuthn24T0T8NSKWRcTR7VFjmdWw/86MiAcj4r6IuDUiarpfTmdSwz7854i4PyLmRsTtEbFre9RZZuvah1X9jo6IjIhOcWuDWtXwGzwpIpqK3+DciPh0e9RZZrX8BiPiY8Xfh/Mi4jdtXWPdZWanf1G5AOFh4F3A5sC9wK6r9fkc8LNi+jjgivauu0yvGvfhQGAYcClwdHvXXKZXjftvHLBlMf1Zf4MbtA+3rpo+HLipvesu06uWfVj06wnMAO4CGtq77rK8avwNngT8pL1rLeurxn04GJgDbFvMb9fedW/qlyNnFbU8OuoIYEoxfTWwX0S0dOPczmqd+zAzF2XmfcCK9iiw5GrZf9My85Vi9i4q9wLUW2rZhy9Wzb6DFm503cnV+hi97wLnAq+1ZXEdgI8h3Hi17MPPAD/NzOcBMvOZNq6x7gxnFS09Oqp/a30ycxmwBOjTJtV1DLXsQ7VuffffKcB/1rWijqemfRgRp0bEw1TCxWltVFtHsc59GBEjgJ0y84a2LKyDqPW/448WpydcHRE7tbC8M6tlH74HeE9E3BERd0XEQW1WXRsxnFXU8uiomh4v1Ym5fzZOzfsvIj4BNAD/t64VdTy1PgLup5n5buAs4Ot1r6pjWes+jIguwPnAF9usoo6llt/gH4CBmTkMuIW3jsioopZ92I3Koc2xwPHAxRGxTZ3ralOGs4paHh21sk9EdAN6Ac+1SXUdg4/f2jg17b+I2B/4GnB4Zr7eRrV1FOv7G/wdML6uFXU869qHPYHdgOkRsQgYBUz1ooCV1vkbzMzFVf/t/gJ4fxvV1lHU+u/x9Zn5ZmY+CvyNSlh72zCcVdTy6KipwIRi+mjgtizORBTg47c21jr3X3E46edUgtnb7hyLTaCWfVj9F/ihwENtWF9HsNZ9mJlLMrNvZg7MzIFUzn08PDNnt0+5pVPLb3CHqtnDgfltWF9HUMu/JddRuUCKiOhL5TDnI21aZZ2V5vFN7SlbeXRURHwHmJ2ZU4FfApdFxEIqI2bHtV/F5VPLPoyIPYBrgW2BwyLi25n5vnYsuzRq/A3+X2Ar4KriWpTHMvPwdiu6ZGrch58vRh/fBJ7nrf/hEjXvQ7Wixv13WkQcDiyj8m/JSe1WcAnVuA//CzgwIh4ElgNfzszF7Vf1pufjmyRJkkrEw5qSJEklYjiTJEkqEcOZJElSiRjOJEmSSsRwJkmSVCKGM0mlExFL2/nzx0bEGo8nioiGiLhgPbc1OSK+tOmqk/R2533OJHV6EdE1M5evq19xs1VvuCqprhw5k1RKEbFVRNwaEX+NiPsj4oiqZScWD46+NyIuK9q2j4hri7Z7I2J00X5dRNwTEfMiYmLVNpZGxHci4i/AXhFxUEQsiIjbgaNaqWnliFoxIvariJgeEY9ExGlV/b4WEX+LiFuAXara3x0RNxX1/DkihhTt10fEicX0P0XE5ZtwV0rqYBw5k1RWrwFHZuaLxSNa7oqIqcCuVJ4vOiYzn42I3kX/C4A/ZeaREdGVytMUAD6Vmc9FxBbA3RFxTXE38XcAD2TmNyOiB5VHOe0LLASuqLHGIVQeI9MT+FtEXAgMo/IEkRFU/o79K3BP0f8i4J8z86GI+ADw/4rPnAjcERGPUnmo+Kj13VmS3j4MZ5LKKoB/i4h9gBVAf2B7KmHm6sx8FiAznyv67wucWLQtB5YU7adFxJHF9E5UHpC8mMpjX64p2ocAj2bmQwAR8R9UAtO63Fg8xPr1iHimqO+DwLWZ+UqxranF+1bAaN56/BZA96LepyPim8A0KoH0OSR1WoYzSWV1AtAPeH9mvhkRi4AeVEJbTc+di4ixwP7AXpn5SkRML7YB8Npq55ltyLPsXq+aXs5bf6e2tK0uwAuZObyVbQ2lEhp33IA6JL2NeM6ZpLLqBTxTBLNxwDuL9luBj0VEH4Cqw5q3Ap8t2rpGxNbFNp4vgtkQWj9cuAAYFBHvLuaP34i6ZwBHRsQWEdETOAwgM18EHo2IY4oaIyJ2L6b3BA6mcij0SxExaCM+X1IHZziTVCoR0Y3KiNTlQENEzKYyirYAIDPnAd8D/hQR9wI/LFb9AjAuIu6nco7X+4CbgG4RcR/wXeCulj4zM1+jchjzxuKCgL9vaP2Z+Vcq56zNpXLY9M9Vi08ATinqngccERHdgV9QOTfuSSrnnP0qqo59SupcInNDRvIlqT6K0aRfZOae7V2LJLUHR84klUZE/DPwW+Dr7V2LJLUXR84kSZJKxJEzSZKkEjGcSZIklYjhTJIkqUQMZ5IkSSViOJMkSSqR/w9+hxR+fI2xOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "\n",
    "plt.figure(figsize = (10,7))\n",
    "plt.title('Mean Jaccard index, N = {}'.format(N))\n",
    "plt.hist(jaccard_means, label='Within  a playlist')\n",
    "plt.hist(jaccard_rand_means, label='From random sample')\n",
    "plt.xlabel('Jaccard index')\n",
    "plt.ylabel('Frequency')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this histogram, we see that there's clearly two separate populations-- the Jaccard index is much larger for songs that are 'related' (they're together in playlists)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Content analysis: Can we predict a Jaccard index from audio features?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, we can pick a random song, and explore what playlists it is in. Instead of just looking at the Jaccard index between it and a single playlist it is in, it would be more interesting to explore songs from ALL the playlists it is in, and then compare it to random songs.\n",
    "\n",
    "Essentially, we will try to see if there's a linear connection between the audio features of a song, and the Jaccard index between it and a random other song.\n",
    "\n",
    "We are picking an equal distribution of random songs, and songs that populate the same playlists as the reference song, in order to ensure that there's a distribution of Jaccard indices (given that most will hae an index of zero).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, pick a random track\n",
    "\n",
    "track_EDA = pd.read_sql_query(\"SELECT * from tracks WHERE num_member > 3 ORDER BY RANDOM() LIMIT 1;\", conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parse all the playlists for a single track:\n",
    "\n",
    "track_EDA_playlists = track_EDA[['playlist_member']].values[0]\n",
    "track_EDA_playlists = [x.strip() for x in track_EDA_playlists[0].split(',')]\n",
    "track_EDA_playlists = [int(x) for x in track_EDA_playlists]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find all tracks that should have a higher Jaccard index from the playlists\n",
    "EDA_ids = []\n",
    "for playlist_id in track_EDA_playlists:\n",
    "    query = 'SELECT tracks FROM playlists WHERE playlist_id = {};'.format(playlist_id)\n",
    "    all_tracks = pd.read_sql_query(query, conn)\n",
    "    ids = all_tracks['tracks'].values[0].split(',')\n",
    "    EDA_ids.append(ids)\n",
    "\n",
    "# Flatten the list\n",
    "EDA_ids = [item for sublist in EDA_ids for item in sublist]\n",
    "\n",
    "# Make unique and drop track\n",
    "EDA_ids = list(set(EDA_ids))\n",
    "\n",
    "track1 = track_EDA[['track']].values[0][0]\n",
    "\n",
    "# Ensure that the tracks do NOT contain the track itself\n",
    "EDA_ids.remove(track1)\n",
    "\n",
    "\n",
    "# Make tuples\n",
    "\n",
    "track_pairs_EDA = [(track1, track2) for track2 in EDA_ids]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Between this song and random songs:\n",
    "\n",
    "len_EDA = len(EDA_ids)\n",
    "\n",
    "# Then pick the same number of songs\n",
    "\n",
    "tracks_rand_EDA = pd.read_sql_query(\"SELECT * from tracks ORDER BY RANDOM() LIMIT {};\".format(len_EDA), conn)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Place the IDs of the random song into a list, and ensure only unique values\n",
    "\n",
    "EDA_rand_ids = tracks_rand_EDA['track'].values.tolist()\n",
    "EDA_rand_ids = list(set(EDA_rand_ids))\n",
    "\n",
    "# Make sure the reference song itself is not in the playlist\n",
    "try:\n",
    "    EDA_rand_ids.remove(track1)\n",
    "except:\n",
    "    pass\n",
    "\n",
    "# Concatenate lists\n",
    "EDA_ids = EDA_ids+EDA_rand_ids\n",
    "\n",
    "# Create tuples\n",
    "track_pairs_EDA = [(track1, track2) for track2 in EDA_ids]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jaccard\n",
    "jaccard_EDA = []\n",
    "\n",
    "for pair_tuple in track_pairs_EDA:\n",
    "    j_EDA,_,_,_,_ = similarity_btw_pair(pair_tuple)\n",
    "    jaccard_EDA.append(j_EDA)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect audio features:\n",
    "import spotipy as sp\n",
    "import spotipy.oauth2 as oauth2\n",
    "import sqlite3\n",
    "import pandas as pd\n",
    "\n",
    "# set up authorization token\n",
    "credentials = oauth2.SpotifyClientCredentials(\n",
    "        client_id='153369a05314402294db1a574caaff2a',\n",
    "        client_secret='c6fff0923a0c44c5851fc4415038e8fa')\n",
    "\n",
    "token = credentials.get_access_token()\n",
    "spotify = sp.Spotify(auth=token)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n",
      "retrying ...1secs\n"
     ]
    }
   ],
   "source": [
    "# All audio features -- for regression purposes\n",
    "\n",
    "danceability = []\n",
    "energy = []\n",
    "key = []\n",
    "loudness = []\n",
    "mode = []\n",
    "speechiness = []\n",
    "acousticness = []\n",
    "instrumentalness = []\n",
    "liveness = []\n",
    "valence = []\n",
    "tempo = []\n",
    "duration_ms = []\n",
    "time_signature = []\n",
    "\n",
    "# Collect into lists:\n",
    "\n",
    "for track_id in EDA_ids:\n",
    "    results = spotify.audio_features(tracks=track_id)\n",
    "    \n",
    "    # Add into lists\n",
    "    danceability.append(results[0]['danceability'])\n",
    "    energy.append(results[0]['energy'])\n",
    "    key.append(results[0]['key'])\n",
    "    loudness.append(results[0]['loudness'])\n",
    "    mode.append(results[0]['mode'])\n",
    "    speechiness.append(results[0]['speechiness'])\n",
    "    acousticness.append(results[0]['acousticness'])\n",
    "    instrumentalness.append(results[0]['instrumentalness'])\n",
    "    liveness.append(results[0]['liveness'])\n",
    "    valence.append(results[0]['valence'])\n",
    "    tempo.append(results[0]['tempo'])\n",
    "    duration_ms.append(results[0]['duration_ms'])\n",
    "    time_signature.append(results[0]['time_signature'])\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dict, in order to make a pandas DataFrame\n",
    "EDA_data = {}\n",
    "\n",
    "EDA_data['ids']= EDA_ids\n",
    "EDA_data['jaccard_EDA'] = jaccard_EDA\n",
    "\n",
    "# Audio features\n",
    "EDA_data['danceability']= danceability\n",
    "EDA_data['energy']= energy\n",
    "EDA_data['key']= key\n",
    "EDA_data['loudness']= loudness\n",
    "EDA_data['mode']= mode\n",
    "EDA_data['speechiness']= speechiness\n",
    "EDA_data['acousticness']= acousticness\n",
    "EDA_data['instrumentalness']= instrumentalness\n",
    "EDA_data['liveness']= liveness\n",
    "EDA_data['valence']= valence\n",
    "EDA_data['tempo']= tempo\n",
    "EDA_data['duration_ms']= duration_ms\n",
    "EDA_data['time_signature']= time_signature\n",
    "\n",
    "EDA_results = pd.DataFrame.from_dict(EDA_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>jaccard_EDA</th>\n",
       "      <th>danceability</th>\n",
       "      <th>energy</th>\n",
       "      <th>key</th>\n",
       "      <th>loudness</th>\n",
       "      <th>mode</th>\n",
       "      <th>speechiness</th>\n",
       "      <th>acousticness</th>\n",
       "      <th>instrumentalness</th>\n",
       "      <th>liveness</th>\n",
       "      <th>valence</th>\n",
       "      <th>tempo</th>\n",
       "      <th>duration_ms</th>\n",
       "      <th>time_signature</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6pK6XXI3JHudVtl77QAfhh</td>\n",
       "      <td>0.024390</td>\n",
       "      <td>0.547</td>\n",
       "      <td>0.395</td>\n",
       "      <td>0</td>\n",
       "      <td>-9.654</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0396</td>\n",
       "      <td>0.7850</td>\n",
       "      <td>0.000016</td>\n",
       "      <td>0.1240</td>\n",
       "      <td>0.697</td>\n",
       "      <td>140.992</td>\n",
       "      <td>197293</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5BsOb3PMGwsMKTRpkabonh</td>\n",
       "      <td>0.073171</td>\n",
       "      <td>0.655</td>\n",
       "      <td>0.692</td>\n",
       "      <td>6</td>\n",
       "      <td>-11.728</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0303</td>\n",
       "      <td>0.1640</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.8210</td>\n",
       "      <td>0.668</td>\n",
       "      <td>111.253</td>\n",
       "      <td>183827</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1dFITmhIit06htQ6dBRQOS</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.830</td>\n",
       "      <td>0.649</td>\n",
       "      <td>8</td>\n",
       "      <td>-6.789</td>\n",
       "      <td>0</td>\n",
       "      <td>0.1510</td>\n",
       "      <td>0.0829</td>\n",
       "      <td>0.000220</td>\n",
       "      <td>0.0792</td>\n",
       "      <td>0.456</td>\n",
       "      <td>123.979</td>\n",
       "      <td>367960</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7Lg5osturAd2VUP1rf0tEU</td>\n",
       "      <td>0.024390</td>\n",
       "      <td>0.369</td>\n",
       "      <td>0.477</td>\n",
       "      <td>7</td>\n",
       "      <td>-8.836</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0354</td>\n",
       "      <td>0.7180</td>\n",
       "      <td>0.000037</td>\n",
       "      <td>0.1040</td>\n",
       "      <td>0.524</td>\n",
       "      <td>94.904</td>\n",
       "      <td>189213</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>27wLc0cxnw9IbhRMyMeBmR</td>\n",
       "      <td>0.333333</td>\n",
       "      <td>0.775</td>\n",
       "      <td>0.553</td>\n",
       "      <td>1</td>\n",
       "      <td>-10.512</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0308</td>\n",
       "      <td>0.4850</td>\n",
       "      <td>0.842000</td>\n",
       "      <td>0.0453</td>\n",
       "      <td>0.820</td>\n",
       "      <td>120.005</td>\n",
       "      <td>286093</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      ids  jaccard_EDA  danceability  energy  key  loudness  \\\n",
       "0  6pK6XXI3JHudVtl77QAfhh     0.024390         0.547   0.395    0    -9.654   \n",
       "1  5BsOb3PMGwsMKTRpkabonh     0.073171         0.655   0.692    6   -11.728   \n",
       "2  1dFITmhIit06htQ6dBRQOS     0.500000         0.830   0.649    8    -6.789   \n",
       "3  7Lg5osturAd2VUP1rf0tEU     0.024390         0.369   0.477    7    -8.836   \n",
       "4  27wLc0cxnw9IbhRMyMeBmR     0.333333         0.775   0.553    1   -10.512   \n",
       "\n",
       "   mode  speechiness  acousticness  instrumentalness  liveness  valence  \\\n",
       "0     0       0.0396        0.7850          0.000016    0.1240    0.697   \n",
       "1     1       0.0303        0.1640          0.000000    0.8210    0.668   \n",
       "2     0       0.1510        0.0829          0.000220    0.0792    0.456   \n",
       "3     1       0.0354        0.7180          0.000037    0.1040    0.524   \n",
       "4     1       0.0308        0.4850          0.842000    0.0453    0.820   \n",
       "\n",
       "     tempo  duration_ms  time_signature  \n",
       "0  140.992       197293               4  \n",
       "1  111.253       183827               4  \n",
       "2  123.979       367960               4  \n",
       "3   94.904       189213               4  \n",
       "4  120.005       286093               4  "
      ]
     },
     "execution_count": 384,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Resulting pd:\n",
    "EDA_results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split into train and test sets\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Response and independent variables:\n",
    "X_vars = EDA_results[['danceability', 'energy', 'key', 'loudness', 'mode','speechiness', \\\n",
    "                      'acousticness','instrumentalness','liveness', 'valence', 'tempo', \\\n",
    "                      'duration_ms', 'time_signature']]\n",
    "\n",
    "response = EDA_results[['jaccard_EDA']]\n",
    "\n",
    "# Split:\n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X_vars,response,test_size=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Random Forest Regression using OOB (lab 9)\n",
    "\n",
    "In order to do the regression, we can start with a random forest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "odict_values([[200, 400, 600, 800], [0.2, 0.4, 0.6, 0.8]])"
      ]
     },
     "execution_count": 331,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "# Code from lab 9, and:\n",
    "# Adventures in scikit-learn's Random Forest by Gregory Saunders\n",
    "from itertools import product\n",
    "from collections import OrderedDict\n",
    "\n",
    "param_dict = OrderedDict(\n",
    "    n_estimators = [200, 400, 600, 800],\n",
    "    max_features = [0.2, 0.4, 0.6, 0.8]\n",
    ")\n",
    "\n",
    "param_dict.values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(600, 0.2)"
      ]
     },
     "execution_count": 333,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use Random Forest, and pick the best combination:\n",
    "ytrain = np.ravel(ytrain)\n",
    "RF_results = {}\n",
    "RF_estimators = {}\n",
    "for Ntrees, maxfeat in product(*param_dict.values()):\n",
    "    params = (Ntrees, maxfeat)\n",
    "    RF = RandomForestRegressor(oob_score=True, \n",
    "                                n_estimators=Ntrees, max_features=maxfeat, max_depth=50)\n",
    "    RF.fit(Xtrain, ytrain)\n",
    "    RF_results[params] = RF.oob_score_\n",
    "    RF_estimators[params] = RF\n",
    "best_params = max(RF_results, key = RF_results.get)\n",
    "best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The resulting accuracy from the RF is 7.116939558688562 %\n"
     ]
    }
   ],
   "source": [
    "# Check accuracy\n",
    "RF_final = RF_estimators[best_params]\n",
    "print('The resulting accuracy from the RF is', RF_final.score(Xtest, ytest)*100, '%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'pd' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-39fa6f540894>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Plot based on significance\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSeries\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mRF_final\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfeature_importances_\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mXtrain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msort_values\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkind\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"barh\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'pd' is not defined"
     ]
    }
   ],
   "source": [
    "# Plot based on significance\n",
    "pd.Series(RF_final.feature_importances_,index=list(Xtrain)).sort_values().plot(kind=\"barh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is very low accuracy, so regression does not seem too appealing. In any case, we can try a different method before rejecting the regression based on linearity.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boosting\n",
    "\n",
    "Last effort before rejecting this method will be boosting, using the \"best predictors\" from the previous analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the worst predictors\n",
    "ind_todrop = [index for index, value in enumerate(RF_final.feature_importances_) \\\n",
    "              if value < RF_final.score(Xtest, ytest)*0.5]\n",
    "\n",
    "Xtrain_new = Xtrain.copy()\n",
    "\n",
    "for ind in ind_todrop:\n",
    "    col_todrop = Xtrain.columns[ind]\n",
    "    Xtrain_new.drop(col_todrop, axis = 1)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to use the best parameters from the random forest eda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 372,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fit an Adaboost Model\n",
    "\n",
    "from sklearn.ensemble import AdaBoostRegressor\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "boosts = {}\n",
    "boost_results = {}\n",
    "\n",
    "for n in range(2,20):\n",
    "    boost_model = AdaBoostRegressor(base_estimator=DecisionTreeRegressor(max_depth=n),\\\n",
    "                                    n_estimators=best_params[0],\\\n",
    "                                    learning_rate=0.05)\n",
    "    boost_model.fit(Xtrain_new, ytrain)\n",
    "    boosts[str(n)] = boost_model\n",
    "    \n",
    "    boost_results[str(n)] = boost_model.score(Xtest, ytest)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best boosting model has a max_depth of 11 with a score of 8.28111106007271%\n"
     ]
    }
   ],
   "source": [
    "#Best depth accuracy\n",
    "\n",
    "best_depth = max(boost_results, key = boost_results.get)\n",
    "boost_final_score = boost_results[best_depth]\n",
    "\n",
    "print('The best boosting model has a max_depth of {0} with a score of {1}%'.format(best_depth,100*boost_final_score))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is stillnot a satisfactory accuracy score, which shows the need for both a different content-based analysis, and a collaborative filtering approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note! always close the connection at the end\n",
    "cur.close()\n",
    "conn.close()"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}