{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# BINF6112 Tute/Lab Questions 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Last revision: Fri 25 Sep 13:18:08 AEST 2020.*" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from IPython.core.display import display, HTML\n", "display(HTML(\"\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial/lab activity, you will answer some questions on topics from lectures and learn how to apply some of the techniques we covered to representative datasets. We will be working with a number of microarray datasets and working on data preparation and cluster analysis. You need to work through the notebook and complete the answers to questions in the notebook cells.\n", "\n", "#### Please ensure all answers are completed in a single cell of the notebook and that your answers are in ```markdown``` format !!!\n", "\n", "Once you have completed the notebook, make sure you save it with the filename \"Lab1_Solutions.ipynb\" and submit it via give. For example, using the command-line version of ```give``` on a CSE Linux machine:\n", "\n", "```\n", "$ give bi6112 lab1 Lab1_Solutions.ipynb\n", "```\n", "\n", "Alternatively, you can submit using the web-based interface to ```give```.\n", "\n", "There are a total of *10 marks* available.\n", "Each notebook mark is worth *0.5 course mark*, i.e., notebook marks will be scaled\n", "to a **course mark out of 5** to contribute to the course total.\n", "\n", "All questions are worth 1 mark, except for question 7 which is worth 2 marks.\n", "\n", "Deadline: 23:59:59, Sunday October 11, 2020." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Your environment should have already been set up to enable to run this notebook if you successfully completed the lab \"Lab1_Intro_Python.ipynb\". The main data structure we will be using is the ```DataFrame``` (something you are probably familiar with from R). The principal toolkit for handling DataFrames in Python is a library called ```Pandas```. You can find some useful information in compact form on basic use of [Pandas](https://towardsdatascience.com/how-to-master-pandas-8514f33f00f6) and another one on some more [advanced uses](https://towardsdatascience.com/learn-advanced-features-for-pythons-main-data-analysis-library-in-20-minutes-d0eedd90d086)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Warm-up questions\n", "\n", "Refer to the lecture slides \"Gene Expression and Functional Genomics (1) and (2)\" and add your answer in the cell after the question." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q1. At which stage of the systems biology methodology should an experiment to obtain data from high-thoughput genome-wide assays be carried out ?**\n", "\n", "1. Define components (``parts list'') and interactions of the system\n", "2. Systematically perturb the real system and monitor results\n", "3. Reconcile experimental data with output from model simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q1 answer.** _Your answer goes here._ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q2. Why is systems biology needed ? (refer to lecture \"Gene Expression & Functional Genomics 1\")**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q2 answer.** _Your answer goes here._" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q3. What is a typical trade-off made in systems biology modelling ? (refer to lecture \"Gene Expression & Functional Genomics 1\")**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q3 answer.** _Your answer goes here._" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data Exploration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a first step, whenever we get some data given to us we will start by doing some basic exploration. This will give us an idea of what is in the data and what types of analysis we may be able to do. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(7129, 78)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns; sns.set()\n", "\n", "df_train = pd.read_csv(\"./datasets/ALL_AML_train.csv\")\n", "df_train.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see above, there are over 7,000 rows and 78 columns. From our knowledge of microrarray data we would assume that the shape of this data frame suggests the rows correspond to genes and the columns to samples. Let's take a peek." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01234
Gene DescriptionAFFX-BioB-5_at (endogenous control)AFFX-BioB-M_at (endogenous control)AFFX-BioB-3_at (endogenous control)AFFX-BioC-5_at (endogenous control)AFFX-BioC-3_at (endogenous control)
Gene Accession NumberAFFX-BioB-5_atAFFX-BioB-M_atAFFX-BioB-3_atAFFX-BioC-5_atAFFX-BioC-3_at
1-214-153-5888-295
callAAAAA
2-139-73-1283-264
call.1AAAAA
3-76-49-307309-376
call.2AAAAA
4-135-11426512-419
call.3AAAAA
5-106-125-76168-230
call.4AAAAA
6-138-8521571-272
call.5AAAAA
7-72-14423855-399
call.6AAAAA
8-413-2607-2-541
call.7AAAAA
95-127106268-210
call.8AAAAA
10-88-10542219-178
call.9AAAMA
11-165-155-7182-163
call.10AAAAA
12-67-938425-179
call.11AAAAA
13-92-119-31173-233
call.12AAAAA
14-113-147-118243-127
call.13AAAMA
..................
24-120-263-114255-342
call.23AAAAA
25-81-150-85316-418
call.24AAAAA
26-112-233-7854-244
call.25AAAAA
27-273-327-7681-439
call.26AAAAA
34-20-207-50101-369
call.27AAAAA
357-100-57132-377
call.28AAAAA
36-213-252136318-209
call.29AAAAA
37-25-20124325-396
call.30AAAAA
38-72-139-1392-324
call.31AAAPA
28-4-116-125241-191
call.32AAAAA
2915-1142193-51
call.33AAAAA
30-318-192-95312-139
call.34AAAAA
31-32-4949230-367
call.35AAAPA
32-124-79-37330-188
call.36AAAAA
33-135-186-70337-407
call.37AAAAA
\n", "

78 rows × 5 columns

\n", "
" ], "text/plain": [ " 0 \\\n", "Gene Description AFFX-BioB-5_at (endogenous control) \n", "Gene Accession Number AFFX-BioB-5_at \n", "1 -214 \n", "call A \n", "2 -139 \n", "call.1 A \n", "3 -76 \n", "call.2 A \n", "4 -135 \n", "call.3 A \n", "5 -106 \n", "call.4 A \n", "6 -138 \n", "call.5 A \n", "7 -72 \n", "call.6 A \n", "8 -413 \n", "call.7 A \n", "9 5 \n", "call.8 A \n", "10 -88 \n", "call.9 A \n", "11 -165 \n", "call.10 A \n", "12 -67 \n", "call.11 A \n", "13 -92 \n", "call.12 A \n", "14 -113 \n", "call.13 A \n", "... ... \n", "24 -120 \n", "call.23 A \n", "25 -81 \n", "call.24 A \n", "26 -112 \n", "call.25 A \n", "27 -273 \n", "call.26 A \n", "34 -20 \n", "call.27 A \n", "35 7 \n", "call.28 A \n", "36 -213 \n", "call.29 A \n", "37 -25 \n", "call.30 A \n", "38 -72 \n", "call.31 A \n", "28 -4 \n", "call.32 A \n", "29 15 \n", "call.33 A \n", "30 -318 \n", "call.34 A \n", "31 -32 \n", "call.35 A \n", "32 -124 \n", "call.36 A \n", "33 -135 \n", "call.37 A \n", "\n", " 1 \\\n", "Gene Description AFFX-BioB-M_at (endogenous control) \n", "Gene Accession Number AFFX-BioB-M_at \n", "1 -153 \n", "call A \n", "2 -73 \n", "call.1 A \n", "3 -49 \n", "call.2 A \n", "4 -114 \n", "call.3 A \n", "5 -125 \n", "call.4 A \n", "6 -85 \n", "call.5 A \n", "7 -144 \n", "call.6 A \n", "8 -260 \n", "call.7 A \n", "9 -127 \n", "call.8 A \n", "10 -105 \n", "call.9 A \n", "11 -155 \n", "call.10 A \n", "12 -93 \n", "call.11 A \n", "13 -119 \n", "call.12 A \n", "14 -147 \n", "call.13 A \n", "... ... \n", "24 -263 \n", "call.23 A \n", "25 -150 \n", "call.24 A \n", "26 -233 \n", "call.25 A \n", "27 -327 \n", "call.26 A \n", "34 -207 \n", "call.27 A \n", "35 -100 \n", "call.28 A \n", "36 -252 \n", "call.29 A \n", "37 -20 \n", "call.30 A \n", "38 -139 \n", "call.31 A \n", "28 -116 \n", "call.32 A \n", "29 -114 \n", "call.33 A \n", "30 -192 \n", "call.34 A \n", "31 -49 \n", "call.35 A \n", "32 -79 \n", "call.36 A \n", "33 -186 \n", "call.37 A \n", "\n", " 2 \\\n", "Gene Description AFFX-BioB-3_at (endogenous control) \n", "Gene Accession Number AFFX-BioB-3_at \n", "1 -58 \n", "call A \n", "2 -1 \n", "call.1 A \n", "3 -307 \n", "call.2 A \n", "4 265 \n", "call.3 A \n", "5 -76 \n", "call.4 A \n", "6 215 \n", "call.5 A \n", "7 238 \n", "call.6 A \n", "8 7 \n", "call.7 A \n", "9 106 \n", "call.8 A \n", "10 42 \n", "call.9 A \n", "11 -71 \n", "call.10 A \n", "12 84 \n", "call.11 A \n", "13 -31 \n", "call.12 A \n", "14 -118 \n", "call.13 A \n", "... ... \n", "24 -114 \n", "call.23 A \n", "25 -85 \n", "call.24 A \n", "26 -78 \n", "call.25 A \n", "27 -76 \n", "call.26 A \n", "34 -50 \n", "call.27 A \n", "35 -57 \n", "call.28 A \n", "36 136 \n", "call.29 A \n", "37 124 \n", "call.30 A \n", "38 -1 \n", "call.31 A \n", "28 -125 \n", "call.32 A \n", "29 2 \n", "call.33 A \n", "30 -95 \n", "call.34 A \n", "31 49 \n", "call.35 A \n", "32 -37 \n", "call.36 A \n", "33 -70 \n", "call.37 A \n", "\n", " 3 \\\n", "Gene Description AFFX-BioC-5_at (endogenous control) \n", "Gene Accession Number AFFX-BioC-5_at \n", "1 88 \n", "call A \n", "2 283 \n", "call.1 A \n", "3 309 \n", "call.2 A \n", "4 12 \n", "call.3 A \n", "5 168 \n", "call.4 A \n", "6 71 \n", "call.5 A \n", "7 55 \n", "call.6 A \n", "8 -2 \n", "call.7 A \n", "9 268 \n", "call.8 A \n", "10 219 \n", "call.9 M \n", "11 82 \n", "call.10 A \n", "12 25 \n", "call.11 A \n", "13 173 \n", "call.12 A \n", "14 243 \n", "call.13 M \n", "... ... \n", "24 255 \n", "call.23 A \n", "25 316 \n", "call.24 A \n", "26 54 \n", "call.25 A \n", "27 81 \n", "call.26 A \n", "34 101 \n", "call.27 A \n", "35 132 \n", "call.28 A \n", "36 318 \n", "call.29 A \n", "37 325 \n", "call.30 A \n", "38 392 \n", "call.31 P \n", "28 241 \n", "call.32 A \n", "29 193 \n", "call.33 A \n", "30 312 \n", "call.34 A \n", "31 230 \n", "call.35 P \n", "32 330 \n", "call.36 A \n", "33 337 \n", "call.37 A \n", "\n", " 4 \n", "Gene Description AFFX-BioC-3_at (endogenous control) \n", "Gene Accession Number AFFX-BioC-3_at \n", "1 -295 \n", "call A \n", "2 -264 \n", "call.1 A \n", "3 -376 \n", "call.2 A \n", "4 -419 \n", "call.3 A \n", "5 -230 \n", "call.4 A \n", "6 -272 \n", "call.5 A \n", "7 -399 \n", "call.6 A \n", "8 -541 \n", "call.7 A \n", "9 -210 \n", "call.8 A \n", "10 -178 \n", "call.9 A \n", "11 -163 \n", "call.10 A \n", "12 -179 \n", "call.11 A \n", "13 -233 \n", "call.12 A \n", "14 -127 \n", "call.13 A \n", "... ... \n", "24 -342 \n", "call.23 A \n", "25 -418 \n", "call.24 A \n", "26 -244 \n", "call.25 A \n", "27 -439 \n", "call.26 A \n", "34 -369 \n", "call.27 A \n", "35 -377 \n", "call.28 A \n", "36 -209 \n", "call.29 A \n", "37 -396 \n", "call.30 A \n", "38 -324 \n", "call.31 A \n", "28 -191 \n", "call.32 A \n", "29 -51 \n", "call.33 A \n", "30 -139 \n", "call.34 A \n", "31 -367 \n", "call.35 A \n", "32 -188 \n", "call.36 A \n", "33 -407 \n", "call.37 A \n", "\n", "[78 rows x 5 columns]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_train.head(5).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This appears to be Affymetrix data. \n", "However, it also looks like this dataset needs to be \"cleaned\" before starting any analysis, so that is what we will do next." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset sizes: \n", "Training set: (38, 7129)\n", "Test set: (34, 7129)\n", "Class distribution: \n", "ALL 47\n", "AML 25\n", "Name: cancer, dtype: int64\n" ] } ], "source": [ "# Read in the training and test datasets\n", "train_dataset = pd.read_csv('./datasets/ALL_AML_train.csv')\n", "test_dataset = pd.read_csv('./datasets/ALL_AML_test.csv')\n", "\n", "train_dataset1 = [col for col in train_dataset.columns if \"call\" not in col]\n", "train_dataset = train_dataset[train_dataset1]\n", "\n", "test_dataset1 = [col for col in test_dataset.columns if \"call\" not in col]\n", "test_dataset = test_dataset[test_dataset1]\n", "\n", "train_dataset = train_dataset.T\n", "test_dataset = test_dataset.T\n", "\n", "train_dataset2 = train_dataset.drop(['Gene Description','Gene Accession Number'],axis=0)\n", "test_dataset2 = test_dataset.drop(['Gene Description','Gene Accession Number'],axis=0)\n", "\n", "# Ensure data frame entries are numeric and indexing has the same order\n", "train_dataset2.index = pd.to_numeric(train_dataset2.index)\n", "train_dataset2.sort_index(inplace=True)\n", "\n", "test_dataset2.index = pd.to_numeric(test_dataset2.index)\n", "test_dataset2.sort_index(inplace=True)\n", "\n", "print(f'Dataset sizes: ')\n", "print(f'Training set: {train_dataset2.shape}')\n", "print(f'Test set: {test_dataset2.shape}')\n", "\n", "# Read in the classes\n", "y = pd.read_csv('actual.csv')\n", "print(f'Class distribution: ')\n", "print(y['cancer'].value_counts())\n", "y = y.replace({'ALL':0,'AML':1})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to do some pre-processing to make sure the data is in a suitable format for the sklearn algorithms, and scale the gene expression values as a normalization step." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# Training set\n", "X_train = train_dataset2.reset_index(drop=True)\n", "Y_train = y[y.patient <= 38].reset_index(drop=True)\n", "\n", "# Test set\n", "X_test = test_dataset2.reset_index(drop=True)\n", "Y_test = y[y.patient > 38].reset_index(drop=True)\n", "\n", "Y_test = Y_test.iloc[:,1].values\n", "Y_train = Y_train.iloc[:,1].values\n", "\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "sc = StandardScaler()\n", "X_train = sc.fit_transform(X_train)\n", "X_test = sc.fit_transform(X_test)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the data is very high-dimensional and we only have a relatively small number of samples we will apply _dimensionality reduction_ using PCA." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.09735759e+03 8.76976485e+02 4.83272186e+02 3.57658145e+02\n", " 3.39170426e+02 2.72510227e+02 2.55597771e+02 2.40859971e+02\n", " 2.18576598e+02 1.93681882e+02 1.83741623e+02 1.72634428e+02\n", " 1.61466953e+02 1.52866654e+02 1.42171050e+02 1.38629132e+02\n", " 1.35241393e+02 1.25441555e+02 1.24851480e+02 1.20423734e+02\n", " 1.12430908e+02 1.11931494e+02 1.06610477e+02 1.04779826e+02\n", " 1.03617441e+02 1.00934278e+02 9.86526997e+01 9.56118071e+01\n", " 9.54129156e+01 9.14541425e+01 8.49250449e+01 8.17690665e+01\n", " 7.56789591e+01 7.30684864e+01 6.76502725e+01 6.28995436e+01\n", " 6.11190279e+01 2.16415626e-28]\n", "[1.49877930e-01 1.19778111e-01 6.60056806e-02 4.88492199e-02\n", " 4.63241532e-02 3.72196529e-02 3.49097368e-02 3.28968370e-02\n", " 2.98533570e-02 2.64532179e-02 2.50955698e-02 2.35785408e-02\n", " 2.20532785e-02 2.08786432e-02 1.94178295e-02 1.89340717e-02\n", " 1.84713718e-02 1.71329025e-02 1.70523096e-02 1.64475646e-02\n", " 1.53558984e-02 1.52876881e-02 1.45609395e-02 1.43109079e-02\n", " 1.41521484e-02 1.37856801e-02 1.34740603e-02 1.30587329e-02\n", " 1.30315682e-02 1.24908759e-02 1.15991269e-02 1.11680809e-02\n", " 1.03362895e-02 9.97974912e-03 9.23972538e-03 8.59086723e-03\n", " 8.34768304e-03 2.95582098e-32]\n" ] } ], "source": [ "from sklearn.decomposition import PCA\n", "\n", "pca = PCA(n_components = None)\n", "X_train_pca = pca.fit_transform(X_train)\n", "# X_train_pca\n", "\n", "# Eigenvalues (sum of squares of the distance between the projected data points and the origin along the eigenvector)\n", "print(pca.explained_variance_)\n", "\n", "# Explained variance ratio (i.e. how much of the change in the variables is explained by change in the respective principal component): eigenvalue/(n variables)\n", "print(pca.explained_variance_ratio_)\n" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PCA Dim: 28; PCA Var: [14.99 26.97 33.57 38.45 43.08 46.8 50.29 53.58 56.57 59.22 61.73 64.09\n", " 66.3 68.39 70.33 72.22 74.07 75.78 77.49 79.13 80.67 82.2 83.66 85.09\n", " 86.51 87.89 89.24 90.55]\n" ] } ], "source": [ "# Calculating the \"explained variance\" of each principal component, up to 90% of the total variance\n", "total = sum(pca.explained_variance_)\n", "k = 0\n", "current_variance = 0\n", "while current_variance/total < 0.90:\n", " current_variance += pca.explained_variance_[k]\n", " k = k + 1\n", "k\n", "\n", "# Applying PCA to select the top k components capturing up to 90% of the variance\n", "pca = PCA(n_components = k)\n", "X_train_pca = pca.fit_transform(X_train)\n", "X_test_pca = pca.transform(X_test)\n", "\n", "var_exp = pca.explained_variance_ratio_.cumsum()\n", "var_exp = var_exp*100\n", "\n", "var1=np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)\n", "print(f'PCA Dim: {k}; PCA Var: {var1}')" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'cumulative explained variance')" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deVxU9f748dewgxsKDLjgiooLuCtg6nUJryhu4ZJeyTJbvf6itKysXMprtlh5y3ure+1bWZqlllZIZmoJprhA4oIrKusAsm+znN8f2txIcQBhhhnez8fDx8MzM+fM+z1H5z3n8zmfz0elKIqCEEKIRs/O0gEIIYRoGKQgCCGEAKQgCCGEuEEKghBCCEAKghBCiBukIAghhACkIAghhLjBwdIB3Ilr14oxGGo+jMLDoyk5OUX1EFHDYes52np+YPs5Sn7mZ2enomXLJlU+b9UFwWBQalUQft/X1tl6jraeH9h+jpJfwyJNRkIIIQApCEIIIW6QgiCEEAKQgiCEEOKGei0IO3bsICwsjNDQUDZu3HjT8/v27SM8PJzw8HCeeuopiouL6zMcIYQQt1FvBSEzM5O1a9fy2WefsX37djZv3sy5c+eMzxcUFLBkyRLWrl3Ljh078Pf3Z+3atfUVjhBC2IT6XLGg3gpCbGwsQUFBuLu74+bmxtixY4mOjjY+f+nSJdq0aYOfnx8AI0eOZPfu3fUVjhBCWK384gp+Tkhj3VeJPL52P3uPpdbL+9TbOISsrCy8vLyM22q1msTERON2x44dycjI4PTp0/j7+/P999+TnZ1dX+EIIYTVUBSF1OxiEs5lc/xsNhfSClCAVs2dCentQx8/z3p533orCAaDAZVKZdxWFKXSdvPmzXn11Vd54YUXMBgMTJ8+HUdHxxq9h4dH01rH5+XVrNb7Wgtbz9HW8wPbz1Hy+x+d3kDS+RwOnczg16QMMnNLAPDzdWfWX/0Z0suHjq2bV/oerWv1VhB8fHyIj483bms0GtRqtXFbr9fj4+PDli1bAEhMTMTX17dG75GTU1SrkYBeXs3QaAprvJ81sfUcbT0/sP0cJT8oKdPx24Ucjp3V8NuFXErLdTg62NGjQ0vGDvalTxdPWjZzNr4+O/vOpsKws1Pd9od0vRWEkJAQ1q1bR25uLq6ursTExLBy5Urj8yqVigceeIAtW7agVqv56KOPCAsLq69whBCiQcgrKufY2WyOJWs4lXINvUGhuZsjA7p70c/Pk54dW+HsZG+R2OqtIHh7exMVFUVkZCRarZaIiAgCAwOZP38+CxcuJCAggBUrVvDggw9SUVFBcHAw8+bNq69whBDCYtJzio1F4HxaAQBqd1fuHuhLv26edGnTAju7+msKqi6VUp/3MNUzaTKqmq3naOv5ge3naMv5GRSFvFIdew5d5thZDek51/sDOvo0o183L/p39aSNZ5N67Q+4FYs1GQkhRGOiNxhIvpzHkWQNR5M15BVVYG+nont7d0b1b0e/rp60au5i6TBvSwqCEELUklanJ+nSNY6e0XDsrIbiMh1ODnYEdPZgxEBfOqmb0MSlZndPWpIUBCGEqIHS8ut3Bh05oyHxQg7lFXpcnR3o6+dB/25qenduhbOjvVU2iUlBEEIIE4pKtRw/m82RM1kkXbqGTm+guZsjQT29GdDNC/8OLXGwt/65QqUgCCHELRSWVHDsbDbxp7OMt4d6NHdhVP+29O/mhV/bhnFnUF2SgiCEEDfkF1dwNFlD/OkszlzOw6AoeLm7EDrYl4Hd1XT0aWb2O4PMSQqCEKJRu1ZYbiwCyVfzUBTwbuVGWHB7BnZX46tuatNF4I+kIAghGp1rheXEn84i/kwW567mowBtPJsQHtKRgd3VtPUy/xiBhkAKghCiUcgtKCP+zPUrgXOp+QC082rCpGGdGNBdTVvPJhaO0PKkIAghbFZuQRnxp7M4fCaL86nXp4xo59WUKcM6MdBfTWsPKQJ/JAVBCGFTblUEfNVNmTK8M4P81fi0crNwhA2XFAQhhNX7vU/g0OlMYxFor27K1OGdGShFoNqkIAghrFJ+UTnxZzQcPpXJ2Rsdw743isAgfzXeUgRqrFoFoaysjJSUFLp160ZZWRmurq71HZcQQtykoKSCIzeKwJnLeShA2xsdw4OkT+COmSwIx48fZ8GCBTg4OLBp0yYmTZrE+vXr6d+/vzniE0I0ckWlWo4mazh0KpPTKdcHi7X2cCN8aEcG+atp61X7pXRFZSYLwpo1a/joo49YtGgRPj4+rFmzhldeeYWvvvrKHPEJIRqhCq2e4+eyOZiUyW8XctAbFNQtXQkLbs9gf+9GO06gvpksCGVlZfj5+Rm3R4wYwdq1a+s1KCFE46M3GDiVco1fkzI5kqyhrEKPe1MnxgxsR1BPH9p7N54Rw5ZisiA4ODiQn59vPBEXLlyo96CEEI2DoihcyigkLimDQ6eyKCiuwNXZgYH+aoJ7etO9fUubm0CuITNZEB599FH+9re/kZ2dzZNPPsmBAwdYsWJFtQ6+Y8cO1q9fj06n47777mP27NmVnk9KSuLFF19Eq9XSunVrXnvtNZo3b167TIQQViMzt4Qfjqay5/BlMq+V4mCvok8XT4J6eRPYxQNHB8ssMt/YVWtN5ZSUFA4cOIDBYCA4OJguXbqYPHBmZib33nsvW7duxcnJiZkzZ/Lmm29Wan6aNWsWDz/8MCNGjGD16tU4OzsTFRVV7eBlTeWq2XqOtp4f2F6O+cUVHDqVycGkTC6mF6BSQXdfd4J6+TCwuxduVrSyWHU0xPN3x2sqZ2RksGHDBpYtW8aFCxd4/fXXWb58OV5eXrfdLzY2lqCgINzd3QEYO3Ys0dHRLFiwwPgag8FAcXExAKWlpbRo0aJaSQkhrENZhY5jydnEnczg5MVrGBSF9uqmTB/px7i7OqNodZYOUfyByYKwZMkSRo0aBUDbtm0ZPHgwzz33HB988MFt98vKyqpUNNRqNYmJiTcd+4EHHmDVqlW4urryxRdf1Cj421U6U7y8mtV6X2th6znaen5gnTnq9AaOJ2v46cgVfk3KoLxCj7qlK/eM8mNE/3Z08Gk8zcLWdv5MFoRr164RGRkJgLOzM3PnzmX79u0mD2wwGCrdEaAoSqXtsrIynn/+eT766CMCAwPZsGEDzzzzDO+//361g5cmo6rZeo62nh9YV46KonAhrcDYOVxUqqWJiwPBvXwI6umNX7sW2N34//97TtaUX200xPzuuMlIr9eTmZmJt7c3ANnZ2VSj2wEfHx/i4+ON2xqNBrVabdxOTk7G2dmZwMBAAGbMmMHbb79t8rhCiIYjv6ic2KQMfklMJz2nBEcHO/r6Xe8cDujsYRPrDDcmJgvC3LlzmTx5MsOGDUOlUhEbG8vTTz9t8sAhISGsW7eO3NxcXF1diYmJYeXKlcbnO3ToQEZGBhcuXKBz5878+OOPBAQE3Fk2Qoh6p9MbSDyfwy+J6SSez8GgKPi1a8Hccf4M8lfj6ixTpFkrk2cuIiKC3r17c/DgQezt7Zk3bx7dunUzeWBvb2+ioqKIjIxEq9USERFBYGAg8+fPZ+HChQQEBPCPf/yDJ554AkVR8PDwYNWqVXWSlBCi7qVqivg5MZ2DSRkUlGhp0dSJvw5pz9AAH5lDyEZU67bT1NRU8vPzKzUV9erVq14Dqw7pQ6iaredo6/lBw8ixpEzHoVOZ/JyYzsX0AuztVPT18+SuwNb07twKe7vaNwk1hPzqU0PM7477EN5++23++9//4uHhYXxMpVLx448/1k2EQogGRVEUzl7NZ9/xNOLPZKHVGWjr1YSZo7sS1Mub5m5Olg5R1BOTBeHrr78mJibG2KkshLBNhSUVxJ7IYH9CGuk5Jbg42TM0oDXDAlvT0aeZzCPUCJgsCK1bt5ZiIISNMigKZ1KusS8hjaPJGnR6hS5tm3N/mD+D/b1xdpIpJBoTkwUhODiYNWvWMHr0aFxcXIyPN4Q+BCFE7eQXV3Dgt3T2J6SRda2UJi4O/KVvW4b3bUM7WV+g0TJZELZu3QpAdHS08THpQxDC+hgUhZMXc9mXkMbxs9noDQrdfN2ZNLQTA7p74eQoVwONncmCsGfPHnPEIYSoJ9cKy/nlt3R+TkgjO7+Mpq6OjBnYjuF92sjtoqISkwUhNzeXb775huLiYhRFwWAwkJKSwhtvvGGO+IQQtWAwKJy4mMu+46kknLs+eKxHh5ZE/KUL/bp64eggI4jFzUwWhCeeeAIXFxfOnTtHSEgIsbGxDBgwwByxCSFqKLegjF8S0/k5MY2cgnKauzkydogvw/u0wbulm6XDEw2cyYKQlpbG7t27WbZsGTNnzuTvf/87jz32mDliE0JUg95g4LcLuew/nkbC+WwUBXp1bMmMUV3p29VT5hMS1WayIHh6egLQsWNHkpOTmThxIjqdzGEuhKUVllSwPyGNn46lkltQTosmToQFdWBYnzao3V0tHZ6wQiYLgoeHBx9++CF9+/Zl3bp1NG3alLKyMnPEJoS4hcuZhew+cpWDSZno9AZ6dGjJvaO70sdPrgbEnTFZEFasWMG3337LwIED6d27N++88w6LFi0yR2xCiBv0BgPHkrPZHX+F5Kv5ODnacVdga0b3b0tbGTcg6ki1JrdrqGRyu6rZeo62nh9cz/FCSk6lZiHPFi6M6t+OYX1a08TK1yC29XPYEPOr9eR29957L59//jn9+vW75RwmR48erZsIhRA3uZxZyOd7zvHTkavGZqHZY7rRx88TOzuZU0jUjyoLwu+rl23YsKHS2shCiPphMCgknMsm5vAVzlzJw8nRnrsCfBg1oJ1MJyHMosqC8Ptyl0uWLKk0bYUQom6Vluv45bd0foy/SlZeKa2aOzNtZBemjupGaXG5pcMTjYjJTuW2bdty9OhR+vbti90dLIYhhKgsO6+U3Ueu8nNiGqXlerq0bc49f+lC/26e2NvZ0dTNSQqCMCuTBeH8+fPMmjULBwcHnJycUBQFlUpVrT6EHTt2sH79enQ6Hffddx+zZ882Pnfq1CmWLFli3M7NzaVFixbs3LmzlqkI0fApisK51HxiDl/haLIGFSoG+ntx9yBfurRpYenwRCNnsiBs3LixVgfOzMxk7dq1bN26FScnJ2bOnMmQIUPw8/MDoEePHnz99dcAlJaWMm3aNJYtW1ar9xKiodPpDcSfziLm8BUuZRTSxMWBvw5pz+j+7WjV3MX0AYQwg2o1GZ08eZKSkhIURUGv13P58mWmT59+2/1iY2MJCgrC3d0dgLFjxxIdHc2CBQtueu2///1vBg0axMCBA2uZhhANU0mZjn0JqeyOv8q1wnJ8WrkxJ7QbIb1by+IzosExWRCWLl3Kjz/+SHl5OWq1msuXLzNgwACTBSErK6vS3UlqtZrExMSbXldYWMgXX3zBjh07ahz87e6nNcXLq1mt97UWtp5jQ84vK7eEb36+QMyvKZSW6wj08+Tv0/sywN+7RreNNuQc64Lk17CYLAixsbH8+OOPLF++nMcff5z09HQ+/PBDkwc2GAyVxi/83vfwZ9988w1jxozBw8OjhqHLwLTbsfUcG2p+F9ML2HXoMvGnNQAM7qlm7KD2dPC5/sWQk1NU7WM11BzriuRnfrUemPY7Ly8v3Nzc6Ny5M8nJyYwZM4aXX37Z5Bv7+PgQHx9v3NZoNMZbWf9o9+7dPPzwwyaPJ0RDZVAUEs/lsOvQZc5cycPV2Z7QQb6MGSj9A8K6mCwIjo6OHD58mC5durB//36GDBlCSUmJyQOHhISwbt06cnNzcXV1JSYmhpUrV1Z6jaIoJCUl0a9fv9pnIISFVGj1xCZlEHPoChm5JbRq7syMUX4M79MGV2eT/7WEaHBM/qtdtGgRn3zyCatXr+b9998nKCiIhx56yOSBvb29iYqKIjIyEq1WS0REBIGBgcyfP5+FCxcSEBBAbm4ujo6OODs710kyQphDabmOvcdT2XXoCgXFFXTwbsZDE3sysLtaZhsVVs3k5HYXLlygc+fOxu3CwkKaNWsYHSXSh1A1W8/REvmVlGnZfeQqPxy+QnGZjl4dWxIW3BH/9u637B+7U3IOrVtDzO+O+xDmzp2Lr68v06ZNY9y4cQ2mGAhhLoUlFfwQf4Ufj1yltFxPXz9Pxod0kIFkwuaYLAh79+7l559/Ztu2bbz++uuEhoYyffp0/P39zRGfEBaTV1TOrkOX+elYKlqtgQHdvZgQ0pH23vKjSNgmkwXBzs6OESNGMGLECM6fP8+zzz7L559/zqlTp8wRnxBml5Nfxve/prA/IR29wUBQT2/CgjvS1rOJpUMTol6ZLAg6nY49e/awdetWEhMTCQsLu+luISFsQXZeKTvjLnHgtwwAhgb4MC6oA94t3SwbmBBmYrIg3HXXXXTt2pWIiAjeeecdnJyczBGXEGaTW1DGzrgUfk5IQ6WCEX3bMG5IBzxayBgC0biYLAibNm2iY8eOZghFCPPKLyrn27gU9h5PQ1EUhvVpw4TgDjKYTDRaJguCFANhawpLKvj+18vsOXIVnV4hJMCHiSEd8XR3tXRoQliUDKcUjUZxmZZdhy7zQ/xVKir0BPXyZuLQTni3kj4CIUAKgmgESst1/HD4CrsOX6G0XMcgfzWT7upEG7lrSIhKqiwIhw8fvu2OgwYNqvNghKhLWp2BPUevsjP2EsVlOvp19WTysM74qmXBeiFupcqCsGLFCuD6amZpaWn4+fnh4OBAcnIyXbp0Ma52JkRDY1AUDp3MZOv+C2Tnl9GrUyvuGdGZjj7NLR2aEA1alQXh9wVrnnjiCdasWUP//v0BSEpK4l//+pd5ohOihk6lXOOLn86RklFIe3VTnprRl16dWlk6LCGsgsk+hIsXLxqLAUCvXr1ISUmp16CEqKmrmiK+3HuexPM5eDR3Zv6Engzp5Y1dPUw6J4StMlkQXFxc2Lp1K5MmTUJRFLZs2ULz5nLpLRqGa4XlbPv5Agd+S8fVyYHpI/0YPaAtjg6yXrEQNWWyIKxatYpFixaxdOlSVCoVvXr14o033jBHbEJUqaRMy1f7zvPD4SsYFIW7B/oyIaQjTV0dLR2aEFbLZEHo0qUL27ZtIy8vDwB3d/d6D0qIqugNBvYfT+PrA5coKK4gqKc3U4Z3xksGlQlxx0wu76TRaHjooYeYMWMGer2eefPmkZWVZY7YhKjkdMo1lm+I55OYZDr4NOfFuQN5aGIvKQZC1BGTBWH58uWMGTMGZ2dnmjdvjr+/P0uXLq3WwXfs2EFYWBihoaFs3LjxpucvXLjAnDlzmDhxIvPmzSM/P7/mGQibl51fynvbfmPN58coLdfx2OTevPJoiNxGKkQdM1kQUlNTmT59OnZ2djg6OrJ48WLS09NNHjgzM5O1a9fy2WefsX37djZv3sy5c+eMzyuKwqOPPsr8+fP55ptv6NGjB++///6dZSNsSrlWz/afL/D8B7+SeD6HycM68cr8IQz0V9fLkpVCNHYm+xBUKhUGg8G4XVRUVGm7KrGxsQQFBRn7HMaOHUt0dDQLFiwAro9ncHNzY/jw4QA88sgjFBQU1CoJYVsUReHw6Sy++OkcuQXlDO6hZtpf/GQ6aiHqmcmCEBoayqJFiygsLGTTpk1s2bKFcePGmTxwVlYWXl5exm21Wk1iYqJx+/Lly3h6evLcc89x6tQpOnfuzAsvvFDLNIStuJxZyGe7z5J8JQ9fdVPmT+hJ9/YtLR2WEI2CyYLwyCOPsH37dgwGA7GxscyYMYNp06aZPLDBYKh0Wa8oSqVtnU7HoUOH+PTTTwkICOCtt95i9erVrF69utrBe3jUfk4aLy/bXxfXmnLMLypnY/Rpdh28RBNXJx6L6EPokA7Y21XdNGRN+dWWreco+TUs1ZrtdPLkyUyePLlGB/bx8SE+Pt64rdFoUKvVxm0vLy86dOhAQEAAABMmTGDhwoU1eo+cnCIMBqVG+1x/72ZoNIU13s+aWEuOBkXhl8R0tvx0jtJyPaP6t2PSsE40cXEkN6eoyv2sJb87Yes5Sn7mZ2enuu0PaZMFYffu3axatYr8/HwU5X9fvkePHr3tfiEhIaxbt47c3FxcXV2JiYmptBZzv379yM3N5fTp0/j7+7Nnzx569epVnZyEjbiSVcTHu05zPrWAbr7u/C20G+28ZCZSISzFZEF47bXXWLJkCT179qzRnR3e3t5ERUURGRmJVqslIiKCwMBA5s+fz8KFCwkICODdd99l6dKllJaW4uPjw5o1a+4oGWEdyip0fP3LRX44fBU3Fwfmje9BSG8fuXNICAtTKX/82X8L06ZNY8uWLeaKp0akyahqDTFHRVE4mqzhs91nuVZYzvA+bYj4S5daTTfREPOra7aeo+RnfnfcZNSnTx/27dvHiBEj6jQw0bho8krZ+EMyiedzaOfVlEcn98avbQtLhyWE+AOTBWHfvn18+umnODo64ujoaLxbyFQfghAAOr2B6F8vszP2Eio7FTNH+TF6YDvs7UyOiRRCmJnJgvDRRx+ZIQxhi85cvsbHu86QnlPCgO5e3Du6K62ay+AyIRqqKgtCXFwcwcHBJCUl3fL5tm3b1ltQwrqVa/V8+dN5fjx6Fc8WLjwxrQ+BXTwsHZYQwoQqC8K3335LcHAwn3zyyU3PqVQqQkND6zUwYZ0uphfwwY6TZOSWcPdAX+4Z0RknR1msRghrUGVBePnllwFuWRCE+DOd3sC3cSnsOHAJ92ZOLJ7Zlx4dZS1jIayJyT6ES5cu8emnn1JSUoKiKBgMBlJSUti0aZM54hNWID2nmA93nuRieiHBvbyZfXc33Fxk5TIhrI3JWz2eeuoptFotx44do23btpw7d45u3bqZIzbRwCmKwo9HrrJ8w2GyrpXy6OTezA/vJcVACCtlsiAUFxezfPly7rrrLoYPH86GDRs4fvy4OWITDdi1wnLe/CKBjT8k0629OyvmDWGQv9r0jkKIBstkk9Hv6xl06NCBs2fPEhgYKFMMNHKHTmXyya4zaPUG5oztzl/6tpF/E0LYAJMFoUOHDrzyyitMmTKF559/npKSEnQ6nTliEw1MSZmWT2OSOXgyk06tmzM/vCc+rdwsHZYQoo6YLAjLli1j//799OzZk2nTpnHgwAFWrFhhjthEA5KSUch7238jJ7+cycM6MT64g4w2FsLGVFkQ8vLyjH8fMmQIeXl5hIWFERYWZpbARMOgKAr7EtL47IezNHNzZMns/vi1kzmIhLBFVRaEoKAgVCoVt5oMVaVScerUqXoNTFheeYWej3edIS4pg16dWjE/vCfN3ZwsHZYQop5UWRBOnz5tzjhEA5OeU8x7206Qll3MpLs6ER7SEbvbLGcphLB+JvsQ9Ho9mzZt4pdffsHe3p5Ro0YxdepUc8QmLOTXk5l89P1pHB3seHJGX3p1khHHQjQGJgvCypUrOX/+PJMmTUJRFL788ktSUlKIiooyR3zCjLQ6A5v3nGXP0VT82rbgkUm9ZHZSIRoRkwUhNjaWb7/9FkfH66NPJ06cyMSJE6tVEHbs2MH69evR6XTcd999zJ49u9Lz//znP/nqq69o3rw5ANOnT7/pNcI8svNKWf/1CS6mFzJ2sC/3jOiCg73cRSREY2KyILRq1Qq9Xm8sCCqVyvgFfjuZmZmsXbuWrVu34uTkxMyZMxkyZAh+fn7G15w4cYI333yTfv363UEK4k4lnMvmw50nMSgKj0/pzYDuMuJYiMbIZEHw9/dn1qxZTJ06FXt7e7777jtatmzJhg0bALj//vtvuV9sbCxBQUHGkc5jx44lOjqaBQsWGF9z4sQJ/v3vf5OamsqgQYN45plncHZ2rou8RDUYFIXtP19gZ2wK7dVNeXRKb7xbykAzIRork20C5eXldO/enaSkJBITE2nXrh0eHh4kJyeTnJxc5X5ZWVl4eXkZt9VqNZmZmcbt4uJievToweLFi9m2bRsFBQW89957d5iOqK6SMh3rvkxkZ2wKdwW25rk5A6QYCNHImbxCWLx4Ma1aVb7L5PTp0/j7+992P4PBUGl+m9/XYv5dkyZN+OCDD4zbDzzwAM8991yNOqs9PJpW+7V/5uXVrNb7WouqckzVFPGPjUfJyCnmkSkBhA3tZJVzETXmc2grJL+GxWRBmDp1Km+88QYDBgwA4OOPP+a9997j4MGDt93Px8eH+Ph447ZGo0Gt/l/bdFpaGrGxsURERADXC4aDg8lwKsnJKcJguHngnCleXs3QaAprvJ81qSrHxPM5/PubJOztVDw1oy/+HVqSnV1kgQjvTGM+h7ZC8jM/OzvVbX9Im/wGXrVqFU8++SQzZ84kISGBwsJCtmzZYvKNQ0JCWLduHbm5ubi6uhITE8PKlSuNz7u4uPDaa68xZMgQ2rVrx8aNG7n77rurmZaoKUVR+P7Xy3y19zy+6qYsuCcAzxaulg5LCNGAmCwIISEhvPjiiyxYsABPT0+++uqrSr/0q+Lt7U1UVBSRkZFotVoiIiIIDAxk/vz5LFy4kICAAFasWMGjjz6KVqulf//+VXZQiztTrtWz4btTHDqVxeAeau4P64GzrHMshPgTlXKryYr+4LXXXuPrr79m5cqVnD17lk8++YQXX3yxQfyalyajqv2eY05+Geu2JnIls4ipIzoTFtTBKvsL/qwxnUNbJfmZ3x03GSUlJbFt2za8vLwYOXIkQUFBPPXUUw2iIIjbO3P5Gu9tP4FOr/D/pgUS2MXT0iEJIRowkwVhw4YNqFQqCgoKaN68OYGBgWzfvt0csYlaUhSF72Iv8v623/Byd+Xv9wTQ2qOJpcMSQjRwJschXLp0ibCwMMaPH09mZibjxo0jIyPDHLGJWtDpDXy86wzrv0qkV6dWLI0cKMVACFEtJgvCypUref755/Hw8MDb25u//e1vvPjii+aITdRQSZmOt7cksO94GtNGd2XhPYG4udTsVl4hRONlsiDk5eUxdOhQ4/bs2bMpKrK++9ZtXU5+Gf/YeITTl/N4IKwHkWE9Zf0CIUSNVOvnY3l5ufHOFI1Gg8FgqNegRM1cyijg7S2JVOgMRE3vQ8+Osn6BEKLmTBaEWbNmMW/ePHJycnjjjTf49ttvefDBB80Rm6iG42ez+VP0yw0AABrySURBVNc3J2jm6sSie/vR1lP6C4QQtWOyIERERNChQwf27t2LTqdj5cqVlZqQhOXsjr/C5z+epYN3M/5fRCAtmspMsUKI2qtWk9GgQYMYNGhQfcciqslgUNi05yy746/Sr6snD4X3wtlJRh4LIe6M3IJiZcor9Ly/I4ljZ7O5e6AvM0b5SeexEKJOSEGwIvlF5bz9ZSIpmYXMvrsbowe0s3RIQggbIgXBSqRqinhrSyKFpRX8fWogfbvKNBRCiLplchyCRqPhoYceYuzYsWRnZzNv3jyysrLMEZu44WJ6Aas+PYpOb2DJ7P5SDIQQ9cJkQVi+fDljxozB2dmZFi1a4O/vz9KlS80RmwBSs4tZ+0UCTVwceD5yAB19mls6JCGEjTJZEFJTU5k+fTp2dnY4OjqyePFi0tPTzRFbo5edV8qbm49jb6di0cy+sqCNEKJemSwIKpWq0sjkoqIiGalsBvnFFby++TjlFXqenNEXdUs3S4ckhLBxJjuVQ0NDWbRoEYWFhWzatIktW7Ywbtw4c8TWaJWUaXlz83HyispZNLMfvuqqF7QQQoi6YvIK4ZFHHmH48OEEBAQQGxvLjBkzePzxx6t18B07dhAWFkZoaCgbN26s8nV79+5l1KhR1Y/ahpVr9bz1ZSJp2cUsmBqAX9sWlg5JCNFImLxC2LRpExMmTGDy5Mk1OnBmZiZr165l69atODk5MXPmTIYMGYKfn1+l12VnZ/Pqq6/WLGobpdMbeG/bCc6n5vPopN707uRh6ZCEEI2IySuEX3/9lTFjxvDcc89x/Pjxah84NjaWoKAg3N3dcXNzY+zYsURHR9/0uqVLl7JgwYKaRW2DDAaFD3ee5LcLOdz3V38G+qstHZIQopExWRDWrl3Lrl276NWrFy+//DITJkzg//7v/0weOCsrCy8vL+O2Wq0mMzOz0ms+/vhjevbsSZ8+fWoRuu1QFIVPf0jm0Kkspo3swvA+bSwdkhCiEarWSOUWLVowY8YM1Go1H3zwAR988AH33XffbfcxGAzGNRTg+pfeH7eTk5OJiYnho48+qvWSnB4ete9s9fJqVut969rH351k77FUIkZ1JXJ8zzo7bkPKsT7Yen5g+zlKfg2LyYJw8uRJvvrqK6Kjo+nZsycPPvhgtTqAfXx8iI+PN25rNBrU6v81g0RHR6PRaLjnnnvQarVkZWUxa9YsPvvss2oHn5NThMGgVPv1v/PyaoZGU1jj/erD97+msOWn84zo24Zxg9rVWVwNKcf6YOv5ge3nKPmZn52d6rY/pE0WhMcee4x77rmHLVu20KZN9ZsyQkJCWLduHbm5ubi6uhITE8PKlSuNzy9cuJCFCxcCcPXqVSIjI2tUDGzB/oQ0tvx0nkH+auaEdq90BSWEEOZmsiD89NNPtfqi8vb2JioqisjISLRaLREREQQGBjJ//nwWLlxIQEBArQK2FUfOZPF/0afp3bkV88Nl/WMhhOWpFEW5ZZvLvffey+eff06/fv1u2Rdw9OhRswVZFWttMjp7NY/XNx2nvbopi2b2q5fFbSydY32z9fzA9nOU/Myv1k1Gb7/9NgA7d+686bkqaoiohvScYt75MpFWzZxZGBEoK50JIRqMKm87/b0D+KWXXqJt27aV/jz55JNmC9CW5BdXsPaLBOzsVERN70MzNydLhySEEEZVXiEsXLiQixcvcuXKFcLDw42P63Q6nJzki6ymyip0vLUlgYKSCp6Z1V8mqxNCNDhVFoSnn36a1NRUXnjhBV544QXj4/b29jdNPyFuT28w8K+vk7icWcjf7wmkU2tZ00AI0fBUWRDatWtHu3btiI6Oxs6ucstSSUlJvQdmKxRF4ZNdySSezyFybHf6+slqZ0KIhsnkbad79uzhnXfeoaSkBEVRMBgM5OXlcezYMXPEZ/V2xqWwPyGN8cEd+Eu/tpYORwghqmSyIKxZs4YnnniCzz//nPnz57N7926aNGlijtis3oHf0tm2/wLBvXyYOryzpcMRQojbMjm5naurK2FhYfTt2xdnZ2eWLVvG3r17zRCadUu6mMtH35+mR4eW3B/mL6OQhRANnsmC4OzsTEVFBe3bt+fUqVPY2dnJl5sJlzMLeXfbb7T2cOPxKQE42Jv8mIUQwuJMNhmNGjWKhx56iFdffZUZM2Zw5MgRWrZsaY7YrFJOfhlvbUnA1dmBJ6b1wc2lWhPKCiGExZn8tnrkkUeYOHEi3t7evPfeexw+fJgJEyaYIzarU1Km5a0tCZRrDTz7t/60au5i6ZCEEKLaqiwIMTExlbZPnDgBQOvWrTly5AihoaH1G5mVURSFf39zkozcEp6c0Zd2XrVfq0EIISyhyoLwySefVLmTSqWSgvAnv/yWzm8Xcph9dzd6dJAmNSGE9alVQRCV5RWVs/nHc3TzdWdkfxlrIISwTib7EF5++eVbPr506dI6D8YaXR+JfAat3sD94/yxkzuwhBBWyuT9kO7u7sY/TZo04dChQ+aIy2ocPp3FsbPZTB7WCe9WMmGdEMJ6mbxCWLBgQaXt+fPn8+ijj9ZbQNaksKSCjT8k09GnGaGDfC0djhBC3JEaj5hq2rQpWVlZ1Xrtjh07CAsLIzQ0lI0bN970/A8//EB4eDjjx49nyZIlVFRU1DQci/r8x7OUlOl4IKwH9nYy+EwIYd1q1IegKApJSUl07mx6Xp7MzEzWrl3L1q1bcXJyYubMmQwZMsQ4dXZJSQkrVqxg27ZteHp6EhUVxbZt25gxY8YdpGM+CeeyOZiUycShHWmnlltMhRDWz2RBcHd3r7Q9ceJEJk6caPLAsbGxBAUFGfcfO3Ys0dHRxiYoNzc39uzZg6OjI6WlpeTk5NC8uXWsE1BSpuPjXWdo69WECSEdLR2OEELUiRr3IVRXVlYWXl5exm21Wk1iYmKl1zg6OrJv3z6efvpp1Go1d911V63ey9y27D1HXlE5C6bKPEVCCNthsiB89913vPPOO+Tn51d6PC4u7rb7GQyGSpPgKYpyy0nxRowYwa+//sqbb77JsmXLeOONN6obOx4etW+q8fJqVqv9Es5q2Hc8jSl/8WNwYMMec1DbHK2FrecHtp+j5NewmCwIr732GkuXLqV9+/Y1OrCPjw/x8fHGbY1Gg1qtNm7n5eVx4sQJ41VBeHg4UVFRNXqPnJwiDAalRvvA9ZOk0RTWeL/yCj1vfX4UdUtXQge0rdUxzKW2OVoLW88PbD9Hyc/87OxUt/0hbbK9o23btowePZquXbtW+mNKSEgIcXFx5ObmUlpaSkxMDMOHDzc+rygKixcvJi0tDYDo6Gj69+9fnZwsZuv+C2Tnl3H/OH+cHe0tHY4QQtQpk1cIkydP5tVXX2X48OE4OPzv5YMGDbrtft7e3kRFRREZGYlWqyUiIoLAwEDmz5/PwoULCQgIYOXKlTz88MOoVCr8/PxYvnz5nWdUT86l5rM7/goj+7Wle3uZq0gIYXtUiqLcts1l8eLF7N+/v1JzD1wfY2Bp5moy0uoMLNtwiHKtnpXzhuDq3PDXOGiIl6t1ydbzA9vPUfIzP1NNRia/2U6ePMn+/ftxdnau08CsyY7Yi6TnlBA1vY9VFAMhhKgNk30Inp6e6HQ6c8TSIF3OLOS7uMuE9PYhoLOHpcMRQoh6Y/Lnrre3N5MmTSIkJAQnJyfj441htlODorDhu9M0dXNk5mjTHelCCGHNTBaE9u3b1/iWU1tx7mo+KZmF3D/On6aujpYORwgh6lW9jVS2BQeTMnBytGNQD7XpFwshhJUzWRDCw8Nv+XhDuMuoPml1Bg6fzqJ/Ny9cnKQjWQhh+0x+073wwgvGv2u1Wr799lt8fW1/7v/fLuRQXKYjuJePpUMRQgizMFkQBg8eXGk7JCSEmTNn2vwiOXFJGTR3c6RnRxmEJoRoHGo8Vee1a9eqvUCOtSop05JwLpvBPb1l4RshRKNR4z6EtLQ0q1nEprbiz2jQ6RVpLhJCNCo16kNQqVS0atWKLl261GtQlhZ3IgOfVm509LGuqWuFEOJOmGwPad++Pd999x2DBw/Gw8ODN954g+zsbHPEZhE5+WWcuZJHcC/vW67fIIQQtspkQViyZIlxDeW2bdsyePBgnn322XoPzFIOnswAYIg0FwkhGhmTBeHatWtERkYC4OzszNy5c9FoNPUemCUoikJcUiZ+7Vqgdne1dDhCCGFWJguCXq8nMzPTuJ2dnY2JGbOt1pWsItKyi6UzWQjRKJnsVJ47dy6TJ09m2LBhqFQqYmNjefrpp80Rm9nFnsjA3k7FIH+ZqkII0fiYLAgRERH07t2bgwcPYm9vz7x58+jWrZs5YjMrg0Hh15OZBHbxkInshBCNUrUm6fH398ff37++Y7GoUynXyC+ukOYiIUSjVa/DcHfs2EFYWBihoaFs3Ljxpud3797NpEmTmDhxIo899hj5+fn1Gc5txSVl4OpsTx8/WQRHCNE41VtByMzMZO3atXz22Wds376dzZs3c+7cOePzRUVFLFu2jPfff59vvvmG7t27s27duvoK57bKK/QcSdYwsLsaRwd7i8QghBCWVm8FITY2lqCgINzd3XFzc2Ps2LFER0cbn9dqtbz00kt4e3sD0L17d9LT0+srnNs6dk5DeYVemouEEI1avU30n5WVhZeXl3FbrVaTmJho3G7ZsiV33303AGVlZbz//vvMmTOnRu/h4dG01vF5ef1vWoqjZ5PwbOHC0P6+2NnZzujkP+Zoi2w9P7D9HCW/hqXeCoLBYKg09YOiKLecCqKwsJDHH38cf39/pkyZUqP3yMkpwmCo+ZgIL69maDSFABQUV3D0dBZjh/iSk1NU42M1VH/M0RbZen5g+zlKfuZnZ6e67Q/pemsy8vHxqTSiWaPRoFZXvr8/KyuLWbNm0b17d1555ZX6CuW2Dp3KxKDIzKZCCFFvBSEkJIS4uDhyc3MpLS0lJiaG4cOHG5/X6/U88sgjjBs3jueff95iE8nFJWXiq25KO6/aNz8JIYQtqLcmI29vb6KiooiMjESr1RIREUFgYCDz589n4cKFZGRkcPLkSfR6Pbt27QKgd+/eZr1SyMgt4WJ6AdNH+pntPYUQoqGq19Xjw8PDb1pg54MPPgAgICCA06dP1+fbm3QwKQMVMKSnt0XjEEKIhqDRrg95fWbTDHp0bEnLZs6WDkcIISyu0RaE82kFaPLKpDNZCCFuaLQFIS4pAycHO/p38zL9YiGEaAQaZUHQ6gwcOplJ366euDrXazeKEEJYjUZZEI6ezqS4TCfNRUII8QeNsiD8dPQqTV0d6dWplaVDEUKIBqPRFYSSMh2HkjIY0sMbB/tGl74QQlSp0X0j5hWVY2+n4q7A1pYORQghGpRG16PaxrMJn60cR961EkuHIoQQDUqju0IAZBEcIYS4hUZZEIQQQtxMCoIQQghACoIQQogbpCAIIYQApCAIIYS4QQqCEEIIwMrHIdjZ1X7ZzTvZ11rYeo62nh/Yfo6Sn3mZikelKIpipliEEEI0YNJkJIQQApCCIIQQ4gYpCEIIIQApCEIIIW6QgiCEEAKQgiCEEOIGKQhCCCEAKQhCCCFukIIghBACaIQFYceOHYSFhREaGsrGjRstHU6dmzNnDuPHj2fSpElMmjSJhIQES4dUJ4qKipgwYQJXr14FIDY2lvDwcEJDQ1m7dq2Fo6sbf87x2WefJTQ01Hguf/jhBwtHWHv//Oc/GT9+POPHj2fNmjWA7Z3DW+VodedQaUQyMjKUkSNHKteuXVOKi4uV8PBw5ezZs5YOq84YDAblrrvuUrRaraVDqVPHjx9XJkyYoPTq1Uu5cuWKUlpaqowYMUK5fPmyotVqlQceeEDZu3evpcO8I3/OUVEUZcKECUpmZqaFI7tzBw4cUGbMmKGUl5crFRUVSmRkpLJjxw6bOoe3yjEmJsbqzmGjukKIjY0lKCgId3d33NzcGDt2LNHR0ZYOq85cuHABgAceeICJEyfy6aefWjiiuvHFF1/w0ksvoVarAUhMTKRDhw74+vri4OBAeHi41Z/HP+dYWlpKWloazz33HOHh4bzzzjsYDAYLR1k7Xl5eLFmyBCcnJxwdHenSpQuXLl2yqXN4qxzT0tKs7hw2qoKQlZWFl5eXcVutVpOZmWnBiOpWQUEBwcHBvPvuu3z00Uds2rSJAwcOWDqsO/bKK68wcOBA47Ytnsc/55idnU1QUBCrVq3iiy++ID4+ni+//NKCEdZe165d6du3LwCXLl3i+++/R6VS2dQ5vFWOw4YNs7pz2KgKgsFgQKX63/SviqJU2rZ2/fr1Y82aNTRr1oxWrVoRERHBvn37LB1WnbP18wjg6+vLu+++i1qtxtXVlTlz5lj9uTx79iwPPPAATz/9NL6+vjZ5Dv+YY+fOna3uHDaqguDj44NGozFuazQa4yW6LYiPjycuLs64rSgKDg5WveTFLdn6eQQ4c+YMu3btMm5b+7k8cuQIc+fO5amnnmLKlCk2eQ7/nKM1nsNGVRBCQkKIi4sjNzeX0tJSYmJiGD58uKXDqjOFhYWsWbOG8vJyioqK2LZtG3fffbelw6pzffr04eLFi6SkpKDX69m5c6dNnUe4/uWxatUq8vPz0Wq1bN682WrPZXp6Oo8//jivv/4648ePB2zvHN4qR2s8hw27XNUxb29voqKiiIyMRKvVEhERQWBgoKXDqjMjR44kISGByZMnYzAYmDVrFv369bN0WHXO2dmZ1atX8/e//53y8nJGjBjBX//6V0uHVaf8/f156KGHuPfee9HpdISGhjJhwgRLh1Ur//nPfygvL2f16tXGx2bOnGlT57CqHK3tHMqKaUIIIYBG1mQkhBCialIQhBBCAFIQhBBC3CAFQQghBCAFQQghxA1SEIRVGzVqFL/99ptZ3quoqIiZM2cyfvx4YmJizPKeDcWWLVtscnZgUVmjGocgxJ04deoUOTk5DX8K43pw5MgRunbtaukwRD2TgiDq1a+//sratWvx9fXl7Nmz6HQ6li9fzoABA1iyZAldu3Zl3rx5AJW2R40axYQJEzh48CD5+fk8+OCDHD16lKSkJBwcHFi/fj3e3t4AfPbZZ5w+fZqKigruv/9+IiIiANizZw/r169Hq9Xi4uLCM888Q79+/Vi3bh3Hjx8nKyuL7t278/rrr1eKeffu3fzzn//EYDDQpEkTnn32WZo2bcpzzz1HZmYmkyZNYvPmzbi4uBj30Wg0vPTSS1y4cAE7OztmzpxJZGQkGRkZLFu2jNTUVBRFYfLkyTz44INcvXqV++67j6FDh3LixAn0ej0LFy5k8+bNXLhwgd69e/Pmm2+SlpbGnDlzGDZsGAkJCSiKwosvvsjAgQPRarWsXr2auLg47O3tCQwMNMY6atQopkyZQlxcHOnp6UyaNIknnnjC5OeSmpqKRqMhNTUVb29vXnvtNRISEtizZw8HDhzAxcWFoKAgnn/+eSoqKlAUhYiICGbPnm2Of06ivllk0m3RaBw8eFDp0aOHcvLkSUVRFOU///mPMnv2bEVRFOWZZ55RPvzwQ+Nr/7g9cuRIZdWqVYqiKMq3336r+Pv7K6dOnVIURVEee+wxZf369cbXvfTSS4qiXF/vIjg4WElOTlYuXryoTJgwQcnNzVUURVGSk5OVoUOHKsXFxco777yjjB079pbrRpw7d04JCQlRLl++rCiKosTGxipDhw5VCgsLlYMHDyrjx4+/ZZ6PP/648uqrryqKoigFBQXK+PHjlUuXLimzZ89W/vvf/xofDw8PV3bu3KlcuXJF6datm7J7925FURTlxRdfVEaOHKkUFhYqZWVlytChQ5UjR44YX/fNN98oiqIoe/fuVYYOHapUVFQob7/9trJgwQKloqJC0ev1ypIlS5QXXnjB+LmsXr3a+LkEBAQoly9fNvm5jB49WiksLFQURVEefvhh5e23377p3Dz77LPKv//9b0VRFCUrK0t54oknFL1ef7t/BsJKyBWCqHdt2rShR48eAPTs2ZNt27ZVa7/Q0FDg+syfnp6e+Pv7A9C+fXvy8/ONr5s5cyZwfWqSoUOHGn8xZ2VlMXfuXOPrVCoVly9fBqBv3763nGjs4MGDBAUF4evrC0BwcDCtWrXixIkTt52NMzY2lsWLFwPQrFkzdu7cSUlJCUePHuW///2v8fGpU6eyf/9++vTpg6OjI6NGjTLm1K9fP5o2bQpcnw46Pz8ftVpNixYtCA8PB2DEiBHY29tz5swZ9u/fT1RUFI6OjsD11fIef/xxY0yjR482fi4eHh7k5+eTkJBw289l8ODBxhh69uxZ6XP+3d13380zzzxDYmIiwcHBLF26FDs76Y60BVIQRL37Y9OKSqVCuTFbyh//DqDVaivt5+TkZPz77196t/LHLyODwYCDgwN6vZ7g4GDeeust43Pp6emo1Wp++OEH3NzcbnmsP0+tDdcnKdPpdLeNwcHBodJ+V65cwd3dvVJ+vx9fp9MZc/rjPlUd397e/qZj2Nvb3xSrwWCo9Bk6Ozsb//77Z20wGG77uVR1rv5o5MiR7Nq1i9jYWOLi4nj33XfZunUrPj4+t/5whNWQsi4spmXLlpw4cQKAzMxMDh06VKvj/H7FkZaWRlxcHMHBwQQHB3PgwAHOnz8PwL59+5g4cSJlZWW3PVZwcDC//PILV65cATC2wffp08fkfl999RVwfdbZ++67j5SUFPr06WO8O6ewsJDt27cTEhJSo/xyc3PZv38/cL3939HRkW7dujFs2DA+//xztFotBoOBjRs3MnToUJNx1uZzsbe3Nxayp556iu+++47x48fz0ksv0bRpU+MVhrBucoUgLGbOnDksWrSIsWPH0q5dO4KCgmp1nPLycqZMmYJWq2Xp0qV06tQJgBUrVvDkk08a56Ffv349TZo0ue2x/Pz8eOmll1iwYAF6vR4XFxf+9a9/0axZs9vu9+KLL7Js2TLCw8NRFIWHH36Y3r178/rrr7NixQq2bt1KRUUF4eHhTJ06ldTU1Grn5+zszNdff83rr7+Oi4sL7777Lvb29jz66KO8+uqrTJ48GZ1OR2BgIC+88ILJ/GrzuQwfPtw4k+djjz3G888/z+bNm7G3t2fMmDEMGjSo2vmIhktmOxWiAbt69Srh4eEcO3bM0qGIRkCajIQQQgByhSCEEOIGuUIQQggBSEEQQghxgxQEIYQQgBQEIYQQN0hBEEIIAUhBEEIIccP/B9w10EUKO+uLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(np.cumsum(pca.explained_variance_ratio_))\n", "plt.xlabel('number of components')\n", "plt.ylabel('cumulative explained variance')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q4. In your opinion, which of the following is the best explanation of the plot of the separation of the two classes when the data are projected only onto the first two principal components ?**\n", "\n", "1. There is good class separation due to the large amount of explained variance from the first two principal components\n", "2. There is only limited class separation although there is a large amount of explained variance from the first two principal components\n", "3. There is good class separation although there is only a limited amount of explained variance from the first two principal components\n", "4. There is only limited class separation due to the limited amount of explained variance from the first two principal components\n", "5. None of the above" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q4 answer.** _Your answer goes here._ " ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEXCAYAAABF40RQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXwT1fo/8E+WJmmalm7pQqFlpywFZG/ZBKGAUJCCiqwuLOIFBGUTEC7IppfvBawKlquXn4pXEVEU2RQEoQVkuUBlXwstXdItbdLsc35/9HYgTdMtbZO2z/v14qVzJjPzdJI8OXNmzjkCxhgDIYSQBkXo7AAIIYTUPkr+hBDSAFHyJ4SQBoiSPyGENECU/AkhpAGi5E8IIQ1Qucn/4sWLmDx5MmJiYjBy5EhMmzYNt27dqo3YynT58mWsWLECAJCUlIS5c+fW+DEtFgtmzZqFoUOH4quvvqrQNkuWLMFnn31Ww5FVzT/+8Q+cPHnSoX2sXr0acXFx1RRR6SZPnoyDBw/alKekpOCpp55yaN8LFizAnTt3Sl3n6HuXlpaGkSNHYvTo0Th37lylP6NPfsZL+u6777Bz584qx+YoR79zo0ePRn5+fpW2rY73va569dVXkZOTUy37Epe10mg0YubMmfj888/RoUMHAMDevXsxffp0HDlyBCKRqFqCqIrbt28jIyMDABAREYEPP/ywxo+ZkZGBkydP4uLFi07926vDxYsXcefOHSxcuNDZoTjVm2++ibfffhvffvstBAJBte77zJkz8Pf3x44dOwAA3bt3r9T2T37GSzp//jxat27taIhV5uh3bu/evdUYTcORkJBQbfsqM/nrdDoUFBSgsLCQLxs1ahQUCgUsFgtEIhGOHj2KrVu3wmQyQSaTYfHixXjqqacQFxeH5ORkpKenQ6VSITw8HGvXroVCocDvv/+OTz/9FEajETk5OXjuuecwb948nDlzBmvXroVcLodWq8X333+PDz74AJcuXYJWqwVjDGvWrEHjxo3x4YcfoqCgAO+88w6ee+45vPfee9i3bx8KCgqwatUqXL9+HQKBAP369cNbb70FsViMiIgIzJgxAwkJCcjMzMS0adMwYcIEm7/73Llz+OCDD6DT6eDm5oZ58+aha9eumDZtGsxmM2JjYxEXF4fQ0FB+G61WizVr1uDChQsQiUQYPHgw5s+fb7Xf3bt349tvv4XJZIJarcb06dMxYcIEqFQqLF68GLm5uQCAAQMGYN68eXbLgaKa33/+8x9wHAdvb2+8++67aNmyJc6dO4cNGzaA4zgAwMyZMzF06FCbvzEuLg6TJk0CgFLP+8mTJ0t9XzUaDZYtW4br168jICAAIpEI3bp1s9q3xWLBoEGD8PHHH6Njx44AgHnz5qFnz57o1asXli1bBqPRCMYYxo0bh4kTJ5b1MQQA/Prrr4iPj4der0dMTAxmzZpl8/fk5ubyNeUnlwsKCrB27VrcvHkTJpMJkZGRWLRoEcRiMZo2bQpPT08cOXIEgwcPtjnu+fPncejQIWg0GvTp0weLFy/G/v378fXXX+Obb74BADx69AgvvPACjh49ColEAgA4ffo0Nm/ejIKCAkyePBmzZ8/mP6NLlixBXl4eHj58iKeffhoDBw60ec86depk9Rlfv3691bk4evQoEhISYDQaER8fj4SEBMjlcqxYsQJ3797lr0yjo6OxdetWcByH1atXIy8vDwKBAK+++iqee+45m7930KBBGDFiBBISElBQUIBXXnkFEyZMsPmMLFq0CO+//z7/9ygUCty4cQPp6elo27Yt3n//fXh4eODSpUtYs2YN/11atGgRIiMj0bZtW5w6dQrHjh3DwYMHwXEcHj16hMDAQGzYsAGBgYG4ePEi/vGPf8BoNEKlUiEqKgrr1q0r83Py+++/Y/PmzeA4DnK5HKtWrUJ4eDh+++03fPTRR+A4Dh4eHnjnnXfQqVMnxMXF4cGDB8jIyIBKpUKHDh3Qq1cv/Pjjj0hJScHChQsxcuTIMvPZrVu3Sj23Z86cwaZNm9C0aVPcunULZrMZq1atQrdu3WA0GrFx40acPXsWFosF7du3x/Lly6FQKDBo0CCMGTMGp06dQlpaGkaPHo158+bhnXfeAQBMnToV8fHx+P333/HNN9/Azc0NUqkUq1evRqtWrcr9LvFYOT7//HPWqVMnNmjQILZgwQL23XffscLCQsYYY/fu3WMjR45kOTk5jDHGbt68yfr06cO0Wi378MMPWf/+/ZlKpWIWi4W99dZbbMOGDYzjODZp0iR27949xhhj6enprF27diw7O5udPn2ahYeHs5SUFMYYYxcuXGBz5sxhFouFMcbYp59+ymbOnMkYY+z7779nM2bMYIwxdvr0aTZixAjGGGOLFi1i7733HuM4jhkMBvbqq6+yTz/9lDHGWJs2bdiXX37JGGMsKSmJdezYken1equ/Nycnh0VGRrKLFy/yf1PPnj3ZgwcP2MOHD1mXLl1KPU/r1q1j8+fPZ2azmRkMBjZx4kR2+vRptnjxYvavf/2LaTQa9sILL/Dn6r///S+/r48++oi9++67jDHGtFotmzdvHsvPz7dbfubMGTZhwgT+fThx4gQbNmwYY4yxKVOmsH379jHGGLt27Rr7+9//bhOrWq1mnTt3ZgaDgT9/T573st7XtWvXskWLFjGO41h2djbr378/+/DDD22OsWXLFrZq1SrGGGN5eXmsZ8+eLD8/n73zzjv8+5GZmcnmzZvHv7/2TJo0ic2cOZOZTCZWUFDAhg0bxo4dO2b1fnz44Yf88UouL1myhH3xxReMMcbMZjNbsGABi4+P51/7+eefs0WLFtkcd/HixWzMmDFMq9Uyg8HAJk2axHbu3MkMBgOLjIxkN2/eZIwxtnnzZrZx40ab7e19RhcvXsymTp3Kv87ee/bk9qXF9q9//YsxxtjkyZPZ0aNHGWOMRUdHs6ioKKbRaNitW7fY8OHDmclkYs888ww7dOgQY6zoO9evXz924cIFm/0OHDiQvfvuu4zjOJaWlsZ69erFrl+/bvMZKfn3vPjii8xgMDCj0ciee+45tnv3bmY0GlmfPn3Y77//zhgr+s6NHDmSWSwW1qZNG5adnc2+//571qVLF3b37l3GGGP/+Mc/2Jw5cxhjjM2fP5+dPn2aMcaYRqNhvXr1YklJSXa/hyqVinXr1o1duXKFMcbYoUOH2GuvvcZu377NoqKi2IMHDxhjjCUmJrI+ffqwgoIC9uGHH7KBAwey/Px8ptPpWI8ePdj69esZY4z9+uuvLDo6mjHG7Oazss7t6dOnWbt27djVq1cZY4x99tlnbOLEiYwxxuLi4vh8yBhj//d//8dWrlzJvwcbNmzg9xcREcHHXnzezGYz69ChA8vIyGCMMfbDDz+wb775ptTPij1l1vwB4JVXXsHzzz+Ps2fP4uzZs9i+fTu2b9+O3bt38zXol19+mX+9QCDAgwcPAADDhg2Dv78/AGDcuHFYt24dFi9ejG3btuHYsWPYt28f7ty5A8YYdDodACA4OBghISEAgKeeegqNGjXCN998g4cPH+LMmTPw8PAoM94//vgD//nPfyAQCCCRSDB+/Hj8v//3/zBjxgwAwDPPPAMA6NChA4xGIwoLCyGVSvntL1++jNDQUHTu3BkA0Lp1a3Tt2hV//vknevXqZfe4iYmJeOeddyASiSASifia1w8//AAA8PDwwLZt23D8+HHcv38f169f56+o+vXrhxkzZiAtLQ1RUVF4++234enpabf82LFjSE5Oxvjx4/nj5+fnIy8vD8OHD8fq1atx9OhRREVF4a233rKJNTk5GUqlkq+lljzvZb2vp06dwtKlSyEQCODr64shQ4aUej7Gjh2LcePGYcmSJdi3bx8GDRoET09PDBkyBIsXL8bly5cRGRmJ5cuXQygs/7mDcePGQSwWQ6FQYOjQoUhMTETLli3L3Q4Ajh07hqSkJOzevRsAoNfrrdY3adIEBw4cKHXb0aNHQy6XAyi66j1+/DgmTJiA559/Ht999x0WL16MH374AV9++WWFYin25NVSRd6zsgwZMgR//PEHQkNDERgYiDZt2uDs2bO4ceMGoqOjcf/+fRgMBkRHRwMAAgMDER0djRMnTpTadj5hwgQIBAIEBQWhX79+SEhIQIcOHaw+IyX169eP/zy1adMGarUaN2/ehFAoxNNPPw0A6NixI37++Webbfv06YPmzZsDAF544QWMHj0aALBhwwb88ccf2LZtG+7evQuDwYDCwkJ4e3uXGsOFCxfQunVrtG/fHkDRVU90dDR27tyJ3r17o2nTpgCAyMhI+Pr64q+//gIAREVFwdPTEwAQEBCAfv36AQBCQ0ORl5fH77+0fDZ27Fi757ZXr15o3Lgx2rVrBwBo3749nw+OHTuGgoICJCYmAgBMJhP8/Pz4YxXnqcDAQPj5+UGtVvPxA4BIJMKwYcMwfvx4PP300+jbty8GDBhQ6nmxp8zkf/78efz3v//FtGnTMHDgQAwcOBBvvfUWRo4ciYSEBHAch8jISGzevJnfJi0tDQEBAfj111+t2sU5joNQKERhYSHGjBmDwYMHo3v37hg7dix+++03sP8NMVT8RSs+QWvXrsUrr7yCZ555Bi1atMBPP/1U5h/EcZxV2y3HcTCbzfxycaIvfg0rMbSRxWKxaftljFntozRisdhqu7S0NMhkMn45PT0dL774Il544QV069YNw4YNw++//w4A6NSpE44cOYJTp07h9OnTeP7557F9+3a75RzHYfTo0Xx7PcdxyMzMRKNGjTB+/HgMHDgQCQkJOHHiBD766CMcPHjQ6gdOIBDwTQzFnjzvZb2vJc+ZvXsfISEhaN++PY4dO4Y9e/Zg6dKlAICBAwfi0KFDSExMxKlTp/Dxxx9jz549CAoKKvP8PnkcxhjEYuuPrkAgsIrLZDJZ/T1btmzhfyzy8/Ot3iuxWGz3B8jeccePH49x48ahZ8+eaN26tdUXsyKePN/23rOKGjJkCCZOnIhmzZqhT58+8PLywsmTJ5GUlIRVq1ZV+jP95Lkt/t6WjLmkJz/rxe+FSCSyOe7NmzfRokULq7KSeaJ4edKkSWjbti369euH4cOH49KlSzbf15L7efJ4jDHcuHHDJicUryv++5+sBJX8+8uKUygUlntuSzsvxdsvXbqUT9harRYGg4F/bcnva2l/98aNG3Hz5k0kJiYiPj4ee/fuxZYtW0qNvTRlVrl8fX2xdetWnDt3ji9TqVTQaDRo06YNIiMjkZCQwD8tcfz4cYwaNYqvWR05cgQFBQXgOA67du3CwIEDkZycDI1Gg3nz5mHQoEE4c+YMjEajTTICimqgAwcOxIQJE9CxY0f89ttvsFgsAIreiNI+vH379sVXX30FxhiMRiN27dqFqKioCp+QLl264O7du7h8+TIA4NatWzh79ix69uxZ5naRkZH44YcfwHEcjEYj5s6di7Nnz/Lr//rrL/j6+uKNN95A3759+cRvsViwceNGfPLJJxg8eDCWLVuGVq1a4datW3bL+/bti19++QWZmZkAgP/85z+YOnUqgKJEcu3aNcTGxuK9995Dfn4+VCqVVayhoaHIzs62+rCV/Fvsva/9+vXD7t27wXEc1Go1jhw5YvecvPDCC9i+fTt0Oh1f03377bexf/9+jBgxAitXroRCoeCvFMvy448/gjEGtVqNAwcO8LWzYj4+Prhy5QoYY9BoNPz5BYo+Ezt27OA/E7NmzbJ6WislJcUmIRX75ZdfYDQaYTAY8MMPP6B///4Aiq6UunTpgnXr1uGll14qN/6y2HvP7H3GAevPf1BQEHx8fPDNN9+gT58+6Nu3Lw4fPoy8vDyEh4ejRYsWEIvFOHz4MICiBxcOHTpk93vx448/Aii6l5GQkMD/zZXVokULCAQC/ibllStXMHXqVJvv+unTp/kb29988w0GDhyI/Px8JCUlYcGCBYiOjkZ6ejoePHhQap4o1rlzZ9y5c4d/GvHIkSNYuHAhIiMjcfLkSTx8+BAA+Lb04qv7iiotn1X23Bbr27cvdu7cyee+d999F//85z/LjaH4fc/JycGAAQPg7e2Nl19+GfPmzUNSUlKl/p4ya/7NmzfHxx9/jE2bNiE9PR1SqRSenp5Yt24d/2VZvXo13nrrLb5WtHXrVr5pxt/fH9OnT0dubi569OiB119/HRKJBE8//TSGDx8OiUSCNm3aoFWrVkhOTrb5BR4/fjzefvttxMTEwGw2o0+fPjh8+DA4jkOXLl3w8ccfY/bs2Zg8eTK/zfLly7FmzRrExMTAZDKhX79+eP311yt8Qnx9fbFlyxa899570Ov1EAgEWL9+PZo3b46UlBS7282ePRtr167F6NGjYbFY8OyzzyI6OhpHjx4FUHRpu3v3bgwbNgwCgQA9e/aEr68vkpOTMXXqVCxZsgQjR46ERCJB27ZtMWLECKjV6lLLJRIJpk+fjldffRUCgQAKhQIfffQRBAIBFixYgHXr1mHz5s0QCASYPXs2mjRpYhWrl5cXunXrhtOnT5d6qdiqVSu77+ucOXOwcuVKDB8+HL6+vmjTpo3dczJo0CCsWrUK06dP58veeOMNLFu2DN9++y1/Y7xHjx7IyMjAjBkzEB8fj8DAQJt9eXp6IjY2Fnq9HpMmTULv3r2t3o9Ro0bhxIkTiI6ORmBgIHr27MnXlpYtW4a1a9fyn4moqChMmzaN3/bEiRP8ze+SmjRpggkTJkCr1WLIkCEYM2YMv644WVf2crske++ZxWLhP+MfffSR1Tb9+/fHhg0bABTdIB4yZAg+//xztG/fHkKhEDKZjL+B7ebmhk8++QRr1qxBXFwcLBYL/va3v6F3796lxpOSksKf6+XLl6NFixY2FYiKkEgkiIuLw7p16/DBBx/Azc0NcXFxNt/zwMBALFy4ECqViv/seXl5YcaMGRgzZgzkcjkCAwPRtWtXJCcn273K8vf3x8aNG7F48WJYLBYoFAps2rQJrVq1wsqVKzF79mxYLBbIZDJs27aNb+qpqNLyWVnn9syZM3b39cYbb+D999/HmDFjYLFY0K5dOyxZsqTcGIYNG4bJkycjLi4Os2bNwssvvwyZTAaRSIQ1a9ZU6u8RsLKuoxxQ8ukL4louXLiAbdu2IT4+3tmhONWDBw+wYMGCSj/qWfz0TOPGjfn7SfXBoEGDsGXLFkRERNTK8fbs2YNDhw7h008/rZXjVVV9zGfUw7eB6tq1K5o3b44//vjD2aE41ebNm7FmzZpKJX6NRoNevXohLS0NU6ZMqcHoCKk5NVbzJ4QQUn00Gg3Gjx+Pbdu22TTlXrt2DcuWLYNWq0X37t2xatUquzeui1HNnxBCXNylS5fw0ksv4f79+6WuX7hwIVasWIFDhw6BMYZdu3aVu09K/oQQ4uJ27dqFlStX8o9bPyk1NRV6vR5dunQBUPQgQkUeFS63kxchhJDql5+fX+rgdl5eXvDy8rIqW7t2rd39ZGZmQqlU8stKpdLumFBPouRPCCHVyGg0QyIpP7VKJBLExsZCrVZblc+ePRtz5syp8PFKdmJjjFXoAYYGkfxzc7XguMrd1/bzUyA7W1NDEVUexVM2V4rHlWIBKJ7yFMcjFArg41P28DEVIZGIkXDsDnQ6k93XuLu7oc/TLbF3716+42qxkrX+8gQFBVn1w8jKyiq1eaikBpH8OY5VOvkXb+dKKJ6yuVI8rhQLQPGUp7rj0elM0BXaT/7FgoODHT5WSEgIpFIpzp8/j27dumHv3r0V6pVNN3wJIaQOmj59Oj+kw8aNG7F+/XoMGzYMhYWFFep/0iBq/oQQUh8UDxcDANu3b+f/Pzw8nB+1tqKo5k8IIQ0QJX9Sr5i1OujSsmExlN/eSkhDRs0+pN7I+fMaCm48BBiDQCyCb49wKFo3KX9DQhogqvmTeqHgfjoKrj8A/jdUFTNbkH36KsxafTlbEtIwUfIn9UJhSinjzTMGfXp27QdDSB1AyZ/UC2KFe+nlHqWXE9LQUfIn9UKj8FCI5FKrMmmAD2RBvk6KiBDXRjd8Sb0glkkQNLw3Cq4/gClfC1mANxRtQ50dFiEui5I/qTfEHjL4dLM/pzAh5DFq9iGEkAaIkj8hhDRAlPwJIaQBouRPCCENECV/QghpgCj5E0JIA0TJnxBCGiB6zp8QQqrZzSwN1AUGu+sbeUoxuBbjKQ0lf0IIqWZB7f3haTDbXe8hdX7qdVqzj0ajwciRI5GSkgIASExMRExMDKKjo7Fp0yb+ddeuXUNsbCyGDh2KZcuWwWy2f0IJIYRUjFOS/6VLl/DSSy/h/v37AAC9Xo+lS5fik08+wf79+/HXX3/h+PHjAICFCxdixYoVOHToEBhj2LVrlzNCJoSQesUpyX/Xrl1YuXIlAgICAACXL19GWFgYmjZtCrFYjJiYGBw8eBCpqanQ6/Xo0qULACA2NhYHDx50RsiEEFKvOKXhae3atVbLmZmZUCqV/HJAQAAyMjJsypVKJTIyMip9PD8/RZXiVCo9q7RdTaF4yuZK8bhSLADFUx5Xi6c2OP+uAwCO4yAQCPhlxhgEAoHd8srKztaA41iltlEqPaFSFVT6WDWF4imbK8XjSrEAFE95iuMRCgVVrijWRS7xnH9QUBBUqsfT8KlUKgQEBNiUZ2Vl8U1FhBBCqs4lkn/nzp1x7949JCcnw2KxYN++fejfvz9CQkIglUpx/vx5AMDevXvRv39/J0dLCCF1n0s0+0ilUmzYsAFz5syBwWDAgAEDMGzYMADAxo0bsXz5cmg0GnTo0AFTpkxxcrSEEFL3OTX5Hz16lP//yMhI/PTTTzavCQ8Px+7du2szLEIIqfdcotmHEEJI7aLkTwghDRAlf0IIaYAo+RNCSANEyZ8QQhogSv6EENIAUfInhJAGiJI/IYQ0QJT8CSHExf3888949tlnER0djZ07d9qsv3LlCsaOHYtRo0Zh5syZyM/PL3eflPwJIcSFZWRkYNOmTfj666/x448/4ttvv8Xt27etXrN27VrMnTsXP/30E5o3b47PPvus3P1S8ieEECdJS0tDSkqK1b+StfbExET07t0b3t7ekMvlGDp0qM2kVhzHQavVAgB0Oh1kMlm5x3aJgd0IIaQ+SdfnQa032l3fiEkAABMnTkRqaqrVutmzZ2POnDn8cmmTXV2+fNlqmyVLluDVV1/FunXr4O7uXqHpbin5E0JINevUxB16s/30KhO7AQB27twJi8Vitc7Ly8tqubxJrfR6PZYtW4YdO3agU6dO+Pe//43FixcjPj6+zBgp+RNCiJMEBweX+5qgoCCcO3eOXy6e7KrYzZs3IZVK0alTJwDAiy++iC1btpS7X2rzJ4QQFxYVFYVTp04hJycHOp0Ohw8ftprUKiwsDOnp6bh79y4A4MiRI4iIiCh3v1TzJ4QQFxYYGIj58+djypQpMJlMGDduHDp16oTp06dj7ty5iIiIwPr16zFv3jwwxuDn54d169aVu19K/oQQ4uJiYmIQExNjVbZ9+3b+/wcMGIABAwZUap/U7EMIIQ2QyyX/vXv3YsSIERgxYgTef/99AMC1a9cQGxuLoUOHYtmyZTCbzU6OkhBC6jaXSv46nQ5r167Fl19+ib179+LcuXNITEzEwoULsWLFChw6dAiMsQo9w0oIIcQ+l0r+FosFHMdBp9PBbDbDbDZDLBZDr9ejS5cuAIDY2Fib3m2EEEIqx6Vu+CoUCrz55psYPnw43N3d0aNHD7i5uVn1blMqlcjIyKjUfv38FFWKR6n0rNJ2NYXiKZsrxeNKsQAUT3lcLZ7a4FLJ//r16/j+++/x+++/w9PTEwsWLEBCQkKZvdsqIjtbA45jldpGqfSESlVQqW1qEsVTNleKx5ViASie8hTHIxQKqlxRrItcqtnn5MmTiIyMhJ+fHyQSCWJjY3HmzBmoVCr+NVlZWVa92wghhFSeSyX/8PBwJCYmorCwEIwxHD16FD179oRUKsX58+cBFD0N9GTvNkIIIZXnUs0+ffv2xdWrVxEbGws3NzdERERgxowZGDJkCJYvXw6NRoMOHTpgypQpzg6VEELqNJdK/gAwY8YMzJgxw6osPDwcu3fvdlJEhBBS/7hUsw8hhJDaQcmfEEIaIEr+hBDSAFHyJ4SQBoiSPyGENEAu97QPIYRUhkmtQd7FOzDm5EPq54VGXVrBzcvD2WG5PEr+hJA6y6I3IuPQWVj0RgCAuaAQ+owcNB7dF0KJm9PiKmS50DG93fUck9ViNKWj5E8IqbMK76fzib+YRWdEYXIGFK2bOCkqoKVHI5g5ud31YqHzfpiKUZs/IaTO4kylT+xkr5w8RsmfEFJnyUMDgBKj/AqEQshDA50UUd1ByZ8QUme5NVLAv28ERHIpAEAkl8KvbwTECncnR+b6qM2fEFKneTQPhjwsEBadESJ3CQRCqtNWBCV/QkidJxAKIfZw/hM0dQn9RBJCSANEyZ8QQhogSv6EENIAUfInhJAGqErJ32g0Ii4uDu+88w6OHz9ute69995zKKCjR48iNjYWw4cPx5o1awAAiYmJiImJQXR0NDZt2uTQ/gkhhFQx+a9atQr3799Hy5YtsXz5csTHx/PrLly4UOVgHj58iJUrV+KTTz7BTz/9hKtXr+L48eNYunQpPvnkE+zfvx9//fWXzQ8OIYSQyqnSo56XL1/GTz/9BIFAgOHDh2PSpElo0qQJnn32WTDGqhzMr7/+imeffRZBQUEAgE2bNiE5ORlhYWFo2rQpACAmJgYHDx7EgAEDqnwcQghp6Kr8nL/gf12qQ0JCEBcXh9deew1hYWF8eVUkJyfDzc0Nr7/+OtLS0vD000+jdevWUCqV/GsCAgKQkZFR5WMQQgipYvKPiIjAkiVLMGvWLISFhaFjx45YsWIFpk+f7lAwFosF586dw5dffgm5XI5Zs2ZBJpNZ/aAwxir9A+Pnp6hSPEqlZ5W2qykUT9lcKR5XigWgeMrjavHUhiol/5UrV2Lr1q1ITU1FWFgYAGDEiBFwd3fHBx98UOVg/P39ERkZCV9fXwDA4MGDcfDgQYhEIv41KpUKAQEBldpvdrYGHFe55iil0hMqVUGltqlJFCQkVbcAACAASURBVE/ZXCkeV4oFoHjKUxyPUCiockWxLqrSDV+pVIp58+YhKirKqnzQoEE4ePBglYMZOHAgTp48ifz8fFgsFpw4cQLDhg3DvXv3kJycDIvFgn379qF///5VPgYhhBAXG9unc+fOmDZtGiZMmACTyYQ+ffrgpZdeQosWLTBnzhwYDAYMGDAAw4YNc3aohBBSa37++Wds3boVZrMZU6dOxcSJE63W3717FytXroRarYZSqcQ///lPNGrUqMx9ulTyB4Bx48Zh3LhxVmWRkZH46aefnBQRIYQ4T0ZGBjZt2oQ9e/ZAIpFg/Pjx6NWrF1q1agWg6D7orFmzsGzZMvTv3x8bN25EfHw8Fi5cWOZ+HerhW9pjnWq12pFdEkIIeUJiYiJ69+4Nb29vyOVyDB061Kp5/cqVK5DL5Xxz+Ouvv25zZVAah5J/bGysTVlFDkoIIfWZUKSGUJRTxr+iSnJaWhpSUlKs/uXn51vtKzMzs8zH3R88eAB/f38sXboUY8aMwcqVKyGX258/uFiVmn2mTp2KpKQk6PV6dO3alS/nOA4RERFV2SUhhNQb3mJfcJzJ7nrh/yZwnzhxIlJTU63WzZ49G3PmzOGXOY4r83F3s9mMP//8E1999RUiIiKwefNmbNiwARs2bCgzxiol/48//hh5eXlYunQp1q9f/3hnYrHVLxQhhBD7du7cCYvFYlXm5eVltRwUFIRz587xyyUfd1cqlQgLC+Mr3iNHjsTcuXPLPXaVkr9CoYBCocAXX3wBo9EInU7Ht//n5+fD29u7KrslhJAGJTg4uNzXREVFIS4uDjk5OXB3d8fhw4etBtB86qmnkJOTg+vXryM8PBxHjx5Fhw4dyt2vQ0/7fPPNN1i3bh1MJhOf/AUCAa5du+bIbgkhhPxPYGAg5s+fjylTpsBkMmHcuHHo1KkTpk+fjrlz5yIiIgIff/wxli9fDp1Oh6CgoAp1thUwB0ZiGzx4MLZs2VKhXxlnoh6+1Y/isc+VYgEonvLURA/f3Pyb5bb5+3i1qZZjVZVDT/v4+/u7fOInhBBiy6Hk37dvX3z99dfIyMhAXl4e/48QQohrc6jNPz4+HkajEatXr+bLqM2fEEJcn0PJ//Lly9UVByGEkFrkULMPx3H47LPPsGTJEmg0Gnz66ac2z6wSQghxPQ4l/w8++AA3btzApUuXwBjDiRMnrDp9EUIIcU0OJf9Tp05hw4YNkEql8PT0xOeff46EhITqio0QQkgNcSj5i8ViCIWPdyGRSCAWu9wo0YQQQkpwKFO3adOGH5vi7t272LFjB8LDw6srNkIIITXEoZr/smXLcOXKFWRnZ+Oll16CVqvF0qVLqys2QgghNcShmr9CocC6deuqKxZCCCG1xKHkf/fuXWzfvh15eXlWs3pt27bNoaDef/995ObmYsOGDbh27RqWLVsGrVaL7t27Y9WqVXRfgRBCHORQs8+SJUvg4eGBIUOGYOjQofw/R5w6dQo//PADv7xw4UKsWLEChw4dAmMMu3btcmj/hBBCHKz563Q6LF++vLpiQV5eHjZt2oTXX38d169fR2pqKvR6Pbp06QKgaNrIDz/8EBMmTKi2YxJCSEPkUM0/LCwMmZmZ1RULVqxYgfnz5/Mz2ZScu1KpVFrNXUkIIaRqHKr5cxyHkSNHokOHDpBKpXx5Vdr8v/vuOwQHByMyMhJ79uzh91/W3JUVVdUxupVKzyptV1MonrK5UjyuFAtA8ZTH1eKpDQ4l/yFDhmDIkCHVEsj+/fuhUqkwevRoqNVqFBYWQiAQQKVS8a/Jysqymruyour7ZC6c0QRzoQFuXnIIhA5dzFVLPM7gSvG4UiwAxVOempjMRYxcMIHB7noBpHbX1RaHkv+YMWOQmpqKP//8E2azGT179kRYWFiV9vXvf/+b//89e/bgzz//xPr16zFy5EicP38e3bp1w969e9G/f39HQq531JfvQP3XPTCzBSK5FL692kPetPI/kISQamRuBJQxkxeEbrUXi70QHNn4xIkTGDt2LH777TccOXIE48aNw2+//VZdsQEANm7ciPXr12PYsGEoLCzElClTqnX/dZnuURbyLt4GMxeNpGopNCDrxGVYDGV86AghBA7W/Lds2YKvvvoKrVq1AgDcunULCxcuxODBgx0KKjY2FrGxsQCA8PBw7N6926H91Ve6h7Y325nZAv2jLHg0D3ZCRISQusKhmr/JZOITPwC0bt2axvOvRUKpxE65G0xqDXSpKnBGugoghNhyqOYvk8mQlJSEiIgIAEBSUhLc3d2rJTBStsIHGTDm5EOflg2JnxeEkqI2RDdvDxTcTIHuQdEjsQKxCH5RHeDRjK4ECCGPOZT8Fy5ciNdff52/yXvv3j1s2bKlWgIj9uWev4n8K/cAAOJGHtA9yoJXeBg8mgdDKJMg77+3+NcyswXZiVfg3tif/4EghBCHkn/37t3xyy+/4NKlS+A4Dl26dIGPj091xUZKYdEbUXAtmV8Wy2UQNwuGexMlfHqEQ3X8ks02zGyBQZUH9xClzTpSs6raN4WQmuZQ8rdYLPjll19w8uRJiEQi5Obm8jdqSc0wa3RgHGdTbsovBACIPWSlbifyoOa42qS5kwr1xdswa/WQBvjAt1c7SHwaXkci4rocuuG7Zs0aHDx4EM888wwGDBiA3bt3Y9OmTdUVGymFxEcBkcz2Rq8syBcA4BkearNeHhYEiXf1dF4h5StMz0F2wl8wa/UAAENmLjKPnC/1R5sQZ3Go5p+QkIBffvkFbm5FbcmjRo3CqFGjMH/+/GoJjtgSiETw7d0eWSeT+Of7pf6N4NW+6L6LWOGOoGd7o+D6A5i1OsiC/KBoHeLMkBuc/FspNmWWQgP0adnU9EZchkPJ39fXFxaLhU/+AoGAH5SN1Bx5aCBCxvpCn5YNkbsEskBfq/VihTt8urd1UnQVV/gwE+rLd2AuKIQ00Bc+3dvCzVPu7LAcZreNn9r+iQtxKPmHh4djwoQJiI2NhUgkwv79++Hj48MP1fDKK69US5DElkjqBo9mQc4Oo8oMWWqojl0E/jcJkO5hJky5BWj8XN9aG5+opjRq2xQP/7zJ/21A0Q9ycdMcIa7AoeRvMBjQtm1bXLlyBQDQpEkTAMDNmzcdj4zUa9o7qVbJESi6ma1Pz4F7Y38nRVU9ZEpvKAd0/t9VjQ7SoKKrmrr+o0bqF4eS//r166srDtLAMGZnlFV75XWMPDQQ8tBAZ4dBiF0OJf8zZ84gPj4earXaqpzG4iHl8WjRGJqb1jdGRXIpZMF+ToqIkIbFoeS/fPlyTJ48GaGhodUVD2kgZAE+8OvT0eZZeGoaIaR2OJT8/fz8aIhlUmWKliFQtAwB4zhK+oTUMoeS/6BBg7Bz507069cPYvHjXTVu3NjhwEjDQYmfkNrnUPLPzc3FP//5T6uRPAUCAS5cuOBwYIQQQor8/PPP2Lp1K8xmM6ZOnYqJEyeW+rpjx45h9erVOHr0aLn7dCj5//777zh58iT8/ev2o3mEEOKqMjIysGnTJuzZswcSiQTjx49Hr169rOZSAYrmOH///fcrvF+Hrrf9/Pzg60sdVwghpCrS0tKQkpJi9S8/P9/qNYmJiejduze8vb0hl8sxdOhQHDx40GZfy5cvx+zZsyt8bIdq/m3atMGECRMwcOBASCSPBxNzpGfvRx99hAMHDgAABgwYgEWLFiExMRHr16+HwWDA8OHDaewgQohLcytMAbPo7K4XiNwB33BMnDgRqampVutmz56NOXPm8MuZmZlQKh+PCRUQEIDLly9bbfPFF1+gffv26Ny5c4VjdCj56/V6NG/eHPfv33dkN7zExEScPHkSP/zwAwQCAaZNm4Z9+/Zh48aN+PLLLxEcHIyZM2fi+PHjGDBgQLUckxBCqpvJIAcz229YEYhlkAHYuXOnzdS3JcdH4zjOaryoknNE3Lx5E4cPH8aOHTuQnp5e4RirpYdvamoqzGYzP6NXVSmVSixZsoS/imjZsiXu37+PsLAwNG3aFAAQExODgwcPUvInhNR5wcHlT68aFBSEc+fO8csqlQoBAQH88sGDB6FSqTB27FiYTCZkZmZiwoQJ+Prrr8vcr0Nt/snJyRgxYgSee+45xMbGYvDgwbhz506V99e6dWt06dIFAHD//n0cOHAAAoHA5pInIyPDkbAJIaTOiIqKwqlTp5CTkwOdTofDhw+jf//+/Pq5c+fi0KFD2Lt3L+Lj4xEQEFBu4gccrPmvXr0a06ZNw5gxYwAA33//PVatWoUvvvjCkd3i1q1bmDlzJhYtWgSRSGTVrFSVafH8/Ko2kYlS6VozL1E8ZXOleFwpFoDiKY+rxfOkwMBAzJ8/H1OmTIHJZMK4cePQqVMnTJ8+HXPnzkVERESV9utQ8s/OzuYTPwCMHTsWO3bscGSXOH/+PObOnYulS5dixIgR+PPPP6FSqfj1JS95KhanBhxXuQHDlEpPqFQFldqmJlE8ZXOleFwpFoDiKU9xPEKhoMoVxZoWExODmJgYq7Lt27fbvK5JkyYVesYfcLDZx2KxIC8vj1/OyclxZHdIS0vD3/72N2zcuBEjRowAAHTu3Bn37t1DcnIyLBYL9u3bZ3XJQwghpPIcqvlPmjQJL774IoYPHw6BQID9+/dj6tSpVd7fZ599BoPBgA0bNvBl48ePx4YNGzBnzhwYDAYMGDAAw4YNcyRsQghp8ATM7sDqFXP69GmcOHECHMehf//+iIyMrK7Yqg01+1Q/isc+V4oFoHjKUxPNPgV3EsHMervrBWIZPFtGVcuxqsqhZp+MjAwcPHgQCxcuxPPPP48vv/zSqn2eEHsYx8GQpYZZY78jDCGk5jiU/BcvXowWLVoAAEJCQtCzZ08sXbq0WgIj9ZchMw+pe04gff9ppO75A6rjl8BKdHQhhNQsh5J/bm4uP56/VCrFyy+/TDV/UibGGLJOXIKl8PElcWFyOgquP3BiVIQ0PA7d8LVYLMjIyEBgYNFcpVlZWfbnZiX1GmMM+VfuQXs3DQKhAIrWTeHZtqnN60y5BTBrbdtCCx+q4NWheW2ESqqB5u4j5CfdhUVngCzYH749wiGSS50dFqkEh5L/yy+/jOeeew79+vWDQCBAYmIiFi1aVF2xkTok7/xN5F+9zy/nnLkKxnHwamc95IdIJgUEApuJ2kXuUjDGwOmMELpLKt2Rj9QeXaoK2SeT+OXC5HSYNYUIHuF6D3sQ+xxK/uPGjUPHjh1x+vRpiEQivPbaa2jTpk11xVardKkqFFx/AM5ggjwsCP79Ozg7pDqDMQbNrRSb8oLrD2yTv1wKjxbB0N55xJcJhEKIPd2R+v1xWAoNEMll8OneFh7Ngmo8dlJ5mtuPbMqM2fkwZudD4udVyhbEFTmU/AEgPDwc4eHh1RGL0+hSVMg8+nj2MUOWGhlCBlE7aoaoEMbAmW1v2Nq7iesX2QESXy/oUlQQySSQNVUi+0QSfzVgKdQj+2QSJH5ecPOU12jopAoYV3oxNfnWKTR5KoD8a8m2ZTcfwmIwOSGaukcgFEIeajvkhjys9Jq7QCiEV7swBA7pDv9+ncBp9TbNQIzjoHtAA/i5Io8WtnN0u3krIPVv5IRoSFVR8gfAGYw2ZYxjYGazE6Kpm3x7tYd7k6LRVwVCITyaB8P7qdYV2lYgLv0C1F45cS55aCB8e7bjb/DKGvshYFBXJ0dFKou+XQDcmwbAmGPd41Dq5wWxh7udLYhZb0RWQlJR0427FF4dmiFgUFdYDCYIBIBQ4lbhfcmbBUF96TYs+sc/wiKZBHJq83dZnuGh8AwPrdIou8Q1UPIH0Khjc5jzC6G9nw4wBjcfTwQ/0xX5JmrDtOfR4XP8TVvOYEJ2wl8QStwgb1q5EVcBQCR1Q+DQHsi7eAfGnHxIfL3g3aUlRNKK/4A4kzFPA92DjKK/v3n5k3PUJ5T46y5K/gAEIhH8+3WCT7e24MxmuHl5QOqtAFxo/BFXYlJroEu3HcFVczu1SskfANwaKaAcUPH5R12F5nYqshP/4pfVSXfgO3GQEyMipGKozf8JIrkUbl4ezg7D5TF7g+RxpT8FUl8xiwW5529YlVl0RmRfuOWkiAipOKr5k0qT+HhC6ueFggLrnrqlPQVSn5kLDeBKeSLMkJ1PX6wGTphxG0xvv+VAIPMEnDyqJ31GSZWERHeHZt8Z6B9lQySTwKtDM3g0sPZusYcMIpnE6kY1AHrkkcCMQHCwPzWkEM7vv0LJn1SJm6ccgYO7g3EcBMKG2XooEArh0yMcWScfd1ATyWXw79YGeXoapZS4Nkr+xCENNfEX82geDIl/I+geZhY97RMWWNQruYxLfkJcASV/Qhzk5imHW/tmzg6DkEqpM9W2n3/+Gc8++yyio6Oxc+dOZ4dTr1j0RuRfS4b6r7swqTXODocQUgvqRM0/IyMDmzZtwp49eyCRSDB+/Hj06tULrVq1cnZodZ4xT4PMw2f5m5Z5/70N/74RDe7mLSENTZ2o+ScmJqJ3797w9vaGXC7H0KFDcfDgQWeHVS+UHFYBjCH3/A2wBvbMPiENTZ2o+WdmZkKpVPLLAQEBuHz5coW39/NTVOm4SqX9R7WcoSbiKTCbIPKU2ZT7ekohltuW13Q8jnCleFwpFoDiKY+rxVMb6kTy5zjOagyRyg4mlZ2tAWevV6odSqUnVC40vENNxaMXS1BYYD1Ug0guQ47GCIHW/pDWDeX8VIUrxQJQPOUpjkcoFFS5olgX1Ylmn6CgIKuJ4VUqFQICqjaGDLHm3aUlRO4SflkgFMK3RzgN2EVIPVcnkn9UVBROnTqFnJwc6HQ6HD58GP3793d2WPWCWyMFGo/uC7/IDvDp3haNn+sLeVigs8MihNSwOtHsExgYiPnz52PKlCkwmUwYN24cOnXq5Oyw6g2hxA2K1k2cHQYhpBbVieQPADExMYiJiXF2GKSBYYzBkJELzmiCLNgPQrc685UhpEz0Sa6DjNn5MObkw83HkwYRq0EWnQEZv52HKbfo5qRQIoby6acgC/J1cmSEOI6Sfx2TffoKNDdT+GWPlo3h3yfCiRHVXRYzh0cpaqjz9JBKxWjcxAsenlJ+vfryHT7xAwBnNCM78S80HtOPboiTOq9O3PAlRfQZOVaJHwC0dx5Bl5Zd5naWQgNMBYVVOiZnMoOx+jmd5Y0rmUh9oIYm34BslRZXLqVDV/i4w5s+I7fovyaGTA2HbC0HfX4hLFq9vV0SUmdQzb8OMajySi/PzIV7sJ9NOWcyIzshCYUPMgEUjTPvP6BzhSamN+YWIOfUFRiy1EXj9Ue0gFe7MMf+ABeiLTAgX22dxDmOIf1RAZq3KjqXYk85Mh/lI1X9uLdzlk6EICaolS+Oxczh/t0cZGdqIRAKEBjsiabNvOmqg1QLqvnXIfammHRrVHrHFPXF23ziBwBDlhrZCX+V+tonMY5D5pELMGSpARQN/JZ79jp0qapytqw7TKbSx9s3mx4nes/2zZBR4oLJLdgPj1Jqp4PSvdvZUKUXdVC0mDk8eqhG6gN1rRybuJbyBrb87bffMHr0aIwaNQpvvPEG1OryPyeU/OsQ96YBNjcbpUpvyENL7/BW+DDTpkyfngPOaL/nLlDU3GEptG3a0N5Lr0S0rs2zkQwikW0N2sfv8VWRwEsBRYfmkAX7QRrgA8/wULiHKFGoNdpsV904jiFbpbUpV6XTqKsNTfHAll9//TV+/PFHfPvtt7h9+za/XqPR4O9//zvi4+Px008/oW3btoiLiyt3v5T86xCBQICAZ7rCr09HeLYLg19URwRGd7c7oYpQ6ma7D7EIAlHZb7u99eVtV5eIREK0aR8AN4kIACAQAEEhXvAPeHwVJZGKIfOSw6NZEBQtG0PiUzT+i4enpNR91gpq8WlwyhvY0mQyYeXKlQgMLOqc2bZtW6SlpZW7X2rzr2MEIhEULUOAluW/1rNdGLJPJlmXtWkKgUhU5nayAB+4+XhaPekCgQCKViFVCdllNfJxx1M9m0CvM8FNIoKbm/V5EQoFCGvph9vXHzd3SSQiNAnzrvHYhEIB/AI8kJVhXftXBjacsWfqMtOdu+AK7De9CD0bAVFAWloaLBbrJkgvLy94eXnxy+UNbOnj44MhQ4YAAPR6PeLj4zF58uRyY6TkX48pWjSGQCiE5lYKmNkCebMgeIaHVmjbwGe6Iff8DejTsyH2cEejTi0hVdZ80qttQqEAcg/7NXn/AA8oPCXIzS6EWCyCr78cInHtXAE1b+UHkVCIrEwNhEIBAoI9ERJK/TrqAqN3S1jcdHbXi/730MXEiRORmppqtW727NmYM2cOv1zRgS0LCgrwt7/9DeHh4RgzZky5MVLyr+c8mgXBo1lQpbcTyaXw70dDaACAzN0NwU1qP+mKREI0b+2H5q1tn+Qi9cPOnTtLrfk/KSgoCOfOneOXSxvYMjMzE6+99hp69+6NpUuXVujYlPwJqWVGgxlpKfnQaAzw8JAguGkjSKX0VWyIgoPLnzEvKioKcXFxyMnJgbu7Ow4fPoz33nuPX2+xWPD6669j+PDheOONNyp8bPrEkVIZstQoTM6AUCKGR4vGEHuUPbELqRiLhcOVS+kw6M0AgAK1ATnZhejULQTiWmpOInWLvYEtp0+fjrlz5yI9PR1Xr16FxWLBoUOHAAAdO3bE2rVry9wvJX9io+DmQ+Scvsov51+5h8DoHpD4epWxVcNh0Juhyih6/t7XXw7FE0NClCc3u5BP/MWMBguyVVoEBje82aRIxZQ2sOX27dsBABEREbh+/Xql90nJn1hhFgvUF29blXFGM9SX70D59FNOisp1aAsMuHI5HZylaMiLRw/VaNHaDwElErdZq0PuuZtFN8wV7mgU0QLy0EAYjaV3LjPZKSekptB1JrFi0RutJ3T/H1OebYej+spktODOzSz8988UXL2cjvy8xx3eUpLz+MRf7MG9XKtpQhljyPztPAqT08EZTDBm50N1/BIMmXnw8S19aA175YTUFEr+VWAxc8jJ0iIvR1fvBj0TyWWltu9LlA3nEcOrSelQpWtg0JuRn6fHtaR0vldvYaFt72izmYP5ieEidBm5MKlL/FgyBs2dVLjLJWjWyhfC//UuFgoFCG3hYzWaKCG1gZp9KilfrceNK5mwmIvGgJG5u6F9p0BI6snTGgKBAL692kN1/CKYpehvFHvI4N2llZMjqx35eXroSkxczxiQkVY04JvCU2rTZi+VivmewkBR01lpGFd0PoMaF/Uk1hWa4C53oxu9xClcJmOdP38e69evh8lkgre3N9atW4eQkBDk5+djwYIFePjwIXx9fbF582ar3m617e7NbD7xA4BeZ8LD+3lo2dbfaTFVN/cmSoTE9ocuNQsCsQjypspyewXXFxYLV2Z50+beKMjXw2goSvBCoQDNW/tadbqRB/tB7CGDucTQzx7NHz/WJxYL4elFtX3iPC5T5Vi4cCHWrFmDvXv3IiYmBmvWrAEAbN68Gd27d8eBAwfw/PPPl/v4Uk0yGszQ62wv+59sE64vRO5SKFqFwKNZUINJ/ADg5S0rtSbu5180oqpM5oYuPZqgdTslWrTxw1O9msDbV271WoFQCOWgrpAG+AAo6jDn26s93BvXnwoCqftcouZvNBrx5ptvIjw8HEDRwERfffUVAODYsWP8EKYjR47E6tWrYTKZ4OZmO2hZTRO7iSASC61q/gAgk7vEaaxTOKMJmtupMBcUQhrgA3mzIJcYp14kEqJtxwDcuZENvc4EkUiAkFBv+Pg9TvBCoQB+ytKH1y4m8fFE0LCeYBZLg/rxBMDfB3OF95PY5xJZSyKRYPTo0QCKxrH46KOPMHjwYADWgxqJxWIoFArk5OTwI9jVJqFQgCZh3ki+k8OXCQRASGj9G/OmJnFGE9L3n4Epv+imaMGNh5A/yIByQBcnR1bE00uGLj1CYDSYIXYTQSisehJrSInfYjAh98xVFD7IhEAkhKJNU3h3bU0/Ai6q1pP/gQMHsH79equyFi1aYMeOHTAajViyZAnMZjNmzpxZ6vaMMQjtDGFsj59f1UZCVCptO90olZ5o0tQHmekFEIuECG7iBQ9F7bTdlhaPM1U1npzLdyBjFsg8n3iqKEcNT1ggc2DwOFc6P64UC1A78aQePgdhdh4U/xsojz1IgyjQC35PtXZKPJXhavHUhlpP/sOHD8fw4cNtyrVaLWbNmgVvb29s3bqVb9YJCAhAVlYWgoKCYDabodVq4e1duQSRna2xeg67IpRKT6hU9mds8gsoagYo1BlRqKv5yT3Ki6e2VTYezmiC5lYKjLka6B6pYNEabOYHSLubAQWqVlN2pfPjSrEAtRMPZzQh7a/kokejnqA/fxtcE+uBBWs6Hr3ehIf38lCg1kPm7oYmYd7w8rY/PElxPEKhoMoVxbrIpW74hoWFYfPmzZBIHg+xO2DAAPz4448AgP3796N79+5Oae8nVceZLUg/cAa5529Ce/cRdA9VKLiebNNHoj4OGd1gCASlN+9U8irdURzHcO1yBrJVWhiNFuSri/pp6AprvoJW17hEm//Vq1dx5MgRtGrVih+HOiAgANu3b8ebb76JJUuWYMSIEfD09MTGjRudHC2pLO3dR1adnqT+jWDKLYBJrYXEu6im5f1Ua7h5yu3twiXpdSbkZBdCKBTAX+kBsZvrt+9zJjOYyQKRvHqbKoVuYsibBUF795FVeW1PAJSXq7Pph8EYoMrQIrS5E2dgc0Eukfzbt2+PGzdulLrO29sb27Ztq+WISHUyF5SY1EIggKJNUyhahkAa5ANpgE+dS/zZKi1uX1fxrRwp9/PQvnNQmRPDOBNjDHnnbqDg5kMwCweJfyP49+kIt0bV18zh27s9hBIxCu+nQyASwjM8FF7twqpt/xXB7DTvcnb6bzRkLpH8Sd2hz8yFLjkDgsBGMPt5Q+xR/pg0QhlmRAAAEY5JREFUsiAf5F+5Z1Pu2T6Mnxe3uuXl6KDVGOAud4OPn7xanzhhjCH5To5V87bZzCElOQ9t2gfY39CJNDceIv9aMr9szFJDdfwSGo/qU23HEIpF8O3ZDr4921XbPiurkY87RCIBLCXGX/ILKPvR3IaIkj+psPxrycg9WzR0LEuRQas3I3Bo+UM9u4cooWjTBJqbKUUFAgEadWpRY4n/1jUVslWPm5m8vGUI7xjo0CObTzIZLaWOzqktcN125cLkdJsyU54GxjwN3/RWH4jFQrTtGIh7t7KhKzTBzU2Ips184OlF81GURMmfVAhntkB9qcRQzyYz1El3K/R8vl/vDvBqFwZjrgZS/0YQK2pmFMt8td4q8QNFPbBzswvL7ZhVUW4SESQSkc0PgFzhmk0+ACBwK/2rLqwD9ykqy6uRDJ27h8BsskAkFlI/Azso+ZMK4fRGcEazTbnN6JVlcGukqNY25tJoNaXXvrUaY7Ulf4FAgLCWvlZt/mKxEE2bue7TSp5tm0KXorIqc28aUKFmu7rKmTfg8y7cgilHbXe9m28jNI6JrsWIbFHyJxUi8pBBJJfBUmg9jpGrPZ7pYeeGq0c118r9lB7wUEiQk1UIoahouAc3F65Fu4cooXy6C/KvJoMzGOHeJACNOrd0dlj1lrlxc5gUhXbXC7yc/4ADJX9SIQKBAH6920P1xyUwc1Fzh5uXB7w7u9ZQz17eMvgpPWza/J8cm6e6yNzd0Lhp3ZnnQB4aCHlo7Q+LQlwTJX9SYfxQz4+yoAzyRqFMBkEtd+KpiNbtlFAGKmrsaR9C6gNK/qRSRDIJFC0aQ6H0hM6FhjAoydvXHd40NSIhdrletY0QQkiNo+RPCCENECV/QghpgKjN30Vo76dBnXQPFq0OsmB/+HRvC7EH9UokhNQMSv4uQJ+eg6w/LvPLhcnpMOVr0TgmyolR1S2FWiOS7+TAbObg6y+vkUc7CalPKPm7AM2tFJsyU24BDKo8l+tE5Yq0BQZcu5QBtbpo9FBVhgZNwrzRJIzOHSH2UJu/K2ClD0NbcrITUrrUh2pYSgzZ++ihGhYzDeNLiD2U/F2AR4vGNmViTznV+itIrzPZlHEcK3XkTUJIEUr+LsC9iRK+PdvxsyvJgnwR8ExX6pVaQV6NbG+MS6QiyNypVZMQe+jb4SI8w0PhGR4KxnEuOWSCKwsJ88aj5HxoNAYAgEgkQIvW/vTjSUgZXC7LXL16FR07duSXjUYjFi5ciOHDh2PMmDG4c+eOE6OreZT4K8/NTYTukaFo3ykIbdor0bVXUxragZByuFTNX6fT4b333oPJ9LgN98svv4S7uzsOHDiAs2fP4p133sGuXbucGCVxVV7e1C+CkIpyqWrmhg0bMHXqVKuyY8eOYdSoUQCAHj16ICcnB48ePXJGeIQQUm+4TPI/cuQI9Ho9hg0bZlWemZkJpVLJLyuVSqSn285HSgghpOJqvdnnwIEDWL9+vVVZixYtoNFosGPHDpvXM8asbtwxxiCsZLu4n1/Vpg5UKmtmgvGqonjK5krxuFIsAMVTHleLp6Sff/4ZW7duhdlsxtSpUzFx4kSr9deuXcOyZcug/f/t3X9MVfX/B/AnF+6FCB2aIGiCDgUXkawIugYyjElwZcivRbDEShGWgm0CFxWcBfEjN7K5NSwHi2AVcPUOpy5MrAQWwWzSlHBThIgAQUSBy/31+vzB1/MFUfjI4F763NfjL84597zP676vvjicc8/rNTwMLy8vHDlyBBYW06d3gyf/4OBgBAcHT1pXUVGBoqKiSW8oLCwMZWVlWL58OXp7e+Hk5AQAuHPnDuzt7Z/qmP39D6DXP90DU3Z2i9C3gOrVczzTW0jxLKRYAI5nJg/jEYnMZn2iOJ96enpQWFgIhUIBiUSCmJgY+Pj4YO3a/++il5qaiuzsbHh6euLAgQP4/vvvERsbO+24C+KGb3R0NKKjo4VlNzc3KJVKAIC/vz+USiW8vLzQ1NQES0tLrFgx9aGo6YhEs/vK32z3my8cz/QWUjwLKRaA45mJSGQ2pzFJFk3/bbOH27u7u6HTTX4YcfHixVi8eLGwXF9fj9deew22tuMPfQYFBeH8+fPYs2cPAKCrqwsqlQqenp4AgIiICHz++ef/juQ/nXfeeQdZWVmQyWSQSCQoKCh46jGWLHl2VsdeaGcBHM/0FlI8CykWgOOZyVzHsz4mYMbXqFQqhIWF4d69e5PW79mzB3v37hWWH73vaW9vj6tXrz5xu52dHXp6emY8/oJM/n/++afws6WlJfLz840YDWOMzT21Wg2FQjFl/cSzfgDQ6/VT7ntOXJ5p+5MsyOTPGGP/6x69vPMkDg4OaGpqEpb7+vom3fd0cHBAX1+fsPzf3hddMF/1ZIwxNtXGjRvR0NCAgYEBjI6O4ocffsCmTZuE7StXroSlpSWam5sBAEqlctL2JzEjrhvMGGMLWnV1NYqKiqDRaBAVFYVdu3Zh165dSE5OhoeHB1pbW3Ho0CE8ePAA7u7uyM3NhUQimXZMTv6MMWaC+LIPY4yZIE7+jDFmgjj5M8aYCeLkzxhjJoiT/yN6e3uRkJCAbdu2ISYmBn/99RcAYGhoCAkJCQgODkZcXNyk79XOt4XS4Ka5uRlRUVEICwtDfHw8urq6ABh3bqqrqxESEoItW7agrKzMYMed6Pjx45DJZJDJZMIT6PX19QgNDcWWLVtQWFho8Jjy8/Mhl8sBjBf9ioiIQFBQEA4ePAitVmuwOC5evIiIiAgEBwcjOzsbgHHnRqlUCp/Vw4dHjTk/RkVskvj4eCovLyciovLyckpJSSEioiNHjlBRUREREZ06dUpYP99GRkYoJiaGXF1dhXVfffUVZWZmEhFRY2MjRUdHGySWgIAAun79OhERVVRUUGJiIhEZb27++ecfCggIoLt379Lw8DCFhobSjRs3DHLsh+rq6uitt96isbExUqvVtH37dqquriZ/f3/q6OggjUZD7733Hl26dMlgMdXX15OPjw+lp6cTEZFMJqMrV64QEVFGRgaVlZUZJI6Ojg7y9fWl7u5uUqvV9Pbbb9OlS5eMNjcjIyP06quvUn9/P2k0GoqKiqK6ujqjzY+x8Zn/BAMDA2htbUVMTAwAIDIyEvv27QMw3lQmNDQUALB161b8/PPPkzqOzZeF0uBGrVYjJSUF69evBzBefK+7u1uIxxhzM7HglbW1tVDwypDs7Owgl8shkUggFovh4uKC9vZ2ODs7Y9WqVbCwsEBoaKjB4hocHERhYSESExMBPL7ol6FiqampQUhICBwcHCAWi1FYWIhnnnnGaHOj0+mg1+sxOjoKrVYLrVYLCwsLo82PsXHyn6CzsxMrVqxAXl4eIiMjkZycDLFYDGBy8SQLCwvY2NhgYGBgXuNZSA1uJBIJwsLCAIzXEjl+/DgCAwOnxGOouXn0uMB4wav/pqDVXFq3bp2QONrb23Hu3DmYmZkZLa6srCx8+OGHQtmA2Rb9mgu3b9+GTqdDYmIiwsLCUF5ebtTPzMbGBikpKQgODoa/vz9WrlwJsVhstPkxNpOt7fO4pjLOzs64du0a9u7di4yMDFRUVEAul6O0tHTK/jSLpjJPE8t8N7iZTTwlJSVQq9WQy+XQarXYvXv3Y/ef63ieZLYFrebDjRs3sHv3bqSlpcHc3Bzt7e0Gj6uiogKOjo6QSqVCwTBjzpFOp0NTUxNKS0thbW2NpKQkWFlZGS2e1tZWVFVVoba2FosWLcL+/ftRV1e3YP4NGZrJJv/HNZXp6OhAeHg4AgLGy7Fu3bpVuEllb2+PO3fuwMHBAVqtFsPDw0J97fmIZb4b3DxtPAAwPDyMpKQk2Nra4osvvhD+KprPuZnOTAWvDKW5uRnJyck4cOAAZDIZGhsbJ930NlRcZ8+eRV9fn1AmeGRkBGZmZrMq+jUXli1bBqlUiqVLlwIAAgMDcf78eZibmwuvMeRndvnyZUilUjz33HMAxi/xnDx50mjzY2x82WcCJycnODg44KeffgIA1NbWwt3dHcB4U5nTp08DGP9P5uXlJSS/+RAdHY0LFy5AqVQKjW2USiVsbGyEBjcAZt3gZjZSU1Ph7OyMzz77bFLdEEPPzUMzFbwyhO7ubnzwwQc4evQoZDIZAGDDhg24deuWcNnjzJkzBomruLgYZ86cgVKpRHJyMjZv3ozc3NxZFf2aCwEBAbh8+TKGhoag0+nwyy+/4M033zTK3ADA+vXrUV9fj5GRERARLl68CG9vb6PNj7FxbZ9H3Lx5E4cPH8bdu3dhY2ODvLw8rF69GoODg5DL5ejs7MSiRYtw9OhRPP/88waLy83NTehzMDY2hqysLPzxxx+QSCTIzs4WfknNl2vXriE8PBxr164VeoPa29vjyy+/NOrcPK7glSFlZ2ejqqpK+CsMAGJiYrB69Wrk5uZibGwM/v7+yMjIMOjlBIVCgcbGRuTl5c2q6NdcqaysRElJCTQaDV5//XUcOnQIv/76q9Hm5sSJE1AoFBCLxfDw8MDhw4dx69Yto82PMXHyZ4wxE8SXfRhjzARx8meMMRPEyZ8xxkwQJ3/GGDNBnPwZY8wEcfJnbI7cv38f27dvn/Y13d3d8PPzM0j5C8amw8mfsTly7949tLS0PHH76dOnERcXh97eXgNGxdjjmWx5B2YclZWVKC4uhkgkwpIlS5Cfnw9HR0d89913KC0thUgkwrJly5CZmYk1a9ZALpfDysoKbW1t6O/vx+bNm2Fra4va2lr09fUhOzsbUqkUcrkclpaWaG1tRX9/v/BAkVgsRlNTEwoKCjA6OgqxWIx9+/Zh06ZNUCgUqKmpgUgkwu3bt2FlZYX8/Hy4uLjg/v37yMnJQVtbGzQaDaRSKdLS0mBhYQEPDw8kJCSgrq4Ovb292LlzJ2JjY5GRkQGVSoWwsDAoFIpJZQx6enpw4cIFnDx5ckqhPsaMwhh1pJlpun79Ovn4+NDff/9NRETFxcWUmZlJ9fX1FBgYSP39/UREVFVVRcHBwaTX6yk9PZ2io6NJrVZTb28vubq60tdff01ERCUlJfTuu+8SEVF6ejpt27aNHjx4QGNjYxQXF0elpaU0MDBAUqmUfv/9dyIiamtrI29vb+ro6KCqqip65ZVXqLu7m4iIPvroI0pLSyMiIrlcLhxHq9XS/v376cSJE0RE5OrqSqWlpURE1NLSQi+++CKpVCrq7OwkT0/PGefB1dVVeK+MGQuf+TODaWhogK+vLxwdHQEAO3bsAAAUFBQgJCREKAAWERGBnJwcoYtaQECAUHrX2toafn5+AMZrMQ0ODgrjh4eH49lnnwUwXgTvxx9/xKpVq+Dk5IQNGzYAGC/B/PLLL6OxsRFmZmZwd3eHg4MDAOCFF15ATU0NgPEeBS0tLaisrAQAqFSqSe/ljTfeAAC4u7tDrVZjZGRkbieLsXnGyZ8ZjLm5+aQaLiqVCl1dXdDr9VNeS0RCO71H66w8rC30uPEn7i8SiaDT6abUjXk4tlgshpWVlbDezMwM9H/VTvR6PY4dOwYXFxcA460qJ45jaWkp7PNwTMb+TfiGLzMYHx8fNDQ0CDc8v/32W3z66afw8/PD2bNnhW/AVFVVwdbWFs7Ozk81/rlz56BWqzE2NoZTp04hICAAnp6euHnzJq5evQpgvO7+b7/9Bm9v72nH8vX1RUlJCYgIarUaSUlJ+Oabb6bdx8LCAjqdjn8RsH8FPvNnBuPm5obU1FTs3LkTwHjXpE8++QTLly/Hjh07EB8fD71ej6VLl6KoqOipG8JYWVkhNjYWQ0NDCAoKQmRkJEQiEY4dO4aPP/4YKpUKZmZmyM3NxZo1a3DlypUnjnXw4EHk5OQgNDQUGo0GGzduFOJ+Ejs7O7z00kuQyWQoKyvDkiVLnip+xgyJq3qy/wlyuRzr1q3D+++/b+xQGPtX4Ms+jDFmgvjMnzHGTBCf+TPGmAni5M8YYyaIkz9jjJkgTv6MMWaCOPkzxpgJ4uTPGGMm6D/Op5xrqza+1wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.cm\n", "\n", "plt.scatter(X_train_pca[:, 0], X_train_pca[:, 1],\n", " c=Y_train, edgecolor='none', alpha=0.4,\n", " cmap=plt.cm.get_cmap('Spectral', 10))\n", "plt.xlabel('component 1')\n", "plt.ylabel('component 2')\n", "plt.title('Separation of classes (red vs. blue) by first two principal components')\n", "plt.colorbar();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q5. what is the reduction in dimensionality if we use the PCA transformed data instead of the original ?**\n", "\n", "1. 28\n", "2. 38\n", "3. 72\n", "4. 7101\n", "5. 7129" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q5 answer.** _Your answer goes here._ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classification" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have two options. On one hand, now that the dimensionality has been reduced, we can try a standard linear classifier. However, since the data has been transformed it is not clear how this relates to the original gene expression dataset. On the other hand, we can try a more powerful \"black-box\" classifier ensemble using boosting on the original data, and see how that turns out." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Logistic Regression Accuracy: 0.9117647058823529\n", "[[19 1]\n", " [ 2 12]]\n" ] } ], "source": [ "# Fitting a logistic regression classifier to the PCA-transformed version of the training set\n", "from sklearn.linear_model import LogisticRegression\n", "\n", "clf=LogisticRegression(penalty='l2',solver='liblinear',random_state=0)\n", "clf.fit(X_train_pca,Y_train)\n", "\n", "Y_pred=clf.predict(X_test_pca)\n", "\n", "# Metrics for evaluation\n", "from sklearn.metrics import accuracy_score\n", "from sklearn.metrics import confusion_matrix\n", "\n", "log_reg_ac = accuracy_score(Y_test, Y_pred)\n", "log_reg_cm = confusion_matrix(Y_test, Y_pred)\n", "\n", "print(f'Logistic Regression Accuracy: {log_reg_ac}')\n", "print(log_reg_cm)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Boosting Accuracy: 0.9117647058823529\n", "[[18 2]\n", " [ 1 13]]\n", "Estimators: [DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=1,\n", " max_features=None, max_leaf_nodes=None,\n", " min_impurity_decrease=0.0, min_impurity_split=None,\n", " min_samples_leaf=1, min_samples_split=2,\n", " min_weight_fraction_leaf=0.0, presort=False,\n", " random_state=209652396, splitter='best')]\n", "Estimator weights: [1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n", " 0. 0.]\n", "Feat. Importances: 7129\n", "Rank: 4846; Feature: 1.0\n" ] } ], "source": [ "# This trains an ensemble classifier of \"decision-stumps\", essentially\n", "# selecting a subset of up to 100 genes that can be used to predict the\n", "# class outcome (cancer type) of the patients in the test set.\n", "#\n", "from sklearn.ensemble import AdaBoostClassifier\n", "\n", "clf = AdaBoostClassifier(n_estimators=50, random_state=0) # 100\n", "clf.fit(X_train, Y_train)\n", "# Y_pred = classifier.predict(X_train)\n", "Y_pred = clf.predict(X_test)\n", "\n", "boosting_ac = accuracy_score(Y_test, Y_pred)\n", "boosting_cm = confusion_matrix(Y_test, Y_pred)\n", "\n", "print(f'Boosting Accuracy: {boosting_ac}')\n", "print(boosting_cm)\n", "\n", "print(f'Estimators: {clf.estimators_}')\n", "\n", "print(f'Estimator weights: {clf.estimator_weights_}')\n", "\n", "print(f'Feat. Importances: {len(clf.feature_importances_)}')\n", "for xrap in range(len(clf.feature_importances_)):\n", " xval = clf.feature_importances_[xrap]\n", " if xval > 0.0:\n", " print(f'Rank: {xrap}; Feature: {xval}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q6. Which classifier achieves the better accuracy?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q6 answer.** _Your answer goes here. [HINT: compare both accuracies and confusion matrices]_ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q7. Can you explain what has happened with the boosting classifier?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q7 answer.** _Your answer goes here. [HINT: you will need to read the SKLearn documentation for the \"AdaBoostClassifier\"]_ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Clustering" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "RangeIndex: 799 entries, 0 to 798\n", "Data columns (total 18 columns):\n", "gene_name 799 non-null object\n", "X 799 non-null float64\n", "X.1 799 non-null float64\n", "X.2 799 non-null float64\n", "X.3 799 non-null float64\n", "X.4 799 non-null float64\n", "X.5 799 non-null float64\n", "X.6 799 non-null float64\n", "X.7 799 non-null float64\n", "X.8 799 non-null float64\n", "X.9 799 non-null float64\n", "X.10 799 non-null float64\n", "X.11 799 non-null float64\n", "X.12 799 non-null float64\n", "X.13 799 non-null float64\n", "X.14 799 non-null float64\n", "X.15 799 non-null float64\n", "X.16 799 non-null float64\n", "dtypes: float64(17), object(1)\n", "memory usage: 112.4+ KB\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
gene_nameYAL022CYAL040C
X-0.63-0.05
X.1-0.66-0.15
X.2-0.24-0.58
X.3-0.82-0.58
X.4-0.840.18
X.5-0.89-0.67
X.60.570.29
X.70.640.2
X.80.540.02
X.90.661.32
X.100.280.27
X.11-0.010
X.12-0.26-0.41
X.130.11-0.05
X.1400.17
X.150.710.14
X.160.81-0.11
\n", "
" ], "text/plain": [ " 0 1\n", "gene_name YAL022C YAL040C\n", "X -0.63 -0.05\n", "X.1 -0.66 -0.15\n", "X.2 -0.24 -0.58\n", "X.3 -0.82 -0.58\n", "X.4 -0.84 0.18\n", "X.5 -0.89 -0.67\n", "X.6 0.57 0.29\n", "X.7 0.64 0.2\n", "X.8 0.54 0.02\n", "X.9 0.66 1.32\n", "X.10 0.28 0.27\n", "X.11 -0.01 0\n", "X.12 -0.26 -0.41\n", "X.13 0.11 -0.05\n", "X.14 0 0.17\n", "X.15 0.71 0.14\n", "X.16 0.81 -0.11" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.cluster import KMeans\n", "\n", "# Load the Spellman dataset\n", "df_sc = pd.read_csv(\"./datasets/gene800_Spellman.csv\")\n", "df_sc.shape\n", "df_sc.info()\n", "df_sc.head(2).T" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0123456789...789790791792793794795796797798
X-0.63-0.05-0.070.53-0.14-0.65-1.04-0.40-2.01...-0.960.73-0.290.71-0.630.12-0.020.11-0.41-0.36
X.1-0.66-0.15-0.58-0.080.210.010.29-0.99-0.57-1.25...0-1.67-0.57-0.34-0.380.550.760.270.180
X.2-0.24-0.58-0.05-0.250.841.521.45-1.06-0.320.37...-2.38-0.560.44-0.54-0.321.551.430.160.510.49
\n", "

3 rows × 799 columns

\n", "
" ], "text/plain": [ " 0 1 2 3 4 5 6 7 8 9 ... 789 \\\n", "X -0.63 -0.05 -0.07 0.53 -0.14 -0.65 -1.04 -0.4 0 -2.01 ... -0.96 \n", "X.1 -0.66 -0.15 -0.58 -0.08 0.21 0.01 0.29 -0.99 -0.57 -1.25 ... 0 \n", "X.2 -0.24 -0.58 -0.05 -0.25 0.84 1.52 1.45 -1.06 -0.32 0.37 ... -2.38 \n", "\n", " 790 791 792 793 794 795 796 797 798 \n", "X 0.73 -0.29 0.71 -0.63 0.12 -0.02 0.11 -0.41 -0.36 \n", "X.1 -1.67 -0.57 -0.34 -0.38 0.55 0.76 0.27 0.18 0 \n", "X.2 -0.56 0.44 -0.54 -0.32 1.55 1.43 0.16 0.51 0.49 \n", "\n", "[3 rows x 799 columns]" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "train_sc = df_sc.T\n", "train_sc2 = train_sc.drop(['gene_name'], axis=0)\n", "train_sc2.head(3)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,\n", " n_clusters=10, n_init=10, n_jobs=None, precompute_distances='auto',\n", " random_state=None, tol=0.0001, verbose=0)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "km_sc = KMeans(n_clusters=10)\n", "km_sc.fit(train_sc2)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "K = 2, Silhouette score = 0.2857984530259937\n", "K = 3, Silhouette score = 0.2687327804042268\n", "K = 4, Silhouette score = 0.3816832030529871\n", "K = 5, Silhouette score = 0.35723272513768023\n", "K = 6, Silhouette score = 0.3508881616379952\n", "K = 7, Silhouette score = 0.35318762639297197\n", "K = 8, Silhouette score = 0.3080070193259968\n", "K = 9, Silhouette score = 0.23401284029555927\n", "K = 10, Silhouette score = 0.24513774506853678\n" ] } ], "source": [ "from sklearn import metrics\n", "\n", "for k in range(2,11):\n", " km_sc = KMeans(n_clusters=k)\n", " km_sc.fit(train_sc2)\n", " lb_sc = km_sc.labels_\n", " ss = metrics.silhouette_score(train_sc2, lb_sc, metric='sqeuclidean')\n", " print(f'K = {k}, Silhouette score = {ss}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q8. Which setting for \"k\" achieves the best clustering ?**" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(64, 2)\n" ] }, { "data": { "text/plain": [ "NSCLC 9\n", "RENAL 9\n", "MELANOMA 8\n", "BREAST 7\n", "COLON 7\n", "OVARIAN 6\n", "LEUKEMIA 6\n", "CNS 5\n", "PROSTATE 2\n", "MCF7D-repro 1\n", "MCF7A-repro 1\n", "K562A-repro 1\n", "K562B-repro 1\n", "UNKNOWN 1\n", "Name: x, dtype: int64" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Load the labels for the NCI60 dataset\n", "df_nci_lb = pd.read_csv(\"nci60_labs.csv\")\n", "print(df_nci_lb.shape)\n", "df_nci_lb.x.value_counts(sort=True)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(64, 6831)\n", " 0 1 2 3 4 5 6 7 8 9 \\\n", "Unnamed: 0 V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 \n", "1 0.3 0.679961 0.94 0.28 0.485 0.31 -0.83 -0.19 0.46 0.76 \n", "2 1.18 1.28996 -0.04 -0.31 -0.465 -0.03 0 -0.87 0 1.49 \n", "\n", " ... 54 55 56 57 58 59 60 61 62 \\\n", "Unnamed: 0 ... V55 V56 V57 V58 V59 V60 V61 V62 V63 \n", "1 ... 0.01 -0.62 -0.38 0.0499805 0.65 -0.03 -0.27 0.21 -0.05 \n", "2 ... -1.28 -0.13 0 -0.720019 0.64 -0.48 0.63 -0.62 0.14 \n", "\n", " 63 \n", "Unnamed: 0 V64 \n", "1 0.35 \n", "2 -0.27 \n", "\n", "[3 rows x 64 columns]\n" ] } ], "source": [ "# Load the microarray data for the NCI60 dataset\n", "df_nci_dat = pd.read_csv(\"datasets/nci60_data.csv\")\n", "print(df_nci_dat.shape)\n", "print(df_nci_dat.T.head(3))" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Silhouette score = 0.20587971878618905\n", "Silhouette score = 0.18975672062168678\n", "Silhouette score = 0.04240091462684446\n", "Silhouette score = 0.19121437351322212\n", "Silhouette score = 0.19922871941255466\n" ] } ], "source": [ "# df_nci_dat = df_nci_dat.T\n", "\n", "# df_nci_X = pd.to_numeric(df_nci_dat.iloc[0:64,1:6831])\n", "df_nci_X = df_nci_dat.iloc[0:64,1:6831]\n", "\n", "from sklearn.cluster import AgglomerativeClustering\n", "\n", "avg = AgglomerativeClustering(linkage='average', affinity='euclidean', n_clusters=20) # 2 = 0.26170 10=0.1897 20 = 0.20587\n", "avg.fit(df_nci_X)\n", "\n", "nci_lbls = avg.labels_\n", "\n", "nci_ss = metrics.silhouette_score(df_nci_X, nci_lbls, metric='sqeuclidean')\n", "print(f'Silhouette score = {nci_ss}')\n", "\n", "linkages = ['average', 'single', 'complete', 'ward']\n", "for link_type in linkages:\n", " avg = AgglomerativeClustering(linkage=link_type, affinity='euclidean', n_clusters=10)\n", " avg.fit(df_nci_X)\n", " nci_lbls = avg.labels_\n", " nci_ss = metrics.silhouette_score(df_nci_X, nci_lbls, metric='sqeuclidean')\n", " print(f'Silhouette score = {nci_ss}') \n", "\n", "# from sklearn.metrics.cluster import adjusted_rand_score\n", "\n", "# act_lbls = df_nci_lb.iloc[0:64,1:2]\n", "#### act_lbls = df_nci_lb['x']\n", "### actual_lbs = pd.to_numeric(act_lbls)\n", "# print(act_lbls)\n", "# actual_lbs = pd.to_numeric(act_lbls['x'])\n", "# act_labels = pd.to_numeric(df_nci_lb.iloc[0:64,1:2], errors='coerce')\n", "# ari = adjusted_rand_score(act_labels, nci_lbls)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q9. Can you find a better clustering by changing the number of clusters (set by the 'n_clusters=' parameter ?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q9 answer.** _Add your code and output for the silhouette scores here._" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [] } ], "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }