"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(2,5, subplot_kw=dict(xticks=[], yticks=[]))\n",
"for ax, digit in zip(axes.flat, model.cluster_centers_):\n",
" ax.imshow(digit.reshape(8,8), cmap=\"gray\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One can recognize these numbers with the exception of maybe number eight. What is the accuracy score of this clustering?"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-24T19:29:57.808213Z",
"start_time": "2020-06-24T19:29:57.802233Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[4, 3, 5, 9, 7, 0, 1, 8, 2, 6]\n",
"Accuracy score is 0.793544796884\n"
]
}
],
"source": [
"permutation3 = find_permutation(10, digits.target, model.labels_)\n",
"print(permutation3)\n",
"acc = accuracy_score(digits.target, [ permutation3[label] for label in model.labels_])\n",
"print(\"Accuracy score is\", acc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is quite a good result for such a simple algorithm!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercise 5 (plant clustering)
\n",
"\n",
"Using the same iris data set that you saw earlier in the classification, apply k-means clustering with 3 clusters.\n",
"Create a function `plant_clustering` that loads the iris data set, clusters the data and returns the accuracy_score.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercise 6 (nonconvex clusters)
\n",
"\n",
"This exercise can give four points at maximum!\n",
"\n",
"Read the tab separated file data.tsv from the `src` folder into a DataFrame. The dataset has two features X1 and X2, and the label y. Cluster the feature matrix using DBSCAN with different values for the eps parameter. Use values in `np.arange(0.05, 0.2, 0.05)` for clustering. For each clustering, collect the accuracy score, the number of clusters, and the number of outliers. Return these values in a DataFrame, where columns and column names are as in the below example.\n",
"\n",
"Note that DBSCAN uses label -1 to denote outliers , that is, those data points that didn't fit well in any cluster. You have to modify the find_permutation function to handle this: ignore the outlier data points from the accuracy score computation. In addition, if the number of clusters is not the same as the number of labels in the original data, set the accuracy score to NaN.\n",
"\n",
" eps Score Clusters Outliers \n",
" 0 0.05 ? ? ?\n",
" 1 0.10 ? ? ?\n",
" 2 0.15 ? ? ?\n",
" 3 0.20 ? ? ?\n",
"\n",
"Before submitting the solution, you can plot the data set (with clusters colored) to see what kind of data we are dealing with.\n",
"\n",
"Points are given for each correct column in the result DataFrame.\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Hierarchical clustering\n",
"\n",
"Hierarchical clustering works by first putting each data point in their own cluster and then merging clusters based on some rule, until there are only the wanted number of clusters remaining. For this to work, there needs to be a distance measure between the data points. With this distance measure `d`, we can define another distance measure between the **clusters** U and V using one of the following methods (*linkages*):\n",
"\n",
"* `single`: $d(U, V) := \\min_{u \\in U, v \\in V} d(u,v)$\n",
"* `complete`: $d(U, V) := \\max_{u \\in U, v \\in V} d(u,v)$\n",
"* `average`: $d(U, V) := \\sum_{u \\in U, v \\in V} \\frac{d(u,v)}{|U||V|}$\n",
"* `ward`: tries to minimize the variance in each cluster\n",
"\n",
"At each iteration of the algorithm two clusters that are closest to each other are merged. After this the distance between the clusters are recomputed, and then it continues to the next iteration.\n",
"\n",
"Below is an example with a botanical dataset with 150 samples from three species. Each species appears in the dataset 50 times. Each sample point has 4 features, which are basically dimensions of the \"leaves\" of the flower.\n",
"\n",
"We use the [seaborn](https://seaborn.pydata.org/index.html) library to both to compute the clustering and to visualize the result. The visualization consists of two parts: the *heatmap*, whose rows and/or columns may be reordered so as to have the elements of the same cluster next to each other; and the *dendrogram*, which shows the way the clusters were merged. The colors give the length of the corresponding features."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-24T19:29:58.300009Z",
"start_time": "2020-06-24T19:29:57.810710Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['setosa' 'versicolor' 'virginica']\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAI/CAYAAAC1XpeNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3X9c1GW+///nMMAo4O8Ef8TuBmp67KTtRyvX0o1CLTXMstLsGFvH2kRSydQ8ZdmSritm22wlax2tVm+brkllpYm/2rSsdsvaslbTE2wILqiI4DAzzPcPv4wgCPhm4M28edxvN27LvGe4rhdrwGuu1/u6Xjafz+cTAAAAGi3E7AAAAACsgsQKAAAgQEisAAAAAoTECgAAIEBIrAAAAAIktK4n3f/+vkkmDbsorknGBQAAMFOdiZW87mYKAwAAIPjVmVj53OXNFUc133//vWbOnOl/nJOTo9TUVN1zzz3+ax9//LEefPBBXXzxxZKkxMREpaSkNHeoAAAAfnWvWHlczRRGdXFxccrKypIkeb1eDRs2TImJiTVeN2jQIK1YsaK5wwMAAKhV3StWXk9zxXFee/bsUWxsrHr27Gl2KAAAAHWqe8XKbc6KVVWbNm3SmDFjan3u888/180336zo6GjNmTNHvXv3buboAAAAzrLV1SvQ9dX7TTKp47KaZb3alJeX69prr9WmTZt00UUXVXuupKRENptNkZGR2rlzp9LT07Vly5amCBcAAKBB6jzHyudxNclHQ+3atUv9+/evkVRJUlRUlCIjIyVJw4cPl8fjUVFR0QV++wAAAIHTokuBmzZt0ujRo2t97ujRo7roootks9m0b98+VVRUqFOnTs0cIQAAwFn17Ao057gFSSotLdXu3bu1cOFC/7W1a9dKkiZOnKjNmzdr7dq1stvtatOmjZYtWyabzWZWuAAAAHXfY1W2a1WTTNp22D1NMi4AAICZ6ikFmrdiBQAAEGxabCkQAAAg2LTKFavy8nKdOHFCXbt2NTsUAABgIXUetyB3edN8mGDmzJk6efKkTp8+rbFjx2r06NF66aWXTIkFAABYU6tJrA4dOqR27dppx44duuqqq7Rz505t3LjRlFgAAIA11VMKdDdTGE3P4znT9/CTTz7R8OHD1bZtW4WE1J1XAgAAXIh6bl63TmIVHx+v++67T99//73S0tJ0+vRps0MCAAAW02JvXi8uLtb//M//6LvvvpPNZtPTTz+tK664wv+8z+dTenq6du7cqTZt2mjx4sXq37//ecf77W9/q7/+9a+69NJLFRERofz8fKWlpTXHtwIAAFqJOhMrn4mlwPT0dF177bX6/e9/r/Ly8horTLt27dLhw4e1ZcsWffHFF3riiSe0bt26Wsfyer2aOHGi3njjDf+1mJgYxcTENOn3AAAAWpe6V6zKzUmsTp48qU8++USLFy+WJIWHhys8PLzaa7KzszVu3DjZbDYNHDhQxcXFKigoUHR0dI3x7Ha7IiIi5HK55HA4muV7AAAArU89pUBPM4VRXW5urjp37qx58+Zp//796t+/v+bPn6+IiAj/a/Lz89WtWzf/427duik/P7/WxEqSLrnkEt11110aOXJktXHuuuuupvtGAABAq1Lntjif290kH/XxeDz6+uuvNXHiRG3cuFFt27ZVZmZmo75Rr9er3r176/vvv9dXX33l/wAAAAiUFrli1a1bN3Xr1k0DBgyQJI0aNapGYhUTE6MjR474Hx85cqTOe6YWLVrUNMECAAD8/+o5INTdNB/16Nq1q7p166bvv/9ekrRnzx7Fx8dXe01CQoI2btwon8+nzz//XO3atTtvGVCSysrKtHz5cv9OwIMHD2rr1q31xgIAANBQde8KLDdnxUqSHnvsMT388MNyu92KjY3VokWLtHbtWknSxIkTNXz4cO3cuVOJiYlq27atnn766TrHe+KJJ9S1a1ft379f0plVsbS0NN1www1N/r0AAIDWoUWWAiWpX79+2rBhQ7VrEydO9H9us9m0YMGCBo/37bff+s+ykqTIyEhVVFQEJlgAAAC14BWrQDv3uAaXyyWfz2dSNAAAwIrqOSDU21xxNLlBgwbpxRdfVHl5uT7++GP97//+rxISEswOCwAAWIjNV8eyzYnkprn/qMP/Nv9N4263WytXrtS2bdvk8/mUkJCg+++/X3a7vdljAQAA1lRnYnX8rqZZ0en4p21NMm5dDh48WGNnYW3XAAAAjKr7uAWPr2k+TPDwww836BoAAIBRdd5jVVEe/LvmioqKVFRUJJfLpYMHD/pvWD958qRKS0tNjg4AAFhJPbsCmyuMpvPWW29p9erVKigo0H//93/7r7dr10733XefiZEBAACrqfMeq6OJw5tk0q7v72yScevy4osv6oEHHmj2eQEAQOtR5z1WFW5bk3w0lNfr1bhx43T//ffXeG7Dhg26+uqrlZSUpKSkJK1bt67OsR544AHt2bNHr732miSpsLBQhw4danAsAAAA9amzFOgtb3gS1BReeeUVxcfHq6SkpNbnb7rpJj3++OMNGiszM1M7d+7U0aNHNXnyZLndbj366KP+NjkAAACNVeeKldcd0iQfDXHkyBHt2LFDt912W0C+0bffflurVq1SRESEpDO9As+XsAEAABhR94pVA5OgpvD0009r9uzZOnXq1Hlfs2XLFn3yySe65JJLNG/ePHXv3v28r23Tpo3CwsKqXbPZzF2RAwAA1lL3ipUnpEk+6rN9+3Z17txZl1122Xlfc91112nbtm1666239Itf/EJz5sypc8xu3brp008/lc1mU0VFhZ5//nn17t273lgAAAAaqs5dgd/2vbFJJr10/7t1Pp+RkaGsrCyFhobK5XKppKREiYmJWrp0aa2v93q9uvLKK/XZZ5+dd8yjR49qzpw52rt3r2w2mwYNGqSlS5eqS5cujfpeAAAAKtWZWH0dP7pJJv2Pg5sa/NqPP/5YL7/8slasWFHtekFBgaKjoyVJ77//vv74xz/q9ddfr3e8srIyVVRUKDIy8sKCBgAAqEed91h5vObdY1WbZ599Vpdddpmuv/56vfrqq9q2bZvsdrs6dOigRYsW1fv1P/zwg3744Qd5vV7/teHDm+asLgAA0PrUuWL1t9ikJpn05zlZTTJuXZYsWaKNGzfqkksuUUjImYTRZrPplVdeafZYAACANdW5YuX22psrjia3detWZWdnq23btmaHAgAALKruUmBFyyoFNkb37t1rHLcAAAAQSHWvWPmsc87T3Llz9cADD2jo0KEKDw/3X7/rrrtMjAoAAFhJ3YlV3cdcBZXMzEwdPXpU33zzjex265Q4AQBAy1FPYmWdFat//OMf2rx5M6etAwCAJlN3YmWzzorVz372M5WWlnJ+FQAAaDL1JFbWWd2JiorS+PHjde2111a7x+qRRx4xMSoAAGAl9SRWzRVGdS6XS3fddZfKy8vl9Xo1cuRIpaamVntNeXm5HnnkEf3jH/9Qx44d9cwzz+jiiy8+75hxcXGKi4tr6tABAEAr1iJXrMLDw7V69WpFRkbK7XZr0qRJGjZsmAYOHOh/zbp169S+fXu9//772rRpk5YuXarly5efd8yUlJTmCB0AALRiLXLFymaz+e+F8ng88ng8NW4637Ztmz9ZGjlypBYuXCifz1fjde+++65uvPFG/elPf6p1Lo5bAAAAgdIiEytJ8nq9Gj9+vH744QdNmjRJAwYMqPZ8fn6+unfvLkkKDQ1Vu3btdOzYMXXu3Lna6/75z3/qxhtv1FdffdVssQMAgNapnuMWzGO325WVlaXi4mJNmzZN3333nfr06XPB41TemzV//nxFRUVVe66kpCQgsQIAAEiq+wTQclvTfFyI9u3b66qrrtIHH3xQ7XpMTIzy8vIknSkXnjx5Up06dTrvOHfffXeDrgEAABhVZ2LltjXNR32KiopUXFwsSTp9+rR2795dY0dfQkKC3njjDUnS5s2bdfXVV9d6+KfH41FZWZkqKip0+vRplZWVqaysTAUFBSorK2vo/08AAAD1qrMU6JWvueKopqCgQHPnzpXX65XP59OoUaN03XXX6dlnn9Vll12m66+/Xrfddptmz56txMREdejQQc8880ytY7344otyOp2y2WzVdhVGRUUpOTm5ub4lAADQCth8Pt95s6cnfto0O+ae+L/ad+g1pYULF+rxxx9v9nkBAEDrUc+uQHNWrJpCZVJVeehopbZt25oVEgAAsJh6dgVaJ7F6//339dRTT6mgoEA2m81/5tU333xjdmgAAMAiWk1itWTJEi1fvlwDBw5USIh1mksDAICWo9UkVh06dNDPf/5zs8MAAAAWVk9iVdFccTS5xMRErVmzRjfddJMcDof/OvdYAQCAQKlzV+CvfnZbk0z68uH1TTJuXfr27ev/nHusAABAU6h7xcpnnRWr/fv3mx0CAACwuDoTK4+FSoGSdOjQIR08eFA33HCDTp06JbfbrY4dO5odFgAAsIi6W9qookk+6jNv3jwNGTJEY8aMqfX5jz/+WP/v//0/JSUlKSkpSU6ns94xN2zYoF//+tdatGiRJCk/P18zZsyo9+sAAAAaqkWWAsePH6/Jkydrzpw5533NoEGDtGLFigaP+corr+gvf/mL7rrrzGnycXFx+ve//93oWAEAACq1yF2BgwcPVm5ubkDHDAsLU2RkZLVrdrs9oHMAAIDWre57rFrwzeuff/65br75ZkVHR2vOnDnq3bt3na/v2LGjDh06JJvNJknKyspSt27dmiNUAADQStRTCvTW9bRp+vfvr23btikyMlI7d+7UtGnTtGXLljq/5tFHH1VaWpoOHTqkhIQEtWnTRi+++GIzRQwAAFqDoFyxioqK8n8+fPhwPfnkkyoqKlLnzp3P+zV5eXlauXKljh07Jp/Pp86dO+vbb7/VT37yk+YIGQAAtAJ17wr0eZvko7GOHj2qynNN9+3bp4qKCnXq1KnOr1myZIk6deqk+Ph49erVSx07dtSSJUsaHQsAAEClelaszCkFzpo1S3v37tWxY8c0bNgwTZ8+XR6PR5I0ceJEbd68WWvXrpXdblebNm20bNky/71T51N50nqlkJAQeb0ts9QJAACCU4u8x2rZsmV1Pj958mRNnjz5gsaMjIzUF198oQEDBkiSvvjiC0VERBiOEQAA4FwtcsWqKcyePVvTpk1Tr169JEkHDhxo0MGiAAAADVVnE+Y+XQc1yaTfHf20Scatz4kTJ/T5559LkgYOHKgOHTqYEgcAALCmOhOruIuuaJJJv//335tkXAAAADPVXQqs8DRXHAAAAEGv7pvXK6xzjxUAAEBTazU3rwMAADS1uhMrznkCAABoMEqBAAAAAVJnYuVtob0CAQAAWqK6EytWrAAAABqsznOsAAAA0HAhZgcAAABgFSRWAAAAAUJiBQAAECAkVgAAAAFCYgUAABAgJFYAAAABQmIFAAAQICRWAAAAAUJiBQAAECAkVgAAAAFCYgUAABAgJFYAAAABQmIFAAAQICRWAAAAAUJiBQAAECAkVgAAAAFCYgUAABAgJFYAAAABQmIFAAAQICRWAAAAAUJiBQAAECAkVgAAAAESanYAAAAAzcHlcumuu+5SeXm5vF6vRo4cqdTUVM2dO1d79+5Vu3btJEmLFy9Wv379DM1BYgUAAFqF8PBwrV69WpGRkXK73Zo0aZKGDRsmSXrkkUc0atSoRs9BYtVKFI4dbnYIMOCtv8WaHQIMet6XY3YIMOiLou/NDgEGuU7X/XNns9kUGRkpSfJ4PPJ4PLLZbAGNgXusLMrpdCojI8PsMAAAaFG8Xq+SkpL0i1/8Qr/4xS80YMAASdIzzzyjsWPH6umnn1Z5ebnh8VmxsiiXy6W0tDT/Y0+JicHAsFLe+gStCF+Y2SHAIG9FhdkhoAnZ7XZlZWWpuLhY06ZN03fffadZs2apa9eucrvdeuyxx5SZmamUlBRD4/NrGwAAtDrt27fXVVddpQ8++EDR0dGy2WwKDw/X+PHj9eWXXxoel8TKwpxOp9khAADQYhQVFam4uFiSdPr0ae3evVtxcXEqKCiQJPl8Pm3dulW9e/c2PAelQAtzuVz+z+3hJgYCAEHEHsKag1UVFBRo7ty58nq98vl8GjVqlK677jr913/9l44dOyafz6e+ffvqySefNDwHiZUFOZ1O5eSwIwkAgKr69u2rjRs31rj+yiuvBGwO0nILcrlcio2N1UcffWR2KAAAtCqsWFlYRZWdLV7jO0cBoFUJDbGbHQKCGCtWFtahQwezQwAAoFVhxcrCOnfubHYIAAC0GHl5eXrkkUdUWFgom82m22+/XVOmTNFzzz2n119/3f93c9asWRo+3FjHEhIri3I4HOrevbv/cWiUicEArVBECFtxg5Wnwmt2CGgidrtdc+fOVf/+/VVSUqJbb71VQ4cOlSTdc889uvfeexs9B6VAi0pJSVFsLH3mAACoFB0drf79+0uSoqKiFBcXp/z8/IDOQWJlYXv37jU7BAAAWqTc3Fx98803/l6Bf/rTnzR27FjNmzdPJ06cMDwupUAL83g8/s997AoMSmE+syOAUeWinBSs6BVofadOnVJqaqoeffRRRUVFaeLEiXrwwQdls9n07LPPavHixVq0aJGhsVmxsjBKgQAAVOd2u5WamqqxY8dqxIgRkqSLLrpIdrtdISEhmjBhAr0CUTsSKwAAzvL5fJo/f77i4uKUnJzsv17ZK1ASvQJxfg6Hw+wQACDo0CvQuj777DNlZWWpT58+SkpKknTmaIW3335b+/fvlyT17NlTCxcuNDwHiZWFOJ1OuVwuf5/AlJQUkyMCAKDlGDRokL799tsa142eWVUb0nILcblcSktLowQIAIBJWLFqJcqK+KcORp28bAsMVh47uwKDFb0C0RisWFmY0+k0OwQAAFqMvLw83X333brppps0evRorV69OuBzsIxhYS6Xy/+5PYxzWYJRSYjN7BBgkNvHz1ywoqWNdZ2vpU2vXr0CNgcrVgAAoFWgpQ0Mczqd2rJli9lhAADQIp3b0iZQKAValMvlUlRUlP+x100OHYzcVAKDVltbmNkhADiPc1vaBBJ/bQEAQKtRW0ubQCKxsrCSkhKzQwAAoMU4X0ubQKIUaGGBXt5E8wvjGKug5RE7y4CW5nwtbQJ58jqJlYVxAjsAAGedr6VNIFEKtCCHw6GcnJyAnssBAADqx4qVBaWkpCgjI6NaE+bwKI+JEcGoqApqgcGqR3g7s0MAYAJWrAAAAAKEFasg4nQ6q7WpOVdOTk6N11ddtQIAoDWbN2+eduzYoS5duujtt9+WJL377rtyOp06ePCg1q1bp//8z/9s1BwkVkHE5XIpLS3tvM9nZGTUeD0Acxz38fMHtDTjx4/X5MmTNWfOHP+1Pn366LnnntOCBQsCMgeJFQAAaBUGDx6s3Nzcatfi4+MDOgf3WFmUw+GgVyAAAM2MFSuLSklJ0fbt2/2PvS5y6GBUEkKzQAAIJvy1tTAOCAUAoHmRWFlYXl6e2SEAANCqUAq0sO7du/s/tzsqTIwERrmpBAJAwMyaNUt79+7VsWPHNGzYME2fPl0dO3bUU089paKiIt1///3q16+fXnrpJcNzkFhZGKVAAADOWrZsWa3XExMTAzYHpUALczgcZocAAECrwoqVhVU9db3sWLiJkcCor0Lp8RisTrhPmx0CDAoNsZsdAoIYiRUAAGgV8vLy9Mgjj6iwsFA2m0233367pkyZot/+9rfavn27wsLC9JOf/ESLFi1S+/btDc1BYmUx5/YTrGyBYw/j5vVgNLLM7AhgVLZKzQ4BBnkqvGaHgCZit9s1d+5c9e/fXyUlJbr11ls1dOhQDR06VGlpaQoNDdXvfvc7rVixQrNnzzY0B4mVxdTXTxAAgNYqOjpa0dHRkqSoqCjFxcUpPz9f11xzjf81AwcO1HvvvWd4DhIrC3E4HDpw4EC1VSuSLAAAasrNzdU333yjAQMGVLv+l7/8RTfeeKPhcUmsLCQlJUUZGRm1rlqdOsEOwWB0OJwf0WAV5uYG6GDlreDWCas7deqUUlNT9eijjyoqKsp//YUXXpDdbtfNN99seGx+awMAgFbD7XYrNTVVY8eO1YgRI/zXN2zYoB07dmjVqlWy2Yyfzsw5VhbjcDiUk5Mj6cyN7AAA4Ayfz6f58+crLi5OycnJ/uu7du3SypUr9cILL6ht27aNmoMVK4upLAdKqrY7MCycXS7BKOKU2RHAqAg75fdgZQ9hzcGqPvvsM2VlZalPnz5KSkqSdKbNzW9+8xuVl5f7k60BAwZo4cKFhuYgsQIAAK3CoEGD9O2339a4Pnz48IDNQVpuQQ6HQxkZGdqyZYvZoQAA0KqwYmVBla1sdu/ebXIkaKxjbCwLWmEVvG8FWiN+8i0sNjbW7BAAAGhVWLGysLy8PLNDAACgxXC5XLrrrrtUXl4ur9erkSNHKjU1VXv27NGSJUtUUVGhiIgILV68WD/96U8NzUFiZWHdu3c3OwQ0ksNndgQwyu3jkEmgpQkPD9fq1asVGRkpt9utSZMmadiwYXriiSf0/PPPKz4+Xn/605/0wgsvaPHixYbmoBRoYZQCAQA4y2azKTIyUpLk8Xjk8Xj8h4GWlJT4/7eyn6ARrFhZmMPBOToAAFTl9Xo1fvx4/fDDD5o0aZIGDBig9PR0TZ06VQ6HQ1FRUXr99dcNj09iZWGVuwMlqeMlp02MBEb9bK/H7BBgUFgYBQGgJbLb7crKylJxcbGmTZum7777TqtWrVJmZqYGDBiglStXatGiRUpPTzc0Pj/5AACg1Wnfvr2uuuoq7dq1S/v379eAAQMkSTfddJP+/ve/Gx6XxMrC6BUIAMBZRUVFKi4uliSdPn1au3fvVnx8vE6ePKlDhw5Jkj788EPFx8cbnoNSoIUdOHDA/7mv3MRAACCIeCvY0WlVBQUFmjt3rrxer3w+n0aNGqXrrrtOv/nNb5SamiqbzaYOHTro6aefNjwHiZWF5eTkmB0CAAAtRt++fbVx48Ya1xMTE5WYmBiQOSgFWlhoKHkzAADNib+8FnbllVf6P7eFmxgIDDthp1kg0NzsIaw5wDgSKwAA0Cp8//33mjlzpv9xTk6OUlNTdc899wRsDkOJldPplMvlClgQaJgLvWdq7969TRQJAADBJy4uTllZWZLOHBQ6bNiwgN1bVclQYuVyuZSWlhbQQFC/jIyMBr2uMvEtLCxs4ojQ1MLoFRi0Qm2UcYGWbM+ePYqNjVXPnj0DOi6FZAuqTHwr+x8BAIDqNm3apDFjxgR8XO6xsrATJ074Pw9pQ5IFAIAklZeXa9u2bU1SfWPFysJC2NkCAEANu3btUv/+/XXRRRcFfGz+8lpYoOvGAABYwaZNmzR69OgmGZtSYCthCyeHDkZhPu5eB4BAKi0t1e7du7Vw4cImGZ/EysJiY2PNDgEAgBYlIiJCH3/8cZONzzKGheXl5ZkdAgAArQorVhbm8XjMDgGN5ObIjKDl8XnNDgGACVixsrDi4mKzQwAAoFUhsbKwkydPmh0CAAAtjtfr1bhx43T//fdLOtMybsKECUpMTNSMGTNUXl5ueGxKgRZW9bgFbzFliWDErsDg5fZVmB0CDAoNoR2R1b3yyiuKj49XSUmJJGnp0qW65557NHr0aD3++ONav369Jk2aZGhsVqwsjF2BAABUd+TIEe3YsUO33XabJMnn8+mjjz7SyJEjJUm33HKLsrOzDY9PYmVh7AoEAKC6p59+WrNnz/Z3Jzl27Jjat2+v0NAzRbxu3bopPz/f8PiUAi3C6XTK5XJJOlMrlqrvCqRXINC8wmy8bw1WngpunbCq7du3q3Pnzrrsssua7CwrEiuLcLlc/maSGRkZkigFAgBQ1d/+9jdt27ZNu3btksvlUklJidLT01VcXCyPx6PQ0FAdOXJEMTExhufgLZWFkVgBAHBWWlqadu3apW3btmnZsmW6+uqrlZGRoauuukqbN2+WJL3xxhtKSEgwPAcrVhbmcDjMDgEAgg67Aluf2bNna+bMmVq+fLn69eunCRMmGB6LxMrCUlJSzA4BAIAW6aqrrtJVV10l6UyFZ/369QEZl1IgAABAgLBi1UpUnOagSaA5RYSEmx0CDGofHmF2CAhiJFYAAKBVyMvL0yOPPKLCwkLZbDbdfvvtmjJlipYvX67s7GyFhISoS5cuWrRokeGdgZQCLczpdJodAgAALYbdbtfcuXP1zjvv6M9//rPWrFmjAwcO6L777tNbb72lrKws/fKXv9Qf/vAHw3OwYmVhlQeGShwQGqxKQ/h3C1ZhvG8NWh4fB4RaVXR0tKKjoyVJUVFRiouLU35+vnr16uV/TVlZmWw24797+cm3oL179yojI0NbtmwxOxQAAFqk3NxcffPNNxowYIAk6ZlnntHw4cP11ltv6aGHHjI8LitWFuTxeJSWlqbdu3f7r3HzenAK458taLlVYXYIMKjc66n/RQhqp06dUmpqqh599FFFRUVJkmbOnKmZM2dqxYoVeu2115SammpobFasLIyT1wEAqM7tdis1NVVjx47ViBEjajw/duzYRlV8SKwsrGrNGACA1s7n82n+/PmKi4tTcnKy//rhw4f9n2dnZysuLs7wHJQCLahyparqyevcvB6con3lZocAgzraaCkFtDSfffaZsrKy1KdPHyUlJUmSZs2apfXr1+vQoUOy2Wzq2bOnnnzyScNzkFhZECVAAABqGjRokL799tsa14cPHx6wOSgFWpTT6VRGRobZYQAA0KqwYmVRLpdLaWlpZocBtFrHfa76X4QW6bSHEjyMY8UKAAAgQFixCiIOh+O85b2cnJxqrztw4ICcTme1G9gBAGjN5s2bpx07dqhLly56++23JUn79+/XggULVFpaqp49e2rp0qX+s62MILEKInUlSVUTrpSUFGVkZFRraQOgeZXTFiVohYbYzQ4BTWT8+PGaPHmy5syZ4782f/58zZkzR1deeaXWr1+vlStXasaMGYbnoBRoUQ6Ho9oqFgAArd3gwYPVoUOHatcOHz6swYMHS5KGDh3a6HZwJFYWlZKSory8PLPDAACgRevdu7eys7MlSe+9914TPyeQAAAgAElEQVSj/3ZSCrSw7t27mx0CGsnt470P0NwoBbYu6enpSk9P1/PPP6+EhASFh4c3ajwSKwvjoFAAAOoWHx+vl19+WZJ06NAh7dixo1Hj8XbYwhwOWmoAAFCXwsJCSVJFRYVeeOEF3XnnnY0ajxUrC6u6i9AWTg4djLq0KTM7BBhVYXYAMOpU+WmzQ0ATmTVrlvbu3atjx45p2LBhmj59ukpLS7VmzRpJUmJiom699dZGzUFiBQAAWoVly5bVen3KlCkBm4NlDAtzOp1mhwAAQKvCipVFOZ1OHThwwP/Y3iXCxGhgVKidelKw8ng5IDRY2UNYc4Bx/NdjUS6Xi12BAABUMW/ePA0ZMkRjxoyp8dzLL7+sSy+9VEVFRY2ag8TKwtgVCADAWePHj9fKlStrXM/Ly9OHH36oHj16NHoOSoEW43Q65XK5lJOTo+XLl/uvh3RuZ2JUMKrMxe6kYFXqO2l2CADOMXjwYOXm5ta4vmjRIs2ePVsPPvhgo+dgxcpiXC6X0tLSKAMCANAAW7duVXR0tPr27RuQ8Vixsqi9e/cqIyNDaWlpkiSfq9zkiGBE506lZocAg9qfaGN2CDCoTWjjWpogeJSVlWnFihX+k9cDgRUri+revbs/qQIAADX98MMPys3NVVJSkhISEnTkyBGNHz9eR48eNTwmiZVF9erVSxkZGWaHAQBAi3XppZdqz5492rZtm7Zt26Zu3bppw4YN6tq1q+ExKQVaVNV2NpJUUVRiUiRojPzCKLNDgEGRYS6zQwBwjtpa2kyYMCGgc5BYAQCAVuF8LW0qbdu2rdFzUAq0CIfDoYyMDOXk5Piv0dIGAIDmxYqVRVSW/qreV+VynS1FhHSmpBSMPD7e+wDNzVNBOyIYx29ti3I6ndVWrwAAQNNjxcqiXC6XevXqZXYYAAC0GHl5eXrkkUdUWFgom82m22+/XVOmTNHx48c1c+ZM/etf/1LPnj21fPlydejQwdAcJFYWVnVnILsCg1OB7SKzQ4BhZWYHAIMoBVqX3W7X3Llz1b9/f5WUlOjWW2/V0KFDtWHDBg0ZMkRTp05VZmamMjMzNXv2bENzUAoEAACtQnR0tPr37y9JioqKUlxcnPLz85Wdna1x48ZJksaNG6etW7canoPEymKq7g5kVyAAALXLzc3VN998owEDBqiwsFDR0dGSpK5du6qwsNDwuJQCLaay/HfHHXeYHAnQup3yuc0OAQbRK9D6Tp06pdTUVD366KOKiqq+a95ms8lmsxkemxUri/J4PIqNjTU7DAAAWhS3263U1FSNHTtWI0aMkCR16dJFBQUFkqSCggJ17tzZ8PgkVhbmcDjMDgEAgBbD5/Np/vz5iouLU3Jysv96QkKCNm7cKEnauHGjrr/+esNzUAq0sHP7BSL49NBps0OAQWUVlAKBluazzz5TVlaW+vTpo6SkJEln+gdOnTpVM2bM0Pr169WjRw8tX77c8BwkVhZFGRAAgOoGDRqkb7/9ttbnVq9eHZA5KAVaVF5eXrX2NgAAoOmxYmVRHo9HaWlp/sc2R5iJ0cCooz7ukwOa22lPudkhIIixYmVRlAIBAKhu3rx5GjJkiMaMGeO/dvz4cSUnJ2vEiBFKTk7WiRMnGjUHiZVFxcbGckAoAABVjB8/XitXrqx2LTMzU0OGDNGWLVs0ZMgQZWZmNmoOSoEW5nK5/J/7XOxQCkZRPnqWBau2IZTfgZZm8ODBys3NrXYtOztbr776qqQz7Wzuvvtuw30CJVasLMvhcCgnJ8fsMAAAaNEC2c5GIrGyrJSUFO6zAgDgAjS2nY1EKdDSqp687iunpAQAwLkq29lER0c3up2NxIqVpXHyOgAAdQtkOxuJFatWI3zcCLNDgAHfv7/P7BBg0MY72pgdAgx6Yv01ZoeAJjJr1izt3btXx44d07BhwzR9+vSAtrORSKwAAEArsWzZslqvB6qdjUQp0NI4xwoAgObFipWFVT3HCsEpqsJndggwKoxfr8EqTI3bFYbWjRUrC+McKwAAmheJlYXl5eWZHQIAAC1Gbb0Cf/vb32rUqFEaO3aspk2bpuLi4mpf8+OPP+qKK67QSy+91KA5WKu2sCuvvNLsEAAg6BxVudkhoImMHz9ekydP1pw5c/zXhg4dqrS0NIWGhup3v/udVqxYUa2lzeLFi3Xttdc2eA4SK4txOp3+e6soBQIAcFZtvQKvuebs8RoDBw7Ue++953+8detW9ezZUxEREQ2eg1KgxbhcLqWlpSktLY1SIAAAF+Avf/mLhg0bJkk6deqU/vjHP17wYdusWFmYx+M5+yCyvXmBwDBXI3tWAbhwkfxpbJVeeOEF2e123XzzzZLOVICmTJmiyMjICxqH/3osjCbMAADUb8OGDdqxY4dWrVrlb8L8xRdfaPPmzVq6dKmKi4sVEhIih8OhyZMn1zkWiZWFUQoEAKBuu3bt0sqVK/Xaa6+pbdu2/utr1qzxf/7cc88pIiKi3qRKIrGytO7du5994DptXiBAK2SLurDyAVoSdgVaVW29AjMzM1VeXq7k5GRJ0oABA7Rw4ULDc5BYWUTlbsCqOwEpBQIAcFZtvQInTJhQ79dNnz69wXOwK9AiKncDVk2mHA6HiREBAND6sGJlYdW2iOb9YF4gMMzho1dgsLKF8+s1WEXIU/+LgPNgxQoAACBASKwsyul0KiMjw+wwAABoMWrrFfjuu+9q9OjR6tu3r7788kv/9dzcXF1++eVKSkpSUlKSHn/88QbNwVq1RVXec+UXFm5eMDCMA0KDl/tryu/BqlQXmR0CmkhtvQL79Omj5557TgsWLKjx+p/85CfKysq6oDlIrCzG4XAoIyODPoEAAJyjtl6B8fHxAZ2DUqDFpKSk+HcHOp1Os8MBACBo5ebmaty4cZo8ebI+/fTTBn0NK1YW5nK5zj5wc+BdMIqqYFdgsLJ34YDQYBUmSvCQoqOjtX37dnXq1ElfffWVpk2bpk2bNikqKqrOr2PFyqIcDoe2bNlidhgAAASl8PBwderUSZJ02WWX6Sc/+YkOHTpU79exYmVRKSkp2r59u/+xz+02MRqg9ak4SRupYOVW2/pfBMsrKipShw4dZLfblZOTo8OHDzeoowmJlYXR0gYAgLNq6xXYsWNHPfXUUyoqKtL999+vfv366aWXXtInn3yi3//+9woNDVVISIiefPJJdezYsd45SKwsrFevXmaHAABAi1Fbr0BJSkxMrHFt5MiRGjly5AXPQWJlYVVb2tgiuJE2GH0bXmR2CDAosdxrdggw6KjPVf+LgPPg5nUAAIAAIbGyKFraAABQXW0tbY4fP67k5GSNGDFCycnJOnHihCRp69atGjt2rJKSkjR+/HjOsWrtarS04RyroBTh4zydYMU5VsGrXJRxraq2ljaZmZkaMmSIpk6dqszMTGVmZmr27NkaMmSIrr/+etlsNu3fv18zZszQe++9V+8cJFYtnNPprH7Q53nQwgYAgLrV1tImOztbr776qiRp3LhxuvvuuzV79mxFRp59c1RWViZbA3u3kli1cDVWns6jatnP6XQqJydHTqez2g3sAACgusLCQkVHR0uSunbtqsLCQv9z77//vjIyMlRUVKQVK1Y0aDwSKwtyuVyKjY2tvtIV2c68gGCYw1dgdggwyMeuwCDG7cetlc1mq7YylZiYqMTERH3yySd69tlntWrVqnrH4L8ei3I4HJQHAQCoR5cuXVRQcOZNbEFBgTp37lzjNYMHD1ZOTo6Kiuo/AofEyqJSUlKUl5dndhgAALRoCQkJ2rhxoyRp48aNuv766yVJ//d//yefzydJ+sc//qHy8nJ/78C6UAq0sO7du599UM6Bd8EozGd2BABgHbW1tJk6dapmzJih9evXq0ePHlq+fLkkafPmzcrKylJoaKjatGmjZ555pkE3sJNYWRi9AgEAOOt8LW1Wr15d49rUqVM1derUC56DxMpCKo9mqLy3yuFwmBwRAACtC4mVhVQezVB59EK1oxY61rwZDy3f96E/mB0CDArpzE7c4HXK7AAQxLh5HQAAIEBIrCxo7969ysjIoFcgAABVXEivwEr79u3Tf/zHfzSonY1EKdCSPB5PzdPaS1naDkaxXt77BC8OCA1WXW1tzA4BTeRCegVKktfr1dKlSzV06NAGz8FvbQtiNyAAADUNHjxYHTp0qHYtOztb48aNk3SmV+DWrVv9z7366qsaOXKkunTp0uA5SKwsqDKxcjqdJkcCAEDLdr5egfn5+dq6dasmTpx4QeNRCrSwar0COSA0KJXw1ido+UpKzQ4BBrkVZnYIMEnVXoHp6el6+OGHFRJyYb+ISawsyOFwKCMjg16BAADUo7JXYHR0dLVegV999ZVmzZolSTp27Jh27typ0NBQ3XDDDXWOR2JlQZXnV7ErEAAuXBh3ybQqlb0Cp06dWq1X4LZt2/yvmTt3rn75y1/Wm1RJ3GNlaZy8DgDAWbNmzdKdd96pQ4cOadiwYVq3bp2mTp2qDz/8UCNGjNDu3bsNtbGpihUri3A4HDpw4EC1a9VOXgcAoJW7kF6BVS1evLjBc5BYWURKSkrdpT9a2gSlqIpcs0OAQa5vj5sdAgz6sYJ2RDCOUiAAAECAkFhZlNPp5OZ1AACqqK2lzbvvvqvRo0erb9+++vLLL6u9fsWKFUpMTNTIkSP1wQcfNGgOSoEW5XK5qre1oaVNUIqoMDsCGBXSxmZ2CDCoYwgbf6yqtpY2ffr00XPPPacFCxZUe+2BAwe0adMmbdq0Sfn5+UpOTtbmzZtlt9vrnIMVKwtyOp2cYQUAwDlqa2kTHx+vuLi4Gq/Nzs7W6NGjFR4ertjYWP30pz/Vvn376p2DxMqCXC6XYmNjaWkDAIBB+fn56tatm/9xTEyM8vPz6/06SoEWVq2lDYJSKW99gpYtnH88oDXiJ9+iHA4H5UAAAAyKiYnRkSNH/I/z8/MVExNT79eRWFlUSkqKevXqZXYYAAAEpYSEBG3atEnl5eXKycnR4cOHdfnll9f7dZQCG8HpdDZ5ua0xq07VTl5v3ykA0QBoKHuXCLNDgEFhYkenVc2aNUt79+7VsWPHNGzYME2fPl0dO3bUU089paKiIt1///3q16+fXnrpJfXu3Vs33nijbrrpJtntdj3++OP17giUSKwapcaRBk3gQs6icjgcysjIoAQIAEAtztfSJjExsdbrv/71r/XrX//6guagFGghKSkpSktLU2xsrCSxKxAAgGbGipVFOZ3O6k2ZHW3MCwaGhfnMjgBG+Vxus0OAYfy+hHGsWFlU5VlWAACg+ZBYWZjDQVsGAAAq1dYrsNLLL7+sSy+9VEVFRZKkgwcP6o477tBll12ml156qcFzUAq0IIfDoQMHDmj58uVmhwK0Wt4iDugNVpFiF7VV1dYrUJLy8vL04YcfqkePHv5rHTt21Pz585WdnX1Bc7BiZUEpKSmUAQEAOEdtvQIladGiRZo9e7ZstrNHbXTp0kWXX365QkMvbA2KxMrC2BUIAEDdtm7dqujoaPXt2zcg41EKtKjKcqDfqWLzgoFhbs4pDFr0CgRavrKyMq1YsUIvv/xywMbkJ9+iKAcCAFC3H374Qbm5uUpKSlJCQoKOHDmi8ePH6+jRo4bHZMXKwtgVCADA+V166aXas2eP/3FCQoLWr1+vzp07Gx6TxMrCqvYKtP30P0yMBEa5bPvMDgEGOa4baHYIMOhffz1Q/4sQlGrrFThhwoRaX3v06FHdeuutKikpUUhIiFavXq133nlHUVFRdc5BYmVhTqezeiNmAABasfP1Cqy0bds2/+ddu3bVrl27LngOEisLc7nOnqNjC6UsGIy4eT14+UpOmR0CDAqX3ewQEMS4ed2inE6ncnJyzA4DAIBWhcTKolwul3r16mV2GAAAtBi1tbQ5fvy4kpOTNWLECCUnJ+vEiROSpDfffFNjx47V2LFjdeedd2r//v0NmoNSoIVVvb/Kd6LAxEhgVFSF2RHAKN+pMrNDgEERNkqBVlVbS5vMzEwNGTJEU6dOVWZmpjIzMzV79mxdfPHFeu2119ShQwft3LlTjz32mNatW1fvHCRWFuF0OqvdU0UZEACA6gYPHqzc3Nxq17Kzs/Xqq69KksaNG6e7775bs2fP1s9//nP/awYOHKgjR440aA4SK4twuVxKS0vzP77jjjuUkZFR7RoAAKiusLBQ0dHRks7sBCwsLKzxmvXr12vYsGENGo/EyqI8Hk/1pCq8jXnBwLAISoFByxbOr9dgddxXbnYIMInNZqvWiFmSPvroI61fv15r1qxp0BjcvG5RtLMBAKB+Xbp0UUHBmfuQCwoKqp26vn//fv3P//yPnn/+eXXq1KlB45FYWVRsbKycTqfZYQAA0KIlJCRo48aNkqSNGzfq+uuvlyT9+OOPmj59upYsWaJLLrmkweOxVm1hVW9mR3Aq5a1P0GJXYPAq9XnMDgFNpLaWNlOnTtWMGTO0fv169ejRQ8uXL5ck/eEPf9Dx48f15JNPSpLsdrs2bNhQ7xwkVhblcDh04AD9rgAAqHS+ljarV6+ucS09PV3p6ekXPAfvhy0qJSWFA0IBAGhmrFhZWNUDQkO6kWQFo4IQn9khwKCQa4abHQIM6r1qh9khIIiRWLVwDodDGRkZ9b6OA0EBADAfiVULV3XVqS7nJl+VJ7FzQCgAAGfMmzdPO3bsUJcuXfT2229LOtMrcObMmfrXv/6lnj17avny5erQoYM+/vhjPfjgg7r44oslSYmJiQ36m0xiZVE1kio3OwSDUXSFrf4XoWU6/m+zI4BBp8SuQKu6kF6BkjRo0CCtWLHiguYgsbKQqv0CKQ0CAFDdhfQKNIpdgRZSuUqVlpbGAaEAADRAXb0CP//8c918882677779M9//rNB47FiZVE1zrEKc5gXDAwr4a1P8AoLNzsCGOQWu3Fbq6q9Avv3769t27YpMjJSO3fu1LRp07Rly5Z6x+DXtkWlpKTQLxAAgHqcr1dgVFSUIiMjJUnDhw+Xx+NRUVFRveORWFmYw8EqFQAAdTlfr8CjR4/K5zuzerlv3z5VVFQ0qBEzpUCLqK2FTUOPakDLdTX9HoFmFyG72SGgiVxIr8DNmzdr7dq1stvtatOmjZYtW+YvE9aFxMoiUlJSaj3LiuQKAIAzLqRX4OTJkzV58uQLnoPEysJcVVY7bG3bmRgJ0Aodr/9eDLRMpfKaHQKCGImVxXCWFQAA5uHmdYupepZVr140XgYAoNK8efM0ZMgQjRkzxn/t+PHjSk5O1ogRI5ScnKwTJ05Ikk6cOKFp06Zp7Nixuu222/Tdd981aA5WrCys6v1V3kN/NzESGFXi4ybaYOU7dszsEGBQD7Gj2qoupKXNiy++qH79+ukPf/iDDh48qIULF9Z6L9a5WLGyEIfDQfkPAIDzGDx4sDp06FDtWnZ2tsaNGyfpTEubrVu3SpIOHjyoq6++WpIUHx+vf/3rX/r3v+vvAUpiZSGVh4I6nU5lZGTU2CUIAACqO19Lm759+/pPWt+3b59+/PFHHTlypN7xKAVaUOV9VlWFdON+q2B0OJwfUaC5larC7BBgkqotbaZOnar09HQlJSWpT58+6tevn+z2+m/P4Le2xdR2UCgAAKhdZUub6OjoGi1tFi1aJEny+Xy6/vrrG9QqjlKgxaSkpKhXr16UAgEAaIDztbQpLi5WeXm5JGndunUaNGiQoqKi6h2PFSsLqu209ZMPppoQCRrrY/tFZocAg3ouppwUrLJ9h80OAU3kQlraHDx4UHPnzpUk9e7dW+np6Q2ag8QKAAC0ChfS0uaKK67Q5s2bL3gOSoEWVbkzEAAANB9WrCzq3J2BIW3q78iNloeDCoPXv8J43xqswtwczAvj+MkHAAAIEFasgljVhsvnqjzU7NzzrAAAaM0SEhIUGRmpkJAQ2e12bdiwQTNmzNChQ4ckSSdPnlS7du2UlZVlaHwSqyBW20GglXbv3l3tOW+xr7nCQgC5xb9bsOrk5d8uWLl9XrNDQBNbvXq1/7wqSf6dgJK0ePHiBh2rcD6UAi2qIYeYAQCAs3w+n959912NGTPG8BgkVhZV2TMQAABUd++992r8+PH685//XO36p59+qi5duuhnP/uZ4bEpBVrUua1tbOEmBgPDKAUCQGCtXbtWMTExKiwsVHJysuLi4jR48GBJ0ttvv92o1SqJFSvLSklJoRwIAMA5YmJiJJ3pEZiYmKh9+/ZJkjwej95//33ddNNNjRqfxMrCHA7OQAIAoFJpaalKSkr8n3/44Yfq3bu3pDObvuLi4tStW7dGzWGoFOhwODjVW1JOTo7ZIdSpas/AtsPiTYwERpV/WWJ2CDCol632o1AAmKewsFDTpk2TJHm9Xo0ZM0bDhg2TJL3zzjsaPXp0o+cwlFjV1uS3NSK5BAAgeMTGxurNN9+s9bnFixcHZA5uXreoysNDK8+y8p0qMzkiGHGxL8zsEGBQWOhps0OAUR6zA0Aw4x4ri6rr8FAAANA0SKwsyOl0Kicnh3OsAAA4x6pVqzR69GiNGTNGs2bN8i9EjBw5UmPGjNG8efPkdrsNj08p0IJcLpdiY2PP20cQQNMrLaeMG6wi7Oyotqr8/Hy98soreuedd9SmTRs99NBD2rRpk26++WYtXbpU0pkeu+vWrdOkSZMMzUFiBQAAWg2v16vTp08rNDRUp0+fVnR0tK655hr/85dffrny8/MNj08p0KIcDoe2bNlidhgAALQYMTEx+tWvfqXrrrtO11xzjaKioqolVW63W1lZWbr22msNz8GKlUWlpKRo+/btZoeBRjJe5YfZwkIqzA4BRtFJyrJOnDih7OxsZWdnq127dnrooYeUlZWlpKQkSdKTTz6pQYMGadCgQYbnYMXKwmhpAwDAWbt379bFF1+szp07KywsTCNGjNDf//53SWc2fhUVFWnevHmNmoMVqxak8uyphqrv5PdevXo1NiQAACyjR48e+uKLL1RWVqY2bdpoz549uuyyy7Ru3Tr99a9/1apVqxQS0rg1JxKrFuRCz56qevJ71aSsMuGqekK+r6Q0QFGiOXXwtTM7BBh01MPOsmAVFmo3OwQ0kQEDBmjkyJG65ZZbFBoaqn79+umOO+7QwIED1aNHD91xxx2SpMTERMNdZkisLKJqUkarHQAAapeamqrU1NRq177++uuAjc89VhbldDpJsAAAaGasWFnQ3r171b17dy1fvvzsxTD+qYPRCRvbk4KV28b71mB1wsOtEzCOn3wL8ng87AgEAMAELGNYlMPhUEZGBo2YAQCoIiEhQZGRkQoJCZHdbteGDRu0f/9+LViwQKWlperZs6eWLl2qqKgoQ+OTWFnUubsZPD8cMykSNEaYupsdAgzqaqNXZ7C6KNTYH1QEj9WrV6tz587+x/Pnz9ecOXN05ZVXav369Vq5cqVmzJhhaGxKgQAAoFU7fPiwBg8eLEkaOnRoo1rCkVhZUGxsLLsCAQA4j3vvvVfjx4/Xn//8Z0lS7969lZ2dLUl67733lJeXZ3hsSoEWFBsbW+OwUXuXCBMjglEFNq/ZIcAgegUGL7ePfzsrW7t2rWJiYlRYWKjk5GTFxcUpPT1d6enpev7555WQkKDw8HDD45NYAQCAViMmJkaS1KVLFyUmJmrfvn2699579fLLL0uSDh06pB07dhgen1KghTmdTrNDAACgxSgtLVVJSYn/8w8//FC9e/dWYWGhJKmiokIvvPCC7rzzTsNzsGJlQQ6HQ2+99Zb69etndihopAje+wBAwBQWFmratGmSJK/XqzFjxmjYsGFavXq11qxZI+lMn8Bbb73V8BwkVhaUkpKi7du3c0goAABVxMbG6s0336xxfcqUKZoyZUpA5uDtsIU5HA6zQwAAoFVhxcqiYmNjqx0Sagu3mxgN0Pq4K3jfGqxKfRzuCuNIrIJYZdsaScrJyan2XOVZVueewA4AQGu2atUqrVu3TjabTX369NGiRYtUUFCgWbNm6fjx4+rfv7+WLFli+MgFEqsgVjVpqu0wUJfr7LsuXznnIQUjt3xmhwCDPD5WrICWJj8/X6+88oreeecdtWnTRg899JA2bdqknTt36p577tHo0aP1+OOPa/369Zo0aZKhOfjJBwAArYbX69Xp06fl8Xh0+vRpde3aVR999JFGjhwpSbrlllv8p7AbQWJlUQ6Ho1G9jgAAsJqYmBj96le/0nXXXadrrrlGUVFR6t+/v9q3b6/Q0DNFvG7duik/P9/wHJQCLaryyAUA5nBTCgxabnHrhFWdOHFC2dnZys7OVrt27fTQQw/pgw8+COgcJFYWxjlWAACctXv3bl188cXq3LmzJGnEiBH629/+puLiYnk8HoWGhurIkSP+tjdG8JbKwnr16mV2CAAAtBg9evTQF198obKyMvl8Pu3Zs0e9evXSVVddpc2bN0uS3njjDSUkJBiegxUrC6u6a9D17UkTI4FRYWpvdggwKL57kdkhwKD4oovMDgFNZMCAARo5cqRuueUWhYaGql+/frrjjjv0y1/+UjNnztTy5cvVr18/TZgwwfAcJFYAAKDVSE1NVWpqarVrsbGxWr9+fUDGJ7GyGKfTWe38qrS0NBOjAQCgdSGxshiXy1VrMtXmPzuaEA0aK2q/zewQYFBJMb06g9W/vaVmh4Agxs3rAAAAAUJiZWFOp9PsEAAAaFFWr16tMWPGaPTo0Vq1apUkaf/+/brjjjs0duxYPfDAAyopKTE8PqVACztw4MDZB2H8UwejMFEKDFYeL+9bgZbmu+++07p167Ru3TqFhYXpvvvu03XXXaf58+drzpw5uvLKK7V+/XqtXLlSM2bMMDQHP/kWlpOTY3YIAAC0GAcPHtTll1+utm3bKjQ0VIMHD9aWLVt0+PBhDR48WJI0dOjQRrWEI7GysMq+RwAAQOrTp48+++wzHTt2TGVlZdq1ay9vB40AACAASURBVJeOHDmi3r17+xsvv/fee8rLyzM8B395LcLhcCgjI6PaKtWVV1559gVujwlRobGiKsyOAIbxthVoceLj43Xffffp3nvvVdu2bdW3b1+FhIQoPT1d6enpev7555WQkKDw8HDDc5BYWUTlKesZGRn+s6woBQIAUN2ECRP8J6svW7ZMMTExio+P18svvyxJOnTokHbs2GF4fN5TWVDlWVaNWcoEAMCKCgsLJUk//vijtmzZorFjx/qvVVRU6IUXXtCdd95peHxWrCzM46lS/mNXINCsQu3UcYOV2+01OwQ0oenTp+v48eMKDQ3VggUL1L59e61evVpr1qyRJCUmJurWW281PD5/bS0sNjbW7BAAAGhRKhOoqqZMmaIpU6YEZHxKgRZGKRAAgObFipWFde/e3ewQ0EguzgcFml2YzW52CAhirFhZWK9evcwOAQCAFqW4uFipqakaNWqUbrzxRv3973+XJL366qsaNWqURo8erSVLlhgenxUri3E4HP5WNpVHMEjiHKsg5fCZHQGMimrvMjsEGNS2MMzsENCE0tPTde211+r3v/+9ysvLdfr0aX300UfKzs7Wm2++qfDwcP8uQSNYsbKYlJQUbloHAKAWJ0+e1CeffKLbbrtNkhQeHq727dtr7dq1mjp1qv9g0C5duhieg8TKwpxOp9khAADQYuTm5qpz586aN2+exo0bp/nz56u0tFSHDx/Wp59+qgkTJmjy5Mnat2+f4TkoBVqYy3W2FFFx8rSJkcC4dmYHAACW4fF49PXXX+uxxx7TgAED9Jvf/EaZmZnyer06ceKEXn/9dX355ZeaMWOGsrOzZbNd+A4iVqwsyul00tIGAIAqunXrpm7dumnAgAGSpFGjRunrr79WTEyMEhMTZbPZdPnllyskJETHjh0zNAeJlUW5XC52BQIAUEXXrl3VrVs3ff/995KkPXv2KD4+XjfccIM+/vhjSWd6BbrdbnXq1MnQHJQCLazqrsCQzlEmRgKjSnjrAzS7sgq32SGgCT322GN6+OGH5Xa7FRsbq0WLFqlt27Z69NFHNWbMGIWFhWnx4sWGyoASiVWjOBwOZWRkBGy8QJXuqh65AAAAzurXr582bNhQ4/rSpUsDMj6JVSNUOycqAAKVpKWkpCgjI0NOpzPgMQIAgPMjsbKwqrsCfac4rDAYlaqN2SHAIK+HOm6wcvvKzQ4BQYyffItyOBzsCgQAoJmRWFlUSkoKuwIBAKhi3rx5GjJkiMaMGeO/9u6772r06NHq27evvvzyyxpf8+OPP+qKK67QSy+91KA5KAVaWPVdgRw0GYwiZGxXCszXuR/l92A17JMeZoeAJjJ+/HhNnjxZc+bM8V/r06ePnnvuOS1YsKDWr1m8eLGuvfbaBs9BYgUAAFqFwYMHKzc3t9q1+Pj4875+69at6tmzpyIiIho8B6VAi3I6nQE9CgIAgNbk1KlT+uMf/3jBu+tZsbIol8ultLQ0/2NfSamJ0cA4SrjBylvsMzsEGHRU7ArEmQWKKVOmKDIy8oK+jsQKAADgHF988YU2b96spUuXqri4WCEhIXI4HJo8eXKdX0diZRFOp7PauVU5OTkcEAoAgEFr1qzxf/7cc88pIiKi3qRKIrGyjHNLfxkZGdUSLYXxTx2Mjtn+v/buPSjK6/wD+He5vSKiBRSDZB1FQFOqEiJtohUVxdgIQuiYNK02sVpsK1IQL0RjjRCvEUviOla0jniJ8YZX1CJIZCqKEjNipkkNXqIIQsPFC+ACy/n94Y9XVsHAuvqyL9/PjDP7vvvuOY97duHhnPec06B0CGQijZ3SERDRo2bNmoWzZ8+ioqICAQEBmDlzJn7yk58gISEB5eXlmD59Ol566aVWL63QHP62VSnuF0hERGRs9erVzZ4PCgp64utmzpzZ6jo4K1ClIiMjUVxcrHQYREREHQp7rFSsvr5efqyROC5hiboILhBqqaw6se0sVZ3gEDyZjj1WKqbVapUOgYiIqN1obkubRps2bUL//v1RXl4OADh48CBCQkIQEhKC3/zmN/j2229bVQcTKxVjYkVERPRQeHg4Nm7c+Nj54uJinDp1Cr16PdzO6MUXX8S2bdtw6NAh/PnPf8bChQtbVQcTKxWTJEnpEIiIiNoNf39/dOvW7bHzy5Ytw5w5c6DRPBzC9/Pzk6/19fXFrVu3WlUHEysV4xpWRERET5aRkQFXV1cMGDCgxWv27NmDgICAVpXHm9ctVHMLgj6JxsH+WYdEz0Ch5p7SIZCJGu5zSxtLVdrALcA6ipqaGqxfvx6bNm1q8ZozZ85gz549RguGPgkTKwvV3IKgRERE1HrXr19HYWEhQkNDAQC3bt1CeHg4du/ejR49euDbb7/FBx98gA0bNsDJyalVZTKxUqnGHq2myRcRERE91L9/f5w+fVo+DgwMxJ49e+Ds7IyioiLMnDkTK1euRN++fVtdJhMrlXosqerStt25qb3gUKClunuDa8cRtTfNbWkzceLEZq9du3YtKisrsXjxYgCAtbU1UlNTf7QOJlYq09hT9WP3XBEREXU0LW1p0+jEiRPy4yVLlmDJkiVtroOzAlWmsadKq9VCp9MpHQ4REVGHwh4rFWs6a1Bja6tgJGQqO/7tY7EM9Ww7S1XTUKd0CGTB+M1XKUmSOBxIRET0nDGxUqnIyEh4enoqHQYREVG70dxegUlJSQgJCUFoaCj+8Ic/oKSkxOg1+fn5+OlPf4pjx461qg4OBapY05XXxb0qBSMhUzkJa6VDIBP1HNdJ6RDIRLafG5QOgZ6R8PBwTJo0CfPmzZPPTZs2DdHR0QCALVu2YO3atYiPjwcAGAwGrFq1CsOGDWt1HeyxUglJkpCYmMjhPyIiohY0t1dgly5d5Mc1NTVG+wVu3boVr7/+OlxcXFpdB3usVKKxdyoxMdFouxsuEEpERPRkf//737F//344Ojpiy5YtAICSkhJkZGRgy5YtuHjxYqvLYmKlQs2tuK7hAqEWSa+5q3QIZCL9fyuVDoFM1MvaWekQ6DmLiYlBTEwM1q9fj23btiEqKgpLlizB7NmzYWXVtsE9JlZEREREAEJCQhAREYGoqCh8/fXXmDVrFgCgoqICJ0+ehI2NDcaMGfPEMphYqZhOpzO6gZ2IiIiMXbt2DX369AEAZGZmwsPDA4DxKuxxcXEYOXLkjyZVABMrVWu6QChZpi5C8+MXUbskapWOgExVC84KVKvm9grMzs7G1atXodFo4O7uLu8NaComViojSRIKCgqg0+k4Q5CIiKiJ5vYKbGkT5qaWL1/e6jq43ILKREZGQqvVQq/Xo7i4WOlwiIiIOhT2WKmYm5ub/FjUce8rS2QLDgVaKgOHAok6JCZWKqbVapUOgYiIqN24cuUKYmJi5OMbN24gKioKJSUlyMrKgq2tLXr37o1ly5aha9euJtXBoUAVatyAWZIkpUMhIiJqNzw8PHDgwAEcOHAAqampsLe3R1BQEIYNG4bDhw/j0KFD6NOnD9avX29yHeyxUqHIyEgkJiYaLbWgsbVVMCIyVYWmQekQyETWdkpHQKaqbOCM6o7g9OnT0Gq1cHd3h7u7u3ze19e31RsuN4c9Viqm0+mUDoGIiKhdSktLQ3Bw8GPn9+7di4CAAJPLZY+VijVdx0pUcHsNS8R1rCyXdVe2naWy01grHQI9Y7W1tThx4sRj27+tW7cO1tbWmDBhgsllM7EiIiKiDiU7Oxs+Pj7o3r27fC41NRVffPEFNm/eDI3G9D+MOBSoUpIkIT09XekwiIiI2p20tDSMHz9ePs7OzsbGjRuxbt062NvbP1XZ7LFSqcjISGRlZcnHGqefKBgNmcoWpUqHQNTh1ApuaaNm1dXVyMnJQXx8vHwuISEBtbW1mDJlCgBg8ODBRs+3BRMrFeM6VkRERMY6d+6M3Nxco3PHjx83W/lMrNoRSZKQmJjYqmtbsw+gp6fn04ZEREREbcDEqh1puu7Uj2lNAmZUXo8XTAmJFOZk4FCgpaq5yVmBRB0Rb14nIiIiMhP2WKmYTqdrUy8YERGR2gUGBsLBwQFWVlawtrZGamoqKisrERMTg5s3b8Ld3R1JSUno1q2bSeUzsVKxpguEolv3li+kdstWKB0BUcdTz1mBqpeSkgJnZ2f5ODk5Ga+99hoiIiKQnJyM5ORkzJkzx6SyORSoQjqdDomJia26wZ2IiKijy8zMRFhYGAAgLCwMGRkZJpfFxEqF9Ho9YmNjUVxcrHQoRERE7c7UqVMRHh6OnTt3AgDKysrg6uoKAOjRowfKyspMLptDgSrm5ub28KC+VrlAyGR1nFhmsewcG5QOgUxko7FTOgR6hnbs2IGePXuirKwMU6ZMgYeHh9HzGo2GW9pQ87hAKBERkbGePXsCAFxcXBAUFIT8/Hy4uLigtPTB8jalpaVG91+1FRMrFZMkSekQiIiI2o3q6mrcu3dPfnzq1Cl4eXkhMDAQ+/fvBwDs378fo0ePNrkODgWqmNFSC1V3lAuEqAO6fbOT0iGQiTgrUL3KysowY8YMAIDBYEBwcDACAgIwcOBAREdHY8+ePejVqxeSkpJMroOJlQrodDqjpRU4G5CIiOhxWq0WBw8efOy8k5MTUlJSzFIHEysVaJwF2KhxuxsuEEpERPR8MbFSKZ1Oh4KCgocn6jgr0BJxViARkWVhYqVSer2eswKJiIiauHLlCmJiYuTjGzduICoqCu+99x62bt2K7du3w9raGiNGjMDcuXNNqoOJlYpxViAREdFDHh4eOHDgAIAHN68HBAQgKCgIZ86cQWZmJg4ePAg7OzsuEErGJElCQUGB8awG+y7KBUQmu8cFUSxWlx73lQ6BTFVpq3QE9BycPn0aWq0W7u7uWLlyJSIiImBn92BxWBcXF5PL5Y9tFYqMjIRWq4VOp1M6FCIionYpLS0NwcHBAIBr164hLy8PEydOxKRJk5Cfn29yueyxUrGmSzDgToVygRB1QLX3+OPVUtlorJUOgZ6x2tpanDhxQp5RbzAYcPv2bezatQsXL15EdHQ0MjMzTdrahj1WRERE1KFkZ2fDx8cH3bt3B/Bgm5ugoCBoNBoMGjQIVlZWqKgwrUOCiZVKSZKE9PR0pcMgIiJqd9LS0jB+/Hj5eMyYMcjNzQUAXL16FXV1dXBycjKpbPZVq1RkZCSysrIenrDlzZiWyFYoHQFRx1PTUKd0CPQMVVdXIycnB/Hx8fK5X//615g/fz6Cg4Nha2uL5cuXmzQMCDCxUjWuY0VERGSsc+fOcu9UIzs7O6xatcos5XMoUMU8PT2VDoGIiKhDYY+VijXdJ1AUcmNmS8QtbSyXtW2D0iGQieoE245Mxx4rIiIiIjNhYqViXCCUiIjooStXriA0NFT+5+fnh82bN+Obb77BW2+9hdDQUISHh3OBUGqe0QKhdZzlYok4K9By1Vbzx6ul6mrNPge1ammvwIULF2LGjBkYMWIETp48iY8//hhbt241qQ5+eoiIiKjDabpXoEajQVVVFQDg7t27cHV1Nblc/klloSRJQmJiIgDgxg3jG9N1Oh30ej3S09Pl5fqJiIjooaZ7Bc6fPx9Tp07FihUr0NDQgM8//9zkcplYWaimM/4aE6xGer0esbGxyMnJeXiSC4RaqFqlAyAT2XWuVzoEMlH9D9wrUO0e3Stwx44deP/99/H666/jyJEjWLBgATZv3mxS2RwKVDEuEEpERPS4R/cK3LdvH8aOHQsA+NWvfvVUN68zsVKx4uJipUMgIiJqdx7dK9DV1RVnz54FAJw5cwZ9+vQxuWwOBaqYm5ub0iHQU5I4K9BiGer4dytRe9TcXoEJCQlYunQp6uvrIUmS0XNtxcRKxTgUSEREZKy5vQKHDBmC1NRUs5TPxEpFGmcDNs4SlCRJ4YiIiIg6FiZWKtI4G7BxlmDTmYMaJyelwqKnoNdUKR0CmajqNv+wsVR14r7SIZAFY2JFREREHUZKSgp2794NIQQmTpyI9957D9HR0bh69SqABwuEOjo6yiu0txUTKxXT6XRGvVZEREQd2aVLl7B7927s3r0btra2mDZtGkaNGoWkpCT5muXLl6NLly4m18HESqV0Oh0KCgrkY1HNISUiIurYLl++jEGDBsHe3h4A4O/vj/T0dPzxj38EAAghcPToUaSkpJhcB+cDq5Rer+esQCIioia8vb3x5ZdfoqKiAjU1NcjOzsatW7fk5/Py8uDi4sJ1rKh5TWcFivLbCkZCpuI6Vparvp5/t1oqeytuAaZW/fr1w7Rp0zB16lTY29tjwIABsLJ6+F09fPiwvH+gqfjNVyFJknDjxg3eX0VERPSIiRMnIjU1Fdu3b0e3bt3k3qn6+nocP34cb7zxxlOVz8RKhSIjI6HVaqHT6ZQOhYiIqF0pKysDABQVFSE9PR0hISEAgJycHHh4eOCFF154qvI5FKhier1efqxxsFcwEjLVPStOOrBUNXoOJ1mqGlQrHQI9QzNnzkRlZSVsbGywaNEidO3aFQBw5MgRo/0DTcXESqUkSTKaFUhERETAZ5991uz55cuXm6V8DgWqgCRJSExMlLeyAR4MB3p6eioYFRERUcfDHisVaLxJvXErm0fPA4CoqnmuMZF5dGlQOgIylefQCqVDIBO5ne2qdAhkwdhjRURERGQmTKxUjLMCiYiIHrpy5QpCQ0Plf35+fti8eTMqKysxZcoUjB07FlOmTMHt26av/cihQBVrOiuQLFNnDgVarNofuLorUXvj4eEhb65sMBgQEBCAoKAgJCcn47XXXkNERASSk5ORnJyMOXPmmFQHe6yIiIiowzl9+jS0Wi3c3d2RmZmJsLAwAEBYWBgyMjJMLpeJlQrpdDokJiYiPT1d6VCIiIjapbS0NHn7mrKyMri6ugIAevToIS8iagoOBaqQXq9HbGwscnJy5HNcINQyVXOBUItlbad0BETUktraWpw4cQKxsbGPPafRaKDRaEwumz1WKqbVapUOgYiIqN3Jzs6Gj48PunfvDgBwcXFBaWkpAKC0tBTOzs4ml83ESsWKi4uVDoGIiKjdSUtLM9q+JjAwEPv37wcA7N+/H6NHjza5bA4Fqpibm9vDg7o65QIh6oAMtUpHQETNqa6uRk5ODuLj4+VzERERiI6Oxp49e9CrVy8kJSWZXD4TKxV5dGsbDgUSEREZ69y5M3Jzc43OOTk5ISUlxSzlM7FSkUe3tpEkSclwiIiIOhwmVirWdK9Askyu9Vxk0lLVVHBaoKWqFvVKh0AWjIkVERERdRiBgYFwcHCAlZUVrK2tkZqaiqNHj0Kn0+Hy5cvYvXs3Bg4caHL5TKxUSqfTyetZERER0UMpKSlGSyp4e3tjzZo1WLRo0VOXzcRKpR5NqjT/v6IsWZZ7VoVKh0DU4dQKg9Ih0HPWr18/s5XFdaxUSJIkeWYgERERGZs6dSrCw8Oxc+dOs5fNxEqFIiMjodVqodPplA6FiIioXdmxYwf27duHDRs2YPv27Th37pxZy+dQoEpJkoSCgoKHJ+q4WqElkgRnBVoqu86cWWap7Mq5t6qa9ezZE8CDbWyCgoKQn58Pf39/s5XPHiuVioyMhKenp9JhEBERtRvV1dW4d++e/PjUqVPw8vIyax3ssVIxrmNFRET0UFlZGWbMmAEAMBgMCA4ORkBAAI4fP46EhASUl5dj+vTpeOmll/DPf/7TpDo0QnCsgYiIiMgcOBRIREREZCZMrIiIiIjMhIkVERERkZkwsSIiIiIyEyZWRERERGbCxIqIiIjITJhYEREREZkJEysiIiIiM2FiRe1ebm4upk+f3uLzqampiI+PN3u9qampKCkpkY8DAwNRXl5u9no6gh9rwx9z8eJFfPTRR80+19gud+7cwfbt281Wp1o9+rluSVxcHI4dO9bi85MnT8bFixfNGRrbsJXM1YY/5pNPPkFOTs5j55u2S25uLs6fP2+2OtWAiRVRC/bt24fS0lKlwyAAAwcOxAcffPDEa+7cuYMdO3Y8p4gsV3v+XLMNW+d5teFf//pXDB069InXnD17Fl999dUzj8WScK9AMovq6mpER0fj1q1baGhowF/+8hf07t0by5cvR3V1NZycnLBs2TK4urpi8uTJ6N+/P86dOweDwYClS5di0KBByM/Px5IlS6DX69GpUycsXboUHh4ebYqjvLwcixYtQlFREQBg/vz5eOWVV7BmzRoUFRWhsLAQRUVFePfdd/H73/8eALB27VocPHgQzs7OcHNzg4+PD9zd3fH1119j9uzZ6NSpE3bu3AkA2LZtG7KyslBfX4+kpCT069fPvG+kgpRsw5CQEGzfvh2Ojo549dVX8f777yMsLAxz585FaGgobGxssGnTJqxfvx4VFRWIjY1FSUkJfH190bgrV2JiIq5fv47Q0FAMHToUI0eORHV1NaKionDp0iX4+Phg1apV0Gg0z/qtfK4KCwsxbdo0+Pj44D//+Q+8vLywYsUKXL58+bG2O3/+/GOf640bNyIrKwt6vR4vv/wy4uPj2/we/fvf/8aaNWtQW1sLrVaLZcuWwcHBAYGBgQgLC3vsO1NeXo7Y2FiUlpbC19cXOTk52Lt3L9vwObZhfn4+kpOTodPpkJGRgVmzZiEvLw9CCLzxxhvIzMxEXFwcRo4ciXHjxiE7OxtLly6Fvb09XnnlFTnuzz//HFZWVjh48CAWLlwIAMjLy8PmzZvxv//9D3PmzMG4ceOe+XvYrggiMzh27JhYsGCBfHznzh3x9ttvi7KyMiGEEGlpaSIuLk4IIcSkSZPka8+ePSvGjx8vhBDi7t27oq6uTgghxKlTp0RkZKQQQogzZ86IiIiIFuveu3evWLx4sRBCiFmzZolz584JIYS4efOmGDdunBBCiE8//VS8/fbbQq/Xi7KyMvHzn/9c1NbWigsXLogJEyaI+/fvi7t374qgoCCxceNGOc78/Hy5nlGjRoktW7YIIYTYtm2bmD9//tO8Ze2Okm24cOFCkZWVJf773/+K8PBwueygoCBRVVVl9PqEhASxZs0aIYQQWVlZwtvbW5SVlYkbN27IcTTW6efnJ4qLi4XBYBBvvfWW/NlQkxs3bghvb2+Rl5cnhBAiLi5ObNiw4Ylt1/RzXVFRIT+ePXu2yMzMFEIIMW/ePHH06NEW620sp6ysTPz2t78VVVVVQggh1q9fL7dPS9+ZxYsXi3/84x9CCCFOnjzJNlSgDevq6kRgYKAQQojly5eL8PBwkZeXJ3Jzc0VMTIzR6+/fvy8CAgLE1atXRUNDg4iKipK/j59++qn8M7PxNTNnzhQGg0F89913YsyYMWZ5jywJe6zILLy9vbFixQp8/PHHGDVqFLp27YpLly5hypQpAICGhgb06NFDvn78+PEAAH9/f9y7dw937txBVVUV5s2bh++//x4ajQZ1dXVtjiMnJwcFBQXy8b1791BVVQUAGDFiBOzs7ODs7AxnZ2eUlZXh/PnzGD16NCRJgiRJGDVq1BPLHzt2LADgZz/7GY4fP97m+NozJdtwyJAhOHfuHHr16oV33nkHu3btQklJCbp27YrOnTsbXXvu3DnodDoAwMiRI9GtW7cWyx00aBBeeOEFAMCAAQNw8+ZNDBkypPVvioVwc3OTexEmTJiA9evXP7HtmsrNzcXGjRtx//59VFZWwsvLC4GBga2u+8KFCygoKMA777wDAKirq4Ovr6/8fHPfmS+//FJuw4CAALYhnn8b2tjYoHfv3rh8+TLy8/MxZcoU5OXlwWAwyHE0unLlCl588UX06dNHjm/Xrl0tlj1mzBhYWVnB09MTP/zwQ2vfAtVgYkVm0bdvX6SmpuLkyZNISkrCq6++Ci8vL3kI7VGPdlNrNBp88skn+MUvfoG1a9eisLBQHqpri4aGBuzatQuSJD32nJ2dnfzY2toa9fX1bS7f1tYWAGBlZQWDwdDm17dnSrahv78/PvvsMxQXFyMmJgYZGRk4duzYU/8CfbTN1dZmjR5tCwcHhye2XSO9Xo/Fixdj7969cHNzw5o1a6DX69tUtxACw4YNw+rVq5t9/mm/M2zDZ9eGQ4YMQXZ2NmxsbDB06FDExcXBYDBg7ty5Jv8/AOM264h48zqZRUlJCezt7REaGoqpU6fiwoULKC8vl29qrKurw3fffSdff+TIEQAPxuIdHR3h6OiIu3fvomfPngAe3Jxpil/+8pfYunWrfPzNN9888Xo/Pz/53oSqqip88cUX8nMODg5yb1dHoGQburm5oaKiAteuXYNWq4Wfnx82bdrUbGLl7++PQ4cOAQBOnjyJ27dvA+h47dVUUVGR3E6HDx/G4MGDW2y7pu9T4y9gJycnVFVV4V//+leb6/b19cX58+fx/fffA3hwr97Vq1ef+Bo/Pz8cPXoUwIP7s9iGyrThkCFDkJKSAl9fXzg7O6OyshJXr16Ft7e30XUeHh64efMmrl+/DgBIS0uTn+vIbdYS9liRWVy6dAkrV66ElZUVbGxs8OGHH8LGxgYfffQR7t69C4PBgHfffRdeXl4AAEmSEBYWhvr6eixduhQAMG3aNMTFxWHdunUYMWKESXEsWLAA8fHxCAkJgcFgwJAhQ564FMOgQYMQGBiICRMmwMXFBd7e3nB0dAQAvPnmm1i0aJHRzetqpnQbDho0CA0NDQAe/MBfvXr1Y0MSADBjxgzExsZi/PjxePnll9GrVy8AD36x+Pn5ITg4GMOHD8fIkSOf4t2wLH379sX27dsxf/58eHp6YvLkyRg+fHizbffo53rixIkIDg5G9+7dMXDgwDbX7ezsjGXLlmHWrFmora0FAERHR6Nv374tviYyMhKzZs3CwYMH4evrix49eqBLly6ws7NjGz7HNhw8eDB++OEH+Pv7AwD69+8PZ2fnx3rPJElCfHw8IiIi5JvXG5OpUaNGISoqCpmZmfLN6x2dRoj/n1JD9JxMnjwZc+fONemH+LNQNymZbwAAANVJREFUVVUFBwcH1NTU4He/+x0SEhLg4+OjdFjtWntrw46ssLAQf/rTn3D48GGlQ2m12tpaOYH/6quv8OGHH+LAgQNKh6UYS2xDahl7rKjD+9vf/oaCggLo9Xq8+eabTKqInrGioiJER0ejoaEBtra2SEhIUDokIrNhjxVZjL1792LLli1G5/z8/LBo0SKFIqK2YhtanhkzZqCwsNDo3OzZszF8+HCFIqK2Yhs+X0ysiIiIiMyEswKJiIiIzISJFREREZGZMLEiIiIiMhMmVkRERERm8n+Hgharariq2gAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import seaborn as sns; sns.set(color_codes=True)\n",
"iris = sns.load_dataset(\"iris\")\n",
"species = iris.pop(\"species\") # Remove the species column\n",
"print(species.unique()) # The samples seems to be from these three species\n",
"sns.clustermap(iris, method=\"ward\", col_cluster=False, cbar_kws={'label': 'centimeters'}); # Cluster only the rows\n",
"#plt.colorbar().ax.set_title('This is a title')\n",
"#plt.gca().images[-1].colorbar.ax.set_title(\"title\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With sharp eye and good will one can discern three clusters in the above heatmap and dendrogram."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercise 7 (binding sites)
\n",
"\n",
"This exercise can give three points at maximum!\n",
"\n",
"A binding site is a piece of DNA where a certain protein prefers to bind. The piece of DNA can be described as a string consisting of letters A, C, G, and T, which correspond to nucleotides Adenine, Cytosine, Guanine, and Thymine.\n",
"In this exercise the length of binding sites is eight nucleotides. They are stored in the file `data.seq`,\n",
"and the binding sites there are classified into two classes.\n",
"\n",
"Part 1. Write function `toint` that converts a nucleotide to an integer. Use the following mapping:\n",
"```\n",
"A -> 0\n",
"C -> 1\n",
"G -> 2\n",
"T -> 3\n",
"```\n",
"\n",
"Write also function `get_features_and_labels` that gets a filename as a parameter. The function should load the contents of the file into a DataFrame. The column `X` contains a string. Convert this column into a feature matrix using the above `toint` function. For example the column `[\"GGATAATA\",\"CGATAACC\"]` should result to the feature matrix\n",
"```\n",
"[[2,2,0,3,0,0,3,0],\n",
"[1,2,0,3,0,0,1,1]]\n",
"```\n",
"The function should return a pair, whose first element is the feature matrix and the second element is the label vector.\n",
"\n",
"Part 2. Create function `cluster_euclidean` that gets a filename as parameter. Get the features and labels using the function from part 1. Perform hierarchical clustering using the function `sklearn.cluster.AgglomerativeClustering`. Get two clusters using `average` linkage and `euclidean` affinity. Fit the model and predict the labels. Note that you may have to use the `find_permutation` function again, because even though the clusters are correct, they may be labeled differently than the real labels given in `data.seq`. The function should return the accuracy score.\n",
"\n",
"Part 3. Create function `cluster_hamming` that works like the function in part 2, except now using the [hamming](https://en.wikipedia.org/wiki/Hamming_distance) affinity. Even though it is possible to pass the function `hamming` to `AgglomerativeClustering`, let us now compute the Hamming distance matrix explicitly. We can achieve this using the function `sklearn.metrics.pairwise_distances`. Use the affinity parameter `precomputed` to `AgglomerativeClustering`. And give the distance matrix you got from `pairwise_distances`, instead of the feature matrix, to the `fit_predict` method of the model. If you want, you can visualize the clustering using the provided `plot` function.\n",
"\n",
"**NB!** When submitting your solution, please comment away all `plot` function calls. This might cause tests to fail on the server.\n",
"\n",
"Which affinity (or distance) do you think is theoretically more correct of these two (Euclidean or Hamming)? Why?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open")
\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}