Skip to content

Instantly share code, notes, and snippets.

@kantale

kantale/lesson_3.ipynb

Created April 8, 2019 16:49
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save kantale/384747eff773b996c50cbad3c5070a78 to your computer and use it in GitHub Desktop.
Save kantale/384747eff773b996c50cbad3c5070a78 to your computer and use it in GitHub Desktop.
Download ZIP
lesson_TEI_8_apr_2019
Display the source blob
Display the rendered blob
Raw
lesson_3.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"!conda install pandas\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"a= '''1 2 id=a;attr=mitsos\n",
"3 4 id=b;attr=kostas\n",
"'''"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"with open('test.csv', 'w') as f:\n",
" f.write(a)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 2 id=a;attr=mitsos\n",
"3 4 id=b;attr=kostas\n"
]
}
],
"source": [
"!type test.csv"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"df = pd.read_csv('test.csv', sep=' ', header=None)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>0</th>\n",
" <th>1</th>\n",
" <th>2</th>\n",
" <th>id</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>id=a;attr=mitsos</td>\n",
" <td>5</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>id=b;attr=kostas</td>\n",
" <td>5</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" 0 1 2 id\n",
"0 1 2 id=a;attr=mitsos 5\n",
"1 3 4 id=b;attr=kostas 5"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Index([0, 1, 2, 'id'], dtype='object')"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.columns"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" print (x)\n",
" return 5"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"df['id'] = df.apply(lambda x: x[2].split(';')[0].split('=')[1], axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"df['attr'] = df.apply(lambda x: x[2].split(';')[1].split('=')[1], axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>0</th>\n",
" <th>1</th>\n",
" <th>2</th>\n",
" <th>id</th>\n",
" <th>attr</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>id=a;attr=mitsos</td>\n",
" <td>a</td>\n",
" <td>mitsos</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>id=b;attr=kostas</td>\n",
" <td>b</td>\n",
" <td>kostas</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" 0 1 2 id attr\n",
"0 1 2 id=a;attr=mitsos a mitsos\n",
"1 3 4 id=b;attr=kostas b kostas"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 2 id=a;attr=mitsos\n",
"3 4 id=b;attr=kostas\n",
"\n"
]
}
],
"source": [
"print (a)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"b = {\n",
" 'col_1': [1,2,3,4,5],\n",
" 'col_2': [10,11,12,13,14]\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>col_1</th>\n",
" <th>col_2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>10</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>2</td>\n",
" <td>11</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>3</td>\n",
" <td>12</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>4</td>\n",
" <td>13</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>5</td>\n",
" <td>14</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" col_1 col_2\n",
"0 1 10\n",
"1 2 11\n",
"2 3 12\n",
"3 4 13\n",
"4 5 14"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame(b)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 2 id=a;attr=mitsos\n",
"3 4 id=b;attr=kostas\n",
"\n"
]
}
],
"source": [
"print (a)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"from collections import defaultdict\n",
"d = defaultdict(list)\n",
"for line in a.split('\\n'):\n",
" if not line:\n",
" continue\n",
" #print (line)\n",
" ls = line.split()\n",
" d['chr'].append(ls[0])\n",
" d['pos'].append(ls[1])\n",
" ids = ls[2]\n",
" ids_s = ids.split(';')\n",
" d['id'].append(ids_s[0].split('=')[1])\n",
" d['attr'].append(ids_s[1].split('=')[1])"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"defaultdict(list,\n",
" {'chr': ['1', '3'],\n",
" 'pos': ['2', '4'],\n",
" 'id': ['a', 'b'],\n",
" 'attr': ['mitsos', 'kostas']})"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>chr</th>\n",
" <th>pos</th>\n",
" <th>id</th>\n",
" <th>attr</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>a</td>\n",
" <td>mitsos</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>b</td>\n",
" <td>kostas</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" chr pos id attr\n",
"0 1 2 a mitsos\n",
"1 3 4 b kostas"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame(d)"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [],
"source": [
"import random"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [],
"source": [
"def get_allele_maf():\n",
" return random.random()/2"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.4266539052580529"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"get_allele_maf()"
]
},
{
"cell_type": "code",
"execution_count": 126,
"metadata": {},
"outputs": [],
"source": [
"def generate_genotypes(N, f):\n",
" allele_1 = int(2*N*f)\n",
" allele_2 = 2*N-allele_1\n",
" \n",
" genotypes = ['A'] * allele_1 + ['B'] * allele_2\n",
" random.shuffle(genotypes)\n",
" genotypes = [(genotypes[x], genotypes[x+1]) for x in range(0, len(genotypes), 2)]\n",
" random.shuffle(genotypes)\n",
" return genotypes"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B')]"
]
},
"execution_count": 127,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"generate_genotypes(10, 0.3)"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [],
"source": [
"def generate_population(N, SNPs):\n",
" return [generate_genotypes(N, get_allele_maf()) for x in range(SNPs)]"
]
},
{
"cell_type": "code",
"execution_count": 158,
"metadata": {},
"outputs": [],
"source": [
"def generate_two_populations(N1, N2, SNPs, pert):\n",
" \n",
" ret = []\n",
" for x in range(SNPs):\n",
" \n",
" #pert = 0.4 # # -0.2 + 0.2 \n",
" random_allele_frequency_1 = get_allele_maf()\n",
" random_perturbation = random.random()/(1/pert) - (pert/2) \n",
" random_allele_frequency_2 = random_allele_frequency_1 + random_perturbation\n",
" \n",
" if random_allele_frequency_2>0.5:\n",
" random_allele_frequency_2 = 0.5\n",
" if random_allele_frequency_2<0.0:\n",
" random_allele_frequency_2 = 0.0\n",
"\n",
" #print (random_allele_frequency_1, random_allele_frequency_2)\n",
" \n",
" gen_1 = generate_genotypes(N1, random_allele_frequency_1)\n",
" gen_2 = generate_genotypes(N2, random_allele_frequency_2)\n",
"\n",
" snp = gen_1 + gen_2\n",
" ret.append(snp)\n",
" \n",
" return ret"
]
},
{
"cell_type": "code",
"execution_count": 159,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[[('B', 'B'), ('B', 'A'), ('B', 'B'), ('B', 'A'), ('B', 'A')],\n",
" [('B', 'B'), ('B', 'B'), ('B', 'B'), ('B', 'A'), ('B', 'B')],\n",
" [('B', 'B'), ('B', 'B'), ('B', 'B'), ('A', 'B'), ('B', 'B')],\n",
" [('A', 'B'), ('B', 'B'), ('A', 'B'), ('A', 'B'), ('B', 'B')]]"
]
},
"execution_count": 159,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"generate_two_populations(2,3,4, 0.2)"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.09994860158815176"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"min([(random.random()/5) - 0.1 for x in range(1000)])"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [],
"source": [
"def numerical_modeling(samples):\n",
" return [[y.count('A') for y in x] for x in samples]"
]
},
{
"cell_type": "code",
"execution_count": 118,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[[0, 1, 0, 0, 2], [0, 1, 0, 0, 2], [0, 1, 1, 0, 0], [1, 0, 0, 2, 0]]"
]
},
"execution_count": 118,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples = generate_two_populations(2,3,4)\n",
"numerical_modeling(samples)"
]
},
{
"cell_type": "code",
"execution_count": 151,
"metadata": {},
"outputs": [],
"source": [
"samples_genotypes = generate_two_populations(100,100,1000)\n",
"samples_numerical = numerical_modeling(samples_genotypes)"
]
},
{
"cell_type": "code",
"execution_count": 129,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'B'),\n",
" ('A', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'B'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B'),\n",
" ('A', 'A'),\n",
" ('B', 'A'),\n",
" ('A', 'A'),\n",
" ('A', 'B')]"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples_genotypes[0]"
]
},
{
"cell_type": "code",
"execution_count": 152,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"samples_np = np.array(samples_numerical)"
]
},
{
"cell_type": "code",
"execution_count": 132,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1000, 200)"
]
},
"execution_count": 132,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples_np.shape"
]
},
{
"cell_type": "code",
"execution_count": 133,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[2, 0, 0, ..., 1, 2, 1],\n",
" [0, 0, 0, ..., 0, 0, 0],\n",
" [0, 0, 0, ..., 1, 1, 0],\n",
" ...,\n",
" [1, 1, 1, ..., 0, 1, 2],\n",
" [0, 2, 0, ..., 0, 1, 1],\n",
" [0, 0, 1, ..., 0, 0, 0]])"
]
},
"execution_count": 133,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples_np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HOW TO INSTALL scikit-learn\n",
"\n",
"```\n",
"!conda install scikit-learn\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 134,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.decomposition import PCA"
]
},
{
"cell_type": "code",
"execution_count": 153,
"metadata": {},
"outputs": [],
"source": [
"pca = PCA(n_components=2)"
]
},
{
"cell_type": "code",
"execution_count": 136,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(200, 1000)"
]
},
"execution_count": 136,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples_np.T.shape"
]
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,\n",
" svd_solver='auto', tol=0.0, whiten=False)"
]
},
"execution_count": 154,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.fit(samples_np.T)"
]
},
{
"cell_type": "code",
"execution_count": 155,
"metadata": {},
"outputs": [],
"source": [
"samples_np_proj = pca.transform(samples_np.T)"
]
},
{
"cell_type": "code",
"execution_count": 142,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(200, 2)"
]
},
"execution_count": 142,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples_np_proj.shape"
]
},
{
"cell_type": "code",
"execution_count": 146,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 156,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHXRJREFUeJzt3W2MHVd5B/D/c/cl4OKQrU3rJOtdY7mxwKYV3lW8KGohjZMm4AbJBEEcEAUFq5KpsEJECZGsylLVVggIElaRa1q1woQ2cSpXFilxJFM1qOt6r0kaO8GJ2WbtdRKwnU2J6jT7cp9+uDvr2dmZO29n3s78f1/I2uO557JnnnnOM2fOEVUFERHZo1F0A4iIyCwGdiIiyzCwExFZhoGdiMgyDOxERJZhYCcisgwDOxGRZRjYiYgsw8BORGSZ7iI+dOXKlbpmzZoiPpqIqLKazeZFVX1X2HGFBPY1a9ZgbGysiI8mIqosEZmIchxLMURElmFgJyKyDAM7EZFlGNiJiCzDwE5EZBkjgV1ErhGRR0XkZyLyvIh8wMR5iYgoPlPTHb8F4F9V9S4R6QWwzNB5K6U5MYXR8UsYWbsCQ4N9RTeHiGoqdWAXkasB/B6APwIAVZ0GMJ32vFXTnJjCPftHMT3bQm93AwfuHWFwJ6JCmCjFrAVwAcDfichPRWS/iPyagfNWyuj4JUzPttBSYGa2hdHxS0U3iYhqykRg7wawCcBfq+r7AfwvgK94DxKRHSIyJiJjFy5cMPCx5TKydgV6uxvoEqCnu4GRtSuKbhIR1ZSoaroTiKwCMKqqa+Z//l0AX1HVjwT9m+HhYbVxSQHW2IkoSyLSVNXhsONS19hV9VUROSci61X1NIBbADyX9rxVNDTYx4BORIUzNSvmTwAcmJ8RMw7gs4bOS0REMRkJ7Kr6NIDQ4QEREWWPb54SEVmGgZ2IyDIM7ERElmFgJyKyDAM7EZFlGNiJiCzDwE5EZBkGdiIiyzCwExFZhoGdiMgyDOxERJZhYCcisgwDOxGRZRjYiYgsw8BORGQZBnYiIsswsBMRWYaBnYjIMgzsRESWYWAnIrIMAzsRkWUY2ImILMPATkRkGQZ2IiLLMLATUWrNiSnsPXoGzYmpoptCALqLbgARVVtzYgr37B/F9GwLvd0NHLh3BEODfUU3q9aMZewi0iUiPxWRw6bOSUTlNzp+CdOzLbQUmJltYXT8UtFNqj2TpZgvAnje4PmIqAJG1q5Ab3cDXQL0dDcwsnZF0U2qPSOlGBHpB/ARAH8O4D4T5ySiahga7MOBe0cwOn4JI2tXdCzDNCemIh1H6ZiqsT8E4MsAlhs6HxFVyNBgX2igZi0+P6lLMSKyFcAvVbUZctwOERkTkbELFy6k/VgiKkjSGTCsxefHRMZ+E4A7ReTDAN4G4GoR+Z6qfsp9kKruA7APAIaHh9XA5xJRztJk3U4tfma2xVp8xlIHdlV9AMADACAiHwJwvzeoE5Ed/LLuqIE9Ti2e0uE8diKKLG3WHaUWT+kZDeyq+mMAPzZ5TiIqD2bd1cCMnYhiYdZdflwrhojIMrUL7FysiIhsV6tSDF+QIKI6qFXGzhckiKgOahXYuVgRS1FUPPbB7NWqFGPDVK00iyixFEWmxe2P7IP5qFVgB6o9VSvtRZHmrUEiL3d/7G4IPj68Gts29XfsU+yD+ahVKabq0j4jYCmKTHL3x+k5xfePncU9+0c7lljYB/NRu4y9yky8zl31UhSVh9Mf35ppQQEowrPwTn2Qa7WbI6r5L7Q4PDysY2NjuX+uDdj5qUyaE1M4eGISjzYnMTfXTjiS1M1Ze49GRJqqOhx2HDP2iqnyMwKyj9MfP7apP1XCwdq7WQzsRJSaO+FIMqrkWu1mMbATkTFJSyp8/mMWAzsRGZN2Iw4GdDMY2InImL5lvWiIAFCWVArEeewVwlexqcyaE1PYc/gU5lqKhgh2b92wJANnH84HM/aSCHvglOQtP6I8OWUYBaCqmLo8vejv2Yfzw4y9BJwO//UnTge+uZfkLT+iPIW9Vco+nB8G9hKIslSAc9HI/M8KYJpLD1OJODNb7rttve9sGL8+zOWzs8HAXgJR1s9wLppb3/ubC3/W0vbDKqKyGBrsw86b1/mWV5w+fPfmAa4XkzHW2DMS5yWNqHN4hwb78Durr8GR534BRfuu7K1jEpWZqTdVqTMG9gwkeUkj6hzekbUrcFUP39CjauOc9WwxsGfAxLoXQRm/iTf0uJAYkd0Y2DOQdt2LsIw/TbbDVfQoL0wgisPAnoG0WXWWK91xFT3KAxOIYqUO7CKyGsA/AFgFoAVgn6p+K+15qy5NVp3lSndcRY/yUEQCwRHCFSYy9lkAX1LVEyKyHEBTRI6o6nMGzl1L3owfAPYePWOkw3IVPcpD3gkERwiLpQ7sqvoKgFfm//sNEXkewPUAShHYq3oXdzL+LDosZyRQ1vJOIFhiXMxojV1E1gB4P4BjJs+blA138UWvYbPDUoXkmUCwxLiYscAuIu8AcBDALlX9lc/f7wCwAwAGBgZMfWxHVb2Lu0cZfct60ZrflralwBtvzhgryxDZgptkL2YksItID9pB/YCqPuZ3jKruA7APaG9mbeJzw1TxLu4dZWzb1I/26taAANj/1H+jpVrZEQhVQ9xg6Hd83gHVb4Rgw6g9CROzYgTAdwE8r6rfSN8kc/Ks85nqxN5RhgALb5oCwFxLFy2eVIdOSvmKGwz9jgdQioBa1VF7WiYy9psAfBrAsyLy9PyffVVVf2jg3KklrfPFCdRpsgLv53hHGds29WPbpn4cPDGJR8bOYW6uPdjpakglRiBUPXGDYdDqpGUIqFUctZtgYlbMU8DCSpxWiBuovQ84H3ryBezackPiG4LfKGN0/BLmWlcqWO+59mozX5bII24wDDq+DAE1bNRua/2db576iJuxOB3b+TdPvXgRx196bckNwduJgj7Hb5Sx8BkzLbQAPHv+f3DP/tHa1AwpH04f3b11A6YuT3fc0cvdl/2CZ1nelwgatdtcf2dg9xE3Y3E69kNPvoCnXrzoWwP360RxPsf9GT85c7HwIS7ZJ2qgCzrOe2zZ35ewuf7OjTZ8hO0EE/Rv7th4LboagobPBgJBnSjO5wwN9mHXlhu4SUHNZbUhdJSdvOIcV3ZRNripKmbsAeJmG+4d2rsaS3doD8rO434OlwSotyzLB1FHkJ2OC5v2CKA0fXdosA+7t27A4ydfwR0bry28PSYxsBsStkO7yYBc9iEuZSfL8kGcnbz8jgub9tjdEEAEs3PlqGk7ydj0bAvHX3oN61ctt+a6YmA3JEq2457h4v45KVuf6FOwrKfvRU0a/I7zmx22+teXXbkRzSmA8ryHYXONnYHdECeL+c6//Ry//NX/4fSrb4S+BRc286ATm5/oU7AsSnGmEgS/2WE93Q10N2ShRAkRzM2VY065zXPcGdgNOv3qGzjy3C8AAM9MPgsA2L75yro4izKamRZ2HzqZeHkAm7MN6sxkKc5kguA3O2xuroVP3jiA6655eylr7LY+r2JgDxEnm3n85CtLfnYHdneGICJoqSYOzFGyDZZqKIzpBMGZuXX8pdcWvT3tPmeZ+qKtz6sY2DuIm83csfFa/PuLFxf97ObOEPqW9WLP4VOJh4FR3qhjqYbCZFGOsC0TrmKCxMDeQdxsxsnOnelT7mzd4c4Q1q9anqrDdMo2WKqhKLIKwrZkwlVNkBjYOxhZuwJdDUFrTgERnH/9zYWXQoIuhO2bB3wDup8sO7/ND4bIrE79MEq2WsWMNqqqJkgM7B2cfvWN+Sla7eVyHz52Fo82JwFVzLbKvSa6bcNhyl+UbLWqGW1UVU2QGNg78D4Mdebfuv+7zHdwW4bDVIwo2WpVM9qoqpogca2YDrwPPwXtNSV6ugRd0l4T/WVXeYbIJlHWUslqvZWs1sNJYmiwDztvXleZoA4AoprLLnWLDA8P69jYWO6fm8T3j53F4ydfwYZrr8byt/csdNyDJybxaHMy9uvRWdcjba53Uv6KqLHbXt5JQ0SaqjocdhxLMSGCHoaOjl/C7Fy8IWiSLcfi7jvJC4JM8ivneftl0qUygvq37eWdPDCwJ5TkoUpYh/Wughc3SPOCoKyZ2t+0UxIS99riKHUpBvaEkjxUCVvu1N3Rt23qjx2kvefvW9aLvUfPsMOTMWH7m07PRNsaslMSEufa4ijVHwN7CklmnWzb1A+Z/99Ow09B/D0j/d5sZYcnk/qW9aIhAkCX7G/qbNv4kzP+W0P6nkfVt39HvbbijILr1P8Z2HPil5G7ebPtbZv6sW1Tf+xO6VwQe4+eYVmGjGlOTOGxE5N4ZOyc72YycbZtdNZBb6mi4bMpTRxxRsF1Sm4qFdiTPEzM+26d9IFQ0PBzaLBvYepXnO8RlhERReUEyLdm2hvJAEs3k/Fb/Cuoz7mvBUH7PM6NQwF8zDOa7aRT2abOz5wqE9iTzCjJ+m7tDeImHwil+R4mMyIi9+5gwJX3OfzKJ0E7K7n/zO9Z0N1/0+7jAPDo2Dk8vOMDsUepXlV9a9SEygT2uHdfE3frThm/X8BN80AoKIAn+R5+GRFRUu4A2dXVwF1D/YFZtTfIBvVr97UwOn5p4Y1uoL3Tkons2vvMydTOZVVQmcAe9+6b9m4dlin7Bdywz+w0J/j862/6BvAk3yNpGaauD5qoszSv1QfNovGeq2d+5yUA6OkSY/sLOMeVpdae1zVWmcAet3OlXeOhU6bcnJjCy6+/ubDllxM8436m++bR3dVYcr4k3yNpGabOD5ooXNJ1h/zKLn797OHPj3Sssafpn2Wpted5jRkJ7CJyO4BvAegCsF9V/9LEeb3idq6knREIzvi9wfgTN65e1BHjvIXn7nDeLcS8D1aTdOI4ZZiydH6yi1/Zxa+fhfXxNP2zLLX2PK+x1IFdRLoA7AVwK4BJAMdF5F9U9bm05y5SUKbsDcbXX/P2SDVFP35THNP+opN24rJ0frJbXv3TW/IowwqNeV5jqRcBE5EPAPgzVf2D+Z8fAABV/Yugf1OlRcC8nMDt/HK8gXvv0TP4+hOn0VKgS4D7bluPnTev63g+0x0u6TlZYyfTgpYgyLJ/lrmsmPYay3MRsOsBnHP9PAlgs4HzllLY3T/uXTlNucj0ObNoC9WbX/kh6RK4fjNuqjZ/Pa9rzERgF58/WzIMEJEdAHYAwMBAtK3jyqrTL6cswz6iMhhZuwLdXVemSppcr/3uff+BmTlFT5csmvfOsqKZwD4JYLXr534AL3sPUtV9APYB7VKMgc8trbJnviy5UK6ccq/BvR8OnpjE9Py2ldNzioMnJiO/M1IHJgL7cQC/JSLvBnAewCcBbDdwXqsEBdO8g2yZ649kn9HxS5htKRTtfYNNlUW8ZQJB8DrxdZQ6sKvqrIh8AcCP0J7u+Leqeip1yyzRnJgK3G3JO3Wy0xt9ppS5/kj2CSuLJE1stm3qxyPNyYXzbrjunZESlrqMVo3MY1fVHwL4oYlz2cRv8SR3MHUH2enZFh4+dhaPnZjMNItm/ZHy5C2LAFhY0A6I90aoNyg//Pnw+fHef1+X0Wpl3jytorDFk5wg6wR+RfZZNOuP5JV1FuuURdJsJhMUlN3HhyUsdRqtMrBnaNHiSQ3Bx4dXL3oJyQmy7nWu88ii615/pCu+f+wsdh86ibmW4qqebLPYuJvJuG84SZe9dqvTaJWBPUNROpsTZMM21ahLbZDy05yYwu5DJzHbmp9dMpM8i3X6Z9+yXkxdno70jseG694JAL7rw3gz9N1bN7R3aZptQUTQt6x3SRvCEpY6jVYZ2DPm19n8gnSnTunu5A0R7PnoRmzfXO13Aah4o+OX0HJNQWw0/FdVDON9ltQQ+Naw3YHVu3Xjxuveueia8GboU5ensXvrhoXRxZ7Dp7B+1fLYwbkuo1UG9pwleYAzOn5p4aJpqWL3oZOJOnVeOLqoBieDdicMSX5f3mdJnWrYTmB1b904PdvC7kMn0VJduCa8a8Cff/1NvPz6m2ip5vIsquqsC+xlDypJHuD0Letd9CqvyfnAptVp5kHVmSpNLNwg5jezboj/Dkt+/2ZmvrQy11ocsHfevG7R86cf/OdZdDcE3V0NzM3ZXyNPy6rAHieoFHUD8HuAE9aWqcvTaO8L39aVcMichzrNPCiKyb5rojThLbEE1dg7/Zs9h08teajplGRmW9peTbWl+MSNq3G9z9LWjrIndnmxKrBHDSqmssqkO7p45/WGtWVk7Qpc1ZN+yJyHOs08KEJZR0RJbhDuf7N+1XLfa8nbnzq9wFfW/2+KYFVgjxpUTO2HmrQTuTu0u9bYqS5Zlaf5VWprFeU9IsorAw66MYT1pzhTIuuksoE9aGZJlKBiIqs01YmitqVKT/Or1NaqyXNEVJYMOKg/BU2J5GixooG9U4eLElRMZJWmLjBmuBRHnv2liAw4zgjBb0okr6W2SgZ2Ex0ubVbpd4G5O6XTzigdzNsWPgCiTvIaEeX9vCTuCMGvfRwttlUysJflAZ27Ey1aqbEhgMiS1RyjKMvwl8qliJt9HqODNDVyjnaDVTKwl/EXuqhTzrWX9EryIgUfAJFXkTf7LDNgEzVyZuj+KhnYgfL9Qhe/cAEoBKrxF/Uqy2iEysPWmz1r5NmpbGAvStCQ2BlFOJtqOCs67t66IVYHLeNohIpVlZt93HIRa+TZYWCPIcqQ+NxrlzE7N7++uiqmLk/H/hx2bnKrws0+SbmoCt+rqhjYY/DuePTQky9g15YbFmbEeFe4S5JdRc16OHOmXsp+s3dfG2/NtBZtLt1J2b9XVTGwx+BeDa+lwFMvXsTxl15byDqcFe4aAG5at3Ih6EcVNesJO45Bn7LQqV+NrF2B7oZgeq49aeDR5mTm+/dSsEbRDagSZ+h407qVC4tyOQ+znKDfJUBvTyN2UAf8H5LFPc4J+l9/4jTu2T+K5sRUim9M1BbWr4YG+/Dx4dWQ+Z/n5oL7L2WPgT2mocE+7NpyA67qaQdx90OfA/eO4L7b1ieejua+OXQq43Q6LurNgSiOKP1q26b+JddFkObEFPYePWMs8Qg7n+nPKzuWYhIIeuiTxduscY+rygwKqpYo/Spq/w0qJSYtIUYpTdbtpT8G9oSyeujT6bzejp9kNTyiJOIkHXHXeHGy/6TBN2yev63vAXRSu8Ce9YPFrM4fJ+vgTAPKgql+5Zf9pwm+YaOJOo5iaxXYsx6SZXn+OmYdVE5+yUuUhMZ9jF/2nzT4ho0m6jiKrVVgzzo4Bp3fRBZfx6yDyscveQHCyyh+/27nzesW/j5t8A0bTdRtFJsqsIvI1wD8IYBpAD8H8FlVfd1Ew7KQZL9RE+cP6vRxPruOWQeVT1B9PCxhipJUVTH4lvWdkbQZ+xEAD6jqrIj8FYAHAPxp+mZlwxscgSuZhrOX6PbNA8bOPzTYF7j1XdJXsMvUeah+gkaOYaNJG0ecZZ5tkyqwq+oTrh9HAdyVrjnZcwdHd9BtqWL3oZNYv2p56imLUaYfsmZORYhbC49arw4bTdo44izzNWyyxv45AP9o8HypROnAI2tXoCGClioAoNVS47+coA5tYwZD5RYlw4xyjN/IMcpo0rYRZ5mv4dDALiJPAljl81cPquqh+WMeBDAL4ECH8+wAsAMABgaSlzuiiDpEGhrsw56PbsTuQyfRail6e7L55QRdCLZlMFRuUTLMJFloWevMWSvzNRwa2FV1S6e/F5HPANgK4BbV+dTX/zz7AOwDgOHh4cDjTIjTObdvHsD6VctT/XKSdmzbMhgqtygZZtws1GSduYo3iLJew2lnxdyO9sPSD6rqZTNNSi9u50zzyynzAxQitygZZtws1FSdmdeRWWlr7N8GcBWAIyICAKOq+sepW5VSnkOkpB27itkJVZ/pWripOnOZH0RWUdpZMevCjypGXkOkJB2b2QmVXdTEw1QS1besFw0RIME+wbRUrd48zUKSjs3shMosbuKRNolqTkxhz+FTaKmikWCfYFqKgd2AuB27zNOkiLJOPLyjAffnCZLtExz1s+qCgb2DrDqF3xuwe4+eqV3no3LKMvHwGw1k9Xl1LnkysAfIulM4WX6dOx+VU5aTD/xGAztvXpfJ59W55MnAHiCvTlHnzkflldXkg6Ds3Pt5XBE1ndoE9rgdJa9OUefOR/UTZTRgahQbdeRhYx2+FoE96UqKecyFL/NryURZCBsNmBzFhn2WraXQWgT2pB0lr7nwZX0tmagIeY5ibS2F1iKws9xBVB15jmJtjQ3SYd2uzAwPD+vY2Fiun2ljHY0oD7ZfO1X6fiLSVNXhsONqkbEDLHcQJZGmBl2VgGljbKhNYE+rKp2UyKQ0i9zZ+FCyKhjYI2AnpbpKWoO29aFkVTCwR8BOSnWV9EGmrQ8lq4KBPQJ2UqqzJDVovp9RrNrMikmLNXYiKhpnxRhm45NzIi8mMHZgYK8IXnCUNU4SsAcDe0J5BlpecJQHThKwBwN7TM2JKTx2YhKPjJ3DbEtzCbS84CgPnCRgDwb2GJzM+a2ZFpxHznkEWl5wlAfOZLEHA3sMTubsBHUBcgm0vOAoL5wkYAcG9hjcmXNXVwN3DfXjY5v6ubQvEZUKA3sMUXd/YWZNREViYI+pU+bM2StEVAaNohtgE7/ZK0REeTMS2EXkfhFREVlp4nxV5dTguySfh6pERH5Sl2JEZDWAWwGcTd+cauPsFSIqAxM19m8C+DKAQwbOVXmcvUJ1xskD5ZAqsIvInQDOq+ozImKoSURURZw8UB6hgV1EngSwyuevHgTwVQC3RfkgEdkBYAcADAwMxGgiEVUBl74oj9DArqpb/P5cRN4H4N0AnGy9H8AJEblRVV/1Oc8+APuA9nrsaRpNROXDpS/KI3EpRlWfBfAbzs8i8hKAYVW9aKBdRFQxnDxQHnxBiYiMKXryAB/ethkL7Kq6xtS5iIji4sPbK/jmKRFZgW9+X8HATkRW4JvfV7DGHgPrd0TlxYe3VzCwR8T6HVH5Ff3wtixYiomI9TsiqgoG9ohYvyOiqmApJiLW74ioKhjYY2D9joiqgKWYkmpOTGHv0TNoTkwV3RQiqhhm7CXEGThElAYz9hLiDBwiSoOBvYQ4A4eI0mAppoQ4A4eI0mBgLynOwCGipFiKISKyDAM7EXXEqbfVw1IMEQXi1NtqYsZORIE49baaGNiJKBCn3lYTSzFEFIhTb6uJgZ2IOuLU2+phKYaIyDIM7ERElmFgJyKyDAM7EZFlGNiJiCzDwE5EZBlR1fw/VOQCgIncP9jfSgAXi25EQer63fm968Wm7z2oqu8KO6iQwF4mIjKmqsNFt6MIdf3u/N71UsfvzVIMEZFlGNiJiCzDwA7sK7oBBarrd+f3rpfafe/a19iJiGzDjJ2IyDIM7C4icr+IqIisLLoteRCRr4nIz0Tkv0Tkn0XkmqLblCURuV1ETovIGRH5StHtyYOIrBaRoyLyvIicEpEvFt2mPIlIl4j8VEQOF92WPDGwzxOR1QBuBXC26Lbk6AiAjar62wBeAPBAwe3JjIh0AdgL4A4A7wVwt4i8t9hW5WIWwJdU9T0ARgDsrMn3dnwRwPNFNyJvDOxXfBPAlwHU5qGDqj6hqrPzP44C6C+yPRm7EcAZVR1X1WkAPwDw0YLblDlVfUVVT8z/9xtoB7nri21VPkSkH8BHAOwvui15Y2AHICJ3Ajivqs8U3ZYCfQ7A40U3IkPXAzjn+nkSNQlwDhFZA+D9AI4V25LcPIR2stYquiF5q80OSiLyJIBVPn/1IICvArgt3xblo9P3VtVD88c8iPaQ/UCebcuZ+PxZbUZnIvIOAAcB7FLVXxXdnqyJyFYAv1TVpoh8qOj25K02gV1Vt/j9uYi8D8C7ATwjIkC7HHFCRG5U1VdzbGImgr63Q0Q+A2ArgFvU7rmvkwBWu37uB/ByQW3JlYj0oB3UD6jqY0W3Jyc3AbhTRD4M4G0ArhaR76nqpwpuVy44j91DRF4CMKyqtiwaFEhEbgfwDQAfVNULRbcnSyLSjfYD4lsAnAdwHMB2VT1VaMMyJu1s5e8BvKaqu4puTxHmM/b7VXVr0W3JC2vs9fZtAMsBHBGRp0XkO0U3KCvzD4m/AOBHaD9A/Cfbg/q8mwB8GsDvz/+On57PYslizNiJiCzDjJ2IyDIM7ERElmFgJyKyDAM7EZFlGNiJiCzDwE5EZBkGdiIiyzCwExFZ5v8B48ZTDN1j3A8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(samples_np_proj[:,0], samples_np_proj[:,1], '.')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 157,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHFdJREFUeJzt3W2MHVd5B/D/s3dthxZQpE2qSHFcg1pVRXUVl1Xaq0hlqU2UQpp8CB8ognUTiSWSU2EJRFmiSJEsxS1IxEhBrZdgC4tItCK8CdGS2GWlVLm8rEOoBYGKRmCcgDCuEKiV7az36Yezk52dnfc5M3POmf9PWq139+7cc71nnvuc55w5I6oKIiIKx1TfDSAiIrsY2ImIAsPATkQUGAZ2IqLAMLATEQWGgZ2IKDAM7EREgWFgJyIKDAM7EVFgpvt40uuuu053797dx1MTEXnrzJkzv1TV64se10tg3717N1ZWVvp4aiIib4nIT8o8jqUYIqLAMLATEQWGgZ2IKDAM7EREgWFgJyIKjJXALiLXisjnROQHIvK8iIxtHJeIiKqztdzx4wD+TVXfLiLbAfyWpeN6ZTIBlpeBuTlgzLc2IupJ48AuIq8F8OcA/gYAVPUKgCtNj+ubyQTYtw+4cgXYvh04fZrBnYj6YaMU83oAFwCcEJHviMhjIvLbFo7rleVlE9SvXjWfl5f7bhERDZWNwD4N4E8A/KOq7gXwvwA+lHyQiCyIyIqIrFy4cMHC07plbs5k6qOR+Tw313eLiGiobAT28wDOq+o317/+HEyg30RVl1R1VlVnr7++cKsD74zHpvxy+DDLMETUr8Y1dlX9uYj8VET+QFV/CGAfgO83b5p/xmMGdCLqn61VMX8L4PH1FTEvALjH0nGJiKgiK4FdVZ8DMGvjWERE1AyvPCUiCgwDOxFRYBjYiYgCw8BORBQYBnYiosAwsBMRBYaBnYgoMAzsRESBYWAnIgoMAzsRUWAY2ImIAsPATkQUGAZ2IqLAMLATEQWGgZ2IKDAM7EREgWFgJyIKDAM7EVFgGNiJiALDwE5EFBgGdiKiwDCwExEFhoGdiCgwDOxERIFhYCei5iYT4MgR85l6N913A4jIc5MJsG8fcOUKsH07cPo0MB733apBs5axi8hIRL4jIl+xdUwi8sDysgnqV6+az8vLfbdo8GyWYt4H4HmLxyMiH8zNmUx9NDKf5+b6btHgWQnsIrITwNsAPGbjeETkkfHYlF8OHy4uw7AW3wlbNfajAD4I4DWWjkdEPhmPi+vqrMV3pnHGLiJ3APiFqp4peNyCiKyIyMqFCxeaPi0R9aVu1s1afGdsZOy3ArhTRN4K4BoArxWRz6jqu+IPUtUlAEsAMDs7qxael4i61iTrjmrx0e+yFt+axhm7qi6q6k5V3Q3gHQD+PRnUiSgQTbLuKrV4aoTr2ImovKZZd5laPDVmNbCr6jKAZZvHJCKHRFn38rIJ6gzSTmLGTkTVMOt2HveKISIKzOACO6+PIKLQDaoUw+sjiGgIBpWx8/oIIhqCQQV27lXEUhQ5gJ2wdYMqxYSwUmsyqd9+lqLIuqodkp2wE4MK7IDfK7WanhNppShf/y/IAfEOORoB994LzM/ndyp2wk4MqhTju6ZzBCxFkVXJDnnsmAn0eSUWdsJODC5j95mNq7l9L0WRQ6IOeekSoGo+irLwvE7YpM5Im4hq9xstzs7O6srKSufPGwL2fXLKZAKcPAmcOAGsrtavm7P2XoqInFHV2aLHMWP3jM9zBBSgqEPOzzfLOFh7t4qBnYiai2ccdYaV3KvdKgZ2IrKnbkmFE0BWMbATkT1NSiqsM1rDwE5E9szMAFNTZoUMSyq94Tp2j/BKbHLaZAIcOmSy9akp4OjRrRk4O3EnmLE7omi+qc5FfkSdisowa2uACHDx4uafsxN3hhm7A6L+/uCD2Rfu1bnIj6hTRVeVshN3hoHdAWW2CojOGRHztSpw+TK3HiaHRCtbDh9OXw2T1om5f3YrGNgdUGb7jOicueuuje+trZm5KiJnjMfA4mJ6eSXqxO99L7BjB/eLaRFr7C2pco1G2SW84zFwyy3Al79sgvrU1NYyJpHTbF2pSrkY2FtQ5xqNskt45+ZMssML9MhrXLPeKpZiWmDjFnxZq8KKyphNjk1EYWDG3oKm214UZfxNkh1uoked4VakvWFgb0HTbS/a3OiOm+hRJ5hB9KpxYBeRmwCcBHADgDUAS6r68abH9V2TrLrNje64iR51oo8MgiOEV9jI2FcBvF9VnxWR1wA4IyJPqer3LRx7kJIZP2Bq4jb6KzfRo050nUFwhLBJ48Cuqj8D8LP1f/9GRJ4HcCMAJwK7r2/iUcbfRn/lggRqXdcZBGuMm1itsYvIbgB7AXzT5nHrCuFNfHnZXGG6trZxpalvr4EGqssMgjXGTawtdxSRVwN4AsAhVf11ys8XRGRFRFYuXLhg62lz2Vh22If4csSZGRPUAfP5V7/iUkWiLfLWAQ9wfa+VjF1EtsEE9cdV9fNpj1HVJQBLgLmZtY3nLeLjm3hylHHggLnCNNow75FHzL99HYGQJ6rWMNMe33UdNG2EEMKwvQYbq2IEwKcAPK+qH2veJHu6LPPZ6sPJUQawcaUpYG4EH987aQB9lLpWNRimPR5wI6AOtPZuI2O/FcC7AZwVkefWv/dhVf2qhWM3VrfMVyVQN0kKks+THGXMz5uPkyeBT33K9E/A7J/kwwiEPFQ1GGbVPF0IqD4O2y2wsSrmPwCIhbY4o2qgjvfry5eBhx4yH3XfENJGGcvLG7V2ALj55iavkChH1WCY9XgXAmrRsN3XZXMFeOVpiqoJS9Svo9Urp04BTz+dPocT70NZz5M2ykg+x8qKeVMYSMmQuhJ10qNHzdahebf0infmtODpygUTWcP2kOvvqtr5xxvf+EZ12TPPqL7qVaqjkfn8zDPlfue221SnplQB87sPP5x/zKrPU/QcRI2U7ZB1ThAXPfyweQ0enUwAVrREjOXujinq7KA4HgN3321q31NTW0efWdl5lecZj02Jh/coGLi2lu+VXR/s6zripDJ3uPEUSzEZqk66Ft2gPasMWfV5uCXAwLVZPihbW897XNGyR8CdzjsemxP1iSdMVtZ3eyxiYLckSmKybtBuMyBzS4ABa3P5XpVbeaU9rmjZ42hkTo7VVTdq2lE2duWKmRTbsyeYE4uB3ZIyyU58hUv867oCndCnPG0v3yubNaQ9Lm152Otfv/G9aFmXKxdiBLzGnYHdkiiJ+chHgJdeAs6eLb4IrmjhQZ6QJ/QpRxu1OFsZQtrysG3bTKYObM3Y+65pB7zGnYHdorNngS9+0fz7W98ynxcWNn6eTGjuv7/+9gABJxtUxGYtzmaGEL3pPPSQCeprayaIv+c9wK5dbtbYA52wYmAvUCWZeeKJrV/HA3s8QRDZGJ3WCcxlkg2WaqiQ7QwhWrr19NObL5+OH9OlzhjohBUDe46qyczddwNPPrn567h4gjAzszFvU2cUWOaCOpZqqFAb5YjQMmEPMyQG9hxVk5koO49WT8Wz9Ug8Qdizp1l/yUs2WKqhUtoKwqFkwp5mSAzsOebmzJr0q1dN6eTcuY1rQrLOg4WF9ICeps2+H/C8ENmW1xHLZKseZrSleZohMbDnOHsWePll8+/VVeDYMeDECbNa6+pVt9/AQxsNUw/KZKueZrSleZohMbDnSE6GRstv4/92+Q08lNEw9aRMtuppRluapxkS94rJkZz8FDFv2tHS3NFoc3mGKChl9lJpa78Vl25nNx4Di4veBHUAELNhWLdmZ2d1ZWWl8+etY2nJZO433wxce+1Gvz150pRlql4d3XY5MuRyJ/Wgjxp76OWdBkTkjKrOFj2OpZgCWZOhy8smqFcZgda541iV84XnA1mXVs9Ldsy6e2VkdfDQyzsdYGCvqc6cSlF/TW6CVzVI83yg1tm6v2leFlL15OIwdQsG9prqzKkU7XYa7+cHDlQP0snjz8yYMiX7O1mTlj0A1e8NmZeFVDm5OExNxcDeQJ1VJwcOmM/Jq6yT/RyoPiLIu7KV/Z2smJkxF3eobr2/adG9IZPHEUm/Kw1Q/uSqMgwe0glQ5jZLtj9cvzVeG4ruJpZ167yHH6535zEP7/pFLnvmGdX77lPdvt3cm3F6WvXYsc0/L3vfxqizT02pbtu2+Th12pV1YoVyC78YlLw1nlcZe53JxK7frOvOB2WNPsfjjZVfVV5HUUJEVFpU7rh0yWTqwNa7yaRt/pXV6aKTIX5XmsnELDUDtg5n8+SVbQY86eRNYK+zoqTt0lsyiNucD2ryOqIbw6ytmeXFydv0EVUSBch4UM8qn2TdWSn+vbTJoDe/2ZRyAOD48WpBOKts4+lVozZ4E9irvvnaeLPOy/jTAm6T+aCsAF7ndaQlRES1xQPk9DRwzz3ZWXUyyGZ17PjJEHXYyMsv28muk5NOVZdjesybwF71zbfpm3VRppwWcIueM29J8Llz6QG8zuuoW4YZ6jwTFWhyWX3WKprksaLJV8Bc2m3rBgPR41xZOdPVSVamEG/7o+7kadXJxLYmH+PzSMl5mSrPGZ/b2bEj/Xh1j1l1XirAeSZyQbJjHTuW3tGik+q++9I7X5MO6spKAgsnGbqcPBWR2wF8HMAIwGOq+vc2jptUdXlhk02wsjLleCY/PW3u+hUflVa5CC+ezACb7yCWnFgt+zrqlmEGPM9EbcoquyQ7WlEnb9JBXam1d3iSNQ7sIjIC8AkAbwFwHsC3ReTLqvr9psfuU9boMxmMd+0qV1JMk+xvVRYDZKnbh13p+xS4rjposuThwg6NHZ5kNjL2WwD8SFVfAAAR+SyAuwB4HdiB9CSi6G9T5U25jf5W95iu9H0KTFqm03YHzcqu+u7UHZ5kNgL7jQB+Gvv6PIA/tXBcJxX9baq+KbfR3+oe04W+T4FJy3TqboGbtuLGt/XrHZ1kNgK7pHxvy17AIrIAYAEAdu3aZeFp+5P3t2HmSxQzN2cmo9bWzGeb+7XPzZmlkdu2bQ7erCtaCeznAdwU+3ongJeSD1LVJQBLgNmP3cLzOsv1zJfLGqlT0YVNNu/9cPLkxtr3K1fM12UvGhkAG4H92wB+X0ReB+BFAO8A8E4Lxw1KVjDtOshyMzzq1PKyKYlENwpusyyStU/8ADUO7Kq6KiL3A/gazHLH46r6vcYtC0S0BUba3ZaSSyfzLuizxeXyIwWoqCxSN7OZnzcnVXTcvXvLZSwDGa5aWceuql8F8FUbxwpJ2t5J8WAaD7JXrwLHjgGf/nS7WTTLj9SpZFkE2NjRDmh2S7Gvf714fXzy9wcyXPVmSwEfRX0ta++kKMhGgV+1/Sya5Ufaou0sNiqLNLmbTJkljEUZy4CGqwzsLYpnx6MRcO+9W69SPX3alGqOHzf9rYsseuDlR4pbWgIOHjSrVnbsaDeLTQZWoHyZpigol8lYBjRcZWBvUZm+FgXZ+fn8xw2kNEhdmkyA++83kz+A2YSrbhYbddCZGbOPRZmLPPbu3fhZcnIpmaEfPbqxUdjUlHmepKKMZUDDVQb2luXt6Ji8oUbRKPTyZZP5P/oosLDQZqtpEKIVK5GpqXpZbLyDrq2Z46Rl//HAmrx34969m0+KZIZ+8aIJ7gcPmu8dOgTs2VM9OA9kuMrA3rE68zfLyxt1+LU107fr9OmucHThibk5E4DjGUOdP1h85znAfM6qYUeB9ciRzTfAPnhw4z6qp09v3QP+3DnzEZ0EgdfImwousLseVOrM38zMbL62o+3lwE0MaOGB/2yVJqIgHM/Yi2rY8cA9NWU6dTxgLy5unoD65CfNm8/0esgKvEbeVFCBvUpQ6esNIG3+pqgtFy+aFTVRcB+N3O3TA1p40B+bnddGaSJZYsmqsef9TrwsE3XuqCQTrQcGsve2jrie2XWlzKbttj/q3mijSNn99G3dVKLujTziv1emLTZv6t423rCjZaH+B2edTFVeb6j/NzHo8kYbrii7mslGVtmk5BBPkuKlxryypC+T+T611UtdD4m6yoCzRg5FHarKksgB8TawZ60sKRNUbCxntdWHyrbFp8l8n9rqnS7XYrsyYZLVobKWRA5gnXoRLwN7Xn8rE1RsZJW2zi9muFRJlx2mjwy4ygghbUkkTyYAngZ2G/2taVaZdn7F+2TUzjL9K9kWzv9Qrq6GRF1fqVl1hJDWPg4XAXga2F25Mjjeh+J9cjQyq1iSuzmW4crolxzTx7t9F6ODJjVyDnczeRnYXfx7xvtkdJ1GnU29OP9DW/T5bt9mBmyjRs4MPZWXgR1w7+8ZH0VEa86jC+mqjChcGY2QQ0J9t2eNvDXeBva+ZI2Io1FEdFONl182F9QdPVqtf7o4GqGe+fJuX7VcxBp5axjYKygzIn7hBRPU19ZM5n7xYvXnYd+mTXx4t69TLvLhdXmKgb2C+Mjx8mXgoYfMR7QiJrnBXZ3kqmzSw5UzA+P6u3385Lh0afPNpfO4/ro8xcBeQXKvo1OngKef3kg6og3upqaA/fs3gn5ZZZOeoscx6FMr8jrW3JxZDhbduPrEifZv4EuZpvpugE+ikeP+/SZ4xzeji4L+aGR2Qq0a1IH0ObKqj4uC/oMPms+TSY0XSpRU1LHGY3OLMBHz9epqdgem1jGwVzQem6C9Y4cJ4vE5n9OngcOH669Gi7855JVx8h5X9s2BqJIyHWt+HrjmmuIODJg3hiNH7GUeRcez/XyuK7NTmO2PtnZ37FLdnR1tHdfGZnhEpZXtWGU6cNaxmmyXmte2gE4KDHF3xy61NeeTd9xkibPOZnhEtZTtWGVOjKzsv+5FWEXr/EO9DiDH4AJ72xOLbR2/ymoyLjSgVtjqWGnr15sE36J1/r5cB2DRoAJ721dmt3n8ASYd5Kq07KVMRhN/TFr2Xzf4Fo0mBjiMHVRgbzs4Zh3fRhY/wKSDXJSWvQDFGU3a7y0ubvy8afAtGk0MbBjbKLCLyEcB/BWAKwD+G8A9qvorGw1rQ537jdo4flafr/LcA0w6yEVZ9fGijKlMVuVj8HX0opGmGftTABZVdVVE/gHAIoC/a96sdiSDI7BxtehoBDz6KLCwYO/443H2re/qXoHtUN+hIcoaOhYNJ0Mccjq8x3ajwK6qT8a+/AaAtzdrTvviwfHIkY2rSNfWgIMHgT17mv1tksE3qz+zZk69qFoLL1uvLhpOhjjkdPgktlljvxfAP1s8XiNl+m90FXS0f/ramv2/TVZ/DjGBIceVyTDLPCZt6FhmOBnakNPhk7gwsIvIKQA3pPzoAVX90vpjHgCwCuDxnOMsAFgAgF27dtVqbFllR0jjsSm/HDxogvqOHe38bbLOg9ASGHJcmQyzThbqaJ25dQ6fxIWBXVX35/1cRA4AuAPAvvUro7KOswRgCQBmZ2czH2dDlb65sGDKL03+NnX7dWgJDDmuTIZZNQu1WWf28Q3C0ZO46aqY22EmS9+kqv9np0nNVe2bTf42Ds+fEG1WJsOsmoXaqjPzRLKqaY39UQA7ADwlZle3b6jqfY1b1VCXI6S6/drH5IQCYLsWbqvO7PBEpI+aror5PVsNsa2rEVKdfs3khJxXNvOwlUXNzJgtf+veoYY2GdSVp22o06+ZnJDTqmYeTbOoyQQ4dMisYBiNqt8omLZgYLegar92eJUUUfuZR3I0ED1fkxsFl32ugWBgz9FWn0i7AvbIkcH1PXJVm5lH2migrecbcM2TgT1D230iyvIH3PfIVW2uPkgbDSwutvN8A655MrBn6KpPDLjvkcvaWn2QlZ0nn49bojYymMBetZ901ScG3PdoiMqMBmwNY8uOPAKsww8isNfdSbGLtfAOX5VM1I6i0YDNYWzRcwVaCx1EYK/bT7paC+/oVclE/ehyGBtoLXQQgZ3lDiKPdDmMDTQ4SM6+Xa2ZnZ3VlZWVTp8zwDIaUTdCP3k8en0ickZVZwsfN5TATkQ1NKlBexQwfVE2sA+iFGMD+ygNUpNd7gKclPQFA3sJ7KM0WHVr0IFOSvqCgb0E9lEarLoTmYFOSvqCgb0E9lEatDrrcXmBRq8Y2EtgHyWqgRdo9IaBvST2URoErhIIAgO7J3i+Ueu4SiAYDOw1dRloeb5RJ7hKIBgM7BVNJsDJk8Dx46b/dxFoeb5RJ7hKIBgM7BVEmfOlS0B0wW4XgZbnG3WCqwSCwcBeQZQ5R0FdpJtAy/ONOsNVAkFgYK8gnjlPTwP33APMz3NrXyJyCwN7BWVv/sLMmoj6xMBeUV7mzNUrROSCqb4bEJK01StERF2zEthF5AMioiJynY3j+SqqwY9GXL1CRP1pXIoRkZsAvAXAuebN8RtXrxCRC2zU2B8B8EEAX7JwLO9x9QoNGlcPOKFRYBeROwG8qKrfFRFLTSIiL3H1gDMKA7uInAJwQ8qPHgDwYQC3lXkiEVkAsAAAu3btqtBEIvIC975wRmFgV9X9ad8XkT0AXgcgytZ3AnhWRG5R1Z+nHGcJwBJgbmbdpNFE5CDufeGM2qUYVT0L4Heir0XkxwBmVfWXFtpFRL7h6gFn8AIlIrKn79UDnLwFYDGwq+puW8ciIqqMk7ev4JWnRBQGXvr9CgZ2IgoDL/1+BWvsFbB8R+QwTt6+goG9JJbviDzQ9+StI1iKKYnlOyLyBQN7SSzfEZEvWIopieU7IvIFA3sFLN8RkQ9YinHUZAIcOWI+ExFVwYzdQVyBQ0RNMGN3EFfgEFETDOwO4gocImqCpRgHcQUOETXBwO4orsAhorpYiiEiCgwDOxHl49pb77AUQ0TZuPbWS8zYiSgb1956iYGdiLJx7a2XWIohomxce+slBnYiyse1t95hKYaIKDAM7EREgWFgJyIKDAM7EVFgGNiJiALDwE5EFBhR1e6fVOQCgJ90/sTprgPwy74b0ZOhvna+7mEJ6XX/rqpeX/SgXgK7S0RkRVVn+25HH4b62vm6h2WIr5ulGCKiwDCwExEFhoEdWOq7AT0a6mvn6x6Wwb3uwdfYiYhCw4ydiCgwDOwxIvIBEVERua7vtnRBRD4qIj8Qkf8UkS+IyLV9t6lNInK7iPxQRH4kIh/quz1dEJGbROTrIvK8iHxPRN7Xd5u6JCIjEfmOiHyl77Z0iYF9nYjcBOAtAM713ZYOPQXgj1T1jwH8F4DFntvTGhEZAfgEgL8E8AYAfy0ib+i3VZ1YBfB+Vf1DAH8G4OBAXnfkfQCe77sRXWNg3/AIgA8CGMykg6o+qaqr619+A8DOPtvTslsA/EhVX1DVKwA+C+CuntvUOlX9mao+u/7v38AEuRv7bVU3RGQngLcBeKzvtnSNgR2AiNwJ4EVV/W7fbenRvQD+te9GtOhGAD+NfX0eAwlwERHZDWAvgG/225LOHIVJ1tb6bkjXBnMHJRE5BeCGlB89AODDAG7rtkXdyHvdqvql9cc8ADNkf7zLtnVMUr43mNGZiLwawBMADqnqr/tuT9tE5A4Av1DVMyIy13d7ujaYwK6q+9O+LyJ7ALwOwHdFBDDliGdF5BZV/XmHTWxF1uuOiMgBAHcA2Kdhr309D+Cm2Nc7AbzUU1s6JSLbYIL646r6+b7b05FbAdwpIm8FcA2A14rIZ1T1XT23qxNcx54gIj8GMKuqoWwalElEbgfwMQBvUtULfbenTSIyDTNBvA/AiwC+DeCdqvq9XhvWMjHZyqcB/I+qHuq7PX1Yz9g/oKp39N2WrrDGPmyPAngNgKdE5DkR+ae+G9SW9Uni+wF8DWYC8V9CD+rrbgXwbgB/sf43fm49i6WAMWMnIgoMM3YiosAwsBMRBYaBnYgoMAzsRESBYWAnIgoMAzsRUWAY2ImIAsPATkQUmP8HcCgrRUBmRN8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(samples_np_proj[:100,0], samples_np_proj[:100,1], 'b.')\n",
"plt.plot(samples_np_proj[100:,0], samples_np_proj[100:,1], 'r.')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"PCA = Principal Component Analysis\n",
"\n",
"Population stratification"
]
},
{
"cell_type": "code",
"execution_count": 189,
"metadata": {},
"outputs": [],
"source": [
"def my_pca(N1, N2, SNPs, pert):\n",
" samples_genotypes = generate_two_populations(N1,N2,SNPs, pert)\n",
" samples_numerical = numerical_modeling(samples_genotypes)\n",
" samples_np = np.array(samples_numerical)\n",
" \n",
" pca = PCA(n_components=2)\n",
" pca.fit(samples_np.T)\n",
" samples_np_proj = pca.transform(samples_np.T)\n",
" #plt.plot(samples_np_proj[:N1,0], samples_np_proj[:N1,1], 'b.')\n",
" #plt.plot(samples_np_proj[N1:,0], samples_np_proj[N1:,1], 'r.')\n",
" plt.plot(samples_np_proj[:,0], samples_np_proj[:,1], '.', c=\"black\")\n",
" plt.show()\n",
" return samples_np_proj"
]
},
{
"cell_type": "code",
"execution_count": 216,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnX+MpEeZ379Pd8+MnbsgS4ORDy+TJQJZOBhsM2oxQlo2GAbDgZnL5qSLDLvxIpax1siru1NzjsUxJ8TOsURhpRhxbR3mdhTfRSftAQkOssHxyBv13MHYZ7CJ4QKnYAi3wTcRISTs7M505Y+ZaldXV9Vb9b71/ui3n4/Usrfnfeutet+3v/XUU089RUIIMAzDMPWhUXYFGIZhmLiwsDMMw9QMFnaGYZiawcLOMAxTM1jYGYZhagYLO8MwTM1gYWcYhqkZLOwMwzA1g4WdYRimZrTKuOjLX/5ycfDgwTIuzTAMM7Y89dRTfy+EuDbpuFKE/eDBg9jc3Czj0gzDMGMLEf3Q5zh2xTAMw9QMFnaGYZiawcLOMAxTM1jYGYZhagYLO8MwTM1gYWcYhqkZ0YSdiJpE9NdE9JVYZTIMU202NjawurqKjY2NsqvCKMSMY78XwPMAXhaxTIapDRsbG1hfX8fhw4exsLBQdnUys7Gxgdtuuw2XL1/G9PQ0Hn/88Vq0qw5EEXYiOgDg1wF8EsBvxyiTYepEHUVwfX0dly9fxu7uLi5fvoz19fWxb1NdiOWKOQugA6AfqTyGqRUmERx3Dh8+jOnpaTSbTUxPT+Pw4cNlV4nZJ7PFTkTvAfBTIcRTRHTYcdwJACcAYG5uLutlGWaskCIoLfY6iODCwgIef/zxWrmX6gIJIbIVQLQK4AMAdgBchT0f+18IId5vO2d+fl5wrhhm0qibj50pHiJ6Sggxn3hcVmHXLnoYwO8KId7jOo6FnWEYJhxfYec4doZhmJoRNW2vEGIdwHrMMhmGYZgw2GJnGIapGSzsDMMwNYOFnWEYpmawsDMMw9QMFnaGYZiawcLOMAxTM1jYGYZhagYLew3gnNgMw6hEXaDEFE8d08EyDJMNttjHnDqmg2WKhUd89YMt9jGnjulgmZfIOyMkj/jqCQv7mMM5setLEaLLuyDVExb2GrCwsMA/xhpShOjyiK+esLAzTEUpQnR5xFdPom604QtvtMEwfvCuS4yK70YbbLEzTIVhNxuThszhjkR0FRF9g4i+RUTfIaI/iFExJgwOWWMYRhLDYt8G8DYhxC+IaArAfyGirwoh/jJC2YwHHLLGMIxKZotd7PGL/X9O7X+Kd9xPMHr0xNraGlvvDDPBRPGxE1ETwFMAXgPgs0KIv4pRLuOHGj3RbDbxhS98ATs7O2y9M8yEEiWlgBBiVwhxM4ADANpE9Hr9GCI6QUSbRLT54osvxrgss48MWfvEJz6B48ePY2dnh1MMMMwEEz3ckYg+DuD/CiH+te0YDnfMD/a3M0x9KSzckYiuBXBFCPEzIroawNsBfCpruUw6eMHJ+MIx60wsYvjYfw3AuX0/ewPAnwshvhKhXCYlHPs8fvBIi4lJZmEXQnwbwC0R6sIwlcVlTcewtDkZFxMTXnnKMAm4rOlYljYn42JiwhttMEwCrs1MYm10okY2sRuGyQpb7AyTgMuajmlp89xIfkzaxDQLO8Mk4Io08o1CqqqwVLVeMZnEiWkWdqYU0ghKLBFKU47Lmk6ytKsqLFWtV2wmcWKahZ2JTpJwphGUWCJUhphVVViqWq/YTOLENE+eVoA6pdyVwvmxj30Mt912m7FNaSYcY01SxionBCkszWbTKixlvAM+9aoDkzgxzRZ7ydRtOOxjBaaxoGJZXWVYb0l++LLegSJXKZfty5+0iWkW9pKp23DYRzjTCEosEcpLzJKEyyUsMd+BUAEtQvDqZryMBUKIwj9vetObBLNHr9cTV199tWg2m+Lqq68WvV6vkGuePn06t2vlXX7VyPoMY70DZbxLPpw+fVo0m00BQDSbTXH69OmyqzS2ANgUHhrLFnvJuCxIX+srxEorwnqatGGvanFfunQJa2trQe2PNYqo6uhvEicvS8dH/WN/2GJPxtf6CrXS0lhPk2aBh9Lr9cT09LTA3s5hYmZmJtq9Crn3VbXYheB3KBZgi3288bW+Qq20UOuJ/aPJLCws4Pjx4+h2uxBCYGdnJ7O1vLGxgbW1NTz00EPY3d31uvdVTtk8aaO4smFhryi+Ahwq1KE/frXj2N7exsrKClZWVvhHqnH06FGcO3cuirtBdqaXLl2C2N8Ix9Vp6664qj6bsiNjJgofsz72Z9JdMb7D0tjHpUEO7xuNhgAgGo1G5Yb5ZaHf91jPQXWXARBENLjn6jV6vZ5YXl4W09PThblf0raxym6icQKerpjMIg3gVQCeAPA8gO8AuDfpnEkW9nF8wXu9nlhcXByIO0c25Psc1bJnZmbE8vLyQMjl99PT02JmZkYQ0aADyPu5ZGkzR8bEwVfYY6w83QHwO0KI1wF4M4CTRHRjhHJrSRkrH7OysLCAlZUVzMzM1H6Voi95Pkd1peQTTzyBz33uc1hYWBi65pUrV3D58uWBq4aIvOdM0q5wzdLmSVnlWhl81D/kA+DLAN7hOoYtdrvVU+XoAZ+6JR1j+nuV22yjjOeoXrPVaolmsykajcaQVZ+lziHXT3v+uD3nqoGiXDFDhQEHAbwA4GWu4yZZ2IWwv+Dj6KZR8RE7/e/j3Gbp49ZFNY82yXem2+2K5eVlMTMzIxqNhpiamhLdbterjBjuEBbncvEV9mhRMUT0qwDOAzglhPi54e8nAJwAgLm5uViXHUtskQtVXWDiS1L9bUN5nzbnvedoWmQkzLlz5wbhiLGfox5y+s53vnPghiEibG1teZUj3SHb29sgIszOzgbXJUbUDUfHFICP+id9AEwBeBTAb/scP+kWu42qWa+h1lleFrvrmKS/5Wld2izg2M9RvU6j0RiKmGk2m94WuxBCdLtdMTU1VVp0k1zMRURienq69Hd83ECBUTEEYA3AWd9zWNjtVGWom1ac8vCxu1wIRYmrrS1FdCq9Xm8QASP96lLYAQS1r+zolOXl5aG6Ly8vF3r9ccdX2GO4Yt4C4AMAniWiZ/a/+1dCiP8UoezKE3tYWZWhbhp3gs91Te1LanOaPUeLcGslbZkX83p7v+m96JdWqzUUERPSPs7bMiH4qH/sT5Ut9nHKzWGzftPWSV/8ElJO3vfC9Vxi34eysLVRt7LlhG3ahUndblcsLi6Kbrdb+AhRHX3EzKkzKaCMqBjfT1WFPVQMyhrWdrtd0W63xdTU1Ehd09ZJtr3RaIhWqxX8o19eXh4slqnKApSsolWk6NnePRl5YxLxNPVTrzMzMxN11apvfaribhxHWNhTECqKvhN/MV/ibrc75KPU65rWUj19+vSQ73ZqairIrz4zMzM4tw6TYkVb/KZ3Txdhn1j1kOsQ0UhnnPZ9HccR0jjiK+ycBEwhdkKtPDIjnj9/fuQ7ta5pM/wdPnwYjUYD/X4fALC7u+ud8Gt9fR07OzsA9nzAx48fH4swNtecgM1Hn1eonundU+sA7IUJZ72mep1Wq4V+v4+dnR20Wi3Mzs6mfl/HPVS3dviof+xPVS12IeJa2FldNaa66Bb70tJS1NGADIWDlvAryccdawegsl0frr/nYZXq8xp6UrE8rGB5nW63OxR6uLy8nPp9ZYu9GMCumPLJOpFpO1ed/HKdn3ZIrSf8Wl5ezt3lVAXXh6lOaptiz6m4/Oo2sY9Fmudsq4vaUbDvPF9Y2D3J20pMW34WEckqkvr5WSw5X0InX8voSGJ3Pkl+9bxi49WJct+RmasTYku9OFjYPchqUefdIagpWkMmzmLnBMn7x+uafDXd5zJdPzGfu2ly1NaJxnwG+krWxcXFxPJs71TZC54mDRZ2D7KGBuZtpfR6vUHCp5B48tBzfOsSy1rUy9EjNeRqRNt9LkNM8nSJqOGMthDEmOGkMUcrtrkIPTlakfMndYaF3YO0Al2ksIRcS21Pq9US7XY7KI9I3oQO533TBeTt2405QjBlgzQtQNInUX03y/YV0JijFX10p9e12+0GL3bjTsAMC7snaV/w0ERGaV/WEBFTBUJav3n7PUPa5bI6TROGLkHwOSZWm2JYyzbB82mDbURjukYR+XGS8vqouzoBEO12O5Vxwj77UVjYc0T6hH2XRfu8rK4fjK+IyesUtV2abRhus+pUP3qr1TLOG4Ra43lHqsiQQB9r2YVJ8NSVw652ut4f9X7nPZL0eTZ6Byafta9rkH32bljYU+BrfYa+fEnH+1opviF6vj52n84k5EfoCpfTRxPNZjPxOJ/VkLEtPP36i4uLXtZyErrgNZvNoD1kTffAJLS+9yLNCFKfdJUdk74qttfriXa7PTTK0d1LrnqxxW6HhT2QkBcq9OVLOt63owito+uHlGQFuiJybCMIV1ik7v+3iVoasYrpk80ilkl1Un3sMVxIPp2gT4fge09tz1B2euq7ksZdabs2+9xfgoU9kFArPPRly2p1yvM7nU7i4iQfXO11+epdw3H5N9uWbVLYlpaWnEmt1DJDLXhfQkYraZ61bZI4pmAlvTe2DjrJZ+9ysenzArp7Sb4r3W43ShZHtuCHYWEPJMlKLeL6ScJvWlCS5XpJFrvJV5/UAXa73YE15xoJmIbvJgvP1JFk3YHHV3ht54a6qOTxaQUq6d0wRdro9TCJrmndgKn+NhebvLZ8Huq7orqwsvjK2ec+DAt7CuSLGiMGPMQiTDpf/4HGnCi0TXRKy1q/F72ee+LYdySg/821s45azxg78IQIr3ptX3E2HZdWoJJGSK731dVBJ0Uo+brY5PF6PaRFH7rhts+9nGR8hT1KdkciegjAewD8VAjx+hhlloHciHhnZydTljpXVkefjI/6MWfPnh1sQtzv99FoNKLsfmPa5Ue/9r333otnnnkGR44cGWQ33Hu/MPivSprdjkLquba2lq6xCXW0bbSt3otjx455ZTC0ZdhM03a1Xtvb2zh58iSEEGg2myAi505Ksh5ra2v4whe+gJ2dnaFry4249fro9XcdK49XOXr06OC7kydPYmdnBydPnsQPfvADXHPNNd6ZMWUmzbNnz2Jra2voPN4QOwEf9U/6ADgE4FYAz/kcX1WLXYg4FkJaq9V1jMn/nAdJkQ8+eWPSjFaSRgKhxyVh8ncnWdk+ydBc+CRvM9VTXlPNvKnmUpf/Do2ACnmn0swRnT59emR/Vl83oo+rcBKteBTtigFwsA7CLkS+k1o+L2WeL25S21ThVEPypHjE3nUnpG6hx2W9vmsSMfTaLpdKUn3UCUvpz5ax4Y1GQzSbzdTpm0PeNVvbbcZKr9cTrVZrSNhdBoFPmUl/qzuVE3YAJwBsAticm5tL1ai8ftB5ENPHHrNOSQtdOp3O4EcjxUP30frGJNcB1wgj5B7YRkIh/nx9tNLpdAZWfFLHEyrKpvuQxlhx5fh3wRa7mcoJu/pJY7FP8sOMhcuyUqNu1GFz1o2TdbJ2eHnj2+mGrmNQJxdNcfw+E5b689MXTx06dGhE6F31de2dq5PUAfg81263G5QcrOrvShnUTtgnefiVhSQrUIjRsDjdOk/rgvBxbajHJo0msrhCfOvrI9hJceC2MmUYrWlxko8/3+TOUcMW9U5Z/kZMv52Qnbhkx5S2c1eFXa8/G2th1E7Yq2Kxl2UpyB9XSHx9ko9Y/8GZLPa08eKma6fxm+rCGDsdsYqvYEuXiLxPrvvjalcaf77+nRq2aOqUpZWsi/Li4uLQ8YuLi8b6q66UNBtqq+3SRyqxYt0niUKFHcCfAfg7AFcA/BjAB13Hp508LXv4laZziVHnXs8/bauKS0hN1t/i4qJRJBAQLy7bu7y8PORXTVqk4zOaUCNBfIQgdCgfIti+GR9D3pnQkYlqSavRMrqImkRZt9i73e5QRy//q05+qqMAX1xzC2yxh1O4xR7yqXpUjI1Qd1CsUcbp06dHQtt8fJyu6yct0tFdM+rmFz5hb1NTUyPC4XO+qw2tVmsQqeMjkjYL2CUoeQu2D2rZNivZdIwUZNU/74pEUcMvZXl6B6G+d1NTU6l84a5ooKR7YxrZTKJvXcLCngO+fmBJrHmBJIvdVi+X+8b2g1P/q8eLJ4mZbl37WHohlqkM72u1WqLT6TjPU2Oo5SSwKR48aSQTU7B90eO/TTHqPqMxVaSTRnmueRbXPfe5X2nvk+kdnXQL31fYo6w8nRRsKwptq0nlCsft7W00Gg3Mzs4GXU9dXbe+vj5Ydamu7AOQuGqy2WwCAJ599tmhFXyyLbOzszh16tRI/W+66aahtq6urg5dZ21tbejv6orOVqsFIQR2d3etKy1N9022R72/6orgfr8PAPjMZz6Dfr9vXb07Ozs7OLbf7+PixYuDugsh0Gg0QESJqy5dqxrlitiNjQ2srq4Ojs+6KlKtO7BnfF26dAlra2uD8lyreNUVp5///Oexs7MjXaZWZHm//OUvh75/73vfi3a7bX1HTO+e3mbTCmcf9LLPnz/vtfKXAVvsMXBZT66kWCZ83AWmc1wTleonyfqTk4a6v1W32G0LlVx+Yv3fPpEgpjbq/mPTZswui913gZAPMaxKdWTlmuvQXSFJlnAa1+Hy8rKYmpoamTT3mdx2tTeN1c4W+yhgV0xxhPqyfcpxuQts55qiLEwCobtG9EnDqampoQk5fQHM8vKyaLfbwfXT75He6S0tLTn923qn51r04vKxxxKEXm8vukS9D+oGE3qnY5tD0HccUs/XI5VCkp6FuJX080JcLj5+8rSCzD72YVjYC8b2wiW91Op5egSBnCxMa51IEdbFQY1Pl6iThnoeElVodf+tayRia5tqnctIiU6n4x2RIstut9vG+qWNeQ853nQfpqenRyaN5f2xWZv6xLgq7gcOHBAHDx7MJOy+IbK+cx1pRDXWXBPDwl4pfEXftIAjRqikLjYy6kV3oehuFpN4652PyQ3i0zZ9ZaVq+eruIJN1aNqLNMtQPdSqNN0HPRpFHX3YYrZNFrvpE5r0LKQ9WSzq2HVh3PgKO0+eFoBt8kifHNra2hqZuJMTcbIcH/RJV3UiDgB2d3cBDKd5NaVqlROrespUddJuZWXFOGGY1DZgOBXszTffjMceewzAnrHxspe9zDixKr9rNBqDdhAR7rrrLmxtbaWeXPOZBFTRJy9XVlaG2tTcT6srU+UeOXIEFy5cGEwsv/DCC9jY2Bjcq7W1NTz99NP45je/uWdxabzjHe8Y3GvTc9brqqf7XVlZGTk/bdtDCZmQZiLho/6xP1Wx2Mv216nD+VarZdxKLsTSkUNv04YHpi3MWq2WaLfbwZsg2Pz56qYKPnXXXTX6hKhu4Sa5qqT167Na1taGEJ+8a0RhOl8+H33SWZ87MLnB0rj4QlxmbFGPB2BXjJuqvMzqkm3dNaJPzPlMvJrix3VBmZ6eFocOHRoSD1XcQzs8Pe661WolCmOSsJrcKqZj9DJ8c7rb4v5NEUD6xHGWbRP1CKSlpaWBADebTXHo0CGxtLQ0iExpNpvGjtfHbx3yDqUxcoo2jMo2xKoAC3sCsSd0Yk4s+Vpb6jVtOUPU5E7q8XqukHa7PTjGd7JX/U7146sdSsiEsknsXRaxz730vedJx8iRkGxjlv1WVZ+6vrpTlp00v+JrmORlwBRtGFXFECsbFvYEYr4oWcoyneszQaknZzJtaCCtQlOd9FwhrVZLLC0tOcMYXe3sdDqDpGFJMce+y/ZDiSl2+jF6x0lE1oljeb5NmNWy5D3Tn5mpbL1M3x2Z0hodSZ1o0lqCGNdRr5fGEKublc/C7kGsh57mpUvyxSZZzS4h161AW52WlpasURimTYiTFqmo55gsXls0S8wfna+7JOnZ6+VIN49+j1z+7ST/t4xAUiNpbGWr91gu7y8zkkU3DIDhBVQhnU5IZxyyOXYdrXwW9gIJfYF8LUYfa0n/yLA/W/5skx/ZFkftu2GDbr2p7iQ1fLLZ3Fs1qsbLh8Rl+xAiFC7/v81VtLy8LF73utc52+Dr/z59+vTIOoPrr7/eOGLS491DM12Gtj/JN297B5eWloyZI22EGEWhq7jrGD/Pwl4wvta/z48mCdXtoQ7ldZHxGQnIXXRMVqPNHaO3s9PpGH/IqnipE4a+P3pf1Dr5iqpL/F1lmCx33d9uKt/2fugiqFricmMO+ZxsoyvbBKvrftnar1rGPqMSvS7tdts717vPs/B9LlnLLoqsXgIW9gri+6NxoYvonXfe6W2hqjlI9PwpMzMzQ5aoyz+vl5uUs1v9gcloD9uxruskTcROT0+LpaWlxJ1+kgTCJcy+KXFNIyNTnfTRjnT/LC0tDYWt2ibHfZ+Tq/3qHra+C9BkG3V3XrfbDbLY9XuVdFyoUFfJxx6joylU2AHcDuB7AL4P4PeSjp9UYQ/50Zjo9Xojw19Zhk8EhS4MS0tLQ/Vpt9tDIZG6n9p0HX1InpSzO22SLJcwqkIrQx1dPvZQV5jeecj0wb6dsx7iqI+q9LKly0od1akpGPQOMnTkZ3OTJU1628qSk+6qgPv62EOpklCHEsM1VJiwA2gC+AGAfwxgGsC3ANzoOmdShT1rj637WX1dGbr46Ra7a6edpLqroxDTIit5TFIYo2/9bZE6+kij3W4njjR862CycKVV7ZuHJSmfvj4a0Hcb0t050t9vG534dPb6NdVwW9kJxwi5ZPYYK4sdwAKAR5V/3wfgPtc5kyrsQmQLPZMWtvStdzodr/JsLiA5yeqTqdHkMrBF9ZiunRSzntT2pCgTmW5W7bxiRdzo11ezS9o6Mx2fEE/TdZIifEz3MURAkjrsLC4tZpSx8bED+OcA/lj59wcAPOA6Z5KFPRQpWqofW06Whf6AVStMHXL7bBKt+011a9KGyQ2Rl690eXl5xPfs8n2HoJ6nu58ajUbi1n++HZyvb96FKyzV1BnbXGw2Sz5E/Jm4FCnsv2kQ9n9rOO4EgE0Am3NzcwXcgjgU7dMz/bBNC1jkMDptyJtrAs12vJp90fe6uhtiamoqOJe7Lz4ujxhC1OuNriOQsdVJIwv1HtvSSUjSRIHY5jFsvnzXHINaRqfTMYYaxvx9jLP/vCjYFROBPCySELeFKQJDCpbck1R+Z1ri7nMtX9+4Wjd1wk1Osrr8sXrkhBoVE9vSkyMck/sipuug2+2OWO2mxGU2V4n+bJOiiUJdKvrz0EdOvht2y47ClDIiJmz9+1GksLcA/C2AVyuTp//Edc64CHtsH2LSy2uyolUxVSfrXJEWPtcSwm0xms7XhV7dYFrWo9lsjuSnMeUbz7oMPQ2xxSMpjYLLcvbJ4CjrbOsY1M4rafSWxmJXOya1vjLJWyx6vfC1HZNq3Rcm7HvXwrsB/A32omPuTzp+XIQ9SdxCSSPGoT5b32slHaO7XpaWlpx+WP0j5wH0azSbo5t8FEmsMDx1xKMucdcFUTcMbPfMJuym6+pJxDqdjleahiS/ut421Q2j1rXT6WS+d/pIMCR8dJKt+0KFPfQzLsIuRJzJLLUsdUhr+gGGdBxJP86kH3u32x0SGT3Xh8nKVof5+t6j+ke6eGwugqKtLtMIJW0dfEZzJveIGqKYZj9TU8ir6lqJkaZBb5tqTWd1w7jcjb6juLySyI0DLOw5kMY1o3cMSX5V1/m272zi5PoByB+YLsYy/vv06dNGsdbjq6VYHTp0yLgy0rXFXehCmCydQK83ukJWdXWFdtS+9ZduE3VFrFwroM89mARZd7vYXFu6tZvlfpk6JDXPfZYtG01hs6HvQZ5J5KoOC3sOpBEj3TrxWaXpup7pR2erk6u+NpeA6i/WBUROtNpWRJrEZmpqKtVEZsyRkryeWsdms5k5Oscmnqq7Ryau0t0u8jx9YxC93SYR6/X21jToHalc/akvXEojxGqHIt8FuetWkn/ehSndQEgnZHMfTorP3VfYa7nnqWsvyCyE7t2o7yUJADMzM9je3kaj0cADDzzgLMO0FyWAoe/Onz+P7e1t9Pt9bG9vD+1XqdZ3dnZ2aO9UuWfnpUuX5DwJgL2OXu5RKvfiBIBbbrkFW1tbmJ2dxalTpwZ7fcr6qPuqEtGgzCtXrqDb7eLcuXN4/PHHB3XT9wyV+6ACe89P3e/02LFjmfbk3NjYwAsvvICpqSlcuXIFRAQiwubmJvr9PhqNxkgdfND3st3Y2MCZM2fwpS99CQDw2GOPodFojNwbea2FhQU88cQTQ/vA6u2+cuXK4FzZ9vvuuw9f/OIX8da3vhVPPvnk4O9XXXUVTp06NfRMt7e3cc8996Df7w/2jvW9d3L/1kajgZ2dHQghsLu7i93d3cF7Evostra2Bvek0Whga2vLuiewCf29OXr06Mj7EtLGulI7Yc/7IWd9CY8ePerdMdjET/1O3QS63+/jZz/7GVZXV4fKf+GFF7CysoLd3d2he/L4449jZWUFX/va1wZCoAuPXseNjQ0cO3YMAHD06FEAewIgOxciQqvVAhHhypUrAwtCFwFXJ6lvxPz000+j2WwO2h4iwOr70Gw2ceLECVy8eBFf/vKXIYRAo9HA/Pw8br31Vu8y9fJlx3nq1Cn88pe/HPq7KuqNRgN33HEHrrvuusF36j1eXV0d6sAuXrw4tGm33vYbb7xxSNgvXbqEy5cvDz1LeX6/3w8SYvUZ9Pv9QZlCCLRaLQghhurja0wdPnwYMzMzxg7dB9N7o9+32JtxjyU+Zn3sT56umKosc5ZDwyz+SLUcm49d37tUJomS4ZFyqA/NDaCWpYbBpUmeJYft6iRYp9MZCof03RxBvY5PDpskTCGk+oKptG4L9X6oLirb59ChQ063mboRuVwrIJ+nGkKqnqP7vtXwWOlGSePG0tumv0Pqsyh6vsRV17pHyWBSfexVeMhp65Dmhb/++uudYqJ+bOldbdfVv3d1mvpy/kajIZaWlsShQ4cGESCh9yJr3npZji0KQyYLUzsk31QJ+v1Qz5Wra2UYIhENOk3T/TN1rrZjTe1LmkhPK6SqceJaoORrTMWYAHdFgbGPvcbCLkT5Dzlt9ExoZ2Cy2HXrSgpYqMXrM3E0wvsWAAAX3klEQVSrlqULu+ljGi24nlPWTto0anJNPusTwz4RS6ZQTv166iSkqT2m98U0sep731z3Im2UjExA57OYzTTqifEsyzbYqsBEC3vZpHkJQzsDfah8/fXXD7I9qlERSe6V0Pq4rHvbPqym0YJaf1enk8XadLk9TFZu2lzxelmuqCVf0evtR8RIaz+N28PUrtD3wfacfO6hulradwRioyou1rJhYS+ZEEHSxTjUDWCyhH02dHbV2WWN2sqW2+wdOnRoSORbrdaIIOi5UpJ2AQoV+LRCkHW0Z1rc4+umUO+rrf4hbg+b/9/UydrabBtNhIbR6nMZbLGng4V9TPC1XF3nmSzhkFzhtvJ0K0yfHHRZkp1ORzSbzSGLU/cl63uH2oQqrZsqhhsnqwCp925qakp0Oh3jytskS1/Ww3dDZ30rPf3Z2URaH1GY7mOSUXH11eaMpEkZRPN6JnXCV9hrF+44bqhhZQAwNzfnFaplCxdcX18fhB72+33cc889uOmmm7xi7tV4+LW1NVy8eBE/+clP8MpXvnIkVh3Yi1GXoWV6zP0zzzwDYM9w2N3dHYnBF0LghhtuwPe+971Buc1mcySe/cyZM9jc3BzEZl+6dAlra2uJ7dHvD4CRMFAbWUJmbc/l5MmT2NnZwZkzZwDshT4SEQAMYtbV+7+1tTVSzsbGBk6dOoV+v49ms4mPfOQjQ2sTJB/96Efx4IMP7llu++hhlzLUUA8tvfvuuwHsrbdQw2L19tjWIADAsWPHcPHiRTzyyCODOPxWq4WjR49mCkMMCTWeeHzUP/aHLfaXiDXEVK1r2+bSLotHXxGoW3gyxwo8LXZ1taJusavJxaBZdGp7bD57U4rimPc4hj9Xvdc29wQU61mfBDeNtGxROLpFb7uOPE/Nx6L68vXjXO02vUu2UZt8ZllDfxl2xYwVvkNM18Sl/oMyJbtyiZu+5F7/3HjjjYOIG1dctepakCFyauqEXm84jFGKuslX6xLCELHNMjGtu6XSuIGSkqVJ37N6P5JcUrYInna77RR100Sy7hLTn1na+6zOL0jXYNV95FV397Cw71P1B+WL74SV6jv1jUFXy7eJjwx3M51vusd6+KO6ubTua7dtnq1b7KowpF1IlCbaJYbFL8u78847R+6rTw4hdURmi+BRO1O1EyQir5w9slPxnZvRn7trnkAd8VU1qmUcJmhZ2MV4PChffCasXKF9PqF8Jmv6wIEDQ6tIbZO1+vem/Udtsd42er29+GmZ3MrUFlPHYHMTpO3gY1j8KnpOeLVTNYmqtKr1mHZX5w1AvOIVr0gUU7WusSbvY4WSFs04hFSysIt8HlRZIwCblatbci6fp6+Y2qwued1OpzMQJle8uz7ED13ZacNkZfq6nNIQUqbreSRdw3a83kna8q3b5jCSImjSvNO6v1/329tchlUePY+DIViIsGNvI+vvAOgDmPc9b1wt9rIffK83nE8kaYuzGLHcehl6DvFOp+McLajx+eqq2CwbNsjnYMqBo84VZL2GyfIM6RSzvB/q9VzCbqqnvjGGFN00bUpqq77zUVIa6SoLuxDVr2NRwv46ADcAWK+isAsR90FVYaim1iFpU+IYQqOXoU/OSdFw3WP5d32LtSzb0+mdnGyfKd931jb73rdY74dp1BSyCtX0vWsk5nKlqHUyjQj1fD62BVllG0V1wVfYM8WxCyGeBzCIx60iMWNfXTnEXbhSmobmjlfr0Gw2QUTY2dmxplANyR9vQo9jfvbZZ/GNb3xj8PcjR44k3mP599XV1ZFc3GmRZeppkNfX14eu8dWvfhXnz5/HkSNHcOLECa+yTXnwQ59NmpS0tuvL3PimNQumevqktj1//rwx178evw8Aa2treOihh0bSPi8sLGBlZQUXLlzA9vY2iAg333wzLly4MHIP0t5TJiU+6p/0gYfFDuAEgE0Am3Nzc3l3bLkROgJImthMY8WodTANr/O2jDqdjnjNa14TvKmx74Sb6bw0ETB6VI1pw+mQevq2MesI0TafkqWePha7PuKQu2O50j4LMbqXbB5JwEztqbLLJC8QyxUD4OsAnjN83icChF39TFIcu2t47pqAyuN6Mcj6Aw3tiNJcT17D5jbyKS+2cKQxCGyRSCp6hI1PHWyRSfq98c3n45pAj+XT19sxqW6daMLuVciEC7vrhfWx2PUJqFjWnupfddUv5McWuzPS5wz0iI8sHZXub5fRPEXHU+fls0/b6Ul/vSm2XR8NyjBLWxy8rR55im8V5rrKgoW9IHxe4CRhTSs2PhNdsV1BsTsjKTRSfE25x7NEWUiLVkbvxOxEfUkTB+8T++2zME3HlDff9d6aJm1ddZbHmNw6NvdhKGyx5yzsAH4DwI8BbAP4nwAe9TmvTsIew3rIIrBJ5/i6gkI7FFtInc+5+o9azURo8+Hq7obQe5aH28vVJv3vaXzhaj5zH0tZzyRpcs+YhN0WFprl3dbnCfQt/3xHlK7y2cees8Ue+qmTsMeyHrK4RFw/utgWu35uiAVsu5465Pe12GOvBjUd7/M8YvvsXYuvkspVY/iBvRQMptGcWr7tOLVtrhWxPpPe+naEaicuJ2hD98WdVFjYCyTG0DKNxRJiBcbysevn+riRTD9wfZLNNuR3TcylnVSNJdZCJI82XHUw1UdeOykSxVRm0t6kEpkBVPrOXWKqR7yonXHaUYhqset5chqNRtAeApNmtbOwF0zal0x/6X1zdcjdikxZFovEJcq93ug2faYoD5f1ndXHnoaQ0ZBrfsB0vE3gXPctyWJX70+n0/HaiCPNCEK9F2nnDfQOrdcbTVzmk1lyUv3svsLOG21EIMvGDOrCjd3dXXS7XZw7d85ZxoMPPogPf/jDg3+/613vKnWxh1ygRkR49tlnsb6+jtnZWZw6dWqwQQawZ0TMz8/j1ltvHdp0IWlhz7FjxwBgZKOGGIvPTAvEfBcara+vY2dnZ9D2u+66y1kf9VnLjS+EECMLdmyLr5LKvHz5Mq655ho8+eSTI+ep7ZTn+Sxam52dRaPRgBBi6F743CP93urPTvLAAw/g7rvvHtwTuTGL772MseApdKFg5fFR/9ifulnsMSaZQobei4uLI/HZSdfIa8iqT0rK1LqmfTblMbZ49bwXteiYylfdGj4RN+pzSHIh+FrsWdvgOmZ6elpMTU0NzWfY3o8QH7teRuizs7l8srTbl3Gy/sEWe3FkWUoul3/ry7ZdZRw5cgSPPfbY0L9tqKOJRqOBW265BR/84Aety+t1y25tbQ3AqLUsUdveaDQG1mij0RikPGi1WnjjG9+Izc1N9Pt9o4Vlsr51q2xtbS2qVaWXf+bMGTzyyCPY3d1Fq9XC8ePHnedvbW0FpUgwbdeXtT22rets7ez3+4MR1Pb2Ns6cOYNHH3108O6ePXsWW1tbQ9vm9ft9ENFI++Qz29jYwOHDh3HlyhVMTU0N0h+EWNQnTpzATTfd5H0/fNrtSy3THfiof+xP3Sx2IeItJfctw3fFoR5lAYd1qVpo0uqWx7vimFUrVw+/842nd1mMMS1cV/mm++TadKKqlp7Letbb2G63h0Zcalpl3xzqpuyTVb03JsapruDJU0aIl15aXbBM7hs9XE792LZrM13PJwLH1hnYXDS2iJqsqOXb2q5P5iVFtsSoT2zXjJyQXVpaGkTEyH1I5fGm7fZ86mNLKzxOUSvjUlcWdmZAr9cb2SzZlMCr14u3gXRSfaSY+GyZFhIlkrVOcrSih+BlCbMMrUOWsn3CQ2X0lZ4zJu0uR72efR0CExdfYWcf+wSwsLCA22+/HRcuXIAQAkSEa665BsBoNMBnP/tZ3HPPPQM/87vf/W5cd911Vh97GnSfr/THm+YW1DmCVquFD33oQ1HrIjGlJ5b3YWZmxpl+Vn5ftK83JKJHLVty6tSpkUiuED+3ZGFhAU888UT0qJLaRaoUiY/6x/6wxT5MyDAwzZDRZvG6hu15DkulhQfF1WGL3y8z4ZNPpE6svTxDRyXqCENdZCTLMSX3UuuZl2srFnIRVZE5fcYBsCsmP2IKX8jwO81Q3TQEl+fFynOT5l74rtisysSWza8eO1eQ7yI1U/oA3/mKpMnssun1el4raCcRFvaciP2DCFmSnkZE0q7q9CHL+aEdWpkTW0kRPXn5xZPqpIufbVs6GyG53IvEJ+eNTtnvSFGwsOdETNdAr5d+SXoai913YZAvWbJD+i4C8iXPH3ZSO/OKZElCd1eEuIWqbrGHJAarcltiw8KeEzFfIlUwTJtM2K5flRSneXQ0RdWjSuXLa+ix574Jy2zn2eYI8gwfjUXIO1vmPEzRsLDnSCyhLMvSSCsaWcqS5PEjLOKHXeRQP8Z7YSpD9+XnGT5aJGyxj34yhTsS0acBvBfAZQA/AHCXEOJnWcocB2Ikn5LluJZF5xHu5ZuwTD9OXWqeJRFXlvQLRZapE+uZ+xBjifva2togAZsalqmGPH7oQx/C3Nzc2IcTxkwvUBt81N/2AbAIoLX//58C8Cmf88bdYi+CvKwQPWmXbSchW3IvV6KpPEM2yyizLGJMaqvhpHJxWdIkcF3uX51B0a4Y7G2T97DPsSzsyeTlXlAnpqRo2yZVfVeH+gjRJApHljYn+cpd6HM37Xbb6bsvYw6BSUcZwv4fAbzf51gW9mTy/LH1emE7H6XZVLmotlSRXm9vkZBpU5E0ZYVGNUmLPSlVsiTvOYpJe/554ivsiT52Ivo6gOsMf7pfCPHl/WPuB7AD4GFHOScAnACAubm5pMtOPHn6DRcWFrCysoILFy44/dKqX1ldag4Aq6urg3ol+bhDfMYh8wplLjm3XVvOTagbjGRJBeu6d675EvGSwWVNlSzJe46ilmlxq46P+rs+AI4B2ADwD3zPYYu9GqQNnQxNQ+BrsaluIle63JAy86DX89ujFUjejNrnWrZ22ixtU6rmpCRu6ugsj/kPttjjgCJcMQBuB/BfAVwbch4L+/iSdVGS60etrzh07X2ZZg1AFtT629LUyuOkiE1PT3vvYeu6rin3i34tvZNV51F8sy7m7f5jH3t2ihL27wP4EYBn9j9/5HMeC3v1sf0Q01jsIWWraYNdOUJskR95oLd5aWnJKuyu9mW9bsiktO88isokLfQZVwoR9rQfFvaXqKIlEzphF2LpJR0bsvdlSJ6dLOiCJydG884/nlVoQy1wdplUHxb2MaCqP6RQQQk53udY386uqPtnuk4RHXKM9oXWs4qGBvMSLOxjQFWHvnlaerHFuCghKkvwWGgZFV9hp71ji2V+fl5sbm4Wft2q4bu8vwxCQwnHJUyRYcYZInpKCDGfeBwLezG44p5Z5BiG8cFX2HnP0wJwWeZFJpcqA+64qgU/j8mAhb0AJnXlXZVdTbEYJ6GchOfB7NEouwKTgFyy3Ww2c0srC+z9cFdXV7GxseH8zvV9TEwdWpUJvSdSKD/2sY/htttu87rvZTJuz4PJgM8Ma+zPJEbF5B3dYAvJc61MLCNMsKqkqastqqmq7S6iXhzFky8oYqMNxp+8fek2a8zkAirKNTROGyCYNqZIqq8teVZVXW95Pw929VQHFvYCydMfaxMZ03dF7DgkGYfJ4Y2NDTz00EN7CzsAtFotr3tiE8oi728oeT6PqnZokwgLe0Hkbc3YRObs2bM4f/48jhw5AuCldLvjYkkXwfr6+mC7OCLCXXfd5X1PTEKpPwtgOM1xXalyhzZx+PhrYn8m0cdexipTPdNgjI0fqkBsP27eWQ2r6G/PC/ax5wvYx14tyrBm1KFxv98HgCAfchXJY+STp+950twT4+B6mwRY2AuijIlEtTNpNpsgIuzs7Iz1MDkvocxLkNg9wZQBC3uBFG3NmHy94+5XHzehHKfIIKY+cK4YZuwYp9WeDBOTQnLFENEnALwPQB/ATwH8SyHET7KUyTBJsB+XYdxkTSnwaSHEG4QQNwP4CoDfj1AnhmEqThVTJjAvkcliF0L8XPnnr2BvD0gmZ9gVwZQJrzCtPpknT4nokwCOAvjfAP6p47gTAE4AwNzcXNbLTiz8o2LKZtJCOMeRRFcMEX2diJ4zfN4HAEKI+4UQrwLwMIB7bOUIIR4UQswLIeavvfbaeC2YMDhDH1M2RWUrZdKTaLELId7uWdafAngEwMcz1YhxMm7hfkz94BDO6pM1Kua1Qoj/tv/POwB8N3uVGBf8o2KqAEcmVZusPvY/JKIbsBfu+EMAy9mrxCTBPyqGYVxkjYo5EqsiDMMUD0dY1RNOKcAwEwpHWNUX3vOUYSYUjrCqLyzsDDOhcNhifWFXDMNMKBxhVV9Y2BlmguEIq3rCrhgmM5wQimGqBVvsTCY4soJhqgdb7EwmOLKCYaoHCzuTCY6sYJjqwa4YJhMcWcEw1YOFnckMR1YwTLVgVwzDMEzNYGFnmAmFw1TrC7tiGGYC4TDVesMWO8NMIBymWm+iCDsR/S4RCSJ6eYzyGIbJFw5TrTeZXTFE9CoA7wDwQvbqMAxTBBymWm9i+Ng/A6AD4MsRymImBN65Jxsx7h+HqdaXrJtZ3wHgfwghvkVEkarE1B2euMsG3z8miUQfOxF9nYieM3zeB+B+AL/vcyEiOkFEm0S0+eKLL2atNzPG8MRdNvj+MUkkWuxCiLebvieimwC8GoC01g8AeJqI2kKIi4ZyHgTwIADMz8+LLJVmxhs5cSctTp64C4PvH5MECRFHY4novwOYF0L8fdKx8/PzYnNzM8p1mfGEfezZ4Ps3mRDRU0KI+cTjWNgZhmHGA19hj7byVAhxMFZZDMMwTHp45SnDMEzNYGFnGIapGSzsDMMwNYOFnWEYpmawsDMMw9SMaOGOQRclehHADwu/sJuXA0gM1awA41DPcagjMB71HIc6AlzPmLjq+I+EENcmFVCKsFcRItr0iQ8tm3Go5zjUERiPeo5DHQGuZ0xi1JFdMQzDMDWDhZ1hGKZmsLC/xINlV8CTcajnONQRGI96jkMdAa5nTDLXkX3sDMMwNYMtdoZhmJrBwq5ARJ8gom8T0TNE9BgRvbLsOpkgok8T0Xf36/pFIrqm7DrpENFvEtF3iKhPRJWKQiCi24noe0T0fSL6vbLrY4KIHiKinxLRc2XXxQURvYqIniCi5/ef971l10mHiK4iom8Q0bf26/gHZdfJBhE1ieiviegrWcphYR/m00KINwghbgbwFXjuDlUCXwPweiHEGwD8DYD7Sq6PiecA/DMAT5ZdERUiagL4LIB3AbgRwL8gohvLrZWRPwFwe9mV8GAHwO8IIV4H4M0ATlbwfm4DeJsQ4o0AbgZwOxG9ueQ62bgXwPNZC2FhVxBC/Fz5568AqOQEhBDiMSHEzv4//xJ7u1dVCiHE80KI75VdDwNtAN8XQvytEOIygH8P4H0l12kEIcSTAP5X2fVIQgjxd0KIp/f///9gT5SuL7dWw4g9frH/z6n9T+V+20R0AMCvA/jjrGWxsGsQ0SeJ6EcA7kR1LXaV4wC+WnYlxojrAfxI+fePUTEhGleI6CCAWwD8Vbk1GWXfxfEMgJ8C+JoQonJ1BHAWQAdAP2tBEyfsCZtzQwhxvxDiVQAeBnBPVeu5f8z92BsKP1zVOlYQMnxXOett3CCiXwVwHsApbeRbCYQQu/su1gMA2kT0+rLrpEJE7wHwUyHEUzHKi7aD0rhg25zbwJ8CeATAx3OsjpWkehLRMQDvAXCbKClmNeBeVokfA3iV8u8DAH5SUl1qARFNYU/UHxZC/EXZ9XEhhPgZEa1jb/6iShPTbwFwBxG9G8BVAF5GRP9OCPH+NIVNnMXugoheq/zzDgDfLasuLojodgAfBXCHEOL/lV2fMeObAF5LRK8momkAvwXgP5Rcp7GFiAjA5wE8L4T4N2XXxwQRXSsjx4joagBvR8V+20KI+4QQB/a3GP0tAP85ragDLOw6f7jvSvg2gEXszVBXkQcA/EMAX9sPzfyjsiukQ0S/QUQ/BrAA4BEierTsOgHA/qTzPQAexd5E358LIb5Tbq1GIaI/A7AB4AYi+jERfbDsOll4C4APAHjb/rv4zL7VWSV+DcAT+7/rb2LPx54pnLDq8MpThmGYmsEWO8MwTM1gYWcYhqkZLOwMwzA1g4WdYRimZrCwMwzD1AwWdoZhmJrBws4wDFMzWNgZhmFqxv8Hforp/K6E1OEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"proj = my_pca(400, 400, 1000, 0.15)"
]
},
{
"cell_type": "code",
"execution_count": 192,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(800, 2)"
]
},
"execution_count": 192,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"proj.shape"
]
},
{
"cell_type": "code",
"execution_count": 193,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.27719121, 0.88676956],\n",
" [ 2.03233192, 1.59504041],\n",
" [ 2.03770526, -0.50666795],\n",
" ...,\n",
" [-1.12542688, 0.21465377],\n",
" [-1.4274677 , 0.05756537],\n",
" [-1.50212695, -0.16900765]])"
]
},
"execution_count": 193,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"proj"
]
},
{
"cell_type": "code",
"execution_count": 194,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.cluster import KMeans"
]
},
{
"cell_type": "code",
"execution_count": 217,
"metadata": {},
"outputs": [],
"source": [
"kmeans = KMeans(n_clusters=2, random_state=0).fit(proj)"
]
},
{
"cell_type": "code",
"execution_count": 199,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,\n",
" 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,\n",
" 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1,\n",
" 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1,\n",
" 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,\n",
" 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,\n",
" 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1, 0,\n",
" 0, 0, 0, 0, 0, 1, 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,\n",
" 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1, 0, 0, 0,\n",
" 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1,\n",
" 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1, 0, 0, 0,\n",
" 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,\n",
" 0, 0, 0, 0, 1, 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, 1,\n",
" 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,\n",
" 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
" 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,\n",
" 0, 0, 0, 0, 0, 0, 0, 0])"
]
},
"execution_count": 199,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kmeans.labels_"
]
},
{
"cell_type": "code",
"execution_count": 218,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztvX+wXFd15/vZ3VeSxwTwy8XGBksjJjFUXAhJtp7hwiM5xJ7ENj8TzVQxMJGwjC+hLIJDKIOG58lNUVPKmF9KjAHLSIpvPWYyUyUwiZGxkeODFWiw9TMGYwI49pWR9Gxung3E8ZVu935/7N599tm9z6/u0z9u9/64XNLtPmeffY6k715nrbXXElJKPB6PxzM6VAY9AY/H4/GUixd2j8fjGTG8sHs8Hs+I4YXd4/F4Rgwv7B6PxzNieGH3eDyeEcMLu8fj8YwYXtg9Ho9nxPDC7vF4PCPGxCAu+pKXvESuXr16EJf2eDyeJcuhQ4d+JqU8N+u4gQj76tWrOXjw4CAu7fF4PEsWIcQTeY7zrhiPx+MZMbywezwez4jhhd3j8XhGDC/sHo/HM2J4Yfd4PJ4Rwwu7x+PxjBilCbsQoiqEOCKEuKusMT0ez7BTA7Y3f/UMC2XmsX8Q+AHwohLH9HhGiBoQAgEwNdCZlEMNuBw4DSwH7mM07mvpU4rFLoS4EHgz8MUyxvN4Rg8tgjc1fx0FCzdEiXq9+Ws4yMl4DMpyxewAbgQaJY3n8YwYIaMnggHKUq82fw0GORmPQdfCLoR4C/CUlPJQxnHTQoiDQoiDTz/9dLeX9XiWGAGjJ4JTKPfLx/FumOFCSCm7G0CI7cAfAIvAWSgf+5ellP856ZwNGzZIXyvGM36Mmo/d02+EEIeklBuyjus6eCql3AZsa140AD6cJuoez/gyhRd0Tz/weewej8czYpRatldKGTIaUSGPx+NZsniL3ePxeEYML+wej8czYnhh93g8nhHDC7vH4/GMGF7YPR6PZ8Twwu7xeDwjhhd2j8fjGTG8sI8AteM1th/YTu34KFQM9Hg83VLqBiVP/6kdr3H57OWcrp9meXU59226j6mVftu6xzPOeIt9iRM+HnK6fpq6rHO6fprw8XDQU/IsOXwXpFHDW+xLnGB1wPLq8pbFHqwOBj0lT6n0uiKk74I0inhhX+JMrZzivk33ET4eEqwOvBtmpOiH6Ia0NwDxf4eWOl7YR4CplVNe0EeSkN6LboBaNPTiEZQ8vmcQeGH3eIaWgN6Lru6CFOIbgIwOXtg9nqGlX6LrG4CMGl7YPZ6hxouupzhlNLM+SwjxoBDimBDi+0KIPytjYp5i+E1KHo9HU4bFvgD8tpTyl0KIZcDfCyHullJ+p4SxPTnwm5Q8Ho9J1xa7VPyy+eOy5v+y23E9+bE3Kc0em/XWu8czxpTiYxdCVIFDwK8Dt0opv1vGuJ58mJuUqpUqe47uYbGx6K13j2dMKaWkgJSyLqVcB1wIXCaEeLV9jBBiWghxUAhx8Omnny7jsp4mepPSx9/0cbas28JiY9GXGPB4xphSa8VIKZ9B5WZd6fhup5Ryg5Ryw7nnnlvmZT0ocd/2xm1sWruJ5dXlVEXVlxjweMaUrl0xQohzgTNSymeEEP8GuAL4713PzNMRvsTAUqbXdWE840IZPvYLgDuafvYK8L+llHeVMK6nQ3yJgaWIL8blKY+uhV1K+Q/A+hLm4vEMMWnWdBmWdogvxuUpC7/z1OPJJM2aLsvSDvDFuDxl4RtteDyZhLRb03m+K4KuC/NxvBvG0y3eYvd4MglItqbTviuKrwvTO8YrMO2F3ePJJK3KYt4KjMMqLMM6rzIZv8C0F3bPQKgdrxVOyezknISRKC5madZ0lqU9rMIyrPMqm5BxC0x7YfeUTpYAd1K0rLxCZ4MQs5DhFJaQ4ZxX2QSMW2DaB0+HgFEquasF+Kb7b+Ly2cud92QXLctT9qCTcxJGopxgZxEClKBUSRaWGrC9+Wu/CMie1ygwfoFpb7EPmFEruesSYPt+zKJlecsedHJOwkj033rL8sMPyiXSz7Z4g/blj1dg2gv7gMkjhEuJPALcSdmD8kol9ErMsoQrTVhCynOJFBXQfgjeuPjyhwcv7AOmPEs0P+UFIdvJK8CdlD0or1RC2WLWrXAFlPMWMawCGjIevvzhwQv7gEkTwrwCXESo++H6Gb9aNSGRcD0PzFJMuMp6izDnMUwCGjBuwctB44V9CHAJYV4BLirUnbh+emnhjwYBKgBZRzUP2wNsori4u44v4loJGE4B7acv3wNe2IeWvAJcVKiLun5GLbjbG6aALcBtKGFfpHtruYay/HejFow8rpVhFtDxCl4OGi/sQ0peAS4q1EWDkObCsbC4wEw4w0ww48W9jU3AHZRjLWtf+fNE7YPTXCu2VT+sfzaDzowZH4SU/e87vWHDBnnw4MG+X3dYKNt33ktXibbYFxYXaNCgIiqsqK7wljvQLlRlCdd24CaUpQ4ggLNQ1jjGNaC4Vd8tnd7jsAZ2lxZCiENSyg1Zx5XRQWkl6m/X+UAD2Cml/Ituxx1Virg28gYhexms1Bb+TDjD/n/aT0M2RiIts3uShKqMZxIQ+congGtQbwQY16yiBP80+az6MuhGnEOGM7A7mpSx83QR+BMp5W8ArwOuF0JcXMK4I0l5Oyj7x9TKKWaCGVZUV/heqi1CereD1dwpeT/w+eZn5jXPEBd1QT4XUDc7XM3rF73ngPHY5ToclNFB6SRwsvn7XwghfgC8HHik27FHkSyf+LBmoOT1zeepE2N/P6z3nE5AegZKt24Zl/VvXlOgRF0Cy4is+l4WIzOvX1SchzmwO3qUGjwVQqxGtcn7bpnjjhJZeevDnIGS5fLJmr/re2Co7zkZLVSzju964U/WC8UO4AgqpfIMygL+S2A6xxgh3blDuhXnYQ7sjhalCbsQ4leAvcANUsqfO76fpvm3b9WqVWVddkmSJJBLvbxA1vyT3FB57jnNqh+sxa8zYe4gEvCQcv3J9kLxu0RuGAHM5xwnaJ6/0DxvsoO5lCHOPjum15Qi7EKIZShR/5KU8suuY6SUO4GdoLJiyrjuqDGI8gJpFBXMrPknfZ91z2lvAlnf9VbwQ9wCHlDuRiHzOgvA3xL51iG/QE+hLP6tzbFuANbQX3GtoZ7HGZQLKezz9ceDMrJiBLAL+IGU8tPdT2l8Ka/QVfd04hbKmn/S91n3nPYmkPRdf9xaAW4BL9ufHKD+qTaIfOuaogI93xynwWCyU2ab16X5a9HyC548lGGxvwH4A+BhIcTR5mf/RUq5r4Sxh56yrcIyUhfLmFOvSg+47i/rntPeBJK+649bK03Ay/Ynm9kvE3Se5hgwnGUHPGVSRlbM36P+to0Ew1ZQK+v6rgyTTudkjjdMpQfS3gSSvuufW6tMAU/yPYdEdWgkKgMG4huTgpzX0O6YvcDG5mfbHdfsFZtQgV+9sGxKP9zTEb6kgEE/CmqVwc5DO9l1eBdHTh2hIRuxuXY6J3OHaaVS4darby3kFpo9Nsvzi88jkT15FmlWfdJbQPdurX4G+ZIyaWrAHCr7BSIxnGr+WnR+NZTr5jTwTdRCUdau1TzPawqVm591nKcbvLAb9KKgVtmump2HdvK+u94X+8yca6eWavh42Cob0Gg02LpvK998zzfZ9sZtmefWjtfYc3QPsukaqFaqAw/8QrdurX5vgQ9xb/7Rc5gAriOeq97J24J5nUbzM2ldM6S3JQN82mOv8cJuUHZBrV64J/Y+srftM3OunVqqweqASqVCo6H+sdcb9dwFv8LHQxYbiwAIBFvWbVkiaZppFmaIO+OlV1Z8QLvv25wDwKoSrmleRwdkF5u/n8SXDBgNvLAbdNqyLem4bl01Lmt/48Ubufexe1vHvONV7+DGN9zY5nfupDvRrVffytZ9W6k36jRosP+x/RyYO9DaSJT0XOwFcdPazvym/c1Hz7IwA9qFtpcbjwLcgdiyA51mwHcS+EDzc4na+NSpOAc9mKunU7ywW5RZUKubAF6StT99qdphuPeRvWy8eGPrZ9f5RUVy+tJp1py3RhX8emw/DVTBr9ljs9xx7I7EN48y/Nn9D0SHpIuYK+Nle8Y5RXEtFNuI6rkEjjmUgR5nBmWtaz87dFYmwdwVO1/yXD2dMPbCPgz9P12kWfvTl04nCjp0J5K64NeBuQOt8yF7d2i3C2LR4Gv3f24B2Ram7QvOc04RQtL96qbYQ3luIL2gLKBcMRWioGxSQDYtuOvL8Q4bYy3sZaUGpp3TqeCZ1n61UmXu2Tlqx2u5xurWBWQvSEDMYi87MJoWfC07pTOik01Evdh4ZPq754g28Lh8+2UJaNgcR4v6FSjr3QzKJp1jzyvpc88gGWth7zY1sNcNoe/bdB+zx2bZc3QPtx++nTuO3ZF5rdrxGnPPzjFRmYAGHQuxvSCVtSPWJdRJwdek51xemmmR7AzTWs7OFMp/fV1IbDdwO+qfpJnaGDR/P0vUUalbAQ2Iv3nM5BjLPidI+Vy39YMoi8fXh+knYy3s3aQG9iN/XYvYYmMx17VMIRRCcOkFl3LtJdcO1Y5Yl1AnBV+TnrN9/OTZk2w/sL2HQdeyrGWX4Gmrt07k574OlQETEInibqKdphMku4Hy5pKX9bZif07zV102YA+q+qTOnc/z/Pwi0C1jLezdpAZWK1Ua9UbunO1OfcJFRMwUQiQ8dOIhHn7qYdact6Zngcgi95XkR3e5frYf2M7k2ZPOhdc8fvLsSW74+g0lvz3ZwlKGtayLX9mCN4/KTjGtXruueki8Td41CdfvdS65Pj60fjbH2o4q8KVZQJWSyuuu8T77MhhrYYfOLVHRrKIgclRTyOO6SRLIIiKmFwEtnr3aBZp2X+BOi7T96EKIWNxA/2+PuePKHcw/N+98LlMrp9h+YHvJb0+2sOwgv7WcRki74G1F+bn1dZIySgLahd+cb0h73nsv/N2uZ2PPOUBVbTxtnHeYSGqygs4h3mffPZVBT2CYqB2vsf3AdmrH09uGafeIRLLYWGzVFU87Pq0dnhazm+6/ictnL2+7/tTKKba9cRvzz82njqMXgfdd+r5cbezS7jfPs7Dva/bYbOJ9mH50ACkltx++3XmcOaYW9fDx0DkXvZiV17IvJC4se8lnLWcRoARPUyVyv5xGCeS2hLG1u+PjxC1YLbQ3NX/Vln+e9nOdtMgLiZcP3tq89puA9zfH0q6ly4hKSOn6Nvb8XQQF7sGTxNhb7JoiAdGivvms4/P67PNcV1uym9ZuSnWR5K1xXq1U2bJuC5vWbop9r98gzPlAclqkOXchBA3ZcDbGdrme0v5cyqkJE3vKxK3jjcAB3NZyGrY7Rwue9rGvJ+53DjLGc7lOQuKL0DztfnCXvzqPu8N1XkC8NZ8uS1AHbkO92WxBPaMdzeN1l6estn3mfea5B08aXtibFAmIFhWTrOOL1Jz5wGs/wNGTR9l48cau0izz1jiv1+vcdui2VkYOkOgqAZUWuVBfoCIqTJ4dNYAws3xO/csp9v1oH/VGPXa/+h7NMZPedsxn2Zk7LUksXMKyJuHYtLGThNMMihYd1yagPSPFXADMeVSJRDckWhCep70mumv+EN+ENIlamHTsQccfbkN1k9pBZLEXLf6adA/e554XL+xNiuaNFxWTtOPz1pzRRboqosKBuQNdBUXz1Dh3+eqBNleJWShsx5U7uH7f9dQbdW74+g1tc9T58BOVCa675LrWm0DteI3gjoAz9TMsqy4j3BwtNLYFn3RcftI224S0pzTaQqOPSbpuSL6NR3kDmFmL0Gz7KW3z0Fa1Ft0JolLAu4lb1Pb8Z4laAJpzX0OUqnmGuMDvJdrVukjnvnJ7Lp2OM154YW/Sad54EkV7dNrCbx6jrdZGsxqfy4XR6f0mzXHz2s2c+uUp7v7x3Sw2FmPiP1GZoFFvMFGZaHu7mH9uHillqxxB0psADVj14lWt72aPzXK6rgJu2l/vyphJOq4YIfmF1xTVpGNsAtILehURqLSAJShh3YMST7PvqjkP26qeR/m8byMqJxA6ztPXJGHu+v9N1jxMF9YCKpTXSX9V11yCtIM9TcrqebobeAvwlJTy1WWMOQiK5o0n0WmPzqTzd1y5g+XV5TGLvYxAoestwr72B1/3wZjrp3a81spska1MkYhOuh0VmefssSTrtAgB+YQX4qK62XFMWrAzpPuCXua8FoDrUWJcRbk40jopmRa9Kbr62qYVbs7Hnn/asfp4E9P6v7553euBnwDnUMylFeLOvvF+9zTKstj/Cvgsye+ES4YyOu/k9V8nLRyuzBAz5dGV/lcW5rUXFhf4TO0zNGSDbz7xTY6cPAKokr4SSb1Rb5t/2ptA2neb1m5iz9E9mRUi8x6XTl7hDYkLueuYtGvYbhOzc1HeP7vAuGaFKGBp1lIHJfJJoqut6pD4/ZpVHkPjeNf80zY02W8V+s9E91fV7pibm/ewgnyblJLejrzfPYtShF1K+YAQYnUZYw2aMrIsurVaXceUsfMTsjcUBauDlqtFCEFd1mnIRiuIury6nGqlmlquICue4PpuauUU92++P/O55z0uG5fwZol9WpGsLMzORQean2VVQrQtVl1mV2eaVJu/F8BbgRtTxnL58/XPeUTS9L9jHRPifpMJUELeMI7N20Q7acys7zwAQsr21+mOBlLCfleSK0YIMQ1MA6xaterSJ554ovA1+luvuzuK+tiLnN/NnNJcROHjIc8sPMOnvv0p6rLORGWCqqhyun46KtAlqlx3yXWsevGqJfHn0D1ZpWrtz5PYjsr5rqPErkq0OSmvPx9Uzrj+7IPAZ5pjriCevZJ3vua8qqhcczNwbJ7fiQW9E5Xvrt80vMXeDUKIQ1LKDVnH9S14KqXcifpTZsOGDYVXk0E3ji5KJ1Zr0WOKkuQGamXd1BdoyMi6ajQavPfS9wKw++juVnqimdNelG4XvN7jyj/Pcj3kEak5on9uZg54EX8+xDNNjhK5ZZ4HPtq8lin0esFwzXcncCfRPsU091LomI/tyw9pXzimidI6J1HNPOxn4zovbcy07zywhLJi+lV4a9QwxTLJDdTKupGN2LlCqPzjTWs3ZW54yrq2q1xA3qCyOY6eb2/EP69gh6TngSeNWUUV+HJtTjLHTPPnm59tRDWk1mmLDxjXXSASX3vsEHgYMHvnvoNkV46roXZgHeNaAF3oIKxOu0wrDpY2Zt7rjSdLRtjLCGqWwaCsytrxWisjJK/F7BLLtGJbOutGI6XktkO3sfvobsLNYa7G1mnX7iSobO+CFYhW+mX5b20h+QQ7ID0PPGlMUBuUTCs2MM7L48+3LdUjRGmLJgIlxjtxi/KMdfxzCfM3XSnLaG+onYW5sAmiNwyd6+595b2grHTH/4n62/ISIcSTwJ9KKXeVMbam/K3jxenEHVTGQqA37+j87T1H93D/5vszx3OJ5bY3bksstrX3kb1847FvtPzpWuSL5Ivr+517dq7l2lmoL7SeQdGgsnkPjbqaTzldlpK2zOcR7CnS88BNAtyWd1bwFsf87PN0nZdlKNdMAyWeleavO4n82rYobwTuNcbdSPRMJomCtTpdEVSgtmhD7ZBIvHVsQWfw2OUaAucInuKUlRXzn8oYJ4te+J2LUNQdVFZcIHw85Ew9qgzourZLxNKE1JVSabbEA9RGIsc9Ze2QPV0/TUVUWq6dhmwwefZkR6mQdo0ZKdVO2Ky3trTKk2955SRrXupyARQR7E2k53a37oz8/mA9B9PKnWjOyV5g7GOmUS6eeZSFfjtRNkoDJc6mKE83f9Xpl2uIt8vTC4Tpnqsm3GdaIDkgvrDZOelZZRXssX3+eh6WjCtmGEgTyqTOQGXEBYLVAcuqy1qC67p20gKyee1moN194yq2ZdZpmTx7kj+6+49i+eJZC1XMujb89RUqzD83D7QvzvZzs5+PvSP4TOMMlUqFD7z2A60SB65nGj4etlxLC4sLsYbc/7ooePV5DYRwpd71UrDzEhIJrFkKwPRBh7hdPKDE7w7iIu0qNzxtnKMbdZv58TpTRqJE/o9pT3fMiktkPae0Z+PadVukYcf44oW9AElWZVZnIFdRrDzEygpsDhN97EmFskzfNMDDTz0c29yUVed9zXlrYvdq1z6fPTYb+95cLCYqE61NTGmFzfLUczd3BDdoQIPWxqmkN6HJsyejEgw0OPUvp1pz/7t/qnDTb1ZZVnVt6unUwt5OeVblJHFLWdLu8w9I3iyl72EW1eRCZ9Kkocf7V+vzt6JK8OqiX7aohmT7yTsNdNpje598XrywF8RlVSZZ5lMrpzKLYtmYJXFtsf38Wz7vPMf1JmFXaPzCoS8AqjHIWRNntcTQFuznF59vLSB6HnqhsIXbTIHU49nBWTsrxvzZfm6mVW0Ltu2S0W8FC/UFZsIZZoKZ2HOdf26+5Q6qiArnv+D81vmHTy7n0Z/tYM1LkzYIFRGiMqxKs2XeepSACeJiLFEibbbUy7KEQ6Kdn2luJX28vRgsI8qU0Ra9axNS0gJj32PSXJOwx/Y++bx4YS+BNBdNWlEsG9tHXW/Uc52X9CZhVmjUSGQrmGmK5kRlgnpdlQrYdWQXu4/uZrG+2KpNs6K6opVVM3tslsMnD3PwxMG2+dkLX9pbzeTZkwghqEhV++bUL085W+fZ96gXPR2c3f/Yfg7MHWhbCFZUV8RcSWbK5pqXlmHp1VDZJdrlcRolirro1kLz+xmSLXm7ZR5Eoq7zy7X1foa41Z61AAUUq0+jx3Nl4ySNlecNp9MNRa6xuy11PB54YS+BtKBgVpqmq4pjXSqBrVQqCClypXe6BFWL8M5DO2NpjIL2tnTXrLuG2w7d1uoKBVGRL7OaZLA6aNVcb9CgQnJBsqR7M63zhlR9Yz/w2g/wF9/5i2iHq6OXrHmPa85bww1fv4EHTzwYW1wgektw/ZmkxziKWJVarLSo64yPI0SWdgPYj7Iykyz5kHjLPJrnC+BlqH+ij2fMJY3NzV+zUhRdG7NM0gQ8a4EJ6dyFYo/dqVtnvPDCXhJpNVCSRD+piqOriUUnQVd9zq4ju2g04puP7LLEm9ZuitVKl8iYxW67eLSL44pXXNHmBslzb0BrHIHg6MmjnGkogRMItqzbAtDWuNt0VR39f4+2rjdRmXB2W8qfe1/UqgyJgo0V4Arg36GyUTQ6bzvNPxzQ3iMUlLg/aY21nGIdnFyFufIcm1YvphNRDfCld/uLF/Y+kCT6aVUcbb+0HicPtqVs7yjVaYy2C8XlH3dVkzQFWot6lv/cvjcg5k9fd8E67n1M5VVLJC8660XOwKrtqgK1EFyz7hpnT9j8C2JIMasyIC5WM83PdUaNzte265Pr9MQ54j1CZ1FNnx/CHej890QuHU3aG4Z5P7ZLyMY8thdBySIBaU8ZjLWwD7o2SStrZnEBIUQr1zvPNnsXenfqnqN7WrsztaVs+9oFAiEEdz56J5NnTzJ96XSif9zE9QZi1pqpiiqfvfqzuSpUmuOEj4dUqLTcO0dPHnVm+iS5qnT53mqlSqPecLpyrKdFXGgC2q3KNOFMEqu0TUZmx6HbiVIYQaUrrkdt9dc+ek0Vt6inWdn6frSrSLuEXNa4697LxrtQ+snYCvswFBXTWTNb922lLuNZM7XjNWbCmVYudt7AqyngC/WFlqU8e2y2lcVSrVR53ctfxwNzD/DgiQd58MSDAExfOt0aK23BswU6fDxsBTIbssH1+67ngfc8kNnuz/5uxUQU7Nx48cbWZinTh5/mqqodryGa/TX1r26SRNEWZVd1RZ29YmanxJ6O9Zn9+5CoVd3zqBrl96AEWABvAH4V+BrK2q8An3NcJyTdytb3M4MS9bRyuZ1a1P3eLOQ3J+VlbIW97KJinVr/88/NtwTRzkHPClCa15w9NttmlTdkgwd/+iDB6oDPv+XzsayQmXAmNtauw7uYvnQ6c8FL2uFaFdVop2mjEXO5uJ6V6xr2QmDn0QOpi4XOddcB4OQ/0xC3KJqibKf3zaK6EC00v99NZy6LAGWB67IFf0uUkgiqkNdyVN+atHrtAdlW9hRK2POkCBa1qPtdOnd8S/V2wtgKe5lFxbqx/tNy0NMClDsP7WxZ+ssqy1ruCZuv/vCr3POTe2J56wAbL97Y8mkDHD51mN/769/jxC9OtKxve8FLus+plVP88dQf88lvfxIpJSsmVjgDmXoccxHKSpV07URNerb5/0wDskXRPgbiAc4zpPut08rRbiEqW+Dyp59BBVvtse0x83Rk6sa/neXD126eLB9+N9cxr9dJHGA8rfyxFfa0bJWidGL9m5ZvUg66HaA0z71+3/WttMSF+kLrO4FobcyRzf9cc5q+dJq7f3w3dz56JwCLjUXu/OGdre8rVNp2yybdZ+14jVu+ewugslO0i8S1S3Xy7El2H93dWoRcDbE7xUzxzDiSfGJnpgpC3GKXJPuts6xLs2zBBMrlYtblcY1tpldWgA8BtzTHOIDy36eJe9G/31n3YO6ObaAKit2PKiM8hSpAphed6S6uowmIYgYV8jXHHl8rf2yFHdKtvyIUtf5dlq+Zlpe16ISPh23pi5rl1eX85VV/yZGTR2I7Q/WczAXlxtffyD0/vqfNhQO0ygGYfv+0Coytkr+Slt/btUtVCNGWzVJ2bENn25jpnMVwpQpOoYRrFiVejxI1uDA3DUF+/3dIVFpX83LgJO0+8ZAoqNoAPtE8XiZcI889hrgXNtfGK3v8edqLhJ1BxQyuIqr1rt8Kk8Q9JJ8lrt9Qrm8eewPpi1mRsUePsRb2sihi/beCognuDnPMpHGeWXgGBAgpYhUPtVDqIKjdHMO1oOy4cge7Du/i0MlDsWqOLms/6T6fWXgmVpfFruQ49+wctx++nbqsI6RoLSISyfoL1nf+4I1nmrQRyv32lGXJhST74EFZ7nohdJX3DcjOsNH/69K6mncTWeJVorTISdpLDJgUqUOUp+2cFvUKbndVgOrSZNeWOYGy1E32kizsAfkzcuaJFrZOUlLTxu4X/XENeWEviTzWfystsGndmht/8vKR/R/h5m/d3Pr5Xa9+F1/+wZdjW+ddc9ILirbOF+ozQV08AAAgAElEQVTxiocTlQle+X+8kkd/9mhLeAXtu17t+6wdr/Hp2qdbP7sqOdaO11rXqYhKK8BpHptFUnDabsRx9UVXZzbbzrbkApKFeY6oPrnGrsOSJ8NGH2tavhXg5ygX0CngbqK0yM2014/RNMhnwWpC2gPDeq76O3Pj1YxjXH2PN6Pa62mubf5q13pPokgMIKB4iYRhyp/vn2uorEYbVwJ/gTIxviil/PMyxh01WkHRZqZL0q7NJGrHa3zy25+Mffb0vzyd+bbgSoVsyAanfhlVPJR1yQuXv5Bl1WWtlMgt67bEKkkmlSY23UJJpQCSKknmWdTSgtN2sbOvPvpVlleXc90l16V0mgpIF4g0Ya6i/tloyzHJok3KsHFVaVxhjL0HtXDY3YYAziJKi6wQVW0s6o4x738C9cZRJypiZj6bmZQxp1BFwkBZ6tcSt8zz+Nj1OHnm3YlQdxJf6BUh/XINdS3sQogqcCtqa9yTwENCiL+RUj7S7dijhu2jLiLqoERMyrjFtvHijZlvC1r8TD96hQrn/8r5sZZ4B08eZFllmVMU00oTr5hYwcLiApVKhc9e/dnMnHVXGmPWfSe5V/Qz1YuWdh8dPnk49vYSJ49AJAkzqE5EoKzq88muwxIQT3Hcg7tK4xzKQnd1G7Jb5UG02UmLsv4c8m+uMq95GvUGob+bbP6KYwx9DdMCXWN8N022oHfCMAl1UQL65Roqw2K/DPixlPIxACHEXwNvB7ywW3STiVM7XuPBnz7YyngRQvDh13+YNeetaaunYmPWhW9I9bawYmJFy799+ORhDp48SEM2WGwssurFq5wBW7txRVpWjznvNGtbP5c00oLTZibMriO7ONM4g0Ty4IkHedMdb0ppIVhEIALi/yB1M2qdobE+Y6wp4imOi7S7bnTmi9nkYweqqFjSnKdwV2LM88qfdE1znKwxQsY1ONkZ/XMNlSHsLweOGz8/Cby2hHFHkqKZOLpMwBePfLGV3lgVVT735s+x5rw1ufLnbVeI7pCkXSLVSpVllWWtMgQu94jduGLXkV2xRhdJxbZMa9us9V4k7z9rQTSfqa47DziDp51tJLP/QYZEwcUG8P7mcdMkW8pJnZns420XkK4IaXdQMueWV3DNa2H83iU2rjGwjgsYvuDksNOfN44yhN21d7stwiOEaL2brVq1qoTL9od+15Mxrwe0+cZB+cfnn5tP3Ojjwl5QzOYaNOC6S65j1YtXJd6n2bhCIFoB0KzrBqsDqpVqrNb74ZOHM7OCsubvYtPaTew+upvTdeWTLtJCMBv7H6SZ6tcAtjZ/n9Rkw2Wt7WyeV0f52c3eq5Dc3CIJsxm1K/hrxgnMAmX3AfbCHFhjPAP8ZvNezbmWaYGO52aiXlCGsD8JrDR+vhAVSYkhpdxJM2F3w4YNWX26hoJe1JPJ2wx6eXU5m9dubvONA60epX/2zT9LrV+edi27ANn6C9a30iRd55uNK8y2d9VKlbln59h5aKezxPDUyimuvujq1kaoM40zPHTiIZUV00FWUBpTK6cKtxDs7M9zChVWej+RuNdxl+aF9jTHWvNc7dsG9QYQkp2dk4Srk5NZksBcJMy+pnnqx0yiFiCdDWTOtSwLdHw3E/WCMoT9IeAiIcQrgJ8C7wTeVcK4A6cX9WTyNoM2rU4tplf9+lWc/yvnt/LTtWtG1y/PE+zUpBUgSzrfLrurK0nuPLyzZclXRIW3vuqt3Pj6G1upjvt+tC/2HHSqY9GsoDyUU3IgD9PAT4BPogRyBe2t2yZxFxK7nPYKjoLs7BzThz7b/L0Ows4aY+ogqGmFB0SLhG2x29fV19DXDYmneFYSzumUPBui0ubnFwCbroVdSrkohNiKKlFXBXZLKb/f9cyGAJcQdOOacfmbzTHs69nt3Ozr2ccWuRbEC5DZx5gVG59ffJ6bv30zl73ssrZUx8XGYqv4l0RSl3XufPRO/vaHf8vn3vw55p+bj+00rYhKa+5li3oWejHb+8jeVjZR59RQG4lA/TPagRJ7s3VbiNuCP012c+nWrIkLl91K73bgT1DZMXrMCfKlcOrf28/Btp4/YM33Q45zimD7+vNsiEqbn7fubUrJY5dS7gP2ZR64xLCDdlAs6GcTrA5Ug4hmwa49R/fE3AVJQcKsgGiSq8X0bdvXAmJ1YLT/Wx8zefZkTLDvfPRO/ubRv2HFxIpWLZjJsydj2TYmdVnn+n3Xc+vVtyaW2oX2Dkm9xCycdmDuQPMNBTqz/EKijTwCZSFDuxCbrpRJVJZLtfmdzoUHJf6zZM8hJN5Kr05UXoDmXK5JGMeVVZN0DXNBOkp8E9U5GXNMwxblzeTbEGViv52EGcePH37naQbmq70ZcOyk2BcQE8AzjTNtY7h2dtribX627Y3bqB2vtQnk1MoptqzbEutjaldqvOHrN8TKCJxpnOGGr9/QEl8bneq4dd/WVkbMjit3cOTkER752SMceOJALB6w2FjkyMkjzgWokyYi3QSxXYXTfjQ/y9RKM1OliOUXkK9s7n0oITqFsnzrqH921zU/u9NxXmzmxN0uAe2t9PQzr6BcQpvozlUREL+3jaj6OKeb155E+ew7GTskXhUS61ozGWPWyH478XhhL0C3xb42r90c22BUFeldflziB/G3hh1X7ojt5DQF0uxjas9Xu2psHjrxEJfPXh7rUaqpUKFSUW8cOqPlyMkjqrn14oKzbPCuI7tacykSyHRlB3UTxLZ3yL7+QsFbX3WYYn5dk7SMELOyIcAXaS9DsAolwHcTLzZmuykCIhHfgypEFqK28n+VuItkA2r3p71xyQ6k5r23WeMzfZ068EfEM2qK/FnYVSHX487FTyIkCjibbyfe524yksLeqxTFohuMXMHQrF2aaefrDT3mZ3sf2RvbOGQKpJ2/bm4IsndsanQK4/xz87EMk/UXrG/Lf481pTYKWQmiQl9nGme47dBtbZUW0xbJpOygToPYteM15p6dY1l1GWfqZ3j9SsH976kwUTlIfr+uC5f/26ydci/tFRD1TtKgea4Wan1t201hul304rMN+ArwW6jGHJqzUOmWZmB2AZXR0qC4EOu3mQpR+QLd/alTN4hdG2eeYpk1AXELXy+G3uduMnLC3uuWd0U2GBUNhmad72oPZzaBbtDgmYVn2twyc8/OMfPNmVYJX7Nj0Uw4wzce+4az8JfrXmvHa2xeq+qU64CtttgbqMyYicoEQgjO1M8k1oRPWyTNBW2hvsDhk4dzFPZyYxcJm750mmvWn6IitLVbQVm6l+Qe07oCUUrgDbRXOzRFvQK8DVWGQGOKmp23fqp5jrZQ7cXnYuLC/jzxwKyuKaNTHIsIcWjMxezwJIlq5ZjzyWsxB0S1cTpdTO03paL5/qPPyAl72SmKnaLfGuy+nOAOhrpIEj/zs4/e99HYOZ/89idb4nrVr1/F3T++O5YLb5fhnQlmWr1FXYW/7HsyF019nN7Sr0vz1ht1PvyGD/Pz53/OnqN7ONM409a0Q9+f6zp2CYSDJw6yrOquYZOF+fdBa+yN39jHvndLllVU/n9FHAUOoSzUIm4L01LURbvS+L9QyWOunaQ1VN0W/U9yApWPIFHB1reiCm6Zc9qEcs9okbwW1QxbF/e6hqj8QVEhDYgsY0FkpdP89TqiVMsiFnOaCysvtoVvzrWTxWL0GDlhLzdXuTM6fWtwuZBc4md+9pN//knsOx2crdfrsY5IkFyGN8lytueTtGhOrZxi9thsKxDboMEnv/1J3vbKt/HaC1/Lt+a+1da0Iw3zbWL/Y/tpkFzDJgv77wPA38/VuXwWfnu1YONvrOeSlx1CCVdRt0VIZClWiPLDKyhBDVALxRlU0PFi4Fu4t/qbu0J1kbHbm3Opokoy2XOxXTlTxNMt9fGuz7IwBVhvUNJuIYmKEeixzOeQZjGbVv22nPNIOt/O7hmm8ryDZ+SEvagfvBd08tbQyWJQO17j13711/jpL37a+myiMkG9Ee9/qgX9mnXXOC3eJJeLPZ8ii2ZDNtoWFlcP1bT6L+bbRCeLtOutCZTr6KGfnubYqeW889Wmlast004bOWhrf5LI6n8HURByPW7LMiRePXJV87s9zbnYmR+uph0al7/aFGAc3yehx6oBb0Y13tabscz5BLSnddpZM936wbPOL+KnH31GTtihmB+8F3Ty1lB0MTCFd6IywUtf8FLe/Zp3845XvYPZY7OtVnRZ7pUi89n2xm2Ji+amtZtihcps7LcFe/6uRaebRTptoTTHXPNS08rVfvJuGjm4tvbbFRttV0+AW/Cl9SuO8dME0vb/6zeCLWSXGnZdT7t39GY4U7xN6/4G4/iriGIK3fjBwy7PHy9GUtgHTVFB0lkbRQKEpvBWqXL9/3m9s8JiEUG3a8OYi9Pk2ZNsP7CdybMnmXt2rpUtY7qMbr36VnYd3sVZy87i28e/3RL5icoE713/3lawdfuB7cw9OxdrkOHKnNHj6tIERTYzpS2U7Qu/ae0VdVvYlmJIXIDs+jH2Vn89hp1eGBL5tc0OTfb4+nMbl/+/0TzvNuI+/rTAp3k9UG8T4F5c7EBmnShDaBnRxqxO/OAB3o+eHy/sPSLvW4NtueYNECa9FZjt9yqVCusvWF94HnZtGJ3iaO8w3X10N+HmsCW8Zhrkh6Y+xKe+/SlV/11UWqJuZqhMVCao15VgpFWK7MRN1XmsxXQ/dLIJJ6B9c88BotrtzxjjQjzVUVv2OojrErJJIj9+msDNEk97rBjfmamKkF48zL6fgPTFRR9v18JZRDW41m6mota296MXwQv7gLGzNvIGCJPeCsyGGI1Gg637tuYKWLoaaZz65SlO/OIEL3vhy1SuulU24Ew92jlrW8hHTx4FaFWBtHPwZUPyqslX8cP5H7bGtStU1o7XuPlbN3PwxMFWvn1S3Zus5wNFyhd04w9OEqDrUeKm+9VWiCpeL0flrC8Q7cicd4xTQ7k5dED1A7j95h9BbZIyhdVOu3SJ9AJRbfms0rxp1rPu2fo1ooDrBPndP0l4P3pevLAPmG6yeMy3Au1GmTx7kkql0tplWZf1lvimBSvtRho7D+2MbTpaVlmGlDL+WXVZa772fWy8eCPhEyGNeiMm2GYao90426xQWTte4zf/6jfbfPYSye6ju3O91ZhunGIWf0h3/lw7WDlP1JtUYwrtaVTDMXNH5iTJbh593Gdoz+DZSbR4uKjQXo+lSuSq0SwQ3bcrQJtUddJcEP8Y+BTqOQpUkNo+x9MLvLAPmCL++CRhtoXrQ1Mf4jO1z1CXdVZUV7SqUqaJm9lIA4gJOMBFv3oR//jP/4hsqDrqb33lW7nxDTcmbjgCkLK5PalZRsGVxghK1M+aOCtWodIsS2yj3wDyxg6KZykFuBtVhBTPb9fujeVEFrlNFZUGqTGLiiXNKymDZ1fKvCooS3zGugdBe7XJKul+bJf1HNJePEyziHpr0RubhnV36GiUJhh5Ye93B6ROyOOPz+odagrXOSvO4Zvv+WbsvrMKmOlGGtodY/PKyVfyw/kftizsy15+WWJRMoD33/V+zjTUa7hZXCzvpqhgdcBEZSIm7hOVCaSUsWBuNzt4k7EtUijmmglpD5jq8b4PfMk49h2orJGdxmeuwlZacMyUSjuDp0a8RypE1vgEagOT7Q4JidexEc1zPptxj+acAtz+eLM+fYXOdsD2k9EpTTDSwt7r8gL9JM3qdAmXnU2iS+wmiVvMmv6n/S3L/cIXXsi7XvMufv78z5moTLRl7eR9xmZxsfnn5p07cu35PPCeB7j5Wzdz4hcnuPaSa1lz3ppYMDdpYbAX887SJk2LtOiW9YB2i98c7zeJioRNEzWU1gFWW1RrwJuM8e5vfm9n8Gwn/kZwHvAzovLA5qYi11zNdMZO4xC2i6bTVNJBEDIqKZUjLey9KC8wqDcAU7x1O7ra8RpAYukCW3TziKm9Keim37opJqLXXXIdLzrrRcyEM2y8eCPzz807n/GmtZvYc3QPC/UFQPnG7ZK/eZpYf+WdX2n7zHz7sFMlwV0Jsru9DQHFWtSFpJcmmG7+37or0jM+ZolK3C4Q1W233SF6ntrl8zRRPZykeXeabRISD7rOEHfxhMb4+rNOdsD2k4BRSansStiFEP8R9af5G8BlUsqDZUyqLMouLzDINwCzJsueo3u4/fDt7D66u9VYOk+rvfnn5p257q7r6MXLzto59ctTfOHQFwC497F7ufENNzqf8dTKKe7ffH9ss5QQURNsuxJlEezKlGaqJNAKzi7UO79Gu5shj/h1m01jW+n6ekXnOQPsp715BcRTOLvZ4h8QX0T2o9wuO0hu6D3sjE5KZbcW+/eA30fteBg6yi4vMOgCYzqtcLGxqGqi16MOR0l+804WNtu6Ncc48YsTsWOPnjya+Iz1OLqi5TMLz3Dzt1TGRoNGW1GwIs/BXOT0whasDnj4qYejALDs9BpJAp31Zx1Szqt8UvBV15zZlHCcnucM8d6rM83j7TFdApwUPHSVMbAXEdeGLP0Mlor/ejRSKrsSdinlDwCEEFmHDowyywt0KpRp7puirh3bJWNa7KbfW4/Z7cJmL44PP/UwD554sPW97h2a5VLRLhSdeVMRFWeXpiLzMhcNs1BZhQoNGlSocPeP7m71OJ2+dDp7YKBzgQ4o51Xevv5887OQuOAmzdNledpxApcAg7v59izxxh3mQjdDtAFLAOuILyr6GSTN1dML+uZjF0K0HIurVq3KOHo46eQNIM1904lrx5VWaM7HNWaW+yXPfZtpjT/5/37Clx/5Mr9/8e8XEMso88a1WzbtmWYVC2vL7plQ1xBCtAqR6Zr105dO51hMAzoT6LJe5c3rV1HlfKHdXZI2zyT/uytjRZ8bEhffWVRQ19xFaovyFMr639o87xaK1cPplNFIS+wVmcIuhNhPvDOA5mNSyq/mvZCUcifNnK4NGzbkbdE+dBR9A0hz39gNJWbCGWaCmdwbb8yf81yvDGrHa9zy3Vs4XT/NLd+9hXe86h25x3ctjFmLW9HFz7zGnY/eGXu72PvIXtactybHeN0IdNKrfBEh0tefRVV4vJ32+u36uB1EGTZp45r3pKtPugTYFF+avzcbd7hEeZ5og5NZD8cuy1CW/3qpuHUGR6awSymv6MdEljJpFmCa+8ZuKLH/sf0cmDvQVVDWdb0yXUGdLkYae1Eyx3OVC+hkodLXmDx7Mibs6y5Yx0w408rVTx+vTF9rJ0I0RZRjnuS+0CUGTqOs7zUZ45rW8xmi3Hbzezt/3ywd7MqD1+PZ1nincYo8hHi3Tjojne7YD7IsyjT3jWsnZhErO6kxh+2qKdMVVPZiFKwOqFaq1Ouqhvyeo3tiOelZvVHTFiXtJtr7yF7WXbCOW757i5p30//ev0YsIcWEyCy3m+a+cI2rPw8SrjHbPBaUuH+B+NuAKb414tZ6Un57Hp/+LPEFI22OWQSMSlpir6h0c7IQ4veEEE+i/nS+JoS4p5xpLR1cFqXN1Moptr1xW2ru+IqJFVRFNbfYaFG+6f6buHz28lZOu329tPnlmbtrvvdtuo8rXnFFK0i5sKgsd3MOafPefmB769iplVNsWbcF0SyItdhYjM1jauUUO67cweWvuLy1czXr/k2mL53mnj+4h3NWnNMqZFYRFa74d1eUmK6qXQ5J9x+gBKhKvjz4y4GbgD8CfhfVUcll5dvjTjbP/b9Rja53ko8FokXBJKS9dHASUyj3i+1Trzb/39O8pzc1v7upOdca2c/Pda37gI/j3TBuus2K+QqqXfrYUkaufCdB2bwuijyuoE7SIfVGJu3W2P9P2ZZ70huC3sykSxfbVrneIHVg7kCrUmVRF419r9nuo7w+8TxuliL+5ZB4PfOvAmcRpTimjRsS5ZU3ULVZbPfMJpTPvm58VsG92AREueoCtXCYpD0jc25zxjX1zlhdOlgHaBeIShnkCciPRlpir/CumC4pI1e+k92seUU5jyuok7nH3EjNEgRp9dTDx8NYcw37WF2DRlrFqJIEvOiiVOxei/jEZ4myRtLcLFqIzIAitAtjQLyeeZ5x9TiTRLVhaP5qnzcFfA4l+nWUBCTVhbEzXm4gWijyLmj6WF1rfoLoDUAHaM3FSJcNziPuPjMmCS/sJWBmuZg/5yFPizgXDz/1MGvOW8PLXviyWJXFpPklfd9Nnr+23MMnQs7Uzzjrqdtt+lz1ZsLHw1afVrtyY5KAd7Io5b/XkHw+8Roqv1svRq7iXfbxZps5U+BMH7fOiDFzx5PGdZXK/TRKJO3epJpp8m/vtzNeQqLAbt64gSsoa/5+F/HFaCvZgWCfGZOGF/YS6KbUgGmRprWIM9l5aCfvu+t9rZ+vuuiqvu6AtdH+cYHg4acejhXq0tv+AWRdsuFlG7jkgktyB0gBNq/dDLS3+Stn85nL6gvIF5wLiSojClQBrbxuFtslERrnaoHf5Jhb2pingXOABxzn2SUKssbVTKJcNbrcrj4/IPsZZTXd1nwWZanrZ2K2AkwipNzMmNGy/r2wl0A3ueNpdU+Sxtj7yN62n9M2CvWycJkucSCRnGmc4fp91yOlRAhBQzZijTQaNDh48iAPP/VwrPZ6kvVtL5jmOeXgsvpA/QNPK+KlmSQu0OszrhcQr6RoWuyB4/g8fmRzTD2OfZ55n9XmdReJKkWCW9R0OmUdJe47iC8+thVu16HJa1Hrv7va5ZP0pmESUF5mzOhZ/17YS6DbLki67ol2WWSNsfHija2dlPrnJExxrIgK689fz7WXXJu4EJiLANBqWp3kHjLvvSIq1Bt1lU4oK62SBxOVCda+dC0HTx5M9MW7rG97wZw9NlvyAhUSt/puRrVz077nLRnnz6METxfbyiqRkOaS6PR+8gRmQ+JvCtp1tIC653uI15DRC5o+r4G7+YfpQw+IatmEFLeoi7iH9LXL2vAUMmp58UJ3t+knGzZskAcPDlUhyK4pwyouMsbOQztz1UDZfmA7N91/k6rOaHDbW25rO89shC2EaFVjBOUT142rk+Zt1km3ywRDej69677t+IP2wZdXWdP2eevNQCYTwK24g3nDaunZbgVznhC/x8uAQ0RWuQ6+apHPU6nx/ah8eM0fotxIw/hsXAzrn2M7QohDUsoNWcd5i70kyvD3Fhlj+tLpXHVatEX9r4v/Gvvc5b4JH48aWtstOs3G1Wnz1s0wXIuTq7m0vRiYgm26aOaeneP2w7eXXCrBTslz5X0v0h7MM4Wz7DKv3fp60xpgzKKaTN+FEvJlqN2kDzePt9vtuRpq56VMi7rXLKW55sML+4ijxfGj932UB554oPX5ugvWtR0brA5ijbBNzMbVWdfLysCx3UPafZPkogHlEqpWqm0ZNd1jp+QtEPUAdQXzXMK5raS5lGE5hiS7FcyUw/eiYgJmzRhXl6M8fv5NqA1I+jwdC8lz7rCwlOaajRf2MWBq5RRX/tqVHHjiABKJQHDOinOAdjfIrVffytZ9W6nLOhOVCa6+6GrOf8H5uVIw82L6zvVO0CruXbe2O+a6S64rdS4RttX2MO5gXkjxbfx5cY2dlfJnXzfAHVQ0x9a43CyddDmaQgVhi56XxWhlqvQTL+xDQBHfeie+/NrxGnPPzrGsuiwWnHWlaU5fOp3qTikD3aha14epiArXrr/WKdh2B6dVL17Vw9RO02pLErmAuHDqbfzd+mdrKHdQtflzVqaHtu4XaN+xubn5q1nbxZ43JNdy72T+ZVu8O1GbqHQ+/vD6vYcRL+wdUGb6YJEc+E7y5dMsXrN3qOkGKeLr7+RZTK2c4pp113DboduQSBqykSjYZbc3LIYpVkl+9ZDuMyrsIO51ZDeUDmkvHwBxK9xMD3Vl5NxBu2U/DNRQb0t6j8ACo5Cp0k+8sBek7L6ns8dmWznsWUHBTvLl0yzebkWzm2exae0m7jh2R1clEfpHll+923zqkLibZBX58tft8gFJbek0tlWdt5Z7vwnJV8vGxLttTLywF6TMRha14zV2H93d2sRjF8Cy6USI087pVjQ7fRbayjfTITstiWCP2RvxD0kWzDIyKgKKLw5TKPeL6a5wdUVKomgt934SoO7HdDNlxRqWRrpiv/DCXpAyXQO6RgqonZnXrLsmU+A6qY+Sdk43aZqdPIuy33h6NWacgHTh7da/7Foc8ligrk095s8Q3w2KMe4cw7spp+hiGTK89zIYvLAXpEzXgC2MebbMdyLE9jl5rdus4zp5Fr1o3dfrdoD9yXO2/fl5LVB7UTHTN13lEkxfft5A7SAoslgG+MYbcboSdiHEJ4C3op7oT4BrpJTPlDGxYaac4lPZwtgL90Je69Y+LsltUvRZ9CIY2p8Aaz/znEO6t0BnaS8nDHFf/nUof37QwfjDRD8W3qVFVyUFhBC/A/ydlHJRCPHfAaSUH8k6bxRLCpRNr9wLZomBCqqTkKvpROw4UaEiKkgpE+fS65TNLHrrY+833fqMa6hORQvNn5cTCXvSuD74uBToS0kBKeW9xo/fAf5DN+N5InrlXmj1LM3ofJRU3Ms1lzyLkC28ZYtvL8bsnk7F0pWaaPvK0wiJlxPWu4yTLNt+BB/9wtFPyvSxbwH+V4njjTW9ci9o909W5yPTTWTXc7HnkrUI9T64OWzUUK6QPUTlcYuKZZqv3BwnaffpBFElx4PNMcxmHiYhvQ0++qyVfpMp7EKI/cD5jq8+JqX8avOYj6H+Bn8pZZxpmlvjVq1a1dFkx4le5m9PrYx6lqYtHEnFvUC5avS8shahIm8fg3bp5CfJAtUipv3b0J1YhiSLbppgmlXc7O5HNgG9DT6G+KyV/pIp7FLKK9K+F0JsBt4CXC5THPZSyp00y+dt2LCh/7WClyC9dC8UXThcBbxM6zttrLxvH2bZ4Eqlwq1X35paN35wbwGu+uP62iFKvPRfcUF3YhmQLLr6WrZghrSXH2GZnvgAAAouSURBVK6mzMF00UwS+ePLep4BPmulv3SbFXMl8BHgt6SUz5UzJU+/6GThSLK+s6o65llEzLLBjUaDrfu2sua8Nc7jzXk8v/g8s8dmeyzspoU+S1Tb/HTzZ9MNokWsivJQZpUHyMJV+8W+limY+vMFoiYZIuMaetxeuEx81kq/6dbH/lnUFrFvCCEAviOl/MOuZ+UZOElujiTrO0/Ou/48bWyzbHBd1hPdNnYhsd1Hd/eo6iO0uzx+N+XYMkXMvq69zyHpWvrzGWA/StwXyXaBhPTOZdLPdFFPt1kxv17WRMaVYUzTS3NzuKzvsgqZ2WWDV1RXpNaQMQuJ1RvJi0D3hMQF73yUPZMmuGXMw75u6Bg36VpTKGHPW2IAvMtkdPA7TwfIsGaLZAU7bbdLkeBo1rFFygbnLSTWPQFxwdvU/D+kt64F+7pBwfOLvj14l8mo4IV9gPR+K3xnFE21LHJ8nmPz+v77V/kxzeXRS8oQ2qJvD95lMgr4ZtYDZFgtdijuIlo6aYoez9Il785TL+x9IknMvMh5PJ689KWkgCcfWQHDURZ0v3ANG35r/zjghb0PDKsvvdcMs6upPJaSUPqt/eNCZdATGAd0wLAqqj3N3qgdr7H9wHZqx2upn6V9XiauBW24qaGKbeV9Jloob2r+WrO+KzJWPwhpT5/0jCLeYu8D/cjecFnHgNNi7pclPdhG1EXpxJoNceeZD6tlHND7PPWl9AYzunhh7xO99qUnWccuF1C/XEP9S0csA1djiqz5BriFMmQ4i171Ok99WBe08cMLex/pZSAxyTp2fdZPS3ppBIdrwG6iwl0T5G8o7RLKgOHdwdnLPPWQ4VzQxg8v7H2i1+6PJOt4x5U72PvIXjZevBGIyu0uHUu6H4RE1RAFcA35BckllLbgQ7FGGUuVgOFd0MYLL+x9oh/uD1fTat0gI3wiRCBYbCy2FpZtb9xW6vX7RflvPgHtJQO6RQv+OLknfEmCYcELe58YRCDRXEwadVUxUSKXdMplb958eilIIePlnvAlCYYBL+x9YhCBRHMxqVaqMYt9uDNUkundm0+vBCnAuyc8/cYLex/pdyDRXkyAJe9XX1oplODdE55B4GvFeJYcvkyBZ1zpS60YIcTHgbejWrQ8BbxHSnmimzE9niyWRgqlxzM4ui0p8Akp5WuklOuAu4D/WsKcPB7P0DOMJRM8mm5b4/3c+PEFRDs8PD3EuyI8g2WcUjiXJl0HT4UQ/w2V+Pss8KaU46aBaYBVq1Z1e9mxZTwqJnqGm5DxSuFcemS6YoQQ+4UQ33P8/3YAKeXHpJQrgS8BW5PGkVLulFJukFJuOPfcc8u7gzFj6VVM9IweAcpSr+JTOIeTTItdSnlFzrH+B/A14E+7mpEnlaWX7ucZPXwK57DTbVbMRVLKHzV/fBvwaPdT8qSxtComekYXv8N0mOnWx/7nQohXodIdnwD+sPspebLw6X4ejyeNbrNiNpY1EY/HMwh8Y4xRxJcU8HjGFp+2OKr4nqcez9gS4nugjiZe2D2esSXApy2OJt4V4/GMLT5tcVTxwu7xjDU+bXEU8a4YT9fUjtfYfmA7teO+IJTHMwx4i93TFb52jcczfHiL3dMVvnaNxzN8eGH3dIWuXVMVVV+7xuMZErwrxtMVvnaNxzN8eGH3dI2vXePxDBfeFePxeDwjhhd2j2ds8X1LRxXvivF4xhJfAGyU8Ra7xzOWhPgCYKNLKcIuhPiwEEIKIV5Sxngej6fXBPgCYKNL164YIcRK4N8Dc91Px+Px9AdfAGyUKcPH/hngRuCrJYzlGRNqx2s+970ryuh85AuAjSrdNrN+G/BTKeUxIURJU/KMOr6+TLf4wKcnnUxhF0LsB853fPUx4L8Av5PnQkKIaWAaYNWqVQWm6Bk1XPVlvLAXIaQ98OmfnyciU9illFe4PhdCrAFeAWhr/ULgsBDiMinlKcc4O4GdABs2bJDdTNqztNH1ZbTF7uvLFCVAWeraYg8GORnPENKxK0ZK+TBwnv5ZCPE4sEFK+bMS5uUZYXx9mW7xgU9POn6Dkmcg+Poy3eIDn55kShN2KeXqssbyeDweT+f4nacej8czYnhh93g8nhHDC7vH4/GMGF7YPR6PZ8Twwu7xeDwjhpCy/3uFhBBPA0/0/cLpvARYCjn4S2GeS2GOsDTmuRTmCH6eZZI2x38rpTw3a4CBCPswIoQ4KKXcMOh5ZLEU5rkU5ghLY55LYY7g51kmZczRu2I8Ho9nxPDC7vF4PCOGF/aInYOeQE6WwjyXwhxhacxzKcwR/DzLpOs5eh+7x+PxjBjeYvd4PJ4Rwwu7gRDi40KIfxBCHBVC3CuEeNmg5+RCCPEJIcSjzbl+RQhxzqDnZCOE+I9CiO8LIRpCiKHKQhBCXCmE+KEQ4sdCiI8Oej4uhBC7hRBPCSG+N+i5pCGEWCmEuF8I8YPmn/cHBz0nGyHEWUKIB4UQx5pz/LNBzykJIURVCHFECHFXN+N4YY/zCSnla6SU64C7gP866Akl8A3g1VLK1wD/CGwb8HxcfA/4feCBQU/ERAhRBW4FrgIuBv6TEOLiwc7KyV8BVw56EjlYBP5ESvkbwOuA64fweS4Avy2lXAusA64UQrxuwHNK4oPAD7odxAu7gZTy58aPLwCGMgAhpbxXSrnY/PE7qO5VQ4WU8gdSyh8Oeh4OLgN+LKV8TEp5Gvhr4O0DnlMbUsoHgH8e9DyykFKelFIebv7+FyhRevlgZxVHKn7Z/HFZ8/+h+7cthLgQeDPwxW7H8sJuIYT4b0KI48C7GV6L3WQLcPegJ7GEeDlw3Pj5SYZMiJYqQojVwHrgu4OdSTtNF8dR4CngG1LKoZsjsAO4EWh0O9DYCbsQYr8Q4nuO/98OIKX8mJRyJfAlYOuwzrN5zMdQr8JfGtY5DiHC8dnQWW9LDSHErwB7gRusN9+hQEpZb7pYLwQuE0K8etBzMhFCvAV4Skp5qIzxxq41XlJzbgf/A/ga8Kc9nE4iWfMUQmwG3gJcLgeUs1rgWQ4TTwIrjZ8vBE4MaC4jgRBiGUrUvySl/PKg55OGlPIZIUSIil8MU2D6DcDbhBBXA2cBLxJC/D9Syv/cyWBjZ7GnIYS4yPjxbcCjg5pLGkKIK4GPAG+TUj436PksMR4CLhJCvEIIsRx4J/A3A57TkkUIIYBdwA+klJ8e9HxcCCHO1ZljQoh/A1zBkP3bllJuk1Je2Gwx+k7g7zoVdfDCbvPnTVfCPwC/g4pQDyOfBV4IfKOZmvmFQU/IRgjxe0KIJ1Edl78mhLhn0HMCaAadtwL3oAJ9/1tK+f3BzqodIcT/BGrAq4QQTwohrh30nBJ4A/AHwG83/y4ebVqdw8QFwP3Nf9cPoXzsXaUTDjt+56nH4/GMGN5i93g8nhHDC7vH4/GMGF7YPR6PZ8Twwu7xeDwjhhd2j8fjGTG8sHs8Hs+I4YXd4/F4Rgwv7B6PxzNi/P/qrmwP8EUubAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"\n",
"plt.plot(proj[kmeans.labels_ == 1][:,0], proj[kmeans.labels_ == 1][:,1], '.', color='green')\n",
"plt.plot(proj[kmeans.labels_ == 0][:,0], proj[kmeans.labels_ == 0][:,1], '.', color='yellow')\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": 219,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"error rate: 0.04875\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnX+QnMV55789MzurwY5LCd6Y8EORdUm5IMKBQyWzZ0dsZCJjB5slJndJZMuOwfIKU5hcfHsoFJgUyRCb1JmUkPDghBwcTs7B2DE4lzM/VjK4Zgksjg3iwDmcS2xOShBcEQew9udzf7zbOz093f32+749M++883yqpqTdfaffft+3328//fTTTwsiAsMwDFMcSv2uAMMwDBMWFnaGYZiCwcLOMAxTMFjYGYZhCgYLO8MwTMFgYWcYhikYLOwMwzAFg4WdYRimYLCwMwzDFIxKP076xje+kTZu3NiPUzMMwwwsTzzxxItENBZ3XF+EfePGjZibm+vHqRmGYQYWIcQ/+hzHrhiGYZiCwcLOMAxTMFjYGYZhCgYLO8MwTMFgYWcYhikYLOwMwzAFI5iwCyHKQoi/FUJ8LVSZDMPknNlZ4MYbo3+Z3BAyjv0TAJ4B8IaAZTJMcZidBQ4dAiYmgPHxftcmO7OzwDvfCSwsANUq8NBDxbiuAhDEYhdCnArglwH8cYjyGKZwSBG89tro3yJYuIcORaK+vBz9e+hQv2vErBLKFXMzgGkAK4HKY5hiUUQRnJiILPVyOfp3YqLfNWJWyeyKEUJcCOAFInpCCDHhOG43gN0AsGHDhqynZZjBQoqgdFsUQQTHxyP3S5HcSwVBEFG2AoS4EcAHASwBWIfIx/5lIvqA7TtbtmwhzhXDDB1F87EzPUcI8QQRbYk9LquwayedAPBJIrrQdRwLO8MwTHJ8hZ3j2BmGYQpG0LS9RHQIwKGQZTIMwzDJYIudYRimYLCwMwzDFAwWdoZhmILBws4wDFMwWNgZhmEKBgs7wzBMwWBhZxiGKRgs7AWAU2IzDKMSdIES03s4JTbDMDpssQ84RcwGy/QYHvIVDrbYB5wiZoNlFLqdEZKHfIWEhX3A4ZTYBaYXomsa8nEjGnhY2AvA+Di/i4WkF6LLQ75CwsLOMHmlF6LLQ75CwsLOMHmlV6LLQ77CwcLOMHmGRZdJQeZwRyHEOiHEY0KI7wghnhZC/G6IijHJ4Ig1hmEkISz2eQDbiegVIcQIgG8KIf6aiB4NUDbjAUesMQyjktlip4hXVn8cWf2E2yGbiUUPnrjzTrbeGWaYCeJjF0KUATwB4GcA7CeivwlRLuOHGjxRLgN/+qfA0hJb7wwzrARJKUBEy0R0FoBTAWwVQmzWjxFC7BZCzAkh5o4dOxbitMwqMnjihhuAj3wkEnVOMcAww0vQXDFE9DKAQwAuMPztNiLaQkRbxsbGQp6WQSTue/cCu3ZFlnq5zOtNGGZYyeyKEUKMAVgkopeFEDUA5wP4dOaaMang9SYDTLfzwjBDQwgf+08BuGPVz14C8BdE9LUA5TIp4dDnAYRDm5iAZBZ2InoSwNkB6sIw+cVlTYewtDkZFxMQXnnKMHG4rOlQljYn42ICwhttMEwcrt1MQu10ooY2sRuGyQhb7AwTh8uaDmlp8+RI9xiyiWkWdoaJwxVq5BuGlFdhyWu9QjKEE9Ms7ExfSKMnwTQoTUEuazrO0s6rsOS1XqEZwolpFnYmOHG6mUZPgmlQP8Qsr8KS13qFZggnpnnyNAcUKeWu1M1rr43+NV1TmvnGUHOU4QpKgBQW13LgfjQCn3oVgSGcmGaLvc8UbTTsYwSmMaCCGV39sN7i/PD9agS9XKbcb1/+kE1Ms7D3maKNhn10M42eBNOgbolZnHC5hCVkI0gqoL0QvKJZLwMAC3uf6YcB2U3jyVc30+hJMA0KLWZZhStUI8irgBbNehkAWNj7jEsIfQU4iVD34t0fslFvu3AdPx7tdJLkBoQaReRVQIdw8rLfsLDnAJMQ+gpwUqFO8+732z2aeyYmognI5WWAKNrpZNeu5OJue8C+Nz+vAsopR3sOC3tO8RXgpEKd9N3P6+g+V4yPRzucNBqRsC8tZbeWZ2cjy//226OH63Pz8yygQzeM6y8s7DnFV4CTCnXSd1/tOObngeuvjz78jmrs2gXccUcYa1n2psePRx0F4O61das+rw+Hh349g4W9D/i07ySTkGkiTHzfK9lxzM8DKyvAgw8CjzzCljuAzgcZylqWvakUdSFanYV6TiC5VZ+VtOLMQ7/eQkSZPgBOA3AQwDMAngbwibjvnHPOOTSsNJtEtRpRuRz922z2u0bxNJtEO3YQlUpEQFT3er3fteoz3XyQatmjo0RTU9Hv1N9Xq9HfhIgeSi8eTJZrrtej73EDygSAOfLQ5RArT5cA/DYRnQ7gXAAfF0KcEaDcQtKPhY9ZGR+P3C+jo8VfpOiN40EePAiceSbwT//U/n9v1JWSBw8Ct94a/U495+Ki3ap3kWWFa5bGOyyrXPOCj/on+QD4KoBfch3DFrvd6Gk2I2Mmj5a8T93ijjH9Pc/XbMXyIGdmiE44gahSIdq+vfX/yy8PfM5KJfq3VGq36lPUOdX5035/4B50voCnxR5a1DcC+D6AN7iOG2ZhJ7K370F006j4dFr6372ueWaGaPNmoqNH2//fb5rNSFAVUd28OdJcIBJ16SUZGwtwrnqdqNGIzjc6Gon6yEj0Ox9CuENYnPuKr7AHSwImhHg9gHsAXEVEPzT8fbcQYk4IMXfs2LFQpx1IxseBvXs7544G0U2jEld/099jr/ngQeDCC4FnnwV27mz9/4Yb1g4xehcy+0M8ueMO4POfX8t4dv/9wLZtwAknAK+9Fh1SqwEHDmQ4h5pZ7aqroutYWIhms1dWgJde8itHukNKpch1c+KJyetia7xJKFLWu7zio/5xHwAjAL4O4D/6HD/sFruNvFnsSY2zrljsMSaw8ftd9YcoGCxgeWpZTSCa59yzJ9B5SqXW/+V5fS12oujYkZGonH40smYzuiFCRP/2u5EPGOiVKwaAAHAngJt9v8PCbicvI920nUxwH/uRIy1xlmJWqxHdfTcRWbwLXfOHGC5Gu0nqqWu1SLsyn7rZbEXASL+62nMkeUD9jk6Zmmqv+9RUb88/4PgKewhXzNsBfBDAdiHEt1c/7wlQ7kAQelSZl5FuGreQb3y+fn3Oa372WeDRR1t+DSCq1MwMAEuwRVf8IZaL0fJ8P/AAsHs3MDYWhZhfemn0/y9+MeO5SIl+GRmJ/pUk8dtxdMpw4KP+oT95ttiTWMz9dp3YrN+0dVLLS1pO1+6FhwnccR+64g/pAbbGp1vZcsK2Wk13wxuNaGFCo9H7IaI6+hgd7f/QdMBAP6JifD95Ffak4tSvUW2jQbR1a+Qq1euatk7y2kulSEeTvvNTU621MkHvxdGjkW98bCxyv+zZE/1/Zsb+nVD+kF6Knq3xycgbk4inqZ+++Clt52Ar26c+efE3DiAs7ClIKoo+HUHoNtxotBuiel3TWs71ervrdmQkmV99dLTdOO7rO5umM9Dp9XDM1PhsK1BDnUeIzt44bYPt9/B1SGBhT0Gatul6D7rR1nfs6BR2vey0hpw0cmXwxY4dfmXoWjEw82GuG2Xr5btlbZoaSzeGhHpnMTLScos0GpwyIOewsKck5Hubta2b6qJb7JOTYUcDMhJOirsaltjtDqynI/SuxGamrIc6saHegG5ZwfI8jUZ76OHUVPoGyxZ7T2BhzwFZJzJt31XnvlzfTzui1hN+TU113+XUc13w6XX1iwptlbr86jaxD0WaB22ri9pRsO+8q7Cwe9JtKzFt+Vk0JKtI6t/PYsj5knTyNfNzS+t3C9n7xPnVQ/jYXNeRZGjm6oTYUu8ZvsI+1PnYs6SI9k1LnXbfA3UDjXIZ+P73o3P6lJV160s9tTgQbg8JE7Oz0W5ytBqqXS63zmG6z0FSe6dNZB9yhyL1IVcq0UO+807zwwuZz1w2kJWVKL3A+ee3756SZKuuvO6zOuz4qH/oT14s9qyhgd02UmSk2+hosnjypN/xrUsoY1Evxzb5arvPfZmn66ZLRA1ntIUghownDTlasc1FaMnRejuBUlzArph40gp0L4UlybnU66lUolj3JGlEuk3S0bwrMEU9PqRr15RI8sX7AvXkJsEzXejUVOckqozFB9wLe3wFNI3QxvnYpajrdU0abcOdgBUWdk/Stu+keYyyTGb6ipiqD9L67bbbM8l1uYxO03yhSw98jkmKLXfYN34ugLVsEzyfi/CNJ+3FUDLugdfr7bs6AZGFkcY6YZ99ByzsXaTZTLYq2qetut4XXxGT5+nVbmm2UbjNqFMXMVUq5vU2Sa3xkKMnU+6wc9Gk4/C0ll2YBE9dOuy6UN8J1W4PJX0ejt6ByYft6xvkeHgnLOwp8LU+k7a9uON9jRTfCD1fH7tPZ5LkHXRFy+mjiXI5/jifxZAhDTxTIslrK3VaLgVYfaULXrmcbBNZ000wCa3vzUgzhFQfjtzkQ84L6P70rVvbRzm6e8lVL7bYrbCwJyRJe0ra9uKO9+0oktbR9R7FGYHyb3Ldio8bxBUWqfv/bZqWRqtCuWRNucN+odKk+UoKobEJsfSxh/Ah+fSCPh2C7021PUTZ6amNJY2/0nZu9rmvwcKekKRWeNK2ltXqlN+fno5fnOSD63pdvnrXaFz+zbZjm9S1yUl3Tiu1TL2ed106Q69s2kx/tPcoHd6XYZs87YHMzLRGOapOAUTvXp/iYcctPDL9nOYafHxzeg8d57N3+dj0eQHdvSQbS6MRJosjW/BtsLAnJM5K7cX544TftJ4ky/niLHaTrz6uA2w0WsacayRgGr2bDDz1O++qztBC9QRaQIUewnZ6BSfQcjnFzkjaxc/ta9IJJ0Q/nnJK6x5v3Bj97ktf0iJkkvqo4hYe+dTX1ThMkTZ6PUyia8va5utjk+eWD05tLDt2hPGVs8+9DRb2FCTxT/uUlcV/rR6jv58hV7TbJjqlZa3fi2bTbYj5jgT0v7k21pH1fGXTZloqRbObryDDzkhaRT77prpxwyU57yf97u8oe7plTCKeVqDihkiuBuvqoeNClHx9bPJ4vR7Sok+64bbPvRxieirsAG4H8AKAwz7H51XYicIYCL7+a5fbxeRnDmmx+9Zdd/3EuU7TXrvXjmlHjtDL52ynVxVRXxqtrW2Tl/Yij93b7Jg0FaJ1v+Xvr0adFuHZOExulzQCpU9YViqtYaXsYV29vU38fdw4uo8trtHqIwc5fBMi+nd6Op3/0hZ9M4R+914L+zYA/7YIwh7CQEhrtbqOcbXxkMQFPvjkjUkzWokbCRBR2J2RlIrYij355Pbfn1dt0mI1Q+Pwyd5mqqdskGrqTTWXuj4REnO9Hb/zaVQ+vkL9vtTrnfuz+lolWa2jgtJzVwyAjUUQdqLuzmmlsdhDttu4a1MFVo3Ik9oRetOdJHXrzk7R9mL1T7VK9JmLUzYOl0vF9R19wlIOl2RseKkUlZk2f3OSxmZ7QDZrpdlsT/KfxI+Y1ToqKLkTdgC7AcwBmNuwYUOqixqk0VdIH3vIOrk6HBl1I98ZqR36aN83JDk4IXZG8ixW7dSMfUjSB2QbCtmiZkwPSx/WTE+3rHg9esXlWrHVyyWSaa0VV5J/F2yxG8mdsKufNBb7ED/LYLgMK9WHr46as+6brJO1w+s6q5V48b6mvQ9J2hh1P7cpkF8v0+Tz0h+gGnkiBNG2bZ1CT476ujbP1YnrAHwebKORLDlY7htL7/EV9oFJ28vZQdOhpr1Vs8Sq6XfVLK4qQkT/7toVfZJmq9VT7royz8b9TU0hHCprrqnCy7/4TtD8An5iXRUfuekhPPymcbzjHcAllwAHDqwed+OhVmM8fjxKtWurjHph5TLw0Y8CZ58NXHVV+4PQGzhgfljq797/fuAb34i+QwQ8/HDrvPPz7tS6Tz0FfOxjreMnJ4HpafN1zM5GKYXL5VYd9NzNvvmpZf7nO+4Abr65/T7oqYhdZabNhz0s+Ki/zwdDYrH3y1AwBR34fMe2zkT9vx51o1rsaRYO2s6dxm2qliMDQbrVBp67rBX1sogyXTcShUF2hMlLl4gtDlzFdWFxUTOuCU/5OzVsUfdjyxWupiGXvnnujh3m+quulDQbaqvXpY9UQsW6DxHocVTMnwM4CmARwPMALnUdn3bytN+jrzSdS4g6N5v+WVtVXEJqmsvbscOsEdbwQ0tdZZi06laNW6Nj+5u+SDL1Lkt6Pt5Nm4j27m2rxM5NTXoVNVpAmV5Fjc5F0z4365sfPe3kpO9EjBRtNVpGF1GTKOub5zYa5slafYfzpOLrmlsImZpzSOipsCf95D0qxkbSyfhQo4x6vTOyzcfF6Tq/LaRSHi//pgt7nNtTjczTdcPn+65rkCHc28UMHRab6fH7jrYLtqUu76rO0NI6JR/vunWtiik35sgRoivOadJ1I/U1Ua/ZwuTTCrYPatk2K9l0jBRk1T/vikRRwy/1iRbZQagNb2QknS/cFQ0Ud29MI5sh9K1LWNi7QJzFqbe3UFFZcRa7rV4u943tfVP/1ePK47RMt659DL0khunoaCTqr+AEWi5pSdM1f4kaQv0kWitWOwLWlQeTOEy+2aTnLqvTzk3Nto050qSusVbeFqPuMxxTRTpumKeWp98fuTDKtMDIp4NLK8amRjrkFj4Le5dIYhmr71eaVdX6yNwnJYjJ+pa5b2zGku190a9VP49pox/ViIyLpPF1K6vnfhKbaQEGkdb8Jaqn4SQcoeffsr1TtbXokaRh8raNOW76lYxWpe4mkeKs51mIWzk6NRU1PJ+JElmeft7JSXcj6WZMuSsKaEh98izsPcTVtl1JsUzEia3tO66JSl0fXMaf1A/d3apb7DbhdrmJfTqKOB/8KaUjdLC0nZZGWyK9UKnRs7/X7i9Rjd7tYoYWRjRRHxkhetvb2k6SNEzetjHHa0gxEaOm87VNduiukDhLOI3v0NYZ+Mxuu643jdXOFnsHLOw9JKkv26ccdU1H4olCpSyTPuiukWazPchjZKR9Pk5f/zI1FYU/J62ffo/0Tm9y0j0f2WxGaXuX1rWL9HFUqVHeYx1FHRaboyyQvqa4J8fubdLnN9VpYrS5Vp3bylO0LC+iVIpE2uUf1v1sAC1D0ApAK6VS9DFNdvjgK7qm7yVxufj4ydMKMvvY22Bh7zG29uYzYpbf0wMI5GRhWuNEirAexiitY7VMNchDT0NicvH4LCS0XZtqnUs31fS0ZwShYibPl2tr29a9gLG2c6kjhsfv8zTFk4hGs0lLozVaVCJozkWTfoTRdiGWN8jlylButvzuIgQ9/WNb6Wmcnk3YfWNkfa49rah2010zZLCw5whf0dff/RAJv5rNzggVGfViW5wo3Swm8dY7H9UgTXJt+sLKHTva5/pUd5Ba/uP3HaW5t11Or71+jP5D5W7ajz30zxijHSMz2UbqSa3KeivmfQFlurZSp2vKSvZHdQbZ5R/WLPYV5bOIEi2uWu/LSLhhRZLryWJRh64L48RX2Adm5ekgY1skpy8IfOmlaPGdvlrz0KFWOT6oKzUPHepcUbq8HP2rruAdH28/t6zfiSdG9VJXeqqLH6+/vr2e8ri4awNaixCrVeCss4D7749+TwS84Q2dK1EB4J3//iQsLOxHqbQfyyvAF3EJrhAH8LFLgdc9dhCP/ehKnI8HsHn+Gfybi68Evv0AcNJJ8TdNq/D3bj+Eyd3jeOAB4JlngCuvBB5Qi5qYQHldFcvzC1iiKiaun8ATTwCLX6miVFpAqVKOlu4uLbVWiT7ySHSOSiVayTk727pZd94JfOtbwOOPQxCBAJSxAgFgCQIv/vwv4aRbr29vBPpNt13P/Hz0oK6/3q8hhl7WrTcuXjHafXzUP/QnLxZ7v911qlujUjFvJZfE0JEjb9N+B6YdzCqVyFeeJVpHvw4Z/eNTd91Vo44QTAauy1V1eF/ke/feXcl2EeUyLY3WaGK02RHtYlyB6loZajqHaRWoPK7RoKXRGi2tWunSgp9HJcoqqZ87zseXxGfGFvVAAHbFuMlLW1ZXbOuuEdU94Tvxaoof1/WkWo3yRalCr4p70g5PD7uuVMy6ptfXFippc0mbjlkrQ/G9t0XAmCZIbQ9/tVI7NzWNOyq9e30Cn7XrZqkhSJOTLQEul+mxE7bRlzFJP8IoLaJE8xihy9DovAwfv3WSRpTGyum1ZdRvSywHsLDHEHo+J+S8kq+xpZ7TljJETdOtHq+nCtm6tXWM72Sv+jvVj692KEkmlE1i7zKI2zhyhDq2QarVzMtGYx6+qajzqk1aGvHMD+NC86nrqztXAFooV+m3Xt+gJ3+jTp+5uGkOufS1TLplwfTaMsqLJdZnWNhjCNlOspRl+q7PBKWem8m0n4E0Ck110tfAVCpRJ+AKY3Rd5/R0K2lYXMixb5qVRCRZNhrzwExFXVOuRxOY6o21zRzLc9h6evUGyJumPzRT2XqZvjsypbU6XN8z+c7SvkQ+9UtriRXMymdh9yDUM0/T5uJcsXFWs0vINSPQWqfJSXMZ8l3VV8vGrVFRv2NbpSo3AZLnSRLo4URxxSyN1mixvHoSgytmZiZK+PXDvXWa29fsSAOweXO0cfXVqNN51SZVq7QWyqha1kso0cq6WmcZvv5vGYJkWklmy6muLu/vZySLaXWsuoAqSaeTZOSRZBl3Aa18FvYekrT9+Bzvayzpn9FRe6ZWvVyTX17tIHz3a9CNN9WdpK9SlXsby3MkCcs2IhO0fOc7RBddRCulMu0Z+TwdwB56AWN0eN9Mx+GmNADqxOiL9zVpvhLFqC9Wa2vukLl9TaKpKfrXnz59zXo3pvf19X/X650LDU45xTxk0uLdE6e6tJ3f1vjifPO2Rjg5ac4caSOJVZR0GXcB4+dZ2HuMr/WfZD7Lhur2UEfypnQi9XoUMUKbo2yI76rO0FPYTBvXHaVms7WJjslotLlj9Oucnja/x6p2qfOFvu98LKpKn302kRC0JMr0FVxET2IznVw62nFvTWkAOgx7TRCeu6y+ZpHP7WstQloB6Eeodqb3NfWAtgaii6BqicskP/JB2YZX5XL6DbJNq0h9J3hqtc66bN3qn+s9ri46aVIk5M1iz+gmYGHPIb7vjAtdRHfujGm7ivj9n01RGOACKrQfl7flZhkdJTr99GQb38trikvZrb5fMg2J7VjXeTreB1Wllc8KQAuo0OfKl3fU3WuOVamwHvp43UhrEdISBB3AlL0MfWhkekj6cEdG3UxOtset2mbHfR+Uis1P5jvBo16j7s9rNJJZ7Pq9ijsuqVDnycceoKPpqbADuADAdwE8B+DquOOHVdiTvDMmms3O0a8sw9p2LWGA/4wxmpxsr8/Wre0hkXpUn+k8+og8LmV3mtWh1vdBqrRF8BbWm/3rXnOsqxXWQx/PhXkzDmd6Xz3E0ZalUd0iSt9tSM3BoPeQSYd+Nj9Z3Ky3rSw5664KuK+PPSl5EuqkBHAN9UzYAZQBfA/AJgBVAN8BcIbrO8Mq7Fk7bN3N6uXKMJior6JGvyrubtMKKe4++zq45vRM9TFNDieNkze+DyaVXv3MlzszPhIlT81rsvDH0aRbMUUHMLU2ueoqg5rN+IT6us9K321Id+c0m+5dxuNusumcaryt7IXjysibqyPPDJLFDmAcwNeVn/cC2Ov6zrAKO1F6g0MaRjLipVSK3DKx5RnET2ZDlJOsPpkaTR4DW1SPXm+fmPW4aze+DxZXjLzGz5X2dJSfNDWvqe/4hUqTXhM1WkSJlsoVumtbw1kGEfnFeOoXKh+Qa0GU6UYmEZC4HttVRgEnJ7vOoPjYAVwC4I+Vnz8I4BbXd4ZZ2JMiDTNVv+Rcmdf7q5moyyOR5fja68eMHgBbObrbVDcmbZi8EMFcpVKly+W1XkcuwZfuJpO/P8l7ZbLwr0adFqH4n0ql+L3/fHs4X9+8C1dcqqk3tvnYbJZ8EvFngtJLYf9Vg7DvMxy3G8AcgLkNGzb04BaEodcuPdN7bVq/IkfRsRFvFhP1rktnrPNnJlSLPUmkne6FGBlJnss9ltVrXFjfmfExrSGrFd12+969vtnK765emOafVuPb1dj5tRPb8klI0kSB2CYybL58V7ytWsb0tDnUMOQLMsj+8x7BrpgAdMMgSeK2MO1JLF20ck9SdQLQZCTGncvXN26bb7Ntu6eiB06oUTGhDT05wvHdQjA1jUb7rLGWuWy51Ipv74iXl5VUH25cOFFSl4r+QPShk0/PrHYUppwRIWHr34teCnsFwN8DeLMyefpzru8MirCHdiHGtV1TFJoqppOT7XnKbYEWPucichuMpu+b5u5GRzs3ttfz02gbBLXpYC/f38P7ZuiwiOLb31WdoVc2Zdx52pFH4TVRo3eUm2vRNPK6373eMgwrl+09sMmHrvdeccO3NBa76n5R6yuzvIWi2Uy+uGNIrftehzu+B8DfrUbHXBN3/KAIe5y4JSWNGCd12fqeK+4Y3fUi9zS2uWH1j5wH0M9RLvfJKFudBV0uR/H8S+tsuXg9UYc86hL31Qd27N6mMV7+yd+w3DSbsJvOqycRm572y9MQ51fXr011w6h1nZ5Od89s9Ui6uGOIrXteoBSIEHNZalnqiNb0/iXpOOLezbh3vdFo1xg91YfJylZH+dJAtaU3kC4em4egp0bX5ta+p+pG2LR+fbpKxAznbNE03zhDCVHUb5xPfgVTzKvqWgmRp0G/NtWazuqGcfkbfYdxXckiNxiwsHeBNK4ZvWOIc6u6vm/7nU0gXe1fvl+6GG/d2irPJNZ6eLWMyNu2zbxOyLbFXdJOMmsn8Pi9R+hgKVp5Kyu3Uh2lpUqVFlGmlXUJe+qYC1Cjac6rNulzpSn6EaJzrS0W0CcfTIKsu11svi3d2s1yw0z++tHR6AHLCZ4sw1Y9bjZpQ+hKFrnBgIW9C6QRI9048Vml6Tqf6Z2z1clVX5sbRXUX6/ohJ1ptCyJNWjMykm4iM+RIiYjorktn2kSdECXwWl79/5JIYfnZxLPRoOPn7aC7zmvQb72+QUvlCi2jtVF1W/jaDTY6AAAdcElEQVSgKpimiQyTiDWbUaeg96Ry9ae+cCmNEKsdimwMctutOP+8C1O6gSSdkM1/OCQ+d19hL+Sep66tILOQdOtGfStJABgdjbagLJWAW25xl2HaihJo/90990TlraxE/6rbVar1PfHE9r1TJyairTiPH4/eMAlRa49SuRUnAJx9dvS7E08ErrqqtReprI+6r6oQrTIXF4FGI9rf9KGHWnWT55flyH1Qgej5qfudfuhD2bbknJ0FJr58JaqIHsJrqKGMZYxiAQRgEWUsUhXP/NgEPnymtrepC30z29lZ4DOfAf7yLzEKYCfux85SqfPmyAseHwcOHmzfCFa/8MXF1nflxe/dC3zlK8B55wEPP9z6+7p10cNRH+r8PHDFFVEd5OaxvjdPbkpbKkV7txJFD2F5udVQkj6Ml16KyltZif596SX7psAm9Iaza1dng0lyjUXFR/1Df7ppsedpXiXr5KuPxW6a19LdNK4Uvjt2JEv8ZfIMqBa7EJGVLg083UjVyzLdC9Uo03PYJH2msn4nl47SraXL6V9Gx+jKk++mW0UU63459tHViHLCyPDRD3+4lQ04NmhGXoQcOpmGLeqQRw1t0tGHMWoyH5PbYWqq02JXj5cWdppFA7plrE+exC20insgWV5Q/VxDtAIWw+qKycszVt/3LCPEOB+7vnepzBElwyPlSN8msOp7Zkr8pdfFtnhSnwObnm4Ph/TdG0E9j08Omzj09vDe99o1V3UfWTewtt0QVUBtn23b3H4zdSdyuVhAPlA1hlT9ju77VuNjpRsljZDq16Y3IvVh9HrCxFXXfltzXWZohT0PzzhtHdK091NOcWuJyX9u0gdbSKWvYaQbj9I43batFQCS9F5kzVsvy1GfxY//ePx9Ut3a1sRe+g1RZ5Xl8loZhihEq9c03UBD7/qD906tpQdeLpXps2+qm0cP+kOKswSS3jyfBUq+1lRWQXd9n33sxRZ2ov4/47TRM0k7A5PFrhtXUtCTWrw+biCXV8D0MY0WXM8paydtGjXdd1/r2cQZ2bWaeS9sawXVOFB9JZfLetYazHOX1WliNEoPvLiaHvgd5WZr9JCmgWd5KeRkra2Htt0H1zFpHma/LbYcMNTC3m/StMGknYE+Uj7llFa2R9WnHudeSVofl3Vv24fVNFpQ6+/qdLIYmy6vR71OtGlTq86mCEJnjnVbBU0il1D0ZA74c9Gk60bq7Ts0JWlcuv8/TYOwPSjbSEF3B8k5BdtoxZe8+Fj7DAt7n0kiSLoY+3QGrnauGolpjTqXMWorW26zt21bu8hXKp16oKdK8Zm0TSLwPjqgJvg67bTWSEe6t2NdMT4nVvLHqBVRJ2dlgrBXPhjdWNuq1bvv9rwwecNs/n9TL2u7uabzuToXWxytnFFniz0TLOwDgq/l6vqeyRJ2JffyLU83wnS3hZp4TP/+9HT0f+le1jsLma7E5aqJu9Y098eGnsXxMxc36YbX1aMNrJNgs9ilD356es3VoiYIe0c52kCbyq3t+IyjB98NnfWt9PSHZxNpfURhupFxVkWtZvYHxqUQ9bm3Q+BHd8HCPiBkGWGa2rm6sE8aSr6uIH1BoNzxTI+6U99VWV8fQ1Wfazz99Pa66hkqpWv31FPTrZjXRyDemhDKuS+/J8V49SKWAVpCieqYXrPKr0adllZzvC+hRL9Tqq9Z6nL0sJZATOkkjBclk5PZ/GJqp6A/FPlxWfRxfi7ZeNQJV9+GyDhhYR8QQo0wVevatrm0S9z0BYG6LsgskL4Wu7pYUbfY1Zh33aBTr8fmszelKA56j0P4c9WbbXBPrKx+LkODAKLPlKfXVqauAHTXeY2OXZ6eu0wTYdNuJ6YkP/qDVfOxNJudiw70hhN3fbYbLYdt8qFljf1lWNgHCV9r0jVxqb9PejreOHHTLX39c8YZLaPTFVatWsjSYFONtWazPYzR5l+v1+11Saq1iXXa5ZdK4wfSfFmqgP81dtC5aNI8Kq2UA+pQyFaubeHR1q1uUTdZ2LpPLK2F7Rq2Sd9g3n3kOXf3sLCvkvPn5I3vfJXqOk2yOE+3pvWP6o7xCVs0LYo0Wfe2IA2Txa7qQqOxmv9lkzIDaVkqmmpUpPdSSQpwPZCdO9es9RWA9lQadE1Z22rPJKrqkEyPPlE7H9X9Ic8vlwPHJe2RnYrv5Eyci0atozrky2tUS6jhcxdhYaeBeE7e+MxXuUL7fBYgmqzpU09tX0Vqm6zVf68Lu5pcTA/1tiF97DK3lXot76pGSb0WUKGXz2nfosjmJUjdwWeJRTXc7JdvatDTp+6g33p9Y22i9jXUaFlYZrylVa0mCzNdlC7SP/mT8WKq1jXU7H1cGGReX8YBCKlkYafuPKd+jQBsVq5uyLlcnr5iajO65Hmnp6MOQN9IQ73HphG+zS2clHqd6ClspoXVTasXRlpxgQvrx8LrRxILwfVA4s5hO17vJW2zx7ZJjLgImjSNWp901f32Np9hnofPA2AJ9kTYEW1k/TSAFQBbfL83qBZ7v597s9luOcftcJa2Y1PfP70MPYX49HR8gIQMX1ZXxWbZr6HZJHrzuiP0ELbTq2gP9r7n1+8OsyeEzfJM0itmaSDq+VzCbqqnvjGGFN001xR3rfrOR3F5pPMs7ES5r2OvhP10AG8BcCiPwk4U9jnlYaSm1iFuT+IQOqOXoc/NSc1w3WP5dz0TpW+MvYnD+2baLHUZeXF42561HycwQy+dHO+Dj71o3xsXqoGYhk16iJGrnqbfu4ZiLleKWifTkFBP6GNZkNV3q6gg+Ap7pnzsRPQMAAghshTTVZKkeo7DlUPchSs/fNLc8WodyuUovffSUnt91DKT5I83oeegf+op4LHHWn9///vj77H8+403dqbiTsvPNa4EaDXJfa22lqx94xNfQql0ANtWDuJruBDrji7gH96+E6cdeRTlpQXghhuA/fvdhZsS4Sd9OEkaSNz5ZXJ8/UHa6mnaOODGGzsT+ZuS/et5zYEoKf/tt0fHqvnOx8eB668HHnkkyvsuBHDWWdHP+j1Ie0+ZdPiof9wHHhY7gN0A5gDMbdiwodsdW9dIOgKIm9hMY8SodTCNrrttGE1PE/3MzyTf09h3vs30vY6/60tFV4O9D++boVqt3QfftnPS2Fj8M8xyE0MMEW0TKlnq6WOx6yMOuT1WXGJ9mf1R3xQ3RGN3XU+OXSbdAqFcMQAeBHDY8LlIOSZW2NXPMMWxu0bnrvmnbpwvBFnfT1NH9K7qDD2FzbRx3VE6vK/dZZLmfM0m0R9dfYSeObl9j9OFSo2e/b27/coLLRxpLAJbKJJKo9Gayfatgy00Sb/Zvgl9XDPooXz6+nUMqVsnmLB7FTLkwu5qrz4Wuz7/FMrY81lbk/RdC90Zyb1IF1Chh7A98psrO1xYOyo9i5buP5+ZoaV12h6n5Srde+qe3odTd8tnn7bXk/56U2y7PhyUYZa2OHhbPbopvnmY7OoTLOw9wqf9xglr2g0lfOa5QruCQndGr2xquUxe1Vwm1jrOzLTi1rdvN293tHnz2gqnhUqNFstVIoBewFiwTtSbFHHwz11Wp0+d3KCVdTVaLpXpNVGjF+/zsJTjempT4nxXwzVN2lrq7FwRpyYAy2q5s8XeXWEHcDGA5wHMA/hnAF/3+V6RhD2E8ZBFYOO+4+sKStqh2CLqfL7b9k4fOUIvn7O9PcKlVmvb4aLD26CIdltuWzXHruaDn3vbHnoBYzSBmWBuL/tFGf6ewBe+NNraYOMP39Kg60bq7Rtt2MrVM0ma3DMmYbfFhWZp3Po8gb7ln++Q0lU++9i7a7En/RRJ2EMZD1lcIq53LrTFrn83ieXutL61sEW5w4XxO0eOkD1hefL7YP2CzwPxLdi3vHp9bUu8BZTpupG6sd8ylluvty8brlTMwznVb247Tr02Ww5onw6tXm/31esxunKCNunGuEMKC3sPCTGyTGOwJAmICOVj17/r40Yyvd9rx27eTMtl6Yqp0XFU21TM2IHFdAZp7kPHgb43V8177ulmcbokmk1aXtey2OXuSc5+S5YTtzepRE0jHCemesSLrGvSBiiP1S32qan2zqhUSraJwJBZ7SzsPSZtG9PbvG+qDrlbkSnLYi9xuWGbzc5t+jqCPI4epbm3XU4vYIzej7vpAPbQK68bi8SbLPqhumJqtdTbHVnnXz2HQ3P7mnQc1bWEXsfFaKcfXL9ZNoFTbtzcvvYt8Zz9ln6Dpqf9NuJIMIIw3osU8wbGDq3Z7Exc5pNZckj97CzsPSRLG1PfD1dEmYqecrufo9dmsz0/lZ58UN/EZ+tWcyCGy13UsRWfJYZddgY+mOZfSyWin/gJohfva/m5V9aZH8bMDNF1Iy23yRIE3SqmOv3gKp7LhhP1W74TqHEjBRtZLHbf8zQanVZ7XEcROjJmQKx/FvYeEmKOKW4NiMqOHe3CvmNH/DkSt9m4cMJV9PBHmVrXtM2myxdv06JuGWW2+VeA6OyzicbRpN8RdbrpV8wn3bw52s7uVdRoQXGbOAcNnhZ7on7LV2DVicyRkfZMkbYGksTHrpeR9OHZOhCfe9mNGOGcwsLeQ7K2C91lEdJiV+s2MtJKges6/q5LV2PAK5W1iJXlshZOaClfnUyVmRxHR6PzJg3pdEXMZcU0/6p+zkWTrkY92o7O8f2J0ei4c9GMm7+NSGs5+5ZpQh8pqBc6OdkZWaNa/b4z9Lo/Lo21k/R+hLp/AxQXz8LeY0K0saQjZJ8Fh7qrx9UZSJF+EubY8oX1ZnNUnbuzaUScu8VlMMYYuKkwzb+qov4qarSIEi2VzZtOpJy/7Q0u61lvDFu3tg+51LzKvjnUTdknB8gKHqS6srAzRNRqs7p4mdw3MlruJByhB9G+HP9VRClxfc7nE4Fj6wxsLhpjRE0GVFfMunXt9+ZqtHY0WjFN5jWb9Nk3RbHlGedvO8lqIdhESg4LJydbETFyH1J5vGm7PZ/62NIKD4jfmogGpq4s7MwazSbRtm3t754pgVezGb3bE5hpT5wF0HFU6ejFYcxRm/vGFS6ZxFXlg+rH/tSnOi32tj1I9R3BazVaKZdpvlKjd69vpp2/NV9oVsvR5lbQhz9TU505Y9LucqTPoOdcHAcZX2HPlLaXGQzGx4ELLoiyqRJF2VXXr4/+pqcN3r8f+HcfuxJVRClxl0ZqwMoyRpcXcNI3vwTgQOb6qBlcZQrfctmc6XZ2tpVJtlIBPvpRYNeu7BlfTzoputb9+4Ezz4zOv7wc/W2uMo6PL+3HflyBkdIyxOhoR/pZsbyMankB/2P6EHAKcMlph3DgqxPZKpY0ta0p57MtdbBatuSqq9pT9I6PRzcjaZ7n8XHg4MFs+aFNJM1pzbTwUf/QH7bY20kyCkwzYrRZvDYD8fH7otjyhfXpwwnj6qNum+fKL9WLeS1bFMrcPsPN1m9aqL08kw5L1IgVdZGRLCcupjS0bys0chFVT5P65B+wK6Z7hHTHJRl9pxmpm0bg8nuh8tykuRe+CzZzM69li2YJnSzId5WaKX2A74RF3Gx2v2k2/VbQDiG+ws6umISorgF1BJuWO+8Ejh+PWnDc6DvNJjT6CHzDhvhRuy9Z7sWuXcAdd8Sf27QZUM8xXejeva2/Z901yfWQbExMRP6jlZXo55UV865Iajn6Vlc33xx9R26DlRcOHWp3GZVK8feV3TZtsLAnJOQOX7Oz0Y5jRNHPlYq7/aYRYtd3sopm2nsh38Gbb452fYs7t8/2hl19r10XGqLnSfNgx8eBW24BPv7xSNRHRyOBNm1LZ2J2tuVjf+SRyLeeF0GcmIiuZ34+6rxuuSV+riGktVUAWNgTEmpbS6DdMBEC+M3fjBe4pBoS950se8KmuRfdeAe7/l7HXWjWjXVND8mnp9q9u3OyU/0ZiPY6VcuQ5X7/+/ndgzRpQ+f9VDvx8deE/rCPvVVOP9ycvvX3OS7pvejGZGhPFg72Ms45RMMwldHNFV/9JM/zBYFBL3zsQoibALwXwAKA7wH4TSJ6OUB/k2uyGmhqOS7DpBvuBV/rVj/O5jZJei9Cjni6WWYHoR66DyEsUNPkDdDuy//oRyN//qD7pXMxEZMzfNTf9gGwA0Bl9f+fBvBpn+8NusXeC7plhOhJu2w7CdmSe2VN5Jf0WF8GZOGgH1kfvh5PKvO3pMnrwOQK9MJiJ6L7lR8fBXBJlvKYFt1yG0rrdn4+mnN78MFo7ky33FUruFRqLSYy1cVnFKCPPkIbVb00qL1JO+TSLVCg01fu4tAhYGkp+r8QwFlnmctV/e7dnnzkqJXe4qP+Ph8A9wH4gM+xbLHH0023YbOZbOejuPDoOB930VygsRmN5SKhjl1FUhB380yWtrrE35UrWdLtSYqiNYA+glAWuxDiQQAnGf50DRF9dfWYawAsAfiCo5zdAHYDwIYNGxJ3QMNGN92G4+PA9dfHR8apVrAr2CLOx51k9JHEsOuHEXjwIHDhhdF17NwJPPpo9P8bbojSE6xZv9K/DWQbcrlunsvSlo4YwD7UknR7koKjVnqPj/q7PgA+BGAWwAm+32GLPR+kTU9gSx7oyuroY7Cpq+RN+zqkKTM0MivkuWjfvm4tq2OaLbFcuC7UZ9s63c/uOo+aECz0EJEt9iCgFykFAFwA4H8BGEvyPRb2wSXtqN2nE9FXybu2vtT3jpCZYrvG6gUcu7dJV5wj87VHOyedV222NthQRaxa9d/E1nVeU+4X/VymbevkzfTNutht/x9PzmamV8L+HIAfAPj26udzPt9jYc8/tvcwjcWepGyZIz0uRYgt8KMrKBe9NFqj28pTa3udLqBM15Tr7Rts9HKhg+tm+kykqAzQTkLDSk+EPe2Hhb1FHg2ZpPN1SQy9uGOTbH3pm0gsM4rgLaJMnxNTyfY6DXDeVBeY1AJnl0nu8RV2TinQR/Ka4iJurksPLUwyNxZ3rGmVvA3fRGKZUSYXyyNV/Ot7duGSmV349HsO4X/+aALf++Y4vvjF7p431QUmnYHnhT6FgYW9j+Q1WCCpniQ53udY35j0numQciIxMYFPjo/jk9EfcCaA/9Sl0wa5wKQB/rlcEMAkRUTWfW/ZsmULzc3N9fy8eSOvFjuQPJQw72GKDFMEhBBPENGW2ONY2HuDTcxY5BiG8cVX2NkV0wNclnnRR77cceUMfiBDAQt7D8irL73b5NnVFIxBEsqheCAMAJT6XYFhQE4Ylsvdjd6YnY2W+s/Oun/n+n1ITB1arkl6U6RQXntt9K/Pje8nA/dAmLSwxd4DehG9YTLGALOB1ivDrSd50gNw8CDwJ5fN4r8deSdoYQHzVMVr9z6EEy9MuXNPXi3jXjyQQRrBFBgW9h7RbV+6zRgz6U6vXEODEBYtk3pddfwQVlYWUMYyKljAN3/vEC6KE3abUObV99btB5LXDm0IYWHvId00ZmwaY/pdLy3pvE8OX3lldB9mViZwDaoYwQIWUUXjuxO4KO7LNqHM81Clmw8krx3aEMLhjj2iX3sZ3HYbcM890Qb2eurdPFvSveLoUeADH4jS7771tVlM4BBmqxO44gvjuCTLtjHqwwCG42azxd51fMMdOVdMj+hHfiU90WCIfR/yQMj8OjMzRCec0EomJhOKtSX1ysKw5V/JY/KjAgHOFZMv+jE6V0fGKyvR74gGe5Qc2iiUrhgAqNVaXoQvfQk4cCBAhYfNPZF339uQwOGOPUK6Y2+4oXcjVDXMcmSkNyGX3SZ0xN4DD0SJx8bGgDvvBC69NPp/sKRevYp1ZRgF9rEXnKK5egfSjcshgEwgOFcMU1hYJ5lhpSe5YoQQNwC4CMAKgBcAfJiIjmQpk2HiYDcuw7jJ6mO/iYjeSkRnAfgagOsC1IlhmLyTx5QJzBqZLHYi+qHy4+sA9N6vM4SwK4LpKwM50TFcZA53FEL8PoBdAP4FwC86jtsNYDcAbNiwIetphxZ+p5i+M2whnANIrCtGCPGgEOKw4XMRABDRNUR0GoAvALjCVg4R3UZEW4hoy9jYWLgrGDI4QR/TdziEM/fEWuxEdL5nWX8G4K8AfCpTjRgneU5DwgwJg5DdbcjJGhXzs0T0v1d/fB+AZ7NXiXHB7xSTCzg0Kddk9bH/gRDiLYjCHf8RwFT2KjFx8DvFMIyLrFEx7w9VEYZh+gCHWBUSTgLGMMMKh1gVFk4CxjDDCodYFRYWdoYZVjhssbCwK4ZhhhUOsSosLOwMM8xwiFUhYVcMkxnOB8Uw+YItdiYTHFjBMPmDLXYmExxYwTD5g4WdyQQHVjBM/mBXDJMJDqxgmPzBws5khgMrGCZfsCuGYRimYLCwM8ywwnGqhYVdMQwzjHCcaqFhi51hhhGOUy00QYRdCPFJIQQJId4YojyGYboMx6kWmsyuGCHEaQB+CcD3s1eHYZiewHGqhSaEj/2zAKYBfDVAWcyQwBv3ZCTEDeQ41cKSdTPr9wH4v0T0HSFEoCoxRYfn7TLCN5CJIdbHLoR4UAhx2PC5CMA1AK7zOZEQYrcQYk4IMXfs2LGs9WYGGJ63ywjfQCaGWIudiM43/V4IcSaANwOQ1vqpAL4lhNhKRP9kKOc2ALcBwJYtWyhLpZnBRs7bSYOT5+0SwjeQiSG1K4aIngLwk/JnIcQ/ANhCRC8GqBdTYHjeLiN8A5kYeIES0xd43i4jfAMZB8GEnYg2hiqLYRiGSQ+vPGUYhikYLOwMwzAFg4WdYRimYLCwMwzDFAwWdoZhmIIhiHq/VkgIcQzAP/b8xG7eCGAQYvAHoZ6DUEdgMOo5CHUEuJ4hcdXxp4loLK6Avgh7HhFCzBHRln7XI45BqOcg1BEYjHoOQh0BrmdIQtSRXTEMwzAFg4WdYRimYLCwt7it3xXwZBDqOQh1BAajnoNQR4DrGZLMdWQfO8MwTMFgi51hGKZgsLArCCFuEEI8KYT4thDifiHEyf2ukwkhxE1CiGdX6/oVIcT6ftdJRwjxq0KIp4UQK0KIXEUhCCEuEEJ8VwjxnBDi6n7Xx4QQ4nYhxAtCiMP9rosLIcRpQoiDQohnVp/3J/pdJx0hxDohxGNCiO+s1vF3+10nG0KIshDib4UQX8tSDgt7OzcR0VuJ6CwAX4Pn7lB94AEAm4norQD+DsDePtfHxGEAvwLg4X5XREUIUQawH8C7AZwB4NeFEGf0t1ZG/iuAC/pdCQ+WAPw2EZ0O4FwAH8/h/ZwHsJ2Ifh7AWQAuEEKc2+c62fgEgGeyFsLCrkBEP1R+fB2AXE5AENH9RLS0+uOjiHavyhVE9AwRfbff9TCwFcBzRPT3RLQA4L8DuKjPdeqAiB4G8P/6XY84iOgoEX1r9f//ikiUTulvrdqhiFdWfxxZ/eTu3RZCnArglwH8cdayWNg1hBC/L4T4AYCdyK/FrvIRAH/d70oMEKcA+IHy8/PImRANKkKIjQDOBvA3/a1JJ6sujm8DeAHAA0SUuzoCuBnANICVrAUNnbDHbM4NIrqGiE4D8AUAV+S1nqvHXINoKPyFvNYxhwjD73JnvQ0aQojXA7gHwFXayDcXENHyqov1VABbhRCb+10nFSHEhQBeIKInQpQ3dFvj2TbnNvBnAP4KwKe6WB0rcfUUQnwIwIUA3kl9illNcC/zxPMATlN+PhXAkT7VpRAIIUYQifoXiOjL/a6PCyJ6WQhxCNH8RZ4mpt8O4H1CiPcAWAfgDUKIu4joA2kKGzqL3YUQ4meVH98H4Nl+1cWFEOICAP8ZwPuI6LV+12fAeBzAzwoh3iyEqAL4NQD39rlOA4sQQgD4EwDPENF/6Xd9TAghxmTkmBCiBuB85OzdJqK9RHTq6hajvwZgJq2oAyzsOn+w6kp4EsAORDPUeeQWAD8G4IHV0MzP9btCOkKIi4UQzwMYB/BXQoiv97tOALA66XwFgK8jmuj7CyJ6ur+16kQI8ecAZgG8RQjxvBDi0n7XycLbAXwQwPbVtvjtVaszT/wUgIOr7/XjiHzsmcIJ8w6vPGUYhikYbLEzDMMUDBZ2hmGYgsHCzjAMUzBY2BmGYQoGCzvDMEzBYGFnGIYpGCzsDMMwBYOFnWEYpmD8fyNoQx08GaZ+AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def plot_all():\n",
" errors = 0\n",
" for i in range(800):\n",
" if kmeans.labels_[i] == 1 and i<400:\n",
" plt.plot(proj[i][0], proj[i][1], '.', color='blue')\n",
" elif kmeans.labels_[i] == 0 and i<400:\n",
" plt.plot(proj[i][0], proj[i][1], 'X', color='blue')\n",
" errors += 1\n",
" if kmeans.labels_[i] == 0 and i>=400:\n",
" plt.plot(proj[i][0], proj[i][1], '.', color='red')\n",
" elif kmeans.labels_[i] == 1 and i>=400:\n",
" plt.plot(proj[i][0], proj[i][1], 'X', color='red')\n",
" errors += 1\n",
" print ('error rate:', errors/800)\n",
" plt.plot()\n",
"plot_all()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Allelic Statistical Testing "
]
},
{
"cell_type": "code",
"execution_count": 223,
"metadata": {},
"outputs": [],
"source": [
"samples = [\n",
" ('A', 'A'), ('A', 'B'), ('A', 'A'), ('A', 'A'), ('A', 'A'), ('A', 'A'), ('A', 'B'), ('A', 'B'), ('B', 'B'), ('B', 'B')\n",
"]\n",
"phenotypes = ['NO', 'NO', 'NO', 'NO', 'NO', 'YES', 'YES', 'YES', 'YES', 'YES']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Contingency table"
]
},
{
"cell_type": "code",
"execution_count": 227,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"NO ('A', 'A')\n",
"NO ('A', 'B')\n",
"NO ('A', 'A')\n",
"NO ('A', 'A')\n",
"NO ('A', 'A')\n",
"YES ('A', 'A')\n",
"YES ('A', 'B')\n",
"YES ('A', 'B')\n",
"YES ('B', 'B')\n",
"YES ('B', 'B')\n",
"defaultdict(<class 'int'>, {('NO', 'A'): 9, ('NO', 'B'): 1, ('YES', 'A'): 4, ('YES', 'B'): 6})\n"
]
}
],
"source": [
"from collections import defaultdict\n",
"\n",
"d = defaultdict(int)\n",
"\n",
"for phenotype, sample in zip(phenotypes, samples):\n",
" print (phenotype, sample)\n",
" for allele in sample:\n",
" d[(phenotype, allele)] += 1\n",
"print (d)"
]
},
{
"cell_type": "code",
"execution_count": 231,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total_NO: 10\n",
"Total_YES: 10\n",
"Total_A 13\n",
"Total_B 7\n"
]
}
],
"source": [
"total_NO = d[('NO', 'A')] + d[('NO', 'B')]\n",
"total_YES = d[('YES', 'A')] + d[('YES', 'B')]\n",
"print ('Total_NO:', total_NO)\n",
"print ('Total_YES:', total_YES)\n",
"total_A = d[('NO', 'A')] + d[('YES', 'A')]\n",
"total_B = d[('NO', 'B')] + d[('YES', 'B')]\n",
"print ('Total_A', total_A)\n",
"print ('Total_B', total_B)"
]
},
{
"cell_type": "code",
"execution_count": 259,
"metadata": {},
"outputs": [],
"source": [
"random_snp = ['A'] * total_A + ['B'] * total_B\n",
"def test():\n",
" random.shuffle(random_snp)\n",
" #print (random_snp)\n",
" random_yes = random_snp[:total_YES]\n",
" random_no = random_snp[-total_NO:]\n",
" #print ('random_yes', random_yes)\n",
" #print ('random_no', random_no)\n",
"\n",
"\n",
" observed_Yes_A = random_yes.count('A')\n",
" observed_Yes_B= random_yes.count('B')\n",
" observed_No_A = random_no.count('A')\n",
" observed_No_B = random_no.count('B')\n",
"\n",
" chi2_Yes_A = (6.5-observed_Yes_A)**2/6.5\n",
" chi2_Yes_B = (3.5-observed_Yes_B)**2/3.5\n",
" chi2_No_A = (6.5-observed_No_A)**2/6.5\n",
" chi2_No_B = (3.5-observed_No_B)**2/3.5\n",
" random_sum_all = chi2_Yes_A + chi2_Yes_B + chi2_No_A + chi2_No_B\n",
" #print (random_sum_all)\n",
" return random_sum_all\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 263,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.05707\n"
]
}
],
"source": [
"c = 0\n",
"for i in range(100000):\n",
" chi_square = test()\n",
" if chi_square>=5.49:\n",
" c+=1\n",
" \n",
"print (c/100000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 236,
"metadata": {},
"outputs": [],
"source": [
"chi2_Yes_A = (6.5-4)**2/6.5\n",
"chi2_Yes_B = (3.5-6)**2/3.5\n",
"chi2_No_A = (9-6.5)**2/6.5\n",
"chi2_No_B = (1-3.5)**2/3.5"
]
},
{
"cell_type": "code",
"execution_count": 237,
"metadata": {},
"outputs": [],
"source": [
"sum_all = chi2_Yes_A + chi2_Yes_B + chi2_No_A + chi2_No_B"
]
},
{
"cell_type": "code",
"execution_count": 238,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.4945054945054945"
]
},
"execution_count": 238,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum_all"
]
},
{
"cell_type": "code",
"execution_count": 239,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats.distributions import chi2"
]
},
{
"cell_type": "code",
"execution_count": 240,
"metadata": {},
"outputs": [],
"source": [
"p_value = chi2.sf(sum_all, 1)"
]
},
{
"cell_type": "code",
"execution_count": 241,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.01907632210177841"
]
},
"execution_count": 241,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_value"
]
}
],
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment